Bilateral filtering prestack AVO inversion method based on precise Zoeppritz equation

Abstract

The invention discloses a bilateral filtering prestack AVO inversion method based on an accurate Zoeppritz equation, which belongs to a seismic exploration method, and is characterized in that CDP gather data, logging data and target horizon data of a target area are acquired, and the accurate Zoeppritz equation is used for inversion, so that the elastic difference of a bottom layer can be increased, the inversion precision can be improved, and the accuracy of the inversion can be improved. And the addition of the second-order difference L2 norm constraint can effectively suppress noise and reduce inversion disturbance, so that the inversion result is more in line with the real situation, and meanwhile, Bilate Filter filtering is added in each iteration to protect the stratum boundary, so that the inversion resolution and precision can be further improved.

Description

Technical Field

[0001] The present invention relates to a seismic exploration method, and more specifically, mainly relates to a bilateral filtering pre-stack AVO inversion method based on an accurate Zoeppritz equation. Background Art

[0002] Seismic exploration methods are an important means for humans to obtain underground spatial information. Seismic inversion solves the spatial structure and physical properties of underground rock formations and converts seismic data into a form that can be directly compared with well logging data. According to the seismic data used for inversion, it can be divided into pre-stack inversion and post-stack inversion. Pre-stack seismic data has more and more accurate information about underground media. Therefore, pre-stack inversion technology using amplitude variation with offset (AVO) or incident angle variation (AVA) has played an increasingly important role in actual production. Traditional pre-stack inversion uses approximate formulas for calculation. The approximate formula assumes that the elastic difference between strata is small and the offset is not too large, which limits its accuracy. Therefore, it is necessary to further study and improve the aforementioned pre-stack AVO inversion method. Summary of the invention

[0003] One of the purposes of the present invention is to address the above-mentioned shortcomings and provide a bilateral filtering pre-stack AVO inversion method based on the precise Zoeppritz equation, in order to solve the technical problems that similar inversion methods in the prior art use approximate formulas for calculation, affecting the inversion accuracy.

[0004] In order to solve the above technical problems, the present invention adopts the following technical solutions:

[0005] The present invention provides a bilateral filtering pre-stack AVO inversion method based on the exact Zoeppritz equation, the method comprising the following steps:

[0006] Step A: obtaining observation data of the target area, wherein the observation data includes pre-stack incident gather data, well logging data and target layer data.

[0007] Step B: Establish the residual L2 norm between the pre-stack incident angle gather data and the simulated data, and the objective function with the L2 norm as the constraint in the second-order difference model:

[0008]

[0009] The above simulation data is obtained by forward modeling of the current model, and in the above formula, S(m) is the above objective function, d obs is the above pre-stack incident angle gather data, d cal is the above simulation data, m is the model parameters P-wave velocity, S-wave velocity and density, B is the second-order difference matrix, and λ is the regularization parameter.

[0010] Step C: Determine the initial model m0 through the above logging data and target horizon data, extract the seismic wavelet of the above pre-stack incident angle gather data through deterministic wavelet extraction technology, and use the precise Zoeppritz equation to forward model the simulated data under the current initial model m0:

[0011] R pp =Z(m0)1

[0012]

[0013] Where Z(m0)1 is the longitudinal wave reflection coefficient calculated by the exact Zoeppritz equation, R pp is the longitudinal wave reflection coefficient, W is the above-mentioned seismic wavelet, is the convolution symbol, d cal For the above simulation data.

[0014] Step D: Using the L-curve method, given the value of the range quantity λ, calculate the appropriate regularization parameter λ to balance the weight between the residual term of the objective function and the regularization term.

[0015] Step E: Apply the LM algorithm to solve the above objective function, calculate the model update Δm0, and calculate the objective function gradient g(m0):

[0016] g(m0)=G T d+λB T Bm0

[0017] Where d is the residual of the pre-stack incident angle gather data and the simulated data, G is the Jacobian matrix, B is the second-order difference matrix, m0 is the initial model, λ is the regularization parameter, and T is the matrix transpose.

[0018] Step F: Derived from the above model update Δm0 through the approximate matrix H of the Hessian matrix:

[0019] H≈G T G+λB T B

[0020] Δm0=-(H+ηI) -1 g(m0)

[0021] Where G is the Jacobian matrix, B is the second-order difference matrix, λ is the regularization parameter, η is the damping factor, I is the identity matrix, and T is the matrix transpose.

[0022] Step G: Update the initial model m0 by the above model update amount Δm0 to obtain the updated model m1.

[0023] m1=m0+Δm0

[0024] Step H: Determine whether the current solution to the objective function has reached the maximum number of iterations.

[0025] If the judgment result is yes, the inversion results of the P-wave velocity, S-wave velocity and density of the target area section obtained by solving the objective function are output, that is, the above m1 of the current iteration is output.

[0026] If the judgment result is no, the updated model m1 is filtered using the Bilateral Filter according to the following formula to protect the stratum edge, and then steps CG are repeated, so as to obtain the m1 after stratum edge protection.

[0027]

[0028] In the formula, p is the center point of the model value to be calculated, q is the input model value, and W p is the filtering parameter, S is the convolution kernel coverage area (for example, S = 5 means that the center point 5*5 range is filtered), m(p) is the model value obtained by filtering, m(p) is the original model value of the center point, m(q) is the model value of a point in the area, G C and G S are the distance weight and similarity weight respectively, t is the time corresponding to the current calculation sample point, i is the time coordinate of point q, σ C and σ S are two standard deviations respectively.

[0029] As a preferred embodiment, a further technical solution is: in the above step D, a suitable regularization parameter λ is calculated by the following formula using the L curve method:

[0030]

[0031] Where K(λ) is the curvature of the point on the L curve, η is the initial model constraint norm, ρ is the residual L2 norm between the pre-stack incident angle gather data and the simulated data, ρ′ and ρ″ represent the first and second order derivatives of ρ with respect to λ, η′ and η″ represent the first and second order derivatives of η with respect to λ, and d obs is the above pre-stack incident angle gather data, d cal is the above simulation data, B is the second-order difference matrix, and m is the model parameters P-wave velocity, S-wave velocity and density.

[0032] A further technical solution is that the above maximum number of iterations is an arbitrarily set value.

[0033] A further technical solution is: the above inversion results are compared with the logging data curve, and appropriate parameters are adjusted, and then each pre-stack incidence angle gather data in the above observation data is inverted to obtain the inversion results of the longitudinal wave velocity, shear wave velocity and density of the entire observation profile of the target area.

[0034] Compared with the prior art, one of the beneficial effects of the present invention is that the use of precise Zoeppritz equation inversion can increase the bottom elastic difference, which helps to improve the inversion accuracy, and the addition of the second-order difference L2 norm constraint can effectively suppress noise, reduce inversion disturbance, and make the inversion result more consistent with the actual situation. At the same time, the addition of Bilateral Filter filtering in each iteration to protect the formation boundary can further improve the inversion resolution and accuracy. BRIEF DESCRIPTION OF THE DRAWINGS

[0035] Figure 1 The present invention is a flowchart for describing an embodiment of the present invention. DETAILED DESCRIPTION

[0036] The present invention will be further described below in conjunction with the accompanying drawings.

[0037] refer to Figure 1 As shown, one embodiment of the present invention is a bilateral filtering pre-stack AVO inversion method based on the exact Zoeppritz equation, which is operated according to the following steps:

[0038] Step 1: Obtain observation data of the target area, which includes pre-stack incident angle gather data, well logging data and target layer data.

[0039] Step 2: Establish the residual L2 norm between the pre-stack incident angle gather data and the simulated data, and the objective function with the L2 norm as the constraint in the second-order difference model:

[0040]

[0041] In the formula of this step, S(m) is the above objective function, d obs is the above pre-stack incident angle gather data, d cal is the above simulation data, m is the model parameters of P-wave velocity, S-wave velocity and density, B is the second-order difference matrix, and λ is the regularization parameter. And in this step, the above simulation data is obtained by forward modeling of the following current model.

[0042] Step 3: Determine the initial model m0 through the above logging data and target horizon data, and extract the seismic wavelet of the above pre-stack incident angle gather data through deterministic wavelet extraction technology, and use the precise Zoeppritz equation to forward model the simulation data under the current initial model m0:

[0043] R pp =Z(m0)1

[0044]

[0045] In the formula of this step, Z(m0)1 is the longitudinal wave reflection coefficient calculated by the exact Zoeppritz equation, R pp is the longitudinal wave reflection coefficient, W is the above-mentioned seismic wavelet, is the convolution symbol, d cal For the above simulation data.

[0046] Step 4: Using the L-curve method, given the value of the range quantity λ, calculate the appropriate regularization parameter λ to balance the weight between the residual term of the objective function and the regularization term.

[0047] In this step, preferably, a suitable regularization parameter λ is calculated by the following formula:

[0048]

[0049] In the formula of this step, K(λ) is the curvature of the point on the L curve, η is the initial model constraint norm, ρ is the residual L2 norm between the prestack incident angle gather data and the simulated data, ρ′ and ρ″ represent the first and second order derivatives of ρ with respect to λ, η′ and η″ represent the first and second order derivatives of η with respect to λ, and d obs is the above pre-stack incident angle gather data, d cal is the above simulation data, B is the second-order difference matrix, and m is the model parameters P-wave velocity, S-wave velocity and density.

[0050] Step 5: Apply the LM algorithm to solve the above objective function, calculate the model update Δm0, and calculate the objective function gradient g(m0).

[0051] g(m0)=G T d+λB T Bm0

[0052] In the formula of this step, d is the residual of the pre-stack incident angle gather data and the simulated data, G is the Jacobian matrix, B is the second-order difference matrix, m0 is the initial model, λ is the regularization parameter, and T is the matrix transpose.

[0053] Step 6: Derived the above model update Δm0 through the approximate matrix H of the Hessian matrix:

[0054] H≈G T G+λB T B

[0055] Δm0=-(H+ηI) -1g(m0)

[0056] In the formula of this step, G is the Jacobian matrix, B is the second-order difference matrix, λ is the regularization parameter, η is the damping factor, I is the identity matrix, and T is the matrix transpose.

[0057] Step 7: Update the initial model m0 by the above model update amount Δm0 to obtain the updated model m1;

[0058] m1=m0+Δm0

[0059] Step 8: Determine whether the current solution to the above objective function has reached the maximum number of iterations;

[0060] If the judgment result is yes, the inversion results of the P-wave velocity, S-wave velocity and density of the target area section obtained by solving the objective function are output, that is, the above m1 of the current iteration is output.

[0061] If the judgment result is no, the updated model m1 is filtered using the bilateral filter according to the following formula to protect the stratum edge, and then steps 3 to 7 are repeated to obtain the above m1 after stratum edge protection.

[0062]

[0063] In the formula, p is the center point of the model value to be calculated, q is the input model value, and W p is the filtering parameter, S is the convolution kernel coverage area (for example, S = 5 means filtering calculation is performed within the 5*5 range of the center point), is the model value obtained by filtering calculation, m(p) is the original model value of the center point, m(q) is the model value of a point in the area, G C and G S are the distance weight and similarity weight respectively, t is the time corresponding to the current calculation sample point, i is the time coordinate of point q, σ C and σ S are two standard deviations respectively.

[0064] In step eight, the number of iterations can be set arbitrarily, the objective function is solved iteratively, the inversion result is compared with the logging data curve, and the appropriate parameters are adjusted. Then, each pre-stack incident angle gather data in the observation data is inverted to obtain the inversion results of the longitudinal wave velocity, shear wave velocity and density of the entire observation section of the target area.

[0065] As mentioned in the background technology, the traditional pre-stack inversion uses an approximate formula for calculation. The approximate formula assumes that the elastic difference between the strata is small and the offset distance is not too large, which limits its accuracy. Therefore, in this embodiment, the use of the precise Zoeppritz equation inversion can avoid these problems. Moreover, the addition of the second-order difference L2 norm constraint can effectively suppress noise, reduce inversion disturbances, and make the inversion results more in line with the actual situation. The general pre-stack inversion method has a strong blurring effect on the stratum boundary, resulting in a decrease in resolution. Therefore, a bilateral filter is added in each iteration to protect the stratum boundary to improve the inversion resolution and accuracy.

[0066] References to "one embodiment", "another embodiment", "embodiment", etc. in this specification refer to specific features, structures or characteristics described in conjunction with the embodiment as included in at least one embodiment generally described in this application. The appearance of the same expression in multiple places in the specification does not necessarily refer to the same embodiment. Further, when a specific feature, structure or characteristic is described in conjunction with any embodiment, it is claimed that the realization of such feature, structure or characteristic in conjunction with other embodiments also falls within the scope of the present invention.

[0067] Although the present invention is described herein with reference to a number of illustrative embodiments of the present invention, it will be appreciated that those skilled in the art may devise many other modifications and implementations that fall within the scope and spirit of the principles disclosed herein. More specifically, within the scope of the present disclosure, drawings, and claims, a variety of variations and improvements may be made to the components and / or layout of the subject combination layout. In addition to variations and improvements made to the components and / or layout, other uses will also be apparent to those skilled in the art.

Claims

1. A bilateral filtering prestack AVO inversion method based on the exact Zoeppritz equation, characterized in that The method comprises the following steps: Acquiring observation data of the target area, wherein the observation data includes pre-stack incident angle gather data, well logging data and target layer data; The residual L2 norm between the pre-stack incident angle gather data and the simulated data of the seismic is established, as well as the objective function with the L2 norm as the constraint in the second-order difference model: Where S(m) is the objective function, d obs is the pre-stack incident angle gather data, d cal is the simulation data, m is the model parameters P-wave velocity, S-wave velocity and density, B is the second-order difference matrix, λ is the regularization parameter; The initial model m0 is determined by the logging data and the target layer data, and the seismic wavelet of the pre-stack incident angle gather data is extracted by the deterministic wavelet extraction technology, and the simulation data under the current initial model m0 is forward modeled by the precise Zoeppritz equation: R pp =Z(m0)1 Where Z(m0)1 is the longitudinal wave reflection coefficient calculated by the exact Zoeppritz equation, R pp is the longitudinal wave reflection coefficient, W is the seismic wavelet, is the convolution symbol, d cal is the simulation data; Using the L-curve method, given the value of the range quantity λ, the appropriate regularization parameter λ is calculated to balance the weight between the residual term of the objective function and the regularization term; The LM algorithm is applied to solve the objective function, calculate the model update Δm0, and calculate the objective function gradient g(m0): g(m0)=G T d+λB T Bm0 Where d is the residual of pre-stack incident angle gather data and simulation data, G is the Jacobian matrix, B is the second-order difference matrix, m0 is the initial model, λ is the regularization parameter, and T is the matrix transpose; The model update Δm0 is derived through the approximate matrix H of the Hessian matrix: H≈G T G+λB T B Δm0=-(H+ηI) -1 g(m0) Where G is the Jacobian matrix, B is the second-order difference matrix, λ is the regularization parameter, η is the damping factor, I is the identity matrix, and T is the matrix transpose; Update the initial model m0 by the model update amount Δm0 to obtain an updated model m1; m1=m0+Δm0 Determine whether the current solution to the objective function has reached a maximum number of iterations; If the judgment result is yes, the inversion results of the longitudinal wave velocity, shear wave velocity and density of the target area section obtained by solving the objective function are output; If the judgment result is no, the updated model m1 is filtered using Bilateral Filter to protect the stratum edge: In the formula, In the formula, p is the center point of the model value to be calculated, q is the input model value, and W p is the filtering parameter, S is the convolution kernel coverage area, is the model value obtained by filtering calculation, m(p) is the original model value of the center point, m(q) is the model value of a point in the area, G C and G S are the distance weight and similarity weight respectively, t is the time corresponding to the current calculation sample point, i is the time coordinate of point q, σ C and σ S are two standard deviations respectively.

2. The bilateral filtering prestack AVO inversion method based on the exact Zoeppritz equation according to claim 1 is characterized in that The L-curve method is used to calculate the appropriate regularization parameter λ by the following formula: Where K(λ) is the curvature of the point on the L curve, η is the initial model constraint norm, ρ is the residual L2 norm between the pre-stack incident angle gather data and the simulated data, ρ′ and ρ″ represent the first and second order derivatives of ρ with respect to λ, η′ and η″ represent the first and second order derivatives of η with respect to λ, and d obs is the pre-stack incident angle gather data, d cal is the simulation data, B is the second-order difference matrix, and m is the model parameters of longitudinal wave velocity, shear wave velocity and density.

3. The bilateral filtering prestack AVO inversion method based on the exact Zoeppritz equation according to claim 1, characterized in that: The maximum number of iterations is an arbitrarily set value.

4. The bilateral filtering prestack AVO inversion method based on the exact Zoeppritz equation according to claim 1, characterized in that: The inversion result is compared with the logging data curve, and appropriate parameters are adjusted. Then, each pre-stack incidence angle gather data in the observation data is inverted to obtain the inversion results of the longitudinal wave velocity, shear wave velocity and density of the entire observation section of the target area.

Images