CN113419278B - Well-seismic joint multi-target simultaneous inversion method based on state space model and support vector regression - Google Patents

Abstract

The invention relates to a well-seismic joint multi-target simultaneous inversion method based on a state space model and support vector regression, which comprises the following steps of: based on measured data of a sample oil well, a state space model and a state transition transformation matrix are established by a discretization method, then a state residual error random variable probability density function is established, an observation function is established based on a support vector regression method, then an observation error random variable probability density function is established, and finally multi-target simultaneous inversion is carried out based on the state space model. The method greatly reduces the difficulty of establishing the state space model, saves manpower, material resources and time, and solves the problems of long period, high cost and low production efficiency of the conventional inversion method.

Description

Well-seismic joint multi-target simultaneous inversion method based on state space model and support vector regression

The technical field is as follows:

the invention relates to the field of seismic data inversion and interpretation, in particular to a well-seismic combined multi-target simultaneous inversion method.

Background art:

the current stochastic inversion method based on logging constraints is a mainstream seismic inversion method, and research and engineering personnel often compare and discuss deterministic inversion and stochastic inversion as a characteristic of the prominent stochastic inversion method. Deterministic inversion assumes that the wave impedance is spatially a definite value, and is usually solved using a minimization criterion based on a convolution model to obtain a smooth estimated wave impedance value. Because seismic data are band-limited, the maximum limitation of deterministic inversion is that the inversion result is lack of both low frequency and high frequency components, the low frequency component is usually compensated by prestack time/depth migration velocity spectrum data or logging low frequency data, and the high frequency component is mainly compensated by logging data. Constrained sparse impulse inversion and model-based inversion are the two most common deterministic inversion methods. Stochastic seismic inversion assumes that the wave impedance is a random variable in space and can be represented by a probability distribution, and any sampling of the probability distribution is realized by one inversion. The random seismic inversion fuses seismic data, well logging data and geostatistical information into posterior probability distribution of the underground wave impedance, the posterior probability distribution is sampled by a Markov chain Monte-Lolo method, and the property of the posterior probability distribution is researched by comprehensively analyzing a plurality of sampling results. It can be seen that the deterministic inversion method focuses on obtaining the model with the highest probability in this distribution, while the stochastic inversion studies the properties of the entire probability distribution and is complementary to the deterministic inversion. Compared with the deterministic inversion method, the stochastic inversion has the following advantages: firstly, the method comprehensively utilizes logging information and geostatistical information, and the inversion process does not adopt local smoothing processing in deterministic inversion, so that more detailed information can be extracted from seismic data; secondly, random inversion starts from a well point, and the coincidence degree of an inversion result and the well is very high; thirdly, quantitative evaluation can be made on the uncertainty of the inversion result by comprehensively analyzing a plurality of implementations; finally, the latest random inversion method is established on the basis of the Betus formula, and multi-scale information can be conveniently fused. A large number of use experiences show that the transverse change of a deterministic inversion result is natural, but the resolution in the longitudinal direction is low, the single sand body and well drilling result is not good, and the actual production requirement cannot be met; the longitudinal resolution of the random inversion based on the logging constraint is greatly improved, the boundary of the thin sandstone is clearer and basically accords with the sandstone explained by drilling, because the initial wave impedance model contains high-frequency logging information exceeding a seismic frequency band and is kept in an inversion result, the high-frequency information corresponds to the response of the thin layer, and because the thin layer is usually small in distribution scale and quick in transverse phase change, enough wells are needed to ensure the reliability of an interpolation result, the random inversion based on the logging constraint provides higher requirements for the quantity of logging information and the spatial distribution of the logging information, and the random inversion is more suitable for reservoir seismic description in the development and production stages.

The state space model is a dynamic time sequence model, which can be expressed as equation (1) and equation (2):

X^t＝f^t(X^t-1)+Q^t (1)

Y^t＝o^t(X^t)+R^t (2)

wherein X^t＝f^t(X^t-1)+Q^tIs a random process representing the process of migration of a state within the model, X, which is not directly known with time sequence^tIs a random variable at time t, representing the state at time t, X^t-1Is a random variable at time t-1, representing the state at time t-1, f^t(X^t-1) Is the state transition function at time t, which may be arbitrary, most commonly a Markov chain state transition matrix, Q^tIs a state noise random variable at time t; y is^t＝o^t(X^t)+R^tIs a state observation at time t, representing a measurement of a state that is not directly known within the model at time t, o^t(X^t) Is a state observation function at time t, R^tIs a random variable representing the measurement error. If the time is discrete, i.e. the value of t can only take non-negative integers, X^t、Q^t、Y^t、R^tAll Gaussian random variables, o^t(X^t) The model is a linear function, and the model is a discrete time linear Gaussian system which is equivalent to Gaussian process regression in the field of machine learning. For convenience of description, let θ ═ X^t、Q^t、Y^t，R^t，f^t，o^t) Are parameters of the model. Based on the state space model, any inversion job can be expressed as: given model parameter theta and batch observation data y ═ y { (y)^t-n+1，y^t-n+2，...，y^t-1，y^tUnder the condition of the model, solving the state sequence x of the model as { x ═ x^t-n+1，x^t-n+2，...，x^t-1，x^tAnd n is the length of the batch of observation data y. Based on Bayesian framework, random variable X can be established^tA posterior probability density function p (X)^t|y，θ)＝c*p(y|X^t，θ)*p(X^tI θ), where a posterior probability density function p (X)^tY, θ) is determined, a random sequence X is established as X^t-n+1，X^{t -n+2}，...，X^t-1，X^tThe joint posterior probability density function of:

p(X|Y，θ)＝p(X^t-n+1，X^t-n+2，...X^t-1，X^t|Y，θ)＝p(X^t-n+1|Y，θ)*p(X^t-n+2|X^t-n+1，Y，θ)*p(X^t-n+3|X_t-n+2，Y，θ)...p(X^t|X^t-1，Y，θ)。

the sequence of states x ═ x can be solved by two methods^t-n+1，x^t-n+2，...，x^t-1，x^tOne is maximum a posteriori, i.e.

Second is Bayesian inference, i.e. solving for the expectation of x, i.e.

x＝{x^t-n+1，x^t-n+2，...，x^t-1，x^t}＝E_p(X|Y，θ)(X)。

At present, state estimation methods built on state space models are mature, and some classical methods and their variants, such as hidden markov chain decoding method, kalman filtering, particle filtering, etc., have been developed. For some complex systems, how to build the state space model, i.e. how to give the model parameter θ ═ X^t、Q^t、Y^t、R^t、f^t、o^t) There are many difficulties, because the small to atomic, large to celestial body, whether it is natural science or human society science, it may be a core target of scientific research to establish and perfect a state space model, and it needs to pay a lot of manpower and material resources. Machine learning is a multi-disciplinary cross specialty, covers probability theory knowledge, statistical knowledge, approximate theoretical knowledge and complex algorithm knowledge, uses a computer as a tool and is dedicated to a real-time simulation human learning mode, and knowledge structure division is carried out on the existing content to effectively improve learning efficiency. By using the machine learning method, the difficulty of establishing the state space model can be greatly reduced, and manpower, material resources and time are saved. At present, most oil fields are in the middle and later stages of exploitation, a large amount of earthquake, well logging, statistical prediction and other data are accumulated, how to mine more information from ever-increasing mass data for decision making and insights is achieved, the period from seismic data acquisition, processing and explanation to well location planning is greatly shortened, the decision making precision is improved, the development cost is saved, and the generation efficiency is improved.

The invention content is as follows:

the invention can solve the problems by the following technical scheme: a well-to-seismic joint multi-target simultaneous inversion method based on a state space model and support vector regression comprises the following steps:

the first step is as follows: establishing a state space model by using a discretization method based on the measured data of the sample oil well, wherein the measured data of the sample oil well comprises the longitudinal wave speed, the transverse wave speed and the density of a measured sample curve;

the second step is that: establishing a state transition transformation matrix, wherein firstly, discretizing measured data of a sample oil well by utilizing discretization transformation;

the third step: establishing a state residual error random variable probability density function;

the fourth step: establishing an observation function based on a support vector regression method;

the fifth step: establishing a probability density function of an observation error random variable;

and a sixth step: and performing multi-target simultaneous inversion based on the state space model.

Further, in the first step, the discretization method is a K-Means discretization method, specifically, setting { x^tThe sample data { X ] is discretized in a cluster K-Means^tBefore, will { x }^tNormalizing continuous attribute values of all the dimensional data to a value range with the mean value of 0 and the variance of 1; the process of discretizing sample data by K-Means is recorded with S (), and X is set^t＝f^t(X^t-1)+Q^tWherein X is^t∈RⁿIs a model state, R is a real number, n is a state X^tFor setting the state space, it is necessary to set the dimensions of (A) for the continuous states X^tDiscretizing, setting S^tE.g. {1, 2, 3.. k }, k being an integer greater than 1, S^t＝S(X^t)，

Wherein S () is the K-Means discretization sample data process, S^-1() For the inverse discretization transform, it is assumed that there is sufficient sample data { x }^t}，X^t＝f^t(X^t-1)+Q^tRewritable to S^t-1＝S(X^t-1)，

{S^tIs regarded as a stationary Markov chain stochastic process, F^t() Is { S^tThe state transition transformation matrix is a time-independent transformation matrix, so F^t() The state space model can be written as F (), and after the state data X of the sample data is discretized by K-Means, the state space model can be adjusted as the following formulas (3), (4), (5) and (6):

S^t＝S(X^t) (3)

S^t＝F^t(S^t-1) (4)

X^t＝S^-1(S^t-1)+Q^t (5)

Y^t＝o^t(X^t)+R^t (6)

further, in the first step, the longitudinal wave velocity, the shear wave velocity and the density of the formation can be simultaneously inverted, wherein n is 3, and X is^t∈R³，

Wherein

As the velocity of the longitudinal wave,

as the velocity of the transverse wave,

is the density.

Further, in the second step, the method for establishing the state transition transformation matrix is specifically to assume that there is sufficient sample data { x }^tFirstly, discretizing the sample by discretizing transform S (), i.e. { S }^t}＝S({x^t}) of which s^tIs the dispersion of sample dataSequence of states, according to equation (3) for successive states X^tDiscretizing to obtain state space {1, 2, 3.. k }, S^tE {1, 2, 3.. k } is a state of the state space, k is an integer greater than 1, and F () is a k-th order square matrix of which each element F_ijRepresents the probability of a state transition from i to j in a Markov chain, in s^tTaking the statistic data sample of the Markov chain state transition probability matrix F as a statistic data sample, and counting each element F of the F_ij＝p(s^t＝j|s^t-1I) where P(s)^t＝j|s^t-1I) denotes the probability that the t-th state is j under the condition that the t-1 th state is i.

Further, in the third step, a state residual error random variable probability density function is established specifically according to X^t＝S^-1(F^t(S^t-1))+Q^tCan obtain the product

Wherein Q^tIs a state residual random variable, due to S^t-1Being discrete random variables, S^t＝F^t(S^t-1) Also referred to as a discrete random variable,

as state S^tCorresponding to a continuous random variable X^tExpected value of (2) X^tViewed as a k-dimensional Gaussian distribution, i.e.

Wherein epsilon²Is X^tThe covariance matrix of (2) is a k-th order square matrix because

So Q^tAnd k-dimensional Gaussian distribution, namely obtaining a state residual error random variable probability density function:

Q^t～N(0，ε²) (7)

where 0 is the mean value,. epsilon²Is a covariance matrix of X^tThe same covariance matrix is used, based on sufficient samplesData { x^tThe state transition transformation matrix F (), the discretization transformation S () and the discretization inverse transformation S established in the first step and the second step^-1() And calculating to obtain a data residual error sequence { q^tWherein q is^t＝x^t-S^-1(s^t)，s^t＝S(x^t) I.e. q^t＝x^t-s^-1(S(x^t) Namely { q) }^t}＝{x^t-S^-1(S(x^t) Q) according to { q }^tConveniently, to statistically infer epsilon²。

Further, in the fourth step, an observation function is established based on the support vector regression method, specifically, the observation function Y is established according to the state space model^t＝o^t(X^t)+R^tWherein o is^t() Representing a state observation function, order

To obtain

{x^tIs the state data sequence in the sample, { y^tIs the corresponding observed data sequence, { x^tAnd { y }^tComposition of State Observation Pair sequences { (x)^t，y^t) Using SVR from { (x)^t，y^t) Learn o^t() The learning process is as follows: to X^tPerforming dimension increasing treatment, namely, increasing the dimension of the original random variable which is obtained n-element for 1 time into n-element random variable for 3 times; if order X^tIs a 2-dimensional random variable, i.e.

Set the upgraded random variable to

Namely, it is

Namely, it is

Is a 9-dimensional vector; is provided with Y^tIs a 4-dimensional random variable, obviously

R^tAnd is also a 4-dimensional vector,

to get o therein^t() Viewed as a linear transformation, i.e. o^t() A matrix of 4 rows and 9 columns, denoted W, i.e.

Called W the observation matrix, so sample data { (x)^t，y^t) Can be expressed as

That is, the objective function of SVR can be established according to the sample data

Wherein

I.e. the objective function is:

wherein | W | is the sum of the modes of the row vectors of the matrix W, and c is a regularization factor, and W is obtained by an SVR learning method.

Further, in the fifth step, the probability density function of the random variable of the observation error is established specifically according to

Can obtain the product

Wherein

R^tIn order to observe the residual random variable,

as a state

Corresponding observed random variable Y^tExpected value of (2) Y^tViewed as an m-dimensional Gaussian distribution, i.e.

Wherein mu²Is Y^tThe covariance matrix is an m-order matrix, and the probability density function of the random variable of the observation error can be obtained:

R^t～N(0，μ²) (8)

where 0 is the mean value, μ²Is a covariance matrix of it, with Y^tThe same covariance matrix, based on sufficient sample data { y^tThe method comprises the steps of (1) and a state transition transformation matrix F (), a discretization transformation S (), a discretization inverse transformation S () established in the first step, the second step, the third step and the fourth step^-1() And an observation matrix W, and calculating to obtain an observation data residual error sequence { r^tTherein of

For the slave sample state data x^tData obtained in ascending dimensions, i.e.

Namely, it is

According to { r^tStatistical inference of μ²。

Further, performing multi-target simultaneous inversion based on the state space model, specifically, the state space model is as follows: x^t＝f^t(X^t-1)+Q^t，Y^t＝o^t(X^t)+R^tWherein X is^t、Q^t、Y^t、R^tSubject to a Gaussian distributed random vector, o^t() Taking the random variable X as a linear transformation matrix, adopting a maximum posterior method and based on a Bayesian framework to establish a random variable X^tA posterior probability density function p (X)^t|y，θ)＝c*p(y|X^t，θ)*p(X^tI θ), where a posterior probability density function p (X)^tI Y, θ) is determined, a random sequence X is established as X ═ X¹，X²，...，X^t-1，X^tThe joint posterior probability density function of: p (X | Y, θ) ═ p (X)¹，X²，...，X^t-1，X^t|Y，θ)＝p(X¹|Y，θ)*p(X²|X¹，Y，θ)*p(X³|X₂，Y，θ)*...*p(X^t|X^t-1Y, θ) solving the sequence of states by maximum a posteriori

Namely, it is

Since X is ═ X¹，X²，...，X^t-1，X^tIs a markov chain, so p (X | Y, θ) is p (X)¹|Y¹，θ)p(X²|X¹，Y²，θ)....p(X^t|X^t-1，Y^tθ), where p (X)^t|X^t-1，Y^t，θ)＝p(X^t|X^t-1，θ)p(Y^t|X^tθ), p (X | Y, θ) can be obtained as p (X)¹|θ)*p(Y¹|X¹，θ)*p(X²|X¹，θ)*p(Y²|X²，θ)*...*p(X^t|X^t-1，θ)*p(Y^t|X^tθ), can be obtained

Due to the fact that

And

equivalence ofIs obtained by

One step derivation

Wherein

J_x，k(X)＝{X^k-S^-1(F(S(X^k-1)))}ε^2-1{X^k-S^-1(F(S(X^{k -1})))}，J_y，k(X)＝{WX^k-Y^t}μ^2-1{WX^k-Y^tGet it out using convex optimization algorithm

In addition, the method also comprises a process of collecting and sorting observation data and sample data for inversion, wherein the data is divided into observation data and state data, and the sample data is a data pair with both observation data and state data. The sample data is used for training and improving an inversion algorithm, the observation data is used as the input of the inversion algorithm, the state data is the data which is not needed yet, and the state data is used as the target of inversion and may have a certain relation with the observation data. For convenience of explanation, the state data is X, and the observation data is Y. The sample data is a data pair of X and Y. X and Y are 2-dimensional arrays, each column represents a timing curve, each row represents a sampling point of the whole timing curve, the rows of X and Y are the same, but the columns can be different, and the number of columns of Y is larger than or equal to that of columns of X. Each row of X and Y represents a feature vector whose values are continuous variables over a range of values. If time domain data and depth domain number exist, time-depth calibration needs to be done in advance, and in addition, normal processing, standardization and normalization work needs to be done. The state data X and the observation data Y correspond to X, Y in expression (1) and expression (2), respectively.

The invention relates to a well-seismic joint multi-target simultaneous inversion method based on a state space model and support vector regression, which has the following characteristics: the method comprises the steps of establishing model parameters by using a machine learning method by taking a state space model as a link, combining an earthquake inversion and geological statistics method, and simultaneously performing inversion and prediction on a plurality of related target parameters, so that all target parameters meet earthquake and geological constraint conditions, and contradictions among the related target parameters are reduced to the maximum extent, so that lithological inversion and physical property prediction are coordinated and verified mutually, and useful information is mined from massive data such as earthquake, well logging, statistical prediction and the like, the inversion and prediction period is greatly shortened, the decision precision is improved, the development cost is saved, and the generation efficiency is improved.

In addition, compared with the background technology, the invention also has the following beneficial effects: the well-seismic joint multi-target simultaneous inversion method based on the state space model and the support vector regression is applied to seismic inversion, so that the difficulty of establishing the state space model can be greatly reduced, and manpower, material resources and time are saved.

Description of the drawings:

FIG. 1 is a flow chart of the method of the present invention;

FIG. 2 is a graph of measured sample curves, wherein FIG. 2(a) is longitudinal wave velocity, FIG. 2(b) is transverse wave velocity, and FIG. 2(c) is density;

FIG. 3 is a diagram of AVO forward results;

FIG. 4 is a graph showing the comparison between the AVO results of well-seismic joint multi-objective simultaneous inversion based on a state space model and support vector regression for a sample well and the actually measured sample curves, wherein FIG. 4(a) is the longitudinal wave velocity, FIG. 4(b) is the transverse wave velocity, and FIG. 4(c) is the density;

FIG. 5 is a cross-sectional view of a line of incident angles superimposed, wherein FIG. 5(a) is 5-10 degrees, FIG. 5(b) is 10-20 degrees, and FIG. 5(c) is 20-35 degrees;

fig. 6 is a diagram of the result of simultaneous inversion of a line AVO by combining multiple targets based on a state space model and a support vector regression to obtain well-seismic, where fig. 6(a) is a longitudinal wave velocity, fig. 6(b) is a transverse wave velocity, fig. 6(c) is a density, and fig. 6(d) is a longitudinal wave velocity ratio.

The specific implementation mode is as follows:

the invention is further illustrated by the following examples in conjunction with the accompanying drawings:

in order to make the purpose, technical scheme and advantages of the invention clearer, the following will take the three-dimensional actual earthquake work area of the Songliao basin in Daqing exploration area as an example, and the embodiment of the invention will be further described in detail with reference to the attached drawings.

As shown in fig. 1, a well-seismic joint multi-objective simultaneous inversion method based on state space model and support vector regression includes the following steps:

the work area has a small number of sample wells, each well has actually measured sample curve longitudinal wave velocity, transverse wave velocity and density, as shown in fig. 2, and fig. 3 is an AVO forward result of the sample curve shown in fig. 2; the work area also has a seismic data volume with incident angles of 5-10 degrees, 10-20 degrees and 20-35 degrees, for example, as shown in fig. 5, the longitudinal wave velocity, the transverse wave velocity and the density volume data of the work area are inverted when the target is inverted. It is assumed that the data has been calibrated to the time domain.

The first step is as follows: and establishing a state space model based on a K-Means discretization method.

In order to establish all state space models for inversion, firstly, discretization processing is carried out on the measured sample curve longitudinal wave velocity, transverse wave velocity and density of a sample well based on a K-Means discretization method, and { x is made^tIs sample data, where x^t∈R³，x^t ₁Is the velocity of longitudinal wave, x^t ₂Is the transverse wave velocity, x^t ₃Is the density. The targets for setting K-Means discretization have 200 types, and the 200 types are used as the state space of the state space model and are identified by 200 natural numbers from 1 to 200.

K-Means discretization sample data { x }^tThe process is as follows: 1) because sheets on object property sets in a datasetBits may be different, so x must be added before clustering^tNormalizing the continuous attribute values of the data of each dimension to mean 0 and variance 1. 2) During initialization, each object is regarded as a cluster; and calculating the distance between the clusters, selecting two clusters with the minimum distance to combine to generate a new object cluster, for example, combining two clusters to form a cluster ab when the distance between the cluster a and the cluster b is minimum. 3) And calculating the distance between the merged cluster and the rest other clusters, and then carrying out the next round of merging. When there is only one object in each cluster of objects, the distance between two clusters is defined as the euclidean distance between the objects. The distance between two clusters and another cluster after they are merged is given by the following formula d_a(bc)＝d_(bc)a＝α_bd_ab+α_cd_ac+βd_bc+γ|d_ab-d_acL wherein

γ is 0. During merging, all the clusters of classes constitute a partition of the domain of discourse U, the elements belonging to the same cluster of classes being considered indistinguishable. Therefore, we should continue the merging of the clusters until the inconsistency of the data generated after merging is higher than that of the original data, so that the clustering process is completed. And in the second stage, the clustering results are sorted and combined, wherein r is the number of generated clusters, and K is a certain cluster generated by clustering. The elements in this class cluster are in attribute A_jThe set of projections on is:

thus, the definition is according to the object at attribute A_jThe value of above defines a class cluster K_jDividing intervals

Comprises the following steps:

for a given oneDefining the domain of the class cluster K may be domain intersections defining other class clusters, in which case the domains contained in the other domains are removed and the remaining domains are sorted so as not to intersect each other. The post-processing procedure involves merging adjacent intervals to reduce the number of object clusters. The process of merging includes calculating the class entropy of all adjacent intervals on all attributes. And selecting two adjacent intervals with the minimum class entropy to be combined until the inconsistency of the data exceeds a certain threshold value. The K-Means discretization sample data process is marked as S (). Let X^t＝f^t(X^t-1)+Q^tWherein X is^t∈RⁿIs a model state, R is a real number, n is a state X^tBecause the dimension of (1) needs to simultaneously invert the longitudinal wave velocity, the transverse wave velocity and the density of the stratum, n is 3, and X is^t∈R³，

Wherein

As the velocity of the longitudinal wave,

as the velocity of the transverse wave,

is the density. To set up the state space, it is necessary to align successive states X^tDiscretizing, setting S^t∈{1、2、3...k}，k＝200，S^t＝S(X^t)，

Wherein S () is a K-Means discretization sample data process called discretization transform, S^-1() Is an inverse discretization transform. Assume sample data { x^tIs sufficient, X^t＝f^t(X^t-1)+Q^tRewritable as s^t-1＝s(X^t-1)，

{S^tIs regarded asStationary Markov chain stochastic Process, F^t() Is { s }^tThe state transition transformation matrix is a time-independent square matrix with 200 rows and 200 columns, so that F^t() Which can be written as F ().

The second step is that: and establishing a state transition transformation matrix.

The state transition transformation matrix F () is a markov chain state transition probability matrix. Assume sample data { x^tIt is sufficient that the samples are first discretized with a discretization transform S (), i.e. { S }^t}＝S({x^t}) of which s^tIs a discrete state sequence of sample data, which can be viewed as sufficient statistical data for a Markov chain state transition probability matrix F (), from { s }^tF () can be counted out in the equation). F () is a 200 th order square matrix with each element F_ijRepresenting the probability of a state transition from i to j in a markov chain.

The third step: and establishing a state residual error random variable probability density function.

According to X^t＝S^-1(F^t(S^t-1))+Q^tCan obtain the product

Wherein Q^tIs a state residual random variable, due to S^t-1Being discrete random variables, S^t＝F^t(S^t-1) Also referred to as a discrete random variable,

as state S^tCorresponding to a continuous random variable X^tExpected value of (2) X^tViewed as a k-dimensional Gaussian distribution, i.e.

Wherein epsilon²Is X^tThe covariance matrix of (2) is a k-th order square matrix because

So Q^tAlso of k-dimensional Gaussian distribution, i.e. Q^t～N(0，ε²) Wherein 0 is thereofMean value of epsilon²Is a covariance matrix of X^tThe covariance matrix of (a) is the same. According to sufficient sample data { x^tThe state transition transformation matrix F (), the discretization transformation S () and the discretization inverse transformation S established in the first step and the second step^-1() And calculating to obtain a data residual error sequence { q^tWherein q is^t＝x^t-S^-1(s^t)，s^t＝S(x^t) I.e. q^t＝x^t-S^-1(S(x^t) Namely { q) }^t}＝{x^t-S^-1(S(x^t))}. Because of q^tIs sufficient sample data due to Q^t～n(0，ε²) Can be according to { q^tConveniently, to statistically infer epsilon²。

The fourth step: and establishing an observation function based on a support vector regression method.

Support Vector Regression (SVR) is an application of SVM (support vector machine) to the regression problem. Observing function Y according to state space model^t＝o^t(X^t)+R^tWherein o is^t() Representing a state observation function, order

To obtain

{x^tIs the sequence of state data in the sample, where x^t∈R³，x^t ₁Is the velocity of longitudinal wave, x^t ₂Is the transverse wave velocity, x^t ₃Is density, { y^tIs the corresponding observed data sequence, where y^t∈R³，y^t ₁Impedance data corresponding to the data are superposed for 5-10 degrees of incidence angle beside the sample well, y^t ₂Impedance data corresponding to the data are superposed for the incidence angle of 10-20 degrees beside the sample well, y^t ₃Superimposing the impedance data corresponding to the data for a 20-35 degree angle of incidence next to the sample well, { x^tAnd { y }^tComposition of State Observation Pair sequences { (x)^t，y^t) From { (x) using SVR^t，y^t) Learn o^t(). The learning process is as follows: 1) to X^tPerforming dimension increasing treatment, namely increasing the dimension of the original random variable with 3 elements and 1 time into 3 elements and 3 times of random variables

Is shown, i.e.

Is a 20-dimensional vector; y is^tIs a 3-dimensional random variable, obviously

R^tAnd is also a 3-dimensional vector,

in the method of the invention^t() Viewed as a linear transformation, i.e. o^t() A matrix of 3 rows and 20 columns, denoted W,

referred to as W as the observation matrix. Therefore sample data { (x)^t，y^t) Can be expressed as

That is, the objective function of SVR can be established according to the sample data

Wherein

I.e. the objective function is:

where | W | is the sum of the modes of the row vectors of the matrix W, and c is the regularization factor. W is obtained by SVR learning.

The fifth step: and establishing a probability density function of the random variable of the observation error.

According to

Can obtain the product

Wherein

R^tIn order to observe the residual random variable,

as a state

Corresponding observed random variable Y^tExpected value of (2) Y^tViewed as a 3-dimensional Gaussian distribution, i.e.

Wherein mu²Is Y^tThe covariance matrix of (2) is a 3 rd order square matrix because

So R^tAlso 3-dimensional Gaussian distribution, i.e. R^t～N(0，μ²) Where 0 is the mean value, ε²Is a covariance matrix of it, with Y^tThe covariance matrix of (a) is the same. According to sufficient sample data y^tCalculating with a state transition transformation matrix F (), a discretization transformation S (), a discretization inverse transformation S-1() and an observation matrix W established in the first, second, third and fourth steps to obtain an observation data residual sequence { r }^tTherein of

For the slave sample state data x^tData obtained in ascending dimensions, i.e.

Namely, it is

Because of r^tIs sufficient sample data due to R^t～N(0，μ²) Can be according to { r^tConveniently, statistically infer mu²。

And a sixth step: and performing multi-target simultaneous inversion based on the state space model.

The state space model is: x^t＝f^t(X^t-1)+Q^t，Y^t＝o^t(X^t)+R^tThe method of the invention comprises^t、Q^t、Y^t、R^tTreated as random vectors obeying a Gaussian distribution, let o^t() The model is regarded as a linear transformation matrix, so the model is a linear Gaussian state model, according to the proof of the predecessor, when the model is used for carrying out prediction statistical inference on the linear Gaussian state model, the maximum posterior method and the Bayesian statistical inference method are equivalent, and for the sake of simplicity, the maximum posterior method is adopted in the method. Based on Bayesian framework, random variable X can be established^tA posterior probability density function p (X)^t|y，θ)＝c*p(y|X^t，θ)*p(X^tI θ), where a posterior probability density function p (X)^tI Y, θ) is determined, a random sequence X is established as X ═ X¹，X²，...，X^t-1，X^tThe joint posterior probability density function of: p (X | Y, θ) ═ p (X)¹，X²，...，X^t-1，X^t|Y，θ)＝p(X¹|Y，θ)*p(X²|X¹，Y，θ)*p(X³|X₂，Y，θ)*...*p(X^t|X^t-1Y, θ). Solving state sequences by maximum a posteriori

Namely, it is

Since X is ═ X¹，X²，...，X^t-1，X^tIs a markov chain, so p (X | Y,θ)＝p(X¹|Y¹，θ)p(X²|X¹，Y²，θ)....p(X^t|X^t-1，Y^tθ), where p (X)^t|X^t-1，Y^t，θ)＝p(X^t|X^t-1，θ)p(Y^t|X^tθ), p (X | Y, θ) can be obtained as p (X)¹|θ)*p(Y¹|X¹，θ)*p(X²|X¹，θ)*p(Y²|X²，θ)*...*p(X^t|X^t-1，θ)*p(Y^t|X^tθ), can be obtained

Due to the fact that

And

equivalent, can obtain

One step derivation

Wherein

J_x，k(X)＝{X^k-S^-1(F(S(X^{k -1})))}ε^2-1{X^k-S^-1(F(S(X^k-1)))}，J_y，k(X)＝{WX^k-Y^t}μ^2-1{WX^k-Y^t}. This is an unconstrained convexSolving the problem by using a convex optimization algorithm

FIG. 4 is a comparative display of the results of inverting AVO based on the method of the present invention on a sample well with a measured sample curve; FIG. 5 is a cross-sectional view of a line of incident angle stack, wherein FIG. 5(a) is 5-10 degrees, FIG. 5(b) is 10-20 degrees, and FIG. 5(c) is 20-35 degrees and FIG. 6 is the result of inverting a line of AVO based on the method of the present invention.

The above description is only a preferred embodiment of the present invention, and is not intended to limit the present invention, and various modifications and changes may be made to the embodiment of the present invention by those skilled in the art. Any modification, equivalent replacement, or improvement made within the spirit and principle of the present invention should be included in the protection scope of the present invention.

Claims (8)

1. A well-to-seismic joint multi-target simultaneous inversion method based on a state space model and support vector regression comprises the following steps:

the first step is as follows: establishing a state space model by using a discretization method based on the measured data of the sample oil well, wherein the measured data of the sample oil well comprises the longitudinal wave speed, the transverse wave speed and the density of a measured sample curve;

the second step is that: establishing a state transition transformation matrix, wherein firstly, discretizing measured data of a sample oil well by utilizing discretization transformation;

the third step: establishing a state residual error random variable probability density function;

the fourth step: establishing an observation function based on a support vector regression method;

the fifth step: establishing a probability density function of an observation error random variable;

and a sixth step: and performing multi-target simultaneous inversion based on the state space model.

2. Method according to claim 1, characterized in that in the first step the discretization method is a K-Means discretization method, in particular given { x [ ]^tOf state data X of sample dataTime sequence representation, discretizing sample data { x ] in clustering K-Means^tBefore, will { x }^tNormalizing continuous attribute values of all the dimensional data to a value range with the mean value of 0 and the variance of 1; the process of discretizing sample data by K-Means is recorded with S (), and X is set^t＝f^t(X^t-1)+Q^tWherein X is^t∈RⁿIs a model state, R is a real number, n is a state X^tFor setting the state space, it is necessary to set the dimensions of (A) for the continuous states X^tDiscretizing, setting S^tE.g. {1, 2, 3.. k }, k being an integer greater than 1, S^t＝S(X^t)，

Wherein S () is the K-Means discretization sample data process, S^-1() For the inverse discretization transform, it is assumed that there is sufficient sample data { x }^t}，X^t＝f^t(X^t-1)+Q^tRewritable to S^t-1＝S(X^t-1)，

{S^tIs regarded as a stationary Markov chain stochastic process, F^t() Is { S^tThe state transition transformation matrix is a time-independent transformation matrix, so F^t() The state space model can be written as F (), and after the state data X of the sample data is discretized by K-Means, the state space model can be adjusted as the following formulas (3), (4), (5) and (6):

S^t＝S(X^t) (3)

S^t＝F^t(S^t-1) (4)

X^t＝S^-1(S^t-1)+Q^t (5)

Y^t＝o^t(X^t)+R^t (6)

wherein o is^t() A matrix of 4 rows and 9 columns, denoted W, i.e.

Is called WAnd observing the matrix.

3. The method of claim 2, wherein in the first step, the compressional velocity, shear velocity, and density of the formation are simultaneously inverted, where n is 3, X^t∈R³，

Wherein

As the velocity of the longitudinal wave,

as the velocity of the transverse wave,

is the density.

4. The method according to claim 2, wherein in the second step, the state transition transformation matrix is established by assuming sufficient sample data { x }^tFirstly, discretizing the sample by discretizing transform S (), i.e. { S }^t}＝S({x^t}) of which s^tIs a discrete state sequence of sample data, for a continuous state X according to equation (3)^tDiscretizing to obtain state space {1, 2, 3.. k }, S^tE {1, 2, 3.. k } is a state of the state space, k is an integer greater than 1, and F () is a k-th order square matrix of which each element F_ijRepresents the probability of a state transition from i to j in a Markov chain, in s^tTaking the statistic data sample of the Markov chain state transition probability matrix F as a statistic data sample, and counting each element F of the F_ij＝P(s^t＝j|s^t-1I) wherein

Which represents the probability that the t-th state is j under the condition that the t-1 th state is i.

5. Method according to claim 2, wherein in the third step a state residual random variable probability density function is established, in particular according to X^t＝S^-1(F^t(S^t-1))+Q^tCan obtain the product

Wherein Q^tIs a state residual random variable, due to S^t-1Being discrete random variables, S^t＝F^t(S^t-1) Also referred to as a discrete random variable,

as state S^tCorresponding to a continuous random variable X^tExpected value of (2) X^tViewed as a k-dimensional Gaussian distribution, i.e.

Wherein epsilon²Is X^tThe covariance matrix of (2) is a k-th order square matrix because

So Q^tAnd k-dimensional Gaussian distribution, namely obtaining a state residual error random variable probability density function:

Q^t～N(0，ε²) (7)

where 0 is the mean value,. epsilon²Is a covariance matrix of X^tThe same covariance matrix, based on sufficient sample data { x^tThe state transition transformation matrix F (), the discretization transformation S () and the discretization inverse transformation S established in the first step and the second step^-1() And calculating to obtain a data residual error sequence { q^tWherein q is^t＝x^t-S^-1(s^t)，s^t＝S(x^t) I.e. q^t＝x^t-S^-1(S(x^t) Namely { q) }^t}＝{x^t-S^-1(S(x^t) Q) according to { q }^tConveniently, to statistically infer epsilon²。

6. Method according to claim 2, wherein in the fourth step an observation function is established based on a support vector regression method, in particular an observation function Y according to a state space model^t＝o^t(X^t)+R^tWherein o is^t() Representing a state observation function, order

To obtain

{x^tIs the sequence of state data in the sample, { y^tIs the corresponding observed data sequence, { x^tAnd { y }^tComposition of State Observation Pair sequences { (x)^t，y^t) Using SVR from { (x)^t，y^t) Learn o^t() The learning process is as follows: to X^tPerforming dimension increasing treatment, namely, increasing the dimension of the original random variable of n units for 1 time into a random variable of n units for 3 times; if order X^tIs a 2-dimensional random variable, i.e.

Set the upgraded random variable to

Namely, it is

Namely, it is

Is a 9-dimensional vector; is provided with Y^tIs a 4-dimensional random variable, obviously

R^tAnd is also a 4-dimensional vector,

to get o therein^t() Viewed as a linear transformation, i.e. o^t() A matrix of 4 rows and 9 columns, denoted W, i.e.

Called W the observation matrix, so sample data { (x)^t，y^t) Can be expressed as

That is, the objective function of SVR can be established according to the sample data

Wherein

I.e. the objective function is:

wherein | W | is the sum of the modes of the row vectors of the matrix W, and c is a regularization factor, and W is obtained by an SVR learning method.

7. The method according to claim 6, wherein in the fifth step the probability density function of the random variables of the observed error is established, in particular, according to

Can obtain the product

Wherein

R^tIn order to observe the residual random variable,

as a state

Corresponding observed random variable Y^tExpected value of (2) Y^tViewed as an m-dimensional Gaussian distribution, i.e.

Wherein mu²Is Y^tThe covariance matrix is an m-order matrix, and the probability density function of the random variable of the observation error can be obtained:

R^t～N(0，μ²) (8)

where 0 is the mean value, μ²Is a covariance matrix of it, with Y^tThe same covariance matrix, based on sufficient sample data { y^tThe method comprises the steps of (1) and a state transition transformation matrix F (), a discretization transformation S (), a discretization inverse transformation S () established in the first step, the second step, the third step and the fourth step^-1() And an observation matrix W, and calculating to obtain an observation data residual error sequence { r^tTherein of

For the slave sample state data x^tData obtained in ascending dimensions, i.e.

Namely, it is

According to { r^tStatistical inference of μ²。

8. The method according to claim 1, wherein the multi-objective simultaneous inversion is performed based on a state space model, in particular the state space model is: x^t＝f^t(X^t-1)+Q^t，Y^t＝o^t(X^t)+R^tWherein X is^t、Q^t、Y^t、R^tSubject to a Gaussian distributed random vector, o^t() Taking the random variable X as a linear transformation matrix, adopting a maximum posterior method and based on a Bayesian framework to establish a random variable X^tA posterior probability density function p (X)^t|y，θ)＝c*p(y|X^t，θ)*p(X^tI θ), where a posterior probability density function p (X)^tI Y, θ) is determined, a random sequence X is established as X ═ X¹，X²，...，X^t-1，X^tThe joint posterior probability density function of:

p(X|Y，θ)＝p(X¹，X²，...，X^t-1，X^t|Y，θ)

＝p(X¹|Y，θ)*p(X²|X¹，Y，θ)*p(X³|X₂，Y，θ)*...*p(X^t|X^t-1y, θ) solving the sequence of states by maximum a posteriori

Namely, it is

Since X is ═ X¹，X²，...，X^t-1，X^tIs a markov chain, so p (X | Y, θ) is p (X)¹|Y¹，θ)p(X²|X¹，Y²，θ)....p(X^t|X^t-1，Y^tθ), where p (X)^t|X^t-1，Y^t，θ)＝p(X^t|X^t-1，θ)p(Y^t|X^tθ), p (X | Y, θ) can be obtained as p (X)¹|θ)*p(Y¹|X¹，θ)*p(X²|X¹，θ)*p(Y²|X²，θ)*...*p(X^t|X^t-1，θ)*p(Y^t|X^tθ), can be obtained

Due to the fact that

And

equivalent, can obtain

Further derivation

Wherein

J_x，k(X)＝{X^k-S^-1(F(S(X^k-1)))}ε^2-1{X^k-S^-1(F(S(X^{k -1})))}，I_y，k(X)＝{WX^k-Y^t}μ^2-1{WX^k-Y^tGet it out using convex optimization algorithm

Images