引言

地震是一种常见的破坏性极大的自然灾害,严重损害着灾区人民的生命和财产安全. 在过去的数十年中, 国内外相关领域的研究人员早已开始关注地震前产生的异常, 以期对地震未来的发展趋势进行预测. 20世纪80年代前苏联学者Gorny等^[1]首次报道了在中亚及东地中海地区,许多中强震前出现大面积卫星热红外辐射增强现象. 其他学者也提出了许多新的震前异常研究成果^[2]. NĚMEC等^[3]使用统计方法研究地震前低频电磁波数据, 发现在地震前短时间内电磁波的强度会出现异常, 这种与多种因素有关. 张元生等^[4]选用静止卫星红外遥感亮温资料证明在大地震发生之前亮温变化存在明显的特征周期和振幅以及热异常分布区域. Kong等^[5]使用几何移动平均鞅算法分析NOAA卫星遥感OLR数据的变化过程, 挖掘汶川地震和芦山地震前OLR数据的异常. Yeha等^[6]采用多个连续GPS基准站建立一个GPS网络实时监测台湾南部地壳形变, 发现GPS信号在地震发生前很短时间内出现了系统性扰动. Kuo等^[7]应用了改进的岩石圈-大气-电离层(LAI)耦合模型计算日本2011年 Tohoku-Oki 9.0级大地震发生前的电离层总电子含量, 并将计算结果与所报道的TEC观测结果进行比较, 以研究地震前电离层TEC异常. 文献[8]中使用涡流场计算和小波变换检测地震前OLR数据的异常. 文献[9]中使用模糊神经算法用于检测震前电磁异常. 文献[10]中提出了基于误差和关键点的自顶向下分段算法以及基于时间邻域的局部异常因子的分析方法对地震前兆观测数据进行异常挖掘. 文献[11]中提出地面持续增温是一种相当普遍的临震前兆增温异常, 这种前兆异常的时空动态,可实时地反映于卫星热红外图象上. 文献[12]中提出使用一种非线性过滤器HMMs来提取GPS数据中可能与地震有关的异常信号, 进而达到预测地震的目的. 文献[13]提出用小波变换来识别并分析数据中的奇异点与地震的关系. Marzocchi等^[14]使用贝叶斯估值方法分析地震数据, 构建地震整体预测. 钱国红^[15]利用中国大陆GPS连续观测站资料,获取了2011年3月11日日本9.0级地震造成的连续站同震位移, 通过对时间序列分析发现,有明显同震位移的连续站,震前水平方向的运动速度都有放缓的趋势. 他们认为这可能是一种地震形变前兆现象. 卫星热红外异常是一种临震预报的新指标. Ouzounov等^[16]通过系统地应用卫星数据分析技术, 将卫星数据分析技术应用于大地震后的图像, 从而验证TIR异常与已知的地震相关联. 从上面可以看出, 地震发生前会产生很多异常信号, 通过对地震前兆异常的研究与分析, 寻找其中蕴含的地震前兆异常信息对预测地震, 防震减灾具有重要意义.

本文的贡献主要在于: 虽然卫星遥感技术已经广泛地应用于地震科学研究, 但由于热红外遥感数据是通过卫星来监测地面的变化, 通常会受到云层等各种因素的影响, 卫星能够监测到的信号往往比较弱. 缺乏有效的热红外遥感数据处理技术来提取地震异常相关信息以及探讨震前热异常内部机制, 导致大多数的热红外卫星遥感数据没有被充分利用. 虽然已经有不少学者开始运用数据挖掘中的异常检测算法来分析地震异常, 但是通过量子漫步算法来分析震前异常的研究相对较少, 这种算法具有鲁棒性, 是一种可行的算法. 因此, 为了能够有效的识别地震热红外(射出长波辐射)异常, 本文基于数据挖掘技术的基本原理, 结合Martingale理论^[17], 提出一种量子漫步概率算法来识别热红外异常特征, 然后在经过一定的信号放大来发现异常. 针对发生在中国比较具有代表性的汶川和芦山地震, 运用该算法, 对地震前后震中附近OLR数据、P值、CD值进行分析研究. 并且, 通过实验将该算法扩展到最近十年全球发生的8.0级及以上地震, 验证该算法的有效性与可靠性. 结果表明, 量子漫步算法是一种具有鲁棒性的OLR异常变化识别方法. 同时, 它对将数据挖掘算法运用于地理科学研究领域具有重要的借鉴意义.

1 相关数据

本文研究的对象是发生在中国的2008年5月12日汶川大地震和2013年4月20日芦山大地震,以及从2005年6月到2014年9月期间全球发生的8级及以上中的16个地震. 具体地震相关数据如表1所示.

2 研究方法

量子漫步作为量子信息领域的一个重要分支, 是经典随机漫步向量子场景转变的概括. 量子漫步大致分为离散量子漫步和连续量子漫步两类. 本文所用的是在直线上的离散量子漫步算法挖掘震前异常. 在直线上的离散量子漫步可以描述为一个维度 ${H_P}$ 无限但却可数的希尔伯特空间^[18]. ${H_P}$ 相应的计算基础是 $\left\{ {{{|n}}\rangle |n \in Z} \right\}$ , 其中Z表示整数集. 硬币空间 ${H_{{C}}}$ 是一个二维的希尔伯特空间, 它的计算基础是 ${\rm{|}}0\rangle $ 和 ${\rm{|}}1\rangle $ , 分别对应量子漫步运动的向左走向右走的方向.

量子漫步异常识别方法具体实现步骤如下:

首先, 由硬币算子 ${U_C}$ 和转换算子 ${U_S}$ 组成演化算子 $U$ , 且有

$U = {U_S} \cdot ({I_P} \otimes {U_C})$

(1)

在每一步中, 硬币粒子都经过一次哈达玛变换, 然后, 运动粒子根据硬币粒子的状态向对应方向运动. 具体漫步规则如图1(a)~(c)所示, 从 ${{t}}$ 时刻开始, 粒子在手性“右”(或“左”), 经过哈达玛变换,粒子现在处于“左”手性态和“右”手性态的相等叠加态, 粒子随之移动到 $t + 1$ 时刻的状态.

其次, 经过 $t$ 步之后, 整个量子漫步的状态为:

$|{\psi _t}\rangle = {(U)^t}|{\psi _0}\rangle $

(2)

其中, ${\rm{|}}{\psi _0}\rangle = |{P_0}\rangle |{C_0}\rangle $ , $|{P_0}\rangle $ 表示漫步的初始状态, ${\rm{|}}{C_0}\rangle $ 表示硬币粒子. 选择当前监测数据点前 $ws$ 个数据点, 通过公式(2), 计算每一个数据点的量子漫步概率特征, 分别计算获得当前数据点根据量子漫步游走的概率分布.

再次, 结合Martingale理论^[17],对前两个步骤获得的特征序列做如下处理,定义数据点的概率变化程度( $CD$ ):

$CD_n^{(\varepsilon )} = \frac{{\sum\limits_{k = 1}^{100} {\prod\limits_{i = 1}^n {(\varepsilon \hat p_{i,k}^{\varepsilon - 1})} } }}{{100}}$

(3)

其中, 通过计算 $\prod\limits_{i = 1}^n {\varepsilon \hat p_{i,k}^{\varepsilon - 1}} $ 得到每个数据点的所对应的Martingale值, 根据文献[19], 令 $\varepsilon = 0.82$ . $\hat p$ 计算如下:

对于无标签的时间序列数据集 $Z = \{ {z_1}, \cdots ,{z_{i - 1}}\} $ 和当前计算的数据点 ${z_i}$ . 异常度评分该点跟其他点的区别定义如下:

${s_i} = s(Z,{z_i}) = \left\| {{z_i} - m} \right\|$

(4)

其中, $\left\| \cdot \right\|$ 定义的距离度量, $m$ 是集合 $Z \cup \left\{ {{z_i}} \right\}$ 中点的均值.

在上面定义的基础上,我们定义 ${\hat p_{i,k}}$ 计算如下:

${\hat p_{i,k}}(Z \cup \{ {z_n}\} ,{\theta _n}) = \frac{{\# \{ j\left| {{s_j} > {s_i}} \right.\} + {\theta _{i,k}}\# \{ j\left| {{s_j} = {s_i}} \right.\} }}{i}$

(5)

其中, ${\rm{\# }}\left\{ {} \right\}$ 是一个用于返回满足给定条件的样本的数量的函数, ${\theta _{i,k}}$ 为(0,1)范围内选取的参数, $i = 1,2, \cdots, n$ , ${s_j}$ 是点 ${z_j},j = 1,2, \cdots, i$ 的异常度评分,通过等式(4)计算得到 ${S_i}$ , 初始化值 $CD_0^{(\varepsilon )} = {\rm{1}}$ .

通过等式(3)的计算可以看出, $CD$ 值越大, 说明地震OLR异常越明显.

最后, 用上一步计算方法计算每一个数据点的异常度. 为了避免 $CD$ 值在短时间内出现大波动, 导致结果不可控, 本文算法中设置了一个阈值 $h$ 作为停止参数, 当某一个数据点满足下面条件时, 说明异常度很大,重新初始化步骤二, 继续计算剩余点的异常度:

$C{D_n} \geqslant h$

(6)

综上所述, 基于量子漫步算法的异常挖掘算法的流程图如图2所示.

3 结果与分析

(1)汶川地震和芦山地震对应的原始射出长波辐射OLR数据分别如图3(a)和图4(a)所示. 量子漫步算法中, 窗口大小 $ws$ 取值为30–45的每个数, 然后求均值,以降低随机因素对结果造成的影响. 实验时, 根据文献[5], 将阈值 ${{h}}$ 设置为1000. 研究的数据时间周期为一年, 即从发生地震上一年的9月1日到第二年的9月1日. 黑色竖线表示地震发生时间.

通过量子漫步算法进行初步处理OLR数据, 计算得到P值, 如图3(b)和图4(b)所示. 进一步处理特征序列, 计算数据点的概率变化程度 ${{CD}}$ , 得到的异常值结果如图3(c)和图4(c)所示.

从图3(c)和图4(c)可以看出, 在距离地震发生时间较远的地方, ${{CD}}$ 值趋近于0, 几乎没有任何变化, 说明OLR数据比较稳定. 而临近地震前后, 两个地方的 ${{CD}}$ 值均出现了显著变化, 说明汶川地震和芦山地震发生前热红外异常就已经出现, 且在地震后仍有异常. 这进一步证实了OLR数据中含有大量地震异常信息.

结合图3(b), 图3(c), 图4(b)和图4(c)可以看出, 当地震的P值图像在短时间内有连续且密集的上下波动的区间内会出现较大的异常值, 即 ${{CD}}$ 值的图像会在该处出现较大幅度的波峰. 而图3(a)和图4(a)的原始OLR数据图像并未出现明显的异常趋势, 这证实了基于量子漫步的异常挖掘算法提取OLR数据异常的有效性.

由于汶川和芦山两个地震的震中位置相对较近, 成因相似, 因此所提取的异常值 ${{CD}}$ 变化规律比较类似, 这从侧面说明本文提出的基于量子漫步算法的异常挖掘算法的可靠性.

(2)为了进一步说明本文所提出的基于量子漫步算法的异常挖掘算法的可靠性和有效性, 现在扩展该算法的适用范围, 将该算法运用于分析2005年6月到2014年9月全球发生的16个8.0级及以上地震的异常数据. 受文章篇幅限制, 本文只展示其中部分地震的实验结果, 实验结果如图5~图10所示.

通过研究近10年来的时间长度为一年的OLR时间序列, 综合以上图5~图10所示图像可以看出, 各个地方CD值图像的变化趋势有相似的规律, 在距离地震发生时间较远的地方, CD值趋近于0, 几乎没有任何变化, 而在地震前后, 各个地方的CD值均出现了显著变化, 即临近地震前后会出现较大异常, 而且越大的地震异常越明显. 然而, 如果仅从地震区域的原始OLR数据图像无法直观地看出是否存在异常. 这也说明本文提出的基于量子漫步算法的异常挖掘算法几乎适用于所有地震, 具有可靠性和有效性.

(3)为了说明P值和CD值在震前出现异常和地震发生之间的因果关系, 本文针对汶川和芦山地震分别设计了对比实验.

首先, 将汶川地震与其他年份(2005年到2017年)汶川地区的P值和CD值图像进行比较. 实验结果表明, 在总共12年数据中, CD值大于100, 且小于200的有8个; 大于200且小于300的有1个; 大于300, 且小于400的有2个; 大于400的有1个. 在这12个数据中, 其中有10个异常发生的时间与汶川或其周边发生的5.0级以上地震发生时间相吻合; 出现误报的仅2个(2007和2016), 即出现异常后, 地震真正发生的概率为83.3%, 误报率为16.7%.

其次, 将芦山地震与其他年份(2005年到2017年)芦山地区的P值和CD值图像进行比较. 实验结果表明, 在总共12年数据中, CD值大于100, 且小于200的7个; 大于200且小于300的有2个; 大于300, 且小于400的有2个; 大于400的有1个. 在这12个数据中, 其中有10个异常发生的时间与汶川或其周边发生的5.0级以上地震发生时间相吻合; 出现误报的仅2个(2007和2016), 即出现异常后, 地震真正发生的概率为83.3%, 误报率为16.7%.

汶川和芦山2007年都出现大于100, 且小于200的异常, 及2016年都出现大于200且小于300, 笔者猜测可能是由于气候、温度或人为活动等原因引起的异常.

(4)为了说明本文提出的算法的优越性, 本文通过实验, 比较本文提出的基于量子漫步算法的异常挖掘算法与经典随机漫步算法提取地震异常特征的方法. 结果如下图所示. 其中, 实线曲线(图中左曲线)为100步量子漫步概率分布, 虚线曲线(图中右曲线)为随机漫步概率分布.

从图11中可以看出, 由于硬币和粒子的量子干涉效应, 在多步以后, 量子漫步算法展现出比经典随机漫步更好的一些性质. 在执行N步之后, 与经典随机漫步不同, 离散量子漫步的概率分布并不是高斯分布, 它的分布偏向左边, 比经典随机漫步收敛性快, 显然本文提出的基于量子漫步算法的异常挖掘算法更具有优势.

4 结束语

本文结合Martingale理论, 提出了一种基于量子漫步算法的震前异常挖掘算法, 并运用该算法分析汶川和芦山两个地震的地震前后, 震中附近OLR数据发生的显著异常变化, 挖掘有用的震前异常信息. 从结果中可以看出, 即使原始OLR数值没有明显异常, 在地震前后一段时间内其对应的 ${{CD}}$ 值也会出现较大的变化. 相对于其他传统的算法, 使用量子漫步算法计算的 ${{CD}}$ 值的变化更能客观形象地刻画出地震前后数据的异常程度. ${{CD}}$ 值越大, 则表明当前数据点的异常变化越大. 这不仅进一步证实了OLR数据中蕴含着震前异常信息, 也为使用量子漫步算法分析OLR数据, 计算 ${{CD}}$ 值预测地震提供了很大的可能. 为了进一步验证该算法的有效性和可靠性, 本文应用该算法分析近10年左右的全球8.0级地震的OLR, 挖掘地震异常, 结果可观, 扩展了该算法的适用范围. 然而, 基于目前的数据量和研究水平, 应用该算法研究地震异常取得的成果很少, 使用量子漫步算法预测地震仍任重道远. 本文提出的基于量子漫步算法的异常挖掘算法所对应的 ${{CD}}$ 值的异常已经找出了一些初步规律, 但还需要收集更多的地震数据, 进行更完善的实验进一步挖掘 ${{CD}}$ 值的变化程度与地震相关参数(震源深度、震级、GPS等)之间的关联性问题, 这也是本文下一步的研究方向.