一種雙極化氣象雷達自適應譜極化濾波方法

安孟昀 殷加鵬 黃建開 龐 晨 李永禎 王雪松

(國防科技大學電子信息系統(tǒng)復雜電磁環(huán)境效應國家重點實驗室 長沙 410073)

1 引言

氣象雷達可實現(xiàn)對大氣的高時空分辨率觀測，是大氣觀測必不可少的工具[1,2]。具有極化測量能力的多普勒氣象雷達獲得的信息被廣泛應用于定量降水估計[2—5]、短期天氣預報[6]和預警[7]等方面。多普勒和極化信息相結(jié)合可以反映采樣體積中分布目標的微物理和動力學特性，顯著提高定量降水估計性能。在氣象觀測中，利用雷達數(shù)據(jù)進行降水估計的前提是數(shù)據(jù)要有足夠的精度，但是由于雷達工作環(huán)境較復雜，雷達數(shù)據(jù)經(jīng)常受到雜波的影響[8]，例如地雜波、海雜波、生物回波、射頻干擾[9]等。這些雜波會影響氣象目標的反射率、多普勒速度和極化參量的測量與估計，進而影響降雨估計和微物理研究。因此，有效的雜波抑制方法研究十分重要[10]。

氣象雷達雜波根據(jù)多普勒速度的不同可以分為靜態(tài)雜波(即地雜波)和動態(tài)雜波。對于地雜波，傳統(tǒng)的方法是通過以零速度為中心的凹口濾波器來抑制，其性能取決于地雜波的譜寬[2,11]。為解決窄帶零速度氣象信號的損失問題，文獻[12]提出了高斯模型自適應處理方法(Gaussian Model Adaptive Processing,GMAP)。但是，當GMAP作用于未被雜波影響的距離單元時，降雨目標反射率的估計值會降低。因此，在應用這些雜波抑制方法之前需要先檢測雜波是否存在。常用的檢測算法有雜波抑制決策(Clutter Mitigation Decision,CMD)算法[13,14]和譜雜波識別(Spectrum Clutter Identification,SCI)[15]算法。在實際中無論是CMD算法還是SCI算法都會造成降雨目標的漏判和誤判。文獻[16]提出了一種模糊算法識別并去除地雜波，這種方法在提高地雜波識別率的同時也會造成降雨目標的損失。為了將雜波檢測和濾波集成到一種算法中，文獻[17]提出自適應雜波環(huán)境分析法(Clutter Environment Analysis Using Adaptive Processing,CLEANAP)，它利用幅度和相位來改變凹口濾波器的凹口寬度，使氣象目標估計的偏差更小。為了在濾波的同時保留較完整的降雨目標，文獻[18]提出一種時域回歸濾波技術(shù)，該技術(shù)提高了對降雨目標的估計性能，但是該方法還未得到實測數(shù)據(jù)的驗證。

除了上面提到的雜波，帶狀雜波也會影響氣象雷達數(shù)據(jù)的精度[10]。帶狀雜波是由雷達系統(tǒng)本身或外部輻射源造成的。大多情況下，這些雜波并不局限于某些距離單元，且在距離-多普勒(Range Doppler,RD)域中是非平穩(wěn)的，這使得傳統(tǒng)的雜波抑制方法無法濾除此類雜波[10,19]。目前，荷蘭代爾夫特理工大學研制的IDRA (IRCTR Drizzlc Radar)雷達的帶狀雜波抑制處理是在RD圖(PPI的每個方位向)上進行的，它由以0 ms—1為中心的窄凹口濾波器和雙譜線性去極化比濾波器(Double spectral Linear Depolarization Ratio,DsLDR)組成[20]。但是由于降雨和雜波的譜線性退極化比(spectral Linear Depolarization Ratio，sLDR)有所重疊，DsLDR濾波器無法將它們完全分開。此外，凹口濾波器在去除噪聲的同時會損失弱降雨目標[10]，導致剩余的降雨區(qū)域不連續(xù)。為了解決上述問題，Yin等人[10]提出移動雙譜線性去極化比(Moving Double spectral Linear Depolarization Ratio,MDsLDR)濾波器，該濾波器利用降雨在RD圖上連續(xù)分布這一特征，用數(shù)學形態(tài)學方法補償缺失的降雨。雖然sLDR可以用來區(qū)分降雨、噪聲與雜波，但是它很容易受到噪聲的影響，引入測量誤差[20,21]，且該方法只適用于全極化雷達。為了解決這一問題，Yin等人[22]提出移動譜去極化比(Moving spectral Depolarization Ratio,MsDR)濾波器，該濾波器引入更穩(wěn)健的sDR參量。MsDR濾波器還可解決無交叉極化測量能力的雙極化雷達的雜波抑制問題。

MDsLDR濾波器和MsDR濾波器都運用了譜極化濾波技術(shù)，利用降雨目標和非降雨目標的譜極化特性在RD圖上的連通特性差異，在盡可能多地保留降雨目標的同時濾除非降雨目標。譜極化濾波技術(shù)是根據(jù)信噪比與對應濾波參數(shù)的關(guān)系，以及雜波和降雨的去除比例來確定濾波閾值，利用此閾值在RD圖上產(chǎn)生一個 {0,1}的二值濾波值，作用于原始的功率譜，實現(xiàn)降雨目標的保留和雜波的濾除[21,22]。但是，在實際PPI掃描中，不同方位向上的降雨情況存在差異，且雜波和噪聲依賴該方位向的地理環(huán)境。為了獲得較好的濾波效果和降雨估計，根據(jù)每個方位向回波情況選取閾值顯得尤為重要。

本文提出的自適應移動譜去極化比(Adaptive Moving spectral Depolarization Ratio,AMsDR)濾波方法可以根據(jù)每個方位向回波情況以及氣象目標和雜波在RD圖上的差異進行氣象目標的保留與雜波的濾除。首先，為了解決雙極化雷達的雜波抑制問題，本文根據(jù)氣象目標的譜極化特性構(gòu)造sDR參數(shù)。其次，本文利用氣象目標和雜波的詹森-香農(nóng)(Jensen-Shannon,JS)散度確定各徑向掃描的濾波閾值，在RD圖上生成一個二元掩模。然后，結(jié)合數(shù)學形態(tài)學方法恢復降雨，最終實現(xiàn)自適應譜極化濾波。本文提出的自適應譜極化濾波器可以廣泛應用于雙極化氣象雷達以及全極化氣象雷達。

本文的結(jié)構(gòu)如下：第2節(jié)主要介紹極化多普勒氣象雷達的典型參量；第3節(jié)闡述AMsDR濾波器的設計步驟，并提出逐個方位向確定閾值的算法；第4節(jié)應用實測數(shù)據(jù)對濾波器的濾波效果進行驗證；最后對本文工作進行總結(jié)。

2 譜極化參量

譜極化技術(shù)利用多普勒和極化信息綜合表征氣象目標的微物理和運動特征，有利于反演氣象目標和抑制非氣象目標回波。下面介紹幾種典型極化氣象雷達測量參量[21]。

假設極化氣象雷達發(fā)射x極化電磁波、接收y極化電磁波，其中x,y ∈{h,v}，h 表示水平極化，v表示垂直極化。基于后向散射的約束，與距離r和多普勒速度v相關(guān)的譜反射功率定義為

其中，C為雷達常數(shù)，sSxy(r,v)表示二維的復距離-多普勒譜，sPxy(r,v)表示譜功率。定義譜差分反射率 sZdr(r,v)、譜共極化相關(guān)系數(shù)sρco(r,v)、譜去極化比s DR(r,v)為[21,22]

其中，sρco(r,v)中的〈〉表示在距離維或多普勒速度維的滑動平均，在本文中為后者，*表示復共軛。利用RD圖進行譜極化濾波之后，雜波和噪聲被濾除的同時保留降雨目標，可以更準確地估計測量參數(shù)。

3 AMsDR濾波器設計

3.1 AMsDR濾波器濾波步驟

利用降雨目標和非降雨目標的譜極化特性和在RD圖上連通特性的差異，譜極化濾波技術(shù)可以在盡可能多保留降雨目標的同時濾除非降雨目標。其中，MsDR濾波器提出的初衷是為了解決大部分同時發(fā)射同時接收(Simultaneous Transmission and Simultaneous Reception,STSR)模式的業(yè)務氣象雷達的雜波濾除問題。該濾波器主要分為4個步驟，第1步利用譜去極化比濾波器濾波，第2步是沿著多普勒維使用一個動態(tài)窗濾除帶狀雜波，第3步是利用二維的動態(tài)窗恢復去除的降雨濾除剩余雜波；最后，用數(shù)學形態(tài)學方法進一步恢復降雨信號。

為了能夠在有效去除雜波的同時保留降雨目標，本文提出一種AMsDR濾波方法。該方法可以根據(jù)每個方位向的數(shù)據(jù)情況自適應地進行譜極化濾波。AMsDR濾波器的設計主要包括以下6個部分，如圖1所示。

圖1 AMsDR濾波器流程圖Fig.1 Flow chart of the AMsDR filter

步驟1 自適應閾值選取。根據(jù)降雨目標與非降雨目標的譜極化參量在統(tǒng)計特性上的差異，利用JS散度確定譜極化濾波的閾值。其中閾值選取方法將在3.2節(jié)用真實的雷達數(shù)據(jù)討論。

步驟2 譜極化濾波。利用步驟1中求得的閾值進行譜極化濾波，得到{ 0,1}二元掩模，其中“0”表示雜波，“1”表示降雨目標。該二元掩模作用于原始的RD功率譜以提取降雨目標并濾除雜波[10]。

步驟3 一維多普勒滑動窗選取降雨目標?；趯铍s波譜寬的分析，設置長度為M的多普勒移動窗口，將中心單元以及其前后的M/2個單元視為一個整體。將移動窗口應用于步驟2中獲得的掩模，當移動窗口在多普勒譜中從第1個單元移動到最后一個單元時，如果移動窗口的任何位置有0，那么將中心單元替換為0，否則保持1。通過此過程，獲得濾波后的掩模[10]。

步驟4 二維滑動窗恢復降雨并濾除雜波。步驟3的處理幾乎會濾除所有的帶狀雜波，但是同時會損失一些降雨。為了進一步濾除剩余的雜波并恢復降雨，對步驟3得到的濾波后的掩模應用一個3×3大小的移動二維窗口[10]。將移動窗口的中心對齊到選擇的RD單元，得到此RD單元鄰近的9個值[M1,M2,...,M9]，然后對這9個元素求和得到

通過設置一個閾值TRD來決定該RD單元的濾波掩模MRD，即

其中，閾值TRD是根據(jù)降水和雜波的去除比例，并結(jié)合降雨目標濾波前后的最小均方誤差決定的，具體選取方法見文獻[10]。

通過此移動二維窗口，可以進一步去除步驟3中殘留的雜波，并恢復降雨區(qū)域附近的點。

步驟5 數(shù)學形態(tài)學法恢復濾除的降雨目標。降雨目標區(qū)域內(nèi)以及邊界的點可能會被濾除，可以利用數(shù)學形態(tài)學的形態(tài)閉合來恢復缺失的降雨[10,23]。將結(jié)構(gòu)單元設置為半徑為5的平面圓盤。通過此過程，獲得恢復降雨目標之后的掩模。

步驟6 圖像填充恢復降雨目標。經(jīng)過步驟5的處理會恢復大部分的降雨邊界，但是降雨區(qū)域內(nèi)部還會存在一些空洞。為了進一步恢復降雨目標，對步驟5得到的二值圖像進行填充。當空洞封閉時，圖像填充會將該空洞填充為周圍背景的值。

下面從香農(nóng)熵的角度闡述自適應閾值選取的原理。

3.2 閾值自適應選取

信息熵(香農(nóng)熵)常被用來作為一個系統(tǒng)信息含量的量化指標，可以用來作為系統(tǒng)方程優(yōu)化的目標或者參數(shù)選擇的判據(jù)[24]。JS散度是一個基于信息熵的用于比較兩個分布之間的相似性或差異性的參數(shù)[25]，其定義如下：

其中，P(X)和Q(X)是兩個不同分布的歸一化概率密度函數(shù)，0≤JS(P||Q)≤1,J S(P||Q)=0表示兩個概率分布完全相同，JS(P||Q)=1表示兩個概率分布完全不同[26]。在譜極化濾波中，X表示譜極化參量之一，即

雜波與無降雨情況下X的概率分布相似，降雨目標與無降雨情況下X的概率分布差異較大。利用這一特點，本文把無降雨情況下X的概率分布作為基準，分別衡量雜波、降雨目標與無降雨情況下X概率分布的差異。

為了防止單一組數(shù)據(jù)計算概率分布造成的特殊性和偶然性，本文利用多組無降雨情況下的數(shù)據(jù)來計算每個方位向上X的平均概率分布，在保持地雜波分布的情況下降低偶然出現(xiàn)的運動雜波的影響。

平均概率分布求解方法為：首先，按照方位角序列依次求解多組數(shù)據(jù)在每個方位向上X的概率分布；其次，根據(jù)多組數(shù)據(jù)的概率分布依次求得每個方位向上平均概率分布為

利用JS散度確定在方位角α上譜極化濾波的閾值。

其中，m ask(X,α)=0 表示在方位角α上非降雨目標部分，JSN(T)表示非降雨目標部分與無降雨目標情況下X分布之間的散度；J SP(T)表示降雨目標部分與無降雨目標情況下X分布之間的散度。

圖2為在方位角α上確定最優(yōu)閾值的算法。該算法輸入是X閾值的搜索范圍[Tmin,Tmax]，此搜索范圍由降雨目標在X上的分布確定。ΔT是閾值的步進值，考慮到計算效率和濾波效果，此步進值取(Tmax-Tmin)/100。在利用JS散度最大限度地提取降雨回波時，對于一個最優(yōu)的閾值T使得J SN(T)較小，JSP(T) 較大，當J SD(T)達到最大值時的閾值T為最優(yōu)的閾值(Topt)。求解步驟具體為：首先，根據(jù)輸入的閾值求J SD(T)；其次，求J SD(T)的最大值點，對JSD(T)求導，如果在搜索范圍內(nèi)存在一個閾值T使得，那么此閾值即為J SD(T)的極大值點，取此閾值為最優(yōu)的閾值；如果不存在使得的點，即在此方位角中沒有降雨目標，為了更好地濾除非降雨目標，取Topt為J SD(T)的最大值點。

圖2 閾值選取算法流程圖Fig.2 Flow chart of the algorithm for threshold selection

4 sDR濾波及結(jié)果分析

4.1 雷達系統(tǒng)與雜波特性分析

將JS散度最優(yōu)濾波閾值選取方法應用于AMsDR濾波器中，并與MsDR濾波器進行比較，使用IDRA雷達的數(shù)據(jù)對所提的濾波器進行性能驗證。IDRA雷達是X波段高分辨極化-多普勒氣象雷達。該雷達的測量模式為交替發(fā)射同時接收(Alternate Transmission and Simultaneous Reception,ATSR)模式，其具體參數(shù)如表1所示。

表1 IDRA雷達參數(shù)Tab.1 IDRA radar specifications

利用測量于2016年3月9日00:00—2016年3月13日12:00(UTC時間，每隔12小時測一組數(shù)據(jù))的無降雨數(shù)據(jù)計算非降雨部分的平均概率分布。以2016年3月9日00:00的數(shù)據(jù)為例，其原始雷達PPI如圖3(a)所示。使用測量于2017年4月26日00:00(UTC時間)含降雨的數(shù)據(jù)進行濾波性能測試，其原始雷達PPI如圖3(b)所示。

使用如圖3(b)所示的數(shù)據(jù)進行處理，其方位角為316.9°(圖3(b)黑色虛線)時譜極化參量如圖4所示。

圖3 原始反射率Fig.3 Raw reflectivity

由圖4(a)可見該方位向數(shù)據(jù)中除了含有氣象目標外還包含帶狀雜波和地雜波。通過多個方位向數(shù)據(jù)的統(tǒng)計分析，地雜波主要集中在零多普勒處，位置相對固定。帶狀雜波主要有以下幾個特點：(1)多普勒速度不為零；(2)在RD圖上隨機出現(xiàn)，帶狀雜波在不同方位向出現(xiàn)的位置都不相同；(3)帶狀雜波強度與降雨目標相當。帶狀雜波的這些特征增加了雜波抑制的難度[21]。

由圖4(b)可見，降雨的 sZdr為正值。帶狀雜波的 sZdr時正時負，并且降雨與周圍環(huán)境的sZdr差別不大。因此，使用 sZdr來區(qū)分降雨和雜波并不合適。大多數(shù)雙極化雷達系統(tǒng)可測量得到譜共極化相關(guān)系數(shù) sρco。大多數(shù)降雨目標的sρco非常接近于1，非降雨目標的 sρco遠低于1。圖4(c)表明，地雜波、帶狀雜波和降雨的sρco相 似。這意味著僅使用sρco不足以區(qū)分降雨，應該結(jié)合其他技術(shù)進一步去除這些雜波[20]。

由圖4(d)可見，降雨與雜波以及噪聲的 sDR值差異較大，這表明 sDR在抑制雜波方面具有較好的能力。對于沒有交叉極化測量能力的雙極化氣象雷達，sDR是有效的濾波極化參量。

基于上述分析，本文選取s DR參數(shù)進行譜極化濾波。下面分別從RD和PPI的角度對AMsDR的性能進行說明。

4.2 譜極化濾波

4.2.1 濾波閾值確定

選取圖4的數(shù)據(jù)，根據(jù)式(10)和式(11)計算JSP散度及 JSN散度和濾波閾值T之間的關(guān)系如圖5(a)所示，根據(jù)式(12)計算 JSD散 度與濾波閾值T的關(guān)系如圖5(b)所示，其中T=-9.6 是J SD的極值點。由圖5可知，當閾值T在[-15，-9.6]區(qū) 間時，J SP散度增長較快，JSN散度增長較慢。說明隨著閾值T的增加，保留的降雨目標變多，這一部分的分布特性與晴空條件下的分布特性差異變大；當閾值T在[-9.6,-5]區(qū)間時，情況相反。因此，在方位角為316.9°時取Topt=-9.6 dB。

圖4 IDRA雷達在方位角為316.9°時的譜極化參數(shù)Fig.4 The spectral polarimetric variables of IDRA radar at an azimuth angle of 316.9°

圖5 JS散度與閾值T的關(guān)系Fig.5 Relationship between the JS divergence and the threshold T

4.2.2 譜極化濾波分析

選取圖4的原始譜功率進行自適應譜極化濾波，得到圖6(a)和圖6(b)所示的結(jié)果。圖6(a)是使用MsDR濾波器濾波后得到的結(jié)果，此時的濾波閾值為—10 dB，該閾值是根據(jù)降雨目標和雜波去除的比例[10,21]，并結(jié)合降雨和雜波的概率密度函數(shù)分布確定的[21,26]；圖6(b)是使用本文提到的自適應閾值選取的算法，得到濾波閾值為—9.6 dB進行濾波的結(jié)果。兩種濾波方法在方位角為344.7°時的對比如圖6(c)和圖6(d)所示，此時濾波閾值為—9.2 dB。

圖6 濾波之后的譜功率Fig.6 Spectral power after filtering

對比圖6(a)和圖6(b)，兩種濾波閾值都可以有效地去除與降雨目標分開較好的地雜波以及與降雨目標重疊的帶狀雜波。在5～8 km處，MsDR濾波會造成氣象目標區(qū)域內(nèi)部和邊緣的降雨目標缺失(圖中圓圈所示)，利用AMsDR方法濾波情況得到改善，但在7～8 km處也會存在一個較小的空洞。這里的空洞是由速度模糊造成的，速度模糊導致降雨區(qū)域不連續(xù)，圖像填充算法無法將降雨目標填充。對比圖6(c)和圖6(d)，方位角為344.7°時AMsDR濾波方法在保留弱降雨目標方面效果優(yōu)于MsDR。

4.3 PPI處理與結(jié)果分析

經(jīng)過以上譜極化濾波的步驟，大部分的帶狀雜波會被濾除。利用降雨的空間連續(xù)性，本文在PPI維度濾除殘留的雜波。此濾波步驟作用于譜極化濾波之后的PPI功率圖，通過產(chǎn)生一個{ 0,1}的二值濾波掩模，使用譜極化濾波步驟4，利用式(5)和式(6)實現(xiàn)降雨目標的保留并進一步濾除雜波。

應用RD譜極化濾波以及PPI濾波處理后的雷達反射率在全掃描角度下的PPI顯示結(jié)果，如圖7所示。

圖7(a)是利用MsDR濾波器處理的PPI結(jié)果，圖7(b)是利用文中提出的AMsDR濾波方法以及PPI維濾波得到的PPI結(jié)果。從圖中可以看出，利用MsDR濾波存在大量散點且降雨區(qū)域不連續(xù)。相比之下，本文的濾波方法散點較少且降雨區(qū)域較連續(xù)。

圖7 濾波處理后的反射率Fig.7 Reflectivity after filtering

利用2016年3月4日00:00與2017年4月26日12:00(UTC時間)的數(shù)據(jù)對自適應濾波器的性能進行驗證。經(jīng)過MsDR以及AMsDR濾波之后的反射率因子如圖8和圖9所示。對比圖8(b)與圖8(c)以及圖9(b)與圖9(c)可知，無論是晴空還是有部分降雨，在濾除雜波和保留降雨方面，AMsDR濾波器的性能優(yōu)于MsDR濾波器。

圖8 2016年3月4日00:00 (UTC時間)的雷達數(shù)據(jù)不同處理后的反射率因子Fig.8 Reflectivity after different processing of the radar data at 00:00,4 March 2016 (UTC time)

圖9 2017年4月26日12:00 (UTC時間)的雷達數(shù)據(jù)不同處理后的反射率因子Fig.9 Reflectivity after different processing of the radar data at 12:00,26 April 2017 (UTC time)

5 結(jié)論

針對雙極化雷達中非氣象回波的濾除問題，本文提出了一種AMsDR濾波方法。該方法根據(jù)氣象目標和雜波在RD圖上的譜極化特征和區(qū)域連通性的特性差異在保留氣象目標的同時濾除雜波。首先根據(jù)徑向掃描數(shù)據(jù)確定該徑向濾波閾值，在RD圖上生成一個二元掩模?；诿嫦?qū)ο蟮乃枷?，將二元掩模標記為氣象對象掩模和雜波對象掩模。然后利用數(shù)學形態(tài)學以及圖像填充方法恢復降水，最終實現(xiàn)自適應譜極化濾波。本文對原始功率譜應用自適應濾波方法，可以更好地在保留氣象目標的同時去除雜波和噪聲。最后，通過實測氣象雷達數(shù)據(jù)的雜波濾除結(jié)果驗證了本文方法的有效性。結(jié)果表明，所提方法與MDsLDR濾波器相比具有更加廣泛的應用，與固定閾值的MsDR濾波器相比有更好的雜波抑制與降雨保留性能。此外，該方法實現(xiàn)簡單，計算復雜度較低。