Fracture detection from Azimuth-dependent seismic inversion in joint time–frequency domain

Detection of fracture properties can be implemented using azimuth-dependent seismic inversion for optimal model parameters in time or frequency domain. Considering the respective potentials for sensitivities of inversion resolution and anti-noise performance in time and frequency domain, we propose a more robust azimuth-dependent seismic inversion method to achieve fracture detection by combining the Bayesian inference and joint time–frequency-domain inversion theory. Both Cauchy Sparse and low-frequency constraint regularizations are introduced to reduce multi-solvability of model space and improve inversion reliability of model parameters. Synthetic data examples demonstrate that the frequency bandwidth of inversion result is almost the same for time, frequency and joint time–frequency domain inversion in seismic dominant frequency band using the noise-free data, but the frequency bandwidth in joint time–frequency domain is larger than that in time and frequency domains using low- signal-to-noise-ratio (SNR) data. The results of cross-correlation coefficients validate that the joint time–frequency-domain inversion retains both the excellent characteristics of high resolution in frequency-domain inversion and the advantage of strong anti-noise ability in time-domain inversion. A field data example further demonstrates that our proposed inversion approach in joint time–frequency domain may provide a more stable technique for fracture detection in fractured reservoirs.

Scientific Reports | (2021) 11:1269 | https://doi.org/10.1038/s41598-020-80021-w www.nature.com/scientificreports/ anisotropic parameters, a weak-anisotropy and linearized PP-wave reflection coefficient can be derived based on the seismic scattering theory and the first-order perturbations in stiffness components of rocks, building the bridge between the microscopic fracture parameters and macroscopic seismic reflection response [12][13][14][15] . As a result, the sensitive weaknesses can be estimated by combining the azimuth-dependent PP-wave seismic data and reflection coefficient equation. Amplitude versus offset and azimuth (AVOA) inversion has been an important method to predict the fracture information via the PP-wave azimuthal seismic data 16 . Gary et al. use the AVOA inversion to estimate the fracture density and fracture strike from the PP-wave seismic data 17 . Bachrach et al. also use the PP-wave seismic data to reconstruct the Thomsen-type anisotropic parameters based on the rock-physic-based Bayesian inversion in time domain 18 . Chen et al. estimate the sensitive fracture weakness parameters based on the difference in PP-wave azimuthal seismic data, and the proposed method is applied to a field data set 19 . Far et al. just use the synthetic data to estimate the sensitive fracture parameters, but they extend the inversion method to an arbitrary anisotropic medium 20 . In addition, Downton and Roure use the azimuthal Fourier coefficients to estimate the fracture weakness parameters 21 . The reflection-coefficient-based or Fourier-coefficient-based AVOA inversion mentioned above is generally performed in time domain. The inversion in time domain behaves better for the noisy seismic data and worse for the resolution of inversion results 22,23 . The seismic inversion in frequency domain has the advantages of high resolution of inversion results, but anti-noise performance is not good. Combining the timedomain seismic inversion, the seismic inversion in joint time-frequency domain can achieve a balance between inversion resolution and anti-noise performance 24,25 .
Moreover, the inverse problem of fracture estimation is ill-conditioned, and the inversion results will be unstable without any constraint to the problem. The estimation of sensitive fracture weaknesses should be implemented under the constraints of regularization terms 15 . Therefore, we attempt to use the PP-wave seismic reflected amplitude data to estimate the normal and shear weaknesses ( δ N and δ T ) of fractures with the regularization constraints in joint time-frequency domain. Integrating the relationship between Thomsen-type anisotropic parameters and fracture weakness parameters 26 , we first construct the forward modelling equation following Rüger's weak-anisotropy PP-wave reflection coefficient 12,13 . Then a Bayesian framework is introduced to the inversion for sensitive fracture weaknesses in joint time-frequency domain, which combines the prior constraint information and sparse-distribution likelihood function to estimate the posterior distribution of fracture parameters. In this paper, we construct the cost function using the assumption of a Cauchy-distribution prior constraint and a Gaussian-distribution likelihood function 15,27 . In addition, a low-frequency smoothing model constraint is also introduced to the cost function to obtain more stable estimation of Bayesian AVOA inversion 15,28 . We finally present a method of azimuth-dependent and azimuthal-seismic-amplitude-difference-based inversion to estimate the fracture parameters in joint time-frequency domain. The iteratively reweighted least-squares (IRLS) algorithm to solve the inversion problem for fracture estimation 29,30 . Synthetic data examples demonstrate that the normal and shear fracture weaknesses can be reasonably and reliably inverted when the PP-wave azimuthal seismic data contains moderate or even relatively high random noises. The real data set acquired over a fractured reservoir further validate that our proposed inversion approach in joint time-frequency domain can achieve the fracture detection from the azimuth-dependent PP-wave seismic data.

Methods
Forward matrix in time-frequency domain. Following Rüger's weak-anisotropy equation for an horizontal transversely isotropic (HTI) medium and the relationship between Thomsen's anisotropic parameters and fracture weaknesses 26 , the azimuth-related linearized PP-wave reflection coefficient R HTI PP (θ, ϕ) for an interface separating two HTI media can be written as in the form 15 , where θ and ϕ are the angles of incidence and azimuth, respectively; R ISO PP (θ) is the azimuth-independent background isotropic reflection coefficient, and R ANI PP (θ, ϕ) is the fracture-induced and azimuth-dependent reflection coefficient in an HTI medium formed by a single set of vertical and rotationally invariant fractures embedded in a homogeneously isotropic background rocks, which is given by where the symbol represents the value changes of normal and shear fracture weaknesses ( δ N and δ T ) between upper and lower layers separated by the reflection interface, and the weighting coefficients of fracture weaknesses ( a δ N and a δ T ) can be expressed as and Here g represents the square of S-to-P-wave velocity ratio of media. In the following paper, g used in the synthetic data examples is the ratio of the square of well log S-wave velocity and P-wave velocity, and g used in the field data example is the ratio of the square of initial S-wave velocity and P-wave velocity model.
To obtain the azimuth-dependent anisotropic parameters, we can just utilize the fracture-induced and azimuth-dependent reflection coefficient R ANI PP (θ, ϕ) to estimate the fracture weaknesses, which can be used for fracture detection. Integrating the estimated seismic wavelets, the vector of azimuth-dependent seismic reflection a δ N (θ, ϕ)= − g cos 2 ϕ sin 2 θ 1 − 2g 1 + sin 2 ϕ tan 2 θ + 1 − g cos 2 ϕ tan 2 θ , (4) a δ T (θ, ϕ)=g cos 2 ϕ sin 2 θ 1 − sin 2 ϕ tan 2 θ . where d t is the azimuth-dependent seismic difference data vector, A is the weight coefficient matrix of model parameters, m is the target model matrix, and G t = WA is the product of the seismic wavelet matrix and the weight coefficient matrix of model parameters, which can be expressed as is the wavelet matrix, and w j denotes the jth term of an extracted seismic wavelet; T is the matrix of reflection coefficient, respectively, and the symbol T represents the transposition of a matrix. In contrast, the seismic data d(ω) in frequency domain can be written as where ω is the angular frequency, W(ω) is the frequency spectrum of seismic wavelets, and R(ω) is the frequency spectrum of the fracture-induced and azimuth-dependent reflection coefficient R ANI PP (θ, ϕ) , which can be expressed as where τ (z) denotes the time-domain depth, and exp (·) is an exponential function. Equation (9) can be rewritten as in the form, where E(ω) represents the Fourier transform operator or the time shift operator. Bayesian inference in time-frequency domain. Bayesian inference in seismic inversion can be used to establish the a posteriori probability density function (PDF) as a product of the a priori PDF and the likelihood function 27 . The likelihood function depends on the PDF of background seismic noises. Assuming that the seismic data in time domain and data in frequency domain are both independent random variables 23,27 , the joint likelihood function in time-frequency domain can be expressed as where p(·) represents a PDF, and the three PDFs denote the degree of matching between inversion results and seismic data in joint time-frequency domain, time domain, and frequency domain, respectively. We further assume that the likelihood functions of seismic data d t and d f in time domain and frequency domain both satisfy the Gaussian PDF with mean zero, and the joint likelihood function can be expressed as where the symbol �·� 2 represents 2-norm function, σ 2 t and σ 2 f are the variances of time-domain and frequencydomain seismic data, respectively. Equation (13) links the seismic response between time-domain and frequencydomain data. The a priori PDF of unknown model parameters is used to describe the prior information of model parameters, and Cauchy distribution is utilized as the a priori PDF, which is given by Scientific Reports | (2021) 11:1269 | https://doi.org/10.1038/s41598-020-80021-w www.nature.com/scientificreports/ where σ 2 m is the variance of model parameter. Based on the Bayesian inference, the joint a posteriori PDF p m|d t , d f can be given by 28 that is, Maximizing the joint a posteriori PDF p m|d t , d f in Eq. (16), we obtain the objective function J(m) , which is given by where J Gauss (m) denotes the measurement of the difference between seismic response in joint time-frequency domain and forwarding synthesized gathers, and J Cauchy (m) denotes the sparse constraint regularization term introduced by the a priori PDF term; χ 1 =σ 2 t σ 2 f and χ 2 = 2σ 2 t are the regularization coefficients of seismic data error in frequency domain and sparse constraint term, respectively. Moreover, we introduce the low-frequencymodel constraint regularization term J mod (m) into the objective function in Eq. (17) which can be written as where χ 3 denotes regularization coefficient of low-frequency-model constraint, ζ and P are the low-frequency smoothing models of unknown model parameters and the integral matrix, respectively. Minimizing the final objective function J ALL (m) , we can get the nonlinear inversion equation. Here we use the iteratively reweighted least-squares (IRLS) algorithm to solve the nonlinear Eq. (18) iteratively 29,30 . After a couple of iterations, the IRLS algorithm can reach the state of convergence. When using the IRLS algorithm, some steps are demanded to implement the nonlinear and iterative inversion, including the construction of initial model parameters, the selection of iteration times, and the setting of convergence threshold. We then calculate the objective function iteratively, and finally obtain the inversion results according to the iteration times or the convergence threshold.

Results and discussions
To validate the proposed approach, we first use the synthetic data generated by PP-wave reflection coefficient (computed with Eq. 1) convoluted with seismic wavelets without noises, and then perform the azimuth-dependent and azimuthal-amplitude-difference-based seismic inversion for normal and shear weaknesses. Figure 1a shows the inversion results in time domain (blue curves), frequency domain (green curves), and joint time-frequency domain (red curves), respectively, and the initial model (dotted black curves) are generated by smoothing the true well log data (solid black curves). We find that the inverted fracture weaknesses are all consistent with the true values in all three domains. Figure 1b shows the corresponding spectra of differences in normal (above) and shear (below) fracture weaknesses, respectively, and the spectra of inverted results in all three domains are still a good match in the seismic band. Therefore, the inversion methods in all three domains perform well when seismic data contains no noises.
To further test the anti-noise ability of inversion methods in different domains, we add moderate Gaussian random white noises into the noise-free data and generate the synthetic data with different signal-to-noise-ratios (SNRs) being 5 and 2, respectively. We then perform the azimuth-dependent and azimuthal-amplitude-difference-based seismic inversion for normal and shear weaknesses in noisy cases. Figure 2a,b show the comparison between original and inverted model parameters and spectra in time domain (blue lines), frequency domain (green lines), and joint time-frequency domain (red lines) with synthetic azimuth-dependent seismic data containing moderate noises (i.e., the SNR of data is 5). From the inversion results in different domains shown in Fig. 2a, we can see that the inversion accuracy of fracture weaknesses in frequency and joint time-frequency domain are better that in time domain, and the inversion results in joint time-frequency domain are more stable than that in frequency domain from the spectra comparison shown in Fig. 2b. Figure 3a,b show the same case but with more noises (i.e., the SNR of data is 2), and we can find that the spectra of inverted fracture weaknesses in joint time-frequency domain are wider than the other inversion results in time or frequency domains. To validate the stability of the proposed inversion approach, we compare the cross-correlation coefficients between the true and inverted results in time, frequency and joint time-frequency domains. Table 1 illustrates the comparison results, and we can obviously find that the cross-correlation coefficients between true and inverted fracture weaknesses in joint time-frequency domain are larger than that in time and frequency domains when the seismic (14) p Cauchy (m) = 1 Scientific Reports | (2021) 11:1269 | https://doi.org/10.1038/s41598-020-80021-w www.nature.com/scientificreports/ data contains moderate or even more noises. Therefore, the inversion method in joint time-frequency domain maintains a balance between the anti-noise ability and resolution effect. A real data acquired from a fractured reservoir in Sichuan Basin, China is also used to further demonstrate the proposed inversion method, which is a wide-azimuth land survey dataset. To quickly estimate the normal and shear fracture weaknesses, we use two azimuth-dependent seismic data with three angles of incidence to implement the azimuth-dependent and azimuthal-amplitude-difference-based seismic inversion for fracture detection in joint time-frequency domain. The azimuth of fracture normal was first calculated using the leastsquares ellipse fitting method, and then the two azimuths were selected based on the estimated fracture normal to obtain the large seismic amplitude differences. Of course, the method can easily be extended to the multi azimuth data easily just changing the azimuth-dependent seismic difference data vector, the weight coefficient matrix of model parameters, and the wavelet matrix. However, the inversion with two azimuths is simple and generally  www.nature.com/scientificreports/ gives acceptable results in practice, and we attempt to implement the fracture detection using only two azimuth data to simplify the inversion processing. The trace spacing is 20 m. Before the seismic inversion, the data are processed to guarantee that the finely processed data is high-quality enough to be used for amplitude versus offset and azimuth (AVOA) inversion. The workflow of data processing is presented in Table 2, and details about the data processing are presented in Dulaijian 32 . Figures 4 and 5 are the azimuth-dependent seismic data with near, middle, and far angles of incidence generated from angle-stacked seismic data. The main frequencies of data vary from 10 to 50 Hz in this work area, and we use three inversion methods to estimate the fracture weaknesses in time, frequency, and joint time-frequency domains, respectively. All three methods are based on the information of seismic amplitude difference in different azimuth, and Fig. 6 shows the corresponding amplitude difference data. Note that the red curves in data profiles are the well log curves of shear weakness, which can be represented as the fracture density. Next we perform the proposed inversion method to characterize the fractured reservoirs.  The interpreted well-log fracture density information in this work area was not available. We estimated the fracture weaknesses based on an azimuthally anisotropic rock-physics model with conventional well log data 15 . Figure 7a,b are the initial models of normal and shear fracture weaknesses, and the curves are the corresponding  www.nature.com/scientificreports/ estimated information of fracture weaknesses. Of course, the estimated fracture weakness parameters should be calibrated by using the interpreted well-log fracture density information and the image logging, especially the micro-resistivity image logging (FMI). In this work area, there is no anisotropic well log information available, such as the interpreted well-log fracture density information, but the FMI information can be interpreted as the initial constraint of the estimation of the fracture development situation. Figure 8a,b are the inverted fracture weaknesses in time domain, Fig. 9a,b are the inverted fracture weaknesses in frequency domain, and Fig. 10a,b are the inverted fracture weaknesses in joint time-frequency domain, respectively. High-value fracture weaknesses illustrate the developed fractures in reservoirs. From the above inversion results of fracture weaknesses in different domains, we find that the inversion results in frequency domain or in joint-time-frequency domain show higher resolution compared with the inverted results in time domain, but the time-domain inversion results exhibit better lateral continuities. However, the joint-time-frequency-domain inversion results may provide more geologically reasonable interpretations for the fractured reservoirs in this area due to the discontinuous reservoirs of fracture development. Figure 11 is the comparison between inversion results in different domain at the well location, and we also find that the time-frequency-domain inversion method can achieve a balance between the anti-noise ability and seismic resolution. Figure 12 illustrates the comparison of histograms between original and inverted normal and shear weaknesses in time, frequency, and joint time-frequency domains, and we find that the a posteriori PDF of inverted fracture weaknesses in all three domains nearly agrees the Gaussian distribution with the a priori PDF as we have assumed.

Conclusions
Motivated by fracture detection in fractured reservoirs based on azimuth-dependent seismic inversion, we establish an inversion method by integrating Bayesian inference and regularization constraints in joint time-frequency domain to estimate the normal and shear fracture weaknesses. Combining the azimuth-dependent seismic amplitude difference information, we express the a posteriori probability distribution as a product of the a priori probability distribution and the likelihood function, and get the objective function by maximizing the a posteriori probability distribution. We finally estimate the characteristic parameters of fracture properties iteratively via a reweighted least-squares algorithm. Compared with the time-and frequency-domain inversion Table 1. Cross-correlation coefficients between true and inverted fracture weaknesses in time, frequency, and joint time-frequency domain.

Cross-correlation coefficients
Noise-free case SNR = 5 case SNR = 2 case  www.nature.com/scientificreports/ can obtain more accurate and robust inversion results than that in a separate time or frequency domain, in which the high-value fracture weaknesses are used to characterize the development areas of fractures. Therefore, the proposed inversion approach may provide a new way to perform the fracture detection by combining the different domain information for the seismic data.

Data availability
The datasets used during the current study are available from the corresponding author on reasonable request.

Appendix A
Following the assumption that a porous fractured medium can be taken as a periodic horizontally stratified layer, the normal and shear fracture compliances ( Z N and Z T ) can be defined as 31     License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.