拉格朗日烟团模型的大气扩散系数自适应修正

谈文姬，王德忠，马元巍，吉志龙

（上海交通大学 机械与动力工程学院，上海 200240）

自1986年发生切尔诺贝利事故以来，国际上对核事故后果评价和决策支持系统的研究有了很大进展，特别是在2011 年福岛核事故之后，放射性核素释放到环境中的影响引起了国际社会的广泛关注。核事故后果决策支持系统能对事故发生后放射性核素泄漏到空气中的浓度分布进行实时和预测性计算，为保护周围公众的安全提供数据依据。

目前，用于核事故应急决策支持系统中的大气扩散模型主要有3类：高斯烟羽模型、拉格朗日烟团模型和蒙特卡罗粒子模型［1］。高斯烟羽模型适用于稳定均匀流动及平坦地形的条件，虽模型简单、运算快捷，但难以适用于复杂流场下的实际核事故。蒙特卡罗粒子模型能较为真实地模拟粒子在大气中的输运过程，但计算时间较长，并存在边界稀疏的问题，往往不能满足实时核事故应急响应［2］。烟团模型能适用于非稳定流动及复杂地形的条件，在局地尺度上有较好的精度，应用较为广泛，如欧洲RODOS等核事故应急决策支持系统均采用该模型［3］。

描述烟团扩散能力的大气扩散系数，是影响拉格朗日烟团模型的重要参数。传统的拉格朗日烟团模型采用的扩散参数是由Pasquill-Gifford（P-G）曲线确定的。P-G 曲线是根据理想条件下的实验结果拟合得到的，与实际情况的不符将给烟团模型带来偏差。本文针对拉格朗日烟团模型不能实时反映局地扩散能力的缺点，提出大气扩散系数自适应的改进方法，并使用实际大气扩散示踪实验的数据，对大气扩散系数自适应修正方法进行检验。

1 烟团模型的大气扩散系数自适应修正

1.1 烟团模型简介

烟团模型是将释放的污染物假想成离散的烟团，对每个烟团中心的输运过程进行模拟。每个烟团在某一特定时间间隔内固定不动，根据此刻固定住的烟团计算浓度。然后烟团在下一个时间步长内继续移动，大小和强度继续变化，直到下次采样时间再次固定。在基本时间步长内，接受点的浓度为周围所有烟团采样时间内的平均浓度总和。有风的条件下，第i个烟团在某一点的浓度［4］为：

式中：Ci（x，y，z）为第i个烟团（x，y，z）点的浓度，g／m3；Q（i）为第i个烟团的强度，g／s；σx、σy、σz分别为沿风向、水平方向、垂直方向的扩散参数，m；（xc（i），yc（i），zc（i））为第i个烟团的中心位置，m。

计算污染物浓度及分布时，必须知道源强、风场、有效源高和大气扩散参数。其中，源强往往是通过测量获得或由设计基准事故给出，气象信息由当地气象测量站获取。大气扩散参数对扩散后果有较大的影响，但目前仍为理想条件下的示踪实验数据计算得到，这些数据在核事故后果评价中与实际情况有较大偏差，应提出新的方法对其进行修正。

1.2 大气扩散系数自适应修正方法

大气扩散参数（σx、σy、σz）表示污染物在随机湍流场中的扩散能力和散布范围，它与下风向距离x、大气稳定度和下垫面粗糙度等相关。大气扩散参数随x 及下垫面粗糙度的增大而增大，随大气稳定度的增大而减小。

大多数模型均通过Pasquill气象稳定度与P-G 曲线来确定大气扩散参数，该曲线是由大量扩散实验和理论分析得到的，具体确定扩散参数的方法为：首先根据云况和日射及地面风速选择大气稳定度，然后根据扩散曲线读出不同下风向距离处的扩散参数。RIMPUFF模型中使用的Karlsruhe-Jülich（K-J）大气扩散系数即采用此方法，K-J大气扩散系数是根据大气稳定度及源项的释放高度将大气扩散系数进行分类，共分为18套大气扩散系数［5］。但由于P-G曲线是依据平坦理想条件下的实验得到的，对于不同源高、下垫面等条件下的扩散曲线还有待进一步修正。

本文提出一种根据当地的实测数据（源项释放数据、观测点位置和测得的浓度数据、气象数据）及烟团模型的浓度计算公式，对大气扩散系数进行实时估计。这种将上一步大气扩散系数的估计结果作为修正的大气扩散系数进行下一个时间步长运算的方法，称为大气扩散系数的自适应修正。将得到的适用于当地的实时大气扩散系数作为修正后的系数继续进行扩散计算，较使用理想条件下的实验得到的扩散参数更符合实际情况，能提高拉格朗日烟团模型的准确性。

假设扩散参数为下风向距离x 的幂函数，即：

式中：py、qy、pz、qz为大气扩散系数。

若确定py、qy、pz、qz，就 可 得 到 适 合 当 时当地地形气象条件的σy和σz。这4个大气扩散系数可利用最小二乘法确定。为了权衡观测数据的精度，引入一标志测量精度的权重g，使地面浓度的计算值Cmi（Cmi=C（xi，yi，0））与实际观测值COi之差的平方和S 最小，则S 可表示为：

式中：n为实际观测的数据个数；gi为采样点的权重，gi=COi／CO，max，CO，max为观测点采集数 据中的最大浓度。

采用梯度迭代法计算求解使S 为最小的非线性方程组，可得到实时的大气扩散系数。本文根据上述方法，针对拉格朗日烟团模型进行修正，开发了相应的大气扩散系数自适应模型。

2 模型验证

2.1 Kincaid扩散实验

MVK 是荷兰国家公共健康与环境研究院开发的大气扩散模型验证对比工具［6］，包括Kincaid、Indianapolis、Copenhagen、Lillestrom 4个实验数据集，定量的统计分析工具（BOOT程序）和定性的图形分析工具（SIGPLOT 程序），以及可用于评价模型的准确性、适用性和揭示模型的错误和数据的坏点。目前，MVK 已广泛用于验证和比较自主开发的大气扩散模型［7］。

Kincaid实验是EPRI烟羽模型验证和改进项目的一部分。EPRI空气质量数据中心提供了实验时的气象数据，包括气压、温度、湿度、风速和风向。Kincaid电站坐落在美国伊利诺斯州（北纬39.59°，西经89.49°），平均海拔约为180m，其周围是平坦的农田及湖泊。图1为以Kincaid电站为中心的80km×80km 范围的等高地形图。图1中为UTM 坐标，x正方向为东，y正方向为北，网格划分为100m×100m。

图1 以Kincaid电站为中心的80km×80km 范围的等高地形图Fig.1 80km×80km wide topographic contour centered at Kincaid Power Plant

实验时，六氟化硫将从电站中的一座高187m、直径9 m 的烟囱顶部释放。与释放点距离为1、2、3、5、7和10km 的圆弧上设有观测点，采样并记录数据。实验数据集包含约171h的数据。

2.2 数据验证方法

BOOT程序是MVK中的统计工具，能将模型预测值与实验观测值进行比较，且可计算出各种统计参数。除了模型预测与实验观测的平均值MEAN、标准差SIGMA 外，其他数据还包括平均偏差BIAS、标准化均方差NMSE、相关系数CORR、比例偏差FB、FA2等，这些统计参数能为模型评价提供依据［8］，其计算公式为：

式中：CM为模型预测浓度，μg／m3；CO为观测浓度，μg／m3；σM为模型预测浓度的标准差，μg／m3；σO为观测浓度的标准差，μg／m3。FA2为0.5≤CM／CO≤2的几率，即模型计算值误差在50%以内的几率。

3 计算结果与讨论

3.1 大气扩散系数最优拟合方法

拟合大气扩散系数所需数据包括：观测站数据（观测站相对于源的位置坐标和测得的浓度数据）、气象数据（当时的风速、风向，由气象站提供）、源项释放率（实验设计提供）。1980年4月20日14：00的测量数据列于表1。其中源项释放处的高度为源高处，源高处风向以正北方向为0°。

表1 Kincaid实验1980年4月20日测量数据Table 1 Measurement data of 1980-04-20in Kincaid experiment

大气扩散系数初值的设置会对计算得到的大气扩散系数带来较大影响。不同的初值求解得到的大气扩散系数不同，需对解得的大气扩散系数进行筛选。本文提出一种大气扩散系数最优拟合方法，根据表1中的数据，当大气扩散系数py、qy、pz、qz初值为（1，1，1，1）时，求解得到大气扩散系数A=（0.739，0.881，0.838，0.959）；当初值为（0.8，0.8，0.8，0.8）时，解得大气扩散系数B=（3.84，0.681，0.553，1.007）。将A 和B代入拉格朗日烟团模型时，两者计算得到的观测点坐标处的浓度相同，无法立即判断这两组大气扩散系数的适用性。在这种情况下，将A、B 两组大气扩散系数作为初值代入下一个时间步长的迭代计算。下一步计算结果表明，使用A 计算的浓度较使用B 计算的浓度更接近观测值，则选择A 代入下一步的计算。重复该过程即可得到大气扩散系数的最优解。

3.2 示踪实验验证结果

用MVK 中的Kincaid数据集对大气扩散系数自适应修正方法进行检验，并对大气扩散系数自适应修正方法的烟团模型与K-J大气扩散系数烟团模型的计算结果进行比较，以判断该方法的可行性。根据Kincaid实验记录，从K-J大气扩散系数中选择与源项释放高度为180m 和D 类大气稳定度相对应的一组大气扩散系数（0.208，0.903，0.307，0.734）［9］。

对实验值（OBS）、K-J大气扩散系数烟团模型的计算值和大气扩散系数自适应修正的烟团模型（APM）的计算值进行对比的结果列于表2～6。

表2 1980年4月20日3h18个有效观测站数据检验结果Table 2 Result of comparison of 18valid experiment data in 3hon 1980-04-20

表3 1980年4月25日6h45个有效观测站数据检验结果Table 3 Result of comparison of 45valid experiment data in 6hon 1980-04-25

表4 1980年5月1日6h20个有效观测站数据检验结果Table 4 Result of comparison of 20valid experiment data in 6hon 1980-05-01

表5 1980年5月4日6h19个有效观测站数据检验结果Table 5 Result of comparison of 19valid experiment data in 6hon 1980-05-04

由表2、3可看出，与大气扩散系数烟团模型相比，大气扩散系数自适应修正方法计算得到的平均值、标准差更接近于观测值，比例偏差FB更接近于0、FA2更接近于1。采取大气扩散系数自适应方法的烟团模型的各项统计参数均优于KJ大气扩散系数烟团模型的。就模型计算值与观测值误差在50%以内的几率而言，准确性提高约20%。表2、3、6的数据修正效果优于表4、5的，这是由于风向变化较大对修正效果产生了负影响。1980年4月20日、4月25日及5月7日3天的风向变化较小（5°～40°），而5月1日、4日两天的风向变化较大，有时达180°。

表6 1980年5月7日9h36个有效观测站数据检验结果Table 6 Result of comparison of 36valid experiment data in 9hon 1980-05-07

图2、3分别为大气扩散系数自适应修正方法模型及K-J大气扩散系数模型的预测值与观测值的分位数图（Q-Q 图）。

图2 APM 模型的Q-Q 图Fig.2 Q-Q plot of APM model

图3 K-J模型的Q-Q 图Fig.3 Q-Q plot of K-J model

从图2和3可看出，两个模型的预测值较实验值均偏大，这种保守的估计，在核事故浓度预测上是可被接受的。APM 模型点的个数多于K-J模型的，这是由于APM 模型预测的观测点处浓度的非零值较多。修正模型的预测值的Q-Q 图更接近y=x 直线，表明APM 模型的预测值较K-J模型的更符合观测值。

3.3 讨论

由于传统烟团模型中的扩散参数是根据大气稳定度及释放高度的条件选取的，可供选取的值是以理想条件下的实验结果为依据，因此，扩散参数在复杂流动和地形条件下，适用性受到了很大局限。为改善扩散参数的适用性，利用当时当地的气象数据、源项数据及观测数据，以最小二乘法求解大气扩散系数，将符合实际情况的大气扩散系数代入模型中，以达到系数自适应修正的效果。通过计算，将传统扩散参数的烟团模型与系数自适应修正的烟团模型的计算结果进行比较，后者与观测数据一致性的统计结果得到明显的提高，证明大气扩散系数自适应修正方法是有效的。

在拟合求解大气扩散系数时，一个时间步长中风速风向的影响会累积到下一个计算时间步长内。缩短计算时间步长可减少风速风向变化带来的影响，更准确地求解大气扩散系数。但在计算过程中要求有与计算时间步长相匹配的气象数据。

由计算过程可发现，大气扩散系数自适应修正方法需要有较为准确的观测值和释放数据。由于对观测值的数量和精度均有一定的要求，因此，在厂区内如何合理有效地布置监测站是需要考虑的问题。事故发生时，破口出现的位置及泄漏的数据应尽可能地详细，因此安全壳内各部件的监测系统需灵敏、准确，且需发展源项反演技术，以便获得准确的释放数据。

4 结论

本文针对拉格朗日烟团模型提出了一种大气扩散系数自适应的方法，建立了动态的大气扩散系数自适应修正的拉格朗日烟团模型。通过Kincaid实验数据验证表明，大气扩散系数自适应修正方法的烟团模型优于传统P-G 曲线扩散参数的烟团模型。基于新方法建立的动态大气扩散系数自适应修正烟团模型可提高大气扩散模型计算结果的准确性。

［1］ 胡二邦，姚仁太，宣仪仁，等.核事故后果评价系统的进展与比较［J］.辐射防护通讯，2003，23（1）：6-12.HU Erbang，YAO Rentai，XUAN Yiren，et al.Comparison and advance of consequences assessment system for nuclear accident［J］.Radiation Protection Bulletin，2003，23（1）：6-12（in Chinese）.

［2］ 曲静原，曹建主，刘磊，等.我国核应急决策支持系统研究开发的现状与展望［J］.原子能科学技术，2001，35（3）：283-288.QU Jingyuan，CAO Jianzhu，LIU Lei，et al.Status and perspective on the research and development of the Chinese decision support system for nuclear emergency management［J］.Atomic Energy Science and Technology，2001，35（3）：283-288（in Chinese）.

［3］ EHRHARDT J，WEIS A.Development of a comprehensive decision support system for nuclear emergency in europe following an accidental release to atmosphere，FZKA-5772［R］.Karlsruhe：Forshungszentrum Karlsruhe，1996.

［4］ 蒋维楣.空气污染气象学［M］.南京：南京大学出版社，2003：102-111.

［5］ 曲静原，曹建主.RODOS4.0与RODOS的未来发展［J］.辐射防护，2002，22（4）：198-206.QU Jingyuan，CAO Jianzhu.The future development of RODOS4.0and RODOS［J］.Radiation Protection，2002，22（4）：198-206（in Chinese）.

［6］ OLESEN H R.Model validation kit-status and outlook［J］.International Journal of Environment and Pollution，2000，14（12）：65-76.

［7］ OLESEN H R.User’s guide to the model validation kit［R］.Denmark：National Environment Research Institute，2005.

［8］ ELEVELD H.Application of a methodology to validate atmospheric dispersion models［C］∥Proceedings of the 7th International Conference on Harmonization Within Atmospheric Dispersion Modelling for Regulatory Purposes.Italy：［s.n.］，2001.

［9］ THYKIER-NIELSEN S，DEME S，MIKKELSEN T.Description of the atmospheric dispersion module RIMPUFF［R］.Denmark：Risoe National Laboratory，1999.