貝葉斯膨脹算法對EnSRF雷達資料同化的影響研究

高士博 閔錦忠 黃丹蓮

?

貝葉斯膨脹算法對EnSRF雷達資料同化的影響研究

高士博1, 2閔錦忠1, 2黃丹蓮1, 2

1南京信息工程大學氣象災害預報預警與評估協同創新中心，南京 210044;2南京信息工程大學氣象災害教育部重點實驗室，南京210044

本文針對2009年6月5日發生在我國華東地區的一次中尺度對流過程（Mesoscale Convective System，簡稱MCS），基于集合均方根濾波（Ensemble Square Root Filter，簡稱EnSRF）方法同化多部多普勒天氣雷達資料，引入具有時空自適應理論優勢的貝葉斯膨脹算法，通過與常數膨脹算法的對比，分析了兩種膨脹算法對EnSRF同化效果的影響。結果表明：貝葉斯膨脹算法同化的雷達組合反射率因子在強對流中心處有所增強，改善了基于常數膨脹算法的EnSRF同化強對流系統偏弱的問題。相比常數膨脹算法，貝葉斯膨脹算法同化的冷池結構更合理，徑向風和反射率因子的均方根誤差均減少。進一步探討貝葉斯膨脹算法對同化效果改善的原因，結果發現：貝葉斯膨脹參數的分布與反射率因子的均方根誤差分布十分吻合，這表明貝葉斯膨脹算法可以在背景場均方根誤差較大，即背景場與觀測差距較大時，給出較大的膨脹參數，進而增加集合的背景場誤差，使得觀測權重增大，從而給出了較大的分析增量。對集合平均分析場進行了1小時的確定性預報發現，貝葉斯膨脹算法提高了預報模式對安徽與江蘇交界處的強對流系統的模擬效果，回波強度更強，冷池強度和范圍更大，且對于不同組合反射率因子的閥值，貝葉斯膨脹算法的評分（Equitable Threat Score，簡稱ETS）均高于常數膨脹算法。這表明貝葉斯膨脹算法有效地改進了基于常數膨脹算法的EnSRF同化雷達資料的效果。

集合卡爾曼濾波 雷達資料同化 貝葉斯膨脹算法 常數膨脹算法

1 引言

中尺度對流系統（Mesoscale Convective System，簡稱MCS）是水平尺度在2～2000 km左右的具有旺盛對流運動的天氣系統，在衛星云圖上表現為一塊連續的高亮的云團，且具有對流核的云結構（Houze, 1993），在雷達回波圖上，回波強度超過40 dB，伸展范圍大于100 km，持續時間為3～24 h（Schumacher and Johnson, 2005）。MCS是產生雷暴、大風、暴雨和冰雹等強對流災害性天氣的重要天氣系統，因此對其進行準確地預報就顯得尤為重要。對于中尺度對流天氣的數值預報和模擬，初值的獲取很重要，而高時空分辨率的多普勒雷達資料包含了豐富的中小尺度信息，為數值天氣預報（Numerical Weather Prediction，簡稱NWP）模式提供了合理和準確的初始場。

多普勒天氣雷達觀測的是非常規氣象數據，需要從探測到的徑向風和反射率因子來估計實際大氣的狀態，得到風、溫、壓、濕等模式變量信息。目前主要的雷達資料同化方法有三維、四維變分法和集合卡爾曼濾波（Ensemble Kalman Filter，簡稱EnKF）等。EnKF結合了傳統卡爾曼濾波和集合預報的優點，通過預報集合樣本估計得到具流依賴的背景誤差協方差，與變分同化方法相比，EnKF算法簡單，不要求編寫切線性和伴隨模式，而且可以顯式地提供集合預報的初始擾動，已經成為國內外熱門的資料同化方法（許小永等, 2006; Evensen, 2003; Zhang et al., 2003; Houtekamer and Mitchell 2005; Szunyogh et al., 2008; Torn and Hakim, 2008; Aksoy et al., 2009, 2010; Bonavita et al., 2010; 秦琰琰等, 2012）。近年來，很多學者將EnKF應用到對流尺度的資料同化上，如Snyder and Zhang（2003）利用EnKF同化模擬雷達資料，初步證明了將EnKF應用于對流尺度上同化雷達資料是可行的。Tong and Xue（2005）和 Xue et al.（2006）在此基礎上進一步考慮了包含復雜冰相作用的微物理過程，利用ARPS（Advanced Regional Prediction System）模式及其EnKF系統同化了一次風暴過程的多普勒雷達資料，取得了較好的效果。Huet al.（2006a, 2006b）發現同化雷達資料能夠更加準確地預報龍卷的位置，并且反射率的作用更加明顯。蘭偉仁等（2010）利用模擬多普勒雷達資料進行一系列風暴尺度的EnKF敏感性試驗，檢驗EnKF在風暴天氣尺度資料同化方面的效果，并分析了集合卡爾曼濾波各參數對同化效果的影響。

需要注意的是，盡管EnKF在對流尺度同化上取得了較好的效果，但它在實際雷達資料同化的應用中，還是不可避免產生濾波發散的問題。EnKF是通過集合來統計觀測增量與集合擾動之間的相關系數的，當集合成員數較小時，相關系數會存在較大的采樣誤差。另外，由于數值模式的不確定性，模式積分預報和參數化的物理過程會引入模式誤差（Miyoshi et al., 2011）。在EnKF 中若不能正確考慮采樣誤差和模式誤差，預報集合的統計協方差將低估真正的預報誤差協方差，集合離散度會迅速減少，使得觀測對同化的影響越來越小，甚至完全失效，造成濾波發散（Houtekamer and Mitchell, 1998; Anderson and Anderson, 1999; Constantinescu, et al. 2007）。為此，一種旨在處理采樣誤差和模式誤差的協方差膨脹技術逐漸發展起來。

最早Anderson and Anderson（1999）提出具有經驗性質的協方差膨脹算法，即通過一個略大于1的膨脹參數來增大背景誤差協方差，但膨脹參數的選擇依賴集合成員數和分析的時間間隔等，需要反復試驗才能得到，在復雜的NWP模式中，計算代價巨大。膨脹參數在時間和空間上都是常數，在觀測稀少的地方會產生多余的離散度。隨后，Zhang et al.（2004）采用了同時保留先驗擾動和后驗擾動的松弛膨脹技術，這種方法避免了對沒有觀測同化更新的集合擾動進行不必要的方差擴張。Whitaker et al.（2008）以及Houtekamer et al.（2009）發現，若對先前的分析成員添加隨機的擾動，能夠較好地處理協方差低估的不均勻性。然而，以上方法均需事先給定最優的參數，如振幅和權重系數等。另外，當天氣系統發展變化快、觀測信息密度比較大時，對所有變量使用空間一致的膨脹參數并不合理，在實際應用中誤差的特性在空間上是各項異性的。近年來，氣象學者致力于發展具有自適應性質的協方差膨脹算法（Miyoshi, 2011; Li et al. 2009; Bishop and Hodyss, 2009a, 2009b）。Anderson（2007）提出一種自適應協方差膨脹算法，膨脹參數的計算公式與貝葉斯濾波更新方程類似，尋找在已知觀測的條件下最有可能出現的期望值，下一時次同化的膨脹參數為上一時次計算得到的期望值。該膨脹參數在時間上具有自適應性，但沒有考慮膨脹參數在模式空間上的不均勻性。為此，Anderson（2009）進一步利用泰勒級數展開，先求解膨脹參數的觀測似然函數，再計算貝葉斯膨脹參數的后驗概率的極大值得到最佳的膨脹參數，從而實現了膨脹參數在時間和空間上的自適應，但該膨脹算法只是在簡單的Lorenz_ 96進行試驗。

由于集合平方根濾波（Ensemble Square Root Filter，簡稱EnSRF）方法能夠避免傳統EnKF由于擾動觀測而造成的采樣誤差（Whitaker and Hamill, 2002），因此本文通過ARPS EnSRF系統，針對一次MCS的發生發展過程，將Anderson（2009）提出的具有時空自適應理論優勢的貝葉斯膨脹算法應用到EnSRF雷達資料同化上，通過與常數膨脹算法的對比，全面評估了貝葉斯膨脹算法在EnSRF雷達資料同化上的適用性。

2 膨脹算法

（1）常數膨脹算法

Anderson andAnderson（1999）提出的常數膨脹算法公式為

其中，表示模式變量，上標inf表示進行了協方差膨脹，下標表示第個集合成員，下標表示第個模式變量，表示模式變量的總數，表示集合成員總數。是根據經驗決定的協方差膨脹參數，不隨時間變化，這對于發展變化比較快的中尺度天氣系統是不合理的。另外，在空間上也是各向均勻的，在觀測稀疏的區域會給出多余的離散度，導致模式積分不穩定，而在觀測稠密的區域又會給出較低的離散度，從而降低同化效果。

（2）貝葉斯膨脹算法

Anderson（2009）進一步提出貝葉斯膨脹算法，基本膨脹公式形式與公式（1）相同，不同的是依據貝葉斯理論來估計膨脹因子的概率分布，通過令其后驗概率分布取得極大值，最終求解得到膨脹參數值。

由貝葉斯公式，膨脹參數的后驗概率為

其中，表示概率密度函數，表示當前時刻，表示由模式分析值計算出的觀測變量，中表示觀測變量，上標o表示實際觀測，下標表示當前時刻經過同化后的值，下標?1表示前一個時刻的分析值。假設先驗概率分布符合正態分布，即其概率密度函數為

在時間自適應的基礎上，進一步考慮膨脹參數在模式空間和觀測空間的關系，引入觀測空間膨脹參數。其中與的關系為：，表示數學期望，為觀測和狀態量的相關系數與局地化因子的乘積。通過可以建立膨脹參數在模式空間和觀測空間的關系，從而實現了膨脹參數在時間和空間上的自適應。

公式（6）可簡化為

其中，

公式（7）、（8）、（9）中：

即. （14）

由上可知，根據公式（6）到（14）求得的貝葉斯膨脹參數能夠根據觀測而改變，是格點，變量和時間的函數，且是在采樣誤差最小化的意義下理論上的最優解，因此具有時空自適應的理論優勢。

3 試驗設計

2009年6月5日，華東地區發生一次典型的MCS過程，大部分地區出現強雷暴天氣，部分地區伴隨大風、冰雹，過程期間有35個站點降水強度超過10 mm h?1，9個站點超過20 mm h?1。這次MCS過程系統移動快，影響范圍廣，造成了非常嚴重的人員傷亡和財產損失。

本文采用非靜力中尺度數值模式ARPS對此次MCS過程進行數值模擬，將NCEP（National Centers for Environmental Prediction）提供的6小時間隔，水平分辨率為0.5°×0.5°的GFS（Global Forecast System）資料作為初始和邊界條件，同時利用EnSRF同化方法對4部S波段天氣雷達，包括南京、鹽城、南通以及合肥進行雷達資料同化。在進行多普勒雷達資料同化之前，需要對雷達資料進行質量控制，包括去除噪聲點、地物回波、二次回波、速度退模糊（Oye et al., 1995）和去距離折疊。最后將質量控制后的雷達資料插值到模式網格上，在插值過程中進行了雷達觀測資料的稀疏化。

為比較不同膨脹算法對雷達資料同化效果的影響，分別設計三組試驗（表1）。三組試驗均采用雙層單向嵌套網格，外層網格水平分辨率為3 km，水平格點數為303×303，內層網格水平分辨率為 1 km，水平格點數為803×803。兩層嵌套的中心均位于（33°N，120°E），垂直方向均取53個不等距分層，垂直平均分辨率為425 m。兩層區域均覆蓋中國江蘇、安徽、山東和浙江大部分地區，地形分辨率為30 s（圖1），且均采用Lin-Tao冰相微物理參數化方案（Hong et al., 2004），1.5階湍流動能方案，NASA大氣輻射傳輸方案和兩層土壤模式等，內外層區域均不使用積云對流參數化方案。ARPS模式只能實現單向嵌套，所以內層區域的模擬是以外層區域的模式輸出資料作為側邊界條件的。

圖1 模擬區域與雷達位置（圖中黑色矩形框代表模擬范圍，倒三角代表雷達位置，大圓圈為半徑為230 km雷達觀測的范圍，其中4部雷達分別為南京（NJRD）、南通（NTRD）、鹽城（YCRD）和合肥（HFRD）雷達

表1 試驗方案

第一組試驗為控制試驗（Experiment of Control，簡稱ExpC），只模擬不同化，外層區域采用6小時間隔的GFS再分析資料作為初始場和邊界條件，從6月5日09:00 UTC（協調世界時，下同）開始模擬，積分2 h，將10:00 UTC的模擬結果插值到內層區域，并利用外層區域提供的側邊界條件，對內層區域進行模擬，從10:00 UTC模擬到11:00 UTC。

第二組試驗為常數膨脹算法同化試驗（Experiment of Multipicative inflation, 簡稱ExpM），首先對6月5日00:00 UTC的外層區域先進行9小時“spin up”積分再同化，然后在5日09:00 UTC的預報場加入初始擾動后產生40個集合成員。參考Putnam et al.（2014）的方法，進行了兩次高斯隨機平滑擾動。表2為兩次高斯隨機平滑擾動的平滑半徑和變量擾動的標準差，可以看出，第一次擾動在是全場范圍內進行的振幅較小的擾動，由于擾動強度較弱，平滑半徑較小。對全場進行擾動是因為MCS在發展過程不斷有新的對流區域生成，可以使每個格點處的集合離散度都大于0。第二次擾動是在對流區域（≥10 dB）范圍進行振幅較大擾動，由于擾動強度較強，所以平滑半徑較大，在對流區域進行較強的擾動是因為在對流區域存在較大的誤差，包括風場，溫度和水物質等。09:06 UTC開始同化雷達資料，每6分鐘同化一次，10:00 UTC結束同化。協方差局地化為五階距離相關函數方案（Gaspari and Cohn, 1999），水平距離為8 km，垂直距離為4 km，協方差膨脹算法為常數膨脹算法，最后在內層區域上對10:00 UTC的集合平均分析場進行1小時的確定性預報。

表2 初始擾動方案

注：、、分別為擾動風場的緯向分量，經向分量，垂直分量p為擾動位溫q為水汽混合比q為云水混合比q為雨水混合比q為冰晶混合比q為雪混合比q為冰雹混合比

第三組試驗為貝葉斯膨脹算法試驗（Experiment of Bayes inflation，簡稱ExpB），設計同第二組試驗，不同的是采用貝葉斯自適應膨脹算法，同樣對EnSRF的集合平均分析場的進行1小時的確定性預報。

4 試驗結果

4.1 觀測與模擬分析

從觀測雷達組合反射率因子隨時間的演變可以看出（圖2），2009年6月5日09:00 UTC，MCS主體部分位于江蘇省和山東南部，內含許多小尺度的對流單體，最強回波強度在60 dB以上。在MCS主體的西側，即江蘇與安徽交界處，存在一條呈西北—東南走向的線狀強對流回波區，南北方向的長度約為400 km，東西方向的寬度約為20 km。在MCS主體的南側，強對流回波中心為弓形，其中在上海附近的最強回波強度達60 dB以上。此外，在合肥附近還存在一條尺度較小的對流系統向東北方向移動，且在移動的過程中，前方不斷有新的對流單體生成。可見，MCS分別在江蘇和合肥附近存在兩個孤立的對流系統。5日09:30 UTC，這兩個對流系統進一步加強，合肥和上海附近的50 dB以上的強對流回波區范圍明顯增大。5日10:00 UTC，合肥附近的對流系統與MCS主體西側的對流單體逐漸靠攏，半小時后，形成一條東西跨度更大的強對流回波區，同時在MCS主體南側的弓形回波面積增大，層云區結構更緊致，但MCS主體中部的回波強度有所減弱。5日11:00 UTC，MCS主體南部的弓形回波發展成熟，具有高度的組織性，而新形成的強對流系統范圍開始減弱，結構變松散。圖3為控制試驗模擬的雷達組合反射率因子，與圖2a–e對比可知，控制試驗模擬的MCS與觀測有較大的偏差，基本沒有模擬出MCS主體的位置和結構。在山東和江蘇交界處模擬出的單體偏強，位于江蘇南部的弓形回波沒有模擬出來。在合肥附近模擬的對流單體組織松散，面積偏小。

圖2 2009年6月5日雷達觀測組合反射率因子（單位：dBZ）：（a）09:00（協調世界時，下同）；（b）09:30；（c）10:00；（d）10:30 ；（e）11:00

圖3 同圖2，但為模式模擬

4.2 同化結果分析

以上述分析表明，控制試驗只能模擬出部分強對流回波區，與觀測有較大的偏差。圖4分別基于常數膨脹算法和貝葉斯膨脹算法同化得到集合平均場。對比圖4a–f與圖3b、c可以發現，兩組同化試驗的組合反射率因子與實況圖在空間上十分相似，明顯優于控制試驗，這說明EnSRF有效地同化了雷達資料，提高了初始場的準確程度。在5日09:30 UTC，兩組試驗均能分析出江蘇境內團狀的中尺度對流系統，其內部包含多個強度不同的對流云。在主對流系統的后部，即在合肥附近也同化出一個正在發展的對流單體向東北方向移動。相比ExpM，ExpB得到的對流單體和實況場更接近，尤其在合肥附近以及江蘇和安徽交界處。在合肥（33°N，117°E）附近，ExpM同化的組合反射率因子偏弱，約55 dB，而ExpB為60 dB。另外，ExpB在江蘇和安徽省交界處（33°N，119°E）線狀結構明顯，而ExpM在該處的對流范圍偏小，強度偏弱。在5日10:00 UTC（圖4d），合肥附近的對流單體并入MCS主體，且團狀MCS的強度繼續加強，雷達回波強度有所增強，水平覆蓋面積也增大。總體來說，兩種膨脹算法基本上同化出MCS的形狀、強度和位置，但在其細節結構，尤其是強對流區域，ExpB得到的最強組合回波的范圍相比ExpM明顯增大，組織性也更強，ExpM得到的對流中心強度偏弱。

圖4 2009年6月5日EnSRF分析的組合反射率因子（單位：dBZ）：（a）09:30觀測；（b）09:30貝葉斯膨脹試驗；（c）09:30常數膨脹試驗；（d）10:00觀測；（e）10:00貝葉斯膨脹試驗；（f）10:00常數膨脹試驗

為了更清楚地看出不同膨脹算法同化效果的差異，圖5a–d為ExpB與ExpM組合反射率因子與觀測的偏差圖。由圖可知，兩組試驗在江蘇和安徽的交界處存在負偏差，在山東和江蘇的交界處存在正偏差。5日09:30 UTC時，相比于ExpM，ExpB與觀測的負偏差的絕對值略小，在10:00 UTC時，ExpB的負偏差的絕對值進一步減小，且正偏差也有不同程度的減少。可見，貝葉斯膨脹算法通過減少采樣誤差，明顯減少了與觀測的偏差，且隨著時間的增加更加明顯。

圖5 2009年6月5日EnSRF分析的組合反射率因子與觀測的差（單位：dBZ）：（a）09:30常數膨脹試驗與觀測的差；（b）09:30貝葉斯膨脹試驗與觀測的差；（c）10:00常數膨脹試驗與觀測的差；（d）10:00貝葉斯膨脹試驗與觀測的差

冷池是由于降水蒸發冷卻導致冷空氣不斷下沉擴展而形成的近地面冷空氣堆（陳明軒等，2012）。冷池的強度和低層風切變是決定MCS強度和結構的主要因素，冷池的冷出流與低層環境風的輻合作用能夠使對流運動長時間維持（張進和談哲敏，2008）。從近地面擾動溫度的模擬結果來看（圖6a、b、d、e），在09:30 UTC，兩組試驗在江蘇的大部以及江蘇與安徽的交界處，近地面的擾動位溫場均存在負中心，表明此時低層已經形成明顯的冷池結構，其中最強冷池中心位于上海（30.3°N，121°E）和江蘇與山東省交界處（34.2°N，119°E），最低擾動位溫能達到－8 K。09:30 UTC時，在江蘇和安徽交界處，ExpB的冷池強度略強于ExpM。由圖6c可以看出，兩者差值的強度和范圍隨時間增加。在10:00 UTC時，ExpB的冷池強度和面積明顯強于ExpM（圖6f）。由此可見，貝葉斯膨脹算法同化的冷池面積更大，強度更強，比常數膨脹算法更加合理。

圖6 2009年6月5日EnSRF分析的地面擾動位溫（單位：K）和水平速度（單位：m s?1）：（a）09:30常數膨脹試驗；（b）09:30貝葉斯膨脹試驗；（c）09:30貝葉斯膨脹試驗與常數膨脹試驗的差；（d）10:00常數膨脹試驗；（e）10:00貝葉斯膨脹試驗；（f）10:30貝葉斯膨脹試驗與常數膨脹試驗的差

為了定量地對比分析不同試驗分析場的誤差，以合肥附近的雷達為例，圖7分別計算了兩組試驗在同化前后徑向風和反射率因子的均方根誤差。由圖可知，無論是哪組試驗，預報誤差在剛開始時均較大，同化雷達資料后，所有變量的分析誤差在各個時刻均有不同程度的下降，這表明雷達資料同化有效地降低了預報誤差。

圖7 2009年6月5日09:06 UTC至10:00 UTC EnSRF合肥雷達處的均方根誤差（橫軸為距離09:00 UTC的時長）：（a）徑向風（單位：m s?1）；（b）反射率因子（單位：dBZ）

對比兩組試驗的預報和分析誤差發現，無論是徑向風還是反射率因子，ExpB的誤差均小于ExpM。對于徑向風的預報和分析誤差，在第一次同化循環中，ExpB和ExpM誤差均隨時間迅速減少，之后趨于平緩。ExpB經過一次循環后，預報誤差一直在4 m s?1以下，分析誤差一直在2 m s?1以下。ExpM在48～54 min時，誤差有個增加的過程，之后開始減少，預報誤差減少至4.2 m s?1，分析誤差減少至2.2 m s?1。對于反射率因子預報誤差，ExpB和ExpM在前兩次循環誤差相差較少，均先隨著 時間迅速減少。24 min后，ExpB和ExpM的誤差趨于穩定，且隨著時間的增加，ExpB誤差明顯 小于ExpM。對于反射率因子分析誤差，ExpB和ExpM在前三次循環中緩慢增加，之后逐漸收斂，其中ExpB和ExpM分別收斂至5 m s?1和8 m s?1左右。

綜上可知，采用自適應的貝葉斯膨脹算法試驗得到的徑向風和反射率因子的誤差始終小于常數膨脹算法試驗，徑向風和反射率因子分別經過1和3次同化循環后，觀測和模式變量場之間的協相關關系開始趨于合理，觀測信息能夠對預報場進行準確的調整，誤差逐漸趨于穩定。

4.3 原因分析

以上分析表明，貝葉斯膨脹算法較常數膨脹算法有更明顯的同化優勢。為了探討貝葉斯膨脹算法對同化效果改善的原因，圖8以冰雹混合比變量為例，給出了貝葉斯膨脹參數的均方根在5 km處的水平結構。由圖可見，在江蘇與安徽的交界處和合肥附近，達到1.8以上，自適應的貝葉斯膨脹算法通過提高協方差膨脹參數，給與雷達觀測更多的權重，改善了EnKF基于常數膨脹算法同化強對流系統偏弱的問題。需要注意的是，在MCS南部，兩種試驗模擬的弓形回波相比觀測均偏強，結構也更緊致，主要體現在回波強度在40～45 dB之間的面積增大。因此該處不高，為1.1左右，這都說明貝葉斯膨脹參數的結構基本上是合理的，在時空上是不均勻的膨脹參數更符合實際情況。

圖8 2009年6月5日EnSRF貝葉斯膨脹試驗5 km處冰雹混合比的膨脹參數的均方根：（a）09:30；（b）10:00

圖9進一步給出了ExpB在5 km處反射率因子的均方根誤差空間分布。由圖9可知，5日09:30 UTC，絕大部分地區MCS的反射率因子的均方根誤差小于4 dB，但在合肥、安徽與江蘇的交界處以及上海西部等地的均方根誤差相對比較大，最大值發生在連云港，可達20 dB。5日10:00 UTC，均方根誤差相對減少，在江蘇南部、安徽與江蘇的交界處大部分均方根誤差在6 dB以下，但在部分區域，均方根誤差依然較大，均超過10 dB。對比圖8發現，ExpB的均方根誤差的高值區對應的均大于1.8，這表明當模式背景場與觀測差距較大時，貝葉斯膨脹算法通過給出較大的膨脹參數，給予觀測更大的權重，使得模式背景場分析的反射率因子增強，從而更加接近觀測。

4.4 模擬結果

以上分析表明，基于貝葉斯膨脹算法的EnSRF雷達資料同化，對MCS回波結構和強度具有明顯地改善，且能夠較準確地分析MCS的熱力場和動力場的細致特征。為了從動力上進一歩驗證貝葉斯膨脹方法的有效性，圖10分別基于兩種不同膨脹算法的EnSRF的集合平均分析場作了1小時的確定性預報。

圖10 基于EnSRF分析場模擬的組合反射率因子（單位：dBZ）2009年6月5日：（a）10:30觀測；（b）10:30貝葉斯膨脹試驗；（c）10:30常數膨脹試驗；（d）11:00觀測；（e）11:00貝葉斯膨脹試驗；（f）11:00常數膨脹試驗

5日10:30 UTC，在安徽與江蘇交界處，已經形成一條有組織的西北東南走向的強對流回波帶，且在外圍有大范圍的層狀回波存在。MCS南部受較強西北風的影響，回波帶呈現彎曲，中部向前突出，與MCS的移動方向一致，為典型的弓狀回波。兩組試驗均較好地再現MCS南部的弓狀回波，也基本模擬出其它對流單體的形狀。在安徽和江蘇交界處，相比雷達觀測，ExpB模擬的對流帶組合回波強度偏弱，約為45 dB，但組合回波已呈現組織化特征。ExpM在此處的組合回波強度約40 dB，對流帶組織結構較差，出現兩個孤立的對流中心。5日11:00 UTC弓狀回波中的對流云有所增強，但強對流回波帶和層狀云稍減弱。兩組試驗基本上模擬出MCS的主要特征，尤其是南部弓狀回波與觀測十分接近，但在江蘇鹽城、南通模擬的組合回波偏強，范圍偏大，結構變松散。在江蘇和安徽處交界處，ExpB雖然模擬的對流帶偏弱，但相比ExpM，組合回波大值區的面積增大，組織性更強。

ETS（Equitable Threat Score）評分是衡量某一等級的預報準確率的標準，ETS值越接近1，表明預報的準確率越高。表3給出了ExpM和ExpB在不同等級的組合反射率因子的ETS評分，從表中可以得到，ExpB在大于15、30和45 dB的回波強度的ETS時間平均值分別為0.6087626、0.4728514和0.2218793，ExpM在大于15、30和45 dB的回波強度的ETS時間平均值分別為0.5670509、0.4594605和0.1858251。在5日10:00 UTC和5日11:00 UTC之間，ExpB在弱回波區，次強回波區和強對流回波區的ETS評分均高于ExpM，表明貝葉斯膨脹算法試驗在不同等級的回波強度的模擬水平均要高于常數膨脹算法試驗。

表3 常數膨脹試驗和貝葉斯膨脹試驗的ETS評分

圖11a、b、d、e給出了基于常數膨脹算法和貝葉斯膨脹算法模擬的擾動位溫和水平速度場，可以看出，兩組試驗在江蘇和安徽的交界處、山東與江蘇的交界處、江蘇的中北部以及上海附近存在冷池中心，其中以上海為中心的弓形回波的冷池強度最強。低層的冷池是對流的產物，比較兩組試驗發現，貝葉斯膨脹算法模擬的冷池狀態與模擬的組合反射率因子趨于一致。10:30 UTC時，在江蘇和安徽交界處，ExpB模擬的冷池強度約為－7 K，最低位溫值的范圍增大。ExpM模擬的冷池強度和范圍均比ExpB小和弱，說明常數膨脹算法在此處模擬的對流偏弱。在11:00 UTC時，冷池進一步加強，在合肥、安徽和江蘇省交界處的兩塊冷池中心已經合并為新的強冷池中心。由ExpB與ExpM擾動位溫差值（圖11c、f）也可以明顯看出，貝葉斯膨脹算法試驗比常數膨脹算法試驗模擬的冷池中心更強。

圖11 基于EnSRF分析場模擬的地面擾動位溫（單位：K）和水平速度（單位：m s?1）2009年6月5日（a）10:30常數膨脹試驗；（b）10:30貝葉斯膨脹試驗；（c）10:30貝葉斯膨脹試驗與常數膨脹試驗的差（d）11:00常數膨脹試驗；（e）11:00貝葉斯膨脹試驗；（f）11:00貝葉斯膨脹試驗與常數膨脹試驗的差

5 結論

本文基于ARPS模式，針對2009年6月5日發生在我國華東地區的一次MCS過程，引入最新的貝葉斯膨脹算法到EnSRF雷達資料同化上，通過與常數膨脹算法的對比，探討了兩種膨脹算法對EnSRF同化效果的影響及其原因，并通過數值模擬，進一步驗證了將貝葉斯膨脹算法應用于EnSRF同化雷達資料的可行性，主要得到以下結論：

（1）貝葉斯膨脹算法同化出的對流單體強于常數膨脹算法，對流中心組合反射率因子增強，尤其在安徽與江蘇交界處，合肥附近。貝葉斯膨脹算法得到的冷池結構更合理，存在兩個明顯的低溫中心，常數膨脹算法給出的冷池強度偏弱，面積偏小。貝葉斯膨脹算法試驗得到的徑向風和反射率因子的誤差始終小于常數膨脹算法試驗。

（2）貝葉斯膨脹參數的空間分布與反射率因子的均方根誤差的空間分布一致，反射率因子均方根誤差的高值區對應的膨脹參數均大于1.8，這表明貝葉斯膨脹算法可以在背景場均方根誤差較大，即背景場與觀測差距較大時，給出較大的膨脹參數，進而增加集合的背景場誤差，使得觀測權重增大，從而給出了較大的分析增量。

（3）基于貝葉斯膨脹算法對中尺度對流系統回波強度、熱力場和動力場的準確反演，對10次同化循環后的集合平均場開展預報試驗。結果表明，與常數膨脹算法相比，雷達回波的強度和范圍更大，冷池的結構更合理，尤其在安徽與江蘇的交界處、合肥附近模擬效果較好。通過計算兩種算法的ETS評分發現，貝葉斯膨脹算法的ETS評分普遍高于常數膨脹算法。

本文的分析表明，貝葉斯膨脹算法有效地改進了基于常數膨脹算法的EnSRF同化雷達資料的效果。這是因為貝葉斯膨脹算法具有時間和空間的自適應性的理論優勢，能夠使膨脹參數的時空分布特征隨觀測而變化，有效地改善了背景誤差協方差的結構，從而改進了EnSRF的同化效果。

參考文獻：（References）

Aksoy A, Dowell D C, Snyder C. 2009. A multicase comparative assessment of the ensemble Kalman filter for assimilation of radar observations. Part I: Storm-scale analyses [J]. Mon. Wea. Rev., 137 (6): 1805–1824, doi:10.1175/2008MWR2691.1.

Aksoy A, Dowell D C, Snyder C. 2010. A multicase comparative assessment of the ensemble Kalman filter for assimilation of radar observations. Part II: Short-range ensemble forecasts [J]. Mon. Wea. Rev., 138 (4): 1273–1292, doi:10.1175/2009MWR3086.1.

Anderson J L. 2007. An adaptive covariance inflation error correction algorithm for ensemble filters [J]. Tellus A, 59 (2): 210–224, doi:10.1111/ j.1600-0870.2006.00216.x.

Anderson J L. 2009. Spatially and temporally varying adaptive covariance inflation for ensemble filters [J]. Tellus A, 61 (1): 72–83, doi:10.1111/ j.1600-0870.2008.00361.x.

Anderson J L, Anderson S L. 1999. A Monte Carlo implementation of the nonlinear filtering problem to produce ensemble assimilations and forecasts [J]. Mon. Wea. Rev., 127 (12): 2741–2758, doi:10.1175/1520- 0493(1999)127<2741:AMCIOT>2.0.CO;2.

Bishop C H, Hodyss D. 2009a. Ensemble covariances adaptively localized with ECO-RAP. Part I: Tests on simple error models [J]. Tellus A, 61 (1): 84–96, doi:10.1111/j.1600-0870.2008.00371.x.

Bishop C H, Hodyss D. 2009b. Ensemble covariances adaptively localized with ECO-RAP. Part II: A strategy for the atmosphere [J]. Tellus A, 61 (1): 97–111, doi:10.1111/j.1600-0870.2008.00372.x.

Bonavita M, Torrisi L, Marcucci F. 2010. Ensemble data assimilation with the CNMCA regional forecasting system [J]. Quart. J. Roy. Meteor. Soc., 136 (646): 132–145, doi:10.1002/qj.553.

陳明軒, 王迎春, 肖現, 等2012. 基于雷達資料四維變分同化和三維云模式對一次超級單體風暴發展維持熱動力機制的模擬分析 [J]. 大氣科學, 36 (5): 929–944. Chen Mingxuan, Wang Yingchun, Xiao Xian, et al. 2012. A case simulation analysis on thermodynamical mechanism of supercell storm development using 3-D cloud model and 4-D variational assimilation on radar data [J]. Chinese Journal of Atmospheric Sciences (in Chinese), 36 (5): 929–944, doi:10.3878/j.issn.1006-9895.2012.11132.

Constantinescu E M, Sandu A, Chai T F, et al. 2007. Ensemble-based chemical data assimilation. Part I: General approach [J]. Quart. J. Roy. Meteor. Soc., 133 (626): 1229–1243, doi:10.1002/qj.76.

Evensen G. 2003. The Ensemble Kalman Filter: Theoretical formulation and practical implementation [J]. Ocean Dyn., 53 (4): 343–367, doi:10.1007/ s10236-003-0036-9.

Gaspari G, Cohn S E. 1999. Construction of correlation functions in two and three dimensions [J]. Quart. J. Roy. Meteor. Soc., 125 (554): 723–757, doi:10.1002/qj.49712555417.

Hong S Y, Dudhia J, Chen S H. 2004. A revised approach to ice microphysical processes for the bulk parameterization of clouds and precipitation [J]. Mon. Wea. Rev., 132 (1): 103–120, doi:10.1175/1520- 0493(2004)132<0103:ARATIM>2.0.CO;2.

Houtekamer P L, Mitchell H L. 1998. Data assimilation using an ensemble Kalman filter technique [J]. Mon. Wea. Rev., 126 (3): 796–811, doi:10.1175/1520-0493(1998)126<0796:DAUAEK>2.0.CO;2.

Houtekamer P L, Mitchell H K. 2005. Ensemble Kalman filtering [J]. Quart. J. Roy. Meteor. Soc., 131(613): 3269–3289, doi:10.1256/qj.05.135.

Houtekamer P L, Mitchell H K, Deng X X. 2009. Model error representation in an operational ensemble Kalman filter [J]. Mon. Wea. Rev., 137 (7): 2126–2143, doi:10.1175/2008MWR2737.1.

Houze R A Jr. 1993. Cloud Dynamics [C]. San Diego: Academic Press, 573pp.

Hu M, Xue M, Brewster K. 2006a. 3DVAR and cloud analysis with WSR- 88D Level-II data for the prediction of the Fort Worth, Texas, tornadic thunderstorms. Part I: Cloud analysis and its impact [J]. Mon. Wea. Rev., 134 (2): 675–698, doi:10.1175/MWR3092.1.

Hu M, Xue M, Gao J D, et al. 2006b. 3DVAR and cloud analysis with WSR-88D Level-II data for the prediction of the Fort Worth, Texas, Tornadic thunderstorms. Part II: Impact of radial velocity analysis via 3DVAR [J]. Mon. Wea. Rev., 134 (2): 699–721, doi:10.1175/MWR3093.1.

蘭偉仁, 朱江, Xue M, 等. 2010. 風暴尺度天氣下利用集合卡爾曼濾波模擬多普勒雷達資料同化試驗I.——不考慮模式誤差的情形 [J]. 大氣科學, 34 (3): 640–652. Lan Weiren, Zhu Jiang, Xue M, et al. 2010. Storm-scale ensemble Kalman filter data assimilation experiments using simulated Doppler radar data. Part I: Perfect model tests [J]. Chinese J.ournal of Atmospheric Sciences (in Chinese), 34 (3): 640–652, doi:10.3878/j.issn.1006-9895.2010.03.15.

Li H, Kalnay E, Miyoshi T. 2009. Simultaneous estimation of covariance inflation and observation errors within an ensemble Kalman filter [J]. Quart. J. Roy. Meteor. Soc., 135 (639): 523–533, doi:10.1002/qj.371.

Miyoshi T. 2011. The Gaussian approach to adaptive covariance inflation and its implementation with the local ensemble transform Kalman filter [J]. Mon. Wea. Rev., 139 (5): 1519–1535, doi:10.1175/2010MWR3570.1.

Oye R, Mueller C, Smith S. 1995. Software for radar translation, visualization, editing, and interpolation [C]// Preprints, 27th Conference on Radar Meteorology. Vail, CO: American Meteorological Society, 359–361.

Putnam B J, Xue Ming, Jung Y, et al. 2014. The analysis and prediction of microphysical states and polarimetric radar variables in a mesoscale convective system using double-moment microphysics, multinetwork radar data, and the ensemble Kalman filter [J]. Mon. Wea. Rev., 142 (1): 141–162, doi:10.1175/MWR-D-13-00042.1.

秦琰琰, 龔建東, 李澤椿. 2012. 集合卡爾曼濾波同化多普勒雷達資料的觀測系統模擬試驗 [J]. 氣象, 38 (5): 513–525. Qin Yanyan, Gong Jiandong, Li Zechun. 2012. Assimilation of Doppler radar observations with an Ensemble Kalman filter: OSS experiments [J]. Meteor. Mon. (in Chinese), 38 (5): 513–525.

Schumacher R S, Johnson R H. 2005. Organization and environmental properties of extreme-rain-producing mesoscale convective systems [J]. Mon. Wea. Rev, 133 (4): 961–976, doi:10.1175/MWR2899.1.

Snyder C, Zhang F Q. 2003. Assimilation of simulated Doppler radar observations with an ensemble Kalman filter [J]. Mon. Wea. Rev., 131(8): 1663–1677, doi:10.1175//2555.1.

Szunyogh I, Kostelich E J, Gyarmati G, et al. 2008. A local ensemble transform Kalman filter data assimilation system for the NCEP global model [J]. Tellus A, 60 (1): 113–130, doi:10.1111/j.1600-0870.2007. 00274.x.

Tong M J, Xue M. 2005. Ensemble Kalman filter assimilation of Doppler radar data with a compressible nonhydrostatic model: OSS experiments [J]. Mon. Wea. Rev., 133 (7): 1789–1807, doi:10.1175/MWR2898.1.

Torn R D, Hakim G J. 2008. Performance characteristics of a pseudo- operational ensemble Kalman filter [J]. Mon. Wea. Rev., 136 (10): 3947–3963, doi:10.1175/2008MWR2443.1.

Whitaker J S, Hamill T M. 2002. Ensemble data assimilation without perturbed observations [J]. Mon. Wea. Rev., 130 (7): 1913–1924, doi:10.1175/1520-0493(2002)130<1913:EDAWPO>2.0.CO;2.

Whitaker J S, Hamill T M, Wei X, et al. 2008. Ensemble data assimilation with the NCEP Global Forecast System [J]. Mon. Wea. Rev., 136 (2): 463–482, doi:10.1175/2007MWR2018.1.

許小永, 劉黎平, 鄭國光. 2006. 集合卡爾曼濾波同化多普勒雷達資料的數值試驗 [J]. 大氣科學, 30 (4): 712–728. Xu Xiaoyong, Liu Liping, Zheng Guoguang. 2006. Numerical experiment of assimilation of Doppler radar data with an ensemble Kalman filter [J]. Chinese Journal of Atmospheric Sciences (in Chinese), 30 (4): 712–728, doi:10.3878/ j.issn.1006-9895.2006.04.16.

Xue M, Tong M J, Droegemeier K K. 2006. An OSSE framework based on the ensemble square root Kalman filter for evaluating the impact of data from radar networks on thunderstorm analysis and forecasting [J]. J. Atmos. Oceanic Technol., 23 (1): 46–66, doi:10.1175/JTECH1835.1.

張進, 談哲敏. 2008. 啟動對流的初始擾動對熱帶颮線模擬的影響 [J]. 大氣科學, 32 (2): 309–322. Zhang Jin, Tan Zhemin. 2008. Impacts of initial convection-triggering perturbations on numerical simulation of a tropical squall line [J]. Chinese Journal of Atmospheric Sciences (in Chinese), 32 (2): 309–322, doi:10.3878/j.issn.1006-9895.2008.02.10.

Zhang F, Snyder C, Rotunno R. 2003. Effects of moist convection on mesoscale predictability [J]. J. Atmos. Sci., 60 (9): 1173–1185, doi:10.1175/ 1520-0469(2003)060<1173:EOMCOM>2.0.CO;2.

Zhang F, Snyder C, Sun J Z. 2004. Impacts of initial estimate and observation availability on convective-scale data assimilation with an ensemble Kalman filter [J]. Mon. Wea. Rev., 132 (5): 1238–1253, doi:10.1175/1520-0493(2004)132<1238:IOIEAO>2.0.CO;2.

Impact of Bayesian Inflation Method on Assimilation of Doppler Radar Data with EnSRF Method

GAO Shibo1, 2, MIN Jinzhong1, 2, and HUANG Danlian1, 2

1,,2100442,,210044

The mesoscale convective system(MCS) occurred on 5 June 2009 in eastern China is simulated using the Advanced Regional Prediction System (ARPS) model and Doppler Radar data is assimilated with EnSRF. Bayesian inflation method is introduced in this study, which allows the inflation parameter to vary in space and time. The impact of Bayesian inflation method on assimilation of radar data with the ensemble square root filter (EnSRF) is investigated by comparing with the simulation using the multiplicative inflation method. Experimental results show that: the simulated composite reflectivity and cold pool from the Bayesian inflation experiment are stronger than that from the multiplicative inflation experiment; Bayes inflation method improves the performance of EnSRF, which always underestimates convection at the storm center. In the convective region, root mean square innovation of radial velocity and reflectivity in the Bayes inflation experiment are lower than that in the multiplicative inflation experiment. Further analysis indicates that the structure of Bayes inflation parameter corresponds very well to the root mean square innovation of reflectivity, which explains why the performance of EnSRF based on Bayes inflation method is improved. It is found that Bayes inflation method can give more weight to radar observations by increasing background error and provides bigger analysis increment when the root mean square innovation (RMSI) of background is bigger. Simulations of the two analysis fields show that the reflectivity near Hefei is stronger and the convective area of MCS is larger in Bayes inflation experiment than in the multiplicative inflation experiment. The simulated cold pool is colder and the area is bigger from Bayes inflation experiment than from the multiplicative inflation experiment, and corresponds well with observed reflectivity. ETS (Equitable Threat Score) of composite reflectivity from Bayes inflation experiment is higher than that from the multiplicative inflation experiment for various thresold. These resulsts suggest that Bayes inflation method improves the performance of EnSRF in radar data assimilation compared to that based on multiplicative inflation method.

EnSRF method, Radar data assimilation, Bayes inflation method, Multiplicative inflation method

1006-9895(2016)05-1033-15

P446

A

10.3878/j.issn.1006-9895.1511.15230

2015-07-12；網絡預出版日期2015-12-01

高士博，男，1987年出生，博士研究生，主要從事雷達資料同化研究。E-mail: shibogao@126.com

閔錦忠，E-mail: minjz@nuist.edu.cn

國家重點基礎研究發展計劃（973計劃）項目2013CB430102，江蘇省普通高校研究生科研創新計劃項目KYLX_0829、KYLX_0844，國家自然科學基金重點項目41430427，江蘇省高校自然科學重大基礎研究項目11KJA170001

Funded by National Basic Research Program of China (973 Program) (Grant 2013CB430102), the Research Innovation Program for College Graduates of Jiangsu Province (Grants KYLX_0829, KYLX_0844), National Natural Science Foundation of China (Grant 41430427), Key University Science Research Project of Jiangsu Province, the Priority Academic Program Development of Jiangsu Higher Education Institutions (Grant 11KJA170001)