刘 建，沈鸿雁，席井昌，李 勤，赵 静，李 萌

(1.山西省地球物理化学勘查院，山西 运城 044004；2.西安石油大学 地球科学与工程学院，陕西 西安 710065；3.陕西省油气成藏地质学重点实验室，陕西 西安 710065)

陷落柱的存在不仅破坏煤层的连续性，减少煤炭可采储量，影响巷道的掘进和煤层的开采，而且还有可能成为特殊的导水通道，对煤矿生产造成严重的危害。由于陷落柱是一种相对孤立的地质体，具有发育规律性差、空间形态不规整、规模悬殊大、隐蔽性强等特征，导致陷落柱探测难度较大。目前，对于长轴直径25 m以上的陷落柱，三维反射地震技术的解释准确率也仅有40%～50%。因此，煤田陷落柱精细探测技术还有待持续加强研究。

由于陷落柱的顶部、内部、围岩3者之间一般存在比较明显的弹性差异，从而会引起地震响应特征发生变化，常以绕射波的形式表现。绕射波的存在，一方面给反射波精确成像带来了极大的挑战；另一方面，绕射波是由地下不均匀地质体引起的地震响应，它们携带了高分辨率甚至超高分辨率的地质信息。因此，绕射波可以作为识别和追踪陷落柱的有效波。自20世纪50年代初，人们就已经感觉到了利用绕射波识别地质异质体(区)的潜力，并开始利用绕射波检测小断层；此外，还利用绕射波分辨率高的特点来识别碳酸岩缝洞型储集体、页岩裂缝和不均匀地质体等。STURZU等的研究结果表明：绕射波成像能有效描绘碳酸盐岩储层中的裂缝；DECKER等认为绕射波成像还能有效改善单个溶洞的水平分辨率，并且可有效识别出低于反射分辨率的异质区域。

然而，在煤田陷落柱的探测中，目前还没有广泛开展类似的研究工作，如果能有效利用绕射波携带的丰富信息进行煤田陷落柱探测，其预测精度将有可能得到进一步提升。为此，笔者开发了一种利用中值阻滤波分离绕射波的方法，并将分离出来的绕射波单独进行偏移成像，其目的是充分利用绕射波分辨率高、能直接指示地质异质体(区)存在的优势，进一步提高煤田陷落柱的预测精度。

1 方法原理

脉冲信号与噪声的主要区别在于它们的统计分布特征，如果要实现这类信号与噪声分离，最好的处理办法是采用基于顺序统计量的滤波器。中值滤波是一种基于排序统计理论来抑制噪声的非线性信号处理方法，其主要特点是运算简单，而且处理速度较快，在滤除噪声(特别是脉冲噪声)的同时，能较好地保持信号的细节信息。最初该方法是针对离散数据平滑问题提出来的，随后广泛应用于数字图像处理和信号处理领域并得到了快速发展。目前在地震数据处理中，中值滤波也得到了应用，主要用于剔除野值、VSP上下行波场分离等。

1.1 中值滤波

假设存在一组信号序列(=1，2，…，)，如果滤波窗口长度为，则对信号系列中第点的中值滤波过程为：取以第点为中心的个样值(=2+1，为正整数，即在窗口左右两边各取个样值)作为输入，并对这个样值按大小顺序重排，取重排后个数据中心位置的样值作为该点滤波的输出。考虑到滤波的边界效应，一般采用式(1)，并通过滑动窗口来实现滤波。

(1)

一般为奇数，若为偶数，输出应为中间2个样值的平均值；为原始信号。

1.2 中值阻滤波

借鉴频域带阻滤波去噪的方法原理，可定义中值阻滤波。中值阻滤波就是把中值滤波压制的噪声当成有效信号被保留，而其他成分则当成噪声被压制掉，其基本原理是原始信号减去中值滤波结果(式(2))。对于一个含噪声的一维信号，如果有效信号不具备顺序统计特征(如脉冲型信号)，而要剔除的噪声又具有良好的统计特性，此时就可采用中值阻滤波来实现信号与噪声的分离。

(2)

式中，为中值阻滤波结果。

1.3 绕射波分离处理流程

在分离绕射波的过程中，彻底去除反射波是绕射波去噪的难点。鉴于一维中值阻滤波可有效分离脉冲型信号，如果通过NMO将反射波校平，此时反射波具有最佳的横向相干性，而绕射波因存在剩余时差致使其横向相干性仍然较差。在不同地震道上相同采样时间的一维地震信号中，绕射波不具备连续性分布特性，其表现特征类似于脉冲型信号，因此它们几乎不存在排序统计特性，而反射波的排序统计特性则最好，通过一维中值阻滤波就可实现绕射波与反射波的分离。在滤波过程中，为了充分满足中值滤波统计特性的要求，即信号采样点要达到一定数量，笔者选用NMO后的炮集记录进行滤波处理。

为了实现绕射波的提取及利用，笔者建立了绕射波分离及成像处理流程(图1)。具体来说，输入地震记录后，需要对原始地震数据进行预处理，具体包括静校正、能量补偿、干扰波压制等；随后进行速度分析和共炮点域NMO，接下来对NMO后的炮记录进行波数谱分析获得中值阻滤波窗口因子，然后依次对不同地震道上相同采样时间的一维地震信号实施中值阻滤波提取绕射波；提取的绕射波地震记录中可能包含有因NMO导致的大偏移距波形畸变，因此需要对这些噪声进行切除；反NMO后就可获得绕射波地震记录；随后以反射波速度分析获得的速度为基础进行成像速度建模，获得地震波成像速度模型，并应用于绕射波叠前偏移成像和反射波偏移成像，就可分别获得用于地质解释的地震成像结果。

图1 绕射波分离及地震成像处理流程

需要强调说明的是：中值阻滤波分离绕射波可在共炮点(CSP)域、共接收点(CRP)域和共中心点(CMP)域中进行，不同处理域分离绕射波的处理流程基本相同，唯一差别是在进行中值阻滤波分离绕射波之前需进行道集选排。不同处理域分离绕射波还需满足批量化处理条件，即每一个道集的道数需相同，笔者选择在共炮点(CSP)域中分离绕射波。此外，对于含多炮的地震测线而言，在地层界面横向变化相对单一的情况下，采用相同的滤波参数就可连续处理整条测线数据，否则需采取分段处理的措施，即不同地层界面特征段的道集采用不同的滤波参数处理。

2 陷落柱模型处理与分析

2.1 正演模拟

为了检验本文方法的有效性，建立了一个包含4层的地质模型，其中第3层是厚度为6 m的煤层，在第2，3，4层中含有尺寸不等的4个陷落柱，模型参数如图2和表1所示，陷落柱几何参数详见表2，采用有限差分解波动方程的技术模拟地震波场。震源子波选用雷克子波，子波主频为60 Hz。采用检波器固定，炮点移动的方式采集数据。地震数据采集参数包括：道距5 m，炮距10 m，201道/炮，共采集101炮，采样时间1 ms，采样长度700 ms。正演模拟获得的原始单炮地震记录如图3(a)所示，由陷落柱引起的绕射波较发育。在数据处理时，为了贴近实际，在原始地震记录中加入了随机噪声，信噪比(SNR)约为5∶1(图3(b))。

图3 单炮地震记录

表1 地质模型参数

表2 陷落柱参数

图2 陷落柱地质模型

2.2 滤波因子确定

在中值阻滤波分离绕射波的过程中，滤波质量的好坏取决于窗口因子(滤波点数)参数的选择。由波数谱分析理论可知：波数表示单位距离内波的个数，如果假设以地震道为单位(不是距离)，通过提取绕射波发育段(目的层)的地震信号进行波数谱分析就可获得绕射波在空间域上的分布道数，以此就可定量确定中值阻滤波的窗口因子参数。

依据图1的绕射波处理流程对该套数据进行处理。从NMO结果(图3(c))来看，反射波已被校正平，绕射波仍然不同程度地没有被校正平，横向相干性差别依然存在较大差距，大偏移距存在因NMO导致的波形畸变现象。由NMO后的单炮地震记录(图3(c))分析可知，绕射波发育于450～700 ms，为了在波数分析时尽可能地减少其他噪声的影响，提取=550 ms的地震数据进行分析(图4(a))；从波数谱分析(图4(b))可知，主波数带为0.049～0.112(其中下切波数=0.049 m，上切波数=0.112 m)，利用式(3)可将波数转换为道数，获得绕射波发育的道数为9～21道。依据最大道数原则，最终确定中值阻滤波窗口因子参数为21道。

图4 地震记录波数谱分析(图3(c)，t=550 ms)

=1

(3)

式中，为中值阻滤波道窗参数(地震道数)；为等效波数。

2.3 绕射波分离

图5(a)为采用沿地震道方向上滤波窗口因子为21的一维中值阻滤波对图3(c)处理后获得的地震记录；对NMO导致的波形畸变切除后进行反NMO，就获得了图5(b)所示的绕射波地震记录，从绕射波分离结果来看，绕射波提取完整，较强的反射波则得到了有效压制；图5(c)是提取绕射波后的反射波地震记录，反射波没有丢失，绕射波分离较彻底。

笔者分别取滤波窗口因子参数为5，11和51对图3(c)做中值阻滤波实验，相对于滤波窗口因子为21的滤波结果(图5(b))，滤波窗口因子为5和11的中值阻滤波结果中(图6(a)，(b))，绕射波提取不完整，尤其低频有效信息丢失严重，而高频噪声又得不到有效压制；在滤波窗口因子为51的中值阻滤波结果(图6(c))中，虽然绕射波分离完整，但同时也提取了部分反射波。因此，滤波窗口因子为5，11和51的中值阻滤波均不利于绕射波的提取。由此可见，在采用中值阻滤波分离绕射波时，滤波窗口因子参数太小，绕射波丢失严重，高频噪声得不到有效压制；滤波窗口因子太大，虽然高频噪声得到了有效压制，但反射波也将随着增强。采用本文方法分离绕射波时，要慎重选择好滤波窗口因子参数。

图5 绕射波分离(第1炮)

图6 不同窗口因子参数分离的绕射波比较(第1炮)

2.4 成像结果对比分析

绕射波的利用只有通过偏移成像处理后才能得以实现，为此笔者采用Kirchhoff叠前时间偏移成像技术对分离出来的绕射波进行了偏移成像处理，并与全波场成像(包括水平叠加、叠后偏移和叠前偏移)结果进行了比较。从成像结果来看，在全波场水平叠加剖面(图7(a))中，陷落柱的地震响应可辨识，表现为抛物线形状特征，然而其纵、横向分辨率均较低，无法准确判断陷落柱的具体位置和规模；对于全波场叠后偏移成像剖面(图7(b))，虽然绕射波得到了有效收敛，但仅有1号、2号和3号陷落柱依稀可辨，纵、横向分辨率也很差；对于全波场叠前偏移成像结果(图7(c))，虽然有效改善了陷落柱的成像质量，1号、2号和3号陷落柱信息均归位较准确，纵、横向上的分辨率也较高，然而受煤层强能量的反射波影响严重，小尺度的4号陷落柱仍然无法识别。绕射波场叠前偏移成像(图7(d))明显优于全波场成像结果，4个陷落柱的地震响应均得到了有效成像，成像位置与地质模型(图2)中的陷落柱位置对应准确。该模型的成像结果证明本文提出的绕射波分离方法是有效的，同时也印证了绕射波成像结果具有较高的分辨率，对煤田陷落柱识别的优势明显高于全波场成像结果。

图7 地震波成像结果比较

3 实际资料应用

利用本文方法处理了山西Z煤田约1.5 km的三维地震资料，采用一炮八线制的方式采集地震数据，采集参数包括：炮线距8 m，炮距2 m，接收线距6 m，道距1 m，1 024道/炮，采样时间1 ms，采样长度700 ms。图8(a)为典型的原始炮记录，由于地形起伏的影响，存在静校正量；此外，该地震记录中面波较发育，能量较强，是主要的干扰噪声。经过静校正和去噪处理后(图8(b))，隐藏在面波下的反射波已经凸现出来了，同时还发育能量较弱的绕射波。采用中值阻滤波方法对该套数据进行处理获得了图8(c)所示的绕射波地震记录，反射波得到了有效压制，并完整地提取了绕射波。

图8 绕射波分离前后比较

图9为该三维地震资料的成像结果，其中图9(a)为全波场成像结果，图9(b)为绕射波场成像结果。从图9(a)可看出，煤层引起的地震同相轴连续性较好，但局部出现同相轴错断现象，从=245 ms的全波场成像结果水平切片(图10(a))来看，地震振幅不均匀，表明煤层界面存在起伏变化情况。绕射波场偏移成像结果(图9(b))进一步明确了煤层起伏变化特征，从=245 ms的绕射波成像结果水平切片(图10(b))来看，存在多个正或负的局部幅值异常能量团，表明煤层存在多个凸起、凹陷区。

图9 地震成像结果比较

图10 水平切片比较

提取=120 m和=200 m的垂直切片(图11)并进一步对比分析发现，图11(c)中的A处存在同相轴错断，然而错断同相轴之下的负幅值同相轴是连续的；此外，在图11(d)中对应于A异常区域表现为一能量较弱的负幅值异常小区域，而在图10(b)中对应于A异常区域表现为正幅值包围的近椭圆状负幅值异常特征，由此说明A异常不是断层；结合图10(b)中的异常棱廓特征，可判定A为一个直径约为12 m的陷落柱。图11(b)中的B处存在一个较强负幅值的异常，在水平切片中(图10(b))，其地震响应特征与A异常相似，结合图11(a)的全波场地震响应特征来看，其由煤层存在倾角所导致，并非陷落柱引起。此外，结合绕射波和全波场成像结果对比分析，图10(b)中的其他正幅值异常均可判定为由煤层凸起(区域)引起的局部异常特征。经后期煤层采掘验证，地震异常A确实由一采空区引起，其长轴长约12.6 m，短轴长约11.7 m，与地震解释结果基本吻合。

图11 垂直切片比较

4 结 论

(1)由于绕射波与反射波的时距规律不相同，NMO后利用基于顺序统计特征差异的中值阻滤波可有效分离绕射波与反射波。

(2)利用中值阻滤波分离绕射波时，滤波窗口因子(滤波点数)参数对绕射波分离质量起着关键性作用，窗口因子参数过小，容易丢失有效绕射波，而窗口因子太大则容易残留反射波，一种可取的办法是通过目的层段的波数谱分析来定量确定中值阻滤波的窗口因子参数。

(3)绕射波单独成像可有效改善陷落柱的成像质量，尤其对改善小尺度陷落柱的成像质量具有绝对的优势。

需要补充说明的是，本文的目的不是完全抛弃反射波，而是通过绕射波单独成像，并与全波场成像结果联合解释，进一步提高煤田陷落柱的识别精度。