汪 玲 田 鳳* 朱岱寅 孟凡旺 吳 迪
①(南京航空航天大學(xué)雷達成像與微波光子技術(shù)教育部重點實驗室 南京 210016)
②(中國航空工業(yè)集團公司雷華電子技術(shù)研究所 無錫 214063)
氣象雷達是一種具有高時空分辨率的降水探測系統(tǒng),是大氣觀測必不可少的工具[1]。具有雙極化功能的多普勒氣象雷達[2]可以同時記錄降水的微物理和動態(tài)特征,極大提高了定量降水估計(Quantitative Precipitation Estimation, QPE)[3,4]性能。一般情況下,QPE應(yīng)用中使用雷達數(shù)據(jù)的前提是要有足夠的測量精度,而雷達數(shù)據(jù)往往會受到雜波的影響。重要的是,與強降水相比,弱降水在雜波去除時可能很容易地被去除。因此,研究雜波抑制方法,在有效濾除雜波的同時保留氣象信息一直是氣象雷達信號處理的關(guān)鍵。
氣象雷達雜波按多普勒速度不同主要分為靜止地物雜波和運動雜波,相應(yīng)地有針對各自特點的雜波抑制濾波器。高斯模型自適應(yīng)處理(Gaussian Model Adaptive Processing, GMAP)[5]是一種典型的地雜波抑制技術(shù),它可以解決氣象目標和地雜波在時域重疊的問題。自適應(yīng)雜波環(huán)境分析法(CLutter Environment ANalysis using Adaptive Processing, CLEAN-AP)[6]將雜波檢測與濾波統(tǒng)一在一起。與GMAP相比,由于使用幅度和相位確定適當?shù)陌伎趯挾萚7],CLEAN-AP具有更好的雜波抑制性能。文獻[8]利用雙極化氣象雷達在晴空模式下接收的地雜波數(shù)據(jù)訓(xùn)練模糊神經(jīng)網(wǎng)絡(luò),自適應(yīng)地計算地雜波各偏振參量隸屬函數(shù)的參數(shù),然后利用訓(xùn)練得到的地雜波隸屬函數(shù)對降水模式下的地雜波進行抑制。但是,以上方法均不能去除運動雜波。
地面氣象雷達探測遇到的運動雜波包括鳥類和昆蟲的回波、海雜波和各種無線電頻率干擾等。文獻[9–12]詳細分析了運動鳥類和昆蟲的回波特征。文獻[13]基于貝葉斯分類器,提出了一種海雜波識別方法。無線電頻率干擾[14]在雷達平面位置指示器(Plan Position Indicator, PPI)呈現(xiàn)為點狀、條帶狀干擾。由雷達系統(tǒng)本身引起的條帶狀雜波也會影響雷達數(shù)據(jù)的使用。大多數(shù)情況下,條帶狀雜波是PPI顯示畫面上沿某個方位向的連續(xù)距離門上的斑點,在頻域觀測時其多普勒速度不固定,傳統(tǒng)的雜波抑制方法無法進行有效抑制。這些條帶狀雜波會影響雷達數(shù)據(jù)質(zhì)量,降低雷達基產(chǎn)品及后續(xù)衍生產(chǎn)品的準確性。文獻[15,16]介紹了條帶狀雜波在國際電信傳輸和雷達研究中心研制的毛毛雨雷達(IRCTR Drizzle RAdar, IDRA)中的影響,由荷蘭Delft理工大學(xué)研制的MESEWI雷達系統(tǒng)[17]也存在這個問題。
為了能夠同時有效去除靜止地物雜波和運動雜波,并保留盡可能多的氣象信息,本文提出了一種基于譜極化參數(shù)(Spectral Polarimetric Parameters,SPP)的雜波濾波方法。根據(jù)氣象和雜波在距離-多普勒(Range Doppler, RD)域內(nèi)的特征不同進行前者的保留和后者的抑制。首先利用氣象目標的頻譜極化特征和連續(xù)性構(gòu)造譜極化參數(shù)并設(shè)置閾值,在RD域內(nèi)生成一個2元掩模,其中“1”對應(yīng)保留氣象,“0”對應(yīng)去除雜波。結(jié)合形態(tài)學(xué)方法,對2元掩模進行內(nèi)部空洞填充和邊界平滑,使被抑制的弱氣象信息得到恢復(fù)。然后采用面向?qū)ο蟮乃季S,運用泛洪填充算法將2元掩模標記為氣象對象掩模和雜波對象掩模。最后引入譜寬作為額外的參數(shù),利用雜波對象譜寬普遍小于氣象對象譜寬的特點,篩選出所有氣象對象掩模并疊加在一起得到SPP濾波器。SPP濾波器適用于同時發(fā)射同時接收(Simultaneous Transmission and Simultaneous Reception, STSR)模式和交替發(fā)射同時接收(Alternate Transmission and Simultaneous Reception, ATSR)模式的雙極化氣象雷達,具有通用性。且本文算法是逐個對雷達徑向數(shù)據(jù)進行處理的,具有較高的計算效率,可以實時實現(xiàn)。
本文各節(jié)內(nèi)容安排如下:第2節(jié)首先引入雙極化氣象雷達譜極化參數(shù),給出詳細的定義;第3節(jié)詳細闡述SPP濾波器的設(shè)計步驟;第4節(jié)討論SPP濾波器的參數(shù)選擇;第5節(jié)采用X波段和C波段雙極化氣象雷達實測數(shù)據(jù)對算法的性能進行驗證;第6節(jié)給出結(jié)論。
假設(shè)雙極化氣象雷達發(fā)射y 極化波、接收x 極化波,x ,y ∈{h,v}, 其中h 表示水平極化, v表示垂直極化。一次徑向掃描得到的回波信號為 Sxy(r,t),其中 r 表示距離向,采樣數(shù)為Nr, t表示駐留時間,采樣數(shù)為 NDwell。對Sxy(r,t)回波序列進行1維離散傅里葉變換得到Sxy(r,v)個距離點的回波強度,而譜反射率因子s Zxy(r,v)還包含了氣象目標的多普勒速度信息。
基于RD譜定義譜線性退極化比s LDR和譜相關(guān)系數(shù)s ρhv[18],如式(3)和式(4)所示
式(4)中 〈·〉表示對信號多普勒維進行1維高斯滑動平均操作。
由上述式(2)—式(4)定義可見,基于RD譜定義的譜極化參數(shù)綜合了目標的極化信息和多普勒頻率信息,還包括了目標的動態(tài)信息,更有利于區(qū)分氣象目標和非氣象目標。
為了能夠有效去除雜波,同時保留盡可能多的氣象信息,本文設(shè)計了一種基于譜極化參數(shù)的雜波濾波方法,SPP濾波器為RD域內(nèi)的一個2元掩模,它可以根據(jù)氣象和雜波在RD域內(nèi)的特征不同進行前者的保留和后者的抑制。SPP濾波器的設(shè)計包含以下4個部分。
(1) RD域掩模生成。設(shè)置閾值,對譜極化參數(shù)進行閾值濾波,得到初步的2元掩模。2元掩模作用于原始回波的RD譜數(shù)據(jù),它可以提取氣象信息,屏蔽雜波信息。雙極化氣象雷達通常為STSR或ATSR工作模式,STSR工作模式下只能得到Vhh和Vvv共極化電壓信號,ATSR工作模式下可以得到Vhh, Vvv共 極化電壓信號和Vhv, Vvh交叉極化電壓信號。所以按照雷達工作模式,分為兩種情況來處理:
當雷達工作模式為STSR時,對譜相關(guān)系數(shù)sρhv進行閾值濾波。假設(shè)獲得的2元掩模為Msρ,Msρ是 大小為Nr×NDwell的 矩陣,矩陣元素msρ∈Msρ取值為
其中,“1”對應(yīng)保留氣象目標,“0”對應(yīng)去除雜波。 Tsρ為譜相關(guān)系數(shù)閾值,可以根據(jù)雜波和氣象目標的特性設(shè)置。氣象目標的s ρhv值一般遠大于除地雜波以外的其他雜波的s ρhv值,而地雜波和氣象目標的s ρhv值分布區(qū)間大致相同,因此需要進一步采用其他方法抑制地雜波。
需要指出的是,掩模生成采用的閾值濾波處理可能會濾除弱氣象信息,并有部分雜波剩余,在后續(xù)處理中需要考慮如何改善。
(2) 基于形態(tài)學(xué)方法的氣象信息恢復(fù)。數(shù)學(xué)形態(tài)學(xué)在二值圖像中有廣泛的應(yīng)用[19],本文采用閉運算來恢復(fù)缺失的氣象目標信息。形態(tài)學(xué)閉運算可以填充二值圖像內(nèi)的小孔,彌合小裂縫,并保持總的位置和形狀不變。經(jīng)過算法第(1)步濾波后,氣象目標內(nèi)部回波強度比較弱的區(qū)域會有一些缺失,基于氣象目標在RD域內(nèi)具有連續(xù)性的特征,所以形態(tài)學(xué)閉運算可以填充這些缺失的氣象信息,而不會過度填充出現(xiàn)虛假氣象。閉運算是利用結(jié)構(gòu)元素對二值圖像先膨脹后腐蝕。由于圓盤對氣象邊界的平滑效果較好,本文將結(jié)構(gòu)元素選為半徑為R 的圓盤。
對 Mstep1進行閉運算后,獲得新的2元掩模Mstep2,該掩模增強了對弱氣象目標的保留。
(3) 標記對象,生成對象掩模。利用泛洪填充算法[20]將Mstep2中具有相同值“1”的連續(xù)區(qū)域標記為獨立的對象。對2元掩模Mstep2進行處理,未被標記的值為“1”的點和其相鄰區(qū)域?qū)⒈粯擞洖橐粋€獨立的對象,基于氣象目標在RD域的連續(xù)性,本文采用8-領(lǐng)域像素填充法進行標記。對RD域標記結(jié)束后,按照對象面積降序排列,得到對象集合O =[O1,O2,···,ON,···], O同時包含了氣象對象與雜波對象。每個對象Oi∈O, 都對應(yīng)一個大小為Nr×NDwell的2元對象掩模MOi。
一般情況下,氣象對象有兩個性質(zhì):對象面積大和對象個數(shù)有限性。因此,選擇對象集合O 中前 N個對象即包含了所有氣象目標,然后再進行下一步處理。
(4) 篩選出所有氣象對象。逐個提取對象集合O中的前 N個對象,經(jīng)過篩選后將其歸類為氣象對象或雜波對象。如果選擇的對象同時包含氣象和雜波,說明氣象目標和雜波此時具有相似的譜極化參數(shù)值和面積,形成的掩模不具有區(qū)分性。為了進一步去除雜波,需要利用譜寬值進一步篩選。
MSPP是作用于原始譜功率的2元掩模,用于去除雜波和噪聲,并保留盡可能多的氣象信息。
圖1給出了SPP雜波抑制濾波器設(shè)計的流程圖。由圖1可知,SPP濾波算法是對每一個徑向數(shù)據(jù)進行處理,在進行下一次徑向掃描的同時可處理上一次徑向掃描的數(shù)據(jù),能夠并行處理,具有較好的實時性。此外,SPP濾波算法一開始就在RD域生成一個初步的濾波掩模,后續(xù)的步驟都是對濾波掩模進行調(diào)整,整個濾波過程只有4步,且每一步都是簡單的“0”和“1”操作,計算復(fù)雜度較低。
在選擇譜極化參數(shù)的閾值時,完全保留氣象信息和去除所有雜波往往會有沖突。閾值 TsLDR和Tsρ是根據(jù)氣象目標和雜波去除的比例[15,21]選取的。除此之外,針對地面氣象雷達,為了抑制地物雜波,使用零多普勒速度附近的凹口濾波器生成MGC。根據(jù)經(jīng)驗數(shù)值,靜止地物雜波的多普勒速度主要分布在 –0.5~0.5 m/s。
利用形態(tài)學(xué)方法進行氣象信息恢復(fù)時,選擇半徑為 R的圓盤作為形態(tài)學(xué)閉運算的結(jié)構(gòu)元素。為了確定 R的大小,首先設(shè)定真實氣象區(qū)域,然后按照第3節(jié)所述的(1)和(2)對原始譜功率 sP(r,v)進行處理。
為了得到所有氣象對象掩模,根據(jù)氣象對象的兩個性質(zhì)(對象面積大和對象個數(shù)有限性),在這里只選擇對象集合 O 中前 N個對象對應(yīng)的掩模進行處理, N的選取準則為前 N個對象應(yīng)包含所有的氣象對象。計算前 N 個對象掩模MO1,MO2,···,MON中值為“1”的像素點總數(shù)Pobject,基于形態(tài)學(xué)進行氣象信息恢復(fù)后的2元掩模Mstep2中值為“1”的像素點總數(shù)為 Pstep2,當Pobject/Pstep2≥90%時就認為得到的前 N個對象包含了所有的氣象對象。根據(jù)這個原則,一般當N =8時即可滿足這個條件,也可以根據(jù)數(shù)據(jù)的不同進行微調(diào)。譜寬閾值Tobject為
其 中,R 為圓盤半徑大小,L 是雜波的最大譜寬。
圖1 SPP濾波器設(shè)計流程圖
為驗證本文所提雜波抑制方法的有效性,采用X波段和C波段氣象雷達數(shù)據(jù)進行處理,其中X波段氣象雷達工作模式為ATSR,有交叉極化測量,一個方位上的駐留脈沖數(shù) NDwell=512, C波段氣象雷達工作模式為STSR,不包含交叉極化測量,駐留脈沖數(shù) NDwell=43。為了說明本文方法在雜波抑制和氣象保留方面的優(yōu)良性能,在ATSR模式下選取MDsLDR濾波方法[15]與本文算法進行比較,在STSR模式下選取傳統(tǒng)的基于時域的門限因子雜波抑 制法[22]進行比較。
5.1.1 譜極化參數(shù)分析
由一個徑向掃描數(shù)據(jù)得到的譜極化參數(shù)如圖2所示,圖2(a)和圖2(b)分別為譜功率 sP和譜線性退極化比s LDRhh。
從圖2(a)可以看出該徑向掃描數(shù)據(jù)包含氣象目標、地雜波和條帶狀雜波。通過對多個徑向掃描數(shù)據(jù)分析,發(fā)現(xiàn)條帶狀雜波具有以下特點:(1)多普勒速度不為0;(2)在RD域內(nèi)隨機出現(xiàn),每個徑向觀察到的多普勒速度都不同;(3)強度和氣象目標強度相當。這些增加了雜波抑制難度。
由圖2(b)可見,氣象目標的s LDRhh值明顯小于雜 波的s LDRhh值。
5.1.2 濾波參數(shù)確定
根據(jù)氣象目標和雜波去除的比例[15,21],并結(jié)合圖3氣象和雜波的s LDR概率密度函數(shù)(Probability Density Function, PDF)分布圖,氣象目標的sLDR值主要分布在–7 dB以下,雜波的s LDR值都比較大主要分布在–7 dB以上,所以取TsLDR=?7 dB。
選取圖2(a)的數(shù)據(jù),根據(jù)式(9)計算得到的Pd和Pfa與圓盤半徑R 之間的關(guān)系如圖4所示。由圖4可知,當半徑 R在[1, 3]區(qū)間時,Pd增長較快,Pfa增長較慢;當 R 在[3, 7]區(qū)間時,情況相反。R=3時, Pd與Pfa之間的差值達到最大值,因此設(shè)置圓盤半徑。
根據(jù)數(shù)據(jù)統(tǒng)計,并結(jié)合圖2(a),該數(shù)據(jù)的雜波最 大譜寬L =5 ,根據(jù)式(11)計算得到Tobject=11。
5.1.3 雜波抑制結(jié)果分析
一個徑向掃描數(shù)據(jù)在RD域生成SPP濾波器的過程如圖5所示。
在圖5(a)中,經(jīng)過 sLDR閾值濾波后殘留了部分雜波點和條帶狀雜波,且氣象區(qū)域內(nèi)部和邊界均存在缺失點。運用形態(tài)學(xué)閉運算,從圖5(b)可以看出氣象區(qū)域得到了恢復(fù)。圖5(c)—圖5(g)為標記的前5個對象,由圖5(g)可知從第5個標記對象開始只包含了雜波對象。圖5(h)為SPP濾波掩模。
對圖2(a)的原始譜功率分別進行MDsLDR濾波和SPP濾波,得到圖6所示的結(jié)果。對比圖6(a)和圖6(b),不管是與氣象目標分開較好的地物雜波,還是與氣象目標重疊的條帶狀雜波,兩種濾波器都可以有效地去除。但是在3~5 km處,MDsLDR濾波算法在抑制條帶狀雜波同時造成了弱氣象目標的抑制,如圖中圓圈所示,所以在保持弱氣象方面SPP濾波器效果更好。
與5.1.2節(jié)類似,可以確定s ρhv閾值為Tsρ=0.6,
圖2 X波段氣象雷達數(shù)據(jù)在方位角為261.9°時的譜極化參數(shù)
圖3 氣象和雜波的s LDR概率密度分布圖
圖4 圓盤半徑大小與P d 和P fa的關(guān)系
圖5 X波段氣象雷達一個徑向掃描數(shù)據(jù)生成的濾波掩模
圖6 濾波處理后的譜功率sP(r,v)
地雜波凹口濾波長度LGC=3 ,圓盤半徑R =1,雜波最大譜寬 L=3,根據(jù)式(11)計算得到譜寬閾值Tobject=5。
針對一個徑向掃描數(shù)據(jù)按照本文算法進行處理得到如圖7所示結(jié)果。圖7(a)中白色區(qū)域代表雜波抑制后的氣象區(qū)域,紅色輪廓包圍區(qū)域代表真實氣象區(qū)域。
圖7(a)為生成的SPP濾波器,圖7(b)為原始譜功率,圖7(c)為SPP濾波后結(jié)果。圖7(c)中零多普勒速度附近的空白值表示抑制的地雜波,該寬度取決于地雜波的多普勒速度分布區(qū)間,因為地物是靜止目標,所以其速度一般分布在–0.5~0.5 m/s。氣象目標(下落的雨滴、雪花等)普遍具有多普勒速度,只有很少一部分分布在零多普勒速度附近,所以對氣象目標強度的影響可以微乎不計。對比圖7(b)和圖7(c)可以看出該濾波器可以有效地去除雜波,最終只保留氣象信息。
C波段和X波段氣象雷達數(shù)據(jù)經(jīng)過SPP濾波后的反射率因子如圖8和圖9所示。
圖8(b)與圖8(c)分別為基于時域的門限因子雜波抑制和SPP濾波處理后的反射率因子。對比圖8(a)的原始反射率,兩種方法均能夠有效地去除雜波。圖8(b)中黑色橢圓圈標注區(qū)域表示時域處理結(jié)果缺失了部分弱氣象信息,在保留弱氣象方面性能比SPP濾波器差。
圖7 C波段氣象雷達數(shù)據(jù)在方位角為244.1°時利用本文算法處理結(jié)果
圖8 C波段氣象雷達數(shù)據(jù)反射率因子估計
圖9 X波段氣象雷達數(shù)據(jù)反射率因子估計
圖9(b)與圖9(c)分別為X波段氣象雷達MDsLDR濾波和SPP濾波處理后的反射率因子,圖9(b)中黑色橢圓圈標注區(qū)域表示MDsLDR濾波結(jié)果缺失了部分弱氣象信息,由此可見,SPP濾波器在保留弱氣象信息方面的性能優(yōu)于MDsLDR濾波器。
針對雙極化氣象雷達中非氣象回波濾除問題,本文提出一種基于譜極化參數(shù)的雜波濾波方法。根據(jù)氣象和雜波在RD域內(nèi)的特征不同進行前者的保留和后者的抑制。首先利用氣象目標的頻譜極化特征和連續(xù)性構(gòu)造譜極化參數(shù),結(jié)合形態(tài)學(xué)方法,在RD域內(nèi)生成一個2元掩模。然后采用面向?qū)ο蟮乃季S,將2元掩模標記為氣象對象掩模和雜波對象掩模。最后引入譜寬值作為額外的參數(shù),篩選出所有氣象對象掩模并疊加得到SPP濾波器。對原始譜功率應(yīng)用SPP濾波器,可以保留氣象,去除雜波和噪聲。氣象雷達實測數(shù)據(jù)雜波抑制結(jié)果驗證了本文方法的有效性。結(jié)果表明,本文方法適用于不同工作模式的雙極化氣象雷達,與MDsLDR濾波器和基于時域的門限因子雜波抑制方法相比,具有較好的雜波抑制效果和氣象保留性能。SPP濾波器實現(xiàn)簡單,計算復(fù)雜度低,可以實時應(yīng)用于雙極化監(jiān)視氣象雷達,且為機載極化氣象雷達雜波抑制算法奠定了基礎(chǔ)。