引言

遥感技术因其快速获取大范围区域影像数据的优势，早在二十世纪八十年代就用于地震热异常的分析研究，并且取得了诸多成果（Tramutoli et al，2001；Choudhury et al，2006；Yao，Qiang，2012），震前热异常也在诸多震例中得到验证（吕月琳等，2009）。强祖基等（1990）通过多个震例探讨了临震卫星热红外异常是否与地表增温异常相关，其结果显示临震出现异常增温是普遍现象。吴立新等（2004）对岩石破裂的热红外辐射规律开展了试验研究，在探索地震热异常的形成机理方面取得了进展。

地震热异常具有强干扰、弱信号的特点，即与温度有关的气象信息容易受到多种复杂因素的干扰，而与地震有关的热异常较为微弱，因此需要排除强烈干扰背景以有效地提取热异常信息（李继业等，2011）。吴立新等（2017）引入热红外背景场，以热异常定量指标为判断依据，排除气象干扰构建较为稳定的背景场，对卫星影像观测值与背景值的偏差程度予以分析，发展了多种方法用于热异常提取。近年来，差值法、涡度法、小波功率谱法等多种方法广泛应用于地震热异常的提取分析，取得了诸多成果（陈顺云等，2006；郭晓等，2010；张元生等，2010；康春丽等，2011；刘军等，2014；魏从信等，2018）。陈顺云等（2006）基于NOAA气象卫星的甚高分辨率扫描辐射计（Advanced Very High Resolution Radiometer，缩写为AVHRR）陆地数据集，运用小波理论将地表热红外辐射分解为三个频率成分，结合气象资料和地震活动资料进行综合分析，其结果表明低频成分的热异常更接近于断层活动所引起的热异常信息，低频成分的高温异常区与构造特征、地震活动密切相关，因此低频成分热异常可为断层或地震活动提供一定的指示信息。马未宇等（2006）关于美国国家环境预报中心（National Centers for Environmental Prediction，缩写为NCEP ）提供的再分析数据温度异常图像与天体引潮力附加构造应力的联合分析结果显示，西藏的三次地震均发生于天体引力潮附加构造应力变化达到相对稳定值的末端。康春丽等（2011）基于NOAA-18提供的分辨率为1°×1°网格化的夜间长波辐射（outgoing longwave radiation，缩写为OLR）数据，使用涡度算法分析了四次大地震震前的热异常，结果显示震前一到两天热辐射异常明显增强，震后逐渐消失。

卫星数据鲁棒分析技术（robust satellite techniques，缩写为RST）自Tramutoli （1998）提出以来，被用于大量的震例研究，相关结果显示多个震例前均有热红外异常现象出现（Filizzola et al，2004；Tramutoli et al，2013；Genzano et al，2015；宋冬梅等，2016）。Ma等（2012）基于欧洲第二代气象卫星（Meteosat Second Generation，缩写为MSG）红外波段数据运用RST算法验证并提取了意大利阿布鲁佐大区M_S6.3地震前的热异常。宋冬梅等（2016）基于地表温度（land surface temperature，缩写为LST）数据使用RST算法再次对汶川地震进行分析，结果显示所提取的热异常与地质断裂更加吻合。Eleftheriou等（2015）也使用RST算法提取了2004—2013年希腊地区的长时序地震热异常，并对其特点予以分析。

上述研究用于热异常提取的遥感数据通常为夜间某一时刻或夜间均值数据，但是某时刻数据并不能准确地反映夜间热环境的整体情况，且均值数据也有可能导致热信息被湮没。为更加准确地判断热异常是否与地震相关，本文拟综合静止气象卫星的高时间分辨率特性，改进热异常与地震相关性的判别方法，运用RST算法提取夜间均值数据热异常后，进一步提取夜间多时刻热异常，以期通过分析多时刻热异常的时空分布来判断热异常是否可能与地震相关，进而为地震热异常判定提供参考。

1.   数据源选取及热异常提取方法

1.1   数据源

本文采用MSG搭载的旋转增强型可见光及红外成像光谱仪（Spining Enhanced Visible and Infrared Imager，缩写为SEVIRI） 10.8 μm波段数据。该型成像仪具有4个可见光—近红外通道和8个红外通道，其数据具有较高的时空分辨率，在空间上，中红外通道分辨率为3 km，比第一代的5 km有较大提高，且图像量化等级提高到10 bit （徐清，2013），极大地提升了数据质量。

1.2   去云

在热异常提取过程中，由于经常受到云层的干扰，可能无法提取出热异常。为去除遥感数据中云层带来的干扰，本文采用一种云检测方法（Cuomo et al，2004），即基于历年同期数据，通过设置迭代相应阈值，识别出云像元，具体公式如下：

This page contains the following errors:

Below is a rendering of the page up to the first error.

1.3   RST算法

RST算法是一种基于多时相数据的分析方法，其优势在于突出时空域上相对空间域的异常，最初被用于云检测、火山、森林火灾监测等方面（Tramutoli，1998；Genzano et al，2009；Tramutoli et al，2013），后经不断改进应用于地震热红外监测。该方法将地震引起的热红外异常看作信号，其它因素引起的异常看作噪声，引入地震热异常信噪比的概念（孙珂，2017）。热异常稳健估计（robust estimator of thermal infrared anomalies，缩写为RETIRA）值R （Tramutoli et al，2001）是以多年同期的遥感数据为基础，运用RST算法构建去除噪声的背景场，提取相对于背景场的偏离值，表达式如下：

式中，T （x，y，t）为（x，y）处的像元值，T （t）为影像区域像元均值，ΔT （x，y，t）为热异常待提取点位的像元值相对于区域均值的差值，μ_ΔT（x，y）为ΔT （x，y，t）多年相同位置、相同时刻的平均值，σ_ΔT（x，y）为对应的标准差。

由于在构建背景场所用到的数据存在干扰噪声，需进一步去除这些噪声。通过

对遥感影像进行处理，构建去除噪声的差值背景场。式（5）和（6）中，A（x，y，t）用来剔除构建背景场的异常值，μ_ΔT（x，y）为剔除背景噪声后的背景场（Zhang，Meng，2019）。构建稳定背景场的具体要求为：对于该点位的时间序列，与均值之差的绝对值大于k （k≥2）倍标准差的A值不可用于构建背景场。将A（x，y，t）赋值为0，迭代直到数值与均值之差均小于k倍中误差（Genzano et al，2009）。由式（6）计算相应的背景场μ_ΔT（x，y），进而得到R值。

2.   多时相地震热异常判定

人类活动及极端天气、大范围洪水、森林火灾等多种因素的影响均可能导致热异常，这使得所提取的地震热异常可能具有假阳性（Filizzola et al，2004），因此需通过筛选来减少非地震因素所造成的干扰。Eleftheriou等（2015）提出的热异常判定条件具体如下：

1） 选取相对偏离阈值为4，根据数据特性选取R≥4可有效检测出热异常；

2） 在时间分布上，第一次出现热异常后七日内应至少再出现一次，且两次热异常的分布位置相同或相邻；

3） 在空间分布上，热异常分布应当是聚集的，检测出的热异常在1°×1°的范围内应当至少覆盖150 km²的面积。

Eleftheriou等（2015）认为满足上述三条即可视为可能与地震相关的热异常。虽然基于影像均值或某时刻数据所获得的热异常提取结果有一定的代表性，可反映夜间一段时间内热辐射的基本情况，但是由于云层及其它干扰因素的不断变化，影像均值可能丢失部分热信息，难以反映夜间热异常的时空变化特征。因此，本文结合Eleftheriou等（2015）的三条热异常判定条件，增加判定条件4）：进一步提取分析夜间多时相数据的热异常，在排除云层干扰的情况下，由多时相数据所提取热异常的时空分布连续。这样，热异常在满足前三条热异常筛选条件的同时，仍需满足条件4），以排除短期天气因素所导致的干扰，满足该条件即可有效地排除异常时空分布强度不足、不连续及空间分布零散等不符合地震热异常的情况。

本文首先选取夜间0时、1时、2时的时刻数据，求均值后提取热异常，对于存在热异常的情况，再选取当日夜间20:30至次日05:00时段每隔30分钟数据，分析夜间18个时刻的遥感数据热异常，以探究其在夜间多个时段的时空演化特征，进而判定热异常与地震的相关性。

3.   热异常提取结果分析

3.1   数据均值热异常

本文基于RST算法探究意大利地区震例的震前热异常，使用2009—2018年SEVIRI遥感10.8 μm数据构建背景场。首先处理夜间0时、1时、2时的均值数据，提取三次热异常，分别为2015年5月6日和5月7日出现在意大利半岛南部和西西里岛北部的热异常、2016年1月9日和12日出现在意大利半岛中部的热异常以及2017年3月21日和23日出现在意大利半岛中部的热异常，其均值热异常如图1所示。根据研究区地形以及热异常序列的分布，本文认为未来地震可能发生于热异常中心附近。事实上，2015年5月6日和5月7日热异常出现后的第三天附近区域发生了M_S4.7地震，2016年1月9日和12日热异常出现后的第七天附近区域发生了M_S4.4地震，2017年3月21日和23日热异常出现后30日内附近区域未发生地震。

图  1  意大利三次夜间均值热异常分布

（a） 2015年5月6日和7日；（b） 2016年1月9日和12日；（c） 2017年3月21日和23日

Figure  1.  Three nightly mean thermal anomalies in Italy

（a） On May 6 and 7，2015；（b） On January 9 and 12，2016；（c） On March 21 and 23，2017

下载: 全尺寸图片 幻灯片

如图1a所示：2015年5月6日在意大利半岛南部和西西里岛北部出现了热异常，7日意大利南部半岛部分区域出现了云层，未能检测出热异常；5月6日热异常分布较7日的分布范围更广。根据热异常的判定条件，西西里岛北部等区域出现的单日热异常不属于有效热异常。由图1b可见：2016年1月9日和1月12日在意大利半岛中部出现了热异常，相关区域的热异常空间分布较为连续；1月12日在意大利半岛中北部和南部出现了热异常。这几处热异常不连续，因此不认为其为有效热异常序列。由图1c可见：2017年3月21日和23日意大利半岛中部区域出现了热异常，两日热异常序列空间分布较为连续，也存在部分区域出现零散热异常的情况。

3.2   多时相数据热异常提取

基于RST算法进一步提取夜间均值多时刻影像热异常，本文选取从当日20:30至次日05:00 （时间均为UTC）每隔30分钟的遥感影像作为基础数据，分析不同时刻的热异常分布特征，并根据文中的判定准则判断该热异常是否可能与地震相关。

通过分析意大利2015年5月6日和5月7日多时刻数据，得到热异常的时空分布序列，如图2所示。可见：5月6日各时刻均出现了热异常，本次热异常序列限于意大利半岛南部和西西里岛部分区域；5月6日凌晨前后多时刻热异常随着时间推移呈现逐渐增多的扩展趋势，其中意大利半岛南部和西西里岛北部的热异常均不断扩展，在震中附近的热异常分布呈现不断增多的趋势；5月7日凌晨前后多时刻热异常时空分布较为连续；由于云层的干扰，5月6日20时30分至5月7日0时0分这期间的7个时刻未出现显著热异常，然而随着时间推移，云层逐渐消散，相关区域热异常逐步显现、增多、扩展，因此认为此次热异常存在的不连续情况是由于云层干扰所致，根据文中的判定方法，判定此次热异常可能与地震相关。

图  2  意大利地区2015年5月6日和7日多时刻热异常时空分布

Figure  2.  Spatio-temporal distribution of multi-time thermal anomalies in Italy on May 6 and 7，2015

下载: 全尺寸图片 幻灯片

基于RST算法提取意大利2016年1月9日和12日的夜间多时相热异常，得到热异常的时空分布序列，如图3所示。可见：1月9日各时刻相关区域均出现了显著热异常，时空连续性较好，云层的变化情况清晰可见，相关区域附近的云量逐步增多，部分出现热异常的区域逐步被云层覆盖，导致无法连续提取出热异常；1月12日意大利中北部出现了较大区域的热异常，根据热异常的分布特性，此部分热异常与9日热异常的空间距离较远且不连续，故认为此部分热异常不属于有效热异常；1月9日和12日夜间意大利半岛中部均出现了连续热异常，时空分布具有较好的连续性，根据文中的判定方法，判定此次热异常可能与地震相关。

图  3  意大利地区2016年1月9日和12日多时刻热异常时空分布

Figure  3.  Spatio-temporal distribution of the multi-time thermal anomalies in Italy on 9 and 12 January 2016

下载: 全尺寸图片 幻灯片

通过分析意大利2017年3月21日和3月23日多时刻数据得到夜间多时相热异常的提取结果，如图4所示。3月20日和22日在21:30和22:00时刻均出现数据缺失，导致这两个时刻无法提取热异常，根据两个时刻前后热异常的提取结果，推断相关区域可能未出现热异常。从图4可见：3月21日凌晨前后多数时刻未出现热异常，仅01:00—04:00出现热异常，而后热异常再次消失，热异常未出现时刻未有云层出现，可排除存在云层的干扰，推测热异常可能由短期天气变化导致；3月23日凌晨前后各时刻数据多数时刻未提取出热异常，仅02:30—05:00少数时刻出现热异常，加上云层在相关区域分布较少，在一定程度上可排除云层对于热异常的干扰，但相关区域仍未出现连续热异常，根据文中的判定方法，认为此次热异常与地震无关。

图  4  意大利地区2017年3月21日和23日多时刻热异常时空分布

Figure  4.  Spatio-temporal distribution of thermal anomalies in Italy at multiple times on March 21 and 23，2017

下载: 全尺寸图片 幻灯片

3.3   多时相热异常分析与统计

根据三次夜间多时相热异常的时空分布可知，云层对于热异常提取的影响较大，云层的扩散与消失可能导致无法提取热异常或引起热异常时空分布变化较大。2015年5月7日和5月9日、2016年1月9和1月12日的热异常提取结果均受到云层干扰，云层的移动可能会引起热异常的空间分布产生变化。由2016年1月9日多时刻热异常时空序列的提取结果（图3），可明显看出1月9日1时起震中北部方向逐渐有云层出现，导致北部区域的热异常逐渐消失，而震中西北部地区逐渐出现了热异常；在2015年5月6日和5月9日、1月12日相关区域热异常（图2）通常因云层边缘的移动而出现变化，故认为云层在相关区域的时空变化对于热异常提取具有重要影响。

对夜间多时相云层面积占成像区域面积的百分比进行统计，结果如图5所示，可见：2015年5月6日凌晨前后云层缓慢增多，5月7日凌晨前后云层所占比例逐渐下降；2016年1月9日凌晨前后多时刻云层呈现逐步扩展增多的趋势，12日凌晨前后多时刻云层分布相对平稳；2017年3月21日和23日凌晨前后云层均呈现逐步增多，但基本保持稳定。

图  5  意大利三次多时刻云层及热异常面积统计

（a） 2015年5月6日和7日；（b） 2016年1月9日和12日；（c） 2017年3月21日和23日

Figure  5.  Statistics of three multi-temporal cloud cover and thermal anomalies in Italy

（a） On May 6 and May 7，2015；（b） On January 9 and January 12，2016；（c） On March 21 and 23，2017

下载: 全尺寸图片 幻灯片

通过统计意大利三次地震热异常的空间分布面积，可进一步了解夜间热异常的分布变化，探究其分布形态。如图5a所示，自5月6日凌晨起各时刻热异常基本呈现逐步增多趋势，其中5月6日3时热异常分布范围最大，面积达15 030 km²；5月7日凌晨前后多时刻热异常面积整体较前一日明显减小；5月6日20:30至24:00时未出现显著热异常。由图2可知，5月6日凌晨前后，因有云层出现无法提取热异常，之后云层逐渐消散，热异常面积逐渐增大。如图5b所示，2016年1月9日凌晨前后各时刻均出现显著热异常，多时刻热异常呈振荡增多趋势；1月8日23时热异常分布范围最大，面积达2 567 km²，随后分布范围逐步缩小，至1月9日1时开始逐步增大，02:30后热异常面积减小，并趋于稳定；1月12日凌晨前后多时刻数据均出现了显著热异常，热异常面积整体呈现先增大后逐渐减小的趋势，在1月12日00:30降至最小值，为486 km²，之后逐渐增大，在05:00达到最大值，为6 228 km²。如图5c所示，3月21日和23日凌晨前后多数时刻未出现热异常，两日在01:30前均未出现显著热异常，其中3月21日03:00热异常分布范围为当日最大，面积为558 km²；3月23日04:00面积达到当日最大值，为459 km²。

4.   讨论与结论

针对热异常与地震相关性的准确判断，本文对Eleftheriou等（2015）提出的热异常与地震相关性的判定条件予以改进，深入地探究基于夜间多时相数据的热异常与地震的相关性，并基于RST算法提取了意大利地区2015—2017年期间的三次热异常，分析了这些热异常的分布及其变化特征。

1） 根据影像均值或某一时刻的影像值判断热异常是否与地震相关可能存在误判情形，具有一定的局限性，本文通过改进热异常与地震相关性的判定条件，探究研究区夜间多时刻连续热环境的变化，可清晰地观察到热异常在夜间的时空变化特征。

2） 根据改进后的热异常与地震相关性判定条件，探究夜间多时相数据的热异常，结果显示：两次夜间多时相热异常序列较为连续，认为可能与地震相关；一次夜间多时相热异常在排除不存在云层干扰后，夜间多数时刻热异常不具有连续性，仅少数时刻存在热异常，认为该次热异常可能与地震无关。本文根据三次热异常后是否发生地震，验证了改进判定条件，提升了热异常与地震相关性判定的准确性。

3） 三次夜间多时刻热异常的时空变化表明，云层对于热异常提取的影响较大。热异常可能因云层的干扰而无法被提取，云层的扩散与消散都可能对热异常提取造成影响，同时云层分布比例与热异常分布面积并无显著相关性。

综上所述，本文基于夜间多时相数据探究了热异常与地震的相关性，进行了意大利地区三次热异常与地震的相关性分析，初步证实了该方法的有效性。为进一步探究夜间多时相数据提取的热异常与地震相关性判定方法的区域适用性，仍需大量震例作进一步研究。