To construct a training set that closely resembles the real subsurface model, we introduce a Gaussian random rough surface (Kobayashi et al., 2002), whose power spectrum is described by P k x , k y = δ 2 l x l y 4 π exp − k x 2 l x 2 + k y 2 l y 2 4 (5) where δ is the root-mean-square height of the rough surface that determines the scale of the rough surface in the vertical direction. k x and k y are the wave numbers of the rough surface in the x and y directions, l x and l y are the correlation lengths along the x and y directions that are related to the scale of the rough surface in the horizontal direction. For rough surfaces with horizontal scales of l x and l y , a larger root mean square height means a rougher surface. Therefore, we simulate the random shape of the subsurface anomalies and the undulating interfaces by the Gaussian random rough surfaces.

To train our network, we use MATLAB to establish certain number of resistivity models and use the staggered finite-difference method to calculate the EM responses. For this purpose, we write the governing equation for the secondary electric field as ∇ × ∇ × E s = − i ω μ σ E s − i ω μ σ − σ b E b (6) where the time harmonic factor e i ω t has been assumed. E b and E s are the primary and secondary field, respectively, σ is the underground conductivity (equal to the reciprocal of the resistivity), σ b is the conductivity of a background half-space, μ is the vacuum magnetic permeability. We apply the spatial differences to replace the derivatives in Equation 6 and discretize the governing equation in the frequency domain, and then we use the quasi-minimal residual (QMR) method to solve the linear equations system and obtain the frequency-domain AEM responses. Finally, we match the calculated responses and the corresponding resistivity models into sample pairs to construct the desired training set.

The 2D resistivity model we establish in this paper is on the x , z -plane. Considering that the footprint of a frequency-domain AEM system is generally very small, we set up the resistivity models with 16 × 32 rectangular grids. The thickness of the first layer is 5m, and the thickness of each layer increases by a factor 1.12. The grid size in the horizontal x direction is 30 m. The parameters of the AEM system are chosen to be consistent with those used for the survey data. Table 1 gives the frequency and transmitter-receiver (T-R) offsets. Assuming that the survey line is aligned in the x direction, the interval between the measuring points is 30 m. We calculate AEM responses of five frequencies at each measuring point, so that we can construct sample pairs with multi-frequency data from all measuring points along a survey line and the corresponding underground resistivity model of 16 × 32 parameters.

To obtain better training performance and accelerate the convergence, it is necessary to preprocess the training set before putting it into the neural network for training. Here we adopt the logarithmic normalization to keep the range of training data within [0, 1]. Specifically, for the resistivity parameters, we have Y i = log Y i log Y i max + a (7)

while for AEM data, we have X i = log X i + b log X i max + c (8) where the parameters a , b , c are introduced to ensure that the normalized parameters and data range from 0 to 1. These parameters should be chosen based on the ranges of the model parameters and the data in the training set. We then put the normalized training set into the training process shown in Figure 1. After that, we use an Adam optimizer (Diederik and Jimmy, 2014) to update the weights of the network and minimize the loss function to obtain the best training performance. In the network training process, we set the initial value of learning rate to 10 − 4 and the number of samples per batch to 32. The loss function used here is defined as the mean square error (MSE) of the difference between the model parameters and the predicted ones, i.e. M S E = 1 n ∑ i = 1 n Y i − Y ^ i 2 (9) where Y i and Y ^ i represent respectively the true and the predicted resistivity model, n denotes the number of samples in each batch.

In the above network training process, we use the TensorFlow (Abadi et al., 2016) to establish the U-net network structure and perform the entire training process on a platform for data modeling and analysis at https://www.kaggle.com/. The GPU version accessed by this platform is Tesla P100-PCIE-16GB. Taking advantage of these supporting conditions, the training time required to loop a training set containing 100,000 samples for one epoch is approximately 60 s.

The U-net has a “U" network structure. It was first used for the medical image segmentation. After that, it has also been successfully used in seismic data processing, inversion and interpretation in recent years (Yu and Ma, 2021). This method improves the traditional full convolution neural network (FCN) (Shelhamer et al., 2015). Its unique convolution operation can not only effectively reduce the large memory consumption caused by too many network layers, but can also largely reduce the number of weights and bias in the network, and thus ease the over fitting to the data.

The U-net architecture used for our training task is shown in Figure 2. Since the network input has a dimension of 32 × 10 (including 32 measuring points and five frequencies with real and imaginary responses), we interpolate the input data into 32 × 16 , and then transpose it to match the output dimension of 16 × 32 . Considering that the sizes of our input and output data are small, too many down-sampling operations may lead to the loss of image information, which will reduce the positioning accuracy in the up-sampling process, we remove two down-samplings and up-samplings operations in the original U-net architecture. However, to solve the problem with insufficient extraction of high-level features due to fewer down sampling and convolutions, we added two convolution operations after each down-sampling.

As shown in Figure 2, when we put AEM responses as input data into the trained U-net, the corresponding geoelectric model will be output after the forward propagation process. The numbers above and below each blue rectangle represent the number of channels and the data dimension of the current layer. The network used in our imaging process takes 64, 128, and 256 channels. The blue arrows represent the convolution operation shown in Figure 3. In each layer of convolution, we select the ELU (Exponential Linear Units) function (Clevert et al., 2015) as the activation function, so that the network has non-linear fitting capability. The red arrows represent the maximum pooling operations, as shown in Figure 4. The green arrows represent the up-sampling operations, which are the operations of data recovery via interpolations. The grey arrows represent the skip connections. This operation introduces the feature information on corresponding scales in the contracting path into the up-sampling process, so as to achieve feature fusion at different levels for a more precise imaging. It should be pointed out that in order to reduce the loss of edge information, we fill in all zeros to outside grids in each convolution operation to ensure that the data dimension will not be changed during the convolution process (see Figure 3).

In the following, we first validate our 2D imaging algorithm by synthetic data. As we know, the training set determines the performance of the network. The more complex the training set is, the more complex geological model the network can predict. Therefore, in this paper we establish two training sets with different complexities. The simple training set contains 50,000 samples, while the complex training set contains 100,000 samples. The model sizes of the two training sets are the same, both have 16 × 32 grids, but the simple training set has uniform background and contains only one abnormal body. For the complex training set, we divide the background into upper and lower layers, with the layer interface being randomly set. In addition, we put up to three (0–3) abnormal bodies randomly in the layered background. As for the resistivities, we assign the background resistivity randomly within 1 Ω⋅m—10,000 Ω⋅m (logarithmically 0–9.21) and the abnormal resistivities within 1 Ω⋅m—1,000 Ω⋅m (logarithmically 0–6.91). Then, we randomly combine them into two cases: a conductive background with resistive anomalous bodies, and a resistive background with conductive anomalous bodies. After that, we take 1,000 samples from each training set to test the imaging effect of the network. The resistivity models used for testing do not participate in the network training.

After numerous tests, we find that 200 training epochs can deliver good imaging results. Figure 5 and Figure 6 respectively show the attenuation process of training errors and imaging results of the two training sets. From Figure 5 it is seen that the curves as a whole show a downward trend, the two curves of training set and validation set fit well with each other at the early training stage. However, as the training epochs increases, the loss error curves gradually get flattened. The minimum loss errors for the complex and simple training set reduce to .0015, .0017, respectively. From the imaging results shown in Figure 6, one sees that ether for a simple or more complex model, the neural network can very accurately establish the mapping relationship between the AEM responses and the underground resistivities, the layer interfaces and the boundaries of the abnormal bodies are clearly revealed, the predicted resistivities are also very close to the true values.

Model1∼Model4 are for a uniform background with a single body; Model5∼Model8 are for a complex background with multiple bodies.

To demonstrate the effectiveness of our imaging algorithm, we compare in Figure 7 the imaging results of theoretical data with the traditional inversions. From the comparison, we can see that the imaging results can clearly reveal the layer interfaces and anomalous boundaries, and deliver an image very close to the true model. The inversion results of Gauss-Newton (GN) method with a uniform half-space as the initial model can only roughly reveal the shape of the anomalous boundaries, the inverted resistivity is not accurate. Furthermore, the time consumptions for the GN inversion and our imaging method are quite different. While a single iteration in GN inversion takes several minutes on DELL workstation of Intel R) Xeon(R) Gold 6256 CPU at 3.60 GHz +3.59 GHz, while our imaging takes only seconds on the same equipment.

In airborne EM, the flight altitude has a very significant impact on EM responses. To analyze the influence of flight altitude on the imaging of our network, we take 50,000 samples from the complex training set in Section 3.1 and calculate EM responses for four flight altitudes of 30, 50, 70, and 100 m for each model, so that the training set of 50,000 samples is expanded to 200,000. During the training process, we take the flight altitude as input in our network. After the network is trained for 200 epochs, the loss decreases to .0011 (see Figure 8).

To check the effectiveness of our imaging method, we test the imaging results at three different altitudes of 40, 50, and 80 m. Among them, the flight altitude of 50 m has been seen by the network, while the altitudes of 40m and 80 m are not included in the training set and thus are unfamiliar to the network. Figure 9 shows the imaging results of AEM data for three flight altitudes. It is seen that no matter whether our network is familiar or unfamiliar with the flight altitude, it can obtain very accurate imaging results. However, when the network is familiar with the flight altitude the imaging results are better than those of the opposite cases. We can hope that with more AEM data for different flight altitudes being added to the training set, we will be able to improve the imaging results even for unfamiliar flight altitudes.

Until now, our imaging results are obtained from the theoretical data. However, AEM survey data are often contaminated by noise, so a stable and practical imaging method must have certain anti-noise capability. As we know, the neural network is sensitive to changes in the input data, and generally the longer the training time of a network is, the less adaptable to the interference of noise and the weaker the stability of the network is. Therefore, in this section we add different levels of Gaussian random noise (1%, 5%, 10%) to the test datasets used in the complex training set in Section 3.1. In addition to the aforementioned 200 epochs of model training, we train the training set for extra 50 epochs (the minimum error of the validation set is .0043). Then, we use the network to image the test set with different noise levels and analyze the influence of noise on the imaging results of the network. From the imaging results shown in Figure 10, we can see that when different levels of Gaussian noise are added to the test set, the network can characterize the boundaries of the anomalous body and the layer interfaces. This indicates that our U-net has certain anti-noise capability and can recover the model at different noise levels. However, analyzing the noise immunity of the network for different training epochs, we find that compared to the network with fewer training epochs, the network adaptability to noise with more training epochs is reduced. Although the network that underwent 200 training epochs achieves good imaging results when no noise was added, with increasing noise level, the depiction to the lower boundaries of the anomalous bodies begins to distort, and the layer interfaces begin to become blurred. In contrast, although the imaging results for 50 training epochs become slightly worse for noise-free data, they are less affected by the noise in the data.

Summarizing the above discussions, we can draw the conclusion that actually there exists a balance between the stability and accuracy of imaging results when using the network to image AEM data. As the noise level in the practical survey data cannot be accurately estimated, numerical experiments need to be done before starting processing AEM data, so that good imaging results can be obtained.