CN103454677A - 基于粒子群与线性加法器结合的地震数据反演方法 - Google Patents

Abstract

基于粒子群与线性加法器结合的地震数据反演方法，包括步骤1.设定粒子空间及其初始参数；步骤2.产生粒子空间中各个粒子的个体极值及全局极值；步骤3.拟合计算出每个粒子的粒子适用值；步骤4.比较并优化粒子空间中各个粒子的个体极值及全局极值；步骤5.将步骤4中得到的新的个体极值与全局极值带回步骤2，重复步骤2-4，完成一次迭代，重复多次迭代，直到满足结束条件。采用本发明所述的基于粒子群与线性加法器结合的地震数据反演方法，使用粒子群算法仅需要调整两个参数即可实现对非线性问题的求解，得到反演结果并应用于地震数据资料解释中，以达到准确预测地层信息的目的。

Description

基于粒子群与线性加法器结合的地震数据反演方法

技术领域

本发明属于地质勘探领域，涉及石油勘探中的地震数据处理，特别是一种基于粒子群与线性加法器结合的地震数据反演方法。

背景技术

 地震反演是利用已经采集的地震资料，将已知地质信息和钻井、测井数据作为反演中的约束条件，对地下地层的空间展布情况和物理性质进行反演成像（求解）的过程。现有的地震反演模型，如测井约束反演、地震岩性模拟反演、广义线性反演、多道反演、地质统计学反演、遗传算法反演、混沌反演和波阻抗多尺度反演等都有各自的使用范围，如地质统计反演使用于各类复杂的地震预测和统计，尤其是钻井资料较多、需要进行精细储层描述的地区；波阻抗多尺度反演使用小波变换把目标函数分成不同尺度的分量，根据不同尺度上目标函数的特征逐步搜索全局最小点。

 测井约束反演测井约束反演是一种基于模型的反演，其反演结果的精度不仅依赖于研究目标的地质特征、钻井数量、井位分布以及地震资料的分辨率和信噪比，还取决于处理工作的精细程度，如测井资料的环境校正、层位的精细标定、子波提取等。

 广义线性反演广义线性反演（General Linear Inversion）也是一种建立在模型基础上的反演技术。它是通过模型正演与实际地震剖面作比较，根据误差的情况，最佳地逼近实际数据，从而迭代反复修改模型，直到满意为止。该算法将模型看作一个线性系统，其反演问题归结为求解一组线性联立方程组（例如用矩阵表示的误差方程组）。

地质统计学反演技术综合利用地震资料、地质知识和测井资料，通过高斯模拟、高斯协模拟和随机反演技术反演出各种储层参数，这一技术适用于各类复杂的地震预测和描述，尤其是钻井资料较多、需要进行精细储层描述的地区。

 波阻抗多尺度反演多尺度反演是一种加快收敛速度、克服局部极值影响、搜索全局最小点的反演策略。使用小波变换把目标函数分成不同尺度的分量，根据不同尺度上目标函数的特征逐步搜索全局最小点。

地震反演模型构造完成后，需要采取各种方法寻求模型参数的最优解。目前广泛采用的方法主要包括数学方法及智能方法。智能方法包括遗传算法反演、模拟退火方法反演等，它们都是基于生物仿真的智能算法，但普遍存在着参数较多、训练慢等缺点。

 地震资料的储层参数反演是一个非线性问题，这是因为地震反演中的目标函数是一个非凸性的复杂的多峰函数。如果使用对初始模型有很大依赖性的线性优化方法、最速下降法和共扼梯度法进行储层参数参数反演，则容易陷入局部极值而难以得到全局极值。并且当反演参数较多时，搜索空间急剧膨胀，计算量非常大，导致搜索效率大大下降。

 粒子群算法(Particle Swarm Optimization, PSO)的基本概念源于对鸟群捕食行为的简化社会模型的模拟，1995年由Kenndy 和Eberhart等人提出，它同遗传算法类似，通过个体间的协作和竞争实现全局搜索。系统初始化为一组随机解，称之为粒子。通过粒子在搜索空间的飞行完成寻优，在数学公式中即为迭代，它没有遗传算法的交叉以及变异算子，而是粒子在解空间追随最优粒子进行搜索。

 PSO算法的数学描述为：假设在一个d维的目标搜索空间中，有m个代表潜在问题解的粒子组成的一个种群                                               

，其中

表示第i个粒子在d维解空间的一个矢量点。将

代入一个与求解问题相关的目标函数可以计算出相应的适应值。用记录第i个粒子自身搜索到的最好点(所谓最好，是指计算得到的适应值为最小，即pbest)。而在这个种群中，至少有一个粒子是最好的，将其编号记为g，则

就是种群搜索到的最好值(即gbest)，其中

。而每个粒子还有一个速度变量，可以用

表示第i个粒子的速度。对每一个粒子，它的第D维(1 ≤D ≤ d)飞行速度和位置根据如下方程进行变化：

       (1)

                          (2)

其中：ω为惯性权值；c1和c2都是正常数，称为加速系数；rand()是在[0,1]范

围内变化的随机数；n为迭代次数。

此外，搜索时，微粒的位置被最大位置和最小位置限制，如果某微粒在某维的位置超出该维的最大位置或最小位置，则该微粒的位置被限制为该维的最大位置或最小位置。同样，微粒的速度也被最大速度和最小速度所限制，即有

                   (3)

线性加法器就是一个线性函数，用直线拟合曲线，整体来看就是把关系复杂的函数图像，分成简单的线性关系，例如如果x是一维，就把曲线划分成很多个小线段，如果x是两维，就把曲面划分成很多个小平面等。通过这样的划分就可以转换成寻找若干个加法器的系数，通过PSO算法和训练数据可以完成这一工作。

发明内容

为克服现有地震数据反演算法反演参数多，搜索空间膨胀，计算量巨大的技术缺陷，本发明公开了一种基于粒子群与线性加法器结合的地震数据反演方法。

本发明所述基于粒子群与线性加法器结合的地震数据反演方法，包括测量一一对应的多个数据点（Xi，Yi），其中Xi，Yi 均为向量，i=1、2、…m，m为数据量个数，还包括如下步骤：

步骤1.设定解空间的惯性权重参数ω，加速系数c1、c2，设定粒子群中各个粒子的初始个体极值和粒子群的初始全局极值，设定各个粒子的初始速度；

步骤2.以Xi为基础，利用①式和②式，随机产生各个粒子的初始位置

与初始速度

，每个粒子的初始个体极值作为其初始位置， 全部个体极值中的最佳值为初始全局极值；

      -----①

                                   ------②

其中

为Xi在第d维的分量，

为Xi对应的个体极值在第d维的分量，

为全局极值在第d维的分量，rand()是在[0,1]范围内变化的随机数；上标n或n+1为迭代次数，n=0时，表示相应参数的初始值； 

步骤3.根据步骤2中的Xi拟合出

，将与Yi比较，计算出每个粒子的粒子适用值；

步骤4.比较粒子适用值与它的个体极值，若前者优于后者，则以粒子适用值作为该粒子新的个体极值；比较各个粒子的粒子适用值与全局极值，若前者中有优于后者的粒子，则以该粒子的粒子适用值作为新的全局极值；

步骤5.将步骤4中得到的新的个体极值与全局极值带回步骤2，重复步骤2-4，完成一次迭代，重复多次迭代，直到满足结束条件。

    具体的，步骤3中，计算Xi粒子的粒子适用值

为利用下式；

      

       

其中

为根据全部Xi拟合出的训练值，i=1、2、…m 。

 进一步的，所属步骤3中利用多层线性计算器对

进行拟合，具体为：

步骤301.利用 ，计算出

，其中

，计算出每一维度的分量，从而计算出

；

其中b为预设常数，下标i表示第i个粒子，d表示粒子的维度。 上标n表示迭代次数；

步骤302.利用

 计算出

。

 优选的，所述步骤5中的结束条件为满足预设的迭代次数或得到的全局极值小于预设值。

 具体的，C1=C2=2、ω=1。

优选的，m取值为20-40。

采用本发明所述的基于粒子群与线性加法器结合的地震数据反演方法，使用粒子群算法仅需要调整两个参数即可实现对非线性问题的求解，从而建立起测井目标曲线与地震波形间的非线性映射关系，得到反演结果并应用于地震数据资料解释中，以达到准确预测地层信息的目的。

附图说明

图1为本发明所述多层线性计算器的步骤结构示意图；

图2为本发明所述多层线性计算器的域门和软最大函数的步骤关系示意图；

图3示出本发明一种具体实施方式的步骤示意图。

具体实施方式

下面结合附图，对本发明的具体实施方式作进一步的详细说明。

本发明所述基于粒子群与线性加法器结合的地震数据反演方法，包括测量一一对应的多个数据点（Xi，Yi），Xi，Yi 均为向量，i=1、2、…m，m为数据点个数，还包括如下步骤：

步骤1.设定解空间的惯性权重参数ω，加速系数c1、c2，设定粒子群中各个粒子的初始个体极值和粒子群的初始全局极值；

步骤2.以Xi为基础，利用①式和②式，随机产生各个粒子的初始位置

与初始速度

，每个粒子的初始个体极值作为其初始位置， 全部初始个体极值中的最佳值为初始全局极值；

      -----①

                                   ------②

其中

为Xi在第d维的分量，

为Xi对应的个体极值在第d维的分量，

为全局极值在第d维的分量，rand()是在[0,1]范围内变化的随机数；上标n或n+1为迭代次数，n=0时，表示相应参数的初始值； 

步骤3.根据Xi拟合出

，将

与Yi比较，计算出每个粒子的粒子适用值；

步骤4.比较粒子适用值与它的个体极值，若前者优于后者，则以粒子适用值作为该粒子新的个体极值；比较各个粒子的粒子适用值与全局极值，若前者中有优于后者的粒子，则以该粒子的粒子适用值作为新的全局极值；

步骤5.将步骤4中得到的新的个体极值与全局极值带回步骤2，重复步骤2至4，完成一次迭代，重复多次迭代，直到满足结束条件。

粒子群优化算法是在一个解空间的多个随机解中，通过叠代找到最优解的算法。在每一次叠代中，每个粒子，也就是随机解通过跟踪两个极值参数来更新自己。第一个极值参数是粒子本身所找到的最优解，叫做个体极值pbest。另一个极值是整个种群目前找到的最优解，这个极值参数是全局极值gbest。

例如，在一个d维的解空间中，有由m个代表随机解的粒子组成的一个种群

，其中

，

表示第i个粒子在d维解空间内的一个矢量点。

是一个向量，包含d个维度，每一维度上的值为

。

将

代入一个与求解问题相关的目标函数可以计算出相应的适应值。用

记录第i个粒子自身搜索到的最好点，所谓最好，是指计算得到的适应值为最小，即个体极值pbest，每个粒子在D维空间的每一维度都有一个个体极值。

对随机粒子的初始化采用随机过程，利用①式和②式，随机产生各个粒子在每一维度的初始位置

与初始速度

，每个粒子的初始个体极值作为其初始位置， 全部个体极值中的最佳值为初始全局极值。

      -----①

                                   ------②

① 和②式中的惯性权重参数ω，加速系数c1、c2根据经验设定，一般

取c1=c2=2，ω=1，解空间中的粒子个数m一般取在20-40，粒子数目越多，算法搜索的空间范围越大，也就越容易发现全局最优解。但是，算法运行的时间会较长。 利用①式和②式，遍历每一个维度，得到一个Xi的全部向量值。

将步骤1中设定的各个参数和初始值代入①式和②式，开始进行迭代，对输入Xi的每一维度分别进行计算更新。

更新后的值代入步骤3进行拟合，计算出

步骤3中，计算Xi粒子的粒子适用值

为利用④式；

      

       

-------④

其中

为根据全部Xi拟合出的训练值，i=1、2、…m 。

所属步骤3中可以利用多层线性计算器对进行拟合。

多层线性计算器是由多个线性计算器和多个域门组成，对于每个输入，都需要代入各个线性计算器，并由域门动态的控制各层的输出权重，最后逐层加权计算得到最终的输出。多层线性计算器模型将输入空间分成一个嵌套的子空间集，在多个以分层方式调整的域门控制下，信息在计算器之间被整合或者重新分配。多层线性计算器模型的体系结构类似树形，域门在树的非终端节点，充当协调者的角色，而多个线性计算器作为树叶终端。如图1-2所示。所谓域门是多层线性计算器中根据输入量得到多层线性计算器中间计算量的第一次计算过程，而软最大是对域门的中间计算量进行权重优化处理的过程。

软最大函数的具体形式为： 

                  

其中K为输入个数，u_k是输入向量X的一个线性加权输出，加权系数向量为a_k，偏置量为c_k，则

。可以看出软最大化函数满足以下要求：

                       

                       

具体在本发明中：

 对输入数据

，其中

，

 利用

 ，计算出

，其中，计算出每一维度

的分量，从而计算出

，这里 b为常数，具体取值意义不大，一般取零即可。下标中的i表示第i个粒子，d表示粒子的维度。 上标n表示迭代次数，由于步骤3为每次迭代都饱含的步骤，因此对不同的迭代次数，使用不同的参数。

            利用

 ；

 计算出

。计算出

后，代入④式求出粒子适用值。

随后进行步骤4.比较粒子适用值与它的个体极值，若前者优于后者，则以粒子适用值作为该粒子新的个体极值；比较各个粒子的粒子适用值与全局极值，若前者中有优于后者的粒子，则以该粒子的粒子适用值作为新的全局极值。

 步骤5.将步骤4中得到的新的个体极值与全局极值带回步骤2，重复步 骤2至4，完成一次迭代，重复多次迭代，直到满足结束条件。

 所述步骤5中的结束条件可以为满足预设的迭代次数或得到的粒全局极值小于预设值。例如将迭代次数设置为10次，或当全部粒子的适用值小于0.001时。由于本方法中每一个粒子意味着一个近似解，全局极值表征着近似解与真实解之间的差，因此全局极值越小，表示解越接近真实值。

 为更清楚简要的总结说明，以图3为例示出本发明一个具体实施方式，首先初始化粒子群空间，利用多层线性计算器拟合出训练值

，并与测量值，即图3中的观测值Y一起计算出粒子适用值和全局极值gbest，将gbest与预先设定的值e比较，若大于e，则更新后重新迭代，直到达到迭代次数或使gbest小于e，或gbest与e的误差在预设范围之内。

 采用本发明所述的基于粒子群与线性加法器结合的地震数据反演方法，使用粒子群算法仅需要调整两个参数即可实现对非线性问题的求解，减小了计算量，并能建立起测井目标曲线与地震波形间的非线性映射关系，得到反演结果并应用于地震数据资料解释中，以达到准确预测地层信息的目的。

利用多层线性计算器模型通过线性计算器的耦合交叉实现对地震资料解释的非线性反演。与粒子群算法结合后，利用多层线性计算其对粒子群算法中的个体极值和全局极值两个参数进行拟合优化，一方面需要优化的参数数量只有两个，减小了计算量，同时采用多层线性计算器使优化精度提高，保证了结果的可靠性。

前文所述的为本发明的各个优选实施例，各个优选实施例中的优选实施方式如果不是明显自相矛盾或以某一优选实施方式为前提，各个优选实施方式都可以任意叠加组合使用，所述实施例以及实施例中的具体参数仅是为了清楚表述发明人的发明验证过程，并非用以限制本发明的专利保护范围，本发明的专利保护范围仍然以其权利要求书为准，凡是运用本发明的说明书及附图内容所作的等同结构变化，同理均应包含在本发明的保护范围内。 

Claims (6)

1.基于粒子群与线性加法器结合的地震数据反演方法，包括测量一一对应的多个数据点（Xi，Yi），其中Xi，Yi均为向量，i=1、2、…m，m为数据量个数，其特征在于，还包括如下步骤：

步骤1.设定解空间的惯性权重参数ω，加速系数c1、c2，设定粒子群中各个粒子的初始个体极值和粒子群的初始全局极值，设定各个粒子的初始速度；

步骤2.以Xi为基础，利用①式和②式，随机产生各个粒子的初始位置                                                

与初始速度

，每个粒子的初始个体极值作为其初始位置， 全部个体极值中的最佳值为初始全局极值；

      -----①

                                   ------②

其中

为Xi在第d维的分量，

为Xi对应的个体极值在第d维的分量，

为全局极值在第d维的分量，rand()是在[0,1]范围内变化的随机数；上标n或n+1为迭代次数，n=0时，表示相应参数的初始值； 

步骤3.根据步骤2中的Xi拟合出

，将

与Yi比较，计算出每个粒子的粒子适用值；

步骤4.比较粒子适用值与它的个体极值，若前者优于后者，则以粒子适用值作为该粒子新的个体极值；比较各个粒子的粒子适用值与全局极值，若前者中有优于后者的粒子，则以该粒子的粒子适用值作为新的全局极值；

步骤5.将步骤4中得到的新的个体极值与全局极值带回步骤2，重复步骤2-4，完成一次迭代，重复多次迭代，直到满足结束条件。

2.如权利要求1所述基于粒子群与线性加法器结合的地震数据反演方法，其特征在于，步骤3中，计算Xi粒子的粒子适用值为利用下式；

              

其中

为根据全部Xi拟合出的训练值，i=1、2、…m 。

3.如权利要求2所述基于粒子群与线性加法器结合的地震数据反演方法，其特征在于，所属步骤3中利用多层线性计算器对

进行拟合，具体为

步骤301.利用

 ，计算出

，其中

，计算出每一维度的分量，从而计算出

；

其中b为预设常数，下标i表示第i个粒子，d表示粒子的维度，上标n表示迭代次数；

步骤302.利用 计算出。

4.如权利要求1所述基于粒子群与线性加法器结合的地震数据反演方法，其特征在于，所述步骤5中的结束条件为满足预设的迭代次数或得到的全局极值小于预设值。

5.如权利要求1所述基于粒子群与线性加法器结合的地震数据反演方法，其特征在于，C1=C2=2、ω=1。

6.如权利要求1所述基于粒子群与线性加法器结合的地震数据反演方法，其特征在于，m取值为20-40。

Images