一种混合全局优化算法的叠前地震多参数反演方法与流程

本发明涉及一种混合全局优化算法的叠前地震多参数反演方法，具体涉及一种结合多维学习粒子群算法和快速模拟退火算法的混合全局优化算法，属于油气地球物理勘探中的地震资料反演技术领域。

背景技术：

叠前地震反演是一种基于褶积模型的反演，它基于振幅随炮检距变化理论，可以同步获得纵波速度、横波速度和密度等多种弹性参数，可以有效识别储层流体类型和地质构造特征。因此，叠前地震反演是目前应用最广泛、发展最成熟的地震反演技术之一，在油气资源的勘探和开发过程中发挥重要作用。

地震反演是一个非线性优化问题，即使得估算的地层参数所对应的合成地震记录与观测地震记录间的误差达到最小。目前，地震反演的最优化算法可以分为两大类：第一类是局部梯度类算法，如最速下降法、共轭梯度法和牛顿类方法等；第二类是全局优化算法，如模拟退火算法、遗传算法和粒子群算法等。梯度类算法计算效率较高，对非线性反演问题可以快速收敛获得最优解，但此类算法对初始模型依赖较高，并且对复杂的非线性反演问题效果不理想。而全局最优化方法可以解决非线性和多峰值反演问题，并且不依赖于初始模型，特别适用于叠前地震多参数反演这类复杂的非线性反演问题。

模拟退火算法和粒子群算法是目前发展较成熟的全局优化算法，已有不少应用于叠前地震反演的成功案例。但模拟退火算法计算量较大，收敛速度较慢，虽然粒子群算法运算效率有显著提高，但易不成熟收敛而陷入局部最小值。因此，结合粒子群算法和模拟退火算法的混合全局最优化算法对两类算法进行优势互补，使得粒子群算法具有模拟退火算法的“跳跃”机制，扩大了粒子搜索和扩展的能力，并且不易陷落局部最优解；此外，叠前地震反演是一种多参数的非线性反演问题，多参数反演结果具有不稳定性，而常规反演算法本身并未有针对该问题的解决方法，考虑到多参数间存在基于地质统计学的约束关系，在粒子群算法中添加多维学习项，可以有效解决叠前地震反演多参数的不稳定问题。

技术实现要素：

本发明所要解决的技术问题是：提供一种混合全局优化算法的叠前地震多参数反演方法，克服叠前地震多参数反演的不稳定问题，并解决传统粒子群算法收敛不成熟的问题，可以同步而准确的获取地震三参数的反演结果。

本发明为解决上述技术问题采用以下技术方案：

一种混合全局优化算法的叠前地震多参数反演方法，包括如下步骤：

步骤1，确定粒子待反演参数的取值范围，设定粒子的数量、最大迭代次数、模拟退火冷却进度表、模拟退火学习项加权系数和多维学习项加权系数；

步骤2，设置各个粒子的初始位置和初始速度，位置包括四个维度：纵波速度、横波速度、密度、前述三个参数的联合概率密度；

步骤3，对第k次迭代，根据粒子的位置计算粒子适应度，其中粒子适应度由观测地震记录与合成地震记录的误差和步骤2所述三个参数的先验约束项构成；

步骤4，根据步骤3的粒子适应度，计算各粒子位置被选择作为模拟退火学习项引导粒子的接收概率，并基于各粒子的接收概率，通过轮盘选择法则确定模拟退火学习项引导粒子；

步骤5，对每个粒子，固定该粒子的横波速度和密度，用所有粒子的纵波速度依次替换该粒子的纵波速度，并利用联合概率密度公式计算所有粒子的纵波速度对应的联合概率密度，找出最大联合概率密度；从每个粒子对应的最大联合概率密度中找出最大值及该最大值对应的纵波速度vpmax；按上述同样的方法，找到横波速度vsmax和密度ρmax；根据vpmax、vsmax和ρmax，得到多维学习项引导粒子；

步骤6，根据模拟退火学习项引导粒子和多维学习项引导粒子更新各粒子的速度和位置；

步骤7，进入第k+1次迭代，重复步骤3至步骤6，直至达到最大迭代次数，且到达模拟退火冷却进度表中的终止退火温度，输出粒子适应度最优的粒子。

作为本发明的一种优选方案，步骤1所述待反演参数包括纵波速度、横波速度和密度。

作为本发明的一种优选方案，步骤1所述模拟退火冷却进度表包括初始退火温度、终止退火温度，温度从初始退火温度开始逐渐降低直至终止退火温度，所有退火温度的个数与最大迭代次数相同，且一次迭代对应一个退火温度。

作为本发明的一种优选方案，步骤2所述联合概率密度表达式为：

其中，fi为第i个粒子纵波速度vpi、横波速度vsi和密度ρi这三个参数的联合概率密度，σ为纵波速度vpi、横波速度vsi和密度ρi这三个参数的协方差矩阵，pi为第i个粒子位置，e(pi)为pi的期望，上标t表示矩阵的转置。

作为本发明的一种优选方案，步骤3所述粒子适应度表达式为：

其中，为第i个粒子适应度，pi为第i个粒子位置，e(pi)为pi的期望，l为观测数据采样长度，θ为入射角，t为地震记录采样时间，w为震源子波，r为纵波反射系数，d为观测地震记录，σ为纵波速度vpi、横波速度vsi和密度ρi这三个参数的协方差矩阵，λ1、λ2均为预设定系数，上标t表示矩阵的转置。

作为本发明的一种优选方案，步骤4所述接收概率表达式为：

其中，piaccept为第i个粒子的接收概率，为第i个粒子适应度，n为粒子的数量，tk为第k次迭代的退火温度。

作为本发明的一种优选方案，所述步骤6更新公式为：

粒子速度更新公式：

粒子位置更新公式：

其中，d表示粒子的维度，vi为第i个粒子的速度，k、k-1分别为第k、k-1次迭代，c1、c2分别为模拟退火学习项加权系数、多维学习项加权系数，rand为[0,1]的随机函数，pmin为模拟退火学习项引导粒子，pmax为多维学习项引导粒子，pi为第i个粒子位置。

本发明采用以上技术方案与现有技术相比，具有以下技术效果：

1、本发明方法将粒子群算法和快速模拟退火算法有效结合，解决传统粒子群算法易不成熟收敛的问题，并在粒子群算法中添加基于三参数联合概率密度择优组合的多维学习项，克服叠前地震多参数同步反演的不稳定，可以同步而准确的获取纵波速度、横波速度和密度三参数反演结果。

2、本发明方法与传统粒子群算法的反演结果相比反演效果改善明显。

附图说明

图1是本发明混合全局优化算法的叠前地震多参数反演方法的流程图。

图2是本发明各反演参数的理论模型，其中，(a)为纵波速度，(b)为横波速度，(c)为密度。

图3是理论模型的合成地震数据(观测地震数据)。

图4是使用本发明方法的各反演参数反演结果与理论模型的对比，其中，(a)为纵波速度，(b)为横波速度，(c)为密度。

图5是使用传统粒子群算法的各反演参数反演结果与理论模型的对比，其中，(a)为纵波速度，(b)为横波速度，(c)为密度。

图6是使用本发明方法和传统粒子群算法进行反演的收敛速度的对比。

具体实施方式

下面详细描述本发明的实施方式，所述实施方式的示例在附图中示出。下面通过参考附图描述的实施方式是示例性的，仅用于解释本发明，而不能解释为对本发明的限制。

如图1所示，为本发明反演方法的流程图，具体步骤如下：

步骤一，初始化设置。待反演模型离散化，设定各位置待反演的纵波速度、横波速度和密度的取值范围；设定粒子的数量n、最大迭代次数k和模拟退火冷却进度表；设置学习项加权系数c1和c2。

加权系数c1和c2分别对应粒子速度更新的模拟退火学习项和多维学习项。对这两种学习项的选择通过加权系数c1和c2进行平衡。模拟退火学习项结合快速模拟退火算法，模拟退火算法的“跳跃”机制使得这种学习有一定的概率接受非最优位置的粒子作为引导粒子，从而扩大粒子的更新范围和搜索能力，不易陷入局部最优解；多维学习项是粒子向其他维度粒子位置的学习，多维学习项的引导粒子是基于三参数间的统计学关系而选择的具有最大概率密度值的组合粒子，可以有效控制三参数同步反演的稳定性。

模拟退火冷却进度表包括初始温度t1和终止温度tend，设定温度个数与最大迭代次数相同，即每个迭代k与一个温度tk对应，随着迭代进行，温度逐渐降低直至终止温度。温度控制模拟退火学习项中各粒子位置被选择作为引导粒子的接收概率。

步骤二，设置粒子初始位置和初始速度。每个粒子位置pi包含四个维度(即纵波速度、横波速度和密度三参数，以及三参数的联合概率密度)，即pi＝[pi¹,pi²,pi³,pi⁴]＝[vpi,vsi,ρi,fi],i＝1,2,…,n。

其中，三参数的联合概率密度的具体表达式为：

其中，e为期望；σ为纵波速度、横波速度和密度三参数的协方差矩阵，由测井数据获取。

在取值范围内，随机设置粒子的初始位置和初始速度。粒子的初始位置组成[4×n]的矩阵，粒子的初始速度设为零。

步骤三，设置目标函数，即粒子适应度表达式。本反演方法的目标函数由观测地震记录与合成地震记录的误差和三参数的先验约束项构成，具体表现形式如下：

其中，f为粒子适应度；r为使用精确佐普里兹方程计算的纵波反射系数；w为震源子波；d为观测地震记录；l为观测数据采样长度；λ1和λ2为预设定系数；θ为入射角；t为地震记录采样时间。

步骤四，设置模拟退火学习项。按公式(2)计算各粒子pi的适应度值并计算各粒子位置被选择作为模拟退火学习项的引导粒子pmin的接收概率，接收概率的表达式为：

其中，piaccpet为位置pi被接收的概率，tk为第k次迭代的退火温度；为最优粒子的适应度值。

基于各粒子的接收概率，使用轮盘选择算法决定模拟退火学习项的引导粒子pmin，显然，适应度最小的粒子具有最大的接收概率，但非最优粒子也有一定概率“跳跃”成为引导粒子，且这种概率在高温时最大，并随着温度降低而减小。这种“跳跃”机制增加了粒子的搜索扩展能力，避免因不成熟收敛而陷入局部最优解。

步骤五，设置多维学习项。

固定横波速度和密度，由公式(1)计算各粒子的纵波速度对应的联合概率密度值，取最大概率密度值对应的纵波速度vpmax；按同样方法，分别获得最大概率密度对应的横波速度vsmax和密度ρmax，即：

多维学习项的引导粒子pmax是具有最大概率密度值的择优组合粒子，即

pmax＝[vpmax,vsmax,ρmax,fmax]

多维学习项的引导粒子是基于三参数间的统计学关系而选择的具有最大概率密度值的组合粒子，该引导粒子不局限于单一粒子所对应的位置，而是多维度、多粒子间择优组合的结果，可以有效控制叠前地震三参数同步反演的稳定性。

步骤六，更新粒子速度和粒子位置。

粒子速度的更新公式为：

其中，为模拟退火学习项，其中的引导粒子pmin由步骤四得到；为多维学习项，其中的引导粒子pmax由步骤五得到，上标d代表粒子的维度，每个粒子pi包含四个维度(即纵波速度、横波速度、密度三参数，以及三参数的联合概率密度)。

粒子位置的更新公式为：

步骤七，进入下一次迭代k+1，并使温度降低至tk+1，重复步骤三至步骤六。直至达到最大迭代数k，温度降至终止温度tend，输出适应度最优的粒子。

下面以一个合成地震数据测试进行具体说明：

合成地震数据测试的理论模型如图2所示，其中，(a)为纵波速度，(b)为横波速度，(c)为密度。使用精确佐普里兹方程计算理论模型中地层的纵波反射系数，计算的角度间隔为5°，角度覆盖范围为0-50°。分别将每层的11组反射系数与主频为50hz的零相位理论ricker子波进行褶积，得到如图3所示的合成地震记录(角度道集)，即观测地震数据，其中有11个角度道，时间采样间隔2ms。

使用本发明方法对观测地震数据进行反演，以获得纵波速度、横波速度和密度三参数的反演结果，使反演结果的合成地震数据与观测地震数据之间误差的二次范数极小，具体实现方式如下：

设置粒子位置的取值范围，每个粒子位置包含四个维度，即纵波速度、横波速度、密度和联合概率密度，对应的最小值为[1.2,0.6,1.2,0]和最大值为[5.6,3.5,3.2,1]，粒子数量n为20，最大迭代次数k为300，初始温度t1为0.9，终止温度tend为0.003，学习项加权系数c1＝c2＝1。

在取值范围内，随机生成粒子的初始位置，并对粒子位置进行对数归一化处理；设定各粒子的初始速度为零。

按步骤三计算各粒子的适应度值，即各粒子对应参数模型的合成地震记录与观测地震数据的误差的二次范数。

按步骤四计算模拟退火学习项，基于各粒子的接收概率，按轮盘选择算法确定该项的引导粒子pmin；按步骤五计算多维学习项，择优组合粒子位置，使其具有最大联合概率密度值，构成该项的引导粒子pmax，这里使用三参数模型的平滑结果作为测井数据，计算三参数的协方差矩阵。

对各粒子按步骤六进行速度更新和位置更新。

降低温度并进入下一次迭代，利用本发明方法对粒子不断迭代更新，直至终止温度和最大迭代次数，输出最优适应度的粒子，作为观测地震数据的反演结果。

图4为使用本发明方法的三参数反演结果与理论模型的对比，其中，(a)为纵波速度，(b)为横波速度，(c)为密度；可以看到反演结果(虚线)与理论结果(实线)吻合很好，特别是一些高频成分得到很好的恢复；图5为使用传统粒子群算法的三参数反演结果与理论模型的对比，其中，(a)为纵波速度，(b)为横波速度，(c)为密度；可以看到反演结果(虚线)与理论结果(实线)吻合较差。图6为分别使用传统粒子群算法与本发明方法进行反演的收敛速度的对比，可以看到使用本算法的收敛速度(实线)比使用传统粒子群算法的收敛速度(虚线)明显较快。

以上实施例仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明保护范围之内。