Simulated annealing nonlinear prestack inversion method based on structural dip angle constraint

Abstract

The present application relates to the technical field of seismic data processing, and relates to a simulated annealing nonlinear prestack inversion method based on structural dip angle constraint. The method comprises the following steps: S1. Dip angle information extraction: obtaining the stratum dip angle θ by means of a plane wave deconstruction method; S2. Inversion objective function construction: inverting an AVO model based on a Bayesian theory framework, introducing the structural stratum dip angle θ as prior information, and adding a penalty term to the objective function to force the inversion result to follow a preset dip angle distribution; S3. Inversion calculation based on the objective function: introducing an adaptive temperature adjustment strategy in the calculation process to enable the inversion calculation to extensively explore the solution space in the early stage and finely search in the later stage; meanwhile, a multi-scale inversion framework is introduced in the calculation process to perform inversion step by step from coarse to fine, and the low-scale level is used to determine the large-scale structure, and the high-scale level is used to reveal the detail features. The present application significantly improves the geological rationality and accuracy of inversion.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of seismic data processing in geophysics, and more particularly to a simulated annealing nonlinear prestack inversion method based on structural dip angle constraint. BACKGROUND

[0002] Seismic inversion is an important geophysical technique for inferring the properties of subsurface rocks, such as wave velocity, density, etc., from surface observed seismic data. Traditional seismic inversion methods often use linear or simple nonlinear methods, which often have difficulty in obtaining accurate results when dealing with complex subsurface structures and noisy data.

[0003] Simulated annealing algorithm as a global optimization method has been widely used in nonlinear seismic inversion, which finds the global optimal solution of the problem by simulating the temperature reduction and particle motion in the solid annealing process. However, the traditional simulated annealing inversion algorithm does not fully consider the dip angle information of geological structure, which limits the accuracy and reliability of the inversion results. SUMMARY

[0004] The purpose of the present application is to overcome the deficiency of low accuracy of the prior art inversion results, and to provide a simulated annealing nonlinear prestack inversion method based on structural dip angle constraint, which effectively improves the accuracy of the inversion results.

[0005] To solve the above technical problems, the technical scheme adopted by the present application is:

[0006] A simulated annealing nonlinear prestack inversion method based on structural dip angle constraint is provided, comprising the following steps:

[0007] S1. Dip angle information extraction: obtain the formation dip angle θ by plane wave deconstruction method;

[0008] S2. Constructing inversion objective function: based on the Bayesian theory framework, the AVO model is inverted, the structural formation dip angle θ is introduced as prior information, and a penalty term is added in the objective function to force the inversion result to follow the preset dip angle distribution;

[0009] S3. Inversion calculation based on objective function: in the calculation process, an adaptive temperature adjustment strategy is introduced to enable the inversion calculation to widely explore the solution space in the early stage and finely search in the later stage; at the same time, a multi-scale inversion framework is introduced in the calculation process, which is performed from coarse to fine, and the low scale level is used to determine the large scale structure, and the high scale level is used to reveal the detail characteristics.

[0010] The application discloses a simulation annealing nonlinear pre-stack inversion method based on structure dip angle constraint.

[0011] Preferably, the step S1 comprises:

[0012] S11. The structure dip angle information is obtained by a plane wave deconstruction method, and a local plane wave model is expressed by the following partial differential equation:

[0013]

[0014] In the formula, wherein, P(x, t) is a local plane wave field, sigma(x, t) is a local slope, and is a variable with respect to space x and time t;

[0015] S12. A plane wave prediction filter is obtained by performing Fourier transformation and Z transformation on the formula (1):

[0016] C(Z x , Z t )=A(Z x , Z t )B(1 / Z x )=B(1 / Z x )-Z t B(Z x ) (2)

[0017] In the formula, Z x is Z transformation with respect to space x, Z t is Z transformation with respect to time t, A(Z t , Z x ) is a two-dimensional prediction filter with respect to time and space domains; B(Z t ), B(1 / Z t ) are filters for adjusting the approximation accuracy, B(Z t ) / B(1 / Z t ) is an approximate form of a phase shift operator e iωσ obtained by Fourier transformation, for a time-varying local slope, that is, time t is also a variable of the local dip angle sigma, and therefore C(Z t , Z x ) is a function with respect to the local dip angle sigma, denoted as C(sigma);

[0018] S13. Convolution of the structure plane wave prediction filter and seismic data d is performed to establish a least square target function:

[0019] C(sigma)d approximately equal to 0 (3)

[0020] S14. Obtain the formation dip angle θ by solving equation (4):

[0021] θ = arctan -1 σ (4).

[0022] Preferably, the step S2 comprises:

[0023] S21. Using the Bayes' theorem, the posterior distribution of the model parameters given the observation data d is given by:

[0024]

[0025] where p(d|m) is the likelihood function of the data, p(m) is the prior distribution, and p(d) is the total probability of d;

[0026] S22. Assuming the noise is Gaussian, the likelihood function of d is given by:

[0027]

[0028] where A is a normalization constant, σ d is the noise variance; and G(m) is the forward operator calculated by the Zoeppritz equation;

[0029] The posterior probability distribution is given by:

[0030]

[0031] where F(m) is the prior term;

[0032] S23. Minimize the following objective function to maximize the solution of the posterior distribution:

[0033]

[0034] S24. In the Bayes' framework, the prior term has the following expression:

[0035]

[0036] where μ and Σm represent the initial model and the covariance matrix of the model parameters, respectively; and m is a model data randomly generated based on the initial model;

[0037] The objective function is given by:

[0038]

[0039] where λ is the regularization parameter;

[0040] S25. Add the lateral constraint operator of the model parameter on the prior constraint, and extend it to multiple channels, and modify the objective function as:

[0041]

[0042] S26. For two-dimensional tilted strata, the structure constraint operator C' of each point l , C' v are obtained by coordinate transformation after the strata dip angle θ is obtained:

[0043]

[0044] Suppose that the underground medium has ny layers and nx channels to be inverted, and the rotation angle coordinate of the point to be inverted is specifically written as:

[0045]

[0046] In the formula, θ1, θ2, …, θ nx×ny are angles corresponding to each layer and each seismic channel position, and θ is an angle vector;

[0047] S27. Then the objective function of inversion becomes:

[0048]

[0049] Preferably, the adaptive temperature adjustment strategy in step S3 comprises:

[0050] The update rule of the temperature parameter T is as follows: define a cooling coefficient α>1, if E(m k+1 )>E(m k ), then T k+1 =αT k , otherwise, T k+1 =T k / α, T k is the temperature of the kth iteration, and E(m k ) is the objective function value of the kth iteration.

[0051] Preferably, the multi-scale inversion strategy in step S3 comprises:

[0052] Inversion is performed at different scale levels l, the initial scale is large scale, and gradually transitions to small scale; the model update equation of each scale level is:

[0053] m (l+1) =m (l) +δ (l) d (l) (15)

[0054] In the formula, d (l)is the model updating direction in the scale, δ (l) is the step size in the scale.

[0055] Preferably, the step S3 specifically comprises:

[0056] S31. Setting an initial search temperature T0, and letting search temperature T i = T0;

[0057] S32. Randomly generating an initial model parameter m0, and calculating the current objective function value E(m0, T0);

[0058] S33. Obtaining new model parameters m n in different scales through a multi-scale inversion strategy, and calculating the objective function value E(m n , T0) of the new model;

[0059] S34. Judging whether to jump from E k to E k+1 according to the probability P, and the judging criterion of P is the Metripolis criterion: if E(m k+1 ) < E(m k ), then P = 1, otherwise P = 0;

[0060] S35. Updating the iteration number in the current temperature;

[0061] S36. If the iteration termination condition is not met, repeating steps S33-S35 until the loop number is met;

[0062] S37. If the iteration requirement is met, reducing the search temperature through an adaptive temperature adjustment strategy;

[0063] S38. Judging whether the search temperature is 0, if not, returning to step S33 to continue the loop, and if yes, ending the loop and outputting the optimal model parameter solution of the objective function.

[0064] The application also provides a simulation annealing nonlinear pre-stack inversion system based on structure dip angle constraint, comprising:

[0065] A dip angle information extraction module: used for obtaining the stratum dip angle θ through a plane wave deconstruction method;

[0066] An inversion objective function construction module: used for inverting an AVO model based on a Bayesian theory framework, introducing the structure stratum dip angle θ as prior information, and forcing the inversion result to follow a preset dip angle distribution by adding a penalty term in the objective function;

[0067] The objective function calculation module is used for introducing an adaptive temperature adjustment strategy in a calculation process, so that the inversion calculation can widely explore the solution space in the early stage and finely search in the later stage; meanwhile, a multi-scale inversion framework is introduced in the calculation process, and the inversion is performed step by step from coarse to fine, the low-scale level is used for determining large-scale structures, and the high-scale level is used for revealing detailed features.

[0068] Preferably, the dip angle information extraction module extracts the formation dip angle θ by the following steps:

[0069] A11. The dip angle information is obtained by a plane wave decomposition method, and a local plane wave model is represented by the following partial differential equation:

[0070]

[0071] In the formula, P(x, t) is a local plane wave field, σ(x, t) is a local slope, and is a variable with respect to space x and time t;

[0072] A12. After Fourier transformation and Z transformation are performed on the formula (1), a plane wave prediction filter is obtained:

[0073] C(Z x ,Z t )=A(Z x ,Z t )B(1 / Z x )=B(1 / Z x )-Z t B(Z x ) (2)

[0074] In the formula, Z x is a Z transformation with respect to space x, Z t is a Z transformation with respect to time t, A(Z t , Z x ) is a two-dimensional prediction filter with respect to time and space domains; B(Z t ), B(1 / Z t ) are filters for adjusting the approximation accuracy, B(Z t ) / B(1 / Z t ) is an approximate form of a phase shift operator e iωσ obtained by Fourier transformation, for a time-varying local slope, that is, time t is also a variable of the local dip angle σ, and therefore C(Z t , Z x ) is a function with respect to the local dip angle σ, denoted as C(σ);

[0075] A13. The convolution of the plane wave prediction filter and seismic data d is constructed, and a least square objective function is established:

[0076] C(σ)d≈0 (3)

[0077] A14. The formation dip angle θ is obtained by solving equation (4):

[0078] θ = arctan -1 σ (4).

[0079] Preferably, the inversion objective function construction module constructs the objective function by the following steps:

[0080] A21. Using the Bayes' theorem, the posterior distribution of the model parameters given the observed data d is given by:

[0081]

[0082] where p(d | m) is the likelihood function of the data, p(m) is the prior distribution, and p(d) is the total probability of d;

[0083] A22. Assuming the noise is Gaussian, the likelihood function of d is expressed as:

[0084]

[0085] where A is a normalization constant, σ d is the noise variance; G(m) is the forward operator calculated by the Zoeppritz equation;

[0086] The posterior probability distribution is expressed as:

[0087]

[0088] where F(m) is the prior term;

[0089] A23. The following objective function is minimized to maximize the solution of the posterior distribution:

[0090]

[0091] A24. In the Bayes' framework, the prior term has the following expression:

[0092]

[0093] where μ and Σm represent the initial model and the covariance matrix of the model parameters, respectively; m is the model data randomly generated based on the initial model;

[0094] The objective function is expressed as:

[0095]

[0096] where λ is the regularization parameter;

[0097] A25. A lateral constraint operator is added to the model parameters on the prior constraint, and is extended to multiple channels, and the objective function is modified as:

[0098]

[0099] A26. For two-dimensional tilted strata, the structural constraint operator C' of each point l , C' v are obtained by coordinate transformation after the strata dip angle θ is obtained:

[0100]

[0101] Suppose that the underground medium has ny layers, and there are nx channels to be inverted, and the rotation angle coordinates of the points to be inverted are specifically written as:

[0102]

[0103] wherein θ1, θ2, …, θ nx×ny are angles corresponding to each layer and each seismic trace position, and θ is an angle vector;

[0104] A27. Then the objective function of the inversion is:

[0105]

[0106] Preferably, the objective function calculation module calculates the template function through the following steps:

[0107] A31. An initial search temperature T0 is set, and the search temperature T i = T0;

[0108] A32. An initial model parameter m0 is randomly generated, and the current objective function value E(m0, T0) is calculated;

[0109] A33. New model parameters m n under different scales are obtained through a multi-scale inversion strategy, and the objective function value E(m n , T0) of the new model is calculated; the multi-scale inversion strategy comprises:

[0110] Inversion is performed at different scale levels l, the initial scale is a large scale, and gradually transitions to a small scale; the model update equation of each scale level is:

[0111] m (l+1) = m (l) + δ (l) d (l) (15)

[0112] wherein d (l) is a model update direction under the scale, and δ (l)is the step size at the scale;

[0113] A34. judging whether to jump from E k to E k+1 according to the probability P, and the judging criterion of P is Metripolis criterion: if E(m k+1 ) < E(m k ), then P = 1, otherwise P = 0.

[0114] A35. updating the iteration number at the current temperature;

[0115] A36. if the iteration termination condition is not met, repeating steps S33-S35 until the iteration number is met;

[0116] A37. if the iteration requirement is met, reducing the search temperature through an adaptive temperature adjustment strategy; the adaptive temperature adjustment strategy comprises:

[0117] The updating rule of the temperature parameter T is as follows: defining a cooling coefficient α > 1, if E(m k+1 ) > E(m k ), then T k+1 = αT k , otherwise, T k+1 = T k / α, T k is the temperature of the kth iteration, and E(m k ) is the objective function value of the kth iteration;

[0118] A38. judging whether the search temperature is 0, if not, returning to step S33 to continue the loop, and if yes, ending the loop and outputting the optimal model parameter solution of the objective function.

[0119] Compared with the prior art, the present application has the beneficial effects that:

[0120] The simulation annealing nonlinear prestack inversion method based on structure dip angle constraint of the present application incorporates the key geological information of structure dip angle into the simulation annealing nonlinear seismic inversion framework, significantly improves the geological rationality and accuracy of inversion, effectively balances exploration and development in the inversion process through the adaptive temperature adjustment strategy and the multi-scale inversion framework, accelerates the convergence speed, and improves the inversion efficiency. The hybrid optimization strategy of the present application combines the advantages of global and local search techniques, further enhances the accuracy and stability of the inversion result. BRIEF DESCRIPTION OF DRAWINGS

[0121] Figure 1 is a flowchart of the simulation annealing nonlinear prestack inversion method based on structure dip angle constraint.

[0122] Figure 2For the seismic data in Example 3, the upper 0-10° partial stack profile, the 11-20° partial stack profile, and the 21-30° partial stack profile are shown.

[0123] Figure 3 For the low-frequency initial model data in Example 3, the upper panel shows the P-wave velocity, the middle panel shows the S-wave velocity, and the lower panel shows the density, which are obtained by smoothing the true model.

[0124] Figure 4 For the stratigraphic dip information in Example 3.

[0125] Figure 5 For the P-wave and S-wave velocities and density obtained by the simulated annealing algorithm without the structural dip constraint in Example 3.

[0126] Figure 6 For the P-wave and S-wave velocities and density obtained by the simulated annealing algorithm with the structural dip constraint in Example 3. DETAILED DESCRIPTION

[0127] The application will be further described below in conjunction with specific embodiments. The accompanying drawings are only used for illustrative purposes, and the representations are only schematic diagrams, not physical drawings, and should not be understood as limiting the patent. In order to better illustrate the embodiments of the application, some components in the drawings may be omitted, enlarged or reduced, and do not represent the actual size of the product. It is understandable to those skilled in the art that some well-known structures and their descriptions in the drawings may be omitted.

[0128] The same or similar reference numerals in the drawings of the embodiments of the application correspond to the same or similar components; in the description of the application, it should be understood that the orientations or positional relationships indicated by terms such as "upper", "lower", "left", "right" are based on the orientations or positional relationships shown in the drawings, and are only for the convenience of describing the application and simplifying the description, and therefore the terms describing the positional relationships in the drawings should not be understood as limiting the patent, and the specific meanings of the above terms can be understood by those skilled in the art according to the specific circumstances.

[0129] Example One

[0130] The embodiment is a simulated annealing nonlinear prestack inversion method based on structural dip constraint, which comprises the following steps:

[0131] Step 1. Dip information extraction: obtain the stratigraphic dip θ by the plane wave decomposition method.

[0132] S11. The structural dip information is obtained by the plane wave decomposition method, and is represented by the following partial differential equation as a local plane wave model:

[0133]

[0134] where P(x,t) is the local plane wave field, σ(x,t) is the local dip, and x and t are variables with respect to space and time;

[0135] S12. After Fourier transformation and Z transformation of formula (1), a plane wave prediction filter is obtained:

[0136] C(Z x ,Z t )=A(Z x ,Z t )B(1 / Z x )=B(1 / Z x )-Z t B(Z x )(2)

[0137] where Z x is the Z transformation with respect to space x, Z t is the Z transformation with respect to time t, A(Z t ,Z x ) is a two-dimensional prediction filter with respect to time and space domains; B(Z t ) and B(1 / Z t ) are filters for adjusting the approximation accuracy, B(Z t ) / B(1 / Z t ) is an approximate form of the phase shift operator e iωσ obtained by Fourier transformation, and for the time-varying local dip, i.e., t is also a variable of the local dip σ, C(Z t ,Z x ) is a function with respect to the local dip σ, denoted as C(σ);

[0138] S13. The convolution of the plane wave prediction filter and the seismic data d is constructed to establish a least square objective function:

[0139] C(σ)d≈0 (3)

[0140] S14. The formation dip θ is obtained by solving formula (4):

[0141] θ=arctan -1 σ (4)。

[0142] Step 2. The inversion objective function is constructed: based on the Bayesian theory framework, the AVO model is inverted, the formation dip θ is introduced as prior information, and a penalty term is added in the objective function to force the inversion result to follow the preset dip distribution.

[0143] S21. It is a common strategy to invert AVO model parameters using the Bayesian theory framework. The inversion process based on the Bayesian framework is as follows. According to the Bayesian principle, the posterior distribution of the model parameters under the condition of known observation data d is given by the following formula:

[0144]

[0145] In the formula, p(d|m) is the likelihood function of the data, p(m) is the prior distribution, and p(d) is the total probability of d;

[0146] S22. Assuming that the noise is Gaussian noise, the likelihood function of d is expressed as:

[0147]

[0148] In the formula, A is a normalization constant, σ d is the noise variance; G(m) is the forward operator calculated by the Zoeppritz equation;

[0149] The posterior probability distribution is expressed as:

[0150]

[0151] In the formula, F(m) is the prior term;

[0152] S23. The following objective function is minimized to maximize the solution of the posterior distribution:

[0153]

[0154] S24. Most seismic inversion methods introduce an initial model to obtain stable inversion results. In the Bayesian framework, the prior term has the following expression:

[0155]

[0156] In the formula, μ and Σm represent the initial model and the covariance matrix of the model parameters, respectively; m is the model data randomly generated based on the initial model;

[0157] Therefore, the objective function is expressed as:

[0158]

[0159] In the formula, λ is a regularization parameter;

[0160] S25. Since the objective function of the above formula represents the case of each trace, it cannot guarantee the inversion accuracy and continuity in the lateral direction, therefore, the equation (10) is modified; the lateral constraint operator of the model parameters is added to the prior constraint, and it is extended to multiple traces, and the objective function is modified as:

[0161]

[0162] When the stratum contains large relief structure, the lateral constraint operator no longer points to the same stratum interface structure direction, resulting in the inability to track lateral continuity information. At this time, it is necessary to find a structure constraint operator consistent with the stratum structure direction.

[0163] S26. For two-dimensional inclined strata, the structure constraint operator C' l , C' v of each point is obtained by coordinate transformation after the stratum dip angle θ is obtained:

[0164]

[0165] Assuming that the underground medium has ny layers and nx channels to be inverted, the rotation angle coordinates of the points to be inverted are specifically written as:

[0166]

[0167] Where θ1, θ2, …, θ nx×ny are angles corresponding to each layer and each seismic trace position, and θ is an angle vector.

[0168] S27. Then the objective function of the inversion is:

[0169]

[0170] Step 3. Inversion calculation based on the objective function: In the calculation process, an adaptive temperature adjustment strategy is introduced to enable the inversion calculation to widely explore the solution space in the early stage and finely search in the later stage; at the same time, a multi-scale inversion framework is introduced in the calculation process, and the inversion is performed step by step from coarse to fine, and the low-scale level is used to determine the large-scale structure, and the high-scale level is used to reveal the detailed features.

[0171] The adaptive temperature adjustment strategy includes:

[0172] The update rule of the temperature parameter T is as follows: define a cooling coefficient α>1, if E(m k+1 )>E(m k ), then T k+1 =αT k , otherwise, T k+1 =T k / α, T k is the temperature of the kth iteration, and E(m k ) is the objective function value of the kth iteration.

[0173] The multi-scale inversion strategy includes:

[0174] The inversion is performed at different scale levels, starting from large scale and gradually over to small scale; the model updating equation at each scale level is:

[0175] m (l+1) = m (l) + δ (l) d (l) (15)

[0176] where d (l) is the model updating direction at scale, and δ (l) is the step size at the scale.

[0177] Step 3 specifically comprises:

[0178] S31. Set an initial search temperature T0, and let search temperature T i = T0;

[0179] S32. Randomly generate an initial model parameter m0, and calculate the current objective function value E(m0, T0);

[0180] S33. Obtain new model parameters m n at different scales by the multi-scale inversion strategy, and calculate the objective function value E(m n , T0) of the new model;

[0181] S34. Determine whether to jump from E k to E k+1 according to the probability P, and the judging criterion of P is the Metripolis criterion: if E(m k+1 ) < E(m k ), then P = 1, otherwise P = 0;

[0182] S35. Update the iteration number at the current temperature;

[0183] S36. If the iteration termination condition is not met, repeat steps S33-S35 until the loop number is met;

[0184] S37. If the iteration requirement is met, reduce the search temperature by the adaptive temperature adjustment strategy;

[0185] S38. Determine whether the search temperature is 0, if not, return to step S33 to continue the loop, if yes, end the loop and output the optimal model parameter solution of the objective function.

[0186] This invention presents a simulated annealing nonlinear pre-stack inversion method based on tectonic dip constraints. This method incorporates the crucial geological information of tectonic dip into the simulated annealing nonlinear seismic inversion framework, significantly improving the geological rationality and accuracy of the inversion. Through an adaptive temperature adjustment strategy and a multi-scale inversion framework, it effectively balances exploration and development during the inversion process, accelerating convergence and improving inversion efficiency. The hybrid optimization strategy of this invention combines the advantages of global and local search techniques, further enhancing the accuracy and stability of the inversion results.

[0187] Example 2

[0188] This embodiment provides a simulated annealing nonlinear pre-stack inversion system based on constructed tilt angle constraints, including:

[0189] Dip angle information extraction module: used to obtain the dip angle θ of the formation through the plane wave deconstruction method.

[0190] The dip angle information extraction module extracts the formation dip angle θ through the following steps:

[0191] A11. The tilt angle information is obtained through the plane wave deconstruction method, and the local plane wave model is represented by the following partial differential equation:

[0192]

[0193] In the formula, P(x,t) is the local plane wave field, σ(x,t) is the local slope, and is a variable with respect to space x and time t;

[0194] A12. A plane wave prediction filter is obtained by performing Fourier transform and Z-transform on formula (1):

[0195] C(Z x Z t )=A(Z x Z t )B(1 / Z x )=B(1 / Z x )-Z t B(Z x (2)

[0196] In the formula, C(Z) x Z t ) is a function of the local tilt angle σ, denoted as C(σ);

[0197] A13. Construct the plane wave prediction filter and the seismic data d, and establish the least squares objective function:

[0198] C(σ)d≈0 (3)

[0199] A14. The dip angle θ of the formation is obtained by solving formula (4):

[0200] θ = arctan -1 σ (4).

[0201] The inversion objective function construction module is used for inverting the AVO model based on a Bayesian theory framework, introduces a constructed stratum dip angle θ as prior information, and forces the inversion result to follow a preset dip angle distribution by adding a penalty term in the objective function.

[0202] The inversion objective function construction module constructs the objective function through the following steps:

[0203] A21. According to the Bayesian principle, the posterior distribution of the model parameter under the condition of known observation data d is given by the following formula:

[0204]

[0205] In the formula, p(d|m) is a likelihood function of data, p(m) is a prior distribution, and p(d) is a total probability of d;

[0206] A22. Assuming that the noise is Gaussian noise, the likelihood function of d is represented as:

[0207]

[0208] In the formula, A is a normalization constant, and σ d is a noise variance;

[0209] The posterior probability distribution is represented as:

[0210]

[0211] In the formula, F(m) is a prior term;

[0212] A23. The following objective function is minimized to maximize the solution of the posterior distribution:

[0213]

[0214] A24. In the Bayesian framework, the prior term has the following expression:

[0215]

[0216] In the formula, μ and Σm respectively represent an initial model and a covariance matrix of the model parameter;

[0217] The objective function is represented as:

[0218]

[0219] In the formula, λ is a regularization parameter;

[0220] A25. Add the lateral constraint operator of the model parameter on the prior constraint, and extend it to multi-channel, modify the objective function as:

[0221]

[0222] A26. For two-dimensional tilted strata, the structure constraint operator C' of each point is l , C' v are obtained by coordinate transformation after the strata dip angle θ is obtained:

[0223]

[0224] Assuming that the underground medium has ny layers and nx channels to be inverted, the rotation angle coordinate of the point to be inverted is specifically written as:

[0225]

[0226] Where θ1, θ2, …, θ nx×ny are the angles corresponding to each layer and each seismic trace position, and θ is the angle vector;

[0227] A27. Then the objective function of inversion is:

[0228]

[0229] The objective function calculation module: used to introduce an adaptive temperature adjustment strategy in the calculation process, so that the inversion calculation can widely explore the solution space in the early stage and finely search in the later stage; at the same time, a multi-scale inversion framework is introduced in the calculation process, and the inversion is performed step by step from coarse to fine, the low scale level is used to determine the large scale structure, and the high scale level is used to reveal the detail features.

[0230] The objective function calculation module calculates the template function through the following steps:

[0231] A31. Set the initial search temperature T0, and let the search temperature T i =T0;

[0232] A32. Randomly generate an initial model parameter m0, and calculate the current objective function value E(m0, T0);

[0233] A33. Obtain new model parameters m n under different scales by a multi-scale inversion strategy, and calculate the objective function value E(m n , T0) of the new model; the multi-scale inversion strategy includes:

[0234] Inversion is performed on different scale levels, the initial scale is large scale, and gradually transitions to small scale; the model update equation of each scale level is:

[0235] m (l+1) =m (l) +δ (l) d (l) (15)

[0236] In the formula, d (l) It is the direction of model update at the scale, δ (l) It is the step size at that scale;

[0237] A34. Determine whether to proceed from E based on probability P. k Leap to E k+1 The criterion for judging P is the Metripolis criterion: if E(m) k+1 ) < E(m k If the answer is yes, then P = 1; otherwise...

[0238] A35. Update the iteration count at the current temperature;

[0239] A36. If the iteration termination condition is not met, repeat steps S33 to S35 until the number of iterations is met;

[0240] A37. If the iteration requirements are met, the search temperature is reduced using an adaptive temperature adjustment strategy; the adaptive temperature adjustment strategy includes:

[0241] The update rule for the temperature parameter T is as follows: Define a cooling coefficient α > 1, if E(m k+1 )>E(m k ), then T k+1 =αT k Otherwise, T k+1 =T k / α,T k Let E(m) be the temperature of the k-th iteration. k Let be the objective function value of the k-th iteration;

[0242] A38. Determine if the search temperature is 0. If it is not 0, return to step S33 to continue the loop. If it is 0, end the loop and output the optimal model parameter solution of the objective function.

[0243] Example 3

[0244] This embodiment is based on the method described in Embodiment 1 for experimental verification, such as... Figure 2 As shown, this is the seismic data of this embodiment, with the 0-10° section overlaid from top to bottom, the 11-20° section overlaid, and the 21-30° section overlaid; as Figure 3 As shown, this is the low-frequency initial model data for this embodiment, obtained by smoothing the actual model. The top represents the P-wave velocity, the middle represents the S-wave velocity, and the bottom represents the density. Figure 4As shown in the figure, it is the stratigraphic dip information of the embodiment.

[0245] Based on the above experimental data, the experiment comparison is carried out, such as Figure 5 As shown in the figure, it is the P-wave and S-wave velocity and density inverted by the simulated annealing algorithm without adding the structural dip constraint, the resolution of the result is not high and there are many discrete points; as Figure 6 As shown in the figure, it is the P-wave and S-wave velocity and density inverted by the simulated annealing algorithm with adding the structural dip constraint, the resolution of the result is obviously improved compared with that without adding the structural constraint, the randomness is reduced, the discrete points are reduced, and the lateral continuity is stronger.

[0246] From the above experimental results, it can be seen that the simulated annealing nonlinear prestack inversion method based on the structural dip constraint provided by the application incorporates the structural dip, a key geological information, into the simulated annealing nonlinear seismic inversion framework, significantly improves the geological rationality and accuracy of the inversion, balances the exploration and development in the inversion process through the adaptive temperature adjustment strategy and the multi-scale inversion framework, speeds up the convergence speed, and improves the inversion efficiency. The hybrid optimization strategy of the application combines the advantages of global and local search techniques, further enhances the accuracy and stability of the inversion result.

[0247] In the specific content of the above specific embodiments, each technical feature can be combined arbitrarily without contradiction, and in order to make the description simple, all possible combinations of the above technical features are not described, however, as long as the combination of these technical features does not exist contradiction, it should be considered as the scope of the present application.

[0248] Obviously, the above embodiments of the application are only examples for clearly illustrating the application, and are not intended to limit the implementation modes of the application. Based on the above description, other different forms of changes or modifications can be made by those skilled in the art. Here, all the implementation modes do not need to be exhausted. Any modification, equivalent replacement and improvement made within the spirit and principle of the application should be included in the protection scope of the claims of the application.

Claims

1. A method of constructing a simulated annealing nonlinear prestack inversion method based on dip angle constraints, characterized in that, Comprising the following steps: S1. Dip information extraction: Obtain the formation dip by the plane wave destruction method ; S2. Constructing the inversion objective function: based on the Bayesian theory framework, the AVO model is inverted, and the structural stratigraphic dip angle is introduced As prior information, a penalty term is added to the objective function to force the inversion result to follow the preset dip angle distribution; S3. Inversion calculation based on the objective function: an adaptive temperature adjustment strategy is introduced in the calculation process to enable the inversion calculation to widely explore the solution space in the early stage and finely search in the later stage; meanwhile, a multi-scale inversion framework is introduced in the calculation process to perform inversion from coarse to fine in stages, with low-scale levels used to determine large-scale structures and high-scale levels used to reveal detailed features; Wherein, the step S2 comprises: S21. Using the Bayesian principle, the posterior distribution of the model parameters is given by the following formula under the condition that the observation data d is known: （5） wherein is the likelihood function of the data, p (m) is the prior distribution, p (d) is the total probability for d; S22. Assuming that the noise is Gaussian noise, the likelihood function of d is expressed as: （6） wherein A is a normalization constant, is the noise variance; G(m) is the forward operator calculated from Zoeppritz equations; The posterior probability distribution is expressed as: （7） In the formula, is a prior term; S23. Minimize the following objective function to maximize the solution of the posterior distribution: （8） S24. In the Bayesian framework, the prior term has the following expression: （9） In the formula, μ and Σm respectively represent the initial model and the covariance matrix of the model parameters; m is the model data randomly generated based on the initial model; Then the objective function is expressed as: （10） In the formula, is a regularization parameter; S25. Add the lateral constraint operator of the model parameters to the prior constraint, and extend it to multiple channels, and the objective function is modified as: （11） S26. For two-dimensional dipping strata, the structural constraint operator for each point , is obtained by coordinate transformation after the dip of the strata is determined: ​ （12） Assuming that there are ny layers of underground medium and nx channels to be inverted, the rotation angle coordinates of the points to be inverted are specifically written as: （13） wherein, θ 1, θ 2, …, θ nx×ny is the angle for each layer, each seismic trace position, θ is the angle vector; S27. Then the inversion objective function becomes: （14）。 2. The method of claim 1, wherein the method is a structure angle constraint based simulated annealing nonlinear prestack inversion method. The adaptive temperature adjustment strategy in the step S3 comprises: Temperature parameter The update rule of the temperature parameter is as follows: define a cooling coefficient > 1, if then , otherwise, , is the temperature of the th iteration, is the objective function value of the th iteration.

3. The method of claim 2, wherein the method is a structure tilt angle constrained simulated annealing nonlinear prestack inversion method. The multi-scale inversion strategy in the step S3 comprises: At different scale levels The inversion is performed on the upper equation, starting from the large scale and gradually overstepping to the small scale; the model updating equation at each scale level is: （15） wherein is the model update direction at the scale, is the step size at the scale.

4. The method of claim 3, wherein the method is a structure angle constraint based simulated annealing nonlinear prestack inversion method. The step S3 specifically comprises: S31. Set initial search temperature , let search temperature ; S32. Randomly generate an initial model parameter and compute the current objective function value ; S33. Obtain new model parameters at different scales by multi-scale inversion strategy , calculate the objective function value of the new model ; S34. In accordance with the probability to determine whether to transition from , The Metropolis criterion is used as the criterion for determining whether to transition: if then = 1, otherwise ;​ S35. Update the number of iterations at the current temperature; S36. If the iteration termination condition is not met, repeat steps S33-S35 until the loop number is met; S37. If the iteration requirement is met, reduce the search temperature through the adaptive temperature adjustment strategy; S38. Determine whether the search temperature is 0, if not, return to step S33 to continue the loop, if yes, end the loop and output the optimal model parameter solution of the objective function.

5. A system for structural dip angle constrained, simulated annealing, nonlinear prestack inversion, comprising: Comprise: Dip information extraction module: used for obtaining the formation dip by the plane wave decomposition method ; Inversion objective function construction module: used for inversion of AVO model based on Bayesian theory framework, introducing construction of stratigraphic dip angle As prior information, by increasing a penalty term in the objective function to force the inversion result to follow the preset dip angle distribution; An objective function calculation module: used to introduce an adaptive temperature adjustment strategy in the calculation process to enable the inversion calculation to widely explore the solution space in the early stage and finely search in the later stage; meanwhile, a multi-scale inversion framework is introduced in the calculation process to perform inversion from coarse to fine in stages, with low-scale levels used to determine large-scale structures and high-scale levels used to reveal detailed features; Wherein, the inversion objective function construction module constructs the objective function through the following steps: A21. Using the Bayesian principle, the posterior distribution of the model parameters is given by the following formula under the condition that the observation data d is known: （5） wherein is the likelihood function of the data, p (m) is the prior distribution, p (d) is the total probability of d; A22. Assuming that the noise is Gaussian noise, the likelihood function of d is expressed as: （6） wherein A is a normalization constant, is the noise variance; G(m) is the forward operator calculated from Zoeppritz equations; The posterior probability distribution is expressed as: （7） In the formula, is a prior term; A23. Minimize the following objective function to maximize the solution of the posterior distribution: （8） A24. In the Bayesian framework, the prior term has the following expression: （9） In the formula, μ and Σm respectively represent the initial model and the covariance matrix of the model parameters; m is the model data randomly generated based on the initial model; Then the objective function is expressed as: （10） In the formula, is a regularization parameter; A25. Add the lateral constraint operator of the model parameters to the prior constraint, and extend it to multiple channels, and the objective function is modified as: （11） A26. For two-dimensional dipping strata, the structural constraint operator for each point , is obtained by coordinate transformation after the dip of the strata is determined: ​ （12） Suppose that the underground medium has ny layers and nx channels to be inverted, and the rotation angle coordinates of the points to be inverted are specifically written as: （13） wherein, wherein, θ 1, θ 2, …, θ nx×ny is the angle corresponding to each layer, each seismic trace position, θ is the angle vector; A27. The objective function of the inversion is: （14）。 6. The construction dip angle constraint based simulated annealing nonlinear prestack inversion system of claim 5, wherein, The objective function calculation module calculates the template function by the following steps: A31. Set initial search temperature , let search temperature ; A32. Randomly generate an initial model parameter and compute the current objective function value ; A33. Obtain new model parameters at different scales by a multi-scale inversion strategy , calculate the objective function value of the new model ; The multi-scale inversion strategy includes: At different scale levels The inversion is performed on the upper equation, starting from the large scale and gradually overstepping to the small scale; the model updating equation at each scale level is: （15） wherein is the model update direction at the scale, is the step size at the scale; A34. In accordance with the probability to determine whether to jump from to , The Metropolis criterion is used as the criterion for determining whether to jump: if , then = 1, otherwise ; A35. Update the iteration number at the current temperature; A36. If the iteration termination condition is not met, repeat steps A33-A35 until the loop number is met; A37. If the iteration requirement is met, reduce the search temperature by an adaptive temperature adjustment strategy; the adaptive temperature adjustment strategy includes: Temperature parameter The update rule of temperature parameter is as follows: define a cooling coefficient >1, if then , otherwise, , is the temperature of the th iteration, is the objective function value of the th iteration; A38. Determine whether the search temperature is 0, if not, return to step S33 to continue the loop, if yes, end the loop and output the optimal model parameter solution of the objective function.

Images