雷達資料同化在可致龍卷的強暴雨 系統分析中的一些應用

傅承浩 郜吉東 王云恒, 唐杰 周長青 葉成志 莊照榮,

（1 湖南省氣象臺，長沙 410118；2 美國國家強風暴實驗室，諾曼 73072，美國；3 美國俄克拉何馬大學中尺度氣象合作研究所，諾曼 OK 73072，美國；4 國家氣象中心，中國氣象局，北京 100081）

0 引言

因為雷達徑向風和數值模式風場之間是線性關系，同化雷達徑向風到對流尺度的數值模式相對來說比較容易，尤其是在變分同化的框架內。在過去30年中，變分資料同化的應用無論是三維變分還是四維變分主要目的是用來初始化數值預報模式[1-11]。但也有很多研究應用在雷達資料的反演過程中[12-20]。在雷達資料同化的研究中，NCAR的Sun等[12, 17]做了很多四維變分同化方面的工作。近年來，集合卡爾曼濾波也在雷達資料同化的應用中取得了很大的進展。本文主要是對作者近10年來在俄克拉荷馬大學強風暴中心和NOAA強風暴實驗室參與的三維變分研究以及在此基礎上發展的集合變分同化系統和混合變分集合同化系統等作些簡要的回顧和展望。

在雷達資料研究中，Gao等[18]較早提出了一個利用兩個以上雷達徑向風場的分析方案。這個方案允許較靈活地同時應用雷達資料和其他與風場有關的資料，例如探空資料、地面站資料和VAD廓線等。這個方案還包括了兩個弱約束：質量連續性方程和平滑算子。無論是在理想試驗還是在實時觀測數據的試驗中，這個方案都得到了較好的結果。但是這個方案也有缺陷，特別是在確定約束項的最優權重過程中存在一些困難。因此，發展了一個更標準化的三維變分同化系統[20]。在這個新系統中，遞歸濾波方法[21]被用于模擬背景場誤差距陣。利用遞歸濾波可以很有效地把每個雷達資料同化到觀測位置附近的模式格點上[22]。同時能降低變分系統的條件數，從而使得數值求解收斂得更快些。對每個弱約束項最優權重的確定是通過對各項相對誤差的估計來實現。后來這個變分系統被應用到很多理想個例研究和實際觀測個例的研究中，取得了較好的效果[23-26]。

以上的研究工作主要集中在對雷達徑向風資料的同化上，但如何有效地利用雷達反射率資料是一個較為困難的課題。首先需要發展一個能把反射率和數值模式水凝物變量聯系起來的較準確的正向算子。因為這個問題和數值模式復雜的微物理過程有關，因而有很多不確定性。其次是反射率和數值模式水凝物變量有很強的非線性關系，這使得在變分問題中數值解的收斂性存在困難。由于這些困難的存在，早期的反射率同化方案選擇利用反射率做云分析。Albers等[27]發展了一個主要利用反射率的“局地分析和預報系統”。這個系統主要利用反射率來估計模式水凝物，包括雨水、雪和冰雹等，同時對云內水汽和溫度做出調整。Zhang等[28]對這個系統做了些改進，使它更適合于對流尺度的暴雨系統。Hu等[23-24]、Schenkman等[26]利用此改進后的云分析系統同化反射率，并利用上述變分同化系統[20]分析雷達徑向風到一個適用于對流尺度的區域預報系統[29-30]。結果顯示該系統對短期強風暴，包含龍卷和冰雹等的超級單體有一定的預報能力。Hu等[23-24]的研究表明反射率的同化對暴雨的發展很重要。但是，云分析方案包含很多關于水凝物和反射率之間的經驗關系，因此存在很多不確定性。這些不確定性使得云分析在對流尺度的天氣數值預報中有一定的局限[31]。

Sun等[6]較早意識到同時定量同化徑向風和反射率的重要性。在他們發展的適合于對流尺度的四維變分同化系統中，對幾種反射率的同化作了系統的比較。在理想實驗的框架里，他們發現最好的反射率同化方案是把反射率轉化為雨水混合比，然后直接把雨水混合比同化到對流尺度數值模式中去。但這一方案的缺點是忽視了反射率中冰相（包括雪、冰雹等）的貢獻，從而在他們的模式只有關于暖雨的微物理過程。他們也嘗試直接同化反射率，發現這樣做會增加問題的非線性，使得四維變分問題的求解更加困難。集合卡爾曼濾波是另外一個很有前途的同化雷達反射率和徑向風的方法[32-39]。一般的結論是用此種方法同化雷達反射率和徑向風資料能提高對流尺度的數值預報模式的預報水平。集合卡爾曼濾波同化雷達反射率資料主要的問題仍然是對流尺度模式的誤差一般較大。因而和反射率有關的模式變量（如雨水、雪、冰雹等的混和比）的誤差會通過協相關傳播到該模式的其他變量（例如溫度、水汽等），從而得到不合理的低層冷池分析結果[38]。

在Gao等[40]的研究中，與冰相微物理過程有關的對流尺度數值模式變量被引入進到三維變分的目標函數中。他們根據對流尺度數值模式的溫度變量來決定雷達反射率的正向算子中水凝物（主要是雨、雪、冰雹等）的分類，這個研究使得在數值模式中對雷達反射率的同化更容易些。三維變分同化方法具有計算效率高、容易并行、編碼簡單等特性，所以比較適用于對流尺度的同化和預報系統。

這篇文章的主要目的是介紹一個基于三維變分同化的適用于對流尺度的雷達資料同化方法，特別是近幾年來這個方法在實時分析和預報中的應用。

1 資料同化方法簡介

1.1 目標函數的定義

正如上一節介紹的，Gao等發展了一個適用于對流尺度雷達資料同化的三維變分方案[18-20, 23-24, 41-44]。不同于大尺度的資料同化系統，這個方案強調利用質量連續性方程[18]、從適用于對流尺度的數值預報方程導出的氣壓診斷方程[44]作為弱約束來建立模式變量之間的聯系。從而在同化中對沒有觀測的變量，例如垂直速度、溫度等也做出一定程度的更新和改善。通過定義一個包括數值模式背景場項、雷達資料的觀測項和弱約束項的目標函數，然后通過應用數學上的下降算法（例如共軛梯度法）使目標函數極小化來求得同化分析結果。同化分析結果可用于初始化對流尺度的數值模式。目前這個三維變分同化系統已經被用來初始化ARPS模式和WRF模式。

三維變分同化的目標函數可寫為以下形式，

右邊第一項定義分析場和模式背景場之間的距離。B是背景誤差協方差距陣。第二項是觀測項，代表分析變量和觀測向量yo之間的距離。在最初的研究中，yo只包含雷達徑向風，現在加入了雷達反射率和一些常規觀測資料及風廓線資料等。關于雷達徑向風和反射率的正向算子將在下文介紹。如果忽略觀測資料之間的相關性，觀測誤差的協方差距陣R可以簡化為對角陣，只包含了與各個觀測相對應的儀器誤差和代表性誤差。第三項，JC(x)可以是包含質量連續性方程和數值模式方程在內的弱約束項。在目前的系統中，分析變量包括風場的3個分量、位溫、氣壓、水汽混合比和模式水凝物（包括云水、云、冰、雨、雪、冰雹的混合比）等。

1.2 徑向風正向算子

在方程（1）右邊第二項的觀測項可包含任何觀測資料，側重雷達資料的應用。根據Doviak等[45]，徑向風算子可用以下公式表示：

式中，vr代表雷達徑向風，r是雷達觀測點與雷達位置的距離，h是雷達觀測點離地球表面的高度，s是雷達射線投影到地面上的距離，是雷達觀測點的徑向角。如果在徑向風算子中也考慮雷達的射線寬度，那么根據Rihan等[46]的研究，較復雜一點的徑向風算子可表述為：

式中，是在某一高度角上的徑向風，是數值模式的垂直格距，G描述了與雷達射線寬度有關的放大率函數[47]。

如果雷達徑向風的觀測誤差為1或2 m/s[48-49]，方程（1）中的第二項可改寫為

式中，H(vr,e)定義為把三維風場變量u，v，w投射到徑向風方向的正向算子（式（2）或（3））和把徑向風從模式格點投影到觀測點上的三維線性插值算子。

1.3 反射率正向算子

早期的三維變分框架內，不包含對反射率的直接同化，反射率的分析是通過云分析來進行的[23-24]。Gao等[40]提出了一個在變分框架內直接同化反射率的方案。在這個方案中，首先發展了一個利用數值模式的溫度背景場來對模式的水凝物進行分類的正向反射率算子。根據Lin等[50]、Gilmore等[51]和Doswell等[38]的研究，反射率因子的正向算子可表述為雨、雪、冰雹混合比的函數如下：

根據Smith等[52]的研究，雨水對反射率的貢獻可表述為

ρ是空氣密度。如果氣溫低于0℃，那么干雪對反射率的貢獻可表述為，

如果氣溫高于0℃，那么濕雪對反射率的貢獻可表述為，

根據Lin 等[50]和 Gilmore 等[51]的研究，冰雹對反射率的貢獻為

用式（5）～（9）把反射率同化到數值模式中去不是很理想，存在著一定的缺陷。這是因為反射率是雨、雪、冰雹混合比三者的函數，直接用以上公式可能引起數值解的欠定性。例如，分析中可能導致在模式的低層產生雪的情況，但實際上模式低層由于溫度較高卻只能有雨。為解決這個問題，式（5）可被改進寫成

式中，α隨背景場溫度的區間（Tb=－5℃和Tb=5℃）在0～1變化，Tb可從數值模式中得到。實際上用式（5）和（10）得到的反射率的因子大致相同。但在式（10）中，隱含著不同種類水凝物對反射率因子的貢獻可以根據模式提供的溫度場有所區分。從而當把反射率同化到數值模式時，使不同種類水凝物（雨、雪和冰雹的混合比）的分析結果更加合理。對反射率因子方程（10）取對數就可計算反射率，

這里，ZdB的單位是dBz。所以目標函數（1）中第二項關于反射率的觀測項可寫為，

類似于式（4），式（12）中的H(ZdB)定義為用方程（5）～（11）把模式變量，雨、雪和冰雹的混合比轉化為反射率的非線性算子和把反射率從模式格點插值到雷達觀測點上線性算子。

這種根據數值模式提供的溫度場來確定水凝物對分射率貢獻的正向算子使得對反射率的同化更精確，一定程度上去掉了解的不確定性。根據以前的對比試驗，用溫度場來確定水凝物的正向算子使得同化的結果更好，特別是對低層雨水的分析更合理。作者正在NOAA的“基于數值模式的預警”研究課題的資助下做進一步試驗，希望將來能把此方法應用到實際業務中。

1.4 目標函數中的弱約束項Jc

在三維變分同化中，一些能把模式變量聯系起來的物理約束（例如質量連續方程）和模式方程以及其導出的相關公式可以作為弱約束結合到目標函數中。在目前的系統中，Jc主要包括兩個弱約束項即：質量連續性方程和從ARPS模式導出的氣壓診斷方程。這里主要介紹前者，后者可以使模式各個變量的分析更平衡[44]。

在這個三維變分同化中，3個風場分量之間的平衡關系是通過一個滯彈性近似的質量連續性方程來實現的，它的定義如下，

式中是給定層的平均空氣密度。是個權重系數，它控制目標函數中這個約束項的相對重要性。一般取值是= 5.0×10-4。這個值也可通過對分析方案的調試得到。JDP定義氣壓診斷方程弱約束項[44]。Gao等[18,20]發現這個弱約束建立了風場3個分量之間的緊密聯系，因而對垂直速度的分析非常有效。正如Gao等[18]在研究中指出的，在三維變分同化方案中把質量連續性方程作為弱約束而不是強約束[53-55]，避免了在垂直方向上積分連續性方程引起的誤差累積。所以垂直速度的分析相對來說能更準確些。

1.5 變量替換和條件化

為了使數值求解更快，或者說使目標函數下降更快，定義控制變量v，使得C=，這樣目標函數（式1）變為如下形式[56-57]：

式中H是H的線性化算子目標函數Jinc關于v的梯度和Hessian算子可從公式（15）導出，

式中，I是單位距陣。從式（17）中可看到，從式（1）到（15）的變換使相對應于新控制變量v的本征值不可能是1，從而降低了條件數，求解過程收斂會更快。

在式（15）中距陣C可分解為

式中，D是背景標準誤差的對角化距陣，L是一個放大的系數。從理論上來說，F是一個對角元素是1，而非對角元素是背景協相關誤差系數的矩陣；但是這個矩陣維數巨大，不可能用在數值求解過程中顯示表示。Purser等[21]提出了一個利用對等空間遞歸濾波的方案來模擬F。一價的遞歸濾波定義如下

Xi是定義在格點i上的初始值，Yi是向前濾波后在格點i的值，Zi是向后濾波后在格點i的值，α是按式（20）定義的系數[58]，

式中，L是水平相關尺度，?是水平格距，N是實施濾波的次數。這個遞歸濾波方法要求必須在正反兩個方向各進行一次，以保證位相誤差為0。這個遞歸濾波方法可分別應用在x、y、z這3個方向上。多尺度濾波可通過多次使用式（19）實現，每次應用不同的系數α。Xie等[59-60]在他們的三維變分研究中證明使用多尺度濾波比單尺度濾波效果會更好。

2 幾個實時個例分析

在上述三維變分的框架上，發展了一個實時的可根據天氣形勢移動分析區域的分析系統。這個系統可用數值預報產品和實時收集的雷達資料進行分析[41]。這個系統有以下優點：1）可用多個雷達的資料，以及任何有關的常規數據；2）自動識別強天氣可能發生的地點并提供有較高空間分辨率（1 km）和時間分辨率（5 min）的分析結果；也可由預報員根據當前天氣形勢確定分析區域；3）這個系統能夠幫助預報員估計超級單體中層中尺度氣旋的強弱從而在發布強天氣預警時提供參考和幫助。這個系統曾經被應用到NOAA災害性天氣試驗基地，并得到預報員的好評[61-62]。下面給出幾個在美國產生強龍卷的災害性天氣試驗個例和一個發生在湖南的強暴雨個例的分析結果。

2.Fig.年在俄克拉荷馬州產生龍卷的一個超級單體

2013年5月20日，一個產生了強龍卷的超級單體襲擊了美國俄克拉荷馬州穆爾市（MOORE）。這個龍卷后來被美國氣象部門歸類于5級——龍卷強度的最高級別。這個龍卷共造成了23個人死亡，377個人受傷，并催毀了很多房屋（圖1）。龍卷產生的瞬時最大風速達到了每小時340 km。龍卷的直徑達2.1 km。根據事后收集的資料，龍卷在當地時間14：56從附近紐爾斯爾市開始形成，在地面穿行了大約27 km，大部分經過穆爾市并于當地時間15：35離開地面。這個產生了龍卷的超級單體被附近的多普勒雷達網監測到，其中包括有4個S波段的WSR-88D雷達，KTLX雷達離風暴系統最近，觀測了系統的中低層的結構；KFDR，KINX和KVNX離風暴較遠，主要觀測到了風暴的高層結構。分析中的背景場是由美國環境預報中心（NCEP）的北美中尺度模式（NAM）提供的。在實時運行中，水平方向的分辨率設為1 km格距，模式區域共有200×200個格點。在垂直方向上有31層，采用逐步遞增的不等距格距，平均格距是400 m。當天系統自動定位的分析中心位于（35.36°N，97.18°W）。分析主要集中在三維風場以及由此導出的垂直渦度和垂直速度上。

圖1 2013年5月20日強龍卷襲擊后的俄克拉荷馬州穆爾市的住宅區 Fig.1 Selected images of tornado damage area in Moore, Oklahoma on 20 May 2013

2013年5月20日早晨，美國天氣服務局暴雨預報中心發布了從密蘇里州東南部到得克薩斯中北部的強天氣預警。當地時間14：00左右，產生強龍卷的超級單體開始在雷達回波中顯現。這個超級單體從出現到成熟發展得很快，在不到一個小時內（14：56），在紐爾斯爾市產生了第一個龍卷（圖2）。當年，作者并沒有在NOAA的強風暴試驗基地每天實時運行分析系統，但在預報的強天氣系統可能性提高時，分析系統會自動啟動。系統在當地時間13：00開始實時接收雷達數據。系統自動定位了靠近穆爾市的分析區域。圖2和圖3顯示了分析得到的位于地面高度以上3 km，從14：00—15：35的水平風場、垂直渦度和雷達反射率、以及相應的垂直剖面圖。在這段時間內，超級單體經歷了從生成、發展到成熟的階段。14：00，在分析區域的西邊界出現了幾個雷達回波較弱的普通暴雨單體（圖2a和圖3a）。產生了穆爾市龍卷的暴雨單體迅速加強，僅僅25 min后，演變成超級單體。在對地面高度以上3 km上空的風場分析表明存在一個很強的中尺度氣旋（圖2b和圖3b）。雖然分析也表明在這個超級單體的北部和南部各存在較強的超級單體，但它們都沒產生龍卷（圖2b和圖2c）。穆爾市的龍卷在14：56下沉接連地面，在14：55的分析表明在雷暴中心存在超強的旋轉和上升運動（圖2e和圖3e）。最強的垂直速度達到20 m/s，最強的垂直渦度達到0.015/s。在近一個小時期間，強大的氣旋一直存在于超級單體內，在垂直方向上延伸幾千米。當地時間15：20，這個超級單體分裂，強度減弱，龍卷于15：35時消散。為了驗證分析是否合理，把從當地時間14：00—17：00分析中得到的垂直渦旋度的軌跡和氣象局事后統計得到的龍卷軌跡一起畫在圖4中。從圖中看出，三維變分分析的垂直渦旋度的軌跡前一段與龍卷造成損害的軌跡基本重合。三維變分的分析區域水平格距是1 km，這樣的分辨率雖不可能分辨龍卷，但與龍卷緊密相關的中尺度氣旋卻能很清晰地被分析出來，并且伴隨龍卷的發生一直非常強大。對這個例中類似的超級單體的預報是個非常困難的事情，因為這種超級單體在很短時間內快速加強并生成龍卷，對預報員提前準確發布強天氣預警具有很大挑戰。但快速三維變分同化的分析產品可以迅速地實時計算出氣流旋轉的強度，從而給預報員發布強天氣預警提供及時幫助。

2.Fig.年在阿拉巴馬州產生龍卷的多個超級單體

圖2 分析得到的俄克拉荷馬穆爾市地面以上3 km的多部雷達合成的反射率（彩色陰影區）、水平風場（風矢量）和垂直渦度（黑色等值線）（紫色實線表明圖3垂直剖面圖的位置） （a） 當地時間14：00，（b）14：20，（c）14：45，（d）14：50，（e）14：55，（f）15：00 Fig.2 The synthesized reflectivity (color shaded) and analyzed horizontal wind fields (vectors), and vertical vorticity (black contours) at 3 km AGL at (a) 14:00, (b) 14:20, (c) 14:45, (d) 14:50, (e) 14:55, and (f) 15:00 CST, 20 May 2013 near Moore, Oklahoma. Maroon line denotes location of cross-sections for vertical slices in Fig.3

圖3 圖例說明同圖2，是圖2中相應位置的垂直剖面圖 Fig.3 Wind vectors and reflectivity (color shaded) through a vertical slice in Fig.2

2011年4月27日，在美國東南部幾個州包括阿拉巴馬、阿肯色、喬治亞、密西西比和田納西等地暴發了劇烈的災害性天氣事件，是美國歷史上死亡人數最多，損失最慘重的天氣事件之一（圖5）。總共有292個龍卷報告，大約造成了60億美元的經濟損失，348個人喪失生命。其中阿拉巴馬州的損失最嚴重，4個龍卷被定性為5級（龍卷所能夠造成破壞的最高級別）。僅在這個州就有238人死于這些龍卷造成的破壞。在該災害事件發生的過程中，實時分析系統雖然在運行，但由于太多的超級單體，分析區域只覆蓋了其中的一小部分。模式格點數為400×300，水平格距為1 km，垂直方向仍為不等距的31層。分析的中心位于（33.9°N，86.4°W），大約在伯明翰市西北80 km的地方。在這個分析區域，一共有7部WSR-88D的雷達資料被應用在分析中。分析區間是從當地時間17：00—18：30。圖6顯示了這段時間超級單體群從西南到東北方向移動。在分析區域中每個超級單體都可能成為一個移動的“小炸彈”，產生龍卷，襲擊沿路的建筑及其他目標。選擇對其中一個超級單體在18時的風場進行較為詳細的分析（圖6中的黑色方框）。這個超級單體產生了級別為4級的龍卷，襲擊了阿拉巴馬州的多拉市（Dora）。在這個時刻，龍卷已經產生。圖7a顯示離地面1.5 km高度上的風場環流疊加在反射率的鉤狀回波上。在中層5 km高度上形成直徑為8～10 km的閉合環流（圖7b），在較高層（10 km高度上）風場的形勢是伴隨著旋轉的向下風方向的輻散。在垂直剖面上存在著大的上升運動區，與強反射率區域重合。最大上升速度達到20 m/s（圖7d）。圖8顯示無論在低層、中層還是高層，暴雨區域的強垂直渦度區和強垂直速度區相重合。表明在這種產生了龍卷的超級單體中存在著極強的與旋轉相伴隨的上升運動。但在垂直方向上最大渦度中心和最強上升中心不重合，前者所在層次低于后者。在很多其他產生龍卷的超級單體的分析中也發現了類似現象。另外，在風暴的上層存在重力波從風暴中心向風暴的下風方向的傳播。把從三維變分分析得到的每個格點垂直方向上最大渦度和最大垂直速度做成一個風暴移動軌跡顯示在圖9中。在分析中得到的對應于超級單體的垂直渦度中心和垂直速度中心的移動與美國天氣服務局暴雨中心得到的龍卷的破壞規跡相重合（圖10）。

圖4 由三維變分分析得到的俄克拉荷馬州穆爾市地面到3 km的當地時間14：00—17：00的垂直渦旋度（Helicity）（圖中深藍線是被龍卷破壞的軌跡圖） Fig.4 3DVAR 0-3 km MSL updraft helicity accumulated over the period from 14:00 to 17:00 CST on 20 May 2013 in Central Oklahoma. The dark blue lines indicate tornado damage path as determined by the National Weather Service Forecast Office at Norman, Oklahoma

圖5 2011年4月27日強龍卷襲擊后的阿拉巴馬州住宅區 Fig.5 Selected images of tornado damage area in Alabama on 27 April 2011

2.Fig.年一個在得克薩斯州產生龍卷的超級單體

2010年4月20日，一個產生了龍卷的超級單體襲擊了得克薩斯州西北部區域，這一區域又被稱為得州刀把區。這個超級單體生成于當地時間17時，終止于當地時間21時，生命期大約為4 h。在這個超級單體的生命期內有2個龍卷生成，美國氣象局暴雨預報中心事后也有收集到關于冰雹和強陣風的報告（圖11），幸運的是并無人員傷亡。這個風暴附近的3個WSR-88D的雷達KAMA、KLBB和KFDX實時觀測到了這個風暴的演變。

這個超級單體是在得州和新墨西哥州邊界形成的。在形成的最初兩個小時內，它的風場旋轉不太強。但到當地時間接近19時，風暴中的中尺度氣旋突然加強。19：04，第一個較弱龍卷形成，之后風暴繼續發展和加強。實時分析系統“跟蹤了”這個風暴的演變。分析得到的水平風場、垂直渦度和反射率都表明風暴從19：30—20：00發展到頂峰（圖12）。在這半小時內，地面以上3 km高度上的最大垂直渦度一直維持在0.01 s-1左右。當地時間19：45和19：50，由垂直渦度表征的中尺度氣旋在垂直方向上延伸達到了10 km，最強的渦度中心和強反射率區域相重合（圖13）。在風暴中心地面上空4 km處存在一個弱回波區和較強的低層入流區（圖13a、13c）。第二個龍卷于19：47生成（NOAA暴雨中心的報告，圖11a）。從19：35到20：00的反射率演變可以看到勾狀云逐漸錮囚或閉合的一個過程（圖14a～14f）。第二個龍卷剛好是在這個錮囚過程中形成的。在這個過程中，這個超級單體也降下了很多冰雹。

2.Fig.年發生在湖南的一場大暴雨過程

北京時間2010年5月6日00：00—06：00，湖南省中南部懷化、婁底等6個地區出現暴雨，其中7個縣級市出現大暴雨，日最大降水出現在婁底新化為147.8 mm；據區域氣象觀測站資料，153鄉鎮自動氣象站出現暴雨，77鄉鎮自動氣象站出現大暴雨。本次降水特點為突發性強、短時雨強大、強降水區呈狹長帶狀分布。此次過程造成婁底、懷化、長沙等8個市（州）27個縣（市、區）336個鄉鎮272.35萬人受災，因災死亡11人，失蹤1人，緊急轉移人口6.01萬人，倒塌房屋9944間，直接經濟損失14.53億元。

圖6 阿拉巴馬州在地面以上3 km的多部雷達合成的反射率（彩色陰影區）、水平風場（風矢量）和垂直渦度（黑色等值線）（圖6e中的黑色小框是圖7和圖8的分析區域） （a）當地時間17：00，（b）17：15，（c）17：30，（d） 17：45，（e）18：00，（f）18：15 Fig.6 The synthesized reflectivity (color shaded) and analyzed horizontal wind fields (vectors), and vertical vorticity (black contours) at 3 km AGL, at (a) 17:00, (b) 17:15, (c) 17:30, (d) 17:45, (e) 18:00, and (f) 18:15 CST, 27 April 2011 near North Alabama and nearby area. The small box in Fig.6e is the domain shown in Fig.7 and Fig.8

鑒于此次暴雨過程從6日當地時間02：00開始加強，至08：00減弱，因此主要分析這個時段的雷達資料。該時段的降水情況如圖15顯示。對湖南境內的兩部多普勒雷達資料進行了三維變分分析，圖16顯示的是2010年5月6日暴雨過程的分析區域和雷達位置，這兩部雷達分別是位于長沙的CINRAD/SA雷達和位于常德的CINRAD/SB雷達。同化的水平分辨率是1 km，平均垂直分辨率為500 m。模式水平格點數為403×403，垂直層數為25。

圖7 從圖6中放大的阿拉巴馬州多拉市在地面以上的多部雷達合成的反射率（彩色陰影區，單位：dBz）和水平風場（風矢量，單位：m/s） （a）地面以上1.5 km，（b）地面以上5 km，（c） 地面以上10 km ，（d）從y=140 km的垂直剖面圖 Fig.7 The synthesized reflectivity (color shaded, unit: dBz) and horizontal wind fields (vectors, unit: m/s) which were zoomed in from the box in Fig.6e, at (a) 1.5 k m AGL, (b) 5 km AGL, (c) 10 km AGL, (d) vertical profile at y=140 km, on 27 April 2011 near Cordova, Alabama

從同化后雷達回波和風場的圖中可以看出（圖17），當地時間6日03時，在湖南中部地區出現從常德南部到益陽西南部和從湘西自治州東南部到懷化中東部兩條回波帶與婁底—湘潭中部—株洲西北部回波帶匯合，呈“人”字型回波帶（圖17a）；04：00時回波北擴，強回波范圍明顯增大，＞45 dBz的回波連成線狀，此時反射率因子達到最強，強回波仍然穩定在湘中一線；05：00—06：00時，回波持續北擴，同時在06：00衡陽北部地區又生成了新的東西向的回波帶，但回波強度明顯減弱（圖17d）。07：00，湘中偏北位置的回波東移減弱，衡陽北部的回波仍然持續（圖17e）；6日08：00后回波逐漸減弱（圖17f）。

多普勒雷達徑向速度無法直接反映風場信息，而同化雷達資料后的風場更為直觀，并且包含更多中尺度信息：03：00出現了一個呈西北—東南向的弱輻合區，這預示著將有一個較強的中尺度對流系統發展（圖17a）。04：00這個輻合區更加明顯，同時，在湘中強回波帶對應的區域，風的流線有明顯的氣旋性彎曲且有風速上的輻合，為強降水的發展提供了有利的動力條件（圖17b）。

需要說明的是，目前中國天氣雷達主要是從美國引進并國產化的新一代多普勒雷達，因此從類型上來說與美國主力雷達系統沒有太大差異。并且WRF，ARPS等模式都能直接處理國內的雷達資料。但國內雷達數據質量還不夠好。造成質量問題的原因很多，首先在客觀上來說，中國地形較為復雜，有些雷達選址不當就會導致擋角影響明顯，往往會出現底層幾個仰角沒有有效觀測數據。比如湖南省懷化雷達的擋角問題就比較明顯。其次，雷達的生產工藝及安裝后的系統參數定標不當，也會影響雷達數據質量，比如雜波偏多，系統顯示不清晰，有時對同一系統在兩個相鄰雷達的回波強度或徑向速度有明顯不同。由于國外雷達數據質量相對較好，因此國外的數值模式對國內的雷達數據往往水土不服，同化結果中，天氣系統顯示不夠清晰，雜波偏多，影響了同化后的預報效果，也不利于天氣系統的分析。因此，同化國內雷達數據時要特別加強雷達質量控制。這些也是今后需要研究的課題。另外做的研究例子還不是很多，將來也應做些關于中國龍卷的例子。

圖8 圖例說明同圖7，圖中彩色陰影區為垂直速度，黑色等值線為垂直渦度（粉色線表明重力波傳播的位置） Fig.8 As same as in Fig.7, color shaded is for vertical velocity and contours are for vertical vorticity (pink curve indicates the propagation of gravity waves)

圖9 阿拉巴馬州中北部及其附近區域，由三維變分分析得到的垂直渦度（a）和垂直速度（b） Fig.9 The composite vertical vorticity and vertical velocity tracks derived from 3DVAR analyses at North and Central Alabama and nearby areas

圖10 2011年4月27日美國東南部的幾個州發生的龍卷、冰雹和大風實況報告 Fig.10 The storm report for tornadoes (red), hails (green) and damaging winds (blue) in southeastern U.S. on 27 April 2011, from NOAA/storm prediction center

3 結論與展望

在過去幾年中，發展了一個基于三維變分的高分辨率（1 km）和高頻率（每5 min）的實時分析系統。這個系統可根據天氣實況確定分析區域，實時讀取分析數據，并做出分析。在過去的一些研究和應用中作者對這個系統不斷地完善。在這個三維變分同化中，對代價函數的定義只包含一個背景場、雷達徑向風和一個弱約束項。其中的背景誤差協方差距陣是由簡單的遞歸濾波來實現。最近幾年在變分同化的框架內增加了對反射率的同化，同時也引進了用集合擾動場來導出背景誤差協方差距陣，和引進更多觀測資料項等，以便使分析更精確穩定，同時有能力對對流尺度的數值預報模式提供初始場。

本文介紹了基于三維變分的實時分析系統。這個系統在NOAA的強災害天氣試驗基地進行了準業務運行。它能根據實時天氣形勢確定分析區域，并利用實時雷達觀測數據和數值預報產品快速做出分析。特別是能實時估計出強風暴中中尺度氣旋環流的強度，從而給預報員發布強天氣預警提供幫助。在試驗時得到很多來試驗基地訪問的預報員的好評。文中典型超級單體和暴雨的分析個例表明這個系統能很好地模擬超級單體的發展和較準確地提供龍卷的預警。但這個分析系統也存在很多需要改進的地方。首先分析產品仍然滯后7～8 min，所以接收數據的時間需縮短，后處理的時間也需縮短。另外這個系統對于較弱的風暴分析不是很合理，還有待提高。

這個三維變分同化系統在其他對流尺度的分析和應用研究方面也發揮了很大作用（如Gao等[63]；Stensrud等[42]；Xue等[64]；Schenkman等[26]）。例如，在美國大氣遙感工程研究中心（The Engineering Research Center for Collaborative Adaptive Sensing of the Atmosphere）建立的X波段雷達網的觀測區域進行的高分辨率低層風場分析和超短期的預報[63, 65]。從2008年以來，這個系統也被應用到NOAA災害性天氣試驗基地春季預報試驗階段來為高分辨率WRF預報模式的確定性和集合預報提供初值。這個試驗的分析預報區域覆蓋全美大陸地區，用到幾乎所有WSR-88D雷達網提供的數據[64, 66]。這個春季試驗主要是為了用計算機和數值模式提高對災害性強風暴天氣的預報能力，從而對NOAA暴雨預報中心的業務提供支撐。

圖11 （a）2010年4月20日得克薩斯州刀把區的龍卷和冰雹實況報告；（b）最內黑框為分析區域（200 km×200 km）（外黑框用來確定所用雷達數目（500×500 km）；有3個WSD-88D雷達（代碼分別是：KAMA、KFDX 和KLBB）被分析采用（美國氣象學會提供）） Fig.11 (a) Reports of the tornadic storm event in Texas Panhandle on 20 April 2010; (b) The 3 radars of WSD-88D: the KAMA, KFDX and KLBB (in black box of 500×500 km) used for analysis, and the analysis area (internal 200 km×200 km) (courtesy of the American Meteorology Society)

圖12 2010年4月20日分析得到的得克薩斯州刀把區域在地面以上3 km的多部雷達合成的反射率（彩色陰影區）、水平風場（風矢量）和垂直渦度（黑色等值線） （a）19：35，（b）19：40，（c）19：45，（d）19：50，（e）19：55，（f）20：00（以上均為當地時間，紫色實線表明圖13中垂直剖面圖的位置；圖12a中小框為圖14的分析區域（美國氣象學會提供）） Fig.12 The synthesized reflectivity (color shaded, dBz), horizontal wind fields (vectors, m/s), and vertical vorticity (contour interval, 2×10-3 s-1) at 3 km AGL observed by 88D radars shown in Fig.6b, at (a) 19:35, (b) 19:40, (c) 19:45, (d) 19:50, (e) 19:55, and (f) 20:00 CST, 20 April 2010 near Umbarger, Texas. Maroon line denotes location of cross-sections in Fig.13. Black box in Fig.12a is the zoom-in area for Fig.14 (courtesy of the American Meteorology Society)

圖13 圖12中的兩個時次的垂直剖面圖 （a）和（b）為垂直渦度（黑色等值線）疊加于反射率； （c）和（d）為風矢量疊加于反射率（美國氣象學會提供） Fig.13 Vertical vorticity (black contours in a and b), wind vectors (in c and d), and reflectivity (color shaded) along a vertical slice indicated in Fig.12, at 19:45 CST (left panel), and 19:50 CST (right panel), on 20 April 2010 near Umbarger, Texas. Between these two time levels, a tornado touched down (courtesy of the American Meteorology Society)

但三維變分系統在理論上不是最優的方法，存在著很多理論上的缺陷。例如不能引入觀測資料的時間分布信息，背景誤差協方差距陣不能隨天氣形勢演變等。由于這些缺陷，不少對流尺度的資料同化和預報的研究工作中也探索了運用集合卡爾曼濾波和四維變分同化的方法。但是這些方法各有優缺點[67-68]。綜合考慮各種因素，例如現階段對流尺度，或對流尺度數值模式的準確性，預報又必須及時發布（時效性），仍然認為用改進的三維變分同化方法對高分辨率數值預報提供初值仍然是一個方便、高效和實用的選擇。首先可以加入一些平衡約束來使分析變量之間保持平衡，同時也可把集合預報的信息引入到三維變分系統用來改進系統的誤差協方差距陣，使之隨流型的改變而變。Lorenc[67]首先提出了將變分同化和集合卡爾曼濾波相結合的方法。在這個方法中，變分同化系統可以是三維變分同化，也可以是四維變分同化。

在最近幾年出現了大量的混和同化方面的研究，大部分研究集中在對大尺度和中尺度天氣的資料同化及數值預報方面[9, 69-75]。對對流尺度的研究相對較少。Gao等[76]在他們的研究中證明混合三維變分和集合卡爾曼濾波方法在對流尺度的數值預報也很有用。把從集合預報中估計的背景誤差協方差距陣用在三維變分同化中可以有效地提高雷達資料同化的精度，甚至在集合數量不太大的情況下也是如此。用理想數據的數值試驗表明這個方法優于三維變分的同化方法。但在這個方法中必須同時使用三維變分和集合卡爾曼濾波。因為雷達資料分辨率高，資料量巨大，使得計算量巨大，特別是對高分辨率的對流尺度的數值模式。為減少計算量，同時利用變分和集合卡爾曼方法的優點，Gao等[77]提出了采用集合變分同化的方法把雷達資料同化到對流尺度的數值模式中。集合變分同化的概念和方法首先在法國氣象局和歐州中期數值預報中心使用[78-79]。在這里Gao等[77]利用增加alpha控制變量的方法[67]來計算集合誤差協方差距陣。如果集合變分資料同化的核心框架用的是四維變分同化，那么這個方法計算量也會很大。但是Gao等采用的是計算效率較高的三維變分同化方法，同時這種方法不需運行集合卡爾曼濾波，所以此方法在集合數量不是太大的情況下計算代價較小，容易并行化。這種方法的優點是能夠應用從集合預報中導出隨時間和流型而變的集合誤差背景場協方差距陣，所以較適用于對流尺度的雷達資料同化。現在作者正在集中精力調試這些方案，希望將來能有效地應用到NOAA的災害性天氣試驗基地的短期天氣預報試驗中，能對基于這些預報的強風暴預警提供支撐。

圖14 圖例說明同圖12，但分析區域為圖12a中的被放達的黑邊小框區域（美國氣象學會提供） Fig.14 As same as in Fig 12, but the analysis area is for the zoom-in area in Fig.12a (courtesy of the American Meteorology Society)

說明和致謝：本文有幾位作者在過去幾年中在俄克拉荷馬大學作訪問學者時參與了此軟件的研究與開發，包括湖南省氣象臺傅承浩、唐杰、周長青，以及國家氣象中心莊照榮等。在研究過程中曾得到俄克拉荷馬大學超級計算中心Henry Neeman博士和他的團隊在計算機資源及應用方面的大力幫助。在此一并表示感謝！

圖15 2010年5月6日湖南省自動雨量站逐小時累積雨量圖 Fig.15 Precipitation observed by the surface rain gauge

圖16 2010年5月6日湖南暴雨過程的分析區域和雷達位置（內黑框為分析區域）Fig.16 Positions of two radars (9731, 9736) and the analysis domain for Hunan heavy rain event on 6 May 2010 (The internal domain is the analysis domain in Hunan Province)

圖17 2010年5月6日同化得到的3 km高度上的雷達回波和風場疊加 （a）03：00，（b）04：00，（c）05：00，（d）06：00，（e）07：00，（f）08：00 Fig.17 The analyzed u-v wind vector and reflectivity (colored) at z=3 km AGL for a local heavy rain and severe convective weather event occurred at central south of Hunan Province in China: at 03:00 (a), 04:00 (b), 05:00 (c), 06:00 (d), 07:00 (e), 08:00 (f) BT on 6 May 2010

[1] Lewis J, Derber J. The use of adjoint equations to solve a variational adjustment problem with advective constraints. Tellus, 1985, 37A: 307-322.

[2] Le Dimet F, Talagrand O. Variational algorithms for analysis and assimilation of meteorological observations: theoretical aspects. Tellus, 1986, 38: 97-110.

[3] Talagrand O, Courtier P. Variational assimilation of meteorological observations with adjoint vorticity equation. I: theory. Quarterly Journal of the Royal Meteorological Society, 1987, 113: 1311-1328.

[4] Thacker C, Long R. Fitting dynamics to data. Journal of Geophysical Research, 1988, 93: 1227-1240.

[5] Daley R. Atmospheric data analysis. Cambridge: Cambridge University Press, 1991: 457.

[6] Sun J, Crook N A. Dynamical and microphysical retrieval from Doppler radar observations using a cloud model and its adjoint. Part I: model development and simulated data experiments. Journal of the Atmospheric Sciences, 1997, 54: 1642-1661.

[7] Rabier F, J?rvinen H, et al. The ECMWF operational implementation of 4D variational assimilation. Part I: experimental results with simplified physics. Quarterly Journal of the Royal Meteorological Society, 2000, 126: 1143-1170.

[8] Barker D, Huang W, Guo Y, et al. A three-dimensional (3DVAR) data assimilation system for use with MM5: implementation and initial results. Monthly Weather Review, 2004, 132: 897-914.

[9] Barker D, Huang X, et al. The Weather research and forecasting model’s community variational/ensemble data assimilation system: WRFDA. Bulletin of the American Meteorological Society, 2012: 93, 831-843.

[10] Xiao Q, Kuo Y H, Sun J, et al. Assimilation of Doppler radar observations with a regional 3DVAR system: impact of Doppler velocities on forecasts of a heavy rainfall case. Journal of Applied Meteorology, 2005, 44: 768-788.

[11] Lewis J, Lakshmivarahan S, Dhall S. Dynamic data assimilation: a least squares approach. Cambridge University Press, 2006: 654.

[12] Sun J, Flicker D, Lilly D. Recovery of three-dimensional wind and temperature fields from simulated Doppler radar data. Journal of the Atmospheric Sciences, 1991, 48: 876-890.

[13] Xu Q, Qiu C. Simple adjoint methods for single-Doppler wind analysis with a strong constraint of mass conservation. Journal of Atmospheric and Oceanic Technology, 1994, 11: 289-298.

[14] Xu Q, Qiu C. Adjoint-method retrievals of low-altitude wind fields from single-Doppler reflectivity and radial-wind. Journal of Atmospheric and Oceanic Technology, 1995, 12: 1111-1119.

[15] Xu Q, Gu H, Yang S. Simple adjoint method for three-dimensional wind retrievals from single-Doppler radar. Quarterly Journal of the Royal Meteorological Society, 2001, 127: 1053-1067.

[16] Qiu C, Xu Q. Least Squares Retrieval of Microburst Winds from Single-Doppler radar data. Monthly Weather Review, 1996, 124: 1132-1144.

[17] Sun, J, Crook N A. Real-time low-level wind and temperature analysis using single WSR-88D data. Weather and Forecasting, 2001, 16: 117-132.

[18] Gao J, Xue M, et al. A variational method for the analysis of threedimensional wind fields from two Doppler radars. Monthly Weather Review, 1999, 127: 2128-2142.

[19] Gao J, Xue M, et al. New development of a 3DVAR system for a nonhydrostatic NWP model. Preprint, 15th Conf. Numerical Weather Prediction and 19th Conf. on Weather Analysis and Forecasting, San Antonio, TX, American Meteorological Society, 2002, 339-341.

[20] Gao J, Xue M, et al. A three-dimensional variational data assimilation method with recursive filter for single-Doppler radar, Journal of Atmospheric and Oceanic Technology, 2004, 21: 457-469.

[21] Purser R J, Wu W, et al. Numerical aspects of the application of recursive filters to variational statistical analysis. Part I: Spatially homogeneous and isotropic Gaussian covariances. Monthly Weather Review, 2003, 131: 1524-1535.

[22] Wu W, Purser R J, Parrish D. Three-dimensional variational analysis with spatially inhomogeneous covariances. Monthly Weather Review, 2002, 130: 2905-2916.

[23] Hu M, Xue M, Brewster K. 3DVAR and cloud analysis with WSR-88D Level-II data for the prediction of the fort worth tornadic thunderstorms. Part I: cloud analysis and its impact. Monthly Weather Review, 2006, 134: 675-698.

[24] Hu M, Xue M, Gao J, et al. 3DVAR and cloud analysis with WSR-88D Level-II data for the prediction of the fort worth tornadic thunderstorms. Part II: impact of radial velocity analysis via 3DVAR. Monthly Weather Review, 2006, 134: 699-721.

[25] Gao J, Stensrud D, et al. Three-dimensional analyses of several thunderstorms observed during VORTEX2 field operations. 34th Conference on Radar Meteorology, Willimsburg, VA, 2009, Online publication.

[26] Schenkman A, Xue M, et al.The analysis and prediction of the 8-9 May 2007 Oklahoma tornadic mesoscale convective system by assimilating WSR-88D and CASA radar data using 3DVAR. Monthly Weather Review, 2011, 139: 224-246.

[27] Albers S, McGinley J, et al. The local analysis and prediction system (LAPS): analysis of clouds, precipitation and temperature. Weather and Forecasting, 1996, 11: 273-287.

[28] Zhang J, Carr F, Brewster K. ADAS cloud analysis. Preprints, 12th Conf on Numerical Weather Prediction, Phoenix, AZ, American Meteorological, Society, 1998: 185-188.

[29] Xue M, Droegemeier K K, Wong V. The advanced regional prediction system (ARPS): a multiscale nonhydrostatic atmospheric simulation and prediction tool. Part I: model dynamics and verification. Meteorology and Atmospheric Physics, 2000, 75: 161-193.

[30] Xue M, Wang D, et al. The advanced regional prediction system (ARPS), storm scale numerical weather prediction and data assimilation. Meteorology and Atmospheric Physics, 2003, 82: 139-170.

[31] Gao J, Ge G, et al. The relative importance of assimilating radial velocity and reflectivity data to storm-scale analysis and forecast. 23nd Conf Wea Anal Forecasting/19th Conf Num Wea Pred, Omaha, NB, American Meteorology Society, 2009, Paper 16A.3.

[32] Snyder C, Zhang F. Assimilation of simulated Doppler radar observations with an ensemble Kalman filter. Monthly Weather Review, 2003, 131: 1663-1677.

[33] Zhang F, Snyder C, Sun J. Impacts of initial estimate and observations on the convective-scale data assimilation with an ensemble Kalman filter. Monthly Weather Review, 2004, 132: 1238-1253.

[34] Tong M, Xue M. Ensemble Kalman filter assimilation of Doppler radar data with a compressible nonhydrostatic model: OSS experiments. Monthly Weather Review, 2005, 133: 1789-1807.

[35] Meng Z, Zhang F. Tests of an ensemble Kalman filter for mesoscale and regional-scale data assimilation. Part II: imperfect model experiments. Monthly Weather Review, 2007, 135: 1403-1423.

[36] Meng Z, Zhang F. Limited-area ensemble-based data assimilation. Monthly Weather Review, 2011, 139: 2025-2045.

[37] Gao J, Xue M, An efficient dual-resolution approach for ensemble data assimilation and tests with assimilated Doppler radar data. Monthly Weather Review, 2008, 136: 945-963.

[38] Dowell, D, Wicker L, et al. Ensemble Kalman filter assimilation of radar observations of the 8 May 2003 Oklahoma City supercell: influences of reflectivity observations on storm-scale analyses. Monthly Weather Review, 2011, 132: 1982-2005.

[39] Yussouf N, Mansell E R, et al. The ensemble Kalman filter analyses and forecasts of the 8 May 2003 Oklahoma city tornadic supercell storm using single-and double-moment microphysics schemes. Monthly Weather Review, 2013, 141: 3388-3412.

[40] Gao J, Stensrud D. Assimilation of reflectivity data in a convectivescale, cycled 3DVAR framework with hydrometeor classification. Journal of the Atmospheric Sciences, 2012, 69: 1054-1065.

[41] Gao J, Smith T, et al. A realtime weather-adaptive 3DVAR analysis system for severe weather detections and warnings with automatic storm positioning capability. Weather and Forecasting, 2013, 28: 727-745.

[42] Stensrud D, Gao J. Importance of horizontally inhomogeneous environmental initial conditions to ensemble storm-scale radar data assimilation and very short range forecasts. Monthly Weather Review, 2010, 138: 1250-1272.

[43] Ge G, Gao J, et al. Impacts of beam broadening and earth curvature in radar data assimilation. Journal of Atmospheric and Oceanic Technology, 2010, 27: 617-636.

[44] Ge G, Gao J, Xue M. Diagnostic pressure equation as a weak constraint in a storm-scale three dimensional variational radar data assimilation system. Journal of Atmospheric and Oceanic Technology, 2012, 29: 1075-1092.

[45] Doviak, R, Zrnic D. Doppler Radar and Weather Observations. 2nd ed. Academic Press, 1993: 562.

[46] Rihan F A, Collier C G, et al. Assimilation of Doppler radial winds into a 3D-Var system: errors and impact of radial velocities on the variational analysis and model forecasts. Quarterly Journal of the Royal Meteorological Society, 2008, 134: 1701-1716.

[47] Wood V T, Brown R A. Effects of radar sampling on single-Doppler velocity signatures of mesocylcones and tornadoes. Weather and Forecasting, 1997, 12: 928-938.

[48] Doviak, R, Ray P, et al. Error estimation in wind fields derived from dual-Doppler radar measurement. Journal of Applied Meteorology, 1976, 15: 868-878.

[49] Miller L J, Sun J. Initialization and forecasting of thunderstorm: specification of radar measurement errors. Preprints, 31st Conf. on Radar Meteorology, Seattle, WA, American Meteorological Society, 2003: 146-149.

[50] Lin Y, Farley R D, Orville H D. Bulk parameterization of the snow field in a cloud model. Journal of Applied Meteorology, 1983, 22: 1065-1092.

[51] Gilmore M, Straka J, Rasmussen E. Precipitation and evolution sensitivity in simulated deep convective storms: comparisons between liquid-only and simple ice and liquid phase microphysics. Monthly Weather Review, 2004, 132: 1897-1916.

[52] Smith P, Myers C, Orville H. Radar reflectivity factor calculations in numerical cloud models using bulk parameterization of precipitation processes. Journal of Applied Meteorology, 1975, 14: 1156-1165.

[53] Sasaki Y. Some basic formalisms in numerical variational analysis. Monthly Weather Review, 1970, 98: 875-883.

[54] Sasaki Y, Numerical variational analysis formulated under the constraints as determined by longwave equations and a lowpass filter. Monthly Weather Review, 1970, 98: 884-898.

[55] Sasaki Y. Numerical variational analysis with weak constraint and application to surface analysis of severe storm gust. Monthly Weather Review, 1970, 98: 899-910.

[56] Courtier P, Thépaut J, et al. A strategy for operational implementation of 4D-Var, using an incremental approach. Quarterly Journal of the Royal Meteorological Society, 1994, 120: 1367-1388.

[57] Courtier P. Dual formulation of four dimensional variational assimilation. Quarterly Journal of the Royal Meteorological Society, 1997, 123: 2449-2461.

[58] Lorenc A. Iterative analysis using covariance functions and filters. Quarterly Journal of the Royal Meteorological Society, 1992, 118: 569-591.

[59] Xie Y, Koch S E, et al. A sequential variational analysis approach for mesoscale data assimilation. Preprints, 21st Conf. on Weather Analysis and Forecasting/17th Conf on Numerical Weather Prediction, Washington, DC, American Meteorological Society, 2005, 15B.7. http://ams.confex.com/ams/pdfpapers/93468.pdf.

[60] Xie Y, Koch S E, et al. A Space and time multiscale analysis system: a sequential variational analysis approach. Monthly Weather Review, 2011, 139: 1224-1240.

[61] Smith T, Gao J, et al. Performance of a real-time 3DVAR analysis system in the hazardous weather testbed. Weather and Forecasting, 2014, 29: 63-77.

[62] Calhoun K, Smith T, et al. Forecaster use and evaluation of realtime 3DVAR analyses during severe thunderstorm and tornado warning operations in the hazardous weather testbed. Weather and Forecasting, 2014, 29: 601-613.

[63] Gao J, Brewster K, et al, Real-time, low-level wind analysis including CASA and WSR-88D radar data using the ARPS 3DVAR. 25th Conf Severe Local Storms, American Meteorological Society, 2010: 7B.4 (online publication).

[64] Xue M, Kong F, et al. CAPS realtime storm scale ensemble and high resolution forecasts for the NOAA hazardous weather testbed 2010 spring experiment. 25th Conf Severe Local Storms, American Meteorological Society, 2010: 7B.3.

[65] Brewster K, Thomas K, et al. A nowcasting system using full physics numerical weather prediction initialized with CASA and NEXRAD radar data. Preprints, 25th Conf Severe Local Storms, Denver, CO, American Meteorological Society, Denver, CO, 2010: 9.4.

[66] Kain J S, Xue M, Coniglio M C, et al. Assessing advances in the assimilation of radar data within a collaborative forecastingresearch environment. Weather and Forecasting, 2010, 25: 1510-1521.

[67] Lorenc A. The potential of the ensemble Kalman filter for NWP: a comparison with 4DVar. Quarterly Journal of the Royal Meteorological Society, 2003, 129: 3183-3204.

[68] Kalnay E, Li H, et al. 4-D-Var or ensemble Kalman filter? Tellus, 2007, 59: 758-773.

[69] Buehner M. Ensemble-derived stationary and flow-dependent background-error covariances: evaluation in a quasi-operational NWP setting. Quarterly Journal of the Royal Meteorological Society, 2005,131: 1013-1043.

[70] Buehner M, Houtekamer P, et al. Intercomparison of variational data assimilation and the ensemble Kalman filter for global deterministic NWP. Part I: description and single-observation experiments. Monthly Weather Review, 2010, 138: 1550-1566.

[71] Buehner M, Houtekamer P, et al. Intercomparison of variational data assimilation and the ensemble Kalman filter for global deterministic NWP. Part II: one-month experiments with real observations. Monthly Weather Review, 2010, 138: 1567-1586.

[72] Wang X, Barker D M, et al. A Hybrid ETKF–3DVAR data assimilation scheme for the WRF model. Part I: observing system simulation experiment. Monthly Weather Review, 2008, 136: 5116-5131.

[73] Wang X, Barker D M, et al. A Hybrid ETKF–3DVAR data assimilation scheme for the WRF model. Part II: real observation experiments. Monthly Weather Review, 2008, 136: 5116-5131.

[74] Wang X, Parrish D, et al. GSI 3DVar-based ensemble-variational hybrid data assimilation for NCEP global forecast system: singleresolution experiments. Monthly Weather Review, 2013, 141: 4098-4117.

[75] Zhang F, Zhang M, Poterjoy J. E3DVar: coupling an ensemble Kalman filter with three-dimensional variational data assimilation in a limited-area weather prediction model and comparison to E4DVar. Monthly Weather Review, 2013, 141: 900-917.

[76] Gao J, Stensrud D. Some observing system simulation experiments with a hybrid 3DEnVAR system for storm scale radar data assimilation. Monthly Weather Review, 2014, 142: 3326-3346.

[77] Gao J, Fu C, et al. OSSEs for an Ensemble 3DVAR data assimilation system with radar observations of convective storms. Journal of the Atmospheric Sciences, 2016, 73: 2403-2426.

[78] Bonavita M, Isaksen L, Holm E. On the use of EDA background error variances in theECMWF 4D-Var. Quarterly Journal of the Royal Meteorological Society, 2012, 138: 1540-1559.

[79] Berre L, Pannekoucke O, Desroziers G, et al. A variational assimilation ensemble and the spatial filtering of its error covariances: increase of sample size by local spatial averaging. Proceedings of the ECMWF Workshop on flow-dependent aspects of data assimilation ECMWF, 2007: 151-168.