Abstract
Based on researches of Baiyao (1995) and Zhang and Xu (2001), this paper proposes an improved 2D MT Zohdy-Oldenburg direct inversion method, in the least-square sense, embodying the features of Zohdy’s ratio method and Oldenburg’s difference method, in the condition of rugged topography, with phase information. It bypasses large calculations of the Jacobian matrix and large sparse linear systems of equations and enables direct modifications and comparisons of the model parameters. According to the calculation and analysis of examples, it shows faster convergence and higher precision. In contrast with the conventional linear inversion, the calculation speed of this new method can be increased by more than 10 times.
1. Introduction
Calculating Jacobian matrix and solving large systems of linear equations are required in conventional magnetotelluric (referred to later as MT) linear or nonlinear optimized inversion method [1–4]. It is imaginable that the calculation is very complicated and the systems of equations are often ill conditioned. Therefore in this paper we think it is necessary to find a quick and effective inversion method to solve practical problems in engineering projects.
Zohdy’s [5] method for the inversion calculation of direct-current resistivity sounding method (1989) converts polar distance (AB/2) to equivalent depth and modifies the model resistivity with the ratio between forward curve and measured curve to eliminate errors caused by zero-order asymptotic inversion (ratio method). This is called the asymptotic iterative inversion. Oldenburg and Ellis [6] in 1991 also proposed an iterative inversion method for MT studies, which modifies the model resistivity and model depth with the difference, calculated from Bostick, of forward result data and measured data of the prediction model inversion (difference method). Zohdy’s method is to modify the model resistivity with ratio method, so the rate of iterative convergence is apparently slower than that of Oldenburg’s method. In Zohdy’s iterative process, sometimes the model modification may be too much that the fitting errors will diverge. To solve this problem, two papers authored, respectively, by Shizhe and Bin (1995) [7, 8] suggest that one correction factor should be added to Zohdy’s method to control the step size which modifies the model resistivity. This inversion method is named as the curve comparison. Baiyao (1995) [9] repeated his study on the methods used by Zohdy and Oldenburg and proposed a least-square-sense-based direct inversion method, which combines Zohdy’s ratio method with Oldenburg’s difference method by filling longitudinal conductivity in the model grid and applies the combination in the inversion of direct-current resistivity sounding method [10] and MT sounding method [11]. On the basis of researches of Shizhe and Bin (1995) [7, 8] and Weidelt’s MT impedance phase and apparent resistivity theory [12], Zhang and Xu [13, 14] added impedance phase information to Zohdy’s ratio method and thus improved the rate of convergence in its iterative process.
This shows that the MT direct inversion, bypassing calculating Jacobian matrix and solving large systems of linear equations, directly modifying and providing comparisons for the model parameters, can fully rival Bostick inversion. Therefore, to solve practical problems like slow speed of inversion and poor fitting precision found in engineering projects, this paper, based on previous studies, aims to improve the direct inversion method, namely, to propose an improved 2D MT Zohdy-Oldenburg direct inversion method with phase information.
2. Improved Zohdy-Oldenburg Direct Inversion in the Least-Square Sense
Suppose that there are parameters of a geoelectric prediction model, and they are . Their relationship with surface observed values, , is as follows: In this equation, means theoretical observed values obtained from the forward computation of model parameters, and means errors between observed values and theoretical values.
To obtain the result , that is to say, the errors between observed values and theoretical values can reach a certain degree of precision, multiple iterations and modifications of , parameters of the geoelectric model, are often needed before a group of optimized are produced. This is the general process of geophysical inversion.
During the th iterative calculation, for the prediction model , we expand (1) to Taylor series and neglect terms which contain higher than second power in , and then we get When there are observed data, in matrix form, the above equation can be put into where is the matrix, the Jacobian matrix (partial derivative matrix).
Establishing the objective function is often needed in conventional linear inversions: where is the matrix norm, and when is 2, it becomes the least-square method. To make the objective function tend to minimum and thus determine the amount of the model parameters’ modification , calculating the Jacobian matrix and solving systems of linear equations are needed. This kind of inversion is burdened with large calculation and often contains ill-conditioned systems of equations.
To avoid this problem, we introduce the Zohdy-Oldenburg direct inversion method, which can also be called the asymptotic iterative method. If we replace the in the Jacobian matrix with , the objective function in (4) becomes And when we use the second-order norm, this turns into a question about the least-square method, and its solution is where and and therefore the following equation can be derived:
To keep the algorithm from converging, a correction factor will be added to (7) to limit the step size of modification. Thus we can obtain the iterative formula of Zohdy-Oldenburg direct inversion method in the least-square sense.
In this case, when it comes to the th model parameter, the inversion formula for th iterative calculation is When , (8) becomes Oldenburg’s asymptotic iterative method, which modifies the model with the difference method.
When , (8) turns into If , (9) turns into Zohdy’s asymptotic iterative method, which modifies the model with the ratio method.
Hence we can conclude that the improved Zohdy-Oldenburg direct inversion has features of asymptotic iterative method of the two scholars, ensuring that the algorithm in the iterative process can achieve faster convergence and higher stability of data.
3. Two-Dimensional Improved Zohdy-Oldenburg Direct Inversion
3.1. Two-Dimensional Direct Inversion
In the 2D inverse modeling of MT, we adopt an improved Zohdy-Oldenburg direct inversion method. The geoelectric prediction model’s parameter, , is divided into two parts, the resistivity and the depth . So (9) can be changed to the MT 2D improved Zohdy-Oldenburg direct inversion iterative formula: Equation (10) is the iterative formula for the improved Zohdy-Oldenburg direct inversion and gives resistivity inversion of the model. In this formula, stands for the true resistivity (or called intrinsic resistivity) of th model at th sounding point, of th layer; the measured apparent resistivity at th sounding point, th frequency point; the theoretical apparent resistivity of th model at th sounding point, th frequency point; in this paper the correction factor .
Equation (11) is the iterative formula for the improved Oldenburg’s iterative inversion and gives model depth inversion. In this formula, represents the depth of th model at th sounding point, of th layer; the apparent depth at th sounding point, of th layer (calculated from Bostick’s inversion formula); the theoretical apparent depth of th model at th sounding point, of th layer (calculated from Bostick’s inversion formula); in this paper the correction factor .
3.2. Two-Dimensional Direct Inversion with Phase Information
Zhang and Xu [14] and Weidelt [12], 1972, derived the approximate relationship between MT apparent resistivity and impedance phase: In this relationship, represents the frequency, the impedance phase, and the apparent resistivity.
If , then Bostick inversion formula is already known, so we get where stands for the true resistivity in Bostick inversion formula.
When we plug (13) into (14), we obtain the inversion formula of impedance phase: When small changes of impedance phase bring changes to resistivity, namely, both sides of (15) take a derivative with respect to , the result is Ignoring the local minimum of differential, we find that (16) turns into If we plug (17) into (10), the iterative formula of Zohdy’s direct inversion, we get Equation (18) is the improved iterative formula with phase information for Zohdy’s direct inversion. In this equation, stands for the measured impedance phase at th sounding point, th frequency point; the theoretical impedance phase of th model at th sounding point, th frequency point.
The model depth inversion still uses (11), the improved iterative formula for Oldenburg direct inversion. If we put (18) and (11) in the simultaneous form, then we will obtain the 2D MT improved Zohdy-Oldenburg direct inversion with phase information.
3.3. Calculation Method of Direct Inversion
The improved Zohdy-Oldenburg direct inversion is based on the direct inversion of Zohdy and Oldenburg, and it constructs the initial model of inversion with Bostick inversion formula and combines the formulas of Zohdy and Oldenburg simultaneously while it inverses the model resistivity and model depth. This method, featuring faster speed, higher precision, and faster convergence, can tackle a series of problems caused by solving large sparse Jacobian matrix.
Process of the improved Zohdy-Oldenburg direct inversion is as follows.
Step 1. Use Bostick inversion formula to calculate Bostick true resistivity and layer depth of each electric horizon, and fill 2D initial prediction model in the mesh model and perform a forward computation is as follows.
Step 2. Compare the measured apparent resistivity with the theoretical apparent resistivity obtained from the forward computation to check whether the precision can meet requirements. If the answer is yes, then it is certain that the prediction model is the result of the inversion.
Step 3. If not, modify the model parameters with Zohdy-Oldenburg iterative formula and obtain a new 2D prediction model. Conduct Bostick inversion, and get true resistivity and layer depth of each electric horizon. Refill the new 2D prediction model in the mesh model and again perform a 2D forward computation.
Step 4. Return to Step 2.
The fitting error between the measured apparent resistivity with the theoretical apparent resistivity will adopt mean square error; that is to say, where stands for the th prediction model’s mean square error between the measured apparent resistivity and the theoretical apparent resistivity, the number of measuring points, and the number of frequencies.
4. Analysis of Examples
4.1. One-Dimensional Layered MT Model
Figure 1 is the H-type MT cross-sectional model of horizontally layered earth (3 layers), which is nonproportional and for illustration only. The resistivity and the layer thickness of the first layer are, respectively, 100 Ω·m and 500 m; for the second layer, they are 10 Ω·m and 500 m; the third layer’s resistivity is 1000 Ω·m. The computation domain is 0~4 km, the dot pitch is 100 m, and it totaled 41 measuring points. The range of frequency is 103~10−3 Hz, it totaled 61 frequency points, and the sampling interval is base-10 logarithms.
Figure 2(a) is the inversion result of apparent resistivity, in which the horizontal axis represents depth (depth, unit: m) while the vertical axis represents resistivity (, unit: Ω·m); Figure 2(b) is the curve graph of inversion fitting error, in which the horizontal axis represents iteration times and the vertical axis represents fitting error.
(a) Inversion result of apparent resistivity
(b) Inversion fitting error
As can be seen from the inversion result of TE mode and TM mode, the anomalies of H-type stratum are revealed clearly. In the TE mode inversion result, the iteration times are 9, fitting errors percentage is 3.9%, and it takes 623 s; as for the TM mode inversion result, the iteration times are 9, fitting errors percentage is 6.2%, and it takes 410 s. The maximum relative error between these 2 modes is 5.3% and the mean square error is less than 2%.
4.2. Rugged Topography 2D Structural Model
Figure 3 is the rugged topography 2D structural model. The largest drop in the rugged topography is 500 m, which contains faults and high-resistivity and low-resistivity layers; the measuring points cover 0~20 km, the dot pitch is 200 m, and it totaled 101 measuring points; the frequency is within the common frequency spectrum of a V8 electrical prospecting apparatus, 102.5~10−2.96 Hz, the sampling interval is base-10 logarithms, and it totaled 38 frequency points.
Figure 4 shows the cross-sectional view of inversion results, in which the horizontal axis represents measuring point location (, unit: km) and the vertical axis depth (Depth, unit: km); Figures 4(a) and 4(b) are apparent resistivity inversion results of TE mode and TM mode, respectively, while Figures 4(c) and 4(d) are apparent resistivity inversion results of TE mode and TM mode with phase. Figure 5 is the fitting error curve of inversion iteration, in which the horizontal axis represents number of iterations or iteration times, while the vertical axis represents fitting errors. Figure 5(a) is the comparison curve of fitting error of TE mode inversion iteration with and without phase while Figure 5(b) is the comparison curve of fitting error of TM mode inversion iteration with and without phase.
(a) Cross-sectional view of TE mode inversion result
(b) Cross-sectional view of TM mode inversion result
(c) Cross-sectional view of TE mode inversion result (with phase)
(d) Cross-sectional view of TM mode inversion result (with phase)
(a) Comparison curves of fitting error of TE mode inversion iteration with and without phase
(b) Comparison curves of fitting error of TM mode inversion iteration with and without phase
From the above cross-sectional view of inversion results, it can be seen that the impact of rugged topography is eliminated, and inversions with topography of TE mode and TM mode are able to reflect locations and anomalies of faults and high-resistivity and low-resistivity layers. But it fails to reflect the second low-resistivity layer in the middle area of the figure, due to the equivalence phenomenon in geophysics.
In the TE mode inversion without phase, the iteration times are 15 and the fitting error is 6.2%. In the TE mode inversion with phase, the iteration times are 17 and the fitting error is 5.7%. In the TM mode inversion without phase, the iteration times are 20 and the fitting error is 5.3%. In the TM mode inversion with phase, the iteration times are 25 and the fitting error is 2.8%.
4.3. Model with Abnormal Body under Rugged Topography
Figure 6 shows the model with abnormal body under rugged topography, which is nonproportional and just for the convenience of easier understanding. In this figure, the drop is 60 m, the background resistivity is 100 Ω·m, and it contains a high-resistivity abnormal body of 500 Ω·m, which is 1600 m long and 150 m high, and it is 500 m away from its upper top surface to the ground surface. The measuring points cover 0~4 km, the dot pitch is 100 m, and it totaled 41 measuring points; the frequency is 104~10−2 Hz, the sampling interval is base-10 logarithms, and it totaled 61 frequency points.
Figure 7 is the cross section of the inversion result. It is a resistivity (or intrinsic resistivity) contour map (, unit: Ω·m), in which the horizontal axis stands for measuring points location (, unit: km) and the vertical axis depth (Depth, unit: km). Figures 7(a) and 7(b) are apparent resistivity inversion results of TE mode and TM mode, respectively, while Figures 7(c) and 7(d) are apparent resistivity inversion results of TE mode and TM mode with phase. Figure 8 is the fitting error curve of inversion iteration, in which the horizontal axis represents number of iteration times, while the vertical axis represents fitting errors. Figure 8(a) is the comparison curve of fitting error of TE mode inversion iteration with and without phase while Figure 8(b) is the comparison curve of fitting error of TM mode inversion iteration with and without phase.
(a) Cross-sectional view of TE mode inversion result
(b) Cross-sectional view of TM mode inversion result
(c) Cross-sectional view of TE mode inversion result (with phase)
(d) Cross-sectional view of TM mode inversion result (with phase)
(a) Comparison curves of fitting error of TE mode inversion iteration with and without phase
(b) Comparison curves of fitting error of TM mode inversion iteration with and without phase
The key of inversions with topography lies in the reliability of numerical modeling, which is influenced by topography, in forward calculation. As can be observed from the cross-sectional view of the inversion results, the impact of rugged topography is significantly reduced in the TE mode and the TM mode. And since the impact of topography on TE mode is relatively small, the inversion result of abnormal body of TE mode is better than that of TM mode. In the TE mode inversion without phase, the iteration times are 2 and the fitting error is 6.6%; in the TE mode inversion with phase, the iteration times are 22 and the fitting error is 0.25%. In the TM mode inversion without phase, the iteration times are 20 and the fitting error is 0.58%; in the TM mode inversion with phase, the iteration times are 25 and the fitting error is 4.3%.
5. Conclusions
In this paper, the proposed rugged topography 2D MT improved Zohdy-Oldenburg direct inversion is a asymptotic iteration method based on the least square, so the modification direction of this model always points to the descent direction of the objective function. This method embodies the features of both Zohdy’s asymptotic iteration method (ratio method) and Oldenburg’s asymptotic iteration method (difference method). The advantages of faster convergence and higher precision bypass large calculations of the Jacobian matrix and large sparse linear systems of equations and enable direct modifications and comparisons of the model parameters. In contrast with the conventional linear inversion, the calculation speed of the new method in this paper can be increased by more than 10 times.