Direct Inversion Method of Brittleness Parameters Based on Reweighted Lp-Norm

Abstract

Brittleness is an important factor that indicates shale properties, as well as fracturability, and it can be well-represented using the elasticity parameter Eρ. Seismic inversion allows direct access to the Eρ ofbrittleness parameters. Seismic inversion is a typical ill-posed problem that has an enormous multiplicity of solutions. In order to invert reservoir brittleness parameters more stably and reliably, a direct inversion method for determining brittleness parameters based on a reweighted Lp-norm is proposed, and the reweighted Lp method is introduced to brittleness parameter inversion for the first time. The alternating direction method of multipliers (ADMM) is used to establish the inversion structure and to optimize the objective function in blocks, which effectively improves the convergence speed. We first introduce a reweighted Lp method and establish a pre-stack inversion objective function based on the reweighted Lp method. Then, theoretical simulation data are applied to compare the inversion outcomes of the new method with those of the traditional method, and the effect of the method in this manuscript is verified. Finally, the feasibility of this method is further verified using actual data for experimental analysis. Through an analysis of the experimental results, we find that this method can be well-applied to seismic pre-stack inversion calculation and provides a new direct inversion method for the determination of brittleness parameters for exploration geophysics.

1. Introduction

The brittleness index is a critical factor for indicating the properties of shale and the effect of fracturing. There are two main methods to evaluate brittleness: one is to evaluate brittleness using mineral content or experimental stress results, and the other is to express brittleness with elastic parameters. In terms of elastic parameters, many scholars have conducted studies to evaluate brittleness and have derived brittleness representation formulas [1,2,3]. Seismic pre-stack inversion is an effective method of obtaining elastic parameters. The Zoeppritz equation is the basis of pre-stack inversion, but its equations are too complicated and inconvenient to use, so the majority of scholars have approximated it [4,5]. At the same time, the Zoeppritz approximation equation provides an opportunity to directly invert brittleness parameters. Zong et al. [6] derived a YPD equation based on the Aki–Richards approximate formula; established a linear relationship between the reflection coefficient and Young’s modulus, Poisson’s ratio, and density; and provided a direct inversion method for brittleness parameters. Zhang et al. [7] derived approximate equations for compression and transition wave reflection coefficients based on Young’s modulus multiplied by density ($E\rho$), as well as Poisson’s ratio and density, providing a new way to obtain direct inversion. After quantitative studies on the lithology and brittleness of shale gas, it was concluded that the $E\rho$ was a better indicator of shale properties [8]. Wang [9] found that the $E\rho$ could be used as a more sensitive elastic parameter to calculate the brittleness index and proposed a new formula to calculate brittleness using the $E\rho$ instead of Young’s modulus. Zhang et al. [10] proposed a binomial equation based on the relationship between Young’s modulus, Poisson’s ratio, and density. Ge et al. [11] derived a new three-property parametric reflection coefficient approximation equation for the prediction of the brittleness index.

Seismic inversion is a typical ill-posed problem. In 1943, Tikhonov explained that this problem could be solved, and in 1963, the Tikhonov regularization was proposed [12]. The Tikhonov regularization is an additional term. The regularization term is a priori information that transforms an ill-posed problem into a conditionally posed problem to be solved, and the result of its solution is biased toward the prior information, thus mitigating the awkward posing of inverse problems. Rudin et al. [13] proposed total variation (TV) regularization and applied it to image denoising. Wang [14] proposed a $P\text{-}wave$ impedance inversion method based on L1 regularization. Zhang et al. [15] proposed a sparse $P\text{-}wave$ impedance inversion method with TV regularization. Liu et al. [16] proposed a $P\text{-}wave$ impedance inversion method with L1 regularization. Li et al. [17] introduced anisotropic total variational regularization based on the Lp-norm (ATpV) to achieve impedance inversion, which led to smaller inversion errors and improved inversion results. Chen et al. [18] used the Lp-norm to express reflection coefficients and found that the Lp-norm improved the accuracy more than the L1-norm. Wu et al. [19] introduced mixed second-order fractional anisotropic total p-variation regularization (MS_ATpV) to $P\text{-}wave$ impedance inversion, which reduced the number of multiple solutions and improved the inversion accuracy. Liu et al. [20] proposed hybrid total variational (HTV) regularization, which achieved better results for inversion. Currently, regularization methods are often used in seismic inversions to mitigate ill-posed problems [21,22,23,24,25,26,27]. They can effectively improve the accuracy and stability of seismic inversion.

To better recover the sparsity of seismic inversions, Candes et al. [28] proposed a reweighted L1 method, which could better penalize non-zero coefficients. Cristian et al. [29] described a new method of sparse filtering based on a reweighted L1 algorithm. Ma et al. [30] applied a reweighted L1 method to synthetic aperture radar imaging. Feng et al. [31] proposed a reweighting method to improve the precision of inverse synthetic aperture radar imaging. Song et al. [32] proposed a new edge-preserving image-smoothing algorithm based on a reweighted L1 algorithm. He et al. [33] proposed an impedance inversion method based on a reweighted L1-norm constraint, which could obtain more accurate impedance boundaries and attenuate the pseudo-layer phenomenon. A reweighting algorithm constrains a regularization term by weight, which eliminates the effect of the amplitude of the regularized parametric term on the sparsity, yielding results with better sparsity than results with no weighting.

In this manuscript, we provide a direct inversion method for brittleness parameters based on a reweighted Lp-norm that introduces a reweighted Lp-norm to brittleness parameter inversion for the first time and recovers better result sparsity. The alternating direction method of multipliers (ADMM) is also used for establishing the inversion structure to optimize the objective function in blocks, which accelerates the convergence speed. In the following sections, the theory of the algorithm and the derivation process of the formula are introduced, and the computational effectiveness of the method is tested with theoretical and field data.

2. Materials and Methods Used for Experimentation

2.1. Reweighted Lp Method

First, the Zoeppritz approximate equation [7] used in this manuscript is introduced, and its expression is as follows:

$$R(\theta )=a(\theta )R_{\_E\rho }+b(\theta )R_{\_\sigma }+c(\theta )R_{\_\rho }$$

(1)

$$a(\theta )=\frac{1}{4}{\sec }^2\theta -2k^2{\sin }^2\theta$$

(2)

$$b(\theta )=\frac{1}{4}{\sec }^2\theta \frac{(2k^2-3){(2k^2-1)}^2}{k^2(4k^2-3)}+2k^2{\sin }^2\theta \frac{1-2k^2}{3-4k^2}$$

(3)

$$c(\theta )=\frac{1}{2}+2k^2{\sin }^2\theta -\frac{1}{2}{\sec }^2\theta$$

(4)

where $\theta$ is the angle of incidence; $E\rho$ is Young’s modulus multiplied by density; $\sigma$ is Poisson’s ratio; $\rho$ is density; k is the ratio of shear to compressional wave velocity ($k=\frac{V_S}{V_P}$); and $R_{\_E\rho }$, $R_{\_\sigma }$, and $R_{\_\rho }$ are the reflectivity of the $E\rho$, Poisson’s ratio, and density, respectively. To obtain the required parameters from the seismic data, it is necessary to establish the relationship between the seismic data and the parameters. Taking $\rho$ as an example, the reflection coefficient of $\rho$ can be approximated as follows:

$$R_{\_\rho }=\frac{{\rho }_{{}_{i+1}}-{\rho }_{{}_i}}{{\rho }_{{}_{i+1}}+{\rho }_{{}_i}}\approx \frac{1}{2}\Delta \ln (\rho )=\frac{1}{2}\left[\ln ({\rho }_{{}_{i+1}})-\ln ({\rho }_{{}_i})\right]=\frac{1}{2}\left(L_{i+1}-L_i\right)$$

(5)

where $L_i=\ln ({\rho }_{{}_i})$.

$R_{\_\rho }$ is represented by the following matrix:

$$R_{\_\rho }=\frac{1}{2}D{L}_{\rho }$$

(6)

where $D$ is the difference matrix, and $L_{\rho }$ is the natural logarithm of density.

According to Equation (1) with the use of a wavelet kernel matrix $W(\theta )$, the seismic record $S(\theta )$ can be shown as follows:

$$\left[\begin{array}{c}S(\theta )\\ ⋮\\ S({\theta }_M)\end{array}\right]=\frac{1}{2}\left[\begin{array}{ccc}W({\theta }_1)a({\theta }_1)D& W({\theta }_1)b({\theta }_1)D& W({\theta }_1)c({\theta }_1)D\\ ⋮& ⋮& ⋮\\ W({\theta }_M)a({\theta }_M)D& W({\theta }_M)b({\theta }_M)D& W({\theta }_M)c({\theta }_M)D\end{array}\right]\left[\begin{array}{c}{L}_{E\rho }\\ L_{\sigma }\\ L_{\rho }\end{array}\right]$$

(7)

where $L_{E\rho }$ is the natural logarithm of the $E\rho$, $L_{\sigma }$ is the natural logarithm of Poisson’s ratio, and $L_{\rho }$ is the natural logarithm of density.

The forward equation based on an approximate linear formula can be abbreviated as follows:

$$S=WBDL=AL$$

(8)

where $B$ is the matrix consisting of Equations (2)–(4), and $L$ is the matrix of the natural logarithms of the $E\rho$, Poisson’s ratio, and density.

As shown in the introduction, we first introduced a weighted Lp-norm into the inversion of brittleness parameters. This method can solve the shortcomings of L1 regularization in inverse objective function. The regularization term is weighted and constrained by the Lp-norm, which yields results with better sparsity and improves the inversion accuracy. The reweighted Lp method for inversion parameters can be defined as follows:

$$Lp(L)=\lambda \left(w_i{\Vert L\Vert }_1\right)$$

(9)

where ${\Vert \Vert }_1$ is the $L1-norm$, $\lambda$ is the regularization factor, and $w_i$ is the weight value.

$${(w_i)}^{l+1}=\frac{1}{{\Vert {\left(L_i\right)}^{l+1}\Vert }_p^p+\epsilon }$$

(10)

where ${\Vert \Vert }_p$ is the $Lp-norm(0\le p<1)$, and $\epsilon$ is a constant greater than 0. $\epsilon$ can ensure the stability of the algorithm, that is, if $L$ has a value of zero in a round of iterations, there is a non-zero value in the corresponding position in the next iteration.

2.2. Direct Inversion Method of Brittleness Parameters Based on Reweighted Lp

In constructing the objective function, a reweighted Lp method was introduced for the first time to constrain the inverse objective function. Meanwhile, the ADMM method was used to minimize the objective function. According to Equation (8), the traditional objective function can be expressed as follows:

$$J(L)=\underset{L}{\min }{\Vert AL-S^{obs}\Vert }_2^2$$

(11)

where ${\Vert \Vert }_2$ is the $L2-norm$, $S^{obs}$ is the observed data, adding the initial model $L_0$ constraint and the reweighted Lp method constraint to the objective function. The regularization term is a priori information. The ill-posed problem is transformed into a conditionally, well-posed problem for solution. The solution results are biased toward the prior information to reduce the uncertainty of the inversion. The reweighted Lp method contains amplitude information, and the constrained objective function can better restore sparsity. Objective function can be expressed as follows:

$$J(L)=\underset{L}{\min }{\Vert AL-S^{obs}\Vert }_2^2+\mu {\Vert L-L_0\Vert }_2^2+\lambda (w_i{\Vert DL\Vert }_1)$$

(12)

where $\mu$ is the initial model regularization parameter.

The matrix $R$ is introduced to replace the matrix $DL$. The objective function is further shown as follows:

$$\begin{array}{l}J(L)=\underset{L}{\min }{\Vert AL-S^{obs}\Vert }_2^2+\mu {\Vert L-L_0\Vert }_2^2+\lambda (w_i{\Vert R\Vert }_1)\\ s.t.R=DL\end{array}$$

(13)

Then, the secondary penalty term and dual variable are introduced to change the above objective function into a constrained optimization problem:

$$J(L)=\underset{L}{\min }{\Vert AL-S^{obs}\Vert }_2^2+\mu {\Vert L-L_0\Vert }_2^2+\lambda (w_i{\Vert R\Vert }_1)+\eta {\Vert R-DL+C\Vert }_2^2$$

(14)

where $\eta$ is the difference of the Lagrangian constraint parameters, and $C$ is the dual term.

We used the ADMM algorithm to minimize the objective function. Based on the ADMM algorithm, the objective function can be decomposed into the sub-objective functions of $L$, $R$, and $C$, as follows:

$$J(L)=\underset{L}{\min }{\Vert AL-S^{obs}\Vert }_2^2+\mu {\Vert L-L_0\Vert }_2^2+\eta {\Vert R-DL+C\Vert }_2^2$$

(15)

$$J(R)=\underset{R}{\min }\lambda (W_i{\Vert R\Vert }_1)+\eta {\Vert R-DL+C\Vert }_2^2$$

(16)

$$J(C)=\underset{C}{\min }\eta {\Vert R-DL+C\Vert }_2^2$$

(17)

We used the reweighted Lp method to construct the inversion objective function. The reweighted Lp algorithm is shown in Equations (9) and (10), and the inversion objective function is constructed in Equation (14). The ADMM algorithm decomposes the objective function into three sub-objective functions, namely, Equations (15)–(17). The inversion results are obtained by solving the sub-objective function. Equations (7) and (8) are forward problems, and Equation (14) is an inverse problem. Figure 1 shows the inversion flow according to which the theoretical and actual data inversion experiments were carried out.

3. Experiments

In order to verify the feasibility of the reweighted Lp method, in this section, the reweighted Lp method and the traditional method were each applied to the inversion of the theoretical model and actual data, respectively, for a series of comparative analyses.

3.1. Theoretical Model Experiment

The model used in the theoretical experiments was part of the Marmousi2 model. To carry controlled tests, we generated a synthetic dataset that was used as observational data in the present work. Such a dataset is generated by using the same mechanism to solve the forward problem. Figure 2 shows the angle stack data volume composed of reflection coefficient and wavelet convolution. The observed data made a synthesis according to Equation (7). Figure 2a–c show the angle superposition data volumes for the incident angles of 10°, 20°, and 30°, respectively. The theoretical models are shown in Figure 3, which includes three parameters: Figure 3a shows the $E\rho$ theoretical model, Figure 3b shows the Poisson’s ratio model, and Figure 3c shows the density theoretical model. The models have 800 traces in the cross-sections, and each trace includes 500 sampling points.

We used the traditional method and the weighted Lp method to invert the theoretical models. Figure 4a, Figure 5a and Figure 6a show the inversion results based on the traditional method for the $E\rho$, Poisson’s ratio, and density. Figure 4b, Figure 5b and Figure 6b show the inversion results based on the reweighted Lp method for the $E\rho$, Poisson’s ratio, and density. Compared with the theoretical models, it was demonstrated that the inversion results of the reweighted Lp method had good continuity in the horizontal orientation and clear stratigraphic boundaries in the vertical orientation. The new method could adequately reflect the geological characteristics of the models and was highly consistent with the theoretical models. Compared with the traditional method, the reweighted Lp method showed approximately the same characteristics. However, through careful observation, it was found that the reweighted Lp method had higher accuracy and better stability. In the black rectangular boxes, the reweighted Lp method had better continuity and was closer to the theoretical models, which further proved the significance of the method’s improvement.

In order to reflect the difference between the two methods more clearly, we calculated the absolute error values of inversion for the two methods. Figure 7, Figure 8 and Figure 9 show the absolute error values of inversion for the two methods. Figure 7a, Figure 8a and Figure 9a show the absolute error values of inversion for the traditional approach for the $E\rho$, Poisson’s ratio, and density, respectively. Figure 7b, Figure 8b and Figure 9b show the absolute error values of inversion for the reweighted Lp method for the $E\rho$, Poisson’s ratio, and density, respectively. It can be seen from the three figures that the errors of the inversion results for the proposed method were significantly reduced, and the results were more stable, especially those in the black rectangular boxes. The above experiments show that this research method could improve inversion accuracy and had important significance for improvement.

Then, we analyzed the inversion results of the 150th trace and the 500th trace. Figure 10 shows the inversion results of the 150th trace, and Figure 11 shows the inversion results of the 500th trace, where the black, green, blue, and red curves represent the actual model, initial model, inversion outcome of the traditional approach, and inversion outcome of the reweighted Lp method, respectively. Figure 12 shows the absolute error values of inversion for the 150th trace, and Figure 13 shows the absolute error values of inversion for the 500th trace, where the blue and red curves represent the absolute error of inversion for the traditional approach and the absolute error of inversion for the reweighted Lp method, respectively. From the individual trace inversion outcomes and absolute errors, it can be seen that the line of the inversion outcome of the reweighted Lp method was consistent with the line of the actual model. The method used was close to the actual inversion results, with high accuracy and stability, which further proves the validity of the method.

In order to carry out the noise resistance experiment, Gaussian white noise was added to seismic records, and seismic records with a signal-to-noise ratio of 5 were synthesized. Figure 14 shows the angle stack data volume after adding noise, and Figure 14a–c express the angle stack data volume with incident angles of 10°, 20°, and 30°, respectively. Figure 15, Figure 16 and Figure 17 show the inversion results of the weighted Lp method and the traditional method after adding noise. Figure 15a,b express the $E\rho$ inversion results of the two methods. Figure 16a,b show the Poisson’s ratio inversion results of the two methods. Figure 17a,b show the density inversion results of the two methods. It can be seen from the figures that the inversion results of the traditional approach had more outliers, which was more obvious in the black rectangular boxes. The experiments show that, compared with the traditional method, the reweighted Lp method had stronger antinoise ability.

We then calculated the absolute error values of inversion for the two methods separately, as shown in Figure 18, Figure 19 and Figure 20. Figure 18a,b show the $E\rho$ absolute error values of inversion for the two methods. Figure 19a,b show the Poisson’s ratio absolute error values of inversion for the two methods. Figure 20a,b show the density absolute error values of inversion for the two methods. From the absolute error plots, it can be seen that the error values of the new method were smaller, with relatively few anomalies near the fault area and higher accuracy. After that, in order to analyze the experimental results for antinoise ability from a microscopic perspective, data for the 150th trace and the 500th trace were extracted to compare the two methods. The inversion results of the two methods are expressed in Figure 21 and Figure 22, and the absolute error values of inversion are expressed in Figure 23 and Figure 24. In Figure 21 and Figure 22, the black, green, blue, and red curves represent the actual model, initial model, inversion outcomes of the traditional approach, and inversion outcomes of the reweighted Lp method, respectively. In Figure 21 and Figure 22, the inversion outcomes of the reweighted Lp method were in better agreement with the actual model. In Figure 23 and Figure 24, the blue and red curves represent the absolute error of inversion for the traditional approach and the absolute error of inversion for the reweighted Lp method, respectively. The accuracy of the two methods is more clearly shown in Figure 23 and Figure 24, and the inversion results of the new method could maintain the minimum error value and were more stable at most sampling points. The antinoise ability experiments show that the reweighted Lp method had stronger noise resistance and better precision, which further illustrates the significance of the method’s improvement.

To consider the effect of the regularization factor $\lambda$ on the experimental results, different values were adjusted for the inversion experiments. The root mean square error (RMSE) of the inversion results was also calculated, and the value of $\lambda$ was determined when the root mean square error was smallest. In Figure 25, a graph of the variation in the RMSE error with A is shown. Figure 26 shows the convergence curve of the objective function.

.

3.2. Field Data Experiment

The field data were collected from Sichuan, China. The angle stack data volume is shown in Figure 27, with incidence angles of 9°, 18°, and 27° in Figure 27a–c, respectively. The section had 800 traces with a vertical time length of 750 ms and a sampling interval of 2 ms. The two red lines represent the two layers of geological interpretation, and the black vertical line is the well trajectory. The inversion method was applied to the field data, and Figure 28, Figure 29 and Figure 30 show the inversion results of the reweighted Lp and the traditional method, respectively. Figure 28a,b express the $E\rho$ inversion results of the two methods. Figure 29a,b show the Poisson’s ratio inversion results of the two methods. Figure 30a,b show the density inversion results of the two methods. For the sake of observation, the black rectangular boxes are enlarged in the figures. In these areas, compared with the traditional approach, the reweighted Lp method had better transverse continuity, while the vertical resolution was effectively improved. In addition, the overall effect in all the other regions was greatly improved and could effectively portrayed stratigraphic features.

For the purpose of verifying the reliability of the inversion method, we compared the logged data with the inversion data of the well location trace. The inversion comparison chart is shown in Figure 31. In order to more obviously reflect the difference between the two methods, we calculated the absolute error values of inversion for the two methods, as shown in Figure 32. As can be seen in the two figures mentioned above, both methods were in good agreement with the logged data. However, when carefully observed, the inversion results of the reweighted Lp method were closer overall to the logged data and were more stable. The black rectangular boxes indicate the reservoir location. In this area, compared with the traditional method, the reweighted Lp method could invert the characteristics of the reservoir more accurately, with less error and higher accuracy. This further authenticated the validity of the method and the significance of the improvement.

4. Conclusions

In this research, we proposed a direct inversion method for brittleness parameters based on a reweighted Lp-norm. The regularization constraint in this method could reduce the number of multiple solutions to the inversion problem and could improve inversion stability. The reweighted Lp method effectively restored the sparsity of the inversion results and improved inversion accuracy by constraining the regularization term. In addition, the inversion framework established by the ADMM method improved the convergence speed.

The method was validated by theoretical models and field data. The inversion results expressed that the reweighted Lp method could reduce the inversion error, improve the vertical resolution, enhance the lateral continuity, and describe the reservoir information more accurately. The proposed methodology can be applied in different contexts, both in post-stacked data inversion processes [34] and also in inverting pre-stacked data [35], and it provides a new direct inversion method for brittleness parameters for exploration geophysics.