何煒琨，高 麗，王曉亮

（中國(guó)民航大學(xué)天津市智能信號(hào)與圖像處理重點(diǎn)實(shí)驗(yàn)室，天津 300300）

作為可再生能源，風(fēng)能已受到各國(guó)青睞。十多年來(lái)，世界風(fēng)力累計(jì)裝機(jī)容量呈現(xiàn)指數(shù)化增長(zhǎng)趨勢(shì)[1]。雷達(dá)波束覆蓋區(qū)域內(nèi)風(fēng)輪機(jī)后向散射引起的雜波與氣象目標(biāo)回波具有相似特性，進(jìn)而影響氣象目標(biāo)參數(shù)估計(jì)的穩(wěn)定性。研究氣象雷達(dá)風(fēng)輪機(jī)雜波抑制技術(shù)，對(duì)于降低氣象目標(biāo)的誤檢測(cè)和誤識(shí)別概率、提高氣象觀測(cè)預(yù)報(bào)準(zhǔn)確性、保障飛行安全具有重要意義。

目前氣象雷達(dá)風(fēng)電場(chǎng)雜波抑制技術(shù)主要圍繞凝視和掃描兩種模式展開(kāi)，凝視模式下的抑制方法主要包括Radon變換法、中值濾波法、匹配追蹤算法以及利用風(fēng)輪機(jī)雜波時(shí)域周期特性來(lái)實(shí)現(xiàn)風(fēng)電場(chǎng)雜波的抑制[2-4]。上述抑制方法要求雷達(dá)駐留時(shí)間較長(zhǎng)，對(duì)于相干脈沖間隔較少的氣象雷達(dá)而言，性能受到一定的局限。對(duì)于掃描模式下的雜波抑制方法：Kong等[5]提出了自適應(yīng)譜處理算法；Feng等[6]提出了距離多普勒回歸算法（RDR,range Doppler regression）；奧克拉荷馬大學(xué)氣象研究中心利用散射儀遙測(cè)系統(tǒng)基于自適應(yīng)算法抑制風(fēng)電場(chǎng)雜波[7]；Kong等[8-9]提出了基于自適應(yīng)濾波器（如維納濾波器等）的風(fēng)電場(chǎng)雜波抑制方法和最大后驗(yàn)概率（MAP,maximum a posterior）算法；Frank等[10]提出了信號(hào)分離法。上述算法或是未污染區(qū)域的數(shù)據(jù)樣本有限，影響運(yùn)算性能，或基于二次回波數(shù)據(jù)（Level II數(shù)據(jù)）進(jìn)行，運(yùn)算復(fù)雜度較高。

在傳統(tǒng)利用風(fēng)輪機(jī)雜波周期特性實(shí)現(xiàn)航管監(jiān)視雷達(dá)風(fēng)電場(chǎng)雜波抑制方法基礎(chǔ)上展開(kāi)研究[11]。該方法利用風(fēng)輪機(jī)雜波的周期特性得到隨頻率變化的校正曲線，修正回波功率譜，從而在保留飛機(jī)目標(biāo)的前提下能有效抑制風(fēng)輪機(jī)雜波。然而與飛機(jī)目標(biāo)不同，氣象目標(biāo)屬于分布式目標(biāo)，具有一定頻帶寬度，傳統(tǒng)基于周期特性的風(fēng)電場(chǎng)雜波抑制算法對(duì)于與氣象目標(biāo)處于相同頻帶范圍內(nèi)的風(fēng)輪機(jī)雜波抑制效果有限。利用氣象目標(biāo)參數(shù)隨距離連續(xù)分布特性，利用臨近未受風(fēng)輪機(jī)雜波污染區(qū)域的回波數(shù)據(jù)作為先驗(yàn)信息，對(duì)傳統(tǒng)基于周期特性的風(fēng)輪機(jī)雜波抑制后的結(jié)果進(jìn)行處理，以此提高算法性能。

1 數(shù)據(jù)獲取

1.1 氣象目標(biāo)回波仿真

氣象目標(biāo)的無(wú)規(guī)則運(yùn)動(dòng)使得其頻譜分布在一定連續(xù)頻譜區(qū)域內(nèi)，即頻譜擴(kuò)展。氣象目標(biāo)回波功率譜密度近似高斯分布為

其中：Pr、fd、σf分別為氣象目標(biāo)回波的平均功率、平均多普勒頻率及譜寬。

在一定信噪比（SNR）下，接收到的氣象目標(biāo)回波功率譜密度[12]為

其中：S（f）為式（1）求得的氣象目標(biāo)功率譜密度；N 為噪聲的功率譜密度，X為服從[0，1]均勻分布的隨機(jī)數(shù)。

接收目標(biāo)回波信號(hào)為

其中：θ（f）為模擬回波在[0，2π]均勻分布的相位譜；IFT為傅立葉逆變換。

1.2 氣象雷達(dá)回波信號(hào)模型

基于散射點(diǎn)疊加模型得到風(fēng)輪機(jī)回波模型[13]。該模型中風(fēng)輪機(jī)葉片回波等效為大量圓片回波的疊加，且輪機(jī)艙、桅桿回波為點(diǎn)目標(biāo)回波的疊加。其中，風(fēng)輪機(jī)葉片與雷達(dá)波束的位置關(guān)系[13]如圖1所示。

由此，風(fēng)輪機(jī)雜波回波信號(hào)為

其中：ρtower、ρblade分別為桅桿和葉片上各散射點(diǎn)的雷達(dá)散射強(qiáng)度；c為光速；λ、T分別表示雷達(dá)波長(zhǎng)與脈沖重復(fù)時(shí)間；Rn，i（t）為第n張葉片上第i個(gè)薄片和雷達(dá)間的徑向距離；u（t）為矩形脈沖信號(hào)；R為風(fēng)輪機(jī)軸心和雷達(dá)間的徑向距離；M為一個(gè)距離門中雷達(dá)回波數(shù)據(jù)的長(zhǎng)度。進(jìn)而，氣象雷達(dá)回波為

圖1 風(fēng)輪機(jī)葉片和雷達(dá)波束的位置關(guān)系Fig.1 Position relation between wind turbine blade and radar beam

2 基于周期特性的氣象雷達(dá)風(fēng)輪機(jī)雜波抑制

風(fēng)電場(chǎng)雜波對(duì)氣象雷達(dá)會(huì)產(chǎn)生顯著影響，氣象雷達(dá)風(fēng)輪機(jī)雜波抑制技術(shù)的研究很有意義。美國(guó)Lone Star風(fēng)電場(chǎng)對(duì)KDYX產(chǎn)生多徑干擾[1]，如圖2所示，產(chǎn)生的多徑散射可能使周圍氣象雷達(dá)信號(hào)受到影響，使氣象雷達(dá)產(chǎn)生大量的虛假目標(biāo)。圖3為平均反射率因子在Lone Star建立前后的對(duì)比圖。從圖3可看出，風(fēng)電場(chǎng)建立后其區(qū)域附近反射率因子數(shù)值變大很多，這會(huì)嚴(yán)重影響對(duì)氣象目標(biāo)的準(zhǔn)確觀測(cè)。

圖2 美國(guó)德克薩斯州Lone Star風(fēng)電場(chǎng)對(duì)附近氣象雷達(dá)產(chǎn)生的多徑干擾Fig.2 Multi-path interference caused by Lone Star Wind Farm on nearby weather radars in Texas

圖3 美國(guó)某地氣象雷達(dá)在附近風(fēng)電場(chǎng)建立前后所獲取的平均反射率Fig.3 Average reflectivities obtained by one meteorological radar before and after establishment of a nearby wind farm in US

掃描模式下的氣象雷達(dá)回波數(shù)據(jù)相當(dāng)于聚束模式下的雷達(dá)回波數(shù)據(jù)采樣。氣象雷達(dá)的相干脈沖數(shù)較少，為此每次掃描只能獲得很少的風(fēng)輪機(jī)回波，為實(shí)現(xiàn)基于風(fēng)輪機(jī)周期特性的風(fēng)輪機(jī)雜波抑制，首先需要對(duì)掃描模式下的風(fēng)輪機(jī)回波缺省數(shù)據(jù)進(jìn)行填充，而后利用周期特性抑制風(fēng)輪機(jī)雜波。

2.1 基于GAPES的回波缺省數(shù)據(jù)重構(gòu)

為實(shí)現(xiàn)掃描模式下風(fēng)輪機(jī)回波缺省脈沖的重構(gòu)，首先需從獲得的多次風(fēng)輪機(jī)回波掃描數(shù)據(jù)中搜索出兩個(gè)與多普勒頻率相反的峰值數(shù)據(jù)，然后進(jìn)行風(fēng)輪機(jī)回波周期的估計(jì)，最后基于GAPES算法實(shí)現(xiàn)兩個(gè)峰值之間缺省數(shù)據(jù)的重構(gòu)。

GAPES方法重構(gòu)風(fēng)輪機(jī)回波信號(hào)問(wèn)題即是估計(jì)雷達(dá)回波信號(hào)幅度及其相位的問(wèn)題。GAPES算法是APES算法在數(shù)據(jù)缺損情況下的一個(gè)擴(kuò)展[14]，其基本思想是利用已知數(shù)據(jù)，采用非參數(shù)估計(jì)方法對(duì)缺省數(shù)據(jù)進(jìn)行幅度和相位估計(jì)，且該算法是在缺省數(shù)據(jù)的譜內(nèi)容類似于已有數(shù)據(jù)譜內(nèi)容的假設(shè)前提下進(jìn)行的，其實(shí)現(xiàn)框圖[11]如圖4所示。

圖4 基于GAPES的風(fēng)輪機(jī)回波數(shù)據(jù)重構(gòu)算法實(shí)現(xiàn)框圖Fig.4 Reconstruction flowchart of echo data based on GAPES wind turbine

2.2 基于周期特性的風(fēng)電場(chǎng)雜波抑制

獲得風(fēng)輪機(jī)完整周期的回波數(shù)據(jù)后，即可基于風(fēng)輪機(jī)雜波的周期特性進(jìn)行雜波抑制。風(fēng)輪機(jī)回波有周期性，而氣象目標(biāo)回波不具有周期性。該抑制方法正是利用這個(gè)特性在時(shí)頻域內(nèi)實(shí)現(xiàn)抑制的。

首先提取風(fēng)輪機(jī)所處距離單元的雷達(dá)回波數(shù)據(jù)并對(duì)其進(jìn)行時(shí)頻分析；然后以風(fēng)輪機(jī)掃描周期為分段周期將頻域內(nèi)的回波數(shù)據(jù)分成若干段，針對(duì)每個(gè)頻點(diǎn)依次計(jì)算任意兩子數(shù)據(jù)的相關(guān)系數(shù)并求其平均值，循環(huán)多次得到隨頻率變化的調(diào)制函數(shù)C（f）。最后用C（f）對(duì)雷達(dá)回波功率譜進(jìn)行校正，即可達(dá)到抑制風(fēng)機(jī)雜波、保留氣象目標(biāo)的目的?；谥芷谔匦缘娘L(fēng)電場(chǎng)雜波抑制算法流程[11]如圖5所示。

圖5 風(fēng)輪機(jī)雜波抑制流程Fig.5 Flowchart of wind turbine clutter suppression

基于周期特性的風(fēng)電場(chǎng)雜波抑制算法，其關(guān)鍵步驟是頻譜調(diào)制函數(shù)C（f）的獲取。由于風(fēng)輪機(jī)雜波具有周期特性，因而按照風(fēng)輪機(jī)回波的周期對(duì)回波數(shù)據(jù)進(jìn)行分段，并求任意兩段子數(shù)據(jù)的相關(guān)系數(shù)，在風(fēng)輪機(jī)頻譜范圍內(nèi)，其調(diào)制函數(shù)C（f）取值較大，接近于1。而氣象目標(biāo)頻譜范圍內(nèi)所表現(xiàn)出的C（f）值則較小，利用調(diào)制函數(shù)C（f）對(duì)回波功率譜進(jìn)行校正[11]，可得

其中：X（f）、X′（f）分別為調(diào)制前后的頻域內(nèi)幅值；α 為需通過(guò)經(jīng)驗(yàn)確定的調(diào)制因子；10-αC(f)為調(diào)制因子，呈遞減特性，使C（f）對(duì)由風(fēng)電場(chǎng)雜波周期性所引起的峰值經(jīng)校正后的幅度譜出現(xiàn)谷值，而C（f）對(duì)噪聲和目標(biāo)回波引起的谷值經(jīng)校正后的幅度譜出現(xiàn)峰值，凸顯出氣象目標(biāo)信號(hào)而使風(fēng)輪機(jī)雜波盡可能小，從而實(shí)現(xiàn)氣象雷達(dá)風(fēng)輪機(jī)雜波抑制。

2.3 改進(jìn)算法

與飛機(jī)目標(biāo)不同，氣象目標(biāo)屬于分布式目標(biāo)，具有群聚性和彌散性等特性。這使得在頻域內(nèi)，氣象目標(biāo)的頻率分布范圍比飛機(jī)目標(biāo)更廣，可能會(huì)出現(xiàn)風(fēng)輪機(jī)回波頻譜與氣象目標(biāo)頻譜部分重疊的問(wèn)題，影響該頻段范圍內(nèi)的風(fēng)輪機(jī)雜波抑制效果。為此需對(duì)傳統(tǒng)的基于周期特性的風(fēng)電場(chǎng)雜波抑制算法進(jìn)行改進(jìn)。

改進(jìn)算法中，首先基于周期特性進(jìn)行風(fēng)輪機(jī)雜波抑制；然后根據(jù)氣象目標(biāo)參數(shù)隨距離均勻分布的特性，以風(fēng)輪機(jī)附近未污染距離單元的氣象回波數(shù)據(jù)作為先驗(yàn)信息對(duì)風(fēng)輪機(jī)雜波與氣象目標(biāo)頻譜重疊的頻段范圍內(nèi)的抑制結(jié)果進(jìn)行處理，以提升算法性能，改進(jìn)算法對(duì)應(yīng)的流程如圖6所示。

圖6 改進(jìn)算法流程圖Fig.6 Flow chart of improved algorithm

風(fēng)輪機(jī)雜波與氣象目標(biāo)頻譜重疊區(qū)域數(shù)據(jù)處理過(guò)程如下所述。首先，對(duì)鄰近風(fēng)輪機(jī)距離單元（可取相比風(fēng)輪機(jī)距雷達(dá)更近的單元）且方位相同的一個(gè)未污染單元回波數(shù)據(jù)做時(shí)頻分析；然后對(duì)風(fēng)輪機(jī)單元經(jīng)傳統(tǒng)基于周期特性的風(fēng)電場(chǎng)雜波抑制后的時(shí)頻譜進(jìn)行處理，在時(shí)頻域內(nèi)分別求取風(fēng)輪機(jī)單元每個(gè)頻點(diǎn)在時(shí)間頻率采樣點(diǎn)數(shù)值與相應(yīng)頻點(diǎn)未污染距離單元在時(shí)間頻率采樣點(diǎn)數(shù)值的差值；最后對(duì)每個(gè)頻點(diǎn)在相同時(shí)間頻率采樣點(diǎn)的差值與閾值進(jìn)行比較。如果其差值大于閾值，則抑制后每個(gè)頻點(diǎn)在時(shí)間頻率采樣點(diǎn)的數(shù)值被該頻點(diǎn)鄰近未污染距離單元在時(shí)間頻率采樣點(diǎn)的數(shù)值取代；反之，則抑制后該頻點(diǎn)在時(shí)間頻率采樣點(diǎn)的數(shù)值不變。需要說(shuō)明的是，閾值的選取需要通過(guò)未污染區(qū)域時(shí)頻譜的變化作為先驗(yàn)信息，選取具體過(guò)程為：對(duì)與風(fēng)輪機(jī)方位相同、距離稍遠(yuǎn)和稍近的若干個(gè)鄰近單元作時(shí)頻分析，可求得任意兩個(gè)相鄰距離單元時(shí)頻譜在各時(shí)頻采樣點(diǎn)的差值，各時(shí)頻采樣點(diǎn)的閾值取該時(shí)頻采樣點(diǎn)多組差值中的最大值。

3 算法性能的定量評(píng)價(jià)指標(biāo)

從抑制前后回波數(shù)據(jù)時(shí)頻譜圖像的熵值和雜波抑制前后氣象目標(biāo)參數(shù)估計(jì)偏差兩方面來(lái)定量評(píng)估風(fēng)輪機(jī)雜波抑制效果。熵值代表圖像的能量信息。由于氣象目標(biāo)的反射率因子、徑向速度和譜寬等參量隨距離均勻分布[15-16]，因而認(rèn)為未受風(fēng)電場(chǎng)污染區(qū)域參量隨距離變化的圖像熵值較小，污染區(qū)域相應(yīng)參量所對(duì)應(yīng)的圖像熵值將有所增加。

4 實(shí)驗(yàn)結(jié)果與性能分析

仿真實(shí)驗(yàn)中，氣象目標(biāo)及風(fēng)輪機(jī)相關(guān)參數(shù)設(shè)置如表1所示。氣象目標(biāo)回波功率譜如圖7所示。信噪比SNR=25 dB，在第190個(gè)距離單元包含風(fēng)輪機(jī)雜波。從該距離單元得到風(fēng)輪機(jī)與氣象回波的距離多普勒譜，如圖8所示。提取風(fēng)輪機(jī)與氣象目標(biāo)同時(shí)存在的距離單元，對(duì)其進(jìn)行時(shí)頻分析，如圖9所示。按照風(fēng)輪機(jī)周期進(jìn)行分段處理，對(duì)任意兩段數(shù)據(jù)求相關(guān)，并對(duì)所有相關(guān)系數(shù)取平均，得到相關(guān)系數(shù)隨頻率變化的回波功率譜調(diào)制函數(shù)C（f），如圖10所示。由圖10可以看出，調(diào)制函數(shù)在氣象目標(biāo)處于-100～110 Hz頻率范圍內(nèi)相關(guān)系數(shù)值較小，在只含風(fēng)輪機(jī)回波的110～310 Hz和-310～-110 Hz頻率范圍內(nèi)相關(guān)系數(shù)值與1接近。

表1 仿真參數(shù)Tab.1 Simulation parameters

圖7 氣象目標(biāo)與噪聲回波的功率譜Fig.7 Power spectrum of weather target and noise echo

圖8 風(fēng)輪機(jī)與氣象回波的距離多普勒譜Fig.8 Distance Doppler spectrum of wind turbine and weather echo

調(diào)制函數(shù)C（f）對(duì)回波信號(hào)的頻譜進(jìn)行校正以抑制風(fēng)輪機(jī)雜波，抑制后的時(shí)頻分析如圖11所示?？梢钥闯?，處于氣象目標(biāo)頻率外的風(fēng)輪機(jī)雜波已得到完全抑制，但氣象目標(biāo)回波頻率范圍內(nèi)的風(fēng)輪機(jī)雜波仍存在。為抑制殘余雜波，利用氣象目標(biāo)均勻分布特性，用其鄰近距離單元的氣象目標(biāo)數(shù)據(jù)作為先驗(yàn)信息，對(duì)圖11抑制后的數(shù)據(jù)進(jìn)行處理，結(jié)果如圖12所示。

圖9 雜波抑制前的時(shí)頻分析Fig.9 Time-frequency analysis before clutter suppression

圖10 調(diào)制函數(shù)C（f）Fig.10 Modulation function C（f）

圖11 傳統(tǒng)雜波抑制算法處理后結(jié)果Fig.11 Processing result of conventional clutter suppression algorithm

圖12 改進(jìn)雜波抑制算法處理后的結(jié)果Fig.12 Processing result of improved clutter suppression algorithm

為進(jìn)一步驗(yàn)證算法的有效性，對(duì)抑制前后的回波數(shù)據(jù)利用脈沖對(duì)法分別進(jìn)行氣象目標(biāo)的參數(shù)（平均徑向速度和譜寬）估計(jì)，如表2所示。從表2可看出抑制后數(shù)據(jù)所估計(jì)的平均徑向速度為10.26 m/s、譜寬2.04 m/s更接近于真實(shí)氣象目標(biāo)的平均徑向速度和譜寬?；夭〞r(shí)頻譜熵值由抑制前的8.23 bit減小到1.76 bit。綜合來(lái)看，改進(jìn)后的基于風(fēng)輪機(jī)回波周期特性抑制算法能夠在保留氣象目標(biāo)信號(hào)的前提下有效抑制了風(fēng)輪機(jī)雜波。

表2 定量評(píng)估抑制效果Tab.2 Quantitative evaluation of suppression performence

5 結(jié)語(yǔ)

利用風(fēng)輪機(jī)雜波的周期特性以及氣象目標(biāo)參數(shù)隨距離連續(xù)分布的特性，實(shí)現(xiàn)氣象雷達(dá)風(fēng)輪機(jī)雜波的有效抑制，并利用抑制前后氣象目標(biāo)參數(shù)估計(jì)誤差及回波時(shí)頻譜圖像熵值的變化定量評(píng)價(jià)算法的性能，驗(yàn)證其有效性。主要討論了氣象雷達(dá)風(fēng)輪機(jī)雜波抑制，是風(fēng)電場(chǎng)雜波抑制后續(xù)研究的方向之一。風(fēng)電場(chǎng)雜波抑制方法的頻譜校正參數(shù)是依經(jīng)驗(yàn)取值的：取較小值時(shí)，風(fēng)輪機(jī)雜波有殘留；取較大值時(shí)，噪聲及部分氣象目標(biāo)也受影響。因此，后期需要找到更科學(xué)合理的參數(shù)值選取方法。