"decompose_insar_velocities" (DIV) is an open-source set of matlab scripts for performing a velocity decomposition (e.g. Watson et al. 2022, and references therein) on multiple overlapping InSAR velocity fields. The code has been written to use InSAR LOS velocities generated by LiCSBAS (https://github.com/yumorishita/LiCSBAS), using interferograms from the COMET-LiCS project (https://comet.nerc.ac.uk/COMET-LiCS-portal/), although any LOS velocities may be used if they are correctly formatted as geotifs.
Written for Matlab 2020a and should work with all newer versions.
For older versions, the main change is replacing tiledlayout with subplot.
This is research code provided to you "as is" with no warranties of correctness.
Andrew Watson, 2022
Watson, A. R., Elliott, J. R., & Walters, R. J. (2022). Interseismic Strain Accumulation Across the Main Recent Fault, SW Iran, From Sentinel‐1 InSAR Observations. Journal of Geophysical Research: Solid Earth, 127(2), e2021JB022674.
The example directory contains velocities for four frames (two ascending and two descending) over the North Anatolian Fault in Turkey.
These have been generated at 1 km resolution using LiCSBAS.
To run the tutorial after cloning/downloading this repository, simply run decompose_insar_velocities('example/north_anatolian_fault.conf') in the Matlab command window with the main repository directory set as your path.
DIV works through the following processing steps, many of which can be altered using within the config file.
- Read the input parameter file that defines how the rest of the script will function.
- Load the line-of-sight velocities, uncertainties, and look vector components for all frames. These are stored in a Matlab cell array, as the dimensions may vary. Optionally, also load a mask, interpolated GNSS velocities, and fault and border polygons for plotting.
- Check and downsample the look vector components if required.
- Optionally perform additional downsampling of all inputs.
- Optionally scale the velocity uncertainties using a semivariogram to mitigate the impact of the reference point (Ou et al. 2022).
- Interpolate inputs onto a common grid, required so that we can perform calculations using data from multiple inputs.
- Optionally merge adajcent frames along-track. If adjacent frames do not overlap, either because they are spatially seperated or because of masking, then the track will be split into two subtracks.
- Optionally correct for the "reference frame effect" (Stephenson et al. 2022), which can produce velocity ramps in the range direction. Requires ITRF2014 plate velocities in No-Net-Rotation. These can be generated using the Unavco Plate Motion Calculator (https://www.unavco.org/software/geodetic-utilities/plate-motion-calculator/plate-motion-calculator.html).
- Shift the relative los-of-sight InSAR velocities into a common reference frame, by tying them to GNSS velocities. This is required to perform the decomposition.
- Optionally generate frame overlap statistics, useful for assessing noise in the InSAR velocities.
- Perform the velocity decomposition to estimate East and Vertical velocities.
- Optionally plot and save outputs.
Ou, Q., Daout, S., Weiss, J. R., Shen, L., Lazecký, M., Wright, T. J., & Parsons, B. E. (2022). Large‐Scale Interseismic Strain Mapping of the NE Tibetan Plateau From Sentinel‐1 Interferometry. Journal of Geophysical Research: Solid Earth, 127(6), e2022JB024176.
Stephenson, O. L., Liu, Y. K., Yunjun, Z., Simons, M., & Rosen, P. (2022). The Impact of Plate Motions on Long-Wavelength InSAR-Derived Velocity Fields. Unpublished.
The example config file provides short descriptions of each option. Below, I give further details.
The number of parallel processes to run for the decomposition, which is performed pixel-by-pixel.
Setting to zero use a for loop as opposed to a parfor loop.
LiCSBAS generates uncertainties using bootstrapping, with a trend towards lower uncertainties close to the reference point due to a reduction in the scatter of the displacement series. This reference point bias can be mitigated by calculating the difference between all points in each uncertainty map and fitting either a spherical or exponential model. Each uncertainty is then scaled by the ratio between the sill and the model value at that distance from the reference point. See Ou et al. (2022) for a full breakdown.
Semivariogram model used to scale the uncertainties.
- exp - Use an exponential model.
- sph - Use a spherical model.
Shift the relative InSAR velocities for each frame/track into the same reference frame as a set of GNSS velocities. Can be performed using either GNSS station velocities (e.g. Hussain et al. 2016) or interpolated GNSS velocities (e.g. Weiss et al. 2020). See Ou et al. (2022) for another example (slightly different method).
Broad steps:
- Project East and North GNSS velocities into line-of-sight.
- Calculate residual between LOS GNSS and LOS InSAR velocities.
- Apply a mask to the residuals using a linear deramp and upper and lower limits, to mitigate the impact of large, short-wavelength signals (e.g. subsidece).
- Smooth the residuals in some manner, either through fitting a polynomial function, or by applying a spatial filter.
- Subtract the smoothed residuals from the InSAR.
In the case of GNSS stations, the residual is calculated using the mean of the InSAR velocities within a given radius of each station (defined below). In the case of interpolated GNSS velocities, the interpolation must be performed before running DIV.
Three options are currently included:
- 0 - disable referencing.
- 1 - use GNSS station velocities.
- 2 - use GNSS interpolated velocities.
Select how the GNSS-InSAR residual is smoothed to create the "referencing surface". Options are:
- 1 - use a polynomial function (based on Weiss et al. 2020).
- 2 - use a median filter, only works with interpolated GNSS velocities (ref2gnss=2) (inspired by Xu et al. 2021).
When using a polynomial to create the referencing surface (ref_type=1), set the order of the 2-D polynomial function. Options are:
- 1 - first order (ax + by + c).
- 2 - second order (ax^2 + b^y2 + cxy + dx + ey + f).
When using a median filter to create the referencing surface (ref_type=2), set the width of the filter window in terms of the number of pixels. Must be an odd number of pixels (so the map area covered by the filter will change if the resolution changes).
When referencing to GNSS station velocities (ref2gnss=1), set the radius around each station that is used to calculate the mean InSAR LOS velocity. Units are the same as the input coordinate system (e.g. in degrees for lon-lat).
Applies downsampling to the input velocities, taking either the mean or the median of a given window size (ds_factor x ds_factor). This can improve the computational requirements of later steps (useful for testing). For ds_method, the options are:
- mean
- median
Using pre-masked velocities can lead to smearing or shrinking of the masked area if additional downsapling is performed within DIV. Hence, I recommend providing unmasked velocities and a matching mask file, which is then optionally applied after both regridding and downsampling.
Merge overlapping frames within each track. This requires that the frame directories have been given in order for each track, although gaps are allowed. Options are:
- 0 - disables.
- 1 - apply the shift to each frame, but leave the frames seperate (i.e. don't take the mean to fully merge them).
- 2 - take the weighted mean of overlapping points, combining multiple frames into a single velocity field. Weights are 1/vstd^2.
Function used to perform the track merging. This is what is fit to the overlapping pixels and then removed from the entire frame to perform the "merge".
- 0 - mean difference.
- 1 - first order polynomial plane (ax + by + c).
- 2 - median difference.
- 3 - mode difference (rounded to the nearest 0.1)
Merge adjacent tracks into a continuous velocity field. This is an experimental method which is useful for comparing velocities with and without GNSS referencing. The LOS velocities for each track are projected into a fixed average LOS, and then the difference is minimised with a static offset. Finally, we take the mean of the overlap. Options are:
- 0 - disables
- 1 - enables
Apply a correction for the "reference frame bias" caused by absolute rigid plate motions in an ITRF reference frame (Stephenson et al. 2022). The "relative" LOS velocities produce by e.g. LiCSBAS are technically in the reference frame of the satellite orbit, which is ITRF 2014. This means that any rigid translation of a plate (i.e. not the deformation that we normally want to measure) will be captured. While this velocity tends to be nearly constant across the area of a frame, the varying LOS makes it appear as a ramp in the range direction. For more details see "https://www.essoar.org/doi/10.1002/essoar.10511538.1". This bias can be mitigated by taking rigid no-net-rotation plate velocities in ITRF (see the UNAVCO plate motion calculator) and projecting them into LOS. Options are:
- 0 - disables
- 1 - enables
Propagate uncertainties on the interpolated GNSS velocities through the decomposition. Options are:
- 0 - disables
- 1 - enables
Set the method for decomposing the line-of-sight velocities into East and Vertical velocities.
- 0 - project the North GNSS velocities into line-of-sight for each frame, and subtract them from the InSAR. Then, decompose into East and Vertical (e.g. Weiss et al., 2020).
- 1 - include the North GNSS velocities as an input to the decomposition, and solve for East, North, and Vertical.
- 2 - decompose into East and a joint Vertical-North plane, and then split the latter into Vertical and North components (e.g. Ou et al., 2022).
- 3 - assume North is zero.
This is a threshold on the value of cond(G), where G is the design matrix for the velocity decomposition. As long as the decomposition is ran for pixels with at least one ascending and descending frame (which is currently forced), then a poorly conditioned G is unlikely to be a problem. Options are:
- 0 - disables
- ~0 - enables with that value
Threshold on the model variance (see Qm), that removes pixels for which either the East or Vertical variances are above the given value. Useful for removing particularly noisy pixels. Options are:
- 0 - disables
- ~0 - enables with that value
Calculate histograms of the differences between overlapping frames, both along and across track. For across-track frames, we project the velocities into a shared LOS. This is a useful test for quantifying the effectiveness of both the along-track merge and the GNSS referencing. Options are:
- 0 - disables
- 1 - enables
For the following parameters, the only options are 0 (disables) or 1 (enables).
Write the decomposed East and vertical velocities to geotiffs, using the same projection as the inputs.
Write the decomposed East and vertical velocities to GMT grd files, using the same projection as the inputs.
Write the regridded and modified frame velocities to geotiffs. Currently, this doesn't work with merged tracks.
Write the along-track overlaps, after subtraction of a given function, to text files. Used for plotting histograms to assess the consistency between adjacent frames.
Plot fault traces on the decomposed velocity maps.
Plot country borders on the decomposed velocity maps.
Plot the input vels as a mosaic of all frames. Useful for checking that the vels have loaded in correctly.
For each frame, plot the original uncertainties, the scaled uncertianties, and the semivariogram.
Plot all scaled uncertainties, split into ascending and descending.
Plot the merged tracks, split into ascending and descending.
For each track, plot the original frame velocities, and the new merged velocity(s). Giving 2 as an input makes this pause on each overlap until the figures are closed.
Plot the residual overlap velocities for all frames.
Plot the merged ascending and descending masks
Plot the InSAR velocities after the plate motion correction has been applied.
Plot individual plate motion corrections, showing the original velocities, the correction, and the corrected velocities.
Plot the referencing surfaces for all frames/tracks, split into ascending and descending.
Plot the original relative velocities, the referencing surface, and the referrenced velocities for each frame/track.
Plot the model uncertainties associated with the decomposed East and Vertical velocities.
Plot the masks for the cond(G) and model variance thresholds.
The following parameters are all paths to inputs files. I recommend using absolute paths, but relative paths will also work.
Path to a .mat file containing the interpolated GNSS velocities.
Path to a .csv file containing GNSS stations velocities.
Path to a text file containing coordinates for faults (see plotting/misc/).
Path to a .mat file containing country border polygons (see plotting/misc/).
Path to a text file containing plate velocities used to apply the plate motion correction.
Path to save output geotifs too.
Prefix for outputs, which then has "vE" etc. appended to it.
These are file ID's used to select which inputs to load from every given framedir.
The file is selected by searching for framedir/*id*.
This allows for multiple versions of the velocities to be toggled between without having to change the framedir paths.
Weiss, J. R., Walters, R. J., Morishita, Y., Wright, T. J., Lazecky, M., Wang, H., ... & Parsons, B. (2020). High‐resolution surface velocities and strain for Anatolia from Sentinel‐1 InSAR and GNSS data. Geophysical Research Letters, 47(17), e2020GL087376.
Hussain, E., Hooper, A., Wright, T. J., Walters, R. J., & Bekaert, D. P. (2016). Interseismic strain accumulation across the central North Anatolian Fault from iteratively unwrapped InSAR measurements. Journal of Geophysical Research: Solid Earth, 121(12), 9000-9019.
Xu, X., Sandwell, D. T., Klein, E., & Bock, Y. (2021). Integrated Sentinel‐1 InSAR and GNSS Time‐Series Along the San Andreas Fault System. Journal of Geophysical Research: Solid Earth, 126(11), e2021JB022579.
GNSS fields file should be a .mat file containing structure with the following:
gnss_field.x - vector of x axis coords (n)
gnss_field.y - vector of y axis coords (m)
gnss_field.N - grid of north gnss vels (mxn)
gnss_field.E - grid of east gnss vels (mxn)
Optional extras:
gnss_field.sN - grid of north 1-sigma uncertainties (mxn)
gnss_field.sE - grid of east 1-sigma uncertainties (mxn)
GNSS stations file should be a .csv file containing the following columns, with one row per stations:
Lon Lat vE vN sE sN cov
where vE and vN are the East and North velocities, sE and sN are the one-sigma uncertainties, and cov is the covariance between vE and vN.
Text file containing fault lines for plotting.
Two column (x y), with each fault segement ended by a ">".
See the GEM faults in plotting/misc/ for an example.
A .mat file containing a structure that defines country outline polygons.
See borderdata.mat in plotting/misc/ for an example.
This contains polygons for most countries.
borders.lon = cell array of vectors containing the longitudes for each polygon
borders.lat = cell array of vectors containing the latitudes for each polygon
borders.places = cell array of strings giving the name of each polygon
Text file of plate motion velocities. Format is [lon lat East North]. No uncertainties are required.
All velocities are expected to be in geotif format.
These can be easily generated from LiCSBAS outputs using the LiCSBAS_flt2geotif.py function.
The spatial extent and pixel spacing of the velocities is read from the geotif metadata.
DIV assumes that positive velocities show motion towards the satellite (default for LiCSBAS).
One sigma uncertainties for each velocity.
0 or 1 value for all velocities. 0 = mask pixel. 1 = keep pixel. Areas outside the velocity field itself can be set to NaN.
Unit vector components (0-1) calculated from the line-of-sight and azimuth of the satellite, for each point.
These are used to project ENU motion into LOS, and vice versa.
Incidence angle and azimuth can be calculated from them (see vel_decomp_vE_vUN).
This work was created while completing a PhD at the University of Leeds, School of Earth and Environment, funded by the Royal Society.
These scripts were designed to work with outputs from LiCSBAS and LiCSAR, projects of COMET, the UK Natural Environment Research Council's Centre for the Observation and Modelling of Earthquakes, Volcanoes and Tectonics.
DIV uses open-access scientific colour maps from https://www.fabiocrameri.ch/colourmaps/ (Crameri et al. 2020), and fault traces from the GEM Active Fault Database (Styron and Pagani, 2020).
I would like to thank Jack McGrath for code advice and debugging, John Elliott for his guidance, and both Dehua Wang and Jessica Payne for their feedback as users.
Andrew Watson, 2022.
Crameri, F., Shephard, G. E., & Heron, P. J. (2020). The misuse of colour in science communication. Nature communications, 11(1), 1-10.
Styron, Richard, and Marco Pagani. “The GEM Global Active Faults Database.” Earthquake Spectra, vol. 36, no. 1_suppl, Oct. 2020, pp. 160–180, doi:10.1177/8755293020944182.