Module 4

Earthquake Monitoring with InSAR

Measuring co-seismic ground deformation using Synthetic Aperture Radar interferometry, from elastic dislocation theory to operational processing pipelines

Okada Elastic Dislocation Model

The Okada (1985) model provides analytical expressions for the surface displacement caused by a rectangular fault in an elastic half-space. It is the fundamental forward model used to relate earthquake source parameters to observable surface deformation.

General Representation Theorem

The displacement at a point $\mathbf{x}$ on the surface due to slip on a fault surface $\Sigma$ is:

$$\mathbf{u}(\mathbf{x}) = \int_{\Sigma} \Delta u_i \, c_{ijkl} \, n_l \, G_{jk} \, d\Sigma$$

where $\Delta u_i$ is the slip vector,$c_{ijkl}$ is the elastic stiffness tensor,$n_l$ is the fault normal, and$G_{jk}$ is the Green's function for the elastic half-space.

Okada Fault Parameters

Strike (phi):Azimuth of fault trace, clockwise from North

Dip (delta):Angle of fault plane from horizontal (0-90 deg)

Rake (lambda):Direction of slip on fault plane (-180 to 180 deg)

Length (L):Along-strike dimension of the fault

Width (W):Down-dip dimension of the fault

Depth (d):Depth to top edge of the fault

Slip (U):Magnitude of displacement on the fault

Location (x0, y0):Surface projection of fault center

Seismic Moment & Magnitude

The seismic moment relates fault geometry to earthquake magnitude:

$$M_0 = \mu \cdot L \cdot W \cdot \bar{U}$$

where $\mu \approx 30$ GPa is the shear modulus and $\bar{U}$ is the average slip. The moment magnitude is:

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

LOS Displacement Projection

InSAR measures only the line-of-sight (LOS) component of the 3D displacement vector. The projection from East, North, Up components to LOS is:

$$d_{LOS} = u_E \sin\theta \cos\alpha - u_N \sin\theta \sin\alpha + u_U \cos\theta$$

where $\theta$ is the incidence angle (typically 20-45 deg for Sentinel-1) and $\alpha$ is the satellite heading angle (azimuth of the satellite ground track, ~-12 deg for ascending, ~-168 deg for descending passes).

Key insight: InSAR is most sensitive to vertical displacement (cos theta factor ~0.7-0.9) and east-west motion, but nearly blind to north-south displacement. Combining ascending and descending passes can decompose the 2D (E, U) displacement field.

DInSAR Processing Pipeline

The standard DInSAR processing chain for earthquake deformation mapping using SNAP (Sentinel Application Platform) or similar tools:

Co-registration

Align the secondary (post-event) SLC to the primary (pre-event) SLC with sub-pixel accuracy using orbital data and a DEM.

Interferogram Formation

Multiply the primary by the complex conjugate of the secondary to form the interferogram. The phase encodes path length differences.

Topographic Phase Removal

Simulate and subtract the topographic phase contribution using a reference DEM (e.g., Copernicus 30m DEM).

Phase Filtering

Apply Goldstein adaptive filter to reduce phase noise while preserving fringe continuity. Window size and alpha parameter control smoothing.

Phase Unwrapping

Convert wrapped phase (modulo 2pi) to continuous phase using SNAPHU or MCF algorithms. Most error-prone step.

Phase to Displacement

Convert unwrapped phase to line-of-sight displacement using the radar wavelength.

Geocoding

Transform from radar geometry (range/azimuth) to geographic coordinates using DEM and orbital information.

Phase to Displacement Conversion

After unwrapping, the phase is converted to LOS displacement:

$$\Delta r = -\frac{\lambda}{4\pi}\phi_{unwrapped}$$

For Sentinel-1 C-band ($\lambda = 5.547$ cm),$\lambda / 4\pi = 0.442$ cm/radian. One complete fringe cycle ($2\pi$) corresponds to$\lambda/2 = 2.774$ cm of LOS displacement. The negative sign convention means positive phase = motion toward the satellite.

PSInSAR & Time Series Analysis

While DInSAR uses a single interferometric pair, Persistent Scatterer InSAR (PSInSAR) exploits long time series of SAR acquisitions to identify stable radar targets and estimate their deformation history.

PS Velocity Estimation

The PS velocity is estimated by maximizing the temporal coherence over the stack of interferograms:

$$v_{PS} = \arg\min_v \sum_{k=1}^{N} \left| e^{i(\phi_k - \frac{4\pi}{\lambda} v t_k)} \right|^2$$

where $\phi_k$ is the differential phase of the $k$-th interferogram,$v$ is the linear velocity, and$t_k$ is the temporal baseline. The temporal coherence $\gamma_t$ measures how well the linear velocity model fits the observed phases.

PSInSAR vs SBAS

PSInSAR (Ferretti et al., 2001)

• Single primary image
• Point-like stable targets (buildings, rocks)
• Best for urban areas and infrastructure
• mm/yr velocity precision
• Tools: StaMPS, SARPROZ, ISCE

SBAS (Berardino et al., 2002)

• Multiple primary images (short baselines)
• Distributed scatterers (fields, desert)
• Better spatial coverage in rural areas
• Slightly lower precision than PSInSAR
• Tools: MintPy, LiCSBAS, ISCE

Key applications: PSInSAR time series are used to monitor urban subsidence (e.g., Mexico City sinking at 30 cm/yr), volcanic inflation/deflation, post-seismic relaxation, landslide creep, and infrastructure stability (bridges, dams, tunnels).

Okada Fault Model: Surface Displacement Simulation

The following code implements a simplified Okada rectangular fault model for a strike-slip earthquake. It computes the 3D surface displacement field and projects it to SAR line-of-sight to simulate interferometric fringes.

Simulation

Python

script.py228 lines

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import Normalize

# ============================================================
# Simplified Okada (1985) dislocation model
# Point-source approximation for a rectangular fault
# ============================================================

def okada_point(x, y, depth, strike_deg, dip_deg, rake_deg, slip, nu=0.25):
    """
    Simplified Okada point-source surface displacement.
    Returns (ux, uy, uz) in East, North, Up.

Parameters:
        x, y: observation coordinates relative to fault center (m)
        depth: depth to fault center (m)
        strike_deg: fault strike (degrees CW from North)
        dip_deg: fault dip (degrees from horizontal)
        rake_deg: slip direction on fault plane (degrees)
        slip: total slip (m)
        nu: Poisson's ratio (default 0.25)
    """
    strike = np.radians(strike_deg)
    dip = np.radians(dip_deg)
    rake = np.radians(rake_deg)

# Rotate coordinates to fault-aligned system
    xs = np.cos(strike) * x + np.sin(strike) * y
    ys = -np.sin(strike) * x + np.cos(strike) * y

# Strike-slip and dip-slip components
    Us = slip * np.cos(rake)  # strike-slip
    Ud = slip * np.sin(rake)  # dip-slip

d = depth
    R = np.sqrt(xs**2 + ys**2 + d**2)
    R3 = R**3
    R5 = R**5

# Mogi-like + dislocation terms (simplified)
    # Using point-source Volterra dislocation in half-space
    factor = 1.0 / (4.0 * np.pi * (1.0 - nu))

# Strike-slip contribution
    ux_s = -factor * Us * d * (
        (2*nu - 1) * ys / R3 +
        3 * d * ys * (xs**2 + d**2) / R5
    ) * np.cos(dip)

uy_s = factor * Us * d * (
        (2*nu - 1) * xs / R3 +
        3 * d * xs * (ys**2 + d**2) / R5
    ) * np.cos(dip)

uz_s = -factor * Us * d * (
        3 * d**2 * ys / R5
    ) * np.sin(dip)

# Dip-slip contribution
    ux_d = factor * Ud * (
        (1 - 2*nu) * xs * np.sin(dip) / R3 +
        3 * d * xs * ys * np.cos(dip) / R5
    )

uy_d = factor * Ud * (
        (1 - 2*nu) * ys * np.sin(dip) / R3 -
        3 * d * (xs**2 + d**2) * np.cos(dip) / R5
    )

uz_d = factor * Ud * (
        (1 - 2*nu) * d * np.cos(dip) / R3 +
        3 * d**2 * ys * np.sin(dip) / R5
    )

# Rotate back to geographic coordinates
    ue = ux_s + ux_d
    un = uy_s + uy_d
    uu = uz_s + uz_d

# Rotate from fault-aligned back to ENU
    ue_geo = np.cos(strike) * ue - np.sin(strike) * un
    un_geo = np.sin(strike) * ue + np.cos(strike) * un

return ue_geo, un_geo, uu

# ============================================================
# Fault parameters (example: M6.5 strike-slip earthquake)
# ============================================================
fault_params = {
    'depth': 8000,       # 8 km depth
    'strike_deg': 45,    # NE-SW strike
    'dip_deg': 80,       # near-vertical
    'rake_deg': 0,       # pure left-lateral strike-slip
    'slip': 1.5,         # 1.5 m of slip
}

# Observation grid (40 km x 40 km)
n = 300
x1d = np.linspace(-20000, 20000, n)
y1d = np.linspace(-20000, 20000, n)
X, Y = np.meshgrid(x1d, y1d)

# Compute displacement field
# Sum contributions from multiple point sources along fault
fault_length = 15000  # 15 km
fault_width = 10000   # 10 km
n_sources = 20

ue_total = np.zeros_like(X)
un_total = np.zeros_like(X)
uu_total = np.zeros_like(X)

for i in range(n_sources):
    frac = (i + 0.5) / n_sources - 0.5
    dx = frac * fault_length * np.cos(np.radians(fault_params['strike_deg']))
    dy = frac * fault_length * np.sin(np.radians(fault_params['strike_deg']))

ue, un, uu = okada_point(
        X - dx, Y - dy,
        fault_params['depth'],
        fault_params['strike_deg'],
        fault_params['dip_deg'],
        fault_params['rake_deg'],
        fault_params['slip'] * fault_length * fault_width / n_sources * 1e-9
    )
    ue_total += ue
    un_total += un
    uu_total += uu

# Scale to get realistic displacement magnitudes
scale = 5e4
ue_total *= scale
un_total *= scale
uu_total *= scale

# ============================================================
# Project to SAR line-of-sight
# ============================================================
theta_inc = np.radians(34)   # Sentinel-1 incidence angle
alpha_head = np.radians(-12) # ascending heading

d_los = (ue_total * np.sin(theta_inc) * np.cos(alpha_head) -
         un_total * np.sin(theta_inc) * np.sin(alpha_head) +
         uu_total * np.cos(theta_inc))

# Wrap to simulate interferometric fringes
wavelength = 0.05547  # Sentinel-1 C-band wavelength (m)
phase_wrapped = np.angle(np.exp(1j * 4 * np.pi * d_los / wavelength))

# ============================================================
# Plot results
# ============================================================
fig, axes = plt.subplots(2, 3, figsize=(15, 10))
fig.suptitle('Okada Fault Model: M~6.5 Strike-Slip Earthquake', fontsize=16, fontweight='bold', color='white')

extent = [-20, 20, -20, 20]

# East displacement
im0 = axes[0, 0].imshow(ue_total * 100, extent=extent, cmap='RdBu_r', origin='lower',
                          vmin=-15, vmax=15)
axes[0, 0].set_title('East Displacement (cm)', fontsize=11, color='white')
plt.colorbar(im0, ax=axes[0, 0], shrink=0.8)

# North displacement
im1 = axes[0, 1].imshow(un_total * 100, extent=extent, cmap='RdBu_r', origin='lower',
                          vmin=-15, vmax=15)
axes[0, 1].set_title('North Displacement (cm)', fontsize=11, color='white')
plt.colorbar(im1, ax=axes[0, 1], shrink=0.8)

# Vertical displacement
im2 = axes[0, 2].imshow(uu_total * 100, extent=extent, cmap='RdBu_r', origin='lower',
                          vmin=-5, vmax=5)
axes[0, 2].set_title('Vertical Displacement (cm)', fontsize=11, color='white')
plt.colorbar(im2, ax=axes[0, 2], shrink=0.8)

# LOS displacement
im3 = axes[1, 0].imshow(d_los * 100, extent=extent, cmap='RdBu_r', origin='lower',
                          vmin=-15, vmax=15)
axes[1, 0].set_title('LOS Displacement (cm)', fontsize=11, color='white')
plt.colorbar(im3, ax=axes[1, 0], shrink=0.8)

# Wrapped interferogram (fringes)
im4 = axes[1, 1].imshow(phase_wrapped, extent=extent, cmap='hsv', origin='lower',
                          vmin=-np.pi, vmax=np.pi)
axes[1, 1].set_title('Wrapped Interferogram (fringes)', fontsize=11, color='white')
plt.colorbar(im4, ax=axes[1, 1], shrink=0.8, label='Phase (rad)')

# Displacement magnitude
disp_mag = np.sqrt(ue_total**2 + un_total**2 + uu_total**2)
im5 = axes[1, 2].imshow(disp_mag * 100, extent=extent, cmap='hot_r', origin='lower',
                          vmin=0, vmax=20)
axes[1, 2].set_title('3D Displacement Magnitude (cm)', fontsize=11, color='white')
plt.colorbar(im5, ax=axes[1, 2], shrink=0.8)

# Draw fault trace on all panels
for ax in axes.flat:
    fx = [-7.5 * np.cos(np.radians(45)), 7.5 * np.cos(np.radians(45))]
    fy = [-7.5 * np.sin(np.radians(45)), 7.5 * np.sin(np.radians(45))]
    ax.plot(fx, fy, 'k-', linewidth=2, label='Fault trace')
    ax.set_xlabel('East (km)', fontsize=9, color='white')
    ax.set_ylabel('North (km)', fontsize=9, color='white')
    ax.tick_params(colors='white', labelsize=8)

fig.patch.set_facecolor('#0a0a1a')
for ax in axes.flat:
    ax.set_facecolor('#0a0a1a')
    ax.title.set_color('white')

plt.tight_layout()
plt.savefig('output.png', dpi=130, bbox_inches='tight', facecolor='#0a0a1a', edgecolor='none')
plt.close()

# Print summary
M0 = 30e9 * 15000 * 10000 * 1.5  # mu * L * W * slip
Mw = 2/3 * (np.log10(M0) - 9.1)
print(f"Fault parameters:")
print(f"  Strike: {fault_params['strike_deg']} deg, Dip: {fault_params['dip_deg']} deg, Rake: {fault_params['rake_deg']} deg")
print(f"  Length: 15 km, Width: 10 km, Depth: {fault_params['depth']/1000} km")
print(f"  Slip: {fault_params['slip']} m")
print(f"  Seismic moment: M0 = {M0:.2e} N m")
print(f"  Moment magnitude: Mw = {Mw:.1f}")
print(f"\nMax surface displacement: {disp_mag.max()*100:.1f} cm")
print(f"Max LOS displacement: {np.abs(d_los).max()*100:.1f} cm")
print(f"Number of fringes (approx): {int(2 * np.abs(d_los).max() / (wavelength/2))}")
print(f"\nSentinel-1 C-band wavelength: {wavelength*100:.3f} cm")
print(f"One fringe = {wavelength/2*100:.3f} cm LOS displacement")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Notable InSAR Earthquake Studies

2023 Turkey-Syria Earthquakes (Mw 7.8 + 7.5)

Sentinel-1 DInSAR revealed >3 m of horizontal displacement along the East Anatolian Fault. The two-event sequence ruptured ~350 km, making it the most destructive earthquake measured by modern InSAR.

2016 Amatrice, Italy (Mw 6.2)

Sentinel-1 captured the 20 cm co-seismic subsidence within 6 days. PSInSAR time series later showed post-seismic afterslip continuing for months.

2015 Gorkha, Nepal (Mw 7.8)

ALOS-2 L-band InSAR (better coherence than C-band in vegetated terrain) measured up to 1.5 m of surface uplift south of the Himalayan range.

2019 Ridgecrest, California (Mw 7.1)

Sentinel-1 + ALOS-2 ascending/descending decomposition resolved the complex conjugate fault system with cm-level precision.

InSAR Error Sources & Mitigation

Error Source	Magnitude	Mitigation
Atmospheric delay	1-15 cm (troposphere)	GACOS correction, ERA5 integration, stacking
DEM error	Proportional to baseline	Use high-quality DEM (Copernicus 30m), short baselines
Temporal decorrelation	Loss of phase info	Short temporal baselines, L-band SAR for vegetation
Unwrapping errors	Multiples of 2pi	Multilooking, careful masking, multiple algorithms
Orbital errors	Long-wavelength ramps	Precise orbit files (POD), polynomial deramping

← M3: Data Access M5: Drought & Vegetation →