Surface-Wave Tomography#

Open In Colab

In this notebook, we will apply CoFI to measurements of surface-wave velocity collected across the USArray from the ambient seismic noise. Specifically, we will retrieve, through CoFI and SeisLib, a Rayleigh-wave phase velocity map of the Conterminous United States at the surface-wave period of 5 s. The employed velocity measurements belong to the data set compiled by Magrini et al. (2022).

If you are running this notebook locally, make sure you’ve followed these steps to set up the environment. (This environment.yml file specifies a list of packages required to run the notebooks)

Theoretical Background#

To map lateral variations in surface-wave velocity, SeisLib implements a least-squares inversion scheme based on ray theory. This method rests on the assumption that surface waves propagate, from a given point on the Earth’s surface to another, without deviating from the great-circle path connecting them. Under this assumption, the traveltime along the great-circle path can be written \(t = \int_{\mathrm{path}}{s(\phi(l), \theta(l)) dl}\), where \(\phi\) and \(\theta\) denote longitude and latitude, and \(s\) the sought Earth’s slowness.

Let us consider a discrete parameterization of the Earth’s surface, and assume each block (or grid cell) of such parameterization has constant slowness. The above integral expression can then be reformulated in the discrete form

\[\label{eq:forward_problem}\tag{1} s = \frac{1}{L} \sum_{n}{s_n l_n},\]

where \(L\) is the length of the great-circle path and \(l\) the distance traveled by the surface wave through the \(n\)th block. Equation (\(\ref{eq:forward_problem}\)) represents the forward calculation that allows for retrieving the average velocity of propagation between two points on the Earth’s surface (i.e., the quantity which is typically measured in ambient-noise seismology), provided that the (discrete) spatial variations in velocity (or slowness) are known.

If we now define the \(m \times n\) matrix such that \(A_{ij} = \frac{l_j}{L_i}\), where \(L_i\) is the length of the great circle associated with \(i\)th observation, we can switch to matrix notation and write

\[\label{eq:forward_matrix}\tag{2} {\bf A \cdot x} = {\bf d},\]

where \(\bf d\) is an \(m\)-vector whose \(k\)th element corresponds to the measured slowness, and \(\bf x\) the sought \(n\)-vector whose \(k\)th element corresponds to the model coefficient \(s_k\). Matrix \(\bf A\), also known as “data kernel” or “Jacobian”, is computed numerically in a relatively simple fashion. For each pair of receivers for which a velocity measurement is available, its \(i\)th entries is found by calculating the fraction of great-circle path connecting them through each of the \(n\) blocks associated with the parameterization.

In geophysical applications, the system of linear equations (\(\ref{eq:forward_matrix}\)) is usually ill-conditioned, meaning that it is not possible to find an exact solution for \(\bf x\). (In our case, it is strongly overdetermined, i.e. \(m \gg n\).) We overcome this issue by first assuming that the target slowness model is approximately known, i.e. \({\bf x}_0 \sim \bf{x}\). We then invert for the regularized least-squares solution

\[\label{eq_inverse_prob}\tag{3} {\bf x} = {\bf x}_0 + \left( {\bf A}^T \cdot {\bf A} + \mu^2 {\bf R}^T \cdot {\bf R} \right)^{-1} \cdot {\bf A}^T \cdot ({\bf d} - {\bf A} \cdot {\bf x}_0),\]

where the roughness of the final model is determined by the scalar weight \(\mu\) and the roughness operator \(\bf R\) is dependent on the parameterization (for technical details on its computation, see Magrini et al. (2022)).

1. Data and Parameterization#

As mentioned earlier, the data used in this notebook consist of inter-station measurements of Rayleigh-wave phase velocity at 5 s period. We parameterize the Earth’s surface through equal-area blocks of \(1^{\circ} \times 1^{\circ}\).

from seislib.tomography import SeismicTomography

tomo = SeismicTomography(cell_size=1) # Parameterization

# To reproduce the results locally, download data.txt and change the below path

tomo.add_data(src='./data.txt')
-------------------------------------
Optimal grid found in 91 iterations
-------------------------------------
-------------------------------------
GRID PARAMETERS
Lonmin - Lonmax : -180.000 - 180.000
Latmin - Latmax : -90.000 - 90.000
Number of cells : 41252
Grid cells of 1.000° : 41252
-------------------------------------
DATA PARAMETERS
Lonmin - Lonmax data : -124.566 - -67.312
Latmin - Latmax data : 24.727 - 49.098
Number of measurements : 171353
Source : ./data.txt
-------------------------------------

Overall, 171,353 velocity measurements are available (check tomo.velocity), each associated with a different pair of receveirs (check tomo.data_coords, consisting of a matrix of 171,353 rows and 4 columns: \(\theta_1\), \(\phi_1\), \(\theta_2\), and \(\phi_2\)).

2. Jacobian#

We use the information about the data coordinates to calculate the matrix \(\bf A\) (i.e. the Jacobian). In doing so, we will discard all blocks parameterizing the Earth’s surface that are not intersected by at least one inter-station great-circle path. These model parameters (referred to as “grid cells” in the below output) have no sensitivity to our data.

# This discards all blocks that are far away from the study area

tomo.grid.set_boundaries(latmin=tomo.latmin_data,
                         latmax=tomo.latmax_data,
                         lonmin=tomo.lonmin_data,
                         lonmax=tomo.lonmax_data)
*** GRID UPDATED ***
-------------------------------------
GRID PARAMETERS
Lonmin - Lonmax : -125.774 - -65.926
Latmin - Latmax : 23.999 - 50.002
Number of cells : 1201
Grid cells of 1.000° : 1201
-------------------------------------
# Computes the coefficients of the A matrix, while discarding all model parameters that are not constrained by our data.
tomo.compile_coefficients(keep_empty_cells=False)
*** GRID UPDATED ***
-------------------------------------
GRID PARAMETERS
Lonmin - Lonmax : -125.774 - -66.537
Latmin - Latmax : 23.999 - 50.002
Number of cells : 775
Grid cells of 1.000° : 775
-------------------------------------

The Jacobian can now be accessed by typing tomo.A, and the associated parameterization can be visualized by typing

tomo.grid.plot()
sw tomography
<GeoAxes: >

3. Inversion – SeisLib style#

The lateral variations in phase velocity can now simply be retrieved, via SeisLib, through

mu = 5e-2 # Roughness damping coefficient, arbitrarily chosen

# The output of tomo.solve is slowness, hence we take the reciprocal

c = 1 / tomo.solve(rdamp=mu) # in km/s

Let’s have a look at the results.

from seislib.plotting import plot_map
import seislib.colormaps as scm

img, cb = plot_map(tomo.grid.mesh, c, cmap=scm.roma, show=False)
cb.set_label('Phase velocity [km/s]')
sw tomography

4. Inversion – CoFI style#

Let’s now reproduce the above results through CoFI. First, we need to define a starting model \({\bf x}_0\) to compute the residuals \({\bf r} = {\bf d} - {\bf A} \cdot {\bf x}_0\), as in equation (3).

import numpy as np

A = tomo.A # Jacobian
x0 = np.full(A.shape[1], 1 / tomo.refvel) # tomo.refvel is the average inter-station phase velocity
d = 1 / tomo.velocity # measurements of (average) inter-station slowness
r = d - A @ x0 # residuals

We now need to define the roughness operator \(\bf R\). This is done under the hood by SeisLib through the “private” method _derivatives_lat_lon.

from seislib.tomography._ray_theory._tomography import _derivatives_lat_lon

# coordinates of each parameterization block, tomo.grid.mesh, should be in radians

R_lat, R_lon = _derivatives_lat_lon(np.radians(tomo.grid.mesh))
R = np.row_stack((R_lat, R_lon))

Almost everything is ready to carry out the inversion through CoFI. Before doing so, we need to define our inverse problem (through BaseProblem) and pass to it the data and the Jacobian (through set_data and set_jacobian). Finally, we will specify the regularizazion criterion (through set_regularization).

from cofi import BaseProblem
from cofi.utils import QuadraticReg

problem = BaseProblem()
problem.set_data(r) # our data are now the residuals defined above
problem.set_jacobian(A)

# As opposed to SeisLib, CoFI does not square the damping coefficient.
problem.set_regularization(mu**2 * QuadraticReg(R, (A.shape[1],)))   # L2 norm of R, i.e. R.T @ R

problem.summary()
=====================================================================
Summary for inversion problem: BaseProblem
=====================================================================
Model shape: Unknown
---------------------------------------------------------------------
List of functions/properties set by you:
['jacobian', 'regularization', 'data']
---------------------------------------------------------------------
List of functions/properties created based on what you have provided:
['jacobian_times_vector']
---------------------------------------------------------------------
List of functions/properties that can be further set for the problem:
( not all of these may be relevant to your inversion workflow )
['objective', 'log_posterior', 'log_posterior_with_blobs', 'log_likelihood', 'log_prior', 'gradient', 'hessian', 'hessian_times_vector', 'residual', 'jacobian_times_vector', 'data_misfit', 'regularization_matrix', 'forward', 'data_covariance', 'data_covariance_inv', 'initial_model', 'model_shape', 'blobs_dtype', 'bounds', 'constraints']

We now carry out the inversion through scipy.linalg.lstsq.

from cofi import Inversion, InversionOptions

options = InversionOptions()
options.set_tool("scipy.linalg.lstsq")

inv = Inversion(problem, options)
inv_results = inv.run()
inv.summary()
=======================================
Summary for Inversion
=======================================
Completed with the following result:

Summary for inversion result
SUCCESS
model: [ 1.18179454e-06 -1.25177071e-05 -1.09358993e-05 -4.91566292e-06
  3.30435896e-06  5.65442303e-06  4.76237564e-06  8.81902823e-06
  1.08410009e-05  1.51420559e-05  2.25037965e-05  2.36951885e-05
  2.04676927e-05  1.50074366e-05  7.02904303e-06  5.04035200e-08
 -8.04395501e-08  3.98518542e-05  1.39452438e-05 -1.54276927e-05
 -1.49697038e-05 -8.59247577e-06 -1.36969843e-05  5.07935541e-06
  7.80266500e-06 -1.16544344e-05  4.66301865e-06  1.03449957e-05
  1.01949660e-05  1.90675486e-05  4.40184003e-05  4.34500804e-05
  8.43876299e-06 -3.07874077e-06 -2.18481778e-05 -2.57935726e-05
 -2.04246725e-05 -2.34561998e-05 -2.85126659e-05 -6.86848259e-06
 -3.66257498e-08  4.26711646e-05  6.04471855e-05  5.06158158e-06
 -4.78646988e-06 -1.25434180e-05 -1.24568555e-05 -9.42127235e-06
 -1.82659875e-06  5.55900889e-06 -1.46681175e-05 -4.53840079e-06
  8.30662686e-06  2.17424298e-05  3.14617201e-05  4.49551156e-05
  5.28191023e-05  2.11956148e-05 -2.62068784e-06 -2.31795897e-05
 -2.77380865e-05 -2.49902336e-05 -3.00944276e-05 -1.95982340e-05
  1.74828943e-05  1.17275126e-05  1.81609410e-05  1.32380112e-05
  1.11335859e-06 -1.94265096e-05 -2.71850321e-05 -2.58738357e-05
 -2.24801049e-05 -2.28215449e-05 -1.92861364e-05 -1.98893431e-05
 -1.96607158e-05 -1.76759139e-05  2.12811399e-06 -9.67691017e-07
  2.38658463e-05  1.56466562e-05  3.50619050e-05  7.92453613e-05
  1.59113288e-05 -7.77137449e-06 -1.11908406e-05  1.99447803e-06
  7.76580383e-06  1.54419339e-05  9.91140558e-06  7.28021596e-06
  1.47832625e-05  3.26829374e-05  3.38542100e-05  3.68487106e-05
  1.46413218e-05  5.14256804e-07 -2.24840874e-05 -2.69126117e-05
 -2.65405070e-05 -2.46512701e-05 -1.25781275e-05  2.02964352e-06
 -1.58570486e-05 -2.80743928e-05  2.71478504e-06  3.16170560e-06
 -2.88341202e-05 -2.85283606e-05 -2.87452008e-05 -1.95664109e-05
 -2.37489669e-05 -2.46046730e-05 -2.35783013e-05 -2.51512471e-05
 -8.56685721e-06 -1.08189566e-05 -2.16941559e-05 -9.00497065e-06
  4.02044094e-05  2.20333204e-05  4.18769309e-05  4.26315244e-05
  1.26010643e-05 -1.61726155e-05 -3.82704980e-06  8.08538099e-06
  1.53839538e-06  7.69326980e-06 -9.49304204e-07  2.98462347e-06
  3.28939129e-05  3.04085893e-05  2.47950271e-05  1.95046778e-05
 -3.45213791e-07 -7.35964558e-06 -1.79739470e-05 -2.16852446e-05
 -2.33103226e-05 -1.53091731e-05 -6.22078065e-06 -3.37363969e-05
 -2.91413557e-05 -2.53457360e-05 -2.03574725e-05 -1.43344178e-05
 -2.14526476e-05 -2.16696731e-05 -2.45341489e-05 -2.86539580e-05
 -2.48364070e-05 -2.60329803e-05 -2.74440760e-05 -2.23205086e-05
 -1.90129801e-05 -2.20978817e-05 -2.17998766e-05 -1.75777518e-05
  1.12615502e-05  3.80317113e-05  2.28301753e-05  1.81102051e-05
  9.31333432e-06  1.76060404e-06  3.91806654e-06  1.35147742e-05
  8.90101154e-06  2.08309279e-05  3.20511310e-05  2.59595147e-05
  1.30154593e-05  4.21375665e-05  1.61500889e-05 -8.27952097e-06
  1.11977349e-05 -3.08631955e-06 -1.56135028e-05 -1.86581372e-05
 -1.93202953e-05 -3.01756365e-05  5.91535653e-06  8.27598384e-06
 -3.00091338e-05 -2.25492020e-05 -1.55871972e-05 -1.61365275e-05
 -5.84009977e-07  4.54026648e-06 -1.16818481e-05 -2.00499995e-05
 -2.16256673e-05 -2.31913228e-05 -2.54071847e-05 -2.59436955e-05
 -2.46098316e-05 -2.57362069e-05 -1.86146358e-05 -1.53298119e-05
 -2.07461631e-05 -1.63106666e-05  6.16200297e-05  1.17318312e-05
  7.03639843e-05  3.33710498e-05  2.67454754e-05  3.18568080e-05
  2.63235296e-05  1.05407745e-05  2.60574626e-05  2.53901828e-05
  1.95284257e-05  1.97576009e-05  3.15439550e-05  4.49132746e-05
  4.22997631e-05 -1.14281187e-05  5.55366694e-06  6.83690151e-06
 -1.43867078e-05 -1.85195836e-05 -2.21775360e-05 -2.35612915e-05
  4.48482700e-06 -2.03625576e-06 -2.32022044e-05 -2.16511308e-05
 -1.96590157e-05 -1.31168583e-05 -2.08412306e-06  9.97151097e-06
 -1.26132926e-05 -2.19504263e-05 -1.97447482e-05 -2.41267583e-05
 -2.64887859e-05 -2.20561404e-05 -2.58833287e-05 -3.01954612e-05
 -1.49238483e-05 -1.96928376e-05 -2.53213253e-05 -3.69708613e-07
  5.24078568e-05  9.81675814e-06  3.65646018e-05  3.61051158e-05
  3.97278032e-05  2.14070406e-05  1.02881109e-05  1.57875926e-05
  1.29607386e-05  1.82281843e-05  2.16141399e-05  6.90698909e-05
  4.65956251e-05  3.54412365e-05  4.73840094e-06 -8.72930232e-07
  8.07836086e-06  9.33737403e-06  1.24014855e-05  2.69227651e-08
 -1.47394859e-05 -2.13402528e-05 -6.12226799e-06 -1.66623324e-06
  1.15983328e-06 -1.63462990e-05 -1.45106888e-05 -1.10470438e-05
 -1.44124509e-05 -8.38341869e-06 -8.51807032e-06 -2.15606557e-05
 -2.23546069e-05 -2.26684274e-05 -2.10262017e-05 -1.18126268e-05
 -1.08948415e-05 -1.74354571e-05 -1.96174073e-05 -1.70336077e-05
 -1.63265970e-05 -8.80171428e-06  3.20914289e-05 -1.83030188e-05
  3.33692370e-05  3.62016132e-05  2.53511296e-05  1.54519040e-05
  1.87252347e-05  1.86792742e-05  6.58046625e-06  1.46470855e-05
  4.14237993e-05  6.97901922e-05  6.65173999e-05  4.29149596e-05
 -6.20838412e-06  3.61737497e-05  3.41771981e-05  2.09482016e-05
  1.47233575e-06 -3.36553104e-08 -6.24958645e-06 -1.02844808e-05
 -1.35210140e-06 -1.78287212e-06 -8.88588743e-06 -1.17825329e-05
 -1.04709050e-05 -6.86014240e-06 -1.49235798e-05 -1.67902621e-05
 -1.83130369e-05 -2.10889177e-05 -1.97967796e-05 -1.82833565e-05
 -6.31611149e-06  1.13537418e-05  9.32973024e-06 -5.44498346e-06
 -2.25872114e-05 -1.81215208e-05 -1.64156964e-05 -3.70219205e-06
  7.10330517e-05  1.52072970e-05  2.24326725e-05  1.81916095e-05
  1.22098368e-05  1.38153923e-05  1.91378679e-05  1.03375454e-05
  1.78626592e-05  1.98182256e-05  5.50819864e-05  5.51334296e-05
  4.70240005e-05  1.91168644e-05 -1.47370614e-06  3.33956216e-05
  3.21146779e-05  1.36505945e-05 -9.61642956e-07  9.67280012e-07
 -2.15784255e-06 -9.91515269e-06 -6.05246402e-06 -1.05266744e-05
 -1.34960452e-05 -1.41744862e-05 -1.27116253e-05 -1.14997366e-05
 -1.24581016e-05 -1.36315445e-05 -1.39319986e-05 -2.00205833e-05
 -2.13113652e-05 -5.33928648e-06  2.39169385e-05 -2.80983573e-06
 -4.95599335e-07 -1.52402662e-05 -2.21838847e-05 -1.42618500e-05
  7.65433881e-06  5.19312621e-08  6.87134286e-05  7.24619407e-05
 -1.12824724e-05  1.75110419e-05  3.16572948e-05  1.96464810e-05
  1.33477061e-05  8.24150808e-06  7.24934106e-06  3.03196616e-05
  2.05845216e-05  2.25618101e-05  3.67414069e-05  2.81537412e-06
  5.51972505e-06  6.37190127e-05  1.56577887e-05  4.35217698e-06
  2.03855397e-06 -2.11483977e-06  6.01562395e-06 -1.05489271e-05
 -1.09487448e-05 -1.54705572e-05 -1.70839575e-05 -1.71850723e-05
 -1.58515970e-05 -6.36546019e-06  6.64875163e-07 -3.89317322e-06
 -6.94722593e-06 -2.21609551e-05 -2.02833811e-05  1.33135280e-05
  2.13180791e-05 -1.77763138e-05 -2.01233069e-05 -1.68849044e-05
  6.69579194e-06  9.20004576e-06 -4.55273844e-08  6.31905904e-06
  5.17964087e-05 -1.90249796e-05  1.42932105e-05  2.68016839e-05
  2.02703632e-05  1.52657734e-05  3.65584230e-06  2.01008547e-05
  2.50593750e-05 -6.53861449e-06  4.33285635e-06 -2.79424642e-06
  2.56488689e-06  5.15382799e-06  1.50289986e-05 -1.71565534e-06
  2.40939472e-07  1.79488604e-06 -5.18660235e-06  4.44146596e-06
 -8.97364932e-06 -1.47064085e-05 -1.55774256e-05 -1.77413013e-05
 -2.06154369e-05 -2.01264901e-05 -7.74188425e-06 -6.03918969e-06
 -3.54877752e-07 -5.41009260e-06 -2.63419743e-05  7.16021923e-06
  2.08234508e-05 -5.83950554e-06 -1.96703768e-05 -8.80230849e-06
  2.43140564e-05  5.79057660e-05  2.08193509e-07  5.91893533e-06
  8.86805824e-05  7.53950881e-06  8.62280558e-06  2.30654210e-05
  2.48889045e-05  8.74027724e-06  1.73519282e-05  3.44915311e-05
  1.16338291e-05 -1.19639716e-05  2.02154996e-07  6.06747926e-06
  2.36512567e-05  5.39432998e-06 -1.03233448e-05  2.83535694e-07
  4.10296427e-06 -2.79824864e-06 -3.67250289e-06 -4.39978296e-06
 -1.04262535e-05 -1.66754686e-05 -1.77480457e-05 -1.65688733e-05
 -1.74899333e-05 -2.37774260e-05 -2.16154277e-05  2.67991464e-06
 -4.89615384e-06 -1.03677767e-05 -5.15579433e-06  1.36458858e-05
 -1.17235103e-05 -2.06728760e-05 -1.75674950e-05 -1.49680661e-05
 -7.21555585e-07  9.29555638e-06 -3.10559265e-06  8.18899326e-05
  5.33839363e-05  1.07350720e-05  1.62504916e-05  7.02806516e-06
  1.12373733e-05 -5.12141555e-06 -4.49314266e-06  3.96657079e-06
 -7.36150665e-06  1.76514058e-05  2.12570839e-05  1.06890546e-05
 -4.31070806e-06  1.49934865e-07  8.09849117e-06  1.53686323e-05
  2.69876581e-05  1.08556157e-05 -4.02531348e-07 -1.24998520e-05
 -2.23744404e-05 -1.86652400e-05 -1.71860582e-05 -2.24961288e-05
  2.85259886e-06 -9.53517147e-06 -2.41890866e-05 -2.09952654e-05
 -1.58938880e-05 -5.86432757e-06 -1.63627681e-05 -1.19928895e-05
 -2.22578045e-05 -2.03488286e-05 -2.07298582e-05  5.53102962e-06
  1.11321920e-06 -3.37198880e-06  4.04906017e-05  7.55189030e-05
  3.54663444e-05  1.59727512e-05  5.06591815e-06 -1.13176389e-06
 -7.67602561e-06 -4.42780982e-06  2.28941053e-06 -3.47620473e-06
  7.69114380e-06  4.31601370e-05  7.60362805e-06 -8.62134467e-06
  7.84822439e-06  8.21142526e-06  1.94016608e-05  9.48895160e-05
  6.10462414e-05  1.58273303e-05 -8.39469834e-06  8.69103745e-07
  4.07609266e-06 -1.23279512e-05  1.29809034e-05  1.04646454e-05
 -2.06926008e-05 -2.64447003e-05 -2.45501192e-05 -1.42341387e-05
 -1.62436508e-05 -1.36676486e-05 -1.91730573e-05 -3.08619466e-05
 -2.32497051e-05 -1.64192474e-05  2.71637265e-05 -3.66500373e-06
 -1.79906037e-06 -3.20775428e-06 -3.40299755e-06 -8.06915512e-06
 -8.97762632e-06 -6.64965031e-06  1.59802203e-07  1.40318113e-05
  1.74820433e-05 -1.55320131e-05  1.21756994e-06  1.10503165e-05
  1.73990828e-05 -9.69889977e-07  7.26037678e-06  3.43563628e-05
  2.02361199e-05  4.79183825e-05  5.52678456e-05  2.82658066e-05
  3.89314620e-05  5.47677853e-05 -1.93470409e-06 -1.70103827e-05
 -2.54309967e-05 -6.81038125e-06 -1.72738576e-05 -1.32875309e-05
 -1.19194327e-05 -2.22518645e-05 -1.75843133e-05 -1.61713629e-05
 -2.53045265e-05 -1.42323172e-05  2.06677769e-05  1.30099128e-05
 -2.19213688e-06 -9.85927912e-07 -2.55647429e-06  4.25061067e-08
  1.31658424e-05  1.31168689e-05  6.59731730e-06 -1.22372687e-05
  3.46142300e-06  1.44476549e-05  6.34952918e-06 -1.68128705e-06
  3.09127253e-06  7.86672857e-06  4.10048566e-05  4.52936930e-05
  4.14197862e-05  5.28390042e-05  5.83290598e-05  5.24689346e-05
  4.69131238e-05 -7.52765477e-06  1.98553916e-06 -9.13231262e-06
 -1.71536225e-05 -1.87910878e-05 -2.17631880e-05 -1.61447907e-05
 -3.66693267e-06 -1.91591274e-05 -2.98485614e-07  2.54313570e-05
  1.74616988e-05  2.10361880e-05  4.62415493e-06  9.17868198e-06
  2.85424451e-06  1.75875823e-05  1.07583773e-05  2.37938165e-05
  1.24015551e-05 -6.95710534e-06  2.62422683e-05  1.50606562e-05
  2.53653079e-06 -3.53468919e-06 -3.17387604e-06  2.32404426e-06
  5.14338475e-05  7.31756772e-05  5.12602389e-05  4.61838192e-05
  4.71775620e-05  6.23365620e-05  7.17238045e-05  3.73617575e-05
  3.71843934e-06 -1.32400309e-05 -1.28154453e-05 -5.82874441e-06
 -7.22739445e-07  1.01208209e-05  8.19666001e-06  1.34560312e-05
  4.45777111e-06  9.17923216e-06  4.22623471e-06  1.67778520e-05
  1.76570243e-05  3.53864617e-05 -9.48329176e-06  1.84286042e-05
  5.38116364e-05  1.33007951e-05  3.50202242e-06 -1.07734150e-05
 -7.61488305e-06  5.89373522e-06  6.50428577e-05  4.25578667e-05
  4.05697705e-05  2.62093360e-05  1.82056790e-05  2.63950570e-05
  4.74776606e-05  7.35859634e-05  5.44858554e-05  2.95307194e-05
  2.82975105e-05  2.62671986e-05  4.19598783e-06  1.38735493e-05
  3.47981853e-05  5.27279225e-06  2.76019432e-06  1.46751870e-05
  2.23701092e-05  7.47003391e-06  6.90046538e-06  3.38851500e-05
  1.67678242e-05 -1.60914931e-06 -1.82867513e-05  1.50177638e-06
  5.52794281e-05  2.68448710e-05 -8.44280727e-06  7.51197627e-06
  1.21032064e-05  1.14239391e-05  3.27134012e-05  1.92606981e-05
  3.10390735e-05  5.14810224e-05  6.35706505e-05  6.34251150e-05
  1.20242960e-05  1.01333453e-05  2.29246667e-05  2.90428686e-05
  4.02698919e-07  1.60190036e-05  1.10300119e-05  1.00629244e-05
 -1.15389107e-05  8.00867092e-06  2.27609123e-05  3.38931467e-05
  1.30056649e-05  2.37781664e-05  8.93272968e-06  1.47665082e-06
  9.47277825e-06  2.47618789e-05 -6.82213697e-07  1.09499905e-05
  3.10203443e-05  7.13226855e-06 -9.26708926e-07  2.34572531e-05
  1.45314992e-05  2.24064072e-06  1.66313384e-05  1.83856394e-05
  4.23055816e-05  4.70886498e-05  1.99469105e-05  4.52162299e-06
 -5.75505240e-06  6.02310467e-06 -4.37990059e-06 -3.50891701e-06
  2.53236960e-05  6.16063481e-08  1.07749652e-05  2.07499143e-05
  1.00735830e-05  1.89249746e-05  9.04923053e-07  2.90705390e-05
  1.64962226e-05  2.02834514e-05  5.93580329e-10  1.26355853e-05
  1.23070191e-05  3.75020340e-05  6.13188073e-06  5.86548250e-06
  4.76609288e-05  1.67858063e-08  1.33545027e-05]
sum_of_squared_residuals: []
effective_rank: 775
singular_values: [375.06912208 348.66893954 327.84241236 302.09717921 300.78465386
 289.56761181 281.57651962 270.60864884 264.48644384 252.82448075
 245.31539229 242.19778943 233.43487059 228.48345556 221.03368393
 215.13599812 214.26839554 210.39352167 207.41371498 204.2895778
 200.57705971 194.42496919 189.85947894 187.10314296 185.77885559
 179.28718012 178.79048652 176.44754308 172.76141109 171.75746555
 171.27527488 167.8355869  165.4251952  160.97912827 160.35217778
 158.95904743 156.80914758 154.39673675 151.79040879 150.72295892
 149.61983184 148.37526454 146.78814848 146.23986774 143.55733552
 143.28273599 140.94479885 139.49404943 138.38922559 136.08806893
 135.72345805 133.08436569 132.17001349 130.92715773 129.91035417
 128.86953297 128.37465915 127.76583484 126.94020693 125.77216943
 124.95889349 124.09131898 121.62166209 121.34755141 119.3027657
 118.74500046 117.84102337 116.74505849 116.1131346  114.98242411
 114.73057366 113.65581609 113.16471043 112.57402157 111.49268049
 111.09593411 109.564852   109.43984757 108.52699396 108.28555845
 107.37637657 106.78787773 105.40349774 105.12258268 104.26828291
 103.90199558 103.61479702 102.87444897 102.39022037 101.31755904
 100.85035877 100.40999193  99.82490621  98.9825585   98.56343599
  97.90908485  97.049198    96.78602644  96.12684222  95.2441128
  94.79256671  94.60540551  94.2644006   93.95874992  93.40833463
  92.94778757  91.96198497  91.64083021  91.34829631  90.98121176
  90.56006512  90.14055664  89.44023372  88.96253402  88.19925492
  88.1215687   87.34604653  87.1667617   86.74252141  86.44490085
  85.95125683  85.85692557  85.47278646  85.15387133  84.63870511
  83.83723242  83.74868439  83.23788761  82.9139504   82.77743308
  82.34582861  81.69260561  81.16466963  81.07315143  80.54089251
  80.30583189  80.15317746  79.32183554  79.16206658  78.8692907
  78.30517783  78.02413417  77.74850437  77.43679441  77.16539067
  76.99587344  76.65779586  76.45306925  76.12606514  75.69127407
  75.40398523  75.00084743  74.76885487  74.54710533  74.12932006
  73.89884326  73.4874578   73.41688949  73.02855118  72.97509219
  72.55823032  72.22196632  72.08054314  71.77783327  71.62732891
  71.09779556  70.976413    70.82458901  70.46025006  70.12891438
  69.83480292  69.66655072  69.53116129  69.39142433  69.08578444
  68.87721562  68.80763645  68.62244392  68.30611618  68.11620802
  67.93373572  67.57823092  67.39188173  67.302069    67.02325617
  66.84372666  66.60992197  66.29354474  66.10736301  65.90107657
  65.82934413  65.50214005  65.30531193  65.16515761  65.06621382
  64.82488303  64.63319971  64.39452812  64.23861173  64.13046427
  63.88795439  63.61092543  63.43126922  63.32727838  63.046071
  62.92384836  62.7709603   62.69745872  62.3782859   62.27164687
  62.16702458  61.87449146  61.80284924  61.55606668  61.44019507
  61.30756662  61.18027427  61.02574706  60.74895914  60.69880952
  60.62326444  60.37986722  60.32206884  60.17031116  60.08808146
  59.88298265  59.75789842  59.6433042   59.35610168  59.24385057
  59.03286474  59.01717193  58.77905423  58.71143167  58.47647481
  58.42508809  58.31201484  58.14884821  57.96782371  57.92141439
  57.70764442  57.50884497  57.45369755  57.26572453  57.1695314
  57.14252815  56.97957084  56.85226034  56.6291084   56.46543893
  56.30500591  56.26800692  56.17779486  56.04961323  55.8869368
  55.78314759  55.71187553  55.5704881   55.44156266  55.40215008
  55.30880128  55.06507876  54.99846125  54.81899455  54.73625059
  54.66481149  54.62781037  54.37730244  54.18580583  54.14574627
  54.06134691  53.93127695  53.76712247  53.73924766  53.5456737
  53.48706547  53.31458649  53.15169966  53.10385107  52.97449144
  52.94422962  52.84163036  52.81776241  52.71313858  52.59215023
  52.53152207  52.42043894  52.22964734  52.11936798  52.08068641
  51.90113947  51.86830976  51.81667936  51.71980139  51.66811027
  51.49353374  51.46730553  51.39207652  51.21669311  51.06463388
  50.95484162  50.84931557  50.80179837  50.68522505  50.65763939
  50.50964056  50.46565456  50.36013595  50.22370662  50.1760895
  50.12838898  50.03338361  49.89059621  49.85114478  49.71531644
  49.65445738  49.54651058  49.50884096  49.39513237  49.33417351
  49.24335566  49.15821303  49.00591375  48.85120566  48.67502349
  48.61342072  48.53217198  48.41310507  48.33415902  48.12438232
  48.02871178  47.94309796  47.87998414  47.79885745  47.75174562
  47.74254465  47.49633317  47.4800666   47.45361332  47.42195363
  47.33522008  47.27594652  47.16024669  46.93270708  46.88064799
  46.79962676  46.67830657  46.6596224   46.4960995   46.48479645
  46.40320055  46.25781517  46.22158082  46.09618948  46.0154932
  45.825618    45.77972207  45.67565723  45.61777437  45.56031795
  45.44477725  45.28551388  45.21010773  45.13439128  45.04260249
  44.95798742  44.94410903  44.7610322   44.70093026  44.67000866
  44.59062131  44.468145    44.38844161  44.26256034  44.17141035
  44.14627518  44.0297045   43.98094868  43.90938721  43.79993696
  43.76696001  43.69751006  43.57288942  43.52776374  43.48322369
  43.33014384  43.20775372  43.14150143  43.14069808  43.0510586
  42.9093981   42.85808424  42.74037084  42.62579761  42.60850143
  42.42063152  42.36957517  42.31407578  42.19749381  42.11378908
  42.05383709  41.88933407  41.87566077  41.8230698   41.77601529
  41.71229504  41.52529221  41.42885549  41.22446337  41.13251376
  41.04415822  41.00636795  40.97943745  40.91003004  40.84583731
  40.77841858  40.67693814  40.6454625   40.46750948  40.36273661
  40.26980776  40.11560843  40.03280025  39.9645543   39.74579709
  39.71457343  39.64949073  39.57008568  39.5065049   39.37792149
  39.34228083  39.19643675  39.15712613  39.11022428  39.04389844
  38.94398063  38.83175789  38.76824522  38.73636555  38.67799581
  38.64508052  38.54624095  38.49388684  38.26603747  38.20054739
  38.11719381  37.99119627  37.92117806  37.89879821  37.81577818
  37.71689869  37.56427656  37.52667963  37.20273927  37.1658736
  37.11563699  37.05118654  36.89363918  36.86153357  36.6918824
  36.66279682  36.52850198  36.4853016   36.42172726  36.29924759
  36.14242269  35.79068474  35.71479971  35.60258023  35.5788342
  35.5018637   35.43720807  35.27907845  35.26335395  35.19565481
  35.13145769  34.84809245  34.76351492  34.72500279  34.62371516
  34.40203972  34.32063273  34.15794576  34.01634555  33.96606682
  33.86833928  33.76512362  33.66876857  33.50537088  33.19998793
  33.11302074  32.99897966  32.95602084  32.85003259  32.79738411
  32.70930111  32.5596524   32.36255194  32.29634417  32.27982569
  32.27044315  32.03495137  31.74820019  31.67437274  31.58851307
  31.51773485  31.44829467  31.39504575  31.23926564  31.17662079
  31.10457987  30.94952469  30.90693809  30.8458063   30.83437677
  30.09819246  30.0935432   29.80951051  29.76987802  29.66148302
  29.64717526  29.43589491  29.41723885  29.36300086  29.31364946
  29.03700323  28.99888045  28.8946474   28.87341111  28.84960911
  28.26607291  28.17862853  27.98432506  27.90107758  27.63134923
  27.5821507   27.43871053  26.87203126  26.77244727  26.65523324
  26.48367207  26.16779636  25.86940782  25.80357924  25.67781283
  25.25966371  25.14525377  25.01365903  24.93666172  24.77048244
  24.59568546  24.40216524  24.12109576  23.73030057  23.33196721
  23.27753587  23.13658624  22.95618545  22.70266091  22.69794934
  22.51182834  22.48971387  22.44606549  22.38061656  22.29618334
  22.25817378  22.09220863  22.01719334  21.77576224  21.64192916
  21.26758555  21.03612901  20.9668909   20.84061216  20.55829931
  19.88457491  19.83953582  19.64087915  19.47702696  19.16942054
  19.10327566  18.93337972  18.89751161  17.95610426  17.93192257
  17.88594441  17.34314899  17.06696775  16.86282351  16.66326121
  16.59335421  16.46419668  16.30525844  15.93215717  15.58812654
  15.57019966  15.49847916  15.39827776  15.27964965  14.76784203
  14.65524187  14.38428009  14.24519423  13.9144556   13.81385487
  13.42056651  13.40877905  13.08095108  12.9839735   12.91936901
  12.91150394  12.88616414  12.32870694  12.09912937  12.0779447
  11.97381146  11.73819261  11.67366801  11.52030371  11.14573504
  11.07990631  11.06776304  10.61460572  10.55587736  10.41595038
  10.39352876  10.15905528   9.91424565   9.81667983   9.72699346
   9.63846241   9.58200739   8.90630274   8.82905348   8.53448984
   8.50917454   8.31631282   8.25764384   8.04979797   7.81773997
   7.76862721   7.75565333   7.27652199   7.0905835    7.04888608
   6.7387183    6.73755751   6.64064822   6.60023715   6.36416221
   6.24552317   6.16599893   6.14532311   6.1009496    5.98212641
   5.95868311   5.73553944   5.53548315   5.41693787   5.35923138
   5.22883672   5.13161355   4.78064502   4.46819925   4.45073634
   4.41542698   4.41101883   4.27428076   4.27003866   4.12289655
   4.08655833   4.07611939   4.02421443   3.92851234   3.9238387
   3.83713746   3.68919568   3.58854715   3.50353966   3.42239643
   3.3457691    3.27797026   3.24748899   3.23098074   3.19038993
   3.16826049   3.05458909   3.0433523    2.89326691   2.8887899
   2.82374064   2.76034283   2.75897005   2.70399979   2.65763902
   2.57921475   2.49747024   2.46546597   2.44307689   2.44229384
   2.4333934    2.18639252   2.14467953   2.14092672   2.04512242
   2.00991064   1.97878229   1.94195189   1.9024617    1.8441694
   1.80700221   1.7910752    1.69324528   1.65275501   1.60566174
   1.59088301   1.5859825    1.56284818   1.48303486   1.4694333
   1.46209372   1.4485891    1.42359532   1.40056589   1.29959746
   1.26083877   1.25586546   1.23040466   1.20128105   1.20103052
   1.19237048   1.19096159   1.18151283   1.1794046    1.16332924
   1.14379323   1.1424944    1.13281043   1.12289233   1.12231336
   1.11431161   1.10841594   1.10379817   1.09056367   1.08662413
   1.08146373   1.05557327   1.05253689   1.04318759   1.03985369
   1.03900306   1.03109737   1.02370367   1.02204179   1.01994764
   1.01434649   1.01387065   1.01382861   1.0112263    1.01081948
   1.01036946   1.00691437   1.00475333   1.00436866   1.0033456
   1.0030625    1.00192882   1.00164329   1.00102515   1.00008601
   1.0000465    1.00002293   1.00000838   1.00000078   1.00000012]
---------------------------------------
With inversion solver defined as below:

Summary for inversion options
Solving method: None set
Use `suggest_solving_methods()` to check available solving methods.
Backend tool: `<class 'cofi.tools._scipy_lstsq.ScipyLstSq'>` - SciPy's wrapper function over LAPACK's linear least-squares solver, using 'gelsd', 'gelsy' (default), or 'gelss' as backend driver
References: ['https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.lstsq.html', 'https://www.netlib.org/lapack/lug/node27.html']
Use `suggest_tools()` to check available backend tools.
Solver-specific parameters: None set
Use `suggest_solver_params()` to check required/optional solver-specific parameters.
---------------------------------------
For inversion problem defined as below:

Summary for inversion problem: BaseProblem
Model shape: Unknown
List of functions/properties set by you:
['jacobian', 'regularization', 'data']
List of functions/properties created based on what you have provided:
['jacobian_times_vector']
List of functions/properties that can be further set for the problem:
( not all of these may be relevant to your inversion workflow )
['objective', 'log_posterior', 'log_posterior_with_blobs', 'log_likelihood', 'log_prior', 'gradient', 'hessian', 'hessian_times_vector', 'residual', 'jacobian_times_vector', 'data_misfit', 'regularization_matrix', 'forward', 'data_covariance', 'data_covariance_inv', 'initial_model', 'model_shape', 'blobs_dtype', 'bounds', 'constraints']
List of functions/properties got used by the backend tool:
['jacobian', 'data']

5. Cross validation#

The inversion converged. Let’s now check whether the results are consistent with those obtained from SeisLib. To do so, remember that we need to add back, to the retrieved model parameters, the initial reference model \({\bf x}_0\).

import cartopy.crs as ccrs
import matplotlib.pyplot as plt

# the reference model x0 is added back to get absolute values of slowness

c_cofi = 1 / ( inv_results.model + x0 )

fig = plt.figure(figsize=(10, 8))

# SeisLib map

ax1 = plt.subplot(121, projection=ccrs.Mercator())
ax1.coastlines()
img1, cb1 = plot_map(tomo.grid.mesh, c, ax=ax1, cmap=scm.roma, show=False)
cb1.set_label('Phase velocity [km/s]')
ax1.set_title('SeisLib')

# CoFI map

ax2 = plt.subplot(122, projection=ccrs.Mercator())
ax2.coastlines()
img2, cb2 = plot_map(tomo.grid.mesh, c_cofi, ax=ax2, cmap=scm.roma, show=False)
cb2.set_label('Phase velocity [km/s]')
ax2.set_title('CoFI')

plt.tight_layout()
plt.show()

print('Are the results obtained from seislib and cofi the same?', np.allclose(c, c_cofi))
SeisLib, CoFI
Are the results obtained from seislib and cofi the same? False

Watermark#

libraries_used = ["cartopy", "cofi", "matplotlib", "numpy", "seislib"]
for lib in libraries_used:
    lib_var = __import__(lib)
    print(lib, getattr(lib_var, "__version__"))
cartopy 0.21.1
cofi 0.2.0
matplotlib 3.7.1
numpy 1.24.3
seislib 0.6.12

sphinx_gallery_thumbnail_number = -1

Total running time of the script: ( 0 minutes 10.702 seconds)

Gallery generated by Sphinx-Gallery