## Abstract

We investigate the practice of regularization (also termed damping) in inverse problems, meaning the use of prior information to supplement observations, in order to suppress instabilities in the solution caused by noisy and incomplete data. Our focus is on forms of regularization that create smooth solutions, for smoothness is often considered a desirable—or at least acceptable—attribute of inverse theory solutions (and especially tomographic images). We consider the general inverse problem, in its continuum limit. By deconstruction into the part controlled by the regularization and the part controlled by the data kernel, we show the general solution depends on a smoothed version of the back-projected data as well as a smoothed version of the generalized inverse. Crucially, the smoothing function that controls both is the solution to the simple data smoothing problem. We then consider how the choice of regularization shapes the smoothing function, in particular exploring the dichotomy between expressing prior information either as a constraint equation (such as a spatial derivative of the solution being small) or as a covariance matrix (such as spatial correlation falling off at a specified rate). By analyzing the data smoothing problem in its continuum limit, we derive analytic solutions for different choices of regularization. We consider four separate cases: (1) the first derivative of the solution is close to zero, (2) the prior covariance is a two-sided declining exponential, (3) the second derivative of the solution is close to zero, and (4) the solution is close to its localized average. First-derivative regularization is put forward as having several attractive properties and few, if any, drawbacks.

This is a preview of subscription content, access via your institution.

## References

Abers, G. (1994).

*Three-dimensional inversion of regional P and S arrival times in the East Aleutians and sources of subduction zone gravity highs*, J. Geophys. Res.*99*, 4395–4412.Aki, K., A. Christoffersson and E. Husebye (1976).

*Three-dimensional seismic structure under the Montana LASA*, Bull. Seismol. Soc. Am.*66*, 501–524.Backus G.E. and Gilbert, J.F. (1968).

*The resolving power of gross earth data*, Geophys. J. R. Astron. Soc.*16*, 169–205.Backus G.E. and Gilbert, J.F. (1970).

*Uniqueness in the inversion of gross Earth data*, Phil. Trans. R. Soc. London Ser. A*266*, 123–192.Boschi, L. and A.M. Dziewonski (1999),

*High- and low-resolution images of the earth’s mantle: implications of different approaches to tomographic modeling*, J. Geophys. Res.*104*, 25567–25594.Ekstrom, G., J. Tromp and E.W.F. Larson (1997).

*Measurements and global models of surface wave propagation*, J. Geophys. Res.*102*, 8127–8157.Gradshteyn, I.S. and Ryzhik, I.M., Tables of Integrals, Series and Products, Corrected and Enlarged Edition (Academic, New York, 1980).

Hetenyi, M., (1979). Beams on Elastic Foundation (University of Michigan Press, Ann Arbor, 1979).

Humphreys, E., R.W. Clayton and B.H. Hager (1984).

*A tomographic image of mantle structure beneath southern California*, Geophys. Res. Lett.*11*, 625–627.Laske, G. and G. Masters (1996).

*Constraints on global phase velocity maps from long-period polarization data*, J. Geophys. Res.*101*, 16059–16075.Lawson, C. and Hanson, R., Solving Least Squares Problems (Prentice-Hall, New York, 1974).

Levenberg, K., (1944).

*A method for the solution of certain non-linear problems in least-squares*, Q. Appl. Math.*2*, 164–168.Menke, W., Geophysical Data Analysis: Discrete Inverse Theory (First Edition). (Academic, New York, 1984).

Menke, W. (2005). Case studies of seismic tomography and earthquake location in a regional context, in Seismic Earth: Array Analysis of Broadband Seismograms, A. Levander and G. Nolet, eds., Geophysical Monograph Series 157. American Geophysical Union, 7–36.

Menke, W., Geophysical Data Analysis: Discrete Inverse Theory (MATLAB Edition), (Elsevier, New York, 2012).

Menke, W. (2014). Resolution and Covariance in Generalized Least Squares Inversion, in press in Surveys of Geophysics.

Menke, W. and Menke, J., Environmental Data Analysis with MATLAB (Elsevier, New York, 2011).

Menke W. and Abbott, D., Geophysical Theory (Columbia University Press, New York, 1989).

Nettles, M., and A.M. Dziewonski (2008).

*Radially anisotropic shear-velocity structure of the upper mantle globally and beneath North America*, J. Geophys. Res.,*113*, doi:10.1029/2006JB004819, 2008.Smith, W. and Wessel, P. (1990),

*Gridding with continuous curvature splines in tension*, Geophysics*55*, 293–305.Trampert, J. and J.H. Woodhouse (1995),

*Global phase velocity maps of Love and Rayleigh waves between 40 and 130 seconds*, Geophys. J. Int.,*122*, 675–690.Tarantola A. and Valette B. (1982a),

*Generalized non-linear inverse problems solved using the least squares criterion*, Rev. Geophys. Space Phys.*20*, 219–232.Tarantola A, Valette B. (1982b),

*Inverse problems**=**quest for information*, J. Geophys.*50*, 159–170.Tromp, J., Tape, C. and Liu, Q. (2005),

*Seismic tomography, adjoint methods, time reversal and banana-doughnut kernels*. Geophys. J. Int.,*160*: 195–216. doi:10.1111/j.1365-246X.2004.02453.x.Wiggins, R.A. (1972),

*The general linear inverse problem: Implication of surface waves and free oscillations for Earth structure*, Rev. Geophys. Space Phys.*10*, 251–285.Wikipedia (2014), Screened Poisson Equation, en.wikipedia.org.

Zha, Y., S.C. Webb, S.S. Wei, D.A. Wiens, D.K. Blackman, W.H. Menke, R.A. Dunn, J.A. Conders (2014),

*Upper mantle shear velocity structure beneath the Eastern Lau Spreading Center from OBS ambient noise tomography*, Earth Planet. Sci. Lett.*408*, 194–206.

## Acknowledgments

This research was supported by the US National Science Foundation under grants OCE-0426369 and EAR 11-47742.

## Author information

### Affiliations

### Corresponding author

## Appendix: Derivations of Smoothing Kernels and Covariances for Case Studies

### Appendix: Derivations of Smoothing Kernels and Covariances for Case Studies

### Case 1: First-Derivative Minimization

The operator \({ \mathcal{L}} = \varepsilon {\text{d}}/{\text{d}}x\) has translational invariance, so we expect that the smoothing kernel *a*(*x*, *x*′) = *a*(*x* − *x*′) will depend only upon the separation distance (*x* − *x*′) (as also will *C*
_{
h
}, *P*
_{
h
}, *Q*
_{
h
}, and *R*). Without loss of generality, we can set *x*′ = 0, so that Eq. (13) becomes

Here, we utilize the relationship that \(\left( {{\text{d}}/{\text{d}}x} \right)^{\dag } = - {\text{d}}/{\text{d}}x\). The solution to this well-known 1D screened Poisson equation is

This solution can be verified by substituting it into the differential equation:

Here, we have relied on the fact that (d/d*x*)|*x*| = sgn(*x*) and \(\left( {{\text{d}}/{\text{d}}x} \right) {\text{sgn}}(x) = 2\delta \left( x \right)\). Note that *a*(*x*) is a two-sided declining exponential with unit area and decay rate *ɛ*
^{−1}. Because of the translational invariance, the integral in Eq. (11) has the interpretation of a convolution. The solution is the observed data *d*(*x*) convolved with this smoothing kernel:

The variance *C*
_{
h
} of the prior information satisfies Eq. (15a):

This is a 1D Poisson equation, with solution

This solution can be verified by substituting it into the differential equation:

The covariance *C*
_{
h
}(*x* − *x*′) implies that the errors associated with neighboring points of the prior information equation *m*(*x*) = 0 are highly and positively correlated, and that the degree of correlation declines with separation distance, becoming negative at large separation.

Finally, we note that the operator \({\mathcal{L}} = \varepsilon {\text{d}}/{\text{d}}x\) is not self-adjoint, so that it is not the continuous analogue of the symmetric matrix **C**
^{−1/2}_{
h
}
. As described earlier, we can construct a symmetric operator by introducing a unitary transformation. \({\mathcal{L}}\) is antisymmetric in *x*, but we seek a symmetric operator, so the correct transformation is the Hilbert transform, \({\mathcal{H}}\), that is, the linear operator that phase-shifts a function by *π*/2. It obeys the rules \({\mathcal{H}}^{\dag } = - {\mathcal{H}}\), \({\mathcal{H}}^{\dag } {\mathcal{H}} = 1\), and \({\mathcal{H}}\left( {{\text{d}}/{\text{d}}x} \right) = \left( {{\text{d}}/{\text{d}}x} \right){\mathcal{H}}\). The modified operator \({\mathcal{L}}_{\text{sa}} = \varepsilon { \mathcal{H}}{\text{d}}/{\text{d}}x\)
*is* self-adjoint and satisfies \({\mathcal{L}}_{\text{sa}}^{\dag } {\mathcal{L}}_{\text{sa}} = {\mathcal{L}}^{\dag } {\mathcal{L}}.\)

### Case 2: Exponentially Decaying Covariance

For a covariance described by a two-sided declining exponential function,

By comparing Eqs. (42) and (43), we find that this prior covariance is the inverse of the operator

The smoothing kernel solves the equation

By analogy to Eqs. (42) and (43), the smoothing kernel is

An operator \({\mathcal{L}}\) that reproduces the form of \({\mathcal{L}}^{\dag } {\mathcal{L}}\) given in Eq. (50) is

The function *P*
_{
h
} solves Eq. (15b), \({\mathcal{L}}^{\dag } P_{h} = \delta (x)\), which for the operator in (30) has the form of a one-sided exponential,

Here, *H*(*x*) is the Heaviside step function. Because of the translational invariance, the inner product in Eq. (14) relating *P*
_{
h
} to *C*
_{
h
} is a convolution. That, together with the rule that the adjoint of a convolution is the convolution backwards in time, implies that \(C_{h} \left( t \right) = P_{h} \left( { - t} \right)*P_{h} \left( t \right) = P_{h} \left( t \right){ \star }P_{h} (t)\), where \({ \star }\) signifies cross-correlation. The reader may easily verify that the autocorrelation of Eq. (31) reproduces the formula for *C*
_{
h
} given in (25). Unfortunately, its Hilbert transform cannot be written as a closed-form expression, so no simple formula for the symmetrized form of *P*
_{
h
}, analogous to **C**
^{1/2}_{h}
, can be given.

### Case 3: Second-Derivative Minimization

The smoothing kernel *a*(*x*) satisfies the differential equation

This well-known differential equation has solution (Hetenyi 1979; see also Menke and Abbott 1989; Smith and Wessel 1990; Menke 2014)

The variance *C*
_{
h
} of the prior information satisfies Eq. (15a):

and by analogy to (47) has solution

This solution can be verified by substituting it into the differential equation:

This function implies a steep drop off in covariance between neighboring points and increasingly great anticorrelation with distance.

### Case 4: Damping towards Localized Average

From Eq. (37), we find that the smoothing kernel *a*(*x*) satisfies

We now make use of the fact that the operator \({\mathcal{L}}_{s} = 1 - \eta^{ - 2} \text{d}^{2} /\text{d}x^{2}\) is the inverse to convolution by *s*(*x*). Applying \({\mathcal{L}}_{s}\) twice to (37) yields the differential equation

We solve this equation by finding its Green function [that is, solving (39) with *f*(*x*) = *δ*(*x*)] and then by convolving this Green function by the actual *f*(*x*). This Green function can be found using Fourier transforms, with the relevant integral given by Equation 3.728.1 of Gradshteyn and Ryzhik (1980) (which needs to be corrected by dividing their stated result by a factor of 2). The result is

where

We determine the area under the smoothing kernel by taking the Fourier transform of (61):

and evaluating it at zero wavenumber. Thus, *a*(*k* = 0) = 1; that is, the area is unity.

## Rights and permissions

## About this article

### Cite this article

Menke, W., Eilon, Z. Relationship Between Data Smoothing and the Regularization of Inverse Problems.
*Pure Appl. Geophys.* **172, **2711–2726 (2015). https://doi.org/10.1007/s00024-015-1059-0

Received:

Revised:

Accepted:

Published:

Issue Date:

### Keywords

- Inverse theory
- tomography
- spatial analysis
- damping
- smoothing
- regularization