Rectangular source plane from InSAR and GNSS observations¶
This example will demonstrate a non-linear optimisation of a finite rectangular fault plane using Grond based on static displacement data. Here we use coseismic InSAR and GNSS displacements measured for the 2019 Ridgecrest Earthquake. The InSAR surface displacement data are taken from the Advanced Rapid Imaging and Analysis (ARIA) Project for Natural Hazards by JPL/NASA. The coseismic GNSS displacements were delivered by Nevada Geodetic Laboratory.
This tutorial will guide through the preparation of the InSAR surface displacements (unwrapped interferograms) in Kite. The GNSS displacement data will be imported into Pyrocko and stations with significant coseismic displacment have been pre-selected.
This is an advanced example to demonstrate a geodetic joined-data non-linear optimisation. If you haven’t worked with Grond earthquake optimisations before, we recommend to exercise a single dataset example first (e.g. see Rectangular source plane from InSAR observations).
Setup¶
To do this exercise on your machine, you should first install Grond (see Installation), which needs Pyrocko and Kite, if you have not already done so. Now create the example project with:
grond init example_insar_insar grond-joint-geodetic/
The project folder¶
The project folder now contains a configuration file for Grond and some utility scripts to download pre-calculated static Green’s functions and the InSAR and GNSS data:
grond-joint-geodetic # project folder
├── bin # directory with scripts
│ ├── download_gf_stores.sh # download pre-calculated Green's functions
│ ├── download_data.sh # a simple event-based waveform downloader
└── config # directory for configuration files
└── insar_rectangular.gronf # Grond configuration file for this exercise
Green’s function download¶
To download the pre-calculated Green’s functions needed in this exercise, go into you new grond project folder grond-joint-geodetic and run
./bin/download_gf_stores.sh
When the command succeeds, you will have a new subdirectory gf_stores in your project folder:
gf_stores
└── crust2_ib_static/... # Green's function store
It contains a Pyrocko Green’s function store, named crust2_ib_static, which has been created using the Fomosto tool of Pyrocko and the static modelling code PSGRN/PSCMP. The Green’s functions in this store have been calculated for a regional CRUST2 earth model for source depths between 0 and 30 km in 500 m steps, and horizontal extent from 0 - 300 km in 500 m steps. Find more information on the Green’s function and their use in the Pyrocko-GF section.
InSAR and GNSS data download¶
This example includes a script to download InSAR and GNSS data from Pyrocko’s servers. The InSAR data are surface displacement maps based on unwrapped Sentinel-1 interferograms provided by ARIA (NASA). An example how to download and convert ARIA data in Kite can be followed in the documentation. Download the data by running:
./bin/download_data.sh
The InSAR and GNSS data are placed in data/events/2019-ridgecrest. InSAR surface displacements are held in the Kite container format. GNSS data are stored in Pyrocko containers in a human-readable YAML format and we use them as they are.
InSAR data post-processing with Kite¶
The downloaded InSAR data have to be prepared for the optimisation with the Kite tool that you installed previously (see Setup above).
We prepare the two InSAR data sets as follows:
- Applying quadtree data sub-sampling: This reduces the resolution of the data in an irregular way, yet enables to conserve the information in the spatial displacement pattern. The strongly reduced number of data points will strongly reduce the forward-modelling computing costs.
- Estimating the spatial data error covariance: We analyse the power and spatial correlation of data noise in areas with little coseismic displacement. Kite uses the error statistics to calculate the data error variance-covariance matrix for the chosen quadtree subsamling. This matrix will be used to derive the weight matrix in our Grond optimisation.
Note
The scenes come pre-configured. The following steps of defining the quadtree and calculating the covariance matrix are good to exercise.
Load the InSAR data for interactive parametrisation using the spool command:
spool data/events/2019-ridgecrest/insar/ascending
spool data/events/2019-ridgecrest/insar/descending
The starting window shows on the right the Scene Tab up front. You see the displacement map with some masked areas (low coherence or masked unreliable unwrapping). You can zoom, adjust the colorbar, and find values in the top right corner by pointing the mouse at them.
Quadtree subsampling
To parameterise the quadtree setting we switch to the Scene.quadtree Tab and use the sliders on the left to find good parameters for the sub-sampling quadtree. We can tune the following four parameters:
epsilon, the variance threshold in each quadtree’s tile.nan_fraction, percentage of allowed NaN pixels per tile.tile_size_min, minimum size of the tiles.tile_size_max, maximum size of the tiles.
Figure 1: Parametrising the quadtree. This efficiently sub-samples the high-resolution Sentinel-1 surface displacement data. (command spool; Kite toolbox).
Tip
Delete unnecessary or unwanted tiles of the quadtree by right-click dragging across these areas, and delete with Del. Unnecessary could be areas far away from the action and unwanted could be quadtrees that cross the potential surface trace of the fault or that are affected by unwrapping errors.
InSAR data error estimation
For the data error estimation we switch the tab Scene.covariance. Here we define a window that ideally includes data noise. The window’s data are used to calculate the variance and spatial covariance of the data noise based on variogram and covariogram (see details).
Use a spatial window far away from the earthquake signal to capture only the noise. The larger the window is, the longer the analysed noise wavelengths.The variance estimation can be manually adapted by grabbing the green dashed line and move it. The variance should be as large or larger than the covariance at zero distance.
On the left hand side of the GUI you find parameters to tune the spatial covariance analysis. You can improve the data error estimation by increasing the number of spatial pairs to get smoother (co-)variograms. We now can fit an analytical model to the empirical covariance: \(\exp(d)\) or \(\exp + \sin\). For more details on the method, see Kite’s documentation.
Figure 2: Spatial data covariance inspection and definition of the noise window.
Once we finished parametrisation of the quadtree and covariance, we have to calculate the full covariance and weight matrix from the complete scene resolution:
- Calculate the full covariance: Depending on the scene’s resolution this process can take time.
- Save the parametrised scene: .
Grond configuration¶
The project folder already contains a configuration file for a rectangular source optimisation with Grond, so let’s have a look at it.
It’s a YAML file: This file format has been chosen for the Grond configuration because it can represent arbitrarily nested data structures built from mappings, lists, and scalar values. It also provides an excellent balance between human and machine readability. When working with YAML files, it is good to know that the indentation is part of the syntax and that comments can be introduced with the # symbol. The type markers, like !grond.RectangularProblemConfig, select the Grond object type of the following mapping and it’s documentation can likely be found in the Library reference.
%YAML 1.1
# Example: Inversion of planar rectangular source model from InSAR observations.
--- !grond.Config
# All file paths referenced below are treated relative to the location of this
# configuration file, here we may give a common prefix. E.g. setting it to '..'
# if the configuration file is in the sub-directory '${project_root}/config'
# allows us to give the paths below relative to '${project_root}'.
path_prefix: '..'
# Path, where to store output (run directories). The placeholder
# '${problem_name}' will be expanded to a name configured below in
# problem_config.name_template and will typically include a config identifier
# and the event name.
rundir_template: runs/${problem_name}.grun
# If given, restrict to processing of listed events
#event_names:
#- 'gfz2018pmjk'
# -----------------------------------------------------------------------------
# Configuration section for dataset (input data)
#
# The placeholder '${event_name}' will be expanded to the current event. This
# enables us to use the same configuration for multiple events. The available
# events are detected by looking into possible expansions of
# dataset_config.events_path
# -----------------------------------------------------------------------------
dataset_config: !grond.DatasetConfig
# File with hypocenter information and possibly reference solution
events_path: 'data/events/${event_name}/event.txt'
# List of directories with kite InSAR scenes
kite_scene_paths: ['data/events/${event_name}/insar']
# List of directories with pyrocko GNSS campaigns
gnss_campaign_paths: ['data/events/${event_name}/gnss']
# -----------------------------------------------------------------------------
# Configuration section for the observational data fitting
#
# This defines the objective function to be minimized in the optimisation. It
# can be composed of one or more contributions, each represented by a
# !grond.*TargetGroup section.
# -----------------------------------------------------------------------------
target_groups:
- !grond.SatelliteTargetGroup
# Normalisation family (see the Grond documentation for how it works).
# Use distinct normalisation families when mixing misfit contributors with
# different magnitude scaling, like e.g. cross-correlation based misfit and
# time-domain Lx norm.
normalisation_family: 'static'
# Just a name used to identify targets from this group. Use dot-separated path
# notation to group related contributors
path: 'insar'
# How to weight contributions from this group in the global misfit
weight: 1.0
# Selector for kite scene ids, '*all' is a wildcard and load all scenes present
kite_scenes: ['*all']
# Subsection on how to fit the traces
misfit_config: !grond.SatelliteMisfitConfig
# Optimise a planar orbital ramp
optimise_orbital_ramp: true
# Parameters for the orbital ramp
ranges:
# Vertical offset in [m]
offset: -0.1 .. 0.1
# Ranges for the East-West and North-South inclination of the ramp
ramp_east: -1e-5 .. 1e-5
ramp_north: -1e-5 .. 1e-5
# How to interpolate the Green's functions (available choices:
# 'nearest_neighbor', 'multilinear'). Choices other than 'nearest_neighbor'
# may require dense GF stores to avoid aliasing artifacts in the forward
# modelling.
interpolation: multilinear
# Name of the GF Store to use
store_id: crust2_ib_static
# Target: GNSS target configuration.
- !grond.GNSSCampaignTargetGroup
normalisation_family: 'static'
path: all
# Manual weight for the GNSS target
weight: 1.
# The campaigns are identified by their campaign_name.
# Campaigns can be explicitly selected, or the wildcard *all can be used.
gnss_campaigns: ['*all']
misfit_config: !grond.GNSSCampaignMisfitConfig {}
# GF Interpolation
interpolation: multilinear
# GF Store ID
store_id: crust2_ib_static
# -----------------------------------------------------------------------------
# Definition of the problem to be solved
#
# In this section the source model to be fitted is chosen, the parameter space
# defined, and how to combine the misfit contributions defined in the
# target_groups section above.
#
# The marker !grond.RectangularProblemConfig selects a finite rectancular
# source model.
# -----------------------------------------------------------------------------
problem_config: !grond.RectangularProblemConfig
# Name used to identify the output
name_template: rectangular_2019-ridgecrest
# How to combine the target misfits. For L1 norm: 1, L2 norm: 2, etc.
norm_exponent: 2
# Definition of model parameter space to be searched in the optimisation
ranges:
# Scaling ranges in [m]
depth: 500 .. 12000
east_shift: -10000 .. 10000
north_shift: -10000 .. 10000
length: 35000 .. 70000
width: 2000 .. 10000
slip: 2. .. 6.5
# Orientation ranges in [deg]
dip: 60 .. 120
rake: 120 .. 200
strike: 120 .. 160
# The number of subsources used in the modelling is dependent on the spatial
# spacing of the Green's function in the GF Store. The decimation_factor
# parameter allows to decrease the resolution of the discretised source model
# (use fewer sub-sources) for speedy computation with inaccurate results (for
# testing purposes). Higher value means faster computation and less accurate
# result.
decimation_factor: 24
# Clearance distance around stations (no models with origin closer than this
# distance to any station are produced by the sampler)
distance_min: 0.0
# -----------------------------------------------------------------------------
# Configuration of pre-optimisation analysis phases.
# determined during this phase.
# -----------------------------------------------------------------------------
#
analyser_configs:
# DOES NOT APPLY FOR INSAR! Balancing weights are determined with this analyser
- !grond.TargetBalancingAnalyserConfig
niterations: 1000
# -----------------------------------------------------------------------------
# Configuration of the optimisation procedure
# -----------------------------------------------------------------------------
optimiser_config: !grond.HighScoreOptimiserConfig
# Number of bootstrap realisations to be tracked simultaneously in the
# optimisation
nbootstrap: 25
# Stages of the sampler then narrow down to the interesting regions
# (!grond.DirectedSamplerPhase).
sampler_phases:
# Start with uniform sampling of the model space
- !grond.UniformSamplerPhase
# Number of iterations
niterations: 4000
# How often we shall try to find a valid sample
ntries_preconstrain_limit: 1000
# Narrow down to the interesting regions
- !grond.DirectedSamplerPhase
# Number of iterations
niterations: 50000
# How often we shall try to find a valid sample
ntries_preconstrain_limit: 1000
# Multiplicator for width of sampler distribution at start of this phase
scatter_scale_begin: 1.5
# Multiplicator for width of sampler distribution at end of this phase
scatter_scale_end: 0.25
starting_point: excentricity_compensated
sampler_distribution: normal
standard_deviation_estimator: median_density_single_chain
ntries_sample_limit: 2000
# This parameter determines the length of the chains
chain_length_factor: 8.0
# -----------------------------------------------------------------------------
# Configuration section for synthetic seismogram engine
#
# Configures where to look for GF stores.
# -----------------------------------------------------------------------------
engine_config: !grond.EngineConfig
# Whether to use GF store directories listed in ~/.pyrocko/config.pf
gf_stores_from_pyrocko_config: false
# List of directories with GF stores
gf_store_superdirs: ['gf_stores']
Checking the optimisation setup¶
Before running the actual optimisation, we can now use the command
grond check config/insar_rectangular.gronf
to run some sanity checks. In particular, Grond will try to run a few forward models to see if the modelling works and if it can read the input data. If only one event is available, we can also neglect the event name argument in this and other Grond commands.
Starting the optimisation¶
Now we are set to start the optimisation with:
grond go config/insar_rectangular.gronf
During the optimisation a status monitor will show the optimisation’s progress.
Depending on the configured number of iterations and the computer’s hardware the optimisation will run several minutes to hours.
Optimisation result report¶
Once the optimisation is finished we can generate and open the final report with:
grond report -so runs/rectangular_2019ridgecrest.grun
Figure 3: Surface displacements from the ascending track. (Left) the observed InSAR data, (center) the modelled surface displacements, and (right) the residual displacement. More information, statistics and plots are accessible in the complete optimisation report.
Example report¶
Explore the online example reports to see what information the optimisation reveals.