混合粒子群算法在地震波阻抗反演中的应用

刘玉敏, 高松岩

(东北石油大学 电气信息工程学院, 黑龙江 大庆 163318)

0 引 言

随着油田储量的日渐枯竭以及复杂的油气藏很难被识别, 因而从地震资料中了解地下构造的基本情况, 对地震资料进行综合处理解释就显得尤为重要。地震反演技术就是按照一定的已知的地质规律, 将地表观测到的地震数据进行反推, 从而推知地球内部的结构形态和物质成分[1]。波阻抗是重要的岩石物理参数, 能较好地提供地层速度的信息, 可直接与钻井资料对比进行储层岩性解释和物性分析, 因而波阻抗反演技术在油藏描述、 储层参数估算、 油气横向预测等研究工作中发挥着越来越重要的作用[2]。

地震波阻抗反演是典型的非线性化问题, 传统的地震波阻抗反演常用线性化算法, 如最快下降法, 共轭梯度法, 牛顿法等[3]。线性优化方法存在的弊端就是其目标函数大部分需要可导或可微, 这就使一些函数难以被优化、 收敛速度慢、 容易陷入局部最优解等问题[4,5]。因此目前许多学者都在尝试将智能非线性优化算法应用到地震波阻抗反演中。如文献[6]将蚁群算法应用到地震波阻抗反演中； 文献[7]提出了层状模型约束下的粒子群算法, 并将其应用到波阻抗反演中； 文献[8]将免疫算法与遗传算法相结合应用于波阻抗反演； 文献[9]将遗传算法和模拟退火相结合, 应用到波阻抗反演中, 取得了良好的效果； 文献[10]将粒子群算法用于火山岩的波阻抗反演中, 并将反演结果与神经网络方法相比较, 证实了粒子群方法的有效性。

粒子群优化算法是一种基于群体的随机优化技术, 是典型的非线性算法[11]。该算法具有实现方便、 收敛速度快、 参数设置少等优点。目前该算法已经广泛应用于各个领域, 并且许多学者在不断从优化算法结构、 参数等方面对传统粒子群进行改进。文献[12]将粒子群算法成功用于组合优化问题； 文献[13]针对粒子群优化算法在搜索的初期收敛较快, 而在后期容易陷入局部最优的缺点, 对粒子群算法进行了优化和仿真； 文献[14]对传统的粒子群算法的惯性权重进行了优化, 使权值能自适应改变, 从而提高了粒子群搜索的速度； 文献[15]采用了一种时变学习因子的方式设计了一种自组织粒子群优化器； 文献[16]提出了一种基于改进粒子群算法的弹性阻抗非线性反演方法, 获得了良好的应用效果。也有很多学者将其他优化算法与粒子群算法相结合, 互相取长补短, 优化其性能。如与神经网络、 模拟退火、 混沌、 量子等相结合[17-20]。笔者在上述文献的基础上, 考虑将粒子群算法与混沌和遗传算法相结合。混沌序列具有遍历性, 混沌优化方法具有全局渐进收敛、 易跳出局部极小点和收敛速度快等特点[18], 因此考虑将混沌思想引入到粒子群中。但是混沌粒子群优化算法又具有不彻底性, 盲目随机性, 因此将遗传算法嵌入到粒子群算法中, 可以增强种群之间的互相交流, 提高种群的多样性[21]。将笔者算法应用于波阻抗反演中, 通过对模型和实际数据的反演, 结果表明笔者算法具有较快的收敛速度和较高的反演精度。

1 算法描述

1.1 粒子群算法基本原理

粒子群算法(PSO: Particle Swarm Optimization)与其他进化算法一样, 也采用群体和进化的概念。它通过个体间的协作与进化, 实现在复杂空间中搜索最优解的过程。算法将每个个体看作是n维搜索空间中的一个没有体积和重量的粒子, 每个粒子的位置都为优化问题的一个潜在解, 粒子在解空间中飞行, 由速度决定其运动方向和所在位置, 并由目标函数为之确定一个适应值, 适应值的大小衡量粒子的优劣。在搜索过程中, 粒子通过跟踪本身迄今找到的最优位置和全种群迄今找到的最优位置不断寻优, 经逐代搜索直到满足结束条件。具体过程如下[11,13]。

设在一个d维的目标搜索空间中, 有m个粒子组成一个群体。其中在第t次迭代时粒子i的位置表示为xi(t)=(xi1(t),xi2(t),…,xid(t))； 相应的飞行速度表示为vi(t)=(vi1(t),vi2(t),…,vid(t))； 粒子本身迄今搜索到的最优解(个体极值)表示为pi(t)=(pi1(t),pi2(t),…,pid(t))； 整个粒子群到目前为止找到的最优解(全局极值)表示为pg(t)=(pg1(t),pg2(t),…,pgd(t))。在每次迭代中, 粒子通过跟踪这两个极值更新自己的速度和位置。在第t+1次迭代计算时, 粒子i的第k维分量根据下列规则更新自己的速度和位置

vik(t+1)=wvik(t)+c1r1(pik(t)-xik(t))+c2r2(pgk(t)-xik(t))(1)

xik(t+1)=xik(t)+vik(t+1)i=1,2,…,m,k=1,2,…,d(2)

式(1)中,w为惯性权重, 体现了前一次迭代的速度对下一次迭代速度的影响。w大全局搜索能力强,w小局部搜索能力强。c1和c2为两个学习因子, 体现了粒子向个体最优解和全局最优解学习的步长。r1和r2为两个均匀分布在(0,1)之间的随机数。

1.2 遗传算法基本原理

遗传算法(GA: Genetic Algorithms)是一种以自然选择和基因遗传理论为基础的随机优化搜索算法[22]。算法结合了生物进化过程中的适者生存规则与群体内部染色体的随机信息交换机制。具体过程为： 首先用一定长度和编排规则的数字或字符将每个可能的解编码为一个染色体个体。在搜索之初, 从搜索空间中随机地选择一组染色体作为初始群体； 选定目标函数计算群体中所有个体的适应度值, 以度量个体的优劣； 根据适应度值从中选取优秀个体作为父本, 这体现了适者生存的进化理论； 由交换概率挑选的每两个父本通过将相异的部分基因进行交换, 得到新的个体, 体现了信息交换的思想； 为了在群体中能够引入新的个体, 用极小的概率对个体进行变异操作, 为新个体的产生提供机会。算法通过不断地使用选择、 交叉、 变异等遗传操作, 使获得的染色体一代更比一代好, 最终得到适应度最佳的解[23]。

1.3 混合粒子群算法

通过对粒子群算法和遗传算法进行分析、 比较, 发现这两种算法都是在自然特性的基础上模拟个体种群的适应性, 然后采用一定的变换规则通过搜索空间求解的。主要不同之处如下。 1) PSO算法在计算适应度后, 所有个体都被更新, 形成下一代个体, 使算法具有了良好的记忆性, 但搜索速度受到影响。而遗传算法是通过某种选择机制, 选出若干个样本, 两两杂交生成下一代的新个体, 以前的知识随着种群的改变被破坏, 记忆性差, 但杂交操作体现了算法信息交换的思想。2) 遗传算法有变异操作, 虽然变异概率很低, 但这为新个体的产生提供了机会, 更容易跳出局部最优。而粒子群算法没有变异操作。可见两种算法基本流程类似, 但各有利弊。笔者以PSO算法为主框架, 在原粒子群更新操作的基础上, 加入了遗传算法中的选择、 交叉、 变异思想。采用这种结合方式, 粒子群优化算法可以通过学习借鉴以前的搜索经验, 保存所有粒子好的解的知识, 减少无效迭代的次数, 提高收敛速度。也可以发挥遗传算法在保持种群多样性和全局搜索的优势。此外, 由于粒子群算法的初值是随机生成的, 不能保证种群的多样性和粒子搜索的遍历性。因此在初值选取时, 还引入了混沌序列的概念。混沌是一种看上去没有规律的运动, 它是存在于确定性的非线性系统中的类似随机的行为。它具有遍历性、 随机性和规律性等特点, 能根据其自身的规律, 在一定范围内不重复遍历整个状态。因此具有易跳出局部最优、 搜索速度快、 全局渐进收敛等优点[24]。所以笔者在初值选取时, 引入了混沌序列代替随机初始化。粒子的位置初值和速度初值以混沌序列形式排列, 既具有随机性的本质, 又具有混沌特性, 从而提高了种群的多样性和粒子搜索的遍历性。笔者混合粒子群算法具体步骤如下。

1) 参数初始化。确定粒子群算法的种群规模m、 种群进化次数以及学习因子等参数； 利用logistic混沌映射模型产生d个混沌变量

(3)

(4)

2) 计算粒子群的每个粒子的适应度值, 将粒子的个体极值记为pbest, 全局极值记为gbest。

3) 根据粒子群优化算法式(1)和式(2)更新粒子的速度和位置, 重新计算粒子的适应度值, 并更新pbest和gbest。

4) 对适应度进行排序, 并择优选取样本作为下一步遗传操作的父本。

5) 执行遗传算法的交叉, 变异操作, 生成新一代种群。

6) 判断是否满足算法的停止条件, 若满足, 算法结束, 输出最优解； 否则返回2)。

混合粒子群算法的具体流程如图1所示。

图1 混合粒子群算法流程图Fig.1 Flow chart of hybrid particle swarm optimization

2 粒子群波阻抗反演

2.1 地震波阻抗反演

目前常用的波阻抗反演方法基本上都是建立在褶积模型的基础上。地震记录用褶积可表示为

S(t)=R(t)*ω(t)+N(t)(5)

2.2 混合粒子群波阻抗反演思路

由褶积模型(5)可以看出, 地震记录与波阻抗为非线性关系, 因此求取波阻抗的过程就是一个求解非线性问题的过程。粒子群算法是一个典型的求解非线性问题的算法, 将其应用于地震道的波阻抗反演, 问题的关键在于如何构建一个能够对反演结果进行评价的目标函数[4]。笔者在简化模型的基础上, 不考虑平衡白噪声, 直接利用具有一定抗噪能力的最小二乘原理建立目标函数, 其数学表达式如下

(6)

其中S(t)为实际地震记录,N为采样点数,ω(t)为提取的地震子波,R(t)为所求的反射系数。

标准PSO波阻抗反演过程中, 每个粒子的位置, 代表着每个波阻抗模型。每个粒子的速度, 则代表着波阻抗模型的修正量。PSO的任意一个粒子在波阻抗空间里进行搜索, 每搜索到一个新位置就相当于构建了一个新的波阻抗模型。将波阻抗根据式(5)与子波褶积, 即可求出模型响应S。在每步迭代过程中依据粒子群速度和位置更新公式、 模型以及模型的修正量, 并最终达到目标函数最小。

笔者混合粒子群算法反演波阻抗过程具体描述如下。

1) 确定求解模型空间。在分析地质、 测井等先验信息的基础上, 建立反演初始模型。以该模型为依据, 随机产生一组波阻抗合理解, 构成模型求解空间。并利用混沌发生器Logistic映射使之成为一个混沌变量, 也即粒子。粒子的维数对应地下介质的层数。

2) 设定算法参数。为每个粒子赋一个初始速度v, 作为波阻抗模型的修正量初值。同时为惯性权重、 学习因子、 交叉变异概率等参数设置初值。

3) 计算每个粒子的适应度。目标函数值越小, 模型响应越接近实际地震记录。记录每次迭代粒子个体最优位置pbest和所有粒子的最优位置gbest。

4) 迭代寻优。采用粒子群位置和速度公式对粒子进行更新, 并对粒子进行选择、 交叉、 变异这3个遗传操作。粒子通过反复追随自身历史最优位置和种群最优位置, 并以一定的概率进行信息交换, 从而逐步趋近全局最优解模型。

5) 波阻抗的输出。当满足迭代终止条件(达到理想精度或者最大迭代次数时)时, 停止迭代, 输出种群最优位置gbest即为最终反演的波阻抗模型值。否则返回3)。

3 算例分析

3.1 标准粒子群波阻抗反演

笔者设计了一个一维的五层层状介质波阻抗模型, 如图2所示。该模型各层介质波阻抗为1.5,1.6,2.0,1.7,1.8, 单位为106g·m-2·s-1。模型空间采样间隔为2 ms, 有60个采样点。选择雷克(Ricker)子波合成地震记录, 褶积建立合成地震记录如图3所示。雷克(Ricker)子波是最常用的一种子波。它有一个幅值较大的主瓣, 两个幅值较小、 对称的旁瓣, 非常接近真实地震子波。雷克子波的数学表达式为f(t)=[1-2(∂ft)2]e-(πft)2。笔者选择的Ricker子波频率为50 Hz, 采样率为0.002 s。

图2 波阻抗模型 图3 合成地震记录 Fig.2 Wave impedance model Fig.3 Synthetic seismogram

将标准粒子群算法用于上述模型的波阻抗反演中。反演波阻抗结果如图4所示, 反演地震记录如图5所示。从反演的波阻抗和地震记录中, 可以看出反演的波阻抗模型趋势与实际模型一致, 基本能反映出模型的各个层位波阻抗的变化。但是反演结果与实际模型相比还是存在一定的误差, 误差精度在104, 误差曲线如图6所示。

图4 标准粒子群波阻抗反演 图5 标准粒子群反演地震记录 Fig.4 Inversive wave impedance with Fig.5 Inversive seismogram with ordinary ordinary particle swarm optimization particle swarm optimization

3.2 混合粒子群算法波阻抗反演

由于标准粒子群算法具有后期收敛速度较慢, 容易早熟且易陷于局部极小值的局限性, 因此采用笔者的混合粒子群算法对上述模型进行反演。设计学习因子大小为1.494 45, 种群为60, 惯性权重采用线性递减惯性权重ω(t)=0.9-0.5t/T, 迭代次数T为1 000次, 交叉概率为0.5, 变异概率为0.05。反演结果如图7～图9所示。

图6 标准粒子群波阻抗反演误差曲线 图7 混合粒子群波阻抗反演结果 Fig.6 Error curve with ordinary Fig.7 Inversive wave impedance with particle swarm optimization hybrid particle swarm optimization

图8 混合粒子群反演地震记录 图9 混合粒子群反演误差曲线 Fig.8 Inversive seismogram with hybrid Fig.9 Error curve with hybrid particle swarm optimization particle swarm optimization

从图7～图9可见, 在相同参数设置的前提下, 混合粒子群算法反演出的波阻抗、 地震记录与模型的吻合度更高, 虽然在部分区域存在误差, 但在误差量级上有了大幅度的减小。这足以说明, 笔者的混合粒子群算法对波阻抗反演有较好的效果, 精度较高, 具有一定的可信度。

3.3 算法抗噪性分析

考虑到实际环境中噪声的存在, 笔者在理论地震记录中加入15%的随机噪声, 对算法进行了抗噪性试验与分析。

图10 反演波阻抗(15%噪声) 图11 反演地震记录(15%噪声) Fig.10 Inversive wave impedance with 15% noise Fig.11 Inversive seismogram with 15% noise

由图10, 图11可以看出, 当地震记录中加入噪声后, 反演结果仍然能跟踪模型参数, 可以较真实的反映有效信息。噪声显然对结果造成了影响, 使某些区域反演误差有所增加, 但还在可接受的范围之内。这说明笔者混合粒子群算法的波阻抗反演具有一定的抗噪能力。

4 应用实例

笔者利用混合粒子群算法对位于松辽盆地某区块的实际地震资料进行波阻抗反演。该数据采样间隔0.5 ms, 地震道为60道, 每道有100个采样点, 采用40 Hz的雷克子波作为提取的子波。反演结果如图12所示。

图12 实际地震资料反演波阻抗剖面Fig.12 Inversive wave impedance section with seismic data

从该剖面的储层预测结果可以看出, 目标层段内储层纵向各薄层砂体在波阻抗反演剖面上显示效果良好, 预测结果未受到地震资料宏观模型影响, 井间薄层砂体连续性关系得到了清晰刻画, 对于油田范围内薄层储层预测具有很好的参考价值。

5 结 语

笔者针对传统粒子群算法在波阻抗反演中效率低, 误差偏大的缺点, 将混沌和遗传算法引入到粒子群中, 并对地震波阻抗模型进行了反演。结果表明, 笔者采用的混合粒子群算法在收敛精度上, 明显优于传统粒子群反演。通过在模型中加入噪声, 验证了笔者的混合粒子群算法也具有良好的抗噪声效果。利用混合粒子群算法对实际地震资料进行波阻抗反演, 获得了残差较小的反演地震记录, 说明笔者方法具有一定的应用价值。