基于壓縮傳感的水下目標(biāo)方位估計(jì)方法*

（杭州應(yīng)用聲學(xué)研究所聲納技術(shù)重點(diǎn)實(shí)驗(yàn)室 杭州 310023）

1 引言

聲信號(hào)的波達(dá)方向（direction-of-arrival，DOA）估計(jì)是聲納在目標(biāo)探測(cè)過程中非常重要的一個(gè)環(huán)節(jié)，它可以估計(jì)水下目標(biāo)的空間方位。常規(guī)波束形成（conventional beamforming，CBF）［1］是 DOA 估計(jì)的傳統(tǒng)方法，其分辨力受到瑞利限的限制。為了突破這一限制，發(fā)展起了一些高分辨的空間譜估計(jì)方法，如基于子空間分解的多重信號(hào)分類（multiple signal classification，MUSIC）方法［2］。除此之外，為了自適應(yīng)地抑制干擾，基于最大信噪比準(zhǔn)則的最小方差無失真響應(yīng)（Minimum Variance Distortionless Response，MVDR）方法［3］可以自適應(yīng)地估計(jì)DOA。無論是MUSIC方法還是MVDR方法，都需要計(jì)算信號(hào)的協(xié)方差矩陣，這導(dǎo)致了低快拍數(shù)、相干源、低信噪比條件下的性能下降。

2006年前后，壓縮傳感（compressive sensing or compressed sensing）又稱為壓縮采樣（compressive sampling）發(fā)展起來。E.Candès和 T.Tao等［4～6］證明了，利用信號(hào)的稀疏性，能夠用比Nyquist-Shannon采樣定理所要求的少得多的采樣實(shí)現(xiàn)稀疏信號(hào)的完全重構(gòu)，因此取名壓縮采樣。同時(shí)期，D.Donoho等［7］提出了一種新的信息獲取指導(dǎo)理論，即壓縮傳感。CS用于求解欠定線性逆問題，致力于發(fā)展信號(hào)的稀疏性或可壓縮性，用比Nyquist采樣率更低的采樣率做采樣而實(shí)現(xiàn)信號(hào)的重構(gòu)。稀疏性或可壓縮性是信號(hào)的自然性質(zhì)，聲信號(hào)的稀疏性或可壓縮性表現(xiàn)在時(shí)域、頻率域、空間域、波數(shù)域、空間波束域等域中。傳統(tǒng)DOA估計(jì)方法的目標(biāo)是使重構(gòu)信號(hào)的能量最小，因而得到的是低分辨力、非稀疏的解。將CS應(yīng)用到聲信號(hào)的DOA估計(jì)中正是利用了空間波束域的稀疏性，因此得到的是高分辨的稀疏解。李璇等［8］分別利用對(duì)角加載最小二乘、λ1正則化和正交匹配追蹤方法求解單快拍情況下的DOA估計(jì)問題，得到了高分辨的方位估計(jì)結(jié)果。房云飛等［9］提出了一種基于波束域的多重測(cè)量向量欠定系統(tǒng)正則化聚焦求解DOA的算法，結(jié)果表明有較好的估計(jì)精度和角度分辨力。A.Xenaki和P.Gerstoft等［10～11］把單快拍和多快拍情況下的 DOA估計(jì)分別建模為單一測(cè)量向量和多重測(cè)量向量的問題，結(jié)果表明CS方法的性能優(yōu)于傳統(tǒng)的波束形成和高分辨力方法。

本文將針對(duì)陣元間距大于半波長(zhǎng)的線陣，即空間采樣率低于Nyquist采樣率進(jìn)行研究。首先把DOA估計(jì)建模為一個(gè)標(biāo)準(zhǔn)的CS問題，并分析傳感矩陣的相干性（coherence），然后利用凸優(yōu)化方法［12～13］和基于閾值的方法［14］對(duì)DOA估計(jì)這一欠定線性逆問題進(jìn)行求解。數(shù)值仿真結(jié)果和實(shí)驗(yàn)結(jié)果均表明了基于CS的波束形成方法的可行性，其分辨力要高于傳統(tǒng)的DOA估計(jì)方法，可在單快拍、相干源、低信噪比條件下估計(jì)目標(biāo)的空間方位。

2 壓縮傳感理論

經(jīng)典的Shannon采樣定理以帶寬為先驗(yàn)信息，是一個(gè)充分非必要條件，對(duì)于頻率低于最高頻率的分量，一個(gè)周期上采的樣本多于兩個(gè)，因此存在冗余。空間采樣亦是如此，陣元間距為最高頻率的半波長(zhǎng)時(shí)，那么小于最高頻率的波長(zhǎng)，一個(gè)周期上有兩個(gè)以上的樣本，存在冗余。對(duì)于傳統(tǒng)的信號(hào)處理，在時(shí)域按Nyquist-Shannon采樣較容易實(shí)現(xiàn)，但由于實(shí)際中信號(hào)帶寬較大或最高頻率較高使得采樣頻率也較高，因而采集到的數(shù)據(jù)量較大；在空間域上按Nyquist-Shannon采樣需要在空間內(nèi)均勻地布放大量傳感器，成本太高。壓縮傳感從根本上顛覆了Nyquist-Shannon采樣定理，在信號(hào)滿足稀疏性或可壓縮性的條件下，直接對(duì)信號(hào)進(jìn)行壓縮采樣，通過復(fù)雜度較高的重構(gòu)算法來重構(gòu)原始信號(hào)。

壓縮傳感有三個(gè)基本要素，分別是稀疏性和可壓縮性、壓縮采樣和重構(gòu)算法。

1）稀疏性和可壓縮性

一個(gè)信號(hào)如果在某一域中的分量大多數(shù)為零，則稱這個(gè)信號(hào)是稀疏的。數(shù)學(xué)上，用λ0范數(shù)表示非零元素的個(gè)數(shù)，如果向量x∈CN最多有s個(gè)元素是非零的，就稱 x是s-稀疏的，表示為‖x‖0≤s。稀疏性是一個(gè)強(qiáng)約束，實(shí)際中的許多信號(hào)是可壓縮的。可壓縮性指信號(hào)可以很好地由稀疏信號(hào)近似表示，這種表示常常是在經(jīng)過適當(dāng)?shù)幕儞Q之后，如Fourier變換、小波變換等。假設(shè)信號(hào)z∈CN可由稀疏向量x∈CN表示為

其中Φ∈CN×N稱為變換基矩陣。

2）壓縮采樣

利用測(cè)量矩陣Ψ∈CM×N對(duì)信號(hào)z進(jìn)行采樣，得到測(cè)量數(shù)據(jù)y∈CM為

其中 A∈CM×N稱為傳感矩陣，是變換基矩陣和測(cè)量矩陣的結(jié)合。測(cè)量矩陣Ψ可以是隨機(jī)的，也可以是確定的，其行數(shù)M一定小于列數(shù)N，故該采樣過程把一個(gè)高維信號(hào)變換為低維數(shù)據(jù)，因而稱為壓縮采樣。

3）重構(gòu)算法

經(jīng)過壓縮采樣后得到了一個(gè)欠定線性方程

壓縮傳感的目的是利用稀疏性，從少量的測(cè)量數(shù)據(jù)y中重構(gòu)x或z，本質(zhì)是求解該欠定線性方程的稀疏解。重構(gòu)算法能夠在適當(dāng)?shù)募僭O(shè)下，求得上述欠定線性方程的稀疏解。要想保證重構(gòu)算法能準(zhǔn)確地和有效地重構(gòu)信號(hào)，傳感矩陣A需要滿足一定的性質(zhì)，如零空間特性（null space property，NSP）、約束等距特性（restricted isometry property，RIP）和相干性等，其中前兩者在計(jì)算上都是NP-hard的，因此論文中考慮相干性。重構(gòu)算法可以粗略地分成三類：優(yōu)化方法，貪婪方法和基于閾值的方法。

3 基于壓縮傳感的DOA估計(jì)方法

3.1 線陣DOA估計(jì)的CS模型

假設(shè)聲源位于線陣的遠(yuǎn)場(chǎng)，即線陣接收到的聲波可以考慮為平面波，考慮聲源輻射的聲信號(hào)是窄帶的。線陣示意圖如圖1所示。

根據(jù)圖1示意圖，目標(biāo)聲源的DOA由相對(duì)于陣列方向的 θ∈[-90°，90°]決定。第i個(gè)聲源對(duì)每個(gè)傳感器的傳播時(shí)延可以由駕駛向量表示：

圖1 線陣示意圖

令向量x∈CN由感興趣方向θ格點(diǎn)上的源信號(hào)幅度組成，y表示聲場(chǎng)在M個(gè)傳感器處的測(cè)量向量。由于我們感興趣的是角度格點(diǎn)上的高分辨力，因此有M＜N。傳感矩陣則由所有可能方向的駕駛向量組成，即

它是測(cè)量矩陣Ψ和基變換矩陣Φ的乘積，其中Ψ表示聲場(chǎng)在傳感器位置處的空間采樣，是一個(gè)確定的矩陣，Φ表示一個(gè)IDFT變換基。考慮存在噪聲n∈CM的情況，測(cè)量向量就可以表示為

這是一個(gè)標(biāo)準(zhǔn)的欠定線性系統(tǒng)，我們的目標(biāo)就是利用x的稀疏性，從y中重構(gòu)x。

3.2 優(yōu)化方法和基于閾值的方法

CS中求解式（1）的問題被稱為（P0）問題，它屬于非凸優(yōu)化問題，且一般是NP-hard的，因此通常把（P0）問題凸松弛為（P1）問題，表示為

（P1）問題是一個(gè)凸優(yōu)化問題，常被稱作λ1最小化或基追蹤（basis pursuit，BP）。對(duì)于存在噪聲的情況，可以利用Lagrange乘子法求解式（4），被稱為問題，表示為

其中正則化參數(shù)λ＞0控制著稀疏解和測(cè)量值擬合之間的相對(duì)重要性。求解問題的方法被稱作基追蹤去噪（basis pursuit denoising，BPD）。在Matlab上可以利用CVX工具箱求解上述凸優(yōu)化問題，它利用內(nèi)點(diǎn)求解器得到優(yōu)化問題的全局解。

除了優(yōu)化方法，還有一種基于閾值的方法。由于λ0最小化是非線性操作，因此基于閾值的方法需要定義一個(gè)非線性算子Hs，它稱為稀疏性為s的硬閾值算子。對(duì)于一個(gè)向量x，Hs（x）保持x的s個(gè)絕對(duì)值最大的元素不變，其他元素置為零。這里介紹一種已知解是s-稀疏情況下的迭代算法，稱為迭代硬閾值（iterative hard thresholding，IHT）。該算法是通過迭代地求解定點(diǎn)方程xn+1=(I-AHA)xn+AHy得到的，其中I為單位陣，這里的上標(biāo)“H”表示矩陣的共軛轉(zhuǎn)置。由于已知所求的解是s-稀疏的，所以只需在每次迭代時(shí)保留s個(gè)絕對(duì)值最大的元素，最終可以得到一個(gè)s-稀疏的向量。

3.3 傳感矩陣的相干性分析

傳感矩陣A的特性會(huì)影響重構(gòu)算法的性能。如果不考慮噪聲，NSP是穩(wěn)健的信號(hào)重構(gòu)的充要條件。如果考慮存在噪聲，Candès和Tao引入了RIP［15］，它對(duì)于各種算法成功地從有噪測(cè)量中重構(gòu)一個(gè)稀疏信號(hào)是充分非必要的。無論是RIP還是NSP，它們?cè)谟?jì)算上都是NP-hard的，因此這里考慮易于計(jì)算的矩陣A的相干性［16］。相干性度量了傳感矩陣任意兩列之間的相關(guān)性（correlation），首先給出一個(gè)矩陣A的相干性μ(A)的定義式為

對(duì)于式（3）所表示的傳感矩陣，有

在DOA估計(jì)中，為了實(shí)現(xiàn)高分辨力，通常需要精細(xì)的格點(diǎn)，因此欠定方程中的傳感矩陣A具有線性相關(guān)的列，其相關(guān)程度反映在G的非對(duì)角元素上。下面給出了陣元間距大于半波長(zhǎng)的線陣的Gram矩陣示意圖。

從圖2中可以看出當(dāng)d=4λ5時(shí)，不出現(xiàn)柵瓣的范圍大約是[-15°，15°]。但是由于CS重構(gòu)算法是非線性的，且稀疏解具有唯一性，因而在一定程度上可以準(zhǔn)確重構(gòu)出稀疏信號(hào)。

圖2 十六元均勻線陣的Gram矩陣示意圖，d/λ=4/5

4 數(shù)值仿真和實(shí)驗(yàn)結(jié)果

4.1 數(shù)值仿真結(jié)果

考慮一個(gè)16元均勻線陣，陣元間距d=4λ/5。首先考慮遠(yuǎn)場(chǎng)存在一個(gè)目標(biāo)，入射方向?yàn)?0°，信噪比為10dB，分別采用CBF、MVDR、MUSIC、BPD和IHT算法進(jìn)行DOA估計(jì)。仿真結(jié)果如圖3所示，其中CBF、MVDR和MUSIC算法的快拍數(shù)為50，BPD和IHT算法采用單快拍。從圖3中可以看出，CBF和MVDR方法會(huì)在-65.3°出現(xiàn)柵瓣，雖然IHT算法仍然存在柵瓣，但BPD算法在一定程度上可以準(zhǔn)確估計(jì)目標(biāo)方位。因此基于CS的DOA估計(jì)可以在單快拍的條件下準(zhǔn)確重構(gòu)出目標(biāo)方位，且分辨力高于CBF和MVDR。

圖3 單目標(biāo)DOA估計(jì)

然后考慮存在兩個(gè)相干源，入射方位分別為0°和5°，信噪比為-5dB，其余條件與之前一致，結(jié)果如圖4所示。從圖4中可以看出基于CS的DOA估計(jì)方法可以在相干源、低信噪比的條件下基本準(zhǔn)確估計(jì)目標(biāo)方位，且分辨力高于經(jīng)典方法。

圖4 相干源DOA估計(jì)

4.2 實(shí)驗(yàn)結(jié)果

在消聲水池中進(jìn)行水下目標(biāo)方位估計(jì)的實(shí)驗(yàn)，目標(biāo)是一塊長(zhǎng)2.35m，寬0.5m的鐵板。利用收發(fā)合置形式的主動(dòng)聲納對(duì)目標(biāo)進(jìn)行探測(cè)，發(fā)射中心頻率為12kHz的PCW信號(hào)，脈沖寬度為10ms，用16元垂直線陣接收回波信號(hào)，陣元間距為0.1m，靜止目標(biāo)與接收陣之間的水平距離為25m。實(shí)驗(yàn)的目的是從接收到的回波信號(hào)中估計(jì)目標(biāo)的空間方位，即目標(biāo)的俯仰角。分別利用CBF、MVDR、BPD和IHT算法對(duì)目標(biāo)的俯仰角進(jìn)行估計(jì)，結(jié)果如圖5所示。其中CBF和MVDR用了100個(gè)快拍，兩者估計(jì)得到的俯仰角均為-1.3°，BPD和IHT計(jì)算了50個(gè)快拍，BPD估計(jì)的俯仰角均值為-1.3°，IHT估計(jì)的俯仰角為-1.5°，這些方法估計(jì)得到的結(jié)果與實(shí)際相符。可以看出，基于CS的DOA估計(jì)方法可以在陣元間距大于半波長(zhǎng)的情況下有效估計(jì)目標(biāo)的空間方位，尤其是在低快拍數(shù)的情況下，且分辨力高于經(jīng)典的CBF和MVDR方法。

圖5 實(shí)驗(yàn)結(jié)果

5 結(jié)語

本文將CS應(yīng)用到線陣的DOA估計(jì)上，尤其是陣元間距大于半波長(zhǎng)的線陣，把DOA估計(jì)轉(zhuǎn)變?yōu)橐粋€(gè)欠定線性方程，闡述了基本原理和基本的重構(gòu)算法BPD和IHT，并分析了傳感矩陣的相干性。數(shù)值仿真結(jié)果和實(shí)驗(yàn)處理結(jié)果均表明了BPD和IHT方法的有效性，尤其是在相干源、低信噪比、低快拍數(shù)的情況下，基于CS的DOA估計(jì)方法相較傳統(tǒng)的CBF和MVDR等方法具有更高的分辨力，能較好地估計(jì)目標(biāo)空間方位。