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

(1 中國(guó)科學(xué)院聲學(xué)研究所 聲場(chǎng)聲信息國(guó)家重點(diǎn)實(shí)驗(yàn)室 北京 100190)

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

0 引言

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

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

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

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

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

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

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

其中，

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

圖1 簡(jiǎn)正波陣列信號(hào)模型Fig.1 Array signal model of normal mode

進(jìn)一步，可將式(2)重寫(xiě)為

其中，θim=arccos(cosθisin?m)，對(duì)方位角θi入射的點(diǎn)源，將會(huì)有M個(gè)等效的模態(tài)方位角{θim}與之對(duì)應(yīng)。將方位角空間{θ}進(jìn)行均勻離散化，假設(shè)離散網(wǎng)格數(shù)為G，不考慮離格問(wèn)題，認(rèn)為網(wǎng)格中包含I個(gè)入射點(diǎn)源方位角，除有限I個(gè)入射點(diǎn)源方位角外(I一般較小)，其他方位上均沒(méi)有聲源，因此信號(hào)是稀疏的，稀疏度為I；又因?yàn)榈趃個(gè)方位角對(duì)應(yīng)有M個(gè)等效的模態(tài)方位角，因此信號(hào)是塊稀疏的。為方便討論，將式(3)寫(xiě)成矩陣的形式：

其中，y=[P(1),P(2),···,P(L)]T為測(cè)量聲場(chǎng)向量，上標(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為測(cè)量噪聲向量。此時(shí)，對(duì)淺海低頻聲源的方位估計(jì)問(wèn)題，轉(zhuǎn)化為已知測(cè)量y和模態(tài)方位導(dǎo)向矢量矩陣A，對(duì)權(quán)值向量x進(jìn)行塊稀疏恢復(fù)的問(wèn)題。針對(duì)以上簡(jiǎn)正波陣列信號(hào)模型中的塊稀疏恢復(fù)問(wèn)題，給出BOMP聲源方位估計(jì)算法，算法偽代碼如表1所示。

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

2 數(shù)值仿真

根據(jù)2011年北黃海聲學(xué)實(shí)驗(yàn)中一次實(shí)測(cè)的聲速剖面，并以介質(zhì)均勻無(wú)限大液態(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]計(jì)算波導(dǎo)模態(tài)，得到4階簡(jiǎn)正模深度函數(shù)。

2.1 兩個(gè)點(diǎn)聲源的時(shí)間方位歷程圖

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

圖2 仿真波導(dǎo)環(huán)境及100 Hz點(diǎn)源激發(fā)的模態(tài)深度函數(shù)Fig.2 The simulation waveguide environment and modal depth functions excited by a point source at 100 Hz

圖3 兩個(gè)點(diǎn)聲源的運(yùn)動(dòng)軌跡和水平陣的相對(duì)位置Fig.3 The tracks of two point-sources and the relative location of HLA

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

2.2 信噪比的影響

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

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

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

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

圖6 實(shí)驗(yàn)氣槍聲源方位譜估計(jì)Fig.6 The experimental air-gun source azimuth spectra estimation by CBF and BOMP,respectively

圖7 實(shí)驗(yàn)氣槍50～200 Hz寬度信號(hào)聲源方位估計(jì)Fig.7 The experimental air-gun 50～200 Hz broadband source azimuth spectra estimation

表2 統(tǒng)計(jì)451個(gè)頻點(diǎn)BOMP估計(jì)的方位離散分布Table 2 Discrete distribution of source azimuth by BOMP with 451 frequency points

4 結(jié)論與討論

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

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

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