稀疏貝葉斯學(xué)習(xí)用于噪聲互相關(guān)提取格林函數(shù)

郭啟超，李風(fēng)華，楊習(xí)山

(1.中國科學(xué)院聲學(xué)研究所，北京 100190；2.中國科學(xué)院大學(xué)，北京 100190)

0 引言

兩點(diǎn)之間的時(shí)域格林函數(shù)(Time Domain Green’s Function,TDGF)可以通過此兩點(diǎn)之間環(huán)境噪聲的互相關(guān)函數(shù)(Noise Cross-correlation Function,NCF)累積得到[1-3]，這為海洋被動聲層析[4-6]提供了一種實(shí)用方法。該方法的主要問題是NCF累積時(shí)間較長。目前，頻域白化等預(yù)處理方法可以一定程度上改善TDGF的提取質(zhì)量[7]。空域?yàn)V波方法[8-10]是一種可以顯著加快TDGF提取的有效技術(shù)。此外實(shí)驗(yàn)數(shù)據(jù)表明縮小帶寬有助于濾除非擴(kuò)散噪聲干擾，從而提高NCF累積效果。但是，當(dāng)頻率帶寬太窄時(shí)，傳統(tǒng)的反傅里葉變換(Inverse Fast Fourier Transform,IFFT)無法區(qū)分TDGF中相鄰的到達(dá)時(shí)間，影響了NCF方法在海洋聲被動層析中的應(yīng)用。

由于NCF方法從環(huán)境噪聲中估計(jì)出少量的TDGF到達(dá)時(shí)間，該過程可被視為壓縮感知(Compressed Sensing,CS)理論中的稀疏重構(gòu)問題[11-13]。并且由于短時(shí)間內(nèi)TDGF是近似穩(wěn)態(tài)的，因此NCF累積過程可以看做CS中的多觀測向量模型(Multiple Measurement Vector,MMV)問題。稀疏貝葉斯學(xué)習(xí)(Sparse Bayesian Learning,SBL)[14-15]是一種 CS 算法，其通過貝葉斯框架尋找欠定線性系統(tǒng)的稀疏解。SBL與人的思考模式類似，首先假設(shè)未知變量服從某種先驗(yàn)分布，然后通過貝葉斯公式對樣本數(shù)據(jù)進(jìn)行“學(xué)習(xí)”，不斷更新已有認(rèn)知，最終作出對未知參數(shù)的推斷。與基于L1懲罰項(xiàng)的CS算法(比如Lasso等)相比，SBL在MMV問題中能自動尋找稀疏度(如TDGF到達(dá)時(shí)間個(gè)數(shù))[16]，保持較高的計(jì)算效率[17]。在水聲學(xué)領(lǐng)域，SBL算法已成功應(yīng)用于波達(dá)方向估計(jì)[18]、聲源定位[19]、水平波數(shù)[20]分離方面。

本文利用SBL從較窄頻帶(帶寬為70 Hz)的海洋環(huán)境噪聲中提取兩接收器之間的TDGF到達(dá)結(jié)構(gòu)，首先構(gòu)造了NCF方法的多觀測向量模型，其中字典矩陣由傅里葉變換算子構(gòu)成，觀測矩陣由預(yù)累積后的頻域NCF組成。預(yù)累積是指對于給定總快拍數(shù)的頻域NCF，首先以某個(gè)固定的快拍數(shù)間隔累積，依次生成每一個(gè)觀測向量。預(yù)累積目的是增強(qiáng)每一個(gè)觀測向量的信噪比，從而提高SBL的分辨率。為了保證SBL的穩(wěn)健性，同時(shí)預(yù)累積快拍數(shù)也不能過大。實(shí)際使用中應(yīng)該選取滿足分辨率門限的最小預(yù)累積快拍數(shù)。

聯(lián)合預(yù)累積處理的SBL方法兼顧了噪聲提取TDGF的分辨率與穩(wěn)健性，能夠有效獲取較窄頻帶下傳統(tǒng)方法無法分辨的TDGF到達(dá)時(shí)間。仿真與海試實(shí)驗(yàn)數(shù)據(jù)驗(yàn)證了該方法的性能，并且分離得到的到達(dá)時(shí)間可以用于海水平均聲速反演，突破了傳統(tǒng)方法提取TDGF到達(dá)時(shí)間精度的限制。該方法拓展了SBL的應(yīng)用場景，并且可以與空域?yàn)V波方法結(jié)合，為自適應(yīng)選取噪聲頻段，實(shí)現(xiàn)快速海洋被動聲層析提供了一種可行方法。

1 理論

1.1 被動TDGF提取

本小節(jié)介紹時(shí)域格林函數(shù)(TDGF)被動提取的基本方法。兩點(diǎn)之間的歸一化噪聲互相關(guān)函數(shù)(NCF)在頻域可以表示為[8]

式中：yl(ω)表示一個(gè)快拍的頻域NCF；L為快拍數(shù)；p1(ω)與p2(ω)分別表示兩個(gè)接收器記錄的噪聲信號的頻域形式；(·)*表示求共軛復(fù)數(shù)，歸一化的目的是通過頻域白化的預(yù)處理方法消除強(qiáng)頻率能量的干擾。利用兩條陣列指向端射方向的波束輸出互相關(guān)代替原本的兩個(gè)單接收器信號互相關(guān)，可以顯著地提高每一快拍NCF的質(zhì)量[18]。對多個(gè)快拍的頻域NCF進(jìn)行累積，然后求反傅里葉變換(IFFT)，即可獲得兩接收器位置之間的TDGF[21]：

式中：x(t)與x(-t)分別表示因果TDGF與非因果TDGF，它們關(guān)于原點(diǎn)對稱；jω為相移因子。IFFT方法存在的問題是當(dāng)頻帶較窄時(shí)，TDGF分辨率低而無法直接應(yīng)用于聲速反演。

1.2 TDGF稀疏模型

本小節(jié)將建立頻域NCF提取TDGF的稀疏表示模型。給定NCF的總快拍數(shù)L，首先將NCF以每ΔL(1≤ΔL≤L)快拍數(shù)的間隔累積，重新得到 L′=L/ΔL快拍的NCF。

式(3)的作用是調(diào)節(jié)稀疏貝葉斯學(xué)習(xí)(SBL)獲取的TDGF的分辨率與穩(wěn)健性(在3.3節(jié)中證明)。只考慮因果TDGF，假設(shè)每一快拍新的NCF可以看作TDGF與加性高斯噪聲的形式，將式(3)代入式(2)，并且對式(2)兩邊進(jìn)行傅里葉變換，可得到：

式中：nl′(ω)表示加性高斯白噪聲。為了方便表述，用 yl′，xl′與 nl′分別表示 -jωyl′(ω)，xl′(t)與 nl′(ω)。令Y=[y1… yL′]∈ CN×L′為觀測矩陣，其每一列的觀測向量等于預(yù)累積后的頻域NCF乘上因子-jω，N為NCF包含的頻率點(diǎn)數(shù)。令 X=[x1… xL′]∈CM×L′，其每一項(xiàng)xm，l′為第l′快拍的TDGF在第m個(gè)時(shí)間點(diǎn)的復(fù)數(shù)幅度，并且m=1，…，M,M為TDGF包含的時(shí)間點(diǎn)數(shù)。此時(shí)，將式(4)中的傅里葉變換離散化后，可以得到：

式中：A∈CN×M，為字典矩陣，其每一項(xiàng) an，m=e-j2π[m(n-1)]/N聯(lián)系著觀測向量中第n個(gè)頻點(diǎn)與TDGF中第m個(gè)時(shí)間點(diǎn)之間的傅里葉變換。Z=[n1… nL′]∈ CN×L′， 為 加 性 復(fù) 高 斯 白 噪 聲 。 假 設(shè)N?M，則式(5)是欠定的。由于穩(wěn)態(tài)的TDGF中只有少量的K個(gè)時(shí)間點(diǎn)的能量不為0，因此每一快拍TDGF是K-稀疏的，即K?M。模型(5)為壓縮感知框架中經(jīng)典的多觀測向量模型。

1.3 SBL求解模型

本小節(jié)簡要概括利用SBL求解TDGF模型的方法。假設(shè)每一快拍噪聲nl′符合均值為0、方差為σ2的 復(fù) 高 斯 分 布 ， 即 nl′～ CN(0，σ2I)， 則 似 然 分布為[17]

選擇我院接受的進(jìn)行超聲檢查的300例甲狀腺結(jié)節(jié)患者作為案例，對患者采用本院現(xiàn)有的超聲量化分級系統(tǒng)實(shí)施診斷和評估。本次研究的研究案例均符合要求，經(jīng)過病理診斷后具備知情權(quán)，簽署知情同意書。

式中：I表示單位向量。假設(shè)每一快拍TDGF的幅度服從復(fù)高斯先驗(yàn)分布，且不同快拍之間相互獨(dú)立，則多快拍信號X的先驗(yàn)分布為[17]

式中：Γ=diag(γ)是一個(gè)對角協(xié)方差矩陣，其對角元素γ=[γ1…γM]T表示TDGF在每個(gè)時(shí)間點(diǎn)處的能量。由于TDGF的稀疏性，在算法運(yùn)行中，大部分時(shí)間點(diǎn)的能量γ會趨于0。根據(jù)高斯先驗(yàn)與高斯似然分布，得到數(shù)據(jù)的邊緣分布P(Y)也是高斯分布[17]：

式中：Σy=AΓAH+σ2I表示模型方差。SBL通過最大化數(shù)據(jù)邊緣分布得到TDGF的能量γ[17]：

令式(9)中目標(biāo)函數(shù)的導(dǎo)數(shù)為0，可得到能量γ的迭代更新解[17]：

式(10)中的每次迭代都需要知道噪聲方差的估計(jì)值σ2，該變量可以通過每次更新的信號能量計(jì)算[17]：

式中：AM為非0位置對應(yīng)的字典向量的集合，Ky=YYH/L′表示樣本協(xié)方差矩陣，(·)+表示矩陣的偽逆操作。實(shí)際使用中，給定信號能量γ與噪聲方差σ2的初始值(例如設(shè)置γ=1，σ2=0.1)，直到前后兩次迭代γ的改善值收斂到一個(gè)預(yù)設(shè)的閾值。此時(shí)，信號X的最大后驗(yàn)估計(jì)為[17]

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

2.1 仿真環(huán)境介紹

圖1 實(shí)驗(yàn)布放圖及聲速剖面Fig.1 Layout of the experiment and sound speed profile

2.2 仿真結(jié)果

圖2給出了不同頻帶下基于IFFT方法提取TDGF結(jié)果，圖中粗實(shí)線是基于噪聲互相關(guān)提取的兩個(gè)接收器之間的TDGF，細(xì)虛線是根據(jù)將一個(gè)接收器作為聲源，計(jì)算得到的相應(yīng)接收器之間的信道響應(yīng)，黑色圓點(diǎn)表示BELLHOP理論計(jì)算得到的各條聲線到達(dá)時(shí)間。從圖2可以看出，被動TDGF波形與相應(yīng)的信道響應(yīng)基本一致，均存在四個(gè)明顯的到達(dá)波包(在圖中用1、2、3、4標(biāo)注)，分別對應(yīng)4條主要聲波路徑的到達(dá)時(shí)間，但兩者對應(yīng)波包(尤其是第1、第2個(gè)波包)的幅度有區(qū)別，這與海面噪聲源提取TDGF的理論相符[22]。仿真結(jié)果同時(shí)也表明，在較寬的頻帶情況下(20～400 Hz)，基于IFFT方法可以獲得較為準(zhǔn)確的TDGF峰值到達(dá)時(shí)間，而在較窄頻帶條件下(110～180 Hz)，基于IFFT方法無法分辨TDGF的時(shí)間結(jié)構(gòu)。

圖2 IFFT法提取TDGF與通過將一個(gè)接收器作為聲源計(jì)算的信道響應(yīng)對比Fig.2 Comparison between the simulated TDGF extracted by IFFT and the channel impulse response obtained by setting one receiver as sound source.

圖3給出了基于本文提出的預(yù)累積SBL方法得到的TDGF時(shí)間結(jié)構(gòu)。在110～180 Hz頻段下，給定NCF的總快拍數(shù)(1 440)，圖3中顯示了經(jīng)過不同快拍數(shù)ΔL預(yù)累積后，基于IFFT(細(xì)虛線)與SBL(粗實(shí)線)的提取結(jié)果，黑色圓點(diǎn)表示計(jì)算得到的理論值。從圖3可以看出，基于IFFT方法的計(jì)算結(jié)果不受ΔL影響。對于SBL，隨著ΔL增加，TDGF的包絡(luò)逐漸變窄，當(dāng)ΔL等于60、90或120時(shí)，其第二個(gè)包絡(luò)成功分辨TDGF的第二個(gè)到達(dá)時(shí)間。但是，也不是ΔL越大越好，當(dāng)ΔL過大時(shí)(比如本文中ΔL大于720)，基于SBL的估計(jì)結(jié)果可能與真值產(chǎn)生偏差。

圖3 不同預(yù)累積量下利用IFFT與SBL從NCFs中提取的TDGF仿真Fig.3 Simulated TDGF extracted from NCFs under different prestacking snapshots by IFFT and SBL

3 海試實(shí)驗(yàn)

3.1 海試實(shí)驗(yàn)數(shù)據(jù)介紹

本節(jié)利用海試實(shí)驗(yàn)數(shù)據(jù)驗(yàn)證本文提出的預(yù)累積SBL算法的有效性。實(shí)驗(yàn)于2018年5月在中國南海某海域開展。有兩套接收系統(tǒng)同步記錄海洋環(huán)境噪聲，在實(shí)驗(yàn)過程中利用溫深儀(Temperature Depth sensor,TD)每30 s同步測量一次海水聲速剖面。數(shù)據(jù)處理中，首先將噪聲信號分割成若干快拍，每一快拍長度為10 s。利用1.1節(jié)中的方法產(chǎn)生每一快拍頻域NCF。

3.2 不同頻帶下提取的TDGF對比

對比不同頻帶下IFFT與SBL的TDGF提取結(jié)果。圖4(a)、4(b)分別為20～400 Hz與20～90 Hz頻段時(shí)、使用IFFT與SBL，從720快拍(2 h)NCF中提取的TDGF，散點(diǎn)為使用BELLHOP根據(jù)實(shí)測水文仿真的TDGF到達(dá)時(shí)間理論值。當(dāng)頻段為20～400 Hz時(shí)，IFFT與SBL中均存在4個(gè)明顯的獨(dú)立波包(在圖中分別用1、2、3、4標(biāo)注)，分別與Bellhop計(jì)算得到的4個(gè)到達(dá)時(shí)間相吻合。當(dāng)頻段為20～90 Hz時(shí)，IFFT無法分辨各到達(dá)時(shí)間，而SBL仍然能夠較好分辨各到達(dá)時(shí)間。

圖4 不同頻段的海洋環(huán)境噪聲中提取的TDGFFig.4 TDGF extracted from ocean ambient noise in different frequency bands

3.3 預(yù)累積處理對SBL提取TDGF的影響

本節(jié)討論預(yù)累積量ΔL與基于SBL方法提取TDGF的分辨率與穩(wěn)健性之間的關(guān)系。圖5分別給出了經(jīng)過不同ΔL預(yù)累積之后，基于IFFT(細(xì)虛線)與SBL(粗實(shí)線)方法得到TDGF，信號頻段為110～180 Hz，散點(diǎn)表示BELLHOP理論值。數(shù)據(jù)處理結(jié)果與仿真結(jié)果一致，IFFT均無法分辨TDGF的前兩個(gè)到達(dá)時(shí)間，SBL在ΔL為30或120時(shí)，可以較好地估計(jì)第二個(gè)到達(dá)時(shí)間(包絡(luò)的波谷位置)。而當(dāng)ΔL=1(不進(jìn)行預(yù)累積)或ΔL=720(全部預(yù)累積)時(shí)，估計(jì)結(jié)果分別出現(xiàn)了分辨率低與偏差大的問題。

圖5 不同預(yù)累積量時(shí)，利用IFFT與SBL從海洋環(huán)境噪聲中提取TDGF的結(jié)果Fig.5 The results of TDGF extracted from ocean ambient noise with different ΔL by IFFT and SBL

定義TDGF的包絡(luò)寬度(Band Width,BW)為其Hilbert包絡(luò)絕對值下降3 dB時(shí)兩個(gè)到達(dá)時(shí)延的差值：

圖6為經(jīng)過70 h平均后，SBL估計(jì)的第二號到達(dá)時(shí)間在包絡(luò)寬度內(nèi)的最大誤差隨ΔL的變化情況。隨著ΔL的增大，最大誤差先減小后增大。原因是SBL估計(jì)的TDGF的分辨率與穩(wěn)健性是一對相反關(guān)系，即隨著ΔL的增大，TDGF的分辨率提高但穩(wěn)健性降低。當(dāng)ΔL較小時(shí)，分辨率的提高起主要作用，也就是說，包絡(luò)變窄使得TDGF的“隨機(jī)誤差”范圍變??；而當(dāng)ΔL較大時(shí)，穩(wěn)健性的損失起主要影響，即包絡(luò)中心位置與真值之間偏差使得TDGF的“系統(tǒng)誤差”增大。圖6同時(shí)也說明，通過選取合適的ΔL，SBL可以在分辨率與穩(wěn)健性間進(jìn)行折中。

圖6 70 h平均的SBL估計(jì)的TDGF第二號到達(dá)時(shí)間的最大誤差隨預(yù)累積量的變化Fig.6 Variation of the 70 h averaged maximum error of the second arrival time of TDGF estimated by SBL with differentΔL

進(jìn)一步討論預(yù)累積量ΔL對基于SBL方法提取TDGF的分辨率與穩(wěn)健性的影響。SBL假設(shè)待估計(jì)的稀疏向量中的每個(gè)元素都服從一個(gè)均值為0、方差為γ的高斯分布。在算法運(yùn)行中，γ中的絕大部分值會變成0(無噪情況下)或者趨于0(有噪情況下)。在給定總快拍數(shù)條件下，一方面，預(yù)累積快拍數(shù)越多，則每一個(gè)觀測向量的信噪比越高，則非真值位置的能量越容易趨于0，也可以理解為數(shù)據(jù)信噪比越高，則SBL的“學(xué)習(xí)”效果越好，越能確切地將靠近真值附近的非真值點(diǎn)的能量置0，表現(xiàn)在波形上就是包絡(luò)窄，分辨率高；另一方面，若預(yù)累積后的各觀測向量已經(jīng)具有較好的信噪比，那么將多個(gè)觀測向量組成矩陣輸入SBL，比簡單地累積為一維的向量效果好。原因是：(1)在穩(wěn)態(tài)情況下，不同觀測向量之間有共同的稀疏性[23]，簡單累積為一維向量會“浪費(fèi)”共同稀疏性的約束，導(dǎo)致估計(jì)的非0位置不穩(wěn)定。(2)預(yù)累積會導(dǎo)致樣本協(xié)方差矩陣維度減小，估計(jì)的噪聲方差不穩(wěn)定(式(11))，而噪聲方差又會影響SBL對非0位置的估計(jì)(式(10))。因此在實(shí)際使用中，應(yīng)該首先根據(jù)頻帶范圍確定分辨率的最低門限，然后選取達(dá)到分辨率門限的最小ΔL。

3.4 結(jié)果

圖7給出了使用IFFT與SBL從每2 h NCFs中提取一次TDGF，共計(jì)70 h的提取結(jié)果，其中黑色點(diǎn)線為利用實(shí)測聲速剖面通過BELLHOP軟件計(jì)算的TDGF到達(dá)時(shí)間理論值，噪聲數(shù)據(jù)的頻段為110～180 Hz，SBL的預(yù)先累積量ΔL設(shè)置為12，為了盡可能提高SBL的穩(wěn)定性，減小聲速反演誤差，ΔL的選取小于圖6中的最優(yōu)解60。從圖7可以看出，基于SBL方法可以清晰看到4條聲線的到達(dá)時(shí)間結(jié)構(gòu)，而IFFT只能看到較寬的三個(gè)區(qū)域，無法準(zhǔn)確分辨各條聲線的到達(dá)時(shí)間，其原因是基于IFFT的提取方法對應(yīng)的不同聲線信號之間發(fā)生“混疊”。

圖7 使用IFFT與SBL從海洋環(huán)境噪聲中提取的共計(jì)70 h的TDGF結(jié)果Fig.7 The results of TDGF extracted from 70 h ocean ambient noise by IFFT and SBL

圖8為利用SBL估計(jì)TDGF的第二條聲線到達(dá)時(shí)間反演得到的海水平均聲速與TD實(shí)測的海水平均聲速的比較，二者之間的均方根誤差為0.6 m·s-1。實(shí)驗(yàn)結(jié)果證實(shí)了相比于傳統(tǒng)方法，預(yù)累積SBL提取的TDGF可以更準(zhǔn)確地獲取較窄頻帶下的聲線到達(dá)結(jié)構(gòu)。

圖8 利用SBL提取的TDGF第二號到達(dá)時(shí)間反演海水平均聲速與溫深儀實(shí)測聲速對比Fig.8 Comparison between the averaged sound speed inverted by the second TDGF arrival time estimated by SBL and the measured sound speed by bathy thermograph

4 結(jié)論

為解決較窄頻帶下傳統(tǒng)的噪聲互相關(guān)提取TDGF方法時(shí)域分辨率低的問題，提出了一種基于SBL的噪聲互相關(guān)TDGF提取方法。首先建立了TDGF的稀疏模型，其中字典矩陣由傅里葉變換算子組成，觀測矩陣由頻域噪聲互相關(guān)函數(shù)組成。然后提出利用預(yù)累積處理來折中SBL估計(jì)TDGF的分辨率與穩(wěn)定性。仿真與實(shí)驗(yàn)數(shù)據(jù)處理結(jié)果表明：通過選取合適的預(yù)累積量，基于SBL方法可以準(zhǔn)確、穩(wěn)健地從較窄頻帶的海洋環(huán)境噪聲中獲取TDGF的時(shí)間到達(dá)結(jié)構(gòu)，有效克服傳統(tǒng)方法時(shí)間精度低的問題，為之后開展自適應(yīng)挑選噪聲頻段，實(shí)現(xiàn)快速被動聲層析提供了一種可行方法。