非线性集合四维变分同化方法NLS-4DVar之局地化改进❋

张洪芹， 田向军， 张承明

(1.山东农业大学，山东 泰安 271018; 2.中国科学院大气物理研究所国际气候与环境科学中心，北京 100029)



非线性集合四维变分同化方法NLS-4DVar之局地化改进❋

张洪芹1,2， 田向军2❋❋， 张承明1

(1.山东农业大学，山东 泰安 271018; 2.中国科学院大气物理研究所国际气候与环境科学中心，北京 100029)

四维变分同化可利用同化窗口内所有可能的观测信息优化大气、海洋模式的初始场，从而极大地提高大气、海洋模式模拟性能，而作为4DVar标准算法的伴随方法始终无法避免繁琐与复杂的预报模式伴随方程的编程、维护以及更新。为避免伴随模式的使用，集合四维变分方法，4DEnVar方法被逐渐开发，为4DVar的求解提供了一种便捷的途径。4DEnVar一般通过局地化过程消除样本不足所造成的虚假相关，而局地化方案的不同也必然会影响到其最终的同化效果。本文将一种集合样本扩展的局地化方案引入到基于Gaussian-Newton迭代算法的非线性集合四维变分同化方法NLS-4DVar中，从而避免了原算法中为进行局地化过程而额外需要的线性化假设，使得算法收敛更稳定。另外，通过将原Gaussian-Newton迭代序列进行变形、避免了矩阵的直接求逆，极大地提高了同化算法的计算效率。利用非线性动力模型Lorenz-96所开展的观测系统模拟试验表明：采用新的样本扩展型局地化方案的NLS-4DVar算法，其同化精度略优于NLS-4DVar原始算法，由于避免了矩阵的直接求逆，其计算效率反而有所提高，同化所需时间有所降低，对于大气与海洋数据同化领域的应用具有极大的潜力。

样本扩展型局地化方案； NLS-4Dvar； 共轭梯度法

引用格式：张洪芹， 田向军， 张承明, 等. 非线性集合四维变分同化方法NLS-4DVar之局地化改进[J]. 中国海洋大学学报(自然科学版), 2016, 46(10)： 10-15.

ZHANG Hong-Qin， TIAN Xiang-Jun, ZHANG Cheng-Ming. An improved localization scheme to the NLS-4DVar method[J]. Periodical of Ocean University of China, 2016, 46(10)： 10-15.

四维变分同化(4DVar)本质上属于最优控制问题，根据同化窗口内所有可能的观测信息优化大气、海洋模式的初始场，从而使得到的分析场与预报、观测之间距离最小，其优势明显: (1) 利用完整的模式方程作为动力约束; (2) 可以同时同化多个时刻的观测资料。但该方法也存在弊端: 4DVar给出的代价函数中, 控制变量(初始场)是以隐函数的形式出现的, 为了求得代价函数相对初始场的梯度, 需要给出预报模式的切线性模式及相应的伴随模式; 极小化代价函数也需要反复积分模式方程和伴随方程, 计算效率极低。自从4DVar的概念引入到数据同化领域以来，很多研究都致力于开发高效、稳健的算法用以求解4DVar[1-8]。直到现在，伴随法几乎被认为是求解4DVar问题的标准算法[1-3]；然而，伴随法始终无法避免繁琐与复杂的预报模式伴随方程的编程、维护以及更新。

为避免伴随模式的使用，所谓的4DEnVar方法被逐渐开发[9-11,14,17]，为4DVar的求解提供了一种便捷的途径。4DEnVar方法的理论基础是利用模拟样本所构成的线性空间近似逼近模式状态变量空间；基于这样的理论基础，4DVar的分析解应属于该线性样本空间且可表征为该空间内所有样本的线性组合；进一步又基于模式状态变量增量与观测向量增量线性相关的假设、可相对容易地求解4DVar的分析解[14]。另外，在4DEnVar的极小化过程也不需要使用伴随模式，大大降低了同化难度。基于4DEnVar方法的数据同化应用都显示了4DEnVar方法在各个领域的巨大潜力[12-13,16-18,21]。然而正如Tian and Feng[19]所指出，4DEnVar中所假设的模式状态变量增量与观测变量增量的线性相关(4DEnVar的理论基础之一)在预报模式或者观测算子高度非线性的情形下可能会出现问题。为了解决这个问题，Tian and Feng[19]提出来一种非线性最小二乘4DEnVar方法NLS-4DVar，采用Gaussian-Newton迭代策略应对数据同化问题中高度非线性问题。观测系统模拟试验[19]以及真实雷达资料的同化试验都表明NLS-4DVar确实可以很好地应对同化问题中高度非线性问题，其同化性能在一定程度要优于基于上述模式增量与观测增量线性假设的4DEnVar方法[22,25]。

另外，由于4DEnVar方法所构造的样本空间的维数要远远小于模式状态变量的总维数mx(对于复杂的非线性天气、气候模式而言尤为突出)，而通常需要采用所谓的“局地化”过程用以减小由于较少的样本个数所造成的采样误差与虚假相关[26]。为了解决这个问题，Houtekamer and Mitchell[27]通过Schür积B∘C对背景误差矩阵B进行修正，其中矩阵C的每个元素来自于认为给定的随某一个距离递减的相关函数。当中，Gaspari and Cohn[28]提出了正态分布的多项式估计递减函数。Houtekamer and Mitchell[29]和Hamill等[30]利用Gaspari and Cohn提出的多项式估计递减函数阐述了一种更一般的局地化，它是观测值和状态变量距离的函数，从而减少虚假相关。2010年，Wang等[17]将Schur积应用到矩阵中，以过滤掉观测站点和模式格点之间的虚假遥相关并提高了相关性。2011年，Tian等[14]将局地化方案引入到了PODEn4DVar同化方法中，并通过观测系统模拟试验进行了分析评估，结果显示引入局地化技术以后的同化结果要明显好于原来的同化结果。很自然地，不同的局地化方案的选择也必然会对相应的4DEnVar算法的同化效果产生影响。本文通过将一种样本扩展型局地化方案引入到NLS-4DVar算法当中，由此实现其同化性能的提升；另外，通过将NLS-4DVar的Gaussian-Newton迭代序列进行变形、避免了矩阵的直接求逆，从而同化效率比NLS-4DVar算法进一步优化。

1　方法

(1)

(2)

将式(2)带入到公式(1)同时经过一系列的数学变换(详见Tian and Feng[19])，式(1)可以转化为如下的非线性最小二乘的形式[22]:

(3)

这里

(4)

(5)

(6)

以及

(7)

(8)

基于公式(3、4、8), NLS-4DVar采用Gauss-Newton格式迭代求βi+1(i=0,1,2,…,Imax)[24]

(9)

(10)

因此公式(9)被重写为

(11)

(12)

因此，

(13)

(14)

(15)

(16)

矩阵ρ的定义如下：

ρi,j=C0(dh,i,j/dh,0)·C0(dv,i,j/dv,0) 。

(17)

(18)

其中dh,0与dv,0分别是水平和垂直的协方差局地化半径。

实际上，式(10)的导出是基于模式增量Px与观测增量Py的线性相关假设,意味着NLS-4DVar所具有的非线性优势也必然会在一定程度上被这种线性假设所影响。

鉴于此,本文将一种集合样本扩展型局地化方案引进NLS-4DVar[19,23]，同时对比了采用原始局地化方案产生对非线性特征的影响。在NLS-4DVar当中,背景误差协方差矩阵B可以利用模式状态变量增量集合Px进行如下的近似计算

(19)

类似地, 相关矩阵C(是相关长度[19]或局地化半径ρ的函数)也可以进行如下分解[23]

C=C′C'T。

(20)

(21)

其中“∘”代表2个矩阵的Schür积[19]，而局地化后的模式状态变量增量Px,ρ∈mx×(N×r)可表示为

Px,ρ=PxC′。

(22)

这里2个矩阵Px与ρ′之间的算子为以下操作

(23)

其中P*x,i(i=1,…,N)为一个r列矩阵，其每列即为Px的第i列。

理想地, 我们需要利用扩展后的初始样本(也就是xb+Px,ρ)反复积分预报模式Mt0→tk从而计算其相应的观测增量样本Py,ρ∈my×(N×r)(其中my为观测变量的总维数)。然而这几乎是不可能的，因为扩展后的样本数N×r急剧变大，其计算代价无法承受。取而代之，本文采用以下的方式近似求解相应的局地化观测变量增量矩阵Py，ρ∈my×(N×r)

(24)

其中

Py,ρ,k=Hk(xb,k+Px,ρ,k)-Hk(xb,k),

(25)

xb,k=Mt0→tk(xb),

(26)

Px,ρ,k=Px,kC′,

(27)

以及

Px,k=Mt0→tk(xb+Px)-xb,k。

(28)

经过这种样本扩展局地化过程之后,公式(9)可进一步地改写为

(29)

其中I为(N×r)×(N×r)单位矩阵。

很显然, 由于此时样本从N扩展到了N×r，直接求解以下的矩阵的逆的难度颇大

为了避免如下矩阵的直接求逆

本文进一步将(29) 转化为如下的形式

(30)

(31)

当然, C=C′C'T的直接分解的计算代价也很昂贵, 而这个过程又可以通过利用对一组Gaussian随机样本进行SVD分解后直接对C′近似来避免[19]事实上, 当C由广泛使用的Gaspari and Cohn[29]方程求得时，作者也已开发了一种高效算法(将另文介绍)来进行C=C′C'T分解。

2　数值试验

本文在此采用Lorenz-96模型[27]作为预报模式来对新发展的NLS-4DVar进行验证，主要与采用原始局地化方案的算法进行对比。Lorenz-96尽管简单，但其基本特征与真实的大气模式非常相似，是一个非常著名的非线性动力模式。作者将从同化精度与计算效率两个方面设计观测系统模拟试验对新发展NLS-4Dvar局地化算法进行验证。

2.1 试验设计

具有周期边界条件的Lorenz[27]模型由以下的方程所控制

(32)

2.2 试验结果

图1阐述了两种(采用新的和原始的)NLS-4DVar局地化算法在模式具有误差(F=10)的情形下的同化效果的对比图。总体而言，就均方根误差而言，两种局地化算法都表现相当不错，其同化结果的均方根误差都相当低。进一步可以看出采用新的局地化方案的迭代算法优势更为明显一点，要略优于采用原始方案的NLS-4DVar算法，尤其在初始的12个同化日阶段。基本上由于采用了新的局地化方案，因此现在的迭代算法不需要额外的线性假设，应该在一定程度上有助于缓解原始方案对于NLS-4DVar方法非线性特征的损害。

图1　新采用的(New)和原始的(Old) NLS-4DVar局地化算法在模式具有误差 (F=10)的情形下的日平均均方根误差的对比图Fig.1　Time series of the daily averaged root mean square (RMS) error for the new proposed (New) and the original (Old) localization schemes to the NLS-4DVar method with model error F=10, respectively

为了检验新采用局地化算法对于保留的相关矩阵特征根数目r的影响,本文进一步对比了新的局地化方案采用r=5, 10, 以及 15同化结果的日平均误差均方根误差。图2表明当r≥5的时候，新的局地化算法对于不同的特征根数目并不是特别明显，这说明所采用的局地化方案有望具有较高的计算效率。

图2　NLS-4DVar算法中新提出的局地化 方案对选取的特征根个数(r=5，10 和15) 的敏感性的日平均均方根误差序列Fig.2　Time series of the daily averaged root mean square (RMS) error for the new proposed localization scheme to the NLS-4DVar method with r=5, 10 and 15, respectively

表1　针对不同特征根数目r两种不同方式 (直接、间接求逆)的迭代方案同化的CPU时间表Table 1　The CPU times (s) for the two (direct and indirect) ways to the iterative scheme to conduct the whole assimilation as for the choice of the number r

另外，本文还对比了新的局地化方案(不直接求逆)与最初NLS-4DVar[19]的局地化方案的收敛速度与迭代方案，对于同样的收敛标准，两个方法都需要迭代三次实现收敛，同时其相应的CUP时间为2.308 649 s(新的局地化方案)和2.318 443 s(原始的局地化方案)。总体而言，新旧两种局地化方法的局地化方案计算效率相当，但新的局地化方案的同化效果更好一些。

3　结语

作为一种4DEnVar方法，非线性最小二乘集合四维变分同化方法NLS-4DVar也需要局地化过程以消除由于样本不足所造成的采样误差以及虚假相关。因而，其同化性能也必然随着其所采用的局地化方案的不同而变化。为了执行其原始的局地化方案，NLS-4DVar需要额外的关于模式变量增量与观测变量增量线性相关的假设，这必然会对NLS-4DVar算法所具有的非线性特征带来一定的冲击，因而也会对该算法的广泛应用带来某种不确定性。为了消除这种不确定性，本文将一种样本扩展型局地化方案引入到NLS-4DVar算法，进一步又对NLS-4DVar算法的Gaussian-Newton迭代序列进行变形，采用共轭梯度法以避免对矩阵的直接求逆，从而使得同化的计算效率得以大幅度提升。

本文采用几组观测系统模拟试验对新改正的NLS-4DVar进行检验：首先利用一组试验对新旧两种不同的局地化方案进行对比，发现就同化精度而言，两种方案都表现不错，但新的方法表现略优，说明不需要额外的线性假设对于NLS-4DVar算法的优化确实具有不错的效果；另外，对于所保留特征根数目r的敏感性试验表明，并不需要保留所有的特征根数据就可以得到不错的同化效果，说明该局地化方案的引入对于NLS-4DVar算法同化计算效率的提升具有不错的潜力，而采用共轭梯度法取代矩阵的直接求逆使得新算法的计算效率更为高效。

[1]Lewis J M, Derber J C. The use of the adjoint equation to solve a variational adjustment problem with advective constraints[J]. Tellus, 1985， 37A: 309-322.

[2]Le Dimet F X, Talagrand O. Variational algorithms for analysis and assimilation of meteorological observations: Theoretical aspects[J]. Tellus， 1986： 38A: 97-110.

[3]Courtier P and Talagrand O. Variational assimilation of meteorological observations with the adjoint vorticity equation II: Numerical results[J]. Q J R Meteorol Soc, 1987， 113: 1329-1347.

[4]Courtier P, Thepaut J N, Hollingsworth A. A strategy for operational implementation of 4DVar using an incremental approach[J]. Q J R Meteorol Soc, 1994， 120: 1367-1387.

[5]Bormann N, Thepaut J N. Impact of MODIS polar winds in ECMWF’s 4DVAR data assimilation system[J]. Mon Wea Rev, 2004, 132: 929-940.

[6]Bauer P, Lopez P, Benedetti A, et al. Implementation of 1D+4D-Var assimilation of precipitation affected microwave radiances at ECMWF II: 4D-Var[J]. Q J R Meteorol Soc, 2006, 132: 2307-2332.

[7]Rosmond T， Xu L. Development of NAVDAS-AR: nonlinear formulation and outer loop tests[J]. TellusA, 2006, 58(1): 45-58.

[8]Gauthier P, Tanguay M, Laroche S, et al. Extension of 3DVAR to 4DVAR: Implementation of 4DVAR at the Meteorological Service of Canada[J]. Mon Weather Rev, 2007, 135(6): 2339-2354.

[9]Qiu C, Shao A, Xu Q, et al. Fitting model fields to observations by using singular value decomposition: An ensemble-based 4DVar approach[J]. J Geophys Res, 2007, 112:11105. doi:10.1029/2006JD007994.

[10]Liu C, Xiao Q, Wang B. An ensemble-based four-dimensional variational data assimilation scheme. Part I: Technique formulation and preliminary test[J]. Mon Weather Rev, 2008, 136: 3363-3373.

[11]Tian X, Xie Z, Dai A. An ensemble-based explicit four-dimensional variational assimilation method[J]. J Geophys Res， 2008, 113: 21124. doi:10.1029/2008JD010358.

[12]Tian X, Xie Z, Dai A, et al. A dual-pass variational data assimilation framework for estimating soil moisture profiles from AMSR-E microwave brightness temperature[J]. J Geophys Res, 2009, 114: 16102. doi: 10.1029/2008JD011600.

[13]Tian X, Xie Z, Dai A, et al. A microwave land data assimilation system: Scheme and preliminary evaluation over China[J]. J Geophys Res, 2010, 115: 21113. doi: 10.1029/2010JD0 14370.

[14]Tian X, Xie Z, Sun Q. A POD-based ensemble four-dimensional variational assimilation method[J]. Tellus, 2011, 63A: 805-816.

[15]Tian X, Xie Z. Implementations of a square-root ensemble analysis and a hybrid localisation into the POD-based ensemble 4DVar[J]. Tellus A, 2012, 64, http://dx.doi.org/10.3402/tellusa.v64i0.18375.

[16]Tian X, Xie Z, Liu Y, et al. A joint data assimilation system(Tan-Tracker) to simultaneously estimate surface CO2fiuxesand 3-D atmospheric concentrations from observations[J]. Atmos Chem Phys, 2014, 14: 13281-13293.

[17]Wang B, Liu J, Wang S, et al. An economical approach to four-dimensional variational data assimilation[J]. Adv Atmos Sci, 2010, 27(4): 715-727, doi: 10.1007/s00376-009-9122-3.

[18]Pan X, Tian X, Li X, et al. Assimilating Doppler radar radial velocity and reflectivity observations in the weather research and forecasting model by a proper orthogonal decomposition based ensemble three-dimensional variational assimilation method[J]. J Geophys Res, 2012, 117, D17113, doi:10.1029/2012JD017684.

[19]Tian X, Feng X. A non-linear least squares enhanced POD-4DVaralgorithm for data assimilation[J]. Tellus A, 2015, 67, 25340, http://dx.doi.org/10.3402/tellusa.v67.25340.

[20]Tian X, Feng X, Zhang H. et al. An enhanced ensemble-based method for computing the CNOPs using an efficient localization implementation scheme and a two-step optimization strategy: formulation and preliminary tests[J]. Q J R Meteorol Soc, 142:1007-1016， DOI:10.1002v.2703.

[21]Zhang B, Tian X, Sun J, et all. PODEn4DVar-based radar data assimilation scheme: Formulation and preliminary results from real-data experiments with advanced research WRF (ARW)[J]. Tellus A. 2015a, 67, 26045, http://dx.doi.org/10.3402/tellusa.v67.26045.

[22]Dennis J E Jr, Schnabel R B. Numerical Methods for Unconstrained Optimization and Nonlinear Equations (Classics in Applied Mathematics)[M]. SIAM, Philadelphia, 1996, 378. ISBN:0-89871-364-1.

[23]Liu C, Xiao Q, and Wang B. An ensemble-based four dimensional variational data assimilation scheme. part II: Observing system simulation experiments with advanced research WRF (ARW)[J]. Mon Wea Rev, 2009, 137: 1687-1704.

[24]Gaspari G, Cohn S E. Construction of correlation functions in two and three dimensions[J]. Quart J R Meteorol Soc, 1999, 125: 723-757.

[25]Lorenz E. Predictability: A Problem Partly Solved[R]. Proc. Seminaron Predictability. Volume 1. United Kingdom: ECMWF, 1996: 119.

[26]Houtekamer P L, Mitchell H L. A sequential ensemble Kalman filter for atmospheric data assimilation[J]. Mon Wea Rev, 2001,129: 123-137.

[27]Houtekamer P L, Mitchell H L. Data assimilation using an ensemble Kalman filter technique[J]. Mon Wea Rev, 1998, 126: 796-811.

[28]Gaspari G, Cohn S E. Construction of correlation functions in two and three dimensions[J]. Quart J Roy Meteor Soc, 1999, 125: 723-757.

[29]Houtekamer P L, Mitchell H L. A sequential ensemble Kalman filter for atmospheric data assimilation[J]. Mon Wea Rev, 2001, 129: 123-137.

[30]Hamill T M, Whitaker J S, Snyder C. Distancedependent filtering of background error covariance estimates in an ensemble Kalman filter[J]. Mon Wea Rev, 2001, 129: 2776-2790.

责任编辑庞旻

An Improved Localization Scheme to the NLS-4DVar Method

ZHANG Hong-Qin1,2， TIAN Xiang-Jun2, ZHANG Cheng-Ming1

(1.College of Information Science and Engineering, Shandong Agriculture University， Taian 271018, China; 2.ICCES, Institute of Atmospheric Physics, Chinese Academy of Sciences, Beijing 100029, China)

The 4DEnVar methods often contain a localization implementation scheme and its final algorithm after the localization implementation will be certainly changed with the choice of its localization scheme. As a step in the improvement of the NLS-4DVar method, we implement an expanding sample localization scheme into NLS-4DVar with an aim at removing the additional linear assumption adopted in the original NLS-4DVar algorithm. And the new proposed iterative scheme can give a fair good performance and its computational costs can be further reduced through an indirect way to avoid the computation of a matrix inverse. Numerical experiments with a nonlinear model of the Lorenz-96 equation show that the proposed new iterative scheme behaves slightly better than the original one in terms of the assimilation precision presumably due to the additional linear assumption adopted in the original one and its computational costs can be further reduced.

expanding example localization scheme； NLS-4DVar； the conjugate gradient method

国家高技术研究发展计划项目(2013AA122002);国家自然科学基金项目(41575100,91437220);山东省省级水利科研与技术推广项目(SDSLKY201503)资助

2016-01-19；

2016-03-09

张洪芹(1988-)，女，硕士生。E-mail:hqzhang1112@163.com

❋❋通讯作者：E-mail:tianxj@mail.iap.ac.cn

P40

A

1672-5174(2016)10-010-06

10.16441/j.cnki.hdxb.20160010

Supported by the National High-Tech R&D Program(2013AA122002); National Natural Science Foundation of China(41575100, 91437220); Shandong Provincial Water Conservancy Scientific Research and Technology Program(SDSLKY201503)