高飛， 潘長(zhǎng)明， 孫磊,2

(1.海軍海洋測(cè)繪研究所, 天津 300061; 2.哈爾濱工程大學(xué) 水聲工程學(xué)院, 黑龍江 哈爾濱 150001)

Bayes匹配場(chǎng)地聲參數(shù)反演:多步退火Gibbs采樣算法

高飛1， 潘長(zhǎng)明1， 孫磊1,2

(1.海軍海洋測(cè)繪研究所, 天津 300061; 2.哈爾濱工程大學(xué) 水聲工程學(xué)院, 黑龍江 哈爾濱 150001)

為實(shí)現(xiàn)海洋地聲環(huán)境參數(shù)的快速高效反演，提出多步退火Gibbs(Multi-AG)采樣算法，可消除因參數(shù)搜索空間設(shè)置對(duì)參數(shù)反演結(jié)果的影響，并有效解決Bayes匹配場(chǎng)高維參數(shù)反演過程中常見的運(yùn)算量大、旁瓣高等問題。分析待反演地聲參數(shù)對(duì)匹配場(chǎng)處理器的敏感性，用以制定多步反演與退火方案，利用Gibbs采樣算法反演敏感性級(jí)別最高的參數(shù)，計(jì)算其均值并代入后續(xù)反演步驟，進(jìn)而采用退火Gibbs采樣算法逐步反演后續(xù)參數(shù)；利用數(shù)值仿真實(shí)驗(yàn)對(duì)比Metropolis-Hastings算法、Gibbs采樣算法、快速Gibbs采樣算法和Multi-AG采樣算法的反演效果。實(shí)驗(yàn)結(jié)果表明，與其他3種算法相比，Multi-AG采樣算法可通過最小的計(jì)算量得到均方差最小、精度最高的參數(shù)反演結(jié)果。

聲學(xué)； 多步退火Gibbs采樣算法； 地聲參數(shù)反演； Gibbs采樣； Bayes匹配場(chǎng)

0 引言

海底地聲參數(shù)信息是海洋水聲環(huán)境的重要參數(shù)組成之一，在海洋資源勘探、水下工程及聲納系統(tǒng)探測(cè)領(lǐng)域應(yīng)用廣泛，極具軍事價(jià)值和民用價(jià)值。當(dāng)前水聲調(diào)查通常包括海洋環(huán)境噪聲、混響、傳播損失、聲速剖面、水深及海面氣象等要素，且相關(guān)理論研究及調(diào)測(cè)設(shè)備發(fā)展相對(duì)成熟[1-2]。然而，受作業(yè)難度與技術(shù)要求制約，當(dāng)前淺地層地聲參數(shù)實(shí)測(cè)手段仍無法滿足大規(guī)模工程作業(yè)的要求，通過間接反演方法獲取地聲參數(shù)信息仍是當(dāng)前及未來較長(zhǎng)時(shí)間內(nèi)的主要方式。

海洋環(huán)境、聲源、垂線陣(VLA)聲信號(hào)是Bayes匹配場(chǎng)3要素[3]，由第1、第3項(xiàng)得到第2項(xiàng)為聲源定位，由后兩項(xiàng)得到第1項(xiàng)稱為海洋環(huán)境參數(shù)反演。當(dāng)前高維地聲參數(shù)Bayes匹配場(chǎng)反演技術(shù)主要面臨計(jì)算量大、旁瓣高等問題，也一直是研究改進(jìn)的焦點(diǎn)之一。

Bayes匹配場(chǎng)反演主要可從匹配處理器、隨機(jī)數(shù)據(jù)統(tǒng)計(jì)算法等方面進(jìn)行優(yōu)化。常見的Bartlett線性處理器[4]最早用于度量線列陣測(cè)量聲場(chǎng)與拷貝聲場(chǎng)的匹配程度，具有較好的穩(wěn)健性；自適應(yīng)場(chǎng)處理器[5]與多約束場(chǎng)處理器[6]通過權(quán)值矢量/后驗(yàn)概率密度(PDF)約束，有效融合了陣列信號(hào)處理的優(yōu)點(diǎn)，可抑制旁瓣與海洋環(huán)境不確定性，提供主瓣保護(hù)；Green函數(shù)處理器較好地處理了海洋環(huán)境不確定性，可削弱環(huán)境參數(shù)失配對(duì)匹配場(chǎng)反演精度的影響[7]。馬爾可夫鏈蒙特卡洛(MCMC)隨機(jī)統(tǒng)計(jì)思想[8]可反演得到參數(shù)PDF，但該方法在處理高維數(shù)據(jù)時(shí)計(jì)算量過大；遺傳算法[9-10]和最優(yōu)算法[11]在參數(shù)搜索空間中尋優(yōu)，可得到高精度反演結(jié)果；Gibbs采樣算法[12]是高維MCMC反演算法的一種改進(jìn)，可極大地減少運(yùn)算量；聚焦法[13]常用于聲源位置與參數(shù)共同反演，首先采用重要度采樣獲取局部最優(yōu)解，進(jìn)而通過Metropolis Gibbs采樣算法確定出聲源位置與參數(shù)PDF。

綜上所述可知，已有的研究成果未考慮各海洋環(huán)境參數(shù)對(duì)匹配場(chǎng)處理器的敏感性差異，導(dǎo)致敏感性較低的參數(shù)反演精度不高，且受搜索空間限制影響顯著。本文提出多步退火Gibbs(Multi-AG)采樣算法，以有效解決上述問題，進(jìn)而得到快速高精度反演結(jié)果。

1 Bayes匹配場(chǎng)反演理論

以m表示地聲參數(shù)信息(如沉積層密度、聲速、厚度等)，s為復(fù)數(shù)形式聲壓信號(hào)，則依據(jù)Bayes計(jì)算法可得

P(m|s)∝P(s|m)P(m)，

(1)

式中：P(m)為地聲參數(shù)m的先驗(yàn)概率密度，通常假設(shè)其在各自取值區(qū)間內(nèi)服從均勻分布，重要度采樣時(shí)與具體的密度分布函數(shù)g(m)相關(guān)；P(s|m)為已知m時(shí)s的條件概率密度，當(dāng)水聲環(huán)境參數(shù)m已知時(shí)，聲信號(hào)s可利用水聲數(shù)值模型計(jì)算得到；P(m|s)為PDF函數(shù)。

Bayes匹配場(chǎng)理論中，拷貝場(chǎng)聲信號(hào)s(在地聲參數(shù)m條件下，利用水聲數(shù)值模型計(jì)算得到)通常難以與觀測(cè)聲壓so完全相同，可利用似然函數(shù)描述兩者的相似程度L(so|m)，即

P(m|so)∝L(so|m)P(m).

(2)

于是Bayes地聲參數(shù)反演問題轉(zhuǎn)化為求取地聲參數(shù)m的PDFP(m|so)。反演結(jié)果通常以參數(shù)的一維邊緣概率密度表示：

P(mi|so)=∫δ(m′i-mi)P(m′|so)dm′，

(3)

(4)

(5)

(6)

(7)

式中：N為VLA陣元數(shù)；F為數(shù)據(jù)的頻率點(diǎn)數(shù)；vf與VLA信噪比(SNR)相關(guān)；Bf為常用的Bartlett處理器，

(8)

2 Multi-AG采樣算法

MCMC隨機(jī)模擬思想主要包括蒙特卡洛積分、馬爾可夫鏈及Metropolis-Hastings抽樣算法，Gibbs采樣算法可替代Metropolis-Hastings算法以提高在高維參數(shù)空間下計(jì)效率。

2.1 MCMC算法

2.1.1 蒙特卡洛積分

MCMC算法通過在各參數(shù)取值區(qū)間內(nèi)隨機(jī)均勻采樣，則任一組樣本參數(shù)mi的PDF為

(9)

式中：Q為樣本數(shù)。通常樣本參數(shù)mi的維度高(為多個(gè)地聲參數(shù)反演的聯(lián)合后驗(yàn)概率)，指示反演結(jié)果不夠直觀，可利用(3)式、(4)式、(5)式得到其一維邊緣概率密度、最大PDF、平均PDF.

當(dāng)隨機(jī)抽取樣本達(dá)到一定數(shù)量時(shí)，蒙特卡洛積分可得到地聲參數(shù)的一種無偏、非對(duì)稱估計(jì)結(jié)果。

2.1.2 馬爾可夫鏈

假設(shè)已生成的一連串隨機(jī)參數(shù)樣本{m0,m1,…,mt}，則t+1時(shí)刻的狀態(tài)概率為

P(mt+1|m0,m1,…,mt)=P(mt+1|mt).

(10)

(10)式表示了馬爾可夫鏈的基本原理，即任一時(shí)刻t+1的狀態(tài)轉(zhuǎn)移概率只與前一時(shí)刻t相關(guān)，而與已生成的樣本{m0,m1,…,mt}無關(guān)。

隨機(jī)抽取先驗(yàn)概率密度分布服從π(·)的地聲參數(shù)變量，生成易于實(shí)現(xiàn)且不可約遍歷鏈{mt}，使抽樣變量平穩(wěn)地服從π(·)分布。則在Bayes匹配場(chǎng)地聲參數(shù)反演過程中，隨著樣本數(shù)量的增加，馬爾可夫鏈逐漸忘記樣本初始化狀態(tài)m0，PDF的反演結(jié)果將最終服從某種分布π(·).

2.1.3 Metropolis-Hastings算法

Metropolis-Hastings算法是MCMC算法中隨機(jī)模擬與決定性算法相結(jié)合使用的一種常用解決方案，其算法的主要流程是：

1)初始化待反演地聲參數(shù)樣本m0；

α=min (1,exp (-ΔE))，

(11)

4)重復(fù)上述步驟，直至產(chǎn)生預(yù)設(shè)的T個(gè)樣本為止。

2.2Multi-AG采樣算法

Gibbs采樣將高維樣本化為一系列一維取樣，可避免Metropolis-Hastings算法在計(jì)算轉(zhuǎn)移概率α(通常α<1)偏小所導(dǎo)致的拒絕率過大的問題，進(jìn)而有效地提高運(yùn)行效率。

(12)

P(mi|so)=

∫…∫P(m|so)dm1…dmi-1dmi+1…dmP.

(13)

由(13)式可知，Metropolis-Hastings算法需經(jīng)P-1重積分方可得到參數(shù)的一維邊緣概率密度，而Gibbs采樣只需一重積分即可，即

(14)

Gibbs采樣按待反演參數(shù)維度逐個(gè)采樣判別，由于各參數(shù)對(duì)處理器敏感性存在差異，不同參數(shù)樣本拒絕率差別顯著。圖1為Gibbs采樣前180次參數(shù)判別跳轉(zhuǎn)軌跡，表1為Workshop’97淺海環(huán)境基準(zhǔn)模型參數(shù)[15]。受樣本拒絕率控制，c1的采樣判別次數(shù)遠(yuǎn)大于ρs、αs的采樣判別次數(shù)，導(dǎo)致前者的反演效果明顯優(yōu)于后二者。為解決各參數(shù)反演效果差異問題，通常的做法是擴(kuò)大采樣數(shù)量，這也極大地增加了計(jì)算量。

事實(shí)上，增加采樣數(shù)量對(duì)c1反演結(jié)果的影響較小，c1收斂速率快，在采樣反演前期即可完成收斂。而將所有參數(shù)共同進(jìn)行反演，c1既增加了計(jì)算量，也限制了其他參數(shù)的反演采樣次數(shù)。多步Gibbs采樣首先完成對(duì)敏感性最高的參數(shù)反演，利用反演結(jié)果計(jì)算其均值并視為常量，進(jìn)而逐級(jí)反演敏感性更低的參數(shù)，可極大地減少計(jì)算量。

圖1 Gibbs采樣參數(shù)跳轉(zhuǎn)軌跡Fig.1 Parameters skipping track of Gibbs sampling method

參數(shù)實(shí)際值聲源深度zs/m20聲源水平距離rs/km1水體深度D1/m100水體吸收系數(shù)αw/(dB·λ-1)1×10-4沉積層厚度D2/m30.7沉積層表層聲速c1/(m·s-1)1530.4沉積層底層聲速c2/(m·s-1)1604.1沉積層密度ρs/(kg·m-3)1500沉積層衰減系數(shù)αs/(dB·λ-1)0.2基底聲速c3/(m·s-1)1689基底密度ρb/(kg·m-3)1700基底衰減系數(shù)αb/(dB·λ-1)0.2

參數(shù)搜索空間對(duì)Gibbs采樣反演結(jié)果的影響較大，特別是在敏感性較低的參數(shù)反演時(shí)，不同的搜索空間得到的參數(shù)均值、均方差差異顯著。圖2為ρs一維PDF反演結(jié)果，分別設(shè)置參數(shù)采樣空間為[1 300 kg/m3,1 600 kg/m3](見圖2(a))、[1 400 kg/m3, 1 700 kg/m3](見圖2(b))時(shí)，Gibbs采樣算法反演均值(±均方差)分別為1 498.6 kg/m3±53.97 kg/m3、1 511.3 kg/m3±62.07 kg/m3，二者差異較大，且與實(shí)際值1 500 kg/m3偏離較遠(yuǎn)。此現(xiàn)象主要源于敏感度較低的參數(shù)對(duì)處理器影響微弱，造成其整個(gè)搜索空間內(nèi)接受率高，進(jìn)而導(dǎo)致其一維PDF較均勻所致。

圖2 參數(shù)ρs一維邊緣PDF分布Fig.2 Marginal posterior probability function of parameter ρs

模擬退火(SA)算法是概率算法的一種。該方法結(jié)合熱力學(xué)與統(tǒng)計(jì)學(xué)理論，將參數(shù)空間內(nèi)各點(diǎn)的概率視為分子動(dòng)能，并計(jì)算前后兩種狀態(tài)轉(zhuǎn)化的目標(biāo)函數(shù)差，進(jìn)而依據(jù)判定準(zhǔn)則判決是否接受該樣本。

定義退火Gibbs采樣算法的目標(biāo)函數(shù)差為

(15)

同樣的搜索空間，當(dāng)采用退火Gibbs采樣算法時(shí)，ρs的反演結(jié)果分別為1 501.1 kg/m3±28.13 kg/m3(見圖2(c))、1 500.7 kg/m3±28.31 kg/m3(見圖2(d))，二者間的微小差異源于隨機(jī)采樣所致，而與搜索空間無關(guān)。

Multi-AG采樣算法結(jié)合退火Gibbs算法，逐步完成高維地聲參數(shù)反演，Multi-AG采樣算法參數(shù)反演算法流程如圖3所示。圖3中，參數(shù)敏感性分析①以Bartlett處理器為量化依據(jù)，Multi-AG反演根據(jù)參數(shù)敏感性等級(jí)劃分為n步，第m(2≤m≤n)步將前m-1步中敏感性較高的參數(shù)反演均值作為常量輸入，即②和③，反演第1步利用Gibbs采樣算法，后續(xù)步驟使用退火Gibbs采樣算法。

圖3 Multi-AG采樣算法參數(shù)反演流程Fig.3 Flow chart of parameters inversion by multi-annealing Gibbs sampling method

3 地聲參數(shù)反演

水聲調(diào)查數(shù)據(jù)的系統(tǒng)性(缺實(shí)測(cè)地聲環(huán)境數(shù)據(jù))限制了對(duì)地聲參數(shù)反演結(jié)果的驗(yàn)證，實(shí)際海洋環(huán)境的不確定性降低了Bayes匹配場(chǎng)的反演精度。為更好地實(shí)現(xiàn)與檢驗(yàn)Multi-AG采樣算法的地聲參數(shù)反演效果，本文采用Workshop’97[15]提供的淺海環(huán)境基準(zhǔn)模型進(jìn)行地聲參數(shù)反演(見圖4)，該模型在水聲學(xué)理論研究中具有重要的參考價(jià)值[12,16]，利用其進(jìn)行反演得到的結(jié)果可與已有的研究成果進(jìn)行對(duì)比驗(yàn)證。

圖4 Workshop’97淺海環(huán)境基準(zhǔn)模型示意圖[15]Fig.4 Sketch of the benchmark of shallow water environment by Workshop’97[15]

[15]中測(cè)試樣例1，相關(guān)參數(shù)物理含義及取值如表1所示。Workshop’97淺海環(huán)境基準(zhǔn)模型為水平不變的兩層海底海洋環(huán)境，考慮沉積層聲速垂向變化，無限大基底層聲速恒等于1 700 m/s，聲速剖面同文獻(xiàn)[15]中的圖2.

VLA含20個(gè)垂向間距4 m的陣元，最淺單元位于3 m，無指向性聲源頻率為100 Hz. 利用Kraken簡(jiǎn)正波數(shù)值模型[17]計(jì)算聲場(chǎng)(其中模型代碼下載地址為：http:∥oalib.saic.com)，垂向計(jì)算點(diǎn)分辨率1 m，水平10 m，SNR為20 dB.

3.1 地聲參數(shù)敏感性等級(jí)劃分

表1中含8個(gè)地聲參數(shù)，分別為αs、ρs、c1、c2、D2、αb、ρb、c3，地聲參數(shù)反演過程中，聲場(chǎng)對(duì)各參數(shù)的敏感性不同，進(jìn)而導(dǎo)致Bartlett處理器變化尺度各異，Bayes匹配場(chǎng)對(duì)敏感性較低的參數(shù)反演效果較差，且增大了運(yùn)算量，故進(jìn)行地聲參數(shù)敏感性分析尤為重要。

基于單因子變量原則，分別計(jì)算各地聲參數(shù)變化時(shí)聲場(chǎng)與參數(shù)實(shí)際值下聲場(chǎng)1 km處的Bartlett處理器PDF，PDF分布特征如圖5所示。

圖5 各地聲參數(shù)變化對(duì)Bartlett處理器PDF分布Fig.5 Distribution of posteriori probability density of Bartlett mismatch function varying with geoacoustic parameters

分析易知：c1敏感性最高(見圖5(c))；c2、D2次之；αs、ρs、c3較??；αb、ρb的PDF幾乎服從均勻分布，對(duì)Bartlett影響可近似忽略。沉積層覆蓋于基底之上，與聲波相互作用顯著，其聲速大小影響聲阻抗造成的海底反射損失是底邊界對(duì)聲場(chǎng)作用的主要來源，故c1對(duì)Bartlett處理器的敏感性最高。

趙航芳等[18]利用Workshop’93典型淺海環(huán)境，分析了各水聲環(huán)境參量不確定性對(duì)Bayes匹配場(chǎng)性能的影響，其研究結(jié)論與本文符合較好。

圖6 各地聲參數(shù)敏感性變化曲線Fig.6 Sensitivity curves of geoacoustic parameters

經(jīng)分析，可將8個(gè)地聲參數(shù)劃分為4個(gè)敏感性等級(jí)。第1級(jí)(敏感性等級(jí)最高)含參數(shù)c1，比例在3%以內(nèi)；第2級(jí)含c2、D2，比例分別為9%、17.5%；第3級(jí)含αs、ρs、c3，比例分別為59%、75%、42%；第4級(jí)含αb、ρb，對(duì)Bartlett處理器幾乎無影響，比例大于99%.

參考已有的研究成果[10,12,19]反演得到的αb、ρb參數(shù)PDF在其搜索空間內(nèi)幾乎服從均勻分布，效果差，故下文不對(duì)αb、ρb進(jìn)行反演，按Workshop’97中的定義相應(yīng)地將其設(shè)置為常數(shù)。

3.2 退火方案

退火算法的目的在于提高參數(shù)收斂性能，消除敏感性較弱的參數(shù)因搜索空間設(shè)置所導(dǎo)致的反演誤差。但T太小或退火速率過快都會(huì)導(dǎo)致樣本拒絕率過大而無法采樣。

本文依據(jù)各參數(shù)的PDF及反演結(jié)果所預(yù)期的參數(shù)均方差范圍，分別制定退火方案。對(duì)于參數(shù)D2、αs、ρs、c3，T0分別取0.7 ℃、0.2 ℃、0.1 ℃、0.5 ℃. 通過數(shù)值實(shí)驗(yàn)并分析其PDF發(fā)現(xiàn)，隨著采樣數(shù)的增長(zhǎng)，在采樣前期(1/2倍總樣本數(shù))，控制T從1 ℃以總采樣數(shù)/20的速率降低至T0，可將參數(shù)反演均方差控制在相應(yīng)搜索空間寬度的15%左右，既提高了參數(shù)的反演精度及收斂性，又保證了退火Gibbs采樣算法較高的樣本接受概率。

3.3 Multi-AG法地聲參數(shù)反演

基于表1中的環(huán)境參數(shù)配置，通過Kraken模型計(jì)算1 km處對(duì)應(yīng)VLA深度的聲壓分布，并假定其為實(shí)測(cè)聲壓數(shù)據(jù)。依據(jù)3.1節(jié)和3.2節(jié)分析得到的參數(shù)敏感性等級(jí)劃分及退火方案，在兼顧水聲環(huán)境參數(shù)模型整體性的前提下，使地聲參數(shù)的采樣范圍最大化，設(shè)定參數(shù)c1、c2、D2、αs、ρs、c3的采樣空間分別為[1 500 m/s2,1 600 m/s2]、[1 500 m/s2,1 700 m/s2]、[10 m, 50 m]、[0.001 dB/λ, 0.500 dB/λ]、[1 300 kg/m3, 1 700 kg/m3]、[1 600 m/s2, 1 800 m/s2]，在采樣空間內(nèi)進(jìn)行隨機(jī)均勻采樣。

通過106次Metropolis-Hastings算法、105次Gibbs采樣算法、7×104次Multi-AG采樣算法(第1步、第2步、第3步分別計(jì)算104次、2×104次、4×104次)計(jì)算得到地聲參數(shù)PDF如圖7～圖9所示，圖中藍(lán)色實(shí)線為參數(shù)實(shí)際值位置。表2總結(jié)了Metropolis-Hastings算法、Gibbs采樣算法、快速Gibbs采樣(FGS)算法、Multi-AG采樣算法4種算法的反演結(jié)果及運(yùn)算量。

分析表2可知，105次Gibbs采樣算法反演得到的各參數(shù)均值及均方差與105次Metropolis-Hastings算法相當(dāng)，表明Gibbs采樣算法可極大地提高運(yùn)算效率。然而，各參數(shù)的反演效果隨各參數(shù)敏感性逐漸變差，敏感性最高的參數(shù)c1(見圖7(a)和圖8(a))反演結(jié)果分別為1 535 m/s±4.67 m/s、1 535 m/s±4.82 m/s，接近實(shí)際值1 530.4 m/s，且收斂性較好；敏感性低的參數(shù)αs(見圖7(d)和圖8(d))、ρs(見圖7(e)和圖8(e))、c3(見圖7(f)和圖8(f))的PDF幾乎接近均勻分布，偏離實(shí)際值較遠(yuǎn)，旁瓣高，且均方差過大。

圖7 Metropolis-Hastings方法地聲參數(shù)反演結(jié)果Fig.7 Inversion results of geoacoustic parameters by Metropolis-Hastings method

圖8 Gibbs采樣算法地聲參數(shù)反演結(jié)果Fig.8 Inversion results of geoacoustic parameters by Gibbs sampling method

對(duì)于多維參數(shù)共同反演，各參數(shù)對(duì)Bartlett處理器的敏感性差異較大，在參數(shù)采樣判別循環(huán)過程中，敏感度較高的參數(shù)循環(huán)判別次數(shù)多，限制了敏感度較低的參數(shù)隨機(jī)采樣次數(shù)。且較低的敏感性使得Bartlett處理器的效能降低，隨機(jī)誤判出現(xiàn)頻繁，無法有效地向參數(shù)真實(shí)值收斂，進(jìn)而出現(xiàn)近似均勻分布的PDF. 故可知，在敏感性較低的參數(shù)整個(gè)搜索空間內(nèi)都具有相當(dāng)?shù)臉颖窘邮芨怕?，?dǎo)致參數(shù)反演均值接近其搜索空間中值，使得參數(shù)搜索空間設(shè)置的差異影響到參數(shù)的反演結(jié)果。參數(shù)αs搜索空間為[0.001 dB/λ,0.500 dB/λ]，反演均值為0.280 dB/λ，接近搜索空間中值0.250 dB/λ，而與實(shí)際值0.200 dB/λ差異較大。

Multi-AG采樣算法分3步分別對(duì)6個(gè)參數(shù)逐步進(jìn)行反演，第1步循環(huán)104次，接受樣本622個(gè)，得到c1的反演結(jié)果為1 535±5.13 m/s，其一維邊緣概率密度如圖9(a)所示。對(duì)比106次的METROPOLIS-HASTINGS算法與105次的Gibbs采樣算法，三者反演效果差異微小，說明第1步Multi-AG算法可快速高精度地反演敏感性最高的參數(shù)c1.

利用第1步得到的c1PDF計(jì)算其均值(1 535 m/s)，并將c1視為海洋環(huán)境模型常量參數(shù)，進(jìn)而對(duì)αs、ρs、c2、D2、c3進(jìn)行退火Gibbs采樣反演。第2步采樣2×104次，接受1 553個(gè)樣本，可得到敏感性較低的參數(shù)c2(見圖9(b))、D2(見圖9(c))的PDF(均值±均方差)分別為1 597.3 m/s±13.5 m/s、30.1 m±3.2 m，為明顯的單峰PDF，較Metropolis-Hastings算法、Gibbs采樣算法更接近實(shí)際值。

圖9 多步退火Gibbs采樣算法地聲參數(shù)反演結(jié)果Fig.9 Inversion results of geoacoustic parameters by multi-step Gibbs sampling method

算法c1/(m·s-1)c2/(m·s-1)D2/mαs/(dB·λ-1)ρs/(kg·m-3)c3/(m·s-1)總樣本數(shù)接受樣本數(shù)實(shí)際值1530.41604.130.700.2015001689Metropolis-Hastings采樣算法1535.0±4.671592.0±13.528.12±8.90.28±0.1001452±101.01699±53.01×106132456Gibbs采樣算法1535.0±4.821591.0±11.729.46±7.10.28±0.0891454±91.31698±42.81×1058233FGS采樣算法[12]1534.0±7.001590.0±13.029.00±2.00.27±0.0761580±90.11687±11.01×105≈1×104Multi-AG采樣算法(第1步)1535.0±5.131585.0±19.725.27±10.20.33±0.1461444±99.51718±55.11×104622Multi-AG采樣算法(第2步)1535.0±5.131597.3±13.530.10±3.20.34±0.1521456±81.81715±41.62×1041553Multi-AG采樣算法(第3步)1535.0±5.131597.3±13.530.10±3.20.22±0.0371490±59.21692±10.64×1044980

依照第2步Multi-AG算法，對(duì)αs、ρs、c3進(jìn)行反演。第3步采樣4×104次，接受4 980個(gè)樣本，各參數(shù)的一維PDF如圖9(d)～圖9(f)所示。圖9中，αs、ρs、c3的PDF(均值±均方差)分別為0.22 dB/λ±0.037 dB/λ、1 490 kg/m3±59.2 kg/m3、1 692 m/s±10.6 m/s，單峰現(xiàn)象明顯，且搜索空間邊界附近的PDF值近似為0，有效地消除了旁瓣，表明參數(shù)反演精度顯著提高。

值得一提的是，對(duì)于多維地聲參數(shù)反演，Metropolis-Hastings算法、Gibbs采樣算法、Multi-AG 3種方法對(duì)敏感度最高的參數(shù)(c1)反演效果相近。

對(duì)比Multi-AG算法各步驟的總樣本數(shù)與接受樣本數(shù)，第2步(5個(gè)反演參數(shù))采樣數(shù)是第1步(6個(gè)反演參數(shù))的兩倍，接受樣本數(shù)大于其2倍；第1步、第3步(3個(gè)反演參數(shù))采樣數(shù)相同，第3步接受到的樣本數(shù)大于第2步?；谒晹?shù)值模型計(jì)算聲場(chǎng)用于反演地聲參數(shù)，是一個(gè)強(qiáng)非線性問題[20]，Stan等[12]研究了各參數(shù)間相關(guān)性對(duì)反演結(jié)果及運(yùn)算量的影響，指出不確定參數(shù)越多，反演所需計(jì)算量越大。本文通過敏感性分析，采用多步反演策略，減少了計(jì)算過程中的參數(shù)數(shù)目，提高了精度，較常規(guī)的Gibbs采樣算法進(jìn)一步減少了計(jì)算量。

退火算法可有效增強(qiáng)參數(shù)對(duì)Bartlett處理器的敏感性，提高其收斂性能。對(duì)于敏感性較高的參數(shù)c1、c2，無需退火，直接進(jìn)行Gibbs采樣判別即可得到較好的反演效果。從圖10(a)和圖10(b)采樣軌跡可看出，參數(shù)c1、c2接受的樣本收斂性較好，所接受的樣本值接近實(shí)際值。

Multi-AG第1步采樣含6個(gè)反演參數(shù)，敏感性弱的參數(shù)αs(見圖10(c))、ρs(見圖10(d))、c3(見圖10(e))接受的樣本軌跡分散。利用分布反演策略后，Multi-AG第3步只含3個(gè)參數(shù)，進(jìn)而對(duì)其進(jìn)行退火Gibbs采樣，隨者T的降低，接受的樣本逐漸向各自的實(shí)際值聚集，可有效提高反演精度并減小結(jié)果的均方差(見圖10(f)、圖10(g)和圖10(h))。

3.4 結(jié)果驗(yàn)證

Stan等[12,16]基于Workshop’97淺海環(huán)境基準(zhǔn)模型，利用FGS算法反演地聲參數(shù)，相關(guān)參數(shù)的反演結(jié)果如表2(引自文獻(xiàn)[12]中的表1)所示。對(duì)比分析參數(shù)實(shí)際值、FGS算法及本文的Multi-AG算法可知，Multi-AG算法經(jīng)過多步反演得到的參數(shù)αs、ρs、c2、D2反演精度更高。其中，F(xiàn)GS算法反演得到的αs、ρs均方差分別為0.076 dB/λ、90.1 kg/m3，Multi-AG算法反演得到的αs、ρs均方差分別為0.037 dB/λ、59.2 kg/m3，表明Multi-AG算法反演得到的參數(shù)收斂性更好。

4 結(jié)論

利用Gibbs采樣算法代替常規(guī)的Metropolis-Hastings算法進(jìn)行高維Bayes匹配場(chǎng)地聲參數(shù)反演，可在一定程度上減小運(yùn)算量。然而，各參數(shù)對(duì)反演處理器敏感性的不同，使得反演獲取各參數(shù)穩(wěn)定的一維PDF所需的計(jì)算量各異，且敏感性較小的參數(shù)搜索空間差異會(huì)在一定程度上影響參數(shù)的反演結(jié)果。

本文綜合Gibbs采樣算法、退火算法及地聲參數(shù)對(duì)Bartlett處理器的敏感性差異，提出Multi-AG采樣算法。通過采用Workshop’97淺海環(huán)境基準(zhǔn)模型，對(duì)比分析Metropolis-Hastings、Gibbs采樣算法、FGS算法、Multi-AG采樣算法4種算法地聲參數(shù)反演所需計(jì)算量、均值、均方差，結(jié)果表明：

1)Multi-AG參數(shù)反演所需計(jì)算量最小，說明分布反演方案進(jìn)一步提高了Gibbs采樣算法的反演效率。

2)Multi-AG參數(shù)反演均值與實(shí)測(cè)值最為接近，均方差最小，敏感性較低的參數(shù)反演改進(jìn)效果尤為明顯，證實(shí)了退火方案可有效提高參數(shù)反演精度。

本文Multi-AG算法用于地聲參數(shù)反演取得了較好效果，該算法也可用于Bayes匹配場(chǎng)其他海洋環(huán)境參數(shù)(水深、聲速剖面等)的反演。由于缺少實(shí)測(cè)地聲參數(shù)數(shù)據(jù)，本文將反演結(jié)果與已有的研究文獻(xiàn)[12]進(jìn)行了對(duì)比驗(yàn)證，這也是下一步將需要深化的研究方向。

參考文獻(xiàn)(References)

[1] 潘長(zhǎng)明, 高飛, 孫磊, 等. 淺海溫躍層對(duì)水聲傳播損失場(chǎng)的影響 [J]. 哈爾濱工程大學(xué)學(xué)報(bào), 2014, 35(4): 401-407. PAN Chang-ming, GAO Fei, SUN Lei, et al. The effects of shallow water thermocline on water acoustic transmission loss field[J]. Journal of Harbin Engineering University, 2014, 35(4): 401-407. (in Chinese)

[2] Peter F W, Matthew A D, James A M, et al. The North Pacific Acoustic Laboratory deep-water acoustic propagation experiments in the Philippine Sea[J]. Journal of the Acoustical Society of America, 2013, 134(6): 3359-3375.

[3] 李倩倩. 不確知海洋環(huán)境下的貝葉斯匹配場(chǎng)定位性能研究[D]. 北京: 中國(guó)科學(xué)院大學(xué), 2013: 1-2. LI Qian-qian. Source localization via Bayesian matched-field processing in an uncertain ocean environment[D].Beijing: University of Chinese Academy of Sciences, 2013: 1-2. (in Chinese)

[4] Bucker H P. Use of calculated sound fields and matched field detection to location sound sources in shallow water[J]. Journal of the Acoustical Society of America, 1976, 59(2): 368-372.

[5] Richardson A M, Nolte L W. A posteriori probability source localization in an uncertain sound speed, deep ocean environment[J]. Journal of the Acoustical Society of America, 1991, 89(5): 2280-2284.

[6] 王奇, 王英民, 茍艷妮. 不確定環(huán)境下后驗(yàn)概率約束的匹配場(chǎng)處理 [J]. 兵工學(xué)報(bào), 2014, 35(9): 1473-1480. WANG Qi, WANG Ying-min, GOU Yan-ni. Posterior probability constraint matched field processing with environmental uncertainty[J]. Acta Armamentarii, 2014, 35(9): 1473-1480. (in Chinese)

[7] Yann L G, Stan E D, Francois X S, et al. Bayesian source localization with uncertain Green’s function in an uncertain shallow water ocean[J]. Journal of the Acoustical Society of America, 2016, 139(3): 993-1004.

[8] Oh S H, Kwon B D. Geostatistical approach to Bayesian inversion of geophysical data: Markov chain Monte Carlo method[J]. Earth Planets Space, 2001, 53: 777-791.

[9] Perter G. Inversion of seismoacoustic data using genetic algorithms and a posteriori probability distributions[J]. Journal of the Acoustical Society of America, 1994, 95(2): 770-782.

[10] 李倩倩, 李整林, 張仁和. 不確知海洋環(huán)境下的貝葉斯聲源定位法 [J]. 聲學(xué)學(xué)報(bào), 2014, 39(5): 535-543. LI Qian-qian, LI Zheng-lin, ZHANG Ren-he. Bayesian localization in an uncertain ocean environment[J]. Acta Acustica, 2014, 39(5): 535-543. (in Chinese)

[11] Li Z L, Li F H. Geoacoustic inversion for sediments in the South China Sea based on a hybrid inversion scheme[J]. Chinese Journal of Oceanology and Limnology, 2010, 28(5): 990-995.

[12] Stan E D. Quantifying uncertainty in geoacoustic inversion I. A fast Gibbs sampler approach[J]. Journal of the Acoustical Society of America, 2002, 111(1): 129-142.

[13] Stan E D, Michael J W. Comparison of focalization and marginalization for Bayesian tracking in an uncertain ocean environment[J]. Journal of the Acoustical Society of America, 2009, 125(2): 717-722.

[14] Chen F H, Perter G, William S H. Validation of statistical estimation of transmission loss in the presence of geoacoustic inversion uncertainty[J]. Journal of the Acoustical Society of America, 2006, 120(4): 1932-1941.

[15] Tolstoy A, Chapman N R, Brooke G. Workshop’97: benchmarking for geoacoustic inversion in shallow water[J]. The Journal of Compute Acoustics, 1998, 6(4): 1-28.

[16] Stan E B, Peter L N. Quantifying uncertainty in geoacoustic inversionⅡ. Application to broadband, shallow-water data [J]. Journal of the Acoustical Society of America, 2002, 111(1): 143-159.

[17] Porter M B, Reiss E L. A numerical method for ocean-acoustic normal modes [J]. Journal of the Acoustical Society of America, 1984, 76(3): 244-252.

[18] 趙航芳，李建龍，宮先儀．不確實(shí)海洋中最小方差匹配場(chǎng)波束形成對(duì)環(huán)境參量失配的靈敏性分析 [J]．哈爾濱工程大學(xué)學(xué)報(bào), 2011, 32(2): 200-208． ZHAO Hang-fang, LI Jian-long， GONG Xian-yi. Sensitivity of minimum variance matched-field beam forming to an environmental parameter mismatch in an uncertain ocean channel[J]．Journal of Harbin Engineering University, 2011, 32(2): 200-208. (in Chinese)

[19] Liu Z, Sun C, Yang Y, et al. Robust source localization using predictable mode subspace in uncertain shallow water environment [C]∥OCEANS 2013 MTS/IEEE Conference. San Diego, CA, US: IEEE, 2013.

[20] Nattapol A, Zoi-Heleni M. Sequential filtering for dispersion for tracking and sediment sound speed inversion [J]. Journal of the Acoustical Society of America, 2014, 136(5): 2655-2674.

Geoacoustic Parameters Inversion of Bayes Matched-field: A Multi-annealing Gibbs Sampling Algorithm

GAO Fei1, PAN Chang-ming1, SUN Lei1,2

(1.Naval Institute of Hydrographic Surveying and Charting, Tianjin 300061, China; 2.College of Underwater Acoustic Engineering， Harbin Engineering University, Harbin 150001, Heilongjiang, China)

A multi-annealing Gibbs (multi-AG) sampling algorithm is developed to obtain a fast, accurate inversion result of ocean geoacoustic parameters. The proposed algorithm can deal well with huge computation load and high side lobe in multi-dimensional inversion of Bayes matched-field, and also eliminate the effects from the sampling bounds. The sensitivity of geoacoustic parameters to the matched-field processor is analyzed, which contributes to establish the multi-step inversion and annealing scheme. The Gibbs sampling algorithm is used to invert the highest sensitive parameters, which mean value is necessary to the following inversion steps. The inversion of remain parameters can be operated with annealing Gibbs sampling algorithm step by step. The inversion effects of Metropolis-Hastings, Gibbs, FGS, and multi-AG algorithms are compared through numerical experiment, and the research shows that multi-AG sampling algorithm can be used to obtain the inversion results with the smallest mean square deviation and the highest precision, while costing the least algorithm computation.

acoustics; multi-AG sampling algorithm; geoacoustic parameters inversion; Gibbs sampling; Bayes matched-field

2016-10-28

國(guó)家自然科學(xué)基金項(xiàng)目(41406004)

高飛(1988—)，男，助理工程師。E-mail:gfei88_lgdx@163.com； 潘長(zhǎng)明(1962—)，男，高級(jí)工程師。E-mail: pcming62@163.com

P733.23

A

1000-1093(2017)07-1385-10

10.3969/j.issn.1000-1093.2017.07.017