基于变分同化的化学危害扩散预测研究∗

蒋金利 诸雪征 顾 进 韩朝帅， 李 廷

（1.陆军防化学院 北京 102205）（2.军事科学院防化研究院 北京 102205）

1 引言

当战场发生化学袭击或者化工厂发生泄漏，爆炸事故后，毒剂被分散成为蒸气、气溶胶进入大气中，形成毒剂云团。而毒剂云团将随气流一起运动，同时向四周扩散，形成化学危害区域［1］。指导部队进行有效化学防护的前提是对战场化学危害区域进行准确的界定、检测及预测，只有掌握毒剂浓度的时空分布规律，准确描绘毒剂云团的危害地域，才能预测未来一段时间窗口内的毒剂浓度分布，进而分析出该段时间内的毒害作用效应。当前对于化学危害的预测主要通过大气扩散理论和数学扩散模型［2］来推导未来一段时间的毒剂浓度分布，从而对化学危害的发展态势做出评估，但是基于大气扩散模式的化学云团浓度预测具有一定的偏差，且战场气象、地形条件复杂多变，计算结果往往会有较大的误差。

数据同化是一种最初来源于数值天气预报的数据处理技术，现在已成功用于对大气污染物浓度的预测中来。在化学危害预测领域引入数据同化算法，以更快的速度和更高的精准度得出危害场未来一段时间的浓度走势及分布图，在保证速度的同时尽最大可能地为防护和救援提供准确的决策依据。当前关于污染物扩散数据同化主要有逐步迭代法［3］，最优插值法［4］、卡尔曼滤波法［5］、集合卡尔曼滤波法［6］、遗传算法［7］、神经网络算法［8］和机器学习［9］等，但逐步迭代法和最优插值法精度偏低，属于早期的同化算法；卡尔曼滤波系列算法对集合数目的要求偏高，难以快速应用于危害环境；遗传算法，神经网络和机器学习是新兴智能优化算法，用于数据同化领域还处于初期探索阶段。所以，本文选用变分同化［10～11］方法，既避免了精度不够的问题，也没有集合卡尔曼滤波对数据样本这么高的要求。

因此，本文结合变分同化理论，结合经典的高斯扩散模式，对化学危害扩散中数据同化算法进行构建和案例分析。

2 化学危害扩散模型

高斯扩散模式［12］是经典的污染物扩散模型。在大气中，几乎处处时时都存在着各种不同尺度的湍流运动，常见的湍流理论主要有梯度输送理论和统计理论两种。但是无论从哪一种理论出发，高斯扩散模型是可推导的。由湍流扩散梯度输送理论建立起来的扩散方程，当扩散系数为常数时，可以得出正态分布形式的解；当气流平稳均匀时，扩散物质的粒子分布也遵循正态分布。尽管实际大气中，并不能严格满足相关条件，但是在不知道实际物质浓度分布的情况下，高斯分布起码是一种较好的近似。在实际操作中，高斯扩散模式所需要的气象条件少，计算速度快，也可以较好地反映浓度的时空分布。

2.1 瞬时扩散模型

瞬时释放又称为高斯烟团模式［13～14］，可以看成是将高斯烟羽模式分成数个烟团。瞬时点源排放的迁移模式的基本方程:

同时，假定在原点（0，0，0），即x=y=z=0处有一个瞬时点源排放；当时间趋近零时，浓度也趋近零，当时间趋向无穷大时，浓度也趋近于零；整个扩散过程，满足质量守恒定律，且假定时间为零时，排放物的总质量为M。

在满足上述条件之后，微分方程式（1）的解为

式中:C（x，y，z，t）为 t时刻，空间中任一点（x，y，z）的浓度，单位为（g/m3）；Qi，t为单个烟团的浓度；σx，i，t，σy，i，t，σz，i，t为t时刻各方向上的扩散参数；xi，t，yi，t，zi，t为某个毒剂烟团的位置参数。

2.2 连续扩散模型

点源扩散是高斯扩散模式的主要研究对象之一，它是指在扩散进行了一段持续时间且气象状况平稳（风向恒定）的条件下所进行的扩散。在这个过程中，我们把释放源的位置作为坐标原点，且源强连续、均衡地释放毒剂气体或气溶胶，释放出来的气体延风向不断进行扩散运动。如果以风向方向为x轴，垂直方向为y轴建立坐标系的话，不难预见，毒剂气体将在下风方向形成圆锥形的云团，并且在x轴的任一截面处都符合二维正态分布。连续扩散下，高斯扩散模式又叫高斯烟羽模型［14］。

地面点源连续排放的粒子迁移工程为

上式中各参数含义与瞬时扩散时相同。在满足质量守恒定律的前提下，该微分方程的通解为

按照所排放毒剂的性质的不同，标准解一般可分为两种情况，即地面反射与地面不反射。

在扩散过程中，当毒剂沉降到地面时，存在反射回大气层或者不沉降到地面的情况，数学上的表达为满足积分:

将式（5）代入式（6），可得:

其中，h为源项与地面的高度，如果是地面源，那么其值为零。其他项的含义和瞬时扩散方程相同。

2.3 扩散参数

在上述扩散过程中，确定扩散参数σx，σy，σz其实是很困难的，一般来说，需要特定条件下的气象观测，并进行大量的计算工作。但是在实际操作中，往往只有常规的观测资料。为了能方便高效地根据已有资料推算出扩散参数，研究人员进行了深入的研究，得出了一些行之有效的方法，包括帕斯奎尔扩散曲线法、GB/T 公式以及 Briggs［15］公式等方法。这些方法将大气划分A-F为六个稳定度级别，稳定度划分如表1所示。

表1 大气稳定度分级

Pasquill和Gifford也给出的不同稳定度下的σy与σz随下风距离变化的经验曲线。在P-G曲线的基础上，我国根据实际情况制订了国家标准《制定地方大气污染物排放标准的技术方法》［16］，确定大气稳定度后，通过下面两个公式来计算扩散参数:

式中，x为下风方向距离；a1，a2及γ1，γ2为扩散参数，都是无量纲常数。

3 化学危害扩散数据同化模型

根据变分同化理论，结合高斯扩散模式，对战场化学危害问题进行数学建模。变分同化方案的基本功能:在化学危害事件发生之后，结合当时的气象条件，观测数据等，输入变分同化的工作流程。当误差统计减小到可接受的范围后，输出危害场的浓度分析值，反映当时状态下的危害程度，为后续的救援行动以及事后评估提供决策依据。下面依据变分理论，进行方案设计，构建基于高斯扩散模型的战场化学危害同化模型。

3.1 变分算法原理

四维变分方法（Four—Dimensional Variational Algorithm:4D-Var）最早由我国学者顾震潮［17］提出，是以三维变分法方法（3D-Var）为基础发展起来的，是3D-Var在时间维上的扩展，实现了整个同化时间窗口观测数据的同化，预测精度更高，适用范围更广。

四维变分同化的公式:

式中Hi表示第i时刻的观测算子，表示第i时刻的观测值，Mt0，ti表示将预报模式从时刻积分到ti时刻的模式预报算子，故:

指由模式算子积分得到的第i时刻的模式变量。同理可得:

指将第i时刻的模式变量经第i时刻的观测算子Hi映射到观测场。

1994年，Courtier等为了避免大量的计算，降低运算成本，提出增量方法［18］，它不直接求解分析变量，而是求解分析增量并修正预测场。

定义分析增量:

同时，对预报模式算子进行泰勒展开:

代表观测算子的切线模式或其他线性近似模式。舍去式中的高阶无穷小项，并令向量:

从而得到增量变分公式:

3.2 化学危害扩散变分同化方案设计和模型建立

3.2.1 化学危害变分同化方案设计

依据变分同化原理，结合高斯扩散模式以及化学危害初生云形式，设计了变分数据同化方案，对工作流程进行了梳理。

该方案主要包括六个部分:外部条件获取模块，源项信息估算模块，观测数据处理模块，高斯扩散预测模块，变分同化模块以及结果分析模块。

外部条件获取模块:获取相关的风场信息，包括风向、风速、温度等；同时确定危害地域的地形地貌信息，有无大型山体、农田等。

源项信息估算模块:主要是为高斯扩散模型提供源项信息，生成同化背景场。

观测数据处理模块:处理各类装备收集到的各种观测资料，主要对浓度以及风场的数据进行处理，在同化时间窗口内生成观测向量以及数据文件。

高斯扩散预测模块:将所获得的信息条件输入，进行预测，生成背景场信息，传输给同化模块。并不断更新时间域，不断向前预测。

变分同化模块:综合处理背景资料，观测资料，通过所建立的变分同化模型得出分析场。

结果分析模块:输出分析结果，并分析过程误差。

3.2.2 化学危害变分数据同化模型构建

1）变分同化梯度函数

变分同化算法的方案设计部分，是在空间场内求解最优分析解的数据同化算法。下式是上一章节提到的变分法目标函数的基础公式:

对标量函数J关于自变量x求导，可以得到其梯度:

求解梯度公式，可得分析值:

式中，x为需要求解的状态参量，xa称为分析值（analysis value）；xb是模型的背景场，在本文中，它是由高斯扩散模型输出的参量；B是背景误差协方差矩阵；y为观测场数据，H是观测算子，它的作用是将观测场数据的类型映射到与状态量相一致；R为观测误差协方差矩阵。

误差是整个同化过程中不可避免的话题，从数学的角度来说，可以把同化过程理解为最小化误差的过程。误差在同化开始前主要体现在背景误差协方差矩阵B上，它是整个同化过程的重要组成部分，和观测误差协方差共同决定了分析场的误差。因此，在同化过程开始之前，必须对误差进行修正，尽可能使分析值误差降低，提高准确度。

一般情况下，确定误差协方差矩阵前要有一个假设，即观测误差与背景误差不相关，即:

式中:εb为背景误差；εo为观测误差，上标T表示矩阵的转置。在实际分析中，这也是完全合理的假设，背景场参量与观测场参量相互独立。

2）背景误差协方差矩阵B

背景误差协方差矩阵表示了预报误差的概率分布函数（Probability Distribution Function，PDF），是一个高斯型的函数分布。它的准确与否直接决定，其同化系统所产生的分析场参量大概只有15%来自观测资料，而剩余的85%均来自背景场［19］。在物理意义上，背景场是同化前的该地域的状态场，它将同化循环中的观测数据不断融合进最新的分析场。所以说，如何对分析场进行准确有效的更新，不仅取决于模型的质量，也取决于背景场质量。其定义如下:

式中xtrue代表的是真值，εb为背景误差，由于很多时刻真值未知，所以，很多时候B矩阵也表示为

式中，上横线代表均值，其余含义与上式相同。由此可以看出，B矩阵是正定对称矩阵，而且特征值均为非负值。

假设背景误差为向量(e1,e2,e3)，那么其背景误差协方差矩阵可以写成:

从上面这种形式下的B矩阵可以看出，其中的每个元素反映的都是两个位置上背景误差的相关性，和两个变量之间的相关性。对角线元素或者对角矩阵代表的是背景误差方差，而非对角元素和矩阵代表的则是协方差。

定义x，x0两点间的背景场相关系数b(x,x0)为两点间距离dist(x-x0)的函数:

3）观测误差协方差矩阵R

R是观测误差协方差矩阵，它包含了关于观测误差的统计信息。对于很多观测系统来说，不同位置上的观测误差在统计上是相互独立的。同样地，在确定观测误差协方差矩阵之前，假定观测误差是无偏的，εo的数学期望为零。

观测误差一般包括观测算子误差、仪器误差和代表性误差。观测算子误差主要是因为对两种不同种类的数据之间联系的理解不足或者理解过于片面导致的，主要出现在非直接观测情况中。代表性误差一般来源于两个方面:一是模型不完美，得不到完全准确的测量；二是来自有限的模型分辨率，分辨率3km的观测系统肯定无法反映10km范围的区域状态。这类误差的减小可以利用高密度观测，提高模式分辨滤等手段。仪器误差是出厂时带有的，一般会带有数据参数，考虑设备老化，使用时间长等实际情况可以对参数进行适当调整。除此之外，就目前的研究水平而言，国内外也鲜有好的解决办法。一般而言，对R的处理相当简化，在观测点之间相互独立的情况下，可以将其简化为对角矩阵处理，方便下一步的计算。

3.2.3 最优化算法

共轭梯度法是介于最速下降法与牛顿法之间的一个方法［20］，它仅需利用一阶导数信息，但克服了最速下降法收敛慢的缺点，又避免了牛顿法需要存储和计算Hessen矩阵并求逆的缺点。本文结合化学危害问题的研究背景以及变分同化的计算要求，选定共轭梯度法为最优算法。

对于一般二次型目标函数来说:

其中G是正定矩阵。由扩展子空间原理可知，若果矩阵G有n个共轭方向，那么从任意的初始点出发，至多做n次精确一维搜索，就可以得到目标函数唯一的极小点。根据这一理论，可以导出对于正定二次矩阵型函数的共轭梯度法。用g表示梯度，d表示下降方向，算法步骤如下:

4 案例分析

为了验证所建模型的同化能力以及各参数对同化质量的影响，本节利用数值模拟实验进行分析。

首先运行危害物质的扩散模型，得到危害区域内各个网格点的浓度，将其作为背景浓度场，并将此时即为记为同化窗口的0时刻；继续向前运行扩散模型，简单起见，将所得到的模型预测值作为“真值”；其次，在危害区域内设定观测点位置，将观测点位置处预测浓度输出并记为观测数据；同时，保持参数不变，运行另一组扩散模型，并给初始预测结果加上10%随机误差，用上述得到观测值同化带有误差的预测值，检验同化效果以及收敛性。

条件设定:假设某平原地区化工厂发生烟囱内单一化学物质泄漏，并且已获取当时的源项信息以及气象信息，厂内为预防事故已经配备了具有输出功能的观测设备，其他条件如表2所示。

表2 实验条件设定

为了使实验尽量简单，设定释放方式为均匀释放，并且理想化外部条件，不考虑降雨、雷电等恶劣天气影响。除同化效果分析之外，本节还将分析收敛性以及过程误差。

首先，给出背景场，真实场以及模型预测场的扩散图。

图3 浓度场比较图

从图中可以分析出来，刚刚开始扩散时所预测出来的浓度图，本试验的背景场，在整个扩散区域内都有着或多或少的浓度分布，但是30min扩散后，除危害物质初生云之外，其他地域的物质浓度似乎为零；同时，加上随机误差后，可以明显的看出预测值与真实值之间有了差异。

取迭代次数为1000次，下图是试验的同化效果图，收敛性以及误差变化。

图4 同化实验结果图

从图中可以看出，通过同化，分析值比观测值跟更接近真值，而背景值在初始状态下与其他三条曲线有着明显偏差，在一定迭代次数后，逐渐与其他几条曲线拟合；由于观测误差协方差矩阵R为对角矩阵，且对角线元素为均方差，为固定值，所以不管迭代次数是多少次，同化时间窗口有多长，观测场方差都不变；同化过程的单次误差虽然各次不同，但是在迭代200次左右开始，便趋于收敛，变的稳定；同样地，同化的背景场方差也是随着迭代次数的增多，开始逐渐收敛，次数也大概在200次左右。

图5是迭代次数20，50，100，200次时的同化图。

图5 不同迭代次数下的同化效果

从图5中对比分析可以看出，在200次迭代之后，同化基本达到效果，各曲线之间拟合程度加大，开始收敛。这也与图5中分析200次迭代后同化过程开始收敛的结果相一致。危害区域的浓度信息经过同化后，危害物质浓度逐步逼近真值，围绕其上下波动。并且随着时间的推进和观测数据的增加，同化结果越来越稳定，逐渐消除了初始时刻所加的10%的随机误差，达到收敛状态。

5 结语

本文主要对四维变分数据同化在化学危害扩散预测中的应用进行了系统研究和验证。首先，分别对高斯扩散模式的瞬时释放模式和连续释放模式进行了研究，提出适用于化学危害预测的扩散参数；其次，构建了适用于化学危害预测的变分同化算法，对同化方案、背景误差协方差矩阵、观测误差协方差矩阵以及最优化算法进行了建模，明确了寻优步骤；最后，通过Matlab仿真，对算法的同化质量和收敛性进行了试验和对比分析。但是，由于为证明参数的影响，示例分析部分应用了理想化数值实验，简化了实验条件，接下来的研究还要更加细化实验条件，做更细致的分析。