用Python+Matplotlib复现断层‘沙滩球’：从参数计算到一键成图实战

用Python+Matplotlib复现断层‘沙滩球’：从参数计算到一键成图实战

地质构造的动力学特征往往通过震源机制解来呈现，而沙滩球图（Beach Ball）正是这种三维空间信息的二维投影表达。本文将带领读者用Python实现从断层参数到可视化图形的全流程，不仅包含数学推导和代码实现，还会深入解析每个环节的工程化细节。

1. 理解沙滩球图的核心参数

沙滩球图的本质是断层运动在球面上的立体投影。要准确绘制图形，必须掌握三个核心参数：

• 走向（Strike）：断层线与正北方向的夹角（0-360°）
• 倾角（Dip）：断层面与水平面的夹角（0-90°）
• 滑移角（Rake）：断层滑动方向与走向线的夹角（-180°到180°）

这三个参数共同定义了断层类型判别矩阵：

注意：当dip=90°时，需要特别处理上/下盘的定义问题，这是许多开源库容易出错的边界条件。

2. 极射赤平投影的数学实现

赤平投影的核心是将三维球面坐标转换为二维平面坐标。我们采用下半球等角度投影（Lower Hemisphere Equal-Angle Projection），其数学变换过程如下：

import numpy as np def stereographic_projection(trend, plunge): """ 将线理/面理要素投影到赤平面 :param trend: 倾向方位角（度） :param plunge: 倾伏角（度） :return: (x, y) 投影坐标 """ theta = np.radians(90 - plunge) r = np.sqrt(2) * np.sin(theta/2) x = r * np.cos(np.radians(trend)) y = r * np.sin(np.radians(trend)) return x, y

这个基础函数可以处理：

• 断层面法线的投影
• 滑动矢量的投影
• 辅助构造线的绘制

对于断层面的投影，需要先计算其法向量：

def fault_normal_vector(strike, dip): """ 计算断层面单位法向量 """ strike_rad = np.radians(strike) dip_rad = np.radians(dip) nx = -np.sin(dip_rad) * np.sin(strike_rad) ny = np.sin(dip_rad) * np.cos(strike_rad) nz = -np.cos(dip_rad) return np.array([nx, ny, nz])

3. 沙滩球绘制的完整流程

3.1 数据预处理模块

创建参数校验函数确保输入合法：

def validate_focal_mechanism(strike, dip, rake): assert 0 <= strike < 360, "走向角超出范围" assert 0 <= dip <= 90, "倾角超出范围" assert -180 <= rake <= 180, "滑移角超出范围" # 处理360度边界条件 strike = strike % 360 rake = (rake + 180) % 360 - 180 return strike, dip, rake

3.2 图形绘制引擎

使用Matplotlib的Patch集合实现高效绘图：

from matplotlib.patches import Polygon, Wedge from matplotlib.collections import PatchCollection def draw_beachball(ax, strike, dip, rake, size=1): # 创建压缩区和膨胀区多边形 patches = [] # 计算断层节面（两个正交平面） P_plane = compute_principal_plane(strike, dip) T_plane = compute_aux_plane(strike, dip, rake) # 生成压缩区多边形 for segment in compress_zones(P_plane, T_plane): polygon = Polygon(segment, closed=True) patches.append(polygon) # 设置颜色方案 colors = ['white', 'black'] collection = PatchCollection(patches, alpha=0.8) collection.set_facecolor(colors) # 添加比例圆 border = Wedge((0,0), size, 0, 360, width=0.05) ax.add_patch(border) ax.add_collection(collection) ax.set_aspect('equal') ax.axis('off')

3.3 断层类型自动判别

基于参数组合的智能分类：

def classify_fault(strike, dip, rake): if abs(dip - 90) < 1e-6: # 垂直断层 if abs(abs(rake) - 90) < 1e-6: return "纯倾滑断层" else: return "垂直走滑断层（{}旋）".format( "左" if rake > 0 else "右") elif abs(rake) < 20 or abs(rake-180) < 20: return "走滑断层" elif -180 < rake < 0: return "正断层" else: return "逆断层"

4. 工程化扩展功能

4.1 批量处理接口

为科研场景设计的数据批处理功能：

def batch_plot_focal_mechanisms(data_frame, output_dir): """ :param data_frame: 包含strike/dip/rake列的DataFrame :param output_dir: 输出目录路径 """ os.makedirs(output_dir, exist_ok=True) for idx, row in data_frame.iterrows(): fig, ax = plt.subplots(figsize=(6,6)) draw_beachball(ax, row['strike'], row['dip'], row['rake']) # 添加标注信息 title = f"Event {idx}\n{classify_fault(**row)}" ax.set_title(title, pad=20) # 保存矢量图形 fig.savefig(f"{output_dir}/event_{idx}.svg", bbox_inches='tight') plt.close(fig)

4.2 交互式调试工具

结合IPython控件实现参数实时调整：

from ipywidgets import interact, FloatSlider @interact( strike=FloatSlider(min=0, max=359, step=1, value=45), dip=FloatSlider(min=0, max=90, step=1, value=30), rake=FloatSlider(min=-180, max=180, step=1, value=90) ) def interactive_beachball(strike, dip, rake): fig, ax = plt.subplots(figsize=(8,8)) draw_beachball(ax, strike, dip, rake) ax.set_title(classify_fault(strike, dip, rake), fontsize=14) plt.show()

5. 性能优化技巧

5.1 向量化计算

使用NumPy广播机制加速批量投影计算：

def batch_projection(trends, plunges): """ 向量化处理多个投影点 """ trends = np.asarray(trends) plunges = np.asarray(plunges) theta = np.radians(90 - plunges) r = np.sqrt(2) * np.sin(theta/2) x = r * np.cos(np.radians(trends)) y = r * np.sin(np.radians(trends)) return np.column_stack((x, y))

5.2 缓存机制

对频繁使用的三角函数结果进行缓存：

from functools import lru_cache @lru_cache(maxsize=360) def cached_cos(degree): return np.cos(np.radians(degree)) @lru_cache(maxsize=360) def cached_sin(degree): return np.sin(np.radians(degree))

6. 常见问题解决方案

图形锯齿问题：

• 增加极射投影的分段数（建议≥100）
• 使用矢量格式输出（SVG/PDF）

def smooth_projection(vertices, n_segments=100): """ 对投影边界进行平滑处理 """ from scipy.interpolate import splprep, splev tck, u = splprep(vertices.T, u=None, s=0.0, per=1) u_new = np.linspace(u.min(), u.max(), n_segments) x_new, y_new = splev(u_new, tck, der=0) return np.column_stack((x_new, y_new))

多图形对比场景： 建议使用Subplot网格配合标准化尺寸：

def plot_mechanism_comparison(events, ncols=3): nrows = (len(events) + ncols - 1) // ncols fig, axes = plt.subplots(nrows, ncols, figsize=(3*ncols, 3*nrows)) for ax, (strike, dip, rake) in zip(axes.flat, events): draw_beachball(ax, strike, dip, rake, size=0.9) ax.set_title(f"{strike}/{dip}/{rake}", fontsize=9) for ax in axes.flat[len(events):]: ax.axis('off') plt.tight_layout() return fig