Changes in seismic amplitude indirectly reflect the spatial variation of underground reservoir parameters. Since there exists a certain geophysical mapping relationship between observed seismic data and reservoir parameters (Li et al., 2022), the interpolation of unknown reservoir parameters in the subsurface can be completed by guiding known well information based on known observed seismic data. The stochastic inversion method based on compressed sensing frequency-division waveform indication prior draws on the traditional geostatistical prior information construction method. The core approach for building a priori information using seismic waveform indication is as follows. The first point is to collect statistics on the well-side seismic trace data of all wells and the prediction parameter curves on the wells based on the data of all wells in the work area. The second point is to solve the low-rank equation group through the compressed sensing method to obtain the correlation coefficient between the seismic waveform at the point to be predicted and the seismic waveform at all sample points of the well-side seismic trace. The third point is the unbiased optimal estimation of well data with high-frequency components using pseudo-Kriging interpolation methods. The final inverted seismic waveforms are adjusted to match the original seismic characteristics, including the more deterministic high-frequency elastic parameters.

Waveform indication modeling is based on the geophysical mapping relationship between observed seismic data and parameters to be inverted, utilizing the similarity of seismic waveforms to simulate unknown underground model parameters from selected well samples. Zhou et al. (2021) described the modeling process of the parameters to be inverted using the ordinary pseudo-Kriging interpolation formula: Z x 0 = ∑ i = 1 n λ i Z x i (1) where, Z x 0 is the value at the interpolation parameter x 0 to be simulated. Z x i is the corresponding value at the predicted well sample point x i selected by the well-side seismic data. λ i is the weight of Z x i at the point to be interpolated. n is the number of participating wells.

Assume that the expected c and variance σ 2 for any point p x , y in a region are the same, and for any point p x , y , we have: E p x , y = E p = c Δ p x , y = σ 2

The unbiased estimation condition is: E z ^ 0 − z 0 = 0 (2)

The ordinary Kriging interpolation formula is: z ^ 0 = ∑ i = 1 n λ i z i (3) where, z 0 is the seismic data at the point to be predicted. z ^ 0 is the seismic data at the point to be predicted after Kriging interpolation. z i is the seismic data at known points. Substituting Equation 3 into Equation 2, there is: E ∑ i = 1 n λ i z i − z 0 = 0 ∑ i = 1 n λ i = 1

The estimated variance is: J = Δ z ^ 0 − z 0 = ∑ i = 1 n ∑ j = 0 n λ i λ j Cov z i , z j − 2 ∑ i = 1 n λ i Cov z i , z 0 + σ 2

Let C i j = Cov z i , z j , then we have: J = ∑ i = 1 n ∑ j = 0 n λ i λ j C i j − 2 ∑ i = 1 n λ i C i 0 + σ 2 (4)

The correlation coefficient is: r i j = C i j Δ z i Δ z j (5) where, r i j is the correlation coefficient between the seismic waveform of seismic data point z i and the seismic waveform of seismic data point z j .Substituting Equation 5 into Equation 4, we have: J = ∑ i = 1 n ∑ j = 0 n λ i λ j r i j Δ z i Δ z j − 2 ∑ i = 1 n λ i r i 0 Δ z i Δ z 0 + σ 2 = ∑ i = 1 n ∑ j = 0 n λ i λ j r i j σ 2 − 2 ∑ i = 1 n λ i r i 0 σ 2 + σ 2

Our goal is to find the set of λ i that minimizes J , and J is a function of λ i . Therefore, the partial derivative of J with respect to λ i is calculated and set to 0, that is: ∂ J ∂ λ i = 0

The optimal λ i obtained from the solution needs to satisfy the ∑ i = 1 n λ i = 1 . Therefore, the new objective function can be reconstructed: W = J + 2 ϕ ∑ i = 1 n λ i − 1

We solve the optimization problem through the above formula. Where ϕ is the Lagrange multiplier. Solving the minimum parameter set ϕ and λ 1 , λ 2 , ⋯ , λ n of W , we can get the minimized J under the ∑ i = 1 n λ i = 1 constraint, that is: ∂ J + 2 ϕ ∑ i = 1 n λ i − 1 ∂ λ k = 0   k = 1 , 2 , ⋯ , n ∂ J + 2 ϕ ∑ i = 1 n λ i − 1 ∂ ϕ = 0 ∂ ∑ i = 1 n ∑ j = 0 n λ i λ j r i j σ 2 − 2 ∑ i = 1 n λ i r i 0 σ 2 + σ 2 + 2 ϕ ∑ i = 1 n λ i − 1 ∂ λ k = 0 ∂ ∑ i = 1 n ∑ j = 0 n λ i λ j r i j σ 2 − 2 ∑ i = 1 n λ i r i 0 σ 2 + σ 2 + 2 ϕ ∑ i = 1 n λ i − 1 ∂ ϕ = 0 ⇒ ∑ j = 1 n λ j r k j + r j k σ 2 − 2 r k 0 σ 2 + 2 ϕ = 0 ∑ i = 1 n λ i = 1 (6)

Due to r k j = r j k , Then Equation 6 becomes: ∑ j = 1 n λ j r k j − r k 0 + ϕ σ 2 = 0   k = 1 , 2 , ⋯ , n ∑ i = 1 n λ i = 1 (7)

From Equation 7, we can get the equation group for solving the weight coefficient: r 11 λ 1 + r 12 λ 2 + ⋯ + r 1 n λ n + ϕ σ 2 = r 10 r 21 λ 1 + r 22 λ 2 + ⋯ + r 2 n λ n + ϕ σ 2 = r 20 ⋮ r n 1 λ 1 + r n 2 λ 2 + ⋯ + r n n λ n + ϕ σ 2 = r n 0

Written in matrix form: r 11 r 12 ⋯ r 1 n 1 r 21 r 22 ⋯ r 2 n 1 ⋮ r n 1 r n 2 ⋯ r n n 1 1 1 ⋯ 1 0 λ 1 λ 2 ⋮ λ n ϕ σ 2 = r 10 r 20 ⋮ r n 0 1 (8)

By performing an inverse operation on the matrix, the weight of the seismic data at the known point at the point to be predicted can be obtained. According to the geophysical mapping relationship between the observed seismic data and the parameters to be inverted, the weight λ j obtained from this matrix replaces λ i in Equation 1, thereby completing the establishment of the initial model.

Seismic reflections from reservoirs of varying thicknesses correspond to specific dominant frequencies in the frequency domain. Generally, high-frequency tuning energy reflects the tuning response of thin layers, while low-frequency energy corresponds to the tuning response of thick layers. The tuning response of the wavelet varies with reservoir thicknesses, as geological bodies of different thicknesses are associated with distinct dominant frequencies. Therefore, when modeling waveform indication, we first model by dividing the frequency bands and establish elastic parameter models under different frequency bands. These models in different frequency bands are integrated in the frequency domain to obtain the final required elastic parameter model.

Next, we simplify Equation 8, let A = r 11 r 12 ⋯ r 1 n 1 r 21 r 22 ⋯ r 2 n 1 ⋮ r n 1 r n 2 ⋯ r n n 1 1 1 ⋯ 1 0 , x = λ 1 λ 2 ⋮ λ n ϕ σ 2 , b = r 10 r 20 ⋮ r n 0 1 . Then Equation 8 can be written as: A x = b

However, when obtaining the matrix x containing weight coefficients, A is a matrix with insufficient rank and very small rank. It is extremely unstable when performing the inverse operation, resulting in the obtained matrix x containing the weight coefficients being inaccurate and unstable. As a result, when modeling waveform indications, the modeling results deviate significantly from the real data and the modeling results are inaccurate. To address these modeling discrepancies caused by an unreliable weight coefficient matrix, we apply the compressed sensing method to obtain the weight coefficient matrix. L1 regularized linear regression is expressed as follows (Chen et al., 2022): l x = 1 / 2 ∥ A x − b ∥ 2 2 (9)

When we solve Equation 9, we are also solving the following equation: minimize  1 / 2 ∥ A x − b ∥ 2 2 + λ ∥ x ∥ 1

At this point, we introduce the alternating direction method of multipliers to solve: minimize  f x + g z subject to  x − z = 0 (10) where, f x = 1 / 2 ∥ A x − b ∥ 2 2 , g z = λ ∥ z ∥ 1 . Based on Equation 9 to Equation 10, we can derive: x k + 1 = arg ⁡ ⁡ min x 1 / 2 ∥ A x − b ∥ 2 2 + ρ 2 x − z k + u k 2 2 z k + 1 = arg ⁡ ⁡ min z ⁡ λ ∥ z ∥ 1 + ρ 2 x k − z + u k 2 2 u k + 1 = u k + x k − z k

Derivation to zero for x update: A T A x − b + ρ x − z k + u k = 0 A T A + ρ I x − A T b + ρ z k − u k = 0 x k + 1 = A T A + ρ I − 1 A T b + ρ z k − u k

A soft thresholding is introduced to constrain the update of z : x k + 1 : = A T A + ρ I − 1 A T b + ρ z k − u k z k + 1 : = S λ / ρ x k + 1 + u k u k + 1 : = u k + x k + 1 − z k + 1

Here A T A + ρ I is always reversible and since ρ > 0 . S is the soft threshold operator with the following formula: S λ / ρ x k + 1 + u k = x k + 1 + u k − λ / ρ , x k + 1 + u k > λ / ρ 0 , x k + 1 + u k ⩽ λ / ρ x k + 1 + u k + λ / ρ , x k + 1 + u k < − λ / ρ

Therefore, a more accurate weight coefficient matrix can be obtained by using the aforementioned compressed sensing method, resulting in more reliable prior information.

To validate the effectiveness of the method in this paper, the marmousi2 partial model is selected for testing. We take near angle elastic impedance as an example for simulation. Figure 1 presents the small-angle elastic impedance model data (Figure 1A), the corresponding synthetic seismic data (Figure 1B), and the extracted pseudo-well data. The pseudo-well data include one pseudo-well (Figure 1C), three pseudo-wells (Figure 1D), and five pseudo-wells (Figure 1E). Figure 2 shows the Kriging simulation results (Figures 2A, C, E) and the simulation results (Figures 2B, D, F) of this paper’s method obtained using one pseudo-well, three pseudo-wells, and five pseudo-wells, respectively. The red arrows indicate the locations of oil-bearing reservoirs in the model data. From the figures, it is evident that the simulation results improve as the number of pseudo-wells increases, for both the Kriging method and the method proposed in this paper. In the case of fewer pseudo-wells, the Kriging method fails to capture the spatial variation of the elastic parameters. However, even with a small number of pseudo wells, the method in this paper is still able to effectively utilize the seismic waveform information, leading to a more accurate simulation of elastic parameters at oil-bearing locations. Moreover, the simulated small-angle elastic impedance follows the trend of the seismic waveform, which aligns with the spatial structural characteristics of the underground medium. The relative errors of Kriging’s method are 0.25, 0.24, and 0.22 for the simulation results of the oil-bearing layers at the red arrows under one pseudo-well, three pseudo-wells, and five pseudo-wells, respectively. For the method proposed in this paper, the relative errors are 0.21, 0.14, and 0.10, respectively. These results demonstrate that the method in this paper yields relatively smaller relative errors in comparison to the Kriging method.

Based on the above method, the modeling effect of the compressed sensing frequency division waveform indication modeling method is further tested by taking the small-angle elastic impedance model data and synthetic seismic data as examples. We set the low frequency range to 0–20Hz, the medium frequency range to 20–45Hz, and the frequency greater than 45Hz as the high frequency range. Figure 3 presents the seismic data corresponding to the near-angle seismic data for the low-frequency (Figure 3A), medium-frequency (Figure 3C), and high-frequency (Figure 3E) bands, along with the simulation results for the near-angle elastic impedance in these frequency bands: the low-frequency (Figure 3B), medium-frequency (Figure 3D), and high-frequency (Figure 3F). As shown in the figure, the elastic impedance simulation results of different frequency bands established by the frequency division method are consistent with the spatial variation trends of seismic data in different frequency bands. And it can better simulate geological bodies of different thicknesses according to the dominant frequencies corresponding to geological bodies of different thicknesses. The results from different frequency bands are integrated in the frequency domain to obtain the final required elastic impedance simulation results. Next, we further utilize the seismic data from different angles to obtain the elastic impedance simulation results at those angles. Figure 4 displays the simulation results for the near angle (Figure 4A), middle angle (Figure 4C) and far angle (Figure 4E) elastic impedance, following the integration process described in this paper. The elastic impedance inversion results (Figures 4B, D, F) are obtained by combining the obtained elastic impedance simulation results with seismic data at different angles. Comparing the inversion results with the real model, the overall trend of the elastic impedance is better inverted and the location of the oil-bearing layer is clearly depicted. Figure 5 shows a comparison of the elastic impedance inversion results for three different angles at a single trace, with the red curve representing the real data, the green curve representing the modeled results, and the blue curve representing the inversion results. Although there are some perturbations in the inversion results, which do not exactly match the model, they still fulfill the requirements for accurate reservoir characterization. Therefore, the model test confirms the validity of the method proposed in this paper.

To further validate the advantages of this paper’s method, we perform near angle elastic impedance inversion using the Kriging simulation results from Figure 2E, and compare them with the inversion results obtained using the method outlined in this paper (Figure 4B). Figure 6A shows the inversion results based on the Kriging simulation, and Figure 6B presents box plots of the relative errors for both the inversion results from this paper and those based on the Kriging simulation. As shown in the figure, the relative errors of the red solid line (inversion of the method of this paper) are all smaller than those of the blue solid line (inversion based on the Kriging simulation). The mean relative error for the Kriging-based inversion is 0.040, while the mean relative error for the proposed method is 0.015. This quantitative analysis demonstrates that the inversion results from the method in this paper exhibit improved accuracy, which further confirming its effectiveness.

To further validate the applicability of this method to complex models, we use the Stanford model for testing. This model is a fluvial channel system containing deltaic deposits and meandering channels. A section in the Stanford model is taken here for testing. Figure 7A presents the near angle elastic impedance of the real model, in which the oil-bearing reservoir shows a low elastic impedance. Figure 7B displays the near angle seismic data and Figure 7C shows the extracted pseudo-well data. We set the low frequency range to 0–20Hz, the medium frequency range to 20–45Hz, and frequencies greater than 45Hz for the high frequency range. The simulation results for different frequency bands are further established. Figure 8 shows the seismic data for the low (Figure 8A), medium (Figure 8C) and high (Figure 8E) frequency bands, along with the near-angle elastic impedance simulations for each frequency range: low (Figure 8B), medium (Figure 8D) and high (Figure 8F). The simulation results for different frequencies are consistent with the trend of seismic waveform changes at the corresponding frequencies. These results are integrated to obtain the full frequency band simulation, which is then compared with the full-frequency results obtained using the Kriging method. Figure 9A shows the full frequency band elastic impedance simulation results from the Kriging’s method and Figure 9B presents the corresponding results from the proposed method. We calculate the relative errors for the oil-bearing layers based on the two inversion results. For the simulation of oil-bearing reservoir segments, the mean relative error for the Kriging’s method is 0.21, while the mean relative error for the proposed method is 0.13.

The inversion is further performed based on the two simulation results, respectively. Thus, the inversion results based on the Kriging method (Figure 9C) and the inversion results based on the proposed method (Figure 9D) are obtained. A comparison of the two inversion results shows that the proposed method more clearly delineates the boundaries and distribution of the reservoir. Additionally, when comparing the mean relative errors between the calculated inversion results and the real model, the relative error of the proposed method (0.055) is smaller than that of the Kriging method (0.083). Therefore, this complex model test further verifies the feasibility and advantages of this method.

We conduct the test using field data from a work area in southern China. Figure 10 displays the near angle (Figure 10A), middle angle (Figure 10B) and far angle stacked seismic data (Figure 10C) of the field data. The method proposed in this paper is evaluated using far angle elastic impedance simulations as an example. For this field data test, we also set the low frequency range as 0–20Hz, the medium frequency range as 20–45Hz, and the frequency greater than 45Hz as the high frequency range. Figure 11 shows the low- (Figure 11A), medium- (Figure 11C), and high-frequency (Figure 11E) seismic data corresponding to the far angle seismic data, as well as the results of the far angle elastic impedance simulations for low- (Figure 11B), medium- (Figure 11D), and high-frequency (Figure 11F). It can be observed that the trends of the far angle elastic impedance simulation results at different frequencies are consistent with the trends of the seismic data at different frequencies. The far angle elastic impedance at low, medium and high frequencies is well simulated. However, the low-frequency seismic data at shallow depths exhibit a low signal-to-noise ratio, leading to less continuity in the elastic impedance simulation results at shallow locations. Therefore, the continuity of the simulation results under waveform indication is also limited by the continuity and signal-to-noise ratio of seismic data. Improving the continuity of the simulation results should be considered as the first step to improve the signal-to-noise ratio and continuity of the seismic data.

The elastic impedance simulation results at different frequencies are integrated in the frequency domain to obtain the final elastic impedance simulation results at three angles. Figure 12 presents the elastic impedance simulation results for the near angle (Figure 12A), middle angle (Figure 12C) and far angle (Figure 12E), respectively. Among them, the elastic impedance simulation results at different angles are all calculated from the corresponding seismic data at different angles. The final elastic impedance inversion results are obtained by further inversion using seismic data from various angles. Figure 12 also shows the elastic impedance inversion results at near angle (Figure 12B), middle angle (Figure 12D) and far angle (Figure 12F). In these inversion results, the black curve represents the elastic impedance curve from the logging data, and the black elliptical box indicates the location of the gas-bearing reservoir. According to the previous statistical analysis of the logging data, the elastic impedance at the gas-bearing reservoir location shows a low value. From the inversion results, it is evident that the elastic impedance inversion results from three different angles exhibit low values at the gas-bearing reservoir location. And the inversion results match well with the logging curve and exhibit good lateral continuity. Figure 13 displays the well side trace inversion results and logging data histograms at different angles. The histograms reveal that both the well side trace inversion results and the logging curves are consistent with the Gaussian distribution, and the distribution trend and distribution pattern are consistent. At the same time, we also calculate the relative errors between the well data and the inversion results. Among them, the relative error average for the small-angle, medium-angle, and large-angle elastic impedance are 0.020, 0.018 and 0.026, respectively. The relative errors between the three-angle elastic impedance inversion results and the well data are all small. This also indirectly illustrates the accuracy of the inversion results. Therefore, the field data test effectively demonstrates the accuracy and applicability of the proposed method.