三维弹性波全波形反演在小保当煤矿勘探中的应用研究

刘 超，刘 鹏，张慧峰，白 坤，李志勇，潘文勇

(1.陕西小保当矿业有限公司，陕西 神木 719399；2.北京中矿大地地球探测工程技术有限公司，北京 100038；3.中国科学院地质与地球物理研究所 油气资源研究院重点实验室，北京 100029)

为实现煤炭资源高效、安全、智能化开采，提高资源利用率，获取精准的三维地下地质模型尤为重要，是采煤工作面信息透明化管理的重要基础[1-4]。然而，常规地震探测技术无法满足当前煤矿开采的成像精度和分辨率需求，亟需全新的地震反演与成像技术。全波形反演(Full-waveform inversion，即FWI)是近年来兴起的“新一代”高分辨率地震成像技术，通过对弹性波动方程的数值求解，充分利用全波形信息进行反演成像，可构建地下复杂构造的弹性参数模型，实现三维地质模型的精细刻画[5-8]。地震全波形反演技术在油气勘探领域已获得大量成功应用，但在煤炭行业应用较少，国外仅基于2D地质模型做理论数据的FWI，帮助煤矿勘探[9]，国内研究机构只针对二维工作面进行合成数据的FWI拟合[10]，受限于精度要求高，计算成本过大的实际操作问题，至今国内外还没有三维FWI在煤矿上的实际应用研究。而为了精细刻画井下煤矿模型，对陷落柱、断层以及煤岩分界面清晰成像，使用三维FWI算法可同时降低计算成本，并确保反演结果精确度。

本研究，首先需要通过求解三维弹性波动方程进行正演模拟获得模拟数据，在计算梯度时还需要求解伴随波场。常规的地震波正演模拟方法包括有限差分法、有限元法、伪谱法和边界元法等[11-13]，与这些方法相比，谱元法(Spectral-Element Method，即SEM法)，具有伪谱法高精度优点，同时可灵活处理不规则边界，非常适合GPU并行计算。本研究采用SEM法模拟地震波在复杂三维弹性介质中的传播[14]，使用CPU/GPU异构并行技术提高计算效率。

以陕西小保当煤矿为例，采用三维弹性波全波形反演技术，对矿区地下三维速度结构进行高分辨反演成像，预测煤层顶底界面，实现采煤工作面的透明化，达到整个矿区的信息透明化管理和安全高效采煤。根据煤矿巷道勘探的特点，发展了多尺度、多目标函数和多波反演方法与技术。由于地震数据缺失低频信息，FWI具有非常强的非线性，容易陷入局部最小值，即周期跳跃问题。采用多尺度反演策略[15]，首先使用低频信息反演速度的长波长部分，然后使用高频信息反演速度结构短波长部分[15]；传统的波形差目标函数非常容易受周期跳跃问题的影响，发展多目标函数反演策略，依次使用互相关旅行时[16]、包络[17]和波形差目标函数进行反演，可有效克服周期跳跃问题；在煤矿巷道采集的地震数据中，包含P波、S波和槽波，可使用P波和S波反演背景速度结构，槽波在煤层的顶底界面经过多次全反射，在煤层中具有较强的穿透性，对煤层顶底界面成像有重要作用。采用多波反演策略，可更好地构建煤层的三维空间结构。

为实现三维弹性波FWI在实际地震数据处理中的成功应用，地震数据预处理至关重要，对小保当矿区采集的实际地震数据进行了一些机预处理，具体包括：重新采样、去除坏道、带通滤波和振幅归一化等。使用已知的测井数据(包括煤层顶底界面位置和煤层厚度)作为FWI的先验信息进行正则化约束，可减弱反演的多解性。从每炮地震数据的直达波中提取子波函数，作为正演模拟的震源函数。在反演中，对梯度进行异常值去除和三维高斯平滑，降低反演的病态性。最终，获得高分辨率的三维弹性速度结构模型，预测煤层的三维空间展布和顶底界面，使用测井数据验证反演结果的准确性。

1 三维弹性波全波形反演基本理论

在FWI中，传统波形差目标函数测量了观测数据di与模拟数据ui之间的直接差值，通过求取目标函数的最小值，使用最优化算法对模型参数进行迭代反演，最终使得模拟数据与观测数据达到最佳匹配，反演的模型与地下真实模型最大程度接近。波形差目标函数可以写为[18]：

式中，i，j，k，l∈[1，2，3]或[x，y，z]；ui为位移场；ρ为密度式中，fi为震源式中，cijkl和ekl分别为弹性刚度张量和应变张量。根据Voigt标记法，cijkl可简化为CIJ(I，J∈[1，2，3，4，5，6])。在三维各向同性弹性介质中，独立的弹性系数简化为C33和C44，与速度的关系可表示为：

式中，VP和VS分别为P波和S波速度。基于谱元法模拟地震波在复杂三维弹性介质中的传播[14]，与传统有限差分法相比，精度更高[18，19]。使用国家“天河”超算平台的GPU节点(显卡型号为NVIDIA A100，40 GB显存，6192 CUDA核)加速计算，与常规多线程CPU节点(18核)相比，计算效率提高50～60倍。

FWI使用局部最优化算法(如最速下降法和牛顿法等)求取目标函数的最小值，需要计算目标函数对模型参数的一阶偏导，即梯度(或敏感核)。为避免大型雅克比矩阵的直接求取，基于伴随状态法，弹性速度的敏感核可以通过正演模拟波场和伴随波场的互相关计算[20-22]：

mn+1=mn+αnΔmn(6)

式中，αn为第n次迭代的步长，可使用线性搜索方法计算获得；Δmn为搜索方向，使用拟牛顿L-BFGS优化算法计算[23，24]，并实现迭代反演。本研究使用mask函数去除梯度异常值，通过三维高斯函数对梯度进行平滑，提高反演的稳定性。

由于地震波形与模型参数之间具有非常强的非线性[25-27]，使用传统波形差目标函数容易陷入局部最小值，得到不准确的反演结果。针对该问题，发展了多尺度、多目标函数和多波反演以及正则化约束等方法：

1)在多尺度反演方法中，对数据进行带通滤波，使用低频段数据构建模型的长波长部分，拓展频率范围，基于低频反演结果，使用高频信息构建短波长部分，获得最终反演结果。

2)采用多目标函数反演策略，依次使用互相关旅行时、包络和波形差目标函数进行反演。旅行时目标函数为：

式中，Δti为旅行时差。包络目标函数为：

3)发展了煤矿巷道多波反演方法，首先使用P波和S波分别反演目标区域的VP和VS结构。由于槽波在煤层顶底界面经过多次全反射，引入槽波使用整个炮集记录进行反演，更加有助于对煤层顶底界面进行成像。

4)基于已有测井数据(包括煤层顶底界面位置和煤层厚度)构建先验模型mprior，在FWI目标函数中加入惩罚项，建立正则化目标函数：

式中，β为权重算子；W为正则化参数；χ(m)为常规目标函数。使用该目标函数可压制不符合先验信息的解，降低反演的多解性。

2 小保当矿区概况及地质背景

小保当矿区位于陕西省神木市西南部，侏罗纪煤田榆神矿区西部鄂尔多斯台地向斜东翼，面积约220 km2。小保当一号煤矿开采2-2煤层，112201和112202为倾角较小的单斜构造，三维地震未发现较大断层和明显的褶皱构造。煤层存在宽缓的波状起伏，其轴向大致平行于勘探线方向，主要表现为西部两翼波长较短，东部波长较长，由西向东呈过渡状态，在平面上形成中部为喇叭状凹型地势。此构造形态在各个煤层底板等高线中都有表现，说明该构造形态与含煤地层沉积基底有关。

FWI的目标区域为112204工作面，位于井田西南侧，长约650 m，宽约350 m，如图1所示，图1中红色和蓝色圆球分别为震源和接收器。煤层赋存于延安组第四段顶部，埋深为305～387 m，煤层底板标高924～990 m，煤厚6.30～7.24 m，平均煤厚6.50 m。煤层由北向南逐渐变厚，地质构造简单，存在宽缓的波状起伏。为实现采煤工作面的地质信息透明化，达到整个矿区的透明化管理和安全高效采煤，需通过三维弹性波全波形反演，构建目标地层的精细三维地质结构模型。

图1 FWI工区平面

3 实际地震数据应用

3.1 实际地震数据分析与预处理

首先对采集的实际地震数据进行了分析和预处理。巷道地震数据包括三分量透射接收数据和反射接收数据，震源类型为炸药，共248炮，采样间隔为0.25 ms，最大记录时间为2 s，采样点数为8000。在透射接收的原始地震数据中(如图2所示)，地震波形(包括P波、S波和槽波)连续且清晰，信噪比较高，但也存在一定的随机噪音和一定数量的坏道。图3为原始地震数据的振幅谱，可看到有效信号的频率范围为40 Hz至200 Hz。由于反射接收的原始地震数据存在非常强的高频和低频噪音，有效反射信号非常弱，几乎无法识别。因此，只使用透射接收的地震数据进行反演和成像。对透射接收的地震数据进行了预处理，包括去坏道、滤波、去噪和能量均衡，如图4所示。可以看到，x和z分量质量较高，y分量数据质量较差。总体而言，处理后的数据集可作为FWI的输入，并达到预期反演效果。

图2 工作面运输巷透射接收的原始地震数据

图3 工作面运输巷透射接收数据频谱

图4 处理后的工作面运输巷透射接收数据

为实现有效反演，设计了专门针对巷道全波形反演的地震数据预处理流程，首先根据正演模拟得到相应的采样率，地震数据再按此采样率进行重新采样，选取合适的时间窗口分别提取P波、S波和槽波，去除数据中的坏道，通过带通滤波控制反演的频率范围，并使用能量均衡将振幅归一化。具体数据处理流程为：输入数据—重新采样——时间窗口选取—去除坏道—带通滤波—能量均衡—输出数据。经数据处理后获得的直达P波数据如图5所示。在迭代反演中，需要对地震进行多次处理，为提高效率，开发了地震数据批量处理脚本程序，并与正演模拟和反演程序融合，形成了一套高效的反演算法流程。

图5 处理后的地震直达P波数据

3.2 三维初始模型和网格剖分

根据勘探地区的实际情况以及震源和接收器的三维空间绝对位置，经过坐标旋转和相对坐标转换，建立了正演模拟和反演所需的相对坐标系统。在工作面运输巷和工作面回风巷中震源和接收器的分布如图6所示。该模型X方向长度为650 m，Y方向长度为350 m，Z方向的深度为-30 m(对应海拔为915 m至945 m)。

图6 旋转后的相对坐标系统以及观测系统分布

FWI需要初始弹性参数模型作为输入，然后进行迭代更新。根据地震数据中P波和S波的旅行时，获得了初始模型的P波速度为4000 m/s，S波速度为1750 m/s，密度为2465 kg/m3，作为FWI的初始输入模型。在X、Y和Z方向的网格数量分别为120，72和60，其中垂直Z向的网格大小为0.5 m。

3.3 三维弹性速度反演

使用全波形反演算法与流程以及反演参数对小保当目标区域的三维速度模型进行构建。首先，在频率段[30 Hz，60 Hz]内，使用旅行时目标函数，对三维P波和S波速度结构进行反演。由于旅行时与速度参数的非线性关系较低，可避免周期跳跃现象。将P波和S波分开反演，可减弱由多参数串扰导致的噪音。随后，在该频率段内分别使用包络和波形差目标函数，进行反演，可以进一步提高反演的分辨率。最后使用P波、S波和槽波对速度结构进行反演。加入槽波信息对煤层的顶底界面的反演更加有效。根据多尺度反演策略，继续在频率段[30 Hz，90 Hz]、[30 Hz，120 Hz]和[30 Hz，150 Hz]，重复使用同样的反演流程与算法，从低频到高频逐步构建小保当矿区地下的速度结构。低频信息可反演速度的长波长背景结构，高频信息可构建地层的高频短波长精细结构。

最终反演得到三维P波速度结构模型如图7所示，P波速度结构切片如图8所示。最终反演得到的三维S波速度模型如图9所示，S波速度结构切片如图10所示。从速度反演结果可以看出，煤层的三维空间结构和顶底界面可以清晰成像，地下煤层分布可清晰识别，总体较平缓，无复杂构造(如断层和褶皱等)。而基于传统槽波的探测方法，只能给出工作面内的二维能量成像图或煤层厚度的等值线图[28，29]，无法对煤层的三维空间变化进行成像和预测。由此说明，全波形反演对地下三维空间结构的成像结果可达到煤层透明工作面的实际需求。

图7 反演得到的三维VP速度模型

图8 VP速度模型切片

图9 反演得到的三维VS速度模型

图10 VS速度模型切片

3.4 煤层顶底界面空间分布预测

基于三维全波形反演得到的三维弹性速度结构反演结果，对煤层的顶底界面的空间分布进行预测，结果如图11所示。从图11中可看出，煤层的总体厚度约为6～7 m，空间变化较为平缓，无明显的特殊构造。为验证预测结果的准确性，将在测井数据中获得煤层顶底板界面位置与全波形反演预测结果进行了对比，见表1。测井位于工作面运输巷内，随机选取了5口测井，在X方向的位置分别为100，230，450，575和625 m。可见，全波形反演结果获得煤层顶板误差较小，而地板误差较大，预测精度达到了煤层透明工作面的实际需求，为高效和安全的煤矿开采和生产提供技术支撑和重要信息。

表1 测井数据煤层顶底界面与全波形反演结果对比 m

图11 煤层顶面(蓝色)和底面(红色)分布预测

4 结 语

针对我国煤矿开采的实际需求，研发了煤矿巷道三维弹性波全波形反演算法与技术流程，并应用于陕西小保当一号煤矿112204工作面的实际三维地震数据处理。在反演中，采用了多尺度、多目标函数和多波的反演策略，并使用测井先验信息进行正则化约束，与标准的全波形反演方法相比，可减弱反演的非线性和多解性，获得更加准确的结果。最终构建了112204采煤工作面精细三维弹性速度结构模型，对煤层的顶底板进行了高精度成像，与传统槽波探测方法相比，可对煤层三维空间分布进行精细刻画，满足了采煤工作面的透明化的实际需求。