地声参数及传播损失不确定性估计与建模∗

过武宏　笪良龙　赵建昕

（海军潜艇学院　青岛　266071）

地声参数及传播损失不确定性估计与建模∗

过武宏†笪良龙赵建昕

（海军潜艇学院青岛266071）

地声参数的不确定性对水声传播具有重要的影响。通过贝叶斯理论建立水声环境不确定性推理模型，理论推导了地声参数的似然函数以及地声参数和传播损失的后验概率密度，并采用MCMC（Markov Chain Monte Carlo）进行了仿真计算，给出了地声参数的二维后验联合概率密度和一维边缘概率密度，在此基础上对传播损失的不确定性进行了估计，得到了传播损失80%的可信区间。仿真和实验结果表明，该方法适用于地声参数反演和不确定性估计，并能获取因地声参数不确定性导致的传播损失不确定性估计。

贝叶斯，后验概率密度，地声参数，马尔可夫链蒙特卡洛

1　引言

海洋水体环境参数与海底地声参数共同构成了水下声传播计算的环境参数，而地声参数反演技术一直是水声学中的热门课题，然而现有的参数反演算法一般只给出参数的点估计，并且很少对参数的不确定性进行进一步描述。相比较而言贝叶斯推理能通过后验概率密度给出整个模型参数的概率分布，所以采用优化算法容易丢失一些地声参数反演问题解的重要信息，而对于具有多参数的模型，贝叶斯方法能用来对任何特定参数或参数的函数进行推理分析。有必要采用贝叶斯理论对水声传播模型中地声参数估计与不确定性进行研究，获取地声参数的后验概率密度，在此基础上研究地声参数不确定性对传播损失的影响。

2　模型与算法

2.1水声环境不确定性推理模型

将水声环境观测数据、地声参数以及传播损失划分为三个域，如图1［1-2］所示，它反映了从数据的观测域d∈D到环境域m∈M，再到应用域u∈U的映射关系。这也说明了从海洋水声观测（观测域）到传播损失（应用域）估计的总体思路。地声参数的估计作为中间步骤给出地声参数的后验概率分布p（m|d）（环境域）。在此基础上，传播损失的概率分布p（u|d）可由Monte Carlo求积获得。由图中可以看出，观测数据d和应用域u都与m相关，并且由映射关系D（m）和U（m）相联系，如果映射关系是唯一的，那么u=U（D-1（d））［3］。

图1　观测域d∈D到环境域m∈M到应用域u∈U的映射关系图［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］

在贝叶斯统计方法中，认为总体分布中未知的环境参数向量m为随机向量，在对m进行统计推断时，需要对m设定一个先验分布p（d|m），它是在获得观测数据前对于未知环境参数信息的概率表述。获得观测数据d之后，使用贝叶斯公式就可将参数的先验分布更新为后验分布。根据Bayesian理论，地声参数的后验分布、观测信息、先验分布之间有如下关系［1-2］：

式中：p（d）为归一化因子，由于观测数据已经给出，它是一个与m无关的常数，可忽略；其中条件概率密度p（d|m）又称似然函数，也可写作L（m），表示估计的环境参数和实际的观测数据的拟合程度，值越大表明拟合程度越高，反之越低。在此基础上将先验信息和观测信息相结合，就可以得到反映未知的随机环境参数向量m综合信息的后验概率密度p（m|d），它定义在整个解空间，表示了问题的“完全”解，而构造并计算p（m|d）是贝叶斯推论的核心问题。

2.2MCMC采样方法

一旦获得了参数的后验分布，就可以获得参数的任何特征，如单个参数的边缘分布或均值、方差等［4］。然而除非对于很简单的情况，后验概率分布都不能以解析解的形式给出，所以必须采用数值积分方法。马尔科夫链蒙特卡罗（MCMC）就是一种为了获得参数后验分布而发展起来，对复杂问题在高维空间上的数值积分方案。MCMC方法的基本原理就是基于建立的平稳分布为f（x）的Markov链来获得f（x）的样本。产生若干条独立并行的Markov链来探索模型参数空间，通过不断更新样本信息而使Markov链收敛于高概率密度区，也就是Bayesian方法中的最大后验估计。

根据Markov链所用转移概率的不同，MCMC方法具有多种抽样算法，如：Gibbs抽样、Metropolis抽样、Metropolis-Hastings抽样［5-7］。本文采用目前较为常用的Metropolis-Hastings算法，其算法流程如下：

（1）t=1，将初值赋给不同变量；

（2）while t＜T

（a）t=t+1；

（b）产生服从建议分布的推荐变量值m∗；

（c）计算接受概率

（d）产生服从均匀分布的随机数u～U（0，1）；

（e）若u＜α，则mt=m∗，否则mt=mt-1；

（3）重复第（2）步（a）～（e），直到产生T个样本为止。

3　似然函数与后验概率密度

3.1似然函数

若观测的数据矢量表示为d=D（m）+e，式中D（m）=d（m）s，其中s表示声源信息，且e是服从均值为0，协方差为Ce正态分布的误差项，则似然函数为［2］

式中：N为数据点的数量，上标∗表示共轭转置，转换函数d（m）由水声传播模型获得。用独立同分布的误差Ce=vI来描述数据的不确定性，其中v为误差方差，那么当∂lgL/∂s=0时

将式（3）代入公式（2），似然函数则变为

式中

为目标函数。在此，将误差方差视为冗余参量并通过对公式（4）求积消除v，那么

式中p（v）=1/v［2］。至此，似然函数可写为［2］

3.2地声参数与传播损失的后验概率密度

获得地声参数的似然函数L（m），结合参数相应的先验信息，便可求得环境参数的后验概率密度p（m|d）。而对于利用水声环境来说，传播损失（TL）与声纳探测性能估计直接相关，在模型中用u表示，它可表示为I个离散位置的传播损失，即ui=u（ri，zi），i=1，···，I，其中ri和zi分别表示距离和深度上的离散值。那么传播损失的后验概率密度p（u|d）可由u和m的联合概率密度函数获得［2］

假设所有的不确定性均在数据域d中，并且数据域d中所有的信息都被映射到水声环境域m，那么

条件概率密度函数p（u|m）用于描述环境知识的不完全情况下，前向映射的不确定性。假设水声传播模型是准确的，函数关系表示为ui=Ui（m），i=1，···，I，简记为u=U（m），这样对于每一个环境值m，都能给出一个较为准确的传播损失u。此时概率密度可写为

如图2所示，直角部分可以采用两条Clothoid曲线A0和A1过渡，曲线A0和A1分别起始于点P0与P1，相遇于P，在起始点处与直线曲率相同为0，在点P两条曲线曲率相同并且切线平行反向，这样保证了整条运动轨迹曲率G2连续。

MCMC方法对于Bayesian推论问题较为适用，公式（12）中的积分即δ（U（m）-u）的期望函数，可采用蒙特卡洛积分方法求解，该方法关键在于获得服从复杂目标分布的大量样本，而MCMC方法可对模型参数后验分布p（m|d）进行Metropolis-Hastings抽样{m（t）}来近似。

至此，我们便得到了传播损失，即应用域u的分布，即传播损失的概率分布，它包含了水声环境参数估计的不确定性。为了更好的描述传播损失的不确定性情况，在此采用后验预测概率分布50%处的中值u50%和β1%处与β2%处之间的区域［uβ1%，uβ2%］来描述传播损失分布。通过满足下式进行求解：

根据地声环境参数不确定性程度及其传播损失后验预测概率分布公式，本文β1、β2分别取10与90，并以β2%与β1%之间的范围表示传播损失估计范围，把其称为传播损失（β2-β1）%的可信区间（Credibility interval）。

4　实验与算例

4.1实验概况

2012年6月30日在黄海北部某海域进行了水声传播实验，实验过程及水听器布放如图2所示，接收船静止布放16元垂直阵，水听器之间距离2 m，发射船以10节速度远离接收船，每隔1 min投放一枚手榴弹作为爆炸声源，接收船接收爆炸声信号，实验海区声速剖面如图3所示。实验数据获取与处理参见文献［8］。

图3　实验海区声速剖面Fig.3 Sound speed profile in experimental area

海洋环境模型是进行地声环境参数及其不确定性估计的基础。目前，国内外估计浅海水声环境参数时大都采用分层的海底模型，其中两层海底模型应用最多。根据浅海水声环境参数估计的特点及试验所在海区的底质情况，采用两层海底模型，如图2所示。考虑待估的地声环境参数为7个，分别为：沉积层密度rhos（g/cm3）、岩石层密度rhor（g/cm3）、沉积层声速Csed（m/s）、岩石层声速Crock（m/s）、沉积层厚度Dsed（m）、岩石层厚度Drock（m）、沉积层衰减alphy（dB/λ）。推理模型中的映射关系D（m）和U（m）采用KRAKEN简正波传播模型。各参数的先验分布，根据参数的取值范围设定为截断正态分布，沉积层密度、声速的均值根据海图标注的底质给出经验取值。各参数初始取值分别为：1.83 g/cm3、2.5 g/cm3、1677 m/s、1850 m/s、4 m、15 m、0.1 dB/λ。

图2　实验示意图及浅海环境模型Fig.2　Schematic diagram of experiment and shallow water model

4.3结果分析

仿真计算时采用声源频率100 Hz，各参数的二维后验概率密度分布如图4所示，其中颜色越亮代表后验概率密度值越高，它表达了各参数在二维空间中的取值范围。同时，图5给出了各参数一维的后验概率密度分布图，表1给出了各地声参数后验概率密度的均值与均方差。从结果可以看出，采用MCMC方法获得的地声参数后验概率密度可以对先验知识进行修正，而其中因沉积层衰减对模型计算的敏感程度较高，后验估计对先验信息有较大的改善。从图中也可以看出，与以往地声参数反演问题的点估计相比，本文给出了参数的概率密度分布，包含的参数信息也更全面。

根据水声环境不确定性推理模型，得出地声参数的后验概率密度可作为传播损失不确定性估计的中间步骤，由公式（13）便可求得传播损失的后验概率密度，并根据公式（14）求得某一距离上满足一定概率的传播损失区间。

如图6、7所示，分别是接收深度13 m和31 m时传播损失的不确定性估计，其中（a）、（b）为10 km以内，传播损失的均值（中间红色实线）与80%的可信区间（两条黑色虚线之间的区间）；（b）和（c）分别是3 km和8 km处传播损失的概率密度直方图。从图中可以看出以上两个深度传播损失的实验数据（“*”表示）均落于80%的可信区间，而且从传播趋势以及均方差可以看出，距离声源越远，不确定性越大。

图4　地声环境参数二维后验概率密度Fig.4 2D posterior probability densities of geoacoustic parameters

图5　地声环境参数一维边缘后验概率密度Fig.5 1D marginal posterior probability densities of geoacoustic parameters

图6　接收深度13 m时的传播损失不确定性Fig.6 Uncertainty of transmission loss（TL）at 13 m depth

图7　接收深度31 m时的传播损失不确定性Fig.7 Uncertainty of transmission loss（TL）at 31 m depth

表1　地声参数后验概率密度均值与均方差Table 1　Mean value and mean square deviation of geoacoustic parameters posterior probability densities mean

同理，其他深度上传播损失的均值与80%的可信区间如图8所示，从图中可以看出除表层（7 m、9 m）和跃层中下深度（23 m、25 m、27 m）的实验数据与计算数据吻合程度相对较差，其余深度的实验数据均落在传播损失80%的可信区间之内。这是因为靠近海面的水听器受海面影响较大，数据的质量相对较低，而23～27 m正好是跃层向均匀层过渡的深度，受跃变层不稳定性和海底因素等影响较大，因此计算的传播损失区间与实验数据有一定出入。

图8　不同接收深度时传播损失的80%可信区间与实验数据对比Fig.8 Predicted and measured TL and its 80%credibility interval at different depth

5　结论

本文首先通过贝叶斯理论建立水声环境不确定性推理模型，介绍了MCMC算法，推导了地声参数的似然函数以及地声参数和传播损失的后验概率密度，采用2012年黄海北部某海域水声传播实验数据反演了地声参数的后验概率密度，并在此基础上对传播损失的不确定性进行了估计，得到了传播损失80%的可信区间。经与实验数据对比，验证了该方法对地声参数的不确定性估计具有一定的优势，与已有的地声参数反演与水声不确定性研究相比，本文的主要改进之处在于将先验信息和观测信息相结合，得到反映地声参数向量综合信息的后验概率密度，已有的地声参数反演是在选定的变化范围根据观测数据进行寻优计算使代价函数达到最小值，来确定地声参数，这里计及环境参数的随机因素，应用概率统计方法获取地声参数后验概率密度均值表征反演后的取值概率，比较科学。并应用到与环境参数密切相关的传播损失的不确定性的估计，表示了问题的“完全”解。但是该方法计算效率相对较低，在抽样算法方面还可进行进一步的深入研究。

［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］笪良龙.海洋水声环境效应建模与应用［M］.北京：科学出版社，2012：131-134. DA Lianglong.Ocean acoustic environment modeling and application［M］.Beijing：Science Press，2012：131-134.

［4］朱嵩，毛根海，程伟平，等.基于贝叶斯推理的水环境系统参数识别［J］.江苏大学学报（自然科学版），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，FARHANG-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］过武宏，笪良龙，赵建昕.动态水声环境不确定性的估计与分析［J］.应用声学，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定稿

∗中国博士后科学基金项目（20110491884），总装预研基金项目（9140A03060213JB15039）

过武宏（1980-），男，浙江杭州人，博士后，研究方向：海洋、水声环境预报与不确定性。

E-mail：g1w2h31980@163.com