一种获取海底地声参数与沉积层分层结构的序贯方法

技术领域

本发明涉及海底反演、海底遥感的技术领域，特别涉及一种获取海底地声参数与沉积层分层结构的序贯方法。

背景技术

获取准确的海洋环境信息是海洋资源开发以及国防建设的基础，而局部海域的海底地声参数和沉积层分层结构是研究海底地质变迁和声纳设备应用的重要条件。传统的提取海底地声参数与沉积层分层结构的方法是采用柱状原状取样设备，利用机械压力或者重力定点取样，通过实验室对底质进行分析，进而得到沉积层声学参数与分层结构信息。但是，这种提取方法的时间成本和经济成本昂贵，而且相对于广阔的海洋，取样范围极其有限，不能获得广阔海域的地声参数信息及沉积层分层结构信息。

因而，根据不同海底结构对声传播的影响，应用声学遥感的方法如匹配场反演可以快速、准确的获得海底地声参数与沉积层分层结构。对于海底特性与距离无关的情形，地声参数与沉积层分层结构都不随距离变化，匹配场反演中使用批处理算法具有较好的反演效果。对于海底特性随距离变化而沉积层层数不变的情形，使用序贯方法可以较好反演地声参数与沉积层厚度。但是，对于海底特性与沉积层层数都随距离变化的地声反演问题，使用粒子滤波等序贯方法无法获得沉积层分层信息。为了获取沉积层分层信息，在粒子滤波后需要使用rjMCMC采样等运算量极大的方法反演沉积层层数，但此方法反演速度慢，对于实时性要求较高的系统，此方法较难实现。

发明内容

本发明的目的是为了解决现有的提取方法所存在的上述缺陷，并提供了一种快速获取海底地声参数与沉积层分层结构的序贯方法。此方法是基于粒子滤波的改进，在此方法中沉积层层数也作为未知数可以改变。观测数据可以选择为水听器阵列接收到的声压。首先，取多个快拍计算声压的交叉谱密度；然后，使用沉积层层数可变的粒子滤波反演海底地声参数与沉积层层数。该方法可以快速获取沉积层分层情况与海底地声参数，可以较为准确的估计出沉积层层数与层厚。此方法解决了传统粒子滤波中沉积层层数固定的问题，在粒子滤波中将沉积层层数考虑在内，给不同复杂度的沉积层模型不同的先验概率，从而在粒子滤波中实现沉积层层数的估计。本发明可以即时给出海底特性随距离变化的海域准确的沉积层分层情况，为实时准确反演局部海域的地声参数，获取沉积层分层结构提供了一种简便快速的方法。

为了实现上述目的，本发明提出了针对阵列接收声压数据的一种获取局部海域沉积层的分层结构和地声参数的方法，此发明利用的搜索算法是沉积层层数可变的粒子滤波，与传统粒子滤波不同，沉积层层数可变的粒子滤波在滤波过程中将沉积层层数也设为未知量；对于不同的沉积层层数，则前向模型的复杂度不同，层数越多则前向模型越复杂，根据贝叶斯定理，不同复杂度的模型有不同的先验概率；为了防止模型过拟合与欠拟合，更复杂的模型：正常模型：更简单的模型的先验概率应该为1∶n∶n

本发明提出了一种获取海底地声参数与沉积层分层结构的序贯方法，所述方法首先根据水听器接收到多个快拍的声压，计算所有帧的交叉谱密度作为观测数据；然后根据沉积层层数变化模型的复杂度按比例进行模型采样，并根据状态方程对地声参数进行高斯采样；基于采样的沉积层层数与地声参数计算声压拷贝场，并计算此声压拷贝场与声压观测数据的Bartlett相关，然后计算似然函数，再根据似然函数计算粒子权重；重采样逐帧反演不同距离处的沉积层层数与地声参数；结合所有帧的海底沉积层层数与地声参数信息，绘制出海底分层结构图。

作为上述技术方案的改进之一，获取海底地声参数与沉积层分层结构的序贯方法包括如下步骤：步骤1)给定一帧的距离为L，总帧数M，快拍数K，起始帧数m为1，初始化沉积层层数与地声参数，其中地声参数包括沉积层层厚、声速、密度、衰减；

步骤2)根据水听器接收到的多个快拍的声压p，计算第m帧的交叉谱密度

步骤3)根据步骤2)采样得到的沉积层层数与地声参数，计算声压拷贝场，并计算此声压拷贝场与声压实测数据的Bartlett相关，然后计算似然函数；

步骤4)根据似然函数，计算粒子权重，并重采样，令m＝m+1；

步骤5)重复步骤2)到步骤4)，逐帧反演不同距离处的沉积层层数与地声参数；

步骤6)结合所有帧的海底沉积层的层数、各层的深度厚度信息、以及地声参数信息，绘制出局部海域海底分层结构图。

作为上述技术方案的改进之一，在步骤1)中，具体步骤包括：给定一帧的距离为L，总帧数M，快拍数K，起始帧数m为1，初始化沉积层层数与地声参数，其中地声参数包括沉积层层厚、声速、密度、衰减α的具体步骤包括：

步骤1-1)将一定范围与距离相关的海底环境分为多个距离无关的海底环境，其中每一段距离无关的海底环境称为1帧，共M帧；

步骤1-2)对于第一帧，初始化沉积层层数m

作为上述技术方案的改进之一，在步骤2)中，具体步骤包括：根据多个快拍水听器接收到的声压p，计算此帧的交叉谱密度

步骤2-1)将较小的距离尺度L看作一帧，此帧中地声参数与沉积层分层结构可以近似看成不变的，即此帧具有相同的沉积层层数m、沉积层层厚h、声速c、密度ρ、衰减α；对此帧k个快拍得到的声压，计算

其中，f

步骤2-2)由于这里一帧的距离尺度较小，海底特性随距离变化通常不会特别迅速，所以这里沉积层层数变化只考虑层数增多一层、层数不变、层数减少一层这三种情况，这三种情况对应的模型复杂度依次减小，根据贝叶斯理论，对于均匀先验，这三种不同复杂度的模型应按照1∶n∶n

步骤2-3)层数增多时，设定最小层间距为0.2m，并设定沉积层层厚、声速、密度和衰减的建议分布方差，之后在满足最小层间距的情况下在沉积层最大深度之间随机插入一层，同时此新生层的地声参数以其两边任一层地声参数作为采样中心，根据建议分布方差进行高斯采样；层数减少时，设定沉积层层厚、声速、密度和衰减的建议分布方差，之后随机合并两层，同时，新层的地声参数以原来这两层中任一层为采样中心，根据建议分布方差进行高斯采样；层数不变时，此步骤中地声参数也不改变；

步骤2-4)粒子滤波状态方程

x

其中x

x

设置状态方程沉积层层厚、声速、密度和衰减的方差，对步骤2-3)中根据建议分布采样之后得到的地声参数，再次根据状态方程进行高斯采样；

作为上述技术方案的改进之一，在步骤3)中，具体步骤包括：根据步骤2采样得到的沉积层层数与地声参数，计算声压拷贝场，计算此声压拷贝场与声压实测数据的Bartlett相关，然后计算似然函数的具体步骤包括：

步骤3-1)根据步骤2采样得到的地声参数x

步骤3-2)根据

计算标准化的Bartlett不匹配度，其中

步骤3-3)使用观测方程计算声压观测值，观测方程为

y

其中v

计算似然函数L(x

其中

作为上述技术方案的改进之一，在步骤4)中，具体步骤包括：根据似然函数，计算粒子权重，重采样，m＝m+1的具体步骤包括：

步骤4-1)由于每一帧粒子都经过重采样过程，所以每一帧结束后粒子权重都相同，则粒子权重正比于似然函数的值；

步骤4-2)粒子权重越大，重采样后粒子数量越多，经过重采样后所有粒子权重相等；令m＝m+1，选取下一帧。

本发明的优点在于：本发明的方法可以自动提取局部海域的沉积层分层结构和地声参数，包括海底沉积层分层数、各沉积层厚度以及地声参数，同步提取，序贯反演沉积层分层情况及地声参数，反演速度快，可以有效应用到利用水听器声压数据进行局部海域综合海底环境测量及反演等领域。另外，本发明的方法将沉积层层数作为未知量，并给不同复杂度的模型不同的先验概率，在粒子滤波中实现了模型的选择，有效解决了粒子滤波无法反演沉积层层数的问题。本发明的方法可以应用于海底沉积层分层结构的快速获取，也可给出局部海域沉积层分层结构信息，为海底考古，声纳设备应用等提供直观、准确的海底结构信息。

附图说明

图1是本发明的一种获取海底地声参数与沉积层分层结构的序贯方法的流程图；

图2是本发明的一种获取海底地声参数与沉积层分层结构的序贯方法待处理的仿真环境；其中，图2(a)为声速随沉积层深度变化仿真图，图2(b)为海水密度随沉积层深度变化仿真图，图2(c)为衰减随沉积层深度变化仿真图；

图3是本发明的一种获取海底地声参数与沉积层分层结构的序贯方法对地声参数的敏感度分析；其中，图3(a)为不同频率下的Bartlett相关与水深的关系曲线图，图3(b)为不同频率下的Bartlett相关与沉积层厚度的关系曲线图，图3(c)为不同频率下的Bartlett相关与沉积层声速的关系曲线图，图3(d)为不同频率下的Bartlett相关与基底声速的关系曲线图，图3(e)为不同频率下的Bartlett相关与沉积层密度的关系曲线图，图3(f)为不同频率下的Bartlett相关与基底密度的关系曲线图，图3(g)为不同频率下的Bartlett相关与沉积层衰减的关系曲线图，图3(h)为不同频率下的Bartlett相关与基底衰减的关系曲线图；

图4(1)是本发明的一种获取海底地声参数与沉积层分层结构的序贯方法的反演结果图像；其中，图4(1-a)为海水密度随沉积层深度变化反演图，图4(1-b)为声速随沉积层深度变化反演图，图4(1-c)为衰减随沉积层深度变化反演图；图4(2)是本发明的一种获取海底地声参数与沉积层分层结构的序贯方法的层数反演结果图像；

图5是本发明的一种获取海底地声参数与沉积层分层结构的序贯方法的概率密度分布图像；其中，图5(a)为沉积层1声速概率密度分布图，图5(b)为沉积层2声速概率密度分布图，图5(c)为基底声速概率密度分布图，图5(d)为沉积层1海水密度概率密度分布图，图5(e)为沉积层2海水密度概率密度分布图，图5(f)为基底海水密度概率密度分布图，图5(g)为沉积层1衰减概率密度分布图，图5(h)为沉积层2衰减概率密度分布图，图5(i)为基底衰减概率密度分布图，图5(j)为沉积层分层概率密度分布图。

具体实施方式

以下结合附图对本发明作进一步的详细说明。

本发明的方法给出交叉谱数据处理的技术流程，和具体的反演与距离相关海域沉积层分层情况与地声参数的实施步骤。本发明的方法在利用粒子滤波反演地声参数过程中，基于贝叶斯原理，增加了对沉积层层数的选择。此方法避免了层数反演错误的情况，获得更准确的沉积层分层情况；而且沉积层层数的正确反演有利于更准确的地声参数反演。

本发明提供了一种在海底沉积层层数未知时获取沉积层分层结构与地声参数的方法，包括：首先根据水听器接收到多个快拍的声压，计算所有帧的交叉谱密度作为观测数据；然后根据沉积层层数变化模型的复杂度按比例进行模型采样，并根据状态方程对地声参数进行高斯采样；基于采样的沉积层层数与地声参数计算声压拷贝场，并计算此声压拷贝场与声压观测数据的Bartlett相关，然后计算似然函数，再根据似然函数计算粒子权重；然后对粒子重采样，逐帧反演不同距离处的沉积层层数与地声参数；结合所有帧的海底沉积层层数与地声参数信息，绘制出海底分层结构图。

如图1所示，本发明提供了一种获取海底地声参数与沉积层分层结构的序贯方法，其具体步骤包括：

步骤1)给定一帧的距离为L＝180m，总帧数M＝13，快拍数k＝10，起始帧数m为1，初始化沉积层层数m

步骤2)根据水听器接收到的多个快拍的声压p，计算第一帧的交叉谱密度

步骤3)根据步骤2)采样得到的沉积层层数与地声参数，计算声压拷贝场，并计算此声压拷贝场与声压实测数据的Bartlett相关，然后计算似然函数；

步骤4)根据似然函数，计算粒子权重，并重采样，令m＝m+1；

步骤5)重复步骤2)到步骤4)，逐帧反演不同距离处的沉积层层数与地声参数，地声参数声速密度衰减的反演结果如图4(1-a)～(1-c)所示，沉积层层数的反演结果如图4(2)所示，各参数的后验概率密度如图5(a)～(j)所示，反演结果显示沉积层层数在第4帧变多，说明此方法能够较好跟踪沉积层层数的变化，各层地声参数均分布在真值附近；

步骤6)结合所有帧的海底沉积层的层数、各层的深度层厚信息、以及地声参数信息，绘制出局部海域海底分层结构图如图4(1-a)～(1-c)所示；

在步骤1)中，给定一帧的距离为L＝180m，总帧数M＝13，快拍数k＝10，起始帧数m为1，初始化沉积层层数m

步骤1-1)将一定范围与距离相关的海底环境分为多个距离无关的海底环境，其中每一段距离无关的海底环境称为1帧，共M＝13帧，水听器阵列为位于距船尾150m处，阵列为10元阵，阵元间距为等间隔的3m，声源深度设置为7m，阵列深度设置为17m，水声速剖面设置为0-15m等声速1530m/s，15-20m声速从1530m/s均匀变为1505m/s，20m到海底声速为1505m/s，海水密度为1g/cm

步骤1-2)对于第一帧，初始化水深h

在步骤2)中，根据多个快拍水听器接收到的声压p，计算第一帧的交叉谱密度

步骤2-1)将较小的距离尺度L＝180m看作一帧，此帧中地声参数与沉积层分层结构可以近似看成不变的，即此帧范围内具有相同的沉积层层数m、沉积层层厚h、声速c、密度ρ、衰减α；对此帧k个快拍得到的声压，计算

其中，f

步骤2-2)由于这里一帧的距离尺度较小，海底特性随距离变化通常不会特别迅速，所以这里沉积层层数变化只考虑层数增多一层、层数不变、层数减少一层这三种情况，这三种情况对应的模型复杂度依次减小，根据贝叶斯理论，对于均匀先验，这三种不同复杂度的模型应按照1∶n∶n

步骤2-3)层数增多时，设定最小层间距为0.2m，并设定各参数建议分布方差分别为

步骤2-4)粒子滤波状态方程

x

其中x

x

将状态方程方差分别设置为σ

在步骤3)中，根据步骤2采样得到的沉积层层数与地声参数，计算声压拷贝场，计算此声压拷贝场与声压实测数据的Bartlett相关，然后计算似然函数的具体步骤包括：

步骤3-1)根据步骤2采样得到的海底信息x

步骤3-2)根据

计算标准化的Bartlett不匹配度，其中

步骤3-3)观测方程为

y

其中v

计算似然函数，其中N＝10为传声器阵列阵元个数，J＝5为频点个数，

其中

在步骤4)中，根据似然函数，计算粒子权重，重采样，m＝m+1的具体步骤包括：

步骤4-1)由于每一帧粒子都经过重采样过程，所以每一帧结束后粒子权重都相同，则粒子权重正比于似然函数的值；

步骤4-2)粒子权重越大，重采样后粒子数量越多，经过重采样后所有粒子权重相等。

从上述对本发明的具体描述可以看出，本发明可以即时给出海底特性随距离变化的海域准确的沉积层分层情况，为实时准确反演局部海域的地声参数，获取沉积层分层结构。

最后应说明的是，以上实施例仅用以说明本发明的技术方案而非限制。尽管参照实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，对本发明的技术方案进行修改或者等同替换，都不脱离本发明技术方案的精神和范围，其均应涵盖在本发明的权利要求范围当中。