可控震源的激发频率及激发时间可根据勘探目标人为设置, 因而具有高效、经济、环保、灵活等优点, 被广泛应用于陆上地震资料采集^[1-2]。但是, 大吨位或超大吨位P波可控震源, 只能在地表使用, 激发垂直向下的地面压缩力^[3-5]。可控震源产生的地面力(可控震源子波)在半空间地球弹性介质中传播, 不仅会产生有效反射P波, 还会产生被视为噪声的各种相干波, 如频散面波等。尤其对于面波, 当遇到地表不规则地形(如建筑、沙丘、山地等障碍物)和近地表横向剧烈变化的不均匀体时, 会产生强散射面波或背散射波^[6-8], 在共炮点道集上形成以直达面波为边界的强能量三角形区, 其能量约占可控震源激发总能量的80%~90%。此外, 可控震源机械系统非线性和震板刚度不够等因素引起的谐波畸变^{[4, 9]}, 以及可控震源震板与地表耦合不佳造成的较强能量谐振^[10]以及声波干扰^[11]等, 使得三角形区内的强能量噪声更加复杂。

压制面波的方法很多。其中, 检波器组合方法已被证明是一种消除面波的有效方法, 但却以牺牲高频成分为代价^[12]。f-k、f-x以及f-p等转换域滤波方法, 虽然对衰减直达面波有效, 但对于陆上可控震源地震采集而言, 由于面波速度低, 对空间采样要求高, 往往产生空间假频; 散射面波频带较宽, 与有效反射波频带重叠, 不满足转换域滤波方法的假设。基于Born近似的迭代反演预测-消除方法^[13], 不仅需要昂贵的计算资源, 而且, 由于面波依赖于近地表高度变化的介质特性, 难以模拟频散、低频、强振幅有关的面波特征, 尤其在沙漠等地区, 散射波能量比反射波要强得多, 因此, 实际应用是沙漠区可控震源地震资料处理难点之一。

DONG等^[7]率先提出面波干涉预测与消除方法, 基本思路是: 通过接收点或炮点间地震干涉法^[14-16](互相关型)从实测地震资料中恢复一对接收点或炮点间的面波响应, 然后采用最小平方匹配滤波方法消除实际地震道集上的面波。由于面波走时呈线性, 且所有震源均为稳相位震源, 因此, 叠加所有共炮点道集上两个地震道之间的互相关可以增强恢复的面波响应; 而对于来自地下反射面的反射波, 因其走时呈非线性, 且只有少数的稳相位震源, 使得叠加后恢复的反射波响应欠估计^[17]。因此, 相比于转换域滤波等方法, 基于数据驱动的面波干涉预测与消除方法优势明显。但是, 该方法的精度受制于稳相位震源位置分布, 仅适用于2D测线。XUE等^[18]通过非线性时移校正, 将3D面波干涉预测问题转换为相应的2D问题。HALLIDAY等^[19]展示了衰减介质层状模型散射面波干涉预测的一个例证, 并根据面波干涉积分的稳相位分析, 讨论地震干涉、面波和震源分布之间的关系^[20], 提出衰减介质散射面波地震干涉中“串扰”形成的原因^[21]。根据震源位于一对接收点之间或一端, 分别采用互相关或互褶积形式的干涉方程, 进行可控震源高密度单点接收实际资料的地滚波消除^[22-23]。但是, 该方法要求用于构建虚震源的接收位置处必须有震源激发, 这意味着每一个接收点处都必须有震源激发, 对采集的要求过于苛刻, 成本过高, 实际上难以满足和应用。

本文利用在新疆某工区采集的可控震源资料, 定性分析黑三角噪声的形成原因。在此基础上, 针对散射面波, 采用基于CURTIS等^[24-26]提出的炮-检地震干涉方法, 研发炮-检散射面波干涉预测技术, 可适应任意的陆上观测方式; 同时, 研发基于模式的匹配相减方法, 以预测的散射面波为模式, 消除可控震源资料中的散射面波, 形成了可控震源面波干涉预测与消除技术。将该技术与异常噪声压制技术相结合, 形成了可控震源黑三角综合去噪方法及其模块, 并以新疆某工区的实际可控震源资料为测试资料, 对比了该综合去噪技术与某商业软件的黑三角去噪应用效果。

1 可控震源黑三角噪声分析

2016年, 在新疆某沙漠探区实施了可控震源3D地震数据采集, 工区覆盖面积约为500 km², 覆盖次数为900次, CDP面元大小为25 m×25 m, 记录长度为8 s, 采样率为2 ms; 可控震源3台1次, 采用4~84 Hz的20 s线性升频扫描, 斜坡起始800 ms, 终了500 ms, 驱动70%。工区地表主要为盐碱浮土区和垄状沙丘区(图 1)。

首先, 图 2为对应于图 1中蓝色炮点①、②和③的原始道集及其频谱, 由图 2可见: 炮点能量差异不大, 60 Hz以上高频吸收严重; 浮土区主频低, 频宽约为5~48 Hz, 主要是因为浮土区吸收衰减较为严重, 而沙丘区主频较高, 频宽约为5~57 Hz。但是, 随着地表地形变化, 即从较为平坦的浮土区过渡到有一定高程起伏的沙垄状沙丘区, 黑三角区内噪声愈发严重, 这主要由地表沙丘或近地表横向不均匀体引起的散射面波所致, 且不能清晰地观察到直达面波同相轴。

其次, 对某个典型的原始单炮道集进行线性动校正(速度为500 m/s)及频谱分析, 如图 3所示。线性动校正后, 黑三角区内的不同类型噪声在一定程度上得以显现。主频约为12 Hz的强能量相干噪声, 具有频散的特征, 是沿地表传播的Rayleigh面波, 以及遇到地表不规则地形或近地表横向不均匀体时引起的散射面波; 而有效波区域的能量弱, 肉眼难于观察到反射波, 其主频约为16 Hz。此外, 在可控震源下方, 还有传播速度比面波更低的线性噪声及异常噪声, 特征是强能量, 主频很高, 约为20 Hz, 与地表介质固有频率相近。因此, 推测是可控震源震板与地表介质耦合发生共振(谐振)所致。可见, 可控震源黑三角区内噪声成分复杂, 特征也不同, 且相互干涉叠加在一起, 虽然主频有一定差异, 但频带宽度却相互重叠。因此, 单靠已有的如f-k等转换域滤波技术, 难以有效压制这些噪声。

最后, 根据HERMAN等^[13]的简化近地表散射模型(图 4a), 可以获得2D测线上频散直达面波及其测线上散射点①和③, 测线外散射点②和近地表散射点④的散射面波时距曲线(图 4b)。图 4a中, 炮点和检波器位置分别用SRC和REC表示, 红色箭头表示传播中的直达面波, 绿色箭头表示传播中的散射面波和来自近地表中不均匀体的散射波, 蓝色箭头表示来自地下反射界面的反射波。图 4b中, 来自测线上单散射①或二次散射④①、③①和②①的时距曲线是线性的, 而来自测线外单散射点②或来自近地表的多次散射②③④的散射时距曲线则是非线性的。

2 方法原理

根据CURTIS等^[24-26]提出的炮-检间地震干涉原理(图 5), 利用已知的在震源x₁上激发, 在接收边界S′上接收的波场以及在震源边界S上所有的震源x激发, 在接收点x₂和接收边界上接收点x′接收的波场, 预测(重建或恢复)任意炮检(x₁~x₂)之间的地震响应。重建的炮-检对地震响应不仅可用于补缺失道(如果该道是坏道)或评价该地震道的质量, 更为重要的是, 还可用于消除该地震道中的面波噪声。

但由于实际地震观测资料通常不能满足炮-检地震干涉的前提条件, 预测的面波响应在振幅上存在一定的误差, 因此, 不能从实际地震资料中直接减去预测的面波响应, 而应该基于最小平方或者模式的匹配滤波技术消除面波。可见, 可控震源散射面波干涉预测与消除, 主要包含两项技术, 即散射面波的干涉预测技术和基于模式的匹配相减技术。

2.1 散射面波干涉预测技术

炮-检地震干涉(图 5)可表述为: 假设声波介质中有一个很大半径的接收点封闭边界S′, 其上接收点x′间距不超过Nyquist空间采样准则所要求的间距, 接收来自炮点x₁的格林函数G(x′, x₁), 封闭边界S′上及其外部的介质是均匀的; 接收点边界内另有一个较小半径的震源封闭边界S, 其上震源x间距不超过Nyquist空间采样准则所要求的间距, 其内的接收点x₂和接收点边界上接收点x′分别接收来自震源的格林函数G(x₂, x)和G(x′, x), 那么, 基于相关型和褶积型互换方程, 在高频近似条件下, 可获得频率域炮检(x₁~x₂)间的干涉积分, 即:

$ \begin{array}{c} G({x_2}, {x_1}) - {G^ * }({x_2}, {x_1}) \approx \frac{{4{k^2}}}{{{{(\omega \rho )}^2}}}\int_S {\int_{S'} {G(x', {x_1}) \cdot } } \\ {G^ * }(x', {x})G({x_2}, x){\rm{d}}S'{\rm{d}}S \end{array} $

(1)

式中: G(x₂, x₁)表示重建的炮-检(x₁~x₂)之间的格林函数, k、ω和ρ分别为接收点边界外均匀介质的波数、圆频率和密度, *表示复共轭。

对于干涉积分方程(1), 其左边表示希望重建的炮检(x₁~x₂)间地震响应, 右边实际上由两个干涉积分构成^[24-26]。一是利用已知的实测资料, 即炮点x₁和震源边界S上震源x至接收点边界上所有接收点x′的地震响应, 进行炮点间的相关型地震干涉, 它是对接收点边界S′求积分的, 获得炮点x₁(虚接收点)与震源边界上某个震源x之间的格林函数及其逆时格林函数, 即

$ \begin{array}{c} G({x_1}, x) - {G^ * }({x_1}, x) \approx \frac{{2k}}{{\omega \rho }}\int_{S'} {G(x', {x_1})} {G^ * }\\ (x', x){\rm{d}}S' \end{array} $

(2)

二是利用震源边界上重建的地震响应G(x₁, x)与已知的地震响应G(x₂, x), 进行接收点间的褶积型地震干涉, 它是对震源边界S求积分的, 从而获得从炮点x₁(虚震源)至检波点x₂之间的格林函数及其逆时格林函数, 即:

$ \begin{array}{c} G({x_2}, {x_1}) - {G^ * }({x_2}, {x_1}) \approx \frac{{2k}}{{\omega \rho }}\int_{S} {[G(x, {x_1})} - \\ {G^ * }(x, {x_1})]G({x_2}, x){\rm{d}}S \end{array} $

(3)

2.2 基于模式的匹配相减技术

基于模式的匹配滤波方法, 主要优势在于能够充分利用模式运动学特征及相对振幅关系, 避免中浅部有效信号及直达波的影响。首先, 从原始资料和干涉预测的散射面波中, 利用子波反褶积分别获得其散射面波噪声(N)预测误差算子P_N和原始资料信号(D)预测误差算子P_D, 通过反褶积求取有效信号预测算子P_S, 即在频率域可表示为:

$ {P_S} = \frac{{{P_D}}}{{{P_N}}} $

(4)

其次, 通过最小二乘目标函数, 即:

$ \min {\left\| {{\mathit{\boldsymbol{P}}_S}(\mathit{\boldsymbol{D}} - \mathit{\boldsymbol{\tilde fN}})} \right\|_2} $

(5)

最终求取频率域的最佳匹配滤波算子${\mathit{\boldsymbol{\tilde f}}} $, 可表示为:

$ \mathit{\boldsymbol{\tilde f}} = {({\mathit{\boldsymbol{N}}^\mathit{\boldsymbol{T}}}\mathit{\boldsymbol{P}}_S^\mathit{\boldsymbol{T}}{\mathit{\boldsymbol{P}}_S}\mathit{\boldsymbol{N}})^{ - 1}}\mathit{\boldsymbol{NP}}_S^\mathit{\boldsymbol{T}}{\mathit{\boldsymbol{P}}_S}\mathit{\boldsymbol{D}} $

(6)

利用获得的最佳匹配滤波算子, 就可得到消除面波后的地震道集, 即:

$ \mathit{\boldsymbol{D'}} = \mathit{\boldsymbol{D}} - \mathit{\boldsymbol{\tilde fN}} $

(7)

3 可控震源实际资料应用 3.1 单炮记录试验

根据上述介绍的新疆某工区采集的可控震源资料, 对本文提出的可控震源散射面波干涉和消除技术进行测试, 结果如图 6所示, 图中的红色十字线仅用于指示道集之间的相同位置。图 6b表示散射面波干涉预测结果, 与原始炮道集相比(图 6a), 除了直达面波外, 两者在形态上基本相同, 但细节上稍有不同。图 6c表示将预测的散射面波从原始道集中匹配相减的结果, 可以看到: 消除了大量散射面波, 但仍有少量残余, 主要原因是用于预测散射面波的炮集数量较少。为了能更清晰地进行比较, 对图 6中的炮道集进行了线性动校正, 结果如图 7所示。图 7c清楚地显示已将散射面波从原始炮道集上消除, 这也说明本文方法可有效应用于可控震源黑三角区内的散射面波的预测与消除。

但是, 可控震源黑三角区还有如谐振等其它强能量噪声, 其特点是同相轴不规则、散碎(图 7c), 可视为异常噪声, 为此, 提出了基于奇异谱分析的异常噪声压制方法技术, 将其与本文的可控震源散射面波干涉预测与消除技术相结合, 形成可控震源黑三角综合去噪方法技术。

3.2 实际应用及效果分析

与某商业软件进行了对比测试, 测试流程如图 8所示。测试流程中, 除了“黑三角去噪”所用方法不同外, 分别采用某商业软件的去噪模块和本文所述的可控震源黑三角综合去噪方法技术(简称本文方法), 其它处理步骤相同, 均采用该商业软件的相同模块, 且所用的模块参数也相同。

图 9对比了黑三角区内噪声消除后的CMP道集, 可以看出: 商业软件仅压制了黑三角区内的噪声强度, 被其掩盖的反射波并未被有效增强或恢复; 而本文方法则有效剔除了黑三角区内的强散射面波及谐振噪声, 突显了被其掩盖的有效反射波, 且与黑三角区外的有效反射波一致。

图 10、图 11和图 12对比了黑三角噪声压制后CMP道集的粗叠加结果。图 10对比了近偏移距(0~2 000 m)的叠加结果, 由于商业软件去噪处理后获得的CMP道集上, 在黑三角区范围内看不到任何有效反射波, 因此, CMP叠加后虽然能显现出少量的有效反射波, 但这很可能是因为叠加次数高或者局部地段地震记录近偏移距数据信噪比较高的缘故。由于本文方法有效消除了黑三角区内的大量强散射面波和谐振噪声(图 9), 与商业软件的相比, 本文方法的近偏移距叠加结果揭示了大量与该工区地震反射标志层相一致的反射波同相轴。

图 11对比了中等偏移距(2 000~4 000 m)的叠加结果, 与近偏移距剖面(图 10)相比, 本文方法略有改善, 尤其是叠加剖面中部的绕射波得到突显。因为该工区是深层缝洞型储层, 绕射波是对缝洞进行成像的关键有效波。

此外, 可以从可控震源原始资料的角度出发, 解释本文方法对中等偏移距的叠加剖面仅略有改善的原因。可控震源原始炮集(图 2)表明: 黑三角区内强噪声通常只影响如0~2 000 m的近偏移距地震道, 且随接收时间的增加, 影响范围也不断增大; 而对中偏移距(2 000~4 000 m)影响甚微, 对远偏移距地震道没有影响。这是因为在地表传播的面波和遇到不规则地表地形或近地表横向不均匀体时产生的散射面波传播速度低所致。

图 12对比了全偏移距叠加结果, 本文方法除了突显了剖面中部的绕射波外, 与商业软件的处理结果基本相同。与近偏移距的叠加结果(图 10)形成了鲜明的反差, 其主要原因在于近偏移距对全偏移距叠加的贡献较小, 0~2 000 m的近偏移距道数约占全偏移距道数(0~9 000 m)的25%。

两种方法去噪后的全偏移距数据的叠后时间偏移剖面如图 13所示。与商业软件的处理结果相比, 本文方法有如下改善: ①深大断裂成像更加清晰, 连通性好; ②目标层奥陶系内幕信噪比较高, 目的层成像同相轴更加连续, 对后续资料解释更加有利; ③串珠的成像更加聚焦。这是因为成像垂向分辨率的改善, 主要来自于近偏移距道(图 10)信噪比的显著提高。此外, 近偏移距叠加剖面(图 10)也印证了成像分辨率改善的另一个原因, 即本文方法能有效压制黑三角区的强噪声, 进而有效保留或增强了更多的绕射波, 且未损害有效反射波。

图 14对比了近偏移距数据叠前时间偏移的结果。与商业软件的处理结果相比, 本文方法有效恢复了黑三角区内幕信息, 显著提升成像质量, 且揭示了更加丰富的地质现象。这些改善也反映在全偏移距数据叠前时间偏移剖面(图 15)和CRP道集(图 16)上。

4 结论与展望

沙漠区可控震源资料黑三角区去噪难题一直困扰着此类资料处理人员。本文基于CURTIS等提出的炮-检地震干涉方法, 基于模式的匹配相减方法, 提出了散射面波干涉预测与消除技术, 并与异常噪声压制技术相结合, 形成了可控震源黑三角综合去噪方法技术。

利用新疆某工区的可控震源资料, 定性分析了黑三角区噪声的形成机理, 认为主要由两个因素造成: 一是地表地形或近地表横向不均匀体导致的强散射面波; 二是可控震源震板与地表介质耦合造成的强能量共振或谐振。利用新疆某工区的可控震源资料, 进行了散射面波干涉预测与消除技术的测试, 并与某商业软件进行了对比, 结果表明: ①商业软件仅压制了黑三角区内的噪声强度, 反射波并未被有效增强或恢复, 而本文方法能有效去除黑三角内的强能量噪声, 尤其是散射面波, 突显了黑三角内幕的有效反射波, 且与黑三角区外保持一致; ②本文方法去噪结果显著改善了叠前时间偏移成像剖面的质量, 进而恢复地层真实反射特征, 断裂和断溶体成像效果优于某商业软件。

但是, 黑三角区内含有多种不同类型的噪声, 其机理尚待进一步研究, 仍需进一步探索相应的噪声压制技术。