Fault-Inversion-Using-Green-Function

Previously, Teacher Zhang and others developed a method for inverting seismic parameters and fault slip. Here, we present additional examples and explanations, along with corrections to some potentially non-standard aspects. We welcome discussion.

• First, you need to calculate the Green's function and prepare far-field or strong-motion waveform files for the specific field; the final required formats are ob_tele.mat and ob_sm.mat
  
• We have consolidated all the data to be imported into the Data file and provided two functions for data processing, with specific parameters described in the function documentation
  
• Green_tele and Green_sm are tools for calculating Green's functions (already package .exe). Please modify the relevant paths and adjust the parameters accordingly (grn.inp). The computation may take several hours, so we have not provided specific values in the file
  
• An example is the 2016 Mw7.8 New Zealand event

Refer to:

• Yong Zhang, Wanpeng Feng, Yuntai Chen, Lisheng Xu, Zhenhong Li, David Forrest, The 2009 L'Aquila MW 6.3 earthquake: a new technique to locate the hypocentre in the joint inversion of earthquake rupture process, Geophysical Journal International, Volume 191, Issue 3, December 2012, Pages 1417–1426, https://doi.org/10.1111/j.1365-246X.2012.05694.x
  
• Initial respiratory： https://github.com/ZhangYongGroupPKU/Finite-fault-joint-inversion.git
  

Test Results:

# [1] Parameters are ready

# Notice: Teleseismic observations and Green's functions have been integrated once (velocity -> displacement).
# Teleseismic Green's functions have been recalculated and saved:
  Data\obg_tele.mat
# Notice: Strong-motion observations and Green's functions have been integrated once (velocity -> displacement).
# Strong-motion Green's functions have been recalculated and saved:
  Data\obg_sm.mat
# [2] Data-Input are ready

# [Now] 2026-07-27 22:46:21 
#   Iteration: 80 / 80
#   Residual: 42.68 / 41.80
#   Stop: Maximum iteration reached
# [Process] 00:00:21 time 
#   Iteration: 75 / 80
#   Residual: 42.63 / 41.80
#   Stop: Maximum iteration reached
# [Re-Process] 00:00:24 time 
# [Final] 2026-07-27 22:47:07 

# Moment magnitude= 8.09
# [3] Finished

Finite-Fault Rupture Inversion Workflow

This MATLAB workflow performs a joint finite-fault rupture inversion using one or more of the following datasets:

1. Teleseismic P-wave records
2. Strong-motion waveforms
3. GPS static displacement
4. InSAR displacement

The inversion estimates the spatial and temporal evolution of fault slip, including slip distribution, rupture timing, source time function, and waveform/static-data fits.

1. Main Script Overview

The main script is organized into four major parts:

• Input parameter definition
• Data and Green's-function preparation
• Joint inversion
• Post-processing and plotting

The main output includes:

• Total slip on each subfault;
• Slip-rate/source-time function for each subfault;
• Overall source time function;
• Rupture propagation;
• Synthetic waveform fitting for teleseismic and strong-motion data;
• Optional GPS and InSAR fitting results;
• Estimated seismic moment and moment magnitude.

---

2. Important Input Parameters

2.1 Fault Geometry

EveLoc = [-42.725,173.065];  % Event location: [latitude, longitude]
depth  = 19;                 % Hypocentral depth (km)
Focal  = [219,38,128];       % [Strike, Dip, Rake] in degrees
NGrid  = [8,25];             % Number of subfaults: [dip, strike]
LGrid  = [10,10];            % Subfault size: [dip_length, strike_length] in km
SGrid  = [3,21];             % Hypocenter grid index: [dip, strike]

• NGrid defines the discretized fault plane.
• LGrid defines the size of each subfault.
• SGrid identifies the reference subfault containing the hypocenter.
• Focal = [strike, dip, rake] defines the fault-plane mechanism.

Important: SGrid, NGrid, and LGrid jointly control the relative position of the hypocenter on the fault plane.

---

2.2 Teleseismic Data Parameters

UseTe     = 1;               % Enable teleseismic data
WayTe     = 0;               % GF method: approximate = 0, exact = 1
DoSTF     = 0;               % Enable apparent STF inversion
BandTe    = [0.01,0.2];      % Teleseismic band-pass frequency range (Hz)
DataTe    = 'Data\ob_tele.mat';
DatagTe   = 'Data\obg_tele.mat';
PathTe    = '...\Green_tele\green_func\';
DoGreenTe = 0;               % Recalculate Green's functions if set to 1
ShiftTe   = 1;               % Allow waveform shifting after initial inversion

Teleseismic Input File: ob_tele.mat

The raw teleseismic data file should contain at least:

ob
loca
da
mm
srate
leng
leng_before

Variable	Description
ob	Teleseismic observed waveforms, usually vertical/P-wave component
loca	Station locations, typically [latitude, longitude]
da	Event-station geometry, e.g., distance and azimuth
mm	Station metadata
srate	Sampling rate in Hz
leng	Length of waveform time window in seconds
leng_before	Time before the P-wave arrival in seconds

---

2.3 Strong-Motion Data Parameters

UseSm     = 1;
BandSm    = [0.02,0.2];
DataSm    = 'Data\ob_sm.mat';
DatagSm   = 'Data\obg_sm.mat';
PathSm    = '...\Green_sm\green_func\';
DoGreenSm = 0;
ShiftSm   = 1;

Strong-motion observations are generally stored as:

ob(Nsample, 3, Nstation)

where the second dimension corresponds to three components, normally East, North, and Up.

---

2.4 Inversion Parameters

Tspan     = 110;     % Source duration for each subfault (s)
Tmax      = 30;      % Maximum rupture time (s)
Vmax      = 3;       % Maximum rupture velocity (km/s)
lambda1   = 50;      % Spatial smoothing weight
lambda2   = 100;     % Temporal smoothing weight
lambda3   = 0.01;    % Minimum scalar moment constraint
Rakevary  = 1;       % Allow rake variation within approximately +/-45 degrees
IsEdge    = 1;       % Edge-slip constraint

Variable	Description
Tspan	Duration of the local source time function on each subfault
Tmax	Maximum allowed rupture-initiation time
Vmax	Maximum rupture velocity constraint
lambda1	Controls spatial smoothness of slip
lambda2	Controls temporal smoothness of slip rate
lambda3	Stabilizes the inversion through a minimum-moment constraint
Rakevary	Allows rake variation; set to 0 during STF inversion
IsEdge	Controls zero-slip conditions at fault edges

---

3. Teleseismic Green's-Function Methods

Three teleseismic processing functions are available.

Fun_getobg_tele3

[obt,gt,locat,mmt,s1,leng] = Fun_getobg_tele3(...);

• Approximate Green's-function method;
• Recommended for standard inversion without STF inversion;
• Efficient because reference Green's functions are shifted according to relative travel times.

Fun_getobg_tele2

[obt,gt,locat,mmt,s1,leng] = Fun_getobg_tele2(...);

• Exact Green's-function method;
• Calculates Green's functions for each subfault-station geometry;
• More computationally expensive than the approximate method.

Fun_getobg_tele4

[obt,gt,locat,mmt,s1,leng,t,t0] = Fun_getobg_tele4(...);

• Approximate Green's-function method for STF-related processing;
• Returns:
  • t: P-wave travel time from each subfault to each station;
  • t0: P-wave travel time from the hypocenter to each station.
• The relative travel-time matrix is calculated as:

tdt = t - repmat(t0',size(t,1),1);

---

4. Green's-Function Cache File

The shared teleseismic cache file is:

DatagTe = 'Data\obg_tele.mat';

Most teleseismic functions save the following variables:

obt
gt
locat
mmt
s1
leng

When Fun_getobg_tele4 is used, it additionally saves:

t
t0

Important cache behavior:
All teleseismic methods use the same DatagTe file. Therefore, the most recently calculated method may overwrite the previous cached Green's functions.

If Fun_getobg_tele4 is called with:

DoGreenTe = 0;

but the cache does not contain t and t0, set:

DoGreenTe = 1;

and run Fun_getobg_tele4 once to regenerate the cache.

Similarly, Green's functions should be recalculated after changing important source parameters, including:

EveLoc, depth, Focal, NGrid, LGrid, SGrid, miun, PathTe, IsDisp

---

5. Processing Workflow

Step 1 — Calculate Shear Modulus

load('earth.mat');
[miu1,miun] = getmiu(depth,Focal,NGrid,LGrid,SGrid,EveLoc,earth);

miun is the shear modulus for each subfault and is used to scale the Green's functions and calculate seismic moment.

---

Step 2 — Prepare Observations and Green's Functions

Teleseismic and strong-motion records are loaded, optionally integrated from velocity to displacement, filtered, padded, and combined.

[obt,gt] = length_unification(obt,gt,obsm);
srate = set_srate(s1,s2);

ob = [obt,obsm];
g  = cat(4,gt,gsm);

The observations and Green's functions are then normalized:

[ob,g,mob] = normalize_obg(ob,g,obt,obsm,UseTe,UseSm);

---

Step 3 — Construct the Inversion Matrix

[G,sGDT0,hdata,meangreen,sumob,trup_con,locasub,locadep,trup] = ...
    get_Gmatrix(...);

This step creates:

• The waveform Green's-function matrix;
• Rupture-time constraints;
• Spatial and temporal smoothing constraints;
• Fault/subfault geometry;
• Data normalization factors.

GPS and InSAR constraints, if enabled, are appended to the inversion system.

---

Step 4 — Perform Joint Inversion

[Xsolve,rsq] = Fun_MainInverse(G_jt,sGDT,ObsUsing,1e-170,1200);

The inversion solution Xsolve contains the fault-slip model parameters.

After the first inversion, optional waveform time shifts are applied:

[Xsolve,ob0,syn,astftele2,synastf2] = shift_ob(...);

---

Step 5 — Extract Results and Plot Figures

Main rupture results are obtained using:

[slip,stf,slipx,slipy,slipe,slipn,zslip,stfsub_out] = ...
    get_result2(...);

The workflow then produces figures for:

• Subfault source time functions;
• Teleseismic station distribution;
• Strong-motion station distribution;
• Static slip distribution;
• Global source time function;
• Teleseismic waveform fitting;
• Strong-motion waveform fitting;
• GPS/InSAR fitting, if enabled;
• Cumulative rupture/slip snapshots;
• STF inversion fitting, if enabled.

---

6. Improvements in the Updated Teleseismic Functions

Compared with the original implementations, the updated Fun_getobg_tele2 and Fun_getobg_tele4 include the following improvements:

• Explicit MAT-file variable loading:
  Raw data are loaded into a structure first, then assigned explicitly. This avoids hidden workspace variables caused by load(obfile).

• Unified naming convention:
  Function outputs use:

  obt, gt, locat, mmt, s1, leng

  which matches the main inversion script.

• Improved cache management:
  Fun_getobg_tele4 stores t and t0 in DatagTe, ensuring that STF-related travel-time information can be reused.

• Clearer error messages:
  Green-function extraction indices, sampling dimensions, frequency bands, and shear-modulus dimensions are checked before processing.

• Explicit time-dimension processing:
  Integration and filtering are applied along the first dimension, which is the waveform time axis.

• Improved comments and readability:
  The Green-function generation, subfault geometry construction, travel-time calculations, and scaling procedures are documented in English.

---