一種雙高斯模型濾波器的研究與應用

(安徽四創(chuàng)電子股份有限公司，安徽合肥 230088)

0 引 言

通常情況下，地物雜波是指雷達波束在正常傳播情況下探測到的地物回波。地物雜波存在于雷達的低仰角掃描數(shù)據(jù)或雷達較近的范圍，對于任一特定仰角，典型的地物雜波污染從一個體掃到下一個體掃很少有變化，并且大多數(shù)時間都會出現(xiàn)。這嚴重影響雷達的數(shù)據(jù)估算，噪聲則會影響弱信號的檢測，從而導致參量估算的問題，因此，在信號處理系統(tǒng)中，濾除地物雜波能夠更好地進行氣象觀測。

在實際項目中可選擇全程濾波和動態(tài)雜波圖濾波。通過CMD(Clutter Mitigation Decision)算法識別雜波圖，根據(jù)雜波圖對雜波位置進行濾波，減少計算量，實現(xiàn)實時處理。在實際項目中設計了多個濾波器以適應復雜環(huán)境條件。

本文描述了IIR(Infinite Impulse Response) 橢圓濾波器、自適應高斯頻域濾波器(GMAP)和雙高斯濾波器(BGMAP)的原理，對GMAP濾波器和BGMAP濾波器進行了詳細的分析和比較，并使用實際天氣雷達數(shù)據(jù)對IIR橢圓濾波器和GMAP濾波器進行了測試和比較，結果表明：BGMAP濾波器在分析數(shù)據(jù)時效果比GMAP濾波器效果好，但耗費時間長，不能滿足實時處理要求；GMAP濾波器濾波效果優(yōu)于IIR濾波器，能滿足實時處理要求。

1 IIR橢圓濾波器

通常情況下，采用的是IIR橢圓濾波器[1]進行地物濾波，因為不需要很大的存儲器和計算量，在工程上比較容易實現(xiàn)。本文選擇的是3階極點和3階零點的IIR濾波器，濾波器系數(shù)按照重復頻率Fs= 300∶50∶2 000等35個頻率，凹口寬度按0.2 m/s步進進行濾波器設置，以滿足不同重復頻率、不同凹口寬度要求的應用需要。

IIR濾波器存在的主要問題是：

1) IIR濾波器有暫態(tài)響應時間，濾波器凹口寬度越窄，輸出暫態(tài)響應時間就越長；反之就比較短。暫態(tài)響應時間長意味著在一個波束寬度內(nèi)需要有比較多的回波脈沖，否則難以達到所設計的穩(wěn)態(tài)響應。

2) 濾波器凹口寬度如果設置凹口過寬，會對氣象回波造成損失；如果設置凹口過窄，濾波后地物剩余仍很多。

3) 會對速度較小的氣象回波或者折疊到零頻附近的氣象回波造成影響。

2 CMD算法

GMAP濾波和BGMAP濾波算法擬合需要一定時間，應當配合動態(tài)雜波圖進行濾波，動態(tài)雜波圖使用CMD雜波識別方法[2-3]，處理流程如下：

1) 檢查信噪比SNR，如果SNR<3 dB,該距離庫為噪聲，不作任何處理。

2) 計算3個特征量TDBZ，SPIN，CPA，其中TDBZ選擇以當前距離庫為中心的9個距離庫，SPIN選擇11個庫。

3)CPA5點中值濾波。

4) 通過隸屬函數(shù)將3個特征量映射到0-1區(qū)間，如式(7)、式(8)、式(9)所示。

5) 使用模糊邏輯組合SPIN和TDBZ。

6) 計算地物概率CP，如式(10)所示。

7) 將地物概率超過0.5的距離庫標識為地物。

8) 對地物標識CF進行平滑和填充。

9) 對標識的距離庫進行濾波。

10) 對濾波后的距離庫重新計算譜數(shù)據(jù)，包括反射率、速度和譜寬。

反射率紋理(TDBZ)是相鄰距離庫的反射率的差平方均值，測量相鄰距離反射率的變化，反映了反射率的平滑程度。天氣信號是大范圍的平滑信號，按式(1)計算：

(1)

雜波相位陣列校準值CPA(Clutter Phase Alignment)用于判斷由一組雷達脈沖回波的相位(I和Q)的穩(wěn)定性。CPA定義所有回波IQ時序的矢量和的幅度與所有回波IQ幅度的計算和的比值。

(2)

式中，xi=Ii+jQi。

SPIN反映了反射率因子在徑向梯度的變化，表示一定距離內(nèi)反射率符號變化的頻率。當符號相反，dBZi距離點SPIN數(shù)值的增加需要滿足

sign{dBZi+1-dBZi}=-sign{dBZi-dBZi-1}

(3)

(4)

時：

MSPINi=1

(5)

(6)

式中，spin_thres表示反射率因子門限，典型值為5 dBz。最后，將SPIN數(shù)值除以計算單元數(shù)，再乘以100。

(7)

(8)

(9)

(10)

3 GMAP濾波算法

當前國內(nèi)的主流濾地物雜波算法有IIR濾波及自適應譜濾波，這兩種方案對零頻附近存在天氣信號時，都會造成譜矩估算的偏差。采用有限積累點數(shù)進行譜矩估計，相當于對信號加矩形窗，當?shù)匚镫s波強度較強時，矩形窗對旁瓣的抑制小(-13 dB)，從雜波中泄露出的旁瓣功率會把弱的天氣信號淹沒。

GMAP算法[4-5]是在2004年由SIGMET公司的兩位工程師Siggia和Passarelli提出的，該算法有以下優(yōu)點：

1) 在零頻附近存在天氣信號與雜波疊加的情況下，可在濾除雜波的基礎上，基本恢復天氣信號，使得譜矩估計更為準確；

2) 在對地物雜波信號的譜寬估計更為準確，可以更好且自適應地估量出剔除地物雜波時所需凹口的寬度(傳統(tǒng)的IIR濾波需手動切換選擇凹口大小)；

3) 可以根據(jù)不同的雜信比(CSR)，自適應選擇所加窗函數(shù)，從而能對旁瓣達到很好的抑制，又不至于使得主瓣過寬影響氣象目標的譜寬估計。

4 BGMAP濾波算法

BGMAP算法在頻域內(nèi)去除地物回波，其基本假設是氣象回波和地物回波的功率譜為高斯分布[5]。雷達信號的功率譜由3個部分組成：氣象回波、地物雜波和噪聲，其模型表示如下：

(11)

式中：第一項為氣象回波模型，Pw為氣象回波功率，vrw為平均多普勒速度，σvw為譜寬；第二項為地物雜波模型，Pc為地物雜波功率，vrc為平均多普勒速度，σvc為譜寬，并且σvc<σvw；第三項為噪聲模型，Pn為噪聲功率，vN為Nyquist速度。

如果假定功率譜的每條譜線是指數(shù)分布的，則有

(12)

(13)

因此，需要求出滿足高斯模型的回波功率P、平均多普勒速度vr、譜寬σ使得式(13)最小。

在BGMAP中，使用標準的非線性擬合算法來擬合功率譜。考慮到氣象回波頻譜不一定都滿足高斯分布，因此重建的頻譜盡可能保留原來的頻譜而只去掉地物的部分。偏振雷達的參量不能夠直接從BGMAP的結果中計算出來。重新構建不含地物回波的IQ數(shù)據(jù)按照式(18)～式(22)進行處理。BGMAP算法流程圖[8]如圖1所示。

圖1 BGMAP流程圖

1) 選擇窗函數(shù)、FFT求頻譜以及估算噪聲

首先分別對IQ數(shù)據(jù)加海明窗，經(jīng)過FFT變換計算功率譜。如式(18)所示，原始頻譜的幅度為A(v)、相位為φ(v)。利用Hildebrand和Sekhon于1974年提出的方法計算噪聲功率。通過式(15)計算噪聲信號均值的平方，式(16)計算信號平方均值減去信號均值的平方。通過式(17)求出噪聲比值R，從而確定噪聲線，低于噪聲線的為噪聲功率譜，高于噪聲線的為信號和雜波譜。

(14)

var(Sn)=〈Sn〉2

(15)

var(Sn)=〈Sn2〉-〈Sn〉2

(16)

(17)

2) 擬合高斯氣象譜和雜波譜

按照雙高斯模型去掉原始頻譜的地物成分，得到新的幅度A′(v)和相位φ′(v)。當氣象回波的速度小于地物回波的速度時，氣象回波的相位會被地物回波干擾，導致與相位相關的參量估算問題，使用隨機相位來代替被地物干擾的那些相位，如式(21)所示。這種方式可以有效減少地物干擾。當新的幅度和相位確定后，通過傅里葉逆變換(IFFT)得到新的IQ數(shù)據(jù)。

FFT[s(n)]=Fs(v)=A(v)ejφ(v)

(18)

F′s(v)=A′(v)ejφ′(v)

(19)

(20)

(21)

s′(n)=IFFT[F′s(v)]

(22)

根據(jù)雙高斯擬合的結果，vc的值對應地物功率小于氣象回波和噪聲功率之和的位置。

BGMAP的實時性會受到兩個問題影響：首先，在收斂之前非線性的遞歸運算需要大量的循環(huán)；其次，如果頻譜不能用雙高斯模型來表示，遞歸可能不收斂，使用單高斯氣象模型擬合。單高斯模型和雙高斯模型的主要區(qū)別就在于是否對存在地物干擾的那些頻點進行擬合，如式(23)所示：

(23)

單高斯擬合方法基于S′(v)，該頻譜只包含氣象回波和噪聲成分。

不含地物回波的IQ數(shù)據(jù)重建過程與雙高斯濾波的重建過程類似，只是新的幅度A′(v)按照式(24)計算：

(24)

3) 檢驗雜信比(CSR)是否滿足窗函數(shù)標準

如果CSR>40 dB，加布萊克曼窗重復BGMAP濾波器的濾波過程和動態(tài)噪聲計算；如果CSR>20 dB，加布萊克曼窗重復BGMAP濾波器的濾波過程；如果CSR>25 dB，就取布萊克曼窗的結果；如果CSR<2.5 dB，加矩形窗重復BGMAP濾波器的濾波過程；如果CSR<1 dB，就取矩形窗的結果；其他情況輸出海明窗的結果。

5 濾波效果分析

圖2 GMAP分析結果

圖3 BGMAP分析結果

分析氣象信號，對比GMAP濾波器和BGMAP濾波器效果。在有氣象和強地物的信號時，GMAP濾波器識別如圖2所示，BGMAP濾波器識別如圖3所示。GMAP濾波器能夠識別出雜波點并擬合出氣象信號，并計算出速度和譜寬。BGMAP濾波器能擬合出雜波信號和氣象信號，并計算出速度和譜寬等信息。二者計算出的速度和譜寬相差不大，但時間差200倍左右。

只有氣象信號時，GMAP分析結果如圖4所示，BGMAP分析結果如圖5所示，二者雖然錯誤地識別了雜波信號，但都能擬合出氣象信號，并計算出速度和譜寬，時間同樣相差200倍左右。

圖4 只有氣象信號時GMAP分析結果

圖5 只有氣象信號時BGMAP分析結果

圖6是對GMAP和BGMAP進行1 000個距離庫運算統(tǒng)計的時間圖，紅線顯示BGMAP與GMAP耗費時間比值，可看出BGMAP計算花費的時間是GMAP的2個數(shù)量級，在工程應用中實時性差。

本文的實際回波數(shù)據(jù)是從安徽四創(chuàng)電子股份有限公司的多普勒天氣雷達獲取。

圖7顯示的是未濾波的多普勒四屏圖，顯示了濾波前后反射率因子、濾波后速度V和譜寬W。圖8顯示了使用IIR濾波器后的反射率因子速度譜寬圖。從圖中可以看出，在零速線上，天氣信號與雜波疊加的情況下，IIR濾波氣象回波損失較大。圖9中GMAP濾波器能識別氣象信號和雜波信號，恢復氣象信號，濾波效果優(yōu)于IIR濾波器。

圖6 GMAP與BGMAP花費時間統(tǒng)計

圖7 未濾波的強度速度譜寬圖

圖8 IIR濾波的強度速度譜寬圖

圖9 GMAP算法強度速度譜寬圖

6 結束語

本文對IIR濾波器、GMAP濾波器和BGMAP濾波器進行了分析和處理，結果表明：在非實時分析時，BGMAP的效果要優(yōu)于GMAP，能更好地識別出氣象信號和雜波信號；BGMAP濾波器和GMAP濾波器的效果優(yōu)于IIR濾波器，由于BGMAP濾波器計算時間要比GMAP計算時間多很多，在工程實踐中，不適合使用BGMAP，會造成資源不夠用。

鑒于本文只是對雙高斯地物濾波器算法的初步研究，還未通過長期實際回波驗證，在今后的項目中應作進一步的分析和研究。