Note
Go to the end to download the full example code
Surface-Wave Tomography#
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
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
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
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()
<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]')
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).
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))
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)