CN109782354B - 基于方向引导的协同差分进化算法及其在射线追踪中的应用 - Google Patents

Abstract

本发明公开了一种基于方向引导的协同差分进化算法及其在射线追踪中的应用，在变异策略中引入了一个符号函数来定义变异的方向以提高算法收敛性，且提出两种灵活有效的迭代终止准则以改进计算效率。不同于局部优化波形反演算法和其他协同进化算法，本发明提出的方法不依赖于初值，对高维可分模型空间有较好的可行性，计算速度快，收敛性好，对维度较大的参数依然能有效反演；并且，在变异策略中引入S函数，S函数的值的正负特征决定变异的方向；引入S函数使得变异沿着缩小误差能量的方向进行，使算法加速收敛，能快速收敛到全局最优解。

Description

基于方向引导的协同差分进化算法及其在射线追踪中的应用

技术领域

本发明属于地震勘探技术领域，涉及射线追踪，尤其是一种基于方向引导的协同差分进化算法及其在射线追踪中的应用。

技术背景

在地球物理中，我们研究地震波在地下介质中的传播来研究地层结构。射线是地震波穿过地下介质的路径，是波动方程的高频近似解。射线追踪是层析成像最关键的步骤，用于迭代反演的初始阶段来模拟正演模型。快速有效的射线追踪算法是提高反演效率的关键。常用的射线追踪算法有波前法、最短路径法、走时法、有序波前重构法等。上述算法计算效率低，收敛性差，分辨率低。最短走时射线追踪方法将地下介质网格化，然后遍历所有结点。网格不能太小，否则追踪速度将很低。但是网格太大将降低分辨率。步进迭代法必须预先假设一条射线，这将引入人为因素。射线方程方法描述地震能量随高频的轨迹。因为射线方程有解析表达式，所以射线追踪将会比较快且更加精确。

本发明基于叠前地震记录研究射线方程方法。炮点和接收点的位置是已知的，如果反射点和出射角也已知，我们就可以确定一条射线。本发明利用迭代方法计算射线的出射角。学者们研究了多种反演方法，如局部最优波形反演方法和全局最优方法等。局部最优方法中，目标函数的最小值可由梯度法、共轭梯度法、牛顿法、高斯牛顿法和伪牛顿法等确定。在这些方法中，基于梯度的波形反演取得了一些成功，但是这类方法具有非线性性，且依赖于初值。Virieux指出波形反演方法的困难是建立精确的初始模型、确定新的最小化准则及改进多参数反演能力等。

此外，有许多全局优化算法也取得了进展，如进化算法、模拟退火法、蒙特卡洛方法等。相对于基于梯度的波形反演方法，全局优化算法不依赖于初值且不需要梯度，但是不足是计算量大。基于此，随着计算机运算能力的大幅提升，全局优化算法得到广泛应用。

另外，认知计算能处理类脑问题。它在学习、推理和解决特殊问题方面有强大的能力。认知系统能充分理解问题，且将思考的过程建模成计算模型。他们研究特征提取的能力类似大脑的学习过程。因此认知计算能加强人工智能和决策能力。

进化算法类似认知计算，它用计算机模拟自然物种进化过程，使个体成为最优。但是，传统的进化算法解决高维问题时效率较低。潘志斌提出CRsADE方法，该方法采用个体交叉率和子分量交叉率自适应地改进交叉运算地效率。汪超提出协同选择差分进化算法(DE-CCS)和协同变异差分进化算法(DE-CCM)。这些方法将高维空间问题分解为几个子空间并用局部适应汗水评价子问题。R.Chandra提出协同差分多任务学习方法用于多步混沌时间序列预测。高照奇提出高效进化算法，该方法使用一种新的物种进化策略来降低种群大小，且提出一种新的多变异流程来加速收敛。崔晓飞将改进了粒子群算法用于叠后波阻抗反演。该方法将群体智能与概率论相结合用于全局优化。S.Mahdavip综述了大尺度全局连续优化的启发式演算法，包括大尺度全局优化方法(LSGO)、进化算法[EAs]、协同进化算法[CC]等。许多差分进化算法，如协同进化算法、ND-CC及协同多任务学习算法等，均可用于解决高维优化问题。同时，大多数现有的方法均有“维数瓶颈”问题。这些方法的优化能力随着维数的增加而急剧降低。为了对地下介质进行精细刻画，我们需要将地下介质剖分成细小的网格，这样就导致模型参数数量巨大。反演高维参数将导致现有的协同进化算法失效或者需要大量迭代。

发明内容

为了解决现有技术中存在的问题，本发明的目的是提出一种方向引导的差分进化算法，用于射线追踪的角度估计。估计的角度可用于射线追踪并最终得到合成叠前地震反射记录。另一个创新点是目标函数中包含了一个符号函数用于引导变异的方向，且提出一个迭代终止准则以加速计算。

为了实现上述目的，本发明采用的技术方案是，基于方向引导的协同差分进化算法，包括以下步骤：

步骤1)：首先计算射线的初始出射角，并设定出射角的上下限；

步骤2)：在步骤1)给出的出射角的上下限范围内随机抽取角度参数，给出一个包含M个初始个体的群体P；

步骤3)：用观测射线与计算射线在地表上的终点的距离的误差能量定义第k个个体的局部适应度函数；

步骤4)：定义第k个个体的第j个分量的变异策略：

其中符号函数S定义为：

步骤5)：用步骤4)定义的变异策略得到变异后群体P的新值进行射线追踪，计算每条射线的局部适应度值，直至满足终止条件，则退出迭代；否则就返回步骤3)；

步骤6)：根据步骤5)得到的局部适应度值选择最优参数，最小的局部适应度值所对应的参数即为最优参数。

进一步的，所述步骤5)中，局部适应度函数定义为观测射线与计算射线在地表上的终点的距离的误差能量。

进一步的，所述步骤4)中，F是控制因子，其值为F＝F₁·2-j,j＝L,L-1,L-2,…,1,0，其中，F₁是分布于0～1之间的随机数，并且当误差能量Ea_j的趋势是上升时，F₁将乘以增益以扩大变异值的搜索边界，否则，无增益；L是一个常数，其值由所要求的精度决定，L越大，取值间隔2^-j越小，则精度越高。

进一步的，所述步骤4)中，变异策略中引入S函数，S函数的值的正负特征决定变异的方向。

进一步的，角度初值调整为：

其中，h为地层深度，x_s为炮点到中点的距离，x_r为中点到检波点的距离，θ_s为炮点与反射点之间射线的出射角，θ_r为检波点与反射点之间射线的出射角。

进一步的，所述步骤5)中局部适应度值迭代计算的终止条件有两个，满足任意终止条件时即终止迭代：一是硬阈值条件，即迭代次数达到给定的阈值即终止迭代；二是软阈值条件，若误差能量满足Ea_j＜E_min则终止迭代，其中，Ea_j是合成地震记录的误差能量，E_min给定的误差精度。

进一步的，所述步骤1)中出射角的上下限范围为±20°。

本发明还提供了基于方向引导的协同差分进化算法在射线追踪中的应用，采用本发明提供的协同差分进化算法，将步骤6)得到的最优参数用于射线追踪，用最优参数根据射线追踪方法求取合成地震资料。

与现有技术相比，本发明至少具有以下有益效果，本发明在变异策略中引入了一个符号函数来定义变异的方向以提高算法收敛性，且提出两种灵活有效的迭代终止准则以改进计算效率。不同于局部优化波形反演算法和其他协同进化算法，本发明提出的方法不依赖于初值，对高维可分模型空间有较好的可行性，计算速度快，收敛性好，对维度较大的参数依然能有效反演；并且，在变异策略中引入S函数，S函数的值的正负特征决定变异的方向；引入S函数使得变异沿着缩小误差能量的方向进行，使算法加速收敛，能快速收敛到全局最优解。

进一步的，所述步骤4)中，控制因子F的值是非均匀的，因此，在最小局部适应度值对应的变量附近有较多的变异值，而在远离最小值的地方变异后的值较少，根据视网膜非均匀感受野机理，本发明采取的这种策略可以更好地分配计算资源、改善计算效率，类似于视感细胞和中间细胞的多对一连接模式，获取最优参数，只要得到一个最小的局部适应度值，就可以在如图3中圆环所示的范围内搜索到最优参数。

进一步的，本发明中是在初始角度上叠加±20°，±20°是根据经验选定的，实际中可以根据不同的观测系统选择不同的数值，如果边界太窄，可能无法得到全局最优解，DE算法仅使用初值的低频分量而不是使用初值本身，因此进化算法不依赖于初值。

进一步的，根据认知计算的交互学习模式，通过对本发明提出的基于方向引导的协同差分进化算法(DG-DECCM)方法的错误控制功能的研究，该功能基于交互原理，如果误差能量Ea_j的趋势是上升的，F₁将乘以一个增益以扩大变异值的搜索边界，否则，无增益。

在实际中，若误差能量满足不等式，则终止迭代；若始终无法满足不等式，为了节省计算资源避免死循环，则根据迭代次数选择硬阈值条件退出循环。

附图说明

图1(a)为叠前CSP资料的观测系统；图1(b)为初始出射角及射线终点示意图。

图2为控制因子F示意图。

图3为计算最优参数示意图。

图4为符号函数S定义示意图。

图5为共中心点道集(CMP)的观测系统。

图6(a)为五层介质速度模型；图6(b)为出射角边界范围。

图7(a)第一层介质中的射线追踪结果；图7(b)为五层介质中的射线追踪结果。

图8(a)误差能量随时间的收敛曲线；图8(b)为合成CMP资料。

图9(a)和图9(b)为用局部最优算法自适应角度步长射线追踪方法得到的误差能量的收敛曲线；图9(c)误差能量函数示意图。

图10(a)为最短路径法的射线追踪结果；图10(b)为试射法的射线追踪结果。

图11(a)为基于CSP数据5层介质的射线追踪结果；图11(b)为第4层介质的射线追踪结果。

图12(a)为高速体模型；图12(b)为高速体模型的射线追踪结果。

图13(a)为低速体模型及其射线追踪结果；图13(b)为误差能量的收敛曲线。

图14为DG-DECCM算法的流程示意图。

具体实施方式

下面结合附图和具体实施方式对本发明做进一步说明。

地震反射波法操作简单，相较于VSP及井间地震来说花费较小。其检波器位于地表，这样就减小了钻井的花费。相对于叠后反射数据，叠前地震反射资料不受正常时差(NMO)的影响，有精确的层位信息，丰富的频率和走时信息。我们可以从叠前数据中观察到细小的地层特征如薄互层中的波形变化等。叠前共炮点道集的观测系统如图1(a)所示。震源位于地表炮点处，检波器沿着地表的测线排列。正演是反问题的关键步骤。我们研究基于叠前地震数据的射线追踪正演模拟方法。由于射线方程计算的射线无法反射，我们将一条完整的射线分为炮点到反射点及检波点到反射点两部分，如图1(b)所示。p_x和p_z分别表示沿x轴和z轴的慢度，即速度倒数的分量。T_s和T_r分别表示两条射线从反射点以出射角θ_s和θ_r射出后到达地表所用的走时。二维射线方程为：

其中，t为走时，v为某空间点的速度，X＝(x,z)为射线的坐标。射线方程的积分变量是时间，微分方程(1)的数值解是每个时间采样点的位置和梯度方向。因此，射线的端点位置是出射角的函数。如果反射点和出射角θ_s及θ_r已知，射线方程可用四阶龙格库塔方法求解。我们用提出的DG-DECCM方法计算出射角。

DE算法模拟自然界中物种进化的过程并寻找最优解。认知科学可用来最小化进化过程。首先，我们给出一个初始种群，然后对它的某些个体进行变异。我们提出的算法可以自主学习，独立思考并判断该变异是否为进化的最优方向。最后，通过迭代循环使得所有个体成为最优。基于进化算法得到的结果，我们可以找到最优的参数并作出快速和有效的决定。提出的DE算法首先给出一个包含M个初始个体的群体P：

P＝[p₁,p₂,…,p_M] (2)

第k个个体的表达式如下：

p_k＝[p_k,1,p_k,2,…,p_k,N],k＝1,…,M (3)

其中，N是估计的参数的数量。对于包含N/2个地震道的叠前地震记录，有N个待估角度，所以待估参数的维度是N。局部适应度函数由观测数据和计算数据的误差能量定义，可以用来控制变异方向。第k个个体的局部适应度函数定义为：

其中，[x_o(j),z_o(j)]是第j条观测射线的终点，[x_c(j),z_c(j)]是第j条计算的射线的终点。x和z是水平和垂直坐标轴。局部适应度函数是一个M×N的矩阵，

一个好的变异策略能有效地修正个体的基因(即角度)，几次迭代后即可得到最优解。一个群体包含M个个体，M可以取任意值，但若M太小，将不能保证得到的是全局最优解；若M太大，运算时间将会很长。实际中，我们取M为100。如果有N个待估参数，个体基因数将会是N。每个群体的M个初始个体需进行变异和更新，最终我们选择最优个体作为所需。

选择变量p_r,j作为排序后的最小的局部适应度函数值对应的参数，再随机选两个不同的变量p_r1,j和p_r2,j，则第k个个体的第j个分量的变异策略定义为：

其中，F是控制因子，其值为F＝F₁·2^-j,j＝L,L-1,L-2,…,1,0，F₁是分布于0～1之间的随机数。L是一个常数，其值由所要求的精度决定。如图2所示，F的值是非均匀的，因此，在最小局部适应度值对应的变量附近有较多的变异值，而在远离最小值的地方变异后的值较少。根据视网膜非均匀感受野机理，我们采取的这种策略可以更好地分配计算资源、改善计算效率。类似于视感细胞和中间细胞的多对一连接模式，获取最优参数的示意图如图3所示。只要我们得到一个最小的局部适应度值，我们就可以在如图中圆环所示的范围内搜索到最优参数。S是符号函数，定义为：

其中，x_o和x_c是观测射线和计算所得射线的终点的水平坐标。图4是定义符号函数的示意图。若计算射线的终点的横坐标x_c1大于观测终点的横坐标，则出射角θ₁应增大。此时，最优参数可能在变量p_r,j的右侧，所以S是正的。反之，若x_c2在xo的左边，则出射角θ₂应减小，因此S是负的。变异后的基因将在最优值附近扰动，这种方向选择性将加速进化速度。

入图14所示，DG-DECCM算法流程图，其中，gen_max为给定的硬阈值。我们设置了两种迭代终止条件，一是硬阈值法，即如果迭代次数达到给定的的阈值则程序退出循环。另一种是软阈值法：

Ea_j＜E_min (7)

其中，

是平均误差能量，E_min是给定的最小精度阈值。若满足式(7)的不等式，则程序退出循环。根据认知计算的交互学习模式，我们研究了DG-DECCM方法的错误控制功能。该功能基于交互原理，如果误差能量Ea_j的趋势是上升的，F₁将乘以一个增益以扩大变异值的搜索边界，否则，无增益。

从流程图中我们可以看出，待估的可分模型参数(即角度)有N个，最终可以得到N条射线。所有的子分量在局部适应度函数的引导下，将同时独立地进行变异和选择，成为新的个体并与原始的个体随机交叉。出射角的初值根据几何关系来定义，其示意图如图1(b)所示：

其中，x_s是炮点和中点的距离，x_r是接收点和中点的距离，h是反射点的深度。射线的出射角定义为z轴与射线的夹角，因此我们将角度初值调整为：

算例分析

1、叠前共中心点道集(CMP)模型

我们将提出的方法用于CMP数据，该数据基于水平层状速度模型得到。观测系统如图5所示。五层介质的速度分别是2000m/s，2200m/s，2400m/s，2600m/s和2800m/s，如图6(a)所示。深度分别为200m，400m，600m，800m和1000m。震源位于地表距离原点240m处，以5m为间隔。检波器以5m为间隔位于地表260m到500m处。

群体中有100个初始个体，每一个个体的变量值均随机产生。为了快速有效地最优化个体，我们定义的搜索区间的上下限。角度初值根据几何关系确定。如图6(b)所示，上下限的定义方法是在初始角度上叠加±20°，±20°是根据经验选定的，实际中可以根据不同的观测系统选择不同的数值。如果边界太窄，可能无法得到全局最优解。DE算法仅使用初值的低频分量而不是使用初值本身，因此进化算法不依赖于初值。

第一层中的射线如图7(a)所示。从图中可以看出，射线从震源出发进入均匀介质，以直线传播。5层正演模拟射线如图7(b)所示，当射线在速度随深度变化的介质中传播时，其路径为曲线。图8(a)是误差能量的收敛曲线。误差能量定义为观测曲线和计算曲线的终点的距离。有49个检波器，所以总射线数为49×2×5＝490个。50次迭代后，误差小于50m，这意味着每条射线的平均误差能量是0.1m左右，这在实际中是可以接受的。图8(b)是基于计算射线的合成单炮集记录。

图9(a)和(b)是使用角度步长射线追踪算法求取的误差能量收敛曲线，该算法是一种局部优化算法。从图中可以看出，误差能量以常振幅的形式振动，这是因为误差能量是如图9(c)所示的出射角的二次函数，寻找最优出射角就是寻找误差能量曲线的“碗底”的过程。但是当图9(c)中虚线的两端离最小值点的距离相等时，误差能量无法收敛到最小值点且会在虚线两端不断震荡。在个人PC机上该方法得到490条射线需运行6个小时。实际上，每条射线的追踪时间从0.1s到80s不等。实验结果表明，越深层的同相轴花费的射线追踪时间越多。这是因为地层越深，地下介质的结构越复杂，因此误差能量越难收敛。因此，局部优化方法对水平层状介质非常有效，但不适用于复杂地下地质情况。

图10(a)和(b)是由最短路径法和试射法得到的射线路径。最短路径法遍历所有节点，然后反向追踪最合适的射线路径。试射法逐步调整初始出射角。这两种方法原理简单，但耗时长且效率低。

2、叠前共炮点道集(CSP)模型

我们用一个CSP模型验证了提出方法的有效性。速度模型与图4所示相同，观测系统如图1(a)所示。深度为200m，400m，600m，800m和1000m，震源位于地表原点处，检波器以10m为间隔位于20m到500m的测线上。图11(a)是5层模型的射线追踪结果，图11(b)是第4层的射线追踪结果。该实验的结论与第1小节类似：在均匀速度模型中射线是直线，在水平层状速度模型中射线是折线。

3、高速体模型

我们用一个高速体模型验证了提出的方法的有效性。速度模型如图12(a)所示，速度定义为：v＝2+0.4sin(0.25πz)sin(0.5πx)km/s。高速体的速度存在连续的变化，反射点位于中点(1000,2000)m处。震源从原点到960m以20m为间隔进行排列，接收器以同样的间隔位于地表1060m到2000m的测线上。图12(b)是射线追踪结果，其射线向异常高速区域外侧弯曲。

4、低速体模型

最后一个例子是低速体模型。速度模型如图13(a)所示，其速度定义为：v＝2-0.4sin(0.25πz)sin(0.5πx)km/s。低速体中也存在连续变化的边界，反射点位于中点(1000,2000)m处，震源和检波点的分布与第3节相同。图13(a)中，黑线是射线追踪的结果，其射线弯向异常低速体内部，图13(b)是单调递减的收敛曲线。

Claims (5)

1.基于方向引导的协同差分进化算法，其特征在于，包括以下步骤：

步骤1)：首先计算射线的初始出射角，并设定出射角的上下限；

步骤2)：在步骤1)给出的出射角的上下限范围内随机抽取角度参数，给出一个包含M个初始个体的群体P；

步骤3)：用观测射线与计算射线在地表上的终点的距离的误差能量定义第k个个体的局部适应度函数；

步骤4)：定义第k个个体的第j个分量的变异策略：

其中符号函数S定义为：

步骤5)：用步骤4)定义的变异策略得到变异后群体P的新值进行射线追踪，计算每条射线的局部适应度值，直至满足终止条件，则退出迭代；否则就返回步骤3)；

步骤6)：根据步骤5)得到的局部适应度值选择最优参数，最小的局部适应度值所对应的参数即为最优参数；

所述步骤4)中，F是控制因子，其值为F＝F₁·2^-j,j＝L,L-1,L-2,…,1,0，其中，F₁是分布于0～1之间的随机数，并且当误差能量Ea_j的趋势是上升时，F₁将乘以增益以扩大变异值的搜索边界，否则，无增益；L是一个常数，其值由所要求的精度决定，L越大，取值间隔2^-j越小，则精度越高；

所述步骤4)中，变异策略中引入S函数，S函数的值的正负特征决定变异的方向；

所述步骤5)中局部适应度值迭代计算的终止条件为，一是硬阈值条件，即迭代次数达到给定的阈值即终止迭代；二是软阈值条件，若误差能量满足Ea_j＜E_min则终止迭代，其中，Ea_j是合成地震记录的误差能量，E_min给定的误差精度。

2.根据权利要求1所述的基于方向引导的协同差分进化算法，其特征在于，所述步骤5)中，局部适应度函数定义为观测射线与计算射线在地表上的终点的距离的误差能量。

3.根据权利要求1所述的基于方向引导的协同差分进化算法，其特征在于，角度初值调整为：

h为地层深度，x_s为炮点到中点的距离，x_r为中点到检波点的距离，θ_s为炮点与反射点之间射线的出射角，θ_r为检波点与反射点之间射线的出射角。

4.根据权利要求1所述的基于方向引导的协同差分进化算法，其特征在于，所述步骤1)中出射角的上下限范围为±20°。

5.基于方向引导的协同差分进化算法在射线追踪中的应用，其特征在于，采用权利要求1至4中任一项所述的协同差分进化算法，将步骤6)得到的最优参数用于射线追踪，用最优参数根据射线追踪方法求取合成地震资料。

Images