简介本资源是一份面向地球物理、地质勘探及地震信号处理初学者与科研人员的MATLAB轻量级工具脚本聚焦于地震记录波形可视化与振幅分析核心需求。压缩包仅含1个关键文件——wigb.m大小仅1KB为纯MATLAB函数脚本可直接加载SEED/ASCII格式地震数据实现Wiggle Traces with BackgroundWIGB图绘制在平滑背景上叠加振幅缩放的波形抖动线直观凸显P波、S波等特征相位及振幅变化趋势。脚本内嵌数据读取、时窗截取、振幅归一化及专业绘图逻辑无需额外工具箱即可运行适合教学演示、野外数据快速质检或科研预处理环节。目前已有275人学习下载读者可即刻获得一套结构清晰、注释明确、开箱即用的地震波形可视化解决方案显著降低MATLAB地震绘图入门门槛。1. wigb 到底是干什么的从地震剖面显示说起先说一个大家可能都经历过的场景。你在处理地震数据的时候手里是一炮野外采集回来的共炮点道集或者是做完速度分析之后的动校正道集又或者是偏移之后的地震剖面。数据本身是个二维矩阵行是时间采样点列是地震道。你想快速看看到底信噪比怎么样、同相轴有没有拉平、断层在哪个位置于是很自然地想到了 MATLAB 里最经典的波形显示函数wigb。标题里的 wigb.zip_matlab_wigb_地震wigb_地震记录_振幅 其实已经把关键信息都点出来了这是一个和 MATLAB 环境下绘制地震记录振幅剖面相关的函数包最核心的东西就是wigb。wigb是 wiggle plot 的缩写直译就是波形摆动图江湖上也叫变面积波形剖面。它的本职工作是把二维地震数据矩阵按道依次排列每一道用一条随振幅左右摆动的曲线绘制出来同时把曲线与零线围成的面积填充颜色从而同时呈现振幅强弱和波形形态。这种显示方式在地震勘探行业里用了快几十年至今仍是野外质量监控和成果展示的主力工具。我遇到很多刚接触地震数据处理的同学第一反应是我直接用plot画出来不就行了吗。用plot直观但是有个致命问题一炮数据几百上千道每道上千个采样点全画一起就是一坨黑线谁也看不清楚。wigb做的事情本质上是按道对齐、按方向展开、把波形姿态立体化把一道压成一条窄带状区域波形往左就是负振幅往右就是正振幅黑色填充部分是正的波形面积。这样几十上百道平铺开同相轴的延续性、振幅强弱的变化、断点的位置一眼就能看个大概。这篇文章就围绕wigb这个函数从函数本体、参数含义、实际调用、常见报错到进阶扩展完整讲一遍。适合正在学 MATLAB 地球物理编程的研究生、刚入职的物探处理工程师以及任何需要快速可视化地震记录的科研人员。即便你是完全的新手文末也有可以直接复制运行的完整脚本。2. 拿到的 wigb.zip 里装了些什么函数定位与配套依赖分析既然标题是wigb.zip你手上应该已经有一个压缩包里面大概率就是一个wigb.m文件也可能连带wigb2.m、seisplot.m之类的变体。先说清楚这个压缩包的来源和版本问题。2.1 wigb 函数的老家CWP 工具包在 MATLAB 地球物理圈子里wigb最能查到出处的地方是 Colorado School of Mines 的 Center for Wave PhenomenaCWP发布的 Seismic Unix / MATLAB 工具包。这个工具包是几十年来学术界和工业界通用的地震数据处理工具集里面包含大量经典的可视化函数。wigb在 CWP 文件里属于 SU 到 MATLAB 的桥接函数之一代码写得非常简洁但功力极深——它用最笨的办法双层循环逐道画图实现了最耐用的功能。如果你从网上或师兄师姐那里拿到的wigb.m打开后看到的函数声明可能是这样function wigb(seismic, t, x, scale, ornt)也可能是function [v, t, x] wigb(v, t, x, scale, ornt)这两种形式分别对应只画图和画图并返回坐标数据的用途。CWP 原版更接近第二种——它允许你传入数据矩阵v如果返回时带上变量v里可能存放了变换后的坐标方便你在图上叠加标注。不同的发行版本之间参数含义有细微差别这一点后面专门讲。2.2 wigb 的前后端依赖风险单独一个wigb.m文件其实是不完整的它对环境有一些隐性要求MATLAB 的fill函数用于填充正波形的面积所有 MATLAB 版本都支持默认调用时如果没传坐标向量t和x函数内部会按采样点数自动生成1:n的自变量但如果没有指定时间采样间隔横轴的物理含义就丢了绘图的 Figure 对象坐标轴方向需要调整通常地震剖面要求时间轴向下也就是纵轴要set(gca, YDir, reverse)否则默认的纵轴向上会把剖面颠倒过来。这些细节不是函数本身能解决的而是使用者的常识。我见过太多人把wigb跑通了但画出来的图上下颠倒以为程序坏了其实只是坐标轴方向的问题。2.3 如果压缩包里只有 wigb.m不要慌有些版本的压缩包还会带上plotimage.m、seisplot.m这些辅助函数。如果没有也不要紧——wigb本身不依赖它们。你只需要保证你的数据是一个二维矩阵行方向是时间/深度列方向是道号就可以直接开画。最基础的一次调用是这样figure; wigb(data);其中data是nt × nx的矩阵nt是每道的采样点数nx是总道数。如果跑出来是一条条波浪线恭喜你基础流程已经通了。3. 核心调用语法与参数逐项拆解从 t、x、scale 到 ornt这部分是 wigb 用得顺不顺手的关键。很多教程只是说这个函数可以画地震剖面但真到自己处理数据时不把参数吃透画出来的图永远差那么点意思。下面逐个拆解。3.1 基础调用形式wigb最完整的函数签名通常是wigb(v, t, x, scale, ornt)参数含义如下参数含义典型值/默认行为v地震数据矩阵nt × nx无默认值必传t时间/深度向量每道的纵坐标默认1:ntx距离/道号向量横坐标默认1:nxscale振幅缩放比例默认 1.0即每个道内自动归一化到最大振幅ornt剖面方向默认 plot即时间轴垂直向下道向水平第 1 个参数v不用多说就是你的数据本体。第 2 个参数t比较有意思。如果你传了行向量t长度必须等于nt如果你传的是列向量函数内部会做转置处理。这里容易出问题因为 MATLAB 对行列向量很敏感代码里如果没有写t t(:)之类的统一操作你传个size(t) [nt, 1]的列向量可能在第 4 维的计算上出错。实测中最稳的办法是先在外部把t统一成行向量传进去之前确认一下length(t) nt。第 3 个参数x同理长度必须等于nx。通常代表检波点桩号或 CDP 号直接传1:nx也能画但横坐标失去物理意义。第 4 个参数scale这是最实用也最容易被忽略的参数。它控制的是每个道内振幅的缩放方式。默认情况下wigb会遍历每一道找到该道的最大绝对值振幅然后把整道除以这个最大值再乘以scale。这样做的效果是不管是振幅为 0.001 的微震记录还是振幅为 1000 的强反射记录每道内部的波形幅度都会大致布满道宽的一半画出来的剖面不会因为某一道特别大或特别小而导致其他道看起来像一条直线。但这里有个副作用道的最大振幅被归一化后真正的振幅相对关系就丢失了。同一道内浅层强反射和深层弱反射的相对大小还在但是不同道之间的绝对能量差异被抹平了。做振幅保真显示时你要么把scale设成很小的数值避免截断要么自己在外部处理数据后再传入。第 5 个参数ornt控制剖面朝哪个方向画。原版wigb默认是时间轴垂直向下道轴水平向右。如果你想旋转 90 度变成测线剖面即道方向垂直向上传入ornt为其他值时函数内部会交换x和t的角色。这在对比不同测线剖面时挺有用。3.2 变体 wigb2同时绘制两组数据叠加显示CWP 工具包里还有一个非常实用的变体wigb2它允许你在同一坐标系下绘制两组地震数据适合做合成记录和实际地震记录的对比或者做动校正前后对比。用法类似wigb2(data1, data2, t1, t2, x, scale1, scale2)data1和data2的尺寸可以不一样但横坐标x必须一致。这里有个简化实现思路wigb2内部其实是先画第一组然后把第二组数据按道偏移叠加到第一组上。注意两组数据的时间采样率必须一致否则对比没有意义。这个函数在检查你的合成地震记录是否和原始剖面吻合这个场景下非常好用。3.3 t、x、scale 三者配合的经典控制技巧我个人的习惯是拿到数据后先计算这三个量dt 0.002; % 采样间隔单位秒 t (0:nt-1)*dt; % 时间向量从0开始 x 1:nx; % 道号 scale 1.0; % 振幅缩放然后调用figure; wigb(data, t, x, scale); set(gca, YDir, reverse); xlabel(Trace number); ylabel(Time (s)); title(Raw shot gather);跑完这一步你就得到一张最标准的地震剖面图了。加粗提示一下set(gca, YDir, reverse)这个操作不是可选项是必选项否则剖面在 MATLAB 默认坐标系下是上下颠倒的。4. 实测示例从合成记录到真实数据的完整可视化流程这一节我把从造数据到成品图的完整流程走一遍。你可以直接复制这段代码跑一跑体会 wigb 在真实工作流中的手感。4.1 构造一个简单的合成地震道集真实野外数据文件不可能人人都能随手拿到所以先用合成数据演示逻辑。假设我们要构造一个包含三个反射层的零偏地震道集每个反射层是一个 Ricker 子波到达时间分别在 0.2 秒、0.45 秒、0.7 秒处振幅分别为 1.0、0.8、0.5同时加入一点随机噪声。dt 0.002; t 0:dt:1.0; nt length(t); nx 60; data zeros(nt, nx); % 三个反射层 f0 30; % Ricker子波主频 w (1 - 2*(pi*f0*(t-0.2)).^2) .* exp(-(pi*f0*(t-0.2)).^2); data data w * linspace(1, 1, nx); w2 (1 - 2*(pi*f0*(t-0.45)).^2) .* exp(-(pi*f0*(t-0.45)).^2); data data 0.8 * w2 * linspace(1, 1, nx); w3 (1 - 2*(pi*f0*(t-0.7)).^2) .* exp(-(pi*f0*(t-0.7)).^2); data data 0.5 * w3 * linspace(1, 1, nx); % 加噪声 rng(42); data data 0.05 * randn(nt, nx);这里w * linspace(1, 1, nx)的作用是把一个行向量波形w扩展到每一道都大小相等。实际数据里每一道的波形当然不可能完全一样这里只是为了演示基础调用。4.2 用 wigb 显示并叠加解释标识画图并加注释figure(Color, w); wigb(data, t, 1:nx, 1.0); set(gca, YDir, reverse, FontSize, 10); xlabel(Trace number); ylabel(Time (s)); title(Synthetic shot gather with wigb); hold on; % 在0.45秒处画一条水平虚线表示解释层位 line([0 nx1], [0.45 0.45], Color, r, LineStyle, --, LineWidth, 1.5);跑完这段你会看到三条明显的波浪带加上红色的解释层位虚线。这个场景在工作中非常常见——处理解释的时候经常要在 wigb 剖面图上叠加层位标注、断层位置、井曲线等。4.3 对比 scale 参数对显示效果的影响这才是 wigb 最容易踩坑的地方。试着把scale从 1.0 调到 0.1 再调到 5figure; subplot(1,3,1); wigb(data, t, 1:nx, 0.1); set(gca, YDir, reverse); title(scale0.1); subplot(1,3,2); wigb(data, t, 1:nx, 1.0); set(gca, YDir, reverse); title(scale1.0); subplot(1,3,3); wigb(data, t, 1:nx, 5.0); set(gca, YDir, reverse); title(scale5.0);你会看到scale 0.1时波形的摆动范围很小噪声基本不可见振幅形态被压缩成细线适合看强反射同相轴的连续性scale 1.0时每道最大振幅被归一化到道宽的一半整体均衡噪声看得清楚常用于常规质量监控scale 5.0时波形剧烈摆动每道之间的波形互相交叠画面基本糊成一团只有在故意夸大微弱信号时才用。所以scale的实用建议是常规显示用 1.0快速巡检用 0.50.8精细解释强反射层位用 1.52.0 都可以但要随时切换对比不要一次定死。4.4 真实数据的格式适配经验拿到野外观测的 SEG-Y 文件时一般是用ReadSEGY函数或SegyMAT工具包读入读入后的数据矩阵维度可能是nt × nx或nt × nx × ns多分量记录。要显示单分量先切片data_2d data3D(:, :, 1); % 第一个分量再取一炮shot1 data_2d(:, 1:120); nx size(shot1, 2);然后正常调用wigb。实测下来2000 道 × 2000 采样点的大矩阵跑wigb会有点慢因为函数内部是逐道fill的。这时候建议先用decimate或resample降采样或者只画其中一部分道比如shot1(:, 1:10:end)先把整体轮廓看清再局部放大。5. 常见坑位与排查经验为什么你的 wigb 画出来是花的、反的、或者空白讲了这么多正面用法下面专门把我在实际使用中遇到过的坑集中列一遍。这些坑如果没踩过看代码几乎不会意识到但只要踩过一次就会长记性。5.1 坑一数据矩阵方向反了剖面看起来像噪声墙wigb内部处理时是把每一列当作一道来画的。如果你的数据矩阵是nx × nt也就是行是道号、列是采样点那么传给wigb后函数会把每一行当一道画出来的图上横轴实际是采样点而不是道号整个剖面形态会完全错乱本应该横向连续的同相轴变成了纵向的条纹纵向应该是随深度变化的时间波形全被压缩到了横轴上。看起来就是一团黑压压的噪声墙。排查方法非常简单先打印size(data)确认第一维是不是采样点数、第二维是不是道数。如果不是用transpose转置一下再传。5.2 坑二时间轴上下颠倒剖面头朝下这是最经典的问题。MATLAB 默认的坐标 y 轴方向向上时间轴的数值从 0 到 1画出来后 0 在底部、1 在顶部意味着浅层信息在屏幕下方——这和地震剖面的常规显示习惯完全相反。地球物理剖面的标准展示方式是浅层在上时间 0 在顶部、深层在下大时间在底部。处理办法就是我反复强调的那一句set(gca, YDir, reverse);如果画多个子图每个子图都要设置一次。另一个办法是在调用wigb前直接对数据取反即wigb(flipud(data))但这时纵轴标签就乱了不推荐。5.3 坑三个别道振幅特别大整条剖面被压扁野外数据经常出现个别道有突发强噪声比如工业干扰、感应雷等幅值可能是正常信号的几百倍。wigb默认的自动归一化是逐道计算的理论上不应该因为一道的问题影响其他道。但如果你用的是自定义的全局缩放即先data data / max(abs(data(:)))再传入那么极大值的那一道会把整条剖面的振幅都压得极其微小剖面上看起来除了那道亮点其他信息全都没了。解决办法是先用clip方式把异常道拉回合理范围data(data threshold) threshold; data(data -threshold) -threshold;或者直接用wigb的默认scale让它逐道归一化如果确实需要全局振幅保存建议把异常道剔除后显示。5.4 坑四空白剖面没有任何波形调用wigb(data)后图窗弹出但是一片空白或者只有坐标轴。这种情况绝大多数是数据本身传成了全零矩阵或 NaN 矩阵。fill函数遇到全 NaN 时会直接跳过。先检查数据是否包含 NaN用sum(isnan(data(:)))统计很多 SEG-Y 读取在道头缺失时会把数据区填成NaN。处理办法data(isnan(data)) 0;5.5 坑五速度太慢数据量大时转圈圈wigb内部是逐道填充2000 道 × 3000 采样点的大数据画一次可能要十几秒到几十秒。这时候可以通过三步优化先用downsample降采样时间方向比如从3000点降到1000点视觉上几乎无损抽稀显示部分道比如每隔 5 道显示一道作为快速预览显示核心区域比如先xlim到某 200 道再画这 200 道。另外wigb函数内部如果用的是patch而不是fill大矩阵下渲染效率会高一些。新版本 MATLAB 的图形引擎对fill也做了优化但仍是逐对象绘制无法与pcolor或imagesc这种栅格化绘制比速度。5.6 坑六填充颜色不对画出来是白底黑线或黑底白线wigb默认的填充色通常是蓝色或黑色填充正振幅部分。如果你打开代码会看到里面有一行fill(xp, yp, k)或fill(xp, yp, b)这就是颜色定义。有的版本还会用copper色系循环。如果你希望统一改成其他颜色比如把所有正振幅填充改为默认红色但保留黑色轮廓线可以直接修改wigb.m里面的fill行把颜色字符串改成[0.8 0 0]之类。改完记得保存到当前工作目录或路径上覆盖原函数。6. 进阶玩法在 wigb 基础上做批量交付图件的实用扩展到这里wigb 的核心用法已经讲透了。但实际生产中我们往往不是画一张图就完事而是需要对一个工区几百炮、几百条测线的地震数据进行批量可视化并按照统一格式输出。这节分享我平时项目中直接在 wigb 之上做二次扩展的三个实用方向。6.1 批量绘制多炮记录的自动排版脚本模板假设你有一个三维的野外采集数据矩阵data_all维度是nt × nx × nshot也就是每个时间采样点、每个接收道、每炮一个值。要批量把所有炮画到同一张多子图图件中nshot size(data_all, 3); figure(Color, w); for ishot 1:nshot subplot(ceil(sqrt(nshot)), ceil(sqrt(nshot)), ishot); wigb(data_all(:, :, ishot), t, 1:nx, 1.0); set(gca, YDir, reverse, FontSize, 8); xlabel(Trace); ylabel(Time (s)); title([Shot , num2str(ishot)], FontSize, 9); end注意这里的subplot网格如果太多每张子图会被压得很小波形细节看不出来。批量预览的目的本来就是快速查看整体质量所以推荐nshot比较大时分批次画每批 9 炮或 16 炮。6.2 叠加振幅曲线到 wigb 剖面上有些场景下我们不仅想看波形剖面还想把各道的均方根振幅叠加上去直观地看能量沿测线的变化。可以先计算每道的 RMS 振幅再用plot叠加rms_amp sqrt(mean(data.^2, 1)); figure; wigb(data, t, 1:nx, 1.0); set(gca, YDir, reverse); hold on; % 将RMS振幅归一化到右侧合适范围后叠加 rms_norm rms_amp / max(rms_amp) * 0.2 * nx; plot(1:nx, rms_norm * 0 0.5, r, LineWidth, 1.5);不过更常见的是把振幅曲线单独放在剖面下方用两个叠置的坐标轴实现。这里只是抛个思路实际项目里大家按需调整叠加位置。6.3 基于 wigb 的自定义剖面绘制函数因为 wigb 是一个成熟的开源函数而我们在具体项目里往往有一些统一格式要求所以很多工程团队会把它封装成自己的函数。我的习惯是写一个plot_seismic_customfunction plot_seismic_custom(data, t, x, title_str, scale, cmap_dir) figure(Color, w); wigb(data, t, x, scale); set(gca, YDir, reverse, FontSize, 12, LineWidth, 1.0); xlabel(CDP, FontSize, 12); ylabel(Time (s), FontSize, 12); title(title_str, FontSize, 13); grid on; colormap(cmap_dir); colorbar; end这里cmap_dir是可选的地图颜色映射参数画完wigb后手动设置colormap和colorbar不影响已绘制的波形但可以给后续叠加的 color 区域一些辅助配色。这个封装的好处是项目中所有成员都调用同一个函数出图格式自然统一。以后要改标题字体大小、加网格、换配色只改一个函数就行不用追着每个人的脚本去改。6.4 与机器学习的联动批量生成训练集预览图还有一个近年用得越来越多的场景做地震数据质量评估模型时需要人工标注大量剖面图像。wigb 画出来的剖面图可以作为自动标注工具的底图。你只需要对批量读取的数据调用 wigb 并导出高清图片for i 1:n h figure(Visible, off); wigb(data_all(:, :, i), t, 1:nx, 1.0); set(gca, YDir, reverse); print(h, [preview_, num2str(i, %03d), .png], -dpng, -r150); close(h); endVisible, off让图窗隐藏print直接输出图片-r150表示 150 DPI批量几千张图几分钟就能跑完。这里有个经验如果只需要预览质量DIP 设成 100 就够设太高会拖慢速度且占据大量磁盘空间。7. 要不要自己再造一个 wigb谈性能边界与替代方案既然 wigb 这么经典是不是所有场景都该无脑用它我坦率讲它并不是万能的。在某些特定需求下你可能需要考虑替代方案。7.1 wigb 的优点和不足优点显示效果专业波形变面积显示在地球物理领域内是默认标准调用简单一款代码复用几十年跨 MATLAB 版本兼容性极好填充色让同相轴强弱一目了然对解释工作帮助很大。不足数据显示量有限大矩阵时性能差没有内置的色标、道头标注、地震道编号标签等现代显示功能交互功能几乎为零不能手动缩放某个局部时自动更新波形填充基于fill的对象图片体积偏大导出 PDF 或 EPS 时经常几百 MB。7.2 替代方案对比功能需求推荐方案说明快速预览大数据剖面imagesc 自定义 colormap栅格化绘图速度极快但不显示波形摆动细节需要波形 色彩双显示pcolor或contourf比 wigb 快但视觉上更像热像图交互式拾取同相轴plot 鼠标响应可以手动点选层位功能自由但开发量大专业出版级剖面Seismic Unix 的suxwigb功能更强大但需要 linux 环境与 Python 混用matplotlib自定义 wiggle 绘制适合 python 技术栈团队7.3 我个人的选型建议如果项目时间允许、数据量不超过 1000 道就放心用 wigb它画出来的图件在答辩和报告中都拿得出手。但如果数据量巨大且要求出图速度建议先用imagesc做快速质检视图再用 wigb 对局部重点区域出正式图。至于自研一个 wigb 替代品我的观点是不必。除非你的项目有大量定制化交互需求否则在成熟工具上做二次封装比重新发明轮子高效得多。最后分享一个我自己工作时的习惯写任何地震数据可视化脚本时我都习惯把 wigb 和 imagesc 的代码放一起保存方便随机切换。因为同样的数据用 wigb 看的是波形细节用 imagesc 看的是能量分布两者互补缺一不可。刚开始接触 wigb 的朋友不用着急追求最花哨的效果先把基本调用跑通再对照本文的参数说明逐步调整你很快就能够画出专业图件。本文还有配套的精品资源点击获取