基于新型觀測(cè)算子的雙偏振雷達(dá)雨滴譜變分反演

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

(1. 南京信息工程大學(xué)大氣物理學(xué)院，江蘇 南京 210044；2. 福建省平潭綜合實(shí)驗(yàn)區(qū)氣象局，福建 平潭 350400；3. 南京信息工程大學(xué)氣象災(zāi)害預(yù)警與評(píng)估協(xié)同創(chuàng)新中心，江蘇 南京 210044；4. 平潭海洋氣象野外科學(xué)觀測(cè)研究站，福建 平潭 350400；5. 福建省災(zāi)害天氣重點(diǎn)實(shí)驗(yàn)室，福建 福州 350001）

1 引 言

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

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

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

無(wú)論是常規(guī)反演還是變分反演，都需要在反演前建立雨滴譜模型，利用模型參數(shù)表示雨滴譜。但在實(shí)際情況中，雨滴譜是復(fù)雜且多變的，某一種雨滴譜模型并不能代表所有的雨滴譜分布。前人在變分反演中所用的觀測(cè)算子多為建立雨滴譜模型后模擬得到，本研究不建立雨滴譜模型，使用滴譜儀實(shí)測(cè)數(shù)據(jù)計(jì)算更能反映當(dāng)?shù)亟邓攸c(diǎn)的觀測(cè)算子，將其應(yīng)用于雙偏振雷達(dá)雨滴譜最優(yōu)化反演中，進(jìn)一步提高了雨滴譜反演的準(zhǔn)確性。

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

2 數(shù)據(jù)介紹

2.1 C波段雙偏振多普勒天氣雷達(dá)

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

2.2 S波段多普勒天氣雷達(dá)

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

2.3 OTT Parsivel滴譜儀

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

2.4 雨量計(jì)

本文所用的雨量計(jì)數(shù)據(jù)來(lái)源于南京市地面雨量自動(dòng)站，共計(jì)85 個(gè)雨量站點(diǎn)。雨量計(jì)數(shù)據(jù)時(shí)間分辨率為1 小時(shí)，雨量分辨率為0.1 mm。本文以雨量計(jì)數(shù)據(jù)為實(shí)際雨強(qiáng)來(lái)驗(yàn)證最優(yōu)雨滴譜計(jì)算降雨率的準(zhǔn)確性。NUIST-CDP、SA 雷達(dá)、三部滴譜儀以及雨量計(jì)站點(diǎn)的空間分布如圖1所示。

圖1 NUIST-CDP、SA雷達(dá)、三部滴譜儀和雨量計(jì)站點(diǎn)分布圖

3 雨滴譜反演方法

3.1 代價(jià)函數(shù)

基于變分理論的雙偏振雷達(dá)雨滴譜反演本質(zhì)上是求狀態(tài)變量的先驗(yàn)估計(jì)、背景誤差、觀測(cè)誤差、雙偏振觀測(cè)量和觀測(cè)算子預(yù)測(cè)結(jié)果組成的非線性方程組最優(yōu)解的過(guò)程。因?yàn)橛甑巫V參數(shù)和Dm、LWC之間可相互轉(zhuǎn)化，且Dm和LWC符合高斯分布[21]，本研究使用Dm和LWC 作為狀態(tài)變量，建立雙偏振觀測(cè)量和雨滴譜之間的關(guān)系。將變分反演雨滴譜的代價(jià)函數(shù)定義為：

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

式(2)～(4)中，Dm(n)和LWC(n)代表第n個(gè)距離庫(kù)上的狀態(tài)變量，當(dāng)狀態(tài)變量確定時(shí)，可用觀測(cè)算子預(yù)測(cè)雙偏振觀測(cè)量。xb為狀態(tài)變量的先驗(yàn)估計(jì)(背景場(chǎng))，由狀態(tài)變量和雙偏振參量之間的經(jīng)驗(yàn)關(guān)系計(jì)算得到。B 為背景場(chǎng)的誤差協(xié)方差矩陣，通過(guò)對(duì)狀態(tài)變量進(jìn)行南京地區(qū)氣候統(tǒng)計(jì)分析后，將背景誤 差 標(biāo)準(zhǔn)差σDm和σLWC分 別 設(shè)置為1 mm 和0.707 g/cm3(方差分別為1 mm2和0.5 g2/cm6)，本研究假設(shè)了一個(gè)高斯相關(guān)模型來(lái)計(jì)算協(xié)方差矩陣[28]：

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

式(6)～(9)中，分號(hào)代表矩陣的縱向合并，n表示距離庫(kù)序號(hào)，Z'H和Z'DR表示衰減訂正后的ZH和ZDR?？紤]衰減效應(yīng)，當(dāng)由Dm和LWC得到雨滴譜并計(jì)算AH和ADP之后，各距離庫(kù)上ZH和ZDR可表示為：

式(10)、(11)中，δr代表雷達(dá)徑向上每個(gè)距離庫(kù)的庫(kù)長(zhǎng)度(即雷達(dá)徑向分辨率)。H(x)為觀測(cè)算子預(yù)測(cè)得到的雙偏振參量，觀測(cè)算子的建立在后文中會(huì)詳細(xì)說(shuō)明。R 為觀測(cè)誤差協(xié)方差矩陣，由三個(gè)雙偏振觀測(cè)量的觀測(cè)誤差標(biāo)準(zhǔn)差組成。對(duì)雙偏振雷達(dá)觀測(cè)量進(jìn)行統(tǒng)計(jì)分析后，將ZH、ZDR和KDP的觀測(cè)誤差標(biāo)準(zhǔn)差分別設(shè)置為1 dB、0.2 dB 和0.6 °/km，假定ZH、ZDR和KDP的誤差相互獨(dú)立時(shí)，

式(12)中，O 代表零矩陣，RZH、RZDR和RKDP分別對(duì)應(yīng)于三個(gè)雙偏振參量的觀測(cè)誤差協(xié)方差矩陣。由于雷達(dá)觀測(cè)量的隨機(jī)誤差在各距離庫(kù)上相互獨(dú)立，因此RZH、RZDR和RKDP副對(duì)角線上的元素為0。

因?yàn)橛^測(cè)算子是非線性的，所以代價(jià)函數(shù)最小值的求解是一個(gè)非線性最優(yōu)化問(wèn)題。該類問(wèn)題(例如變分資料同化問(wèn)題)通常用高斯-牛頓迭代方法求解，該方法的求解公式為[29]：

式(13)、(14)中，W 是最優(yōu)化權(quán)重矩陣，H 是觀測(cè)算子H的雅可比式，即觀測(cè)算子的線性近似，由每個(gè)觀測(cè)算子分別對(duì)兩個(gè)狀態(tài)變量求偏導(dǎo)數(shù)得到。

3.2 觀測(cè)算子的建立

降雨參數(shù)和雙偏振參量都由雨滴譜決定，因此降雨參數(shù)和雙偏振參量之間可定量關(guān)聯(lián)，這也是由雙偏振參量反演雨滴譜并計(jì)算降雨參數(shù)的基本原理。將降雨參數(shù)中的Dm、LWC(即狀態(tài)變量)和雙偏振參數(shù)之間的函數(shù)關(guān)系定義為觀測(cè)算子。

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

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

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

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

狀態(tài)變量用階矩法計(jì)算，可用雨滴譜的積分量表示，一般用Mn來(lái)表示雨滴譜的n階矩：

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

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

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

圖2 新型觀測(cè)算子、模擬觀測(cè)算子和滴譜儀數(shù)據(jù)計(jì)算值對(duì)比

3.3 反演步驟

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

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

首先，導(dǎo)入雙偏振雷達(dá)數(shù)據(jù)，對(duì)其進(jìn)行質(zhì)量控制。再利用狀態(tài)變量和雙偏振參量間的經(jīng)驗(yàn)關(guān)系計(jì)算狀態(tài)變量的先驗(yàn)估計(jì)。然后，如3.1 節(jié)所述，建立背景誤差協(xié)方差矩陣和觀測(cè)誤差協(xié)方差矩陣。隨后，導(dǎo)入預(yù)處理過(guò)后的滴譜儀數(shù)據(jù)，利用3.2 節(jié)中的公式計(jì)算雙偏振參量和狀態(tài)變量，推導(dǎo)新型觀測(cè)算子并計(jì)算觀測(cè)算子關(guān)于狀態(tài)變量的微分，進(jìn)而構(gòu)建代價(jià)函數(shù)，將觀測(cè)量、狀態(tài)變量的先驗(yàn)估計(jì)和誤差協(xié)方差矩陣代入其中，再由新型觀測(cè)算子將狀態(tài)變量轉(zhuǎn)化為觀測(cè)量。最后利用高斯-牛頓迭代法求解代價(jià)函數(shù)，當(dāng)結(jié)果收斂或達(dá)到迭代次數(shù)后得到最優(yōu)雨滴譜，結(jié)束迭代，否則更新?tīng)顟B(tài)變量，繼續(xù)迭代。得到最優(yōu)雨滴譜后由3.2 節(jié)中相關(guān)公式計(jì)算雙偏振參量和降雨參數(shù)。

4 實(shí)驗(yàn)驗(yàn)證

為驗(yàn)證上述基于新型觀測(cè)算子的變分反演方法的可行性，本節(jié)利用一個(gè)理想模擬試驗(yàn)和NUIST-CDP實(shí)測(cè)個(gè)例最優(yōu)化反演雨滴譜。

4.1 理想模擬試驗(yàn)

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

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

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

4.2 實(shí)測(cè)個(gè)例驗(yàn)證

除了理想模擬試驗(yàn)，本研究還利用NUISTCDP觀測(cè)到的一次對(duì)流云降水天氣個(gè)例和一次層狀云降水天氣個(gè)例驗(yàn)證該算法。

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

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

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

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

本研究提出的最優(yōu)雨滴譜反演方法能得到很好的滴譜反演效果，但從反演結(jié)果可看出，圖像并不完整，這是因?yàn)楦咚?牛頓迭代方法要求進(jìn)行迭代的向量必須是連續(xù)的，因此本研究在變分反演時(shí)對(duì)每條雷達(dá)徑向數(shù)據(jù)只選取第一段連續(xù)的部分。另外，在實(shí)際應(yīng)用當(dāng)中，如果有非雨滴的水凝物(例如冰雹)出現(xiàn)，可能會(huì)出現(xiàn)沒(méi)有解或異常解的情況，這是因?yàn)楸狙芯恐杏^測(cè)算子推導(dǎo)的前提是粒子為雨滴，這也是后續(xù)研究的主要內(nèi)容之一。

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

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

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

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

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

圖11 最優(yōu)化雨滴譜計(jì)算結(jié)果和雨量計(jì)輸出結(jié)果散點(diǎn)圖 a. 2015年8月10日11時(shí)；b. 2016年7月2日20時(shí)。

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

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

分析反演結(jié)果及結(jié)果優(yōu)化原因，與常規(guī)反演方法相比，變分反演方法考慮了觀測(cè)數(shù)據(jù)的誤差和雨滴譜的不確定性，能夠最大化地降低誤差。與基于模擬觀測(cè)算子的變分反演方法相比，基于新型觀測(cè)算子的變分反演方法在推導(dǎo)觀測(cè)算子時(shí)，不需要建立雨滴譜模型，使用滴譜儀實(shí)測(cè)數(shù)據(jù)計(jì)算更能反映當(dāng)?shù)亟邓攸c(diǎn)的觀測(cè)算子。前人在變分反演中所用的觀測(cè)算子多為建立雨滴譜模型后模擬得到，而某一種雨滴譜模型并不能代表所有的雨滴譜分布。滴譜儀數(shù)據(jù)提供了更能代表所在地區(qū)降水特性的雨滴譜信息，用其計(jì)算的觀測(cè)算子參與最優(yōu)化反演，求解出的雨滴譜也更接近真實(shí)值。

5 結(jié)論與討論

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

(1) 由于變分反演在反演過(guò)程中考慮了觀測(cè)數(shù)據(jù)的誤差和雨滴譜的不確定性，且新型觀測(cè)算子能夠反映本地區(qū)降水的實(shí)際滴譜分布，解決了單一雨滴譜模型不能代表本地區(qū)雨滴譜分布特征的問(wèn)題，基于新型觀測(cè)算子的變分反演方法可實(shí)現(xiàn)雨滴譜的最優(yōu)化反演?；谛滦陀^測(cè)算子的變分反演方法得到的LWC 和Dm與滴譜儀數(shù)據(jù)計(jì)算結(jié)果的相關(guān)系數(shù)達(dá)到了0.96 和0.80，利用最優(yōu)反演的LWC和Dm得到雨滴譜，從而計(jì)算微物理參數(shù)和降雨量，可實(shí)現(xiàn)降雨參數(shù)的最優(yōu)估計(jì)。

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