CN112649860B - Layer velocity inversion method and system based on continuous ant colony algorithm - Google Patents
Layer velocity inversion method and system based on continuous ant colony algorithm Download PDFInfo
- Publication number
- CN112649860B CN112649860B CN201910966262.6A CN201910966262A CN112649860B CN 112649860 B CN112649860 B CN 112649860B CN 201910966262 A CN201910966262 A CN 201910966262A CN 112649860 B CN112649860 B CN 112649860B
- Authority
- CN
- China
- Prior art keywords
- layer
- speed
- target layer
- grid
- umin
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 39
- 238000012545 processing Methods 0.000 claims abstract description 5
- 239000003016 pheromone Substances 0.000 claims description 37
- 241000257303 Hymenoptera Species 0.000 claims description 14
- 230000015572 biosynthetic process Effects 0.000 claims description 11
- 230000009191 jumping Effects 0.000 claims 1
- 230000000694 effects Effects 0.000 abstract description 2
- 238000005553 drilling Methods 0.000 description 11
- 238000005755 formation reaction Methods 0.000 description 9
- 238000010586 diagram Methods 0.000 description 3
- 238000004364 calculation method Methods 0.000 description 2
- 238000005070 sampling Methods 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000007796 conventional method Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/30—Analysis
- G01V1/303—Analysis for determining velocity profiles or travel times
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/60—Analysis
- G01V2210/62—Physical property of subsurface
- G01V2210/622—Velocity, density or impedance
- G01V2210/6222—Velocity; travel time
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Acoustics & Sound (AREA)
- Environmental & Geological Engineering (AREA)
- Geology (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
- Earth Drilling (AREA)
Abstract
The invention provides a continuous ant colony algorithm-based interval velocity inversion method and system, and belongs to the field of seismic data processing. The method comprises the following steps: setting a maximum layer speed Umax and a minimum layer speed Umin of a work area; setting the number n of target layers to be inverted; mesh subdivision is carried out on the speed interval [ Umin, umax ] to obtain a speed mesh; and obtaining the optimal solution of the layer velocity of each target layer through multiple iterations by utilizing the velocity grid and the continuous ant colony algorithm. By using the method, the layer velocity can be rapidly solved, the global optimal solution is found by crossing the local optimal solution, the accuracy of the inversion of the layer velocity is high, and the effect is good.
Description
Technical Field
The invention belongs to the field of seismic data processing, and particularly relates to a continuous ant colony algorithm-based interval velocity inversion method and system.
Background
Vertical seismic data (VSP) contains abundant seismic wave information, and the pre-drilling prediction by utilizing the VSP data is an important geophysical means for obtaining information of a target layer in front of a drill bit. The prediction content before drilling comprises the longitudinal wave speed, the depth, the formation pressure and the like below the well bottom, and can provide assistance for a drilling process and reduce the drilling risk. The conventional method for predicting the velocity of the VSP front layer before drilling is to invert the reflection coefficient before drilling by adopting a proper inversion method on the basis of zero-offset VSP data, and then calculate the layer-by-layer velocity of a stratum by the reflection coefficient so as to achieve the purpose of predicting the velocity of the front layer before drilling.
The domestic scholars put forward different researches on the VSP data inversion layer speed. Chen Xinping (1992) an algorithm for inverting the layer velocity according to VSP first arrival time is proposed based on Redshaw modified gaussian-newton algorithm, and the error of the fitting result is analyzed. Zhang Houzhu et al (1995) discuss the stability of genetic algorithm and its sensitivity to noise, etc., and it is considered that the genetic algorithm has strong anti-interference ability, and it has potential ability to invert the horizon speed of actual data by using the genetic algorithm. Li Wenjie, etc. (2004) utilize VSP data to extract the reflection coefficient of the underground formation, and then utilize the reflection coefficient to invert the layer-to-layer velocities of the various formations, by utilizing the characteristics of the VSP data that the interference degree is small, the signal-to-noise ratio is high, the high frequency components are rich, and the up-going wave and the down-going wave can be accurately separated out. Huang Handong et al (2007) introduce a method and principle of an ant colony algorithm for layer velocity inversion, start with the application of the ant colony algorithm to solve an objective function optimization problem, and utilize a model to test the inversion result of the method. Zhou (2012) introduced a method for inverting the stratum velocity layer by layer according to VSP data, which is mainly characterized by dividing the propagation path of the direct wave ray into straight line propagation and broken line propagation, and inverting the stratum velocity of each stratum from top to bottom by using the VSP data. Yellow light south et al (2016) can reconstruct an underground stratum velocity model by combining a three-dimensional VSP seismic ray tracing algorithm and a digital inversion method.
The prediction of the velocity, the depth and the formation pressure of the longitudinal wave and the transverse wave before drilling of the VSP data can provide help for a drilling process and reduce the drilling risk. However, the above layer velocity prediction method has certain limitations, such as relying on an initial model, and easily falling into a locally optimal solution.
Disclosure of Invention
The invention aims to solve the problems in the prior art, provides a layer velocity inversion method and system based on a continuous ant colony algorithm, overcomes the limitation that the conventional ant colony algorithm is easy to fall into a local optimal solution, and improves the stability of layer velocity inversion.
The invention is realized by the following technical scheme:
a layer velocity inversion method based on a continuous ant colony algorithm comprises the following steps:
setting a maximum layer speed Umax and a minimum layer speed Umin of a work area;
setting the number n of target layers to be inverted;
mesh subdivision is carried out on the speed interval [ Umin, umax ] to obtain a speed mesh;
and obtaining the optimal solution of the layer velocity of each target layer through multiple iterations by using the velocity grid and the continuous ant colony algorithm.
The operation of meshing the speed interval [ Umin, umax ] to obtain the speed mesh comprises the following steps:
setting the maximum layer speed of each target layer to be Umax, and setting the minimum layer speed of each target layer to be Umin;
setting the number of equally dividing the speed interval as a speed division grid number M;
calculating a subdivision grid variable h by using the following formula:
h＝(Umax-Umin)/M；
dividing a speed interval by using h to obtain a speed grid, wherein each row in the speed grid corresponds to one speed, and the speeds are Umin, umin + h, umin +2h to Umax-h and Umax from bottom to top in sequence;
each column in the velocity grid corresponds to a target layer to be inverted, and the target layer is 1,2 to n from left to right in sequence;
and intersecting rows and columns in the speed grid to obtain a plurality of grid nodes, wherein each grid node corresponds to a speed value on a destination layer.
The operation of obtaining the optimal solution of the layer velocity of each target layer through multiple iterations by using the velocity grid and the continuous ant colony algorithm comprises the following steps:
calculating to obtain an optimal objective function value by using the maximum layer speed Umax of each target layer and an objective function value formula;
step two: setting the number Nvar of ants and the total iteration times D; setting an initial value of iteration number i: i =1;
step three: assigning the initial value of local pheromone of each grid node in the speed grid to be 1/M, assigning the initial value of global pheromone to be 0, and assigning the initial value of grid node probability of each grid node to be 0;
step four: randomly generating the position of a grid node of each target layer of an ant, and obtaining the speed value of the ant on each target layer according to the position of the grid node;
step five: calculating the objective function value of the ant by using the speed value of the ant on each target layer and an objective function value formula;
step six: judging whether the objective function value of the ant is smaller than the optimal objective function value, if so, updating the local pheromone, the global pheromone and the probability of each grid node of the grid node on each target layer where the ant is positioned by using the objective function value of the ant, and if not, not updating;
step seven: sequentially carrying out the treatment of the fourth step, the fifth step and the sixth step on each ant until Nvar ants are completely treated;
step eight: obtaining the optimal solution of the layer speed of each target layer, and updating the optimal objective function value;
step nine: updating the grid nodes on each destination layer, wherein the grid nodes on all the destination layers form a new speed grid;
step ten: judging whether i = D is true, if so, turning to the step eleven, otherwise, turning to the step III, and returning to the step III, wherein i = i + 1;
and step eleven, outputting the layer velocity optimal solution of each target layer obtained by the last iteration.
The operation of updating the local pheromone, the global pheromone and the probability of each grid node of the grid node on each destination layer where the ant is positioned by using the objective function value of the ant in the sixth step comprises the following steps:
and updating local pheromones of the grid nodes on each destination layer where the ant is positioned by using the following formula:
Tnew＝(1-Rou)*Told+Q/y
wherein, told is the current local information content of the grid node, tnew is the updated local information content of the grid node, rou is the information content residual factor, and Q is the information content intensity;
the global pheromone Ts is updated using the following equation:
Ts＝sum(Tnew)；
and updating the grid node probability P of the grid node on each destination layer where the ant is positioned by using the following formula:
P＝Tnew/Ts。
the operation of the step eight comprises the following steps:
and respectively carrying out the following treatment on each target layer: comparing the probability of each grid node on the target layer to obtain the position max _ ind of the grid node corresponding to the maximum grid node probability, and obtaining the layer speed optimal solution local of the target layer by using the position of the grid node:
local＝Umin+(Umax-Umin)*(max_ind-1)/M；
after all the target layers are processed, a new objective function value is calculated by using the layer speed optimal solution and the objective function value formula of each target layer, and the new objective function value is assigned to the optimal objective function value.
The objective function value formula is as follows:
y＝sum((R*W-S) 2 )
wherein y is an objective function value, R is a formation reflection coefficient, R is a product of a layer velocity V and a formation density, W is a seismic source wavelet, and S is an actual seismic trace record.
The operation of the ninth step comprises the following steps:
and respectively carrying out the following treatment on each target layer: and updating the maximum layer speed Umax and the minimum layer speed Umin of the target layer by using the max _ ind and the floating index delt of the target layer, and then calculating a subdivision grid variable h of the target layer: h = (Umax-Umin)/M, and finally, carrying out speed mesh division on the speed of the target layer by using h to obtain each mesh node on the target layer.
Further, in order to avoid falling into the locally optimal solution, the method further includes, in step two: setting a cycle upper limit value count;
further comprising between step nine and step ten:
avoiding falling into the local optimal solution step: and respectively carrying out the following treatment on each target layer: judging whether the count minimum layer speeds Umin obtained by the target layer after the continuous count iterations are the same or not, if so, judging that the Umin obtained by the target layer after the continuous count iterations is the same, and if not, judging that the Umin obtained by the target layer after the continuous count iterations is different;
after all the target layers are processed, judging whether at least one target layer with different Umin obtained by continuous count iterations exists, if so, judging that the target layer does not fall into the local optimal solution, then, turning to the step ten, if not, judging that the target layer falls into the local optimal solution, and then, respectively processing each target layer as follows: randomly placing max _ ind of the target layer, then updating the maximum layer speed Umax and the minimum layer speed Umin of the target layer by using the max _ ind, and then calculating a mesh variable h of the target layer: h = (Umax-Umin)/M, the speed mesh division is carried out on the speed of the target layer by using h to obtain each mesh node of the target layer, all mesh nodes on the target layer form a new speed mesh, and then the step ten is carried out.
The operation of updating the maximum layer speed Umax and the minimum layer speed Umin of the destination layer comprises:
the maximum layer velocity for the destination layer is updated using:
Umax_new＝Umin_old+(max_ind+delt)*(Umax_old-Umin_old)/M；
the minimum layer velocity for the destination layer is updated using the following equation:
Umin_new＝Umin_old+(max_ind-delt)*(Umax_old-Umin_old)/M；
wherein Umax _ new and Umin _ new are the maximum layer speed and the minimum layer speed of each target layer after updating, and Umax _ old and Umin _ old are the maximum layer speed and the minimum layer speed of each target layer before updating.
The invention also provides a layer velocity inversion system based on the continuous ant colony algorithm, which comprises:
the layer speed setting unit is used for setting the maximum layer speed Umax and the minimum layer speed Umin of the work area;
the target layer number setting unit is used for setting the number n of target layers to be inverted;
the mesh generation unit is respectively connected with the layer speed setting unit and the target layer number setting unit and is used for carrying out mesh generation on a speed interval [ Umin, umax ] to obtain a speed mesh;
and the layer velocity inversion unit is connected with the mesh generation unit and is used for obtaining the layer velocity optimal solution of each target layer through multiple iterations by utilizing the velocity mesh and the continuous ant colony algorithm.
Compared with the prior art, the invention has the beneficial effects that: by using the method, the layer velocity can be rapidly solved, the global optimal solution is found by crossing the local optimal solution, the accuracy of the inversion of the layer velocity is high, and the effect is good.
Drawings
FIG. 1 is a block diagram of the steps of the method of the present invention
FIG. 2 is a schematic diagram of a velocity mesh generation;
FIG. 3 shows the results of a first inversion in an embodiment of the method of the invention;
FIG. 4 shows the results of a second inversion in an embodiment of the method of the invention;
FIG. 5 shows the result of a third inversion in an embodiment of the method of the invention;
FIG. 6 is a block diagram of the components of the system of the present invention.
Detailed Description
The invention is described in further detail below with reference to the accompanying drawings:
the method improves the global information element updating and the local information element updating, sets the cycle upper limit, and adopts the newly generated variable group to jump out of the cycle when the optimal solution is continuously and repeatedly equal, thereby avoiding falling into the local optimal solution. The method can be used for the layer velocity inversion in VSP data to predict the layer velocity before drilling.
Specifically, the method roughly determines the maximum speed and the minimum speed in a work area range, and carries out mesh generation on the speed; the artificial ants move among the speed space grids, and information content is left, so that the moving direction of the next batch of ants is influenced; and when iteration is circulated to the given total iteration times, finding out the grid point with the maximum information quantity according to the information quantity, thereby determining the speed value of each speed sampling point. In order to avoid falling into the local optimal solution, the rules of global pheromone updating and local pheromone updating are changed, the upper limit of the circulation is set, and when the optimal solution is continuously and repeatedly identical, the newly generated variable group is adopted to jump out of the circulation
The basic idea of the method of the invention is as follows: and setting a speed range to be inverted according to past experience information of a past work area. And then, carrying out spatial meshing subdivision on the speed, wherein each grid point corresponds to a state, ants are randomly placed on the grid points of each stratum, pheromone and probability at the node are calculated, the grid position with the maximum node probability is obtained, the maximum and minimum stratum speed ranges of each stratum are updated according to the grid position, the grid is reduced, and next iteration is carried out. In order to avoid falling into a local optimal solution, a loop upper limit value is set, and when the optimal solutions are the same within the number of times of the loop upper limit value, a new generation variable group is adopted and the loop is jumped out.
Specifically, as shown in fig. 1, the method of the present invention comprises:
setting a maximum layer speed Umax and a minimum layer speed Umin of a work area and setting the number n of target layers to be inverted according to past experience information of the past work area; the interval velocity of each destination layer to be inverted represents a velocity variable. All velocity variables that need inversion fall within the [ Umin, umax ] velocity interval.
The method comprises the following steps: setting the maximum layer speed of each target layer to be Umax, and setting the minimum layer speed of each target layer to be Umin; dividing the speed interval equally, setting the number of the divided intervals as the speed division grid number M, and setting the variable size of the division grid as h, h = (Umax-Umin)/M; and (2) performing speed grid subdivision on the speed to obtain an initial speed grid, as shown in fig. 2, each row in the initial speed grid corresponds to one speed, sequentially from bottom to top, umin + h, umin +2h to Umax-h and Umax, each column in the initial speed grid corresponds to a target layer to be inverted, sequentially from left to right, 1 and 2 … … n target layers, the rows and the columns intersect to obtain a plurality of grid nodes, each grid node represents a possible layer speed on the target layer corresponding to the column where the grid node is located, and the grid nodes are represented by V1 and V2 … … Vn in fig. 2.
And calculating to obtain an optimal objective function value by using the maximum layer speed Umax of each stratum and an objective function value formula, wherein the initial optimal objective function value can be updated along with the updating of the range of the maximum speed values of each stratum in the later iteration. The objective function value is calculated by the standard deviation of the synthetic seismic trace R x W synthesized by the maximum velocity of each stratum and the actual seismic trace S, and the objective function value formula is as follows:
y＝sum((R*W-S) 2 )
wherein y is the objective function value, R is the formation reflection coefficient (the product of the layer velocity V and the formation density), W is the seismic source wavelet, and S is the actual seismic trace record. The formation density, source wavelet, and actual seismic trace recordings are all known. In the first step of calculation, because the maximum interval velocity of each target interval is Umax, the formation reflection coefficient R is obtained by multiplying the maximum interval velocity and the formation density.
Step two: the number Nvar of ants is set and can be generally selected to be 1/4 of the mesh number M of the speed subdivision. Setting the total iteration times D; setting a loop upper limit value count, wherein the count is far smaller than the total iteration number D and can be generally 5; setting an initial value of iteration times i: i =1;
step three: and respectively assigning initial values to the local pheromone, the global pheromone and the grid probability. The local pheromone refers to the information quantity on each grid node, the global pheromone refers to the sum of the information quantities on all grid nodes of each stratum, and the grid probability refers to the probability value on each grid point. In this step, the initial values of the local pheromones of the respective mesh nodes are all given 1/M, the initial value of the global pheromone is given 0, and the initial value of the mesh probability is given 0, and the values of these three parameters are updated in the subsequent steps.
Step four: randomly generating the position of a grid node of an ant on each destination layer by using a vector (P) 1 ,P 2 ,…P n ) Represents the element P j (j =1,2.. N) represents the location of the ant at the mesh node corresponding to the jth stratum. And obtaining the speed value V of the ant on each destination layer from the positions of the grid nodes (the speed value corresponding to the ordinate of the grid node on each row where the ant is positioned is the speed value V of the ant on the destination layer corresponding to the row).
Step five: and calculating the objective function value of the ant by using the speed value V of the ant at each destination layer and the objective function value formula.
Step six: and judging whether the objective function value of the ant is smaller than the optimal objective function value or not. If the objective function value is smaller than the optimal objective function value, updating local pheromones of the grid nodes on each destination layer where the ants are located by using the objective function values: tnew = (1-Rou) × Told + Q/y, where Told is the current local traffic, tnew is the updated local traffic, rou is the traffic residual factor, rou =0.3, Q is the traffic intensity, Q =1.0 × 1O 8 . Calculating a global pheromone Ts, wherein the global pheromone is the sum of local pheromones of all nodes: ts = sum (Tnew). Calculating the probability P of each grid node, P = Tnew/Ts。
And if the objective function value is larger than the optimal target value, not updating the local pheromone, the global pheromone and the grid node probability of the ant.
Step seven: and (4) sequentially carrying out the treatment of the fourth step, the treatment of the fifth step and the treatment of the sixth step on each ant to obtain the local pheromone, the global pheromone and the grid node probability after all the ants (namely Nvar ants) are placed.
Step eight: obtaining a speed optimal solution of each target layer: and respectively carrying out the following treatment on each target layer: comparing the probability of each grid node on the target layer to obtain the position max _ ind of the grid node corresponding to the maximum grid node probability of the target layer, and obtaining the layer speed optimal solution corresponding to the maximum probability grid node of the target layer by using the position of the grid node: local = Umin + (Umax-Umin) × (max _ ind-1)/M. And calculating a new objective function value by using the speed optimal solution of each target layer and the objective function value formula, and assigning the new objective function value to the optimal objective function value, namely updating the optimal objective function value into the new objective function value.
Step nine: and respectively carrying out the following treatment on each target layer: and updating the maximum layer speed Umax and the minimum layer speed Umin of the target layer by using the max _ ind and the floating index delt of the target layer, and calculating a subdivision grid variable h of the target layer: h = (Umax-Umin)/M, the speed mesh division is carried out on the speed of the target layer by using h to obtain each mesh node of the target layer, and all the mesh nodes form a new speed mesh; the delt index is a delt grid node which is floated up and down by taking the optimal solution as a center, the delt is a natural number, and the value can be adjusted according to actual needs. After each iteration, the maximum layer speed and the minimum layer speed of each destination layer may not be the same any more, the grids of each destination layer are not the same any more, and each destination layer obtains each grid node of the destination layer according to h of the destination layer, that is, each grid node is on each vertical line in fig. 2, but a horizontal line is not intersected with all vertical lines to obtain grid nodes on the same horizontal line.
Specifically, the maximum layer velocity Umax of the destination layer is updated using the following equation:
Umax_new＝Umin_old+(max_ind+delt)*(Umax_old-Umin_old)/M，
the minimum layer velocity for the destination layer is updated using the following equation:
Umin_new＝Umin_old+(max_ind-delt)*(Umax_old-Umin_old)/M，
wherein Umax _ new and Umin _ new are the maximum and minimum layer speeds of each stratum after updating, and Umax _ old and Umin _ old are the maximum and minimum layer speeds of each stratum before updating.
Step ten: and judging whether i = D is established, if so, shifting to the step eleven, otherwise, returning to the step three, wherein i = i + 1.
And step eleven, outputting the layer velocity optimal solution of each target layer obtained by the last iteration. With the increase of the iteration times, the speed range is continuously reduced, and finally the layer speed values of all the target layers at the inversion position in the given iteration times, namely the optimal solution of the layer speeds of the n target layers, is obtained.
Further, in order to avoid falling into the local optimal solution, the method further comprises, between step nine and step ten:
and sequentially performing the following treatment on each target layer respectively: judging whether count minimum layer speeds Umin obtained by the target layer after continuous count iterations are the same or not, if yes, judging that the Umin obtained by the target layer through the continuous count iterations is the same, and if not, judging that the Umin obtained by the target layer through the continuous count iterations is different;
judging whether at least one target layer with different Umin obtained by continuous count iteration exists, if so, judging that the local optimal solution does not fall into, then turning to the step ten, if not, judging that the local optimal solution falls into (namely judging that the local optimal solution falls into only under the condition that the Umin obtained by continuous count iteration of all the target layers is the same), and then respectively carrying out the following processing on each target layer: randomly placing the max _ ind of the target layer, namely reassigning the max _ ind of the target layer to a random position, then updating the maximum layer speed and the minimum layer speed of each target layer by using the randomly placed max _ ind of each target layer, and then calculating a subdivision grid variable h of the target layer: h = (Umax-Umin)/M, the speed grid of the target layer is obtained by carrying out speed grid division on the speed of the target layer by h (similar to the equal division in the step one, only each target layer is respectively equally divided, and the number from uman to Umax of the target layer is equally divided into M to obtain each grid node); then step ten is carried out.
Specifically, the maximum layer speed Umax is updated using the following equation:
Umax_new＝Umin_old+(max_ind+delt)*(Umax_old-Umin_old)/M，
the minimum layer velocity is updated using:
Umin_new＝Umin_old+(max_ind-delt)*(Umax_old-Umin_old)/M，
wherein Umax _ new and Umin _ new are the maximum and minimum layer speeds of each stratum after updating, and Umax _ old and Umin _ old are the maximum layer speed and minimum layer speed of each stratum before updating.
As shown in fig. 6, the present invention further provides a continuous ant colony algorithm based layer velocity inversion system, which includes:
a layer speed setting unit 10, configured to set a maximum layer speed Umax and a minimum layer speed Umin of a work area;
a target layer number setting unit 20, configured to set a number n of target layers to be inverted;
the mesh generation unit 30 is respectively connected with the layer speed setting unit 10 and the target layer number setting unit 20 and is used for carrying out mesh generation on a speed interval [ Umin, umax ] to obtain a speed mesh;
and the layer velocity inversion unit 40 is connected with the mesh generation unit 30 and is used for obtaining the layer velocity optimal solution of each target layer through multiple iterations by utilizing the velocity mesh and the continuous ant colony algorithm.
The examples of the invention are as follows:
to test the correctness and feasibility of this inversion method, a theoretical model as shown in table 1 was designed. The model was divided into 7 layers in total.
| |
1 | 2 | 3 | 4 | 5 | 6 | 7 |
| Speed/(m/s) | 700 | 2000 | 1800 | 1700 | 2100 | 1900 | 1300 |
TABLE 1
A zero-phase Rake wavelet with a main frequency of 30Hz and a sampling rate of 1ms is adopted. Basic parameters used in the inversion process: the number of grid points is 200, the number of ants is 50, and the number of iterations is 200. Fig. 3 to 5 are graphs comparing inversion results obtained after three times of inversion with actual layer velocities. As can be seen from the comparison graph, the method has good stability, and the calculation error is also within the allowable range.
The above-described embodiment is only one embodiment of the present invention, and it will be apparent to those skilled in the art that various modifications and variations can be easily made based on the application and principle of the present invention disclosed in the present application, and the present invention is not limited to the method described in the above-described embodiment of the present invention, so that the above-described embodiment is only preferred, and not restrictive.
Claims (8)
1. A layer velocity inversion method based on a continuous ant colony algorithm is characterized by comprising the following steps: the method comprises the following steps:
setting a maximum layer speed Umax and a minimum layer speed Umin of a work area;
setting the number n of target layers to be inverted;
mesh subdivision is carried out on the speed interval [ Umin, umax ] to obtain a speed mesh;
obtaining a layer velocity optimal solution of each target layer through multiple iterations by utilizing the velocity grid and the continuous ant colony algorithm;
setting a cycle upper limit value count, and jumping out of the cycle by adopting a newly generated variable group when the continuous counts of the optimal solution are the same for the times;
the operation of obtaining the optimal solution of the layer velocity of each target layer through multiple iterations by using the velocity grid and the continuous ant colony algorithm comprises the following steps:
step one, calculating to obtain an optimal objective function value by utilizing the maximum layer speed Umax and an objective function value formula of each objective layer;
step two: setting the number Nvar of ants; setting the total iteration times D, and setting the initial value of the iteration times i: i =1;
step three: assigning the initial value of local pheromone of each grid node in the speed grid to be 1/M, assigning the initial value of global pheromone to be 0, and assigning the initial value of grid node probability of each grid node to be 0;
step four: randomly generating the position of a grid node of each target layer of an ant, and obtaining the speed value of the ant on each target layer according to the position of the grid node;
step five: calculating the objective function value of the ant by using the speed value of the ant on each target layer and an objective function value formula;
step six: judging whether the objective function value of the ant is smaller than the optimal objective function value, if so, updating the local pheromone, the global pheromone and the probability of each grid node of the grid node on each target layer where the ant is positioned by using the objective function value of the ant, and if not, not updating;
step seven: sequentially carrying out the treatment of the fourth step, the fifth step and the sixth step on each ant until Nvar ants are completely treated;
step eight: obtaining the optimal solution of the layer speed of each target layer, and updating the optimal objective function value;
step nine: updating the grid nodes on each destination layer, wherein the grid nodes on all the destination layers form a new speed grid;
step ten: judging whether i = D is true, if so, turning to the step eleven, otherwise, turning to the step III, and returning to the step III, wherein i = i + 1;
step eleven, outputting the layer velocity optimal solution of each target layer obtained by the last iteration;
the operation of the ninth step comprises the following steps:
and respectively carrying out the following treatment on each target layer: and updating the maximum layer speed Umax and the minimum layer speed Umin of the target layer by using the max _ ind and the floating index delt of the target layer, and then calculating a subdivision grid variable h of the target layer: h = (Umax-Umin)/M, and finally, carrying out speed mesh subdivision on the speed of the target layer by using h to obtain each mesh node on the target layer; wherein, max _ ind is the position of the grid node corresponding to the maximum grid node probability; delt is delt grid nodes which float up and down by taking the optimal solution as a center, and delt is a natural number; m is the speed subdivision grid number.
2. The continuous ant colony algorithm-based interval velocity inversion method according to claim 1, wherein: the operation of mesh generation on the speed interval [ Umin, umax ] to obtain the speed mesh comprises the following steps:
setting the maximum layer speed of each target layer to be Umax, and setting the minimum layer speed of each target layer to be Umin;
setting the number of equally dividing the speed interval as a speed division grid number M;
calculating a subdivision grid variable h by using the following formula:
h＝(Umax-Umin)/M；
dividing a speed interval by using h to obtain a speed grid, wherein each row in the speed grid corresponds to one speed, and the speeds are Umin, umin + h, umin +2h to Umax-h and Umax from bottom to top in sequence;
each column in the velocity grid corresponds to a target layer to be inverted, and the target layer is 1,2 to n from left to right in sequence;
and intersecting rows and columns in the speed grid to obtain a plurality of grid nodes, wherein each grid node corresponds to a speed value on a destination layer.
3. The continuous ant colony algorithm-based interval velocity inversion method according to claim 2, wherein: the operation of updating the local pheromone, the global pheromone and the probability of each grid node of the grid node on each destination layer where the ant is positioned by using the objective function value of the ant in the sixth step comprises the following steps:
and updating local pheromones of the grid nodes on each destination layer where the ant is positioned by using the following formula:
Tnew＝(1-Rou)*Told+Q/y
wherein, told is the current local information content of the grid node, tnew is the updated local information content of the grid node, rou is the information content residual factor, and Q is the information content intensity; y is the objective function value;
the global pheromone Ts is updated using the following equation:
Ts＝sum(Tnew)；
and updating the grid node probability P of the grid node on each destination layer where the ant is positioned by using the following formula:
P＝Tnew/Ts。
4. the continuous ant colony algorithm-based interval velocity inversion method according to claim 3, characterized in that: the operation of the step eight comprises the following steps:
and respectively carrying out the following treatment on each target layer: comparing the probability of each grid node on the target layer to obtain the position max _ ind of the grid node corresponding to the maximum grid node probability, and obtaining the layer speed optimal solution local of the target layer by using the position of the grid node:
local＝Umin+(Umax-Umin)*(max_ind-1)/M；
after all the target layers are processed, a new objective function value is calculated by using the layer speed optimal solution and the objective function value formula of each target layer, and the new objective function value is assigned to the optimal objective function value.
5. The continuous ant colony algorithm-based interval velocity inversion method according to claim 4, wherein: the objective function value formula is as follows:
y＝sum((R*W-S) 2 )
wherein y is an objective function value, R is a formation reflection coefficient, R is a product of a layer velocity V and a formation density, W is a seismic source wavelet, and S is an actual seismic trace record.
6. The continuous ant colony algorithm-based interval velocity inversion method according to claim 5, wherein: the second step further comprises: setting a cycle upper limit value count;
further comprising between step nine and step ten:
avoiding trapping in a local optimal solution step: and respectively carrying out the following treatment on each target layer: judging whether the count minimum layer speeds Umin obtained by the target layer after the continuous count iterations are the same or not, if so, judging that the Umin obtained by the target layer after the continuous count iterations is the same, and if not, judging that the Umin obtained by the target layer after the continuous count iterations is different;
after all the target layers are processed, judging whether at least one target layer with different Umin obtained by continuous count iterations exists, if so, judging that the target layer does not fall into the local optimal solution, then, turning to the step ten, if not, judging that the target layer falls into the local optimal solution, and then, respectively processing each target layer as follows: randomly placing max _ ind of the target layer, then updating the maximum layer speed Umax and the minimum layer speed Umin of the target layer by using the max _ ind, and then calculating a mesh variable h of the target layer: h = (Umax-Umin)/M, the speed mesh division is carried out on the speed of the target layer by using h to obtain each mesh node of the target layer, all mesh nodes on the target layer form a new speed mesh, and then the step ten is carried out.
7. The continuous ant colony algorithm-based interval velocity inversion method according to claim 6, wherein: the operation of updating the maximum layer speed Umax and the minimum layer speed Umin of the destination layer comprises:
the maximum layer velocity for the destination layer is updated using:
Umax_new＝Umin_old+(max_ind+delt)*(Umax_old-Umin_old)/M；
the minimum layer velocity for the destination layer is updated using the following equation:
Umin_new＝Umin_old+(max_ind-delt)*(Umax_old-Umin_old)/M；
the Umax _ new and the Umin _ new are respectively the maximum layer speed and the minimum layer speed of each target layer after updating, and the Umax _ old and the Umin _ old are respectively the maximum layer speed and the minimum layer speed of each target layer before updating.
8. A layer velocity inversion system based on continuous ant colony algorithm is characterized in that: the system comprises:
the layer speed setting unit is used for setting the maximum layer speed Umax and the minimum layer speed Umin of the work area;
the target layer number setting unit is used for setting the number n of target layers to be inverted;
the mesh generation unit is respectively connected with the layer speed setting unit and the target layer number setting unit and is used for carrying out mesh generation on a speed interval [ Umin, umax ] to obtain a speed mesh;
the layer velocity inversion unit is connected with the mesh generation unit and is used for obtaining a layer velocity optimal solution of each target layer through multiple iterations by utilizing the velocity mesh and a continuous ant colony algorithm;
the layer velocity inversion unit sets a cycle upper limit value count, and when the continuous counts of the optimal solution are the same for the times, a cycle is skipped by using a newly generated variable group;
the operation of obtaining the optimal solution of the layer velocity of each target layer through multiple iterations by using the velocity grid and the continuous ant colony algorithm comprises the following steps:
calculating to obtain an optimal objective function value by using the maximum layer speed Umax of each target layer and an objective function value formula;
step two: setting the number Nvar of ants; setting the total iteration times D, and setting the initial value of the iteration times i: i =1;
step three: assigning the initial value of local pheromone of each grid node in the speed grid to be 1/M, assigning the initial value of global pheromone to be 0, and assigning the initial value of grid node probability of each grid node to be 0;
step four: randomly generating the position of a grid node of each target layer of an ant, and obtaining the speed value of the ant on each target layer according to the position of the grid node;
step five: calculating the objective function value of the ant by using the speed value of the ant on each target layer and an objective function value formula;
step six: judging whether the objective function value of the ant is smaller than the optimal objective function value, if so, updating the local pheromone, the global pheromone and the probability of each grid node of the grid node on each target layer where the ant is positioned by using the objective function value of the ant, and if not, not updating;
step seven: sequentially carrying out the treatment of the fourth step, the fifth step and the sixth step on each ant until Nvar ants are completely treated;
step eight: obtaining the optimal solution of the layer speed of each target layer, and updating the optimal objective function value;
step nine: updating the grid nodes on each destination layer, wherein the grid nodes on all the destination layers form a new speed grid;
step ten: judging whether i = D is true, if so, turning to the step eleven, otherwise, i = i +1, and then returning to the step three;
step eleven, outputting the layer velocity optimal solution of each target layer obtained by the last iteration;
the operation of the step nine comprises the following steps:
and respectively carrying out the following treatment on each target layer: and updating the maximum layer speed Umax and the minimum layer speed Umin of the target layer by using the max _ ind and the floating index delt of the target layer, and then calculating a subdivision grid variable h of the target layer: h = (Umax-Umin)/M, and finally, carrying out speed mesh subdivision on the speed of the target layer by using h to obtain each mesh node on the target layer; wherein max _ ind is the position of the grid node corresponding to the maximum grid node probability; delt is delt grid nodes which float up and down by taking the optimal solution as a center, and is a natural number; m is the speed subdivision grid number.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910966262.6A CN112649860B (en) | 2019-10-12 | 2019-10-12 | Layer velocity inversion method and system based on continuous ant colony algorithm |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910966262.6A CN112649860B (en) | 2019-10-12 | 2019-10-12 | Layer velocity inversion method and system based on continuous ant colony algorithm |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN112649860A CN112649860A (en) | 2021-04-13 |
| CN112649860B true CN112649860B (en) | 2022-10-14 |
Family
ID=75343594
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201910966262.6A Active CN112649860B (en) | 2019-10-12 | 2019-10-12 | Layer velocity inversion method and system based on continuous ant colony algorithm |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN112649860B (en) |
Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103440527A (en) * | 2013-07-29 | 2013-12-11 | 辽宁大学 | Method for improving ant colony algorithm optimization support vector machine parameters |
Family Cites Families (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US8112374B2 (en) * | 2005-11-23 | 2012-02-07 | Henry Van Dyke Parunak | Hierarchical ant clustering and foraging |
| CN104199106B (en) * | 2014-09-11 | 2017-01-18 | 中国石油集团川庆钻探工程有限公司地球物理勘探公司 | Seismic data residual static correction method based on ant colony algorithm |
| WO2019027449A1 (en) * | 2017-08-01 | 2019-02-07 | Halliburton Energy Services, Inc. | Prediction ahead of bit using vertical seismic profile data and global inversion |
| CN109214498A (en) * | 2018-07-10 | 2019-01-15 | 昆明理工大学 | Ant group algorithm optimization method based on search concentration degree and dynamic pheromone updating |
| CN109001801B (en) * | 2018-07-30 | 2021-03-19 | 中国石油化工股份有限公司 | Fault variable-scale identification method based on multiple iteration ant colony algorithm |
-
2019
- 2019-10-12 CN CN201910966262.6A patent/CN112649860B/en active Active
Patent Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103440527A (en) * | 2013-07-29 | 2013-12-11 | 辽宁大学 | Method for improving ant colony algorithm optimization support vector machine parameters |
Also Published As
| Publication number | Publication date |
|---|---|
| CN112649860A (en) | 2021-04-13 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| US11016214B2 (en) | Dolomite reservoir prediction method and system based on well and seismic combination, and storage medium | |
| CN106779148B (en) | A kind of method for forecasting wind speed of high speed railway line of multi-model multiple features fusion | |
| CN105093319B (en) | Ground micro-seismic static correcting method based on 3D seismic data | |
| CN106772577B (en) | Source inversion method based on microseism data and SPSA optimization algorithm | |
| CN109738940B (en) | Acoustic emission/microseismic event positioning method under condition of existing empty zone | |
| CN108508482A (en) | A kind of subterranean fracture seismic scattering response characteristic analogy method | |
| CN104360396B (en) | A kind of three kinds of preliminary wave Zoumaling tunnel methods of TTI medium between offshore well | |
| CN103984996B (en) | Optimization prediction method for taboo searching algorithm and genetic algorithm of algal bloom mechanism time varying model | |
| CN116774292B (en) | Seismic wave travel time determining method, system, electronic equipment and storage medium | |
| CN104422969A (en) | Method for reducing non-uniqueness of electromagnetic sounding inversion result | |
| CN105607119B (en) | Near-surface model construction method and static correction value acquiring method | |
| CN110765665B (en) | Dynamic modeling method and system for geography | |
| CN112649860B (en) | Layer velocity inversion method and system based on continuous ant colony algorithm | |
| CN108680968B (en) | Evaluation method and device for seismic exploration data acquisition observation system in complex structural area | |
| CN110726416A (en) | Reinforced learning path planning method based on obstacle area expansion strategy | |
| CN103778657B (en) | Space partition based sound beam tracing method | |
| CN117331119A (en) | Rapid earthquake wave travel time calculation method for tunnel detection | |
| US20150134308A1 (en) | Method and device for acquiring optimization coefficient, and related method and device for simulating wave field | |
| CN103886129B (en) | By the discrete method and apparatus to reservoir grid model of log data | |
| CN103543478A (en) | Geologic morphological interpolation KM (Kriging and Multiple-point geostatistics) method | |
| CN110907990B (en) | Quantitative prediction method and system for post-stack seismic cracks | |
| CN106443829A (en) | Method and apparatus for constructing near-surface model | |
| CN102262245B (en) | Self-adaptive weighted SIRT inversion method and system in seismic tomography processing | |
| CN112649859B (en) | Method and system for establishing seismic wave velocity self-adaptive gridless field node | |
| CN110568497B (en) | Accurate solving method for seismic first-motion wave travel time under complex medium condition |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |