CN116879946B - Pre-stack inversion method for improving particle filtering algorithm based on multi-objective dayfish algorithm - Google Patents

Abstract

The invention relates to a particle filter algorithm improved pre-stack inversion method based on a multi-objective dayfish algorithm, which comprises the following steps: establishing an initial model; assigning the initial model as a prior probability model to the particles according to preset particle filtering algorithm parameters; carrying out a particle filtering algorithm prediction process from the seismic trace at the first incident angle; performing iteration by adopting a multi-objective dayf algorithm, and finishing the updating process of the particle filtering algorithm during the iteration to obtain a final Pareto set under the incident angle; assigning the position of the f in the final Pareto set obtained by the incidence angle to the particles, and taking the position of the f as a posterior probability model for inverting the incidence angle and taking the position of the f as an prior probability model for inverting the seismic record of the observation of the next incidence angle; and obtaining a global final Pareto set after calculating the observation seismic records of all incidence angles, and assigning the f-positions on the global final Pareto set as particles to obtain an inversion result.

Description

Pre-stack inversion method for improving particle filtering algorithm based on multi-objective dayfish algorithm

Technical Field

The invention relates to the technical fields of geophysics, seismic exploration and oil and gas exploration, in particular to a pre-stack inversion method for improving a particle filtering algorithm based on a multi-target Alfa algorithm.

Background

Inversion of pre-stack AVA (Amplitude versus angle, amplitude with incident angle transformation) is a conventional method for obtaining elastic parameters of an underground medium in the technical field of seismic exploration, and the elastic parameters mainly comprise longitudinal wave speed, transverse wave speed and density. The longitudinal wave speed, the transverse wave speed and the density are important reservoir parameters, the accurate inversion of the parameters is the important research content of conventional prestack AVA inversion, and the accurate three parameters can not only effectively improve the identification rate of a high-quality reservoir, but also reduce the exploration and development cost.

Currently, pre-stack AVA inversion methods are based mainly on various approximate Zoeppritz equations and accurate Zoeppritz equations. Based on the approximate zopriz equation, the derivatives of Shuey, russell and Aki-Richards et al predominate, but are difficult to apply to large angle incidence seismic gathers, thus causing instability in the inversion of the density parameters. Based on an accurate zopriz equation, the longitudinal and transverse wave reflection coefficient and the longitudinal and transverse wave transmission coefficient can be calculated simultaneously, and the nonlinear equation is extremely complex and difficult to understand; in order to solve the zopriz equation more scientifically to obtain an accurate elastic coefficient, attempts are made to solve the accurate zopriz equation by using various unconstrained optimization methods such as local nonlinearity, global nonlinearity and the like.

The local nonlinear solving method which has been used at present mainly comprises the following steps: gauss Newton (Gauss-Newton) iterative method, steepest descent method, conjugate gradient method, and damped least squares (LM), etc.; the global nonlinear solving method which is used at present and has higher inversion precision is mainly an intelligent optimization algorithm which has better performance at present, such as: genetic algorithms (Genetic algorithm, GA), particle swarm optimization algorithms (Particle swarm optimization, PSO), and whale optimization algorithms (Whale optimization algorithm, WOA), and the like.

Both of the above solutions have some problems: for the local nonlinear solving method, although the method does not depend on accurate prior model distribution and inversion results can certainly show the characteristics of reservoirs, the inversion accuracy is low in ubiquitous and noise immunity is not achieved under the condition that a Bayesian framework is not introduced. The damping least square method is used as a representative, statistical correlation among model parameters and prior distribution characteristics obeyed by the parameters are not considered, so that the accuracy of an inversion result is lower than that obtained by a Bayesian inversion method introducing prior information. And the second partial derivative is needed to be calculated when solving the objective function by taking the Gaussian Newton iteration method as a representative, but when the nonlinear forward equation is complex in form and strong in nonlinearity, the second partial derivative can be very difficult and unstable to solve, and omitting the second partial derivative can cause the whole Gaussian Newton iteration method to be not converged.

For the global nonlinear solving method, although the solving precision is very high, the accuracy of the prior model is very dependent, if the prior model is inaccurate, the algorithm is usually easy to fall into a local extremum, and finally the result of the inversion of the severe image is obtained. Aiming at the problems that the global nonlinear solving method has higher dependence on the prior model, a Markov Monte Carlo (Markov chain Monte Carlo, mcMC) method is generally adopted to generate a reasonable prior model, the typical representation of which is Particle Filter (PF), and a plurality of pre-stack AVA inversion methods based on the PF algorithm are developed by a plurality of scholars at present. The particle filtering algorithm has the greatest defect that the particle in the resampling stage is poor and degenerated, the weight of most particles is 0, and the result of serious image inversion is shown, and various improved resampling algorithms are proposed by a plurality of scholars aiming at the problem, such as using a generalized regression neural network to adjust a particle sample, so that the particle sample is closer to the posterior probability density, and the filtering precision is improved; combining the unscented Kalman filtering and the weight optimization method, optimizing the particle state by taking the unscented Kalman filtering as an importance density function, guiding particle sampling, and introducing an influence factor to optimize the weight of the particles; the particle weight in the particle filtering algorithm is optimized by the single-target group intelligent optimization algorithm, the particle weight is taken as an objective function, the process can be analogous to posterior probability distribution that particles are continuously approaching to the real system state, the particles in the particle filtering algorithm are in a high-level likelihood domain through the strong global optimization capability of the group intelligent optimization algorithm, the diversity of the particles in the particle filtering algorithm is ensured, the inversion precision is improved, and the pre-stack AVA inversion method for improving the particle filtering algorithm for the whale optimization algorithm is adopted at present and has better effect.

Disclosure of Invention

The invention aims to realize accurate inversion of three parameters of elasticity based on a particle filtering algorithm improved by a multi-objective dayf algorithm, provide powerful support for accurate identification of a high-quality reservoir, and provide a pre-stack inversion method based on the particle filtering algorithm improved by the multi-objective dayf algorithm.

In order to achieve the above object, the embodiment of the present invention provides the following technical solutions:

a method of pre-stack inversion based on a multi-objective, dayfish algorithm to improve a particle filtering algorithm, comprising the steps of:

step 1, collecting logging data and establishing an initial model;

step 2, according to preset particle filtering algorithm parameters, assigning an initial model as a global prior probability model to particles;

step 3, starting from the seismic trace at the first incident angle, performing a particle filtering algorithm prediction process;

step 4, iteration is carried out by adopting a multi-objective dayf algorithm, the updating process of the particle filtering algorithm is completed in the period, and a final Pareto set under the incident angle is obtained after the iteration is carried out for a plurality of times;

step 5, assigning the f-positions in the final Pareto set obtained by the incidence angle to the particles, and using the f-positions as a posterior probability model for inverting the incidence angle and a priori probability model for inverting the seismic record of the observation of the next incidence angle;

And 6, repeating the steps 3 to 5 until the final global Pareto set is obtained after the observation seismic records of all the incidence angles are calculated, assigning the positions of the head on the final global Pareto set as particles, and calculating expected values of all the head to obtain an inversion result.

The step 1 specifically comprises the following steps: the collected logging data comprises longitudinal wave velocity logging data, transverse wave velocity logging data and density logging data; and filtering out high-frequency signals in the logging data by a linear smoothing method, and establishing an initial model.

The step 2 specifically comprises the following steps: the particle filter algorithm parameters include particle number, observed noise variance, and process noise variance.

The step 3 specifically comprises the following steps:

establishing a nonlinear dynamic equation:

wherein m is _θ-1 Representing the calculated theta-1 deg. angleA posterior probability model after completion; m is m _θ Representing the model m when calculating the angle θ° _θ-1 A priori probability model after adding process noise; f () represents a state transition equation; w (w) _θ Representing a process noise variance of the system; z _θ Representing an observed seismic record of the system at an angle of θ°; h () represents the observation equation; v _θ Representing observation noise of the system when calculating the theta angle;

Model parameters at a given angle θ -1 momentPredicting model parameters at a given angle theta>Is predicted by the following equation:

wherein,a posterior probability model at the time of theta-1 DEG; />The prior probability model at the time theta is adopted;a priori probability from 3 DEG to theta-1 DEG; d, d _obs(3:θ-1) Representing a set of observation seismic records from 3 ° to θ -1 ° time instants; />Is shown in event->Event under Condition->Probability of occurrence; (m) _θ-1 |d _obs(3:θ-1) ) Representing model parameters m _θ-1 In the seismic observation record set d _obs(3:θ-1) Distribution below.

In step 4, before iterating with the multi-objective dayf algorithm, the method further comprises the steps of:

initializing parameters of a multi-objective dayfish algorithm, setting a male dayfish number q=100, a female dayfish number p=100, a male dayfish mating number nc=100, a offspring dayfish mutation number nm=50, a gravity coefficient g=0.8, and a positive attraction constant a ₁ Positive attraction constant a =1 ₂ =1.5, mutation coefficient mu=0.02, modified visibility coefficient beta=2, random flight coefficient fl=0.77, wedding Dance coefficient dance=0.77, wedding Dance coefficient decrease rate dd=0.99, random flight coefficient decrease rate fd=0.99, pareto aggregate capacity np=100, maximum number of iterations t _max ＝100。

The step 4 specifically comprises the following steps:

Determining the number of the dayfish population, and initializing the positions of the male and female dayfles.

M＝{x _i |x _i ＝(x _i1 ,x _i2 ,...,x _id ),i＝1,2,...,q}

F＝{y _i |y _i ＝(y _i1 ,y _i2 ,...,y _id ),i＝1,2,...,p}

Wherein M represents a population of male dayfish; f represents a female population; x is x _i Representing the position of the ith male; y is _i Representing the position of the ith female; q represents the number of male dayfish (q=100); p represents the female dayfish number (p=100); d represents the total number of dimensions (d=375); x is x _id Representing the position of the ith male dimension; y is _id Representing the position of the ith female f iotaeld d dimension;

the prior probability model is known only at the θ timeUpdating posterior probability model->The process of (1) is as follows:

wherein,the posterior probability at time θ; d, d _obs(θ) The observation seismic record corresponding to the time theta is obtained; />The prior probability at the time theta is given; />Likelihood functions for observing seismic records; p (d) _obs(θ) ) To observe seismic record d _obs(θ) Can be expressed as 1 constant;

assuming that the likelihood function satisfies the gaussian distribution, there are:

wherein C is _n A variance matrix representing the actual noise,is C _n In the present algorithm, an identity matrix I may be multiplied by the observed noise v at time θ _θ To represent; />Representing model parameters +.>Is used for the co-variance matrix of (a),is->An inverse matrix of (a); />Representing model parameters +.>Is a mean vector of (a); t represents the transpose of the matrix; / >Representing the synthetic seismic record corresponding to the moment theta, namely z in a nonlinear dynamic equation established by a particle filtering algorithm _θ Expansion by Taylor's formula>And omit higher-order items above second order, update

Minimizing the objective function J ₁ (m _θ )：

Minimizing the objective function J ₂ (m _θ )：

Where Δm represents the model delta, i.e., the current modelIs +.>Is a difference in (2);representing composite seismic records->For model parameters->Taking the first partial derivative, expressed as:

wherein W represents a longitudinal wave matrix, K represents a Accord ratio matrix of longitudinal wave coefficients, and the Accord ratio matrix of reflection coefficients for an incident angle θ° is:

wherein R is _N Representing the reflection coefficient of the N-th dimension, V _PN 、V _SN And ρ _N Representing the transverse wave speed, the longitudinal wave speed and the density corresponding to the Nth dimension;

the reflectance terms in the Accord matrix K can be obtained by the exact Zoeppritz equation, which can be expressed as:

AR＝B

wherein,

wherein V is _P1 、V _P2 Representing longitudinal wave velocity at both ends of the medium, V _S1 、V _S2 Representing transverse wave velocity, ρ, at both ends of the medium ₁ 、ρ ₂ Representing the density of both ends of the medium; r is R _PP 、R _PS 、T _PP 、T _PS Respectively representing the longitudinal wave reflection coefficient, the longitudinal wave transmission coefficient, the transverse wave reflection coefficient and the transverse wave transmission coefficient;

the partial derivative of the reflection coefficient with respect to the model parameters can be expressed as:

Wherein ψ= [ ψ ] ₁ ,ψ ₂ ,ψ ₃ ,ψ ₄ ]Is the proportion of elastic parameters, and and->The specific parameters of (a) are as follows:

the 6 elastic parameters between the single reflection coefficient and the reflection interface can be obtained by the chain rule:

updating the current position of the 1 st male-f to be the initial position, and performing an updating process on the 1 st male-f to obtain a current objective function value of (J) ₁ (m _θ ),J ₂ (m _θ ) A) is provided; and determining the current objective function values of all the remaining male and female dayfles in the same manner;

combining all male and all female individuals, performing one-time Pareto non-dominant sorting, taking Np individuals before sorting as a Pareto set, and selecting the dayf individuals with the rank of 1 in the Pareto set as the head of the male and female individuals;

starting a first iteration:

the current position of female dayfish 1 is calculated from the initial stage, and the current objective function value (J ₁ (m _θ ),J ₂ (m _θ ) Updating the current position and the current objective function value of the female-figet 1; calculating the current positions and the current objective function values of all the remaining female fomes in the same way, and updating the historical optimal positions and the historical optimal objective function values of all the remaining female fomes;

the current position of the male-fim 1 is calculated from the initial stage, and the current objective function value (J) is calculated from the current position of the male-fim 1 ₁ (m _θ ),J ₂ (m _θ ) Updating the historical position and the historical objective function value of the male dayfish 1; and calculating the current positions and the current objective function values of all the remaining male dayfds in the same way, and updating the historical optimal positions and the historical optimal objective function values of all the remaining male dayfds;

randomly selecting every 2 dayf individuals from the Pareto set for mating calculation, and randomly selecting half of the dayf individuals from the Pareto set for mutation calculation; combining the mated and mutated fava populations, and respectively incorporating 1/2 of the fava population into the current male and female fava populations;

respectively carrying out Pareto non-dominant ranking on the male and female populations, and taking q male and p female individuals ranked ahead as a new Pareto set; taking objective function values of all populations in the new Pareto set to obtain a Pareto set under the 1 st iteration;

and iterating for 99 times in the same way, and obtaining the Pareto set under the 100 th iteration as the final Pareto set of the optimization.

The current position of the female-type 1 in the initialization stage is calculated, and the current objective function value (J) is calculated according to the current position of the female-type 1 ₁ (m _θ ),J ₂ (m _θ ) A step of updating the current position and the current objective function value of the female-fipresent 1, comprising:

calculation is performed starting from the female f 1 at the initialization stage, at which time the male f 1 dominates the female f 1, and the European distance between the male f 1 and the female f 1 individuals is calculated.

Wherein r is _mf Representation ofAnd->European distance between->Representing the initial position in the j-th dimension of male-fimbriae 1,/and->An initial position representing a j-th dimension of female maydie 1;

thereafter, the current speed of female dayfish 1 is updated, and the speed calculation result of each dimension of female dayfish 1 is:

wherein,representing the current speed of the j-th dimension of the i-th female; g represents a gravity coefficient; a, a ₂ Representing a positive attraction constant; beta represents the modified visibility coefficient; />Represents the position of the j-th dimension of the ith male at time t,representing the position of the jth dimension of the ith female at time t; f () represents an objective function; LF is a random flight coefficient; r is [ -1,1]Is a random number of (a);

then, update the current position of female dayfish 1:

wherein,indicating the position of the ith female at the current time; />Representing the current speed of the ith female; />Indicating the position of the ith female at time t;

calculating a current objective function value (J) based on the current position of the female dayfish 1 ₁ (m _θ ),J ₂ (m _θ ) Update femaleThe current position and current objective function value of the dayfish 1;

the same calculation is carried out on the rest 99 female dayflets according to the method, the current positions and the current objective function values of all the rest female dayflets are calculated, and the historical optimal positions and the historical optimal objective function values of all the rest female dayflets are updated;

calculation starts from male f 1 at the initialization stage, where male f 1 governs the speed calculation results for each dimension of male f 1 as:

wherein,representing the current speed of the j-th dimension of the i-th male; />Representing the position of the j-th dimension of the ith male at time t; a, a ₁ 、a ₂ Is a positive attraction constant; pbest (p best) _ij Representing the current optimal position of the j-th dimension of the i-th male; beta represents the modified visibility coefficient; r is (r) _p Represents x _i And the current optimum position pbest _i Euclidean distance between them; r is (r) _g Represents x _i And the Euclidean distance between the global optimal position gbest;

then, the current position of stamen 1 is updated:

wherein,representing the position of the ith male at the current time; />Representing the current speed of the ith male; />Representing the position of the ith male at time t;

calculating a current objective function value (J) based on the current position of the male f-1 ₁ (m _θ ),J ₂ (m _θ ) Updating the historical optimal position and the historical optimal objective function value of the male dayfish 1.

Said step 5 comprises the steps of:

and assigning the position of the f in the Pareto set obtained when the incidence angle is theta to the particles, and using the position of the f as a posterior probability model for inversion when the incidence angle is theta and using the position of the f as an priori probability model for observation seismic record inversion when the incidence angle is theta+1 deg.

The step 6 specifically comprises the following steps:

step 3 to step 5 are executed 42 times for incidence angles theta=4 DEG to theta=45 DEG, a global final Pareto set is obtained, 100 f-positions on the global final Pareto set are assigned to 100 particles, expectations of all the particles are calculated, average values are respectively calculated for all the particles according to different dimensions, and final expected vectors are obtainedAs a result of the inversion, N represents the nth dimension of the corresponding longitudinal wave velocity, transverse wave velocity, and density.

Compared with the prior art, the invention has the beneficial effects that:

the invention adopts a multi-objective-type Algorithm to improve a particle filtering algorithm based on a Bayesian framework, adopts a multi-objective-type Algorithm to simultaneously optimize two objective functions, and can simultaneously minimize absolute errors and pearson correlation coefficients of an observed seismic record and a synthesized seismic record compared with a mode that a single objective optimization algorithm can only optimize one objective function. The special Pareto set in the multi-objective optimization algorithm is essentially one resampling of all the dayfies, and then the last sampling result is assigned to the particles, so that the resampling process in the standard particle filtering algorithm can be completely replaced, and the weight of each particle can be balanced and kept in a high-order likelihood domain by the Pareto non-dominant ordering. The multi-objective Algorithm and the particle filter algorithm are combined, so that the gain of 2 algorithms is improved in the pre-stack AVA inversion, a more accurate combined inversion method is obtained, the inversion accuracy is improved, the accurate inversion of three parameters of elasticity is realized based on the particle filter algorithm improved by the multi-objective Algorithm, and a powerful support is provided for accurate identification of a high-quality reservoir.

Drawings

In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings that are needed in the embodiments will be briefly described below, it being understood that the following drawings only illustrate some embodiments of the present invention and therefore should not be considered as limiting the scope, and other related drawings may be obtained according to these drawings without inventive effort for a person skilled in the art.

FIG. 1 is a flow chart of the pre-stack inversion method of the present invention;

FIG. 2 is the original pre-stack AVA longitudinal wave angle gather seismic data of embodiment 1 of the present invention;

FIG. 3 is a Pareto set at iteration 1 for an angle of incidence of 3℃in example 1 of the present invention;

FIG. 4 is the final Pareto set at iteration 100 for an angle of incidence of 3℃in example 1 of the present invention;

FIG. 5 is a final Pareto set for each incident angle in example 1 of the present invention when the incident angle is 3-45;

FIG. 6 is a graph of elastic three-parameter inversion obtained using a conventional local nonlinear method;

FIG. 7 shows the result of elastic three-parameter inversion obtained by the improved particle filter algorithm of example 1 of the present invention using the multi-objective F-algorithm.

Detailed Description

The following description of the embodiments of the present invention will be made clearly and completely with reference to the accompanying drawings, in which it is apparent that the embodiments described are only some embodiments of the present invention, but not all embodiments. The components of the embodiments of the present invention generally described and illustrated in the figures herein may be arranged and designed in a wide variety of different configurations. Thus, the following detailed description of the embodiments of the invention, as presented in the figures, is not intended to limit the scope of the invention, as claimed, but is merely representative of selected embodiments of the invention. All other embodiments, which can be made by a person skilled in the art without making any inventive effort, are intended to be within the scope of the present invention.

It should be noted that: like reference numerals and letters denote like items in the following figures, and thus once an item is defined in one figure, no further definition or explanation thereof is necessary in the following figures. Also, in the description of the present invention, the terms "first," "second," and the like are used merely to distinguish one from another, and are not to be construed as indicating or implying a relative importance or implying any actual such relationship or order between such entities or operations. In addition, the terms "connected," "coupled," and the like may be used to denote a direct connection between elements, or an indirect connection via other elements.

Example 1:

the invention adopts a multi-objective and dayf algorithm to improve a particle filtering algorithm, and realizes inversion of three parameters of pre-stack elasticity. The invention uses the special Pareto non-dominant sorting mode of multi-objective algorithm to obtain Pareto set, and keeps the field of the multi-objective algorithm in high-order likelihood domain, and optimally assigns the field to the particles, thereby replacing the particle resampling process in the traditional particle filtering algorithm; in addition, the traditional intelligent optimized particle swarm algorithm is used for directly solving the inversion equation, so that the large particle swarm number is required, the convergence of the algorithm is very high, and the method is difficult to achieve at present.

The practical pre-stack seismic data angle gather range applied in this embodiment is an integer angle of 3 ° to 45 °, the sampling interval is 0.002s, the time range is 2538 to 2788ms, and 125 data points are sampled in total.

Because the combination inversion algorithm involves more calculation and the range of each dimension parameter is larger, the general intelligent optimized particle swarm algorithm, genetic algorithm, simulated annealing algorithm and the like are used, so that the optimal value is difficult to obtain and the local optimization is easy to fall into. The particle filter algorithm is improved by using the multi-objective dayfish algorithm, the advantages of the multi-objective dayfish algorithm and the particle filter algorithm are effectively combined, and the method possibly applied is as follows.

The invention is realized by the following technical scheme, as shown in fig. 1, a pre-stack inversion method for improving a particle filtering algorithm based on a multi-objective dayf algorithm comprises the following steps:

step 1, collecting logging data and establishing an initial model.

Longitudinal wave velocity log data, transverse wave velocity log data, and density log data from the AM block W1 well are collected. And filtering out high-frequency signals in the logging data by a linear smoothing method, and establishing an initial model. Referring to Table 1, log data is shown, in part, for use in establishing an initial model.

Table 1 logging data to build initial model

And 2, assigning the initial model to the particles as a global prior probability model according to preset particle filtering algorithm parameters.

The parameters of the particle filtering algorithm include: particle count, process noise variance, observed noise variance, etc. The initial values of the parameters of the particle filter algorithm are preset as follows: particle number npf=100, observed noise variance q=2, process noise variance r=2.011, and the initial model is assigned to the particles in the particle filtering algorithm as a prior probability model.

Step 3, starting from the seismic trace at the first incident angle, performing a particle filtering algorithm prediction process;

assuming that the first incident angle is 3 °, a particle filter algorithm prediction step is performed starting from the seismic data at the 3 rd incident angle, i.e. the 100 particles are added with observation noise, and table 2 is a part of the particles added with observation noise.

Table 2 particles (part) with observation noise added thereto

Dimension(s)	Particle 1	Particle 2	Particle 3	…	Particle 100
						1	3726.24	3727.29	3729.85	…	3728.26
2	3840.86	3840.13	3846.28	…	3844.55
						3	4028.06	4030.05	4029.33	…	4028.07
4	3830.69	3824.85	3828.36	…	3827.23
						5	3610.24	3608.29	3606.21	…	3608.50
6	3590.00	3598.67	3595.86	…	3595.91
						7	3615.46	3613.52	3616.64	…	3615.07
8	3469.85	3472.98	3473.51	…	3468.23
						9	3459.64	3456.49	3455.25	…	3458.20
10	3496.52	3496.64	3498.95	…	3496.57

According to the basic theory of a standard particle filter algorithm, the nonlinear dynamic equation at the angle of theta DEG is as follows:

wherein m is _θ-1 Representing a posterior probability model after the angle theta-1 degrees is calculated; m is m _θ Representing the model m when calculating the angle θ° _θ-1 A priori probability model after adding process noise; f () represents a state transition equation; w (w) _θ Representing a process noise variance of the system; z _θ Representing the observed value of the system at an angle of θ ° (i.e., observing seismic records); h () represents the observation equation (i.e., the operation from the prior probability model to the synthetic seismic record); v _θ Representation ofThe system observes noise when calculating the θ° angle.

Note that, when θ=3°, the model parameter m is ₃ The different phases according to the prior and posterior can be defined as respectivelyAnd->Wherein->The model is an initial model (low-frequency model) and also serves as a global prior probability model; prior probability model->After calculation of the observation equation h (), a posterior probability model +.>Whereas a posterior probability model for θ=3° is +.>Calculated by state transition equation f (), then used as a prior probability model for θ=4°>The calculation is continued and reciprocated. When the angle is θ, the time θ is defined.

The kernel of the particle filtering algorithm is based on cyclic operation and probability accumulation of different angles theta (1 DEG increase each time) for a plurality of times, and the particle filtering algorithm passes through a global prior probability modelI.e. solve the final posterior probability model +. >The global implementation is shown at 43 angles43 prediction processes and an update process.

(1) The prediction process can be expressed as model parameters at a given angle θ -1Predicting model parameters at a given angle theta>Is a probability of (2). The predictive equation can be expressed in terms of an integral as:

wherein,a priori probability from 3 DEG to theta-1 DEG; d, d _obs(3:θ-1) Representing a set of observation seismic records from 3 ° to θ -1 ° time instants; />Is shown in event->Event under conditionsProbability of occurrence; (m) _θ-1 |d _obs(3:θ-1) ) Representing model parameters m _θ-1 In the seismic observation record set d _obs(3:θ-1) Distribution below.

(2) The updating process is the inversion of the invention, namely from model parametersUpdating to model parametersIs also a pre-stack inversion process in the field of seismic exploration,this intermediate uses the observed seismic record d at the current angle θ° _obs(θ) Updating the prior probability +.>To distinguish between the prediction process and the update process, the update process is denoted +.>Obtaining posterior probabilityThe update equation is:

/>

wherein,a posterior probability at a time from 3 DEG to theta DEG; p (d) _obs(θ) |d _obs(3:θ) ) Is a normalization constant; />Representing likelihood functions, which are also the specific calculation process for the inversion of the present invention;representing a priori probabilities at time θ °.

According to the full probability formula, the normalization constant is expressed in integral form as:

and solving an objective function established by a subsequent Bayesian theory, and calculating the maximum posterior probability to obtain posterior probability distribution of model parameters, namely an inversion result in the invention.

The updating process based on the multi-objective and the improved particle filtering algorithm is completed in the iteration of the multi-objective and the improved particle filtering algorithm, and the subsequent particle weight updating and the particle resampling steps in the improved particle filtering algorithm are all improved in the manner of a Pareto (Pareto) set solved by the multi-objective and the improved particle filtering algorithm.

And 4, iterating by adopting a multi-objective dayfish algorithm, finishing the updating process of the particle filtering algorithm during the iteration, and obtaining a final Pareto set under the incident angle after iterating for a plurality of times.

Before performing the multi-objective, dayfish algorithm iterations, multiple key parameters of the multi-objective, dayfish algorithm need to be initialized. Setting the number of male and female dayfish q=100, the number of mating father and dayfish p=100, the number of filial mating nc=100, the number of filial mutation nm=50, the gravity coefficient g=0.8, and the positive attraction constant a ₁ Positive attraction constant a =1 ₂ =1.5, mutation coefficient mu=0.02, modified visibility coefficient beta=2, random flight coefficient fl=0.77, wedding Dance coefficient dance=0.77, wedding Dance coefficient decrease rate dd=0.99, random flight coefficient decrease rate fd=0.99, pareto aggregate capacity np=100, maximum number of iterations t _max ＝100。

Determining the number of the dayfish population, and initializing the positions of the male and female dayfles.

M＝{x _i |x _i ＝(x _i1 ,x _i2 ,...,x _id ),i＝1,2,...,q}

F＝{y _i |y _i ＝(y _i1 ,y _i2 ,...,y _id ),i＝1,2,...,p}

Wherein M represents a population of male dayfish; f represents a female population; x is x _i Representing the position of the ith male; y is _i Representing the position of the ith female; q represents the number of male dayfish (q=100); p represents the female dayfish number (p=100); d represents the total number of dimensions (d=375); x is x _id Representing the position of the ith male dimension; y is _id Representing the position of the ith female and the d-th dimension.

In the initialization stage of the multi-objective Algorithm, the average value of the particle values after the prediction process is assigned to male and female Algorithms, and the initial position is taken as the initial position, and the initial speed is set to be 0; and the objective function value for each of the fakes is calculated at the initialization node.

From bayesian theory, the problem of inverting the elastic coefficient of the subsurface medium can be solved into a posterior probability under the condition that pre-stack seismic data are known, and the posterior probability is only calculated for the time theta:

Wherein,the posterior probability at time θ; d, d _obs(θ) The observation seismic record corresponding to the time theta is obtained; />The prior probability at the time theta is given; />Likelihood functions for observing seismic records; p (d) _obs(θ) ) To observe seismic record d _obs(θ) Can be expressed as 1 constant. Assuming that the likelihood function satisfies the gaussian distribution, there are:

wherein C is _n A variance matrix representing the actual noise,is C _n In the present algorithm, an identity matrix I may be multiplied by the observed noise v at time θ _θ To represent；/>Representing model parameters +.>Is used for the co-variance matrix of (a),is->An inverse matrix of (a); />Representing model parameters +.>The mean value vector of (different model parameters need to be calculated respectively); t represents the transpose of the matrix; />Representing the synthetic seismic record corresponding to the moment theta, namely z in a nonlinear dynamic equation established by a particle filtering algorithm _θ 。

Obtaining the maximum posterior probability distribution problem of model parameters, in a single-objective algorithm, the objective function value J can be minimized ₁ (m _θ )：

In the multi-objective algorithm, another objective function J can be additionally minimized ₂ (m)：

Wherein corrcoef (& gt) represents the pearson correlation coefficient for obtaining two sequences; due to observation of seismic recordsThe recorded values are larger, and the original seismic record is normalized to [ -1000,1000 before calculation ]Within the range; because the global optimization algorithm is difficult to directly solve the objective function, taylor first order approximation is generally adopted, so that the calculated amount is reduced, the smoothing effect is achieved at the same time, and the visual recognition of the inversion result is improved. At this time, update Then the update procedure in the nonlinear dynamic equation established by the standard particle filter algorithm is equivalent:

/>

where Δm represents the model delta, i.e., the current modelIs +.>Is a difference in (2);representing composite seismic records->For model parameters->Taking the first partial derivative can be expressed as:

wherein W represents a longitudinal wave matrix, K represents a Accord ratio matrix of longitudinal wave coefficients, and the Accord ratio matrix of reflection coefficients for an incident angle θ° is:

wherein R is _N Representing the reflection coefficient of the N-th dimension, V _PN 、V _SN And ρ _N Representing the shear wave velocity, the longitudinal wave velocity and the density corresponding to the nth dimension.

The reflectance terms in the Accord matrix K can be obtained by the exact Zoeppritz equation, which can be expressed as:

AR＝B

wherein,

/>

wherein V is _P1 、V _P2 Representing longitudinal wave velocity at both ends of the medium, V _S1 、V _S2 Representing transverse wave velocity, ρ, at both ends of the medium ₁ 、ρ ₂ Representing the density of both ends of the medium; r is R _PP 、R _PS 、T _PP 、T _PS The reflection coefficient of the longitudinal wave, the transmission coefficient of the longitudinal wave, the reflection coefficient of the transverse wave and the transmission coefficient of the transverse wave are respectively expressed.

The partial derivative of the reflection coefficient with respect to the model parameters can be expressed as:

wherein ψ= [ ψ ] ₁ ,ψ ₂ ,ψ ₃ ,ψ ₄ ]Is the proportion of elastic parameters, and and->The specific parameters of (a) are as follows:

/>

the 6 elastic parameters between the single reflection coefficient and the reflection interface can be obtained by the chain rule:

updating the current position of the 1 st male-f to be the initial position, and performing an updating process on the 1 st male-f to obtain a current objective function value of (J) ₁ (m _θ ),J ₂ (m _θ ) A) is provided; and the current objective function values of all the remaining male and female dayfis determined in the same manner.

Starting a first iteration:

calculating the current position of the male dayfish 1 from the current position of the male dayfish 1 in the initialization stageObjective function value (J) ₁ (m _θ ),J ₂ (m _θ ) Updating the current position and the current objective function value of stamen's dayfish 1; and calculating the current positions and the current objective function values of all the remaining male dayfds in the same way, and updating the historical optimal positions and the historical optimal objective function values of all the remaining male dayfds;

the current position of female dayfish 1 is calculated from the initial stage, and the current objective function value (J ₁ (m _θ ),J ₂ (m _θ ) Updating the historical position and the historical objective function value of female dayf 1; and calculating the current positions and the current objective function values of all the remaining female fomes in the same way, and updating the historical optimal positions and the historical optimal objective function values of all the remaining female fomes.

Explanation is made by taking a male dayfish 1 as an example:

the current position of male 1 is x ₁ Some of the position parameters are shown in table 3.

TABLE 3 position parameters of F1 (section)

Dimension(s)	VP/(m/s)	Dimension(s)	VS/(m/s)	Dimension(s)	DEN/(g/cm 3 )
						1	3727.22	126	2069.67	251	2420.25
2	3846.07	127	2164.23	252	2442.82
						3	4028.15	128	2315.69	253	2465.03
4	3831.72	129	2155.53	254	2439.81
						5	3611.35	130	1986.22	255	2402.21
6	3591.93	131	1960.97	256	2400.01
						7	3613.74	132	1981.03	257	2402.64
8	3471.89	133	1882.08	258	2377.97
						9	3458.60	134	1868.78	259	2379.66
10	3495.73	135	1899.52	260	2385.28
						...	...	...	...	...	...

Calculating a model increment delta m by the model increment delta m andmultiplying to obtain an increment of the synthetic seismic record; at this time, composite seismic record +.>Adding the observed noise v of the current angle θ _θ =5.5000 e+15, the update procedure in the standard particle filter algorithm is completed. Current model parametersThe calculated value is 2.2192e+03, and the current objective function J is obtained by calculating the residual error and the correlation coefficient between the synthetic seismic record and the observed seismic record ₁ (m _θ ) Has a value of 2.0568e+07, an objective function J ₂ (m _θ ) Is 2.2200e+03, the objective function value of male f 1 in the initialization stage is (2.0568e+07, 2.2200e+03). And the objective function values of all the remaining male and female dayfis determined in the same manner.

The initialization of the multi-objective dayfish algorithm is completed above.

Thereafter, a main loop of the multi-objective dayfish algorithm is entered, iterated with the multi-objective dayfish algorithm, and the positions (model parameters) of the dayfish are updated. The same calculation is performed for the remaining 99 male-falfs as described above, updating the historical optimal positions and the historical optimal objective function values for male-falfids 2-100. Similarly, 100 female dayfs are subjected to the same calculation according to the method, and the historical optimal positions and the historical optimal objective function values of the female dayfs 1-100 are updated.

All male and all female individuals are combined, one-time Pareto non-dominant sorting is carried out, np individuals before sorting are taken as Pareto sets, and the dayf individuals with the rank of 1 in the Pareto sets are selected as the head of the male and female individuals.

100 male and 100 female individuals were combined, and a total of 200 individuals were Pareto non-dominant ordered, with the top 100 (np=100) as the Pareto set. The dayfish individuals in the Pareto collection ranked 1 are selected as the leading dayfish.

Start iteration 1 (t=1):

calculation is performed starting from the female f 1 at the initialization stage, at which time the male f 1 dominates the female f 1, and the European distance between the male f 1 and the female f 1 individuals is calculated.

Wherein r is _mf Representation ofAnd->European distance between->Representing the initial position in the j-th dimension of male-fimbriae 1,/and->Representing the initial position of the j-th dimension of female dayfish 1.

Thereafter, the current speed of female dayfish 1 is updated, and the speed calculation result of each dimension of female dayfish 1 is:

wherein,representing the current speed of the j-th dimension of the i-th female; g represents a gravity coefficient; a, a ₂ Representing a positive attraction constant; beta represents the correctedIs a visible coefficient of (2); />Represents the position of the j-th dimension of the ith male at time t, Representing the position of the jth dimension of the ith female at time t; f () represents an objective function; LF is a random flight coefficient; r is [ -1,1]Is a random number of the random number group.

Then, update the current position of female dayfish 1:

wherein,indicating the position of the ith female at the current time; />Representing the current speed of the ith female; />Indicating the position of the ith female at time t.

Calculating a current objective function value (J) based on the current position of the female dayfish 1 ₁ (m _θ ),J ₂ (m _θ ) The current position and current objective function value of female-fican 1 are updated.

The same calculation is performed on the remaining 99 female dayflets according to the above method, the current positions and current objective function values of all the remaining female dayflets are calculated, and the historical optimal positions and historical optimal objective function values of all the remaining female dayflets are updated.

Calculation starts from male f 1 at the initialization stage, where male f 1 governs the speed calculation results for each dimension of male f 1 as:

/>

wherein,representing the current speed of the j-th dimension of the i-th male; />Representing the position of the j-th dimension of the ith male at time t; a, a ₁ 、a ₂ Is a positive attraction constant; pbest (p best) _ij Representing the current optimal position of the j-th dimension of the i-th male; beta represents the modified visibility coefficient; r is (r) _p Represents x _i And the current optimum position pbest _i Euclidean distance between them; r is (r) _g Represents x _i And the euclidean distance between the global optimum position gbest.

Then, the current position of stamen 1 is updated:

wherein,representing the position of the ith male at the current time; />Representing the current speed of the ith male; />The position of the ith male at time t is indicated.

The current objective function value (J) is calculated from the current position of male f 1 ₁ (m _θ ),J ₂ (m _θ ) At this time, the historical optimal objective function value of male-fiveleaf 1 governs the current objective function value of male-fiveleaf 1, and randomly selects a range of [0,1]Currently, the number of (2)The number of choices is 0.436 (less than 0.5), the history optimal position of male dayf 1 is updated, and the history optimal objective function value of male dayf 1 is updated.

The same calculation is performed on the remaining 99 male dayfs according to the above method, and if the current objective function value of the male dayf 1 dominates the historical optimal objective function value of the male dayf 1, the historical optimal position and the historical optimal objective function value of the male dayf 1 are updated; if the history optimal objective function value of male 1 governs the current objective function value of male 1, then there is a probability of 1/2 to update the history optimal position and history optimal objective function value of male 1.

For example, when the male-dayfish 77 does not satisfy the above conditions, the history optimum position and history optimum objective function values of the male-dayfish 2 to 76 and the male-dayfish 78 to 100 are updated.

Then randomly selecting every 2 f individuals from the Pareto set for mating calculation, and randomly selecting half f individuals from the Pareto set for mutation calculation; the mated and mutated fava populations were combined, 1/2 each was incorporated into the current male and female populations, respectively.

Next, pareto non-dominant ordering is performed on the male and female populations, respectively, and the q male and p female individuals ranked first are taken as a new Pareto set. And obtaining the Pareto set under the 1 st iteration by taking the objective function values of all the populations in the new Pareto set.

For example, 2 individuals from the Pareto set were randomly selected and subjected to the 1 st mating calculation, and the partial calculation results are shown in table 4.

TABLE 4 Efalciform population 1 st mating calculation (part)

The remaining 98 individuals from the Pareto collection were subjected to 49 more mating calculations.

After the completion of the mating calculation for 50 times, mutation calculation is performed on 50 selected 1 st individual of the Pareto set, and the current objective function value of the mutated individual is updated, and the mutation dimension calculation result is shown in table 5.

TABLE 5 calculation of mutation dimension of the 1 st individual

The mated and mutated falciform populations were combined, and a total of 58 falciform individuals were pooled together in 1 set, and then 1/2 of the falciform populations were incorporated into the initializing stage male and female falciform populations, respectively. Thereafter, the male and female dayfish populations were Pareto non-dominant ordered, respectively, with each male dayfish individual with the front 50 of the order and the female dayfish individual with the front 50 of the order combined as a new Pareto collection. Pareto non-dominant ranking is performed on the populations in the new Pareto collection. Finally, the objective function values of the top 100 populations are taken as Pareto sets at iteration 1, see fig. 3.

The parameters of the multi-objective dayfish algorithm, i.e., the wedding Dance coefficients, are modified by dance=0.77×0.99, and the random flight coefficients are modified by fl=0.77×0.99.

And iterating for 99 times in the same way, and obtaining the Pareto set under the 100 th iteration as the final Pareto set of the optimization, see fig. 4.

And 5, assigning the position of the head in the final Pareto set obtained by the incidence angle to the particle, and using the position of the head as a posterior probability model for inverting the incidence angle and using the position of the head as an priori probability model for inverting the seismic record of the observation of the next incidence angle.

And assigning the f-positions in the Pareto set acquired when the incidence angle is theta=3 degrees to the particles, wherein the f-positions are used as a posterior probability model for inversion when the incidence angle is theta=3 degrees, and the f-positions are used as an prior probability model for observation seismic record inversion when the incidence angle is theta=4 degrees.

And 6, repeating the steps 3 to 5 until the final global Pareto set is obtained after the observation seismic records of all the incidence angles are calculated, assigning the positions of the head on the final global Pareto set as particles, and calculating expected values of all the head to obtain an inversion result.

Step 3 to step 5 are executed 42 times for incidence angles theta=4 DEG to theta=45 DEG, a global final Pareto set is obtained, 100 f-positions on the global final Pareto set are assigned to 100 particles, expectations of all the particles are calculated, average values are respectively calculated for all the particles according to different dimensions, and final expected vectors are obtainedAs an inversion result, N represents the nth dimension of the corresponding longitudinal wave velocity, transverse wave velocity and density, and a partial result is shown in table 6, and the final Pareto set at each incident angle is shown in fig. 5.

Watch 6100 expectations of particles (section)

Dimension(s)	VP/(m/s)	Dimension(s)	VS/(m/s)	Dimension(s)	DEN/(g/cm 3 )
						1	3653.002	126	2028.87	251	2406.261
2	3849.892	127	2149.502	252	2442.033
						3	4044.108	128	2318.26	253	2471.697
4	3828.972	129	2105.174	254	2443.106
						5	3625.422	130	1980.983	255	2405.535
6	3627.05	131	1905.8	256	2417.457
						7	3650.477	132	1926.128	257	2420.68
8	3470.682	133	1874.296	258	2379.322
						9	3476.499	134	1868.723	259	2380.72
10	3530.83	135	1894.978	260	2388.637

In summary, the present invention performs elastic three-parameter inversion on the well log data of the investigation region, please refer to fig. 6. Compared with the conventional inversion method (damping least square method), please refer to fig. 7, the inversion algorithm result established by the invention is more fit with the real logging, and the accuracy is higher. The resampling process of the standard particle filter algorithm is replaced by the Pareto non-dominant sorting mode, so that the phenomenon of particle lean degradation is effectively avoided, and the position of a f-in the Pareto set is ensured to be in a high-order paraphrasing domain; meanwhile, compared with a single-target optimization algorithm, the multi-target optimization algorithm additionally constrains correlation coefficient indexes between the seismic observation records and the synthetic seismic records. Therefore, the accuracy of pre-stack AVA inversion is greatly improved by reasonably fusing the particle filtering algorithm and the multi-objective dayf algorithm. The particle filter algorithm based on the multi-objective dayf algorithm improvement established by the invention realizes accurate inversion of three parameters of elasticity and provides powerful support for accurate identification of high-quality reservoirs.

The foregoing is merely illustrative of the present invention, and the present invention is not limited thereto, and any person skilled in the art will readily recognize that variations or substitutions are within the scope of the present invention. Therefore, the protection scope of the present invention shall be subject to the protection scope of the claims.

Claims (7)

1. The prestack inversion method for improving the particle filter algorithm based on the multi-objective dayfish algorithm is characterized by comprising the following steps of: the method comprises the following steps:

step 1, collecting logging data and establishing an initial model;

step 2, according to preset particle filtering algorithm parameters, assigning an initial model as a global prior probability model to particles;

step 3, starting from the seismic trace at the first incident angle, performing a particle filtering algorithm prediction process;

the step 3 specifically comprises the following steps:

establishing a nonlinear dynamic equation:

wherein m is _θ-1 Representation ofCalculating a posterior probability model after the angle theta-1 DEG is finished; m is m _θ Representing the model m when calculating the angle θ° _θ-1 A priori probability model after adding process noise; b () represents a state transition equation; w (w) _θ Representing a process noise variance of the system; z _θ Representing an observed seismic record of the system at an angle of θ°; h () represents the observation equation; v _θ Representing observation noise of the system when calculating the theta angle;

model parameters at a given angle θ -1 momentPredicting model parameters at a given angle theta>Is predicted by the following equation:

wherein,model parameters at the time theta-1; />Model parameters at time theta; />A priori probability from 3 DEG to theta-1 DEG; d, d _obs(3:θ-1) Representing a set of observation seismic records from 3 ° to θ -1 ° time instants;is shown in event->Event under Condition->Probability of occurrence; (m) _θ-1 |d _obs(3:θ-1) ) Representing model parameters m _θ-1 In the seismic observation record set d _obs(3:θ-1) Distribution of the lower part;

step 4, iteration is carried out by adopting a multi-objective dayf algorithm, the updating process of the particle filtering algorithm is completed in the period, and a final Pareto set under the incident angle is obtained after the iteration is carried out for a plurality of times;

the step 4 specifically comprises the following steps:

determining the number of the dayfish population, and initializing the positions of the male and female dayfles.

M＝{x _i |x _i ＝(x _i1 ,x _i2 ,...,x _id ),i＝1,2,...,q}

F＝{y _i |y _i ＝(y _i1 ,y _i2 ,...,y _id ),i＝1,2,...,p}

Wherein M represents a population of male dayfish; f represents a female population; x is x _i Representing the position of the ith male; y is _i Representing the position of the ith female; q represents the number of male dayfish; p represents the number of female dayflower; d represents the total number of dimensions; x is x _id Representing the position of the ith male dimension; y is _id Representing the position of the ith female f iotaeld d dimension;

model parameters are known only at the θ timeUpdating model parameters +.>The process of (1) is as follows:

wherein,the posterior probability at time θ; d, d _obs(θ) The observation seismic record corresponding to the time theta is obtained;the prior probability at the time theta is given; />Likelihood functions for observing seismic records; p (d) _obs(θ) ) To observe seismic record d _obs(θ) Can be expressed as 1 constant;

assuming that the likelihood function satisfies the gaussian distribution, there are:

wherein C is _n A variance matrix representing the actual noise,is C _n Is multiplied by an identity matrix I by the observed noise v at time θ _θ To represent; />Representing model parameters +.>Covariance matrix of>Is->An inverse matrix of (a);representing model parameters +.>Is a mean vector of (a); t represents the transpose of the matrix; />Representing the synthetic seismic record corresponding to the moment theta, namely z in a nonlinear dynamic equation established by a particle filtering algorithm _θ Expansion by Taylor's formula>And omit higher-order terms above second order, update +.>

Minimizing the objective function J ₁ (m _θ )：

Minimizing the objective function J ₂ (m _θ )：

Wherein Δm represents the model increment, i.e. the current model parameterAnd initial model parameters- >Is a difference in (2);representing composite seismic records->For model parameters->Taking the first partial derivative, expressed as:

wherein W represents a longitudinal wave matrix, K represents a Accord ratio matrix of longitudinal wave coefficients, and the Accord ratio matrix of reflection coefficients for an incident angle θ° is:

wherein R is _N Representing the reflection coefficient of the N-th dimension, V _PN 、V _SN And ρ _N Respectively representing the transverse wave speed, the longitudinal wave speed and the density corresponding to the Nth dimension;

the reflectance terms in the Accord matrix K can be obtained by the exact Zoeppritz equation, which can be expressed as:

AR＝B

wherein,

wherein V is _P1 、V _P2 Representing longitudinal wave velocity at both ends of the medium, V _S1 、V _S2 Representing transverse wave velocity, ρ, at both ends of the medium ₁ 、ρ ₂ Representing the density of both ends of the medium; r is R _PP 、R _PS 、T _PP 、T _PS Respectively representing the longitudinal wave reflection coefficient, the longitudinal wave transmission coefficient, the transverse wave reflection coefficient and the transverse wave transmission coefficient;

the partial derivative of the reflection coefficient with respect to the model parameters can be expressed as:

wherein ψ= [ ψ ] ₁ ,ψ ₂ ,ψ ₃ ,ψ ₄ ]Is the proportion of elastic parameters, and and->The specific parameters of (a) are as follows:

the 6 elastic parameters between the single reflection coefficient and the reflection interface can be obtained by the chain rule:

updating the current position of the 1 st male-f to be the initial position, and performing an updating process on the 1 st male-f to obtain a current objective function value of (J) ₁ (m _θ ),J ₂ (m _θ ) A) is provided; and determining the current objective function values of all the remaining male and female dayfles in the same manner;

combining all male and all female individuals, performing one-time Pareto non-dominant sorting, taking Np individuals before sorting as a Pareto set, and selecting the dayf individuals with the rank of 1 in the Pareto set as the head of the male and female individuals;

starting a first iteration:

the current position of female dayfish 1 is calculated from the initial stage, and the current objective function value (J ₁ (m _θ ),J ₂ (m _θ ) Updating the current position and the current objective function value of the female-figet 1; calculating the current positions and the current objective function values of all the remaining female fomes in the same way, and updating the historical optimal positions and the historical optimal objective function values of all the remaining female fomes;

the current position of the male-fim 1 is calculated from the initial stage, and the current objective function value (J) is calculated from the current position of the male-fim 1 ₁ (m _θ ),J ₂ (m _θ ) Updating the history position and the history order of the male dayfish 1Marking a function value; and calculating the current positions and the current objective function values of all the remaining male dayfds in the same way, and updating the historical optimal positions and the historical optimal objective function values of all the remaining male dayfds;

Randomly selecting every 2 dayf individuals from the Pareto set for mating calculation, and randomly selecting half of the dayf individuals from the Pareto set for mutation calculation; combining the mated and mutated fava populations, and respectively incorporating 1/2 of the fava population into the current male and female fava populations;

respectively carrying out Pareto non-dominant ranking on the male and female populations, and taking q male and p female individuals ranked ahead as a new Pareto set; taking objective function values of all populations in the new Pareto set to obtain a Pareto set under the 1 st iteration;

iterating for 99 times in the same way, and obtaining a Pareto set under the 100 th iteration as a final Pareto set of the optimization;

step 5, assigning the f-positions in the final Pareto set obtained by the incidence angle to the particles, and using the f-positions as a posterior probability model for inverting the incidence angle and a priori probability model for inverting the seismic record of the observation of the next incidence angle;

and 6, repeating the steps 3 to 5 until the final global Pareto set is obtained after the observation seismic records of all the incidence angles are calculated, assigning the positions of the head on the final global Pareto set as particles, and calculating expected values of all the head to obtain an inversion result.

2. The method of pre-stack inversion based on the multi-objective dayf algorithm to improve the particle filtering algorithm of claim 1, wherein: the step 1 specifically comprises the following steps: the collected logging data comprises longitudinal wave velocity logging data, transverse wave velocity logging data and density logging data; and filtering out high-frequency signals in the logging data by a linear smoothing method, and establishing an initial model.

3. The method of pre-stack inversion based on the multi-objective dayf algorithm to improve the particle filtering algorithm of claim 1, wherein: the step 2 specifically comprises the following steps: the particle filter algorithm parameters include particle number, observed noise variance, and process noise variance.

4. The method of pre-stack inversion based on the multi-objective dayf algorithm to improve the particle filtering algorithm of claim 1, wherein: in step 4, before iterating with the multi-objective dayf algorithm, the method further comprises the steps of:

initializing parameters of a multi-objective dayfish algorithm, setting a male dayfish number q=100, a female dayfish number p=100, a male dayfish mating number nc=100, a offspring dayfish mutation number nm=50, a gravity coefficient g=0.8, and a positive attraction constant a ₁ Positive attraction constant a =1 ₂ =1.5, mutation coefficient mu=0.02, modified visibility coefficient beta=2, random flight coefficient fl=0.77, wedding Dance coefficient dance=0.77, wedding Dance coefficient decrease rate dd=0.99, random flight coefficient decrease rate fd=0.99, pareto aggregate capacity np=100, maximum number of iterations t _max ＝100。

5. The method of pre-stack inversion based on the multi-objective dayf algorithm to improve the particle filtering algorithm of claim 1, wherein: the current position of the female-type 1 in the initialization stage is calculated, and the current objective function value (J) is calculated according to the current position of the female-type 1 ₁ (m _θ ),J ₂ (m _θ ) A step of updating the current position and the current objective function value of the female-fipresent 1, comprising:

calculation is performed starting from the female f 1 at the initialization stage, at which time the male f 1 dominates the female f 1, and the European distance between the male f 1 and the female f 1 individuals is calculated.

Wherein r is _mf Representation ofAnd->European distance between->Representing the initial position in the j-th dimension of male-fimbriae 1,/and->An initial position representing a j-th dimension of female maydie 1;

thereafter, the current speed of female dayfish 1 is updated, and the speed calculation result of each dimension of female dayfish 1 is:

wherein,representing the current speed of the j-th dimension of the i-th female; g represents a gravity coefficient; a, a ₂ Representing a positive attraction constant; beta represents the modified visibility coefficient; />Represents the position of the jth dimension of the ith male at time t,/f>Representing the position of the jth dimension of the ith female at time t; f () represents an objective function; f (F)LIs a random flight coefficient; r is [ -1,1]Is a random number of (a);

then, update the current position of female dayfish 1:

wherein,indicating the position of the ith female at the current time; />Representing the current speed of the ith female; />Indicating the position of the ith female at time t;

calculating a current objective function value (J) based on the current position of the female dayfish 1 ₁ (m _θ ),J ₂ (m _θ ) Updating the current position and the current objective function value of the female-figet 1;

the same calculation is carried out on the rest 99 female dayflets according to the method, the current positions and the current objective function values of all the rest female dayflets are calculated, and the historical optimal positions and the historical optimal objective function values of all the rest female dayflets are updated;

calculation starts from male f 1 at the initialization stage, where male f 1 governs the speed calculation results for each dimension of male f 1 as:

wherein,representing the current speed of the j-th dimension of the i-th male; />Representing the position of the j-th dimension of the ith male at time t; a, a ₁ 、a ₂ Is a positive attraction constant; pbest (p best) _ij Representing the current optimal position of the j-th dimension of the i-th male; beta represents the modified visibility coefficient; r is (r) _p Represents x _i And the current optimum position pbest _i Euclidean distance between them; r is (r) _g Represents x _i And the Euclidean distance between the global optimal position gbest;

then, the current position of stamen 1 is updated:

wherein,representing the position of the ith male at the current time; />Representing the current speed of the ith male; />Representing the position of the ith male at time t;

calculating a current objective function value (J) based on the current position of the male f-1 ₁ (m _θ ),J ₂ (m _θ ) Updating the historical optimal position and the historical optimal objective function value of the male dayfish 1.

6. The method of pre-stack inversion based on the multi-objective dayf algorithm to improve the particle filtering algorithm of claim 1, wherein: said step 5 comprises the steps of:

and assigning the position of the f in the Pareto set obtained when the incidence angle is theta to the particles, and using the position of the f as a posterior probability model for inversion when the incidence angle is theta and using the position of the f as an priori probability model for observation seismic record inversion when the incidence angle is theta+1 deg.

7. The method of pre-stack inversion based on the multi-objective dayf algorithm to improve the particle filtering algorithm of claim 1, wherein: the step 6 specifically comprises the following steps:

Step 3 to step 5 are executed 42 times for incidence angles theta=4 DEG to theta=45 DEG, a global final Pareto set is obtained, 100 f-positions on the global final Pareto set are assigned to 100 particles, expectations of all the particles are calculated, average values are respectively calculated for all the particles according to different dimensions, and final expected vectors are obtainedAs a result of the inversion, N represents the nth dimension of the corresponding longitudinal wave velocity, transverse wave velocity, and density.