技术领域
[0001] 本发明属于地质工程技术领域，涉及一种采空区探测方法。
相关背景技术
[0002] 煤矿开采造成的采空区是指煤层采掘后形成的空洞区域，主要包括煤层采空、巷道工程及硐室三个部分。采空区的形成演化直接关系到矿区地表沉陷、岩溶塌陷、水资源破坏、生态环境恶化问题，是制约矿区可持续发展的瓶颈。因此，及时探测查明采空区的空间分布和稳定性状况，对保障煤炭开采安全、指导矿区生态修复、促进矿区经济转型具有重要意义。
[0003] 常规的采空区探测方法主要有钻探法、物探法和矿测法，但上述方法存在成本高、分辨率低、实时性差的问题。近年来，物探技术在采空区探测中得到广泛应用和快速发展，如采用瞬变电磁法对某采煤沉陷区开展了探测研究，揭示了沉陷区的三维电性结构；利用浅层地震反射法确定了采空区顶板破坏带的埋深和范围；综合重力和视电阻率数据反演了采空区的密度和电阻率分布。但目前针对采空区探测的联合反演方法多局限于定性描述，缺乏统一的数学模型和反演框架，难以准确刻画复杂条件下采空区的精细结构。
具体实施方式
[0078] 为了使本申请所要解决的技术问题、技术方案及有益效果更加清楚明白，以下结合附图及实施例，对本申请进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本申请，并不用于限定本申请。
[0079] 基于瞬变电磁与地震波场联合反演的采空区探测方法，流程图如图1所示，具体如下所述。
[0080] 首先，根据采空区的埋藏条件，在采空区所处的隧洞穿越区域布设瞬变电磁探测系统和反射地震探测系统。
[0081] 具体地，瞬变电磁探测系统的探测频率为0.1‑100Hz，发射功率不小于500W；反射地震探测系统的激发频率为20‑200Hz，采用可控震源激发。根据采空区埋藏条件，优化两种探测系统的空间布设方式和测线参数。
[0082] 然后，采集瞬变电磁视电阻率和反射地震波形数据，并进行包括数据配准、噪音压制、同步成像在内的预处理，提取采空区电性及弹性差异特征参数。
[0083] 具体地，采空区电性及弹性差异特征参数由瞬变电磁多属性特征参数和反射地震多属性特征参数组成，瞬变电磁多属性特征参数有：视电阻率、视电导率、感应极化系数、瞬变电磁波形衰减系数，反射地震多属性特征参数有：地震波初至走时、振幅、频率、相位、反射系数、衰减系数。
[0084] 然后，基于采空区电性及弹性差异特征参数，提取得到多属性特征数据集，采用机器学习分类算法对多属性特征数据集进行训练、验证及测试，建立采空区分类模型，对采空区类型、规模、稳定性进行判别。
[0085] 具体地，采空区分类模型的建立方法包括：
[0086] S1、基于采空区电性及弹性差异特征参数，提取得到多属性特征数据集；
[0087] S2、采用Pearson相关系数、互信息指标对多属性特征数据集中各特征与采空区类型的相关性，用递归特征消除和序列特征选择方法筛选最佳特征组合，作为对采空区分类贡献度最大的特征子集；
[0088] S3、采用随机分层抽样的方式将特征子集划分为训练集、验证集和测试集三部分；
[0089] S4、机器学习分类算法采用支持向量机，在高维特征空间中构建最优分类超平面，使不同类别的样本间隔最大化，将训练集输入至支持向量机中进行训练，支持向量机的数学模型为:
[0090] min12||w||2+C∑i＝1Nξi，
[0091] s.t.yi(w·xi+b)≥1‑ξi，ξi≥0，i＝1，2，…，N；
[0092] 其中，为分类超平面法向量，b为偏置项，C为惩罚因子，ξi为松弛变量；
[0093] 通过求解上述优化问题，得到分类决策函数:
[0094] f(x)＝sign(∑i＝1Nαiyi(x·xi)+b)，
[0095] 式中，αi为拉格朗日乘子，yi为支持向量，(x·xi)为核函数；
[0096] 训练过程中，采用网格搜索、交叉验证方法优选包括核函数类型、惩罚因子、核函数参数在内的模型超参数，得到训练后模型；
[0097] S5、利用测试集评估训练后模型的采空区分类模型性能，采用精度、召回率、F1值、ROC曲线指标综合评价分类效果；同时学习曲线分析训练集和验证集精度随样本量的变化趋势，判断训练后模型是否存在过拟合或欠拟合；对于表现欠佳的训练后模型，返回步骤S2‑S4，调整特征选择、优化模型参数，直至获得满意的分类性能的模型，将其作为采空区分类模型。
[0098] 然后，基于麦克斯韦方程组和弹性波动方程，采用有限差分和有限元方法进行离散求解，建立采空区物性参数与采空区电性及弹性差异特征参数的定量响应关系，构建采空区瞬变电磁与地震波场联合正演数值模型。
[0099] 具体地，瞬变电磁的理论基础是麦克斯韦方程组，表达式为:
[0100]
[0101] 式中，E为电场强度，H为磁场强度，Js为电流源密度，D为电位移，B为磁感应强度，σ为电导率，ρ为电荷密度；
[0102] 考虑准静态近似条件，忽略位移电流，上述方程组可简化为:
[0103]
[0104] 在采空区瞬变电磁数值模拟中，采用交错网格有限差分方法对上述方程组进行离散求解，首先将研究区域剖分为Nx×Ny×Nz的立方体网格，然后将电场、磁场分量分别配置在网格棱边和面的中心点，电导率、磁导率等参数赋值在网格中心；时间方向采用二阶精度的中心差分格式，空间方向采用二阶精度的交错网格差分算子，建立三维瞬变电磁正演数值方程为：
[0105] En+1‑En‑1，
[0106]
[0107] Hn+1‑Hn，
[0108]
[0109] 将上述方程在整个空间网格和时间步长上进行组装，形成大型稀疏线性方程组，采用直接法或迭代法进行求解，即可得到三维电磁场时空分布。
[0110] 具体地，反射地震数值模拟采用弹性波动方程，表达式为:
[0111]
[0112] 式中，u为位移矢量，σ为应力张量，ρ为密度，f为震源项；
[0113] 应力‑应变关系满足胡克定律:
[0114] σ＝λtr(ε)δ+2με，
[0115]
[0116] 式中，δ为单位张量，λ和μ分别为第一拉梅常数和第二拉梅常数；
[0117] 在数值求解中，采用交错网格有限差分方法，应力分量配置在网格中心，速度分量配置在相应棱边中点；时间和空间方向均采用二阶精度的中心差分格式，建立速度‑应力形式的弹性波动方程差分格式：
[0118]
[0119] 其中，v为速度分量；
[0120] 对全部网格和时间步长组装得到巨型稀疏矩阵方程组，采用显式时间步进方法进行求解，得到三维弹性波场时空演化信息。
[0121] 更具体地，采空区瞬变电磁与地震波场联合正演数值模型的构建方法包括：
[0122] S1、根据先验地质信息，利用建模软件构建采空区初始模型，包括采空区的空间展布形态以及围岩的电阻率、密度、波速物性参数；
[0123] S2、设置包括发射电流、波形在内的瞬变电磁发射源激发参数和包括接收通道数、采样间隔在内的接收系统参数，在采空区初始模型中配置发射源和接收点的空间位置；
[0124] S3、基于三维电磁场时空分布，依据发射源和接收点的空间位置，正演计算采空区初始模型的三维瞬变电磁响应，得到不同时间和空间位置的电场、磁场分量，并合成视电阻率表征参数；
[0125] S4、设置包括震源类型、主频在内的反射地震激发源参数和包括检波器组数、道间距在内的接收系统参数，在初始模型中配置炮点和检波点的空间位置；
[0126] S5、基于三维弹性波场时空演化信息，依据炮点和检波点的空间位置，正演计算采空区初始模型的三维地震波场快照，提取不同道集的地震记录，并拾取包括初至走时、频散、振幅衰减在内的属性参数；
[0127] S6、通过改变采空区初始模型的包括空区范围、厚度在内的几何参数和包括空区充填物电阻率、密度在内的物性参数，重复步骤S2‑S6，开展电磁和地震响应的参数敏感性分析，定量刻画采空区物性参数与采空区电性及弹性差异特征参数的定量响应关系；
[0128] S7、依据不同时间和空间位置的电场、磁场分量，视电阻率表征参数，不同道集的地震记录，包括初至走时、频散、振幅衰减在内的属性参数，以及采空区物性参数与采空区电性及弹性差异特征参数的定量响应关系，构建采空区瞬变电磁与地震波场联合正演数值模型。
[0129] 然后，依据采空区瞬变电磁与地震波场联合正演数值模型，结合采空区包括几何形态、围岩应力状态地质在内的先验信息，构建电磁‑地震联合反演目标函数。
[0130] 具体地，电磁‑地震联合反演目标函数为：
[0131] Φ(m)＝1||Wd(d‑f(m))||2+α||Wm(m‑m0)||2+β||L(m)||2，
[0132] 式中，第一项为数据拟合项，d为观测数据，包含视电阻率、地震走时、波幅信息；f(m)为采空区模型m的理论响应，Wd为数据加权矩阵，用于平衡不同数据类型的贡献；第二项为模型约束项，m0为参考模型，通常基于地质资料和经验设定，Wm为模型加权矩阵，对反演模型的空间分布进行约束，α为模型约束因子，控制先验信息对反演结果的影响程度；第三项为模型光滑项，L为空间二阶差分算子，β为光滑约束因子，控制反演模型的光滑性，其中模型参数m主要包括采空区各网格的电阻率、密度、纵横波速度。
[0133] 然后，依据电磁‑地震联合反演目标函数，采用贝叶斯‑马尔可夫链蒙特卡洛方法进行反演求解，得到采空区电阻率、地震波速三维分布及其不确定性定量评估结果。
[0134] 具体地，电磁‑地震联合反演目标函数的求解过程为：
[0135] 将电磁‑地震联合反演目标函数转化为求解如下极小化问题：
[0136] minΦ(m)，s.t.ml≤m≤mu，
[0137] 式中，ml和mu分别为模型参数的上下限，反映已知的地质约束条件；
[0138] 采用贝叶斯‑马尔可夫链蒙特卡洛方法进行反演求解，给定观测数据d，采空区模型m的后验概率p(m|d)正比于似然函数p(d|m)和先验概率p(m)的乘积，即：p(m|d)∝p(d|m)p(m)；其中，似然函数p(d|m)体现了模型m与观测数据d的拟合程度，先验概率p(m)体现了模型参数的物理约束和统计特征；
[0139] 通过构造马尔可夫链实现后验概率分布的随机抽样，进而得到采空区模型参数的最优估计和置信区间，流程如下：
[0140] S1、基于地质资料和经验判断，构建采空区初始模型，作为马尔可夫链的起点；
[0141] S2、基于Metropolis‑Hastings准则，按照概率p从当前模型mi扰动到下一状态mi+1；其中转移概率p的表达式为:
[0142] p＝min{1，p(mi+1|d)q(mi|mi+1)p(mi|d)q(mi+1|mi)}，
[0143] 式中，q(mi+1|mi)为马尔可夫链的转移核函数，通常选择对称分布如高斯分布等；
[0144] S3、重复步骤S2，不断更新和调整马尔可夫链的状态，直至达到收敛条件：后验概率无明显变化；
[0145] S4、基于马尔可夫链的平稳分布，随机抽取一系列反演模型样本，对其进行统计分析，得到采空区模型参数的最优估计和置信区间，得到采空区电阻率、地震波速三维分布及其不确定性定量评估结果。
[0146] 最后，依据采空区类型、规模、稳定性判别结果，结合采空区电阻率、地震波速三维分布及其不确定性定量评估结果，建立采空区三维精细结构与物性模型。
[0147] 实施例
[0148] 针对某水利枢纽工程输水隧洞穿越老采空区段，采用本发明对该采空区进行探测。
[0149] 根据已知采空区分布，设计瞬变电磁测线8条，测线长度5km，测点间距20m；布设地震测线4条，测线长度5km，检波点间距5m。瞬变电磁发射源采用美国EMIT公司的EMP‑400T发2
射机，最大发射功率800A·m ，频率0.3Hz～300Hz；地震激发源选用美国IDE公司的BIS‑SHD50电火花震源，单点激发能量2000J。瞬变电磁接收设备采用加拿大Phoenix公司V8多通道接收仪，共150个接收通道；地震采集采用Inova公司ARIES‑Ⅱ系统，共1000个检波点。
[0150] 对全测线的瞬变电磁视电阻率与地震三分量波形资料采集完成后，对电磁数据进行仪器响应校正、同步叠加、三维空间插值预处理，信噪比提高5倍以上；对地震数据进行静校正、振幅恢复、面波衰减预处理，有效压制高频随机噪音，优化同相轴响应，获得高质量叠后偏移剖面。提取了地震反射波频散、振幅、视速度等12个属性，协同表征采空区电性、弹性差异特征。
[0151] 选取典型采空区地段，构建了覆盖层、基岩、采空区三层介质地震‑电磁波场耦合模型。根据地质资料，设置空区充填类型：裂隙型、崩落型、充水型，厚度为1m‑20m，埋深为10m‑100m，开展全响应数值模拟，获得不同类型采空区的视电阻率5Ω·m‑500Ω·m和地震波速300m/s‑2500m/s的响应特征，建立了"采空区发育程度‑物性参数"映射关系，如图2所示，为联合反演提供约束。
[0152] 以正演模拟数据为训练样本，融合实测视电阻率、地震波频散等8维多源属性，构建了采空区分类模型。在井孔资料约束下，对采空区类型、规模、稳定性分类标签进行智能标定，对现场采集数据进行分类识别，采空区判别准确率达85％，较传统经验判读方法提高20％以上。
[0153] 以瞬变电磁视电阻率和地震波反射波速度为约束，同时加入采空区埋藏规律、围岩应力状态等先验地质信息，构建了电磁‑地震联合反演目标函数:
[0154] min{||Wd(dobs‑F(m))||2+α||Wm(m‑m0)||2+β||Lm||2},
[0155] 其中,m为采空区空间分布、物性参数；dobs为实测物探数据；F为正演算子；Wd、Wm分别为数据、模型加权矩阵；m0为先验模型；L为光滑约束算子；α、β为先验信息与光滑约束权重。利用马尔可夫链蒙特卡洛方法，对采空区模型进行随机扰动，获得反演结果的后验分布，生成置信概率图，定量评估了反演参数的不确定性。
[0156] 综合瞬变电磁、地震联合反演及不确定性分析结果，结合钻孔揭露，建立了XX输水隧洞沿线采空区三维精细结构模型,如图3所示，圈定了7个高风险采空区段，整体定位精度优于5m，规模估算误差小于10％。采空区最大埋深50余米，体积5万立方米以上。根据采空区稳定性预测，提出了隧洞优化线路方案,如图4所示，局部绕避2处强耦合采空区，采空区失稳风险系数降低80％以上，节约工程投资3000万元。指导现场开展了针对性采空区注浆加固处理，实现了XX引水隧洞的安全高效穿越施工，创造了显著的工程和社会效益。
[0157] 本领域内的技术人员应明白，本申请的实施例可提供为方法、系统、或计算机程序产品。因此，本申请可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本申请可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD‑ROM、光学存储器等)上实施的计算机程序产品的形式。本申请实施例中的方案可以采用各种计算机语言实现，例如，面向对象的程序设计语言Java、C++、Python和直译式脚本语言JavaScript等。
[0158] 本申请是参照根据本申请实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
[0159] 这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
[0160] 这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
[0161] 尽管已描述了本申请的优选实施例，但本领域内的技术人员一旦得知了基本创造性概念，则可对这些实施例作出另外的变更和修改。所以，所附权利要求意欲解释为包括优选实施例以及落入本申请范围的所有变更和修改。
[0162] 显然，本领域的技术人员可以对本申请进行各种改动和变型而不脱离本申请的精神和范围。这样，倘若本申请的这些修改和变型属于本申请权利要求及其等同技术的范围之内，则本申请也意图包含这些改动和变型在内。
