There is a clear demand to increase detection depths in the context of raw material exploration programs. Semi-airborne electromagnetic (semi-AEM) methods can address these demands by combining the advantages of powerful transmitters deployed on the ground with efficient helicopter-borne mapping of the magnetic field response in the air. The penetration depth can exceed those of classic airborne EM systems because low frequencies and large transmitter-receiver offsets can be realized in practice. A novel system has been developed that combines high-moment horizontal electric bipole transmitters on the ground with low-noise three-axis induction coil magnetometers, a three-axis fluxgate magnetometer, and a laser gyro inertial measurement unit integrated within a helicopter-towed airborne platform. The attitude data are used to correct the time series for motional noise and subsequently to rotate into an earth-fixed reference frame. In a second processing step, and as opposed to existing semi-AEM systems, we transform the data into the frequency domain and estimate the complex-valued transfer functions between the received magnetic field components and the synchronously recorded injection current by regression analysis. This approach is similar to the procedure used in controlled-source EM. For typical source bipole moments of 20–40 kAm and for rectangular current waveforms with a fundamental frequency of approximately 10 Hz, we can estimate reliable three-component (3C) transfer functions in the frequency range from 10 to 5000 Hz over a measurement area of 4×5  km2 for a single source installation. The system has the potential to be used for focused exploration of deep targets.

Airborne electromagnetic (AEM) methods measure the electrical conductivity distribution within the subsurface and are widely used for mineral exploration and groundwater studies (Siemon et al., 2009; Smith, 2014). AEM surveys can cover large areas in a short amount of time and are much more efficient than ground-based measurements. Furthermore, they can be performed in rough or even inaccessible terrains. With the exception of the audio-frequency magnetic (AFMAG) technique (Ward, 1959), which uses natural magnetic field variations, AEM systems use man-made transmitters (Fraser, 1972; Thomson et al., 2007). Because the power of an AEM transmitter is practically limited when carried by an aircraft, the penetration depth of airborne EM is also limited, typically to a few hundred meters. But, modern exploration targets for airborne geophysics require exploration depths in excess of 1000 m.

The semi-airborne EM (semi-AEM) approach uses ground-based transmitters and air-towed magnetic field receivers, as illustrated in Figure 1. This approach facilitates the use of very powerful transmitters on the ground, whereas the advantages of survey efficiency and spatial coverage of AEM can be partially maintained, when a passive sensor system is towed by an aircraft over the footprint of the transmitter. Semi-AEM systems of various types have been developed in the past and include the TURAIR (Associates Australasia Pty. Ltd., 1971; Becker, 1979), FLAIRTEM (Elliott, 1996, 1998), GREATEM (Mogi et al., 1998), TerraAir (Smith et al., 2001), and HeliSAM (GAP Geophysics, Cattach and Boggs, 2005) systems. These systems differ in the type of transmitter (grounded electric bipole or loop), the magnetic-field sensor (air coil, induction coil, and total field magnetometer), and the data evaluation domain (time or frequency domain).

AEM systems use towed-loop sources to transmit electromagnetic (EM) fields into the subsurface. On the ground, loop sources and electric bipole sources can be realized. Because contributions from opposed wire segments cancel out at large distances outside of a loop source, fields from this source decay more rapidly in space when compared to a long bipole source. This suggests that bipole sources are more useful for semi-AEM surveys, at least in environments in which ground coupling of electrodes can be established. Furthermore, an electric bipole is a grounded galvanic source and excites vertical and horizontal (or poloidal and toroidal) modes of current flow in the subsurface (as opposed to the primary currents induced by an inductively coupled loop source, which are entirely horizontal). To take advantage of the complementary information contained in the two modes, our concept includes the installation of additional (sparse) ground-based electric field recorders (see Figure 1) to measure the bimodal electric field. This can be achieved at the expense of only a little extra logistic effort.

Within the deep electromagnetic sounding for mineral exploration (DESMEX) project, we aim at establishing semi-AEM for the exploration of deep-seated resources of raw materials. Typically, shallow deposits in Germany, as well as in other densely populated areas of Europe, are exploited, but the potential at greater depths is unknown. The challenge to explore in these environments is related to a high level of background noise associated with a high population density, which makes it difficult to perform natural source audiomagnetotelluric surveys. However, the survey areas are typically accessible on the ground, which makes it logistically feasible to install transmitters (and an auxiliary set of stations) on the ground. Therefore, the concept of semi-AEM is suitable for such areas.

Two airborne receiver systems have been developed in the project. In this paper, we describe in detail the setup of and test flight results obtained with an airborne receiver, which combines fluxgate and purpose-built broadband induction coil magnetometers with a high-precision inertial navigation system (INS) on a rigid platform housed in a so-called bird. Although data from this system have been already used to report the feasibility of 3D semi-AEM inversion (Smirnova et al., 2019), the objective of this work is to present the technical specifications and the characteristics of the measured data in detail. The second airborne receiver developed within the project relies on a novel generation of superconducting quantum interference device-based magnetometers. This system and the results will be reported separately because both systems require technical descriptions and data processing workflows.

The test measurements were conducted in September 2016 near the town of Schleiz in southern Thuringia, Germany (Figure 2; Smirnova et al., 2019). The region hosts several antimony deposits, which were mined until the late 1950s. The abandoned mine “Grube Halber Mond” in the east of the village of Oberböhmsdorf is located centrally within the flight area. Although antimony is electrically not conductive and cannot be targeted directly with EM surveys, the survey area covers a highly heterogeneous and tectonically faulted and folded geology with occurrences of electrically conductive Silurian alum shales embedded within moderately resistive to resistive Ordovician silt shales, Devonian diabase intrusions, and Carboniferous graywackes. Preliminary ground geophysical surveys and petrophysical studies have revealed electrical resistivities of less than 10  Ωm for the alum shales and background resistivities of 300–5000 Ωm for the other sequences, thus providing good parameter contrasts for a test survey area.

Three transmitter locations were deployed during the survey (the blue, red, and green lines in Figure 2) and overflown in separate flights but along coincident flight lines (the black lines in Figure 2). The flight velocity was approximately 33 m/s, and the altitude of the bird was 60 m above ground on average. The helicopter path was following the topography, which is in the range of 400600 m above sea level in the area. The flight lines for each transmitter location were flown in a total time of 2.5 h, covering an area of 4 × 5 km2 each, with 21 parallel lines in the northwest direction and four perpendicular control lines. Flight lines perpendicular to the sources TxB3 and TxB4 have a spacing of 100 m in the central part of the survey area and have been increased to 300 m to the northeast and southwest. The spacing of the four control lines is 1 km. In addition, a ground station has been installed in the vicinity of the three transmitters (the magenta star in Figure 2).

We use the data from the test survey described above to illustrate our processing steps. We first introduce the principle objective of processing semi-AEM data in the frequency domain. Then, we describe the details of the technical components integrated in the airborne platform, followed by a detailed outline of our data processing scheme. Data processing involves two main steps: (1) removing motional noise from the raw recordings and (2) estimating transfer functions in the frequency domain from spectral estimates of the magnetic flux density and the current. The processing results are validated by comparison with ground recordings and with comparisons of half-space responses. A more comprehensive 3D inversion of the data has been presented in a separate publication (Smirnova et al., 2019).

The primary objective of semi-AEM data acquisition and processing in the frequency domain is the estimation of the spatial distribution and the frequency dependency of the magnetic transfer functions Bx,y,z, which relate the observed magnetic flux densities Bx,y,z(ω,r) to the injection current I(ω) as
Bx,y,z(ω,r)=Bx,y,z(ω,r)I(ω).
(1)
Here, ω is the angular frequency and r are the observation points. An equivalent expression holds for electric transfer functions Ex,y, that is,
Ex,y(ω,r)=Ex,y(ω,r)I(ω),
(2)
where Ex,y(ω,r) represents the ground-based electric field measurements. All quantities in equations 1 and 2 are complex-valued. The terms Bx,y,z and Ex,y are dependent on the particular measurement geometry and the subsurface conductivity structure but independent of the particular current waveform (Streich et al., 2011). Quasi-analytical expressions for Ex,y and Bx,y,z can be found for a layered half-space with a conductivity depth profile σ=σ(z) and for dipole sources, straight wire sources of finite length (Streich and Becken, 2011), as well as magnetic dipole sources. More complicated source geometries can be approximated by a superposition of elementary sources. In the presence of 2D and 3D conductivity structures, or in the presence of topography, numerical solutions to Maxwell’s equations must be used (Streich, 2009; Rochlitz et al., 2019).

For this method, current waveforms with sharp switch-on and switch-off ramps are suitable for measurements of Bx,y,z (and Ex,y) over a broad frequency range because a broad spectrum of harmonics will be excited. It is mandatory that measurements of the injection current and of airborne data are precisely time synchronized. We use GPS time as the time reference and the same type of data loggers for the current recordings and the magnetic field recordings to exclude any hardware-specific time delays.

The estimation of Bx,y,z is challenging in airborne applications because the magnetometers are constantly translating and rotating at high velocities and angular rates, respectively. Motional noise compensation is therefore a crucial task because these disturbing signals overlap the signal frequency range under investigation. Because the motion is not stationary, this correction using precise INS data must be applied in the time domain prior to spectral estimation. The INS attitude data are therefore also time-synchronized via GPS. The general workflow for the estimation of transfer functions is schematically laid out in Figure 3; the indicated processing steps are detailed in the following sections.

Ground transmitter

We have tested two different transmitters. For the measurements reported in this paper, the current was supplied by a Zonge GGT30 transmitter with a 30 kVA power generator. The maximum output voltage of the system is 1000  V, and we injected up to 24  A into approximately 1  km long grounded bipoles. Electrodes were built of arrays of metal rods and plates. The Zonge GGT30 transmitter is especially designed for time-domain applications in a long-offset transmitter-receiver configuration (Hördt et al., 1992; Mörbe et al., 2020) and can generate current waveforms with a short turn-off time. A second transmitter, which has been successfully used, has originally been designed for direct-current (DC) resistivity surveys. It uses a 40 kVA power generator, with a maximum output voltage of 1500 V and a maximum output current of 50 A. Both transmitters are programmable in repetition rate as a multiple of 1  ms and in duty cycle. For the data discussed here, we injected a rectangular waveform with 50% duty cycle and a repetition rate of 96  ms (f0=10.416¯  Hz). The output waveform of the current is measured with current clamps during the transmission time and stored as GPS-synchronized time series at a sampling rate of 16,384  Hz. The current recording unit is an ADU-07e 24 bit data logger from Metronix Geophysics. The transfer function of the current clamp is flat and displays a sensitivity of ci=10  mV/A output voltage; the ADU-07e has, however, an analog Butterworth high-pass filter of order one with a cutoff frequency of 1  Hz integrated at the signal input, which can be expressed by the filter impulse response chp(t). Hence, the recorded voltage is
ui(t)=cichp(t)i(t),
(3)
where i(t) is the current and denotes the convolution operator. The voltage ui(t) has the spectrum
Ui(ω)=F{ui(t)}=ciChp(ω)I(ω),
(4)
where F{.} denotes the Fourier transform and the convolution theorem of the Fourier transform has been used.

An example recording of the source current — injected through transmitter TxB4 (see Figure 2) — with a strength of ±20  A is depicted in Figure 4. We display a 1  s time-series window with approximately 10 switching cycles (Figure 4) and the power spectral density (PSD) of the current (Figure 4), obtained from averaging 13 windows of a short-time Fourier transform (STFT) with a window length of 200 cycles approximately 20  s each. For both displays in the time and frequency domains, voltages have been converted to amperes, but the high-pass filter has not been deconvolved. The time-series waveform is therefore slightly distorted. The PSD displays the fundamental frequency f0 and approximately 400 odd harmonics nf0,n=1,3, in the frequency range from 10  to  8000  Hz. These coefficients will be evaluated in the processing scheme devised later.

Sensor platform

The sensor platform is a 10 m long glass-fiber reinforced plastic tube. The outer shell is based on an out-of-service DIGHEM V bird, which was refurbished. The bird has two suspension points with two ropes of slightly different length attached to the main rope of 40 m length. The front rope is shorter to maintain a horizontal orientation of the system during flight. A skirt is mounted at the back of the bird to ensure the in-flight stability of the towed platform. The instruments are attached to four wooden sleds within the tube. A system design that suppresses vibrations is desirable to reduce motion induced noise in the magnetometer systems. We rely on heavy batteries and a low center of mass of the bird to increase the in-flight stability of the system; we have abstained from mechanical suspension systems within the tube because they might impact the overall mechanical stability of the bird.

The bird contains two three-component (3C) magnetometers to detect the magnetic field variations: MFS-11e induction coil magnetometers from Metronix Geophysics and a MAG-03 Fluxgate from Bartington Instruments. The motion and attitude of the bird are recorded with the iNAT-RQT-4003 from iMAR Navigation. The magnetometers are connected to an ADU-07e, and the INS system includes its own separate data logging unit. In addition, an optically pumped GSMP-35U potassium magnetometer from GEM Systems is integrated in the bird. Three GPS antennas are mounted on top of the bird, one antenna for each of the recording systems. The ADU-07e and the INS can be accessed through an Ethernet connector panel; an optional WiFi connection is available to examine the magnetometer recordings and the bird motion during the flight from within the helicopter. A laser altimeter is attached at the bottom of the tube. For the test flights described here, the potassium magnetometer and the laser altimeter have not yet been operational and are therefore not discussed.

The arrangement of the different systems and antennas within and on top of the bird is illustrated in Figure 5. The INS has to be rigidly attached to the magnetometer platform to perform the same motion as the magnetic field sensors. However, the INS is an active measurement system and as such emits electromagnetic noise, which can disturb the magnetic field measurements over a broad frequency range. To evaluate the optimum distance of the sensors within the bird and to avoid mutual interferences of the systems, we performed several distance tests on the ground. We found 60 cm to be a reasonable distance between the INS and the three induction coils. The fluxgate is mounted approximately 2  m away from the INS and 1.4 m away from the coils. More details can be found in Nittinger (2018).

The induction coil sensors are thoroughly designed to be sensitive to the desired frequency range of 1–10 kHz. The MFS-11e induction coil is a unique sensor system manufactured by Metronix solely for this airborne system. The coil has an integrated feedback system resulting in a calibration function with a high-pass cutoff frequency of 32  Hz, a dynamic range of ±78  nT, a noise level of 30  fT/Hz at 10  Hz, and a sensitivity of 135  mV/nT in the frequency range, where the coil response is flat (see Figure 6). It is adapted to the conditions in the moving platform by suppressing the motion-induced voltages from the dominant long-period pendulum motion of the bird (in the range of a few seconds).

With a length of 71.5 cm and a diameter of 7.5 cm, the coils are too long to fit into the bird shell in an orthogonal arrangement. Instead of adding holes to the outer shell, which might impair flight stability, we used an oblique arrangement of the coils (the insets in Figure 5). Test measurements on the ground have not retrieved any noticeable interference between the coils in this oblique geometry (Nittinger, 2018). For subsequent processing, the coil readings must be transformed into the orthogonal body coordinate system with the unit vectors x^b,y^b, and z^b. We refer to the three coils as u, v, and w, arranged in the u^,v^, and w^directions. Then, these oblique sensor axes can be expressed in the body coordinate system as
[u^v^w^]=TC[x^by^bz^b],
(5)
where TC is an oblique 3×3 transformation matrix that is constructed from the orientation angles of the coils out of the body coordinate system as depicted in the insets in Figure 5. In detail, coils u,v,andw are rotated by αu=16°,αv=16°, and αw=30° (positive clockwise) out of the long axis of the bird (the x^b-axis) and dipping by βu=26.57°, βv=26.57°, and βw=0° (positive downward). Then,
TC=[cosαucosβusinαucosβusinβucosαvcosβvsinαvcosβvsinβvcosαwcosβwsinαwcosβwsinβw].
(6)
Thus, the voltage outputs of the coils, uu,uv, and uw, can be converted into output voltages ux,b,uy,b, and uz,b, respectively, of a hypothetical orthogonal arrangement of coils aligned with the sensor platform as
(ux,b,uy,b,uz,b)T=TC1(uu,uv,uw)T.
(7)

The Bartington MAG-03 provides a bandwidth of 3  kHz with a flat response up to 1 kHz (see Figure 6) and a noise level of 4–6 pT/Hz at 1  Hz. The sensitivity is 0.1  mV/nT. The three-axis sensors possess an orthogonality error <0.1°. The frequency-domain response function of the fluxgate was validated in the Metronix calibration facility. Calibration parameters regarding zero-pole offsets, scale factors, and orthogonality of the sensors were determined using the method of Auster et al. (2002).

All magnetometers are sampled at a rate of 16,384 Hz. The fluxgate magnetometer is oversampled, but it is convenient to have a common sampling rate for all field components and for the current. As for the current recordings, an analog Butterworth high-pass filter of order one with a cutoff frequency of 1 Hz is applied to signals from all magnetic field sensors prior to digitization. The filter is effective in suppressing long-period motional noise and removes the DC part of the field: a requirement to set a gain to the fluxgate output without exceeding the input range of the data logger. Here, we use an amplification of 64 to take advantage of the full resolution of the fluxgate sensor; no gain is set for the induction coils.

The maximum sampling rate of the INS is 400  Hz. Attitude data are internally processed, and outputs include attitude (the accuracy is 0.01° according to the technical specification provided by the manufacturer), angular velocity, acceleration, position, and velocity.

During the flight tests presented next, the system is towed by the research helicopter Sikorsky S-76B, operated by the German Federal Institute for Geosciences and Natural Resources (BGR). The helicopter is also equipped with an INS and a laser altimeter to record the position of the helicopter and the height above ground. It also provides a WiFi router and corresponding antennas to access and operate the system during flight and at the airport.

Removal of motional noise

Let us denote the magnetic flux density at some point in space r and time t as b(r,t)=b0(r)+bsig(r,t)+μ(t), where b0(r) is the static geomagnetic field, bsig(r,t) is the signal due to current injection i(t) at the transmitter, and μ(t) is ambient noise superimposed on the data. Measurements of b(r,t) are performed in flight mode with magnetometers attached to the towed platform. The attitude of the platform is changing with time and expressed in terms of special Euler angles ϕ(t),ϑ(t),andψ(t), referred to as roll-pitch-yaw angles as is common in aircraft navigation. In the body axis of the sensor platform, the magnetic flux density can thus be expressed as
bb(r,t)=RT(t)b(r,t),
(8)
where
R(ϕ,ϑ,ψ)=R(ψ)R(ϑ)R(ϕ)=[cosψsinψ0sinψcosψ0001][cosϑ0sinϑ010sinϑ0cosϑ][1000cosϕsinϕ0sinϕcosϕ]
(9)
is the time-dependent rotation matrix (time dependency has been omitted) and subscript b denotes a vector in the body coordinate system.
The dominant contribution in equation 8 comes from the geomagnetic field b0(r) because the sensor rotation translates into temporal flux-changes in the body frame. This gives rise to motional noise as well as to motion-induced noise, the latter term referring to the fact that magnetometers with a nonflat frequency response are also sensitive to temporal flux changes. In our case, we apply a high-pass filter to the voltage output of the magnetometers prior to digitization. Combining the sensor impulse responses cis(t) with the impulse response of the analog high-pass filter of the recording unit chp(t), the registered voltage ub,i for the ith component of a rotating magnetic 3C sensor can therefore be expressed by
ub,i(r,t)=chp(t)cis(t)bb,i(r,t).
(10)
With knowledge of RT(t), b0(r), and the response functions of the sensor components and the high-pass filter, we can predict the motion effects in the geomagnetic field from
ub,ipred(r,t)=chp(t)cis(t)[RT(t)b0(r)]i,
(11)
and subtract this quantity from the measured voltages ub,iobs(r,t),
ub,isig(r,t)=ub,iobs(r,t)ub,ipred(r,t),=chp(t)cis(t)[RT(t)bsig(r,t)]i,
(12)
where the ambient noise term has been omitted and bsig is the magnetic flux signal of interest. The residual voltage as defined in equation 12 is thus considered as the signal.
Recovery of the EM signal bsig from equation 12 involves deconvolution of the sensor and filter response functions followed by reverse rotation. Because deconvolution is a numerically unstable operation, we would like to avoid the explicit removal of the filter functions in the time domain. Direct transformation of equation 12 into the frequency domain would allow us to calibrate the voltage spectrum by division with the spectral representations of the filter functions; unfortunately, however, reverse rotation would correspond to a deconvolution in the frequency domain. To avoid these complications, we tested as an approximate solution to equation 12 a change in the order of multiplication (rotation) and convolution (filter processes) and express the 3C vector of residual signal voltages in body coordinates as
ubsig(r,t)RT(t)[chp(t)cs(t)bsig(r,t)].
(13)
Left multiplication with R(t) then yields the voltages in an earth-centered earth-fixed frame, that is,
usig(r,t)R(t)ubsig(r,t)chp(t)cs(t)bsig(r,t).
(14)
Equation 14 now has the advantage that the filter responses can easily be corrected for by making use of the convolution theorem of the Fourier transform after transformation of the signal voltage into the spectral domain.

The approximation described above is numerically stable but has several shortcomings. First, note that the individual sensor response functions cis(t) had to be substituted by a common response function cs(t), which is possible only under the assumption that the three sensors’ components exhibit identical frequency characteristics. This assumption is very well satisfied for the sensors used.

Second, the rearrangement in equation 13 is only applicable if either (1) the filter functions have a flat frequency response in the frequency range of interest or (2) if the frequencies of motion do not overlap with the signal frequencies. The former condition is approximately satisfied for the high-pass filter and the fluxgate sensor, but it is violated for the induction coils at frequencies below 100  Hz (Figure 6). With respect to the latter condition, we observe dominant pendulum motions of the bird excited in the period range between 1  and20  s, but higher frequency vibrations are also encountered. Removal of motion-induced voltages in the induction coils due to such vibrations is therefore only possible to within the approximation made above. Later, we show that transfer functions derived from the fluxgate and induction coil data compare very well at low frequencies (see Figure 12). This suggests that the errors introduced above are not significant and that the approximation made in equation 13 can be justified in practice.

In Figure 7, we display the raw fluxgate and induction coil time series (the blue line) recorded along a segment of flight line L22. The shown time window corresponds roughly to a flight distance of 2  km; the transmitter has been traversed at a relative time of 50  s. The data in Figure 7 correspond to sensor output voltages, converted to nT with the sensor sensitivities from the calibration procedure. Superimposed on the raw data are the predicted voltages (equation 11, the red lines in Figure 7), based on the attitude data from the INS describing the motion of the sensors in the geomagnetic field. Here, we determined the local geomagnetic field vector from the international geomagnetic reference field (IGRF) (Finlay et al., 2019). Because the area covered by a single flight is only on the order of 20km2, the effect of spatial variations of the reference field (which are in the range of a few nT) on the compensation procedure for EM data is small and can be disregarded. It also can be noted that the effect of geomagnetic field gradients on AEM responses is generally neglectable on these scales, even in cases of strong field gradients due to local magnetic anomalies (Smith, 1999). The IGRF field does not take into account daily variations and local magnetic anomalies; these deviations from the local conditions are therefore also disregarded here.

We observe field variations on the order of ±1000  nT (±100  mV sensor output) for the fluxgate readings and ±10  nT (±1350  mV sensor output) for the induction coils, which can be almost entirely explained by motional noise. The EM signal is superimposed on the motional noise and is clearly visible in the induction coil data when the helicopter is traversing the transmitter; because the sensitivity of the fluxgate data to flux changes is much lower in the excited frequency range, the EM signal is less obvious in the raw fluxgate data.

The prediction of voltages due to motion is limited to 200  Hz, that is, half the sampling rate of the INS. For further processing, we interpolate all data to a common sampling rate of fs=20,000  Hz using spline interpolation. Interpolated magnetic field predictions using the INS data are low-pass filtered with a corner frequency of 200  Hz to ensure that high-frequency contents, artificially introduced by the interpolation scheme, are suppressed. Before subtracting the interpolated time series one from another, we test the time synchronization by cross-correlating both series. Synchronization is not guaranteed because different hardware and software components used in the ADU07e and the INS can lead to system-specific time delays. In our case, we find a lag of 8.5  ms, corresponding to the distance of 3.5 samples at the original INS sampling rate, and we shift the INS data accordingly.

Despite the high precision of the INS recordings, attitude data are afflicted with noise, which propagates further into the predicted sensor voltages via equation 14 and may thus impact the further processing. Simple estimates can show that attitude uncertainties of 0.01° (i.e., the accuracy claimed for the INS deployed in the bird) can be the cause of signal well above the noise level of the magnetic field sensors. Further statistical analysis (i.e., based on coherence measures of observed and predicted voltages) could help to mitigate noise propagation from the INS into the EM signal, but this has not been investigated yet.

Heading errors and misalignments

Figure 8 depicts the magnetic field residuals after subtracting the predicted signal as determined with equation 11. It can be seen that the fluxgate residuals (the red lines in the left panels of Figure 8) are still afflicted with motional noise in the order of 50–100 nT; residuals for the induction coil data are significantly less but still obvious (the red lines in the left panels of Figure 8). These residuals can be further reduced when taking into account undetected misalignments of the sensor components (Auster et al., 2002) as well as so-called heading errors (Leliak, 1961).

Misalignments and offsets are calibration parameters of vector magnetometers and are usually determined from calibration measurements in multiple orientations (Auster et al., 2002). Heading errors describe distortions of the magnetic field at the sensor position due to the influence of nearby permanent magnets, and of statically and dynamically induced magnetic fields in magnetically permeable and in electrically conductive parts of the platform, respectively. These sources of errors can be compensated for, because they can be predicted from the orientation of the aircraft and its rate of change of orientation using calibration data collected in an area of low magnetic relief (Leliak, 1961).

Offsets that are constant in the coordinate frame aligned with the sensor axes can be ignored in our setup because they are removed by the hardware high-pass filter as explained in the previous section. It means that intrinsic sensor offsets as well as the influence of permanent magnets do not play a role for the registered signal. However, due to sensor misalignments, the actual signal readings correspond in every component to a weighted linear combination of the 3Cs of the magnetic field that would be measured with an ideal 3C magnetometer triple. Hence, a total of nine weighting coefficients is required to describe sensor misalignments.

Statically induced magnetic fields in magnetizable parts of the platform are proportional to the inducing magnetic field. Thus, they also can be expressed by a linear combination of the inducing field weighted with constant coefficients, similar to the description of sensor misalignments.

If we take the inducing field to be the magnetic flux in the body frame, bb(r,t)=RT(t)b0(r), express the heading error as bhead(r,t)=Cheadbb(r,t) and take into account sensor misalignments by projecting the field components onto the sensor axes via the projection matrix Cmis, then the predicted field can be expanded to
bbpred(r,t)=Cmis(bb(r,t)+Cheadbb(r,t))=Cbb(r,t),
(15)
where we introduce C=Cmis(I+Chead) and I is the identity matrix. It means that sensor misalignments and heading errors arising from statically induced fields in parts of the platform are indistinguishable and can be captured by a single 3×3 coefficient matrix C.

Eddy currents in electrically conductive parts of the platform can add to further heading errors. Because eddy currents are driven by the rate of change of the inducing field rather than the momentary field, obeying the law of induction, a separate term in equation 15 would be needed that is proportional to the time derivative of the inducing field. We omit this term here for the sake of simplicity, assuming that its contribution is minor compared to sensor misalignments and statically induced fields.

It means that in place of equation 11, the predicted voltage in each sensor component i{x,y,z} can be expanded to
ub,ipred(r,t)=chp(t)cis(t)[CRT(t)b0(r)]i.
(16)
Under the premise that the sensors have identical impulse responses, we can rewrite equation 16 as
ubpred(r,t)=C[chp(t)cs(t)[RT(t)b0(r)]],
(17)
where it is implied that each component of the vector of [RT(t)b0(r)] is convolved with the filter responses. The coefficients of C are estimated row-wise from flight data in a least-squares sense using linear regression techniques to minimize the squared residuals
ub,iobs(r,t)ub,ipred(r,t)22min.
(18)
For the induction coil data, the minimization can be either performed in the orthogonal body coordinate system of the bird, or in the oblique sensor axes of the coils. A corresponding mapping is achieved by making use of equation 7.

This minimization problem is ideally performed on data from a calibration flight, for example, at high altitudes in a homogeneous geomagnetic field and with flights in multiple directions to cover a broad range of sensor orientations. However, for the sake of fitting the predicted voltage to the observed voltages, we found that the coefficients in C can be reliably estimated from single flight line data as well.

Having estimated C, the compensated signal voltages can be recovered in the body frame as C1ubsig(r,t), taking the residuals as the signal of interest (see equation 12). Finally, compensated signal voltages in an earth-fixed frame are given by
usig(r,t)R(t)C1ubsig(r,t),
(19)
and are used for further processing instead of the uncompensated voltages from equation 14.

The residuals, obtained with the compensation procedure described above, are plotted in Figure 8 on top of the residuals without compensation. Compensated residuals as given by equation 19 are depicted in blue, and uncompensated residuals (equation 12) are depicted in red. The compensation can significantly reduce the remaining motional noise in the fluxgate data from approximately ±100 to ±10  nT (Figure 88), but appears to be of minor importance for the induction coil data (Figure 88). Figure 88 and 88 displays a magnification into the time series of the compensated fluxgate and induction residuals, respectively, to emphasize the EM signal contained in the time series.

Estimation scheme

In this section, we describe our scheme to estimate magnetic transfer functions Bi(ω,r) between the spectral components of any of the measured magnetic flux density components Bi(ω,r) and the injection current I(ω). We assume now that all signal components have been rotated and projected in a geostationary reference frame. Let us denote the Fourier transforms of the high-pass-filtered signals as
Y(r,ω)=Chp(ω)Bi(ω,r)=Cs,1(ω)F{uisig(r,t)}
(20)
and
X(ω)=Chp(ω)I(ω)=F{uii(t)},
(21)
respectively. The linear relation given by equation 1 is then equivalent to
Y(r,ω)=Bi(r,ω)X(ω),
(22)
and X(ω) and Y(r,ω) can be conceived as the input and output of a linear system. Equation 22 corresponds to an univariate regression problem, which can be solved for Bi(r,ω) in a straightforward manner (Streich et al., 2011; Smirnova et al., 2019).
It is implied in equations 2022 that the measured signals can be attributed to a single position r=rwin. This is only approximately true, however, because spectral analysis must be carried out over temporal windows of finite duration, in which the receiver is towed through the air; that is, r=r(t). As a consequence, the signal is not stationary (see Figure 8), particularly when the receiver approaches or moves away from the transmitter. Instead, the field can be thought of as an amplitude-modulated signal
y(r(t),t)env[y(r(t))]y(rwin,t),
(23)
where env[.] denotes the envelope over the considered time window and y(rwin,t) corresponds to the signal at a point in the center rwin of the flight path covered during the considered time window. Note that the spectrum of the signal as given by equation 23 contains the spectral lines associated with the signal frequencies as well as satellite frequencies due to the amplitude modulation with the envelope. At distances from the transmitter in which the envelope of the signal varies slowly over a given temporal window, reliable spectral estimates of the EM signal can be obtained and attributed to a central point rwin. However, it is obvious that, in the immediate vicinity of the transmitter, the envelope of the signal exhibits a strong curvature, and the spectrum can be expected to be distorted in this region.

In practice, spectral estimates are obtained from the time series by means of an STFT. The time window length Lwin is chosen in exact multiples Ncyc of the fundamental cycle duration, Lwin=Ncycfs/f0; this ensures that the satellite frequencies due to the amplitude modulation do not affect the signal frequencies. We choose the temporal windows to overlap by 50% and apply a Hanning window as a taper. In our example, f0=10.416¯  Hz=1/0.096  s, the interpolated sampling rate is fs=20,000  Hz and one cycle corresponds exactly to Lwin=1920 samples.

In the spectral domain, the EM signal is contained in the fundamental frequency f0 and a large number of odd harmonics nf0,n=1,3, (see Figure 4). As is common in electromagnetics, the response functions vary smoothly with frequency. Hence, we bin all signal coefficients of the input and output channels into bands of one-half octave bandwidth 2log2fc±1/4 centered at logarithmically equidistant estimation frequencies fc, and collect them for each band and for a few adjacent temporal windows Nw in data vectors X and Y. The components of X and Y are considered as independent realizations of equation 22 and admit a least-squares estimate of the transfer function
Biest(rc,fc)=[XTX]1XTY.
(24)
Direct solution of the regression problem in equation 22 ignores, however, the spatial variation of Bi(rwin,ω) over adjacent windows. To improve our estimate, we approximate Bi(rwin,ω) to be locally linear around the central position rc as
Bi(rwin,ω)Bi(rc,ω)+dBi(r,ω)dr|r=rc(rwinrc),
(25)
and we solve this equation for Bi(rc,ω) and the slope in the flight direction dBi(rc,ω)/dr in place of equation 22. The augmented system of equations corresponds to fitting a straight line through the observations, where Bi(rc,ω) is the intercept and dBi(rc,ω)/dr is the slope (Menke, 2012). The primary effect of the expansion 25 is a reduction of the prediction error (and thus the estimation error) in regions where the slope is large (e.g., near the transmitter), whereas the estimate Biest(rc,ω) is nearly unaffected.

We use standard and robust least-squares regression techniques (Huber, 1981) to solve for the transfer function and the slope. The regression result is then attributed to the central position rc and the central frequency fc. The two lowermost frequency bands bin only one and three signal coefficients, respectively. In these bands, the transfer function is attributed to the fundamental frequency and the logarithmically averaged frequencies of the first three odd harmonics, respectively, rather than the logarithmic center of the frequency band. Given that only a few windows can be averaged, robust regression schemes are not effective at low frequencies, but they might also be less required because the signal is strongest at the lowest frequencies. However, as the number of harmonics within each frequency band increases with frequency, robust schemes may help to improve the estimates from equation 25 at high frequencies.

The window length Ncyc and the number of averaging windows Nwin are important parameters. Short temporal windows yield a good spatial resolution but a low spectral resolution and vice versa. We tested the effect of different window lengths and window numbers on the regression result. Figure 9 depicts the real part of Bz(rc,fc) as a function of the distance from the transmitter and for selected frequencies. Transfer functions are shown for individual Fourier coefficients and windows, that is, the ratio Y(rwin,nf0)/X(nf0) using small symbols, and the large symbols are regression estimates using all Fourier coefficients within each band and from adjacent windows, the latter chosen to overlap by 50%. We show results for window lengths of 8 and 32 cycles and for 2 and 8 windows for each of these cases.

It is seen from all panels in Figure 9 that the real part decays strongly away from the source and with frequency, as expected. The decay is asymmetric, suggesting higher conductivities toward the northwest (negative distances) than to the southeast (positive distances). In Figure 9, we show results for Ncyc=8 and Nwin=2, corresponding to a spatial averaging length of approximately 45  m, which gives a high spatial resolution but yields some scatter of responses with offsets increasing. The robustness of responses is improved with a larger number of windows (Figure 9, Nwin=8, averaging length of approximately 135  m) or with increasing window length (Figure 9, Ncyc=32, Nwin=2, averaging lengths approximately 185  m; Figure 9, Ncyc=32, Nwin=8, averaging lengths approximately 550  m), but at the expense of loss of resolution, particularly near the transmitter. We therefore conceive that the choice of these parameters should be adopted to the frequency and the distance from the transmitter, that is, the signal strength and use Ncyc=8 and Nwin=8 for high spatial resolution and Ncyc=32 and Nwin=8 for improved spectral resolution and robustness.

The vertical magnetic flux density decays rapidly away from the source. However, the horizontal components dominate at large distances, where the field becomes more homogeneous. In Figure 10, we depict the real and imaginary parts for the same line and frequencies as in Figure 9. Note that stable transfer functions can be estimated at high frequencies and large offsets for the horizontal component. The imaginary part of the transfer function is purely inductive, whereas the primary field is entirely real. The primary field dominates the response in the near-field of the transmitter; however, the imaginary part becomes increasingly important with increasing distance and at high frequencies, as expected for inductive responses.

In Figure 11, we display maps of transfer functions, derived for the fundamental frequencies, 10.416¯  Hz (Figure 11 and 11), and for 1448.2  Hz (Figure 11 and 11). We depict the imaginary parts of the vertical and northward components I{Bz} and I{Bx}, respectively. All imaginary components reflect a spatially consistent scattered field, which is in support of good data quality. Noisy estimates are seen in regions where the imaginary part of one component of a transfer function is close to zero, for example, in the eastern corner of the flight map (Figure 11, Bx, f=10.4  Hz). Note also from the data maps that — in accordance with potential field characteristics — the horizontal components are large in regions where the vertical component changes sign and vice versa. This suggests that it is useful to make use of the vertical and horizontal components. We also observe from the maps in Figure 11 that visually indistinguishable values are obtained at crossover points on line and tie-line recordings, suggesting correct timing and no strong altitude dependency of the data.

Finally, to validate our semi-AEM response estimates, we compare fluxgate and induction coil data with the response from a ground station (station A) under operation during the survey, which was overflown by the helicopter-based instrument. The location of the ground station is indicated in Figure 2. The estimation of the response at the ground station has used 20 min of data, as opposed to the approximately 5  s of flight data used above. Frequency-sounding curves, shown in Figure 12, depict the responses of the ground station (the dashed line) and the responses estimated from airborne induction coil data (the open circles) and separate estimates for airborne fluxgate data (crosses). Fluxgate responses are limited to frequencies below 150  Hz because of the low sensitivity of the fluxgate at higher frequencies. We first note that the airborne data are consistent with the ground data. The average differences between the real parts of the horizontal components and the imaginary parts of all three components of the airborne and ground data are below 10  pT/A; the real parts of the vertical component Bz of the airborne estimates are however systematically larger by 1020  pT/A than the estimates from the ground recording. The exact reason for this observation is not entirely clear; we conjecture that this effect is due to the spatial averaging involved in the estimation of airborne transfer functions.

We further note from Figure 12 that the fluxgate and airborne data are consistent at low frequencies. Any systematic discrepancy between data from these sensors would point toward effects that were not correctly accounted for during the processing. We tested the relation between transfer functions obtained from induction coils and the fluxgate statistically. First, we collect the six low-frequency transfer function estimates of all three components for the validation point (site A) discussed above in a crossplot (Figure 13) and perform a linear regression. The regression yields an intercept of 6·104  nT/A and a slope of 0.998. This indicates the equality of induction coil and fluxgate estimates. Second, we construct the same crossplot, but for all estimation points within a distance of less than 1500  m from the transmitter. In total, 619 points, three components, and six frequencies are combined for this analysis. The regression on these data yields a vanishing intercept (2.9·105±3·105  nT/A) and a slope close to identity (1.01±4·104, see Figure 13), where the uncertainties have been estimated from the standard deviation of data residuals propagated into model covariances. The slope estimate is biased by approximately 1%, probably reflecting the degree of inaccuracy of the sensor calibrations. Third, the distance ΔBi between the induction coil and fluxgate transfer functions is binned to a histogram (Figure 13). We observe a balanced distribution centered at zero and with a standard deviation σΔ=6  pT/A for the real part and σΔ=4  pT/A for the imaginary part. These numbers thus indicate the noise estimates of the combined distribution of induction coil and fluxgate data.

All of these results do not indicate any systematic bias or error in our analysis and thus support the technical feasibility of our data acquisition and processing approach. In particular, the coincidence of fluxgate and induction coil responses at low frequencies suggests that the approximations made for handling motional noise are justified (see equation 13).

Comparison to half-space responses

Ultimately, the subsurface conductivity distribution must be extracted from the semi-AEM transfer functions described above by means of inversion. As a first proxy for the information content of the data, EM data and geoelectrical data are often transformed to apparent resistivities (or conductivities), for example, by performing half-space inversions. In the presented case, we found that such data display are not useful because of the 3D nature of the conductivity structure and the large inductive volume that contributes to local responses. For the same reason, we found that 1D inversions are not adequate although they can fit portions of the data locally. Therefore, full 3D inversion is generally necessary and has also been accomplished by Smirnova et al. (2019) for subsets of the data discussed.

Here, we compare estimated responses to half-space responses in terms of transfer functions rather than attempting to transform the data. As an example, the vertical transfer function Bz is shown in Figure 14 for the frequency f=362  Hz along two profile lines, L22 and L26. The data are shown for two independent transmitters, TxB4 and TxB3, which are located at 0  km and at 1  km profile distances, respectively. Profile line L22 crosses both transmitters, whereas profile line L26 passes the transmitter at a distance of approximately 500  m. We only show the vertical component, but we note that similar data quality with redundant information is available for the horizontal components.

The transfer function estimates are shown for two sets of STFT window lengths, using 8 (the gray symbols) and 32 (the colored symbols) signal cycles, respectively. It should be noted that these sets of estimates coincide except for regions where the field exhibits a strong gradient along the profile. Therefore, the effect of spatial averaging using long windows can be clearly observed in terms of smoothed spatial variations associated with increased error estimates.

We compare the estimated responses with the responses simulated for a half-space with resistivity of 10, 100, and 1000  Ωm, respectively (the dotted, solid, and dashed-dotted lines in Figure 14). In the vicinity of the transmitter, the real part is dominated by the primary field. Inductive effects are evident in the real part at some offset to the transmitter and throughout in the imaginary part. The data are largely consistent with the response of a 100  Ωm half-space, but deviations from this simple model are also obvious and indicate lateral and/or vertical resistivity variations. This is consistent with the 3D model presented by Smirnova et al. (2019), which contains several conductivity anomalies embedded within a background of a few hundred Ωm. In conclusion, we state that the transfer function estimates appear plausible in view of the simulated primary and the scattered field components for a homogeneous half-space.

In our novel semi-AEM system, we integrated induction coil magnetometers (in an oblique arrangement) and a fluxgate sensor in an airborne sensor platform to measure the three components of the EM field generated by a transmitter located on the ground. The induction coils are modified versions of induction coils that are normally used in ground-based magnetotelluric or controlled source EM surveys. Motional noise is a major source of noise in AEM surveys. The motion and attitude of the platform are therefore monitored with high-resolution INS.

The induction coils that were deployed in the new system represent a balance between sensor sensitivity and signal resolution. An important aspect in the sensor design is the integrated feedback system of the induction coils. The feedback system flattens the coil sensitivity within a frequency range of approximately 100–5000 Hz. In this range, the output voltage is proportional to the magnetic flux rather than to temporal flux changes. It means that the deployed induction coils are much less sensitive to high-frequency motion in the geomagnetic field as compared to air coils.

With the present design, sensor output voltages are within a few volts and within the dynamic input range of the 24-bit analogue to digital converter. Assuming an effective dynamic range of 21 bits, the least significant bit corresponds to 10  μV, corresponding to a magnetic field resolution below 100  fT (given a coil sensitivity of 135  mV/nT). This value is still higher than the intrinsic noise level of the sensors but is below ambient noise fields and thus does not impose a limitation on the data quality.

The sensitivity of the fluxgate sensor is 0.1  mV/nT. The resolution of the AD converter then corresponds to 100  pT only. Therefore, we amplify the sensor output with a gain factor of 64, resulting in a digital resolution of approximately 2  pT, which is slightly less than the 6  pT noise level of the fluxgate. Consequently, we have to remove the DC component of the magnetic field with an analog high-pass filter prior to digitization to remain within the input range of the AD converter. Then, the amplified fluxgate sensor output voltages are within a few volts, mainly due to motion in the geomagnetic field.

Coupling between the obliquely installed induction coil sensors components was anticipated, due to induced magnetic fields when the coils move in the geomagnetic field, and due to feedback currents within each coil. However, neither ground test measurements nor the results from airborne experiments indicate that coupling plays a role, given that the induction coil responses and the fluxgate responses are compatible at low frequencies.

A technically limiting factor remains the recording of platform motion. Although the INS installed within the bird offers one of the highest attitude resolutions available, the accuracy is still below the sensitivity of the magnetic field sensors. It means that motion-induced voltages, particularly at higher frequencies, can not entirely be removed. Future improvements may therefore consider more sophisticated data fusion processing schemes that allow treatment of the INS data also as noisy observations (e.g., an errors-in-variables approach), as opposed to the correction scheme used in this work. However, we do not observe systematic distortions of the processing result due to motion, as exemplified by comparison of transfer functions with those of an overflown ground station.

Together with improved motion noise removal, a further noise reduction could be achieved with increased suspension lengths of the bird from the helicopter. Hover tests with the bird on the ground were conducted (not shown in this paper) and have revealed that broadband EM signals, which overlap with the signal frequencies regarded here, are visible even at helicopter distances in excess of 100  m.

As opposed to most airborne systems, we process the data in the frequency domain rather than in the time domain. Data processing consists in the straightforward estimation of univariate transfer functions between the magnetic fields and the injection current. It is mandatory to have full waveform recordings of the current and accurate time synchronization. The imaginary part of the transfer function is entirely induced (or scattered) and any synchronization error would distort the information contained in the responses.

Short window lengths are required for airborne data processing as a trade-off between spatial and spectral resolution, thus data quality. We have tested various parameter ranges and find it useful to combine short windows within the vicinity of the transmitter with longer temporal windows at greater distances. Overall, our scheme suggests that good data quality can be achieved with our system for frequencies from approximately 10  to  5000  Hz for the survey geometry described here. Generally, the most challenging part of semi-AEM data acquisition and processing is to obtain good-quality low-frequency responses, that is, at the base frequency, where the largest penetration depths can be achieved.

AEM surveying is a key technology for the exploration of mineral deposits. The challenge to explore deeper than presently possible with AEM surveys may be tackled with the semi-AEM approach. The semi-AEM method is conceptually a variant of controlled source EM, where the receiver is towed through the air and the transmitter is deployed on the ground. This removes any technical limitations on the strength of an EM transmitter mounted on an aircraft, and it facilitates the realization of large transmitter-receiver offsets.

We have laid out the technical details of a new semi-AEM system that has been developed for the purpose of exploration of deeply seated mineral deposits. The primary data products of our system are frequency- and distance-dependent transfer functions between the airborne recorded field components and the synchronously recorded injection current. In a test survey, we were able to obtain high-quality transfer functions over a frequency range of approximately 10–5000 Hz with an areal coverage of approximately 20  km2, centered at a 1  km long bipole transmitter on the ground. The 3C magnetic field recordings from highly sensitive broadband induction coils with an integrated feedback system were processed together with platform orientation measurements from a high-precision INS to remove the distortion from motional noise and to estimate transfer functions. Simultaneous fluxgate measurements were processed as independent recordings and confirmed the induction coil processing results at low frequencies. The induction coils developed for this airborne system are clearly superior to fluxgate sensors across the EM bandwidth used in this experiment.

Transfer function estimates reflect the source field characteristics, but their spatial decay and frequency dependency are controlled by induction phenomena in the conductive ground. These data are the input data for inversion to recover the electrical conductivity distribution within the subsurface. In practice, multiple transmitters should be deployed on the ground for real case studies and should be overflown in multiple overlapping flights to achieve good spatial coverage. This involves the need to access the ground and thus poses additional effort in terms of logistics and manpower when compared to AEM surveys; however, this extra effort may be worthwhile when deep targets are to be investigated.

We are grateful to Y. Li and three anonymous reviewers, and to the associate editor A. Swidinsky, for constructive comments that greatly helped to improve the manuscript. This research has been funded by the Federal Ministry of Education and Research of Germany in the framework of DESMEX (project nos. 003R130A-E).

Data associated with this research can be requested by contacting the corresponding author. Release of data is subject to the approval of all project partners.

Freely available online through the SEG open-access option.