CMA-TYM混合En3DVar方案的設計和初步試驗*

瞿安祥 麻素紅 張 進

1 中國氣象局地球系統數值預報中心，北京 100081 2 災害天氣國家重點實驗室，北京 100081

提 要： 基于CMA-TYM的3DVar系統,發展了利用擴展控制變量引入流依賴背景誤差協方差(集合擾動成員統計表達)的混合En3DVar同化方案。測試顯示,單點臺風中心氣壓數據同化會引起風場非對稱增量形成，以及導致傳統3DVar方案認為不相關的濕度增量出現。臺風個例試驗表明,混合En3DVar方案能提取臺風內零散觀測資料信息，并按實際臺風動力特征和分布區域向四周傳播，從而影響分析場中臺風渦旋強度、結構。同時，與3DVar方案相比，混合En3DVar方案對提高臺風路徑、強度預報效果明顯。

引 言

2012年國家氣象中心將基于中尺度數值模式CMA-MESO開發的區域臺風模式系統CMA-TYM投入業務運行。近些年來檢驗表明，該系統120小時內路徑強度預報性能接近國際水平，在業務應用中具有重要參考價值。但與此同時，CMA-TYM在進一步提高其預報能力方面遇到了諸多技術難題，一個最主要的困難是如何提高初始臺風渦旋質量。雖然存在豐富的、時空密集的衛星探測數據，但受云和降水污染，這些資料并不能被分析系統有效吸收和同化。同時，受限于傳統變分同化技術局限性，部分臺風探測數據同化反而會對模式預報產生負面影響(Aberson,2008)。

進入21世紀，一種引入集合預報擾動表達流依賴背景誤差協方差的混合變分同化方案被發展起來(Hamill and Snyder，2000；Lorenc，2003；Etherton and Bishop，2004；Wang et al，2007；Wang，2010)，相比于傳統方案使用氣候統計(靜態)背景誤差參數，混合方案將集合預報擾動引入變分同化系統，顯著改善背景誤差協方差對實時天氣系統描述能力，因而對觀測數據信息的提取和傳播，更符合實際天氣系統分布特征?？茖W研究(Wang et al，2008a;2008b；Buehner et al，2010a;2010b； Zhang and Zhang，2012;Zhang et al,2013)和業務應用(Clayton et al,2013；Bonavita et al,2012;Wang,2011;Wang et al,2013；Kleist and Ide，2015)都表明，混合變分同化方案會明顯提高分析場質量及模式預報水平。隨著技術成熟，混合變分同化方案也開始應用于臺風模擬研究(Wang,2011；Hamill et al，2011；Li et al,2012)取得不錯預報效果。

因此，在CMA-TYM 3DVar基礎上，本文設計發展了通過三維Alpha控制變量引入集合預報擾動(隱式流依賴背景誤差協方差表達)的混合En3DVar同化方案。然后利用單點臺風中心氣壓數據進行測試，驗證En3DVar方案對觀測信息的提取和傳播，是否與集合預報擾動統計的背景誤差協方差特征相符合，以及形成的分析增量是否與實際臺風環流分布及動力學屬性相匹配。最后基于新建立CMA-TYM En3DVar系統進行臺風個例試驗，來分析其改進臺風初始結構方面的能力，以及提高模式路徑和強度預報方面的效果。

1 CMA-TYM 3DVar同化系統

CMA-TYM衍生于中國氣象局自主研發的新一代全球與區域同化預報系統GRAPES(Global/Regional Assimilation and Prediction System)，是一個全可壓非靜力有限差分格點模式(薛紀善等，2008；陳德輝等，2012)。水平方向為經緯網格點球面坐標，垂直方向為高度地形追隨坐標。主要特點包括：模式預報變量在水平方向采用Arakawa-C 跳點格式，垂直方向采用非均勻Charney-Philips跳層格式，時間上采用半隱式-半拉格朗日差分方案。CMA-TYM是一個覆蓋西北太平洋、南海和部分東亞大陸的區域臺風數值模式系統(張進等，2017；麻素紅等，2018;麻素紅，2019；麻素紅和陳德輝，2018)。水平分辨率為0.12°，垂直方向為50層，模式層頂約為33 000 m，參考大氣采用水平平均廓線。模式物理過程包括(郭云云等，2015；聶皓浩等，2016；鄭曉輝等，2016)：包含水汽、雨、雪、云水、云冰、霰的WSM6顯式微物理方案、YSU邊界層方案、Monin-Obukhov近地面層方案、Goddard短波輻射、RRTM長波輻射，陸面過程采用SLAB熱量擴散方案。另外，考慮0.12°分辨率不能完全描述臺風對流特征，Meso-SAS積云對流參數化方案(Pan et al,2014)也考慮在內。

針對CMA-MESO發展的CMA-TYM 3Dvar系統(馬旭林等，2009；張華等，2004；莊世宇等，2005；莊照榮等，2006)，其變量水平、垂直分布與模式完全一致。雖然CMA-TYM為非靜力模式，預報變量包括水平和垂直風、溫度、比濕、無量綱氣壓(Exner函數)、云水等。但考慮到觀測資料屬性和變分同化可操作性，3DVar選取水平風(u,v)、位溫θ、比濕q、無量綱氣壓π作為同化分析變量，其中質量場變量：無量綱氣壓π或位溫θ只能二選一(本文選取了無量綱氣壓π)。

2 混合En3DVar同化方案

以一種簡單形式，用權重平均(集合預報擾動統計的和氣候統計的背景誤差協方差)后的背景誤差協方差替換3DVar目標函數中相應項，就可以獲得混合En3DVar同化方案直接表達公式：

(1)

其中

x′=x-xb

β1+β2=1

(2)

(3)

滿足梯度為零目標函數極小化條件：

HTR-1[H(x′+xb)-yo]=0

(4)

Hx′=H(x′+xb)-H(xb)

(5)

分析增量場最優解為：

HTR-1[yo-H(xb)]

(6)

通過引入擴展控制變量α=(α1,…,αm)T，Lorenc(2003)以一種巧妙方式獲得了混合En3DVar同化方案間接表達式：

(7)

可以看出，只要選擇合適局地化矩陣模擬方法，基于擴展控制變量技術的混合En3DVar方案很容易在傳統3DVar方案上融入實現(具體實現細節見下節)。

3 混合En3DVar同化方案實現

R-1[H(Uυ+xb)-yo]

(8)

接下來為消除分析狀態變量(π,u,v,q)之間相關，通過物理算子Up將其變換為互不相關分析變量：流函數ψ、勢函數χu(與ψ不平衡部分)、πu(與ψ不平衡部分)、比濕q。然后假設ψ，χu，πu，q自身空間相關結構在水平和垂直方向上具有可分離屬性，利用垂直模態映射算子Uv與水平濾波算子Uh的Kronecker積來共同模擬實現。經過上述系列變換，目標函數及其梯度計算就轉為對新控制變量υ的迭代求解，分析增量x′的計算表達式為：x′=Uυ=UpUvUhυ。

(9)

x′的計算表達式為：

(10)

目標函數及其梯度的計算就轉換為對新控制變量(υ,w)的迭代求解。

3)讀取各種類型觀測資料yo并設定相應觀測誤差；

4)計算新息向量：d=yo-H(xb)；

5)設置分析控制變量(υ,w)迭代初值，其中w由m個向量場w1,w2,…,wm組成；

6)計算基于控制變量(υ,w)表達的目標函數值；

7)計算目標函數關于分析控制變量(υ,w)梯度值；

8)基于獲得的梯度值，計算梯度下降方向和最優步長，并更新控制變量(υ,w)迭代值。后重復計算目標函數和目標函數梯度值(步驟6和步驟7)，直至達到收斂標準；

9)計算最終分析值：xa=xb+x′。

計算資源方面，混合En3DVar方案利用的是三維α控制變量，在傳統3DVar方案分析狀態變量總維數為n前提下，會額外增加n×m維分析狀態變量空間，且隨著集合預報成員數目m增加，空間還在增大，這對計算資源硬件需求是一個比較大的挑戰。不過，在采用并行計算環境下，這個需求會隨著并行核增加而變得相應容易滿足。

4 集合預報擾動數據來源

本試驗所需集合預報擾動場數據，源自96個成員的全球(譜模式)集合預報系統T574-EnKF(1)T574-EnKF集合系統是NCEP基于采用EnSRF技術建立的全球集合預報-同化循環滾動業務系統。(Whitaker et al,2008；Wang et al,2013；Kleist and Ide，2015)的，通過CMA-TYM積分獲得。具體方案為，首先基于T574-EnKF系統輸出的分析場集合，利用降尺度技術獲得CMA-TYM初值集合。然后，考慮En3DVar分析更多關注臺風內中小尺度信息提取和傳播，而降尺度形成的初值僅包含全球模式可分辨的(較粗分辨率)天氣尺度信息，因此需通過CMA-TYM模式短時間積分，來獲得具有中小尺度信息的集合預報場(Pu et al,2016;Lu et al,2017)。

圖1顯示的是1306號臺風溫比亞(Rumbia)形成初期，基于96個集合預報擾動場統計表達的：背景場狀態變量π,u,v,q在850 hPa層上誤差(集合預報離散度)分布信息。從圖中可以看出，各個狀態變量誤差水平分布范圍，具有明顯臺風環流特征：無量綱氣壓 (圖1a)和風場u(圖1b)、v(圖1c)都顯示越靠近臺風中心，誤差值越大，并且呈現閉合環流形狀。很明顯，這些狀態變量誤差分布與模式對臺風系統描述能力吻合，臺風內部模擬準確程度要遠低于外圍環流準確程度。同時，圖1d顯示的比濕q也符合臺風降雨云系分布的螺旋形狀(螺旋雨帶降水模擬也是常見臺風模擬誤差來源)。

從圖2中可以看出，在垂直方向上，臺風中心氣壓誤差與中低層無量綱氣壓π誤差存在正相關，與高層存在負相關，這與臺風環流中低層氣旋輻合、高層反氣旋輻散結構相吻合。同時臺風中心氣壓誤差與風速u,v分量誤差在不同象限區域顯示著有正有負的相關，并且正負值所在區域與臺風風壓動力學平衡關系相匹配。

圖1 由集合預報擾動樣本統計的背景場狀態變量 (a)π,(b)u,(c)v,(d)q在850 hPa層上的誤差分布 (：臺風中心位置，下同)Fig.1 Background error distribution of (a) π, (b) u, (c) v, (d) q at 850 hPa level calculated by the ensemble forecast perturbation (：center of typhoon, same as below)

圖2 由集合預報擾動樣本統計的背景臺風中心氣壓與狀態變量(a)π，(b)u，(c)v，(d)q 在模式垂直層上(沿臺風中心經、緯向剖面)的誤差相關性分布Fig.2 Distribution of model vertical levels (meridional and zonal profile along typhoon center) of background error correlation between typhoon central pressure and (a) π, (b) u, (c) v, (d) q calculated by the ensemble forecast perturbation

5 臺風中心氣壓的單點同化分析

依據前文算法步驟，本文發展實現了混合En-3DVar同化方案。為了驗證其合理性和正確性，選取單點臺風中心氣壓數據進行同化測試，并將分析結果和傳統3DVar方案進行了對比分析。

圖3 同化臺風中心氣壓資料后，3DVar方案形成的(a)π,(c)u增量與En3DVar方案形成的 (b)π,(d)u增量在地面層上的水平分布Fig.3 Surface increment of (a) π and (c) u analyzed by 3DVar and (b) π and (d) u analyzed by En3DVar with typhoon central sea level pressure assimilation

圖4 同圖3，但為在模式垂直層上經向沿臺風中心的分布Fig.4 Same as Fig.3, but for meridional distribution along typhoon center at the model vertical levels

6 實際觀測資料同化效果分析

基于新建立的CMA-TYM En3DVar系統，本文進行了實際觀測資料同化試驗。試驗模式區域范圍為5°～49°N、93°～156°E，覆蓋部分大陸和太平洋、南海的東南亞地區，水平格點數為531×451。試驗時間是2013年6月28日12：00 UTC(1306號臺風溫比亞生成初期)。試驗冷啟動數據來自T1148-GSI(2)為本文試驗提供冷啟動數據來源的NCEP全球EnKF集合預報系統和全球確定性預報系統，是一套采用高低雙分辨率運行的雙向耦合系統，出于業務計算資源和運行時效考慮，EnKF系統采用較低分辨率的T574模式運行，而確定性預報采用高分辨率的T1534模式運行。本文出于試驗資源的考慮，確定性預報選擇了T1148模式配置運行。全球數值預報模式系統(提供初始場和側邊界)和T574-EnKF全球集合預報系統(提供集合預報場)。同化所需觀測資料來自實時數值預報業務常規觀測資料集，包括探空、地面站、洋面浮標和船舶資料、飛機報和云導風、洋面風等。除此之外，由預報員實時分析的臺風中心海平面氣壓也作為常規資料進入了分析同化系統。

γ(k1,k2)=exp(-d2/L2)

(11)

式中：d為模式層k1和k2的幾何距離，L為垂直局地化特征長度。為獲得最優取值，本文對100～1 000 km 范圍內的水平局地化特征長度進行了測試，最終發現在260 km情況下會顯示非常合理的分析增量。L也采取了類似的測試確定取值為36。

圖5、圖6顯示的是En3DVar系統同化觀測資料后，狀態變量π,u,v,q在水平(地面層)和垂直方向上(沿臺風中心經緯向剖面)的增量分布。從π來看，水平方向上(圖5a)在臺風中心周圍形成一個氣壓增量負值區；垂直方向上(圖6a)，氣壓增量呈現中低層有效加深，高層反氣旋弱增強的分布形式，匹配臺風低層輻合高層輻散的質量場屬性。從u，v分析增量來看，圖5b和5c顯示在水平方向上，伴隨著氣壓場加深，圍繞臺風中心形成明顯氣旋性環流。實際數據顯示增量分布并不均勻，東部、北部區域更大一些，呈現明顯非對稱加深結構。圖6b和6c顯示在垂直方向上，u，v分析增量緊密圍繞臺風中心，且遍布整個中低層，同樣非對稱特征也非常明顯。從濕度q分析增量來看,水平方向上(圖5d)其分布區域與臺風環流螺旋狀云系特征緊密匹配,垂直方向上(圖6d)，在臺風中心附近產生的分析增量遍布眼區及周圍環流區域(看似“雜亂”的增量分布實際和臺風中小尺度對流區域十分吻合)。

7 臺風預報試驗

為評估混合En3DVar方案對臺風路徑強度預報影響，基于CMA-TYM模式系統對1306號臺風溫比亞進行了預報試驗。試驗分兩種方案進行，一種為3DVar方案，另一種為En3DVar方案，即分別應用3DVar、En3DVar形成的分析場進行模式積分預報。運行方案為，在“溫比亞”生命史期間，每天選取00 UTC和12 UTC時刻進行同化分析和預報試驗。預報所需數據同樣來自全球模式系統T1148-GSI和全球集合預報系統T574-EnKF。

圖5 同化常規觀測資料后，En3DVar方案形成的(a)π，(b)u，(c)v，(d)q增量在地面層上的水平分布Fig.5 Surface increment of (a) π, (b) u, (c) v, (d) q analyzed by En3DVar with conventional data assimilation

圖6 同圖5，但為在模式垂直層上沿臺風中心的(a,b)經向、(c,d)緯向分布Fig.6 Same as Fig.5, but for (a, b) meridional and (c, d) zonal distribution along typhoon center at the model levels

圖7顯示在“溫比亞”生命史期間，CMA-TYM模式分別應用3DVar系統和混合En3DVar系統產生的預報路徑和實際觀測路徑情況。從圖中可以看出，兩種方案對36 h內臺風路徑預報效果差別不大，但是對長時效的48～96 h預報，混合En3DVar方案表現較好，特別在“溫比亞”生命史早期幾個預報時刻，相比于3DVar方案，混合En3DVar方案產生的預報路徑更趨向于實際觀測，有效糾正其右偏趨勢。而在“溫比亞”生命史后期，盡管混合En3DVar方案優勢表現的并不如前期那么明顯，但是實際檢驗數據顯示，其路徑預報還是好于3DVar方案。

從強度預報檢驗來看(圖8)，無論是中心氣壓還是近地面最大風速，混合En3DVar方案產生的預報結果都好于3DVar方案。特別在“溫比亞”生成初期，混合En3DVar方案就提前72 h準確預報出其在近?？焖僭鰪姳l過程，并且這種表現一直延續貫穿后期所有預報時刻，直至臺風登陸消亡。反觀3DVar方案，盡管也預報描述出臺風后期增強過程，但是這個趨勢緩慢且不明顯，而且在數據量級(氣壓與風速大小)上也和實際觀測相差甚遠。

圖7 CMA-TYM系統分別基于En3DVar和3DVar方案預報的 1306號臺風溫比亞移動和實況路徑對比Fig.7 Tracks predicted by CMA-TYM model with En3DVar and 3DVar scheme for No.1306 Typhoon Rumbia compared with the observed track

圖8 2013年6月28日至7月1日CMA-TYM系統分別基于En3DVar和3DVar方案預報的1306號臺風溫比亞 (a)中心氣壓,(b)最大地面風速和實況對比Fig.8 Central pressure (a) and maximum surface wind (b) predicted by CMA-TYM model with En3DVar and 3DVar scheme for No.1306 Typhoon Rumbia compared with the observed data from 28 June to 1 July 2013

8 結論與討論

基于CMA-TYM 3DVar同化系統，本文設計發展了通過三維Alpha控制變量引入集合預報擾動(隱式流依賴背景誤差協方差表達)的混合En3DVar同化方案。在匹配現有3DVar方案算法流程的基礎上，具有計算高效、易于實現的特點。

單點臺風中心氣壓數據同化測試表明，相比于3DVar方案，混合En3DVar方案形成的氣壓場分析增量明顯分布于臺風環流區域內，風場分析增量結構上更加匹配臺風動力學特征：數值量級上更大且非對稱結構明顯。更重要的是，氣壓分析增量會導致濕度場(比濕)增量的出現(傳統3DVar方案假設兩者是不相關的)。由此可見，在集合預報擾動統計表達的背景誤差協方差能準確表征背景臺風誤差條件下,臺風內部資料的同化就能夠獲得與臺風環流范圍相匹配的分析增量,并且這些分析增量內部之間的相關平衡符合實時臺風動力學特征。

應用常規觀測資料的同化試驗顯示,混合En3DVar方案能有效提取臺風區域內零散觀測信息，并且這部分信息會按實時臺風特征向周圍傳播出去，從而影響分析場臺風結構、強度變化。這使得通過有限觀測資料信息同化改進初始臺風渦旋質量成為可能。同時，臺風個例預報試驗表明，無論是路徑還是強度，混合En3DVar產生的預報結果都要明顯好于3DVar預報結果。

需要指出的是，混合變分同化方案形成的分析增量實質是多個集合成員預報擾動的線性組合(被限制在集合成員擾動所張的子空間)，其同化效果與集合預報擾動統計表達的背景誤差協方差的準確程度密切相關。只要有限集合預報擾動成員能夠抓住實時天氣系統在模式格點尺度體現的主要誤差特征，那么基于其統計所得的誤差協方差矩陣就具有實時天氣系統的流依賴型表達能力。