Python快速可视化SEG-Y地震数据：三步实现剖面绘制与批量抽检

简介这份资源面向地质勘探、地震预测与环境监测领域的技术人员及Qt开发者提供一套基于Qt3.3框架的地震数据快速可视化实现方案重点解决大规模地震记录的读取、处理与图形化展示问题。压缩包共77个文件约2.07MB以h头文件、cpp源文件、ui界面文件、png与xpm图标资源为主另含dsp、pro、dsw等工程配置及可执行程序与依赖库结构完整便于直接编译运行。项目支持SEGY、GRISYS、DSK等常见格式并允许自定义数据格式内置波形、充填、便面积、彩图及插值彩图五种显示方式分别适用于事件识别、异常检测、趋势分析与多维数据关联观察。已有907人学习下载读者可从中获取完整视图组件、数据处理模块与交互逻辑理解海量地震数据下的流畅渲染与界面设计思路适合作为地震数据可视化工具开发与二次扩展的参考。1. 地震数据的快速可视化从 SEG-Y 到屏幕只要三步手头有一批 SEG-Y 格式的地震记录想快速看一眼剖面长什么样结果打开专业软件要等几分钟加载换个工区又得重新配坐标、调增益、设色标。这种场景在野外验收、井位论证、教学演示里太常见了——你只是想确认数据有没有问题却被迫走一遍完整处理流程。地震数据的快速可视化核心诉求就一个把二进制道集在几秒内变成能看的图像不依赖重型软件不写冗长的配置文件。适合谁做地震资料处理的工程师、带学生看实际数据的导师、需要快速抽检数据质量的现场人员。下面这份资源包围绕“快速”二字展开用 Python 生态把 SEG-Y 读取、道头解析、剖面绘制串成一条最短路径同时保留参数调整的入口让你在快和可控之间找到平衡。2. 读懂 SEG-Y 的道头与数据体为什么不能直接 imshow2.1 地震记录的二进制结构决定了读取方式SEG-Y 文件不是一张图片它由 3200 字节的文本头、400 字节的二进制文件头以及若干道组成。每道又分 240 字节的道头和一个浮点数组。很多人第一次用numpy.fromfile直接读整个文件然后 reshape 成二维数组结果画出来的剖面全是噪声——因为道头混进了数据体。常见做法是先用segyio或obspy把道头和数据分离拿到data矩阵和header字典再决定怎么画。这里有个关键参数采样间隔dt和道间距dx它们决定了横纵坐标的物理单位。如果道头里没写你得从观测系统表里手动补否则画出来的图只有像素意义没有时间深度。import segyio import numpy as np # 打开 SEG-Y 文件只读模式不修改原始数据 with segyio.open(line_001.sgy, r, ignore_geometryTrue) as f: # 读取所有道的数据转成 (道数, 采样点数) 的二维数组 data segyio.tools.collect(f.trace[:]) # 提取采样间隔微秒用于计算时间轴 dt f.bin[segyio.BinField.Interval] # 提取道数 n_traces f.tracecount # 提取每道的 CDP 号用于横坐标 cdps [f.header[i][segyio.TraceField.CDP] for i in range(n_traces)] print(f数据体形状: {data.shape}, 采样间隔: {dt} us, 道数: {n_traces})这段代码的逻辑是先建立文件句柄用collect一次性把所有道读进内存。ignore_geometryTrue是为了跳过几何定义加快打开速度——快速可视化不需要完整的观测系统。dt从二进制头里取单位是微秒后面画图时转成毫秒。cdps列表用来做横轴标签如果道头里 CDP 号缺失可以退而用道序号。参数说明segyio.open的strict参数默认 True遇到不规范文件会报错快速查看时可以设strictFalse跳过校验但要注意可能读到错误数据。2.2 数据体转图像增益与色标是两道坎拿到data矩阵后直接plt.imshow(data)往往一片黑或者一片白。原因是地震振幅动态范围大深部弱信号被浅部强能量淹没。常见做法是加一个增益函数比如t^2增益或者 AGC自动增益控制。t^2增益简单每个采样点乘以时间的平方补偿球面扩散。AGC 更灵活用滑动窗口计算均方根振幅做归一化。色标方面地震剖面习惯用红白蓝或灰度对称色标要设vmin-vmax让零振幅对应白色。下面是一个带增益和色标控制的绘图函数。import matplotlib.pyplot as plt def plot_section(data, dt2000, gain_power2, clip_percentile98): data: 二维数组 (道数, 采样点数) dt: 采样间隔微秒 gain_power: 增益指数常用 1~3 clip_percentile: 裁剪百分位控制色标范围 n_traces, n_samples data.shape # 时间轴单位毫秒 time_ms np.arange(n_samples) * dt / 1000.0 # 构建 t^gain_power 增益曲线 gain_curve (time_ms / time_ms.max()) ** gain_power # 对每道应用增益 gained data * gain_curve[:, np.newaxis] # 计算对称色标范围 vmax np.percentile(np.abs(gained), clip_percentile) # 绘制剖面 plt.figure(figsize(10, 6)) plt.imshow(gained.T, aspectauto, cmapseismic, vmin-vmax, vmaxvmax, extent[0, n_traces, time_ms.max(), 0]) plt.xlabel(道号) plt.ylabel(时间 (ms)) plt.colorbar(label振幅) plt.title(地震剖面快速可视化) plt.show() # 调用 plot_section(data, dtdt, gain_power2, clip_percentile98)逻辑说明gain_curve是一个从 0 到 1 的数组乘以数据后浅层增益小、深层增益大补偿能量衰减。clip_percentile取绝对值第 98 百分位避免个别强振幅把色标拉爆。extent参数把时间轴翻转符合地震剖面深度向下的习惯。参数怎么改如果剖面还是太暗把gain_power调到 3如果噪声被放得太大降到 1。clip_percentile一般 95~99数值越小对比度越强但可能丢失弱反射。注意imshow的aspectauto必须加否则道数和采样点数比例失调剖面会被压扁。3. 批量抽检与交互式浏览把单张图变成工作流3.1 多文件批量出图循环里加异常捕获实际工区往往有几十条测线一条条手动打开不现实。常见做法是写一个批处理脚本遍历目录下所有.sgy或.segy文件每读一个就生成一张 PNG 缩略图文件名带上测线号。这里有个坑不同文件的道头格式可能不一致有的用标准 SEG-Y有的用扩展文本头。所以循环里要加try...except遇到读不了的文件跳过并记录不要整个脚本崩掉。下面是一个批量出图的骨架。import os import glob def batch_plot(input_dir, output_dir, dt_default2000): os.makedirs(output_dir, exist_okTrue) sgy_files glob.glob(os.path.join(input_dir, *.sgy)) \ glob.glob(os.path.join(input_dir, *.segy)) for sgy_path in sgy_files: base_name os.path.splitext(os.path.basename(sgy_path))[0] try: with segyio.open(sgy_path, r, ignore_geometryTrue) as f: data segyio.tools.collect(f.trace[:]) dt f.bin[segyio.BinField.Interval] or dt_default # 调用前面的绘图函数保存为文件 plot_section(data, dtdt, gain_power2, clip_percentile98) plt.savefig(os.path.join(output_dir, f{base_name}.png), dpi120) plt.close() print(f[OK] {base_name}) except Exception as e: print(f[FAIL] {base_name}: {e}) # 调用 batch_plot(./sgy_data, ./quicklook)逻辑说明glob同时匹配两种扩展名避免漏文件。dt从文件头取取不到就用默认值 2000 微秒。plt.savefig后必须plt.close()否则内存里堆积几十张图会爆。dpi120是缩略图的合适分辨率再高文件太大。参数说明如果目录里文件很多可以加max_files参数限制数量先抽检前 20 条。异常捕获里打印文件名和错误信息方便后续排查是文件损坏还是格式不兼容。3.2 交互式浏览用 matplotlib 的滑块调增益静态图只能看个大概想细看某段同相轴还是得能实时调参数。matplotlib.widgets里的Slider可以挂在图上拖动时重绘。常见做法是把增益指数和色标百分位做成两个滑块范围分别设 0.5~4 和 90~100。这样不用改代码就能找到最佳显示效果。下面是一个最小交互示例。from matplotlib.widgets import Slider def interactive_section(data, dt2000): n_traces, n_samples data.shape time_ms np.arange(n_samples) * dt / 1000.0 fig, ax plt.subplots(figsize(10, 6)) plt.subplots_adjust(bottom0.25) # 初始增益和裁剪 init_gain 2.0 init_clip 98 gain_curve (time_ms / time_ms.max()) ** init_gain gained data * gain_curve[:, np.newaxis] vmax np.percentile(np.abs(gained), init_clip) img ax.imshow(gained.T, aspectauto, cmapseismic, vmin-vmax, vmaxvmax, extent[0, n_traces, time_ms.max(), 0]) ax.set_xlabel(道号) ax.set_ylabel(时间 (ms)) # 增益滑块 ax_gain plt.axes([0.15, 0.1, 0.65, 0.03]) slider_gain Slider(ax_gain, 增益指数, 0.5, 4.0, valinitinit_gain) # 裁剪滑块 ax_clip plt.axes([0.15, 0.05, 0.65, 0.03]) slider_clip Slider(ax_clip, 裁剪百分位, 90, 100, valinitinit_clip) def update(val): g slider_gain.val c slider_clip.val gc (time_ms / time_ms.max()) ** g gd data * gc[:, np.newaxis] vm np.percentile(np.abs(gd), c) img.set_data(gd.T) img.set_clim(-vm, vm) fig.canvas.draw_idle() slider_gain.on_changed(update) slider_clip.on_changed(update) plt.show() # 调用 interactive_section(data, dtdt)逻辑说明update函数在滑块变化时重新计算增益和色标范围只更新图像数据不重建画布所以响应很快。set_clim控制颜色映射范围draw_idle触发重绘。参数说明滑块范围根据经验设增益指数低于 0.5 几乎没效果高于 4 噪声会盖过信号。裁剪百分位低于 90 会把有效反射裁掉高于 100 没意义。注意plt.subplots_adjust(bottom0.25)给滑块留出空间否则滑块会盖住横轴标签。4. 避坑与排查快速可视化里最容易翻车的五件事4.1 现象剖面上下颠倒浅层跑到下面去了原因imshow默认原点在左上角而地震数据的时间轴是向下增加的。如果extent参数写成[0, n_traces, 0, time_ms.max()]时间轴就是从上到下递增浅层反而在底部。解决把extent的第三、四个参数对调写成[0, n_traces, time_ms.max(), 0]让 y 轴反向。或者用ax.invert_yaxis()但extent方式更直接。4.2 现象读大文件时内存爆了程序被系统杀掉原因segyio.tools.collect会把所有道一次性读进内存。一条 10GB 的 SEG-Y 文件读进来就是 10GB 的float32数组加上增益计算时的临时数组内存翻倍。解决分块读取。用f.trace[i:ichunk]逐块处理每块画完就释放。或者用segyio.tools.cube只读需要的道范围。快速可视化不需要全精度可以隔道抽样比如每 5 道取 1 道数据量直接降到 20%。4.3 现象色标范围被个别异常振幅主导剖面一片灰原因数据里有个别坏道或者尖脉冲振幅是正常值的几百倍。np.percentile取 98% 时如果坏点占比超过 2%vmax会被拉得很大正常信号全被压成灰色。解决先做异常值剔除。用np.median和np.std算每道的稳健统计量把超过 5 倍标准差的采样点置零。或者把clip_percentile降到 95牺牲一点弱信号对比度换取整体可读性。4.4 现象道头里的 CDP 号全是 0横轴没有物理意义原因有些处理系统输出的 SEG-Y 文件不写 CDP 号或者写在非标准字节位置。segyio按标准道头读读出来就是 0。解决先打印道头字典看看哪些字段有值常见的有TraceField.TRACE_SEQUENCE_LINE或TraceField.FieldRecord。如果都没有就用道序号代替并在图注里说明“横轴为道序号非 CDP”。不要硬编一个假 CDP 号否则后续对比会出错。4.5 现象批量出图时中文标题变成方框原因matplotlib默认字体不支持中文plt.title(地震剖面)会显示成豆腐块。解决在脚本开头设置plt.rcParams[font.sans-serif] [SimHei]和plt.rcParams[axes.unicode_minus] False。如果系统没有 SimHei换成Microsoft YaHei或WenQuanYi Micro Hei。注意axes.unicode_minus要设 False否则负号也会变方框。5. 进阶技巧把可视化嵌进处理流程做质量监控5.1 用子图矩阵对比不同增益参数的效果调增益时经常纠结选哪个指数。与其反复改代码不如一次性画出 2x2 子图分别用 1、2、3、4 四个增益指数直观对比。下面是一个快速对比函数。def compare_gains(data, dt2000, gains(1, 2, 3, 4)): n_traces, n_samples data.shape time_ms np.arange(n_samples) * dt / 1000.0 fig, axes plt.subplots(2, 2, figsize(12, 8)) for ax, g in zip(axes.flat, gains): gc (time_ms / time_ms.max()) ** g gd data * gc[:, np.newaxis] vm np.percentile(np.abs(gd), 98) ax.imshow(gd.T, aspectauto, cmapseismic, vmin-vm, vmaxvm, extent[0, n_traces, time_ms.max(), 0]) ax.set_title(f增益指数 {g}) ax.set_xlabel(道号) ax.set_ylabel(时间 (ms)) plt.tight_layout() plt.show() # 调用 compare_gains(data, dtdt)逻辑说明zip(axes.flat, gains)把四个子图和四个增益值配对循环里复用同一套绘图逻辑。tight_layout自动调整间距避免标签重叠。参数说明gains元组可以按需改成(0.5, 1, 2, 3)覆盖更宽的范围。如果数据量大这个函数会开四个图像缓冲内存吃紧时改成逐个画、逐个保存。5.2 把缩略图拼成工区概览图一个工区几十条测线每条的缩略图单独看太碎。常见做法是用PIL或matplotlib把缩略图拼成网格加上测线号标签一张图看全工区。下面是用PIL拼接的示例。from PIL import Image import os def mosaic_thumbnails(thumb_dir, output_path, cols5): files sorted([f for f in os.listdir(thumb_dir) if f.endswith(.png)]) if not files: print(没有找到缩略图) return # 读取第一张图获取尺寸 first Image.open(os.path.join(thumb_dir, files[0])) w, h first.size rows (len(files) cols - 1) // cols # 创建空白画布 mosaic Image.new(RGB, (cols * w, rows * h), white) for idx, fname in enumerate(files): img Image.open(os.path.join(thumb_dir, fname)) r, c divmod(idx, cols) mosaic.paste(img, (c * w, r * h)) mosaic.save(output_path) print(f工区概览已保存: {output_path}) # 调用 mosaic_thumbnails(./quicklook, ./overview.png, cols5)逻辑说明divmod同时算出所在行和列paste把每张缩略图贴到对应位置。cols5表示每行 5 张根据缩略图数量和屏幕宽度调整。参数说明如果缩略图尺寸不一致paste会按原尺寸贴导致错位。所以批量出图时要固定figsize和dpi保证每张缩略图尺寸相同。Image.new的背景色设白色和地震剖面的底色一致拼出来更整洁。5.3 验证可视化结果是否忠实于原始数据快速可视化最怕“好看但失真”。验证方法很简单选一道数据把增益后的振幅曲线和原始振幅曲线画在一起看增益是否只改变了幅度比例没有引入虚假同相轴。另外用np.allclose检查增益计算前后的数据形状和道数是否一致。我一般会在批量出图后随机抽 3 条测线用segyio重新读原始文件对比data矩阵的均方根振幅误差超过 1% 就说明读取环节有问题。从那以后我每次做快速可视化都强制走一遍“读一道、画一道、对一道”的流程确认无误再批量跑。希望帮到你。本文还有配套的精品资源点击获取