基于頻域低秩矩陣補(bǔ)全的氣象雷達(dá)WTC抑制算法

邱存銀， 沈明威， 韓國棟

(1. 河海大學(xué)， 江蘇 南京 211106;2. 中國電子科技集團(tuán)公司第五十四研究所， 河北 石家莊 050081)

0 引 言

隨著能源不斷的被開發(fā)利用，煤炭、石油等不可再生能源正在不斷枯竭，世界各國對可再生能源的需求也更加強(qiáng)烈。近年來，諸多國家對風(fēng)電等的研究逐步深入，技術(shù)也愈發(fā)成熟，這使得風(fēng)電場的數(shù)量和規(guī)模呈指數(shù)倍增加。然而研究表明，風(fēng)力渦輪機(jī)大型葉片高速旋轉(zhuǎn)會造成雷達(dá)回波的多普勒頻譜持續(xù)展寬，氣象目標(biāo)被淹沒于風(fēng)力渦輪機(jī)雜波(Wind Turbine Clutter，WTC)的多普勒頻譜中，嚴(yán)重影響了雷達(dá)系統(tǒng)對氣象目標(biāo)檢測的精確性[1-3]。因此，國內(nèi)外學(xué)者詳細(xì)分析了氣象雷達(dá)不同工作模式下風(fēng)力渦輪機(jī)雜波WTC與氣象回波的時(shí)、頻域分布特性，發(fā)現(xiàn)由于WTC的時(shí)變性，使得針對靜止雜波的時(shí)域?yàn)V波、頻域?yàn)V波等傳統(tǒng)雜波抑制技術(shù)均失效，此時(shí)傳統(tǒng)的雜波抑制技術(shù)無法在抑制WTC的同時(shí)無損恢復(fù)氣象信息[4-9]。

針對上述問題，本文根據(jù)氣象信號的低秩分布特性，將稀疏優(yōu)化理論引入到WTC抑制中，研究了基于矩陣補(bǔ)全(Matrix Completion，MC)的氣象雷達(dá)WTC雜波抑制新方法。MC算法是一種源于壓縮感知的秩最小化理論，其主要作用是在存在數(shù)據(jù)缺失或只能觀測到部分?jǐn)?shù)據(jù)的低秩觀測矩陣中，利用矩陣元素之間的相關(guān)性，有效地恢復(fù)矩陣的缺失數(shù)據(jù)[10]。本文的目的是在氣象信號矩陣受到WTC污染后，利用MC對污染后缺失的數(shù)據(jù)進(jìn)行填充恢復(fù)，在抑制WTC的同時(shí)還能夠精確恢復(fù)氣象信號矩陣。另外將常用的頻域多重二次插值與頻域MC算法進(jìn)行對比，理論分析和實(shí)驗(yàn)結(jié)果表明，頻域MC算法較頻域多重二次插值能更好地恢復(fù)受WTC干擾的氣象信號。

1 雷達(dá)回波數(shù)據(jù)預(yù)處理

在雷達(dá)系統(tǒng)中，假設(shè)距離單元同時(shí)包含WTC和天氣信號。在本文實(shí)驗(yàn)中給定被污染的方位-距離單元在相干積累時(shí)間(Coherent Integration Time， CIT)內(nèi)接收到的信號表示為

xl(m)=sl(m)+cl(m)+wl(m)+nl(m)

(1)

式中：m=1,2，…,M，M是脈沖個數(shù)；l為某一距離單元；sl(m)為氣象信號；cl(m)為地雜波信號；wl(m)為WTC；nl(m)為其他噪聲信號。

由于雜波集中在零頻附近，對氣象信號不造成干擾，因此設(shè)計(jì)脈沖可以將其濾除。此時(shí)，被污染的雷達(dá)方位-距離單元在CIT內(nèi)接收到的信號可表示為

xl(m)=sl(m)+wl(m)+nl(m)

(2)

信噪比(Signal to Noise Ratio，SNR)是影響氣象信號恢復(fù)精度的決定因素。因此，可將實(shí)測數(shù)據(jù)的雷達(dá)回波信號從時(shí)域轉(zhuǎn)換至多普勒頻域，通過相干積累提高SNR，這一過程可通過快速傅里葉變換(Fast Fourier Transform，F(xiàn)FT)實(shí)現(xiàn)。多普勒頻域的氣象信號X(k)為

k=1,2,…,K,l=1,2,…,L

(3)

式中：K為地雜波對消后的脈沖數(shù)；L為距離單元。在頻域進(jìn)行矩陣元素補(bǔ)全，以實(shí)現(xiàn)低SNR情況下，被污染距離單元?dú)庀笮盘柕母呔然謴?fù)。

2 頻域矩陣補(bǔ)全算法

2.1 低秩矩陣重構(gòu)

文獻(xiàn)[11-12]證明了當(dāng)矩陣的奇異值具有稀疏性(即矩陣是低秩的)且采樣數(shù)目滿足一定條件時(shí)，大多數(shù)矩陣可通過求解核范數(shù)最小化問題來精確地恢復(fù)所有元素。所以對于任何抽樣的數(shù)據(jù)，如果將一行或一列作為一個整體進(jìn)行抽樣，那么此時(shí)觀測矩陣不能滿足元素采樣集均勻隨機(jī)采樣的條件，缺失的元素也就不能被高精度恢復(fù)。因此，需選取特殊的矩陣形式重新構(gòu)造觀測矩陣，使信號的采樣點(diǎn)隨機(jī)分布在整個回波數(shù)據(jù)中，滿足低秩重構(gòu)的條件。

對于預(yù)處理后的多普勒域信號X(k)，將其重構(gòu)成k1×k2維的Hankel矩陣，其中k1+k2-1=L。用雜波檢測算法提取出WTC所在距離單元l，并將干擾元素置零，在重構(gòu)的矩陣中，這些零元素將被重新排列以進(jìn)行隨機(jī)采樣[13]。根據(jù)Hankel矩陣模型和式(1)，取氣象信號距離-多普勒頻域矩陣的第k列，構(gòu)造的Hankel矩陣S如下

構(gòu)造一個方陣B=S×SH，根據(jù)矩陣的特征分解理論，可以得到

(4)

式中：Σ=diag(σ1,σ2,…,σr)，σi為矩陣S的全部奇異值；V為特征向量排成列的矩陣。根據(jù)r(SSH)=r(S)，因此，矩陣B秩的數(shù)值可由該矩陣的奇異值數(shù)量推導(dǎo)得出。假設(shè)k=470，則矩陣B的特征值分布如圖1所示。

圖1 仿真數(shù)據(jù)構(gòu)建的Hankel矩陣特征值分布

圖1為氣象信號仿真數(shù)據(jù)下構(gòu)建的Hankel矩陣特征值分布，可以明顯看出，重構(gòu)Hankel矩陣的秩r=2?51，此時(shí)觀測矩陣滿足所需低秩的要求[14]。

2.2 基于NNM-IALM算法的WTC抑制

矩陣補(bǔ)全是一種重建隨機(jī)抽樣低秩矩陣或近似低秩約束的過程。根據(jù)MC理論，當(dāng)雷達(dá)回波信號被WTC污染時(shí)，可去除所有被污染的樣本，只保留剩余的氣象數(shù)據(jù)，即可將對WTC抑制轉(zhuǎn)化為一個矩陣補(bǔ)全問題。

根據(jù)上述Hankel矩陣的構(gòu)造方法，將預(yù)處理后的頻域信號第k個脈沖下的列向量構(gòu)造為Hankel矩陣X。同理，用S，C，W，N分別表示構(gòu)造的關(guān)于Sl(k)，Cl(k)，Wl(k)，Nl(k)的Hankel矩陣。

為了解決上述矩陣補(bǔ)全問題，構(gòu)造一個矩陣補(bǔ)全模型，通過約束秩最小化的問題，根據(jù)已知的矩陣元素補(bǔ)全未知元素。去除WTC后，氣象數(shù)據(jù)可以通過下式中核范數(shù)最小化(Nuclear Norm Minimization，NNM)模型進(jìn)行恢復(fù)

(5)

其中，

(6)

式中：min(·)表示最小化；‖·‖*為矩陣的核范數(shù)；PΩ為正交投影算子，表示投影到僅在指標(biāo)集Ω非零的稀疏矩陣子空間上的映射。

利用文獻(xiàn)[15]中非精確增廣拉格朗日乘子(Inexact Augmented Lagrange Multiplier，IALM)求解式(4)中的優(yōu)化問題時(shí)，將其擴(kuò)展到復(fù)數(shù)域，其對應(yīng)的拉格朗日函數(shù)可表示為

L(S,N,Y,μ)=‖S‖*+〈Y,X-S-N〉+

‖S‖*+tr(YH((X-N)-S))+

(7)

定義一個收縮運(yùn)算符為

(8)

利用IALM交替更新S，N，Y，μ，其迭代步驟如下。

步驟1：固定變量N，Y和μ，更新變量S

(9)

對X-Nk+μk-1Y進(jìn)行奇異值分解

(U,Σ,VH)=svd(X-Nk+μk-1Yk)

(10)

則S的迭代更新公式為

Sk+1=UT(2μk)-1[Σ]VH=

(11)

步驟2：固定變量S，Y和μ，更新變量N

(12)

對式(12)中的N進(jìn)行求導(dǎo)可得：

(13)

步驟3：固定變量S，N和μ，更新變量Y

Yk+1=Yk+μ(X-Sk+1-Nk+1)

(14)

步驟4：更新變量μ

μk+1=ρμk

(15)

其中，Y表示拉格朗日乘子矩陣；μ>0表示懲罰因子；ρ為觀測數(shù)據(jù)率；=tr(XHY)表示矩陣的內(nèi)積；tr(·)表示矩陣的跡。

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

本節(jié)通過計(jì)算機(jī)驗(yàn)證算法的性能，本次實(shí)驗(yàn)使用的實(shí)測數(shù)據(jù)來自于荷蘭代爾夫特理工大學(xué)的IDRA雷達(dá)[16]。實(shí)測數(shù)據(jù)參數(shù)如表1所示。

表1 Hankel矩陣仿真參數(shù)

選取氣象信號實(shí)測數(shù)據(jù)的第30個距離單元作為本次實(shí)驗(yàn)被污染的距離單元，由圖2(a)可知，氣象信號分布在425 Hz～526 Hz，在第30個距離單元中加入WTC，從圖2(b)可以看出，氣象信號淹沒在WTC中。圖2為被WTC污染前后的多普勒圖。

圖2 加入WTC前后氣象信號多普勒圖

從圖3可以看出，實(shí)測數(shù)據(jù)構(gòu)建的Hankel矩陣特征值分布下降的趨勢相較于圖1仿真數(shù)據(jù)構(gòu)建的Hankel矩陣更緩慢，該現(xiàn)象主要是由實(shí)測數(shù)據(jù)中噪聲成分的復(fù)雜性造成，但近似滿足低秩特性，因此可采用MC方法高精度恢復(fù)氣象信號[17]。

圖3 實(shí)測數(shù)據(jù)構(gòu)建的Hankel矩陣特征值分布

圖4為分別使用頻域多重二次插值、時(shí)域NNM-IALM和頻域NNM-IALM算法恢復(fù)的氣象信號多普勒圖。從圖中可以看出，三種算法都可以有效抑制WTC并恢復(fù)氣象信號。但是在氣象信號恢復(fù)的精度上，頻域NNM-IALM算法恢復(fù)的效果明顯優(yōu)于頻域多重二次插值算法和時(shí)域NNM-IALM算法，這是由于在矩陣補(bǔ)全實(shí)驗(yàn)前對實(shí)測數(shù)據(jù)的雷達(dá)回波信號進(jìn)行預(yù)處理，將雷達(dá)回波信號從時(shí)域轉(zhuǎn)換至多普勒頻域，通過相干積累提高SNR，進(jìn)而提高實(shí)驗(yàn)的恢復(fù)精度。另外在時(shí)域的條件下對氣象信號補(bǔ)全，不僅不能抑制噪聲而且氣象信號的旁瓣較高，對氣象信號的恢復(fù)產(chǎn)生了影響，這一問題在多普勒頻域?qū)崿F(xiàn)補(bǔ)全時(shí)得到了解決。在實(shí)際工程應(yīng)用中，并不需要將所有頻點(diǎn)處的信號恢復(fù)，我們只需恢復(fù)氣象信號所在的頻率范圍即可，這樣可極大的減少運(yùn)算量。

圖4 不同算法恢復(fù)出的氣象信號功率譜

使用FFT法對氣象信號的徑向速度與譜寬進(jìn)行估計(jì)，將兩種算法的譜矩結(jié)果與真實(shí)數(shù)據(jù)做對比，采用交叉驗(yàn)證法來評價(jià)兩種算法的性能。

根據(jù)譜矩估計(jì)量的定義，信號功率表示為

(16)

平均頻率表示為

(17)

譜寬表示為

(18)

式中：S(fi)為信號功率譜fi=i/KTr，i=1,2，…,K,Tr為脈沖重復(fù)間隔；K為脈沖個數(shù)。

氣象回波的徑向速度和速度譜寬分別為

(19)

(20)

由實(shí)驗(yàn)估計(jì)出被WTC污染的第30個距離單元?dú)庀笮盘柕膹较蛩俣群妥V寬如表2所示。

表2 仿真數(shù)據(jù)性能對比 m/s

由表2的實(shí)驗(yàn)數(shù)據(jù)可以看出，在相同的實(shí)驗(yàn)環(huán)境下，頻域NNM-IALM算法比頻域多重二次插值能更精確的恢復(fù)出氣象信號，恢復(fù)后的徑向速度和速度譜寬也更接近原始的實(shí)測數(shù)據(jù)且估計(jì)精度比頻域多重二次插值算法分別提高了約22.25%和31.48%。綜上，頻域NNM-IALM算法能夠更好的抑制WTC，并有效恢復(fù)被污染的氣象信號。

4 結(jié) 語

本文將頻域矩陣補(bǔ)全運(yùn)用到氣象雷達(dá)實(shí)測數(shù)據(jù)的WTC抑制中。把預(yù)處理后的多普勒域回波信號重構(gòu)為等效的低秩矩陣，并通過IALM算法來有效地恢復(fù)缺失數(shù)據(jù)，同時(shí)在相同的實(shí)驗(yàn)條件下與頻域多重二次插值進(jìn)行對比。實(shí)驗(yàn)結(jié)果表明，上述兩種方法都達(dá)到了雜波抑制的效果，實(shí)現(xiàn)了氣象信號的恢復(fù)。但基于Hankel矩陣構(gòu)建的高精度頻域補(bǔ)全算法無論在恢復(fù)性能還是降噪效果上都更優(yōu)，并解決了低信噪比下時(shí)域矩陣補(bǔ)全的高旁瓣問題。