Kinematic Slip Inversion

Kinematic Source Model

The kinematic source model studies also the temporal evolution of coseismic slips. The kinematic source model currently implemented in AlTar is formulated as follows.

The fault plane, as a rectangular area, is divided into \(N_{dd} \times N_{as}\) square patches (dd=down dip, as=along strike), and time is divided into \(N_t\) intervals. The kinematic slip function \({\bf M}_b({\vec \xi}, t)\) (big M) can be in general written as

\[{\bf M}_b({\vec \xi}, t) = {\bf D}({\vec \xi}) S(t-T_R({\vec \xi}); T_r({\vec \xi})),\]

where \({\vec \xi}, t\) label the patch and time, respectively, and

\({\bf D}({\vec \xi})\) is the coseismic slip vector, also subject to the static source model inversion;
\(T_R({\vec \xi})\) is the initial rupture time of each patch, determined by solving the Eikonal equation with a Fast Sweeping algorithm, given the location of the hypocenter \(H_0\) and the rupture velocity \(V_r({\vec \xi})\) (assumed to be isotropic and the same within each patch);
\(T_r({\vec \xi})\) is the slip duration (rise time) for each patch;
\(S(t, T_r)\) is the source time function which is finite only for \(t \in [0, T_r]\) and integrated to 1: we choose a triangular source time function.

To produce smooth synthetics, we refine each patch into \(N_{mesh} \times N_{mesh}\) grids when solving the eikonal equation, while dividing each time interval into \(N_{pt}\) points. \({\bf M}_b({\vec \xi}, t)\) are then interpolated and integrated from these finer spatial and temporal meshes.

The predicted observations can be calculated by,

\[d_{prediction} ({\vec x}, t) = \sum_{{\vec \xi}, t'}{\cal G}_b({\vec x}, {\vec \xi}; t-t') {\bf M}_b({\vec \xi}, t')\]

which is a linear model (note that the Eikonal equation is non-linear). Here, \({\cal G}_b({\vec x}, {\vec \xi}; t-t')\) (big G) are the kinematic Green’s functions relating the source \({\bf M}_b({\vec \xi}, t')\) to an observation location at a given time \(({\vec x}, t)\). The kinematic Green’s functions are pre-calculated as inputs to AlTar. There are some existing software packages for the calculation, e.g., the frequency-wavenumber integration code of Zhu and Rivera, or AXITRA.

In summary, the kinematic inversion uses the following source parameters

the two components of coseismic slip \({\bf D}({\vec \xi})\) (\(2\times N_{dd} \times N_{as}\) elements),

the rupture velocity \(V_r({\vec \xi})\) (\(N_{dd} \times N_{as}\) elements),

the rise time \(T_r({\vec \xi})\) (\(N_{dd} \times N_{as}\) elements),

the location of the hypocenter \({\bf H}_0\) (2 elements),

and the forward model

\[d_{pred} = G( {\bf D}({\vec \xi}), V_r({\vec \xi}), T_r({\vec \xi}), {\bf H}_0)\]

is preformed in two steps,

to obtain \({\bf M}_b({\vec \xi}, t')\) from the Eikonal equation sovler,

to calculate \(d_{pred} = {\cal G}_b {\bf M}_b\).

Joint Kinematic-Static Inversion

Because the slips \({\bf D}({\vec \xi})\) are subject to both static and kinematic observations, we in general run static and kinematic models together, by a model ensemble with or without cascading.

The joint Bayesian probability for static and kinematic models can be written as

\[P({\boldsymbol \theta}_c, {\boldsymbol \theta}_s, {\boldsymbol \theta}_k | {\bf d}_s, {\bf d}_k) = P({\boldsymbol \theta}_c)P({\boldsymbol \theta}_s)P({\boldsymbol \theta}_k) P({\bf d}_s| {\boldsymbol \theta}_c, {\boldsymbol \theta}_s) P({\bf d}_k| {\boldsymbol \theta}_c, {\boldsymbol \theta}_k).\]

where \({\boldsymbol \theta}_c = [strikeslip, dipslip]\) as shared (common) parameters, \({\boldsymbol \theta}_s = [ramp]\), and \({\boldsymbol \theta}_k = [risetime, rupturevelocity, hypocenter]\). The data likelihoods are computed from

the static model with observations \({\bf d}_s\) and the forward model \({\bf d}_s^{pred} = G_s({\boldsymbol \theta}_c, {\boldsymbol \theta}_s)\),
the kinematic model with \({\bf d}_k\) and \({\bf d}_k^{pred} = G_k({\boldsymbol \theta}_c, {\boldsymbol \theta}_k)\).

In the annealing schemes such as CATMIP, we can introduce transitioning distributions

\[P_{\beta_s, \beta_k} ({\boldsymbol \theta}_c, {\boldsymbol \theta}_s, {\boldsymbol \theta}_k | {\bf d}_s, {\bf d}_k) = P({\boldsymbol \theta}_c)P({\boldsymbol \theta}_s)P({\boldsymbol \theta}_k) [P({\bf d}_s| {\boldsymbol \theta}_c, {\boldsymbol \theta}_s)]^{\beta_s} [P({\bf d}_k| {\boldsymbol \theta}_c, {\boldsymbol \theta}_k)]^{\beta_k}.\]

where both \(\beta_s\) and \(\beta_k\) vary from 0 to 1, either jointly or independently. Depending on how \(\beta_s, \beta_k\) vary, AlTar supports two schemes:

In the non-cascading scheme, we set \(\beta=\beta_s=\beta_k\), i.e., both increasing at the same pace, while their increment in the COV scheduler is determined by the coefficient of variation of the weights \(w = [P({\bf d}_s| {\boldsymbol \theta}_c, {\boldsymbol \theta}_s) P({\bf d}_k| {\boldsymbol \theta}_c, {\boldsymbol \theta}_k)]^{\beta_{m+1}-\beta_m}\).
In the cascading scheme, we run AlTar twice:
1. In the first run, we perform inversion on the static model only, or producing samples of \({\boldsymbol \theta}_c\) and \({\boldsymbol \theta}_s\) pursuant to the posterior distribution of the static model \(P({\boldsymbol \theta}_c, {\boldsymbol \theta}_s | {\bf d}_s) = P({\boldsymbol \theta}_c)P({\boldsymbol \theta}_s) [P({\bf d}_s| {\boldsymbol \theta}_c, {\boldsymbol \theta}_s)]\).
2. In the second run, we run static and kinematic models together, by setting \(\beta_s=1\) (cascaded = True) and varying \(\beta_k\) from 0 to 1 (cascaded = False). Here, the samples of \({\boldsymbol \theta}_c\) and \({\boldsymbol \theta}_s\) from the static inversion are used as the initial samples for the joint inversion. The increment of \(\beta_k\) in the COV scheduler is determined by the coefficient of variation of the weights \(w = [P({\bf d}_k| {\boldsymbol \theta}_c, {\boldsymbol \theta}_k)]^{\beta_{k,m+1}-\beta_{k,m}}\).

With the optimized and (usually greatly) reduced search range of \({\boldsymbol \theta}_c\) and \({\boldsymbol \theta}_s\) in parameter space from the static inversion, the cascaded joint-static-kinematic inversion runs more efficiently than the non-cascading scheme. In general, the cascading scheme is recommended for all model ensembles if there is a computation-intensive model (such as the kinematic model) present.

Configurations (Kinematic Model only)

We illustrate the settings of the kinematic model by assuming to run it alone (in practice this is rarely adopted).

An example configuration file

The configuration file (kinematicg_only.pfg) for the kinematic model appears as

; application instance name
slipmodel:

    ; model to be sampled
    model = altar.models.seismic.cuda.kinematicg
    model:

        dataobs:
            observations = 14148 ; number of observed data points
            data_file = kinematicG.data.h5
            cd_std = 5.0e-3
            ; or cd_file = kinematicG.cd.h5 if using a file input

        ; fixed model parameters
        ; green's function (2*Ndd*Nas*Nt, observations)
        ; [Nt][2(strike/dip)][Nas][Ndd] with leading dimensions on the right
        green = kinematicG.gf.h5

        Ndd = 3 ; patches along dip
        Nas = 3 ; patches along strike
        Nmesh = 30 ; mesh points for each patch
        dsp = 20.0 ; length for each patch, km
        Nt = 90 ; number of time intervals
        Npt = 2 ; mesh points for each time interval
        dt = 1.0 ; time unit for each interval, second
        ; initial starting time for each patch, in addition to the fast sweeping calculated arrival time
        t0s = [0.0] * {slipmodel.model.patches}

        ; parameters to be simulated
        ; provide a list at first, serving as their orders in theta
        psets_list = [strikeslip, dipslip, risetime, rupturevelocity, hypocenter]

        ; define each parameterset
        psets:
            strikeslip = altar.cuda.models.parameterset
            dipslip = altar.cuda.models.parameterset
            risetime = altar.cuda.models.parameterset
            rupturevelocity = altar.cuda.models.parameterset
            hypocenter = altar.cuda.models.parameterset

            ; variables for patches are arranged along dip direction at first [Nas][Ndd]
            strikeslip:
                count = {slipmodel.model.patches}
                prep = altar.cuda.distributions.preset ; load preset samples
                prep.input_file = theta_cascaded.h5 ; file name
                prep.dataset = ParameterSets/strikeslip ; dataset name in h5
                prior = altar.cuda.distributions.gaussian
                prior.mean = 0
                prior.sigma = 0.5

            dipslip:
                count = {slipmodel.model.patches}
                prep = altar.cuda.distributions.preset
                prep.input_file = theta_cascaded.h5 ; file name
                prep.dataset = ParameterSets/dipslip ; dataset name in h5
                prior = altar.cuda.distributions.uniform
                prior.support = (-0.5, 20.0)

            risetime:
                count = {slipmodel.model.patches}
                prior = altar.cuda.distributions.uniform
                prior.support = (10.0, 30.0)

            rupturevelocity:
                count = {slipmodel.model.patches}
                prior = altar.cuda.distributions.uniform
                prior.support= (1.0, 6.0)

            ; along strike(first), dip directions
            ; could be separated into 2 for dip and strike direction
            hypocenter:
                count = 2
                prior = altar.cuda.distributions.gaussian
                prior.mean = 20.0
                prior.sigma = 5.0

Parameter Sets

The parameter sets or psets for the kinematic models are psets_list = [strikeslip, dipslip, risetime, rupturevelocity, hypocenter].

The names the parameter sets can be changed per your preference, e.g., strike_slip, StrikeSlip. But the order of the parameter sets must be preserved because the forward model uses the order to map appropriate parameters. strikeslip and dipslip may be switched as long as their order is consistent with the Green’s functions.
strikeslip and dipslip are two components of the cumulative slip displacement. If you prefer to load their initial samples from the static inversion results, use the altar.cuda.distributions.preset distribution for prep, see Preset distribution for more details. Only HDF5 format is accepted for Preset prior and therefore, its dataset name prep.dataset=ParameterSets/strikeslip is also required. If you choose to generate samples from a given distribution, e.g., gaussian/moment scale distributions, please follow the Parameter Sets example in static inversion to set their prep and prior distributions. The slips are usually in unit of meters.
risetime (in unit of seconds) and rupturevelocity (in unit of km/s) are rupture duration and velocities for each patch. As they are positive, usually uniform or truncated gaussian distributions are used as their priors.
strikeslip, dipslip, risetime and rupturevelocity are defined for each patch and their counts are the same as the number of patches. The sequence of patches is arranged as, for \(N_{dd} \times N_{as}\) patches, \((as_0, dd_0), (as_0, dd_1), ... (as_0, dd_{Ndd-1}), (as_1, dd_0), ..., (as_{Nas-1}, dd_{Ndd-1})\). Or dd is the leading dimension.
hypocenter (in unit of km) is the location of the hypocenter measured from the CENTER of the \((as_0, dd_0)\) patch (note that it’s not the origin or the corner), in unit of kilometers. If the distances down dip (dd) and along strike (as) directions are different, you may separate them as two parameter sets hypo_as and hypo_dd, with as component being first.

Input files

The kinematic model requires the following input files

green:: the kinematic Green’s functions, with the shape=(2*Ndd*Nas*Nt, observations). The observations is the number of observed data points, and is the leading dimension. [Nt][2(strike/dip)][Nas][Ndd] labels the spatial-temporal source displacements with leading dimensions on the right (or which comes first):

(t=0, strike, as_0, dd_0, obs_0), (t=0, strike, as_0, dd_0, obs_1), ..., (t=0, strike, as_0, dd_0, obs_{Nobs-1})
(t=0, strike, as_0, dd_1, obs_0), (t=0, strike, as_0, dd_1, obs_1), ..., (t=0, strike, as_0, dd_1, obs_{Nobs-1})
... ...
(t=0, strike, as_0, dd_{Ndd-1}, obs_0), (t=0, strike, as_0, dd_{Ndd-1}, obs_1), ...,  (t=0, strike, as_0, dd_{Ndd-1}, obs_{Nobs-1})
(t=0, strike, as_1, dd_0, obs_0), (t=0, strike, as_1, dd_0, obs_1), ..., (t=0, strike, as_1, dd_0, obs_{Nobs-1})
... ...
(t=0, strike, as_{Nas-1}, dd_{Ndd-1}, obs_0), (t=0, strike, as_{Nas-1}, dd_{Ndd-1}, obs_1), ..., (t=0, strike, as_{Nas-1}, dd_{Ndd-1}, obs_{Nobs-1})
(t=0, dip, as_0, dd_0, obs_0), (t=0, dip, as_0, dd_0, obs_1), ..., (t=0, dip, as_0, dd_0, obs_{Nobs-1})
... ...
(t=0, dip, as_{Nas-1}, dd_{Ndd-1}, obs_0), (t=0, dip, as_{Nas-1}, dd_{Ndd-1}, obs_1), ..., (t=0, dip, as_{Nas-1}, dd_{Ndd-1}, obs_{Nobs-1})
(t=1, strike, as_0, dd_0, obs_0), (t=1, strike, as_0, dd_0, obs_1), ..., (t=1, strike, as_0, dd_0, obs_{Nobs-1})
... ...
(t={Nt-1}, dip, as_{Nas-1}, dd_{Ndd-1}, obs_0), (t={Nt-1}, dip, as_{Nas-1}, dd_{Ndd-1}, obs_1), ..., (t={Nt-1}, dip, as_{Nas-1}, dd_{Ndd-1}, obs_{Nobs-1})

You need to follow the above order when preparing the Green's functions as it's the order how big-M is arranged in the forward model.

dataobs.data_file:: 1d vector of observed data.
dataobs.cd_file:: the data covariance matrix with shape=(observations, observations). If not available, a constant dataobs.cd_std may be used instead.

The input files can be a text file (.txt), a raw binary (.bin or .dat) or an HDF5 (.h5) file, with its format recognized by the file suffix.

Other attributes

Ndd:: integer, number of patches down the dip direction
Nas:: integer, number of patches along the strike direction
Nmesh:: integer, number of mesh points for each patch, i.e., each patch is divided into \(N_{mesh} \times N_{mesh}\) grids for solving the Eikonal equation
dsp:: float, the length for each patch, in km
Nt:: integer, number of time intervals, should be long enough to cover the rupture process
Npt:: integer, number of mesh points for each time interval
dt:: float, time unit for each time interval, in second
t0s:: a list of floats with \(N_{patch}\) elements, initial starting time for each patch, in addition to the fast sweeping calculated arrival time. If configured properly, they can reduce the total number of time intervals needed for the computation.

Configurations (Joint inversion)

The configuration for the joint kinematic-static inversion (kinematicg.pfg) appears as

model = altar.models.seismic.cuda.cascaded
model:
    ; parameters to be simulated (priors)
    ; provide a list at first, serving as their orders in theta
    psets_list = [strikeslip, dipslip, ramp, risetime, rupturevelocity, hypocenter]
    ; define parametersets
    psets:
        ; define the prior for each parameter set
        ; use preset prior to load samples from static inversion for cascading scheme
        ; or use regular priors for non-cascading scheme
        strikeslip = ... ...
        dipslip = ... ...
        ... ...

    ; the model ensemble
    models:
        static = altar.models.seismic.cuda.static
        kinematic = altar.models.seismic.cuda.kinematicg

        static:
            cascaded = True ; or False for non-cascading scheme
            psets_list = [strikeslip, dipslip, ramp]
            ; other static model configurations
            ... ...

        kinematic:
            cascaded = False ; default setting for model
            psets_list = [strikeslip, dipslip, risetime, rupturevelocity, hypocenter]
            ; other kinematic model configurations
            ... ...

Here, the main model is a model ensemble altar.models.seismic.cuda.cascaded, while its embedded-models [static, kinematic] listed as elements of the attribute models (a dict).

The parametersets are properties of the main model and are processed by the main model for sample initializations and prior probability computations. Each embedded-model only requires a psets_list attribute to extract a sub set of parameters from model.psets for its own forward modelling, with the data likelihood computed with respect to its own data observations. The main model collects the data likelihood from all embedded models and assembles them into the Bayesian posterior.

The configuration for each embedded model will be the same as when running it independently, except for an extra flag cascaded (default=``False``) to control the cascading scheme.

For the non-cascading scheme with \(\beta_s = \beta_k = \beta\) varying from 0 to 1 simultaneously, set

static:
    cascaded=False
kinematic:
    cascaded=False

while for the cascading scheme with \(\beta_s =1\), and \(\beta_k = \beta\) varying from 0 to 1, after running the static inversion,

static:
    cascaded=True
kinematic:
    cascaded=False

Examples

The examples for the joint static and kinematic inversion are available at models/seismic/examples. Input files for both static and kinematic models are stored under the 9patch directory.

Cascading Scheme

The first step is to run the static inversion only:

$ slipmodel --config=static.pfg

Please refer to Static Slip Inversion for more details.

The results are saved in the directory results/static specified by the config controller.archiver.output_dir, which include HDF5 files for all or selected annealing steps. The final step (\(\beta=1\)) results are saved in step_final.h5. Copy that file to 9patch directory so that the final samples of strike/dip slips serve as initial samples for the joint inversion:

$ cp results/static/step_final.h5 9patch/theta_cascaded.h5

Please also note that the number of chains job.chains in the static inversion should be the same or larger than that of the joint inversion so that there are enough samples available.

We now can run the joint static-kinematic inversion,

$ slipmodel --config=kinematicg.pfg

The results for the jointly inversion will be saved to results/cascaded, or any other directory by changing controller.archiver.output_dir in kinematicg.pfg.

Non-cascading Scheme

For the non-cascading scheme, you don’t need the step to run static inversion.

You may edit the kinematicg.pfg file (or make a copy at first),

change the static.cascaded to False;

change the prep distributions for strikeslip, dipslip, and ramp from preset to appropriate distributions, e.g., copying them from static.pfg file.

change the output directory controller.archiver.output_dir to, e.g., results/non-cascaded.

Then run the joint inversion:

$ slipmodel --config=kinematicg.pfg

In general, the non-cascading takes long iterations to converge and therefore is slower than the cascading scheme.

Please refer to the AlTar Framework for the Bayesian MCMC framework options and job/output controls. For example, the Adaptive Metropolis Sampler in general has better performance than the fixed-length Metropolis Sampler, which can be selected by setting sampler=altar.cuda.bayesian.adapativemetropolis in the configuration file.

Forward Model Application (new version)

It is essentially the same as the static Forward Model Application. The steps are,

prepare a file with a set of parameters in case input directory;
add the forward problem settings to the configuration file kinematicg.pfg and change the job configuration to run with one GPU,

; the model
model = altar.models.seismic.cuda.cascaded
model:

    ; settings for running forward problem only
    ; forward theta input
    theta_input = kinematicG_mean_model.txt
    ; forward output file
    forward_output = forward_prediction.h5

... ...

job:
tasks = 1 ; number of tasks per host
gpus = 1  ; number of gpus per task
gpuprecision = float32 ; double(float64) or single(float32) precision for gpu computations
;gpuids = [0] ; a list gpu device ids for tasks on each host, default range(job.gpus)

run the plexus command,

$ slipmodel.plexus forward --config=kinematicg.pfg

the predicted data from both static and kinematic models, as well as the bigM, will be saved to one forward_prediction.h5 file.

An example is available at examples/kinematicg.pfg.

Note also that the same script can be used for Bayesian simulation, with the commands,

$ slipmodel --config=kinematicg.pfg
# or
$ slipmodel.plexus sample --config=kinematicg.pfg

See Forward Model Application for more details.

Forward Model Application (old version)

Note

This section describes an old implementation, which will be depreciated in the next release.

When analyzing the results, you may need to run the forward model once for the obtained mean-model or any set of parameters, to produce data predictions in comparison with data observations. Since the kinematic forward model is not straightforward, we provide an additional application for running the forward model only, named kinematicForwardModel.

An example configuration file is available as examples/kinematicg_forward.pfg. You may use the model configuration copied from kinematicg.pfg, with extra settings

; theta input
theta_input = kinematicG_synthetic_theta.txt

; output h5 file name
; data prediction is 1d vector with dimension observations
data_output = kinematicG_synthetic_data.h5
; Mb is 1d vector arranged as [Nt][2(strike/dip)][Nas][Ndd] with leading dimensions on the right
mb_output = kinematicG_synthetic_mb.h5

where theta_input is the input of a mean model or any synthetic model, and data_out and mb_output are output file names for the data predictions and the big M (you can create an animation from it to observe the rupture process).

The forward model application may be run as

$ kinematicForwardModel --config=kinematicg_forward.pfg