无需测速的微震源反演数学模型研究
2023-09-03庞聪李查玮马武刚江勇廖成旺
庞聪 李查玮 马武刚 江勇 廖成旺
摘要：微震源反演数学模型、反演算法等对象是微震源定位误差的重要来源。针对P波速度易衰减、首波到时数据拾取误差较大的情况，基于到时差模型（ATD）和P波到时方程组，建立到时差比值法反演数学模型（ATDRM），从目标函数中消除了P波速度和发震时刻等参数；同时，引入范数概念，得到改进的L1、L2、L3、L4范数模型，设计一系列反演模型对比实验和反演算法对比实验综合分析研究这些模型的优越性。工程数据显示：ATDRM模型的定位效果整体比经典ATD模型更好，ATDRM-L1模型震源误差可减少13.437 0 m；反演模型的定位误差受反演算法影响较大，单形替换法相比遗传算法（GA）、非线性最小二乘法（NLSM）、模拟退火算法（SA）等方法定位精度更高，稳定性更好；在应用单形替换法求解模型时，反演误差会随着范数参数N的增加（L1→L4）呈现先下降再上升的现象；ATDRM模型及其范式形式对提高微震源定位精度有一定价值。
关键词：反演模型；微震源定位；范数；单形替换法；随机三角；P波速度
doi：10.16256/j.issn.1001-8956.2023.01.006
在微震源定位中，有许多因素影响着定位精度，这些因素造成的定位影响统称为系统误差，包括监测点坐标测量误差、波速误差、算法误差等［1］，而反演模型将会直接影响到反演算法误差的大小，因此反演模型的正确选择对微震源定位具有重要意义。常见的微震源反演数学模型主要有基于走时差［2］的反演模型、基于到时差［3］的反演模型、基于到时差商［4］的反演模型、基于矢量走时差［5］的反演模型等。其中，走时差模型由单个监测点的理论走时和观测走时之差来获得；到时差模型基于微震观测网中2个有效监测点的理论观测到时差与观测到时差的差值；到时差商模型包括3个监测点数据，以两两组合成分式的形式组成模型目标函数；矢量走时差模型基于走时差模型，引入空间三分量概念，建立空间位移差目标函数，相比于最小路程差，避免了算法的多解性问题。
针对微震震源反演数学模型的研究已有许多显著成效，但较多集中在到时差模型与走时差模型，虽取得了一些实际应用价值，但是对于可消除波速参数的到时差商模型及其改进形式的反演效果讨论存在较多空白。例如，丁恩杰等［4］创新地应用TDOA定位原理确立微震定位数学模型，该方法减少了速度因素带来的误差影响，但是研究内容仅局限于若干监测点布置在单维直线上（二维平面）的情况，代表性不足，且微震源定位目标函数存在主观简化的现象，没有给出相应简化的理由；戴峰［6］等简单描述了走时差模型、到时差模型及到时差商模型的基本内涵，并在分层反演中应用到走时差模型与到时差模型，并未对到时差商模型进一步讨论；李绍红［3］等建立的3参数模型和4参数模型本质上属于到时差模型，而5参数模型为走时误差，对于反演模型的讨论局限于P波速度是否已知而定的模型未知参数数目，以及4参数到时差模型相对于5参数走时差模型，消除了发震时刻。上述研究中的模型基本仍然包含地震波波速或首波初至时刻等未知参数，在一定程度上会将人工测量与数据处理带来的参数误差引入到微震震源定位中，进而降低定位精度；在传播介质较为均匀、地震波能量衰减不明显、小范围内定位等条件下，讨论并寻求一种无需测速的反演模型对提高微震震源定位精度是有必要的。
针对微震震源定位反演模型研究不够全面的情况，本文中在经典到时差模型基础上进行数学变换，通过差值法和比值法，将微震源定位精度的误差来源（P波速度和发震时刻）予以消除。各范数的引入，在一定程度上减小了算法本身带来的误差。然后加入到时差反演模型，对比研究不同的反演模型对单一算法的敏感性和稳健性，从而测试本文中提出的新模型的定位有效性；借助单形替换法、遗传算法、最小二乘法、模拟退火算法等，对比研究不同算法本身的算法误差带给反演模型的影响。
1 微震源反演到时差模型与波速分析
在微震源反演过程中，反演模型的选择对定位结果具有较大影响。TDOA（Time Difference of Arrival，到达时间差）原本属于一种室内无线定位技术，是一种利用时间差进行定位的方法，即到时差反演模型［3］。它不同于走时差模型［2］，其目标函數仅包含检波器空间坐标、P波速度、首波到时、发震时刻及震源等重要参数，这是由于到时差模型利用观测到时差与理论到时差相减而从数学公式上消去了发震时刻参数可以有效减小单个到时数据带来的误差，即时间检测不同步带入的误差。
根据微震源反演模型目标函数未知量的数量，可以将到时差模型的目标函数分为3参数型函数、4参数型函数等［3］，其中3参数型函数将P波速度作为已知量，其未知量是微震源的三分量坐标，4参数型函数的未知量除震源空间坐标外，还包括P波速度参数。3参数型模型虽然不包含波速参数，但是需要预先测速，而测速行为本身就会引入一系列误差，增加不必要的经济成本；4参数模型包含波速参数，虽然不会引入未知的波速误差，但是会增加算法寻优工作量，未知量的增多甚至可能导致算法陷入局部最优点或者算法收敛失败。
因此，预先测定的波速存在有3个方面的误差来源：测速设备误差、岩体误差、时空误差。测速设备误差主要和设备的稳定性、精确度、灵敏度等指标有关；岩体误差的产生是由于岩体的不稳定性、非均质特点等决定的，经过涉笔测得的数值也只能代表测速区域的波速平均值，这种单一的、静态的速度模型无法准确描述区域内P波传播的实际情况；时空误差主要变现为测定波速与实际波速在时间和空间上的不一致性，实际震源位置与测速区域不一定处在同一区域，且波速的传播路径随着岩体应力不断变化也不是固定的。
2 微震源反演数学模型的改进
1） 步骤1
记li为真实微震源（x0，y0，z0）与监测点（xi，yi，zi）的空间距离，并建立如下到时差方程组［3］
上式中，Δti、Δtj、Δto分别为P波从微震源传播至第i、j、o个微震监测点的走时数据；Δti=ti-t0，Δtj=tj-t0，Δto=t0-t0且i≠j≠o，其中t0为发震时刻，ti、tj、to分别为第i、j、o个检波器记录到的P波初至时间。
2） 步骤2
对方程组（1）进行数学形式转换，利用比值法，消去了P波速度参数v和发震时刻参数t0，即
将式（2）的左右两端相减，可得到微震源反演数学模型的到时差比值形式（Arrival-Time-Difference Ratio Model，简称ATDRM），即
上式中，i，j，o=1，2，…，n，i≠j≠o，且n（n≥3）为有效微震监测点的数量；由等式公式定义。按照等式定义，f的值理论上应为零；且区别于其他反演数学模型进行数值优化时的四四组合累加形式，当监测点数量不小于8时，三三组合的数值计算适应度函数将包含更少的累加项，累加次数为6×C3n。
3） 步骤3
为了深入了解到时差比值法模型对定位误差的影响，在优化计算中防止过拟合与提高解的可逆概率，引入范数概念［7］，其中L1范数和L2范数在最佳估计、正态分布统计等数据分析领域应用较广。L1范数表示为各个数据的绝对值之和，能有效实现稀疏，便于优化求解；L2范数被描述为各个数据求取平方和后再计算该累和值的平方根，增强矩阵的求逆能力；同时增加L3范数和L4范数作为对比模型，以对比观察范数与反演误差的关系。基于范数的微震震源定位目标函数（ATDRM-LN模型）的一般形式为
3 工程验证
3.1 监测数据描述
为了验证微震反演模型的有效性和创新性，应用吕进国［6］、王泉棟等［7］公开的柿竹园矿工程实验数据进行验证。该矿一共安装了30个微震监测单元，编号为N1～N30，爆破位置为（8 732.70， 6 570.60， 511.30）及起爆时刻为10：27。爆破成功后观测到8个P波到时数据，重新标记为T1～T8，监测点的空间坐标和首波到时数据如表1所示。
为了了解监测范围内各个三点组合的情况，分别计算各个监测点之间的空间距离以及三角周长，详细数据见表2。
由表2可知，空间距离较小（<55 m）的两点组合为（T1，T2）、（T1，T7）、（T3，T4）、（T6，T7），空间距离较大（>150 m）的两点组合为（T1，T5）、（T5，T6）、（T5，T7）、（T5，T8）、（T6，T8）；同时，就三点形成的三角形周长而言，（T1，T2，T7）、（T1，T6，T7）、（T2，T3，T4）等组合的周长明显较小，分别为164.983 5 m、198.301 2 m、199.316 3 m，它们之间的地震波传输速度大小在较小观测范围内可视作相同，更适合ATDRM模型及其变形的应用。
3.2 不同反演模型下的定位结果对比
将ATD模型的各范数形式分别记为ATD-L1、ATD-L2、ATD-L3、ATD-L4，ATDRM模型的各范数形式分别记为ATDRM-L1、ATDRM-L2、ATDRM-L3、ATDRM-L4。微震源反演算法统一采用单形替换法（又称单纯形法，通过对初始四面体单纯形进行扩张、收缩、镜像等操作，建立起新的四面体，而将极小点包括在内，从而不断逼近最优解），扩张因子为1.8，收缩因子为0.9，反演模型定位结果如表3所示。
根据表3可知：相对于微震源反演模型的标准形式，模型范数形式的定位误差有所减小，到时差L1模型定位误差最小，相对ATD下降3.924 3 m，但是在L2、L3及L4范数模型下突然出现误差反弹上升的情况，这说明微震源定位误差与范数并不是连续的正相关关系；而到时差比值法L1模型相对ATD模型下降13.437 0 m，其他范数形式误差总体保持在8.605 2 m；在实际应用中，微震源反演模型取标准模型的L1范数或L2范数形式即可。
3.3 不同模型求解方法的定位效果对比
经过国内外学者的多年研究，反演方法不断推陈出新，且不同的定位方法对反演结果的精确度具有差异性影响。李健等［8］将微地震定位方法分为迭代算法、非迭代算法、区域定位法及三分量定位法等，李自红等［9］将反演方法分为几何作图法、绝对定位法、相对定位法及非线性定位法等，其中绝对定位法和相对定位法都是线性方法。例如，遗传算法［10］、单形替换法、粒子群算法［11］、模拟退火算法等属于迭代算法或非线性法，INGLADA法、四四组合法［12］等属于非迭代算法，Geiger法［13］、JED法属于绝对定位法，主事件定位法和双差定位法［14］属于相对定位法；其中，单形替换法区别于其他迭代类算法，它是一种直接搜索方法，而模拟退火算法等是随机搜索方法，模拟退火算法与遗传算法虽同属于启发式算法，但是模拟退火算法的局部寻优能力较强，全局寻优能力较差，遗传算法与之相反，全局寻优能力比局部搜索能力较强。总之，针对微震震源反演数学模型的数值求解方法研究，对提高微震震源定位精度，剖析反演过程的可靠性，以及评价反演模型的优劣有重要意义［15-19］。
为了验证分析本文中提出的各类反演模型的有效性和稳健性，以及探索不同范数反演模型与不同类型定位方法的关联，利用遗传算法（Genetic algorithm，GA）、非线性最小二乘法（Nonlinear least square method，NLSM）、单形替换法（Simplex method，SM）、模拟退火算法（Simulated annealing，SA）等经典定位方法进行微震源定位解算。不同算法下的定位结果见表4和图1。
根据图1可知，算法定位结果表现最好的为单形替换法，其得到的定位误差总体最小；遗传算法与模拟退火算法的误差总体较大，且SA的结果不唯一，这是由于GA和SA算法为有约束局部寻优算法，算法结果除了与约束条件相关外，还易出现局部最优解，较难在全域内寻找最优值。这说明，在同一反演模型下，不同反演算法带来的算法误差有较大差距。
4 结论
基于到时差模型建立无需提前测速的反演数学模型及其范数形式，利用实际工程数据验证这些模型的优劣性以及对比在不同反演方法下的定位效果差异，并观察范数参数对反演模型的影响，实验结论为：
1） 在经典到时差模型（ATD）基础上，提出了5种到时差比值法模型：ATDRM标准型及其范数形式（L1、L2、L3、L4），从数学模型上消除对微震源定位精度影响较大的P波速度、发震时刻等难以精准确定数值的参数，克服了传统微震震源定位中波速难以确定的情况，简化现场实际微震震源定位的复杂流程，减少不必要的现场测量任务。
2） 利用单形替换法，求解包括ATD模型及ATDRM在内的10种微震源反演模型目标函数，对比分析本文中反演模型的优劣。结果表明，ATDRM的定位结果整体优于经典到时差模型，并在范数形式下具有一定的稳定性，定位精度较高。
3） 分别采用遗传算法、非线性最小二乘法、单形替换法、模拟退火算法等经典算法对反演模型进行定位解算，结果表明应用单形替换法的实验结果相对于标准模型的定位精度更高，求解稳健性更优，没有出现明显的定位结果偏差现象，这是由于单形替换法在数值计算过程中不需要求解逆矩阵和偏导数，通过及时更新最小误差近似解避免发散现象的发生，提高定位过程稳健性和计算效率，而其他方法易受到局部最小值的干扰；同时，观察到单形替换法外的其他解法具有随范数的增加而定位误差先下降后上升的现象，这是因为L1范数模型不易受极个别异常数据的影响，而L2范数模型相对更适合数据误差服从正态分布的定位计算问题，随着范数幂次的增加，L3与L4范数模型不仅运算量迅速增大，初始监测数据中的异常值或离群点也会加重对定位精度的影响，即受台网布设质量与检波器数量的影响愈加明显。
针对无需测速的反演数学模型的研究，丰富了微震震源定位方法理论体系，通过数学公式变换的形式消除波速参数，使反演仿真实验更加关注单一因素对定位精度的影响，对相关人员采用范数原理分析微震震源定位精度提供了一定参考依据，适用于浅层天然地震定位、矿山微震监测、页岩气资源勘探等重要领域，并要求检波器尽量布设在地震波传播介质相对均匀、岩体属性一致的小范围区域，从而减弱人工测量因素和波速变化对微震震源精准定位的干扰。
参考文献：
［1］ 董陇军， 李夕兵， 唐礼忠. 影响微震震源定位精度的主要因素分析［J］. 科技导报， 2013， 31（24）：26-32.
［2］ 呂进国， 姜耀东， 赵毅鑫，等. 基于稳健模拟退火-单纯形混合算法的微震定位研究［J］. 岩土力学， 2013，（8）：2 195-2 203.
［3］ 李绍红， 吴礼舟， 杨戒，等. 微震源定位的两步反演方法研究［J］. 岩石力学与工程学报， 2017， 36（7）：1 710-1 717.
［4］ 丁恩杰， 吕雅洁， 胡东平，等. 基于Spark平台的微震监测快速定位方法研究［J］. 煤炭科学技术， 2016， 44（7）：22-27.
［5］ 张唤兰， 朱光明， 王保利. 微地震震源定位中目标函数的研究与应用［J］. 煤田地质与勘探， 2015，（6）：105-108.
［6］ 戴峰， 郭亮， 徐奴文，等. 基于异向波速模型的微震定位改进［J］. 地球物理学报， 2016， 59（9）：3 291-3 301.
［7］ 李楠， 王恩元， 孙珍玉，等. 基于L1范数统计的单纯形微震震源定位方法［J］. 煤炭学报， 2014， 39（12）：2 431-2 438.
［8］ 李健， 高永涛， 谢玉玲，等. 基于无需测速的单纯形法微地震定位改进研究［J］. 岩石力学与工程学报， 2014， 33（7）：1 336-1 346.
［9］ 李自红， 陈慧. 地震定位方法研究综述［J］. 山西煤炭， 2018，（1）：35-40.
［10］王泉栋， 李国和， 吴卫江，等. 多种群遗传算法在微震震源定位中的应用［J］. 计算机测量与控制， 2015， 23（4）：1 285-1 288.
［11］陈炳瑞， 冯夏庭， 李庶林，等. 基于粒子群算法的岩体微震源分层定位方法［J］. 岩石力学与工程学报， 2009， 28（4）：740-749.
［12］朱权洁， 姜福兴， 缪华祥，等. 基于聚类分析的微震定位二次优化研究［J］. 采矿与安全工程学报， 2014， 31（2）：196-202.
［13］张唤兰，朱光明，王保利. 微地震震源定位中目标函数的研究与应用［J］.煤田地质与勘探，2015，43（6）：105-108.
［14］柳云龙， 田有， 冯晅，等. 微震技术与应用研究综述［J］. 地球物理学进展， 2013， 28（4）：1 801-1 808.
［15］龐聪，江勇，廖成旺. 地震观测仪器电子地图信息服务平台的设计与实现［J］.内陆地震，2018，32（4）：359-369.
［16］赵石柱，张勇，闫新义. Sg、Sn震相对新疆及相邻地区速报地震定位的影响分析［J］.内陆地震，2021，35（3）：237-244.
［17］庞聪，江勇，廖成旺，等. 基于机器学习的强震动监测环境抗干扰方法对比研究［J］.内陆地震，2020，34（2）：119-124.
［18］曹昌军，李艳永，任林，等. 地震速报工作中S震相对震中位置影响分析［J］.内陆地震，2022，36（2）：121-129.
［19］王桂丹.基于JOPENS/MSDP地震定位方法的定位结果分析［J］.内陆地震，2016，30（1）：74-80.
MATHEMATICAL MODELING OF MICRO-SEISMIC SOURCE
INVERSION IN MINES WITHOUT VELOCIMETRY
PANG Cong1，2， LI Cha-wei1，2， MA Wu-gang1，2， JIANG Yong1，2， LIAO Cheng-wang1，2
（1. Institute of Seismology，CEA，Wuhan 430071，Hubei，China；
2. Hubei Key Laboratory of Earthquake Early Warning，Wuhan 430071，Hubei，China）
Abstract： The micro-seismic source inversion mathematical model， inversion algorithm and other objects are important sources of micro-seismic source localization errors. To address the situation that the P-wave velocity is easily attenuated and the first wave arrival data pickup error is large， the inversion mathematical model of arrival time difference method （ATDRM） is established based on the arrival time difference model （ATD） and the set of P-wave arrival time equations， and the parameters of P-wave velocity and seismic moment are eliminated from the objective function; meanwhile， the concept of parametric number is introduced to obtain the improved parametric models of L1， L2， L3 and L4， and a series of inversion A series of model comparison experiments and inversion algorithm comparison experiments are designed to comprehensively analyze the superiority of these models. The engineering data show that the overall localization effect of the ATDRM model is better than that of the classical ATD model， and the source error of the ATDRM-L1 model can be reduced by 13.437 0 m. The localization error of the inversion model is influenced by the inversion algorithm， and the monomorphic replacement method has higher localization accuracy and better stability than the genetic algorithm （GA）， nonlinear least squares method （NLSM）， simulated annealing algorithm （SA） and other methods. The inversion error decreases and then increases with the increase of the parametric parameter N （L1→L4） when the monomorphic substitution method is applied to solve the model; the ATDRM model and its parametric form are valuable for improving the localization accuracy of micro-seismic sources.
Key words： Inversion model; Localization of micro-seismic source; Norm; Simplex substitution method; Random Triangle；P-wave velocity
