Simulated annealing nonlinear pre-stack inversion method based on tectonic dip angle constraint

Abstract

The invention relates to the technical field of seismic data processing, and relates to a simulated annealing nonlinear pre-stack inversion method based on tectonic dip angle constraint. Comprising the following steps: S1, dip angle information extraction: obtaining a stratigraphic dip angle theta through a plane wave deconstruction method; s2, constructing an inversion objective function: inverting an AVO model based on a Bayesian theory framework, introducing a tectonic stratigraphic dip angle theta as prior information, and forcing an inversion result to follow preset dip angle distribution by adding a penalty term in the objective function; s3, performing inversion calculation based on the objective function: introducing an adaptive temperature adjustment strategy in the calculation process, so that the inversion calculation can widely explore a solution space in the initial stage and can perform refined search in the later stage; meanwhile, a multi-scale inversion framework is introduced in the calculation process, inversion is carried out step by step from coarse to fine, the low-scale level is used for determining a large-scale structure, and the high-scale level is used for revealing detail features. According to the invention, the geological rationality and accuracy of inversion are significantly improved.

Description

Technical Field

[0001] The invention relates to the technical field of seismic data processing in geophysics, and more specifically to a simulated annealing nonlinear prestack inversion method based on structural dip constraints. Background Art

[0002] Seismic inversion is an important geophysical technique used to infer underground rock properties, such as wave velocity and density, based on seismic data observed on the surface. Traditional seismic inversion methods often use linear or simple nonlinear methods, which often have difficulty in obtaining accurate results when dealing with complex underground structures and noisy data.

[0003] As a global optimization method, simulated annealing algorithm has been widely used in nonlinear seismic inversion. This method simulates the temperature drop and particle motion during solid annealing to find the global optimal solution of the problem. However, the traditional simulated annealing inversion algorithm does not fully consider the dip information of geological structures, which limits the accuracy and reliability of the inversion results. Summary of the invention

[0004] The purpose of the present invention is to overcome the disadvantage of low precision of inversion results in the prior art, and to provide a simulated annealing nonlinear prestack inversion method based on structural dip angle constraints, thereby effectively improving the precision of inversion results.

[0005] In order to solve the above technical problems, the technical solution adopted by the present invention is:

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

[0007] S1. Dip information extraction: The formation dip angle θ is obtained by plane wave deconstruction method;

[0008] S2. Constructing the inversion objective function: Based on the Bayesian theoretical framework, the AVO model is inverted, the structural formation dip angle θ is introduced as prior information, and a penalty term is added to 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 explore the solution space extensively in the early stage and to refine the search in the later stage. At the same time, a multi-scale inversion framework is introduced in the calculation process to perform inversion 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 detailed features.

[0010] The simulated annealing nonlinear prestack inversion method based on structural dip constraint of the present invention incorporates the key geological information of structural dip into the simulated annealing nonlinear seismic inversion framework, which significantly improves the geological rationality and accuracy of the inversion; through the adaptive temperature adjustment strategy and multi-scale inversion framework, the exploration and development in the inversion process are effectively balanced, the convergence speed is accelerated, and the inversion efficiency is improved.

[0011] Preferably, the step S1 comprises:

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

[0013]

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

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

[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 the Z transformation about space x, Z t is the Z-transformation about time t, A(Z t , Z x ) is a two-dimensional prediction filter in the time and space domains; B(Z t )、B(1 / Z t ) is a filter that adjusts the approximate accuracy, B(Z t ) / B(1 / Z t ) is the phase shift operator e obtained by Fourier transform iωσ The approximate form of the time-varying local slope, that is, time t is also a variable of the local inclination angle σ, so C(Z t , Z x ) is a function of the local inclination angle σ, denoted by C(σ);

[0018] S13. Construct the convolution of the plane wave prediction filter and the seismic data d, and establish the least squares objective function:

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

[0020] S14. By solving formula (4), the formation inclination angle θ is obtained:

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

[0022] Preferably, the step S2 comprises:

[0023] S21. Using the Bayesian principle, given the known observation data d, the posterior distribution of the model parameters 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 that the noise is Gaussian noise, the likelihood function of d is expressed as:

[0027]

[0028] Where A is the normalization constant, σ d is the noise variance; G(m) is the forward operator for the Zoeppritz equation calculation;

[0029] The posterior probability distribution is expressed as:

[0030]

[0031] In the formula, 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 Bayesian framework, the prior term has the following expression:

[0035]

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

[0037] Then the objective function is expressed as:

[0038]

[0039] Where λ is the regularization parameter;

[0040] S25. Add the lateral constraint operator of the model parameters to the prior constraints and extend it to multiple channels. The objective function is modified as follows:

[0041]

[0042] S26. For two-dimensional inclined strata, the structural constraint operator C' at each point l , C' v After obtaining the formation dip angle θ, they are obtained through coordinate transformation:

[0043]

[0044] 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:

[0045]

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

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

[0048]

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

[0050] The updating 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, E(m k ) is the objective function value of the kth iteration.

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

[0052] The inversion is performed at different scale levels l, with the initial scale being large and gradually transitioning to small scales; the model update equation for each scale level is:

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

[0054] Where, d (l)is the model update direction under scale, δ (l) is the step size at this scale.

[0055] Preferably, the step S3 specifically includes:

[0056] S31. Set the initial search temperature T0, and set the search temperature T i =T0;

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

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

[0059] S34. Determine whether to go from E based on probability P k Jump to E k+1 , the judgment criterion of P is the Metripolis criterion: if E(m k+1 )<E(m k ), then P=1, otherwise

[0060] S35. Update the number of iterations at the current temperature;

[0061] S36. If the iteration termination condition is not met, repeat steps S33 to S35 until the number of cycles is met;

[0062] S37. If the iteration requirements are met, the search temperature is lowered through an adaptive temperature adjustment strategy;

[0063] S38. Determine whether the search temperature is 0. If not, 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.

[0064] The present invention also provides a simulated annealing nonlinear prestack inversion system based on structural dip angle constraints, comprising:

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

[0066] Inversion objective function construction module: used to invert the AVO model based on the Bayesian theory framework, introduce the structural formation dip angle θ as prior information, and add a penalty term in the objective function to force the inversion result to follow the preset dip angle distribution;

[0067] 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 refine the 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 detailed features.

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

[0069] A11. The structural dip information is obtained by plane wave deconstruction method. The local plane wave model is expressed by the following partial differential equation:

[0070]

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

[0072] A12. After Fourier transform and Z transform of 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 the Z transformation about space x, Z t is the Z-transformation about time t, A(Z t , Z x ) is a two-dimensional prediction filter in the time and space domains; B(Z t )、B(1 / Z t ) is a filter that adjusts the approximate accuracy, B(Z t ) / B(1 / Z t ) is the phase shift operator e obtained by Fourier transform iωσ The approximate form of the time-varying local slope, that is, time t is also a variable of the local inclination angle σ, so C(Z t , Z x ) is a function of the local inclination angle σ, denoted by C(σ);

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

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

[0077] A14. By solving formula (4), we can obtain the formation inclination angle θ:

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

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

[0080] A21. Using the Bayesian principle, given the known observation data d, the posterior distribution of the model parameters is given by the following formula:

[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 that the noise is Gaussian, the likelihood function of d is expressed as:

[0084]

[0085] Where A is the normalization constant, σ d is the noise variance; G(m) is the forward operator for the Zoeppritz equation calculation;

[0086] The posterior probability distribution is expressed as:

[0087]

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

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

[0090]

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

[0092]

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

[0094] Then the objective function is expressed as:

[0095]

[0096] Where λ is the regularization parameter;

[0097] A25. Add the lateral constraint operator of the model parameters to the prior constraints and extend it to multiple channels. The objective function is modified as follows:

[0098]

[0099] A26. For two-dimensional inclined strata, the structural constraint operator C' at each point l , C' v After obtaining the formation dip angle θ, they are obtained through coordinate transformation:

[0100]

[0101] 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:

[0102]

[0103] Among them, θ1, θ2, …, θ nx×ny is the angle corresponding to each layer and each seismic trace position, and θ is the angle vector;

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

[0105]

[0106] Preferably, the target function calculation module calculates the template function by the following steps:

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

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

[0109] A33. Obtain new model parameters m at different scales through multi-scale inversion strategy n , calculate the objective function value E(m n , T0); the multi-scale inversion strategy includes:

[0110] The inversion is performed at different scale levels l, with the initial scale being large and gradually transitioning to small scales; the model update equation for each scale level is:

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

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

[0113] A34. Based on the probability P, determine whether to k Jump to E k+1 , the judgment criterion of P is the Metripolis criterion: if E(m k+1 )<E(m k ), then P=1, otherwise

[0114] A35. Update the number of iterations at the current temperature;

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

[0116] A37. If the iteration requirements are met, the search temperature is lowered by an adaptive temperature adjustment strategy; the adaptive temperature adjustment strategy includes:

[0117] The updating 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, E(m k ) is the objective function value of the kth iteration;

[0118] A38. Determine whether the search temperature is 0. If not, 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.

[0119] Compared with the prior art, the present invention has the following beneficial effects:

[0120] The simulated annealing nonlinear prestack inversion method based on structural dip constraint of the present invention incorporates the key geological information of structural dip into the simulated annealing nonlinear seismic inversion framework, which significantly improves the geological rationality and accuracy of the inversion. Through the adaptive temperature adjustment strategy and multi-scale inversion framework, the exploration and development in the inversion process are effectively balanced, the convergence speed is accelerated, and the inversion efficiency is improved. The hybrid optimization strategy of the present invention combines the advantages of global and local search technologies, and further enhances the accuracy and stability of the inversion results. BRIEF DESCRIPTION OF THE DRAWINGS

[0121] Figure 1 This is a flow chart of the simulated annealing nonlinear prestack inversion method based on structural dip constraints.

[0122] Figure 2For the seismic data in Example 3, from top to bottom, the 0-10° section is superimposed, the 11-20° section is superimposed, and the 21-30° section is superimposed.

[0123] Figure 3 The low-frequency initial model data in Example 3 are obtained by smoothing the real model. The upper part is the longitudinal wave velocity, the middle part is the shear wave velocity, and the lower part is the density.

[0124] Figure 4 It is the formation dip information in Example 3.

[0125] Figure 5 Schematic diagram of P- and S-wave velocities and densities inverted by the simulated annealing algorithm without adding structural dip constraints in Example 3.

[0126] Figure 6 Schematic diagram of P- and S-wave velocities and densities inverted by the simulated annealing algorithm with structural dip constraints added in Example 3. DETAILED DESCRIPTION

[0127] The present invention is further described below in conjunction with specific implementation methods. The accompanying drawings are only used for exemplary descriptions and are only schematic diagrams, not actual drawings, and cannot be understood as limiting this patent; in order to better illustrate the embodiments of the present invention, some parts of the accompanying drawings may be omitted, enlarged or reduced, and do not represent the size of the actual product; for those skilled in the art, it is understandable that some well-known structures and their descriptions in the accompanying drawings may be omitted.

[0128] The same or similar numbers in the drawings of the embodiments of the present invention correspond to the same or similar parts; in the description of the present invention, it should be understood that if the terms "upper", "lower", "left", "right" and the like indicate directions or positional relationships based on the directions or positional relationships shown in the drawings, it is only for the convenience of describing the present invention and simplifying the description, rather than indicating or implying that the device or element referred to must have a specific direction, be constructed and operated in a specific direction. Therefore, the terms describing the positional relationship in the drawings are only used for illustrative purposes and cannot be understood as limitations on this patent. For ordinary technicians in this field, the specific meanings of the above terms can be understood according to specific circumstances.

[0129] Embodiment 1

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

[0131] Step 1. Dip information extraction: The formation dip angle θ is obtained by plane wave deconstruction method.

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

[0133]

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

[0135] S12. After Fourier transform and Z transform 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] In the formula, Z x is the Z transformation about space x, Z t is the Z-transformation about time t, A(Z t ,Z x ) is a two-dimensional prediction filter in the time and space domains; B(Z t )、B(1 / Z t ) is a filter that adjusts the approximate accuracy, B(Z t ) / B(1 / Z t ) is the phase shift operator e obtained by Fourier transform iωσ The approximate form of the time-varying local slope, that is, time t is also a variable of the local inclination angle σ, so C(Z t ,Z x ) is a function of the local inclination angle σ, denoted by C(σ);

[0138] S13. Construct the convolution of the plane wave prediction filter and the seismic data d, and establish the least squares objective function:

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

[0140] S14. By solving formula (4), the formation inclination angle θ is obtained:

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

[0142] Step 2. Construct the inversion objective function: Invert the AVO model based on the Bayesian theoretical framework, introduce the structural formation dip angle θ as prior information, and add a penalty term in the objective function to force the inversion result to follow the preset dip angle distribution.

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

[0144]

[0145] 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;

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

[0147]

[0148] Where A is the normalization constant, σ d is the noise variance; G(m) is the forward operator for the Zoeppritz equation calculation;

[0149] The posterior probability distribution is expressed as:

[0150]

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

[0152] S23. Minimize the following objective function 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 covariance matrix of the model parameters respectively; m is the model data randomly generated based on the initial model;

[0157] Then the objective function is expressed as:

[0158]

[0159] Where λ is the regularization parameter;

[0160] S25. Since the objective function of the above formula represents the situation of each channel, the inversion accuracy and continuity in the lateral direction cannot be guaranteed. Therefore, equation (10) is modified; the lateral constraint operator of the model parameters is added to the prior constraints, and it is extended to multiple channels. The objective function is modified as follows:

[0161]

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

[0163] S26. For two-dimensional inclined strata, the structural constraint operator C' at each point l , C' v After obtaining the formation dip angle θ, they are obtained through coordinate transformation:

[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] Among them, θ1, θ2, …, θ nx×ny is the angle corresponding to each layer and each seismic trace position, and θ is the angle vector;

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

[0169]

[0170] Step 3. Perform inversion calculation based on the objective function: 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 refine the search in the later stage; at the same time, introduce a multi-scale inversion framework in the calculation process, and perform inversion 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 detailed features.

[0171] Adaptive temperature adjustment strategies include:

[0172] The updating 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, E(m k ) is the objective function value of the kth iteration.

[0173] The multiscale inversion strategy includes:

[0174] The inversion is performed at different scale levels l, with the initial scale being large and gradually transitioning to small scales; the model update equation for each scale level is:

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

[0176] Where, d (l) is the model update direction under scale, δ (l) is the step size at this scale.

[0177] Step 3 specifically includes:

[0178] S31. Set the initial search temperature T0, and set the 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 at different scales through multi-scale inversion strategy n , calculate the objective function value E(m n , T0);

[0181] S34. Determine whether to go from E based on probability P k Jump to E k+1 , the judgment criterion of P is the Metripolis criterion: if E(m k+1 )<E(m k ), then P=1, otherwise

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

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

[0184] S37. If the iteration requirements are met, the search temperature is lowered through an adaptive temperature adjustment strategy;

[0185] S38. Determine whether the search temperature is 0. If not, 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.

[0186] The simulated annealing nonlinear prestack inversion method based on structural dip constraint of the present invention incorporates the key geological information of structural dip into the simulated annealing nonlinear seismic inversion framework, which significantly improves the geological rationality and accuracy of the inversion. Through the adaptive temperature adjustment strategy and multi-scale inversion framework, the exploration and development in the inversion process are effectively balanced, the convergence speed is accelerated, and the inversion efficiency is improved. The hybrid optimization strategy of the present invention combines the advantages of global and local search technologies, and further enhances the accuracy and stability of the inversion results.

[0187] Embodiment 2

[0188] This embodiment provides a simulated annealing nonlinear prestack inversion system based on structural dip constraint, comprising:

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

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

[0191] A11. The structural dip information is obtained by plane wave deconstruction method. The local plane wave model is expressed by the following partial differential equation:

[0192]

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

[0194] A12. After Fourier transform and Z transform of formula (1), a plane wave prediction filter is obtained:

[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 inclination angle σ, denoted by C(σ);

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

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

[0199] A14. By solving formula (4), we can obtain the formation inclination angle θ:

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

[0201] Inversion objective function construction module: used to invert the AVO model based on the Bayesian theoretical framework, introduce the structural formation dip angle θ as prior information, and add a penalty term in the objective function to force the inversion result to follow the preset dip angle distribution.

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

[0203] A21. Using the Bayesian principle, given the known observation data d, the posterior distribution of the model parameters is given by the following formula:

[0204]

[0205] 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;

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

[0207]

[0208] Where A is the normalization constant, σ d is the noise variance;

[0209] The posterior probability distribution is expressed as:

[0210]

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

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

[0213]

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

[0215]

[0216] Where μ and Σm represent the initial model and covariance matrix of model parameters respectively;

[0217] Then the objective function is expressed as:

[0218]

[0219] Where λ is the regularization parameter;

[0220] A25. Add the lateral constraint operator of the model parameters to the prior constraints and extend it to multiple channels. The objective function is modified as follows:

[0221]

[0222] A26. For two-dimensional inclined strata, the structural constraint operator C' at each point l , C' v After obtaining the formation dip angle θ, they are obtained through coordinate transformation:

[0223]

[0224] 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:

[0225]

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

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

[0228]

[0229] 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 refine the 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 detailed features.

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

[0231] A31. Set the initial search temperature T0, and set 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 at different scales through multi-scale inversion strategy n , calculate the objective function value E(m n , T0); the multi-scale inversion strategy includes:

[0234] The inversion is performed at different scale levels l, with the initial scale being large and gradually transitioning to small scales; the model update equation for each scale level is:

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

[0236] Where, d (l) is the model update direction under scale, δ (l) is the step size at this scale;

[0237] A34. Based on the probability P, determine whether to k Jump to E k+1 , the judgment criterion of P is the Metripolis criterion: if E(m k+1 )<E(m k ), then P=1, otherwise

[0238] A35. Update the number of iterations at the current temperature;

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

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

[0241] The updating 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, E(m k ) is the objective function value of the kth iteration;

[0242] A38. Determine whether the search temperature is 0. If not, 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] Embodiment 3

[0244] This embodiment is based on the method described in Example 1 for experimental verification. Figure 2 As shown, it is the seismic data of this embodiment, from top to bottom 0-10° part superimposed section, 11-20° part superimposed section, 21-30° part superimposed section; Figure 3 As shown, it is the low-frequency initial model data of this embodiment, which is obtained by smoothing the real model. The upper part is the longitudinal wave velocity, the middle part is the shear wave velocity, and the lower part is the density; Figure 4As shown, the formation dip information of this embodiment.

[0245] Based on the above experimental data, experimental comparison is carried out, such as Figure 5 As shown in the figure, the P-wave and S-wave velocities and densities inverted by the simulated annealing algorithm without adding structural dip constraints are shown in the figure. The resolution of the result is not high and there are many discrete points. Figure 6 As shown in the figure, the P-wave and S-wave velocities and densities inverted by the simulated annealing algorithm with structural dip constraints are added. The resolution of the result is significantly improved compared to when no structural constraints are added, and the randomness is reduced, the discrete points are reduced, and the lateral continuity is stronger.

[0246] It can be seen from the above experimental results that the simulated annealing nonlinear prestack inversion method based on structural dip constraint provided by the present invention incorporates the key geological information of structural dip into the simulated annealing nonlinear seismic inversion framework, which significantly improves the geological rationality and accuracy of the inversion. Through the adaptive temperature adjustment strategy and multi-scale inversion framework, the exploration and development in the inversion process are effectively balanced, the convergence speed is accelerated, and the inversion efficiency is improved. The hybrid optimization strategy of the present invention combines the advantages of global and local search technologies to further enhance the accuracy and stability of the inversion results.

[0247] In the specific contents of the above-mentioned specific implementation methods, the various technical features can be combined in any non-contradictory manner. In order to make the description concise, not all possible combinations of the above-mentioned technical features are described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this specification.

[0248] Obviously, the above embodiments of the present invention are merely examples for clearly illustrating the present invention, and are not intended to limit the embodiments of the present invention. For those skilled in the art, other different forms of changes or modifications can be made based on the above description. It is not necessary and impossible to list all the embodiments here. Any modifications, equivalent substitutions and improvements made within the spirit and principles of the present invention should be included in the protection scope of the claims of the present invention.

Claims

1. A simulated annealing nonlinear prestack inversion method based on structural dip constraint, characterized in that: The following steps are involved: S1. Dip information extraction: The formation dip angle θ is obtained by plane wave deconstruction method; S2. Constructing the inversion objective function: Based on the Bayesian theoretical framework, the AVO model is inverted, the structural formation dip angle θ is introduced as prior information, and 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 objective function: In the calculation process, an adaptive temperature adjustment strategy is introduced to enable the inversion calculation to explore the solution space extensively in the early stage and to refine the search in the later stage. At the same time, a multi-scale inversion framework is introduced in the calculation process to perform inversion 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 detailed features.

2. The simulated annealing nonlinear prestack inversion method based on structural dip angle constraint according to claim 1 is characterized in that: The step S1 comprises: S11. The structural dip information is obtained by plane wave deconstruction method, and the local plane wave model is expressed by the following partial differential equation: Where, P(x, t) is the local plane wave field, σ(x, t) is the local slope, which is a variable with respect to space x and time t; S12. After Fourier transform and Z transform of formula (1), a plane wave prediction filter is obtained: C(Z x ,Z t )=A(Z x ,Z t )B(1 / Z x )=B(1 / Z x )-Z t B(Z x ) (2) In the formula, Z x is the Z transformation about space x, Z t is the Z-transformation about time t, A(Z t ,Z x ) is a two-dimensional prediction filter in the time and space domains; B(Z t )、B(1 / Z t ) is a filter that adjusts the approximate accuracy, B(Z t ) / B(1 / Z t ) is the phase shift operator e obtained by Fourier transform iωσ The approximate form of the time-varying local slope, that is, time t is also a variable of the local inclination angle σ, so C(Z t ,Z x ) is a function of the local inclination angle σ, denoted by C(σ); S13. Construct the convolution of the plane wave prediction filter and the seismic data d, and establish the least squares objective function: C(σ)d≈0 (3) S14. By solving formula (4), the formation inclination angle θ is obtained: θ=arctan -1 σ (4)。 3. The simulated annealing nonlinear prestack inversion method based on structural dip constraint according to claim 1, characterized in that: The step S2 comprises: S21. Using the Bayesian principle, given the known observation data d, the posterior distribution of the model parameters is given by: 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; S22. Assuming that the noise is Gaussian noise, the likelihood function of d is expressed as: Where A is the normalization constant, σ d is the noise variance; G(m) is the forward operator for the Zoeppritz equation calculation; The posterior probability distribution is expressed as: In the formula, F(m) is the prior term; S23. Minimize the following objective function to maximize the solution of the posterior distribution: S24. In the Bayesian framework, the prior term has the following expression: In the formula, μ and Σm represent the initial model and covariance matrix of the model parameters respectively; m is the model data randomly generated based on the initial model; Then the objective function is expressed as: Where λ is the regularization parameter; S25. Add the lateral constraint operator of the model parameters to the prior constraints and extend it to multiple channels. The objective function is modified as follows: S26. For two-dimensional inclined strata, the structural constraint operator C' at each point l , C' v After obtaining the formation dip angle θ, they are obtained through coordinate transformation: 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: Among them, θ1, θ2, …, θ nx×ny is the angle corresponding to each layer and each seismic trace position, and θ is the angle vector; S27. Then the objective function of the inversion becomes:

4. The simulated annealing nonlinear prestack inversion method based on structural dip constraint according to claim 3 is characterized in that: The adaptive temperature adjustment strategy in step S3 includes: The updating 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, E(m k ) is the objective function value of the kth iteration.

5. The simulated annealing nonlinear prestack inversion method based on structural dip constraint according to claim 4 is characterized in that: The multi-scale inversion strategy in step S3 includes: The inversion is performed at different scale levels l, with the initial scale being large and gradually transitioning to small scales; the model update equation for each scale level is: m (l+1) =m (l) +δ (l) d (l) (15) Where, d (l) is the model update direction under scale, δ (l) is the step size at this scale.

6. The simulated annealing nonlinear prestack inversion method based on structural dip constraint according to claim 5, characterized in that: The step S3 specifically includes: S31. Setting the initial search temperature T0, setting the search temperature T i =T0; S32. Randomly generate an initial model parameter m0 and calculate the current objective function value E(m0,T0); S33. Obtain new model parameters m at different scales through multi-scale inversion strategy n , calculate the objective function value E(m n ,T0); S34. Determine whether to go from E based on probability P k Jump to E k+1 , the judgment criterion of P is the Metripolis criterion: if E(m k+1 ) <E(m k ), then P=1, otherwise S35. Update the number of iterations at the current temperature; S36. If the iteration termination condition is not met, repeat steps S33 to S35 until the number of cycles is met; S37. If the iteration requirements are met, the search temperature is lowered through an adaptive temperature adjustment strategy; S38. Determine whether the search temperature is 0. If not, 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.

7. A simulated annealing nonlinear prestack inversion system based on structural dip constraint, characterized in that: include: Dip information extraction module: used to obtain the formation dip angle θ through the plane wave deconstruction method; Inversion objective function construction module: used to invert the AVO model based on the Bayesian theory framework, introduce the structural formation dip angle θ as prior information, and add a penalty term in the objective function to force the inversion result to follow the preset dip angle distribution; 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 refine the 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 detailed features.

8. The simulated annealing nonlinear prestack inversion system based on structural dip constraint according to claim 7, characterized in that: The dip information extraction module extracts the formation dip angle θ by the following steps: A11. The structural dip information is obtained by plane wave deconstruction method. The local plane wave model is expressed by the following partial differential equation: Where, P(x, t) is the local plane wave field, σ(x, t) is the local slope, which is a variable with respect to space x and time t; A12. After Fourier transform and Z transform of formula (1), a plane wave prediction filter is obtained: C(Z x ,Z t )=A(Z x ,Z t )B(1 / Z x )=B(1 / Z x )-Z t B(Z x )(2) In the formula, Z x is the Z transformation about space x, Z t is the Z-transformation about time t, A(Z t ,Z x ) is a two-dimensional prediction filter in the time and space domains; B(Z t )、B(1 / Z t ) is a filter that adjusts the approximate accuracy, B(Z t ) / B(1 / Z t ) is the phase shift operator e obtained by Fourier transform iωσ The approximate form of the time-varying local slope, that is, time t is also a variable of the local inclination angle σ, so C(Z t ,Z x ) is a function of the local inclination angle σ, denoted by C(σ); A13. Construct the convolution of the plane wave prediction filter and the seismic data d, and establish the least squares objective function: C(σ)d≈0(3) A14. By solving formula (4), we can obtain the formation inclination angle θ: θ=arctan -1 σ (4)。 9. The simulated annealing nonlinear prestack inversion system based on structural dip constraint according to claim 7, characterized in that: The inversion objective function construction module constructs the objective function through the following steps: A21. Using the Bayesian principle, given the known observation data d, the posterior distribution of the model parameters is given by the following formula: 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; A22. Assuming that the noise is Gaussian, the likelihood function of d is expressed as: Where A is the normalization constant, σ d is the noise variance; G(m) is the forward operator for the Zoeppritz equation calculation; The posterior probability distribution is expressed as: In the formula, F(m) is the prior term; A23. Minimize the following objective function to maximize the solution of the posterior distribution: A24. In the Bayesian framework, the prior term has the following expression: In the formula, μ and Σm represent the initial model and covariance matrix of the model parameters respectively; m is the model data randomly generated based on the initial model; Then the objective function is expressed as: Where λ is the regularization parameter; A25. Add the lateral constraint operator of the model parameters to the prior constraints and extend it to multiple channels. The objective function is modified as follows: A26. For two-dimensional inclined strata, the structural constraint operator C' at each point l , C' v After obtaining the formation dip angle θ, they are obtained through coordinate transformation: 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: In the formula, θ1, θ2, …, θ nx×ny is the angle corresponding to each layer and each seismic trace position, and θ is the angle vector; A27. Then the objective function of the inversion becomes:

10. The simulated annealing nonlinear prestack inversion system based on structural dip constraint according to claim 9, characterized in that: The target function calculation module calculates the template function by the following steps: A31. Set the initial search temperature T0, and set the search temperature T i =T0; A3 2. Randomly generate an initial model parameter m0 and calculate the current objective function value E(m0,T0); A33. Obtain new model parameters m at different scales through multi-scale inversion strategy n , calculate the objective function value E(m n , T0); the multi-scale inversion strategy includes: The inversion is performed at different scale levels l, with the initial scale being large and gradually transitioning to small scales; the model update equation for each scale level is: m (l+1) =m (l) +δ (l) d (l) (15) Where, d (l) is the model update direction under scale, δ (l) is the step size at this scale; A34. Based on the probability P, determine whether to k Jump to E k+1 , the judgment criterion of P is the Metripolis criterion: if E(m k+1 ) <E(m k ), then P=1, otherwise A35. Update the number of iterations at the current temperature; A36. If the iteration termination condition is not met, repeat steps S33 to S35 until the number of cycles is met; A37. If the iteration requirements are met, the search temperature is lowered by an adaptive temperature adjustment strategy; the adaptive temperature adjustment strategy includes: The updating 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, E(m k ) is the objective function value of the kth iteration; A38. Determine whether the search temperature is 0. If not, 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.

Images