用Python绘制地震波时距曲线：从理论公式到动态可视化

地震波时距曲线是地球物理勘探中的核心分析工具，它能将抽象的地下结构转化为可量化的时间-距离关系。对于从事地震数据处理的研究人员和工程师而言，掌握时距曲线的编程实现不仅有助于理解波场传播规律，更能为后续的反演解释工作奠定基础。本文将使用Python科学计算栈，带您逐步构建三种典型地震波（直达波、反射波、折射波）的时距曲线模拟系统，并通过交互式可视化展现地质参数变化对曲线形态的影响。

1. 环境配置与基础理论

1.1 Python科学计算环境搭建

推荐使用Anaconda创建专属的地球物理分析环境：

conda create -n seismology python=3.9
conda activate seismology
conda install numpy matplotlib ipywidgets scipy

关键库的作用说明：

• NumPy ：处理数组运算与数学公式实现
• Matplotlib ：生成出版级质量的二维图表
• IPyWidgets ：创建交互式参数调节面板
• SciPy ：提供高级数学函数和优化工具

提示：Jupyter Notebook/Lab非常适合本类数值模拟实验，能实时查看代码运行结果。

1.2 时距曲线物理模型精要

时距曲线描述地震波传播时间t与炮检距x的函数关系，其形态由波类型和地层结构共同决定：

波类型	典型方程	几何形态	关键参数
直达波	t = x/v	直线	波速v
反射波	t² = t₀² + x²/v²	双曲线	自激自收时间t₀、波速v
折射波	t = x/v₂ + 2hcosθ/v₁	直线	界面深度h、临界角θ

其中折射波的临界角θ满足Snell定律：sinθ = v₁/v₂。这些方程将成为我们后续编程实现的理论基础。

2. 直达波时距曲线实现

2.1 均匀介质中的波传播建模

在均匀各向同性介质中，直达波传播路径最简单，其Python实现仅需基本数组运算：

import numpy as np
import matplotlib.pyplot as plt

def direct_wave(x, v):
    """计算直达波旅行时"""
    return x / v

# 参数设置
v = 2000  # 波速(m/s)
x = np.linspace(0, 5000, 100)  # 炮检距序列(m)

# 计算并绘图
t = direct_wave(x, v)
plt.plot(x, t, 'r-', lw=2, label=f'直达波 (v={v}m/s)')
plt.xlabel('炮检距 (m)'); plt.ylabel('旅行时 (s)')
plt.legend(); plt.grid()

这段代码揭示了两个重要特性：

1. 曲线斜率与波速成反比
2. 零偏移距处时间必然为零

2.2 多速度层对比分析

通过列表循环可快速比较不同速度模型的曲线差异：

velocity_range = [1500, 2000, 2500]  # 典型地层速度(m/s)
for v in velocity_range:
    plt.plot(x, direct_wave(x, v), label=f'v={v}m/s')

实际应用中，我们常通过拟合野外实测直达波数据来估算近地表速度。这种速度信息对后续的静校正处理至关重要。

3. 反射波时距曲线模拟

3.1 水平界面双曲线实现

水平层状介质的反射波时距曲线呈现典型双曲线特征，其实现需要平方根运算：

def reflection_wave(x, v, t0):
    """水平界面反射波旅行时"""
    return np.sqrt(t0**2 + x**2/v**2)

# 示例参数
v = 3000  # 上覆地层速度(m/s)
h = 1500  # 界面深度(m)
t0 = 2*h/v  # 自激自收时间(s)

# 计算曲线
t_ref = reflection_wave(x, v, t0)
plt.plot(x, t_ref, 'b--', label=f'反射波 h={h}m')

关键观察点：

• 曲线极小点对应t₀时间
• 远偏移距渐近线斜率趋近1/v
• 双曲率随深度增加而减小

3.2 倾斜界面修正模型

当反射界面倾斜时，需引入倾角ϕ修正方程：

def tilted_reflection(x, v, h, phi):
    """倾斜界面反射波旅行时"""
    phi_rad = np.deg2rad(phi)
    return np.sqrt(x**2 + 4*h**2 + 4*h*x*np.sin(phi_rad)) / v

# 对比不同倾角效果
for phi in [0, 15, 30]:
    t_tilt = tilted_reflection(x, v, h, phi)
    plt.plot(x, t_tilt, label=f'倾角{phi}°')

倾斜界面会导致双曲线极小点偏移激发点正下方，这种偏移量与界面倾角呈正相关。实际资料处理时需要特别关注此类现象以避免构造解释错误。

4. 折射波时距曲线构建

4.1 临界折射的数学表达

折射波时距曲线实现需要计算临界角并处理截距时间：

def refraction_wave(x, v1, v2, h):
    """折射波旅行时计算"""
    theta_c = np.arcsin(v1/v2)  # 临界角
    intercept = 2*h*np.cos(theta_c)/v1  # 截距时间
    return x/v2 + intercept

# 两层介质模型
v1, v2 = 2000, 3500  # 上下层速度(m/s)
h = 800  # 界面深度(m)

# 生成曲线
t_refr = refraction_wave(x, v1, v2, h)
plt.plot(x, t_refr, 'g-.', label='折射波')

折射波曲线的两个关键特征：

1. 直线斜率对应下层介质速度倒数
2. 与直达波的交点称为临界距离，可用于估算界面深度

4.2 实际应用中的盲区问题

折射波存在一个重要的"盲区"现象——当接收点距震源过近时无法观测到折射波。通过修改代码可以清晰展示这一特性：

x_critical = 2*h*np.tan(np.arcsin(v1/v2))  # 盲区临界距离
x_refr = x[x >= x_critical]  # 只计算盲区外的数据
t_valid = refraction_wave(x_refr, v1, v2, h)
plt.plot(x_refr, t_valid, 'm:', lw=3)

这个特性提醒我们在野外观测系统设计时，需要合理布置检波器排列以避免漏失关键数据。

5. 综合可视化与参数分析

5.1 交互式曲线对比工具

利用IPyWidgets创建动态调节面板，直观理解参数影响：

from ipywidgets import interact

@interact(v=(1000,5000,100), h=(500,3000,100), phi=(0,45,5))
def plot_curves(v=2000, h=1500, phi=0):
    t_dir = direct_wave(x, v)
    t_ref = tilted_reflection(x, v, h, phi)
    t_refr = refraction_wave(x, v, v*1.8, h)
    
    plt.figure(figsize=(10,6))
    plt.plot(x, t_dir, label='直达波')
    plt.plot(x, t_ref, label='反射波')
    plt.plot(x[x>=2*h*0.58], t_refr[x>=2*h*0.58], label='折射波')
    plt.legend(); plt.grid()

这种交互式工具特别适合用于教学演示，学员可通过滑动条实时观察各参数对曲线形态的影响。

5.2 实际数据拟合案例

假设我们有一组野外采集的初至时间数据：

field_x = np.array([0,500,1000,1500,2000,2500])  # 炮检距(m)
field_t = np.array([0,0.28,0.55,0.82,1.08,1.35])  # 观测时间(s)

# 最小二乘速度拟合
from scipy.optimize import curve_fit
popt, _ = curve_fit(direct_wave, field_x, field_t)
v_est = popt[0]  # 估计波速

通过此类拟合分析，我们可以将理论模型与实际观测相结合，验证地下结构的假设是否合理。这种正演模拟能力是地震资料解释人员的重要技能。

在地震勘探实践中，时距曲线分析只是整个工作流程的起点。现代地震处理软件虽然已经高度自动化，但理解这些基础原理仍然至关重要——它们是我们判断处理结果合理性的最后防线。当我在处理某次海上地震数据时，就曾通过时距曲线异常发现了一个被忽略的速度异常体，这再次证明了基础理论不可替代的价值。