用Python动态绘制地震波时距曲线从理论到可视化实战地震勘探的核心在于理解波在不同介质中的传播规律而时距曲线正是这种规律的直观体现。传统教学中学生往往需要死记硬背各种曲线方程却难以建立直观认知。本文将带您用Python的Matplotlib和NumPy库通过代码实现从水平界面到倾斜地层的多种时距曲线动态生成让抽象公式变为可交互的图形。1. 环境配置与基础准备在开始前确保已安装Python 3.8环境及以下库pip install numpy matplotlib ipywidgets基础参数设置是模拟的起点。我们首先定义全局变量import numpy as np import matplotlib.pyplot as plt from ipywidgets import interact # 基础参数 v1 1500 # 上覆地层速度(m/s) v2 3000 # 下伏地层速度(m/s) h 500 # 界面深度(m) theta_c np.arcsin(v1/v2) # 临界角(弧度)注意速度单位需保持一致建议使用m/s。界面深度h可根据实际勘探目标调整典型陆上勘探范围在200-3000米之间。2. 水平地层模型的三类基本曲线2.1 直达波最简单的直线模型直达波时距曲线是理解波传播的基础案例。在均匀介质中其方程为$$ t \frac{x}{v} $$对应Python实现def direct_wave(x, v): 计算直达波旅行时 return np.abs(x) / v # 生成曲线数据 x np.linspace(0, 3000, 100) # 炮检距范围0-3000米 t_direct direct_wave(x, v1) # 可视化 plt.figure(figsize(10,6)) plt.plot(x, t_direct, labelf直达波 (v{v1}m/s)) plt.xlabel(炮检距 (m)); plt.ylabel(旅行时 (s)) plt.legend(); plt.grid()关键观察点曲线斜率直接反映地层速度零偏移距处时间为0符合物理实际双边接收时表现为通过原点的两条直线2.2 反射波双曲线特征与参数敏感性水平界面反射波时距曲线方程为$$ t \frac{\sqrt{x^2 4h^2}}{v_1} $$代码实现中加入参数交互功能def reflection_wave(x, h, v): 计算反射波旅行时 return np.sqrt(x**2 4*h**2) / v interact(h(100, 1000, 50), v(1000, 5000, 100)) def plot_reflection(h500, v2000): t_ref reflection_wave(x, h, v) plt.figure(figsize(10,6)) plt.plot(x, t_ref, labelf反射波 h{h}m, v{v}m/s) plt.xlabel(炮检距 (m)); plt.ylabel(旅行时 (s)) plt.legend(); plt.grid()通过滑块可发现深度h主要影响曲线弯曲程度速度v改变整体时间尺度极小点t0 2h/v是重要解释参数2.3 折射波临界角后的线性传播折射波时距曲线方程较为复杂$$ t \frac{2h\cos\theta_c}{v_1} \frac{x}{v_2} $$其Python实现需要考虑临界距离def refraction_wave(x, h, v1, v2): 计算折射波旅行时 theta_c np.arcsin(v1/v2) x_crit 2 * h * np.tan(theta_c) # 临界距离 t np.where(x x_crit, 2*h*np.cos(theta_c)/v1 x/v2, np.nan) # 小于临界距离时不产生折射波 return t t_refr refraction_wave(x, h, v1, v2) # 绘制三种曲线对比 plt.figure(figsize(12,7)) plt.plot(x, t_direct, label直达波) plt.plot(x, reflection_wave(x,h,v1), label反射波) plt.plot(x, t_refr, label折射波) plt.xlabel(炮检距 (m)); plt.ylabel(旅行时 (s)) plt.legend(); plt.grid()典型特征包括存在最小临界距离后才出现折射波直线斜率反映下层介质速度v2与直达波存在交叉点可用于估算界面深度3. 倾斜界面模型的扩展实现当地层倾斜时反射波时距曲线方程变为$$ t \frac{\sqrt{x^2 4h^2 4hx\sin\phi}}{v_1} $$Python实现需扩展参数def tilted_reflection(x, h, v, phi): 倾斜界面反射波 phi_rad np.deg2rad(phi) # 角度转弧度 return np.sqrt(x**2 4*h**2 4*h*x*np.sin(phi_rad)) / v interact(phi(0, 45, 5)) def compare_tilt(phi10): t_horiz reflection_wave(x, h, v1) t_tilt tilted_reflection(x, h, v1, phi) plt.figure(figsize(12,7)) plt.plot(x, t_horiz, label水平界面) plt.plot(x, t_tilt, labelf倾斜界面 φ{phi}°) plt.xlabel(炮检距 (m)); plt.ylabel(旅行时 (s)) plt.legend(); plt.grid()倾斜效应表现为双曲线极小点偏离零偏移距位置上倾方向与下倾方向不对称倾角越大极小点偏移越明显4. 实际应用与解释技巧4.1 速度分析实用方法利用反射波时距曲线估算速度的常用方法def velocity_analysis(t, x, t0): 速度谱计算 dt t**2 - t0**2 v_est np.sqrt(x**2 / dt) # 由双曲线方程变形 return v_est # 示例从合成数据反演速度 t_syn reflection_wave(x, h, v1) np.random.normal(0, 0.001, len(x)) # 添加噪声 t0 2*h/v1 v_estimated velocity_analysis(t_syn, x, t0) print(f真实速度: {v1}m/s, 估计速度: {np.mean(v_estimated[1:]):.1f}m/s)提示实际数据需进行动校正处理上述简化方法仅适用于理想情况4.2 层状介质建模进阶对于多层介质时距曲线可通过迭代计算实现def multi_layer_reflection(x, layers): layers [(厚度1,速度1), (厚度2,速度2), ...] t_total np.zeros_like(x) h_accum 0 for i, (thickness, v) in enumerate(layers): h_accum thickness t_layer reflection_wave(x, h_accum, v) t_total t_layer if i 0 else t_layer - reflection_wave(x, h_accum-thickness, v) return t_total # 三层模型示例 layers [(300, 1800), (400, 2200), (500, 3000)] t_multi multi_layer_reflection(x, layers) plt.figure(figsize(12,7)) plt.plot(x, t_multi) plt.title(三层介质反射波时距曲线); plt.grid()复杂介质建模要点每层贡献的时差需累加计算深层反射曲线更平缓实际应用中需考虑波阻抗差异5. 交互式模拟系统开发整合所有功能创建交互界面from ipywidgets import FloatSlider, Dropdown, VBox # 控件定义 style {description_width: 150px} h_slider FloatSlider(min100, max2000, step50, value500, description界面深度(m), stylestyle) v1_slider FloatSlider(min1000, max5000, step100, value2000, description上覆层速度(m/s), stylestyle) v2_slider FloatSlider(min1500, max6000, step100, value3500, description下伏层速度(m/s), stylestyle) phi_slider FloatSlider(min0, max45, step5, value0, description界面倾角(度), stylestyle) wave_type Dropdown(options[直达波, 反射波, 折射波, 全显示], value全显示, description波形选择) def update_plot(h, v1, v2, phi, wave_type): x np.linspace(0, 3000, 100) plt.figure(figsize(14,8)) if wave_type in [直达波, 全显示]: t_dir direct_wave(x, v1) plt.plot(x, t_dir, b, labelf直达波 v{v1}m/s) if wave_type in [反射波, 全显示]: t_ref tilted_reflection(x, h, v1, phi) if phi 0 else reflection_wave(x, h, v1) plt.plot(x, t_ref, r, labelf反射波 h{h}m (f, φ{phi}° if phi 0 else )) if wave_type in [折射波, 全显示] and v2 v1: t_refr refraction_wave(x, h, v1, v2) plt.plot(x, t_refr, g, labelf折射波 v2{v2}m/s) plt.xlabel(炮检距 (m)); plt.ylabel(旅行时 (s)) plt.legend(); plt.grid() plt.title(地震波时距曲线模拟系统) controls VBox([h_slider, v1_slider, v2_slider, phi_slider, wave_type]) interact(update_plot, hh_slider, v1v1_slider, v2v2_slider, phiphi_slider, wave_typewave_type);这个交互系统允许实时调整地质参数观察曲线变化单独显示或对比不同波型直观理解各参数的地震响应特征