Example 4: Static finite-fault estimation¶
In this example we will determine a variable slip distribution for the L’aquila 2009 earthquake by using static InSAR data. The data is the exact same from Example 3, where the overall geometry of the fault plane was estimated. It is a requirement to have Example 3 completed in order to follow the instructions and commands given in this example.
Please make sure that you are one level above the Laquila project folder (created earlier).:
In this example you will first make use of the mode argument. The default is: geometry, which is why it was not necessary to specify it in the earlier examples. Now we will always have to set mode to ffi, which is an abbreviation for finite-fault-optimization. The following command will create a configuration file for the ffi mode called config_ffi.yaml right next to the config_geometry.yaml:
beat init Laquila --mode='ffi' --datatypes=geodetic
It will load the config_geometry.yaml and port arguments that have been specified before to ensure consistency and will only use geodetic data.
The main differences in the two configuration files are in the geodetic_config.gf_config and the problem_config. You may want to have a first glance at the new config. The general structure is the same and the argument names are chosen to give the user an initial idea what these are for. A short explanation for each argument is again given in the module API here, where you can use the search function of your browser to find the argument of interest.
Calculate Greens Functions¶
For the distributed slip optimization a reference fault has to be defined that determines the overall geometry. Once this has been done the problem becomes linear as the only unknown parameters are the slips in rake perpendicular and rake parallel direction. The fault geometry needs to be defined in the geodetic.gf_config.reference_sources.:
gf_config: !beat.GeodeticLinearGFConfig store_superdir: /home/vasyurhm/GF/Laquila reference_model_idx: 0 n_variations: [0, 1] earth_model_name: ak135-f-average.m nworkers: 4 reference_sources: - !beat.sources.RectangularSource lat: 42.29 lon: 13.35 north_shift: 5542.073672733207 east_shift: 10698.176839272524 elevation: 0.0 depth: 2926.702988951863 time: 2009-04-06 01:32:49.19 stf: !pf.HalfSinusoidSTF duration: 0.0 anchor: 0.0 stf_mode: post strike: 144.48588106798306 dip: 54.804317242125464 rake: -114.58259929068664 length: 12219.56636799338 width: 9347.802691161543 velocity: 3500.0 slip: 0.5756726498300268 patch_widths: [2.0] patch_lengths: [2.0] extension_widths: [0.6] extension_lengths: [0.4] sample_rate: 1.1574074074074073e-05
The values shown above are parts of the MAP solution from the optimization from Example 3. The results can been imported through the import command specifiying the –results option. We want to import the results from the Laquila project_directory from an optimization in geometry mode and we want to update the geodetic part of the config_ffi.yaml:
beat import Laquila --results=Laquila --mode='ffi' --datatypes=geodetic --import_from_mode=geometry
Of course, these values could be edited manually to whatever the user deems reasonable.
The reference point(s) on the reference_fault(s) are the top, central point of the fault(s)! Ergo the depth parameter(s) relate(s) to the top edge(s) of the fault(s).
Now we need to specify the dimensions of the patches we want to discretize the reference fault(s) into. The patch_widths and patch_lengths arguments have to be lists and the respective entry needs to result in square patches. The parameters extension_withs and extension_lengths specify by how much the reference fault(s)should be extended in each direction. Example: 0.1 means that the fault is extended by 10% of its width/length value in each direction and 0. means no extension.
Once we decided for the discretization and the reference fault values we can create the discretized fault through.:
beat build_gfs Laquila --datatypes=geodetic --mode=ffi
This will create a directory ($linear_gfs): Laquila/ffi/linear_gfs where the fault geometry is saved as a pickle file ‘fault_geometry.pkl’. We can inspect the geometry of the resulting extended discretized fault wrt. the reference fault with:
beat check Laquila --what=discretization --mode=ffi
This will open an interactive 3d plot of the fault geometry, which looks along the lines of
The grey rectangle shows the geometry of the fault specified under reference_sources and the red rectangle(s) show the extended fault with the respective discretization of the sub-patches. The grey and red dots mark the centres of the reference_fault(s) and the extended faults, respectively. The numbers are the indexes of the repsective sub-patch in the Green’s Function matrix we are going to calculate next.
If the upper edge of the fault would intersect the surface (no topography assumed) due to the extension it is truncated at the intersection and not extended further. Which is why the extent of the red fault is assymetric around the grey reference fault in dip-direction.
To repeat the fault discretization after changing some parameters please add the –force option and the previously fault geometry will be overwritten.:
beat build_gfs Laquila --datatypes=geodetic --mode=ffi --force
The next command starts the calculation of the linear Green’s Function matrixes (also called library) using nworkers CPUs in parallel with unit slip in each slip-direction.:
beat build_gfs Laquila --datatypes=geodetic --mode=ffi --execute
The slip components are not dip-slip and strike-slip, but rake-parallel (uparr in config_ffi.yaml priors) and rake-perpendicular (uperp in config_ffi.yaml priors) wrt to reference_fault(s) rake angle(s). This is following the convention of [Minson2013]. In addition to that there is the component utens, which is normal to the previously mentioned components and would be needed to simulate tensile opening or closing. We ignore that here as we want to model a shear-dislocation.
- This will create two files for each GF library in the $linear_gfs directory:
geodetic_uparr_static_0.traces.npy a numpy array containing the linear GFs
geodetic_uparr_static_0.yaml a yaml file with the meta information
Now we are ready to prepare the optimization setup.
Under the problem_config we find the parameters that we need to adjust:
problem_config: !beat.ProblemConfig mode: ffi mode_config: !beat.FFIConfig regularization: none npatches: 121 initialization: random source_type: RectangularSource stf_type: HalfSinusoid decimation_factors: geodetic: 1 seismic: 1 n_sources: 1 datatypes: [geodetic, seismic] hyperparameters: h_SAR: !beat.heart.Parameter name: h_SAR form: Uniform lower: [-20.0] upper: [20.0] testvalue: [0.0] priors: uparr: !beat.heart.Parameter name: uparr form: Uniform lower: [-0.05] upper: [6.0] testvalue: [1.15] uperp: !beat.heart.Parameter name: uperp form: Uniform lower: [-0.3] upper: [4.0] testvalue: [0.5] utens: !beat.heart.Parameter name: utens form: Uniform lower: [0.0] upper: [0.0] testvalue: [0.0] hierarchicals: Laquila_ascxn_offset: !beat.heart.Parameter name: Laquila_ascxn_offset form: Uniform lower: [-0.004496268249748271] upper: [-0.004496268249748271] testvalue: [-0.004496268249748271] Laquila_ascxn_ramp: !beat.heart.Parameter name: Laquila_ascxn_ramp form: Uniform lower: [-0.00043773457168120667, -0.00023808150002277328] upper: [-0.00043773457168120667, -0.00023808150002277328] testvalue: [-0.00043773457168120667, -0.00023808150002277328] Laquila_dscxn_offset: !beat.heart.Parameter name: Laquila_dscxn_offset form: Uniform lower: [-0.003754963750062188] upper: [-0.003754963750062188] testvalue: [-0.003754963750062188] Laquila_dscxn_ramp: !beat.heart.Parameter name: Laquila_dscxn_ramp form: Uniform lower: [4.978325480108451e-05, -0.00025072248953317104] upper: [4.978325480108451e-05, -0.00025072248953317104] testvalue: [4.978325480108451e-05, -0.00025072248953317104]
The npatches parameter should not be manually adjusted. It is automatically set by running the fault discretizeation step during GF calculation(above).
Please notice the hierarchicals parameters! These are the MAP parameters for the orbital ramps for each radar scene that have been optimized in Example 2. These parameters are imported if the fit_plane parameter in the geodetic_config was set to True. The default is to fix these ramp parameters during the static distributed slip optimization, because leaving them open often results in tradeoffs with patches at greater depth and thus artificial slip is optimized at greater depth. Nevertheless, the user may want to try out to free the upper and lower bounds again to include the parameters into the optimization.
The upper and lower bounds of the two prior variables can be adjusted to reduce the solution space (slip parameters [m]). For the L’aquila earthquake it is highly unlikely to have 6 meters of slip, which is simply the default parameter. A maximum slip of 2 meters in slip parallel direction may be more reasonable. In order to be able to sample the zero value at the lower bound it is necessary to allow for some backslip- ergo negative uparr; here 0.1 might be a reasonable choice.
To also allow for variable rake angles accross the fault we may want to allow some rake perpendicular slip. Here the lower and upper bounds should be set to -1. and 1., respectively.
In order to fix a variable at a certain value, the lower and upper bounds as well as the testvlue need to be set to the same value.
The regularization argument should be set to laplacian to introduce a smoothing constraint that penalizes high slip gradients between neighboring patches. Once this is enabled we need to update the configuration file to initialize the slip-smoothing weight as a random variable in the optimization [Fukuda2008]. Adding the –diff option will display the changes to the config to screen instead of applying them to the file.:
beat update Laquila --mode=ffi --diff --parameters=hypers
Once happy with the displayed changes the changes will be applied to the file with:
beat update Laquila --mode=ffi --parameters=hypers
Sample the solution space¶
Please refer to the ‘Sample the solution space section’ of example 3 example for a more detailed description of the sampling and associated parameters.
Firstly, we only optimize for the noise scaling or hyperparameters (HPs) including the laplacian smoothing weight:
beat sample Laquila --hypers --mode=ffi
Checking the $project_directory/config_ffi.yaml, the hyperparameter bounds show something like:
hyperparameters: h_SAR: !beat.heart.Parameter name: h_SAR form: Uniform lower: [-1.0] upper: [5.0] testvalue: [2.0] h_laplacian: !beat.heart.Parameter name: h_laplacian form: Uniform lower: [-5.0] upper: [5.0] testvalue: [0.5]
Markov Chain initialization¶
The initialization argument determines at which point in the solution space to initialize the Markov Chains. The default value random simply draws a random point in the solution space from the prior distributions for each Markov Chain to be sampled. However, as we are using a laplacian smoothing constraint we can use the non-negative least-squares solution as a starting value for a randomly drawn smoothing weight (from the initial guess on the h_laplacian parameter range) [Fukuda2008]. To do, so we need to set the initialization to “lsq”:
mode_config: !beat.FFIConfig regularization: laplacian npatches: 121 initialization: lsq
The ‘n_jobs’ number should be set to as many CPUs as the user can spare under the sampler_config. The number of sampled MarkovChains and the number of steps for each chain of the SMC sampler should be set to high values as we are optimizing now for ca 250 random variables (if the values from the tutorial haven’t been altered by the user); for example to 5000 and 400, respectively.
With these sampler parameters a huge amount of samples are going to be stored to disk! With the values from the tutorial approximately 140GB of samples are created in the course of the sampling. Please see example 0 for an instruction on how to keep only the important samples to reduce the disk usage. Another way to reduce the required disc space is through the buffer_thinning parameter described here.
Finally, we are set to run the full optimization for the static slip-distribution with:
beat sample Laquila --mode=ffi
Summarize and plotting¶
After the sampling successfully finished, the final stage results have to be summarized with:
beat summarize Laquila --stage_number=-1 --mode=ffi
After that several figures illustrating the results can be created.
For the slip-distribution please run:
beat plot Laquila slip_distribution --mode=ffi
To get histograms for the laplacian smoothing, the noise scalings and the posterior likelihood please run:
beat plot LaquilaJointPonlyUPDATE_wide stage_posteriors --stage_number=-1 --mode=ffi --varnames=h_laplacian,h_SAR,like
For a comparison between data, synthetic displacements and residuals for the two InSAR tracks in a local coordinate system please run:
beat plot Laquila scene_fits --mode=ffi
The plot should show something like this. Here the residuals are displayed with an individual color scale according to their minimum and maximum values.
For a plot using the global geographic coordinate system where the residuals have the same color bar as data and synthetics please run:
beat plot Laquila scene_fits --mode=ffi --plot_projection=latlon
Duputel, Z., Agram, P. S., Simons, M., Minson, S. E., and Beck, J. L. (2014). Accounting for prediction uncertainty when inferring subsurface fault slip. Geophysical Journal International, 197(1):464–482
Fukuda, J. and Johnson, K. M. (2008). A fully Bayesian inversion for spatial distribution of fault slip with objective smoothing. Bulletin of the Seismological Society of America, 98(3):1128–1146
Ragon, T., Sladen, A., Simons, M. Accounting for uncertain fault geometry in earthquake source inversions – I: theory and simplified application, Geophysical Journal International, 214(2):1174–1190