Multiscale Constrained Inversion Method for Direct Current Resistivity Tomography Based on Haar Wavelet Transform

Smoothness constrained inversion dominates the field of geoelectrical imaging. However, the resulting smooth images are unrealistic. In the absence of a priori underground information, the actual shape of the target geological structures are not well described. To address these limitations, we propose a multiscale inversion scheme. The inversion uses the wavelet transform method to convert the model parameters into Wavelet parameters. The Wavelet parameters are used as target values for the inversion. On this basis, different smoothness values were applied to different regions by controlling the feature parameter weights. Finally, the resistivity model was developed from the Wavelet parameters using the inverse wavelet transform. This study mainly considers the features of different depths and scales. By reducing the weight of the deep Wavelet parameters in the model objective function, the role of the data fitting objective function is enhanced, thereby improving the imaging of deep targets. Similarly, by reducing the weight of small-scale Wavelet parameters in the model objective function, the accuracy of the small-scale Wavelet parameters is improved. The small-scale Wavelet parameters correspond to the target boundary (fine or local) in the space domain, and hence, the actual shape of the target geological structures is better described. Compared with the traditional smoothness constrained inversion, the new method has a stronger boundary description effect on the target body. Numerical simulations have verified that different weights can improve the inversion performance. The feasibility of the algorithm was verified by using the sandbox test.


I. INTRODUCTION
For geological interpretation purposes, inverse problems have been studied via linear and nonlinear methods. The linear method based on Tikhonov regularization (Tikhonov and Arsenin [1]) has been widely applied in the fields of landslide monitoring, hydrogeophysics, and environmental engineering (Loke et al. [15], Loperte et al. [7], Wilkinson et al. [2], Uhlemann et al. [3], and Liu et al. [4]). The inversion result using the linear method sometimes falls into a local optimal solution (nonuniqueness of the inverse problem). In contrast, The associate editor coordinating the review of this manuscript and approving it for publication was Yi Zhang . nonlinear methods can perform global searches without relying on the initial model. The problems of slow convergence and overfitting in traditional nonlinear methods are overcome by deep learning methods (Liu et al. [5] and Li et al. [6]). However, deep learning methods are rarely applied to field data due to lack of data (paired geological models and survey data).
In general, the inherent nonuniqueness of the resistivity inverse problem can be alleviated by adding regularization constraints, including smoothness constraints (Sasaki [8]), inequality constraints (Kim et al. [9]), structural constraints (Li et al. [10] and Bergmann et al. [11]), level-set method (Li et al. [12]), and image-guided VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ constraints (Zhou et al. [13]). Most constrained inversion methods rely on priori information. However, smoothness constrained inversion is one of the most popular methods that needs no priori information. This method provides a smooth image of the real geoelectric structures, which sometimes seems geologically unrealistic (Portniaguine and Zhdanov [14]). The shape of the target is usually not well described as the discontinuities in the real resistivity distribution are ignored and the effects of other factors such as potential-field and complex environmental noise are not considered.
In response to these problems, multiscale inversion methods have been developed in recent years. In seismic and electromagnetic techniques, a common method is to extract the high-and low-frequency data using wavelet transform and then perform an inversion on models of different scales (Bunks et al. [16], Fichtner et al. [17], and Li et al. [18]). The method of direct current (DC) resistivity tomography is generally formulated as an ill-posed inverse problem and aims to reconstruct the accurate resistivity distribution underground. The data are usually obtained by injecting low-frequency current into the ground and collecting potential information. The number of data points is usually smaller than the number of the model results to be obtained; thus, the nonuniqueness of the solution is inherent. However, the data of DC resistivity tomography consist of electric potential or apparent resistivity. Thus, the wavelet transform cannot be used on DC resistivity tomography data. Chiao and Kuo [19] proposed a model parameterization method using the Haar wavelet transform. Liu et al. [20] studied the multiscale inversion method using the Daubechies wavelet in frequency domain airborne electromagnetics. They verified that low-order wavelets (such as Haar wavelets) can enable the recovery of sharp discontinuities. A multiscale resistivity inversion based on the convolutional wavelet transform was proposed and experiments showed that the method can accurately locate and delineate the boundaries of geological targets (Pang et al. [21]). In these methods, inversion is transformed from the model domain to the feature domain by wavelet transform. These methods aim to produce an image with clear boundaries by applying uneven smoothing constraints. Studies on the applications of multiscale parameterization to resistivity inverse problems are limited.
The Haar wavelet is the simplest orthonormal wavelet with compact support. It can quickly extract features of different scales in the model (low computing requirement) and it is widely used in pattern recognition and image processing (Stankovicá and Falkowski [22]). However, the multiscale inversion of DC resistivity using the Haar wavelet may produce redundant structures, and the imaging effect in deep areas is poor. In this study, we investigate the multiscale resistivity inversion method where the inversion process is guided by changing the weightage; i.e., smoothness is given different weights at different depths and scales. This method allows for sharp contrasts between different zones, and small fake anomalous regions generated by strong oscillations of the wavelet basis function are suppressed. The remainder of this paper is organized as follows. Section II introduces the multiscale inversion method with the weight matrix. Section III discusses the weighting effect through numerical simulation. In Section IV, the performance of the new method is verified using synthetic models.

II. INVERSION METHODOLOGY A. HAAR WAVELET TRANSFORM
Wavelet transform can decompose signals into different scales and preserve a sense of spatial location (Ridsdill-Smith and Dentith [23]). This study uses the simplest Haar wavelet, which is represented by H(·). The calculation process of the resistivity model is shown in Fig. 1. First, the resistivity parameter matrix is expanded to a one-dimensional (1-D) vector. The detail coefficient is obtained by averaging the differences in adjacent parameter values. Similarly, approximate coefficients are obtained by averaging adjacent parameter values. This operation is then continuously applied to the approximation coefficients until only one approximation coefficient remains. Haar wavelet transform is expressed as follows: where m represents the 1-D vector andm represents all the detail coefficients and the last approximation coefficient. Note that the horizontal resolution of the inversion result is usually higher than the resolution in the depth direction. Since discontinuous smoothing constraints are imposed on highconfidence horizontal directions, the resistivity parameter expands to 1-D by row.

B. MULTISCALE INVERSION EQUATION
The DC resistivity multiscale inversion function was reported with the following objective function (Pang et al. [21]): The L1 norm is used in the model term to obtain a sharp resistivity distribution. We further consider the effect of smoothness at different depths and scales. The objective function is then expanded to the following: where W Z and W S are the weight matrices of depth and scale, respectively. The values of the weight matrix elements are proposed by Li and Oldenburg [24] in the following format: where S and Z represent the scale and depth of the corresponding parameter. S 0 , Z 0 , β1 and β2 are empirical parameters. It should be noted that detail coefficients on a large scale often correspond to model parameters at multiple depths. The above formula of the depth weighting matrix only corresponds to small-scale detail coefficients. The weight values corresponding to large-scale detail coefficients are all 1.0 to ensure the integrity of the matrix. The data fit term can be expressed as follows: where G is the forward modelling response, calculated using the finite element method.J represents the sensitivity matrix: This study uses the general form of the Lp norm (Ekblom [25]), which is where ε is a threshold value. In this work, ε = 0.1. After minimizing the objective function (3) and adding a damping term, the inversion equation is obtained.
where λ and µ are the regularization coefficients. In this work, µ = 3 × 10 −3 and λ = 5 × 10 −6 . The elements of matrix R are Using an accurate initial model makes multiscale inversion more efficient. If a uniform resistivity model is used as the initial model for multiscale inversion, the detail coefficients of the first iteration are all 0; i.e., the model terms do not affect the inversion. Therefore, we only use the traditional smoothness constrained inversion in the first iteration.

III. NUMERICAL EXPERIMENTS
To validate our inversion algorithm, we used three geoelectric models (Fig. 2, Fig. 4, and Fig. 6). The first model was to illustrate the effectiveness of weighting the model terms. The last two models were compared with the traditional smooth least-squares constraint.

A. WEIGHTED EFFECT
Since the depth resolution of the surface data is poor, when minimizing the model terms, the effect tends to concentrate near the surface. Therefore, we used a simple weighting matrix (Z 0 = 20, β2 = 0.5) to improve deep imaging. Moreover, the interest of Wavelet parameters at different scales was different. Large-scale Wavelet parameters typically focus on the background, whereas small-scale Wavelet parameters are more concerned with the internal smoothness of each region and thus they display the boundaries more clearly. However, the number of large-scale feature data points is relatively small, and sometimes there is greater noise in the background. Therefore, the weight of the large-scale data in the equation needs to be increased. Our solution was S 0 = 2 and β1 = 0.8.
As shown in Fig. 2, a two-dimensional (2-D) finite element geoelectric model (124 m × 32 m) was constructed using a quadrilateral grid (2 m × 2 m). The background resistivity was 100 ·m. The resistivity of the two rectangular low-resistivity objects was 20 ·m, and their sizes were 6 m × 10 m and 8 m × 14 m, respectively. Their top buried depths are 8m and 10m, respectively. Data were obtained from 64 electrodes on the surface (electrode spacing was 2 m). The complete dataset contained 4925 independent data points, including Schlumberger and Dipole-Dipole arrays.
The relative root mean square (RMS) data misfit error is calculated as follows: where N denotes the number of data points.
The imaging results using the scale and depth weighted constraints are shown in Fig. 3a. As a comparative test, Fig. 3b is the result without weighting. Fig. 3 shows that inversion results successfully characterize the shape of the rectangular body. However, due to the comprehensive weighting of the Wavelet parameters, we can better identify the lower boundary of the deep objects in Fig. 3a. RMS values are relatively small, indicating that the inversion process shows a good convergence.

B. SIMPLE MODEL
To verify the imaging effect of the high resistivity target, a simple model was designed. In this sample, we compared the imaging results of the new and traditional smoothness constrained methods.  As shown in Fig. 4, the resistivity of the two rectangular high-resistivity objects was 500 ·m, and their sizes were 6 m × 6 m and 8 m × 6 m, respectively. Their top buried depths are 10 m and 12 m, respectively. Other conditions for modeling, data acquisition, and inversion parameters were consistent with those described in section III A. Comparing Figs. 3 and 5, it can be seen that the imaging effect of the high-resistivity target body is not as good as that of the low-resistivity target body. Advantageously, both imaging methods can accurately locate the high resistivity target, and the RMS value is also less than 4%. There is a clear transition area around the target boundary in the image based on smoothness constraints (Fig. 5b). But the new method performs better in recognizing the shape of the target's boundary than the traditional smoothness constrained method. The imaging results using the new method were closer to the true shape.

C. COMPLICATED MODEL
To verify the imaging effect of complex targets, two target bodies with complex shapes were designed. The shapes of the objects are shown in Fig. 6. Other conditions for modeling, data acquisition, and inversion parameters were consistent with those described in section III A.    7 shows the imaging results. The center positions of the two targets were identified accurately by using both methods. The RMS difference was small, suggesting that the new method achieves better inversion convergence. The current density and electric field strength gradually decrease with increasing depth, resulting in low sensitivity of potential data to deep targets. As shown in Fig. 7b, imaging based on traditional smoothness constrained method does not correctly describe the boundary of the left target body. But the new method can get the general outline features of the target body.

IV. PHYSICAL MODEL TEST
Physical model tests were designed to detect and image low resistivity objects at a farm in Zhangqiu City, Shandong Province, China. As shown in Fig. 8, the volume of the physical model was set at 6.5 m × 1.0 m × 1.2 m (length × width × depth). The model was supported by wooden baffles on both sides and filled with standard sand grains (sizes ranging from 0.5 mm to 2.0 mm). Sand grains with low water content represent the background of high resistivity in the model. After the sand was compacted, the background resistivity ranges from 200 ·m to 300 ·m. All the targets were made of a mixture of clay, salt, and water and act as low resistivity bodies in the model. The resistivity value was controlled within the range of 20 ·m to 50 ·m through ratio adjustment. Two sets of experiments were carried out. The respective site environments and spatial locations of the targets are shown in Figs. 9 and 10. Other specific parameters such as shape and size are shown in Table 1. Multiple targets were placed in the model, and the surface ERT method was used to collect the data. Finally, the new method and the traditional smoothness constrained inversion method were used for processing.

A. WEIGHTED EFFECT
A total of 64 electrodes (copper nails) were installed in the centerline with a spacing of 0.1 m. Each on and off time was 0.3 s and 0.5 s, respectively, and the data collection time was 1 h. After reciprocal data noise evaluation, the dataset contained 1985 independent data points. The temperature on the day of the test was -3 • C to 5 • C. The basic grid size in the finite element model was 0.1 m × 0.1 m. All inversion parameters were consistent with those presented  in section III. The inversion time cost of the two methods is almost the same. The results of the two tests are shown in Fig. 11 and Fig. 12, respectively.

B. RESULTS AND ANALYSIS
There are many artifacts in both inversion images owing to the uneven spatial compaction, uneven sand moisture content, and the non-uniqueness of inversion. Besides, the imaging areas of the two targets are larger than the actual ones due to the inherent volume effect of the resistivity method and because the imaging profile is affected by lateral anomalies in three-dimensional (3-D) detection.  We first analyzed the results of Test 1, as shown in Fig. 10. For L-shaped objects, smoothness constraints produce dull and blurred images. However, the new method has uneven smoothness, and its imaging results are closer to reality than the smoothness constrained method. For rectangular objects, the imaging results are larger than the actual model, which means that the new method cannot eliminate the volume effect. Nevertheless, the distance between the position of the two targets and the actual placements is less than 10 cm. To analyze whether the physical model fits the data well, we compared the inversion results of forward modeling (new method) with the measured data. According to statistics, the ratio of the inversion results to the measured data within 5% is about 83.6%. It suggests that the new method can be effectively used for the real-world data, but it needs further research for improvement.
The results of Test 2 are shown in Fig. 11. The results of the two methods can roughly get the spatial position of the two target bodies. The left target is offset 10cm to the right. The imaging result of the traditional method has a smooth transition boundary, which is much larger than that of the actual model. But the results of the new method can more accurately describe the boundary. According to the statistical results of the data, the proportion of the data volume in which the error between the inversion results and the measured data is within 5% is 85%.

V. DISCUSSION AND CONCLUSION
This paper develops a multiscale resistivity inversion method with weighted constraints. The new method considers both depth and scale effects. The weighting factor is determined by hundreds of numerical simulations. Compared with the other algorithm, the advantage of this method is that the shape of the target object can be more clearly represented by controlling the utility of smoothing constraints in different scales. This advantage is especially obvious in the numerical simulation. In the physical model test, the performance of multi-scale inversion is slightly better than that of the smooth constraint. This may be caused by factors such as the complexity of the environment and errors in the test process. Therefore, it is necessary to continue to focus on the processing of actual data. This will improve the practicality and reliability of field detection.
The DC multiscale inversion method is specifically prescribed for spatial model parameters with multiple scales. Multigrid technology can be introduced in future research for further improvements (Pessel and Gibert [26]). The inversion process transitions from a coarse grid to a fine grid to achieve true multiscale inversion. This may further improve the imaging effect.
It is difficult to obtain a good boundary effect only by the inversion of the surface potential data. Increasing the effective information in the data may improve the imaging accuracy. By using other observation methods (e.g., seismic monitoring and radar observations) and techniques (e.g., drilling) can be developed to analyze geological conditions comprehensively based on engineering requirements. In short, a new priori information is added to the objective function of DC multiscale inversion, or in other words, a joint inversion of the multiscale inversion method and other geophysical methods is realized. It is expected to solve the problem of engineering detection efficiently.
It was the first time that the multiscale weighted constraint method was used in resistivity inversion. This method can be extended to generalized electrical methods such as induced polarization and complex resistivity inversion. This is a new research direction in electrical inversion, which is expected to gain more influence in practical application. The goal of our work is to improve the resolution and reliability of electrical exploration. The subject of future research could be involved in multigrid technology, the dynamic weighting of specific regions, and joint inversion methods. Those methods as mentioned above contain more prior information that can further improve the imaging resolution.