3.2 Waveform-based Moment Tensor Inversion

This is a follow-up of the Amplitude-based Moment Tensor Inversion tutorial. In this tutorial, we will extend the MT inversion algorithm to work directly with waveforms instead of picked amplitudes. By avoiding any picking, this method can determine the moment tensor of a microseismic source when the source location is not known a-prior. As such, it can be considered a joint location and MT inversion algorithm.

We will start once again from the far-field particle velocity expression from a moment tensor source in a homogeneous full space (from eq. 4.29, Aki and Richards) - see the Amplitude-based Moment Tensor Inversion tutorial for more details.

However, in comparison to the Amplitude-Based Moment Tensor Inversion tutorial, in this waveform-based approach we assume a distributed source within a subsurface area of interest and use the following integral relation to reconstruct the surface data:

\[v_1^P(\mathbf{x_r}, t) = w(t) * \int\limits_V G_{pq}(\mathbf{x_r}, \mathbf{x}, t) M_{pq}(\mathbf{x})\,\mathrm{d}\mathbf{x}\]

where \(M_{pq}\) with \(p,q=1,2,3\) are the so-called MT kernels and \(G_{pq}\) are the so-called Green’s functions, whose high-frequency approximation can be written as:

\[G_{pq}(\mathbf{x_r}, \mathbf{x}, \omega) = a_{pq}(\mathbf{x_r}, \mathbf{x}) e^{j \omega t(\mathbf{x_r}, \mathbf{x})}\]

Here \(a_{pq}\) with \(p,q=1,2,3\) represent the same coefficients used in the Amplitude-based Moment Tensor Inversion tutorial.

To summarize, we will apply the following workflow:

• Load model and data;

• Compute the traveltimes & ray angles

• Compute the Greens functions for the subsurface area of interest

• Define the Kirchhoff-MT operator

• Jointly solve for the location and MT with a least-squares solver

Assumptions: for now, the MTWI procedure assumes a homogeneous velocity model.

import os
import numpy as np
import matplotlib.pyplot as plt

import fracspy

from pylops.utils.wavelets import ricker
from fracspy.utils.sofiutils import read_seis
from fracspy.mtinversion.utils import get_mt_max_locs, get_mt_at_loc
from fracspy.mtinversion.mtwi import *
from fracspy.visualisation.momenttensor_plots import MTMatrix_comparisonplot
# sphinx_gallery_thumbnail_number = 7

Load model and seismic data

For this example, we will use a toy homogenous model with a gridded surface receiver array. The data are modelled using the SOFI3D Finite Difference modelling software. The model is the same that we have used in the FD modelling to generate the data. As such, it contains additional boundaries, which we need to remove prior to performing localisation.

default_cmap = 'seismic'

# Directory containing input data
input_dir = '../data/pyfrac_SOFIModelling'

# Model parameters
abs_bounds = 30
dx = dy = dz = 5
nx = 112
ny = 128
nz = 120

# Modelling parameters
dt = 1e-3  # SOFI3D Time sampling rate
t_shift = 167  # Time shift required to align FD data to zero time for Kirchhoff operators
tdur = 500  # Recording duration

# Load model
mod_w_bounds = np.fromfile(os.path.join(input_dir,'inputs',
                                        'models',
                                        'Homogeneous_xyz.vp'),
                           dtype='float32').reshape([nx, ny, nz])

# Load receiver geometry
recs_xzy = np.loadtxt(os.path.join(input_dir,'inputs',
                                   'griddedarray_xzy_20m.dat')).T
nr = recs_xzy.shape[1]

# Load seismic data (note that Vz is Vy given the SOFI convention)
expname = 'MT-90-90-180_Homogeneous_griddedarray'
vz = read_seis(os.path.join(input_dir, 'outputs', 'su',
                            '%s_vy.txt' % expname),
               nr=nr)
vz = vz[:, t_shift: t_shift + tdur]
efd_scaler = np.max(abs(vz))  # Scaler to make data more friendly
vz /= efd_scaler

# Remove absorbing boundaries for both the model and receiver coordinates
mod = mod_w_bounds[abs_bounds:-abs_bounds, abs_bounds:-abs_bounds, :-abs_bounds] # z has free surface
nx, ny, nz = mod.shape
x, y, z = np.arange(nx)*dx, np.arange(ny)*dy, np.arange(nz)*dz
recs = np.array([recs_xzy[0]-(abs_bounds*dx), recs_xzy[2]-(abs_bounds*dx), recs_xzy[1]])

Let’s now double-check that the data has been loaded correctly. Observe the changes in polarity across the traces, this is the information that we utilise to determine the Moment Tensor.

fig, ax = fracspy.visualisation.traceviz.traceimage(vz, climQ=99.99, figsize=(10, 4))

Create modelling operator

First, we will define a Ricker wavelet with peak frequency of 20Hz. This is the same wavelet that we used in modelling; in real applications, this will need to be estimated from the data.

omega_p = 20
nt = vz.shape[1]
t = np.arange(nt)*dt
wav, wavt, wavc = ricker(t[:81], f0=omega_p)

fig, ax = plt.subplots(figsize=(10, 2))
ax.plot(wavt, wav, 'k', lw=2)
ax.set_xlabel('t [s]')
ax.set_title('Wavelet')
ax.grid('on')
plt.tight_layout()

Second we define an area of interest where we expect the source to be located. In fact, whilst in practice one could consider the entire subsurface, this comes with a computational and storage burden for the Green’s functions.

sx = nx // 2
sy = ny // 2
sz = 2 * nz // 3
sloc_ind = [sx, sy, sz]

hwin_nx_aoi, hwin_ny_aoi, hwin_nz_aoi = 15, 13, 11  # half window lengths in x, y, z
winc_x, winc_y, winc_z = nx // 2, ny // 2, 2 * nz // 3  # Center points of the area of interest

# Defining area of interest
xsi, xfi = winc_x-hwin_nx_aoi, winc_x+hwin_nx_aoi+1   # start/end index of x-region of interest
ysi, yfi = winc_y-hwin_ny_aoi, winc_y+hwin_ny_aoi+1   # start/end index of y-region of interest
zsi, zfi = winc_z-hwin_nz_aoi, winc_z+hwin_nz_aoi+1   # start/end index of z-region of interest
nx_aoi = xfi - xsi
ny_aoi = yfi - ysi
nz_aoi = zfi - zsi

# MT in area of interest
MT_aoi = np.zeros([6, nx_aoi, ny_aoi, nz_aoi])  # MT components as images
MT_selected = -1 * np.array([0, 0, 0, 1, 0, 0])
MT_aoi[:, nx_aoi//2, ny_aoi//2, nz_aoi//2] = MT_selected

Next, we create our Kirchhoff-MT operator

Ms_scaling = 1.92e10
mtw = MTW(x, y, z, recs, mod, sloc_ind,
          2, omega_p, (hwin_nx_aoi, hwin_ny_aoi, hwin_nz_aoi),
          t, wav, wavc, multicomp=False,
          Ms_scaling=Ms_scaling,
          engine='numba')
data = mtw.model(MT_aoi)

Joint localisation and MT inversion

Finally, we are ready to invert our waveform data for the 6 MT kernels.

# Adjoint
mt_adj = mtw.adjoint(vz)

# Inversion
mt_inv = mtw.lsi(vz, niter=100, verbose=True)

LSQR            Least-squares solution of  Ax = b
The matrix A has 71500 rows and 115506 columns
damp = 0.00000000000000e+00   calc_var =        0
atol = 0.00e+00                 conlim = 1.00e+08
btol = 0.00e+00               iter_lim =      100

   Itn      x[0]       r1norm     r2norm   Compatible    LS      Norm A   Cond A
     0  0.00000e+00   2.334e+01  2.334e+01    1.0e+00  4.2e+02
     1 -1.15762e-06   1.355e+01  1.355e+01    5.8e-01  4.5e-01   1.2e+04  1.0e+00
     2 -1.71343e-06   8.313e+00  8.313e+00    3.6e-01  3.0e-01   1.4e+04  2.5e+00
     3 -7.81990e-07   6.895e+00  6.895e+00    3.0e-01  2.1e-01   1.7e+04  3.8e+00
     4  4.13818e-07   6.036e+00  6.036e+00    2.6e-01  2.5e-01   1.9e+04  5.4e+00
     5  1.64581e-06   4.969e+00  4.969e+00    2.1e-01  1.1e-01   2.2e+04  8.0e+00
     6  1.42365e-06   4.358e+00  4.358e+00    1.9e-01  1.2e-01   2.3e+04  1.0e+01
     7  1.02796e-06   4.057e+00  4.057e+00    1.7e-01  1.2e-01   2.5e+04  1.2e+01
     8 -1.57272e-08   3.647e+00  3.647e+00    1.6e-01  1.2e-01   2.7e+04  1.5e+01
     9 -8.50817e-07   3.404e+00  3.404e+00    1.5e-01  6.3e-02   2.9e+04  1.7e+01
    10 -1.20292e-06   3.311e+00  3.311e+00    1.4e-01  6.5e-02   3.1e+04  1.9e+01
    90 -6.78468e-06   8.598e-01  8.598e-01    3.7e-02  4.1e-03   8.6e+04  5.6e+02
    91 -6.48549e-06   8.572e-01  8.572e-01    3.7e-02  5.9e-03   8.6e+04  5.7e+02
    92 -6.31525e-06   8.558e-01  8.558e-01    3.7e-02  5.4e-03   8.7e+04  5.7e+02
    93 -6.01709e-06   8.533e-01  8.533e-01    3.7e-02  4.0e-03   8.8e+04  5.8e+02
    94 -5.46185e-06   8.488e-01  8.488e-01    3.6e-02  8.5e-03   8.8e+04  6.0e+02
    95 -5.13239e-06   8.462e-01  8.462e-01    3.6e-02  4.8e-03   8.9e+04  6.1e+02
    96 -4.70270e-06   8.428e-01  8.428e-01    3.6e-02  5.3e-03   8.9e+04  6.2e+02
    97 -4.54091e-06   8.417e-01  8.417e-01    3.6e-02  5.2e-03   8.9e+04  6.3e+02
    98 -4.18106e-06   8.392e-01  8.392e-01    3.6e-02  5.2e-03   9.0e+04  6.4e+02
    99 -3.97774e-06   8.380e-01  8.380e-01    3.6e-02  5.8e-03   9.1e+04  6.4e+02
   100 -3.61493e-06   8.358e-01  8.358e-01    3.6e-02  5.0e-03   9.1e+04  6.6e+02

LSQR finished
The iteration limit has been reached

istop =       7   r1norm = 8.4e-01   anorm = 9.1e+04   arnorm = 3.8e+02
itn   =     100   r2norm = 8.4e-01   acond = 6.6e+02   xnorm  = 8.6e-03

Let’s now compare the expected and true MT source parameters at the true location. Note that for the expected MT parameters, we display their normalized version since the modelling operator is only accurate up to relative amplitudes.

mt_at_loc = get_mt_at_loc(mt_inv, [sloc_ind[0]-xsi, sloc_ind[1]-ysi, sloc_ind[2]-zsi])

MTMatrix_comparisonplot(MT_selected, mt_at_loc / np.abs(np.array(mt_at_loc)).max())

Finally we visualize the estimated kernels from the inversion.

clim = 1e-4

fracspy.visualisation.eventimages.locimage3d(mt_inv[0], sloc_ind[0]-xsi, sloc_ind[1]-ysi, sloc_ind[2]-zsi,
                                             clipval=[-clim, clim], title='Mxx')
fracspy.visualisation.eventimages.locimage3d(mt_inv[1], sloc_ind[0]-xsi, sloc_ind[1]-ysi, sloc_ind[2]-zsi,
                                             clipval=[-clim, clim], title='Myy')
fracspy.visualisation.eventimages.locimage3d(mt_inv[2], sloc_ind[0]-xsi, sloc_ind[1]-ysi, sloc_ind[2]-zsi,
                                             clipval=[-clim, clim], title='Mzz')
fracspy.visualisation.eventimages.locimage3d(mt_inv[3], sloc_ind[0]-xsi, sloc_ind[1]-ysi, sloc_ind[2]-zsi,
                                             clipval=[-clim, clim], title='Mxy')
fracspy.visualisation.eventimages.locimage3d(mt_inv[4], sloc_ind[0]-xsi, sloc_ind[1]-ysi, sloc_ind[2]-zsi,
                                             clipval=[-clim, clim], title='Mxz')
fracspy.visualisation.eventimages.locimage3d(mt_inv[5], sloc_ind[0]-xsi, sloc_ind[1]-ysi, sloc_ind[2]-zsi,
                                             clipval=[-clim, clim], title='Mzz')

• 
• 
• 
• 
• 
• 

(<Figure size 800x800 with 5 Axes>, (<Axes: xlabel='x samples', ylabel='z samples'>, <Axes: xlabel='y samples'>, <Axes: ylabel='y samples'>))