袁 峰,申 涛,谢晓深,马 丽,汶小岗,2
(1.自然资源部煤炭资源勘查与综合利用重点实验室,陕西 西安 710026; 2.陕西省煤田物探测绘有限公司,陕西 西安 710005; 3.陕西省地质调查院,陕西 西安 710065; 4.西安科技大学 地质与环境学院,陕西 西安 710054)
煤层开采形成的采空塌陷会造成地表含水层水流失、溃沙等灾害[1-5]。采空区中导水裂隙带是否发育到地表和含水层对保水采煤至关重要。因此,有必要对导水裂隙带的发育高度进行精确探测[6-7]。
目前,国内外导水裂隙带的探测方法主要有钻孔冲洗液观测法、井下钻孔注水法、相似材料模拟、钻孔电视法等,探测成本较高,探测结果仅为离散的钻孔资料,孔间推测存在误差[8-10]。三维地震技术具有数据横向连续、纵向分辨率较高的特点,在一定程度上能够弥补钻孔资料的不足[11-12]。
导水裂隙带探测的核心是裂隙识别。目前,常用于裂隙识别的地震属性归纳起来可分为:波形相似类,包括相干体、边缘检测和方差体等;构造几何类,包括曲率属性、倾角属性等;吸收衰减类,包括振幅属性、频率属性和频谱属性等[13-14]。对同一物探手段来说,在实际生产中需要根据其适宜性和有效性进行选择,但无论采用何种解释方法,都存在多解性。综合以往技术实践可以看出,利用多种解释方法进行综合解释,可以提高解释的精度和可靠性。
笔者采用以漏失量为监督数据的深度学习方法融合多种地震属性对导水裂隙带发育情况进行探测。首先,在钻孔冲洗液漏失量观测数据分析的基础上,优选地震属性,结合钻孔冲洗液漏失量观测数据进行有监督的深度学习,融合多种属性信息建立裂隙模型,利用裂隙模型分析确定采动覆岩结构破坏和导水裂隙带的三维空间范围、形态特征、垂向岩石破坏程度定量变化特征及导水裂隙带发育高度,并结合2个工作面的开采时间差异性推断导水裂隙的发育、闭合规律。
研究区位于陕北黄土高原北端,毛乌素沙漠东南缘地带,地形相对平坦,地貌类型主要为风沙滩地地貌。主采3号煤层,3号煤层全区可采,区内可采煤层厚度9.8~10.2 m,煤层由东南向西北缓倾,倾角约0.5°。该煤层属厚煤层,采用复合假顶综合机械化采煤法开采,全垮落式管理顶板,分层开采上分层煤层留底煤,开采厚度为5 m。3号煤层底板标高+1 100~+1 110 m,煤层埋深250~260 m,顶板岩性以中粒砂岩为主。研究区位于30101工作面中段,根据采掘资料,30101工作面煤层已经全部采空,该工作面宽300 m,与其西侧相邻的30102工作面,也已开采完毕,探测范围的选择思路是一方面保证完全覆盖30101,30102工作面,另一方面保证完全覆盖垮落、裂隙及变形带,最终确定探测区为规则矩形,横向长1 000 m,纵向宽500 m,探测区面积0.50 km2,如图1所示,图中蓝色范围为三维地震工作范围。30101,30102工作面开采时间不同,30101工作面煤层采掘距本次野外数据采集时间超过18个月,30102工作面采掘距本次野外数据采集时间仅4个月。三维地震工作的同时,在沿30101工作面走向中心位置布置一条倾向剖面线A,在该剖面线上布置4个钻孔(H1,H2,H3,H4号钻孔)。对照钻孔(H1,H2)布置在30101工作面加风巷东未采动区;为了查明导水裂隙带沿工作面倾向的高度变化特征,在距30101加风巷外边界以内25 m处布置1个H3号钻孔;为了查明导水裂隙带最大高度,在30101工作面中心位置布置1个H4号钻孔。所有钻孔均进行抽水试验、冲洗液漏失量观测。
图1 研究区位置Fig.1 Position of study area
图2为A勘探线地震时间剖面。经人工合成记录标定,Tq为一组可连续追踪反射波,该波组标定为基岩面反射波,在地震时间剖面中用黄线表示。Tz为一组较连续反射波,该波组标定为直罗底反射波,在地震时间剖面中用橙线表示。T3为一组连续性好,能量较强的反射波,该波组标定为3号煤层底板反射波,在地震时间剖面中用蓝线表示。
图2 地震时间剖面(变密度显示)Fig.2 Seismic time profile(variable density display)
研究区地质构造简单,采动之前地震时间剖面主要由Tq,Tz,T3等几组反射波构成,剖面以近水平同相轴为主。煤层采空后,其上覆岩石失去支撑而导致平衡破坏,应力重新分布,使上覆岩体产生变形、位移和破坏,形成垮落带、断裂带和弯曲变形带,采空区内往往会表现出有别于正常地层表征的反射波场特征。
垮落带,岩层破碎、地层完全性被破坏,该带内岩层的波组动力学性质发生变化,波组错断、波形凌乱。垮落带中T3反射波同相轴消失并伴有绕射波和各种散射。
断裂带,地层破坏程度小于垮落带,地层完整性未被完全破坏,地层中发育的大小不一裂隙使地层连续性和地震反射波能量(如Tz波)均受到影响。地震波在隙裂带中传播后,各频率成分的能量分布将发生变化,主要表现为反射波频率降低,为研究断裂带发育提供了依据。断裂带上部约110 ms位置处反射波同相轴断续出现,虽有一定的能量,但连续性差,可以用来判断断裂带发育最大高度。30102工作面内断裂带反射波呈杂乱状,该位置裂隙发育较为剧烈,对反射波的吸收也强烈;同一位置30101工作面有部分反射波发育,说明经历一定时间后,部分裂隙发生了闭合,减少了对反射波的吸收。
变形带属于韧性变形,波组连续性较好和能量较强(如Tq波)。从图2可以看出,变形带底部相对煤柱位置出现反射波同相轴弯曲特征。
经以上分析,可以对采空区三带进行大致的划分,如图2所示。但单从剖面特征来判断导水裂隙带发育高度精度有限,不能满足生产的要求。
正常地层、30102工作面(开采4个月)、30101工作面(开采18个月)对应地层的反射波组特征是不一样的,其频率特征如图3所示。
图3 采后不同时间断裂带地震响应Fig.3 Seismic response of fissure zone at different mining time
从图3可以看出,正常地层反射波主频62 Hz;当煤层采空后,断裂带发育,对反射波高频成分吸收严重,开采4个月后反射波的主频变为26 Hz;采动18个月后,随着地层沉降、压实,部分裂隙闭合,减少了对反射波的吸收,反射波主频有所提高,变为38 Hz。
研究区测井资料受断裂带影响,数据采集并不理想,导水裂隙带发育高度主要靠钻孔冲洗液漏失量观测确定,此为判断导水裂隙带发育高度的直接依据。从三维地震反射波特征来看,单从剖面特征来判断导水裂隙带发育高度较为困难。笔者在钻孔、抽水试验、冲洗液漏失量观测等数据的基础上选择对裂隙变化敏感的地震属性,以冲洗液漏失量观测数据为监督样本,对三维地震频谱分解、相干、蚂蚁追踪等属性数据进行有监督的深度神经网络学习,得到基于三维地震的精细裂隙模型对导水裂隙带发育情况进行预测,预测流程如图4所示。
图4 导水裂隙带预测流程Fig.4 Flow chart of water diversion fracture zone prediction
目前,常用于裂隙预测的地震属性归纳起来可分为:波形相似类,包括相干体、边缘检测和方差体等;几何特征类,包括曲率属性、倾角属性等;吸收衰减类,包括振幅属性、频率属性和频谱属性等。
为了从各种地震属性中优选出与漏失量观测数据密切相关的地震属性,提取了井旁地震道的15种地震属性数据。利用式(1)求取地震属性S与漏失量q之间相关系数,选择与漏失量相关性较好的属性。
(1)
表1为地震属性与漏失量相关系数。从表1可见地震属性与漏失量相关系数最大者为相干属性,相关系数达0.54,而最小者为方位倾角,相关系数仅为0.04。
表1 地震属性与漏失量相关系数Table 1 Correlation coefficient between seismic attributes and leakage
根据不同属性的特点及与漏失量的相关性,提取相干体属性、地层倾角属性、最大曲率属性、谱分解30 Hz数据体、蚂蚁体、瞬时振幅属性、瞬时频率进行裂隙预测,地震属性如图5所示,从图5可以看出,不同地震属性对裂隙所表现的敏感程度是不同的。
图5 地震属性剖面Fig.5 Seismic attribute profile
由表1可以看出,单一地震属性与钻孔冲洗液漏失量的相关性较低,为了更好地研究裂隙,本次采用深度前馈神经网络(DFNN)技术进行属性融合。
深度前馈神经网络模型是一种典型的深度学习模型[15-17]。其目的是当某个近似函数f信息经过x的函数,定义f期间的计算过程,最终到达输出y。模型的输出与模型自身没有反馈连接,如图6所示。输入训练样本为7维列向量T[x1,x2,x3,x4,x5,x6,x7],其表示为地震属性。训练过程中,对输入训练数据和经过神经元输出的结果,采用Xavier方法初始化,随机产生对应的权重系数w和偏置项b。
图6 DFNN的工作流程示意Fig.6 DFNN schematic workflow
神经元由一个非线性Sigmoid逻辑函数构成,其表达形式为
(2)
在向前传播的过程中,数据以加权平均和的形式作为逻辑函数输入到第1个隐藏层中的神经元中。第1个隐藏层中7个神经元的输入数据分别为
(3)
经过非线性逻辑函数处理,会得到3个输出结果,分别为f1(z1),f2(z2),f3(z3),它们加权平均的结果将作为第2个隐藏层中神经元的输入。第2个隐藏层中2个神经元的输入数据分别为
同理,可以得到第2隐藏层中每个神经元的输出结果,分别为f4(z4),f5(z5),它们的加权平均之和会作为输出层神经元的输入,即
z6=w(4,6)f4(z4)+w(5,6)f5(z5)+b6
经过输出层之后,将得到最终的预测结果f6(z6)。由于采用的代价函数是非线性,通过其求解的方程不能实现期望结果。所以,必须利用数值优化进行求解。本文采用梯度下降法进行优化求解。梯度下降法从初始点采用一阶线性逼近,沿着负梯度方向移动,后回到原函数,反复迭代至收敛[18]。利用梯度下降法对w和b进行优化,为了满足梯度为0的一阶最优条件,需要使目标函数为凸函数,但是,目标函数实际上为非线性函数,不属于凸函数。因此,在计算中利用正则项为L1范数来缓解。当达到1 000次迭代时,终止数值优化。
从图7可以看出,通过深度学习地震属性融合,融合后的地震属性与漏失量的相关系数r达到95%,预测结果如图8所示。
图7 相似系数Fig.7 Similarity coefficient
图8 多属性融合剖面(过孔垂直工作面剖面)Fig.8 Multi-attribute fusion profile
通过单属性和多属性融合剖面比较(图5,8)可以看出,多属性融合是单一属性的综合,集合了各个属性的优点,如瞬时频率、相干体对大裂隙反映较好,瞬时振幅对断裂带整体形态反映较好,蚂蚁追踪对细微裂隙反映较好,30 Hz谱分解数据对保留煤柱反映较好。多属性融合数据通过多种属性的融合降低了地震解释的多解性。深度学习结果与多元回归、概率神经网络相比(图9),深度学习预测结果与漏失量数据吻合度最高,预测结果也符合地质规律,而概率神经网络计算结果在垮落带附近发生连续跳跃(如图9蓝色方框内),不符合力学特征。
多属性融合体值域与漏失量相同,值在0~4,值越大代表裂隙越发育。从图9可以看出,从垂直方向看裂隙发育呈高角度状,断裂带自上而下分为紫色部分(值在0~1)、蓝色部分(值在1~2)、绿色部分(值在2~4)、黄、红色部分(值在3~4),代表裂隙从弱到强,越靠近煤层,裂隙越发育。
图9 H3孔多方法裂隙预测比较Fig.9 Comparison of multi-method fracture prediction
通过多属性融合一方面大大提高了三维地震的成像精度,消除多解性;另一方面使得地震数据转换成冲洗液漏失量观测数据,漏失量观测数据大小和裂缝发育强度是呈正相关的[5,19-20],根据这一原理可以实现导水裂隙带的半定量预测。多属性融合结合了不同地震属性的优点,其解释结果较为客观。通过与漏失量数据对比分析认为,融合数据中数值在0~0.5为裂隙不发育区,划分为变形带和正常区;数值在0.5~3.5划分为断裂带;数值>3.5划分为垮落带。
利用深度学习融合的数据进行断裂带、垮落带解释,结合钻孔拟合时深转换公式:
Hd=0.006 9t2-3.784 8t+464.1
(6)
将地震时间域数据转为深度域数据。经计算,30101工作面断裂带发育最大高度为120 m,30102工作面断裂带发育最大高度为133 m。对单个工作面来说,工作面中部断裂带发育高度达到最大,向两侧采空边界处断裂带发育高度逐渐减小,保护煤柱附近断裂带发育高度约为70 m。从图10可以看出,随着埋深增加,裂隙发育程度不断加大,在垮落带附近裂隙发育程度达到最大。单个工作面断裂带发育形态为“拱形”,裂隙特征表现为密集的网状分布,同一深度下30102工作面裂隙发育程度大于30101工作面。
图10 裂隙发育平面特征Fig.10 Plane characteristics of fissure development
通过对融合后的三维地震数据进行裂隙自动提取,共解释出裂隙200多个,如图11(a)所示。裂隙饼状图如图11(b)所示,图11中半径大小表示倾角,径线方向表示方位角。由于地震分辨率、信噪比和裂隙组合的原因,追踪出的裂隙数量不代表实际裂隙的数量,但可反映区内裂隙的整体发育程度。
图11 提取裂隙及裂隙统计特征Fig.11 Predict cracks and their statistical characteristics
从图11(b)可以看出裂隙走向主要以平行和垂直巷道两个方向为主;裂隙倾角多为高角度,一般大于70°。图11(c)为采动裂隙倾角分布图,倾角小于30°的裂隙占8%,倾角为30°~50°的裂隙占12%,倾角为50°~60°的裂隙占9%,倾角为60°~70°的裂隙占17%,倾角为70°~80°的裂隙占18%,倾角为80°~90°的裂隙占36%。煤层开采后在岩层中形成两类裂隙,一类为离层裂隙,另一类为竖向破裂裂隙[21]。限于裂隙自动提取精度,本文将倾角小于30°的裂隙划分成离层裂隙。可以看出采动裂隙以高角度甚至垂直岩层层面的裂隙为主。
4.2.1断裂带发育高度预测评价
将过井地震融合数据与钻探数据进行比较,如图12所示。
图12 钻孔综合成果图(部分层段)Fig.12 Comprehensive results of drilling(part of the layer)
H1号钻孔:在基岩段单位时间冲洗液消耗量变化为0.062~0.180 L/s,平均值0.329 L/s。在整个观测过程中,冲洗液消耗量没有随着钻孔孔深增加而增加,冲洗液循环正常,没有出现中断或全部漏失的现象。过井地震融合数据数值在0~0.3。根据3.3节判断依据认为没有采动裂隙发育,与钻探判断结果一致。
H3号钻孔:孔深135.16 m时冲洗液消耗量突然增大,从135.16 m以浅的0.062 L/s,增大至2.622 L/s,增大了42倍,比对比钻孔同层位消耗量(0.077 1)增大了34倍,钻孔水位呈缓慢下降趋势;结合岩芯编录判定导水裂隙带顶界位置孔深为135.16 m处。过井地震融合数据在深度138 m时数值增大,达到0.5,根据3.3节判断依据认为导水裂隙带顶界位置为138 m处。
H4号钻孔:孔深为125.20 m时钻孔冲洗液突然漏失,循环中断,继续钻进2.0 m时冲洗液又开始循环,并且漏失量呈忽大忽小的宽幅度震荡。结合岩芯编录判定导水断裂带顶界位置孔深为125.20 m处。过井地震融合数据在深度129 m时数值增大,达到0.5,根据3.3节判断依据认为导水裂隙带顶界位置为129 m处。
4.2.2断裂带倾角预测评价
本区进行了部分电视测井,将电视测井与过井地震融合数据进行对比,如图13,14所示。从图13可以看出,H3孔电视测井中在134.50 m开始出现垂向裂隙,延伸较短,随着孔深增大,裂隙密度增加,且裂隙延伸增大;144.10~146.10 m为1条长2.0 m的斜向破坏裂隙,裂隙倾角约87°。截取相同深度地震融合数据可以看出,在135~145 m间有裂隙,裂隙倾角约80°与电视测井结论接近。
图13 H3号钻孔电视测井图像Fig.13 H3 borehole TV logging image
从图14可以看出,H4孔电视测井中在128.60 m开始出现斜向裂隙,延伸较短,随着孔深增大,裂隙密度增加,以斜向裂隙为主,且裂隙延伸较短,裂隙倾角约45°。截取相同深度地震融合数据可以看出,在125~130 m间有裂隙,裂隙倾角20°~60°角度较小与电视测井结论接近。
图14 H4号钻孔电视测井图像Fig.14 H4 borehole TV logging image
从以上分析可以看出基于深度学习的地震属性融合技术导水裂隙带预测效果较好,精度较高。
研究区煤层稳定,工作面布设及开采方式相同,因此可以将30102工作面采动情况看作是30101工作面14个月前的状态。本次研究采用类比法将采集的地震数据按照采动时间进行分类,研究裂隙生成、发育、闭合情况。
为了比较2个工作面裂隙发育情况,以深度学习后融合的地震数据为基础,分别提取30101,30102工作面范围内裂隙,并统计其规律如图15所示。
图15 不同时间开采后裂隙饼状图Fig.15 Pie chart of cracks after mining at different times
从图15可以看出,煤层开采18个月后(30101工作面)裂隙比煤层开采4个月(30102工作面)有所减少,经比较认为减少了约21%,其中低角度离层裂隙(倾角0°~30°)减少了50%。
随着采后时间推移,除裂隙数量有所减少之外,断裂带发育最大高度也有所降低,经计算认为30101工作面断裂带最大高度比30102工作面断裂带最大高度减小了约13 m。通过2个工作面比较可以看出距煤层顶板较远岩体中裂隙发育强度较小,裂隙闭合较好;距煤层顶板较近岩体中裂隙发育强度大,裂隙未能完全闭合。工作面中部裂隙发育强度变化较大,裂隙闭合较好;工作面边缘裂隙发育强度变化较小,裂隙闭合较差。
(1)煤层采动过程中其上覆岩层变形与破坏受到多种因素影响,其中煤层采高、采厚、采速、覆岩类型及埋深等参数对导水裂隙带高度发育起到主导作用。岩层破坏在地震剖面中根据其振幅、频率的时空变化可以分辨出结构破坏及裂隙发育。探测应用表明,基于深度学习的地震属性融合技术对煤层采动引起的岩层破坏规律探测具有针对性,其精度能满足生产需要,费用低,是一种有效探测技术,具有推广应用价值。
(2)研究表明30101工作面断裂带发育最大高度为120 m,30102工作面断裂带发育最大高度为133 m,以高角度裂隙为主,主要沿垂直、平行工作面方向发育。导水裂隙带的裂隙发育是先增大后降低,断裂带上部裂隙闭合较好,断裂带下部和工作面边缘裂隙闭合较差,采动后18个月裂隙比采动后4个月减少了21%,离层裂隙减少了50%。
(3)本次研究受原始资料所限,采用漏失量数据进行深度学习,今后研究中可采用其他更好表征裂隙特征的数据如裂隙密度数据进行深度学习,可以提高预测精度。