Bayes匹配场地声参数反演:多步退火Gibbs采样算法

高飞， 潘长明， 孙磊,2

(1.海军海洋测绘研究所, 天津 300061; 2.哈尔滨工程大学 水声工程学院, 黑龙江 哈尔滨 150001)

Bayes匹配场地声参数反演:多步退火Gibbs采样算法

高飞1， 潘长明1， 孙磊1,2

(1.海军海洋测绘研究所, 天津 300061; 2.哈尔滨工程大学 水声工程学院, 黑龙江 哈尔滨 150001)

为实现海洋地声环境参数的快速高效反演，提出多步退火Gibbs(Multi-AG)采样算法，可消除因参数搜索空间设置对参数反演结果的影响，并有效解决Bayes匹配场高维参数反演过程中常见的运算量大、旁瓣高等问题。分析待反演地声参数对匹配场处理器的敏感性，用以制定多步反演与退火方案，利用Gibbs采样算法反演敏感性级别最高的参数，计算其均值并代入后续反演步骤，进而采用退火Gibbs采样算法逐步反演后续参数；利用数值仿真实验对比Metropolis-Hastings算法、Gibbs采样算法、快速Gibbs采样算法和Multi-AG采样算法的反演效果。实验结果表明，与其他3种算法相比，Multi-AG采样算法可通过最小的计算量得到均方差最小、精度最高的参数反演结果。

声学； 多步退火Gibbs采样算法； 地声参数反演； Gibbs采样； Bayes匹配场

0 引言

海底地声参数信息是海洋水声环境的重要参数组成之一，在海洋资源勘探、水下工程及声纳系统探测领域应用广泛，极具军事价值和民用价值。当前水声调查通常包括海洋环境噪声、混响、传播损失、声速剖面、水深及海面气象等要素，且相关理论研究及调测设备发展相对成熟[1-2]。然而，受作业难度与技术要求制约，当前浅地层地声参数实测手段仍无法满足大规模工程作业的要求，通过间接反演方法获取地声参数信息仍是当前及未来较长时间内的主要方式。

海洋环境、声源、垂线阵(VLA)声信号是Bayes匹配场3要素[3]，由第1、第3项得到第2项为声源定位，由后两项得到第1项称为海洋环境参数反演。当前高维地声参数Bayes匹配场反演技术主要面临计算量大、旁瓣高等问题，也一直是研究改进的焦点之一。

Bayes匹配场反演主要可从匹配处理器、随机数据统计算法等方面进行优化。常见的Bartlett线性处理器[4]最早用于度量线列阵测量声场与拷贝声场的匹配程度，具有较好的稳健性；自适应场处理器[5]与多约束场处理器[6]通过权值矢量/后验概率密度(PDF)约束，有效融合了阵列信号处理的优点，可抑制旁瓣与海洋环境不确定性，提供主瓣保护；Green函数处理器较好地处理了海洋环境不确定性，可削弱环境参数失配对匹配场反演精度的影响[7]。马尔可夫链蒙特卡洛(MCMC)随机统计思想[8]可反演得到参数PDF，但该方法在处理高维数据时计算量过大；遗传算法[9-10]和最优算法[11]在参数搜索空间中寻优，可得到高精度反演结果；Gibbs采样算法[12]是高维MCMC反演算法的一种改进，可极大地减少运算量；聚焦法[13]常用于声源位置与参数共同反演，首先采用重要度采样获取局部最优解，进而通过Metropolis Gibbs采样算法确定出声源位置与参数PDF。

综上所述可知，已有的研究成果未考虑各海洋环境参数对匹配场处理器的敏感性差异，导致敏感性较低的参数反演精度不高，且受搜索空间限制影响显著。本文提出多步退火Gibbs(Multi-AG)采样算法，以有效解决上述问题，进而得到快速高精度反演结果。

1 Bayes匹配场反演理论

以m表示地声参数信息(如沉积层密度、声速、厚度等)，s为复数形式声压信号，则依据Bayes计算法可得

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

(1)

式中：P(m)为地声参数m的先验概率密度，通常假设其在各自取值区间内服从均匀分布，重要度采样时与具体的密度分布函数g(m)相关；P(s|m)为已知m时s的条件概率密度，当水声环境参数m已知时，声信号s可利用水声数值模型计算得到；P(m|s)为PDF函数。

Bayes匹配场理论中，拷贝场声信号s(在地声参数m条件下，利用水声数值模型计算得到)通常难以与观测声压so完全相同，可利用似然函数描述两者的相似程度L(so|m)，即

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

(2)

于是Bayes地声参数反演问题转化为求取地声参数m的PDFP(m|so)。反演结果通常以参数的一维边缘概率密度表示：

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

(3)

(4)

(5)

(6)

(7)

式中：N为VLA阵元数；F为数据的频率点数；vf与VLA信噪比(SNR)相关；Bf为常用的Bartlett处理器，

(8)

2 Multi-AG采样算法

MCMC随机模拟思想主要包括蒙特卡洛积分、马尔可夫链及Metropolis-Hastings抽样算法，Gibbs采样算法可替代Metropolis-Hastings算法以提高在高维参数空间下计效率。

2.1 MCMC算法

2.1.1 蒙特卡洛积分

MCMC算法通过在各参数取值区间内随机均匀采样，则任一组样本参数mi的PDF为

(9)

式中：Q为样本数。通常样本参数mi的维度高(为多个地声参数反演的联合后验概率)，指示反演结果不够直观，可利用(3)式、(4)式、(5)式得到其一维边缘概率密度、最大PDF、平均PDF.

当随机抽取样本达到一定数量时，蒙特卡洛积分可得到地声参数的一种无偏、非对称估计结果。

2.1.2 马尔可夫链

假设已生成的一连串随机参数样本{m0,m1,…,mt}，则t+1时刻的状态概率为

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

(10)

(10)式表示了马尔可夫链的基本原理，即任一时刻t+1的状态转移概率只与前一时刻t相关，而与已生成的样本{m0,m1,…,mt}无关。

随机抽取先验概率密度分布服从π(·)的地声参数变量，生成易于实现且不可约遍历链{mt}，使抽样变量平稳地服从π(·)分布。则在Bayes匹配场地声参数反演过程中，随着样本数量的增加，马尔可夫链逐渐忘记样本初始化状态m0，PDF的反演结果将最终服从某种分布π(·).

2.1.3 Metropolis-Hastings算法

Metropolis-Hastings算法是MCMC算法中随机模拟与决定性算法相结合使用的一种常用解决方案，其算法的主要流程是：

1)初始化待反演地声参数样本m0；

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

(11)

4)重复上述步骤，直至产生预设的T个样本为止。

2.2Multi-AG采样算法

Gibbs采样将高维样本化为一系列一维取样，可避免Metropolis-Hastings算法在计算转移概率α(通常α<1)偏小所导致的拒绝率过大的问题，进而有效地提高运行效率。

(12)

P(mi|so)=

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

(13)

由(13)式可知，Metropolis-Hastings算法需经P-1重积分方可得到参数的一维边缘概率密度，而Gibbs采样只需一重积分即可，即

(14)

Gibbs采样按待反演参数维度逐个采样判别，由于各参数对处理器敏感性存在差异，不同参数样本拒绝率差别显著。图1为Gibbs采样前180次参数判别跳转轨迹，表1为Workshop’97浅海环境基准模型参数[15]。受样本拒绝率控制，c1的采样判别次数远大于ρs、αs的采样判别次数，导致前者的反演效果明显优于后二者。为解决各参数反演效果差异问题，通常的做法是扩大采样数量，这也极大地增加了计算量。

事实上，增加采样数量对c1反演结果的影响较小，c1收敛速率快，在采样反演前期即可完成收敛。而将所有参数共同进行反演，c1既增加了计算量，也限制了其他参数的反演采样次数。多步Gibbs采样首先完成对敏感性最高的参数反演，利用反演结果计算其均值并视为常量，进而逐级反演敏感性更低的参数，可极大地减少计算量。

图1 Gibbs采样参数跳转轨迹Fig.1 Parameters skipping track of Gibbs sampling method

参数实际值声源深度zs/m20声源水平距离rs/km1水体深度D1/m100水体吸收系数α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沉积层衰减系数αs/(dB·λ-1)0.2基底声速c3/(m·s-1)1689基底密度ρb/(kg·m-3)1700基底衰减系数αb/(dB·λ-1)0.2

参数搜索空间对Gibbs采样反演结果的影响较大，特别是在敏感性较低的参数反演时，不同的搜索空间得到的参数均值、均方差差异显著。图2为ρs一维PDF反演结果，分别设置参数采样空间为[1 300 kg/m3,1 600 kg/m3](见图2(a))、[1 400 kg/m3, 1 700 kg/m3](见图2(b))时，Gibbs采样算法反演均值(±均方差)分别为1 498.6 kg/m3±53.97 kg/m3、1 511.3 kg/m3±62.07 kg/m3，二者差异较大，且与实际值1 500 kg/m3偏离较远。此现象主要源于敏感度较低的参数对处理器影响微弱，造成其整个搜索空间内接受率高，进而导致其一维PDF较均匀所致。

图2 参数ρs一维边缘PDF分布Fig.2 Marginal posterior probability function of parameter ρs

模拟退火(SA)算法是概率算法的一种。该方法结合热力学与统计学理论，将参数空间内各点的概率视为分子动能，并计算前后两种状态转化的目标函数差，进而依据判定准则判决是否接受该样本。

定义退火Gibbs采样算法的目标函数差为

(15)

同样的搜索空间，当采用退火Gibbs采样算法时，ρs的反演结果分别为1 501.1 kg/m3±28.13 kg/m3(见图2(c))、1 500.7 kg/m3±28.31 kg/m3(见图2(d))，二者间的微小差异源于随机采样所致，而与搜索空间无关。

Multi-AG采样算法结合退火Gibbs算法，逐步完成高维地声参数反演，Multi-AG采样算法参数反演算法流程如图3所示。图3中，参数敏感性分析①以Bartlett处理器为量化依据，Multi-AG反演根据参数敏感性等级划分为n步，第m(2≤m≤n)步将前m-1步中敏感性较高的参数反演均值作为常量输入，即②和③，反演第1步利用Gibbs采样算法，后续步骤使用退火Gibbs采样算法。

图3 Multi-AG采样算法参数反演流程Fig.3 Flow chart of parameters inversion by multi-annealing Gibbs sampling method

3 地声参数反演

水声调查数据的系统性(缺实测地声环境数据)限制了对地声参数反演结果的验证，实际海洋环境的不确定性降低了Bayes匹配场的反演精度。为更好地实现与检验Multi-AG采样算法的地声参数反演效果，本文采用Workshop’97[15]提供的浅海环境基准模型进行地声参数反演(见图4)，该模型在水声学理论研究中具有重要的参考价值[12,16]，利用其进行反演得到的结果可与已有的研究成果进行对比验证。

图4 Workshop’97浅海环境基准模型示意图[15]Fig.4 Sketch of the benchmark of shallow water environment by Workshop’97[15]

[15]中测试样例1，相关参数物理含义及取值如表1所示。Workshop’97浅海环境基准模型为水平不变的两层海底海洋环境，考虑沉积层声速垂向变化，无限大基底层声速恒等于1 700 m/s，声速剖面同文献[15]中的图2.

VLA含20个垂向间距4 m的阵元，最浅单元位于3 m，无指向性声源频率为100 Hz. 利用Kraken简正波数值模型[17]计算声场(其中模型代码下载地址为：http:∥oalib.saic.com)，垂向计算点分辨率1 m，水平10 m，SNR为20 dB.

3.1 地声参数敏感性等级划分

表1中含8个地声参数，分别为αs、ρs、c1、c2、D2、αb、ρb、c3，地声参数反演过程中，声场对各参数的敏感性不同，进而导致Bartlett处理器变化尺度各异，Bayes匹配场对敏感性较低的参数反演效果较差，且增大了运算量，故进行地声参数敏感性分析尤为重要。

基于单因子变量原则，分别计算各地声参数变化时声场与参数实际值下声场1 km处的Bartlett处理器PDF，PDF分布特征如图5所示。

图5 各地声参数变化对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几乎服从均匀分布，对Bartlett影响可近似忽略。沉积层覆盖于基底之上，与声波相互作用显著，其声速大小影响声阻抗造成的海底反射损失是底边界对声场作用的主要来源，故c1对Bartlett处理器的敏感性最高。

赵航芳等[18]利用Workshop’93典型浅海环境，分析了各水声环境参量不确定性对Bayes匹配场性能的影响，其研究结论与本文符合较好。

图6 各地声参数敏感性变化曲线Fig.6 Sensitivity curves of geoacoustic parameters

经分析，可将8个地声参数划分为4个敏感性等级。第1级(敏感性等级最高)含参数c1，比例在3%以内；第2级含c2、D2，比例分别为9%、17.5%；第3级含αs、ρs、c3，比例分别为59%、75%、42%；第4级含αb、ρb，对Bartlett处理器几乎无影响，比例大于99%.

参考已有的研究成果[10,12,19]反演得到的αb、ρb参数PDF在其搜索空间内几乎服从均匀分布，效果差，故下文不对αb、ρb进行反演，按Workshop’97中的定义相应地将其设置为常数。

3.2 退火方案

退火算法的目的在于提高参数收敛性能，消除敏感性较弱的参数因搜索空间设置所导致的反演误差。但T太小或退火速率过快都会导致样本拒绝率过大而无法采样。

本文依据各参数的PDF及反演结果所预期的参数均方差范围，分别制定退火方案。对于参数D2、αs、ρs、c3，T0分别取0.7 ℃、0.2 ℃、0.1 ℃、0.5 ℃. 通过数值实验并分析其PDF发现，随着采样数的增长，在采样前期(1/2倍总样本数)，控制T从1 ℃以总采样数/20的速率降低至T0，可将参数反演均方差控制在相应搜索空间宽度的15%左右，既提高了参数的反演精度及收敛性，又保证了退火Gibbs采样算法较高的样本接受概率。

3.3 Multi-AG法地声参数反演

基于表1中的环境参数配置，通过Kraken模型计算1 km处对应VLA深度的声压分布，并假定其为实测声压数据。依据3.1节和3.2节分析得到的参数敏感性等级划分及退火方案，在兼顾水声环境参数模型整体性的前提下，使地声参数的采样范围最大化，设定参数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]，在采样空间内进行随机均匀采样。

通过106次Metropolis-Hastings算法、105次Gibbs采样算法、7×104次Multi-AG采样算法(第1步、第2步、第3步分别计算104次、2×104次、4×104次)计算得到地声参数PDF如图7～图9所示，图中蓝色实线为参数实际值位置。表2总结了Metropolis-Hastings算法、Gibbs采样算法、快速Gibbs采样(FGS)算法、Multi-AG采样算法4种算法的反演结果及运算量。

分析表2可知，105次Gibbs采样算法反演得到的各参数均值及均方差与105次Metropolis-Hastings算法相当，表明Gibbs采样算法可极大地提高运算效率。然而，各参数的反演效果随各参数敏感性逐渐变差，敏感性最高的参数c1(见图7(a)和图8(a))反演结果分别为1 535 m/s±4.67 m/s、1 535 m/s±4.82 m/s，接近实际值1 530.4 m/s，且收敛性较好；敏感性低的参数αs(见图7(d)和图8(d))、ρs(见图7(e)和图8(e))、c3(见图7(f)和图8(f))的PDF几乎接近均匀分布，偏离实际值较远，旁瓣高，且均方差过大。

图7 Metropolis-Hastings方法地声参数反演结果Fig.7 Inversion results of geoacoustic parameters by Metropolis-Hastings method

图8 Gibbs采样算法地声参数反演结果Fig.8 Inversion results of geoacoustic parameters by Gibbs sampling method

对于多维参数共同反演，各参数对Bartlett处理器的敏感性差异较大，在参数采样判别循环过程中，敏感度较高的参数循环判别次数多，限制了敏感度较低的参数随机采样次数。且较低的敏感性使得Bartlett处理器的效能降低，随机误判出现频繁，无法有效地向参数真实值收敛，进而出现近似均匀分布的PDF. 故可知，在敏感性较低的参数整个搜索空间内都具有相当的样本接受概率，导致参数反演均值接近其搜索空间中值，使得参数搜索空间设置的差异影响到参数的反演结果。参数αs搜索空间为[0.001 dB/λ,0.500 dB/λ]，反演均值为0.280 dB/λ，接近搜索空间中值0.250 dB/λ，而与实际值0.200 dB/λ差异较大。

Multi-AG采样算法分3步分别对6个参数逐步进行反演，第1步循环104次，接受样本622个，得到c1的反演结果为1 535±5.13 m/s，其一维边缘概率密度如图9(a)所示。对比106次的METROPOLIS-HASTINGS算法与105次的Gibbs采样算法，三者反演效果差异微小，说明第1步Multi-AG算法可快速高精度地反演敏感性最高的参数c1.

利用第1步得到的c1PDF计算其均值(1 535 m/s)，并将c1视为海洋环境模型常量参数，进而对αs、ρs、c2、D2、c3进行退火Gibbs采样反演。第2步采样2×104次，接受1 553个样本，可得到敏感性较低的参数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采样算法更接近实际值。

图9 多步退火Gibbs采样算法地声参数反演结果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)总样本数接受样本数实际值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算法，对αs、ρs、c3进行反演。第3步采样4×104次，接受4 980个样本，各参数的一维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，单峰现象明显，且搜索空间边界附近的PDF值近似为0，有效地消除了旁瓣，表明参数反演精度显著提高。

值得一提的是，对于多维地声参数反演，Metropolis-Hastings算法、Gibbs采样算法、Multi-AG 3种方法对敏感度最高的参数(c1)反演效果相近。

对比Multi-AG算法各步骤的总样本数与接受样本数，第2步(5个反演参数)采样数是第1步(6个反演参数)的两倍，接受样本数大于其2倍；第1步、第3步(3个反演参数)采样数相同，第3步接受到的样本数大于第2步。基于水声数值模型计算声场用于反演地声参数，是一个强非线性问题[20]，Stan等[12]研究了各参数间相关性对反演结果及运算量的影响，指出不确定参数越多，反演所需计算量越大。本文通过敏感性分析，采用多步反演策略，减少了计算过程中的参数数目，提高了精度，较常规的Gibbs采样算法进一步减少了计算量。

退火算法可有效增强参数对Bartlett处理器的敏感性，提高其收敛性能。对于敏感性较高的参数c1、c2，无需退火，直接进行Gibbs采样判别即可得到较好的反演效果。从图10(a)和图10(b)采样轨迹可看出，参数c1、c2接受的样本收敛性较好，所接受的样本值接近实际值。

Multi-AG第1步采样含6个反演参数，敏感性弱的参数αs(见图10(c))、ρs(见图10(d))、c3(见图10(e))接受的样本轨迹分散。利用分布反演策略后，Multi-AG第3步只含3个参数，进而对其进行退火Gibbs采样，随者T的降低，接受的样本逐渐向各自的实际值聚集，可有效提高反演精度并减小结果的均方差(见图10(f)、图10(g)和图10(h))。

3.4 结果验证

Stan等[12,16]基于Workshop’97浅海环境基准模型，利用FGS算法反演地声参数，相关参数的反演结果如表2(引自文献[12]中的表1)所示。对比分析参数实际值、FGS算法及本文的Multi-AG算法可知，Multi-AG算法经过多步反演得到的参数αs、ρs、c2、D2反演精度更高。其中，FGS算法反演得到的αs、ρs均方差分别为0.076 dB/λ、90.1 kg/m3，Multi-AG算法反演得到的αs、ρs均方差分别为0.037 dB/λ、59.2 kg/m3，表明Multi-AG算法反演得到的参数收敛性更好。

4 结论

利用Gibbs采样算法代替常规的Metropolis-Hastings算法进行高维Bayes匹配场地声参数反演，可在一定程度上减小运算量。然而，各参数对反演处理器敏感性的不同，使得反演获取各参数稳定的一维PDF所需的计算量各异，且敏感性较小的参数搜索空间差异会在一定程度上影响参数的反演结果。

本文综合Gibbs采样算法、退火算法及地声参数对Bartlett处理器的敏感性差异，提出Multi-AG采样算法。通过采用Workshop’97浅海环境基准模型，对比分析Metropolis-Hastings、Gibbs采样算法、FGS算法、Multi-AG采样算法4种算法地声参数反演所需计算量、均值、均方差，结果表明：

1)Multi-AG参数反演所需计算量最小，说明分布反演方案进一步提高了Gibbs采样算法的反演效率。

2)Multi-AG参数反演均值与实测值最为接近，均方差最小，敏感性较低的参数反演改进效果尤为明显，证实了退火方案可有效提高参数反演精度。

本文Multi-AG算法用于地声参数反演取得了较好效果，该算法也可用于Bayes匹配场其他海洋环境参数(水深、声速剖面等)的反演。由于缺少实测地声参数数据，本文将反演结果与已有的研究文献[12]进行了对比验证，这也是下一步将需要深化的研究方向。

参考文献(References)

[1] 潘长明, 高飞, 孙磊, 等. 浅海温跃层对水声传播损失场的影响 [J]. 哈尔滨工程大学学报, 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] 李倩倩. 不确知海洋环境下的贝叶斯匹配场定位性能研究[D]. 北京: 中国科学院大学, 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] 王奇, 王英民, 苟艳妮. 不确定环境下后验概率约束的匹配场处理 [J]. 兵工学报, 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] 李倩倩, 李整林, 张仁和. 不确知海洋环境下的贝叶斯声源定位法 [J]. 声学学报, 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] 赵航芳，李建龙，宫先仪．不确实海洋中最小方差匹配场波束形成对环境参量失配的灵敏性分析 [J]．哈尔滨工程大学学报, 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

国家自然科学基金项目(41406004)

高飞(1988—)，男，助理工程师。E-mail:gfei88_lgdx@163.com； 潘长明(1962—)，男，高级工程师。E-mail: pcming62@163.com

P733.23

A

1000-1093(2017)07-1385-10

10.3969/j.issn.1000-1093.2017.07.017