利用均线差值振幅波动水平分析汶川、雅安芦山、九寨沟地震前热红外异常＊

廖洪月 董 娜 王 刚

1) 西安市地震局，西安 710007

2) 西安市地震监测中心，西安 710007

引言

在地震监测预报中，卫星热红外观测资料与其他地震前兆资料相比，其突出的优点是覆盖面广，平面分辨率高，时间连续好，且不受地理环境限制[1]。在多山的四川省，于2008年、2013年和2017年分别发生了汶川MS8.0地震（31°N，103.4°E）、雅安芦山MS7.0地震（30.3°N，103°E）和九寨沟MS7.0地震（33.2°N，103.8°E）。地震发生后，国内学者相继进行了与热红外相关的研究并取得了一定进展。目前，热红外异常提取方法主要包括4大类：①目视解译；②基于差值分析的异常提取算法；③基于信号分析的异常提取算法；④基于背景场分析的异常提取算法[2]。本文主要运用第②类方法，较为成功的提取出了汶川、雅安芦山、九寨沟地震前热红外异常的时空序列。

当前学术界主要研究方法为小波功率谱分析法、距平法等，主要侧重于短时间内绝对亮温升温异常的研究。本文的差值振幅分析法则侧重于中短期相对升温异常的研究，且在不需要剔除错误数据和去云干扰影响等人工干预手段的情况下，提取的震前震后异常时空特征较为明显，异常时间和异常区域的研究结果与有关文献有相同之处（在后面具体震例分析中有描述）。

1 数据选取

本文热红外数据是由四川省地震局提供的TERRA卫星MODIS遥感传感器亮温数据。数据时间范围为2004年1月1日—2020年12月31日，数据空间区域范围为（30°N—45°N，100°E—130°E），数据分辨率为1 km，由450万个像元组成，单个像元每天一组数据，数据覆盖中国大陆西南、西北、华北、华东等部分地区。

在数据分析过程中，将分辨率为1 km数据分别反演成10 km、20 km、30 km、50 km和100 km分辨率数据后进行相关处理分析。对比发现，数据分辨率并非越高越好，虽然不同分辨率数据的时间曲线、空间图形态总体相似，但是分辨率较低时，空间图相关指数等值线平滑度更佳，且分辨率越低后续数据处理效率越高。但分辨率过低，存在影响异常空间定位精度和异常峰值时间滞后等缺点。本文以100 km分辨率热红外数据作为研究的基础数据。

2 数据处理及研究方法

前人利用不同的卫星遥感热红外数据，采用不同方法研究大震前热红外异常的结果显示，获取的异常信息几乎都是升温异常[3]。本文研究也是以升温异常为研究对象，不同之处在于侧重研究升温过程中均线差值振幅的变化水平异常，而非绝对的亮温升温异常。

2.1 数据处理

本文的数据处理流程如下：

（1）数据反演：将空间分辨率为1 km数据分别反演为10 km、20 km、30 km、50 km和100 km数据。以100 km分辨率为例，即将100 km×100 km区域内10000组1 km分辨率网格数据进行均值化处理，转化为一组数据，将450万个1 km网格单元转换为450个100 km网格单元。

（2）滑动均线计算：以365天作为年周期对反演后的分辨率为100 km亮温数据进行滑动均线计算。将数据区域内450个网格单元每天的亮温数据转化为年均线数据。

（3）差值分析：对每个单元网格年均线差值的振幅变化特征进行分析。

（4）产出异常区域单元网格时间曲线、区域空间异常叠加图和剖面图等。

实现流程（3）和流程（4）相关功能的计算公式为：

式（1）中，Ai和Ai−t分 别为第i天、第（i−t）天的亮温年均线值。在研究中发现，数据分辨率越低，代表网格区域空间范围越大，能观察到区域内的亮温异常现象就可能越多，观察到异常的时间也可能越长。当用均线分析异常时，均线周期应大于热红外异常发展周期，这样的均线值就可能包含整个异常发展周期的信息。为避免均线值受到季节变化的影响，本文特将均线周期设定为365天，默认为1年，在理想情况下，年均线任何时间点的均线值相等。当出现亮温升温异常时，年均线值会伴随亮温升温发展过程不断增大。△Bi是红外亮温均线值Ai的 差值，t为差值分析的时间间隔，可根据需要设置。文献[4-6]表明，热红外异常通常只持续10天到几十天之间。当前均值与t天前均线值进行差值对比时，假设t天前处于正常状态。t值不宜过小，当t值过小时，则可能出现异常值与异常值间的差值比较，因异常对冲而导致提取的异常减弱或异常被屏蔽现象；t值也不宜过大，t值过大时可能导致本次周期异常与上一个周期异常进行差值对比，也可能出现异常屏蔽现象。本研究经过反复测试，认为t值约为热红外异常的发展周期为最佳，但判断异常发展周期非常复杂。不同的t值，产生的差值系列不一致，最后提取的异常结果也有所不同（表1）。本文研究上述3次地震时，t值统一设定为20天。式（2）中，△Bmin为数据起始时间到当前时间段内的最小差值，可能随时间变化而变化，是一个变量，△Ti是当前亮温均线差值与最小均线差值的差值，也就是当前差值的振幅变化值，振幅变化值越小，表示热红外活跃性越低，振幅变化值越大，表示活跃性越高。式（3）中，是当前的平均值，n为数据起始时间到当前的总天数，随着时间变化而变化，也是一个变量，值高意味着对应时段内热红外活跃性高。本研究认为，n的最大值也就是选择的数据时段不宜过长，要避免数据时段内存在多个热红外异常时段，包括已经对应发生地震的异常和其他离当前时间较长的热红外异常，数据时段内只存在一个异常时段最为理想。在本文中，选择连续数据时长一般不超过2年，不小于1年。式（4）中，用于计算均线差值振幅变化值与平均变化值的变化比，比值为增大倍数，将其作为亮温异常参数，其基准值为零，大于零时表示处于升温异常状态，小于零时表示升温活跃性水平下降，其活跃性低于历史平均水平，也可理解为处于降温异常状态。在研究中，为突显异常观测效果，在不改变异常参数正负属性条件下，对增大倍数进行3次幂运算放大，当存在异常时，放大后的幂值在时间曲线中有一目了然的视觉效果。空间剖面图等值线Yi则为原始值。

表1 地震热红外异常特征相似点对比Table 1 Comparison of the thermal infrared anomaly similarity

2.2 异常分析原理

异常分析原理和异常判断标准：①时间异常分析。理论上均线差值应该是一个稳定值，包括两种情形：当均线为水平线时，差值应该为0；当均线为一定斜率的直线时，差值为一个稳定值，此时差值振幅变化值应为0。当均线不稳定时，差值数据系列的振幅也发生变化。振幅变化越大，出现异常的可能性越大。本文设定的异常判别条件是：当振幅变化大于平均值的2倍时，则判断为增温异常，值越大，表示时间异常越高。②空间异常分析。当空间平面某区域的异常值高于2且在整个空间中处于相对高位时，认为该区域存在空间异常。当时间和空间异常同时存在且相互对应时，将该现象判定为震前异常。

3 结果分析

3.1 汶川8.0级地震

以汶川地震为例，实际工作中快速异常分析分为3步：①通过时空扫描绘制2006年8月1日—2008年5月12日即汶川地震前这段时间的异常叠加图，通过异常叠加图找出异常区域；②提取异常核心区时间曲线，分析该区域震前热红外活动特征，判断异常发生时间等；③以5天为步长，绘制连续的异常平面图，分析震前与震后热红外活动的空间变化特征。

本文的异常叠加图是指一段时间内产生的异常在同一平面上的投影，它能显示异常发生的位置但不能显示异常发生的时间。本文3个震例分析时选取空间的范围为（30°N—45°N，100°E—115°E）。

由2006年8月1日—2008年5月12日时间段内时间扫描异常叠加图发现，以（37°N，102°E）为核心的较大区域出现异常（图1）。异常核心区网格单元的时间曲线显示异常峰值发生于3月18日，距离汶川地震55天，异常放大值约为55，折算值≈3.8（图2），即最大峰值倍数为3.8倍，表示当天差值振幅变化水平高于历史平均水平3.8倍。从异常叠加图可以看出，异常区域在震中的北部偏西，地理上异常主体区域位于柴达木块体，该结果与路茜等[3]的研究有相同之处，异常发生时段的结果与Xie等[4]的研究有相同之处。

图1 2006年8月1日—2008年5月12日异常叠加图（红圆为震中）Fig.1 Abnormal stacking diagram from August 1，2006 to May 12，2008（the red circle is the epicenter）

图2 异常区均线差值振幅随时间变化关系（纵轴为异常放大系数）Fig.2 The relation between the amplitude of mean difference in the abnormal area and time（the vertical axis is the abnormal amplification index）

对步长为5天的连续空间剖面分析发现：从2008年1月下旬开始，汶川地震中北部区域出现大范围的低值异常，表明该区域热红外活跃性降低，该状态持续约1个多月。从3月初开始，震中北部区域异常值开始缓慢上升，3月中旬高值异常范围达到最大，之后开始减弱。至汶川地震时，异常值较低，震后约20天之后的整个6月份，整个数据区域内大范围呈现低值异常状态，说明这个时段热红外活跃性水平非常弱（图3）。

图3 汶川地震前后热红外活跃性水平时空变化过程Fig.3 Temporal and spatial variation of thermal infrared activity level before and after Wenchuan earthquake

3.2 雅安芦山7.0级地震

雅安芦山地震的分析方法步骤与汶川地震相同。2012年1月1月—2013年4月20日异常叠加图显示存在两个高异常区域，一个异常位于（40°N，112°E）区域（图4右上），异常区呈短条带状。该区域时间曲线显示于2012年5月出现过一次高值异常，经查询发现同年5月9日和6月18日在山西与陕西交界处各发生一次4.3级地震（数据来源于USGS），本文不作详细描述。另一个异常位于雅安芦山震中的西北方向（图4左下），由于受热红外原始数据空间范围局限，只能显示部分异常区域而不能看到全貌。

图4 2012年1月1月—2013年4月20日异常叠加图（红圆为震中）Fig.4 Abnormal stacking diagram from January 1，2012 to April 20，2013（the red circle is the epicenter）

异常核心区域时间曲线显示，异常峰值发生于3月16日，距雅安芦山地震35天，最大幅度变化高于历史平均水平3.1倍（图5）。

图5 异常区均线差值振幅随时间变化关系（纵轴为异常放大指数）Fig.5 The relation between the amplitude of mean difference in the abnormal area and time（the vertical axis is the abnormal amplification index）

步长为5天的连续异常图（图6）显示：从2013年2月上旬至2月下旬，震中的偏西区域，热红外活跃性水平低值异常持续了20多天，之后活跃性强劲上升，3月中旬达到峰值后迅速下降。至4月20日地震时，震中及附近区域无明显异常。震后约20天之后，震中附近区域热红外活跃性水平开始下降，低活跃性持续了约20天，该异常现象与汶川地震后低值异常相似。本次地震前活跃性高值异常区域的研究结果与魏乐军等[5]的研究结果相近。

图6 雅安芦山地震前后热红外活跃性水平时空变化过程Fig.6 Temporal and spatial variation of thermal infrared activity level before and after Lushan earthquake in Yaan

3.3 九寨沟7.0级地震

分析方法与汶川地震相同。2016年1月1日—2017年8月8日时空异常叠加图显示，以（35°N，101°E）为核心的较大区域存在高值异常（图7），以（35°N，101°E）为中心的异常区域时间曲线显示，分别在1月、5月和7月出现过3次较大异常，最大异常峰值发生于7月17日，距离九寨沟地震22天，当日最大振幅变化高于历史平均水平3.1倍（图8）。

图7 2016年1月1日—2017年8月8日异常叠加图（红圆为震中）Fig.7 Anomalous stacking diagram from January 1，2016 to August 8，2017（the red circle is the epicenter）

图8 异常区均线差值振幅随时间变化关系（纵轴为异常放大系数）Fig.8 The relation between the amplitude of mean difference in the abnormal area and time（the vertical axis is the abnormal amplification index）

异常剖面图变化过程显示，从2017年4月中旬至5月初，震中西部区域有一个明显的异常指数下降过程，说明该区域热红外活动性水平在这段时间呈下降状态，之后活动性呈现宽幅的振荡变化，分别于5月20日和7月20日在该区域产生高值异常后迅速下降，第2次异常强度和异常区域面积均大于第1次，8月8日发震后震中南部区域活跃性指数有下降现象（图9）。

图9 九寨沟地震前后热红外活跃性水平时空变化过程Fig.9 Temporal and spatial variation of thermal infrared activity level before and after Jiuzhaigou earthquake

本次地震的震前异常时间与异常区域研究结果与戴勇等[6]、杨星等[7]研究结果有相同之处，但异常区域研究结果与张丽峰等[8]的研究结果不一致。

3.4 亮温异常的时空变化特征

不同的时间间隔，均线差值振幅变化结果有所不同。汶川地震、雅安芦山地震、九寨沟地震前与震后相关异常特征简要对比如表1。

3.5 低分辨率反演与年均线差值分析的潜在优点

在处理亮温数据时，通常需要人工干预：剔除错值、云干扰值及不符合黑体辐射公式的高值，计算剩余数据的均值，从而得到扣除部分云影响的亮温日值[9]，本文数据处理过程中省略人工剔除错误数据和去云干扰工作等环节，也达到几乎相同异常提取的效果，也就是说，本文方法具有一定的容错功能和自动去云干扰影响功能，其原理如下：

容错原理：将高分辨率由1 km×1 km数据反演为100 km×100 km数据后，单个错误值的影响稀释到原来的1/10000，将分辨率为100 km×100 km网格数据进行年均线转换后，可再将影响再稀释为现有的1/365，总影响量降为最初的1/3650000，此时错误数据负面影响非常有限，说明低分辨率反演和年均值转换处理具有非常理想的容错功能。

自动去云干扰影响原理：我们可以简单地认为，任意一组年均值由受云干扰影响的亮温值和正常陆地亮温值组成，同时也可以简单地认为，每组年均值受云干扰影响的量值是相同的。这样，通过差值计算，直接将双方均拥有的云干扰影响去除。本文中时间相差20天的年均线数据包含了相差无多的云干扰影响，在差值计算过程中云干扰影响被自动全部对冲或绝大部分被对冲掉，因此，无需人工处理。本研究认为，当均线周期较短时，计算差值、时间间隔中可能存在季节差，尤其是存在阴雨季节差别时，不能忽略云干扰影响，仍有必要做去云干扰工作。但当均线周期为年或年的整数倍时，均线中每组数据中包含相同的季节，滑动差值计算时不存在季节差，因此，可以忽略云干扰影响。

由于不需人工干预，且不同的数据处理流程均有统一的算法，因而数据处理的全流程均可实现计算机全自动分析完成。这样，不但提高了数据处理效率，同时也降低了数据处理和分析的技术难度，从技术角度为研究成果转化和应用提供了便利条件。

4 结论与讨论

4.1 结论

通过以上的系统分析，得出主要结论如下：

（1）本文热红外异常参数Yi的实质是变化率，其基准值为零，其异常活跃性水平是针对升温异常而言，其数值的大小仅代表升温异常活跃性水平的高低，其内在的意义如下：异常值越趋近于零，表示热红外活动越稳定；当异常值大于零时，异常值越高表示热红外升温异常活跃性水平越强；当异常值小于零时，异常值越低则表示热红外降温异常活跃性水平越强。

（2）汶川、雅安芦山和九寨沟地震前均出现热红外亮温均线差值振幅变化正异常，也就是活跃性上升异常，且地震均在峰值后较长时间发生，该研究结果具有地震短临预测的意义。

（3）上述3次地震中，汶川地震级别最高，其震前异常值大于2倍的异常区域也最大。3次地震震中均在异常区域外围而不在异常优势区域，震中均位于异常区东南方向。

（4）上述3次地震前伴随有活跃性低值异常和活跃性高值异常。先有低活跃性异常，再演变为高活跃性异常，低值时低至历史平均水平50%以下，且低值状态存续时间较长，高值时高于平均水平3倍以上。震前热红外存在两种性质相反的异常现象，在文献中较为罕见。

（5）本方法在提取异常过程中虽然缺少了剔除错误数据和去云干扰影响两个看似不可或缺的环节，但仍能提取到与前人相近或相同的研究结果，说明本方法可行性较高。

（6）通过低分辨率反演原始热红外数据和对反演后日均值数据实行年周期均值转换，有效地稀释了原始数据中可能存在的错误数据带来的负面影响，使本方法具有一定的容错功能。

（7）本方法巧妙利用年周期均线差值对冲功能，从而自动地实现了最大限度去云干扰影响。换句话说，通过均线年周期的设定，使本方法具有抗云干扰功能。

（8）本研究方法是对全时段、大区域的卫星原始数据进行系统的整体的分析，不对原始数据做任何的删改，确保在百分之百真实的原始数据基础上提取异常，所有数据参与处理的流程统一，分析方法统一，每组原始数据均同权参与全流程计算，最后得出的数据结果不含任何主观因素。由于实现了高速高效的自动化数据处理流程，本文认为该方法适合大规模、大范围数据甚至全球热红外数据同步或准同步处理，可用于全国性大区域热红外实时异常动态监测。

4.2 讨论

在研究过程中，同时也发现了另一些现象，值得进一步讨论和研究：

（1）与前人通常侧重于绝对升温异常不同，本文是以相对异常作为研究方向。亮温升温绝对高值异常与相对高值异常并不一定同步，绝对亮温升温高值异常时并不一定产生相对亮温高值异常，同理，相对高值异常也不一定会同时出现绝对升温高值异常，因此，本文关于异常区域和异常时间的研究结论与其他学者研究结论存在出入属于正常现象。

（2）汶川、雅安芦山、九寨沟3个地震前均存在相同的异常发展过程，先有活跃性水平大幅下降，后又转为活跃性大幅上升，高值异常峰值过后发震。震前两种性质相反的异常现象是否为大地震前共有现象以及其产生的原因，值得深入研究。

（3）对冲的作用与副作用。本文差值分析时应用对冲特性去云干扰影响，从而省略了去云干扰数据处理人工干预工作，但差值间隔时间设置不当可能会对冲掉部分真实异常，影响最终的异常提取效果。因此，差值间隔时间参数设置是一个值得深入研究的问题。

（4）3次地震震前异常值在3.1—3.8倍之间，以本文的2倍可作为判断异常的参考指标是否恰当，需更多震例验证和深入研究讨论。

（5）汶川地震前，距震中500—600 km的大面积异常区域是不是该震前兆值得讨论。本文认为，虽然异常区域距震中较远，但从图2的时间曲线分析，此区域在超过10年的时间内，仅在汶川地震前出现过一次这种罕见异常现象，因此，从时间和概率上对应分析，它是汶川地震前兆的可能性较大。根据马瑾院士[10]的研究，一个地震的发生可对周边不同构造区造成升温或降温影响，距离有远有近，500—600 km在该文献的震例中不算是远距离。因此，这种现象在空间上认为是汶川地震前兆也是可能的。

（6）本研究方法的缺点，在特定条件下可能出现伪异常。假设热红外历史数据非常平稳，则其年线差值振幅变化的历史平均值趋近于零，此时如果热红外出现小幅异常变化，则可能出现数值特别高的异常值，因此，要慎重对待高值异常的计算结果。

本文采用与众不同的热红外分析方法，发现了一些有趣的现象，这些现象是否可作为地震前兆还需进一步研究。由于本文所用到的数据范围较小，震例仅3例，虽然本文研究结果显示，震前异常较为明显且异常峰值出现在地震发生前，但由于3个震例的震中均不在异常区内，因此，本研究方法用于实际地震预测还有很远的距离，研究方法的科学性也有待更大范围的数据累积和更多震例的验证，希望获得更多学者的批判指正。

致谢

作为一名从事行政工作的业余研究者，感谢西安市地震局领导对本人和本研究提供宽松的研究环境与相关支持，同时感谢同事田晖先生为本文图片合成的辛苦付出。