如何使震源定位更精确、快速，是学者研究微震监测的基础问题之一。目前，研究人员主要通过拾取P波初至时间来进行反演定位［9］。王辉等［10］基于单纯形-最短路径射线追踪法建立了微震震源混合定位算法，为复杂地区微震监测提供了理论和数据支撑；贾宝新等［11］基于高密度台阵和粒子群算法提出了新的微震源定位方法，为中小尺度微震的研究提供了参考；Li等［12］通过研究传感器阵列的空间布局与放置方向，提出了传感器网络优化原则，使定位精度得到一定程度的提高。这些定位方法均要预先设定一个波速模型，而实际矿山中的岩体波速难以真实测量，因此基于预先建立波速模型下的定位方法具有一定局限性，导致定位误差较大。
董陇军等［13］提出了3种无需预先建立波速模型的震源定位方法，节约了预先测速的成本，避免了测速区与实际微震区的波速差异导致的定位误差，但是震源位置和波速耦合在一起反演时，对以微震震源位置为主要目的的反演是不利的［14］。
高明仕等［15］的研究结果表明，微震信号的能量沿传播距离呈乘幂关系衰减。为使微震信号能量衰减规律在震源定位上得到有效应用，笔者采用颗粒流程序［16-18］建立微震信号在花岗岩中传播的仿真模型，采用经验模态分解（empirical mode decomposition，EMD）处理信号，使用幂函数拟合、反函数变换建立震源距离-能量关系式进行震源反演，通过室内声发射试验考察反演精度，验证该数值计算模型的有效性，探究微震源定位新思路。
1 颗粒流数值模型构建
1.1 几何模型构建
建立尺寸为2.0 m×0.5 m的计算区域，花岗岩试样模型尺寸为1.8 m×0.3 m（如图1所示），颗粒半径取值为0.5 cm，试样密度设为2 900 kg/m3。
<G:\武汉工程大学\2023\第3期\李 贲-1.tif>
图1 花岗岩试样模型
Fig. 1 Model of granite sample
1.2 粘结模型选择及参数设定
模型构建所采用的粘结模型为平直节理模型。其原因为平直节理模型将圆形颗粒构造成多边形，可抑制接触破坏后颗粒的旋转，能够模拟出大压拉比［19］。根据张权等［20］得到的平直节理模型宏细观参数对应关系式［式（1）］，求解得到平直节理模型的细观参数如表1所示。
[E=1.03Ec-3.9kn/ks+14.6 v=0.146ln(kn/ks)+0.075 σf=1.724C1+80.1μb+0.97σc-13σt=0.692σc-1.181ln(kn/ks)+1.405C2=16.07+0.28C1-15.54μbφ=arcsin(1-40/(49.8+254μb+0.16C1))]
（1）
其中，E、υ、σf、σt、C2、φ分别为岩石材料宏观力学参数弹性模量、泊松比、抗压强度、抗拉强度、黏聚力、摩擦角，Ec、kn/ks、μb、C1、σc分别为颗粒流程序中平直节理模型的细观参数。
表1 平直节理模型细观参数计算结果
Tab. 1 Calculation results of fine view parameters of flat and straight nodal model
[Ec / MPa kn/ks μb C1 / MPa σc / MPa 88.90 5.30 0.26 88.70 18.60 ]
1.3 震源与监测点的布置
以试样模型左下角为原点建立直角坐标系，如图2所示。设置震源位置为（0.00 m，0.15 m），6个监测点位置为（0.28 m，0.15 m）、（0.56 m，0.15 m）、（0.84 m，0.15 m）、（1.12 m，0.15 m）、（1.40 m，0.15 m）、（1.68 m，0.15 m）。以实测振动波信号为震源模拟声发射，监测各测点波形。
<G:\武汉工程大学\2023\第3期\李 贲-2.tif>[震源][监测点1 监测点2 监测点3 监测点4 监测点5 监测点6]
图2 震源与监测点布置图
Fig. 2 Location of hypocenter and monitoring points
2 仿真结果分析
2.1 原始信号处理
由于监测到的信号中存在反射波以及其他噪声信号，因此需要对原始信号进行分解、滤波、降噪。图3所示为监测点1原始信号，将其进行EMD［21］得到固有模态函数（intrinsic mode function，IMF）分量，如图4所示。
通过分析可知，IMF0为高频率的噪声信号，IMF6、IMF7为零漂现象所导致的信号。将监测点1原始信号分解后的IMF信号进行筛选重构，作为后续能量计算的依据，如图5所示。
2.2 能量计算
要精确计算出各个监测点的微震信号能量较为困难，通常采用某一波形参数来表征测点能量大小，如最大峰值振动速度、峰值振动速度平方、持续时间内各采样点振幅的平方和（E）等，一般认为微震波持续时间内各采样点振幅平方和表示能量更准确［21］。笔者采用一段时间内E的最大值Emax作为表征能量EZ的值。Emax与EZ相差一个比例系数，只能表征能量的大小，并不能代表真正的能量值，它的单位不能用国际单位的焦耳（J）表示，应该为m2/s4。EZ的真实值应为：
[EZ=ZEmax] （2）
式中，Z=1 J·s4/m2。
根据式（2）对原始信号处理后得到的信号进行计算，在0~6 ms持续时间内各监测点的表征能量如图6所示。
2.3 能量衰减规律分析
从图6中可以看到，各测点的能量随时间的推移迅速衰减，并随着与震源距离的增加呈现递减趋势。通过能量计算，得到6个监测点的微震信号的表征能量值，如表2所示。有学者的研究表明幂函数、指数函数都能较好地描述能量随传播距离的增大而衰减的规律［22］，基于Levenberg-Marquardt最优算法，使用幂函数、指数函数将表2中的数据进行拟合，拟合结果如图7所示。通过对比，幂函数的拟合结果方差更小，验证了王进强等［22］所述的幂函数更适宜描述微震信号随传播距离的增大而呈现出前期衰减迅速，后期衰减缓慢的特性。
<G:\武汉工程大学\2023\第3期\李 贲-6.tif>[0 1 2 3 4 5 6
t / ms][30
24
18
12
6
0][能量 / kJ][测点1][测点2][测点3][测点4][测点5][测点6]
图6 不同监测点表征能量图
Fig. 6 Characterization energy map of different
monitoring points
表2 各个监测点的表征能量值
Tab. 2 Characterization energy values for each
monitoring point
[距离 / m 0.28 0.56 0.84 1.12 1.40 1.68 表征能量 / kJ 29.80 17.58 11.03 9.30 8.51 6.92 ]
为定量研究震源距离与表征能量之间的关系，通过幂函数拟合得到震源距离与表征能量的关系式：
[EZ=E0x-η] （3）
式（3）中：EZ为表征能量；x为测点与震源之间的距离；η=0.824 5，为能量衰减系数，与岩石的宏观力学参数有关；E0=10 461.104，与震源能量与衰减系数有关。由此，通过数值模拟及函数拟合，得到了微震信号能量在花岗岩中随传播距离的增大呈现出幂函数衰减的规律。
<G:\武汉工程大学\2023\第3期\李 贲-7.tif>
图7 能量-距离幂函数、指数函数拟合曲线
Fig. 7 Fitting curves of power function and exponential function of energy-distance
在实际工程应用中，一个区域内能量衰减系数[η]可通过一次爆破试验测得。而E0的值会随震源能量发生改变，因此需要通过2个传感器所得能量进行计算。1号、2号传感器与震源的距离分别为x1，x2，2个传感器的间距为[?x]，E0的计算推导公式如下：
[EZ1=E0x1-ηEZ2=E0x2-η] （4）
[EZ11-η=E01-ηx1EZ21-η=E01-ηx2] （5）
[EZ21-η-EZ11-η=E01-η(x2-x1)] （6）
[E0=(EZ21-η-EZ11-η?x)-η] （7）
2.4 距离反演
在原模型上重新设置5个监测点，其坐标分别为（0.10 m，0.15 m）、（0.50 m，0.15 m）、（0.90 m，0.15 m）、（1.30 m，0.15 m）、（1.70 m，0.15 m）。给震源相同的初始振动波，使振动信号在模型中传播，并对这5个监测点的信号进行记录。通过将数值模拟记录的结果进行筛选重构，根据能量计算理论得出0~6 ms持续时间内各监测点的距离-能量，如表3所示。
通过求式（3）能量-距离关系式的反函数可以得到距离-能量关系式：
[x=E0EZη] （8）
式（8）中：E0=10 461.104，η=0.824 5。
将表3中的能量值代入式（8）中，计算出各监测点与震源之间的距离，并与真实距离进行比较，采用相对误差来衡量计算误差。相对误差的计算公式为：
[δ=Lc-LrLr×100%] （9）
式（9）中：δ表示相对误差，Lc表示计算距离，Lr表示真实距离。所得结果如表4所示。
表4 计算结果与误差
Tab. 4 Calculation results and errors
[测点 计算距离 / m 真实距离 / m 相对误差 / % 1 0.100 13 0.10 0.129 2 0.498 51 0.50 0.297 3 0.897 28 0.90 0.302 4 1.295 00 1.30 0.381 5 1.706 40 1.70 0.379 ]
从表4的计算结果可以看到，各个测点的计算距离与实际距离的相对误差均在0.5%的范围之内，因此根据能量-距离关系式所推出的距离-能量关系在原模型上的不同监测点都具有很好的适用性。由式（3）可知，随着距离增加，能量不断衰减，能量变化对距离影响愈敏感，当距离增大到280 m以上，能量值会随之衰减到100 J以下，此时每0.5 J的能量计算误差就会造成0.6%的距离计算误差。而将原模型的能量计算精度控制在0.5 J以内是相当困难的，因此，在实际的工程应用中，需将原模型尺寸控制在280 m以内。
3 室内试验验证
3.1 试验过程
在实验室中，采用尺寸为0.60 m×0.15 m×0.15 m的花岗岩作为实验材料。首先使用石膏作为耦合剂，将1、4、5、6、7、8号传感器如图8所示分布在试件的上表面，间距均为0.10 m，以试件左下角为原点建立空间直角坐标系，震源位置为（0.000 m，0.075 m，0.075 m）。通过断铅模拟震源，使用江苏东华测试技术公司制造的DH5922D动态信号测试分析系统对各个传感器接收到的信号进行记录。
<G:\武汉工程大学\2023\第3期\李 贲-8.tif>[1#][4#][5#][6#][7#][8#]
图8 室内试验岩样
Fig. 8 Rock sample for laboratory tests
3.2 试验结果与分析
将室内声发射试验6个监测点记录到的数据通过分解、滤波、降噪及能量计算，得到0~6 ms持续时间内各监测点的表征能量图，如图9所示。
<G:\武汉工程大学\2023\第3期\李 贲-9.tif>[0 2 4 6
t / ms][210
168
126
84
42
0][能量 / kJ][测点1][测点2][测点3][测点4][测点5][测点6]
图9 室内试验不同监测点表征能量图
Fig. 9 Characterization energy diagram of different
monitoring points in laboratory tests
记录各个监测点的能量如表5所示。将表5中测点1~测点6的能量-距离按照式（3）函数进行拟合，得到室内试验的能量-距离关系式：
[EZ=43 837.780x-0.822 9] （10）
据式（10）可知室内花岗岩岩样的能量衰减系数为0.822 9，与数值试验测得的能量衰减系数相对误差为0.20%，证明了在颗粒流程序中所建立的数值模型能较准确地模拟微震信号在花岗岩中的能量衰减规律。
3.3 距离反演结果验证
为进一步探究距离反演方法的准确性，将震源位置改为（0.000 m，0.075 m，0.150 m），重复室内试验过程，计算得到5个监测点的能量，如表6所示。
将测点2和测点3的能量值代入式（7）计算得到E0=39 561.838，η=0.822 9，再代入式（8）得到室内试验距离-能量关系式：
[x=39 561.838EZ0.822 9] （11）
将测点4、测点5和测点6的能量值分别代入式（11），得到各测点的反演距离分别为0.229、0.399、0.501 m，根据式（9）得到与真实距离的相对误差分别为0.17%、0.16%、0.12%，均在0.5%以内，证明该方法可以实现距离精准反演。
表6 室内试验各测点能量反演结果
Tab. 6 Energy inversion results for each measurement
point in laboratory tests
[测点 真实距离 / m 表征能量 / kJ 2 0.10 263.142 3 0.20 148.753 4 0.30 106.615 5 0.40 841.363 6 0.50 69.864 ]
4 结 论
（1）应用平直节理模型构建花岗岩仿真模型模拟声发射试验，通过室内声发射试验与模拟声发射试验对比，计算得到能量衰减系数误差率为0.2%，证明该模型能有效模拟出微震信号在花岗岩中的衰减规律。
（2）采用EMD处理原始信号，使用幂函数拟合微震信号能量衰减趋势，构建震源距离反演数学模型，室内试验距离反演结果与真实距离误差均在0.5%以内，证明了基于能量衰减的微震源距离反演的可行性。
（3）所建立的花岗岩震动传播颗粒流仿真模型以及距离反演数学模型结构简单，计算方便，在一定距离范围内反演精度较高，可为震源定位提供一种新思路。