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

梁玉權(quán) 周士弘 宮在曉

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

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

0 引言

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

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

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

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

1 塊稀疏陣列信號模型與塊稀疏恢復(fù)算法

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

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

其中，

為第i個點源激發(fā)的第m階模態(tài)，在HLA參考陣元1的位置的復(fù)數(shù)權(quán)值。公式(2)表明接收聲場是I個聲源激發(fā)的M階模態(tài)線性組合的結(jié)果，共有I×M項線性求和，其中聲源i激發(fā)的聲信號，傳播到參考陣元1位置的模態(tài)間頻散體現(xiàn)在Wim的ejkrmri項中，聲信號從陣元1傳播到陣元l的模態(tài)間頻散體現(xiàn)在ejkrm(l?1)dcosθi項中。假設(shè)已知接收深度的水介質(zhì)聲速為c，則波數(shù)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個等效的模態(tài)方位角{θim}與之對應(yīng)。將方位角空間{θ}進行均勻離散化，假設(shè)離散網(wǎng)格數(shù)為G，不考慮離格問題，認(rèn)為網(wǎng)格中包含I個入射點源方位角，除有限I個入射點源方位角外(I一般較小)，其他方位上均沒有聲源，因此信號是稀疏的，稀疏度為I；又因為第g個方位角對應(yīng)有M個等效的模態(tài)方位角，因此信號是塊稀疏的。為方便討論，將式(3)寫成矩陣的形式：

其中，y=[P(1),P(2),···,P(L)]T為測量聲場向量，上標(biāo)T表示轉(zhuǎn)置，A∈CL×GM為按塊結(jié)構(gòu)組合的模態(tài)方位導(dǎo)向矢量矩陣(字典)：

其每一列為模態(tài)方位導(dǎo)向矢量(原子)

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

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

2 數(shù)值仿真

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

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

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

圖2 仿真波導(dǎo)環(huán)境及100 Hz點源激發(fā)的模態(tài)深度函數(shù)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參數(shù)與2.1小節(jié)相同，仿真聲源方位分別為0°、10°、20°和30°的陣列信號，加不同的陣元SNR的復(fù)高斯白噪聲，并進行1000次蒙特卡洛仿真實驗，角度搜索方位為0°～90°，統(tǒng)計CBF和所提BOMP估計的聲源方位的均值和標(biāo)準(zhǔn)差，對應(yīng)誤差棒曲線如圖5(a)～圖5(d)所示?？梢钥闯觯诼曉捶较蚩拷嚵卸松浞较驎r，CBF估計的方位角與真實的方位角存在較大的偏差，隨著SNR變大，標(biāo)準(zhǔn)差變小，偏差并沒有變??；隨著聲源方位角從陣列端射方向過渡到30°方向，CBF估計的聲源方位偏差逐漸減小。通過所提BOMP方法與CBF估計結(jié)果對比發(fā)現(xiàn)，在SNR大于?10 dB且聲源靠近陣列端射的情況下，BOMP估計的方位角與真實方位角存在較小的偏差，且隨著SNR的增大，標(biāo)準(zhǔn)差減小，在高信噪比的情況下，所提BOMP方法可以實現(xiàn)聲源方位的準(zhǔn)確估計。在SNR小于?10 dB時，隨著SNR的降低，CBF和BOMP估計的聲源方位逐漸趨于0°～90°均勻分布，均值收斂于45°方向，估計的聲源方位沒有參考意義。由于陣元數(shù)為201，即陣增益為10lg201≈23(dB)，?10 dB的陣元域信噪比對應(yīng)陣列波束輸出信噪比約為13 dB。實際應(yīng)用時，當(dāng)陣元數(shù)L已知，通過不等式SNR≥13 dB?10lgL可以判斷聲吶探測的海洋環(huán)境中陣元域信噪比是否符合本文方法的要求；當(dāng)陣元域信噪比SNR(即輸入信噪比)已知，通過不等式10lgL≥13 dB?SNR以檢驗聲吶的線列陣陣元數(shù)L是否滿足使用本文方法的要求。

圖5 蒙特卡洛仿真實驗1000次，統(tǒng)計CBF和BOMP估計的聲源方位均值和標(biāo)準(zhǔn)差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 實驗數(shù)據(jù)處理

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

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

圖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 統(tǒng)計451個頻點BOMP估計的方位離散分布Table 2 Discrete distribution of source azimuth by BOMP with 451 frequency points

4 結(jié)論與討論

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

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

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