本发明属于核与辐射突发事件后果大尺度模拟领域，具体涉及一种气载放射性核素长距离迁移拉格朗日粒子扩散计算方法，用于对核设施在发生核事故时释放或可能释放到环境中的物质所造成的后果进行预测或评价。
背景技术：
一直以来，世界各国均非常重视和关注放射性核素大尺度大气迁移问题，不仅是军事、国防方面的需要，而且是核事故应急、环境危害问题等的需要。放射性物质在大气中的长距离迁移和后果评价与分析技术涉及众多的科学领域，属于多学科交叉的技术，模拟技术和评价方法也在不断地更新与发展，因此，包括美国、欧盟、日本在内的发达国家、地区也在不断加强该技术的研发。另外，放射性核素大尺度迁移数值模拟技术已经作为一种国家的技术资源，随着相关领域技术的不断发展与更新，从战略角度考虑，一个国家必须有技术储备和可持续发展的能力。当前国际形势变化多端，需要我们不断建立与完善大气输运模拟技术，并提供技术支撑。
近些年来，在朝核试验、日本福岛核电泄漏等事件接连发生的情况下，建立并不断完善一套可应对我国周边数百至数千公里范围内气载放射性物质释放情形应急响应系统的必要性已经成为共识。该系统的建立将为决策人员能够面对类似事件，快速、科学、有效地提出决策建议与方案，最终达到保护我国公众与环境安全的目标。
为此，中国辐射防护研究院自主开发、建立了核与辐射突发事件后果大尺度模拟系统，包括放射性核素大尺度大气迁移数值模拟技术和latmes1.0计算程序，以及相关的模拟结果统计分析程序，可应对全球范围的核与辐射突发事件下，气载放射性物质的迁移、扩散及其环境安全后果的评价。针对朝鲜核试验，根据总装的要求，中辐院定期提供阶段报告预测和分析了朝鲜核试验场址以及宁边核场址放射性核素释放大气迁移对我国的影响。目前已成功用于对日本福岛核事故、朝鲜第三次核试验等事件的有效响应。
境外核爆辐射后果评价系统是通过将区域气象数值预报产品和中长距离放射性污染物迁移模拟和后果评价集成为一体。该系统设计为一个连接多种模型、数据库和地理信息的综合系统，用以模拟、预测和评估指定区域范围内突发事件产生的影响。
技术实现要素：
针对现有技术中存在的问题，本发明提供一种气载放射性核素长距离迁移轨迹拉格朗日粒子扩散计算方法，根据该方法能够估算所有粒子在时间和空间上的总体分布，便于在发生突发事件时，决策人员能够面对该类事件，快速、科学、有效的提出决策建议与对策，从而保护公众与环境的安全。
为达到以上目的，本发明采用的技术方案是：提供一种气载放射性核素长距离迁移拉格朗日粒子扩散计算方法，包括如下步骤：
1)获取天气预报数据、地形数据及模式参数；天气预报数据两个水平风分量u和v、一个垂直风分量w、温度t、湿度的三维场；模式参数包括释放源项参数和设定区域范围；
2)空间坐标和时间坐标变换；
3)调用风场；
4)通过积分粒子运动的拉格朗日方程，计算粒子运行轨迹；xi(t+δt)＝xi(t)+vi(xi(t),t)δt+vi′(xi(t),t)δt+vi″(xi(t),t)δti＝1,2,3
其中，xi为粒子的三维坐标分量；vi为平均风速分量
5)计算放射性衰变和干、湿沉积；
6)统计粒子密度、计算网格浓度；
7)考虑粒子分裂条件，跟踪所有粒子，检查是否有新粒子产生，如果有新粒子产生，返回步骤3)，继续计算；如果没有新粒子产生，结束计算。
进一步，步骤2)中，所述空间坐标和时间坐标变换是将等压层气象数据转换为随地笛卡尔坐标。
进一步，步骤4)中，对运行轨迹中湍流参数和中尺度风脉动参数确定；湍流参数，采用hanna提出的参数化方法，即根据边界层参数混合层高度h、莫宁-奥布霍夫长度l、对流速度尺度w*、粗糙长度z0和摩擦速度u*来计算；中尺度风脉动参数，根据maryon提出的方法，通过假定中尺度风速脉动与hanna参数化方法覆盖的湍流脉动无关，解一个独立的拉格朗日方程来计算，所用到的时间尺度和速度的方差是通过对一个测站的风观测时间序列进行谱分析得到。
进一步，步骤5)中，所述放射性衰变通过下列公式计算：n(1)＝n(0)*0.5δt*λ，
其中，n(1)为衰减后的放射性物质量；n(0)为衰减前的放射性物质量；△t为衰变时间；λ为衰变常数。
进一步，步骤5)中，所述干、湿沉积通过下列公式计算：cd＝c*vd，
其中，cd为沉降后的浓度；c为空气浓度；vd为沉降速度。
进一步，步骤6)中，所述网格浓度通过下列公式计算：
其中，q，释放的放射性活度；n，释放的粒子总数；t，所有粒子扩散时间的总和；nik，第i个粒子在k网格中穿行时间内的时间步数；δt，时间步长；vk，网格k的体积。
本发明的有益技术效果在于：本发明采用粒子随机游走的方法模拟大气扩散，把每个污染质点当成有标志的质点，通过释放大量粒子，计算粒子的轨迹，而这些粒子描述了气载污染物在大气中的迁移扩散；本发明粒子在流场中按平均风输送，同时又用一系列随机位移来模拟湍流扩散，这样就表达了平流和湍流扩散两种作用，对于中长距离迁移问题还考虑了中尺度风脉动的影响，最后由这些粒子在空间和时间上的总体分布估算出污染物的分布，从而保证数据的准确性；由此，在发生突发事件时，决策人员能够面对该类事件，快速、科学、有效的提出决策建议与对策，从而保护公众与环境的安全。
附图说明
图1是本发明气载放射性核素长距离迁移轨迹拉格朗日粒子扩散计算方法的流程图。
具体实施方式
下面结合附图，对本发明的具体实施方式作进一步详细的描述。
如图1所示，是本发明提供的气载放射性核素长距离迁移轨迹拉格朗日粒子扩散计算方法，该方法包括如下步骤：
1)获取天气预报数据、地形数据、下垫面特征数据、模式参数；
天气预报数据包括：(a)空间三维场数据：两个水平风分量u和v，一个垂直风分量w，温度t，特征湿度；(b)空间二维场数据：地表压力，雪厚度，海面压力，云状，10m高u、v风速，2m高温度，2m高露点温度，大尺度降水，地面感应热通量，太阳辐射，地面应力。
模式参数包括：(a)释放源项参数，即爆炸当量、核素种类及放射性活度；(b)模拟区域范围，并对该区域范围进行网格划分；(c)坐标系类型、时间长度、起止时间。
2)空间坐标和时间坐标变换，是通过已知的方法将等压层气象数据转换为随地笛卡尔坐标。
3)调用风场；按照文件给定顺序读取风场文件名称、各风场时刻、以及各时刻各节点气象数据信息。
4)通过积分粒子运动的拉格朗日方程，计算粒子运行轨迹；利用风场插值，获得精确的风场数据。
拉格朗日粒子弥散模式是把每个污染质点当成有标志的质点，通过释放大量粒子，计算粒子的轨迹，而这些粒子描述了气载污染物在大气中的迁移扩散。粒子在流场中按平均风输送，同时又用一系列随机位移来模拟湍流扩散，这样就表达了平流和湍流扩散两种作用，最后由这些粒子在空间和时间上的总体分布估算出污染物的分布。
积分粒子运动的拉格朗日方程，将粒子的运行轨迹写成下列形式：
xi(t+δt)＝xi(t)+vi(xi(t),t)δt+vi′(xi(t),t)δti＝1,2,3(1)
其中，xi为粒子的三维坐标分量；vi为平均风速分量
u′(t+δt)＝u′(t)r+(1-r2)1/2σuξ
v′(t+δt)＝v′(t)r+(1-r2)1/2σvξ(2)
其中，式中右边第二项代表速度涨落量中的随机部分，ξ为符合高斯分布(平均值为0、标准差为σi)的随机数；
方程(1)考虑了平均风输送影响和大气湍流风脉动的影响，湍流脉动反映了时间尺度小于1小时，对应于较短的长度尺度。中尺度运动可以使弥散的烟羽显著增大(gupta等，1997)，对于大尺度模拟问题，需要考虑中尺度风脉动影响。因此，考虑中尺度风脉动影响的粒子运动方程为：
xi(t+δt)＝xi(t)+vi(xi(t),t)δt+vi′(xi(t),t)δt+vi″(xi(t),t)δti＝1,2,3(3)
其中，vi″为中尺度风脉动速度分量(u″,v″,w″)。
本发明运用拉格朗日粒子弥散模拟方法的关键在于确定决定粒子运动平均轨迹的平均风场、三个速度分量的拉格朗日时间尺度和涨落的标准差。
对平均风场，通过解二阶差分形式的轨迹运动方程，
对于湍流参数，采用hanna(1982)提出的一种参数化方法，即根据边界层参数混合层高度h,莫宁-奥布霍夫长度l,对流速度尺度w*,粗糙长度z0和摩擦速度u*来计算湍流参数。由于hanna的方法在整个行星边界层并不总能得到光滑的σw廓线，导致粒子不能充分混合，因此采用ryall和maryon(1997)提出的一种修正方法确定σw。
对于边界层参数(如l，u*)，利用模式第一层高度、地面10m和2m高度的风温数据然后用廓线方法计算上述参数，采用迭代的方法解下列方程：
其中，κ，karman常数；zi，模式第一层高度；δu，模式第一层高度和10m高之间的风速差；δθ，模式第一层高度和2m高之间的位温差；g，重力加速度；θ*，温度尺度；l，莫宁-奥布霍夫长度；
其中，φ1和φ2分别为风速和温度的廓线函数，其形式分别为
对于中尺度风脉动参数，采用maryon(1998)提出的方法：通过假定中尺度风速脉动与hanna参数化方法覆盖的湍流脉动无关，解一个独立的拉格朗日方程来解决，所用到的时间尺度和速度的方差是通过对一个测站的风观测时间序列进行谱分析得到的。假定在网格尺度上观测的风的方差也为亚网格风速标准差提供一些信息，这样，用一种简化方法确定风速标准差的取值方法是：计算粒子所在位置周围16个网格点(在时间和空间上)的速度的标准差，然后取该标准差的一半作为解中尺度拉格朗日方程时所用的风速标准差。
5)采用源耗减的概念来计算干、湿沉积
给定物质的干沉积用沉积速度vd(m/s)来描述，可根据核素类型和下垫面特征来确定。雨水冲洗效应是造成地面高放射性污染水平的最重要的因素之一。湿沉积可以用类似干沉积的方法计算，不同之处仅在于用冲洗系数λ(s-1)代替干沉积速度，冲洗系数的大小取决于降雨强度。放射性衰变的估算：根据公式
n(1)＝n(0)*0.5δt*λ
其中，n(1)为衰减后的放射性物质量；n(0)为衰减前的放射性物质量；△t为衰变时间，s；λ为衰变常数，s-1。
干、湿沉积量的估算。只针对气溶胶进行沉积计算。
cd＝c*vd
cd为沉降后的浓度，bq/m2；
c为空气浓度，bq/m3；
vd为沉降速度，m-1。
当近地层有降水发生时，同时计算干沉积与湿沉积量；
当近地层无降水发生时，只计算干沉积。
6)统计粒子密度、计算网格浓度
为了提高扩散计算精度和效率，当采用粒子分裂技术和核函数方法。利用核函数方法对空间关注点浓度进行计算。核函数方法认为离散化后的烟团在迁移过程中，其自身也呈现高斯扩散的趋势。对本模式而言，由于空间网格间距较大，而单纯的粒子模式如果不使用网格嵌套技术，则网格内部呈现不合理的均匀浓度，并引起监测点浓度的系统误差。所以，在模式中引入核函数方法，这样不仅可以避免复杂的嵌套计算，而且可以快速合理地得到有限关注点的浓度。同时，对于一般意义的核函数方法，不强调单个质子团的空间扩散方向性，但由于大尺度模式中，质子团稀疏，对称分布的核函数仍然会引起不真实的空间浓度分布。对此，模式根据气象条件设定不同参数，使质子团的核函数扩散呈现更为真实的三维不对称分布。
放射性核素空气浓度的计算方法：每个网格中的浓度正比于质点通过该网格时所需时间的总和，因此每个网格的浓度ck(bq/m3)用下列公式计算：
其中，q，释放的放射性活度，bq；n，释放的粒子总数；t，所有粒子扩散时间的总和，s；nik，第i个粒子在k网格中穿行时间内的时间步数；δt，时间步长，s；vk，网格k的体积，m3。在浓度计算中，根据放射性核素的半衰期和粒子迁移的总时间对放射性衰变进行修正。
7)输出时间间隔
如果达到设定的浓度输出时间，计算该时刻网格的浓度，并输出空气浓度场；如果没有达到设定的浓度输出时间，考虑粒子分裂条件，跟踪所有粒子、检查是否有新粒子释放；如果有新粒子释放返回步骤3)，如果没有新粒子释放，结束程序。
其中：粒子分裂条件，是当某粒子所在空间的12个相邻空间均无粒子分布时，对粒子进行分裂。即，分裂后的粒子质量减半，总模拟粒子数增加1。跟踪所有的粒子，检查是否有新粒子释放，是判断模拟时间是否处于释放时间段内，如果处于释放段内，则新增n个粒子的初始迁移位置，并与其它已有粒子一同进行迁移扩散模拟。
对于方程(3)求解，当给定初始条件，通过迭代计算就可以确定各随机游走粒子的运动轨迹。
以下是本发明在计算过程中所涉及到的数据：
(1)坐标系
模式采用混合坐标系，即x、y、η坐标，其中η是一种对气压坐标进行变换的垂直坐标，η与气压坐标的变换方法为：
pk＝ak+bkps
ηk＝ak/p0+bk(8)
其中，ηk为模式第k层的η值；ps为地表压力；p0为压力常数(101325pa)。ak和bk为系数，由最接近地表(随地坐标)的值和最大压力高度层的值确定，中间高度层的系数值则根据随地表层和压力高度层之间的压力梯度确定。
此外，如果纬度超过75°时考虑极地立体投影。
(2)时间、空间分辨率
模式的时间、空间分辨率取决于两方面因素：一是输入资料的时间、空间分辨率；二是模式计算的时间、空间分辨率，在输入资料的时空分辨率一定的情况下，模式计算的时空分辨率主要取决于时间步长。
①输入资料的时间、空间分辨率
在数值模式中，地形和气象数据的空间分辨率是一致的，可接受的空间分辨率从0.25°到2.5°，时间分辨率从3h到12h。
数值计算分析结果表明，随着风场网格分辨率的降低，数值模拟结果的准确程度也随之降低。数值模拟实验表明：风特征大约可以被0.5°(约45km)空间分辨、6小时时间分辨水平的预报风场解析，而利用2.5°(约225km)与12小时的时空分辨预报风场则无法满足需要。因此，可以认为6小时的风场时间间隔是大尺度模式运行的基本条件。综合考虑模拟精度和计算资源与计算速度，在本系统中采用0.5°的空间分辨率、3～6小时的时间分辨率水平完全满足实际应用的需要。
②时间步长
时间步长的确定取决于网格的格距、风速和离散风场的时间间隔，从空间角度遵循：
其中，δxi为空间格距，i代表三个方向；vi分别代表三个方向的风分量。从时间角度遵循：
其中，δt风场为输入风场的时间间隔，即风观测或数值天气预报风场的周期。粒子迁移计算模式的时间步长最大值要小于δt空间和δt时间的最小值，最小时间步长可为1秒。
本发明的气载放射性核素长距离迁移轨迹拉格朗日粒子扩散计算方法并不限于上述具体实施方式，本领域技术人员根据本发明的技术方案得出其他的实施方式，同样属于本发明的技术创新范围。