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:

  1. 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.
  2. 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:

  1. epsilon, the variance threshold in each quadtree’s tile.
  2. nan_fraction, percentage of allowed NaN pixels per tile.
  3. tile_size_min, minimum size of the tiles.
  4. tile_size_max, maximum size of the tiles.
../../_images/example_spool-ridgecrest-quadtree.png

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.

../../_images/example_spool-ridgecrest-covariance.png

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:

  1. Calculate the full covariance: Tools ‣ Calculate Full Matrix Depending on the scene’s resolution this process can take time.
  2. Save the parametrised scene: File ‣ Save 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.

config/insar_rectangular.gronf (in project folder)
%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
../../_images/example_spool-ridgecrest-result.png

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.