非線性集合四維變分同化方法NLS-4DVar之局地化改進?

張洪芹， 田向軍， 張承明

(1.山東農(nóng)業(yè)大學(xué)，山東 泰安 271018; 2.中國科學(xué)院大氣物理研究所國際氣候與環(huán)境科學(xué)中心，北京 100029)

?

非線性集合四維變分同化方法NLS-4DVar之局地化改進?

張洪芹1,2， 田向軍2??， 張承明1

(1.山東農(nóng)業(yè)大學(xué)，山東 泰安 271018; 2.中國科學(xué)院大氣物理研究所國際氣候與環(huán)境科學(xué)中心，北京 100029)

四維變分同化可利用同化窗口內(nèi)所有可能的觀測信息優(yōu)化大氣、海洋模式的初始場，從而極大地提高大氣、海洋模式模擬性能，而作為4DVar標(biāo)準(zhǔn)算法的伴隨方法始終無法避免繁瑣與復(fù)雜的預(yù)報模式伴隨方程的編程、維護以及更新。為避免伴隨模式的使用，集合四維變分方法，4DEnVar方法被逐漸開發(fā)，為4DVar的求解提供了一種便捷的途徑。4DEnVar一般通過局地化過程消除樣本不足所造成的虛假相關(guān)，而局地化方案的不同也必然會影響到其最終的同化效果。本文將一種集合樣本擴展的局地化方案引入到基于Gaussian-Newton迭代算法的非線性集合四維變分同化方法NLS-4DVar中，從而避免了原算法中為進行局地化過程而額外需要的線性化假設(shè)，使得算法收斂更穩(wěn)定。另外，通過將原Gaussian-Newton迭代序列進行變形、避免了矩陣的直接求逆，極大地提高了同化算法的計算效率。利用非線性動力模型Lorenz-96所開展的觀測系統(tǒng)模擬試驗表明：采用新的樣本擴展型局地化方案的NLS-4DVar算法，其同化精度略優(yōu)于NLS-4DVar原始算法，由于避免了矩陣的直接求逆，其計算效率反而有所提高，同化所需時間有所降低，對于大氣與海洋數(shù)據(jù)同化領(lǐng)域的應(yīng)用具有極大的潛力。

樣本擴展型局地化方案； NLS-4Dvar； 共軛梯度法

引用格式：張洪芹， 田向軍， 張承明, 等. 非線性集合四維變分同化方法NLS-4DVar之局地化改進[J]. 中國海洋大學(xué)學(xué)報(自然科學(xué)版), 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)本質(zhì)上屬于最優(yōu)控制問題，根據(jù)同化窗口內(nèi)所有可能的觀測信息優(yōu)化大氣、海洋模式的初始場，從而使得到的分析場與預(yù)報、觀測之間距離最小，其優(yōu)勢明顯: (1) 利用完整的模式方程作為動力約束; (2) 可以同時同化多個時刻的觀測資料。但該方法也存在弊端: 4DVar給出的代價函數(shù)中, 控制變量(初始場)是以隱函數(shù)的形式出現(xiàn)的, 為了求得代價函數(shù)相對初始場的梯度, 需要給出預(yù)報模式的切線性模式及相應(yīng)的伴隨模式; 極小化代價函數(shù)也需要反復(fù)積分模式方程和伴隨方程, 計算效率極低。自從4DVar的概念引入到數(shù)據(jù)同化領(lǐng)域以來，很多研究都致力于開發(fā)高效、穩(wěn)健的算法用以求解4DVar[1-8]。直到現(xiàn)在，伴隨法幾乎被認(rèn)為是求解4DVar問題的標(biāo)準(zhǔn)算法[1-3]；然而，伴隨法始終無法避免繁瑣與復(fù)雜的預(yù)報模式伴隨方程的編程、維護以及更新。

為避免伴隨模式的使用，所謂的4DEnVar方法被逐漸開發(fā)[9-11,14,17]，為4DVar的求解提供了一種便捷的途徑。4DEnVar方法的理論基礎(chǔ)是利用模擬樣本所構(gòu)成的線性空間近似逼近模式狀態(tài)變量空間；基于這樣的理論基礎(chǔ)，4DVar的分析解應(yīng)屬于該線性樣本空間且可表征為該空間內(nèi)所有樣本的線性組合；進一步又基于模式狀態(tài)變量增量與觀測向量增量線性相關(guān)的假設(shè)、可相對容易地求解4DVar的分析解[14]。另外，在4DEnVar的極小化過程也不需要使用伴隨模式，大大降低了同化難度。基于4DEnVar方法的數(shù)據(jù)同化應(yīng)用都顯示了4DEnVar方法在各個領(lǐng)域的巨大潛力[12-13,16-18,21]。然而正如Tian and Feng[19]所指出，4DEnVar中所假設(shè)的模式狀態(tài)變量增量與觀測變量增量的線性相關(guān)(4DEnVar的理論基礎(chǔ)之一)在預(yù)報模式或者觀測算子高度非線性的情形下可能會出現(xiàn)問題。為了解決這個問題，Tian and Feng[19]提出來一種非線性最小二乘4DEnVar方法NLS-4DVar，采用Gaussian-Newton迭代策略應(yīng)對數(shù)據(jù)同化問題中高度非線性問題。觀測系統(tǒng)模擬試驗[19]以及真實雷達資料的同化試驗都表明NLS-4DVar確實可以很好地應(yīng)對同化問題中高度非線性問題，其同化性能在一定程度要優(yōu)于基于上述模式增量與觀測增量線性假設(shè)的4DEnVar方法[22,25]。

另外，由于4DEnVar方法所構(gòu)造的樣本空間的維數(shù)要遠遠小于模式狀態(tài)變量的總維數(shù)mx(對于復(fù)雜的非線性天氣、氣候模式而言尤為突出)，而通常需要采用所謂的“局地化”過程用以減小由于較少的樣本個數(shù)所造成的采樣誤差與虛假相關(guān)[26]。為了解決這個問題，Houtekamer and Mitchell[27]通過Schür積B°C對背景誤差矩陣B進行修正，其中矩陣C的每個元素來自于認(rèn)為給定的隨某一個距離遞減的相關(guān)函數(shù)。當(dāng)中，Gaspari and Cohn[28]提出了正態(tài)分布的多項式估計遞減函數(shù)。Houtekamer and Mitchell[29]和Hamill等[30]利用Gaspari and Cohn提出的多項式估計遞減函數(shù)闡述了一種更一般的局地化，它是觀測值和狀態(tài)變量距離的函數(shù)，從而減少虛假相關(guān)。2010年，Wang等[17]將Schur積應(yīng)用到矩陣中，以過濾掉觀測站點和模式格點之間的虛假遙相關(guān)并提高了相關(guān)性。2011年，Tian等[14]將局地化方案引入到了PODEn4DVar同化方法中，并通過觀測系統(tǒng)模擬試驗進行了分析評估，結(jié)果顯示引入局地化技術(shù)以后的同化結(jié)果要明顯好于原來的同化結(jié)果。很自然地，不同的局地化方案的選擇也必然會對相應(yīng)的4DEnVar算法的同化效果產(chǎn)生影響。本文通過將一種樣本擴展型局地化方案引入到NLS-4DVar算法當(dāng)中，由此實現(xiàn)其同化性能的提升；另外，通過將NLS-4DVar的Gaussian-Newton迭代序列進行變形、避免了矩陣的直接求逆，從而同化效率比NLS-4DVar算法進一步優(yōu)化。

1　方法

(1)

(2)

將式(2)帶入到公式(1)同時經(jīng)過一系列的數(shù)學(xué)變換(詳見Tian and Feng[19])，式(1)可以轉(zhuǎn)化為如下的非線性最小二乘的形式[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分別是水平和垂直的協(xié)方差局地化半徑。

實際上，式(10)的導(dǎo)出是基于模式增量Px與觀測增量Py的線性相關(guān)假設(shè),意味著NLS-4DVar所具有的非線性優(yōu)勢也必然會在一定程度上被這種線性假設(shè)所影響。

鑒于此,本文將一種集合樣本擴展型局地化方案引進NLS-4DVar[19,23]，同時對比了采用原始局地化方案產(chǎn)生對非線性特征的影響。在NLS-4DVar當(dāng)中,背景誤差協(xié)方差矩陣B可以利用模式狀態(tài)變量增量集合Px進行如下的近似計算

(19)

類似地, 相關(guān)矩陣C(是相關(guān)長度[19]或局地化半徑ρ的函數(shù))也可以進行如下分解[23]

C=C′C'T。

(20)

(21)

其中“°”代表2個矩陣的Schür積[19]，而局地化后的模式狀態(tài)變量增量Px,ρ∈mx×(N×r)可表示為

Px,ρ=PxC′。

(22)

這里2個矩陣Px與ρ′之間的算子為以下操作

(23)

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

理想地, 我們需要利用擴展后的初始樣本(也就是xb+Px,ρ)反復(fù)積分預(yù)報模式Mt0→tk從而計算其相應(yīng)的觀測增量樣本Py,ρ∈my×(N×r)(其中my為觀測變量的總維數(shù))。然而這幾乎是不可能的，因為擴展后的樣本數(shù)N×r急劇變大，其計算代價無法承受。取而代之，本文采用以下的方式近似求解相應(yīng)的局地化觀測變量增量矩陣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)

經(jīng)過這種樣本擴展局地化過程之后,公式(9)可進一步地改寫為

(29)

其中I為(N×r)×(N×r)單位矩陣。

很顯然, 由于此時樣本從N擴展到了N×r，直接求解以下的矩陣的逆的難度頗大

為了避免如下矩陣的直接求逆

本文進一步將(29) 轉(zhuǎn)化為如下的形式

(30)

(31)

當(dāng)然, C=C′C'T的直接分解的計算代價也很昂貴, 而這個過程又可以通過利用對一組Gaussian隨機樣本進行SVD分解后直接對C′近似來避免[19]事實上, 當(dāng)C由廣泛使用的Gaspari and Cohn[29]方程求得時，作者也已開發(fā)了一種高效算法(將另文介紹)來進行C=C′C'T分解。

2　數(shù)值試驗

本文在此采用Lorenz-96模型[27]作為預(yù)報模式來對新發(fā)展的NLS-4DVar進行驗證，主要與采用原始局地化方案的算法進行對比。Lorenz-96盡管簡單，但其基本特征與真實的大氣模式非常相似，是一個非常著名的非線性動力模式。作者將從同化精度與計算效率兩個方面設(shè)計觀測系統(tǒng)模擬試驗對新發(fā)展NLS-4Dvar局地化算法進行驗證。

2.1 試驗設(shè)計

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

(32)

2.2 試驗結(jié)果

圖1闡述了兩種(采用新的和原始的)NLS-4DVar局地化算法在模式具有誤差(F=10)的情形下的同化效果的對比圖。總體而言，就均方根誤差而言，兩種局地化算法都表現(xiàn)相當(dāng)不錯，其同化結(jié)果的均方根誤差都相當(dāng)?shù)?。進一步可以看出采用新的局地化方案的迭代算法優(yōu)勢更為明顯一點，要略優(yōu)于采用原始方案的NLS-4DVar算法，尤其在初始的12個同化日階段?；旧嫌捎诓捎昧诵碌木值鼗桨?，因此現(xiàn)在的迭代算法不需要額外的線性假設(shè)，應(yīng)該在一定程度上有助于緩解原始方案對于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

為了檢驗新采用局地化算法對于保留的相關(guān)矩陣特征根數(shù)目r的影響,本文進一步對比了新的局地化方案采用r=5, 10, 以及 15同化結(jié)果的日平均誤差均方根誤差。圖2表明當(dāng)r≥5的時候，新的局地化算法對于不同的特征根數(shù)目并不是特別明顯，這說明所采用的局地化方案有望具有較高的計算效率。

圖2　NLS-4DVar算法中新提出的局地化 方案對選取的特征根個數(shù)(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　針對不同特征根數(shù)目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]的局地化方案的收斂速度與迭代方案，對于同樣的收斂標(biāo)準(zhǔn)，兩個方法都需要迭代三次實現(xiàn)收斂，同時其相應(yīng)的CUP時間為2.308 649 s(新的局地化方案)和2.318 443 s(原始的局地化方案)?？傮w而言，新舊兩種局地化方法的局地化方案計算效率相當(dāng)，但新的局地化方案的同化效果更好一些。

3　結(jié)語

作為一種4DEnVar方法，非線性最小二乘集合四維變分同化方法NLS-4DVar也需要局地化過程以消除由于樣本不足所造成的采樣誤差以及虛假相關(guān)。因而，其同化性能也必然隨著其所采用的局地化方案的不同而變化。為了執(zhí)行其原始的局地化方案，NLS-4DVar需要額外的關(guān)于模式變量增量與觀測變量增量線性相關(guān)的假設(shè)，這必然會對NLS-4DVar算法所具有的非線性特征帶來一定的沖擊，因而也會對該算法的廣泛應(yīng)用帶來某種不確定性。為了消除這種不確定性，本文將一種樣本擴展型局地化方案引入到NLS-4DVar算法，進一步又對NLS-4DVar算法的Gaussian-Newton迭代序列進行變形，采用共軛梯度法以避免對矩陣的直接求逆，從而使得同化的計算效率得以大幅度提升。

本文采用幾組觀測系統(tǒng)模擬試驗對新改正的NLS-4DVar進行檢驗：首先利用一組試驗對新舊兩種不同的局地化方案進行對比，發(fā)現(xiàn)就同化精度而言，兩種方案都表現(xiàn)不錯，但新的方法表現(xiàn)略優(yōu)，說明不需要額外的線性假設(shè)對于NLS-4DVar算法的優(yōu)化確實具有不錯的效果；另外，對于所保留特征根數(shù)目r的敏感性試驗表明，并不需要保留所有的特征根數(shù)據(jù)就可以得到不錯的同化效果，說明該局地化方案的引入對于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.

責(zé)任編輯龐旻

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

國家高技術(shù)研究發(fā)展計劃項目(2013AA122002);國家自然科學(xué)基金項目(41575100,91437220);山東省省級水利科研與技術(shù)推廣項目(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)