Full Waveform Inversion of Cross-hole Radar Data Using Envelope Objective Function Xintong Liu, Sixin Liu*, Lei Fu College of Geo-exploration Science & Technology, Jilin University No.938, Ximinzhu Street, 130026, Changchun, China
[email protected] Abstract—The full waveform inversion (FWI) has been used as a high resolution imaging method for cross-hole radar inversion implementation, however, when FWI is used to process field data, it encounters a variety of problems and causes the inversion to fall into a local minimum. One problem is that the GPR data lacks low frequency information. An effective way is to provide an accurate initial model for FWI. The fact is that the underground media is unknown which results in difficult to obtain ideal result in FWI. By means of the derivation of the objective function, the gradient formula of the envelope waveform inversion (EWI) is derived by taking the derivation of the EWI misfit function with respect to the model parameter. By comparing the inversion results of the EWI with that of the traditional FWI without low-frequency information in the observed GPR data, we found that the EWI can effectively restore missing low-frequency information and has better inversion ability for low-frequency missing data. Keywords—full waveform inversion; cross-hole radar; envelope; low frequency component; initial model
I.
INTRODUCTION
Ground penetrating radar tomography technology has been widely used in environmental, hydrological and engineering geophysical prospecting. The cross-hole radar collects the direct wave by placing the transmitting antenna and receiving antenna respectively in two adjacent boreholes, and the transmission wave contains the spatial distribution information of dielectric constant and conductivity. Traditional cross-hole radar imaging methods are based on ray theory, such as traveltime inversion and maximum amplitude of the first cycle signal inversion. Due to the high frequency assumption in the ray-based method, the resolution of these imaging methods is around the diameter of first Fresnel zone. In order to improve the resolution, methods based on full waveform inversion (FWI) techniques have been introduced into the radar data processing recently. These methods include many kinds of Born approximation iteration methods for seismic data [1-4]. Kuroda et al. [5] and Ernst et al. [6] used the FWI method for cross-hole radar data imaging, and achieved good results. The difference is that Kuroda et al. [5] only inverted the permittivity, but Ernst et al. [6] used the cascade way to update the permittivity and conductivity simultaneously. Meles et al. [7] proposed vector FWI technology, and simultaneously inverted the permittivity and electrical conductivity. Wu et al. [8] presented the detailed
derivation of the theoretical foundation of the FWI, and used the steepest descent method to update the model parameters. Meng et al. [9] extended the cross-hole radar waveform inversion to the frequency domain, and used the logarithmic objective function and conjugate gradient method to update the model. Babcock and Bradford [10] implemented a targeted reflection waveform inversion algorithm to quantify thin-layer permittivity, thickness, and conductivity for NAPL thin (≤1⁄2 dominant wavelength) and ultrathin (≤1⁄8) layers by using GPR reflection data. The FWI is widely used in the inversion of the local optimization method. This method is highly dependent on the initial model, however, it can not start from the accurate initial model in the underground media when there is no priori information. In addition, with the absence of low frequency information and the presence of noise in the actual collected data, the local optimization method is easy to fall into the local extremum [11]. In the field of full waveform inversion research, it is common to choose a more accurate initial model to reduce the degree of non-linearity for better results. One of the more common methods in the existing methods is to use the results of ray as an initial model for full waveform inversion. If the waveform inversion method is used to build the initial model, the use of low frequency information is particularly important. However, for the real field data, lowfrequency information is often difficult to obtain because of the antenna bandwidth or data acquisition filter design, etc. Taking borehole radar with center frequency of 100 MHz as an example, it is difficult to obtain effective information below 30 MHz. Therefore, reconstruction of effective low frequency information is very important for waveform inversion of cross hole ground penetrating radar. The envelope of radar data can reconstruct the large-scale information of the subsurface model and it is appropriate to use it to invert the long wavelength components. As early as 2011, Bozdag [12] carried out the study of the envelope objective function, but he did not use the envelope objective function for FWI. Subsequently, Chi et al. [13] and Wu et al. [14] used the FWI based on acoustic wave equations to process seismic data respectively. The envelope objective function was used to effectively invert the velocity information of underground media, and use the result obtained by this method as the initial model of the conventional FWI to further obtain a more accurate results with high resolution. However, the above studies for the envelope-based FWI are
This work was supported in part by the Natural Science Foundation of China under Grant (41504085) and China Postdoctoral Science Foundation (801151010423)
978-1-5386-5777-5/18/$31.00 ©2018 IEEE
currently focused on the inversion of seismic data. Research on this aspect of radar is still in its infancy. In this paper, based on the time domain transmission radar waveform inversion, the envelope objective function and the derivation method are used to derive the basic theory of envelope waveform inversion (EWI) in detail. The EWI is successfully implemented to the field of GPR. Conjugate gradient method is used to invert two kinds of physical parameters (dielectric constant and conductivity) simultaneously. By simulating the absence of low-frequency information in the original data, the inversion performance of this method for low-frequency missing cases is verified. The inversion results can not only evaluate the underground media, but also can be used as an initial model to reduce the degree of nonlinearity and further participated in traditional FWI to obtain higher resolution results. The forward modelling of the FWI method uses a high-order finite-difference time-domain (FDTD) method with CPML absorption boundary. Due to the huge amount of data in the inversion process, Matlab parallel computing mode of MDCE is adopted to improve the program efficiency with ordinary microcomputer. II.
THEORY
A. Traditional FWI The main purpose of the FWI is to invert the space distribution of the permittivity and conductivity by minimizing the objective function:
S ( , )
1 s s ]d , [Es ( , ) Eobs ]d , [Es ( , ) Eobs 2 s d (1)
where the E ( , ) and E are the forward and observed data, respectively. Equation (1) is the L2 norm of the difference between the forward and observed data in three ranges including the source positions s, receiving position d and the observation time . s
s obs
In the time domain, the Maxwell equations can be expressed as (Meles, et al, 2010),
(x) t (x)
E s ( x, t ) J s ( x, t ) (2) 0 t H s (x, t ) 0
(2) is abbreviated as follows,
E s J s M ( , ) s H 0
(3)
Ignore the magnetic field, yields to,
ME J or E = M J (4) In the time domain, the electric field of arbitrary position can be considered as the convolution between the source wavelet and the corresponding Green's function g, Es g J s (5) where g is the Green’s operator of M. It is clear that the M-1 must contain the discrete approximations to the Green's functions (Pratt, 1999; Pratt and Shipp, 1998). Thus, s
s
s
-1 s
M 1 g
(6) The gradient corresponding to the objective function is obtained by the derivation of equation (1) with respect to the model parameters. For simplicity, we use p to represent model parameters and ,
S Es s Es Eobs p p s d
(7)
S v t g t r d dt p s d
(8)
The gradient of the objective function is obtained from the zero-lag cross-correlation between the forwarded wavefields and the backprogated wavefields.
where,
s r Es Eobs
Since the backward propagated field
(9)
g t r d is
inverted in time. B. Envelope Waveform Inversion In this paper, we introduce the envelope objective function into the vector waveform inversion technique of Meles et al. (2010) and reconstruct the basic theory of envelope waveform inversion. Different from Meles' method of deriving the gradient by using the first-order approximation of perturbed objective function, this paper uses the method of derivation of the envelope objective function to derive the gradient. (1) is the most widely used traditional objective function in inversion, this paper presents an envelope objective function based on the objective function, that is, the L2 norm of the difference between the synthetic and observed envelope : S ( , )
1 obs [E env ( , ) Eobs env ]s , d , [ E env ( , ) E env ]s , d , 2 s d
(10) Here, E env ( , ) is the envelope of E( , ) , defined as
( , ) is a Eenv ( , ) E 2 ( , ) E 2 ( , ) . E Hilbert transform of E ( , ) . Obtain the gradient by deriving the model parameters through the objective function (8). For brief, we use p to represent model parameters and .
E S env Eenv ( , ) Eobs env p s d p
(11)
Eenv ( , ) Eobs env E( , ) S Eenv ( , ) E( , ) obs p p Eenv ( , ) Eenv s d H E( , ) Eenv ( , ) (12) An expression for the partial derivatives in eq. (12) in terms of the forward-modelling matrix eq. (4) can now be obtained by
taking the partial derivative of both sides of eq. (4) with respect to the parameter p:
E M s Ε M 1 v , M E or p p p s
(13)
,k
where we have introduced the virtual source term
M s v E p
(14)
where,
v is the virtual source vector, the expression for different model parameters and are different. Substituting (13) into (12) we can obtain
S E s radjoint = v*[M 1 ] renv (15) p p s d s d where,
radjoint
Insert equation (6) into (15), yields to the integral form,
d
s
d
s
d
s
d
γ r γ γ γ r γ γ
env
env
(22)
(23)
renv [Εenv ( , ) Eobs env ]d ,
(24)
γ [Eenv ( , C ,k ) Eenv ( , )]Td ,
(25)
γ [Eenv ( , C , k ) Eenv ( , )]Td ,
(26)
and
are stability factors, which should be
III.
(16)
S v t g t r d dt (17) p s d
Since the backward propagated field
Here,
s
carefully selected in inversion.
obs Eenv ( , ) Eenv E( , ) Eenv ( , )
E ( , ) Eobs env H env E( , ) Eenv ( , )
,k
s
g t r d is
SYNTHETIC DATA
Unlike the traditional FWI method that compares complete waveforms, our method compares the envelope between the synthesized and observed data. Figure 1a shows a Ricker wavelet (red line) and its envelope (blue line). It can be seen from the comparison that the Ricker wavelet signal itself changes drastically while the envelope waveform changes relatively slowly. Figure 1b shows the comparison between the wavelet spectrum (red line) and the envelope spectrum (blue line). It can be clearly seen that the low frequency information is created in the envelope spectrum by non-linear conversion.
inverted in time, the equation (17) can be considered as the zero-lag crosscorrelation between the backward propagated field and the source wavefield. It is consistent with the gradient results obtained using the perturbation of the objective function by Meles et al.(2010). According to (17), the gradient can be expressed as,
ˆ T Rs t Es G S S s T s ˆ R s E G
(18)
ˆ R represent backward propagated process of the where, G s sum of residual R , and T
s
Rs r
(19)
d
In order to increase the convergence, we use the conjugate gradient method (Scales, L., 1985) to update the model parameters:
x k 1 x k ,k S x k x k 1 x k ,k S x k
(20) (21)
where, S x k and S x k are gradients for permittivity
and conductivity at the k-th iteration. ,k and ,k are the iterative step length of the permittivity and conductivity.
Fig 1. (a) Ricker wavelet and its envelope; (b) spectra of the wavelet and its envelope.
The traditional cross-borehole radar FWI can obtain good inversion results when the low-frequency information of radar data is sufficient and the initial model is accurate enough, but FWI often can not recover the large-scale background information when the low-frequency information is missing. The update of the model will be in the wrong direction, and eventually the inversion will not converge. The lack of lowfrequency information in radar data is a condition that exists in many radar examples. The traditional full waveform inversion method (FWI) can not work well under similar conditions. However, the envelope-based waveform inversion (EWI) has the ability to update the long wavelength information of the model in the absence of large amounts of low-frequency data, the reason is that its adjoint wavefield contains low-frequency information. As shown in Fig. 2, we remove the low frequency components below 30MHz for the Ricker wavelet. In Fig. 2a, the Ricker wavelet and its envelope after removing the low frequency components show that the Ricker wavelet changes
more fiercely after removing the low frequency components. As shown in Fig. 2, we remove the low frequency components below 30MHz for the Ricker wavelet. In Fig. 2a, the Ricker wavelet and its envelope after removing the low frequency components show that the Ricker wavelet changes more fiercely after removing the low frequency components, while the envelope curve is still stable. Figure 2b shows the comparison of the spectra. Compared with Figure 1b, it can be clearly seen that after removing the low frequency components, the low frequency part of the Ricker wavelet has been removed. However, the low frequency portion of the envelope is almost unaffected, and significant low frequency information is still added.
wavelength information, both inversion methods can obtain good inversion results. Among them, the resolution obtained by the traditional FWI method is slightly higher than that based on the envelope waveform inversion. The reason is that the information used in the FWI are more complete but there is also a higher degree of non-linearity. From the inversion results in Fig. 4, when the waveform information is complete, the inversion result is better than the envelope waveform inversion. Especially for the two pipe-shaped anomaly bodies in the middle layer, the conductivity of the conventional waveform inversion is more accurate. The result of EWI has a very high resolution, which is consistent with the conclusion obtained by the seismic EWI.
Fig 2. (a) Ricker wavelet without low frequency data and its envelope; (b) spectra of Ricker wavelet without low frequency data and its envelope.
In order to verify the capability of FWI tomography in the time domain, the relative permittivity and conductivity model are built in Figure 3a and 3b, respectively. The model size is 6 m in both x- and y- directions, and contains two pipeline anomalies buried in the three-layer medium with diameter of 0.5 m. The relative permittivity and conductivity of the two anomalous with depth of 3 m are 6.5 and 0.008 S/m, respectively. The relative permittivity and conductivity of the top and the bottom-layer medium are 5 and 0.001 S/m, respectively. The relative permittivity and conductivity of the mid-layer medium are 5.5 and 0.0028 S/m, respectively. The initial permittivity and electrical conductivity model are the same as the mid-layer medium with the relative permittivity and conductivity of 5.5 and 0.0028 S/m. Other parameters are the same as the first example.
Fig 3. Real model (a) Relative permittivity; (b) Conductivity.
First, we compare the results of the two methods when using the complete waveform information. The uniform initial model has a relative dielectric constant of 5.5 and a conductivity of 0.0028 S/m, which are consistent with the background physical parameters of the original model. Since the original data contains complete and accurate long-
Fig 4. (a) Relative permittivity result of traditional FWI; (b) Conductivity result of traditional FWI; (c) Relative permittivity result of EWI; (d) Conductivity result of EWI.
Secondly, we simulate the low-frequency missing condition. In full waveform inversion, the low-frequency component of the data is very important. If the low-frequency component is sufficient, it can be used in FWI to restore the large-scale background information of the formation through a series of frequency strategies. It is different from highfrequency information missing only affects the resolution. Once the low-frequency information is missing, it will not only affect the resolution of the inversion result, but also affect the stability of full-waveform inversion. In Figure 5, it can be seen that the loss of low-frequency information has a very strong effect because we simulate a situation where a large amount of low-frequency information is missing and causes a center frequency shift. During the numerical simulation, it was found that the shift of the center frequency has a great influence on the inversion of the conventional waveform inversion. Also, we only use the uniform model as the initial model, which greatly increases the nonlinearity of the inversion. By comparison, the envelope waveform inversion can still get good results and provide the initial model for traditional waveform inversion. In Figure 5, the number of
iterations of the traditional FWI result is 20 times. Due to the lack of low-frequency information, when the number of iterations is increased, the result of the inversion can not converge. Instead, the inversion falls into a local minima and updates the model in the wrong direction, with the objective function failing to converge to the minimum. When the lowfrequency information is missing, the results of traditional waveform inversion can not be obtained consistent with the FWI using the complete waveform information, no matter the dielectric constant or the conductivity. When the low frequencies are missing, the results of the envelope inversion are better than the complete information. This is because when the envelope information is complete, the contribution of lowfrequency information to inversion is much greater than that of high-frequency information, making the inversion focused on large-scale stratigraphic information. When filtering out a part of the low-frequency information, the weight of the highfrequency information is increased in the inversion, so that the situation of the middle two abnormal bodies is restored, but the resolution is still lower than the complete information FWI. It can be clearly seen that the inversion resolution of the upper and lower strata has decreased, especially at the corners. We can use this to build a frequency strategy for envelope waveforms inversion.
waveforms, and in particular verifies its ability to invert the initial wave. Through the numerical simulation of cross-bore radars, it is proved that this method has high application value, especially for the low-frequency missing data. Even if traditional waveform inversion does not work, envelope waveform inversion can provide good inversion results. Regardless of dielectric constant or conductivity, envelope waveform inversion is significantly superior to FWI with missing data of low frequency information. Envelope waveform inversion is a complementary method to the waveform inversion of the cross-bore radar, which makes the cross-borehole radar waveform inversion method can be applied to more examples.
References [1] [2]
[3] [4] [5] [6]
[7]
[8] [9]
[10]
Fig 5. Inversion result without low-frequency information. (a) Relative permittivity result of traditional FWI; (b) Conductivity result of traditional FWI; (c) Relative permittivity result of EWI; (d) Conductivity result of EWI.
IV.
[11] [12]
CONCLUSION
The research of envelope waveform inversion mainly faces the problem of low frequency information missing in GPR data. Envelope inversion is a waveform inversion method which is more stable and less nonlinear than the traditional full waveform inversion. This paper successfully introduces the envelope objective function into the inversion of gpr
[13]
[14]
A. Tarantola, “Inversion of seismic-reflection data in the acoustic approximation,” Geophysics, vol. 49, no. 8, pp. 1259–1266, Aug 1984. R. G. Pratt, C. Shin, and G. J. Hicks, “Gauss-Newton and full Newton methods in frequency domain seismic waveform inversion,” Geophysical Journal International, vol. 133, no. 2, pp. 341–362, Apr 1998. R. G. Pratt, “Seismic waveform inversion in the frequency domain, Part 1: Theory and verification in a physical scale model,” Geophysics, vol. 64, no. 3, pp. 888–901, May 1999. R. G. Pratt, and R. M. Shipp, “Seismic waveform inversion in the frequency domain, Part 2: Fault delineation in sediments using crosshole data,” Geophysics, vol. 64, no. 3, pp. 902–914, May 1999. S. Kuroda, M. Takeuchi, and H. J. Kim, “Full waveform inversion algorithm for interpreting cross-borehole GPR data,” SEG Technical Program Expanded Abstracts 2005, pp. 1176–1180. J. R. Ernst, H. Maurer, A. G. Green, and K. Holliger, “Full-waveform inversion of crosshole radar data based on 2-D finite-difference timedomain solutions of maxwell’s equations,” IEEE Transactions on Geoscience and Remote Sensing, vol. 45, no. 9, pp. 2807–2828, Sept 2007. G. A. Meles, J. Van der Kruk, S. A. Greenhalgh, J. R. Ernst, H. Maurer, and A. G. Green, “A new vector waveform inversion algorithm for simultaneous updating of conductivity and permittivity parameters from combination crosshole/borehole-to-surface GPR data,” IEEE Trans. Geosci. Remote Sens, vol. 48, no. 9, pp. 3391–3407, Sep 2010. J. Wu, S. Liu, Y. Li, et al., “Study of cross-hole radar tomography using full-waveform inversion,” Chinese Journal of Geophysics, vol. 57, no. 5, pp. 1623–1635, May 2014. X. Meng, S. Liu, L. Fu, et al., “Frequency domain waveform inversion of cross-hole GPR data based on a logarithmic objective function,” Chinese Journal of Geophysics, vol. 59, no. 5, pp. 1875-1887, May 2016. E. Babcock, and J. H. Bradford, “Reflection waveform inversion of ground-penetrating radar data for characterizing thin and ultrathin layers of nonaqueous phase liquid contaminants in stratified media,” Geophysics, vol. 80, no. 2, pp. H1-H11, Mar 2015. J. Virieux, and S. Operto, “An overview of full-waveform inversion in exploration geophysics,” Geophysics, vol. 74, no. 6, pp. WCC1– WCC26, Nov/Dec 2009. E. Bozdağ, J. Trampert, and J. Tromp, “Misfit functions for full waveform inversion based on instantaneous phase and envelope measurements,” Geophys. J. Int., vol. 185, no. 2, pp. 845–870, May 2011. B. Chi, L. Dong, and Y. Liu, “Full waveform inversion method using envelope objective function without low frequency data,” Journal of Applied Geophysics, vol. 109, pp. 36–46, Oct 2014. R. S. Wu, J. Luo, and B. Wu, “Seismic envelope inversion and modulation signal model,” Geophysics, vol. 79, no. 3, pp. WA13– 21WA, May-Jun 2014.
