Generation of synthetic wide-band electromagnetic time series

The estimation of the earth transfer functions in MT prospecting method poses the greatest difficulty. As in the seismic prospecting method, this task requires the development of advanced processing techniques. In order to assess the performance of each technique, controlled synthetic data and different noise types, which simulate the observed signals, are required. This paper presents a procedure to generate a wide-band noise-free electromagnetic field to be used both for magnetotelluric and audio-magnetotelluric studies. Furthermore, an effort was made to extend the simulation procedures to generally stratified and simple inhomogeneous earth structures. The discretetime magnetic field values are generated through the inverse Fourier transform of a continuous amplitude spectrum and a sampling procedure. The electric field time series are obtained by the convolution of the magnetic field time series, calculated in the interested frequency band, with a non-causal impedance impulse response. Polarized fields, which are important when inhomogeneous media are considered, are also generated.


Introduction
Computer simulation of electromagnetic field components has frequently been used for testing estimation techniques of the magnetotelluric impedance tensor (Goubau et al., 1978;Mc-Mechan and Barrodale, 1985;Yee et al., 1988;Larsen et al., 1996) and for assessing preprocessing methods of extraction of the stationary and coherent part of signals corrupted by noise (San Filipo and Hohmann, 1983;Lamarque, 1999).The testing and assessing tasks are generally accomplished in the frequency domain through the power spectra of the measured signals.Hence, the simulation is carried out in the frequency domain.
Some authors have pointed out the numerous problems arising in the frequency-domain approach and have proposed estimation methods in the time-domain (Kunetz, 1972;Ernst, 1981;McMechan and Barrodale, 1985;Yee et al., 1988;Spagnolini, 1994).Moreover, for some purposes (e.g., noise identification and removal, filling missing data, search for outliers) the time-domain methods are also useful in a frequency-domain context (Egbert, 1992;Egbert et al., 1992).Therefore, these considerations justify the interest in generating synthetic time series for simulation purposes.
In the above-mentioned papers, the synthetic data are constructed by means of different techniques.Larsen et al. (1996) used measured magnetic time series and an estimated 1D transfer function, modified by a distortion function, to generate the electric time series in a test area.Simple pseudo-random numbers were used to represent the noise-free random components of the time series (Lamarque, 1999) or the real and imaginary parts of the incident magnetic field components (Goubau et al., 1978).The latter authors calculated the electric field spectra using a bi-dimensional impedance tensor made by simple complex relationships.San Filipo and Hohmann (1983) and Yee et al. (1988) assume a function which tends to zero exponentially with increasing values of frequency to represent the main spectral characteristics of the natural magnetic field, simulating this behaviour through recursive digital filters.These same authors generated the magnetic time series by convolving a pseudo-random number sequence with the filter coefficients.Yee et al. (1988) computed the noise-free electric field components by convolving in time the magnetic time series with the impedance impulse responses modelled as rational forms.
Uncorrelated noise is simulated by adding random numbers either to the spectra (Goubau et al., 1978) or to the time series (Yee et al., 1988;Larsen et al., 1996).Some authors introduce the noise by adding sinusoids to the signal.These are characterized by frequencies and amplitudes which are typical of industrial and atmospheric noise (Lamarque, 1999).
In recent years, the English version of a systematic study on the MT synthetic data generation, carried out by Russian researchers, was published (Varentsov and Sokolova, 1995).The study was conducted within the framework of the project to compare MT data processing techniques, using synthetic data sets (COMDAT Project).The discrete values of the spatial components of the magnetic field spectrum, considered in the MT frequency band, were constructed by opportunely multiplying selected random functions with a simple realistic field structure (the form was used with < 0).The electric field spectra were calculated via arbitrary one-dimensional scalar impedances.Finally, the time series were obtained by a discrete inverse Fourier transform of the spectral data approximated by a piecewise constant function.Noise components were added both in the frequency and in the time domain.
The present paper describes an approach to generate wide-band electromagnetic time series.The procedure adopted is represented in a flow chart (fig.1).The aim was to simulate real acquisition techniques to be used both in magnetotelluric and audio-magnetotelluric studies.In this first work only noisefree signals have been considered.The results were then applied to 1D and simple non 1D structures.

EM field simulation
The noise-free horizontal electromagnetic field components are related by means of the earth transfer functions in a two-input/twooutput linear system (Cantwell, 1960).In the time-domain, the relation between the timevarying input magnetic field h(t) and the timevarying output electric field e(t) can be written in terms of convolution equations among the field components where z xx , z xy , z yx , z yy are the elements of the impedance impulse response tensor.The input and output channels in the system are perturbed by added noise (Swift, 1967;Sims et al., 1971)  For application purposes, the continuous-time convolutions of eq.(2.1) must be replaced by discrete forms.This task can be accomplished by applying the well-known classical sampling theorem (Freeman, 1965) to lowpass (anti-alias) filtered continuous signals.
In the following, the procedure adopted to Construction of a continuous function representing the magnetic field spectrum H(ω ) in the 10 -4 -10 4 Hz frequency range (eq.(2.5) and Fig. 2) Construction of the analytical expression for the frequency domain earth impedance Z (ω) for a stratified model (eq.(2.9))

Input a sequence of Random Numbers
Noise-Free Magnetic Field Time Series Noise-Free Electric Field Time Series DFT DFT H(∆ω ) Fig. 4 E(∆ω ) Fig. 7 generate the noise-free electromagnetic field components and the impedance impulse response is illustrated.

H-field generation
The natural EM field is due to a variety of causes.At frequencies above 6 Hz, roughly, it is primarily due to spherics, which are EM transients generated by lightning strokes.Below 6 Hz, the field is of geomagnetic and ionospheric origin.The general features of the natural field in different frequency bands have been illustrated by many authors (Campbell, 1967;Filloux, 1973;Serson, 1973;Mcnae et al., 1984;Labson et al., 1985;Malergue et al., 1986).
The horizontal components of the field present average amplitudes, generally increasing with longer periods and reaching a value of about 10 nT at 10 4 -10 5 Hz.Moreover, two distinct minima in the 1 Hz to 6 Hz and 1 kHz to 2 kHz frequency ranges are observed.These minima represent the separation between ionospheric pulsations and spherics and between lowfrequency and high-frequency spherics, respectively.
A smoothed, mean spectrum of the magnetic field is shown in fig. 2. To obtain a relationship which describes the field, the plotted behaviour can be modelled through a linear time invariant system.
Depending on the goal of the study (geomagnetic,magnetotelluric, audiomagnetotellurics prospecting), the wide-band field is used in a limited frequency range.This operating characteristic influences the procedure for the generation of the magnetic field time series.In fact, in accordance with some authors (San Filipo and Hohmann, 1983;Yee et al., 1988), the approximation of the frequency characteristics of the field could be performed directly in the discrete domain by writing the transfer function of the system in terms of poles and zeros in the z domain.An alternative is to transform the transfer function from the continuous to the discrete domain through a bilinear transformation (Oppenheim and Schafer, 1975).In the present study, these methods cannot be adopted for the following reasons: 1) For each sub-band, the power of the H field at frequencies greater than the Nyquist one is not negligible, thus resulting in the aliasing problem when the normal z transform is used.
2) With regards to the bilinear z, this method maps the complex continuous s plane into the complex sampled z plane without aliasing, but with a frequency distortion (warping) (Oppenheim and Schafer, 1975).This distortion, resulting from the non-linear character of the transformation, is compensated when the frequency-magnitude curve of the system can be approximated by a piecewise constant function in which only two or three characteristic frequencies need to be transformed without distortion (prewarping procedure) (Golden and Kaiser, 1964).This requirement is not fulfilled in the simulation of a physical field which is characterized by a complex morphology spectrum.
Considering the previous arguments against the use of a discrete-time relationship, a continuous function for the H field was constructed.The sampled impulse response was obtained by a windowed integral Fourier transform, computed with the time interval dictated by the Nyquist frequency of the studied band.
In particular, by using the poles and zeros form in the continuous (Laplace) domain, a system function can be expressed (Seshu and Balabanian, 1959) as where, s = σ + iω is the Laplace complex variable with σ the damping factor of the system and ω the frequency.The complex numbers s 0k are the m zeros while the s ph are the n poles of the function, respectively; K is a scale factor.
The analytic properties of the function are determined by its poles and zeros.In order for the system to be stable, the poles must be in the left half of the complex s-plane and any pole on the imaginary iω axis must be simple.Moreover, since the network function must be real for a real variable s, it follows that the poles and zeros must either be real or in complex conjugate pairs.
The simple zeros and poles and the complex conjugate pairs in the numerator and denominator of (2.3) can be factored in the following form (Gatti et al., 1966) (2.4) where the indexes i, j, l, k indicate the generic simple zero, simple pole, complex conjugate zero and complex conjugate pole, respectively.The scale factor K 1 and the correction factors ζ l and χ k are defined by the following equations: If s = iω is inserted into eq.(2.4), the magnetic field relationship in the frequency domain can be obtained.The frequency and characteristics of the poles and zeros are selected by considering that simple poles and zeros introduce a 20 dB/decade slope change in the function curve on a cartesian decibel versus the log (ω) representation (Bode's plot), while the complex conjugate pairs double these changes.Moreover, the ζ l and χ k quantities act as damping factors at each critical frequency.
Considering the magnetic field spectrum shown in fig.2, the downward change of the slope (energy decrease) can be modelled by poles, while the upward changes (energy increase) by zeros.In particular, four poles and four zeros, with values as shown in table I, were sufficient to represent the changes in the spectrum.The simple poles and zeros were chosen to be real, while for the complex conjugate pairs, correction factors ζ l = 1 and χ k = 1 were adopted.These conditions are equivalent in considering real double poles and double zeros.For the sake of system stability, all the poles have negative values.Finally, the scale factor K 1 was purposely selected in order to obtain a low frequency asymptote of 10 nT for the field spectrum.
With the above conditions eq. (2.4) becomes (2.5) The continuous-time impulse response of the EM field can be obtained from the inverse Fourier transform of eq.(2.5) If ω l and ω h are the low-frequency and highfrequency cut-off of the band, and H r (ω) and H i (ω) the real and imaginary parts of the complex field H(ω), respectively, then eq.(2.6) becomes (2.7) If ∆t is the sampling period and T the maximum acquisition time, being H r (ω) and H i (ω) even and odd functions, respectively, the sampled impulse response can be written as (2.8)Both for the numerical testing of the procedure, discussed at the end of Section 2.3, and for the generation of polarized fields in Section 3, the magnetic field in the frequency range 0.0005-5 Hz has been computed.The latter range is adequate for investigating the earth models considered.The resulting sampling period and maximum period were 0.1 s and 2000 s, respectively.
The sampled magnetic field impulse response, calculated from eq. (2.8), is shown in fig.3a.The Discrete Fourier Transform (DFT) of the impulse response (fig.3b) is in perfect agreement with the continuous magnetic field amplitude spectrum, within the frequency range of interest.
In order to generate a magnetic field time series, the impulse response of fig.3a was convolved with a sequence of pseudo-random numbers with a zero mean and unity-variance.The DFT of the resulting time series is shown in fig. 4.
The generation of the electric field could be performed in the frequency domain as the product of H(ω) and Z(ω).However, if eq. (2.5) for H(ω) and, for example, eq.(2.9) for Z(ω) are considered, their product is a very complex function; thus, its inverse Fourier transform is performed with difficulty.Hence, it was decided to inverse transform the two functions separately, thereby obtaining the electric field by a convolution procedure in the time domain.
The determination of the discrete-time impedance impulse response which produces the same outputs as the continuous-time counterpart remains a problem which is thoroughly and elegantly discussed by Egbert (1992).As the continuous-time impedance impulse response is unbounded, the non-causality of the corresponding discrete form generally follows.Hence, the discrete-time impulse response can be estimated from the frequency-domain impedance by band limiting and sampling it, as well as, by taking into consideration the negative time values.This procedure, applied to a stratified model, is presented in the following paragraph.

Impedance impulse response
The analytical expression of the frequency domain MT impedance, relating magnetic and electric field components, can be written for a stratified earth as (Kunetz, 1972) (2.9)  where the first term on the right hand side represents the impedance over a half space of constant resistivity 1 , the second one results from the reflections of the q m images on the boundaries between the layers and , h 1 being the thickness of the first layer (see Kunetz, 1972, for the determination of the images from the reflection coefficients).
The corresponding continuous time impulse response can be written as (2.10)If the low-frequency and high-frequency limits of the band are indicated, as previously, with l and h , respectively, then eq.(2.10) becomes (2.11)By integrating by parts, eq.(2.11) yields (2.12)By placing l = 2 /T and h = / t, with t and T the sampling period and the longest period, respectively, the solution of the integrals on the right side of (2.12) leads to the following algorithm for the discrete-time impedance impulse response for a stratified earth: (2.13) The above expression can be easily calculated numerically for a given sequence of resistivities and thicknesses.Figure 5 shows an example of the impulse response for an H-type three-layer earth (the model is shown in the same figure).
The impulse response shows a structure which is mainly localized around the origin, with oscillating decreasing values.Moreover, a different behaviour between positive and negative time lags is observed.
The procedure previously adopted for layered earth models could be extended and applied to more complex situations, the analytical solutions of which are available, e.g., vertical discontinuity (d'Erceville and Kunetz, 1962), vertical dyke (Rankin, 1962), anisotropic stratified sequence (Douglas et al., 1967).

E-field generation
The electric field time series is calculated by convolving in time the magnetic field time series with the earth impulse response (eq.(2.13)).The procedure adopted for the convolution operation requires further discussion.
If M values for h( t) and N values for z( t) are considered, with N M, the e( t) resulting from the convolution will contain M + N -1 values.However, there are (N -1) values both at the beginning and at the end of the convolved series which are not fully constrained by the data.A solution to this problem is to strip off these edge regions from the data, thus obtaining M -N + 1 e( t) values.Moreover, in order to generate magnetic and electric field time series the same length as well as a causal electric field, the (N -1)/2 values at the and at the end of h( t) must be discarded.Thus, the value e(N) will be generated by h[(N + 1)/2].
Figures 6 and 7 show the impulse response and the spectral amplitudes of the noise-free electric field for the 1D structure of fig. 5.For the generation of the electric field impulse response 20 000 samples of the magnetic field impulse response, together with 4001 samples of the impedance impulse response were used.
For the purpose of numerically testing the procedure, the impedance transfer function Z( ), obtained as E( )/H( ) for a single time series, has been computed.The computed discrete impedance transfer function (fig.8) fits the theoretical curve well in all but low frequency values, due to the number of samples used.The impulse response shows a structure which is mainly localized around the origin, with oscillating decreasing values.Moreover, a different behaviour between positive and negative time lags is observed.
The procedure previously adopted for layered earth models could be extended and applied to more complex situations, the analytical solutions of which are available, e.g., vertical discontinuity (d'Erceville and Kunetz, 1962), vertical dyke (Rankin, 1962), anisotropic stratified sequence (Douglas et al., 1967).

E-field generation
The electric field time series is calculated by convolving in time the magnetic field time series with the earth impulse response (eq.(2.13)).The procedure adopted for the convolution operation requires further discussion.
If M values for h( t) and N values for z( t) are considered, with N M, the e( t) resulting from the convolution will contain M + N -1 values.However, there are (N -1) values both at the beginning and at the end of the convolved series which are not fully constrained by the data.A solution to this problem is to strip off these edge regions from the data, thus obtaining M -N + 1 e( t) values.Moreover, in order to generate magnetic and electric field time series of the same length as well as a causal electric field, the (N -1)/2 values at the beginning and at the end of h( t) must be discarded.Thus, the value e(N) will be generated by h[(N + 1)/2].
Figures 6 and 7 show the impulse response and the spectral amplitudes of the noise-free electric field for the 1D structure of fig. 5.For the generation of the electric field impulse response 20 000 samples of the magnetic field impulse response, together with 4001 samples of the impedance impulse response were used.
For the purpose of numerically testing the procedure, the impedance transfer function Z( ), obtained as E( )/H( ) for a single time series, has been computed.The computed discrete impedance transfer function (fig.8) fits the theoretical curve well in all but low frequency values, due to the number of samples used.

Field polarization
The fields simulated through the previously described procedure were generated by considering a horizontally stratified isotropic half-space.In this case, it is well-known that the polarization of the waves is not important (Porstendorfer, 1975).On the contrary, in an earth having a complicated electrical structure, the measured electric and magnetic fields depend strongly on the orientation of the primary electric and magnetic fields.The effect of an electrical inhomogeneity varies as the orientation of the primary field varies.Thus, the total field polarization is a function of the current sources in the ionosphere.Since, in general, the location of the sources changes with time, or several sources may contribute simultaneously to the field, the total electromagnetic field observed along arbitrary directions often changes with time at a given frequency.
The representation of the impedance as a tensor instead of as a scalar quantity was developed to take into account the orientation of the primary field and the structural complexities.
In considering the most general polarization of the magnetic field, that is the elliptical one, an electric field, itself elliptically polarized (Porstendorfer, 1975), results.
A 2D structure is now considered and measuring axes which coincide with the principal axes of the anisotropy tensor, are assumed.Remembering that any elliptically polarized magnetic wave can be synthesized from two linearly polarized waves, the field may be divided into two parts, parallel and perpendicular to the strike of the structure, respectively.In this particular case, each of the electric field components relates to the perpendicular magnetic field component alone.Thus, the magnetotelluric eqs.(2.1) become where z xy and z yx are the impedance impulse responses of two different 1D models.In summary, to preserve the real field behaviour it must be reminded that it presents (Grillot, 1975): 1) An impulsive nature, occurring in short bursts.
2) A distinct, randomly variable, polarization in each burst.
Together with the previous characteristics of the fields, the field components must be correlated to each other.
In order to satisfy the previous requirements, the magnetic time series have been modulated with a pseudo-random sequence, uniformly distributed in the interval (-1,1).This modified time series was assumed to represent one component of the field.The other component was obtained by modulating the same initial field with a sequence of numbers, complementary to the pseudo-random ones.
As earth model, a stratified anisotropic structure, characterized by the two 1D models shown in table II has been selected.The electric field components were generated by the procedure discussed in the previous sections.
In accordance with Grillot (1975), the occurrence of the field characteristics was tested through a band-pass filtering of the time series.Figure 9 shows a 1600 s specimen of the horizontal components of the field time series filtered with a band-pass centered at 5 s and a bandwidthcentre frequency ratio of 0.3.In the figure the polarization ellipses of the total horizontal field vectors for short sections of the data have also been plotted.
The figure shows ellipses well defined over 2-3 cycles with polarizations changing from section to section.As expected, the magnetic and electric ellipses are perpendicular.
Considering independent sets of field components at a given frequency, which occurs when the sections have different polarizations, is the basis for any statistical estimation of the earth impedances.

Conclusions
A procedure to generate synthetic wide-band electromagnetic time series has been developed and described in detail.It attempts to take into account and to reproduce the main characteristics of the actual electromagnetic field.The main intent was to carry out a simple yet accurate MT experiment simulation, consisting in the interaction of a synthetically generated magnetic field with assigned resistivity models of the underground.These models were complex with regards to their analytical representation.The electric field time series were calculated in the timedomain through a non-causal impedance impulse response.The procedure was extended to simple inhomogeneous structures considering the field polarizations.
Further studies will include: 1) More complex earth structures.
2) The effects of different noise types.
3) The relative performances of different statistical approaches for a reliable impedance estimation.

Fig. 1 .
Fig. 1.Flow chart outlining the procedure for generating the electromagnetic field time series.

Fig. 2 .
Fig. 2. Smoothed, mean amplitude of the horizontal electromagnetic field component variations simulated through a linear time invariant system.
Fig. 3a,b.a) Impulse response of the electromagnetic field in the frequency range 0.0005-5 Hz, sampled every 100 ms; b) theoretical amplitude spectrum of the magnetic field in the frequency range considered (continuous line) and DFT of the impulse response in a) (dashed line).

Fig. 4 .
Fig. 4. Amplitude spectrum of the magnetic time series generated by the convolution of the impulse response of fig.3a with a sequence of pseudo-random numbers.

Fig. 5 .
Fig. 5. Non-causal discrete-time impedance impulse response for a three-layer earth model (shown in the upper part of the

Fig. 6 .
Fig. 6.Electric field impulse response resulting from the convolution in time of the magnetic field impulse response (fig.3a) with the non-causal impedance impulse response of fig. 5.

Fig. 7 .
Fig. 7. Amplitude spectrum of the electric field generated by the magnetic field time series.

Fig. 6 .
Fig. 6.Electric field impulse response resulting from the convolution in time of the magnetic field impulse response (fig.3a) with the non-causal impedance impulse response of fig. 5.

Fig. 7 .
Fig. 7. Amplitude spectrum of the electric field generated by the magnetic field time series.

Fig. 8 .
Fig. 8. Fit of the impedance transfer function generated by synthetic data (solid line) to the theoretical curve (dashed line) for the three-layer model in fig. 5.

Fig. 9 .
Fig. 9. Horizontal components of the time series filtered at 5 s and total field polarization of short sections.The ellipses are centered at each of the sections and are in the same scale as the corresponding field.

Table I .
Characteristics of the poles and zeros of the system function used in the generation of the magnetic field variations.All poles and zeros are real and shown as absolute values.

Table II .
Two 1D layered earth sections representing a weakly anisotropic stratified structure, used to generate the impedance impulse response along the anisotropy directions.