Inverse Theory (Geophysics 529/Math 519) Class Homepage

Richard Aster
Geophysics Program
Department of Earth and Environmental Science

Brian Borchers and Oleg Makhnin
Mathematics Department

New Mexico Institute of Mining and Technology
Socorro, New Mexico

Cliff Thurber
Department of Geology and Geophysics
University of Wisconsin-Madison

2010 Class Information

Spring 2010 Syllabus and General Class Guidelines

Ancillary Reading List (pdf)

Matlab Primer, c/o Kermit Sigmon (pdf)

Example Past Class Projects (2009)

Example Past Class Projects (2007)

Example Past Class Projects (2005)

Example Past Class Projects (2004)

Class Bulkmail

Parameter Estimation and Inverse Problems Textbook Homepage

2010 Class Projects

Image Processing

Dapo Oduba

An image is a matrix from a real vector space R(m x n x d) where m= the number of rows, n= number of columns and d=1, the dimension of a color space. This project models the image deblurring process characterized by a blurring matrix and an observed image. Methods of recovering the image are considered while avoiding noise amplification. An effective regularization technique is examined which involves comparing methods such as the least square method and the Tikhonov method. Also, positivity is restricted for the purpose of additional regularization benefits.

Mean Free Path Length and Cosmic Ray Flux

Michelle Hinz

 A parameter estimation of the absorption mean free path length, Λ, for the cosmic-ray flux will be recalculated using a multiplicity modified neutron flux network. Λ increases with increasing altitude as does the uncertainty associated with generally used values. Depending on geomagnetic latitude and altitude, Λ values range from 121 g/cm2 to 149 g/cm2. Dunai (2000) has reduced the overall uncertainty of calculations based on Λ by using Λ values which are dependent on paleo-inclination, I. Also, along with other considerations, he normalized neutron flux data so neutron monitor counting rates would be comparable. Despite these improvements to the scaling scheme, cosmogenic nuclide production rates using scaling factors published by Dunai (2000) can still have overall uncertainties of 20% for high altitudes (e.g. 7000m). Most recently, neutron monitor results are indicative of neutron multiplicity (a higher energy neutron counted as multiple neutrons rather than a singular neutron) and modeling efforts are addressing this multiplicity effect.  After a neutron flux network with multiplicity corrections have been made and following the approach used in Dunai 2000, Λ(I) can be calculated with non linear regression and a five parameter sigmoidal function.   

Lightning Location Using LEFA Data

Jeff Lapierre

Many people think that lightning is simply the flash they see. In reality this is just one of many processes which occur during the lifetime of a lightning discharge. The main processes of a flash consist of leaders and return strokes with smaller components before, after, and in between these two. The initial leader is charges propagating from the origin of the flash while breaking down and ionizing the air. If the ash is a cloud-to-ground (CG) ash this leader will propagate downward. The leader propagates downward in steps rather than in a smooth motion, thus earning the name stepped leader. Once the leader reaches the ground, a 'short-circuit' is formed and a current can flow to negate the charge on the channel created by the leader. This is what causes the bright ash and thunder. The Langmuir Electric Field Array (LEFA) is an 8 instrument array that measures electric field changes. Using highly time resolved data (50,000 sample/s) the propagation of the leader as it approaches ground (for CGs) can be approximated. This model assumes that an equal amount of charge is being deposited onto the channel for each step that the leader takes until it reaches the ground and a return stroke occurs. Error estimations will be obtained using Monte Carlo techniques. Because the model is non-linear, the error propagation techniques for linear systems will not be sufficient. By introducing random error into the data and running the model many times the sensitivity of the model can be found. Finally, results obtained by the Lightning Mapping Array (LMA) can be compared to the model results. The LMA uses the emission of radio waves during electrical breakdown of air to locate leaders.

The Metropolis Algorithm: A Parameter Search Method for Models of Arbitrarily Complexity

Michael Herman

In many situations, models designed to simulate complicated physical behavior reach a level of complexity such that popular inverse methods cannot be used to determine a desirable set of model parameters. The Navier-Stokes equations, for example, are a set of three non-linear equations describing fluid dynamics. Coupled with expressions describing the continuity of mass and conservation of energy, the resulting system of equations is often used to model planetary atmospheres and in practice becomes an intricate set of iterated computer instructions. In this way, sufficiently complicated models can appear like a black box to the user and even to the model’s author when the need for parameter tuning arises. The Metropolis algorithm is presented as a tool apropos to this circumstance. The algorithm is described and it is shown how this flexible, simple numerical method can be used to locate an optimal parameter set for a model of arbitrary complexity based on the slow approach to minimizing a cost function. An interesting analogy with the Metropolis algorithm, described by Kirkpatrick et al 1983, related the minimization process to statistical mechanics. The benefits of this slow annealing approach presented in terms of Kirkpatrick’s analogy will be discussed and a brief illustration of how this approach can be applied to atmospheric modeling will provide a hint of the method’s broad research utility.


Source Location with a Small-Aperture Array at Santiaguito Volcano, Guatemala

Aída Quezada-Reyes

During a volcanic eruption, infrasonic pressure waves generate low-frequency acoustic signals that reflect the impulsivity of the degassing source at the vent. In this work, I will determine the directions of acoustic sources using infrasound data recorded with a 3-microphone array deployed at Santiaguito’s lava dome complex, Guatemala, from January to April 2009. To achieve this, I will perform an exact inverse technique defined by the distance between the array sensors and the differences in time of arrival of the infrasound signals at the microphones. The later are related to the slowness vector and will be used to calculate the direction of plane waves given by the array back azimuth. Source locations are strongly influenced by the arrival time measurements of the infrasound signals and the number of microphones used. Error propagation will also be performed to determine the accuracy of my results. Source locations using this technique will provide valuable information to discriminate sources originated by eruptive events from those of non-volcanic origin.

Validation of a Distributed Parameter Watershed Model to Observations of Soil Moisture Using Inverse Methods

Todd Umstot

A distributed parameter watershed model has been developed that estimates the water balance at a location including precipitation, runoff, runon, soil water storage, net infiltration (e.g., recharge), evapotranspiration, sublimation, and snow melt. Although unsaturated water movement is a non-linear process, numerical methods cannot be practically applied at the watershed scale and so the vertical flow of water in the root zone is represented using a lumped empirical approach referred to as the “bucket model” method. Additionally, the watershed model has a large number of parameters that may or may not be significant to the estimation of root zone soil moisture. Using 27 years of soil moisture data collected in the Jornada Basin, the appropriateness of the bucket model approach and the sensitivity of model parameters will be tested using inverse methods.


Multi-station Deconvolution of Earthquake Source Time Functions

Katherine Anderson

The methodology of multi-station deconvolution of earthquake source time functions for body wave seismograms is presented in detail and applied to seismograms obtained from the IRIS DMC (via SOD) for five aftershocks of the 2010 Maule 8.8 Mw earthquake. Synthetic P and SH Green’s functions are computed for each event using a best double couple of the GCMT solution. Thereafter the synthetics were deconvolved by Green’s functions generated for a range of source depths using a model with a water layer over a uniform half space with a P-wave velocity of 6.0 km/s. The above method allows for the simultaneous inversion of several seismograms to obtain a ‘best’ source time function and source depth. Although it is the most stable linear inverse procedure in that it accounts for the inadequacies in the Green’s functions (spectral holes and a limited band of frequencies) by obtaining a more reliable source time function that fits the observed waveshapes remarkably well, ‘unexplained amplitude scatter’ occurs such that the synthetics must be scaled by arbitrary trace factors to account for the degradation of the solution.


Borehole Temperature Inversion

Holly Rotman

Borehole temperature reconstructions represent an inverse problem with nonunique solutions.   Theoretically, borehole temperatures can be used to resolve air temperatures over ~1000 years or more, depending on the depth of the boreholes sampled.  Researchers working on the problem have used multiple approaches: least squares inversion, singular value decomposition, and Markov Chain Monte Carlo.  Depending on the scale of the reconstruction, seasonal variations and regional lack of borehole data may further bias the solutions.  I will present methods and resultant solutions from work in Utah and in Alaska/Canada.  I will also present at least one regional to global climate reconstruction, with comments on its potential bias and other pitfalls.


Source Location and Sound Speed from Inversion of Infrasonic Data from Kilauea Volcano, Hawaii

Rebecca L. Johnson

Arrival time data from four near-field infrasound stations will be used to invert for the x and y coordinates of the infrasonic signal source and determine sound speed for events in Hawaii. The inversion will be treated as analogous to an iterative linear inversion for source location of an
earthquake whose signal travels through a homogeneous medium with an unknown velocity structure. A network of thirty-eight infrasound sensors was deployed near Halema'uma'u crater on Kilauea volcano in Hawaii between 2 and 4 August 2010. Continuous signal was recorded at 200 Hz for 72 hours. Infrasound instruments record acoustic signals traveling through the atmosphere with frequencies between 0.05 and 20 Hz. The initial inversion will use a subset of data from only four sensors. Reflections and refractions in a layered atmosphere will be neglected and linear ray paths assumed due to proximity between the presumed source and station locations. The z coordinates will be fixed by topography with the source elevation assumed to equal the elevation of Halema'uma'u crater. Sound speed is influenced by ambient temperature and wind velocities. The effects of wind speed will be neglected in the inversion but the first inversion iteration will use a seed sound speed value of 343 m/s, the speed of sound at 20 degrees Celsius. Error analysis will be performed using variance-covariance matrices. The inversion process will be repeated with data from additional stations to briefly investigate the effects on parameters for an overdetermined problem that is inconsistent due to noise. Tasks will be
completed in MATLAB. 


Seismic Scatterer Inversion

Robert Anthony

The Tramelli et al (2009) study of the spatial distribution and strength of seismic wave scatters within Mt. Vesuvius, Italy, is a cutting-edge example of the application of inverse techniques to solve complex geophysical problems. The study uses the coda of local earthquakes to obtain the energy residuals, which are a result of the spatial distribution and strength of scatters within the volcano. Thus, the coda waves produced by local earthquakes are manipulated to form the observables (d ) for the inversion of the spatial distribution and strength of the scatters (m). For this particular inverse problem, the G matrix is quite complex, requiring both the location of the earthquakes as well as a 3D velocity structure of the volcano. Ultimately, this inverse problem is carried through twice in order to obtain higher resolution locations. First the inversion is applied to a large volume with large cells to obtain a baseline estimate of the model parameters. This initial estimate is then used as starting model for a smaller volume divided into smaller cells. The results as well as some of the resolution and stability tests for this method will be presented.


Distance Effects on Frequency of Acoustic Signals Generated by Lightning

Elias Badillo

The goal of this work is to determine the effects of distance on the frequencies of acoustic signals produced by lightning. Multiple arrays of infra-sound sensors were deployed at different distances from a central location in the Magdalena mountains of New Mexico. The lightning mapping array (LMA) will be used as a reference for the location of RF, and acoustic sources. It is expected that there will be a 1/r relationship between the amplitude of the acoustic signals and the distance (r) from the acoustic sources to the sensors; and as it is already known, low frequencies are less attenuated than high frequencies, and will travel longer distances. Knowledge of the differences in frequency propagation could be used to perform tasks such as locating sources more accurately by selecting a frequency or band of frequencies that reach the sensors with higher signal to noise ratio. Assuming an isotropic medium and no multi-path effects due to topography, a point source will be used as a first order approximation, where the acoustic waves start as concentric circles of radius r (plane waves when they arrive to the acoustic arrays). The model can then be extended to line sources in order to take into account the main channels of a lightning flash assuming that each LMA source or groups of LMA sources emit an acoustic signal. Forward and inverse models will be considered to explore the propagation of acoustic signals due to lightning, and the effects of distance on the frequency of these signals. Since the function 1/r is nonlinear, linearization would be performed so that a linear model can be used for the relationship of frequency versus distance.


Large Parameter Estimation Problems and CAT Scans

Jason Mattax

An analysis of the computational demand in solving large parameter estimation problems is presented. Medical imaging, particularly Computed X-Ray Tomography (CT scan) data and models are considered. The number of parameters and constraints involved in this standard test case and solution difficulties are addressed. The use of standard ``direct" solving and iterative approaches is considered, in addition to model simplifications that may be possible in systems with many more constraints than parameters. Finally, a method that is capable of better parameter estimation is considered, as well as the computational costs it incurs.


Hydrus 1D Hydrological Modeling

Alicia Candelaria

Hydrus 1D is a critical zone or vadose zone modeling program which in the forward model takes meteorological data a soil properties and outputs the soil moisture profile and propagation curves. Hydrus 1D is also equipped with an inverse solution method in which the user inputs the soil moisture data metrological data and receives the soil hydraulic properties. This project uses an extensive data set, from the Jornada transect, to explore the inverse method of Hydrus 1D. This exploration will be done in two parts: First a solution will be found for weighted data and compared to the solution of non-weighted data, while the number of iterations and optimized parameters are held constant. In the second portion, the number of parameters optimized by Hydrus will be varied, while the data is not weighted and iterations held constant.


2-D Tomographic Inversion of Ross Island, Antarctica

Shoba Maraj

Seismic tomography is used to estimate certain properties of the Earth. More specifically, P-wave refraction tomography is used to to determine the Earth's velocity structure. Ross Island, Antarctica, is located in an active rift environment and is the site of the active Erebus volcano. During the 2008-2009 austral summer field season, 21 three-component 4.5Hz seismic recorders were deployed along a 76 km east-west profile along the island. The observed data from this experiment are the source-receiver traveltimes, while the model parameters are the slowness (reciprocal of velocities) values, as defined by an input forward model. Using seismic tomography on the dataset in Antarctica, we expect to determine the gross velocity structure below Ross Island and Erebus volcano. As this is a rift environment, it is of interest to determine the boundary between the crust and mantle (Moho), and possibly any velocity contrasts due to the volcanic conduit.


Estimating Pole/Zero Errors in GSN-IU Network Calibration Metadata

Adam Ringler

Mapping the voltage output of a seismometer of the digital record of a seismomograph into true ground motion requires correction of the data using a description of the instrument’s response.  For the Global Seismographic Network (GSN; Butler et al., 2004), as well as many other networks, this instrument response is represented as a Laplace pole/zero model and published in the Standard for the Exchange of Earthquake Data (SEED) format.  This Laplace representation assumes that the seismometer behaves as a linear system, with any abrupt changes described adequately via multiple time-invariant epochs.  The SEED format allows for published instrument response errors as well, but these typically have not been estimated or provided to users. We developed an iterative three-step method to estimate instrument response parameters (poles, zeros, and sensitivity and gain) and their associated errors using random calibration signals.  First, we solve a coarse non-linear inverse problem using a least squares grid search to yield a first approximation to the solution.  This approach reduces the likelihood of poorly estimated parameters (a local-minimum solution) caused by noise in the calibration records and enhances algorithm convergence.  Second, we iterative solve a non-linear parameter estimation problem to obtain the least squares best-fit Laplace pole/zero/gain model.  Third, by applying the central limit theorem we estimate the errors in this pole/zero model by solving the inverse problem at each frequency in a 2/3rds-octave band centered at each best-fit pole/zero frequency.  This procedure yields error estimates of the 99% confidence interval. We demonstrate this method by applying it to a number of recent IRIS/USGS (IU) GSN calibrations.


Path Delay Correction for the Magdalena Ridge Observatory Interferometer

Tyler McCracken

Stellar interferometers are capable of spatial resolutions on the sky unmatched by any other conventional telescopes. To gain this resolution, light from multiple telescopes is combined creating a fringe pattern that simulates one large aperture. The path length travelled by the light from each telescope to combination plane must be within a coherence length, defined by the mean observing wavelength squared over the bandwidth being observed, in order for a fringe pattern to be present. At the Magdalena Ridge Observatory, a method termed group delay fringe tracking is used to track the position of fringes and keep the interferometer within the imposed tolerances. Group delay fringe tracking is commonly built around an FFT based algorithm that is used to nd the peak in a power spectrum of several integrated data frames. Of particular interest is the effect of the atmosphere. Its constant fluctuations will inhibit the effectiveness of the algorithm and limits the number of data frames that can be integrated together. This paper provides an overview of group delay fringe tracking and the role the atmosphere plays in limiting its performance.


Earthquake b-Values and The Best Constraints to Calculate Them

Emily Morton

Earthquakes are hazardous processes and scientists as well as the general population have wanted to be able to predict them for a long time.  A relative way to “predict” earthquakes is to look at the distribution of earthquake magnitude to the log of earthquake frequency per year, called the Gutenberg-Richter Distribution.  This produces a linear trend that can potentially tell us when to expect large, dangerous earthquakes in relation to small, harmless earthquakes using the b-value, or slope of the line.  However, in this distribution the smaller earthquakes tend to be under-sampled due to instrument sensitivity and the larger earthquakes tend to be incomplete due to a lack of recorded earthquakes.  Therefore, to calculate a b-value, we must limit the magnitudes used to solve for the line parameters.  Looking at a global catalog of earthquakes from 1973 to the present, this paper will explore how to calculate the global b-value, what it means, and find the ideal constraints with which to calculate the b-value.  I am predicting that a plot relating the lower limit of magnitude and the upper limit of magnitude will show an area of stable b-value and areas of unstable b-values, with the stable area corresponding with the best b-value constraints.


Separation of Karst Spring Discharge into Baseflow and Quickflow Sources

Katrina Koski

The chemistry and discharge measurements from springs at the distal end of underwater karst caves (conduits) can provide information about sources of recharge, flow paths, and aquifer water rock interactions. Knowledge of these processes can be used to describe dissolution of minerals, and contaminant and nutrient transfer in the aquifer. The goal of this research is to separate spring discharge during storm events into two components: baseflow due to groundwater drainage from the rock, and quickflow from recent surface water inputs due to a storm event. Using a well studied system where percent baseflow and quickflow have been determined for a specific storm event by other means, the spring chemistry, physical water properties, and discharge rate (e.g. conductivity, temperature, and flow rate respectively) will be used individually and in combination in search of an inverse model which can be used to specify baseflow and quickflow for the storm event. An attempt will be made to generalize the model for other storm events in the system. Discussion includes the possibility of generalization to other systems with similar geometry, recharge conditions, and surrounding rock.


Blurring and Deblurring

Gunter Leguy

When taking a picture with a digital camera, one of the primary causes of image blurring comes from camera shake. Usually this happens during low light condition and slow shutter speed. When the blur kernel of the picture is known, there exists different models to solve this problem such as Weiner filtering or Richardson-Lucy deconvolution. However, it is not necessary trivial to know information about the noise of the blurred image. In this paper we will first show an example of filtering using deconvolution principle assuming the blur kernel known. We will apply this technique to a picture to which we add the noise ourselves. Such a practice is called a non blind deconvolution. Secondly, we will focus on a mathematical algorithm which deals with blind deconvolution. It is known
that convolution methods work pretty well if the blurry image does not contain noise and the blur kernel does not contain error. Based on these facts, a probabilistic approach combining both, blur kernel estimation and un-blurred image restoration can be developed and lead to a Maximum a posteriori problem to solve.


2010 Homework Assignments


Note that the MATLAB data files for relevant problems can be found on the cd volume (

Homework 1: Due 9/13 (RA):

Exercises: 1.1 (15 pts), 1.2 (10 pts), A3 (10 pts), A5 (10 pts), A6 (10 pts), A7 (10 pts), A18 (10 pts), A19 (10 pts) from the textbook.

New Problem (15 pts):

a) Write a MATLAB subroutine, quadform(K), that diagonalizes a general 2 by 2 positive definite matrix, K, and plots ts associated ellipsoidal bounding region (Equation A.84) using Delta = 1.
You may find the plotting code noted in Exercise 2.1(c), and using the "axis square" and "grid" commands to be helpful here.

b) Produce the appropriate plot for the matrix:

K= [ 2.5000 -0.5000 ; -0.5000 2.5000 ];

c) What are the exact semiaxis directions and lengths for your plotted ellipsoid in this case, and describe in words their relationship to the elements of the eigenvalue-eigenvector decomposition.


Homework 2: Due 9/27 (RA):

Exercises: B.1 (10 pts), B.8 (20 pts), 2.2 (20 pts), 2.3 (30 pts), 2.5 (20 pts) from the textbook.


Homework 3: Due 10/11 (RA):

Exercises: 3.1 (20 pts), 4.1 (20 pts), 4.3 (40 pts), 4.4 (20 pts) from the textbook (you may find a Piccard plot to be useful in the analysis of problem 4.4).


Homework 4: Due 11/1 (OM):

Exercises: 5.2 (35 pts), 5.3 (35 pts), 6.1 (20pts), 6.3 (10pts) from the textbook.


Homework 5: Due 11/22 (OM):

Each problem is worth 25 points.

1) A 30 by 30 by 30 m cube is parameterized into 1 cubic meter cubes and imaged with seismic tomography using 10,000 raypaths traversing each of the three parallel side pairs.

Solve and plot (show horizontal 30 by 30 slices for z=1 to 30 in compact figures using subplot and with a color axis range of +- 0.03 slowness units) the slowness variations of the block using Kazmarz' algorithm, then again using the ART method. To map the slowness vector into a 3-dimensional 30 by 30 by 30 matrix, use the reshape command.

Extra credit: Use the cputime command to find the necessary cpu time to perform the inversion on two computers.

The data file (ttn; travel time residual relative to a uniform slowness model in s) and the 30,000 by 27,000 sparse G matrix for the problem can be found as a .mat file here

2) Using the data and forward model from Exercise 3.1, estimate m(x) using least squares and the following Fourier basis on the unit interval:

h_0(x) = 1
h_1(x) = cos(pi*x)
h_2(x) = sin(pi*x)
h_3(x) = cos(2*pi*x)
h_4(x) = sin(2*pi*x)

No regularization is necessary (why?) You may numerical compute the integrals necessary to set up the necessary matrix G for calculating the basis expansion coefficients.

3) Problem 9.1 from the textbook.

4) Problem 9.5 from the textbook.


Homework 6: Due 12/6 (OM): Emailed to on 11/25.




EES Home | Geology | Geochemistry | Hydrology | Geophysics | Site Map | Search

Last Updated: January 16, 2012
Please contact Rick Aster regarding content on this page.