0 引言
获取震源准确位置对于动态灾害的早期预警具有重要作用[1-3]。通过震源传感器接收震波并根据震波的初至时刻反演出震源位置称为震源定位[4]。根据震源定位公式，可将震源定位问题转换成数学中未知量的求解问题。数学中未知量求解方法包括线性方法和非线性方法。常用的线性方法有梯度下降法、牛顿法、拟牛顿法等，非线性方法有模拟退火法，遗传算法(Genetic Algorithm，GA)、单纯形法等。其中，梯度下降法存在收敛速度慢等问题[5]。拟牛顿法需寻找替代矩阵，收敛速度快，但易陷入局部极小值点[6]。单纯形法从几何搜索的角度迭代，易收敛[7]，是目前的主流算法[8]。GA算法可全局搜索[9]，前期搜索收敛速度快，后期收敛速度慢[11-12]，且相较于单纯形法，具有天然的并行性，很容易应用到并行计算环境[13]中，顺应了目前物联网发展研究的大趋势[14]。
在震源置于远场时，定位精度急剧下降[15]。针对该问题，本文提出了基于改进GA算法的远场震源定位方法，综合利用拟牛顿法收敛速度快的特点和GA算法可全局搜索的特点实现震源定位，提高了震源定位的精度和速度。
1 震源定位原理
震源定位公式为
(1)
式中：(x,y,z)为震源坐标；(xi,yi,zi)是编号为i的检波器的坐标；vi是编号为i的检波器接收到的震波在介质中传播的平均波速；ti是编号为i的检波器接收到震动的时刻；ts为震源的起震时刻。
当检波器位于同一条直线上时，二维远场震源定位模型如图1所示。
图1 二维远场震源定位模型
Fig.1 Two-dimensional far-field seismic source positioning model
由图1可得
(2)
式中：v为波速；t12=t1-t2；t23=t2-t3；t34=t3-t4。
对式(2)进行处理可得到目标函数：
f(x,y)=
(3)
震源定位的问题转化为上述目标函数的优化求解问题，找到使目标函数最优的x和y值，就认为找到了震源位置。
2 改进GA算法原理及性能分析
2.1 GA算法
GA算法原理：首先确定搜索空间的范围，组成遗传空间的初始种群。为挑选出符合预期的个体，需要设定合适的适应度函数，即本文的目标函数。模拟种群交配原理，选择具有更高适应性的个体相互交配，以产生优良个体；同时随机修改待求解变量中的某些部分，变异出新个体，继续选择适应性更强的个体；循环上述过程，直到满足终止条件，得到最优解。
为了确保震源位置在所设定的范围内，刚开始会给出一个较大的空间范围。GA算法需要将该空间范围离散为很多的点逐个搜索，到接近真实值的时候收敛速度很慢。为了提高收敛速度，本文对GA算法进行了改进。
2.2 改进GA算法
改进GA算法将GA算法和拟牛顿法相结合，通过优势互补提高震源定位精度和速度。
改进GA算法流程如图2所示。首先使用GA算法进行初始寻优，当搜索范围缩小到一定程度时结束寻优过程；再将GA算法得到的最优解作为拟牛顿法寻优的初值进行精确寻优；当拟牛顿法寻优过程结束时，判断所得的最优解是否符合预期，如果符合则将其作为最终结果输出，否则使用GA算法进行一定次数迭代，将得到的最优解作为最终结果输出。由于GA算法不需要寻找精确解，所以，只要个体数量足够多、能够较好地覆盖最优解可能存在的位置即可。
图2 改进GA算法流程
Fig.2 Flow of improved GA algorithm
2.3 定位速度对比分析
为了验证改进GA算法的有效性，分别采用GA算法和改进GA算法进行仿真分析。设定如下参数：检波器坐标为(0，0)，(1 m，0)，(2 m，0)，(3 m，0)；介质波速v=230 m/s；在横轴[-50,50] m、纵轴[9,50] m范围内，分别以1 m为间隔设置震源，对共计(50+50+1)×(50-9+1)=4 242个震源进行定位。其他参数设置见表1。
表1 改进GA算法参数设置
Table 1 Parameters setting for improved GA algorithm
GA算法相关参数设置：最大遗传代数为200，变异概率为0.2，其他参数与改进GA算法相同。
采用2种算法定位得到的纵、横坐标的定位误差e1—e4如图3、图4所示。
(a) 横坐标误差e1
(b) 纵坐标误差e2
图3 遗传算法定位误差
Fig.3 Positioning error of GA algorithm
(a) 横坐标误差e3
(b) 纵坐标误差e4
图4 改进GA算法定位误差
Fig.4 Positioning error of improved GA algorithm
重复进行10次上述仿真，取10次仿真所用时间的平均值，GA算法用时1 885.155 s，改进GA 算法用时588.75 s，改进GA算法运行速度约为GA算法的3倍。这是由于拟牛顿法提高了GA算法的收敛速度，从而提高了震源定位的速度。
在图3、图4中的绝大部分区域，改进GA算法的误差曲线要比GA算法光滑，误差更小，定位速度更快。但是一旦出现误差，误差值会比较大。
2.4 改进GA算法定位误差分析
假设震源真实坐标为(30 m，25 m)，取搜索区域中的4组坐标(10 m，20 m)，(20 m，25 m)，(25 m，35 m)，(40 m，50 m)，每组重复进行15次反演定位，通过控制变量法将反演出的定位坐标和真实坐标作对比，结果如图5所示。
(a) 坐标(10 m，20 m)
(b) 坐标(20 m，25 m)
(c) 坐标(25 m，35 m)
(d) 坐标(40 m，50 m)
图5 定位坐标和真实坐标对比
Fig.5 Contrast of positioning coordinates and real coordinates
从图5可看出，x轴定位误差明显小于y轴，这是由于检波器的位置在同一条直线上，导致震源的到达时间差有变动。图5(c)和图5(d)中x轴定位误差波动较大，图5(b)中y轴定位误差波动较大。震源真实坐标为(30 m，25 m)，在算法搜索到越接近震源真实位置的时候，误差越大。这是由于在搜索过程中，目标函数越接近真实值，其变化越小，以致算法最终搜索到错误的位置。本文通过构建一种新的目标函数来解决该问题。
3 新目标函数的构建
3.1 新目标函数
在原目标函数中设置的是x坐标定位，检波器的单轴排列导致y轴的定位精度无法保证。为了改善y轴的定位精度，将检波器的位置改为双轴布置，在双轴上各自设置独立的子目标函数，保证双轴的误差最小，此时的二维远场震源定位模型如图6所示。
图6 双轴二维远场震源定位模型
Fig.6 Biaxial and two-dimensional far-field seismic source positioning model
以坐标分别为(x1,y1)，(x2,y2)的检波器为例，测得震波到达2个检波器的时间差T=Nt0，其中t0为AD转换周期，N为通过人工或机器读取的到达时间，在实际情况下，t0和N为整数，为了讨论方便，假设T>0。震源坐标(x,y)满足如式(4)所示的双曲线方程和式(5)所示的渐近线方程。
(4)
(5)
当震源置于远场的时候，渐近线和双曲线的误差非常小，用双曲线代替渐进线来定位，最终转换得到x关于v的函数：
(6)
当双坐标排列的检波器坐标为(x1,y1)，(x2,y2)，(x3,y3)，(x4,y4)，(x5,y5)，(x6,y6)，(x7,y7)，(x8,y8)时，将在x轴和y轴分别放置的4组检波器两两组合，得f1(x)=|x1-x2|，f2(x)=|x3-x4|。令目标函数f(x)=|f1(x)|+|f2(x)|，由于x是关于v的函数，要使得目标函数最小，只需要找到满足条件的最优解v，再由v反演相应的x，y值，即可对远场震源进行定位。
3.2 收敛速度对比分析
假设检波检波器的坐标为(970 m，0)，(980 m，0)，(990 m，0)，(1 000 m，0)和(0，970 m)，(0，980 m)，(0，990 m)，(0，1 000 m)，波速为3 000 m/s，震源坐标为(600 m，1 000 m)。通过计算可得各个检波检波器的到达时间差分别为t12≈1.170 4 ms，t23≈1.197 6 ms，t34≈1.224 6 ms，t56≈0.138 8 ms，t67≈0.083 3 ms，t78≈0.027 8 ms。
设原目标函数为A，新目标函数为B，其中B=|x1-x2|+|y1-y2|，进行20组实验，对比A和B的性能，结果见表2。
表2 新旧目标函数性能对比
Table 2 Performance comparison of the old and new target functions
新旧目标函数的收敛速度对比如图7所示。从图7可以看出，与A相比，B的收敛速度更加稳定，而且B在迭代20次之前就会收敛，收敛速度快于A。
图7 新旧目标函数的收敛代数对比
Fig.7 The comparison of convergence speed between old and new target functions
新旧目标函数的速度误差对比如图8所示。从图8可看出,B的速度一直很稳定，而A的速度不稳定，波动很大。同时B的速度和实际速度的误差一直稳定为385.3 m/s，而A的速度误差波动很大，刚开始为58.1 m/s，最大值为1 243.7 m/s，都没有收敛到真实速度，这是由于计算机的计时精度等因素引起，与目标函数和算法等无关。
图8 新旧目标函数的速度误差对比
Fig.8 The comparison of speed error between old and new target functions
3.3 新目标函数定位误差分析
假设x1在[-5，94] m范围内移动，y1在[7，106] m范围内移动，一共有100×100个点，通过计算得出真实波速，再通过搜索目标函数，得到满足条件的最优波速值，最后反演得到相应的震源坐标(x2，y2)。记震源横坐标的理论误差值ex=|x1-x2|，纵坐标的理论误差值ey=|y1-y2|，则采用新目标函数得到的震源定位误差如图9所示。
将图9(a),(b)分别与图3(a),(b)比较可知，采用原目标函数时，定位误差散落不均匀，同时y轴定位误差要小于x轴定位误差。而采用新目标函数得到的误差在绝大部分区域都很光滑，震源定位更接近真实值，且y轴定位精度得到明显改善。
(a)x轴
(b)y轴
图9 震源定位误差
Fig.9 Error of source positioning
4 结论
(1) 改进GA算法利用原有GA算法的全局搜索能力，同时结合拟牛顿算法的快速收敛特性，在一定误差范围内，可以提高远场震源定位的精度。仿真分析结果表明，在并行化的环境下，改进GA算法的平均定位速度为原GA算法的3倍。
(2) 针对检波器单轴排列导致y轴定位精度无法保证的问题，将检波器的位置改为双轴布置，在双轴上各自设置独立的子目标函数，保证双轴误差最小。新旧目标函数的性能对比分析表明，新目标函数的收敛速度更快、更稳定。在误差允许的范围内，用该目标函数实现震源定位，可有效提高y轴定位精度。
参考文献(References)：
[1] 张红才. 地震预警系统关键技术研究[D]. 哈尔滨:中国地震局工程力学研究所,2013.
[2] 程久龙,宋广东,刘统玉,等. 煤矿井下微震震源高精度定位研究[J]. 地球物理学报,2016,59(12):4513-4520.
CHENG Jiulong, SONG Guangdong, LIU Tongyu, et al. High precision location of micro-seismic sources in underground coal mine[J]. Chinese Journal of Geophysics,2016,59(12):4513-4520.
[3] 王子珺. 地震动破坏性与地震预警理论方法研究[D]. 北京:北京交通大学,2017.
[4] 李楠. 微震震源定位的关键因素作用机制及可靠性研究[D]. 徐州:中国矿业大学,2014.
[5] 郭跃东,宋旭东. 梯度下降法的分析和改进[J]. 科技展望,2016,26(15):115.
[6] 于杰. 改进的多步拟牛顿法及其收敛性[D]. 南京:南京航空航天大学,2011.
[7] 李健,高永涛,谢玉玲,等. 基于无需测速的单纯形法微地震定位改进研究[J]. 岩石力学与工程学报,2014,33(7):1336-1346.
LI Jian, GAO Yongtao, XIE Yuling, et al. Improvement of microseism locating based on simplex method without velocity measuring[J]. Chinese Journal of Rock Mechanics and Engineering, 2014, 33(7): 1336-1346.
[8] LI J,WU S C,GAO Y T,et al. An improved multidirectional velocity model for micro-seismic monitoring in rock engineering[J]. Journal of Central South University,2015,22(6):2348-2358.
[9] 马永杰,云文霞. 遗传算法研究进展[J]. 计算机应用研究,2012,29(4):1201-1206.
MA Yongjie, YUN Wenxia. Research progress of genetic algorithm[J]. Application Research of Computers, 2012, 29(4): 1201-1206.
[10] 赛琳伟. 基于遗传算法的团簇结构全局搜索[D].大连:大连理工大学,2012.
[11] 朱钰,韩昌佩. 一种种群自适应收敛的快速遗传算法[J]. 计算机科学,2012,39(10):214-217.
ZHU Yu, HAN Changpei. Improved genetic algorithm with adaptive convergence populations[J]. Computer Science, 2012, 39(10): 214-217.
[12] 汪民乐. 遗传算法的收敛性研究[J]. 计算技术与自动化,2015,34(1):58-62.
WANG Minle. Research on convergence of genetic algorithm[J]. Computing Technology and Automation, 2015, 34(1): 58-62.
[13] 余涛,刘泽燊. 基于Spark的并行遗传算法研究[J]. 计算机时代,2017(1):43-46.
YU Tao, LIU Zeshen. Research of massive parallel genetic algorithm based on spark[J]. Computer Era, 2017(1):43-46.
[14] 丁恩杰,刘亚峰,吕雅洁. 感知矿山物联网发展现状及展望[J]. 信息通信技术,2015,9(5):15-20.
DING Enjie, LIU Yafeng, LYU Yajie. Situation and prospect of IoT in mine[J]. Information and Communications Technologies, 2015,9(5):15-20.
[15] 徐逸鹤,徐涛,王敏玲,等.井中震源的远场波场特征研究[J].地球物理学报,2015,58(8):2912-2926.
XU Yihe,XU Tao,WANG Minling,et al. Far-field wavefield characteristics of downhole seismic sources[J]. Chinese Journal of Geophysics,2015,58(8): 2912-2926.