A simple method for determining the spatial resolution of a general inverse problem

Summary

1 Introduction

A class of geophysical problems involves retrieving the covered structure m (an m-vector including: m  ₁, m  ₂, …, m_m) from an outside observation d (an n-vector) on the basis of a forward equation d=g(m). This is a typical inverse problem that yields a solution m=g  ^−g(d), where m (m  ₁, m  ₂, …, m  _m) is the expected solution vector and g  ^−g is the true or generalized inverse of the observation operator g. If the observation error δd is given, the solution uncertainty δm can also be obtained. The determination or inversion of m has been a significant focus of mathematical inverse studies for quite some time. Textbooks are readily available (e.g. Menke 1984; Tarantola 1987, 2004; Parker 1994; Aster 2005), discussing methods for determining the solution efficiently or at least acceptably. For example, for a very large and sparse linear/linearized equation, a solution may be obtained quickly by least square QR-factorization (LSQR) (Paige & Saunders 1982a,b), which has become a standard algorithm for solving tomographic problems (Zhang & Thurber 2007); even global solutions to complex non-linear problems may be efficiently obtained via optimization methods (e.g. Goldberg 1992; Sambridge 1999; An & Assumpção 2004; Lawrence & Wiens 2004).

The determination of both m and δm are insufficient for a real geophysical study. In geophysics, each parameter (e.g. m_i) of a model m represents information at a certain spatial or temporal position, expressed as m_i(x, y, z, t). The spatial and temporal resolutions (or the resolution length), that is, the size of the smallest possible feature that can be detected, as a function of the position (x, y, z, t) are relevant quantities and are informative for appraising the solution m as well as the solution uncertainty δm. The computation of solution's spatial resolution is non-trivial, and is more difficult than solving an inverse problem. Most geophysical studies, except for tomographic studies, almost uniformly neglect the calculation of a practical spatial resolution.

A generally preferable spatial-resolution estimate, widely used in seismic tomography studies, involves visual inspection of the restoration of a synthetic structure (e.g. checkerboard tests; Lévěque 1993; Thurber & Ritsema 2009; Feng & An 2010). The average resolution length in such cases is taken to be equal to half of the recovered checker dimension (Lebedev & Nolet 2003). The observation of small-scale heterogeneities in a visual inspection assessment is too easily mistaken as an indicator of a high resolution, and the creation of synthetic structures required by the tests becomes difficult for models with complex discretization.

An effective strategy for obtaining the model resolution-length scale is the use of Backus–Gilbert resolution kernels (Backus & Gilbert 1968), also referred to as a resolution matrix. Backus–Gilbert kernels can indicate the finiteness of a model parametrization, the incompleteness of data coverage and the effects of damping applied during the inversion (Ritsema 2004). The kernels have a finite spatial extent (Ritsema 2004) corresponding to the resolution length. The resolution matrix can be estimated easily for a moderate inversion (Berryman 2000a,b), and this is also feasible for large-scale problems through Lanczos iteration inversion (Jackson 1972; Yao 1999; Zhang & Thurber 2007), the truncated singular value decomposition (SVD) inversion (Aster 2005; Kalscheuer & Pedersen 2007), the Choleskey factorization inversion (Boschi 2003; Soldati & Boschi 2005) or the resolution-matrix approximation (Nolet 1999; Vasco 2003).

Very large inverse problems can be solved using efficient methods, for example, LSQR (Paige & Saunders 1982b); however, estimating the resolution matrix of a large problem is non-trivial, and not all resolution matrices can provide resolution-length information. The Backus–Gilbert resolution kernels treat essentially underdetermined linear problems (Tarantola 2004), and every estimate requires a full inversion of the system (Nolet 2008). A formal resolution test, therefore, becomes impractical for extremely large problems (Thurber & Ritsema 2009). In reality, the main utility of a resolution matrix is to measure the resolution obtainable from a data set according to the degree to which the resolution matrix approximates the identity matrix (Jackson 1972; Berryman 2000b). If the resolution matrix is an identity matrix, the solution m can be taken as unique, each element is perfectly resolved (Jackson 1972; Tarantola 2004), and the resolution length estimated from the resolution matrix should only be determined by the parameter's spatial extent.

Several researchers have published other skills to calculate the spatial resolution in typical inversions or tomographic studies. Derivative weighted sums measure the sampling of each node; they have long been used as an approximate measure of resolution (Thurber 1983; Toomey & Foulger 1989; Thurber & Eberhart-Phillips 1999; Zhang & Thurber 2007). Yanovskaya & Kozhevnikov (2003) estimated the radius of the averaging area in surface wave tomography to evaluate the resolution of the associated data. Fichtner & Trampert (2011) explained how to retrieve resolution information from the Hessian operator of an inverse problem.

Despite being crucial for solution appraisal, memory-efficiency estimates of the spatial resolution of a general inversion have yet to be developed. Here, a simple method is presented for inverting a solution's spatial/temporal resolution on the basis of a limited number of forward and inversion calculations.

2 Traditional Resolution Matrices and Resolution Length

Spatial-resolution calculations often arise from an analysis of a resolution matrix (Nolet 2008; Thurber & Ritsema 2009). Resolution matrices described in previous publications can be divided into three classes, as introduced below.

Consider a linear inverse problem in which the forward and inverse equations can be expressed, respectively, by eqs (1) and (2):

(1)

 

(2)

where G  ^{−  g} is a true or generalized inverse of an n-by-m observation matrix G. Replacing d in eq. (2) by the expression in the forward eq. (1) yields (Jackson 1972; Menke 1989):

(3)

and

(4)

where R is an m-by-m matrix referred to as the model resolution matrix, resolution operator or resolving kernels. Here, R is referred to as the direct resolution matrix. Eq. (3), which embodies the basic assumption underlying the Backus–Gilbert inversion method (Backus & Gilbert 1968), defines a linear projection of which solution m can be a local average calculated by integrating the model m. If G  ^{−  g} is not calculated, SVD analysis gives the solution (Jackson 1972):

(5)

where V is unitary.

The main utility of a resolution matrix is to provide a measure of the resolution obtainable from the data, and this measure is based on the degree to which the resolution matrix approximates the identity matrix (Jackson 1972; Berryman 2000b). If G is not a full rank matrix, R will not be an identity matrix, which indicates that some parameter (e.g. m_i) in the model m cannot be constrained well by the observation. The inverse problem under such conditions is ill-posed. Figs 2(a) and (b) show the resolution matrix R calculated by SVD for a simple underdetermined inverse problem, shown in Fig. 1. The inverse problem gives a null value for parameters in the range 86–100. These parameters are totally unresolved. If the observation matrix G is a full rank matrix, R should be an identity matrix because it is the product of a full rank matrix G and its inverse. An identity resolution matrix such as that shown in Figs 2(c) and (d) indicates that the model was recovered exactly and that the inversion was perfect (Jackson 1972; Tarantola 2004).

Figure 1

Illustration of a simple 1-D ray-propagation 100-parameter model (a) and its observation matrix (b). All model parameters were separated with a distance interval of 1 km. Only five observations/rays were used, no observational constraints were imposed on the parameters corresponding to distances of 85–100 km in the observations; therefore, the inverse problem was underdetermined and the observation matrix in (b) was not a full rank matrix.

Open in new tabDownload slide

Figure 2

Three resolution matrices for the underdetermined inverse problem presented in Fig. 1. Panels (b), (d) and (f) are the normalized matrices corresponding, respectively, to the panels (a), (c) (e). The matrix in (c,d) is an identity matrix. The matrices are calculated by the SVD of (a,b) the G matrix and its inverse, (c,d) the A matrix and its inverse and (e,f) the G matrix and the inverse of the A matrix. A is a matrix that combines G and the flatness constraints C among all neighbouring parameters. The weight (λ= 1) is used here. In panel (a), all components of the resolution matrix for parameters, m_i (i= 86–100) are 0 because no information is available for the parameters in the observation matrix G. The extent of the square-like pattern indicates that the related solution parameters should have the same value in this inversion problem. The extent can be taken as the resolution length of the solution parameters. The horizontal red bars in (d,f) of each model parameter are the same as in (b), and indicate the resolution lengths from (a), but are centralized along the diagonal element of the matrix.

Open in new tabDownload slide

If a resolution matrix is a nearly diagonal matrix, each element (e.g. m  _i) is in fact the weighted sum of nearby values m_j, where m_j are near m_i based on eq. (3) (Jackson 1972). The resolution matrix can, therefore, provide some information about the resolution length. For a full rank observation matrix G, the identity matrix R (e.g. in Figs 2c and d) indicates that the resolution length at a position m  _i only corresponds to the spatial extent of m  _i. Considering that the half width of the recovered checker extent can be taken as the resolution length (Lebedev & Nolet 2003), the half spatial extent of a parameter, m  _i, can also be taken as the resolution length of the parameter, given an identity resolution matrix. For the example in Figs 2(a) and (b), the shaded squares indicate that the parameters inside each square are constrained by a single observation and should have the same inverted values. The resolution lengths of the parameters located inside each square are half of the square width and they are marked in Fig. 2(b) by horizontal red bars centred on the diagonal elements of the matrix.

For an ill-posed problem, such as the underdetermined problem shown in Fig. 1, regularization (e.g. damping, smoothness/flatness or other a priori constraints, here represented by the k-vector c and by the k-by-m matrix C) must be used to stabilize the inversion. The forward equation then assumes a form similar to eq. (6):

(6)

and

(7)

where A is an (n+k)-by-m matrix, and λ is the weight required to balance the a priori constraints and observations in the inversion. The solution to eq. (6) can be expressed by

(8)

In an approach similar to that used to obtain the resolution matrix R in eq. (4), the forward and inverse eqs (6) and (8) yield a new resolution matrix R (referred to here as a regularized resolution matrix),

(9)

 Many inversion applications use existing inversion programs directly. For example, iterative inversion systems such as LSQR (Paige & Saunders 1982a,b) are often directly used in tomographic studies. Because inversion programs only provide limited options to include a priori constraints (e.g. in LSQR only damping is provided), if one wants to use one's preferred constraints, the simplest way is to take the regularized vector b and matrix A, as expressed in eq. (6), as the input observation vector and observation matrix in the inversion programs, respectively. Subsequently the resolution matrices provided by the programs are examples of a regularized resolution matrix R. Basically, stabilization on the eq. (6) requires, at a minimum, that A is a full rank matrix. R should then be an m-by-m identity matrix (e.g. in Fig. 2c) because R is the product of the full rank A and its inverse. The identity resolution matrix R (e.g. in Fig. 2c) indicates that the resolution length of each parameter is equal to half the spatial extent of the parameter; however, the resolution matrix R from the original forward and inverse equations indicated that some model parameters could not be well solved, and none of the parameters (Fig. 2a) have a resolution length equal to half their spatial extent. The resolution matrix R cannot, therefore, provide any resolution-length information. On the other hand, if we still use the original observation vector d and matrix G in an inversion program without any constraints, the resulting resolution matrix R is meaningful, but the inversion itself is ill-posed.

The vector c is often an empty vector, the components of which are zero. In this case, c can be ignored, and the vector b in the inverse eq. (8) can be replaced by d. In this case, the solution to eq. (6) will be identical to the solution to eq. (10),

(10)

If C is the identity matrix, eq. (10) assumes the form of typical Tikhonov regularizations. Furthermore a third class of resolution matrix () may be obtained based on the forward eq. (1) and inverse eq. (8),

(11)

Here, is referred to as a hybrid resolution matrix because it combines the original observation matrix G and another matrix, the regularized observation matrix A. The resolution matrices in many approaches (e.g. Crosson 1976), especially approaches that offer spatial-resolution information (e.g. Barmin 2001; Boschi 2003; Soldati & Boschi 2005), belong to this class. Obviously, even for an overdetermined problem, the hybrid resolution matrix is seldom an identity matrix because it is the product of the observation matrix G and the inverse of another matrix A, and the difference between A and G often cannot be ignored. By the same reasoning, can provide information about the original observation matrix G, the a priori matrix C, and the balance weight λ. The hybrid resolution matrix (Figs 2e and d) for the underdetermined problem shown in Fig. 1 can be very similar to R (Figs 2a and b).

Consider an overdetermined (n > m) problem based on the problem presented in Fig. 1. A new observation matrix G  _o is defined as the combination of the above G and the m-by-m identity matrix I,

(12)

Using this new observation matrix G  _o, the direct resolution matrix R is certainly an identity matrix (like Fig. 2c) because G  _o is a full rank matrix. Then each solution parameter can be exactly determined, even without a priori constraint. The resolution length at a position at which a parameter located is equal to half of the spatial extent of the parameter. R should also be equal to the identity matrix; however, the identity of with the identity matrix depends on the weight λ. For a small weight λ, the spatial resolution indicated by in Fig. 3(a) is similar to that expected from the identity matrix. For a larger λ, the hybrid resolution matrix (in Fig. 3b) ironically becomes very similar to the resolution matrix of an underdetermined problem, such as that presented in Fig. 2(e). This similarity mistakenly indicates that the resolution length is much larger than the real resolution length. Therefore, a poor value of λ may lead to an incorrect estimate of the resolution length based on a hybrid resolution matrix .

Figure 3

Hybrid resolution matrix calculated by the SVD for an overdetermined inverse problem using the weights (λ) of (a) 1 and (b) 10. The problem's observation matrix G_o is a matrix that combines G from Fig. 1 and the identity matrix. The regularized matrix A combines G  _o and C.

Open in new tabDownload slide

The above discussion suggests that a direct resolution matrix R can provide resolution-length information and the perfectness of the inversion; however, this information is seldom pursued in studies because the inversion is often ill-posed and must be regularized. The second class of resolution matrices, R, can evaluate the inversion perfectness for a given mathematic system, but such matrices cannot provide resolution-length information. The third class resolution matrices, , can be used to extract resolution-length information and to evaluate the perfectness of an inversion; however, these matrices depend on the regularization contribution to the inversion. From this perspective, it may be optimal to simultaneously calculate all three resolution matrices in a given study, although it is often impossible to perform the required matrix operations in the context of very large inverse problems. Below, I will introduce a new class of resolution matrices that are suited for very large inverse problems and do not need matrix operation.

3 Statistical Resolution Matrix and Resolution Length

3.1 Basic equations for a linear problem

Unlike previous treatments, the goal here is to directly invert for a resolution matrix based on the linear projection eq. (3). This treatment is distinct from that applied in the previous three classes of resolution matrices (R, R and ). Here, the inverted resolution matrix is represented by the symbol, ℜ, and the projection definition from the real model to a calculated solution is expressed as eq. (13),

(13)

Using eq. (13), any solution parameter m  _i can be expressed as:

(14)

where r_i,_j is the ith row and jth column element of ℜ. Each row (e.g. ith row) of ℜ is a resolution map that defines the contributions of all model parameters to the ith solution parameter. Alternatively, eq. (14) indicates that each solution parameter is generally derived from the weighted averages of the neighbouring model elements/parameters for a near-diagonal resolution matrix (Jackson 1972). A small r_i,j, indicates that the true model parameter m_j contributes little to the solution parameter m  _i.

A long resolution length indicates that neighbouring model element/parameter (m_j) at a long distance from m_i can contribute significantly to the parameter m  _i. In this case, r_i,i should be small according to eq. (14) or, say, m_i contributes little to m  _i. This implies that the amplitude uncertainty (δm  _i) of the solution parameter m  _i may be large. In general, a long resolution length may correspond to a large uncertainty in the anomaly amplitude, and vice versa.

The sum of the weights can be a constant, c_i, as expressed by:

(15)

 If the summation c_i is equal to 1, the solution parameter m  _i can be considered as unbiased under the condition that eqs (13) or (14) yield a true average (Nolet 2008). Otherwise, if c_i≠ 1, in eq. (14) is not the true average over the model parameters (Nolet 2008). In this case, the weighted average of the model elements/parameters should be written as:

(16)

3.2 Direction-independent or 1-D resolution matrix estimation

For a good inversion, the resolution matrix must be an identity matrix or nearly diagonal (Jackson 1972), that is, a given parameter (e.g. m_i) and its neighbours can provide a large contribution/weight to the average summation for a solution parameter (m  _i), and the weights (r_i,_j, j= 1, m) should decrease rapidly with the distance from m_i to another parameter (m_j); therefore, each row of the resolution matrix (r_i,_j, j= 1, m) can be approximated by either a cone-shaped function (Barmin 2001) or a Gaussian function (e.g. Nolet 2008). Here, the approximation of a Gaussian function is used. In one dimension, the approximate values of r_i,_j, the ith row and jth column of ℜ, in the Gaussian function can be written as

(17)

where x_i and x_j represent dimension (e.g. x, y, z or t) positions for m_i and m_j, L_i,j is the Euclidian distance between the dimensional positions of the parameters m_i and m_j, a_i is a constant, σ_i represents the half width of the peak at ∼60 per cent of its full height, and w_i (=1.17σ_i) is the half width at half maximum (and is often used to represent the shape's width). The spatial form (curve or surface) of a Gaussian function is very similar to the form of a recovered checker in a traditional checkerboard test. Since the half width of the recovered checker dimension can be taken as the resolution length (Lebedev & Nolet 2003), half of the Gaussian width, w_i, can be taken as the resolution length at the position m_i; therefore, w_i will be referred to as the resolution length at position i. If a_i is not considered, a normalized inverted resolution matrix can be constructed according to

(18)

Replacing r_i,_j in eq. (16) by the expression given in eq. (17) yields

(19)

 If we write

(20)

we obtain the equation

(21)

Because a resolution matrix does not depend on the specific model or observation values, but rather on the exclusive properties of the observation matrix G, we can use a random synthetic model to solve eq. (21). Given a random synthetic model m  ¹, we can calculate the synthetic observation from m  ¹ using eq. (1). We can then obtain a solution m  ¹ directly using the inverse eq. (8). For a very large inverse problem, the forward and inversion calculations (e.g. inversion by LSQR) are often much easier than to obtain the traditional resolution matrices. If the model and solution given in eq. (21) are known, w_i becomes a unique unknown in eq. (21). In this case, eq. (21) can be solved; however, the solution w_i may depend on the given m  ¹ and m  ¹. For ns (>1) given different random synthetic models M (=m  ¹, m  ², …, m  ^ns), we can obtain ns respective solutions M(=m  ¹, m  ², …, m  ^ns), and we can construct ns non-linear equations similar to eq. (21). w_i is the only unknown in the equations and can be obtained by minimizing eq. (22),

(22)

In most real inverse problem, w_i has a limited number of possible values. For example, in a 1-D inversion for an equal thickness layered model (e.g. see Fig. 1) with m parameters and thickness (h), w_i can be only one of 0.5h–mh. For the unique unknown eq. (22), an exhausted grid search is the most powerful and easiest approach to obtaining the resolution w_i. As ns increases, the dependence of w_i on the given synthetic model decreases, and the final solution w_i can be obtained by increasing ns until w_i becomes stable. By repeating the above grid-search inversion at parameter positions other than i, we can obtain all resolution lengths, w_i (i= 1, m). The full procedure outlined is indicatively illustrated by the flowchart in Fig. 4.

Figure 4

Flow chart of inversion for the resolution lengths. The procedure can be divided into two steps: (1) generate synthetic models and invert for solutions using your preferred constraints and inversion skill and (2) invert for the resolution lengths.

Open in new tabDownload slide

After obtaining all w_i (i= 1, m), we can directly construct a normalized resolution matrix using eq. (18). Note that we can obtain a_i using eqs (14) and (17), and we can then construct the inverted resolution matrix ℜ using eq. (17); however, because we mainly want to extract the resolution length from a resolution matrix, which w_i can provide directly, is mainly shown in the figures.

3.3 Method validity

The linear projection (eqs 3 and 13) from the real model to the calculated solution to a linear inverse problem is well accepted and widely applied in linear inverse studies. In contrast with the traditional resolution matrices obtained by operating on an observation matrix, the resolution matrix presented here may be solved directly using random synthetic models and their respective solutions. The solution to each synthetic model may be obtained by forward and inversion processes on the basis of a real observation matrix and other parameters (e.g. regularization), similar to the process used to invert real data. The solution is influenced by all parameters and methods used, from the forward calculation to the final inversion; however, the traditional resolution matrices are determined only using the observation matrix and a regularization procedure; therefore, the inverted resolution matrix presented here is not identical to the traditional resolution matrices and may include more information than is present in the traditional resolution matrices.

The flow chart (Fig. 4) and the associated discussion about the inversion procedure indicate that the skill used here may be considered as a simple application of matrix-probing techniques on the basis of randomization (e.g. Halko 2011; Chiu & Demanet 2012; Demanet 2012), while the inversion processes for the matrices ℜ and are more like a statistical analysis based on random samples (random models and their respective solutions). This is so, because the resulting solution (ℜ) depends on the sample number, which is introduced in detail later. Therefore, the inverted resolution matrix (ℜ) is referred to here as a statistical resolution matrix, and is referred to as the normalized statistical resolution matrix.

The above-mentioned approach was used to invert for the statistical resolution matrix (see Fig. 5) of the linear inverse problem presented in Fig. 1. 25 (ns= 25) pairs of random synthetic models and solutions were used to invert for the statistical resolution matrix. Fig. 5(b) showed that the resolution lengths at the positions for the parameters m_i(i= 1–85) were similar to those calculated directly from eq. (4). These values were reasonable at the positions for the parameters m_i(i= 86–100). Because the flatness constraint is associated with the parameters m_i(i= 86–100), similar to the neighbours m_i (i= 81–85), the resolution lengths for the positions of the parameters m_i(i= 86–100) should be 10 km (the half width of the spatial extent, 20 km) similar to the resolution lengths determined using the inversion approach of this study, as shown in Fig. 5(b).

Figure 5

(a) Statistical resolution matrix and (b) normalized matrix for the problem presented in Fig. 1. The statistical resolution matrix is inverted by the SVD using 25 pairs of synthetic models and solutions. The models were constructed by adding ≤10 per cent random deviations onto a reference constant. In panel (b), the red bars are the same as those shown in Fig. 2(b). The red contour marks indicate the half height position of the Gaussian shape, which is taken as the resolution length of the statistical resolution matrix.

Open in new tabDownload slide

During the extraction of resolution-length information from a resolution matrix, it is assumed that the diagonal components of each row should be largest, and that the values of the other components should decrease as the distance from the respective component to the diagonal component increases, as illustrated by the cone-shaped function (Barmin 2001) or the Gaussian function used here, and also by Nolet (2008). However, the real resolution matrix does not always behave in this way. The geometry of the resolution matrices shown in Figs 2(a) and (e) is neither symmetric nor Gaussian and differs from the above geometry variation assumption. Ritzwoller (2001) described a Gaussian smoothness constraint on the solution parameters associated with a surface wave tomography inversion and estimated the resolution lengths. This constraint forces the solution parameters to satisfy an equation similar to m=Cm, but the constraint does not directly influence the linear projection equation, m=  m, which was the assumption used to estimate the resolution length.

In general, the resolution length estimated from a resolution matrix based on a specific geometry shape can often provide an approximate indicator of the spatial resolution. The similarity between the resolution lengths from, respectively, R and ℜ, as shown in Fig. 5(b), indicates that the approximation can adequately represent the real resolution length for a linear inverse problem, at least for the inverse problem in Fig. 1.

3.4 The flatness influences the inverted resolution length

The influence of flatness was tested by estimating the statistical resolution matrices using different weights (λ) in eq. (7), which balanced the flatness constraint and the synthetic observations. Here, the forward eq. (6) was inverted using the LSQR, which is a conjugate gradient method. For λ= 1, the results shown in Figs 6(a) and (b) were identical to those shown in Fig. 5. Fig. 6 showed that the resolution lengths increased with the weights (λ) on the flatness values, as expected. In general, a large weight indicates a large contribution from the flatness constraint in the inversion and results in a long resolution length.

Figure 6

Statistical resolution matrices (a, c, e) and normalized matrices (b, d, f) of the problem presented in Fig. 1, obtained using different weights and the LSQR inversion method. The flatness weights were (a & b) 1, (c & d) 10 and (e & f) 20. The other properties are the same as those described in Fig. 5.

Open in new tabDownload slide

3.5 Non-linear case

For a small model space around a real model (or, e.g. for a given reference model m  ₀ sufficiently close to a real model, m) the non-linear forward equation, d=g(m), can be approximated using a linearization procedure,

(23)

where Δd=d−g(m  ₀) and Δm=m−m  ₀. Eq. (23) is similar to eq. (1). Similarly, the statistical resolution matrix was defined as

(24)

Using the above equation, an approximate statistical resolution matrix for the non-linear problem could easily be determined based on the approach described for estimating the statistical resolution matrix of a linear problem.

On the basis of eq. (24) we obtain:

(25)

In an inversion iteration for a linearized non-linear inversion, the inverse problem is, in fact, a linear inverse problem based on the relevant reference model. If the linear inversion during the iteration is perfect, the resolution matrix ℜ can be the identity matrix I or, say, ℜ=I. In this case we obtain:  

(26)

Eq. (24) is similar to eq. (13). However, the resolution matrices defined by these two equations have different levels of significance. The resolution matrix defined in eq. (24) does not represent the relation between the model and the solution like that in eq. (13) but it is an approximate relation between the model difference and solution difference for a given reference model. Therefore, we write ℜ in eq. (24) as ℜ(m  ₀). Eqs (23)–(26) indicate—but only when the given reference model is close enough to the true model and the linear inversion for the iteration is perfect—that the resolution matrix defined in eq. (24) can have a similar significance as that in eq. (13). However, the above-mentioned conditions are rarely true for non-linear inversion. Therefore, the resolution length estimated from ℜ(m  ₀) is often not our expected absolute model resolution length for an inversion, but an approximate resolution length of the model difference related to a given reference model. The real model resolution length will depend on both the reference model and the non-linearity of observation operator g of the inverse problem. In this case, the primary use of the resolution length retrieved from ℜ(m  ₀) is to measure the relative solution precision between parameters based on the degree to which the row of the matrix approximates a delta function, on the basis of the given reference model, as indicated by the example below.

1-D S-velocity inversions from surface wave dispersions are typical seismological non-linear inverse problems widely applied in crust and lithosphere studies. Although the hybrid resolution matrix may be applicable to linearized inversions (Herrmann & Ammon 2002), no publications have explicitly described the spatial resolution in related studies. Here, a synthetic example is presented for calculating the matrices and ℜ(m  ₀) for a 1-D inversion. The spatial resolution is also estimated: see Fig. 7. The program surf96 (Herrmann & Ammon 2002) was used to perform the forward calculation and linearized inverse iterations, and outputs the hybrid resolution matrix () at each inverse iteration. The vertical S-velocity model included 39 layers (each 5 km thick) and a half-space layer at the bottom (Fig. 7a), and damping in all inversions was 0.01. The value of at the 15th iteration is shown in Fig. 7(b). The last column of the hybrid resolution matrix is related to the bottom half-space layer, and each component in the column differs significantly from the other components in the row. I discarded all matrix components related to the bottom layer in the normalized hybrid resolution matrix (Fig. 7c). On the basis of the output model during the 15th iteration, I made 30 random synthetic models and their respective inverted solutions using the program, finally the statistical resolution matrix (Fig. 7d) was obtained by the method suggested earlier.

Figure 7

Resolution matrices of a Rayleigh-wave dispersion inversion for a 1-D S-velocity model. The synthetic observation includes the Rayleigh-wave group velocities over the periods from 5 to 185 s with intervals of 5 s. The forward calculation and the linearized inversion was performed using surf96 (Herrmann & Ammon 2002). The lines in (a) are the synthetic/true model (black) and the reference model or the output model at each linearized inversion iteration, indicated by the coloured lines, which vary from red to blue with the iteration number. The blue line indicates the outputted model at the 15th iteration and the red line indicates a constant velocity for all depths in the initial model. (b) The hybrid resolution matrix at the 15th iteration was the output of surf96, and its normalization is shown in (c). The matrix components of the half-space were discarded in (c). The resolution matrix in (d) is the normalized statistical resolution matrix based on the above-mentioned output model (blue line in panel a) at the 15th iteration.

Open in new tabDownload slide

As shown in the comparison in Figs 7(c) and (d), the spatial-resolution lengths from the matrices and ℜ(m  ₀) generally increased with the depth. For example, the resolution lengths of increased from 12 km (the spatial extent of a parameter is equal to the layer thickness, 5 km) for the 10th parameter to 25 km for the 20th parameter. The resolution lengths for ℜ(m  ₀) increased from 8 to 20 km. Non-negligible differences between the two resolution matrices were observed. Besides the difference in the resolution-length values of and ℜ(m  ₀), the resolution lengths of the layers at or near the model bottom (Fig. 7d) were very small (<10 km), which differed from the conditions of the hybrid resolution matrix (Fig. 7c). Because the statistical resolution matrix could be obtained directly using synthetic models and their respective solutions, the small resolution length at or near the bottom provided special information about the solutions to the random synthetic models. In fact, a high spatial resolution for a parameter may indicate that the parameter can be easily and well inverted. Fig. 7(a) shows that the parameters at or near the bottom layer could approach the synthetic model (black line) after only three iterations in this synthetic example, indicating a high resolution in the layer at or near the bottom. The small resolution length calculated based on the statistic resolution matrix should, therefore, be reasonable.

3.6 High-dimension or direction-dependent cases

For nd dimensions (nd > 1), a Gaussian-shaped function for an element r_i,_j at the ith row and jth column of ℜ can be written as

(27)

A resolution matrix for a specific direction can be expressed similar to the equation given. For a resolution matrix that does not depend on the dimension or direction, all w_i  ^k (k= 1, nd) are equal and can be written as w_i. Therefore, the statistical resolution matrix can be inverted using the same approach as was described for the 1-D case, above. Otherwise, a resolution matrix may be inverted along each dimension separately using the equation with one unknown. All w_i  ^k (k= 1, nd) also can be inverted simultaneously. Using the same approach described for eqs (19)–(22), we can obtain nd–dimensional equations according to

(28)

We can then obtain eq. (29), which includes nd unknowns, w_i  ^k (k= 1, nd). The approach used to solve eq. (22) was also used to obtain the unknowns by minimizing eq. (29),

(29)

4 Discussion

4.1 Required synthetic model number

The number of synthetic models determines the computational efficiency of the above approach. Each model and its respective solution are analogous to a measurement sample for resolution length, w_i. Consequently, the determination of the minimum number of synthetic models required to provide a stable inversion of w_i or a resolution matrix can be considered as a statistical problem of minimum sample size determination on the basis of a given confidence interval. Because a measurement is only taken for one parameter (w_i), the minimum sample number (ns) will not depend on the total number of parameters in the model but instead on the number of all possible values of w_i. As pointed out, with increasing ns, the mean effect of all samples can result in a stable w_i that is close to the true resolution length. Therefore, this problem can be taken as satisfying the central limit theorem.

The central limit theorem states that when the sample size is large, sample means will fall within ±Z_{α  /2} times the standard errors of the population mean, that is,

(30)

and

(31)

where E_i is called the margin of error, Z_{α  /2} is a constant for the relevant confidence level, and d_i is the standard deviation of all possible w_i. If the margin of error is known, the minimum sample size for the given confidence level can be obtained as:

(32)

For the 1-D example in Fig. 1, the real w_i can only be in the range of 0–20h. If all possible values of w_i have the same probability, the standard errors of the values, d_i, should be ∼5.8h. If the 95 per cent confidence level is used, Z_{α  /2}= 1.96. If we assume a margin of error (E_i) of 2.5 h, we obtain ns= 20. As described, a real test showed that 25 synthetic models (ns= 25) were sufficient to obtain a stable w_i.

The minimum sample number was also tested for by preparing a large inverse example. Here, a 2-D Rayleigh-wave group-velocity tomographic inversion was performed on the Antarctic. The inverted results using real observations of the Antarctic will not be discussed here, rather, the observation matrix will be used as an example for extracting the resolution-length information. The inverse horizontal 2-D model was parametrized using 12 163 equal-area hexagonal and pentagonal cells with widths of ∼1 great-circle degree (see Fig. 8a). The application of traditional checkerboard tests on the tetragonal cells was impossible due to the use of hexagonal and pentagonal cells. The inversion was a linear inverse problem, and LSQR was used here to obtain a solution. For this problem, the possible values of w_i are in the range of 0°–15°, and the standard errors of the values should be ∼4.35°; if the margin of error is ∼0.5°, the minimum sample number (ns) should be ∼290 at the 95 per cent confidence level according to eq. (32). Figs 8(b) and (c) present a random synthetic model and its inverted solution. Statistical resolution matrices were obtained using different ns pairs of random synthetic models and their respective solutions; the resolution lengths from the matrices using 20, 200 and 300 samples are shown in Figs 8(d)–(f). A comparison among Figs 8(c), (d) and (f) reveals that ns= 300, which is similar to the predicted minimum sample number above and sufficient to obtain a stable resolution length for the entire study region. The resolution lengths (Fig. 8f) are consistent with the path coverage (Fig. 8a). For example, within the continental Antarctic, the path density was high (Fig. 8a); the resolution length (Fig. 8f) was also high and approached the cell dimensions.

Figure 8

Resolution-length information for an ∼1°-wide cell 2-D Rayleigh-wave dispersion tomography set over a period of 50 s. Panel (a) shows the cells and the 5384 rays from source (circles) to station (triangles) used in this example. Most of the sources are earthquakes, and a few are seismic stations for ambient-noise cross-correlation rays. Panel (b) shows one random synthetic 2-D model, (c) shows the inverted solution to the synthetic model in (b); (d–f) show the spatial-resolution maps from 20, 100 and 300 pairs of random synthetic models and their solutions.

Open in new tabDownload slide

In the 1-D tests using 100 parameters, as described, 25 synthetic models (ns= 25) were sufficient to obtain a stable w_i. Here, a 2-D example with 12 163 parameters required only 300 models. The tests showed that fewer than several hundred pairs of synthetic models and solutions produced a stable w_i, and the minimum number of synthetic models depends on the possible values of w_i but not on the number of parameters. The minimum sample numbers in the real tests are similar to the prediction by the central limit theorem.

4.2 Does the inversion of a random synthetic model simulate a checkerboard test?

Interestingly, the solution pattern (Fig. 8c) for a random synthetic model (Fig. 8b) submitted to the tests described based on the 2-D Rayleigh-wave group-velocity tomography of the Antarctic resembled a checkerboard test output. The inverted solution for any random synthetic model was found to display notable patterns similar to the output of a traditional checkerboard test using an input model with checker pattern anomalies. The anomalous pattern size in each solution (Fig. 8c) was comparable to the resolution-length scale estimated from the resolution matrix analyses, as shown in Fig. 8(f).

The checker-like patterns in the solution based on a random synthetic model are also shown using ∼4°-wide equally sized cells (see Fig. 9). Tests (not shown here) on a linearized inversion for a 3-D random synthetic model also yielded a checker pattern solution, similar to that shown in Figs 8(c) and 9(c). These tests demonstrated that even though the inverted solutions should depend on the respective random input models, the inverse output solution of any random synthetic input model directly yielded a specific checker pattern that can give some indicative resolution-length information.

Figure 9

Resolution-length information for an ∼4°-width cell 2-D Rayleigh-wave dispersion tomography set for a period of 50 s. Subpart (a) shows the cells and the 5393 rays used in this example. The rays are nearly the same as those in Fig. 8(a). Subpart (b) shows one random synthetic 2-D model, (c) shows the inverted solution to the synthetic model (b); (d) shows the resolution-length map for 300 pairs of random synthetic models and their solutions.

Open in new tabDownload slide

The resolution length corresponds to the size of the smallest possible feature that can be detected. For a given resolution matrix, the resolution length at the position of a parameter is consistent with the number of neighbouring parameters, which make an obvious contribution to the solution parameter; see eq. (14) and Fig. 2(b). The inverted solutions (e.g. Figs 8c or 9c) revealed that the resolution length at the position of a parameter (Figs 8f or 9d) is consistent with the number of neighbouring parameters which have a value similar to that of the parameter in the inverted solution; in contrast, the parameters in the synthetic/real model (Figs 8b or 9b) are not similar to one another. Put simply, the resolution-length estimation can be considered as an analysis of how many neighbouring parameters have similar values to each other around a given position. According to this idea, the resolution length could be estimated directly via visualization of the solution inverted from a random synthetic model. Following this line of reasoning, variation in the model/solution parametrization (e.g. absolute value of the physical quantity, anomaly and relative anomaly) may influence the inversion for a statistical resolution matrix especially for a linearized non-linear inverse problem. It seems that the anomaly or relative anomaly could be superior in the model/solution parametrization when retrieving resolution-length information using the approach presented earlier.

Even though the information obtained from the inverted solution of a random synthetic model depends on said synthetic model, the visualization of the solution inverted from a random synthetic model can not only give some indicative resolution-length information (as explained), but it can also yield the direction dependence of the resolution. Because most of the observation stations (Fig. 8a) are located inside continental Antarctica and earthquakes originate at the plate boundaries, the rays have a better defined cross-section at a position far from the plate boundary, while rays at positions close to the plate boundary are oriented much more parallel to their neighbouring rays. In practice, the resolution length along the parallel direction (nearly parallel to the meridians) of the rays should be longer than that in the normal direction (nearly parallel to the latitude). In the solution (Fig. 8c) inverted from a random synthetic model, most of the anomaly pattern looks like a strip with a long axis along the meridian and a short axis along a line of latitude, especially for the oceanic region where the patterns of the anomaly in the input synthetic model (Fig. 8b) exhibit little similarity to the strip anomaly in the solution. According to the reasoning above, the strip anomaly pattern in the oceanic region should indicate that the resolution length along the meridian is longer than that along the line of latitude in the oceanic region, which is consistent with indications from the ray distribution.

4.3 Resolution lengths of discontinuities between neighbouring parameters

Discontinuities or sharp variations between neighbouring model parameters (e.g. m_i and m_j) provide important information in geophysical studies. Therefore, the resolution lengths of a discontinuity are needed for an inversion, if possible. Using eq. (14), an equation was obtained describing the differences between neighbouring parameters,

(33)

Eq. (33) indicates that the differences between two solution parameters (m  _i and m  _j) depend not only on the entire model, but also on the resolution matrix corresponding to the ith and jth row elements.

For a perfect inversion, each row of the resolution matrix is a delta function, and eq. (33) can be written in terms of eq. (34),

(34)

Eq. (34), corresponding to a perfect inversion, can only indicate the differences between two neighbouring parameters that can be well solved; however, although the width of a discontinuity (such as a Moho) may be infinitely small or thin, the discontinuity cannot be localized to the border between two parameters. In fact, it can only be concluded that the resolution length of a discontinuity does not exceed the sum of the half spatial extents of two parameters, and the real resolution length of a discontinuity can be much shorter. Fig. 2(a) provides an example showing that the discontinuity resolution length at some depth (e.g. 30 or 60 km) can be much higher than the resolution length of a model parameter if the discontinuity is located on a cell's boundary. A simple way to obtain a reasonable resolution length of a discontinuity in a perfect inversion is to decrease the parameter spatial extent and re-calculate the resolution matrix until the matrix rows associated with the discontinuity are not delta function. The parameter spatial extent can then be taken as the largest bound on the discontinuity resolution length.

4.4 Computation requirements and efficiency

Our approach to calculating the statistical resolution matrix can be divided into two steps: the first step is to create tens or hundreds of random synthetic models and solutions; the second step is to directly invert for the resolution lengths using these synthetic models and solutions. The first step is similar to the procedure used in traditional checkerboard tests; however, the input model is easier to create because it is a random model without a regular anomaly pattern distribution. The forward calculation and the solution inversion procedure applied in the first step are basic calculations for an inversion study. The only coding burden for the skill presented here is to write a short program that inverts the synthetic models and their respective solutions to identify the resolution lengths using a grid search. This coding effort is trivial because it is a one-parameter inversion, and the grid search is the simplest of the inversion methods.

The second step, involving a grid search for the resolution-length scale, does not require much memory. The primary memory consumption involves retaining tens or hundreds of pairs of models and their solutions in memory for use in the grid-search inversion. The maximal memory required for the computation approaches 2 ×ns×m, where ns is generally not more than several hundreds and m is the model parametrization size. The advantages of the approach described here are that the maximal memory consumption is much smaller than that required from a calculation of the resolution matrix using matrix operations. The idea presented here is especially suitable for very large inverse problems. If the information of the models and solutions are stored on a disk, the memory burden will approach zero, but the computation time will be much longer due to the need for frequent disk reading actions.

The skill presented here is isolated from the specific forward/inversion procedure and does not directly use any relation between observation and model parameters such as observation matrix/operators. Therefore, it can be applied to general inverse problems. For the same reason, one should not be surprised that the skill may have a worse computational efficiency than other skills using the relation between observation and model parameters of the relevant inverse problem. For example, for a small-scale inversion using SVD, the direct/regularized/hybrid resolution matrices can be obtained easily using decomposition matrices from the observation matrix, but calculation of a statistical resolution matrix needs much more time, because ns forward/inversion computations and m one–unknown grid searches for resolution lengths are needed. On the other hand, since ns is often much smaller than m for large inverse problems, the efficiency of statistical resolution matrix computation is better than that of brute-force calculations of Backus–Gilbert kernels, because the latter need to solve the system m times, and the computation time to solve the system normally is much longer than a one–unknown grid search for the relevant resolution length.

Lanczos iteration inversions, such as LSQR, can be used to quickly obtain the solution and they have been widely applied to large linear inversions. After some modifications to the inversion codes, a resolution matrix can be obtained (e.g. Yao 1999; Zhang & Thurber 2007). Because inversion programs only provide limited options to include a priori constraints, if one wants to use one's preferred constraints, the simplest way is to take the regularized vector b and matrix A, as expressed in eq. (6), as the input observation vector and observation matrix, respectively. In this case, the output resolution matrix is a regularized resolution matrix which cannot give any resolution-length information, as shown earlier. Mostly, researchers engaged in tomographic studies have to perform checkerboard tests, which can give indicative resolution-length information. The skill presented here can be applied directly to these inversions to obtain resolution lengths. Lanczos iteration inversions are efficient, which will ensure the efficiency of the statistical resolution matrix computation. For the reasons outlined, statistical resolution matrix computations may be particularly suitable for applications of Lanczos iteration inversions.

A discussion to improve the efficiency of the first step of the forward/inversion computations is beyond the scope of this paper. As for the second step, the efficiency of grid searches for resolution lengths can be easily improved for modern computers equipped with processors containing more than a single core. The one-parameter grid-search inversion for the relevant resolution length must be repeated for all model parameters. Modern parallel computational approaches (e.g. OpenMP: see http://www.openmp.org) suitable for repeating calculations, can be easily applied to these codes, which decreases the computation time significantly.

In general, the computation efficiency of the skill presented here may be worse than that of some other skills. However, the skill discussed here can be applied to general problems, and needs little coding work and computation memory. The computation efficiency can be improved using modern parallel computing approaches. Statistical resolution matrix computation may be particularly suitable for application of Lanczos iteration inversions, because the application can quickly obtain a solution although it is difficult to obtain resolution lengths.

5 Conclusions

Applications of the Backus–Gilbert inversion method are few; however, the basis of the inversion (resolution matrices and kernels that connect real models to calculated solutions) are widely used to evaluate inversion perfectness and to analyse the spatial resolution of a linear inversion. The resolution matrices described in previous publications fall into three classes: (1) direct resolution matrices, which are the product of an observation matrix and its inverse; (2) regularized resolution matrices, which are the product of a regularized observation matrix and its inverse; (3) and hybrid resolution matrices, which are the product of an observation matrix and the inverse of the regularized observation matrix. The first class of matrices provides resolution-length information and the perfectness of the inversion but is seldom presented in studies because the inversion is often ill-posed and must be regularized. The second class matrices can be used to evaluate the inversion perfectness under a given mathematical system, but these matrices cannot provide any information about the resolution length. The third class of resolution matrices can be used to extract resolution-length information and to evaluate the perfectness of an inversion; however, these matrices depend on the contributions of regularization to the inversion. The simultaneous presentation of all three resolution matrices in a study could fully articulate the inversion perfectness and resolution-length information. Unfortunately, it is often impossible to perform any of these calculations in the context of very large inverse problems. Here, a new class of resolution matrices, called the statistical resolution matrices, is suggested to permit the direct inversion from synthetic models. These matrices are especially suitable for very large linear/linearized inverse problems.

A resolution matrix defines a linear projection in which each solution parameter is derived from the weighted average of the neighbouring solution parameters and the resolution matrix elements are the weights. Ideally, a parameter and its neighbours provide a large contribution/weight to the average summation for a parameter in the solution; therefore, all elements in each row of the resolution matrix are assumed to be Gaussian functions over the parameter distances. The resolution matrix exclusively provides the properties of the observation matrix but does not depend on any particular model or observation values; therefore, random synthetic models and their solutions may be used, and the projection equation becomes a non-linear inversion of only one parameter, the width of the Gaussian shape. The projection can then be directly inverted using a simple grid-search method. The inverted widths directly indicate the resolution length.

The tests described here showed that the statistical resolution matrix not only measures the resolution obtainable from the data, it also provides reasonable spatial/temporal resolution or resolution-length information. Estimates of the statistical resolution matrix can be used for both direction-dependent and direct-independent resolution estimates. The estimates depend on the real forward/inversion processes and are independent of the approach used to determine the solution inversion. Therefore, it is not necessary to modify the preferred forward/inversion codes/methods if a user wishes to obtain the resolution length. Matrix operations are not needed in the estimation, indicating that estimating the statistical resolution matrix is especially suitable for very large linear/linearized inverse problems.

Interestingly, the tests showed that for a random synthetic input model without a specific checker pattern, the inverse output solution directly provided a checkerboard anomaly pattern that indicatively provided information on not only the resolution length but also on the direction dependence of the resolution.

All the codes to calculate statistic resolution matrix with examples for the inverse problem in Fig. 1 are available through: http://www.seismolab.org/people/meijian/.

Acknowledgments

This work was financed by grants 40874021 and 41174039 (NSFC-China), and by CHINARE2012-02-05. I would particularly thank the editor, Jeannot Trampert, and the two reviewers, Andreas Fichtner and Jonathan MacCarthy, for their thoughtful and constructive comments and suggestions. All figures are generated using the Generic Mapping Tools (GMT; Wessel & Smith 1991).

References