稀疏貝葉斯學習用于噪聲互相關提取格林函數

郭啟超，李風華，楊習山

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

0 引言

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

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

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

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

1 理論

1.1 被動TDGF提取

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

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

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

1.2 TDGF稀疏模型

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

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

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

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

1.3 SBL求解模型

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

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

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

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

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

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

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

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

2 仿真實驗

2.1 仿真環境介紹

圖1 實驗布放圖及聲速剖面Fig.1 Layout of the experiment and sound speed profile

2.2 仿真結果

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

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

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

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

3 海試實驗

3.1 海試實驗數據介紹

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

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

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

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

3.3 預累積處理對SBL提取TDGF的影響

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

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

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

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

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

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

3.4 結果

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

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

圖8為利用SBL估計TDGF的第二條聲線到達時間反演得到的海水平均聲速與TD實測的海水平均聲速的比較，二者之間的均方根誤差為0.6 m·s-1。實驗結果證實了相比于傳統方法，預累積SBL提取的TDGF可以更準確地獲取較窄頻帶下的聲線到達結構。

圖8 利用SBL提取的TDGF第二號到達時間反演海水平均聲速與溫深儀實測聲速對比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 結論

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