Two‐dimensional regularized inversion of AMT data based on rotation invariant of Central impedance tensor

Considering the uncertainty of the electrical axis for two‐dimensional audo‐magnetotelluric (AMT) data processing, an AMT inversion method with the Central impedance tensor was presented. First, we present a calculation expression of the Central impedance tensor in AMT, which can be considered as the arithmetic mean of TE‐polarization mode and TM‐polarization mode in the two‐dimensional geo‐electrical model. Second, a least‐squares iterative inversion algorithm is established, based on a smoothness‐constrained model, and an improved L‐curve method is adopted to determine the best regularization parameters. We then test the above inversion method with synthetic data and field data. The test results show that this two‐dimensional AMT inversion scheme for the responses of Central impedance is effective and can reconstruct reasonable two‐dimensional subsurface resistivity structures. We conclude that the Central impedance tensor is a useful tool for two‐dimensional inversion of AMT data.


Introduction
Audio-magnetotelluric (AMT) sounding is used to detect how the electrical resistivity of the earth varies with depth. It is similar to standard magnetotelluric measurement in that it utilizes a passive natural source, but in AMT the frequency band is limited to the audio range, generally from 0.1 Hz to 8000 Hz. AMT has wide applications in mineral and ground water exploration as well as in environmental tasks such as fault-zone mapping, geothermal investigation, and engineering detection (Barcelona et al., 2013;Huang et al., 2017;Share et al., 2014;Tarabees et al., 2017;Tuncer et al., 2006;Yamaguchi et al., 2010). The main application of the AMT method, the modeling of subsurface resistivity structures, requires inversion of the AMT responses.
Two-dimensional inversion methods and technology applicable to AMT responses have become increasingly sophisticated and mature, reaching the stage of routine practical application (deGroot-Hedlin and Constable, 1990;Nittinger and Becken, 2018;Rodi and Mackie, 2001;Smith and Booker, 1991;Pedersen and Engels, 2005). However, for complex geological structures, if the traditional two-dimensional inversion method is applied to the quantitative inversion interpretation of the measured data of AMT, three aspects should be considered (Muñíz et al., 2017). Problem 1: The polarization mode of the AMT field data, that is, the discrimination between the TE-polarization mode and the TM-polarization mode, must be accurately described. Problem 2: Under normal circumstances, strike direction of the two-dimensional AMT field data is a curve that varies with frequency, but sometimes the difference between high frequency and low frequency strike directions is close to 90°. Problem 3: In AMT exploration, there is no case in which the survey line and the geological structure line are completely orthogonal. Pedersen and Engels (2005) proposed the invariants of the impedance tensor, and figured out that it is one solution for the uncertainty of the strike direction for two-dimensional AMT data inversion. All above three questions are caused by the uncertainty of strike direction in two-dimensional magnetotelluric data processing. The fact that the rotational invariance of the Central impedance tensor can make the forward modelling data have the same results in different electrical axes allows us to simplify the inversion process. We employ this idea to do regularized inversion in two-dimensional AMT.
The main difference between the inversion method presented in this paper and those in previous works is that the responses of the Central impedance tensor are applied to two-dimensional AMT regularized inversion. First, we present the calculation expression of the Central impedance tensor in AMT, which avoids the uncertainty of the strike direction for two-dimensional inversion. Second, we develop a two-dimensional AMT regularized inversion algorithm based on smoothness-constraint least-squares technique, adopting an improved L-curve scheme to determine the regularization parameter. Finally, we test the inversion performance on synthetic and field data, which confirms that this inversion scheme can reconstruct reasonable two-dimensional subsurface resistivity structures.

Rotation Invariant of Central Impedance Tensor
AMT rotation tensors are those whose amplitudes remain constant after 360° mathematical rotation. The four-tensor elements in AMT, i.e. Z xx , Z xy , Z yx and Z yy describe ellipses of the same size and ellipticity, not only theoretically but verified by field observations that show the same property (Roy et al., 2004). In the case of 2D structures the AMT impedance tensor connecting horizontal components of the surface electromagnetic field at a particular position and frequency becomes particularly simple in the principal coordinate system, with the x-axis pointing in the strike direction: where subscripts TE and TM stand for transverse electric and transverse magnetic modes, respectively.
The Central impedance tensor is rotationally invariant, and the Central impedance for 2D model can be written as (Lilley, 1993) where r and i, respectively, stand for the real and imaginary components of the impedance tensor element Z.
In order to show the characteristics of the Central impedance tensor, the AMT responses of a two-dimensional geo-electrical model are computed by the finite difference method (Tong XZ et al., 2018). The model consists of a rectangular conductor with a resistivity of 10 Ω·m homogeneous half-space, and the sizes of the ρ Cen surrounding homogeneous space and the conductive block are 8 km × 4 km and 2 km × 2 km, respectively. Figure 1 shows the plots for the AMT responses of the TE-polarization, TM-polarization and the rotation invariant Central impedance tensor over a two-dimensional geo-electrical model. Whatever the discrepancy of separation between the TE-polarization and TM-polarization mode apparent resistivities, the rotationally invariant values are very close to each other at all frequencies ( Figure 1a). Meanwhile, 2D Central impedance data can be considered as the arithmetic mean of TE-polarization and TM-polarization mode data.

Regularized Inversion Method
The AMT two-dimensional forward problem can be generally represented in a discrete form as d where F is a forward modeling operator which is generally nonlinear, m is a model parameter vector, and is a model response (predicted data) vector.
The problem of inverting AMT data is ill-posed and non-linear, but can be solved with various regularization methods. The following parametric functional is minimized in the regularization (Sun JJ and Li YG, 2014;Tikhonov and Arsenin, 1987): where is a misfit functional, is a stabilizing functional and is known as regularization parameter (trade-off parameter), which controls the trade-off between these two contributions in a minimization process. The misfit and stabilizing functionals are written as d obs W d where is the observed data vector, C indicates the model parameter weighting matrix, and is a diagonal data weighting matrix whose elements are reciprocals of the standard deviations where is the standard deviation of the noise in the ith observation.
Many algorithms have been suggested to solve equation (5) (Van Beusekom et al., 2011;Candansayar, 2008;Mehanee and Zhdanov, 2002;Singh et al., 2017;Siripunvaraporn and Egbert, 2007;Yan P et al., 2017). In this paper, the smoothness-constrained leastsquares inversion method is adopted for solving the regularized inversion problems (Negi et al., 2013;Lee et al., 2009). Linearization of equation (4) and the least-squares inversion model can be updated by solving where , denotes the model updates to be obtained, and J is the sensitivity matrix or Jacobian matrix of forward modeling operator F, and C is a Laplacian (second-order) smoothness operator (Kelbert et al., 2008;Lee et al., 2009;Tong XZ et al., 2017). Solution of equation (8) yields the model perturbation for updating the model . The inversion then proceeds to the next iteration.

Determination of Regularization Parameter
Determination of the regularization parameter, which balances the minimization of the data misfit and model roughness, may be a critical procedure to achieve both resolution and stability (Farquharson and Oldenburg, 2004;Haber and Oldenburg, 2000;Roy, 2002). A large regularization parameter, in principle, would apply more constraint or regularization to the solution and give poor resolution of the model parameters. On the other hand, too small a regularization parameter may be harmful to the stability of the inversion. In our implementation, an improved L-curve scheme is adopted, in which the selection of the regularization parameter needs to be given in advance.
Through extensive numerical experiments and Newman's practice in three-dimensional massively parallel electromagnetic inversion (Newman and Alumbaugh, 1997), the minimum and maximum regularization parameters ( , ) can be obtained when the regularization parameters are selected as the minimum and the maximum row sums of the matrix product . Therefore, the and can be given as follows and Here is an indicated element of and M is the number of model parameters.
Therefore, the best regularization parameter can be computed using the formula (Farquharson and Oldenburg, 2004;Hansen, 1997) where , and . The improved L-curve for our inverse problem is shown in Figure 2, and the best regularization parameter is the maximizer of the curvature of the L-curve given by equation (11). Hz and 100 Hz. Synthetic data are generated by the finite difference method (Tong XZ et al., 2018), and are contaminated by a 2% random noise.· Figure 4 shows the Bostick inversion results for TE-polarization, TM-polarization and the rotationally invariant Central responses. The results can be reconstructed both in shape and in spatial anomaly centers in Figure 4. Because the two-dimensional geo-electrical model has different AMT responses in TE-polarization mode and TM-polarization mode, the two-dimensional semi-quantitative inversion results are different from those of the one-dimensional Bostick inversion. Figure 4c shows that the inversion result of model I, based on the rotationally invariant Central parameters, has good recovery and high precision, which can be seen from the arithmetic mean of the TE-polarization and TM-polarization inverted models.

Inversion Tests for Synthetic
The total number of inversion blocks, 32×20(=640), is generated by the model parameterization of the two-dimensional AMT inversion in the computational domain. For an initial model for numerical test, we set resistivity of 100 Ω·m in the homogeneous model. The inversion result is shown in Figure 5; the anomalous low resistivity can be clearly identified in the inverted section. Figure 5c shows clearly that the simple geo-electrical model can be effectively recovered using the smoothness-constrained inversion algorithm for AMT responses of the Central impedance.
A comparison between observed and predicted responses is shown in Figure 6. As can be seen from Figure 6, the apparent resistivity-frequency pseudo-section and impedance phase-frequency pseudo-sections show that the Central impedance responses of the inverted model are in good agreement with the Central impedance responses of the true model.

Model II
Model II is the same as the one tested in Sasaki (1989), as shown in Figure 7a. In this model, a conductive layer of 5 Ω·m is embedded in a host medium of 50 Ω·m, but this layer is disconnected or faulted at the range of 9000-11000 m with a gap of 2000 m in the horizontal coordinate. There also exist one low-resistivity (10 Ω·m) and one high-resistivity (100 Ω·m) anomalous bodies just below the surface. The stations are located along a profile from 0 m to 18000 m at 500 m intervals in the horizontal coordinate, resulting in 35 sites in total. There are a total of nine frequencies, 0.1, 0.22, 0.5, 1, 2.2, 5, 10, 22 Hz, and 50 Hz. Synthetic data are given in the Central impedance tensor of model II, which were also generated by the finite difference method (Tong XZ et al., 2018), and are contaminated by 2% random noise.
The model parameterization for AMT inversion has divided the model into 40×28 (=1120) inversion blocks in the computational domain. For an initial model in the numerical test, a homogeneous model of resistivity of 100 Ω·m was chosen; its inversion res-  (Hz)  35  40  45  50  55  60  65  70  75  80  85  90  95  100  105  110  115   41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  ult is shown in Figures 7b-7c. The conductive layer and a fault gap in the range of 7500-12000 m are recovered very well in the inversion result. From the above model study, it is confirmed that twodimensional AMT smoothness-constrained inversion of the Central impedance tensor can efficiently recover the complex subsurfaces (see Figure 7c).
Moreover, comparison between observed and predicted responses, as shown in Figure 8, shows very good agreement. Therefore, the inversion result for the synthetic model shows that a smoothness-constrained regularized inversion method based on the Central impedance is effective, and can reconstruct reasonable two-dimensional subsurface resistivity structures.

Application to Field Data
The field data set in this paper is one acquired around the area of the Huangshanling Pb-Zn deposit in Chizhou City, Anhui Province, China. The AMT survey was used to obtain an electrical conductivity image to improve current interpretations in the east-south part of the Huangshanling Pb-Zn mine. The profiling line has 2360 m in survey, and the total number of AMT stations is 59. A 5-channel AMT data acquisition system (MTU-5) from Phoenix Geophysics Ltd was used to record the raw data, and impedance was done with SSMT2000 processing software.
The inversion for the rotationally invariant Central responses takes nine iterations to reduce the misfit to 13.2%; the two-dimensional inversion result is shown in Figure 9b. The electric resistivity structure characteristics can be divided into two general domains: (1) the shallow layer on the west side of the electric resistivity structure is high resistance, and the east slope is buried deep; (2) the shallow east side is mainly the low resistivity anomaly, and the high resistivity anomaly is confined to a relatively small locality. Figure 9c shows a sketch of the geological section, which can be explained by the two-dimensional electric resistivity structure characteristics (see Figure 9b) and magnetic data (see Figure 9a). The strata are diagonally inclined to the southeast, with an inclination of about 40° and partial secondary anticlines. The Ordovician strata have been buried deeper in the east, and the deep rock masses are shallow on both sides and buried deep in the middle. Combined with the magnetic exploration data, two magnetic layers were inferred from the top and bottom of the main interface: a pyrrhotization zone rich in S 1 g, seated above the main interface, and a skarn zone below the main interface, of thickness 10 m.

Conclusions
(1) A quantitative inversion interpretation of AMT data was established using responses of the Central impedance tensor. The uncertainty of the strike direction for two-dimensional data processing was successfully avoided, demonstrating a new way to interpret AMT data quantitatively.
(2) A two-dimensional AMT regularized inversion method, based on smoothness-constraint least-squares technique, was shown. Meanwhile, an improved L-curve scheme was adopted for determination of the regularization parameter, which balances minimization of data misfit and model smoothness.
(3) Inversion results for synthetic data and field data show that this two-dimensional AMT inversion scheme for the responses of the Central impedance is effective for reconstructing two-dimensional subsurface resistivity structures. We note that there are obvious limitations in using only two degrees of freedom for the Central impedance, rather than the four degrees of freedom in the TE-polarization and TM-polarization modes.