地聲參數(shù)及傳播損失不確定性估計(jì)與建模?

過(guò)武宏　笪良龍　趙建昕

（海軍潛艇學(xué)院　青島　266071）

地聲參數(shù)及傳播損失不確定性估計(jì)與建模?

過(guò)武宏?笪良龍趙建昕

（海軍潛艇學(xué)院青島266071）

地聲參數(shù)的不確定性對(duì)水聲傳播具有重要的影響。通過(guò)貝葉斯理論建立水聲環(huán)境不確定性推理模型，理論推導(dǎo)了地聲參數(shù)的似然函數(shù)以及地聲參數(shù)和傳播損失的后驗(yàn)概率密度，并采用MCMC（Markov Chain Monte Carlo）進(jìn)行了仿真計(jì)算，給出了地聲參數(shù)的二維后驗(yàn)聯(lián)合概率密度和一維邊緣概率密度，在此基礎(chǔ)上對(duì)傳播損失的不確定性進(jìn)行了估計(jì)，得到了傳播損失80%的可信區(qū)間。仿真和實(shí)驗(yàn)結(jié)果表明，該方法適用于地聲參數(shù)反演和不確定性估計(jì)，并能獲取因地聲參數(shù)不確定性導(dǎo)致的傳播損失不確定性估計(jì)。

貝葉斯，后驗(yàn)概率密度，地聲參數(shù)，馬爾可夫鏈蒙特卡洛

1　引言

海洋水體環(huán)境參數(shù)與海底地聲參數(shù)共同構(gòu)成了水下聲傳播計(jì)算的環(huán)境參數(shù)，而地聲參數(shù)反演技術(shù)一直是水聲學(xué)中的熱門(mén)課題，然而現(xiàn)有的參數(shù)反演算法一般只給出參數(shù)的點(diǎn)估計(jì)，并且很少對(duì)參數(shù)的不確定性進(jìn)行進(jìn)一步描述。相比較而言貝葉斯推理能通過(guò)后驗(yàn)概率密度給出整個(gè)模型參數(shù)的概率分布，所以采用優(yōu)化算法容易丟失一些地聲參數(shù)反演問(wèn)題解的重要信息，而對(duì)于具有多參數(shù)的模型，貝葉斯方法能用來(lái)對(duì)任何特定參數(shù)或參數(shù)的函數(shù)進(jìn)行推理分析。有必要采用貝葉斯理論對(duì)水聲傳播模型中地聲參數(shù)估計(jì)與不確定性進(jìn)行研究，獲取地聲參數(shù)的后驗(yàn)概率密度，在此基礎(chǔ)上研究地聲參數(shù)不確定性對(duì)傳播損失的影響。

2　模型與算法

2.1水聲環(huán)境不確定性推理模型

將水聲環(huán)境觀測(cè)數(shù)據(jù)、地聲參數(shù)以及傳播損失劃分為三個(gè)域，如圖1［1-2］所示，它反映了從數(shù)據(jù)的觀測(cè)域d∈D到環(huán)境域m∈M，再到應(yīng)用域u∈U的映射關(guān)系。這也說(shuō)明了從海洋水聲觀測(cè)（觀測(cè)域）到傳播損失（應(yīng)用域）估計(jì)的總體思路。地聲參數(shù)的估計(jì)作為中間步驟給出地聲參數(shù)的后驗(yàn)概率分布p（m|d）（環(huán)境域）。在此基礎(chǔ)上，傳播損失的概率分布p（u|d）可由Monte Carlo求積獲得。由圖中可以看出，觀測(cè)數(shù)據(jù)d和應(yīng)用域u都與m相關(guān)，并且由映射關(guān)系D（m）和U（m）相聯(lián)系，如果映射關(guān)系是唯一的，那么u=U（D-1（d））［3］。

圖1　觀測(cè)域d∈D到環(huán)境域m∈M到應(yīng)用域u∈U的映射關(guān)系圖［1-2］Fig.1　An observation d　∈　D is mapped into a distribution of environmental parameters m∈M that potentially could have generated it.　These environmental parameters are then mapped into the usage domain u∈U［1-2］

在貝葉斯統(tǒng)計(jì)方法中，認(rèn)為總體分布中未知的環(huán)境參數(shù)向量m為隨機(jī)向量，在對(duì)m進(jìn)行統(tǒng)計(jì)推斷時(shí)，需要對(duì)m設(shè)定一個(gè)先驗(yàn)分布p（d|m），它是在獲得觀測(cè)數(shù)據(jù)前對(duì)于未知環(huán)境參數(shù)信息的概率表述。獲得觀測(cè)數(shù)據(jù)d之后，使用貝葉斯公式就可將參數(shù)的先驗(yàn)分布更新為后驗(yàn)分布。根據(jù)Bayesian理論，地聲參數(shù)的后驗(yàn)分布、觀測(cè)信息、先驗(yàn)分布之間有如下關(guān)系［1-2］：

式中：p（d）為歸一化因子，由于觀測(cè)數(shù)據(jù)已經(jīng)給出，它是一個(gè)與m無(wú)關(guān)的常數(shù)，可忽略；其中條件概率密度p（d|m）又稱似然函數(shù)，也可寫(xiě)作L（m），表示估計(jì)的環(huán)境參數(shù)和實(shí)際的觀測(cè)數(shù)據(jù)的擬合程度，值越大表明擬合程度越高，反之越低。在此基礎(chǔ)上將先驗(yàn)信息和觀測(cè)信息相結(jié)合，就可以得到反映未知的隨機(jī)環(huán)境參數(shù)向量m綜合信息的后驗(yàn)概率密度p（m|d），它定義在整個(gè)解空間，表示了問(wèn)題的“完全”解，而構(gòu)造并計(jì)算p（m|d）是貝葉斯推論的核心問(wèn)題。

2.2MCMC采樣方法

一旦獲得了參數(shù)的后驗(yàn)分布，就可以獲得參數(shù)的任何特征，如單個(gè)參數(shù)的邊緣分布或均值、方差等［4］。然而除非對(duì)于很簡(jiǎn)單的情況，后驗(yàn)概率分布都不能以解析解的形式給出，所以必須采用數(shù)值積分方法。馬爾科夫鏈蒙特卡羅（MCMC）就是一種為了獲得參數(shù)后驗(yàn)分布而發(fā)展起來(lái)，對(duì)復(fù)雜問(wèn)題在高維空間上的數(shù)值積分方案。MCMC方法的基本原理就是基于建立的平穩(wěn)分布為f（x）的Markov鏈來(lái)獲得f（x）的樣本。產(chǎn)生若干條獨(dú)立并行的Markov鏈來(lái)探索模型參數(shù)空間，通過(guò)不斷更新樣本信息而使Markov鏈?zhǔn)諗坑诟吒怕拭芏葏^(qū)，也就是Bayesian方法中的最大后驗(yàn)估計(jì)。

根據(jù)Markov鏈所用轉(zhuǎn)移概率的不同，MCMC方法具有多種抽樣算法，如：Gibbs抽樣、Metropolis抽樣、Metropolis-Hastings抽樣［5-7］。本文采用目前較為常用的Metropolis-Hastings算法，其算法流程如下：

（1）t=1，將初值賦給不同變量；

（2）while t＜T

（a）t=t+1；

（b）產(chǎn)生服從建議分布的推薦變量值m?；

（c）計(jì)算接受概率

（d）產(chǎn)生服從均勻分布的隨機(jī)數(shù)u～U（0，1）；

（e）若u＜α，則mt=m?，否則mt=mt-1；

（3）重復(fù)第（2）步（a）～（e），直到產(chǎn)生T個(gè)樣本為止。

3　似然函數(shù)與后驗(yàn)概率密度

3.1似然函數(shù)

若觀測(cè)的數(shù)據(jù)矢量表示為d=D（m）+e，式中D（m）=d（m）s，其中s表示聲源信息，且e是服從均值為0，協(xié)方差為Ce正態(tài)分布的誤差項(xiàng)，則似然函數(shù)為［2］

式中：N為數(shù)據(jù)點(diǎn)的數(shù)量，上標(biāo)?表示共軛轉(zhuǎn)置，轉(zhuǎn)換函數(shù)d（m）由水聲傳播模型獲得。用獨(dú)立同分布的誤差Ce=vI來(lái)描述數(shù)據(jù)的不確定性，其中v為誤差方差，那么當(dāng)?lgL/?s=0時(shí)

將式（3）代入公式（2），似然函數(shù)則變?yōu)?/p>

式中

為目標(biāo)函數(shù)。在此，將誤差方差視為冗余參量并通過(guò)對(duì)公式（4）求積消除v，那么

式中p（v）=1/v［2］。至此，似然函數(shù)可寫(xiě)為［2］

3.2地聲參數(shù)與傳播損失的后驗(yàn)概率密度

獲得地聲參數(shù)的似然函數(shù)L（m），結(jié)合參數(shù)相應(yīng)的先驗(yàn)信息，便可求得環(huán)境參數(shù)的后驗(yàn)概率密度p（m|d）。而對(duì)于利用水聲環(huán)境來(lái)說(shuō)，傳播損失（TL）與聲納探測(cè)性能估計(jì)直接相關(guān)，在模型中用u表示，它可表示為I個(gè)離散位置的傳播損失，即ui=u（ri，zi），i=1，···，I，其中ri和zi分別表示距離和深度上的離散值。那么傳播損失的后驗(yàn)概率密度p（u|d）可由u和m的聯(lián)合概率密度函數(shù)獲得［2］

假設(shè)所有的不確定性均在數(shù)據(jù)域d中，并且數(shù)據(jù)域d中所有的信息都被映射到水聲環(huán)境域m，那么

條件概率密度函數(shù)p（u|m）用于描述環(huán)境知識(shí)的不完全情況下，前向映射的不確定性。假設(shè)水聲傳播模型是準(zhǔn)確的，函數(shù)關(guān)系表示為ui=Ui（m），i=1，···，I，簡(jiǎn)記為u=U（m），這樣對(duì)于每一個(gè)環(huán)境值m，都能給出一個(gè)較為準(zhǔn)確的傳播損失u。此時(shí)概率密度可寫(xiě)為

如圖2所示，直角部分可以采用兩條Clothoid曲線A0和A1過(guò)渡，曲線A0和A1分別起始于點(diǎn)P0與P1，相遇于P，在起始點(diǎn)處與直線曲率相同為0，在點(diǎn)P兩條曲線曲率相同并且切線平行反向，這樣保證了整條運(yùn)動(dòng)軌跡曲率G2連續(xù)。

MCMC方法對(duì)于Bayesian推論問(wèn)題較為適用，公式（12）中的積分即δ（U（m）-u）的期望函數(shù)，可采用蒙特卡洛積分方法求解，該方法關(guān)鍵在于獲得服從復(fù)雜目標(biāo)分布的大量樣本，而MCMC方法可對(duì)模型參數(shù)后驗(yàn)分布p（m|d）進(jìn)行Metropolis-Hastings抽樣{m（t）}來(lái)近似。

至此，我們便得到了傳播損失，即應(yīng)用域u的分布，即傳播損失的概率分布，它包含了水聲環(huán)境參數(shù)估計(jì)的不確定性。為了更好的描述傳播損失的不確定性情況，在此采用后驗(yàn)預(yù)測(cè)概率分布50%處的中值u50%和β1%處與β2%處之間的區(qū)域［uβ1%，uβ2%］來(lái)描述傳播損失分布。通過(guò)滿足下式進(jìn)行求解：

根據(jù)地聲環(huán)境參數(shù)不確定性程度及其傳播損失后驗(yàn)預(yù)測(cè)概率分布公式，本文β1、β2分別取10與90，并以β2%與β1%之間的范圍表示傳播損失估計(jì)范圍，把其稱為傳播損失（β2-β1）%的可信區(qū)間（Credibility interval）。

4　實(shí)驗(yàn)與算例

4.1實(shí)驗(yàn)概況

2012年6月30日在黃海北部某海域進(jìn)行了水聲傳播實(shí)驗(yàn)，實(shí)驗(yàn)過(guò)程及水聽(tīng)器布放如圖2所示，接收船靜止布放16元垂直陣，水聽(tīng)器之間距離2 m，發(fā)射船以10節(jié)速度遠(yuǎn)離接收船，每隔1 min投放一枚手榴彈作為爆炸聲源，接收船接收爆炸聲信號(hào)，實(shí)驗(yàn)海區(qū)聲速剖面如圖3所示。實(shí)驗(yàn)數(shù)據(jù)獲取與處理參見(jiàn)文獻(xiàn)［8］。

圖3　實(shí)驗(yàn)海區(qū)聲速剖面Fig.3 Sound speed profile in experimental area

海洋環(huán)境模型是進(jìn)行地聲環(huán)境參數(shù)及其不確定性估計(jì)的基礎(chǔ)。目前，國(guó)內(nèi)外估計(jì)淺海水聲環(huán)境參數(shù)時(shí)大都采用分層的海底模型，其中兩層海底模型應(yīng)用最多。根據(jù)淺海水聲環(huán)境參數(shù)估計(jì)的特點(diǎn)及試驗(yàn)所在海區(qū)的底質(zhì)情況，采用兩層海底模型，如圖2所示?？紤]待估的地聲環(huán)境參數(shù)為7個(gè)，分別為：沉積層密度rhos（g/cm3）、巖石層密度rhor（g/cm3）、沉積層聲速Csed（m/s）、巖石層聲速Crock（m/s）、沉積層厚度Dsed（m）、巖石層厚度Drock（m）、沉積層衰減alphy（dB/λ）。推理模型中的映射關(guān)系D（m）和U（m）采用KRAKEN簡(jiǎn)正波傳播模型。各參數(shù)的先驗(yàn)分布，根據(jù)參數(shù)的取值范圍設(shè)定為截?cái)嗾龖B(tài)分布，沉積層密度、聲速的均值根據(jù)海圖標(biāo)注的底質(zhì)給出經(jīng)驗(yàn)取值。各參數(shù)初始取值分別為：1.83 g/cm3、2.5 g/cm3、1677 m/s、1850 m/s、4 m、15 m、0.1 dB/λ。

圖2　實(shí)驗(yàn)示意圖及淺海環(huán)境模型Fig.2　Schematic diagram of experiment and shallow water model

4.3結(jié)果分析

仿真計(jì)算時(shí)采用聲源頻率100 Hz，各參數(shù)的二維后驗(yàn)概率密度分布如圖4所示，其中顏色越亮代表后驗(yàn)概率密度值越高，它表達(dá)了各參數(shù)在二維空間中的取值范圍。同時(shí)，圖5給出了各參數(shù)一維的后驗(yàn)概率密度分布圖，表1給出了各地聲參數(shù)后驗(yàn)概率密度的均值與均方差。從結(jié)果可以看出，采用MCMC方法獲得的地聲參數(shù)后驗(yàn)概率密度可以對(duì)先驗(yàn)知識(shí)進(jìn)行修正，而其中因沉積層衰減對(duì)模型計(jì)算的敏感程度較高，后驗(yàn)估計(jì)對(duì)先驗(yàn)信息有較大的改善。從圖中也可以看出，與以往地聲參數(shù)反演問(wèn)題的點(diǎn)估計(jì)相比，本文給出了參數(shù)的概率密度分布，包含的參數(shù)信息也更全面。

根據(jù)水聲環(huán)境不確定性推理模型，得出地聲參數(shù)的后驗(yàn)概率密度可作為傳播損失不確定性估計(jì)的中間步驟，由公式（13）便可求得傳播損失的后驗(yàn)概率密度，并根據(jù)公式（14）求得某一距離上滿足一定概率的傳播損失區(qū)間。

如圖6、7所示，分別是接收深度13 m和31 m時(shí)傳播損失的不確定性估計(jì)，其中（a）、（b）為10 km以內(nèi)，傳播損失的均值（中間紅色實(shí)線）與80%的可信區(qū)間（兩條黑色虛線之間的區(qū)間）；（b）和（c）分別是3 km和8 km處傳播損失的概率密度直方圖。從圖中可以看出以上兩個(gè)深度傳播損失的實(shí)驗(yàn)數(shù)據(jù)（“*”表示）均落于80%的可信區(qū)間，而且從傳播趨勢(shì)以及均方差可以看出，距離聲源越遠(yuǎn)，不確定性越大。

圖4　地聲環(huán)境參數(shù)二維后驗(yàn)概率密度Fig.4 2D posterior probability densities of geoacoustic parameters

圖5　地聲環(huán)境參數(shù)一維邊緣后驗(yàn)概率密度Fig.5 1D marginal posterior probability densities of geoacoustic parameters

圖6　接收深度13 m時(shí)的傳播損失不確定性Fig.6 Uncertainty of transmission loss（TL）at 13 m depth

圖7　接收深度31 m時(shí)的傳播損失不確定性Fig.7 Uncertainty of transmission loss（TL）at 31 m depth

表1　地聲參數(shù)后驗(yàn)概率密度均值與均方差Table 1　Mean value and mean square deviation of geoacoustic parameters posterior probability densities mean

同理，其他深度上傳播損失的均值與80%的可信區(qū)間如圖8所示，從圖中可以看出除表層（7 m、9 m）和躍層中下深度（23 m、25 m、27 m）的實(shí)驗(yàn)數(shù)據(jù)與計(jì)算數(shù)據(jù)吻合程度相對(duì)較差，其余深度的實(shí)驗(yàn)數(shù)據(jù)均落在傳播損失80%的可信區(qū)間之內(nèi)。這是因?yàn)榭拷Ｃ娴乃?tīng)器受海面影響較大，數(shù)據(jù)的質(zhì)量相對(duì)較低，而23～27 m正好是躍層向均勻?qū)舆^(guò)渡的深度，受躍變層不穩(wěn)定性和海底因素等影響較大，因此計(jì)算的傳播損失區(qū)間與實(shí)驗(yàn)數(shù)據(jù)有一定出入。

圖8　不同接收深度時(shí)傳播損失的80%可信區(qū)間與實(shí)驗(yàn)數(shù)據(jù)對(duì)比Fig.8 Predicted and measured TL and its 80%credibility interval at different depth

5　結(jié)論

本文首先通過(guò)貝葉斯理論建立水聲環(huán)境不確定性推理模型，介紹了MCMC算法，推導(dǎo)了地聲參數(shù)的似然函數(shù)以及地聲參數(shù)和傳播損失的后驗(yàn)概率密度，采用2012年黃海北部某海域水聲傳播實(shí)驗(yàn)數(shù)據(jù)反演了地聲參數(shù)的后驗(yàn)概率密度，并在此基礎(chǔ)上對(duì)傳播損失的不確定性進(jìn)行了估計(jì)，得到了傳播損失80%的可信區(qū)間。經(jīng)與實(shí)驗(yàn)數(shù)據(jù)對(duì)比，驗(yàn)證了該方法對(duì)地聲參數(shù)的不確定性估計(jì)具有一定的優(yōu)勢(shì)，與已有的地聲參數(shù)反演與水聲不確定性研究相比，本文的主要改進(jìn)之處在于將先驗(yàn)信息和觀測(cè)信息相結(jié)合，得到反映地聲參數(shù)向量綜合信息的后驗(yàn)概率密度，已有的地聲參數(shù)反演是在選定的變化范圍根據(jù)觀測(cè)數(shù)據(jù)進(jìn)行尋優(yōu)計(jì)算使代價(jià)函數(shù)達(dá)到最小值，來(lái)確定地聲參數(shù)，這里計(jì)及環(huán)境參數(shù)的隨機(jī)因素，應(yīng)用概率統(tǒng)計(jì)方法獲取地聲參數(shù)后驗(yàn)概率密度均值表征反演后的取值概率，比較科學(xué)。并應(yīng)用到與環(huán)境參數(shù)密切相關(guān)的傳播損失的不確定性的估計(jì)，表示了問(wèn)題的“完全”解。但是該方法計(jì)算效率相對(duì)較低，在抽樣算法方面還可進(jìn)行進(jìn)一步的深入研究。

［1］GERSTOFT P，HUANG C F，HODGKISS W S.Estimation of transmission loss in the presence of geoacoustic inversion uncertainty［J］.IEEE Journal of Ocean Engineering，2006，31（2）：299-307.

［2］HUANG C F，GERSTOFT P，HODGKISS W S.Validation of statistical estimation of transmission loss in the presence of geoacoustic inversion uncertainty［J］.Journal of Acoustical Society of America，2006，120（4）：1932-1941.

［3］笪良龍.海洋水聲環(huán)境效應(yīng)建模與應(yīng)用［M］.北京：科學(xué)出版社，2012：131-134. DA Lianglong.Ocean acoustic environment modeling and application［M］.Beijing：Science Press，2012：131-134.

［4］朱嵩，毛根海，程偉平，等.基于貝葉斯推理的水環(huán)境系統(tǒng)參數(shù)識(shí)別［J］.江蘇大學(xué)學(xué)報(bào)（自然科學(xué)版），2007，28（3）：237-240. ZHU Song，MAO Genhai，CHENG Weiping，et al.Parameters identification for water environmental system based on Bayesian inference［J］.Journal of Jiangsu University（Natural Science Edition），2007，28（3）：237-240.

［5］PENG R H，CHEN R R，F(xiàn)ARHANG-BOROUJENY B. Markov Chain Monte Carlo detectors for channels with intersymbol interference［J］.IEEE Transactions on Signal Processing，2010，58（4）：2206-2217.

［6］SPALL J C.Estimation via Markov Chain Monte Carlo［J］. IEEE Control Systems Magazine，2003，23（2）：34-45.

［7］JOSEPH L S，RECAI M Y.Computational strategies for multivariate linear mixed-effects models with missing values［J］.Journal of Computational and Graphical Statistics，2002，11（2）：437-457.

［8］過(guò)武宏，笪良龍，趙建昕.動(dòng)態(tài)水聲環(huán)境不確定性的估計(jì)與分析［J］.應(yīng)用聲學(xué)，2013，32（6）：464-472. GUO Wuhong，DA Lianglong，ZHAO Jianxin.Uncertainty estimation and analysis for dynamic underwater acoustic environment［J］.Applied Acoustics，2013，32（6）：464-472.

Estimation and modeling of geoacoustic parameters and transmission loss uncertainty?

GUO Wuhong?DA LianglongZHAO Jianxin
（Navy Submarine Academy，Qingdao 266071，China）

Uncertainty of geoacoustic parameters has a great effect on acoustic propagation.In this paper,the uncertainty inference model of the acoustic environment is established via Bayes theorem.Then the likelihood function of the geoacoustic parameters and the posterior probability distribution of the geoacoustic parameters and transmission loss are deduced.Furthermore,a simulation is made using Markov Chain Monte Carlo（MCMC）.Two-dimensional and one-dimensional posterior probability densities of the geoacoustic parameters are given.Meanwhile,the uncertainty of transmission loss is estimated.The mean and 80%credibility interval of transmission loss in different depths are given.The results of the simulation and experiments show that the method is good at geoacoustic parameters inversion and uncertainty estimation,and it can get the uncertainty estimation of the transmission loss that is caused by uncertainty of the geoacoustic parameters.

Bayesian，Posterior probability density，Geoacoustic parameters，Markov Chain Monte Carlo

TP391

A

1000-310X（2015）01-0071-08

10.11684/j.issn.1000-310X.2015.01.011

2014-01-18收稿；2014-05-10定稿

?中國(guó)博士后科學(xué)基金項(xiàng)目（20110491884），總裝預(yù)研基金項(xiàng)目（9140A03060213JB15039）

過(guò)武宏（1980-），男，浙江杭州人，博士后，研究方向：海洋、水聲環(huán)境預(yù)報(bào)與不確定性。

E-mail：g1w2h31980@163.com