一種基于高階累積量的波達(dá)估計(jì)方法
【專利摘要】本發(fā)明公開了一種基于高階累積量的波達(dá)估計(jì)方法，屬于信號(hào)處理技術(shù)領(lǐng)域。該方法利用等間距直線傳感器陣列所接收到的觀測(cè)信號(hào)，估計(jì)出信號(hào)源的波達(dá)方向及波達(dá)時(shí)間；包括以下步驟：步驟1、對(duì)觀測(cè)信號(hào)做傅里葉變換后進(jìn)行空域?頻域平滑處理；步驟2、構(gòu)造出空域?頻域平滑處理后信號(hào)的四階累積量矩陣；步驟3、利用迭代地部分SVD方法，根據(jù)四階累積量矩陣構(gòu)建觀測(cè)信號(hào)的信號(hào)子空間和噪聲子空間；步驟4、根據(jù)觀測(cè)信號(hào)的信號(hào)子空間和噪聲子空間之間的正交性，估計(jì)出信號(hào)源的波達(dá)方向及波達(dá)時(shí)間。本發(fā)明還公開了一種基于聲線傳播時(shí)間層析的海洋聲層析方法及一種定位方法。本發(fā)明可大幅降低算法的計(jì)算復(fù)雜度，提高算法實(shí)時(shí)性并降低硬件資源消耗。
【專利說明】
一種基于高階累積量的波達(dá)估計(jì)方法
技術(shù)領(lǐng)域
[0001] 本發(fā)明涉及信號(hào)處理技術(shù)領(lǐng)域，尤其涉及一種基于高階累積量的波達(dá)估計(jì)方法。
【背景技術(shù)】
[0002] 波達(dá)方向（D0A)估計(jì)指的是要確定同時(shí)處在空間某一區(qū)域內(nèi)多個(gè)感興趣信號(hào)的空 間位置，即各個(gè)信號(hào)到達(dá)陣列參考陣元的方向角。波達(dá)方向技術(shù)是陣列信號(hào)處理中的重要 研究課題。波達(dá)方向估計(jì)技術(shù)的關(guān)鍵在于利用處于空間不同位置的天線陣列，接收多個(gè)不 同方向的信號(hào)源發(fā)出的信號(hào)，運(yùn)用現(xiàn)代處理方法快速且高精度地估計(jì)出信號(hào)源的方向。波 達(dá)方向估計(jì)技術(shù)目前廣泛應(yīng)用于雷達(dá)、聲納、地震等領(lǐng)域。
[0003] 八十年代以后，學(xué)術(shù)界提出一類基于矩陣特征值分解的譜估計(jì)方法，以Schmidt等 人提出的MUSIC(Multipie Signal Classification)算法為代表，它是基于信號(hào)子空間和 噪聲子空間的正交性，具有很好的角度分辨能力。但是當(dāng)信號(hào)相關(guān)或者相干時(shí)，樣本的協(xié)方 差矩陣會(huì)不滿秩，MUSIC算法則不能給出正確的分離結(jié)果。為了解決這個(gè)問題，學(xué)者們給出 了許多方法。特別地，Jiang 等人在"Raypath separation with high resolution processing" 一文中，提出在點(diǎn)到陣列結(jié)構(gòu)中，結(jié)合空域-頻域平滑方法與主動(dòng)寬帶多信號(hào) 分離算法，給出平滑的主動(dòng)寬帶多信號(hào)分離算法smoothing-MUSICAL(smoothing Multiple Signal Classification Active Large-band)。該算法仍然存在著一些問題：（1)分辨能力 并不能滿足實(shí)際應(yīng)用的需求；（2)該算法基于二階統(tǒng)計(jì)量，需要假設(shè)信號(hào)為高斯信號(hào)；（3)傳 感器陣列的陣元數(shù)目必須大于待分離信號(hào)源的數(shù)目。
[0004] 為了解決這些問題，一篇中國(guó)發(fā)明專利申請(qǐng)(發(fā)明名稱為"基于高階累積量的多路 徑傳播聲信號(hào)分離方法"，申請(qǐng)?zhí)枮?01610006173.3,申請(qǐng)日為2016-1-5)提出了采用基于 高階累積量的陣列信號(hào)處理方法。該方法的具體步驟如下：
[0005] 步驟1、對(duì)傳感器所接收聲信號(hào)做傅里葉變換后進(jìn)行空域-頻域平滑處理，獲得K = (2Ks+l)(2Kf+l)個(gè)窄帶估計(jì)組成矩陣砂^，1(3、1^分別為空域平滑階次、頻域平滑階次;且1( 大于聲源個(gè)數(shù)P;
[0006] 步驟2、利用下式計(jì)算矩陣lks,kf的四階累積量& ^ - ￡ f few ^
[0007] 、
[0008] 其中，E表示求期望，*表示求共輒，H表示求共輒轉(zhuǎn)置，#表示克羅內(nèi)克積；
[0009] 步驟3、對(duì)四階累積量f進(jìn)行EVD特征分解，并將所得到的(MF)2個(gè)特征值從大到小 排列為…>2(#廠其中，F(xiàn)為頻域平滑時(shí)選取的頻率數(shù);M為聲傳感器陣列中的傳 感器數(shù)量；
[0010] 步驟4、對(duì)f進(jìn)行特征分解后得到(MF)2個(gè)特征值?為、…、心、…、，以
[0011 ]其中P2個(gè)較大的特征值斗…、為^對(duì)應(yīng)的特征向量構(gòu)造信號(hào)子空間，以剩下的(MF )2-p2個(gè)較小特征值士。、…、2(A/f卩對(duì)應(yīng)的特征向量構(gòu)造噪聲子空間；
[0012] 步驟5、為聲線路徑構(gòu)造以下的副本矢量a(0，T):
[0013] a{OJ')^h{0：r)?h{OJ')^ 玢(…2卿知§)
[0014] bUt} … e(外知'"淑
[0015] = [1,廠2咖內(nèi).辦),…2坎
[0016] 其中，0表示聲線路徑的到達(dá)角度，T表示聲線路徑的到達(dá)時(shí)間，e(Vi)是聲信號(hào)在 頻率Vl處的幅值，i = l，2,…，F(xiàn)，發(fā)表示克羅內(nèi)克積，Tld(0)表示聲線路徑到達(dá)第j個(gè)傳感器 相對(duì)于到達(dá)作為參考傳感器的第1個(gè)傳感器的時(shí)間延遲，j = 2,…，M-1;
[0017] 步驟6、搜索與所述噪聲子空間正交的副本矢量，這些副本矢量即為最終分離出的 聲線路徑，其0和T即為聲線路徑相應(yīng)的到達(dá)方向和到達(dá)時(shí)間。
[0018] 由于高階累積量對(duì)任意的高斯噪聲有自然盲性，基于高階累積量的算法使原有的 D0A估計(jì)所適應(yīng)的觀測(cè)噪聲擴(kuò)展到高斯空間有色噪聲或者對(duì)稱的非高斯空間有色噪聲和白 噪聲。
[0019] 然而，上述基于高階累積量的算法，由于引入了樣本的高階累積量矩陣，這將導(dǎo)致 矩陣規(guī)模極其龐大，隨之而來的問題是對(duì)更大的存儲(chǔ)空間和內(nèi)存的需求，同時(shí)基于特征值 分解的這樣一個(gè)計(jì)算過程的時(shí)間消耗也會(huì)指數(shù)級(jí)增長(zhǎng)。

【發(fā)明內(nèi)容】

[0020] 本發(fā)明所要解決的技術(shù)問題在于克服現(xiàn)有技術(shù)不足，提供一種基于高階累積量的 波達(dá)估計(jì)方法，一方面可獲得更高的信號(hào)估計(jì)精度并能適用于傳感器數(shù)目小于射線路徑的 情況，同時(shí)還可抑制高斯噪聲和非高斯噪聲的影響；另一方面可大幅降低算法的計(jì)算復(fù)雜 度，提高算法實(shí)時(shí)性并降低硬件資源消耗。
[0021] 本發(fā)明具體采用以下技術(shù)方案解決上述技術(shù)問題：
[0022] 一種基于高階累積量的波達(dá)估計(jì)方法，利用等間距直線傳感器陣列所接收到的觀 測(cè)信號(hào)，估計(jì)出信號(hào)源的波達(dá)方向及波達(dá)時(shí)間;包括以下步驟：
[0023] 步驟1、對(duì)所述觀測(cè)信號(hào)做傅里葉變換后進(jìn)行空域_頻域平滑處理，得到K=(2KS+ l)(2Kf+l)個(gè)窄帶估計(jì)作為列向量所組成的矩陣X，Ks、Kf分別為空域平滑階次、頻域平滑階 次;且K大于信號(hào)源個(gè)數(shù)N;
[0024]步驟2、根據(jù)下式構(gòu)造出空域-頻域平滑處理后信號(hào)的四階累積量矩陣C: C = ^{(X ? X* )(X <8) X* )H }- ^{x <g) X* }^{(X (E) X* }
[0025] -￡{XXH}0￡{XX^}
[0026] X表示空域-頻域平滑處理之后得到的矩陣，E表示期望，上標(biāo)*表示共輒，上標(biāo)H表 示共輒轉(zhuǎn)置，?表示克羅內(nèi)克積；
[0027] 步驟3、根據(jù)所述四階累積量矩陣構(gòu)建觀測(cè)信號(hào)的信號(hào)子空間和噪聲子空間；
[0028] 步驟4、根據(jù)觀測(cè)信號(hào)的信號(hào)子空間和噪聲子空間之間的正交性，估計(jì)出信號(hào)源的 波達(dá)方向及波達(dá)時(shí)間；
[0029]所述步驟3具體包括以下子步驟：
[0030]步驟301、對(duì)兩個(gè)規(guī)模均為(MF)2 X p的正交矩陣U、V進(jìn)行隨機(jī)初始化，記為lAV13，矩 陣右上角的數(shù)字標(biāo)號(hào)表示迭代的次數(shù)，M為所述等間距直線傳感器陣列的陣元個(gè)數(shù)，F(xiàn)為空 域-頻域平滑處理過程中選取的頻率通道數(shù)，P=K 2且p<<(MF)2;
[0031] 步驟302、對(duì)所述兩個(gè)正交矩陣U、V進(jìn)行迭代處理，并在每一次迭代后判斷預(yù)設(shè)條 件是否滿足，如是，則停止迭代并轉(zhuǎn)至下一步驟，否則，繼續(xù)下一次迭代:對(duì)于第k次迭代：
[0032] 首先對(duì)矩陣做部分SVD如下：
[0033] l^C = Vi+1S_V_ 廣
[0034] 再對(duì)矩陣CVk做部分SVD如下：
[0035] CVk = Uk+1Stemp2Vtemp2H
[0036] 其中：
[0037]矩陣右上角的數(shù)字標(biāo)號(hào)k表示第k次迭代的初值，每一次的迭代更新該值；
[0038] Uk，Vk表示第k次迭代的初值；
[0039] Uk+1，Vk+1表示第k次迭代的更新之后的結(jié)果，亦作為下一次迭代的初值；
[0040] U和V的規(guī)模均為：（MF)2Xp;
[00411下標(biāo)tempi和temp2表示該矩陣為臨時(shí)變量，不需要保存；
[0042]所述預(yù)設(shè)條件具體如下：
[0045] 范數(shù)2: [(V~0) -I]C#(UiifCV*)s + )~IJC^U"(U^C)； F -F
[0046] 其中，|.|a|表示計(jì)算矩陣A最大的奇異值，| |A| |f表示矩陣A的Frobenius范數(shù)，e為 預(yù)設(shè)閾值；
[0047] 步驟303、截取當(dāng)前矩陣U的前N2列，用Un2表示，即為觀測(cè)信號(hào)的信號(hào)子空間，。 [0048]優(yōu)選地，所述根據(jù)觀測(cè)信號(hào)的信號(hào)子空間和噪聲子空間之間的正交性，估計(jì)出信 號(hào)源的波達(dá)方向及波達(dá)時(shí)間，具體如下：
[0049] 求取使得如下估計(jì)函數(shù)取最大值時(shí)的0、T，即為信號(hào)源的波達(dá)方向及波達(dá)時(shí)間：
[0050] P((l T) = a(6?, TY Ux, Ux；a(^, 7')
[0051] 其中：
[0052] a{0.T) = h{0,T)?h(0J)' ~ eio^e^div^&f
[0053] h(0,T)= : _e(uF)e_27*fI：d;(%,6 ,)
[0054] d(^,^) = [l,e
[0055] 0表示信號(hào)路徑的波達(dá)角度;T表示信號(hào)路徑的波達(dá)時(shí)間；e(Ui)表示信號(hào)在頻率Ui 處的幅值，i = l，2, . . .，F(xiàn);Tu(0)表示信號(hào)路徑到達(dá)第j個(gè)傳感器相對(duì)于到達(dá)參考傳感器的 時(shí)間延遲，j = 2，"、M-l。
[0056] 優(yōu)選地，閾值e的取值范圍為1 e -3~1 e -6。
[0057]優(yōu)選地，所述空域-頻域平滑中的頻域平滑使用頻域子帶平均方法。
[0058] 上述波達(dá)估計(jì)方法可廣泛用于雷達(dá)、聲納、地震等領(lǐng)域，以下為兩個(gè)具體應(yīng)用方 案：
[0059] -種基于聲線傳播時(shí)間層析的海洋聲層析方法，利用聲音在海洋中傳播速度的變 化來反演海洋環(huán)境參數(shù)，首先利用以上任一技術(shù)方案所述基于高階累積量的波達(dá)估計(jì)方法 對(duì)從聲傳感器陣列所接收到的多路徑傳播聲信號(hào)進(jìn)行波達(dá)估計(jì)，從而分離出每一條聲線路 徑;然后根據(jù)聲線路徑的到達(dá)時(shí)間反演出海洋環(huán)境參數(shù)。
[0060] -種定位方法，首先利用以上任一技術(shù)方案所述基于高階累積量的波達(dá)估計(jì)方法 進(jìn)行波達(dá)方向估計(jì)，然后利用估計(jì)出的波達(dá)方向確定信號(hào)源的位置。
[0061] 相比現(xiàn)有技術(shù)，本發(fā)明具有以下有益效果：
[0062] 首先，和基于二階統(tǒng)計(jì)量的技術(shù)相比，本發(fā)明具有更高的信號(hào)路線分離精度，可以 分離出間隔較小的路徑，且本發(fā)明可以在傳感器陣列陣元數(shù)目小于待分離信號(hào)源數(shù)目的情 況下正確分離出信號(hào)的路徑。同時(shí)本發(fā)明中假設(shè)的加性噪聲是彩色噪聲，而傳統(tǒng)的技術(shù)都 是假設(shè)噪聲為高斯噪聲，而彩色噪聲會(huì)更符合實(shí)際情況，本發(fā)明可以更好地抑制各類環(huán)境 噪聲，包括高斯噪聲和非高斯噪聲。
[0063]其次，和基于特征值分解的高階累積量D0A算法相比，本發(fā)明利用迭代地部分SVD 的方法，設(shè)定收斂標(biāo)準(zhǔn)，可以極大地加快算法的計(jì)算速度，減少耗時(shí)。
【附圖說明】
[0064]圖1為空域平滑方法的一個(gè)原理示例；
[0065]圖2為頻域平滑方法的一個(gè)原理示例；
[0066]圖3為本發(fā)明方法的流程示意圖；
[0067]圖4為本發(fā)明方法與Smoothing-MUSICAL算法以及四階累積量基于特征值分解方 法的對(duì)比實(shí)驗(yàn)結(jié)果。
【具體實(shí)施方式】
[0068] 下面結(jié)合附圖對(duì)本發(fā)明的技術(shù)方案進(jìn)行詳細(xì)說明：
[0069] 本發(fā)明的發(fā)明思路是針對(duì)基于高階累積量的陣列信號(hào)處理方法 (CN201610006173.3)所存在的計(jì)算復(fù)雜度過高的問題，在根據(jù)四階累積量矩陣構(gòu)建觀測(cè)信 號(hào)的信號(hào)子空間和噪聲子空間過程中，利用迭代地部分奇異值分解的方法可以大幅減少計(jì) 算復(fù)雜度，從而大幅節(jié)約算法的耗時(shí)。
[0070] 為了便于公眾理解，下面以一個(gè)具體實(shí)施例來對(duì)本發(fā)明技術(shù)方案進(jìn)行詳細(xì)說明。 [0071]首先建立信號(hào)模型：
[0072]各待測(cè)信號(hào)源是相關(guān)的或相干的，設(shè)其數(shù)目為N。傳感器陣列中陣元數(shù)目為M，并不 需要限制1>1傳感器陣列是等間距直線陣，各陣元的特性相同，間隔為d，且陣元間隔不大 于最高頻率信號(hào)的1/2波長(zhǎng)。信號(hào)源的傳播過程中有彩色加性噪聲。
[0073] 第m個(gè)傳感器上接收到的信號(hào)可以表示為： N.
[0074] xm (J) = V anc(t -Tmn) + noise m (〇 (1) n=\
[0075] 其中an為第n個(gè)射線路徑對(duì)第m傳感器的影響，e(t)為信號(hào)源發(fā)射的信號(hào)，^,"為第 n個(gè)射線路徑到達(dá)第m個(gè)傳感器的時(shí)間，noisem(t)表示第m個(gè)傳感器上的加性彩色噪聲。 [0076] 到達(dá)時(shí)間可以表示為：
[0077] Xm,n = Tn+tm(0n) (2)
[0078] 其中Tn表示第n個(gè)射線路徑到達(dá)參考傳感器的時(shí)間，^(0")則表示第n個(gè)射線路徑 到達(dá)第m個(gè)傳感器相比參考傳感器的時(shí)間延遲，角度0"則表示第n個(gè)射線路徑到達(dá)傳感器的 方向。
[0079]式（1)的頻域表示為： M
[0080] xm(u) = ^?"4：u)exp(-/2^：m J + KOM-e,,,!!；) C3) n=l
[0081] 在以上信號(hào)模型的基礎(chǔ)上，本發(fā)明方法具體包括以下步驟：
[0082]步驟1.對(duì)傳感器所接收信號(hào)進(jìn)行空域-頻域的平滑處理，獲得K=(2Ks+l)(2Kf+l) 個(gè)窄帶估計(jì)作為列向量所組成的矩陣X，其中K4PKf分別為空域平滑階次和頻域平滑階次， 且K大于信號(hào)源個(gè)數(shù)N。
[0083]空域-頻域平滑為現(xiàn)有技術(shù)，下面對(duì)其進(jìn)行簡(jiǎn)要說明。
[0084] 空域平滑通過對(duì)空域子陣列的平均來實(shí)現(xiàn)。如圖1所示為空域平滑的原理圖。M個(gè) 傳感器被分成若干尺寸相同，部分重疊的子陣列，假設(shè)子陣列是線性一致的，則子陣列傳感 器上信號(hào)的強(qiáng)度不會(huì)發(fā)生劇烈變化。當(dāng)子陣列數(shù)目大于或等于射線路徑數(shù)時(shí)空間譜矩陣為 非單秩矩陣。假設(shè)空間平滑階次為Ks，則每個(gè)子陣列的尺寸為M-2K S，子陣列的數(shù)目為2KS+1。
[0085] 頻域平滑原理圖如圖2所示，根據(jù)操作對(duì)象是時(shí)域數(shù)據(jù)還是頻域數(shù)據(jù)可分為如下 兩種：（1)加權(quán)相關(guān)矩陣；（2)頻域子帶平均。本發(fā)明優(yōu)選采用頻域子帶平均法。假設(shè)頻域的 平滑階次為K f，則平滑后可以得到2Kf+l個(gè)尺寸為M-2Kf的子帶。
[0086] 對(duì)空域或頻域的過度平滑均會(huì)導(dǎo)致互譜矩陣產(chǎn)生嚴(yán)重的誤差，為了減小估計(jì)誤 差，需要將空域和頻域的平滑相結(jié)合。
[0087] 傳感器陣列的陣元個(gè)數(shù)用M表示，選取的頻率通道數(shù)用F表示，則平滑后以K個(gè)窄帶 估計(jì)作為列向量所組成的矩陣X的規(guī)模為：（MF) XK。
[0088] 步驟2.構(gòu)造平滑后的矩陣X的四階累積量矩陣如下： r n c = ^{(x ? x* )(x ? x* }- e\x ? x* k{(x 0 x* f}
[0089] lv A (( 1 (4) -ir{xxH}祕(mì){xx*}
[0090] 其中：X表示平滑處理之后的觀測(cè)矩陣;C表示四階累積量矩陣;E表示期望；上標(biāo)* 表示共輒;上標(biāo)H表示共輒轉(zhuǎn)置;0表示克羅內(nèi)克積。構(gòu)造出的四階累積量的規(guī)模為：（MF) 2 X(MF)2。
[0091 ]步驟3.該步驟是一個(gè)迭代的過程：
[0092] (1)迭代前先初始化，生成兩個(gè)隨機(jī)的正交矩陣作為U、V的初始值，這意味 著U?HU a = I，VaHVQ = I。這兩個(gè)正交矩陣的規(guī)模均為：（MF)2 Xp。
[0093] 其中：
[0094] p=K2 且 p<<(MF)2。
[0095] I表示單位矩陣。
[0096]矩陣右上角的數(shù)字標(biāo)號(hào)表示迭代的次數(shù)，下面沿用該表示方式。
[0097] (2)迭代的更新（以第k次迭代為例，迭代從第0次開始）：
[0098] 首先對(duì)矩陣做部分SVD如下：
[0099] 細(xì)/ (5、
[0100] 再對(duì)矩陣CVk做部分SVD如下：
[0101] CVk = Uk+1Stemp2Vtemp2H (6)
[0102] 其中：
[0103]矩陣右上角的數(shù)字標(biāo)號(hào)k表示第k次迭代的初值，每一次的迭代更新該值;Uk，Vk表 示第k次迭代的初值;Uk+1，Vk+1表示第k次迭代的更新之后的結(jié)果，亦作為下一次迭代的初 值;U和V的規(guī)模均為：（MF) 2Xp。下標(biāo)tempi和temp2表示該矩陣我們并不關(guān)注。
[0104] (3)終止條件:設(shè)定作為收斂判定標(biāo)準(zhǔn)的閾值e的取值，經(jīng)大量實(shí)驗(yàn)發(fā)現(xiàn)，e的量級(jí) 選取在le-3至le-6左右時(shí)精確度以及計(jì)算復(fù)雜度都會(huì)比較好，由于越接近最優(yōu)值收斂的速 度越慢，如果標(biāo)準(zhǔn)設(shè)定在le-7及更小的數(shù)量級(jí)，迭代的次數(shù)會(huì)大幅增加，相比直接計(jì)算四階 累積量矩陣的特征值其計(jì)算復(fù)雜度的優(yōu)勢(shì)不會(huì)那么明顯，同樣的，如果收斂標(biāo)準(zhǔn)的值過大， 雖然會(huì)極快地地達(dá)到收斂標(biāo)準(zhǔn)停止迭代，但是精確度也會(huì)降低，因此本發(fā)明優(yōu)選在該范圍 內(nèi)取值，例如可取e = 1 e -5。當(dāng)然本領(lǐng)域技術(shù)人員也可根據(jù)需要確定其取值，如需要突出精 確度，則可取較小的值，反之，如需要快速計(jì)算，則可取較大的值。
[0105]計(jì)算兩個(gè)范數(shù)： _6]范數(shù) 2 (7)
[0107] 范數(shù) 2:
[0108] [(VAV^) _ I]CV" (U^CV" )H + [(U*U^) - I]C*U：i (U^^CV4) - S ' F
[0109] 其中；|A.|計(jì)算的是矩陣A最大的奇異值；I |A| If是Frobenius范數(shù)，計(jì)算的是矩陣A 協(xié)方差矩陣的跡的二次方根。
[0110] 每一次迭代我們都作一次判定，若
：則停止迭代。
[0111] 最終得到的(MF)2Xp矩陣U和V則分別是對(duì)應(yīng)于四階累積量矩陣的前p個(gè)最大的奇 異值的左奇異向量子空間和右奇異向量子空間。
[0112] 信號(hào)源個(gè)數(shù)用N來表示，截取矩陣U的前N2列，用UN2表示，即為我們需要的信號(hào)子 空間。
[0113]現(xiàn)有方法是直接進(jìn)行累積量矩陣C的特征值分解計(jì)算，從而來獲得信號(hào)子空間和 噪聲子空間，然而應(yīng)用了高階累積量之后，C的規(guī)模會(huì)變得很大，四階累積量的情況下其規(guī) 模為(MF)2X(MF)2。直接對(duì)其做特征值分解的計(jì)算量耗時(shí)相當(dāng)嚴(yán)重，而在本發(fā)明方法中，并 不直接對(duì)該規(guī)模的矩陣做特征值分解，而是在有限次的迭代過程中做規(guī)模為(MF) 2Xp的矩 陣的部分SVD，由于p<<(MF)2,因此計(jì)算復(fù)雜度可以極大地降低，從而節(jié)約計(jì)算時(shí)的時(shí)間 消耗。
[0114] 步驟4.求取如下估計(jì)函數(shù)的值：
[0115] P(OJ) = mTTVx：Vx：a(OJ) (9)
[0116] 其中：
[0117] a(0S) = h(0,T)?hW,TY 、UJ)
[0118] b(0, T)= : (11) e(v,)e-2j^Td(oF,8)_
[0119] d(u, 0) = [1, ,..., e-2JWir'M-'m ] (12)
[0120] 0表示信號(hào)路徑的波達(dá)角度;T表示信號(hào)路徑的波達(dá)時(shí)間；e(Ui)表示信號(hào)在頻率Ui 處的幅值，i = l，2, . . .，F(xiàn);Tu(0)表示信號(hào)路徑到達(dá)第j個(gè)傳感器相對(duì)于到達(dá)作為參考傳感 器的時(shí)間延遲，j = 2，"_，M-l。
[0121] 估計(jì)函數(shù)的最大值所對(duì)應(yīng)的0和T即信號(hào)源的達(dá)波方向和達(dá)波時(shí)間。整個(gè)處理過程 的流程參見圖3。
[0122]為了驗(yàn)證本發(fā)明的效果，將其與現(xiàn)有的smoothing-MUSICAL算法進(jìn)行仿真對(duì)比實(shí) 驗(yàn)，同時(shí)與基于特征值分解的四階累積量算法進(jìn)行對(duì)比。實(shí)驗(yàn)使用的射線路徑的相關(guān)信息 為:五條射線路徑在各個(gè)傳感器間的延遲時(shí)間分別為:〇、1、-1、1.6、-1.6;五條射線路徑到 達(dá)傳感器的時(shí)間分別為10、12、12.5、14、14;選取第一個(gè)傳感器為參考傳感器;傳感器米樣 數(shù)據(jù)長(zhǎng)度為129;兩種方法所采用的空域和頻域的平滑階次均為1。且該實(shí)驗(yàn)中收斂標(biāo)準(zhǔn)設(shè) 定為le-3,可以更突出其在計(jì)算復(fù)雜度方面的提升。算法運(yùn)行的操作系統(tǒng)是Red Hat Enterprise Linux Server release 5.6(Tikanga)，內(nèi)存24GB，處理器是Intel X56500 2.67GHz。
[0123] 對(duì)比實(shí)驗(yàn)的結(jié)果如圖4所示，圖中每行的3個(gè)結(jié)果分別是smoo th i ng-MUS I CAL (左側(cè) 的結(jié)果)和四階累積量基于特征值分解（中間的結(jié)果)和本發(fā)明方法(右側(cè)的結(jié)果）的方法使 用同樣的觀測(cè)數(shù)據(jù)得出的結(jié)果。同樣的實(shí)驗(yàn)對(duì)象是5個(gè)信號(hào)源：傳感器間的延遲時(shí)間分別 為:0、1、-1、1 ? 6、-1 ? 6，信號(hào)到達(dá)傳感器的時(shí)間分別為10、12、12 ? 5、14、14的五條射線路徑， 圖中從最上面一行至最后一行的實(shí)驗(yàn)條件依次如下：
[0124] (1)不加噪聲，傳感器數(shù)目為4，選取第一個(gè)傳感器為參考傳感器；
[0125] (2)加彩色噪聲，信噪比為20dB，傳感器數(shù)目為4,選取第一個(gè)傳感器為參考傳感 器；
[0126] (3)加彩色噪聲，信噪比為_5dB，傳感器數(shù)目為4，選取第一個(gè)傳感器為參考傳感 器；
[0127] (4)加彩色噪聲，信噪比為-10dB，傳感器數(shù)目為4,選取第一個(gè)傳感器為參考傳感 器；
[0128] (5)不加噪聲，傳感器數(shù)目為6，選取第一個(gè)傳感器為參考傳感器；
[0129] (6)加彩色噪聲，信噪比為20dB，傳感器數(shù)目為6,選取第一個(gè)傳感器為參考傳感 器；
[0130] (7)加彩色噪聲，信噪比為_5dB，傳感器數(shù)目為6,選取第一個(gè)傳感器為參考傳感 器；
[0131] (8)加彩色噪聲，信噪比為-10dB，傳感器數(shù)目為6,選取第一個(gè)傳感器為參考傳感 器；
[0132] 表1為該實(shí)驗(yàn)中四階累積量D0A方法分別基于特征值分解（即CN201610006173.3中 的方法)和基于部分SVD迭代算法(即本發(fā)明方法)獲取信號(hào)子空間過程的時(shí)間消耗。
[0133] 表1.四階累積量D0A算法分別基于特征值分解和部分SVD迭代的時(shí)間消耗對(duì)比
[0135] 根據(jù)以上對(duì)比實(shí)驗(yàn)可知，與基于二階統(tǒng)計(jì)量的技術(shù)相比，本發(fā)明具有更高的信號(hào) 路線分離精度，可以分離出間隔較小的路徑，且本發(fā)明可以在傳感器陣列陣元數(shù)目小于待 分離信號(hào)源數(shù)目的情況下正確分離出信號(hào)的路徑。同時(shí)本發(fā)明可以抑制非高斯噪聲的影 響。
[0136] 和基于特征值分解的高階累積量D0A算法相比，本發(fā)明基于部分SVD迭代的方法， 可以極大地加快算法的計(jì)算速度，減少耗時(shí)。
[0137] 本發(fā)明波達(dá)估計(jì)方法可可廣泛用于雷達(dá)、聲納、地震等領(lǐng)域，以下為兩個(gè)具體應(yīng)用 方案：
[0138] -種基于聲線傳播時(shí)間層析的海洋聲層析方法，利用聲音在海洋中傳播速度的變 化來反演海洋環(huán)境參數(shù)，首先利用以上任一技術(shù)方案所述基于高階累積量的波達(dá)估計(jì)方法 對(duì)從聲傳感器陣列所接收到的多路徑傳播聲信號(hào)進(jìn)行波達(dá)估計(jì)，從而分離出每一條聲線路 徑;然后根據(jù)聲線路徑的到達(dá)時(shí)間反演出海洋環(huán)境參數(shù)。
[0139] -種定位方法，首先利用以上任一技術(shù)方案所述基于高階累積量的波達(dá)估計(jì)方法 進(jìn)行波達(dá)方向估計(jì)，然后利用估計(jì)出的波達(dá)方向確定信號(hào)源的位置。
【主權(quán)項(xiàng)】
1. 一種基于高階累積量的波達(dá)估計(jì)方法，利用等間距直線傳感器陣列所接收到的觀測(cè) 信號(hào)，估計(jì)出信號(hào)源的波達(dá)方向及波達(dá)時(shí)間;包括以下步驟： 步驟1、對(duì)所述觀測(cè)信號(hào)做傅里葉變換后進(jìn)行空域-頻域平滑處理，得到K=(2Ks+l)(2Kf +1)個(gè)窄帶估計(jì)作為列向量所組成的矩陣X，Ks、Kf分別為空域平滑階次、頻域平滑階次;且K 大于信號(hào)源個(gè)數(shù)N; 步驟2、根據(jù)下式構(gòu)造出空域-頻域平滑處理后信號(hào)的四階累積量矩陣C:X表示空域-頻域平滑處理之后得到的矩陣，E表示期望，上標(biāo)*表示共輒，上標(biāo)Η表示共 輒轉(zhuǎn)置，發(fā)表示克羅內(nèi)克積； 步驟3、根據(jù)所述四階累積量矩陣構(gòu)建觀測(cè)信號(hào)的信號(hào)子空間和噪聲子空間； 步驟4、根據(jù)觀測(cè)信號(hào)的信號(hào)子空間和噪聲子空間之間的正交性，估計(jì)出信號(hào)源的波達(dá) 方向及波達(dá)時(shí)間； 其特征在于，所述步驟3具體包括以下子步驟： 步驟301、對(duì)兩個(gè)規(guī)模均為(MF)2 X ρ的正交矩陣U、V進(jìn)行隨機(jī)初始化，記為U*3、V*3，矩陣右 上角的數(shù)字標(biāo)號(hào)表示迭代的次數(shù)，Μ為所述等間距直線傳感器陣列的陣元個(gè)數(shù)，F(xiàn)為空域-頻 域平滑處理過程中選取的頻率通道數(shù)，Ρ=Κ 2且p<<(MF)2; 步驟302、對(duì)所述兩個(gè)正交矩陣U、V進(jìn)行迭代處理，并在每一次迭代后判斷預(yù)設(shè)條件是 否滿足，如是，則停止迭代并轉(zhuǎn)至下一步驟，否則，繼續(xù)下一次迭代:對(duì)于第k次迭代： 首先對(duì)矩陣If 做部分SVD如下：再對(duì)矩陣CVk做部分SVD如下： 其中：矩陣右上角的數(shù)字標(biāo)號(hào)k表示第k次迭代的初值，每一次的迭代更新該值； Uk，Vk表示第k次迭代的初值； Uk+1，Vk+1表示第k次迭代的更新之后的結(jié)果，亦作為下一次迭代的初值； U和V的規(guī)模均為：（MF)2Xp; 下標(biāo)tempi和temp2表示該矩陣為臨時(shí)變量，不需要保存； 所述預(yù)設(shè)條件具體如下：其中，|Ag表示計(jì)算矩陣A最大的奇異值，|^47表示矩陣A的Frobenius范數(shù)，ε為預(yù)設(shè)閾 值； 步驟303、截取當(dāng)前矩陣U的前N2列，用Un2表示，即為觀測(cè)信號(hào)的信號(hào)子空間。2. 如權(quán)利要求1所述波達(dá)估計(jì)方法，其特征在于，所述根據(jù)觀測(cè)信號(hào)的信號(hào)子空間和噪 聲子空間之間的正交性，估計(jì)出信號(hào)源的波達(dá)方向及波達(dá)時(shí)間，具體如下： 求取使得如下估計(jì)函數(shù)取最大值時(shí)的Θ、Τ，即為信號(hào)源的波達(dá)方向及波達(dá)時(shí)間：Θ表示信號(hào)路徑的波達(dá)角度;Τ表示信號(hào)路徑的波達(dá)時(shí)間；e(Ui)表示信號(hào)在頻率Ui處的 幅值，i = l，2, . . .，F(xiàn);t^(0)表示信號(hào)路徑到達(dá)第j個(gè)傳感器相對(duì)于到達(dá)參考傳感器的時(shí)間 延遲，j = 2，"_，M-l。3. 如權(quán)利要求1所述波達(dá)估計(jì)方法，其特征在于，閾值ε的取值范圍為1 e-3~1 e-6。4. 如權(quán)利要求1所述波達(dá)估計(jì)方法，其特征在于，所述空域-頻域平滑中的頻域平滑使 用頻域子帶平均方法。5. -種基于聲線傳播時(shí)間層析的海洋聲層析方法，利用聲音在海洋中傳播速度的變化 來反演海洋環(huán)境參數(shù)，其特征在于，首先利用權(quán)利要求1~4任一項(xiàng)所述基于高階累積量的 波達(dá)估計(jì)方法對(duì)從聲傳感器陣列所接收到的多路徑傳播聲信號(hào)進(jìn)行波達(dá)估計(jì)，從而分離出 每一條聲線路徑;然后根據(jù)聲線路徑的到達(dá)時(shí)間反演出海洋環(huán)境參數(shù)。6. -種定位方法，其特征在于，首先利用權(quán)利要求1~4任一項(xiàng)所述基于高階累積量的 波達(dá)估計(jì)方法進(jìn)行波達(dá)方向估計(jì)，然后利用估計(jì)出的波達(dá)方向確定信號(hào)源的位置。
【文檔編號(hào)】G01S7/539GK105929386SQ201610230971
【公開日】2016年9月7日
【申請(qǐng)日】2016年4月14日
【發(fā)明人】劉杰, 姜龍玉, 洪亞萍, 伍家松, 舒華忠
【申請(qǐng)人】東南大學(xué)