1. 引言

利用资料同化的方法提取观测资料的有效信息是改进数值预报初始场，进而提高数值预报水平的一个重要手段 [1]。在变分同化下，资料同化系统将观测资料和背景信息进行优化组合，为数值天气预报系统提供大气初始状态的最优估计值 [2]，就需要定义背景误差协方差矩阵B和观测误差协方差矩阵R。承载着背景场误差信息的B矩阵是资料同化成败的关键因素之一 [3]，合理估计B矩阵的量值和结构是提高同化分析和数值预报水平的基本前提 [4]。背景误差协方差的估计方法有NMC (National Meteorology Center)方法 [5]、Monte Carlo方法 [6] 以及集合资料同化(EDA, Ensemble Data Assimilation)方法 [7] [8] [9] [10] 等。通常构建均匀、各向同性的背景误差协方差矩阵B，作为B矩阵真实值的一个替代品 [11]，这种假设能够简化运算，但无法反映大气的时空变化特征 [12]。为了使背景误差协方差矩阵B具有流依赖(flow-dependent)性质(即能够反映天气形势的演变)，采用EDA方法进行估计。

集合四维变分资料同化(Ensemble Four-Dimensional Variational Data Assimilation, En4DVAR) [13] [14] [15] [16] [17] 能在充分利用变分同化方法的优势的基础上，估计出流依赖背景误差协方差，是解决资料同化中背景误差协方差矩阵静态不变问题的一种有效途径。在实际的变分同化系统中，一般会将分析变量划分为与平衡变量相关的平衡部分和剩余的非平衡部分。对于平衡部分，较为合理的做法是只引入其中一种变量的流依赖信息，通过平衡约束关系将流依赖信息传播到其他变量上。如ECMWF (European Centre For Medium Range Weather Forecasts，欧洲中期天气预报中心)通过更新涡度的流依赖背景误差，引入流依赖信息。ECMWF已在Cycle37R2版本(2011.5)采用EDA方法估计涡度、湿度、臭氧变量的流依赖背景误差方差 [18]，通过平衡算子M、N、P将流依赖信息传播到散度、温度和地面气压的平衡项部分(η_b, (T,p_s)_b)，GRAPES-4DVAR [19] [20] [21] 则采用流函数作为平衡变量。对于非平衡部分，由于此部分与其他变量无关，需要估计出其他所有变量的非平衡部分的流依赖背景误差信息。陈耀登等 [22] [23] 使用WRFDA (Weather Research and Forecasting Data Assimilation) [24] [25] [26] 通过分析非平衡变量(速度势、温度等)的分布特点，研究背景误差协方差的特征以及背景误差协方差对台风同化和预报的影响。

已有研究表明，随着模式分辨率的日益提升，非平衡方差在总方差的比重已显著增加 [27] [28]。因此如何估计和应用流依赖的非平衡项方差，实现真正全流依赖的背景误差协方差，对提升流依赖背景误差信息的应用水平具有十分重要的研究意义。

2. 实验

2.1. 背景误差协方差模型

YH4DVAR [29] [30] 是我国自主研发的全球四维变分资料同化系统，使用基于球面小波的背景误差协方差模型 [8] [31] [32]。在变分资料同化中，为了估计大气系统的状态，需要寻找合适的分析场使得目标泛函最小化。3DVAR (Three-Dimensional Variational Data Assimilation，三维变分资料同化) [5] [33] [34] [35] [36] [37] 的目标泛函表示如下：

$J\left(x\right)=J_{\text{b}}+J_{\text{o}}=\frac{1}{2}{\left(x-x_{\text{b}}\right)}^{\text{T}}{B}^{-1}\left(x-x_{\text{b}}\right)+\frac{1}{2}{\left(y-Hx\right)}^{\text{T}}{R}^{-1}\left(y-Hx\right)$(1)

其中x_b为背景场，x为模式变量，B为背景误差协方差矩阵，y是观测向量，H为观测算子，R为观测误差协方差矩阵。

使用控制变量变换方法 [38]，引入算子L，将x转化为χ，算子L满足LL^T= B，从而对B的定义转化为对L的定义。则目标泛函背景项写作：

$\begin{array}{c}{J}_b\left(x\right)=\frac{1}{2}{\left(x-x_{\text{b}}\right)}^{\text{T}}{B}^{-1}\left(x-x_{\text{b}}\right)=\frac{1}{2}{\left(L\chi \right)}^{\text{T}}{B}^{-1}\left(L\chi \right)\\ =\frac{1}{2}{\chi }^{\text{T}}{L}^{\text{T}}{B}^{-1}L\chi =\frac{1}{2}{\chi }^{\text{T}}{L}^{\text{T}}{\left(L{L}^{\text{T}}\right)}^{-1}L\chi =\frac{1}{2}{\chi }^{\text{T}}\chi \end{array}$(2)

使用涡度(ζ)、散度(η)、温度(T)、表面气压(p_s)、比湿(q)，作为模式变量。涡度被视为平衡变量(即自身不存在非平衡项)，散度、温度和表面气压依据平衡关系能够被分解为平衡项和非平衡项。比湿是独立的变量。将模式变量(ζ,η, (T,p_s),q)转化控制变量(ζ,η_u, (T,p_s)_u,q)。

将L表示为两个算子的乘积：

$L=K{B}_{\text{u}}^{1/2}$(3)

其中K是将非平衡变量转化为全变量的平衡算子。对称矩阵 $B_{\text{u}}^{1/2}$是非平衡变量背景误差协方差矩阵B_u的平方根，即：

$B_{\text{u}}^{\text{T}/2}{B}_{\text{u}}^{1/2}=B_{\text{u}}$(4)

通常需要对平衡算子K和非平衡变量背景误差协方差矩阵B_u进行简化处理，如假定B_u是一个块对角矩阵，变量之间不存在相关性。B_u的矩阵形式为：

$B_{\text{u}}=\left[\begin{array}{cccc}C\left(\zeta \right)& 0& 0& 0\\ 0& C\left({\eta }_{\text{u}}\right)& 0& 0\\ 0& 0& C{\left(T,p_{\text{b}}\right)}_{\text{u}}& 0\\ 0& 0& 0& C\left(q\right)\end{array}\right]$(5)

其中“C( )”表示每个变量的自协方差矩阵。引入算子N、M、P，全变量与非平衡量之间的关系为：

$\{\begin{array}{l}\eta =M\zeta +{\eta }_{\text{u}}\\ \left(T,p_{\text{s}}\right)=N\zeta +P{\eta }_{\text{u}}+{\left(T,p_{\text{s}}\right)}_{\text{b}}\end{array}$(6)

其中下标“u”表示变量的非平衡部分。根据关系(2.4)和(2.5)得到K的矩阵形式为：

$K=\left[\begin{array}{cccc}I& 0& 0& 0\\ M& I& 0& 0\\ N& P& I& 0\\ 0& 0& 0& I\end{array}\right]$(7)

最终背景误差协方差的矩阵形式为：

$B=L{L}^{\text{T}}=K{B}_{\text{u}}{K}^{\text{T}}=\left[\begin{array}{cccc}C\left(\zeta \right)& C\left(\zeta \right)M^{\text{T}}& C\left(\zeta \right)N^{\text{T}}& 0\\ MC\left(\zeta \right)& MC\left(\zeta \right)M^{\text{T}}+C\left({\eta }_{\text{u}}\right)& MC\left(\zeta \right)N^{\text{T}}+C\left({\eta }_{\text{u}}\right)P^{\text{T}}& 0\\ NC\left(\zeta \right)& NC\left(\zeta \right)M^{\text{T}}+PC\left({\eta }_{\text{u}}\right)& NC\left(\zeta \right)N^{\text{T}}+PC\left({\eta }_{\text{u}}\right)P^{\text{T}}+C\left(T,p_{\text{s}}\right)& 0\\ 0& 0& 0& C\left(q\right)\end{array}\right]$(8)

对非平衡项，每个自相关矩阵具有相同的形式：

$C=S^{\text{T}}{V}^{\text{T}/2}{S}^{\text{T}}ED{E}^{\text{T}}{S}^{-1}{V}^{1/2}S$(9)

其中，S表示谱到格点的变换，V是模式高斯网格上非平衡变量的方差对角矩阵，E是特征向量矩阵，D是格点空间中方差归一化所得的非平衡变量在谱空间中的垂直协方差矩阵的特征值构成的对角矩阵。

算子M、N又可划分为水平平衡算子和垂直平衡算子的乘积：

$\{\begin{array}{l}M=MH\\ N=NH\end{array}$(10)

其中H是块对角阵水平算子，作用是将涡度的谱系数转化为中间量P_b，P_b为质量场，定义为 [6] [27]：

$P_{\text{b}}\left(n,m\right)={\beta }_1\left(n,m\right)\zeta \left(n,m+1\right)+{\beta }_2\left(n,m\right)\zeta \left(n,m-1\right)$(11)

其中n为总波数，m为纬向波数，β₁、β₂是H的水平平衡系数。H中的水平算子具有与球面上标准分析的线性平衡算子相同结构，使得谱空间中J_b结构随纬度变化。M、N和P算子具有相同的结构：块对角、一个谱分量对应一个全垂直矩阵，该矩阵仅与总波数n有关。M、N、P矩阵、β₁、β₂矩阵采用NMC或集合方法计算得到。NMC方法是使用相同时刻不同预报时效的预报场之间的差异统计得到估计结果。集合方法的统计样本既可以是样本之间的差也可以是样本与样本集合平均之间的差。

2.2. 非平衡项方差的流依赖估计及后处理

采用集合四维变分资料同化系统得到用于统计流依赖背景误差方差的集合样本，系统架构见图1。

图1给出了集合四维变分资料同化的系统架构。在EDA循环中，输入扰动的背景场、观测场和SST (Sea Surface Temperature，海表面温度)，得到扰动分析，经过预报步骤后得到扰动预报量。在方差后处理过程中，原始方差经过方差尺度调整和方差滤波之后，得到EDA方差。在4DVAR循环过程中，从“方差后处理”过程中得到的EDA方差与背景场成为输入量，得到分析场，经过预报步骤得到预报场。

为了表征模式中的不确定性，模式积分过程中采用SPPT (Stochastic Perturbed Parameterized Tendency，参数化倾向随机扰动) [39] 方案现实背景场的隐式扰动。非平衡项方差为总方差减去方差的平衡部分，为了减少非平衡项方差估计值中的随机误差项，采用了谱滤波器 [40]：

$\rho \left(n\right)={\left[\cos \left(0.5\pi \frac{\min \left(n,N_{\text{trunc}}\right)}{N_{\text{trunc}}}\right)\right]}^2$(12)

其中截断波数固定设为1000。

2.3. 效果检验

基于2.1、2.2节的原理和方法，使用数值预报模式进行计算，垂直分层总共137层，第1层为最高层，第137层为最低层，第114层对应850 hPa等压面。

解释方差能够体现假设中的平衡关系反映的程度 [6]，方差解释率的垂直廓线能够较为直观地反映出不同的垂直高度上非平衡项方差对于变量总方差的贡献。解释方差占重比(方差解释率，Explained Variance Ratio) [4] 表示为：

$1-\frac{v_{\text{unbal}}}{v_{\text{total}}}$(13)

其中v_total为总方差，v_unbal为变量的非平衡项方差。

图2给出了温度平衡项方差在总方差中的解释方差占重比，非平衡分量的误差方差对变量的总误差方差有显着贡献。在图1中，模式层60以上，解释方差比重不低于0.3，最大是在模式层20的0.6；低于第80层的高度，解释方差比重小于0.1。可见方差的平衡分量对总误差方差贡献，高层明显大于低层。

图3给出诊断误差结构的一个例子，其中显示了温度非平衡分量的纬度和月时间平均垂直剖面，可以看出在各个垂直层上，热带地区的非平衡误差往往相对较大。在垂直方向非平衡误差中，模式层60 (约100 hPa)附近的上对流层与下平流层(UTLS, Upper Troposphere Lower Stratosphere)，在靠近平层顶层(模式层10到20，约1~3 hPa)以及行星边界层(模式层115至137，约850 hPa至表面)中具有极大值，并且在模式层底层60˚S~80˚S的区域出现最大值。在这些区域，4DVAR中的质量—风平衡关系做不到对总变量的误差的完全解释。

图4显示非平衡表面气压误差的月平均值。黑色等值线为平均海平面气压，蓝色区域为非平衡表面气压误差较大区域。高原山地地形(例如青藏高原、落基山脉、安第斯山脉等)，其非平衡部分误差较大，可能是受地形影响，本身不存在实际的海平面气压。在热带地区尤其是赤道附近，非平衡误差相对较大。在南半球高纬度地区的低压系统周围，也存在着非平衡误差相对较大的区域。总体看来，非平衡表面气压误差较大值更倾向于沿着纬圈分布(例如0˚N、40˚S、70˚S)，同时会受到地形的影响。

选取第114层(对应850 hPa等气压面)，对温度非平衡项方差及其校正、滤波的结果的全球分布进行分析。图5显示温度的非平衡项方差的全球分布，能够看到温度的非平衡项方差的高值区域主要位于南、北半球的中高纬度，形状呈条带状，并且更多分布海洋上空。在北大西洋中部、新西兰东部、非洲好望角、南美洲西海岸和南太平洋东部附近的海域，温度的非平衡项方差最大。图6显示的是校正后温度的非平衡项方差的全球分布。经过校正，消除系统误差的影响。图7显示的是校正和滤波后的温度的非平衡项方差的全球分布，通过滤波消除掉噪声，起到一定平滑作用。

3. 结论和讨论

3.1. 结论

本文简要介绍了YH4DVAR的背景误差协方差模型和集合四维变分资料同化系统架构，重点分析了平衡算子，通过集合方法估计得到了散度、温度和表面气压的非平衡项流依赖方差，为了减少样本对方差估计的影响，对非平衡项方差进行了校正和滤波处理，经过校正消除系统误差，再经过滤波消除噪声，保留有价值的信息。主要工作及结论总结如下：

1) 使用EDA方法，通过计算解释方差分析，对非平衡项方差对总方差的贡献，对月平均的非平衡温度EDA标准差的经向平均的垂直剖面以及非平衡项表面气压的标准差进行分析，看出方差解释率在第80层以下空间区域的远小于第80层以上，在各个垂直层上，热带地区的非平衡误差往往相对较大，而表面气压非平衡项标准差高值区域呈条带状，大致沿着纬圈分布，也存在于南半球高纬度地区的低压系统周围。

2) 探究模式第114层(850 hPa)的温度非平衡项方差的分布特点，并作校正、滤波处理。从2017年7月2日的温度非平衡项方差的全球分布看出，温度非平衡项方差具有条带状分布特点，高值聚集于南、北半球中高纬度，并且在海洋上空分布多于陆地。

3.2. 讨论

本文对于非平衡项流依赖方差进行分析，未来可以结合台风、暴雨的个例来检验流依赖非平衡项方差对同化和分析的影响。

致谢

感谢国防科技大学气象海洋学院曹小群、刘柏年老师的指导。

基金项目

国家重点研发计划(批准号：2018YFC1506704)和国家自然科学基金(批准号：41475094)资助的课题。

参考文献