The Fourier transform of controlled-source time-domain electromagnetic data by smooth spectrum inversion

SUMMARY In controlled-source electromagnetic measurements in the near zone or at low frequencies, the real (in-phase) frequency-domain component is dominated by the primary ﬁeld. However, it is the imaginary (quadrature) component that contains the signal related to a target deeper than the source–receiver separation. In practice, it is difﬁcult to measure the imaginary component because of the dominance of the primary ﬁeld. In contrast, data acquired in the time domain are more sensitive to the deeper target owing to the absence of the primary ﬁeld. To estimate the frequency-domain responses reliably from the time-domain data, we have developed a Fourier transform algorithm using a least-squares inversion with a smoothness constraint (smooth spectrum inversion). In implementing the smoothness constraint as a priori information, we estimate the frequency response by maximizing the a posteriori distribution based on Bayes’ rule. The adjustment of the weighting between the data misﬁt and the smoothness constraint is accomplished by minimizing Akaike’s Bayesian Information Criterion (ABIC). Tests of the algorithm on synthetic and ﬁeld data for the long-offset transient electromagnetic method provide reasonable results. The algorithm can handle time-domain data with a wide range of delay times, and is effective for analysing noisy data.


I N T R O D U C T I O N
Controlled-source electromagnetic (CSEM) methods are generally divided into frequency-domain electromagnetic (FEM or FDEM) and time-domain (or transient) electromagnetic (TEM or TDEM) methods, depending on the waveform of the transmitted electrical current. In the FEM method, we measure the electromagnetic (EM) response caused by a sinusoidally alternating source current. The measured response consists of both primary and secondary ®elds. The primary ®eld is the freespace EM ®eld produced directly by the source current, so it contains no information about the resistivity of the Earth. The secondary ®eld, by contrast, is caused by electrical currents induced in the Earth by the primary ®eld and depends on the resistivity of the Earth. In the far zone, where the source± receiver separation, r, is much greater than the skin depth d (m), equal to 503(r/f ) 1/2 , in which r is the resistivity (V m) and f is the frequency (Hz), the secondary ®eld is comparable to the primary ®eld. In contrast, the secondary ®eld becomes extremely small compared to the primary ®eld in the near zone, r%d (West & Macnae 1991). This property imposes a major limitation on the FEM measurements. In the TEM method, we record a transient EM response after the source current is abruptly turned off. Therefore, the secondary ®eld due to induced eddy currents can be measured in the absence of the primary ®eld (Nabighian & Macnae 1991). The TEM and FEM ®elds are related by the Fourier transform, and the FEM secondary ®eld can thus be obtained through the Fourier transform of TEM data.
There have been several studies of Fourier transform techniques based on inversion. Sacchi & Ulrych (1996) proposed an inversion method using Bayes' rule. They used a Cauchy a priori distribution to stabilize the inversion and applied it to the 2-D spectral analysis of seismic records. Duijndam et al. (1999) used a diagonal weighting scheme corresponding to sampling intervals and employed a diagonal regularization to stabilize the inversion of irregularly sampled data. Oldenburg (1976) applied Backus±Gilbert linear inverse theory to the Fourier transform of data contaminated by noise and gaps. These approaches are more useful than the conventional Fast Fourier Transform (FFT) when time-series data are contaminated by noise or sampled irregularly. We desire the FEM response sampled logarithmically and need to deal with noisy data. Thus, the development of a new numerical Fourier transform technique based on inversion is necessary. A difference from these previous studies is that we impose a smoothness constraint on the FEM response. This is a reasonable assumption because our FEM response is usually continuous with respect to frequency. Egbert (1992) estimated the discrete time impulse response of magnetotelluric (MT) impedance from the frequency-domain MT impedance data using a least-squares inversion with a smoothness constraint. This is the reciprocal to estimating the FEM response. However, this study considered the natural EM ®elds of magnetotellurics, not the CSEM ®eld analysed here.
In this paper, we ®rst compare the FEM response with the TEM response in the near zone and show the importance of the imaginary (quadrature) component of the FEM response and the superiority of the TEM measurements in the near zone. Second, we will describe our least-squares inversion theory, which uses a smoothness constraint for the estimation of the FEM response. Finally, we apply our technique to synthetic data with Gaussian noise and real TEM data.
As a result of transmitting an alternating current cos(v 0 t), with an angular frequency v 0 , we measure a signal f (t), where A and w are the amplitude and the phase difference, respectively. The Fourier transform of f (t) is given by where d D is the Dirac delta function. In the near zone, r%d, w becomes extremely small, meaning that the imaginary part of F(v) becomes a small fraction of the real part. In the received signal, the second term on the right-hand side of eq. (1) corresponds to the imaginary part and it becomes much smaller than the ®rst term, corresponding to the real part. This means that almost the whole waveform of the received signal consists of the in-phase component. Consider a horizontally two-layered earth model ( Fig. 1) consisting of a conductive 1 V m basement underneath a 1000 m thick resistive overburden of 100 V m. A horizontal electric dipole source is used to excite the earth and a receiver is located on the surface 1000 m from the source. Fig. 2 shows the frequency-domain vertical magnetic ®eld calculated with the 1-D modelling code EM1D (Kim et al. 1997). For comparison, the FEM response for a 100 V m homogeneous half-space is also shown. One can see that the primary ®eld, 1/4pr 2 =7.96r10 x8 A m x1 , calculated from the Biot±Savart law, is dominant in the real component at low frequencies, and that it masks signals from the lower layer. On the other hand, the response of the conductive basement appears in the imaginary component at low frequencies. In ®eld practice, accurate measurements of the imaginary component are dif®cult in the near zone because the signal corresponding to the imaginary component is very small, as shown in Fig. 2 and eq. (1).

Impulse and step responses
When the electric current of the source is suddenly turned off, the step h off (t) and impulse dh(t)/dt responses are represented by the inverse Fourier transforms of the FEM response H(v) as  follows (Newman et al. 1986; Appendix A): where Re[ ] and Im[ ] denote the real and imaginary parts of the complex argument, respectively. In addition, when the current is suddenly turned on, the step response h on (t) is given by These equations are derived from the causality of the transient responses and the measurement system response, as discussed in Appendix A. We calculated these responses for the model of Fig. 1 using the Gaver±Stehfest inverse Laplace ®lter (Knight & Raiche 1982;Hasegawa 1985) as shown in Figs 3 and 4. As Spies & Frischknecht (1991) pointed out, the greatest advantage of TEM soundings is their ability to locate interfaces at depths that are large compared to the source±receiver separation. The impulse response (Fig. 3) and the step response for the turn-off of source current (Fig. 4) are sensitive to the presence of the conductive basement. By contrast, the step response for the current turn-on shows no evidence of the conductive basement. This is because the primary ®eld of the ®rst term on the right-hand side of eq. (8) obscures the secondary ®eld of the second term caused by the imaginary part.
The FEM response is generally convenient for the deconvolution of a measurement system response, and is easier to model for 1-D or multidimensional interpretations in comparison with the TEM response. However, the high-quality acquisition of the imaginary component in the near zone is dif®cult because the observed signal in the FEM measurements is dominated by the primary ®eld corresponding to the real component. Eqs (4) and (6), however, suggest that we can evaluate the imaginary part of the FEM data from the TEM signal through the Fourier transform.
In practical measurements, acquired TEM signals tend to be distorted by low-pass ®ltering within instruments to avoid the aliasing, which makes the FEM response band-limited. Usually it appears as a signi®cant depression of the early-time TEM response (Strack 1992). However, it is a bene®t for the accurate acquisition of the late-time TEM signal because the dynamic range required for the measurement of the whole waveform of the TEM signal is reduced and we can easily increase the gain for accurate measurement of the late-time TEM signal.

Discretization of frequency response
The Fast Fourier Transform (FFT) is widely used for the spectral analysis of time-series data. However, our TEM data decay slowly and the FEM responses are usually sampled as a function of the logarithm of frequency. In addition, the TEM data are contaminated by noise. Therefore, a technique is required to handle such noisy data in terms of logarithmic time and frequency sampling. Discretizing the integrals of eqs (3)±(6), we can, in principle, solve for the frequency response by a linear inversion. As Fig. 2 shows, the FEM response is not a line spectrum but a smoothly varying continuous function. Therefore, we distribute nodal points with a constant interval  in the logarithm of frequency, and approximate the FEM response with local linear interpolation functions as shown in Fig. 5.
Integrating the product of the interpolation function and the cosine or sine function over each division of v, we can approximate eqs (3)±(6) with the matrix equation and A is an NrM linear operator matrix calculated through the integration over each division. The derivation of the elements of A for eqs (3)^(6) is described in Appendix B. Since the FEM response is a continuous function, we desire the response to be represented by the nodal points as closely as possible. Accordingly, we consider the overparametrized case of MiN.
Maximum a posteriori estimation Throughout this paper, we assume that our data contain Gaussian noise and that we only know relative values of data covariance in the form of the relative covariance matrix C 0d of the data. This assumption is demanded for the use of the least-squares method and the analytical calculation of Akaike's Bayesian Information Criterion (ABIC). The reliability of the assumption of Gaussian noise is beyond the scope this paper, but if the data contain outliers, we can easily identify them in the time domain and remove them, as suggested for magnetotelluric measurements (Egbert 1992). In addition, the data acquisition of controlled-source electromagnetic methods can be easily modi¢ed so that the noise can be described with a Gauss distribution on the basis of the central limit theorem. That is, the distribution of the mean value of randomly sampled data approaches the Gauss distribution as sample size increases, even if the data are not distributed with a Gauss distribution (Menke 1989). Under these assumptions, the probability distribution is given by where s 0 2 is a correction factor to compensate C 0d , det denotes the matrix determinant and T signi¢es the transpose. It is necessary to employ appropriate a priori information for the model parameters to solve our underdetermined problem. Assuming a smoothly varying FEM response, we de®ne the prior distribution of m as where the MrM matrix D is the second-order di¡erence operator whose rank is Mx2. Although the parameter m adjusts the reliability of the a priori information and has to be speci®ed to stabilize the inversion, we assume that m has already been given subjectively here. The determination of m will be considered later.
According to Bayes' rule (Tarantola 1987), we write the posterior distribution where the numerator of the right-hand side is the joint probability density and the denominator p(d), which is the marginal probability density and independent of m, is assumed to be a constant (Sen & Stoffa 1995). Accordingly, where The model existing at the peak of the posterior probability distribution is the most probable model. This estimation, based on Bayes' rule, is called the maximum a posteriori (MAP) estimation (Sacchi & Ulrych 1996). Searching for the maximum of the distribution of eq. (14) corresponds to an estimation of the model that minimizes the objective function U(m). Since the covariance C 0d is a symmetric positive de¢nite matrix, it can be decomposed to C 0d =LL T by the Cholesky factorization (Sakamoto et al. 1986), where L is a lower-triangular matrix. Thus, eq. (15) can be rewritten as follows: and the weighting matrix W is W=L x1 . Minimization of eq. (16) is achieved via the least-squares method. Applying hU(m)/hm=0, we obtain the normal equation In this paper, we factorize the matrix G= QR using the modi¢ed Gram^Schmidt method (Saito 1983), where R is an MrM non-singular upper-triangular matrix and Q is an (N+M)rM orthogonal matrix. That is, Q T Q=I, where I is the unit MrM matrix. By substituting G= QR into eq. (19), the normal equation can be rearranged into Consequently, the estimator m* for our least-squares problem is obtained by solving eq. (20). The covariance matrix of m*, corrected by s 0 2 , is derived in a similar way to Menke (1989) and is given by In the following numerical experiments, we show the standard deviation of m* as error bars that are the square roots of the diagonal elements of C m* .

ABIC minimization method
The prior distribution, eq. (11), is strongly affected by m, that is, the estimation of m depends on m. Therefore, the determination of the value of m is particularly important. Regarding m as the parameter that controls the estimation of the model parameters m, we call m a hyperparameter, which is often called a trade-off parameter or a regularization parameter. In addition, since s 0 2 is also an unknown parameter, we regard it as one of the hyperparameters. Taking these parameters into account, we rewrite p(d|m)[p(d|s 0 2 , m) and p(m)[p(m|s 0 2 , m) and integrate the joint probability density with respect to m in order to evaluate the marginal probability density p(d|s 0 2 , m) for d as follows: The resulting distribution indicates that the probability of occurrence of d is a function of m and s 0 2 and is no longer a function of m. Now our problem becomes that of the estimation of m and s 0 2 from the observed data d. We need to choose parameters m and s 0 2 such that they maximize p(d|s 0 2 , m). Akaike (1980) called this probability the Bayesian likelihood and developed the Akaike Bayesian Information Criterion based on the entropy maximization principle, which is represented by q in this paper and is de®ned as where N H is the number of hyperparameters. We determine the optimum values of the hyperparameters and N H by minimizing the ABIC and call this scheme the ABIC minimization method. This method has been applied to many geophysical problems: by Noro (1989) for curve ®tting, by Murata (1993) for estimating sur®cial density from gravity data, by Mitsuhata (1994) and Uchida (1993) for 1-D and 2-D inversions of magnetotelluric data, respectively, and by Yokota et al. (1997) for seismic tomography. Substituting eqs (10) and (11) into eq. (22) and calculating the integral (Appendix B), we obtain q N À 2 ln2np 2 0 À 2 lndet W À M À 2 ln k lndet G T G Uk min ap 2 0 4 , where U min is the minimum of U, which is evaluated by substituting m* into eq. (16). Considering hq/hs 0 2 =0 to estimate the optimum s 0 2 , we obtain Consequently, by using the above relationship, it can be shown that In eq. (26), since det(G T G)=det(R T R)={det(R)} 2 , the third term on the right-hand side is computed in the following manner: To search for the optimum m that minimizes q, we carry out a 1-D search using the Golden section method (Press et al. 1992). Fig. 6 shows a¯ow chart of the inversion algorithm. Hereafter, we call our smoothness-constrained least-squares inversion a smooth spectrum inversion.

A P P L I C A T I O N S A N D R E S U L T S Synthetic data example
To check the validity of the smooth spectrum inversion, we tested it using synthetic data sets consisting of impulse and step responses. These responses were calculated for the two-layered model of Fig. 1. We added 1 per cent independent Gaussian noise to each datum, which may not be realistic, but the objective of this experiment was to check the validity of the code when the noise level is low. These synthetic data sets were composed of 52 data points sampled logarithmically from t=10 x4 to 10 s and at t=0. The number of model parameters representing the FEM response was 61 for the cosine transform and 60 for the sine transform. The nodal points of each parameter were logarithmically arranged between the limits of f=0.01 and 7943 Hz. In the cosine transform, the static component was also estimated. The smooth spectrum inversion allows us to extrapolate the frequency response beyond the Nyquist frequency and interpolates between the static ®eld and the fundamental frequency estimated from the record. Figs 7(a) and (b) show the inversion result and the conventional FFT result from the impulse response time-domain data, respectively, compared with the values obtained from EM1D. The inversion result agrees well with the correct values in the frequency range lower than 1 kHz. The data predicted by the inversion ®t the synthetic data perfectly, as shown in Fig. 7(c). In the application of the FFT, we were forced to use 131 072 data points to estimate the FEM response over four decades. The FEM response is estimated precisely in the frequency range lower than 100 Hz. However, in the higher range, the real part is strongly affected by the noise. Fig. 8(a) shows the inversion result from the step response time-domain data. The result from the step-off data is also accurate, but the range is limited to lower than 300 Hz. This is because the step-off response has the role of high-cut ®ltering and decreases the high-frequency components, as seen in eqs (3)  and (4). The conventional FFT is not available for the step-off time-domain data sampled within a boxcar sampling window. This is because the step response has a static ®eld that is not zero for t<0, whereas the FFT assumes that the data continue periodically beyond the boxcar window. This produces a large difference from the real step response. Fig. 8(b) shows the result of the impulse response when the Gaussian noise is increased to 10 per cent. In contrast to the result of Fig. 7(a), the frequency response at frequencies higher than 300 Hz becomes inaccurate, although it agrees well with the correct values in the lower range.
In these experiments the smooth spectrum inversion dealt well with the wide range of data over ®ve decades of delay time. The results show that our method provides accurate estimations at low frequencies, i.e. in the near zone.

Field data example
We applied our method to the long-offset transient EM (LOTEM) data that were acquired in the Yurihara oil and gas ®eld in northeastern Japan. This survey was carried out by the Japan National Oil Corporation. Approximately 45 A of electric current was transmitted through an earthed wire with a length of 2.5 km. Vertical induction coil sensors were used to measure the time derivative of the vertical magnetic ®eld at 23 sites. At each site, at least 180 repeat time records were made, with each record consisting of 1024 data points sampled with an equal spacing of Dt=8 ms. In order to test our method, two data sets were selected: one with good data quality and the other with poor data quality. They were recorded through a 50 Hz low-pass ®lter and 50 and 150 Hz notch ®lters to avoid aliasing and to remove arti®cial noise caused by power lines.
To apply and test our method, we resampled the linearly sampled raw data into 31 logarithmically spaced delay-time windows, and an average of the raw data was computed in each window. However, to maintain the causality and the main characteristics of the waveform, we did not average the 10 data for tj0.08 s. Fig. 9 shows the good data obtained at the site offset 3 km from the source and a comparison between the data and the predicted values by the smooth spectrum inversion. We stacked the averages of 159 records, evaluated the ®nal averages in the windows and estimated the relative covariance Figure 8. FEM responses of the model of Fig. 1 estimated from the step response with 1 per cent Gaussian noise (a) and from the impulse response with 10 per cent Gaussian noise (b) using the smooth spectrum inversion. The correct response calculated with the 1-D modelling code EM1D is also shown. The smooth spectrum inversion can still estimate the accurate response at frequencies lower than 300 Hz from the impulse response with 10 per cent independent Gaussian noise.  matrix. Fig. 9(a) also shows the data obtained by stacking the raw data recorded at an equal spacing of 8 ms; it suggests that the quality of data is very good because the stacked data do not have any oscillations, allowing us to recognize the waveform easily. Fig. 10 shows the FEM response estimated from the data in the 31 windows by our inversion method and the result computed from the stacked raw data of 1024 points by conventional FFT. The agreement between them is very good; these two methods provide almost identical results when the quality of data is excellent. In addition, we would like to   emphasize that conventional LOTEM measurements suppress the high-frequency components by using a low-pass ®lter, as the result of Fig. 10 indicates. In practice, the low-pass ®lter is indispensable in recording the impulse response, which rises quickly and then decays slowly. Fig. 9(b) shows the values predicted from the FEM response estimated by the inversion method, which agree very well with the data.
We then applied our method to the poor-quality data shown in Fig. 11(a). The data were acquired at a site offset 8 km from the source using a 20 Hz low-pass ®lter and a 50 Hz notch ®lter. Using 189 records we evaluated the ®nal averages in the windows and also estimated the averages by simple stacking over 189 records. Fig. 12 shows the results of the smooth spectrum inversion and the conventional FFT. The FFT result shows large oscillations in the frequency range higher than 2 Hz and a notch around 0.3 Hz. In contrast, the ABIC method provides a reasonable response with the oscillations and the notch smoothed. This result suggests that our method is useful in obtaining the FEM response for data of poor quality in comparison with the conventional FFT. Fig. 11(b) shows the impulse response data predicted from the resultant FEM response via eqs (5) and (6) compared to the real data. The predicted impulse response data agree with the observed data for t<0.05 s. However, the predicted values appear to be lower than the observed data for later delay times.
Neglecting the off-diagonal elements of the relative covariance matrix C 0d for the data of Fig. 11, we applied the smooth spectrum inversion again (Fig. 13). This approach assumes the noise is uncorrelated. The data ¢t (Fig. 13a) is better than that of Fig. 11(b). In the estimated FEM response (Figs 13b and c), the discrepancy with the result of Fig. 12 becomes visible in the frequencies lower than 0.3 Hz. This approach is not su¤cient because of band-limited noise (Duijndam et al. 1999;Appendix A). Measurement systems are designed so that the noise is band-limited to avoid aliasing and is distributed on the data of other delay times. Consequently, the distribution of the noise becomes a correlated distribution. Fig. 14 shows the correlation coef®cients, r ij =C 0dij /(C 0dii C 0djj ) 1/2 for the data of Fig. 11, which represents the degree of correlation between the noises, and |r ij |j1. For the data of t=0.04 s, the correlation coef®cients of two points ajoining of the peak are large and cannot be neglected. Moreover, at t=2.27 s, the correlation increases with the delay time. This correlation should be included as the off-diagonal elements of C 0d in the inversion. In eq. (17), the observed data, d, are transformed through Wd. Therefore, we should compare the predicted data with the transformed data. Fig. 15 shows the comparison between them for the results shown in Fig. 11. The predicted data are not always lower than the transformed data.

C O N C L U S I O N S
In this paper we have demonstrated the importance of the imaginary component of FEM data in the near zone, and the advantage of TEM measurements in the near zone owing to the absence of the primary ®eld. To estimate correctly the FEM response from the TEM data, we have developed a numerical Fourier transform method using the least-squares approach with a smoothness constraint (smooth spectrum inversion). We have shown that our method can handle data with a wide range of delay times, and that it is capable of analysing noisy TEM data to provide an FEM response without unreasonable oscillations and gaps.
In controlled-source magnetotelluric (CSMT) measurements, Murakami (1988) proposed the time-domain CSMT system and displayed ef®cient and accurate measurements of the near-®eld data by comparison with the conventional frequency-domain CSMT measurements. In light of the recent developments in CSAMT inversion techniques (Routh & Oldenburg 1999 for 1-D; Lu et al. 1999 for 2-D) for analysing the near-®eld CSMT data, our method will be an ef®cient tool in obtaining accurate FEM responses.

A C K N O W L E D G M E N T S
We would like to thank N. Shiga and H. Ishikawa of MINDECO for their efforts in the LOTEM data acquisition, and we also thank Ki Ha Lee for providing his 1-D modelling code EM1D and for his comments, which helped to improve the manuscript. Comments and suggestions from M. J. Unsworth, D. W. Oldenburg, C. G. Farquharson and an anonymous reviewer were very helpful in improving the paper. Figure 15. Comparisons between the transformed data (Wd) and the values predicted by the inverse-cosine transform (a) and the inversesine transform (b) for the results shown in Figs 11 and 12. The weighting matrix Wis di¡erent in the inverse-cosine and the inverse-sine transforms because the datum at t= 0 s is not used in the inverse-sine transform.