Method and system for improving seismic inversion resolution through reflection coefficient protection

Abstract

The invention provides a method and system for improving seismic inversion resolution through reflection coefficient protection, and belongs to the technical field of oil-gas exploration. The method comprises the steps that seismic pure wave data, interval transit time, a density logging curve and target layer position data of a work area are acquired; building an earthquake model by utilizing the multi-well comprehensive wavelet extracted from the target work area and three adjacent earthquake channels of the earthquake pure wave data; inverting a reflection coefficient value by using a multi-task relevance vector machine Bayesian machine learning algorithm; establishing an impedance low-frequency model of a target layer by using a logging impedance curve of a target layer series of the target work area, and performing channel integration and low-frequency compensation on a reflection coefficient to obtain a high-resolution impedance inversion result; and thin dessert body tracking and reservoir quantitative depiction. According to the method, the signal-to-noise ratio of the inversion reflection coefficient is high, the weak reflection coefficient of the thin layer is well protected, particularly, the resolution ratio of the inversion result is far higher than that of an original earthquake, and the recognition effect on the shale oil thin dessert is good.

Description

Technical Field

[0001] The present invention relates to the technical field of oil and gas exploration, and in particular to a method and system for improving seismic inversion resolution by reflection coefficient protection. Background Art

[0002] Seismic inversion is one of the important geophysical technologies in oil exploration and development, and is a key means to predict the lateral distribution of reservoirs and fluids. Wave impedance is the most reliable and accurate parameter for direct seismic inversion. Wave impedance is the product of sound waves and density, and has important reference significance in predicting reservoir properties and characterizing lithologic trap boundaries.

[0003] At present, the most commonly used wave impedance inversion methods are least squares inversion and post-stack constrained sparse pulse inversion. Conventional least squares inversion is simple and efficient, but the noise resistance of the reflection coefficient is poor. The inverted reflection coefficient contains a lot of noise, the signal-to-noise ratio of the wave impedance inversion result is poor, the accuracy is low, the lateral tracking interpretation of the sand body is difficult and the multi-solution is large; another conventional post-stack constrained sparse pulse inversion has the advantages of being independent of the model, fast calculation speed, and high signal-to-noise ratio, but in order to suppress noise, the weak reflection coefficients of many thin layers will be suppressed as noise, and only the reflection coefficients with strong amplitude will be retained. However, the resolution of this inversion method depends on the resolution of the seismic, and the resolution is generally low. It cannot identify thin sand layers (thickness less than 10m), and it is less applicable to thin sweet spots (thickness less than 10m) of continental shale oil.

[0004] The key to the resolution of seismic inversion results lies in the degree of wavelet removal. When there is no noise in the seismic record, the traditional method can obtain accurate reflection coefficients. However, when the seismic record contains noise, the greater the degree of wavelet removal by traditional methods such as deconvolution, the worse the signal-to-noise ratio of the seismic record. In order to obtain a high signal-to-noise ratio and high-precision reflection coefficient, the high-frequency weak reflection coefficient will be considered to be filtered out together with the noise as invalid information, but this will greatly reduce the inversion accuracy and resolution, which is very unfavorable for the prediction of thin sweet spots.

[0005] The existing technology has at least the following deficiencies:

[0006] 1. The traditional post-stack constrained sparse pulse inversion method suppresses the weak reflection coefficients of many thin layers as noise in order to suppress noise, and only retains the reflection coefficients with strong amplitude. The resolution of the inversion depends on the resolution of the seismic, which is low and cannot identify the problem of thin sweet spots. Summary of the invention

[0007] In order to solve the problems existing in the prior art, the present invention provides a method and system for improving the resolution of seismic inversion by protecting weak reflection coefficients, and the method includes the following steps: step 1, obtaining seismic pure wave data, acoustic wave time difference and density logging curves and target layer layer data of the target work area; step 2, establishing a seismic model using three adjacent seismic traces of multi-well comprehensive sub-waves and seismic pure wave data extracted from the target work area; step 3, inverting the reflection coefficient value using a multi-task correlation vector machine Bayesian machine learning algorithm; step 4, establishing an impedance low-frequency model of the target layer using the logging impedance curve of the target layer system in the target work area, performing trace integration and low-frequency compensation on the reflection coefficient to obtain a high-resolution impedance inversion result; step 5, thin sweet spot body tracking and reservoir quantitative characterization. The reflection coefficient inverted by the present invention has a high signal-to-noise ratio, wherein the weak reflection coefficient of the thin layer is well protected, especially the resolution of the inversion result is much higher than the resolution of the original seismic, and the recognition effect of the shale oil thin sweet spot is better.

[0008] The present invention provides a method for improving seismic inversion resolution by protecting weak reflection coefficient, comprising the following steps:

[0009] Step 1: Obtain seismic pure wave data, acoustic wave time difference and density logging curves and target layer position data of the target work area;

[0010] Step 2: Establish a seismic model using three adjacent seismic traces of multi-well integrated wavelets and seismic pure wave data extracted from the target work area;

[0011] Step 3: Invert the reflection coefficient value using the multi-task relevance vector machine Bayesian machine learning algorithm;

[0012] Step 4: Use the logging impedance curve of the target layer in the target work area to establish the impedance low-frequency model of the target layer, perform trace integration and low-frequency compensation on the reflection coefficient to obtain high-resolution impedance inversion results;

[0013] Step 5: Thin sweet spot object tracking and reservoir quantitative characterization. The thin sweet spot has a thickness of less than 10m.

[0014] Preferably, in step 1, based on the seismic data collected from the target work area, a seismic imaging method is used for stacking to obtain a seismic pure wave data volume that has not been subjected to filtering and trace equalization seismic processing.

[0015] Preferably, in step one, the acoustic time difference logging curve is obtained by logging the well logging project in the target work area, the velocity logging curve is obtained by taking the inverse of the acoustic time difference logging curve, and the impedance curve is obtained by multiplying the velocity logging curve with the density logging curve.

[0016] Preferably, in step 1, the seismic record is synthesized by Ricker wavelet and impedance curve, and the comprehensive wavelet of the seismic data volume is obtained through well seismic calibration and wavelet extraction.

[0017] Preferably, the Ricker wavelet and the impedance curve are convolved to synthesize the seismic record.

[0018] Preferably, in step 1, the seismic reflection event axis corresponding to the target layer is determined, and the seismic reflection event axis tracking interpretation of the target layer is performed to obtain the time domain stratigraphic data of the target layer.

[0019] Preferably, the seismic model established in step 2 using three adjacent seismic traces of multi-well integrated wavelets and seismic pure wave data extracted from the target work area is specifically:

[0020] S=WR+e

[0021] in,

[0022] S is the seismic record vector of the three adjacent tracks of the ith track, S = {S i ,S i+1 ,S i+2};

[0023] W is the seismic sub-wave vector of the three adjacent traces of the ith trace, W = {W i ,W i+1 ,W i+2};

[0024] R is the reflection coefficient vector of the three adjacent tracks of the ith track, R = {R i ,R i+1 ,R i+2};

[0025] e is the cell array vector of the interference noise of the three adjacent channels of the ith channel, e = {e i ,e i+1 ,e i+2}.

[0026] The correlation between three adjacent seismic traces is high, and the three reflection coefficients inverted from the three adjacent seismic traces are similar. The three reflection coefficients are superimposed as the reflection coefficient of the middle seismic trace, and the reflection coefficients of all seismic traces can be obtained by analogy. Compared with the single-channel inversion reflection coefficient, the present invention can improve the signal-to-noise ratio of the inverted reflection coefficient and reduce the influence of noise on the reflection coefficient, thereby achieving the purpose of protecting the reflection coefficient.

[0027] Preferably, in step 2, the seismic model is established using three adjacent seismic traces of multi-well integrated wavelets and seismic pure wave data extracted from the target work area. The specific method is as follows:

[0028] Based on the convolution theory, seismic records can be viewed as the convolution of the wavelet w and the reflection coefficient r:

[0029] s=w*r (1.1)

[0030] Where s is the synthetic seismic record.

[0031] Write equation (1.1) in matrix form:

[0032] S i =W i R i +e i (1.2)

[0033] Among them, S i , W i , R i and e i represent the i-th seismic record trace, seismic wavelet, reflection coefficient and interference noise matrix respectively.

[0034] According to formula (1.2), the reflection coefficient inversion can be transformed into a linear regression problem under Bayesian theory. Through sparse Bayesian learning based on relevant vector machine, a priori model of unknown parameters is constructed under the method of probability theory, and sampling information is added to modify the unknown parameter model to obtain the parameter estimation with maximum probability.

[0035] However, in the actual reflection coefficient inversion problem, the sampling data of only a single channel is often incomplete, resulting in poor accuracy of the high-frequency components of the inverted reflection coefficient. Taking into account the correlation between adjacent seismic channel signals, multi-task Bayesian learning is to establish a model with strong generalization ability by integrating the existing connections between incomplete data, so as to improve the learning ability of the model and the high-precision estimation of parameters.

[0036] The earthquake model established by three adjacent seismic traces is as follows:

[0037] S=WR+e (1.3)

[0038] Preferably, step three of inverting the reflection coefficient value using a multi-task relevance vector machine Bayesian machine learning algorithm specifically includes:

[0039] Taking the technology of a single seismic trace as a task, the Gaussian likelihood model of multi-task reflection coefficient inversion in the Bayesian framework is:

[0040]

[0041] Among them, 1 / α 0 is the covariance of the cell array vector e of the interference noise;

[0042] Assume that the reflection coefficient vector is affected by the hyperparameter α vector, then R = {R i ,R i+1 ,R i+2 The prior of} follows the following zero-mean Gaussian distribution:

[0043]

[0044] in,

[0045] j is the task number in the multi-task, ranging from 1 to 3;

[0046] R i，j represents the reflection coefficient of the i-th seismic trace in the j-th task;

[0047] α and α 0 They obey Gamma distribution respectively;

[0048] α j is the element of the hyperparameter α vector for the jth task;

[0049] The calculation of a single seismic channel is one task, while the present invention calculates three adjacent channels simultaneously, thus forming three tasks, and utilizes the correlation between adjacent channels to improve the calculation accuracy.

[0050] The hyperparameters α and α are obtained according to the following formula 0 :

[0051]

[0052] in,

[0053]

[0054]

[0055]

[0056] represents the kth element on the diagonal, C i,-j It is C i The coefficient matrix after removing the elements on the diagonal;

[0057] I is the identity matrix;

[0058] A=diag(α 1 ,α 2 ,α 3 ) is α 1 ,α 2 and α 3 The diagonal matrix composed of;

[0059] According to the Bayesian formula, the posterior distribution of the reflection coefficient R is:

[0060]

[0061] in,

[0062] μ=α 0 ΣW T S (2.6)

[0063] Σ=(α 0 W T W+A) -1 (2.7)

[0064] Formulas (2.4) and (2.5) are used to iteratively solve the problem and finally the reflection coefficient R that meets the convergence condition is obtained.

[0065] Preferably, step three of inverting the reflection coefficient value using a multi-task relevance vector machine Bayesian machine learning algorithm specifically includes:

[0066] Taking the technology of a single seismic trace as a task, the Gaussian likelihood model of multi-task reflection coefficient inversion in the Bayesian framework is:

[0067]

[0068] Among them, 1 / α 0 is the covariance of the cell array vector e of the interference noise;

[0069] Assume that the reflection coefficient vector is affected by the hyperparameter α vector, then R = {R i ,R i+1 ,R i+2 The prior of} follows the following zero-mean Gaussian distribution:

[0070]

[0071] in,

[0072] j is the task number in the multi-task, ranging from 1 to 3;

[0073] R i，j represents the reflection coefficient of the i-th seismic trace in the j-th task;

[0074] α j is the element of the hyperparameter α vector for the jth task;

[0075] α and α 0 They obey the following Gamma distribution respectively:

[0076]

[0077] P(α 0 |c,d)=Γ(α 0 |c,d) (1.7)

[0078] According to the Bayesian formula, we get the hyperparameters α and α 0 The posterior distribution of is:

[0079]

[0080] When the hyperparameters α and α 0 The posterior probability is the largest, which is the hyperparameters α and α 0 The optimal estimate of , considering the complexity of formula (1.8), is obtained by taking the limit a, b, c, d → 0 and taking the logarithm of both sides of formula (1.8):

[0081]

[0082] in, A=diag(α 1 ,α 2 ,α 3 ) is α 1 ,α 2 and α 3 A diagonal matrix composed of .

[0083] Considering that the hyperparameters are independent of each other, C i Rewritten as:

[0084]

[0085] in, represents the kth element on the diagonal, C i,-j It is C i The coefficient matrix with the diagonal elements removed.

[0086] According to formula (2.0), we get:

[0087]

[0088] Substituting formula (2.1) into formula (1.9), we get:

[0089]

[0090] in,

[0091] Using formula (2.2) j Taking the partial derivative and setting it equal to 0, we get:

[0092]

[0093] Assume α j ≈s i,j , we get α j The approximate value of is as follows:

[0094]

[0095] The hyperparameters α and α can be obtained by formula (2.4) 0 , when the hyperparameters are known, the posterior distribution of the reflection coefficient R is obtained according to the Bayesian formula:

[0096]

[0097] μ=α 0 ΣW T S (2.6)

[0098] Σ=(α 0 W T W+A) -1 (2.7)

[0099] Where A = diag(α 1 ,α 2 ,α 3 ).

[0100] Formulas (2.4), (2.6) and (2.7) are used to iteratively solve the problem and finally the reflection coefficient R that meets the convergence condition is obtained.

[0101] Preferably, in step 4, the impedance curve and target layer data are calculated based on the acoustic time difference curve and density curve obtained in step 1, and the impedance model of the target layer in the target work area is established using a local weighted interpolation algorithm.

[0102] Preferably, the impedance curve is calculated by using the acoustic wave time difference curve and the density curve obtained by logging, impedance=acoustic wave velocity×density, and the acoustic wave velocity is the inverse of the acoustic wave time difference.

[0103] The stratigraphic data of the target layer are data interpreted manually or automatically by a machine. Generally, after well-seismic calibration, the seismic isochronous axis corresponding to the target layer is determined and then continuously tracked. This is a basic technology in geophysical prospecting and does not require detailed introduction.

[0104] Based on the layer data and impedance curve, spatial interpolation can be performed to obtain a three-dimensional impedance data body, which is the impedance model. The interpolation method is the local weighted method.

[0105] Preferably, in step 4, the reflection coefficient inversion result is integrated by using the established impedance model below 10 Hz as a low-frequency model to obtain the impedance inversion result after low-frequency compensation.

[0106] Preferably, in step five, based on the characteristics of different wave impedance value ranges corresponding to different lithologies in logging intersection analysis, the thin sweet spots in the target work area are predicted and laterally tracked and characterized on the inverted wave impedance data volume to obtain the distribution range of the sweet spots.

[0107] The inversion results in a three-dimensional wave impedance data body. Through the previous logging intersection analysis, it is determined that different lithologies correspond to different wave impedance value ranges. According to the wave impedance value range, the data corresponding to this range can be found on the wave impedance three-dimensional data body. The wave impedance data in this value range represents the sweet spot. This process is the prediction of the sweet spot. The lateral tracking characterization is interpreted based on the prediction results, and the distribution range of the sweet spot is obtained by extracting plane attributes or three-dimensional carving.

[0108] The present invention provides a system for improving seismic inversion resolution by using weak reflection coefficient protection, comprising a processor, wherein the processor executes the steps of any of the above-mentioned methods for improving seismic inversion resolution by using weak reflection coefficient protection.

[0109] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0110] (1) The present invention adopts a method of simultaneously obtaining the reflection coefficients of three adjacent seismic traces, and then superimposing the three reflection coefficients as the reflection coefficient of the middle seismic trace. Compared with the single-channel calculation method of the traditional sparse pulse inversion method, the calculation is more stable, and the signal-to-noise ratio of the inverted reflection coefficient is high. The weak reflection coefficient of the thin layer is better protected. In particular, the resolution of the inversion result is much higher than the resolution of the original seismic trace, and the identification effect of the thin sweet spot of shale oil is better.

[0111] (2) The present invention uses a multi-task relevance vector machine Bayesian machine learning algorithm to solve the inversion equation and obtain a high signal-to-noise ratio and high-precision reflection coefficient series, effectively protecting the weak reflection coefficient of the thin layer in the reflection coefficient series, improving the resolution of the inversion, and solving the shortcomings of the traditional inversion method in thin reservoir identification, providing more accurate impedance parameters for thin sweet spot prediction. BRIEF DESCRIPTION OF THE DRAWINGS

[0112] Figure 1 This is a flow chart of a method for improving seismic inversion resolution by weak reflection coefficient protection according to an embodiment of the present invention.

[0113] Figure 2 In the figure, (a) is Gaussian random noise at a signal-to-noise ratio of 50 dB according to an embodiment of the present invention, and (b) is a synthetic seismic record containing Gaussian random noise at a signal-to-noise ratio of 50 dB according to an embodiment of the present invention.

[0114] Figure 3 In the figure, (a) is the original reflection coefficient of an embodiment of the present invention, (b) is the reflection coefficient inverted by the method of the present invention, and (c) is the reflection coefficient inverted by the traditional least squares method.

[0115] Figure 4 In the figure, (a) is a post-stack constrained sparse pulse impedance inversion profile of an embodiment of the present invention, and (b) is an impedance inversion profile of an embodiment of the present invention. DETAILED DESCRIPTION

[0116] The following is combined with Figure 1-4 , the specific implementation methods of the present invention are described in detail.

[0117] The present invention provides a method for improving seismic inversion resolution by protecting weak reflection coefficient, comprising the following steps:

[0118] Step 1: Obtain seismic pure wave data, acoustic wave time difference and density logging curves and target layer position data of the target work area;

[0119] Step 2: Establish a seismic model using three adjacent seismic traces of multi-well integrated wavelets and seismic pure wave data extracted from the target work area;

[0120] Step 3: Invert the reflection coefficient value using the multi-task relevance vector machine Bayesian machine learning algorithm;

[0121] Step 4: Use the logging impedance curve of the target layer in the target work area to establish the impedance low-frequency model of the target layer, perform trace integration and low-frequency compensation on the reflection coefficient to obtain high-resolution impedance inversion results;

[0122] Step 5: Thin sweet spot object tracking and reservoir quantitative characterization. The thin sweet spot has a thickness of less than 10m.

[0123] According to a specific implementation of the present invention, in step 1, based on the seismic data collected from the target work area, a seismic pure wave data volume that has not been subjected to filtering and trace equalization seismic processing is obtained by stacking using a seismic imaging method.

[0124] According to a specific implementation scheme of the present invention, in step one, an acoustic time difference logging curve is obtained by logging the well logging project in the target work area, the inverse of the acoustic time difference logging curve is obtained to obtain a velocity logging curve, and the velocity logging curve is multiplied by the density logging curve to obtain an impedance curve.

[0125] According to a specific embodiment of the present invention, in step 1, the seismic record is synthesized by Ricker wavelet and impedance curve, and the comprehensive wavelet of the seismic data volume is obtained through well seismic calibration and wavelet extraction.

[0126] According to a specific embodiment of the present invention, the Ricker wavelet and the impedance curve are convolved to synthesize the seismic record.

[0127] According to a specific implementation of the present invention, in step 1, the seismic reflection event axis corresponding to the target layer is determined, and the seismic reflection event axis tracking and interpretation of the target layer is performed to obtain the time domain stratigraphic data of the target layer.

[0128] According to a specific embodiment of the present invention, the seismic model established in step 2 using three adjacent seismic traces of multi-well integrated wavelets and seismic pure wave data extracted from the target work area is specifically:

[0129] S=WR+e

[0130] in,

[0131] S is the seismic record vector of the three adjacent tracks of the ith track, S = {S i ,S i+1 ,S i+2};

[0132] W is the seismic sub-wave vector of the three adjacent traces of the ith trace, W = {W i ,W i+1 ,W i+2};

[0133] R is the reflection coefficient vector of the three adjacent tracks of the ith track, R = {R i ,R i+1 ,R i+2};

[0134] e is the cell array vector of the interference noise of the three adjacent channels of the ith channel, e = {e i ,e i+1 ,e i+2}.

[0135] The correlation between three adjacent seismic traces is high, and the three reflection coefficients inverted from the three adjacent seismic traces are similar. The three reflection coefficients are superimposed as the reflection coefficient of the middle seismic trace, and the reflection coefficients of all seismic traces can be obtained by analogy. Compared with the single-channel inversion reflection coefficient, the present invention can improve the signal-to-noise ratio of the inverted reflection coefficient and reduce the influence of noise on the reflection coefficient, thereby achieving the purpose of protecting the reflection coefficient.

[0136] According to a specific implementation of the present invention, in step 2, a seismic model is established using three adjacent seismic traces of multi-well integrated wavelets and seismic pure wave data extracted from the target work area. The specific method is as follows:

[0137] Based on the convolution theory, seismic records can be viewed as the convolution of the wavelet w and the reflection coefficient r:

[0138] s=w*r (1.1)

[0139] Where s is the synthetic seismic record.

[0140] Write equation (1.1) in matrix form:

[0141] S i =W i R i +e i(1.2)

[0142] Among them, S i , W i , R i and e i represent the i-th seismic record trace, seismic wavelet, reflection coefficient and interference noise matrix respectively.

[0143] According to formula (1.2), the reflection coefficient inversion can be transformed into a linear regression problem under Bayesian theory. Through sparse Bayesian learning based on relevant vector machine, a priori model of unknown parameters is constructed under the method of probability theory, and sampling information is added to modify the unknown parameter model to obtain the parameter estimation with maximum probability.

[0144] However, in the actual reflection coefficient inversion problem, the sampling data of only a single channel is often incomplete, resulting in poor accuracy of the high-frequency components of the inverted reflection coefficient. Taking into account the correlation between adjacent seismic channel signals, multi-task Bayesian learning is to establish a model with strong generalization ability by integrating the existing connections between incomplete data, so as to improve the learning ability of the model and the high-precision estimation of parameters.

[0145] The earthquake model established by three adjacent seismic traces is as follows:

[0146] S=WR+e (1.3)

[0147] According to a specific embodiment of the present invention, step three uses the multi-task relevance vector machine Bayesian machine learning algorithm to invert the reflection coefficient value, which specifically includes:

[0148] Taking the technology of a single seismic trace as a task, the Gaussian likelihood model of multi-task reflection coefficient inversion in the Bayesian framework is:

[0149]

[0150] Among them, 1 / α 0 is the covariance of the cell array vector e of the interference noise;

[0151] Assume that the reflection coefficient vector is affected by the hyperparameter α vector, then R = {R i ,R i+1 ,R i+2 The prior of} follows the following zero-mean Gaussian distribution:

[0152]

[0153] in,

[0154] j is the task number in the multi-task, ranging from 1 to 3;

[0155] R i，jrepresents the reflection coefficient of the i-th seismic trace in the j-th task;

[0156] α and α 0 They obey Gamma distribution respectively;

[0157] α j is the element of the hyperparameter α vector for the jth task;

[0158] The hyperparameters α and α are obtained according to the following formula 0 :

[0159]

[0160] in,

[0161]

[0162]

[0163]

[0164] represents the kth element on the diagonal, C i,-j It is C i The coefficient matrix after removing the elements on the diagonal;

[0165] I is the identity matrix;

[0166] A=diag(α 1 ,α 2 ,α 3 ) is α 1 ,α 2 and α 3 The diagonal matrix composed of;

[0167] According to the Bayesian formula, the posterior distribution of the reflection coefficient R is:

[0168]

[0169] in,

[0170] μ=α 0 ΣW T S (2.6)

[0171] Σ=(α 0 W T W+A) -1 (2.7)

[0172] Formulas (2.4) and (2.5) are used to iteratively solve the problem and finally the reflection coefficient R that meets the convergence condition is obtained.

[0173] According to a specific embodiment of the present invention, step three uses the multi-task relevance vector machine Bayesian machine learning algorithm to invert the reflection coefficient value, which specifically includes:

[0174] Taking the technology of a single seismic trace as a task, the Gaussian likelihood model of multi-task reflection coefficient inversion in the Bayesian framework is:

[0175]

[0176] Among them, 1 / α 0 is the covariance of the cell array vector e of the interference noise;

[0177] Assume that the reflection coefficient vector is affected by the hyperparameter α vector, then R = {R i ,R i+1 ,R i+2 The prior of} follows the following zero-mean Gaussian distribution:

[0178]

[0179] in,

[0180] j is the task number in the multi-task, ranging from 1 to 3;

[0181] R i，j represents the reflection coefficient of the i-th seismic trace in the j-th task;

[0182] α j is the element of the hyperparameter α vector for the jth task;

[0183] α and α 0 They obey the following Gamma distribution respectively:

[0184]

[0185] P(α 0 |c,d)=Γ(α 0 |c,d) (1.7)

[0186] According to the Bayesian formula, we get the hyperparameters α and α 0 The posterior distribution of is:

[0187]

[0188] When the hyperparameters α and α 0 The posterior probability is the largest, which is the hyperparameters α and α 0 The optimal estimate of , considering the complexity of formula (1.8), is obtained by taking the limit a, b, c, d → 0 and taking the logarithm of both sides of formula (1.8):

[0189]

[0190] in, A=diag(α 1 ,α 2 ,α 3 ) is α 1 ,α 2 and α 3 A diagonal matrix composed of .

[0191] Considering that the hyperparameters are independent of each other, C i Rewritten as:

[0192]

[0193] in, represents the kth element on the diagonal, C i,-j It is C i The coefficient matrix with the diagonal elements removed.

[0194] According to formula (2.0), we get:

[0195]

[0196] Substituting formula (2.1) into formula (1.9), we get:

[0197]

[0198] in,

[0199] Using formula (2.2) j Taking the partial derivative and setting it equal to 0, we get:

[0200]

[0201] Assume α j ≈s i,j , we get α j The approximate value of is as follows:

[0202]

[0203] The hyperparameters α and α can be obtained by formula (2.4) 0 , when the hyperparameters are known, the posterior distribution of the reflection coefficient R is obtained according to the Bayesian formula:

[0204]

[0205] μ=α 0 ΣW T S (2.6)

[0206] Σ=(α0 W T W+A) -1 (2.7)

[0207] Where A = diag(α 1 ,α 2 ,α 3 ).

[0208] Formulas (2.4), (2.6) and (2.7) are used to iteratively solve the problem and finally the reflection coefficient R that meets the convergence condition is obtained.

[0209] According to a specific implementation scheme of the present invention, in step 4, the impedance curve and target layer position data are calculated based on the acoustic wave time difference curve and density curve obtained in step 1, and the impedance model of the target layer in the target work area is established using a local weighted interpolation algorithm.

[0210] According to a specific embodiment of the present invention, the impedance curve is calculated by the acoustic time difference curve and the density curve obtained by logging, impedance = acoustic wave velocity x density, and the acoustic wave velocity is the inverse of the acoustic wave time difference.

[0211] The stratigraphic data of the target layer are data interpreted manually or automatically by a machine. Generally, after well-seismic calibration, the seismic isochronous axis corresponding to the target layer is determined and then continuously tracked. This is a basic technology in geophysical prospecting and does not require detailed introduction.

[0212] Based on the layer data and impedance curve, spatial interpolation can be performed to obtain a three-dimensional impedance data body, which is the impedance model. The interpolation method is the local weighted method.

[0213] According to a specific implementation of the present invention, in step 4, the reflection coefficient inversion result is integrated using the established impedance model below 10 Hz as a low-frequency model to obtain the impedance inversion result after low-frequency compensation.

[0214] According to a specific implementation scheme of the present invention, in step five, based on the characteristics of different wave impedance value ranges corresponding to different lithologies in logging intersection analysis, the thin sweet spots in the target work area are predicted and laterally tracked and characterized on the inverted wave impedance data volume to obtain the distribution range of the sweet spots.

[0215] The inversion results in a three-dimensional wave impedance data body. Through the previous logging intersection analysis, it is determined that different lithologies correspond to different wave impedance value ranges. According to the wave impedance value range, the data corresponding to this range can be found on the wave impedance three-dimensional data body. The wave impedance data in this value range represents the sweet spot. This process is the prediction of the sweet spot. The lateral tracking characterization is interpreted based on the prediction results, and the distribution range of the sweet spot is obtained by extracting plane attributes or three-dimensional carving.

[0216] The present invention provides a system for improving seismic inversion resolution by using weak reflection coefficient protection, comprising a processor, wherein the processor executes the steps of any of the above-mentioned methods for improving seismic inversion resolution by using weak reflection coefficient protection.

[0217] Example 1

[0218] According to a specific embodiment of the present invention, Figure 1 The method for improving seismic inversion resolution by weak reflection coefficient protection of the present invention is described in detail below.

[0219] The present invention provides a method for improving seismic inversion resolution by protecting weak reflection coefficient, comprising the following steps:

[0220] Step 1: Obtain seismic pure wave data, acoustic wave time difference and density logging curves and target layer position data of the target work area;

[0221] Step 2: Establish a seismic model using three adjacent seismic traces of multi-well integrated wavelets and seismic pure wave data extracted from the target work area;

[0222] Step 3: Invert the reflection coefficient value using the multi-task relevance vector machine Bayesian machine learning algorithm;

[0223] Step 4: Use the logging impedance curve of the target layer in the target work area to establish the impedance low-frequency model of the target layer, perform trace integration and low-frequency compensation on the reflection coefficient to obtain high-resolution impedance inversion results;

[0224] Step 5: Thin sweet spot object tracking and reservoir quantitative characterization. The thin sweet spot has a thickness of less than 10m.

[0225] Example 2

[0226] According to a specific embodiment of the present invention, the method for improving seismic inversion resolution by weak reflection coefficient protection of the present invention is described in detail below.

[0227] The present invention provides a method for improving seismic inversion resolution by protecting weak reflection coefficient, comprising the following steps:

[0228] Step 1: Obtain seismic pure wave data, acoustic wave time difference and density logging curves and target layer position data of the target work area;

[0229] Step 2: Establish a seismic model using three adjacent seismic traces of multi-well integrated wavelets and seismic pure wave data extracted from the target work area;

[0230] Step 3: Invert the reflection coefficient value using the multi-task relevance vector machine Bayesian machine learning algorithm;

[0231] Step 4: Use the logging impedance curve of the target layer in the target work area to establish the impedance low-frequency model of the target layer, perform trace integration and low-frequency compensation on the reflection coefficient to obtain high-resolution impedance inversion results;

[0232] Step 5: Thin sweet spot object tracking and reservoir quantitative characterization. The thin sweet spot has a thickness of less than 10m.

[0233] Among them, in step one, based on the seismic data collected from the target work area, the seismic imaging method is used for stacking to obtain a seismic pure wave data body that has not been filtered and trace equalized seismic processing.

[0234] Among them, in step one, the acoustic time difference logging curve is obtained by logging the well logging project in the target work area, the reciprocal of the acoustic time difference logging curve is obtained to obtain the velocity logging curve, and the velocity logging curve is multiplied by the density logging curve to obtain the impedance curve.

[0235] The seismic model established in step 2 using three adjacent seismic traces of multi-well integrated wavelets and seismic pure wave data extracted from the target work area is specifically:

[0236] S=WR+e

[0237] in,

[0238] S is the seismic record vector of the three adjacent tracks of the ith track, S = {S i ,S i+1 ,S i+2};

[0239] W is the seismic sub-wave vector of the three adjacent traces of the ith trace, W = {W i ,W i+1 ,W i+2};

[0240] R is the reflection coefficient vector of the three adjacent tracks of the ith track, R = {R i ,R i+1 ,R i+2};

[0241] e is the cell array vector of the interference noise of the three adjacent channels of the ith channel, e = {e i ,e i+1 ,e i+2}.

[0242] Among them, in step five, according to the characteristics of different wave impedance value ranges corresponding to different lithologies in the logging intersection analysis, the thin sweet spots in the target work area are predicted and laterally tracked and characterized on the inverted wave impedance data body to obtain the distribution range of the sweet spots.

[0243] Example 3

[0244] According to a specific embodiment of the present invention, the method for improving seismic inversion resolution by weak reflection coefficient protection of the present invention is described in detail below.

[0245] The present invention provides a method for improving seismic inversion resolution by protecting weak reflection coefficient, comprising the following steps:

[0246] Step 1: Obtain seismic pure wave data, acoustic wave time difference and density logging curves and target layer position data of the target work area;

[0247] Step 2: Establish a seismic model using three adjacent seismic traces of multi-well integrated wavelets and seismic pure wave data extracted from the target work area;

[0248] Step 3: Invert the reflection coefficient value using the multi-task relevance vector machine Bayesian machine learning algorithm;

[0249] Step 4: Use the logging impedance curve of the target layer in the target work area to establish the impedance low-frequency model of the target layer, perform trace integration and low-frequency compensation on the reflection coefficient to obtain high-resolution impedance inversion results;

[0250] Step 5: Thin sweet spot object tracking and reservoir quantitative characterization. The thin sweet spot has a thickness of less than 10m.

[0251] Among them, in step one, based on the seismic data collected from the target work area, the seismic imaging method is used for stacking to obtain a seismic pure wave data body that has not been filtered and trace equalized seismic processing.

[0252] Among them, in step one, the acoustic time difference logging curve is obtained by logging the well logging project in the target work area, the reciprocal of the acoustic time difference logging curve is obtained to obtain the velocity logging curve, and the velocity logging curve is multiplied by the density logging curve to obtain the impedance curve.

[0253] Among them, in step 1, the seismic record is synthesized by Ricker wavelet and impedance curve, and the comprehensive wavelet of the seismic data body is obtained through well seismic calibration and wavelet extraction.

[0254] Among them, the Ricker wavelet and the impedance curve are convolved to synthesize the seismic record.

[0255] Among them, in step 1, the seismic reflection event axis corresponding to the target layer is determined, and the seismic reflection event axis tracking and interpretation of the target layer is performed to obtain the time domain stratigraphic data of the target layer.

[0256] The seismic model established in step 2 using three adjacent seismic traces of multi-well integrated wavelets and seismic pure wave data extracted from the target work area is specifically:

[0257] S=WR+e

[0258] in,

[0259] S is the seismic record vector of the three adjacent tracks of the ith track, S = {S i ,S i+1 ,S i+2};

[0260] W is the seismic sub-wave vector of the three adjacent traces of the ith trace, W = {W i ,W i+1 ,W i+2};

[0261] R is the reflection coefficient vector of the three adjacent tracks of the ith track, R = {R i ,R i+1 ,R i+2};

[0262] e is the cell array vector of the interference noise of the three adjacent channels of the ith channel, e = {e i ,e i+1 ,e i+2}.

[0263] Among them, step three uses the multi-task relevance vector machine Bayesian machine learning algorithm to invert the reflection coefficient value, which specifically includes:

[0264] Taking the technology of a single seismic trace as a task, the Gaussian likelihood model of multi-task reflection coefficient inversion in the Bayesian framework is:

[0265]

[0266] Among them, 1 / α 0 is the covariance of the cell array vector e of the interference noise;

[0267] Assume that the reflection coefficient vector is affected by the hyperparameter α vector, then R = {R i ,R i+1 ,R i+2 The prior of} follows the following zero-mean Gaussian distribution:

[0268]

[0269] in,

[0270] j is the task number in the multi-task, ranging from 1 to 3;

[0271] R i，j represents the reflection coefficient of the i-th seismic trace in the j-th task;

[0272] α and α 0 They obey Gamma distribution respectively;

[0273] α j is the element of the hyperparameter α vector for the jth task;

[0274] The hyperparameters α and α are obtained according to the following formula 0 :

[0275]

[0276] in,

[0277]

[0278]

[0279]

[0280] represents the kth element on the diagonal, C i,-j It is C i The coefficient matrix after removing the elements on the diagonal;

[0281] I is the identity matrix;

[0282] A=diag(α 1 ,α 2 ,α 3 ) is α 1 ,α 2 and α 3 The diagonal matrix composed of;

[0283] According to the Bayesian formula, the posterior distribution of the reflection coefficient R is:

[0284]

[0285] in,

[0286] μ=α 0 ΣW T S (2.6)

[0287] Σ=(α 0 W T W+A) -1 (2.7)

[0288] Formulas (2.4) and (2.5) are used to iteratively solve the problem and finally the reflection coefficient R that meets the convergence condition is obtained.

[0289] Among them, in step 4, the impedance curve and target layer position data are calculated based on the acoustic wave time difference curve and density curve obtained in step 1, and the impedance model of the target layer in the target work area is established using the local weighted interpolation algorithm.

[0290] Among them, in step 4, the established impedance model below 10 Hz is used as a low-frequency model to perform trace integration on the reflection coefficient inversion result to obtain the impedance inversion result after low-frequency compensation.

[0291] Among them, in step five, according to the characteristics of different wave impedance value ranges corresponding to different lithologies in the logging intersection analysis, the thin sweet spots in the target work area are predicted and laterally tracked and characterized on the inverted wave impedance data body to obtain the distribution range of the sweet spots.

[0292] Example 4

[0293] According to a specific embodiment of the present invention, the method for improving seismic inversion resolution by weak reflection coefficient protection of the present invention is described in detail below.

[0294] The present invention provides a method for improving seismic inversion resolution by protecting weak reflection coefficient, comprising the following steps:

[0295] Step 1: Obtain seismic pure wave data, acoustic wave time difference and density logging curves and target layer position data of the target work area; the specific method is as follows:

[0296] 1) Based on the seismic data collected from the target work area, the seismic imaging method is used to stack and obtain the seismic pure wave data volume that has not been filtered and trace equalized seismic processing;

[0297] 2) Obtain the acoustic time difference logging curve by logging the well logging project in the target work area, obtain the velocity logging curve by taking the inverse of the acoustic time difference logging curve, and obtain the impedance curve by multiplying the velocity logging curve with the density logging curve;

[0298] 3) Synthesize seismic records by using Ricker wavelets and impedance curves, and obtain comprehensive wavelets of seismic data through well-seismic calibration and wavelet extraction;

[0299] 4) Determine the seismic reflection event axis corresponding to the target layer, perform tracking and interpretation of the seismic reflection event axis of the target layer, and obtain the time domain stratigraphic data of the target layer.

[0300] Step 2: Establish a seismic model using three adjacent seismic traces of multi-well integrated wavelets and seismic pure wave data extracted from the target work area; the specific method is as follows:

[0301] Based on the convolution theory, seismic records can be viewed as the convolution of the wavelet w and the reflection coefficient r:

[0302] s=w*r (1.1)

[0303] Among them, s is the synthetic record of earthquake.

[0304] Formula (1.1) can be written in matrix form:

[0305] S i =W i R i +e i (1.2)

[0306] Among them, Si , W i , R i and e i represent the i-th seismic record trace, seismic wavelet, reflection coefficient and interference noise matrix respectively.

[0307] According to formula (1.2), reflection coefficient inversion can be transformed into a linear regression problem under Bayesian theory. Through sparse Bayesian learning based on correlation vector machine, a priori model of unknown parameters is constructed under the method of probability theory, and sampling information is added to modify the unknown parameter model to obtain the parameter estimation with maximum probability. In the actual reflection coefficient inversion problem, the sampling data of a single channel is often incomplete, resulting in poor accuracy of high-frequency components in the inverted reflection coefficient. Considering the correlation between adjacent seismic channel signals, multi-task Bayesian learning is to establish a model with strong generalization ability by taking into account the relationship between incomplete data, so as to improve the learning ability of the model and the high-precision estimation of parameters. To this end, the present invention uses three adjacent channels of the seismic data body to establish a mathematical model for multi-task reflection coefficient inversion:

[0308] S=WR+e (1.3)

[0309] Where S = {S i ,S i+1 ,S i+2},W={W i ,W i+1 ,W i+2}, R = {R i ,R i+1 ,R i+2} and e={e i ,e i+1 ,e i+2} respectively represent the cell arrays of seismic record traces, seismic wavelets, reflection coefficients and interference noise of the three adjacent traces of the i-th trace.

[0310] Step 3: Use the multi-task relevance vector machine Bayesian machine learning algorithm to invert the reflection coefficient value; the specific method is as follows:

[0311] Under the multi-task Bayesian model, the technology of a single seismic trace is regarded as a task. The number of tasks performed by the present invention is 3. The seismic observation record vector used for each inversion of reflection coefficient is S={S i ,S i+1 ,S i+2}, the corresponding reflection coefficient vector is R = {R i ,R i+1 ,R i+2}, the purpose of multi-task learning is to establish a generalization mechanism so that each task shares all information. In the Bayesian framework, the Gaussian likelihood model for multi-task reflection coefficient inversion is:

[0312]

[0313] Among them, 1 / α 0 is the covariance of the interference noise e.

[0314] Assuming that the reflection coefficient vector is affected by the hyperparameter α vector, then R = {R i ,R i+1 ,R i+2 The prior follows a zero-mean Gaussian distribution:

[0315]

[0316] Among them, α j are the elements of the hyperparameter α vector.

[0317] Assume that α and α 0 They obey the Gamma distribution respectively, that is:

[0318]

[0319] P(α 0 |c,d)=Γ(α 0 |c,d) (1.7)

[0320] According to the Bayesian formula, we get the hyperparameters α and α 0 Posterior distribution:

[0321]

[0322] When the hyperparameters α and α 0 The posterior probability is the largest, which is the hyperparameters α and α 0 The optimal estimate of , considering the complexity of formula (1.8), is obtained by taking the limit a, b, c, d → 0 and taking the logarithm of both sides of formula (1.8):

[0323]

[0324] in, A=diag(α 1 ,α 2 ,α 3 ) is α 1 ,α 2 and α 3 A diagonal matrix composed of .

[0325] Further considering the mutual independence between hyperparameters, C i Rewritten as:

[0326]

[0327] in, represents the kth element on the diagonal, C i,-j It is C i The coefficient matrix with the diagonal elements removed.

[0328] According to formula (2.0), we get:

[0329]

[0330] Substituting formula (2.1) into formula (1.9), we get:

[0331]

[0332] in,

[0333] Using formula (2.2) j Find the partial derivative and set it equal to 0, and we get:

[0334]

[0335] The above formula is complex and difficult to solve. In order to simplify the calculation, assume that α j ≈s i,j , then the denominator of the above formula can be regarded as a constant, and we get α j Approximate value of :

[0336]

[0337] The hyperparameters α and α can be obtained by formula (2.4) 0 , when the hyperparameters are known, the posterior distribution of R is obtained according to the Bayesian formula:

[0338]

[0339] μ=α 0 ΣW T S (2.6)

[0340] Σ=(α 0 W T W+A) -1 (2.7)

[0341] Where A = diag(α 1 ,α 2 ,α 3 ).

[0342] Formulas (2.4), (2.6) and (2.7) are used to iteratively solve the problem and finally the reflection coefficient that meets the convergence conditions is obtained.

[0343] Step 4: Use the logging impedance curve of the target layer in the target work area to establish the impedance low-frequency model of the target layer, perform trace integration and low-frequency compensation on the reflection coefficient to obtain high-resolution impedance inversion results; the specific method is as follows:

[0344] The impedance curve and target layer position data are calculated based on the acoustic wave time difference curve and density curve obtained in step 1, and the impedance model of the target layer in the target work area is established using the local weighted interpolation algorithm;

[0345] The established impedance model below 10 Hz is used as the low-frequency model to perform trace integration on the reflection coefficient inversion results to obtain the impedance inversion results after low-frequency compensation.

[0346] Step 5: Thin sweet spot object tracking and reservoir quantitative characterization. The thin sweet spot has a thickness of less than 10m.

[0347] According to the different impedance range characteristics corresponding to different lithologies in logging intersection analysis, the reservoir in the work area can be predicted and laterally tracked on the inverted impedance profile to obtain the reservoir distribution range.

[0348] The seismic records are synthesized by using the P-wave velocity and density data of real well logging in the actual work area.

[0349] Figure 2 (a) is a Gaussian random noise at a signal-to-noise ratio of 50 dB according to an embodiment of the present invention, Figure 2 (b) is a synthetic seismic record containing Gaussian random noise at a signal-to-noise ratio of 50 dB according to an embodiment of the present invention. It can be seen that noise has a great influence on the weak amplitude reflection coefficient.

[0350] Figure 3 (a) is the original reflection coefficient of an embodiment of the present invention, Figure 3 (b) is the reflection coefficient inverted by the method of the present invention. It can be seen that the overall reflection coefficient is well maintained with the original reflection coefficient, and the weak reflection coefficient of the thin layer is also well protected. Figure 3 (c) is the reflection coefficient inverted by the traditional least squares method. Due to the poor noise resistance of the traditional method, the inverted reflection coefficient has been distorted.

[0351] Figure 4 (a) is a post-stack constrained sparse pulse impedance inversion profile of an embodiment of the present invention. It can be seen that the resolution of the inverted profile is low and can only distinguish sand bodies above 15 m. Figure 4 (b) is an impedance inversion profile of an embodiment of the present invention. It can be seen that the impedance inversion resolution of the method of the present invention is much higher than the post-stack constrained sparse pulse impedance inversion result, and thin sand bodies of about 8m can also be well identified.

[0352] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. For those skilled in the art, the present invention may have various modifications and variations. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present invention are included in the protection scope of the present invention.

Claims

1. A method for improving seismic inversion resolution by reflection coefficient protection, characterized in that: The steps include: Step 1: Obtain seismic pure wave data, acoustic wave time difference and density logging curves and target layer position data of the target work area; Step 2: Establish a seismic model using three adjacent seismic traces of multi-well integrated wavelets and seismic pure wave data extracted from the target work area; Step 3: Invert the reflection coefficient value using the multi-task relevance vector machine Bayesian machine learning algorithm; Step 4: Use the logging impedance curve of the target layer in the target work area to establish the impedance low-frequency model of the target layer, perform trace integration and low-frequency compensation on the reflection coefficient to obtain high-resolution impedance inversion results; Step 5: Thin sweet spot object tracking and reservoir quantitative characterization. Thin sweet spot refers to a sweet spot with a thickness of less than 10m.

2. The method for improving seismic inversion resolution by reflection coefficient protection according to claim 1, characterized in that: In step 1, based on the seismic data collected from the target work area, the seismic imaging method is used for stacking to obtain a seismic pure wave data volume that has not been filtered or subjected to trace equalization seismic processing.

3. The method for improving seismic inversion resolution by reflection coefficient protection according to claim 1, characterized in that: In step 1, the acoustic time difference logging curve is obtained by logging the well logging project in the target work area, the velocity logging curve is obtained by taking the inverse of the acoustic time difference logging curve, and the impedance curve is obtained by multiplying the velocity logging curve with the density logging curve.

4. The method for improving seismic inversion resolution by reflection coefficient protection according to claim 3, characterized in that: In step 1, the seismic record is synthesized by Ricker wavelet and impedance curve, and the comprehensive wavelet of the seismic data volume is obtained through well seismic calibration and wavelet extraction.

5. The method for improving seismic inversion resolution by reflection coefficient protection according to claim 4, characterized in that: The seismic records are synthesized by convolution of Ricker wavelet and impedance curve.

6. The method for improving seismic inversion resolution by reflection coefficient protection according to claim 4, characterized in that: In step 1, the seismic reflection event axis corresponding to the target layer is determined, and the seismic reflection event axis tracking and interpretation of the target layer is performed to obtain the time domain stratigraphic data of the target layer.

7. The method for improving seismic inversion resolution by reflection coefficient protection according to claim 1, characterized in that: The seismic model established in step 2 using three adjacent seismic traces of multi-well integrated wavelets and seismic pure wave data extracted from the target work area is as follows: S=WR+e in, S is the seismic record vector of the three adjacent tracks of the ith track, S = {S i ,S i+1 ,S i+2 }; W is the seismic sub-wave vector of the three adjacent traces of the ith trace, W = {W i ,W i+1 ,W i+2 }; R is the reflection coefficient vector of the three adjacent tracks of the ith track, R = {R i ,R i+1 ,R i+2 }; e is the cell array vector of the interference noise of the three adjacent channels of the ith channel, e = {e i ,e i+1 ,e i+2 }.

8. The method for improving seismic inversion resolution by reflection coefficient protection according to claim 7, characterized in that: Step 3 uses the multi-task relevance vector machine Bayesian machine learning algorithm to invert the reflection coefficient value, which specifically includes: Taking the technology of a single seismic trace as a task, the Gaussian likelihood model of multi-task reflection coefficient inversion in the Bayesian framework is: Among them, 1 / α0 is the covariance of the cell array vector e of the interference noise; Assume that the reflection coefficient vector is affected by the hyperparameter α vector, then R = {R i ,R i+1 ,R i+2 The prior for} follows the following zero-mean Gaussian distribution: in, j is the task number in the multi-task, ranging from 1 to 3; R i，j represents the reflection coefficient of the i-th seismic trace in the j-th task; α and α0 follow Gamma distribution respectively; α j is the element of the hyperparameter α vector for the jth task; The hyperparameters α and α0 are obtained according to the following formula: in, represents the kth element on the diagonal, C i,-j It is C i The coefficient matrix after removing the elements on the diagonal; I is the identity matrix; A=diag(α1,α2,α3) is a diagonal matrix consisting of α1,α2 and α3; According to the Bayesian formula, the posterior distribution of the reflection coefficient R is: in, μ=α0ΣW T S (2.6) Σ=(α0W T W+A) -1 (2.7) Formulas (2.4) and (2.5) are used to iteratively solve the problem and finally the reflection coefficient R that meets the convergence condition is obtained.

9. The method for improving seismic inversion resolution by reflection coefficient protection according to any one of claims 1 to 8, characterized in that: In step 4, the impedance curve and target layer data are calculated based on the acoustic time difference curve and density curve obtained in step 1, and the impedance model of the target layer in the target work area is established using the local weighted interpolation algorithm.

10. The method for improving seismic inversion resolution by reflection coefficient protection according to claim 9, characterized in that: In step 4, the established impedance model below 10 Hz is used as a low-frequency model to perform trace integration on the reflection coefficient inversion result to obtain the impedance inversion result after low-frequency compensation.

11. The method for improving seismic inversion resolution by reflection coefficient protection according to claim 1, characterized in that: In step five, based on the characteristics of different wave impedance value ranges corresponding to different lithologies in logging intersection analysis, the thin sweet spots in the target work area are predicted and laterally tracked and characterized on the inverted wave impedance data volume to obtain the distribution range of the sweet spots.

12. A system for improving seismic inversion resolution by reflection coefficient protection, characterized in that: The method comprises a processor, wherein the processor executes the steps of the method for improving seismic inversion resolution by protecting reflection coefficient as described in any one of claims 1 to 11.

Images