数学建模竞赛论文-基于日本核辐射扩散污染问题的研究.docx

北京信息科技大学2011年数学建模竞赛论文论文名称： 基于日本核辐射扩散污染问题的研究 学院、班级：理学院 信计0902 信息管理学院 信管0903机电工程学院 车辆0901 参赛组别： b题 组 员： 电 话： 电子邮箱： 提交日期： 2010-5-11 核辐射扩散问题 摘要 本文通过建立数学模型研究了日本核电站泄漏所放出的放射性粉尘对中国地区的影响问题。针对远距离以及风向风速的复杂变化等因素，要想得到精确的放射性粉尘浓度显然不太可能做到。对这种求解动态情况下的连续空间以及连续时间的放射性粉尘浓度，我们基于经典的高斯点源扩散模式而建立了动态高斯核辐射扩散模型，这种模型考虑到了距离、风向、风速等复杂因素，而且还将受损核电站的处理进程和现状考虑在内，使得该模型更加具有可信度以及实用性。在建模以及求解过程中，我们抓住主要矛盾，暂时放弃一些次要矛盾，着重解决了如下几个问题：（1）、利用我们的搜索技术知识在各大网站以及图书资料上查询到充足的数据，用excel对所有数据进行分类汇总、运算，筛选出所有有用信息。包括各大相关城市的经纬度、风向、风速以及扩散系数等，并用matlab作出相关图形以便我们进行直观分析；（2）、从整体思想出发，通过层层推进建立动态高斯核辐射扩散模型并编写c#语言程序计算出最终的放射性粉尘浓度；（3）、对模型进行检验，并改变时间以及高空距离等变量以计算出多组不同的数据值，从横向以及纵向比较分析；（4）、为了把该模型应用到更加广泛的领域，我们最后对模型进行了必要而恰当的改进分析。 在对该问题建立模型时，我们考虑了问题的实际性和模型的适用性。我们的模型能够做到有理可依，有据可查，实用而不失创新。关键字：放射性粉尘、扩散系数、风向、风速、源强、动态高斯核辐射扩散模型 一、问题重述 b题：核辐射扩散问题2011年3月11日, 日本近海发生9.0级地震并引发了大海啸，沿海的核电站受到破坏，开始释放出具有放射性物质。很多人担心这些放射性物质会危害自己的健康，因此急切希望了解：地震中损坏的日本核电站散发出的放射性物质，究竟会在什么时候到达自己身边，以及什么时候会达到对人体有害的程度。专家们认为，对日本之外的国家和地区而言，会随空气移动的放射性粉尘可能是主要的威胁。若对此进行预测，需要考虑到风向、风速以及距离受损核电站的距离的远近。截止到2011年3月30日，在我国上海、天津、重庆、河北、山西、内蒙古、吉林、黑龙江、江苏、安徽、浙江、福建、河南、广东、广西、四川、陕西、宁夏部分地区空气中监测到来自日本核事故释放出的极微量人工放射性核素碘131.请在上述背景下,请尝试建立数学模型，解决如下问题:(1) 评估来自日本的放射性粉尘是否会到将覆盖到我国全国范围，以及是否会危害到我国居民的人身安全；(2) 如果会达到危害人体的程度，请预测会在什么时候达到；(3) 如果不会达到危害人体的程度，请估计我国哪个地区受到的核污染是最重的, 并尝试给出在什么时候达到最严重程度。注1：建立模型时，请注意将受损核电站的处理进程和现状考虑在内。注2：根据自己掌握的知识和技术，采用数学模型或者计算机仿真方法均可。 二、符号说明放射性粉尘浓度；离核电站爆炸泄漏的时间；碘-131的半衰期【8.3天】；有效高度；城市i的高空距离；水平扩散系数；垂直扩散系数；扩散系数方程中的参数；地区的风速值；下风向距离；横截风向距离；福岛核电站源强；球面上两点的弦长；球面上两点的弧长；城市i的经度；城市i的纬度；旋转角度；地球半径； 三、模型假设1、从福岛核辐射扩散点源到目的地的距离假设为直线距离；2、全部高度风速均匀稳定；3、在某一个地区假设风速固定不变；4、放射性粉尘浓度量等于来自点源连续不断的量加上放射性元素衰变之后剩余量；5、由于核电站处理进程以及现状难以定量估摸，所以我们假设源强q符合关于时间t的标准正体分布；6、 污染物浓度在y、z风向上分布为正态分布；7、建立类似于高斯点源扩散模式的动态核辐射扩散模型进行求解连续空间的放射性粉尘浓度； 四、问题分析与模型建立4.1问题分析 首先，我们通过查询资料统计出中国的31个城市以及日本福岛的经纬度【见附件1】，并利用matlab作出二维平面概图如下： 从图中可以看到，我们主要选择的城市均离散分布于中国的大江南北，能够整体性地代表中国地区。同时，日本福岛离我们所选择的中国城市的距离较远并且各异，因此，若要考虑到中国地区所受放射性粉尘的威胁，我们有必要考虑到其与受损核电站的距离。 其次，由于放射性粉尘主要是由于空气移动所造成的无边界扩散，显然，风也是不得不考虑的因素之一。 此外，放射性粉尘中含有的主要对人体有害的物质乃是碘-131，它是一种人工放射性元素，那么这就涉及到衰变周期，因此也可以将其纳入我们的考虑因素之中。 最后，由于日本在核电站发生泄漏后采取了一定的措施，也就使得源强在不断改变，所以我们也有必要考虑处理进程及现状。 综上，针对该问题，我们决定在建模过程中考虑到距离、风向、风速、衰变期、源强的变化等因素。4.2模型一：通过经纬度计算距离4.2.1模型的建立与推导各个城市的经纬度值【见附件1】已据实统计.已知地球上两地经纬度计算两地距离公式的推导如下：设有空间直角坐标系o-xyz，为空间中一点，它在xoy面上的射影为，ox轴到的夹角为，与z轴正方向的夹角为，这时p的位置可用数组确定，数组称为点p的球面坐标。显然，点p的空间直角坐标与球面坐标间的关系为： 设地球上任一点的经度为，纬度为，地球半径为r，则的直角坐标为：设为地球上的任意两地，则a、b的直角坐标分别为：a： b： 那么 将代入并化简得： 下面求，即ab的球面距离设,则球面上任意两地的距离为： 4.3模型二：扩散参数的选定4.3.1模型的估算以及确定-【扩散参数的确定pg曲线法】 扩散系数包括水平扩散系数和垂直扩散系数，其值与大气稳定度和风速密切相关。而大气稳定度又与太阳辐射等级相关，通过顺层次递进分析，我们查询出实际的资料并统计【表1、表2、表3】如下： 表1 太阳辐射等级总/低云量夜间h1515h3535h654/4-2-1+1+2+35-7/8/5/5-70000+18/800000 表2 大气稳定度等级地面风速（m/s） 大气辐射等级+3 +2 +1 0 -1 -26dddddd 表3 风速1m/s的扩散系数大气稳定度abcda0.240.930.240.93a-b0.230.920.230.92b0.230.910.230.91b-c0.350.910.2350.855c0.470.910.240.80c-d0.380.780.420.85d0.490.760.130.76d-e0.370.780.100.78e0.250.800.070.80e-f0.340.720.410.57f0.170.800.410.500.76 通过分析比较整个中国大部分地区的夜间时间长短、云量的多少以及中国3-6月份平均风速【见附件2】等，我们最终确定了中国大部分地区大气稳定度可以取为d级，故扩散系数公式中的参数取为：a=0.49,b=0.76,c=0.13,d=0.76扩散系数公式为： 故有： 4.4模型三：通过旋转变换计算下风向距离与横截风向距离4.4.1模型的建立与推导 所谓下风向距离即是两地的距离沿着风向的分距离，而横截风向距离乃是两地距离沿着垂直于风向的分距离。以下利用经纬度通过旋转变换来推导并确定下风向距离与横截风向距离。设为福岛的球面坐标，为地球上任意一点，则与a同经度、与b同纬度的点的球面坐标为以福岛为中心平行于纬线所在直线为x轴，平行于经线所在直线为y轴建立平面直角坐标系a-xy.设b点在a-xy坐标系下的坐标为（x,y），于是对应的坐标分别为 : : : 则有： 将-代入并化简可得： 那么b在a-xy坐标系下的坐标为: 【说明：其中的正负号由b点的经纬度与a点的经纬度的大小关系决定：若b的经度大于a的经度，则x取正，反之取负；若b点的纬度大于a的纬度，则y取正，反之取负】下面以a为原点，风的方向【通过查询资料我们选定为西南方向】为x轴，横截风向为y轴建立平面直角坐标系.设b在坐标系下的坐标为（）则 下风向距离： 横截风向距离为：坐标系可看作原坐标系绕a点旋转所得假设坐标系o-xyz（旧系）旋转到o-（新系）它们的坐标向量分别为和，设新系的坐标向量是，在旧系中的坐标为： m是空间中的任一点，且在新系和旧系中的坐标分别为即 将代入并与式比较可得： 或 且满足正交条件： ， 和 ， 设旋转角为，使用xoy上的旋转变换可得： 和 以上旋转坐标变换公式应用到坐标系a-xy和可得： 将代入最终得出下风向距离与横截风向距离分别为： 4.5模型四：动态高斯核辐射扩散模型4.5.1模型的引用及变形 复杂的大气扩散模式很多，这些模式成为支持环境空气质量评价、大气环境容量核定、污染物浓度等的基本工具。同时，现有扩散模式不断向大型化和复杂化发展，我们所了解到的高斯点源扩散模式也正在向各个方向进行变形发展，以至于能够解决更加实际复杂的问题，在此基础上，我们分析了日本核辐射扩散问题，发现该模式有一定的利用价值；经典的高斯点源扩散模式公式为：它是针对于固定污染物的浓度计算模式，而我们如今要解决的是动态的、衰变的、源强不断变化的核辐射扩散问题，故我们对其进行合理的变形：（1） 由于核电站从发生泄漏之后一直都在进行必要的处理措施，使得当地的核辐射源强也在不断变化，但是要想根据实际数据定量地找出其与时间的函数关系，似乎较为困难，于是根据我们的假设有：源强q是关于时间t的标准正态分布函数关系：；（2） 作为放射性粉尘中对人体有害的碘-131乃是一种人工放射性元素，其随着时间的增加必然有规律地进行衰减。通过物理学知识了解到：放射性元素衰变公式为： 放射性元素的密度为： 通过上述两个公式，我们可以推导出： 而在连续的时间内，根据假设：放射性粉尘浓度量等于来自点源连续不断的量加上放射性元素衰变之后剩余量。于是利用高等数学中的积分可得：综上（1）（2），我们可以得出最终的动态高斯核辐射扩散模型公式为： 此为我们计算放射性粉尘浓度的最终函数关系式，其中的多项变量在前几个模型中已建立了相关函数关系式以待直接使用。 五、模型检验与分析我们建立的最终模型乃是基于经典的高斯点源扩散模式。通过有理有据地分析与测试，最终对其进行向着我们所需要研究核辐射扩散问题的方向有效地变形，使其适用于研究我们的问题。 5.1、城市的选择 日本核辐射属于无边界扩散，作为日本近邻的中国，幅员辽阔，其受扩散影响不容忽视，但是为了定性定量地分析解决所提出的问题，我们必须有的放矢地选择部分城市作为研究对象。针对于扩散的连续性，我们选择的城市属于离散型地区，尽管如此，该方案还是比较可行，具有一定的代表性、合理性、整体性。5.2、扩散系数的估算 核辐射既为扩散性问题，那么就少不了扩散系数的分析以及选择问题，针对整个中国地区来说，如果我们根据每一个小地区的气象问题来查询其扩散系数并带入最终模型中进行计算，不仅复杂，而且量大，根本没有此必要。所以，我们选择了数学中的模糊计算方法通过总体数据估算出了一个整体扩散系数。5.3、多项数据的比较与分析 毋庸置疑，一个地区的放射性粉尘浓度与其高空距离也是有着不容忽视的联系，所以我们在进行模型求解时会选择不同的高空距离进行代入求解，并进行比较分析。 由于本问题研究的不仅仅是一般物质的扩散现象，而是一种人工放射性物质碘131，所以我们在模型建立以及模型求解时得考虑到其随着时间的推移会进行有规律地衰变以及累加。由此，我们在模型求解时会选择多组不同的时间值带入模型公式进行计算比较。 六、模型求解 根据我们所建立的动态高斯核辐射扩散模型编写c#语言程序【见附件3】对其进行求解并输出多组结果【见附件4】，并且利用matlab作出浓度变化曲线。例如以下几个城市： 黑龙江 北京 重庆 山东6.1问题一的求解 通过对模型求解出来的多组数据进行比较分析，我们可以得出：来自日本的放射性粉尘会覆盖到我国全部地区。相对于碘-131危害人体健康时的浓度标准【国家标准(gb18871-2002)规定限值为24.3贝克/立方米】，我们通过数学模型所得出的浓度值远远小于此标准，故不会危害到我国居民的人身安全。6.2问题三的求解 既然不会危害到我国居民的人身安全，不管是纵向比较还是横向比较，我们都可以通过数据得出，相对于其它所选择的研究对象城市，沿海地区城市的放射性粉尘浓度值始终保持较高，故我们可以估计出沿海地区受到的核辐射较为严重，例如：上海、福建、广东、山东、天津等；而这些地区在福岛核电站发生泄漏之后的大概20多天左右，即4月上旬这段时间。 七、模型评价7.1模型的合理性 为了回答日本核辐射是否会覆盖到整个中国地区，我们将放射性粉尘浓度作为标准来进行主要判断，但是由于中国幅员辽阔，地形复杂，天气变化难以捉摸，要想精确地得到每个地区的放射性粉尘浓度，显然不太可能，而且计算量极大，根本没有必要，对于我们来说，也没有那么多的精力。所以我们所建立的动态高斯核辐射扩散模型主要针对于中国地区的具有代表性质的31个城市进行浓度计算分析，比较合理，且相对于其它方法计算较易实现。在此模型中，我们考虑了主要的影响因素，包括距离、风向、风速、衰变期等，比较全面而且合理。7.2模型的创新性由于受到各种因素的影响，核辐射扩散本是比较难以捉摸的，但是我们秉承在解决实际问题中抓住主要矛盾进行分析的方法，利用我们所学的数学知识以及前人在类似问题上的研究结果基础上，进行创新建立了“动态高斯核辐射扩散模型”，以此能够解决该日本核辐射扩散问题，特别是应用于像我们这种邻近核辐射源头的国家地区。该模型具有一定的创新性。7.3根据模型所求出数据的合理性 我们将中国整体地区划分为各个局部地区，利用该模型以小求大，我们通过编写c#程序对其进行求解，与此同时，为了得到的数据更加具有说服力，更加具有代表性，更加具有合理性，我们选择了不同的时间、不同的高空距离进行不同的浓度计算，不仅从纵向而且从横向进行对比分析，然后得出比较合理而真实的结果，以便解答我们所需要回答的问题。八、模型改进及应用该模型中的变量因素距离、风向、风速等尚欠更加具体的数据，由于日本福岛离中国地区距离较远，风向也每日不同，风速也并不稳定，这些因素复杂多变，虽然我们的假设具有一定的合理性，所建立的模型也具有一定的创新性及实用性，但是由于空间的广度性、因素的复杂性、时间的有限性，我们的模型难免会有一些不足之处，得出的数据虽然足以解决我们所要解决的问题。不得不承认，要将该模型推广到其它对象作为实际应用，尚需改进。因此我们认为改进的地方有以下几点：（1） 对中国地形的相关数据进行实际分析归类，将我们所选择的31个城市推广并深入到更加具有代表性的不同地区，并不以大城市为主要的对象，普遍性、代表性还需深入研究。（2） 中国各个地区的大量而具体的气象数据还需认真核对统计并进行研究分类，以确立更加精确的风向与风速。（3） 该模型是基于经典的高斯点源扩散模式建立起来的，涉及到的高等数学知识以及其它知识较多，所以如果可能的话，进行必要的简化改进，以使更加多的读者能够深入理解。就中国地区而言，可以利用其对中国受核辐射问题进行研究。该模型具有一定的普遍适用性，所以可以将其推广到日本的所有邻近国家进行类似研究。凡是涉及到放射性粉尘扩散问题，我们均可以利用该动态高斯核辐射扩散模型进行预测。因此，我们的模型具有足够的实用性以及应用性。九、参考资料1、 环境影响评价技术导则 hj/t2.1-2.3-93；2、 e.a.bender. 数学模型引论 朱尧辰、徐伟宣译 科学普及出版社19823、 基于mpi的大气污染扩散模型的并行计算研究 （华东师范大学地理信息科学教育部开放实验室，上海，200062 姚丽萍 王远飞）；4、 王凤林、毕彤、刘振山 突发性污染事故大气扩散数学模型初探（沈阳市环境监测中心站 沈阳11015）； 5、百度百科知识；6、中国气象局网站、日本气象厅网站、中国环境保护部网站、中国气象科学数据共享服务网；十、附录附录1主要城市经纬度城市经度纬度日本福岛141.03 37.04 北京116.46 39.92 上海121.48 31.22 天津117.20 39.13 重庆106.54 29.59 昆明102.73 25.04 呼和浩特111.65 40.82 长春125.35 43.88 成都104.06 30.67 银川106.27 38.47 合肥117.27 31.86 济南117.00 36.65 太原112.53 37.87 广州113.23 23.16 南宁108.33 22.84 乌鲁木齐87.68 43.77 南京118.78 32.04 南昌115.89 28.68 石家庄114.48 38.03 郑州113.65 34.76 杭州120.19 30.26 海口110.35 20.02 武汉114.31 30.52 长沙113.00 28.21 兰州103.73 36.03 福州119.30 26.08 拉萨91.11 29.97 贵阳106.71 26.57 沈阳123.38 41.80 西安108.95 34.27 西宁101.74 36.56 哈尔滨126.63 45.75 附件21971-2000年中国城市平均风速值（0.1m/s）城市三月四月五月六月平均风速北京31 32 28 24 28.75 上海33 32 33 32 32.50 天津29 33 30 25 29.25 重庆16 15 14 14 14.75 昆明29 27 25 21 25.50 呼和浩特22 36 25 22 26.25 长春44 52 48 37 45.25 成都33 33 32 28 31.50 银川24 28 26 23 25.25 合肥34 31 29 28 30.50 济南35 38 36 32 35.25 太原21 29 26 21 24.25 广州16 16 17 17 16.50 南宁30 24 21 17 23.00 乌鲁木齐22 29 30 28 27.25 南京29 27 26 26 27.00 南昌26 24 22 21 23.25 石家庄14 17 15 12 14.50 郑州30 30 27 25 28.00 杭州24 24 23 23 23.50 海口27 27 24 24 25.50 武汉20 20 19 18 19.25 长沙24 23 21 20 22.00 兰州25 26 24 20 23.75 福州22 24 23 26 23.75 拉萨24 24 22 20 22.50 贵阳24 23 22 20 22.25 沈阳32 36 35 29 33.00 西安18 18 18 19 18.25 西宁20 21 19 16 19.00 哈尔滨40 47 45 36 42.00 附件3c#语言程序代码：using system;using system.collections.generic;using system.linq;using system.text;namespace 动态高斯核辐射扩散模型公式求解 class program static void main(string args) city head = new city(福岛, 141.03, 37.04, 8.3); /福岛作为污染源头，设为距离起点 city china_city = new city31; china_city0 = new city(北京, 116.46, 39.92, 2.875); china_city1 = new city(上海, 121.48, 31.22, 3.250); china_city2 = new city(天津, 117.20, 39.13, 2.925); china_city3 = new city(重庆, 106.54, 29.59, 1.475); china_city4 = new city(昆明, 102.73, 25.04, 2.550); china_city5 = new city(呼和浩特, 111.65, 40.82, 2.625); china_city6 = new city(长春, 125.35, 43.88, 4.525); china_city7 = new city(成都, 104.06, 30.67, 3.150); china_city8 = new city(银川, 106.87, 38.47, 2.525); china_city9 = new city(合肥, 117.27, 31.86, 3.050); china_city10 = new city(济南, 117.00, 36.65, 3.525); china_city11 = new city(太原, 112.53, 37.87, 2.425); china_city12 = new city(广州, 113.23, 23.16, 1.650); china_city13 = new city(南宁, 108.33, 22.84, 2.300); china_city14 = new city(乌鲁木齐, 87.68, 43.77, 2.725); china_city15 = new city(南京, 118.78, 32.04, 2.700); china_city16 = new city(南昌, 115.89, 28.68, 2.325); china_city17 = new city(石家庄, 114.48, 38.03, 1.450); china_city18 = new city(郑州, 113.65, 34.76, 2.800); china_city19 = new city(杭州, 120.09, 30.26, 2.350); china_city20 = new city(海口, 110.35, 20.02, 2.550); china_city21 = new city(武汉, 114.31, 30.52, 1.925); china_city22 = new city(长沙, 113.00, 28.21, 2.200); china_city23 = new city(兰州, 103.73, 36.03, 2.375); china_city24 = new city(福州, 119.30, 26.08, 2.375); china_city25 = new city(拉萨, 91.11, 29.97, 2.250); china_city26 = new city(贵阳, 106.71, 26.57, 2.225); china_city27 = new city(沈阳, 123.38, 41.80, 3.300); china_city28 = new city(西安, 108.95, 34.27, 1.825); china_city29 = new city(西宁, 101.74, 36.51, 1.900); china_city30 = new city(哈尔滨, 126.63, 45.75, 4.200); double distance = new double31; for (int i = 0; i 31; i+) console.write(福岛到0的距离为：, china_cityi.name); distancei = city.distance(head, china_cityi); console.writeline(distancei); for (int i = 0; i = b.latitude) return math.abs(-r * math.cos(b.latitude * math.pi / 180) * math.sqrt(2 - 2 * math.cos(a.longitude * math.pi / 180 - b.longitude * math.pi / 180) * math.cos(135 * math.pi / 180) - r * math.sqrt(2 - 2 * math.cos(a.latitude * math.pi / 180 - b.latitude * math.pi / 180) * math.sin(135 * math.pi / 180); else return math.abs(-r * math.cos(b.latitude * math.pi / 180) * math.sqrt(2 - 2 * math.cos(a.longitude * math.pi / 180 - b.longitude * math.pi / 180) * math.cos(135 * math.pi / 180) + r * math.sqrt(2 - 2 * math.cos(a.latitude * math.pi / 180 - b.latitude * math.pi / 180) * math.sin(135 * math.pi / 180); public double juli_y(city a, city b) /返回横截风向距离 if (a.latitude = b.latitude) return math.abs(r * math.cos(b.latitude * math.pi / 180) * math.sqrt(2 - 2 * math.cos(a.longitude * math.pi / 180 - b.longitude * math.pi / 180) * math.sin(135 * math.pi / 180) - r * math.sqrt(2 - 2 * math.cos(a.latitude * math.pi / 180 - b.latitude * math.pi / 180) * math.cos(135 * math.pi / 180); else return math.abs(r * math.cos(b.latitude * math.pi / 180) * math.sqrt(2 - 2 * math.cos(a.longitude * math.pi / 180 - b.longitude * math.pi / 180) * math.sin(135 * math.pi / 180) + r * math.sqrt(2 - 2 * math.cos(a.latitude * math.pi / 180 - b.latitude * math.pi / 180) * math.cos(135 * math.pi / 180); class diffusivity /扩散系数 public double diffusivity_y(double x) return 0.49 * math.pow(x, 0.76); public double diffusivity_z(double x) return 0.13 * math.pow(x, 0.76); public class city private string name; /城市名称 private double longitude; /城市经度 private double latitude; /城市纬度 private double wind_speed; /城市平均风速 public const double r = 6378.040; /地球半径 public const double t = 8.3; /碘131半衰期 public string name get return name; set name = value; public double longitude get return longitude; set longitude = value; public double latitude get return latitude; set latitude = value; public double wind_speed get return wind_speed; set wind_speed = value; public city() public city(string name, double longitude, double latitu