, 苏国锋 清华大学 工程物理系, 公共安全研究院, 北京 100084
收稿日期：2018-05-23
基金项目：国家重点研发计划(2016YFF0103901)
作者简介：黎岢(1989-), 男, 博士研究生
通信作者：梁漫春, 副研究员, E-mail:lmc@tsinghua.edu.cn
摘要：发生核事故后，可以通过模型迅速预测泄漏的放射性物质的大气扩散情况，但释放源项、气象等参数的不确定性导致大气扩散模型计算的准确性受到影响。通常可以采用数据同化来改善模型预测结果。该文提出一种基于Gauss烟团模型的大气扩散数据同化方法，可以结合观测数据改善模型的预测结果。该方法在迭代搜索中对烟团参数进行线性变换以简化Gauss烟团模型，利用粒子群优化算法对泄漏速率、释放高度、风向、平均风速这4个模型参数进行校正。该方法适用于平坦地形和均匀稳定气象条件下的中尺度扩散。采用稳态条件下双生子实验进行验证，观测点位的模拟值与真实观测值的相关系数达到0.99。在非稳态条件下对源项估计结果进行了测试，结果略优于集合Kalman滤波方法，相关系数达到0.68。该同化方法计算速度快，能有效提升模型的预测，可用于大气扩散数据同化。
关键词：放射性核事故大气扩散数据同化粒子群优化
Data assimilation method for atmospheric dispersion based on a Gaussian puff model
LI Ke, LIANG Manchun
, SU GuofengInstitute of Public Safety Research, Department of Engineering Physics, Tsinghua University, Beijing 100084, China
Abstract: Models are needed to quickly predict the atmospheric dispersion of radioactive material released in a nuclear accident. However, the uncertainties in the source term, meteorological data, and other conditions reduce the dispersion model prediction reliability. Data assimilation (DA) is usually introduced to improve the model predictions. The paper presents a DA method based on a Gaussian puff model to improve the predictions using some observed data. The method modifies the puff parameters to approximate the observed data in an iterative search. The four model parameters modified using particle swarm optimization in this study are the release rate, release height, wind direction, and mean wind speed. The method is applicable to mesoscale atmospheric dispersion with uniform and stable conditions over a flat area. Twin experiments are used to verify this DA method. The correlation coefficient between the experimental group and the control group at the observation points is 0.99. The source estimation in the non-steady condition is also tested with the correlation coefficient of 0.68, slightly better than the ensemble Kalman filter method. The method converges rapidly with good model predictions; thus, this method is useful for data assimilation of atmospheric dispersion.
Key words: radioactive nuclear accidentatmospheric dispersiondata assimilationparticle swarm optimization
大气扩散模型是对污染物在大气中弥散状况的模拟，是核事故后果评价的重要工具。它可以对泄漏的有害物质的分布进行预测，从而为应急决策提供技术支持。现今一些研究机构已经开发了很多较为成熟的大气扩散模型，并应用于核事故应急决策支持系统中，如ARAC/NARAC[1]、RIMPUFF[2]、CALPUFF[3]等。然而，由于输入模型的释放源项、气象数据以及模型本身等存在不确定性[4]，模型计算结果可能与实际情况存在很大差距，不利于做出正确的防护行动决策。因此，提高模型预测结果的准确性、减小模型输入参数的不确定性对模型预测结果的影响是当前研究的热点。
引入观测数据对大气扩散模型进行校准的方法被称为数据同化。数据同化可以将不同来源、不同误差的观测数据与模型的计算结果相结合，从而使模型的计算结果接近实际观测值[5]。目前，在大气扩散模拟的数据同化方法研究中，较多采用集合Kalman滤波(ensemble Kalman filter, EnKF)算法、变分方法等。
欧洲核应急决策支持系统(realtime online decision support system, RODOS[6])基于Kalman滤波算法对大气扩散、沉降、食物链及水文模型开发了数据同化算法[7]。一些****也基于EnKF算法开展了很多大气扩散数据同化的研究[8-11]。一些****还基于Gauss烟羽模型与Gauss烟团模型，应用变分方法[12-14]或者粒子滤波算法[15-16]来改善模型参数，以及进行源项反演。
本文提出一种基于Gauss烟团模型的大气扩散数据同化方法，优化大气扩散模型的计算结果。该方法计算速度快，可以适用于中尺度(小于50 km)平坦地形条件下稳态的大气扩散数据同化。本文采用双生子模拟实验对数据同化方法进行验证，并与EnKF算法进行了比较。
1 研究方法大气扩散系统的状态与观测值可以表述为：
| ${\mathit{\boldsymbol{X}}_{k + 1}} = M\left( {{\mathit{\boldsymbol{X}}_k}} \right),$ | (1) |
| ${\mathit{\boldsymbol{Y}}_{k + 1}} = H\left( {{\mathit{\boldsymbol{X}}_{k + 1}}} \right).$ | (2) |
本文所使用的大气扩散模型是Gauss烟团模型，这里的状态量包括了所有烟团的水平位置、高度以及所包含的放射性物质的活度。
1.1 大气扩散模型Gauss烟团模型被广泛用于核事故后果评价系统中。事故释放出的污染物被模拟成多个Gauss烟团，每个烟团含有一定质量的污染物，烟团在大气中输运与扩散。每个烟团对某一位置的污染物活度浓度贡献为
| $\begin{array}{*{20}{c}}{{C_i}\left( {{x_{\rm{g}}},{y_{\rm{g}}},{z_{\rm{g}}}} \right) = \frac{{A\left( i \right)}}{{{{\left( {2\pi } \right)}^{3/2}}{\sigma _x}\left( i \right){\sigma _y}\left( i \right){\sigma _z}\left( i \right)}} \cdot }\\{\exp \left[ { - \frac{{{{\left( {{x_{\rm{g}}} - {x_{\rm{c}}}\left( i \right)} \right)}^2}}}{{2\sigma _x^2\left( i \right)}} - \frac{{{{\left( {{y_{\rm{g}}} - {y_{\rm{c}}}\left( i \right)} \right)}^2}}}{{2\sigma _y^2\left( i \right)}}} \right] \cdot }\\{\left\{ {\exp \left[ { - \frac{{{{\left( {{{\rm{z}}_{\rm{g}}} - {z_{\rm{c}}}\left( i \right)} \right)}^2}}}{{2\sigma _z^2\left( i \right)}}} \right] + \exp \left[ { - \frac{{{{\left( {{z_{\rm{g}}} + {z_{\rm{c}}}\left( i \right)} \right)}^2}}}{{2\sigma _z^2\left( i \right)}}} \right] + } \right.}\\{\left. {\exp \left[ { - \frac{{{{\left( {2{z_{{\rm{inv}}}} - {z_{\rm{c}}}\left( i \right)} \right)}^2}}}{{2\sigma _z^2\left( i \right)}}} \right]} \right\}.}\end{array}$ | (3) |
对所有烟团的污染物活度浓度贡献进行叠加，即可得到(xg, yg, zg)位置的污染物活度浓度值，
| $C\left( {{x_{\rm{g}}},{y_{\rm{g}}},{z_{\rm{g}}}} \right) = \sum\limits_i {{C_i}\left( {{x_{\rm{g}}},{y_{\rm{g}}},{z_{\rm{g}}}} \right)} .$ | (4) |
表 1 开阔平原田野地区的Briggs经验公式[4]
| 稳定度 | σy/m | σz/m |
| A | 0.20x | |
| B | 0.12x | |
| C | ||
| D | ||
| E | 0.03x/(1+0.003x) | |
| F | 0.016x/(1+0.003x) | |
| 注：x表示放射性物质从释放开始输运的距离。 | ||
表选项
第i个烟团在大气中的输运过程满足
| ${\mathit{\boldsymbol{X}}_{i,k + 1}} = {\mathit{\boldsymbol{X}}_{i,k}} + {\mathit{\boldsymbol{V}}_{i,k}}\Delta t.$ | (5) |
1.2 基于Gauss烟团模型的数据同化方法大气扩散模型实际是对污染物在大气中弥散现象的模拟，但往往由于模型本身对物理过程的简化以及模型参数的不确定性(如源项、气象数据)等原因，导致预测结果与实际情况不符。数据同化方法将获取的大气中放射性物质活度浓度的监测数据与模型的预测结果相结合，找到最优的模型输入参数，并计算得到更符合真实观测值的模型预测结果。
常用的数据同化方法普遍存在计算耗时过长的问题，如EnKF算法在每一时间步同化时，通常需要进行数十次的大气扩散模型计算。本文基于Gauss烟团模型，提出一种快速简便的大气扩散数据同化方法，通过对烟团位置的仿射变换，来模拟扩散模型计算结果的修正，从而避免多次代入扩散模型计算过程。
一般认为源项与气象数据是大气扩散模型误差的主要来源，因此本文中考虑的影响模型准确性的参数为释放速率、释放高度、风向与平均风速。
1.2.1 基本原理在Gauss烟团模型中，活度浓度由各个烟团的贡献叠加计算。当模型参数存在不确定性时，烟团的状态(烟团的位置、烟团中的污染物含量)也会受到影响。
| ${\mathit{\boldsymbol{X}}_i} = M\left( {Q,H,\phi ,v,{t_i}} \right).$ | (6) |
| ${\mathit{\boldsymbol{X}}_{i,{\rm{ dist}}}} = M\left( {{Q_{{\rm{dist}}}},{H_{{\rm{dist}}}},{\phi _{{\rm{dist}}}},{v_{{\rm{dist}}}},{t_i}} \right).$ | (7) |
该状态对应的模拟观测数据可由状态量计算得到，
| ${\mathit{\boldsymbol{Y}}_{{\rm{dist}}}} = H\left( {{\mathit{\boldsymbol{X}}_{{\rm{dist}}}}} \right).$ | (8) |
| $J = F\left( {{\mathit{\boldsymbol{Y}}_{{\rm{dist}}}},{\mathit{\boldsymbol{Y}}_{\rm{o}}}} \right).$ | (9) |
模型参数会对模型状态量产生影响。释放高度的不确定性会影响烟团的初始高度，释放速率的不确定性会影响烟团所携带污染物的质量，风向与风速的不确定性则影响烟团输运的方向与距离。根据Gauss烟团模型的原理与计算公式，模型参数对烟团状态{Xi}(xi, yi, zi, Ai)的影响可近似归纳如下：
当释放高度参数与实际释放高度存在偏差ΔH时，即Hdist=Htrue+ΔH，烟团的初始高度会在实际高度的基础上有一个增量(不考虑烟羽抬升高度的影响)。下标true代表准确的参数或状态量。
| ${z_{i{\rm{, dist}}}} = {z_{i{\rm{,true}}}} + \Delta H.$ | (10) |
| ${A_{i{\rm{, dist}}}} = \gamma {A_{i{\rm{,true}}}}.$ | (11) |
| $\left[ {\begin{array}{*{20}{c}}{{d_{x{\rm{, dist}}}}}\\{{d_{y{\rm{, dist}}}}}\\1\end{array}} \right] = \left[ {\begin{array}{*{20}{c}}{\cos \Delta \phi }&{ - \sin \Delta \phi }&0\\{\sin \Delta \phi }&{\cos \Delta \phi }&0\\0&0&1\end{array}} \right] \cdot \left[ {\begin{array}{*{20}{c}}{{d_{x{\rm{, true}}}}}\\{{d_{y{\rm{, true}}}}}\\1\end{array}} \right].$ | (12) |
| $\left[ {\begin{array}{*{20}{c}}{{d_{x{\rm{, dist}}}}}\\{{d_{y{\rm{, dist}}}}}\\1\end{array}} \right] = \left[ {\begin{array}{*{20}{c}}\lambda &0&0\\0&\lambda &0\\0&0&1\end{array}} \right] \cdot \left[ {\begin{array}{*{20}{c}}{{d_{x{\rm{, true}}}}}\\{{d_{y{\rm{, true}}}}}\\1\end{array}} \right],$ | (13) |
| ${x_{{\rm{dist}}}} = {x_{{\rm{true}}}} + {d_{x{\rm{, dist}}}},$ | (14) |
| ${y_{{\rm{dist}}}} = {y_{{\rm{true}}}} + {d_{y{\rm{, dist}}}}.$ | (15) |
1.2.2 适应度与代价函数模型参数的最优估计与适应度函数有关，这里采用观测数据与模拟数据的相关系数来确定, 相关系数越大，则适应度越好。
| ${\rm{Corr = }}\frac{{{\rm{Cov}}\left( {{\mathit{\boldsymbol{Y}}_{{\rm{dist}}}},{\mathit{\boldsymbol{Y}}_{\rm{o}}}} \right)}}{{\sigma \left( {{\mathit{\boldsymbol{Y}}_{{\rm{dist}}}}} \right) \cdot \sigma \left( {{\mathit{\boldsymbol{Y}}_{\rm{o}}}} \right)}}.$ | (16) |
代价函数取相关系数加1的倒数，
| $J = 1/\left( {1 + {\rm{Corr}}} \right).$ | (17) |
1.2.3 搜索算法数据同化的重要问题是寻找最优参数，从而最小化代价函数，这里需要优化的参数包括泄漏速率、释放高度、风向与平均风速。因为适应度是相关系数的函数，泄漏速率并不影响适应度，所以本文对泄漏速率的优化是在对其他3个参数优化之后，利用最小二乘法确定。因此，优先考虑对释放高度、风向与平均风速这3个参数的优化。
寻找最优参数可以采取遍历方法，但这种方法在精度要求高且搜索范围较大时十分耗费时间，因此本文采用粒子群算法寻找最优参数。
粒子群算法模拟鸟类群体的觅食行为，将鸟类飞行的范围比作优化问题的解空间。每只鸟被抽象为一个粒子，代表问题的一个可能解。通过个体之间的协作与竞争实现多维空间内最优解的搜索[18]。
在D维搜索空间中，由N个粒子组成粒子群，每个粒子表示为空间中的一个向量：
| ${\mathit{\boldsymbol{X}}_i} = \left( {{x_{i1}},{x_{i2}}, \cdots ,{x_{iD}}} \right),i = 1,2, \cdots ,N.$ | (18) |
| ${\mathit{\boldsymbol{V}}_i} = \left( {{v_{i1}},{v_{i2}}, \cdots ,{v_{iD}}} \right),i = 1,2, \cdots ,N.$ | (19) |
| $\begin{array}{*{20}{l}}{{\mathit{\boldsymbol{p}}_{i{\rm{, best}}}}}\end{array} = \left( {{p_{\begin{array}{*{20}{l}}{i{\rm{1, best}}}\end{array}}},{p_{\begin{array}{*{20}{l}}{i{\rm{2, best}}}\end{array}}}, \cdots ,{p_{\begin{array}{*{20}{l}}{iD{\rm{, best}}}\end{array}}}} \right),\\~~~~~~~~~~i = 1,2, \cdots ,N.$ | (20) |
| ${\mathit{\boldsymbol{g}}_{{\rm{best}}}} = \left( {{g_1},{g_2}, \cdots ,{g_D}} \right).$ | (21) |
| $\begin{array}{*{20}{c}}{{v_{ij}}\left( {t + 1} \right) = w \cdot {v_{ij}}\left( t \right) + {c_1}{r_1}\left( t \right)\left[ {{p_{ij}}\left( t \right) - {x_{ij}}\left( t \right)} \right] + }\\{{c_2}{r_2}\left( t \right)\left[ {{g_j}\left( t \right) - {x_{ij}}\left( t \right)} \right],}\end{array}$ | (22) |
| ${x_{ij}}\left( {t + 1} \right) = {x_{ij}}\left( t \right) + {v_{ij}}\left( {t + 1} \right).$ | (23) |
应用粒子群算法进行优化搜索的基本流程如下：
1) 初始化粒子群，包括群体规模、迭代总次数、每个粒子的位置和速度等；
2) 计算所有粒子的适应度，得到粒子的个体极值与全局极值；
3) 计算动态权重，更新粒子的速度与位置，更新个体最优位置与全局最优位置；
4) 判断是否满足终止条件，满足则结束搜索输出最优解，否则继续迭代优化。
2 模拟实验本文采用双生子模拟实验对所提出的数据同化方法进行验证。每组实验分为两部分，采用不同的初始参数。大气扩散模型均采用Gauss烟团模型。其中一组作为参考实验，利用“真实”的输入参数来模拟实际放射性物质扩散情况。另一组作为模拟实验，利用具有不确定性的输入参数作为模型输入，产生与“真实”情况不符的结果。一共进行40次双生子模拟实验，并进行统计分析。
2.1 模拟实验基本条件假定的实验区域是一个10 km×10 km的平坦地区，释放点位于区域中心处。监测数据来源于布置在计算区域内的虚拟观测点位。假设观测得到的数据就是放射性污染物的活度值。
参考实验的释放高度设为20 m，从泄漏源释放气态137Cs，泄漏速率为1×1010 Bq/s。参考实验的风向为西风，风速2.0 m/s。
2.2 输入参数误差设置参考实验采用“真实”的泄漏速率、释放高度、风向与平均风速来模拟，而模拟实验则在“真实”参数基础上添加随机误差，误差大小在一定范围内均匀分布，如表 2所示。
表 2 模型参数输入误差设置
| 输入参数 | 单位 | 表达式 | 误差范围 |
| 泄漏速率Q | Bq/s | Q=γ·Qtrue | γ∈[0.1, 10] |
| 释放高度H | m | H=Htrue+ΔH | ΔH∈[-30, 30] |
| 风向? | (°) | ?=?true+Δ? | Δ?∈[-30, 30] |
| 平均风速v | m/s | v=λ·vtrue | λ∈[0.5, 2] |
表选项
2.3 观测点位与观测误差双生子实验的虚拟观测点位布置在泄漏源(0, 0)的下风向处(东侧)，探测器距离地面1.5 m高。观测点位分为两组：一组为同化组，用于输入数据同化的观测数据；另一组为验证组，用于对数据同化的结果进行验证。
同化组采用弧线形布置，以(0, 0)为中心在下风向半径100、200、500、1 000、2 000、3 000、5 000 m的半圆弧线上，每隔2°设置一个观测点, 一共637个观测点(只有小部分可以捕获到烟羽)。验证组则在泄漏源东侧5 km×2 km范围内，每隔200 m设置一个观测点，一共286个观测点。观测点位分布如图 1所示。
 |
| 图 1 (网络版彩图)两组观测点位设置 |
| 图选项 |
假设各探测器没有明显差别，可以测量泄漏的放射性物质137Cs的活度浓度。探测器测量数据具有30%的误差，其误差随机均匀分布。图 2所示为包含30%误差的观测数据与真实活度浓度的对比。
 |
| 图 2 含有30%随机误差的观测结果与真实活度浓度的对比 |
| 图选项 |
2.4 粒子群搜索算法粒子群搜索算法初始化群体粒子个数为N=200，粒子维数D=3(即释放高度、平均风速、风向)，最大迭代次数T=50。学习因子c1=c2=1.5。采用动态惯性权重，Wmin=0.4, Wmax=0.8，线性减小。3个维度的速度范围分别为：释放高度Vh_min=-1, Vh_max=1；风向Vd_min=-1, Vd_max=1；风速Vs_min=-0.1, Vs_max=0.1。
2.5 结果与讨论图 3-6展示了40组模拟实验结果，将未同化的模型输入参数与同化后的模型参数进行了对比，包括泄漏速率、释放高度、风向与平均风速。其中：星号代表同化前含有随机误差的输入参数，方块代表同化后寻找到的最优估计参数，直线代表参考实验中设置的“真实参数”。
 |
| 图 3 泄漏速率的数据同化结果 |
| 图选项 |
 |
| 图 4 释放高度的数据同化结果 |
| 图选项 |
 |
| 图 5 风向的数据同化结果 |
| 图选项 |
 |
| 图 6 平均风速的数据同化结果 |
| 图选项 |
由图 3-6可以看出，对于泄漏速率、释放高度与风向参数，同化后的估计结果明显与“真实”值更为接近；而对于平均风速，在少数实验的数据同化中对原始参数的改善效果并不理想(实验8、14、20、21、25、26)，相对应的释放速率的同化结果也较差。
根据同为Gauss条件假设的直线Gauss烟羽模型计算公式，
| $\begin{array}{*{20}{c}}{C\left( {x,y,z} \right) = \frac{Q}{{2\pi {\sigma _y}{\sigma _z}v}}\exp \left( { - \frac{{{y^2}}}{{2\sigma _y^2}}} \right) \cdot }\\{\left\{ {\exp \left[ { - \frac{{{{\left( {z - {H_e}} \right)}^2}}}{{2\sigma _z^2}}} \right] + \exp \left[ { - \frac{{{{\left( {z + {H_{\rm{e}}}} \right)}^2}}}{{2\sigma _z^2}}} \right]} \right\}.}\end{array}$ | (24) |
活度浓度与泄漏速率Q成正比，与平均风速v成反比，但当泄漏速率与平均风速的比值一定时，单个参数的取值并不影响活度浓度的计算结果。因此，在同化过程中，即使与“真实”参数不符的粒子，也可能有较高的适应度。如果以泄漏速率与平均风速的比值作为参数，同化的结果普遍比同化前更为理想，都在可接受的范围内，如图 7所示。
 |
| 图 7 泄漏速率与平均风速的比值 |
| 图选项 |
表 3展示了40组模拟实验利用验证组的观测点得到的同化后活度浓度数据与“真实”活度浓度数据的对比。这里采用3种统计参数：相关系数(Corr)、归一化均方差(NMSE)、偏差分数(FB)。它们的定义如下：
表 3 40组双生子模拟实验的数据同化结果
| 实验序号 | Corr | NMSE | FB |
| 1 | 0.998 | 0.12 | -0.09 |
| 2 | 0.998 | 0.16 | -0.11 |
| 3 | 0.998 | 0.10 | -0.05 |
| 4 | 1.000 | 0.03 | -0.05 |
| 5 | 0.995 | 0.35 | -0.09 |
| 6 | 0.988 | 0.60 | -0.11 |
| 7 | 0.997 | 0.10 | 0.00 |
| 8 | 0.991 | 0.39 | -0.05 |
| 9 | 0.995 | 0.23 | -0.07 |
| 10 | 0.999 | 0.30 | 0.18 |
| 11 | 1.000 | 0.11 | -0.10 |
| 12 | 1.000 | 0.03 | -0.06 |
| 13 | 0.998 | 0.11 | 0.03 |
| 14 | 0.996 | 0.39 | -0.09 |
| 15 | 0.999 | 0.05 | -0.05 |
| 16 | 1.000 | 0.05 | 0.06 |
| 17 | 0.995 | 0.24 | -0.07 |
| 18 | 0.999 | 0.11 | 0.02 |
| 19 | 0.997 | 0.15 | 0.04 |
| 20 | 0.999 | 0.06 | 0.06 |
| 21 | 0.999 | 0.06 | 0.09 |
| 22 | 0.996 | 0.24 | 0.19 |
| 23 | 0.998 | 0.11 | 0.04 |
| 24 | 0.997 | 0.19 | 0.05 |
| 25 | 0.997 | 0.17 | 0.08 |
| 26 | 0.996 | 0.15 | 0.00 |
| 27 | 0.998 | 0.25 | -0.14 |
| 28 | 1.000 | 0.02 | 0.02 |
| 29 | 1.000 | 0.10 | -0.09 |
| 30 | 0.999 | 0.06 | 0.07 |
| 31 | 1.000 | 0.09 | 0.08 |
| 32 | 1.000 | 0.02 | 0.05 |
| 33 | 1.000 | 0.06 | -0.08 |
| 34 | 0.995 | 0.23 | -0.09 |
| 35 | 0.998 | 0.10 | 0.04 |
| 36 | 0.995 | 0.22 | -0.04 |
| 37 | 0.997 | 0.12 | 0.02 |
| 38 | 0.998 | 0.17 | -0.03 |
| 39 | 0.998 | 0.15 | -0.10 |
| 40 | 1.000 | 0.21 | -0.11 |
表选项
| ${\rm{Corr = }}\frac{{{\rm{Cov}}\left( {{O_{{\rm{ass}}}},{O_{{\mathop{\rm ref}\nolimits} }}} \right)}}{{\sigma \left( {{O_{{\rm{ass}}}}} \right) \cdot \sigma \left( {{Q_{{\mathop{\rm ref}\nolimits} }}} \right)}},$ | (25) |
| ${\rm{NMSE = }}\frac{{\frac{1}{N}\sum\limits_i {{{\left( {{O_{{\rm{ass,}}i}} - {O_{{\mathop{\rm ref}\nolimits} ,i}}} \right)}^2}} }}{{{{\bar O}_{{\rm{ass}}}}{{\bar O}_{{\mathop{\rm ref}\nolimits} }}}},$ | (26) |
| ${\rm{FB}} = \left( {{{\bar O}_{{\rm{ass}}}} - {{\bar O}_{{\mathop{\rm ref}\nolimits} }}} \right)/\left[ {0.5\left( {{{\bar O}_{{\rm{ass}}}} + {{\bar O}_{{\mathop{\rm ref}\nolimits} }}} \right)} \right].$ | (27) |
由表 3结果可以看出，同化后的结果较为理想，在观测数据含有30%相对误差的情况下，也与参考实验符合得很好：Corr~0.99，NMSE~10-1，大部分FB~10-2。本文提出的同化方法有效地改善了Gauss烟团模型对污染物活度浓度分布的预测结果，同时也能够对泄漏速率、释放高度、风向与平均风速作出更准确的估计。
3 与EnKF算法的对比为了进一步验证本文同化方法的性能，将它与应用广泛的EnKF算法作了简单的比较。
EnKF算法的计算假设条件与结果基本采用文[10]提供的数据，包含同化后的源项估计结果。稍有区别的是文[10]中的实验为600 min的非稳态实验，而本文提出的同化算法以稳态为前提，这一差别会对同化方法的性能有所影响。
源项与气象数据随时间的变化按文[10]中的条件给定(因为包含随机参数，并作了离散处理，故与文[10]不完全相同)。137Cs的释放持续约500 min，在100 min左右有一个释放速率上的突变，并在250 min以后逐渐衰减。风向逐渐从西风转向北风，风速从开始的3 m/s逐渐提升，在375 min达到5.5 m/s后，逐渐减弱到4 m/s。
源项释放位置与监测点位布置如图 8所示，释放点在(0，0)，监测点位布置在距地表 2 m高处，具有30%的观测误差。
 |
| 图 8 释放位置与监测点位布置 |
| 图选项 |
最终的源项估计结果如图 9所示。可以看出，本文同化方法基本能够得到源项较为准确的估计，但是相对真实的变化情况具有一定的滞后性。
 |
| 图 9 本文提出的同化算法对源项的估计结果 |
| 图选项 |
文[10]中采用4种统计参数对结果进行评估，本文方法计算结果与文[10]中的EnKF算法源项估计结果的对比如表 4所示。
表 4 与EnKF算法的源项估计结果对比
| 统计参数 | EnKF算法 | 本文的同化算法 |
| Corr | 0.54 | 0.68 |
| FAC2 | 0.55 | 0.70 |
| NMSE | 0.34 | 0.57 |
| FB | 0.66 | 0.11 |
表选项
Corr、NMSE、FB的定义见式(25)-(27)，FAC2为满足条件0.5≤Oass/Oref≤2.0的数据比例。
由表 4的结果可以看出，对于文[10]中的算例，本文提出的数据同化算法对源项的估计略好于EnKF算法，相关系数达到0.68。但是，笔者改变初始条件，也会出现效果变差的情况。另外，由于代价函数与观测点位模拟值与实测值的相关系数有关，该参数对释放高度并不敏感，导致本文方法对释放高度的估计效果较差，这一问题可以通过改变代价函数的方式进行改进。
4 总结本文提出一种基于Gauss烟团模型的大气扩散数据同化方法，它可以在中尺度平坦地形与均匀稳定的气象条件下，结合观测数据，有效地改善模型的预测结果。对模型的4个输入参数(释放速率、释放高度、风向与平均风速)进行了校正。采用双生子模拟实验验证了本文方法的有效性，并与EnKF算法作了简单比较。本文方法相比于常用的数据同化算法，主要具有计算量小、计算参数需求简单的特点。
本文提出的同化方法目前还存在一定的局限性：1)适用于平坦地形条件与均匀稳定的气象条件，当地形或气象条件复杂时，同化的效果会受到影响；2)本文方法是在Gauss烟团模型的基础上提出的，暂无法应用到其他种类的大气扩散模型中；3)目前方法考虑存在不确定性的模型参数较为简单，未考虑诸如大气边界层高度、摩擦速度、扩散参数等。后续的研究会围绕上述方面对算法进行改进, 并考虑用真实的大气示踪实验结果进行验证。
参考文献
| [1] | NASSTROM J S, SUGIYAMA G, BASKETT R L, et al. The national atmospheric release advisory center modelling and decision-support system for radiological and nuclear emergency preparedness and response[J]. International Journal of Emergency Management, 2007, 4(3): 524-550. DOI:10.1504/IJEM.2007.014301 |
| [2] | THYKIER-NIELSEN S, DEME S, MIKKELSEN T. Description of the atmospheric dispersion module RIMPUFF[R]. Roskilde, Denmark: Ris? National Laboratory, 1999. |
| [3] | SCIRE J S, STRIMAITIS D G, YAMARTINO R J. A user's guide for the CALPUFF dispersion model[M]. Concord, USA: Earth Tech, 2000. |
| [4] | 王醒宇, 康凌. 核事故后果评价方法及其新发展[M]. 北京: 原子能出版社, 2003. WANG X Y, KANG L. The method and new development of nuclear accident consequences assessment[M]. Beijing: Atomic Energy Press, 2003. (in Chinese) |
| [5] | 王跃山. 数据同化:它的缘起、含义和主要方法[J]. 海洋预报, 1999, 16(1): 11-20. WANG Y S. Data assimilation:Its cause, its meaning and main procedures[J]. Marine Forecasts, 1999, 16(1): 11-20. (in Chinese) |
| [6] | EHRHARDT J. The RODOS system:Decision support for off-site emergency management in Europe[J]. Radiation Protection Dosimetry, 1997, 73(1-4): 35-40. |
| [7] | ROJAS-PALMA C, MADSEN H, GERING F, et al. Data assimilation in the decision support system RODOS[J]. Radiation Protection Dosimetry, 2003, 104(1): 31-40. DOI:10.1093/oxfordjournals.rpd.a006160 |
| [8] | ZHENG D Q, LEUNG J K C, LEE B Y, et al. Data assimilation in the atmospheric dispersion model for nuclear accident assessments[J]. Atmospheric Environment, 2007, 41(11): 2438-2446. DOI:10.1016/j.atmosenv.2006.05.076 |
| [9] | ZHENG D Q, LEUNG J K C, LEE B Y. Online update of model state and parameters of a Monte Carlo atmospheric dispersion model by using ensemble Kalman filter[J]. Atmospheric Environment, 2009, 43(12): 2005-2011. DOI:10.1016/j.atmosenv.2009.01.014 |
| [10] | ZHANG X L, SU G F, YUAN H Y, et al. Modified ensemble Kalman filter for nuclear accident atmospheric dispersion:Prediction improved and source estimated[J]. Journal of Hazardous Materials, 2014, 280: 143-155. DOI:10.1016/j.jhazmat.2014.07.064 |
| [11] | ZHANG X L, SU G F, CHEN J G, et al. Iterative ensemble Kalman filter for atmospheric dispersion in nuclear accidents:An application to Kincaid tracer experiment[J]. Journal of Hazardous Materials, 2015, 297: 329-339. DOI:10.1016/j.jhazmat.2015.05.035 |
| [12] | QUéLO D, SPORTISSE B, ISNARD O. Data assimilation for short range atmospheric dispersion of radionuclides:A case study of second-order sensitivity[J]. Journal of Environmental Radioactivity, 2005, 84(3): 393-408. DOI:10.1016/j.jenvrad.2005.04.011 |
| [13] | KRYSTA M, BOCQUET M, SPORTISSE B, et al. Data assimilation for short-range dispersion of radionuclides:An application to wind tunnel data[J]. Atmospheric Environment, 2006, 40(38): 7267-7279. DOI:10.1016/j.atmosenv.2006.06.043 |
| [14] | 刘蕴, 方晟, 李红, 等. 基于四维变分资料同化的核事故源项反演[J]. 清华大学学报(自然科学版), 2015, 55(1): 98-104. LIU Y, FANG S, LI H, et al. Source inversion in nuclear accidents based on 4D variational data assimilation[J]. Journal of Tsinghua University (Science and Technology), 2015, 55(1): 98-104. (in Chinese) |
| [15] | HIEMSTRA P H, KARSSENBERG D, VAN DIJK A. Assimilation of observations of radiation level into an atmospheric transport model:A case study with the particle filter and the ETEX tracer dataset[J]. Atmospheric Environment, 2011, 45(34): 6149-6157. DOI:10.1016/j.atmosenv.2011.08.024 |
| [16] | HIEMSTRA P H, KARSSENBERG D, VAN DIJK A, et al. Using the particle filter for nuclear decision support[J]. Environmental Modelling & Software, 2012, 37: 78-89. |
| [17] | BRIGGS G A. Diffusion estimation for small emissions. Preliminary report[R]. Washington, DC, USA: Department of Energy, USA, 1973. |
| [18] | 包子阳, 余继周. 智能优化算法及其MATLAB实例[M]. 北京: 电子工业出版社, 2016. BAO Z Y, YU J Z. Intelligent optimization algorithms and the MATLAB examples[M]. Beijing: Electronic Industry Press, 2016. (in Chinese) |
