Geophysical inversion method and system based on multiple constraints

Abstract

The invention discloses a geophysical inversion method and system based on multiple constraints, is applied to the field of geophysical exploration, and aims to solve the problem that an existing Bayesian inversion theory framework is easily influenced by various interferences such as random noise. According to the method, multiple mixed prior information such as the characteristics of the underground geologic structure and the spatial extension relation of the seismic data is utilized at the same time, the compressive sensing theory is combined, the transverse continuity and extension relation constraint of the seismic data channel are considered, and the multiplicity of solutions of the inversion result is reduced while the reliability of the inversion result is improved.

Description

Technical Field

[0001] The invention relates to a geophysical inversion method and system based on multiple constraints, belonging to the technical field of seismic information processing in geophysical exploration. Background Art

[0002] Inversion is one of the important technologies for seismic data processing and has been widely used in high-resolution seismic exploration. Considering that the bandwidth of seismic data is limited, the low-frequency and high-frequency components outside the effective frequency band are missing, resulting in strong multi-solutions in the inversion results. To this end, various mathematical assumptions are required, and the Bayesian inversion theory framework is used to introduce prior information into the inversion algorithm to obtain an inversion result that conforms to the prior information assumption. At present, sparse models are usually used to make assumptions about the reflection coefficient sequence, and mathematical constraints are integrated into the inversion system to reduce the multi-solutions of the inversion results.

[0003] In addition, conventional inversion methods generally only consider the constraint relationship of the vertical stratum reflection coefficient, and do not consider the lateral influence between different seismic traces. The inversion method is easily affected by various interferences such as random noise, and the inversion results have large multi-solution and uncertainty. Summary of the invention

[0004] In view of the deficiencies in the prior art, the present invention provides a geophysical inversion method and system based on multiple constraints to alleviate the multi-solution and uncertainty of the existing geophysical inversion results.

[0005] In a first aspect, the present invention provides a geophysical inversion method based on multiple constraints, comprising the following steps:

[0006] S1: obtaining pre-stack seismic data records to be processed and pre-processing the seismic data, wherein the pre-stack seismic record data includes a plurality of continuous pre-stack seismic record data;

[0007] S2: Determine a geophysical forward model based on the preprocessed seismic data;

[0008] Abbreviated as y = Gx, where d i and r i represent the i-th seismic data and the i-th reflection coefficient, respectively, and W i represents the wavelet matrix of the ith channel, K represents the number of seismic data channels; y and x represent the multi-channel seismic data column vector and the multi-channel reflection coefficient column vector respectively, and G is the multi-channel wavelet matrix operator;

[0009] S3: Based on the geophysical forward model in step S2, a hybrid priori inversion system is constructed by acquiring multiple priori information;

[0010]

[0011] Where λ1, λ2 and λ3 are adjustment constraint parameters, μ is the non-convex norm adjustment factor, ‖*‖ p represents the p-norm, D n is the geosteering representation matrix; R is the seismic information spatial representation matrix.

[0012] S4: Obtain the data set and use the low-frequency seismic data to obtain the geosteering matrix D n ,At the same time, the inversion algorithm is used to obtain the seismic information space representation operator, and the Helix mathematical transformation is used to construct the seismic information space representation matrix R;

[0013] S5: inputting the pre-stack seismic data set to be processed into the inversion model, solving the target functional by using the DCA algorithm and the alternating direction multiplier method, and outputting the seismic data inversion result.

[0014] Furthermore, in the above step S5, the DCA algorithm and the alternating direction multiplier method are used to solve the above target functional to obtain the seismic data inversion result, and the specific process is as follows:

[0015] Considering that the hybrid prior inversion system is a nonlinear target functional, it cannot be directly solved by gradient algorithms. Therefore, the DCA algorithm is first used to decompose it.

[0016] φ(x)=H(x)-F(x)

[0017] In the formula,

[0018] Furthermore, the solution of the inversion system can be obtained by iteratively solving the following two sub-problems:

[0019]

[0020] Where u is the intermediate variable that relates the two sub-problems.

[0021] The first sub-problem is the linear inversion problem, which can be solved by the gradient method.

[0022]

[0023] The second sub-problem is the L1 norm constraint problem, which can be solved by the alternating direction multiplier method. Its iterative formula can be expressed as:

[0024]

[0025] In this solution, the DCA algorithm is first used to decompose the complex inversion system into two commonly used inversion systems, and then the alternating direction multiplier method is used to solve the sub-problems, which improves the efficiency and accuracy of the inversion and obtains good output seismic data inversion results.

[0026] The present invention also provides a geophysical inversion system based on multiple constraints, characterized in that the inversion system is implemented by the above-mentioned geophysical inversion method based on multiple constraints.

[0027] The beneficial effects of the present invention are as follows: the present invention simultaneously utilizes the two prior information of geology and seismology, and combines with the compressed sensing theory, and not only has stronger sparse constraint conditions, but also has better lateral continuity and extension relationship constraints of seismic data, which can improve the reliability of the inversion results while reducing the multi-solution of the inversion results, can provide key technical support for high-resolution seismic exploration, and has broad application prospects in precise seismic exploration. BRIEF DESCRIPTION OF THE DRAWINGS

[0028] In order to more clearly illustrate the specific implementation methods of the present invention or the technical solutions in the prior art, the drawings required for use in the specific implementation methods or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are some implementation methods of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative work.

[0029] Figure 1 A schematic diagram of constructing a non-convex constraint term of a reflection coefficient according to the present invention;

[0030] Figure 2 A schematic diagram of geosteering prior information of the present invention;

[0031] Figure 3 It is a schematic diagram of the seismic information space characterization operator of the present invention. DETAILED DESCRIPTION

[0032] The following will be combined with the drawings in the embodiments of the present application to clearly and completely describe the technical solutions in the embodiments of the present application. Obviously, the described embodiments are only part of the embodiments of the present application, not all of the embodiments. Based on the embodiments in the present application, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of this application.

[0033] In a first aspect, the present invention provides a geophysical inversion method based on multiple constraints, comprising the following steps:

[0034] S1: obtaining pre-stack seismic data records to be processed and pre-processing the seismic data, wherein the pre-stack seismic record data includes a plurality of continuous pre-stack seismic record data;

[0035] First, the pre-stack seismic data is subjected to high-precision denoising processing to protect the characteristic information of the effective signal as much as possible while attenuating the seismic noise. Then, the pre-stack seismic data is subjected to true amplitude recovery and trace equalization processing to enhance the energy of deep seismic signals and equalize the energy between trace sets. Next, the seismic data is subjected to static correction processing to eliminate the influence of surface factors and restore the original characteristics of the seismic signal. Subsequently, the seismic data is subjected to high-precision dynamic correction processing to minimize the influence of dynamic correction stretching and excision at far offsets. Finally, the pre-stack data is subjected to migration processing.

[0036] S2: Determine a geophysical forward model based on the preprocessed seismic data;

[0037] Abbreviated as y = Gx, where d i and r i represent the i-th seismic data and the i-th reflection coefficient, respectively, and W i represents the wavelet matrix of the ith channel, K represents the number of seismic data channels; y and x represent the multi-channel seismic data column vector and the multi-channel reflection coefficient column vector respectively, and G is the multi-channel wavelet matrix operator;

[0038] Among them, seismic data can refer to data and information collected in seismic exploration, and the seismic data includes but is not limited to seismic records at multiple angles, wavelet data, well logging data, and stratigraphic data, so as to facilitate the establishment of a seismic forward model. In the embodiment of the present invention, a geophysical forward model is constructed based on seismic data, and geophysical inversion is performed using seismic data.

[0039] S3: Based on the geophysical forward model in step S2, a hybrid priori inversion model is constructed by acquiring multiple priori information;

[0040] Conventional inversion methods generally only consider the constraint relationship of vertical stratum reflection coefficients, and do not consider the lateral influence between different seismic traces. Inversion methods are easily affected by various interferences such as random noise, and the inversion results have large multiple solutions and uncertainty. Therefore, in order to jointly use low-rank priors and inversion network priors to remove noise, the following hybrid prior optimization model is proposed:

[0041]

[0042] Where λ1, λ2 and λ3 are adjustment constraint parameters, μ is the non-convex norm adjustment factor, ||*|| p represents the p-norm, D nis the horizontal representation matrix of seismic information; R is the spatial representation matrix of seismic information; using Helix mathematical transformation, the spatial representation matrix of seismic information can be constructed.

[0043] In order to ensure the stability and accuracy of the reflection coefficient inversion process in seismic data processing, it is necessary to introduce a priori information constraints on the reflection coefficient. This involves the use of non-convex norm regularization methods in compressed sensing theory, which controls the sparsity of the reflection coefficient by imposing specific constraints during the optimization process.

[0044] Specifically, by constructing a non-convex constraint term, the sparsity of the reflection coefficient can be effectively controlled to ensure the well-posedness of the reflection coefficient inversion process. Figure 1 As shown, the non-convex norm in the field of compressed sensing is used to construct the non-convex constraint term of the reflection coefficient:

[0045] φ1(x)=||x||1-μ||x||2

[0046] In the formula, ||*|| p represents the p-norm, μ represents the non-convex norm adjustment factor;

[0047] Furthermore, the non-convex norm adjustment factor is set to μ = 0.5.

[0048] Reference Figure 2 As shown, considering that seismic data records the characteristics of underground geological structures, such as faults, folds, lithology changes, etc. The extraction of these features usually relies on the preprocessing and interpretation of seismic data, such as the analysis of attributes such as amplitude, frequency, and phase.

[0049] S4: Obtain the data set and use the low-frequency seismic data to obtain the geosteering matrix D n ,At the same time, the inversion algorithm is used to obtain the seismic information space representation operator, and the Helix mathematical transformation is used to construct the seismic information space representation matrix R;

[0050] The data set is artificially synthesized seismic data or actual collected seismic data. The data set is divided into a training set, a validation set and a test set, and the geophysical model is trained and validated based on the data set.

[0051] On the basis of extracting geological structure characteristics, a spatial gradient matrix is ​​constructed to describe the changes of geological structure in time and space. The spatial gradient matrix is ​​used to construct the geological steering regularization constraint term. The spatial gradient matrix is ​​used to construct the geological steering regularization constraint term.

[0052]

[0053] Where D nIt is the geosteering matrix parallel to the structural direction;

[0054] Furthermore, the spatial gradient matrix is ​​calculated.

[0055]

[0056] Where L is the spatial gradient operator, dt and dx represent the gradients along the time direction and spatial direction respectively.

[0057] According to matrix analysis theory, the eigenvalues ​​and eigenvectors of the spatial gradient matrix L can be expressed as,

[0058]

[0059] Where m and n represent the eigenvectors of the spatial gradient matrix L, μ m and μ n represents the eigenvalues ​​corresponding to eigenvectors m and n, m1 and m2 represent the first and second elements of eigenvector m, and n1 and n2 represent the first and second elements of eigenvector n. Assume that μ m ≥μ n ≥0, then m represents the vector perpendicular to the construction direction, and n represents the vector parallel to the construction direction.

[0060] Using the eigenvector n parallel to the structural direction, a geosteering matrix parallel to the geological trend is constructed.

[0061] D n =N1D t +N2D x

[0062] In the formula, represents the diagonal matrix constructed by the first element n1 of the eigenvector n, J and K represent the number of time sampling points and spatial seismic traces of geological data, represents the diagonal matrix constructed from the second element n2 of the eigenvector n, represents the first-order derivative matrix in the time direction, Represents the first-order derivative matrix with respect to the spatial direction.

[0063] Reference Figure 3 As shown in the figure, in order to ensure the continuity and integrity of the data during the acquisition and processing of seismic data, the spatial extension relationship of seismic waves propagating in the underground medium must be considered. This relationship not only involves the propagation of seismic waves in the vertical direction, but also includes the continuity in the horizontal direction. In order to enhance the continuity of seismic data, the information of the adjacent seismic data channels on the left and right sides of the target seismic data can be introduced, and the seismic information space characterization operator can be used to construct the seismic reflection regularization constraint term.

[0064]

[0065] Where R is the seismic information spatial representation matrix; based on the seismic information spatial representation operator b i,j , using Helix mathematical transformation, the seismic information space representation matrix R can be constructed.

[0066] Calculate the seismic information space representation operator b i,j ,

[0067] b -2,-2 b -1,-2 × b 1,-2 b 2,-2

[0068] b -2,-1 b -1,-1 × b 1,-1 b 2,-1

[0069] b -2,0 b -1,0 0 b 1,0 b 2,0

[0070] b -2,1 b -1,1 × b 1,1 b 2,1

[0071] b -2,2 b -1,2 × b 1,2 b 2,2

[0072] In the formula, the middle "0" indicates the position of the point to be predicted, "×" indicates the point not involved in the calculation, and b i,j Represents the coefficient of the spatial characterization operator of the seismic signal, i represents the time direction, and j represents the space direction.

[0073] Apply the seismic signal space characterization operator to the seismic data to obtain the seismic data prediction results.

[0074]

[0075] In the formula, s J,K represents the original 2D seismic data, represents the predicted seismic data, J and K represent the number of time sampling points and the number of spatial seismic channels.

[0076] During the inversion process, the adjustment direction and step size of the model parameters are continuously updated, and the inversion model is verified and optimized by minimizing the residual functional of the original seismic data and the predicted seismic data. The residual functional of minimizing the original seismic data and the predicted seismic data can be expressed as:

[0077]

[0078] Among them, b i,j is the spatial representation operator of the seismic signal, s J,K represents the original 2D seismic data, represents the predicted seismic data, J and K represent the number of time sampling points and the number of spatial seismic channels.

[0079] Furthermore, the Adam optimization algorithm can be used to adjust the network weights.

[0080] S5: inputting the pre-stack seismic data set to be processed into the inversion model, solving the target functional by using the alternating multiplier direction method, and outputting the seismic data inversion result.

[0081] Furthermore, in the above step S5, the DCA algorithm and the alternating direction multiplier method are used to solve the above target functional to obtain the seismic data inversion result, and the specific process is as follows:

[0082] Considering that the hybrid prior inversion system is a nonlinear target functional, it cannot be directly solved by gradient algorithms. Therefore, the DCA algorithm is first used to decompose it.

[0083] φ(x)=H(x)-F(x)

[0084] In the formula,

[0085] Furthermore, the solution of the inversion system can be obtained by iteratively solving the following two sub-problems:

[0086]

[0087] Where u is the intermediate variable that relates the two sub-problems.

[0088] The first sub-problem is the linear inversion problem, which can be solved by the gradient method.

[0089]

[0090] The second sub-problem is the L1 norm constraint problem, which can be solved by the alternating direction multiplier method. Its iterative formula can be expressed as:

[0091]

[0092] In this solution, the DCA algorithm is first used to decompose the complex inversion system into two commonly used inversion systems, and then the alternating direction multiplier method is used to solve the sub-problems, which improves the efficiency and accuracy of the inversion and obtains good output seismic data inversion results.

[0093] The method of the present invention not only has stronger sparse reflection coefficient constraints, but also has better spatial continuity and extension relationship constraints. It can improve the reliability of the inversion results while reducing the multi-solution of the inversion results, and can provide key technical support for high-resolution seismic exploration.

[0094] The present invention also provides a geophysical inversion system based on multiple constraints, characterized in that the inversion system is implemented by the above-mentioned geophysical inversion method based on multiple constraints.

[0095] The present invention provides a computer storage medium, wherein the computer storage medium stores at least one instruction, wherein the instruction is suitable for being loaded by a processor and executing the method steps of one or more embodiments of the present specification.

[0096] The present invention provides a computer program product, wherein the computer program product stores at least one instruction, wherein the instruction is suitable for being loaded by a processor and executing the method steps of one or more embodiments of the present specification.

[0097] The present invention provides an electronic device, which may include: a processor and a memory; wherein the memory stores a computer program, and the computer program is suitable for being loaded by the processor and executing the method steps of one or more embodiments of the present specification.

[0098] The above embodiments only express several implementation methods of the present invention, and the descriptions thereof are relatively specific and detailed, but they cannot be understood as limiting the scope of the present invention. It should be pointed out that, for those of ordinary skill in the art, several modifications and improvements can be made without departing from the concept of the present invention, and these all belong to the protection scope of the present invention. Therefore, the protection scope of the present invention shall be subject to the attached claims.

Claims

1. A geophysical inversion method based on multiple constraints, characterized in that: The following steps are involved: S1: obtaining pre-stack seismic data records to be processed and pre-processing the seismic data, wherein the pre-stack seismic record data includes a plurality of continuous pre-stack seismic record data; S2: Determine a geophysical forward model based on the preprocessed seismic data; Abbreviated as y=Gx, Where, d i and r i represent the i-th seismic data and the i-th reflection coefficient, respectively, and W i represents the wavelet matrix of the ith channel, K represents the number of seismic data channels; y and x represent the multi-channel seismic data column vector and the multi-channel reflection coefficient column vector respectively, and G is the multi-channel wavelet matrix operator; S3: Based on the geophysical forward model in step S2, by acquiring multiple prior information, an inversion system constrained by mixed prior information is constructed: Where λ1, λ2 and λ3 are adjustment constraint parameters, μ is the non-convex norm adjustment factor, ||*|| p represents the p-norm, D n is the geosteering representation matrix; R is the seismic information spatial representation matrix; S4: Obtain the data set and use the low-frequency seismic data to obtain the geosteering matrix D n ,At the same time, the inversion algorithm is used to obtain the seismic information space representation operator, and the Helix mathematical transformation is used to construct the seismic information space representation matrix R; S5: inputting the pre-stack seismic data set to be processed into the inversion model, solving the target functional by using the DCA algorithm and the alternating direction multiplier method, and outputting the seismic data inversion result.

2. The geophysical inversion method according to claim 1, characterized in that: The non-convex norm adjustment factor is set to μ = 0.

5.

3. The geophysical inversion method according to claim 1, characterized in that: Using low-frequency seismic data to obtain the geosteering matrix D in S4 n , specifically including the following contents: on the basis of extracting geological structure characteristics, constructing a spatial gradient matrix to describe the changes of geological structure in time and space, and using the spatial gradient matrix to construct the geological steering regularization constraint term, Among them, the spatial gradient matrix is ​​calculated, Where L is the spatial gradient operator, dt and dx represent the gradients along the time direction and the spatial direction respectively; According to matrix analysis theory, the eigenvalues ​​and eigenvectors of the spatial gradient matrix L can be expressed as, Where m and n represent the eigenvectors of the spatial gradient matrix L, μ m and μ n represents the eigenvalues ​​corresponding to eigenvectors m and n, m1 and m2 represent the first and second elements of eigenvector m, and n1 and n2 represent the first and second elements of eigenvector n; assuming μ m ≥μ n ≥0, then m represents the vector perpendicular to the construction direction, and n represents the vector parallel to the construction direction; Using the characteristic vector n parallel to the structural direction, a geosteering matrix parallel to the geological trend is constructed. <h2 style=";text-align:left;direction:ltr">D<h2 style=";text-align:left;direction:ltr"> n <h2 style=";text-align:left;direction:ltr"> =N1D<h2 style=";text-align:left;direction:ltr"> t <h2 style=";text-align:left;direction:ltr"> +N2D<h2 style=";text-align:left;direction:ltr"> x In the formula, represents the diagonal matrix constructed by the first element n1 of the eigenvector n, J and K represent the number of time sampling points and spatial seismic traces of the input data, represents the diagonal matrix constructed from the second element n2 of the eigenvector n, represents the first-order derivative matrix in the time direction, Represents the first-order derivative matrix with respect to the spatial direction.

4. The geophysical inversion method according to claim 1, characterized in that: In S4, the inversion algorithm is used to obtain the seismic information spatial representation operator, and the Helix mathematical transformation is used to construct the seismic information spatial representation matrix R, which specifically includes the following contents: considering the spatial continuity prior information of the seismic data, the information of the adjacent seismic data channels on the left and right sides of the target seismic data, and constructing the seismic reflection regularization constraint term. Calculate the seismic information space representation operator b i,j , In the formula, the middle "0" indicates the position of the point to be predicted, "×" indicates the point not involved in the calculation, and b i,j represents the coefficient of the spatial representation operator of seismic information, i represents the time direction, and j represents the space direction; Apply the seismic signal space characterization operator to the seismic data to obtain the seismic data prediction results. In the formula, s J,K represents the original 2D seismic data, represents the predicted seismic data, J and K represent the number of time sampling points and spatial seismic channels of the input data.

5. The geophysical inversion method according to claim 4, characterized in that: During the inversion process, the adjustment direction and step size of the model parameters are continuously updated, and the inversion model is verified and optimized by minimizing the residual functional of the original seismic data and the predicted seismic data. The residual functional of minimizing the original seismic data and the predicted seismic data can be expressed as: Among them, b i,j is the spatial representation operator of the seismic signal, s J,K represents the original 2D seismic data, represents the predicted seismic data, J and K represent the number of time sampling points and the number of spatial seismic channels.

6. The geophysical inversion method according to claim 1, characterized in that: In S5, the DCA algorithm and the alternating direction multiplier method are used to solve the above target functional, which specifically includes the following contents: Firstly, the DCA algorithm is used to decompose the mixed prior inversion system. In the formula, F(x)=λ1μ||x||2 The solution of the objective functional can be obtained by iteratively solving the following two sub-problems: In the formula, u is the intermediate variable that relates the two sub-problems; The first subproblem can be solved by the gradient method, The second subproblem can be solved by the alternating direction multiplication method. By introducing the intermediate variable z, the second subproblem can be expressed as: The Lagrange augmented expression of the above formula can be written as: In the formula, v is the Lagrange multiplier, τ is the Lagrange penalty factor; The alternating direction multiplier method can be used to iterate and solve the problem. The iterative formula can be expressed as:

7. The geophysical inversion method according to claim 1, characterized in that: In S1, seismic data is preprocessed, including the following contents: First, the pre-stack seismic data is subjected to high-precision denoising processing to protect the characteristic information of the effective signal as much as possible while attenuating the seismic noise. Then, the pre-stack seismic data is subjected to true amplitude recovery and trace equalization processing to enhance the energy of deep seismic signals and equalize the energy between trace sets. Next, the seismic data is subjected to static correction processing to eliminate the influence of surface factors and restore the original characteristics of the seismic signal. Subsequently, the seismic data is subjected to high-precision dynamic correction processing to minimize the influence of dynamic correction stretching and excision at far offsets. Finally, the pre-stack data is subjected to migration processing.

8. A geophysical inversion system based on multiple constraints, characterized in that: The inversion system is implemented by the multi-constrained geophysical inversion method according to any one of claims 1-7.

Images