Seismic Reflection Coefficient Inversion Using Basis Pursuit Denoising in the Joint Time-Frequency Domain

: Seismic reﬂection coe ﬃ cient inversion in the joint time-frequency domain is a method for inverting reﬂection coe ﬃ cients using time domain and frequency domain information simultaneously. It can e ﬀ ectively improve the time-frequency resolution of seismic data. However, existing research lacks an analysis of the factors that a ﬀ ect the resolution of inversion results. In this paper, we analyze the inﬂuence of parameters, such as the length of the time window, the size of the sliding step, the dominant frequency band, and the regularization factor of the objective function on inversion results. The SPGL1 algorithm for basis pursuit denoising was used to solve our proposed objective function. The applied geological model and experimental ﬁeld results show that our method can obtain a high-resolution seismic reﬂection coe ﬃ cient section, thus providing a potential avenue for high-resolution seismic data processing and seismic inversion, especially for thin reservoir inversion and prediction.


Introduction
Seismic inversion is a popular research topic in the field of seismic exploration [1][2][3]. With the development of seismic exploration, the identification of thin reservoirs has become increasingly important. Due to the limitations of seismic vertical resolution, traditional seismic inversion methods cannot accurately identify and predict thin reservoirs when the reservoir thickness is less than 1/4 wavelet length [4][5][6].
In order to overcome the shortcomings of traditional seismic inversion methods in thin reservoir identification and prediction, researchers have proposed a thin layer inversion method based on spectral decomposition [7,8]. Subsequently, the concept of using a parity decomposition was also proposed to enhance the seismic resolution [9,10]. These studies show that the even component of seismic data helps to improve the seismic vertical resolution, but the odd component reduces the resolution of seismic data. Therefore, the resolution of seismic data can be improved by increasing the proportion of even components in the reflection coefficient [11][12][13][14].
In recent years, a large number of sparse optimization methods based on the compressed sensing theory have been proposed. These methods can also be used for solving sparse constraint problems in seismic inversion [15][16][17]. Based on the sparsity of reflection coefficients, many seismic sparse inversion methods have been proposed, such as sparse deconvolution [18], sparse acoustic impedance [19][20][21], sparse de-noising [22][23][24], sparse elastic parameters inversion [25], etc. In the field of seismic sparse inversion, the basis pursuit denoising (BPDN) algorithm is widely used due to its strong anti-noise performance [26]. Chai et al. (2014) proposed a sparse reflection coefficient inversion method for non-stationary seismic data [27]. To enhance the resolution of the inversion result, Zhang et al. (2013) suggested a prestack basis pursuit seismic inversion method. They also used the BPDN algorithm in reflection coefficient inversion [28]. Yin et al. (2015) utilized the BPDN technique in seismic inversion to determine the brittleness of shale and obtained a high-resolution result [29]. Chao et al. (2017) applied the curvelet transform in seismic inversion and obtained an inversion method based on the BPDN algorithm [30]. Chai et al. (2019) proposed a sparse seismic deconvolution model by dictionary learning and used the BPDN method to solve the model [31].
Furthermore, joint time-frequency analysis can obtain information from seismic data and this method is widely used in geophysics, e.g., Hilbert Huang transform [32], fractional S transform [33], Wigner-Ville transform [34], fractional Fourier transform [35] and so on. In seismic inversion, Zhang (2017) found that the wavelet transform improved the correlation between the seismic record and synthetic seismograms. They used this transform in impedance inversion and improved the resolution of the results [36]. Hao (2018) used the Gaber transform in acoustic impedance inversion and improved the resolution of the results [37]. Hu (2019) utilized the multi-scale time-frequency domain in full waveform inversion and reduced the model dependence of results [38]. Chen (2020) applied the sparse time-frequency analysis in seismic inversion and improved the resolution of the results [39].
Inspired by these works, we propose a joint time-frequency reflection coefficient inversion method based on the BPDN algorithm, which combines the advantages of spectral decomposition and the parity decomposition theory. We used parity decomposition to directly link the time domain and the frequency domain of the reflection coefficient. This approach is able to obtain the reflection coefficient in the time domain by optimizing the objective function in the frequency domain. The vertical resolution of seismic data can be improved by properly increasing the proportion of even components in the objective function. Due to the existence of noise in seismic data, we used the BPDN algorithm to enhance the anti-noise ability [30], thus ensuring the quality and stability of the inversion results. Although there are many factors that affect the inversion result of the reflection coefficient, few studies have analyzed the effects of these factors on the inversion results. Chai et al. (2014) studied the influence of noise and the Q value on the nonstationary seismic reflection inversion [27]. Zhang and Castagna (2011) analyzed the effect of the regularization factor of the sparse constraint term and noise on the inversion result [40].  discussed the influence of noise on reflection coefficient inversion [41]. In fact, the length of the time window, the size of the sliding step, and the dominant frequency band also affect the inversion result. In this paper, we also tested and analyzed the influence of these factors on the inversion result.

The Parity Decomposition of Reflection Coefficients in the Joint Time and Frequency Domain
A simple strata model with three layers can be represented by two reflection coefficients [10] where r 1 and t 1 are the value and position of the upper reflection coefficient and r 2 and t 2 are the value and position of the bottom reflection coefficient, respectively. For convenience, the analysis point can be moved to the middle time of these two coefficients. Figure 1 shows this simple model. Now, Equation (1) can be reformulated as where T is the time interval between t 1 and t 2 .
Energies 2020, 13, x FOR PEER REVIEW 3 of 15 where T is the time interval between 1 t and 2 t . By using the Fourier transform, we can express Equation (2) in the frequency domain as where Rmodel is the expression of the reflection coefficients in the frequency domain, and f denotes the frequency. Applying the Euler equation to Equation (3) with the trigonometric function, the frequency domain equation can be written as As shown in Figure 2, the reflection coefficient series can be decomposed into an even part and an odd part by using the parity decomposition principle. A scalar representation has the following form  Equivalently, on the basis of Equations (4) and (5), the even and odd components can be expressed as By using the Fourier transform, we can express Equation (2) in the frequency domain as where R model is the expression of the reflection coefficients in the frequency domain, and f denotes the frequency. Applying the Euler equation to Equation (3) with the trigonometric function, the frequency domain equation can be written as As shown in Figure 2, the reflection coefficient series can be decomposed into an even part and an odd part by using the parity decomposition principle. A scalar representation has the following form where r e1 = r e2 = r e , and r o1 = −r o2 = r o .
Energies 2020, 13, x FOR PEER REVIEW 3 of 15 where T is the time interval between 1 t and 2 t . By using the Fourier transform, we can express Equation (2) in the frequency domain as where Rmodel is the expression of the reflection coefficients in the frequency domain, and f denotes the frequency. Applying the Euler equation to Equation (3) with the trigonometric function, the frequency domain equation can be written as As shown in Figure 2, the reflection coefficient series can be decomposed into an even part and an odd part by using the parity decomposition principle. A scalar representation has the following form  Equivalently, on the basis of Equations (4) and (5), the even and odd components can be expressed as  Equivalently, on the basis of Equations (4) and (5), the even and odd components can be expressed as 2r e = r 1 + r 2 (6) 2r o = r 1 − r 2 (7) When we consider the coefficients of the trigonometric function in Equation (4), we find that the coefficient of the cosine term is identical to the right-hand side of Equation (6) and the coefficient of the sine term is equal to the right-hand side of Equation (7). Therefore, a simpler expression can be obtained by substituting Equations (6) and (7) into Equation (4).
When we move the analysis point to the surface to conform to the actual situation, we obtain the following equation where,∆t = t 1 + T/2. The relationship between the frequency domain representation of the reflection coefficients and its time domain representation is established.
Finally, the simple model can be extended to a more complicated version with multiple layers, which is shown in Figure 3. In this case, the reflection coefficient can be expressed in the joint time and frequency domain as where ∆t = t k + T k /2, r k and r N−k+1 are the pairs of reflection coefficients, and N is the number of layers.
When we consider the coefficients of the trigonometric function in Equation (4), we find that the coefficient of the cosine term is identical to the right-hand side of Equation (6) and the coefficient of the sine term is equal to the right-hand side of Equation (7). Therefore, a simpler expression can be obtained by substituting Equations (6) and (7) into Equation (4).
When we move the analysis point to the surface to conform to the actual situation, we obtain the following equation . The relationship between the frequency domain representation of the reflection coefficients and its time domain representation is established.
Finally, the simple model can be extended to a more complicated version with multiple layers, which is shown in Figure 3. In this case, the reflection coefficient can be expressed in the joint time and frequency domain as where

The Objective Function and its Solution
According to the convolution model, a seismic trace is the convolution of a reflection coefficient series and a wavelet. Consequently, we can express this relationship in the frequency domain by using the convolution theorem, which can be expressed as

The Objective Function and Its Solution
According to the convolution model, a seismic trace is the convolution of a reflection coefficient series and a wavelet. Consequently, we can express this relationship in the frequency domain by using the convolution theorem, which can be expressed as where R real ( f ) is the reflection coefficient in the frequency domain, S( f ) and W( f ) are the Fourier transform of the seismic trace and wavelet, respectively. Based on the representation of the reflection coefficient in the frequency domain, we can construct an objective function in the frequency domain as According to the rule of the subtraction of complex numbers, we can rewrite the objective function as A discrete version of the objective function can be obtained by spreading the sum formula, which is shown as where a o and a e are the weight of the odd part and even part, respectively. A 11 and A 21 are the parameters of a e . The equations are represented as A 12 and A 22 are the parameters of a o . The equations are shown as where a o and a e are the weight of the odd part and even part, respectively. According to previous studies on spectral inversion [10], a larger even weight can improve the resolution of the spectral inversion. M is the number of sampling points in the valid frequency band. The L1 norm constraint is based on the consideration that the reflection coefficient should be sparse. λ is a regularization factor related to the noise level, which increases with the noise level. In general, given a higher λ, we can obtain a more robust inversion result.
To get the optimal solution for Equation (14), the spectral projected gradient for L1 (SPGL1), which is a type of basis pursuit denoising algorithm, was used in this paper [26,42]. Finally, the inversion result was obtained by the combination of the even component and odd component obtained from Equation (14).

Synthetic Examples
We constructed a series of reflection coefficient models to verify the validity of the seismic reflection coefficient inversion algorithm, and then we analyzed the effect of the main influencing factors on the inversion results. Figure 4a shows a single-trace reflection coefficient model with ten layers. Figure 4b shows a Ricker wavelet with a dominant frequency of 50 Hz, with a sampling frequency of 1 ms and length of 80 ms. We defined the amplitude value of the Ricker wavelet, which is used to synthesize seismic traces. The amplitude value is limited in the range of (−1, 1). Figure 4c shows the seismogram obtained from the reflection coefficient model (shown in Figure 4a) and the wavelet (shown in Figure 4b) by using the convolution operation. We applied the seismic reflection coefficient inversion algorithm to the model data to get the reflection coefficient series, which is shown in Figure 4d.
In Figure 4d, the circle markers represent the reflection coefficient shown in Figure 4a as true values and the asterisks represent the inversion result. From this figure, we can see that the inversion result agrees well with the real value. The root means square error (RMSE) of the inversion result is 1.6426 × 10 −4 , which indicates that the algorithm is able to get a very good inversion result for single-trace model data.
To evaluate the effectiveness of the seismic reflection coefficient inversion algorithm for complicated seismic data, we constructed a seismic model, which is shown in Figure 5a and the corresponding synthetic seismogram is shown in Figure 5b. The seismogram shown in Figure 5b includes 480 traces and the length of each trace is 380 ms and the sampling rate is 1 ms.
We applied the seismic reflection coefficient inversion algorithm to this seismogram and obtained the inversion result, which is shown in Figure 6a. The wavelet is a Ricker wavelet with a dominant frequency of 38 Hz, a sampling rate of 1 ms and length of 81 ms. Then we calculated the synthetic seismogram using the inversion result and wavelet. The amplitude spectrum of the synthetic seismogram is shown in Figure 6b.
All the strata shown in Figure 5a can be found in the inversion result, and the lateral continuity of the reflection coefficient section is high. Besides, the amplitude spectrum of the synthetic seismogram is significantly widened compared with that of the original seismic data.
The seismic data used in the above numerical test do not contain any noise whereas the actual seismic data is inevitably affected by noise. In order to test the anti-noise ability of the proposed approach, we constructed a model with two reflection coefficients, which is shown in Figure 7a. Then, white noise with an amplitude equal to 5% of the seismic amplitude was added to the synthetic seismic data. We set λ (the regularization factor) to be 0.01, 0.04 and 0.08. The corresponding inversion results are showed in Figure 7b-d. Note that the type and parameters of wavelet used for the synthetic seismic data is consistent with those depicted in Figure 4.
As can be seen in Figure 7b-d, because of the noise, when λ is a small value, the inversion result is greatly influenced by the noise, and there are many error points in the obtained reflection coefficient. As the regularization factor λ increases, the inversion result improves. For example, when λ = 0.08, the inversion result is better than that when λ = 0.01 and λ = 0.04.
Energies 2020, 13, 5025 7 of 15 of 80 ms. We defined the amplitude value of the Ricker wavelet, which is used to synthesize seismic traces. The amplitude value is limited in the range of (−1, 1). Figure 4c shows the seismogram obtained from the reflection coefficient model (shown in Figure 4a) and the wavelet (shown in Figure  4b) by using the convolution operation. We applied the seismic reflection coefficient inversion algorithm to the model data to get the reflection coefficient series, which is shown in Figure 4d. In Figure 4d, the circle markers represent the reflection coefficient shown in Figure 4a as true values and the asterisks represent the inversion result. From this figure, we can see that the inversion result agrees well with the real value. The root means square error (RMSE) of the inversion result is 1.6426 × 10 −4 , which indicates that the algorithm is able to get a very good inversion result for singletrace model data.
To evaluate the effectiveness of the seismic reflection coefficient inversion algorithm for complicated seismic data, we constructed a seismic model, which is shown in Figure 5a and the corresponding synthetic seismogram is shown in Figure 5b. The seismogram shown in Figure 5b includes 480 traces and the length of each trace is 380 ms and the sampling rate is 1 ms. We applied the seismic reflection coefficient inversion algorithm to this seismogram and obtained the inversion result, which is shown in Figure 6a. The wavelet is a Ricker wavelet with a dominant frequency of 38 Hz, a sampling rate of 1 ms and length of 81 ms. Then we calculated the synthetic seismogram using the inversion result and wavelet. The amplitude spectrum of the synthetic seismogram is shown in Figure 6b.  We applied the seismic reflection coefficient inversion algorithm to this seismogram and obtained the inversion result, which is shown in Figure 6a. The wavelet is a Ricker wavelet with a dominant frequency of 38 Hz, a sampling rate of 1 ms and length of 81 ms. Then we calculated the synthetic seismogram using the inversion result and wavelet. The amplitude spectrum of the synthetic seismogram is shown in Figure 6b. All the strata shown in Figure 5a can be found in the inversion result, and the lateral continuity of the reflection coefficient section is high. Besides, the amplitude spectrum of the synthetic seismogram is significantly widened compared with that of the original seismic data. The seismic data used in the above numerical test do not contain any noise whereas the actual seismic data is inevitably affected by noise. In order to test the anti-noise ability of the proposed approach, we constructed a model with two reflection coefficients, which is shown in Figure 7a. Then, white noise with an amplitude equal to 5% of the seismic amplitude was added to the synthetic seismic data. We set λ (the regularization factor) to be 0.01, 0.04 and 0.08. The corresponding inversion results are showed in Figure 7b-d. Note that the type and parameters of wavelet used for the synthetic seismic data is consistent with those depicted in Figure 4.
Since there are many sampling points in the actual seismic data, in order to get better inversion results, we used the windowed Fourier transform to implement the time-frequency conversion. The time window in our scheme is a rectangle window. In order to test the effect of window length and sliding step size on the inversion results, we constructed a reflection coefficient model, which is Since there are many sampling points in the actual seismic data, in order to get better inversion results, we used the windowed Fourier transform to implement the time-frequency conversion. The time window in our scheme is a rectangle window. In order to test the effect of window length and Energies 2020, 13, 5025 9 of 15 sliding step size on the inversion results, we constructed a reflection coefficient model, which is shown in Figure 8a. Similarly, we added 5% white noise to the synthetic seismic record. The relationship between the time window length WL and the sliding step size SP was set as SP = WL, SP = 0.5 WL, and SP = 0.25 WL. From these inversion results, it can be seen that the method that obtains the best result is when the sliding step size is equal to 1/2 of the time window length, and that an increase or decrease in the size of the sliding step leads to inaccuracy or even erroneous reconstruction of the reflection coefficients at some sampling points.
Another important factor that affects the inversion result is the selection of the dominant frequency band (f1 through fM in Equation (14)). In the proposed inversion algorithm, we selected part of the frequency band of the seismic data and that of the wavelet instead of the full frequency band to improve the anti-noise ability. The different selection of the frequency band for the same seismic data will cause differences in the performance of the inversion. To test the influence of the dominant frequency band on the proposed algorithm, we constructed the reflection coefficient model as shown in Figure 9a  From these inversion results, it can be seen that the method that obtains the best result is when the sliding step size is equal to 1/2 of the time window length, and that an increase or decrease in the size of the sliding step leads to inaccuracy or even erroneous reconstruction of the reflection coefficients at some sampling points.
Another important factor that affects the inversion result is the selection of the dominant frequency band (f 1 through f M in Equation (14)). In the proposed inversion algorithm, we selected part of the frequency band of the seismic data and that of the wavelet instead of the full frequency band to improve the anti-noise ability. The different selection of the frequency band for the same seismic data will cause differences in the performance of the inversion. To test the influence of the dominant frequency band on the proposed algorithm, we constructed the reflection coefficient model as shown in Figure 9a. The model contains white noise with 5% of the maximum amplitude of the seismic record. The dominant frequency band of the seismic data is about 5-80 Hz. Due to the existence of a large amount of noise in the high-frequency component of the signal, when the selected frequency band is too wide, the high-frequency noise remains in the seismic data, resulting in the degradation of the inversion result. On the contrary, if the frequency band is too narrow, a lot of frequency information will be missed, thus causing the obtained reflection coefficient to deviate from its true value. Only when the selected frequency band agrees with the dominant frequency band of the seismic data and the wavelet, can we obtain an inversion result of high quality. As shown in Figure 9c, the amplitude of the obtained reflection coefficient is exactly in accordance with the real value, and the inversion error is substantially reduced compared to Figure 9b. In addition, we can further reduce the inversion error by properly adjusting the regularization factor λ .
As shown in the numerical examples above, the influence of noise, the dominate frequency band, the length of the time window, the size of the sliding step and other factors must be considered in order to obtain a high-quality inversion result. The improper setting of any of these factors will lead to an inaccurate inversion result. In practice, we need to select reasonable inversion parameters according to the characteristics of the seismic data, such as the dominant frequency band range, the noise level, the wavelet length, and so on, to obtain a high-quality inversion result. Due to the existence of a large amount of noise in the high-frequency component of the signal, when the selected frequency band is too wide, the high-frequency noise remains in the seismic data, resulting in the degradation of the inversion result. On the contrary, if the frequency band is too narrow, a lot of frequency information will be missed, thus causing the obtained reflection coefficient to deviate from its true value. Only when the selected frequency band agrees with the dominant frequency band of the seismic data and the wavelet, can we obtain an inversion result of high quality. As shown in Figure 9c, the amplitude of the obtained reflection coefficient is exactly in accordance with the real value, and the inversion error is substantially reduced compared to Figure 9b. In addition, we can further reduce the inversion error by properly adjusting the regularization factor λ.
As shown in the numerical examples above, the influence of noise, the dominate frequency band, the length of the time window, the size of the sliding step and other factors must be considered in order to obtain a high-quality inversion result. The improper setting of any of these factors will lead to an inaccurate inversion result. In practice, we need to select reasonable inversion parameters according to the characteristics of the seismic data, such as the dominant frequency band range, the noise level, the wavelet length, and so on, to obtain a high-quality inversion result.

Field Data Example
In order to evaluate the practical use of our algorithm, we applied the proposed method to 2-D field seismic data, which are shown in Figure 10a. The area inside the red ellipse represents the gas reservoir (target zone), and it is also the area that the interpreter focused on. The seismic data were acquired from the Sichuan basin, Northeast Sichuan, China, and consist of a common depth point (CDP) gather of 200 traces. The length of each seismic trace and wavelet are 400 ms and 81 ms respectively, and the sampling interval is 1 ms. The parameters are as follows: the regularization factor λ = 2, the length of the time window is 120 ms, the sliding step is 250 ms, and the dominant frequency band is 5-90 Hz. The inversion result is shown in Figure 10b. It can be seen that the resolution of the obtained reflection coefficient section is very high and the lateral continuity is very good.

Field Data Example
In order to evaluate the practical use of our algorithm, we applied the proposed method to 2-D field seismic data, which are shown in Figure 10a. The area inside the red ellipse represents the gas reservoir (target zone), and it is also the area that the interpreter focused on. The seismic data were acquired from the Sichuan basin, Northeast Sichuan, China, and consist of a common depth point (CDP) gather of 200 traces. The length of each seismic trace and wavelet are 400 ms and 81 ms respectively, and the sampling interval is 1 ms. The parameters are as follows: the regularization factor 2 λ = , the length of the time window is 120 ms, the sliding step is 250 ms, and the dominant frequency band is 5-90 Hz. The inversion result is shown in Figure 10b. It can be seen that the resolution of the obtained reflection coefficient section is very high and the lateral continuity is very good. We used the reflection coefficient section shown in Figure 10b and a high-frequency wavelet to perform a convolution operation, to obtain a reconstructed synthetic seismic section, which is shown in Figure 11a. Note that in order to keep the amplitude consistent with the original seismic data in Figure 10a, we calibrated the range of color bar values. As can be seen from the ellipse region, the blurry layers in the field seismic section can be distinguished very clearly in the reconstructed seismic section. This illustrates that the proposed algorithm effectively improves the vertical resolution of the seismic data. It is of great significance to the subsequent seismic interpretation.  We used the reflection coefficient section shown in Figure 10b and a high-frequency wavelet to perform a convolution operation, to obtain a reconstructed synthetic seismic section, which is shown in Figure 11a. Note that in order to keep the amplitude consistent with the original seismic data in Figure 10a, we calibrated the range of color bar values. As can be seen from the ellipse region, the blurry layers in the field seismic section can be distinguished very clearly in the reconstructed seismic section. This illustrates that the proposed algorithm effectively improves the vertical resolution of the seismic data. It is of great significance to the subsequent seismic interpretation.
Moreover, we compared the amplitude spectra of the field seismic section in Figure 10a and that of the reconstructed seismic section in Figure 11a, and plotted these amplitude spectra lines in Figure 11b. The high-frequency component of the reconstructed seismic section was significantly increased compared to the original field seismic section. This means that by using spectral inversion, we can enhance the high-frequency component of the reconstructed seismic section, and the high-frequency component helps to improve the resolution of the seismic section.
The experiments using the field seismic data demonstrate that by using the approach proposed in this paper, we can obtain a reflection coefficient section of high resolution and good lateral continuity. By combining our algorithm with seismic forward modeling, we can reconstruct a high-resolution seismic section so as to provide high-quality data for seismic interpretation.
in Figure 11a. Note that in order to keep the amplitude consistent with the original seismic data in Figure 10a, we calibrated the range of color bar values. As can be seen from the ellipse region, the blurry layers in the field seismic section can be distinguished very clearly in the reconstructed seismic section. This illustrates that the proposed algorithm effectively improves the vertical resolution of the seismic data. It is of great significance to the subsequent seismic interpretation.

Conclusions and Discussion
In this paper, we proposed a joint time-frequency reflection coefficient inversion method based on the basis pursuit de-noising algorithm. The experimental results show that if the length of the time window, regularization factor, dominate frequency band are optimized, we can obtain a reflection coefficient section of high quality. We evaluated and discussed the influences of these parameters on the inversion result, and found that the regularization factor and dominant frequency band have an important influence on inversion results. The obtained reflection coefficient section can then be used to reconstruct the synthetic seismic section. Compared with the original seismic section, the frequency bandwidth of the reconstructed seismic section was widened. Thus, the vertical resolution of the seismic data is effectively enhanced.
The advantage of our proposed method is that it makes full use of the time-frequency characteristics and sparsity of the signal to carry out the inversion of the reflection coefficient, which further reduces the multiple solutions of the inversion. Although our method obtains high-quality inversion results, the feasibility depends on the settings of the parameters. If these parameters are not optimized, it is not possible to obtain good results. The optimal setting of these inversion parameters should be studied further. Additionally, the reflection coefficient is related to physical properties of rock, so we can generalize our method to inverse other physical properties of rock. For example, seismic impedance can be represented by reflection coefficients through the convolution model, and we can inverse the impedance by using our method [20]. The velocity of the P-wave and S-wave can be obtained from the reflection coefficient by the Aki-Richards equation and we can extend our method to inverse velocity [43]. The Poisson's ratio and Young's modulus can be linearly represented by the Young s modulus, Poisson ratio and density (YPD) equation and we can obtain Young's modulus by our method and the YPD equation [44]. Moreover, we only utilized the sparsity of the reflection coefficient as prior information. In seismic inversion, there is other prior information that is favorable for inversion. The initial model can improve the stability and efficiency of the inversion [45]. The sparse information of the seismic record can improve the resolution of the inversion and enhance the anti-noise ability [46,47]. We can use this information to improve the resolution of the inversion. Finally, we constructed sparse constraints by using the L1 norm in this work. However, the L1 norm is not the best expression to describe the sparsity of gradient. The Lp quasi-norm has better sparsity than the L1 norm [48], therefore we can replace the L1 norm with the Lp norm in our method. Also, we only represented the vertical sparsity in our method. The vertical difference is sparse. We can replace the L1 norm to total variation and introduce the sparse information of the vertical difference [49]. Finally, we ignore the structured