煤矿隐蔽致灾地质体地震绕射波探测方法

赵惊涛彭苏萍陈宗南柳倩男

1. 煤炭资源与安全开采国家重点实验室,北京 100083；2. 中国矿业大学(北京)地球科学与测绘工程学院,北京 100083；3. 空军研究院工程设计研究所,北京 100068

在煤矿开采过程中,隐蔽小尺度断层和陷落柱等不连续地质体破坏了煤层的连续性,易诱发透水与瓦斯突出等事故,严重威胁煤矿安全生产[1]。断层易使顶板岩层破碎,支架失稳,陷落柱迫使采煤机械停产搬家,由此造成的储量损失高达15% ～30%。

长期以来,我国地质、采矿界科技人员对此类问题十分重视,进行了大量的研究工作,但由于小尺度断层和陷落柱隐伏性强、分布规律性差,仅靠地质和钻井资料无法满足煤矿生产需求,矿井设计存在盲目性,制约矿井安全高效生产。三维地震勘探方法广泛应用在断层和陷落柱等不连续地质体空间分布的研究上。但传统三维地震勘探以反射波为主。在此基础上开展的不连续地质体检测方法,如相干技术[2]、地震属性分析方法[3]和蚂蚁体追踪算法[4]等,主要适用于大尺度断层和陷落柱,分辨率难以满足煤矿安全高效生产需求。一般地,煤矿隐蔽性小尺度地质异常体主要指长轴直径小于20 m 的陷落柱和落差小于5 m 的断层。据统计,传统三维地震勘探技术,对于长轴直径20 m 以上陷落柱解释的准确率约为40% ～50% ,在探明落差小于5 m 断层上仍存在极大挑战[5]。因此,如何提高隐蔽小尺度断层和陷落柱探测精度,一直都是地震勘探的难题。

从地震波场看,绕射波是地质体尺度小于一个菲涅尔带的地震响应,携带了关于断层、尖灭点、陷落柱和岩溶洞穴等小尺度地质体信息,因此有效利用地震数据中的绕射波,对提高三维地震成像分辨率至关重要。但由于绕射波能量较反射波弱1 ～2个数量级,因此绕射波利用难度较大,需单独分离绕射波。近年来,利用绕射波识别小尺度不连续地质体(断层、岩溶洞穴和陷落柱等)的研究日益增多,研究热点集中在三维叠前绕射波分离与成像方法。绕射波为弱信号,通常淹没在强反射波背景中,因此分离绕射波是研究难点,在方法上主要分为数据域和成像域两类。

数据域分离方法研究较早,是指在地震数据偏移之前通过信号变换等方法分离绕射波,主要利用数据域反射波和绕射波几何形态差异,通过压制平滑连续的反射波,分离曲率大、动力学特征连续性差的双曲绕射波。Asgedom 等[6]研究基于相似度和多重信号分类准则的绕射波分离方法,用于增强绕射波能量。Figueiredo 等[7]利用模式识别技术研究绕射波自动成像方法,通过k最近邻值技术区分反射波与绕射波。为克服地震采集数据有限孔径和带宽限制,Gelius 等[8]利用共反射面元叠加和多重信号分类方法,探讨了窗口导向的高分辨率绕射波成像方法。朱生旺等[9]利用局部倾角滤波和预测反演方法实现了炮集数据绕射波分离。蒋波等[10]利用反射层拉平方法研究了绕射波分离问题。刘太臣等[11]研究了基于奇异值分解的共偏移距绕射波分离方法。Zhao 等[12-13]利用一致渐近绕射理论和双指数拟合最优化方法,研究了炮集数据域绕射波分离与成像方法,以及基于马氏距离和绕射波振幅衰减的绕射波自动分离与成像方法。Yu等[14]通过正则化平面波破坏方法,研究了共偏移距数据绕射波分离方法,解决了在噪声背景下反射波倾角估计稳定性问题。

成像域分离方法需要提供偏移速度模型,利用反射波和绕射波在成像域的运动学和动力学特征差异来分离绕射波。孙赞东和白英哲[15]提出了一种基于倾角域反射波广义拉东谱的绕射波分离方法。Zhang 等[16]通过在角度域压制菲涅尔带,实现了叠前时间域绕射波成像。刘斌等[17]在倾角域研究了基于局部倾角估计的绕射波分离与成像方法。徐辉[18]利用张角和地层倾角成像道集中反射波和绕射波特征差异,通过倾角信息分离了绕射波。Silvestrov 等[19]利用逆时偏移方法产生共成像点道集,并采用Radon 变换分离绕射波,但由于逆时偏移对速度模型依赖性极强,该方法应用难度较大。

综上所述,数据域分离方法难以有效处理复杂波场情况下绕射波分离问题,当地震波相交和相切时,很难保证绕射波极性和振幅特征完整性,致使在后续绕射波成像过程中小尺度地质信息极易缺失；成像域分离方法虽能够利用偏移算子处理地震波场相交情况下绕射波分离问题,但对速度模型依赖性强,复杂地质条件适用性不强。

本文基于曲波稀疏变换和平面波分解,研究三维叠前炮域多参数稀疏优化绕射波分离方法,解决在地震波场相交和相切情况下绕射波分离的难题,使得分离的绕射波具有波形一致性和完整性,进而实现绕射波单独成像,提高不连续地质体成像分辨率。该研究一定程度上解决了煤矿小尺度断层和陷落柱等隐蔽致灾地质体探测难题,为我国煤矿安全高效开采提供技术支撑。

1 绕射波弱信号提取方法

传统绕射波分离方法一般利用平面波破坏滤波器,在共偏移距数据分离绕射波。该方法要求观测系统规则化,但一般近偏移距和远偏移距中地震数据分布不规则,绕射波难以有效分离,且在地震波场交叉或相切情况下,该方法也不能保证绕射波分离的完整性。

本研究以曲波变换表示绕射波,结合平面波分解,构建三维叠前炮域多参数稀疏优化绕射波分离方法,解决在地震波场相交和相切情况下绕射波分离难题。首先,对三维叠前炮集数据动校正,消除偏移距造成的时差,使得反射波表现为线性特性,而绕射波为双曲线特征；其次,利用曲波变换和平面波分解,构建多参数稀疏最小二乘目标函数,其中曲波变换系数为稀疏约束,对应于绕射波,而平面波分解算子对应的倾角场为平滑约束,对应于反射波；最后,通过反动校正算子,完成三维叠前炮域绕射波分离。

1.1 多尺度多方向曲波变换方法

基于Wrapping 算法的曲波变换,在每个象限每个角度的网格是一样的,每个曲波变换都有正确的方向,定义笛卡尔坐标下的曲波系数为

式中,˜(ω)为信号f的傅里叶变换；为剪切矩阵Sθl的逆矩阵；j、l、k分别为尺度、角度和角度平移量；ω为角频率。

定义原点周围L1,j×L2,j的区域:

基于Wrapping 算法的curvelet 变换方法,改变了数据的符号,数据d[n1,n2]的Wrapping 可以简化为Wd[n1modL1,j,n2modL2,j]=d[n1,n2],(n1,n2)∈Pj,l。

这种表示将原来的点(n1,n2)映射到原点附近新的位置。对于角度θ∈(π/4,3π/4),Wrapping算法是类似的,需要交换两个坐标的位置。

基于Wrapping 算法的快速曲波变换如下:

第一步:对f[t1,t2]∈L2(R)采用二维快速傅里叶变换得到˜[n1,n2],其中t1,t2为空间变量；n1,n2为频率变量,-n/2＜n1,n2＜n/2。

第二步:对每个尺度j和方向参数l,构造乘积,从而实现局部化。

这个算法的计算复杂度为O(n2logn)。在实际计算中,计算量不会超出二维快速傅里叶变换的6 ～10 倍。该变换的逆变换为其转置,属于一种连续变换。

1.2 平面波分解方法

平面波分解利用局部平面波方程表征地震数据构造特征,是一种基于地震局部倾角估计的预测方法。该算法依次从相邻地震道预测当前地震道,预测算子按照局部倾角模式对原始输入地震数据进行平滑估计,因此可以预测地震反射波。平面波分解方法可定义为线性算子,表达式如下:

式中,r为预测的地震反射波数据；s为动校正后的炮集地震道；L2(σ)为线性平面波算子；σ为地震倾角斜率。

式中,I为单位矢量；σ为局部倾角斜率矢量；Pi,j(σi)表示由第i道地震数据预测第j道数据。

1.3 多参数稀疏优化绕射波分离方法

基于曲波稀疏变换和平面波分解构建的多参数稀疏优化绕射波分离模型如下:

式中,‖·‖1为l1稀疏约束范数；d为三维叠前炮集动校正数据；L1为曲波变换算子；m1为其变换系数；m2为其倾角场；α为正则化参数。

虽然基于l1范数误差函数具有非平滑性质,但数值极小值难以计算,并且l1范数误差函数在小误差处也会按照大误差处进行“严重对待”,因此适用性有限。Huber 范数收敛更快,且由于与较少参数有关而更容易操作,可用于近似表示l1范数:

式中,ξ为可调节的正则化参数；hξ为分段函数,稀疏约束项。

引入拉格朗日算子后,多参数稀疏优化绕射波分离模型等价于求解如下无约束问题:

式中,β为正则化因子。

该模型可通过交替方向迭代法求解,步骤如下:

由于三维叠前炮集地震数据规模较大,因此采用拟牛顿法求解子问题,拟牛顿法只需计算一阶偏导数,计算简单且运算量小,比解析求得Hessian矩阵更加高效稳定。

基于curvelet 稀疏变换的三维叠前炮域多参数稀疏优化绕射波分离方法的流程如图1 所示。

图1 三维叠前炮域多参数稀疏优化绕射波分离流程Fig.1 Flow chart of three-dimensional pre-stack diffraction separation with multi-parameter sparse optimization

2 数值模型实验

三维地质模型由水平反射界面、断层和散射点组成,如图2 所示,图中A 和B 标记为散射点,通常描述地质上的岩溶洞穴或小尺度陷落柱。水平反射界面将产生较强能量的反射波,而断层边界和散射点产生的地震响应分别为边缘绕射波与散射波。采用三维Kirchhoff 正演算法合成地震记录,射线追踪路径如图3 所示,用于模拟地震波传播规律。激发震源采用主频为30 Hz 的Ricker 子波,其中全波场(反射波与绕射波)单炮记录的一条接收线如图4(a)所示,在该图中强反射波对应于水平反射界面,1 和2 为断层上下边界产生的边缘绕射波。通常,断层产生的绕射波能量较弱于反射波,但强于散射点,因此在图4(a)中,更为微弱的散射波几乎不可见。利用本文提出的三维叠前炮域多参数稀疏优化绕射波分离方法去除强反射后,得到的绕射波炮集记录如图4(b)所示,该图中强反射波去除干净,1 和2 对应的绕射波特征清晰,A 和B对应的散射波弱信号突显。在反射波与绕射波相切或交叉情况下,传统平面波滤波器的绕射波分离方法[9]由于波场估计不准,使得分离的绕射波不完整,出现绕射波顶部残缺或者交叉位置波场不连续等问题。本文提出的方法在绕射波分离结果图4(b)中较好地解决了此类问题。利用三维叠前时间偏移算法分别处理全波场和绕射波地震数据,获得成像结果如图5 所示。图5(a)能够清晰刻画水平反射界面和断层边界,但散射点无法表征；图5(b)强反射去除干净,散射点和断层断点成像清楚。由于图5(b)中的断层断点成像存在左右极性反转,而散射点无此特征,因此该现象可作为区分断层和陷落柱两类小尺度隐蔽地质体的关键因素。

图2 三维地质模型Fig.2 3D geological model

图3 三维模型单炮射线追踪路径Fig.3 3D ray tracing path exciated by one shot

图4 模拟数据炮集Fig.4 Synthetic shot gather

图5 模拟数据偏移剖面Fig.5 Profile of synthetic migation result

数值模型测试表明,本文提出的基于曲波和平面波分解的三维叠前炮域多参数稀疏优化绕射波分离方法,能够有效去除反射波,突显弱绕射信息,有助于实现小尺度不连续地质体的高精度探测。

3 实际资料应用

3.1 研究区基本地震地质条件

在煤炭生产中,断层、陷落柱及采空区等不良地质体,是导致煤矿事故发生的主要灾害源。因此,在矿井采掘活动前,需要查明开采区域地质构造情况,以确保矿井安全生产。研究区为山西阳泉煤业二矿九采扩区,含煤地层从老到新依次为石炭系上统太原组、二叠系下统山西组；下伏地层为石炭系中统本溪组,上覆地层为二叠系下统下石盒子组。构造基本形态为一走向北西、向南西倾伏的单斜构造,以褶皱构造为主,断裂构造次之。陷落柱剖面形状多为上小下大的不规则柱体,局部呈葫芦状,柱体内岩性杂乱。在陷落柱边缘,煤层及顶板裂隙密集,有时伴有小断层或煤层底板牵引下弯现象。本次主要勘探目的层为8 号和15 号煤层,煤层发育较稳定、结构较简单,与围岩物性差异大,煤岩之间能构成较好的波阻抗界面,可形成能量较强、连续性较好的反射波,具有较好的地震地质条件。

3.2 实际资料绕射成像结果

野外三维地震采集为中点激发12 线16 炮线束状观测系统,接收道数840,道间距10 m,接收线距40 m,最小炮检距5 m,最大炮检距458 m。经过观测系统加载、静校正、噪声去除和反褶积等一系列常规处理后,得到的全波场炮集记录的部分数据如图6(a)所示,其中强振幅连续好的双曲线同相轴对应于较为平缓的强波阻抗反射界面,如煤层,而断层和陷落柱则表现为弱绕射波信号,通常在炮集中很难被直接观察到。利用本文提出的三维叠前炮域多参数稀疏优化绕射波分离方法,消除强反射信号屏蔽作用后,得到的绕射波炮集记录部分如图6(b)所示,该记录去除了双曲形态的强反射信号。为进一步对分离的绕射波信号进行质量控制,将绕射波叠加剖面图与常规全波场叠加剖面对比(图7),绕射波叠加剖面中呈现为线性的反射波同相轴去除较好,而双曲形态的绕射波弱信号得到了突显。分别对全波场和绕射波炮集记录进行三维叠前偏移处理,成像结果如图8 所示。

图6 绕射和全波场炮集Fig.6 Diffraction and full-wave shot gather

图7 绕射和全波场叠加Fig.7 Diffraction and full-wave stacked section

图8 绕射和全波场叠前偏移Fig.8 Diffraction and full-wave migration

陷落柱塌陷引起地层不连续性,在三维地震数据体上通常表现为同相轴扭曲和分叉、能量变弱、连续性变差。常规三维叠前偏移剖面解释出的DX6、DX7 和DX8 三个陷落柱,发育在8 号和15 号煤层中,形状近椭圆状,但陷落柱边界形态不清楚,尤其是DX7 陷落柱。在绕射波三维叠前偏移剖面中,3 个陷落柱边界形态断点清晰,且陷落柱内部特征丰富,陷落柱边缘伴有小断裂,说明隐蔽地质体探测精度得以提高。三维叠前绕射波偏移剖面显示DX6 陷落柱在8 号煤相交位置处存在短弱同相轴错断,在8 号煤层以下还存在短反射同相轴；在15 号煤层以下发现了断裂。常规三维叠前偏移剖面解释出的DX7 陷落柱顶点位置在绕射波成像剖面中存在一段弱同相轴。绕射波成像剖面中给出的DX8 陷落柱与8 号和15 号煤的相交位置更为清晰,内部的同相轴下弯变形严重,一定程度上反映了内部的塌陷情况。

3.3 绕射波沿层属性分析

地震数据体中包含大量的地质信息,不同的地震属性可能与地质情况存在关联,能够增加对断层、陷落柱和挠曲等地质构造的识别精度。沿着15 号煤层,对常规三维叠前偏移数据进行方差属性分析,计算参数选取为3 道×3 道乘法模式,时窗长度为15 ms,方差属性平面如图9 所示。图9 中断层和陷落柱等构造引起的地层出现不连续,大断裂形态清晰,但隐蔽陷落柱和断裂刻画精度不够。三维叠前绕射波偏移数据方差属性平面如图10 所示。图中隐蔽断裂清晰,小尺度陷落柱刻画精度提高；主测线201 和联络线251 交叉位置显示出若干环形区带,该特征符合陷落柱平面展布规律,而附近存在的小尺度线性断裂,在常规偏移数据方差平面图上并无显示。同时,还出现了一系列与主断层具有相同走向的线性小尺度断裂,呈现北东向间断分布特征。这种线性断裂进一步揭示了岩溶发育与断层活动的伴生关系。煤矿应用效果表明,本文提出的三维叠前炮域多参数稀疏优化绕射波分离方法,能够有效去除煤层强反射信息屏蔽作用,分离出完整的绕射弱信号,叠加剖面中绕射波的双曲形态进一步表明分离出绕射波的可靠性。在绕射波地震解释上,通常需要全波场成像结果提供宏观的地质背景,在此基础上进一步利用绕射波成像结果揭示隐蔽地质体,如陷落柱和断层等。从成像剖面上看,绕射波成像结果对陷落柱的边界形态刻画更为清晰,内部细节描述更为丰富。从沿层属性上看,绕射波成像数据对隐蔽断裂的揭示更为细腻,对陷落柱的探测能力有很大提升。

图9 15 号煤常规偏移数据方差平面Fig .9 Variance map of conventional migration data of No.15 coal

图10 15 号煤三维叠前绕射波偏移数据方差平面Fig.10 Variance map of 3d pre-stack diffraction migration data of Coal No.15

4 结 论

(1) 常规叠前地震绕射波分离方法利用平面波破坏滤波器在共偏移距进行,该方法要求观测系统规则化,并且在地震波场交叉或相切情况下很难保证绕射波分离的完整性。基于曲波变换和平面波分解,本文构建了三维叠前炮域多参数稀疏优化绕射波分离方法,解决了地震波场相交和相切情况下绕射波分离问题。

(2) 在地震高分辨率探测上,本文提出的方法能够有效消除强反射波,分离出完整的绕射波/散射波弱信号,提高断层断点和散射点成像精度,有助于实现小尺度不连续地质体的高精度探测。

(3) 在绕射波成像结果解释上,需要结合全波场成像结果的宏观地质背景,进一步揭示隐蔽致灾地质体,其中极性反转特征可用于区分断层断点和散射点。

(4) 绕射波成像虽然具有高分辨率探测能力,但在地质认识上还存在很多研究难点,如绕射波高分辨率成像结果的地质可靠性分析、绕射波成像智能解释等。利用绕射波多方位性波场特征,结合人工智能技术建立多属性模型是未来的研究方向,将有助于进一步提升煤矿小尺度不连续隐蔽致灾体探测能力。