基于S波段雙極化雷達的變分法的定量降水估計算法

劉陳帥，張阿思，陳生

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

1 引 言

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

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

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

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

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

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

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

2.1 質(zhì)量控制

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

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

2.2 物理約束

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

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

圖1中的黑色星點為一條ΦDP徑向觀測數(shù)據(jù),紅色實線是基于Scarchilli給出的經(jīng)驗參數(shù)進行的ΦDP的重構(gòu)，藍色實線是基于本地化參數(shù)進行的ΦDP的重構(gòu)。基于經(jīng)驗參數(shù)重構(gòu)的ΦDP整體上高于觀測數(shù)據(jù)，而基于本地化參數(shù)重構(gòu)的ΦDP位于觀測數(shù)據(jù)之間，較好地擬合出觀測數(shù)據(jù)的大致走勢。

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

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

2.3 邊界條件

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

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

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

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

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

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

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

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

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

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

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

圖3 邊界條件計算

2.4 代價函數(shù)

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

然后給出的φ的定義：

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

其中?r是雷達距離分辨率。因為ΦDP是遞增的，KDP是非負的，為了保證KDP非負，引入了k2，其表示如下：

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

通過引入后項算子H2,φ'可以寫成與式（11）相同的形式：

觀測的差分相位與邊界條件的差值為：

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

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

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

圖4 變分擬合結(jié)果

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

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

3 實證分析

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

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

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

其中，ZH(mm6/m3)是水平反射率因子，ZDRl=10ZDR/10是線性尺度上的差分反射率因子。本文使用經(jīng)過變分重構(gòu)后的KDP數(shù)據(jù)，基于Zhang等[9]的RQPE方法和經(jīng)驗公式進行降水估計，整個算法（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ǎn)品，同時將基于V-RQPE估測的降水量也進行重采樣轉(zhuǎn)換到0.1°×0.1°的網(wǎng)格場上，以方便比較降雨事件的空間分布。圖7為基于KRIG插值、VRQPE的24小時累計降水估測值圖像，V-RQPE在雷達中心位置可以很好估計出降水量，但是在雷達的北部估測效果與實際值相比偏低。

圖7 24小時累計降水觀測值與估測值

3.2 評估指標

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

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

3.3 累計降水估測效果

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

圖8 累計降水散點圖

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

3.4 多種QPE方法比較

為了比較基于變分擬合的降水估計方法的性能，本文采用了六種不同的降水估計方案進行對比，具體方案如表1所示。

表1 六種QPE方案配置表

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

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

圖9 六種QPE方法的1小時累計降水估計對比

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

表2 六種QPE方法24小時累計降水估計評估表

4 結(jié) 論

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

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

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

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