瞿安祥 麻素紅 張 進

1 中國氣象局地球系統(tǒng)數(shù)值預(yù)報中心，北京 100081 2 災(zāi)害天氣國家重點實驗室，北京 100081

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

引 言

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

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

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

1 CMA-TYM 3DVar同化系統(tǒng)

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

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

2 混合En3DVar同化方案

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

(1)

其中

x′=x-xb

β1+β2=1

(2)

(3)

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

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

(4)

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

(5)

分析增量場最優(yōu)解為：

HTR-1[yo-H(xb)]

(6)

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

(7)

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

3 混合En3DVar同化方案實現(xiàn)

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

(8)

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

(9)

x′的計算表達式為：

(10)

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

3)讀取各種類型觀測資料yo并設(shè)定相應(yīng)觀測誤差；

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

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

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

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

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

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

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

4 集合預(yù)報擾動數(shù)據(jù)來源

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

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

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

圖1 由集合預(yù)報擾動樣本統(tǒng)計的背景場狀態(tài)變量 (a)π,(b)u,(c)v,(d)q在850 hPa層上的誤差分布 (：臺風(fēng)中心位置，下同)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 由集合預(yù)報擾動樣本統(tǒng)計的背景臺風(fēng)中心氣壓與狀態(tài)變量(a)π，(b)u，(c)v，(d)q 在模式垂直層上(沿臺風(fēng)中心經(jīng)、緯向剖面)的誤差相關(guān)性分布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 臺風(fēng)中心氣壓的單點同化分析

依據(jù)前文算法步驟，本文發(fā)展實現(xiàn)了混合En-3DVar同化方案。為了驗證其合理性和正確性，選取單點臺風(fēng)中心氣壓數(shù)據(jù)進行同化測試，并將分析結(jié)果和傳統(tǒng)3DVar方案進行了對比分析。

圖3 同化臺風(fēng)中心氣壓資料后，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，但為在模式垂直層上經(jīng)向沿臺風(fēng)中心的分布Fig.4 Same as Fig.3, but for meridional distribution along typhoon center at the model vertical levels

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

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

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

(11)

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

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

7 臺風(fēng)預(yù)報試驗

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

圖5 同化常規(guī)觀測資料后，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，但為在模式垂直層上沿臺風(fēng)中心的(a,b)經(jīng)向、(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模式分別應(yīng)用3DVar系統(tǒng)和混合En3DVar系統(tǒng)產(chǎn)生的預(yù)報路徑和實際觀測路徑情況。從圖中可以看出，兩種方案對36 h內(nèi)臺風(fēng)路徑預(yù)報效果差別不大，但是對長時效的48～96 h預(yù)報，混合En3DVar方案表現(xiàn)較好，特別在“溫比亞”生命史早期幾個預(yù)報時刻，相比于3DVar方案，混合En3DVar方案產(chǎn)生的預(yù)報路徑更趨向于實際觀測，有效糾正其右偏趨勢。而在“溫比亞”生命史后期，盡管混合En3DVar方案優(yōu)勢表現(xiàn)的并不如前期那么明顯，但是實際檢驗數(shù)據(jù)顯示，其路徑預(yù)報還是好于3DVar方案。

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

圖7 CMA-TYM系統(tǒng)分別基于En3DVar和3DVar方案預(yù)報的 1306號臺風(fēng)溫比亞移動和實況路徑對比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系統(tǒng)分別基于En3DVar和3DVar方案預(yù)報的1306號臺風(fēng)溫比亞 (a)中心氣壓,(b)最大地面風(fēng)速和實況對比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 結(jié)論與討論

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

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

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

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