从地震波旅行时间到地下构造：手把手教你解读时距曲线图（附Python可视化代码）

当地震波穿越地层时，它们的传播时间记录着地下结构的秘密。对于地球物理和地质工程领域的研究者来说，时距曲线图就像是一把打开地下世界大门的钥匙。本文将带你从零开始，通过Python代码实现直达波和反射波时距曲线的可视化，并教会你如何从这些曲线中解读出有价值的地下信息。

1. 时距曲线基础与Python环境准备

时距曲线（Time-Distance Curve）描述的是地震波从震源出发，传播到不同距离检波器所需的时间关系。这种曲线能够直观反映地下介质的物理特性，是地震勘探中不可或缺的分析工具。

1.1 安装必要的Python库

在开始之前，我们需要准备Python环境。推荐使用Anaconda创建虚拟环境：

conda create -n seismic python=3.8
conda activate seismic
pip install numpy matplotlib scipy

这三个库将构成我们分析的基础：

• NumPy ：处理数值计算
• Matplotlib ：数据可视化
• SciPy ：科学计算工具

1.2 基本物理概念与数学模型

直达波的时距关系最为简单，其数学表达式为：

t = x/v

其中：

• t ：旅行时间（秒）
• x ：炮检距（米）
• v ：地层波速（米/秒）

在Python中，我们可以这样实现：

import numpy as np
import matplotlib.pyplot as plt

def direct_wave(x, v):
    return x / v

# 示例参数
v = 2000  # 波速 2000 m/s
x = np.linspace(0, 5000, 100)  # 炮检距 0-5000米
t_direct = direct_wave(x, v)

plt.plot(x, t_direct)
plt.xlabel('Offset (m)')
plt.ylabel('Travel Time (s)')
plt.title('Direct Wave Time-Distance Curve')
plt.grid()
plt.show()

2. 水平界面反射波时距曲线分析

反射波时距曲线比直达波复杂得多，它能够揭示地下界面的深度和形态。对于单一水平界面，反射波时距曲线呈现双曲线特征。

2.1 双曲线特征与物理意义

水平界面反射波时距曲线方程为：

t = √(t₀² + x²/v²)

其中：

• t₀ = 2h/v 为零偏移距时间
• h 为界面深度
• v 为上覆地层速度

关键特征点：

1. 顶点 ：位于x=0处，时间值为t₀
2. 渐近线 ：斜率±1/v，反映地层速度
3. 曲率 ：与界面深度相关

def reflection_wave(x, v, h):
    t0 = 2 * h / v
    return np.sqrt(t0**2 + (x/v)**2)

# 示例参数
h = 1000  # 界面深度 1000米
t_reflect = reflection_wave(x, v, h)

plt.plot(x, t_direct, label='Direct Wave')
plt.plot(x, t_reflect, label='Reflection Wave')
plt.legend()
plt.xlabel('Offset (m)')
plt.ylabel('Travel Time (s)')
plt.title('Comparison of Wave Time-Distance Curves')
plt.grid()
plt.show()

2.2 从曲线解读地下信息

通过分析生成的时距曲线，我们可以提取以下关键参数：

曲线特征	物理意义	计算方法
顶点时间t₀	界面深度	h = v*t₀/2
渐近线斜率	地层速度	v = Δx/Δt
曲率半径	界面起伏	二阶导数分析

注意：实际解释时需要考虑波速随深度的变化，上述公式适用于均匀介质假设。

3. 倾斜界面反射波时距曲线

当地层界面存在倾斜时，时距曲线会表现出不对称性，这种不对称性正反映了地层的倾角信息。

3.1 数学模型与Python实现

倾斜界面反射波时距曲线方程为：

t = √(t₀² + x²/v² - 2xt₀sinξ/v)

其中ξ为界面倾角。Python实现如下：

def tilted_reflection(x, v, h, xi):
    xi_rad = np.deg2rad(xi)  # 角度转弧度
    t0 = 2 * h / v
    return np.sqrt(t0**2 + (x/v)**2 - 2*x*t0*np.sin(xi_rad)/v)

# 示例参数
xi = 15  # 倾角15度
t_tilted = tilted_reflection(x, v, h, xi)

plt.plot(x, t_reflect, label='Horizontal Interface')
plt.plot(x, t_tilted, label=f'Tilted Interface ({xi}°)')
plt.legend()
plt.xlabel('Offset (m)')
plt.ylabel('Travel Time (s)')
plt.title('Effect of Interface Dip on Time-Distance Curve')
plt.grid()
plt.show()

3.2 倾角提取方法

从倾斜界面时距曲线中提取倾角的主要步骤：

1. 确定极小点位置xₘ
2. 计算虚震源位置
3. 根据几何关系求解倾角ξ

# 寻找极小点位置
min_idx = np.argmin(t_tilted)
x_min = x[min_idx]
t_min = t_tilted[min_idx]

# 计算视倾角
apparent_dip = np.arcsin(x_min * v / (2 * h**2 + x_min**2))
true_dip = np.arcsin(np.tan(apparent_dip))

print(f"Estimated dip angle: {np.rad2deg(true_dip):.1f}°")

4. 实际应用与案例分析

时距曲线分析在地震勘探中有广泛应用，从简单的构造解释到复杂的储层描述都离不开它。

4.1 速度分析

通过分析时距曲线的形状可以估算地层速度，常用的方法包括：

• 正常时差分析 ：利用不同偏移距的时间差
• 速度谱 ：通过扫描速度参数寻找最佳拟合

# 速度扫描示例
velocity_range = np.linspace(1500, 2500, 11)
plt.figure(figsize=(10,6))

for v_test in velocity_range:
    t_test = reflection_wave(x, v_test, h)
    plt.plot(x, t_test, label=f'v={v_test:.0f} m/s')

plt.legend(bbox_to_anchor=(1.05, 1))
plt.title('Velocity Analysis via Curve Fitting')
plt.xlabel('Offset (m)')
plt.ylabel('Travel Time (s)')
plt.grid()
plt.tight_layout()
plt.show()

4.2 实际数据处理流程

典型的地震数据处理流程中，时距曲线分析是关键环节：

1. 数据采集 ：获取野外地震记录
2. 初至拾取 ：确定各道初至时间
3. 曲线拟合 ：用数学模型拟合观测数据
4. 参数反演 ：提取地下介质参数
5. 地质解释 ：将物理参数转化为地质认识

提示：实际工作中需要考虑噪声、多次波等干扰因素，通常需要结合其他处理方法。

5. 高级主题与扩展应用

当时距曲线分析结合现代计算技术，可以解决更复杂的地球物理问题。

5.1 各向异性介质中的时距曲线

当地层表现出各向异性时，时距曲线方程需要引入更多参数：

def anisotropic_reflection(x, v, h, epsilon, delta):
    t0 = 2 * h / v
    return np.sqrt(t0**2 + (x/v)**2 * (1 + 2*delta - 8*epsilon*(x/(v*t0))**2))

# Thomsen参数示例
eps, delta = 0.2, 0.1
t_aniso = anisotropic_reflection(x, v, h, eps, delta)

plt.plot(x, t_reflect, label='Isotropic')
plt.plot(x, t_aniso, label='Anisotropic')
plt.legend()
plt.title('Anisotropy Effect on Time-Distance Curve')
plt.xlabel('Offset (m)')
plt.ylabel('Travel Time (s)')
plt.grid()
plt.show()

5.2 三维时距曲面

对于三维勘探，时距关系扩展为曲面：

from mpl_toolkits.mplot3d import Axes3D

X, Y = np.meshgrid(np.linspace(-5000,5000,50), np.linspace(-5000,5000,50))
R = np.sqrt(X**2 + Y**2)
T = reflection_wave(R, v, h)

fig = plt.figure(figsize=(10,7))
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, Y, T, cmap='viridis')
ax.set_xlabel('X Offset (m)')
ax.set_ylabel('Y Offset (m)')
ax.set_zlabel('Travel Time (s)')
ax.set_title('3D Time-Distance Surface')
plt.show()

6. 常见问题与解决技巧

在实际工作中，时距曲线分析常会遇到各种挑战。以下是一些经验分享：

数据质量问题 ：

• 初至时间拾取不准：尝试不同的拾取算法
• 噪声干扰：应用滤波处理
• 静校正问题：检查高程和低速带校正

解释陷阱 ：

• 速度-深度歧义：结合其他资料约束
• 各向异性影响：考虑更复杂的模型
• 多次波干扰：识别和压制多次波

代码优化技巧 ：

# 使用向量化运算加速计算
def fast_reflection(x, v, h):
    x = np.asarray(x)  # 确保输入为数组
    t0 = 2 * h / v
    return np.sqrt(t0**2 + (x/v)**2)

# 对大炮检距使用渐近近似
def approx_reflection(x, v, h, threshold=3000):
    x = np.asarray(x)
    t0 = 2 * h / v
    exact = np.sqrt(t0**2 + (x/v)**2)
    approx = np.abs(x)/v  # 大炮检距渐近线
    return np.where(np.abs(x) > threshold, approx, exact)

在地震解释项目中，我发现将时距曲线分析与地质露头资料结合，能显著提高解释的可靠性。特别是在复杂构造区，单纯依靠地震数据往往难以确定真实的地层产状。