ARTICLE DETAIL

资讯详情

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

NMO动校正实战:从CMP道集到叠加剖面的Python实现与避坑指南

NMO动校正实战:从CMP道集到叠加剖面的Python实现与避坑指南 简介这份资源面向地震勘探数据处理方向的学习者与工程人员聚焦动校正NMO这一地震处理基础环节帮助理解如何消除因地下速度差异造成的时间错位使同一反射界面的地震波在时间轴上对齐从而提升剖面分辨率与对比度。压缩包为rar格式共4个文件包含3个m脚本与1份readme说明整体仅2KB属于轻量级代码示例便于快速阅读与移植到Matlab环境中运行调试。内容涉及速度模型估算、时间域校正及与叠前深度偏移等技术的配合思路适合作为动校正算法入门与代码实现的参考。目前已有450人学习下载读者可借此梳理NMO处理流程、对照脚本理解校正步骤并在此基础上开展地震资料处理实验与地质解释练习。1. NMO 动校正从 CMP 道集到叠加剖面的关键一步如果你手头有一份 CMP 道集反射同相轴弯得像一把弓远偏移距的反射时间明显大于近偏移距那说明动校正还没做。NMONormal Moveout正常时差动校正就是要把这些因为炮检距不同而产生的时间延迟拉平让同一反射界面的能量在叠加时能对齐。没有这一步后面叠加出来的剖面就是一团模糊分辨率根本谈不上。地震勘探里NMO 动校正属于常规处理流程中承上启下的环节——上面接着速度分析下面直接影响叠加和偏移成像的质量。做地震处理的人不管你是刚入行的处理员还是做了几年的老手NMO 都是绕不过去的基本功。这篇文章不讲教科书定义只讲怎么在代码里把它跑通、参数怎么调、哪些地方容易翻车。2. NMO 动校正的原理与速度场依赖关系2.1 时距曲线方程决定了校正量的计算方式NMO 动校正的核心公式并不复杂。对于水平层状介质反射波旅行时 t 与炮检距 x 的关系可以写成t² t₀² x²/v²_nmo其中 t₀ 是零炮检距双程旅行时v_nmo 是动校正速度。动校正量 Δt t - t₀。实际操作中我们拿到的是离散采样后的地震道数据每个样点对应一个时间值需要根据上述公式计算每个样点应该移动到哪个时间位置。这里有个容易被忽略的点动校正速度不是介质真实速度而是一个使同相轴拉平的等效速度。它通常大于均方根速度在存在各向异性或倾斜地层时偏差更大。很多新手直接把叠加速度当成层速度用结果校正后同相轴还是弯的这就是原因。从实现角度看NMO 校正有两种常见做法一种是直接对每个样点做时间映射把原始时间 t 上的振幅搬到 t₀ 位置另一种是在频率域做相位移动。时域方法直观、容易实现但需要处理拉伸畸变频域方法精度高但计算量大。实际生产中用得最多的还是时域映射加拉伸切除。2.2 速度分析精度直接决定校正效果NMO 校正对速度的敏感度非常高。速度偏大同相轴校正不足远道仍然向下弯速度偏小同相轴被过度校正远道反而向上翘。这两种情况在速度谱上都表现为能量团不聚焦。我一般会先做速度分析生成速度谱拾取一系列 t₀ 对应的 v_nmo 值形成速度函数。然后把这个速度函数应用到 CMP 道集上做 NMO 校正。校正完了之后再看道集是否拉平如果不平回去修改速度拾取迭代两到三次基本就能收敛。这里有个经验速度谱的等速度线间隔不要太粗一般 50 m/s 一档比较合适。太粗了拾取精度不够太细了计算量大且拾取效率低。另外速度分析之前最好先做一下道集均衡或者振幅补偿不然浅层能量太强会压制深层弱反射速度谱上深层的能量团根本看不出来。2.3 用 Python 实现一个最小可跑的 NMO 校正下面这段代码演示了如何对一个合成的 CMP 道集做 NMO 校正。输入是一个二维数组行是时间采样点列是炮检距。速度用一个简单的常速实际生产中应该从速度分析结果里读取。import numpy as np def nmo_correction(cmp_gather, dt, offsets, v_nmo): cmp_gather: 2D array, shape (n_samples, n_offsets) dt: 采样间隔秒 offsets: 炮检距数组米 v_nmo: 动校正速度米/秒常速简化版 返回校正后的道集同 shape n_samples, n_offsets cmp_gather.shape corrected np.zeros_like(cmp_gather) for j in range(n_offsets): x offsets[j] for i in range(n_samples): t i * dt # 当前样点的旅行时 t0_sq t**2 - (x**2) / (v_nmo**2) if t0_sq 0: continue # 切除拉伸部分 t0 np.sqrt(t0_sq) idx t0 / dt i0 int(np.floor(idx)) if i0 0 or i0 n_samples - 1: continue frac idx - i0 # 线性插值 corrected[i, j] (1 - frac) * cmp_gather[i0, j] frac * cmp_gather[i0 1, j] return corrected这段代码的逻辑是对每个炮检距的每一道遍历时间样点根据时距方程反算 t₀然后把原始振幅通过线性插值搬到 t₀ 对应的样点位置。参数说明dt是采样间隔常见值 0.001 s 或 0.002 soffsets是炮检距数组单位米v_nmo是动校正速度这里用常速简化实际应该随 t₀ 变化。注意t0_sq 0的判断这是为了切除远道浅层的拉伸畸变不做切除的话叠加剖面上会出现高频假象。跑完这段代码你可以把校正前后的道集画出来对比。校正前同相轴是弯的校正后应该基本拉平。如果没拉平先检查速度是不是给错了再检查炮检距单位是不是和速度单位一致——我见过有人炮检距用千米、速度用米每秒结果校正量差了三个数量级。3. 从 CMP 道集到叠加NMO 校正的完整处理链路3.1 数据准备与道集加载的常见格式实际地震数据通常以 SEG-Y 格式存储。做 NMO 之前需要先把 SEG-Y 读进来按 CMP 道集组织数据。常用的工具是segyio或者obspy。下面是一个读取 SEG-Y 并抽取 CMP 道集的示例。import segyio import numpy as np def load_cmp_gather(segy_path, cmp_key21, offset_key37): 从 SEG-Y 文件中按 CMP 号抽取道集 cmp_key: 道头中 CMP 号的字节位置 offset_key: 道头中炮检距的字节位置 with segyio.open(segy_path, r, ignore_geometryTrue) as f: cmp_nos f.attributes(cmp_key)[:] offsets f.attributes(offset_key)[:] traces f.trace.raw[:] # 所有道 unique_cmps np.unique(cmp_nos) gathers {} for c in unique_cmps: mask cmp_nos c gathers[c] { data: traces[mask].T, # 转置为 (时间, 炮检距) offsets: offsets[mask] } return gathers这里cmp_key21和offset_key37是 SEG-Y 标准道头中 CMP 号和炮检距的字节位置不同处理系统可能用不同的字节位置需要根据实际数据调整。读进来之后每个 CMP 道集是一个二维数组行是时间样点列是炮检距。注意traces[mask].T的转置操作segyio 读出来是 (道数, 样点数)转置后才是我们需要的 (样点数, 道数)。3.2 速度函数插值与逐道集校正实际速度分析得到的是离散的 (t₀, v_nmo) 对做 NMO 时需要插值成每个时间样点对应的速度值。常用线性插值或样条插值。from scipy.interpolate import interp1d def apply_nmo_to_gather(gather, dt, t0_picks, v_picks): gather: dict with data and offsets t0_picks: 速度分析拾取的 t0 数组 v_picks: 对应的 v_nmo 数组 data gather[data] offsets gather[offsets] n_samples data.shape[0] t_axis np.arange(n_samples) * dt # 插值速度函数 v_func interp1d(t0_picks, v_picks, kindlinear, fill_value(v_picks[0], v_picks[-1]), bounds_errorFalse) v_at_samples v_func(t_axis) corrected np.zeros_like(data) for j, x in enumerate(offsets): for i, t in enumerate(t_axis): v v_at_samples[i] t0_sq t**2 - (x**2) / (v**2) if t0_sq 0: continue t0 np.sqrt(t0_sq) idx t0 / dt i0 int(np.floor(idx)) if i0 0 or i0 n_samples - 1: continue frac idx - i0 corrected[i, j] (1 - frac) * data[i0, j] frac * data[i0 1, j] return corrected这段代码和前面常速版本的区别在于速度随 t₀ 变化。interp1d的fill_value参数确保插值时超出拾取范围的速度用边界值填充避免出现 NaN。bounds_errorFalse也是同样的目的。实际跑的时候如果速度函数在浅层变化剧烈线性插值可能不够平滑可以换成kindcubic但要注意样条插值在数据点稀疏时可能产生振荡。3.3 叠加与质量监控NMO 校正做完之后把同一个 CMP 道集内的所有道加起来就得到叠加道。所有 CMP 的叠加道按位置排列就是叠加剖面。def stack_gather(corrected_gather): 对校正后的道集做水平叠加 return np.mean(corrected_gather, axis1) # 假设 gathers 是前面加载的 CMP 道集字典 stacked_section [] cmp_positions sorted(gathers.keys()) for c in cmp_positions: g gathers[c] corrected apply_nmo_to_gather(g, dt0.002, t0_pickst0_picks, v_picksv_picks) stacked_section.append(stack_gather(corrected)) stacked_section np.array(stacked_section).T # (时间, CMP)叠加的时候用np.mean而不是np.sum这样叠加道的振幅不会随道数变化方便后续做振幅归一化。质量监控主要看两点一是叠加剖面上同相轴是否连续、能量是否聚焦二是道集上同相轴是否拉平。如果叠加剖面上出现明显的画弧现象说明速度场横向变化剧烈常速或一维速度函数不够用需要考虑横向变速或者各向异性参数。4. NMO 动校正避坑指南五个血泪教训4.1 拉伸畸变不切除导致浅层假频现象叠加剖面上浅层出现强烈的高频噪声同相轴附近有毛刺状干扰。原因远炮检距浅层的动校正量很大校正后波形被拉伸高频成分被畸变放大。如果不做切除这些畸变能量会参与叠加。解决在 NMO 校正后加拉伸切除。切除准则一般用拉伸系数常见阈值是 1.2 到 1.5即校正后波形长度超过原始长度的 20% 到 50% 就切除。代码里可以在t0_sq 0的基础上再加一个判断如果t / t0 stretch_factor则该样点置零。4.2 速度拾取太粗导致同相轴拉不平现象道集上近道拉平了远道还弯着或者反过来远道拉平了近道翘起来。原因速度谱拾取间隔太大速度函数不能精确描述实际动校正速度随 t₀ 的变化。解决加密速度谱的等速度线从 100 m/s 一档改成 50 m/s 甚至 25 m/s 一档。拾取的时候不要只拾能量团中心要沿着能量团的趋势拾取保证速度函数光滑。拾取完了做一次 NMO 校正看道集拉平情况不平就微调。4.3 炮检距单位与速度单位不一致现象校正后同相轴完全乱掉甚至反向弯曲。原因炮检距用千米速度用米每秒公式里 x²/v² 的量纲错了校正量差了 10⁶ 倍。解决统一单位。炮检距和速度都用米和米每秒或者都用千米和千米每秒。在代码里加一个单位检查如果offsets.max() 100大概率是千米单位需要乘以 1000。4.4 各向异性介质中忽略 η 参数现象在页岩或薄互层发育区常规 NMO 校正后远道总是校正不足速度怎么调都拉不平。原因各向异性介质中反射波时距曲线不是标准的双曲线远道旅行时比双曲线预测的更大。需要引入各向异性参数 η 做高阶校正。解决使用各向异性 NMO 公式或者在做速度分析时同时拾取 η 参数。常见做法是在常规 NMO 之后加一个残余时差校正把远道的剩余时差拉平。4.5 叠加前不做道集均衡导致弱反射被压制现象叠加剖面上浅层能量很强深层弱反射几乎看不见。原因地震波传播过程中振幅随时间和炮检距衰减浅层强反射的能量在叠加时压制了深层弱反射。解决在 NMO 之前做振幅补偿常用球面扩散补偿和吸收衰减补偿。也可以在叠加时用加权叠加给深层道更大的权重。我一般会在速度分析之前先做一次道集均衡这样速度谱上深层的能量团也能看清楚。5. 用残余时差校正验证 NMO 质量并迭代优化NMO 校正做完之后怎么判断做得好不好最直接的方法是看道集是否拉平但肉眼判断有主观性。更客观的做法是做残余时差分析。残余时差分析的思路是对 NMO 校正后的道集在一定的时窗内计算不同炮检距之间的互相关如果互相关峰值不在零延迟说明还有剩余时差。把剩余时差拾取出来转换成速度修正量反馈到速度函数里再做一次 NMO 校正。这个过程迭代两到三次道集基本就能拉平。下面是一个简单的残余时差计算示例def residual_moveout(corrected_gather, dt, window0.1): 计算相邻炮检距道之间的残余时差 window: 时窗长度秒 n_samples, n_offsets corrected_gather.shape n_win int(window / dt) residuals [] for j in range(n_offsets - 1): trace1 corrected_gather[:, j] trace2 corrected_gather[:, j 1] # 分时窗做互相关 for k in range(0, n_samples - n_win, n_win // 2): seg1 trace1[k:k n_win] seg2 trace2[k:k n_win] if np.std(seg1) 1e-10 or np.std(seg2) 1e-10: continue corr np.correlate(seg1, seg2, modefull) lags np.arange(-n_win 1, n_win) peak_lag lags[np.argmax(corr)] residuals.append((k * dt, j, peak_lag * dt)) return residuals这段代码对相邻两道分时窗做互相关找到互相关峰值对应的延迟。如果峰值延迟不为零说明这两道之间还有剩余时差。window参数控制时窗长度一般取 0.1 到 0.2 秒太短了互相关不稳定太长了会平均掉时差变化。n_win // 2是时窗滑动步长取一半重叠是为了保证时窗之间的连续性。拿到残余时差之后可以把它转换成速度修正量。近似关系是Δv/v ≈ Δt_residual / (t * x²/(v² * t²))具体推导这里不展开实际处理软件里都有现成的模块。手工做的话可以简单地根据残余时差的符号和大小微调速度函数对应 t₀ 位置的速度值然后重新做 NMO。我自己的习惯是每做完一次 NMO先看道集再看叠加剖面最后跑一次残余时差分析。如果残余时差在 2 ms 以内基本可以接受超过 5 ms说明速度函数还需要调整。这个迭代过程听起来繁琐但做熟了之后两三轮就能收敛。最怕的是速度拾取阶段就偷懒后面怎么迭代都救不回来。希望帮到你。本文还有配套的精品资源点击获取
返回列表