过去几十年, 通过研究大量震例的遥感数据, 研究者们发现大地震发生前常有各种异常现象出现. 如地表潜热通量 (Surface Latent Heat Flux, SLHF)^[1], 热红外 (ThermaL Infrared, TIR) 异常^[2], 电离层总电荷量(TOtal Electron Content, TEC)^[3]和射出长波辐射(Outgoing Longwave Radiation, OLR)^[4]等, 在地震发生之前都会出现明显的变化. Qin 等^[5]分析2010年新西兰7.1级地震发生前后3个月内的SLHF数据的变化, 发现在地震发生前一个月左右, 震中的东北方向出现了一个孤立的高达160 W/m²的SLHF正异常. 结合GPS位移观测和构造背景, 确定SLHF异常可能与地壳运动时上升的地下热物质有关, 从而解释了震中东北部, 北岛中心和南岛西南部的局部温度升高现象. 张璇等^[6]分析了尼泊尔8.1级地震前热红外遥感信息, 发现震前出现明显的热红外亮温异常变化. Aggarwal^[7]分析了2010年4月13–14日青海地区发生的地震中电离层中总电荷数的变化情况, 发现电离层总电荷数的异常变化或许可以作为这次地震相关前兆特征. Ouzounov等^[8]研究OLR数据, 并发现震前会有连续的异常变化发生.

OLR是地球-大气系统透过大气顶部向宇宙空间发射的主要波长在4～120 μm的辐射. 目前很多研究者研究OLR数据与地震前兆的关系, 并取得了一定的成果. Kong等^[4]采用基于Martingale理论的方法分析了汶川和庐山两次地震OLR数据中的异常信息, 发现在汶川地震和芦山地震发生前后OLR数据都会显著增大, 并在震后一段时间后恢复正常. Xiong和Shen^[9]提出了一种统计学方法, 分析了从2007年9月1日至2015年5月23日, 约3376个震例的OLR异常信息, 统计结果表明震前靠近地震中心的地方有异常现象出现, 并且这些异常与该地区即将发生的地震有关. 康春丽等^[10]提出用涡度场的方法分析昆仑山口西8.1级地震前后长波辐射强度OLR 的变化, 结果显示, 地震活动过程中, 不论是卫星遥感长波辐射场强度, 还是亮度温度都出现了较明显的异常变化. 孙珂等^[11]采用RST (Robust Satellite Technique) 算法和异常指数算法, 分别对2015年4月25日和5月12日尼泊尔发生的两次大地震前后, 卫星长波辐射变化特征进行分析, 发现运用RST算法, 两次地震前后都没有检测到明显异常变化, 而运用异常指数算法, 两次地震前都检测到明显的热红外异常. Xiong等^{[12, 13]}使用小波变换和小波最大值的空间/时间连续性分析来分析连续OLR数据. 他们从OLR数据中发现了与地震前兆相对应的奇点. 这些研究成果都表明了在OLR数据中隐藏着重要的前兆信息, 地震前后, OLR数据会出现明显的异常变化. 分析OLR的变化与地震之间的关系有助于我们进一步了解地震相关机理, 给未来的防震减灾工作提供启示.

然而目前基于OLR数据的相关工作, 大多只是基于地学领域的方法分析, 算法上有一定的局限性. 并且大多只是考虑地震中心所在区域产生的异常, 而对周边区域的影响和历年背景环境产生的影响, 以及季节性影响并未全面考虑, 而这些因素对结果都会产生较大影响. 基于以上两点考虑, 本文提出一种基于OLR数据的震前异常分析方法-时空分析法, 结合统计学中的Martingale理论^[14-16], 聚类分析, 随机漫步等计算机相关算法, 并通过汶川地震, 日本地震, 以及芦山地震3个震例详细阐述该算法的分析方法及实验结果, 探索OLR数据中隐藏的前兆信息.

1 相关数据

由于2008年5月12日汶川8级地震是10年来中国国内发生的最大一次地震, 而2013年4月20日芦山7级地震的震中区域与汶川地震在同个研究区域中, 各种地理特性比较相似. 2011年3月11日发生的日本9级大地震是近年来全球范围内的一次特大地震. 这3个震例比较具有代表性, 故本文采用时空分析法分析这3个震例的OLR数据中隐藏的前兆信息. OLR数据由美国国家海洋和大气管理局(National Oceanic and Atmosphere Administration, NOAA)卫星监测^[17]. 具体震例信息如表1所示.

2 研究方法

本文提出的基于OLR数据的时空分析法分为6个部分: 1) 计算相对于背景数据的时间维度变化增量; 2) 计算空间维度相关变化增量; 3) 提取随机漫步概率特征, 4) 获得概率特征序列; 5) 计算各个点的异常度, 统计异常概率; 6) 计算当前点的异常值, 与预设阈值进行比较, 判断异常是否发生, 若发生, 则异常标记为1, 并重启算法. OLR数据集是将全球按照 $2.5^\circ \times 2.5^\circ $ 分辨率划分为 $72 \times 144$ 个单元. 由于OLR数据分为日间数据和夜间数据两类, 故全球每一天的数据可用144×144单元格表示. 若用 $Data\_All$ 表示一天的数据, 则有数组:

$Data\_Date = \left\{ \begin{array}{l} data\_dat{e_{1,1}}, \cdots ,data\_dat{e_{1,144}},\\ \quad\quad \cdots ,data\_dat{e_{x,y}}, \cdots ,\\ data\_dat{e_{144,1}}, \cdots ,data\_dat{e_{1,144}}, \end{array} \right\}$

其中, 下标 $x$ 和 $y$ 分别为每个数据点的对应的单元格横纵坐标. 研究区域的横纵坐标一旦确定, 从时间维度看, 每个区域的数据 $Data\_Dat{e_{x,y,n}}$ 也是一个数组, 存储每天该区域的OLR值, 则有 $Data\_Dat{e_{x,y,n}} = Z \cup {{\textit{z}}_i}$ , 其中 $Z = \{ {{\textit{z}}_1}, \cdots ,{{\textit{z}}_{i - 1}}\} $ 是同一个单元格中, 随时间变化的连续数据点, 代表该区域的历史数据, 而 ${{\textit{z}}_i}$ 表示该区域的当前数据点. 数据子集 $Data\_Dat{e_{x,y,n}} = Z \cup {{\textit{z}}_i}$ 可以视为一个连续变化的时间序列. 下面详细描述具体的算法流程.

在计算相对于背景数据的时间维度变化增量时, 先以地震中心所在单元格为中心, 获取 $5 \times 5$ 个单元格的震前历年数据, 根据式(1)计算地震当年的OLR序列相对于背景数据的增量值 ${\nabla _{i,t}}$ , 获取地震发生年的时间维度上的增量序列, 称为背景增量序列.

${\nabla _{i,num}} = {{\textit{z}}_{i,num}} - \frac{{\displaystyle\sum\limits_{num = 1}^{num - 1} {{{\textit{z}}_{i,num}}} }}{{num - 1}}$

(1)

其中, $num$ 表示从可获取数据的起始年份到地震发生年的总年数.

由于地震发生时, 其周边环境也可能相互影响, 出现异常现象. 故本文将空间因素考虑在内, 研究地震中心区域相对于周围区域的变化情况. 根据式(2), 处理背景增量序列, 计算各个数据点的空间增量值 $\nabla ({\nabla _{x,y,i}})$ .

$\nabla ({\nabla _{x,y,i}}) = \frac{{{\nabla _{x,y,i}} \times 4 - ({\nabla _{x,y - 1,i}} + {\nabla _{x,y + 1,i}} + {\nabla _{x - 1,y,i}} + {\nabla _{x + 1,y,i}})}}{4}$

(2)

得到具有时空相关性的新增量序列, 记为 ${T_n}$ , 则有 ${T_n} = T \cup {t_i}$ , 其中T表示历史数据点, ${t_i}$ 表示当前数据点.

在提取随机漫步概率特征时, 设置一个滑动窗口, 窗口大小记为 $ws$ , 判断当前数据点的值是否大于或等于前一个数据点的值, 若是, 则当前监测数据点的值定义为向右走, 否则, 当前监测数据点的值定义为向左走. 并根据随机漫步概率式(3)计算出随机漫步概率特征序列.

$rw(ws) = {q^{\frac{{ws + i}}{2}}}{(1 - q)^{\frac{{ws - i}}{2}}} \times C\left(ws,\frac{{ws + i}}{2}\right)$

(3)

此处, $i \in [ - ws,ws]$ , $q = 0.5$ , 而 $C\left(ws,\dfrac{{ws + i}}{2}\right)$ 是数学中的组合公式.

计算异常度时, 采用K-means聚类方法, 计算各个点到聚类中心的距离, 记为 ${s_i}$ . 再根据式(4)统计前i个点的异常概率, 记为 ${\hat p_i}$ .

${\hat p_i}({Z_n},\theta ) = \dfrac{{\# \{ j\left| {{s_j} > {s_i}} \right.\} + \theta \cdot \# \{ j\left| {{s_j} = {s_i}} \right.\} }}{i}$

(4)

此处, ${\rm{\# }}\left\{ {} \right\}$ 是一个用于统计满足条件的数据点的下标个数. 如 ${\rm{\# }}\left\{ {j|{s_j} > {s_i}} \right\}$ 统计满足 ${s_j} > {s_i}$ ( $j = 1,2, \cdots ,i,\;$ $\;i = 1,2, \cdots, n$ )的 $j$ 的个数, ${\theta _{}}$ 是[0, 1]区间的一个随机值.

最后结合Martingale理论, 并且为了减小滑动窗口大小设置不当带来的影响, 设置滑动窗口值取值范围, 求多组滑动窗口计算得到的均值作为最终异常指标. 即根据式(5)计算异常值, 记为 $CD$ .

$CD_n^{(\varepsilon )} = \dfrac{{\displaystyle\sum\limits_{k = {\rm{1}}}^{{{winend - winbegin}}} {\prod\limits_{i = 1}^n {\left( {\varepsilon \hat p_i^{\varepsilon - 1}} \right)} } }}{{{{winend - winbegin}} + 1}}$

(5)

其中, $\varepsilon $ 是[0, 1]区间的一个任意值. $winbegin$ 和 $winend$ 是窗口的起止取值, 本文设置滑动窗口取值范围为[30, 45]. 同时预设一个阈值 $h$ , 将各个数据点的异常值 $CD$ 与 $h$ 进行比较. 若从某点开始, $CD$ 值突然显著增大并超过 $h$ , 则可视为出现异常. 此时, 需标记异常, 并重启算法, 重新计算该点之后的数据点异常值. 即算法重启规则如表达式(6)所示.

$C{D_n} \ge h$

(6)

综上所述, 震前异常时空分析法的算法流程图如图1所示.

3 结果与分析

实验时, 根据文献[4], 将阈值 ${{h}}$ 设置为1000. 参考文献[18], 取 $\varepsilon = 0.82$ . 研究的数据时间周期为一年, 研究数据从发生地震前一年的9月1日到地震发生当年的8月31日. 图2中的竖线表示地震发生时间.

3.1 汶川地震

2008年5月12日在中国四川省的汶川县发生了8级大地震. 地震中心所在位置的经纬度是(31.0°N, 103.3°E). 故地震中心所在单元格的经纬度为(28.75°N–31.25°N, 102.5°E–105°E), 实验研究区域为(26.25°N–33.75°N, 100°E–107.5°E), 总共9个单元格, 地震中心位于第5单元格中. 图2显示了9格单元格中的CD值变化情况. 从中可以看出, 在地震发生前, 第2, 3, 4单元格出现明显的波峰. 并且这些波峰都出现在汶川地震发生前3个月内. 这可能是由于地震前, 大量能量积累引起的. 第1, 7和9单元格在震前未发现明显变化, 却在震后出现了明显的波峰. 这可能是由于地震后大量能量释放导致的. 而第6和8单元格在地震前后都没有明显变化.

从图3可以看到, 第5单元格, 即地震中心所在单元格也在地震前一个月也有小波峰, 但是显然第5单元格的异常趋势没有第1, 2, 3, 4, 7和9单元格的异常明显.

综合图2和图3, 可以看出, 地震前后, 震中区域OLR值会发生显著变化, 并且地震中心周边区域也会出现较大异常现象. 而且较大的异常现象可能并不是出现在震中所在区域, 而是出现在其周边区域. OLR的变化在时间和空间上都有很大的关联性.

3.2 日本地震

2011年3月11日在日本发生了9级大地震, 给日本东北部造成了毁灭性破坏. 震中经纬度是(38.1°N, 142.6°E). 故地震中心所在单元格的经纬度为(36.25°N–38.75°N, 142.5°E–145°E), 实验研究区域为(33.75°N–41.25°N, 140°E–147.5°E), 地震中心位于第5单元格中. 从图4可以看出, 9个单元格除了第1和3单元格外, 都出现了较明显的波峰. 其中第2, 5, 7, 8和9都在地震发生前3个月内出现明显异常. 第4单元格, 虽然波峰出现在地震后, 但在地震发生前两个月内CD值也有明显波动的趋势, 并且3月11日, 即地震当天, 至3月15日的CD值均大于100. 第6单元格, 在地震后一个月内出现了明显的波峰. 从图5可以看出, 震中所在区域, 即第5单元格, 在地震发生前一个月内急剧增大, 并达到峰值. 在2月10日出现最大值, 为444.4.

3.3 芦山地震

为了探索异常大小与地震震级的关系, 本文也分析了2013年4月20日的芦山7级地震. 地震中心所在位置的经纬度是(30.3°N, 103.0°E). 故地震中心所在单元格的经纬度为(28.75°N–31.25°N, 102.5°E–105°E), 与汶川地震的震中位于同个单元格内. 研究区域与汶川地震相同. 从图6可以看出, 第1, 2, 3, 5和6单元格都在芦山地震发生前3个月内出现明显的波峰, 这与前面两个震例发现的时间规律吻合. 而第8单元格的异常出现时间更早, 在地震发生前5个月已经出现明显的异常波峰.

而从图7可以看出, 芦山地震发生前4个月已经出现明显, 并在地震发生前一个月内CD值再次急剧增大, 并在震前半个月达到最大值. 第4, 7和9单元格震前均无明显变化趋势.

从图2, 图4和图6整体分析, 可以发现, 芦山地震的震级明显小于汶川地震和日本3.11地震, 但震中所在区域的异常值却并不比汶川地震和日本地震的异常小. 这也可以看出, 并非震级越大, 异常值就越大. 异常与震级并不是简单的正比关系, 还与其他因素有关.

此外, 结合图1, 图3和图5, 从空间角度看, 发现3个震例的正北方向, 即第2单元格所在方向上, 都在震前3个月出现了明显的波峰. 其他震例的异常现象是否也具有类似的方位规律, 还有待验证.

3.4 对比实验

吴立新等^[19]采用概率统计学中的 $\mu + 2\sigma $ 作为地震异常的衡量指标, 其中 $\mu $ 表示当地震前多年的OLR均值, $\sigma $ 表示方差. 他们分析以震中区域的地震前后3个月的OLR序列, 发现很多震例的地震前有些点的值会大于 $\mu + 2\sigma $ . 为了验证本文所提出的震前异常时空分析法的先进性, 本文采用与他们相同的方法及相似的实验设置分析本文3个震例的OLR震前异常. 受篇幅影响, 仅以日本地震的结果为例作一说明. 我们分析了第5单元格2011年1月1日到4月1日共3个月的OLR序列. 实验结果如图8所示. 震前仅有一个点的值大于 $\mu + 2\sigma $ , 视为异常点. 汶川和芦山两个震例的结果中震前也只能看到分散的几个异常点出现, 但是这些零星分散的异常点不易被发现和注意. 相对而言, 本文所提出的时空分析法也同样可以检测到震前异常值, 并且异常持续时间较长, 可以很直观地反映出了异常值的变化趋势. 因此, 我们的方法可以更直观有效的反映变化过程, 呈现当前数据点与它之前的数据之间的关系.

3.5 参数讨论

为了比较不同震例, 不同年份, 不同单元格的实验结果, 我们在实验中设置了相同的参数初值. 本文中设 $h = 1000$ , $\varepsilon = 0.82$ , $ws \in [30,45]$ . 为了说明本文所提出的方法对参数初值敏感性, 以及参数初值对实验结果的影响, 本文设计了对比实验, 讨论了3个参数的初值对3个震例的结果影响, 受篇幅影响, 以图9–图11的日本震例结果为例, 做一说明. 首先, 固定 $\varepsilon = 0.82$ , $ws \in [30,45]$ , 分别令 $h = 250$ , $h = 500$ , $h = 750$ , $h = 1000$ , $h = 1500$ 和 $h = 2000$ . 得到的结果如图9所示, 可以看出, $h$ 的初值大小仅改变波峰极大值出现的时间以及异常值的幅度, 如果 $h$ 设置太小, 则会较早出现波峰极大值, 如果 $h$ 设置过大, 则会导致极大值幅度很大, 但是并不会改变异常变化的趋势. 本文设 $h = 1000$ , 大小较为适中, 能较好地反映异常变化趋势. 然后, 固定, $h = 1000$ , $ws \in [30,45]$ , 分别令 $\varepsilon = 0.1$ , $\varepsilon = 0.2$ , $\varepsilon = 0.4$ , $\varepsilon = 0.6$ , $\varepsilon = 0.8$ 和 $\varepsilon = 0.9$ , 得到的结果如图10所示, 可以看出, $\varepsilon $ 的初值大小对异常值的变化幅度有影响, 若 $\varepsilon $ 太小, 可能检测不到异常, 而从 $\varepsilon = 0.6$ 开始, 异常变化明显; 当 $\varepsilon = 0.9$ 时, 异常值有所减小, 但是整体的变化趋势并没有太大影响, 不同 $\varepsilon $ 得到的波峰极大值出现的时间很接近. 本文根据适中原则以及实验结果, 设置 $\varepsilon = 0.82$ , 可以检测到较明显的变化趋势. 最后固定 $h = 1000$ , $\varepsilon = 0.82$ , 分别讨论 $ws = 25$ , $ws = 30$ , $ws = 35$ , $ws = 40$ , $ws = 45$ 和 $ws = 60$ . 得到的结果如图11所示. 从中可以看出, 不同的滑动窗口也会影响波峰出现的时间以及变化幅度. $ws$ 小, CD值的变化幅度也较小, 变化趋势不明显, 随着 $ws$ 的增大, CD值的变化幅度也会更大, 趋势更明显. 但是如果滑动窗口太大, 会导致窗口起始值前面很多数据被浪费. 一般来说, 选取60天以内的窗口尺度比较合适. 本文为了减少不同窗口初值对结果的影响, 分别设置窗口值为30到45, 计算异常值, 再取多窗口尺度的平均异常值作为最终的异常指标.

综合图9–图11可以发现, 不同的参数初值可能对实验结果的变化幅度以及异常出现的时间先后有细微影响, 但是不会改变异常的变化过程, 整体变化趋势是类似的.

4 结论

本文提出一种基于OLR数据的震前异常分析方法-时空分析法, 并结合汶川地震, 日本3.11地震和芦山地震3个震例阐述该算法的适用性. 从实验结果可以看出, 时空分析法由于既考虑了不同年份, 同一个地区背景数据的影响, 又考虑了同一时间, 周边区域的影响, 还考虑了同个季节内, 历史数据对当天数据的影响. 相对于只考虑其中一种因素影响的研究来说, 本文提出的时空分析法考虑得更加全面, 分析更合理.

通过分析以震中区域为中心的9个单元格, 可以发现几个规律: 1) 地震发生前后, 地震周边区域的OLR也会出现明显的变化; 2) 震中区域出现的异常不一定会比周边区域的异常明显; 3) 并且一般在震前3个月就已经出现明显的异常变化趋势. 4) 震级大小与异常大小并不是简单的正比关系, 还与其他因素有关; 5) 从方位角度分析, 发现3个震例的正北方向, 即第2单元格所在方向, 都在震前3个月出现了明显的波峰. 方位角度的这个规律是否具有普遍性, 还需要后续进一步验证. 受限于当前研究的数据量和研究水平, 目前只得到了一些初步的规律, 要得到更精确的, 更具有普遍性的结论还需要收集更多的震例数据, 进行更完善的实验, 进一步探索与震前异常相关的参数, 及各参数之间的关联性.