麻雀搜索算法在微地震定位中的应用研究

Application of Sparrow Search Algorithm in Microseismic Location

摘要: 为了更加准确地获取微地震震源所在的位置，基于麻雀搜索算法(SSA)，提出了利用纵波传播时差来定位的方法。该方法首先利用纵波时间差信息构建目标函数，并通过迭代搜索走时目标函数的全局最小值，从而求得定位误差相对最小的震源位置。此外，将该定位算法与传统Geiger算法进行定位仿真实验对比，进一步验证了SSA算法在微地震震源定位中的可行性及高精度特性。

Abstract: To obtain the location of microseismic sources more accurately, a method based on the Sparrow Search Algorithm (SSA) is proposed to locate by using the time difference of primary wave propa-gation. The method first constructs an objective function using the longitudinal wave time differ-ence information and iteratively searches for the global minimum of the travel time objective function to find the location of the source with relatively minimal positioning error. In addition, the localisation algorithm is compared with the traditional Geiger algorithm in a localisation simulation to further verify the feasibility and high accuracy of the SSA algorithm in microseismic source localisation.

1. 引言

微地震监测技术已经在盐矿地面塌陷监测 [1] 、地下储气库安全监测 [2] 、智能家居 [3] 等行业得到了广泛应用。在石油采集过程中，可以利用井下或地表布置的监测传感器接收非常规油气生产过程中的微小的破裂或错断信号，并通过一系列处理方法反推出信号发生的位置、延伸方向等信息，从而优化非常规压裂设计方案，以降低成本提高出产率 [4] 。震源定位是微震监测技术中最重要的优势功能，它能够描述不稳定区域的活动范围和时间特征。微震监测区域中所有传感器台站以及它们之间的空间形态构成了一个完整的微震监测台网，震源定位所需要的输入数据主要包括微震监测台网中各传感器空间坐标、各传感器的到达波到时和微震波在监测介质中的传播速度等数据。所以，如何更好更快速地反演出微地震的震源位置，是微地震监测的关键问题之一。

很多学者对微地震震源位置反演方法进行了大量有效的研究。1912年Geiger提出的经典方法 [5] ，将目标函数线性化，降低求解难度。赵仲和将Geiger的改进方法用于北京台网 [6] ，定位精度大幅度提高；2021年常亚琼 [7] 结合Chan算法和Geiger算法相结合的联合定位方法对微震震源进行定位，其结果显示定位精度更高。随着计算机的飞速发展与跨学科研究的不断深入，陈炳瑞 [8] 和Shishay [9] 等基于粒子群算法，在某种程度上提高了震源定位精度；王泉栋 [10] 发现遗传算法在定位方面具有较高的性能和精度。综上所诉，Geiger定位方法存在迭代初值选择不合适而造成求解效率低、定位准确性差甚至无法定位等问题，群智能优化算法具有计算方法简单，容易实现等优点被广泛应用于各行各业。本文以L2范数作为目标函数，通过麻雀搜索算法计算微地震震源发生的位置，最后通过数值实验对比证明本文算法的可行性。本文结构如下：1) 引言部分介绍了微地震定位的国内外研究现状，微地震定位存在的问题及本文提出的解决方法；2) 方法原理部分介绍了麻雀搜索算法的基本原理及在该算法上对L2范数的改进；3) 仿真实验部分通过在二维平面采集雷克子波震源的初至数据来验证本文方法的有效性与准确性；4) 结论部分对本文工作进行总结。

2. 方法原理

2.1. 麻雀搜索算法原理

麻雀搜索算法(Sparrow Search Algorithm, SSA) [11] 是由Xue等人，总结麻雀捕食行为中的优点得来，该算法具有寻优能力强，收敛速度快的优点。在该算法中麻雀扮演三种角色(捕食者，警戒者，加入者)。警戒者负责提供觅食的范围和方向，并且可以与加入者交换位置，二者保持动态统一。能量越低的加入者越容易成为警戒者。在捕食过程中，加入者总是可以找到优质的警戒者，并且驱使其去寻找食物资源，以获得最高的利益。当麻雀发现捕食者时会发出预警信号。并判断是否转移到其它安全位置进行捕食。当感觉到危险到来，离群体最远的麻雀会马上向群体靠拢，并且找到食物多的位置，重新捕食。种群内部的麻雀会为了捕食而随机走动。在算法中食物储备的多少取决于适应度值的优劣。其算法流程如图1所示。

Figure 1. Algorithm flowchart

图1. 算法流程图

麻雀中三者的位置变化公式如下。捕食者：

$X_{i,j}^{t+1}=\{\begin{array}{ll}{X}_{i,j}\cdot \exp \left(-\frac{i}{\alpha \cdot ite{r}_{\max }}\right),\hfill & \text{if}{R}_2<ST\hfill \\ X_{i,j}+Q\cdot L,\hfill & \text{if}{R}_2\ge ST\hfill \end{array}$ (1)

式中，t表示迭代次数， $\alpha$ 和Q是随机数， $j=1,2,3,\cdots ,d$ ， $ite{r}_{\max }$ 表示最大的迭代次数。 $X_{i,j}$ 表示第j维中第i个麻雀的位置信息。 $\alpha \in \left(0,1\right]$ ， $R_2\left(R_2\in \left[0,1\right]\right)$ 表示预警值 $ST\left(ST\in \left[0.5,1\right]\right)$ 表示安全值。Q是服从正态分布。L表示一个元素全为1的 $1\times d$ 矩阵。

当 $R_2<ST$ 时，捕食区域周围没有捕食者，警戒者可以在大范围内寻找食物。如果 $R_2\ge ST$ ，说明该区域内有捕食者，并给其他麻雀发出警报，前往其他安全区域觅食。

加入者：

$X_{i,j}^{t+1}=\{\begin{array}{ll}Q\cdot \exp \left(\frac{X_{worst}-X_{i,j}^t}{i^2}\right),\hfill & \text{if}i>n/2\hfill \\ X_P^{t+1}+\left|X_{i,j}-X_P^{t+1}\right|\cdot A^{+}\cdot L,\hfill & \text{otherwise}\hfill \end{array}$ (2)

式中， $X_{worst}$ 表示当前全局最差位置， $X_p$ 表示最优位置。A表示一个元素随机为1或−1的 $1\times d$ 矩阵，并且A满足 $A^{+}=A^T{\left(A{A}^T\right)}^{-1}$ 。当 $i>\frac{n}{2}$ 说明，第i个加入者没有食物，需要前往其它地方觅食，并更多的食物。

警戒者：

$X_{i,j}^{t+1}=\{\begin{array}{ll}{X}_{best}^t+\beta \cdot \left|X_{i,j}^t-X_{best}^t\right|,\hfill & \text{if}{f}_i>f_g\hfill \\ X_{i,j}^t+K\cdot \frac{\left|X_{i,j}^t-X_{worst}^t\right|}{\left(f_i-f_g\right)+\epsilon },\hfill & \text{if}{f}_i=f_g\text{}\text{}\text{}\text{}\hfill \end{array}$ (3)

式中 $X_{best}$ 是当前迭代次数中全局最优的位置。β和K是随机数，且 $\beta \sim N\left(0,1\right)$ ， $K\in \left[-1,1\right)$ 。 $f_i$ 是第i个麻雀的适应度值。 $f_g$ 和 $f_w$ 分别是截至目前为止迭代中的全局最佳和最差适应度值。ε的常数。K表示方向也表示步长。 $X_{best}$ 表示最好的位置。

2.2. 走时目标函数

地震震源定位研究，实质上就是通过数学方程得到地震发生的震源位置，主要通过检波器记录到的地震实际走时和地震理论走时，联系检波器的位置，建立一个目标函数，通过计算目标函数的最小值得到震源位置，也即是目标函数的最优解。以二维平面为例，过程就是某个 $t_0$ 时刻，某点 $s\left(x,y\right)$ 处发生了岩石破裂产生破裂力，以波形的方式传播出去，然后被某个检波器k接受到波形信息，记录时间为 $t_k^{obs}$ ，通常我们选择初至波。然后利用记录到的时间反演地震震源点的位置 $s\left(x,y\right)$ 和发震时刻 $t_0$ ，得到了如下目标函数：

$\phi \left(t_0,x,y,z\right)={{\displaystyle \underset{k=1}{\overset{Ng}{\sum }}\left|t_k^{cal}-t_k^{obs}\right|}}^p$ (4)

其中 $t_k^{cal}$ 表示地震波理论到达的时间，Ng表示接收到这个波形的检波器个数，p是范数。P一般为1或者2，本文选取2。又因为 $t_k^{cal}=t_0+T_k\left(m,x,y\right)$ ，式(4)化简为：

$\phi \left(t_0,x,y\right)={\displaystyle \underset{k=1}{\overset{Ng}{\sum }}{\left|t_0+T_k\left(m,x,y\right)-t_k^{obs}\right|}^p}$ (5)

其中m是速度结构模型，T表示算子。如果计算到的 $\left(t_0,x,y\right)$ 是真实解时， $\phi \left(t_0,x,y\right)$ 得到了最小值，就到得到了目标函数：

$\min \phi \left(t_0,x,y\right)={\displaystyle \underset{k=1}{\overset{Ng}{\sum }}{\left|t_0+T_k\left(m,x,y\right)-t_k^{obs}\right|}^p}$ (6)

式(6)就是所求的目标函数其中 $p=2$ 。

3. 仿真实验

本文为了验证本文方法的有效性及准确性，在二维空间同时布设了3个检波器，设起震时刻未知，震源所在位置为(200, 200)，各检波器空间所在位置及初至到时如表1所示：

Table 1. Geophone coordinates and first arrival time

表1. 检波器坐标及初至到时

理论模拟数据采用雷克子波模拟震源，得到如下的地震数据。

Figure 2. Simulated test waveform

图2. 模拟测试波形图

如图2所示，共计3个传感器获取模拟波形数据，通过STA/LTA方法自动拾取初至，分别得到3个传感器观测到时为2826.0 ms，2825.3 ms，2845.0 ms。将表1数据代入到SSA中，可以得出震源坐标，与真实震源坐标比较见表2。

Table 2. Solves the coordinate values

表2. 求解坐标值

从表2上看SSA方法比Geiger定位精度高，结果好。在本次实验进程中，迭代次数对定位精度也有着一定的影响，迭代次数过大会导致算法的寻优时间较长，过小又会致使算法陷入局部最优解。所以，一个合适的迭代次数既能减少优化时间，提高搜索效率，还能节省不必要的能量消耗。迭代次数对定位影响如图3所示。

由于SSA算法在迭代13次后就能够求得目标函数的最小值，这说明随着目标函数值的减小，定位结果也会越来越精确。

Figure 3. SSA iteration diagram

图3. SSA迭代图

4. 结论

SSA通过迭代搜索得到走时目标函数的全局最小值，既是定位误差最小的微地震震源坐标，从而达到工程中对定位精度的需求。一般而言，传感器越多，拾取到的时间精度越高，定位效果与精度就越好，从而削弱个别传感器的到时误差。Geiger方法在远距离定位中，如果迭代初值的选择不佳会影响定位精度，还会计算不到震源发生的位置。而SSA方法可以避免这种情况，因为它具有快速的收敛速度和较强的寻优能力，可以有效提高微地震定位的速度和精度，从而实现快速定位微地震事件的位置，预测地震发展趋势等方面。

然而，在实际微定位过程中还有很多方面的问题，如地质构造速度的精确程度，纵横波的初至拾取准确度等方面，在今后进行微地震定位相关的工作时，如何做到全自动，高精准度，更加快速地定位，从而构建地质地层结构，以达到地球可视化的程度。