Python实战用代码绘制地震波时距曲线的完整指南地震勘探数据的可视化是理解地下结构的关键步骤。本文将带您从零开始使用Python的Matplotlib和NumPy库完整实现直达波、反射波和折射波时距曲线的绘制过程。不同于传统理论讲解我们将通过可运行的代码示例让抽象的地震波传播规律变得直观可见。1. 环境准备与基础概念在开始编码前我们需要明确几个核心概念。时距曲线描述的是地震波传播时间(t)与炮检距(x)之间的关系不同类型的波会形成不同形状的曲线直达波直线反射波双曲线折射波直线不同斜率首先设置Python环境import numpy as np import matplotlib.pyplot as plt plt.style.use(seaborn) # 使用更美观的绘图样式基础参数设置后续可调整v 2000 # 波速(m/s) h 1000 # 反射界面深度(m) theta_c 30 # 临界角(度) x_max 3000 # 最大炮检距(m)2. 直达波时距曲线实现直达波是最简单的时距曲线其方程为t x / v代码实现def direct_wave(x, v): 计算直达波传播时间 return x / v x_values np.linspace(0, x_max, 100) t_direct direct_wave(x_values, v) plt.figure(figsize(10,6)) plt.plot(x_values, t_direct, label直达波) plt.xlabel(炮检距(m)) plt.ylabel(传播时间(s)) plt.title(直达波时距曲线) plt.legend() plt.grid(True)关键观察点曲线斜率为1/v通过改变v值可直观观察斜率变化实际应用中v值可通过近地表测量获得3. 反射波时距曲线绘制水平界面反射波时距曲线方程为双曲线t √(x² 4h²) / v代码实现def reflection_wave(x, v, h): 计算反射波传播时间 return np.sqrt(x**2 4*h**2) / v t_reflect reflection_wave(x_values, v, h) plt.figure(figsize(10,6)) plt.plot(x_values, t_direct, label直达波) plt.plot(x_values, t_reflect, label反射波) plt.xlabel(炮检距(m)) plt.ylabel(传播时间(s)) plt.title(反射波与直达波时距曲线对比) plt.legend() plt.grid(True)交互实验尝试修改以下参数观察曲线变化# 尝试不同深度和速度 h_values [800, 1000, 1200] v_values [1800, 2000, 2200] plt.figure(figsize(12,8)) for h_val in h_values: t reflection_wave(x_values, v, h_val) plt.plot(x_values, t, labelfh{h_val}m) plt.xlabel(炮检距(m)) plt.ylabel(传播时间(s)) plt.title(不同深度反射波时距曲线) plt.legend() plt.grid(True)4. 折射波时距曲线实现折射波时距曲线方程为直线t x/v2 (2h cosθc)/v1代码实现def refraction_wave(x, v1, v2, h, theta_c): 计算折射波传播时间 theta_rad np.deg2rad(theta_c) intercept (2 * h * np.cos(theta_rad)) / v1 return x/v2 intercept v2 3000 # 下层介质波速 t_refract refraction_wave(x_values, v, v2, h, theta_c) plt.figure(figsize(10,6)) plt.plot(x_values, t_direct, label直达波) plt.plot(x_values, t_reflect, label反射波) plt.plot(x_values, t_refract, label折射波) plt.xlabel(炮检距(m)) plt.ylabel(传播时间(s)) plt.title(三种波型时距曲线对比) plt.legend() plt.grid(True)参数敏感性分析参数影响典型值范围v1上层介质速度1500-2500 m/sv2下层介质速度2500-4500 m/sh界面深度500-3000 mθc临界角20-40度5. 高级可视化与交互功能为了让图形更具专业性和交互性我们添加以下功能from ipywidgets import interact def plot_curves(v2000, h1000, v23000, theta_c30): x np.linspace(0, 3000, 100) plt.figure(figsize(12,8)) plt.plot(x, direct_wave(x, v), label直达波) plt.plot(x, reflection_wave(x, v, h), label反射波) plt.plot(x, refraction_wave(x, v, v2, h, theta_c), label折射波) plt.xlabel(炮检距(m), fontsize12) plt.ylabel(传播时间(s), fontsize12) plt.title(地震波时距曲线综合展示, fontsize14) plt.legend(fontsize12) plt.grid(True) plt.tight_layout() interact(plot_curves, v(1500, 2500, 50), h(500, 2000, 50), v2(2500, 4500, 50), theta_c(20, 40, 1))图形美化技巧# 专业论文级图形输出示例 plt.figure(figsize(12,8), dpi300) plt.plot(x_values, t_direct, linewidth2, label直达波) plt.plot(x_values, t_reflect, linewidth2, label反射波) plt.plot(x_values, t_refract, linewidth2, label折射波) # 添加标注 plt.annotate(交叉时, xy(1500, 0.8), xytext(1800,0.7), arrowpropsdict(facecolorblack, shrink0.05)) plt.xlabel(炮检距 (m), fontsize14) plt.ylabel(传播时间 (s), fontsize14) plt.title(地震波时距曲线特征对比, fontsize16) plt.legend(fontsize12, frameonTrue) plt.grid(True, linestyle--, alpha0.7) plt.tight_layout()6. 实际应用案例分析通过一个完整的地层模型展示实际应用# 定义三层地层模型 layers [ {v: 1800, h: 500}, # 第一层 {v: 2200, h: 800}, # 第二层 {v: 3000, h: 1200} # 第三层 ] def multi_layer_reflection(x, layers): 计算多层反射波时距曲线 t0 0 v_rms 0 total_depth 0 # 计算均方根速度 for layer in layers: t0 2 * layer[h] / layer[v] v_rms layer[v]**2 * (2*layer[h]/layer[v]) total_depth layer[h] v_rms np.sqrt(v_rms / t0) return np.sqrt(t0**2 x**2 / v_rms**2) x np.linspace(0, 5000, 200) t_multi multi_layer_reflection(x, layers) plt.figure(figsize(12,8)) plt.plot(x, t_multi, label多层反射波) plt.xlabel(炮检距(m)) plt.ylabel(传播时间(s)) plt.title(复杂地层时距曲线示例) plt.legend() plt.grid(True)数据处理技巧使用scipy.optimize进行曲线拟合添加随机噪声模拟实际数据使用移动平均平滑曲线from scipy.optimize import curve_fit # 添加噪声的模拟数据 noise np.random.normal(0, 0.01, len(x)) t_noisy t_multi noise # 定义拟合函数 def fit_func(x, v_rms, t0): return np.sqrt(t0**2 x**2 / v_rms**2) params, _ curve_fit(fit_func, x, t_noisy) v_rms_est, t0_est params print(f估计参数: v_rms{v_rms_est:.2f} m/s, t0{t0_est:.4f} s)通过这个完整的实现过程我们不仅掌握了时距曲线的绘制方法更重要的是理解了参数变化对曲线形态的影响规律。在实际项目中这些代码可以直接应用于地震数据解释和速度分析工作。