CN108802814B - 一种隧道围岩微震波速的获取方法 - Google Patents

Abstract

本发明公开了一种隧道围岩微震波速的获取方法，该方法基于启发式类算法中的引力搜索法GSA，在隧道掌子面后方布置至少4个微震监测传感器，建立隧道空间坐标系并准确测定传感器空间坐标，利用人工爆破事件触发传感器震动响应，拾取各个传感器上波形信号的观测到时；并设定观测到时与计算到时的累积绝对差值为微震波速反演的目标函数，最终通过计算目标函数满足终止准则来搜索微震波波速，本发明能够利用隧道施工的爆破事件实时准确获取隧道围岩的微震波速，保证了微震定位精度；该方法不需要专门的爆破活动确定微震波速，经济实用，操作简单，适用于各种交通、水利水电等隧道工程微震监测。

Description

一种隧道围岩微震波速的获取方法

技术领域

本发明涉及微震监测领域，尤其是一种基于启发式算法GSA的隧道围岩微震波速获取方法，适用于各种交通、水利水电等隧道工程微震监测。

背景技术

微震是指岩体应力超过自身强度时，岩体内部产生裂纹错动、开裂和破坏，累积的能量以弹性波形式释放并产生震动。微震震源定位是微震监测技术的核心，它是利用微震传感器记录的微震波形信号、到时数据和微震波波速反演微震事件的空间坐标和发震时刻。准确的微震波波速对于提高震源定位精度和稳定性至关重要，同时也是较难精确确定的震源定位输入参数。

目前，微震震源定位过程中，主要还是通过室内试验和现场勘探来获取微震波波速。室内试验是指在室内通过对岩块试样进行声波试验而获取岩块的波速；现场勘探主要指在微震监测区域内利用现场试验反演微震波波速。由于微震监测区域地质构造和介质特性的复杂多变性，现场勘探只能反映微震监测区域中某些区域的微震波波速，并不能完全代表监测区域内所有路径下的微震波波速；而室内声波试验中所用的岩块并不能代表现场岩体的性质，将所测试的波速作为微震波速用于微震监测必然造成较大的定位误差。

现有技术中，中国专利CN103697999A公开了一种高应力硬岩TBM施工隧道微震波速实时获取方法，在高应力硬岩TBM施工隧道掌子面后方岩体内布置至少4个微震传感器，采集掌子面附近岩体破裂产生的震动信号，将掌子面附近发生的岩体破坏事件作为微震源，测量微震源及微震传感器坐标，提取微震传感器监测到的微震源震动信号到时，反演获取微震波速。该方法反演的微震波速基于天然微震事件，由于岩体破裂位置未知，可能造成反演的微震波速存在较大误差。

现有技术中，中国专利CN104406681B公开的一种实时确定微震波波速的测试方法，利用微震监测系统自带的微震传感器和已知震源的人工爆破实时确定微震波波速。在微震监测系统已有的空间坐标系内准确测定各微震传感器和人工爆破事件的空间坐标；拾取微震监测系统中各微震传感器监测到的微震波到时，将第一个被触发的微震传感器到时作为比较对象，分别采用其它微震传感器和第一个被触发微震传感器的到时差以及它们相对于人工爆破震源的距离差实时反演微震波波速。随着采掘工程的进行。但是，该方法仅能获取单一的微震波速，不能获取微震波速的向量或不同传播路径上的波速，可能出现应用该微震波速进行定位时的较大误差。

现有技术中，中国专利CN104502964A公开了一种基于空间几何关系的获得微震波速的方法，根据震源点和基站的空间几何关系，构建微震等效波速V的不等式，构建形如开口向上的二次函数的变化规律的时差适应函数，提出一个单调单向搜索算法，求出最佳波速。但是，该方法仅能获取单一的微震波速，不能获取微震波速的向量或不同传播路径上的波速，可能出现应用该微震波速进行定位时的较大误差。

目前的现有技术还存在获取难度高、工作量较大、误差较大和效率不高的缺陷。因此，研究一种能够综合利用现场已知震源的人工爆破事件和已有的微震监测系统实时测试微震波波速的方法十分有意义，能够很大程度上保证微震震源定位精度。

发明内容

本发明为解决上述技术问题是提供一种能够综合利用现场已知震源的人工爆破事件和已有的微震监测系统实时测试微震波波速的方法。

本发明所采用的技术方案是，一种隧道围岩微震波速的获取方法，其特征在于，包括如下步骤：

a.在隧道施工掌子面后方布置n个微震检测传感器，n≥4；

b.采用人工爆破事件为震源位置，收集微震检测传感器接收范围内的岩体破裂的波形信号，定义所有传感器的观测到时与计算到时的累积绝对差值为震源波速反演的目标函数；

c.假定微震定位空间有多个粒子，并假定各个粒子的质量与多维波速向量；

d.以各粒子的多维波速向量计算目标函数，记录得到的最小目标函数值及对应最小目标函数值的多维波速向量；

e.判断最小目标函数值是否小于规定量值ε，若最小目标函数值小于规定量值ε，则对应最小目标函数值的多维波速向量为波速，若最小目标函数值不小于规定量值ε，则更新各粒子的多维波速向量，重新计算目标函数值，至最小目标函数值小于规定量值ε为止。

所述步骤e中各粒子的多维波速向量的更新采用引力搜索算法GSA(Gravitational SearchAlgorithm)进行。

引力搜索算法GSA(Gravitational SearchAlgorithm)是基于万有引力定律和牛顿第二定律的种群优化算法，万有引力是自然界4种基础力之一。在自然界中，万有引力的作用无处不在，使得任意一个粒子都会与其它的粒子相互吸引而不断的靠近，即较大质量的粒子能吸引较小质量的粒子，较大质量粒子可代表较准确的微震波速模型，从而最终获取到最准确的微震波速。因此，利用基于自然界物理法则的引力搜索算法能实现对微震源的准确定位与获取。

进一步，所述的规定量值ε，其范围为1e-4到1e-5。

作为优选，如上所述的一种隧道围岩微震波速的获取方法，包括如下步骤：

a在隧道施工掌子面后方布置至少4个微震监测传感器；

b采用人工爆破事件，以人工爆破事件位置为震源，建立隧道空间坐标系，定义所有传感器的观测到时与计算到时的累积绝对差值为震源波速反演的目标函数，所述目标函数的计算公式如下：

式中，fit为到时的累积绝对差值，n为传感器数量；t_i为第i个传感器的观测到时，上标p,s为P波或S波，t₀为震源的初始发震时刻，R_i/V为计算走时，R_i代表人工爆破事件位置(x₀，y₀，z₀)与第i个传感器位置(x_i，y_i，z_i)的距离；v_i代表震源与第i个传感器传播路径上的微震波速，对于所有的传感器则有微震波速向量V；

c假定隧道空间有N个粒子，每个粒子代表一种可能的微震波速模型，在初始时刻，每个粒子有质量M_i和n维波速向量V_i；

M＝(M₁,...,M_i,...,M_N)，(i＝1,2,...,N) (2)

式2代表有N个具有质量的粒子，式3代表每个粒子具有n维的波速。其中，表示第i个粒子在第d维的波速，并有上下限值，即

d将式(2)带入式(1)中计算各粒子的目标函数值fit，并记录历史循环的最小目标函数值Fbest及其对应粒子波速向量Vbest；

e判断当前最小目标函数值Fbest是否满足终止准则，所述终止准则为最小目标函数值Fbest是否小于规定量值ε，即是否Fbest<ε。如果Fbest<ε则输出最优微震波波速Vbest；如果Fbest＞ε则执行下一步骤，所述ε的范围为1e-4到1e-5，ε取值越小，获取的微震波波速精度就越高，但需要运算时间也会越久，因此ε的具体取值与工程实际情况相关；

f计算粒子间相互作用的引力，在第k次迭代，定义为在d维度上粒子i受到粒子j作用的引力，计算公式如下：

式中，M_aj(k)为主动粒子j的惯性质量，M_pi(k)为被动粒子i的惯性质量，为粒子j在第d维度上的波速，为粒子i在第d维度上的波速，ε为规定量值，R_ij(k)为粒子i和粒子j的欧式距离，G(k)为引力系数函数，G(k)和R_ij(k)满足如下公式：

式中，G₀与α为确定值，k为当前迭代次数，K为迭代总次数；随着当前迭代次数逐渐增加到总迭代次数，引力系数逐渐衰减；

R_ij(k)＝||X_i(k),X_j(k)||₂ (6)

第d维度上第i个粒子受到其他所有粒子引力作用的总和为：

式中，rand_j为[0，1]之间的随机数，为引力的总和添加随机；

根据以下公式更新粒子的惯性质量M_i：

M_ai＝M_pi＝M_ii＝M_i,(i＝1,2,...,N) (8)

式中，M_ai表示主动粒子i的惯性质量，M_pi表示被动粒子i的惯性质量，M_ii表示粒子i的惯性质量，fit_i(k)表示粒子i在第k次迭代的目标函数值大小；

对于求解最小目标函数值问题，best(k)为N个粒子在第k次迭代的目标函数的最小值，worst(k)为N个粒子在第k次迭代的目标函数的最大值，worst(t)和best(t)的定义如下：

g计算每个粒子的加速度。根据牛顿第二定理，第d维上粒子i的加速度为：

式中，F_i ^d(k)为当前迭代次数粒子i受到的引力总和，M_i(k)为粒子i的惯性质量；

h在每一次迭代中，每个粒子都会根据以下公式更新粒子速度和微震波速

式中，为第k+1次迭代更新后的粒子速度，为第k次迭代的粒子速度；

i将步骤h中各粒子更新的微震波速带入式(1)计算目标函数值fit，并判断是否Fbest<ε。若满足则退出循环，输出全体粒子目标函数的最小值Fbest和对应的最优微震波波速Vbest；若不满足，则继续执行迭代循环。

进一步，G₀可取100，α可取20。

本发明的有益效果是：

(1)采用施工活动爆破事件作为微震震源反演微震波波速，人工爆破事件的空间坐标能够采用全站仪精确测量，从而避免了因震源空间坐标误差导致的微震波波速误差，提高了微震波波速的精度。

(2)本方法不需要进行专门爆破获取微震波速，经济实用，操作简单，参数设置较少，收敛速度快，且能够很好的和各种优化问题相结合，且通用性比较强，对问题中不确定的信息具有一定的适应能力。

(3)能够实时获取微震波波速，对微震震源定位中采用的微震波波速进行实时调整，从而解决了由于隧道施工导致的监测区域中微震波波速可能发生变化的问题，保证了微震波波速的实时准确性，提高了微震震源定位精度。

附图说明

图1是万有引力作用示意图。

图2是本发明隧道围岩微震波速获取方法流程图。

具体实施方式

下面结合附图与实施例对本发明做进一步说明。

以某高地应力高速公路双洞隧道安装微震监测系统为例，依照图2流程获取微震波速。在该高地应力高速公路双洞隧道安装微震监测系统中，该监测系统含有9个微震传感器，分别布置在左右双洞掌子面后方，其中滞后洞布置3个，超前洞布置6个，分2个断面布置。由此开展针对隧道围岩内部破裂发生的实时监测，捕捉破裂触发的微震波形信号，拾取微震波形信号在每个传感器上的观测到时。在所建立的隧道空间坐标系内采用全站仪精确测得各传感器空间坐标(见表1)。2017年3月9日15时37分30秒在微震监测范围内进行了人工爆破试验，采用全站仪准确测量该人工爆破事件空间坐标为M(1008，770，1010)。同时，提取对各微震传感器所接受的微震波形到时(见表1)。

表1各微震传感器坐标

定义所有传感器的观测到时与计算到时的累积绝对差值为震源波速反演的目标函数，其计算公式如下：

在上述基础上，基于启发式算法GSA获取隧道围岩微震波速，其参数选取如下：K＝1000，G₀＝100，α＝20。最终根据上述发明，计算的最小目标函数Fbest＝5.44E-4，对应的粒子波速向量Vbest＝(5601，5566，4178)，即该实例反演获取的最优微震波波速v¹＝5601(对应传感器S1、S2和S3)，v²＝5566(对应传感器S4、S5和S6)，v³＝4178(对应传感器S7、S8和S9)。

以上实例仅用以说明本发明的技术方案而非限制，本领域的普通技术人员应当理解，本发明的技术方案进行修改或同等替换，而不脱离本发明方案的精神和范围，均应覆盖在本发明中。

Claims (4)

1.一种隧道围岩微震波速的获取方法，其特征在于，包括如下步骤：

a.在隧道施工掌子面后方布置n个微震检测传感器，n≥4；

b.采用人工爆破事件为震源位置，收集微震检测传感器接收范围内的岩体破裂的波形信号，定义所有传感器的观测到时与计算到时的累积绝对差值为震源波速反演的目标函数；

c.假定微震定位空间有多个粒子，并假定各个粒子的质量与多维波速向量；

d.以各粒子的多维波速向量计算目标函数，记录得到的最小目标函数值及对应最小目标函数值的多维波速向量；

e.判断最小目标函数值是否小于规定的规定量值ε，若最小目标函数值小于规定量值ε，则对应最小目标函数值的多维波速向量为最终获取的波速，若最小目标函数值不小于规定量值ε，则更新各粒子的多维波速向量，重新计算目标函数值，至最小目标函数值小于规定量值ε为止；

所述步骤e中各粒子的多维波速向量的更新采用引力搜索算法进行。

2.根据权利要求1所述的一种隧道围岩微震波速的获取方法，其特征在于：所述的规定量值ε，其范围为1e-4到1e-5。

3.根据权利要求1或2所述的隧道围岩微震波速的获取方法，其特征在于，包括以下步骤：

a在隧道施工掌子面后方布置n个微震监测传感器，n≥4；

b采用人工爆破事件，以人工爆破事件位置为震源，建立隧道空间坐标系，定义所有传感器的观测到时与计算到时的累积绝对差值为震源波速反演的目标函数，所述目标函数的计算公式如下：

式1中，fit为到时的累积绝对差值，n为传感器数量；t_i为第i个传感器的观测到时，上标p,s为P波或S波，t₀为震源的初始发震时刻，R_i/V为计算走时，R_i代表震源位置与第i个传感器位置的距离；v_i代表震源与第i个传感器传播路径上的微震波速，对于所有的传感器则有微震波速向量V；

c假定隧道空间有N个粒子，每个粒子代表一种可能的微震波速模型，在初始时刻，每个粒子有质量M_i和n维波速向量V_i；

M＝(M₁,...,M_i,...,M_N)，(i＝1,2,...,N) (2)

式2代表有N个具有质量的粒子，式3代表每个粒子具有n维的波速，其中，表示第i个粒子在第d维的波速，并有上下限值，即

d将式(3)带入式(1)中计算各粒子的目标函数值fit，并记录历史循环的最小目标函数值Fbest及其对应最优微震波波速Vbest；

e判断当前最小目标函数值Fbest是否满足终止准则，所述终止准则为最小目标函数值Fbest是否小于规定量值ε，如果Fbest<ε则输出最优微震波波速Vbest；如果Fbest＞ε则执行下一步骤；

f计算粒子间相互作用的引力，在第k次迭代，定义为在d维度上粒子i受到粒子j作用的引力，计算公式如下：

式中，M_aj(k)为主动粒子j的惯性质量，M_pi(k)为被动粒子i的惯性质量，为粒子j在第d维度上的波速，为粒子i在第d维度上的波速，ε为规定量值，R_ij(k)为粒子i和粒子j的欧式距离，G(k)为引力系数函数，G(k)和R_ij(k)满足如下公式：

式中，G₀与α为确定值，k为当前迭代次数，K为迭代总次数；

R_ij(k)＝||X_i(k),X_j(k)||₂ (6)

第d维度上第i个粒子受到其他所有粒子引力作用的总和为：

式中，rand_j为[0，1]之间的随机数；

根据以下公式更新粒子的惯性质量M_i：

M_ai＝M_pi＝M_ii＝M_i,(i＝1,2,...,N) (8)

式中，M_ai表示主动粒子i的惯性质量，M_pi表示被动粒子i的惯性质量，M_ii表示粒子i的惯性质量，fit_i(k)为粒子i在第k次迭代的目标函数值大小；

对于求解最小目标函数值问题，best(k)为N个粒子在第k此迭代的目标函数的最小值，worst(k)为N个粒子在第k此迭代的目标函数的最大值，worst(t)和best(t)的定义如下：

g计算每个粒子的加速度，根据牛顿第二定理，第d维上粒子i的加速度为：

式中，为当前迭代次数粒子i受到的引力总和，M_i(k)为粒子i的惯性质量；

h在每一次迭代中，每个粒子都会根据以下公式更新粒子速度和微震波速

式中，为第k+1次迭代更新后的粒子速度，为第k次迭代的粒子速度；

i将步骤h中各粒子更新的微震波速带入式(1)计算目标函数值fit，并判断是否Fbest<ε，若满足则退出循环，输出全体粒子目标函数的最小值Fbest和对应的最优微震波波速Vbest；若不满足，则继续执行迭代循环。

4.根据权利要求3所述的一种隧道围岩微震波速的获取方法，其特征在于：G₀为100，α为20。