雙偏振雷達差分傳播相移的五種濾波方法對比分析

魏 慶，胡志群，劉黎平

(1.成都信息工程學院電子工程學院，四川成都610225;2.中國氣象科學研究院災害天氣國家重點實驗室，北京100081)

0 引言

近年來，偏振技術在天氣雷達方面取得了巨大發展。雙線偏振多普勒天氣雷達能發射和接收水平與垂直的線偏振波，相對于普通的多普勒天氣雷達，除了能探測反射率因子Z外，還能測得差分反射率因子ZDR、差分傳播相移ΦDP、差分傳播相移率KDP及相關系數ρhv等偏振參量。這對提高雷達定量估測降水的精度以及識別降水粒子的相態和空間取向具有非常重要的作用[1]。

差分傳播相移率KDP在降水估測以及粒子相態識別等方面有著巨大潛力[2]，而KDP又是由ΦDP估算得到，因此研究ΦDP的預處理方法對降水估測的影響有重要意義。ΦDP定義為電磁波在非球形降水粒子區中傳播時，水平偏振和垂直偏振的傳播常數不同造成的相位差。Mueller算法定義:ΨDP=ΦDP+δ，ΨDP為雷達的總差分相移，δ為差分散射相移[3]。對于δ只與降水粒子的特性相關，滿足Rayleigh散射條件的降水粒子，δ約為0。但是，當有Mie散射影響時，δ就不能忽略，而且δ可能存在較大值，這就是δ效應。對于較強的雨區，降水粒子可能嚴重不滿足Rayleigh散射條件，ΨDP包含ΦDP和δ 2個分量。如果ΦDP為信號，那么δ則為噪聲，δ的值會隨著降水粒子增大而增大。因此在使用雙線偏振雷達資料之前，減少δ噪聲對ΦDP信號測量結果的影響是十分必要的[4]。除此之外，在實際探測過程中，氣象目標環境和氣象目標自身的擾動以及雷達系統本身的各種噪聲，都會使ΦDP受到干擾[5]。因此，對噪聲污染的ΦDP資料進行前期的數據預處理對于后期數據應用于雷達定量降水估測尤為重要。

Hubbert等[6]提出采用有限沖激響應(Finite Impulse Response，FIR)和無限沖激響應(Infinite Impulse Response，IIR)低通濾波的方法處理ΦDP徑向數據。Hubbert和Bringi[3]在該方法的基礎上，又提出一種迭代濾波方法，此方法比低通濾波方法的應用范圍更廣。在中國，對ΦDP數據，一般做滑動平均和中值濾波等數據預處理工作。曹俊武[7]采用平滑的方法處理ΦDP數據資料。何宇翔等[4]根據差傳播相移ΦDP的特點引入卡爾曼濾波方法。胡志群等[8]提出小波分析法處理差傳播相移ΦDP數據，即利用db5小波對ΦDP進行濾波處理。這些ΦDP濾波方法應用于降水估測方面的效果究竟如何?文中對比分析13點滑動平均、13點中值濾波、FIR濾波、卡爾曼濾波以及小波分析5種不同ΦDP數據處理方法的效果，及其在降水估測中的應用效果進行對比分析討論。

1 濾波方法介紹

對滑動平均、中值濾波、FIR濾波、卡爾曼濾波和小波分析5種不同ΦDP徑向數據處理方法進行簡單介紹。

1.1 滑動平均法

假設動態測試數據y(t)由確定性成分f(t)和隨機性成分e(t)組成，且f(t)為有效信號，e(t)為隨機測試誤差(即噪聲)。經離散采樣后，可以將動態測試數據表示為[9]

為了抑制隨機測試誤差{ei}的影響，對非平穩的數據{yi}作局部平均處理，以減小{ei}所造成的隨機起伏。這樣沿全長N個數據逐一小區間上進行不斷的平均處理，即可得出較平滑的有效信號{fi}，而去除隨機起伏的噪聲。

POLC雷達一個徑向上共有1000個距離庫，庫長為150m。根據曹俊武[10]分析雙線偏振雷達徑向數據平滑方法表明，隨著平滑庫數N的增大，平滑后曲線的變化更加平穩，資料隨機起伏的點減少，使雷達徑向數據資料更加連續、可靠;但是，隨著滑動平均庫數N的增加，平滑的效果也逐漸減弱，且會平滑掉部分有效的回波信號，從而再次造成誤差。因此，選用N=13，即采用13點滑動平均的方法對POLC雷達徑向數據ΦDP進行處理。

1.2 中值濾波法

中值濾波是一種非線性數字濾波技術，能夠有效地處理脈沖噪聲。中值濾波器的原理就是利用一個奇數點的移動窗口，將窗口的中心點值用窗口內的中值代替。假設一個序列a1，a2，…，an，取窗口長度為M(M為奇數)，對其進行中值濾波，即從序列中相繼抽出M個數ai－v，…，ai－1，ai，ai+1，…，ai+v(ai為窗口中心值，M=2v+1)，再將這M個數按照從小到大的順序排列，排在數列中心點的那個數即為中值濾波器的輸出值[11]。

基于與選擇平滑庫數N的同樣考慮，取中值濾波窗口M=13，即13點中值濾波方法處理POLC雷達徑向數據。

1.3 FIR 濾波法

設有限沖激響應FIR濾波器的單位沖激響應為h(n)，那么一個長度為M的FIR濾波器的系統函數為

它是次數為M－1的z－1的一個多項式。一個輸入信號x(n)，通過FIR濾波器輸出信號y(n)，可以表示為

在處理POLC雷達ΦDP徑向數據時，k為雷達徑向距離庫序數，M為FIR濾波窗口的大小。

1.4 卡爾曼濾波法

卡爾曼濾波器是一種最優化自回歸數據處理算法。以最小均方誤差為最佳估計準則，采用信號與噪聲的狀態空間模型，利用前一時刻的估計值和現在時刻的觀測值更新對狀態變量的估計，求出現在時刻的估計值。根據ΦDP信號的特性，r表示距離，假定a(r)為動態白噪聲向量，均值為零、方差為Q1的高斯白噪聲序列，p(a)～N(0，Q1)，且a(i)和a(j)(i≠j)互不相關，E{a(i)}=0，E{a(i)a(j)}=Q1ε(i－1)，其中 ε(i)滿足:

V(r)為觀測白噪聲向量，均值為零、方差為Q2的高斯白噪聲序列，p(V)～N(0，Q2)，a(r)與V(r)不相關。抽象為離散時間過程的狀態變量估計，定義觀測變量為Z(r)，用離散隨機差分方程表示[4]，則有

1.5 小波濾波法

胡志群等[8]根據雙線偏振雷達ΦDP的特性，提出利用小波分析的方法處理雙偏振雷達ΦDP數據。方法可以使ΦDP數據更加平滑和連續，使有效的氣象信息得到保留。小波分析處理過程大致分為以下3個步驟:

(1)信號分解。利用db5小波函數對原始ΦDP數據進行5層分解，每層均可分解成兩個信號分量，即ΦDP的近似信號(低頻分量)和ΦDP的細節信號(高頻分量)。

(2)信號去噪。根據固定閾值原則計算出閾值thr，并利用該閾值對細節信號進行量化處理，閾值函數如式6所示

式中cj，k為處理后的小波系數，sgn(·)為符號函數。

(3)信號重構。重構經小波分析后的ΦDP數據，采用小波分解的第5層近似信號和經過量化處理的第1～5層的細節信號。

2 個例分析

2.1 雷達資料介紹

中國氣象科學研究院與安徽四創電子股份有限公司聯合研發的車載移動式C波段雙線偏振雷達(POLC)，參加了2013年的“華南季風強降水”外場觀測試驗，觀測地點位于廣東鶴山氣象觀測新站(112.981°E，22.737°N)。表1為POLC雷達的主要技術指標。

3.2 ΦDP徑向距離廓線分析

理論上，ΦDP是一個積分效應量，其徑向距離廓線具有隨著距離的增加而緩慢上升的特點[12]。當雷達電磁波穿過相對較弱且分布均勻的降水區時，ΦDP通常隨著距離而緩慢的遞增;但是，相對于較強的降水回波而言，ΦDP在強回波區域呈現明顯的遞增現象。雷達在實際探測中，ΦDP資料受到非氣象目標、雷達系統噪聲和誤差的影響，使其徑向廓線存在波動起伏現象。因此，有必要對ΦDP資料做前期預處理工作。下面對13點滑動平均、13點中值濾波、FIR濾波、卡爾曼濾波以及小波分析5種不同方法對ΦDP徑向資料數據處理的效果進行對比分析。

圖1為POLC雷達于2013年5月8日17時12分(北京時)0.5°仰角、133°方位角差傳播相移ΦDP的原始資料和經5種不同濾波方法處理后的徑向距離廓線圖。從ΦDP原始數據(圖1a)不難發現，徑向ΦDP數據的連續性較好，且整體上具有隨距離增加而上升的變化趨勢，這與理論上ΦDP的變化特征相一致。但也明顯存在大的波動和毛刺現象。經濾波處理后(圖1b～圖1f)，ΦDP徑向距離廓線的整體變化趨勢與原始數據基本一致，對原始數據的毛刺現象都能起到不同程度的抑制作用。其中小波分析處理后的ΦDP徑向距離廓線(圖1f)平滑程度最好，卡爾曼濾波(圖1e)和FIR濾波(圖1d)次之，中值濾波(圖1c)和滑動平均(圖1b)處理的ΦDP徑向距離廓線波動較大。

表1 POLC雷達主要技術指標

下面定義波動指數FIX反應ΦDP徑向數據的波動情況，表達式如下

其中，i為距離庫序數;ΦiDP為對應庫的ΦDP值。FIX的值越大，說明徑向距離廓線的波動性越大。對圖1的ΦDP連續徑向距離庫(共719個)求出FIX值，如表2所示。由表2統計的結果可知:小波分析處理后的ΦDP數據資料效果最好，其他不同方法處理的效果由好到差的排列順序是:卡爾曼濾波、FIR濾波、滑動濾波、中值濾波。

圖1 ΦDP徑向距離廓線圖

表2 ΦDP徑向距離廓線波動指數統計

2.3 KDP距離廓線分析

差傳播相移率KDP可以反映橢球形雨滴的形態，且幾乎不受衰減的影響。POLC雷達直接觀測到的參量是ΦDP，根據KDP的定義，若降水區中相鄰距離rn與rn+1處測得的差分傳播相移分別為ΦDP(rn)和ΦDP(rn+1)，則[13]

在實際應用中，采用變距離法對KDP進行最小二乘法擬合，擬合的距離r(r=150*N，N為距離庫數，雷達探測距離庫長度為150m)根據對應距離庫的ZH平均值選擇，N取值如下[14]

這樣處理可以有效減小觀測資料擾動產生的誤差;對于較強的回波區域，短距離的擬合(即10×150 m=1.5 km)，能夠有效降低周圍的弱回波區對其影響，使強回波區的特點得以保留;對于弱降水，較長的距離平均能夠減少因信噪比(SNR)減弱引起的ΦDP噪聲增加的影響，從而提高KDP估測弱降水的穩定性與精度。

圖2為差傳播相移ΦDP的原始資料以及分別經上述5種方法數據處理后經變距離法擬合得到的KDP徑向距離廓線圖。理論上，KDP的值為正，由于ΦDP徑向數據的波動使估算的KDP存在負值。由圖2(a)可以看出，ΦDP原始數據經變距離法擬合得到的KDP存在大量的負值，且KDP徑向數據的波動性較大，存在較多的峰值(或谷值)。由圖2(b)～圖2(e)可以看出，KDP徑向數據的波動性得到一定的改善，但是仍然有大量的負值存在。再看圖2(f)，KDP徑向數據的整體波動趨勢基本一致，KDP負值基本不存在。這說明小波分析處理方法對存在較大波動的ΦDP數據處理效果較好，能夠有效抑制KDP的負值，且能保留數據的真實性。

圖2 KDP徑向距離廓線圖

3 不同濾波方法對QPE的影響

采用劉黎平等[15]根據數值模擬方法得到C波段雙線偏振雷達基于偏振參量KDP測量降水強度R的關系式

由于ΦDP數據存在波動起伏現象，尤其是在較低信噪比時，KDP往往存在負值，因此在采用式(7)反演降水時，KDP為負值的時次不參加計算與對比分析。

2013年5月8日POLC雷達探測范圍內發生持續性層狀云降水，且持續時間較長。選取09:00～20:00共12小時雷達降水資料定量估測降水。地面雨量計資料為雷達探測半徑70km內的小時降水量站點資料。以地面雨量計小時降水量作為實況降水量，以評估雷達降水估測的效果。由于每個雷達體掃時間在5～7分鐘，為了得到小時降水量，采用每個體掃的數據分別反演降水強度，再根據體掃時間權重后相加得到每小時降水量(即降水強度)。

統計雷達-雨量計對的相對誤差RE和相關系數ρ對雷達降水估測的效果進行評估，表達式為

其中Ri和Gi分別為雷達反演降水量和雨量計實測降水量，n為有效雷達-雨量計對的數量。相對誤差表示數據的可信度，RE值越小數據的可信度越高;相關系數表示數據變化趨勢的一致程度，ρ值越大表明線性一致程度越高。

表3統計了ΦDP分別經不同濾波方法處理后擬合的KDP應用于降水估測中雷達-雨量計對的相對誤差RE和相關系數ρ。由于KDP在降水較小時的測量誤差比較大，為避免其影響降水估測的結果，只對自動雨量計小時降水量大于1mm的雷達-雨量計對有效數據進行誤差統計。表中時次表示雷達-雨量計對有效數量。根據統計結果看，在每個雨強階段，ΦDP小波分析方法處理后擬合的KDP應用于降水估測的相對誤差都最小，相關系數均在0.7以上;隨著平均雨強的增大，相對誤差值減小，相關系數增大，說明小波分析方法對ΦDP數據的處理效果最好，應用于降水估測的精度最高;從表3可以看出，其他ΦDP數據處理方法應用于降水估測中的效果，由好到差依次為:卡爾曼濾波、FIR濾波、滑動平均、中值濾波。

表3 ΦDP的不同濾波處理方法應用于QPE的誤差統計

4 結束語

利用中國氣象科學研究院災害天氣國家重點實驗室的C波段雙線偏振雷達外場試驗觀測資料，對差分傳播相移ΦDP徑向數據采用5種不同濾波方法進行了徑向距離廓線、擬合的KDP徑向距離廓線以及應用于降水估測的效果進行了對比分析，并得到以下結論

(1)對比分析了5種不同ΦDP徑向數據濾波方法，發現FIR濾波、卡爾曼濾波和小波分析較傳統的滑動平均和中值濾波方法對波動的ΦDP徑向數據更能起到抑制作用，且能有效地保留氣象回波信息;進一步對比FIR濾波、卡爾曼濾波和小波分析3種較好效果的濾波方法，小波分析方法處理ΦDP徑向數據的效果最好。

(2)采用變距離法對5種不同濾波方法處理后的ΦDP數據進行最小二乘法擬合得到差傳播相移率KDP，其中小波分析處理后的ΦDP數據擬合的KDP負值最少。

(3)采用R(KDP)對5種方法得到的KDP進行降水估測評估，發現小波分析方法應用在降水估測中的精度最高;其他方法應用在降水估測中的效果，由好到差的順序依次為卡爾曼濾波、FIR濾波、滑動平均、中值濾波。

根據國內外的許多科研結論，雙線偏振雷達使用KDP進行降水估測能夠有效提高雷達定量估測降水的精度。而ΦDP數據的測量誤差直接影響KDP的估算值，因此，研究ΦDP數據的濾波方法對后期KDP數據的應用，如相態粒子識別和降水估測，起重要作用。

[1] Seliga T A，Bringi V N.Potential use of radar differential reflectivity measurements at orthogonal polarizations for measuring precipitation[J].J Appl Meteor，1976，15(1):69－76.

[2] Zrnic D S，Ryzhkov A V .Advantages of rain measurements using specific differential phase[J].J Atmos Oceanic Technol，1996，13(2):454－464.

[3] Hubbert J，Bringi V N.An iterative filtering technique for the analysis of copolar differential phase and dual-frequency radar measurements[J].J Atmos Ocean Technol，1995，12(3):643－648.

[4] 何宇翔，呂達仁，肖輝，等.X波段雙線極化雷達反射率的衰減訂正[J].大氣科學，2009，33(5):1027－1037.

[5] 杜牧云，劉黎平，胡志群，等.雙線偏振雷達差分傳播相移的小波濾波初探[J].暴雨災害，2012，31(3):248－254.

[6] Hubbert J，Chandrasekar V，Bringi V N，et al.Processing and interpretation of coherent dual-polarized radar measurements[J].J Atmos Ocean Technol，1993，10(2):155－164.

[7] 曹俊武，胡志群，陳曉輝，等.影響雙線偏振雷達相位探測精度的分析[J].高原氣象，2011，30(3):817－822.

[8] Hu Z Q，Liu L P.Applications of wavelet analysis in differential propagation phase shift data denoising[J].Adv Atmos Sci，2014，31(4):824－834.

[9] 胡廣書.數字信號處理理論算法與實現[M].北京:清華大學出版社，2003:1－649.

[10] 曹俊武.雙線偏振多普勒雷達的資料質量控制及降水粒子相態識別方法的研究[D].北京:中國氣象科學研究院，2006:1－120.

[11] 張明謙，李雷.改進的中值濾波算法[J].兵工測量與控制，2007，26(5):77－80.

[12] Bringi V N，Chandrasekar V，Balakrishnan N，et al.An examination of propagation effects in rainfall on radar measurements at microwave frequencies[J].J Atmos Oceanic Technol，1990，7:829－840.

[13] 畢永恒，劉錦麗，段樹，等.X波段雙線偏振氣象雷達反射率的衰減訂正[J].大氣科學，2012，36(3):495－506.

[14] Wang Y，Chandrasekar V.Algorithm for estimation of the specific differential phase[J].J Atmos Oceanic Technol，2009，26(12):2565－2578.

[15] 劉黎平，葛潤生，張沛源.雙線偏振多普勒天氣雷達遙測降水強度和液態含水量的方法和精度研究[J].大氣科學，2002，26(5):709－719.