The magnetotelluric (MT) method uses the earth’s natural electromagnetic signals to image the resistivity distribution of the subsurface. Its application is widespread in studies to improve our understanding of the lithospheric architecture, earthquake zones, and for the exploration of natural resources, such as critical minerals and geothermal energy. Inversion of MT data is a fundamental step in the data analysis routine to retrieve a subsurface geo-electrical model that can be used to inform geological interpretations. In order to reduce the effect of non-uniqueness and the local minimum trapping problems of conventional deterministic inversion, a data-driven mathematical method with a deep neural network can be used to estimate the subsurface properties in many geophysical inversion techniques. In this study, we investigate a deep learning (DL) inversion method composed of a multi-head convolutional neural network (CNN) architecture for 1-D MT data analysis. We created synthetic datasets consisting of 100,000 random samples of resistivity layers to train the CNN model parameters. The trained model was validated with independent synthetic datasets, and the predicted resistivity distribution displayed an acceptable resolution and reliability, which demonstrate the potential application of DL inversion for MT data. The trained CNN model was used to analyze real MT data collected in the southwestern Athabasca Basin, Canada. Predicted results from the DL method displayed a similar subsurface resistivity distribution with traditional iterative inversion. Because this approach can predict a resistivity model without multiple forward modelling operations after the trained CNN model is created, this framework is suitable to speed up multidimensional MT inversion for predicting subsurface resistivity.
Research Article
Deep learning based one-dimensional inversion of magnetotelluric data, and an application in the southwestern Athabasca Basin, Canada
https://doi.org/10.21203/rs.3.rs-1599373/v1
This work is licensed under a CC BY 4.0 License
Version 1
posted
You are reading this latest preprint version
magnetotelluric
deep learning
resistivity
CNN
inversion
An inverse problem is the process of predicting the causal factor from the outcome of measurements, given a partial description of a physical law-based system. Application of inverse theory is widely applied in the science and engineering disciplines, such as geophysics and computer vision (Kim and Nakata, 2018). Inversion of magnetotelluric (MT) data is one such application that uses surface measurements of the earth’s natural electric and magnetic fields as inputs to retrieve the subsurface resistivity structure for the interpretation of subsurface geology. The forward modeling procedure is based on the calculation of electromagnetic (EM) fields from a subsurface geo-electric model using numerical solutions of Maxwell’s equations (e.g. Ward and Hohmann, 1987; Zhdanov, 2002; Craven et al., 2006). Different inversion optimization methods, such as nonlinear conjugate gradient and Gauss-Newton have been studied previously (e.g. Constable et al., 1987; Zhdanov, 2002; Egbert and Kelbert, 2012). The main principle of these approaches is the use of an iterative optimization technique to minimize the parametric function, including a data misfit function and a stabilizing function. Some useful software packages have been developed successfully for use by the MT community in academia and industry, such as Occam and ModEM (Constable et al., 1987; Egbert and Kelbert, 2012; Kelbert et al., 2014). Although there are some advantages inherent to traditional inversion methods, such as delivering details of the electrical structure of the subsurface, these inversion methods still face a few challenges. The first issue is that an initial model is needed to avoid the objective function trapping the search in a local minimum. Additionally, forward modeling, especially for 3D geo-electric structure, is numerically challenging, requiring significant computational resources despite recent advances (Ansari et al., 2020) to avoid a lengthy computation time. This issue is exacerbated for iterative inversion schemes that creep toward a final model using multiple forward computations. It must also be borne in mind that non-linear inversion of band-limited non-precise data is unstable and non-unique (Zhdanov, 2002; Zhdanov et al., 2011). Although regularization methods and physical constraints may be employed to stabilize the inverse solution, an alternative technique such as deep learning (DL) that leverages modern GPU hardware and algorithms should be investigated to assess its ability to replace traditional inversion techniques.
The DL method has matured in the past few decades (Fukushima, 1980; Hochreiter, 1998), and more recently, increased in popularity in different disciplines (Chen et al., 2020; Long et al., 2014; Ronneberger et al., 2015) along with leaps in computer hardware development ideally suited for the underlying DL algorithms. As an increasingly popular method for data-driven inferences in geoscience disciplines, artificial neural networks (NN) have demonstrated a considerable number of successful applications in different areas. For example, similar to image segmentation in self-driving cars (Paszke et al, 2016) and medical image analysis (Ronneberger et al., 2015), a NN framework has been used to aid in the extraction of clay mineralogy and fracture mapping from scanning electron microscope images with U-net and Deeplab architecture (Chen et al., 2020). For optimization problems in geophysics, previous studies have shown DL and machine learning can be used to solve inverse problems. Kim and Nakata (2018) and Russell (2019) compared machine learning and geophysical inversion, and showed machine learning yields results with high spatial resolution. Das et al. (2019) deployed convolutional neural networks (CNN) for seismic impedance inversion and demonstrated promise in the performance of more accurate and faster seismic reservoir characterizations in comparison to traditional inversion techniques. The DL inversion procedure is computationally instant after the model is created, which can improve the inversion efficiency. To overcome the drawbacks of traditional EM inversion methods, it is worth investigating a DL architecture to retrieve the unknown geo-electric structure of the subsurface. Moghadas (2020), Puzyrev (2019), and Puzyrev and Swidinsky (2021) presented 1D inversion of EM induction data, controlled source EM and time domain EM (TEM) data using CNN. The two studies displayed impressive inversion results using a CNN architecture, suggesting the potential of CNN inversion for more widespread use in EM exploration. Liu et al. (2020) presented a sample-compressed neural network algorithm for MT inversion and an adaptive-clustering analysis algorithm for resistivity boundary demarcation, and Guo et al. (2020) employed a supervised descent method (SDM) for 2D MT data inversion to reduce uncertainty. Conway et al. (2019) applied a non-convolutional neural network model to approximate the 3D MT forward method instead of resolving the expensive forward functions, speeding up the inversion procedure and producing reasonable resistivity models. In spite of these advances, further study of MT inversion using CNN has yet to be performed.
Because a CNN may fail to respond optimally to the entirety of scales inherent to the training data, the task of applying CNN model training can be challenging. In order to reduce the number of datasets required for CNN training and to improve predictive capabilities and flexibility, Long et al. (2014) and Roy et al. (2018) have developed fully convolutional networks by re-architecting and fine-tuning classification nets to use direct dense prediction of semantic segmentation, and Ronneberger et al. (2015) have modified the fully convolutional network by adding more feature channels in the upsampling part, and applied the U-net architecture for medical image analysis. Based on previous study of CNN architectures for semantic image segmentation, the DeepLabv3 + architecture with decoder modules were applied by Chen et al. (2017) using an atrous convolutional operator with encoder-decoder network structure for solving the segmentation task.
In this study, we implement a multi-head CNN with residual block, which permits the use of measured MT resistivity and phase data and facilitates working in a multi-scale context. We introduce a workflow of dataset generation, model training, and validation, and then assess this CNN inversion method with a study of synthetic and real data from Athabasca Basin. We also test the sensitivity of our method to data noise and compare the results with traditional inversion method, which shows the advantage of the proposed CNN inversion method.
The MT method retrieves the subsurface distribution of electrical resistivity by measuring the variations of the EM field at the earth’s surface, over which a source signal is assumed to be a naturally-occurring planar EM wave. The recorded electric and magnetic fields at the surface of the earth contain resistivity information on the subsurface, and their ratio in frequency is termed an impedance tensor,
$${Z}\left(\omega \right)=\frac{{E}\left(\omega \right)}{{H}\left(\omega \right)}$$
1
where \({E}\) is the electric field strength, \({H}\) is defined as the magnetic field strength, and \(\omega\) is the angular frequency (Ward and Hohmann, 1987; Zhdanov, 2002). The frequency dependent tensor can be written in Cartesian coordinates:
$$\left[\begin{array}{c}{E}_{x}\\ {E}_{y}\end{array}\right]=\left[\begin{array}{cc}{Z}_{xx}& {Z}_{xy}\\ {Z}_{yx}& {Z}_{yy}\end{array}\right]\left[\begin{array}{c}{H}_{x}\\ {H}_{y}\end{array}\right]$$
2
Assuming natural EM fields travel in terms of 1D resistivity media, the components of impedance tensor can be simplified as:
$$\left\{\begin{array}{c}{Z}_{xx}={Z}_{yy}=0\\ {Z}_{xy}=-{Z}_{yx}\end{array}\right.$$
3
The apparent resistivity and phase are commonly used in 1-D data analysis. The phase is the lag between the\({E}\) and \({H}\) fields, whereas the apparent resistivity (\({\rho }_{a})\) is a normalization of the field amplitudes and can be expressed as:
$${\rho }_{a}\left(\omega \right)=\frac{1}{\omega \mu }{\left|{Z}_{xy}\left(\omega \right)\right|}^{2}$$
4
where \(\mu\) is the magnetic permeability.
It therefore represents an average resistivity between the earth’s surface and skin depth,
which is function of conductivity\(\sigma\), magnetic permeability \(\mu\), and angular frequency\(\omega\). Figure 1 shows the apparent resistivity and phase for various 1D models. The Occam1D forward modelling code is used to calculate the response of the three layer geo-electrical models (Constable and Weiss, 2006). We can see that when the resistivity decreases, the phase is higher than 45o; and when the apparent resistivity increases, the phase switches to below 45o. Based on this relation between apparent resistivity and phase, a consideration of both responses as input will help to achieve high generalized DL model. For a model generating a 1D CNN inversion, the trained DL neural parameters will be calculated based on a group of a randomly generated models. The simulated apparent resistivity and phase are gathered as the input and the resistivity model as the output for calculating the weight parameters during model training.
Geophysical inversion attempts to estimate the distribution of physical properties of the subsurface from surface observations. The proposed workflow implements a significant concept of DL with CNN to calculate the resistivity distribution from the measured data directly. Similar to the usual DL procedure, the flowchart schematic in Fig. 2 shows the major components and CNN architecture of a MT inversion consisting of three modules: data generation; neural network training and validation; and prediction.
3.1 Training dataset generation
The quantity and quality of the dataset for training the neural network model parameters is important for prediction accuracy. As the left box in Fig. 2 shows, the dataset (apparent resistivity and phase) were generated using this workflow and split into training set and validation set by an 80/20 ratio. The random resistivity models were created for a 50-layers model with a limited resistivity range in [1–1E4 Ωm]. The resistivity distribution in the model was randomly generated with a maximum of 10 major anomalous continuous resistive or conductive zones. The thickness of the layers was fixed for all models and started at 20 meters increasing in thickness to 10,000 meters for the bottom-most layer. The maximum depth is 100 km. 128 frequency points were calculated within a range of 100–0.001 Hz. The forward module in the Occam inversion package was used for generating the apparent resistivity and phase (Constable and Weiss, 2006). The large amount of variably shaped resistivity models will cover most of the possible real world subsurface geo-electric models without using any data augmentation methods. For this study, we generated 100,000 layered resistivity samples for model training and validation, which took about 40 minutes using a 56-cores CPU desktop. The base 10 logarithm of resistivity \(m=log10\left(\rho \right)\) was used as the targeted model parameters, and the input data were standardized using mean \(\mu\) and standard deviation \(\sigma\) values on the training processing (Puzyrev and Swidinsky, 2021):
$$\tilde{x}=\frac{x-\mu }{\sigma }$$
6
3.2 CNN training and validation framework
The most important part of this inversion procedure is calculating the weight parameters of the network model that includes training and validation. Since apparent resistivity and phase are different features and need to be input for training, the popular CNN approach named multi-head model architecture is employed for this study (Fig. 2) (Canizo et al., 2019; Khan and Ahmad, 2021). As only 128 frequency points were used for this model framework, a feature map with a convolution kernel size of 3 was used. The framework includes two parallel heads for feature analysis of standardized apparent resistivity and phase. Each head was designed using similar 1-D CNN layers but different weights, to extract features from the corresponding input, then the interpretations of both heads were concatenated into one interpretation as a fully-connected layer (Canizo et al., 2019; Khan and Ahmad, 2021). For each convolution layer, the number of filters increased from 32 to 128 with a rectified linear (ReLU) activation function and a stride of 2, which is followed by 0.5 dropout and max-pooling layers. The ReLU is a linear function, \(f\left(x\right)=\text{m}\text{a}\text{x}(0,x)\), which has the advantage of overcoming the vanishing gradient problem in comparison to the sigmoid function. The stride is a parameter to control the filter movement; the max-pooling layer is a method to down sampling feature maps by calculating the maximum value for each patch of feature map, and the dropout layer is a regularization method that sets some hidden layer neurons to zero randomly during each epoch. The main function of the dropout layer is to prevent model overfitting in order to get a better convergence of validation. The details regarding the definition and function of these neural network layers are based on the Keras API (Chollet et al., 2015). The increased number of filters allows the neural network to propagate context information to complex layers for a better understanding the features. More details on the network model can be found in Long et al. (2014), Ronneberger et al. (2015), Chen et al. (2020), and Puzyrev (2019). During the training and validation epochs, the network parameters are updated by minimizing the objective function, which is specified with the root mean squared error (RMSE) as the loss function metrics,
$$rmse=\sqrt{\frac{1}{n}\sum _{1}^{n}{({y}_{i}^{pre}-{y}_{i}^{mod})}^{2}}$$
7
\({y}_{i}^{pre}\) is predicted value and \({y}_{i}^{mod}\) is input model. The objective function of known resistivity models and their corresponding forward results is used to train the network with the Adam (adaptive moment estimation) optimization algorithm, which is an extension of the stochastic gradient descent for first-order gradient-based optimization that has recently seen broader adoption for DL applications (Kingma and Ba, 2017). Adam combines the advantages of both the Adaptive Gradient Algorithm and the Root Mean Square Propagation, such that it is able to deal with sparse gradients and non-stationary objectives. One important property of Adam’s update rule is to choose bounded step sizes using the step size hyperparameter, as it renders Adam well-suited to a wide range of objective function optimizations in the machine learning field. Python and some relevant DL open-source libraries, e.g. Tensorflow, Keras, are used for realizing the procedure (Chollet et al., 2015; Abadi et al., 2016), and MTpy (Krieger and Peacock, 2014) and Occam1D package are used for model generation and output display.
3.3 Model test and elements selection
The network model training and validation is a key step to calculate the model parameters and reduce overfitting. In this study, the model training took about 5 hours of computation time with a batch size 64 using one NVIDIA Quadro-P5000 GPU in order to reach a stable and acceptable accuracy; the model with best validation loss is saved during maximum 150 epochs. Figure 3 shows an example of convergence during model training and validation—we can see the convergence stability for both the training set and validation set. Because 1D CNN MT inversion is a data-driven technique, we need to deepen the test to the created model to make sure it is accurate and stable. The extra random models are generated with different parameters, such as layer thickness and kernel size to test the CNN model. If the RMSE and the correlation coefficient of the test result are not satisfied, different sample models will be added for training the network. Although the core skeleton is similar for all neural network frameworks, there are some difference in shapes and structures. Before starting model training, the elements selection are a key element for high accuracy prediction. In order to improve the RMSE accuracy, the different hyper-parameters (e.g. network structure, number of hidden layers, batch size, filter number and kernel sizes, activation function, and optimization function) were adjusted during model training period. In addition, we also checked different learning rates for model updating. We chose four hidden layers with a fixed kernel size 3 for the final model. After the tests were verified, the same network architecture and weight parameters were saved and used in the rest of the program.
4.1 Results of synthetic model
In order to evaluate the functionality and efficiency of the trained network model, we created synthetic models and calculated apparent resistivities and phases for model testing. The purpose of model testing is to demonstrate if an accurate result can be obtained for a general resistivity distribution. The MTpy and Occam1D inversion packages were used to create the input files and conduct the forward modeling, respectively (Constable et al., 1987; Krieger and Peacock, 2014). In Fig. 4, the resistivity model samples were generated with a maximum 10 major anomalous resistive or conductive layers, selected randomly and not used during model training/validation period. One can see the variation of resistivity models, which cover different distributions and are diversified in order to create a generalized network model. Figure 5 shows a comparison of resistivity model and predicted inversion results with the trained neural network model. We can see most of the predicted results identified the deep anomalies and displayed the desired trend of resistivity variation, however a few predicted areas show large bias, such as in the complex resistivity distribution zone (e.g. fourth model on the top row). In general, the calculated RMSE and correlation coefficient (\(R\)) from the model testing are in acceptable range. Most of the misfits reach 0.05 and \(R\) are larger than 0.97 (Table 1).
| ID | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15 | 16 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| RMSE | 0.0316 | 0.0539 | 0.0708 | 0.1236 | 0.0193 | 0.0574 | 0.0499 | 0.0363 | 0.0291 | 0.0474 | 0.0519 | 0.0518 | 0.0474 | 0.0414 | 0.0341 | 0.0330 |
| \({R}\) | 0.9961 | 0.9917 | 0.9770 | 0.9139 | 0.9981 | 0.9823 | 0.9895 | 0.9921 | 0.9960 | 0.9814 | 0.9860 | 0.9904 | 0.9869 | 0.9934 | 0.9961 | 0.9943 |
During the CNN model training and validation, we used layered models with continuous variable resistivity at the boundary between layers. In order to further testify reliability of the model prediction, we created a few simplified models with discontinuous layer boundaries and compared with the results of CNN inversion (Fig. 6). We can see that all of the predicted resistivity distribution are consistent with the true models and show the trend of spatial variation, which verifies the trained CNN model works well for both continuous and discontinuous structures.
4.2 Comparison with traditional iterative inversion
In order to verify the accuracy of the proposed method, we compared the CNN predicted results with a traditional deterministic inversion method, which takes an initial model and updates the resistivity model until a best-fit is achieved with the observed data (Constable et al., 1987, 2006). 300 synthetic data were randomly generated using 1D forward simulation. During the inversion, the same layer thickness were used and the number of frequencies remained the same.
We calculated and plotted the statistic distribution of misfit and \(R\) based on testing with the synthetic models (Fig. 7, Table 2). The figure shows that the misfit RMSE distribution peak (Fig. 7 (left)) of the traditional inversion is higher than the CNN inversion in some cases, and the correlation coefficient peak (Fig. 7 (right)) is lower as well. We also can see the results of traditional inversion displays a larger variation than the CNN, and quartile lines is wider than CNN inversion, which indicates the trained CNN model produces geo-electric structure with less uncertainty. In general, the \(R\) and RMSE in Table 2 demonstrate that the resistivity models with CNN inversion are better recovered for this specific set of examples in comparison to its traditional inversion. The peak values of calculated \(R\) are 0.981 and 0.968, and the RMSE are 0.0455 and 0.0752, respectively.
| Inversion methods | Correlation coefficient (peak) | Misfit RMSE (peak) |
|---|---|---|
| CNN | 0.981 | 0.0455 |
| Traditional Inversion | 0.967 | 0.0752 |
4.3 Stability test for performance
In the study so far, the CNN architecture has been trained and validated with a noise-free dataset. In order to check the stability of this network in the presence of noise, the original resistivity and phase are noise distorted using Eq. (8). We added random noise at different levels (3% and 5%) to the simulated apparent resistivity and phase. Noise was added according to the following formula:
$${D}_{n}=D+s\bullet ({D}_{max}-{D}_{min})\bullet \left(r-0.5\right)$$
8
where \({D}_{n}\) is noised data, \(D\) is raw data, \({D}_{max}\) and \({D}_{min}\) are maximum and minimum of \(D\), \(s\) is used to adjust noise magnitude, \(r\) is a generated pseudo-random number. The noised datasets were loaded into the trained neural network model directly to calculate resistivity. In comparison to the true model, the eight inverted results show the correct trend of the resistivity variations (Fig. 8). As expected, the dataset with 5% noise exhibits a larger bias than the data with 3% noise level. E.g. the fourth testing model is noised with [0%, 3%, 5%], the \(R\) values are [0.9842, 0.9744, 0.9692], and the RMSE are [0.0549, 0.0722, 0.1037] respectively. In general, the output is quite stable in the presence of noise. In high resistivity zones, the noise had larger effect however.
5.1 Geological setting of the study area
The Athabasca Basin (Fig. 9a) straddles the border of northern Saskatchewan and Alberta and is the premier region for uranium exploration in Canada, hosting the highest-grade deposits in the world. The study area is in the western Athabasca Basin (Fig. 9b) where flat-lying Paleoproterozoic sedimentary rocks of the Athabasca Supergroup unconformably overlies Taltson domain basement rocks. There is a large range in the electrical resistivity properties of the rocks in and below the Athabasca Basin: the sandstones and underlying orthogneisses and intrusive rocks are characterized by high resistivity, with resistivities greater than 5000 Ωm on average in the eastern Athabasca Basin; silicified fault zones with high resistivities; and hydrothermal alteration haloes with elevated conductivities (Craven et al., 2007; Farquharson and Craven, 2009).
5.2. Real data inversion results
To visualize the lower crustal electrical properties, 34 broadband MT stations were collected using Phoenix Geophysics MTU-5A recorders along two northwest-southeast profiles transecting the PLC (Tschirhart et al., 2022). To test the accuracy of the trained 1D CNN model proposed in this study four sites (PLC013, PLC014, PLC015, and PLC016; Fig. 9a) were selected from the Athabasca dataset. In the upper crust, these data have roughly symmetric phase tensors suggesting 1D conductivity distributions (Caldwell et al., 2004). Geologically, at these sites, upwards of ~ 750 m of resistive Athabasca Supergroup sedimentary rocks unconformably overlies Taltson Domain basement rocks, with no major interpreted intersecting faults and therefore broadly approximates a 1-D Earth (Tschirhart et al., 2021). The determinant average of the impedance tensor at different frequencies was calculated for use in CNN model testing (Eq. 9) (Smirnov, 2009).
$${Z}_{det}=\left|{Z}_{det}\right|{e}^{i{\varnothing }_{det}}=\sqrt{{Z}_{xx}{Z}_{yy}-{Z}_{xy}{Z}_{yx}}$$
9
The predicted subsurface resistivity with the trained CNN model are presented in Fig. 10. The dotted magenta line is the forward apparent resistivity and phase of the predicted resistivity model using the CNN inversion. We can see most of the results predict the trend of the measured data. Comparing with the results of iterative inversion, both methods identified the high resistivity zone around 1 km, and displayed a similar trend in the resistivity distribution. There are some differences in the resistivity between the two methods in the near surface (< 300 m) and mantle (> 40 km) because of effect of high noise levels in the data. Both CNN and iterative inversion show similar trend of resistivity variation from 300 m to ~ 40 km. The prediction in the high frequency zone shows a larger difference potentially owing to static shift effects arising due to a resistive overburden as these sites. The deeper parts of the models is consistent across all four sites, where the resistivity decreases from 25 km to 40 km coinciding with the boundary between lower crust and mantle, which is located at ~ 35 km depth in the study area (Schetselaar and Snyder, 2017). Based on inversion results in Fig. 10, the inline classified resistivity section with site to site correlation is shown in Fig. 11. The resistivity values are displayed in logarithm space, the result of iterative inversion (Fig. 11b) is as a reference. On the whole, the section shows a reasonable correlation in the resistivity variation between sites with both methods. The shallow layers shows the resistivity increasing to a depth of ~ 1 km corresponding to the thickness of the Athabasca Supergroup cover (Tschirhart et al., 2021). The predicted response using the CNN inversion matches the resistivity variant trend shown by a previous 2D inversions (Potter et al., 2020) at some points, particularly at depths less than 20 km, which is reasonable for the potential accuracy of a one dimensional study.
In this study, we investigated the possible use of DL for 1D MT inversion. Compared to traditional inversion methods, the computation time required for this 1D MT CNN inversion is mainly spent on the synthetic dataset generation and training/validation the neural network model parameters. Based on the computational resources we had available for this study (i.e. NVIDIA Quadro-P5000 GPU), the dataset generation and CNN model training/validation stage is the most time consuming part, taking approximately five hours. The prediction stage is much faster with a single calculating step once an accepted network model is created. In contrast, traditional iterative convergence approaches with regularizing modeling constraints require updating the response after each iteration with computationally intensive forward modeling. This study displayed the potential of CNN inversion to efficiently reconstruct subsurface resistivities with less computational requirements, which is important for solving two/three dimensional problems.
Despite the advantages of using the CNN model framework and its positive potential in MT inversion from the above analysis of synthetic and real data, there are a few limitations for CNN inversion where improvement could be addressed. The generated dataset for training is not unlimited, which might lead to a bias of the related parameters in the network model. During model training period, we considered a wide range of resistivity distributions, so that better model parameters will be calculated and the predicted models are more close to real world data. To realize a more generalized DL model, one approach is to make a model generator that covers even wider resistivity ranges so that the training datasets will be better suited to a wider variety of input data. During the synthetic model study, the trained model can predict more complex subsurface structure accurately after having added more layers in the training dataset. In addition, a more complex network architecture could improve the accuracy of simulating the connection between physical models and MT observed data. Because of limitations in computer power, we only considered four convolutional layers in each head for this 1-D CNN framework. By increasing hidden layers in the network architecture and the number of filters, it will improve the noise effect of distorted data to the model. Finally, a different objective function with a stabilizer can be added to the training procedure. In general, optimizing the training dataset, the weights of layers, and a tuned objective function will make the DL inversion algorithm more powerful and efficient based on a given optimization strategy.
As noise-tolerance testing showed, the predicted results show larger bias at high noise levels. Updating the framework could be a solution to make the DL network process real world noisy data more efficiently. Based on the results of the current study, it is possible that application of this neural network framework to anisotropic or 3-D resistivity model inversion studies will be fruitful. Since a large number of datasets are required for training and validation, the training period for 3-D CNN inversion will require significant computational resources for both data generation and the calculation of weighting parameters. However, the prediction processing with a trained model will be accomplished instantaneously and be more efficient for skipping the optimization parameters adjusting of traditional iterative inversion. In this study, we demonstrated a proposed algorithm to implement the optimization problem in MT 1D inversion with a CNN framework. Continuous research is still needed to improve DL algorithms for solving high dimensional MT data optimization problems.
In this study, 1D MT inversion using a CNN architecture has been developed successfully. The merits of this technique are demonstrated by synthetic and real data testing. As a data-driven mathematical method, a trained CNN architecture can discover the relationship between observed MT data and subsurface resistivity distributions and also speed up the inversion procedure to derive reliable subsurface physical property distributions. The application of the trained network model to a noised synthetic test dataset produced satisfactory results, which suggest that developing a network for complex geo-electric model inversion is feasible. Although containing some 2D and 3D geo-electrical structure and real-world noise, the predicted resistivity model using the test stations from the western Athabasca Basin displayed reasonable geo-electric structure in comparison to results using traditional inversion algorithm. Despite there being some obvious processing advantages when using DL for geophysical inversion, there are still some drawbacks, including needing large amounts of datasets, model framework limitations, and noise tolerance bias all of which are continuously improving as the DL research space matures. In future work, we will investigate a more suitable neural network framework and workflow for multidimensional MT inversion. It will also be necessary to test the integrated objective function as a metric to improve the reliability of the predicted results, which would be a solution for solving stability and testing the sensitivity of real data applications.
Statements & Declarations
Funding
This work is fully funded by Targeted Geoscience Initiative (TGI-6) program of Natural Resources Canada (Grant number: 340333).
Conflict of interest
All co-authors declare that there is no conflict of interest in this study.
Author contributions
Xiaojun Liu is responsible deep learning method including algorithms, coding, and synthetic and real data processing, writes initial version of manuscript.
James A Craven helps collect the data and contributed to final manuscript development.
Victoria Tschirhart is responsible for MT data preparation, data QC, geological interpretation, edits and reviews manuscript.
Acknowledgements
The authors are grateful to Research Scientist Dr. SeyedMasoud Ansari at the Geological Survey of Canada for a review and discussion in an early period of this study. The NRCan Contribution number is 20210009.
- Abadi, M., Barham, P., Chen, J., Chen, Z., & Davis, A. (2016). Tensorflow: a system for large-scale machine learning. Operating Systems Design and Implementation (OSDI), 16, 265–283
- Ansari, S., Schetselaar, E., Craven, J., & Farquharson, C. (2020). Three-dimensional magnetotelluric numerical simulation of realistic geologic models. Geophysics, 85, E171–E190. https://doi.org/10.1190/geo2019-0214.1
- Caldwell, T. G., Bibby, H. M., & Brown, C. (2004). The magnetotelluric phase tensor. Geophysical journal international, 158(2), 457–469
- Canizo, M., Triguero, I., Conde, A., & Onieva, E. (2019). Multi-head CNN–RNN for multi-time series anomaly detection: An industrial case study. Neurocomputing (Amsterdam), 363, 246–260
- Chen, L. C., Papandreou, G., Schroff, F., & Adam, H. (2017). Rethinking Atrous Convolution for Semantic Image Segmentation. CoRR, abs/1706.05587. URL: https://arxiv.org/abs/1706.05587v3
- Chen, Z., Liu, X., Yang, J., Little, E., & Zhou, Y. (2020). Deep learning-based method for SEM image segmentation in mineral characterization, an example from Duvernay Shale samples in Western Canada Sedimentary Basin. Computers & Geosciences, 138(5), 104450. DOI: 10.1016/j.cageo.2020.104450
- Chollet, F., et al. (2015). Keras (webpage), URL: https://keras.io
- Constable, S. C., Parker, R. L., & Constable, C. G. (1987). Occam’s inversion: A practical algorithm for generating smooth models from electromagnetic sounding data. Geophysics, 52(3), 289–300
- Constable, S. C., & Weiss, C. J. (2006). Mapping thin resistors and hydrocarbons with marine EM methods: Insights from 1D modeling. Geophysics, 71(2), G43–G51
- Conway, D., Alexander, B., King, M., Heinson, G., & Kee, Y. (2019). Inverting magnetotelluric responses in a three-dimensional earth using fast forward approximations based on artificial neural networks. Computers & Geosciences, 127, 44–52
- Craven, J. A., Farquharson, C. G., Mackie, R. L., Siripunvaraporn, W., Tuncer, V., & Unsworth, M. J. (2006). A comparison of two- and three-dimensional modelling of audiomagnetotelluric data collected at the world’s richest uranium mine, Saskatchewan, Canada. In: 18th International Workshop on Electromagnetic Induction in the Earth, El Vendrell, Spain, pp. 1
- Craven, J. A., McNeice, G., Powell, B., Koch, R., Annesley, I. R., Wood, G. … Xiao, W. (2007). Audio-magnetotelluric studies at the McArthur River mining camp and Shea Creek area, northern Saskatchewan. Bulletin-Geological Survey of Canada, 588, 413–424
- Das, V., Pollack, A., Wollner, U., & Mukerji, T. (2019). Convolutional neural network for seismic impedance inversion. Geophysics, 84(6), R869–R880
- Egbert, G. D., & Kelbert, A. (2012). Computational recipes for electromagnetic inverse problems. Geophysical journal international, 189(1), 251–267
- Farquharson, C. G., & Craven, J. A. (2009). Three-dimensional inversion of magnetotelluric data for mineral exploration: An example from the McArthur River uranium deposit, Saskatchewan, Canada. Journal of Applied Geophysics, 68(4), 450–458
- Fukushima, N. (1980). A self-organizing neural network model for a mechanism of pattern recognition unaffected by shift in position. Biological Cybernetics, 36(4), 193–202. doi:10.1007/bf00344251
- Guo, R., Li, M., Yang, F., Xu, S., & Abubakar, A. (2020). Application of supervised descent method for 2D magnetotelluric data inversion. Geophysics, 85(4), WA53–WA65. https://doi.org/10.1190/geo2019-0409.1
- Hochreiter, S. (1998). The Vanishing Gradient Problem During Learning Recurrent Neural Nets and Problem Solutions. International Journal of Uncertainty, Fuzziness and Knowledge-Based Systems, 06(2), 107–116
- Jefferson, C. W., Thomas, D. J., Gandhi, S. S., Ramaekers, P., Delaney, G., Brisbin, D. … Olson, R. A. (2007). Unconformity-associated uranium deposits of the Athabasca Basin, Saskatchewan and Alberta. In: Jefferson, C.W., Delaney, G. (Eds.), EXTECH IV: geology and uranium exploration technology of the Proterozoic Athabasca Basin, Saskatchewan and Alberta: Geological Survey of Canada Bulletin, 588, p. 23–67
- Kelbert, A., Meqbel, N., Egbert, G. D., & Tandon, K. (2014). ModEM: A modular system for inversion of electromagnetic geophysical data. Computers & Geosciences, 66, 40–53
- Khan, Z. N., & Ahmad, J. (2021). Attention induced multi-head convolutional neural network for human activity recognition. Applied soft computing 110, 107671. URL: https://doi.org/10.1016/j.asoc.2021.107671
- Kim, Y., & Nakata, N. (2018). Geophysical inversion versus machine learning in inverse problems. Leading edge (Tulsa, Okla.), 37(12), 894–901
- Kingma, D. P., & Ba, J. (2017). Adam: A method for stochastic optimization: arXiv preprint arXiv:1412.6980. URL: https://arxiv.org/abs/1412.6980
- Krieger, L., & Peacock, J. R. (2014). MTpy: A Python toolbox for magnetotellurics. Computers and Geosciences, 72, 167–175
- Liu, W., Lü, Q., Yang, L., Lin, P., & Wang, Z. (2020). Application of Sample-Compressed Neural Network and Adaptive-Clustering Algorithm for Magnetotelluric Inverse Modeling. IEEE Geoscience and Remote Sensing Letters, 1–5. doi: 10.1109/LGRS.2020.3005796
- Long, J., Shelhamer, E., & Darrell, T. (2014). Fully convolutional networks for semantic segmentation. arXiv:1411.4038. URL: https://people.eecs.berkeley.edu/~jonlong/long_shelhamer_fcn.pdf
- Moghadas, D. (2020). One-dimensional deep learning inversion of electromagnetic induction data using convolutional neural network. Geophysical journal international, 222(1), 247–259
- Paszke, A., Chaurasia, A., Kim, S., & Culurciello, E. (2016). Enet: A deep neural network architecture for real-time semantic segmentation. arXiv:1606.02147. URL: https://arxiv.org/abs/1606.02147v1
- Potter, E. G., Tschirhart, V., Powell, J. W., Kelly, C. J., Rabiei, M., Johnstone, D. … Bethune, K. M. (2020). Targeted Geoscience Initiative 5: integrated multidisciplinary studies of unconformity-related uranium deposits from the Patterson Lake corridor, northern Saskatchewan. Geological Survey of Canada, Bulletin 615, 37 pp
- Puzyrev, V. (2019). Deep learning electromagnetic inversion with convolutional neural networks. Geophysical journal international, 218(2), 817–832
- Puzyrev, V., & Swidinsky, A. (2021). Inversion of 1D frequency- and time-domain electromagnetic data with convolutional neural networks. Computers & geosciences, 149, 104681
- Ronneberger, O., Fischer, P., & Brox, T. (2015). U-net: convolutional networks for biomedical image segmentation. Medical Image Computing and Computer-Assisted Intervention (MICCAI) (9351 vol., pp. 234–241). LNCS: Springer
- Roy, A. G., Navab, N., & Wachinger, C. (2018). Concurrent Spatial and Channel ‘Squeeze & Excitation’ in Fully Convolutional Networks. In A. Frangi, J. Schnabel, C. Davatzikos, C. Alberola-López, & G. Fichtinger (Eds.), Medical Image Computing and Computer Assisted Intervention – MICCAI 2018. MICCAI 2018 (11070 vol.). Cham: Springer. Lecture Notes in Computer Sciencehttps://doi.org/10.1007/978-3-030-00928-1_48
- Russell, B. (2019). Machine learning and geophysical inversion; a numerical study. Leading edge (Tulsa, Okla.), 38(7), 512–519
- Schetselaar, E. M., & Snyder, D. B. (2017). National database of MOHO depth estimates from seismic refraction and teleseismic surveys. Geological Survey of Canada, Open File, 8243, 14. https://doi.org/10.4095/305396
- Smirnov, M. Y., & Pedersen, L. B. (2009). Magnetotelluric measurements across the Sorgenfrei-Tornquist Zone in southern Sweden and Denmark. Geophysical Journal International, 176(2), 443–456
- Tschirhart, V., Pehrsson, S., Card, C., Potter, E. G., Powell, J., & Pană, D. (2021). Interpretation of buried basement in the southwestern Athabasca Basin, Canada, from integrated geophysical and geological datasets. Geochemistry: Exploration, Environment, Analysis, 21(1), geochem2019–geochem2061
- Tschirhart, V., Potter, E., Powell, J., Roots, E., & Craven, J. (2022). Deep geological controls on formation of the highest-grade uranium deposits in the world: magnetotelluric imaging of unconformity-related systems from the Athabasca Basin, Canada. Geophysical Research Letters, in press
- Ward, S. H., & Hohmann, G. W. (1987). Electromagnetic Theory for Geophysical Applications, Electromagnetic Methods in Applied Geophysics. Volume 1, pp. 201
- Zhdanov, M. S. (2002). Geophysical Inverse Theory and Regularization Problems (p. 214). Elsevier Science
- Zhdanov, M. S., Wan, L., Gribenko, A., Cuma, M., Key, K., & Constable, S. (2011). Large-scale 3D inversion of marine magnetotelluric data: Case study from the Gemini prospect, Gulf of Mexico. Geophysics, 76(1), F77–F87