微震定位是地下空间目标定位领域中的核心技术, 也是解决高价值弹药地下炸点定位以及侵彻轨迹测量等军事问题的关键; 同时, 它也是实现煤矿勘察, 文物防盗、工程爆破、地质监测等民用问题的重要手段^[1–3].

传统的地下浅层微震定位主要借鉴深层次、大深度的天然类地震定位方法. 在这个阶段, 主要使用以下3种定位模型: (1)德国物理学家Geiger提出的走时反定位模型^[4], 国内外学者在此基础上不断优化和完善, 如Waldhauser和Ellswotrh提出双重残差定位法^[5], Zhang和Thurber在双差定位的基础上提出了双差层析成像法^[6], 该方法能够在水平层状速度模型对震源群重定位的同时可以得到改进的速度结构. 但是该类方法主要应用于水平层状速度模型, 但群波混叠严重, 时差数据提取误差大, 其定位精度有限; (2)以波场偏振理论为核心的DOA (Direction Of Arrival)定位模型, 蒋鑫将上述模型成功应用于某煤矿试验区的爆破及检测工作, 该方法无需求取初至时间, 采用少量观测节点即可实现震源的快速定位^[7], 但是在比较复杂的地质结构中尤其存在一些具有强反射特性的地质界面时, 表层入射角度无法表征震源和传感器接收点之间的真实射线路径, 因此存在定位假象的问题; (3)针对未知地质结构条件下的定点群定位问题, Cross提出多震源与速度结构联合迭代反演定位(SSH)模型^[8], Aki等在SSH基础上, 将地层介质横向非均匀速度结构网格化, 提出了三维速度结构与震源联合反演理论^[9]. Kissling等通过研究震源位置和一维速度模型耦合问题, 提出了确定最佳一维速度模型的Velest方法, 定位精度得到显著提高^[10]. 在国内, 赵仲和采用SSH算法建立新的地震波速度模型MDBJ81, 并将该模型用于SSH方法, 提高了北京台网的测定能力^[11]. 刘福田引入正交投影算子实现参数分离, 提高联合反演的效能, 减少了运算量^[12]. 但是由于震源的定位精度受速度场建模精度、震源次数的影响, 该方法主要适用于未知地质结构条件下的定点群定位问题. 而地下浅层微震定位数量较少, 炮射投影数据有限, 定位精度难以保证.

针对以上问题, 本文通过将浅层偏振信息与深层走时技术相结合, 建立基于时差和偏振角的混合定位方程, 提高由单纯时差定位引起的定位误差, 以多事件相关信息为基础, 运用改进后的粒子群算法对所建立的混合定位模型进行解算, 提高波场信息的冗余度, 在实现地下空间震动目标的高精度定位的同时也为其提供夯实的理论支持.

1 基于走时-偏振混合定位模型的构建方法

定位模型的构建对定位结果的影响至关重要, 本文通过利用传统走时定位模型并在此基础上, 联合深层偏振信息, 建立走时-偏振混合定位模型^[13,14]. 三维定位示意图如图1所示.

根据时差定位原理^[15], 我们可以建立以下方程组:

$\left\{ {\begin{array}{*{20}{l}} {r_0^2 = {{(x - {x_0})}^2} + {{(y - {y_0})}^2} + {{(z - {z_0})}^2}} \\ {r_i^2 = {{(x - {x_i})}^2} + {{(y - {y_i})}^2} + {{(z - {z_i})}^2}} \\ {{r_{0i}} = {r_0} - {r_i} = c{t_{0i}}} \end{array}} \right. $

(1)

其中, $\left( {x,y,z} \right)$ 为地下震动源坐标, 各传感器节点坐标分别为 $\left( {{x_i},{y_i},{z_i}} \right)\;\left( {i = 0,1,2,3, \cdots ,n} \right)$ , ${r_i}(i = 0,1,2,3, \cdots ,n)$ , ${t_{0i}}$ 分别为震源到各传感器节点的距离和时差数据.

如图1所示, 基于传感器 $i$ 与震源俯仰角 ${\gamma _i}$ 之间的几何约束关系, 我们可以得到:

$\tan {\gamma _i} = \frac{{z - {z_i}}}{{{r_i}}}$

(2)

其中, ${r_i} = \sqrt {{{(x - {x_i})}^2} + {{(y - {y_i})}^2}} ,\;i = 1,\;2,\;\cdots, n$ .

将式(1)和式(2)联立方程组可得:

$\left\{ {\begin{aligned} & {\tan {\gamma _i} = \frac{{z - {z_i}}}{{r{}_i}}} \\ & {r_0^2 = {{(x - {x_0})}^2} + {{(y - {y_0})}^2} + {{(z - {z_0})}^2}} \\ & {r_i^2 = {{(x - {x_i})}^2} + {{(y - {y_i})}^2} + {{(z - {z_i})}^2}} \\ & {{r_{oi}} = {r_0} - {r_i} = c{t_{0i}}} \end{aligned}} \right.$

(3)

其中, $i $ =1, 2, 可改写为:

$ \left\{ {\begin{array}{*{20}{l}} {z = ({r_0} - {r_{01}})\tan {\gamma _1} + {z_1}} \\ {({x_0} - {x_1})x + ({y_0} - y{}_1)y + ({z_0} - {z_1})z = {k_1} + {r_0}{r_{01}}} \\ {({x_0} - {x_2})x + ({y_0} - {y_2})y + ({z_0} - {z_2})z = {k_2} + {r_0}{r_{02}}} \end{array}} \right. $

(4)

其中, ${k_i} = \dfrac{1}{2}\left[ {r_{0i}^2 + (x_0^2 + y_0^2 + z_0^2) - (x_i^2 + y_i^2 + z_i^2)} \right]$ , $i = 1,2$ , 可改成矩阵形式为:

$AX = B$

(5)

其中, A是可逆矩阵.

$A = \left\{ {\begin{array}{*{20}{c}} 0&0&1 \\ {{x_0} - {x_1}}&{{y_0} - {y_1}}&{{z_0} - {z_1}} \\ {{x_0} - {x_2}}&{{y_0} - {y_2}}&{{z_0} - {z_2}} \end{array}} \right\}$ ,

$B = \left\{ {\begin{array}{*{20}{c}} {({r_0} - {r_{01}})\tan {\gamma _1} + {z_1}} \\ {{k_1} + {r_0}{r_{01}}} \\ {{k_2} + {r_0}{r_{02}}} \end{array}} \right\}$ .

求解方程组可得

$X = {A^{ - 1}}B$

(6)

2 基于粒子群优化的混合定位解算方法 2.1 粒子群算法原理

在粒子群算法中, 将群体中每个飞行个体视为 $D$ 维搜索空间中的一个点, 称之为“粒子”. PSO初始化为一群随机粒子并按照一定的速度在搜索空间中飞行, 然后粒子根据当前的最优粒子在解空间中搜索最优解^[16,17]. 在每一次迭代中, 粒子都是通过跟随两个“极值”来更新自己, 一个是粒子自身找到的最优解, 称为个体极值 ${p_{\rm{best}}}$ ; 另一个极值是整个群体找到的最优解, 称为全局极值 ${g_{\rm{best}}}$ .

假设粒子种群规模大小为 $m$ , 目标搜索空间维数为 $D$ , 第 $i$ 个飞行粒子的位置向量表示为 ${x_i} = [{x_{i1}}, {x_{i2}}, \cdots ,$ ${x_{iD}}]$ , 飞行速度表示为 ${v_i} = [{v_{i1}},{v_{i2}}, \cdots ,{v_{iD}}]$ , 第 $i$ 个飞行粒子搜索到的最优位置为, 整个粒子群搜索到的最优位置为 ${p_g} = [{p_{g1}},{p_{g2}}, \cdots ,{p_{gD}}]$ , 则粒子群速度和位置更新过程可表示为:

$ {v_i}(n + 1) = \omega {v_i}(n) + {c_1}{r_1}({p_i} - {x_i}(n)) + {c_2}{r_2}({p_g} - {x_i}(n)) $

(7)

${x_i}(n + 1) = {x_i}(n) + {v_i}(n + 1)$

(8)

其中, $i = 1,2, \cdots ,m$ , $\omega $ 为惯性权重, ${c_1}$ 和 ${c_2}$ 分别是学习因子, ${r_1}$ 和 ${r_2}$ 分别是均匀分布在 $[0,1]$ 区间的随机数; 粒子在目标搜索空间内不断更新粒子的位置和速度, 搜索个体极值和全局极值, 直到达到设定的迭代终止条件为止.

相较于其他传统定位算法, 粒子群优化PSO (Particle Swarm Optimization)算法最大的特征是收敛速度快且具有更强的全局优化能力. 在一定程度上提高了收敛速度并且能够保证粒子可以处于全局最优的状态. 但是, 由于种群规模不大以及种群质量不好等问题可能造成PSO算法在搜索后期容易出现局部最优, 甚至收敛停滞.

针对这一不足之处, 本文在走时定位的基础上, 一方面结合偏振角信息, 充分利用传感器测量信息, 通过混合参数测量对传统的PSO算法进行种群策略方面的改进. 另一方面, 通过在传统PSO的基础上引入遗传算法中的交叉变异机制, 可以防止PSO算法在搜索的后期陷入局部最优的状态.

2.2 种群改进策略

如图2所示通过极化分析, 利用极化度最高的两个高匹配节点的信息构建目标函数, 其余样本信息利用式(6)构建目标解空间, 并加入变异粒子群, 扩展初始种群的规模并增加其多样性. 这种改进方式与传统PSO所产生的初始粒子群方法相比, 大大提高了收敛速度并提高了解的正确性.

2.3 基于交叉变异的混合粒子群优化算法

通过引入遗传算法中的交叉变异机制, 可以防止PSO算法在搜索的后期陷入局部最优的状态, 导致收敛停滞并且在一定程度上提高粒子群算法的全局收敛性.

(1)自适应交叉策略

根据式(6)所构建的解空间, 随机选取5个解, 通过粒子范围扩展成种群数量为10的新生粒子群 $\bar Y$ , 并将解空间里其他的解构成初始粒子群 $\bar X$ , 对所构建的粒子种群 $\bar X$ 和 $\bar Y$ 通过计算惯性权重以及判别条件按照一定的概率进行交叉操作, 从而更新粒子种群, 产生新的搜索空间^[18]. 则交叉后得到的子代个体为:

$\left\{ {\begin{array}{*{20}{c}} {{{\bar X'}} = {P_c}\bar X + (1 - {P_c})\bar X} \\ {{{\bar Y'}} = (1 - {P_c})\bar X + {P_c}\bar Y} \end{array}} \right.$

(9)

式中, 交叉概率 ${P_c} = U(0,1)$ . 当它不是一个固定值时, 此时我们称之为为非一致性交叉机制, 这种机制不仅有助于增加种群规模的多样性, 而且还增强了种群的全局搜索能力.

当适应度值在算法进程中发生变化的时候, ${P_c}$ 可以根据适应度的变化自动调节自己的参数值, 当种群适应度趋于相同时, ${P_c}$ 增加, 而当种群适应度比较分散时, ${P_c}$ 减小; 交叉概率越小, 说明个体适应度高于种群平均适应度, 应该得到保护; 否则, 个体趋于集中, 应该被淘汰. 根据上述指导思想, 我们设计的自适应交叉概率 ${P_c}$ 可表示为:

$ {P_c} = \left\{ {\begin{aligned} & {{P_{c\max }} - \frac{{({P_{c\max }} - {P_{c\min }})(f - {f_{\rm{avg}}})}}{{{f_{\max }} - {f_{\rm{avg}}}}},f \ge {f_{\rm{avg}}}} \\ & {{P_{c\max }},f < {f_{\rm{avg}}}} \end{aligned}} \right. $

(10)

式中, ${P_{c\max }}$ 和 ${P_{c{\rm{min}}}}$ 分别为交叉概率的最大值和最小值, $f$ 为适应度值, ${f_{\rm{avg}}}$ 为整个种群的平均适应度值, ${f_{\max }}$ 为种群适应度的最大值.

(2)混合纵向变异机制(VMAPSO)

混合纵向变异机制(VMAPSO)包含均匀分布变异和高斯分布变异^[19]. 相对于传统的变异机制而言, 通过引入VMAPSO可以在算法出现收敛停止的现象时确保各粒子能够朝着新的方向继续搜索, 通过对全局极值进行变异来提高粒子群优化算法跳出局部最优的能力.

由于在该自适应混合纵向变异机制中, 两者是以一定的比例相互交替进行, 因此在混合纵向变异机制中假定均匀分布变异的概率为 ${p_u}$ , ${p_u}$ $ \in \left[ {0,1} \right]$ , 并且高斯分布变异的概率为 $1-{p_u}$ 因此, 对进行交叉后的粒子种群进行混合纵向变异的过程可以描述为:

$ \begin{aligned}[b] & {\rm{if}}\;\;rand \le {p_u}\\ & \;\;\;\;\;\;\;{x_{id}} = {x_{\min }} + rand \times ({x_{\max }} - {x_{\min }})\\ & {\rm{else}}\\ & \;\;\;\;\;\;\;{x_{id}} = {x_{id}} \times (1 + Gaussian\;(\sigma ))\\ & {\rm{end}} \end{aligned} $

(11)

其中, ${x_{\max }}$ 为均匀变异的最大边界值, ${x_{\min }}$ 最小边界值, $\sigma $ 为高斯变异的方差.

3 算法实现流程

通过在传统的PSO算法基础上引入交叉变异机制, 对构建的走时-偏振混合定位模型进行定位解算以及迭代寻优.

在图3中, 算法的主要步骤为:

Step 1. 首先对粒子群算法的初始参数进行设定, 包括(种群规模、学习因子、循环次数、最大迭代次数、均匀变异概率等), 将粒子解空间分为两部分粒子群, 即初始粒子群和新生粒子群;

Step 2. 计算粒子的适应度并实时更新粒子位置信息;

Step 3. 判断算法是否达到收敛条件, 如果达到直接输出定位结果, 否则判断是否达到最大迭代次数, 如达到则但会重新计算粒子群, 否则向下执行Step 4;

Step 4. 依据自适应交叉概率将初始粒子群 $\bar X$ 和新生粒子群 $\bar Y$ 进行交叉并更新粒子群;

Step 5. 重新计算更新后的粒子群的惯性权重, 并实时更新粒子的个体极值和全局极值;

Step 6. 在此判别是否达到收敛条件, 如果不满足执行下一步, 否则输出;

Step 7. 将更新后的粒子群进行混合纵向变异操作, 引入变异粒子群;

Step 8. 重新计算引入变异粒子群的粒子种群惯性权重以及实时更新各极值位置;

Step 9. 判断算法是否达到收敛条件, 如果不达到条件执行Step 7, 否则执行Step 10;

Step 10. 将解的位置进行结果输出, 即目标位置并结束操作.

4 实验分析

为对震源仿真, 本文选择5个分布式传感器节点构成的阵列对算法进行试验, 定位区域为, 其中选取传感器阵列中的一个节点设为主站, 其坐标为(4.648, 4.528, –3.59), 其他4个传感器节点作为基站, 坐标分别为(7.634, 1.0315, –3.5), (7.112, –7.018, –3.65), (7.650, 7.530, –3.28), (–12.813, 12.565, –3.00), 创建震源位置坐标为(0, 0, –5), 各基站及震源坐标分布图如图4所示. 在小区域范围内, 对地下震动目标定位, 介质可看做是均匀的, 地震波传播方式以纵波为主, 根据试验的实际爆破数据分析, 地下介质的传播速度可以设定为640 m/s; 为保证仿真数据的真实性, 在仿真实验中对所获取的各传感器基站间的时差、偏振角信息加入随机的干扰噪声, 且 $\delta $ =0.5.

对粒子群算法的参数进行设置: 种群大小规模设置为20, 维数为4, 最大迭代次数为200, 粒子飞行速度范围为, 粒子搜索范围是5, 学习因子和为2, 惯性权重的最大值和最小值分别为0.6和0.3.

通过频谱分析, 计算出相应的时差信息, 并结合传感器偏振角信息, 通过Taylor法计算出相应的坐标位置信息构成目标解空间, 在产生的解空间中随机选择五组解, 将其扩展到10个粒子作为新生粒子群, 并将解空间中剩余的解作为初始粒子群, 交叉概率的最大值和最小值设为0.9和0.4, 均匀变异的概率选择为0.03; 结合上述参数, 运用本文所改进的粒子群算法对混合定位模型进行解算并寻优, 结合其他寻优算法并将各算法得出的适应度值进行对比, 结果如图5所示.

由图5可知, 引入偏振角信息并用改进后的粒子群算法进行解算时, 在50代的时候开始收敛, 适应度值接近于0.031; 而运用传统的PSO定位时, 迭代到大约70次时, 发生早熟并且进入局部收敛. 传统的时差定位, 在120次开始收敛; 引入偏振角信息后, 在100次收敛.

由表1对比可知, 结合走时和偏振信息后, 运用本文所提出的定位算法解算定位模型时结果明显比传统的粒子群定位精度大幅度提高并且误差降低了30%.

5 结论与展望

为了解决传统的走时定位模型在地下浅层定位空间中误差大的问题, 本文创新性地提出一种高精度的地下震源定位算法. 分析实验结果, 得到以下结论.

(1)由表1对比可知, 结合走时和偏振信息后, 运用改进后的定位算法解算结果明显比传统的粒子群定位精度大幅度提高并且误差降低了30%.

(2)利用多事件相关信息, 建立走时-偏振混合定位模型, 运用改进的PSO寻优算法对混合定位模型进行迭代寻优, 这样不仅大大地减少了算法陷入局部最优的风险而且提高收敛速度. 经过试验和仿真, 证明本文所提出的算法可有效地提高定位精度.