深海聲影區時頻譜干涉結構與聲源定位?

劉與涵 郭良浩 章偉裕 閆 超 董 閣

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

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

0 引言

深海聲源定位是海洋聲學的重要問題之一。對于近水面深度(1000 m 以淺)的聲源和接收點，由于水體聲速剖面的折射效應，在幾千米至幾十千米的水平距離范圍內通常會形成海底反射區。該區域內只有經海底反射的聲線到達，由于缺少直達聲線和水體反轉聲線，聲傳播損失較大，也稱為聲影區。聲影區的水平覆蓋范圍較大，是中近程水聲探測的主要區域。

匹配場處理方法較早被應用于深海聲源定位。Fizell 等[1]利用垂直陣定位了深海中的低頻聲源。周士弘等[2]將廣義相積分簡正波方法用于匹配場定位，以提高深?？截悎龅挠嬎阈?。陳連榮等[3]將高斯射線束方法用于深海匹配場定位。但是，匹配場處理方法對于海底底質等環境參數的失配問題比較敏感，且通常需要較大的垂直接收孔徑以提高定位精度。

為了提高定位穩健性以及減小接收孔徑，國內外學者開展了大量深海聲源定位方法研究。Tiemann 等[4]利用單水聽器接收信號提取多途時延差，對抹香鯨發出的寬帶脈沖聲源進行定位。楊士莪[5]通過小型矢量立體陣接收的多途到達角，聯合解算聲源方位、距離和深度。McCargar 等[6]通過垂直陣的窄帶信號波束輸出能量隨距離變化的明暗干涉周期，估計直達聲區內的聲源深度。王文博等[7]進一步利用寬帶聲場的頻率-掠射角干涉結構，估計直達聲區內的聲源深度。

在聲影區中，聲場由經一次海底反射的多途聲線主導，這些聲線的到達俯仰角和時延差等多途特征可利用小孔徑陣列獲取。此外，上述多途參數受海底底質的影響較小，且由于海底反射聲線的掠射角較大，聲速剖面的不確定性對其影響也較小。因此，利用聲影區的多途特征可實現比較穩健的聲源定位。Duan 等[8]利用近水面深度的垂直短陣進行窄帶觀測，通過加權子空間匹配方法進行聲源定位。吳禹沈等[9]利用水下滑翔機接收脈沖聲信號，通過多途到達時延差進行聲源定位。謝亮等[10]通過垂直陣提取脈沖簇到達結構，進行水下聲源定位。劉與涵等[11]利用短間距垂直雙水聽器接收窄帶寬信號，通過頻譜搬移和稀疏時延估計方法提取多途到達時延差，估計聲源距離。

此外，聲影區的聲場強度在距離-頻率維通常會形成比較穩健的干涉結構。Li 等[12]指出，聲影區中的波導不變量(β)值接近1，聲場的干涉條紋斜率與淺海比較相近。翁晉寶等[13]研究了聲影區中的距離-頻率干涉結構，發現存在兩種干涉條紋，第一種條紋的頻率周期與水平距離和聲源深度有關，第二種條紋的頻率周期與水平距離和接收深度有關。吳俊楠等[14]進一步指出，水面船輻射噪聲形成的第一種條紋在較窄的頻帶內可忽略，通過水平陣波束輸出信號的自相關函數，可有效提取第二種條紋的頻率周期，估計水面聲源的距離。翁晉寶等[15-16]也利用單水聽器接收的時頻譜干涉結構，驗證了上述定位方法的有效性。

深海聲影區的傳播損失較大，在實際的水聲探測中，隨著接收信噪比(Signal-to-noise ratio,SNR)降低，現有的聲影區定位方法的性能逐漸下降或失效。理論上，對接收信號時頻譜進行二維傅里葉變換(Two-dimensional Fourier transform,2-D FT)后處理，可將變換前能量分散的干涉條紋聚焦在變換后能量集中的較小區域內，以提高低SNR 條件下的干涉結構檢測能力[17]。但是，現有2-D FT 研究主要關注于淺海環境，針對深海聲影區的有關研究較少。淺海聲場通常具有多組斜率近似相同的干涉條紋，對應的β值在1 附近，因此2-D FT 幅度譜通常呈現一條“脊”狀分布[18]?；谶@一特性，Cockrell等[19]利用2-D FT估計淺海聲場的干涉條紋斜率，從而解算聲源距離。相比之下，深海聲影區中的β分布復雜多變，影區聲場的特點在于干涉條紋種類較少，對應的2-D FT 幅度譜通常呈“點”狀分布，不能直接沿用淺海聲源定位方法。因此，有必要進一步研究影區干涉條紋及其2-D FT 幅度譜的分布規律與聲源位置的關系，再通過2-D FT 方法提取上述幅度譜分布結構，以實現低SNR條件下的深海聲源定位。

本文基于射線理論，在研究深海聲影區中水面運動聲源形成的時頻譜干涉結構的基礎上，建立了頻率干涉周期與聲源距離的關系、時間干涉周期與聲源徑向速度的關系。通過時頻譜的2-D FT 幅度譜提取上述兩種時頻干涉周期，再通過多頻帶處理方法提高低SNR 下的檢測能力。數值仿真及海試數據處理結果表明，相比于現有的利用接收信號自相關的聲源測距方法，所提出的利用時頻譜2-D FT的定位方法得到的距離和徑向速度估計精度較高，且測距連續性較好，比較適用于低SNR條件下的聲源定位。

1 聲影區時頻譜干涉結構

1.1 干涉條紋形成機理

根據射線聲學理論，對于圖1(a)所示的典型深海聲速剖面，直達聲線和水體反轉聲線無法到達聲影區。在第一影區中，經二次及以上海底反射的聲線傳播損失較大，對聲場的影響通?？珊雎?。經一次海底反射的4 條多途聲線對聲場能量起主要貢獻，分別為海底反射(BR)聲線、海面-海底反射(SBR)聲線、海底-海面反射(BSR)聲線和海面-海底-海面反射(SBSR)聲線，其傳播軌跡如圖1(b) 所示。

圖1 典型深海聲速剖面和第一影區中的聲線傳播軌跡Fig.1 A typical sound speed profile in deep water and the propagation paths of the acoustic rays in the first shadow zone

位于接收點(r,zr)、頻率為f的聲壓可表示為上述多途聲線的相干疊加：

其中，r表示水平距離，zr表示接收深度，ti(i=1,2,3,4)分別表示BR、SBR、BSR 和SBSR 聲線的到達時延。

因此，聲場強度可近似表示為[13]

其中，τ12≡t2-t1為SBR聲線與BR聲線的到達時延差，τ13≡t3-t1為BSR聲線與BR聲線的到達時延差。需要說明的是，SBSR 聲線與BSR 聲線的到達時延差τ34≡t4-t3與τ12近似相等，在推導過程中由τ12替換，因此式(2)省略了與SBSR 聲線有關的時延項。

對于水面船輻射噪聲，聲源深度通常僅為幾米，到達接收點的SBR 聲線與BR 聲線的傳播軌跡十分接近，導致τ12通常小于幾毫秒。對于500 Hz以下的低頻聲場，f與τ12的乘積通常小于5。因此，在較窄的頻帶和較短的觀測時間內，式(2)中的cos 2πfτ12項可近似為常數。與之相比，水聽器接收深度通常數倍于聲源深度，導致τ13也數倍于τ12。因此，式(2)中cos 2πfτ13項的周期變化會形成聲場強度的干涉現象。當τ13=n/f，n∈Z+時，聲強出現干涉相消；當τ13=(m+1/2)/f，m∈Z+時，聲強出現干涉相長。

在實際的被動觀測中，接收深度zr通常在一定時間內保持恒定，聲源與接收點之間的水平距離r隨時間t變化，時延差τ13(t)也隨時間變化。接收信號的時頻譜可近似表示為

1.2 頻率干涉周期

由式(3)可知，在t時刻，接收信號時頻譜沿頻率軸的干涉周期為?f=1/τ13(t)。文獻[13]推導了到達時延差τ13的近似表達式

其中，c0為聲源處的聲速，n(z)為深度z處的折射率，α3為BSR 聲線的出射角。例如，在等聲速近似條件下，α3(t)=arctan{(2H+zr)/r(t)}，其中H為海深?？梢钥闯靓?隨水平距離r的增大而單調遞減，因此在接收深度zr恒定的前提下，τ13隨r的增大而單調遞減，頻率干涉周期?f隨r的增大而單調遞增。

圖2 為聲影區內頻率干涉周期?f的值隨水平距離的變化曲線。由于?f隨水平距離單調變化，該特征量可用于聲源距離估計。

圖2 頻率干涉周期隨水平距離的變化曲線Fig.2 The frequency interference period as a function of the horizontal range

1.3 時間干涉周期

由式(3)可知，對于頻率f，接收信號時頻譜沿時間軸的干涉周期?t可表示為

假設在時間間隔?t內，聲源運動的徑向速度近似恒定為vr，BR 聲線和BSR 聲線的到達俯仰角分別近似恒定為θ1(t),θ3(t)，則有：

其中，cr為接收點處的海水聲速。

將式(7)代入式(6)，得到頻率f處的時間干涉周期?t的表達式

圖3 為200 Hz 頻點處?t的值隨水平距離的變化曲線，其中vr=10 m/s。盡管式(8)中的cosθ1(t)-cosθ3(t)項導致?t與水平距離呈非單調關系，但在距離已知的條件下可計算cosθ1(t)-cosθ3(t)的值，從而估計聲源的徑向速度vr。

圖3 時間干涉周期隨水平距離的變化曲線Fig.3 The time interference period as a function of the horizontal range

1.4 時頻譜的2-D FT

取一段時頻譜窗Iwin(f,t)，中心時刻為t0，時長為T，中心頻率為f0，帶寬為B。如圖4(a)所示，當選取的時長和帶寬較小時，窗內的頻率干涉周期?f和時間干涉周期?t近似恒定。經去直流處理后，該時頻譜窗的2-D FT幅度譜為

圖4 時頻譜及其2-D FT 示意圖Fig.4 The spectrogram and its 2-D FT illustration

其中，?表示卷積，sinc(·)表示sinc 函數，δ(·)表示Dirac 函數式(9)中，橫坐標kf的單位為秒(s)，縱坐標kt的單位為赫茲(Hz)。

由于式(9)關于原點對稱，以下僅討論kf>0對應的I?win(kf,kt)第I、第IV象限。式(9)表明，時頻譜窗的2-D FT 幅度譜中，在位置將出現幅度最大值，如圖4(b)中白色十字標記位置所示。

2 聲源距離和徑向速度估計

由第1 節分析可知，水面聲源在聲影區激發的聲場干涉結構與水平距離和徑向速度有關。對一段較小的接收信號時頻譜窗進行2-D FT，可將變換前能量分散的干涉條紋聚焦在變換后能量集中的較小區域內。檢測2-D FT 幅度譜中最大值點對應的橫縱坐標，即可提取頻率干涉周期?f和時間干涉周期?t，從而解算聲源距離和徑向速度。上述分析的前提條件是時頻譜窗的處理時長和帶寬較小，窗內的干涉條紋近似為相互平行的直線，而對于較大帶寬的接收信號，可進行多頻帶處理，處理方法如下。

考慮一段接收信號，中心時刻為t0，時長為T。將接收信號時頻譜分割為Q個具有一定重疊率的子頻帶，其帶寬均為B，中心頻率分別為f(1),f(2),···,f(Q)。對上述Q個子頻帶分別進行2-D FT 處理，相應的幅度譜最大值點橫坐標記為，縱坐標記為

由式(4)可知，上述Q個橫坐標相等，且與聲源距離有關：

由式(8)可知，上述Q個縱坐標與對應子頻帶的中心頻率成正比，且與聲源距離和徑向速度有關：

然后，對上述Q個子頻帶的2-D FT 幅度譜沿縱軸(kt方向)進行線性縮放并累加，得到多頻帶處理的2-D FT幅度譜：

式(12)中，橫坐標kf的單位為秒(s)，縱坐標kt/fcenter無量綱。

本文將橫坐標的界限設定為

其中，τ13,min、τ13,max分別表示τ13的最小值和最大值，由射線模型計算得到；引入±1/B項的作用是包含sinc函數沿橫軸的主瓣寬度。

本文將縱坐標的界限設定為

其中，vr,max為目標船的最大徑向速度，(cosθ1-cosθ3)max由射線模型得到；引入±1/T項的作用是包含sinc函數沿縱軸的主瓣寬度。

最后，通過Bellhop 射線追蹤模型計算不同聲源距離下τ13(r)和cosθ1(r)-cosθ3(r)的拷貝值。根據式(10)，將實際提取的與拷貝值τ13(r)進行匹配，估計聲源的水平距離r(t0)；再根據式(11)，計算聲源的徑向速度vr。

3 數值仿真

本節通過數值仿真，驗證所提出方法在低SNR條件下的有效性。

仿真環境中的海水聲速剖面如圖1(a)所示，海底密度為1.8 g/cm3，海底聲速為1600 m/s。聲源深度為7 m，以10 m/s的徑向速度從水平距離3.6 km處逐漸遠離至54 km 處，運動范圍基本覆蓋了第一影區的全部距離。接收深度為150 m，觀測頻帶為50～500 Hz，接收信號混有加性高斯白噪聲。下面分別在接收SNR 0 dB 和-15 dB 條件下進行數值仿真。

圖5(a)和圖5(b)分別為接收SNR 0 dB 和-15 dB 條件下的時頻譜，在-15 dB 條件下難以直接觀察到干涉條紋，而在0 dB 條件下仿真結果呈現出比較明顯的兩種干涉條紋。其中，頻率周期較大的條紋與式(2)中的τ12項有關，對于后續2-D FT 處理所選用的時頻譜窗理論上可忽略，在現有的水面船輻射噪聲數據中通常也觀察不到此類條紋[14-16]。頻率周期較小的條紋與式(2)中的τ13項有關，可用作聲源距離和徑向速度估計。

圖5 仿真數據接收信號時頻譜Fig.5 The spectrograms of the received signals from the simulations

利用所提出的2-D FT 方法分別對上述兩種SNR條件的時頻譜進行干涉結構檢測。時頻譜窗的帶寬為100 Hz，頻帶重疊率為95%，時長為3 min，時間重疊率為90%。

圖6 仿真數據中在最大值點沿自變量kf 和kt/fcenter 切片的時間歷程Fig.6 The time records of slices along the independent variables kf and kt/fcenter centered on the maximum points of from the simulations

圖7 仿真數據中水平距離和徑向速度估計結果Fig.7 The estimated results of the horizontal range and the radial velocity from the simulations

在接收SNR -15 dB 條件下，盡管難以直接觀察到時頻譜干涉條紋，但利用所提出的2-D FT 方法可有效提取所需的干涉結構。沿kf方向切片的時間歷程如圖6(b)所示，對應的聲源距離估計結果如圖7(b)所示，測距結果出現了少量野點，但在大部分時間正確估計了聲源距離。沿kt/fcenter方向切片的時間歷程如圖6(d)所示，對應的徑向速度估計結果如圖7(d)所示，除少量野點外，徑向速度估計結果在前80 min 內的誤差小于3 m/s，在第80 min后偏離真實距離。

在低SNR 的仿真實驗中，所提方法的檢測增益一方面來自2-D FT對干涉條紋能量的聚焦作用，另一方面來自較大的頻帶重疊率。假設背景噪聲在整個處理頻帶內服從高斯白噪聲分布，那么對于每個子頻帶，2-D FT 前后的背景噪聲也均滿足高斯白噪聲特性。理論上，頻帶重疊率越高，子頻帶數目N就越多，公式(12)中的求和項也越多，在一定SNR 條件下可有效積累這些求和項的譜峰能量，從而提高檢測增益。但事實上，由于仿真中的重疊頻帶取自于同一個樣本源，這些重疊頻帶內的背景噪聲不具有獨立性。因此，增大頻帶重疊率不一定能有效增加獨立觀測次數，對檢測增益的影響比較有限。

下面通過蒙特卡洛仿真實驗，進一步分析接收SNR 和頻帶重疊率對測距性能的影響。仿真環境、收發深度、徑向速度、時頻譜窗的帶寬和頻帶重疊率均與以上實驗保持一致，聲源從水平距離9.1 km 處遠離至10.9 km 處(水平距離近似固定為r0=10 km、時長為3 min)，進行100 次蒙特卡洛仿真，然后統計測距結果的均方根誤差(Root mean squared error,RMSE)。RMSE定義為

對于-14～-10 dB 的接收SNR 范圍，不同頻帶重疊率條件下RMSE 隨接收SNR 變化的關系曲線如圖8(a)所示。首先，在相同的接收SNR 條件下，頻帶重疊率越高，測距結果的RMSE越低，不同頻帶重疊率對應的3條曲線在-14 dB處趨于接近。其次，對于相同的頻帶重疊率，測距結果的RMSE與SNR 呈負相關，且RMSE 均小于0.3 km，單次測量結果即可有效實現測距。對于-16～-10 dB 的接收SNR范圍，圖8(a)擴展為圖8(b)。當SNR小于-14 dB 時，隨著SNR 繼續降低，不同頻帶重疊率對應的測距結果RMSE 同步顯著增大，的單次估計結果不再能穩健地實現聲源測距。但本節在-15 dB 仿真條件下，仍可通過的時間歷程(圖6(b)所示)或測距結果的點跡(圖7(b)所示)觀察出聲源距離的變化軌跡。上述結果表明，增大頻帶重疊率能帶來一定的處理增益，但對測距精度的影響比較有限，影響測距精度的主要因素是SNR。

圖8 不同頻帶重疊率條件下RMSE 隨接收SNR變化的關系曲線Fig.8 The relationship between the RMSE and the received SNR at difference overlapping ratios of the frequency band

以上分析中，重疊頻帶內的背景噪聲不具有獨立性。如果通過其他方法增加獨立觀測次數或增大等效帶寬，則有機會進一步提高檢測增益。此外，以上結論要求處理頻帶內的背景噪聲滿足高斯白噪聲特性。對于實際的海試接收信號時頻譜，背景噪聲具有一定的空間分布特性，使得檢測增益和測距性能通常低于仿真結果。

4 海試數據處理與分析

2022 年4 月，在菲律賓海開展了一次深海聲學實驗。實驗站點的海深約4800 m，根據World Ocean Atlas(WOA)數據庫的歷史溫度[20]、鹽度[21]數據計算得到的海水聲速剖面如圖9所示。

圖9 實驗期間的海水聲速剖面Fig.9 The sound speed profile during the experiment

實驗期間，一艘非合作的水面船逐漸靠近接收點而后駛離，在接收點利用近水面深度布放的水聽器對其輻射噪聲進行觀測。其中，較高SNR 的接收信號時頻譜如圖10(a)所示，較低SNR 的接收信號時頻譜如圖10(b)所示。由于目標船的徑向速度發生了多次改變，時頻譜中的干涉條紋出現了多個拐點。此外，時頻譜中不同區域的干涉條紋清晰程度不同，背景噪聲不能簡單假定為加性高斯白噪聲。

圖10 海試數據接收信號時頻譜Fig.10 The spectrograms of the received signals from the experimental data

在圖10(a)中，分別取兩段時頻譜窗，如圖11(a)和圖11(b)所示，對應的2-D FT 幅度譜分別如圖11(c)和圖11(d)所示。其中，圖11(a)中的背景噪聲更加接近高斯白噪聲，因此圖11(c)中的噪聲分布也比較均勻，由干涉結構形成的幅度譜最大值出現在紅色虛線框(第2節提出的橫縱坐標上下界)之內，絕大部分噪聲則分布在紅框之外。與之相比，圖11(b)中的背景噪聲和條紋能量具有一定的空間分布特性，0.625～0.75 歸一化頻帶內的條紋能量明顯高于0.5～0.625 頻帶，導致圖11(d)在坐標原點附近形成了幅度較大的偽峰，其幅度大于紅框內由干涉結構形成的譜峰。因此，直接搜索2-D FT幅度譜的最大值點可能會錯誤地提取所需干涉周期。所幸的是，通過第2 節提出的橫縱坐標上下界，類似于對時頻譜窗進行二維圖像帶通濾波，利用橫坐標下界可濾除圖11(d)中的偽峰。

圖11 圖10(a)中兩段時頻譜窗和2-D FT 幅度譜Fig.11 Two segments of the spectrogram in Fig.10(a) and their 2-D FT magnitude spectrums

事實上，由背景噪聲形成的譜峰有可能出現在分析的上下界內，但仍可通過譜峰的相對位置以及多頻帶處理后的譜峰能量加以判別。首先，受背景噪聲和條紋能量分布的影響，觀測到的條紋數量相比于原有數量通常減少而不是增加，因此2-D FT幅度譜中的偽峰通常比干涉結構形成的譜峰更加靠近原點。其次，式(12)將各個頻帶內的2-D FT幅度譜進行累加，不同頻帶形成的偽峰通常會出現在不同位置，因此偽峰能量在累加過程中逐漸被平均掉。綜上所述，當2-D FT幅度譜中出現兩個或以上譜峰時，可初步判定靠近原點的為偽峰，在多頻帶累加后的2-DFT 幅度譜中，可進一步判定能量較低的為偽峰。以上分析針對目標船輻射噪聲的連續譜成分，當存在強線譜成分時，2-D FT 幅度譜中可能會出現其他性質的譜峰，有待進一步研究。

利用所提出的2-D FT方法分別檢測圖10中兩張時頻譜中的干涉結構，時頻譜窗所用的帶寬、時長和重疊率參數均與仿真實驗相同。

圖12 海試數據中在最大值點沿自變量kf 和kt/fcenter 切片的時間歷程Fig.12 The time records of slices along the independent variables kf and kt/fcenter centered on the maximum points of from the experimental data

圖13 海試數據的水平距離和徑向速度估計結果Fig.13 The estimated results of the horizontal range and the radial velocity from the experimental data

其次，為了對比所提出的2-D FT 方法在低SNR 下的優越性，利用單個時刻的接收信號自相關方法提取?f，得到的測距結果如圖13(a)和圖13(b)中的藍色點跡所示。在較高SNR 條件下，自相關方法的測距結果能夠體現距離變化的總體趨勢，但野點較多，實際中需要后處理方法以消除野點。相比之下，2-D FT 方法的野點較少，測距結果的穩健性較高。在較低SNR條件下，自相關方法的測距野點過多，難以有效提取聲源距離，而2-D FT方法仍然在大部分時間段實現了測距。

5 結論

利用射線聲學理論，分析了深海聲影區中接收信號時頻譜的干涉結構，給出了頻率干涉周期與聲源距離的近似表達式，以及時間干涉周期與徑向速度的近似表達式。利用水面船輻射噪聲信號的時頻譜2-D FT 幅度譜，提高時頻干涉周期的檢測增益，實現了低SNR 條件下聲源距離和徑向速度的聯合估計。最后通過數值仿真和海試數據處理，驗證了所提方法的有效性。結果表明，在低SNR 條件下，所提方法相比于自相關方法具有較好的穩健性。下一步工作包括研究水面船輻射噪聲中的強線譜成分所形成的干涉結構。

致謝感謝參與海上實驗的全體人員，他們為研究工作提供了寶貴的實驗數據。感謝中國科學院聲學研究所的任歲玲副研究員為實驗數據處理提供了重要幫助。