《1. 引言》

1. 引言

随着浅部矿产资源的逐渐枯竭和地下空间的利用，地下工程逐渐向深部进发。在深部区域，巨大的地应力加上构造应力，使得岩体积蓄了大量的能量。地下开挖如同卸载，使得岩体积蓄的能量瞬间猛烈释放，形成岩爆等灾害。这种高应力引起的灾害往往对岩体工程造成很大破坏，且常造成人员的伤亡。岩爆灾害涉及诸多复杂因素的岩体力学行为过程，很难用传统的岩石力学理论来进行分析和解释，是目前的世界性难题之一。

微震监测技术作为一种有效手段，可以对岩爆进行监测预警[1]。目前，微震监测技术已经广泛应用于深、高应力矿山的地压安全监测[2–6]。此外，加拿大、南非等国家也建立了全国性的矿山微震监测网络系统 [7,8]。并且，微震监测技术还在航空、军工、桥梁结构、隧道工程、精密制造等各个领域得到了迅速发展和应用 [9–13]。

震源定位是微震/声发射监测中最经典和最基本的问题之一[13]。1912年，Geiger [14]首次提出了定位方法。此后，提高微震震源定位的精度和效率，一直是学者们在微震监测技术领域的重要研究内容。随着大量定位方法的提出，震源定位的误差和所需时间不断降低。经典的源定位方法有Inglada方法[15]、USBM方法[16,17]、 Thueber方法[18]、Powell方法[19]、单纯型法[20]、双重残差法[21]等。近些年，Dong和Li等[22,23]研究了波速对定位的影响，提出了无需预先测速的定位方法。该方法避免了测量速度区域与实际微震区域的差异，大幅提高了定位精度。该方法还节约了前期测量速度造成的人员、时间和经济成本，比传统方法更加方便实用。

Dong等在此基础上发展出了三维综合解析定位 [24,25]、多步定位[26]、数值解析协同定位[27]等方法，取得了良好的效果。随着研究的深入，一些学者开始关注物体的几何形状对定位的影响。Park等[9]、Baxter等 [28]、Eaton等[29]通过训练数据建立数据集合来评估震源位置。这类方法原理上可以适用于任何发生在物体表面震源的定位，但对于块体结构内部震源的定位，其适应性较弱。目前关于此类研究的试验对象均为平面结构，且定位方法本身非常耗时。Gollob等[30]考虑了对象的几何情况对定位的影响，但没有考虑预测波速与实际波速之间的误差。

本文针对含孔洞的三维结构，提出了一种基于A* 搜索算法的无需测速震源定位方法（VFH）。该方法通过等距网格点来搜索路径的方式避免了人工重复训练，引入A* 算法[31,32]并采用网格点来适应含有空区的复杂结构，并利用了无需预先测速的定位方法的优点。该方法能有效地适应于三维复杂结构，并显著提高了震源的定位精度。

《2. 理论与方法》

2. 理论与方法

VFH使用改进的A* 搜索算法和无需预先测速定位方法的思想来确定声发射源的位置。该方法可分为四个部分：首先，划分需要定位物体的网格并用0和1表示物体的形状；其次，采集传感器接收到的声发射事件产生的P波到时信号；然后，使用A* 算法搜索每个传感器与每个网格点之间的最快波形路径；最后，根据引入最小偏差量D 来确定声发射源的位置。图1为VFH的定位流程图。本节将描述该方法的每一部分。

《图1》

图1. VFH的定位过程流程图。

《2.1. 确定初始环境》

2.1. 确定初始环境

在定位区域上，确定空区的几何形状和具体位置。根据空区的情况和定位的精度要求，确定单位立方体网格的尺寸大小。一般来说，网格划得越密，定位精度越高，计算量也会成倍地增加。但当网格足够密时，继续增加网格密度，定位精度将不会有明显变化。建立和网格节点相同尺寸的零矩阵M，将矩阵索引位置(x, y, z )与网格节点位置一一对应，并将对应空区的M (x, y, z )值更改为1。网格节点形成一个集合，在后续节点间搜索最快波形路径时它们被作为起始点。假定P 波在周围非空区域的传播速度为一个定值，为未知数，用C 来表示。对于未知震源P0，设其位置坐标为(x0, y0, z0)，激发的初始时间为t0。

《2.2. 搜索最快波形路径》

2.2. 搜索最快波形路径

将集合内的每个网格节点Pxyz 当作潜在震源的激发位置，追踪最短路径，得到网格节点到第k 个台站的理论最短路径。若网格节点Pxyz位于空区内，则认为 。VFH采用改进的A* 算法来追踪最短路径。

2.2.1. A* 搜索算法

A* 搜索算法由斯坦福研究所（现为SRI国际）的 Peter Hart、Nils Nilsson和Bertram Raphael于1968年首次提出[31,32]。在计算机科学中，A* 算法作为Dijkstra 算法的扩展，因其高效性而被广泛应用于搜寻最短路径及图的遍历。A* 算法通过使用启发式的方法来指导搜索路径，从而实现更好的性能。在其主循环的每次迭代中， A* 算法根据当前节点到终点的估价代价来确定扩展路径。具体来说，A* 算法始终选择最小长度的路径。

式中，m 表示选择路径的当前节点；表示从起始节点到节点m 的实际最快路径长度；而 表示一个启发式函数，用于估算从节点m 到目标的最快路径长度[33]。在本文中估值函数采用欧氏距离。当从起始点拓展到目标点或者没有路径可以拓展时，A* 算法停止搜索。

在实际三维应用中，传统的A* 算法采用中心点，一般只考虑相邻层的26个节点来选择下一个节点，如图 2 （a）、（b）所示。在L形的区域内，利用传统的A* 算法在内追踪最短路径，搜索到如图2（c）所示的一条路径。从图中可以发现，追踪的最短路径存在两处不合理的地方：①搜索得到的路径有明显的锯齿状，这是传统的A* 算法的自身限制造成的；②路径的节点均为立方体网格的中心，这意味着传感器也要安放在立方体网格的中心，与实际不符合。

《图2》

图2. 传统的A* 搜索算法。（a）当前节点连接到相关的26个节点；（b）图2（a）的主要视图；（c）传统A* 搜索算法搜索的路径。

2.2.2. 改进的 A* 算法

为了更有效地追踪最短路径，我们改进了A* 算法，采用网格点进行搜索，见图3（a）。这样可以使得搜索得到的路径的节点均在网格节点上。这也意味着传感器可以贴在物体的表面节点上，从而更加符合实际情况。

为了让搜索得到的路径不具有明显的锯齿状，我们让节点与周围更多层的节点建立有效联系。传统的A* 算法中，一个节点向相邻一层的26个节点进行拓展。这意味着当前节点向周围拓展的方向只有26个可以选择。

让节点与周围更多相邻层的节点之间进行联系，可以使得当前节点搜索路径时选择的方向更多，这使得搜索得到的路径更精确。根据节点拓展模型的对称性，只画出模型的八分之一来进行解释说明。图3（b）、（c）分别显示了当前网格节点与周围两层（124个节点）、三层（342个节点）建立了联系。节点与层数i 之间的关系表示为：

《图3》

图3. 当前网格节点与关联网格节点间形成的方向。（a）当前节点连接到相邻的一层；（b）当前节点连接到相邻的两层；（c）当前节点连接到相邻的三层。

在向外拓展的过程中，部分节点间形成的方向重复，故可以不用考虑（图1）。去除这些方向对应的节点间连接可以减少计算量。随着每个节点与周围更多的网格节点建立联系，搜索得到的路径的误差越小，但同时带来计算量的增加。

设计一个块体模型和一个长条状模型来探讨层数 i的合理取值。将模型划分网格，如图4所示。假设在 Source点触发震源，波从Source点到达K 点（K = A,B,C, D,E,F,G,H,A′,B′,C′,D′,E′,F′,G′,H′），形成路径LK。将波的实际最快路径距离DR 记录到表1中。分别使用i 层（i = 1,2,3,4,5,6,…）的模型追踪路径LK 。搜索得到的路径距离DSi 与DR 之间的相对误差E 可以表示为：

将得到的路径的相对误差E 记录到表1中，并选出每个模型的最大路径误差Ei –max。算法的计算量O 与拓展的节点数Z 呈正相关关系。因而，可以用节点数Z 来近似地表示计算量O，即。图（5）显示层数i 分别与Ei –max和关联节点数之间的关系。

《图4》

图4. 节点间的理想路径。

《表1》

表1 真实路径距离与追踪路径距离之间的相对误差

从图5中可以看出，当拓展的层数达到3层时，理论路径的相对误差小于3%，能基本满足定位要求，且计算量增速相对较小。因此，本文中的A* 算法考虑当前节点与周围3层的节点建立联系。为了确定路径的具体位置，A* （或Dijkstra）算法在确定最短路径的长度后，需再次进行反向搜索。本文在搜索路径时增加一个数组，在当前节点对应的位置上记录前一个节点的坐标。这样就避免了反向搜索，提高了运算效率。

《图5》

图5. 相对误差和相关节点数与模型层数的关系。

《2.3. 采集到时数据》

2.3. 采集到时数据

在待测物体的不同位置分别装上m 个传感器。各个传感器均位于网格点上，其位置均为已知。对于三维模型，未知数有5个[P 波的波速、声发射源坐标(x0, y0, z0)、激发的初始时间t0]，因而m 需为大于或等于5的整数。对于接收到信号的第k 个传感器Sk，记录其位置坐标为(xk, yk, zk )，接收到声发射P波信号的初至到时为 。计算两个不同的传感器和的实际到时差，用表示。

《2.4. 震源定位》

2.4. 震源定位

对于点Pxyz 激发的震源，理论旅行时间等于震源到第k个传感器之间的最短传播路径除以波速C。两个不同的台站和m 的到时差等于旅行时间之差。

根据与之差的平方，引入Dxyz 来描述点 Pxyz 与未知声发射源P0的偏离程度，表示为：

式中，当采样点落入空区内，则令Dxyz = 。

每个网格点都将得到对应的Dxyz 值。Dxyz 的值越大，表示点Pxyz 与未知震源P0 的偏离程度越大。因此，最小的Dxyz  值对应的坐标便可认为是微震/AE定位的坐标。

《3. 结果与讨论》

3. 结果与讨论

为了评价新的微震震源定位算法性能，对空心柱体砂浆结构件进行了断铅试验。将10 cm×10 cm×10 cm 的立方体试件中间挖去了尺寸为Φ 6 cm×10 cm的圆柱体，如图6（a）所示。为了更快地计算，将该试件划分成25×25×25个相同尺寸的立方体小网格方块，形成定位模型，如图6（b）所示。

将空区所在方块标记为不通过，其他方块标记为通过。将6个声发射传感器固定在试件上，在传感器与试件之间加入耦合剂来获得更好的声耦合。传感器的位置均位于网格节点上[图6（c）]，坐标位置见表2。

《图6》

图6. VFH定位过程。（a）试件；（b）试件建模和网格划分；（c）确定传感器的坐标（单位：cm）；（d）使用改进的A* 搜索算法进行的从节点到传感器的路径搜索（单位：cm）。

《表2》

表2 传感器布置的坐标

为方便计算，将模型坐标进行转化，来使得节点矩阵M 的索引与坐标位置一一对应。AE数据采集时采用 40 dB的阈值和5 MHz的采样率。在试件不同位置进行断铅试验，每个位置进行两次，记录传感器接收声发射事件产生的P波到时，记录在表3中。搜索每个潜在震源到传感器的路径LS，并计算距离，如图6（d）所示。再根据传感器接收到的到时，计算每个潜在震源的偏差值Dxyz 。确定最小Dxyz 值对应在试件上的位置坐标(x, y, z )，并将位置坐标进行转化，结果见表3。同时，用传统的Geiger法进行定位，定位结果记录在表3中。

《表3》

表3 断铅点的定位结果和误差

分析定位结果，将新方法和传统方法的定位结果与实际断铅点进行对比，将误差记录在表3内。从表中可以看出，传统的Geiger法的最大定位误差为4.5 cm，远大于VFH的定位误差。为了便于观察，图7（a）为两种方法的定位结果和定位误差的可视化图形。圆圈尺寸代表震源定位结果的误差大小。可以明显地看出传统定位方法的圆圈相较于VFH的要大很多。图7（b）为两种方法定位误差的箱线图。在图7（b）中，使用VFH得到的震源定位误差的中位数约为1.0 cm，而使用传统方法得到的震源定位误差的中位数约为1.9 cm。根据表中每个断铅试验的定位误差，可以容易地计算出新方法的平均值为1.20 cm，而传统方法的定位误差为2.02 cm。新方法的平均定位精度较传统方法提高了近40%。因此，在复杂的三维结构中，新方法的定位精度较传统方法有了较大的提高。

《图7》

图7. 两种方法得到的震源定位结果及误差。（a）定位结果及误差可视化；（b）定位误差箱线图。

《4. 结论》

4. 结论

在本文中，我们提出了一种基于A* 算法的无需测速定位方法（VFH）。我们通过等距网格点来搜索路径的方式避免了人工重复训练。对于含有不规则空区的三维复杂结构，我们改进了A* 算法并采用网格点来搜索符合弹性波实际传播规律的路径。并且，我们还在计算中将弹性波的波速当作未知数，来减小因监测过程中物体波速的变化带来的定位误差。

为了评价新方法的准确性和有效性，在空心立方体混凝土结构件上进行了断铅试验。定位结果显示，新方法得到的平均定位误差为1.20 cm，它的平均定位精度较传统的Geiger法提高了40%。这表明新方法能有效地适应几何不规则的三维复杂结构。

该方法克服了传统定位方法中波速难以确定的缺陷，考虑了几何形状对定位的影响，且适用于三维定位，完善了微震定位方法。本文的方法亦适用于涉及声发射定位的无损检测等其他领域。

《致谢》

致谢

感谢国家自然科学基金（51822407、51774327）、湖南省自然科学基金（2018JJ1037）、中南大学创新驱动计划项目（2020CX014）的资助。

《Compliance with ethics guidelines》

Compliance with ethics guidelines

Longjun Dong, Qingchun Hu, Xiaojie Tong, and Youfang Liu declare that they have no conflict of interest or financial conflicts to disclose.

《Appendix A. Supplementary data》

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.eng.2019.12.016.