二维非测速条件下声发射震源定位方法数值验证

吴顺川，张光，张诗淮，郭沛，储超群

(1.昆明理工大学国土资源工程学院，云南昆明，650093；2.北京科技大学金属矿山高效开采与安全教育部重点实验室，北京，100083)

岩石类材料在受内力或外力作用下发生变形或破裂，在破裂过程中会激发应力波。利用波在传播过程中的时间信息，由时空坐标关系可定位震源。吴顺川等[1]提出了一种二维非测速条件下声发射震源定位方法，并采用花岗岩平板验证其适用性。声发射定位技术在结构检测、矿山破岩、油井致裂以及隧道掘进等诸多方面得到了广泛应用。定位方法是确定岩石内部缺陷以及破裂损伤位置的关键，国内外已有大量学者研究了定位方法。传统定位方法包括GEIGER 法[2]、单纯形法[3]和网格搜索法[4]等。TOBIAS[5]提出了三角时差定位法；CIAMPA 等[6]提出了无需测量波速的震源定位方法；DONG 等[7-9]提出无需测速的微地震/声发射定位方法，并采用复杂结构模型和原位微震监测验证了定位效果；KUNDU 等[10-11]总结了不同定位方法在各向同性与各向异性材料中的定位差异，并提出了一种特定传感器布设方式的无需测速震源定位方法；YIN等[12]基于新三角时差定位方法原理，提出了一种新的无需测速震源定位方法，其传感器布设方式为“Z”形。在定位方法的研究基础上，众多学者研究定位精度，LI等[13]采用花岗岩和大理岩研究速度偏差对线定位和平面定位精度的影响；许江等[14]采用砂岩试样研究了岩石尺寸效应与声发射定位精度影响因素。由于应力波的传播规律极其复杂，室内试验中波的传播衰减规律难以满足理论条件，数值模拟对应力波进行研究具有独特优势[15-16]。目前，模拟应力波方法主要包含连续介质理论和非连续介质理论。岩体包含众多断层、节理和裂隙等地质结构，其主要特性为非连续性、非均匀性和各向异性。采用连续介质理论难以准确、客观地描述波在岩体中的传播规律。CUNDALL等[17]提出的离散元法能够有效模拟岩体的力学特性，反映其非连续性。离散元法适用于计算岩体等非连续性介质的动力学问题，在岩土工程领域应用广泛，国内外众多学者采用离散元模型研究了岩石类材料应力波传播问题。TRENT 等[18]采用离散元法，从细观尺度上研究了介质的波动过程；HAZZARD等[19]基于颗粒流理论构建砂岩模型，研究了荷载损伤下的波速演化规律；NABIPOUR[20]应用颗粒流理论建立了六角形排布和方形排布颗粒模型，研究了颗粒施加应力波传播特性，验证了颗粒流模型中应力波传播的合理性；TOOMEY 等[21]应用颗粒流理论构建六角形排布颗粒模型，研究了应力波在大尺度规则排列颗粒模型中的传播规律；SADD等[22]基于离散元理论研究了枝矢量、接触法向和轴向矢量分布与波速和振幅衰减的关系；张国凯等[23]采用颗粒流程序建立9种不同各向异性模型，研究了波的传播规律，揭示孔隙率、枝矢量、配位数张量和刚度张量等细观参数对波的传播和衰减规律；张诗淮等[24]基于颗粒流理论研究应力波在大尺度一维和二维模型中传播，在震源激发方式、颗粒黏结模型、数值弥散和边界条件等方面研究了P波在数值模型中的传播规律；郭易圆等[25]采用面-面接触模型模拟纵波在一维岩体中的传播，得出了阻尼比、软弱夹层以及节理对波传播规律的影响；徐小敏等[26]针对砂土等颗粒材料建立离散元模型，分析了激发频率、激发幅值、激发源和接收源尺寸和阻尼等因素对剪切波速的影响；吴顺川等[27]基于颗粒流理论构建层状岩体模型，在震源定位计算中考虑了层状结构对应力波传播的影响。前人多基于颗粒流理论采用均匀颗粒构建一维和二维模型，在相对简易模型中研究应力波的传播规律，且未进行定量分析。传统定位方法多适用于各向同性介质，对于岩石类非连续性介质，由于内部结构较为复杂，采用传统单一波速模型进行定位会产生误差。吴顺川等[1]提出的非测速条件下声发射震源定位方法无需预先测速，避免了波速误差对定位结果的影响。室内试验会存在较多的外界干扰因素和人为因素引起的试验误差，数值试验条件下应力波的传播规律更接近于理论条件下的传播规律，采用离散元方法研究应力波具有独特优势。本文基于颗粒流程序构建平板模型，在平板上布设震源点和监测点，采用互相关理论分析监测点的波形衰减规律，确定监测点时延，依据二维非测速震源定位方法确定震源位置并进行误差分析，并基于平板模型探究波动理论在平节理模型(FJM)的适用性。

1 二维非测速定位方法

1.1 二维非测速条件下震源定位方法

根据文献[1]可知，定位监测所用的簇内传感器布设方式可为任意三角形，如图1所示。当在平板上布设多簇传感器时，可对震源进行定位，图1中在平板上布设4簇传感器，传感器簇编号为Ci(i=1，2，3，4)，簇内传感器编号为Si-j(j=1，2，3)，每簇传感器将Si-1设定为参考传感器，计算簇内其他2 个传感器与该参考传感器的时延分别为ti-21和ti-31。

由文献[1]可知，震源A与参考传感器Si-1连线的坐标方位角θi以及超声波信号在该路径的传播速度c(θi)计算结果如下：

图1 四簇传感器布设示意图[1]Fig.1 Four sets of sensors on a plate[1]

式中：α为Si-1Si-2与x轴正方向夹角；β为Si-1Si-3与x轴正方向夹角；a为Si-1Si-2的距离；b为Si-1Si-3的距离。

由式(1)和(2)可知，坐标方位角θi以及该传播路径上波速c(θi)的计算结果准确性取决于ti-21和ti-31。依据式(1)求出一簇传感器的震源坐标方位角，2簇传感器联合求解即可确定一个震源坐标。

1.2 时延测量和波形互相关系数分析

数值试验中采用正弦波和雷克子波作为震源输入波形，在平板中设置不同监测点对震源施加的波动信息进行监测，通过计算同一簇监测点的时延确定震源坐标，确定时延的方法主要包括初达波法、峰-峰值法和互相关分析法[26](crosscorrelation，CCR)。采用初达波法和峰-峰值法确定监测点间时延，测量准确性会受到波形扭曲和弥散的影响；采用互相关方法不仅可以确定各监测点时延，还可利用互相关系数评估不同监测点波的传播质量及衰减规律。

对所监测得到的2组波形信号x(t)和y(t)进行处理，其互相关函数[28-29]Rxy(τ)定义为

式中：T为波形时长；τ为时延。

在对其进行处理中，波形时长会选取固定值，互相关函数估计值(τ)为

当τ=τ0时，(τ)函数取最大值，可知波形信号x(t)和y(t)在此处互相关程度最高，其时延为τ0。

对随机信号进行相关性分析时，如果分析整个时间序列信号，较大样本量会造成差异增大，导致互相关系数计算不准确，可采用时间序列延迟相关分析法提高相关系数计算的准确性[30]。信号x(t)={x1，…，xt，…，xn}和y(t)={y1，…，yt，…，yn}在时延τ处相关系数r(τ)计算公式为

相关系数r(τ)的取值范围为[-1，1]，|r(τ)|越接近于1 时，表明信号x(t)和y(t)线性关系越密切；|r(τ)|接近于0 时，表明信号x(t)和y(t)线性关系越弱。

2 震源定位颗粒流模拟

2.1 离散元模型构建

文献[1]采用花岗岩平板进行震源定位试验，并对非测速条件下声发射震源定位方法进行了验证，为便于与室内试验进行对比，基于颗粒流程序建立的平节理模型[31-32](flat-joint model,FJM)与室内试验所用花岗岩平板宏观力学参数均一致。花岗岩力学性质如表1所示[33]。在PFC3D中生成直径为50 mm和高为100 mm的圆柱试样进行单轴压缩试验与直接拉伸试验，调整模型细观参数，将模型的宏观力学性质与实验室结果相匹配，当计算结果与实验室结果相近时，可将该组细观参数应用于模型计算。

表1 花岗岩力学性质[33]Table 1 Mechanical properties of granite[33]

本文采用FJM 模型构建花岗岩平板，验证定位方法，表2所示为平节理模型细观参数，当采用表2所示细观参数匹配结果时，所得力学性质与室内试验中花岗岩试样基本一致(见表1)。采用该细观参数构建长为500 mm、宽为500 mm、厚为18 mm 的平板模型(如图2所示)，模型中颗粒最小直径为3.0 mm，最大直径为4.5 mm，包含颗粒总数为101 764个。在平板上分别布设震源点和监测点颗粒，对震源颗粒施加脉冲信号，采集监测点应力波信号并进行分析，研究波形衰减以及震源定位结果。

在数值试验中，震源布置位置和监测点布设方式与室内试验布设方式均一致，监测点采用图3所示的I，II 和III 类3 种布设方式，通过数值试验验证非测速条件下声发射震源定位方法的适用性与准确性。

2.2 震源激发方式

对震源颗粒以集中力形式分别施加正弦波(Sine wave)和雷克子波(Ricker wave)脉冲信号，模拟声发射信号。脉冲激发频率为100 kHz[19]，为保证模型的稳定，激发幅值不宜过大，在该模型中设定振幅为1.0×10-8N，每一步运行时长Δt为5× 10-8s。在震源点(A)的z方向上分别施加正弦波和雷克子波脉冲，在平板上布设监测点记录应力波信息，分析其传播规律，计算各监测点时延进行定位计算。

表2 平节理模型细观参数Table 2 Micro-parameters used in FJM

图2 颗粒流平板模型Fig.2 Plate model in particle flow code

图3 监测点布设方式示意图[1]Fig.3 Different sensor arrangements on a plate[1]

2.3 应力波振幅衰减及互相关系数分析

在平板上布设多个监测点，研究正弦波和雷克子波在不同震中距的振幅衰减以及不同监测点间的相关性规律。在震源点A与监测点S3-1(此为第10 衰减监测点，Satt-10)之间等间距布设9 个监测点(监测点编号为Satt-i，i=1～9)，记录各监测点z方向位移变化信息，如图4所示。

图5所示为AS3-1传播路径上监测点应力波振幅与互相关系数变化规律。为便于对比分析，以第1 个衰减监测点(Satt-1)为基准对各监测点振幅归一化处理，互相关系数变化规律为各监测点波形与第1 个监测点波形进行互相关处理所得。由图5可见：正弦波和雷克子波初峰振幅随着与震源距离增大而不断降低，且在较近距离范围内会急剧下降，第3个监测点(Satt-3)振幅降低为第1个监测点的50%，第4～10 监测点振幅逐渐减小，最终在第10监测点减小至第1监测点的20%以下。

图4 平板模型监测点分布示意图Fig.4 receivers placed on the plate

对各监测点波形与第1个监测点波形进行互相关处理得到各监测点的互相关系数。正弦波在第2～6 监测点与第1 监测点的波形互相关系数缓慢降低，第7～10监测点互相关系数骤降，在第2～10监测点互相关系均大于0.7，可知其与第1 监测点波形高度相关。雷克子波在第6～8监测点比第2～5监测点的互相关系数的降低速率快，第9 和10 监测点的互相关系数与第8监测点相比有所增大，在第2～10 监测点互相关系数值均大于0.5，与第1 监测点波形显著相关。

图5 AS3-1传播路径上监测点应力波振幅与互相关系数变化规律Fig.5 Variations of the stress wave amplitude and crosscorrelation coefficient of the receivers on the AS3-1 propagation path

对于2种不同震源激发方式，振幅的衰减规律存在差异。随监测点与震源距离增大，正弦波和雷克子波的振幅均在不断减小，但雷克子波振幅衰减快于正弦波振幅衰减，在远离震源位置的各监测点上正弦波互相关系数大于雷克子波互相关系数。

2.4 震源定位计算及误差分析

2.4.1 时延测量

由式(1)可知：声发射震源定位计算中，定位结果的准确性主要取决于不同监测点间时延计算。图6所示为一次试验中C1簇监测点(S1-1，S1-2和S1-3)所记录的应力波信息，3个监测点所采集的波形高度相似，可采用互相关方法对应力波进行处理计算监测点间时延。图7所示为S1-1，S1-2和S1-3监测点所采集的应力波进行互相关处理所得结果，S1-2和S1-3分别对S1-1监测点应力波进行计算可得各监测点间时延，互相关函数峰值处横坐标即为监测点间时延。

图6 不同监测点z方向位移波形Fig.6 Waveforms of z-direction displacement in monitoring points

2.4.2 震源定位结果

平板模型上布设4 簇监测点(C1，C2，C3和C4)分别采集1 号～4 号震源应力波波形信息，依据式(4)对监测点波形信号与参考监测点(Si-1)进行互相关计算，得到时延并根据式(1)求解震源定位结果。图8所示为采用花岗岩平板进行定位试验中一次定位结果计算示意图，在试验中，传感器布设方式为II 类布设，对4 号震源进行定位验证[1]。在该定位试验中，每2簇传感器定位直线交点可确定1个计算定位点，采用4簇传感器进行震源定位，可计算得到6个定位结果。

图7 波形信号互相关处理结果Fig.7 Cross-correlation of the acoustic emission waveform

图8 花岗岩平板4号震源定位计算示意图Fig.8 Location calculation of 4号source on granite plate

图9和图10分别为对震源施加正弦波脉冲和雷克子波脉冲模拟声发射信号进行试验所得定位结果，其中，监测点布设方式采用图3所示的3种布设方式(I，II 和III 类)。分析定位结果精度时，引入较优点概念[14]，即定位结果与实际震源距离相对误差小于20%的点为较优点，定位结果与实际震源距离相对误差和定位直线与实际震源距离相对误差统计见表3。

由图9可见：1～4 号震源定位试验中定位结果大部分集中在实际震源附近，落在较优点范围内，但仍有部分定位结果误差较大。采用III 类布设方式时，1 号和4 号震源定位试验中C2-C3簇监测点定位结果误差较大，2 号震源定位试验中存在1 个较大定位偏离点(图9(b)虚线椭圆标记处)。2 号震源与C1簇和C4簇监测点连线线段之间距离较小，C1簇传感器与2号震源连线的坐标方位角与2号震源与C4簇传感器连线的坐标方位角相近，在定位计算所得的坐标方位角的微小误差会导致定位点的较大误差，该误差产生原因见文献[1]。

由图10可见：1 号震源采用II 类布设方式和4号震源采用III 类布设方式所得定位结果偏离实际震源较远，超出较优点范围(即相对误差大于20%)，其余定位结果较为集中，均在较优点范围内。

由图9和10所示的震源定位结果和表3中定位结果相对误差统计可知，2 种震源激发方式(正弦波和雷克子波)均较好地验证了该震源定位方法在非连续介质的适用性与准确性。

表4所示为文献[1]中采用花岗岩平板进行室内试验所得定位结果相对误差统计结果，表5所示为数值模拟定位结果相对误差统计结果。在数值试验定位结果中，采用正弦波激发震源，3种监测点布设方式(I 类、II 类和III 类)所得震源定位结果较优点比例分别为100.0%，100.0%和87.5%；采用雷克子波激发震源，3 种监测点布设方式(I 类、II类和III 类)所得震源定位结果较优点比例分别为100.0%，95.8%和95.8%。与室内试验定位结果相比，数值试验定位结果相对误差较小，定位准确性较高。在定位试验中，室内试验中影响因素更为复杂，花岗岩平板中存在的固有缺陷会对定位结果造成误差，而数值试验干扰因素较少，应力波传播更接近理论条件，故数值试验定位结果较优点比例明显高于室内试验较优点比例。

图9 正弦波激发方式震源定位结果Fig.9 Source location results by the type of sine wave excitation

图10 雷克子波激发方式震源定位结果Fig.10 Source location results by the type of ricker wave excitation

数值试验中颗粒粒径及位置为随机生成，颗粒在空间上非均匀排布。定位计算时，监测颗粒z方向坐标差异会对定位结果带来一定误差，且该定位误差无法通过计算消除。由于该二维定位方法并未考虑监测点z方向影响，定位方法理论计算中仅包含x及y方向，监测颗粒在实际z方向上会存在差异，导致时延测量出现误差。

3 P波辐射花样定量分析

文献[24]定性验证了应力波在接触模型中的合理性，由于接触模型无法有效反映岩石力学特性，基于FJM 模型，进一步定量对比应力波辐射花样与理论解的差异。

模型中弹性波由点震源激发，根据弹性波动力学理论，其三维波动矢量方程及其格林函数解[24,34]如下：

表3 数值试验平板定位结果相对误差统计Table 3 Relative error of source location on the plate from the numerical simulation

表4 室内试验平板定位结果较优点统计[1]Table 4 The better points counts from the laboratory test[1]

式中：ψ(r，t)为质点位移；r为辐射距离矢量；t为时间；c为波速；F(r，t)为源项；uij为位移分量；i为位移方向(i=x，y和z)；j为作用力方向(j=x，y和z)；γi和γj分别为震源与监测点的连线矢量与坐标轴的方向余弦；r为震源与监测点距离；T(t)为震源激励函数；vP为P波波速；δij为狄拉克函数；vS为S波波速；ρ为密度。

表5 数值试验平板定位结果较优点统计Table 5 The better points counts from the numerical simulation

在平板模型震源点以力形式在z方向施加脉冲波形。式(7)中为近场项，为P 波远场项，为S 波远场项。在实践中震源机制研究中主要取决于位移场远场项，因此，在计算中可忽略近场项[35-36]，并将远场项中的常数项取1，则有

式中：方向余弦为γx=sinθcosφ，γy=sinθsinφ，γz=cosθ。

图11所示为P波和S波三维辐射花样。由图11可知：P波在平行于施加力方向上(z方向)，位移最大，而在垂直于施加力方向上(x和y方向)，位移为0；S 波在垂直于施加力方向上(x和y方向)，位移最大，而在平行于施加力方向上(z方向)，位移为0。在x-z二维坐标系中，P波和S波辐射花样如图12所示，其位移分量如下：

图11 z方向集中力激发震源P波和S波震动位移场远场项三维空间辐射花样分布Fig.11 Three-dimensional representations of the far-field P-and S-wave radiation patterns generated by a point source in the z direction

P波和S波的振幅表达式为

图12 z方向集中力激发震源P波和S波震动位移场远场项x-z平面辐射花样分布[24]Fig.12 Two-dimensional representations of the far-field P-wave and S-wave radiation patterns generated by a point source in the z direction[24]

图13所示为平板模型运行400 时步(2×10-5s)颗粒位移矢量图。在震源z方向施加雷克子波，观察震源周边颗粒位移变化，颗粒位移振动均为z方向。图13中A范围为P波波峰，B范围为波谷，波形由震源点向四周传播时，颗粒位移方向均与施加力方向(z方向)一致，在x和y方向上无位移。应力波在x-y平面上以圆周形式向外传播，平板边界为自由边界，未设置吸收边界，C范围为平板边界，故该范围处颗粒z方向位移明显大于平板内部颗粒位移，应力波传播至平板边界会进行反射，入射波与反射波在此处位移矢量叠加，该处位移明显大于其他传播方向位移，该模拟结果与张诗淮等[24]研究的边界条件对应力波传播振幅影响规律一致。

图13 2×10-5 s时刻平板模型颗粒位移矢量图Fig.13 Displacementvectors intheplatemodelat 2×10-5s

图14所示为震源A与衰减监测点Satt-5的位移监测波形信息图，为便于比较，将震源点与监测点位移波形信息振幅进行归一化处理，其中z方向位移分量显著大于x和y方向位移分量，与图11和图12中辐射花样相符。P波传播位移方向均平行于施加力方向(z方向)，在垂直于施加力方向上无位移。

在AS3-1传播路径上布设的各监测点与震源之间距离(r)和波长(λ)之间关系均满足r>>λ，可知在该研究中位移场主要取决于远场项[36]，忽略近场项影响，由式(7)可得，各监测点P 波在z方向位移如下：

式中：vP为由10 个监测点所得的P 波速度，取4 335.3 m/s，ρ=2 800 kg/m3，如表2所示。震源激励函数T(t)如下：

图14 震源点与监测点(Satt-5)x，y和z方向振幅对比Fig.14 The amplitude of x,y,and z direction displacement from the source and receiver(Satt-5)

式中：A0取1.0×10-8N，f取100 kHz。

若以集中力作为激发震源，式(13)的理论解如图15(a)所示。由图15(a)可见：在距离震源较近时，位移振幅随r增大迅速衰减；在距离震源较远时，衰减速率显著变缓，两者呈反比关系。进一步，将AS3-1传播路径上布设的监测点位移振幅与式(13)理论解进行对比，如图15(b)所示，监测点实际振幅与理论解基本一致，证明了波动理论在FJM模型中的适用性与准确性。

图15 z方向位移振幅衰减性质Fig.15 Attenuation property of z-direction displacement

4 结论

1)采用互相关技术处理声发射信号，可准确得到相邻监测点时延，根据时延反演震源定位结果；随着监测点与震源点间距离增大，各监测点振幅及互相关系数均不断减小。

2)计算定位点均分布在实际震源周边，相对误差较小，验证了二维非测速下震源定位方法在非连续介质中的适用性。

3)在正弦波激发方式下，等腰直角三角形，一般直角三角形和等边三角形3种监测点布设方式定位结果较优点比例分别为100.0%，100.0%和87.5%；在雷克子波激发方式下，3 种监测点布设方式定位结果较优点比例分别为100.0%，95.8%和95.8%。

4)应力波信号在平板模型x-y平面上以圆周形式向外传播，各监测点P 波位移在z方向上最大，在x方向和y方向无位移。平板各监测点实际振幅与理论解基本一致，证明了波动理论在FJM 模型中的适用性与准确性。