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

劉陳帥，張阿思，陳生

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

1 引 言

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

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

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

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

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

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

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

2.1 質量控制

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

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

2.2 物理約束

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

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

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

圖1 基于物理約束的ΦDP重構

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

2.3 邊界條件

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

圖2 ΦDP觀測數據質量控制示意圖

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

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

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

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

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

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

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

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

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

圖3 邊界條件計算

2.4 代價函數

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

然后給出的φ的定義：

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

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

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

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

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

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

通過代價函數和關于k的偏導數，使用牛頓迭代法進行求解，最終重構ΦDP[15]。

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

圖4 變分擬合結果

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

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

3 實證分析

3.1 數據與降水估計方法

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

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

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

圖6 V-RQPE算法流程圖

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

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

3.2 評估指標

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

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

3.3 累計降水估測效果

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

圖8 累計降水散點圖

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

3.4 多種QPE方法比較

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

表1 六種QPE方案配置表

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

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

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

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

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

4 結 論

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

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

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

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