改进粒子滤波在海底地声参数反演中的应用

吴伟文，任群言，鹿力成，马 力

(1.中国科学院声学研究所中国科学院水声环境特性重点实验室，北京 100190；2.中国科学院大学，北京 100049)

0 引言

声信号在水中具有良好的传播性能，因此常用来进行水下目标探测或信息传递。浅海环境中，声传播受海底的影响较大，准确的地声参数有助于分析浅海环境中声传播特性，因此地声参数估计是海洋声学中的一个基础研究问题。

地声参数可以通过直接测量法[1]获得，但此方法时间成本和经济代价都很高。反演算法如匹配场处理[2]是比直接测量法更有效率的方式。匹配场处理将地声参数反演问题作为优化问题，其目的是找到一个模型，此模型包含一组地声参数(假设模型的真实地声参数向量包含在参数集合中)，使前向模型产生的数据与观测数据达到最优匹配。

匹配场处理可分为两步。(1)选择目标函数，对观测数据与前向模型生成数据之间的误差进行衡量；(2)基于最小化目标函数的准则搜索最优的参数向量。目标函数中待匹配的数据，一般可选择为对地声参数变化较敏感的物理量。例如，海底反射系数随沉积层分层特性、声速、密度的变化而明显变化，因此海底反射系数[3]可以作为目标函数中待匹配的物理量。本文将海底反射系数作为待匹配的物理量，利用优化算法搜索地声参数向量。由于地声反演问题具有维度高，局部最优多，前向模型非线性等特点，所以参数搜索方式一般选遗传算法[4]、模拟退火[4]、MCMC[5]采样等全局算法。上述优化算法一般称之为批处理反演算法，不仅搜索过程费时而且一批数据只能生成一个反演结果，反演结果之间互相独立，无法高效处理环境随距离变化时的参数估计问题。

为了利用相邻两段距离之间地声参数的联系，减少运算量，可以将序贯寻优算法如粒子滤波应用到地声参数的反演中。在地声参数随距离变化缓慢的情况下粒子滤波算法已经取得了良好的效果[6]。但是在地声参数随距离变化剧烈，比如沉积层分层个数改变的情况下，传统粒子滤波可能会失效。为了解决此问题，Dettmer等[7]将带桥接的粒子滤波应用到地声参数反演中，通过桥接过渡剧烈变化的目标分布，以实现地声参数随距离变化剧烈时的反演，并通过仿真数据证明了此方法的有效性。本文在粒子滤波中增加桥接的同时改进了粒子的采样方式，由传统的固定方差采样改进为随机方差采样，同时对沉积层厚度的采样采用高维与低维混合采样的方式。将此方法应用到地声参数反演中，通过仿真分析得出，此方法在地声参数随距离变化剧烈的情况下仍可得到较稳定的反演结果。

本文利用反射系数构造似然函数，分析其对沉积层厚度、声速、密度的敏感性；给出了改进粒子滤波的算法流程；然后进行仿真数据处理；得到改进粒子滤波在地声参数随距离变化剧烈的情况下仍可得到较稳定的反演结果的结论。

1 海底反射系数

1.1 反射系数计算

海底反射系数对沉积层厚度以及声速都较为敏感，因此常作为地声参数反演中待匹配的物理量。反射系数R为反射到海水的所有反射波与入射波的比值，基于射线声学理论，平面波反射系数计算原理如图1所示[8]，R的计划公式为

图1 沉积层数为1时的海底反射图Fig.1 Bottom reflection diagram when the number of sediment layer is 1

其中：ρ1、ρ2、ρ3分别为海水、沉积层和基底的密度；c1、c2、c3分别为海水、沉积层和基底的声速(当考虑海底声衰减系数α2、α3时，海底声速为复数)；h2为沉积层厚度；R12为声波由海水入射到沉积层的反射系数，T12为声波由海水入射到沉积层的透射系数，R21为声波由沉积层入射到海水的反射系数，T21为声波由沉积层入射到海水的透射系数，R23为声波由沉积层入射到基底的反射系数，T23为声波由沉积层入射到基底的透射系数，R32为声波由基底入射到沉积层的反射系数，T23为声波由基底入射到沉积层的透射系数。相位延迟为

其中：k2=2πf/c2，f为声源频率,θ2为折射波方向(已知入射波方向θ1，声速c1和c2，可根据斯涅尔定律计算θ2)。根据R12=-R21和T12T21=1-R212，反射系数R进一步可表示为

多层沉积层的反射系数可以依公式(3)递推获得，也可通过OASR工具箱[9]计算获得。

1.2 敏感性分析

反射系数随沉积层层厚，声速，密度的变化而变化，本文用反射系数构造如式(4)所示的似然函数：

其中：xt为第t帧的地声参数向量，包含沉积层深度、声速、密度、衰减系数；yt为第t时刻的观测值，对应于反演问题中的第t帧的反射系数观测向量(包含N个掠射角)；g为测量函数。研究似然函数对沉积层声速、密度、衰减系数、厚度的敏感性，对反演的可行性以及反演结果分析具有重要作用。下面在四层沉积层的情况下，对似然函数进行敏感性分析。

令声源频率分别为500 Hz和1 000 Hz，仿真参数如表1所示。

表1 仿真信道参数Table 1 Parameters of the simulation channel

图2 Φ随ci、ρi、αi、hi各参数变化的曲线图Fig.2Variation curves of Φ with parameters ci、ρi、αi、hi

可以清楚地看到，Φ对声速和沉积层深度敏感，对密度、衰减系数不敏感。同时可以看到Φ随h变化时，局部最优明显。为了防止算法在对沉积层深度的搜索时陷入局部最优，本文提出了沉积层厚度混合搜索方式，此方法将在第2.2节介绍。

2 改进粒子滤波

2.1 带桥接的粒子滤波

在随距离变化的地声反演问题中，可以将地声参数分为很多帧。本文的帧是指一小段范围，在此范围内地声参数可以近似认为是不变的。由于每一帧之间地声参数具有相关性，所以可以将粒子滤波等序贯寻优方式应用到地声反演中。根据粒子滤波原理建立状态方程和观测方程分别为

式中：xt是第t时刻的状态，在地声反演问题中，即为第t帧的地声参数向量，包含沉积层厚度、声速、密度、衰减系数；xt-1为上一时刻的状态，对应于反演问题中上一帧的地声参数向量；yt为第t时刻的观测值，对应于反演问题中的第t帧的反射系数观测向量(包含N个掠射角)；f为状态转移函数，g为测量函数；g(xt)为以xt为地声参数向量，代入反射系数计算公式，所求得的反射系数向量；nt和vt分别是状态方程与观测方程的噪声，vt一般为高斯噪声，则yt即为以g(xt)为均值的高斯分布。假设各掠射角的高斯噪声独立，似然函数L可以表示为

粒子滤波可能会由于似然变化较大而出现粒子枯竭的情况。为了应对粒子枯竭，可以采用重采样的方式，重采样是指维持总粒子数量不变的情况下，复制权重大的粒子，同时舍去权重非常小的粒子，是一种有效防止粒子枯竭的方法。

将粒子滤波应用到地声参数反演中，在地声参数随距离缓慢变化时，传统粒子滤波具有良好的效果。但当地声参数随距离剧烈变化时，似然函数变化剧烈，即使使用了重采样，传统粒子滤波仍可能会出现粒子枯竭[11]的情况，无法正确收敛到目标分布。为了使粒子滤波可以在地声参数剧烈变化的情况下收敛到目标分布，有文献使用了带桥接的粒子滤波的方法。使用桥接可以使粒子滤波结果缓慢收敛到目标分布，实现目标分布之间的平滑过渡。带桥接的粒子滤波原理为[12]

其中，1＞α1＞α2…＞αM＞0，M为桥接层数。桥接的作用是使目标分布从逐渐收敛到，即为重要性采样概率密度，当m=1时，αM+1=1，。当m=0时，，即为最终的目标分布。通过桥接，目标分布从q(x1:t)逐渐收敛到，即最终收敛到参数后验概率密度。桥接原理如图3所示，每一步桥接中都使用重采样从q(x1:t)平缓过渡到，以防止粒子枯竭。

图3 桥接原理图Fig.3 Bridging priciple diagram

2.2 方差随机、单多层层厚混合采样的粒子滤波

粒子滤波解决地声反演问题时，状态方程式(5)可以表示为式(15)形式：

其中：nt一般选择为高斯分布，对于地声反演这种局部最优多的优化问题，选择高斯分布可能会因随机性的不足而导致反演结果易陷入局部最优中。尤其是对于地声参数随距离剧烈变化的情况，即xt与xt-1相距甚远的情况，这时候采用高斯分布噪声，更可能由于采样随机性的不足而导致反演陷入局部最优。容易联想到增大方差来增加粒子随机性的方式，但如此一来由于搜索精度不足，也容易使反演问题陷入局部最优。有文献使用在粒子滤波后增加马尔可夫链蒙特卡洛采样[9]的方法，利用马尔科夫链蒙特卡洛方法可以在目标分布附近采样，从而增加粒子多样性。但是马尔可夫链蒙特卡洛方法计算量巨大，代价高昂。

为解决此问题，在每一步中令nt的方差随机，使用随机方差采样。令z=0+nt，用方差固定和方差高斯这两种高斯采样方式，分别进行100 000个点的采样。方差固定采样是在方差固定为1的零均值高斯分布中采样，随机方差采样相对于固定方差采样，其方差由固定的1变为随机。随机的方差可以在方差为1的零均值高斯分布中采样获得。如图4所示为两种采样时，变量z的样本分布情况。

图4 两种采样方式(方差固定和方差高斯)时，变量z的样本分布Fig.4 Sample distributions of two sampling methods(fixed variance and Gaussian variance)

从图4中可以看到，在方差为高斯分布的情况下，样本在中心更集中，同时，其样本又有较好的两端扩展能力，即在距离中心较远的地方仍然具有相当数量的样本点。在距离相关的地声参数反演问题中，每帧之间的地声参数既可以缓慢变化，也可以剧烈变化，所以方差随机的高斯采样方式十分适合距离相关的地声参数反演问题。本文中的方差随机采样方法，同样达到了增加粒子多样性的效果，进一步防止了粒子枯竭，有助于目标分布收敛。相较于固定方差的采样，3.3节中的模拟结果验证了方差随机采样可以使反演结果具有更好的稳定性。

对于采样的维度，传统粒子滤波中一般使用多维高斯分布采样，采样维度与待反演地声参数的维度相同。但是根据敏感度曲线，沉积层厚度有局部最优，使用高维高斯分布采样虽然保证了全局性能，但由于搜索精度不够，不容易找到最优解。如果在搜索中混合低维采样，相当于提高了搜索精度，但同时全局性能不足，容易陷入局部最优。本文提出了单层层厚搜索与多层层厚搜索混合的方式，层厚混合搜索的粒子滤波采样时使用了沉积层厚度的单层多层混合改变的搜索方式，单层搜索指所有沉积层的层厚只改变任意一个，多层搜索指所有沉积层的层厚全部改变，采用50%概率所有层层厚都变化，50%概率只改变任一层层厚的方式。此方法在增加了全局搜索能力的同时，增加了搜索结果的精度。3.2节的仿真结果验证了此方法的有效性。

基于地声参数反演的特点，本文使用了带桥接、方差随机、单层层厚与多层层厚混合搜索的粒子滤波方法，在提高搜索精度的同时保证了解的全局性，更好地处理了地声参数随环境剧烈变化的问题。

3 仿真结果分析

3.1 仿真环境及传统粒子滤波反演结果

采用的仿真环境设置如下：地声参数随距离变化的100帧海底环境，在第10帧和第11帧之间沉积层分层数目改变，同时沉积层声速和密度也有较大变化。在第20帧和第21帧之间以及第80帧和第81帧之间沉积层分层数目也发生改变。第21帧沉积层由4层变为3层，第81帧沉积层变为2层。第1帧共4层沉积层，地声参数参见表1中。在所有帧中，沉积层总厚度设置为固定的10 m。粒子数量设置为10 000个。水听器数量设置为13个，收集在36°～79°之间共13个掠射角的反射系数。发射声源发射两个单频信号，频率分别为1 000 Hz和500 Hz。实验装置如图5所示。

图5 实验布设示意图Fig.5 Layout diagram of experiments

模拟的真实地声参数分布如图6～8所示。利用传统粒子滤波算法反演结果如图9～11所示。图9中，每帧声速用10 000个粒子的均值表示。图10中，每帧密度用10 000个粒子的均值表示。图11中，每帧衰减系数用10 000个粒子的均值表示。

图6 模拟真实环境声速Fig.6 Simulated sound speed in real environment

图7 模拟真实环境密度Fig.7 Simulated density in real environment

图8 模拟真实沉积层衰减系数Fig.8 Simulated attenuation in real environment

图9 传统粒子滤波沉积层声速c反演结果Fig.9 Acoustic velocity c inversion results of sedimentary layer obtained by inversion with traditional particle filtering

图10 传统粒子滤波沉积层密度ρ反演结果Fig.10 Density ρ of sediment layer obtained by inversion with traditional particle filtering

图11 传统粒子滤波沉积层衰减系数α反演结果Fig.11 Attenuation coefficient α of sediment layer obtained by inversion with traditional particle filtering

可见传统粒子滤波在地声参数剧烈变化时，无法正确估计地声参数，获得的地声参数有较大误差。

3.2 层数混合采样

使用带桥接的粒子滤波，一定程度上能够处理地声参数随距离剧烈变化的反演问题。但是仅仅使用桥接并不够稳定，粒子滤波会出现无法收敛或收敛到局部最优的情况。基于地声反演的特点，本节将2.2节提出的混合层厚搜索的方式，与只有桥接但没有混合层厚搜索的粒子滤波做比较。

由于仿真环境设置为第11帧变化最剧烈，第11帧将是粒子滤波最具挑战性的1帧。接下来对比这三种采样方式在第11帧时的性能差异。随机生成30个不同的仿真环境分别用单层层厚改变、多层层厚均改变、层厚混合改变的采样方式进行带桥接的粒子滤波，记录通过这三种粒子滤波后，第11帧的反射系数最小均方误差(反射系数最小均方误差即指通过多层桥接后得到的最小均方误差)。为了使均方误差不至于太小，这里的均方误差没有除以样本数量(26)，而是直接用|R-R̂|2表示，其中R为观测反射系数，R̂为用所有粒子均值求得的反射系数。各掠射角的反射系数的标准差设置为0.01。将这三种方式所得的30次实验的最小均方误差按照从大到小排列，得到如图12所示曲线。

图12 三种采样方式得出的第11帧的最小均方误差Fig.12 The minimum mean square errors of the 11th frame obtained by the three sampling methods

从图12可以看出，层厚的混合改变采样方式相较于层厚均改变的采样方式有更小的均方误差，前者更容易获得全局最小的反演结果；相较于层厚单层改变的采样方式反演效果也有小幅度提升。粒子滤波采用层厚混合采样的方式时，第11帧地声参数有更大的概率收敛到全局最优。

3.3 方差随机搜索

状态方程中的噪声一般设为方差较大的高斯扰动，设置此方差较大的原因是为了增加粒子滤波的搜索范围，使地声参数至少有从某一帧采样到下一帧真值的可能性。为了对比随机方差的粒子滤波与固定方差的粒子滤波对反演效果的影响，进行了17次实验，分别使用随机方差(均值为0，方差为1的高斯分布)采样、固定方差0.6(状态方程噪声方差的0.6倍，下同)采样、固定方差0.8采样、固定方差1采样、固定方差1.2采样、固定方差1.4采样。沉积层厚度采用层厚均改变的采样方式。得到的最小均方误差结果如图13所示。

由图13可以看出，随机方差的采样方式大大提高了第11帧地声参数收敛到全局最优的可能性。注意到方差过大(1.4)或过小(0.6)将会出现第11帧最小均方误差大于1的情况，说明此时粒子滤波效果已十分不佳。同时还可以看出，所有固定方差的采样方式，最小均方误差曲线都不足17个点，这是由于某些实验中，固定方差的采样方式已经不能使粒子滤波稳定运行，即没有搜索到一个最优或局部最优点，这将导致粒子滤波失效。而随机方差采样的方式则可大幅提高粒子滤波的稳定性，几乎每次粒子滤波都能正常收敛。

图13 不同方差采样得出的第11帧的最小均方误差Fig.13 The minimum mean square errors of the 11th frame obtained by different variance sampling methods

为了减少运算量，可以令桥接数量可变，桥接数量最大为9层，令反射系数均方误差基本收敛时停止桥接。实验采用3.2中的沉积层层厚混合改变的采样方式，统计所用桥接数量。图14是方差固定和方差随机时桥接数量随帧数的变化，图15是将这些桥接数量按照从大到小排列的结果。

图14 方差固定和方差随机时桥接数量随帧数的变化Fig.14 Variation of the number of bridges with the number of frames when variance is fixed or random

图15 当帧数增加到100时桥的数量按降序排列图Fig.15 Diagram of the number of bridges arranged in descending order when the number of frames is increased to 100

从图15可以看出，方差随机时平均桥接数量小于方差为固定值时的平均桥接数量。随机方差采样平均桥接数量为4.65层，固定方差采样平均桥接数量为5.43层，方差随机时算法收敛速度快。可以看到在地声参数随距离变化剧烈时需要的桥接数量较多，变化缓慢时需要的桥接数量较少，这是可以预见的。

3.4 改进粒子滤波反演结果

在粒子滤波中使用桥接的同时，同时使用状态方程方差随机以及沉积层层厚混合采样的方式，得到如图16～18的反演结果。图16～18中，每一帧的声速、粒子、密度均用粒子均值代替。

图16 改进粒子滤波声速c反演结果Fig.16 Acoustic velocity c of sediment layer obtained by inversion with improved particle filtering

图17 改进粒子滤波密度ρ反演结果Fig.17 Density ρ of sediment layer obtained by inversion with improved particle filtering

图18 改进粒子滤波声衰减系数α反演结果Fig.18 Attenuation coefficient α of sediment layer obtained by inversion with improved particle filtering

由反演结果可见，改进粒子滤波在沉积层参数随距离剧烈变化的情况下，仍可较准确地反演地声参数。真实环境中，沉积层层数会发生改变。但在改进粒子滤波反演结果中，沉积层层数一直使用的都是4层。可以看到在20～80帧，以及80～100帧中，会出现某些层极薄且和上层或下层声速、密度相近的情况，此时反演结果中这些层将不含明显的分层显示，所以一直使用4层沉积层仍能取得较好的反演结果。改进粒子滤波在第11帧、第21帧、第81帧声速和密度都有较好的反演结果。改进粒子滤波在第11帧时反射系数拟合效果如图19所示，反射系数曲线拟合效果较好。

图19 第11帧反射系数R拟合效果Fig.19 Fitting results of the 11th frame reflection coefficient R

3.5 反演结果的边缘概率密度

为了获得某个地声参数的分布情况，分析了反演结果的边缘概率密度。图20～22为第11帧声速、密度和衰减系数反演结果的边缘概率密度图。图20～22中粉色虚线为真值。

从图20～22中可见声速密度和衰减系数均在真值附近，其中声速的不确定性较小，密度和衰减系数的不确定性大，这是由于似然函数对声速的敏感性高，对密度和衰减系数的敏感性低。不确定性分析对应于第1节敏感性的分析，证实了反演结果的可靠性。

图20 声速c边缘概率密度图Fig.20 Marginal probability density diagram of sound velocityc

图21 密度ρ边缘概率密度图Fig.21 Marginal probability density diagram of density ρ

图22 衰减系数α边缘概率密度图Fig.22 Marginal probability density diagram of attenuation coefficient α

4 结论

针对海底地声参数随距离变化的地声参数反演问题，本文提出了一种改进的粒子滤波序贯估计方法。在带桥接的粒子滤波基础上，使用了沉积层厚度混合变化采样的方式，增加全局搜索能力的同时使搜索结果更加精确。在状态方程中使用噪声方差随机采样，增加了粒子的多样性，同时减少了所需的桥接系数数量。相较于传统粒子滤波，该方法在地声参数随距离剧烈变化的情况下，依旧能较稳定地反演地声参数。相较于只有桥接的粒子滤波，本文方法的反演结果更加稳定，计算效率更高。

利用不同掠射角和不同频率的反射系数作为待拟合的观测数据，采用改进粒子滤波的方法进行地声参数反演。仿真结果证明了在沉积层参数剧烈变化的环境中，本文方法仍然具有良好的效果。