基于解卷積的自適應(yīng)空間譜估計方法

邱 嵐，謝于晨

(江西科技學(xué)院信息工程學(xué)院，江西 南昌 330098)

0 引言

自適應(yīng)空間譜估計方法可按最大信干噪比(minimum signal to interference plus noise ratio，SINR)、最小二乘(least squares，LS)、最小方差信號無畸變響應(yīng)(minimum variance distortion response，MVDR)等準則，根據(jù)噪聲(含干擾)自適應(yīng)調(diào)整導(dǎo)向權(quán)向量，實現(xiàn)期望功率最大輸出，所得空間譜具有較高的方位分辨率和檢測能力，被廣泛應(yīng)用于聲納、雷達、通信等領(lǐng)域[1-3]。但這類方法需要傳感器陣列拾取數(shù)據(jù)具有一定信噪比，當拾取數(shù)據(jù)中信噪比含有量較低時，該方法所得空間譜的背景級將顯著升高，甚至不能準確估計目標空間方位角分布[4-5]。為了降低自適應(yīng)空間譜估計中背景級，研究學(xué)者提出了Chebyshev濾波[6]、凸優(yōu)化(convex optimization)[7]、分子陣預(yù)處理[8]等方法。但還是存在如下問題：Chebyshev濾波方法存在背景級設(shè)置和主瓣寬度控制折中選擇問題；凸優(yōu)化方法針對空間譜背景級、主瓣寬度等多個性能指標進行優(yōu)化，但這些參數(shù)需要綜合考慮陣列模型結(jié)構(gòu)，如果設(shè)置不當有可能達不到預(yù)期優(yōu)化結(jié)果；分子陣預(yù)處理方法提升背景級有限，且在傳感器陣列數(shù)量有限情況下，還存在主瓣變寬問題。針對MVDR空間譜估計方法存在的旁瓣較高且抑制干擾性能不穩(wěn)健的情況，文獻[9]提出一種旁瓣級可控的自適應(yīng)波束形成后置處理方法。該方法在MVDR 基礎(chǔ)上進行峰值搜索，將獲得的峰值點從大到小進行排序，取次大值作為最高旁瓣的值，然后在其方位添加虛擬干擾加以抑制，從而得到新的波束圖，實現(xiàn)對旁瓣級控制，該方法與凹槽噪聲場法類似，同樣存在旁瓣抑制過程增加了主瓣寬度問題。針對自適應(yīng)空間譜估計方法背景級對目標空間方位角分布值估計的影響問題，本文提出了基于解卷積的自適應(yīng)空間譜估計方法，稱之為DMVDR方法。

1 MVDR方法

1.1 信號模型

假設(shè)K個相互獨立的遠場目標所輻射的頻率f信號(之間互不相關(guān))被一個N元傳感器組成的均勻線列陣拾取，如圖 1所示。

圖1 均勻線列陣示意圖Fig.1 The schematic diagram of uniform linear array

令K個遠場目標相對該線列陣所在空間方位角為ΘK=[θ1,θ2,…,θK]，則圖1所示的均勻線列陣各傳感器拾取數(shù)據(jù)可表示為：

X(f)=A(ΘK)S(f)+V(f)

(1)

式(1)中，S(f)=[S1(f),S2(f),…,SK(f)]T為K個目標輻射頻率f信號，V(f)=[V1(f),V2(f),…,VK(f)]T為線列陣拾取背景噪聲在頻率f上分量，A(ΘK)為K個遠場目標相對該線列陣陣列流形矩陣，(·)T為矩陣轉(zhuǎn)置運算符。

為了后續(xù)論述的方便性，令u=sinθ，即u∈[-1,1]為方向正弦，則陣列流形矩陣A(ΘK)可表示：

(2)

式(2)中，c為聲速。

1.2 MVDR方法空間譜估計

令線列陣各傳感器拾取的背景噪聲屬于高斯白噪聲，之間相互獨立，K個遠場目標輻射信號與背景噪聲之間互不相關(guān)，線列陣拾取數(shù)據(jù)的協(xié)方差矩陣RX(f)可表示為：

RX(f)=X(f)HX(f)=
A(ΘK)HRS(f)A(ΘK)+RV(f)

(3)

式(3)中，RS(f)和RV(f)分別為目標信號和背景噪聲協(xié)方差矩陣，(·)H為矩陣共軛轉(zhuǎn)置運算符。

此時，由RX(f)可獲得權(quán)向量最優(yōu)解為：

(4)

式(4)中，a(u)=[1,ej2πfdu/c,…,ej2πf(N-1)du/c]T。

根據(jù)獲得的權(quán)向量最優(yōu)解，可得搜索方位角θ對應(yīng)的空間譜估計值為：

(5)

對協(xié)方差矩陣RX(f)進行特征分解[10]，可得：

(6)

式(6)中，λk和λm分別為信號和噪聲子空間對應(yīng)特征值，Ek和Em分別為信號和噪聲子空間對應(yīng)特征向量。

此時，式(5)可進一步表示為：

(7)

當搜索方位角θ屬于空間方位角ΘK，由于噪 聲子空間Em與a(uk)正交性，式(7)可進一步表示為：

(8)

由式(8)可知，當搜索方位角θ屬于空間方位角ΘK，MVDR方法空間譜估計過程中自動截取了線列陣拾取數(shù)據(jù)的協(xié)方差矩陣，保留了K個遠場目標信號子空間信息。由于信號子空間可更好地作用于目標信號檢測和方位角估計，理想情況下，MVDR方法空間譜估計結(jié)果近似為狄利克函數(shù)。

2 DMVDR方法

2.1 點擴展函數(shù)

自適應(yīng)空間譜估計的最終目的是對K個遠場目標方位角ΘK進行估計，假設(shè)相對線列陣，K個遠場目標在空間方位角分布為Ω，則Ω中的一個元素Ω(u)可表示為：

Ω(u)=Akδ(u-uk)

(9)

式(9)中，Ak為第k個目標信號幅度，δ為狄利克函數(shù)。

通過空間譜實現(xiàn)目標空間方位角ΘK估計等同于通過空間譜實現(xiàn)Ω估計，即MVDR方法輸出空間譜PMVDP(f,u)可比作是對目標空間方位角分布Ω的估計。為了后續(xù)分析方便，本文借鑒圖像復(fù)原處理中點擴展函數(shù)描述空間譜PMVDR(f,u)和空間方位角分布Ω之間關(guān)系，即MVDR方法輸出空間譜PMVDR(f,u)可表示為點擴展函數(shù)h(u)與空間方位角分布Ω的卷積。

PMVDR(f,u)=Ω(u)*h(u)+v(u)

(10)

式(10)中，*為一維卷積運算符，v(u)為空間譜形成過程所含噪聲分量。

由式(10)可知，在已知點擴展函數(shù)情況下，可通過解卷積處理，實現(xiàn)對目標空間方位角分布Ω估計，進而獲得目標方位ΘK估計值。

另外，由式(7)可知，當搜索方位角θ屬于空間方位角ΘK，式(7)變?yōu)槭?8)，空間譜取得最大值。而當搜索方位角θ不屬于空間方位角ΘK時，式(7)中分母由兩部分組成，空間譜取得最小值。因此，可利用一個“類狄利克函數(shù)”構(gòu)造MVDR方法輸出空間譜PMVDR(f,u)的點擴展函數(shù)，即：

(11)

式(11)中，ξ∈(0,1)的常數(shù)。由式(11)可知，MVDR空間譜的點擴展函數(shù)對于入射信號方向正弦具有平移不變性。

將式(7)中MVDR空間譜表示成卷積形式為：

RMVDR(f,u)=∏(u)*h(u)+v(u)

(12)

式(12)中，∏(u)為k目標信號所在方位a(uk)到信號子空間Ek倒數(shù)的真實分布，即：

(13)

2.2 解卷積技術(shù)

圖像復(fù)原處理[11-14]常用解卷積技術(shù)主要有Wiener濾波、clean方法、Richardson-Lucy方法[15-16]。Wiener濾波是在高斯分布模型基礎(chǔ)上，使均方差最小，隸屬于線性方法；clean方法是在多個脈沖形函數(shù)加權(quán)組合模型基礎(chǔ)上對輸入信號進行組合，對信號要求嚴格；Richardson-Lucy方法基于貝葉斯理論，采用最小鑒別準則在非負實函數(shù)空間使鑒別函數(shù)最小，隸屬于非線性方法。本文解卷積任務(wù)是在已知含噪聲模糊圖像和其點擴展函數(shù)來估計清晰圖像，使其滿足輸出模糊圖像和對應(yīng)的清晰圖像距離最小，采用Richardson-Lucy方法進行解卷積迭代計算，經(jīng)過一定次數(shù)的迭代后可得到清晰圖像的最大似然估計[17]。

由于篇幅限制，本文只給出Richardson-Lucy方法的迭代公式，具體原理和推導(dǎo)可參考文獻[18—19]。

(14)

式(14)中，i為迭代次數(shù)，h(x,y)為已知點擴展函數(shù)，g(x,y)為已知模糊圖像，fi(x,y)為第i次迭代計算所得清晰圖像估計值。在迭代過程中可采用同一常數(shù)作為迭代初始值獲得解卷積結(jié)果，但無法獲得各目標在空間譜中的貢獻相對大小，即強弱目標相應(yīng)幅度差信息將會丟失，為了保留強弱目標相應(yīng)幅度信息，本文采用f0(x,y)=g(x,y)對迭代計算進行初始化[17]。

由式(14)可見，Richardson-Lucy方法只需在給定點擴展函數(shù)的條件下，通過迭代方法便可從輸入模糊圖像中估算原始清晰圖像。由于該方法沒有迭代終止條件，隨著迭代次數(shù)的增加，fi(x,y)將逐漸收斂于原始清晰圖像。

2.3 DMVDR方法空間譜估計

本文依據(jù)圖像處理中利用點擴展函數(shù)和模糊圖像解卷積恢復(fù)清晰圖像的原理，針對自適應(yīng)空間譜估計中背景級較高問題，提出基于解卷積的自適應(yīng)空間譜估計方法。該方法可作為現(xiàn)有自適應(yīng)空間譜估計方法的一種后置處理方法，用于實際應(yīng)用中，即利用MVDR方法輸出空間譜歸一化平均背景級作為點擴展函數(shù)非目標方位處值實現(xiàn)點擴展函數(shù)設(shè)計，然后基于該點擴展函數(shù)對MVDR方法所得空間譜PMVDR(f,u)進行解卷積，逐步恢復(fù)K個遠場目標在空間方位角分布圖，最終估計出目標空間方位角。

利用Richardson-Lucy方法對式(12)進行迭代解卷積，可估計出∏(u)。

(15)

采用式(15)進行迭代估計時，迭代次數(shù)與式(11)中常數(shù)ξ有關(guān)，且ξ取值大小將影響解卷積空間譜估計的性能。本文中，ξ取為歸一化空間譜平均值，相當于歸一化平均背景級，即：

(16)

式(16)中，mean(·)和max(·)分別為取均值和最大值函數(shù)，在采用mean(·)函數(shù)取均值時應(yīng)扣除最大值周圍的峰高值。

設(shè)式(15)中最終的迭代次數(shù)為I次，則經(jīng)過第I次迭代后的∏(u)即為解卷積所得空間譜，即：

PDMVDR(f,u)=∏(u)I

(17)

依據(jù)上面所述數(shù)據(jù)處理過程，可將本文方法實現(xiàn)過程分為如下步驟：

步驟1 對線列陣拾取數(shù)據(jù)進行濾波、時頻變換，采用時頻變換結(jié)果的單頻帶數(shù)據(jù)進行協(xié)方差矩陣估計，得到協(xié)方差矩陣估計結(jié)果RX(f)；

步驟2 利用MVDR方法進行空間譜估計，獲得含有噪聲的空間譜PMVDR(f,u)；

步驟3 利用空間譜PMVDR(f,u)，根據(jù)式(11)和式(16)構(gòu)造點擴展函數(shù)h(u)；

步驟4 將PMVDR(f,u)設(shè)置為初始化值，并設(shè)定迭代次數(shù)，根據(jù)式(15)進行Richardson-Lucy迭代解卷積，得到目標空間方位角分布估計值∏(u)即為DMVDR方法空間譜PDMVDR(f,u)；

步驟5 根據(jù)解卷積所得DMVDR方法空間譜，按式(18)進行累積求和，得到寬帶空間譜：

(18)

式(18)中，flow和fhigh分別為空間譜估計處理頻帶下限和上限。

2.4 DMVDR方法運算量分析說明

DMVDR方法作為一種后置處理方法，是在MVDR方法基礎(chǔ)上進行解卷積實現(xiàn)的，相比MVDR方法，該方法增加了2Θ2+Θ乘法和Θ2加法運算，Θ為空間譜估計所含方位角個數(shù)。下例直觀給出了DMVDR方法解卷積增加時間，采用 Inter Corei7 2核處理器(Windows10 操作系統(tǒng))，在Matlab2018a編程環(huán)境下，對MVDR輸出空間譜(Θ=180)進行100次迭代解卷積所需運算時間為0.013 s。由此可見，DMVDR方法可滿足實際應(yīng)用中的實時性要求，且實際應(yīng)用中也可將解卷積部分單獨作為并行處理的一部分，進一步降低DMVDR方法串行運算時間。

3 數(shù)據(jù)處理分析

3.1 數(shù)值仿真分析

為了驗證本文方法可以很好地降低自適應(yīng)空間譜估計中的背景級。接下來對DMVDR方法輸出空間譜進行如下數(shù)值仿真分析。數(shù)值仿真中采用圖 1形式的32元均勻線列陣作為數(shù)據(jù)拾取平臺，相鄰傳感器間距為2 m，采樣率為5 kHz，一次拾取數(shù)據(jù)長度為1 s。仿真中采用MVDR方法，MVDR分子陣方法(文獻[8]方法)，MVDR旁瓣級迭代處理方法(文獻[9]方法)和DMVDR方法(本文方法)對不同目標數(shù)情況進行對比分析。

1) 目標數(shù)K=1情況

目標相對線列陣方位角為0°(正橫方向)，目標輻射信號為375 Hz單頻連續(xù)信號，其與噪聲譜級比為SNR，噪聲為加性高斯白噪聲。圖 2—圖4分別給出了SNR為10、0、-10 dB情況下，由4種方法所得空間譜。其中，DMVDR方法給出了不同迭代次數(shù)I情況下所得空間譜，文獻[9]方法設(shè)置期望旁瓣級為-40 dB。

圖2 SNR=10 dB時4種方法所得空間譜Fig.2 Spatial spectrum of three methods vs SNR=10 dB

圖3 SNR=0 dB時4種方法所得空間譜Fig.3 Spatial spectrum of three methods vs SNR=0 dB

圖4 SNR=-10 dB時4種方法所得空間譜Fig.4 Spatial spectrum of three methods vs SNR=-10 dB

由圖 2—圖4顯示結(jié)果可得到如下結(jié)論：

當線列陣拾取數(shù)據(jù)所含信號譜級比SNR為10、0和-10 dB時，MVDR方法所得空間譜背景級約為-25、-15和-5 dB；文獻[8]方法通過對線列陣拾取數(shù)據(jù)進行分子陣處理，使其輸出信噪比空間譜背景級達到-35、-20和-10 dB，降低了空間譜背景噪聲級；文獻[9]方法通過在最高旁瓣的方位添加虛擬干擾源和50次以上的迭代以實現(xiàn)對旁瓣級的控制，降低了空間譜背景噪聲級，達到了設(shè)置期望旁瓣級，但主瓣寬度被加寬；而DMVDR方法分別通過一次和四次迭代解卷積實現(xiàn)了與文獻[8]、文獻[9]一樣的效果，且隨著解卷積迭代次數(shù)的增多，所得空間譜背景級得到大幅度的降低，譜峰也變得更加“尖銳”；在解卷積迭代次數(shù)為4、6和8時，DMVDR方法所得空間譜的背景級達到-80 dB以下，遠低于MVDR方法、文獻[8]方法、文獻[9]方法，能更好地估計目標空間方位角分布值。

由上可知，在進行空間譜估計時，DMVDR方法可通過增加解卷積迭代次數(shù)使其輸出空間譜背景級遠低于MVDR方法、文獻[8]方法、文獻[9]方法；在處理數(shù)據(jù)所含信號譜級比較低時，DMVDR方法可保持較好的背景級，便于對目標空間方位角分布值實現(xiàn)估計。

2) 目標數(shù)K=2情況

為了進一步驗證本文方法可以降低背景級對目標方位角估計的影響。接下來采用兩個目標進行仿真分析，兩個目標相對線列陣方位角分別為0°(正橫方向)和-3°，目標輻射信號均為375 Hz單頻連續(xù)信號，背景噪聲為加性高斯白噪聲，目標1與背景噪聲譜級比為-12 dB，目標2與背景噪聲譜級比為-15 dB。圖 5給出了由4種方法所得空間譜，其中，DMVDR方法給出了不同迭代次數(shù)I情況下所得空間譜，文獻[9]方法設(shè)置期望旁瓣級為-40 dB。

圖5 4種方法所得空間譜Fig.5 Spatial spectrum of three methods

由目標數(shù)K=1仿真所得結(jié)論及圖 5所結(jié)果可知：當線列陣拾取數(shù)據(jù)所含信號譜級比SNR=0 dB時，MVDR方法所得空間譜背景級約為-15 dB，致使在目標信號譜級比SNR=-15 dB時，受背景級影響，目標方位角分布值估計效果較差；文獻[8]方法通過對線列陣拾取數(shù)據(jù)進行分子陣處理，降低了空間譜背景噪聲級，輸出空間譜能夠顯示出目標方位分布，但是相對背景級還是較高；文獻[9]方法通過在最高旁瓣的方位添加虛擬干擾源和100次的迭代以實現(xiàn)對旁瓣級的控制，降低了空間譜背景噪聲級，達到了設(shè)置期望旁瓣級，但主瓣寬度被加寬，無法實現(xiàn)對兩目標辨識；而DMVDR方法在解卷積迭代次數(shù)為4次時，其所得空間譜的背景級達到-40 dB，已遠低于MVDR方法和文獻[8]方法，能較好地顯示目標空間方位角分布，且兩目標之間背景級更低、空間譜峰更加“尖銳”，具有更好的目標空間方位角高分辨估計能力。該結(jié)果進一步說明了DMVDR方法作為一種后置處理方法，繼承了MVDR方法，在線列陣拾取數(shù)據(jù)所含信號譜級較低時，能夠通過多次迭代解卷積使其輸出空間譜保持較低的背景級，能夠更好地估計出目標空間方位角分布值。

3.2 實測數(shù)據(jù)處理分析

為了進一步驗證本文方法對目標檢測和方位估計性能，接下來利用某次試驗數(shù)據(jù)對上述4種方法進行處理分析。

試驗中，接收陣為圖 1所示結(jié)構(gòu)的32元線列陣，線列陣采樣率為5 kHz，介質(zhì)平均聲速為1 496 m/s，背景噪聲級為90 dB@100 Hz，F(xiàn)IR濾波器頻帶為200～400 Hz，該段數(shù)據(jù)在-45°附近存在一個合作目標，其他為非合作目標。一次處理數(shù)據(jù)長度為10 000個采樣數(shù)據(jù)，MVDR方法分200個子帶進行空間譜估計，單個子帶協(xié)方差矩陣估計中由40個快拍實現(xiàn)，每個快拍包含256個數(shù)據(jù)樣本，相鄰快拍之間重疊128個數(shù)據(jù)樣本。圖 6—圖9給出了由3種方法所得空間譜方位歷程圖，3種方法所得方位歷程圖已進行歸一化顯示，并將動態(tài)范圍約束在-50～0 dB之間(如圖中右側(cè)色棒所示)，以便于對不同方法所得空間譜背景級壓制能力進行對比分析。其中，DMVDR方法給出了不同迭代次數(shù)I=4情況下所得空間譜方位歷程圖，文獻[9]方法設(shè)置期望旁瓣級為-40 dB。圖 10給出了t=1 s 時刻，4種方法所得空間譜。

圖6 MVDR方法輸出空間譜方位歷程圖Fig.6 Spatial spectrum BTR of MVDR method

圖7 文獻[8]方法輸出空間譜方位歷程圖Fig.7 Spatial spectrum BTR of Ref.8 method

圖8 文獻[9]方法輸出空間譜方位歷程圖Fig.8 Spatial spectrum BTR of Ref.9 method

圖9 DMVDR方法輸出空間譜方位歷程圖(I=4)Fig.9 Spatial spectrum BTR of DMVDR method(I=4)

由圖6—圖9給出的空間譜方位歷程圖可知，在該時段內(nèi)，由于試驗時存在多個目標，線列陣拾取數(shù)據(jù)所含目標信號較多、目標信號相對背景噪聲信噪比較復(fù)雜，MVDR方法和文獻[8]方法空間譜估計性能發(fā)生了退化，輸出空間譜背景級被抬高，損失了動態(tài)范圍(小于20 dB)；而DMVDR方法通過多次解卷積迭代處理，有效降低了空間譜背景級，可獲得良好的動態(tài)范圍(優(yōu)于50 dB)，具有很強的背景級壓制能力。由圖 6和圖 10可知，在合作目標空間方位角附近，MVDR方法輸出空間譜背景級約為-13 dB，合作目標所在空間方位空間譜峰約為-13 dB，較難估計出合作目標空間方位角；由圖 7和圖 10可知，文獻[8]方法通過對線列陣拾取數(shù)據(jù)進行分子陣處理，使其輸出信噪比空間譜背景級達到-20 dB，降低了空間譜背景噪聲級，能夠估計出合作目標空間方位角，但方位歷程圖顯示效果不夠清晰；文獻[9]方法通過在最高旁瓣的方位添加虛擬干擾源和70次的迭代以實現(xiàn)對旁瓣級的控制，使其輸出信噪比空間譜背景級達到-40 dB，降低了空間譜背景噪聲級，能夠估計出合作目標空間方位角，但主瓣被加寬；由圖 9和圖 10可知，DMVDR方法通過4次迭代解卷積處理后輸出空間譜背景級低于-50 dB，能夠穩(wěn)定地估計出合作目標空間方位角，且方位歷程圖顯示效果清晰可辨。

圖10 4種方法所得空間譜(t=1 s)Fig.10 Spatial spectrum of three methods (t=1 s)

實測數(shù)據(jù)處理結(jié)果進一步驗證了，在線列陣拾取數(shù)據(jù)所含目標信號信噪比較低時，MVDR方法輸出空間譜背景級較高，較難估計出目標空間方位角分布。文獻[8]和文獻[9]的方法雖然使MVDR方法輸出信噪比空間譜背景級得到降低，能夠估計出合作目標空間方位角，但均存在一定待改善的地方，而利用DMVDR方法所得空間譜方位歷程圖顯示效果清晰可辨，能夠更好地估計出合作目標空間方位角，可作為現(xiàn)有自適應(yīng)空間譜估計方法的一種后置處理方法用于實際應(yīng)用中。

4 結(jié)論

本文提出了一種基于解卷積的自適應(yīng)空間譜估計方法(DMVDR)。該方法首先在自適應(yīng)空間譜估計方法輸出空間譜基礎(chǔ)上，采用類狄利克函數(shù)和空間譜歸一化平均背景級實現(xiàn)點擴展函數(shù)設(shè)計，然后利用點擴展函數(shù)和圖像復(fù)原中Richardson-Lucy迭代方法對自適應(yīng)空間譜估計結(jié)果進行解卷積，降低背景級及其對目標空間方位角分布值估計的影響。數(shù)值仿真和實測數(shù)據(jù)處理結(jié)果表明，該方法作為一種后置處理方法，繼承了自適應(yīng)空間譜估計方法高分辨估計能力，空間譜峰值更加“尖銳”；信號譜級比較低時，通過增加迭代次數(shù)使其輸出空間譜背景級遠低于自適應(yīng)空間譜方法和分子陣處理方法，具有更好的目標空間方位角分布估計能力。