基于解卷積的自適應空間譜估計方法

邱 嵐，謝于晨

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

0 引言

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

1 MVDR方法

1.1 信號模型

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

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

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

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為矩陣轉置運算符。

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

(2)

式(2)中，c為聲速。

1.2 MVDR方法空間譜估計

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

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

(3)

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

此時，由RX(f)可獲得權向量最優解為：

(4)

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

根據獲得的權向量最優解，可得搜索方位角θ對應的空間譜估計值為：

(5)

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

(6)

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

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

(7)

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

(8)

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

2 DMVDR方法

2.1 點擴展函數

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

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

(9)

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

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

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

(10)

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

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

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

(11)

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

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

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

(12)

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

(13)

2.2 解卷積技術

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

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

(14)

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

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

2.3 DMVDR方法空間譜估計

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

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

(15)

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

(16)

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

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

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

(17)

依據上面所述數據處理過程，可將本文方法實現過程分為如下步驟：

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

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

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

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

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

(18)

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

2.4 DMVDR方法運算量分析說明

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

3 數據處理分析

3.1 數值仿真分析

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

1) 目標數K=1情況

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

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

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

2) 目標數K=2情況

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

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

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

3.2 實測數據處理分析

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

試驗中，接收陣為圖 1所示結構的32元線列陣，線列陣采樣率為5 kHz，介質平均聲速為1 496 m/s，背景噪聲級為90 dB@100 Hz，FIR濾波器頻帶為200～400 Hz，該段數據在-45°附近存在一個合作目標，其他為非合作目標。一次處理數據長度為10 000個采樣數據，MVDR方法分200個子帶進行空間譜估計，單個子帶協方差矩陣估計中由40個快拍實現，每個快拍包含256個數據樣本，相鄰快拍之間重疊128個數據樣本。圖 6—圖9給出了由3種方法所得空間譜方位歷程圖，3種方法所得方位歷程圖已進行歸一化顯示，并將動態范圍約束在-50～0 dB之間(如圖中右側色棒所示)，以便于對不同方法所得空間譜背景級壓制能力進行對比分析。其中，DMVDR方法給出了不同迭代次數I=4情況下所得空間譜方位歷程圖，文獻[9]方法設置期望旁瓣級為-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給出的空間譜方位歷程圖可知，在該時段內，由于試驗時存在多個目標，線列陣拾取數據所含目標信號較多、目標信號相對背景噪聲信噪比較復雜，MVDR方法和文獻[8]方法空間譜估計性能發生了退化，輸出空間譜背景級被抬高，損失了動態范圍(小于20 dB)；而DMVDR方法通過多次解卷積迭代處理，有效降低了空間譜背景級，可獲得良好的動態范圍(優于50 dB)，具有很強的背景級壓制能力。由圖 6和圖 10可知，在合作目標空間方位角附近，MVDR方法輸出空間譜背景級約為-13 dB，合作目標所在空間方位空間譜峰約為-13 dB，較難估計出合作目標空間方位角；由圖 7和圖 10可知，文獻[8]方法通過對線列陣拾取數據進行分子陣處理，使其輸出信噪比空間譜背景級達到-20 dB，降低了空間譜背景噪聲級，能夠估計出合作目標空間方位角，但方位歷程圖顯示效果不夠清晰；文獻[9]方法通過在最高旁瓣的方位添加虛擬干擾源和70次的迭代以實現對旁瓣級的控制，使其輸出信噪比空間譜背景級達到-40 dB，降低了空間譜背景噪聲級，能夠估計出合作目標空間方位角，但主瓣被加寬；由圖 9和圖 10可知，DMVDR方法通過4次迭代解卷積處理后輸出空間譜背景級低于-50 dB，能夠穩定地估計出合作目標空間方位角，且方位歷程圖顯示效果清晰可辨。

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

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

4 結論

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