李海鵬 孫大軍 鄭翠娥

(哈爾濱工程大學(xué)水聲技術(shù)重點實驗室 哈爾濱 150001)

(海洋信息獲取與安全工信部重點實驗室(哈爾濱工程大學(xué)) 哈爾濱 150001)

(哈爾濱工程大學(xué)水聲工程學(xué)院 哈爾濱 150001)

1 引言

海洋作為地球最大的生態(tài)系統(tǒng)影響著全球能量流動、氣候變化與生態(tài)安全，將地球連結(jié)為一個命運共同體。我國是一個海洋大國，海洋面積遼闊，認識海洋、經(jīng)略海洋、建設(shè)海洋強國具有重要的戰(zhàn)略地位。深海面積超過海洋總面積的90%，走向深海是海洋強國的必由之路。對深海環(huán)境特性的精確認知和對深海資源的科學(xué)開發(fā)利用是建設(shè)海洋強國的基礎(chǔ)[1–3]。隨著對海洋特別是深海探索和開發(fā)的深入，對各類水下潛器、平臺的高精度定位導(dǎo)航需求越來越強烈。水聲定位系統(tǒng)是現(xiàn)代深海作業(yè)必備的高精度水下定位裝備，針對復(fù)雜多變的海洋環(huán)境開展深海高精度水下聲學(xué)定位技術(shù)研究，將成為推動海洋強國建設(shè)不斷取得新成就的必要手段[3]。

基于多元傳感器陣列的水聲定位系統(tǒng)通過水面聲學(xué)基陣與水下聲學(xué)應(yīng)答器進行聲波交互，通過估計聲波從應(yīng)答器到達基陣中心的傳播時延和聲波在水中的傳播速度估計基陣與目標的距離，通過聲波到達各基元的時延差(或相位差)估計目標方位，從而獲得目標相對基陣中心的位置，再結(jié)合羅經(jīng)、全球定位系統(tǒng)(Global Positioning System, GPS)等外接輔助設(shè)備轉(zhuǎn)換得到目標的絕對位置。目前，高精度水聲定位面臨著定位信號長距離傳播導(dǎo)致能量衰減，水面作業(yè)船、海洋生物以及海洋環(huán)境等因素產(chǎn)生的干擾(噪聲)，導(dǎo)致接收數(shù)據(jù)的信干噪比(Signal to Interference and Noise Radio, SINR)較低，而SINR是決定時延估計精度的重要因素，從而制約定位精度進一步提高，因此抑制干擾(噪聲)的影響是高精水聲定位中不可避免且亟待解決的問題之一。

目前的干擾抑制算法主要分為兩類，一類是針對單通道的降噪方法，如最小均方誤差(Least Mean Square, LMS)自適應(yīng)濾波器法[4–6]、短時傅里葉變換(Short Time Fourier Transform, STFT)[7,8]等。LMS 算法基于最小均方誤差準則，通過輸入量與期望響應(yīng)的差值對權(quán)值進行迭代更新，以獲取最優(yōu)權(quán)值。LMS 算法具有計算量小、穩(wěn)健性強、易于實現(xiàn)等優(yōu)點，但該方法的收斂過程慢，而且對于隨機干擾的適應(yīng)性較差，而在復(fù)雜的海洋環(huán)境中干擾的統(tǒng)計特性往往是復(fù)雜且隨機的?；诙虝r傅里葉變換的干擾抑制算法通過對接收數(shù)據(jù)進行短時傅里葉變換，根據(jù)期望信號和干擾在時頻域的能量分布對期望信號進行重構(gòu)，從而達到抑制干擾的效果，但該方法對于SINR要求較高。與此同時，LMS和STFT算法的共同問題是會影響期望信號的相位，這會嚴重影響定位系統(tǒng)的定位精度。另一類算法是針對陣列信號的算法，包括波束形成類方法和子空間類算法。波束形成類算法[9–12]通過對各基元的接收數(shù)據(jù)進行加權(quán)，從而在期望方向形成波束，抑制非期望方向的信號，可以視為一種空域濾波器。但是在進行波束設(shè)計時往往需要已知陣列流型，從而針對性地設(shè)計波束，而在水聲定位作業(yè)中目標的方位通常是未知的。子空間類算法是依靠對數(shù)據(jù)矩陣的奇異值分解或?qū)f(xié)方差矩陣的特征值分解估計信號子空間和噪聲子空間，子空間類算法有3個重要理論基礎(chǔ)：(1)信號子空間與噪聲子空間垂直；(2)對于窄帶信號模型，陣列流型所張成的子空間與信號子空間相等；(3)對于窄帶信號模型，信號子空間維度等于信源個數(shù)。子空間類算法的典型應(yīng)用為多信號分類(MUltiple SIgnal Classification, MUSIC)算法，利用信號子空間與噪聲子空間垂直特性進行頻率估計、方位估計等?；谧涌臻g的特性，還有學(xué)者提出了基于子空間理論的干擾抑制方法。Bose等人[13]提出了通過子空間類算法進行語音降噪，通過對信號協(xié)方差矩陣進行特征值分解估計信號子空間，并利用信號子空間與噪聲子空間的垂直特性將帶噪數(shù)據(jù)線性投影到信號子空間中以實現(xiàn)數(shù)據(jù)降噪，該算法只考慮了信號中僅包含高斯噪聲的情況，當(dāng)干擾存在時，信號子空間與干擾子空間將發(fā)生空間糾纏，導(dǎo)致無法抑制干擾。針對干擾和噪聲同時存在的情況，目前的研究內(nèi)容主要集中在對窄帶信號的抑制。張春海等人[14]提出基于子空間跟蹤的直接序列擴頻 (Direct-Sequence Spread-Spectrum, DSSS)通信系統(tǒng)抗窄帶干擾研究，通過跟蹤接收信號自相關(guān)矩陣大特征值對應(yīng)特征矢量構(gòu)成的干擾子空間，實現(xiàn)對窄帶干擾的有效抑制。周峰等人[15]提出了一種用于合成孔徑雷達的基于回波數(shù)據(jù)特征子空間濾波的干擾抑制方法，首先在頻域?qū)φ瓗Ц蓴_進行識別，然后在時域?qū)Ω蓴_進行抑制處理。張小飛等人[16]提出一種基于斜投影的波束形成算法，算法通過構(gòu)造斜投影矩陣先對接收數(shù)據(jù)進行斜投影抑制干擾和噪聲的影響，然后進行波束設(shè)計，進而提高了波束形成的魯棒性，但算法只適用于窄帶模型，且要求陣列流型已知，但在水聲定位過程中，信號和干擾均為寬帶且目標的方位是未知的。

綜上，基于子空間的干擾抑制算法目前還面臨如下挑戰(zhàn)：(1)接收數(shù)據(jù)中同時包含期望信號、干擾和噪聲；(2)期望信號和干擾均為寬帶；(3)期望信號入射方位未知。基于上述挑戰(zhàn)，本文提出一種基于子空間理論的寬帶強干擾抑制方法，首先通過貝葉斯信息量準則(Bayesian Information Criterion,BIC)估計信號子空間和干擾子空間的維度，然后推導(dǎo)不同信號假設(shè)下的概率密度函數(shù)，求解未知參數(shù)的最大似然估計，構(gòu)造廣義似然比并通過最優(yōu)匹配廣義似然比檢測法估計與期望信號最匹配的子空間，然后以此構(gòu)造空間投影算子對接收數(shù)據(jù)進行線性投影，實現(xiàn)對干擾和噪聲的抑制。本方法的優(yōu)點在于，適用于寬帶陣列信號，且無需已知陣列流型和信源個數(shù)，同時不影響期望信號的相位信息。仿真結(jié)果表明本文所提方法能夠在低SINR條件下有效抑制干擾的影響，提高水聲定位系統(tǒng)時延估計精度。

2 陣列信號模型與子空間理論

2.1 子空間框架下的陣列信號模型

假設(shè)定位系統(tǒng)接收基陣由 N 個無指向性的基元構(gòu)成，每個基元進行時域均勻采樣接收 Lx個快拍(即采樣點個數(shù))，位于遠場的目標發(fā)射快拍數(shù)為Ls的波形已 知 的 定 位 信 號s0(t) , t ∈?s≡{K+1,K+2,···,K+Ls}且 K+Ls

即在 H0假設(shè)下，接收信號只包含干擾信號和噪聲；在 H1假設(shè)下，接收信號包含期望信號、干擾信號和噪聲。其中，期望信號可以表示為

其中，[ ·]T表示轉(zhuǎn)置運算。τk表示期望信號到達第k個基元和參考點之間的時延差，該時延差由入射信號方位和基元位置關(guān)系決定。

因此，包含 Lx個快拍的接收數(shù)據(jù)矩陣X ∈?N×Lx可以表示為

在子空間理論框架下，假設(shè)期望信號和干擾分別位于兩個獨立的子空間< Q >和< F>中。因此，期望信號可以表示為一個列滿秩矩陣 Q ∈?N×r所有列的線性組合[17,18]，r為期望信號子空間的維度。于是s (t)可以表示為

其中， A(t)∈?r×1為 期望信號在子空間中的坐標。

同理，干擾信號可以表示為一個列滿秩矩陣F ∈?N×q所有列的線性組合，q為干擾子空間的維度，于是i (t)可以表示為

其中，B (t)∈?q×1為 干擾信號在子空間< F>中的坐標。

將式(4)和式(5)代入式(1)可得子空間理論框架下基陣在t 時刻的接收數(shù)據(jù)

進一步假設(shè)干擾為隨機干擾且協(xié)方差為Ri，即

其 中，σ2I為噪聲的方差。

2.2 線性子空間投影

根據(jù)2.1節(jié)所述的接收信號模型，假設(shè)期望信號、干擾信號和噪聲相互獨立，則接收數(shù)據(jù)的協(xié)方差矩陣可表示為

其中，RA=EA(t)AT(t)。

令子空間由列滿秩矩陣[ Q F]的列張成，即

其中， span{M} 表示由矩陣 M的列張成子空間。因此

其中，對角陣 P=diag{RARB}表示復(fù)合向量[A(t) B(t)] 的協(xié)方差矩陣。因此Rx可以簡化為

由于 M為列滿秩矩陣且期望信號和干擾信號是相互獨立的，所以矩陣 P為非奇異的埃爾米特矩陣，因此可以推斷 Rx一定是非負正定埃爾米特矩陣[19]。根據(jù)埃爾米特矩陣的特性，可以對 Rx進行特征值分解，則有

其 中， W=diag(λ1λ2··· λN)為 特 征 值 矩 陣，U=[u1u2··· uN]為 特征向量矩陣，當(dāng)r+q

當(dāng)i r+q時，{λiui} 對應(yīng)噪聲的特征值和特征向量。因此，Rx的特征分解可以進一步表示為

需要注意的是， UM和 M的列分別為子空間< M>的一組基，在數(shù)值上并不一定相等，這就意味著無法直接從UM中將Q 和F 分離。

如果能正確估計期望信號和干擾信號的子空間，可通過線性投影的方式實現(xiàn)抑制干擾和噪聲的影響。一般地，子空間和< F>的正交投影算子可以表示為

綜上，若已知張成信號子空間和干擾子空間的矩陣 Q和 F，可以構(gòu)造對應(yīng)的正交投影算子或斜投影算子，通過線性投影運算實現(xiàn)干擾和噪聲抑制，主要有兩種實現(xiàn)方式：

圖1 子空間及其投影關(guān)系示意圖

方式1 將接收數(shù)據(jù)投影在信號子空間上從而在保留期望信號的前提下，抑制部分干擾和噪聲的影響

方式2 將接收數(shù)據(jù)投影在斜投影空間上從而在保留期望信號的前提下，抑制全部干擾和部分噪聲的影響

為了便于討論，下文分別將兩種投影方式命名為 SP ?I和 SP ?II。

3 基于最大似然估計的信號子空間估計

根據(jù)式(13)—式(15)可知，決定子空間的有兩個參數(shù)，即子空間維度和接收數(shù)據(jù)協(xié)方差矩陣。當(dāng)干擾存在時，信號子空間和干擾子空間會出現(xiàn)糾纏現(xiàn)象，通過對接收數(shù)據(jù)協(xié)方差矩陣進行特征值分解僅能得到信號和干擾的聯(lián)合子空間，為了抑制干擾的影響則需要單獨估計信號子空間，本節(jié)的目的是分別估計子空間維度和信號協(xié)方差矩陣，然后構(gòu)造投影算子抑制干擾和噪聲。

3.1 子空間維度估計

根據(jù)文獻[21]可知，入射信號確定時，信號子空間的維度與入射角度無關(guān)。因此對于定位信號已知的水聲定位系統(tǒng)而言，可以假定任意方位的入射信號，然后直接估計信號子空間的維度。假設(shè)長度為Ls接收數(shù)據(jù)xs(t),t ∈{1,2,···,Ls}中只包含期望信號和高斯白噪聲，即

其中，‖ ·‖表示矩陣(向量)的二范數(shù)。對式(29)取對數(shù)，則

其中， Tr(·)表 示矩陣的跡。因此，通過遍歷r =1,2,···,N并使BIC最小即可獲得信號子空間的維度。同理，在水聲定位過程中，通常可以在每次接收期望信號前，采集一段僅包含干擾和噪聲的測試信號xc(t)，同樣可以通過上述方法估計干擾信號子空間維 度。

3.2 期望信號協(xié)方差矩陣估計

廣義似然比作為最大的不變統(tǒng)計量，廣泛用于信號檢測、方位估計等領(lǐng)域[23]。本節(jié)推導(dǎo)了不同信號假設(shè)下的概率密度函數(shù)，求解出未知參數(shù)的最大似然估計，構(gòu)造廣義似然比并用匹配廣義似然比檢測法估計與期望信號最匹配的子空間。廣義似然比定義為有約束條件下的概率密度函數(shù)最大值與無約束條件下概率密度函數(shù)最大值之比，首先假設(shè)信號和干擾的子空間均已知，根據(jù)式(6)所定義的信號模型，構(gòu)造廣義似然比檢測的表達式為

根據(jù)第2節(jié)所述內(nèi)容可知，在信號波形確定的情況下，信號子空間只與期望信號的入射方位有關(guān)，因此可以將Q 寫成Qα,β，其中α 和 β 分別表示信號入射的俯仰角和方位角。期望信號的協(xié)方差可以表示為

4 仿真驗證

本節(jié)通過仿真分析所提方法的性能并與現(xiàn)有方法進行對比。如圖2所示，仿真采用一個均勻分布的平面陣列，基元個數(shù)為30個，相鄰陣元間距為5 cm。每個定位周期的采樣信號長度為25 ms，采樣頻率為100 kHz。定位信號為線性調(diào)頻信號，頻帶寬度為9～kHz，信號長度為5 ms。信號入射的俯仰角和方位角分別為75°和30°。干擾為隨機干擾，入射的俯仰角和方位角分別為10°和60°。噪聲為高斯白噪聲。

互相關(guān)時延估計法是水聲定位中最常用且最有效的時延估計算法，然而互相關(guān)時延估計法的精度受到干擾和噪聲的影響，本節(jié)首先通過仿真SINR對互相關(guān)時延估計精度的影響。圖3展示了不同的干擾噪聲比(Interference to Noise Radio, INR)條件下互相關(guān)時延估計精度隨SINR的變化。

根據(jù)圖3可以發(fā)現(xiàn)，互相關(guān)時延估計精度會隨著SINR的減小而降低；在相同的SINR條件下，INR越小時延估計誤差越大，原因是干擾與期望信號是非相干的，當(dāng)干擾為主要分量時，互相關(guān)時延估計精度要優(yōu)于噪聲為主要分量的情況。接下來驗證本文所提出的方法對于時延估計精度的提升效果。按上述仿真條件，得到SINR為0 dB，INR為0 dB的基陣接收數(shù)據(jù)，圖4展示了1號基元的接收數(shù)據(jù)。

根據(jù)圖4所示的接收數(shù)據(jù)估計期望信號子空間，圖5展示了在角度空間內(nèi)的搜索結(jié)果，其中顏色對應(yīng)廣義似然比的大小。為了便于分析，分別提取峰值位置對應(yīng)的俯仰角和方位角方向的切面，如圖6所示，結(jié)果表明本文所提方法能夠精確估計信號子空間。

接下來，根據(jù)式(47)的估計結(jié)果，分別構(gòu)造正交投影算子和斜投影算子，將接收數(shù)據(jù)進行對應(yīng)的線性投影，然后通過互相關(guān)法估計各基元接收信號的時延，并與LMS和STFT算法的結(jié)果進行對比。當(dāng)INR=0 dB時，經(jīng)過1000次蒙特卡羅仿真，各種方法輸出數(shù)據(jù)的時延估計均方根誤差(Root Mean Squared Error, RMSE)隨SINR的變化如圖7所示。

為了進一步驗證算法在不同INR情況下的性能，分別對INR=20 dB和INR=50 dB兩種情況進行仿真，經(jīng)過1000次蒙特卡羅仿真，上述兩種情況的時延估計的RMSE隨SINR的變化分別如圖8和圖9所示。

圖2 定位系統(tǒng)接收基陣陣型圖

圖3 相關(guān)時延估計隨SINR的變化

圖4 基陣1號基元接收信號

圖5 廣義似然比在角度空間內(nèi)的搜索結(jié)果

圖6 廣義似然比隨入射角度變化圖

圖7 時延估計誤差隨SINR變化圖，INR=0 dB

圖8 時延估計誤差隨SINR變化圖，INR=20 dB

圖9 時延估計誤差隨SINR變化圖，INR=50 dB

根據(jù)圖7—圖9的仿真結(jié)果可以發(fā)現(xiàn)，當(dāng)INR較大即非期望信號中干擾信號為主要分量時，基于斜投影的降噪方法( SP ?II)的性能明顯優(yōu)于其他方法；當(dāng)INR較小即非期望信號中高斯噪聲為主要分量時，基于正交投影的降噪方法(S P ?I)明顯優(yōu)于其他方法；隨著INR的降低， SP ?I 和S P ?II的性能均會下降，但是 SP ?I對于噪聲的敏感程度低于SP ?II。上述現(xiàn)象產(chǎn)生的原因是正交投影和斜投影的投影方式導(dǎo)致的，斜投影沿著干擾子空間的方向?qū)⑵谕盘柾队暗叫盘栕涌臻g，因此能夠最大限度地消除干擾，而正交投影僅將期望信號進行正交投影，當(dāng)干擾子空間與信號子空間非正交時，會有部分干擾信號分量投影到信號子空間中。當(dāng)噪聲分量較高時，會對干擾子空間的估計產(chǎn)生影響，進而影響斜投影的性能，而正交投影不受干擾子空間的影響。

5 結(jié)論

本文針對強干擾降低水聲定位系統(tǒng)時延估計精度的問題，提出一種基于子空間理論的寬帶強干擾抑制方法，通過估計期望信號子空間和干擾子空間，構(gòu)造投影算子并對接收數(shù)據(jù)進行線性投影，從而抑制干擾和噪聲對定位系統(tǒng)時延估計精度的影響。相比傳統(tǒng)方法，本文所提方法可適用于寬帶陣列信號，且無需已知陣列流型和信源個數(shù)，同時不影響期望信號的相位信息。仿真結(jié)果顯示，本文所述的方法能有效抑制寬帶強干擾的影響，提高系統(tǒng)時延估計精度。時延估計誤差的仿真結(jié)果顯示，當(dāng)非期望信號中干擾信號為主要分量時，基于斜投影的降噪方法性能最優(yōu)；當(dāng)非期望信號中高斯噪聲為主要分量時，基于正交投影的降噪方法最優(yōu)。