煤层陷落柱多波地震数值模拟与成像研究

我国具有丰富的煤炭资源，但煤层构造比较复杂。煤炭的形成与分布受到各种地质作用的影响及控制，因此正确了解煤系地层的各种地质因素，才能保障煤矿生产的顺利进行。影响煤矿生产的地质因素很多，如地质构造、煤层的厚度变化、岩浆侵入体、陷落柱、煤层顶底板岩性、地层压力、地层水及瓦斯含量等。其中，陷落柱是我国华北、华南等地区煤田开采中的一种典型地质体，其形成主要与构造破碎带、裂隙、地下水的活动和岩溶等因素有关(张茂林等，2007)。在含煤地层的下部由于地下水的常年活动，导致周围可溶性岩石溶蚀形成空洞，随着地下水不断的溶蚀，溶洞范围越来越大，在上覆岩层重力以及地质构造应力的长期作用下，溶洞发生坍塌，上部含煤岩层随之塌陷，形成陷落柱(曹志勇等，2012)。陷落柱柱体充填物成分复杂、松散, 密度小, 速度低, 成层性差, 而陷落柱附近及顶部围岩多为砂岩、泥岩或煤层, 其沉积稳定,因此与围岩相比，在速度、密度以及地层的连续性、产状、岩性等方面均有很大的差别，导致正常反射波组不连续或能量变弱(尹尚先等，2004)。陷落柱对煤炭与煤层气资源的开采具有重大危害，不仅破坏煤层的完整性和连续性，影响巷道的掘进和煤层的开采，而且有可能成为特殊的导水通道，给煤矿高效安全生产带来巨大的隐患(尹奇峰等，2012)。因此，准确的探测煤系地层中陷落柱，对煤炭生产中工作面的合理布置、巷道维护和及时修正，提高开采效率，保障煤矿安全生产具有重要的意义。

近年来，煤田地震勘探方法在煤田勘探上取得了良好的效果(Gochioco, 2000; Walton et al., 2000; Zuo Jianping et al., 2009)，已经成为查明煤矿生产中断层、褶曲、陷落柱、采空区、冲刷带等地质信息的主要技术手段，在地震地质条件较好的地区可以查明30 m左右的陷落柱，对解决煤矿生产中的诸多地质问题提供了重要的依据。但随着大规模机械化采煤的展开，需要解决的地质任务越来越精细，对煤层探测精度的要求越来越高。由于陷落柱等构造在分布上具有复杂性以及其对地震响应的特殊性(李飞等，2009), 常规的地震勘探技术已无法满足现阶段煤田精细勘探的需求。多分量地震勘探能够获得纵波与转换波地震资料，提供更为丰富的地震波场信息。此外，由于横波(S波)速度小于纵波(P波)，在浅层S波的波长相对很小，而地震垂直分辨率是波长的1/4，因此转换PS波在浅层具有更高的分辨率，对埋深较浅的煤系地层中小幅度构造能够获得更好的成像，更清晰的刻画出构造形态及内部变形等特征(Stewart et al., 2002, 2003)。虽然在实际地震勘探中由于S波的衰减较快，转换PS波在深层的分辨率通常低于PP波，但充分利用多波地震资料不仅能更好的揭示地下构造信息，而且能有效的反映裂缝等各向异性信息以及流体特性，有利于煤田勘探精度的提高。

目前，多波地震勘探技术正处于迅速发展时期，在多分量地震技术理论研究以及多分量地震数据采集、处理和解释反演等关键技术上取得了一系列的成果(Gu Bingluo et al., 2015, 2017, 2018；Han Jianguang et al., 2014, 2016, 2017；Lu Jun et al., 2015, 2016, 2018)。多波地震勘探技术的发展为探测煤田陷落柱提供了新的思路，本文将多波地震技术引入到煤层陷落柱研究中，利用数值方法对煤层陷落柱进行多波地震勘探模拟研究，采用弹性波有限差分方法进行正演模拟，得到多分量地震记录，然后分别采用适用于PP波以及转换PS波的叠前深度偏移方法对其进行偏移成像，对PP波和PS波的成像结果进行对比分析。通过数值模拟研究验证多波地震技术对于煤层陷落柱勘探的可行性和有效性，为进一步将多波地震勘探技术应用于实际煤层陷落柱的探测提供理论基础。

1 原理与方法

1.1 弹性波正演模拟基本原理

弹性波动方程是多波地震波场正演模拟的理论基础，地震波数值计算本质上是利用数值方法求解弹性波动方程。二维各向同性介质中的弹性波动方程为

式中，u和w分别表示位移的水平分量和垂直分量，ρ为介质密度(g/cm3)，λ和μ为Lame常数。

有限差分法是最重要的一种地震波数值模拟方法(Moczo et al., 2007)，利用差分算子近似相应的微分算子，建立差分方程来求波动方程的数值解。差分系数是影响有限差分地震波场数值模拟精度的重要因素，本文采用规则网格有限差分进行数值模拟，这里直接给出二阶偏导数2N阶有限差分系数:

利用如上系数可以得到规则网格二阶导数中心差分算子:

地震波数值模拟中震源项为初始条件，波场在震源驱动下，通过波动方程向外传播。震源作用具有一定的时间延续和空间范围，我们利用震源子波函数f(t)来表征时间延续，通过一个随着到震源点距离呈指数衰减的函数g(x,y,z)来表示空间作用范围，则震源作用力可以表示为:

S(x,y,z,t)=g(x,y,z) · f(t)

震源有多种形式，其中爆炸震源的震源力均匀作用于震源点周围一个球腔壁上，方向与球面垂直，在各项同性介质中，只产生无旋场。本文采用爆炸震源，二维情况下弹性波方程爆炸震源的加载方式如下：

此外，地震波场正演模拟的另外一个关键问题是边界条件，目前已经发展了衰减吸收边界(Cerjan et al., 1985；Sochacki et al., 1987)、旁轴近似吸收边界(Clayton et al., 1977; Higdon, 1987, 1991)以及完美匹配层(PML)吸收边界条件等多种类型，其中有限差分数值模拟中应用最为广泛的是PML吸收边界条件。PML吸收边界通过在计算边界处引入一个吸收层，使波场模拟时在模型边界上不产生反射。

为了更好地说明PML边界条件，这里对空间坐标p(p=x或p=z)引入复数坐标系：

其中，

为复数空间坐标，i为虚数单位，dp(s)为边界衰减系数，它是随着坐标p衰减的实函数，ω为角频率。关于p和

的一、二阶偏导数满足如下关系：

基于上述复数坐标系，弹性波方程可以变化为(以左边界为例)：

其中，

和

分别是频率域的水平和垂直位移波场；dx为x方向的吸收系数，其表达式如下：

式中，x表示到计算域的距离，

为dx在x方向的导数。

同理，可获得垂直方向的包含吸收边界的弹性波方程。将包含吸收边界的弹性波方程变换至时间—空间域，利用有限差分即可实现弹性波数值模拟。

1.2 多波叠前深度偏移成像方法

在多分量地震勘探中，X分量与Z分量均可以接收到纵波与横波的能量，我们首先对原始的多分量地震数据进行叠前分离，获得相应的PP波和转换PS波地震记录，然后分别对PP波和转换PS波进行叠前深度偏移成像。本文采用高斯束偏移成像方法进行测试试验，高斯束偏移是一种优秀的偏移方法，同时兼具波动方程偏移方法的成像精度以及Kirchhoff偏移方法灵活、高效的优点，并且能有效的解决多值走时问题。在二维各向同性介质中，PP波和PS波共炮点道集高斯束叠前深度偏移成像见Han Jianguang 等(2013)。

2 数值模拟分析

煤田地质中典型的陷落柱按平面形态可分为圆锥状、筒状、斜塔状、漏斗状以及不规则状(林建东等，2012)。本文以最为常见的圆锥状陷落柱为例，结合实际陷落柱的发育特征，建立一个基本的陷落柱模型，陷落柱底部直径50 m，高110 m(如图1)。设计模型长为1600 m，深800 m，纵横向网格大小均为5 m。模型共4层，第1层为第四系沉积物，第2层为砂岩，第3层为煤层，第4层为灰岩，具体参数如表1所示。下面对上述陷落柱模型进行多波地震数值计算试验，采用弹性波有限差分法进行正演模拟，震源为50 Hz的Ricker子波，全排列接收，共模拟51炮，炮间距为30 m，每炮321道接收，道间距为5 m，采样时间1.2 s，采样间隔为1 ms。图2a和图2b分别为第26炮Z分量和X分量地震记录，从图中可以看到两个分量上同时含有PP波和转换PS波信号，而且X分量记录存在极性反转现象。对多分量地震记录进行波场分离得到纯的PP波和PS波炮集记录，分别采用相应的高斯束叠前深度偏移方法进行偏移成像，偏移过程中校正PS波极性反转现象且消除直达波的影响，得到最终的PP波和PS波偏移成像结果，如图3a和图3b所示。从图中可以看到，对于两种波地震数据均取得了准确的偏移结果，两个偏移剖面上对应的反射界面均归位到准确的层位深度，而且对比两个剖面可以看到PS波的成像结果具有更高的分辨率。此外，由于模型中陷落柱的尺寸较小，且陷落柱的柱体为高陡构造，在短排列采集条件下无法获得柱体成像，但从图中可以看到煤层在陷落柱位置处出现明显的中断和缺失(如箭头所示)，可以直观地辨别陷落柱的真实位置。

为了进一步研究多波地震技术对于更小尺寸陷落柱构造的勘探效果，将上述模型中陷落柱尺寸修改为底部直径30 m，高70 m，模型其他参数不变，如图4所示。采用相同的正演参数对小尺寸陷落柱模型进行弹性波地震波场模拟，图5a和图5b为第26炮Z分量和X分量单炮地震记录，图6a和图6b分别为波场分离后PP波和PS波地震记录的高斯束叠前深度偏移成像结果，偏移过程中校正了PS波极性反转现象且消除了直达波的影响，可以看到两个剖面同样均获得了准确的偏移成像结果，反射界面归位准确，煤层在陷落柱处存在清晰可见的中断(如箭头所示)，可以有效的识别出陷落柱。

3 结论

通过利用数值方法对煤层陷落柱进行多波地震勘探模拟研究，获得以下认识：

(1)多波地震勘探技术是一种有效的煤层陷落柱探测方法，联合PP波和PS波地震数据可以提供更多的地下介质信息，能够对煤层陷落柱取得更好的勘探效果。

(2)通过对基于实际陷落柱发育特征建立的两个小尺寸陷落柱模型进行多波地震数值试算，可以看到PP波和PS波地震数据均获得了准确的偏移成像，煤层在陷落柱位置处出现明显的中断和缺失，可以有效的识别陷落柱的真实位置，且PS波成像结果相对于PP波具有更高的分辨率，充分利用多波地震资料有利于查明煤层陷落柱构造。

(3)本文通过数值模拟研究验证了多波地震技术对于煤层陷落柱探测的可行性和有效性，为多波地震技术应用于实际煤层陷落柱勘探提供了理论基础。

(The literature whose publishing year followed by a “&” is in Chinese with English abstract; The literature whose publishing year followed by a “#” is in Chinese without English abstract)

曹志勇，王伟，王赟．2012．煤层陷落柱散射波数值模拟与成像．地球物理学报，55(5)：1749～1756.

李飞，张智，曹志勇，李晶．2009．陷落柱地震波场特征分析．地球物理学进展，24(3)：886～892．

林建东，杨振，颜中辉，赵显宗．2012．煤矿陷落柱的地震响应特征与识别．地球物理学进展，27(6)：2631～2638．

尹奇峰，潘冬明，于景邨，刘盛东．2012．煤矿陷落柱地震识别技术研究．地球物理学进展，27(5)：2168～2174．

尹尚先，武强，王尚旭．2004．华北煤矿区岩溶陷落柱特征及成因探讨．岩石力学与工程学报，23(1)：120～123．

张茂林，尹尚先． 2007．华北型煤田陷落柱形成过程研究．煤田地质与勘探，35(6)：26～29．

Cao Zhiyong， Wang Wei， Wang Yun. 2012&. Numerical simulation and imaging of scattered wave of sunk pillar in coal seam. Chinese Journal of Geophysics, 55(5): 1749～1756.

Cerjan C, Kosloff D, Kosloff R, and Reshef M. 1985. A nonreflecting boundary condition for discrete acoustic and elastic wave equation. Geophysics, 50(4): 705～708.

Clayton R, and Engquist B. 1977. Absorbing boundary conditions for acoustic and elastic wave equations. Bulletin of the Seismological Society of America, 67(6): 1529～1540.

Gochioco L M. 2000. High-resolution 3-D seismic survey over a coal mine reserve area in the U.S.——A case study. Geophysics, 65: 712～718.

Gu Bingluo， Li Zhenchun， Han Jianguang. 2018. A wavefield-separation-based elastic least-squares reverse time migration. Geophysics, 83(3): S279～S297.

Gu Bingluo， Li Zhenchun， Yang Peng， Xu Wencai， Han Jianguang. 2017. Elastic least-squares reverse time migration with hybrid l1/l2 misfit function. Geophysics, 82(3): S271～S291.

Gu Bingluo， Li Zhiyuan， Ma Xiaona， Liang Guanghe. 2015. Multi-component elastic reverse time migration based on the P and S separating elastic velocity—stress equation. Journal of Applied Geophysics, 112C: 62～78.

Han Jianguang， Wang Yun. 2017. A method of directly extracting multiwave angle-domain common-image gathers. Acta Geophysica, 65: 1009～1017.

Han Jianguang， Wang Yun， Yu Changqing. 2016. Multiwave Gaussian beam prestack depth migration of exploration scale seismic data with complex near-surface effects. Near Surface Geophysics, 14: 307～313.

Han Jianguang， Wang Yun， Lu Jun. 2013. Multi-component Gaussian beam prestack depth migration. Journal of Geophysics and Engineering, 10: 055008.

Han Jianguang， Wang Yun， Han Ning， Xing Zhantao， Lu Jun. 2014. Multiwave velocity analysis based on Gaussian beam prestack depth migration. Applied Geophysics, 11(2): 186～196.

Higdon R. 1987.Numerical absorbing boundary conditions for the wave equation. Mathematics of Computation, 49(179): 65～90.

Higdon R. 1991.Absorbing boundary conditions for elastic waves. Geophysics, 56(2): 231～241.

Li Fei， Zhang Zhi， Cao Zhiyong， Li Jing. 2009&. Analysis on seismic features of subsided column. Progress in Geophysics, 24(3): 886～892.

Lin Jiandong， Yang Zhen， Yan Zhonghui， Zhao Xianzong. 2012&. Seismic response and identification of coal mine collapse columns. Progress in Geophysics, 27(6): 2631～2638.

Lu Jun， Wang Yun， Chen Jingyi， An Ying. 2018. Joint anisotropic amplitude variation with offset inversion of PP and PS seismic data. Geophysics, 83(2): N31～N50.

Lu Jun， Meng Xinghun， Wang Yun， Yang Zhen. 2016. Prediction of coal seam details and mining safety using multicomponent seismic data: A case history from China. Geophysics, 81(5): B149～B165.

Lu Jun， Yang Zhen， Wang Yun， Shi Ying. 2015. Joint PP and PS AVA seismic inversion using exact Zoeppritz equations. Geophysics, 80(5): R239～R250.

Moczo P, Kristek J, Galis M, Pazak P, and Balazovjech M. 2007. The finitedifference and finite-element modeling of seismic wave propagation and earthquake motion. Acta Physica Slovaca, 57(2): 177～406.

Sochacki J, Kubichek R, George J, Fletcher W, and Smithson S. 1987. Absorbing boundary conditions and surface waves. Geophysics, 52(1): 60～71.

Stewart R R, Gaiser J E, Brown R J, and Lawton D C. 2002. Converted wave seismic exploration: Methods. Geophysics, 67: 1348～1363.

Stewart R R, Gaiser J E, Brown R J, and Lawton D C. 2003. Converted wave seismic exploration: Applications. Geophysics, 68: 40～57.

Walton C, Evans B, and Urosevic M. 2000. Imaging coal seam structure using 3-D seismic methods. Exploration Geophysics, 31: 509～514.

Yin Qifeng， Pan Dongming， Yu Jingchun， Liu Shengdong. 2012&. Research on seismic identification technique of coal mine collapse column. Progress in Geophysics, 27(5): 2168～2174.

Yin Shangxian， Wu Qiang， Wang Shangxu. 2004&. Studies on characters and forming mechanism of karstic collapse columns at mine area of north China. Chinese Journal of Rock Mechanics and Engineering, 23(1): 120～123.

Zhang Maolin， Yin Shangxian. 2007&. Forming process of subsided column in coalfields of North China. Coal Geology & Exploration, 35(6): 26～29.

Zuo Jianping， Peng Suping， Li Yongjun， Chen Zhonghui， Xie Heping. 2009. Investigation of Karst collapse based on 3-D seismic technique and DDA method at Xieqiao coal mine, China. International Journal of Coal Geology, 78: 276～287.