
简介本资源是一份面向地球物理、地质勘探及地震信号处理初学者与科研人员的MATLAB轻量级工具脚本聚焦地震记录中波形可视化与振幅分析核心需求。压缩包仅含1个关键文件——wigb.m体积仅1KB为纯MATLAB函数脚本可直接加载SEED/ASCII格式地震数据实现Wiggle Traces with BackgroundWIGB图绘制在平滑背景上叠加抖动波形直观凸显P波、S波等特征振幅变化。脚本封装了数据读取、时域滤波、绝对振幅计算及二维振幅分布绘图等关键流程调用plot、imagesc等基础绘图命令完成专业级地震图件生成。已有275人学习下载适用于高校地球科学实验教学、科研项目快速原型验证及地震信号处理入门实践无需额外依赖工具箱开箱即用便于理解WIGB图原理并拓展自定义分析逻辑。 做地震数据处理的大概没有几个人能绕开wigb。这个来自加拿大卡尔加里大学 CREWES 实验室 MATLAB 工具箱的函数几乎成了地震剖面的“标准画法”。你给它一个二维振幅矩阵它给你画出变面积波形图——反射同相轴连不连续、振幅强不强、断层怎么就错了位一眼就能读出来。这篇文章就是围绕wigb展开的我会讲清楚它的原理、参数怎么设、实际绘图流程以及我这么多年里踩过的坑。适合刚开始接触地震记录可视化、或者想把手头数据画得更规范的人参考。wigb这个东西看着不起眼但如果你只是第一次拿到wigb.zip解压扔进 MATLAB 搜索路径就开始画大概率会画出一张让人摸不着头脑的图。因为它的输入输出、数据方向、缩放系数都有讲究一个地方没弄对图就是反的、糊的甚至什么都看不见。下面我把它拆开讲。1. 从“画一条曲线”到“画一个剖面”wigb 到底解决了什么问题1.1 wigb 的出身与定位wigb全名是 wiggle trace plot后面那个 b 代表 black 或者 banded也就是“变面积黑白波形显示”。它最早是作为 CREWES 研究工具箱的一部分发布的专门用来画地震道集和地震剖面。和 MATLAB 自带那些通用绘图函数不一样wigb是真正面向地震数据组织方式的数据按道存放、每道是一条时间序列横向是道号或者测线位置纵向是时间或深度。很多人拿到的wigb.zip实际上就是提取出来的单个wigb.m文件。这个文件虽然不大但想自己从零写一个同样效果的函数还真得费点功夫。它涉及波形抽稀、正半周闭合多边形、批量填充绘图这些细节。所以我的建议是能用现成的就用现成的但一定要搞清楚它内部干了什么不然出了问题你都不知道往哪查。1.2 变面积波形图到底好在哪先回答一个基础问题为什么地震剖面要用变面积波形图而不是直接用plot(data(:, 10))一类的裸曲线原因很简单——一张剖面上有几十上百道如果每道都画成细线远看就是一团乱麻但如果把正振幅部分涂黑负振幅留白那么反射强的位置自然形成又宽又黑的“带子”弱反射则变成细细的灰色痕迹。这种明暗相间的纹理恰好对应地质层位的反射特征解释人员可以像看照片一样快速识别地层界面。变面积显示的另一个好处是振幅信息可视化。波形图的横向偏移量本质就是振幅大小正半周填充后振幅越强填充面积越宽视觉上“黑度”就越高。这也是为什么很多处理报告、答辩 PPT 里都用这种显示方式——它能让振幅强弱直接形成图像对比且不依赖颜色映射。2. wigb 函数语法与参数逐项拆解2.1 调用格式与数据组织方式标准调用是wigb(data, scale, xcoord, tcoord)四个输入参数前两个基本是必须的后两个看情况data二维矩阵大小是nt × nx。行是时间采样点列是地震道。这是最核心的约定很多人在这一步就栽了跟头把矩阵存成了nx × nt画出来整个剖面横竖颠倒却还以为是函数的问题。scale振幅缩放系数默认是 1。它控制波形横向展宽的大小。xcoord每道对应的横向坐标向量缺省时是1:nx。如果你有实际道头里的 CDP 号、炮检距、测线桩号等都可以传进来。tcoord每个采样点对应的时间或深度向量缺省时是1:nt。注意它是纵轴坐标和data的行数必须严格相等。提示wigb不是在 MATLAB 基础工具箱里自带的使用前需要先把wigb.m文件所在目录添加到搜索路径。如果用的是完整的 CREWES 工具箱那直接调用即可。2.2 scale 参数画地震剖面最需要调的东西scale这个参数经常被忽略但其实最关键。它的物理含义可以理解成每个采样点的振幅值乘以scale之后相当于多少道间距。scale越大波形横向偏移越夸张看起来越“胖”scale太小波形挤在道中心附近几乎看不出变化。我在实际项目里的习惯是先做一次归一化data data ./ max(abs(data(:))); wigb(data, 1.0, 1:nx, t);这样振幅范围控制在[-1, 1]scale1时最大振幅正好对应大概一个道间距的偏移量显示效果比较均衡。如果归一化之后还觉得波形太瘦再慢慢加到1.2、1.5不要一开始就瞎调一个很大的数否则整张图就是一片纯黑毫无信息量。2.3 坐标向量与方向检查xcoord和tcoord是两个很容易搞错长度的参数最典型的报错是Vectors must be the same length或者画出来坐标对不上。我每次写代码前都会强制自己检查一遍矩阵维度[nt, nx] size(data); length(tcoord) nt % 必须为 true length(xcoord) nx % 必须为 true还有一个方向问题。地震剖面通常约定时间向下增加也就是说道剖面显示的时候纵轴要从上往下是t0, tdt, t2dt...。有些版本的wigb内部已经做了set(gca, YDir, reverse)有些则没有。保险起见我在调用wigb之后一定会手动加一句set(gca, YDir, reverse);如果画完之后发现深层反射跑到了图上面就是这里没设置对。3. 从合成记录到真实数据的完整绘制流程3.1 准备数据从雷克子波构建合成地震记录为了让大家能直接跑通流程我先用合成数据做一个完整示例。先写一个生成雷克子波的函数function w ricker_wavelet(t, fdom, delay) % 零相位雷克子波 % t : 时间向量 % fdom : 主频单位 Hz % delay : 子波延迟单位 s w (1 - 2*pi^2*fdom^2*(t-delay).^2) .* exp(-pi^2*fdom^2*(t-delay).^2); end接下来生成一个 50 道的合成地震剖面包含三个反射层振幅各不相同dt 0.001; % 1ms 采样率 t 0:dt:1; % 记录长度 1 秒 nt length(t); nx 50; % 50 道 % 先构造单道反射序列三个子波振幅递减 ref zeros(1, nt); ref ref 1.0 * ricker_wavelet(t, 30, 0.1); ref ref 0.6 * ricker_wavelet(t, 30, 0.3); ref ref 0.3 * ricker_wavelet(t, 30, 0.55); % 复制成多道再加上横向振幅变化模拟透镜体 data repmat(ref(:), 1, nx); taper exp(-((1:nx) - 25).^2 / 200); data data .* repmat(taper, nt, 1); % 画图 wigb(data, 1.0, 1:nx, t); xlabel(道号); ylabel(时间 (s)); set(gca, YDir, reverse); title(合成地震记录剖面);这段代码跑通之后你会看到三组比较明显的同相轴中间一组因为加了横向衰减左右两端振幅弱、中间强这种效果是imagesc很难直接看出层次来的但wigb能很直观地反映出来。3.2 振幅归一化与显示预处理真实地震记录的振幅范围非常夸张。有的道能量特别强有的道基本是死道如果不做归一化直接送进wigb结果往往是那个强能量道把整张图横向撑开其余道全被压成一条细线。我的建议是显示之前先做两个处理第一按道做或者按全数据做振幅归一化。比如全数据归一化data_norm data / max(abs(data(:))); wigb(data_norm, 1.0, ...);如果某些道能量差异太大可以按道归一化env max(abs(data), [], 1); data_norm data ./ env; % 每道单独归一化 wigb(data_norm, 0.8, ...);但要注意按道归一化会抹掉道间振幅差异看相对振幅用全数据归一化更合适看同相轴形态用按道归一化更合适。具体选哪种取决于你的目的。第二如果数据里有明显的直流成分或低频漂移先detrend或者做一次带通滤波再显示。因为wigb对基线偏移很敏感基线只要偏一点正半周填充区域就会整体异常看起来整道都是灰的。3.3 坐标轴、图例与图形修饰wigb画完之后很多人就直接截图完事但作为要放进报告或者论文里的图还有几个修饰动作必不可少。纵轴如果是时间单位一定要带上。横轴如果是道数也要说明是道号还是 CDP 号。更推荐的做法是把实际坐标传进去比如cdp header(:, 1); % 从道头读出 CDP 号 time_axis (0:nt-1) * dt; % 实际时间轴 wigb(data, 1.0, cdp, time_axis); xlabel(CDP号); ylabel(时间 (s));字体的调整也不可忽略。wigb内部用的图形句柄操作较多有时会干扰后续的set(gca, ...)。但一般情况下在wigb之后设置坐标轴属性和字体是没问题的set(gca, FontName, Times New Roman, FontSize, 10);另外wigb画出来的图形对象是patch而不是普通线条这意味着你可以在顶部自由叠加解释线、井位、断层标记这些对象会保持在波形上面不会相互遮挡错乱。这个特性非常实用我在第 6 部分还会专门说。3.4 真实数据与大矩阵的绘图策略真实工业场景下一个三维地震数据体可能非常大。以二维测线为例常见规模是时间采样点数nt1000~3000道数nx500~2000。直接用wigb画 2000 道数据MATLAB 会创建上万个patch对象绘制速度会明显变慢甚至会卡死。我遇到这种情况一般是先抽稀再画。比如每两道取一道或者每四道取一道idx 1:4:nx; wigb(data(:, idx), 1.0, xcoord(idx), t);抽稀之后道间距变小波形密度自然增加显示效果通常反而更清晰。如果你主要目的是快速浏览那直接imagesc是最快的imagesc(xcoord, t, data); set(gca, YDir, reverse); colormap(gray); xlabel(道号); ylabel(时间 (s));我自己的工作流是先用imagesc全局看一遍确定感兴趣的时窗和道范围再截取局部数据交给wigb做精细出图。这样既省时间又能保住wigb的专业效果。4. 结合 SEG-Y 数据画实际剖面4.1 读取与组织地震数据现实中地震数据大多存在 SEG-Y 文件里。MATLAB 读取 SEG-Y 有现成工具箱比如 CREWES 的ReadSEGY、SEG-Y 官方发行的SegyMAT或者 MATLAB 自带fopen加底层读二进制。不同读取工具返回的数据格式略有差异但最终组织方式是一样的得到一个nt × nx的振幅矩阵外加道头信息。以比较常见的读取方式为例读出来之后通常是data segy.data; % nt × nx dt segy.dt; % 采样间隔 t (0:size(data,1)-1) * dt; cdp segy.cdp; % 每道的 CDP 号然后就是把data直接交给wigb。注意 SEG-Y 文件里的道有可能是按 CDP 排序也可能按炮集排列先确认顺序再画不然剖面会出现道顺序错乱的问题。4.2 快速浏览与精细出图结合对一个完整的二维测线我一般分三步走第一步全局浏览。用imagesc全测线扫一遍判断数据整体质量、有没有坏道、时窗范围同时锁定我们要展示的目标区域。第二步局部抽稀。在目标区域附近每隔几道抽样用wigb做一次预览调整scale让波形幅度看起来舒服。第三步最终出图。使用完整目标区域数据通常抽到不超过 200 道加上精细的坐标轴、字体、叠加标记输出高分辨率图片。这三步看起来简单但能帮你省下大量时间尤其是当数据规模较大的时候。我见过不少同事直接一把wigb画全测线结果等了十几分钟出一张没法用的黑图这就是没做抽稀和归一化。5. 常见问题与排查技巧实录5.1 图里一条线都没有或者整张图全黑全黑是wigb使用中最常见的问题原因基本是scale设置过大或者数据没有归一化。如果数据的振幅量级是10^4scale1意味着波形横向偏移达到上万道间距那整个图当然就是一片黑色。解决办法先做一次全数据归一化data data / max(abs(data(:)))然后从scale0.5开始一点点试。如果归一化后波形太瘦看不出同相轴那就加大 scale如果黑成一片就减小。这个参数本质上是“以道间距为单位的振幅展宽系数”记住这一点就很好调了。5.2 时间方向颠倒与坐标错位地震剖面时间向下这是行业惯例。但wigb有的版本内部默认YDir为normal也就是时间向上很多新手一画完就觉得“我的地层怎么反了”。另外如果tcoord向量是从大往小排也会导致类似问题。修正是统一的调用后固定加set(gca, YDir, reverse);坐标错位的问题多半是xcoord、tcoord长度与矩阵维度不匹配。比如xcoord长度不等于列数wigb内部会出错。这个没有捷径每次作图前写一个assert检查即可。5.3 道数过多画出来黑糊糊一片当横向道数很多比如超过 300 道变面积波形图很容易因为道间距太窄导致相邻道的黑色填充区域连在一起看起来不是剖面而是一块黑板。这其实不是振幅参数的问题而是道距太密了。解决办法有两个方向一是增大图幅宽度用set(gcf, Position, [100 100 1600 800])或者让坐标轴平铺占满更多空间二是抽稀显示比如只画奇数道或者三分之二的道。出论文图的时候我通常控制在 120~200 道这样既保留波形细节也能在印刷尺寸下清楚显示同相轴。5.4 绘制太慢与导出文件巨大wigb本质是逐道执行patch绘制对象数量大时性能会明显下降。如果只是浏览我强烈建议先用imagesc。如果需要wigb出图可以先把数据裁剪到目标时窗和道范围再做抽稀然后调用。另外要注意导出矢量图时的文件大小问题。wigb产生的patch对象在矢量图里会保留所有折线细节导出的 EPS、PDF、SVG 可能非常大有的甚至几十上百 MB插入论文会自动造成编译慢。我的经验是PPT 演示用 PNG 300dpi 就够正式论文里如果对矢量图要求很高可以把时窗截短一点再导出别一口气导整个长剖面。5.5 常见问题速查表现象可能原因处理方式整张图全黑scale 过大或数据振幅量级过大归一化后从小到大试 scale波形几乎看不见scale 过小增大 scale 到 1~2深层显示在图上方YDir 未设置加set(gca, YDir, reverse)坐标轴长度报错x/t 向量和矩阵维度不匹配size检查assert拦截道密后糊成一片道距太窄抽稀或拉宽图幅绘制极慢patch 对象太多抽稀、裁剪时窗、换 imagesc背景有灰蒙蒙的基线偏置数据含直流分量先 detrend 或滤波再显示6. 几个我平时最常用的实操经验6.1 wigb 之后叠加解释线wigb画出来的是带着patch对象的地震剖面完全可以用hold on继续往上添加内容。这个特性在做层位拾取、断层解释时非常有用。wigb(data, 1.0, x, t); hold on; plot(x, horizon_time, r-, LineWidth, 1.5); % 叠地层位线 scatter(x(10:10:nx), well_pos, 20, b, filled); % 标井位要注意的是因为纵轴通常已经设置了YDirreverse用plot叠加时坐标是自动跟随的不需要额外处理。这样一张带解释结果的地震剖面就出来了。6.2 颜色与背景的调整技巧经典 wigb 是白底黑波形。但有时候深色背景更能突出强反射。虽然 wigb 不接收颜色参数但可以在画完后通过设置gca的Color来改背景色或者用colormap影响后续叠加的正半周颜色取决于版本。如果你发现 wigb 画出来的填充是灰色的试着在调用前设置colormap(gray)或者colormap(flipud(gray))通常能恢复成黑正白负的标准样式。我自己的习惯是保持白底黑波形这样打印或用 contrast 方式看最清楚。想强调强振幅时把 scale 稍微调大一点点而不是靠改颜色因为颜色在这种图里反而会引入视觉干扰。6.3 把 wigb 封装成自己的绘图函数最后分享一个工程化技巧wigb参数虽然不多但每次都要写归一化、坐标轴、字体、标题还是太绕。我会把它封装成一个自己的绘图函数把常用设置一次性处理掉。function plot_seismic_profile(data, t, x, sc, title_str, cdp_label) % 快速绘制统一风格的地震剖面 if nargin 4 || isempty(sc), sc 1.0; end if nargin 5, title_str ; end if nargin 6, cdp_label 道号; end data data ./ max(abs(data(:))); wigb(data, sc, x, t); set(gca, YDir, reverse, FontName, Times New Roman, FontSize, 10); xlabel(cdp_label); ylabel(时间 (s)); title(title_str); end这种做法能保证同一个处理流程里所有出图风格一致也方便团队其他人复用。你完全可以在wigb基础上继续扩展比如支持按道归一化、自动截取时窗、自动保存 PNG 等功能。最后再碎碎念几句我在实际项目里用 wigb 的频率非常高但也越用越知道它的边界在哪。它适合出成果图、展示同相轴、突出振幅差异却不是一个适合做快速浏览或者大数据交互查看的工具。很多人一上来就指着 wigb 说“这软件怎么这么慢”其实是把它的适用场景理解错了。正确的做法是把它放在出图的最后一步前面用 imagesc 和各类质量监控手段把数据看清楚、调好最后再交给 wigb 出一张漂亮的剖面。另外也真心建议有条件的话把 wigb.m 的源码读懂一遍不长但里面用到的逐道 patch 绘制逻辑、振幅偏移计算、归一化方法都是很有启发的地震可视化思路。读懂之后你再调参数、再改造成自己的工具会有底气得多的。本文还有配套的精品资源点击获取