一种用于高维参数空间波形反演的协同变异差分进化算法

摘要

本发明公开了一种用于高维参数空间波形反演的协同变异差分进化算法，该算法将分解——协调的思想引入差分进化算，将高维个体分解为一系列的子成分，并引入局部适应度函数对每个子成分进行评价。然后在变异操作中利用局部适应度引导各子成分的变异方向，而在选择操作中利用全局适应度协调各子成分达到共同进化。与常用的快速模拟退火法和遗传算法相比，协同变异差分进化算法更加适合于高维参数空间波形反演；当层较薄、待反演参数多的情况下，协同变异差分进化算法能搜索到更加接近真实值的解。另外，协同变异差分进化算法的收敛速度对维数增加不敏感，维数很高时仍用很快的收敛速度。

说明书

技术领域

本发明属于地球物理勘探领域，涉及一种地球模型物性参数反演方 法，特别涉及一种用于高维参数空间波形反演的协同变异差分进化算法。

背景技术

随着油气田勘探开发程度的加深，油气勘探对象日趋复杂，逐渐由简单构 造油气藏向复杂构造油气藏转移，从构造油气藏向地层－岩性等隐蔽油气藏 转移。如何精细描述复杂构造和地层岩性是当前我国油气勘探面临的关键问 题。波形反演方法直接应用地震波形反演地下介质参数，可包含各种类型的 波，如多次波、转换波、面波等，充分地利用了地震波的旅行时、振幅、相 位和频率等信息。所以，波形反演是准确了解地下复杂结构和介质物理性质 （如速度、密度等）的潜在有效方法，近年来在勘探地球物理领域备受关注。

已有的波形反演方法大体上可分为两大类：一类是基于梯度的局部优化方 法，如梯度法、牛顿法和共轭梯度法等，这类方法容易陷入局部极值，结果 依赖于初始模型；另一类是全局优化方法，如模拟退火法、邻域法、遗传算 法^[4]、粒子群优化法等，这类方法不依赖于初始模型和函数梯度值，但计算量 较大。随着计算机计算能力的提高，全局优化方法得到越来越多的应用，但 计算时间仍是制约全局优化方法在波形反演中普遍应用的一个因素。另一方 面，现有的全局优化方法遭受到“维数瓶颈”问题，这些算法用于低维（维 数100以下）问题时具有卓越的寻优能力，但它们的搜索能力随着搜索空间 维数的增加会急剧下降。要精细刻画地下介质，必需对介质进行细小的划分， 这就使得模型参数非常庞大，反问题的维数高达几百几千，常用的全局优化 方法不再有效或需要大量的迭代。

协同进化算法是一种新的进化算法框架，该算法采用分解——协调的思 想将高维问题分解为一系列相互关联的子问题，每个子问题轮流交替进行优 化。协同进化算法对高维可分问题较为有效，但不适用于强非线性的高维波 形反演。差分进化（DE）算法是由Store和Price于1996年首先提出的基于 实数编码的连续空间全局优化算法，是一种原理简单易于实现而又能力强大 的模拟进化算法，之后，一些改进的DE算法被相继提出。差分进化算法成功 地应用到很多领域，如信号处理和地球物理反演^[。但差分进化算法同样受到 “维数瓶颈”问题，在高维参数空间波形反演应用中仍需改进。

发明内容

本发明的目的在于克服上述现有技术的缺点，提供一种用于高维参 数空间波形反演的协同变异差分进化算法，该算法利用局部适应度和全 局适应度同时引导算法的搜索方向，从而大大提高了收敛速度和全局寻 优能力。

本发明的目的是通过以下技术方案来解决的：

该种用于高维参数空间波形反演的协同变异差分进化算法，包括以 下步骤：

1）采集原始地震资料，然后对采集到的地震资料进行预处理，处理 后得到叠前共中心点道集或叠后地震数据，称这个地震数据为观测地震 数据，记为Seis(r,t)，其中r表示检波器接收点位置，t表示时间轴；

2）构建水平层状地质模型，给定地质模型的层数N和层的厚度， 每层的介质模型参数包纵波波速度V_p、横波波速度V_s和密度ρ；

3）确定地质模型参数V_p、V_s、ρ的搜索空间，并指定待优化的目标 函数；

4）在地质模型参数搜索空间内随机地生成NP个地质模型，并进行 实数编码得到含NP个个体的初始群体；

5）计算群体中每个随机模型的合成地震数据，然后根据目标函数估 计各随机模型对应个体的全局适应度值F_i；

6）进行变异操作，对群体中每个个体产生一个变异个体首先 将高维个体按照基因间相互依赖的强弱程度分解为一系列子成分，给每 个子成分指定一个局部适应度函数，根据局部适应度函数估计各子成分 的局部适应度值；然后随机选出三个不同个体，根据其相应子成分的局 部适应度值变化情况得到梯度信息，并以此梯度的负方向作为各子成分 的变异方向；

7）进行交叉操作，随机地从变异个体和与其相应的当代个体中 抽取基因，组合成一个试验个体

8）进行选择操作，根据当代个体和试验个体的全局适应度值选取下 一代；

9）进化代数g=g+1，判断是否满足终止条件，如果进化代数g小于 或等于设定的值G，则返回到步骤6）；否则执行步骤10）；

10）选出第G代群体中全局适应度值最小的个体，将该个体解码后 即得到最终搜索到的最优地质模型。

进一步，上述步骤3）中，参数V_p、V_s、ρ的搜索空间是根据测井 资料中相应参数的测井记录确定，即以测井得到的地质模型的低频分量 为背景，背景值增加设定值后为搜索空间上界，背景值减小设定值后为 搜索空间下界；或者纵波波速度V_p和横波波速度V_s的搜索空间也能够 根据高精度速度分析确定；所述目标函数是刻画最优解的标准，以计算 波场与观测波场之间的拟合程度或误差大小为标准。

上述步骤5）中，以观测地震数据与由模型合成的计算地震数据之 间误差绝对之和为目标函数，即目标函数为：

$F_i={\Sigma }_{r=1}^{\mathrm{NR}}\int |\mathrm{Seis}(r,t)-{\mathrm{Syn}}_i(r,t)|\mathrm{dt}.$

式中NR为检波器个数，适应度值F_i越小，表示观测地震数据与计 算地震数据之间误差越小，则其相应的个体就越优秀；反之，适应度 值F_i越大，则表示其相应的个体就越差。

上述步骤6）中，个体的变异方向是由子成分的局部适应度值来引 导：随机选出三个不同个体，将子成分的基因由小到大排列，如果其对 应的局部适应度值是递增的case1，则在最小基因的左方寻找新基 因；如果其对应的局部适应度值是递减的case2，则在最大基因的右 方寻找新基因；如果其对应的局部适应度值是先减后增case3，则在中 间基因的附近寻找新基因；如果其对应的局部适应度值是先增后减 case4，则在最小基因的左方或在最大基因的右方寻找新基因：

$v_{i,j}^g=\left(\begin{array}{cc}{x}_{r_{11},j}^g-\alpha |x_{r_2,j}^g-x_{r_3,j}^g|,& \mathrm{case}1\\ x_{r_{31},j}^g+\alpha |x_{r_2,j}^g-x_{r_3,j}^g|,& \mathrm{case}2\\ x_{r_{21},j}^g+\alpha (x_{r_2,j}^g-x_{r_3,j}^g),& \mathrm{case}3\\ p_1{x}_{r_{11},j}^g+p_2{x}_{r_{21},j}^g+p_3{x}_{r_{31},j}^g& \mathrm{case}4\end{array}\right)j=\mathrm{1,2},...,J$

$p_1=p_4/{\mathrm{LF}}_j\left(x_{r_{11},j}^g\right),p_2=p_4/{\mathrm{LF}}_j\left(x_{r_{21},j}^g\right),p_3=p_4/{\mathrm{LF}}_j\left(x_{r_{31},j}^g\right),$

$p_4=1/(L{F}_j{\left(x_{r_{11},j}^g\right)}^{-1}+{{\mathrm{LF}}_j\left(x_{r_{2}1,j}^g\right)}^{-1}+{\mathrm{LF}}_j{\left(x_{r_{31},j}^g\right)}^{-1}).$

其中α为变异尺度因子，一般选取0到1之间的实数，J为子成分 的个数，LF_j为各基因对应的局部适应度值。

进一步的，在上述步骤7）中，

$u_{i,n}^g=\left(\begin{array}{cc}{v}_{i,n}^g& \mathrm{ifrand}\le \mathrm{CRorn}=n_{\mathrm{rand}},\\ x_{i,n}^g& \mathrm{otherwise},\end{array}\right)n=\mathrm{1,2},...,D$

其中CR∈[0,1)是交叉率；rand为[0，1]均匀分布随机数；n_rand是1到 D之间的一随机整数。

上述步骤8）中，比较当代个体的全局适应度值和试验个体的全局 适应度值，全局适应度值小的个体被选入下一代，具体操作为：

$x_i^{g+1}=\left(\begin{array}{cc}{u}_i^g& \mathrm{ifF}\left(u_i^g\right)<F\left(x_i^g\right),\\ x_i^G& \mathrm{otherwise},\end{array}\right)i=\mathrm{1,2},...,\mathrm{NP}..$

本发明具有以下有益效果：

本发明将分解—协调的思想引入差分进化算，将高维个体分解为一 系列的子成分，并引入局部适应度函数对每个子成分进行评价。然后在 变异操作中利用局部适应度引导各子成分的变异方向，而在选择操作中 利用全局适应度协调各子成分达到共同进化。与常用的快速模拟退火法 和遗传算法相比，协同变异差分进化算法更加适合于高维参数空间波形 反演；当层较薄、待反演参数多的情况下，协同变异差分进化算法能搜 索到更加接近真实值的解。另外，协同变异差分进化算法的收敛速度对 维数增加不敏感，维数很高时仍用很快的收敛速度。

附图说明

图1是本发明流程示意图；

图2是本发明对高维个体的分解和局部适应度函数示意图；

图3是本发明中层状介质的分解示意图和计算局部适应度函数的时 窗示意图。

具体实施方式

下面结合附图对本发明做进一步详细描述：

地震波形反演的目的是获取一个最优的地层模型，使由该地层模型 计算的地震数据与实测地震数据拟合最佳。本发明采用一种改进的差分 进化算法来搜索最优地层模型参数。基于差分进化算法框架下，借鉴协 同进化算法中分解——协调的思想，把高维问题按参数相互依赖程度适 当地分解为一系列子问题，并引入局部适应度函数对各子问题进行评价， 利用局部适应度和全局适应度同时引导算法的搜索方向。

本发明的物质基础是通过野外高分辨率地震采集设备采集到的大量 地震数据。本发明的基于改进差分进化算法的波形反演框架如图1所示， 具体步骤分别为：

1）采集原始地震资料，然后对采集到的地震资料做常规预处理，包 括静校正处理、噪音压制处理、真振幅恢复处理等。处理后得到叠前共 中心点道集或叠后地震数据，称这个地震数据为观测地震数据，记为 Seis(r,t)，其中r表示检波器接收点位置，t表示时间轴。

2）构建水平层状地质模型，给定地质模型的层数N和层的厚度， 每层的介质模型参数包含V_p、V_s、ρ。这里，V_p、V_s、ρ分别为纵波波速 度、横波波速度和密度参数。

3）确定地质模型参数V_p、V_s、ρ的搜索空间，并指定待优化的目标 函数。参数搜索空间可根据测井资料中相应参数的测井记录确定，即以 测井得到的地质模型的低频分量为背景，背景值增加一定值后为搜索空 间上界，背景值减小一定值后为搜索空间下界。纵波波速度和横波波速 度的搜索空间也可根据高精度速度分析确定。目标函数是刻画最优解的 标准，一般以计算波场与观测波场之间的拟合程度或误差大小为标准。

4）在地质模型参数搜索空间内随机地生成NP个地质模型，并进行 实数编码得到含NP个个体的初始群体。生成的初始群体应尽可能地 均匀分布在整个搜索空间。每个个体可用下面的向量表示：

$x_i^g=\{x_{i,1}^g,x_{i,2}^g,...,x_{i,D}^g\}i=\mathrm{1,2},...,\mathrm{NP},$

其中上标g表示繁殖代数，初始群体时g=0，下标D为个体的维数，也 即待优化问题的参数个数。

5）计算群体中每个随机模型的合成地震数据Syn_i(r,t)，然后根据 目标函数估计各随机模型对应个体的全局适应度值F_i。

本发明以观测地震数据与由模型合成的计算地震数据之间误差绝对 之和为目标函数，即目标函数为：

$F_i={\Sigma }_{r=1}^{\mathrm{NR}}\int |\mathrm{Seis}(r,t)-{\mathrm{Syn}}_i(r,t)|\mathrm{dt}.$

式中NR为检波器个数。适应度值F_i越小，表示观测地震数据与计 算地震数据之间误差越小，则其相应的个体就越优秀；反之，适应度 值F_i越大，则表示其相应的个体就越差。

6）进行变异操作，对群体中每个个体产生一个变异个体首先 根据分解——协调的思想将高维个体按照适当条件分解为一系列子成 分，给每个子成分指定一个局部适应度函数，根据局部适应度函数估计 各子成分的局部适应度值，然后利用局部适应度值引导各子成分的变异 方向。

目前所有进化算法均使用一个统一的适应度函数来评价个体中的所 有基因（变量），这样的评价准则对高维问题不是很有效。本发明借鉴协 同进化算法和人工培育优良物种的思想，将高维个体分解为一系列子成 分，并引入局部适应度函数对每个子成分进行评价。为了便于区分，我 们将评价整个个体的适应度函数称为全局适应度函数。如图2所示，将 D维的个体分解为J个子成分，每个子成分记为为了便于描述， 假定每个子成分均包含M个基因。个体的表达式可重写为：

$x_i^g=\{x_{i,1}^g,x_{i,2}^g,...,x_{i,J}^g\}i=\mathrm{1,2},...,\mathrm{NP},$

$x_{i,j}^g=\{x_{i,j1}^g,x_{i,j2}^g,...,x_{i,\mathrm{jM}}^g\}j=\mathrm{1,2},...,J.$

本发明中，对于层状介质波形反演问题，可以将连续几层的模型参 数划分到同一个子成分里。如图3左栏所示，给出了将连续两层划分为 一个子成分的示意图。

然后，为每个子成分指定一个局部适应度函数LF_j，此局部适应度函 数必须满足以下两个条件：第j个局部适应度函数LF_j的变量至少包含第 j个子分成所有局部适应度函数在全局最优点处必须为最优适应度 值。对于波形反演，由于介质是根据旅行时分层，各层介质参数肯定与 其相应旅行时处的反射地震波形有必然联系，从而我们指定每个子成分 的局部适应度函数如下：

$L{F}_i={\Sigma }_{r=1}^{\mathrm{NR}}\int \left|\right[\mathrm{Seis}(r,t)-{\mathrm{Syn}}_i(r,t)]·{\mathrm{win}}_j\left(t\right)|\mathrm{dt}.$

其中win_j(t)是与第j个子成分相对应的时窗函数，如图3中间栏所 示，对于叠后数据，时窗的中心应位于子成分的中心，对于叠前数据， 随着偏移距的增加，时窗的中心应沿着图3右栏中红线所示的双曲线滑 行（该曲线根据经典Dix公式求取）。上面定义的局部适应度函数有一很 好的优点，只需一次正演即可计算一个个体的全局适应度值和所有局部 适应度值。因此，局部适应度的计算不需增加波场正演次数，这对正演 非常耗时的波形反演问题是非常重要的。

实际进行波形反演时也可将一个层看成一个子成分，即每个子成分 的参数只包含一层的速度和密度。如此分解的一个原因是使得算法易于 实现。局部适应度函数利用一个时窗截取部分反射波来评价各子成分， 要使局部适应度函数能正确地评价各子成分，必须选择适合的时窗宽度。 地面记录的地震数据中每个界面的反射波不是单个的脉冲而是具有一定 时宽的子波，且来自相邻界面的反射子波常常是相互重叠的。所以时窗 宽度的最佳值应为一个子波的宽度加上子成分的厚度（层按旅行时划 分）。地震子波的长度往往远大于层的厚度，这使得相邻时窗大部分是相 互重叠的。因此，虽然每个子成分只包含一层，但相邻的连续几层之间 是通过局部适应度函数关联起来的，相邻层之间的参数是协调进化的， 并且计算局部适应度值的时窗是逐层滑动的，从而层间的关联性也是逐 渐过渡的。

完成了分解和指定局部适应度函数后，即可利用局部适应度来引导 各子成分的变异方向。为了给目标个体产生一个变异个体，首先从其 余NP-1个个体中随机地选出3个个体i,r₁,r₂,r₃是[1，NP]之 间互不相等的随机整数。常用的变异策略是对个体进行整体操作，新变 异策略则以子成分为单位进行操作，即一个子成分接一个子成分地生成 变异个体。为了便于描述，假设每个子成分仅含一个变量，对第j个子 成分，先将3个随机个体对应子成分中的变量和按升序进 行排列：

$\left(\begin{array}{ccc}{x}_{r_{11},j}^g& x_{r_{21},j}^g& x_{r_{31},j}^g\end{array}\right)=\mathrm{sort}\left(\left(\begin{array}{ccc}{x}_{r_1,j}^g& x_{r_2,j}^g& x_{r_3,j}^g\end{array}\right)\right).$

排序后三个变量对应的局部适应度值记为这三个局部适应度值只可能出现四种排列情况，分别为：

Case1：${\mathrm{LF}}_j\left(x_{r_{11},j}^g\right)\le {\mathrm{LF}}_j\left(x_{r_{21},j}^g\right)\le {\mathrm{LF}}_j\left(x_{r_{31},j}^g\right).$

这种情况可以用一个增函数来拟合，增函数指示出在的左边很可能 存在一个更好的变量。

Case2：${\mathrm{LF}}_j\left(x_{r_{11},j}^g\right)\ge {\mathrm{LF}}_j\left(x_{r_{21},j}^g\right)\ge {\mathrm{LF}}_j\left(x_{r_{31},j}^g\right).$

此种情况可以由一个减函数来拟合，而减函数指示出在的右边 很可能存在一个更好的变量。

Case3：${\mathrm{LF}}_j\left(x_{r_{11},j}^g\right)<{\mathrm{LF}}_j\left(x_{r_{21},j}^g\right)\mathrm{and}{\mathrm{LF}}_j\left(x_{r_{21},j}^g\right)>{\mathrm{LF}}_j\left(x_{r_{31},j}^g\right).$

此种情况可以由一个凸函数来拟合，而凸函数意味着在的左边 或右边均有可能存在一个更好的变量。

$\mathrm{Case}4:{\mathrm{LF}}_j\left(x_{r_{11},j}^g\right)>{\mathrm{LF}}_j\left(x_{r_{21},j}^g\right)\mathrm{and}{\mathrm{LF}}_j\left(x_{r_{21},j}^g\right)<{\mathrm{LF}}_j\left(x_{r_{31},j}^g\right)..$

此种情况可以由一个凹函数来拟合，凹函数有一个极小值点，所以 它意味着其极小值点可能是一个更好的变量。但是拟合一个凹函数并精 确找到其极值点是非常麻烦的，与其相比，取的加权平均 是一个更好的选择。

总结上面四种情况，新的变异策略可以写为：

$v_{i,j}^g=\left(\begin{array}{cc}{x}_{r_{11},j}^g-\alpha |x_{r_2,j}^g-x_{r_3,j}^g|,& \mathrm{case}1\\ x_{r_{31},j}^g+\alpha |x_{r_2,j}^g-x_{r_3,j}^g|,& \mathrm{case}2\\ x_{r_{21},j}^g+\alpha (x_{r_2,j}^g-x_{r_3,j}^g),& \mathrm{case}3\\ p_1{x}_{r_{11},j}^g+p_2{x}_{r_{21},j}^g+p_3{x}_{r_{31},j}^g& \mathrm{case}4\end{array}\right)j=\mathrm{1,2},...,J$

$p_1=p_4/{\mathrm{LF}}_j\left(x_{r_{11},j}^g\right),p_2=p_4/{\mathrm{LF}}_j\left(x_{r_{21},j}^g\right),p_3=p_4/{\mathrm{LF}}_j\left(x_{r_{31},j}^g\right),$

$p_4=1/(L{F}_j{\left(x_{r_{11},j}^g\right)}^{-1}+{{\mathrm{LF}}_j\left(x_{r_{2}1,j}^g\right)}^{-1}+{\mathrm{LF}}_j{\left(x_{r_{31},j}^g\right)}^{-1}).$

其中α为变异尺度因子，一般选取0到1之间的实数。值得注意的 是，上面的变异操作在实施时并不需要进行任何曲线拟合。这个变异算 子同时拥有随机性和贪婪性。

7）进行交叉操作，随机地从变异个体和与其相应的当代个体中 抽取基因，组合成一个试验个体

$u_{i,n}^g=\left(\begin{array}{cc}{v}_{i,n}^g& \mathrm{ifrand}\le \mathrm{CRorn}=n_{\mathrm{rand}},\\ x_{i,n}^g& \mathrm{otherwise},\end{array}\right)n=\mathrm{1,2},...,D$

其中CR∈[0,1)是交叉率；rand为[0，1]均匀分布随机数；n_rand是1到 D之间的一随机整数，它保证了试验个体中至少有一个基因是来自于变 异个体。

8）进行选择操作，根据当代个体和试验个体的全局适应度值选取下 一代。比较当代个体的全局适应度值和试验个体的全局适应度值，全局 适应度值小的个体被选入下一代，具体操作为：

$x_i^{g+1}=\left(\begin{array}{cc}{u}_i^g& \mathrm{ifF}\left(u_i^g\right)<F\left(x_i^g\right),\\ x_i^G& \mathrm{otherwise},\end{array}\right)i=\mathrm{1,2},...,\mathrm{NP}.$

9）g=g+1，判断是否满足终止条件。如果进化代数g小于或等于规 定的值G，则返回到步骤6）；否则执行步骤10）。

10）选出第G代群体中全局适应度值最小的个体，将该个体解码后 即得到最终搜索到的最优地质模型。