Guang-Sen Cheng ,Xing-Yao Yin ,Zhao-Yun Zong ,*,Ya-Ming Yang

a School of Geosciences,China University of Petroleum(East China),Qingdao,Shandong 266580,China

b Pilot National Laboratory for Marine Science and Technology(Qingdao),Qingdao,Shandong 266580,China

c Shandong Provincial Key Laboratory of Deep Oil and Gas,Qingdao,Shandong 266580,China

Keywords:

ABSTRACT

1.Introduction

Linear(Aki and Richards,1980;Shuey,1985;Goodway et al.,1997;Zong et al.,2015;Li et al.,2020),nonlinear(Wang,1999;Stovas and Ursin,2003;Yin et al.,2013b;Cheng et al.,2018;Liu et al.,2020;Zhou et al.,2020),and exact plane-wave reflection coefficients(Ursin and Tj?land,1996;Pan et al.,2017;Yin et al.,2018;Zhou et al.,2021)have taken a vital role in the reservoir prediction of pre-stack seismic exploration(Zong et al.,2012;Yin et al.2013a,2014;Li et al.2017,2020).The plane-wave reflection coefficient(PRC)is derived from the planar-wavefront assumption,which is an approximation of spherical-wave reflection coefficient(SRC)in far field.At the pre-critical incident angle,PRC is always real value and there is no phase shift.However,the critical angle is often reached in field data acquisition,such as in the case of strong properties contrast,salt body,carbonate rock and so on.At the critical and post-critical incident angles,PRC becomes inapplicable(O'Brien,1963).

SRC can be expressed as the integral of PRC(Aki and Richards,1980;Haase,2004;Ursenbach et al.,2007;Skopintseva et al.,2011),which is more accurate than PRC to describe the seismic reflection wave excited by point source,especially at the critical and post-critical incident angles.Due to the attenuation of viscoelastic medium,the complex-valued PRC(Innanen,2011;Bird,2012;Zong et al.,2015)is also generated.No matter what the incident angle is,SRC is always a complex value,which includes the spherical-wave amplitude and phase reflection information.Compared with the conventional AVO inversion,phase variation with offset/angle(PVO/PVA)inversion mainly uses the phase-shift information of reflected seismic waves,which provides the potential of accurate density estimation(Zhu and McMechan,2012).However,the complexity of the calculation of SRC integral formula brings difficulties to practical application.The use of imaginary part of SRC further exacerbates this problem.How to make full use of the amplitude and phase information of complex-valued SRC has become an important subject.

Complex seismic traces have been used to estimate seismic attributes(Barnes,2007),such as delineating thin lenses in seismic sections(Robertson and Nogami,1984).Complex seismic traces are also used for AVO inversion.Zong et al.(2015)utilized the complexvalued PRC to estimate P-and S-wave quality factors simultaneously.The frequency-dependent complex-valued SRC in acoustic media was also analyzed,and the spherical-wave AVO inversion of single-reflection was implemented in the synthetic data example(Li et al.,2017).Cheng et al.(2020)applied the SRC to the real seismic traces near borehole,but only the real part of seismic data is used.

We first investigate the amplitude and phase characteristics of SRC,and obtain the complex-valued spherical-wave seismic traces by using a simple harmonic superposition model.It is assumed that geophone can only record the real part of complex seismic traces.To utilize the complex-valued SRC,a novel complex spherical-wave EI inversion approach is proposed and the elastic parameters are further extracted.Our inversion approach is split into two steps:the Bayesian framework is used to estimate the complex sphericalwave elastic impedance(EI)from seismic data with different frequency components and incident angles,and extracting elastic parameters from the complex-valued EI.EI is a generalization of acoustic impedance(AI)and was first put forward by Connolly(1999).Whitcombe(2002)further normalized the EI in terms of P-and S-wave velocities and density.Subsequently,EI is widely used for elastic parameters estimation and reservoir prediction(Ma,2003;Martins,2006;Yin et al.,2013a;Zong et al.,2013;Su et al.,2014;Chen et al.,2018;Cheng et al.,2019).In the complex spherical-wave EI inversion section,the Bayesian scheme(Buland and More,2003)is used to estimate the real and imaginary parts of complex-valued EI simultaneously.We consider that the prior probability obeys the Cauchy probability distribution and the likelihood function obeys the Gaussian probability distribution(Zong et al.,2017;Chen et al.,2018).To extract the elastic parameters,the complex spherical-wave EI equation is derived based on the SRC.The accuracy of the EI equation is basically consistent with that of the exact complex-valued SRC.The P-and S-wave velocities and density are further obtained by combining the complex EI equation and the inverted complex spherical-wave EI.Synthetic and field data examples show that our approach is valid,and the inversion results of complex spherical-wave EI and velocities are in good agreement with the corresponding true value.

2.Spherical-wave forward modeling

2.1.Amplitude and phase of complex spherical-wave reflection coefficient

In our study,only the PP-wave reflection is considered.The exact three-parameter SRC(Cheng et al.,2020)that has compensated for the geometrical spreading can be written ascosθpi,where Rppis the exact three-parameter PRC(Yin et al.,2018),which is the function of P-wave velocity reflectivityS-wave velocity reflectivitydensity reflectivityand the P-wave incident angleθpi.The P-and S-wave velocities and density in the upper medium are denoted vp1,vs1andρ1,and in the lower medium are denoted vp2,vs2,andρ2.vp,vs,andρa(bǔ)re the corresponding average.ωis the angular frequency of harmonic wave.i is the imaginary unit and x is the integral variable.J0and J1are the zero-and first-order Bessel function.h and z are the vertical distances from the reflected interface to the source and geophone,respectively.r is the horizontal offset.To calculate the complicated integral equation(1),a powerful and stable algorithm,the adaptive Gauss-Kronrod quadrature(Shampine,2008),is used to solve the integrand function.In the synthetic data example,vp1is assumed to be known.In field data example,vp1can be obtained by tomographic velocity(Zhu and McMechan,2012).

At the angular frequency ofωn,the SRC in equation(1)at different reflector depths can be expressed by a time-continuous function and written as

where SRCR(t,ωn)and SRCI(t,ωn)are the real and imaginary parts of SRC(t,ωn),respectively.

The elastic parameters of Model 1 are shown in Table 1(Haase,2004),which are used to calculate the complex-valued SRC(equation(1))at the frequency of 15 Hz and the reflected interface depth of 600 m.The PRC is given for comparison.Fig.1 displays the comparisons of amplitude and phase between SRC and PRC.The black and red solid curves denote the SRC and PRC,respectively.We observe that the amplitude and phase of PRC is discontinuous at the critical incident angle,and there is no phase change in PRC before the critical angle.Unlike the characteristics of PRC,the amplitude and phase of SRC are smooth around the critical angle,and there is always phase shifts of SRC at the critical,pre-and post-critical incident angles.Fig.2 displays the comparisons of real and imaginary parts of SRC(black solid line)and PRC(red solid line).From Fig.2 we can see that the imaginary part of PRC is always zero value at the pre-critical incident angles,and the amplitude and real part of SRC are basically consistent with PRC at the small pre-critical incident angles.However,the imaginary part of SRC cannot be ignored.

Table1 Model 1.

Fig.1.Comparisons of(a)amplitude and(b)phase between SRC(black solid curves)and PRC(red solid curves).

Fig.2.Comparisons of(a)real and(b)imaginary parts of SRC(black solid curves)and PRC(red solid curves).

2.2.Modeling the complex-valued spherical-wave synthetic seismogram

To fully exploit the amplitude and phase reflection information,it is necessary to use the complex-valued SRC to construct the spherical-wave synthetic seismogram.Seismic wavelet can be decomposed into harmonic waves with different amplitudes,phases,and frequencies(Krebes,2019),which are usually described by complex exponentials and written as

where|An|andφnare the amplitude and phase of harmonic wave.t is the time.ωmin≤ωn≤ωmax,andωmax-ωminis the bandwidth of seismic wavelet.The real part of equation(3)is

which is the frequency component of the real part of wavelet and describes the physical properties of harmonic wave(Krebes,2019).The imaginary part of equation(3)is

which is the frequency component of the imaginary part of wavelet.Since ei·θ=cosθ+i·sinθ,i=√————1,equation(3)can be written as

We use the convolutional model(Robinson,1985)to convolute equation(6)with equation(2),the spherical-wave syntheticseismogram at the angular frequency ofωncan be obtained and expressed as

Stratumvp,m/svs,m/sρ,kg/m3 Upper medium289812902.425 Lower medium285716662.275

where asterisk denotes the convolution operation.Based on the superposition principle of wave,the synthetic seismogram of fullfrequency band is obtained by using the simple harmonic superposition model

equation(8)is further written as

where

and

Fig.3.Comparisons between equations(1)and(27)calculated by Model 2.(a)real and(c)imaginary parts of equation(1),(b)real and(d)imaginary parts of equation(27),(e)the relative error between(a)and(c),and(f)the relative error between(b)and(d).

where H[·]denotes the Hilbert transform operator.From equation(10)and equation(11)we can see that the imaginary part of wavelet can be got from the real part based on Hilbert transform.Substituting equations(10)and(11)into dr(t),we can obtain

Equation(9)describes the relation between the complex seismic data and SRC at the angular frequency ofωn,and its matrix form at the incident angleθpiis

Fig.5.Comparisons between true values(red solid curves)and inverted complex spherical-wave EI(blue dotted curves)at different frequencies and incident angles without noise.Green curves denote the initial model constraints.

Fig.6.Comparisons between true values(red solid curves)and inverted plane-wave EI(blue solid curves)at different frequencies and incident angles without noise.Green curves denote the initial model constraints.

3.Elastic inversion of complex spherical wave

3.1.Complex EI inversion with the Bayesian scheme

Connolly(1999)initially proposed the relation between PRC and elastic impedance

where EI1and EI2denote the elastic impedance in upper medium and lower medium,respectively.EI is the corresponding average andΔEI=EI2-EI1.The plane-wave elastic impedance is frequency independent.To make full use of the real and imaginary parts of complex seismic signal,we proposed the concept of spherical-wave elastic impedance,which is also expressed by EI in this paper

where EIR and EII are the real and imaginary parts of the complex spherical-wave EI,respectively.Similarly,

and

whereΔln(EIR)=ln(EIR2)-ln(EIR1)andΔln(EII)=ln(EII2)-ln(EII1),ln(-)is the natural logarithm.The subscript 1 and 2 denote the medium 1 and medium 2,respectively.After integrating equations(16)and(17),we can further obtain

and

where t0and t are the start and end travel time of seismic data.To estimate the complex spherical-wave EI at the frequency component of fn(fn=)and the incident angle ofθpi,the forward solver dωn=WωnRωnis used.

The complex spherical-wave EI inversion is implemented in a Bayesian framework(Buland and More,2003).The posterior probability density function p(Rωn|dωn)is

Fig.7.Comparisons between true values(red solid curves)and inverted complex spherical-wave EI(blue dotted curves)at different frequencies and incident angles.Green curves denote the initial model constraints and the signal to noise ratios(S/N)is 5:1.

Fig.8.Comparisons between true values(red solid curves)and inverted plane-wave EI(blue solid curves)at different frequencies and incident angles.Green curves denote the initial model constraints and the signal to noise ratios(SNR)is 5:1.

Fig.9.Comparisons between true values(red solid curves)and inverted complex spherical-wave EI(blue solid curves)at different frequencies and incident angles.Green curves denote the initial model constraints and the signal to noise ratios(S/N)is 2:1.

To improve the inversion resolution,the prior probability p(Rωn)is considered to obey the Cauchy probability distribution(Alemie and Sacchi,2011)and written as

The likelihood function p(dωn|Rωn)obeys the Gaussian probability distribution and written as

Maximizing the posterior distribution function(Zong et al.,2017),the objective function can be obtained and written as whereinitial model constraint of the real and imaginary parts of complex spherical-wave EI.λrandλirespectively are the corresponding constraint coefficients of the real and imaginary parts.The introduction of the initial model constraint can make the inversion results stable and without serious distortion.The larger the constraint coefficients are,the closer the inverted elastic impedance are to the initial model value.Iteratively Re-weighted Least Squares(Daubechies et al.,2010)is used to solve equation(24),the complex-valued spherical-wave EI of different frequencies are estimated from the observed seismic data(dω1,dω2???)at the corresponding frequency components.

Fig.10.Comparisons between true values(red solid curves)and inverted plane-wave EI(blue solid curves)at different frequencies and incident angles.Green curves denote the initial model constraints and the signal to noise ratios(S/N)is 2:1.

Fig.11.Comparisons of P-wave velocity,S-wave velocity,and density between the spherical-wave inversion results(black solid curves),plane-wave inversion results(blue solid curves),and true values(red solid curves)without noise.Green curves denote the initial model constraints.

3.2.Elastic parameters extraction from spherical-wave EI with different frequencies

The complex EI is estimated firstly,and then the elastic parameters can be further extracted from the complex EI.The complex EI equations are derived based on equations(16)and(17)and written as

Fig.12.Comparisons of P-wave velocity,S-wave velocity,and density between the spherical-wave inversion results(black solid curves),plane-wave inversion results(blue solid curves),and true values(red solid curves)in the case of noise.Green curves denote the initial model constraints and the SNR is 5.

Equations(25)and(26)establish the relation between the complex spherical-wave EI andθpi,andωn,which are utilized to estimate the elastic parameters(P-and S-wave velocities,and density)from the inverted complex spherical-wave EI.To verify the accuracy of equations(25)and(26),Model 2(Goodway et al.,1997)and Model 3(Ostrander,1984)as shown in Table 2 and Table 3 are used.A novel SRC is derived and expressed as

Table3 Model 3.

We implement the comparisons between equations(1)and(27)at different frequencies and incident angles using Model 2 and Model 3,as shown in Figs.3 and 4 respectively.Fig.3(a)and(b)displays the real and imaginary parts of exact SRC(equation(1))using Model 2.The real(Fig.3(c))and imaginary(Fig.3(d))parts of novel SRC(equation(27))are given for comparison.The relative error of real and imaginary parts between equations(1)and(27)are displayed in Fig.3(e)and(f),respectively.Fig.4 displays the accuracy comparisons using Model 3.From Figs.3 and 4 we can see that the novel SRC is in good agreement with the exact SRC.

The inverted complex spherical-wave EI with different frequency components and incident angles are preserved as the observed datasets(EIobs)to extract the elastic parameter vector m,Given a model vector m,the spherical-wave EI at different frequencies and incident angles can be obtained and expressed as EImod.With equations(25)and(26),the parameterextraction objective function can be constructed and expressed as

Fig.13.Comparisons of P-wave velocity,S-wave velocity,and density between the spherical-wave inversion results(black solid curves),plane-wave inversion results(blue solid curves),and true values(red solid curves)in the case of noise.Green curves denote the initial model constraints and the SNR is 2.

Fig.14.Workflow of complex spherical-wave elastic inversion for P-and S-wave velocities and density based on the complex spherical-wave elastic impedance.

4.Synthetic data examples

A well log is utilized to test the feasibility and stability of our inversion approach.Given f1=10 Hz,f2=20 Hz and f3=30 Hz,the corresponding SRC can be computed using equation(1).Given the real and imaginary parts of Ricker wavelets with the dominant frequencies of 10 Hz,20 Hz,and 30 Hz,the corresponding complexvalued spherical-wave synthetic seismograms can be obtained.The Gaussian random noises are also added into the synthetic seismograms to test the stability of the complex spherical-wave inversion approach.The signal to noise ratios(S/N)are 5:1 and 2:1.The complex-valued spherical-wave synthetic seismograms with the incident angle ranges of 8?-16?,17?-25?and 26?-34?are respectively stacked to obtain the small-angle,middle-angle,and large-angle observed seismic data,and their dominant incident angles areθp1=12?,θp2=21?,andθp3=30?respectively.To show the advantages of using SRC instead of PRC,the conventional planewave EI inversion is also implemented by using the spherical-wave synthetic seismograms with different frequencies as the observed data.

Figs.5 and 6 display the complex spherical-wave EI inversion results and plane-wave EI inversion results when there is no noise.In Figs.5 and 6,it can be observed that the complex spherical-wave EI inversion results are good agreement with the true EI value,and the differences between the plane-wave EI inversion results and true EI value decrease with the increase of frequency and propagation distance,which indicates that the SRC cannot be ignored in the case of near field and the spherical-wave effect is not obvious in far field.Figs.7-10 display the complex spherical-wave EI inversion results and plane-wave EI inversion results in noise situation.In Figs.7 and 8,SNR is 5.In Figs.9 and 10,SNR is 2.Considering that the weight of the imaginary part to the spherical-wave reflection coefficient is less than the real part(Fig.2),the imaginary part is more sensitive to noise than the real part.In Figs.7-10,we can observe that the real parts of spherical-wave EI inversion results are better than the imaginary parts and plane-wave EI inversion results.Even in the case of noise,the real parts of spherical-wave EI and plane-wave EI can be estimated stably.Similarly,the differences between the plane-wave EI inversion results and true EI value decrease with the increase of frequency and propagation distance.However,due to the influence of noise,the spherical wave effect become weaker.

Fig.15.The(a)real part and(b)imaginary part inversion results of the complex spherical-wave EI at the incident angles and frequencies of 12?and 30 Hz,12?and 40 Hz,21?and 30 Hz,and 21?and 40 Hz,respectively.

Fig.16.The plane-wave EI inversion results at the incident angles of(a)12?,(b)21?,and(c)30?,respectively.

With the inverted complex-valued spherical-wave EI of different incident angles and frequencies,we can further extract the P-and S-wave velocities and density using equation(28).The inverted plane-wave EI and the exact Zoeppritz equation of PPwave are also used for plane-wave elastic parameters estimation.One of the advantages of our spherical-wave inversion approach is that SRC is frequency dependent and we can use more reflection information of seismic data with different frequencies and incident angles.In the far field,the available incident angle is very small.The more seismic reflection information can be provided from the seismic data with different frequency components.When the seismic frequency band is narrow,the more seismic reflection information can be provided from the seismic data with different incident angles/offsets.

Figs.11-13 display the comparisons between true values and the inversion results of P-wave velocity,S-wave velocity,and density when there is no noise,SNR=5,and SNR=3,respectively.The red,black,blue,and green curves denote the true values,spherical-wave inversion results,plane-wave inversion results,and the initial model,respectively.From Fig.11 we can see that the inversion results of P-and S-wave velocities,and density show a good agreement with the true values when there is no noise.Compared with the inverted elastic parameters using our sphericalwave inversion approach,the difference between the true values and the plane-wave inversion results is larger.With the increase of noise,the differences between the inversion results and true values increases.Under the influence of noise,the spherical wave effect is no longer obvious.The P-and S-wave velocities and density estimated from the spherical-wave inversion approach are better than those obtained by plane-wave approach.Especially for the density term,the advantage of our approach is more obvious.Even in the case of noise,no matter the signal-to-noise ratio is 2 or 5,the inversion results can also be well estimated using our sphericalwave inversion approach.

5.Field data example

We employ a field data example to verify our complex sphericalwave inversion approach.The field data are acquired from the seismic traces near borehole(known CDP)which are located in an oilfield of eastern China and have been amplitude-preserved processed.The seismic data are partially stacked,and the incident angle ranges are 8?-16?and 17?-25?.So,the dominant incident angles of the observed seismic data areθp1=12?andθp2=21?,respectively.Next,the different frequency components of partially stacked seismic data at the frequency range of 20-40 Hz and 30-50 Hz are preserved by the continuous wavelet transform,and their dominant frequencies respectively are 30 Hz and 40 Hz.There is no specific rule on how to determine the frequencies used in the inversion.Since the energy of seismic reflection data is concentrated near the dominant frequency(35 Hz),the frequency components near the dominant frequency are utilized for inversion.With the real-valued seismic data of different incident angles and frequencies,the corresponding real-valued wavelets can be extracted.The imaginary parts of seismic data and wavelets can be further obtained by Hilbert transform.To illustrate the advantages of the proposed approach,the conventional plane-wave inversion is also implemented using the seismic datasets with the incident angles of 12?,21?,and 30?,respectively.Based on the inverted complex spherical-wave EI,we implement the complex sphericalwave elastic inversion for P-and S-wave velocities and density using the field datasets,the workflow is shown in Fig.14.

Fig.15 displays the inversion results of complex spherical-wave EI of different incident angles and frequencies.Different from the synthetic data examples,the dominant frequencies of observed field seismic data respectively are f1=30 Hz and f2=40 Hz,and the dominant incident angles of observed seismic data respectively are θp1=12?andθp2=21?.The plane-wave EI inversion results are displayed in Fig.16 for comparison.As can be seen from Figs.15 and 16,compared with the imaginary part of the spherical-wave EI,the plane-wave EI inversion results are much the same as the real part of the spherical-wave EI inversion results.

With the inverted plane-wave EI and complex spherical-wave EI,the corresponding P-and S-wave velocities and density can be further extracted,as shown in Fig.17a and b.It can be seen from Fig.17 that the inversion results estimated by our spherical-wave approach have higher resolution and continuity than that estimated by the plane-wave approach.Fig.18 further displays the inversion results and the well logs near the borehole.In Fig.18,the P-and S-wave velocities and density estimated from our approach can match the filtered well-logging data,and are better than that estimated by plane-wave approach.The field data example verifies the feasibility and practicability of our complex spherical-wave inversion approach using spherical-wave amplitude,phase,and frequency information.

Fig.17.Comparisons of P-wave velocity,S-wave velocity,and density between the(a)plane-wave inversion results and(b)spherical-wave inversion results.

6.Conclusions

Based on the theory of wave decomposition,the simple harmonic wave is convoluted with SRC to obtain the spherical-wave synthetic seismic seismogram at a certain angular frequency.It is assumed that geophone can only record the real part of complex seismic trace,which describes the physical properties of seismic motion.We further demonstrate that the imaginary parts of seismogram and wavelet can be got by the Hilbert transform of their corresponding real parts.We then propose the concept of complex spherical-wave elastic impedance(EI)and derive a complex spherical-wave EI equation.To fully exploit the amplitude and phase information of spherical-wave reflection coefficient(SRC),the complex spherical-wave EI inversion and elastic parameters extraction are implemented.Our inversion approach consists of two steps:estimating complex spherical-wave EI from seismic data of different frequency components and incident angles with Bayesian framework,and extracting P-and S-wave velocities and density from the inverted complex spherical-wave EI.The synthetic data and field data examples show that our spherical-wave inversion approach can reasonably estimate the velocities and density,which demonstrates the feasibility and practicability of using the amplitude and phase information of spherical wave reflection with different offsets and frequencies to estimate the elastic parameters.It is worth noting that the P-and S-wave velocities and density estimated from our spherical-wave inversion approach are better than that estimated by plane-wave inversion approach.

It must be emphasized that we have made some assumptions in the real world,only the real part of seismic signal can be recorded,the P-wave velocity in upper medium is assumed to be known a priori and it can be obtained from tomographic velocity approximately.Actually,seismic wavelets with different dominant frequencies have a certain bandwidth.The spherical-wave reflection coefficient with the angular frequencyωnis used to approximately describe the seismic reflection of the wavelet with the dominant frequency ofωn.

Fig.18.Comparisons of P-wave velocity,S-wave velocity,and density between true values(red solid curves),spherical-wave inversion results(black solid curves),and plane-wave inversion results(blue solid curves)near the borehole.

Acknowledgements

We would acknowledge the sponsorship of the Marine S＆T Fund of Shandong Province for Pilot National Laboratory for Marine Science and Technology(Qingdao)(Grant No.2021QNLM020001-6)and National Natural Science Foundation of China(42030103,41974119)and Science Foundation from Innovation and Technology Support Program for Young Scientists in Colleges of Shandong province and Ministry of Science and Technology of China(2019RA2136).