日本核泄漏放射性污染物定量分析模型.doc

日本核泄漏放射性污染物定量 分析 定性数据统计分析pdf销售业绩分析模板建筑结构震害分析销售进度分析表京东商城竞争战略分析 模型.doc 日本核泄漏放射性污染物定量分析模型 摘要 2011年3月11日下午1点46分，日本本州岛的仙台港以东130公里处的海域发生了里氏8.8级的强烈地震。此次地震以及其所引发的核泄漏事故对日本的政治、经济、人民的生产生活以及周边的地区产生了巨大的影响。 本文针对日本地震对本国及全球影响的评估问 题 快递公司问题件快递公司问题件货款处理关于圆的周长面积重点题型关于解方程组的题及答案关于南海问题 ，在收集处理大量数据的基础上，建立放射性污染物在空气中扩散的一类基于一阶传递偏微分方程的空气污染模型，并进行数学求解。再利用c-free数学软件以及matlab软件对解中的未知常数进行数据计算与拟合验证。并进一步完善模型，增强其科学性与适用性。依据污染物在空气中的衰减时间长短，将其影响力划分为当前影响力和长远影响力，并建立当前污染程度评价模型和污染趋势预测模型。最后运用层次分析法建立各种危害物综合影响力评价体系模型;并具体计算各评价指标的权重值，进而对如何防御核辐射对环境造成危害提出建议。 第一，我们对日本地震所引发的核泄漏的影响力概念作出了合理的解释。核泄漏的影响力主要表现在核泄漏污染物在空气及水中的传播，本文主要研究其中污染物碘-131在空气中的扩散传播。评估日本核泄漏影响力最重要的是建立污染物随时间、空间扩散的科学、全面、合理的数学模型，同时在借鉴国内外重大事件影响评估的 方法 快递客服问题件处理详细方法山木方法pdf计算方法pdf华与华方法下载八字理论方法下载 ，以及相关数据挖掘的基础之上，我们建立了一类基于一阶传递偏微分方程的空气污染模型: 2,,,NNN,,,DukN20,,,txx 第二，根据以上偏微分方程模型，我们在查找大量科学资料的情况下，根据2,ND日本岛国的气候特点以及其周边地区的空气流动规律，合理省去 项，并最2,x终求得: 0,0,xut,,,0, kxNxt(,),,, ， u0,Pexut,,,,00, 并运用c-free软件的大量实际数据的对比计算以及matlab的图像拟合检验求k,0.1193得 。 u0 2,ND第三，针对更加一般的情况下，考虑 的影响，建立更具有代表性的2,x模型，并通过大量数据的拟合，使模型适用性增强。结论更具有科学性与适用性。 第四，我们以空气中的碘(I-131)含量为方案层，并依据这项指标的内在属性，以日本核泄漏影响力为目标，运用层次分析法建立各种危害物综合影响力评价体系模型。 【关键字】:日本核泄漏 空气污染 一阶传递偏微分方程 数据拟合 层次分析法 目录 1、问题的重述与分析 2 2、模型的假设与符号说明 2.1 模型假设. . . . . . . . . . . . . . . . . . . . . . . . . . . . .2 2.2 符号说明. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .3 3、模型的建立与求解 3.1 引言. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 3.2 建立数学模型 3.2.1 模型导出的原因分析. . . . . . . . . . . . . . . . . . . .4 3.2.2 函数 excel方差函数excelsd函数已知函数     2 f x m x mx m      2 1 4 2拉格朗日函数pdf函数公式下载 的约束简化与建立. . . . . . . . . . . . . . . . . . . .4 3.3 函数的数学求解. . . . . . . . . . . . . . . . . . . . . . . . . . . .4 3.4常数的确定与验证 . . . . . . . . . . . . . . . . . . . . . . . .6 3.4.1常数的确定 . . . . . . . . . . . . . . . . . . . . . 7 3.4.2常数的matlab验证 4、模型的修正 8 5、放射性污染物的层次分析 5.1模型的假设. . . . . . . . . . . . . . . . . . . . . . . .......10 5.2符号的说明. . . . . . . . . . . . . . . . . . . . . . . .......10 5.3 建立层次结构模型. . . . . . . . . . . . . . . . . . . . . . . 11 5.4 计算组合权向量和组合一致性检验. . . . . . . . . . . . . . . . 13 5.5模型的求解. . . . . . . . . . . . . . . . . . . . . . . .......14 6、模型的评价 6.1 模型的优点. . . . . . . . . . . . . . . . . . . . . . . . .. 16 6.2 模型的缺点. . . . . . . . . . . . . . . . . . . . . . . . . . 16 7、参考文献 17 8、附录 18 1 一、问题的重述与分析 2011年3月11日下午1点46分，日本本州岛的仙台港以东130公里处的海域发生了里氏8.8级的强烈地震。震源深度为24公里，引发了10米高的海啸，可以说裹胁了大量的汽车、房屋还有轮船，冲入了日本本州东北部的部分地区。这次地震到目前为止，新华网快讯:日本警察厅1日宣布，截至当天10时，日本东北部海域发生的强震及其引发的海啸已确认造成11578人死亡，16451人失踪.核辐射的阴云已从陆地飘向空中，截至目前，除日本外，已有11个国家相继检测到微量放射性物质;更严重的是，反应堆中的高辐射积水已经快要超标，海水也面临着难以想象的污染危机。 作为日本有仪器记录以来的最高震级的地震，这次的地震包括所引发的核泄漏对日本的政治、经济、人民的生活以及周边的地区产生了巨大的影响。由于本次地震牵涉面广，影响巨大，我们无法一一做出定量评估，因此我们希望通过数学建模，在广泛收集各种有价值的数据资源的情况下，建立核泄漏中核污染物在空气中扩散的偏微分方程，并求得一定精确度的解，从而对地震所引发的核泄漏问题进行定量评估分析。 为了定量评估核泄漏的影响情况，我们主要从以下几个方面进行研究: (1)分析情况，对问题进行合理的简化，最终确定研究重点，展开更加深入的分析，建立合适的数学模型。 (2)根据情景，建立污染物在空气中扩散的一类基于一阶传递偏微分方程的空气污染模型，并进行数学求解。再利用c-free数学软件以及matlab软件对解中的未知常数进行数据计算与拟合验证。使该方程能很好的反映本次核泄漏对日本各城市以及周边地区在近几日以及未来一段时间的污染扩散情况，有利于及时对污染情况做出初步的预计估测。 2,N(3)针对更加一般的情况下，考虑 的影响，建立更具有代表性的模型，D2,x并通过大量数据matlab拟合，使模型适用性增强。结论更具有科学性与适用性。 (4)分析核泄漏所产生的各种危害物，确定各因素的影响力，运用层次分析法 建立综合影响力评价体系模型;并具体计算各评价指标的权重值，进而对如何防御核辐射对环境造成危害提出建议。 二、模型的假设与符号说明 2.1 模型的假设 (1)假设文中所用到的所有数据都是真实可靠的，能反映真实情况。 (2)假设在近一段时间内，日本及周边地区的气候不发生显著变化。 (3)假设日本已经较好的控制住核辐射的泄露,且核泄漏危害物质以其自身的衰减速率递减。 (4)由日本岛国的气候特点，假设日本及周边地区空气流动速率远大于放射性污染物分子热运动的扩散速度。 (5)初始时刻探测空间区域中污染物碘元素的浓度为0mg/L。 2 2.2 符号说明 N(x,t) 污染物碘-131的浓度, 单位为mg/L; t 距核泄漏污染时间, 单位为h; x 检测地距核泄漏事发地的距离, 单位为km; u 纵向空气流速, 单位为km /h; 02kmh/ D 污染物碘-131的扩散系数, 单位为 ; ,1 k 污染物碘的衰减速率系数, 单位为 ; h P 城市空气中碘元素的浓度值，单位为mg/L; 三、模型一的建立与求解 3.1 引言 本大地震的爆发，以及其所引起的核泄漏、核辐射问题，环境污染事 随着日 其中由核辐射的主要载体空气所引发的污染事故是日本及日本故已经日趋严重, 周边地区环境污染事故的主要类型. 采用数学模型进行空气质量模拟计算具有灵活、快速、可操作性强等优点, 有助于决策部门了解污染带的迁移状况和污染 从而对事故物在时间、空间上的变化, 掌握污染物对各个城市造成的污染影响, 的发展做出及时、准确的反应. 这类数学模型一般是以反应扩散方程的形式进行模拟. 反应扩散方程(Reaction - diffusion equation)的基本形式是[2]: ,P,,，DPftxP(,,),x, 这里的P是日本核泄漏地碘元素密度函数; D 是扩散系数; DΔP是扩散项; 是反 ftxP(,,)应项; 在这里 2n,P(),,,,,PdivP,. 2,xi1i 在对空气污染扩散问题的研究中发现, 在空气流动速度迅速的地区或者污染物扩散受空气温度影响的情况下, 有时污染物自身的扩散速度相较于空气流动速度微不足道, 这时没有扩散项的一阶传递方程( Firstorder transport equation) 数学模型更能精确的模拟出空气污染的真实情况. 因此我们所要建立的数学模型就是一类基于一阶传递偏微分方程的空气污染模型. 选取近期日本大地震所引起的核泄漏污染事故为研究对象, 建立一阶传递方程的空气污染模型, 并对其进行求解. 3.2 建立数学模型 2011年3月11日日本大地震导致核泄漏引发日本及其周边地区空气污染. 选取以核泄漏地一定半径的区域为研究对象,由国家环保局检测出微量I-131， 具有较可靠的具体实验数据。 3 3.2.1 模型导出的原因分析 (1)在以核泄漏地为原点的一定半径区域内碘元素在整个地区空气的横断面上已经基本混合均匀[ 6 ] ; (2)由化学知识可知碘元素在空气中具有一定的衰减速率[ 7, 8 ] ; (3)假设日本已经较好的控制住核辐射的泄漏,空气中没有另外的碘元素增加。 (4)由于在以核泄漏地为原点的一定半径区域内碘元素在整个地区空气的横断面上已经基本混合均匀，其之间的微小差异比之日本核泄漏地到检测地的空间球面距离可以忽略不计; 由以上各种因素分析可知，我们可以用一维模型进行碘元素含量预测: 2,,,NNN,,,DukN20,,,txx 3.2.2 函数的约束简化与建立 查找资料可知，日本列岛受有规律的季风影响，冬季寒流南下，夏季热带气流北上。在这两种气流影响下，冬夏两季各长达数月，而处于季风交替之际的春秋季，各为两个月左右。因此，在日本及其周边地区空气流动速度较快的空间范2,N围内，我们不考虑 项,可以建立没有扩散项的一阶传递方程的模型进行D2,x空气质量模拟[ 1 ] : ,,NN,,,,,,ukNxt,0,0,0,,,tx,,Nxx(,0)0,0,,,(1),, ,NtPt(0,)0,0,,,,,Ntt(,)0,0,,,,,, ,N其中符号表示的意义见符号说明.由Fick定律知流量( flux) 与梯度( ) 成正比,指向相反的方向(因为流动总是从高密度区域流向低密度区域) ,另外还要,x减掉污染物碘元素自身的衰减( kN ) ,这样就得到了一阶传递方程. 初始时刻这段空间区域中碘元素的浓度为0mg/L, 这确定了模型的初值条件.任意时刻, 城市空间污染物的浓度值为P并且距离探测城市很远的区域空气中不含有碘元素, 因此确定了模型的边值条件. 3.3 函数求解 求解一阶传递方程的一维模型(1) [ 3,5 ] . 求解过程: zsNxusts()(,),，，令 (向量函数) , 那么可以得到 0 ,，，,，，NxustsNxusts(,)(,)'00zsu(),,，0,,xt，. 4 'zskNxustskzs()(,),,，，,,，，由(1) 第一个方程得 解得 0,kszsze()(0),, (2) zNxtztNxut(0)(,),()(,0),,,,由定义的向量函数得 ，又由(2)得 0 ,,ktkt，，，，ktztzezeNxte()(0)(0)(,),,,,,,,。 ktNxutNxte(,0)(,),,,这样就得到 ,即 0 ,ktNxtNxute(,)(,0),,,0。 Nxut(,0)0,,Nxt(,)0,xut,,0 当 时,由(1) 第二个方程得 , 因此 。 00 xxzNt()(0,),,,xut,,0 当 时, 由定义的向量函数得 . 0uu00 同理可得 xkxkx,,()kxuuu000()(0)(0)(,),,,,,,,zzezeNxteu0。 这样就得到 kxxu0(0,)(,)NtNxte,,,u0, 即 kxxu0(,)(0,)NxtNte,,,u0。 再由(1) 第三个方程得 kx,u0NxtPe(,),,。 因此方程组(1) 的解为 0,0,xut,,,0,kxNxt(,)(3),。,,u0,Pexut,,,,00, xx,xxut(), 由(3) 式可以看出, 碘元素的浓度以 点为临界点, 在 时1110 碘元素的浓度在随着x的增大而减小, 即空气中的污染物随着空气流动距离的增 xx,加而不断减少; 在 时碘元素的浓度变成0mg/L,即在x1点之后的空气中不1 含有污染物. xx,因此决策部门应该采取措施大力控制 这一距离的污染物.另外在结果 1 x1x中也可以看出对于一个固定的地点 ,当 时碘元素的浓度为0, 即空 0t, x1u0气的这一横截面在时间 之前不含有污染物. t, u0 5 3.4常数的确定与验证 3.4.1常数的C-Free确定 环境保护部(国家核安全局)于3月31日在距离核泄漏中心不同距离的监测点进行气溶胶取样中检测到了微量的人工放射性核素碘-131，浓度均在10-4贝克/立方米量级及以下;另还在部分地区检测到了极微量放射性核素铯-137和铯-134，浓度均在10-5贝克/立方米量级及以下。 结合近年来辐射环境监测数据分析，初步确认各地所检测到的人工放射性核素来自日本福岛核事故。因此我们查找了同一时间，不同位置空气中所含的碘浓度. 我们由国家环境保护部公布的数据中，选取了人工放射性元素碘-131在距离核泄漏中心不同球面距离的碘-131浓度分布情况，其具体数据见下表一: 核泄漏中心附近范围的I-311浓度分布 距核泄漏中心的球面距离(km) 空气中131I浓度(mBq/m3) 200 1 400 0.9 600 0.8 800 0.72 1000 0.68 1200 0.63 1600 0.56 1800 0.45 2000 0.36 2200 0.3 2400 0.2653 2600 0.21 2800 0.19 3000 0.15 3200 0.139 3400 0.09 3600 0.0565 3800 0.056 4000 0.04 表一 6 k ,0.1193根据以上数据，我们运用C-Free数学软件计算出了 ，P=1.15390u(具体算法见附录一)，故此我们得到了空气中碘元素扩散趋势模型，其公式,0.1194xNxte(,)1.154,，为: ,运用matlab可画出其三维空间中的浓度变 化曲线大致为(如图一): 图一 3.4.2常数的matlab验证 利用matlab软件进行数据拟合。 选用指数趋势模型对日本大地震引发的核泄漏中污染物碘元素在空气中扩散情k况进行拟合，得到了碘元素扩散的指数趋势模型，验证有关参数(P, )，以及0u基于一阶传递偏微分方程的空气污染趋势模型。模型拟合效果见图二: 空气中碘元素浓度拟合曲线 图二 7 同时，我们由matlab可得模型拟合函数形式如下: ,0.1194xNxte(,)1.154,，， 即由matlab拟合所得数据与c-free软件求得数据极为近似，因此，我们得出结论:不管是我们用C-Free软件得到的碘元素扩散趋势模型，还是用matlab软件进行数据拟合得到的碘元素扩散趋势模型，都很符合实际空气碘元素的扩散情况。 四、模型的修正 在模型一中，我们根据日本岛国的气候特点，分析建立了适用于类似日本气候规律地区的空气中污染物扩散的一类基于一阶传递偏微分方程的约束型空气污染模型(特点为受有规律的季风影响，在日本及其周边地区空气流动速度相对2,N较快，因此我们不考虑 项)，然而若对模型一进一步推广，使得其科学D22,N,x性、普适性更好,我们必须考虑 项的影响.故我们建立模型二: D 22,x,,,NNN,,,DukN20,,,txx ， 在收集大量不同时间(t)，不同距核泄漏中心的球面距离(x)检测数据的情况下，运用计算机拟合的方法，绘制出关于N(x,t),x,t之间的三维污染物-碘131空气中扩散的图形。 下面是我们收集的实际观测数据(见表二): 空气中131I浓度(mg/l) 距核泄漏中心的球面距离 (km) 31日 27日 23日 19日 15日 12日 200 1 1.31 1.80 2.62 3.43 4.49 400 0.9 1.21 1.64 2.35 3.1 4.05 600 0.8 1.1 1.45 2.08 2.72 3.6 800 0.72 1 1.295 1.875 2.448 3.3 1000 0.68 0.87 1.224 1.768 2.312 3.1 1200 0.63 0.816 1.133 1.638 2.142 2.835 1600 0.56 0.728 1.007 1.453 1.906 2.52 1800 0.45 0.586 0.86 1.17 1.56 2.02 2000 0.36 0.469 0.649 0.938 1.224 1.63 2200 0.3 0.39 0.545 0.785 1.06 1.35 2400 0.265 0.3449 0.477 0.689 0.903 1.194 2600 0.21 0.274 0.385 0.546 0.716 0.946 2800 0.19 0.246 0.342 0.495 0.645 0.855 3000 0.15 0.199 0.273 0.395 0.53 0.674 3200 0.139 0.182 0.250 0.361 0.475 0.626 3400 0.09 0.117 0.162 0.234 0.306 0.405 3600 0.057 0.0735 0.102 0.147 0.193 0.254 3800 0.056 0.0726 0.101 0.145 0.191 0.252 4000 0.04 0.051 0.076 0.14 0.136 0.185 表二 8 由于考虑扩散项的一阶传递方程( Firstorder transport equation) 数学2,,,NNN模型 ，所以我们不能轻易求得其简便的数学解。 ,,,DukN02,,,txx 因此我们只能求助于matlab这一数学软件进行数据拟合，从而由拟合后的图形对污染趋势进行预测，得到较精确的预测数据，给人们的灾后防治工作带来帮助。下面就是我们借由数据所绘制的拟合图形: (1)根据上述观测数据，我们先由matlab拟合出众多分立的曲线，每条曲线代表一天的碘元素浓度分布情况:(见图三) 图三 9 (2)在适当修改程序，绘制出三维坐标下的关于不同时间t，不同距离x的 碘-131浓度曲面，该曲面具有一定的真实性，是对污染物碘-131扩散趋势的一种较好表现:(见图四) 图四 五(模型一的建立与求解 5.1模型的假设 (1)假设正常情况下的人群在核泄漏之前都是健康的; (2)假设文中所用到的所有数据都是真实可靠的，能反映真实情况; (3)假设这些放射性元素只存在于空气，水体和食物当中，其它地方的含量忽略不计; (4)假设日本已较好地控制住核辐射的泄露,空气中没有另外的放射性元素增加; 5.2符号的说明 A:碘-131，铯-137，锶-90，和钚对人体健康的影响的判断矩阵; :碘-131对空气，水体和食物的污染程度的比较矩阵; C1 :铯-137对空气，水体和食物的污染程度的比较矩阵; C2 :锶-90对空气，水体和食物的污染程度的比较矩阵; C3 10 :钚对空气，水体和食物的污染程度的比较矩阵; C4 w :权重向量; i aijaa :元素与相对于准则C的重要性的比例标度; ij CI:一致性指标; RI:平均随机一致性指标; CR:一致性比例; 5.3 建立层次结构模型 5..31准则层的确定和模型建立 作为核泄漏主要的四种污染物碘-131 ，锶-90，铯-137和钚，这些放射性元素对人体健康有比较大的影响。当这些元素的剂量超过一定浓度时会对人体造成伤害，严重则直接导致死亡。 据此，我们确定了以下三个准则。 1. 碘-131。碘-131的半衰期是八天，意味着它需要数月时间才会完全消失。美 国环保局称，放射性的碘可通过摄入或吸入进入人体。它可能落在草上被牛 食入，并最终通过牛奶进入人体，也可能落在带叶蔬菜上或聚集到海鱼和淡 水鱼体内，被人食用。 核反应堆完全熔解后，会释放出另一种危险物质铯-137。铯-137的半衰期为 30年。人会通过食物和水将其摄入，或类似灰尘被人吸入。与其接触，会增 加癌症的患病率。 2. 锶-90。锶-90的半衰期长达29.1年。锶-90主要通过食物和水进入人体，摄 入锶-90会导致骨癌(骨骼附近软组织的癌症)和白血病。 3. 铯-137。铯-137的半衰期为30年。人会通过食物和水将其摄入，或类似灰尘 被人吸入。与其接触，会增加癌症的患病率。 4. 钚。钚的半衰期最长可达2万年，毒性非常大。 在前面分析的基础上，建立 2011 年核泄漏污染物对人体健康的影响的层次结构，利用层次分析法，构造了层次分析模型。(图) 目核泄漏污染物对人体健康的影响程度 标 层 准碘-131 锶-90 铯-137 钚 则 层 子 准 则 空气 水体 食物 层 11 5.32.构造判断(成对比较)矩阵 我们有比例标度值(如表) 表三: 比例标度值 标度 含 义 aij 与的影响相同 1 CCij 3 比的影响稍强 CCij 5 比的影响强 CCij 7 比的影响明显地强 CCij 9 比的影响绝对地强 CCij 2,4,6,8 与的影响之比在上述两个相邻等级之间 CCij 11a 与的影响之比为上面的互反数 CCiji,？,j29 根据表三，我们构造出成对比较矩阵。 1a,,a,1,(i,j,1,2,？,n)jiiiaij aa？a,,11121n,, aa？a,,21222nAa,,，，ij,, n，n？？？？,,,,aa？a12,,nnnn 我们结合表三，再根据碘-131，铯-137锶-90和钚对人体健康的影响以及CCC243这些放射性污染物空气，水体和食物的污染程度构造出A，，，，。 CCCC1243 1,,13,,59,,,,111,,,,,,131,,11,,112,,125,,,,,,51673,, 111,,,,,,,,,C,1124C2,713C,1C,1334,,1,,A11,,33,,2,,1111,,,,131,,,,111，,,131,,22,,3,,1,,,,36853,,,,9481,, 5.33 .判断矩阵的一致性检验及权重的计算 我们根据个元素对于准则的判断矩阵，求出他们对准则C的相对排序权重 Tw,w,w？ww,(w,w,w？w)，写成权向量形式，即。这里我们首先要通123n123n 12 过判断矩阵的一致性检验和计算权重。 一致性检验我们通过计算: ,n,max(1) 一致性指标:。 CI,n,1 随机一致性指标:，可查找到相应的平均随机一致性指标(如表4)。 (2) RI 1 2 3 4 5 6 7 8 9 10 11 12 n 13 14 15 0 0 0.58 0.90 1.12 1.24 1.32 1.41 1.45 1.49 1.51 RI 1.54 1.56 1.58 1.59 表4:随机一致性指标 CI(3) 一致性比率指标:，当时，认为判断矩阵的一致性是CR,0.10CR,RI ,可以接受的，则对应的特征向量可以作为排序的权重向量。此时 max n aw,ijjnn，，,AW1,1ji,,,,,maxnwnw,1,1iiii (A,W)其中表示的第个分量。 iA,Wi 在权重计算上，我们取判断矩阵个列向量归一化后的算术平均值，近似作n na1ijw,(i,1,2,？,n)为权重，即 ,inn,1ja,kj,1k 类似地，也可以对按行求和所得向量作归一化，得到相应的权重向量。 5.4 计算组合权向量和组合一致性检验 5.41 计算组合权向量 n设第层上个元素对总目标(最高层)的排序权重向量为 k,1k,1 T(k,1)(k,1)(k,1)(k,1)，，W,w,w,？,w 12nk,1 n第层上个元素对上一层(层)上第个元素的权重向量为 jkk,1k T(k,1)(k)(k)(k)，，P,p,p,？,p,j,1,2,？,n j1j2jnjk,1k则矩阵 (k)(k)(k)(k),,P,P,P,？,P 12nk,1n，n是阶矩阵，表示第层上的元素对第层各元素的排序权向量。那么kk,1kk,1 第层上的元素对目标层(最高层)总排序权重向量为 k 13 ,,(k)(k)(k1)(k)(k)(k)(k1),,WPWPPPW,,,,,？,,n12,1k T(k)(k)(k)，，,w,w,？,w12nk 或 nk,1kkk()()(,1) w,pw,i,1,2,？,n,iijjkj,1 对任意的有一般公式 k,2 (k)(k)(k,1)(3)(2)W,P,P,？,P,W(k,2) (2)W其中是第二层上各元素对目标层的总排序向量。 5.42 组合一致性 (k)(k)(k)CI,CI,？,CI 设层的一致性指标为，随机一致性指标为 k12nk,1(k)(k)(k)RI,RI,？,RI 12nk,1则第层对目标层的(最高层)的组合一致性指标为 k (k)(k)(k)(k)(k,1)，，CI,CI,CI,？,CI,W 12nk,1组合随机一致性指标为 (k)(k)(k)(k)(k,1)，，RI,RI,RI,？,RI,W 12nk,1组合一致性比率指标为 (k)CI(k)(k1), CR,CR，(k,3)(k)RI (k)CR,0.10当时，则认为整个层次的比较判断矩阵通过一致性检验。 5.5模型的求解 通过Matlab软件，我们计算得到了准则层判断矩阵(表五)，碘131的C1 层次单排序(表六)，铯137的层次单排序(表七)，锶-90的层次单排序CC23 (表八)，钚的层次单排序(表九)以及层次总排序结果(表十)。 C4 准则层判断矩阵 CCCC污染物影响 核辐射伤害程度排序 3124 碘-131 C1 3 0.0774 111 59 锶-90 C5 1 6 0.2518 12 4铯-137 C1 0.0439 1113 368 钚 C9 4 8 1 0.6269 4 表五 14 碘131的层次单排序 C1 SSS碘131影响程度 排序结果 312 空气 S1 2 5 0.5816 1 水体 S1 3 0.3090 12 2 食物 S1 0.1095 113 53 ,,3.0037CI,0.00185CR,0.0036,0.1 max 表六 铯137的层次单排序 C2 SSS铯137影响程度 排序结果 312 空气 S1 3 1 0.4286 1 水体 S1 0.1428 112 33食物 S1 3 1 0.4286 3 ,,2.9999CI,0CR,0.1 max 表七 锶-90的层次单排序 C3 SSS锶-90影响程度 排序结果 312 空气 S1 1 2 0.4 1 水体 S1 1 2 0.4 2 食物 S1 0.2 113 22 ,,3CI,0CR,0,0.1 max 表八 钚的层次单排序 C4 SSS钚影响程度 排序结果 312 空气 S1 0.088 111 73水体 S7 1 3 0.6694 2 食物 S3 1 0.2426 13 3 ,,3.006CI,0.003CR,0.006,0.1 max 表九 15 层次总排序结果 CCCC影响程度 子准则层对目标1111 层的权重 权重系数 0.0774 0.2518 0.0439 0.6269 S0.5816 0.4286 0.4 0.088 0.2556 1 S0.3090 0.1428 0.4 0.6694 0.4968 2 S0.1095 0.4286 0.2 0.2426 0.2775 3 CI,0.006,0CR,0.003,0.1 总 表十 从这些表格的数据中，我们可以知道目前的辐射影响程度从大到小分别是锶-90，钚，碘-131，铯-137。而三个载体中，受这些放射性污染物影响最大的是水体，其次是食物，最后是空气。这个结果与实际情况很符合，因为碘-131的挥发性较强，在空气中会迅速挥发，并且它的半衰期相对较短，只有8天，所以它的危害程度不是最高的。在三个载体中，空气由于体积庞大，对污染物的稀释作用较好，随着时间的推移，它受污染的程度会减弱。 六、模型的评价 6.1 模型的优点 (1)借助于matlab数学软件， 本文建立了一类基于一阶传递偏微分方程的空气污染模型并得到了它的精确解。这些精确解不仅包括了以前的 论文 政研论文下载论文大学下载论文大学下载关于长拳的论文浙大论文封面下载 中求得的类孤立波解, 而且较之更具有一般性, 因而扩展了该类反应扩散方程解的范围。 )该模型不仅能够求解污染物-碘的污染扩散情况，还能够运用于其他污染物(2 针对不同时间、距离扩散情况的预测与分析，简答明了。 (3)通过计算机数据拟合，绘出的图像能够更好的验证建模求解结果的正确性，有一定科学性和说服力。 6.2 模型的缺点 (1)本次建模所收集到的数据量不是特别充足，对所建模型的求解具有一定的影响，使得评价具有一定的偏差; (2)现实自然是复杂的,其复杂性就在于影响因素的众多和不确定性，并且由于各影响因 素的相互作用及内在的有机联系，在对核辐射污染物扩散的预测时，只对空气中碘元素进行分析，有些片面，且评估技术存在误差; (3)运用传递偏微分方程分析的方法得出的预测函数只能对短期的发展进行预测，若要对长期发展进行预测，还需考虑气候等其他因素的影响。 它的的局限性主要表现在: (4)从建立层次结构模型到给出成对比较矩阵，人主观因素对整个过程的影响 很大，这就使得结果难以让所有的决策者接受。该模型中的比较、判断以及结果的计算过程都是粗糙的，精确度不是很高。 16 七、参考文献 [ 1 ] Robert B. Banks. Growth and Diffusion Phenomena: mathe2 matical frameworks and app lications, volume 14 of Texts in ap2 p lied mathematics. Sp ringer - Verlag, Berlin Heidelberg, 1994. [ 2 ] Junp ing Shi. Reaction diffusion equations and mathematical bi2 ology. Undergraduate Lecture Notes, 2004. [ 3 ] Lawrence C. Evans. Partial differential equations, volume 19 of Graduate Studies in Mathematics. American Mathematical Society, Providence R I, 1998. [ 4 ] 伍卓群, 尹景学, 王春朋. 椭圆与抛物型方程引论. 北京: 科学出版社, 2003. [ 5 ] 陈恕行. 现代偏微分方程导论. 北京: 科学出版社, 2005. [ 6 ] 张波, 王桥, 李顺, 孙强, 王利强, 富尔江. 基于系统动力 学模型的松花江水污染事故水质模拟. 中国环境科学, 2007, 27 (6) : 811,815. [ 7 ] 袁星, 郎佩珍. 硝基苯在江水中生物降解动力学模拟研究环境化 学, 1991, 10 (6) : 24 - 30. [ 8 ] 刘玉萍, 王丽红, 刘薇, 刘侨博. 松花江有机污染物的污 染特征研究. 环境科学与管理, 2006, 31 (3) : 73,75. 17 附录一: #include #include int main(){ double C[25],X[25],fa[500]; int i,j,k; freopen("data.txt","r",stdin); for(i=1;i<=23;i++) scanf("%lf %lf",&C[i],&X[i]); for(i=1,k=0;i<23;i++) for(j=i+1;j<=23;j++){ fa[k++]=log(C[j]/C[i])/(X[i]-X[j]); } //freopen("res.txt","w",stdout); freopen("ave.txt","w",stdout); double sum=0; int cnt=0; for(i=0;i0)sum+=fa[i],cnt++; } printf("ave = %lf\n",sum/cnt); return 0; } 解得 : ave = 0.1193 附录二: 本例程序采用MATLAB 程序源代码 maxlmta: function [w,lmta,CI,RIn,flag]=maxlmta(A) RI=[0,0,0.58,0.90,1.12,1.24,1.32,1.41,1.45,1.49,1.51]; 18 %将A的每一列向量归一化 Asum=sum(A); n=length(A); for j=1:n for i=1:n Wij(i,j)=A(i,j)./Asum(1,j); end end %将Wij按行求和 Wij=Wij'; Wi=sum(Wij); Wi=Wi'; %将Wi归一化 W=sum(Wi); w=Wi./W; %计算lmta lmta=sum(1/n*(A*w)./w); %计算CR RIn=RI(1,n); CI=(lmta-n)/(n-1); CR=CI/RIn; %判断一致性检验 if CR<0.1 disp('通过一致性检验') flag=1; else disp('不能通过一致性检验') flag=0; end AHP: function AHP() %参数赋予初值 numprilay = 0; numevepri = 0; m = 0; n = 0; 19 i = 1; lmta2 = 0; CI2 = 0; RI2 = 0; lmta3 = 0; CI3 = 0; RI3 = 0; lmta4 = 0; CI4 = 0; RI4 = 0; lmtap = 0; CIp = 0; RIp = 0; %初始化各个参数 numprilay=input('输入准则层层数:'); n=numprilay; while n>0 m=input(['输入第',num2str(i),'准则层准则数:']); numevepri=[numevepri,m]; n=n-1; i=i+1; end numevepri=numevepri(1,2:numprilay+1); %除去第一个无意义元素 numpla=input('输入方案层数方案数:'); %初始化第1准则层对目标层的成对比较矩阵并检验一致性，若不通过需重新输 入 flag1=0; while flag1~=1 if numprilay>=1 A1=input('输入第1准则层对目标层的成对比较矩阵:'); [w1,lmta1,CI1,RI1,flag1]=maxlmta(A1); if flag1==0 disp('重新输入') end end end %若存在第2准则层，计算其各个参数 if numprilay>=2 flag2=zeros(1,numevepri(1)); %设置标志位，检测第几个矩阵不能通过一致 性检验 %判断是否全部通过一致性检验，若不通过这重新输 20 while sum(flag2)~=numevepri(1) A2=input('输入第2准则层对第1准则层的成对比较矩阵:'); A2size=size(A2,2); w2=zeros(A2size,1); %初始化权向量 for i=1:numevepri(1); a2=A2((i-1)*A2size+1:i*A2size,:); %提取每个成对比较矩阵 [w2temp,lmta2temp,CI2temp,RI2temp,flag2temp]=maxlmta(a2); w2=[w2,w2temp]; %生成权向量 lmta2=[lmta2,lmta2temp]; %生成最大特征根 CI2=[CI2,CI2temp]; %生成CI RI2=[RI2,RI2temp]; %生成RI flag2(i)=flag2temp; %生成标志位 end w2=w2(1:A2size,2:numevepri(1)+1); %除去第一个无意义元素 lmta2=lmta2(1,2:numevepri(1)+1); %除去第一个无意义元素 CI2=CI2(1,2:numevepri(1)+1); %除去第一个无意义元素 RI2=RI2(1,2:numevepri(1)+1); %除去第一个无意义元素 %显示哪些不能通过一致性检验的矩阵 for i=1:numevepri(1) if flag2(i)==0 disp(['需要重新输入第',num2str(i),'个成对比较矩阵']); end end end end %若存在第3准则层，计算其各个参数 if numprilay>=3 flag3=zeros(1,numevepri(2)); %设置标志位，检测第几个矩阵不能通过一致性检验 %判断是否全部通过一致性检验，若不通过这重新输入 while sum(flag3)~=numevepri(2) A3=input('输入第3准则层对第2准则层的成对比较矩阵:'); A3size=size(A3,2); w3=zeros(A3size,1); %初始化权向量 for i=1:numevepri(2); a3=A3((i-1)*A3size+1:i*A3size,:); %提取每个成对比较矩阵 [w3temp,lmta3temp,CI3temp,RI3temp,flag3temp]=maxlmta(a3); w3=[w3,w3temp]; %生成权向量 lmta3=[lmta3,lmta3temp]; %生成最大特征根 CI3=[CI3,CI3temp]; %生成CI RI3=[RI3,RI3temp]; %生成RI 21 flag3(i)=flag3temp; %生成标志位 end w3=w3(1:A3size,2:numevepri(2)+1); %除去第一个无意义元素 lmta3=lmta3(1,2:numevepri(2)+1); %除去第一个无意义元素 CI3=CI3(1,2:numevepri(2)+1); %除去第一个无意义元素 RI3=RI3(1,2:numevepri(2)+1); %除去第一个无意义元素 %显示哪些不能通过一致性检验的矩阵 for i=1:numevepri(2) if flag3(i)==0 disp(['需要重新输入第',num2str(i),'个成对比较矩阵']); end end end end %若存在第4准则层，计算其各个参数 if numprilay>=4 flag4=zeros(1,numevepri(3)); %设置标志位，检测第几个矩阵不能通过一致性检验 %判断是否全部通过一致性检验，若不通过这重新输入 while sum(flag4)~=numevepri(3) A4=input('输入第4准则层对第3准则层的成对比较矩阵:'); A4size=size(A4,2); w4=zeros(A4size,1); %初始化权向量 for i=1:numevepri(3); a4=A4((i-1)*A4size+1:i*A4size,:); %提取每个成对比较矩阵 [w4temp,lmta4temp,CI4temp,RI4temp,flag4temp]=maxlmta(a4); w4=[w4,w4temp]; %生成权向量 lmta4=[lmta4,lmta4temp]; %生成最大特征根 CI4=[CI4,CI4temp]; %生成CI RI4=[RI4,RI4temp]; %生成RI flag4(i)=flag4temp; %生成标志位 end w4=w4(1:A4size,2:numevepri(3)+1); %除去第一个无意义元素 lmta4=lmta4(1,2:numevepri(3)+1); %除去第一个无意义元素 CI4=CI4(1,2:numevepri(3)+1); %除去第一个无意义元素 RI4=RI4(1,2:numevepri(3)+1); %除去第一个无意义元素 %显示哪些不能通过一致性检验的矩阵 for i=1:numevepri(3) if flag4(i)==0 22 disp(['需要重新输入第',num2str(i),'个成对比较矩阵']); end end end end %初始化方案层对最后一个准则层的成对比较矩阵，计算其各个参数 flagp=zeros(1,numevepri(numprilay)); %设置标志位，检测第几个矩阵不能通过一致性检验 %判断是否全部通过一致性检验，若不通过这重新输入 while sum(flagp)~=numevepri(numprilay) P=input(['输入方案层对第',num2str(numprilay),'准则层的成对比较矩阵:']); psize=size(P,2); wp=zeros(psize,1); %初始化权向量 for i=1:numevepri(numprilay); p=P((i-1)*psize+1:i*psize,:); %提取每个成对比较矩阵 [wptemp,lmtaptemp,CIptemp,RIptemp,flagptemp]=maxlmta(p); wp=[wp,wptemp]; %生成权向量 lmtap=[lmtap,lmtaptemp]; %生成最大特征根 CIp=[CIp,CIptemp]; %生成CI RIp=[RIp,RIptemp]; %生成RI flagp(i)=flagptemp; %生成标志位 end wp=wp(1:psize,2:numevepri(numprilay)+1); %除去第一个无意义元素 lmtap=lmtap(1,2:numevepri(numprilay)+1); %除去第一个无意义元素 CIp=CIp(1,2:numevepri(numprilay)+1); %除去第一个无意义元素 RIp=RIp(1,2:numevepri(numprilay)+1); %除去第一个无意义元素 %显示哪些不能通过一致性检验的矩阵 for i=1:numevepri(numprilay) if flagp(i)==0 disp(['需要重新输入第',num2str(i),'个成对比较矩阵']); end end end %计算组合权向量和组合一致性检验 switch numprilay case 1 w=w1'*wp'; %计算组合权向量 CRfinal=CI1/RI1+(CIp*w1)/(RIp*w1); %计算组合一致性检验 case 2 23 w=w1'*w2'*wp'; %计算组合权向量 CRfinal=CI1/RI1+(CI2*w1)/(RI2*w1)+(CIp*w2)/(RIp*w2); %计算组合一 致性检验 case 3 w=w1'*w2'*w3'*wp'; %计算组合权向量 CRfinal=CI1/RI1+(CI2*w1)/(RI2*w1)+(CI3*w2)/(RI3*w2)+(CIp*w3)/(RIp*w3); % 计算组合一致性检验 case 4 w=w1'*w2'*w3'*w4'*wp'; %计算组合权向量 CRfinal=CI1/RI1+(CI2*w1)/(RI2*w1)+(CI3*w2)/(RI3*w2)+(CI4*w3)/(RI4*w3)+(C Ip*w4)/(RIp*w4); %计算组合一致性检验 otherwise disp('错误～'); end if CRfinal<0.1 disp('通过组合一致性检验') disp('组合权向量:'); w else disp('不能通过组合一致性检验') end 24