岩石在外界应力作用下，其内部将产生局部弹塑性能集中现象，当能量积聚到某一临界值之后，就会引起岩体微裂隙的产生与扩展，微裂隙的产生与扩展伴随着弹性波或者应力波的释放并在周围岩体内快速传播，这种弹性波在地质上称为微震(microseismic/acoustic emission)(Xu et al.，2011).微震现象是由美国矿业局Obert(1977)发现，后来Kaiser(1950)对该现象命名为Kaiser效应.微震监测技术是通过在监测区域岩体布置检波器或者传感器来收集微震信号，然后分析该信号以圈定岩体损伤区域，从而为区域岩体稳定性评价提供有效依据.如今，微震监测分析技术作为工程建设安全预测预警方法被广泛应用到矿山开采(姜福兴等，2006；陆菜平等，2010)、石油开采(宋维琪等，2015)、水利水电工程(Xu et al.，2011；张伯虎等，2012)等行业.

为了使微震监测技术更好地发挥预测预警作用，满足工程需求，定位精度的提高至关重要.微震定位问题属于地球物理反演问题，它根据一次微震事件在各传感器的震相到时来反演微震的发生时刻和空间位置(Spencer and Gubbins，1980；Pavlis and Booker，1980；刘福田，1984；刘杰等，2003；高永涛等，2015).但是，在实际工程应用中，由于区域岩体是非均匀各向异性介质，同时岩体中存在不同走向和倾角的节理断层，微震震源发出的应力波在传播过程中能量、频率和波速改变，这都给微震定位精度带来了较大的影响.为了实现高精度微震定位，需要采用合理的传感器空间布置方案、有效的微震波形识别技术以及合适的微震定位方法.

微震定位方法主要有几何作图法、相对定位法、空间域定位法、线性定位法和非线性定位法.几何作图法是一种直观的定位方法，依据走时方程，通过二维或者三维作图确定震源位置(廉超等，2006).相对定位法又称为主事件定位法，是根据一个震源坐标较精确的主事件来计算其附近其他微震事件的坐标(Spence，1980).空间域定位法是由距离残差代替到时残差的一种定位方法.微震震相到时和震源参数的关系是一种非线性关系.线性定位法是通过求导、级数展开方法将该非线性问题转化为线性问题的定位方法，例如Inglada(1928)法和Usbm(1970)法.这类方法在到时拾取精度不高的情况下定位效果较差.非线性定位法是利用迭代或优化求解方式的定位方法统称，例如Thurber(1985)用非线性牛顿迭代法进行定位.随着计算机计算能力不断发展，不同的非线性优化算法工具被应用到定位中，用于算法优化的目标函数也出现了不同的形式.赵珠和曾融生(1987)在定位中采用了单纯形法.但是单纯形法采用单向搜索，形式简单，容易陷入局部最优解.针对这一缺点，全局最优化算法被应用到定位中，王洪体等(2006)介绍了一种基于浮点遗传算法的定位方法，具有稳定收敛、精度高、时效性强的特点；张华等(2014)基于遗传算法对厂矿区内的人工爆破进行定位，得到了很好的效果；针对速度模型给不准和联合定位法震源位置、发震时刻和微震传播速度相互关联，解不唯一的问题，陈炳瑞等(2009)提出微震震源定位分层处理方法，并采用粒子群算法进行优化求解；董陇军等(2011)针对单一速度模型提出了无需预先测速的微震震源定位方法，利用非线性拟合工具进行求解；李健等(2014)针对单一速度模型，提出了无需测速的定位方法，采用单纯形法进行最优化求解.上述定位算法多是建立在单一速度模型的基础上.虽然实际岩体中不同区域和方向上的波速值不尽相同，但是单一速度模型因不考虑太多未知量和关系，计算稳定性强，是最常见的一种速度模型.

近几年，一些专家和学者在寻找适应各种区域岩体微震定位的速度模型上做出了不同探索.一般情形下，某个微震定位算法速度模型的优劣取决于两点：(1)是否可以较好地考虑区域岩体的波速特性，且能够给出恰当的应力波传播路径；(2)是否简洁，从而利于算法稳定实现(Feng et al.，2015).射线追踪理论能够找到复杂岩体中应力波传播路径，从而逐渐被应用到波速差异大、带空洞的区域岩体微震定位中(Bai and Greenhalgh，2006；Bai et al.，2010；Collins et al.，2014；Hosseini et al.，2015).但是，复杂岩体的速度模型建立及射线追踪的实现对计算机要求相当高，计算耗时大且不利于算法稳定地实现.为此，一些专家和学者针对不同波速分布特征的区域岩体分别对应力波传播特性作了相应简化处理，在提高定位精度的同时保证了计算的稳定性.其中，Feng等(2015)针对隧洞开挖工程中围岩波速特征，提出按隧洞断面对传感器进行波速值分组的定位速度模型，考虑了开挖隧洞中各断面波速值差异，同时采用粒子群算法对所建立的目标函数实现稳定求解；巩思园等(2012)针对煤矿上覆岩层多以层状形式赋存，P波在垂向和水平向的速度及传播路径差异大的情况，提出建立“异向波速模型”，模型求解选用遗传算法与CMAES(Covariance Matrix Adaptation Evolution Strategy，协方差矩阵自适应进化策略)算法结合的混合算法，经验证定位效果优于单一速度模型.

本文综合分析各种非线性定位方法，在“异向波速模型”基础之上，将垂向上的应力波按各岩层波速值大小进行分段处理，提出一种简化的分层速度定位模型SV.该定位模型在考虑了不同岩层应力波波速值大小不同的基础上，对复杂的层间应力波传播路径进行了简化处理，采用具有全局寻优特性的遗传算法对建立的分层速度目标函数稳定求解.相比单一速度定位模型，在波速分层的区域岩体微震定位中有利于定位精度的提高.

微震发生时向周围岩体释放纵波(P波)和横波(S波)，前者传播速度比后者快.当传播距离足够远时，在波形图中能够将P波与S波起跳时刻区分开来.然而，对于微震监测区域来说，由于监测区域不大，往往S波的起跳点叠加在P波中，要精确拾取S波初至时刻很困难.因此，通常仅拾取P波初至时刻用于定位.

对于单一速度模型，第i个传感器接收P波初至时刻t_i与震源参数(x₀，y₀，z₀，t₀)之间的关系由(1)式表述：

(1)

式中：l_i为震源与第i个传感器之间的空间距离；(x_i，y_i，z_i)为第i传感器空间坐标，m是监测区域内接收到P波的传感器数量；V为监测区域P波波速值.

理论上，微震事件发生时刻与走时之和t₀+Δt_i应同t_i相等，即

(2)

但是，在实际工程应用中，由于受到监测仪器以及人为因素的影响，微震信号起跳时刻的拾取存在误差，则ξ_i值不为零.在这样的情况下，为了实现微震震源定位，需使用非线性优化工具进行最优化求解来逼近真值.通常单一速度定位模型采用的非线性优化目标函数有如下三种形式：

因变量为到时的目标函数：

(3a)

因变量为到时差的目标函数：

(3b)

因变量为到时差商的目标函数：

或

(3c)

运用非线性优化方法，得到使上述目标函数尽量小的震源参数逼近解作为微震震源参数的最终计算结果.其中，第一种目标函数(3a)能将震源发生时刻一起反演；第二、三种目标函数(3b)、(3c)不能反演震源发生时刻；第一、二种目标函数(3a)、(3b)含有岩体波速参数，在优化求解过程中可以与震源参数一起反演，从而实现均匀波速模型下无需预先测速的微震定位(董陇军等，2011)；第三种目标函数(3c)不含岩体波速参数，从而能够实现均匀波速模型下的无需测速微震定位(李健等，2014).

在实际工程中，微震监测区域岩体并不是均匀介质，不同区域的应力波波速不同.因此，应力波射线是按最小走时路径传播，而并非几何最短路径.在许多工程中，岩体呈层状分布，例如层状矿层、层状边坡.

在层状煤矿微震监测中，震源应力波传播到矿区地面的检波器比到井下巷道经过的岩体介质更为复杂，导致P波在垂向上速度变化很大.针对煤矿中这种层状赋存和离层带的特点，巩思园等(2012)提出对安装在井下巷道和矿区地面的检波器进行区别对待，采用两种速度建立异向波速模型：针对井下的检波器，岩体波速较为均匀，P波传播最小走时路径近似为最短路径传播，速度为V_i=V_pu；针对矿区地表的检波器，将实际最小走时路径(折线)简化为最短路径，速度为V_i=V_pg.本文在此异向波速模型基础之上，将垂向上的应力波按各岩层波速值大小进行分段处理，提出一种简化的分层速度定位模型SV，以期使微震定位精度得到提升.

如图 1所示，假定微震监测区域岩体由四种不同波速的平行岩层组成，自下而上依次编号为①、②、③和④类，对应的层速度值分别表示为V₁、V₂、V₃和V₄.实际地震波在不同介质间传播时是遵从费曼原理按折线传播，巩思园等(2012)将矿区地表检波器接收到微震应力波(层面垂向上的应力波)的传播路径近似为直线，并将垂向上应力波波速值视作区别于水平向的均值.本文考虑到垂向应力波路径穿过的各岩层波速值不同这一特征，将对应在各层岩体中的走时由相应的空间距离以及波速值表示出来，然后求和得到应力波总走时.具体地，以3号传感器接收到的在①类波速岩层发生的微震事件应力波为例，推导其走时关系式如下：

图 1

Fig. 1

图 1 速度分层岩体中应力波传播特征 Fig. 1 Propagation characteristics of stress waves in rock mass with stratified velocities

通过震源与编号3传感器的空间直线PS可用参数式方程组表达为：

(4)

式中：k为直线参数；(x₀，y₀，z₀)为震源(P)的空间坐标；(x₃，y₃，z₃)为编号3传感器(S)的空间坐标.

由于①类和②类波速岩体交界面、②类和③类波速岩体交界面以及③类和④类波速岩体交界面互相平行，其空间平面一般解析式分别表示为：

(5)

式中：A，B，C是岩体交界面法向向量的三个分量值.D_1，2、D_2，3、D_3，4分别为相应交界面一般式方程中常数项.

经过推导，编号3传感器接收到发生在①类波速岩层中的微震事件应力波走时Δt₃：

(6)

式中，α₃和β₃定义为空间距离参数，由(7)式计算得到：

(7)

在实际工程中，为了达到良好的监测效果，监测网络呈网状布置，各传感器位于不同波速的岩层中.这样，同一个微震事件的应力波传播到各个传感器所经过的岩层个数和岩层类别不一样.由前文可以类比得到，某微震事件应力波传播到编号i传感器的走时Δt_i可以由(8)式表示为：

(8)

其中，空间距离参数α_i和β_i为：

(9)

式(8)和(9)中：从震源发出的应力波传播到编号i传感器所经过的波速层依次编号为1，2，…，n；层速度依次为V_i1，V_i2，…，V_in；应力波所穿过不同波速层交界面的数学解析式中常数项依次为D_i(1，2)，D_i(2，3)，…，D_i(n-1，n).

特别的，如果震源与编号i传感器在同一波速层，则对应的Δt_i为：

(10)

基于前文速度分层区域岩体的应力波走时关系，可以构建因变量为到时和因变量为到时差两种形式的分层速度定位目标函数如下：

因变量为到时的分层速度定位目标函数：

(11a)

因变量为到时差的分层速度定位目标函数：

(11b)

上述等式右边含有参数：传感器空间坐标(x_i，y_i，z_i)、分层模型空间几何参数(A，B，C，D)、微震信号初至时刻t_i、层速度值V_i、微震震源参数(x₀，y₀，z₀，t₀).本文采用因变量为到时差的目标函数(11b).

通过上述推导，可以看出，如果微震震源发生在不同的岩层，那么对应目标函数的具体表达形式也不同.因此，在利用遗传算法求解微震空间坐标的过程中，需要分波速层进行搜索，具体波速层η对应的搜索可行域Ω_η按如下线性规划不等式加以限定：

(12)

式中，

其中：η为波速层编号；A_η，B_η，C_η，D_η为边界面解析式方程参数；γ为边界面个数.

在分层速度定位模型中，采用遗传算法进行分层速度定位目标函数最优化求解，属于目标函数最小化问题.但是，在遗传算法中，种群是向适应度值高的方向进化.因此，需要建立目标函数与适应度函数值之间的映射关系.本文采用的适应度函数H与目标函数f的映射关系如下：

(13)

式中：c_max是当前所有或最近k代种群中f的最大值，并且c_max随遗传代数不断变化.图 2是目标函数的遗传算法优化流程.在优化求解过程中，当遗传算法达到以下两个条件之一则终止算法：达到最大遗传代数N(条件A)；遗传代种群距离小于阈值Δd(条件B).对每个波速层求解时，如果因达到条件A而终止算法，则认为该波速层的震源优化解不存在；如果因达到条件B而终止，则认为该波速层的解存在，同时如果终止种群的最优个体适应度值和平均适应度值分别大于其阈值H和h，则认为解是可行的.终止种群中的最优个体作为该波速层的震源优化解，同时记录其对应的目标函数值.从解存在的各个波速层中，选取目标函数值最小者对应的震源优化解作为最终定位结果，同时其目标函数值作为误差值来评价该定位结果的可靠性.

图 2

Fig. 2

图 2 遗传算法非线性目标函数优化流程图 Fig. 2 Flow diagram of the nonlinear objective function optimization using GA

在此基础上，如图 3所示，本文构建分层速度定位模型分为三个模块：参数准备、层速度反演、微震定位：

图 3

Fig. 3

图 3 分层速度定位模型流程图 Fig. 3 Flow chart of stratified velocity location model

(1) 模块Ⅰ：参数准备.确定监测区域岩体范围，建立坐标系和分层速度几何模型.然后划定波速层，并计算不同波速层交界面解析式方程；确定各岩层遗传算法搜索可行域Ω；根据监测区域岩体范围，优化布置传感器网络，得到传感器坐标；综合计算速度和精度要求，设定遗传算法优化终止参数.

(2) 模块Ⅱ：层速度反演.收集现场岩爆或施工爆破事件数据，拾取相应各通道传感器采集到的波形起跳时刻.根据岩爆或者施工爆破事件空间坐标(x_b，y_b，z_b)判断其所在波速层编号为k，由岩爆或者施工爆破事件空间坐标、分层模型边界面解析式参数和传感器坐标，按公式(11b)生成岩层k(波速层k)目标函数并采用遗传算法计算层速度组合优化解，将其用于模块Ⅲ微震定位.

(3) 模块Ⅲ：微震定位.微震系统采集现场微震事件波形信号，拾取波形起跳时刻.由分层模型边界面解析式参数、传感器坐标和模块Ⅱ中求解得到的层速度组合值，按公式(11b)依次生成各岩层(波速层)目标函数并采用遗传算法在可行域Ω_η中搜索得到相应微震震源参数优化解(x_η⁰，y_η⁰，z_η⁰)及相应目标函数值f_η.在满足解存在条件的各波速层中选取目标函数值最小者对应的微震震源优化解作为最终微震震源参数计算值.

实际工程应用中，微震监测初期阶段，根据地质资料对区域岩体进行岩层划分，并在建议的岩体波速值附近，利用模块Ⅱ进行层速度组合最优化求解.随着对区域岩体的施工开挖和支护开展，岩体质量不断改变、应力不断调整，实际层速度组合也随之改变.为了保证微震源定位精度，在监测期间，应定期对波速层重新划分并进行现场人工爆破试验对速度组合值进行更新计算.

白鹤滩水电站位于金沙江下游，左岸边坡为顺向坡，倾向河流上游南东方向，倾角较缓；岩体内部发育主要断层F14、F17及多个小断层、层内层间错动带.左岸拱肩槽建基面开挖高程为600~834 m，规模巨大.岩体开挖卸荷必然引起断层及层间、层内错动带高应力区裂隙萌生、发育及扩展，内部微破裂往往是岩质边坡外观变形和失稳破坏的前兆.

如图 4所示，为了实现对坝基边坡及其坡内廊道开挖卸荷过程围岩损伤的实时监测，左岸边坡安装了加拿大ESG微震监测系统，系统由HNAS主机处理系统、Paladin(台站)采集系统、单轴加速度传感器和三维可视化软件四部分组成.其中，HNAS通过网络控制各个台站(Paladin)，设置信号增益、采样频率和触发阈值等参数，从而实现全天候连续采集；Paladin为广泛使用的24位数据模数转换系统；单轴加速度传感器灵敏度为0.7Volts、频率响应范围为50 Hz~5 kHz.白鹤滩左岸边坡安装了3套Paladin采集系统，共18通道单轴加速度传感器阵列分别布置在610 m高程的1#灌排洞及施工连接洞、660 m高程的2#灌排洞和750 m高程的4#灌排洞及抗力体排水洞内.

图 4

Fig. 4

图 4 左岸边坡传感器阵列布置图 Fig. 4 Layout of sensor array on slope of left bank

在实际工程中，按照地层岩石物性将区域岩体分成若干个厚度不等的岩层，各层之间的波速值各不相同，这种波速的分层同地层的地质年代、岩性上的分层基本上是一致的.但地质分层比速度分层细，在通常情况下地质年代不相同但岩性相同的一些地层可以成为一个速度层，这在地震地层学中被称为层序“穿时现象”.

如图 5所示，白鹤滩左岸边坡为顺层岩质边坡，各岩层产状近似为N45°E，SE∠20°.开挖过后的坝基拱肩槽岩体岩质坚硬，呈微新状态，无卸荷或弱卸荷.根据坝基勘探资料，其岩体波速分类列于表 1中.如图 6所示，根据传感器网络布置范围和边坡坡面形状，确定模型尺寸和边界面；结合岩性和岩石类别，将模型区域岩体进行波速分层处理：共分为11层，4种类别，各波速交替出现，相应厚度由地质勘探图纸资料换算得到.

图 5

Fig. 5

图 5 左岸边坡拱轴线地质剖面图 Fig. 5 Geologic cross section along the arch axis at slope of left bank

表 1(Table 1)

表 1 岩体波速分类 Table 1 Classification of rock mass velocities 亚类 岩性 声波波速 Vps(m·s-1) 地震波波速 Vpd(m·s-1) 波速编号 Ⅱ1 斜斑玄武岩、隐晶质玄武岩、杏仁状玄武岩 >4700 >4100 V1 Ⅱ2 柱状节理玄武岩 >5100 >4700 V2 Ⅲ1 斜斑玄武岩、隐晶质玄武岩、杏仁状玄武岩、角砾熔岩 4200~4700 >3500 V3 Ⅲ1 柱状节理玄武岩 4700~5100 4500~4700 V4

表 1 岩体波速分类 Table 1 Classification of rock mass velocities

图 6

Fig. 6

图 6 分层速度定位几何模型 Fig. 6 Geometric model of SV

现场收集到610 m高程灌排廊道的4个开挖爆破数据，采用模块Ⅱ在波速建议值附近用遗传算法搜寻层速度组合最优解，计算结果列于表 2中.取由4个爆破事件计算得到的层速度优化值的近似平均值V₁=4500 m·s^-1，V₂=4900 m·s^-1，V₃=3750 m·s^-1，V₄=5000 m·s^-1作为分层速度定位模型的层速度组合值并用于模块Ⅲ进行微震定位.

表 2(Table 2)

表 2 层速度反演结果 Table 2 Inversion results of layer velocity 时间 爆破孔底坐标(m) 层速度优化值(m·s-1) N E U V1 V2 V3 V4 2015-04-18 15∶25 3012794.600 589164.126 613.343 4563.16 4857.44 3641.86 5032.97 2015-04-19 03∶02 3012750.942 589220.404 607.045 4533.67 4708.46 3863.04 4957.35 2015-04-19 04∶23 3012793.056 589166.741 612.979 4512.95 4881.00 3616.94 5050.29 2015-04-19 13∶01 3012848.356 589028.321 617.949 4450.42 5053.86 3938.78 4942.95

表 2 层速度反演结果 Table 2 Inversion results of layer velocity

本文利用ESG微震监测系统对左岸边坡微震事件进行实时采集，拾取各有效通道信号起跳时刻，然后采用Matlab语言编程来分别实现单一速度已知模型和分层速度定位模型定位.

2 015年4月1—18日期间，白鹤滩坝址区边坡坡外部分未进行开挖，工程施工集中在750 m高程置换洞和610 m高程灌排廊道的开挖.边坡层内错动带LS331和LS337受到坡内廊道开挖卸荷扰动，内部微破裂沿其分布带附近萌生和发展，是主要的岩体损伤区.单一速度已知模型采用4个波速区波速的近似均值4600 m·s^-1，分层速度模型采用前文所述模型.首先，利用单一速度已知模型对2015年4月1—18日期间微震监测系统采集的微震事件进行定位，圈定层内错动带LS331和LS337附近区域的微震事件；然后，利用分层速度模型对上述圈定事件进行重新定位.图 7给出了2015年4月6日18∶59∶55.234″事件的波形图，其中有12个通道监测到了有效的微震波形(DA_001等为通道编号).

图 7

Fig. 7

图 7 2015年4月6日18∶59∶55.234″事件波形 Fig. 7 Waveforms of an event at 18∶59∶55.234″ on April 6，2015

单一速度模型和分层速度模型的定位结果以小球的形式分别画在图 8a和8b中，其中小球的颜色表示微震事件发生所在的时间段.从图 8中定位结果图对比来看，单一速度已知模型定位的微震事件分散在层内错动带LS331和LS337附近，无明显的分布规律；分层速度模型定位的微震事件更接近层内错动带LS331和LS337，并呈条带状分布.通过计算某个微震事件到两个层内错动带的距离，然后选择较小的距离值作为评价该微震事件偏离层内错动带程度大小的指标.该指标值越小，表明微震事件的空间位置越靠近层内错动带；反之，则表明越远离层内错动带.计算结果表明，单一速度模型定位的微震事件中，最大偏离程度指标为32.505 m，最小为1.364 m，平均为9.134 m；分层速度模型定位的微震事件中，最大偏离程度指标为13.923 m，最小为0.866 m，平均为3.895 m.由此可以看出，采用分层速度定位模型，最大、最小和平均偏离程度指标分别降低了57.17%、36.51%和57.35%，这说明白鹤滩分层速度定位模型的微震事件定位结果较单一速度模型更加合理可靠.需要指出的是，精确的速度层划分是提高该算法定位精度的关键.白鹤滩分层速度定位模型中，对层速度组合进行更加精细的划分能够进一步提高定位精度.

图 8

Fig. 8

图 8 不同定位模型下微震事件定位结果的空间分布特征 (a)单一速度定位模型；(b)分层速度定位模型. Fig. 8 Spatial distribution characteristics of MS events located using different models

(1) 本文针对层状区域岩体，在异向波速模型的基础上，对垂向上应力波按波速值大小进行分段区别，推导了分层速度应力波走时关系式，建立了分层定位目标函数，并基于此提出了简化的分层速度定位模型SV.该模型包括参数准备、层速度反演、微震定位三个模块.

(2) 在进行分层速度目标函数优化求解时，采用了具有全局最优、不依赖优化初值和强鲁棒性特点的遗传算法作为最优化工具，使目标函数优化求解过程稳定性和全局最优的要求得到了保证.

(3) 将该速度模型应用到白鹤滩水电站左岸边坡，通过定位效果对比，分层速度定位模型定位微震事件的最大、最小和平均偏离层内错动带程度指标较单一速度模型分别降低了57.17%、36.51%和57.35%.

(4) 本文提出的改进定位模型考虑了因各岩层应力波波速值大小不同而对走时计算的影响，同时能采用遗传算法对建立的目标函数实现稳定求解.相比单一速度模型，定位精度得到了提高.然而，受到岩体应力松弛和风化作用影响，同层岩体各部位波速存在差异，囿于层速度划分无法达到绝对精细，势必在一定程度上影响微震事件的定位精度.为了提高定位精度，后期需要考虑对层内速度进行分区，在寻找能够稳定实现、更加符合应力波实际走时计算的定位模型方面进行改进和完善.