Summary of the invention
Technical problem to be solved by the invention is to provide the use for carrying out Inversion Calculation to elastic parameter using earthquake recordIn the improved adaptive GA-IAGA of earthquake data before superposition parametric inversion.
The technical scheme to solve the above technical problems is that
A kind of improved adaptive GA-IAGA for earthquake data before superposition parametric inversion, includes the following steps:
S1, oil field subsurface reservoir well logging region is sampled, obtains at least two well logging sample point datas, obtains at leastTwo groups of elastic parameters；Earthquake record corresponding at least two well logging sampled point to oil field subsurface reservoir well logging regionSampled point carries out seismic prospecting, obtains at least one set of real seismic record；
S2, it is iterated operation with the elastic parameter, updates elastic parameter；
S3, the elastic parameter being calculated by S2 calculate corresponding earthquake record calculated value, and with the reality in the S1Border earthquake record compares, and is verified calculated result；
S4, when verify calculated result be less than threshold value when, be proved to be successful to obtain the inversion result of the elastic parameter；Work as verifyingWhen calculated result is more than or equal to threshold value, authentication failed then returns to S2.
Further, operation is iterated to the elastic parameter by the genetic algorithm that Real-valued encodes in the S2,The numerical value of the gene position of individual in the genetic algorithm in initial population is according to the several groups elastic parameter obtained in the S1Float up and down generation within a preset range；Individual adaptation degree is counted by the verifying calculated result in the genetic algorithmIt calculates.
Further, elastic parameter described in each group includes seismic wave velocity of longitudinal wave Vp, seismic wave shear wave velocity VsAnd rockStone density p.
Further, oil field subsurface reservoir well logging region is sampled in the S1, obtains n+1 well logging samplingThe total n+1 group elastic parameter of point:
[Vpiwell, Vsiwell, ρiwell], (wherein i=1,2 ..., n+1, n are positive integer)；
It is corresponding with the n+1 well logging sampled point to seismic prospecting acquisition is carried out in oil field subsurface reservoir well logging regionN earthquake record sampled point m kind different angle real seismic record it is as follows:
[s(θi,j)] (wherein i=1,2 ..., n；J=1,2 ... m, m are positive integer).
Further, each of described genetic algorithm individual is the Real-valued one-dimension array that length is 3n+3；Initial kindGroup in individual in each element value according to the total 3n+3 elastic parameter of the n+1 group obtained in the S1 within a preset range onLower floating randomly selects.
Further, the value of first three element in the individual of the initial population is according to the of acquisition of logging well in the S1One group of elastic parameter floats up and down and randomly selects；
The value of the 4th element to the last one element in the individual of the initial population is randomly selected by following constraint:
Further, it is calculated and is corresponded to Zoeppritz equation by the S2 elastic parameter being calculated and angleEarthquake record calculated value:
[s'(θi,j)] (wherein i=1,2 ..., n；J=1,2 ... is m).
Further, the calculation of the verifying calculated result is as follows:
Further, by the S2 elastic parameter being calculated and angle with the Aki& in Zoeppritz equationRichard approximate equation calculates corresponding earthquake record calculated value.
The present invention utilizes the theoretical model of pre-stack seismic, establishes the correlativity between elastic parameter and earthquake record, bulletEarthquake record is calculated by forward modeling in property parameter；Well logging obtains initial elastic parameter, with genetic algorithm to elasticityThe elastic parameter newly obtained is constantly carried out forward modeling and compared with earthquake record by the continuous iteration of parameter, the earthquake note for obtaining forward modelingRecord infinitely approaches real seismic record, finally to obtain accurate elastic parameter value.
Specific embodiment
The principle and features of the present invention will be described below with reference to the accompanying drawings, and the given examples are served only to explain the present invention, andIt is non-to be used to limit the scope of the invention.
As shown in Fig. 2, the present invention includes the following steps,
S1, oil field subsurface reservoir well logging region is sampled, obtains at least two well logging sample point datas, obtains at leastTwo groups of elastic parameters；Earthquake record corresponding at least two well logging sampled point to oil field subsurface reservoir well logging regionSampled point carries out seismic prospecting, obtains at least one set of real seismic record；
S2, it is iterated operation with the elastic parameter, updates elastic parameter；
S3, the elastic parameter being calculated by S2 calculate corresponding earthquake record calculated value, and with the reality in the S1Border earthquake record compares, and is verified calculated result；
S4, when verify calculated result be less than threshold value when, be proved to be successful to obtain the inversion result of the elastic parameter；Work as verifyingWhen calculated result is more than or equal to threshold value, authentication failed then returns to S2.
Each step of the invention is described in detail below:
The inversion problem of earthquake data before superposition:
As shown in Figure 1, obtaining n+1 survey by sampling to oil field subsurface reservoir well logging region first in the present inventionThe total n+1 group elastic parameter of well sampled point, wherein can also be calculated by the exploration result or empirical value of related geologic provinceThe initial estimate (and being modified by initial estimate of the real seismic record to elastic parameter) of elastic parameter:
[Vpi, Vsi, ρi], (wherein i=1,2 ..., n+1)；
Seismic prospecting is carried out to oil field subsurface reservoir well logging region and obtains n corresponding with the n+1 group elastic parameterThe real seismic record of a earthquake record sampled point m kind different angle is as follows:
[s(θi,j)] (wherein i=1,2 ..., n；J=1,2 ... is m).
Elastic parameter and real seismic record described in S1 are obtained through the above way.
In several sampled points, (sampled point is generally arranged at different elevations, sampled point to adjacent upper layer and lower layer well logging sampled pointDensity is higher, and data are more accurate) an earthquake record sampled point is corresponded to, therefore n+1 well logging sampled point corresponds to n earthquake recordSampled point, and several different angles can be used in each earthquake record sampled point, when using m kind different angle, then obtain total nmGroup real seismic record.
Spirit of the invention, which is that, is modified elastic parameter by real seismic record, that is, passes through earthquake record pairElastic parameter carries out inverting.Using elastic parameter, there is correspond to each other in Zoeppritz equation with earthquake record for invertingRelationship；The essence of inverting is forward modeling many times, i.e., corresponding earthquake record is calculated by the elastic parameter that fitting obtains, with groundShake record is compared.
Since Zoeppritz equation is complex, the present invention is using the Aki& in Zoeppritz equationRichard approximate equation:
Establish inverting convolution model:
Establishing inverting convolution model is one of the key step for carrying out prestack AVO elastic parameter inversion, establishes convolution modelMainly there are following steps:
First, calculate reflection Rpp.The present invention calculates R using Aki&Rechard approximate equationpp:
Wherein, Δ Vp, Δ Vs, Δ ρ respectively indicates adjacent upper layer and lower layer Vp、VsDifference between ρ, andWithIndicate upper and lower level Vp、VsWith the average value of ρ, θ is angle,It is calculated according to real data, it can according to the formulaTo obtain Rpp, one-component as earthquake record convolution operation.
Second, obtain seismic wavelet.Seismic wavelet is another component of seismic convolution model, by by wavelet withReflection coefficient does convolution operation and obtains earthquake record, suitable for establishing forward model and production Synthetic seismic record, the present inventionRicker wavelet is used, is a kind of seismic wavelet of zero phase, expression formula is as follows:
WhereinVmIndicate dominant frequency, t indicates the time, can be with manual setting.Usually extracted according to real seismic record, it is one group of definite value.
Reflection coefficient and Ricker wavelet are carried out convolution operation by third, and equation is as follows:
S (θ)=Rpp(θ)*f(t)+n(t)
Wherein f (t) indicates that seismic wavelet, n (t) indicate noise.Noise factor is not considered in the present invention.Calculated s(θ) is used to construct objective function.
Here RppDuring carrying out convolution operation with Ricker wavelet, all angles are individually calculated, not associated.As shown in Figure 3, it can be seen that the earthquake record of 8 angles.
Since the core of the calculating of entire earthquake record is by calculating reflection Rpp, and every layer of sampled point of each angleCorresponding RppCalculating be that the difference for subtracting upper layer data by the average value and lower layer of three elastic parameters of upper layer and lower layer is asked, and average value is equal to the difference that upper layer numerical value adds 1/2, it can thus be seen that RppIt can be joined completely by upper layer elasticitySeveral and difference indicates, and can be seen that difference occupies leading position in formula.
By the above-mentioned process for establishing inverting convolution model, forward modeling is completed, that is, completes and passes through bullet described in S3The process of corresponding earthquake record calculated value can be calculated in property parameter:
[s'(θi,j)] (wherein i=1,2 ..., n；J=1,2 ... is m)
The elementary tactics of genetic algorithm:
Interative computation in S2 generallys use genetic algorithm and is iterated update to elastic parameter.Each of genetic algorithmIndividual is adopted as the Real-valued one-dimension array that length is 3n+3；The value of each element is according in individual in initial populationThe total 3n+3 elastic parameter of the n+1 group obtained in S1, which floats up and down, to be randomly selected.And individual adaptation degree is then logical in genetic algorithmIt crosses Zoeppritz equation and calculates the identical property of new elastic parameter and real seismic record that genetic algorithm iteration obtains to countIt calculates.
With earthquake record sampled point n=240, for angle m=8, existing genetic algorithm carries out just in the following wayThe value of beginning population at individual, each individual is in following form:
Gi=(Vp1,Vs1,ρ1,Vp2,Vs2,ρ2,L,Vpn+1,Vsn+1,ρn+1), n=1,2, L, 243
The conventional initialization strategy of genetic algorithm:
Each element, that is, each gene position is using following constraint in initial individuals:
0.8gVpwell≤Vp≤1.2gVpwell
0.8gVswell≤Vs≤1.2gVswell
0.9gρwell≤ρ≤1.1gρwell
Due to generating elastic parameter in the threshold range at random, it is possible that, upper layer elastic parameter is most left in sectionSide (is minimized), and lower layer's elastic parameter (is maximized) in the section rightmost side, it is thus possible to will lead to bilevelThe difference error that elastic parameter generates is larger, and bigger by left and right shock range, leads to calculated RppIt is inaccurate,And then cause calculated earthquake record error larger.
The improved initialization strategy of genetic algorithm:
Since traditional initialization strategy there is a problem of inaccurate, received by improving initialization operation to accelerate algorithmIt holds back, initialization operation is the threshold value that a series of elastic parameters obtained according to log data acquire, and is selected at random in threshold rangeIt takes, according to a series of difference of the available upper and lower level elastic parameters of log data, then the difference of elastic parameter is setThreshold value randomly selects within the scope of difference threshold, as follows:
The value of first three element in genetic algorithm in the individual of initial population is according to logging well the first of acquisition in the S1Group elastic parameter, which floats up and down, to be randomly selected；
The 4th element in the individual of initial population introduces the bullet of adjacent well logging sampled point to the value of the last one elementDifference DELTA Vp between property parameteri、ΔVsi, Δ ρi, and randomly selected by following constraint:
Fitness calculates:
To the fitness of each individual in genetic algorithm, by by the value of the elastic parameter of the representative of each element in the individualForward modeling is carried out, calculates corresponding earthquake record, and to several groups earthquake record compared with real seismic record, verifying is calculated and calculatesAs a result it is used as fitness, while being also verified meter to be compared earthquake record calculated value with real seismic record in S3Calculate the process of result:
The calculation for verifying calculated result is as follows:
Crossover operation:
Because the present invention is to use single point crossing method merely, the search capability of algorithm is weaker, institute based on real codingWith the present invention using a kind of crossover operator for being more suitable real coding, using arithmetic crossover strategy, the strategy is by twoThe linear combination of body, generates two new individuals, and λ is the random distribution number in [0,1] section.Expression formula is as follows:
Child1=λ × parent2+ (1- λ) × parent1
Child2=λ × parent1+ (1- λ) × parent2
Assuming that λ is 0.3, a gene position in parent individuality intersected, then the corresponding gene position of offspring individualValue is as shown in Figure 4.
Mutation operation:
A kind of non-homogeneous, adaptive mutation operation suitable for the research topic is used.Select needs at random firstThe gene v of variation, value range are [bound1,bound2], when carrying out mutation operation, it is assumed that current evolutionary generation is pi, kindThe maximum evolutionary generation of group is p, which is given before algorithm starts, it is assumed that v1, v2It is respectively as follows:
v1=v-bound1
v2=bound2-v
Variation starts, and firstly generates a random number λ:
As λ > 0.5, following operation is executed to the gene to make a variation:
vnew=v+ Δ v
As λ < 0.5, following operation is executed to the gene to make a variation:
vnew=v- Δ v
Genetic algorithm brief summary:
In experiment, the initialization of population of most basic genetic algorithm (abbreviation BGA here) selects conventional measures, selection operationUsing roulette selection, crossover operation is using single point crossing, and mutation operation is using non-homogeneous, adaptive changeETTHER-OR operation.The initialization of population of Revised genetic algorithum selects improvement strategy, and selection operation is intersected using algorithm of tournament selectionOperation is using arithmetic crossover, and mutation operation is using non-homogeneous, adaptive mutation operation (abbreviation IGA here).ThisIn two different initialization operations be respectively adopted to two kinds of algorithms test, be briefly referred to as BGA and IGA, algorithm hereRelative parameters setting is as shown in the table:
1 genetic algorithm general parameter of table
| Parameter | Parameter declaration | Parameter size |
| Pc | Crossover probability | 0.7 |
| Pm | Mutation probability | 0.05 |
| popsize | Population Size | 40 |
| IterNum | The number of iterations | 5000 |
Borehole log data in data set is the data for having 241 sampled points, including velocity of longitudinal wave Vp, shear wave velocity VsAnd density p.The corresponding 8 different angles of each sampled point: [0 °, 6 °, 11 °, 17 °, 23 °, 29 °, 34 °, 40 °], each dataCollection all uses this 8 angles.Forward modeling is carried out to log theoretical model using Aki&Richard formula, utilizes logModel calculates reflection coefficient, then reflection coefficient and wavelet are carried out convolution, needs to utilize upper and lower two because generating earthquake recordRelationship between group sampled point, therefore include 240*8 data in earthquake record.Using it is proposed that Revised genetic algorithum with mostThe genetic algorithm on basis carries out earthquake data before superposition parametric inversion respectively, is tested, can must be tested according to algorithm parameter settingAs a result such as following table.
The comparison of 2 three parameter average fitness value numerical value of table
| Algorithm title | BGA | IGA |
| Average fitness value | 0.006332 | 0.001231 |
The comparison of 3 three parameter average correlation coefficient numerical value of table
| BGA | IGA |
| Vp | 0.556364 | 0.941373 |
| Vs | 0.650805 | 0.915569 |
| ρ | 0.462346 | 0.949102 |
Proof of algorithm:
Due to reflection coefficientRppCalculating be by elastic parameter Vp、VsIt obtains, and is finally calculated same with ρGroup inverting seismic data is by unlimited kind of Vp、VsIt combines to obtain with ρ.Accordingly, there exist three elastic parameters, to have error to be calculated anti-Drilling seismic data and two parameter does not have error, and the inverting seismic data that a parameter has that error calculation obtains is identical.In order to preferably evaluate optimization algorithm to the quality of prestack AVO elastic parameter inversion result, present invention employs Pearson product-momentsRelated coefficient (Pearson product-moment correlation coefficient, also referred to as PPMCC or PCCs) comesMeasure the correlation circumstance of three elastic parameters and actual three elastic parameter that are finally inversed by.Here the correlation coefficient function that we establishIt is as follows:
Wherein, XiFor some parameter standard value of three elastic parameters, YiFor corresponding inverting value,Respectively oneThe average value of class value.
Since seismic data solution procedure is sufficiently complex, target function value is smaller, and the related coefficient of three elastic parameters is differentIt is fixed higher, also, the related coefficient of three elastic parameters is higher, and target function value is equally also not necessarily smaller, still, when target letterWhen number reaches theoretially optimum value and is 0, Vp、VsTheoretially optimum value 1 can be reached with the related coefficient of ρ.Therefore, of the inventionJudge the superiority and inferiority of inversion result jointly by way of combining target function value and related coefficient.Final goal is can to makeThe obtained target function value of inverting is small, while also allowing the related coefficient of three elastic parameters high.
From correlation data as can be seen that genetic algorithm gradually improves after a series of improvement, the elasticity being not only finally inversed byParameter more levels off to practical logging data, and the seismic data being finally inversed by also more coincide with actual seismic data, while elasticityThe related coefficient of parameter is also higher.
The foregoing is merely presently preferred embodiments of the present invention, is not intended to limit the invention, it is all in spirit of the invention andWithin principle, any modification, equivalent replacement, improvement and so on be should all be included in the protection scope of the present invention.