Geophysical Inversion Method and System Based on Multiple Constraints

Abstract

The present invention discloses a geophysical inversion method and system based on multiple constraints, which are applied to the field of geophysical exploration. Aiming at the problem that the existing Bayesian inversion theoretical framework is easily affected by various interferences such as random noise, the present invention will simultaneously utilize multiple mixed prior information such as the characteristics of the underground geological structure and the spatial extension relationship of seismic data, and combine the compressed sensing theory. Considering the lateral continuity and extension relationship constraints of seismic data traces, it reduces the non-uniqueness of the inversion result while improving the reliability of the inversion result.

Description

Technical Field

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

[0002] Inversion is one of the important techniques in seismic data processing and has been widely used in high-resolution seismic exploration. Considering that the bandwidth of seismic data is limited, resulting in the lack of low-frequency and high-frequency components outside the effective frequency band, which leads to strong non-uniqueness in the inversion results. Therefore, various mathematical assumptions are needed, 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. Currently, a sparse model is usually used to make assumptions about the reflection coefficient sequence, and mathematical constraints are incorporated into the inversion system to reduce the non-uniqueness of the inversion results.

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

[0004] Aiming at the deficiencies in the prior art, the present invention provides a geophysical inversion method and system based on multiple constraints to alleviate the non-uniqueness and uncertainty of existing geophysical inversion results.

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

[0006] S1: Obtain the pre-stack seismic data record to be processed and preprocess the seismic data, where the pre-stack seismic data record includes multiple consecutive pre-stack seismic data records;

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

[0008] Briefly written as y = Gx,

[0009] In the formula, d i and r i respectively represent the i-th trace seismic data and the i-th trace reflection coefficient, W i represents the wavelet matrix of the i-th trace, and K represents the number of seismic data traces; y and x respectively represent the multi-trace seismic data column vector and the multi-trace reflection coefficient column vector, and G is the multi-trace wavelet matrix operator;

[0010] S3: Based on the geophysical forward model in step S2, construct a hybrid prior inversion system by obtaining multiple prior information;

[0011]

[0012] where λ 1 , λ 2 and λ 3 are adjustment constraint parameters, μ is a non-convex norm adjustment factor, ||*|| p represents the p-norm, D n is the geosteering characterization matrix; R is the seismic information space characterization matrix.

[0013] S4: Obtain seismic data records, and use low-frequency seismic data to calculate the geosteering matrix D n , and at the same time, use the inversion algorithm to calculate the seismic information space characterization operator, and use the Helix mathematical transform to construct the seismic information space characterization matrix R;

[0014] S5: Input the to-be-processed prestack seismic data records into the inversion system, and use the DCA algorithm and the alternating direction multiplier method to solve the above-mentioned objective functional, and output the seismic data inversion result.

[0015] Furthermore, in the above step S5, using the DCA algorithm and the alternating direction multiplier method to solve the above-mentioned objective functional to obtain the seismic data inversion result, the specific process is as follows:

[0016] Considering that the hybrid prior inversion system is a non-linear objective functional and cannot be directly solved by gradient-based algorithms, therefore, first use the DCA algorithm to decompose it,

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

[0018] where F(x) = λ 1 μ||x|| 2 .

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

[0020]

[0021] where u is an intermediate variable that associates the two sub-problems.

[0022] The first sub-problem is a linear inversion problem and can be solved by the gradient method,

[0023]

[0024] The second sub-problem is an L1 norm constraint problem and can be solved by the alternating direction multiplier method, and its iterative formula can be expressed as:

[0025]

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

[0027] 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.

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

[0029] In order to more clearly illustrate the specific embodiments of the present invention or the technical solutions in the prior art, the following will briefly introduce the drawings required for the description of the specific embodiments or the prior art. Obviously, the drawings in the following description are some embodiments of the present invention. For those of ordinary skill in the art, other drawings can be obtained based on these drawings without creative efforts.

[0030] Figure 1 Schematic diagram for constructing the non-convex constraint term of the reflection coefficient of the present invention;

[0031] Figure 2 Schematic diagram of the geological steering prior information of the present invention;

[0032] Figure 3 Schematic diagram of the seismic information spatial representation operator of the present invention. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0033] The following will clearly and completely describe the technical solutions in the embodiments of the present application with reference to the drawings in the embodiments of the present application. Obviously, the described embodiments are only a part of the embodiments of the present application, rather than all the embodiments. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present application without creative efforts belong to the scope of protection of the present application.

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

[0035] S1: Obtain the pre-stack seismic data records to be processed and preprocess the seismic data. Among them, the pre-stack seismic data records include multiple consecutive pre-stack seismic data records;

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

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

[0038] Abbreviated as y = Gx,

[0039] where d i and r i respectively represent the seismic data of the i-th trace and the reflection coefficient of the i-th trace, W i represents the wavelet matrix of the i-th trace, and K represents the number of traces of the seismic data; y and x respectively represent the multi-trace seismic data column vector and the multi-trace reflection coefficient column vector, and G is the multi-trace wavelet matrix operator;

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

[0041] S3: Based on the geophysical forward model in step S2, construct a hybrid prior inversion system by obtaining multiple prior information;

[0042] Conventional inversion methods generally only consider the constraint relationship of the vertical formation reflection coefficient and do not consider the influence between different seismic traces horizontally. The inversion method is easily affected by various interferences such as random noise, and the inversion result has large multi-solution and uncertainty. Therefore, in order to jointly use the low-rank prior and the inversion network prior to remove noise, the following hybrid prior optimization model is proposed:

[0043]

[0044] where λ 1 、λ 2 and λ 3To adjust the constraint parameter, μ is the non-convex norm adjustment factor, and ||*|| p represents the p-norm, and D n is the matrix representing the horizontal seismic information; R is the matrix representing the spatial seismic information; using the Helix mathematical transformation, the matrix representing the spatial seismic information can be constructed.

[0045] To ensure the stability and accuracy of the reflection coefficient inversion process in seismic data processing, it is necessary to introduce prior information constraints on the reflection coefficient. This involves using the non-convex norm regularization method in compressive sensing theory, which controls the sparsity of the reflection coefficient by imposing specific constraints during the optimization process.

[0046] Specifically, by constructing a non-convex constraint term that can effectively control the sparsity of the reflection coefficient and ensure the well-posedness of the reflection coefficient inversion process, it is necessary to introduce prior information constraints on the reflection coefficient for reference Figure 1 as shown, a non-convex constraint term for the reflection coefficient is constructed using the non-convex norm in the field of compressive sensing:

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

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

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

[0050] For reference Figure 2 as shown, considering the characteristics of the underground geological structure recorded by seismic data, such as faults, folds, lithological changes, etc. The extraction of these characteristics usually depends on the preprocessing and interpretation of seismic data, such as the analysis of attributes such as amplitude, frequency, and phase.

[0051] S4: Obtain the data record, and use the low-frequency seismic data to calculate the geological steering matrix D n , and at the same time, use the inversion algorithm to calculate the operator representing the spatial seismic information, and use the Helix mathematical transformation to construct the matrix R representing the spatial seismic information;

[0052] The data set is synthetic seismic data or actual seismic data collected. 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.

[0053] On the basis of extracting geological structure features, a spatial gradient matrix is constructed to describe the changes of geological structure in time and space. Using the spatial gradient matrix, a geological steering regularization constraint term is constructed. Using the spatial gradient matrix, a geological steering regularization constraint term is constructed.

[0054]

[0055] In the formula, D n is the geological steering matrix parallel to the tectonic direction;

[0056] Furthermore, calculate the spatial gradient matrix,

[0057]

[0058] In the formula, L is the spatial gradient operator, and dt and dx represent the gradients along the time direction and the spatial direction respectively.

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

[0060]

[0061] In the formula, m and n represent the eigenvectors of the spatial gradient matrix L, and μ m and μ n represent the eigenvalues corresponding to the eigenvectors m and n, m 1 and m 2 represent the first and second elements of the eigenvector m, n 1 and n 2 represent the first and second elements of the eigenvector n. Assume that μ m ≥μ n ≥0, then m represents the vector perpendicular to the tectonic direction, and n represents the vector parallel to the tectonic direction.

[0062] Using the eigenvector n parallel to the tectonic direction, construct a geological steering matrix parallel to the geological strike

[0063] D n =N 1 D t +N 2 D x

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

[0065] Refer to Figure 3 As shown, during the acquisition and processing of seismic data, in order to ensure the continuity and integrity of the data, it is necessary to consider the spatial extension relationship of seismic waves propagating in the subsurface medium. This relationship not only involves the propagation of seismic waves in the vertical direction but also includes the continuity in the horizontal direction. To enhance the continuity of seismic data, the information of adjacent seismic data traces on both the left and right sides of the target seismic data can be introduced, and using the seismic information spatial representation operator, a seismic reflection regularization constraint term can be constructed.

[0066]

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

[0068] Calculate the seismic information spatial representation operator b i,j ,

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

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

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

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

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

[0074] In the formula, the position of the middle "0" represents the position of the point to be predicted, "×" represents the points not participating in the calculation, b i,j represents the coefficient of the seismic signal spatial representation operator, i represents the time direction, and j represents the spatial direction.

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

[0076]

[0077] 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.

[0078] During the inversion process, the adjustment direction and step size of the model parameters are continuously updated, and the inversion system 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:

[0079]

[0080] 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.

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

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

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

[0084] 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.

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

[0086] In the formula, F(x)=λ 1 μ||x|| 2 .

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

[0088]

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

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

[0091]

[0092] The second sub - problem is an L1 - norm constraint problem, which can be solved by the alternating direction method of multipliers, and its iterative formula can be expressed as:

[0093]

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

[0095] The method of the present invention not only has stronger sparse reflection coefficient constraint conditions, but also has better spatial continuity and extension relationship constraints. It can reduce the non - uniqueness of the inversion results while improving the reliability of the inversion results, and can provide key technical support for high - resolution seismic exploration.

[0096] 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.

[0097] The present invention provides a computer storage medium, which stores at least one instruction, and the instruction is suitable for being loaded and executed by a processor to perform the method steps of one or more embodiments of this specification.

[0098] The present invention provides a computer program product, which stores at least one instruction, and the instruction is suitable for being loaded and executed by a processor to perform the method steps of one or more embodiments of this specification.

[0099] 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 and executed by the processor to perform the method steps of one or more embodiments of this specification.

[0100] The above embodiments only represent several implementation manners of the present invention, and their descriptions are relatively specific and detailed, but should not be construed as limiting the scope of the present invention. It should be noted that for those of ordinary skill in the art, without departing from the concept of the present invention, several deformations and improvements can still be made, and these all belong to the protection scope of the present invention. Therefore, the protection scope of the present invention should be subject to the appended claims.

Claims

1. A geophysical inversion method based on multiple constraints, characterized in that: The following steps are involved: S1: obtaining a pre-stack seismic data record to be processed and pre-processing the seismic data, wherein the pre-stack seismic data record includes a plurality of continuous pre-stack seismic data records; 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, a hybrid priori inversion system is constructed by obtaining multiple priori information: 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 seismic data records and use 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 record to be processed into the inversion system, 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 system 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