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

0 引 言

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

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

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

1 IIR橢圓濾波器

通常情況下，采用的是IIR橢圓濾波器[1]進(jìn)行地物濾波，因?yàn)椴恍枰艽蟮拇鎯ζ骱陀?jì)算量，在工程上比較容易實(shí)現(xiàn)。本文選擇的是3階極點(diǎn)和3階零點(diǎn)的IIR濾波器，濾波器系數(shù)按照重復(fù)頻率Fs= 300∶50∶2 000等35個頻率，凹口寬度按0.2 m/s步進(jìn)進(jìn)行濾波器設(shè)置，以滿足不同重復(fù)頻率、不同凹口寬度要求的應(yīng)用需要。

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

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

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

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

2 CMD算法

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

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

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

3)CPA5點(diǎn)中值濾波。

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

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

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

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

8) 對地物標(biāo)識CF進(jìn)行平滑和填充。

9) 對標(biāo)識的距離庫進(jìn)行濾波。

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

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

(1)

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

(2)

式中，xi=Ii+jQi。

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

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

(3)

(4)

時：

MSPINi=1

(5)

(6)

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

(7)

(8)

(9)

(10)

3 GMAP濾波算法

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

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

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

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

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

4 BGMAP濾波算法

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

(11)

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

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

(12)

(13)

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

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

圖1 BGMAP流程圖

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

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

(14)

var(Sn)=〈Sn〉2

(15)

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

(16)

(17)

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

按照雙高斯模型去掉原始頻譜的地物成分，得到新的幅度A′(v)和相位φ′(v)。當(dāng)氣象回波的速度小于地物回波的速度時，氣象回波的相位會被地物回波干擾，導(dǎo)致與相位相關(guān)的參量估算問題，使用隨機(jī)相位來代替被地物干擾的那些相位，如式(21)所示。這種方式可以有效減少地物干擾。當(dāng)新的幅度和相位確定后，通過傅里葉逆變換(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ù)雙高斯擬合的結(jié)果，vc的值對應(yīng)地物功率小于氣象回波和噪聲功率之和的位置。

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

(23)

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

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

(24)

3) 檢驗(yàn)雜信比(CSR)是否滿足窗函數(shù)標(biāo)準(zhǔn)

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

5 濾波效果分析

圖2 GMAP分析結(jié)果

圖3 BGMAP分析結(jié)果

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

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

圖4 只有氣象信號時GMAP分析結(jié)果

圖5 只有氣象信號時BGMAP分析結(jié)果

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

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

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

圖6 GMAP與BGMAP花費(fèi)時間統(tǒng)計(jì)

圖7 未濾波的強(qiáng)度速度譜寬圖

圖8 IIR濾波的強(qiáng)度速度譜寬圖

圖9 GMAP算法強(qiáng)度速度譜寬圖

6 結(jié)束語

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

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