解析四維集合變分參數(shù)優(yōu)化方法研究

賈彬鶴，李威*，梁康壯*

(1.天津大學(xué) 海洋科學(xué)與技術(shù)學(xué)院，天津 300072)

1 引言

進(jìn)行數(shù)值預(yù)報(bào)時(shí)，預(yù)報(bào)的準(zhǔn)確率會隨著預(yù)報(bào)時(shí)間的增長而迅速下降。造成這種現(xiàn)象的原因主要有3點(diǎn)：一是模式方程對物理過程的參數(shù)化方案有一定誤差；二是數(shù)值預(yù)報(bào)模式的初始場有誤差；三是數(shù)值預(yù)報(bào)的計(jì)算方法有一定的誤差[1]。

為提高預(yù)報(bào)準(zhǔn)確率，如何改進(jìn)模式物理參數(shù)是需要解決的問題。數(shù)據(jù)同化方法可以從觀測數(shù)據(jù)中提取有效信息，并修正數(shù)值模式誤差。四維變分?jǐn)?shù)據(jù)同化 （Four Dimensional Variational Data Assimilation，4DVar）方法能夠把各個不同時(shí)刻的觀測資料納入統(tǒng)一的分析預(yù)報(bào)系統(tǒng)中，通過相應(yīng)的切線性模式和伴隨模式優(yōu)化初值條件，并使各要素自然滿足動力熱力條件[2–3]。由于其較好的同化能力，眾多學(xué)者也采用4DVar方法進(jìn)行了模式參數(shù)優(yōu)化的研究[4–7]。例如，Inazu等[8]將進(jìn)化算法應(yīng)用于區(qū)域海潮模型，優(yōu)化了模型的邊界條件和物理參數(shù)。Wang等[9]提出通過遺傳算法和神經(jīng)網(wǎng)絡(luò)交替組合的方式來糾正系統(tǒng)參數(shù)誤差。王云峰等[10]提出了一種新的利用觀測資料來同時(shí)優(yōu)化模式初始場和物理參數(shù)的擴(kuò)展四維變分同化方法，并以Ekman邊界層模式和Lorenz模型為例進(jìn)行了數(shù)值試驗(yàn)。魏敏[1]采用4DVar方法，利用Lorenz模型和Hamilton方程，結(jié)合觀測數(shù)據(jù)對該模型的參數(shù)進(jìn)行了優(yōu)化。然而，4DVar方法需要采用伴隨模式，對于原數(shù)值模式需要編寫與之相對應(yīng)的特定的伴隨模式，并且伴隨復(fù)雜程度高、代碼編輯量大，使得這種同化方法可移植性很差。另一方面，4DVar采用的靜態(tài)背景誤差協(xié)方差矩陣不具有流依賴特性，限制了四維變分的同化能力。

序列同化方法如集合卡爾曼濾波（Ensemble Kalman Filter，EnKF）能在避免伴隨的同時(shí)產(chǎn)生流依賴的背景誤差協(xié)方差矩陣。為結(jié)合兩種方法的優(yōu)點(diǎn)，可以在變分框架下引入集合狀態(tài)變量，從而產(chǎn)生四維集合變 分 （Four Dimensional Ensemble Variational Data Assimilation，4DEnVar）。最初，Liu 等[11]將四維變分增量優(yōu)化方法與集合方法結(jié)合提出了該方法，但將其命名為集合四維變分?jǐn)?shù)據(jù)同化方法（Ensemble Four Dimensional Variational Data Assimilation， En4DVar） ， 并采用一維淺水模型對該方法進(jìn)行了驗(yàn)證。隨后，Liu等[12]設(shè)計(jì)了局地化方法，采用高等天氣研究與預(yù)報(bào)（ARW-WRF）模型，設(shè)計(jì)了觀測系統(tǒng)模擬試驗(yàn)，并檢驗(yàn)了它們在實(shí)際數(shù)據(jù)同化中的性能，證實(shí)了4DEnVar局地化技術(shù)可以有效地減輕采樣誤差對分析的影響。根據(jù)世界氣象組織（WMO）會議建議，Liu和Xiao[13將En4DVar更名為4DEnVar，以便將4DEnVar和使用伴隨模型的En4DVar區(qū)分開來。4DEnVar 方法在業(yè)務(wù)化系統(tǒng)中也有所應(yīng)用。Arbogast等[14]對四維軌跡的輸入和存儲等問題，提出了一種具有分布式輸入和存儲擾動的四維集合變分集成并行實(shí)現(xiàn)辦法。Yang和Mémin[15]探索了一個動力學(xué)公式，該公式允許在一個大規(guī)模的流動模型中同化高分辨率的數(shù)據(jù)，并將該原理結(jié)合4DEnVar技術(shù)來估計(jì)隨機(jī)淺水模型的流動初始條件和非均勻時(shí)變子網(wǎng)格參數(shù)。Song和Kang[16]采用混合四維集合變分方法（hybrid-4DEnVar），改進(jìn)了北半球500 hPa夏季月位勢高度預(yù)報(bào)。楊雨軒[17]利用WRF模型研究了POD-4DEnVar和傳統(tǒng)三維變分方 法 （Three Dimensional Variational Data Assimilation 3DVar）對華南暴雨的模擬效果，并研究了物理集合方案對4DEnVar的影響。

4DEnVar可以構(gòu)造多個能反映出背景誤差協(xié)方差分布特征的樣本形成集合，為變分同化系統(tǒng)提供流依賴背景誤差協(xié)方差估計(jì)。然而，需要指明的是Liu等[11]將集合成員擾動關(guān)于集合均值展開來構(gòu)造的協(xié)方差矩陣，反映了集合的統(tǒng)計(jì)性質(zhì)，而非誤差的傳播性質(zhì)。在線性模型下，這兩種協(xié)方差矩陣等價(jià)，但對于非線性模型則不然。此外，目前4DEnVar方法研究集中在初始場的優(yōu)化，尚未有進(jìn)行參數(shù)優(yōu)化的研究。為此，本文提出了專門進(jìn)行模式參數(shù)優(yōu)化的解析四維集合變分?jǐn)?shù)據(jù)同化（Analytical Four Dimensional Ensemble Variational Data Assimilation， A-4DEnVar） 方法，該方法在4DVar框架下，將集合樣本擾動關(guān)于迭代生成的模式參數(shù)處展開，有效刻畫了誤差在非線性模型下的傳播方式，避免了伴隨模式的使用。在此基礎(chǔ)上，進(jìn)一步求得代價(jià)函數(shù)梯度為0的解析解，能迅速得到代價(jià)函數(shù)極小值，完成模式參數(shù)優(yōu)化。本文采用 Lorenz-63 模型[18]通過 4DVar與 A-4DEnVar 方法進(jìn)行參數(shù)優(yōu)化孿生試驗(yàn)和敏感性試驗(yàn)，證實(shí)算法的有效性。

2 解析四維集合變分參數(shù)優(yōu)化方法

假設(shè)動力系統(tǒng)為

式中，Xi?1為第 (i?1)時(shí)刻的狀態(tài)變量；Xi為 第i時(shí)刻的狀態(tài)變量；M(i?1)→i為從時(shí)刻 (i?1)到時(shí)刻i的演化算符；M(i?1)→i中包含了有偏差的模式參數(shù) Λ。按照式（1）可以遞推得到如下方程

假設(shè)初始場X0與驅(qū)動場Fi都是準(zhǔn)確的，唯一不準(zhǔn)確的是模式參數(shù) Λ，在 Λ上疊加一個小擾動 λ?得到另外的一個模式參數(shù)為可以預(yù)期，只要積分步長不是很長，受新的模式參數(shù) λ影響，各個時(shí)刻狀態(tài)變量可以表示為原狀態(tài)變量Xi與擾動量相加的形式Xi+則可以得到擾動量所滿足的方程

和

于是演化公式可改為

注意到，此處假設(shè)初始場和外界強(qiáng)迫場都是準(zhǔn)確的，即初始場和外界強(qiáng)迫場都不擾動，在這種情況下，狀態(tài)變量的擾動都來自于模式參數(shù)的擾動。

定義廣義背景場誤差協(xié)方差矩陣為

式中，矩陣實(shí)際上包含了切線性演化算符，其轉(zhuǎn)置矩陣則包含了伴隨算符信息，于是，

如果只通過擾動參數(shù)構(gòu)造集合，則可以在最小二乘意義下，顯式地求解出式（7）矩陣的具體形式，即

如果只優(yōu)化模式參數(shù)，則四維變分的目標(biāo)函數(shù)為

式中，Yi為第i時(shí)刻的觀測；Hi為第i時(shí)刻的觀測投影算符；Ri為第i時(shí)刻的觀測誤差協(xié)方差矩陣。

傳統(tǒng)四維變分對模式參數(shù)進(jìn)行最優(yōu)估計(jì)的時(shí)候，就是將目標(biāo)函數(shù)對模式參數(shù) Λ求梯度，并將目標(biāo)函數(shù)值和梯度值代入到最優(yōu)化算法中，通過線性搜索和逐步迭代得到最優(yōu)的模式參數(shù)。仿照這一過程，在當(dāng)前模式參數(shù) Λ附近進(jìn)行擾動展開，推導(dǎo) Λ附近使得目標(biāo)函數(shù)取極小值的最優(yōu)擾動量，設(shè)擾動量為δΛ，將在當(dāng)前模式參數(shù) Λ 附近進(jìn)行展開，由于擾動量是小量，忽略高階項(xiàng)后，可以得到

將式（12）代入式（10），則四維變分目標(biāo)函數(shù)為

接著求出式（13）中目標(biāo)函數(shù)相對于擾動量δΛ的梯度，

令式（14）梯度為0，即?J(δΛ)=0，則可以求出δΛ的解析表達(dá)式

需要指出的是，這一最優(yōu)增量是在 Λ為初猜場的情況下得到的，對于非線性動力系統(tǒng)，只要稍微有一點(diǎn)改動，上述最優(yōu)增量的值就必然會發(fā)生變化。因此，為了讓這一方法具有一定的穩(wěn)定性，需要仿照傳統(tǒng)四維變分的做法，引入線性搜索過程，即在相空間的δΛ方向上，以一定的步長 β線性搜索（例如采用二分法）使得目標(biāo)函數(shù)在該方向上取極小的最優(yōu)值 β ·δΛ作為訂正量的最優(yōu)取值，將 Λ +β·δΛ作為下一次迭代的猜測值代入目標(biāo)函數(shù)，繼續(xù)計(jì)算解析解，循環(huán)往復(fù)上述過程從而獲得最優(yōu)的模式參數(shù)。

3 數(shù)值試驗(yàn)

1963年美國麻省理工學(xué)院的氣象學(xué)家Lorenz在對天氣預(yù)報(bào)動力學(xué)模型進(jìn)行數(shù)值計(jì)算時(shí)，發(fā)現(xiàn)了一個由非線性微分方程組所描述的著名的 Lorenz方程[18]該方程是一個3個變量模型，是對流非周期模式的簡化，同時(shí)也是混沌理論的原型個例[10]。由于Lorenz方程具有較強(qiáng)的非線性，且模型較為簡潔，所以常采用orenz-63模型來檢驗(yàn)數(shù)據(jù)同化算法的性能[19–20]。

Lorenz方程為

式中， σ >0，b>0； 參數(shù) σ 、r、b、t分別 表示Prantl數(shù)、Rayleigh數(shù)、對流尺度聯(lián)系的參數(shù)和時(shí)間；x是對流強(qiáng)度；y是 最大溫度差；z是對流引起的層變化[21]。模型標(biāo)準(zhǔn)參數(shù)一般為σ =10,r=28,b=8/3，從而產(chǎn)生混沌解[22]。

本文采用Lorenz-63模型進(jìn)行數(shù)值試驗(yàn)。首先，對標(biāo)準(zhǔn)參數(shù)模型采用四階榮格庫塔（Runge-Kutta）方法積分生成真實(shí)場（積分步長設(shè)置為 dt=0.01）；隨后在一個同化時(shí)間窗口內(nèi)以一定的采樣間隔進(jìn)行觀測采集，并在觀測數(shù)據(jù)上加入方差為1的白噪聲，生成帶噪聲的觀測數(shù)據(jù)。具體優(yōu)化步驟為：根據(jù)式（15）計(jì)算每次優(yōu)化迭代之后的參數(shù)最優(yōu)擾動方向；確定最優(yōu)擾動步長，將計(jì)算出的最優(yōu)擾動疊加到當(dāng)前參數(shù)值上得到下一次優(yōu)化迭代的參數(shù)初值；經(jīng)過多次迭代，得到收斂的模型參數(shù)值，作為最終參數(shù)優(yōu)化結(jié)果。本文設(shè)置單參數(shù)、雙參數(shù)和三參數(shù)試驗(yàn)來檢測新方法的模型參數(shù)優(yōu)化能力，并與傳統(tǒng)4DVar方法進(jìn)行模型參數(shù)優(yōu)化對比試驗(yàn)；此外，通過設(shè)置不同的模型參數(shù)真值、集合成員個數(shù)、同化時(shí)間窗口長度及觀測采樣間隔 等敏感性試驗(yàn)，考察新方法的參數(shù)估計(jì)能力。

3.1 新方法參數(shù)優(yōu)化結(jié)果

本節(jié)試驗(yàn)積分窗口取500步，集合成員個數(shù)取500個，將 σ =10,r=28,b=8/3這些特殊參數(shù)值作為本次試驗(yàn)的標(biāo)準(zhǔn)參數(shù)值，取背景場狀態(tài)變量為x0=1,y0=2,z0=3對模式進(jìn)行積分構(gòu)造真實(shí)場，以10步為采樣間隔采集觀測并加入白噪聲構(gòu)造觀測場。假設(shè)初猜的參數(shù)值為σ1=9,r1=25,b1=5，并采用相同的背景場狀態(tài)變量向前積分500步。傳統(tǒng)4DVar方法采用了同樣的參數(shù)設(shè)置，編輯伴隨模式計(jì)算代價(jià)函數(shù)梯度并采用同樣的線性搜索方法優(yōu)化迭代計(jì)算。表1為本試驗(yàn)的主要參數(shù)配置。

表1 A-4DEnVar和傳統(tǒng)4DVar對Lorenz-63模型的參數(shù)優(yōu)化試驗(yàn)設(shè)計(jì)Table 1 Experimental design of Lorenz-63 model parameter optimization by A-4DEnVar and traditional 4DVar

本節(jié)設(shè)計(jì)7種試驗(yàn)方案，包括單參數(shù)試驗(yàn)、雙參數(shù)試驗(yàn)以及三參數(shù)試驗(yàn)用以檢驗(yàn)新方法的模型參數(shù)優(yōu)化能力。結(jié)果表明，針對不同個數(shù)的參數(shù)，新方法均能收斂到真值，對三參數(shù)優(yōu)化后的模型分析場與真實(shí)場如圖1所示（其余結(jié)果未列出）。由圖1可知，A-4DEnVar參數(shù)優(yōu)化方法和傳統(tǒng)的4DVar參數(shù)優(yōu)化方法均可達(dá)到一個理想的優(yōu)化效果，分析場擬合真實(shí)場程度高。各試驗(yàn)中代價(jià)函數(shù)隨迭代變化情況如圖2所示。由圖2可知，A-4DEnVar在優(yōu)化迭代過程中代價(jià)函數(shù)值在100步左右可以實(shí)現(xiàn)收斂，收斂結(jié)果較為穩(wěn)定，因而具有良好的模型參數(shù)優(yōu)化能力。

圖1 三參數(shù)同時(shí)優(yōu)化試驗(yàn)結(jié)果Fig.1 Results of three parameters optimization experiment

圖2 A-4DEnVar代價(jià)函數(shù)值隨迭代次數(shù)的變化Fig.2 Changes of value of cost function of A-4DEnVar with the number of iterations

3.2 同化時(shí)間窗口長度和觀測采樣間隔影響

為了測試同化時(shí)間窗口長度和觀測采樣間隔對新方法在模型參數(shù)優(yōu)化方面能力的影響，本節(jié)采用與3.1節(jié)相同的初始場和初猜參數(shù)，選取4組不同的同化時(shí)間窗口長度進(jìn)行試驗(yàn)，并計(jì)算分析場與真實(shí)場的均方根誤差（RMSE）。同化時(shí)間窗口長度為：100步、200步、300步和400步，在這4個時(shí)間窗口下再分別取5步、10步和20步的觀測采樣間隔（Observation Sampling Interval，OSI）。為方便對比對比，傳統(tǒng) 4DVar方法采用了同樣的參數(shù)設(shè)置進(jìn)行試驗(yàn)。

圖3給出了改變積分時(shí)間窗口長度和觀測采樣間隔后的模型變量的平均RMSE變化。試驗(yàn)結(jié)果表明，針對不同的同化窗口和采樣間隔，A-4DEnVar方法均可實(shí)現(xiàn)模型參數(shù)的理想優(yōu)化，效果與傳統(tǒng)4DVar方法相當(dāng)；新方法的模型參數(shù)優(yōu)化速度較傳統(tǒng)的4DVar方法快，這是由于A-4DEnVar方法采用了解析解而不是梯度下降方向作為搜索方向。圖4描繪了不同的同化時(shí)間窗口和觀測采樣間隔下優(yōu)化迭代收斂后的RMSE。結(jié)果表明，采用新方法與傳統(tǒng)4DVar方法均能使分析場RMSE低于觀測誤差標(biāo)準(zhǔn)差，可實(shí)現(xiàn)模型參數(shù)的收斂，同一種方法在同一采樣間隔的情況下，隨著積分窗口長度的增長，最終收斂的平均RMSE會增大；同一種方法在同一積分窗口長度的情況下，隨著采樣間隔的增大，最終收斂的平均RMSE會增大。

圖3 改變積分時(shí)間窗口長度和觀測采樣間隔之后的模型均方根誤差變化Fig.3 Model RMSE changes after changing integral time window length and observation sampling interval

圖4 不同積分時(shí)間窗口長度和觀測采樣間隔組合下收斂后的平均均方根誤差Fig.4 Convergent average RMSE under the combination of different integral time window length and observation sampling interval

3.3 集合成員數(shù)量影響

為了測試集合成員數(shù)量的大小對新方法參數(shù)優(yōu)化能力的影響，本節(jié)采用與3.1節(jié)相同的初始場和初猜參數(shù)及同化時(shí)間窗口和參數(shù)標(biāo)準(zhǔn)值，選取4組不同的集合成員數(shù)量進(jìn)行試驗(yàn)。4組集合成員數(shù)量分別是：10、100、200和300。試驗(yàn)參數(shù)如表3所示。

表2 改變積分窗口長度和觀測采樣間隔條件下Lorenz-63模型參數(shù)優(yōu)化試驗(yàn)Table 2 Lorenz-63 model parameter optimization experiment under the changes of integral window length and observation sampling interval

表3 不同集合成員數(shù)量情況下解析4DEnVar對Lorenz-63模型的參數(shù)優(yōu)化Table 3 Lorenz-63 parameter optimization experiment by A-4DEnVar under different number of members of the set

圖5展示了不同集合成員數(shù)量情況下模型分析場與真實(shí)場的擬合程度。由圖5可知，在不同的集合成員數(shù)量下，新方法均可實(shí)現(xiàn)參數(shù)的理想收斂，模型分析場擬合真實(shí)場程度高，這說明新方法對集合成員數(shù)量不敏感；本節(jié)最小集合成員個數(shù)取的是10，若再取更小的集合成員個數(shù)（例如3或4），則需要更多的迭代次數(shù)來實(shí)現(xiàn)模型參數(shù)的收斂。

圖5 改變集合成員數(shù)量后模型分析場擬合真實(shí)場軌跡Fig.5 Trajectory chart of analysis field fitting real field after changing the number of members of the set

3.4 參數(shù)標(biāo)準(zhǔn)值影響

為了測試不同參數(shù)標(biāo)準(zhǔn)值對提出方法的影響，本節(jié)采用與3.1節(jié)相同的初始場和初猜參數(shù)及同化時(shí)間窗口和集合成員個數(shù)，在模式原標(biāo)準(zhǔn)參數(shù)上加上一定的擾動倍數(shù)，形成3組不同的參數(shù)標(biāo)準(zhǔn)值進(jìn)行組合試驗(yàn)[10]。3種參數(shù)所屬區(qū)間分別是： σ ∈(0.9×10,1.1×10)，r∈(0.9×28,1.1×28)和b∈(0.9×8/3,1.1×8/3)，每個區(qū)間生成10個任意值，一共1 000組試驗(yàn)。試驗(yàn)參數(shù)如表4所示。

表4 不同參數(shù)標(biāo)準(zhǔn)值情況下解析4DEnVar對Lorenz-63模型的參數(shù)優(yōu)化Table 4 Lorenz-63 parameter optimization experiment by A-En4DVar under different standard values of parameters

圖6展示了改變模型參數(shù)標(biāo)準(zhǔn)值后的各組代價(jià)函數(shù)值隨著迭代次數(shù)的變化，橫坐標(biāo)表示組別，縱坐標(biāo)表示代價(jià)函數(shù)值。圖7展示了改變模型參數(shù)標(biāo)準(zhǔn)值后平均代價(jià)函數(shù)值隨著迭代次數(shù)的變化，橫坐標(biāo)表示迭代次數(shù)，縱坐標(biāo)表示代價(jià)函數(shù)值。由圖6和圖7可知，對選取的1 000組參數(shù)標(biāo)準(zhǔn)值，在迭代50步之內(nèi)，代價(jià)函數(shù)值即可收斂至最小值，此時(shí)模型參數(shù)收斂至標(biāo)準(zhǔn)值，由此可知新方法優(yōu)化能力較強(qiáng)，且對參數(shù)的選取不敏感。

圖6 改變模型參數(shù)標(biāo)準(zhǔn)值后各組代價(jià)函數(shù)值變化Fig.6 Changes of each group's cost functions after changing the standard values of model parameters

圖7 改變模型參數(shù)標(biāo)準(zhǔn)值后平均代價(jià)函數(shù)值變化Fig.7 Changes of average cost function value after changing the standard values of model parameters

4 結(jié)論

本文將模型參數(shù)視為一種特殊的控制變量，在傳統(tǒng)4DVar框架下，以迭代得到的模式參數(shù)為基準(zhǔn)展開集合擾動，計(jì)算誤差協(xié)方差矩陣。在此基礎(chǔ)上，使用代價(jià)函數(shù)最小值的解析解構(gòu)造線性搜索最優(yōu)化方法，提出了A-4DEnVar參數(shù)優(yōu)化方法，避免了傳統(tǒng)4DVar方法中的伴隨模式的使用。為驗(yàn)證有效性，采用傳統(tǒng)的4DVar方法和新的A-4DEnVar參數(shù)優(yōu)化方法對Lorenz-63模型進(jìn)行參數(shù)優(yōu)化對比試驗(yàn)，針對不同個數(shù)的參數(shù)，新方法同化效果與傳統(tǒng)方法4DVar相當(dāng)。針對不同的同化時(shí)間窗口、觀測采樣間隔、集合樣本個數(shù)及不同標(biāo)準(zhǔn)參數(shù)的敏感性試驗(yàn)結(jié)果顯示，新方法能達(dá)到一個準(zhǔn)確的收斂效果，在較長的同化時(shí)間窗口和觀測采樣間隔時(shí)也可以實(shí)現(xiàn)理想的模型參數(shù)優(yōu)化效果，并且該方法對集合樣本個數(shù)和標(biāo)準(zhǔn)參數(shù)不敏感，采用較少的集合樣本即可達(dá)到同化效果，節(jié)約了工作量，這在實(shí)際模型的參數(shù)優(yōu)化工作中，將具有重要意義。