根据非线性贝叶斯理论的浅海多层海底地声参数反演方法

技术领域

本发明属于地声参数反演技术领域，具体涉及一种根据非线性贝叶斯理论的浅海多层海底地声参数反演方法。

背景技术

海底地声参数是组成海洋水声环境的重要参数之一，海底的声速、密度、声速衰减等声学参数对海洋环境中，尤其是浅海环境中的声传播有着重要影响。对上述海底地声参数的掌握程度，将直接影响对水声设备性能的预测与评估、海洋声场的数值预报、海洋声场特征的利用等。如何高效、准确获取海底地声参数信息一直是水声领域的研究热点。

对海底地声参数的获取方法目前分为直接测量与间接测量两类。相比通过钻孔取样等方式获取海底底质样品进行鉴别的直接测量方法，以声学反演技术为代表的地声参数间接测量方法因其具有实时、快速、高效的技术优势，而被广泛应用于海底地声参数的获取。由于常规声呐多依托中/高频段声波开展工作，因而在以往对海底地声参数的反演研究中多只关注海底表层声学特性，且假设海底为液态介质。随着近年来声呐设备向低频/甚低频的发展，以往只对海底表层声学特性的了解和掌握已不能满足对当前声传播问题的解析和验证，开展包含海底结构在内的、深层海底地声参数反演技术的研究愈发迫切。且已有研究成果已经证明，在研究低频/甚低频水声传播问题时，海底横波声速的影响不可忽略，因而在当前海底地声参数反演问题的研究中，将海底视为分层弹性介质，对包含分层结构、横波声速及其衰减在内、深层地声参数进行准确反演是当前海底地声参数的发展目标，相关研究工作亟待开展。

此外，考虑海底地声参数反演核心是对一个非线性、多参数且参数间存在相关性的复杂函数的求解，在以往研究中多用遗传算法、模拟退火算法等[6]全局寻优方法进行求解。但上述算法寻优结果仅是一组能反映声场数据信息“最佳”参数。由于海上实测数据不可避免的存在各类误差，这些因素将导致反演结果存在一定的不确定性。

非线性贝叶斯反演方法是基于概率论、应用数学和最优化理论等交叉学科的全局优化算法，其不仅可以有效估计最大后验概率(MAP)模型参数，而且可以从统计的角度对参数反演结果的不确定性做出解释，更符合当下地声参数反演的研究需要。

因此，基于上述技术问题需要设计一种新的根据非线性贝叶斯理论的浅海多层海底地声参数反演方法。

发明内容

本发明的目的是提供一种根据非线性贝叶斯理论的浅海多层海底地声参数反演方法。

为了解决上述技术问题，本发明提供了一种浅海多层海底地声参数反演方法，包括：

根据非线性贝叶斯反演方法对海底地声参数进行反演。

进一步，所述根据非线性贝叶斯反演方法对海底地声参数进行反演的方法包括：

构建海底模型，并获取海底模型的海底地声参数；

获取缩比实验中的实验数据；

缩比实验中的实验数据对应的向量d和海底模型参数对应的向量m满足贝叶斯定理：

P(m|d)＝P(d|m)P(m)/P(d)；

其中，P(m|d)为后验概率密度；d的条件概率P(d|m)用似然函数L(m)来表示；P(m)是向量m的先验概率密度函数；P(d)是向量d的概率密度函数；

P(d)与m无关，则

P(m|d)∝P(d|m)P(m)；

似然函数L(m)为：

L(m)＝P(d|m)∝exp[-E(m)]；

其中，E(m)为误差函数；

将E(m)进行归一化：

其中，积分域跨越M维参数空间，M为待反演参数的个数；m′为表示任意模型参数向量m的未知数。

进一步，所述根据非线性贝叶斯反演方法对海底地声参数进行反演的方法还包括：

获取待反演参数的MAP值：

获取待反演参数的均值：

获取待反演参数的一维概率密度分布：P(m

其中，δ狄拉克Delta函数；

获取待反演参数的相互关系：

其中，相关值R

进一步，所述根据非线性贝叶斯反演方法对海底地声参数进行反演的方法还包括：

当数据误差为独立同分布的随机变量时似然函数为：

其中，p

预测声压p

其中，p

当

其中，*为共轭转置；

将角线协方差处理为C

其中，方差v

当

根据方差v

进一步，所述根据非线性贝叶斯反演方法对海底地声参数进行反演的方法还包括：

当测量数据d时，相应模型的似然函数为：

其中，M是相应模型中参数的个数，N为数据参数个数；

用误差函数代替似然函数以获取：

其中，BIC为贝叶斯信息准则，BIC值最小的模型即为最优模型。

进一步，所述根据非线性贝叶斯反演方法对海底地声参数进行反演的方法还包括：

基于误差函数在先验区间搜索MAP值，并基于声阻抗的物理规律进行微正向扰动，即

初始化参数m

计算误差函数；

当满足m

判断随机扰动计算的结果是否满足Metropolis准则，满足时保留参数m

判断保留的参数m

m

所述声阻抗的物理规律为：c

其中，c表示声速，ρ表示密度。

本发明的有益效果是，本发明通过根据非线性贝叶斯反演方法对海底地声参数进行反演，实现了从实测数据中准确判定实际海底结构；采用拨正模拟退火算法可提高反演结果的准确度及辨析度；与两类声速衰减及密度相比，纵波声速和横波声速的不确定性更小；相比半无限基底，沉积层中各参数的不确定性更小。

本发明的其他特征和优点将在随后的说明书中阐述，并且，部分地从说明书中变得显而易见，或者通过实施本发明而了解。本发明的目的和其他优点在说明书以及附图中所特别指出的结构来实现和获得。

为使本发明的上述目的、特征和优点能更明显易懂，下文特举较佳实施例，并配合所附附图，作详细说明如下。

附图说明

为了更清楚地说明本发明具体实施方式或现有技术中的技术方案，下面将对具体实施方式或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图是本发明的一些实施方式，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1是本发明所涉及的浅海多层海底地声参数反演方法的流程图；

图2是本发明所涉及的海底模型的示意图；

图3是本发明所涉及的实测数据与反演数据对比示意图；

图4是本发明所涉及的模型选择结果示意图；

图5是本发明所涉及的纵波声速、横波声速以及海底密度随深度变化刨面结构图；

图6是本发明所涉及的沉积层海底参数一维边缘概率分布示意图；

图7是本发明所涉及的半无限海底参数一维边缘概率分布示意图；

图8是本发明所涉及的反演参数间相关矩阵图；

图9是本发明所涉及的纵波声速、横波声速、密度随深度变化的二维边缘概率分布示意图。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合附图对本发明的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

如图1所示，本实施例提供了一种浅海多层海底地声参数反演方法，包括：根据非线性贝叶斯反演方法对海底地声参数进行反演；应用声压场为研究对象，进而采用贝叶斯方法估算海底各层密度、纵波声速、横波声速、两类声速衰减等参数及其不确定性；各参数的最大后验概率密度(MAP)估计值和边缘概率分布分别通过拨正模拟退火算法和Metropolis-Hastings采样法在各参数先验区间内搜索获得，并依据贝叶斯信息准则(BIC)从不同参数化模型中选择最优模型；可从实测数据中准确判定实际海底结构；采用拨正模拟退火算法可提高反演结果的准确度及辨析度；与两类声速衰减及密度相比，纵波声速和横波声速的不确定性更小；相比半无限基底，沉积层中各参数的不确定性更小。

在本实施例中，所述根据非线性贝叶斯反演方法对海底地声参数进行反演的方法包括：构建海底模型(如图2所示，其中c

P(m|d)＝P(d|m)P(m)/P(d)；

其中，P(m|d)为后验概率密度(PPD)；d的条件概率P(d|m)通常用似然函数L(m)来表示；P(m)是向量m的先验概率密度函数，表示独立于数据的可用模型参数先验信息；P(d)是向量d的概率密度函数；

P(d)与m无关可以看作一个常数，则

P(m|d)∝P(d|m)P(m)；

似然函数由数据形式和数据误差的统计分布决定。考虑到在实际应用过程中，误差的统计特征很难独立获得，在处理过程中采用无偏高斯误差的假设，似然函数L(m)为：

L(m)＝P(d|m)∝exp[-E(m)]；

其中，E(m)为误差函数；

将E(m)进行归一化：

其中，积分域跨越M维参数空间，M为待反演参数的个数；m′表示可以表示任意模型参数向量m的未知数，m和m′实际上相同，只是在积分项中习惯性写法；在贝叶斯理论中，后验概率密度(PPD)可作为反演问题的解。

在本实施例中，所述根据非线性贝叶斯反演方法对海底地声参数进行反演的方法还包括：由于反演中存在对维参数问题，为更合理的解释参数反演结果还需对模型参数间的相关特性进行研究；

获取待反演参数的MAP值：

获取待反演参数的均值：

获取待反演参数的一维概率密度分布：P(m

其中，δ狄拉克Delta函数；

待反演参数间的相互关系可通过标准化协方差产生的相关矩阵来定量描述：

其中，相关值R

在本实施例中，所述根据非线性贝叶斯反演方法对海底地声参数进行反演的方法还包括：选取合适的先验参数范围是反演成功与否的关键，为保证采样的完全，本文先验信息考虑每个参数在其足够宽范围内均匀分布；MAP估计值采用拨正模拟退火算法(Correction Simulated Annealing，CSA)求解，CSA法是一种改进的模拟退火算法，可依据上下层海底物理特性在寻优求解过程中对随机参数进行拨正；参数的一维、二维边缘概率分布通过Metropolis-Hasting采样法(MHS)得到的参数PPD的数值积分求解，虽然MHS采样法的采样结果不依赖初始模型，但是如果选用了正确的初始值，可提高反演的精度和速度，本文采用CSA法搜索到的MAP估计值作为MHS采样的初始模型；

在贝叶斯反演理论中，求解参数PPD需要获得似然函数L(m)，似然函数与数据误差(包括测量误差和理论误差)的统计分布有关，是定量描述参数不确定性的重要指标；本实施例中假设数据误差是独立同分布的随机变量，则似然函数可以表示为：

其中，p

预测声压p

其中，p

当

其中，*为共轭转置；

忽略数据的空间相关性将角线协方差处理为C

其中，方差v

当

根据方差v

获取满足最大似然函数估计值时对应的误差函数E(m)：

在本实施例中，所述根据非线性贝叶斯反演方法对海底地声参数进行反演的方法还包括：合理参数化模型是贝叶斯反演的关键，欠参数化模型使得结构无法完整解析，导致模型的不确定性偏低；过参数化模型对参数的约束不够，导致模型的不确定性增加，欠参数化和过参数化模型都会对反演结果造成一定的影响。本文应用贝叶斯信息准则(BIC)选择与实测数据最符合的参数化模型。BIC值是从多维变量的正态分布中得到的，并不是一个精确值，是模型I的贝叶斯定理P(d|I)的渐进近似，也就是假定测量数据d(缩比实验中的实验数据)，模型I(模型I指的是不同的反演模型，即海底模型，BIC就是从不同的反演模型中找出最符合实验数据的反演模型)的似然函数，其表达式为：

其中，M是模型I(相应模型)中参数(海底模型的参数，海底分层的时候，不同的海底分层，海底模型中包含的反演参数是不桐的，比如一层海底反演参数M＝5,二层海底模型M＝11)的个数，N为数据参数(缩比实验中的实验数据)个数；

用误差函数代替似然函数以获取：

其中，BIC为贝叶斯信息准则，BIC值最小的模型(BIC是从不同的模型I中选择出最符合实验或者仿真参数的模型，这个模型是从模型I中选择出的一个模型)即为最优模型，BIC数值的大小由误差函数、模型参数个数和数据个数共同决定，从而避免了欠参数化和过参数化模型，更有效的选择出最优参数化模型。

在本实施例中，所述根据非线性贝叶斯反演方法对海底地声参数进行反演的方法还包括：基于误差函数在先验区间搜索MAP值，是贝叶斯反演方法的重要一步，由于误差函数是非线性，多维度且可能存在多个峰值，给搜索过程带了极大的挑战，且参数间的相关性加剧了搜索的困难度；目前，绝大多数学者运用局部寻优或全局寻优算法应用到地声参数反演研究中；局部寻优方法通过向下移动局部梯度来迭代改善初始模型，这些方法可以有效下坡移动并适应相关参数，但容易陷入局部最小；全局寻优相对于局部寻优算法，通过适时向上移动的方式避免了陷入局部最小，可定向随机的搜索整个参数空间；

在进行MAP值搜索时，参数区间的选择是尤为重要的，而参数区间选择随着参数化模型逐步复杂而变得困难；以往针对半无限海底反演时，在误差函数与寻优算法选择合理的条件下，逐步寻优可得到MAP值，但随着海底模型层数的增加，沉积层厚度的加入使得寻优过程变得更为复杂，参数空间的选择更为困难，参数空间范围太大增加寻优难度以及准确度，范围太小不能保证实测数据真值包含在寻优区间内，因此本文基于模拟退火算法(SA)进行改进提出了拨正模拟退火算法(CSA)，并应用到贝叶斯反演方法中：

初始化参数m

计算误差函数；

当满足m

判断随机扰动计算的结果是否满足Metropolis准则(Metropolis接受准则)，满足时保留参数m

判断保留的参数m

m

在水声学中声阻抗随着海底深度的增加而增大，根据这一海底物理规律利用上一层寻优结果拨正下一层寻优，所述声阻抗的物理规律为：c

其中，c表示声速，ρ表示密度；根据声阻抗的物理规律，第二层的声阻抗大于第一层的声阻抗，所以在寻优时加入微正向扰动，使得反演过程更加符合物理规律。

在本实施例中，采用尺寸为1.5m×1.1m×0.105m(长×宽×高，测量误差±1mm)的PVC板(Polyvinyl Chloride Polymer，聚氯乙烯，测得密度为(1.20g·cm-3)模拟弹性海底，通过在PVC板上铺放细沙的方法模拟具有弹性沉积层与弹性半无限海底的浅海波导环境，声源深度z

表1：反演结果

模型研究考虑3组不同的模型，分别为1～3层各向同性的均匀海底分布模型，并给出了每种模型经CSA求解的MAP值对应的数据误差值、参数个数和BIC值。

从图4中(a)中可以看出(m1、m2、m3、分别表示1层、2层、3层均匀海底分布模型，即海底模型)，随着讨论模型(海底模型)的层数不断增加其误差值逐渐降低，整体呈现下降的趋势，这与实际物理情况是相符的；从图4中(b)中可以看出随着层数的增加，所需考虑的参数也在增加，这会增加反演结果的不确定性；BIC充分考虑了模型的误差值、参数个数以及数据个数之间的平衡，选择出与实测数据最符合的参数化模型。通过图4中(c)可以看出，通过BIC准则最终确定的最优模型为2均匀海底分布模型。

图5给出了反演得到的对应的纵波声速、横波声速以及密度随深度变化的剖面结构(m1、m2和m3表示1层、2层、3层均匀海底分布反演模型，即海底模型)。通常海底声阻抗(ρc，其中，ρ为密度c为声速)与深度为正比例关系，即声阻抗随深度的增大而增大，从结果中可以看出1，2和3层海底分布模型趋势均符合这一物理规律，随着海底分层的增加，参数的分布并未出现“折回”现象，进而验证了拨正模拟退火算法在参数先验区间中搜寻参数时的关键作用。从剖面结构中可以看出，以1000mm为分界线，在该值以内纵波声速、横波声速以及密度有着较为明显的结构划分，以外均趋向于某一值，其中纵波声速趋势更为明显。

非线性贝叶斯反演方法可对反演结果的不确定性进行分析，由于地声参数反演是一个非线性、多参数优化的复杂问题，只是通过一组最优解作为反演结果显然是不合理的，因此，反演结果的不确定分析尤为重要。

选择BIC准则确定的3层均匀分布的最优模型，利用CSA法搜索得到的估计值作为采样的初始模型，并通过Metropolis-Hasting采样法得到参数的PPD数值积分求解，图6、图7给出了2层均匀海底分布模型对应的参数一维概率分布图，其中红色实线表示各参数的MAP值；从图中可以看出，沉积层以及半无限海底中各参数的一维边缘概率分布在先验区间内均符合正态分布，相比半无限海底沉积层各参数分布较窄，表明沉积层参数更敏感；从单个参数来看，c

考虑到地声参数反演是多参数优化问题，参数间是相互耦合的，参数间的相关性会直接影响反演结果，增加待反演参数的不确定性，为此本文分析了反演参数间的相关性，如图8所示。从结果中可以看出，α

为更直观、准确的表示纵波声速、横波声速和密度随深度的变化关系，图9给出了纵波声速、横波声速以及密度与深度的二维边缘概率分布图，由于反演参数为沉积层及下层海底参数，即二维边缘概率密度分布从沉积层开始。从结果中可以看出：相对于半无限海底，沉积层的纵波声速、横波声速以及密度的概率分布更窄，不确定性更小，进一步表明沉积层各参数更敏感。

利用声压场进行贝叶斯反演可以获取海底参数，通过BIC准则确定的最优模型为3层均匀分布模型，并且与实验结果相符合；针对多层海底模型的非线性贝叶斯反演方法，本文所提出的拨正模拟退火算法可高效准确的获取MAP估计值，避免了在多层海底模型中出现的参数“折回”现象。通过对反演结果的不确定分析可以得出，在进行声压场反演时，与纵波衰减、横波衰减以及密度相比，纵波声速和横波声速的不确定性更小，其敏感性更强；相比半无限海底沉积层各参数更敏感，不确定性更小。

综上所述，本发明通过根据非线性贝叶斯反演方法对海底地声参数进行反演，实现了从实测数据中准确判定实际海底结构；采用拨正模拟退火算法可提高反演结果的准确度及辨析度；与两类声速衰减及密度相比，纵波声速和横波声速的不确定性更小；相比半无限基底，沉积层中各参数的不确定性更小。

在本申请所提供的几个实施例中，应该理解到，所揭露的装置和方法，也可以通过其它的方式实现。以上所描述的装置实施例仅仅是示意性的，例如，附图中的流程图和框图显示了根据本发明的多个实施例的装置、方法和计算机程序产品的可能实现的体系架构、功能和操作。在这点上，流程图或框图中的每个方框可以代表一个模块、程序段或代码的一部分，所述模块、程序段或代码的一部分包含一个或多个用于实现规定的逻辑功能的可执行指令。也应当注意，在有些作为替换的实现方式中，方框中所标注的功能也可以以不同于附图中所标注的顺序发生。例如，两个连续的方框实际上可以基本并行地执行，它们有时也可以按相反的顺序执行，这依所涉及的功能而定。也要注意的是，框图和/或流程图中的每个方框、以及框图和/或流程图中的方框的组合，可以用执行规定的功能或动作的专用的基于硬件的系统来实现，或者可以用专用硬件与计算机指令的组合来实现。

另外，在本发明各个实施例中的各功能模块可以集成在一起形成一个独立的部分，也可以是各个模块单独存在，也可以两个或两个以上模块集成形成一个独立的部分。

所述功能如果以软件功能模块的形式实现并作为独立的产品销售或使用时，可以存储在一个计算机可读取存储介质中。基于这样的理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分或者该技术方案的部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质中，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本发明各个实施例所述方法的全部或部分步骤。而前述的存储介质包括：U盘、移动硬盘、只读存储器(ROM，Read-Only Memory)、随机存取存储器(RAM，Random Access Memory)、磁碟或者光盘等各种可以存储程序代码的介质。

以上述依据本发明的理想实施例为启示，通过上述的说明内容，相关工作人员完全可以在不偏离本项发明技术思想的范围内，进行多样的变更以及修改。本项发明的技术性范围并不局限于说明书上的内容，必须要根据权利要求范围来确定其技术性范围。