Envelope-based constrained model for multiparameter pre-stack seismic inversion

Pre-stack seismic inversion is an effective method to estimate multiparameter from seismic reflection data in the oil–gas reservoir exploration. In order to get a reliable inverted parameter result, a constrained model is commonly used to reduce the multi-solution in the seismic inversion. However, it is hard to get an appropriate constrained model by conventional well-log data interpolation or velocity field in offshore oil–gas exploration with sparse wells. Envelope inversion has been successfully applied to get a good initial model with low-frequency components in full-waveform inversion. In order to fully exploit the rich information contained in seismic data, we introduce envelope inversion into background modeling based on convolution theory and propose a multiparameter pre-stack seismic inversion with envelope prior constraint. Both inversion results of the model and field data indicate the effectiveness and adaptability of the proposed method in multiparameter estimation.

As the theoretical basis of pre-stack seismic inversion, the relationship between seismic reflection wave amplitude and elastic parameters is characterized by the expression of the seismic reflection coefficient to quantitatively describe subsurface information. For the simplest assumption, the underground formation can be set as homogeneous isotropic media based on equivalent medium theory, and the seismic wave amplitude only varies with offset (AVO). And the exact form of reflection and transmission can be expressed as the Zoeppritz equation (Zoeppritz and Erdbebnenwellen, 1919). However, the nonlinearity and complex expression of the Zoeppritz equation limit its application to pre-stack seismic inversion. Corresponding simplification of the exact form with a well-posed assumption is the main approach to reducing the nonlinearity. Bortfeld (1961) proposed an approximation equation with the assumption of distinguishing the fluid and solid. Aki and Richards (1980) constructed an approximate linear equation of P-wave reflection coefficient by P-wave velocity, S-wave velocity, and density; the relationship between media parameters and reflection amplitude on the interface was directly described with linearized perturbation of corresponding parameters. Based on such linear approximation theory, several linearized analytical formulas were proposed for the expression between rock physical characteristics and elastic parameters (Shuey, 1985;Smith and Gidlow, 1987;Fatti, 1994;Goodway et al., 1997;Gray et al., 1999;Russell et al., 2011;Zong et al., 2012Zong et al., , 2013Zong et al., , 2015Yin et al., 2014Yin et al., , 2015Li et al., 2020;Wu et al., 2022).
Pre-stack seismic inversion infers the subsurface parameters by mapping the connection of seismic data and inverted model parameters. According to the analytical expression of forward modeling, the objective function for seismic inversion is deduced to calculate relevant parameters that are applied in reservoir prediction and fluid identification. Unfortunately, pre-stack inversion is an ill-posed problem that always results in multiple solutions while affecting the accuracy of reservoir characterization. To reduce the multisolutions, the regularization term is introduced in seismic inversion as a constraint, which aims to enhance the stability of the inverse equation (Tarantola, 1987;Downton, 2005;Theune et al, 2010;Karimi et al, 2010;Alemie and Sacchi, 2011;Zhang et al, 2013;Yin et al, 2016;Pérez et al, 2017).
Another factor that should be considered when introducing the constrained term in seismic inversion is the quality of observed data for the inverse application. Reliable inverted results also depend on the seismic data quality, and the corresponding low-frequency component affects the stability for seismic inversion (Zhang et al., 2015), which is always missed in actual observed data. Therefore, the solution of to this problem is how to build a good initial model that contains low-frequency information to constrain the inversion process. For the traditional initial model construction, logging data interpolation (Doyen et al., 1997;Liu et al., 2008) or velocity information (Al-Yahya and Kamal, 1989;Li et al., 2014;Zhao et al., 2019a, b) is usually used for modeling, while the inversion quality affects by well number or velocity field's information. To solve the problem of low-frequency missing, a seismic envelope inversion based on signal modulation theory is proposed for the initial model construction in full-waveform inversion (Wu et al., 2014;Luo et al., 2015;Chen et al., 2018;Luo et al., 2019Luo et al., , 2021Yang et al., 2022).
In this paper, we propose a constrained multiparameter seismic inversion method by set the envelope information as the initial model. Following the medium decomposition theory, we decompose the inverted parameters as background terms and perturbation terms. To inverse the information of background terms, we deduce the envelope form based on the convolution model, and the corresponding multiparameter inversion method is constructed via the envelope expression. In the inverse process of perturbation term, we update the former envelope inversion result as the model constraint to calculate the reliable corresponding parameters. To demonstrate the feasibility and stability of the proposed method, we discuss the results obtained from synthetic tests and an actual exploration case. The test results show that the proposed inversion approach in this study is easily extended to more complex geologic conditions because there is no need to implement interpolation or extrapolation processing to estimate the low-frequency components of the elastic parameters.

Principle of modulation and demodulation for seismic signal
The principle of modulation and demodulation contains two aspects: signal modulation and demodulation. For the signal 1 3 modulation, a low-frequency baseband signal loads into the high-frequency carrier wave, and a relevant spectrum shift can be realized through signal modulation. As the former modulated signal is the product of the baseband signal and carrier signal, the modulated result can be equivalent to resampling the baseband signal, while the corresponding frequency relates to the carrier one. For signal demodulation, it can be seen as the reverse process of signal modulation.
Since we set the signal modulation into the seismic wave propagation, the source information can be approached as the modulated signal, while the seismic wave is equivalent to the carrier one. Therefore, seismic wave propagation can be expressed as follows (Wu et al, 2014): with where u x s , x r ;t is the seismic wavefield transmitted from x s to x r , G sr x s , x r ;x i , t is the Green function that represents the modulated signal, r (t) is the convolution result of source wavelet s(t) and reflection coefficient sequence r i x i , t that considered as carrier signal. And the low-frequency baseband signal E x s , x r ;t can be extracted by demodulating received seismic signal u x s , x r ;t , that is, the envelope of the seismic data. Figure 1 plots the process of signal modulation and demodulation. As the figure shows, Fig. 1a is the modulation signal (seismic signal envelope) with low-frequency, Fig. 1b is the carrier signal (seismic wave) with high-frequency, Fig. 1c is the product of the modulation signal ( Fig. 1a) and carrier signal (Fig. 1b) modulated signal called modulated signal (Wu et al., 2014), and Fig. 1d is the envelope signal obtained by demodulation. The relevant spectrums for these signals show in Fig. 1e-h, respectively. By the spectrum analysis, we can see that the frequency band of the modulated signal is consistent with the carrier signal one; the modulation frequency does not evolve visible on the modulated signal spectrum (Fig. 1f). Therefore, the low-frequency information can be extracted by envelope calculation. Compared with the modulation signal (Fig. 1e) and demodulation envelope (Fig. 1h), the demodulation envelope spectrum is consistent with modulation signal one, which illustrates the reasonability of supplementing the low-frequency information by the modulation envelope.
According to the principle of signal demodulation, the seismic record d(t) can be demodulated by Hilbert transform to obtain an analytical signal d (t) with non-negative frequency (Wu et al, 2014), which is given as follows: Signal modulation and demodulation. a modulation signal; b carrier signal; c modulated signal and its envelope; d demodulation envelope; e modulation signal spectrum; f carrier signal spectrum; g modulated signal spectrum; h demodulation envelope spectrum where the real part of the analytical signal is seismic data, and the corresponding imaginary part is transformed seismic data; then, the envelope E(t) can be expressed as follows (Wu et al, 2014): With the assumption of convolution theory, the relationship of source wavelet s(t) , reflection coefficient r(t) , and seismic record d(t) in a linear time-invariant system is as follows (Robinson 1957): Setting Eqs. (4) into (3), the new form of envelope E(t) can be characterized as follows: Using the properties of convolution calculation and Hilbert transform, Eq. (5) can be rewritten in the following form: As Eq. (6) shows, the envelope expression for a plane seismic wave is deduced. Figure 2 shows the envelope information of seismic, and synthetic data, in which the black curve is the seismic record, the red line is the seismic envelope by signal demodulation, and the blue dash line is the envelope calculated by Eq. (6). As the figure shows, the calculation result via Eq. (6) is consistent with that directly extracted in seismic data, which proves the reliability of the modulation envelope based on the convolution theory.
Due to the complexity of plane wave envelope expression, the expression of reflection coefficient with envelope information cannot be simply decomposed from the convolution model. If we assume envelope seismic is satisfied with Eq. (4), the quasi-reflection coefficient operator for convolution only represents the reflection coefficient's low-frequency information. Although obtaining the exact expression of the envelope reflection coefficient is difficult, the relevant low-frequency component is useful to update the inverted parameters in the inversion process.

Envelope-based constrained seismic inversion method
Under the Bayesian inversion framework (Zong et al., 2018), the inverse problem is expressed in the form of the probability distribution. As a constraint of the objective function in the inversion process, the background model is used to reduce the degree of freedom for model parameter estimation, which is widely used in the geophysical inversion. In this section, we propose an envelope-based constrained multiparameter seismic inversion method under the combination of envelope inversion and Bayesian inversion framework. Assuming that parameters for inversion are P-wave velocity, S-wave velocity, and density, thus the corresponding reflection operator m can be represented by reflectivity for P-wave velocity m P , S-wave velocity m S , and density m D , m = m P , m S , m D (Aki and Richards 1980), and relevant observed data set can be written as follows (Tarantola 2005): where d(t i , i ) represents the observed seismic record with incident angle i at time t i , and the time series under the same incident angle i can be written as d( i ). If we take the observed seismic data as a modulated signal set, the seismic envelope data set can be characterized as follows: According to the convolution theory, the seismic envelope can be expressed in a form similar to Eq. (4). Assuming that quasi-envelope reflection sequence R E consists of the linear perturbation for P-wave velocity, S-wave velocity, and density (Aki and Richards 1980). The matrix of quasi-envelope reflection sequence R E with incident angles can be shown as follows: where R E i is quasi-reflection sequence with incident angle i , A( i , B( i and C( i represent the coefficient terms of P-wave velocity, S-wave velocity, and density, respectively, while the corresponding envelope forms for the parameters are characterized as m pE , m sE and m dE . With the combination of wavelet matrix e and coefficient matrix, the forward operator is rewritten as follows (Buland and Omre 2003): And the expression of the Robinson convolution model for the seismic envelope is deduced as follows: where G e is the forward operator, e is the parameter for inversion, and is random noise series.
Under the Bayesian inversion framework (Buland and Omre 2003;Alemie and Sacchi 2011;Zong et al. 2012), the posterior probability density function (PPDF) in the data-model parameter space can be characterized by prior distribution p e and likelihood function p | e , which is given as follows: where p( ) represents the seismic data acquisition probability, independent of model parameters. This term is a constant term for the inversion problem with known seismic data to ensure that the posterior probability is integrated into 1 in the model space.
Assuming that the noise of the observed seismic data obeys the Gaussian distribution, and prior distribution follows the Cauchy distribution, the PPDF can be then expressed as follows (Tarantola 2005): where 2 n represents the noise variance in observed data, 2 m is the variance of model parameters, and m ei is the reflection coefficient to be inverted.
Introducing the smooth background model as the constraint for the Bayesian envelope inversion, the fitting error between model parameters and smooth model follows Gaussian distribution. Accordingly, the new form for prior distribution contains the joint information of model parameter and smooth background model. The PPDF then can be rewritten as follows (Yin et al. 2008;Zong et al. 2012): where is the integral matrix, represents the background model data, and 2 is the variance between background model and inversion parameters. Characterizing Eq. (13) in logarithmic domain, the objective function is obtained as follows: where is the Cauchy constraint weight, i is the model constraint weight for multiparameter, respectively, and the accuracy of inverted parameters is controlled by relevant weight.
Taking the envelope inversion results as model input , the corresponding objective function for envelope-based constrained inversion can be written as follows: To maximize the posterior probability of obtaining model parameters, we deduce the first derivative form of Eq. (16), and the final objective equation is as follows (Yin et al., 2008;Zong et al. 2012): with The multiparameter results can be calculated by solving Eq. (17) with the iteratively reweighted least square (IRLS) method. Since the envelope inversion provides a more accurate low-frequency trend, the pre-stack seismic inversion based on envelope constraint can obtain more precise parameter results. (14)

Model test
To demonstrate the feasibility and accuracy of the proposed method, the SEG/EAGE overthrust model is used for the model test. The grid size of the model is 186 × 801 , and a Ricker wavelet with a central frequency of 25 Hz is used to calculate synthetic data (plot in Fig. 3). And model parameters are shown in Fig. 4a, 5a, 6a.
For the model test of envelope pre-stack seismic inversion, the input seismic data are the corresponding seismic record envelope information. We first extract the seismic envelope from the seismic data with the incident angle of 10 • ,20 • and 30 • , and the linear model is introduced as a background constraint (Fig. 4b, 5b, 6b). The envelope parameters inverted by the former algorithm are shown in Figs. 4c, 5c, 6c, and a single trace's corresponding multiparameter results are plotted in Fig. 7. It can be seen that the parameter results of envelope inversion mainly update the background information of the overthrust model, which clearly restored main river channel and background structure as the slice shown (Fig. 8b, e and h).
For the model test of envelope-constrained pre-stack seismic inversion, we then set the envelope inversion results as model constraints, and seismic data replace the input data of envelope information for former inversion with three angles. The final inverted parameters are shown in Fig. 4d, 5d, 6d,   Fig. 4, 5, 6 shows that the inverted results are consistent with real model parameters, and the envelope-constrained inversion effectively updates relevant detail information (Fig. 8c, f, i).
To compare the parameter estimation accuracy of this approach and the conventional inversion methods, we calculate the relative error between inversion results and model data, respectively. Figure 9a, 10a. Figure 10a shows the inversion results of conventional parameters, and the relative error is shown in Fig. 9c, 10c; Fig. 9b, 10b plots the inversion results of conventional parameters, and the relative error is shown in Fig. 9d, 10d. We can see that since the low-frequency information of the initial model is updated using envelope inversion, the envelope-constrained inversion results are more accurate than conventional inversion (Fig. 11).
To test the stability of the proposed method, we add random noise into synthetic seismic with a signal-to-noise ratio equal to 3:1. Figure 12 shows the multiparameter inversion results, while the comparison of the single trace is shown in Fig. 13. We can see that the inversion based on proposed method can still obtain available parameter results. However, since the linear approximation is chosen to construct the objective function for our inversion, the accuracy of inverted density is lower than velocities. If the objective function is deduced by nonlinear approximation or exact equation, the accuracy of corresponding inverted parameters could be further improved.

Field data application
A set of seismic data is used to test the feasibility of the proposed multiparameter estimation approach in the field data. There is only one well in the test field, and the partial angle stack data sets of three angles ( 8 • ,16 • and 24 • ) are obtained by angle gather.
We first calculate the seismic envelope by Hilbert transform, the mid-angle stack data and its envelope information which are shown in Fig. 14. For the procedure of envelope inversion, we set envelope information as input seismic data, and the multiparameter inversion results are shown in Fig. 15a, b, c. As the results shown, the inverted parameters supplement the large set of background information.
For the procedure of constrained inversion, relevant inverted parameters with envelope information are then updated as model constraints to following seismic inversion, while the final inverted parameters are shown in Fig. 15d, e, f. Contrast to envelope results, the final results mainly update the detail characteristic based on envelope constraints. For the comparison with well-log curves and relevant inverted parameters at well location, we can see that the final inversion results effectively characterize the parameter variation, illustrating the good inversion accuracy of the proposed method.
To further confirm the adaptability of the proposed approach, we select another line for the field data application. The partial angle stack data sets are same as the first field test. As Fig. 16 shows, the horizontal characteristics of seismic amplitude data change in a wide range, while the corresponding signal-to-noise ratios vary in different regions. There are two wells in the test field, but the far-spaced location and incomplete well-log data make it challenging to construct a reliable model for seismic inversion.
We use the conventional and proposed methods to estimate the elastic parameters under the same initial model. The inversion results are shown in Fig. 17. Since envelope inversion results provide a more accurate ultralow-frequency trend and abundant low-frequency information reduces the multi-solution in the inversion process, the envelope-based constrained inversion can improve the parameter estimation accuracy. The inversion results indicate the applicability and calculation accuracy of the proposed method.

Conclusions
A multiparameter pre-stack seismic inversion approach based on envelope constraint is proposed in this study. By updating the data input and model constraint, the proposed method can obtain more precise elastic parameters inversion results rich in low-frequency under the Bayesian theory framework. The conclusions can be drawn as follows: 1. Compared with the classical pre-stack inversion method, this study estimates the low-frequency components of the initial model by envelope inversion. It uses the updated model as the model constraint input to improve the accuracy of inverted parameters. 2. The main advantage of this approach is that it reduces the dependence on the well-log data, effectively enhances the stability of elastic parameter inversion in the area with sparse wells, and decreases the multisolution of inversion results. 3. In the proposed method, only the first-order perturbations of parameters consist in the expression of the forward modeling equation, which may result in parameter estimation with low resolution in the case of a complex medium. 4. In future research, this study can be improved mainly from both forward modeling and inverse process aspects. For the expression of forwarding equations, making full use of high-order perturbation information is the key to improving the representation accuracy of parameters, and deducing forward operators with high-order perturbation can improve the representation accuracy of equations. On the other hand, for the inversion algorithm, a nonlinear inversion algorithm is needed to match the high-order disturbance equation, which is an effective