基于新型觀測算子的雙偏振雷達雨滴譜變分反演

陳垚 ，寇蕾蕾 ，蔣銀豐，楊春生，林正健，楚志剛

(1. 南京信息工程大學大氣物理學院，江蘇 南京 210044；2. 福建省平潭綜合實驗區氣象局，福建 平潭 350400；3. 南京信息工程大學氣象災害預警與評估協同創新中心，江蘇 南京 210044；4. 平潭海洋氣象野外科學觀測研究站，福建 平潭 350400；5. 福建省災害天氣重點實驗室，福建 福州 350001）

1 引 言

雨是最常見的降水形式。雨滴譜(Raindrop Size Distribution，DSD)是描述降雨的一個重要參量，它與雨滴形狀、取向以及下落末速度一起，能夠提供較完整的降雨信息，并可以計算雨滴散射振幅、雙偏振雷達參量和降雨參數。降雨參數主要包括降雨率(R)、液態水含量(LWC)、質量權重平均直徑(Dm)、體積中值直徑(D0)、歸一化截距參數(NW)和總數密度(Nt)等，可用于雷達定量降水估測、全球水循環分析和數值預報雷達資料同化等研究。通過分析不同類型降水的降雨參數，可更加了解降水微物理過程，改善降水微物理參數化，因此雨滴譜和降雨參數的精確反演對準確監測和預報降水十分重要。

雙偏振雷達通過發射垂直和水平方向的電磁波來實現探測，可獲得兩個通道的回波信息。它不僅能獲得常規多普勒雷達的參數如水平反射率因子(ZH)、徑向速度(v)和譜寬(SW)，還可得到雙偏振觀測量如差分反射率因子(ZDR)、比差分傳播相移(KDP)、總差分相移(ΦDP)、協相關系數(ρhv)等，這些觀測量能夠反映降水系統內部水凝物大小、形狀、相態、濃度等信息，所以雙偏振天氣雷達的觀測量可用于雨滴譜反演，進而計算降雨參數，有助于更好地研究我國的降水微物理特性并定量估測降水[1-4]。

目前，雨滴譜反演主要有兩類方法。第一類是常規的雨滴譜反演，主要通過先建立雨滴譜模型，再推導雙偏振參量和雨滴譜參數的關系式，實現雨滴譜反演。Seliga 等[5]在研究中建立了指數雨滴譜模型，并由ZH、ZDR和KDP確定模型參數后得到雨滴譜。Ulbrich 等[6]提出了gamma 分布模型，并將雙偏振雷達參量和降雨特性聯系在一起。Zhang 等[7]提出利用μ-Λ關系約束gamma 分布，得到兩參數模型(constrained-gamma model，簡稱CG模型)[8-12]，由ZH和ZDR計算C-G模型參數，成功反演雨滴譜。國內也有許多學者利用常規方法反演雨滴譜，如張鴻發等[13]通過統計分析得到ZH、ZDR與雨滴譜參數和雨強的關系，用于反演雨滴譜。李宗飛等[14]通過對雨滴模型的散射模擬以及gamma 分布的擬合，建立雨滴譜參數與雙偏振參量之間的關系反演雨滴譜。Liu 等[15]用廣東惠州龍門站滴譜儀數據計算μ-Λ關系，并基于適用于華南地區的C-G 模型，用ZH和ZDR反演雨滴譜。第二類方法是基于最優化理論的雨滴譜反演，主要將最優化理論和常規的雨滴譜反演方法相結合，或通過建立觀測算子，基于最優化算法求解包含觀測誤差在內的非線性方程組反演雨滴譜。Vulpiani等[16]基于神經網絡算法由S波段雙偏振雷達數據計算雨滴譜。Cao 等[17-18]先利用貝葉斯理論反演雨滴譜，用雨滴譜的先驗估計降低觀測誤差，之后又基于C-G 模型建立觀測算子，利用二維變分理論反演雨滴譜。Yoshikawa 等[19]基于極大似然法由X 波段雷達雙偏振觀測量反演雨滴譜。Wen 等[20]提出了一個由雙偏振觀測量反演雨滴譜的模型，首先基于C-G 模型建立觀測算子，其次構造一個偽訓練集，利用最鄰近算法反演雨滴譜。Vivek 等[21]和Huang 等[22]分別基于指數模型和C-G模型建立觀測算子，利用S波段和C波段雷達數據變分反演得到最優雨滴譜。與觀測算子同理，Sun等[23]建立了雨滴譜模型參數和雷達觀測量之間的逆映射表，由X波段偏振雷達數據反演雨滴譜。

無論是常規反演還是變分反演，都需要在反演前建立雨滴譜模型，利用模型參數表示雨滴譜。但在實際情況中，雨滴譜是復雜且多變的，某一種雨滴譜模型并不能代表所有的雨滴譜分布。前人在變分反演中所用的觀測算子多為建立雨滴譜模型后模擬得到，本研究不建立雨滴譜模型，使用滴譜儀實測數據計算更能反映當地降水特點的觀測算子，將其應用于雙偏振雷達雨滴譜最優化反演中，進一步提高了雨滴譜反演的準確性。

為反演得到更精確的雨滴譜，本文提出了一種基于新型觀測算子的雙偏振雷達變分反演雨滴譜方法。首先，介紹本文所用到的C、S波段雙偏振多普勒天氣雷達數據和Parsivel滴譜儀數據；其次，描述基于新型觀測算子的變分反演方法，并闡述觀測算子的建立；隨后，闡述基于新型觀測算子變分反演雨滴譜的詳細過程；最后，利用一個理想模擬試驗和NUIST-CDP觀測到的實測個例驗證該方法的可行性，并在文末分析了結果優化的原因。

2 數據介紹

2.1 C波段雙偏振多普勒天氣雷達

本研究基于南京信息工程大學C 波段雙偏振多普勒天氣雷達(簡稱NUIST-CDP)數據最優反演雨滴譜。雷達站點坐標為118.717 2 °E，32.206 9 °N，天線直徑11 m，掃描半徑150 km，脈沖寬度0.5 μs，距離庫長75 m，波束寬度0.54 °，完成一次體掃所需時間為7～8 min，可探測得到水平反射率因子(ZH)、徑向速度(v)、速度譜寬(SW)、差分反射率因子(ZDR)、線性退偏振比(LDR)、差分傳播相移(φDP)、差分傳播相移率(KDP)、協相關系數(ρhv)等參量。本研究在反演之前對雷達數據進行如下預處理：基于地物雜波識別算法[24]去除雷達數據中的地物雜波。利用21 點中值濾波器平滑ZH、ZDR、φDP，消除原始數據中出現的異常跳變值。對φDP進行退模糊，當相鄰兩點的φDP值之差大于320 °，若前一點大于0，則將后一點加上360 °，反之減360 °，模糊處理后的φDP再用21 點中值濾波器平滑，最后用最小二乘法計算KDP。因為信噪比和ρhv較小時雷達數據質量較差[25]，所以在實際應用中剔除了信噪比小于20 dB、ρhv小于0.9、ZDR和KDP為負值的雷達數據。本研究在使用常規方法反演雨滴譜時，對雷達數據進行了如下衰減訂正[26]：當0.1 °/km≤KDP≤3 °/km 時，使用AH和KDP的關系進行訂正(該關系由滴譜儀數據擬合得到，下同)，訂正公式為AH=0.092×KDP；當KDP＜0.1 °/km或KDP＞3 °/km 時，使用AH和ZH的關系訂正，訂正公式為AH=1.436×10-5×ZH0.855(ZH的單位為mm6/m3)。

2.2 S波段多普勒天氣雷達

由于本研究的最優雨滴譜反演可同時實現衰減訂正，南京龍王山S 波段多普勒天氣雷達(簡稱SA 雷達)反射率因子PPI(Plane Position Indicator)圖像主要用于檢驗衰減訂正的結果。SA 雷達站點坐標為118.696 9 °E，32.190 8 °N，距離庫長1 km，波束寬度約為1 °，完成一次體掃所需時間為6 min，其距離NUIST-CDP 約2.5 km，二者高度相差約50 m，對比時選取二者體掃起始時刻最接近的數據。

2.3 OTT Parsivel滴譜儀

本研究使用滴譜儀數據計算結果來建立觀測算子、模擬雷達徑向廓線和檢驗反演結果，數據來自于三部德國OTT 公司生產的Parsivel 滴譜儀。第一部滴譜儀位于南京信息工程大學大氣綜合觀測基地，站點坐標為118 °42 ＇E，32 °12 ＇N，距離NUIST-CDP 約1 800 m，相對雷達的方位角為250 °。第二部滴譜儀位于高淳站點，坐標為118.903 9 °E，31.333 3 °N，距離NUIST-CDP 約98km，相對雷達的方位角為170 °。第三部滴譜儀位于江寧站點，坐標為118.899 7 °E，31.931 4 °N，距離NUIST-CDP 約35 km，相對雷達的方位角為151 °。在反演之前對滴譜儀數據進行如下預處理：用6 點中值濾波器平滑數據，當觀測通道的雨滴個數小于10 個或計算出的降雨率小于0.5 mm/h，則剔除該點數據。本研究從滴譜儀數據計算雙偏振雷達觀測量和降雨參數，先從滴譜儀數據中讀取出滴譜信息即N(D)，再利用T 矩陣方法[27]求得C 波段粒子散射系數，進而計算雙偏振雷達觀測量和降雨參數。

2.4 雨量計

本文所用的雨量計數據來源于南京市地面雨量自動站，共計85 個雨量站點。雨量計數據時間分辨率為1 小時，雨量分辨率為0.1 mm。本文以雨量計數據為實際雨強來驗證最優雨滴譜計算降雨率的準確性。NUIST-CDP、SA 雷達、三部滴譜儀以及雨量計站點的空間分布如圖1所示。

圖1 NUIST-CDP、SA雷達、三部滴譜儀和雨量計站點分布圖

3 雨滴譜反演方法

3.1 代價函數

基于變分理論的雙偏振雷達雨滴譜反演本質上是求狀態變量的先驗估計、背景誤差、觀測誤差、雙偏振觀測量和觀測算子預測結果組成的非線性方程組最優解的過程。因為雨滴譜參數和Dm、LWC之間可相互轉化，且Dm和LWC符合高斯分布[21]，本研究使用Dm和LWC 作為狀態變量，建立雙偏振觀測量和雨滴譜之間的關系。將變分反演雨滴譜的代價函數定義為：

式(1)中，x 為每條雷達徑向上n個距離庫Dm和LWC的向量形式：

式(2)～(4)中，Dm(n)和LWC(n)代表第n個距離庫上的狀態變量，當狀態變量確定時，可用觀測算子預測雙偏振觀測量。xb為狀態變量的先驗估計(背景場)，由狀態變量和雙偏振參量之間的經驗關系計算得到。B 為背景場的誤差協方差矩陣，通過對狀態變量進行南京地區氣候統計分析后，將背景誤 差 標準差σDm和σLWC分 別 設置為1 mm 和0.707 g/cm3(方差分別為1 mm2和0.5 g2/cm6)，本研究假設了一個高斯相關模型來計算協方差矩陣[28]：

式(5)中，σb2是背景誤差協方差(σDm2和σLWC2)，rij是第i個和第j個雷達距離庫之間的距離，rL是空間去相關長度，會影響最后解的平滑程度，通過多次模擬測試后設置為1 000 m。y為上文提到的雙偏振觀測量，表示為向量形式：

式(6)～(9)中，分號代表矩陣的縱向合并，n表示距離庫序號，Z'H和Z'DR表示衰減訂正后的ZH和ZDR。考慮衰減效應，當由Dm和LWC得到雨滴譜并計算AH和ADP之后，各距離庫上ZH和ZDR可表示為：

式(10)、(11)中，δr代表雷達徑向上每個距離庫的庫長度(即雷達徑向分辨率)。H(x)為觀測算子預測得到的雙偏振參量，觀測算子的建立在后文中會詳細說明。R 為觀測誤差協方差矩陣，由三個雙偏振觀測量的觀測誤差標準差組成。對雙偏振雷達觀測量進行統計分析后，將ZH、ZDR和KDP的觀測誤差標準差分別設置為1 dB、0.2 dB 和0.6 °/km，假定ZH、ZDR和KDP的誤差相互獨立時，

式(12)中，O 代表零矩陣，RZH、RZDR和RKDP分別對應于三個雙偏振參量的觀測誤差協方差矩陣。由于雷達觀測量的隨機誤差在各距離庫上相互獨立，因此RZH、RZDR和RKDP副對角線上的元素為0。

因為觀測算子是非線性的，所以代價函數最小值的求解是一個非線性最優化問題。該類問題(例如變分資料同化問題)通常用高斯-牛頓迭代方法求解，該方法的求解公式為[29]：

式(13)、(14)中，W 是最優化權重矩陣，H 是觀測算子H的雅可比式，即觀測算子的線性近似，由每個觀測算子分別對兩個狀態變量求偏導數得到。

3.2 觀測算子的建立

降雨參數和雙偏振參量都由雨滴譜決定，因此降雨參數和雙偏振參量之間可定量關聯，這也是由雙偏振參量反演雨滴譜并計算降雨參數的基本原理。將降雨參數中的Dm、LWC(即狀態變量)和雙偏振參數之間的函數關系定義為觀測算子。

本研究提出了一種由滴譜儀實測數據計算的新型觀測算子，首先選取位于高淳站點的2015 年和2016 年7 月、8 月滴譜儀觀測數據得到雨滴譜，假設入射波波段為C，波長為5 cm，溫度為20 ℃，再利用T 矩陣方法計算粒子的C 波段散射系數，由散射系數和雨滴譜估計雙偏振參量ZH、ZDR、KDP、AH和ADP，具體計算公式如下：

式(15)中，λ是電磁波波長，D是雨滴等效直徑，K是介電因子，定義為是雨滴的相對介電常數，Dmax和Dmin分別代表最大和最小雨滴等效直徑，shh,vv(π,D)代表后向散射系數。

式(18)中，Re 代表復數的實部部分，shh(0,D)代表水平入射水平散射的前向散射系數，svv(0,D)代表垂直入射垂直散射的前向散射系數。

式(19)中，Im代表復數的虛部部分。

狀態變量用階矩法計算，可用雨滴譜的積分量表示，一般用Mn來表示雨滴譜的n階矩：

LWC正比于雨滴譜的三階矩：

Dm為四階矩和三階矩之比：

最后將計算出的雙偏振觀測量和狀態變量擬合成曲線，即為新型觀測算子。之前學者使用的觀測算子多為模擬得到，先假設雨滴譜分布模型，設置參數模擬出狀態變量和雨滴譜后，由T 矩陣方法模擬的散射系數計算雙偏振參量，再擬合雙偏振參量和狀態變量的曲線得到模擬觀測算子。將新型觀測算子、模擬觀測算子以及南京信息工程大學綜合觀測基地2015 年和2016 年7 月、8 月的滴譜儀觀測數據計算的雙偏振參量和狀態變量散點對比。圖2a～2e 分別為Zh/LWC、ZDR、KDP/LWC、AH/LWC和ADP/LWC關于Dm的新型觀測算子(紅色實線)、模擬觀測算子(黑色實線)和滴譜儀數據計算值(藍色散點)對比結果，與模擬觀測算子相比，基于實測數據的新型觀測算子和滴譜儀數據計算值有更高的一致性。從圖2a 和圖2b 可看出，在雨滴較大時，新型觀測算子和模擬觀測算子存在較大差異，這可能是因為，雨滴越大越偏離球形，而計算雙偏振參量和狀態變量時，將雨滴看作近似球形，所以模擬觀測算子時，雨滴越大計算出的雙偏振參量和狀態變量越不準確。在以后的研究中，考慮基于非球形粒子，利用T 矩陣方法模擬散射系數，從而計算更準確的雙偏振參量和狀態變量。

圖2 新型觀測算子、模擬觀測算子和滴譜儀數據計算值對比

3.3 反演步驟

基于新型觀測算子由雙偏振雷達觀測量變分反演雨滴譜的流程圖如圖3所示。

圖3 基于新型觀測算子的變分雨滴譜反演流程圖

首先，導入雙偏振雷達數據，對其進行質量控制。再利用狀態變量和雙偏振參量間的經驗關系計算狀態變量的先驗估計。然后，如3.1 節所述，建立背景誤差協方差矩陣和觀測誤差協方差矩陣。隨后，導入預處理過后的滴譜儀數據，利用3.2 節中的公式計算雙偏振參量和狀態變量，推導新型觀測算子并計算觀測算子關于狀態變量的微分，進而構建代價函數，將觀測量、狀態變量的先驗估計和誤差協方差矩陣代入其中，再由新型觀測算子將狀態變量轉化為觀測量。最后利用高斯-牛頓迭代法求解代價函數，當結果收斂或達到迭代次數后得到最優雨滴譜，結束迭代，否則更新狀態變量，繼續迭代。得到最優雨滴譜后由3.2 節中相關公式計算雙偏振參量和降雨參數。

4 實驗驗證

為驗證上述基于新型觀測算子的變分反演方法的可行性，本節利用一個理想模擬試驗和NUIST-CDP實測個例最優化反演雨滴譜。

4.1 理想模擬試驗

理想模擬試驗的設計如下：滴譜儀數據觀測誤差比較小且不受衰減效應影響，因此將滴譜儀數據計算結果作為理想模擬試驗的真值。首先，假定雷達為C 波段，波長為5 cm，溫度為20 ℃。根據南京信息工程大學觀測基地2015 年5 月8 日的滴譜儀觀測數據，由前文所述方法計算Dm、LWC、ZH、ZDR和KDP，并使用中值濾波和立方插值模擬出尺度約為35 km(500 個雷達距離庫)且徑向分辨率為75 m 的一次天氣過程的徑向廓線，如圖4a～4e 中藍色實線所示，以此作為真值。考慮到電磁波的衰減效應和隨機觀測誤差(ZH、ZDR和KDP的隨機觀測誤差標準差分別為1 dB、0.2 dB 和0.6 °/km)，模擬的觀測值廓線如圖4c～4e 中黑色實線所示。在距離雷達250～400個距離庫的地方出現較強降雨區，ZH和ZDR明顯受到衰減效應的影響，觀測值與真值的差距逐漸增大。

利用滴譜儀數據模擬的雙偏振雷達觀測量根據基于新型觀測算子的變分方法反演，迭代求解得到最優雨滴譜(由狀態變量Dm和LWC 得到)，再由公式(15)～(18)計算得到ZH、ZDR和KDP，結果如圖4a～4e 中紅色虛線所示。反演值和最優化雨滴譜計算結果都十分接近真值。

圖4 Dm(a)、LWC(b)、ZH(c)、ZDR(d)、KDP(e)的理想模擬試驗

4.2 實測個例驗證

除了理想模擬試驗，本研究還利用NUISTCDP觀測到的一次對流云降水天氣個例和一次層狀云降水天氣個例驗證該算法。

第一個個例是2015 年8 月10 日03—14 時(北京時間，下同)經過南京的一次對流云降水天氣過程，選取11:15 的NUIST-CDP PPI 體掃數據(0.5 °仰角，下同)反演。圖5a、圖5b 分別為最優反演得到的狀態變量Dm、LWC，由狀態變量和雨滴譜參數的關系可得到最優雨滴譜，再根據3.2 節中的公式便可計算Nt、R(圖5c、圖5d)和雙偏振觀測量(圖6c、圖7b、圖7d)。圖6a、圖6b、圖6d分別為2015年8 月10 日11:15 的0.5 °仰角NUIST-CDPZHPPI圖、常規反演雨滴譜計算的ZH和SA 雷達ZHPPI圖。圖7a、圖7c 分別為2015 年8 月10 日11:15 的0.5 °仰角NUIST-CDPZDR、KDPPPI 圖。由圖6a 可看出，此次天氣過程造成了大范圍降雨，雷達北部和東南部分布著塊狀回波單體，結構緊密，邊緣清晰，ZH大于40 dBZ。總體上，Dm與反演雨滴譜計算出的ZDR具有相似的空間分布形式，Nt、LWC、R和最優化雨滴譜計算的ZH也具有相似的空間分布形式。當雷達波束穿過降雨區時，ZH和ZDR衰減明顯。一般認為S 波段雷達信號不受衰減效應的影響[30]，因此可利用起始時刻最接近的SA 雷達反射率因子PPI 圖像(0.5 °仰角，下同)作為真值檢驗反演結果的好壞。將圖6a、圖6c 和圖6d 對比可得，受降雨區影響被衰減的回波強度得到了訂正，并且接近真實值。對比圖6b、圖6c和圖6d三張圖可發現，最優化雨滴譜計算的ZH比常規方法反演雨滴譜計算的ZH更接近真值，整體上更連續，在距離雷達較遠處回波更完整。原因是變分反演方法中AH和ADP直接由雨滴譜計算得到，在降低雨滴譜不確定性對衰減訂正影響的同時，避免了線性關系衰減訂正方法中誤差的多步傳播。再對比圖7a和圖7b、圖7c 和圖7d，會發現距離雷達較遠處原本為負值的ZDR被訂正為正值，由于穿過降雨區受衰減效應影響的ZDR值也得到了修正，最優化雨滴譜計算的KDP與原始值具有幾乎相同的空間分布和強度。

圖5 基于新型觀測算子的變分方法2015年8月10日11:15的0.5 °仰角NUIST-CDP數據反演雨滴譜計算的Dm(a)、LWC(b)、Nt(c)、R(d) Nt以對數尺度顯示。

圖6 2015年8月10日11:15的0.5 °仰角NUIST-CDP ZH PPI圖(a)、常規反演雨滴譜計算的ZH(b)、最優化雨滴譜計算的ZH(c)、SA雷達ZH PPI圖(d)

圖7 2015年8月10日11:15的0.5 °仰角NUIST-CDP ZDR PPI圖(a)、最優化雨滴譜計算的ZDR(b)、NUIST-CDP KDP PPI圖(c)、最優化雨滴譜計算的KDP(d)

本研究提出的最優雨滴譜反演方法能得到很好的滴譜反演效果，但從反演結果可看出，圖像并不完整，這是因為高斯-牛頓迭代方法要求進行迭代的向量必須是連續的，因此本研究在變分反演時對每條雷達徑向數據只選取第一段連續的部分。另外，在實際應用當中，如果有非雨滴的水凝物(例如冰雹)出現，可能會出現沒有解或異常解的情況，這是因為本研究中觀測算子推導的前提是粒子為雨滴，這也是后續研究的主要內容之一。

第二個個例是2016 年7 月2 日19—23 時發生在南京的一次層狀云降水天氣過程，選取20:26 的NUIST-CDP PPI體掃數據反演。圖8a～8d分別為最優反演得到的Dm、LWC、Nt和R。圖9a～9d分別為2016 年7 月2 日20:26 的0.5 °仰角NUIST-CDPZHPPI 圖、常規反演雨滴譜計算的ZH、最優化雨滴譜計算的ZH和SA 雷達ZHPPI 圖。圖10a～10d 分別 為2016 年7 月2 日20:26 的0.5 °仰 角NUISTCDPZDR、KDPPPI 圖和最優化雨滴譜計算的ZDR、KDP。由圖9a 可看出，雷達回波呈片狀均勻分布，回波整體強度不大，雷達南面存在30 dBZ 以上較強回波，最大可達40 dBZ。同樣地，總體上Dm與反演雨滴譜計算出的ZDR具有相似的空間分布形式，Nt、LWC、R和最優化雨滴譜計算的ZH也具有相似的空間分布形式。由于降雨強度不大，衰減效應并不明顯，但依然存在。對比圖9a 與9c，圖10a 與10b，發現雖然ZH和ZDR受衰減效應影響較小，但在反演結果中也得到了修正。對比圖9b、9c和9d，可看出常規方法反演雨滴譜計算的ZH同樣出現不連續的情況，在雷達南側區域，基于新型觀測算子的變分反演雨滴譜計算的ZH更接近真值。對比圖10c 和10d 可得，最優化雨滴譜計算的KDP與原始值圖像具有相似的空間分布形式和數值大小。

圖8 基于新型觀測算子的變分方法2016年7月2日20:26的0.5 °仰角NUIST-CDP數據反演雨滴譜計算的Dm(a)、LWC(b)、Nt(c)、R(d) Nt以對數尺度顯示。

圖9 2016年7月2日20:26的0.5 °仰角NUIST-CDP ZH PPI圖(a)、常規反演雨滴譜計算的ZH(b)、最優化雨滴譜計算的ZH(c)、SA雷達ZH PPI圖(d)

圖10 2016年7月2日20:26的0.5 °仰角NUIST-CDP ZDR PPI圖(a)、最優化雨滴譜計算的ZDR(b)、NUIST-CDP KDP PPI圖(c)、最優化雨滴譜計算的KDP(d)

為了定量分析基于新型觀測算子的變分反演方法效果，先將2015 年8 月10 日11 時、2016 年7月2 日20 時的最優化雨滴譜計算結果和雨量計輸出結果對比。結果如圖11a、11b 所示，相關系數CC、均方根誤差RMSE、相對偏差RB 分別為0.870 7 和0.888 6、4.833 6 和2.094 4、29.523 2%和14.783 6%，可見最優化雨滴譜計算的R和實際雨強具有很高的相關性。圖11a的RMSE和RB比圖11b偏大較多的原因可能是：總體上對流云降水過程對應的ZH比層狀云降水對應的ZH更大，因此反演時更嚴重的衰減造成更大的誤差。

圖11 最優化雨滴譜計算結果和雨量計輸出結果散點圖 a. 2015年8月10日11時；b. 2016年7月2日20時。

再將第一個個例即2015 年8 月10 日03—14時12 小時的反演結果繪制成時序圖，并且將位于江寧站點的滴譜儀數據計算結果作為真值檢驗反演結果。同時，將基于C-G 模型的常規反演結果和基于模擬觀測算子的變分反演結果一同加入對比，繪制了圖12a～12d所示的四張時序圖，分別為R、LWC、Nt和Dm的反演結果。從圖12a、圖12b 可看出，總體上三種反演方法的結果隨時間變化趨勢都和真值一致，但在反演強降水時低估降水且存在較大誤差，這是因為降水越強時ZH受衰減效應影響越大，反演結果越不準確。基于新型觀測算子的變分反演結果在反演強降水時比另外兩種方法更加接近真值，并且在弱降水區幾乎貼合真值，說明新算法不僅反演更準確且衰減訂正更有優勢。將三種反演方法結果與真值的相關系數CC、均方根誤差RMSE、相對偏差RB對比，基于新型觀測算子的變分反演結果與真值有最高的相關性，CC 分別為0.973 9 和0.962 1，且RMSE 和RB最小，分別為21.553 6 和0.635 6，41.761 8%和25.189 6%。由圖12c 可得，常規方法和基于模擬觀測算子的反演結果出現了更多的波動，基于新型觀測算子的Nt反演結果更加穩定，并且有最大的CC(0.831 3) 和 最 小 的RMSE(0.324 0)、RB(4.6279%)。但Nt反演的準確度并沒有R和LWC那么高，這可能是因為使用高階矩量反演低階矩量時存在誤差。由圖12d可見，基于新型觀測算子Dm的反演結果依然最接近真值(CC=0.793 5，RMSE=0.472 4，RB=10.625 6%)，但準確度同樣不如R和LWC，原因可能是在階矩法計算Dm時，階矩組合為4 階矩和3 階矩之比，兩個階矩量分別計算時本就存在誤差，相比時造成了更大的誤差。

圖12 2015年8月10日03—14時反演雨滴譜計算結果時序圖(時間分辨率為7～8 min)

分析反演結果及結果優化原因，與常規反演方法相比，變分反演方法考慮了觀測數據的誤差和雨滴譜的不確定性，能夠最大化地降低誤差。與基于模擬觀測算子的變分反演方法相比，基于新型觀測算子的變分反演方法在推導觀測算子時，不需要建立雨滴譜模型，使用滴譜儀實測數據計算更能反映當地降水特點的觀測算子。前人在變分反演中所用的觀測算子多為建立雨滴譜模型后模擬得到，而某一種雨滴譜模型并不能代表所有的雨滴譜分布。滴譜儀數據提供了更能代表所在地區降水特性的雨滴譜信息，用其計算的觀測算子參與最優化反演，求解出的雨滴譜也更接近真實值。

5 結論與討論

本研究提出了一種基于新型觀測算子的雙偏振多普勒天氣雷達變分反演雨滴譜方法，并且利用理想模擬試驗和NUIST-CDP 實測個例驗證了該方法的可行性，主要結論如下。

(1) 由于變分反演在反演過程中考慮了觀測數據的誤差和雨滴譜的不確定性，且新型觀測算子能夠反映本地區降水的實際滴譜分布，解決了單一雨滴譜模型不能代表本地區雨滴譜分布特征的問題，基于新型觀測算子的變分反演方法可實現雨滴譜的最優化反演。基于新型觀測算子的變分反演方法得到的LWC 和Dm與滴譜儀數據計算結果的相關系數達到了0.96 和0.80，利用最優反演的LWC和Dm得到雨滴譜，從而計算微物理參數和降雨量，可實現降雨參數的最優估計。

(2) 本研究在改進觀測算子后提升了反演的精度，可看出變分反演雨滴譜的效果與觀測算子關系很大。由于不同水凝物的滴譜模型及滴譜參數變化較大，在未來的工作中將考慮針對不同相態多種類型水凝物粒子的成分、尺度、形狀分布等特點建立觀測算子，并應用到變分反演方法中以得到不同水凝物的最優滴譜反演結果，從而更好地了解不同水凝物粒子的微物理特性。