本發(fā)明涉及一種波達(dá)估計(jì)方法，尤其涉及一種基于高階累計(jì)量的波達(dá)快速估計(jì)方法，屬于信號(hào)處理技術(shù)領(lǐng)域。

背景技術(shù)：

在淺水中，聲波的傳播過(guò)程總是伴隨著在海洋表面或水中暗礁處的反射或折射現(xiàn)象。因此聲波是在多路徑模型中傳播的。多路徑傳播特性在淺水中有廣泛地應(yīng)用，例如主動(dòng)聲納，海洋聲層析術(shù)與聲源定位等。多路徑傳播可以覆蓋更多的海洋環(huán)境，但同時(shí)產(chǎn)生的相干信號(hào)會(huì)對(duì)處理造成干擾。因此需要對(duì)接收信號(hào)進(jìn)行分離處理。

這個(gè)問(wèn)題吸引了相關(guān)研究者極大的關(guān)注并提出了很多有價(jià)值的算法。在這些算法中，基于子空間類算法由于高分辨能力而受到研究者的重視。特別是由Schmidt[Schmidt R O.Multiple emitter location and signal parameter estimation[J].Antennas and Propagation,IEEE Transactions on,1986,34(3):276-280.]提出的多重信號(hào)分離算法(Multiple Signal Classification:MUSIC)，這是十分典型的基于子空間方法。該算法的基本思想是通過(guò)對(duì)接收數(shù)據(jù)的協(xié)方差矩陣進(jìn)行特征分解，構(gòu)成與代表信號(hào)部分的信號(hào)子空間正交的噪聲子空間。該算法可以用來(lái)分離非相干的窄帶信號(hào)線路徑。但是在多路徑傳播模型中由于接收的線路徑是發(fā)射信號(hào)的不同時(shí)延版本，因此它們是高度相關(guān)甚至相干的。

為了可以有效地分離高度相關(guān)甚至相干的寬帶信號(hào)，Jiang等人在“Raypath separation with high resolution processing”一文[Jiang L,Aulanier F,Le TouzéG,et al.Raypath separation with high resolution processing,2011]中提出一種高分辨率方法：平滑的主動(dòng)寬帶信號(hào)分離算法(smoothing-MUltiple SIgnal Classification Active Large band：s-MUSICAL)。它結(jié)合了主動(dòng)寬帶算法與空域-頻域平滑算法，并且可以同時(shí)在到達(dá)角-到達(dá)時(shí)間域內(nèi)分離信號(hào)。實(shí)驗(yàn)表明s-MUSICAL算法極大地改善了信號(hào)分離表現(xiàn)，尤其是到達(dá)角度相差極小的兩個(gè)線路徑。s-MUSICAL 算法建立在環(huán)境中僅存在高斯白噪聲的假設(shè)之上，但在實(shí)際海洋環(huán)境中存在譜未知的高斯有色噪聲。

近期，Jiang提出了一種高階主動(dòng)寬帶MUSIC算法(參見公開日為2016/6/8的中國(guó)專利CN 105652264A)，可以分離受高斯色噪聲污染的線路徑。它利用四階累積量，表達(dá)出非高斯過(guò)程的有價(jià)值信息。該方法主要步驟為：步驟A、對(duì)觀測(cè)數(shù)據(jù)做傅里葉變換后進(jìn)行空域-頻域平滑處理，并構(gòu)造出空域-頻域平滑處理后信號(hào)的四階累積量矩陣步驟B、根據(jù)所述四階累積量矩陣?yán)闷娈愔捣纸獾姆椒?gòu)建觀測(cè)數(shù)據(jù)的信號(hào)子空間；步驟C、根據(jù)觀測(cè)數(shù)據(jù)的信號(hào)子空間，估計(jì)出信號(hào)源的波達(dá)方向及波達(dá)時(shí)間。實(shí)驗(yàn)表明該算法可以更有效地抑制噪聲進(jìn)而獲得更準(zhǔn)確的結(jié)果。但是在處理過(guò)程中需要對(duì)四階累積量矩陣進(jìn)行奇異值分解。奇異值分解(Singular Value Decomposition:SVD)處理步驟需要大量的計(jì)算成本，包括運(yùn)算時(shí)間與內(nèi)存空間。對(duì)矩陣進(jìn)行SVD處理過(guò)程需要的時(shí)間和空間與矩陣規(guī)模成超線性關(guān)系。這無(wú)法滿足實(shí)際應(yīng)用中的實(shí)時(shí)性要求。唐建紅等人在“改進(jìn)的基于四階累積量的MUSIC算法”(唐建紅、司錫才與初萍,改進(jìn)的基于四階累積量的MUSIC算法.系統(tǒng)工程與電子技術(shù),2010.32(2):第256-259頁(yè))一文中提出一種改進(jìn)方法，朱敏等人在“一種新的基于四階累積量的DOA估計(jì)算法”(朱敏與何培宇,一種新的基于四階累積量的DOA估計(jì)算法.四川大學(xué)學(xué)報(bào)(自然科學(xué)版),2011.48(2):第343-348頁(yè).)一文中也提出一種類似的改進(jìn)方法，均利用四階累積量矩陣自身的結(jié)構(gòu)特性，根據(jù)接收端為均勻線陣的情況，通過(guò)觀察矩陣內(nèi)部的冗余進(jìn)而縮小矩陣規(guī)模以降低運(yùn)算量。但這樣的方法與陣列情況等因素有關(guān)，不同陣列類型的情況下有失效的可能。因此需要一種適用性較強(qiáng)的改進(jìn)策略，增強(qiáng)基于高階累積量的信號(hào)分離算法的實(shí)時(shí)性。

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

本發(fā)明所要解決的技術(shù)問(wèn)題在于克服現(xiàn)有技術(shù)不足，提供一種基于高階累計(jì)量的波達(dá)快速估計(jì)方法，在保證較高分辨能力的同時(shí)，可大幅降低現(xiàn)有技術(shù)的時(shí)間復(fù)雜度與空間復(fù)雜度，減少硬件資源的消耗并提高信號(hào)處理的實(shí)時(shí)性。

本發(fā)明具體采用以下技術(shù)方案解決上述技術(shù)問(wèn)題：

一種基于高階累計(jì)量的波達(dá)快速估計(jì)方法，利用等間距直線傳感器陣列所接 收到的觀測(cè)數(shù)據(jù)，估計(jì)出信號(hào)源的波達(dá)方向及波達(dá)時(shí)間；該方法包括以下步驟：

步驟A、對(duì)所述觀測(cè)數(shù)據(jù)做傅里葉變換后進(jìn)行空域-頻域平滑處理，并構(gòu)造出空域-頻域平滑處理后信號(hào)的四階累積量矩陣

步驟B、根據(jù)所述四階累積量矩陣構(gòu)建觀測(cè)數(shù)據(jù)的信號(hào)子空間；

步驟C、根據(jù)觀測(cè)數(shù)據(jù)的信號(hào)子空間，估計(jì)出信號(hào)源的波達(dá)方向及波達(dá)時(shí)間；

所述步驟B具體如下：

步驟B1、對(duì)矩陣的全部n個(gè)列向量分別按照概率進(jìn)行隨機(jī)抽樣，得到c個(gè)列向量，其中表示矩陣的第i列，||||_F表示Frobenius-范數(shù)，具體抽取過(guò)程：產(chǎn)生0～1的隨機(jī)數(shù)，若該隨機(jī)數(shù)小于前i列概率之和且大于前i-1列概率之和，則第i列被抽中，按此方法進(jìn)行c次，抽取出c個(gè)列向量，1<c<n；之后對(duì)抽取出的c個(gè)列向量進(jìn)行縮放調(diào)整，第t次抽樣抽取出的列向量的縮放因子為：將調(diào)整后的列向量組成矩陣

步驟B2、對(duì)矩陣進(jìn)行特征值分解；

步驟B3、利用矩陣較大的前k個(gè)特征值對(duì)應(yīng)的特征向量y，計(jì)算得到矩陣 的k個(gè)特征向量h；其中特征向量y與特征向量h之間的關(guān)系為： 其中，h^t表示第t個(gè)特征向量h，y^t表示第t個(gè)特征向量y，表示矩陣的第t個(gè)特征值，k為傳感器陣列所接收到的信號(hào)個(gè)數(shù)的平方；

步驟B4、以這k個(gè)特征向量h張成的空間近似作為觀測(cè)數(shù)據(jù)的信號(hào)子空間。

上述波達(dá)快速估計(jì)方法可廣泛用于雷達(dá)、聲納、地震監(jiān)測(cè)等領(lǐng)域，以下為兩個(gè)具體應(yīng)用方案：

一種基于聲線傳播時(shí)間層析的海洋聲層析方法，利用聲音在海洋中傳播速度的變化來(lái)反演海洋環(huán)境參數(shù)，首先利用以上技術(shù)方案所述方法對(duì)從聲傳感器陣列所接收到的多路徑傳播聲信號(hào)進(jìn)行波達(dá)估計(jì)，從而分離出每一條聲線路徑；然后根據(jù)聲線路徑的到達(dá)時(shí)間反演出海洋環(huán)境參數(shù)。

一種定位方法，首先利用以上技術(shù)方案所述方法進(jìn)行波達(dá)方向估計(jì)，然后利用估計(jì)出的波達(dá)方向確定信號(hào)源的位置。

相比現(xiàn)有技術(shù)，本發(fā)明具有以下有益效果：

和基于二階統(tǒng)計(jì)量的傳統(tǒng)技術(shù)相比，本發(fā)明對(duì)實(shí)際應(yīng)用場(chǎng)景中成分復(fù)雜的環(huán)境噪聲有明顯的抑制效果，同時(shí)可以在傳感器陣列陣元數(shù)目小于入射線路徑數(shù)目的情況下正確分離線路徑。更重要的是，本發(fā)明與現(xiàn)有基于四階累積量的方法相比，僅需占用極少的計(jì)算時(shí)間與內(nèi)存空間，為基于高階統(tǒng)計(jì)量原理的信號(hào)分離技術(shù)應(yīng)用于實(shí)時(shí)系統(tǒng)提供了可行性。

附圖說(shuō)明

圖1為本發(fā)明方法的流程示意圖；

圖2a～圖2c為本發(fā)明方法與現(xiàn)有方法的對(duì)比實(shí)驗(yàn)結(jié)果。

具體實(shí)施方式

下面結(jié)合附圖對(duì)本發(fā)明的技術(shù)方案進(jìn)行詳細(xì)說(shuō)明：

本發(fā)明的發(fā)明思路是針對(duì)現(xiàn)有基于高階累積量的陣列信號(hào)處理方法所存在的計(jì)算復(fù)雜度過(guò)高的問(wèn)題，利用線性時(shí)間的奇異值分解方法進(jìn)行信號(hào)子空間的近似，從而大幅縮短算法的執(zhí)行時(shí)間并大幅節(jié)約算法的占用內(nèi)存。

本發(fā)明基于高階累計(jì)量的波達(dá)快速估計(jì)方法的流程如圖1所示，包括以下步驟：

步驟A、對(duì)所述觀測(cè)數(shù)據(jù)做傅里葉變換后進(jìn)行空域-頻域平滑處理，并構(gòu)造出空域-頻域平滑處理后信號(hào)的四階累積量矩陣

首先建立信號(hào)模型：

考慮P個(gè)線路徑被由M個(gè)傳感器組成的均勻垂直陣列接收。則第m個(gè)傳感器接收到的信號(hào)在時(shí)域內(nèi)可以如下建模：

其中e(t)是發(fā)射信號(hào)，a_p表示第p個(gè)線路徑的幅值，τ_m,p是第p個(gè)線路徑的 傳播時(shí)延，n_m(t)是第m個(gè)傳感器接收到的加性噪聲。另外，第p個(gè)線路徑的傳播時(shí)延τ_m,p可以表示為：

τ_m,p＝T_p+t_m(θ_p) (2)

其中T_p代表第p個(gè)線路徑到達(dá)參考傳感器的傳播時(shí)間，t_m(θ_p)代表參考傳感器與第m個(gè)傳感器之間的時(shí)延。t_m(θ_p)是第p個(gè)線路徑到達(dá)接收陣列的角度θ_p的函數(shù)。

在頻域內(nèi)，(1)式可以表示為：

將(2)式與(3)式結(jié)合可得到：

其中：

Ψ_p＝2πT_p Φ_p＝2πt_m(θ_p)

在(3)式中，e(v)項(xiàng)是確定的發(fā)射信號(hào)在頻率v上的幅值；a_p項(xiàng)是每個(gè)線路徑的幅值，可以看做隨機(jī)的，不相干的。

假定信號(hào)被分為F個(gè)頻點(diǎn)，因此(3)是可以被寫成矩陣形式：

x_g＝H.A+n_g (5)

其中：

·x_g＝[x⁺(v₁)，x⁺(v₂)，...，x⁺(v_F)]⁺是一個(gè)(M×F)維的長(zhǎng)向量，包含所有傳感器的所有頻點(diǎn)信息，其中x(v_i)＝[x₁(v_i)，x₂(v_i)，...，x_M(v_i)]⁺

·n_g＝[n⁺(v₁)，n⁺(v₂)，...，n⁺(v_F)]⁺是一個(gè)(M×F)維的長(zhǎng)向量，是每個(gè)頻點(diǎn)上接收到的噪聲向量的串聯(lián)，其中n(v_i)＝[n₁(v_i)，n₂(v_i)，...，n_M(v_i)]⁺

·H＝[h₁，h₂，...，h_P]⁺是一個(gè)(M×F，P)維矩陣，將表征發(fā)射端與接收端之間傳播函數(shù)的項(xiàng)整合到一起，其中e(v_i)表示發(fā)射信號(hào)

·A＝[a₁，a₂，...，a_P]⁺是一個(gè)P維向量

·上標(biāo)+表示轉(zhuǎn)置

因此接收數(shù)據(jù)的四階累積量矩陣可以表示為：

·表示克羅內(nèi)克積

·上標(biāo)*表示共軛

·H表示共軛轉(zhuǎn)置

鑒于到達(dá)接收陣列的線路徑是發(fā)射信號(hào)不同時(shí)延或振幅增益的副本，它們是高度相關(guān)甚至相干的。因此矩陣C的秩虧損。通過(guò)空間-頻域平滑預(yù)處理將矩陣C的秩恢復(fù)為P²。所以估計(jì)的平滑后的四階累積量矩陣為：

則此時(shí)矩陣的秩最大可為K＝(2k_s+1)(2k_f+1)，k_s、k_f分別為空間平滑與頻域平滑的階。因此為了分離P個(gè)相干的線路徑，需保證K≥P²。

步驟B、根據(jù)所述四階累積量矩陣構(gòu)建觀測(cè)數(shù)據(jù)的信號(hào)子空間；

接下來(lái)需要對(duì)四階累積量矩陣進(jìn)行特征分解，進(jìn)而構(gòu)造出信號(hào)子空間與噪聲子空間。由于現(xiàn)有基于高階統(tǒng)計(jì)量的信號(hào)分離方法中的特征分解步驟需要耗費(fèi)大量的計(jì)算資源，因此，本發(fā)明采用隨機(jī)SVD方法，使分解過(guò)程所消耗的時(shí)間與占用的內(nèi)存空間極大地減少；具體地，

首先對(duì)矩陣的列向量以一定的概率分布進(jìn)行隨機(jī)抽取，即產(chǎn)生0～1的隨機(jī)數(shù)，若該隨機(jī)數(shù)小于前i列概率之和且大于前i+1列概率之和，則第i列被抽中，按此方法進(jìn)行c次，抽取出c個(gè)列向量。概率為即第a列被選中的概率為：Q[i_t＝a]＝q_a,a＝1,…,n。

為了使估計(jì)更準(zhǔn)確，將抽選出的列向量進(jìn)行放縮并構(gòu)成矩陣縮放因子為： 即

對(duì)矩陣進(jìn)行SVD操作，令

通過(guò)矩陣較大的P²個(gè)特征值對(duì)應(yīng)的特征向量，利用計(jì)算得到矩陣較大的P²個(gè)特征值對(duì)應(yīng)的特征向量，組成矩陣即則矩陣的列空間近似為信號(hào)子空間。

對(duì)于抽樣個(gè)數(shù)c，可取c≥4η²/∈²或c≥4P²η²/ε²。其中若使抽樣個(gè)數(shù)滿足c≥4η²/∈²，則利用衡量近似信號(hào)子空間與信號(hào)子空間之間的誤差時(shí)，該誤差D₂滿足的概率為1-δ；若使抽樣個(gè)數(shù)滿足c≥4P²η²/ε²，則利用衡量近似信號(hào)子空間與信號(hào)子空間之間的誤差，該誤差D_F滿足的概率為1-δ。||||_F表示Frobenius-范數(shù)，表示矩陣的最優(yōu)P²階近似矩陣，是一個(gè)與δ相關(guān)的量，表示矩陣的第i列。

步驟C、根據(jù)觀測(cè)數(shù)據(jù)的信號(hào)子空間，估計(jì)出信號(hào)源的波達(dá)方向及波達(dá)時(shí)間:

利用陣列流型張成空間與信號(hào)子空間的重合性，構(gòu)造估計(jì)函數(shù)：

θ表示信號(hào)路徑的波達(dá)角度；T表示信號(hào)路徑的波達(dá)時(shí)間；e(v_i)表示信號(hào)在頻率v_i處的幅值，i＝1,2,..,F；τ_1,j(θ)表示信號(hào)路徑到達(dá)第j個(gè)傳感器相對(duì)于到達(dá)作為參考傳感器的時(shí)間延遲，j＝2,3,…,M-1。

估計(jì)函數(shù)的最大值所對(duì)應(yīng)的θ和T即信號(hào)源的波達(dá)方向和波達(dá)時(shí)間。

為了驗(yàn)證本發(fā)明的效果，將其應(yīng)用于仿真實(shí)驗(yàn)數(shù)據(jù)中，并與傳統(tǒng)基于高階累積量的信號(hào)分離算法的運(yùn)行結(jié)果與運(yùn)行時(shí)情況進(jìn)行比較。實(shí)驗(yàn)使用的射線路徑的 相關(guān)信息為：四條射線路徑在各個(gè)傳感器間的延遲時(shí)間分別為：0s，2×10^-3s，-1×10^-3s，1×10^-3s；四條射線路徑到達(dá)傳感器的時(shí)間分別為5×10^-3s，15×10^-3s，15×10^-3s，20×10^-3s；接收端是由六個(gè)傳感器組成的均勻陣列，選取第一個(gè)傳感器為參考傳感器，傳感器采樣數(shù)據(jù)長(zhǎng)度為128；兩種方法所采用的空域和頻域的平滑階次均為1；添加高斯彩色噪聲，信噪比為20dB；失敗概率為0.01，期望誤差系數(shù)ε＝2。

算法運(yùn)行在Intel(R)Core(TM)i5-4590CPU@3.30GHz計(jì)算機(jī)上，操作系統(tǒng)為Windows 7Ultimate 64-bit。

實(shí)驗(yàn)中，四階累積量矩陣是一個(gè)規(guī)模為14400×14400的方陣，占用了3.16Gb

表1為該實(shí)驗(yàn)中分別將使用經(jīng)典的四階累積量信號(hào)分離方法與使用本發(fā)明方法獲得信號(hào)子空間過(guò)程的時(shí)間與空間消耗，定義M為L(zhǎng)TSVD方法與傳統(tǒng)方法時(shí)間或所占空間的比值。

表1

根據(jù)以上實(shí)驗(yàn)結(jié)果可知，與直接將四階累積量矩陣進(jìn)行特征分解相比，本發(fā)明方法在時(shí)間和空間方面均有顯著地性能提升：將原方法的執(zhí)行時(shí)間縮短至近千分之一，申請(qǐng)的內(nèi)存空間也縮小至不到百分之一。

圖2a～圖2c為本發(fā)明方法與方法結(jié)果的對(duì)比：圖2a、圖2b、圖2c依次為現(xiàn)有方法、F-范數(shù)約束的本發(fā)明方法和2-范數(shù)約束的本發(fā)明方法的信號(hào)分離結(jié)果?？梢?，本發(fā)明方法沒(méi)有影響原方法的準(zhǔn)確度。

本發(fā)明方法可廣泛用于雷達(dá)、聲納、地震監(jiān)測(cè)等領(lǐng)域，以下為兩個(gè)具體應(yīng)用方案：

一種基于聲線傳播時(shí)間層析的海洋聲層析方法，利用聲音在海洋中傳播速度的變化來(lái)反演海洋環(huán)境參數(shù)，首先利用以上技術(shù)方案所述方法對(duì)從聲傳感器陣列所接收到的多路徑傳播聲信號(hào)進(jìn)行波達(dá)估計(jì)，從而分離出每一條聲線路徑；然后根據(jù)聲線路徑的到達(dá)時(shí)間反演出海洋環(huán)境參數(shù)。

一種定位方法，首先利用以上技術(shù)方案所述方法進(jìn)行波達(dá)方向估計(jì)，然后利用估計(jì)出的波達(dá)方向確定信號(hào)源的位置。