一種雙高斯模型濾波器的研究與應(yīng)用

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

0 引 言

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

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

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

1 IIR橢圓濾波器

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

IIR濾波器存在的主要問(wèn)題是：

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

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

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

2 CMD算法

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

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

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

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

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

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

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

7) 將地物概率超過(guò)0.5的距離庫(kù)標(biāo)識(shí)為地物。

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

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

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

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

(1)

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

(2)

式中，xi=Ii+jQi。

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

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

(3)

(4)

時(shí)：

MSPINi=1

(5)

(6)

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

(7)

(8)

(9)

(10)

3 GMAP濾波算法

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

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

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

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

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

4 BGMAP濾波算法

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

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

圖1 BGMAP流程圖

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

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

(14)

var(Sn)=〈Sn〉2

(15)

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

(16)

(17)

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

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

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

(23)

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

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

(24)

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

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

5 濾波效果分析

圖2 GMAP分析結(jié)果

圖3 BGMAP分析結(jié)果

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

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

圖4 只有氣象信號(hào)時(shí)GMAP分析結(jié)果

圖5 只有氣象信號(hào)時(shí)BGMAP分析結(jié)果

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

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

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

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

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

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

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

6 結(jié)束語(yǔ)

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

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