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

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

0 引言

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

由于NCF方法從環(huán)境噪聲中估計出少量的TDGF到達時間，該過程可被視為壓縮感知(Compressed Sensing,CS)理論中的稀疏重構(gòu)問題[11-13]。并且由于短時間內(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è)未知變量服從某種先驗分布，然后通過貝葉斯公式對樣本數(shù)據(jù)進行“學(xué)習(xí)”，不斷更新已有認知，最終作出對未知參數(shù)的推斷。與基于L1懲罰項的CS算法(比如Lasso等)相比，SBL在MMV問題中能自動尋找稀疏度(如TDGF到達時間個數(shù))[16]，保持較高的計算效率[17]。在水聲學(xué)領(lǐng)域，SBL算法已成功應(yīng)用于波達方向估計[18]、聲源定位[19]、水平波數(shù)[20]分離方面。

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

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

1 理論

1.1 被動TDGF提取

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

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

式中：x(t)與x(-t)分別表示因果TDGF與非因果TDGF，它們關(guān)于原點對稱；jω為相移因子。IFFT方法存在的問題是當(dāng)頻帶較窄時，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)兩邊進行傅里葉變換，可得到：

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

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

1.3 SBL求解模型

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

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

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

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

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

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

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

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

2 仿真實驗

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

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

2.2 仿真結(jié)果

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

圖2 IFFT法提取TDGF與通過將一個接收器作為聲源計算的信道響應(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時間結(jié)構(gòu)。在110～180 Hz頻段下，給定NCF的總快拍數(shù)(1 440)，圖3中顯示了經(jīng)過不同快拍數(shù)ΔL預(yù)累積后，基于IFFT(細虛線)與SBL(粗實線)的提取結(jié)果，黑色圓點表示計算得到的理論值。從圖3可以看出，基于IFFT方法的計算結(jié)果不受ΔL影響。對于SBL，隨著ΔL增加，TDGF的包絡(luò)逐漸變窄，當(dāng)ΔL等于60、90或120時，其第二個包絡(luò)成功分辨TDGF的第二個到達時間。但是，也不是ΔL越大越好，當(dāng)ΔL過大時(比如本文中ΔL大于720)，基于SBL的估計結(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 海試實驗

3.1 海試實驗數(shù)據(jù)介紹

本節(jié)利用海試實驗數(shù)據(jù)驗證本文提出的預(yù)累積SBL算法的有效性。實驗于2018年5月在中國南海某海域開展。有兩套接收系統(tǒng)同步記錄海洋環(huá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頻段時、使用IFFT與SBL，從720快拍(2 h)NCF中提取的TDGF，散點為使用BELLHOP根據(jù)實測水文仿真的TDGF到達時間理論值。當(dāng)頻段為20～400 Hz時，IFFT與SBL中均存在4個明顯的獨立波包(在圖中分別用1、2、3、4標注)，分別與Bellhop計算得到的4個到達時間相吻合。當(dāng)頻段為20～90 Hz時，IFFT無法分辨各到達時間，而SBL仍然能夠較好分辨各到達時間。

圖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(細虛線)與SBL(粗實線)方法得到TDGF，信號頻段為110～180 Hz，散點表示BELLHOP理論值。數(shù)據(jù)處理結(jié)果與仿真結(jié)果一致，IFFT均無法分辨TDGF的前兩個到達時間，SBL在ΔL為30或120時，可以較好地估計第二個到達時間(包絡(luò)的波谷位置)。而當(dāng)ΔL=1(不進行預(yù)累積)或ΔL=720(全部預(yù)累積)時，估計結(jié)果分別出現(xiàn)了分辨率低與偏差大的問題。

圖5 不同預(yù)累積量時，利用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時兩個到達時延的差值：

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

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

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

3.4 結(jié)果

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

圖7 使用IFFT與SBL從海洋環(huán)境噪聲中提取的共計70 h的TDGF結(jié)果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。實驗結(jié)果證實了相比于傳統(tǒng)方法，預(yù)累積SBL提取的TDGF可以更準確地獲取較窄頻帶下的聲線到達結(jié)構(gòu)。

圖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 結(jié)論

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