从‘沙滩球’反推断层运动：手把手用Python绘制震源机制解

从‘沙滩球’反推断层运动：手把手用Python绘制震源机制解

地震学研究中最直观的工具莫过于震源机制解图示——那个黑白相间、形似沙滩球的图案。这种专业图表不仅能展示断层的三维运动特征，还能帮助研究者快速判断地震类型。本文将带您用Python从零实现沙滩球的可视化与逆向解析，让抽象的地震参数变成可交互的代码实践。

1. 环境准备与基础概念

在开始编码前，我们需要明确几个核心参数的定义。断层的几何特征由三个关键角度决定：

• 走向(strike)：断层线与正北方向的夹角（0-360°）
• 倾角(dip)：断层面与水平面的夹角（0-90°）
• 滑移角(rake)：断层滑动方向与走向线的夹角（-180°到180°）

安装必要的Python库（建议使用conda环境）：

conda create -n focal_mech python=3.8 conda install -c conda-forge obspy matplotlib numpy

注意：Obspy库是地震学分析的瑞士军刀，其focal_mechanism模块已内置沙滩球绘制功能

2. 正向工程：参数转沙滩球

让我们从一个实际案例开始。假设某次地震的断层参数为：走向45°，倾角30°，滑移角80°。用matplotlib绘制基础沙滩球：

import matplotlib.pyplot as plt from obspy.imaging.beachball import beachball # 输入参数：走向、倾角、滑移角 focal_mechanism = [45, 30, 80] fig = plt.figure(figsize=(6,6)) ax = fig.add_subplot(111, projection='3d') beachball(focal_mechanism, size=200, linewidth=2, ax=ax) plt.title('震源机制解示例', pad=20) plt.tight_layout() plt.show()

执行后会生成标准的赤平投影图，其中：

• 白色区域：压缩象限（P波初动向上）
• 黑色区域：膨胀象限（P波初动向下）
• 边界线：代表断层面的交线

3. 逆向工程：图像反演参数

实际研究中更常见的情况是：我们获得了一个沙滩球图像或地震矩张量，需要反推断层参数。以下代码演示如何从矩张量反解：

from obspy.imaging.beachball import MT2Plane # 假设已知矩张量分量（单位：N·m） moment_tensor = [1.23e+18, -2.45e+17, -9.87e+17, 5.32e+17, 3.21e+17, -7.65e+16] # 获取可能的断层参数解（通常有两个解） fault_planes = MT2Plane(moment_tensor) for i, plane in enumerate(fault_planes, 1): print(f"解{i}: 走向={plane[0]:.1f}°, 倾角={plane[1]:.1f}°, 滑移角={plane[2]:.1f}°")

典型输出示例：

解1: 走向=125.3°, 倾角=42.7°, 滑移角=78.5° 解2: 走向=302.6°, 倾角=48.3°, 滑移角=102.1°

提示：两个解在数学上等效，需结合地质背景判断哪个更合理

4. 断层类型自动判别系统

基于反演得到的参数，我们可以编写自动分类程序：

def classify_fault(strike, dip, rake): if dip == 90: # 垂直断层 if abs(rake) == 90: return "垂直倾滑断层" elif rake == 0: return "左旋走滑断层" elif abs(rake) == 180: return "右旋走滑断层" elif 0 < rake < 180: return "逆断层" elif -180 < rake < 0: return "正断层" elif abs(rake) == 90: return "倾滑断层" else: return "混合型断层" # 使用前文反演结果测试 for plane in fault_planes: fault_type = classify_fault(*plane) print(f"走向{plane[0]:.0f}°/{plane[1]:.0f}°/{plane[2]:.0f}° → {fault_type}")

5. 高级可视化技巧

为了让结果更专业，我们可以增强可视化效果：

import numpy as np from mpl_toolkits.mplot3d import Axes3D fig = plt.figure(figsize=(12,5)) # 子图1：标准沙滩球 ax1 = fig.add_subplot(121, projection='3d') beachball(focal_mechanism, size=100, ax=ax1) ax1.set_title('赤平投影', y=1.05) # 子图2：三维断层示意 ax2 = fig.add_subplot(122, projection='3d') strike, dip, rake = np.radians(focal_mechanism) # 绘制断层面 x = np.linspace(0, 1, 10) y = np.linspace(0, np.cos(dip), 10) X, Y = np.meshgrid(x, y) Z = Y * np.tan(dip) ax2.plot_surface(X, Y, Z, alpha=0.5, color='blue') # 添加滑动方向箭头 slip_x = 0.5 * np.cos(rake) slip_y = 0.5 * np.sin(rake) * np.cos(dip) ax2.quiver(0.5, 0, 0, slip_x, slip_y, 0, color='red', length=0.3, arrow_length_ratio=0.1) ax2.set_title('三维断层模型') ax2.set_xlim(0,1); ax2.set_ylim(0,1); ax2.set_zlim(0,1) plt.tight_layout()

6. 实战案例：真实地震分析

以2023年土耳其地震为例（USGS提供的矩张量）：

# 土耳其地震矩张量（M7.8） turkey_mt = [1.97, -0.57, -1.40, 0.71, -0.13, -0.47] # ×1e20 N·m # 反演断层参数 planes = MT2Plane([x*1e20 for x in turkey_mt]) for i, (strike, dip, rake) in enumerate(planes, 1): print(f"解{i}: {strike:.1f}°/{dip:.1f}°/{rake:.1f}° → {classify_fault(strike, dip, rake)}") # 可视化比较 fig, axes = plt.subplots(1,2, figsize=(10,5), subplot_kw={'projection':'3d'}) for ax, plane in zip(axes, planes): beachball(plane, size=200, ax=ax) ax.set_title(f"{classify_fault(*plane)}\n{plane[0]:.0f}°/{plane[1]:.0f}°/{plane[2]:.0f}°")

输出结果将显示这是典型的左旋走滑断层，与实际地质调查结论一致。通过调整代码中的矩张量参数，可以分析任何历史地震事件。