Fault-Inversion-Auto

Zhang et al. updated the inversion method in 2020, making it more efficient and faster; We preserved the overall logic of the code and made minor modifications to suit our individual work requirements and preferences.

Workflow

Load Data Results

[Now] 2026-08-02 13:20:42
[1] Folders are ready 
[2] Earthquake-Info are ready 
[3] Stations are ready 
     Counts: 50 
   Processed station 1/50: CHBH04, samples = 1223
   Processed station 2/50: CHBH10, samples = 367
   ...
   Processed station 49/50: TCGH17, samples = 1251
   Processed station 50/50: TKYH02, samples = 391
[4] Successfully loaded observations (Three-components) 
[5] Successfully clean observations (Three-components) 
[S1] Successfully save stations location 
[S2] Successfully save wave data 
[S3] Successfully save stations Info 

Main Calculation Results

[Loading] Waveforms loaded successfully
   Observed data: 756x3x50
   Origin Sampling: 4.0 Hz
   Target Sampling: 2.0 Hz
[Loading] Stations: 50


[Setting] Fault model parameters:
   Rupture length : 111.17 km
   Rupture width  : 35.24 km
   Subfault size  : [10.0, 10.0] km
   Grid size      : [5, 13]
   Initial source : [3, 7]
   Frequency band : 0.02-0.10 Hz

[Setting] Station screening:
   Initial      : 50
   Removed Dist : 0
   Removed Amp  : 0
   Final        : 50

[Setting] Subfault-Depth correction:
   Source index    : 33
   Shift           : 30.48 km
   Hypocenter      : 42.60 km
   Minimum subfault: 32.90 km
   Maximum subfault: 52.30 km

[Setting] Manually removed 0 stations
    90%-Interval reduces the wavelength to 687/756 

[Do] Preliminary Station Screening
# Source duration (s): 115
# Num.Params: 14101
[Iter 1]  Res=0.811
[Iter 2]  Res=0.784
[Iter 3]  Res=0.772
# Source duration (s): 122.5
# Num.Params: 15076
[Iter 1]  Res=0.771
[Iter 2]  Res=0.700
[Iter 3]  Res=0.683
# Source duration (s): 108
# Num.Params: 13191
[Iter 1]  Res=0.655
[Iter 2]  Res=0.556
[Iter 3]  Res=0.525
# Source duration (s): 123.5
# Num.Params: 15206
[Iter 1]  Res=0.705
[Iter 2]  Res=0.673
[Iter 3]  Res=0.654
[Check] Retained 46 stations after preliminary screening
[Check] Fault-Plane Pre-Processing
   Removed 20 redundant stations
   Remained 26 stations

[Do] Fault-Plane Processing [1/5]
# Source duration (s): 130
# Num.Params: 16051
[Iter 1]  Res=0.558
[Iter 2]  Res=0.481
[Iter 3]  Res=0.448
[Iter 4]  Res=0.438
[Iter 5]  Res=0.434
[Iter 6]  Res=0.432
[Iter 7]  Res=0.431
[Result] Current Moment Magnitude: Mw=7.85
   Zero-slip margins [Top Bottom Left Right] = [0 1 3 2]
   New Grid: [14, 17]; New Source: [9, 8]; New Index: 107

[Do] Fault-Plane Processing [2/5]
# Source duration (s): 130
# Num.Params: 56791
[Iter 1]  Res=0.555
[Iter 2]  Res=0.456
[Iter 3]  Res=0.411
[Iter 4]  Res=0.395
[Iter 5]  Res=0.390
[Iter 6]  Res=0.386
[Iter 7]  Res=0.379
[Iter 8]  Res=0.368
[Result] Current Moment Magnitude: Mw=7.90
   Zero-slip margins [Top Bottom Left Right] = [3 4 3 4]
[Result] The fault boundary is stable---Stop updating

[Result] Final Moment-Magnitude: Mw = 7.89
[Result] Final Grid: [10, 13]
[Result] Final Source: [7, 6]

Plots

---

1. Overview

This project implements an automatic finite-fault inversion workflow for near-field three-component waveform data, with strong-motion records as the primary application. The workflow estimates a kinematic finite-fault rupture model from observed waveforms and precomputed Green’s functions. The final products include:

• final fault geometry and subfault grid;
• slip distribution;
• subfault source-time functions (STFs);
• total STF;
• synthetic waveforms and waveform residuals;
• scalar seismic moment and moment magnitude;
• rupture-related plots and rupture-information output files.

The inversion framework follows the automatic strategy of Zheng et al. (2020), which is based on the Iterative Deconvolution and Stacking (IDS) method. The central idea is to begin from an empirically estimated fault model, invert the near-field waveforms, inspect whether significant slip reaches the fault boundaries, and automatically update the fault dimensions when necessary. This avoids manually prescribing a final fault size before inversion.

The present implementation uses QSSP-based Green’s functions and is designed for earthquakes for which the hypocenter, focal mechanism, and near-field waveform records are available.

---

2. Directory Structure

A recommended project organization is:

FaultInverseAuto/
│
├── Input/
│   ├── Info.txt
│   │
│   ├── Datafunc/
│   │   └── other utility and plotting functions
│   │
│   ├── DataWave/
│   │   └── Original or intermediate waveform data, if used by the
│   │       preprocessing workflow
│   │
│   ├── DataVel/
│   │   ├── DataEarth.mat
│   │   └── Earth-structure, velocity, rigidity, or auxiliary model data
│   │
│   ├── OutSta/
│   │   └── Station coordinate and station-selection files
│   │
│   ├── OutWave/
│   │   └── Processed observed waveform files
│   │
│   └── OutRup/
│       └── Final rupture-model output files, e.g. RupInfo.txt
│
├── Greenlib/
│   └── green_func/
│       └── Precomputed QSSP Green's-function database
│
└── Main.m
└── Other Calculating functions with Fun_*.m

2.1 Input folders and files

Item	Role in the workflow
Info.txt	Earthquake information, including event identification, epicenter, hypocentral depth, magnitude, focal mechanism, and type
OutSta/	Station information used by Fun_LoadWaveStation, including station locations and related metadata.
OutWave/	Processed observed waveforms used by the inversion. The waveform channels are arranged in E/N/U component order.
DataVel/	Auxiliary Earth-model data. This folder is also used in Green’s-function correction and later rupture-information output.
Greenlib/green_func/	QSSP Green’s-function library for subfault-station combinations.
Datafunc/ and func/	MATLAB functions required by the main script. This folder is added to the MATLAB search path using addpath.

2.2 Output products

The calculation produces variables and files describing the final rupture model:

Output	Description
slip	Final slip distribution on the trimmed fault grid, in meters.
substf	Final subfault source-time functions.
syn	Synthetic waveforms calculated from the inversion result.
obr	Observation data returned by the final inversion procedure.
res	Residual history from the IDS inversion.
locasub	Final subfault-center coordinates, stored as latitude and longitude.
locadep	Final subfault-center depths, in kilometers.
grid	Final fault grid, stored as [Ndip, Nstrike].
source	Final hypocentral subfault index, stored as [idip, istrike].
SeisMoment	Final scalar seismic moment, in N·m.
Mw	Final moment magnitude.
RupInfo.txt	Rupture-information file generated by Slip_Info.

---

3. Main Workflow

The automatic inversion is divided into four major stages after the initial earthquake, waveform, and Green’s-function data have been prepared.

Earthquake information + waveform data + Green's functions
                         │
                         ▼
Initial fault-model estimation
                         │
                         ▼
Station selection and preprocessing
                         │
                         ▼
IDS finite-fault inversion
                         │
                         ▼
Slip distribution and boundary inspection
                         │
              ┌──────────┴──────────┐
              │                     │
        Fault is stable       Fault requires update
              │                     │
              ▼                     ▼
         Final trimming      Rebuild geometry and Green's functions
              │                     │
              └───────────── repeat ┘
                         │
                         ▼
Final rupture model, plots, and output files

3.1 Basic loading and initial model setup

The focal mechanism is represented by SDR, generally containing strike, dip, and rake. The variable IndexNode selects one of the two nodal planes:

IndexNode = 1;   % Nodal Plane I
IndexNode = 2;   % Nodal Plane II

Observed waveform data and station information are then loaded; The main waveform array ob is organized by time sample and waveform channel. The three-component channel order follows:

[E1, E2, ..., EN, N1, N2, ..., NN, U1, U2, ..., UN]

The initial fault dimensions, subfault size, grid number, source location within the grid, and frequency band are automatically estimated by: Fun_SetFaultModel according to a magnitude and type

Following Zheng et al. (2020), the initial fault model should be sufficiently large to contain the expected rupture area, while the subfault size should remain compatible with the selected upper cutoff frequency, typical magnitude-dependent settings described are:

Magnitude range	Typical subfault size	Frequency band
$$5.5 \leq M_W < 6.5$$	2 km	0.02–0.5 Hz
$$6.5 \leq M_W < 7.5$$	5 km	0.02–0.2 Hz
$$7.5 \leq M_W < 8.5$$	10 km	0.02–0.1 Hz
$$8.5 \leq M_W < 9.5$$	20 km	0.02–0.05 Hz

---

4. Methods

4.1 Waveform model and IDS inversion

For a fixed fault geometry and focal mechanism, the observed waveform at station/component channel $$j$$ can be represented schematically as:

$$ d_j(t) = \sum_{i=1}^{N_{\rm sub}} G_{ij}(t) * s_i(t) + e_j(t), $$

where:

• $$d_j(t)$$ is the observed displacement waveform;
• $$G_{ij}(t)$$ is the Green’s function between subfault $$i$$ and waveform channel $$j$$;
• $$s_i(t)$$ is the apparent source-time function of subfault $$i$$;
• $$*$$ denotes convolution;
• $$e_j(t)$$ represents observational noise, modelling error, and site-effect residuals;
• $$N_{\rm sub}$$ is the number of subfaults.

The IDS procedure estimates subfault STFs by iterative deconvolution and stacking of waveform information. Because the focal mechanism is fixed during the inversion, the procedure is computationally efficient and appropriate for automatic processing. The resulting source model is rake-fixed; therefore, it primarily resolves the spatial and temporal variation of slip amplitude rather than a fully variable rake-angle field.

The observed waveforms are prepared before inversion by the loading and Green’s-function preprocessing routines. The workflow is designed for near-field displacement waveforms. For strong-motion data, acceleration records are commonly integrated to displacement, with filtering used to suppress low-frequency drift and high-frequency components inconsistent with the point-source approximation of individual subfaults.

4.2 Automatic station selection

Reliable finite-fault inversion requires sufficient station coverage and waveform quality. The project performs two station-selection stages through: Fun_SetInvStations

• The first stage is a preliminary station-screening procedure. When the number of stations is greater than or equal to:

opts.MinStationScreen = 12;

the stations can be divided into candidate groups and evaluated through preliminary inversions. Zheng et al. (2020) used approximately 12 stations as a practical minimum for finite-fault inversion and paired stations with approximately opposite azimuths to improve source coverage.

The partner-station selection can be expressed with a two-dimensional Gaussian weighting based on takeoff angle $$\theta$$ and azimuth $$\phi$$:

$$ f(x,y)=\frac{1}{2\pi\sigma_\theta\sigma_\phi} \exp\left[ -\frac{1}{2} \left( \frac{(x-\theta)^2}{\sigma_\theta^2} + \frac{(y-\phi-180^\circ)^2}{\sigma_\phi^2} \right) \right]. $$

This design helps construct station groups with more balanced azimuthal coverage and reduces the influence of isolated poor-quality records.

• The second stage reduces redundancy in very dense networks:

opts.DenseStationLimit = 40;
opts.MinRemainStation  = 20;

If the station number exceeds the dense-network threshold, azimuthally redundant stations may be removed while retaining a sufficiently large station set. This is consistent with the goal of reducing computation time without strongly degrading spatial coverage.

4.3 Slip estimation from subfault STFs

After IDS inversion, the slip for each subfault is calculated from the time integral of the corresponding STF. In the current implementation:

Conceptually, the slip of subfault $$i$$ is:

$$ D_i = \frac{\sum_t s_i(t)} {\mu_{\rm ref} A}, $$

where:

• $$D_i$$ is the estimated slip;
• $$s_i(t)$$ is the subfault STF;
• $$\mu_{\rm ref} = 3\times10^{10}$$ Pa is the reference rigidity used in the STF-to-slip conversion;
• $$A$$ is the subfault area in square meters.

Because gridsize is given in kilometers, the code includes the factor:

$$ 3\times10^{16}= 3\times10^{10} \times 10^6, $$

which converts the subfault area from km² to m².

4.4 Automatic fault-plane update

A fault plane that is too small artificially confines seismic moment near its boundaries and can produce unrealistically large edge slip. Conversely, a fault plane that is unnecessarily large increases the number of unknown subfaults and the computational cost.

After each inversion, the slip distribution is normalized:

$$ \widetilde{D}_i = \frac{D_i}{D_{\max}}, $$

where:

$$ D_{\max} = \max_i(D_i). $$

Only the main rupture region is retained for fault-boundary inspection:

$$ \widetilde{D}_i = 0, \qquad \text{if } \widetilde{D}_i \leq 0.2. $$

The algorithm then checks whether significant slip reaches the top, bottom, left, or right boundaries of the fault. If significant edge slip remains, the current fault plane is not regarded as sufficiently stable. The grid is revised using:

[newGrid, newSource] = subfaults_making(NormSlip, grid, source);

The updated geometry is shifted so that the hypocentral subfault remains at the prescribed hypocentral depth. Rows above the free surface are removed automatically. New Green’s functions are then assembled for the updated fault model:

[g, locasub, locadep, muSubfault] = Fun_RebuildGreenModel(...);

The outer-loop update process is controlled by:

opts.MaxOuterIteration = 5;
opts.InnerIteration    = 20;

Thus, up to five fault-geometry updates may be considered, while each IDS calculation uses the specified number of internal iterations.

4.5 Final trimming and seismic moment

After the final inversion, subfaults outside the main rupture region are removed through: Fun_SetInvFaultTrim

The final scalar seismic moment is calculated using the local subfault rigidity:

$$ M_0 = \sum_{i=1}^{N_{\rm sub}} \mu_i D_i A, $$

Under the usual SI convention, the moment-magnitude relation is:

$$ M_W = \frac{2}{3} \left( \log_{10} M_0 - 9.1 \right), $$

when $$M_0$$ is expressed in N·m.

---

5. Special Features of This Implementation

1. Automatic initial fault model.
   Fun_SetFaultModel estimates rupture dimensions, subfault size, grid number, frequency band, and initial source position from earthquake magnitude and fault type. This reduces dependence on manual trial-and-error fault construction.

2. Depth-consistent fault geometry.
   Fun_SetSourceDep adjusts the designed fault plane so that the designated hypocentral subfault matches the reported earthquake depth. During subsequent grid updates, the new geometry is again shifted to preserve this depth constraint.

3. Automatic station management.
   The workflow combines preliminary waveform-based screening with dense-network redundancy reduction. This is important because poor-quality records, strong local site effects, or highly clustered station distributions can degrade automatic inversion results.

4. Iterative fault-boundary control.
   Instead of assuming that the empirical initial fault size is final, the workflow uses the inferred slip distribution to determine whether the fault should be expanded, reduced, or accepted.

5. Consistent final trimming.
   Green’s functions, subfault locations, depths, rigidity values, STFs, and slip values are trimmed together, ensuring that the final fault model remains internally consistent.

6. Three-dimensional visualization.
   The final slip model can be plotted both on the actual subsurface fault plane and on its surface projection using: Fun_PlotFaultSlip

---

Notes and Recommendations

1. Focal mechanism and fault-plane ambiguity
   The IDS inversion uses a fixed focal mechanism. Therefore, the selected nodal plane strongly affects the resulting fault orientation and rupture-directivity interpretation. When the causative plane is uncertain, both nodal planes should be tested:

2. Station coverage is critical
   The inversion is most reliable when stations surround the source with limited azimuthal gaps and when the nearest source-to-station distance is not excessively large compared with rupture dimensions. Sparse or one-sided station distributions may weaken the resolution of slip location, rupture direction, and rupture speed.

3. Green’s-function quality controls model quality
   The method depends on the QSSP Green’s-function database and the assumed velocity structure. Inaccurate structural models, inappropriate source depths, or inconsistent waveform preprocessing can increase waveform residuals and bias the inferred rupture model.

---

Reference

Zheng, X., Zhang, Y., Wang, R., Zhao, L., Li, W., & Huang, Q. (2020). Automatic inversions of strong-motion records for finite-fault models of significant earthquakes in and around Japan. Journal of Geophysical Research: Solid Earth, 125, e2020JB019992.

Strong-Motion Data Sources

Database / Organization	Data Code	Website
Center for Engineering Strong Motion Data	CESMD	strongmotioncenter.org
Pacific Earthquake Engineering Research Center Strong Ground Motion Databases	PEER	peer.berkeley.edu
European Strong-Motion Database	ESD	isesd.hi.is
USGS National Strong-Motion Project	NSMP	earthquake.usgs.gov
Canadian National Seismograph Network	CNSN	earthquakescanada.nrcan.gc.ca
Italian Accelerometric Archive	ITACA	itaca.mi.ingv.it
Iran Strong Motion Network	ISMN	ismn.bhrc.ac.ir
Institute of Engineering Seismology and Earthquake Engineering Research, Research and Technical Institute	ITSAK	itsak.gr
Strong-Motion Database of Turkey	TR-NSMN	kyhdata.deprem.gov.tr
GeoNet Strong-Motion Data	—	geonet.org.nz

Note: Availability, download procedures, registration requirements, and data formats may vary by database. Check each official website for the latest access instructions and metadata.