CN102854528B - Pre-stack non-linear inversion method based on particle swarm optimization algorithm - Google Patents

Abstract

The invention relates to a pre-stack non-linear inversion method based on a particle swarm optimization algorithm. Along with the sustained and steady growth of the economy of the country, the need for energy sources is also sustainedly and quickly increased and a serious test of shortage of oil and natural gas resources can inevitably appear. Reservoir prediction is an important research link in oil exploration production; information (contained in a seismic data) about structure, lithology and the like of subsurface formation needs to be fully utilized to find high-quality reservoirs; the information is frequently obtained through seismic inversion; but the conventional inversion is based on a linear and simplified model and is disturbed by a plurality of human factors, and high-precision reservoir parameters frequently cannot be obtained. According to the method, a research on the topic of the particle swarm optimization algorithm is carried out; the Zoeppritz equation is deduced and simplified on the basis of the idea and the principle of the particle swarm optimization algorithm to ensure that the equation is suitable for the solution of the particle swarm optimization algorithm; the particle swarm optimization algorithm is used for the pre-stack inversion; and the pre-stack inversion is respectively carried out on the two-dimensional theoretical model and a offshore seismic data so as to obtain a satisfactory result.

Description

Particle swarm optimization algorithm prestack nonlinear inversion

Technical field

The present invention relates to exploration of oil and gas field technical field, belong to seismic data inverting category, by particle swarm optimization algorithm, solve nonlinear seismic inversion problem specifically.

Background technology

Utilize seismic data to carry out reservoir lithology identification and fluid prediction is the target that geophysicist lays siege to always.Earthquake prestack inversion is that the theoretical foundation of AVO inverting is famous left Puli thatch (Zoeppritz) system of equations.Ostrander (1984) is in research " bright spot type " sandstone reservoir seismic amplitude characteristic procedure, found the phenomenon of " gas sand reflection amplitude increases with offset distance; water bearing sand reflection amplitude increases and reduces with offset distance ", this phenomenon has been improved the ability that hydro carbons detects greatly, from poststack, guide people's sight line into prestack, indicate the appearance of practical AVO technology.But left Puli thatch equation is very complicated, physical meaning is also indefinite.Geophysicist has carried out a lot of forms to it both at home and abroad research and simplification, obtained a series of approximate formulas.The method that Bortfeld (1961) utilizes zone thickness to go to zero to approach single interface has been calculated the reflection coefficient of plane compressional wave and transmitted wave, has provided the formula of reduction of distinguishing fluid and solid.The impact of the variation that the approximate equation that Aki and Richards proposed in 1980 has stressed to describe P-wave And S speed and density on reflection coefficient.Shuey (1985) approximate equation reflection coefficient is expressed as normal incident near, in, the relational expression of three sums of different incidents far away, reflected more intuitively the relation of amplitude and incident angle.Smith and Gidlow (1987), when gas sand being carried out to AVO analysis, have proposed the weighted overlap-add procedure method that AVO analyzes, and have introduced " fluid factor " concept.Fatti etc. (1994), from Aki and Richards equation, are weakening under the prerequisite of density of earth formations item, and the formula that has proposed comparatively to simplify is for predicting the AVO response of gas sand.Elastic impedance (EI) concept and mathematic(al) representation that Connolly (1999) proposes, greatly enriched the intension of AVO analytical technology, expanded the extension of AVO inversion technique, traditional p-wave impedance inverting has been generalized to a point angle superposition of data.Ozdemir etc. (2001) have proposed to utilize AVO information to carry out elastic parameter Simultaneous Retrieving, are called " joint inversion ".Horse strong wind etc. (2003) have proposed the concept of broad sense elastic wave impedance.In recent years, AVO analyzed and stochastic inverse technical method development, and prestack Simultaneous Retrieving method is increasingly mature, the non-linear inversion technology that adopt in Simultaneous Retrieving method more.The proposition of AVO is only used to improve hydro carbons detectability at first, and today, the development of AVO exceeded this category already, was penetrated into the every field of seismic prospecting.At Crack Detection, pressure prediction, oil reservoir detection of dynamic, petroleum-gas prediction, reservoir heterogeneity, be widely applied aspect describing.

Particle swarm optimization algorithm (Particle Swarm Optimization, PSO), claims again particle swarm optimization, is the one of swarm intelligence algorithm (Swarm Intelligence, SI).To be proposed in nineteen ninety-five by American scholar Kennedy and Eberhart at first.Search strategy in defence, the predation of particle cluster algorithm Shi Shou birds colony inspires and forms.The many biologies of occurring in nature have certain group behavior, and as flock of birds, the shoal of fish etc., although single individuality only has simple rule of conduct in colony, the group behavior of composition is but very complicated.A lot of scientists are studied the group behavior of flock of birds or the shoal of fish, comprise Computer Simulation.Particle swarm optimization algorithm takes a hint and for solving optimization problem, regards the solution of optimization problem as point in search volume from this model, is referred to as " particle ".The motion of each particle is according to " flying experience " optimizing of own and other particle, thereby reaches the object of total space search optimum solution.Many scholars and researchist mutually aspect fusion, have proposed a lot of improved PSO algorithms in parameter selection, topological structure and with other optimized algorithms on the basis of basic PSO algorithm.Kennedy and Eberhart have proposed scale-of-two population (Binary Particle Swarm Optimization in 1997, BPSO) algorithm, this is discrete particle cluster (Discrete PSO, the DPSO) algorithm based on continuous space.Clerc (2000) has proposed TSP-DPSO algorithm for travelling entropy problem (Traveling Salesman Problem, TSP), is the DPSO based on discrete space.Jun Sun etc. (2003) are incorporated into quantum behavior in particle cluster algorithm, have proposed quanta particle swarm optimization (Quantum PSO, QPSO).Within 2004, high hawk has proposed the particle cluster algorithm (SA-PSO) based on simulated annealing.The hybrid algorithm that the people such as P.S.Shelokar proposed based on ant group and particle cluster algorithm in 2007, i.e. PSACO (Particle Swarm Ant Colony Optimization) algorithm.Zhou Yalan etc. (2008) are incorporated into Estimation of Distribution Algorithm thought in PSO algorithm, propose discrete particle cluster (Estimation of Distribution PSO, the EDPSO) algorithm of estimating based on distributing.Lin Dongyi etc. (2008) are incorporated into immunologic mechanism in BPSO algorithm, certain Importance of Attributes based on decision table differential matrix is measured as vaccine pattern, propose a kind of minimal attributes reductions algorithm of optimizing based on immunity particle cluster, be called immune discrete particle cluster algorithm (IPSO algorithm).Particle cluster algorithm development in recent years is very fast, is successfully applied to that function optimizing, neural metwork training, pattern recognition classifier, fuzzy system are controlled and the various fields such as engineering, and a large amount of practical applications prove that it is effective.It has good development prospect, is worth doing further research.

Summary of the invention

The present invention utilizes the high-resolution seismic data of relative amplitude preserved processing to carry out non-linear inversion, has used the reflection coefficient formula upgrading, and used particle cluster algorithm to carry out optimizing in the non-linear process of asking of solution.This measure can reduce personal error, effectively improves the precision of inverting.

First be defined as follows:

In geophysics, the reflection coefficient of continuous function x (t) is defined as,

$\underset{\mathrm{\Δt}\→0}{R_{x\left(t\right)}}=\frac{\mathrm{\Δx}}{2\stackrel{\‾}{x\left(t\right)}}$

Sampling interval with dt is sampled into discrete function by continuous function x (t), so,

$R_x=\frac{x_2-x_1}{x_2+x_1}$

Wherein x can be velocity of longitudinal wave α, shear wave velocity β, density p.Subscript 1 represents the parameter of underlying formation, and subscript 2 represents the parameter of overlying strata.Definition transverse and longitudinal wave velocity ratio

$K=\frac{{\mathrm{\β}}_1}{{\mathrm{\α}}_1}$

This patent is used the reflection angular domain CRP gather that obtains of relative amplitude preserved processing as input, realizes the technical scheme that above-mentioned purpose takes as follows:

Step 1: for underground optional position.The span of first given velocity of longitudinal wave reflection coefficient, shear wave velocity reflection coefficient, density reflection coefficient and transverse and longitudinal wave velocity ratio, the size of setting maximum iteration time and population.Given inertia weight and the study factor, set tolerance scope.

Wherein inertia weight, the study factor, Population Size have following meaning:

Inertia weight: the balance of exploration ability and development ability is an importance that affects optimized algorithm performance.For particle swarm optimization algorithm, these two kinds of energy equilibrium of forces are realized by inertia weight w.Larger inertia weight makes particle have larger speed in own original direction, thereby it is farther in former direction, to fly, and has better exploration ability; Less inertia weight makes particle inherit the speed of less former direction, thereby flight is nearer, has better development ability.

The study factor: the study factor is nonnegative constant, represents the weights of particle preference, make particle have the ability that oneself sums up and learn to excellent individual in colony, thereby in colony or in neighborhood, optimum point is close.Kennedy thinks, learns factor sums for two and should be 4.0 left and right, and search effect is now relatively good.Common way is that they are made as to 2.05.

Group size: when group size set very little time, the possibility that is absorbed in local better solution is very high; When group size is for the moment, particle swarm optimization algorithm becomes the method based on individuality search, once it is excellent to be absorbed in office, can not jump out; When group size is very large, the optimization ability of particle swarm optimization algorithm is better, but can cause significantly increase computing time, and when number of groups rises to certain level, continues growth and will no longer include significant effect.As for how many particles actually, participate in search and can obtain desirable effect, forefathers, by using multiple reference functions to calculate its average fitness to the different population of quantity, think that when population quantity remains on 30 left and right, search efficiency is better.

Step 2: velocity of longitudinal wave reflection coefficient, shear wave velocity reflection coefficient, density reflection coefficient and the transverse and longitudinal wave velocity to each particle is than carrying out initialization by random number, and to set initial velocity be zero, make the individual optimal value of particle identical with initial value, and calculate the fitness of each particle, selecting the particle of fitness minimum is current globally optimal solution.

The computing method of particle fitness are:

First define matrix

$M^{+}=\left[\begin{array}{cccc}D& \sqrt{1-{\left(\mathrm{KD}\right)}^2}& -\mathrm{AD}& \sqrt{1-{\left(\mathrm{KBD}\right)}^2}\\ \sqrt{1-D^2}& -\mathrm{KD}& \sqrt{1-{\left(\mathrm{AD}\right)}^2}& \mathrm{KBD}\\ 2D\sqrt{1-{\mathrm{<figure-callout\; id="2"\; label="D"\; filenames="HSA00000748356000012.png"\; state="\{\{state\}\}">D</figure-callout>}}^2}& \frac{1-2{\left(\mathrm{KD}\right)}^2}{K}& -2{\mathrm{CB}}^{2}D\sqrt{1-{\left(\mathrm{AD}\right)}^2}& \frac{\mathrm{BC}[1-2{\left(\mathrm{KBD}\right)}^2]}{K}\\ 1-2{\left(\mathrm{KD}\right)}^2& -2K^{2}D\sqrt{1-{\left(\mathrm{KD}\right)}^2}& -\mathrm{AC}[1-{\left(\mathrm{KD}\right)}^2]& -2{\mathrm{<figure-callout\; id="2"\; label="CB"\; filenames="HSA00000748356000012.png"\; state="\{\{state\}\}">CB</figure-callout>}}^2{\mathrm{<figure-callout\; id="2"\; label="DK"\; filenames="HSA00000748356000012.png"\; state="\{\{state\}\}">DK</figure-callout>}}^2\sqrt{1-{\left(\mathrm{KBD}\right)}^2}\end{array}\right]$

A, B, C, D and K in formula define with the following methods

$A=\frac{{\mathrm{\&alpha;}}_2}{{\mathrm{\&alpha;}}_1}=\frac{1+R_{\mathrm{\&alpha;}}}{1-R_{\mathrm{\&alpha;}}}$

$B=\frac{{\mathrm{\&beta;}}_2}{{\mathrm{\&beta;}}_1}=\frac{1+R_{\mathrm{\&beta;}}}{1-R_{\mathrm{\&beta;}}}$

$C=\frac{{\mathrm{\&rho;}}_2}{{\mathrm{\&rho;}}_1}=\frac{1+R_{\mathrm{\&rho;}}}{1-R_{\mathrm{\&rho;}}}$

D＝sini ₁

$K=\frac{{\mathrm{\&beta;}}_1}{{\mathrm{\&alpha;}}_1}$

Wherein i ₁for incident angle.Define again

$M^{-}=\left[\begin{array}{cccc}-D& \sqrt{1-{\left(\mathrm{KD}\right)}^2}& -\mathrm{AD}& \sqrt{1-{\left(\mathrm{KBD}\right)}^2}\\ \sqrt{1-D^2}& -\mathrm{KD}& \sqrt{1-{\left(\mathrm{AD}\right)}^2}& \mathrm{KBD}\\ 2D\sqrt{{1-\mathrm{<figure-callout\; id="2"\; label="D"\; filenames="HSA00000748356000012.png"\; state="\{\{state\}\}">D</figure-callout>}}^2}& \frac{1-2{\left(\mathrm{KD}\right)}^2}{K}& -2{\mathrm{CB}}^{2}D\sqrt{1-{\left(\mathrm{AD}\right)}^2}& \frac{\mathrm{BC}[1-2{\left(\mathrm{KBD}\right)}^2]}{K}\\ 2{\left(\mathrm{KD}\right)}^2-1& -2K^{2}D\sqrt{1-{\left(\mathrm{KD}\right)}^2}& -\mathrm{AC}[1-{\left(\mathrm{KD}\right)}^2]& -2{\mathrm{<figure-callout\; id="2"\; label="CB"\; filenames="HSA00000748356000012.png"\; state="\{\{state\}\}">CB</figure-callout>}}^2{\mathrm{<figure-callout\; id="2"\; label="DK"\; filenames="HSA00000748356000012.png"\; state="\{\{state\}\}">DK</figure-callout>}}^2\sqrt{1-{\left(\mathrm{KBD}\right)}^2}\end{array}\right]$

By left Puli thatch equation, can be obtained the solution of following form

$R_{\mathrm{pp}}=\frac{\left|M^{-}\right|}{\left|M^{+}\right|}$

The expression formula that can be obtained particle fitness by this expression-form is

E＝||w(θ)*R _pp(R _α，R _β，R _ρ，K，θ)-d(θ)||

In formula, E is particle fitness, and w (θ) is seismic wavelet, and d (θ) is seismological observation data.

Step 3: whether the fitness of more current global optimum has reached error permissible range.If reach error permissible range, just stop calculating, output global optimum is inversion result.Otherwise enter step 4.

Step 4: the more speed of new particle, and with the more position of new particle of new particle rapidity.Iterations increase once, if iterations exceedes maximum iteration time, stops calculating output global optimum.

The update method of particle rapidity and position is:

v _id(t+1)＝wv _id(t)+c ₁r ₁(p _id(t)-x _id(t))+c ₂r ₂(g _d(t)-x _id(t))

x _id(t+1)＝x _id(t)+v _id(t+1)

In formula, v is speed, and x is particle position.M is Population Size, the dimension that n is particle, and 1≤i≤m, 1≤d≤n, w is called the inertia weight factor; c ₁, c ₂be called the study factor or speedup factor, wherein, c ₁for regulating particle to fly to the step-length of self desired positions direction, c ₂for regulating particle to fly to the step-length of overall desired positions; r ₁, r ₂for the random number in [0,1]; T is current iteration algebraically.

Step 5: recalculate the fitness of particle, whether the fitness that detects each particle is less than the fitness before renewal.If be less than, more the individuality of new particle is optimum is the position after particle renewal.

Step 6: newer individuality optimum, select the individuality optimum of fitness minimum as global optimum.Judge whether global optimum reaches error permissible range, reaching error permissible range, just to export global optimum be inversion result, otherwise get back to step 4.Until handle underground all calculation levels.

Above embodiment is to illustrate the invention and not to limit the present invention.

Accompanying drawing explanation

Fig. 1 is the schematic diagram that particle cluster algorithm particle position upgrades.

Fig. 2 is the CRP gather of a synthetic reflection angular domain.

Fig. 3 is that road shown in Fig. 2 collection is used longitudinal wave reflection coefficient that the inverting of this method obtains and the contrast of model longitudinal wave reflection coefficient.

Fig. 4 (a) is the velocity of longitudinal wave that obtains of the inverting of the collection this method of road shown in Fig. 2 and the contrast of model velocity of longitudinal wave;

Fig. 4 (b) is the shear wave velocity that obtains of the inverting of the collection this method of road shown in Fig. 2 and the contrast of model shear wave velocity;

Fig. 4 (c) is the density that obtains of the inverting of the collection this method of road shown in Fig. 2 and the contrast of model density;

Fig. 4 (d) be the transverse and longitudinal wave velocity that obtains of the inverting of the collection this method of road shown in Fig. 2 than with the contrast of model transverse and longitudinal wave velocity ratio.

Fig. 5 is that the p-and s-wave velocity that obtains than (right side) and business software of this method p-and s-wave velocity that inverting obtains in the real data of somewhere is than the contrast on (left side).

Embodiment

By a synthetic road collection, illustrate:

First, utilize left Puli thatch equation just to drill this model, every 3 °, calculate the reflection coefficient on each stratum within the scope of 1 ° to 30 °.10 reflection coefficient sequences that calculate carry out convolution with the zero phase Ricker wavelet of 40Hz respectively, the theogram obtaining, and Ji Jiao road collection data are as shown in Figure 2.

Zhe Shigejiao road collection data are carried out to inverting as input.To solve the longitudinal wave reflection coefficients R calculating _ppand between real seismic record, error function E minimum is as objective function.Adopt the step of particle swarm optimization algorithm prestack inversion as follows:

(1) will obtain Jiao road collection by forward modeling records as observed reading (true value).

(2) utilize the reflection coefficient sequence of particle swarm optimization algorithm generation Ge Jiao road collection as r (t) candidate solution.The position vector of each particle is r (t) candidate solution, and the dimension of each particle is 4, i.e. 4 independent variable longitudinal wave reflection coefficients R _α, transverse wave reflection coefficients R _β, density reflection R _ρ, and transverse and longitudinal wave velocity compares K.

(3) according to step 2, calculate fitness function E.To calculate R _ppthe collection data comparison of Yu Jiao road, calculates fitness function E.

(4) by the fitness function E substitution globally optimal solution G asking for _best, and according to step 4, each particle is carried out to position renewal.

(5) utilize particle cluster algorithm loop iteration, optimum solution is constantly optimized, until meet end condition (as E is less than certain truncation error, or reach maximum iteration time), now individual optimum consistent with global optimum, all particles converge on a bit, and corresponding position vector is the longitudinal wave reflection coefficients R of requirement _α, transverse wave reflection coefficients R _β, density reflection R _ρ, and transverse and longitudinal wave velocity compares K.

While utilizing in this example the reflection coefficient of road, particle cluster algorithm inverting numerical model angle collection data, parameter value is as follows: error function is less than 1e-6, Population Size m=30, inertial factor w=0.1, study factor c ₁=c ₂=2.04, the maximal value of particle position value is [0.1,0.1,0.1,0.75], and minimum value is [0.1 ,-0.1 ,-0.1,0.35], respectively corresponding longitudinal wave reflection coefficients R _α, transverse wave reflection coefficients R _β, density reflection R _ρ, and transverse and longitudinal wave velocity is than the limited field of K.

The reflection coefficient being finally inversed by and the contrast of former reflection coefficient are as shown in Figure 3.Can find out identical fine of the reflection coefficient (dotted line) being finally inversed by and real reflection coefficient (circle).The reflection R obtaining by inverting _α, R _β, R _ρ, and K, can calculate the elastic parameters (solid line) such as the velocity of longitudinal wave, shear wave velocity, density of each layer, contrasts as shown in Figure 4 with the elastic parameter (dotted line) of model.By Fig. 4, found out, P wave data and density data and truth fit like a glove, just shear wave data and p-and s-wave velocity are than there being minimum error, the reason one that produces these errors is because in refutation process, the levels of precision of shear component is subject to the impact of the factors such as initial value design, rock physics relation, and stability is not as compressional component; The 2nd, because the truncation error of refutation process Computer causes.For error is reduced as far as possible, possible in the situation that, not only use the observed reading of P wave data as inverting, can also use the information of converted shear wave.Even in the situation that having VSP data, can be using the data of reflected P-wave, reflection wave, transmitted P-wave, transmitted shear wave all as the observed reading of inverting, the elastic parameter of the various reflection coefficients that inverting obtains so and the rock of deriving thus will be more accurate.

List of references

Zoeppritz?K，Erdbebenwellen?VIIIB.On?the?reflection?and?propagation?of?seismic?waves.Gottinger?Naehriehten，1919，66-84

Ostrander?W?J.Plane?waves?reflection?coefficient?s?for?gas?sands?at?normal?angles?of?incidence.Geophysics，1984，49(10)：1637-1648.

Bortfeld?R.Approximation?to?the?reflection?and?transmission?coefficients?of?plane?longitudinal?and?transverse?wave.Geophysics?Prospecting，1961，9(4)：485-502.

Aki?K，Richards?P?G.Quantitative?seismology：Theory?and?methods.W.H.Freeman?and?Co，1980，100-170.

Shuey.A?simplification?of?the?Zoeppritz?equations.Geophysics，1985，50(4)：609-634.

Smith?G?C，Gidlow?P?M.Weighted?stacking?for?rock?property?estimation?and?detection?of?gas.Geophysical?Prospecting，1987，35：993-1014.

Fatti?J?L，Vail?P?J，Smith?G?C，et?al.Detection?of?gas?in?sandstone?reservoirs?using?AVO?analysis：A?case?seismic?case?history?using?the?Geostack?technique.Geophysics，1994，59(5)：1362～1376.

Connolly?P.Elastic?Impedance.The?Leading?Edge，1999，18(4)：438～452.

Ostrander?W?J.Plane-wave?reflection?coefficients?for?gas?sands?at?non-normal?angels?of?incidence.Exploration?Geophysics，1984，49：1637-1648.

Horse strong wind. the FORWARD AND INVERSE PROBLEMS of broad sense elastic impedance in seismic prospecting. Chinese Journal of Geophysics, 2003,46 (1): 118-124.

Du Yan .AVO technology and the application in Soviet Union's Sulige gas field reservoir gas-bearing property detects: (master thesis). Chengdu. Chengdu University of Technology, 2009.

Kennedy，Eberhart?R?C.Particle?Swarm?Optimization.International?conference?on?Neural?Networks，1995，1942-1948.

Liu Bo. particle swarm optimization algorithm and engineering application thereof. Beijing: Electronic Industry Press, 2010.

Kennedy?J，Eberhart?R?C.A?discrete?binary?version?of?the?Particle?Swarm?Algorithm.The?1997?Conference?on?System，Cybernetics?and?Informatics，1997，4104-4109.

Clerc?M.Discrete?Particle?Swarm?Optimization，Illustrated?by?the?Traveling?Salesman?Problem.http://www.mau-ricecierc.net，2000.

Sun?J，Feng?B，Xu?W.Particle?Swarm?Optimization?with?Particle?Having?Quantum?Behavior.IEEE?Proc.Of?Congress?on?Evolutionar?Computation，2004，325-331.

Gao Ying, thanks to triumph. based on the particle swarm optimization algorithm of simulated annealing. and computer engineering and application, 2004 (1): 47-50.

Shelokar?P?S，Siarry?P，Jayaraman?V?K，Kulkarni?B?D.Particle?swarm?and?ant?colony?algorithms?hybridized?for?improved?continuous?optimization.Applied?Mathematics?and?Computation，2007，188：129-142.

Zhou Yalan, Wang Jiahai etc. a kind of Discrete Particle Swarm Optimization Algorithm based on distributing and estimating. electronic letters, vol, 2008,36 (6): 1242-1248.

Lin Dongyi, Liao Jiankun. an immune discrete particle cluster algorithm of minimum reduction problem. small-sized microcomputer system, 2008,29 (6): 1089-1092.

He Yi. the inverting of Seismic Reservoir parametrical nonlinearity and Study on Forecasting Method: (doctorate paper). Qingdao. Chinese Marine University, 2008.

Claims (1)

1. a particle swarm optimization algorithm prestack nonlinear inversion, utilize the high-resolution seismic data of relative amplitude preserved processing to carry out non-linear inversion, in the process that solves nonlinear problem, use the reflection coefficient formula upgrading, use particle cluster algorithm to carry out optimizing, the reflection angular domain CRP gather that use relative amplitude preserved processing obtains is as input, and concrete steps are:

Step 1: for underground optional position, the span of first given velocity of longitudinal wave reflection coefficient, shear wave velocity reflection coefficient, density reflection coefficient and transverse and longitudinal wave velocity ratio, set the size of maximum iteration time and population, given weight factor and the study factor, set tolerance scope;

Step 2: velocity of longitudinal wave reflection coefficient, shear wave velocity reflection coefficient, density reflection coefficient and the transverse and longitudinal wave velocity to each particle is than carrying out initialization by random number, and to set initial velocity be zero, make the individual optimal value of particle identical with initial value, and calculate the fitness of each particle, selecting the particle of fitness minimum is current globally optimal solution;

Step 3: whether the fitness of more current global optimum has reached error permissible range, just stops calculating if reach error permissible range, and output global optimum is inversion result, otherwise enters step 4;

Step 4: the more speed of new particle, and with the more position of new particle of new particle rapidity, iterations increase once, if iterations exceedes maximum iteration time, stops calculating output global optimum;

Step 5: recalculate the fitness of particle, whether the fitness that detects each particle is less than the fitness before renewal, if be less than, more the individuality of new particle is optimum is the position after particle renewal;

Step 6: newer individuality optimum, select the individuality optimum of fitness minimum as global optimum, judge whether global optimum reaches error permissible range, reaching error permissible range, just to export global optimum be inversion result, otherwise get back to step 4, until handle underground all calculation levels;

Velocity of longitudinal wave reflection coefficient in each described step, shear wave velocity reflection coefficient, density reflection coefficient, the definition in the following manner of transverse and longitudinal wave velocity ratio:

In geophysics, the reflection coefficient of continuous function x (t) is defined as:

Sampling interval with dt is sampled into discrete function by continuous function x (t):

Wherein x can be velocity of longitudinal wave α, shear wave velocity β, density p; Subscript 1 represents the parameter of underlying formation, and subscript 2 represents the parameter of overlying strata; Definition transverse and longitudinal wave velocity ratio:

It is characterized in that, the computing method of the particle fitness in described step 2,3,5 are as follows:

Definition matrix:

A, B, C, D and K in formula define with the following methods:

Wherein i ₁for incident angle;

Definition matrix:

By left Puli thatch equation, can be obtained the solution of following form:

The expression formula that can be obtained particle fitness by this expression-form is:

In formula, E is particle fitness, and w (θ) is seismic wavelet, and d (θ) is seismological observation data.

Images