2. 中国气象局数值预报中心,北京 100081
2. CMA Numerical Prediction Centre, Beijing 100081
随着为数值天气预报(numerical weather prediction, NWP)模式提供初始场的复杂性和重要性不断提升,大气资料同化本身已成为一门快速发展的学科(Daley, 1991; Kalnay, 2003; 邹晓蕾, 2009; 朱国富, 2015a;2015b;2015c)。自20世纪90年代以来,变分同化(包括三维变分3DVAR和四维变分4DVAR)逐步成为各先进NWP中心的主流资料同化方案(Rabier, 2005)。在同化方法的不断发展创新中,核心问题之一是如何更好地保证分析场的平衡。如果分析场不平衡,后续模式积分会激发虚假重力波,不仅使得同化引入的观测信息被快速频散,更会导致实际大气演变信息被噪音所掩盖,造成预报失败(Daley, 1991)。
由于分析场的平衡十分重要,人们还专门发展了各种初始化方法对同化分析场做进一步处理,以便为模式提供更加平衡的初值(Lynch et al, 2010)。不过,独立的初始化操作会使得分析场远离观测,造成分析阶段引入的观测信息流失(Bloom et al, 1996)。为减少初始化对分析场的修正,通过在同化框架中构建合理的动力平衡约束,使得分析向观测靠拢的同时保证变量间的平衡显得尤为重要(Parrish et al, 1992)。
从运动学角度来看,“动力平衡”是指作用于空气块或空气质点上的不同作用力相互抵消,空气块加速度为零时的一种状态,例如地转平衡是科氏力与水平气压梯度力相抵消时的平衡状态(Holton, 1992)。如果将动力平衡引入到同化分析中去,约束不同分析变量的协同变化,使资料同化变为一个有约束的最优化问题,那么此时这些动力平衡就充当起了“动力平衡约束”的角色(Daley, 1991)。
在变分同化框架中引入动力平衡约束的途径主要包括三种:一是背景误差协方差,二是平衡弱约束,三是4DVAR中的模式约束。其中背景误差协方差是构建动力平衡约束的基础性手段,几乎所有业务运行的变分同化系统均采用了该方案。平衡弱约束以及模式约束是针对该方案的进一步补充和强化。通过背景误差协方差引入动力平衡约束主要可以带来三方面的好处:(1) 将分析场的相空间限定于缓变模态,滤除难于准确分析的非平衡运动,保证分析场的平衡(Bannister, 2008a);(2) 提高观测的使用效率,弥补观测系统存在的不足,例如可以从质量场观测中提取有用的风场信息(Thepaut et al, 2010);(3) 改善变分极小化问题的性状,是变分预处理的重要环节(Lorenc, 1997; Talagrand, 2010)。值得注意的是,与早期单变量分析结合静力初始化保证分析场平衡的方案(例如,同化系统只分析高度场,而风场由高度场根据平衡关系导出,分析场中不存在非平衡)不同,采用前述背景误差协方差等三种形式引入动力平衡约束时,仍会允许非平衡运动在一定合理范围内出现(Daley, 1991),这在下文通过背景误差协方差构造动力平衡约束的实施细节中可以明显看出。
正因为通过背景误差协方差构造动力平衡约束十分重要,有关这方面的科学研究和业务应用发表了大量的论文。然而,由于不同系统在分析变量选取、离散方案设计以及主要关注区域等问题上存在差异,讨论和侧重的科学问题、技术细节各不相同。为了总结这些具体研究中的共性问题和共同认识,理清该领域研究中尚存在的主要不足以及未来的发展方向,有必要整理综述这方面的研究进展。由于不同天气尺度下的平衡约束差异很大,本文主要关注全球和网格较粗的有限区域同化框架中平衡约束构造的研究进展,而对于刚刚起步的面向对流尺度系统的高分辨率同化框架中的平衡约束构造,本文将在未来研究展望中给出讨论。
1 动力平衡约束的构建 1.1 变分框架背景变分资料同化通过极小化一个目标函数得到所要求解的分析场,这个目标函数一般包含以下两个部分(Lorenc, 1986; Courtier et al, 1994):
$J\left( \mathit{\boldsymbol{x}} \right) = {J_{\rm{b}}}\left( \mathit{\boldsymbol{x}} \right) + {J_{\rm{o}}}\left( \mathit{\boldsymbol{x}} \right)$ | (1) |
这里采用了增量形式的目标函数,x表示分析场相对于参考场的增量。式(1) 中,Jb度量了分析场与背景场之间的距离,Jo度量了分析场与观测场之间的距离,变分同化通过极小化该目标函数得到所需要的分析场。如果采用背景场作为参考场,Jb可以简单表达为:
${J_{\rm{b}}}\left( \mathit{\boldsymbol{x}} \right) = {\mathit{\boldsymbol{x}}^{\rm{T}}}{\mathit{\boldsymbol{B}}^{ - 1}}\mathit{\boldsymbol{x}}$ | (2) |
式中,B是背景误差协方差矩阵,反映了背景场的误差结构。B矩阵在变分资料同化中起着十分重要的作用,是决定资料同化成败的核心因素。
在实际情形中,B矩阵是一个超大矩阵,元素值难以一一确定,并接近病态。因此,需对其进行变分预处理,采用模块化形式表达它。预处理一般采用控制变量变换进行,是各业务NWP中心的通用方案(Parrish et al, 1992; Derber et al, 1999; Bannister, 2008b; 薛纪善等, 2008)。控制变量变换的思路是寻找一组背景误差间完全不相关的新分析变量,称为控制变量w,替代原有分析变量x:
$\mathit{\boldsymbol{x}} = \mathit{\boldsymbol{Uw}}$ | (3) |
式中,算子U满足条件B=UUT,此时式(2) 给出的Jb可以进一步简化为:
${J_{\rm{b}}}\left( \mathit{\boldsymbol{w}} \right) = {\mathit{\boldsymbol{w}}^{\rm{T}}}\mathit{\boldsymbol{w}}$ | (4) |
通过控制变量变换,背景场项的Hessian矩阵由B-1简化为单位阵,条件数大为减小,极小化问题的性状得到了极大改善(Lorenc, 1997)。
在具体实施中,U算子被分成Up、Uh、Uv三个算子逐步实施(Barker et al, 2004; 张华等, 2004)。Up算子称为物理变换,它抽取背景误差之间的交叉协相关,表达不同变量间的动力平衡约束:
$\mathit{\boldsymbol{x}} = {\mathit{\boldsymbol{U}}_{\rm{p}}}{\mathit{\boldsymbol{x}}_{\rm{u}}}$ | (5) |
通过Up变换,xu中不同变量间相互独立,被称为独立分析变量,其对应的背景误差协方差矩阵简化为块对角阵。经Up变换以后,同一变量在不同空间位置上的背景误差自协方差仍然存在,采用水平变换Uh和垂直变换Uv做进一步处理,最终获得满足条件的控制变量w(朱宗申等, 2002; 庄世宇等, 2005)。
通过控制变量变换方案进行预处理,采用U算子模块化构造B矩阵,避免了直接表达和求解B矩阵的困难,而将问题转移到了U算子的准确构造上去。本文主要关注的是物理变换算子Up,它承载着动力平衡约束的构造。下面给出物理变换的基本方案和相关研究进展。
1.2 物理变换的实施在理想情形下,经Up变换得到的xu的不同变量间应完全独立,这在实际框架中难以实现。可行的方案是找一组变量使它们之间的背景误差交叉协相关尽可能小,并忽略剩余的相关。为了达到上述目标,Up变换先选择一个能表征大气主要平衡模态的变量作为xu的第一个独立分析变量,然后依次将其余分析变量投影到已选择的独立分析变量上去,并将其剩余部分作为下一个独立分析变量。
对于中高纬大尺度大气运动而言,最主要的平衡约束为准地转平衡和静力平衡(Holton, 1992)。基于这两个平衡,大气运动具有准水平无辐散的特征。因此,可以选择质量(气压或温度,两者通过静力平衡相互转换)或者旋转风作为第一独立分析变量,表征大气的平衡模态(Lorenc et al, 2003)。不过,旋转风(涡度ζ或流函数ψ)是绝大部分系统的共同选择(见表 1第2列),其合理性可从地转适应理论出发理解。根据该理论,水平尺度大、垂直尺度小的系统,旋转风场向质量场调整,宜选用质量表征平衡模态;而对于水平尺度小、垂直尺度大以及在赤道地区的系统,质量场向旋转风场调整,宜选用旋转风表征平衡模态。针对不同尺度选择不同的变量是难以实现的,但如果假设背景场已经能比较好地表达大尺度的信息,更多需要通过观测改善较小尺度上的信息,那么使用旋转风场作为平衡模态会更具优势(Lorenc et al, 2003)。同样基于地转适应理论,在不同尺度下均保持守恒的变量——位涡应是平衡分析变量的理想候选者,但由于涉及三维椭圆方程的求解,计算代价高昂,未进入业务系统中去(Cullen, 2003; Wlasak et al, 2006)。
选择了旋转风作为第一独立分析变量以后,将下一个分析变量分解为与已选独立分析变量相平衡的部分和非平衡部分,其中非平衡部分作为下一个独立分析变量,以此类推,直到完成所有独立分析变量的构造。若选择流函数ψ表征旋转风、势函数χ表征散度风、气压p表征质量、比湿q表征水汽,那么独立分析变量的构造可以表达如下:
$\begin{array}{l} \quad \quad \quad \psi = \psi \\ \quad \quad {\chi _{\rm{u}}} = \chi - \mathit{\boldsymbol{M}}\psi \\ \quad {p_{\rm{u}}} = p - \mathit{\boldsymbol{N}}\psi - \mathit{\boldsymbol{L}}{\chi _{\rm{u}}}\\ {q_{\rm{u}}} = q - \mathit{\boldsymbol{P}}\psi - \mathit{\boldsymbol{Q}}{\chi _{\rm{u}}} - \mathit{\boldsymbol{R}}{p_{\rm{u}}} \end{array}$ | (6) |
式中, M、N、L、P、Q和R为平衡算子,表达变量之间的动力平衡约束。与风压场平衡有关的M、N、L三个算子中,N算子最为重要,它表达了旋转风场和质量场之间的动力平衡约束,在赤道外自由大气中接近地转平衡关系,在所有变分同化系统中均予以考虑(见表 1第3列)。而M和L算子只在部分变分同化系统中有所考虑,其中M算子表达了旋转风场和散度风场之间的动力平衡约束,该平衡算子在边界层以及高空急流区影响较大,而在其他地区影响很小(Wu et al, 2002; 王瑞春等, 2012);L算子则表达了非平衡散度风和质量场之间的平衡约束,该约束一般较小,但在赤道地区以及平流层有一定影响(Chen et al, 2013)。与水汽有关的P、Q和R三个算子将水汽与动力学变量联系到了一起。但由于水汽本身的观测以及预报都存在很大误差,将其与动力学变量联系到一起还存在许多问题,只在部分区域系统中对这三个算子有所研究(Derber et al, 1999; Berre, 2000)。不过,目前欧洲中期天气预报中心(ECMWF)、英国气象局(Met Office)以及我国GRAPES等全球同化系统已开始研究引入湿度与温度之间的平衡约束,构造新的有关湿度的独立分析变量(Andersson et al, 2005; Ingleby et al, 2013; 龚建东等, 2016)。
对于式(6),还有两点需要注意。第一,该式在变分同化框架中一般无需显式出现,引入动力平衡约束的物理变换采用的是它的逆运算,也即式(5),此时可以展开为如下形式:
$\left[ {\begin{array}{*{20}{c}} \psi \\ \chi \\ p\\ q \end{array}} \right]{\rm{ = }}\left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{I}} & 0 & 0 & 0\\ \mathit{\boldsymbol{M}} & \mathit{\boldsymbol{I}} & 0 & 0\\ \mathit{\boldsymbol{N}} & \mathit{\boldsymbol{L}} & \mathit{\boldsymbol{I}} & 0\\ \mathit{\boldsymbol{P}} & \mathit{\boldsymbol{Q}} & \mathit{\boldsymbol{R}} & \mathit{\boldsymbol{I}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} \psi \\ {{\chi _{\rm{u}}}}\\ {{p_{\rm{u}}}}\\ {{q_{\rm{u}}}} \end{array}} \right]$ | (7) |
第二,式(6) 中我们采用了气压p作为质量变量,因而N算子主要体现的是风压场在水平方向的准地转约束,而静力平衡会独立地体现在由气压p导出温度T的环节中。而如果我们采用T作为质量变量,那么N就会涵盖准地转平衡和静力平衡两个约束:先由旋转风计算平衡的气压,再由平衡气压计算与之相平衡的温度。
考虑到准地转平衡和静力平衡是变分同化所需要考虑的主要动力平衡约束,下面两节就分别给出有关这两者求解方案的研究进展。其他平衡约束算子的求解方法与之类似,我们会在叙述过程中指出。
2 准地转平衡约束的求解准地转平衡约束的求解方案总结起来可以分为三类:第一类是采用动力平衡方程直接表达;第二类是采用统计方法,根据一定的样本统计出平衡算子;第三类是将平衡方程和统计方法结合起来使用,下面分别予以讨论。
2.1 动力平衡方程方案该方案采用动力气象学中推导得到的平衡(诊断)方程表达平衡约束,常用的平衡方程是线性平衡方程(linear balance equation, LBE)(吕美仲等, 2004):
${p_{\rm{b}}} = \mathit{\boldsymbol{N}}\left( \psi \right) = \nabla _{\rm{h}}^{ - 2}[\nabla {_{\rm{h}}} \cdot (\rho f\nabla {_{\rm{h}}}\psi )]$ | (8) |
式中,pb表示与旋转风相平衡的气压。LBE考虑了科氏参数随纬度的变化,是地转平衡方程在球面上的自然延拓(Daley, 1996)。有时为更加接近大气实际情形,可以采用更加复杂的非线性平衡方程(NLBE)替代LBE求解N算子(Fisher, 2003; Barker et al, 2004)。NLBE与梯度风关系相近,在处理具有强旋转的天气系统时比LBE要有一定优势(庄照荣等, 2006; 万齐林等, 2007)。
然而平衡方程在具体应用过程中会面临一些问题。Parrish等(1997)针对SSI的研究表明,式(8) 逐层求解平衡气压后,进一步由静力平衡差分导出温度分析时误差较大。其原因在于:一方面,逐层二维求解未考虑大气的垂直连续性,造成温度分析的垂直廓线出现锯齿状分布(Lorenc et al, 2003);另一方面,这也与SSI采用的Lorenz垂直离散方案有关,这将在第3节作进一步分析。此外,线性和非线性平衡方程均是基于中高纬大气运动特征简化得到,在热带地区并不适用。由于变分同化是全局求解的,平衡方程在热带地区会造成虚假平衡,严重影响分析效果(Daley, 1991; 王瑞春等, 2015a;2015b)。
由于难以给出简单的诊断方程,其他平衡算子的求解很少采用动力平衡方程方案。不过,Fisher (2003)在ECMWF变分同化系统中采用简化的准地转ω方程表达了旋转风和散度风之间的平衡约束(M算子)。而在其他系统中,M算子一般采用下文统计方案构造。
2.2 统计方案统计方案根据变量背景误差的估计样本,采用适当的线性回归方案直接统计求解不同变量之间的平衡算子(Parrish et al, 1997; Derber et al, 1999)。由于大气真值无法得知,背景误差的样本无法直接得到,因而需要采用一些替代样本。目前,替代样本获取的方法主要有两类:NMC方法以及集合方法(Fisher, 2003)。
统计方法既可以逐层构建统计关系,也可以在变量垂直廓线间建立统计关系,考虑不同层次上变量间的相关(Wu et al, 2002)。一般为考虑大气的三维特征,可以构建如下的多层回归方程(Wang et al, 2015):
${p_{{k_0}}} = \sum\limits_{k = 1}^K {N_{{k_0}}^k} {\psi _k} + {({p_{\rm{u}}})_{{k_0}}}$ | (9) |
式中,K是模式层总层数。给定p和ψ的误差样本,采用线性回归统计得到N算子。在地转平衡成立较好的地区,统计得到的平衡算子与LBE的效果类似;而在地转平衡不适用的地区,统计得到的平衡约束远小于LBE,风压场的分析更为独立,可以避免直接使用LBE作为平衡约束带来的虚假平衡问题。此外,式(9) 采用整层的旋转风场去求解每一层的平衡质量场,这样的处理考虑了大气在垂直方向的连续性,可以提高温度分析效果(Parrish et al, 1997; Wang et al, 2015)。
由于统计方法简单易行,除用于求解N算子以外,物理变换中其他平衡算子的研究和求解基本均采用了该方案(Derber et al, 1999; Berre, 2000)。但统计方法也有其自身的不足:一方面,它只能考虑变量之间的气候态的线性关系,难于将天气形势的具体变化考虑进去,也难于向非线性方向扩展;另一方面,该方法的实际效果依赖于所选用的预报误差估计样本,统计过程中需小心控制噪音的影响。
2.3 动力-统计相结合方案为了尽可能将动力平衡方程方案与统计方案的优点结合起来,并尽量减少各自的不足,动力-统计相结合方案被引入进来。其基本思路是:首先在水平方向根据平衡方程求解出与风场相平衡的质量场,然后再用统计得到的平衡算子对其垂直廓线做进一步处理,以修正平衡方程不适用地区的虚假平衡, 并改善变量的垂直分布(Lorenc et al, 2000; 王瑞春等, 2015a)。具体的,该方案首先根据式(8) 由旋转风场逐层求解与之相平衡的pb,然后再构造如下的回归方程:
${p_{b1,{k_0}}} = \sum\limits_{k = 1}^K {R_{{k_0}}^k} {p_{b,k}} + {p_{u,{k_0}}}$ | (10) |
在平衡方程基础上引入式(10) 后,大气在垂直方向的连续性被考虑进来,可以帮助提高温度分析效果。此外,在LBE不适用的地区,引入式(10) 可以减少这些地区风压场的耦合程度,避免虚假平衡。与单纯统计方法相比,动力与统计相结合的混合方案可以引入变量间的非线性平衡关系。例如,可以采用NLBE计算式(10) 中的pb,使得物理变换模型更为接近大气真实状况(Fisher, 2003)。
综上,动力平衡方程方案通过线性(或非线性)平衡方程给定准地转平衡约束,在中高纬具有较高的稳健性,但在热带地区会造成虚假平衡,并且逐层二维求解过程未考虑大气的三维连续性。统计以及动力-统计相结合的方案给定的约束,在中高纬与线性平衡方程相近(见图 1a~1c的对比),而在热带地区要远弱于线性平衡方程(见图 1d~1f的对比),风压场分析在该地区解耦,可以避免虚假平衡对同化分析的影响。此外,这两个方案考虑了不同层次上变量间的平衡约束,可以帮助改进温度同化分析效果,这在下一节将进一步分析。不过,统计的具体效果在很大程度上依赖于样本以及统计方案的选取,并且统计系数需要跟随模式和同化系统的变动及时更新,维护成本相对较高。
在1.2节中已经提到,物理变换中静力平衡约束的引入可以分为两种类型:一种是如式(7) 给出的那样,采用气压p作为描述质量的分析变量,温度T的分析由p根据静力平衡关系差分导出;另一种方式是选择T作为质量场分析变量(需要结合地面气压ps),p的分析由(T, ps)根据静力平衡关系积分得到。在这两种方案中,非静力平衡的信息都被同化系统滤除了。这一方面是因为大气运动以静力平衡为主,另一方面也是考虑到数值模式和观测系统均还不能准确地描述非静力平衡运动(Lorenc et al, 2000)。
前文提到,逐层采用动力平衡方程或统计方案求解准地转平衡约束会影响T分析效果。对于第一种引入静力平衡约束的方案而言,这易于理解,因为T的分析需要由p垂直差分得到,逐层求解易导致相邻层次上p的分析变动过大(Wang et al, 2015)。而对于采用T作为质量分析变量的第二种方案,情况同样如此。这是因为在引入质量场和风场间平衡约束时,还是先计算与风场相平衡的pb,再差分导出Tb(Parrish et al, 1997)。因而采用多层统计回归或者动力-统计相结合的方案构造准地转平衡约束,考虑大气在三维空间的连续性,对于两种静力平衡约束引入方式同样重要。
3.2 垂直离散方案的影响引入静力平衡的差分计算除了受准地转平衡约束求解方式影响以外,还会受模式和同化系统所采用的垂直离散方案的影响,出现欠定问题,严重影响分析质量。该问题在目前广泛使用的Lorenz以及Charney-Phillips离散方案中均有所表现,但形式有所不同。下面就分别以NCEP的SSI(Lorenz离散,见图 2a)以及Met Office的变分框架(Charney-Phillips离散,见图 2b)为例分别予以说明,并给出其他系统的情况。
SSI垂直坐标采用σ坐标,垂直离散采用Lorenz跳层方案,并采用涡度ζ表征旋转风,(T, ps)表征质量(Parrish et al, 1992)。在σ坐标中,线性平衡方程LBE可以表示为(Temperton et al, 1981):
${\Phi _{\rm{b}}} = \nabla {^{ - 2}}\{ \nabla \cdot f\nabla ({\nabla ^{ - 2}}\zeta )\} $ | (11) |
式中,Φ由温度T和ps结合静力平衡组合得到:(Φ1,Φ2,…,ΦN)NT=G(ps,T1,T2,…,TN)N+1T,G是线性变换算子。注意到,这里G算子并不是一个方阵,无法求逆,因而在物理变换中由Φb根据静力平衡导出(ps,T1,T2,…,TN)N+1T时就出现了欠定问题。为克服该问题,SSI在计算温度廓线时增加了一个约束条件:温度在垂直方向二阶导数的累加和最小。该约束的引入只能部分解决欠定问题,温度分析误差仍然较大。为此,Parrish等(1997)在针对SSI进一步优化中,根据误差样本构建(T, ps)和Φ间的统计约束关系,采用统计系数表达G-1,取得了较好的结果。而在SSI的升级版本GSI中,Wu等(2002)直接构建了(T, ps)和ζ间的统计约束关系,取消了LBE的显式计算,其形式类似于式(9)。这样的方案也被WRF系统所采用(Huang et al, 2009)。
上述欠定问题在其他采用Lorenz离散的系统中也有所出现,例如ECMWF的变分系统。ECMWF为解决该问题,起初也是引入约束条件,从G的广义逆矩阵中寻找一个符合条件的特解(Courtier et al, 1998)。之后,也沿用SSI的经验,采用统计方案求解G-1(Derber et al, 1999)。
Arakawa等(1988;1996)从动力框架设计角度指出,欠定问题的出现与Lorenz离散方案中变量分布密切相关。根据图 2a,Lorenz离散方案中风场和温度场变量被设定在同一层次上。而根据热成风关系,旋转风场(ζ1,ζ2,…,ζN)NT仅能确定
Charney-Phillips跳层离散方案中,风场变量放置于两层温度变量之间的层次上(见图 2b),Met Office变分同化系统就采用该设置(Lorenc et al, 2000)。该系统中位温有N+1层,比气压和风多一层,因而在由平衡气压导出平衡温度时也存有欠定问题。然而,由于温度和风场不在同一层次上,当风场廓线确定时,根据热成风关系可以唯一确定除上下边界以外的所有层次上的温度值,这与Lorenz网格下的欠定问题并不一样。进一步的,Met Office预报模式在第0层和第1层之间做了等温大气假设,同化框架不分析第0层的位温,并且设定顶层的位温是辐射边界条件,也不参与同化分析。通过这样的边界设置,引入静力平衡求解温度的欠定问题得以解决。
GRAPES变分同化系统也采用的是Charney-Phillips离散方案,因此欠定问题带来的温度分析的不确定性主要集中在上下边界处。之前的系统中是采用外插的形式确定上下边界的值,最新的GRAPES全球变分系统采用统计约束补充边界上的欠定条件。
综上,由于风场和温度之间满足三维热成风关系,变分框架在引入静力平衡时需要考虑两方面的影响:一方面,在求解风压场准地转平衡时需要考虑不同层次上的变量相关,保证垂直差分计算的温度在垂直方向的连续性;另一方面,要考虑由于垂直离散方案造成的求解欠定问题,Lorenz离散方案中整层温度的求解均是欠定的,Charney-Phillips离散方案中上下边界处的温度求解也是欠定的。在业务应用中,常通过整层统计求解(或动力-统计相结合)减缓这两方面的影响:一方面,整层统计考虑了不同层次上变量相关,保证了分析结果在垂直方向的连续性;另一方面,统计方案将气候态约束作为一类定解条件引入到求解中去,避免欠定问题的出现。
4 未来研究展望本文前面总结的动力平衡约束的求解方案究其本质都是气候态的平衡约束(Talagrand, 1997)。这些气候态约束抓取了中高纬大尺度大气运动中最主要的平衡约束信息——准地转平衡与静力平衡。这些平衡约束的构建在全球以及较粗网格的区域同化系统中取得了很大的成功,对NWP水平的稳步提升起到了重要作用。但这些气候态的平衡约束与真实大气变量间的复杂非线性关系相比,仍过于简单。随着NWP的不断发展,利用背景误差协方差准确构建变量间的动力平衡约束面临的挑战和机遇可以总结为以下几个方面。
4.1 赤道等特殊区域的平衡约束构建在赤道地区,准地转理论不再适用(Gill, 1982)。并且由于该地区风压场关系十分复杂,难于简化得到合适的平衡方程充当动力平衡约束。采用第2节介绍的统计方案或者动力-统计相结合方案给出的该区域的风压场平衡约束远弱于准地转约束,可以避免直接使用LBE带来的虚假平衡问题。不过,由于约束很弱,质量场中包含的风场信息难于提取出来。目前,有研究基于赤道波动理论,利用不同波动的正规模态构造B矩阵,以考虑不同赤道波动对于动力平衡约束的影响。该方法最早由Phillips (1986)提出,并在一些研究模型中得以应用(Daley, 1993; Žagar et al, 2004; Körnich et al, 2008)。但是由于在不同波动所占比例的确定,以及如何与中高纬衔接等方面仍有待进一步细化,该方案仍处于研究阶段。此外,随着模式层顶的不断提高,平流层甚至中间层等更高层次上的资料同化的重要性日益显现,而这些地区平衡约束的准确构建尚未实现(Polavarapu et al, 2005b)。
另一方面,与中高纬对流层相比,赤道等特殊区域的同化分析对于动力平衡约束的需求更为迫切。这些地区资料相对稀疏,并且主要集中于质量场观测,风场资料极度匮乏(Bauer et al, 2015)。因此,在赤道等特殊地区,通过合理的平衡约束帮助从质量观测中提取有效的风场信息,对于同化分析效果的提升十分关键。
4.2 高分辨率系统中平衡约束构建在高分辨率(特别是对流尺度的)同化系统中,前文总结的气候态平衡约束均已不再适用,而新的平衡约束尚未建立。Bannister等(2011)研究指出,当NWP模式水平分辨率<35 km时,静力平衡的重要性下降;而当水平分辨率<75 km时,地转平衡的重要性就会出现下降。Vetra-Carvalho等(2012)研究指出,当模式分辨率<20 km时,静力平衡在对流区不适用,但在非对流区仍然适用(1.5 km时仍是如此)。虽然,气候态平衡约束准确的适用范围仍有待进一步研究,但显然已不适合水平分辨率达到千米级的高分辨率同化系统。而在高分辨率变分同化系统中,继续采用气候态平衡约束会产生虚假平衡,阻碍资料的有效同化与应用(陆续等, 2015)。
为避免使用气候态平衡约束造成的虚假平衡,现有许多高分辨率同化框架直接不考虑变量间的误差协方差。而为减少由此带来的非平衡的影响,这些系统在目标函数中增加有关平衡的弱约束项(Liang et al, 2007; Gao et al, 2013)。弱约束项的引入虽然可以帮助改善分析场的平衡特征,但不足以补偿变量间误差协方差缺失造成的影响,并且还会影响极小化收敛的速度和精度。为此,如何以更为恰当的形式引入对流尺度下的动力平衡约束是高分辨率同化系统的重要挑战(Sun et al, 2014)。
4.3 混合变分同化带来的机遇与挑战最后,近十多年逐步发展起来的集合-变分混合同化系统为利用背景误差协方差构建动力平衡约束的发展提供了新的机遇。混合变分同化系统在原有气候态背景误差协方差的基础上通过集合预报样本引入具备流依赖特征的背景误差信息(Lorenc, 2003; 马旭林等, 2014)。在混合变分同化系统中,背景误差协方差矩阵可以表示为如下形式(Wang et al, 2007):
$\mathit{\boldsymbol{B}} = \alpha {\mathit{\boldsymbol{B}}_{\rm{C}}} + \left( {1 - \alpha } \right){\mathit{\boldsymbol{B}}_{\rm{E}}}$ | (12) |
式中,BC表示气候态的背景误差协方差矩阵,BE表示集合估计的背景误差协方差矩阵,α是标量参数,用于调节BE在整体B中所占的比重。与不随时间变化的BC不同,BE是由实时集合预报样本生成,考虑到了天气形势对背景误差的影响,也即BE具有流依赖特征。根据第1节有关B矩阵构造和求解的介绍,BE中流依赖特征可以来源于两个部分:一方面,由于集合样本中包含了完整的(至少是非线性模式给定的)变量间的平衡约束,合理应用可以带来更为全面的动力平衡约束,并可以随具体天气形势变化(秦琰琰等, 2012);另一方面独立分析变量的空间相关尺度以及背景误差方差也可以随流型演变(Clayton et al, 2013; 张卫民等, 2016)。图 4给出了这样的示例,当单点观测位于槽后脊前时,同化分析增量会沿着流线延展,并且风压场的配合也出现了一定的非地转特征。
Bauer等(2015)在讨论资料同化未来发展方向时指出,混合变分同化或者纯粹的集合资料同化在未来一段时间将成为主流资料同化框架。在新的框架中,基于集合样本信息给定的流依赖约束,不仅可以帮助改善全球和区域同化系统原有气候态约束的不足(Clayton et al, 2013),也为高分辨率同化系统弥补对流尺度平衡约束缺失的问题提供了可行解决方案(Sun et al, 2014)。不过,这一过程会受集合样本数不足、局地化方案破坏原有平衡等多方面因素的影响。研究、解决这些问题将会为动力平衡约束的不断改进提供较大的上升空间。
5 结论通过适当的动力平衡约束条件引导不同变量协同变化,保证分析场平衡,对于提高观测使用效率以及保证整个NWP系统的同化、预报效果具有十分重要的意义。这其中,基于背景误差协方差构造动力平衡约束是变分框架构造中的基础性通用方案。本文总结了这方面的研究进展,分析讨论了其中的关键问题,并展望了未来发展趋势。
变分同化框架在变分预处理的物理变换中构造不同变量的背景误差协方差,引入所考虑的动力平衡约束。在全球以及较粗网格的有限区域系统中主要考虑准地转平衡和静力平衡,因而常选用旋转风表达大气运动的平衡部分。在构建旋转风和质量场的地转平衡时,可以采用三类求解方案:平衡方程方案、统计方案以及动力-统计相结合的方案。动力平衡方程方案通过线性(或非线性)平衡方程给定准地转平衡约束,在中高纬具有较高的稳健性,但在热带地区会造成虚假平衡。统计以及动力-统计相结合的方案给定的约束,在中高纬与线性平衡方程相近,而在热带地区要远弱于线性平衡方程,帮助减少虚假平衡的影响。
在采用气压差分求解温度引入静力平衡时,计算精度受两方面的影响。一方面,在求解地转平衡约束时如采用整层统计或者动力-统计相结合的方案,考虑不同层次上变量相关,能更好保证气压在垂直方向的连续性,进而提高温度分析效果。另一方面,会受系统所采用的垂直离散方案的影响,Lorenz离散方案中整层温度的差分计算均是欠定的,Charney-Phillips离散方案中上下边界处的温度求解也是欠定的。为避免欠定问题的影响,可以通过统计得到的气候态约束补足定解条件,例如采用统计方案直接构造温度和旋转风间的平衡约束,也可以引入特定的边条件假设。
而在未来发展中,该领域研究面临的挑战和机遇仍然很多,本文简要讨论了其中几个较为突出的方面:赤道、平流层等特殊区域平衡约束的构造;高分辨率同化框架中准地转和静力平衡不再适用时,新的平衡约束如何构建;混合变分框架中如何更好发挥集合样本信息的作用,构造具备流依赖的平衡约束等。
龚建东, 王瑞春, 郝民, 2016. 温湿统计平衡约束关系对GRAPES全球湿度分析的作用[J]. 气象学报, 74(3): 380-396. DOI:10.11676/qxxb2016.035 |
吕美仲, 侯志明, 周毅, 2004. 动力气象学[M]. 北京: 气象出版社, 419.
|
陆续, 马旭林, 王旭光, 2015. 三维变分同化机载雷达资料对飓风预报的影响研究——2012年Isaac试验[J]. 大气科学, 39(6): 1111-1122. |
马旭林, 陆续, 于月明, 等, 2014. 数值天气预报中集合-变分混合资料同化及其研究进展[J]. 热带气象学报, 30(6): 1188-1195. |
秦琰琰, 龚建东, 李泽椿, 2012. 集合卡尔曼滤波同化多普勒雷达资料的观测系统模拟试验[J]. 气象, 38(5): 513-525. DOI:10.7519/j.issn.1000-0526.2012.05.001 |
万齐林, 薛纪善, 2007. 曲率修正线性平衡方程及其在变分同化风压约束中的应用[J]. 热带气象学报, 23(5): 417-423. |
王瑞春, 龚建东, 张林, 2012. GRAPES变分同化系统中动力平衡约束的统计求解[J]. 应用气象学报, 23(2): 129-138. DOI:10.11898/1001-7313.20120201 |
王瑞春, 龚建东, 张林, 等, 2015a. 热带风压场平衡特征及其对GRAPES同化预报的影响研究Ⅱ[J]. 动力与统计混合平衡约束方案的应用.大气科学, 39(6): 1225-1236. |
王瑞春, 龚建东, 张林, 等, 2015b. 热带风压场平衡特征及其对GRAPES同化预报的影响研究Ⅰ[J]. 平衡特征分析.大气科学, 39(5): 953-966. |
薛纪善, 庄世宇, 朱国富, 等, 2008. GRAPES新一代全球/区域变分同化系统研究[J]. 科学通报, 53(20): 2408-2417. DOI:10.3321/j.issn:0023-074X.2008.20.003 |
张华, 薛纪善, 庄世宇, 等, 2004. GRAPES三维变分同化系统的理想试验[J]. 气象学报, 62(1): 31-41. DOI:10.11676/qxxb2004.004 |
张卫民, 刘柏年, 曹小群, 等, 2016. 流依赖球面小波背景误差协方差模型设计和初步试验[J]. 气象学报, 74(3): 410-420. DOI:10.11676/qxxb2016.024 |
朱国富, 2015. 理解大气资料同化的根本概念[J]. 气象, 41(4): 456-463. DOI:10.7519/j.issn.1000-0526.2015.04.008 |
朱国富, 2015. 数值天气预报中分析同化基本方法的历史发展脉络和评述[J]. 气象, 41(8): 986-996. DOI:10.7519/j.issn.1000-0526.2015.08.008 |
朱国富, 2015. 理解大气资料同化的内在逻辑和若干共性特征[J]. 气象, 41(8): 997-1006. DOI:10.7519/j.issn.1000-0526.2015.08.009 |
朱宗申, 胡铭, 2002. 一种区域格点三维变分分析方案——基本框架和初步试验[J]. 大气科学, 26(5): 684-694. |
庄世宇, 薛纪善, 朱国富, 等, 2005. GRAPES全球三维变分同化系统——基本设计方案与理想试验[J]. 大气科学, 29(6): 872-884. |
庄照荣, 薛纪善, 朱宗申, 等, 2006. 非线性平衡方案在三维变分同化系统中的应用[J]. 气象学报, 64(2): 137-148. DOI:10.11676/qxxb2006.014 |
邹晓蕾, 2009. 资料同化理论和应用(上册)[M]. 北京: 气象出版社, 120.
|
Andersson E, Bauer P, Beljaars A, et al, 2005. Assimilation and modeling of the atmospheric hydrological cycle in the ECMWF forecasting system[J]. Bull Amer Meteor Soc, 86(3): 387-402. DOI:10.1175/BAMS-86-3-387 |
Arakawa A, Moorthi S, 1988. Baroclinic instability in vertically discrete systems[J]. J Atmos Sci, 45(11): 1688-1708. DOI:10.1175/1520-0469(1988)045<1688:BIIVDS>2.0.CO;2 |
Arakawa A, Konor C S, 1996. Vertical differencing of the primitive equations based on the Charney-Phillips grid in hybrid σ-p vertical coordinates[J]. Mon Wea Rev, 124(3): 511-528. DOI:10.1175/1520-0493(1996)124<0511:VDOTPE>2.0.CO;2 |
Bannister R N, 2008a. A review of forecast error covariance statistics in atmospheric variational data assimilation[J]. Ⅰ:Characteristics and measurements of forecast error covariances.Quart J Roy Meteor Soc, 134(637): 1951-1970. |
Bannister R N, 2008b. A review of forecast error covariance statistics in atmospheric variational data assimilation[J]. Ⅱ:Modelling the forecast error covariance statistics.Quart J Roy Meteor Soc, 134(637): 1971-1996. |
Bannister R N, Migliorini S, Dixon M A G, 2011. Ensemble prediction for nowcasting with a convection-permitting model -Ⅱ:Forecast error statistics[J]. Tellus A, 63(3): 497-512. DOI:10.1111/j.1600-0870.2010.00500.x |
Barker D M, Huang W, Guo Y R, et al, 2004. A three-dimensional variational data assimilation system for MM5:Implementation and initial results[J]. Mon Wea Rev, 132(4): 897-914. DOI:10.1175/1520-0493(2004)132<0897:ATVDAS>2.0.CO;2 |
Bauer P, Thorpe A, Brunet G, 2015. The quiet revolution of numerical weather prediction[J]. Nature, 525(7567): 47-55. DOI:10.1038/nature14956 |
Berre L, 2000. Estimation of synoptic and mesoscale forecast error covariances in a limited-area model[J]. Mon Wea Rev, 128(3): 644-667. DOI:10.1175/1520-0493(2000)128<0644:EOSAMF>2.0.CO;2 |
Bloom S C, Takacs L L, Da Silva A M, et al, 1996. Data assimilation using incremental analysis updates[J]. Mon Wea Rev, 124(6): 1256-1271. DOI:10.1175/1520-0493(1996)124<1256:DAUIAU>2.0.CO;2 |
Chen Y, Rizvi S R H, Huang X Y, et al, 2013. Balance characteristics of multivariate background error covariances and their impact on analyses and forecasts in tropical and Arctic regions[J]. Meteorol Atmos Phys, 121(1-2): 79-98. DOI:10.1007/s00703-013-0251-y |
Clayton A M, Lorenc A C, Barker D M, 2013. Operational implementation of a hybrid ensemble/4D-Var global data assimilation system at the Met Office[J]. Quart J Roy Meteor Soc, 139(675): 1445-1461. DOI:10.1002/qj.v139.675 |
Courtier P, Thépaut J N, Hollingsworth A, 1994. A strategy for operational implementation of 4D-Var, using an incremental approach[J]. Quart J Roy Meteor Soc, 120(519): 1367-1387. DOI:10.1002/(ISSN)1477-870X |
Courtier P, Andersson E, Heckley W, et al, 1998. The ECMWF implementation of three-dimensional variational assimilation (3D-Var)[J]. Ⅰ:Formulation.Quart J Roy Meteor Soc, 124(550): 1783-1807. |
Cullen M J P, 2003. Four-dimensional variational data assimilation:A new formulation of the background-error covariance matrix based on a potential-vorticity representation[J]. Quart J Roy Meteor Soc, 129(593): 2777-2796. DOI:10.1256/qj.02.10 |
Daley R, 1991. Atmospheric data analysis[M].
Cambridge: Cambridge University Press, 457.
|
Daley R, 1993. Atmospheric data assimilation on the equatorial beta plane[J]. Atmos Ocean, 31(4): 421-450. DOI:10.1080/07055900.1993.9649479 |
Daley R, 1996. Generation of global multivariate error covariances by singular-value decomposition of the linear balance equation[J]. Mon Wea Rev, 124(11): 2574-2587. DOI:10.1175/1520-0493(1996)124<2574:GOGMEC>2.0.CO;2 |
Derber J, Bouttier F, 1999. A reformulation of the background error covariance in the ECMWF global data assimilation system[J]. Tellus A, 51(2): 195-221. DOI:10.3402/tellusa.v51i2.12316 |
Fisher M.2003.Background error covariance modelling//ECMWF Seminar on Recent developments in data assimilation for atmosphere and ocean, 8-12 September 2003, 63.
|
Gao J, Smith T M, Stensrud D J, et al, 2013. A real-time weather-adaptive 3DVAR analysis system for severe weather detections and warnings[J]. Wea Forecasting, 28(3): 727-745. DOI:10.1175/WAF-D-12-00093.1 |
Gill A E, 1982. Atmosphere-ocean dynamics[M].
New York: Academic Press, 662.
|
Holton J R, 1992. An introduction to dynamic meteorology[M].
Vol.3rd, New York: Academic Press511.
|
Huang X Y, Xiao Q, Barker D M, et al, 2009. Four-dimensional variational data assimilation for WRF:formulation and preliminary results[J]. Mon Wea Rev, 137(1): 299-314. DOI:10.1175/2008MWR2577.1 |
Ingleby N B, Lorenc A C, Ngan K, et al, 2013. Improved variational analyses using a nonlinear humidity control variable[J]. Quart J Roy Meteor Soc, 139(676): 1875-1887. DOI:10.1002/qj.2073 |
Körnich H, Källén E, 2008. Combining the mid-latitudinal and equatorial mass/wind balance relationships in global data assimilation[J]. Tellus A, 60(2): 261-272. DOI:10.1111/j.1600-0870.2007.00286.x |
Kalnay E, 2003. Atmospheric modeling, data assimilation and predictability[J]. Cambridge:Cambridge University Press: 341. |
Kleist D T, Parrish D F, Derber J C, et al, 2009. Introduction of the GSI into the NCEP global data assimilation system[J]. Wea Forecasting, 24(6): 1691-1705. DOI:10.1175/2009WAF2222201.1 |
Krysta M, Rizvi S R H, Huang X Y.2009.A new formulation of WRFDA analysis control variables//10th annual WRF Users' Workshop, June 23-26 2009.
|
Liang X, Wang B, Chan J C L, et al, 2007. Tropical cyclone forecasting with model-constrained 3D-Var.Ⅰ:Description[J]. Quart J Roy Meteor Soc, 133(622): 147-153. DOI:10.1002/(ISSN)1477-870X |
Lorenc A C, 1986. Analysis methods for numerical weather prediction[J]. Quart J Roy Meteor Soc, 112(474): 1177-1194. DOI:10.1002/(ISSN)1477-870X |
Lorenc A C, 1997. Development of an operational variational assimilation scheme[J]. J Meteor Soc Japan, 75(1B): 339-346. DOI:10.2151/jmsj1965.75.1B_339 |
Lorenc A C, 2003. The potential of the ensemble Kalman filter for NWP-a comparison with 4D-Var[J]. Quart J Roy Meteor Soc, 129(595): 3183-3203. DOI:10.1256/qj.02.132 |
Lorenc A C, Ballard S P, Bell R S, et al, 2000. The Met.Office global three-dimensional variational data assimilation scheme[J]. Quart J Roy Meteor Soc, 126(570): 2991-3012. DOI:10.1002/(ISSN)1477-870X |
Lorenc A C, Roulstone I, White A.2003.On the choice of control fields in Var.Met Office Forecasting Research Technical Report No.419, 35.
|
Lynch P, Huang X Y.2010.Initialization//Data assimilation:Making Sense of Observations.Berlin Heidelberg:Springer, 241-260.
|
Parrish D F, Derber J C, 1992. The national meteorological center's spectral statistical interpolation analysis system[J]. Mon Wea Rev, 120(8): 1747-1763. DOI:10.1175/1520-0493(1992)120<1747:TNMCSS>2.0.CO;2 |
Parrish D F, Derber J C, Puser R J, et al, 1997. The NCEP global analysis system:recent improvements and future plans[J]. J Meteor Soc Japan, 75(1B): 359-365. DOI:10.2151/jmsj1965.75.1B_359 |
Phillips N A, 1986. The spatial statistics of random geostrophic modes and first-guess errors[J]. Tellus A, 38(4): 314-332. DOI:10.3402/tellusa.v38i4.11721 |
Polavarapu S, Ren S Z, Rochon Y, et al, 2005a. Data assimilation with the Canadian middle atmosphere model[J]. Atmos Ocean, 43(1): 77-100. DOI:10.3137/ao.430105 |
Polavarapu S, Shepherd T G, Rochon Y, et al, 2005b. Some challenges of middle atmosphere data assimilation[J]. Quart J Roy Meteor Soc, 131(613): 3513-3527. DOI:10.1256/qj.05.87 |
Rabier F, 2005. Overview of global data assimilation developments in numerical weather-prediction centres[J]. Quart J Roy Meteor Soc, 131(613): 3215-3233. DOI:10.1256/qj.05.129 |
Sun J, Xue M, Wilson J W, et al, 2014. Use of NWP for nowcasting convective precipitation:Recent progress and challenges[J]. Bull Amer Meteor Soc, 95(3): 409-426. DOI:10.1175/BAMS-D-11-00263.1 |
Talagrand O, 1997. Assimilation of observations, an introduction[J]. J Meteor Soc Japan, 75(1B): 191-209. DOI:10.2151/jmsj1965.75.1B_191 |
Talagrand O, 2010. Variational Assimilation//Data assimilation:Making Sense of Observations[M].
Berlin Heidelberg: Springer, 41-67.
|
Temperton C, Williamson D L, 1981. Normal mode initialization for a multilevel grid-point model[J]. Part Ⅰ:Linear aspects.Mon Wea Rev, 109(4): 729-743. |
Thepaut J-N, Andersson E.2010.The global observing system//Data assimilation:Making Sense of Observations.Berlin Heidelberg:Springer, 263-281.
|
Vetra-Carvalho S, Dixon M, Migliorini S, et al, 2012. Breakdown of hydrostatic balance at convective scales in the forecast errors in the Met Office Unified Model[J]. Quart J Roy Meteor Soc, 138(668): 1709-1720. DOI:10.1002/qj.v138.668 |
Wang R, Gong J, Zhang L, et al, 2015. Numerical experiments on multi-level statistical estimation of dynamic balance constraints in GRAPES-3DVAR[J]. J Trop Meteor, 21(4): 417-427. |
Wang X, Snyder C, Hamill T M, 2007. On the theoretical equivalence of differently proposed ensemble-3DVAR hybrid analysis schemes[J]. Mon Wea Rev, 135(1): 222-227. DOI:10.1175/MWR3282.1 |
Wlasak M, Nichols N K, Roulstone I, 2006. Use of potential vorticity for incremental data assimilation[J]. Quart J Roy Meteor Soc, 132(621C): 2867-2886. DOI:10.1256/qj.06.02 |
Wu W S, Purser R J, Parrish D F, 2002. Three-dimensional variational analysis with spatially inhomogeneous covariances[J]. Mon Wea Rev, 130(12): 2905-2916. DOI:10.1175/1520-0493(2002)130<2905:TDVAWS>2.0.CO;2 |
Žagar N, Gustafsson N, Köllén E, 2004. Variational data assimilation in the tropics:The impact of a background-error constraint[J]. Quart J Roy Meteor Soc, 130(596): 103-125. DOI:10.1256/qj.03.13 |