基于S波段雙極化雷達(dá)的變分法的定量降水估計(jì)算法

劉陳帥，張阿思，陳生

(1.中山大學(xué)大氣科學(xué)學(xué)院，廣東 珠海519082；2.廣東省氣候變化與自然災(zāi)害研究重點(diǎn)實(shí)驗(yàn)室，廣東 珠海519082；3.熱帶大氣海洋系統(tǒng)科學(xué)教育部重點(diǎn)實(shí)驗(yàn)室，廣東 珠海519082；4.南方海洋實(shí)驗(yàn)室（珠海），廣東 珠海519082；5.廣東省氣象臺(tái)，廣東 廣州510641；6.中國(guó)科學(xué)院西北生態(tài)環(huán)境資源研究院黑河遙感站和甘肅省遙感重點(diǎn)實(shí)驗(yàn)室，甘肅 蘭州730000)

1 引 言

上世紀(jì)七十年代，美國(guó)科學(xué)家提出了雙極化雷達(dá)理論。經(jīng)過(guò)數(shù)十年的發(fā)展，雙極化雷達(dá)已經(jīng)步入業(yè)務(wù)實(shí)用階段[1-2]。定量降水估計(jì)（QPE）是雷達(dá)氣象學(xué)的主要任務(wù)之一，與傳統(tǒng)的多普勒雷達(dá)相比，雙極化天氣雷達(dá)在QPE上有著優(yōu)異的表現(xiàn)。雷達(dá)測(cè)量誤差是雷達(dá)QPE的主要誤差之一，在進(jìn)行QPE之前，應(yīng)當(dāng)對(duì)測(cè)量的極化變量進(jìn)行重構(gòu)以消除測(cè)量中的隨機(jī)誤差。

傳統(tǒng)的QPE一般是通過(guò)水平反射率因子與降雨率（R）之間的冪律關(guān)系來(lái)反演降雨率。而雙極化雷達(dá)相較于傳統(tǒng)雷達(dá)，可以提供更多的信息如差分反射率因子（ZDR），比差分傳播相移（KDP）和差分相位（ΦDP），極化變量受到雨滴譜數(shù)據(jù)的影響更小，因此用極化變量估計(jì)降雨率有更好的表現(xiàn)[3]。

差分相位（ΦDP）為水平和垂直極化信號(hào)傳播相移的差分，比差分傳播相移（KDP）是ΦDP關(guān)于距離的導(dǎo)數(shù)。與基于水平反射率因子ZH降水估計(jì)算法R(ZH)相比，基于KDP的降水估計(jì)算法R(KDP)受雨滴譜影響更小。KDP是雷達(dá)的間接產(chǎn)品，在進(jìn)行定量降水估計(jì)之前需要對(duì)KDP進(jìn)行估計(jì)，如下面的公式（1）：

在實(shí)際應(yīng)用中，傳播中的隨機(jī)誤差和后向散射相位（backscattering phase）可能導(dǎo)致出現(xiàn)負(fù)值的KDP。Maesaka提出一種基于非負(fù)KDP的變分方法來(lái)估計(jì)ΦDP，通過(guò)計(jì)算邊界條件，以測(cè)量值和理論值之間的誤差為代價(jià)函數(shù)，通過(guò)逐次迭代來(lái)求解最優(yōu)的KDP[4]。當(dāng)徑向數(shù)據(jù)中存在異常值或者大量連續(xù)缺失值時(shí)，Maesaka的方法不能準(zhǔn)確求解出邊界條件，進(jìn)而導(dǎo)致錯(cuò)估KDP，而且計(jì)算邊界條件的方法會(huì)損失較多的信息，本文設(shè)計(jì)了一種比較穩(wěn)健的邊界條件求解方法，在保留更多的信息情況下，得到準(zhǔn)確的邊界條件。

變分方法的優(yōu)點(diǎn)就是在確定遠(yuǎn)-近邊界條件后，可以有效地減小隨機(jī)誤差的影響，擁有較好的擬合效果，在加入低通濾波器的情況下，擬合曲線也十分光滑。不足之處就是不能展現(xiàn)出小尺度的變化。解決方法就是以ZDR、ZH和KDP三個(gè)極化變量構(gòu)建的物理約束，通過(guò)ZDR和ZH兩個(gè)極化變量反演出具有較好精細(xì)尺度的KDP。

雖然基于KDP的降水估計(jì)算法R(KDP)相較于R(ZH)更少地受雨滴譜影響，但是對(duì)于較小的降雨率，當(dāng)ZH在20～30 dBZ時(shí)，基于KDP計(jì)算的降雨率估計(jì)并不理想[4-5]。Chen等[6-8]針對(duì)不同的水凝物采用不同的關(guān)系式。Zhang等[9]針對(duì)不同的KDP的值選取不同的降水估計(jì)方法進(jìn)行定量降水估計(jì)，獲得比較好的結(jié)果，本文中采用該方法進(jìn)行降水估計(jì)。

2 基于物理約束的變分方法

2.1 質(zhì)量控制

質(zhì)量控制算法可以有效地剔除非氣象回波的數(shù)據(jù)和地物雜波。其中由于地物雜波和異常傳播造成高ZH值可以通過(guò)異常值檢測(cè)的方法消除，但是由飛禽類等生物目標(biāo)造成的低ZH值很難被消除。首先剔除信噪較小（SNR<20 dB）的觀測(cè)值[10]。Valliappa Lakshmanan(2014)[11]提出了一種天氣雷達(dá)的極化變量的質(zhì)量控制方法，其主要是對(duì)水平反射率因子(ZH)、極化相關(guān)系數(shù)(ρHV)和差分反射率因子(ZDR)三個(gè)變量進(jìn)行閾值控制，具體方法為：

對(duì)ZH在徑向上使用1 km的滑動(dòng)窗口進(jìn)行滑動(dòng)平均，對(duì)ZDR和ρHV在徑向上使用2 km的滑動(dòng)窗口。對(duì)于ΦDP的質(zhì)量將在邊界條件部分詳細(xì)描述。由于本文應(yīng)用的變分法對(duì)純降水估測(cè)有效，因此選取融化層以下的數(shù)據(jù)。以0.5°仰角為例，選取150 km以內(nèi)的數(shù)據(jù)，最大高度約為1.3 km，根據(jù)一些華南降水的研究[12-13]，可以確保所選擇的數(shù)據(jù)在融化層高度以下。

2.2 物理約束

Scarchilli等[14]研究了水平反射率因子，差分反射率因子和比差分傳播相移三者之間的關(guān)系，確定了三者之間的物理約束，即可以通過(guò)任意ZDR和ZH兩個(gè)極化變量通過(guò)物理約束計(jì)算出KDP。三者之間的關(guān)系式如（5）式所示。

其中，ZH的單位是mm6/m3，ZDR的單位是dB，由ZH和ZDR兩個(gè)變量估計(jì)得到,φi為根據(jù)物理約束得到的反演得到的差分相位，φ0為徑向上的初始差分相位，本文中將近邊界條件作為初始差分相位，具體計(jì)算過(guò)程在2.3節(jié)。Scarchilli等[14]使用大量S波段的觀測(cè)數(shù)據(jù)進(jìn)行擬合，得到了公式（5）的經(jīng)驗(yàn)參數(shù)，其 中C=1.05×10-4、α=0.96、β=0.26。本文基于Scarchilli給出的經(jīng)驗(yàn)參數(shù)，結(jié)合本地觀測(cè)數(shù)據(jù)給出了本地化參數(shù)，其中C=1.37×10-4、α=0.88、β=0.20。通過(guò)公式（5）和（6）可以通過(guò)估計(jì)的KDP,進(jìn)而對(duì)ΦDP進(jìn)行重構(gòu)。兩組經(jīng)驗(yàn)參數(shù)組合實(shí)際應(yīng)用情況如圖1所示，本地化參數(shù)可以很好反映出觀測(cè)數(shù)據(jù)的趨勢(shì)。

圖1中的黑色星點(diǎn)為一條ΦDP徑向觀測(cè)數(shù)據(jù),紅色實(shí)線是基于Scarchilli給出的經(jīng)驗(yàn)參數(shù)進(jìn)行的ΦDP的重構(gòu)，藍(lán)色實(shí)線是基于本地化參數(shù)進(jìn)行的ΦDP的重構(gòu)?；诮?jīng)驗(yàn)參數(shù)重構(gòu)的ΦDP整體上高于觀測(cè)數(shù)據(jù)，而基于本地化參數(shù)重構(gòu)的ΦDP位于觀測(cè)數(shù)據(jù)之間，較好地?cái)M合出觀測(cè)數(shù)據(jù)的大致走勢(shì)。

圖1 基于物理約束的ΦDP重構(gòu)

在經(jīng)過(guò)異常值檢測(cè)并剔除后，ΦDP徑向觀測(cè)數(shù)據(jù)可能會(huì)變得不連續(xù)，因此需基于物理約束重構(gòu)的ΦDP對(duì)徑向觀測(cè)數(shù)據(jù)的缺失值進(jìn)行填補(bǔ)。

2.3 邊界條件

變分方法需要計(jì)算近端和遠(yuǎn)端的邊界條件Φnear和Φfar，以保證重構(gòu)的ΦDP在兩個(gè)邊界條件之間。由于Maesaka(2012)計(jì)算邊界條件的步驟較為簡(jiǎn)單，對(duì)異常數(shù)據(jù)的影響較敏感，且計(jì)算邊界條件需選取連續(xù)20個(gè)有效樣本進(jìn)行邊界條件的計(jì)算，信息損失較大。基于上述缺陷，我們進(jìn)行了一定的修正，在盡可能保留足夠多信息的情況下，使邊界條件的計(jì)算更加穩(wěn)健。首先對(duì)ΦDP觀測(cè)數(shù)據(jù)進(jìn)行質(zhì)量控制，以一條徑向數(shù)據(jù)為例，在剔除ρHV<0.9的數(shù)據(jù)之后，通過(guò)遍歷所有數(shù)據(jù)計(jì)算連續(xù)數(shù)組的切片索引。之后對(duì)索引進(jìn)行如下處理。(1)若兩索引之間間隔小于5，且前一個(gè)索引最后的數(shù)據(jù)和后一個(gè)索引的首個(gè)數(shù)據(jù)相差小于30°，將兩個(gè)索引合并為一個(gè)索引。(2)若索引內(nèi)的數(shù)據(jù)個(gè)數(shù)小于3個(gè)，將被視為孤立數(shù)據(jù)被刪除。(3)遍歷所有數(shù)據(jù)如果某個(gè)距離庫(kù)的ΦDP的值與相鄰的數(shù)據(jù)差值大于35°，則將該數(shù)據(jù)視為被雜波污染的數(shù)據(jù)并使用兩個(gè)相鄰數(shù)據(jù)的平均值進(jìn)行替換。進(jìn)行質(zhì)量控制示意圖如圖2所示。橙色的距離庫(kù)表示缺失值，紅色的距離庫(kù)表示異常值徑向數(shù)據(jù)（a）有兩個(gè)連續(xù)數(shù)組的切片索引，分別為（1，6）和（10，16），這兩個(gè)索引對(duì)應(yīng)的距離庫(kù)間隔為3，小于5，則將兩個(gè)索引合并為一個(gè)索引，即（1，16）。徑向數(shù)據(jù)（b）有兩個(gè)連續(xù)數(shù)組的切片索引，分別為（4，5）和（11，11），這兩個(gè)索引對(duì)應(yīng)的距離庫(kù)長(zhǎng)度均小于3，因此視為孤立點(diǎn)并刪除。徑向數(shù)據(jù)（c）有一個(gè)異常值，即 |φDP7-φDP8|>35°, |φDP6-φDP7|>35°，因此視為被雜波污染的數(shù)據(jù)并使用相鄰的距離庫(kù)的平均值進(jìn)行替換。

圖2 ΦDP觀測(cè)數(shù)據(jù)質(zhì)量控制示意圖

然后計(jì)算近邊界條件Φnear，步驟如下。

(1)從頭到尾遍歷徑向數(shù)據(jù)的索引，若索引所代表的數(shù)組元素大于20個(gè)，則該數(shù)組前20個(gè)數(shù)據(jù)作為近邊界條件計(jì)算的樣本。

(2)計(jì)算選取的樣本(X1,X2,……,X20)與對(duì)應(yīng)距離(r1,r2,……,r20)的線性回歸的斜率K。

(3)如果K>0，則近端邊界條件Φnear為最近距離r1的預(yù)測(cè)值，否則，近端邊界條件為所取樣本的中位數(shù)。

近邊界條件可以等效視為系統(tǒng)初始差分相位φ0，在計(jì)算完近邊界條件之后，可以使用物理約束重構(gòu)ΦDP，然后對(duì)徑向數(shù)據(jù)中的缺失值進(jìn)行填補(bǔ)。遠(yuǎn)端邊界條件Φfar的計(jì)算方法與近端邊界條件相似。

(1)從尾到頭遍歷徑向數(shù)據(jù)的索引，若索引所代表的數(shù)組元素大于20個(gè)，則該數(shù)組后20個(gè)數(shù)據(jù)作為遠(yuǎn)邊界條件計(jì)算的樣本。

(2)計(jì)算選取的樣本(Xend-19,Xend-18,……,Xend)與對(duì)應(yīng)距離(rend-19,rend-18,……,rend)的線性回歸的斜率K。

(3)如果K>0，則遠(yuǎn)端邊界條件Φfar為最遠(yuǎn)距離rend的預(yù)測(cè)值，否則，遠(yuǎn)端邊界條件為所取樣本的中位數(shù)。

圖3為遠(yuǎn)-近邊界條件計(jì)算方法的示意圖，近邊界條件為選取最近的20個(gè)有效樣本點(diǎn)進(jìn)行線性擬合，遠(yuǎn)邊界條件為選取最遠(yuǎn)的20個(gè)有效樣本點(diǎn)進(jìn)行線性擬合，遠(yuǎn)近邊界條件受到噪聲的影響較小，在本文中，將近邊界條件視為該徑向方向上的初始相位?0，初始相位可以在填補(bǔ)數(shù)據(jù)起到作用。

圖3 邊界條件計(jì)算

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

在構(gòu)造代價(jià)函數(shù)之前，需要構(gòu)造幾個(gè)中間變量，對(duì)于一條有N個(gè)距離庫(kù)的雷達(dá)徑向觀測(cè)數(shù)據(jù)，本文定義差分相位觀測(cè)值(ΦDP)i，理論的差分 相 位為(?DP)i，比差分傳播相移(KDP)i，其中i=1,2,……,N。

然后給出的φ的定義：

根據(jù)（1）式中所給出的K DP的定義，這里引入一個(gè)前向算子H1，φ可以表示為：

其中?r是雷達(dá)距離分辨率。因?yàn)棣礑P是遞增的，KDP是非負(fù)的，為了保證KDP非負(fù)，引入了k2，其表示如下：

然后將式（10）帶入式（9），φ的表達(dá)式可以寫為：

通過(guò)引入后項(xiàng)算子H2,φ'可以寫成與式（11）相同的形式：

觀測(cè)的差分相位與邊界條件的差值為：

Jobs項(xiàng)是觀測(cè)項(xiàng)和重構(gòu)的理論項(xiàng)的均方誤差，是代價(jià)函數(shù)中的主要部分，使重構(gòu)的差分相位更好地?cái)M合觀測(cè)的差分相位。Jlpf是參數(shù)k的拉普拉斯算子，相當(dāng)于一個(gè)低通濾波器，Clpf是濾波器的參數(shù)，一般取值較大。將當(dāng)代價(jià)函數(shù)取最小值時(shí)的k作為最終的解，然后通過(guò)k計(jì)算KDP，本文選取的迭代方法為擬牛頓法，該方法需要代價(jià)函數(shù)的偏導(dǎo)數(shù)作為迭代方向。代價(jià)函數(shù)關(guān)于k的偏導(dǎo)數(shù)如下所示：

通過(guò)代價(jià)函數(shù)和關(guān)于k的偏導(dǎo)數(shù)，使用牛頓迭代法進(jìn)行求解，最終重構(gòu)ΦDP[15]。

圖4為變分方法擬合和物理約束的效果圖，其中黑實(shí)線為徑向數(shù)據(jù)觀測(cè)值ΦDP，藍(lán)實(shí)線為變分?jǐn)M合的?DP，紅實(shí)線為變分?jǐn)M合的KDP，紫色實(shí)線為物理約束KDP，綠色實(shí)線為物理約束基于物理約束KDP得到的?D（P簡(jiǎn)稱物理約束?DP）。從圖中可以看出變分?jǐn)M合的?DP可以很好地反映出觀測(cè)值ΦDP的趨勢(shì)，并且消除了大部分的隨機(jī)誤差，物理約束KDP在30 km處有一個(gè)較大的異常值，相比于物理約束KDP，變分?jǐn)M合KDP表現(xiàn)得更加平滑，且符合實(shí)際情況。

圖4 變分?jǐn)M合結(jié)果

2016年5月7日06點(diǎn)54分0.5°仰角?DP、KDP的掃描數(shù)據(jù)圖如圖2a~2b所示，經(jīng)過(guò)質(zhì)控和填補(bǔ)后的結(jié)果如圖5c~5d所示。在?DP、KDP的原始數(shù)據(jù)里存在局部范圍的缺失值，如雷達(dá)的東北方向和西北方向，在經(jīng)過(guò)質(zhì)量控制和物理約束的填補(bǔ)后數(shù)據(jù)在空間變得連續(xù)，過(guò)濾掉大部分的噪聲，并且保留了主要的回波部分。圖5e~5f是經(jīng)過(guò)2.4節(jié)所描述的變分方法擬合后的數(shù)據(jù)圖像，可以看出經(jīng)過(guò)變分?jǐn)M合后，?DP和KDP的數(shù)據(jù)在空間上變得光滑，可以更好地反映出雷達(dá)KDP的數(shù)據(jù)。

圖5 2016年5月7日06點(diǎn)54分0.5°仰角?DP、K DP的數(shù)據(jù)圖

3 實(shí)證分析

3.1 數(shù)據(jù)與降水估計(jì)方法

本文選取廣東省2 400多個(gè)雨量站作為地面實(shí)際觀測(cè)數(shù)據(jù)，主要研究2017年5月7日00—24時(shí)的降水事件，當(dāng)天最大累計(jì)降水超過(guò)300 mm。雷達(dá)數(shù)據(jù)選取廣東省廣州市的一個(gè)S波段的極化雷達(dá)。

KDP對(duì)小雨滴的降水估計(jì)效果較差[8]。Zhang等[9]運(yùn)用了一種組合降水估計(jì)方法（RQPE），該方法使用ZH、ZDR和KDP三個(gè)極化變量進(jìn)行估測(cè)，并給出了基于2014年4月—2015年5月的位于廣東省西南部的2DVD觀測(cè)到的雨滴譜數(shù)據(jù)擬合的經(jīng)驗(yàn)公式，經(jīng)驗(yàn)公式如下所示：

其中，ZH(mm6/m3)是水平反射率因子，ZDRl=10ZDR/10是線性尺度上的差分反射率因子。本文使用經(jīng)過(guò)變分重構(gòu)后的KDP數(shù)據(jù)，基于Zhang等[9]的RQPE方法和經(jīng)驗(yàn)公式進(jìn)行降水估計(jì)，整個(gè)算法（V-RQPE）的流程圖如圖6所示。其中，T1=38 dBZ，T2=0.5 dB，T3=0.1°/km。

圖6 V-RQPE算法流程圖

首先用克里金插值方法（KRIG）將雨量站插值到0.1°×0.1°的網(wǎng)格場(chǎng)上，獲得逐小時(shí)的降水產(chǎn)品，同時(shí)將基于V-RQPE估測(cè)的降水量也進(jìn)行重采樣轉(zhuǎn)換到0.1°×0.1°的網(wǎng)格場(chǎng)上，以方便比較降雨事件的空間分布。圖7為基于KRIG插值、VRQPE的24小時(shí)累計(jì)降水估測(cè)值圖像，V-RQPE在雷達(dá)中心位置可以很好估計(jì)出降水量，但是在雷達(dá)的北部估測(cè)效果與實(shí)際值相比偏低。

圖7 24小時(shí)累計(jì)降水觀測(cè)值與估測(cè)值

3.2 評(píng)估指標(biāo)

在該部分，本文選取了2017年5月7號(hào)廣州市00—24時(shí)的降水過(guò)程進(jìn)行試驗(yàn)，其中最大的降水量超過(guò)100 mm/h。在評(píng)估過(guò)程中，采用相關(guān)系數(shù)（CC）、均方誤差（RMSE）、相對(duì)偏置（RB）和分?jǐn)?shù)均方差（FRMSE）四個(gè)評(píng)估標(biāo)準(zhǔn)來(lái)評(píng)估該變分方法的性能：

公式(25)~(28)分別是相關(guān)系數(shù)（CC）、均方誤差（RMSE）和相對(duì)偏置（RB）的計(jì)算公式，其中Rguage表示雨量站的觀測(cè)值，guage為Rguage的對(duì)應(yīng)0.1°×0.1°的網(wǎng)格點(diǎn)的空間平均值，RQPE為估測(cè)的降雨率，QPE為RQPE的平均值。相關(guān)系數(shù)（CC）表示QPE估測(cè)值和實(shí)際觀測(cè)值的相關(guān)關(guān)系，相關(guān)關(guān)系越大，估測(cè)效果越好。均方誤差（RMSE）表示觀測(cè)值與QPE估測(cè)值的差距，RMSE越大，QPE的估測(cè)效果越差。相對(duì)偏置（RB）表示QPE估測(cè)值與觀測(cè)值的相對(duì)誤差的程度，如果大于0則表示降雨率被高估，小于0則表示降雨率被低估。分?jǐn)?shù)均方差（FRSME）為均方誤差（RMSE）與觀測(cè)值均值之比，是一個(gè)無(wú)量綱的統(tǒng)計(jì)量，在針對(duì)不同量級(jí)的降水事件時(shí)，可以更好評(píng)估估測(cè)誤差。

3.3 累計(jì)降水估測(cè)效果

本文使用2017年5月7日00—24時(shí)的廣州雷達(dá)第一仰角（0.5°）觀測(cè)數(shù)據(jù)，圖5a~5b為06點(diǎn)54分0.5°仰角?DP、KDP的PPI掃描圖像。經(jīng)過(guò)質(zhì)量控制之后，使用基于物理約束的變分方法，對(duì)ΦDP進(jìn)行變分?jǐn)M合，然后利用3.1節(jié)所描述的RQPE算法進(jìn)行降水估計(jì)，24小時(shí)累計(jì)降水估計(jì)如圖7b所示。之后將雨量站進(jìn)行重采樣轉(zhuǎn)換到0.1°×0.1°的網(wǎng)格場(chǎng)上的數(shù)據(jù)與RQPE算法估計(jì)的降水量數(shù)據(jù)繪制成散點(diǎn)圖，并計(jì)算相應(yīng)的評(píng)估指標(biāo)，圖8分別為1、3、6、24小時(shí)累計(jì)降水的散點(diǎn)圖，其色標(biāo)為散點(diǎn)相對(duì)密度，紅色相對(duì)密度最大，藍(lán)色最小。

圖8 累計(jì)降水散點(diǎn)圖

從圖8中可以看出，對(duì)于不同時(shí)長(zhǎng)的累計(jì)降水估計(jì)，RQPE的性能表現(xiàn)也有所差異。一般來(lái)說(shuō)，隨著時(shí)長(zhǎng)的增加，相關(guān)系數(shù)會(huì)略微增大（CC從0.778到0.838），均方誤差也會(huì)增大（RMSE從4.06到23.73），但是分?jǐn)?shù)均方差（FRMSE）卻有一個(gè)起伏的過(guò)程，從不同時(shí)間的累積降水對(duì)比圖可以看出，小降雨量存在高估的現(xiàn)象，大降雨量存在低估的現(xiàn)象，產(chǎn)生這種誤差的可能是雷達(dá)數(shù)據(jù)和雨量站數(shù)據(jù)的質(zhì)量，RQPE算法中的不確定性造成的。

3.4 多種QPE方法比較

為了比較基于變分?jǐn)M合的降水估計(jì)方法的性能，本文采用了六種不同的降水估計(jì)方案進(jìn)行對(duì)比，具體方案如表1所示。

表1 六種QPE方案配置表

表1顯示了六種不同的降水方案的配置，分別V-RQPE采用本文的變分?jǐn)M合方法，使用RQPE的方法進(jìn)行降水估計(jì)；RQPE采用Zhang等[9]的質(zhì)控方法進(jìn)行降水估計(jì)；R-Z采用公式（21）進(jìn)行降水估計(jì)；R-Z-ZDR采用公式（23）進(jìn)行降水估計(jì)；RVKDP采用公式（22），基于變分?jǐn)M合的KDP數(shù)據(jù)進(jìn)行降水估計(jì)；R-OKDP采用公式（22），基于滑動(dòng)平均和物理約束填補(bǔ)的KDP數(shù)據(jù)進(jìn)行降水估計(jì)。圖9為上述六種QPE算法一小時(shí)累計(jì)降水估計(jì)的散點(diǎn)圖。

從圖9可以看出各種QPE方法的差異比較明顯，除了R-OKDP方法外，其他五個(gè)QPE方法的相關(guān)系數(shù)均較高（CC>0.7）。兩個(gè)組合算法V-RQPE和RQPE中，V-RQPE的CC（0.77）、RMSE(4.06)和FRMSE（2.769）相對(duì)較小，但是存在低估的現(xiàn)象。從圖9b中看出，雖然RQPE的RB絕對(duì)值較小，但是對(duì)小降雨量存在高估現(xiàn)象，大降雨量存在低估的現(xiàn)象。其他四個(gè)單一算法里R-Z的相關(guān)系數(shù)最高(0.841)，但其低估嚴(yán)重（-72.77%），估測(cè)效果較差。R-Z-ZDR算法也擁有較高的相關(guān)系數(shù)（0.801 1），但是其在降雨率較大時(shí)估計(jì)誤差明顯增大。R-VKDP和R-OKDP兩個(gè)算法相比，但是R-VKDP擁有更好的統(tǒng)計(jì)性能。就單一算法而言，R-VKDP估測(cè)效果較好。

圖9 六種QPE方法的1小時(shí)累計(jì)降水估計(jì)對(duì)比

表2為6種QPE方法計(jì)算的24小時(shí)累計(jì)降水的評(píng)估結(jié)果。與1小時(shí)累計(jì)降水估計(jì)相比，24小時(shí)累計(jì)降水的RMSE均有不同程度的增大，24小時(shí)累計(jì)降水由于誤差累積，可能出現(xiàn)誤差抵消，所以六種方法的相關(guān)系數(shù)都有一定升高。就四個(gè)評(píng)估指標(biāo)而言，V-RQPE和R-VKDP的估測(cè)效果是最好的，RQPE次之，R-OKDP估測(cè)效果最差。

表2 六種QPE方法24小時(shí)累計(jì)降水估計(jì)評(píng)估表

4 結(jié) 論

本文將物理約束應(yīng)用到ΦDP觀測(cè)值的缺失填補(bǔ)中，使用了更加穩(wěn)健的方法求解邊界條件，這種方法可以在穩(wěn)健地求解遠(yuǎn)近邊界條件的同時(shí)，保留更多的原始信息。之后使用變分方法重構(gòu)ΦDP，進(jìn)而反演出非負(fù)的KDP，并將其應(yīng)用到定量降水估計(jì)中。以2017年5月7日的廣州S波段雷達(dá)的回波數(shù)據(jù)為例，使用了Zhang等[9]的RQPE方法進(jìn)行降水估計(jì)，并使用了其他五種不同的方法進(jìn)行對(duì)比，采用相關(guān)系數(shù)、均方誤差、均方誤差、分?jǐn)?shù)均方差等指標(biāo)進(jìn)行評(píng)估，驗(yàn)證了基于變分方法的降水估計(jì)的性能。

（1）變分?jǐn)M合的?DP保留了觀測(cè)值ΦDP的大多數(shù)的信息，并消除了隨機(jī)誤差。對(duì)于24小時(shí)累計(jì)降水，RQPE算法可以較準(zhǔn)確地估計(jì)降水中心的降水。不同時(shí)長(zhǎng)的累計(jì)降水估計(jì)性能表現(xiàn)不同，但存在一定低估的現(xiàn)象。

（2）對(duì)于六種不同的降水估計(jì)算法，在1小時(shí)累計(jì)降水估計(jì)的評(píng)估中，V-RQPE擁有較高的相關(guān)系數(shù)（CC=0.77），較低的估計(jì)誤差（RMSE=4.06，F(xiàn)RMSE=2.76）；R-Z相關(guān)系數(shù)最高（CC=0.84），但其低估嚴(yán)重（RB=-72%），估計(jì)誤差也較大；R-ZZDR在降雨率大估計(jì)誤差會(huì)變大；R-VKDP和RSKDP估計(jì)性能相似，但是R-SKDP會(huì)出現(xiàn)負(fù)的降雨率；R-OKDP表現(xiàn)最差。在24小時(shí)累計(jì)降水估計(jì)的評(píng)估中，隨著時(shí)間的增加，降水量也會(huì)增加，RMSE在增加，但是無(wú)量綱的FRMSE在減小。六種QPE方法中，V-RQPE和R-VKDP的評(píng)估指標(biāo)是最好的，擁有較高的相關(guān)系數(shù)、較低的估測(cè)誤差。

本文所設(shè)計(jì)的變分方法可以很好消除KDP的誤差，不論實(shí)在組合算法中還是單一算法中，都可以有效地降低定量降水估計(jì)中的不確定性，但是依然存在一定的低估現(xiàn)象，可能是由于雷達(dá)數(shù)據(jù)和雨量站數(shù)據(jù)的質(zhì)量，RQPE算法中的不確定性造成的，也可能是由于平滑濾波器的選擇問(wèn)題，雖然最大限度地減少了隨機(jī)誤差，但可能會(huì)使KDP變小導(dǎo)致低估。消除誤差會(huì)伴隨著損失信息，兩者如何才能達(dá)到均衡，是本文以后的主要研究方向。