二维平流扩散模型下的集合卡尔曼滤波模拟同化研究

摆玉龙，尤元红，邵 宇，徐宝兄

西北师范大学物理与电子工程学院，甘肃 兰州 730070

二维平流扩散模型下的集合卡尔曼滤波模拟同化研究

摆玉龙，尤元红，邵 宇，徐宝兄

西北师范大学物理与电子工程学院，甘肃 兰州 730070

为提高环境数值预报水平，构建了一个针对污染物扩散的模拟数据同化系统。采用集合卡尔曼滤波方法对二维平流扩散模型的状态变量进行了实时校正，实现污染物浓度的实时模拟预报，完成了敏感性实验中集合数目变化、观测方差变化和同化窗口长度变化研究。比较考察观测点位置与污染源距离不同时的预报效果，探讨了优化条件下的同化策略，提出一种根据距离远近动态调节卡尔曼增益权重的方法。在集合数目较小时，可降低计算代价，得到优化的同化效果。

数据同化；平流扩散模型；集合卡尔曼滤波

海洋与大气中的污染物扩散、迁移和转化问题是环境监测中的焦点问题，用科学的方法监测预测污染物浓度，开展实时、全面、准确的环境质量预报工作具有重要意义[1-3]。数据同化技术可用来构建海洋和空气质量预报系统[4-6]，基于变分方法和集合卡尔曼滤波的混合同化方法在空气质量预报领域受到了广泛关注[7-9](如将卡尔曼滤波与T106气象因子相结合对空气质量进行预报，得出卡尔曼滤波比传统的多元回归法更具物理意义、准确率更高[10])。然而，不同的数据同化方法各有优缺点，且计算量大，发展更为有效的数据同化方法，提高海洋、大气环境观测成效，仍是当前面临的严峻问题[11]。作为海洋和大气领域数值预报中的核心方法[12-13]，集合数据同化方法因集合样本数过少和模型误差等问题出现滤波发散,降低了预报精度；另外，当存在多个污染源，点位分布相对稀疏的情况下，从少量的观测中给出区域空间内污染物浓度空间分布较为困难。因此，借助污染物扩散传输模式采用数据同化方法来解决这一问题值得研究与探索。

本文在常用的二维平流扩散模型基础上[14]，利用集合卡尔曼滤波方法建立能够对污染物浓度进行实时预报的模拟同化系统，将5个同种污染源随机置于模拟风速场中来模拟实际污染物在空气中的扩散过程，利用集合卡尔曼滤波对模型的状态变量进行更新和校正，通过分析比较集合数目变化、观测方差变化和同化窗口长度变化对同化效果的影响，比较研究不同观测点位置情况下的预报效果，提出优化条件下的同化策略。

1 同化方案设计

1.1 集合卡尔曼滤波

集合卡尔曼滤波是由Evensen[15]提出的一种基于蒙特卡罗方法的顺序同化算法，其特点是通过统计有限数量的样本来估计模型状态向量的统计特征。与传统的卡尔曼滤波相比较，集合卡尔曼滤波引入了集合思想，扩展其在非线性系统中的应用，解决了传统方法在实际应用中的诸多问题。集合卡尔曼滤波的实现过程包含预报和更新2部分，在预报时，将状态向量的所有样本在平流扩散模型中计算获得状态向量的预报值；更新过程中，利用状态向量协方差矩阵以及观测误差协方差矩阵权衡预报值与观测数据，从而得到状态向量的最优估计。

集合卡尔曼滤波算法的执行过程如下：

1)随机产生状态变量样本

以任意一组随机数作为初始样本，并以这组样本作为数据同化的初始值，同时导入观测数据对样本进行实时更新。

2)模型预报

通过求解样本状态的控制方程得到状态向量在第k个同化步的预报值，假设集合的个数为N，则系统的状态方程为

(1)

3)同化污染物浓度观测资料，更新模式预报污染物浓度场。在每一个同化步通过公式(2)对状态向量进行更新，以获取状态变量的分析值。

(2)

(3)

(4)

系统运行过程中，在有观测的时刻，利用观测信息在观测和模型误差分别加权的基础上对模型状态进行更新，从而获得模型状态的后验优化估计；状态更新后，模型利用新的状态继续向前积分，直到获得新的观测信息。图1概念性的反映了在顺序同化中的滤波过程(图中在i+1时刻，背景场通过观测场修正，得到分析场，向前演进)。

图1 顺序同化滤波的概念示意图

1.2 局地化技术

集合同化过程中通常采用局地化方法解决虚假相关问题[16-17]。针对某个模式点，将观测与该模式点上变量的协方差乘以一个随距离增加而减小的连续函数。该函数是五阶分段有理函数，并且随着距离的增大而单调下降，其表达式如下：

(5)

在计算同化点与观测点距离的基础上，引入一个简单的权重因子(w)。在各点同化的过程中，根据观测距离远近，动态调节卡尔曼增益。

1.3 实验设计

通常，可以通过设计观测系统仿真实验(OSSE)来检验数据同化系统的性能[18]。OSSE实验一般需要包括系统状态空间模型、观测模型、模拟“真实”数据和模拟“观测”数据等。在一般对流扩散理论中,常常认为扩散质的存在不会改变流体质点的流动特性；同时，在整个运动过程中,流体质点之间不发生扩散质的转移，扩散质的扩散完全是由带有扩散质的流体发生掺混的结果[14]。

二维平流扩散模型可用于研究流体扩散、平流过程，已被成功应用于西湖水质总磷浓度预测[19]，计算水体交换能力，分析水体交换的主要因素[20]。本文在计算中忽略污染物在空气中的物理、化学、生物降解,仅考虑污染物受风作用的输移和扩散过程。简化后的平面二维平流扩散方程为

(6)

式中：c为污染物的浓度；u，v为风速沿x，y方向的分量；Dx，Dy为污染物沿x，y方向的扩散系数。污染物在空气中向四周扩散的速率相同，所以可认为污染物沿x，y方向的扩散系数相同，上面的模型可以进一步简化为如下形式：

(7)

式中：σ为污染物的扩散系数。

按照Heemink等[21]和Dimitriu[22]提出的系统设计方案，模拟过程与分析过程均采用无量纲量，假定在一个30×30的方形区域，污染源在区域里随机分布，流入边界污染物浓度为零，并且外界没有污染物流入。σ=0.2，风速沿x，y方向的分量由公式(8)和公式(9)决定：

(8)

(9)

式中：(xgrid,ygrid)为当前格点坐标；(xcvortex,ycvortex)为对应的漩涡中心点坐标；vel为速度比例系数。在模拟的风速场中，假定气旋中心点位于方形区域的几何中心处，即其坐标为(15,15)。设定vel=0.06。图2为模拟区域内的风矢量图。

图2 风矢量图

通过数值实验验证集合卡尔曼滤波的可行性，使用拉格朗日法将该模型方程离散在30×30网格的风速场上，实验中假定每一格点风速恒定且不随时间变化。在模拟区域中随机放置5个污染源，其坐标为{(5,8),(7,8),(15,5),(3,18),(20,15)}，在每一积分步后各污染源释放污染物浓度的增长量分别为{0.2,0.5,0.4,0.2,0.25}。按照经典数据同化方法研究中的做法，模拟“观测”数据通过在模拟真实数据上添加高斯白噪声获得。假定观测信息符合自回归过程,每一步观测信息由公式(10)给出：

(10)

式中:wj为满足均值为1，方差为0的高斯白噪声；j为观测的位置；λj为对应观测点的衰减系数，每一步后的衰减系数{λ1,…,λj}={0.9,0.8,0.8,0.9,0.9,0.8,0.9,0.9,0.8}。

1.4 性能指标

实验采用均方根误差(RMSE)来评价数据同化效果，RMSE用于评价污染物预报浓度与真实浓度之间的差异，RMSE值越小，表明预报浓度与真实浓度越接近，同化效果越好。

(11)

图3给出了整个同化系统的框图。二维平流扩散模型和集合卡尔曼滤波算法结合可以用来研究不同要求下的同化策略。

图3 二维平流扩散模型同化系统框图

2 结果与讨论

2.1 基本同化实验

基本同化实验将观测点位置随机设置{(4,11),(13,3),(26,15),(15,10),(23,2),(11,9),(15,20),(23,10),(5,25)}，进行不同集合数下的同化实验。基本同化实验假设模式无误差，背景场通过模型向前演进获得，同化的时间窗口长度设计为10。图4给出了第100步时浓度真值等值线图，以图4为基本参照，图5给出了选择20个集合数时基本的同化结果。

图4 第100步时浓度真值等值线图

图5 第100步时同化效果图

分析比较图5中“同化后估计值”曲线和“未同化估计值”曲线，可看出“同化后估计值”曲线与“未同化估计值”曲线相比更接近“真值”曲线，表明经同化方法改善的模型预报效果优于未同化的模型预报效果。此外，图中实线和虚线的大致走向相同，说明集合卡尔曼滤波结合观测信息可以用来研究海洋或者大气污染物扩散问题。图中实线与虚线距离越近，则预测值与真值之间的误差越小，即同化效果越好。

2.2 同化策略研究

不同敏感参数下的RMSE变化情况见图6。

图6 不同敏感参数下的RMSE变化情况

从图6可以看出，①观测方差小于0.1时对同化效果的影响明显，当观测方差取值大于0.05，RMSE值显著增大。②同化窗口长度的取值为1～20时，RMSE值为0.58～0.89，变化幅度较小，同化窗口长度为10时，RMSE值为0.59，大于10时，RMSE值显著增大。③随着集合卡尔曼滤波集合样本数目的增大，系统RMSE变化明显，集合样本数为2～10时，整体呈下降趋势但局部RMSE异常，集合样本数为10～20范围内变化时RMSE值显著减小，当集合样本数增大到20以后，RMSE值基本没有变化而是始终稳定在0.7左右，集合样本数增大到40时，其RMSE值与集合样本数为20时的RMSE值十分接近，即表明在20个集合样本时系统能达到最好的同化效果。

通常情况下，在大气或海洋数据同化系统中应用集合卡尔曼滤波时，采用100个集合样本时能取得较好的同化效果[23]。在理想实验中，20个样本即可达到较好的数据同化效果，这与污染物同化系统模式自由度相关。通常情况下，样本选取方法合理，所需要的样本容量要远小于模式自由度[24]，理想化的二维平流扩散模型自由度较小，所以实验在20个样本数时就取得了较好的同化效果。

2.3 观测策略研究

为了研究观测位置对同化效果的影响，假设污染源水平分布，研究观测位置与污染源距离不同时的同化效果。

1)将污染源位置理想化的设置为{(5,10),(10,10),(15,10),(20,10),(25,10)}，观测点位置设置为{(4,25),(13,25),(26,25),(15,25),(23,25),(11,25),(5,25)}，集合数为20时其RMSE为0.808，实验效果图如图7所示。

图7 观测点与源点间隔为15的同化效果图

2)将污染源位置设置为{(5,10),(10,10),(15,10),(20,10)(25,10)}，观测点位置设置为{(4,15),(13,15),(26,15),(15,15),(23,15),(11,15),(5,15)}，其RMSE为0.686,实验效果图如图8所示。

图8 观测点与源点间隔为5的同化效果图

3)将污染源位置与观测点位置重合进行实验，实验结果得到的RMSE为0.096，实验结果如图9所示。

图9 观测点与源点重合的同化效果图

通过分析比较图7～图9，可以看出观测点与污染物排放点间距离越小误差就越小，即同化效果就越好，观测点与污染源重合情况下同化效果最好。理想情况下，当观测点与污染源重合时，其他污染源对该观测点的影响远小于与观测点重合污染源的影响，这样该观测点测定的数据更准确，可大大提高同化效果。

2.4 局地化策略研究

污染源位置与观测点位置平行分布时，观测点与污染源位置距离越近，同化效果越好。然而，当污染源随机分布，即使观测位置与污染源完全重合，同样取20个集合，基于集合卡尔曼滤波的同化系统依然得出较大的RMSE。图10给出了同化效果图，RMSE为0.583。

图10 污染源随机分布观测位置与污染源位置重合同化效果图

为了解决这个问题，引入局地化技术，采用权重因子(w)按照同化点与观测点的距离限制其卡尔曼增益的大小。图11给出了随着集合数(N)变化改变w时的同化效果。

研究发现，①不同的w带来不同的同化效果，w可以按照观测点和同化点的距离远近，赋予集合卡尔曼增益不同的权重值，获得更好的同化效果；②随着N的增大，同化效果逐步变优，在同一N值下，呈现出多峰多谷的现象，优选w可以获得较低计算代价下的同化效果。③不同N值情况下的最优w值不同，且随着N值的增大，最优w值也会增大，即对背景误差协方差的依赖越大。④N值较少时滤波容易发散，此时同化对观测场的依赖要远大于对背景场的依赖。

图11 w值变化对同化效果的影响

3 结论

基于二维平流扩散模型，建立了一个理想的海洋或者大气环境污染物监测模拟数据同化系统。该系统可以在模拟风速场中，设计污染源位置，模拟实际污染物的扩散过程。利用集合卡尔曼滤波对模型的状态变量进行更新和校正，分析比较了集合数目变化、观测方差变化和同化窗口长度变化对同化效果的影响。为了研究环境监测中观测位置对于数据同化系统的影响，设计了理想状态的平行分布和重合分布等，讨论了观测信息对于数据同化系统的重要性。得到的主要结论有：①基于集合卡尔曼滤波在模拟的风速场上建立了一个环境污染预报系统，利用该系统进行同化实验后，发现集合卡尔曼滤波融合模拟观测信息对二维平流扩散模型的状态变量能够实时校正，实验结果表明，集合卡尔曼滤波用于环境污染预报在理论上是可行的，在实际环境问题中的研究，还需在此基础上进一步开展。②数据同化方法敏感性实验发现：理想情况下，集合样本数为20时，系统的同化性能较好；同化窗口长度和观测方差越小，误差越小。③实验中调整观测点与污染点之间距离，结果显示观测的位置距离污染源较近的同化效果明显优于距离较远的同化效果，观测点与污染点位置重合，同化效果最好。污染源随机分布时，在同化过程中引入局地化技术，利用w依据同化点和观测点间的距离限制了卡尔曼增益的大小，讨论了不同集合数情况下，权重因子对同化效果的影响。但研究仍存在一定的局限性，主要表现在：①实验中使用的二维平流扩散模型比较简单，并且没有考虑模式误差，观测资料由模型模拟生成，实际问题中结果将会受到一定的影响。②实验是基于模拟的风速场上进行的，虽然与污染物的扩散过程比较接近，但与实际环境相比较还是存在差异，所以存在一定的局限性。进一步研究中需要将观测点间的相关性融入进来，进一步提高集合卡尔曼滤波的同化效果。③污染源与观测点的平行分布过于理想化，与实际情形相差太远，未来研究中，将进一步考虑实际监测工作的需要。

[1] 李小飞，张明军，王圣杰，等.中国空气污染指数变化特征及影响因素分析[J].环境科学，2012，33(6)：1 936-1 943.

[2] 高庆先，师华定，张时煌，等.空气污染对气候变化的影响与反馈研究[J].资源科学，2012，34(8)：1 384-1 391.

[3] 师华定，高庆先，张时煌，等.空气污染对气候变化影响与反馈的研究评述[J].环境科学研究，2012，25(9)：974-980.

[4] 李宏，许建平.资料同化技术的发展及其在海洋科学中的应用[J].海洋通报，2011，30(4)：463-472.

[5] 白晓平，李红，方栋，等.资料同化在空气质量预报中的应用[J].地球科学进展，2007，22(1)：66-73.

[6] 王晓彦，刘冰，李健军.区域环境空气质量预报的一般方法和基本原则[J].中国环境监测，2015，31(1)：134-138.

[7] CONSTANTINESCU E M,SANDU A,CHAI T F,et al.Ensemble-based chemical data assimilation.II: Covariance localization[J].Q J R Meteorol Soc,2007,133:1 245-1 256.

[8] HANEA R G,VELDERS G J M,HEEMINK A.Data assimilation of ground-level ozone in Europe with a Kalman filter and chemistry transport model[J].J Geophys Res,2004,109(D10302):doi:10.1029/2003JD004283.

[9] ELBERN H,STRUNK A,SCHMIDT H,et al.Emission rate and chemical state estimation by 4-dimensional variational inversion[J].Atmos Chem Phys,2007,7(14):3 749-3 769.

[10] 周势俊，宋煜，吴士杰.Kalman滤波法在城市空气污染预报中的应用[J].中国环境监测，2000，16(4)：50-56.

[11] 王辉，刘桂梅，万莉颖.数据同化在海洋生态模型中的应用和研究进展[J].地球科学进展，2007，22(10)：989-996.

[12] HOUTEKAMER P L,MITCHELL H L.A sequential ensemble Kalman filter for atmospheric data assimilation[J]. Monthly Weather Review,2001,129:123-137.

[13] 刘成思.集合卡尔曼滤波资料同化方案的设计和研究[D].北京：中国气象科学研究院,2005，6-20.

[14] 刘浩.溢油逆时追踪问题的数值摸拟[D].大连：大连海事大学，2003：5-32.

[15] EVENSEN G.The ensemble Kalman filter theoretical formulation and practical implementation[J].Ocean Dynamics,2003,53(4):343-367.

[16] ANDERSON J L.A local least squares framework for ensemble filtering[J].Monthly weather review,2003,131(4):634-642.

[17] SAKOV P,BERTINO L.Relation between two common localisation methods for the EnKF[J].Computational Geosciences,2011,15(2):225-237.

[18] HAMID M.Hydrologic remote sensing and land surface data assimilation[J].Sensors,2008,8:2 986-3 004.

[19] 裴洪平，何金土，王维维.杭州西湖总磷动态变化预测[J].湖泊科学，1996，8(2)：139-143.

[20] 张玮，王国超，刘燃，等.环抱式港池水体交换与改善措施研究[J].水运工程，2013，4：37-41.

[21] HEEMINK A W,VERLAAN M,SEGERS A J.Variance reduced ensemble Kalman filtering[J].Monthly Weather Review,2001,129:1 718-1 728.

[22] DIMITRIU G.A numerical comparative study on data assimilation using Kalman filters[J].Computers and Mathematics with Applications，2008,55:2 247-2 265.

[23] 张学峰，黄大吉，章本照，等.集合数据同化方法的发展与应用概述[J].海洋学研究，2007，25(1)：88-95.

[24] 王金成，李建平.4DSVD分析误差与样本选取方法和样本容量的关系初探[J].气候与环境研究，2010，15(6)：729-742.

Assimilation Research with Ensemble Kalman Filter Using Two-dimensional Advection Diffusion Model

BAI Yulong，You Yuanhong，SHAO Yu，XU Baoxiong

College of Physics and Electrical Engineering, Northwest Normal University, Lanzhou 730070, China

A simulation data assimilation system aiming at pollutant concentration was built for further improving environmental numerical predication level. Ensemble Kalman filter had been applied to real-time adjust state variables of the two-dimensional advection diffusion model, which could simulate and forecast pollutant concentration. Sensitivity tests were conducted with the change of ensemble numbers, observation variance and assimilation window’s length. Comparing the forecast effects by the different distances between the locations of observation point with pollution source, the assimilation strategies under optimum conditions were discussed, and then a novel method to adjust Kalman gain weight by distances was presented. In addition, the proposed methods could reduce computational cost and obtain better assimilation effect with the smaller ensemble sizes.

data assimilation；advection diffusion model；ensemble Kalman filter

2015-09-22；

2015-11-03

国家自然科学基金(41461078)；兰州市科技计划项目(2015-3-34)

摆玉龙(1973-)，男，甘肃会宁人，博士，教授。

X830.3

A

1002-6002(2016)06- 0123- 07

10.19316/j.issn.1002-6002.2016.06.20