集合預(yù)報誤差在GRAPES 全球四維變分同化中的應(yīng)用研究Ⅰ：局地化方案設(shè)計與動力平衡特征分析*

龔建東 王瑞春

GONG Jiandong WANG Ruichun

1. 國家氣象中心，北京， 100081

2. 中國氣象局?jǐn)?shù)值預(yù)報中心，北京， 100081

1. National Meteorological Centre，Beijing 100081，China

2. Numerical Weather Prediction Center of CMA，Beijing 100081，China

1 引 言

背景誤差協(xié)方差簡稱背景誤差，它的準(zhǔn)確估計直接影響資料同化的分析性能。背景誤差一般使用氣候歷史樣本進(jìn)行統(tǒng)計，并考慮到誤差協(xié)方差結(jié)構(gòu)的復(fù)雜性，采用高度模型化的誤差自相關(guān)、協(xié)相關(guān)等模型和各向同性均質(zhì)誤差假設(shè)來簡化協(xié)方差結(jié)構(gòu)，使得背景誤差可以顯式表達(dá)并用于資料同化系統(tǒng)中（Lorenc，1986；Bannister，2008）。這樣處理的不足是在變分資料同化的每個分析循環(huán)中，雖然數(shù)值預(yù)報背景場在不斷更新，誤差特征也在不斷變化，但同化框架中背景誤差的特征卻固定不變。導(dǎo)致背景誤差在不同時刻不斷變化的因素較多，但其明顯受到所在地天氣尺度系統(tǒng)的影響，天氣系統(tǒng)的斜壓不穩(wěn)定成為背景誤差增長的主導(dǎo)因素，呈現(xiàn)出隨不同天氣流型變化的“流依賴”特征（Kalnay，2003，第5.6 節(jié)）。如果能在同化框架中引入具備流依賴特征的背景誤差，將會幫助同化框架更加有效地確定背景場和觀測的相對權(quán)重，也能更加合理地將觀測信息在空間上傳遞出去，從而提高分析質(zhì)量和預(yù)報技巧。

隨著集合卡爾曼濾波（EnKF，ensemble Kalman filtering）等集合資料同化技術(shù)的快速發(fā)展與逐步成熟（Evensen，1994；Anderson，2001；Houtekamer，et al，2005），利用集合樣本的短期預(yù)報來估計集合預(yù)報誤差協(xié)方差（簡稱集合預(yù)報誤差）作為背景誤差更為合理。從使用集合預(yù)報誤差角度分析，在現(xiàn)有資料同化方法中，三維變分資料同化（3DVar）使用氣候背景誤差協(xié)方差。4DVar 在同化時間窗起始時刻使用氣候背景誤差協(xié)方差，并利用切線性和伴隨模式在同化時間窗內(nèi)實現(xiàn)背景誤差的時間演變?；?DVar，通過擾動觀測資料和海溫、疊加隨機物理過程等來表示觀測資料、海洋下邊界條件及模式的不確定性，產(chǎn)生多組分析和預(yù)報的集合樣本來估計集合預(yù)報誤差，這種方法被稱為集合方式的四維變分資料同化方法，簡稱集合資料同化（EDA，ensemble of data assimilations）方 法。在4DVar 同化時間窗的起始時刻，同時使用氣候背景誤差和EnKF 或EDA 方法生成集合樣本估計的集合預(yù)報誤差，使得在同化時間窗的起始時刻就具有隨天氣形勢演變特征的流依賴動態(tài)背景誤差，這種方法被稱為集合四維變分同化（En4DVar，ensemble four dimensional variational data assimilation），其與4DVar 的主要差別在于是否使用集合預(yù)報誤差。而在同化時間窗內(nèi)，使用集合預(yù)報在時間序列上的多組集合樣本來直接估計整個同化時間窗內(nèi)不同時刻的背景誤差，這種方法被稱為四維集合變分 同 化（ 4DEnVar， four dimensional ensemble variational data assimilation），其與4DVar 和En4DVar的主要區(qū)別是不再引入切線性和伴隨模式。熊春暉等（2013）、Fairbairn 等（2014）對這些方法的差異做了詳細(xì)介紹。

目前，國際上領(lǐng)先的業(yè)務(wù)變分資料同化系統(tǒng)大多已實現(xiàn)從靜態(tài)的、表示氣候特征的背景誤差向基于集合樣本估計的、能表示隨天氣形勢演變特征的流依賴動態(tài)背景誤差的升級（熊春暉等，2013）。英國氣象局采用在目標(biāo)函數(shù)中增加擴展控制變量的方式將集合轉(zhuǎn)置卡爾曼濾波（ETKF）估計的集合預(yù)報誤差與氣候背景誤差結(jié)合，發(fā)展了En4DVar 系統(tǒng)（Clayton，et al，2013），加拿大氣象局也發(fā)展了EnKF與4DVar 結(jié)合的En4DVar 同化系統(tǒng)（Buehner，et al，2010a，2010b）。業(yè)務(wù)數(shù)值預(yù)報中心的試驗結(jié)果表明，考慮集合預(yù)報誤差改善了4DVar 的分析質(zhì)量，特別是在南半球、熱帶地區(qū)，并對臺風(fēng)路徑預(yù)報也有顯著改進(jìn)。歐洲中期天氣預(yù)報中心與法國氣象局使用基于4DVar 的集合資料同化方法來生成集合樣本，采用譜濾波來濾除噪聲以保留有用的天氣系統(tǒng)變化的集合預(yù)報誤差并用于4DVar，顯著地改善了分析和模式預(yù)報質(zhì)量，且正貢獻(xiàn)在全球范圍10 d 預(yù)報中都有效（Bonavita，et al，2012）。此外，英國氣象局（Lorenc，et al，2015；Bowler，et al，2017）、美國環(huán)境預(yù)報中心（NCEP）還進(jìn)一步采用擴展控制變量方式將集合預(yù)報誤差用于全球4DEnVar（Wang，et al，2013；Kleist，et al，2015a，2015b），加拿大氣象局也發(fā)展并業(yè)務(wù)應(yīng)用了4DEnVar 系統(tǒng)（Buehner，et al，2013）。

應(yīng)用集合預(yù)報誤差雖然可以通過不同方式實現(xiàn)，但事實證明這些方式之間是等價的（Buehner，2005；Wang，et al，2007）。可以看到，在上述業(yè)務(wù)實踐中，Lorenc（2003a，2003b）提出的增加擴展控制變量方法是引入集合預(yù)報誤差的常用方案。增加擴展控制變量可以使得分析增量有更大的空間自由度，從而幫助更好地引入觀測信息。使用集合樣本來估計集合預(yù)報誤差有優(yōu)勢，但其不足也很明顯。Lorenc（2003a）指出，當(dāng)集合樣本數(shù)較少時預(yù)報誤差容易被低估，且存在虛假遙相關(guān)，需要進(jìn)行誤差協(xié)方差的局地化，而局地化容易破壞集合樣本中所包含的動力平衡關(guān)系（Kepert，2009）。

GRAPES 于2018 年6 月實現(xiàn)了全球4DVar 業(yè)務(wù)運行（Zhang，et al，2019），中國成為世界上為數(shù)不多的采用四維變分同化的全球業(yè)務(wù)預(yù)報中心之一。4DVar 方法的基本原理表明，雖然其能在同化時間窗內(nèi)通過切線性與伴隨模式積分來隱式計算流依賴背景誤差，但其初始時刻的背景誤差一直是氣候態(tài)的。要克服這一缺陷，需要在4DVar 中引入集合預(yù)報誤差，發(fā)展GRAPES 全球En4DVar 技術(shù)，從而實現(xiàn)在同化時間窗起始時刻就應(yīng)用流依賴背景誤差。

在GRAPES 全球4DVar 基礎(chǔ)上進(jìn)行拓展，有效使用集合樣本估計的集合預(yù)報誤差來改進(jìn)全球分析質(zhì)量是本研究的重點。文中主要研究這幾個方面的問題：首先，GRAPES 全球4DVar 采用增量分析方案（Courtier，et al，1994），在具體實施上有自己的特殊性（薛紀(jì)善等，2008），有必要研究在該系統(tǒng)框架中增加擴展控制變量的科學(xué)方案。雖然已有一些研究基于GRAPES 區(qū)域3DVar 開展了使用集合預(yù)報誤差的工作（Chen，et al，2015），但全球4DVar 中尚未開展相關(guān)研究。其次，通過增加集合樣本數(shù)來降低抽樣誤差，并對集合預(yù)報誤差給予較大權(quán)重是各業(yè)務(wù)中心同化技術(shù)發(fā)展的趨勢?？紤]到擴展控制變量的空間維數(shù)與集合樣本數(shù)有關(guān)，集合樣本數(shù)越大對應(yīng)的擴展控制變量的空間維數(shù)越大。在大于100 個集合樣本數(shù)的條件下，如何在GRAPES 全球4DVar 中既能有效使用集合預(yù)報誤差，又使得擴展控制變量的維數(shù)增加受到限制，以減少4DVar 系統(tǒng)的計算與存儲消耗。最后在GRAPES全球4DVar 引入擴展控制變量并對其進(jìn)行局地化時，需要研究局地化對分析場平衡造成的影響，并發(fā)展減緩局地化影響的辦法。

2 GRAPES 全球4DVar 背景誤差公式表達(dá)

GRAPES 全球模式為水平經(jīng)緯度格點模式。GRAPES 全球4DVar 的氣候背景誤差矩陣（B）可以表示為 B=B1/2BT/2≈UUT。B1/2是B 的平方根矩陣，在具體實現(xiàn)時通過一系列誤差結(jié)構(gòu)函數(shù)（或變量變換矩陣）來完成。從極小化求解的控制變量（ v）到分析變量（ δX1）的變換可以表示為

式中， δX1=(δu,δv,δπ,δq)T分別是緯向風(fēng)與經(jīng)向風(fēng)、無量綱氣壓及比濕的分析增量，為 n維空間向量。v=(vψ,vχu,vπu,vq)T為控制變量，包括流函數(shù)、非平衡速度勢函數(shù)、非平衡無量綱氣壓及比濕。 Uh是水平變量變換，在全球模式中采用水平譜濾波形式來模擬背景誤差水平相關(guān)的各向同性分布特征，實現(xiàn)分析增量信息由觀測點向周圍模式格點的傳播。水平譜濾波需要采用水平高斯網(wǎng)格，與水平經(jīng)緯度網(wǎng)格不同，因而在GRAPES 全球4DVar 中包含水平經(jīng)緯度網(wǎng)格與水平高斯網(wǎng)格兩套網(wǎng)格。 Uv為垂直變量變換，采用經(jīng)驗正交函數(shù)分解（EOF）方式表達(dá)背景誤差垂直相關(guān)結(jié)構(gòu)。 εb為對角矩陣，表示非平衡部分的背景誤差均方根誤差。式（1）中水平譜濾波和經(jīng)驗正交函數(shù)分解的非平衡分析變量的具體實現(xiàn)過程可由下式表示

式中， δXu是 非平衡分析變量， E、 Λ是背景誤差垂直相關(guān)矩陣（ Rv）通過經(jīng)驗正交函數(shù)分解求得的特征向量與特征值，其表達(dá)與模式垂直層次有關(guān)，而與格點空間或譜空間無關(guān)，可作用在每個譜系數(shù)上。S 、 S?1分別是從格點空間到譜空間及從譜空間到格點空間的變換。 Bs是譜空間的水平相關(guān)系數(shù)矩陣。Wg是 勒讓德變換權(quán)重系數(shù)。式（2）表明， v是垂直特征模態(tài)上的投影系數(shù)。 G1s→d是從控制變量（ v）所在的水平高斯格點到水平經(jīng)緯度格點的插值算子。

式（1）中 UK表示背景誤差中變量間的協(xié)相關(guān)部分，也即平衡算子，可表示為

式中， N是流函數(shù)對無量綱氣壓解釋部分的平衡算子，通過風(fēng)壓線性平衡方程或平衡方程與統(tǒng)計回歸求解兩種方法結(jié)合的方式給出（王瑞春等，2015）。M是流函數(shù)對勢函數(shù)解釋部分的平衡算子，通過統(tǒng)計方式給出。通過 UK變量變換將相互獨立非平衡分 析 變 量 (δψ,δχu,δπu,δq)T變 換 為 分 析 變 量 全 量(δψ,δχ,δπ,δq)T。考慮到實際觀測及模式預(yù)報變量是風(fēng)場，因而還需進(jìn)行物理變換 UP，從流函數(shù)、速度勢函數(shù)變換為緯向風(fēng)與經(jīng)向風(fēng)。此外， UP還包括由靜力平衡關(guān)系從無量綱氣壓（ δπ ）診斷得出位溫（ δθ）

式中， cp是摩爾定壓熱容，g 是重力加速度。

式（2）中經(jīng)驗正交函數(shù)分解和水平譜濾波的參數(shù)由背景誤差垂直相關(guān)模型和水平相關(guān)模型決定。誤差垂直相關(guān)矩陣（ Rv）的構(gòu)造由垂直相關(guān)模型（ ρl,k）決定，與垂直層次所在的高度有關(guān)，表示為

式中，zl、zk是任意模式層l、k 的平均高度。K 是經(jīng)驗系數(shù)，當(dāng)K=0 時，所有垂直層次間的相關(guān)均為1，當(dāng)其數(shù)值增大時，垂直相關(guān)結(jié)構(gòu)逐步變得更加局地。垂直相關(guān)結(jié)構(gòu)也可以通過誤差樣本統(tǒng)計給出（王金成等，2014）。龔建東等（2020）在GRAPES全球4DVar 基礎(chǔ)上發(fā)展了EDA 方法生成多組集合樣本，也可以用于統(tǒng)計垂直相關(guān)結(jié)構(gòu)。

定義新息向量 dn=?HnM(Xb)， 其中是 tn時的觀測值，在高分辨率外循環(huán)進(jìn)行計算。定義 H為從模式狀態(tài)變量到觀測空間的觀測算子，其線性化算子表示為H。定義 M為非線性預(yù)報模式，其切線性模式表示為 M。于是GRAPES 全球4DVar 的目標(biāo)函數(shù)可以表示為

式中，R 是觀測誤差協(xié)方差。 t0與 tL分別為同化時間窗起始與終止時刻。通過極小化求解式（8）可以求得 t0時 刻分析增量 δX1。 Jc是目標(biāo)函數(shù)約束項，采用數(shù)字濾波方案約束重力波引發(fā)的不平衡（劉艷等，2019）。后文為公式表述簡潔而略去 Jc項。

由式（1）與（8）可知，在4DVar 的 t0時刻，分析增量（ δX1）利用氣候背景誤差（B）或變量變化矩陣（ U）只能獲得反映氣候背景誤差的分析增量，而在tn時 刻的分析增量表示為 δXn= MδX1，通過模式動力約束可以部分反映流依賴特征。

3 GRAPES 全球4DVar 引入集合預(yù)報誤差

3.1 擴展控制變量的引入

這時 αi是 δ對應(yīng)的局地化分析變量。受限于高性能計算機的存儲能力，引入局地化分析變量后需要考慮對其維數(shù)進(jìn)行控制。GRAPES 全球模式水平采用Arakawa C 網(wǎng)格，垂直方向采用Chaney-Philips模式分層，不同分析變量 (δu 、 δv 、 δπ 、 δq)所在的模式網(wǎng)格空間位置不同，如果對同一集合樣本 δ的不同分析變量都取同一個局地化分析變量 αi，并視質(zhì)量點δ π所在位置為共同位置，可以降低 αi的長度。

式中， C 為局地化矩陣， C1/2為其平方根矩陣，可通過變量變換矩陣 Uα構(gòu) 造。 Uα可進(jìn)一步表示為垂直局地化和水平局地化矩陣 Uα=。的作用與UhUv的作用一致，同樣可以采用水平譜濾波和垂直經(jīng)驗正交函數(shù)分解實現(xiàn)。于是式（9）表示為

3.2 局地化分析變量實現(xiàn)方式

在GRAPES 全球4DVar 中實現(xiàn)局地化分析變量有兩種方式，一是只進(jìn)行水平局地化而不進(jìn)行垂直局地化。此時，擴展控制變量（）蛻變?yōu)樗蕉S變量，的維數(shù)顯著減少，因而可以不考慮水平局地化的降維處理。也即，水平局地化所用的譜濾波的波數(shù)與4DVar 內(nèi)循環(huán)水平譜濾波的波數(shù)相同， vαi 定義在控制變量（ v）所在的水平高斯網(wǎng)格點上， αi定義在內(nèi)循環(huán)低分辨率水平經(jīng)緯度網(wǎng)格點上。這在GRAPES 全球4DVar 中表示為

式中， Bsα表示水平局地化變量在譜空間的水平相關(guān)系數(shù)矩陣。

另一種辦法是進(jìn)行三維局地化—同時考慮水平和垂直方向的局地化。此時，通過在式（13）的水平譜濾波中取更少波數(shù)，垂直經(jīng)驗正交函數(shù)分解中取前幾個主導(dǎo)特征向量來大幅度縮減控制變量（）的維數(shù)。在GRAPES 全球4DVar 中表示為

在實際應(yīng)用時，當(dāng)局地化所采用網(wǎng)格的水平分辨率與內(nèi)循環(huán)的分辨率不一致時，需要重新定義一組網(wǎng)格用于局地化變量，這會使得GRAPES 的網(wǎng)格定義更為復(fù)雜。一種替代的辦法是在水平局地化網(wǎng)格分辨率不變的情況下，僅在垂直局地化時使用較少的主導(dǎo)特征向量。Clayton 等（2013）的研究表明，當(dāng)經(jīng)驗正交函數(shù)分解前若干個特征值對總誤差的解釋比例達(dá)到95%時可以獲得滿意的結(jié)果。目前的GRAPES 全球4DVar 中，水平局地化與內(nèi)循環(huán)的網(wǎng)格一致，垂直方向按照解釋比達(dá)到90%來截取主導(dǎo)特征向量。試驗表明前8 個特征值占比已經(jīng)超過90%，后續(xù)試驗中取前8 個特征模態(tài)進(jìn)行垂直局地化。

3.3 GRAPES 全球En4DVar 公式表達(dá)

采用增加擴展控制變量，是在具有氣候背景誤差特征的分析增量上增加擴展控制變量給出的具備隨流型演變的分析增量，最終獲得的分析增量是這兩部分增量的加權(quán)之和（Lorenc，2003b；Wang，et al，2007）

由新息向量定義，分析值在觀測空間投影（ yn）表示為

于是 yn與觀測（）的偏差表示為

考慮了擴展控制變量的目標(biāo)函數(shù)表示為

由式（20）及式（16）—（17），對控制變量（ v）和擴展控制變量（）的梯度表示為

對式（21）和（22）再次求偏導(dǎo)數(shù)，計算目標(biāo)泛函相對于控制變量的二階偏導(dǎo)數(shù)（公式略），構(gòu)成海森矩陣可以發(fā)現(xiàn)，將權(quán)重系數(shù) βc、 βe放在分析增量的計算中（式（16）、（17）），而不是放在目標(biāo)函數(shù)（式（20））背景誤差與集合預(yù)報誤差項上，主要是改進(jìn)海森矩陣的條件數(shù)，加快變分同化系統(tǒng)收斂。試驗表明，在GRAPES 全球4DVar 中若將權(quán)重放在背景誤差與集合預(yù)報誤差項上，當(dāng) βc、 βe取值差異大時將顯著影響極小化迭代收斂。

最后，在GRAPES 全球4DVar 中具體實現(xiàn)的步驟是：首先由外循環(huán)計算新息向量 dn，并按照下列步驟計算

（2）由控制變量（ v）通過式（16）計算 δX1，由擴展控制變量（,···,）通過式（17）計算 δX2；

（3）由式（19）計算分析值在觀測空間投影（ y）與觀測（ yo） 的偏差（?yn）；

（5）由式（21）計算 ?vJ 并由式（22）計算 ?vαiJ；從而獲得目標(biāo)函數(shù)對所有控制變量的梯度（ ?v,vα1,...,vαNJ）；

（7）回到第（2）步，進(jìn)入下一次循環(huán)。如精度滿足需求，則退出循環(huán)。

采用以上的計算策略時，通過循環(huán)集合樣本，全部樣本只需要存儲一個 αi即可，大幅度減少了對變量存儲的需求。然而這種處理需要采用循環(huán)計算的方式反復(fù)計算并依次求得每個 αi，在程序設(shè)計上更為復(fù)雜，這與二維空間變量的處理方式有明顯不同。

4 水平局地化對地轉(zhuǎn)平衡約束關(guān)系的影響

4.1 水平局地化尺度影響

GRAPES 全球4DVar 的氣候背景誤差的水平相關(guān)模型采用融合水平相關(guān)模型（龔建東等，2020），水平局地化模型采用高斯型相關(guān)函數(shù)。借鑒Xu 等（2008a，2008b）提出的在分析時刻前后利用集合樣本在時間序列上的擴展來增加集合樣本數(shù)的方法，試驗利用20 個NCEP 全球集合預(yù)報的6 與9 h 預(yù)報的集合樣本，并采用扣除對應(yīng)時間集合平均的擾動樣本來減緩因預(yù)報時長不一致帶來的大尺度場存在的相位誤差，總計40 個樣本。任選2018 年7 月9 日太平洋北部地區(qū)一次鋒面天氣過程進(jìn)行分析。由集合樣本估計出來的集合預(yù)報誤差如圖1 所示，在300 hPa 高度水平面上（圖1a、b，約為模式35 層）集合樣本估計的氣壓與風(fēng)場誤差最大的位置位于冷鋒槽及槽后位置，在沿54°N 的垂直剖面上（圖1c、d），氣壓誤差隨高度升高而遞減，槽后明顯大于槽前，緯向風(fēng)場的誤差主要位于風(fēng)速急流區(qū)附近。

將單個觀測資料放置在冷鋒鋒面后的西北平直氣流上，位于（54°N，173°W）的模式35 層（約250 hPa）位置（圖1）。取氣候背景誤差與集合預(yù)報誤差的權(quán)重分別為0.1 與0.9，這時分析增量主要體現(xiàn)集合預(yù)報誤差的作用，并考慮不同水平局地化尺度參數(shù)來研究其對分析增量的影響，試驗結(jié)果如圖2 所示。水平局地化采用高斯相關(guān)分布特征，在水平局地化尺度的4 倍距離的位置，水平相關(guān)系數(shù)基本降為0。當(dāng)水平相關(guān)尺度參數(shù)取值較大時（如28000 km），地球上最遠(yuǎn)大圓距離（20000 km）處的分析增量受局地化作用將會衰減到原值的70%，水平局地化基本不起作用，可以作為不進(jìn)行水平局地化的試驗結(jié)果。由于集合樣本數(shù)有限，直接由集合樣本估計會產(chǎn)生虛假的遙相關(guān)，在距觀測位置較遠(yuǎn)的距離上會產(chǎn)生虛假的分析增量（圖略）。在單點試驗中分析增量不僅分布在觀測點周圍，在遠(yuǎn)離觀測點的位置也有較大的分析增量（圖2a）。試驗發(fā)現(xiàn)當(dāng)局地化水平相關(guān)尺度在1500 km 以下時距離觀測點較遠(yuǎn)處的虛假分析增量明顯受到抑制。例如：圖2b給出的水平相關(guān)尺度在850 km 時的分析增量，由于所選單點在鋒面后且集合預(yù)報誤差的權(quán)重占主導(dǎo)，因而分析增量呈現(xiàn)出隨天氣流型演變分布的特征—即分析增量沿氣流方向延展較長，而在垂直氣流方向分析增量范圍受到壓縮。當(dāng)不考慮集合預(yù)報誤差，僅有氣候背景誤差時在同化時間窗起始（ t0）分析增量的特點如圖9a 所示。

圖1 40 個集合樣本估計的 t0時刻集合預(yù)報誤差 （色階）（a. 氣壓誤差，單位：hPa；b. 風(fēng)場誤差，單位：m/s；c. 氣壓誤差垂直剖面，單位：hPa；d. 風(fēng)速誤差垂直剖面，單位：m/s；a、b 中流線為模式背景風(fēng)場，c、d 中等值線分別為模式背景氣壓與緯向風(fēng)速；圖中“+”為單點觀測位置，水平位置為（54°N，173°W），位于模式面35 層）Fig. 1 Estimated ensemble forecast errors at time t0 with 40 ensemble samples （shaded）（a. pressure error，unit：hPa；b. wind error，unit：m/s；c. vertical cross section of pressure error，unit：hPa；d. vertical cross section of wind error，unit：m/s；streamlines in （a） and （b） are background wind field，contours in （c） and （d） are background pressure and zonal wind，respectively，the cross "+" denotes observation site，which is located at（54°N，173°W）and on model level 35）

圖2 局地化水平相關(guān)尺度（a. 28000 km，b. 850 km）對單點氣壓觀測產(chǎn)生的分析增量的影響（新息向量大小為1 hPa，觀測時間位于 t0時刻，觀測位置如圖1 所示，圖中等值線為氣壓，單位：hPa）Fig. 2 Impacts of correlation scale of horizontal localization on analysis increment when assimilating single-point pressure observations （a. 28000 km，b. 850 km；the innovation is 1 hPa at time t 0，the observation site is the same as that shown in Fig. 1；the contours are for pressure，unit：hPa）

4.2 非平衡變量空間局地化

Mitchell 等（2002）、Clayton 等（2013）研究表明，直接對水平風(fēng)場及氣壓進(jìn)行局地化時，由于集合樣本中的水平風(fēng)與氣壓變量之間已經(jīng)包含大尺度準(zhǔn)地轉(zhuǎn)平衡約束關(guān)系，進(jìn)行水平局地化容易破壞這種關(guān)系。以緯向風(fēng)（ δu） 為例，其地轉(zhuǎn)分量（ δug）與氣壓（ δp）的關(guān)系可以表示為

式中，f 為科里奧利參數(shù)， ρ為空氣密度。當(dāng)對無量綱氣壓和風(fēng)場分量進(jìn)行局地化時，等式右側(cè)包含了對局地化函數(shù)的偏導(dǎo)數(shù)項而左側(cè)沒有，從而破壞對準(zhǔn)地轉(zhuǎn)平衡關(guān)系的表述。Clayton 等（2013）研究表明，如果對相互獨立的分析變量或非平衡變量進(jìn)行局地化時，將減緩局地化造成的影響。在GRAPES全球4DVar 中，先對 δ進(jìn)行物理逆變換，再將其中所包含的準(zhǔn)地轉(zhuǎn)平衡部分抽取出來，僅保留變量間不相關(guān)的非平衡部分。

TP(δu,δv,δπ,δq)T主要實現(xiàn)從分析變量 到分析變量(δψ,δχ,δπ,δq)T的 物 理 逆 變 換。 TK變 換 是 抽 取(δψ,δχ,δπ,δq)T中流函數(shù)與勢函數(shù)平衡關(guān)系，以及流函數(shù)與無量綱氣壓平衡關(guān)系，給出相互獨立的非平衡分析變量 (δψ,δχu,δπu,δq)T。此時對非平衡分析變量進(jìn)行水平局地化，由于平衡部分已經(jīng)被抽取出來，就不會影響平衡特征。在水平局地化后，重新使用 UK變量變換，計算流函數(shù)對應(yīng)的氣壓平衡部分。于是式（11）可以表示為

圖3 水平局地化對分析增量準(zhǔn)地轉(zhuǎn)平衡的影響 （a、d. 分析變量空間分布；b、e. 非平衡分析變量空間分布，局地化水平尺度為850 km；c、f. 分析變量空間分布，水平局地化尺度為28000 km；a、b、c. t 0 ，d、e、f. tL）Fig. 3 Impacts of horizontal localization on geostrophic balance of analysis increments （a，d. localization on analysis variables；b， e. unbalanced analysis variables，and the horizontal correlation scale of localization is 850 km；c，f. analysis variables，and the scale of localization is 28000 km；a，b，c. t0 ，d，e，f. tL）

圖3 給出了分別在分析變量空間與非平衡分析變量空間進(jìn)行局地化時，對分析平衡影響的差異。為方便比較起見，同時給出水平局地化相關(guān)尺度28000 km 的結(jié)果，作為不進(jìn)行水平局地化的對比試驗。由圖3c、f 可見，當(dāng)不進(jìn)行局地化時，同化時間窗起始時刻（ t0）的地面氣壓分析增量與Z=35 層（約250 hPa）的分析增量對應(yīng)，但隨著切線性模式積分，在同化時間窗終止時刻（ tL）分析增量大值區(qū)向東北移動，同時在鋒面前端出現(xiàn)接近?0.6 hPa 的氣壓增量。由圖1c 可知，集合樣本估計的集合預(yù)報誤差在模式低層的數(shù)值遠(yuǎn)大于35 層，單點分析增量在模式低層出現(xiàn)比觀測位置更強的分析增量。而當(dāng)水平局地化尺度減小時，以水平局地化相關(guān)尺度850 km 的結(jié)果為例，在同化時間窗起始時刻（ t0），模式第一層同樣出現(xiàn)與高層分析增量對應(yīng)的分析增量大值區(qū)，但在同化時間窗終止時刻（ tL），在鋒面前端出現(xiàn)?3 hPa 的氣壓增量，并在遠(yuǎn)離觀測位置區(qū)域出現(xiàn)大片區(qū)域的氣壓增量。該氣壓增量圍繞觀測位置以環(huán)形分布，其與觀測位置的距離隨水平局地化尺度變化。進(jìn)一步分析發(fā)現(xiàn)其位置與高斯水平相關(guān)模型二階導(dǎo)數(shù)的最大值處對應(yīng)。高斯相關(guān)模型的二階導(dǎo)數(shù)是表征氣壓與風(fēng)場的準(zhǔn)地轉(zhuǎn)耦合關(guān)系，分析圖3d 表明，在分析變量空間進(jìn)行局地化將產(chǎn)生虛假的分析增量。而當(dāng)在非平衡分析變量空間進(jìn)行局地化時（圖3b、e），可以發(fā)現(xiàn)其結(jié)果與不進(jìn)行水平局地化的結(jié)果相似，分析增量正值區(qū)向東北方向移動，同樣在鋒面前出現(xiàn)?0.6 hPa 的分析增量。上述結(jié)果表明將集合預(yù)報誤差引入GRAPES 全球4DVar 時，采用在非平衡變量空間進(jìn)行局地化來保證氣壓與風(fēng)場的準(zhǔn)地轉(zhuǎn)耦合關(guān)系非常重要。

5 三維局地化對分析平衡的影響

5.1 垂直局地化特征尺度

GRAPES 全球4DVar 采用靜力平衡分析方案，根據(jù)靜力平衡關(guān)系通過式（5）由無量綱氣壓增量計算虛位溫增量。垂直局地化的相關(guān)結(jié)構(gòu)可以采用式（6）和（7）所給出的垂直相關(guān)模型來實現(xiàn)，也可以直接利用模式預(yù)報的統(tǒng)計樣本來估計。王金成等（2014）比較了垂直相關(guān)模型與同一檢驗時刻的不同預(yù)報時效統(tǒng)計樣本估計出的垂直相關(guān)結(jié)構(gòu)的差異，結(jié)果表明統(tǒng)計樣本估計的垂直相關(guān)結(jié)構(gòu)更加合理。文中采用EDA 方法生成的集合樣本來估計垂直相關(guān)結(jié)構(gòu)（龔建東等，2020）。由集合樣本估計的垂直相關(guān)結(jié)構(gòu)（）存在因樣本數(shù)有限造成的虛假遙相關(guān)，文中利用相關(guān)范圍較寬的垂直相關(guān)模型（ Rv）對距離遠(yuǎn)的統(tǒng)計噪音進(jìn)行衰減限制，表示為Rv?， “ ?”表示“Schur”乘積，結(jié)果如圖4b 所示。采用垂直相關(guān)模型給出的垂直相關(guān)結(jié)構(gòu)（ Rv）在不同層次變化較為平滑，而集合樣本統(tǒng)計獲得的垂直相關(guān)結(jié)構(gòu)（）在邊界層相關(guān)較強，在邊界層之上對流層中層垂直相關(guān)迅速減弱，到對流層中高層垂直相關(guān)又進(jìn)一步變強?？傮w上，統(tǒng)計的垂直相關(guān)結(jié)構(gòu)明顯窄于相關(guān)模型。借鑒Clayton 等（2013）的方法，使用流函數(shù)的垂直相關(guān)結(jié)構(gòu)進(jìn)行所有變量的垂直局地化，垂直局地化結(jié)構(gòu)表示為 Cs=(Rv?)P，這里使用參數(shù)P 控制垂直局地化的幅度，當(dāng)參數(shù)P＜1時，垂直相關(guān)變寬，試驗中取為0.5。這時，局地化垂直相關(guān)結(jié)構(gòu)如圖4c 所示。比較垂直相關(guān)模型（圖4a）給出的垂直局地化結(jié)構(gòu)，可以發(fā)現(xiàn)兩者在邊界層比較類似，而在對流層中層與高層差異較大。在實際應(yīng)用垂直局地化矩陣（ Cs）時，為了避免模式層頂附近的層次存在較大噪音的干擾，不考慮模式層頂最高的5 層。

圖4 集合樣本估計的垂直局地化結(jié)構(gòu)（a. 統(tǒng)計模型，b. 集合樣本估計的流函數(shù)垂直相關(guān)結(jié)構(gòu)，c. 集合樣本估計的垂直局地化結(jié)構(gòu)）Fig. 4 Vertical localization structure estimated by ensemble samples （a. statistic model，b. stream-function vertical correlation structure estimated by ensemble samples，c. vertical localization structure estimated by ensemble samples）

進(jìn)行垂直局地化時，要對無量綱氣壓增量（ δπ）進(jìn)行局地化，由于式（5）右端涉及到 δπ在垂直方向的梯度，對 δπ進(jìn)行局地化會顯著改變梯度的計算，進(jìn)而影響到虛位溫的計算。但由于虛位溫是通過靜力平衡關(guān)系計算得到的，因而垂直局地化后的分析場仍然滿足靜力平衡關(guān)系。使用單點試驗來分析垂直局地化及局地化相關(guān)尺度對分析的影響。由圖5 可見，如果垂直局地化相關(guān)尺度較大（式（7）中經(jīng)驗系數(shù)K 取為1）時，單點觀測產(chǎn)生的氣壓分析增量分布在整個模式層次。對比圖1c 可以發(fā)現(xiàn)，由于集合預(yù)報氣壓誤差從觀測位置向地面呈現(xiàn)出逐步增大的趨勢，垂直局地化對氣壓分析增量的影響非常大。對比圖1d 可以發(fā)現(xiàn)，緯向風(fēng)場誤差大值區(qū)在垂直方向位于模式30—40 層（400—200 hPa），相應(yīng)的緯向風(fēng)場分析增量也在30—40 層，隨垂直相關(guān)尺度的影響比較小。當(dāng)垂直局地化相關(guān)尺度適當(dāng)時（K 值為7—20），分析增量被約束在觀測點附近，氣壓分析增量與風(fēng)場分析增量對應(yīng)較好。采用集合樣本估計的垂直局地化結(jié)構(gòu)對氣壓分析增量在垂直方向的限制作用更強。垂直局地化的影響表明，影響觀測信息在空間傳播的重要因素來自集合預(yù)報誤差本身，垂直局地化主要起到約束虛假遙相關(guān)的作用。這也表明當(dāng)考慮集合預(yù)報誤差在GRAPES 全球4DVar 中應(yīng)用時，集合樣本質(zhì)量及其估計的預(yù)報誤差的合理性，對集合預(yù)報誤差的應(yīng)用起到關(guān)鍵作用。

圖5 垂直局地化相關(guān)尺度對單點氣壓觀測分析增量的影響 （a. K=1，b. K=7，c.K=20，d. 集合樣本估計的局地化尺度；等值線為緯向風(fēng)分析增量，單位：m/s；色階為氣壓分析增量，單位：hPa）Fig. 5 Impacts of correlation scale of vertical localization on analysis increment when assimilating single-point pressure observations （a. K=1，b. K=7，c.K=20，d. correlation scale of vertical localization estimated by ensemble samples；contours are for zonal wind analysis increment，unit：m/s；shaded are for pressure analysis increment，unit：hPa）

進(jìn)一步考察垂直局地化對位溫分析增量的影響，結(jié)果如圖6 所示。由式（5）可知，位溫分析增量一般位于無量綱氣壓增量（ δπ）的垂直梯度最大的位置。當(dāng)相關(guān)尺度取不同數(shù)值時， δπ的垂直梯度發(fā)生了明顯變化，由靜力平衡關(guān)系導(dǎo)出的位溫分析增量同樣也發(fā)生明顯變化。當(dāng)垂直局地化相關(guān)尺度較大時位溫分析增量最大在1 K 左右，而當(dāng)垂直局地化相關(guān)尺度較小時位溫分析增量可以達(dá)到1.4 K左右。采用集合樣本估計的垂直相關(guān)模型，其對無量綱氣壓與虛位溫分析增量的影響，與垂直相關(guān)模型的影響較為一致。這表明垂直局地化結(jié)構(gòu)對分析增量的影響較小，影響較大的是垂直相關(guān)尺度，需要合理給定。

5.2 分析場平衡特征分析

4DVar 中低分辨率切線性模式積分時存在由諸多因素引起的高頻重力振蕩，當(dāng)集合預(yù)報誤差引入4DVar 后，與之相隨的水平與垂直局地化也是造成分析不平衡、激發(fā)高頻重力振蕩的原因之一。劉艷等（2019）研究表明，通過在4DVar 中引入基于數(shù)字濾波方案的初始化過程，可以很好地抑制高頻振蕩過程，使得分析場在動力上更為平衡。由式（8），Jc是目標(biāo)函數(shù)約束項，其數(shù)值反映了對重力波的控制幅度。以單點試驗為例比較4DVar 與En4DVar，分析不同局地化方案對 Jc項收斂的影響，結(jié)果如圖7 所示。4DVar 在單點觀測的目標(biāo)函數(shù)迭代下降27 步后， Jc項數(shù)值減小到0.14。而在En4DVar中， Jc項收斂過程要慢于4DVar，最終收斂數(shù)值在0.2 左右，也要明顯大于4DVar 的結(jié)果。采用不同局地化相關(guān)尺度的結(jié)果類似，不過在非平衡分析變量空間進(jìn)行局地化的結(jié)果要略好于在分析變量空間進(jìn)行局地化。這表明采用集合預(yù)報誤差后，切線性模式積分產(chǎn)生的高頻重力振蕩要強于4DVar。不過由于目標(biāo)函數(shù)中 Jc項的存在，分析場仍能保證動力平衡約束關(guān)系。

圖6 同圖5，但為虛位溫分析 （等值線為虛位溫分析增量，單位：K；色階為無量綱氣壓分析增量）Fig. 6 Same as Fig. 5 but for virtual potential temperature （Contours are for virtual potential temperature analysis increment，unit：K；shaded are for non-dimensional pressure analysis increment）

圖7 集合預(yù)報誤差引入4DVar 對目標(biāo)函數(shù)約束項 Jc的影響Fig. 7 Ensemble forecast errors introduced into 4DVar and their impact on Jc constraint term in cost-function

這里將4DVar 退化為3DVar，不考慮切線性模式積分的影響，以進(jìn)一步評估水平與垂直局地化對分析場平衡特征的影響。利用引入集合預(yù)報誤差的3DVar 同化實際常規(guī)資料與衛(wèi)星資料，并用分析場驅(qū)動GRAPES 全球模式，通過模式啟動初期地面氣壓的時間傾向來分析局地化對分析場平衡特征的影響（圖8）。由于背景場中未引入新的同化增量，由其驅(qū)動的模式軌跡平衡特征較好，這里作為對照試驗。從圖8 中可以看到，當(dāng)在分析變量空間進(jìn)行水平局地化時，模式積分初始時刻的地面氣壓傾向明顯大于背景場啟動的結(jié)果，且在5 h 積分后仍然高于背景場結(jié)果。而在非平衡分析變量空間進(jìn)行水平局地化時，平衡特征明顯改善。這進(jìn)一步說明，在非平衡分析變量空間進(jìn)行水平局地化更有優(yōu)勢。當(dāng)在水平局地化基礎(chǔ)上進(jìn)一步考慮垂直局地化時，地面氣壓傾向增加十分顯著。這表明相較于水平局地化，垂直局地化是造成分析不平衡的主要原因。對不同大小的垂直局地化相關(guān)尺度，隨局地化特征尺度變寬，不平衡程度會略有所減緩，但改善不大。這表明垂直局地化造成的不平衡與垂直局地化的相關(guān)尺度有關(guān)，但其影響較小。此外，采用集合樣本統(tǒng)計給出的垂直局地化結(jié)構(gòu)對分析平衡的改善也有限（圖略）。

圖8 背景與分析地面氣壓傾向隨預(yù)報時間的演變 （水平局地化尺度為1500 km，background 為背景，uv-nonloc 為引入集合預(yù)報誤差但不進(jìn)行局地化，uv-2dloc 為對分析變量進(jìn)行二維局地化，psi-2dloc 為對非平衡分析變量進(jìn)行二維局地化，psi-3dloc 為對非平衡分析變量進(jìn)行三維局地化）Fig. 8 Temporal evolution of background and analysis surface pressure tendency （the correlation scale of horizontal localization is 1500 km；"background" indicates background field，"uv-nonloc" indicates the introduced ensemble forecast errors without non-localization，"uv-2dloc" indicates horizontal localization on analysis variables，"psi-2dloc" is for horizontal localization on unbalanced analysis variables，"psi-3dloc" is for three-dimensional localization on unbalanced analysis variables）

6 集合預(yù)報誤差對分析的影響

6.1 集合預(yù)報誤差權(quán)重影響

氣候背景誤差與集合預(yù)報誤差所起的作用由其權(quán)重系數(shù) βc、 βe來 決定。當(dāng) βe逐步增大時，集合預(yù)報誤差逐步占主導(dǎo)。試驗利用40 個集合樣本，當(dāng)βc=1而 βe=0時，同化時間窗起始時刻的分析增量呈現(xiàn)各向同性（圖9a）。隨著 βe的值逐步增大，分析增量明顯沿著風(fēng)速較大的方向進(jìn)行延伸，而與風(fēng)速垂直的方向分析增量受到抑制（圖9c）。起始時刻分析增量呈現(xiàn)出隨集合背景誤差變化的特征，也即背景誤差與背景場的動力演變機制有關(guān)，具備隨天氣系統(tǒng)流型而變的特征。對于僅采用氣候背景誤差的分析增量，由于切線性模式動力演變機制的存在，在同化時間窗終止時刻分析增量也呈現(xiàn)出隨流型變化的特點（圖9b），而對于集合預(yù)報誤差占主導(dǎo)作用的分析增量（圖9c、d），終止時刻的分析增量隨流型變化的特點更加明顯。對比圖9c 和d，終止時刻的分析增量與起始時刻的分析增量在特征上更加一致，流依賴分析增量的中心位置隨氣流向下游傳播，這是使用集合預(yù)報誤差對4DVar 改進(jìn)的一個主要特征。

圖9 4DVar 同化時間窗內(nèi)分析增量的變化（a、b. 氣候背景誤差在 t0 與 tL 時的分析增量，c、d. 引入集合預(yù)報誤差后t0 與tL 時的分析增量；等值線為氣壓，單位：hPa，箭矢為風(fēng)）Fig. 9 Evolution of analysis increment in the 4DVar assimilation time window （a，b. analysis increments at t0 and tL with climatic background error；c，d. increments at t0 and tL after introducing ensemble forecast errors；the contours are for pressure，unit：hPa，arrows are wind speed）

6.2 計算效率影響

圖10 增加擴展控制變量與集合樣本數(shù)對計算效率的影響 （采用GRAPES 3DVar 參數(shù)配置的數(shù)值試驗；右側(cè)縱軸為集合樣本數(shù)對應(yīng)的抽樣誤差）Fig. 10 Impacts of adding extended control variables and the number of ensemble samples on computational efficiency （numerical experiment with GRAPES 3DVar configuration；the vertical axis on the right shows the sampling error corresponding to the number of ensemble samples）

7 總結(jié)與討論

GRAPES 模式是中國自主發(fā)展的非靜力格點模式，全球模式采用三/四維變分資料同化作為分析方法。文中研究在現(xiàn)有四維變分資料同化的基礎(chǔ)上，有效使用集合樣本估計預(yù)報誤差擴展變分資料同化的分析能力。圍繞GRAPES 全球4DVar 在具體實施上的特殊性，探討了在引入擴展控制變量過程中，如何在大于100 個集合樣本數(shù)條件下，增加擴展控制變量的同時能顯著降低擴展控制變量的維數(shù)增加，并減少對計算機存儲和計算的壓力，提高算法的實用性。此外，重點研究了水平局地化與垂直局地化對分析場平衡性的影響，通過單點試驗、常規(guī)與衛(wèi)星觀測資料試驗，分析了局地化對地面氣壓傾向的影響。

研究表明，可以充分利用GRAPES 全球4DVar氣候背景誤差結(jié)構(gòu)模型在水平方向采用譜濾波和在垂直方向采用經(jīng)驗正交分解的特點，利用譜濾波實現(xiàn)水平局地化，利用經(jīng)驗正交分解實現(xiàn)垂直局地化?？紤]到水平局地化所采用的相關(guān)模型更寬，對譜濾波可以采用更少的波數(shù)、更粗的高斯網(wǎng)格來降低擴展控制變量的空間維數(shù)?？紤]到需要獨立定義擴展控制變量所依托的網(wǎng)格，在實施上較為復(fù)雜，目前GRAPES 變分同化系統(tǒng)只用了較寬的水平特征長度模型，沒有再定義獨立的高斯網(wǎng)格。對垂直局地化，可以通過僅采用前8 個主導(dǎo)特征模態(tài)來降低空間的維數(shù)，此時主導(dǎo)特征模態(tài)方差解釋占比超過90%。在引入20—180 個集合樣本的情形下，對二維局地化，控制變量增加至1.1—1.8 倍，墻鐘時間增加至1.1—1.6 倍。對三維局地化，僅考慮垂直局地化降低空間維數(shù)，擴展控制變量增加至1.7—7.1 倍，墻鐘時間增加至1.4—4.1 倍。對100—200 個集合樣本，引入擴展控制變量的計算機墻鐘時間增加可控，文中所發(fā)展的引入擴展控制變量的方式具備實用性。

為了減少局地化操作對分析平衡的負(fù)面影響，發(fā)展水平局地化方案時，先將風(fēng)壓準(zhǔn)地轉(zhuǎn)平衡關(guān)系抽取出來，然后在非平衡分析變量上進(jìn)行局地化操作，最后再將風(fēng)壓準(zhǔn)地轉(zhuǎn)平衡關(guān)系加回到非平衡分析變量上。試驗結(jié)果表明，相比直接在分析變量空間進(jìn)行水平局地化，本研究發(fā)展的在非平衡分析變量空間進(jìn)行局地化的方案能更好保持分析場的平衡特征。研究還發(fā)現(xiàn)，引入垂直局地化會對分析場平衡特征產(chǎn)生較大影響，但依靠GRAPES 全球4DVar目標(biāo)函數(shù)中重力波弱約束控制，分析場仍能保證動力平衡約束關(guān)系。

氣候背景誤差與集合預(yù)報誤差所起的作用由其權(quán)重系數(shù) βc、 βe來 決定。隨著 βe的值逐步增大，單點試驗表明分析增量明顯沿著風(fēng)速較大的方向延伸，而在垂直風(fēng)速方向分析增量受到抑制。分析增量呈現(xiàn)出隨天氣形勢變化的特征。將權(quán)重系數(shù)βc、βe放在分析增量的計算中，而不是放在目標(biāo)函數(shù)的背景誤差與集合預(yù)報誤差項上，主要是改進(jìn)海森矩陣的條件數(shù)，加快變分同化系統(tǒng)收斂。

致 謝：感謝數(shù)值預(yù)報中心蘇勇、孫健、王金成博士對本研究給予的支持。