技术领域
[0001] 本发明涉及一种包含多层速度结构的微震震源定位方法。
相关背景技术
[0002] 自然环境中的矿产资源对于人类社会具有非常重要的意义，在整个国民经济发展中扮演着重要的角色。随着资源需求的增加，矿产深度开采已成为一种趋势。而在矿产开采过程中，随着采空区的增加以及施工的影响，岩层地质结构会发生改变，这时矿区就可能发生微小的破裂，而随着破裂的增多就可能引发大的震动，危害人员的生命安全。如果能够对微震进行实时监测定位，通过对微震震源的监测分析，及早发现隐患，就可以预防灾害发生，将生命财产损失降到最低。
[0003] 现有的微震震源定位算法，诸如网格搜索算法、牛顿迭代法等，在监测区域较小时可以高效的确定微震事件的发生位置。而在实际情况中，待监测的矿区范围往往较为庞大，由于部分算法对初值的选取较为敏感，这就导致了定位速度慢、定位不准确的问题。另外，将待监测区域等效为介质均一，波速单一的环境时，使用单一波速模型进行震源定位往往可以得到较为准确的微震事件发生位置。然而，在实际环境中，矿区往往存在着复杂的工程地质环境，地震波在不同的地质构造中的传播速度往往有着较大差异，采用单一波速会极大的干扰微震震源定位结果，造成较大的误差。
具体实施方式
[0038] 以下结合附图，对本发明上述的和另外的技术特征和优点进行清楚、完整地描述，对本发明做进一步详细说明。
[0039] 如图1所示为本申请一种包含多层速度结构的微震震源定位方法，主要包含以下步骤：通过检波器采集微震信号数据。拾取微震信号初至时刻。在同一时间段内拾取到微震信号初至时刻的检波器的数量大于预设数量时根据待定位区域介质层结构构建不同速度介质交界面方程。通过粒子群算法初步获取震源定位结果。通过牛顿迭代法对初步获取到的震源定位结果进行优化得到最终震源位置结果。本申请的包含多层速度结构的微震震源定位方法将矿区不同介质层假设为波速不同的结构，构建介质层交界面平面方程，对检波器采集到的微震信号进行处理并尝试拾取微震初至时刻，将拾取到初至时刻的微震采用以粒子群算法为辅、牛顿迭代法为主的搜索算法进行微震震源定位，完成微震事件的震源定位。通过上述步骤，本申请采用多层速度模型更加符合真实的矿区地质环境，且通过粒子群与牛顿迭代法相结合的震源定位算法更加快速准确。如图2所示为本申请所采用的技术方案框图。其中，振动信号采集模块用于进行信号采集。初至时刻采集模块用于拾取微震信号初至时刻。速度层构件模块用于根据待定位区域介质层结构构建不同速度介质交界面方程。微震震源定位模块用于通过粒子群算法和牛顿迭代法完成震源定位。以下具体介绍上述步骤。
[0040] 对于步骤：通过检波器采集微震信号数据。
[0041] 在待预测区域提前排布检波器用以采集震动数据。排布检波器时，各个检波器分布位置具有明显的差异，同时避免将数个检波器排布在同一个球面，从而造成震源定位结果出现多个解或无解的情况。可以理解的是，有效检波器数量增加可以提高震源定位结果的准确性。在本申请中，有效检波器数量需大于4个。
[0042] 对于步骤：拾取微震信号初至时刻。
[0043] 微震波初至时刻的拾取是震源定位的数据基础。而在矿区采集到的震动信号数据往往不仅包含着微震信号，还有大量各种各样的噪声信号。如果不去除这些噪声信号，会对微震初至时刻的拾取产生较大的干扰，从而影响震源定位的准确性。
[0044] 因此，在上一步骤的通过若干检波器采集微震信号数据之后，包含多层速度结构的微震震源定位方法还包含：对采集到的微震信号数据进行预处理。对采集到的波形数据进行降噪，以降低噪声水平。
[0045] 在本申请中，对采集到的微震信号数据进行预处理的具体方法为：对采集到的微震信号数据进行低通滤波和陷波降噪处理。由于微震信号通常为200Hz以下的低频信号，本发明采用200Hz的低通滤波器对采集到的原始信号进行滤波。另外，工频及其谐波的干扰也是微震信号噪声的重要来源，因此在进行低通滤波后，采用陷波器针对50Hz的工频信号及其谐波进行滤波，即完成了降噪处理。
[0046] 完成上述的微震信号数据预处理之后，进行微震信号初至时刻的拾取。
[0047] 在本申请中，通过STA/LTA算法拾取微震信号初至时刻。STA/LTA算法是根据两个一长一短、首尾相连的时窗在微震波形上不断滑动向后，计算短时窗内数据的均值与长时窗内数据均值的比值(ratio)。由于当微震事件到达时，短时窗的信号均值会发生迅速的变化。因此，设定一定的阈值，当比值大于该阈值时，即判定该时刻为微震事件的到时时刻。STA/LTA计算公式如公式(1)‑(3)所示，
[0048]
[0049]
[0050]
[0051] 其中lw为长时窗长度、sw为短时窗长度、ratio为长短时窗比值、X(i)为预处理后的微震信号数据。当比值超过预设的阈值时，即可判定为微震事件的到达时刻。
[0052] 但由于噪声水平、微震发生位置及频率等的影响，会导致STA/LTA算法在一些情况下的拾取结果出现较大的误差。因此，为了使得比值对信号中的振幅和频率变化更加敏感，在本申请中，对预处理后的微震信号数据采用如公式(4)所示的特征函数进行处理，[0053] F(i)＝X(i)2‑X(i‑1)X(i+1) (4)
[0054] 这样，当同一时间范围内拾取到初至时刻的检波器数量大于预设数量时，进行下一步的微震事件震源定位。
[0055] 对于步骤：在同一时间段内拾取到微震信号初至时刻的检波器的数量大于预设数量时根据待定位区域介质层结构构建不同速度介质交界面方程。
[0056] 在本申请中，将预设数量设定为4。
[0057] 具体地，根据待定位区域介质层结构构建不同速度介质交界面方程的具体方法为：采用空间平面方程构建介质层交界面，空间平面方程的一般解析式如公式(5)所示，[0058] Ax+By+Cz+D(i,i+1)＝0 (5)
[0059] 其中，A,B,C是介质层i与介质层(i+1)的交界面法向量的三个分量，D(i,i+1)为常数项。将平面上的三个点带入方程，即可求得各个速度层交界面方程的解析式。在本申请中，设定介质层交界面互相平行，因此，各个平面的交界面方程的A,B,C值相等。
[0060] 地震波在不同岩体介质下的传播遵从费曼原理按照折射路线进行传播。为简化模型，本申请中，仅考虑垂直方向上波传播路径穿过各个介质层波速不同的情况，将地震波在不同介质层上的路径近似为震源位置垂直经过介质层时两介质分割平面的距离，通过计算不同介质层的空间距离及波速值，求得地震波从震源到检波器的总走时。假设震源(x0,y0,z0)处于j介质层，检波器i(xi,yi,zi)处于k介质层。V为不同介质层对应的不同的速度。如公式(6)‑(8)所示为震源到检波器i的走时计算公式，
[0061]
[0062] αi＝|Axi+Byi+Czi‑(Ax0+By0+Cz0)| (7)
[0063]
[0064] 其中，αi,βi分别为空间距离的调节系数。
[0065] 对于步骤：通过粒子群算法初步获取震源定位结果以及通过牛顿迭代法对初步获取到的震源定位结果进行优化得到最终震源位置结果。
[0066] 具体地，根据前述的微震信号初至时刻以及分层速度微震波走时关系，可构建如公式(9)的目标函数，
[0067]
[0068] 其中，t0为微震信号的起始时刻，ti为第i个检波器拾取到的微震信号到时时刻，Δti为微震信号从发生位置传播至第i个检波器的理论时间。
[0069] 本申请以目标函数(9)最小为求解目标，本质是非线性方程的求解问题。在非线性方程求解中，粒子群算法由于其收敛速度快得到广泛应用，但该算法在局部搜索时精度会有所下降，容易陷入局部最优解。而牛顿迭代法虽然搜索精度高，但其依赖于初始值的选取，不恰当的初始值容易造成算法不收敛的问题。鉴于以上两种算法的优劣，本发明首先采用粒子群算法得到初步的震源位置，以该位置作为牛顿迭代法的初值进一步进行搜索，得到最终的震源定位结果。将粒子群算法与牛顿迭代法相结合进行震源定位，既避免了粒子群算法容易陷入局部最优而导致结果不够精确的问题，又避免了牛顿迭代法对初始值选取敏感的问题。
[0070] 如图3所示为粒子群算法的流程图。粒子群算法将树林比作一搜索空间，鸟群中的个体对应单个的粒子。给予每个个体一个随机的初始位置及初始速度，个体从初始位置以一定的速度在树林中飞行，在鸟群中选举出最优的个体，其他个体会向最优个体‘学习’，对自身的速度及位置进行动态更新，以趋近于最优的个体。更新方法如公式(10)、(11)所示。每次更新后再选举出新的优秀个体，重复以上过程进行迭代，直到所有个体均趋于最优解附近，即完成搜索。最优解即为最终结果，
[0071] vij(n+1)＝ωivij(n)+c1r1[pij(n)‑xij(n)]+c2r2[pbestj(n)‑xij(n)] (10)[0072] xij(n+1)＝xij(n)+vij(n+1) (11)
[0073] 其中，c1,c2为学习因子，通常取c1＝c2＝2。r1,r2为(0,1)区间的随机数。i为第i个粒子，j为解空间维度。Pij为第i个粒子在第j维度上的个体最优位置，Pbestj为种群在第j维度的最优位置，ωi为惯性权重。本发明中惯性权重根据粒子群算法的迭代次数进行线性更新。惯性权重的更新公式如公式(12)所示。
[0074]
[0075] 其中，i为当前迭代次数，max_step为最大迭代次数，ω范围取[0.2,1]。可以看出，粒子群的每次更新都包含了粒子个体的更新与种群的全局更新。
[0076] 由于粒子群算法在局部搜索时精度会有所下降，因此采用牛顿迭代法进一步精确震源位置。
[0077] 将粒子群算法求得的解作为初始值带入求解方程(9)得到残差值，沿着目标函数值减小的方向进行迭代，当残差小于某一阈值时结束迭代，即可得到最优解。迭代公式如公式(13)所示，
[0078] fn+1＝fn‑Hess(fn)‑1Grad(fn) (13)
[0079] 其中，Hess为fn上的黑塞矩阵，Grad为fn上的梯度，两者表达式如公式(14)、(15)所示，
[0080]
[0081]
[0082] 其中，m为未知数个数,f为目标函数表达式，在本申请中如公式(9)所示。xi为第i个未知数。
[0083] 采用上述粒子群算法得到的震源定位结果作为牛顿迭代法的初始值，带入牛顿迭代法中，利用牛顿迭代法进行更为精确的震源定位，即可获取最终震源位置结果。采用该算法可以快速准确的获取震源定位结果。
[0084] 以上显示和描述了本发明的基本原理、主要特征和优点。本行业的技术人员应该了解，上述实施例不以任何形式限制本发明，凡采用等同替换或等效变换的方式所获得的技术方案，均落在本发明的保护范围内。
