Evaluation of Grounding Resistance and Inversion Method to Estimate Soil Electrical Grounding Parameters

Soil resistivity plays a key role in designing grounding systems for highvoltage transmission lines and substations. The objectives of this paper are to determine the best estimated value of the apparent resistivity or electrode grounding resistance of N-layer soil and to use a new inversion method to precisely determine earth parameters. The inversion of electrical sounding data does not yield a unique solution, and a single model to interpret the observations is sought. This paper presents a new inversion method to statistically estimate soil parameters from Schlumberger and Wenner measurements. To validate the method and test the inversion scheme, four soundings were selected: two theoretical and two in the field. The procedure was applied using test data and a satisfactory soil model was obtained.


I. INTRODUCTION
Ground electrode resistance plays an important role in designing power substations and transmission towers.Predicting this resistance is very complicated in some situations, especially when the soil is represented by an N-layer horizontal model with different resistivities.
Soil conductivity can vary across seasons and the depth of each layer is usually unknown, which makes determining the grounding resistance of an electrode particularly difficult.
The N-layer soil model is a realistic representation of the actual soil.Soil parameters can be estimated based on a set of soil resistivity measurements.
Several authors have contributed to solving the problem of soil parameter estimation to achieve the best fit between measured and computed resistivity values.
The Soil Measurements Interpretation Program (SOMIP) developed by Meliopoulos and Papalexopoulos [1] statistically estimates soil parameters based on resistivity measurements, but applied to field tests using four pin or three pin measurements.The program is designed to reject erroneous measurements resulting from instrument inaccuracies and changes in local soil resistivities.
Dawalibi and Barbeito [2] used imaging methods to compute the resistivity of multilayer soils by subdividing the soil into layers with thicknesses equal to a multiple of a base value.
The authors also used a method based on convolution and multiple filter theories to examine soil resistivity measurements in multilayer soil.
Takahashi and Kawase [3] developed a very compact formula to analyze changes in calculated apparent resistivity in multilayer soils, using a comparison of ρ-a curves.However, they did not provide an analytical solution to the problem of parameter estimation in such soils.
Del Alamo [4] compared eight different techniques to develop an optimum estimation of electrical grounding parameters for a two-layered soil model.
Lagacé et al. [5] used combined electrostatic images to evaluate N-layer soil resistivities and interpret sounding measurements from Wenner and Schlumberger electrode configurations.
Slaoui et al. [6] applied linear electric filter theory to derive the resistivity transform and proposed a method to estimate N-layer soil parameters from Schlumberger measurements using a ridge regression estimator.
Slaoui et al. [7] developed an innovative method to calculate the apparent resistivity of horizontally multilayered soils using an inversion method to determine the electrical grounding parameters of N-layered earth (resistivities and thicknesses).
Lagace and Vuong [8] proposed a method with a graphical user interface to estimate soil parameters of multilayered horizontal soil.The N-layer soil model illustrated in Figure 1 is a simplification of the actual soil.
This paper attempts to better estimate soil parameters from resistivity measurements.A statistical method to estimate the parameters of an N-layer soil model from Schlumberger measurements is presented.
This interpretation method works well for field data.Two field cases are presented.The two-sounding method was selected, not because the analysis was easy using the ridge regression method, but because the interpretation was very difficult using any method.

II. APPARENT RESISTIVITY EXPRESSIONS
Resistivity data is conventionally expressed in terms of horizontal layers having an isotropic resistivity.A model consisting of n layers would be parameterized by a set of layer resistivities ρ 1 , ρ 2 , ..., ρ n , and thicknesses h 1 , h 2 , ..., h n − 1 .The nth layer is the deepest and extends to infinite depth.Reducing Slaoui's [7] calculation of theoretical apparent resistivity curves, the following general recurrence expression is obtained: (j = n − 1, n − 2, ..., 2, 1). (2) where ρ j is the resistivity of the jth layer, h j is the thickness of the jth layer, ρ jaS is the Schlumberger apparent resistivity at the jth layer, n is the number of layers in the soil model, a is half the distance between current electrodes, and b is half the distance between voltage electrodes.The apparent resistivity ρ aS of the N-layer model is obtained at j = 1.

III. INVERSION METHOD WITH RIDGE REGRESSION
The problem of Schlumberger sounding over a plane-layered earth is nonlinear in the unknown parameters, namely, the resistivity and thickness of each layer.
As a starting point, we have a series i (i = 1, 2,..., N) of resistivity measurements ρ mS , and N is the number of measurements.
In order to estimate soil parameters in multilayer soils from resistivity measurements, the least squares fitting error between the measured values must be minimized and apparent values must be analytically computed for the same spacing.
The objective function to minimize is (5)

LEVENBERG-MARQUARDT METHOD
The Levenberg-Marquardt method is used to minimize multivariate functions.In N-layer soil, the values X = (ρ 1 , ρ 2 , ... ρ n, h 1 , h 2 , ... h n−1 ) T must be determined in order to minimize the distance function F. The F gradient is equal to zero at the absolute minimum.The Levenberg-Marquardt method consists of iteratively computing the F gradient by adjusting X from the Jacobien J such that it reduces the values of the gradient: From [4,9], the gradient of F can be evaluated at X +∆X: where J is a Jacobien matrix with partial derivatives of weighted differences with respect to each parameter.The apparent resistivity can be extended in the vicinity of X: (10) Equation ( 10) represents a linear system of N equations with M unknowns, where the number of unknowns is 2n − 1.Since this is a nonlinear problem, by substituting (10) in (6), the application of the least-squares inverse method yields (11) The method begins by calculating the eigenvalues of the matrix J T J. Small eigenvalues indicate a near-singular system, which is unstable in the presence of noisy data, and the average difference between estimated ∆X and true ∆X becomes very large.The ridge regression method [9,11] seeks to reduce this difference during the iteration process by damping the diagonal terms of (J T J).
The ridge regression estimate of ∆X rr is (12 where I is the unit identity matrix and .In the above equation, α = 0 and ∆X approaches the Taylor series direction.The eigenvalues of (J T J + αI) −1 are (λ i 2 + α), whereas λ i 2 are the eigenvalues of J T J. Increasing the size of all eigenvalues results in significant decreases in a) the mean of the squared length between ∆X and ∆X rr and b) the variance of the estimated solution.Therefore, in some cases, the solution ∆X rr is much closer to ∆X than the standard least-squares solution.The residual sum of squares for the ridge regression solution is given by (13) where ∆ρ* = ρ aS − ρ mS is using the values X* = X + ∆X rr .

IV. DATA ERROR AND WEIGHTING
Generally, field curves contain noise due to measuring errors, superficial inhomogeneities, telluric noise, and limited instrumentation precision.
When data are weighted, a relative degree of importance is assigned to each value.Such weighting may be used to remove inherent bias in the data or to bias the least squares fit so that it is more accurate in one area of the curve than another.
If there is a large numerical difference between the data values in different regions of the curve, an undesirable bias may be introduced into the final solution.The bias is such as to  = ( ) aS mS ρ cause the ridge regression estimator to be biased toward the large values while neglecting the smaller values, which may be as accurate and may contain some very important information.
In general, it is desirable to weight each data point according to the noise in that data point, and so as not to give it a false degree of importance due to its large or small value compared with other data points.
We assume that each data point has the same percentage standard deviation unless it is known or suspected that certain points are significantly noisier.It is further assumed, initially, that each point has a standard deviation equal to one percent of its measured value.The problem standard deviation σ is then adjusted to the estimated noise level of the survey.Most resistivity surveys yield data that is accurate within five per cent or less.
The solution consists of incorporating the weights within the residual sum of squares for the ridge regression, as given by ( 14) The residual error is now defined as (15) If the error in each measurement is independent of the error in other measurements, as is usually assumed, then W reduces to a diagonal matrix of data variance σ i 2 : (16) The choice of a weighting scheme does not appreciably affect minimization speed.However, it can drastically affect the position of the minimum in parameter space and the parameter statistics.
The weighted ridge regression estimator is given as (17)

V. PARAMETER STATISTICS
The next most important requirement is to obtain some idea about parameter uncertainty.The statistical parameters for our models are a) parameter standard errors and b) the parameter correlation coefficient.In addition to these statistical parameters, the parameter and data eigenvectors with the associated eigenvalues can yield great insight into the relations between individual model parameters and specific data.
Parameter standard errors and correlations are derived from the covariance matrix cov(X) evaluated at the minimum, as shown in [1,11]: When is greater than σ 2 , σ 2 has been underestimated.Although the value of σ 2 varies between surveys, it is assumed that the error at any data point is five percent or less.If is found to be smaller than that of the weighting, either the estimated variance of the observations have been overestimated or the curve calculated from the hypothesized model fits the noise in the data.
Thus, the residual variance can be used as an indicator of goodness-of-fit.The residual variance is independent of the linearity or nonlinearity of the problem with respect to the model parameters.
The parameter standard errors are defined by the square root of the diagonal term of cov(X) (e.g., [cov(X) ii ] 1/2 equals the standard error for numbered parameter i).
The correlation matrix is an indicator of the linear dependence between parameters.The correlation matrix elements are given in [2,11]. (20) If the value of [cor(X)] ij is near unity, then the parameters X i and X j are strongly correlated and nearly linearly dependent.For example, if i represents the thickness h and j, the resistivity ρ of a layer (i.e., [cor(X)] ij represents the correlation between the thickness and resistivity of a layer), only the ratio h/ρ is well determined by the data if [cor(X)] ij = 1.This is true for layers that are highly conductive relative to their surroundings.If [cor(X)] ij = −1, only the product hρ is well determined, as is the case for relatively resistive layers.This is the familiar equivalence problem, as discussed by Sunde [10].
The diagonal elements of the covariance matrix are the variance terms for each parameter.If the correlations are small, then the standard deviation is a good measure of the uncertainty of each parameter.If two parameters are strongly correlated, then the standard deviation given by the square roots of the diagonal terms of (18) will be larger than the actual uncertainties.
The eigenvectors and their associated eigenvalues are also very useful in defining the relations between soil parameters and their effects on data generated from a particular model.
The generalized inverse of J is defined in terms of these eigenvectors and eigenvalues, as described by Sunde [10]. (21) The matrix C consists of r (r is a rank of J) eigenvectors C i of length N associated with the columns (data) of J. B comprises the r eigenvectors B i of length M associated with the rows (parameters) of J.The matrix A −1 is the inverse of the diagonal matrix comprised of the eigenvalues of J.

VI. EXAMPLES OF DATA INTERPRETATION
In this section, four Schlumberger sounding curves and their associated models illustrate the ridge regression method and the estimation method for all parameter statistics.
In Figure 2, the first column shows the eigenvectors (columns of B), the second column shows the eigenvalues (diagonal elements of A), and the third column shows the data eigenvectors (columns of C), where the abscissas represents the number of measurements.The correlation coefficients and the model parameter values with estimated standard errors, calculated using data errors as weights in equation ( 14), are also shown.
The solid curve in Figure 3 is the final fit to the data points, which appears to be very good because there is little data noise.The final estimated model is close to the original model.
The first three eigenvalues have high magnitude.The linear combination of parameters represented by the first three eigenvectors has the greatest effect on the sounding curve.In these three eigenvectors, the elements ρ 1 and h 1 have opposite signs, whereas the elements ρ 2 and h 2 have the same sign.This indicates that if ρ 2 and h 2 are both changed in the same direction, the effect on the sounding curve (Figure 3) will be larger than the effect of similar   changes on other parameters.In addition, if ρ 1 increases while h 1 decreases, or vice versa, the sounding curve will change accordingly.The eigenvector associated with λ 4 = 1.17 indicates that increasing or decreasing ρ 1 and h 1 together will have little effect on the sounding curve (Figure 3).In other words, the ratio h 1 /ρ 1 is the combination of these parameters that affects the sounding curve the most (note that the correlation coefficient between ρ 1 and h 1 is +0.94).The eigenvector associated with λ 5 = 0.043 indicates that increasing ρ 2 while decreasing h 2 , or vice versa, has little effect on the sounding curve.Again, this is also indicated by considering the correlation coefficient between ρ 2 and h 2 , which is −0.91.
Considering the parameter eigenvector and parameter correlation coefficients of this model, it can be seen that the products ρ 2 h 2 are better determined by, and have greater effect on, this data than either ρ 2 or h 2 separately.
From the analysis of the eigenvectors, we see a wide range of resistivities and thicknesses for a very resistive layer that will yield nearly identical sounding curves as long as the thickness-resistivity ratio changes only slightly.In addition, for a very conductive layer, there is large range of resistivities and thicknesses that yield nearly the same sounding curve as long as the resistivity-thickness product is relatively constant.Thus, the parameters of highly resistive or highly conductive layers are difficult to determine.The correlation matrix must be considered in any interpretation of the standard deviations given by equation ( 18).
For the second example, we consider four-layer earth model #2: X(150, 700, 15, 200, 3, 20, 40). Figure 4 shows the theoretical data points for this model.The solid curve in Figure 4 is the final fit to the data points, and the final estimated model is close to the original model.
Figure 5 illustrates the parameter eigenvectors, eigenvalues, data eigenvectors, parameter correlation coefficients, and model parameter values with estimated standard errors.
In this example, we note five large eigenvalues and two small eigenvalues.Due to the several high eigenvalues, the inverse interpretation for this model is not very sensitive to noise, but it is very sensitive to model parameter changes.The eigenvectors may  be interpreted as for the three-layer problem.The eigenvectors associated with the small eigenvalues indicate the expected combination of parameters to which the sounding curve is insensitive, and the eigenvectors associated with the large eigenvalues indicate the combinations of parameters to which the problem is highly sensitive.For example, the eigenvector of λ 6 = 0.38 indicates that if the ρ 2 value increases while h 2 decreases, or vice versa, there will be almost no effect on the sounding curve.

B. FIELD EXAMPLE
The field data were taken from the Bryson site by Hydro-Québec.Consider the interpretation of the three-layer earth model #3. Figure 6 contains the data for the Schlumberger sounding and its computed sounding curve.The curve fits the data quite well, and the estimated noise level based on this fit is approximately 4 percent of the value at several data points.Since this is considered the approximate accuracy of the data, a closer fit is not justified.Table 1 illustrates the parameter correlation coefficients, and the model parameter value with estimated standard errors.All calculated parameters may be interpreted as for the three-layer problem in model #1.The data do not appear to provide much information about the third layer: we expected a much higher standard deviation for resistivity ρ 3 .Another interesting observation is the strong correlation between parameter ρ 3 and the thickness of the second layer.The final example, an interpretation of resistivity sounding measurements, was analyzed for a four-layer case using the Schlumberger configuration in Table 2.
Figure 7 illustrates the measured and calculated apparent resistivities, showing better correlation.
Table 3 illustrates the parameter eigenvectors, eigenvalues, and data eigenvectors.Table 4 illustrates the parameter correlation coefficients and the model parameter value with estimated standard errors.
All calculated parameters may be interpreted as for the four-layer problem in model #2.The data do not appear to provide much information about the fourth layer: we expected a higher standard deviation for the resistivity ρ 4. Another interesting observation is the strong correlation between parameter ρ 3 and the thickness of the second and third layers.

VII. ELECTRODE RESISTANCE EXPRESSIONS
Using an N-layer earth model, as shown in Figure 8, the electrode resistance can be approximated using Wenner's arrangement.The general recurrence expression for Wenner apparent resistivity is where ρ j is the resistivity of the jth layer, h j is the thickness of the jth layer, ρ jaW is the Wenner apparent resistivity at jth layer, n is the number of layers in the soil model, and C is the distance between any two electrodes.The apparent resistivity ρ aw of the N-layer model is obtained at j = 1: where R e is the ground electrode resistance Ω, and r is the electrode radius m.The height of the top layer (i.e., h 1 ) varies.When h 1 increases, the approximation of the electrode resistance becomes more accurate.When h 1 approaches infinity (homogeneous soil), the expression reduces to ρ ρ This is the expression for the resistance of hemisphere buried in a homogeneous soil.The computed ground electrode resistance and power loss are tabulated in Table 5 for varying values of h 1 .

VIII. CONCLUSION
The ridge regression estimator is a powerful method for interpreting resistivity soundings over plane-layered earth structures.Accordingly, it is possible to find a model that fits the data, to indicate the accuracy of the fit relative to the noise level in the data, and to predict the accuracy with which each parameter is estimated.The inversion method is based on a statistical estimation of the parameters of an N-layer soil model from Schlumberger measurements.This method provides the best estimate of soil parameters.
Finally, computation of the ground electrode resistance of hemisphere buried in N-layer soil is easier and faster using this method.

IX. APPENDIX
In order to obtain a recursive relationship for the partial derivatives, we differentiate Eq. ( 1): Using these expressions for j = 2, 3, ..., n, recursive application gives the partial derivatives of r aS = r naS Introducing G i (a, b, d j , i) = G i , we have:

Figure 1 N
Figure 1 N-layer earth model with Schlumberger configuration.

Figure 2
Figure 2 Parameters and data eigenvectors with associated eigenvalues, parameter correlations, and best fit model parameters for minimization with standard error weighted data.

Figure 3
Figure 3 Theoretical data and best fit curve for three-layer resistivity model #1 with Schlumberger configuration.

Figure 4
Figure 4 Theoretical data and best fit curve for a four-layer resistivity model #2 with Schlumberger configuration.

Figure 6
Figure 6 Field data from the Bryson site (Hydro-Québec) and interpreted best fit curve for three-layer resistivity model #3.

Figure 8
Figure 8  Hemisphere buried in N-Layer soil.

Table 1
Parameter correlations and best fit model parameters for minimization with standard error weighted data

Table 2
Numerical values of the 18 Schlumberger array measurements Figure 7 Field data and interpreted best fit curve for four-layer resistivity model #4.

Table 3
Parameter eigenvectors, eigenvalues, and data eigenvectors

Table 4
Parameter correlations and best fit model parameters for minimization with standard error weighted data

Table 5
Electrode resistance and power loss vs. height