This introduction is a brief overview of the moment tensor inversion problem as handled by the MTUQ package. It is directly copied from the Google Colab notebook that I used for the SCOPED 2024 Workshop session on moment tensor inversion. The notebook can be found here.

Please reach out to me if you have any issues or questions about the content of the Colab notebook.


Introduction to Moment tensor inversion with the Cut-and-Paste method

In this brief introduction, we will see what the moment tensor inversion problem is about. This should cover all the theory you should be familiar with before handling MTUQ, starting with:

The seismic moment tensor

The seismic moment tensor $\mathbf{M}$ is a mathematical representation of the motion applied to a point-sized seismic source. It is a $3 \times 3$ second-order symmetric tensor with 6 independent parameters.

Representation of the double couple on a fault
Representation of the moment tensor $\mathbf{M}$ components in terms of linear vector dipoles and force couples. From Stein & Wysession (2002).

In the ideal case of a motion applied onto a perfect fault-plane, it can be represented by a supperposition of force-couples or dipoles (with zero net-torque):

Representation of the double couple on a fault
Motion on an ideal fault plane represented as double-couples. From Stein & Wysession (2002).

Therefore, by estimating the moment tensor parameters for seismic sources that fall into the point-source approximation we can learn about their geometrical properties, and the general motion that was applied to the subsurface - generating seismic waves.


A moment tensor can be conveniently represented as the so-called beachball, which is a representation of the P-wave first motions in 3D space (binary representation of the P-waves radiation pattern). The color of the beachball reflects the direction of motion applied in different directions. It is generally white for motion pointing inward and black (or color shaded) for motion pointing outward of the unit sphere.

Gif of a moment tensor in 3D with vector field for the motion
Vector field representation of the motion associated with double-couple point source. represented as a 'beachball'. Courtesy of Stuart Tierney.

Most of the time, moment tensor “beachballs” are visualized as 2D diagrams. These diagrams are stereonet projections of the lower hemisphere of the unit sphere, as seen from above in a map view. This representation allows for easier interpretation and analysis, as it conveys the general motion on a fault at a glance, or the source type in the case of non double-couple sources.

A few key moment tensor and fault block-models.
Common double couple moment tensor and associated fault block-models USGS.

Moment tensor parameterization and the Lune diagram

Moment tensors attributes can be broadly defined as source type and orientation, which are both informed by the geometrical properties of the tensor. Given a moment tensor $\mathbf{M}$, whose eigvenvalues $(\lambda_1, \lambda_2, \lambda_3)$ satisfy $\lambda_1 \geq \lambda_2 \geq \lambda_3$, we can represent a moment tensor as its eigendecomposition:

$\mathbf{M} = \mathbf{U \Lambda U}^T$

where $\mathbf{U}$ contains the eigenvectors that satisfies $\mathbf{UU}^T = \mathbf{I}$, and $\mathbf{\Lambda}$ is the diagonal matrix containing the eigenvalue triple in descending order.

With this decomposition in mind, we see that $\mathbf{U}$ encodes the source orientation (the tensor’s eigenframe), and $\mathbf{\Lambda}$ encodes the source type. From this decomposition stems a very natural way of representing source-types, by projecting the tensor position on the eigenvalue lune:

Geometry of the eigenvalue lune
Definition of the eigenvalue lune as a set of ordered and normalized eigenvalues. From Tape & Tape 2019.

In the following example, we show for 4 fixed eigenframes, how the moment-tensor source type changes as a function of $\mathbf{\Lambda}$ only:

4 moment tensor lunes
Moment-tensor source type plot for a set of 4 fixed source orientations. For a fixed orientation, the location on the eigenvalue lune controls the source type. From Carl Tape figure gallery.

In MTUQ, we adopt a moment-tensor parameterization that exploits the geometrical properties of the moment tensor so that we explore the moment-tensor space in a uniform manner. Instead of exploring the 6 parameters of the moment tensor directly, we will express moment tensors with the following set of coordinates:

such that we have $\mathbf{M}(ρ, γ, δ, κ, σ, θ)$.

For more details, you can read Tape & Tape (2015).

Moment tensor inversion

In MTUQ we solve the moment tensor inversion problem by finding the source that minimizes the following cost function:

$C(\mathbf{M}) = \sum_{r=1}^{N} \sum_{i=1}^{3} \int_0^T \left[d_i^r(t,\mathbf{M}) - d_i^{r,obs}(t)\right]^2 dt$

which translates to the sum of least-squares waveform difference between the synthetic waveforms $d(t)$ and the observed data $d^{obs}(t)$ for N stations with 3 components.

In order to speed up the optimization process, we rely on a set of pre-computed Green’s functions $\mathbf{G}^r$ between each source-receiver pair $r$, such that we can recast the problem as

$C(\mathbf{M}) = \sum_{r=1}^{N} \sum_{i=1}^{3} \int_0^T \left[\mathbf{G}_i^r \mathbf{m} - d_i^{r,obs}\right]^2 dt$

where the time dependency has been dropped for convenience, and $\mathbf{m}$ is a vector containing the 6 independent parameters of the Moment tensor $\mathbf{M}$. Using pre-computed Green’s functions allows the rapid evaluation of the residual term, which we simplify as:

$\mathbf{r} = \mathbf{Gm} - \mathbf{d}$

which is less expensive than explicitly modeling the waveform for each source in a given earth model.

The Cut-and-paste method

The algorithm underlying MTUQ is known as the “cut-and-paste” (CAP)method of Zhu & Helmberger (1996).

The main goal of the CAP method is to mitigate the effects of 3-D velocity structures that are unaccounted for with a 1D reference model (or an inexact 3D model).

It does so by splitting all the data into specific window groups (typically P-waves and Surface waves) and allowing a time-shift window for the synthetics. Within this allowable time-shift window, a cross-correlation is performed between $d_{syn}$ and $d_{obs}$, and the residual is computed only for the maximum cross-correlation value (when signals are supposedly aligned).

Illustration of the cut and paste method.
Illustration of the cut and paste method.

We do this for each data group (P, Rayleigh, Love waves), for each station and each component (Z-R-T), and return one misfit value corresponding to the maximum cross-correlation value.