范文濤，章新華，康春玉，蔣 飚

（1.海軍大連艦艇學(xué)院 博士生隊(duì)，大連 116018；2.海軍大連艦艇學(xué)院 信息與通信工程系，大連 116018）

方位估計(jì)，或稱波達(dá)方向估計(jì)（Direction of Arrival，DOA）是被動(dòng)聲納表現(xiàn)目標(biāo)態(tài)勢(shì)的方式之一［1］，其中又以方位估計(jì)精度和方位分辨率作為衡量DOA估計(jì)性能的重要指標(biāo)。但是在淺海環(huán)境下，尤其是航道附近，往往存在各種強(qiáng)干擾源，嚴(yán)重影響了被動(dòng)聲納對(duì)于弱目標(biāo)的方位分辨性能和估計(jì)精度［1］，因此改善現(xiàn)有被動(dòng)聲納方位估計(jì)性能顯得尤為迫切。

常規(guī)波束形成［2］是被動(dòng)聲納中用于目標(biāo)方位估計(jì)的常用方法，其優(yōu)勢(shì)在于與環(huán)境無(wú)關(guān)，性能穩(wěn)健。但是受到孔徑的限制，使得常規(guī)波束形成的方位分辨率始終無(wú)法突破瑞利限［2］。Capon［2］和 Schmidt［3］先后提出的基于陣列觀測(cè)數(shù)據(jù)協(xié)方差矩陣的高分辨方位估計(jì)算法，即 MVDR［2］和 MUSIC［3］算法，彌補(bǔ)了常規(guī)波束形成在方位估計(jì)方面的缺點(diǎn)。然而受到低信噪比、陣形失配或者源數(shù)目的限制，使得高分辨類算法在處理實(shí)際被動(dòng)聲納數(shù)據(jù)時(shí)表現(xiàn)不夠穩(wěn)健［4］。

近些年發(fā)展起來(lái)的盲信號(hào)處理［5］技術(shù)得到了研究者的廣泛關(guān)注，其中尤以對(duì)盲源分離技術(shù)［6-8］的研究居多。由于盲源分離算法所具有的特殊優(yōu)點(diǎn)［6］，一些研究者將其引入到了水聲信號(hào)處理中［7］。文獻(xiàn)［6-9］成功地從多個(gè)輻射噪聲混合后的信號(hào)中分離出了所需源信號(hào)，驗(yàn)證了盲源分離用于分離水下目標(biāo)源的有效性，同時(shí)指出盲源分離算法可以分離空間上幾乎重疊的多目標(biāo)源［9］。然而，這些方法同時(shí)存在對(duì)于海洋多途信道影響下的寬帶卷積混合源分離［7］研究較少的局限性。日本Sawatari課題組［10］針對(duì)室內(nèi)多途混響下的源分離問(wèn)題，提出了多種波束形成與盲源分離融合的處理結(jié)構(gòu)，但是這些方法中的目標(biāo)DOA估計(jì)均是為解決寬帶頻點(diǎn)間次序模糊而服務(wù)的。范文濤［11］、康春玉［12］針對(duì)被動(dòng)聲納DOA估計(jì)與波形恢復(fù)，分別提出了將高分辨DOA估計(jì)算法與基于卷積混合的頻域盲源分離算法進(jìn)行融合，改善了傳統(tǒng)高分辨算法的DOA估計(jì)性能?？墒菍?duì)于寬帶頻域盲源分離算法存在的頻點(diǎn)間次序模糊和頻點(diǎn)內(nèi)尺度模糊問(wèn)題［9］，依舊沒(méi)有較好的解決辦法。

基于時(shí)間結(jié)構(gòu)的盲源分離算法［11］使用多個(gè)時(shí)間延遲二階相關(guān)矩陣來(lái)完成目標(biāo)源的提取。但是這類算法同時(shí)存在時(shí)滯量選取難，對(duì)沒(méi)有時(shí)間結(jié)構(gòu)特征的聲源的分離效果差等問(wèn)題［11］。近期提出的非參數(shù)化獨(dú)立分量分析算法［13］可以對(duì)目標(biāo)源的統(tǒng)計(jì)特性沒(méi)有任何先驗(yàn)知識(shí)的前提下完成目標(biāo)源的提取，是一種真正意義上的盲處理算法?？紤]將兩者融合，既利用了目標(biāo)源的時(shí)間結(jié)構(gòu)特性，又可以在不需預(yù)知目標(biāo)源統(tǒng)計(jì)特性的前提下完成源信號(hào)的分離，為了仿真比較方便，本文算法簡(jiǎn)稱為時(shí)間結(jié)構(gòu)非參數(shù)盲分離算法（Blind Temporal Structure Non-parametric Separation，BTSNPSEP）。本文根據(jù)寬帶卷積混合［14］，在每個(gè)子帶內(nèi)將基于時(shí)間結(jié)構(gòu)與非參數(shù)化特性的兩種盲源分離算法進(jìn)行融合，建立新的解混矩陣代價(jià)函數(shù)，然后根據(jù)非線性最優(yōu)化得到的解混矩陣得到子帶內(nèi)方位能量譜，最后仿照寬帶頻域波束形成對(duì)所有子帶方位能量譜累加，得到總的寬帶方位譜。通過(guò)計(jì)算機(jī)仿真與實(shí)際海試數(shù)據(jù)，分別同經(jīng)典MUSIC算法和MVDR方法進(jìn)行了比較，結(jié)果表明本文方法能夠用于寬帶被動(dòng)方位估計(jì)，尤其是強(qiáng)干擾背景下的弱目標(biāo)檢測(cè)。

1 寬帶盲源分離算法

1.1 問(wèn)題描述

考慮滿足半波長(zhǎng)且陣元間距為d的均勻線列陣（Uniform Iinear Array，ULA），陣元數(shù)為N，各陣元之間無(wú)互耦。遠(yuǎn)場(chǎng)寬帶目標(biāo)源數(shù)目小于陣列陣元數(shù)目，即（Q＜N）。由此可以得到第i個(gè)陣元的卷積混合為［12］：

其中：htq（l）表示從目標(biāo)源q到陣元i之間的沖激響應(yīng)；P表示沖激響應(yīng)的長(zhǎng)度，ni（t）為疊加在每個(gè)陣元上的噪聲。

由于大部分盲源分離算法是基于瞬時(shí)混合模型的［7］，因此通過(guò)短時(shí)傅里葉變換 （Short-time Fourier transform，STFT）將時(shí)域卷積混合x(chóng)i（t）變換為頻域信號(hào)xi（k，fj），可以表示為：

其中：fj表示第j個(gè)頻點(diǎn)，j=1，…，J；k=1，…，K表示頻域快拍數(shù)。aiq（fj）為從目標(biāo)源q到陣元i之間的頻率響應(yīng)。sq（k，fj），ni（k，fj）分別為sq（t）和ni（t）的短時(shí)傅里葉變換。由此可以得到第j個(gè)頻點(diǎn)的陣列頻域觀測(cè)數(shù)據(jù)快拍矩陣［12］：

其中：X（k，fj）=［x1（k，fj），…，xN（k，fj）］T為陣列觀測(cè)數(shù)據(jù)矩陣，S（k，fj）=［s1（k，fj），…，sQ（k，fj）］T為源信號(hào)矩陣，N（k，fj）=［n1（k，fj），…，nN（k，fj）］T為加性噪聲矩陣。A（fj，Θ）=［aθ1（fj），…，aθQ（fj）］T表示陣列流形矩陣，其中aθ（fj）為：

1.2 基于時(shí)間結(jié)構(gòu)的非參數(shù)盲源分離算法

根據(jù)文獻(xiàn)［13］和式（3）可知，傳統(tǒng)瞬時(shí)盲源分離算法通過(guò)求取解混矩陣來(lái)完成對(duì)目標(biāo)源的分離和混合系統(tǒng)的辨識(shí)，即：

其中：Y（k，fj）=［y1（k，fj），…，yN（k，fj）］T為恢復(fù)出的目標(biāo)源信號(hào)。W（fj）表示解混矩陣。

圖1為本文所提算法框架。如圖1所示，利用每個(gè)子帶內(nèi)估計(jì)的解混矩陣W（fj）和目標(biāo)源信號(hào)Y（k，fj），求取目標(biāo)方位能量譜，最后將所有子帶內(nèi)方位能量譜累加，即可得到總的寬帶方位能量譜。

圖1 算法框架Fig.1 The framework of the algorithm

文獻(xiàn)［11］介紹了利用觀測(cè)數(shù)據(jù)的時(shí)間結(jié)構(gòu)特性進(jìn)行盲源分離的思想，其代價(jià)函數(shù)：

其中：τb，b=1，…，B表示頻點(diǎn)fj處的第b個(gè)時(shí)間延遲量，又稱時(shí)滯量。

根據(jù)文獻(xiàn)［13］，利用得到基于信息論準(zhǔn)則的盲源分離代價(jià)函數(shù)：

其中：H（·）表示信息熵，H（X）為關(guān)于解混矩陣W的常量。

將式（6）與式（7）所表示的代價(jià)函數(shù)融合，可以得到基于時(shí)間延遲結(jié)構(gòu)的盲源分離最優(yōu)化問(wèn)題：

其中：I（yi（k，fj）；yi（k+τb，fj））為：

由式（7）和（8）易得式（8）所對(duì)應(yīng)的代價(jià)函數(shù)：

其中：h為內(nèi)核帶寬［13］，φ（·）為高斯核［13］：

Yim為內(nèi)核質(zhì)心［13］：

其中：X（m）（fj）表示每個(gè)子帶內(nèi)觀測(cè)數(shù)據(jù)矩陣的第m個(gè)頻域快拍。Wi（fj）表示解混矩陣W（fj）的第i行向量。

根據(jù)式（5），則yi（k，fj）=Wi（fj）X（fj），將式（13）代入式（11），可得：

目標(biāo)函數(shù)（15）是解混矩陣W的非線性函數(shù)，約束（16）將最優(yōu)化問(wèn)題的可能解空間限定到一個(gè)有限集內(nèi)。目標(biāo)函數(shù)（15）可以通過(guò)牛頓法對(duì)代價(jià)函數(shù)（10）的非線性最優(yōu)化解混矩陣Wopt，最優(yōu)化求解過(guò)程同文獻(xiàn)［12］類似。由此估計(jì)出的源信號(hào)Y（k，fj）=WoptX（k，fj）。

對(duì)于盲源分離恢復(fù)的結(jié)果，每一路輸出的方位都能估計(jì)出來(lái)：

由式（17）在每個(gè)頻點(diǎn)fj處估計(jì)出的N個(gè)θi（fj）值，將每個(gè)θi（fj）值在θ=-90°，…，90°方位空間中進(jìn)行搜索匹配，找到其對(duì)應(yīng)的方位位置，該方位位置的能量值為該路恢復(fù)出的源信號(hào)的能量。因此可以在每個(gè)子帶內(nèi)構(gòu)造出類似于波束形成的方位能量譜：

對(duì)于式（18）需要說(shuō)明的是，每個(gè)子帶內(nèi)盲源分離算法估計(jì)的N個(gè)θi（fj）值有可能會(huì)出現(xiàn)重復(fù)的情況，本文采取累加的辦法。頻域卷積盲源分離算法普遍存在頻點(diǎn)間次序和頻點(diǎn)尺度模糊的問(wèn)題。然而式（18）并不受這兩個(gè)問(wèn)題的影響，即使是相鄰頻點(diǎn)間盲源分離輸出次序不一致，對(duì)于最后的方位能量譜的構(gòu)建不構(gòu)成影響。仿照頻域?qū)拵Рㄊ纬傻姆窍喔衫奂拥乃枷?，最后總的方位能量譜為：

2 數(shù)據(jù)檢驗(yàn)

2.1 仿真數(shù)據(jù)

寬帶仿真實(shí)驗(yàn)考慮均勻線列陣位于兩寬帶目標(biāo)源的遠(yuǎn)場(chǎng)，兩目標(biāo)源均用海上實(shí)錄艦船輻射噪聲代替。兩目標(biāo)源的原始時(shí)域波形與頻譜如圖2所示。寬帶分析頻段為800～1 600 Hz，陣元數(shù)為32，陣元間距滿足分析頻段最高頻率對(duì)應(yīng)的半波長(zhǎng)。采樣頻率為25 kHz。環(huán)境噪聲添加方式使用 MATLAB中 Awgn函數(shù)［12］對(duì)寬帶陣列數(shù)據(jù)整體添加加性高斯白噪聲。

圖2 兩目標(biāo)源原始時(shí)域波形與頻譜Fig.2 The original time domain waveform and frequency spectrum of the two targets

2.1.1 不同信噪比下方位分辨率比較

兩個(gè)目標(biāo)等功率，分析數(shù)據(jù)長(zhǎng)度為1 s，變換不同環(huán)境噪聲信噪比（-20 dB～10 dB），采用兩種方位間隔，測(cè)試BTSNPSEP與MVDR的方位分辨性能，其中每個(gè)信噪比處的估計(jì)運(yùn)行100次蒙特卡洛仿真［12］。圖3為8°目標(biāo)與11°目標(biāo)（3°間隔）在不同信噪比條件下的方位分辨性能。圖3（a）為BTSNPSEP方法結(jié)果，圖3（b）為MVDR方法結(jié)果。比較圖2中兩個(gè)子圖可以看出，BTSNPSEP方法在低信噪比時(shí)的分辨率要好于MVDR方法。

圖4為8°目標(biāo)與13°目標(biāo)（5°間隔）在不同信噪比條件下的方位分辨性能。圖4（a）為BTSNPSEP方法結(jié)果，圖4（b）為MVDR方法結(jié)果。比較圖3中兩個(gè)子圖可見(jiàn)當(dāng)目標(biāo)間距拉大以后，MVDR的性能變得較好，BTSNPSEP能夠在低信噪比時(shí)能夠區(qū)分兩個(gè)目標(biāo)，但對(duì)于13°目標(biāo)的估計(jì)誤差了1°。

2.1.2 不同目標(biāo)源數(shù)目下方位分辨率比較

考慮到陣元數(shù)較多時(shí)的BTSNPSEP方法運(yùn)算速度問(wèn)題，通過(guò)加入白化預(yù)處理步驟［7］，使得BTSNPSEP方法也涉及到源數(shù)目設(shè)置的問(wèn)題，由此同MUSIC方法比較不同源數(shù)目設(shè)置時(shí)的方位分辨率。兩個(gè)目標(biāo)依舊等功率，同時(shí)環(huán)境噪聲信噪比固定為-5 dB，其中每個(gè)源數(shù)目設(shè)置處的估計(jì)運(yùn)行100次蒙特卡洛仿真［12］。圖5為8°目標(biāo)與11°目標(biāo)（3°間隔）在不同源數(shù)目設(shè)置下的方位分辨性能。圖5（a）為BTSNPSEP方法結(jié)果，圖5（b）為MVDR方法結(jié)果。圖6為8°目標(biāo)與13°目標(biāo)（5°間隔）在不同源數(shù)目設(shè)置下的方位分辨性能。圖6（a）為BTSNPSEP方法結(jié)果，圖6（b）為MVDR方法結(jié)果。從圖5（a）可以看出，當(dāng)方位間隔3°時(shí)，源數(shù)目設(shè)置的變化對(duì)于BTSNPSEP方法影響還是非常大的，對(duì)11°目標(biāo)的估計(jì)不如MUSIC算法。從圖6可以看出當(dāng)方位間距擴(kuò)大（5°間隔），BTSNPSEP方法能夠在各種源數(shù)目設(shè)置下區(qū)分兩個(gè)目標(biāo)，尤其是當(dāng)源數(shù)目為1個(gè)時(shí)（欠定情況）要好于MUSIC。

2.1.3 信干比為-6 dB下弱目標(biāo)檢測(cè)能力

設(shè)定8°聲源為干擾源，11°聲源為所需目標(biāo)源，信干比［12］為 -6 dB，環(huán)境噪聲信噪比固定為 -5 dB，BTSNPSEP采用預(yù)白化處理源數(shù)目設(shè)置為3，MUSIC方法同之。共進(jìn)行20次仿真，每次為100次蒙特卡洛仿真結(jié)果。圖7為分別為BTSNPSEP、MVDR、MUSIC三種方法估計(jì)結(jié)果。比較圖7三個(gè)圖可以看出BTSNPSEP對(duì)于弱目標(biāo)的檢測(cè)效果要好于MVDR與MUSIC。

圖7 信干比為-6 dB下弱目標(biāo)檢測(cè)能力比較Fig.7 The comparison of weak target detection capability under SIR=-6dB.

2.2 海上實(shí)測(cè)數(shù)據(jù)

試驗(yàn)海區(qū)位于東經(jīng)121°33＇～121°42＇，北 38°46＇～38°52＇的主航道附近，海區(qū)平均水深46 m，聲傳播速度約1 500 m/s。海區(qū)地形復(fù)雜，水流急，過(guò)往船只多，離岸近、干擾源多，屬典型的復(fù)雜海區(qū)［11］。接收陣列由28個(gè)水聽(tīng)器等間隔組成，陣元間距0.225 m，接收陣列深度約為20 m。水聽(tīng)器接收信號(hào)經(jīng)信號(hào)調(diào)理機(jī)后到Sony sir1000i錄音機(jī)磁帶記錄，整個(gè)記錄數(shù)據(jù)期間，接收船一直輔機(jī)發(fā)電，有較大噪聲。試驗(yàn)時(shí)目標(biāo)船在大約5.3 km的距離沿正橫經(jīng)過(guò)接收船［11］。同時(shí)，目標(biāo)船運(yùn)動(dòng)時(shí)在視覺(jué)范圍內(nèi)還發(fā)現(xiàn)有漁船目標(biāo)運(yùn)動(dòng)，而且漁船先于目標(biāo)船經(jīng)過(guò)接收船正橫。

圖8為三種方法的估計(jì)得方位歷程圖比較，其中BTSNPSEP與MUSIC的源數(shù)目均設(shè)置為20。從圖中可以看出對(duì)于-80°～-40°區(qū)域的弱干弱目標(biāo)BTSNPSEP的估計(jì)結(jié)果相對(duì)清晰。第0～100 s時(shí)間，對(duì)于20°～40°區(qū)域的目標(biāo)船和漁船的估計(jì)BTSNPSEP方法也比MUSIC、MVDR方法要清晰，空間增益提高比較明顯。

圖8 三種方法海試數(shù)據(jù)比較Fig.8 The comparison of sea trial result during three methods

3 結(jié)論

本文針對(duì)被動(dòng)聲納信號(hào)的特殊性，結(jié)合兩種盲源分離算法的特點(diǎn)，提出了一種融合時(shí)間延遲結(jié)構(gòu)和非參化特性的盲源分離算法（BTSNPSEP），并將其用于被動(dòng)聲納方位估計(jì)。分別利用寬帶仿真和海試數(shù)據(jù)與經(jīng)典MVDR和MUSIC高分辨方位估計(jì)方法進(jìn)行了比較，得到以下結(jié)論：

（1）BTSNPSEP方法能夠突破瑞利限的限制分辨空間接近的多目標(biāo)，且在低信噪下的方位分辨性能要好于MVDR方法。

（2）BTSNPSEP方法能夠在欠定條件下分辨出多目標(biāo)源。

（3）BTSNPSEP方法比起MVDR和MUSIC方法，對(duì)于強(qiáng)干擾背景下的弱目標(biāo)檢測(cè)具有一定的優(yōu)勢(shì)。

［1］Meng T Z，Buck J R.Rate distortion bounds on passive sonar performance［J］.IEEE Trans.on Acoustics.Speech and Signal Processing，2010，58（1）：326-336.

［2］ Capon J.High resolution frequency wavenumber spectrum analysis［C］//Proc.of the IEEE，1969，57（8）：1408-1418.

［3］Stoica P，Nehorai A.Music，maximum likelihood，and cram r-rao bound［J］.IEEE Trans.on Acoustics.Speech and Signal Processing，1989，37（5）：720-741.

［4］ Bao C Y.Adaptive beamforming for sonar audio［J］.Proc.of ACOUSTICS 2005.Australia：Australian Acoustical Society Press，2005：483-485.

［5］蔡艷平，李艾華，石林鎖，等.基于盲解卷積的柴油機(jī)振動(dòng)信號(hào)分離研究［J］.振動(dòng)與沖擊，2010，29（9）：38-41.

［6］Benchekroun N，Mansour A.A general structure for the separation of underwater acoustic signals［C］//Proc.of Symposium on Communications，Control and Signal Processing.Marrakech：IEEE Press，2006：273-278.

［7］De Moura，Seixas N N，F(xiàn)ilho J M.Independent component analysis for optimal passive sonar signal detection［C］//Seventh International Conference on Intelligent Systems Design and Applications.Rio de Janeiro：IEEE Press，2007：671-678.

［8］雷衍斌，李舜酩，門(mén)秀花，等.基于自相關(guān)降噪的混疊轉(zhuǎn)子振動(dòng)信號(hào)分離［J］.振動(dòng)與沖擊，2011，30（1）：218-222.

［9］張安清.盲分離技術(shù)及其在水聲信號(hào)中的應(yīng)用研究［D］.大連：海軍大連艦艇學(xué)院，2006.

［10］Saruwatari H，Kurita S，Takeda K.Blind source separation combining independent component analysis and beamforming［J］.EURASIP Journal on Applied Signal Processing，2003，1135-1146.

［11］康春玉，章新華，韓 東.盲源分離與高分辨融合的DOA估計(jì)與信號(hào)恢復(fù)方法［J］.自動(dòng)化學(xué)報(bào)，2010，36（3）：443-445.

［12］章新華，范文濤.波束形成與獨(dú)立分量分析融合的寬帶高分辨方位估計(jì)方法［J］.聲學(xué)學(xué)報(bào)，2009，34（4）：175-180.

［13］ Boscolo R，Pan H，Roychowdhury V P.Independent component analysis based on nonparametric density estimation［J］.IEEE Trans.on Neural Networks，2004，15（1）：55-65.

［14］范 濤，李志農(nóng)，岳秀廷.基于變分貝葉斯算法和MLP網(wǎng)絡(luò)的后非線性混合盲源分離方法研究［J］.振動(dòng)與沖擊，2010，29（6）：21-24.