淺海低頻長線陣聲源方位估計?

梁玉權 周士弘 宮在曉

(1 中國科學院聲學研究所 聲場聲信息國家重點實驗室 北京 100190)

(2 中國科學院大學 北京 100049)

0 引言

根據淺海波導低頻聲傳播的簡正波理論和平面波假設，水平陣接收點源激發的聲場是由不同的模態疊加而成，每階模態對應一個來波方向，每個來波方向是由聲源實際方位角和聲源激發模態的俯仰角耦合而成。如果采用常規波束形成(Conventional beamforming,CBF)的聲源方位估計方法，將會產生估計方位相對真實方位的偏移和分裂現象，特別是聲源靠近長線陣端射方向時，這種現象更加明顯。解決該問題的基本思路是將聲傳播模型引入到方位估計中[1]，為此，很多學者提出了基于水平陣和淺海簡正模態理論的聲源方位估計方法。

1998年，Yang[2]將水平直線陣的陣元域-場匹配處理，擴展為波束域-波束匹配處理，實現與匹配場處理類似的淺海聲源方位的無偏估計。然而，采用全模態場的陣元域或波束域匹配處理，并不能改善固有的環境失配的影響，且進行深度、距離和方位角三維參數空間搜索運算量較大。2004年，Lakshmipathi等[3]提出了子空間相交(Subspace intersection,SI)方法，將淺海目標聲源方位估計與簡正波模型結合起來，利用信號子空間與模態子空間的相交特征，對方位角進行一維的空間優化搜索，縮減了計算量。2006年，張愛民等[4]提出約束最小二乘子空間相交方法和總體最小二乘子空間相交方法，改善SI方法在進行矩陣奇異值分解時的數值穩定性。2013年，易峰等[5]基于正交投影子空間分解原理進行淺海聲源方位估計。由于SI方法估計信號子空間準確性依賴于樣本協方差矩陣奇異值分解，因此并不適用于快拍數較少的低信噪比接收信號，在存在相干源和未知聲源個數的情況下，子空間估計精度受限。

2006年以來，壓縮感知理論得到了各個研究領域的關注，并在聲學方面得到廣泛的應用和發展[6]。其中，塊稀疏壓縮感知是一個重要的發展方向，其恢復算法可大體分為3類：基于凸優化的塊稀疏恢復算法[7]、基于塊稀疏貝葉斯學習算法[8]和貪婪匹配追蹤類算法[9]。

本文基于簡正波理論，將淺海水平陣接收聲場模型表示為聲源方位空間的塊稀疏信號模型，將聲源方位估計問題表示為塊稀疏恢復問題，基于陣元域的接收聲場信號，提出一種塊正交匹配追蹤(Block orthogonal matching pursuit,BOMP)的淺海聲源的方位估計方法，該方法屬于貪婪匹配追蹤類算法的一種，區別于SI方法，所提方法無需構建協方差矩陣以及特征分解估計子空間，適用于單快拍接收聲場信號。BOMP方法是3類塊稀疏恢復算法中原理最簡單、收斂速度最快的算法[10]。通過仿真對比分析CBF和BOMP淺海聲源方位估計的結果，來說明所提方法是否能準確估計靠近陣列端射方向的聲源方位，并討論信噪比的影響，最后通過海上實驗數據驗證所提方法的有效性。

1 塊稀疏陣列信號模型與塊稀疏恢復算法

考慮水平分層柱對稱平坦海底的波導環境，建立水平均勻直線陣(Horizontal line array,HLA)淺海聲場模型。對分析頻率f，假設僅有I個遠程點聲源激發的聲場被HLA接收，其中，第i個點源場可近似為離散的M階簡正模線性疊加[11]：

其中，Si(f)為第i個點源頻譜，ri和zsi分別為對應的收發距離和聲源深度，krm為第m階簡正模水平波數，zr為接收器深度，ρ(zsi)為聲源處水介質密度，Ψm為模深度函數。假設HLA有L個間隔為d的陣元，第i個點源方位角為θi，以第1個陣元為參考陣元，ri表示第i個點聲源與參考陣元1的水平距離，第i個點聲源與參考陣元l的水平距離可表示為ril=ri+(l?1)dcosθi，假設其滿足遠場近似條件|ril?ri|?ri，則第l個陣元位置的聲場可近似為

其中，

為第i個點源激發的第m階模態，在HLA參考陣元1的位置的復數權值。公式(2)表明接收聲場是I個聲源激發的M階模態線性組合的結果，共有I×M項線性求和，其中聲源i激發的聲信號，傳播到參考陣元1位置的模態間頻散體現在Wim的ejkrmri項中，聲信號從陣元1傳播到陣元l的模態間頻散體現在ejkrm(l?1)dcosθi項中。假設已知接收深度的水介質聲速為c，則波數k=2πf/c。定義第m階簡正模的俯仰角?m，則sin?m=krm/k。簡正波陣列信號模型如圖1所示。

圖1 簡正波陣列信號模型Fig.1 Array signal model of normal mode

進一步，可將式(2)重寫為

其中，θim=arccos(cosθisin?m)，對方位角θi入射的點源，將會有M個等效的模態方位角{θim}與之對應。將方位角空間{θ}進行均勻離散化，假設離散網格數為G，不考慮離格問題，認為網格中包含I個入射點源方位角，除有限I個入射點源方位角外(I一般較小)，其他方位上均沒有聲源，因此信號是稀疏的，稀疏度為I；又因為第g個方位角對應有M個等效的模態方位角，因此信號是塊稀疏的。為方便討論，將式(3)寫成矩陣的形式：

其中，y=[P(1),P(2),···,P(L)]T為測量聲場向量，上標T表示轉置，A∈CL×GM為按塊結構組合的模態方位導向矢量矩陣(字典)：

其每一列為模態方位導向矢量(原子)

其G塊中僅有I塊元素非零，其余塊元素均為0，n∈CL×1為測量噪聲向量。此時，對淺海低頻聲源的方位估計問題，轉化為已知測量y和模態方位導向矢量矩陣A，對權值向量x進行塊稀疏恢復的問題。針對以上簡正波陣列信號模型中的塊稀疏恢復問題，給出BOMP聲源方位估計算法，算法偽代碼如表1所示。

表1 BOMP估計聲源方位算法Table 1 BOMP algorithm for azimuth estimation

2 數值仿真

根據2011年北黃海聲學實驗中一次實測的聲速剖面，并以介質均勻無限大液態海底作為仿真波導模型，如圖2所示。水介質中聲速接近等聲速c=1480 m/s，海底聲速設為cb=1650 m/s，海底深度設為H=69 m，水介質密度和海底密度分別設為ρw=1 g/cm3和ρb=1.8 g/cm3，水介質聲衰減系數和海底聲衰減系數分別設為αw=0 dB/λ和αb=0.1 dB/λ。仿真分析頻率取f=100 Hz，使用Kraken程序[12]計算波導模態，得到4階簡正模深度函數。

2.1 兩個點聲源的時間方位歷程圖

仿真兩個點聲源同時做5 m/s水平勻速直線運動，運動軌跡如圖3所示。聲源深度均為5.7 m，源譜Si(f=100 Hz)幅度均設為1，由于聲源運動速度相對水介質中聲速很小，所以不考慮多普勒頻移效應。同時接收陣元信噪比固定設為20 dB，即不隨聲源位置變化，根據圖2波導環境參數，仿真計算海底HLA接收聲場信號，HLA孔徑設為1 km，均勻陣元間隔5 m，陣元個數為201個，使用CBF和BOMP方法分別進行聲源的方位歷程估計，結果如圖4(a)、圖4(b)所示，其中點虛線和點劃線分別表示聲源1和聲源2的真實時間-方位歷程圖。圖4的結果表明，CBF在對淺海低頻聲源方位估計時會產生波束偏移和分裂現象，特別是在靠近端射方向上，而基于塊稀疏模型的BOMP算法可準確估計聲源方位。

圖2 仿真波導環境及100 Hz點源激發的模態深度函數Fig.2 The simulation waveguide environment and modal depth functions excited by a point source at 100 Hz

圖3 兩個點聲源的運動軌跡和水平陣的相對位置Fig.3 The tracks of two point-sources and the relative location of HLA

圖4 聲源方位歷程估計Fig.4 The estimation of source tracks by CBF and BOMP,respectively

2.2 信噪比的影響

為考察信噪比(Signal to noise ratio,SNR)對所提BOMP淺海聲源方位估計方法的影響，仿真只考慮一個聲源的情況，聲源與參考陣元的水平距離為10 km，聲源深度5.7 m，HLA參數與2.1小節相同，仿真聲源方位分別為0°、10°、20°和30°的陣列信號，加不同的陣元SNR的復高斯白噪聲，并進行1000次蒙特卡洛仿真實驗，角度搜索方位為0°～90°，統計CBF和所提BOMP估計的聲源方位的均值和標準差，對應誤差棒曲線如圖5(a)～圖5(d)所示。可以看出，在聲源方向靠近陣列端射方向時，CBF估計的方位角與真實的方位角存在較大的偏差，隨著SNR變大，標準差變小，偏差并沒有變小；隨著聲源方位角從陣列端射方向過渡到30°方向，CBF估計的聲源方位偏差逐漸減小。通過所提BOMP方法與CBF估計結果對比發現，在SNR大于?10 dB且聲源靠近陣列端射的情況下，BOMP估計的方位角與真實方位角存在較小的偏差，且隨著SNR的增大，標準差減小，在高信噪比的情況下，所提BOMP方法可以實現聲源方位的準確估計。在SNR小于?10 dB時，隨著SNR的降低，CBF和BOMP估計的聲源方位逐漸趨于0°～90°均勻分布，均值收斂于45°方向，估計的聲源方位沒有參考意義。由于陣元數為201，即陣增益為10lg201≈23(dB)，?10 dB的陣元域信噪比對應陣列波束輸出信噪比約為13 dB。實際應用時，當陣元數L已知，通過不等式SNR≥13 dB?10lgL可以判斷聲吶探測的海洋環境中陣元域信噪比是否符合本文方法的要求；當陣元域信噪比SNR(即輸入信噪比)已知，通過不等式10lgL≥13 dB?SNR以檢驗聲吶的線列陣陣元數L是否滿足使用本文方法的要求。

圖5 蒙特卡洛仿真實驗1000次，統計CBF和BOMP估計的聲源方位均值和標準差Fig.5 The mean and the standard deviation of the estimated source azimuth by CBF and BOMP,with 1000 times of Monte Carlo simulation experiment for each SNR

3 實驗數據處理

實驗數據源自2011年北黃海冬季聲學實驗，實驗船拖曳著深度5.7 m氣槍聲源在準南北方向航行，航速約為4 kn，氣槍間隔為1 min發射一次脈沖信號，接收陣列為海底32元準南北放置的HLA，有效陣元數為28，陣列水平孔徑約250 m，通過8組鄰近的氣槍聲信號進行孔徑擴展處理[13?14]，得到水平跨度為1 km、均勻陣元間隔5 m的201元HLA接收的聲場信號。

對頻率間隔1/3 Hz的50～200 Hz寬帶聲場信號(離散的頻點數為451)，進行逐頻點聲場的CBF和BOMP方位估計，其中，角度離散網格間隔取0.5°，觀察角度范圍取0°～40°，最大值歸一化后波束輸出結果如圖6所示。圖6(a)為CBF估計的方位譜結果，可觀察到明顯的波束偏移和分裂現象，每個波束峰值點并不表示真實的聲源方位，由于真實的聲源方位接近HLA端射方向，即方位角約為0°，分裂的波束峰對應的角度與簡正模俯仰角近似互為余角關系。基于圖2所示的波導模型，使用BOMP進行逐頻點的聲源方位估計，結果如圖6(b)所示，可以看出，在大部分頻點上，BOMP估計的波束輸出在靠近端射方向有較強的能量。表2給出了對應451個頻點估計方位角的離散統計分布。估計聲源方位為0°的占比大于90%，在某些頻點上出現較大的方位估計偏差，其原因包括某些頻點上信噪比不足、存在干擾源影響和模型失配。圖7給出了50～200 Hz寬帶內波束能量求和，最大值歸一化輸出的結果，其中CBF波束輸出旁瓣較寬，在17.5°方向上能量最大，該方向與端射方向存在較大的偏差；BOMP波束輸出在0°方向能量最大，其他方位的能量只有很小的波動，與表2統計451個頻點的方位估計為0°的占比大于90%結果相符，估計方位與實際聲源方位一致，實驗數據驗證了所提BOMP方法能準確估計淺海低頻寬帶聲源的方位。

圖6 實驗氣槍聲源方位譜估計Fig.6 The experimental air-gun source azimuth spectra estimation by CBF and BOMP,respectively

圖7 實驗氣槍50～200 Hz寬度信號聲源方位估計Fig.7 The experimental air-gun 50～200 Hz broadband source azimuth spectra estimation

表2 統計451個頻點BOMP估計的方位離散分布Table 2 Discrete distribution of source azimuth by BOMP with 451 frequency points

4 結論與討論

本文基于淺海水平分層波導的聲傳播簡正波理論，建立了塊稀疏信號模型，提出一種塊正交匹配追蹤(BOMP)的淺海水平陣低頻聲源方位估計方法。仿真結果表明，所提方法在陣列波束輸出信噪比大于13 dB和無失配環境情況下，能規避常規波束形成(CBF)方法在估計近陣列端射方向的淺海低頻聲源方位時所產生的波束偏移和分裂現象，從而克服較大方位估計偏差的問題。對2011年北黃海聲學實驗水平長線陣接收低頻聲信號，進行逐頻點方位估計以及能量求和波束輸出。CBF估計的方位譜中可以觀察到明顯的波束偏移和分裂現象，寬帶能量求和波束輸出旁瓣較寬，且波束峰值方位與真實聲源方位存在較大的偏差。本文所提方法估計的聲源方位除部分頻點估計結果與端射方向偏離較大外，大部分頻點上方位估計的結果均靠近端射方向；BOMP波束寬帶能量輸出，具有超低旁瓣的特點，且在0°方向取得能量值最大，與真實聲源方位一致。仿真和實驗數據驗證了所提方法的有效性和較CBF的優越性。

本文所提BOMP算法需要已知波導環境參數，用于計算模態方位導向矢量矩陣字典，當環境參數存在失配，特別是在較高頻段，聲源方位的估計偏差會變大，在實際應用中需要選擇合適的等效波導環境參數。另外，本文所提方法針對大孔徑基陣近端射聲源方位估計中存在的波束分裂和方位偏差問題，只適用于較高信噪比的水平陣信號，當信噪比不足時，估計的聲源方位會出現較大的偏差。

致謝感謝參加2011年北黃海冬季聲學實驗的全體工作人員，他們的辛苦工作為本文提供了寶貴的實驗數據。