1 引言

资料同化方法经过多年的发展，从基础的松弛逼近方法和逐步订正法为代表的经验资料同化方法，发展到最优插值、三维变分同化和四维变分同化(Cao et al, 2011, 2015; Zhu, et al, 2014)为代表的统计资料同化方法，以及现在的以集合卡尔曼滤波和混合资料同化方法为代表的集合资料同化方法。在资料同化方法发展的过程中，变分资料同化方法和集合卡尔曼滤波方法均包含有线性假设和高斯假设，粒子滤波不需要线性和高斯分布假设，是一种相对先进的资料同化方法。

粒子滤波采用蒙特-卡洛方法采样的样本集合来表示和计算后验概率分布，摒除了线性化模式和线性化观测算子的假设，同时也不需要先验概率的高斯分布假设，粒子滤波可以模拟非高斯分布的后验概率函数分布，这成为粒子滤波最大的优势。粒子滤波第一次由Doucet等(2001)引入并成功应用到低维系统当中。由于地球物理模式或者大气模式的状态量的维数大致在10⁹维，受现在超级计算机的计算和存储能力限制，无法使用大规模的粒子集合，因此，粒子集合的成员数量级约为10—100。Snyder等(2008, 2015)和Snyder(2011)研究表明，粒子滤波粒子数目随着模式维数的增加而呈指数增长。粒子滤波通常分为两步：预报步，模式在时间上的演进；分析步(权重步)，给出集合中每个粒子的相对权重。在大多数情况下，粒子滤波的粒子权重和每个粒子的似然是成比例的。因此，越接近所有观测的粒子，其所获得的权重越大。在粒子滤波模式不断向前演进的情况下，大多数的粒子会渐渐远离观测，以至于其粒子权重变得越来越小直至可以忽略。在整个粒子滤波运行过程的最后阶段，绝大多数的权重都会集中到一个粒子之上，使得这个粒子的权重接近于1，即所谓的“维数灾难(Curse of Dimensionality)”(Doucet et al, 2000, 2001；van Leeuwen, 2009)。

在高维系统中，粒子滤波采用一个有限大小粒子集合来模拟后验概率分布，为了解决“维数灾难”问题，对于粒子权重本身构成的集合的方差取最小值，此时粒子权重在权重均值周围的波动最小，权重之间的差异也就最小，最优提议密度(Doucet, et al, 2000)的思想来源于此。Bocquet等(2010)指出，最优提议密度在小的集合成员数情况下，在Lorenz96系统中的表现要优于连续重要性重采样方法(SIR)。提议密度通过从接近所有观测的高概率区域采样，并且通过引导粒子向观测区域演进，部分解决了“维数灾难”问题。Chorin等(2009)引入了隐式采样，这个方法在复杂状态空间和原始状态空间之间建立了隐式映射关系。Ades等(2013)研究发现，当每个时间步都存在观测时，隐式粒子滤波会衰退成为最优提议密度方法。

隐式等权重粒子滤波(Implicit Equal Weights Particle Filter, IEWPF)(Zhu, et al, 2016)，结合了提议密度和隐式采样的两种方法的优点，通过假设一个从复杂提议密度空间到标准高斯分布的隐式映射关系，将一个非常复杂的等权重问题简化成为一个求解伸缩参数的问题。IEWPF方法在1000维的线性模型和Lorenz96模型上进行了测试，两个系统的试验结果显示了非常好的性质和一致性。每一个粒子都被引导至最差表现的粒子权重的区域，即目标权重被设置成表现最差的粒子的权重值，使得每个粒子的权重均等于目标权重。通过使得粒子提议密度依赖于前一时刻的所有粒子，IEWPF解决了困扰粒子滤波领域多年的“维数灾难”问题。

IEWPF最终目标是进行业务化应用，但是这个过程是一个循序渐进的进程，Zhu等(2016)展示了IEWPF的理论框架和在Lorenz96系统中的应用特性。文中将研究IEWPF在较为复杂的准地转模式中的基础原理和所表现的特性。准地转理论由Charney(1948)在动力学分析中首先应用。随后，Burger(1958)、Phillips (1963)等发展和扩充了准地转理论，给出了多种准地转近似运动的原始方程组。准地转模式显示了静力平衡和地转平衡限制并以逼近真实的方式简化了大气运动，其提供了一个简单的学习框架可用以理解和诊断二维或者三维的天气尺度的大气或者海洋运动。准地转运动可帮助理解质量场(通过横向温度平流)和动量场(通过水平涡度平流)交互产生垂直循环，进而产生真实的天气尺度大气或者海洋模式。准地转运动提供了垂直运动的强迫场以及和中尺度热带气旋相关的云/降水模式的物理视角。

本研究将针对IEWPF的基础理论及其在二层高维准地转模式中的应用特性进行描述，IEWPF通过构建后验概率分布函数的采样样本，使得所有的粒子在同化观测后均具有相等的粒子权重。因此，几乎目前所有粒子滤波方法中使用的重采样就变成了不必要的步骤。更重要的是，本研究将会详细展示在高维系统中粒子权重计算是如何非常明显地被简化的。IEWPF的提出为解决困扰粒子滤波领域多年的“顽疾”——“维数灾难”问题提供了一个非常有效的途径，同时也提供了一类崭新的解决高维粒子滤波资料同化问题的思路和方法。

2 IEWPF基础理论 2.1 分析时间步

粒子滤波是两个理论的结合：贝叶斯理论和蒙特-卡洛方法。贝叶斯理论显示了后验概率密度函数p(x|y)可以通过先验概率p(x)和似然p(y|x)的乘积来表示，观测的概率密度函数p(y)可以被认为是一个正则化常数

(1)

蒙特-卡洛方法是指一个概率密度函数可以通过一个随机采样点或者一个粒子集合来模拟，这些随机采样点或者粒子可以通过δ函数来代表指定的粒子状态

(2)

式(1)和(2)的结合使得后验概率密度函数可以通过特定粒子状态的权重和来进行表示

(3)

式中，x为模式状态，y为观测，p(y|x)作为似然，代表着给定的模式状态和观测的相似程度，p(x)为先验概率，p(x|y)为给定观测后的后验概率密度分布。粒子权重和似然成正比，并且决定于粒子本身和观测的远近程度。在有限的集合成员情况下，可能不会存在足够的粒子和观测非常接近，因此，造成了滤波衰退问题。最后，只剩下唯一一个粒子，却贡献了绝大多数的权重值，使其相对权重接近于1，这意味着由一个粒子来代表后验概率密度函数，从而失去了粒子滤波的统计学意义，这就是前文提到的“维数灾难”。

将粒子滤波理论扩展到一个完整的时间序列，其中模式状态序列为{x⁰，x¹，…，xⁿ}，观测序列为{y¹，y²，…，yⁿ}。这样，后验概率密度函数可以表示为

(4)

提议密度是可以同时乘到式(4)上面的一个因子，在IEWPF理论中选择提议密度为q(xⁿ|X^n-1，yⁿ)。在q(xⁿ|X^n-1，yⁿ)中X^n-1表示来自于前一个时间步的所有粒子的一个集合。那么，模式状态从原始的概率密度函数p(xⁿ|x_i^n-1)中采样，便转变成了从提议密度函数中采样。这样后验概率密度函数可以表示为

(5)

式中，w_i为粒子权重

(6)

从提议密度函数q(xⁿ|X^n-1，yⁿ)中进行采样是IEWPF中的隐式部分。在IEWPF中，并没有从原始的提议密度函数q(xⁿ|X^n-1，yⁿ)中进行采样，而是从一个简单的标准高斯分布提议密度q(ξ)中进行采样(Chorin, et al, 2009)。通过随机映射，这两个概率密度分布函数的关系为

(7)

式中，表示变量变换x_i=g(ξ_i)的雅可比矩阵的行列式的绝对值。

在IEWPF中，定义函数g(·)为

(8)

式中，x_i^a是提议密度函数q(x_iⁿ|X^n-1，yⁿ)的取峰值点，P是概率密度函数的协方差矩阵，α_i是一个标量参数。IEWPF中α_i值的选取使得所有的粒子可以得到相等的粒子权重，称这个相等的粒子权重为目标粒子权重w_target。因此，由式(9)求取每个粒子的α_i值

(9)

式(9)为IEWPF的思想核心，是IEWPF等权重的基本原理部分。该公式保证了IEWPF在拥有很大数量独立观测的高维系统中，集合粒子权重不会产生粒子衰退现象，从而解决了困扰粒子滤波多年的“维数灾难”问题。

令J_i和J_i^prev分别代表分析时刻和分析时刻前一时刻的粒子权重的两倍对数值，则粒子的权重公式可以写为

(10)

在观测误差和模式误差为高斯分布的前提假设下，利用Sylvester's行列式定理可以将公式简化为

(11)

在高维假设下，经过化简(Zhu, et al, 2016)，可得

(12)

式中，存在如下的关系：

(13)

通过隐式提议密度和目标权重的设置，可以保证所有的粒子的权重是相等的，进而解决了“维数灾难”问题。因此，修改后的粒子i的模式状态可以表示为

(14)

式中，K=QH^T(HQH^T+R)^-1，P=(Q^-1+H^TR^-1H)^-1，α_i为一个伸缩参数，可以将粒子权重伸缩至目标权重。Q是模式误差协方差矩阵，H是观测算子，R是观测误差协方差矩阵。f(·)为n-1时刻至n时刻的正演模式。在高维系统模式的前提假设下，α_i的解析解可由兰伯特W函数求得(Zhu, et al, 2016)

(15)

式中，N_x为模式状态量的维数，c=C-ϕ_i-J_i^prev，C=-2lgw_target为常数，取之等于目标权重的两倍对数值，α_i的解析解存在特定的结构和特性(Zhu, et al, 2016)。

2.2 松弛逼近项

在两个观测时刻之间会存在若干个模式时间步，通常情况下可以扩展上面描述的IEWPF框架至这一系列的时间步中，因此x_i^a就变成了一条模式轨迹。在本试验中的非观测时间步，采用较为简单的松弛逼近方法对粒子集合进行处理，并指定在观测时间步之前的非观测时间步的收敛到观测位置的松弛提议密度函数为

(16)

式中，B(τ)[yⁿ-hM(x_i^j-1)]为松弛项，在时间步n迫使模式状态收敛至观测位置。为修改后模式的随机强迫量的协方差矩阵，简化处理选择该协方差矩阵等于Q矩阵。从提议密度中进行采样代替从原始模式方程中进行采样，可得

(17)

式中，。

松弛强度控制量B(τ)可表示为

(18)

式中，τ随着时间步的增加而从0至1呈线性增加，b是一个可以定义的常数，B(τ)控制着松弛收敛的强度。式(17)的表达式在每一个时间步通过简单的p/q乘以粒子权重得到w_i^prev。

3 试验及分析 3.1 两层准地转模式

采用一个高维两层准地转模式做理想试验，以研究新资料同化方法的基本特性。其中，资料同化系统代码采用的是MatLab代码，而模式系统代码采用的是Fortran代码，这样整个系统兼具了速度和便捷性，同时保证了精度和守恒性。两层准地转模式的方程(Fandry, et al, 1984; Pedlosky, 1979)采用给定的无量纲变量表示为

(19)

式中，D/Dt表示个别微商(全导数，对应个别变化)，q₁和q₂分别表示上下两层的准地转位势涡度，其中下标1表示上层

(20)

式中，β为科里奥利参数(无量纲)的北向微分值，R_s表示地形强迫或者热强迫(净的热收支)。模式区域纬向采用循环方式，并且假设经向速度在距离区域南部和北部边界一个网格点的空间距离处消失。

无量纲化是一种标准的表示形式，为了方程的完整性给出无量纲化的过程。定义一个典型的长度量L，一个典型的速度量U，上、下两层的厚度定义为D₁和D₂，南部边界的科里奥利参数为f₀，其北向微分值为β₀，重力加速度为g，两层之间的位温差异为Δθ，平均位温为θ。变量头上的波浪号表示量纲的时间、空间坐标和速度，则有

(21)

式中，罗斯贝数为β₀=U/(f₀L)。

模式的预报量是定义在方形网格[n_x×n_y]上的流函数，格点数目的索引从南至北，从西至东顺序增大。该模式的时间步进算法的设计以运行速度而不是精确度为目的，其仅对于时间步长Δt的一阶导数是精确的。这样设计出于实用性的考虑，在给定仅一个单独的同一时间尺度上的信息情况下可以运行一个时间步。通常情况下，模式运行一个时间步，流函数从时刻t积分到时刻t+Δt。然而，为了可以在分析层中获得风和位势涡度从而可以同化风的观测，对时间步的离散方法进行了特殊处理，鉴于本研究的目的是测试新的IEWPF方法，因此，对时间步的离散处理不再赘述。

3.2 集合初始化方案和参数设置

Zhu等(2016)的高维Lorenz96试验中采用的集合初始化方案为从背景误差协方差矩阵B中进行采样。这种集合初始化方案的优点是简单易行，针对简易的理想混沌模式试验是可行的，但是对于模拟真实大尺度大气或者海洋模式的准地转模式或者其他的一般大气/海洋环流模式并不适用，因为其背景误差协方差矩阵并不是真实的背景误差协方差矩阵，会随着一些因素(比如季节、日间变化等)发生变化。因此，产生的随机格点上的资料同化状态量也非常有可能不满足模式系统平衡条件，极易造成模式系统运行的崩溃。因此，文中选用了另外一种集合成员初始化方案：运行一段较长时间的准地转模式，选取运行过程中的随机时间点的模式状态量，作为IEWPF的初始集合成员。

准地转模式采用确定性模式。本节设计的试验为资料同化试验，进行资料同化的变量有两个，分别是纬向风u和经向风v。由于风速和风向是直接观测的变量，其与纬向风和经向风只需要一个简单转换关系，为了简化观测预处理过程，采用直接对纬向风和经向风进行观测。试验中的参数设置分别为：观测误差协方差矩阵为对角矩阵，其主对角元素为R=0.8，观测算子矩阵H为单位矩阵，其主对角元素为1，模式误差协方差矩阵也是对角矩阵，在试验中设定的模式误差协方差矩阵Q主对角元素的值为0.2。所有资料同化试验采用每5个时间步为一个同化周期，即每隔4个非观测时间步有一个观测时间步(同化分析时间步)。非观测时间步(非同化分析步)采用松弛逼近框架，其参数设置为b=0.15，τ=0.3。由于计算资源和代码开发周期的限制，采用MatLab和Fortran的混合编程方式。由于计算资源有限，采用占用计算资源和存储资源较少的对角矩阵作为模式误差协方差矩阵，值得注意的是，实际的模式误差各元素之间应该是存在相关性的。

模式水平格点数n_x=128，垂直格点数n_y=64，模式的纬向区域长度为12×10⁶ m，模式的经向区域长度为6.3×10⁶ m。模式分为两层，高度分别为6000和4000 m，进行资料同化时，模式高度分别设置为5500和4500 m。模式的时间步长设置为Δt=0.036 s，南部边界的科里奥利参数f₀=10^-4，f的经向梯度β₀=1.5×10^-11。边界条件的处理在模式的运行过程中非常重要，其涉及到模式和外界的能量和质量交换，影响到整个模式的运行效果，所以模式运行过程中设计了边界条件更新条件。

粒子滤波集合成员数设置为20，每个粒子(集合成员)的位势涡度和流函数的初始化状态如图 1所示。

3.3 水平动能谱

水平动能谱是在每一层上对风场(u，v)进行二维离散余弦变换得到。由于离散余弦变换可以有效避免离散傅里叶变换在小尺度范围内的虚假上翘(Peng, et al, 2014)，得到更宽的谱，因此，该方法适用于有限区域的物理量的谱展开，并且优于基于二维离散傅里叶变换进行的谱分析，其计算方法如下：

记任意变量q的二维离散余弦变换为，其中k=(k_x, k_y)表示水平波数矢量，k_x= ，k_y= ,N_i、N_j分别表示纬向和经向格点总数，Δ表示格点距离，m=1, 2, 3, …，N_i-1, n=1, 2, 3, …，N_j-1, 总的水平波数为。这样在给定的时刻t，给定的高度z上，单位质量的水平动能谱为

(22)

一维水平波数k_h谱的构造是通过在k_x-k_y平面上等|k|圆环截断求和得到

(23)

式中，，l=1，2，3，…，N-1，N=min(N_i，N_j)。

3.4 资料同化试验 3.4.1 资料同化效果

观测系统试验分别截取第11个时间步和第96个时间步的位势涡度以及流函数的真值场、IEWPF同化分析场以及等权重粒子滤波(Equivalent Weights Particle Filter, EWPF，见3.4.3节)同化分析场(图 2a—f)。虽然进行同化的变量并不是直接的位势涡度和流函数，而是纬向和经向风，但从两个时刻的对比可以看出，IEWPF资料同化框架可以很好地对不同的初始场进行资料同化，同化分析场和真值场的主要特征完全相同，只是在位势涡度和流函数的强度上有一定的减弱。这表明，IEWPF对于具有复杂关系的多变量模式，仍然能够在后验概率密度函数的高概率区采样。

从图 3a、b显示的第11个时间步、第96个时间步的位势涡度及流函数的同化分析场与真实场的差值可以看出，同化分析场和真值场的差异非常小，而且在区域内分布比较平均，在边界处相对略大一些，并且随着资料同化的时间增长，两者的差值越来越小。

3.4.2 水平动能谱

水平动能谱和一维水平波数谱是评判模式的一个标准，同化前后水平动能谱和一维水平波数谱的变化可以反映出同化系统本身对于变量平衡性的一种度量标准。

观测的生成是在真值场从模式状态空间映射到观测空间之后，添加分布为n(0, R)的高斯噪声得到，可以看出在图 4b和5b中，噪声的分布影响了整个区域的能量谱的分布。由于高斯噪声属于全频域范围内的影响，因此在图 4c和d以及图 5c和d的对比中，可以明显看出高斯噪声对水平动能谱和一维水平波数谱的影响，主要差异在于加在模拟观测之上的高斯噪声影响了频谱的高频部分的分布。

3.4.3 和EWPF方法进行对比

EWPF (Ades et al, 2013, 2015a, 2015b; van Leeuwen, 2010, 2011)同样利用提议密度引导粒子成员演进到后验概率密度函数的高概率区域，和IEWPF的基础思想相似。EWPF已经经过了真实的海洋气候模式HadCM3试验验证，是一种有效的资料同化方法(Browne, et al, 2015)。因此，EWPF和IEWPF作为两种资料同化方法进行比较。

EWPF试验中，采用和IEWPF完全相同的参数设置和收敛松弛项，不采用重采样步骤，保留100%的粒子成员。EWPF资料同化分析场结果如图 2c和f所示。经过对比发现，两者的资料同化分析场非常接近，因此两者的资料同化效果对于两层准地转模式是近乎相同的。但是，从资料同化系统运行的时间上，运行相同时间步的试验，IEWPF所花的时间为628.4426 s，而EWPF运行的时间为937.0624 s，可见在运行效率上，IEWPF要远好于EWPF。

4 结论

IEWPF方法采用等权重隐式采样方法，使得粒子成员向着观测集中的位置运动，使得所有粒子成员的权重等于目标权重，从而解决了“维数灾难”问题。论文通过在较大规模的两层准地转模式上验证了IEWPF的一致性和优异性。对于静力平衡的两层准地转模式，IEWPF通过资料同化得到的分析场很好地得到了真实场的主要特征，并且随着时间的增加推移，真实场和同化分析场之间的差值逐渐减小，其水平能量谱也和真实场的水平能量谱相似，其主要差异在于加在模拟观测之上的高斯噪声影响了频谱的高频部分的分布。

IEWPF方法为高维非线性资料同化问题提供了有效的解决方案，经过文中的试验以及与EWPF的对比试验证明了该资料同化方法的高效及可靠性。同时，IEWPF也为资料同化领域的其他问题提供了优秀的解决方案，比如利用资料同化方法进行参数估计等。当然，IEWPF距离实际应用仍然存在许多问题，如粒子的初始化、参数的选取、并行算法的设计等。而最主要的问题在于，IEWPF需要模式误差协方差矩阵构成移动粒子的增益矩阵，但目前国际上对于模式误差本身知之甚少，因此针对模式误差的研究将作为后续研究的主要内容。

致谢: 感谢英国雷丁大学Peter Jan van Leeuwen教授的悉心指导，感谢国防科大海洋科学与工程研究院吴建平研究员、彭军博士关于动能谱的讨论。