孫大軍， 馬超， 梅繼丹， 楊宛珊， 魏秋雨

(1.哈爾濱工程大學 水聲技術重點實驗室， 黑龍江 哈爾濱 150001；2.海洋信息獲取與安全工信部重點實驗室(哈爾濱工程大學) 工業(yè)和信息化部，黑龍江 哈爾濱 150001；3. 哈爾濱工程大學 水聲工程學院， 黑龍江 哈爾濱 150001)

陣列信號處理是現代信號處理領域的重要分支，被廣泛的應用在雷達、聲吶中，高分辨波束形成技術作為提高聲吶目標探測性能的重要手段，幾十年來一直是國內外學者的研究熱點。但是由于聲吶陣列規(guī)模的逐漸增大和本艇的機動等特殊性的應用環(huán)境，使得理想的陣列流型和長時間信號平穩(wěn)性難以保證，導致了許多高分辨陣列信號處理方法的實際應用效果并不理想。近年來具有高穩(wěn)健性和實時性(小快拍、短樣本)特征的高分辨陣列處理方法已成為了陣列信號處理中新的研究熱點。

由于經典的常規(guī)波束形成算法受到空域傅里葉變換極限分辨能力的限制，如何突破瑞利限曾一度成為廣大學者的重要研究方向。自Capon波束形成技術被提出以來，產生了一批高分辨空間譜估計方法，比如諧波分析法[1]、最大熵法[2]、最小方差無畸變響應法(minimum variance distortionless response, MVDR)[3]、基于ARMA (autoregressive moving average model)的線性預測法，MUSIC (multiple signal classification)[4]算法，ESPRIT (estimating signal parameter via rotational invariance techniques)[5]算法等。由于MVDR和MUSIC以及ESPRIT這類經典的高分辨處理方法無法分辨相干源，之后又出現了一系列的解相干處理方法。空間平滑技術其基本思想是利用均勻線陣的平移不變性，將陣列劃分為多個重疊子陣，通過對所有子陣協(xié)方差矩陣的平均來解相干。但在水聲信號處理領域，由于經典類高分辨處理算法對陣列流形精度要求高、低信噪比性能下降嚴重等問題，實用性一直受到質疑。學者們又針對具有良好干擾抑制效果的MVDR算法開展了穩(wěn)健性算法研究，典型的方法有線性約束最小方差算法[6]、對角加載方法[7]，基于特征空間波束形成等方法[8]。但是根據經典的陣列信號處理理論，穩(wěn)健性、高分辨與處理增益往往不可兼得，高分辨方法穩(wěn)健性差，而為了提高穩(wěn)健性對權值進行對角加載這類方法往往又會導致實際處理增益下降，現有方法大多是在這三者之間尋找平衡點，所以每種方法都會有應用局限性。

到了20世紀，人們開始利用其他學科的成果來解決陣處理中的高分辨問題，凸優(yōu)化理論也被應用到陣列信號處理中，Eldar等[9]運用凸優(yōu)化理論，將波束形成問題轉化為二階錐優(yōu)化問題，突破了傳統(tǒng)波束形成方法中存在的未能求取全局最優(yōu)解的局限。隨后，壓縮感知理論被推廣于陣列信號處理領域[10]，相比于前幾種高分辨算法，該方法的優(yōu)點是可應用于小快拍數據，無需預先估計聲源個數，可應用于陣元數少于聲源數的情形等，但受噪聲功率估計誤差影響較大，且計算量大，使工程應用受限。

Yang[11-12]將反卷積算法與波束形成相結合，提出一種高增益穩(wěn)健反卷積波束形成技術，并應用到聲壓陣和圓陣的處理中，對反卷積波束形成的穩(wěn)健性、分辨率和增益方面進行了驗證，改善效果顯著。但文獻中所提出的基于Richardson-Lucy(RL)的反卷積波束形成方法只適用于PSF具有移不變特性的陣列中。針對上述問題，文獻[13]針對矢量陣PSF的移變性，將非負最小二乘(nonnegative least-squares, NNLS)和DAMAS (deconvolution approach for the mapping of acoustic sources)算法應用到矢量陣中，并提出了一種改進的R-L算法[13-14]，該方法在低信噪比的情況下與NNLS和DAMAS 相比，具有更窄的主瓣寬度，更低的旁瓣和更高的穩(wěn)健性。多目標情況下的常規(guī)波束形成空間譜輸出可以看成每一個角度的指向性圖和該角度的源強度乘積之和，在數學上可以用一個疊加積分來表示。如果陣列的PSF具有移不變性，此時疊加積分可以被表示為卷積過程。因此，對常規(guī)波束輸出空間譜和自然指向性函數反卷積可以得到目標函數。

在實際陣列信號處理中，由于應用場景和需求的不同，陣列的形式多種多樣，除了等間距直線陣、圓陣、矢量陣等典型陣列形式外，不等間距陣、十字陣、面陣等陣列形式也是常采用的陣列形式。本文主要分析了水聲常用陣列結構PSF的特點，將多種典型陣列按照卷積模型類型進行了分類綜述，將其劃分為波束圖移變陣列和波束圖移不變陣列。然后分別總結了移不變和移變模型陣列的反卷積波束形成典型求解方法，以期為反卷積技術在水聲陣列信號處理中的應用提供一定的指導。

1 移不變點擴散函數的反卷積波束形成問題

1.1 移不變的點擴散函數定義

在空間域處理中，有一大類陣列的自然指向性函數直接或者通過簡單的變換就能夠具有PSF移不變的特點。在這里移不變的定義就是指向某個角度的波束輸出是和角度無關的，其可以看做自然指向性函數的移位，即R(θ|?)=R(θ-?)?；蛘咴谀骋粋€變換域具有波束輸出具有移不變性，不同位置的波束輸出可以視指向性函數的移位。這種移不變性可能是在一維，也可能是在二維或高維空間。

1.2 典型移不變模型陣列

1.2.1 圓陣

對于某些陣列，PSF是移不變條件是可以直接滿足的，例如，一個等間距圓陣，如果僅考慮二維平面的情況下其自然指向性函數為:

(1)

式中：N為陣元數；r為圓陣半徑；λ為波長。則其在角度?的指向性函數可表示為：

(2)

當陣元數N滿足N>4πr/λ+2時，指向性函數可以近似為：

(3)

其中，J0(x)為零階貝塞爾函數：

(4)

此時，圓陣?方向的指向性函數R(θ|?)是θ-?的函數，即R(θ|?)=R(θ-?)，近似具有關于方位?的移不變性質。

常規(guī)波束形成功率譜P(θ)可以看成每一個角度的指向性函數R(θ|?)和該角度的信源強度S(?)乘積之和，在數學上可以用疊加積分表示：

(5)

由于圓陣的指向性函數具有移不變特性，波束輸出疊加積分的過程可以等效為卷積過程：

(6)

窄帶信號頻率3 kHz，圓陣的陣元數20，半徑0.5 m，改變信號的方向，源1的方位為0°，源2為30°，源3為60°。指向不同角度的圓陣指向性函數如圖1所示，可以看出圓陣的指向性函數是自然指向性函數的圓周移位，因此是移不變的。

圖1 指向不同角度的圓陣的指向性函數Fig.1 The PSF of a circular array steering to different angles

.2.2 均勻聲壓線列陣

另外還有一類陣列的指向性函數是需要進行一定的映射變換才能夠滿足移不變性質的，例如均勻等間距聲壓線陣、非均勻線陣。以均勻聲壓線陣為例，均勻線陣幾何模型圖如圖2所示。其自然指向性函數為：

圖2 均勻線陣模型Fig.2 Array geometry for uniform line array

(7)

其在角度?方向的指向性函數為：

(8)

這里取平方是為了和CBF的功率譜對應，顯然R(θ|?)≠R(θ-?)。因此，疊加過程不能直接用卷積過程代替。但是經過簡單的變換，單位線陣的自然指向性函數具有關于cosθ移不變的性質，即有R(cosθ|cos ?)=R(cosθ-cos ?)。

因此均勻線陣波束輸出可以等效為卷積過程：

R(cosθ)*S(cosθ)

(9)

窄帶信號，信號頻率500 Hz，直線陣的陣元數20，陣元間距半波間距，信號來波方位的余弦值分別為0、0.25和0.5。指向不同角度的聲壓線陣的指向性函數如圖3所示，可以看出指向不同角度的指向性函數和角度無關，具有移不變性，可以看做自然指向性函數的移位。

圖3 指向不同角度的聲壓線陣的指向性函數Fig.3 The PSF of ULA steering to different angles

1.2.3 非均勻聲壓線列陣

假設N個聲壓傳感器非均勻的分布在同一條直線上，信號從陣列的遠場入射，如圖4所示。

圖4 非均勻線陣模型Fig.4 Array geometry for non-uniformly linear array

假設入射信號是單頻信號，第1個陣元H1接收到的信號是Acos(2πft)，把1號陣元當做參考陣元，那么第i個傳感器接收到的信號為：

si(t)=Acos(2πft+φi)，i=1,2,…，N

(10)

其中φi為聲程差HiPi導致的相位差:

(11)

將陣元接收到的信號相加，整個陣列的輸出為：

(12)

將s(t)取平方并歸一化得到的指向性函數為：

(13)

相似地，其在角度?方向的指向性函數為：

(14)

顯然，非均勻線陣的指向性函數是sinθ-sin ?的函數，因此非均勻線陣的PSF是關于角度的正弦值的移不變函數，R(θ|?)=R(sinθ-sin ?)。非均勻線陣常規(guī)波束形成功率可以寫做如下卷積過程:

R(sinθ)*S(sinθ)

(15)

為了說明非均勻線陣移不變性，以20元非均勻陣為例，各陣元的位置分別選取半波間距陣列的[1,2,3,4,5,6,8,10,12,14,16,18,20,22,24,26,27,30,31,35]號陣元構成非等間距陣列。信號是500 Hz單頻信號，其中λ對應信號頻率的波長。

圖5 指向不同角度的非均勻聲壓線陣的指向性函數Fig.5 The PSF of NLA steering to different angles

1.2.3 遠場平面陣

一個任意分布的N元平面陣，如圖6所示放置于xoy平面，假設信號從(θ0,φ0)方位以遠場平面波入射，N個陣元位置坐標分別為[xn,yn]，n=1,2,…,N。則第n個陣元與參考陣元之間的聲程差為：

圖6 任意面陣模型Fig.6 Array geometry for planar array

βn=k(xnsinφcosθ+ynsinφsinθ)

(16)

其中，k為波數，k=ω/c=2π/λ，λ為信號的波長，定義ux=sinφcosθ，uy=sinφsinθ，那么基陣的陣列流形向量為a(k)=(e(-jk(uxx1+uyy1)),…,e(-jk(uxxN+uyyN)))T，常規(guī)波束形成器的加權向量為wc=a(θ0,φ0)/M。因此常規(guī)波束形成的指向性函數為：

Rp(ux,ux0,uy,uy0)=|wHa(θ,φ)|2=

Rp(ux-ux0,uy-uy0)

(17)

故面陣的波束可以看做是一個關于ux=sinφ·cosθ和uy=sinφsinθ的二維移不變函數，且ux和uy都在[-1,1]范圍內，因此可以將矩形面陣歸為PSF移不變陣的二維反卷積求解問題。以矩形陣為例說明平面陣的二維反卷積問題，矩形陣仿真條件：X軸陣元數8，Y軸陣元數8，信號形式為頻率500 Hz單頻信號，陣元間距為1.5 m。

圖7是矩形面陣在不同角度所對應的單目標指向性函數。由圖可見，2個指向性函數形狀相同，只是進行了平移，因此面陣的指向性函數可以看成是移不變的。

圖7 矩形陣指向不同角度的三維指向性函數Fig.7 PSF of rectangle array steering to different source angles

.3 PSF移不變模型反卷積求解方法

對于移不變的點擴散函數，常用的卷積方法有CLEAN[15]、DAMAS[16]、NNLS[17]、R-L[18-19]算法等。在其基礎上又衍生出適用于相干信號的CLEAN-C算法，基于傅里葉變換的快速DAMAS2[20]和DAMAS3[20]算法，基于稀疏信號重構的SC-DAMAS[21]算法，基于小波變換計算網格的快速DAMAS算法[22]，基于傅里葉變換的FFT-NNLS算法等。Klaus在文獻[17]中比較了DAMAS，NNLS,基于傅里葉變化的DAMAS2, FFT-NNLS和R-L算法，結果表明DAMAS和NNLS不要求PSF是移不變的，而DAMAS2,FFT-NNLS和RL算法適用于PSF移不變的反卷積問題中，且RL算法和DAMAS,NNLS相比具有更低的計算量。文獻[23]將CLEAN算法，DAMAS,NNLS和R-L算法應用到三維空氣聲源定位中，結果表明幾種方法都能顯著提高聲源分辨能力，但在處理相干源時R-L的效果最好。文獻[24]分析比較了DAMAS,NNLS和ISCA 3種反卷積算法，結果表明DAMAS的定位誤差精度最差， ISCA有最高的定位精度和最高的穩(wěn)健性。在這里簡要介紹幾種常用的反卷積算法。

1.3.1 DAMAS

將系統(tǒng)的輸入輸出描述成線性方程組的形式：

P=RS+N

(18)

式中：P代表常規(guī)波束形成的輸出；R代表基陣的指向性函數PSF；S代表源的分布函數；N代表噪聲。DAMAS是在忽略噪聲的條件下用高斯賽德爾方法求解線性方程組。對于線性移不變的PSF，矩陣R是Toeplitz 矩陣且是一個半正定對稱陣。

高斯賽德爾方法通過引入正約束來解決式(18)，迭代過程：

(19)

(20)

1.3.2 NNLS算法

最小二乘問題，其基本原理可以推廣到陣列反卷積處理，具體方法是在常規(guī)波束輸出，矢量陣點擴散函數字典，聲源目標函數之間建立差函數方程組，通過最小化差函數的原則來實現對目標函數S的求解，獲取聲源的分布，矩陣為：

φ=‖RS-P‖2

(21)

式中：P是矢量陣常規(guī)波束形成輸出空間譜；R是所有角度矢量陣的指向性函數形成的指向性函數字典矩陣；S是反映聲源方位和強度信息的目標函數。對于上述方程組求解可以采用梯度投影法，該方法的核心是負梯度方向指向標量場下降最快的方向，通過在φ關于S的負梯度方向上按特定步長反復迭代搜索來求得S的結果。計算的具體步驟可以參考文獻[13]。

NNLS算法并不局限于移不變反卷積模型求解，只需要獲得陣列所有方位的指向性函數矩陣即可。需要注意的是，對于移不變模型采用NNLS算法時可采用FFT變換域計算的方式實現快速計算，但這并不適用于移變模型陣列。該方法并不局限于矢量陣。對結構固定不變的陣列該方法都適用。

1.3.3 R-L算法

R-L算法是一種基于貝葉斯理論的迭代方法，該方法被廣泛用于圖像復原中。R-L算法中要求PSF是移不變的，在陣列信號處理中，就是指基陣的波束圖是和角度無關的(指向某一個角度的波束圖都可以看做自然指向性函數的移位)，R-L算法的迭代方程為：

R(sinθ)*Si(sinθ)

(22)

(23)

dsinθ≤ε

(24)

1.3.4 仿真對比

以均勻分布聲壓線陣為例，分析陣列的PSF是移不變的反卷積問題求解。由式(9)可知，ULA的波束輸出為：

P(cosθ)=R(cosθ)*S(cosθ)

(25)

基于陣列的參數可以得到陣列的自然指向性函數，對接收到的數據做常規(guī)波束形成可以得到波束輸出的方位譜，對常規(guī)波束輸出空間譜P(θ)和自然指向性函數R(θ)反卷積可以得到目標函數S(θ)，其理想情況下是δ函數，必然是對目標方位的高分辨估計結果。關于PSF是移不變的反卷積問題求解引用文獻[14]中所采用的R-L方法。

為了說明反卷積波束形成算法在PSF是移不變的陣列的應用效果，以圓陣為例進行仿真分析。

仿真條件同圖1，為理想無噪聲條件下，雙目標的仿真結果圖。圖8說明常規(guī)波束輸出的空間譜是陣列的PSF和信號的目標強度函數的卷積。與CBF相比，反卷積波束形成方法有更窄的主瓣和更低的旁瓣。

圖8 雙目標圓陣仿真結果Fig.8 Simulation results of two target for circular array

2 移變點擴散函數的反卷積波束形成問題

2.1 移變的點擴散函數

在空間域處理中，有一大類陣列的自然指向性函數具有PSF移變的特點，移變的定義就是指向某個角度的波束輸出是和角度有關的，即R(θ|?)≠R(θ-?)，不可以看做自然指向性函數的簡單移位。但與一般的時域系統(tǒng)相比，空間域上的移變問題有一個典型特點：在水聲中一般是將具有固定拓撲結構指的陣列去采集聲信號，對于任意形狀的陣列，當陣列拓撲結構固定后，PSF無論是移不變還是移變都是可以預先測定得到的，如果采用測量獲得的PSF函數由于其包含了對陣型誤差等因素等影響的校正會有更好的反卷積求解效果。

2.1.1 矢量陣

以均勻矢量線陣為例來說明這一問題，角度定義同圖1。與聲壓陣相似，直線聲矢量陣的波束輸出可以寫為[16]：

(26)

式中：Rv(θ|?)矢量陣指向?的指向性函數；S(?)是以方位為變量的目標信源強度分布函數；θ和?分別代表掃描角度和目標信號方位。

矢量陣在90°方向指向性函數為：

(27)

式中：Rp(θ)為聲壓陣的自然指向性函數，可以看出，矢量陣的自然指向性函數Rv(θ)為聲壓陣自然指向性函數Rp(θ)與一個角度指向因子u(θ)=(1+cosθ)/2的乘積。

指向?方向的矢量陣指向性函數為：

Rv(θ|?)=

(28)

比較Rv(θ)和Rv(θ|?)，不難發(fā)現Rv(θ|?)≠Rv(θ-?)，也不可通過簡單的正弦，余弦變換將其變?yōu)橐撇蛔兒瘮怠Ｒ驗?，u(θ)=(1+cosθ)/2是關于θ移不變的函數，而Rp(θ)是關于cosθ的移不變函數，二者相乘后無法通過簡單的變量代換將其等效為某個統(tǒng)一變量的移不變模型。

矢量陣仿真：20元均勻分布矢量陣，陣元間距3.75 m，為信號頻率半波間距，信號頻率200 Hz，采樣頻率20 kHz。改變信號的方向，指向不同角度的矢量陣的指向性函數如圖9所示，可以看出不同角度的矢量陣的指向性函數是移變的，不是自然指向性函數的移位。而且盡管矢量有一定的抗左右模糊作用，但是其效果與目標方位有關，在端射方向附近效果仍然較差，會在對稱方向形成偽峰。

圖9 指向不同角度的矢量陣指向性函數Fig.9 PSF of VSLA steering to different source angles

2.1.2 共形陣

共形陣是傳感器附著在非平面載體表面上的陣列，簡稱共形陣。由于其具有良好的空氣動力學特性，可以充分利用導流罩的縱向空間，增加基陣有效孔徑，共形陣正被當做一種新型布陣方式應用在新型聲吶系統(tǒng)中。

假設各陣元的位置pn是隨機分布，如圖10所示，第nth個傳感器的位置坐標為pn=[xn,yn,zn]。K個平面波從陣列的遠場入射，簡單起見，假設K=1。θ和φ分別是方位角和俯仰角，定義信號傳播方向的單位向量可以表示為：u=-[cosθcosφ,sinθcosφ,sinφ]T，把原點當做參考點時，第nth個傳感器相對于參考點的時延差為：

圖10 任意分布陣列模型Fig.10 Arbitrary conformal array geometry

(29)

式中c為聲速。

若入射信號為Acos(2πft)，那么第nth個傳感器接收到的信號為：

sn(t)=Acos[2πf(t+τn)]

(30)

將接收到的信號求和，并將幅度歸一化并且取平方，得到陣列的指向性函數為：

相似的，其在角度?方向的指向性函數為：

R(θ,φ)=

(32)

顯然，共形陣的指向性函數不具有規(guī)律可言，是移變函數。為了對該問題簡要說明，這里采用半橢圓陣列來說明共形陣反卷積波束形成問題。對于共形陣聲吶(非成像聲吶)，探測目標為陣所在平面內的遠距離目標，即只考慮目標方位方向上的一維求解問題，此時是一維移變反卷積問題。

將41個傳感器如圖11所示，放置在橢圓軌跡上，橢圓方程是 4y2+x2=100; 每個傳感器在x軸的位置是xn,其中xn=[0:0.5:10]∪[9.5:-0.5:0]。

圖11 橢圓形共形陣模型Fig.11 Oval conformal array geometry

信號頻率為1 kHz，改變信號的方向，指向不同角度的共形陣的指向性函數如圖12所示，可以看出不同角度的共形陣的指向性函數是移變的，不是自然指向性函數的移位。

圖12 指向不同角度的共形陣指向性函數Fig.12 PSF of conformal array steering to different source angles

2.2 PSF是移變模型的反卷積求解方法

文獻[14]針對矢量陣分析比較了3種反卷積波束形成方法，DAMAS，NNLS和改進的R-L算法。改進的R-L算法在低信噪比的情況下與NNLS和DAMAS 相比，具有更窄的主瓣寬度，更低的旁瓣和更高的穩(wěn)健性。因此針對空間域反卷積問題點擴散函數移變情況的特點——移變但可預知的特點，應用文獻[14]中提出的改進R-L反卷積算法。

(33)

上述方法在空間域處理中并不限于之前提到的

幾種陣列，凡是指向性函數移變但可預測的情況都適用，而對于具有固定拓撲結構的實際陣列PSF都是可以預先測定得到的，因此適用的陣列形式非常廣泛。該方法亦可推廣到非空間域陣列信號處理領域的其他單位沖擊響應移變但可預測的系統(tǒng)中。

為了證明反卷積波束形成方法在PSF是移變陣列波束形成中的使用性能，以橢圓形共形陣為例進行仿真分析，橢圓形共形陣的參數見2.1.2節(jié)。圖13展現了預存的不同角度指向性函數(PSF)。

圖13 預存的不同角度共形陣指向性函數Fig.13 The predefined beam pattern

圖14是對單目標及雙目標的共形陣做CBF與dCv的仿真對比結果，陣列結構同圖12，圖14(a)是單目標，信號方位90°，圖14(b)是等強度的雙目標，信號方位是100°和240°。

圖14 共形陣CBF和dCv的波束輸出對比Fig.14 Comparison of dCv and CBF in conformal array

通過上圖可以看出，新方法可以很好地將反卷積波束形成算法應用到具有任意結構的共形陣這種PSF具有移變性的陣列中。與CBF相比，dCv方法有更高的角度分辨能力和低的旁瓣，同時能準確估計出多目標信號的功率值。

通過上述具有橢圓形狀的共形陣的仿真結果也證明了文中所提出的PSF移變模型的反卷積波束形成方法是具有可行性的。

3 結論

1)文章介紹了水下陣列信號處理中的反卷積問題。對多種水下常用陣列形式的波束形成卷積模型進行了分析推導，將其按照卷積模型類型進行了分類綜述，劃分為波束圖移變陣列和波束圖移不變陣列。

2)圓陣、均勻聲壓線列陣、非均勻聲壓線列陣均是一維移不變陣列，平面陣是二維移不變陣列，矢量陣和共形陣是移變陣列。

本文總結了移不變和移變模型陣列的反卷積波束形成典型求解方法，是對各類典型水下陣列卷積模型的基礎總結歸納。