Resolution analysis of geophysical images: Comparison between point spread function and region of data influence measures

Practical decisions are often made based on the subsurface images obtained by inverting geophysical data. Therefore it is important to understand the resolution of the image, which is a function of several factors, including the underlying geophysical experiment, noise in the data, prior information and the ability to model the physics appropriately. An important step towards interpreting the image is to quantify how much of the solution is required to satisfy the data observations and how much exists solely due to the prior information used to stabilize the solution. A procedure to identify the regions that are not constrained by the data would help when interpreting the image. For linear inverse problems this procedure is well established, but for nonlinear problems the procedure is more complicated. In this paper we compare two different approaches to resolution analysis of geophysical images: the region of data influence index and a resolution spread computed using point spread functions. The region of data influence method is a fully non-linear approach, while the point spread function analysis is a linearized approach. An approximate relationship between the region of data influence and the resolution matrix is derived, which suggests that the region of data influence is connected with the rows of the resolution matrix. The point-spread-function spread measure is connected with the columns of the resolution matrix, and therefore the point-spread-function spread and the region of data influence are fundamentally different resolution measures. From a practical point of view, if two different approaches indicate similar interpretations on post-inversion images, the confidence in the interpretation is enhanced. We demonstrate the use of the two approaches on a linear synthetic example and a non-linear synthetic example, and apply them to a non-linear electromagnetic field data example. I N T R O D U C T I O N Geophysical inversion is inherently ill-posed and non-unique. To solve ill-posed geophysical inverse problems, we routinely incorporate prior information into our solution. This stabilizes the solution, but also complicates the appraisal process because of the bias introduced into the solution. If the end goal of our data analysis involves decisions based on our solution, it is of great importance that we take steps to quantify how ∗ E-mail: camiller@cgiss.boisestate.edu much information came from data and how much came from prior information. For different applications it is difficult to produce a unified interpretation goal; however, an important question is to determine which regions of the image are influenced by a priori information and which regions are influenced by the data. Although different ways of incorporating a priori information are sometimes considered to be subjective (Friedel 2003), regularization an essential step in obtaining a stable solution that is physically meaningful. Therefore, we usually cannot neglect imposing a priori information through the regularization procedure in the inverse problem. However, if we can C © 2007 European Association of Geoscientists & Engineers 835 836 C.R. Miller and P.S. Routh develop some procedure to separate out the regions that are not constrained by the data, then it would greatly simplify the interpretation and provide a robust way of connecting the interpretation with the physical experiment that generated the data. Oldenburg and Li (1999) used two reference models in a non-linear DC resistivity inversion to determine the regions that are honoured by the data and generate a filter function called depth-of-investigation curves to interpret the inverted image. Alumbaugh and Newman (2000) used 50% spread criteria from the point spread functions to highlight the regions of the model that are better resolved. In seismic tomography, some commonly applied approaches with a similar goal are the impulse test and the checkerboard test. In these approaches, the raypaths of the real data are used on a synthetic impulse model or checkerboard model to generate the data. Subsequently, the inverted images are compared with the synthetic model in order to understand the model resolution (Shearer 1999). In the depth-of-investigation approach of Oldenburg and Li (1999), the index is computed using the models obtained from two inversion runs with different reference models. The image region that is most influenced by the data has a low depth-of-investigation index and we refer to this region as the region of data influence. The advantage of this approach is that it can be applied to non-linear problems without any linearization assumption. Alternatively, the resolving capability of the geophysical data can be obtained by computing columns of the resolution matrix, called point spread functions. We examine spread measures computed using point spread functions and compare them with the region of data influence index. Both of these methods have previously been used in the context of appraisal analysis (Oldenburg and Li 1999; Alumbaugh and Newman 2000). Oldenborger (2006) applied the region of data influence approach to the 3D electrical-resistivity-tomography experiment and showed that it is a better approach to examining resolution, compared to sensitivity measures. In this paper, we show the similarity of the two approaches. From a practical point of view, if two different approaches indicate similar interpretations on post-inversion images, the confidence in the interpretation is enhanced. The examples presented in this paper are one-dimensional (1D). We recognize that in reality the earth is three-dimensional (3D), but the generic procedure of the approaches presented here will be valid for problems in any dimensions. For computational ease, we use 1D examples throughout the paper. The paper is presented as follows: Using a linear problem we examine the relationship and the mathematical connections that exist between the two appraisal tools. We illustrate the relationship with a linear synthetic example and a non-linear synthetic example, and we apply it to a non-linear field example using controlled source electromagnetics. R E S O L U T I O N M E A S U R E S Linearized resolution analysis Insight into the resolution analysis problem can be gained by considering a linear inverse problem, dobs = Gm + ε. The linear problem has been studied in detail in previous literature (Backus and Gilbert 1968, 1970; Oldenburg 1983; Snieder 1990; Parker 1994; Alumbaugh and Newman 2000) and here we present the necessary equations to illustrate the various components of the resolution measure. In a linear problem, the objective function to be minimized is φ = ∥∥∥Wd(dobs − Gm)∥∥∥2 + β‖Wm(m − m0)‖, (1) where G ∈RN×Mis forward modelling operator, Wd ∈ RN×N is the data weighting matrix, Wm ∈ RM×M is the model weighting matrix used to impose constraints in the solution such as smoothing, m ∈ RM×1 is the model vector, m0 ∈ RM×1 is the reference model vector, dobs ∈ RN×1 is the data vector and β is the regularization parameter that controls the relative contribution of the data misfit and the model norm. Denoting the Hessian by H = (GWd WdG + βWmWm), the recovered model is given by m̂ = HGWd WdG m + HGWd Wd ε + HWmWm m0. (2) In a more compact notation we can represent the estimated model in equation (2) by m̂ = R m + η + γ , where R is the resolution matrix, R = (GTWTd WdG + βWTmWm)−1 GWd WdG. (3) R provides a linear mapping between the estimated model parameters, m̂ and the true model parameters, m. For a nonlinear problem given by dobs = F(m) + ε, the estimated solution for the model perturbation is given by δm̂ = HGWd WdG δm + HGWd Wd(Q(δm) − P(δm̂) + ε) + HWmWm δm0, (4) where G is the sensitivity matrix, δm = m(k+1) − m(k) is the model perturbation, δm0 = m0 − m(k) is difference between the reference model and the model estimate at the kth iteration and Q(δm), P(δm̂)are the higher-order terms, typically neglected in the linearized analysis. When the higher-order C © 2007 European Association of Geoscientists & Engineers, Geophysical Prospecting, 55, 835–852 Resolution analysis of geophysical images 837 terms are not negligible, bias is introduced into the solution of the estimated model, in addition to the bias introduced by the reference model. For our analysis we will neglect the higherorder terms and carry out the interpretation using linearized analysis, using R in equation (3) as resolution mapping. Point spread functions The columns of the resolution matrix are called the point spread functions. A point spread function describes how a delta-like perturbation in the model manifests itself in the inverted result. If a particular point spread function is wide and/or has large side lobes, then the corresponding model is poorly resolved at that location and the resolution ‘width’ is broader than the cell dimension. This introduces the idea that spread of these point spread functions can be used as a measure to assess the resolving capability of the data and the prior information. We note that the prior information can effectively be treated as data in the inverse problem, so the point spread function and hence the spread will be influenced by the prior information. Spread criterion for the point spread function To compare the resolution capability of the data and prior information in different regions of the model, point spread functions can be compared directly with each other. This can be useful especially in large-scale problems where the question about resolvability is limited to few specific regions of the model. Although the point spread function comparison is a valid procedure and has the capability to convey a significant amount of detail, the process can be tedious when we are faced with comparing thousands of cells throughout the entire model instead of a few cells within a specific region of interest. Here, we develop a spread criterion to summarize the information con