1. Introduction
Transient electromagnetics (TEM), a time-domain electromagnetic remote sensing technique, has proven to be an effective approach for near-surface geophysical surveys in the past few decades. Its application scenarios include mineral exploration, engineering, hydrology, and environmental investigations [
1,
2,
3,
4,
5]. With the improvement of the instrument and interpretation, the application scenarios of the TEM approach have been expanded [
6,
7].
Airborne and semi-airborne TEM are two emerging branches of the TEM method. In airborne TEM (ATEM), the transmitter loop and the receiver coils fly together with the aircraft; the ATEM system transmits electromagnetic signals while measuring the earth’s responses in the magnetic field continuously during the flight. In semi-airborne TEM (SATEM), the transmitter is a long wire, or a large closed loop fixed to the ground, whereas the receiver coil flies in the air taking magnetic field measurements at a series of locations. It is important to note that electric field measurements on airborne survey platforms are impossible, so the receiver can only sense the magnetic field via a coil or magnetometer. Because of the airborne platform, ATEM and SATEM have high efficiency and mobility compared to ground TEM. They have been shown particularly useful in rapidly assessing extensive or inaccessible areas.
The efficiency of airborne systems comes at the cost of data quality. The bird in ATEM and SATEM is suspended and towed tens of meters above the surface. Non-constant movement of the aircraft and external forces may cause unpredictable rotation, displacement, and tilt of the bird away from its nominal position and orientation. A large ground loop or long wire is intended to increase the transmitter moment, but long transmitter wires cannot always be deployed along straight lines as planned in practice because of surface obstacles. As a result, the inevitable and sometimes unpredictable deviation of survey configuration in ATEM and SATEM from the nominal parameters can be a significant source of error. This paper refers to such deviation as inexact transmitter and receiver geometry.
An error from inexact geometries is a kind of systematic error or coherent noise, which cannot be easily suppressed by stacking or other statistical methods. In many cases, such errors are considered negligible. However, the demand for high-resolution imaging and quantitative interpretations continues to increase nowadays. Not taking into account the effect of inexact geometry in the inversion may undermine confidence in the interpretation or even produce misleading results. Previous researchers have noticed this problem and proposed solutions. The authors of [
8] developed a novel “agree-to-disagree” strategy to identify and correct misrecorded flight height in ATEM. In [
9], the authors proposed to record the orientation of the ATEM bird and decompose the tilted loop as a linear combination of three orthogonal magnetic dipoles. The authors of [
10] modeled the actual geometry of a large ground loop using a series of electrical dipoles. Generally, it is required to record the actual size, position, and orientation of the transmitter and receiver, and model them by decomposing them into a number of elementary dipoles.
Previous works studied the inexact geometry problem of ATEM and ground large loop TEM and treated them using different dipole source models or data correction methods. Such a problem in SATEM still needs further investigation, as the combination of a grounded long wire source and a drone-borne coil receiver has gained popularity, but the study of such a system is lacking. This paper presents a universal algorithm for both ATEM and SATEM. Our method models an arbitrary 3D source as a weighted sum of small electrical dipoles oriented in the x-, y-, and z-directions; the tilt of the ATEM and SATEM receiver coil is modeled by projecting three-component dH/dt data to the normal direction of the coil; and the forward modeling and sensitivity calculations are carried out using the principle of superposition. Although based on the layered earth model, our improvement advances the technology by providing a unified, convenient, and accurate modeling and inversion framework for ATEM and SATEM data, which are often used by the same group of practitioners but traditionally require separate numerical programs.
In the following, we first set up the physical simulation problem for ATEM and SATEM. Then, the numerical formulas are derived for the semi-analytical solution of the loop source in ATEM and the grounded wire source in SATEM. Using the algorithms, we investigate the effect of inexact geometries in data under different circumstances. In addition, we formulate an inverse problem for ATEM and SATEM based on the Gauss–Newton method and our accurate modeling algorithms. The importance of considering inexact survey geometry in ATEM and SATEM is demonstrated by some synthetic examples. In the end, we show that our methods can make a difference to the inversion and interpretation of field data acquired with inexact survey geometries.
2. Inexact Survey Geometry
Benefitting from the high efficiency and simplicity of 1D modeling, data collected from TEM surveys are frequently forward modeled using a 1D layered earth model. In particular, this kind of model is a sufficient representation of the sedimentary environment and is favored in hydrological and other near-surface investigations. Our development in this work is for the modeling and inversion of the layered earth model in a right-handed Cartesian coordinate system (
Figure 1).
In this paper, we study two airborne variants of TEM methods—ATEM and SATEM surveys. For ATEM, we pay particular attention to systems towed by a helicopter, for example, the SkyTEM and VTEM systems. These systems feature a bird of rigid frame for mounting the transmitter loop and receiver coil (
Figure 2). Such a design has been proven advantageous in near-surface and environmental applications that demand high-resolution imaging. For SATEM, we consider a current-carrying transmitter wire with its two ends grounded on the surface through a pair of electrodes. Once the ground source is set up, the aircraft tows the receiver coil and flies above the survey area (
Figure 2). The transmitters in ATEM and SATEM are sometimes termed the inductive source and galvanic source, respectively; in both cases, the receiver is a horizontal coil that measures the vertical component of the magnetic field or its time derivative.
The configurations in
Figure 2 are nominal. In practice, ATEM and SATEM always work with inexact survey geometries. Due to various factors such as terrain and wind, the rigid frame of ATEM cannot be exactly leveled as designed (
Figure 3a). Modern ATEM systems have been equipped with measuring devices for displacement, tilt angles, and frame clearance. Although this information is often available and delivered to the clients, it is ignored in most cases. SATEM also has the same inexact geometry problem for its receiver coil. In addition, its transmitter wire laid on the surface cannot be guaranteed to be straight when it has to bypass obstacles such as buildings, lakes, roads, etc. (
Figure 3b). As a result, the sources of systematic errors in SATEM data can be crooked transmitter wires and tilted receiver coils.
Fortunately, modern technology such as differential GPS, LiDAR, and the tiltmeter has enabled continuous monitoring of actual survey geometries for airborne systems. The actual way path of the SATEM transmitter wire can also be precisely logged by RTK or other surveying methods. This paper intends to utilize this positioning information to accurately model ATEM and SATEM data and subsequent high-resolution subsurface imaging. Here, we propose a universal solution for both ATEM and SATEM. An arbitrary transmitter configuration can be described by a series of waypoints in 3D space and can be simulated using a set of moment-scaled 3D unit electrical dipoles. Each electrical dipole is further decomposed to a linear combination of
x-,
y-, and
z-oriented dipole sources, which has semi-analytical solutions for layered earth.
Figure 3 provides a conceptual diagram using a set of end-to-end electric dipoles to simulate actual transmitter geometry. We note that, in general, the transmitter wire path must be divided finely enough that the distance between the receiver and the dipole is greater than five times the length of the segment [
11]. Once the transmitter is decomposed, the original modeling problem becomes many parallel subproblems, each of which is excited by an electrical dipole source in either
x-,
y-, or
z-direction. By using the principle of superposition, the modeled data or sensitivity are the sum of these dipole solutions, with the orthogonally projected dipole moments as weights.
The inexact geometry of an airborne receiver is mostly from its rotation or tilt. This paper assumes the receiver coil may have an arbitrary orientation specified by a zenith angle and an azimuthal angle of its normal direction. Our forward-modeled data are always three-component in the x-, y-, and z-directions; then, the desired data (or sensitivity) are obtained by projecting the modeled vector field data to the actual direction.
3. Forward Modeling and Analysis
The core of our modeling method is the computation of the time-domain magnetic field (or its time derivative) responses of layered earth excited by an
x-,
y-, or
z-oriented electrical dipole source on or above the surface (
Figure 3). We first obtained the three-component magnetic field data calculated in the frequency domain. Then, they were frequency-time transformed to the time-domain magnetic field or its time derivative. Finally, convolution with the real transmitter waveform was carried out to obtain the time-domain solution of a specific ATEM or SATEM system.
For a layered earth model, the frequency-domain solution of the magnetic field induced by an
x-oriented electric dipole (XED) located at (0, 0, −
h) is [
12]
where the subscript of the magnetic field presents the field component, and the superscript of the magnetic field denotes the dipole direction. (
x,
y,
z) are the coordinates of the receiver.
h is the height of the dipole above the surface.
ρ is the distance between the receiver and dipole.
Ids is the transmitter moment.
J0 and
J1 are the Bessel functions of the zeroth and first order, respectively.
u0 is the wavenumber of the air space.
λ is the integration variable of the Hankel transformation.
rTE is the reflection coefficient for TE mode, and
rTM is the reflection coefficient for TM mode.
The reflection coefficient is defined as
And
with
where
n is the index of the layer.
un is the wavenumber.
hn is the thickness.
i is the imaginary unit.
ω is the angular frequency and
f is the frequency.
kn is the vertical wavenumber.
μ0 and
μn are the magnetic permeability of free space and the
nth layer, respectively.
The magnetic field induced by a
y-oriented electric dipole (YED) located at (0, 0, −
h) is [
12]
The magnetic field induced by a
z-oriented electric dipole (ZED) located at (0, 0, −
h) is unavailable from the publicly available literature, so we directly derived them using the Schelkunoff potentials method [
11] as follows:
Once the frequency-domain responses were obtained, a cosine transform was employed to calculate the time-domain magnetic field and its time derivative induced by a step-off waveform [
13].
For efficient computing, we utilized a digital linear filter to calculate the Hankel transform and cosine transform involved in the equations (1, 5, 6, 7). Details about the digital linear filter can be found in a number of publications [
14]. When the transmitter waveform is not step off, the program utilizes the convolution of the response induced by a step-off waveform and the derivative of the arbitrary waveform to calculate the response [
15].
Our first modeling example concerned a typical central-loop ATEM configuration with a transmitter moment of about 531 A∙m
2 and a step-off waveform (
Figure 3a). The circular transmitter was approximated by 40 electric dipoles. The zenith angle of the frame varied from 5° to 20°, and the azimuth angle was set to 90°. The flight height was set to 10 m. We compared our accurate modeling with the common practice of simulation using the analytical formula of the central loop configuration for the vertical magnetic field, assuming that the transmitter was a standard circle and was perfectly leveled (the zenith angle is 0°).
Figure 4 presents the responses over a 0.01 S/m half-space. If actual survey geometry is not modeled, the magnitude of data can be reduced by 0.03~11%.
Figure 4b also demonstrates that a larger zenith angle, corresponding to a more seriously tilted frame, results in a more pronounced mismatch between the magnitudes of the actual case and the nominal situation. We also found that the discrepancy is most significant in the range of commonly measured time channels, and the late-time responses are more affected.
In our second modeling example, we simulated a SATEM survey on a three-layered earth model with conductivity [0.01, 0.1, 0.001] S/m and thickness [100, 100, ∞] m. The planned transmitter was a 1000 m grounded long straight wire, and the waveform was step off. The receiver was located near the center of the source wire, with a 200 m offset and a height of 10 m above the surface. In practice, the crooked transmitter was simulated by adding random disturbances to the waypoints of the wire; the maximum deviation was not more than 100 m (see
Figure 3b for a plan view). The receiver’s zenith and azimuth angles were 10° and 45°, respectively. Errors due to inexact geometries were more complex and more significant for SATEM surveys. If only the transmitter wire was crooked, but the receiver coil was leveled, a large discrepancy could be seen at early times (
Figure 5a,b); the errors from a tilted receiver when the transmitter wire was straight were much less significant, but could still reach more than 10% (
Figure 5c,d); the compound effect of inexact transmitter and receiver geometries was mostly dominated by the transmitter, with a maximum error of 300% at early times (
Figure 5e,f). In field surveys, it is preferred to place the receiver as close to the source as possible to maximize the received secondary field excited by the transmitter. Our simulation showed that the effect of inexact survey geometries in SATEM can be very strong, especially when the transmitter wire cannot be laid out according to the original plan, but the nominal position is used in the modeling. Proper care is required in the modeling and inversion to prevent the bias from propagating to the inversion models and the geological interpretation.
4. Inversion
In this section, we construct the inverse problem based on the forward modeling approach proposed in the previous section. Using the inversion, we further investigate how the effect of inexact survey geometry in ATEM and SATEM can bias the inverted models. In this work, we employ a classical objective function
where the first term of the right-hand side (RHS) is the data misfit, and the second term is the model norm, which is introduced to stabilize the inverse problem and to add constraints. Specifically,
d represents the measured data,
m is the vector of model parameters, and
F is the nonlinear forward modeling operator presented in
Section 3.
Wd is the data weighting matrix, which provides the estimation of data uncertainty. If the data error is Gaussian,
Wd can be a diagonal matrix whose elements are the reciprocals of the standard deviations of each datum.
Wm denotes the model weighting matrix, which can be a thickness weighting operator, first-order difference operator, or focusing operator [
16]. A combination of the thickness weighting operator and the first-order difference operator is employed in our implementation.
µ is the regularization parameter, balancing the relative weighting between the data misfit and the model norm during the inversion process. In Equation (8), we use an
L2 norm to evaluate the data misfit and model norm, but other forms, such as
L1 norm and
Lp norm, are alternatives.
We use a Gauss–Newton iterative method to minimize the Equation (8). Because the forward modeling operator
F is nonlinear with respect to
m, a Taylor series expansion is utilized to linearize the
F(
m) locally as
where the
mk+1 and
mk are the inverted models at
k + 1 and
kth iteration, respectively.
Jk represents the Jacobian matrix, the derivative of
d with respect to
mk. The iterative equation is found by substituting the Equation (9) into Equation (8), differentiating the Equation (8) with respect to
mk+1, and equating it to zeros. Finally, we get the following iterative equation
The regularization parameter plays an essential role in balancing the data misfit and model norm in the inversion process. Here, we use the following equation to calculate the initial regularization parameter [
17]
where
smax (
A) denotes the largest singular value of matrix
A. In the iterations, a cooling strategy is adopted to update the regularization parameter
µ. The iterative process continues until the data misfit has decreased to a specified target (e.g., when the data are assumed to be Gaussian and uncorrelated, an appropriate target is the number of measured data). Alternatively, the inversion can be terminated when the gradient of the objective function falls below a given tolerance [
18]. During the inversion, we employed an analytical approach to calculate the Jacobian matrix, which characterizes the sensitivity of the EM fields to changes in layer conductivity [
19,
20,
21].
The first ATEM synthetic inversion example employed a circular transmitter with an ideal step-off waveform over a three-layered earth model. The receiver was located in the center of the loop, measuring the z component of the dH/dt data. The conductivity and thickness of each layer were set as [0.01, 0.002, 0.01] S/m and [100, 200, ∞] m, respectively. The transmitter moment was 314.16 A·m2; the frame’s height was 10 m, and there were 20 time channels from 10−5 s to 3 × 10−3 s. In the actual geometry, the zenith and azimuth angles of the receiver coil were 20° and 90°, respectively.
Synthetic data were generated by forward modeling with the tilted ATEM frame simulated by our electrical dipole methods. Then the data were inverted by assuming the frame was still perfectly leveled and accurately modeling the actual geometries. Our synthetic test showed that although both inversions achieved a similar level of data fitting and produced a similar structure in the models, noticeable discrepancies in the models were observed in the near-surface and deep layers. Recalling the forward modeling response analysis in
Section 3, the discrepancies in both shallow and deep can be explained as the effect of inexact geometry existing in all time channels from early to late. Model bias due to survey geometry does not change qualitative interpretation but may be a problem when high-accuracy results are required in near-surface applications.
The consequences of not modeling actual survey geometry are more prominent for SATEM. Our synthetic test found that inversion with an accurate transmitter wire path and receiver coil orientation can greatly improve the characterization of near-surface structure and enhance the clarity of deep structure compared to an inversion case in which the survey geometry follows the nominal parameters. Although the data fit equally well, the large discrepancy of the models in
Figure 6c reminds us that recording and modeling actual survey geometry can be a critical issue in the processing and interpreting of SATEM data.
Finally, we tested our method in a realistic scenario. This example aimed to demonstrate inversion involving a survey geometry deviating from the nominal. The data were from an ATEM survey of the SkyTEM system collected in 2006 over the Bookpurnong irrigation district in southern Australia. The Murray River and adjacent floodplains in the Bookpurnong area have been extensively salinized. This survey was carried out to map the salinity of the groundwater. We refer readers to [
22] for a complete background on the geology and hydrogeology of the study area. The dataset is publicly available in [
23]; previous processing and inversions can be found in a number of publications [
24,
25,
26], in which the tilt of the frame was not considered.
The SkyTEM operates in high-moment and low-moment modes that are supposed to detect both the deep and shallow. Here we used the low-moment data, which operated at a peak current of 40 A and a base frequency of 222.5 Hz [
27]. The survey flowed at an altitude of about 60 m [
25]. For the purpose of comparison, we chose one sounding at location (462,585 m, 6,200,782 m), and performed two inversions: one with the actual frame tilt recorded during the survey and the other with the nominal leveled-frame geometry (
Figure 7). At that sounding, the zenith and azimuth angles of the frame were about 25° and 0°, respectively. In general, there was good agreement in the two inversion models, i.e., both of them supported a resistive layer at a depth of about 50 m. However, the near-surface conductivities recovered in the two inversions were evidently different. Not considering the actual tilt of the ATEM frame tended to underestimate the conductivity, and this feature was consistent with that of the synthetic data example. Such a discrepancy in conductivity may lead to non-negligible errors in the interpretation of groundwater salinity.
SATEM stands out for its cost effectiveness and efficiency and has received more attention in recent years; however, most SATEM systems are still in the experimental stage. To the best of our knowledge, publicly available, reliable SATEM field data are not accessible. Furthermore, in many SATEM surveys, the importance of transmitter and receiver geometry is not realized, and only nominal layouts are reported. Therefore, inversion of SATEM field data taking into account the geometry of the transmitter and receiver is absent in this work. Nonetheless, our work provides theoretical guidance for future SATEM surveys.