CN103116703B

CN103116703B - A kind of covariation differential evolution algorithm for higher-dimension parameter space waveform inversion

Info

Publication number: CN103116703B
Application number: CN201310044383.8A
Authority: CN
Inventors: 高静怀; 汪超; 王大兴
Original assignee: Xian Jiaotong University
Current assignee: Xian Jiaotong University
Priority date: 2013-02-04
Filing date: 2013-02-04
Publication date: 2016-03-30
Anticipated expiration: 2033-02-04
Also published as: CN103116703A

Abstract

The invention discloses a collaborative variation differential evolution algorithm for high-dimensional parameter space waveform inversion. The algorithm introduces the idea of decomposition-coordination into the differential evolution algorithm, decomposes high-dimensional individuals into a series of subcomponents, and introduces The local fitness function evaluates each subcomponent. Then in the mutation operation, the local fitness is used to guide the mutation direction of each sub-component, and in the selection operation, the global fitness is used to coordinate the sub-components to achieve co-evolution. Compared with the commonly used fast simulated annealing method and genetic algorithm, the co-variation differential evolution algorithm is more suitable for high-dimensional parameter space waveform inversion; when the layer is thin and there are many parameters to be inverted, the co-variation differential evolution algorithm can search to a solution closer to the true value. In addition, the convergence speed of the co-variation differential evolution algorithm is not sensitive to the increase of dimension, and it still uses a fast convergence speed when the dimension is high.

Description

A Co-variational Differential Evolution Algorithm for Waveform Inversion in High-Dimensional Parameter Spaces

技术领域technical field

本发明属于地球物理勘探领域，涉及一种地球模型物性参数反演方法，特别涉及一种用于高维参数空间波形反演的协同变异差分进化算法。The invention belongs to the field of geophysical exploration, and relates to an inversion method for physical property parameters of an earth model, in particular to a collaborative variation differential evolution algorithm for inversion of high-dimensional parameter space waveforms.

背景技术Background technique

随着油气田勘探开发程度的加深，油气勘探对象日趋复杂，逐渐由简单构造油气藏向复杂构造油气藏转移，从构造油气藏向地层－岩性等隐蔽油气藏转移。如何精细描述复杂构造和地层岩性是当前我国油气勘探面临的关键问题。波形反演方法直接应用地震波形反演地下介质参数，可包含各种类型的波，如多次波、转换波、面波等，充分地利用了地震波的旅行时、振幅、相位和频率等信息。所以，波形反演是准确了解地下复杂结构和介质物理性质（如速度、密度等）的潜在有效方法，近年来在勘探地球物理领域备受关注。With the deepening of exploration and development of oil and gas fields, the objects of oil and gas exploration are becoming more and more complex, gradually shifting from simple structural oil and gas reservoirs to complex structural oil and gas reservoirs, and from structural oil and gas reservoirs to stratigraphic-lithological and other subtle oil and gas reservoirs. How to finely describe complex structures and formation lithology is a key issue facing my country's current oil and gas exploration. The waveform inversion method directly uses the seismic waveform to invert the parameters of the underground medium, which can include various types of waves, such as multiple waves, converted waves, surface waves, etc., and fully utilizes the travel time, amplitude, phase and frequency of seismic waves. . Therefore, waveform inversion is a potentially effective method to accurately understand the complex subsurface structure and physical properties of media (such as velocity, density, etc.), and has attracted much attention in the field of exploration geophysics in recent years.

已有的波形反演方法大体上可分为两大类：一类是基于梯度的局部优化方法，如梯度法、牛顿法和共轭梯度法等，这类方法容易陷入局部极值，结果依赖于初始模型；另一类是全局优化方法，如模拟退火法、邻域法、遗传算法^[4]、粒子群优化法等，这类方法不依赖于初始模型和函数梯度值，但计算量较大。随着计算机计算能力的提高，全局优化方法得到越来越多的应用，但计算时间仍是制约全局优化方法在波形反演中普遍应用的一个因素。另一方面，现有的全局优化方法遭受到“维数瓶颈”问题，这些算法用于低维（维数100以下）问题时具有卓越的寻优能力，但它们的搜索能力随着搜索空间维数的增加会急剧下降。要精细刻画地下介质，必需对介质进行细小的划分，这就使得模型参数非常庞大，反问题的维数高达几百几千，常用的全局优化方法不再有效或需要大量的迭代。The existing waveform inversion methods can be roughly divided into two categories: one is gradient-based local optimization methods, such as gradient method, Newton method, and conjugate gradient method. The other is the global optimization method, such as simulated annealing method, neighborhood method, genetic algorithm ^[4] , particle swarm optimization method, etc. This kind of method does not depend on the initial model and function gradient value, but the calculation amount is relatively large. Big. With the improvement of computer computing power, the global optimization method has been used more and more, but the calculation time is still a factor that restricts the general application of the global optimization method in waveform inversion. On the other hand, the existing global optimization methods suffer from the "dimensionality bottleneck" problem. These algorithms have excellent optimization capabilities for low-dimensional (dimensions below 100) problems, but their search capabilities increase with the search space dimension. The increase in the number will drop sharply. In order to describe the subsurface medium finely, it is necessary to divide the medium finely, which makes the model parameters very large, and the dimension of the inverse problem is as high as hundreds or thousands. The commonly used global optimization methods are no longer effective or require a large number of iterations.

协同进化算法是一种新的进化算法框架，该算法采用分解——协调的思想将高维问题分解为一系列相互关联的子问题，每个子问题轮流交替进行优化。协同进化算法对高维可分问题较为有效，但不适用于强非线性的高维波形反演。差分进化（DE）算法是由Store和Price于1996年首先提出的基于实数编码的连续空间全局优化算法，是一种原理简单易于实现而又能力强大的模拟进化算法，之后，一些改进的DE算法被相继提出。差分进化算法成功地应用到很多领域，如信号处理和地球物理反演^[。但差分进化算法同样受到“维数瓶颈”问题，在高维参数空间波形反演应用中仍需改进。Coevolutionary Algorithm is a new evolutionary algorithm framework, which uses the idea of decomposition-coordination to decompose a high-dimensional problem into a series of interrelated sub-problems, and each sub-problem is optimized in turn. The coevolutionary algorithm is more effective for high-dimensional separable problems, but it is not suitable for high-dimensional waveform inversion with strong nonlinearity. The differential evolution (DE) algorithm is a continuous space global optimization algorithm based on real number coding first proposed by Store and Price in 1996. It is a simple and easy-to-implement and powerful simulated evolutionary algorithm. were proposed successively. Differential evolution algorithm has been successfully applied to many fields, such as signal processing and geophysical inversion ^[ . However, the differential evolution algorithm also suffers from the "dimensionality bottleneck" problem, and it still needs to be improved in the application of high-dimensional parameter space waveform inversion.

发明内容Contents of the invention

本发明的目的在于克服上述现有技术的缺点，提供一种用于高维参数空间波形反演的协同变异差分进化算法，该算法利用局部适应度和全局适应度同时引导算法的搜索方向，从而大大提高了收敛速度和全局寻优能力。The purpose of the present invention is to overcome the shortcoming of above-mentioned prior art, provide a kind of collaborative variation differential evolution algorithm that is used for high-dimensional parameter space waveform inversion, this algorithm utilizes local fitness and global fitness to guide the search direction of algorithm simultaneously, thereby The convergence speed and global optimization ability are greatly improved.

本发明的目的是通过以下技术方案来解决的：The purpose of the present invention is solved by the following technical solutions:

该种用于高维参数空间波形反演的协同变异差分进化算法，包括以下步骤：The collaborative mutation differential evolution algorithm for high-dimensional parameter space waveform inversion includes the following steps:

1）采集原始地震资料，然后对采集到的地震资料进行预处理，处理后得到叠前共中心点道集或叠后地震数据，称这个地震数据为观测地震数据，记为Seis(r,t)，其中r表示检波器接收点位置，t表示时间轴；1) Collect the original seismic data, then preprocess the collected seismic data, and obtain the pre-stack common midpoint gather or post-stack seismic data after processing, which is called the observed seismic data, denoted as Seis(r,t ), where r represents the receiving point position of the geophone, and t represents the time axis;

2）构建水平层状地质模型，给定地质模型的层数N和层的厚度，每层的介质模型参数包纵波波速度V_p、横波波速度V_s和密度ρ；2) Construct a horizontal layered geological model, given the number of layers N and layer thickness of the geological model, and the medium model parameters of each layer include P-wave velocity V _p , S-wave velocity V _s and density ρ;

3）确定地质模型参数V_p、V_s、ρ的搜索空间，并指定待优化的目标函数；3) Determine the search space of geological model parameters V _p , V _s , and ρ, and specify the objective function to be optimized;

4）在地质模型参数搜索空间内随机地生成NP个地质模型，并进行实数编码得到含NP个个体的初始群体；4) Randomly generate NP geological models in the geological model parameter search space, and encode them with real numbers to obtain NP individuals the initial group;

5）计算群体中每个随机模型的合成地震数据，然后根据目标函数估计各随机模型对应个体的全局适应度值F_i；5) Calculate the synthetic seismic data of each stochastic model in the population, and then estimate the global fitness value F _i of each stochastic model corresponding to the individual according to the objective function;

6）进行变异操作，对群体中每个个体产生一个变异个体首先将高维个体按照基因间相互依赖的强弱程度分解为一系列子成分，给每个子成分指定一个局部适应度函数，根据局部适应度函数估计各子成分的局部适应度值；然后随机选出三个不同个体，根据其相应子成分的局部适应度值变化情况得到梯度信息，并以此梯度的负方向作为各子成分的变异方向；6) Perform a mutation operation to generate a mutant individual for each individual in the group First, high-dimensional individuals are decomposed into a series of subcomponents according to the degree of interdependence between genes, and a local fitness function is assigned to each subcomponent, and the local fitness value of each subcomponent is estimated according to the local fitness function; Three different individuals are selected, and the gradient information is obtained according to the change of the local fitness value of the corresponding subcomponent, and the negative direction of the gradient is used as the variation direction of each subcomponent;

7）进行交叉操作，随机地从变异个体和与其相应的当代个体中抽取基因，组合成一个试验个体 7) Perform a crossover operation, randomly from the mutant individual and corresponding contemporary individuals Extract genes from the group and combine them into a test individual

8）进行选择操作，根据当代个体和试验个体的全局适应度值选取下一代；8) Perform a selection operation, and select the next generation according to the global fitness value of the current individual and the test individual;

9）进化代数g=g+1，判断是否满足终止条件，如果进化代数g小于或等于设定的值G，则返回到步骤6）；否则执行步骤10）；9) Evolutionary algebra g=g+1, judge whether the termination condition is satisfied, if the evolutionary algebra g is less than or equal to the set value G, return to step 6); otherwise, execute step 10);

10）选出第G代群体中全局适应度值最小的个体，将该个体解码后即得到最终搜索到的最优地质模型。10) Select the individual with the smallest global fitness value in the G-th generation group, and decode the individual to obtain the optimal geological model that is finally searched.

进一步，上述步骤3）中，参数V_p、V_s、ρ的搜索空间是根据测井资料中相应参数的测井记录确定，即以测井得到的地质模型的低频分量为背景，背景值增加设定值后为搜索空间上界，背景值减小设定值后为搜索空间下界；或者纵波波速度V_p和横波波速度V_s的搜索空间也能够根据高精度速度分析确定；所述目标函数是刻画最优解的标准，以计算波场与观测波场之间的拟合程度或误差大小为标准。Further, in the above step 3), the search space of the parameters V _p , V _s , and ρ is determined according to the logging records of the corresponding parameters in the logging data, that is, the low-frequency component of the geological model obtained from the logging is used as the background, and the background value increases The upper bound of the search space is after the set value, and the lower bound of the search space is after the background value decreases the set value; or the search space of the longitudinal wave velocity V _p and the shear wave velocity V _s can also be determined according to high-precision velocity analysis; the target The function is the standard for describing the optimal solution, taking the degree of fit or the size of the error between the calculated wave field and the observed wave field as the standard.

上述步骤5）中，以观测地震数据与由模型合成的计算地震数据之间误差绝对之和为目标函数，即目标函数为：In the above step 5), the absolute sum of the errors between the observed seismic data and the calculated seismic data synthesized by the model is used as the objective function, that is, the objective function is:

${F f}_{i i} = = {Σ Σ}_{r r = = 11}^{NR NR} &Integral; &Integral; | | Seis Seis ((r r,, t t)) - - {Syn Syn}_{i i} ((r r,, t t)) | | dt dt . .$

式中NR为检波器个数，适应度值F_i越小，表示观测地震数据与计算地震数据之间误差越小，则其相应的个体就越优秀；反之，适应度值F_i越大，则表示其相应的个体就越差。where NR is the number of geophones, and the smaller the fitness value F _i is, the smaller the error between the observed seismic data and the calculated seismic data is, and the corresponding individual The more excellent; on the contrary, the larger the fitness value F _i is, it means that the corresponding individual worse.

上述步骤6）中，个体的变异方向是由子成分的局部适应度值来引导：随机选出三个不同个体，将子成分的基因由小到大排列，如果其对应的局部适应度值是递增的case1，则在最小基因的左方寻找新基因；如果其对应的局部适应度值是递减的case2，则在最大基因的右方寻找新基因；如果其对应的局部适应度值是先减后增case3，则在中间基因的附近寻找新基因；如果其对应的局部适应度值是先增后减case4，则在最小基因的左方或在最大基因的右方寻找新基因：In the above step 6), the variation direction of the individual is guided by the local fitness value of the subcomponent: randomly select three different individuals, and arrange the genes of the subcomponent from small to large, if the corresponding local fitness value is increasing For case1, the minimum gene Find a new gene on the left of ; if its corresponding local fitness value is decreasing case2, then the maximum Find new genes on the right side of Find a new gene near ; if its corresponding local fitness value increases first and then decreases case4, then find a new gene on the left of the smallest gene or on the right of the largest gene:

${v v}_{i i,, j j}^{g g} = = \{\begin{matrix} {x x}_{{r r}_{1111},, j j}^{g g} - - α α | | {x x}_{{r r}_{22},, j j}^{g g} - - {x x}_{{r r}_{33},, j j}^{g g} | |,, & case case 11 \\ {x x}_{{r r}_{3131},, j j}^{g g} + + α α | | {x x}_{{r r}_{22},, j j}^{g g} - - {x x}_{{r r}_{33},, j j}^{g g} | |,, & case case 22 \\ {x x}_{{r r}_{21 twenty one},, j j}^{g g} + + α α (({x x}_{{r r}_{22},, j j}^{g g} - - {x x}_{{r r}_{33},, j j}^{g g})),, & case case 33 \\ {p p}_{11} {x x}_{{r r}_{1111},, j j}^{g g} + + {p p}_{22} {x x}_{{r r}_{21 twenty one},, j j}^{g g} + + {p p}_{33} {x x}_{{r r}_{3131},, j j}^{g g} & case case 44 \end{matrix} j j = = 1,2 1,2,, . . . . . .,, J J$

${p p}_{11} = = {p p}_{44} / / {LF LF}_{j j} (({x x}_{{r r}_{1111},, j j}^{g g})),, {p p}_{22} = = {p p}_{44} / / {LF LF}_{j j} (({x x}_{{r r}_{21 twenty one},, j j}^{g g})),, {p p}_{33} = = {p p}_{44} / / {LF LF}_{j j} (({x x}_{{r r}_{3131},, j j}^{g g})),,$

${p p}_{44} = = 11 / / ((L L {F f}_{j j} {(({x x}_{{r r}_{1111},, j j}^{g g}))}^{- - 11} + + {{LF LF}_{j j} (({x x}_{{r r}_{22} 11,, j j}^{g g}))}^{- - 11} + + {LF LF}_{j j} {(({x x}_{{r r}_{3131},, j j}^{g g}))}^{- - 11})) . .$

其中α为变异尺度因子，一般选取0到1之间的实数，J为子成分的个数，LF_j为各基因对应的局部适应度值。Among them, α is the variation scale factor, generally a real number between 0 and 1 is selected, J is the number of subcomponents, and LF _j is the local fitness value corresponding to each gene.

进一步的，在上述步骤7）中，Further, in the above step 7),

${u u}_{i i,, n no}^{g g} = = \{\begin{matrix} {v v}_{i i,, n no}^{g g} & ifrand ifrand \leq \leq CRorn CRorn = = {n no}_{rand rand},, \\ {x x}_{i i,, n no}^{g g} & otherwise otherwise,, \end{matrix} n no = = 1,2 1,2,, . . . . . .,, D D.$

其中CR∈[0,1)是交叉率；rand为[0，1]均匀分布随机数；n_rand是1到D之间的一随机整数。Among them, CR∈[0,1) is the crossover rate; rand is a uniformly distributed random number in [0,1]; n _rand is a random integer between 1 and D.

上述步骤8）中，比较当代个体的全局适应度值和试验个体的全局适应度值，全局适应度值小的个体被选入下一代，具体操作为：In the above step 8), the global fitness value of the current individual is compared with the global fitness value of the test individual, and the individual with a small global fitness value is selected into the next generation. The specific operation is as follows:

${x x}_{i i}^{g g + + 11} = = \{\begin{matrix} {u u}_{i i}^{g g} & ifF ifF (({u u}_{i i}^{g g})) < < F f (({x x}_{i i}^{g g})),, \\ {x x}_{i i}^{G G} & otherwise otherwise,, \end{matrix} i i = = 1,2 1,2,, . . . . . .,, NP NP . . . .$

本发明具有以下有益效果：The present invention has the following beneficial effects:

本发明将分解—协调的思想引入差分进化算，将高维个体分解为一系列的子成分，并引入局部适应度函数对每个子成分进行评价。然后在变异操作中利用局部适应度引导各子成分的变异方向，而在选择操作中利用全局适应度协调各子成分达到共同进化。与常用的快速模拟退火法和遗传算法相比，协同变异差分进化算法更加适合于高维参数空间波形反演；当层较薄、待反演参数多的情况下，协同变异差分进化算法能搜索到更加接近真实值的解。另外，协同变异差分进化算法的收敛速度对维数增加不敏感，维数很高时仍用很快的收敛速度。The invention introduces the idea of decomposition-coordination into the differential evolution algorithm, decomposes the high-dimensional individual into a series of subcomponents, and introduces a local fitness function to evaluate each subcomponent. Then in the mutation operation, the local fitness is used to guide the mutation direction of each sub-component, and in the selection operation, the global fitness is used to coordinate the sub-components to achieve co-evolution. Compared with the commonly used fast simulated annealing method and genetic algorithm, the co-variation differential evolution algorithm is more suitable for high-dimensional parameter space waveform inversion; when the layer is thin and there are many parameters to be inverted, the co-variation differential evolution algorithm can search to a solution closer to the true value. In addition, the convergence speed of the co-variation differential evolution algorithm is not sensitive to the increase of dimension, and it still uses a fast convergence speed when the dimension is high.

附图说明Description of drawings

图1是本发明流程示意图；Fig. 1 is a schematic flow chart of the present invention;

图2是本发明对高维个体的分解和局部适应度函数示意图；Fig. 2 is a schematic diagram of decomposition and local fitness functions of high-dimensional individuals in the present invention;

图3是本发明中层状介质的分解示意图和计算局部适应度函数的时窗示意图。Fig. 3 is a schematic diagram of decomposition of a layered medium and a schematic diagram of a time window for calculating a local fitness function in the present invention.

具体实施方式detailed description

下面结合附图对本发明做进一步详细描述：The present invention is described in further detail below in conjunction with accompanying drawing:

地震波形反演的目的是获取一个最优的地层模型，使由该地层模型计算的地震数据与实测地震数据拟合最佳。本发明采用一种改进的差分进化算法来搜索最优地层模型参数。基于差分进化算法框架下，借鉴协同进化算法中分解——协调的思想，把高维问题按参数相互依赖程度适当地分解为一系列子问题，并引入局部适应度函数对各子问题进行评价，利用局部适应度和全局适应度同时引导算法的搜索方向。The purpose of seismic waveform inversion is to obtain an optimal stratigraphic model, so that the seismic data calculated by the stratigraphic model can best fit the measured seismic data. The invention adopts an improved differential evolution algorithm to search for optimal formation model parameters. Based on the framework of differential evolution algorithm, referring to the idea of decomposition-coordination in co-evolution algorithm, the high-dimensional problem is properly decomposed into a series of sub-problems according to the degree of interdependence of parameters, and a local fitness function is introduced to evaluate each sub-problem. The search direction of the algorithm is guided simultaneously by local fitness and global fitness.

本发明的物质基础是通过野外高分辨率地震采集设备采集到的大量地震数据。本发明的基于改进差分进化算法的波形反演框架如图1所示，具体步骤分别为：The material basis of the invention is a large amount of seismic data collected by field high-resolution seismic acquisition equipment. The waveform inversion framework based on the improved differential evolution algorithm of the present invention is shown in Figure 1, and the specific steps are respectively:

1）采集原始地震资料，然后对采集到的地震资料做常规预处理，包括静校正处理、噪音压制处理、真振幅恢复处理等。处理后得到叠前共中心点道集或叠后地震数据，称这个地震数据为观测地震数据，记为Seis(r,t)，其中r表示检波器接收点位置，t表示时间轴。1) Collect the original seismic data, and then perform conventional preprocessing on the collected seismic data, including static correction processing, noise suppression processing, true amplitude restoration processing, etc. After processing, pre-stack common center point gathers or post-stack seismic data are obtained, which are called observed seismic data and denoted as Seis(r, t), where r represents the location of the receiving point of the geophone, and t represents the time axis.

2）构建水平层状地质模型，给定地质模型的层数N和层的厚度，每层的介质模型参数包含V_p、V_s、ρ。这里，V_p、V_s、ρ分别为纵波波速度、横波波速度和密度参数。2) Construct a horizontal layered geological model, given the layer number N and layer thickness of the geological model, the medium model parameters of each layer include V _p , V _s , and ρ. Here, V _p , V _s , and ρ are the longitudinal wave velocity, shear wave velocity, and density parameters, respectively.

3）确定地质模型参数V_p、V_s、ρ的搜索空间，并指定待优化的目标函数。参数搜索空间可根据测井资料中相应参数的测井记录确定，即以测井得到的地质模型的低频分量为背景，背景值增加一定值后为搜索空间上界，背景值减小一定值后为搜索空间下界。纵波波速度和横波波速度的搜索空间也可根据高精度速度分析确定。目标函数是刻画最优解的标准，一般以计算波场与观测波场之间的拟合程度或误差大小为标准。3) Determine the search space for geological model parameters V _p , V _s , and ρ, and specify the objective function to be optimized. The parameter search space can be determined according to the logging records of the corresponding parameters in the logging data, that is, the low-frequency component of the geological model obtained from the logging is used as the background, and the upper limit of the search space is set after the background value increases by a certain value; is the lower bound of the search space. The search spaces for P-wave and S-wave velocities can also be determined from high-precision velocity analysis. The objective function is the standard for describing the optimal solution, and generally takes the degree of fit or the size of the error between the calculated wave field and the observed wave field as the standard.

4）在地质模型参数搜索空间内随机地生成NP个地质模型，并进行实数编码得到含NP个个体的初始群体。生成的初始群体应尽可能地均匀分布在整个搜索空间。每个个体可用下面的向量表示：4) Randomly generate NP geological models in the geological model parameter search space, and encode them with real numbers to obtain NP individuals the initial group. The generated initial population should be as evenly distributed as possible across the search space. Each individual can be represented by the following vector:

${x x}_{i i}^{g g} = = {{{x x}_{i i,, 11}^{g g},, {x x}_{i i,, 22}^{g g},, . . . . . .,, {x x}_{i i,, D D.}^{g g}}} i i = = 1,2 1,2,, . . . . . .,, NP NP,,$

其中上标g表示繁殖代数，初始群体时g=0，下标D为个体的维数，也即待优化问题的参数个数。Among them, the superscript g represents the reproduction algebra, g=0 in the initial population, and the subscript D is the dimension of the individual, that is, the number of parameters of the problem to be optimized.

5）计算群体中每个随机模型的合成地震数据Syn_i(r,t)，然后根据目标函数估计各随机模型对应个体的全局适应度值F_i。5) Calculate each random model in the population Synthetic seismic data Syn _i (r, t) of Syn i (r, t), and then estimate the global fitness value F _i of each stochastic model corresponding to the individual according to the objective function.

本发明以观测地震数据与由模型合成的计算地震数据之间误差绝对之和为目标函数，即目标函数为：The present invention takes the absolute sum of errors between the observed seismic data and the calculated seismic data synthesized by the model as the objective function, that is, the objective function is:

式中NR为检波器个数。适应度值F_i越小，表示观测地震数据与计算地震数据之间误差越小，则其相应的个体就越优秀；反之，适应度值F_i越大，则表示其相应的个体就越差。where NR is the number of detectors. The smaller the fitness value F _i is, the smaller the error between the observed seismic data and the calculated seismic data is, and the corresponding individual The more excellent; on the contrary, the larger the fitness value F _i is, it means that the corresponding individual worse.

6）进行变异操作，对群体中每个个体产生一个变异个体首先根据分解——协调的思想将高维个体按照适当条件分解为一系列子成分，给每个子成分指定一个局部适应度函数，根据局部适应度函数估计各子成分的局部适应度值，然后利用局部适应度值引导各子成分的变异方向。6) Perform a mutation operation to generate a mutant individual for each individual in the group First, according to the idea of decomposition-coordination, the high-dimensional individual is decomposed into a series of subcomponents according to appropriate conditions, and a local fitness function is assigned to each subcomponent, and the local fitness value of each subcomponent is estimated according to the local fitness function, and then used The local fitness value guides the variation direction of each subcomponent.

目前所有进化算法均使用一个统一的适应度函数来评价个体中的所有基因（变量），这样的评价准则对高维问题不是很有效。本发明借鉴协同进化算法和人工培育优良物种的思想，将高维个体分解为一系列子成分，并引入局部适应度函数对每个子成分进行评价。为了便于区分，我们将评价整个个体的适应度函数称为全局适应度函数。如图2所示，将D维的个体分解为J个子成分，每个子成分记为为了便于描述，假定每个子成分均包含M个基因。个体的表达式可重写为：All current evolutionary algorithms use a unified fitness function to evaluate all genes (variables) in an individual, and such evaluation criteria are not very effective for high-dimensional problems. The present invention learns from the cooperative evolution algorithm and the idea of artificially cultivating excellent species, decomposes high-dimensional individuals into a series of subcomponents, and introduces a local fitness function to evaluate each subcomponent. In order to facilitate the distinction, we will evaluate the fitness function of the whole individual as the global fitness function. As shown in Figure 2, the D-dimensional individual Decomposed into J subcomponents, each subcomponent is recorded as For the convenience of description, it is assumed that each subcomponent contains M genes. Individual expressions can be rewritten as:

${x x}_{i i}^{g g} = = {{{x x}_{i i,, 11}^{g g},, {x x}_{i i,, 22}^{g g},, . . . . . .,, {x x}_{i i,, J J}^{g g}}} i i = = 1,2 1,2,, . . . . . .,, NP NP,,$

${x x}_{i i,, j j}^{g g} = = {{{x x}_{i i,, j j 11}^{g g},, {x x}_{i i,, j j 22}^{g g},, . . . . . .,, {x x}_{i i,, jM jM}^{g g}}} j j = = 1,2 1,2,, . . . . . .,, J J . .$

本发明中，对于层状介质波形反演问题，可以将连续几层的模型参数划分到同一个子成分里。如图3左栏所示，给出了将连续两层划分为一个子成分的示意图。In the present invention, for the layered medium waveform inversion problem, the model parameters of several consecutive layers can be divided into the same subcomponent. As shown in the left column of Fig. 3, a schematic diagram of dividing two consecutive layers into a subcomponent is given.

然后，为每个子成分指定一个局部适应度函数LF_j，此局部适应度函数必须满足以下两个条件：第j个局部适应度函数LF_j的变量至少包含第j个子分成所有局部适应度函数在全局最优点处必须为最优适应度值。对于波形反演，由于介质是根据旅行时分层，各层介质参数肯定与其相应旅行时处的反射地震波形有必然联系，从而我们指定每个子成分的局部适应度函数如下：Then, specify a local fitness function LF _j for each subcomponent, and this local fitness function must meet the following two conditions: the variable of the jth local fitness function LF _j contains at least the jth subcomponent All local fitness functions must have the optimal fitness value at the global optimum. For waveform inversion, since the medium is layered according to the travel time, the parameters of each layer of the medium must be related to the reflected seismic waveform at the corresponding travel time, so we specify the local fitness function of each subcomponent as follows:

$L L {F f}_{i i} = = {Σ Σ}_{r r = = 11}^{NR NR} &Integral; &Integral; | | [[Seis Seis ((r r,, t t)) - - {Syn Syn}_{i i} ((r r,, t t))]] \cdot &Center Dot; {win win}_{j j} ((t t)) | | dt dt . .$

其中win_j(t)是与第j个子成分相对应的时窗函数，如图3中间栏所示，对于叠后数据，时窗的中心应位于子成分的中心，对于叠前数据，随着偏移距的增加，时窗的中心应沿着图3右栏中红线所示的双曲线滑行（该曲线根据经典Dix公式求取）。上面定义的局部适应度函数有一很好的优点，只需一次正演即可计算一个个体的全局适应度值和所有局部适应度值。因此，局部适应度的计算不需增加波场正演次数，这对正演非常耗时的波形反演问题是非常重要的。where win _j (t) is the time window function corresponding to the jth subcomponent, as shown in the middle column of Fig. 3, for the post-stack data, the center of the time window should be located at the center of the subcomponent, for the pre-stack data, along with As the offset increases, the center of the time window should slide along the hyperbola shown by the red line in the right column of Figure 3 (the curve is calculated according to the classic Dix formula). The local fitness function defined above has a very good advantage. It only needs one forward modeling to calculate the global fitness value and all local fitness values of an individual. Therefore, the calculation of local fitness does not need to increase the number of wave field forward modeling, which is very important for the time-consuming waveform inversion problem of forward modeling.

实际进行波形反演时也可将一个层看成一个子成分，即每个子成分的参数只包含一层的速度和密度。如此分解的一个原因是使得算法易于实现。局部适应度函数利用一个时窗截取部分反射波来评价各子成分，要使局部适应度函数能正确地评价各子成分，必须选择适合的时窗宽度。地面记录的地震数据中每个界面的反射波不是单个的脉冲而是具有一定时宽的子波，且来自相邻界面的反射子波常常是相互重叠的。所以时窗宽度的最佳值应为一个子波的宽度加上子成分的厚度（层按旅行时划分）。地震子波的长度往往远大于层的厚度，这使得相邻时窗大部分是相互重叠的。因此，虽然每个子成分只包含一层，但相邻的连续几层之间是通过局部适应度函数关联起来的，相邻层之间的参数是协调进化的，并且计算局部适应度值的时窗是逐层滑动的，从而层间的关联性也是逐渐过渡的。In the actual waveform inversion, a layer can also be regarded as a sub-component, that is, the parameters of each sub-component only include the velocity and density of a layer. One reason for such a decomposition is to make the algorithm easy to implement. The local fitness function uses a time window to intercept part of the reflected wave to evaluate each sub-component. To make the local fitness function evaluate each sub-component correctly, a suitable time window width must be selected. In the seismic data recorded on the ground, the reflected wave of each interface is not a single pulse but a wavelet with a certain time width, and the reflected wavelets from adjacent interfaces often overlap each other. So the optimal value of the window width should be the width of a wavelet plus the thickness of the subcomponents (layers are divided by travel time). The length of the seismic wavelet is often much longer than the thickness of the layer, which makes most of the adjacent time windows overlap each other. Therefore, although each subcomponent contains only one layer, the adjacent consecutive layers are related by the local fitness function, the parameters between adjacent layers are coordinated evolution, and the time to calculate the local fitness value The windows slide layer by layer, so the correlation between layers is gradually transitioned.

完成了分解和指定局部适应度函数后，即可利用局部适应度来引导各子成分的变异方向。为了给目标个体产生一个变异个体，首先从其余NP-1个个体中随机地选出3个个体i,r₁,r₂,r₃是[1，NP]之间互不相等的随机整数。常用的变异策略是对个体进行整体操作，新变异策略则以子成分为单位进行操作，即一个子成分接一个子成分地生成变异个体。为了便于描述，假设每个子成分仅含一个变量，对第j个子成分，先将3个随机个体对应子成分中的变量和按升序进行排列：After completing the decomposition and specifying the local fitness function, the local fitness can be used to guide the variation direction of each sub-component. for target individuals To generate a mutant individual, first randomly select 3 individuals from the remaining NP-1 individuals i, r ₁ , r ₂ , r ₃ are random integers not equal to each other among [1, NP]. The commonly used mutation strategy is to operate on the individual as a whole, while the new mutation strategy operates on sub-components, that is, sub-components generate mutant individuals one by one. For the convenience of description, it is assumed that each subcomponent contains only one variable, and for the jth subcomponent, firstly, 3 random individuals correspond to the variable in the subcomponent and Sort in ascending order:

$[\begin{matrix} {x x}_{{r r}_{1111},, j j}^{g g} & {x x}_{{r r}_{21 twenty one},, j j}^{g g} & {x x}_{{r r}_{3131},, j j}^{g g} \end{matrix}] = = sort sort (([\begin{matrix} {x x}_{{r r}_{11},, j j}^{g g} & {x x}_{{r r}_{22},, j j}^{g g} & {x x}_{{r r}_{33},, j j}^{g g} \end{matrix}])) . .$

排序后三个变量对应的局部适应度值记为这三个局部适应度值只可能出现四种排列情况，分别为：The local fitness values corresponding to the three variables after sorting are recorded as There are only four possible arrangements of these three local fitness values, which are:

Case1： ${LF}_{j} (x_{r_{11}, j}^{g}) \leq {LF}_{j} (x_{r_{21}, j}^{g}) \leq {LF}_{j} (x_{r_{31}, j}^{g}) .$ Case1: ${LF}_{j} (x_{r_{11}, j}^{g}) \leq {LF}_{j} (x_{r_{twenty one}, j}^{g}) \leq {LF}_{j} (x_{r_{31}, j}^{g}) .$

这种情况可以用一个增函数来拟合，增函数指示出在的左边很可能存在一个更好的变量。This situation can be fitted with an increasing function that indicates that in There is likely to be a better variable to the left of .

Case2： ${LF}_{j} (x_{r_{11}, j}^{g}) &GreaterEqual; {LF}_{j} (x_{r_{21}, j}^{g}) &GreaterEqual; {LF}_{j} (x_{r_{31}, j}^{g}) .$ Case2: ${LF}_{j} (x_{r_{11}, j}^{g}) &Greater Equal; {LF}_{j} (x_{r_{twenty one}, j}^{g}) &Greater Equal; {LF}_{j} (x_{r_{31}, j}^{g}) .$

此种情况可以由一个减函数来拟合，而减函数指示出在的右边很可能存在一个更好的变量。This situation can be fitted by a subtraction function, which indicates that in There is likely to be a better variable on the right-hand side of .

Case3： ${LF}_{j} (x_{r_{11}, j}^{g}) < {LF}_{j} (x_{r_{21}, j}^{g}) and {LF}_{j} (x_{r_{21}, j}^{g}) > {LF}_{j} (x_{r_{31}, j}^{g}) .$ Case3: ${LF}_{j} (x_{r_{11}, j}^{g}) < {LF}_{j} (x_{r_{twenty one}, j}^{g}) and {LF}_{j} (x_{r_{twenty one}, j}^{g}) > {LF}_{j} (x_{r_{31}, j}^{g}) .$

此种情况可以由一个凸函数来拟合，而凸函数意味着在的左边或右边均有可能存在一个更好的变量。This situation can be fitted by a convex function, and a convex function means that in There may be a better variable to the left or right of .

$Case case 44 : : {LF LF}_{j j} (({x x}_{{r r}_{1111},, j j}^{g g})) > > {LF LF}_{j j} (({x x}_{{r r}_{21 twenty one},, j j}^{g g})) and and {LF LF}_{j j} (({x x}_{{r r}_{21 twenty one},, j j}^{g g})) < < {LF LF}_{j j} (({x x}_{{r r}_{3131},, j j}^{g g})) . . . .$

此种情况可以由一个凹函数来拟合，凹函数有一个极小值点，所以它意味着其极小值点可能是一个更好的变量。但是拟合一个凹函数并精确找到其极值点是非常麻烦的，与其相比，取的加权平均是一个更好的选择。This situation can be fitted by a concave function, which has a minimum point, so it means that its minimum point may be a better variable. But it is very troublesome to fit a concave function and find its extremum point precisely. Compared with it, take The weighted average of is a better choice.

总结上面四种情况，新的变异策略可以写为：Summarizing the above four situations, the new mutation strategy can be written as:

其中α为变异尺度因子，一般选取0到1之间的实数。值得注意的是，上面的变异操作在实施时并不需要进行任何曲线拟合。这个变异算子同时拥有随机性和贪婪性。Among them, α is the variation scale factor, generally a real number between 0 and 1 is selected. It is worth noting that the mutation operation above does not require any curve fitting to be implemented. This mutation operator is both random and greedy.

其中CR∈[0,1)是交叉率；rand为[0，1]均匀分布随机数；n_rand是1到D之间的一随机整数，它保证了试验个体中至少有一个基因是来自于变异个体。Among them, CR∈[0,1) is the crossover rate; rand is a uniformly distributed random number in [0,1]; n _rand is a random integer between 1 and D, which ensures that at least one gene in the test individual is from mutant individual.

8）进行选择操作，根据当代个体和试验个体的全局适应度值选取下一代。比较当代个体的全局适应度值和试验个体的全局适应度值，全局适应度值小的个体被选入下一代，具体操作为：8) Carry out the selection operation, and select the next generation according to the global fitness value of the current individual and the test individual. Comparing the global fitness value of the current individual with the global fitness value of the test individual, the individual with a small global fitness value is selected into the next generation, and the specific operation is as follows:

${x x}_{i i}^{g g + + 11} = = \{\begin{matrix} {u u}_{i i}^{g g} & ifF ifF (({u u}_{i i}^{g g})) < < F f (({x x}_{i i}^{g g})),, \\ {x x}_{i i}^{G G} & otherwise otherwise,, \end{matrix} i i = = 1,2 1,2,, . . . . . .,, NP NP . .$

9）g=g+1，判断是否满足终止条件。如果进化代数g小于或等于规定的值G，则返回到步骤6）；否则执行步骤10）。9) g=g+1, judging whether the termination condition is satisfied. If the evolution algebra g is less than or equal to the specified value G, return to step 6); otherwise, perform step 10).

10）选出第G代群体中全局适应度值最小的个体，将该个体解码后即得到最终搜索到的最优地质模型。10) Select the individual with the smallest global fitness value in the G-th generation population, and decode the individual to obtain the optimal geological model that is finally searched.

Claims

1. A collaborative mutation differential evolution algorithm for high-dimensional parameter space waveform inversion, characterized in that, comprising the following steps:

1) Collect the original seismic data, then preprocess the collected seismic data, and obtain the pre-stack common midpoint gather or post-stack seismic data after processing, which is called the observed seismic data, denoted as Seis(r,t ), where r represents the receiving point position of the geophone, and t represents the time axis;

2) Construct a horizontal layered geological model, given the number of layers N and the thickness of the layer in the geological model, the medium model parameters of each layer include P-wave velocity V _p , shear-wave velocity V _s and density ρ;

3) Determine the search space for geological model parameters V _p , V _s , and ρ, and specify the objective function to be optimized;

4) Randomly generate NP geological models in the geological model parameter search space, and perform real number encoding to obtain NP individuals the initial group;

5) Calculate the synthetic seismic data of each stochastic model in the population, and then estimate the global fitness value F _i of each stochastic model corresponding to the individual according to the objective function;

6) Perform a mutation operation to generate a mutant individual for each individual in the population First, high-dimensional individuals are decomposed into a series of subcomponents according to the degree of interdependence between genes, and a local fitness function is assigned to each subcomponent, and the local fitness value of each subcomponent is estimated according to the local fitness function; Three different individuals are selected, and the gradient information is obtained according to the change of the local fitness value of the corresponding subcomponent, and the negative direction of the gradient is used as the variation direction of each subcomponent;

The direction of individual variation is guided by the local fitness value of the subcomponent: randomly select three different individuals, and arrange the genes of the subcomponent from small to large. If the corresponding local fitness value is increasing case1, then the minimum Gene Find a new gene on the left of ; if its corresponding local fitness value is decreasing case2, then the maximum Find new genes on the right side of Find a new gene near ; if its corresponding local fitness value increases first and then decreases case4, then find a new gene on the left of the smallest gene or on the right of the largest gene:

{v v}_{i i,, j j}^{g g} = = {{\begin{matrix} {x x}_{{r r}_{1111},, j j}^{g g} - - α α | | {x x}_{{r r}_{22},, j j}^{g g} - - {x x}_{{r r}_{33},, j j}^{g g} | |,, & c c a a s the s e e 11 \\ {x x}_{{r r}_{3131},, j j}^{g g} + + α α | | {x x}_{{r r}_{22},, j j}^{g g} - - {x x}_{{r r}_{33},, j j}^{g g} | |,, & c c a a s the s e e 22 \\ {x x}_{{r r}_{21 twenty one},, j j}^{g g} + + α α (({x x}_{{r r}_{22},, j j}^{g g} - - {x x}_{{r r}_{33},, j j}^{g g})),, & c c a a s the s e e 33 \\ {p p}_{11} {x x}_{{r r}_{1111},, j j}^{g g} + + {p p}_{22} {x x}_{{r r}_{21 twenty one},, j j}^{g g} + + {p p}_{33} {x x}_{{r r}_{3131},, j j}^{g g} & c c a a s the s e e 44 \end{matrix},, j j = = 11,, 22,, ... ...,, J J

{p p}_{11} = = {p p}_{44} / / {LF LF}_{j j} (({x x}_{{r r}_{1111},, j j}^{g g})),, {p p}_{22} = = {p p}_{44} / / {LF LF}_{j j} (({x x}_{{r r}_{21 twenty one},, j j}^{g g})),, {p p}_{33} = = {p p}_{44} / / {LF LF}_{j j} (({x x}_{{r r}_{3131},, j j}^{g g})),,

{p p}_{44} = = 11 / / (({LF LF}_{j j} {(({x x}_{{r r}_{1111},, j j}^{g g}))}^{- - 11} + + {LF LF}_{j j} {(({x x}_{{r r}_{22} 11,, j j}^{g g}))}^{- - 11} + + {LF LF}_{j j} {(({x x}_{{r r}_{3131},, j j}^{g g}))}^{- - 11}))

Among them, α is the variation scale factor, generally a real number between 0 and 1 is selected, J is the number of subcomponents, and LF _j is the local fitness value corresponding to each gene;

7) Perform a crossover operation, randomly from the mutant individual and corresponding contemporary individuals Extract genes from the group and combine them into a test individual

8) Carry out the selection operation, and select the next generation according to the global fitness value of the current individual and the test individual;

9) Evolutionary algebra g=g+1, judging whether the termination condition is satisfied, if the evolutionary algebra g is less than or equal to the set value G, then return to step 6); otherwise, perform step 10);

10) Select the individual with the smallest global fitness value in the G-th generation population, and decode the individual to obtain the final optimal geological model.

2. According to claim 1, the collaborative variation differential evolution algorithm for high-dimensional parameter space waveform inversion is characterized in that, in step 3), the search space of parameters V _p , V _s , ρ is based on the well logging data The logging record of the corresponding parameters is determined, that is, the low-frequency component of the geological model obtained from the logging is used as the background value, the background value is increased by the set value to be the upper bound of the search space, and the background value is reduced to the set value to be the lower bound of the search space; or The search space of the longitudinal wave velocity V _p and the shear wave velocity V _s can also be determined according to the high-precision velocity analysis; the objective function is a standard for describing the optimal solution to calculate the degree of fit between the wave field and the observed wave field or The size of the error is standard.

3. According to claim 1, the collaborative variation differential evolution algorithm for high-dimensional parameter space waveform inversion is characterized in that, in step 5), the absolute difference between the error between the observed seismic data and the calculated seismic data synthesized by the model is and is the objective function, that is, the objective function is:

{F f}_{i i} = = {Σ Σ}_{r r = = 11}^{N N R R} &Integral; &Integral; | | S S e e i i s the s ((r r,, t t)) - - {Syn Syn}_{i i} ((r r,, t t)) | | d d t t . .

where NR is the number of geophones, and the smaller the fitness value F _i is, the smaller the error between the observed seismic data and the calculated seismic data is, and the corresponding individual The more excellent; on the contrary, the larger the fitness value F _i is, it means that the corresponding individual worse.

4. The co-variation differential evolution algorithm for high-dimensional parameter space waveform inversion according to claim 1, wherein, in step 7),

{u u}_{i i,, n no}^{g g} = = \{\begin{matrix} {v v}_{i i,, n no}^{g g} & i i f f r r a a n no d d \leq \leq C C R R o o r r n no = = {n no}_{r r a a n no d d},, \\ {x x}_{i i,, n no}^{g g} & o o t t h h e e r r w w i i s the s e e \end{matrix},, n no = = 11,, 22,, ... ...,, D D.

Among them, CR∈[0,1) is the crossover rate; rand is a uniformly distributed random number in [0,1]; n _rand is a random integer between 1 and D; D is the dimension of the individual.

5. the co-variation differential evolution algorithm that is used for high-dimensional parameter space waveform inversion according to claim 1, is characterized in that, in step 8), compare the global fitness value of contemporary individual and the global fitness value of test individual, Individuals with small global fitness values are selected into the next generation, and the specific operations are as follows:

{x x}_{i i}^{g g + + 11} = = \{\begin{matrix} {u u}_{i i}^{g g} & i i f f F f (({u u}_{i i}^{g g})) < < F f (({x x}_{i i}^{g g})),, \\ {x x}_{i i}^{G G} & o o t t h h e e r r w w i i s the s e e,, \end{matrix},, i i = = 11,, 22,, ... ...,, N N P P .. ..