CN120214898B - Wave impedance inversion method based on improved Bayesian algorithm - Google Patents
Wave impedance inversion method based on improved Bayesian algorithmInfo
- Publication number
- CN120214898B CN120214898B CN202510685228.7A CN202510685228A CN120214898B CN 120214898 B CN120214898 B CN 120214898B CN 202510685228 A CN202510685228 A CN 202510685228A CN 120214898 B CN120214898 B CN 120214898B
- Authority
- CN
- China
- Prior art keywords
- wave impedance
- data
- representing
- posterior distribution
- disturbance
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/282—Application of seismic models, synthetic seismograms
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/30—Analysis
- G01V1/301—Analysis for determining seismic cross-sections or geostructures
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/30—Analysis
- G01V1/306—Analysis for determining physical properties of the subsurface, e.g. impedance, porosity or attenuation profiles
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/40—Seismology; Seismic or acoustic prospecting or detecting specially adapted for well-logging
-
- G—PHYSICS
- G06—COMPUTING OR CALCULATING; COUNTING
- G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
- G06N7/00—Computing arrangements based on specific mathematical models
- G06N7/01—Probabilistic graphical models, e.g. probabilistic networks
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Remote Sensing (AREA)
- Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- General Life Sciences & Earth Sciences (AREA)
- Geology (AREA)
- Environmental & Geological Engineering (AREA)
- Acoustics & Sound (AREA)
- Theoretical Computer Science (AREA)
- Data Mining & Analysis (AREA)
- Mathematical Analysis (AREA)
- General Engineering & Computer Science (AREA)
- Mathematical Physics (AREA)
- Software Systems (AREA)
- Pure & Applied Mathematics (AREA)
- Mathematical Optimization (AREA)
- Computing Systems (AREA)
- Evolutionary Computation (AREA)
- Computational Mathematics (AREA)
- Artificial Intelligence (AREA)
- Algebra (AREA)
- Probability & Statistics with Applications (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
The application belongs to the technical field of geophysical exploration, and particularly discloses a wave impedance inversion method based on an improved Bayesian algorithm, which comprises the steps of obtaining pre-stack seismic records and logging wave impedance data of a target area; performing time-depth conversion on logging wave impedance data to obtain an actual wave impedance model and a low-frequency wave impedance model, introducing random disturbance processed by Kalman filtering into the low-frequency wave impedance model to construct an initial wave impedance model, inverting based on the initial wave impedance model, sampling inversion results by using a posterior distribution sampling mode to obtain an accepted initial wave impedance model, repeating the steps to obtain a large number of posterior distribution samples, generating inversion wave impedance sections based on the posterior distribution samples, and discussing and evaluating inversion results. The inversion method has high inversion precision, can obtain the inversion wave impedance profile with higher precision, and has good application prospect in the aspects of improving inversion resolution and uncertainty analysis.
Description
Technical Field
The application relates to the technical field of geophysical exploration, in particular to a wave impedance inversion method based on an improved Bayesian algorithm.
Background
The Bayesian algorithm is an inversion method based on the principle of probability statistics, updates prior probability according to observation data, finally obtains the distribution condition of posterior data through iteration, the method utilizes likelihood functions formed by priori information and observation data to form constraint conditions, fuses forward calculation and posterior distribution sampling methods, presumes posterior probability distribution of model parameters, and extracts statistical parameters for uncertainty analysis. The bayesian formula can be generally expressed as:
,
wherein X represents model parameters, D represents observed data; representing the prior probability distribution of the model X, and determining by prior information; representing a full probability model, representing the occurrence probability of the observed data D; representing likelihood functions representing the size of the probability of obtaining the observed data D when the model X is known; The posterior probability is expressed as the magnitude of probability of occurrence of the model parameter X under the common constraint of the prior-examination information and the observed data when the observed data D is known.
In practical inversion, due toIrrespective of the model parameters, only regularization is performed in the inversion, and therefore can be regarded as a normalization constant of 1/h, so that the sum of the integrals of the probability values over the entire region remains 1. In summary, the above formula can be rewritten as:
,
according to the prior information and the observation data, the prior distribution can be optimized and iterated, and the posterior distribution situation is finally obtained through solving, which is also the basic principle of the Bayesian inversion method.
However, the traditional Bayesian wave impedance inversion method has the problems of excessive random noise components in an inversion initial model, insufficient constraint capability of likelihood functions during inversion and the like in the actual seismic data inversion, influences the inversion effect,
Therefore, there is a need for improving the bayesian algorithm, improving the seismic data inversion effect by using the improved bayesian algorithm, and performing uncertainty analysis on the seismic data inversion result to evaluate the quality of the seismic data inversion result.
Disclosure of Invention
In view of the above, the application provides a wave impedance inversion method based on an improved Bayesian algorithm, which is used for realizing smooth analysis and processing of seismic data and logging data, and finally obtaining real underground wave impedance information by carrying out Bayesian wave impedance inversion on each item of data.
In a first aspect, an embodiment of the present application provides a wave impedance inversion method based on an improved bayesian algorithm, including:
S01, acquiring pre-stack seismic records and logging wave impedance data of a target area;
s02, performing time-depth conversion on logging wave impedance data to obtain an actual wave impedance model;
S03, introducing random disturbance processed by Kalman filtering into the low-frequency wave impedance model, and constructing an initial wave impedance model;
s04, generating prior probability distribution and likelihood function based on the initial wave impedance model to perform inversion to obtain a posterior distribution sample;
S05, sampling the posterior distribution sample by using a posterior distribution sampling mode to obtain an acceptable posterior distribution sample;
s06, repeating the steps S03-S05 until a large number of posterior distribution samples are obtained;
S07, extracting a posterior mean solution or a maximum posterior probability solution from the posterior distribution sample, generating an inversion wave impedance section, and carrying out uncertainty analysis by using a posterior standard deviation to discuss and evaluate an inversion result.
In one possible implementation, after obtaining pre-stack seismic records and log impedance data for the target region, normal distribution verification is performed on the log impedance data to ensure that the data points are located substantially near normal distribution reference lines.
In one possible implementation manner, the time-depth conversion is performed on the logging wave impedance data to obtain an actual wave impedance model, specifically, a relation is established between a depth domain and a time domain, the logging wave impedance data in the depth domain is converted to the time domain to obtain wave impedance information in the time domain, and an actual wave impedance model is obtained, wherein the calculation formula is as follows:
,
Wherein, the Representing time; Representing depth; Representing velocity information over the depth field.
In one possible implementation manner, step S03 specifically includes generating a random disturbance, performing Kalman filtering processing on the random disturbance, storing a filtering value, and obtaining a filtered disturbance after the filtering processing is completed, where the initial wave impedance model is:
,
Wherein, the Representing an initial wave impedance model; Representing a low frequency wave impedance model; Is a constant; Representing the variance of the actual wave impedance model, Representing the post-filter perturbation.
In one possible implementation manner, the Kalman filtering process is performed on the random disturbance, specifically:
step S321, a state equation and an observation equation are established;
the state equation is:
,
Wherein, the Representing random disturbance data inA state of time; Representing a state transition matrix; Representing random disturbance data in A state of time; Representing a control input matrix; Representation of A control input of time; Representation of Process noise at time;
the observation equation can be expressed as:
,
Wherein, the Representation ofThe observation data of the moment of time,Representing the observation matrix of the image of the object,Representation ofObservation noise at the moment;
step S322, predicting the current state and updating the uncertainty of the state, wherein the calculation formula is as follows:
,
,
Wherein, the Representation ofA predicted state of time; Representation of Filtering and disturbing at moment; Representation of A prediction covariance matrix of the moment; Representing covariance of process noise; Representation of The covariance matrix of the time of day,Representing state transition matricesIs a transposed matrix of (a);
step S323, updating a prediction result to obtain a filtered disturbance, wherein the calculation formula is as follows:
,
,
,
Wherein, the Representing Kalman coefficients; Representing an observation matrix Is a transposed matrix of (a); representing covariance of observed noise; Representation of A control input of time; Representation of Filtering and disturbing at moment; Representation of A covariance matrix of the time; Representing the identity matrix.
In one possible implementation, the likelihood function is:
,
Wherein, the The likelihood function is represented as a function of the likelihood,Which represents the adjustment parameters of the device,The observation data is represented by a graph of the observation data,The forward data of the model is represented,The noise weighting operator is represented by a set of values,Representation matrixIs a transposed matrix of (a);
Can be expressed as:
,
Wherein, the Representation extractionTo create a diagonal matrix of elements in the matrix,Representing the post-filter perturbation.
In one possible implementation, the energy difference corresponding to the initial wave impedance model parameter before and after iteration is calculated;
If it isThe posterior distribution sample is an acceptable posterior distribution sample;
If it is Generating a first random number between 0 and 1 by using the random function to determine the initial acceptance rateWhether the first number is larger than the first random number, if the first acceptance rate isIf the posterior distribution sample is larger than the first random number, the posterior distribution sample is an acceptable posterior distribution sample, if the initial acceptance rate is thatIf the number is smaller than the first random number, a second random number between 0 and 1 is generated again by using a random function, and the acceptance rate is judged againWhether or not the received signal is greater than the second random number, if so, the receiving rate is againGreater than a second random number, the posterior distribution sample is an acceptable posterior distribution sample, and if the posterior distribution sample is accepted againLess than a second random number, the posterior distribution sample is an unacceptable posterior distribution sample;
the accepted posterior distribution samples are stored and the unacceptable posterior distribution samples are discarded.
In one possible implementation, the energy differenceThe calculation formula of (2) is as follows:
,
Wherein, the Representing the initial wave impedance model numberEnergy after the iteration; representing the initial wave impedance model number The energy after the iteration is calculated as follows:
Wherein, the Representing a constant; representing observed data; Representation model number Forward data after the iteration; Representing a noise weighting operator; Representation matrix Is a transposed matrix of (a).
In one possible implementation, the first acceptance rate is:
,
Wherein, the Is a constant;
The re-acceptance rate is as follows:
,
Wherein, the Is a constant; And (3) representing Pearson correlation coefficients between the synthetic seismic records and the actual seismic records.
In one possible implementation, pearson correlation coefficients between a synthetic seismic record and an actual seismic recordThe calculation formula of (2) is as follows:
,
Wherein, the A data length representing a synthetic seismic record; Representing model forward data; Representing the observed data.
Compared with the prior art, the application achieves the following remarkable technical effects.
(1) Compared with the prior art that the random disturbance component in the initial inversion wave impedance model is subjected to Kalman filtering, the method has the advantages that compared with the fact that the Butterworth filtering or the Chebyshev filtering is adopted in the prior art, the method simply carries out attenuation processing on signals in a specific frequency range, the Kalman filtering is utilized to carry out filtering processing on the signals based on Bayesian theory, observation and prediction are combined, the interference of the random disturbance component is reduced to a certain extent, meanwhile, the built initial impedance model contains more effective information, and the built initial impedance model is more similar to the actual wave impedance model.
(2) Compared with the traditional method, the likelihood function only considers the model forward dataAnd actual seismic recordingThe present application considers the influence of weighting error in likelihood function and introduces noise weighting operatorThe likelihood function is enabled to fully reflect the relationship between the data and the model.
(3) In the traditional bayesian inversion method, MH MCMC sampling is often used when sampling the posterior distribution. The application expands on the basis, introduces a twice acceptance-rejection strategy when posterior distribution sampling is carried out, thereby increasing the searching range of the algorithm, realizing more comprehensive exploration of the global probability space of model parameters and improving the accuracy of inversion.
Drawings
In order to more clearly illustrate the technical solutions of the embodiments of the present application, the drawings that are needed in the embodiments will be briefly described below, it being obvious that the drawings in the following description are only some embodiments of the present application, and that other drawings may be obtained according to these drawings without inventive effort for a person skilled in the art.
FIG. 1 is a schematic flow chart of a wave impedance inversion method based on an improved Bayesian algorithm provided by an embodiment of the application;
FIG. 2 is an actual prestack seismic data provided by an embodiment of the present application;
FIG. 3 is a schematic diagram of a normal distribution test of log impedance data according to an embodiment of the present application;
FIG. 4 is a schematic diagram of a low-frequency wave impedance model according to an embodiment of the present application;
FIG. 5 is an inverted wave impedance profile generated by a posterior mean solution provided by an embodiment of the present invention;
FIG. 6 is a plot of uncertainty analysis of the inversion result of a well bypass according to an embodiment of the present invention.
Detailed Description
For a better understanding of the technical solution of the present application, the following detailed description of the embodiments of the present application refers to the accompanying drawings.
It should be understood that the described embodiments are merely some, but not all, embodiments of the application. All other embodiments, which can be made by those skilled in the art based on the embodiments of the application without making any inventive effort, are intended to be within the scope of the application.
The terminology used in the embodiments of the application is for the purpose of describing particular embodiments only and is not intended to be limiting of the application. As used in this application and the appended claims, the singular forms "a," "an," "the," and "the" are intended to include the plural forms as well, unless the context clearly indicates otherwise.
It should be understood that the term "and/or" as used herein is merely an association relationship describing the associated object, and means that three relationships may exist, for example, a and/or b, and may mean that a single first exists while a single first and a single second exist. In addition, the character "/" herein generally indicates that the front and rear associated objects are an "or" relationship.
Referring to fig. 1, a flowchart of a wave impedance inversion method based on an improved bayesian algorithm is provided in an embodiment of the present application. As shown in fig. 1, it mainly includes the following steps.
And S01, acquiring pre-stack seismic records and logging wave impedance data of a target area, and carrying out normal distribution inspection on the logging wave impedance data to ensure that data points are basically positioned near normal distribution reference lines.
Preferably, pre-stack seismic records and log impedance data for the target region are obtained through seismic exploration and logging operations.
It should be noted in particular that ensuring that the data points lie substantially near the normal distribution reference line is to verify the rationality and feasibility of applying the bayesian algorithm to the target region.
Referring to fig. 2, as shown in fig. 2, a red dashed box is an inversion target area, the longitudinal time is 1660 to 1960ms, the time sampling interval is 4ms, the number of sampling points per channel is 76, the channel interval is 25m, and the total number of channels is 31, wherein the well side channel is 1573 th channel, the main frequency of the pre-stack seismic record is 33Hz through spectrum analysis, and a rake wavelet (RICKER WAVELET) with the main frequency of 33Hz is selected in inversion to generate a synthetic seismic record. Referring to fig. 3, a schematic diagram of normal distribution inspection of log impedance data provided by an embodiment of the present application is shown in fig. 3, in which time-depth conversion is used to convert log impedance data in a depth domain into time domain, and normal distribution inspection is performed, so that except for a very small number of discrete points, the data are substantially located near a standard normal distribution curve (red dashed line in fig. 3), that is, the method of the present application is applicable to the target area.
Step S02, carrying out time-depth conversion on logging wave impedance data to obtain an actual wave impedance model, carrying out horizon interpolation and extrapolation on wave impedance information on a corresponding time sequence in the actual wave impedance model, and carrying out smoothing treatment to synthesize relevant information such as the position, the form and the like of a phase axis in the pre-stack seismic record to obtain a low-frequency wave impedance model corresponding to a target area. Referring to fig. 4, a schematic diagram of a low-frequency wave impedance model provided by an embodiment of the present application, as shown in fig. 4, is shown that amplitude information in a pre-stack seismic record has been converted into wave impedance information, and the model conforms to a basic trend of the pre-stack seismic record and contains necessary information in logging wave impedance data.
The method comprises the steps of carrying out time-depth conversion on logging wave impedance data to obtain an actual wave impedance model, and specifically, establishing a relation between a depth domain and a time domain, converting the logging wave impedance data of the depth domain to the time domain to obtain wave impedance information of the time domain, and obtaining the actual wave impedance model. The calculation formula is as follows:
,
Wherein, the Representing time; Representing depth; Representing velocity information over the depth field, may be obtained from logging data or using empirical formulas, Representation pairTaking the reciprocal.
And S03, introducing random disturbance processed by Kalman filtering into the low-frequency wave impedance model, and constructing an initial wave impedance model. The method comprises the following steps:
Step S31, generating random disturbance ,, wherein,The number of data points representing random perturbations. Taking a model consisting of 120 data points as an example,Can be expressed as:。
Step S32, performing Kalman filtering processing on the random disturbance to obtain a filtered disturbance Initial wave impedance modelExpressed as:
,
Wherein, the Representing an initial wave impedance model; Representing a low frequency wave impedance model; Is a constant; Representing the variance of the actual wave impedance model, Representing the post-filter perturbation.
Specifically, the Kalman filtering processing is performed on the random disturbance, specifically:
step S321, a state equation and an observation equation are established.
The state equation is:
,
Wherein, the Representing random disturbancesData inA state of time; Representing a state transition matrix; Representing random disturbances Data inA state of time; Representing a control input matrix; Representation of A control input of time; Representation of Process noise at time subject to a desired 0, covarianceIs a normal distribution of (c).
The observation equation can be expressed as:
,
Wherein, the Representation ofThe observation data of the moment of time,Representing the observation matrix of the image of the object,Representation ofObservation noise at time subject to expectations of 0, covariance ofIs a normal distribution of (c).
Step S322, predicting the current state and updating the uncertainty of the state.
The calculation formula is as follows:
,
,
Wherein, the Representation ofA predicted state of time; Representation of Time-of-day post-filter perturbation;Representation ofA prediction covariance matrix of the moment; Representing covariance of process noise; Representation of The covariance matrix of the time of day,Representing state transition matricesIs a transposed matrix of (a).
Step S323, updating the prediction result to obtain the disturbance after filtering. The calculation formula is as follows:
,
,
,
Wherein, the Representing Kalman coefficients; Representing an observation matrix Is a transposed matrix of (a); representing covariance of observed noise; Representation of A control input of time; Representation of Time-of-day post-filter perturbation;Representation ofA covariance matrix of the time; Representing the identity matrix.
In the above process, the values of the parameters can be specifically set according to the actual problem.
And step S04, generating prior probability distribution and likelihood function based on the initial wave impedance model to perform inversion, and obtaining a posterior distribution sample.
During inversion, a likelihood function is improved by introducing a noise weighting operator, so that the relationship between data and a model can be fully and accurately reflected. The calculation formula of the likelihood function is as follows:
,
Wherein, the The likelihood function is represented as a function of the likelihood,Which represents the adjustment parameters of the device,The observation data is represented by a graph of the observation data,The forward data of the model is represented,The noise weighting operator is represented by a set of values,Representation matrixIs a transposed matrix of (a).
Can be expressed as:
,
Wherein, the Representation extractionTo create a diagonal matrix of elements in the matrix,Representing the post-filter perturbation.
And step S05, sampling the posterior distribution sample by using a posterior distribution sampling mode to obtain an acceptable posterior distribution sample. The method comprises the following steps:
Step S51, when Bayesian inversion is performed, the initial wave impedance is modeled The energy after a number of iterations is defined as:
,
Wherein, the Representing a constant; representing observed data; Representation model number Forward data after the iteration; Representing a noise weighting operator; Representation matrix Is a transposed matrix of (a).
The energy difference corresponding to the initial wave impedance model before and after iteration is recorded as,Can be expressed as:
,
Wherein, the Representing the initial wave impedance model numberEnergy after a number of iterations.
If it isIndicating that the energy of the initial wave impedance model is reduced after modification, and the posterior distribution sample is an acceptable posterior distribution sample;
If it is Generating a first random number between 0 and 1 by using the random function to determine the initial acceptance rateWhether the first number is larger than the first random number, if the first acceptance rate isIf the posterior distribution sample is larger than the first random number, the posterior distribution sample is an acceptable posterior distribution sample, if the initial acceptance rate is thatIf the number is smaller than the first random number, a second random number between 0 and 1 is generated again by using a random function, and the acceptance rate is judged againWhether or not the received signal is greater than the second random number, if so, the receiving rate is againGreater than a second random number, the posterior distribution sample is an acceptable posterior distribution sample, and if the posterior distribution sample is accepted againLess than the second random number, the posterior distribution sample is an unacceptable posterior distribution sample.
Further, in the primary acceptance-rejection strategy, the primary acceptance rate is:
,
Wherein, the Is a constant;
further, in the secondary acceptance-rejection strategy, the re-acceptance rate is:
,
Wherein, the Is a constant; And (3) representing Pearson correlation coefficients between the synthetic seismic records and the pre-stack seismic records.
Further, the Pearson correlation coefficient between the synthetic and pre-stack seismic recordsThe calculation formula of (2) is as follows:
,
Wherein, the A data length representing a synthetic seismic record; Representing model forward data, namely, synthetic seismic records; and representing the pre-stack seismic record, namely, the observation data.
The synthetic seismic record is generated by convolution of a reflection coefficient obtained by posterior distribution samples and a Rake wavelet, and can be expressed as follows:
,
Wherein, the Representing model forward data, namely, synthetic seismic records; a data length representing a synthetic seismic record; Representing a Rake wavelet; A data length representing a Rake wavelet; representing the reflection coefficient obtained from the posterior distribution sample; representing the data length of the reflection coefficient.
Rake waveletThe expression in the time domain is:
,
Rake wavelet The expression in the frequency domain is:
,
Wherein, the A time variable representing a wavelet; a frequency variable representing a wavelet; the dominant frequency of the Rake wavelet can be set and modified as appropriate according to the information such as the dominant frequency recorded in the prestack.
Preferably, the Rake wavelet in this embodiment is a Rake wavelet with a dominant frequency of 33 Hz.
Step 52, storing the accepted posterior distribution samples and discarding the unacceptable posterior distribution samples.
Step S06, repeating steps S03-S05 until a plurality of posterior distribution samples are obtained.
And S07, extracting a posterior mean solution or a maximum posterior probability solution from the posterior distribution sample, generating an inversion wave impedance section, and carrying out uncertainty analysis by using a posterior standard deviation to discuss and evaluate an inversion result.
In this embodiment, a posterior mean solution is extracted from a posterior distribution sample, whether the variation trend between the posterior mean solution and an actual wave impedance model is the same or not is observed, the coincidence degree of the posterior mean solution and the actual wave impedance model is compared, and whether the actual wave impedance model is basically located within a 99% confidence interval of the posterior distribution sample is observed, wherein the 99% confidence interval is withinWithin the interval, where isPosterior standard deviation.
Referring to fig. 5, an inversion wave impedance section generated by a posterior mean solution provided by the embodiment of the present invention, as shown in fig. 5, is a simultaneous logging wave impedance curve, which is drawn near a well bypass, and is found to be substantially the same as the trend between the inversion results of the well bypass wave impedance. The inversion wave impedance profile is consistent with the characteristics of the logging wave impedance curve, the resolution is higher, and the inversion wave impedance profile shows better continuity and consistency in the transverse direction. From this, it is clear that the result is reasonable and reliable. In the embodiment, the inversion result can be discussed and evaluated through uncertainty analysis, so that the quality of the inversion result can be better known. Referring to fig. 6, in the uncertainty analysis curve of the inversion result of the well bypass provided by the embodiment of the invention, as shown in fig. 6, the difference between the inversion result of the 1573 rd channel and the logging impedance curve is small as a whole, the fitting degree at the position marked by the arrow is high, the trend is consistent, and most of the logging impedance curve values fall within 99% of the inversion result (i.e.Within the interval, where isPosterior standard deviation) from which it can be determined that the inversion result is reasonable and reliable.
The invention provides a wave impedance inversion technology based on an improved Bayesian algorithm, which is an important method for analyzing and processing seismic data and logging data. However, the constraint capability of the conventional Bayesian inversion algorithm is limited, and the inversion result is affected to a certain extent. The invention aims to further improve inversion precision, and improves a Bayesian wave impedance inversion algorithm, wherein the main improvement points comprise that Kalman filtering processing is carried out on random disturbance components in an initial inversion wave impedance model to enable the constructed initial impedance to be more approximate to an actual wave impedance model, a noise weighting operator is introduced into a likelihood function to strengthen inversion constraint, expansion is carried out on the basis of MH MCMC sampling in a traditional method, a twice acceptance-rejection strategy is introduced when posterior distribution data are sampled, and an learned search space is enlarged to a certain extent. For verification effect, the algorithm is applied to actual seismic data, an inversion wave impedance section with higher precision is obtained, and good prospects in the aspects of improving inversion resolution and uncertainty analysis are revealed.
The foregoing embodiments are merely for illustrating the technical solution of the present invention, but not for limiting the same, and although the present invention has been described in detail with reference to the foregoing embodiments, it will be understood by those skilled in the art that modifications may be made to the technical solution described in the foregoing embodiments, or equivalents may be substituted for some or all of the technical features thereof, without departing from the spirit of the corresponding technical solution from the scope of the technical solution of the embodiments of the present invention.
Claims (9)
1. A wave impedance inversion method based on an improved bayesian algorithm, comprising:
S01, acquiring pre-stack seismic records and logging wave impedance data of a target area;
S02, performing time-depth conversion on logging wave impedance data to obtain an actual wave impedance model;
S03, introducing random disturbance processed by Kalman filtering into the low-frequency wave impedance model, and constructing an initial wave impedance model;
S04, generating prior probability distribution and likelihood function based on the initial wave impedance model to invert to obtain a posterior distribution sample, wherein the likelihood function is as follows:
p(D|X)∝exp[p_adjust(dobs-dsyn)T·W·(dobs-dsyn))]
Where p (d|x) represents likelihood functions, p_adjust represents adjustment parameters, D obs represents observation data, D syn represents model forward data, W represents noise weighting operators, (D obs-dsyn)T represents a transposed matrix of matrix (D obs-dsyn);
W can be expressed as:
wherein diag represents extraction To create a diagonal matrix, filtered disturbance representing the filtered perturbation;
S05, sampling the posterior distribution sample by using a posterior distribution sampling mode to obtain an acceptable posterior distribution sample;
s06, repeating the steps S03-S05 until a large number of posterior distribution samples are obtained;
S07, extracting a posterior mean solution or a maximum posterior probability solution from the posterior distribution sample, generating an inversion wave impedance section, and carrying out uncertainty analysis by using a posterior standard deviation to discuss and evaluate an inversion result.
2. The improved bayesian algorithm-based wave impedance inversion method according to claim 1, wherein after obtaining pre-stack seismic records and log wave impedance data of the target area, normal distribution inspection is performed on the log wave impedance data to ensure that the data points are located substantially near normal distribution reference lines.
3. The improved Bayesian algorithm-based wave impedance inversion method as set forth in claim 1, wherein the step of performing deep conversion on the logging wave impedance data to obtain an actual wave impedance model comprises establishing a relation between a depth domain and a time domain, converting the logging wave impedance data of the depth domain to the time domain to obtain wave impedance information of the time domain, and obtaining the actual wave impedance model, wherein the calculation formula is as follows:
Where T (y) represents time, y represents depth, and V (y) represents speed information over the depth field.
4. The method of claim 1, wherein step S03 comprises generating random disturbance, performing Kalman filtering on the random disturbance, storing the filtered value, and obtaining filtered disturbance disturbance after the filtering, wherein the initial wave impedance model is:
Where c2 represents the initial wave impedance model, b represents the low frequency wave impedance model, β 2 is a constant, var (a) represents the variance of the actual wave impedance model, and filtered_ disturbance represents the post-filter perturbation.
5. The improved bayesian algorithm-based wave impedance inversion method according to claim 4, wherein said performing Kalman filtering on random disturbances is specifically:
step S321, a state equation and an observation equation are established;
the state equation is:
xi＝Axi-1+Bui-1+wi-1
Wherein x i represents the state of random disturbance random_ distufbance data at i time, A represents a state transition matrix, x i-1 represents the state of random disturbance random_ disturbance data at i-1 time, B represents a control input matrix, u i-1 represents a control input at i-1 time, and w i-1 represents process noise at i-1 time;
the observation equation can be expressed as:
zi＝Hxi+vi
Wherein Z i represents the observation data at the moment i, H represents the observation matrix, and v i represents the observation noise at the moment i;
step S322, predicting the current state and updating the uncertainty of the state, wherein the calculation formula is as follows:
Wherein, the A predicted state at the time i is represented; Filtered_ disturbance representing the filtered disturbance at time i-1; The method comprises the following steps of (1) representing a prediction covariance matrix at the moment i, Q representing covariance of process noise, P i-1 representing a covariance matrix at the moment i-1, and A T representing a transpose matrix of a state transition matrix A;
step S323, updating the prediction result to obtain filtered disturbance filtered_ disturbance, where the calculation formula is:
Wherein, C k represents Kalman coefficient, H T represents transpose matrix of observation matrix H, R represents covariance of observation noise, u i represents control input of i time; The filtered disturbance at time I is represented by filtered_ disturbance, P i is represented by the covariance matrix at time I, and I is represented by the identity matrix.
6. The improved bayesian algorithm-based wave impedance inversion method according to claim 1, wherein the energy difference deltae corresponding to the initial wave impedance model before and after parameter iteration is calculated;
if delta E is less than or equal to 0, the posterior distribution sample is an acceptable posterior distribution sample;
If delta E >0, a first random number between 0 and 1 is generated by using a random function, whether the primary acceptance rate rho (delta E) is larger than the first random number is judged, if the primary acceptance rate rho (delta E) is larger than the first random number, the posterior distribution sample is an acceptable posterior distribution sample, if the primary acceptance rate rho (delta E) is smaller than the first random number, a second random number between 0 and 1 is generated again by using the random function, whether the re-acceptance rate rho (r d) is larger than the second random number is judged, if the re-acceptance rate rho (r d) is larger than the second random number, the posterior distribution sample is an acceptable posterior distribution sample, and if the re-acceptance rate rho (r d) is smaller than the second random number, the posterior distribution sample is an unacceptable posterior distribution sample;
the accepted posterior distribution samples are stored and the unacceptable posterior distribution samples are discarded.
7. The improved bayesian algorithm based wave impedance inversion method according to claim 6 wherein the energy difference Δe is calculated as:
ΔE=E(k+1)-E(k)
Wherein E (k+1) represents the energy after the k+1th iteration of the initial wave impedance model, E (k) represents the energy after the k-th iteration of the initial wave impedance model, and the calculation formula is as follows:
Wherein, the Representing constants, d obs representing observed data, d syn (k) representing forward data after the kth iteration of the model, W representing noise weighting operators, and d obs-dsyn(k))T representing the transposed matrix of the matrix (d obs-dsyn (k)).
8. The improved bayesian algorithm based wave impedance inversion method according to claim 6, wherein said first acceptance rate is:
ρ(ΔE)=exp(μΔE)
Wherein μ is a constant;
The re-acceptance rate is as follows:
Wherein, the And r d represents the Pearson correlation coefficient between the synthetic seismic record and the actual seismic record.
9. The improved bayesian algorithm-based wave impedance inversion method according to claim 8, wherein the Pearson correlation r d between the synthetic seismic record and the actual seismic record is calculated as:
where N represents the data length of the synthetic seismic record, d syn represents model forward data, and d obs represents observation data.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202510685228.7A CN120214898B (en) | 2025-05-27 | 2025-05-27 | Wave impedance inversion method based on improved Bayesian algorithm |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202510685228.7A CN120214898B (en) | 2025-05-27 | 2025-05-27 | Wave impedance inversion method based on improved Bayesian algorithm |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN120214898A CN120214898A (en) | 2025-06-27 |
| CN120214898B true CN120214898B (en) | 2025-08-01 |
Family
ID=96099276
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202510685228.7A Active CN120214898B (en) | 2025-05-27 | 2025-05-27 | Wave impedance inversion method based on improved Bayesian algorithm |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN120214898B (en) |
Family Cites Families (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| GB2463242B (en) * | 2008-09-03 | 2012-11-07 | Statoilhydro Asa | Method of modelling a subterranean region of the earth |
| CN103760600A (en) * | 2014-01-07 | 2014-04-30 | 中国石油天然气股份有限公司 | Gas saturation inversion method |
| US10353099B2 (en) * | 2014-10-24 | 2019-07-16 | Ion Geophysical Corporation | Methods and systems for seismic inversion and related seismic data processing |
| CN107783183B (en) * | 2016-08-31 | 2019-10-29 | 中国石油化工股份有限公司 | Depth Domain seismic impedance inversion and system |
| CN106772604B (en) * | 2016-12-28 | 2019-02-15 | 中国石油化工股份有限公司 | Prestack seismic inversion method based on the fluid volume compressed coefficient |
| CN117471536A (en) * | 2022-07-19 | 2024-01-30 | 中国石油化工股份有限公司 | Pre-stack multi-parameter inversion optimization solving method and system based on MCMC |
| CN117805930A (en) * | 2022-09-30 | 2024-04-02 | 中国石油化工股份有限公司 | High-resolution inversion method based on Bayesian theory |
| CN117310802B (en) * | 2023-09-13 | 2024-06-07 | 成都捷科思石油天然气技术发展有限公司 | Depth domain reservoir seismic inversion method |
-
2025
- 2025-05-27 CN CN202510685228.7A patent/CN120214898B/en active Active
Also Published As
| Publication number | Publication date |
|---|---|
| CN120214898A (en) | 2025-06-27 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Godfrey et al. | Zero memory non‐linear deconvolution | |
| RU2300123C2 (en) | Radon's transformation of high resolution for processing of seismic data | |
| CN113887398A (en) | GPR signal denoising method based on variational modal decomposition and singular spectrum analysis | |
| CN110221256B (en) | SAR interference suppression method based on deep residual error network | |
| CN114785379B (en) | Method and system for estimating parameters of underwater sound JANUS signals | |
| CN113030880B (en) | Signal estimation method for intermittent sampling repeater interference based on ADMM | |
| CN116068619B (en) | Self-adaptive multi-order frequency dispersion surface wave pressing method, device and equipment | |
| CN112526599B (en) | Wavelet phase estimation method and system based on weighted L1 norm sparse criterion | |
| CN120214898B (en) | Wave impedance inversion method based on improved Bayesian algorithm | |
| Zhang et al. | Combining CEEMD and recursive least square for the extraction of time-varying seismic wavelets | |
| CN110749923A (en) | Deconvolution method for improving resolution based on norm equation | |
| CN118070648B (en) | Electromagnetic backscattering method, equipment and storage medium combining SOM and complex phase | |
| CN113589384A (en) | Pre-stack gather amplitude-preserving and denoising method based on signal characteristic changing along with offset distance | |
| CN111856559A (en) | Multi-channel seismic spectrum inversion method and system based on sparse Bayesian learning theory | |
| CN114866110B (en) | A Parameter Estimation Method of Frequency Hopping Signal Based on Elastic Network Model and Generalized ADMM | |
| CN120335012B (en) | Stable deconvolution method based on hybrid norm constraint and frequency band optimization | |
| Neri et al. | Bayesian iterative method for blind deconvolution | |
| Nsiri et al. | Blind submarine seismic deconvolution for long source wavelets | |
| CN119596265B (en) | Two-dimensional concave bag feature target detection method based on Bayes estimation | |
| CN118169756B (en) | Irregular grid high-dimensional seismic data reconstruction method based on rapid low-rank prediction | |
| CN119902281B (en) | Method and system for determining nonlinear isochronous stratum slice | |
| CN119575463B (en) | A prestack seismic stochastic inversion method based on frequency domain co-simulation | |
| CN116405349B (en) | An orthogonal matching pursuit channel estimation method based on correlation peak search | |
| CN118915142B (en) | A spatiotemporal joint sparse deconvolution method and system for enhancing geological body recognition accuracy | |
| CN119291778B (en) | An iterative optimization method for improving the signal-to-noise ratio of seismic data in extremely shallow layers of sandstone-type uranium deposits |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |