基于多線譜干擾抑制的水下對空運動聲源線譜探測

劉凱悅，彭朝暉，張靈珊，王光旭，李晟昊

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

0 引言

自20世紀50年代起，國內外學者開展了一系列有關空中運動聲源透射入水問題的研究[1-2]，證明了水面附近的空氣聲源可等效為聲源正下方十分接近水面的水下聲源。因此，通常的水下目標探測方法理論上可用于水下對空探測。

但是，由于空氣、海水介質特性阻抗存在巨大差異，空氣中的聲能量只有極少部分能夠透射入水，加之海洋環境背景噪聲的影響，傳統的寬帶能量檢測器從水下探測空中聲源相對困難。

空中聲源的輻射噪聲譜成分包括寬帶連續譜和離散線譜，其中線譜強度高且穩定性好[3-5]，因此利用線譜探測是一種可能的途徑。但是，由于空中聲源的運動速度普遍較高，導致目標方位變化快、接收聲場多普勒頻偏大，因此水下目標探測技術相繼產生的多種對線譜目標的探測方法[6-9]，如基于瞬時頻率方差加權、基于線譜方位穩定性的目標檢測方法等均難以應用于水下對空運動目標的線譜探測。

陳敬軍等[10]提出一種基于人工智能的線譜自動檢測算法，充分利用線譜形狀特點及其在多個時刻的頻率信息，根據最近時刻線譜頻率外推得到下一時刻線譜搜索中心，有效解決了存在較大多普勒頻移的線譜提取問題。文獻[11]借助上述線譜檢測算法，初步實現了水下對空中高速運動目標的遠程探測。

實際淺海環境中，從水下探測空中目標會不可避免地受到來自多個水下源的干擾，其中大型水面船占主要成分。通常情況下，這些水面船的線譜基頻很小[12]，在低頻范圍內具有較高的線譜能量，其線譜數量遠多于空中目標激發的水下線譜，甚至存在同頻干擾即目標與干擾的頻率相同或十分相近的情況。因此，單純依靠線譜檢測方法難以有效地探測到空中弱線譜目標，還需要進一步進行干擾抑制。

國內外學者在對強干擾抑制的研究方面已總結出許多成熟的理論和技術，包括基于波束形成的干擾抑制算法[13-15]、基于信號分離技術的干擾抑制算法[16]以及基于子空間類的自適應干擾抑制算法[17-18]等。前兩種算法分別在干擾與目標位于同一波束內、要求實時處理的情況下性能受損，相比之下，基于子空間類的干擾抑制算法具有更廣泛的適用范圍。任歲玲等[18]提出了基于特征分析的自適應干擾抑制算法，首先對接收信號的互譜密度矩陣進行特征分析，其次利用目標的方位范圍信息剔除非目標信號子空間，重構互譜密度矩陣，實現自適應干擾抑制。該算法具有更穩健的干擾抑制能力，然而需要目標方位的先驗信息，且在強干擾與目標軌跡交叉的情況下可能會對交叉后的探測結果產生誤導。此外，文獻[11]通過提取并剔除強干擾的線譜信息達到了干擾抑制的目的，但對于存在同頻干擾的情況，該方法性能受限。

針對上述問題，本文提出一種基于多線譜干擾抑制的水下對空運動聲源的線譜探測方法，以實現對同頻干擾的抑制。通過2019年5月南中國海北部海域進行的水下探測空中直升機實驗結果，驗證了所提方法的有效性。

1 基于最小方差無失真響應的線譜探測方法

1.1 最小方差無失真響應

常規波束形成因其計算簡單、穩健性高而廣泛應用于陣列信號處理中。然而，由于其波束旁瓣較高，可能會導致弱目標信號被強干擾的“能量泄露”所淹沒，抗干擾能力較差。基于此，采用最小方差無失真響應(MVDR)波束形成器[19]處理陣列數據。理論上，MVDR是多目標干擾情況下的最佳陣列處理器，具有較強的旁瓣抑制能力。

首先，給出頻率分量ωk對應的MVDR波束形成器的加權向量(k=1,2,…,K，K為頻點數)：

(1)

式中：ωk表示第k個頻率值；a(ωk,θ)表示頻率ωk、期望方向θ對應的陣列流形向量；上標H表示共軛轉置；R(ωk)是由N個快拍估計得到的頻率分量ωk的互譜密度矩陣；R-1(ωk)表示對R(ωk)的求逆運算。

對于寬帶信號，經頻率非相干處理后的寬帶MVDR空間譜可表示如下：

(2)

當進行線譜檢測時，需要利用MVDR波束形成器對每條線譜進行空間譜估計，并對所有線譜的估計結果進行綜合處理。

假設分析時間長度為L，經線譜檢測后得到I條線譜，設第i(i=1,2,…,I)條線譜在l(l=1,2,…,L)時刻共有K個頻點，則第i條線譜在l時刻的空間功率譜估計可表示為

(3)

由于不同線譜的能量有所差異，為了提高對低信噪比目標的探測能力，下面對各條線譜能量進行歸一化，即對各條線譜的空間功率譜估計結果進行歸一化。由此，l時刻綜合所有線譜貢獻的歸一化空間功率譜表示如下：

(4)

1.2 線譜提取

參照文獻[10]中的算法進行線譜提取，通過模擬聲納兵在低頻分析與記錄(LOFAR)圖中檢測和提取線譜的過程，將線譜識別、檢測和跟蹤相結合，可以在低信噪比、線譜數量未知的情況下，以較高的檢測概率實現線譜的自動提取。方法流程包括背景均衡[20]、線譜識別、線譜跟蹤、噪聲點剔除和丟失線譜產生。

下面以一次水下探測直升機實驗為例，給出上述線譜檢測算法的執行結果，實驗對應的水文環境在3.1節中給出(但與3.1節中所給實驗數據屬不同批次)。圖1(a)與圖1(b)是利用上述線譜檢測方法處理得到的線譜識別結果以及最終的線譜提取結果。圖1(b)中不同目標的線譜分別用不同顏色的圓點表示，各干擾源線譜對應其他顏色圓點。當選取第6個時間點作為線譜搜索起始時刻時，該方法在前6個時間點內未能檢測到直升機譜線。這是因為該方法在進行線譜跟蹤時，只選取搜索中心時刻上存在的線譜點作為跟蹤起始點，而未在搜索中心時刻處檢測到線譜點的譜線可能會被忽略。

圖1 文獻[10]線譜檢測各流程結果Fig.1 Detected results of line spectrum in each process in Ref.[10]

針對上述問題，為了降低線譜的漏檢概率，修改線譜跟蹤過程如下：以搜索中心時刻及其左右相鄰時刻包含的所有線譜點(去除重復元素)作為線譜搜索的起始點，進行線譜跟蹤，具體跟蹤方法以及后續剔除噪聲點、產生丟失線譜等步驟保持不變。

經上述步驟檢測得到的線譜中可能會存在兩條或兩條以上重合度較高的譜線，它們實際對應同一條線譜，此外，在檢測過程中對目標線譜頻點的識別還有可能會存在偏差。基于此，在方法執行結束后增加線譜分離、線譜加密兩步操作：

1) 線譜分離。對線譜跟蹤結果進行歸類，合并相同譜線。分為如下兩步：①計算每條線譜在搜索時間區間內的平均頻率值，將平均值之差在一定門限(可參考1 Hz)內的線譜劃分至疑似合并線譜集合。②以每個集合內的第1條線譜作為參考線譜，選取其他線譜與參考線譜同時存在線譜點時對應的線譜信息，求取兩組數據差值的標準差；若標準差小于一定門限，則與參考線譜視為同一條譜線，并將該線譜信息合并至參考線譜中，依次類推。

2) 線譜加密。對分離出的每條線譜在每個時刻分別進行加密，將線譜點左右一定頻率范圍內能量大于線譜檢測門限的頻點一并納入該線譜。

利用改進后的線譜檢測方法提取到的線譜如圖2所示。從圖2中可以看出，改進后的方法通過借助左右相鄰時刻(第5、7個時間點)處的目標線譜點信息進行搜索，可跟蹤得到完整的直升機譜線，并通過插值產生第6個時間點處的丟失線譜點。

圖2 改進后線譜提取結果Fig.2 Extracted results of improved line spectrum

利用MVDR，分別對圖1(b)、圖2中的每條線譜做波束形成處理，并綜合所有線譜貢獻繪制方位歷程圖，如圖3所示。圖3(a)、圖3(b)分別對應方法改進前后的檢測結果。與線譜提取結果相對應，圖3(a)中前6個時間點內均未檢測到目標方位歷程，而圖3(b)中第2～11個時間點內直升機目標方位歷程明顯，與真實方位軌跡較為吻合。

圖3 方位歷程圖Fig.3 Detected target bearings of two methods

此外，由于改進后方法對線譜跟蹤結果進行了線譜點加密(見圖2)，使得圖3(b)中經加密處理后的目標方位更接近真實方位歷程。因此，改進后線譜檢測方法在一定情況下可相對提高線譜檢測概率，增強線譜檢測方法的探測性能。

2 基于多線譜干擾抑制的線譜探測

2.1 算法背景及思路

線譜檢測對寬帶背景噪聲、寬帶強干擾起到了很好的抑制作用，但是，當目標為空中聲源時，近水面船因存在線譜成分而成為窄帶強干擾。圖4給出了直升機噪聲入水實驗中常見的水下接收聲場功率譜密度曲線。從圖4中可以看出，分析頻段內干擾1、干擾2對應3條線譜，而目標只對應1條線譜；此外，3個干擾源均存在能量高于目標的線譜。即使對每條線譜的能量進行歸一化，在綜合所有線譜貢獻后，仍存在干擾1、干擾2的空間譜能量大于空中目標，從而影響對目標的探測能力。因此，有必要在線譜探測的基礎上繼續進行干擾抑制。

圖4 常見接收聲場功率譜密度曲線Fig.4 Common power spectral density curve of receiving sound field

文獻[11]與文獻[18]所提出的干擾抑制算法分別在抑制同頻干擾和同向干擾的情況下性能受損。基于此，本節提出一種針對多線譜干擾的抑制方法，以有效抑制同頻干擾，且對于存在同向干擾的情況避免因干擾抑制而造成的目標方位信息損失。

首先，根據線譜提取結果，將各條線譜的順序按照頻率由高到低排列，利用MVDR波束形成器，依次對各條線譜進行空間功率譜估計。其次，對于每一時刻，當存在兩條或兩條以上線譜的空間譜最大值方位相同或相近時(認為這些線譜對同一信源起主要貢獻)，只允許接收來自第1條線譜的全部貢獻，抑制其他線譜在該信源方位附近的空間譜估計結果。最后，按照這一準則累積所有線譜貢獻，得到方位歷程圖。

下面仍然以圖4為例進行具體說明，圖4中干擾2存在3條線譜，按頻率軸從右到左依次估計每條線譜的空間譜。由于干擾2第2、3條線譜空間譜的最大值方位與第1條線譜空間譜最大值方位相同或相近，在累積線譜貢獻時，只加入第1條線譜的完整空間譜以及第2、3條線譜空間譜經抑制后的結果。因此，最終的空間譜估計中相當于只有第1條線譜對干擾2起主要貢獻。

2.2 算法具體流程

步驟1對第1個時刻，計算第1條線譜的空間譜并進行歸一化，根據空間譜最大值對應的方位信息確定向量Θ1,1.

步驟2計算第2條線譜的歸一化空間譜，若其空間譜最大值所對應的方位與Θ1,1中的元素相同，則找出空間譜中與該最大值相鄰的左右兩個極小值點(對應方位分別為θ1和θ2)，將介于θ1和θ2之間的空間譜替換為這兩點之間的插值結果。

步驟3按照步驟2的方法，依次計算其他線譜的歸一化空間譜，通過比較空間譜最大值方位與Θi-1,l中的元素，決定是否更新空間譜。

步驟4將所有線譜在第1個時刻處的歸一化空間譜相加，即為第1個時刻干擾抑制處理后的空間譜P_IS1(θ).

步驟5按照上述步驟依次計算第2，3，4，…，L時刻的空間譜。

經過上述干擾抑制，l時刻第i條線譜的歸一化空間譜可表示如下：

P_ISi,l(θ)=

(5)

綜合所有線譜貢獻，l時刻最終的歸一化空間功率譜為

(6)

在每個信源所對應的空間譜中，該方法只允許一條線譜起主要貢獻，從而抑制了多線譜干擾的能量。由于此處計算的是每條線譜的歸一化空間譜，不同線譜之間的能量差異不予考慮。因此，選取哪條線譜作為干擾源的主要貢獻者，對結果的影響不大。為了盡量提高線譜空間譜的方位分辨率，此處以同一信源的最高頻線譜作為具有主要貢獻的線譜。

為了提高方法的穩健性，可根據實際情況采用以下兩種方法對參數進行微調。

1) 單條線譜空間譜方位分辨率較差，導致對空間譜最大值方位的估計存在一定誤差，因此在判斷第i條線譜對應的目標是否為多線譜目標時，可適當放寬條件：若第i條線譜對應空間譜的最大值方位與數組Θi-1,l中某一個元素值之差的絕對值小于一定的角度容限，則將二者視為同一信源，并抑制此線譜在最大值方位處的空間譜能量。

2) 當存在兩個或兩個以上強干擾源對應的線譜頻率相同或相近時，線譜空間譜將同時出現多個較大極值。如果只抑制能量最大的干擾源，即只將最大值方位寫入數組Θi,l而忽略其他極大值方位，則其他干擾源在其他頻率處仍會貢獻線譜，干擾抑制不完全。因此，記錄入數組Θi,l中的每條線譜空間譜極大值的個數可根據實際情況在1～3之間進行選擇。

綜上所述，給出基于多線譜干擾抑制的水下對空運動聲源線譜探測流程圖，如圖5所示。

圖5 方法流程圖Fig.5 Flowchart of the proposed method

3 海試數據驗證

將本文多線譜干擾抑制方法應用于水下對空探測的海試數據中，并與文獻[11]、文獻[18]提出的兩種干擾抑制方法相比較，驗證該方法在實際工程應用中的有效性。分別選取兩組來自不同實驗的海試數據進行處理分析，兩組數據中空中聲源均以較高速度飛行，采用布放于海底的水平陣接收聲信號，接收聲場中均伴隨有較多的水面船只干擾。

3.1 存在同頻干擾，且目標與干擾的方位歷程存在交叉

3.1.1 實驗概況

2019年5月，在中國南海北部海域進行了一次海底水平陣探測空中直升機的實驗。海底水平陣由56個陣元組成。實驗期間直升機繞海底陣運動，運動軌跡如圖6所示，運動范圍始終保持在與陣相距2.5 km以內，聲源的速度與高度隨時間波動(此處直升機的位置及工況與1.2節中的實驗數據有所差別)。實驗海區海深約80 m，聲速剖面如圖7所示，水深5 m以內、5～50 m之間存在兩個溫躍層，呈現出不同的負梯度，水深50 m以下接近等聲速分布。

圖6 空中聲源相對參考陣元運動軌跡Fig.6 Trajectory of airborne sound source relative to reference array element

圖7 實驗海區聲速剖面Fig.7 Sound velocity profile of experimental sea area

3.1.2 線譜檢測

選取低頻范圍作為分析頻段，歸一化后的頻率范圍為[6×10-3，12.4×10-3]。接收信號總時長為440 s，其中1次空間譜估計的積分時間為40 s，將其分成互不重疊的10個頻域快拍，每個快拍對應的傅里葉變換點數為65 536.相鄰時段空間譜估計的滑動步長為20 s，整段信號被分為21個時間段。

圖8給出了參考陣元接收聲場時頻圖，橫坐標第1～21個時間點，分別對應1～21個分析時段的中間時刻。觀察圖8可以看出，分析頻段內存在多條線譜，不同線譜之間的能量相差較大。按照1.2節所給方法分別進行背景均衡、線譜識別與跟蹤、線譜分離，最終分離出的10條線譜如圖9所示，相鄰兩條譜線用不同顏色的圓點進行區分。圖9中黑色實線為根據聲源全球定位系統(GPS)信息計算得到的接收聲場瞬時頻率曲線，由于聲源的運動速度與運動方向都在不斷變化，因此得到的瞬時頻率曲線具有明顯的波動。

圖8 參考陣元時頻圖Fig.8 Spectrogram of reference array element

圖9 線譜分離結果Fig.9 Result of line spectrum separation

對每條線譜進行加密，并利用MVDR波束形成器進行線譜檢測，綜合所有線譜的貢獻，累積時間得到方位歷程圖，如圖10所示。由圖10可以看出，在190°和230°方向附近均存在較強的干擾，干擾的方位歷程較為清晰、連續，從而掩蓋了第7～11個時間點內的目標軌跡。此外，在第15～18個時間點處目標方位歷程出現間斷。總體而言，由于多個強干擾的存在，僅依靠線譜探測無法得到較清晰完整的目標方位歷程。

圖10 線譜探測結果Fig.10 Detection result of line spectrum

為了明確多個水面干擾源的存在對水下探測空中目標的影響，圖11給出了各個時間點處空中目標與各干擾源之間的信干比隨時間的變化曲線。選取圖10中3個強度較大的干擾源，對應方位190°、227°和236°附近，分別計算目標與3個干擾源之間的信干比。由圖11可以看出，目標與單個干擾源的信干比基本維持在0 dB以下，最低可達到-9 dB.因此，為了提高對空中目標的探測能力，需要尋求有效的干擾抑制方法。

圖11 空中目標與水面干擾的信干比估計Fig.11 Estimation of signal-to-interference ratio between airborne target and surface interferences

3.1.3 干擾抑制

分別利用強干擾線譜剔除、基于特征分析的自適應干擾抑制以及多線譜干擾抑制3種方法進行數據處理，比較3種方法的探測效果。

1) 強干擾線譜剔除：由于對230°附近2個干擾的方位分辨力較差，此處選取190°干擾作為抑制對象，借助文獻[11]中的方法對干擾進行線譜提取并予以剔除。

2) 基于特征分析的自適應干擾抑制：從圖10中大致確定目標在初始時刻的方位為105°，后續時刻的目標方位由前一時刻經干擾抑制后的方位估計結果給出。空中運動目標方位變化較快，此處設置目標的方位搜索范圍為前一時刻方位估計結果的左右15°區間內。

3) 多線譜干擾抑制：基于每個信源所在方位的空間譜主要接收一條線譜貢獻的準則，依次更新每條線譜的歸一化空間譜，并綜合所有線譜貢獻得到最終方位歷程圖。其中，角度容限取5°，記錄入數組Θi,l中每條線譜空間譜極大值的個數(取為2)。

圖12(a)、圖12(b)、圖12(c)分別給出了3種方法處理得到的方位歷程圖。從圖12(a)中可以看出，190°方位附近的干擾已被抑制，目標方位歷程在第1～8個時間點內清晰可見，但在第15～18個時間點處仍未顯現。圖13為190°干擾的低頻率分析與記錄(LOFAR)譜圖，由圖13可以發現，在第15～18個時間點內，190°干擾線譜與目標線譜的頻率十分接近，因此在對190°干擾進行線譜剔除的同時損失了目標線譜信息，導致目標在該段時間內方位缺失。

圖12 干擾抑制結果Fig.12 Results of interference suppression

圖13 190°干擾LOFAR譜圖Fig.13 LOFAR spectrogram of 190° interference

觀察圖12(b)發現，在第1～6個時間點內，由于算法對非目標信號子空間進行了抑制，使得190°和230°附近的干擾并未出現在方位歷程圖中。然而在第7個時間點處，目標與190°干擾發生軌跡交叉，由于干擾能量遠大于目標，下一時刻得到的方位估計值為干擾所對應的方位，以至于自第8個時刻起，對目標的方位估計發生偏差，轉變為對190°干擾的方位跟蹤。此外，由于所給的目標方位搜索范圍較大，對非干擾子空間的選取不夠精確，導致方位歷程圖中的峰值旁瓣比較低，背景相對模糊。

相比圖12(a)與圖12(b)，圖12(c)中的結果既顯示出較完整的目標、干擾運動軌跡，又保持了較高的峰值旁瓣比，便于同時對目標與干擾進行分析。圖12對應的海試數據處理結果表明，本文所提方法對于解決存在同頻干擾或目標與強干擾軌跡交叉的線譜探測問題具有明顯的優勢。

3.2 不存在同頻干擾

為了證明本文所提方法同樣適用于抑制非同頻干擾，下面用該方法處理文獻[11]中直升機探測實驗的海試數據。本次實驗與3.1節中所述實驗毫無關聯，具體實驗細節可參考文獻[11]。圖14(a)是利用所提方法處理得到的方位歷程結果，為了進行對比，在圖14(b)中給出了文獻[11]中的處理結果。對比發現，圖14(a)中目標的方位歷程在前11個時間點內更為連續，能量更加集中；在第12～13個時間點內，由于150°干擾與目標的運動軌跡存在交叉，且干擾的能量相對較強，使得目標的方位歷程在后續時間段內銜接性較差。總體而言，利用所提方法可有效處理文獻[11]中的實驗數據，實現非同頻干擾下水下對空運動聲源的線譜探測。

圖14 目標探測結果對比Fig.14 Results of target detection

4 結論

為了解決水下對空中高速運動聲源的探測問題，本文在線譜探測的基礎上提出了一種對多線譜干擾的抑制方法。針對分析頻段內干擾源線譜數量遠大于目標線譜且干擾線譜能量較強的特征，限定每個信源的空間譜主要接收來自一條能量歸一化的線譜的貢獻，最后綜合所有貢獻線譜的空間譜估計得到最終的方位估計結果。由于分析頻段內目標與干擾的空間譜各自只對應一條線譜的貢獻，且對線譜的能量進行了歸一化，因此目標的方位歷程得到增強。本文方法可同時有效抗多個干擾源，且在干擾與目標同頻的情況下依舊保持較好的性能。海試實驗結果表明，在多個近水面船干擾情況下，基于多線譜干擾抑制的線譜探測方法可成功實現水下對空中運動聲源的探測，相比其他方法具有更廣泛的適用范圍。