The MT method is a passive-source geophysical technique used to investigate the resistivity distribution of the Earth’s subsurface by analyzing the natural electromagnetic fields at the surface (Tikhonov, 1950; Cagniard, 1953). These fields mainly originate from regional and global thunderstorm activity, as well as the interaction of solar wind with the Earth’s magnetosphere. As electromagnetic waves propagate through the subsurface, their penetration depth is affected by the diffusion of waves, which varies with frequency. Specifically, waves with lower frequencies demonstrate a greater capacity for deeper penetration into geological structures. Since MT sources are sufficiently distant and the Earth possesses a higher refractive index than air, electromagnetic waves are assumed to propagate as planar waves perpendicular to the Earth’s surface. However, the diffusion process through various subsurface structures can exhibit arbitrary polarization. The MT method typically measures the orthogonal components of the natural electric and magnetic fields at the Earth’s surface, and a complex impedance tensor Z is used to describe the mathematical relationship (Equation 1) between the electric and magnetic fields (Berdichevsky, 1976): E = Z H , (1) where E = [ E x , E y ] T and H = [ H x , H y ] T are the electric and magnetic fields, respectively.

The apparent resistivity and phase are defined as Equations 2, 3, respectively. ρ = 1 ω μ 0 Z , (2) ϕ = a r c t a n I m Z R e Z , (3) where R e and I m are real and imaginary components of Z , μ 0 is the magnetic permeability of free space, and ω = 2 π f . f is the frequency.

Building on the AE framework, we replace the decoder with the MT forward operator to create a customized ANN (Figure 1). The loss function Φ d (Equation 5) in PGAE, representing the difference between input features d obs and output predictions d pre , corresponds to the data misfit typically employed in traditional inversions. Moreover, a regularization term Φ m (Equation 6) is employed to ensure the smoothness of the predicted geoelectric model m . Φ = Φ d + λ Φ m , (4) Φ d = 1 2 N ∑ i = 1 N d i obs − d i pre 2 2 , (5) Φ m = 1 2 N ∑ i = 1 N ∇ m i 2 2 , (6) d i pre = F m i = F A N N w , d i obs , (7) where N is the batch size used in training, F is the forward operator, λ is the regularization parameter, and w is a set of network parameters.

Algorithm 1 outlines the workflow for implementing PGAE. The process begins by initializing the encoder parameters using a truncated normal distribution, which helps mitigate the exploding or vanishing gradient problems during training Li H. et al. (2020). After initialization, the training samples are shuffled randomly, and a batch of samples is selected and fed into the encoder to generate the corresponding geoelectric models. The MT forward operator is then applied to produce the predicted MT data. Subsequently, the loss function is evaluated using Equation 4. These steps constitute the forward propagation phase. Next, the gradient g of the loss function is computed with respect to the network parameters w , and the parameters are updated using a suitable optimization algorithm, such as the Adam optimizer introduced by Kingma and Ba (2017). These steps constitute the back propagation phase. If the value of the loss function, Φ , falls below a predefined threshold, ε , the network is saved. Otherwise, the training loop continues until convergence.

Algorithm 1

The PGAE training algorithm.

Input: Initialize ANN weights w , target loss ε , N e max epochs, k = 0 ;

Scaling laws provide a crucial theoretical framework for simulating the electromagnetic response of real-world resistivity structures within a laboratory setting (Ward, 1967; Nabighian, 1987; Wong et al., 2009). Building on the foundational work of Ward (1967), we utilize the coordinates ( x , y , z ) to represent any point in the real-world system and the coordinates ( x ′ , y ′ , z ′ ) to denote the corresponding point in the model-scale system. The two systems are related through a length scale factor, denoted as p , and the transformation between them is governed by Equation 8. x = p x ′ , y = p y ′ , z = p z ′ . (8)

Ignoring the displacement current, the MT field is described by the diffusive forms of Maxwell’s equations with a time-dependent component e iwt . For the real-world system, we have ∇ × E x , y , z , f = − i ω μ 0 H x , y , z , f , (9) ∇ × H x , y , z , f = σ x , y , z , f E x , y , z , f , (10) where σ is the electrical conductivity.

For the model-scale system, we have ∇ × E ′ x ′ , y ′ , z ′ , f ′ = − i ω ′ μ 0 H ′ x ′ , y ′ , z ′ , f ′ , (11) ∇ × H ′ x ′ , y ′ , z ′ , f ′ = σ ′ x ′ , y ′ , z ′ , f ′ E ′ x ′ , y ′ , z ′ , f ′ . (12)

Generally, we want the fields from the real-world system to be linearly transformed from the model-scale system, that is E x , y , z , f = a E ′ x ′ , y ′ , z ′ , f ′ , (13) H x , y , z , f = b H ′ x ′ , y ′ , z ′ , f ′ . (14) where a and b are scale factors for the electric and magnetic fields, respectively. Therefore, ∇ × E x , y , z , f = a p ∇ × E ′ x , y , z , f , (15) ∇ × H x , y , z , f = b p ∇ × H ′ x ′ , y ′ , z ′ , f ′ . (16)

Putting Equations 9–16 together, we have f = a p b f ′ , (17) σ = b p a σ ′ . (18)

According to Equations 17, 18, to simulate the MT response of real-world resistivity structures in a model-scale system, we should ensure that the model-scale system has the same induction number as the real-world system: σ f l 2 = σ ′ f ′ l ′ 2 , (19) where l is any characteristic dimension in the real-world system and l ′ is the corresponding dimension in the model-scale system.

Scaling laws enable us to model large-scale MT responses using a small-size model. While researchers generally employ these principles for physical modeling, they are also applicable to numerical modeling. When we have an ANN trained for audio-frequency MT tasks, we can leverage the trained network to invert long-period MT data by applying scaling laws, rather than training a new network from scratch. Additionally, discrepancies between the number of frequencies in the real-world data and those used in training present another challenge. A practical solution is to interpolate the frequency sampling of the real-world data to match that of the training set. The interpolation method will perform effectively since MT apparent resistivity curves are typically smooth.

To prepare the training dataset, we analytically compute the apparent resistivity curves for 1D layered Earth models consisting of 31 layers and 25 frequencies, ranging from 10,000 Hz–1 Hz. The thicknesses of the layers begin at 20 m and increase logarithmically, as described in Equation 20, with the final layer always being a homogeneous half-space. For the resistivity of each layer, we uniformly sample values within the range of 0.1–100,000 Ω ⋅ m in the logarithmic domain and subsequently smooth the geoelectric model using a Gaussian filter (Liu L. et al., 2023). A total of 120,000 samples are generated, with 100,000 allocated for training and 20,000 for evaluation. Some random smooth geoelectric models can be found in Figure 2. d z i = 20 + 1 0 0.115 × i − 1 i = 2,3 , … , 30 , (20) where i represents the number of layers and d z i denotes the thickness of the i t h layer.

The primary objective of this study is to enhance the generalization capabilities of the trained ANN rather than to design a high-performance ANN. Accordingly, we employ an ANN architecture that is as simple as possible. Specifically, we construct an ANN with a single hidden layer containing 500 neurons. The number of neurons in the input (output) layer corresponds to the MT data (geoelectric model) used in the training dataset. The regularization parameter λ in Equation 4 is set to 4.5, determined through trial and error. We code the whole PGAE framework using TensorFlow (an open-source Python library by Google) so that automatic differentiation (Liu L. et al., 2024) can be employed to calculate the gradient of the loss function.

We use the Adam optimizer to train the established network, as illustrated in Figure 3A, the data misfit decreases rapidly as the number of training epochs increases. A significant reduction in data misfit is observed at the end of the first epoch, indicating that the training of the network converges quickly. The gradual increase in model roughness shown in Figure 3B suggests that the well-trained PGAE possesses the capability to produce models with complex structures. After 20 epochs, the loss function, data misfit, and model roughness exhibit minimal changes, indicating the completion of the training stage. The network was trained using an NVIDIA Tesla K80 GPU with 12 GB of onboard memory, taking approximately 40 min.

To evaluate the performance of the trained PGAE, we simulate numerous synthetic 1D MT data as test samples, which share the same observation frequencies as the training samples. For comparative analysis, we also invert these synthetic data using the classical Occam’s inversion (Constable et al., 1987). To further assess the noise adaptability of the trained PGAE, we introduce 5% Gaussian noise into the test samples. As illustrated in Figures 4, 5, the well-trained PGAE produces models comparable to those generated by the Occam’s inversion, with both methods yielding a satisfactory data fit.

Although no noise was intentionally added to the training samples, the network demonstrates remarkable adaptability to noise. Given that noise is unavoidable in real-world data, we further evaluate the noise tolerance of the trained PGAE by testing it against three higher levels of Gaussian noise: 10%, 20%, and 30%. As illustrated in Figures 6, 7, the models predicted by the PGAE are comparable to those derived from the Occam’s inversion, indicating that the PGAE exhibits significant noise tolerance even at high levels of Gaussian noise.

We select 11 MT sites recorded in the Tarim Basin for our real-world examples (Figure 8). The phase tensor analysis (Caldwell et al., 2004) shows that the resistivity distribution in the shallow area closely resembles a 1D structure (Figure 9), suitable for 1D inversion. Data from these sites were recorded using the Phoenix MTU-V5 system over approximately 22 h during the summer of 2019 and subsequently processed with the robust MT time series processing code, EMTF (Egbert and Booker, 1986), to estimate the impedance and vertical transfer functions. We utilize 41 frequencies ranging from 515 Hz to 0.0439 Hz, which share a similar bandwidth with the frequencies of the training samples but employ distinct frequency values. This enables the application of the scaling laws described in Section 2.3.

As illustrated in Figure 10, we first transform the frequencies of the real-world data to align with the frequency range of the training data. Subsequently, we employ the linear interpolation algorithm to match the number of real-world data to the neurons of the input layer of the trained ANN. Once the trained ANN generates an inverted model, the final step is transforming it back to its original scale, ensuring consistency with the real-world system. In this scenario, we scale the real-world frequencies by 10000/515, resistivity by 10, and model size by 515 / 100000 , meeting the conditions specified in Equation 19. The models predicted by PGAE are comparable to those inverted by Occam’s method (Figure 11). Both methods yield excellent data fitting (Figure 12). Additionally, we present the pseudo-sections of the 1D models derived from both methods in Figure 13, revealing consistent 2D structures. The prediction stage of PGAE requires only 0.113 s for all 11 sites, significantly faster than the Occam’s inversion, which takes approximately 5 s per site.