ARTICLE DETAIL

资讯详情

深耕网站视觉设计与运营推广的一线实战洞察。

MATLAB wigb.m 地震信号 WIGGLE 剖面显示与避坑指南

MATLAB wigb.m 地震信号 WIGGLE 剖面显示与避坑指南 简介这份资源面向地震学与地球物理方向的学习者和研究人员提供基于MATLAB的地震信号剖面可视化工具核心解决将地震记录转换为WIGGLE波形图以便观察P波、S波到达时间及振幅相位特征的问题。压缩包内共1个文件为单个.m脚本整体约1KB属于轻量级代码资源可直接在MATLAB环境中读取地震数据并生成剖面图适合作为信号处理与可视化练习的入门参考。目前已有331人学习下载说明该工具在相关教学与科研场景中具有一定实用价值。通过运行和理解脚本中的算法逻辑读者可以掌握地震数据预处理、动态范围调整与波形绘制的基本思路进而用于分析地震波传播特性、震源性质及地壳结构特征为地震活动模式研究和灾害风险评估提供直观的图形依据。1. 地震信号剖面生成从 wigb.m 看 WIGGLE 显示为什么值得单独封装如果你处理过主动源地震记录大概率遇到过这种场景采集回来的 SAC 或 SEG-Y 数据想快速看一眼同相轴连续性、初至是否拾取干净、面波有没有压制到位结果用 imagesc 画出来一片糊用 plot 逐道画又慢又乱。wigb 这个函数就是解决这个问题的——它把多道地震信号按 WIGGLE 方式叠成剖面图正负振幅用填充色区分道间距可调显示增益可控。压缩包里只有一个 wigb.mMATLAB 脚本不依赖工具箱拿来就能用。适合做近地表折射、反射数据处理的人也适合做被动源噪声成像时快速检查互相关函数的对称性。它不负责数据处理只负责把数据画对但恰恰是这一步很多人翻车。2. wigb.m 的输入输出约定读懂参数才能改对图2.1 函数签名与数据维度要求wigb.m 的典型调用形式是wigb(data, scale, x, z)其中 data 是二维矩阵行对应时间采样点列对应地震道。这个维度约定和 MATLAB 里 imagesc 的默认行为一致但和某些地震处理软件“道在行、时间在列”的习惯相反。如果你从 SEG-Y 读出来的数据是 nt × ntr直接传进去没问题如果是 ntr × nt必须先转置否则画出来的剖面时间轴和道轴会互换同相轴变成竖条纹。scale 参数控制波形幅度相对于道间距的比例。设道间距为 1scale 取 0.5 表示最大振幅占半个道间距取 1.0 表示波形可以碰到相邻道。实际用的时候如果数据动态范围大强振幅道会盖住弱道这时候要么先做 AGC要么把 scale 压到 0.3 以下。x 和 z 是可选的道号和时间轴向量不传就默认用 1:N 和 1:nt。% 假设 data 是 nt × ntr 的地震数据矩阵 nt size(data, 1); ntr size(data, 2); dt 0.002; % 采样间隔 2ms x 1:ntr; % 道号 z (0:nt-1) * dt; % 时间轴单位秒 % 基本调用scale 取 0.8 figure; wigb(data, 0.8, x, z); xlabel(Trace Number); ylabel(Time (s)); title(WIGGLE Display);这段代码里z 向量决定了纵轴的物理意义。如果你传的是采样点序号而不是时间纵轴标签就得手动改成“Sample Index”否则读图的人会误以为单位是秒。x 向量一般用道号就行但如果是变偏移距数据也可以传偏移距值这样横轴直接反映空间位置。scale 的选取没有公式我一般先试 0.6看强振幅有没有溢出再微调。2.2 填充逻辑与颜色映射WIGGLE 显示的核心是正振幅填充黑色、负振幅填充白色或反过来中间用一条竖线表示零线。wigb.m 内部通常用fill或patch实现对每一道循环把大于零的部分和小于零的部分分别填充。这里有个细节如果数据里有 NaNfill 会报错或者画出断裂的色块。常见做法是在调用前把 NaN 替换成 0或者用data(isnan(data)) 0处理掉。另一个坑是道数太多时循环填充会慢。比如 1000 道以上MATLAB 的 fill 循环可能跑几十秒。这时候可以改用imagesc加colormap(gray)做近似但会丢失 WIGGLE 的波形轮廓感。折中方案是每隔几道抽一道画 WIGGLE或者把数据分块显示。我一般超过 500 道就考虑抽道除非必须看全部道的连续性。% 处理 NaN 并限制显示道数 data(isnan(data)) 0; max_traces 500; if size(data, 2) max_traces idx round(linspace(1, size(data, 2), max_traces)); data data(:, idx); x x(idx); end这段逻辑说明NaN 替换是必须的否则 fill 会出问题抽道用 linspace 均匀取避免只显示前 500 道导致空间采样不均。抽道后 x 向量也要同步截取不然横轴标签和实际道号对不上。2.3 增益控制与道间距调整scale 参数本质上是一个全局增益但地震数据里浅层强、深层弱是常态。如果只想看深层弱信号全局 scale 调大会让浅层溢出。这时候有两个选择一是分时窗做 AGC二是用 wigb 的变体允许传入每道的 scale 向量。wigb.m 原版不一定支持向量 scale但改起来不难——把标量乘法改成逐道乘法即可。道间距的调整更直接x 向量如果传的是实际物理距离比如偏移距道间距就自动反映了采集几何。如果传的是道号道间距就是 1。想让剖面看起来更紧凑可以把 x 乘以一个小于 1 的系数但这样横轴刻度就失真了。我一般不动 x而是调 figure 的宽高比用axis tight让显示区域贴合数据范围。提示如果画出来的剖面上下颠倒检查 z 向量是不是递减的。wigb 默认按 z 的顺序画z 递减会导致时间轴反向。3. 从 SAC 到 WIGGLE完整数据准备与绘图流程3.1 读取 SAC 数据并组装矩阵SAC 格式在地震学里很常见MATLAB 里可以用readsac或rdsac读取。假设你有一组 SAC 文件按道号命名需要先读进内存再拼成矩阵。这里的关键是确保所有道的采样率和时间窗一致否则拼出来的矩阵行数不匹配。% 读取一组 SAC 文件并组装成数据矩阵 files dir(*.sac); ntr length(files); data []; for i 1:ntr [t, amp, hdr] readsac(files(i).name); if isempty(data) nt length(amp); data zeros(nt, ntr); dt hdr.delta; end if length(amp) ~ nt error(采样点数不一致%s, files(i).name); end data(:, i) amp; end z (0:nt-1) * dt; x 1:ntr;这段代码先读第一个文件确定 nt 和 dt然后预分配 data 矩阵。循环里检查每个文件的采样点数不一致就报错。实际项目中SAC 头里的delta可能因为浮点精度有微小差异建议用第一个文件的 dt 统一或者用round对齐。如果道数多readsac循环读可能慢可以改用readsegy批量读 SEG-Y但 SAC 没有批量接口只能循环。3.2 预处理去均值、带通滤波与 AGC原始 SAC 数据通常有直流偏移直接画 WIGGLE 会看到零线偏移。去均值是第一步。带通滤波看需求如果只想看波形连续性不滤波也行如果要突出特定频段用 Butterworth 滤波器。AGC 是可选项但做 WIGGLE 显示时AGC 能让深浅层振幅更均衡。% 去均值 data data - mean(data, 1); % 带通滤波 1-20 Hz fs 1/dt; [b, a] butter(4, [1 20]/(fs/2), bandpass); data filtfilt(b, a, data); % AGC窗口 0.5 秒 win round(0.5 / dt); for i 1:ntr env movmax(abs(data(:, i)), win); env(env 0) 1; data(:, i) data(:, i) ./ env; end去均值用mean(data, 1)对每道单独做避免道间均值污染。滤波用filtfilt实现零相位不会引入时移。AGC 用滑动最大值做包络窗口 0.5 秒适合看区域地震信号如果是反射数据窗口可以缩到 0.1 秒。注意movmax在旧版 MATLAB 里可能没有可以用循环或filter替代。3.3 调用 wigb 并调整显示参数数据准备好后调用 wigb 就是一行的事但显示参数需要根据数据特点调。比如折射数据初至强scale 可以取 0.5反射数据弱scale 取 1.0 甚至 1.5。如果道数少比如 48 道道间距可以设大一点让波形不重叠。figure(Position, [100 100 800 600]); wigb(data, 0.7, x, z); colormap(gray); xlabel(Trace Number); ylabel(Time (s)); title(WIGGLE Profile); set(gca, YDir, reverse); % 时间轴向下 axis tight;set(gca, YDir, reverse)是地震剖面的标准做法时间向下增加。wigb 内部可能已经处理了但手动设一次更保险。axis tight去掉多余白边。如果剖面看起来太挤可以调 figure 的 Position 宽度或者用pbaspect设宽高比。注意如果 wigb.m 内部用了axis image或类似命令可能会覆盖你的 YDir 设置。调用后检查一下纵轴方向不对就再设一次。4. 避坑与排查wigb 用错场景的五个血泪教训4.1 现象剖面显示为空白或只有零线原因数据矩阵全零或者 scale 太小导致波形缩成一个点。常见于读 SAC 时头文件里的scale因子没乘或者滤波后数据被归一化到极小值。解决先max(abs(data(:)))看数据范围如果是 1e-10 量级说明没做单位换算。SAC 头里的scale因子要乘到 amp 上。scale 参数至少取数据最大值的 0.5 倍道间距。4.2 现象同相轴断裂或错位原因道间采样率不一致或者拼矩阵时某些道少采样点用零填充导致相位跳变。也可能是 SAC 头里的b值不同时间零点没对齐。解决读数据时检查每道的delta和npts不一致的道要么重采样要么剔除。时间零点用b值对齐t b (0:nt-1)*dt。4.3 现象填充色块溢出到相邻道原因scale 太大强振幅波形超过了道间距的一半fill 时覆盖了相邻道的区域。这在浅层折射或面波干扰严重时常见。解决降低 scale 到 0.3-0.5或者先做 AGC。如果只想看弱信号可以对数据做时变增益浅层压制、深层放大。4.4 现象绘图速度极慢几百道要等几分钟原因wigb.m 用循环 fill每道一次 patch道数多时 MATLAB 图形渲染开销大。解决抽道显示或者改用imagesc加colormap(gray)做快速预览。如果必须看全部道把数据分块每块 200 道画多个子图。4.5 现象保存为 PDF 或 PNG 后波形变模糊原因MATLAB 默认的渲染器是 paintersfill 对象在矢量输出时可能被简化。或者分辨率设太低。解决保存时用print -dpdf -r300指定 300 dpi或者用exportgraphics函数。如果还是模糊把 figure 的 Renderer 设为 opengl但矢量输出会变成位图。5. 进阶技巧用 wigb 做被动源噪声成像的互相关函数质量检查被动源噪声成像里互相关函数NCF的对称性和信噪比直接决定频散提取的可靠性。我一般把台站对按间距排序把 NCF 按道号排成矩阵用 wigb 画出来。正负延迟对称的话WIGGLE 图里应该看到左右对称的同相轴如果只有一边有信号说明互相关没做对称化或者某方向噪声源占优。具体操作读 NCF 数据每个台站对一个文件时间轴从 -T 到 T。拼矩阵时按台站间距从小到大排列这样横轴就是间距纵轴是延迟时间。scale 取 0.5 左右因为 NCF 振幅通常归一化过。% NCF 矩阵组装与 wigb 显示 ncf_files dir(NCF_*.mat); n_pairs length(ncf_files); ncf_data []; dist []; for i 1:n_pairs load(ncf_files(i).name); % 包含变量 ncf 和 t if isempty(ncf_data) nt length(ncf); ncf_data zeros(nt, n_pairs); end ncf_data(:, i) ncf; dist(i) header.dist; % 台站间距 end [~, idx] sort(dist); ncf_data ncf_data(:, idx); dist dist(idx); figure; wigb(ncf_data, 0.5, dist, t); xlabel(Inter-station Distance (km)); ylabel(Lag Time (s)); title(NCF Symmetry Check);这段代码按间距排序后画图横轴是物理距离纵轴是延迟时间。如果看到对称的“V”形或“X”形同相轴说明 NCF 质量不错如果一边强一边弱就得回去检查互相关前的单台数据处理比如是否做了时域归一化和谱白化。这个检查步骤我每次做噪声成像都强制走一遍比看单条 NCF 波形高效得多。提示t 向量要包含负延迟从 -T 到 T且零延迟在中间。如果 t 是从 0 开始的画出来的图只有正半轴看不出对称性。从那以后我每次拿到新的地震数据不管多急都先用 wigb 画一张全道剖面确认同相轴没断、时间轴没反、振幅没溢出再往下做处理。这个习惯帮我省了至少三次返工。希望帮到你。本文还有配套的精品资源点击获取
返回列表