郝望 段睿? 楊坤德

1) (西北工業(yè)大學航海學院,西安 710072)

2) (西北工業(yè)大學海洋研究院,太倉 215400)

3) (西北工業(yè)大學,海洋聲學信息感知工業(yè)和信息化部重點實驗室,西安 710072)

大多數(shù)基于淺海簡正波模態(tài)頻散數(shù)據(jù)的地聲參數(shù)反演方法無法對深層底質(zhì)聲學參數(shù)進行可靠估計,究其原因是僅利用了簡正波水波頻散特征,忽略了與深層底質(zhì)聲學參數(shù)密切相關(guān)的底波頻散特征,因此,本文在分析了包含水波和底波的淺海寬帶數(shù)據(jù)的基礎(chǔ)上,結(jié)合底波頻散特征對深層底質(zhì)聲學參數(shù)變化更加敏感的物理特性,實現(xiàn)了基于完整簡正波頻散特性的貝葉斯地聲參數(shù)反演,并針對簡正波寬帶聲場模型計算復雜度較高的現(xiàn)實問題,利用變分貝葉斯蒙特卡羅方法的推斷優(yōu)勢,完成了未知參數(shù)的可靠估計和快速后驗分析.仿真和海上實驗結(jié)果表明: 聯(lián)合簡正波水波和底波頻散特征數(shù)據(jù)的貝葉斯地聲參數(shù)反演,不僅可以有效估計深層底質(zhì)聲學參數(shù),而且降低了其他相關(guān)環(huán)境參數(shù)的估計不確定性.

1 引言

淺海波導環(huán)境復雜多樣,其聲學特性參數(shù)主要包括水體和海底參數(shù),例如聲速和密度等,這些參數(shù)對于聲場預報和目標定位等研究有著重要意義.水體聲學參數(shù)可以通過直接測量的方式快速準確地獲得,但是對于海底聲學參數(shù),也稱地聲參數(shù),直接測量起來較為困難,需要耗費較大的人力和物力,且通常只能獲得實驗站位處的淺層測量結(jié)果,因此,利用水聲觀測數(shù)據(jù)進行地聲參數(shù)反演的研究就顯得尤為重要[1],地聲參數(shù)反演具有成本低、效率高和范圍廣等眾多優(yōu)勢,現(xiàn)有的地聲參數(shù)反演方法有傳統(tǒng)的匹配反演和線性擾動反演,以及近些年研究的熱點,如貝葉斯反演和神經(jīng)網(wǎng)絡(luò)等[2].

在淺海地聲參數(shù)反演研究中,經(jīng)常使用的水聲觀測數(shù)據(jù)有聲壓場數(shù)據(jù)[3,4]、海底反射數(shù)據(jù)[5,6]和模態(tài)數(shù)據(jù)[7-10]等,其中,利用模態(tài)頻散數(shù)據(jù)進行地聲參數(shù)反演的研究已經(jīng)發(fā)展較為成熟,該方法的主要特點是僅利用單個接收水聽器就可以實現(xiàn)相關(guān)環(huán)境參數(shù)的可靠估計,無需已知聲源級,例如: 文獻[7]提出利用Warping 變換對一定距離上的寬帶脈沖接收信號進行模態(tài)分離,基于分離的水波模態(tài)頻散特征數(shù)據(jù)進行地聲參數(shù)反演;另外,針對遠距離海底參數(shù)反演性能下降的問題,文獻[8]提出采用雙聲源模態(tài)相對到達時間差數(shù)據(jù)進行地聲參數(shù)反演,消除了距離累積帶來的誤差,實現(xiàn)了遠距離地聲參數(shù)的可靠估計.雖然,上述研究已經(jīng)充分證明利用淺海模態(tài)頻散數(shù)據(jù)可以有效地估計相關(guān)地聲參數(shù),但是不難發(fā)現(xiàn),幾乎所有的研究都是基于模態(tài)的水波頻散特征,而對于底波頻散特征的研究極少,所能利用的實驗數(shù)據(jù)也較少,導致反演結(jié)果僅能給出海底淺層地聲參數(shù)的有效估計,對于深層地聲參數(shù)基本無效,因此,本文聯(lián)合簡正波的水波和底波頻散特性進行地聲參數(shù)反演,一方面實現(xiàn)對深層地聲參數(shù)的可靠估計,另一方面降低其他模型輸入?yún)?shù)的反演不確定性.

值得注意的是,2004 年和2018 年,Lin 等[11]和Wan 等[12]分別提出利用包含Airy 相結(jié)構(gòu)的模態(tài)頻散特征數(shù)據(jù)對海底的壓縮波聲速進行反演,但是對于其他類型的環(huán)境參數(shù)文中沒有進行過多討論,除此之外,前者沒有進行參數(shù)反演的不確定性分析,后者采用效率較低的基于采樣策略的推斷方法進行了簡單的不確定性分析.隨著非線性貝葉斯近似推斷理論的發(fā)展,地聲參數(shù)反演不僅可以給出模型輸入?yún)?shù)的后驗估計結(jié)果,而且可以定量地描述參數(shù)估計的不確定性以及參數(shù)之間的相關(guān)性,對于大多數(shù)地聲參數(shù)反演問題,前向模型的巨大計算代價導致基于采樣策略的近似推斷方法無法在短時間內(nèi)完成參數(shù)的后驗分析,比如遺傳算法(genetic algorithm,GA)和馬爾科夫鏈蒙特卡羅(Markov chain Monto Carlo,MCMC)采樣方法等,相比之下,變分推斷(variational inference,VI)[13]作為另一類近似推斷方法,因其計算效率更高,且計算結(jié)果較為可靠,將有望代替MCMC 采樣,成為一種更加實用的后驗分析工具.本文使用一種數(shù)據(jù)高效的近似推斷方法,也稱變分貝葉斯蒙特卡羅(variational Bayesian Monte Carlo,VBMC)[14,15],完成地聲參數(shù)反演的后驗分析工作.

本文第2 節(jié)為基本理論,分別介紹了包含Airy相結(jié)構(gòu)的頻散曲線、貝葉斯地聲參數(shù)反演理論、VBMC 方法以及地聲參數(shù)反演流程;第3 節(jié)為仿真數(shù)據(jù)驗證;第4 節(jié)為實驗數(shù)據(jù)驗證;第5 節(jié)為結(jié)論.

2 基本理論

2.1 包含Airy 相結(jié)構(gòu)的頻散曲線

根據(jù)簡正波聲場理論,當聲源發(fā)射深度為zs,水聽器接收深度為zr,收發(fā)距離為r時,接收復聲壓譜Y(f) 可以表示為[16]

其中,f為信號頻率,A為常數(shù)因子,S(f) 為聲源的發(fā)射譜,M為總的傳播模態(tài)數(shù),ψ m和k rm分別為第m階簡正波模態(tài)的深度函數(shù)和水平波數(shù).

由(1)式可得,復聲壓譜Y是由各階簡正波疊加組成,對于第m階簡正波,不同頻率的信號能量具有不同的傳播速度,也稱群速度其表達式為

各階簡正波模態(tài)信號具有隨頻率變化的群速度,因此,當收發(fā)距離r足夠遠時,寬帶脈沖信號會表現(xiàn)出明顯的模內(nèi)頻散和模間頻散現(xiàn)象,假設(shè)在t0時刻發(fā)射一個寬帶脈沖信號,第m階模態(tài)信號到達接收點的絕對時間可以表示為

然而,在實際情況中,因為收發(fā)系統(tǒng)非嚴格同步,所以通常獲得第m階模態(tài)信號的相對到達時間,本文稱為模態(tài)頻散曲線,可以表示為

2.2 貝葉斯地聲參數(shù)反演理論

在貝葉斯反演框架內(nèi),未知的模型輸入?yún)?shù)被描述為多維隨機變量,并從后驗概率的角度對其進行估計和不確定性分析.結(jié)合本文的具體問題,采用隨機變量θ表示未知的環(huán)境參數(shù),包含部分水體參數(shù)和地聲參數(shù),d=[d1,d2,···,dM] 表示頻散曲線的觀測向量集合,其中,d m(m=1,2,···,M) 為第m階模態(tài)的觀測結(jié)果,頻點數(shù)為N m,通常情況下,d m可以表示為正演模型預測項g m(θ) 和誤差項εm之和,即

gm(θ) 由(5)式給出,ε m包含模型計算誤差和實際測量誤差.

貝葉斯地聲參數(shù)反演的核心是貝葉斯公式,隨機變量θ的后驗概率密度函數(shù)p(θ|d) 可以寫作[3]:

N(0,Cm),因此,似然函數(shù)L(θ) 可以表示為式中,p(θ) 表示先驗概率密度函數(shù),p(d) 為歸一化因子,與參數(shù)θ無關(guān),似然函數(shù)p(d|θ) 或L(θ) 表示誤差項的統(tǒng)計分布,本文假設(shè)不同階模態(tài)的頻散曲線數(shù)據(jù)誤差ε m之間相互獨立,服從高斯分布

其中,?表示向量轉(zhuǎn)置,C m為N m×Nm維的數(shù)據(jù)誤差協(xié)方差矩陣,通常假設(shè)在地聲參數(shù)反演中,可以利用多次迭代的數(shù)據(jù)殘差分析計算得到[17].

獲得環(huán)境參數(shù)的后驗概率密度函數(shù)p(θ|d) 后,其他統(tǒng)計結(jié)果可以利用(10)—(12)式計算得到,包括一維邊緣后驗概率密度函數(shù)p(θi|d)、后驗均值估計θ0和參數(shù)后驗協(xié)方差矩陣D:

其中,δ為狄拉克函數(shù),除上式外,還可以采用最高概率密度可信區(qū)間來定量描述參數(shù)的不確定性,例如,95%可信區(qū)間表示包含邊緣概率密度函數(shù)95%覆蓋面積的最小寬度區(qū)間.

2.3 變分貝葉斯蒙特卡羅方法

對于非線性貝葉斯推斷問題,通常采用MCMC或VI 方法近似求解,當前向模型的計算復雜度較高時,相比于MCMC 方法(通常需要數(shù)以萬計次的前向模型計算),基于尋優(yōu)策略的VI 方法可以更加高效地完成后驗分析工作,本文不再贅述其基本理論,感興趣的讀者可以參考相關(guān)文獻[13],這里僅給出散度(Kullback-Leibler divergence,KL)和證據(jù)下界(evidence lower bound,ELBO)的計算表達式,分別為

其中,q ?(θ) 為參數(shù)化的分布族,也稱變分族,?表示變分參數(shù),E ?表示關(guān)于q?求期望,p(d,θ) 為聯(lián)合概率密度函數(shù).

VI 方法的核心思想是利用優(yōu)化算法求出使ELBO 取得最大值時的變分參數(shù)?*,此時后驗概率密度函數(shù)p(θ|d) 可以被近似為為了計算(14)式中非線性函數(shù)的復雜期望,Acerbi[14]在2018 年提出了一種數(shù)據(jù)高效的近似推斷方法,也稱VBMC 方法,該方法將參數(shù)化的分布族q ?(θ) 建立為高斯混合模型(Gaussian mixed model,GMM),并使用高斯過程模型(Gaussian process,GP)擬合 l ogp(d,θ),同時,將基于主動采樣的貝葉斯正交(Bayesian quadrature and active sampling)以及蒙特卡羅采樣(Monte Carlo sampling)和VI 方法相結(jié)合,利用隨機梯度下降法對證據(jù)下界函數(shù)進行尋優(yōu),使其非常適合具有復雜似然函數(shù)的非線性貝葉斯推斷問題,大量仿真和實驗結(jié)果表明,該方法使用較小的計算代價可以獲得和MCMC 類方法幾乎相同的后驗分析準確度[15],主要步驟如表1所示.

表1 VBMC 算法步驟Table 1.Steps of the VBMC algorithm.

2.4 地聲參數(shù)反演流程

地聲參數(shù)反演流程如圖1 所示,假設(shè)模型的輸入?yún)?shù)已知,利用Kraken 聲場模型[18]計算出模態(tài)群速度后,代入(5)式求出模態(tài)頻散曲線,即為地聲參數(shù)反演中的正演模型g m(θ) .對于實際接收信號,利用Warping 變換[19]和同步壓縮小波變換[20]分別對水波和底波信號的頻散曲線進行分離和提取,即得到地聲參數(shù)反演中的觀測數(shù)據(jù)d m.

在貝葉斯反演框架內(nèi)(圖1 中虛線框所示),首先,設(shè)定待反演參數(shù)θ的先驗區(qū)間,假設(shè)在區(qū)間上服從均勻分布;然后,基于正演模型g m(θ) 和觀測數(shù)據(jù)d m建立貝葉斯反演模型,利用多次迭代的數(shù)據(jù)殘差分析方法估計誤差協(xié)方差矩陣C m,圖中θML為參數(shù)的最大似然估計值,在本文中使用自適應單純形模擬退火算法求得[21];最后,基于穩(wěn)定的Cm估計結(jié)果,采用VBMC 方法推斷后驗概率密度函數(shù)p(θ|d),完成參數(shù)后驗分析.在實際應用中,對于最終獲得的反演結(jié)果,通常還需要采用多種直接或間接的方式對其進行有效性驗證,例如原位測量對比和傳播損失對比等.值得注意的是,變分推斷類方法對初始狀態(tài)的選取較為敏感,因此,本文將θML作為VBMC 算法的初始值.

圖1 地聲參數(shù)反演流程圖Fig.1.Flow chart of geoacoustic inversion.

3 仿真數(shù)據(jù)驗證

3.1 仿真環(huán)境

圖2 給出了一個具有3 層水平分層介質(zhì)的波導模型,包括水體層、沉積層和半空間基底層: 在水體層中,海水聲速cw和海深dw均為常數(shù);在沉積層中,壓縮波聲速隨著沉積厚度的增大線性變化,在 [ 0,h1] 的 范圍內(nèi)從c1U變化到c1L,其密度ρ1為一個常數(shù);在基底半空間中,壓縮波聲速c2和密度ρ2均為常數(shù),值得注意的是,對于大多數(shù)淺海粉砂質(zhì)砂或砂質(zhì)粉砂海底,忽略剪切波的影響是合理的,所有的模型輸入?yún)?shù)總結(jié)如表2 所示.

圖2 3 層水平分層介質(zhì)波導模型Fig.2.Three-layer waveguide model.

將表2 中的仿真參數(shù)值代入Kraken 聲場模型,計算1 階模態(tài)群速度曲線,如圖3 所示,圖中的圓圈給出Airy 相頻率的位置,從整體上看,群速度曲線被Airy 相頻率分成了左右兩部分,分別稱為簡正波傳播中的底波和水波部分.當h1,c1L,c2和ρ2取不同值時,群速度曲線在小于Airy 相頻率的部分變化明顯,而在大于Airy 相頻率的部分基本不變,說明底波的群速度特征對深層的地聲參數(shù)變化較為敏感,充分利用這一特性,可以實現(xiàn)相關(guān)地聲參數(shù)的可靠反演,另外,從圖3 的結(jié)果可以看出,h1和c1L的變化會引起Airy 相頻率的變化,而c2和ρ2的變化對Airy 相頻率無影響.

圖3 1 階模態(tài)群速度曲線 (a)不同的 h1 ;(b)不同的 c1L ;(c)不同的 c2 ;(d)不同的 ρ2,圓圈給出了艾里相頻率的位置Fig.3.Group velocity curves of Mode 1: (a) Different values of h1 ;(b) different values of c1L ;(c) different values of c2 ;(d) different values of ρ2 .The circles indicate the Airy phase frequencies.

表2 反演參數(shù)空間和先驗區(qū)間Table 2.Parameter spaces and prior bounds for inversion.

3.2 底波頻散曲線提取誤差分析

在仿真接收信號的頻散特性時,假設(shè)聲源深度zs=6 m,接收深度zr=11.5 m,收發(fā)距離r=6.24 km,海底的衰減系數(shù)為0 d B/m,在不考慮環(huán)境噪聲的情況下,寬帶脈沖信號的接收時域波形和時頻分析結(jié)果如圖4 所示,在0—500 Hz 頻帶范圍內(nèi),接收信號包含前3 階簡正波模態(tài),其中,1 階和3 階模態(tài)信號能量較強,2 階模態(tài)信號能量較弱.利用圖1 所示的方法對1 階和3 階模態(tài)頻散曲線進行分離和提取,結(jié)果如圖4(b)中的“×”符號所示,可以看出,在理想情況下,實驗提取的頻散曲線與時頻譜圖顯示的頻散特征基本一致,并且與理論計算出的頻散曲線(圖中的實線)基本吻合.

圖4 仿真信號 (a)歸一化時域波形;(b)時頻圖,其中1 階模態(tài)的底波、水波和Airy 相結(jié)構(gòu)清晰可辯Fig.4.Simulation signal: (a) Normalized time domain waveform;(b) time-frequency diagram.The ground wave,water wave,and Airy phase component of Mode 1 are well-defined.

在實際中,底波信號的幅度通常遠小于水波信號的幅度,因此,實驗中聲源級大小、沉積物吸收以及環(huán)境噪聲等因素都會影響底波的觀測,導致其頻散曲線的提取存在誤差,綜合考慮上述因素,在仿真的底波信號中加入不同量級的高斯白噪聲后,進行頻散曲線的提取,分析不同信噪比(signal to noise ratio,SNR)條件下的提取誤差,結(jié)果如圖5所示,圖5(a)對比了不同信噪比條件下的單次提取結(jié)果和理論模型計算結(jié)果;假設(shè)提取誤差在不同頻點上服從獨立的零均值高斯分布,圖5(b)給出了不同信噪比條件下誤差標準差的估計結(jié)果,可以看出,隨著信噪比逐漸降低,誤差標準差逐漸增大.

圖5 不同信噪比條件下底波頻散曲線的單次提取結(jié)果(a)和提取誤差標準差(b)Fig.5.Single extraction results (a) and standard deviations of extraction errors (b) of ground wave dispersion curve under different SNRs.

3.3 地聲參數(shù)反演仿真驗證

為了討論引入簡正波的底波頻散特性對地聲參數(shù)反演的影響,本文針對兩種數(shù)據(jù)條件,分別進行仿真研究: 數(shù)據(jù)條件1(Data 1),包括1 階模態(tài)的底波和1 階,3 階模態(tài)的水波頻散特征數(shù)據(jù);數(shù)據(jù)條件2(Data 2),包括1 階,3 階模態(tài)的水波頻散特征數(shù)據(jù),特別地,對于數(shù)據(jù)條件1,考慮不同的底波信噪比情況,即SNR=5,0,—5 和—10 dB.

根據(jù)圖1 所示的流程進行地聲參數(shù)反演,部分深層參數(shù)的一維邊緣后驗概率密度函數(shù)如圖6 所示,從圖6(a)—(d)分別為h1,c1L,c2,ρ2,圖中的點劃線表示參數(shù)真值,可以發(fā)現(xiàn): 1)當?shù)撞ㄐ盘栃旁氡容^高時(SNR=5 dB,0 dB),數(shù)據(jù)條件1 的邊緣后驗概率密度函數(shù)較為尖銳,即反演的不確定性較小,且最大后驗概率估計值(概率密度函數(shù)的最大值點)與參數(shù)真值較為接近,相比之下,數(shù)據(jù)條件2 的邊緣后驗概率密度函數(shù)非常平坦,幾乎等于先驗概率密度函數(shù),因此,無法對該參數(shù)進行可靠估計;2)隨著底波信噪比的降低,其頻散曲線的提取誤差標準差逐漸增大(如圖5(b)所示),包含的有用信息越來越少,數(shù)據(jù)條件1 的地聲參數(shù)反演效果將逐漸退化為數(shù)據(jù)條件2 的地聲參數(shù)反演效果.總之,當?shù)撞ㄐ盘柕男旁氡茸銐蚋邥r,基于簡正波水波和底波完整頻散特性的貝葉斯地聲參數(shù)反演可以準確估計出深層地聲參數(shù).

圖6 兩種數(shù)據(jù)條件下部分參數(shù)的一維邊緣后驗概率密度函數(shù) (a) h1 ;(b) c1L ;(c) c2 ;(d)ρ2Fig.6.1 D marginal posterior probability densities of some parameters for Data 1 and Data 2: (a) h1 ;(b) c1L ;(c) c2 ;(d) ρ2 .

4 實驗數(shù)據(jù)驗證

4.1 實驗環(huán)境

圖7 給出某次渤海實驗的實驗地點和實測的水體聲速剖面(sound speed profile,SSP),在圖7(a)中,符號▲表示11 號自容式水聽器的具體位置,深度傳感器顯示其布放深度約為11.5 m,其靈敏度約為-170 dB re 1 V/μPa,采集系統(tǒng)的采樣率為30 kHz,實驗過程中,寬帶脈沖聲源在圖中★符號所示位置處發(fā)射信號,發(fā)射深度約為6 m,發(fā)射距離約為6.07 km,利用溫鹽深儀對水體聲速剖面進行測量,結(jié)果如圖7(b)所示,測點處的海深約為20 m,海面附近的聲速約為1512.5 m/s,海底附近的聲速約為1511 m/s,海水混合程度較高,介質(zhì)較均勻.

圖7 實驗描述 (a)實驗地點;(b)水體聲速剖面Fig.7.Description of the experiment: (a) Experimental site;(b) sound speed profile in water.

4.2 水波和底波頻散曲線提取

自容式水聽器接收信號的時域波形和時頻分析結(jié)果如圖8 所示,信號包含兩階簡正波模態(tài),通過仿真對比可以確定分別為1 階和3 階模態(tài),其中,1 階模態(tài)包含水波和底波兩組成分,而3 階模態(tài)僅包含水波成分,1 階模態(tài)的底波成分到達接收點的時間最早,準正弦的底波信號能量主要集中在35—74 Hz 之間,大幅值的水波信號緊隨其后,與底波信號相比,水波信號的波形更加尖銳,頻散特征更弱,1 階模態(tài)信號中到達時間最晚的部分(圖8(b)中方框區(qū)域)即為Airy 相結(jié)構(gòu),對應的Airy 相頻率約為74 Hz(點劃線).實驗中,接收信號的信噪比較高,有利于進行模態(tài)頻散特性分析.

利用Warping 變換和同步壓縮小波變換對水波信號和底波信號的頻散曲線分別進行提取,結(jié)果如圖8(b)中的“×”符號所示,頻散曲線和時頻圖中的頻散特征一致性較好,值得注意的是,寬帶聲源信號在頻譜上存在弱干涉,導致接收信號的模態(tài)能量在感興趣的頻段范圍內(nèi)出現(xiàn)強弱交替的現(xiàn)象,但是,其對頻散特征的影響較小,實際處理時可以忽略.

圖8 實驗信號 (a)歸一化時域波形;(b)時頻圖,其中1 階模態(tài)的底波、水波和Airy 相結(jié)構(gòu)清晰可辯Fig.8.Experimental signal: (a) Normalized time domain waveform;(b) time-frequency diagram.The ground wave,water wave,and Airy phase component of Mode 1 are well-defined.

4.3 地聲參數(shù)反演

與仿真研究相同,針對兩種實驗數(shù)據(jù)條件分別進行地聲參數(shù)反演,如圖9 所示,圖中使用誤差條給出了實驗數(shù)據(jù)誤差標準差σ m的估計結(jié)果.

圖9 頻散曲線的實測結(jié)果和后驗預測結(jié)果 (a)數(shù)據(jù)條件1,水波和底波頻散數(shù)據(jù),插圖為低頻部分的局部放大;(b)數(shù)據(jù)條件2,水波頻散數(shù)據(jù)Fig.9.Measurements and posterior predictions of dispersion curves: (a) Data 1,water waves and ground waves,the enlarged view of the low frequency part is inset;(b) Data 2,water waves.

待反演參數(shù)的邊緣后驗概率密度函數(shù)如圖10所示,對角線圖形給出了一維估計結(jié)果,實線和虛線分別表示數(shù)據(jù)條件1 和條件2,在先驗區(qū)間內(nèi),不同數(shù)據(jù)條件下的10 個參數(shù)結(jié)果均表現(xiàn)為單峰結(jié)構(gòu),但是,因為引入了底波頻散特征,數(shù)據(jù)條件1 的參數(shù)估計不確定性(概率密度函數(shù)的分布范圍)明顯小于數(shù)據(jù)條件2,特別是對于深層地聲參數(shù),這一結(jié)果尤為明顯,除此之外,聯(lián)合水波和底波完整頻散特性進行地聲參數(shù)反演,可以準確地分辨沉積層和基底半空間層,即h1的估計不確定性非常小,綜上所述,底波頻散特征數(shù)據(jù)可以提供更多關(guān)于海底聲學特性的信息,使相關(guān)參數(shù)的反演更加準確.

圖10 待反演參數(shù)的一維(數(shù)據(jù)條件1 和條件2)和二維(數(shù)據(jù)條件1)邊緣后驗概率密度函數(shù)Fig.10.1D (Data 1 and Data 2) and 2D (Data 1) marginal posterior probability densities of unknown parameters.

利用待反演參數(shù)的后驗估計結(jié)果對模態(tài)頻散曲線進行預測,如圖9 所示,虛線給出后驗均值樣本的預測值,細實線給出100 個后驗隨機樣本的預測值,在不確定性范圍內(nèi),預測結(jié)果和實驗數(shù)據(jù)的一致性較好,充分說明圖2 所示的淺海波導模型可以有效模擬此次海上實驗的真實環(huán)境.

針對數(shù)據(jù)條件1 的地聲參數(shù)反演給出了待反演參數(shù)的可靠估計,后面重點對其進行討論.參數(shù)的后驗統(tǒng)計結(jié)果總結(jié)在表3 中,包括均值、方差和95%可信區(qū)間,其中,沉積層參數(shù)結(jié)果表明,實驗海域的海底覆蓋一層中高聲速沉積物,平均厚度約為(46.65±3.77) m,上部的壓縮波聲速約為(1675.89±9.40) m/s,下部的壓縮波聲速約為(1673.52±8.04) m/s,平均密度約為(1.52±0.09) g/cm3,沉積層中介質(zhì)的低頻聲學特性隨著沉積厚度的增大變化不大.在基底層中,壓縮波聲速約為(1765.03±15.78)m/s,其95%的可信區(qū)間為[1731.16,1789.77] m/s,密度約為(2.12±0.19) g/cm3,其95%可信區(qū)間為[1.75,2.44] g/cm3.本文假設(shè)環(huán)境參數(shù)與距離無關(guān),因此海水深度dw和聲速cw代表了實驗海域的平均估計結(jié)果,分別為(20.69±0.45) m (實測結(jié)果為19—22 m)和(1510.98±4.13) m/s (實測結(jié)果為1511—1513 m/s),收發(fā)距離r的估計結(jié)果 為(6.24±0.3) km(實測結(jié)果為6.07 km),時間因子δt作為一個輔助變量,在反演過程中快速收斂,并且估計不 確定性較小,其均 值為0.280 s,方差為0.002 s.

表3 針對數(shù)據(jù)條件1 的后驗統(tǒng)計結(jié)果Table 3.Posterior statistical results for Data 1.

待反演參數(shù)之間的相關(guān)性可以通過二維邊緣后驗概率密度函數(shù)(圖10 中的非對角線圖形)和歸一化后驗參數(shù)協(xié)方差矩陣(圖11 中的熱圖)進行描述,結(jié)果顯示,不同參數(shù)之間存在一定的相關(guān)性,例如: 水平距離r和平均水深dw之間存在明顯的正相關(guān);水平距離r和時間因子δt之間存在明顯的負相關(guān),因此,當其中一個參數(shù)的估計出現(xiàn)偏差時,與之相關(guān)的其他參數(shù)的估計也會出現(xiàn)偏差.

圖11 針對數(shù)據(jù)條件1 的歸一化后驗參數(shù)協(xié)方差矩陣Fig.11.Normalized posterior covariance matrix of unknown parameters for Data 1.

圖12(a),(b)分別給出VBMC 推斷過程中ELBO 和KL 散度的收斂情況,可以看出,在經(jīng)過約80 次迭代之后,該方法已經(jīng)基本收斂,此時前向模型的計算次數(shù)僅為400 余次,與MCMC 方法相比,計算次數(shù)顯著減少,雖然VBMC 推斷還包括有梯度下降算法等額外的計算損耗,但是其總的后驗分析時間是遠小于MCMC 方法的.利用VBMC 方法求出近似的變分后驗概率密度函數(shù)后,通過簡單的采樣算法獲得其有效樣本,并進行統(tǒng)計分析.

圖12 針對數(shù)據(jù)條件1 的VBMC 收斂情況 (a)證據(jù)下界ELBO;(b) KL 散度Fig.12.Convergences of the VBMC method for Data 1:(a) ELBO;(b) KL divergence.

進一步討論海底的分層情況,針對數(shù)據(jù)條件1,假設(shè)在沉積層和基底層之間還存在一個中間層,該層的厚度、壓縮波聲速和密度均為未知常數(shù),同樣,利用2.4 節(jié)所述方法,可以獲得如圖13(b)所示的波導環(huán)境聲學參數(shù)剖面,包括壓縮波聲速(左)和密度(右),圖中粗實線表示后驗均值結(jié)果,細實線表示若干個后驗隨機結(jié)果,圖13(a)給出沒有中間層的波導環(huán)境聲學參數(shù)剖面,對比發(fā)現(xiàn),兩種海底分層結(jié)構(gòu)下的波導環(huán)境幾乎相同,海水的實測聲速剖面(虛線)和平均估計結(jié)果(實線)均擬合較好,海底中間層的出現(xiàn)等效于在50 m 深度附近將沉積層分成上下兩部分,除了下部的密度外,兩部分的底質(zhì)聲學特性基本一致.

圖13 針對數(shù)據(jù)條件1 的壓縮波聲速和密度剖面后驗估計結(jié)果 (a) 2 層海底模型;(b) 3 層海底模型Fig.13.Posteriori estimates of compressed-wave sound speed profiles and density profiles for Data 1: (a) Two-layer seabed model;(b) three-layer seabed model.

4.4 反演結(jié)果驗證

將表3 給出的參數(shù)均值代入Kraken 聲場模型,計算模態(tài)群速度,結(jié)合另一組已知收發(fā)距離的寬帶脈沖信號進行頻散曲線預測,圖14 給出了不同距離上接收信號的時頻分析結(jié)果,從圖14(a)—(c)分別為3.42,5.35 和6.51 km,實線為頻散曲線的預測結(jié)果,兩者的一致性較好,充分證明上述地聲參數(shù)反演結(jié)果的有效性和可靠性.

圖14 另一組寬帶脈沖信號時頻圖和頻散曲線預測結(jié)果(實線) (a)接收距離3.42 km;(b)接收距離5.35 km;(c)接收距離6.51 kmFig.14.Time-frequency diagrams and dispersion curve predictions (solid lines) of other broadband pulse signals: (a) The range is 3.42 km;(b) the range is 5.35 km;(c) the range is 6.51 km.

歷史調(diào)查和研究顯示,實驗海域的海深在15—25 m 之間,海底沉積物多為粉砂質(zhì)砂或砂質(zhì)粉砂[22,23],平均粒徑的變化范圍為[2Φ,5Φ],文獻[10,24]給出了相似淺海環(huán)境下粉砂質(zhì)沉積物低頻衰減系數(shù)α的取值,當頻率為200 Hz 時,α ≈0.02 dB/m,當頻率為315 Hz 時,α ≈0.06 dB/m .聲源深度為6 m,接收深度為11.5 m,結(jié)合其他地聲參數(shù)反演結(jié)果進行傳播損失預報,如圖15 所示,傳播損失的距離變化范圍為[0,10] km,圖中圓圈為實測傳播損失,距離分別為3.42,5.35 和6.51 km,可以看出,利用反演得到的地聲參數(shù)(圖中實線)可以較好地預測實驗海域的傳播損失.

圖15 傳播損失的理論計算結(jié)果和實測結(jié)果比較,聲源深度6 m,接收深度11.5 m (a)頻率200 Hz;(b)頻率315 HzFig.15.Comparisons of the theoretical and experimental transmission loss at source depth 6 m and receiver depth 11.5 m: (a) The frequency is 200 Hz;(b) the frequency is 315 Hz.

5 結(jié)論

針對基于簡正波水波頻散特性的地聲參數(shù)反演方法無法有效估計深層地聲參數(shù)的現(xiàn)實問題,本文聯(lián)合底波和水波頻散特性進行貝葉斯地聲參數(shù)反演,創(chuàng)新地采用變分貝葉斯蒙特卡羅算法對反演結(jié)果進行快速后驗分析,不僅實現(xiàn)了深層地聲參數(shù)的可靠估計,而且大大降低了淺層地聲參數(shù)和水體參數(shù)的估計不確定性,仿真和實驗結(jié)果顯示: 待反演參數(shù)的邊緣后驗概率密度函數(shù)表現(xiàn)為先驗區(qū)間內(nèi)的尖峰結(jié)構(gòu),頻散曲線的預測結(jié)果與測量結(jié)果擬合較好,說明通過此方法反演出的波導環(huán)境模型比較可靠.在實驗數(shù)據(jù)驗證中,收發(fā)距離和海水平均聲速兩個參數(shù)的估計結(jié)果與實際測量結(jié)果基本一致,反演得到的地聲參數(shù)不但可以用于解釋本組實驗數(shù)據(jù),而且可以較好地用于不同距離上另一組實驗數(shù)據(jù)的模態(tài)頻散特征預測,除此之外,傳播損失的準確預報也充分證明了聯(lián)合反演方法的有效性.

在本文的地聲參數(shù)反演研究中,假設(shè)海洋環(huán)境參數(shù)不隨距離發(fā)生變化,然而,實際的淺海環(huán)境復雜多樣,在水平和垂直方向上均表現(xiàn)出非均勻性,因此,如何實現(xiàn)距離相關(guān)淺海波導環(huán)境下地聲參數(shù)的可靠反演是下一步需要研究的工作.