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

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

1 引言

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

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

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

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

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

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

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

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

和

于是演化公式可改為

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

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

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

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

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

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

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

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

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

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

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

3 數(shù)值試驗

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

Lorenz方程為

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

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

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

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

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

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

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

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

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

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

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

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

圖4 不同積分時間窗口長度和觀測采樣間隔組合下收斂后的平均均方根誤差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ù)標準值，選取4組不同的集合成員數(shù)量進行試驗。4組集合成員數(shù)量分別是：10、100、200和300。試驗參數(shù)如表3所示。

表2 改變積分窗口長度和觀測采樣間隔條件下Lorenz-63模型參數(shù)優(yōu)化試驗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ù)量情況下模型分析場與真實場的擬合程度。由圖5可知，在不同的集合成員數(shù)量下，新方法均可實現(xiàn)參數(shù)的理想收斂，模型分析場擬合真實場程度高，這說明新方法對集合成員數(shù)量不敏感；本節(jié)最小集合成員個數(shù)取的是10，若再取更小的集合成員個數(shù)（例如3或4），則需要更多的迭代次數(shù)來實現(xiàn)模型參數(shù)的收斂。

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

3.4 參數(shù)標準值影響

為了測試不同參數(shù)標準值對提出方法的影響，本節(jié)采用與3.1節(jié)相同的初始場和初猜參數(shù)及同化時間窗口和集合成員個數(shù)，在模式原標準參數(shù)上加上一定的擾動倍數(shù)，形成3組不同的參數(shù)標準值進行組合試驗[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組試驗。試驗參數(shù)如表4所示。

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

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

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

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

4 結論

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