本發(fā)明涉及的是一種傳感器陣列信號(hào)處理方法，具體地說是一種聲矢量圓陣寬帶相干目標(biāo)方位估計(jì)方法。

背景技術(shù)：

水下目標(biāo)輻射噪聲可分為線譜及寬帶連續(xù)譜，與線譜相比，寬帶連續(xù)譜具有攜帶信息量大、背景噪聲相關(guān)性弱，有利于目標(biāo)檢測、參數(shù)估計(jì)及目標(biāo)特征提取等優(yōu)點(diǎn)，但當(dāng)不同方向目標(biāo)相關(guān)性較強(qiáng)、測量距離較大時(shí)，將導(dǎo)致方位估計(jì)算法性能下降。因此，開展水下遠(yuǎn)程寬帶相干源方位估計(jì)算法研究具有重要的工程應(yīng)用價(jià)值。

聲矢量傳感器較聲壓傳感器可獲取更多的聲場信息，為水下目標(biāo)探測提供了強(qiáng)有力的支撐。聲矢量傳感器陣列信號(hào)處理可有效提高水下目標(biāo)的探測能力，其中寬帶信號(hào)方位估計(jì)問題受到了廣泛關(guān)注。陳偉華等人(陳偉華，趙俊渭.聲矢量傳感器陣寬帶相干信號(hào)子空間最優(yōu)波束形成[j].聲學(xué)學(xué)報(bào)，2005，30(1):76-82.)將寬帶聚焦思想引入了聲矢量陣寬帶處理中，為聲矢量線陣寬帶聚焦矩陣的構(gòu)造提供了一種方法，提高了聲壓線陣的目標(biāo)方位估計(jì)性能。文獻(xiàn)“widebandmvdrbeamformingforacousticvectorsensorlineararray[j].ieeeproc.radarsonarnaving，2004，151(1):339-343.”采用寬帶相干信號(hào)子空間處理方法構(gòu)造聚焦矩陣，通過最小方差無畸變響應(yīng)(mvdr)算法實(shí)現(xiàn)了相干聲源的方位估計(jì)。但上述方法均將振速信息視為與聲壓相同的獨(dú)立陣元信息來處理，沒有充分利用聲矢量陣中聲壓與振速聯(lián)合處理的抗噪能力，在低信噪比條件下方位估計(jì)性能有限，無法滿足水下目標(biāo)遠(yuǎn)程探測的要求。國內(nèi)外學(xué)者h(yuǎn)awkes等(hawkesm,nehoraia.acousticvector-sensorcorrelationsinambientnoise[j].ieeejournalofoceanicengineering，2001，26(3):337-47.)以及惠俊英等人(惠俊英，劉宏，余華兵，范敏毅.聲壓振速聯(lián)合處理及其物理基礎(chǔ)初探[j].聲學(xué)學(xué)報(bào)，2000，25(4):303-307.)研究表明對(duì)于各向同性噪聲場，空間同點(diǎn)測得的聲壓振速不相關(guān)；對(duì)于有限尺度的信號(hào)源，空間同點(diǎn)測得的聲壓振速是完全相關(guān)的，從理論上證明了聲壓振速聯(lián)合處理的抗各向同性噪聲原理。在此基礎(chǔ)上，白興宇等人(白興宇，楊德森，趙春暉.基于聲壓振速聯(lián)合處理的聲矢量陣相干信號(hào)子空間方法[j].聲學(xué)學(xué)報(bào)，2006，31(5):410-417.)將相干信號(hào)子空間方法的寬帶高分辨能力與聲壓振速聯(lián)合處理的抗噪能力有機(jī)的結(jié)合起來，實(shí)現(xiàn)了對(duì)寬帶聲源的遠(yuǎn)程、高分辨方位估計(jì)。

與上述聲矢量線陣的豐碩研究成果相比，聲矢量圓陣的研究較少。均勻圓陣可獲360°全方位無模糊的方位角信息，且對(duì)于任意方位，均勻圓陣的分辨力幾乎相同，被廣泛應(yīng)用于雷達(dá)、聲納系統(tǒng)。但圓陣的陣列流形不具有vandermonde結(jié)構(gòu)，解相干算法不能直接應(yīng)用于圓陣，模式空間變換解決了該問題?；谀Ｊ娇臻g變換，wax等(waxm，sheinvaldj.directionfindingofcoherentsignalsviaspatialsmoothingforuniformcirculararrays[j].ieeetransactionsonantennas&propagation，1994，42(5):613-620.)引入空間平滑技術(shù)，實(shí)現(xiàn)了相干目標(biāo)的方位估計(jì)；張薇等人(張薇,韓勇.基于均勻圓陣的矢量重構(gòu)解相干算法[j].哈爾濱工業(yè)大學(xué)學(xué)報(bào).2016，48(5):62-66.)提出了一種模式空間矢量重構(gòu)算法，實(shí)現(xiàn)了相干目標(biāo)的方位估計(jì)。但上述研究方法可處理的信噪比門限高，不適合遠(yuǎn)程目標(biāo)的方位估計(jì)。文獻(xiàn)“circularacousticvector-sensorarrayformodebeamforming[j].ieeetransactiononsignalprocessing，2009，57(8):3041-3052.”中將模式空間變換應(yīng)用于聲矢量圓陣的振速分量，提高了方位估計(jì)性能，但未充分利用聲矢量陣中聲壓振速聯(lián)合處理的抗噪能力，低信噪比的探測能力有限。朱中銳等(楊德森，朱中銳，時(shí)勝國，莫世奇.聲矢量圓陣相位模態(tài)域目標(biāo)方位估計(jì)[j].聲學(xué)學(xué)報(bào)，2014，39(1):19-26.)利用聲場分解理論，將陣元域信號(hào)分解成為若干正交的相位模態(tài)，實(shí)現(xiàn)了窄帶遠(yuǎn)程目標(biāo)方位估計(jì)。針對(duì)寬帶方位估計(jì)問題，申請(qǐng)?zhí)枮?01610532431.1、名稱為“一種聲矢量圓陣寬帶相干源方位估計(jì)方法”的專利文件中，提出來了聲矢量圓陣寬帶聚焦變換方法，實(shí)現(xiàn)了遠(yuǎn)程寬帶目標(biāo)方位估計(jì)，但沒有引入解相干算法，在信號(hào)相干性較強(qiáng)的情況下不能有效估計(jì)目標(biāo)方位。

技術(shù)實(shí)現(xiàn)要素：

本發(fā)明的目的在于提供一種能提高低信噪比條件下寬帶強(qiáng)相干目標(biāo)的探測性能，可實(shí)現(xiàn)水下目標(biāo)遠(yuǎn)程探測的基于矢量奇異值分解的聲矢量圓陣寬帶相干源方位估計(jì)方法。

本發(fā)明的目的是這樣實(shí)現(xiàn)的：

(1)利用子帶分解，得到各個(gè)子帶的聲壓通道的頻域信號(hào)pe(fj)、徑向振速通道的頻域信號(hào)ver(fj)、切向振速通道的頻域信號(hào)

(2)構(gòu)建聲壓的模式空間變換矩陣tp(fj)、徑向振速的模式空間變換矩陣tvr(fj)、切向振速的模式空間變換矩陣將頻域信號(hào)變換為模態(tài)域聲壓通道的頻域信號(hào)p(fj)、模態(tài)域徑向振速通道的頻域信號(hào)vr(fj)、模態(tài)域切向振速通道的頻域信號(hào)

(3)采用聲壓p與振速聯(lián)合處理方式構(gòu)建每個(gè)子帶的互協(xié)方差矩陣，通過求和平均實(shí)現(xiàn)聲矢量圓陣寬帶接收信號(hào)的互協(xié)方差矩陣估計(jì)，vr為徑向振速、為切向振速；

(4)采用修正矢量奇異值方法，即將步驟(3)中所述的互協(xié)方差矩陣進(jìn)行特征值分解，尋找最大特征值對(duì)應(yīng)的特征矢量e1；并利用所述特征矢量e1構(gòu)建重構(gòu)矩陣y；然后，增加反向平滑項(xiàng)jd×dy^*jc×c得到修正重構(gòu)矩陣r；r為非方陣，對(duì)修正重構(gòu)矩陣r進(jìn)行奇異值分解，獲得修正重構(gòu)矩陣的信號(hào)子空間和噪聲子空間，jd×d、jc×c分別為d×d、c×c的反對(duì)角矩陣；

(5)采用music算法得到空間譜，根據(jù)譜峰位置輸出相干目標(biāo)方位。

本發(fā)明還可以包括：

1、步驟(2)具體包括：

所述聲壓的模式空間變換矩陣tp(fj)、徑向振速的模式空間變換矩陣tvr(fj)、切向振速的模式空間變換矩陣表示為：

其中，jp(fj)＝diag[i^-njn(kjr),…,iⁿjn(kjr),…,iⁿjn(kjr)]是聲壓通道對(duì)應(yīng)的由貝塞爾函數(shù)組成的矩陣，是振速通道對(duì)應(yīng)的由貝塞爾函數(shù)組成的矩陣，是振速通道對(duì)應(yīng)的由貝塞爾函數(shù)組成的矩陣，jn(kjr)為第一類n階貝塞爾函數(shù)，j'n(kjr)為jn(kjr)對(duì)r的導(dǎo)數(shù)，diag[·]表示對(duì)角矩陣，ρ為流體密度，kj＝2πfj/c表示波數(shù)，c表示聲速，kmax＝2πfmax/c為信號(hào)上限頻率fmax對(duì)應(yīng)的波數(shù)，-n≤n≤n，為圓陣可激發(fā)的最大相位模態(tài)數(shù)、為向下取整符號(hào)，i表示虛數(shù)單位；f＝[w-n,w-n+1,…wn,…,wn]第n列計(jì)算公式為wn＝[1,e^i2πn/m,…,e^i2πn(m-1)/m]^t；(·)^h表示(·)的共軛轉(zhuǎn)置運(yùn)算，(·)^-1為求逆運(yùn)算，m為陣元個(gè)數(shù)；

所述模態(tài)域聲壓通道的頻域信號(hào)p(fj)、模態(tài)域徑向振速通道的頻域信號(hào)vr(fj)、模態(tài)域切向振速通道的頻域信號(hào)如下：

2、步驟(3)具體包括：

(1)子帶的互協(xié)方差矩陣為：

(2)聲矢量圓陣寬帶接收信號(hào)的互協(xié)方差矩陣估計(jì)為：

3、步驟(4)具體包括：

(1)所述重構(gòu)矩陣y為：

其中，en為的特征矢量e1的第n個(gè)元素，c＝2k+1-d+1,d＞h,c＞h；

(2)所述修正重構(gòu)矩陣r為：

r＝]y,jd×dy^*jc×c]

(·)^*為(·)的共軛。

本發(fā)明將聲壓振速聯(lián)合處理、解相干算法及子空間算法進(jìn)行有效結(jié)合，比同陣型的聲壓陣、質(zhì)點(diǎn)振速獨(dú)立處理方法及p*vc聲壓振速聯(lián)合處理方式具有更強(qiáng)的抑制噪聲能力，且在低信噪比情況下該方法比未修正的矢量奇異值分解方法和前后向空間平滑具有更強(qiáng)的空間分辨能力。解決了相干源分辨困難、可處理的門限較高等問題。

本發(fā)明的有益效果是：針對(duì)聲矢量圓陣寬帶相干源方位估計(jì)問題，將聲壓振速聯(lián)合處理、解相干算法及子空間算法進(jìn)行有效結(jié)合，比同陣型的聲壓陣、質(zhì)點(diǎn)振速獨(dú)立處理方法及p*vc聲壓振速聯(lián)合處理方法具有更強(qiáng)的抑制噪聲能力，在低信噪比情況下該方法比矢量奇異值方法和前后向空間平滑具有更強(qiáng)的空間分辨能力。解決了相干源分辨困難、可處理的門限較高等問題，為聲矢量圓陣寬帶相干目標(biāo)的遠(yuǎn)程被動(dòng)探測問題提供了有效方法，具有良好的應(yīng)用前景。

附圖說明

圖1基于奇異值分解的聲矢量圓陣寬帶相干源方位估計(jì)方法流程圖。

圖2聲矢量圓陣布放示意圖。

圖3聲壓、矢量獨(dú)立、p*vc以及兩種聲壓振速聯(lián)合處理方式的空間譜對(duì)比分析。

圖4修正的矢量奇異值分解與矢量奇異值分解算法空間譜對(duì)比分析。

圖5不同算法的相干目標(biāo)方位估計(jì)性能對(duì)比分析。

圖6聲壓、矢量獨(dú)立、p*vc以及兩種聲壓振速聯(lián)合處理方式的試驗(yàn)結(jié)果。

圖7修正的矢量奇異值分解與矢量奇異值分解方法的試驗(yàn)結(jié)果。

圖8不同算法的相干目標(biāo)方位估計(jì)試驗(yàn)結(jié)果。

具體實(shí)施方式

下面結(jié)合附圖舉例對(duì)本發(fā)明作進(jìn)一步說明。

結(jié)合圖1，本發(fā)明的處理流程具體包括以下步驟：

第一步，基于子帶分解原理，得到各個(gè)子帶的頻域信號(hào)pe(fj)、ver(fj)、(1≤j≤j)。

結(jié)合圖2，陣元數(shù)為m的聲矢量圓陣位于xoy平面內(nèi)，半徑為r，聲矢量傳感器的振速分量x、y通道正軸方向分別沿著圓陣的徑向和切向，假設(shè)h個(gè)相干寬帶信號(hào)s1(t),s2(t),…,sh(t)分別從θ1,θ2，…,θh入射至聲矢量圓陣，信號(hào)與噪聲相互獨(dú)立。將接收到的寬帶信號(hào)進(jìn)行j點(diǎn)離散傅里葉變換，頻率fj的快拍數(shù)據(jù)可寫為：

式中，j＝1,2,…j；s(fj)＝[s1(fj),…,sh(fj)]^t為h個(gè)寬帶信號(hào)s1(t),s2(t),…,sh(t)對(duì)應(yīng)的頻率fj的信號(hào)向量，(·)^τ表示(·)的轉(zhuǎn)置運(yùn)算；pe(fj)、ver(fj)、分別為聲壓通道、徑向振速通道及切向振速通道的頻域信號(hào)；nep(fj)、ner(fj)、為聲壓通道、徑向振速通道和切向振速通道頻域的噪聲向量；aep(fj,θ)、aer(fj,θ)、分別為聲壓通道、徑向振速通道和切向振速通道對(duì)應(yīng)的陣列流型矩陣，具體形式為：

其中，θh表示第h個(gè)入射信號(hào)sh(t)的入射角度；aep(fj,θh)、aer(fj,θh)、表示對(duì)應(yīng)于入射角度θh的聲壓通道、徑向通道及切向通道的導(dǎo)向矢量，具體表達(dá)式如下：

其中，聲壓通道導(dǎo)向矢量中的第m個(gè)元素aepm(fj,θh)、徑向通道導(dǎo)向矢量中的第m個(gè)元素aerm(fj,θh)、切向通道導(dǎo)向矢量中的第m個(gè)元素的計(jì)算方法如下：

式中，i為虛數(shù)單位，i²＝-1，kj＝2πfj/c表示波數(shù)，c為聲速，為第m陣元與x正軸的夾角，cos(·)表示余弦運(yùn)算符號(hào)，sin(·)表示正弦運(yùn)算符號(hào)。

第二步，構(gòu)建模式變換矩陣tp(fj)、tvr(fj)、將頻域信號(hào)變換為模態(tài)域聲壓、振速通道的頻域信號(hào)p(fj)、vr(fj)、

(1)構(gòu)造聲壓pe(fj)、徑向振速ver(fj)和切向振速的模式空間變換矩陣tp(fj)、tvr(fj)、

式中，jp(fj)＝diag[i^-njn(kjr),…,iⁿjn(kjr),…,iⁿjn(kjr)]是聲壓通道對(duì)應(yīng)的由貝塞爾函數(shù)組成的矩陣，是振速通道對(duì)應(yīng)的由貝塞爾函數(shù)組成的矩陣，是振速通道對(duì)應(yīng)的由貝塞爾函數(shù)組成的矩陣，jn(kjr)為第一類n階貝塞爾函數(shù)，j'n(kjr)為jn(kjr)對(duì)r的導(dǎo)數(shù)，diag[·]表示對(duì)角矩陣，ρ為流體密度，-n≤n≤n，為圓陣可激發(fā)的最大相位模態(tài)數(shù)、為向下取整符號(hào)，kmax＝2πfmax/c為信號(hào)上限頻率fmax對(duì)應(yīng)的波數(shù)；f＝[w-n,w-n+1,…wn,…,wn]其中第n列計(jì)算公式為wn＝[1,e^i2πn/m,…,e^i2πn(m-1)/m]^t；(·)^h表示(·)的共軛轉(zhuǎn)置運(yùn)算，(·)^-1為求逆運(yùn)算，m為陣元個(gè)數(shù)。

模態(tài)域聲壓、徑向振速及切向振速通道的頻域信號(hào)為：

式中，表示模態(tài)域的陣列流型矩陣。

第三步，采用聲壓p與振速聯(lián)合處理方式構(gòu)建每個(gè)子帶的互協(xié)方差矩陣，通過求和平均實(shí)現(xiàn)聲矢量圓陣寬帶接收信號(hào)的互協(xié)方差矩陣估計(jì)，具體步驟如下：

(1)采用聲壓p與振速聯(lián)合處理方式構(gòu)建各個(gè)子帶的互協(xié)方差矩陣：

(2)通過求和平均實(shí)現(xiàn)聲矢量圓陣寬帶接收信號(hào)的互協(xié)方差矩陣估計(jì)：

第四步，采用修正的矢量奇異值方法，將第三步中的互協(xié)方差矩陣特征值分解得到最大特征值對(duì)應(yīng)的特征矢量e1；并利用特征矢量e1構(gòu)建重構(gòu)矩陣y；然后，增加反向平滑項(xiàng)jd×dy^*jc×c得到修正重構(gòu)矩陣r；對(duì)矩陣r(非方陣)進(jìn)行奇異值分解，獲得修正重構(gòu)矩陣的信號(hào)子空間和噪聲子空間，具體步驟如下：

(1)將第三步中的聲壓振速聯(lián)合處理的互協(xié)方差矩陣進(jìn)行特征值分解

(2)尋找最大特征值對(duì)應(yīng)的特征矢量e1。

(3)將特征矢量e1按下面規(guī)律排列得到重構(gòu)矩陣，方法如下：

其中，en為特征矢量e1的第n個(gè)元素，c＝2k+1-d+1,d＞h,c＞h。

(4)為了提高分辨能力，對(duì)矩陣y增加反向平滑項(xiàng)jd×dy^*jc×c，得到修正的重構(gòu)矩陣r(r為d×2c維矩陣)：

r＝[y,jd×dy^*jc×c](10)式中，jd×d、jc×c分別為d×d、c×c的反對(duì)角矩陣，(·)^*為(·)的共軛。

(5)對(duì)矩陣r(非方陣)進(jìn)行奇異值分解。

(6)獲得修正的重構(gòu)信號(hào)子空間us和噪聲子空間un。

第五步，采用music算得到空間譜，根據(jù)譜峰位置確定寬帶相干目標(biāo)方位：

式中，pmusic表示空間譜，θ為掃描角度，為對(duì)應(yīng)重構(gòu)矩陣的導(dǎo)向矢量。

上面對(duì)發(fā)明內(nèi)容各部分的具體實(shí)施方式進(jìn)行了說明。下面通過仿真實(shí)例和試驗(yàn)實(shí)例對(duì)本發(fā)明具體實(shí)施做進(jìn)一步描述。

仿真實(shí)例：

仿真采用8陣元均勻布放的聲矢量圓陣，圓陣半徑為0.35m；兩相干聲源的信號(hào)分別為s1(t),s2(t)，其中s1(t)為500-2500hz的寬帶信號(hào)，且s2(t)＝αs1(t)，α為不隨時(shí)間變化的常數(shù)。兩相干信號(hào)分別從θ1＝100°、θ2＝100°入射至聲矢量圓陣，采樣頻率為51.2×2.56khz，積分時(shí)間為t＝1s，寬帶時(shí)域信號(hào)分為128段，每段重疊點(diǎn)數(shù)為0，每段作4096點(diǎn)離散傅里葉變換，將工作頻帶分為j＝63份。

圖3為聲壓、矢量獨(dú)立、p*vc以及兩種聲壓振速聯(lián)合處理方式的空間譜結(jié)果。對(duì)比分析可知，聲壓圓陣在-5db時(shí)已不能分辨兩相干源的方位；矢量獨(dú)立處理方法在-5db時(shí)可估計(jì)相干雙聲源位置，然而該處理方法沒有利用聲壓振速聯(lián)合處理的抗噪性能，低信噪比條件下估計(jì)性能差；p*vc和充分利用了聲壓振速聯(lián)合處理的抗噪能力，較聲壓和矢量獨(dú)立處理方法有較好的抑制噪聲能力及較高的空間分辨能力，而且隨著信噪比的降低，聯(lián)合處理方法較p*vc聯(lián)合處理方法具有更高的噪聲抑制能力及較高的角分辨能力。

圖4為修正的矢量奇異值分解與矢量奇異值分解空間譜對(duì)比分析。其中，svd表示矢量奇異值分解的空間譜，msvd表示修正的矢量奇異值分解的空間譜。對(duì)比分析可知，在信噪比較高的條件下，svd及msvd均能分辨兩相干源，當(dāng)信噪比降低至-20db時(shí)，svd方法已不能有效估計(jì)相干目標(biāo)方位，而msvd方法仍能有效估計(jì)兩相干源方位。可見，msvd方法的相干源方位估計(jì)性能優(yōu)于svd方法。

圖5為不同算法的相干目標(biāo)方位估計(jì)性能對(duì)比分析。其中，msvd為聯(lián)合處理的修正矢量奇異值方法空間譜，fbss為聯(lián)合處理的前后向空間平滑方法空間譜。對(duì)比分析可知，當(dāng)snr＝-5db時(shí)，fbss方法及msvd方法均可有效估計(jì)相干源方位，但信噪比降低，fbss方法空間增益降低，當(dāng)snr＝-20db時(shí)，fbss方法已失效，而msvd方法仍可有效分辨兩相干源方位。

試驗(yàn)實(shí)例：

在消聲水池開展聲矢量圓陣寬帶相干目標(biāo)方位估計(jì)實(shí)驗(yàn)測試研究，驗(yàn)證本算法的可行性及有效性。實(shí)驗(yàn)采用8元均勻布放的生矢量圓陣，基陣半徑為0.35m，聲矢量傳感器的振速x、y分別與圓陣的徑向、切向重合。實(shí)驗(yàn)的相干目標(biāo)信號(hào)是由發(fā)射換能器發(fā)射的頻率為600-2400hz寬帶信號(hào)，即由信號(hào)源產(chǎn)生的寬帶白噪聲輸出給帶通濾波器，并將帶通濾波器的輸出信號(hào)分別通過兩套發(fā)射系統(tǒng)(包括功率放大器、發(fā)射換能器)發(fā)射的寬帶信號(hào)。兩只發(fā)射換能器相對(duì)于聲矢量圓陣的入射角度分別約為198°，260°。實(shí)驗(yàn)采集記錄了背景噪聲和目標(biāo)信號(hào)，通過調(diào)整背景噪聲大小并將其混入測量的目標(biāo)信號(hào)中，獲得不同信噪比條件下聲矢量圓陣接收數(shù)據(jù)。定義信噪比snr為：

式中，分別為處理帶寬內(nèi)信號(hào)功率及噪聲功率。

圖6為聲壓、矢量獨(dú)立、p*vc以及兩種聲壓振速聯(lián)合處理方式的試驗(yàn)結(jié)果。由圖可知，在幾種處理方式中聲壓振速聯(lián)合處理方法比其他方式具有更高的空間分辨力及抑制噪聲能力，實(shí)現(xiàn)了寬帶相干目標(biāo)的遠(yuǎn)程探測。

圖7為為修正的矢量奇異值分解與矢量奇異值分解空間譜對(duì)比分析。由圖可知修正的矢量奇異值方法較矢量奇異值方法的相干源分辨能力較高，抑制噪聲能力較強(qiáng)。

圖8為不同算法的相干目標(biāo)方位估計(jì)性能對(duì)比分析。由圖可知，矢量奇異值方法比前后向空間平滑方法具有更高的空間分辨力及抗噪能力，更適合寬帶相干目標(biāo)的遠(yuǎn)程探測。