## ABSTRACT

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\xd75\u2009\u2009km2$ for a single source installation. The system has the potential to be used for focused exploration of deep targets.

## INTRODUCTION

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\u2009\u2009\Omega 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 $400\u2212600$ 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 km^{2} 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).

## PRINCIPLE OBJECTIVES OF THE SEMI-AEM METHOD

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.

## SYSTEM DESIGN AND INSTRUMENTATION

### Ground transmitter

An example recording of the source current — injected through transmitter TxB4 (see Figure ^{2}) — with a strength of $\xb120\u2009\u2009A$ is depicted in Figure ^{4}. We display a $1\u2009\u2009s$ 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\u2009\u2009s$ 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,\u2026$ in the frequency range from $10\u2009\u2009to\u2009\u20098000\u2009\u2009Hz$. 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\u2009\u2009m$ 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\u2009\u2009Hz$, a dynamic range of $\xb178\u2009\u2009nT$, a noise level of $30\u2009\u2009fT/Hz$ at $10\u2009\u2009Hz$, and a sensitivity of $135\u2009\u2009mV/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).

^{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,\u2009y^b$, and $z^b$. We refer to the three coils as $u$, $v$, and $w$, arranged in the $u^\u2212,v^\u2212$, and $w^\u2212$directions. Then, these oblique sensor axes can be expressed in the body coordinate system as

^{5}. In detail, coils $u,v,and\u2009w$ are rotated by $\alpha u=16\xb0,\alpha v=16\xb0$, and $\alpha w=\u221230\xb0$ (positive clockwise) out of the long axis of the bird (the $x^b$-axis) and dipping by $\beta u=26.57\xb0$, $\beta v=\u221226.57\xb0$, and $\beta w=0\xb0$ (positive downward). Then,

The Bartington MAG-03 provides a bandwidth of $3\u2009\u2009kHz$ with a flat response up to 1 kHz (see Figure ^{6}) and a noise level of 4–6 $pT/Hz$ at $1\u2009\u2009Hz$. The sensitivity is $0.1\u2009\u2009mV/nT$. The three-axis sensors possess an orthogonality error $<0.1\xb0$. 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\u2009\u2009Hz$. Attitude data are internally processed, and outputs include attitude (the accuracy is $0.01\xb0$ 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.

## TIME SERIES

### Removal of motional noise

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\u2009\u2009Hz$ (Figure ^{6}). With respect to the latter condition, we observe dominant pendulum motions of the bird excited in the period range between $1\u2009\u2009and\u200920\u2009\u2009s$, 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\u2009\u2009km$; the transmitter has been traversed at a relative time of $50\u2009\u2009s$. 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 $20\u2009km2$, 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 $\xb11000\u2009\u2009nT$ ($\xb1100\u2009\u2009mV$ sensor output) for the fluxgate readings and $\xb110\u2009\u2009nT$ ($\xb11350\u2009\u2009mV$ 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\u2009\u2009Hz$, 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\u2009\u2009Hz$ using spline interpolation. Interpolated magnetic field predictions using the INS data are low-pass filtered with a corner frequency of $200\u2009\u2009Hz$ 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\u2009\u2009ms$, 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\xb0$ (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.

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.

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.

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 $\xb1100$ to $\xb110\u2009\u2009nT$ (Figure ^{8}–^{8}), but appears to be of minor importance for the induction coil data (Figure ^{8}–^{8}). Figure ^{8}–^{8} and ^{8}–^{8} 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.

## TRANSFER FUNCTIONS

### Estimation scheme

^{8}), particularly when the receiver approaches or moves away from the transmitter. Instead, the field can be thought of as an amplitude-modulated signal

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\xaf\u2009\u2009Hz=1/0.096\u2009\u2009s$, the interpolated sampling rate is $fs=20,000\u2009\u2009Hz$ and one cycle corresponds exactly to $Lwin=1920$ samples.

^{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\xb11/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

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.

## RESULTS

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\u2009\u2009m$, 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\u2009\u2009m$) or with increasing window length (Figure ^{9}, $Ncyc=32$, $Nwin=2$, averaging lengths approximately $185\u2009\u2009m$; Figure ^{9}, $Ncyc=32$, $Nwin=8$, averaging lengths approximately $550\u2009\u2009m$), 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\xaf\u2009\u2009Hz$ (Figure ^{11} and ^{11}), and for $1448.2\u2009\u2009Hz$ (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\u2009\u2009Hz$). 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\u2009\u2009s$ 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\u2009\u2009Hz$ 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\u2009\u2009pT/A$; the real parts of the vertical component $Bz$ of the airborne estimates are however systematically larger by $10\u221220\u2009\u2009pT/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 $\u22126\xb710\u22124\u2009\u2009nT/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\u2009\u2009m$ 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\xb710\u22125\xb13\xb710\u22125\u2009\u2009nT/A$) and a slope close to identity ($1.01\xb14\xb710\u22124$, 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 $\Delta 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 $\sigma \Delta =6\u2009\u2009pT/A$ for the real part and $\sigma \Delta =4\u2009\u2009pT/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\u2009\u2009Hz$ along two profile lines, L22 and L26. The data are shown for two independent transmitters, TxB4 and TxB3, which are located at $0\u2009\u2009km$ and at $\u22121\u2009\u2009km$ profile distances, respectively. Profile line L22 crosses both transmitters, whereas profile line L26 passes the transmitter at a distance of approximately $500\u2009\u2009m$. 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\u2009\u2009\Omega 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\u2009\u2009\Omega 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 $\Omega 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.

## DISCUSSION

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\u2009\u2009\mu V$, corresponding to a magnetic field resolution below $100\u2009\u2009fT$ (given a coil sensitivity of $135\u2009\u2009mV/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\u2009\u2009mV/nT$. The resolution of the AD converter then corresponds to $100\u2009\u2009pT$ only. Therefore, we amplify the sensor output with a gain factor of 64, resulting in a digital resolution of approximately $2\u2009\u2009pT$, which is slightly less than the $6\u2009\u2009pT$ 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\u2009\u2009m$.

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\u2009\u2009to\u2009\u20095000\u2009\u2009Hz$ 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.

## CONCLUSION

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\u2009\u2009km2$, centered at a $1\u2009\u2009km$ 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.

## ACKNOWLEDGMENTS

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 AND MATERIALS AVAILABILITY

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.