宋其巖 馬曉川 李 璇 詹 飛

（1.中國科學(xué)院聲學(xué)研究所中科院水下航行器信息技術(shù)重點(diǎn)實(shí)驗(yàn)室，北京 100190；2.中國科學(xué)院大學(xué)，北京 100049）

1 引言

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

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

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

2 信號(hào)模型

假設(shè)有K個(gè)遠(yuǎn)場窄帶信號(hào)以角度θk，k=1，2，...，K入射到均勻直線陣列，陣元間距d=λ2，陣元m的輸出可表示為

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

將式（1）寫成矢量形式

3 方位估計(jì)算法

3.1 常規(guī)波束形成

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

各陣元的權(quán)矢量，

陣列輸出為，

則整個(gè)陣列輸出的平均功率為

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

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

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

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

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

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

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

3.3 最小方差無失真響應(yīng)波束形成

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

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

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

3.4 解卷積MVDR

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

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

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

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

3.4.1 Richardson-Lucy解卷積MVDR

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

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

文獻(xiàn)［21］應(yīng)用矢量外推算法提出了一種快速RL算法，本文仿真計(jì)算均是基于快速RL算法，詳細(xì)執(zhí)行過程請(qǐng)查看文獻(xiàn)［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 算法的使用方式完全一樣，算法的輸入和輸出完全一樣，均是首先計(jì)算方位譜，然后采用解卷積算法利用PSF對(duì)方位譜解卷。

4 仿真實(shí)驗(yàn)

4.1 單次仿真實(shí)驗(yàn)

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

其中，代表信號(hào)能量，代表噪聲能量。

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

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

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

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

4.2 時(shí)間方位歷程譜圖仿真實(shí)驗(yàn)

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

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

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

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

4.3 蒙特卡洛統(tǒng)計(jì)實(shí)驗(yàn)

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

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

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

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

5 結(jié)論

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