李海鵬 孫大軍 鄭翠娥

(哈爾濱工程大學(xué)水聲技術(shù)重點(diǎn)實(shí)驗(yàn)室 哈爾濱 150001)

(海洋信息獲取與安全工信部重點(diǎn)實(shí)驗(yàn)室(哈爾濱工程大學(xué)) 哈爾濱 150001)

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

1 引言

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

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

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

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

2 陣列信號(hào)模型與子空間理論

2.1 子空間框架下的陣列信號(hào)模型

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

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

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

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

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

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

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

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

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

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

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

2.2 線性子空間投影

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 基于最大似然估計(jì)的信號(hào)子空間估計(jì)

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

3.1 子空間維度估計(jì)

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

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

其中， Tr(·)表 示矩陣的跡。因此，通過(guò)遍歷r =1,2,···,N并使BIC最小即可獲得信號(hào)子空間的維度。同理，在水聲定位過(guò)程中，通?？梢栽诿看谓邮掌谕盘?hào)前，采集一段僅包含干擾和噪聲的測(cè)試信號(hào)xc(t)，同樣可以通過(guò)上述方法估計(jì)干擾信號(hào)子空間維 度。

3.2 期望信號(hào)協(xié)方差矩陣估計(jì)

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

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

4 仿真驗(yàn)證

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

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

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

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

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

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

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

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

圖4 基陣1號(hào)基元接收信號(hào)

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

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

圖7 時(shí)延估計(jì)誤差隨SINR變化圖，INR=0 dB

圖8 時(shí)延估計(jì)誤差隨SINR變化圖，INR=20 dB

圖9 時(shí)延估計(jì)誤差隨SINR變化圖，INR=50 dB

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

5 結(jié)論

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