A non-linear truncated SVD variance and resolution analysis of two-dimensional magnetotelluric models

SUMMARY A novel approach to assess variance and resolution properties of 2-D models of electrical resistivity derived from magnetotelluric measurements is presented. Based on a truncated singular value decomposition (TSVD) scheme on a local subspace, it partly takes the non-linearity of the inverse problem into account. The TSVD resolution and variance analysis is performed on a single cell at a time. A variance threshold is selected and the resulting model resolution is determined. As an improvement over existing schemes, non-linear semi-axes are introduced to describe the non-linear conﬁdence surface in the directions of the model eigenvectors and they replace the inverse singular values entering into the standard expression of model variances. The model variance of the cell considered is estimated from the sum of squares of the non-linear semi-axes up to the given variance threshold. This, in turn, gives the truncation level of the TSVD and the row of the model resolution matrix belonging to the considered cell can be computed from the model eigenvectors of the TSVD. The information contained in the resolution matrix is condensed to easily comprehensible measures like the centre of resolution and horizontal and vertical resolution lengths. The validity of our non-linear model variance and resolution estimates is tested with a most-squares technique which gives an improved estimate of model variability. A synthetic model with a conductive block in a homogenous half-space is analysed. TSVD analyses for model cells on the upper edge of the block and outside the block illustrate how the truncation process works. Typically, the linear and non-linear semi-axes are almost equal up to a certain singular value number, after which the non-linear semi-axes increase much less than the linear semi-axes. This important result indicates that the resolution of 2-D magnetotelluric models is signiﬁcantly better than previously suggested by linear schemes for the computation of model variance and resolution properties. A ﬁeld example from the Skediga area (Sweden) shows that the electrical resistivity distribution of sand and gravel formations which are only laterally bounded by conductive clay lenses is relatively well resolved whereas there is little resolution for the transition between the sand and gravel layer and the basement under a clay cover.


I N T RO D U C T I O N
An inverse problem can be split into several parts (Parker 1977): (1) investigating the existence and uniqueness of a solution, (2) constructing a solution, (3) analysing the stability of the problem and (4) assessing the significance of a solution.
Here, we describe a new method for evaluating the significance of a solution with a model resolution and variance analysis in the spirit of Backus & Gilbert (1968).
In order to compute model resolution and covariance properties, one generally needs to know both the sensitivity matrix of the considered model and the corresponding generalized inverse-preferably in the form of a singular value decomposition. The availability of the sensitivity matrix and its generalized inverse depends strongly on the inverse scheme that is used.
Several inverse methods have been devised for constructing a solution for non-unique and ill-posed, that is, unstable problems.  and later Lines & Treitel (1984) describe the singular value truncation technique and the Marquardt method as adequate means for tackling ill-posed inverse problems. More recently, Pedersen (2004) provided a 1-D inverse scheme for magnetotelluric data which determines an optimal truncation level for a singular value decomposition based on the total mean-square error of the estimated model. Due to the important statistical information that is given by the singular value decomposition, model variance and resolution information is readily available from these inverse schemes Lines & Treitel 1984;Pedersen 2004). Wiggins (1972) proposes to construct the generalized inverse for the solution of a linear inverse problem such that the model variances stay below a given threshold. As a consequence of the intimidating computing times of the singular value decomposition as the number of data and model parameters increases, these schemes are not widely used in 2-D or 3-D inverse computations. In both the original variant of the 2-D Occam scheme (de Groot-Hedlin & Constable 1990) and its data subspace formulation (Siripunvaraporn & Egbert 2000) the sensitivity matrix is computed and can be used directly to appraise model resolution and covariance. Mackie & Madden (1993) propose the conjugate gradient relaxation method which avoids explicit construction of the sensitivity matrix and only requires knowledge about the effect of the sensitivity matrix within matrix-vector products. Hence, the use of the conjugate gradient method requires other techniques of model resolution and covariance appraisal (see below).
Despite the fact that model resolution and covariance can be derived relatively easily from inverse schemes that compute the sensitivity matrix and its generalized inverse, this way of interpretation was not commonly taken in the past. This is partly due to the fact that the linear approximation to the inverse problem does not account for all models that fit the data to the same tolerance. Hence, most approaches taken to assess the resolution and covariance properties were restricted to derive a maximal depth of investigation pertaining to the considered data set from either linear or so-called non-linear studies or a combination of both. In this context, non-linear studies involve models that differ significantly from a preferred model obtained by a linearized inverse process. For the simple case of the 1-D inversion of LOTEM data, Kalscheuer et al. (2006) utilized different smoothness constraints of the Occam scheme to derive a maximal depth of exploration from non-linearly related models. Oldenburg & Li (1999) advocate the use of non-linear studies to compute a depth of investigation (DOI) index for DC resistivity and induced polarization data sets. After computing a preferred model m * that minimizes the model objective function subject to the data constraints, they modify the model objective function to obtain other models that fit the data to the same degree and are close but not linearly close to m * . This procedure is supposed to guarantee a more complete exploration of model space than granted by a linear analysis, like a Backus-Gilbert analysis (Parker 1977), which only considers models linearly close to m * . Depending on the model objective function chosen, the DOI index is either a cellwise normalized model difference or derived from a cross-correlation. Model cells above a certain DOI isoline are not considered to be resolved by the data. Schwalenberg et al. (2002) recommend a combination of linear and non-linear sensitivity studies. In the linear approach, columnwise sums of the sensitivity matrix normalized by the corresponding data errors are computed, that is, the sensitivities of all data points for one model cell enter into the sum. Inspired by Oldenburg & Li (1999), an isoline of a sensitivity deemed appropriate is used to exclude cells with lower normalized sensitivities from the interpretation. Because of the shortcoming of the linearized inverse process to find models that differ significantly from the estimated model, Schwalenberg et al. (2002) propose systematic non-linear sensitivity studies based on forward modelling as well as the inclusion of a priori models in forward modelling and inverse computations. Since forward modelling does not account for the interdependence between model parameters, the use of a priori model information in inverse computations is assumed to be more appropriate for deriving depth of investigation information (Schwalenberg et al. 2002).
A more rigorous resolution analysis applied to resistivity tomography data, that is not limited to estimating a depth of investigation, is proposed by Friedel (2003). He presents a single-step linearized inverse routine based on a truncated singular value decomposition. Because of the low-contrast approximation used, model resolution, model covariance and data resolution estimates are computed a priori, that is, from the SVD of the sensitivity matrix for the starting model. Friedel (2003) proposes to calculate a posteriori estimates if a non-linear inverse scheme with several iteration steps is used. Friedel (2003) also describes simple measures of resolution as a radius of resolution and a distortion flag (see below). From the model covariance matrix a percentage of image noise is extracted. Minkoff (1996) uses the conjugate gradient algorithm to solve seismic inverse problems. Since the sensitivity matrix and -more importantly -the normal matrix is never explicitly formed during the inversion, Minkoff (1996) forms an approximate model resolution matrix by estimating eigenvalues and eigenvectors of the normal matrix via the Lanczos method during the inverse procedure. In a similar fashion, Yao et al. (1999) derive both approximate model resolution and covariance matrices for an iterative LSQR solver of seismic tomography data.
For electromagnetic problems, Alumbaugh & Newman (2000) compute model resolution and covariance matrices for both direct and iterative inverse schemes. For both cases, however, resolution and covariance are estimated from the generalized inverse pertaining to the Tikhonov regularized inverse problem. Thus, the approach taken by Alumbaugh & Newman (2000) yields resolution and covariance appraisals which partly reflect the Tikhonov regularization of the inverse scheme.
In none of the above analyses (Minkoff 1996;Yao et al. 1999;Alumbaugh & Newman 2000;Friedel 2003) is the trade-off between model resolution and covariance (Menke 1989) or the influence of the non-linearity of the inverse problem on estimating model variance and resolution properties accounted for.
Here, we present a model variance and resolution analysis of resistivity models from a 2-D inverse scheme by means of a non-linear truncated SVD analysis in a local subspace. The analysis is carried out after the last iteration of the inverse process. No assumptions about smoothness constraints or other means of reducing the illposedness used during the inverse process are made, and the analysis is based on the SVD of the Jacobian as well as the semi-axes along the model eigenvectors to a given level of the non-linear confidence surface. Non-linear semi-axes were previously used by Johansen (1977) and Pedersen & Rasmussen (1989) to estimate non-linear parameter confidence limits of 1-D resistivity models, that were obtained from DC resistivity and magnetotelluric data, respectively. For the first time, we apply this technique to 2-D resistivity models and, more importantly, use the non-linear model parameter variance to determine a truncation level of the SVD for which the resolution kernel of the considered cell is computed. This approach allows to study efficiently the influence of non-linearity on model resolution, an aspect of geophysical inverse theory which largely has been neglected. The computation of the non-linear variance up to a pre-determined model variance threshold gives the truncation level of the SVD and fixes a point on the model resolution-covariance trade-off curve at the same time. From the resolution matrix corresponding to the thus determined truncation level, representative resolution measures in the form of a centre of resolution as well as horizontal and vertical resolution lengths are estimated. We have added our resolution and variance analysis to the prominent 2-D magnetotelluric inverse code REBOCC by Siripunvaraporn & Egbert (2000).
In the following, our non-linear analysis is introduced by first reviewing the theory of linear model variance and resolution estimates before the non-linear aspects are discussed. Short descriptions of model subspace techniques, resolution measures and a general scheme round off the theoretical part. Thereafter, a synthetic and a real world example are presented.

T H E O RY
In inverse problems, a vector of N measurement readings d = (d 1 , . . . , d N ) T with its corresponding standard errors σ = (σ 1 , . . . , σ N ) T is related to a vector of M model parameters m = (m 1 , . . . , m M ) T through a generally non-linear forward operator F: where the notation a W = Wa denotes a multiplication with the diagonal weighting matrix W = diag{σ −1 1 , . . . , σ −1 N } which contains the standard errors.
The error vector e describes the quality of the fit of the forward response F[m] of the current model m to the measured data d. The basic strategy of improving the fit is to minimize the cumulative normalized squared error Q = e T W e W with respect to the model parameters (Lines & Treitel 1984;Menke 1989), that is, to demand that the gradient of Q with respect to m vanishes, For this purpose, the forward operator needs to be linearized by expanding it in a Taylor series in the vicinity of the current model m (Lines & Treitel 1984): where J W is the normalized Jacobian matrix of partial derivatives atm. The parameter change vector m = m −m describes the change of model parameters away fromm. With the definition of the data discrepancy vector d = d − F[m] = e(m), the original problem of minimizing the non-linear quantity Q = e T W e W is substituted by minimizing the linearized quantity Q lin = (e lin W ) T e lin W = ( d W − J W m) T ( d W − J W m). According to Golub & van Loan (1996), a truncated singular value decomposition (TSVD) gives an approximate factorization of the normalized sensitivity matrix as where the columns of U p ∈ R N * p and V p ∈ R M * p consist of the p data and p model eigenvectors, respectively, of the p largest nonzero singular values. The singular values are positive and are ordered in the diagonal matrix p ∈ R p * p such that λ 1 ≥ · · · ≥ λ p > 0. The inclusion of only the p largest singular values in the solution process is often desirable in order to grant the stability of the solution. The factorization given by eq.
(3) is exact if all non-zero singular values are included. The truncated generalized inverse estimate of the minimum point of Q lin is (Lines & Treitel 1984;Menke 1989) where J −g With the definition of the generalized parameters α p = V T p m and generalized data β p = U T p d, the solution of the inverse problem (eq. 4) can be written as

Model resolution and covariance
The model resolution matrix R relates the model m est obtained from the TSVD inversion to the unknown true model m (Menke 1989): The kth row of the resolution matrix can be considered as a set of coefficients that describe with which weights the parameters of the true model m enter into the kth parameter of the estimated model m est . If, for instance, the kth element of the kth row is the only non-zero element, then only the kth parameter of the true model m will contribute to the kth parameter of the estimated model m est and the kth parameter is said to be perfectly resolved. In contrast to this, if there is a large number of non-zero elements in the kth row with approximately equal amplitudes, then the kth parameter of the estimated model m est will be an average over a wide range of model parameters and is said to be badly resolved. A resolution matrix identical to the identity matrix, that is, R = I, implies perfect resolution of all model parameters and can only be obtained for overdetermined least-squares problems (Menke 1989). The elements of the resolution matrix scale with the staggered grid chosen for solving the forward problem as can be seen from the relation to the continuous resolving kernel or averaging kernel r (Backus & Gilbert 1968) as defined below. In an integral formulation, the estimated model parameter m est at a position (y, z) is related to the continuous true model m on an area A (Backus & Gilbert 1968;Parker 1977) via m est (y, z) = A r ((y, z), (y , z ))m(y , z ) dy dz .
In the transition to a discrete case where the cell index k is assigned to the cell at the coordinate (y, z) and the cell index i is assigned to the cell at (y , z ), eq. (7) becomes where y i and z i are the width and height of cell number i, respectively. For cell number k, however, we have from eq. (6) Hence, the elements of the model resolution matrix R are related to the elements of the discretized resolving kernel by that is, the resolution matrix R contains samples of the resolving kernel r scaled by the cell area.
In the resolution plots in Sections 3 and 4, we depict the elements r ki of the resolving kernel of the kth cell scaled by its maximum modulus.
The model covariance matrix describes how the data errors given by the data covariance matrix [cov d] propagate into the errors of the estimated model m est (Menke 1989): In practice, [cov d] contains diagonal elements only, that is, the squares of the total data errors. However, since we normalized the Jacobian by the total data errors the model covariance matrix simplifies to and the standard errors ste of the model parameters are the square roots of the diagonal elements of eq. (12), that is,

Iterative computation
Due to the orthogonality of the model eigenvectors, the model resolution and covariance matrices can be computed iteratively with increasing singular value number p. The resolution matrix containing the first p + 1 model eigenvectors can then be expressed as the resolution matrix of the first p model eigenvectors plus the outer vector product of the ( p + 1)th model eigenvector: Increasing the number of model eigenvectors improves the resolution.
Similarly, the model covariance can be computed in an iterative fashion by: We are, however, only interested in the variance of each cell. For model cell number k we then obtain: where v k p+1 is the kth component of th ( p + 1)th model eigenvector.
Increasing the number of eigenvectors and singular values means to increase the variance.

Resolution-covariance trade-off
According to eqs (14) and (16), resolution improves at the same time as variance rises if the number of included eigenvalues and eigenvectors is increased. This behaviour indicates a trade-off between model resolution and covariance (Backus & Gilbert 1968;Parker 1977;Menke 1989). The trade-off justifies two different approaches. One could fix a resolution for the cells under consideration and then determine the variances of the cells or set a variance threshold for all cells and compute the resolution properties. In our approach, we keep the model parameter variances under a pre-determined threshold. According to eq. (16), this corresponds to computing a truncation level p cut−off of the SVD. For this truncation level, we then compute resolution properties from eq. (14).
In the REBOCC 2-D magnetotelluric inverse code (Siripunvaraporn & Egbert 2000), the model parameters are the logarithms of the resistivities ρ with respect to the base 10. This means that the considered variance estimate for the kth cell is the square of the logarithmic resistivity error, that is, var (m est k ) = ( log 10 ρ k ) 2 . The range of resistivities that we permit in the resolution analysis by defining a variance threshold is then given by log 10 ρ ± ( log 10 ρ) thresh . The error log 10 ρ can be understood as a factor f on the resistivity ρ by letting log 10 ρ = log 10 f . It follows log 10 ρ ± ( log 10 ρ) thresh = log 10 ρ ± log 10 f thresh = log 10 (ρ · f ±1 thresh ). Hence, resistivity values within the range ρ/ f thresh , . . ., f thresh · ρ given by a factor f thresh are allowed for in our model variance and resolution analysis.

Non-linear analysis
A more adequate choice of model variance and resolution estimates reflects the non-linearity of the underlying inverse problem at least partly. A simplified description of the non-linear confidence surface in parameter space is a pseudo-hyperellipsoid (Johansen 1977) whith semi-axes that reflect at which distance from the optimum model a given misfit is reached in direction of the model eigenvectors. In this context, we refer to the model with minimum misfit for a certain truncation level as the optimum model. The pseudohyperellipsoid can have different lengths of the semi-axes in the positive and negative direction of a particular model eigenvector. It thus resembles the hyperellipsoid of the linear confidence surface which has semi-axes of equal length.
For the linear case, we consider the misfit Q (see Section 2, 2nd paragraph) in the vicinity of the set of model parameters m 0 in which Q attains its minimum and let m =m 0 + m, (Johansen 1977). This corresponds to U T p e 0 W = 0 since all eigenvalues are nonzero. From eq. (2), we have Together with the compatibility condition, this yields In terms of the generalized parameters (see eq. 5), Q = m T V p Λ 2 p V T p m can be expressed as This is a hyperellipsoid with semi-axes s lin i = √ Q/λ i and axes parallel to the model eigenvectors (Johansen 1977). Hence, the linear confidence surface for a given Q has the shape of a hyperellipsoid in parameter space (Johansen 1977).
A simplified description of the non-linear confidence surface can be obtained by determining the lengths of the actual semi-axes s ± i in both positive and negative directions of an eigenvector v i (Johansen 1977;Pedersen 1979;Pedersen & Rasmussen 1989). The actual shape of the non-linear confidence surface in directions other than those of the eigenvectors is neglected and, therefore, approximated as a pseudo-hyperellipsoid. Fig. 1 depicts a pseudo-hyperellipsoid for the simple example of two model parameters. The actual semiaxes s ± 1 and s ± 2 correspond to the points in model space where the pseudo-hyperellipsoidal surface (thin line) intersects the true nonlinear confidence surface (thick line) for a given Q. The pseudohyperellipsoid consists of four ellipsoidal sections with semi-axes . The problem of computing the actual semi-axes can be formulated as finding the zeros of for all eigenvectors v i and m = m 0 ± s ± i v i . The desired deviation Q from the minimum misfit Q 0 needs to be specified. Eq. (22) is a non-linear equation in one unknown (either +s + i or −s − i ) and is solved with a safeguarded method that uses the secant method as a rapidly convergent method but keeps a bracket of the bisection method around the solution in case that the next approximate solution of the secant method falls outside the bracketing interval (Heath 2002). For an eigenvector v i with i > 1, the initial brackets of the bisection method are found by evaluating eq. (22) at ±s ± i = ±s ± i−1 and then either increasing or decreasing ±s ± i by a constant factor untilQ[m] changes its sign. For the first eigenvector v 1 , the search for the initial brackets is started from ±s lin 1 . Due to the non-linearity, the lengths of the non-linear and linear semiaxes can be different especially for the smallest singular values (Johansen 1977;Pedersen & Rasmussen 1989).
After determining the linear/non-linear semi-axes, the confidence limits of the model parameters can be computed by finding their extreme points on the surface of the (pseudo-)hyperellipsoid that is characterized by the specific choice of Q. The extreme points of a parameter m k are found by requiring that the gradient of Q is parallel to the normal n of the hyperplane m k = const (Johansen 1977;Pedersen 1979) since the hyperplane corresponding to n is tangential to the (pseudo-)hyperellipsoidal surface in the extreme points (Fig. 1). The kth component of n is set equal to +1 for the maximum point or −1 for the minimum point whereas all other components are set equal to zero. Furthermore, m needs to be on the (pseudo-)hyperellipsoidal surface characterized by Q. For the linear case, this gives (Johansen 1977;Pedersen 1979;Pedersen & Rasmussen 1989) For the choice Q = 1, eq. (23) simplifies to the standard error of the model parameters as in eq. (13) and corresponds to a confidence level of 68 per cent. The computation of the non-linear extremal points demands that the s ± i are substituted for the s lin i in eq. (23) according to the following rule (Johansen 1977). Given that the maximum point of m k on the pseudo-hyperellipsoid is to be determined, one has to substitute s + i for s lin the semi-axes s + i of the positive direction of v i must be substituted if n · v i > 0 and vice versa. A similar rule has to be applied for the minimum point of m k .
Only in the case where the estimated model m est equals the optimum model m 0 can it be expected that the non-linear and linear estimates of the semi-axes have equal values for the largest singular values. For the widely used smoothed inverse operators, the inverse process does not hit the optimum model of Q and m est is shifted away from m 0 . Hence, also the s + i and s − i will be systematically shifted away from s lin i in a way that if s + i < s lin i then s − i > s lin i and vice versa. Consequently, non-linear effects cannot be distinguished from a displacement from the optimum model m 0 . In our approach, we therefore, use a combination of mean-square error (MSE) inversion (Pedersen 2004) and TSVD inversion  to move to the optimum model before a final model variance and resolution analysis is carried out. During the construction of the inverse solution, it is not known how many pairs of non-linear semi-axes have to be included in the following non-linear analysis in order to reach the permitted variance threshold. For this reason, we perform an initial inversion with the MSE method, which usually includes less eigenvectors than required by the following initial non-linear analysis. Including more eigenvectors in the non-linear analysis means that certain eigenvectors are included that were considered part of the null-space (Menke 1989) during the preceding MSE inversion. Therefore, a final TSVD inversion is performed with a truncation level that equals that of the initial non-linear analysis and, after that, a final non-linear analysis is performed. Typically, the estimated truncation level does not change strongly from the initial to the final non-linear analysis and two iterations of model improvement and following non-linear analysis suffice in most cases.

Most-squares inversion
In our non-linear analysis, we investigate non-linearity only in the directions of the model eigenvectors. The trustworthiness of our variance estimates can be tested with a most-squares inversion (Jackson 1976;Meju & Hutton 1992) that computes the extreme resistivities of a cell for a given misfit threshold and truncation level.
If the true variance of the considered cell is larger than that from our non-linear analysis, then the given variance threshold would already have been reached with a smaller truncation level. Hence, too many model eigenvectors were included in the computation of the resolving kernel, and our resolution estimates are too optimistic.
The aim of the most-squares inversion is to to find the extreme values of a functional m T n under the constraint that the misfit Q does not exceed some threshold Q t . If n is chosen as a normal of the hyperplane m k = const, the kth model parameter will be maximized or minimized. The cost-functional of a stable iterative most-squares scheme can then be written as (Meju & Hutton 1992) where m = m 0 + m, 1 2μ is a Lagrange multiplier, β is a Marquardt damping factor, and L 2 constrains the energy of the model correction step m. In our case, the threshold misfit is chosen as Q t = Q[m 0 ]+ Q. The extremum of the cost functional (eq. 24) corresponds to a model correction step (Meju & Hutton 1992) which consists of a least-squares model correction m ls and a compensating term related to the maximized or minimized parameter. The Lagrange multiplier μ is determined by the quadratic constraint that the linearized misfit Q lin = ( d W − J W m) T ( d W − J W m) equals the threshold misfit, that is, Q lin = Q t . This constraint yields (Meju & Hutton 1992) where Q lin ls is the linearized misfit function pertaining to the leastsquares model correction m ls .
For a consistent comparison of the variance estimates derived from our non-linear analysis with the non-linear variances given by the most-squares inversion, the generalized inverse used in the most-squares inversion is constructed with the truncation level of the final non-linear analysis.
The damping factors must be chosen carefully since they influence the directions of the model correction steps. Although the damping factor can be chosen such that the confidence surface corresponding to Q 0 + Q is reached with one iteration, there is no guarantee that the resistivity of the considered cell is maximized or minimized by doing so. Therefore, an iterative implementation of the most-squares algorithm has to be applied. We limit the range of tested damping factors to those, which do not result in misfits exceeding Q 0 + Q, and pick the one which gives the largest ratio in resistivity change of the considered cell to average resistivity change of the model.

Model subspace
The use of a subspace technique is motivated by the fact that a SVD is very costly in terms of CPU time. In order to solve for singular values and model eigenvectors, the computing time scales as M 3 (Golub & van Loan 1996) where M is the number of model parameters. Hence, we restrict ourselves to considering a subspace around the model cell of interest. The extent of the subspace is determined by the shape of the subspace and the horizontal and vertical subspace radii s y and s z . In practice, subspaces of half-elliptical form, centred at the surface, and subspace radii chosen as the depth to the cell considered turned out to be most adequate. This choice is justified by the diffusive behaviour of the electromagnetic fields. With increasing depth the electromagnetic fields spread over a larger volume and resolution is diminished. An overview of different subspace choices is given in Appendix A.

Resolution measures
In order to obtain a quick overview of the resolution properties, it is desirable to condense the resolving kernel into a reduced set of resolution measures. Alumbaugh & Newman (2000) propose the use of a 50 per cent spread width that is the horizontal or vertical width of the resolving kernel between points at which the amplitude is 50 per cent of its maximum. Friedel (2003) reduces the information contained in the model resolution matrix to a few representative measures such as a radius of resolution derived from the diagonal elements of the resolution matrix and a distortion flag that describes whether the position of maximum resolution for a certain cell is offset from the cell itself. Both approaches deliver, however, limited resolution measures since the possible existence of side-lobes in the resolving kernel is not allowed for. Also, the position of the centre of resolution with respect to the considered cell is not incorporated. We follow the definition of resolution centres and lengths as integrals over the resolving kernel by Pedersen (1977). This choice takes account of the shape of the complete resolving kernel on the subspace. Details can be found in Appendix B.

General scheme
The general scheme of our TSVD resolution and variance analysis can be summarized as follows.
(1) Select the mesh cells for which TSVD resolution/variance analyses are to be carried out.
(2) Select a model variance threshold in form of a factor f thresh on the model resistivities (see end of Section 2.1.2).
(3) For each selected cell: (i) Determine the radii of the subspace around the selected cell, for example, on the basis of the depth to the cell (see Section 2.5 and Appendix A).
(ii) For analysis step = 1,2 (1) If analysis step = 1, perform a MSE inversion on the subspace giving a truncation level p MSE (see end of Section 2.2). If analysis step = 2, perform a truncated SVD inversion on the subspace with a truncation level p TSVD = p 1 of analysis step 1 (see end of Section 2.2).
(2) Scale subspace Jacobian of last model by total data errors.
(4) Determine the non-linear semi-axes in direction of the model eigenvectors with a combination of the secant and bisection methods (Section 2.2) and substitute the non-linear semi-axes for the linear semi-axes in eq. (23) for the model parameter variance until the given variance threshold is reached. This determines the truncation level p 1,2 of analysis step 1 or 2, respectively, and typically p 1 > p MSE .
(iii) For the final truncation level p 2 of analysis step 2, compute the subspace resolution matrix (eq. 14) and the resolution measures (Section 2.5). (iv) Check the model variances with a most-squares inversion (Section 2.3).

A S Y N T H E T I C E X A M P L E
Synthetic forward data of a simple 2-D model with a 10 m block in a homogeneous half-space of 1000 m (Fig. 2) were computed for the transverse magnetic (TM) and transverse electric (TE) modes (see e.g. Brewitt-Taylor & Weaver 1976; Aprea et al. 1997). The responses were computed at 10 receiver sites for 12 periods ranging from 0.0031 to 6.4 s giving a total of 480 data points. Gaussian white noise corresponding to 2.5 per cent of the modulus of the computed impedances was added to the forward responses of both polarizations.
After that, an inversion of the synthetic data set was performed with the REBOCC inverse scheme. The error floor was assumed to  correspond to 2.5 per cent of the modulus of the impedances, and the starting model was a homogeneous half-space of 300 m. At the end of the seventh iteration a model was obtained (Fig. 3) which fits the data to a rms misfit of 0.88. Additional iterations with REBOCC did not decrease the rms misfit further.
On the basis of this REBOCC model, the TSVD resolution and variance analysis was performed on a set of three cells where the first cell is located just on top of the conductive block (marked as cell I in Fig. 3), the second cell is located aside the conductor (cell II in Fig. 3) and the last cell lies just under the conductive block (cell III in Fig. 3). The shape of the subspaces was in all cases chosen to be a half-ellipse (see eq. A1). The horizontal and vertical subspace radii s y and s z , respectively, were in all cases chosen as the depth to the centre of the cell. The results of the analysis are summarized in Table 1.
First, the variance and resolution properties of cell I on top of the conductive block are discussed. A total of 168 cells enters into the subspace. In Fig. 4, the resolving kernel for cell I is shown. Cell I is marked with a white diamond, and the border of the subspace is illustrated with a dashed line. The values of the resolving kernel of cells pertaining to the subspace are plotted according to the normalized colour scale to the right of the plot. The initial MSE inversion gives a conservative truncation level of 71. For the resolution and variance analysis, the relative model variance threshold was set to f thresh = 3.0. In both the initial and final variance analysis, the total of 168 eigenvectors and semi-axes is exploited without exceeding  the given variance threshold. The linear and non-linear semi-axes are depicted in Fig. 5. The linear semi-axes (red curve) have a dynamic range between ≈2 × 10 −2 and ≈3 × 10 +3 whereas the non-linear semi-axes (blue and green curves) do not exceed 1 in magnitude. Up to singular value number 50, the linear and non-linear semi-axes are almost equal indicating that non-linearity becomes important only for higher singular value numbers. Since the non-linear semi-axes increase very slowly for singular value numbers higher than 50, it is evident that a higher number of eigenvectors and semi-axes can be included during the computation of the variance if non-linearity is allowed for. Consequently, more model eigenvectors enter into non-linear resolution estimates than into linear resolution estimates and non-linearity has the property of focusing resolution. As a rule of thumb, we found that the SVD needs to be truncated when the semi-axes are on the order of 1 if variances corresponding to factors f thresh = 3 are demanded. As the non-linear semi-axes are smaller than 1, this threshold cannot be reached even if all 168 eigenvectors and singular values are included as is the case here. The resulting factors on the resistivity of the investigated cell are f − = 1.50 in direction of decreasing and f + = 1.49 in direction of increasing resistivity. Since all eigenvectors are included, the cell is uniquely  After singular value number 50, the non-linear semi-axes increase very little, whereas the linear semi-axes continue to increase. Hence, more semi-axes and eigenvectors can be included in the computation of non-linear model variances until a given variance threshold is reached, and non-linearity tends to focus resolution. resolved and the resolving kernel is different from zero only in cell I itself (Fig. 4). This result is confirmed by the most-squares inversion which gives factors f − MSQ = 2.06 and f + MSQ = 1.47, respectively. Although the former factor is is somewhat higher than its non-linear equivalent f − , it clearly stays below the given threshold. For comparative purposes, the results of a linear model variance and resolution analysis of cell I are described. With only 83 of a total of 168 model eigenvectors and linear semi-axes, a factor f lin = 2.66 on the resistivity of cell I comes closest to f thresh = 3.0. The linear semi-axes computed during the linear analysis and the linear semi-axes computed during the non-linear analysis (Fig. 5) are rather similar. The length of the linear semi-axis for the truncation level p cut-off = 83 in Fig. 5 is approximately 2, which conforms to the above rule of thumb. According to the linear analysis, the elements of the resolving kernel in Fig. 6 are non-zero outside cell I, and the estimated resistivity of cell I is predominantly an average of resistivities of cells close to cell I and of cells close to the surface. However, also the linear resolution estimate shows that cell I is rather well resolved. In Fig. 6, the red bars outline the estimated vertical and horizontal resolution lengths, which are computed according to eq. (B2). The point where the bars cross is the centre of resolution calculated from eq. (B1). The centre of resolution coincides with cell I, and the resolution lengths are approximately a factor 3 larger than the dimensions of cell I. In comparison to non-linear resolution estimates, the linear resolution estimates of relatively badly resolved cells tend to deteriorate more drastically than those of rather well resolved cells. This is basically due to the fact, that the elements in the model eigenvectors of rather well resolved cells are large in the first model eigenvectors, which already enter in a linear analysis, and small in later model eigenvectors (see examples in . In contrast to this, the elements of badly resolved cells appear to be more concentrated in model eigenvectors with a high singular value number, which perhaps do not enter in a linear analysis. Next, cell II within the resistive host to the left of the conductive block is analysed. The corresponding subspace contains 167 cells  Figure 6. Subspace resolution kernel for cell I as obtained from a linear model variance and resolution analysis. In contrast to the findings of the non-linear analysis, the resistivity of cell I is not perfectly resolved (see Fig. 4) but an average of resistivities of cells in the vicinity of cell I and of cells close to the surface. The values of the resolving kernel represent the weights with which the resistivity of each cell enters into the average. However, also the linear resolution estimate shows that cell I is rather well resolved. The red bars depict the estimated vertical and horizontal resolution lengths (see eq. B2). The point where the bars cross defines the centre of resolution according to eq. (B1).  (Fig. 7). The truncation level of the initial MSE inversion is 49. The following initial variance analysis predicts a truncation level of 133 for resistivity factors slightly smaller than 3 (see Table 1). This estimated truncation level remains remarkably constant after the final truncated SVD inversion and final variance analysis. The final variance analysis yields factors of f ± = 2.99 in directions of both decreasing and increasing resistivity. Both vertical and horizontal resolution lengths are about a factor of 4 larger than the corresponding cell dimensions. The centre of resolution coincides with the considered cell. The variances are, however, underestimated as is proven by a most-squares inversion, which reports the variability as f − MSQ = 4.70 and f + MSQ = 18.4. Hence, our resolution estimates are too optimistic, and the resolving kernel corresponding to a factor f ± ≈ 3 is less concentrated on the considered cell than shown in Fig. 7.
For EM induction processes, it is generally accepted that the resolution below a thick conductor is relatively poor due to the strong  . 2-D model obtained from the final truncated SVD inversion in the subspace pertaining to cell III. The model fits the synthetic data to a rms misfit of 0.86. As an effect of the reduced misfit obtained by the TSVD inversion, this model is slightly rougher than the REBOCC model in Fig. 3. The cell on top with a black filling has a resistivity of more than 1800 m as indicated by the legend.
damping of the electromagnetic fields in the conductor itself. We investigate the resolution properties of cell III centred below the conductive block. The resolving kernel is shown in Fig. 8. For a maximal period of 6.4 s and a half-space resistivity of 10 m the skin depth is about 4000 m. Since the conductive block has a thickness of about 5800 m, it can be expected that there is relatively little resolving power for the investigated cell. A total of 284 cells is included in the subspace. The extension of the subspace incorporates the conductive block completely (see Fig. 8). The model of the final TSVD subspace inversion is displayed in Fig. 9. In comparison with the smooth REBOCC model (Fig. 3), the subspace model has a slightly improved rms data fit of 0.86. The upper and lateral boundaries of the block appear to be equally well determined in both models. Also, in both models the lower edge is smeared out towards greater depth. There is no mode or period range in which one of the models shows a systematic improvement of data fit over the other model (not shown). The number of included eigenvectors is the same for the initial and final variance analysis (Table 1), that is, 245. The final variance analysis yields factors of f − = 2.86 in direction of decreasing and f + = 2.88 in direction of increasing resistivity. In accordance with our expectation, the resolving kernel (Fig. 8) is strongly distorted. Many cells close to the surface contribute to the resolving kernel indicating that the estimated resistivity of cell III is an average over a large area. Fig. 4 suggests that the cells in the upper part of the conductive block are very well resolved, while resolution deteriorates for cells deeper in the block. As the resolution matrix is symmetric, the contribution of the cells in the upper part of the block to the resistivity of cell III is negligible. This results in small elements of such cells for the resolving kernel of cell III (see Fig. 8). Cells in the deeper part of the conductive block are less resolved than those in the upper part meaning that their resolving kernels have non-zero values spread over cells in their vicinity including those below the block. Consequently, cells of the estimated model below the conductive block will exhibit lower resistivity than the corresponding cells of the true model, and thick conductors smear out towards depth in the estimated model (see Fig. 8). Due to the large resolving kernel values of the cells close to the surface, the centre of resolution is displaced from the investigated cell, and the horizontal and vertical resolution lengths exceed the dimensions of the cell. The most-squares inversion gives factors f − MSQ = 2.73 and f + MSQ = 2.75 (Table 1) which confirm our non-linear variance estimates.  report on the problem of determining the geometry of a glacio-fluvial aquifer system composed of sand and gravel formations in Skediga (Sweden) with the controlled source radiomagnetotellurics method. The lower boundary of the aquifer system is the crystalline basement, and it is overlain by a system of clay lenses. Due to the conductive clay lenses, it is difficult to determine the structure of the aquifer and especially the transition into the resistive basement.

A F I E L D D ATA E X A M P L E F RO M T H E S K E D I G A A R E A ( S W E D E N )
The model, that is arrived at after four inverse steps using the determinant mode (Pedersen & Engels 2005) in the REBOCC inverse code (Siripunvaraporn & Egbert 2000), is shown in Fig. 10. The number of data points is 528, and the rms fit is 1.47. Between profile metres 10 and 100 as well as 130 and 220, there are two conductive clay lenses with average resistivities of 10 m (in yellow colours).    (i.e. the transition between the two greenish colours as log 10 30 ∼ = 1.5) as the lower bound of the clay lenses. The average resistivity of the underlying aquifer is 100 m. At about profile metre 115, the sand and gravel formations of the aquifer stretch up to the surface. According to boreholes in the vicinity of the profile, the transition from the aquifer to the underlying crystalline basement occurs at about 30 m depth. At a few boreholes, the depth to the basement correlates very well with the 300 m contour line in Fig. 10 and suggests a basement resistivity larger than 300 m. The variance and resolution properties of the four cells labelled A-D in Fig. 10 are investigated in subspaces which have the shape of half-ellipses (see Appendix A). The results of the variance analyses and most-squares inversions are summarized in Table 2. In order to clearly distinguish the resistivity ranges permitted for the clay lenes and the aquifer, a resistivity factor threshold of f thresh = 3.0 is used. For this choice, the resistivity range of the clay lenses is {3.3, . . . , 30 m} at an average resistivity of 10 m whereas the resistivity range of the aquifer is {33, . . . , 300 m} at an average resistivity of 100 m. Hence, the permitted resistivity ranges do not overlap. Fig. 11 shows the resolution properties of cell A which is located inside a clay lens. The number of cells in the subspace is 184. After including all non-linear semi-axes in the initial analysis, perfect cell resolution and a variance corresponding to a factor of f + = 1.75 in  Figure 11. Subspace resolution kernel of cell A (marked with a white diamond), which is situated inside a clay lens, overlain by log 10 ρ(y, z)-isolines of the model shown in Fig. 10. The resolution values outside cell A are zero, and cell A is perfectly resolved. For cell B, which is deemed to lie in the upper few metres of the basement and underneath a clay lens, the resolving kernel is shown in Fig. 12. The initial variance analysis gives factors f ± = 2.93 for 214 of a total of 260 eigenvectors. The number of eigen pairs remains remarkably stable in the final variance analysis yielding 215 included eigenvectors and factors f − = 2.81 and f + = 2.84, respectively. The resulting resolving kernel (Fig. 12) spreads over five cells in both vertical and horizontal direction and is centred around the cell of interest thereby suggesting moderate resolution. However, the most-squares inversion gives a factor of f − MSQ = 24.04 in direction of decreasing resistivity and f + MSQ = 46.59 in direction of increasing resistivity meaning that the resolution is grossly overestimated. Less severe damping and hence better resolution of the transition between aquifer and basement can be expected in the area where the sand and gravel formations stretch up to the surface. We investigate cell C in the top section of the window between the clay lenses (Fig. 13). The initial variance analysis gives a truncation level of 142 of a total of 192 eigenvectors which is augmented to 153 by the final variance analysis. Factors of f − = 2.59 and f + = 2.95 are obtained in direction of decreasing and increasing resistivity, respectively. The most-squares inversion confirms our non-linear estimates by giving a factor f − MSQ = 3.88 in direction of decreasing resistivity and of only f + MSQ = 1.53 in direction of increasing resistivity. Hence, the resistivity is better determined in the latter direction which can be attributed to the TM mode which enters into the determinant.
A similar behaviour is observed for cell D on top of the basement below the window between the clay lenses (Fig. 14). The final variance analyses yields a truncation level of 222 with final factors f − = 2.85 and f + = 2.80. The following most-squares inversion reveals that the resistivity of this cell is also better constrained in direction of increasing resistivity with factors f − MSQ = 2.03 and f + MSQ = 1.45. In fact, the given variance threshold is not reached as indicated by the results of the most-squares inversion and the resolution of cell  Figure 14. Subspace resolution kernel of cell D under a window between the clay lenses and within the upper few metres of the basement. The spread of the resolution kernel resembles that for cell B (Fig. 12). In contrast to cell B, the variance estimates from the most-squares inversion are closer to the non-linear estimates (see Table 2) indicating that cell D is better resolved than cell B.  Figure 15. Subspace resolution kernel of cell D as obtained from a linear analysis. The resolution kernel is more spread out than that of its non-linear analogue (Fig. 14) since model eigenvectors with essential contributions to cell D are not included. D is possibly underestimated. A comparison with the results of a linear model variance and resolution analysis underlines the importance of allowing for non-linearity. During the construction of the subspace model and the linear analysis, a truncation level p cut-off = 95 gives a linear factor f lin = 2.97 on the resistivity of cell D. According to the resolving kernel ( Fig. 15) from this linear analysis, the estimated resistivity of cell D is no longer an average of cells that belong solely to the aquifer and upper part of the basement but also includes cells from the clay lenses and cells close to the surface. As explained in Section 3, the exclusion of model eigenvectors with essential elements for cell D has lead to a significantly deteriorated resolution estimate of cell D.

D I S C U S S I O N A N D C O N C L U S I O N S
The method described here is a valuable means of analysing model resolution and variance properties of a class of non-linear inverse problems, which do not explicitly include a priori models and model covariances. With a few minor modifications, our method can be applied to a wide range of non-linear inverse problems of this class. A point on the resolution-covariance trade-off curve for a given model parameter is selected by fixing the model variance. The resolution and variance estimates are only influenced by the truncation of the SVD, and the truncation is solely governed by the data, the prescribed model variance and the non-linearity of the inverse problem. Our variance estimates partly account for the non-linearity of the inverse problem by systematically replacing the inverse singular values with the non-linear semi-axes in the direction of the model eigenvectors. The rotation of the model eigenvectors with position in parameter space is, however, not allowed for.
Inverse schemes based on Tikhonov regularization, that includes a priori information, are commonly used. This preference is partly due to the important computational savings in comparison with inverse schemes based on the SVD. In this context, the question must be raised whether smoothing inverse operators and a priori models are required by the operating geophysicist or whether they are used as a convenient way to reduce ill-posedness and facilitate a solution of the problem.
In the former case, that a priori models and a priori model covariances are integral part of the stochastic formulation of the inverse problem, the a posteriori model covariance and resolution matrices are directly given by well known formulae (Menke 1989).
In the latter case, when a smoothing inverse operator is used merely to construct a solution, that fits the data to some degree, and model smoothness is not a main aim, our model variance and resolution analysis can be a valuable alternative.
In assessing the resolution and variance properties of an estimated model of electrical resistivity, it is necessary to include non-linearity. As can be seen from the spectra in Fig. 5, the linear and non-linear semi-axes deviate strongly for small singular values and the volume of the non-linear pseudo-hyperellipsoid is much smaller than that of the corresponding linear hyperellipsoid. Hence, more eigenvectors can be included in a non-linear variance analysis until a given variance threshold is reached and non-linearity has the property of focusing resolution.
It is very often a critical task to assess the resolution of resistive structures in a model, especially if they are located under a conductive cover. The choice of truncation level tends to affect the resolution estimates of conductive and resistive structures differently. Conductive near surface structures have large elements in model eigenvectors with small singular value numbers and, hence, are estimated to have good resolution even if only a relatively small number of model eigenvectors is included in, for instance, a linear analysis. As we have demonstrated, the resolution of resistive structures is significantly underestimated if a linear analysis is employed, as such structures tend to have large elements in model eigenvectors with intermediate to high singular value numbers.
The reliability of our variance analysis can be tested with the most-squares inversion, which delivers an independent and improved estimate of the model variability. According to the results of the most-squares analysis, the variance threshold, that is used in our non-linear analysis, can be adjusted to obtain a most-squares estimate that is closer to some variance deemed appropriate for the problem. In the case of cell B of the Skediga model, for instance, the variance threshold f thresh should be decreased below 3.0 to approach f ± MSQ ≈ 3.0 more closely. By doing so, the truncation level will be decreased and an improved model resolution estimate will be obtained.
For a small-scale inverse problem like our RMT field example with about 2800 model blocks and 500 data points, a model variance analysis for a subspace of approximately 200 cells takes about one hour on a modern 32-bit processor with a CPU frequency of 2 GHz. The main contributions to the computing time are the calculation of the Jacobian matrix and its subspace SVD during the subspace inversion and the large number of forward computations during the search for the non-linear semi-axes. For 3-D or large 2-D problems, the time to perform a forward calculation and to compute the Jacobian matrix becomes very long and memory limitations restrict the ability to store the complete Jacobian matrix in memory. In order to make our method more efficient it would be desirable to compute a Jacobian for the subspace only.

A C K N O W L E D G M E N T S
Weerachai Siripunvaraporn of Mahidol University, Thailand and Gary Egbert of Oregon State University made the REBOCC code available for modifications. We would also like to thank three reviewers and the editor for their comments which helped to improve the manuscript.