吳仁彪 丁 紅 何煒琨 王曉亮 張 鑫

(中國民航大學(xué)天津市智能信號與圖像處理重點實驗室，天津 300300)

1 引言

隨著化石能源不斷耗盡，風(fēng)能作為一種清潔環(huán)保、可再生的新能源呈現(xiàn)出最迅猛的發(fā)展勢頭，據(jù)全球風(fēng)能理事會(Global Wind Energy Council， GWEC)報告顯示，近十多年以來，全球風(fēng)力發(fā)電累計裝機容量呈現(xiàn)飛速增長趨勢，2017年全球新增裝機容量超過50 GW[1]。根據(jù)近五年數(shù)據(jù)預(yù)測，到2021年全球新增裝機容量每年將達到75 GW左右的水平，到2021年底累計裝機容量將達到800 GW[2]。我國風(fēng)力資源極為豐富，風(fēng)力發(fā)電很可能在今后能源產(chǎn)業(yè)中起到領(lǐng)軍作用。然而，隨著風(fēng)電場規(guī)模在不斷擴建，其產(chǎn)生的雜波影響了氣象目標(biāo)參數(shù)估計的準(zhǔn)確性和穩(wěn)健性，造成對雷雨風(fēng)暴、湍流、風(fēng)切變等惡劣天氣的誤識別和誤判，對氣象雷達系統(tǒng)性能帶來愈發(fā)嚴(yán)重的影響。

關(guān)于氣象雷達風(fēng)電場雜波抑制問題，國外相關(guān)學(xué)者已關(guān)注多年。Kong Fanxing等人根據(jù)氣象目標(biāo)與風(fēng)輪機的頻譜特性差異利用自適應(yīng)譜處理方法抑制風(fēng)電場雜波[3]，該方法基于估計臨近非污染區(qū)域的氣象目標(biāo)參數(shù)自適應(yīng)地生成帶通濾波器，濾除污染區(qū)域的風(fēng)電場雜波，同時保留有用的氣象信息。Nai Feng等人提出的距離多普勒自回歸算法[4]則是首先檢測出風(fēng)電場雜波，用未受污染的數(shù)據(jù)來估計污染區(qū)數(shù)據(jù)，該方法是將抑制問題看成一個估計問題，而不是直接抑制污染區(qū)的風(fēng)輪機雜波，這類方法的有效性依賴于未受污染的數(shù)據(jù)及其與污染區(qū)的歐式距離。Beauchamp R M等人利用雷達回波的時域相關(guān)性，將風(fēng)輪機雜波近似為周期循環(huán)平穩(wěn)過程，地雜波和氣象目標(biāo)近似為廣義平穩(wěn)過程，利用自適應(yīng)滑動平均濾波器濾除風(fēng)電場雜波[5]，該方法適用于凝式模式的氣象雷達風(fēng)輪機雜波抑制。2016年, faruk U.等人根據(jù)氣象目標(biāo)和風(fēng)輪機在不同域中具有不同的稀疏性，采用信號分離理論從受污染的雷達回波數(shù)據(jù)中分離出風(fēng)電場雜波實現(xiàn)雜波抑制[6]，實驗中某些參數(shù)根據(jù)經(jīng)驗選取，算法較為復(fù)雜且在信噪比較低時抑制性能有限。國內(nèi)關(guān)于風(fēng)輪機雜波抑制，主要還是在航管監(jiān)視雷達和沿海超高頻雷達等方面研究，針對氣象雷達風(fēng)輪機雜波抑制我國目前研究相對較少，更多集中在風(fēng)電場雜波檢測及其對氣象雷達的影響評估等方面研究[7- 8]。因此，在我國風(fēng)能飛速發(fā)展的情況下，研究掃描模式的氣象雷達風(fēng)電場雜波抑制技術(shù)將有利于解決其對氣象雷達的影響，降低擬建的氣象觀測設(shè)備及風(fēng)電場的選址要求。

本文研究了氣象雷達風(fēng)電場雜波抑制的匹配追蹤(Matching Pursuit， MP)算法。該方法首先需要構(gòu)造超完備字典，往往規(guī)模較大的字典可以顯著提高算法性能，但同時也會極大地增加存儲量和計算量。為了解決這個問題，本文提出改進的匹配追蹤算法，首先根據(jù)大地主題反解方法獲得風(fēng)輪機與雷達的相對距離和方位先驗信息，提取出每個風(fēng)輪機對應(yīng)位置的雷達回波數(shù)據(jù)，從而可以降低待處理數(shù)據(jù)量及字典規(guī)模。實驗結(jié)果表明，改進的匹配追蹤算法在保證算法性能的前提下可以有效地降低計算量，提高運算效率。

2 數(shù)據(jù)模型

氣象目標(biāo)是由大量散射粒子組成，氣象粒子在一定范圍內(nèi)作無規(guī)則運動，因而在某一頻譜范圍內(nèi)會產(chǎn)生多普勒頻移。理論研究和實驗觀測表明氣象目標(biāo)的功率譜可以近似為高斯分布[9-11]，相位譜服從[0，2π]上的均勻分布，結(jié)合功率譜和相位譜信息，并對其進行離散傅里葉逆變換，其對應(yīng)的時域信號可以寫為[9]

(1)

其中I、Q分別是雷達回波的同向通道和正交通道時域序列，Sk為功率譜，θk為相位譜，N為頻域序列的長度。

風(fēng)輪機信號模型主要是依據(jù)散射點疊加原理，遠場條件下，可以將風(fēng)輪機看作是一系列散射點疊加而成的幾何體[12-13]。首先以風(fēng)輪機的輪機艙為坐標(biāo)中心，風(fēng)輪機的旋轉(zhuǎn)面為yoz平面，則x軸垂直于旋轉(zhuǎn)面，如圖1建立風(fēng)輪機與氣象雷達幾何關(guān)系圖。

圖1 風(fēng)輪機和氣象雷達的幾何關(guān)系圖

圖1中，φ為雷達視線(Line of Sight， LOS)與風(fēng)輪機葉片夾角，β為雷達視線與z軸夾角的俯仰角，α為雷達視線在xoy平面上的投影與x軸的夾角，即雷達視線相對于垂直旋轉(zhuǎn)面的方位角。P為風(fēng)輪機葉片上任意散射點，lk表示輪機艙中心到風(fēng)輪機葉片散射點P的距離，R和RP分別表示輪機艙中心和風(fēng)輪機葉片散射點P到氣象雷達的距離，由于(lk/R)2→0，所以RP近似為R-lkcosφ(t)，再去掉載波后，那么由散射點疊加理論可知風(fēng)輪機三張葉片總的回波信號表達式為[13]

(2)

其中λ為雷達波長，設(shè)葉片長度為L，每個葉片有K個散射點，桅桿和輪機艙為靜止目標(biāo)，設(shè)有M個散射點，桅桿高度為H， hm為桅桿上散射點到雷達的距離，同樣利用散射點疊加理論可得回波表達式[13]

(3)

最后整個風(fēng)輪機的回波即為二者之和

sWT(t)=sblade(t)+smast(t)

(4)

對于氣象雷達系統(tǒng)，接收到的氣象雷達回波可以表示為

s=sWT+sweather+n

(5)

其中，sWT表示風(fēng)輪機回波信號， sweather表示氣象目標(biāo)回波信號，n表示高斯白噪聲。

3 傳統(tǒng)匹配追蹤算法

基于傳統(tǒng)匹配追蹤算法抑制氣象雷達風(fēng)電場雜波，就是從一組超完備字典中，找到若干個向量(也稱為原子)盡可能的表示回波中的風(fēng)輪機雜波。把氣象雷達回波中最匹配原子所對應(yīng)的成分從中減去，獲得殘余信號，然后對殘余信號繼續(xù)匹配，重復(fù)迭代，直到雜波殘差值可以忽略為止。對應(yīng)的流程圖如圖2所示。

依據(jù)傳統(tǒng)匹配追蹤算法，字典規(guī)模越大，原子越多，越有可能匹配到最合適原子，算法性能也就越好。然而每次匹配都需要做一次內(nèi)積計算，字典規(guī)模越大，計算量和存儲量自然就越大。當(dāng)風(fēng)輪機三個參數(shù)葉片轉(zhuǎn)速frot、初始夾角θ0、方位角α步長設(shè)置較小時，構(gòu)造字典總的原子數(shù)可達到上百萬個，字典的建立及匹配過程運算量較大，在僅有四臺風(fēng)輪機的情況下就需要幾個小時完成，隨著風(fēng)電場中風(fēng)輪機個數(shù)的增加，運算量將呈現(xiàn)指數(shù)增長。為了解決這個問題，本文利用風(fēng)輪機與雷達位置相關(guān)先驗信息，基于大地主題坐標(biāo)反解算法提取被風(fēng)電場污染的回波數(shù)據(jù)，再利用傳統(tǒng)匹配追蹤方法抑制風(fēng)輪機雜波，進而能夠在保證該方法的良好抑制效果的前提下，能有效降低匹配原子的冗余度，提高運算效率和減少存儲空間。

圖2 傳統(tǒng)匹配追蹤算法流程圖

4 改進的匹配追蹤算法

為了降低傳統(tǒng)匹配追蹤算法計算量，利用風(fēng)輪機距離和方位先驗信息提取氣象雷達回波中被污染的距離單元數(shù)據(jù)。本文首先通過谷歌地球讀取雷達和風(fēng)輪機所在位置的經(jīng)緯度，利用白塞爾大地主題反解方法[15]獲取到每個風(fēng)輪機的精確距離和方位信息，白塞爾方法是解算方法中計算精度較高的一種，能夠滿足對距離和方位信息精度的要求。設(shè)氣象雷達所在位置的經(jīng)緯度分別為λ1、φ1，風(fēng)輪機所在位置經(jīng)緯度分別為λ2、φ2，λ為兩點的經(jīng)差。兩點間的方位角可表示為[15]

(6)

兩點的距離可表示為[15]

tanr=

(7)

由此可知，通過谷歌地球獲取氣象雷達和每個風(fēng)輪機具體的經(jīng)緯度，便可利用白塞爾大地主題反解方法求得其兩者的距離和方位。

利用白塞爾大地主題反解算法獲取氣象雷達與風(fēng)輪機之間精確的位置信息，從而可以提取受風(fēng)電場污染的氣象雷達回波數(shù)據(jù)，即利用距離信息確定風(fēng)輪機所在距離單元，利用方位信息確定氣象雷達掃描到風(fēng)輪機時的相干脈沖回波數(shù)據(jù)。提取字典中有效數(shù)據(jù)，根據(jù)匹配追蹤方法，將接收數(shù)據(jù)與字典原子進行匹配，即投影到不同原子上尋找最大投影值，確定最匹配原子并從回波數(shù)據(jù)中減去其所對應(yīng)分量，計算殘余信號，對殘余信號繼續(xù)匹配，同樣進行反復(fù)迭代，直到風(fēng)輪機雜波殘差值可以忽略為止。算法流程圖如圖3所示。

圖3 改進后的風(fēng)電場雜波抑制流程圖

具體過程如下：

(1)在谷歌地球上，定位每臺風(fēng)輪機以及雷達所在位置，獲取其經(jīng)緯度信息。

(2)根據(jù)白塞爾大地主題坐標(biāo)反解算法，計算每臺風(fēng)輪機相對雷達的距離和方位信息。

(3)根據(jù)距離和方位先驗信息提取雷達回波中污染區(qū)數(shù)據(jù)。

(5)利用匹配追蹤算法抑制風(fēng)電場雜波，設(shè)迭代次數(shù)為n，依據(jù)先驗信息，提取雷達回波中風(fēng)輪機所在距離單元的相干脈沖回波yn，將該距離單元回波投影到字典中的不同原子上，每個原子對應(yīng)的參數(shù)值不同，即

(8)

(9)

(7)求殘余信號，從回波中減去該原子對應(yīng)的分量即為殘余信號

(10)

(8)用殘余信號yn+1代替yn，從(5)重復(fù)該過程進行迭代，直到殘余信號中風(fēng)輪機雜波剩余量值足夠小為止。依次處理M臺風(fēng)輪機所在距離單元的相干脈沖回波，直到所有風(fēng)輪機雜波被抑制。

5 實驗結(jié)果及性能分析

仿真實驗中，氣象目標(biāo)參數(shù)、氣象雷達參數(shù)及風(fēng)輪機參數(shù)如表1～表3所示。

表1 氣象目標(biāo)參數(shù)

表2 氣象雷達參數(shù)

表3 風(fēng)輪機仿真參數(shù)

實際風(fēng)電場具有幾十甚至上百臺風(fēng)輪機，本實驗考慮四臺風(fēng)輪機情況，相關(guān)參數(shù)如表4所示。

表4 四個風(fēng)輪機參數(shù)

5.1 實驗結(jié)果

為了模擬真實天氣回波的功率譜，信噪比(SNR)一定的情況下，在氣象目標(biāo)中加入平均功率為1的加性高斯白噪聲。如圖4所示，只含有噪聲的氣象目標(biāo)回波功率譜圖和時頻圖可以看出，氣象目標(biāo)的頻譜分布在一定的頻譜范圍內(nèi)，具有一定的頻譜擴展，其中心頻率和譜寬由氣象目標(biāo)相對于氣象雷達的平均徑向速度以及速度譜寬所決定。

整個風(fēng)電場的氣象雷達回波如圖5所示，圖5(a)是氣象雷達掃描范圍的距離方位圖，四臺風(fēng)輪機分布在相對氣象雷達的不同距離和方位上，圖中四個較亮的點就是四臺風(fēng)輪機所對應(yīng)的位置，由圖可以看出主要分布在方位30°～60°和距離20～30km之間。

圖4 含有噪聲的氣象目標(biāo)回波

圖5 氣象雷達回波的距離方位和距離多普勒譜

圖6 不同初始夾角的氣象雷達回波功率譜和時頻譜

圖5(b)是氣象雷達掃描范圍的距離多普勒，可以看出氣象目標(biāo)頻譜具有隨距離連續(xù)分布特性，另外四條明顯的頻譜是風(fēng)輪機所對應(yīng)的頻譜。提取風(fēng)電場中第四臺風(fēng)輪機所在距離單元的回波，其功率譜和時頻譜如圖6所示。葉片的初始夾角θ0不同，風(fēng)輪機的頻譜分布情況就不同。圖6(a)當(dāng)初始夾角為0°時，可以看出在零頻處和正負頻率處風(fēng)輪機能量較強，零頻處是參考葉片的回波能量，而另外兩個葉片關(guān)于參考葉片對稱，所以其頻率也是對稱分布。圖6(b)初始夾角為30°時，由于雷達視線與參考葉片垂直，葉片上所有的點到氣象雷達距離相同，回波較強，另外兩個葉片回波能量相對較弱，所以能量較強的部分主要在零頻到最大多普勒頻率分布。在功率譜和時頻譜中可以看出兩種信號的不同特征，本文根據(jù)風(fēng)輪機雜波不同于氣象目標(biāo)特征，采用匹配追蹤算法抑制風(fēng)輪機雜波并保留真實的氣象目標(biāo)。

在谷歌地球上獲得四臺風(fēng)輪機的經(jīng)緯度如表5所示。利用該先驗信息獲得每臺風(fēng)輪機相對于雷達的距離和方位，然后提取有效相干脈沖對應(yīng)的回波數(shù)據(jù)，基于匹配追蹤方法抑制風(fēng)電場雜波。

表5 四個風(fēng)輪機和雷達的經(jīng)緯度

采用改進的匹配追蹤方法實現(xiàn)風(fēng)電場雜波抑制的結(jié)果如圖7所示，圖7與圖5對比分析可以看出，圖5(a)中四臺風(fēng)輪機所在距離方位單元的回波已經(jīng)剔除，同時圖5(b)中四臺風(fēng)輪機對應(yīng)的四條展寬的頻譜也基本剔除。圖8是第四臺風(fēng)輪機所在距離單元的雜波抑制后的功率譜及時頻譜，與真實的氣象目標(biāo)結(jié)果(如圖4所示)一致，由圖8和圖4對比可以看出，風(fēng)電場雜波得到了有效的抑制。

圖7 風(fēng)電場雜波抑制結(jié)果

圖8 風(fēng)輪機雜波抑制結(jié)果

5.2 性能分析

(1)算法改進前后性能分析

在參數(shù)設(shè)置、回波數(shù)據(jù)以及字典原子的數(shù)量相同的條件下，改進前后的算法運算時間對比如表6所示(該實驗是在一個中央處理器為英特E5-2640的服務(wù)器上完成)。由表6數(shù)據(jù)對比可以看出，改進前抑制處理時間約三個小時，改進后五分半鐘左右?？梢娕c傳統(tǒng)方法相比，改進算法在保證抑制性能不變的前提下，明顯提高了運算效率。

表6 改進前后抑制時間對比

(2)抑制前后風(fēng)輪機雜波平均功率

為了分析改進后的匹配追蹤算法性能，本文在不同雜噪比(CNR)的條件下，進行了500次的蒙特卡羅實驗，對比抑制前后風(fēng)輪機雜波平均功率變化。抑制前后的風(fēng)輪機雜波平均功率隨著雜噪比(CNR)的變化性能曲線如圖9所示。

圖9 抑制前后風(fēng)輪機雜波平均功率隨雜噪比的變化

圖9(a)和圖9(b)中的兩條曲線分別是抑制前后風(fēng)輪機雜波平均功率的變化情況，隨雜噪比的增大，抑制前的雜波平均功率逐漸變大，抑制后雜波功率逐漸變小并趨于0dB以下的一個穩(wěn)定值。除此之外，圖9(a)和圖9(b)對比可以看出，在不同信噪比的條件下該方法都能保持良好的抑制效果，信噪比的變化并不影響該方法的性能。需要說明的是，當(dāng)雜噪比小于信噪比(如圖9(a)中CNR小于SNR=10dB)時，抑制效果不理想，這是因為此時雜波平均功率也就等于甚至低于氣象目標(biāo)平均功率，依據(jù)匹配追蹤算法原理找最大投影值，可能會出現(xiàn)匹配錯誤的原子，導(dǎo)致抑制效果不理想。也就是在信雜比(SCR)較高時，該方法抑制效果不明顯。事實上，這種情況并不需要抑制雜波，因為在信雜比(SCR)較高時，也就是氣象目標(biāo)回波較強，此時風(fēng)輪機對氣象目標(biāo)影響較小，無需進行雜波抑制。例如，當(dāng)信噪比是10dB時，雜噪比(CNR)分別為10dB和0dB兩種情況下的風(fēng)輪機雜波抑制前后結(jié)果如圖10所示。圖10(a)和圖10(c)分別是雜噪比為10dB和0dB時抑制前的雷達回波，由于雜噪比較低甚至低于信噪比，氣象目標(biāo)特性占主導(dǎo)，不能明顯地看出風(fēng)輪機的頻譜分布。圖10(b)和圖10(d)是其相對應(yīng)的抑制結(jié)果，雜波功率較小，抑制前后結(jié)果變化不明顯。

(3)抑制前后氣象目標(biāo)參數(shù)估計結(jié)果

基于脈沖對法對風(fēng)電場雜波抑制前后數(shù)據(jù)分別進行氣象目標(biāo)參數(shù)(平均功率，平均徑向速以及平均譜寬)估計結(jié)果如表7所示。

圖10 抑制前后的雷達回波的功率譜和時頻譜對比

目標(biāo)狀態(tài)平均功率/dB平均徑向速度/(m/s)平均譜寬/(m/s)真實氣象目標(biāo)10.0002.02.0抑制前氣象雷達回波29.925911.08148.6550抑制后氣象雷達回波10.0011.98671.9913

由表7可以看出，抑制前的氣象雷達回波參數(shù)估計值出現(xiàn)較大的偏差，由此可知風(fēng)電場雜波對氣象目標(biāo)參數(shù)估計的影響較大。對氣象雷達回波抑制后，參數(shù)估計值更接近于真值，雜波抑制算法效果較為顯著。

6 結(jié)論

本文利用白塞爾大地主題反解算法獲得雷達和風(fēng)輪機的距離、方位信息，提取出雷達回波中被風(fēng)電場污染的相干脈沖數(shù)據(jù)，再對該數(shù)據(jù)作雜波抑制處理。實驗結(jié)果顯示，改進方法在較好的保留了原有的氣象目標(biāo)信息的前提下抑制了大部分的風(fēng)電場雜波，同時有效地降低了計算量和數(shù)據(jù)存儲量，提高了運算效率。除此之外，從雜波功率去除率及氣象信息保留完好性兩個方面對算法的性能進行了定量的分析驗證算法的有效性。本文算法中，只考慮了風(fēng)輪機直達路徑回波情況，沒有考慮風(fēng)電場中風(fēng)輪機間的多徑散射對抑制效果的影響，這也是后續(xù)需要進一步研究的工作。