宋其巖 馬曉川 李 璇 詹 飛

（1.中國科學院聲學研究所中科院水下航行器信息技術重點實驗室，北京 100190；2.中國科學院大學，北京 100049）

1 引言

來波方向（Direction of Arrival，DOA）估計一直是陣列信號處理的重要內(nèi)容。聲波頻率越高，介質(zhì)吸收損失越大，為了使聲納的作用距離增大，需要盡可能降低探測信號的頻率，根據(jù)半波長布陣原則，頻率越低陣列的陣元間距往往越大，陣元數(shù)量越少，波束主瓣越大，旁瓣造成的能量泄露越嚴重，陣列的分辨能力越差［1］。因此，如何在尺寸固定的聲納平臺，盡可能提高方位估計的分辨力和估計精度，一直是陣列信號處理的難點問題。另外，隨著主動傳感器陣列作用范圍的增大，搜索范圍內(nèi)很可能會出現(xiàn)多個目標，例如海上航母護航編隊、拖曳誘餌、水下航行器集群，如何對多目標進行準確、快速的方位估計亟待解決。在空氣聲領域，基于麥克風陣列的室內(nèi)說話人聲源定位、語音增強同樣需要應用到多目標方位估計技術［2-6］。

經(jīng)過幾十年的發(fā)展，眾多學者提出了眾多方位估計算法［7-10］。1948 年，Bartlett 提出傳統(tǒng)波束形成技術（Conventional Beamforming，CBF）［11］，CBF算法將各個陣元的接收數(shù)據(jù)進行延時求和（Delay and Sum，DAS），對于不同的波達方向，陣元的加權值不同，然后求波束輸出功率，功率譜的峰值對應的角度就是真實的波達方向。CBF 算法計算速度快，在工程中得到了普遍的應用，但是由于CBF算法的波束主瓣較寬，分辨能力較差，旁瓣較高引起能量泄露，功率強的信號容易掩蔽功率低的信號，最終導致模糊效應［12-13］。1969 年Capon 提出了高分辨的最小方差無失真響應（Minimum Variance Distortionless Response，MVDR）波束形成器［14］。1986 年Schmidt 提出了高分辨多重信號分類算法（Multiple Signal Classification，MUSIC）［15］，MUSIC算法需要對數(shù)據(jù)采樣協(xié)方差矩陣執(zhí)行特征分解，大特征值對應的特征向量是信號主特征向量，通過小特征值對應的特征向量構造噪聲子空間，利用信號子空間與噪聲子空間的正交性實現(xiàn)波達方向高分辨，但是MUSIC 估計的空間功率譜屬于偽譜，無法準確估計出聲源的功率。1986 年Roy 等人利用規(guī)則陣列的旋轉不變性提出了旋轉不變參數(shù)估計技術（Estimating Signal Parameter via Rotational Invariance Techniques，ESPRIT）［16］，計算復雜度比MUSIC 算法低，但是分辨能力低于MUSIC。繼解卷積后處理技術在圖像恢復領域取得了成功的應用后，有學者將解卷積技術擴展應用到了參數(shù)估計領域，2018 年，T.C.Yang 將常規(guī)波束形成表示為卷積形式，并將觀察方向為0°的波束圖作為點擴散函數(shù)（Point Spreading Function，PSF），將Richardson-Lucy（RL）［17］算法應用于CBF 算法進行解卷積后處理，提高了分辨能力，2019 年謝磊等人對MUSIC 的方位譜進行后處理，獲得了背景級更低的方位譜［18］。Ma 等人提出了聲源圖像的快速解卷積方法（Deconvolution Approach for the Mapping of Acoustic Source，DAMAS）［19］，2017 年，Lylloff 將快速迭代收縮閾值算法（Fast Iterative Shrinkage-Thresholding Algorithm，F(xiàn)ISTA）解卷積技術應用于平面麥克風陣列聲源定位中［20］，并指出如果陣列具有平移不變性則可以采用快速傅里葉變換（Fast Fourier Transform，F(xiàn)FT）執(zhí)行解卷積計算，提高了解卷過程中的計算效率，這類后處理算法的顯著優(yōu)勢是能夠繼承被解卷積算法的優(yōu)勢，在付出少量的計算代價下，能夠提高原始算法的分辨性能。

MVDR 算法分辨能力高，估計出的方位譜為功率真實譜，然而低信噪比（Signal-to-Noise Ratio，SNR）下，方位譜背景級較高，估計精度較差。本文對MVDR 算法的方位譜重新表述為卷積模型，并運用兩種解卷積算法對MVDR 方位譜進行解卷積，算法在付出少量的計算代價下，能夠提高MVDR 算法的方位分辨能力。仿真實驗顯示提出的算法具有高分辨能力。

2 信號模型

假設有K個遠場窄帶信號以角度θk，k=1，2，...，K入射到均勻直線陣列，陣元間距d=λ2，陣元m的輸出可表示為

其中，陣元數(shù)為M，陣元號m=0，1，...，M-1，快拍數(shù)為N，n=1，2，...，N。

將式（1）寫成矢量形式

3 方位估計算法

3.1 常規(guī)波束形成

經(jīng)典波束形成法（Conventional Beamforming，CBF）又稱延遲相加法（Delay and Sum，DAS），是DOA 領域最基礎的方法，通過權系數(shù)對各個傳感器記錄的數(shù)據(jù)執(zhí)行加權求和，從而達到調(diào)整陣列的波束指向的目的，令波束主瓣指向所期望的角度方向。

各陣元的權矢量，

陣列輸出為，

則整個陣列輸出的平均功率為

對波達方向在[-90°，90°]范圍內(nèi)進行角度掃描，計算陣列輸出功率，最大功率點所對應的角度θs即為聲源的DOA。

3.2 Richardson-Lucy解卷積常規(guī)波束形成

由于傳統(tǒng)方位估計CBF 的分辨率受到陣列陣元數(shù)量的限制，即存在瑞利限。因此對CBF 的波束輸出進行后處理是提高CBF 分辨能力的有效技術途徑，文獻［17］采用Richardson-Lucy 算法對CBF 方位譜進行解卷積后處理（CBF Richardson-Lucy，CBF-RL），提高了CBF算法的性能。

重新考慮常規(guī)波束形成方位譜公式（6），

將式（7）重新表示為卷積的形式，

PSFCBF(sinθ)是CBF 算法的PSF，PSF 的物理意義是空間中的點光源經(jīng)過成像系統(tǒng)沖激響應函數(shù)卷積后的擴散后的圖像，如圖1 所示，點光源經(jīng)過采樣函數(shù)卷積后形成漣漪狀的擴散圓環(huán)圖像。在圖像恢復領域，Richardson-Lucy（RL）是一種最流行的解卷積技術，RL 迭代解卷積過程如式（10）所示［21］，

其中，S=S(sinθ)，sinθ∈[-1，1]，Si代表第i次估計 結果，PCBF=P(sinθ)，sinθ∈[-1，1]，PSFCBF是CBF 算法的波束圖，波束形狀與sinθ=0 處的單目標的功率譜圖相同。式（10）的物理意義是PSFCBF?Si是第i次迭代對PCBF的擬合，即預測的模糊圖像，Perror=是PCBF與其預測圖像之間的誤差的度量，是以點擴散函數(shù)PSFCBF為加權值對估計誤差Perror的加權求和，當前估計Si小于最優(yōu)解時，下一次迭代Si+1會變大，當前估計Si大于最優(yōu)解時，下一次迭代Si+1會變小，以趨近最優(yōu)解。

3.3 最小方差無失真響應波束形成

傳統(tǒng)波束形成聚焦于目標信號的相干累加，沒有考慮其他方位信號（視為干擾）和噪聲的信息，當不存在其他方位信號且環(huán)境噪聲為白噪聲時，則DAS 波束形成具有最大的輸出信噪比。當存在其他方位信號或環(huán)境噪聲非白噪聲時，如能有效利用干擾和噪聲的統(tǒng)計信息則會實現(xiàn)比DAS更高的輸出信干噪比。Capon 法或MVDR 的角度分辨率好于CBF，該算法將陣列自由度分成兩部分，其中一部分用于在期望角度處形成波束，另一部分自由度用于排除干擾的影響，即在非期望信號方向設計較深的凹陷。Capon 法的準則是讓陣列的輸出功率達到最小同時又保證在期望角度方向上的增益維持不變，因此可最大程度地排除干擾的影響［14］。即

運用拉格朗日算子可將上式轉化為無約束最優(yōu)化模型，其解為：

從上式可以看出，權值W依賴于數(shù)據(jù)協(xié)方差矩陣R的逆，因此MVDR屬于數(shù)據(jù)依賴的自適應算法。

3.4 解卷積MVDR

對現(xiàn)有算法進行解卷積后處理是提升算法性能的有效技術，現(xiàn)把式（12）表示為

將式（14）代入式（13），

其中，σ2是噪聲功率δ(sinθ-sinθk)是多目標角度的真實空間分布，PSFMVDR=|WH(sinθ)a(0)|2是MVDR的PSF。

根據(jù)PSF 的物理意義可知，一種可行的確定PSFMVDR的方法是：假設感興趣區(qū)間的中心位置有一個點目標，目標的功率為MVDR 方位譜譜峰的平均值，即根據(jù)式PMVDR(sinθ)=計 算PSFMVDR，其 中R=P0a(0)a(0)H+I，P0是MVDR 估計出的聲源的功率的平均值，?·?1表示歸一化操作。圖3直觀顯示了CBF算法與MVDR算法的PSF，顯然，CBF 的方位譜主瓣更寬，旁瓣更高，分辨能力差，MVDR 算法的PSF 更加尖銳，分辨能力更高。

3.4.1 Richardson-Lucy解卷積MVDR

在圖像處理領域RL 算法是一種已被理論證實和工程實際證明有效的后處理解卷積算法，RL解卷積MVDR迭代公式如下所示，

其中，PMVDR=PMVDR(sinθ)，sinθ∈[-1，1]。

文獻［21］應用矢量外推算法提出了一種快速RL算法，本文仿真計算均是基于快速RL算法，詳細執(zhí)行過程請查看文獻［21］。

RL解卷積MVDR（MVDR-RL）的執(zhí)行流程框圖如下所示：

3.4.2 FISTA解卷積MVDR

Oliver 等人提出可以采用快速迭代收縮閾值算法（Fast Iterative Shrinkage-Thresholding Algorithm，F(xiàn)ISTA）完成解卷積過程［20］，F(xiàn)ISTA 解卷積MVDR（MVDR-FISTA）算法執(zhí)行流程圖如下所示。注意RL算法與FISTA 算法的使用方式完全一樣，算法的輸入和輸出完全一樣，均是首先計算方位譜，然后采用解卷積算法利用PSF對方位譜解卷。

4 仿真實驗

4.1 單次仿真實驗

本節(jié)通過仿真分析提出的兩種解卷積MVDR算法的性能。假定陣列為半波長布陣的均勻直線陣列，陣元數(shù)為M=10，兩個不等強度的非相干聲源波達方向為sinθ1=-0.08，sinθ2=0.1，SNR 分別為5 dB 和3 dB?？炫臄?shù)為100（對于窄帶信號一個快拍數(shù)據(jù)是一個時域采樣點，采樣率為100 kHz時，100 個快拍數(shù)據(jù)對應1 ms 時長），為防止卷繞誤差，將sinθ∈[-1，1]擴展為[-2，2]。其中SNR 的定義為

其中，代表信號能量，代表噪聲能量。

CBF 算法是陣列信號處理領域的常用算法，其波束主瓣的3 dB 寬度約為110o，由于CBF 算法主瓣寬、旁瓣高，當多個目標的角度間隔小于等于波束主瓣時，CBF 算法無法分辨［1］。通過圖4 可以看出CBF 算法無法準確估計這兩個聲源信號的DOA。CBF-RL 算法通過解卷積運算，降低了CBF算法寬主瓣，高旁瓣的不利影響，能夠成功分辨出兩個信號，但是估計的聲源位置有偏移。MVDR 算法能夠分辨出兩個聲源，但是由于功率譜圖的譜峰和譜谷的差值為6 dB，沒有將兩個譜峰完全分離。CBF-RL，MVDR-RL，MVDR-FISTA 三種算法能夠準確分辨兩個聲源，由于采用解卷積后處理，兩個聲源的譜峰完全分離。其中MVDR-RL 和MVDRFISTA算法的方位譜主峰寬度一致，為CBF-RL算法主峰寬度的16。

將上述實驗進行1000次仿真，統(tǒng)計不同算法的迭代次數(shù)和計算時間如表1所示。

表1 不同算法計算效率（Intel（R）CPU：i5-1035G1@3.4 GHz）Tab.1 Computational efficiency of different algorithms

表1 表明，CBF-RL 算法的計算時間最短，但是對比圖4，CBF-RL 算法的并沒有準確估計目標方位，因此沒有尋找到全局最優(yōu)值。MVDR-RL 的迭代次數(shù)少于MVDR-FISTA 算法，兩種算法的迭代策略不同，RL 算法采用矢量加速的算法進行迭代，更容易尋找到最優(yōu)值，F(xiàn)ISTA 算法采用梯度下降思想進行迭代計算，收斂速度較慢。但是由于RL 算法需要計算加速因子，因此平均每次迭代的計算時間比FISTA算法的計算時間更長。

4.2 時間方位歷程譜圖仿真實驗

假設有兩個不等強度的非相干目標，其中目標1 的角度sinθ=-0.05，SNR=5 dB，目標2 的角度sinθ從-0.5 到0.5 變化，SNR=3 dB，每次方位估計取100 個快拍數(shù)據(jù)，以下是不同算法的時間方位譜估計結果。

從圖5可以看出，與其他算法相比較，CBF 算法的功率譜圖的主瓣寬，當運動目標角度在sinθ=±0.2范圍內(nèi)時，兩個目標不可分，而且旁瓣較高較多。圖6 表明CBF-RL 算法能夠提高CBF 算法的分辨能力，這是由于解卷積后處理可以降低CBF 算法波束主瓣寬，旁瓣高的不利影響，CBF-RL 主瓣更窄，能夠去除大部分旁瓣，功率譜的背景能量更低。當運動目標角度在sinθ=±0.15 范圍內(nèi)，CBF-RL 算法難以分辨兩個目標。圖7表明MVDR算法的估計結果好于CBF 算法，這是由于MVDR 的優(yōu)勢是抑制干擾，屬于高分辨算法。圖8 和圖9 是分別采用RL 算法和FISTA 算法對MVDR 的功率譜進行解卷積后處理，兩者差別不大，都能夠提高MVDR 的性能，其分辨能力高于CBF 算法、CBF-RL 算法和MVDR 算法。解卷積MVDR 的背景級更低，幾乎沒有旁瓣。另外，在方位交匯附近，由于兩個目標太近，此時所有算法均無法準確分辨兩個聲源，因此會形成模糊，并最終兩個譜峰合并為一個譜峰。

通過上述仿真實驗可以發(fā)現(xiàn)，CBF 算法通過旁瓣產(chǎn)生能量泄露，因此當兩個目標相互接近時，方位估計結果會相互干擾，影響分辨性能。MVDR 算法對干擾會自適應形成凹陷，即形成空域陷波器，因此即使是兩個相近的目標，也不會通過能量泄露，彼此干擾，本文從MVDR 方位估計的結果著手分析，實際上隱含地利用了MVDR 抑制干擾這一優(yōu)勢。仿真結果顯示，解卷積MVDR 優(yōu)于解卷積CBF算法。

毋庸置疑，解卷積CBF 算法理論非常完美，提高了CBF 算法的分辨能力，因此在水聲方位估計領域，解卷積思想得到越來越多的重視。解卷積MVDR 盡管只是一種近似表示，但是仿真結果表明解卷積MVDR 算法提高了MVDR 分辨能力，這是由于MVDR 會抑制干擾，屬于高分辨算法，對其進行解卷積，能夠進一步提升分辨能力，其性能好于解卷積CBF。以下通過蒙特卡洛仿真實驗進行詳細說明。

4.3 蒙特卡洛統(tǒng)計實驗

為進一步定量統(tǒng)計不同算法的均方根誤差（Root Mean Square Error，RMSE），假設兩個等強度的非相干信號，波達方向為sinθ1=-0.08，sinθ2=0.1，SNR 從-10 dB～10 dB 范圍內(nèi)變化，對于每個SNR 取值，進行100 次蒙特卡洛仿真實驗，其中RMSE的定義為，

當方位譜的峰值與波谷相差3 dB 時認為成功分辨出了兩個聲源。在低SNR 條件下，算法存在無法以概率1 分辨兩個聲源的情況，本文將無法分辨定義為：1.兩個較近的目標的波峰合為一個主峰，無法通過主峰的位置確定兩個目標的波達方向，如圖4 中的CBF 算法的方位譜；2.波峰與兩個目標的波峰之間的波谷的差值小于3 dB。當?shù)趐次蒙特卡洛仿真實驗無法分辨兩個目標時，此時方位譜最高峰對應的角度為θpeak，此時令目標1，2 的估計角度

圖10 是不同算法的分辨概率與SNR 的關系曲線，所有算法的分辨概率隨SNR 的增大而增大。MVDR 算法的分辨概率最差，主要是由于在低SNR環(huán)境下，波谷與波峰之間的差值小于3 dB，分辨概率低。CBF-RL、MVDR-RL、MVDR-FISTA 算法由于采用解卷積后處理，消除模糊效應，因此分辨率更高。MVDR-RL 算法的分辨概率在低SNR 條件下，高于MVDR-FISTA 算法。當SNR 大于-6 dB 時，MVDR 解卷積算法的分辨概率達到1，當SNR 大于-4 dB時，CBF-RL算法的分辨概率達到1。

圖11 顯示了不同算法在不同SNR 條件下的RMSE，由于CBF 算法無法分辨角度間隔小于瑞利限的聲源，因此沒有展示CBF 的RMSE 曲線。通過圖11 可以看出，本文提出的兩種解卷積MVDR 算法均具有比傳統(tǒng)MVDR 更低的RMSE，MVDR-RL和MVDR-FISTA 算法的估計性能非常接近。在低SNR 環(huán)境下，MVDR-RL 的分辨概率大于MVDRFISTA 算法，因此MVDR-RL 的估計精度明顯好于MVDR-FISTA 算法。當信噪比大于2 dB 時，MVDR算法和解卷積MVDR 算法的估計精度相同。圖12是不同算法的估計偏差，偏差的定義為：偏差=，由于偏差有正有負，因此正偏差和負偏差分開統(tǒng)計，圖12縱軸的正半軸是偏差為正的部分，負半軸是偏差為負的部分，圖12 結果與圖11 基本一致，MVDR-RL 和MVDR-FISTA 算法的效果最好，低SNR 下，MVDR 的性能最差，估計偏差最大。

5 結論

為提高MVDR 算法的方位分辨性能，本文將MVDR 算法的輸出功率譜重新表示為卷積的形式，并運用兩種解卷技術對MVDR 的方位譜進行后處理，提升了算法的性能。該方法將角度空間中心位置的單個聲源的MVDR 方位譜作為PSF，并利用RL算法和FISTA 算法分別對MVDR 的方位譜進行解卷積后處理，獲得具有更低背景級的MVDR-RL 和MVDR-FISTA 方位譜，同時提高了分辨能力和估計精度。仿真實驗顯示了所提算法的良好性能。