馮業(yè)榮 薛紀(jì)善 李夢婕 戴光豐

中國氣象局廣州熱帶海洋氣象研究所/廣東省區(qū)域數(shù)值天氣預(yù)報(bào)重點(diǎn)實(shí)驗(yàn)室，廣州，510641

1 引言

初值是數(shù)值天氣預(yù)報(bào)的關(guān)鍵科學(xué)技術(shù)問題。已有研究和業(yè)務(wù)實(shí)踐表明，在變分資料同化方法中由于考慮不同時(shí)刻觀測資料的影響并受預(yù)報(bào)模式的約束，四維變分資料同化（4DVar）相比于三維變分資料同化（3DVar）可為數(shù)值模式提供更好的初值，進(jìn)而改善模式預(yù)報(bào)效果（Lorenc，et al，2005；Rawlins，et al，2007）。故此，歐洲中期天氣預(yù)報(bào)中心（ECMWF）等國際上的先進(jìn)數(shù)值預(yù)報(bào)中心從20世紀(jì)末開始相繼采用4DVar作為基本的業(yè)務(wù)同化系統(tǒng)（Courtier，et al，1994；Gauthier，et al，2007；Rawlins，et al，2007）。尤其是ECMWF，1996年3DVar投入使用（Courtier，et al，1998），但很快便于1997年由4DVar取代（Rabier，et al，2000）。4DVar和衛(wèi)星觀測資料的結(jié)合，提高了其全球模式的 預(yù) 報(bào) 質(zhì) 量（McNally，et al，2006；Bauer，et al，2010）。尤其對南半球，盡管直接觀測資料不足，但ECMWF全球模式的可用預(yù)報(bào)時(shí)效在2000年前后有較大提高，4DVar與衛(wèi)星觀測資料同化應(yīng)用功不可沒（Bauer，et al，2015）。經(jīng)過多年努力，中國氣象局的全球四維變分同化系統(tǒng)（GRAPES-4DVar，現(xiàn)為CMA-GFS 4DVar）也于2018年7月投入業(yè)務(wù)運(yùn)行（Zhang，et al，2019），使中國自主數(shù)值預(yù)報(bào)發(fā)展的步伐又邁出了關(guān)鍵的一步。

迄今為止，4DVar較多的是與全球模式結(jié)合。對于區(qū)域模式而言，初值或直接取自全球模式降尺度分析場，或通過混合使用區(qū)域3DVar與其他方法（Jones，et al，1997；Dixon，et al，2009；Bloom，et al，1996）得到，且同化所用的資料以單時(shí)刻觀測為主。為了科學(xué)利用本地區(qū)多時(shí)刻稠密觀測資料，改善區(qū)域模式的初值質(zhì)量，提高模式預(yù)報(bào)準(zhǔn)確度，建立區(qū)域4DVar同化系統(tǒng)是一個(gè)較好的選擇。中國由于區(qū)域GRAPES模式已被廣泛用于業(yè)務(wù)，開發(fā)與GRAPES模式相配套的區(qū)域4DVar系統(tǒng)更具有現(xiàn)實(shí)意義。

開發(fā)4DVar最復(fù)雜與最繁重的工作在于開發(fā)切線性模式及其伴隨模式。切線性模式通常要求擾動(dòng)趨于比較小時(shí)與非線性模式的偏差保持一致，它的程序可以根據(jù)非線性模式代碼按一定規(guī)則直接編程，這種做法使切線性模式和伴隨模式的建立成為近似機(jī)械的操作，盡管工作量大，難度一般不大，但后續(xù)系統(tǒng)維護(hù)和升級比較困難。而且由于非線性模式中存在各種復(fù)雜物理開關(guān)，連續(xù)的切線性模式其實(shí)并不存在，對物理過程處理不好會(huì)導(dǎo)致極小化迭代失??；另一種做法是不依賴非線性模式代碼，而是直接建立擾動(dòng)模式方程組并通過簡化的線性物理過程使擾動(dòng)模式保持較好的線性特征，以利于極小化計(jì)算平穩(wěn)進(jìn)行。這種所謂擾動(dòng)模式并不追求在擾動(dòng)為0的極限情況下無限逼近切線性模式，僅要求一定隨機(jī)分布的擾動(dòng)期望值的演變與擾動(dòng)在切線性模式中的演變一致。為此需要從非線性模式控制方程組出發(fā)，推導(dǎo)一套線性擾動(dòng)方程組（包括線性物理過程），進(jìn)而開發(fā)擾動(dòng)預(yù)報(bào)模式及其伴隨模式，如英國氣象局的4DVar系統(tǒng)便是基于擾動(dòng)模式建立的（Rawlins，et al，2007）。

馮業(yè)榮等（2020）在GRAPES模式原始方程組基礎(chǔ)上建立了線性擾動(dòng)預(yù)報(bào)模式（GRAPES_PF），試驗(yàn)表明該擾動(dòng)模式是GRAPES切線性模式的一個(gè)合理近似。本研究是上述工作的繼續(xù)和深入，即通過GRAPES_PF搭建適應(yīng)GRAPES區(qū)域模式的4DVar系統(tǒng)框架，為最終建立GRAPES區(qū)域四維變分同化業(yè)務(wù)系統(tǒng)奠定基礎(chǔ)。

2 擾動(dòng)4DVar框架設(shè)計(jì)

2.1 控制變量的設(shè)置完善

四維變分同化問題最終歸結(jié)為求以下目標(biāo)函數(shù)的極小值

式中，w為控制變量，定義在模式空間；Ri為觀測余項(xiàng)，定義在觀測空間；O為觀測誤差協(xié)方差矩陣；上標(biāo)T代表對矩陣或向量進(jìn)行轉(zhuǎn)置，上標(biāo)?1代表對矩陣求逆，下標(biāo)i表示觀測對應(yīng)的時(shí)間序號。觀測余項(xiàng)Ri=di+HM0→iPUw。 其中di=H[xb(ti)] ?yoi稱為新息向量，用于衡量背景場與觀測的偏離程度。xb(ti)是 按非線性模式演變的背景場，yoi是觀測量，函數(shù)H[xb(ti)]稱為由背景場計(jì)算得到的觀測相當(dāng)量；H是H(·)的 線性觀測算子；M0→i代 表從分析時(shí)刻t0到觀測時(shí)刻ti的切線性模式積分；P和U分別稱為物理變換矩陣和預(yù)條件變換矩陣，是為了改善極小化收斂條件而引入的。變換P先將風(fēng)速變量（u,v）變換為流函數(shù)ψ和勢函數(shù)χ，再將勢函數(shù)χ分解為平衡部分χb和非平衡部分χu；位勢高度也分解為平衡（φb） 和非平衡（φu） 兩部分，φb與ψ通過線性平衡方程建立聯(lián)系（薛紀(jì)善等，2008）；變換U將背景誤差協(xié)方差矩陣分解為水平相關(guān)矩陣與垂直相關(guān)矩陣的乘積，計(jì)算時(shí)針對垂直相關(guān)矩陣進(jìn)行經(jīng)驗(yàn)正交函數(shù)（EOF）分解，水平相關(guān)矩陣使用遞歸濾波來逼近。經(jīng)過以上系列變換處理，背景誤差協(xié)方差矩陣的相關(guān)性和規(guī)模得到降低，便于極小化問題的求解。

原GRAPES 3DVar的平衡方程是定義在氣壓坐標(biāo)上的，先在等壓面上求解平衡方程得到φb，然后通過靜力平衡求平衡氣壓Πb。由于GRAPES模式采用了地形追隨垂直坐標(biāo)（），等壓面上的求解結(jié)果最后還需要插值到模式面上，容易增加額外的計(jì)算誤差，因此本研究在模式面上建立平衡方程直接求解Πb。 經(jīng)過推導(dǎo)，坐標(biāo)的線性平衡方程為

u，v用流函數(shù)ψ和勢函數(shù)χ表示為

于是水平渦度

其中拉普拉斯算子為

進(jìn)一步假設(shè)平衡運(yùn)動(dòng)為水平無輻散的渦旋運(yùn)動(dòng)，于是式（2）可寫成

式（3）已用無量綱氣壓Π代替原來使用的位勢高度（φ=gz），但其中的變量均為全量，變分同化使用的是擾動(dòng)量，因此經(jīng)線性化后，得到擾動(dòng)形式的線性平衡方程（泊松方程）

式中，上標(biāo)“′”代表擾動(dòng)量（或稱增量）。

原GRAPES 3DVar只使用4個(gè)控制變量：流函數(shù)ψ，非平衡速度勢χu，質(zhì)量變量m（代表非平衡氣壓Πu或非平衡位溫θu，二擇一）和濕度變量Hum（比濕（q）或相對濕度（RH），二擇一）。本研究控制變量增加了非平衡溫度Tu，即w由5個(gè)分量組成：流函數(shù)ψ，非平衡速度勢χu、 非平衡氣壓Πu、非平衡溫度Tu和相對濕度（RH）。增加溫度作為第5個(gè)控制變量對于擾動(dòng)四維變分同化系統(tǒng)是方便的，不僅使得分析變量與擾動(dòng)模式變量直接對應(yīng)，同時(shí)也為今后開展非靜力分析做準(zhǔn)備，否則所需要的位溫（或氣壓）必須通過靜力平衡關(guān)系由對應(yīng)的氣壓（或位溫）求得，由于GRAPES模式（包括其擾動(dòng)模式）是非靜力的，這樣無疑增加了誤差。

2.2 擾動(dòng)模式和伴隨模式

切線性模式采用擾動(dòng)預(yù)報(bào)模式GRAPES_PF（馮業(yè)榮等，2020）代替。擾動(dòng)模式的積分過程可以表示為

式中，Mi?1→i表示擾動(dòng)模式從ti?1時(shí)刻到ti時(shí)刻的積分。

目標(biāo)函數(shù)極小化過程是一個(gè)迭代過程，采用LBFGS擬牛頓法(Liu，et al，1989)求解。初始控制矢量一般取w=0，在模式反復(fù)約束下通過不斷引入多時(shí)刻觀測資料，得到使目標(biāo)函數(shù)達(dá)到最小值的w。每次迭代需要搜索合適的下降方向，最速下降方向是目標(biāo)函數(shù)梯度的反方向，為此需要反復(fù)計(jì)算目標(biāo)函數(shù)的梯度。

與多數(shù)4DVar相似，采用“增量分析”的策略，因此計(jì)算涉及兩個(gè)循環(huán)（外循環(huán)和內(nèi)循環(huán)）。外循環(huán)過程中高分辨非線性模式從分析時(shí)刻開始連續(xù)積分到同化窗的每一個(gè)觀測時(shí)刻，所得到的預(yù)報(bào)量用于計(jì)算新息向量。內(nèi)循環(huán)在極小化過程內(nèi)部發(fā)生，在搜索目標(biāo)函數(shù)的下降點(diǎn)過程中需要反復(fù)計(jì)算梯度 ?J，因此需要多次反復(fù)積分?jǐn)_動(dòng)模式和伴隨模式，因計(jì)算量巨大，故采用低分辨率的擾動(dòng)模式及其伴隨模式。本研究中外循環(huán)分辨率為12 km，內(nèi)循環(huán)分辨率取36 km。

擾動(dòng)模式的計(jì)算最終體現(xiàn)為通過廣義共軛余差法（GCR）迭代求解擾動(dòng)亥姆霍茲方程（馮業(yè)榮等，2020）。由于GCR迭代是否終止由收斂判據(jù)決定，事先并不能確定迭代次數(shù)，故無法根據(jù)一般伴隨規(guī)則直接編寫伴隨程序。實(shí)際上擾動(dòng)亥姆霍茲方程的伴隨方程也是亥姆霍茲方程，其系數(shù)矩陣不過是原來擾動(dòng)亥姆霍茲方程系數(shù)矩陣的轉(zhuǎn)置。因此，文中仍利用GCR方法迭代求解亥姆霍茲伴隨方程，只需將系數(shù)矩陣“行、列”位置互換即可，當(dāng)然系數(shù)矩陣的“行、列”在亥姆霍茲求解程序中并不顯現(xiàn)，編程者須將下標(biāo)變量的對應(yīng)關(guān)系弄清楚。

編寫伴隨模式需要開展伴隨檢驗(yàn)，若設(shè)擾動(dòng)模式和伴隨模式的初值分別為 dx0和dxi?，則擾動(dòng)模式和伴隨模式構(gòu)成的內(nèi)積應(yīng)滿足檢驗(yàn)判據(jù)

3 構(gòu)造同化試驗(yàn)數(shù)據(jù)

為了開展同化試驗(yàn)，利用GRAPES模式模擬天氣“實(shí)況”，借此形成幾組人造 “探空”觀測數(shù)據(jù)。所選個(gè)例是2015年南海臺(tái)風(fēng)“彩虹”?！安屎纭逼鹪从诜坡少e呂宋島東部，2015年10月2日移入中國南海，隨后快速加強(qiáng)并向西北方向移動(dòng)，4日06時(shí)10分（世界協(xié)調(diào)時(shí)，下同）以強(qiáng)臺(tái)風(fēng)級別在湛江市登陸?！安屎纭庇绊懫陂g，華南區(qū)域氣象中心的GRAPES臺(tái)風(fēng)業(yè)務(wù)模式（現(xiàn)為CMA-TRAMS）曾于2015年10月2日12時(shí)準(zhǔn)確預(yù)報(bào)了“彩虹”的移動(dòng)和登陸，預(yù)報(bào)路徑幾乎與實(shí)況重合（圖1）。鑒于當(dāng)前擾動(dòng)4DVar動(dòng)力框架尚未包括物理過程，因此，文中沒有選擇包含全套物理過程且預(yù)報(bào)效果較好的業(yè)務(wù)模式作為4DVar的非線性模式，僅對GRAPES非線性模式進(jìn)行干模式積分，即關(guān)閉了所有物理選項(xiàng)（輻射、微物理、邊界層和積云對流等）。模擬范圍（12°—24°N，105°—123°E），水平分辨率0.12°，模式垂直67層，模式頂高30 km。取ECMWF全球預(yù)報(bào)模式IFS資料作初邊值，以2日12時(shí)為模式初始時(shí)刻，積分42 h。圖1可見，與濕模式（業(yè)務(wù)模式）相比，干模式預(yù)報(bào)的登陸位置偏東，登陸時(shí)間早6 h，但仍可報(bào)出“彩虹”向西北方向移動(dòng)的趨勢，因此將此干模式積分結(jié)果作為理想 “真實(shí)”大氣，既可對應(yīng)現(xiàn)有擾動(dòng)模式框架，又不失去“彩虹”運(yùn)動(dòng)的基本特征。需要指出，圖1中ECMWF分析場的臺(tái)風(fēng)中心位置（A）比實(shí)際位置（O）偏西32.6 km，文中做干模式積分時(shí)未對初始場臺(tái)風(fēng)中心進(jìn)行重定位。

圖1 臺(tái)風(fēng)“彩虹”的觀測和預(yù)報(bào)路徑：實(shí)況 （虛線）、業(yè)務(wù)模式 （點(diǎn)線） 和干模式 （實(shí)線），干模式臺(tái)風(fēng)初始位置 （A）位于臺(tái)風(fēng)實(shí)際位置 （O） 偏西32.6 kmFig.1 The observed and forecasted tracks of typhoon Mujigae：best track （dashed），simulations of operational model （dotted） and dry model （solid），symbol A indicates Mujigae's initial position in the dry model，which is 32.6 km west of the observed position （O）

在中國南海及其周邊放置25個(gè)人造“探空”觀測站，間距100—300 km不等（圖2）。為更好描述臺(tái)風(fēng)中心特征，有意在“彩虹”中心附近使探空站分布密一些（間距80—100 km）。人為構(gòu)造了間隔3 h的3組“探空”觀測數(shù)據(jù)，即將前面干模式的初始場（ECMWF的2日12時(shí)分析場）和預(yù)報(bào)場（15 和18時(shí)）插值到每個(gè)“探空”站的基本等壓面上，形成人造“探空”觀測數(shù)據(jù)。要素包括氣壓（p）、緯向風(fēng)（u）、經(jīng)向風(fēng)（v）、位勢高度（gz）、溫度（T）和相對濕度（RH），垂直層次包括地面、1000、925、850、700、500、400、300、250、200、150、100、80、70、60、50、40、30和20 hPa。

圖2 25個(gè)人造“探空”站點(diǎn)分布Fig.2 The layout of the 25 pseudo radiosonde stations

為評估4DVar效果，人為制造一個(gè)有誤差的背景場。利用Kurihara等（1993）提出的方法將ECMWF的2日12時(shí)分析場（圖3a）進(jìn)行渦旋分離和重定位，即以不同波長在水平面上連續(xù)多次濾波，得到一個(gè)平滑場，然后通過柱狀濾波器，從“真實(shí)”場分離出臺(tái)風(fēng)渦旋，將分離出來的臺(tái)風(fēng)渦旋向東北方向移位約230 km，同時(shí)將臺(tái)風(fēng)渦旋的u和v風(fēng)場乘以一個(gè)隨高度逐漸減小的系數(shù)，即對低層的風(fēng)速也做一定的改變，使渦旋環(huán)流減弱，最后得到與“真實(shí)”場有較大差異的初始場（圖3b）?；谶@個(gè)“差”的初始場重新積分干模式36 h（對照試驗(yàn)CTL），得到“差”的背景場。圖3d是由對照試驗(yàn)得到的18時(shí)的背景場，同樣與同時(shí)刻的“真實(shí)”場（圖3c）差異明顯。利用上述試驗(yàn)數(shù)據(jù)對所開發(fā)的擾動(dòng)四維變分同化系統(tǒng)進(jìn)行測試。

圖3 2015年10月2日12 （a、b） 和18 （c、d） 時(shí)850 hPa位勢高度場 （黑色線條，單位：dagpm）、風(fēng)矢量和風(fēng)速（色階，單位：m/s）（a、c.“真實(shí)”大氣，b、d.人為改變的背景場）Fig.3 Initial fields of 12：00 （a，b） and 18：00 （c，d） UTC 2 October 2015 at 850 hPa geopotential heights （black contour，unit：dagpm）, wind vectors and wind speed （shaded，unit：m/s）（a，c.idealized "truth" fields；b，d.artificially changed fields）

4 4DVar與3DVar比較試驗(yàn)

4.1 目標(biāo)函數(shù)收斂特征

四維變分同化中目標(biāo)函數(shù)J的極小化解與高斯分布條件下背景場與觀測量的聯(lián)合概率密度函數(shù)的極大似然估計(jì)等價(jià)，因此J值衡量由分析場計(jì)算的模式相空間軌跡相對于初猜場（背景場）預(yù)報(bào)軌跡和觀測量變化軌跡的距離，其極小值反映了分析場和基于分析場的模式預(yù)報(bào)同時(shí)靠近背景場和觀測值的最優(yōu)狀態(tài)。

文中將三維變分同化與四維變分同化設(shè)計(jì)成一體化系統(tǒng)，因而可以很方便地進(jìn)行3DVar或4DVar試驗(yàn)。利用這一系統(tǒng)做了4個(gè)同化試驗(yàn)：2個(gè)三維變分（3DVar_1和3DVar_2）、2個(gè)四維變分（4DVar_1和4DVar_2）。3DVar_1和3DVar_2的分析時(shí)刻分別為10月2日12 時(shí)和18時(shí)。4DVar_1和4DVar_2的分析時(shí)刻均為10月2日12時(shí)，但4DVar_1使用2個(gè)時(shí)次（12和15時(shí)）的觀測，4DVar_2使用3個(gè)時(shí)次（12、15和18時(shí)）的觀測。盡管4DVar_2分析時(shí)刻是12時(shí)，但也用到了18時(shí)的觀測，因此模式有效預(yù)報(bào)時(shí)效是18時(shí)以后，故增加分析時(shí)刻為10月2日18時(shí)的一次三維變分同化試驗(yàn)（3DVar_2），以便公平評估觀測資料對于3DVar和4DVar效果的影響。

三維變分極小化迭代的初始值取w=0，由于此時(shí)目標(biāo)函數(shù)蛻變?yōu)椋蔎完全由t0時(shí)刻的新息向量d0決定。從表1可見，3DVar_1開始時(shí)J=375.73，說明背景場的觀測相當(dāng)量H[xb(t0)]與觀測量有明顯偏離。極小化結(jié)束后J減小到95.56，下降了 74.57%； 3DVar_2類似，J值由417.54減小到160.0，下降了 61.68%。表明通過三維變分同化觀測資料后，分析結(jié)果向觀測靠攏，背景場得到明顯修正。

作為合理選擇，4DVar_1和4DVar_2控制變量初始值不取w=0，而取自3DVar_1得到的w，這時(shí)4DVar_1和4DVar_2的目標(biāo)函數(shù)J分別為281.87和538.66，均比3DVar_1極小化后的J=95.56大很多，說明引入新的觀測資料后，原先由3DVar_1得到的w值已不再是最優(yōu)的。四維變分極小化過程中，擾動(dòng)預(yù)報(bào)模式將分析時(shí)刻t0的 分析增量PUw傳播到新的觀測時(shí)刻ti，得到新的觀測余差Ri=di+HM0→iPUw，如果不能使J值下降，則需要對w進(jìn)行修正。修正的過程是通過伴隨模式將觀測余差傳回分析時(shí)刻，重新計(jì)算目標(biāo)函數(shù)梯度 ?J，尋找使J下降的方向。如此通過擾動(dòng)模式和伴隨模式反復(fù)積分，隱式實(shí)現(xiàn)了背景誤差協(xié)方差的流依賴（Lorenc，2003），最終分析場受多組觀測和模式的反復(fù)約束，在達(dá)到概率意義最優(yōu)的同時(shí)，也因遵循模式內(nèi)在規(guī)律改變，得到物理意義上更好的分析效果。表1中可見，2個(gè)四維變分同化試驗(yàn)均使J值減少超過69%。

表1 不同同化試驗(yàn)?zāi)繕?biāo)函數(shù)的變化Table 1 Variations of cost function value during minimization in different data assimilation experiments

4.2 分析結(jié)果

圖4—6給出了4個(gè)同化試驗(yàn)在850、500和200 hPa的分析增量（風(fēng)場和位勢高度場）。在850 hPa（圖4）和500 hPa（圖5），4個(gè)試驗(yàn)均分析出偶極子型非對稱增量場，即渦旋對。盡管背景場臺(tái)風(fēng)與“實(shí)況”的位置偏離均超過210 km，4個(gè)同化試驗(yàn)得到的臺(tái)風(fēng)分析位置明顯靠近“實(shí)況”，其中2個(gè)4DVar臺(tái)風(fēng)分析位置幾乎與“實(shí)況”重疊。偶極子型的分析增量恰好反映了同化觀測資料之后對背景場的合理修正：背景場臺(tái)風(fēng)中心附近升壓，升壓區(qū)配以反氣旋性環(huán)流增量，而分析場臺(tái)風(fēng)中心附近降壓，降壓區(qū)配以氣旋性環(huán)流增量。500 hPa（圖5）的結(jié)果與850 hPa大致相同。200 hPa（圖6）上的風(fēng)場和高度場增量很小，2個(gè)3DVar試驗(yàn)沒有出現(xiàn)偶極子流型，風(fēng)場和高度場不太配合，但2個(gè)4DVar試驗(yàn)仍給出有風(fēng)壓場配合的弱偶極子環(huán)流。由于200 hPa臺(tái)風(fēng)渦旋較弱，且個(gè)別時(shí)次部分“觀測”資料未通過質(zhì)量控制（最多時(shí)有約24%的氣壓觀測列入黑名單），一定程度對3DVar有影響，故該層的分析有一定不確定性；而4DVar因使用了多時(shí)次觀測，不同時(shí)次的有效觀測可通過擾動(dòng)模式以及伴隨模式傳播分配，對分析場有整體貢獻(xiàn)。由于臺(tái)風(fēng)移動(dòng)主要受大尺度環(huán)流引導(dǎo)，偶極子風(fēng)壓場增量不僅有助于對“差”的背景場進(jìn)行修正，而且渦旋對之間強(qiáng)的通風(fēng)氣流（850 hPa可超過20 m/s）對調(diào)整臺(tái)風(fēng)引導(dǎo)氣流有作用。通風(fēng)氣流的作用可參閱Fiorino等（1989）和Holland（1983）。Mathur（1991）曾將偶極子環(huán)流加到颶風(fēng)渦旋區(qū)域，對改進(jìn)颶風(fēng)路徑預(yù)報(bào)有幫助。

圖4 四個(gè)試驗(yàn)850 hPa分析增量：風(fēng)矢量 （箭頭）、風(fēng)速 （色階） 和位勢高度 （等值線，單位：gpm）（a.3DVar_1，b.3DVar_2，c.4DVar_1，d.4DVar_2；臺(tái)風(fēng)記號黑色為“實(shí)況”位置，藍(lán)色為背景場位置，紫色為分析場位置）Fig.4 Analysis increments at 850 hPa of four experments：wind vectors （arrow），wind speed （shaded） and geopotential height （contour，unit：gpm）（a.3DVar_1，b.3DVar_2，c.4DVar_1，d.4DVar_2；the typhoon symbols in black，blue and purple mark the positions of typhoon center for the "truth"，the background and the analysis，respectively）

圖5 同圖4，但為500 hPaFig.5 Same as Fig.4 but for 500 hPa

圖6 同圖4，但為200 hPaFig.6 Same as Fig.4 but for 200 hPa

圖7給出同化試驗(yàn)得到的850 hPa分析場（風(fēng)和位勢高度）。4個(gè)試驗(yàn)的分析場和背景場（圖3b和d）相比均得到明顯改善，更靠近“真值”（圖3a和c），臺(tái)風(fēng)渦旋的風(fēng)場和高度場配置合理，風(fēng)速不對稱分布也與“真值”比較符合。與3DVar相比，4DVar分析場的臺(tái)風(fēng)渦旋更對稱，146 dagpm閉合線位置和范圍更接近實(shí)況，臺(tái)風(fēng)中心的位置和強(qiáng)度更合理且更接近“實(shí)況”。4個(gè)試驗(yàn)均顯示臺(tái)風(fēng)環(huán)流東偏北處存在曲率較大的倒槽，與背景場臺(tái)風(fēng)中心（藍(lán)色臺(tái)風(fēng)記號）對應(yīng)，反映了分析結(jié)果對背景場的敏感性，可通過多次外循環(huán)即利用分析場更新背景場積分非線性模式，做多次極小化加以改善。

圖7 850 hPa風(fēng)矢、風(fēng)速 （色階，單位：m/s） 和位勢高度 （等值線，單位：dagpm） 的分析結(jié)果（a.3DVar_1，b.3DVar_2，c.4DVar_1，d.4DVar_2；臺(tái)風(fēng)記號同圖4）Fig.7 Analyses of wind voctor，wind speed （shaded，unit：m/s） and geopotential height （contour，unit：dagpm） at 850 hPa（a.3DVar_1，b.3DVar_2，c.4DVar_1，d.4DVar_2；the symbols of typhoon are the same as Fig.4）

圖8是850 hPa分析場與“真值”的差（簡稱“A?O”）。從位勢高度來看，差值大的地方主要集中在背景場臺(tái)風(fēng)中心和臺(tái)風(fēng)“實(shí)況”的位置，背景場臺(tái)風(fēng)中心處為負(fù)偏差，而臺(tái)風(fēng)“實(shí)況”位置處為正偏差。由于分析場是背景場和觀測的加權(quán)平均（權(quán)重分別與背景誤差協(xié)方差和觀測誤差協(xié)方差有關(guān)），因而分析場既不完全靠近背景場也不完全靠近觀測，而是取兩者中間的最佳狀態(tài)，所以圖8出現(xiàn)這樣的偏差結(jié)構(gòu)是自然的（也是圖7中倒槽出現(xiàn)的原因）。4DVar試驗(yàn)相對于3DVar試驗(yàn)，偏差更加集中于臺(tái)風(fēng)中心，反映臺(tái)風(fēng)中心氣壓的影響。3DVar偏差場較發(fā)散，從另一個(gè)角度反映了4DVar的風(fēng)壓場分析效果優(yōu)于3DVar。

圖8 850 hPa 風(fēng)矢 （箭頭）、風(fēng)速 （色階） 和位勢高度 （等值線，單位：gpm） 分析場與“真值”之差（分析?觀測）（a.3DVar_1，b.3DVar_2，c.4DVar_1，d.4DVar_2；臺(tái)風(fēng)記號同圖4）Fig.8 Analysis minus "truth" （A?O） for wind vector （arrow），wind speed （shaded） and geopotential height（contour，unit：gpm） at 850 hPa （a.3DVar_1，b.3DVar_2，c.4DVar_1，d.4DVar_2；the symbols of typhoon are same as Fig.4）

針對2日12時(shí)的分析增量（簡稱“A?B”），在“實(shí)況”臺(tái)風(fēng)中心（17.2°N，117.6°E）附近分別沿17°N和沿117°E做緯向（圖9）和經(jīng)向（圖10）高度剖面。為了和“真值”與背景場的差異（簡稱“O?B”）進(jìn)行對比，圖9和10中增加了“O?B”的剖面。如果分析合理，“A?B”的結(jié)構(gòu)應(yīng)與“O?B”相似。

圖9中，3DVar_1、4DVar_1和4DVar_2 各變量的分析增量“A?B”的主要分布特征是一致的，且與“O?B”相似。 在臺(tái)風(fēng)東側(cè)，u風(fēng)有較強(qiáng)的東風(fēng)增量，v風(fēng)有較強(qiáng)的南風(fēng)增量。在臺(tái)風(fēng)西側(cè)，u風(fēng)有較弱的西風(fēng)增量，v風(fēng)有弱的北風(fēng)增量。對流層中下部，臺(tái)風(fēng)中心附近降壓，東側(cè)升壓。東西向氣壓增量差與強(qiáng)南風(fēng)增量相對應(yīng)。位溫（θ）增量在臺(tái)風(fēng)中心上空20—25層（約400—300 hPa）有暖心結(jié)構(gòu)，比濕（q）在臺(tái)風(fēng)中心附近的15—25層存在弱的增濕。3個(gè)試驗(yàn)結(jié)果的垂直結(jié)構(gòu)很相似，但也存在細(xì)節(jié)上的差異。

圖10與圖9的結(jié)論相同，3DVar 和4DVar的“A?B”垂直結(jié)構(gòu)十分接近，也與“O?B”相似。u風(fēng)在臺(tái)風(fēng)北側(cè)有較強(qiáng)的東風(fēng)增量，南側(cè)有較強(qiáng)的西風(fēng)增量。v風(fēng)在臺(tái)風(fēng)中心及北側(cè)有北風(fēng)增量，南側(cè)低層有弱的南風(fēng)增量。對流層中下部臺(tái)風(fēng)中心附近有較明顯降壓，強(qiáng)南北向氣壓增量差與強(qiáng)東風(fēng)增量相對應(yīng)。位溫在臺(tái)風(fēng)中心上空20—25層有升高，比濕在中低層臺(tái)風(fēng)中心北側(cè)有弱的增大。

圖9 沿17°N穿過臺(tái)風(fēng)中心所做的分析增量緯向高度剖面（a.u，b.v，c.Π，d.θ，e.q； a1—e1.3DVar_1，a2—e2.4DVar_1，a3—e3.4DVar_2，a4—e4.觀測減背景場）Fig.9 Latitude-height cross sections of analysis increments along 17°N（a.u，b.v，c.Π，d.θ，e.q； a1—e1.3DVar_1，a2—e2.4DVar_1，a3—e3.4DVar_2，a4—e4.O?B）

圖10 同圖9，但沿117°EFig.10 Same as Fig.9 but for along 117°E

續(xù)圖9Fig.9 Continued

續(xù)圖9Fig.9 Continued

續(xù)圖10Fig.10 Continued

續(xù)圖10Fig.10 Continued

4.3 同化分析場對預(yù)報(bào)的影響

將4個(gè)試驗(yàn)的分析結(jié)果作為初值進(jìn)行干模式積分，并與對照試驗(yàn)（CTL）和“實(shí)況”（真值）進(jìn)行對比。3DVar_2試驗(yàn)的起報(bào)時(shí)間為2日18時(shí)，其他3個(gè)同化試驗(yàn)（3DVar_1、4DVar_1和4DVar_2）的起報(bào)時(shí)間均為2日12時(shí)。表2是臺(tái)風(fēng)預(yù)報(bào)位置與“實(shí)況”位置的偏差，0 h代表分析時(shí)刻，其余為預(yù)報(bào)時(shí)刻。表中可以看到，3DVar和4DVar均明顯改善了臺(tái)風(fēng)的分析位置。3DVar_1和3DVar_2分析場臺(tái)風(fēng)位置偏差分別為115.1和73.9 km，相比背景場的位置偏差（12時(shí)為230.8 km，18時(shí)為218.1 km），積分初始時(shí)刻的臺(tái)風(fēng)中心位置得到明顯修正。4Dvar_1和4DVar_2分析位置偏差分別只有30.8和54.0 km，明顯好于2個(gè)3DVar。從預(yù)報(bào)位置偏差來看，對照試驗(yàn)路徑預(yù)報(bào)誤差最大，3DVar比對照試驗(yàn)有較大改善，且3DVar_2優(yōu)于3DVar_1，但2個(gè)4DVar試驗(yàn)則明顯優(yōu)于2個(gè)3DVar試驗(yàn)，可見初值質(zhì)量對臺(tái)風(fēng)路徑預(yù)報(bào)有明顯影響。4DVar_2盡管在初始位置上不及4DVar_1，但多個(gè)預(yù)報(bào)位置比4DVar_1好，表明6 h同化時(shí)間窗比3 h同化時(shí)間窗對路徑預(yù)報(bào)有一定改善。

表2 不同同化試驗(yàn)方案的臺(tái)風(fēng)“彩虹”路徑預(yù)報(bào)誤差Table 2 Absolute forecast errors of the path of typhoon Mujigae in different data assimilation experiments

圖11 是不同初值（對照試驗(yàn)和4次同化試驗(yàn)）預(yù)報(bào)的“彩虹”路徑和“真實(shí)”路徑。以同化分析場作為初值，“彩虹”的路徑預(yù)報(bào)均比以背景場為初值有明顯改善，其中2個(gè)4DVar試驗(yàn)的預(yù)報(bào)路徑更接近“實(shí)況”（圖11a），2個(gè)3DVar試驗(yàn)的預(yù)報(bào)路徑與“實(shí)況”有一定偏離。分析時(shí)刻為18時(shí)的3DVar_2路徑預(yù)報(bào)誤差小于分析時(shí)刻為12時(shí)的3DVar_1，但仍明顯不及分析時(shí)刻為12時(shí)的4DVar試驗(yàn)。

圖11 “彩虹”臺(tái)風(fēng)的路徑預(yù)報(bào) （a.起報(bào)時(shí)間2015年10月2日12時(shí)，b.起報(bào)時(shí)間2015年10月2日18時(shí)）Fig.11 Track forecasts of typhoon Mujigae（a.initial time at 12：00 UTC 2 October 2015，b.initial time at 18：00 UTC 2 October 2015）

圖12是4個(gè)同化試驗(yàn)和對照試驗(yàn)給出的3日18時(shí)的預(yù)報(bào)。其中3DVar_1、4DVar_1和4DVar_2是30 h預(yù)報(bào)，3DVar_2是24 h預(yù)報(bào)。對照試驗(yàn)預(yù)報(bào)（圖12e）臺(tái)風(fēng)強(qiáng)度與“真值”（圖12f）相比明顯弱很多，臺(tái)風(fēng)登陸時(shí)間明顯偏早。4個(gè)同化試驗(yàn)均比

對照試驗(yàn)有顯著改善（圖12a—d），臺(tái)風(fēng)強(qiáng)度和風(fēng)速非對稱分布更像“真值”（圖12f）。3DVar_2（圖12b）較3DVar_1（圖12a）渦旋結(jié)構(gòu)更接近實(shí)況，顯示同是3DVar，基于最近時(shí)刻分析場的預(yù)報(bào)優(yōu)于更早時(shí)刻分析場的預(yù)報(bào)。4DVar的預(yù)報(bào)比3DVar試驗(yàn)更好一些，風(fēng)場分布細(xì)節(jié)上4DVar_2略優(yōu)于4DVar_1。

圖12 基于不同初始場預(yù)報(bào)的2015年10月3日18時(shí)850 hPa“彩虹”臺(tái)風(fēng)的風(fēng)場 （風(fēng)矢和色階） 和位勢高度場（等值線，單位：dagpm）（a.3DVar_1，b.3DVar_2，c.4DVar_1，d.4DVar_2，e.CTL，f.truth）Fig.12 Model forecasts of wind （arrow and shaded） and geopotential height （contour，unit：dagpm） at 850 hPa at 18：00 UTC 3 October 2015 for typhoon Mujigae （a.3DVar_1，b.3DVar_2，c.4DVar_1，d.4DVar_2，e.CTL，f.truth）

圖13是不同初值得到的850 hPa風(fēng)場和高度場的均方根誤差隨預(yù)報(bào)時(shí)間變化?？梢钥闯鲆?DVar分析場為初值，預(yù)報(bào)場的均方根誤差整體明顯下降，預(yù)報(bào)效果明顯比對照試驗(yàn)好很多，也優(yōu)于3DVar。除了12 h前的預(yù)報(bào)，3DVar_2整體優(yōu)于3DVar_1，尤其是在高度場的預(yù)報(bào)上。但是計(jì)算結(jié)果也顯示，同化3個(gè)時(shí)次資料的4DVar_2并不比同化2個(gè)時(shí)次資料的4DVar_1均方根誤差小，原因可能與短的同化時(shí)間帶來模式積分的噪音有關(guān)，留待今后更加細(xì)致地分析。

圖13 850 hPa 水平風(fēng)場分量（a.u，b.v）和位勢高度（c）均方根誤差隨預(yù)報(bào)時(shí)效變化Fig.13 Temporal variations of RMSEs for forecasts of horizontal wind components（a.u，b.v）and geopotential height（c）at 850 hPa

5 總結(jié)

在GRAPES區(qū)域擾動(dòng)預(yù)報(bào)模式基礎(chǔ)上開發(fā)了四維變分同化（4DVar）系統(tǒng)，該同化系統(tǒng)采用增量分析策略建立，是一個(gè)三維和四維一體化的區(qū)域同化系統(tǒng)框架。與原有GRAPES 3DVar相比，增加溫度作為第5個(gè)控制變量，同時(shí)采用了地形追隨坐標(biāo)下的平衡方程，即建立氣壓（Π）和流函數(shù) (ψ)有線性約束關(guān)系的泊松方程，并在模式面上求解，相比于原來GRAPES 3DVar減少了變量轉(zhuǎn)換和空間插值帶來的誤差積累。利用與擾動(dòng)模式相同的廣義共軛余差（GCR）方法求解亥姆霍茲方程的伴隨方程，保證了伴隨模式的精度要求。

為評估所建立的四維變分同化系統(tǒng)，針對2015年10月2—3日臺(tái)風(fēng)“彩虹”進(jìn)行了理想試驗(yàn)。利用在中國南海及周邊區(qū)域人為設(shè)置的25個(gè)假想“探空”站，通過運(yùn)行關(guān)掉物理過程選項(xiàng)的GRAPES非線性模式（干模式），產(chǎn)生測試需要的“真實(shí)”大氣，并用于生成理想人造“探空”觀測資料。共設(shè)計(jì)了4個(gè)同化試驗(yàn)，即2個(gè)三維變分同化試驗(yàn)及2個(gè)分別使用3和6 h時(shí)間窗的四維變分同化試驗(yàn)。利用渦旋分離和重定位技術(shù)產(chǎn)生一個(gè)明顯變“差”的背景場。

試驗(yàn)表明，4個(gè)同化試驗(yàn)均在對流層中低層臺(tái)風(fēng)環(huán)流附近分析出偶極子型風(fēng)壓增量場，在背景場臺(tái)風(fēng)中心附近為反氣旋高壓增量，分析場臺(tái)風(fēng)中心附近為氣旋性低壓增量。偶極子增量對初始背景場產(chǎn)生了很好的修正作用，使得分析場的臺(tái)風(fēng)渦旋更接近觀測。垂直剖面圖顯示分析增量符合臺(tái)風(fēng)環(huán)流物理量場基本分布特征。使用同化得到的分析場作為初值，不論是三維變分還是四維變分，GRAPES模式對臺(tái)風(fēng)“彩虹”的路徑和強(qiáng)度預(yù)報(bào)都得到了明顯改進(jìn)。四維變分同化試驗(yàn)比三維變分同化試驗(yàn)有顯著的分析和預(yù)報(bào)效果。試驗(yàn)也表明，6 h同化時(shí)間窗僅在臺(tái)風(fēng)路徑預(yù)報(bào)上比3 h同化時(shí)間窗略優(yōu)。測試表明，開發(fā)的四維變分同化框架得到了符合理論預(yù)期的分析結(jié)果，同化更多資料并反復(fù)受到模式約束的四維變分同化系統(tǒng)確實(shí)能有效改善初值質(zhì)量，進(jìn)而改善區(qū)域數(shù)值預(yù)報(bào)。本項(xiàng)工作為進(jìn)一步發(fā)展有完整物理過程的區(qū)域四維變分同化系統(tǒng)打下科學(xué)基礎(chǔ)。

由于缺乏實(shí)際觀測資料，文中對于4DVar系統(tǒng)框架的測試僅基于“理想”大氣進(jìn)行，人為構(gòu)造的觀測資料也不能完全反映實(shí)際大氣的狀態(tài)。另外，系統(tǒng)框架也沒有包含完整的物理過程，后續(xù)工作還要進(jìn)一步開發(fā)各種主要的線性物理過程方案，并通過實(shí)際觀測資料進(jìn)一步檢驗(yàn)。