《1、 引言》
1、 引言
随着土地资源的日益枯竭，岩土工程重心逐渐向地下建设迈进[1‒2]。然而，岩体中隐伏的异常区域给地下岩土工程的开展带来重大挑战。地下岩土工程的安全受到各类潜在地下灾害的威胁[3‒4]，如岩爆[5‒6]、顶板垮塌[7]、突水[8]等。这些潜在危险区域与岩体往往在应力分布上存在差异，在开挖和勘探等地下工程中可能造成极大的财产损失和人员伤亡[9‒10]。由于岩土工程进行中的应力场时刻处于重分布状态[11‒14]，因此仍有待进一步揭示其时空演化规律。利用现有应力场识别技术来探测复杂岩体结构中隐伏的危险区域极具挑战性。
当岩体受到外部应力而发生破坏时，往往伴随着明显的速度降低[15‒16]，因而波速场成像技术能够实现对于处于裂纹稳定和不稳定扩展阶段中内部损伤演化的监测[17]。该技术可进一步应用于采矿工程中潜在危险源的预警，为风险区域的疏散工作提供充足的时间[18]，并给地下无人驾驶车辆提供引导[19]。波速场成像技术对于保障岩土工程活动安全高效开展起着至关重要的作用。
异常区域识别是地下岩土工程中经典和基础的问题之一[20‒21]。20世纪70年代，Aki等[22‒23]开创性地将医学CT技术引入地球物理科学，并在勘探地层结构、油气存储和地质灾害等领域得到不断发展。其中，两点射线追踪方法[24‒26]曾被广泛应用于天然地震层析成像，但由于该方法存在计算效率低、易跟踪盲区、易陷入局部最优等问题，随着射线追踪技术的发展被逐步淘汰。Nakanish和Yamaguchi [27]提出了最短路径法（shortest-path method, SPM），将区域离散化为若干单元体，然后计算单元边界上节点走时，并将节点最小走时连线作为射线路径。Nasr等[28]提出了一种动态SPM（dynamic shortest-path method, DSPM），通过在SPM拓扑中添加和移除节点，进一步提高了精度和效率。Dong等[29‒30]研发了一种改进的A*最短路径算法，结合匹配方法识别二维结构中的空洞区域。除了射线追踪理论外，程函方程也是层析成像的基本原理之一。快速行进法（fast marching method, FMM）结合窄带技术描述波前传播方式，有效地确定网格点的更新顺序[31‒33]，计算区域内的走时。Dong等[34]采用FSM结合主被动震源探讨模拟实验中复杂结构下声层析成像的影响因素。Jiang等[35]提出了一种改进的FMM方法并实现了岩体中损伤区域裂隙的精确定位。提出了一种利用FMM的主-被动联合层析反演算法并将其应用于砂岩压密实验。然而，该方法受岩体与空洞区域之间过大的速度差限制，难以应用于空洞区域的检测。快速扫描法（fast sweeping method, FSM）通过Gauss-Seidel迭代交替扫描实现了走时场程函方程数值解的快速求解[37‒39]。该方法具有高效的正演模拟效率，对层析成像技术的工程应用具有重要意义。层析成像作为一种实用技术，能在大规模地震和油气储藏勘探中定量揭示地质构造[40‒43]。当地下开口中周围积聚的应变能突然释放，其岩石破裂过程往往伴随着微震活动，故而可利用微震信号监测应力的重分布以及微破裂的形成[44‒46]。层析成像对于探测复杂采矿环境下的隐伏异常区域至关重要。
本文提出利用走时层析成像技术，结合矩阵分解和高斯滤波的阻尼最小二乘法（DLSQR-GF），辨识复杂岩体结构中的异常区域，不仅克服了空洞区域辨识的速度差限制，而且减轻了迭代中孤立速度突变的影响。在走时层析成像计算中使用SPM、DSPM和FSM作为正演方法，以获得计算区域的到时。为了量化评估本文方法的辨识精度和计算效率，进行了数值和室内实验，并通过现场应用以验证其可行性。
《2、 方法》
2、 方法
《2.1 正演建模》
2.1 正演建模
精准快速的到时计算能力是层析成像正演建模的基础，这对于异常区域辨识同样也是必不可少的。常用的正演建模方法包括基于射线理论的射线追踪技术和基于程函方程的波前追踪技术，本文主要采用SPM、DSPM和FSM等作为正演建模方法。下面简要介绍这些方法的基本原理。
《2.1.1. 最短路径法》
2.1.1. 最短路径法
SPM源自图论中的最短路径问题，将区域离散为单元，在每个单元的边界上设置节点。通过源节点计算得到相邻节点的到时，然后将其设置为其他节点的次级源。链接最小到时及其对应的次级源节点，并将其设置为射线追踪路径。
《2.1.2. 动态最短路径法》
2.1.2. 动态最短路径法
为了消除网格均匀慢度和平面波假设带来的误差，Nasr和Giroux [28]提出了DSPM来定义节点的慢度。DSPM首先通过改进的SPM计算所有节点到时，并在迭代拓扑中增加和移除次级节点以动态改变慢度精度。在射线路径跟踪过程中DSPM使用最陡走时梯度作为射线路径，然后利用射线路径与网格相交区域的慢度值重新计算每个传感器的到时。DSPM可在不增加网格精度的情况下提高计算精度，从而减少内存和计算成本，在同等SPM计算成本下能够更精准地正演结果。
《2.1.3. 快速扫描法》
2.1.3. 快速扫描法
FSM通过Gauss-Seidel迭代交替扫描计算区域节点的到时。首先，源节点的到时被初始化为零，其他节点的到时设置为极值。然后，采用Gauss-Seidel迭代算法计算程函方程的解，保留原始值和更新值中的最小值。扫描过程持续进行，直到满足收敛条件。FSM正演计算效率高，远远超过射线追踪法，正演计算一直是层析成像中最耗时的部分，FSM的出现对层析成像的工程应用具有重要意义。
《2.2 反演算法》
2.2 反演算法
层析成像的反演问题可视为求解由观测走时和计算走时所组成的最优化问题
（1）
（2）
（3）
式中，
式（3）为欠定和矛盾方程组，可以通过反投影法、牛顿法、DLSQR等迭代反演方法求解。其中，DLSQR是一种基于迭代的反演方法，尤其适用于病态线性方程组的求解，具有计算量小，能够利用矩阵的稀疏性简化计算，对于病态问题具有较好的收敛性和可靠性的特点，适合求解由射线路径所组成的非线性稀疏算子
走时层析由以下步骤组成：
（1）建立初始模型并划分网格。网格的大小和数量决定结果的分辨率，更精细的网格也伴随着更多的计算成本和信号数量。
（2）采集信号数据，包括震源、传感器和断铅点的走时和坐标。采用降噪算法对所采集的信号数据进行去噪处理，以保证到时的准确性。
（3）联合正演模拟和DLSQR-GF进行层析成像计算。通过DLSQR-GF将正演部分计算得到的走时与观测走时进行反演计算，得到迭代过程中的速度场。
（4）当收敛准则达到时停止迭代。扩展信息准则（extended information criterion, EIC）可利用bootstrap统计量确定最佳迭代次数[47]，避免了使用残差作为收敛标准带来的过度迭代。EIC和最大迭代次数共同作为收敛准则。
（5）从走时层析成像结果中识别异常区域。异常区域由与背景速度明显不同的速度网格聚集而成。算法针对异常区域的识别能力可通过统计方法进行评估。
走时层析成像流程如图1所示。
《图1》
图1 层析成像算法流程图，采用SPM、DSPM和FSM作为正演模拟方法，DLSQR-GF作为反演算法。
《3、 实验》
3、 实验
为量化评价走时层析成像对于复杂岩体结构中异常区域的超前辨识能力，本文开展了数值、室内实验和现场试验。数值实验由包含两个椭圆形低速异常区的50 m × 40 m × 20 m掘进模型构成，传感器等距地分布于4个角落。室内实验包含了二维和三维实验，研究区域为一个500 mm × 200 mm × 160 mm的长方体花岗岩试样，其顶部有5个直径（D）为50 mm的钻孔，代表低速异常区。数值和室内实验的参数如表1所示，传感器和异常区域的位置如图2（a）、（b）所示。室内实验中布置传感器按照尽可能增加岩体中的射线覆盖率的原则配置，源自射线路径的慢度场扰动和足够的射线覆盖确保了岩体中异常区域的识别精度。图2（c）、（d）描述了室内二维和三维实验中的网格分布。在室内实验中脉冲和断铅信号作为主动源，脉冲信号的基频在160 kHz左右。通过建立含坐标的解析方程，利用极大似然估计法确定断铅点的起震时间。室内实验采用AMSY-6多通道设备和VS45-H传感器对信号进行高标准采样，传感器的频率范围为20~450 kHz，频率响应在280 kHz处达到峰值。
《表1》
表1 数值实验和室内实验的参数
| Parameter | Numerical experiment | 2D laboratory experiment | 3D laboratory experiment |
|---|---|---|---|
| Model size | |||
| Abnormal region size | |||
| Grid number | |||
| Receiver number | 36 | 12 | 20 |
| Source number | 58 | 12 | 40 |
| Ray number | 2088 | 144 | 800 |
《图2》
图2 室内实验中传感器和网格的分布。（a）二维传感器网络；（b）三维传感器网络；（c）二维网格；（d）三维网格。
超前辨识要求从岩体中快速识别隐伏地异常区域，本文使用识别率（identify rate, IR）和代价率（cost rate, CR）来定量评估走时层析成像的识别效率。当异常区域的速度明显低于背景速度时，可以认为满足有效识别。本文将IR设定为低于背景速度下四分位数的异常区域速度的比例，下四分位值由非空洞区域的网格计算得到，然后与空洞区域的网格进行比较。IR值可由非空洞区域内速度高于下四分位数速度值的网格数与空洞区域内所有网格数的比值得到，IR值越高代表异常区域辨识度越高。
《3.1 数值实验》
3.1 数值实验
为了评估方法在理论条件下的可行性，本文建立如图3（a）所示的数值实验，模拟实际掘进过程中的异常区域检测。目标岩体区域尺寸为50 m × 40 m × 20 m，周围有其他岩体。在目标岩体中埋设了两个不同的异常区域，可能会对掘进过程产生影响。为探测异常区域的分布，在目标岩体四角布置钻孔，沿钻孔间隔5 m布置传感器。传感器布置为目标岩体构建了监测网络。目标岩体的参数化如图3（b）所示。背景速度设为
《图3》
图3 掘进过程中异常区域探测的传感器布设（a）和其速度模型（b）。SPM（c）、DSPM（d）和FSM（e）成像结果的速度模型。
《图4》
图4 初始模型（a）和真实模型（b）。（c）~（h）成像结果：SPM（c）、DSPM（e）、FSM（g）三维数值实验的成像结果； SPM（d）、DSPM（f）和FSM（h）在5%到时误差下的成像结果。
SPM、DSPM和FSM的层析成像结果如图3（c）~（e）所示。从水平速度剖面结果可以看出，低速体位于实际异常区域周围，其中速度较低的异常区位于目标岩体边界附近，而其他异常区位于掌子面附近，与真实数值模型接近。然而，层析反演需要求解一个待定、多解的问题，其欠定程度与慢度场的信号和网格直接相关。因此，反演模型在异常区域的形状和数值方面可能与真实模型存在差异。异常区域外的速度在层析成像结果中也表现出不同程度的变化。低速区周围的单元被赋予比背景速度更高的速度值，以最小化异常区域形状和速度值的差异带来的到时差异。
图4（c）、（e）、（g）展示了初始模型和层析成像结果的透视图。为了更好地展示异常区域分布，将反演模型沿高度进行切片。从仅包含掘进面的初始模型开始，将正演模拟与DLSQR-GF相结合得到层析成像结果。层析成像结果清晰地显示了掘进面和异常区域的分布情况，所提方法实现了预设异常区域的辨识。SPM、DSPM和FSM的IR值分别为93.96%、94.02%和94.16%。尽管使用这三种方法的数值实验中的IR值差别不大，但FSM表现出了相对更好的性能。以SPM的计算代价为评价标准，DSPM和FSM的CR值分别为28.92%和18.13%。FSM的计算成本明显低于其他两种方法，具有更高的计算效率。为了验证所提方法对噪声的鲁棒性，在理论走时中加入了平均5%的到时误差，结果如图4（d）、（f）和（h）所示。结果表明，虽然在该情况下成像精度受到一定程度的影响，但仍实现了预设异常区域的辨识。在三维层析成像数值实验中，DSPM和FSM在异常区域辨识中表现出优异的识别能力和计算效率。SPM虽然满足识别异常区域的要求，但其计算成本明显高于其他方法。
《3.2 室内实验》
3.2 室内实验
与数值实验中的理想条件不同，在层析成像室内实验中，误差因素是不可避免的，如走时误差、非均匀背景速度、传感器-介质耦合等，层析成像的精度也受到这些因素的影响。因此，有必要对走时层析成像在室内实验中的应用进行验证与量化评估研究。
《3.2.1. 二维室内实验》
3.2.1. 二维室内实验
传感器分布如图2（a）所示。两侧传感器分别用作震源端和接收端。设背景速度
《图5》
图5 初始模型（a）、SPM（b）、DSPM（c）、FSM（d）的成像和射线追踪结果；初始模型（e）、SPM（f）、DSPM（g）、FSM（h）的异常区域分布对比。
结果表明，两个左侧异常空洞区域之间存在一个低速带，在其他异常空洞区域周围也存在扩展的低速带。这可能是由钻孔过程中操作不当造成的，虽然仅从岩样表面没有观测到明显的裂痕，但可能存在内部损伤影响了其速度结构。边界附近的低速伪影可能是由射线覆盖不足和正演方法的稳定性不佳造成的。波速场的反演是一个欠定问题，当反演过程中存在信息不足或误差时，就可能会诱发伪影的产生。
层析成像室内实验受误差因素影响，包括到时误差、背景速度、传感器-介质耦合等，增加了异常区域超前识别的应用难度。其中，SPM受到严重影响，与其他方法相比，反演结果不理想；DSPM和FSM的辨识能力相当，但DSPM在计算效率方面表现出更好的性能。
《3.2.2. 三维室内实验》
3.2.2. 三维室内实验
传感器分布如图2（b）所示。传感器作为接收端，断铅点作为震源端。SPM、DSPM和FSM层析成像结果的水平切片和垂直切片如图6所示。低速区速度大小与背景速度大小明显不同，且与空洞的位置相对应。然而，三维层析成像实验对信号数量的要求较二维层析成像实验更严格，且更容易受到误差因素的影响，这加剧了三维异常区域辨识的难度，并易在射线路径无法覆盖的区域诱发伪影的产生。波速场顶部水平切片上出现较大的低速区，SPM成像结果中无法准确识别空洞的实际分布。DSPM和FSM的层析成像结果的伪影较少，但仅能更新顶部的浅层水平切片。
《图6》
图6 三维层析成像室内实验中的成像结果。（a）初始模型；（b）SPM；（c）DSPM；（D）FSM。（ⅰ）为水平切片；（ⅱ）为垂直切片。
由于三维网格与二维网格相比其网格数量剧增，因此在计算过程中也需要更多有效信号。然而，断铅产生的射线不能完全覆盖空洞，因此在有限信号下很难有效识别空洞。而且三维结构中射线路径变化更复杂，进一步增加了层析成像过程的难度。由于异常区域覆盖不足，低速带模糊了各个异常体之间的区域。SPM、DSPM和FSM等方法的IR均在0.32左右。与二维层析成像室内实验相比，三维层析成像实验虽可以为分析速度场提供更全面的视角，但射线追踪路径更复杂，对用于三维层析成像的信号提出了更为严格的要求。
《3.3 现场试验》
3.3 现场试验
为验证所提方法在实际工程中的可行性，在陕西震奥矿山开展了应用。如图7所示，安装的微震监测系统分布在795~860层，由26个单分量传感器组成。用于层析成像的信号数据为2022年4月1—4日的事件。目标区域可以分为两部分；第一部分包括795-1分层和795-2-810分层（A分层），第二部分包括795-4-842分层和860-1分层（B分层）。根据定位结果和传感器分布情况，建立的模型尺寸约为300 m × 520 m × 200 m。在对比数值实验和室内实验层析成像性能的基础上，将DSPM和FSM应用于现场实验。反演部分采用双差法消除起震时间误差的影响。
《图7》
图7 陕西震奥矿山井下巷道分布及微震监测系统。
A、B分层的初始模型如图8（a）、（b）所示，对应岩体速度为4500 m·s-1。DSPM和FSM的成像结果如图8（c）~（f）所示，其中低速区对应井下巷道分布。结果表明，所提方法能够将井下巷道与岩体区分开来，等同于在掘进过程中进行异常区域检测。成像结果存在高速异常区域，主要分布在传感器和爆破点附近。一方面，这些高速区可能是由大量射线引起的速度扰动造成的，在传感器和爆破点附近的速度有明显变化趋势；另一方面，掘进面高速区可能是由应力条件改变引起的变形和速度变化。DSPM和FSM辨识结果与现场实际情况一致。
《图8》
图8 初始模型在A分层（a）和B分层（b）的水平切片；DSPM在A分层（c）和B分层（d）的水平切片成像结果；FSM在A分层（e）和B分层（f）的水平切片成像结果。
尽管现场试验受到复杂岩体结构和误差因素的影响，所提出的方法仍然实现了从岩体内部对地下巷道的重构。若提供更多样的数据和更符合实际的初始模型，则可以进一步改善反演结果。
《4、 结论》
4、 结论
本文提出了一种用于识别复杂岩体结构中隐伏异常区域的走时层析成像方法。该方法通过SPM、DSPM和FSM等正演模拟方法获得计算区域的走时，利用DLSQR-GF迭代优化观测走时与计算走时之间的残差，克服了空区辨识中速度差限制，缓解了优化过程中孤立速度突变的影响。为了定量评估正演模拟方法与DLSQR-GF结合的识别精度和计算效率，开展了数值和室内实验，并通过对比目标范围内异常区域辨识的IR和CR结果证明了DSPM和FSM的可靠性和有效性。陕西震奥矿山现场试验结果验证了所提方法工程应用的可行性，该方法利用了矿山开采中爆破、微震等多类声源分布广的特点，实现了对井下巷道和高应力区的重建，满足异常区域超前辨识精度和效率要求。
传感器和震源的分布能够影响实验中异常区域辨识的准确性，使用被动源对于三维异常区域的识别更加可靠和有效。此外，通过构建更合理的传感器布局和初始模型，可以进一步提高辨识结果。该方法将为隧道工程、采矿工程等地下岩土工程中潜在危险区域的监测提供理论基础和技术支持。
京公网安备 11010502051620号
