基于STL和GESD算法的于田两次Ms7.3地震热异常时空分布特征识别研究

米日阿依·买土地， 玉素甫江·如素力， 喀迪阿依·阿力木， 艾尔肯·图尔孙

（1.新疆师范大学地理科学与旅游学院，乌鲁木齐 830054；2.新疆干旱区湖泊环境与资源重点实验室，乌鲁木齐 830054）

地震是地壳快速释放能量过程中造成的振动. 地震给人类所带来的损失很大. 为了尽可能地减少地震灾害造成的损失，相关领域的学者多年一直专注于地震前兆信息提取和地震预测研究［1-3］. 地震前发生的，所有与地震有关系的一切不正常现象一般称为地震前兆，分为宏观前兆和微观前兆. 常见的地震前兆现象有：地震活动异常，地震波速度变化，地下水不正常改变，重力异常，动物异常，地声，地光，地表温度（LST，Land Surface Temperate）和地表亮温（Brightness Temperature）异常等［4-5］. 其中，震前亮度温度，地表温度（LST）异常是利用热红外遥感监测地震前兆的主要参量.

最初对地球内部的热红外异常与地震活动的对应关系进行探讨的是苏联Gorny 等，随着空间技术的发展，采用热红外遥感数据监测震前LST异常成为可能［6-7］. 热红外遥感在地震监测以及地震预报中的应用越来越丰富. 地震所涉及的异常中存在很多复杂问题，怎样用有效方法，从地震热异常中获取地震相关的异常信息也是个难点. 最前期的热红外异常现象主要是通过目视判读发现的，现在跟着研究的需求，一个个又复杂又精细的研究思路，研究方法开始出现在地震热红外信息的提取研究当中. 当今，涉及遥感数据检测地震热异常的主要方法有背景场差值法［8］、小波分析法［9-11］、异常点比值时序法［12］、RST（Robust Satellite Techniques）等［13］. 目前为止，不同算法在不同地震活动异常信息监测中的作用和效果都不同，提高异常信息提取中的不确定性. 因此，热异常监测的算法及其应用还需深入研究.

然而，LST的变化过程不仅具有明显的季节性变化，还有年际趋势变化的特征，是由多种因素共同作用的结果，不能直接对其进行监测分析，而要先对其进行分解. 相关研究表明，时间序列分解法（STL，Seasonal-Trend Decomposition using LOESS）［14］可将时间序列分解为趋势项（trend）、季节项（seasonal）和残余项（remainder）.如果STL分解法能用于地震区地表温度的时空数据集，有效去除趋势和季节性变化的背景场信息，其对残余项用GESD检验（Generlized Extreme Studentized Deviate Test）进一步进行异常检测，取出地震热红外异常. 这对防震减灾有着重要的实际应用价值和现实指导意义.

本文拟利用基于MODIS LST 数据，以于田县2008年3月21日7.3级与2014年2月12日7.3级地震为例，采用STL分解法、GESD检验，以震中周围区异常面积时空变化作为主要参量，获取与此次地震相关的地震热红外异常演化特征.

1 研究区概况

于田县位于新疆维吾尔自治区西南部，东临民丰县，北邻塔克拉玛干沙漠与沙雅县接壤，西与策勒县相毗邻，南与西藏自治区改则县、日土县相接. 由南部山区、中部平原、北部沙漠区三大部分组成，南部山区海拔高度平均在4000 m以上.

2008 年3 月21 日新疆于田发生Ms7.3 地震. 震中（35.60°N，81.60°E），震源深度33 km. 地震发生后发生345 次余震. 2014 年2 月12 日新疆于田发生Ms7.3 地震，据中国地震台网中心测定，震中位于（36.1°N，82.5°E），震源深度为12 km. 地震的前1天曾发生Ms5.4地震. Ms7.3主震后余震活动频繁，最大余震为2月14日17时24分发生的Ms5.7地震.

野外地震地质调查表明，2008年于田地震震中位于西昆仑山中的阿什库勒盆地，距阿尔金断裂南缘约50 km. 处于西昆仑地震带与阿尔金地震带的交会区，即塔里木盆地的南部，靠近西昆仑山区的边缘部分.造成了全长约30 km，走向南北—北东东，同时兼有左旋走滑和正断分量的张剪地表破裂，最大左旋位移1.8 m，最大垂直位移约2.0 m，发震构造属于阿尔金断裂西南尾端的张性构造区. 美国地质调查局（USGS）与哈佛大学CMT震源机制解结果表明，此次地震为略带走滑的正断层性质，断层走向为近南北向. 而2014年于田地震也是发生在阿尔金断裂西南段，这次地震因发生在高山无人地区并未造成人员伤亡，其发震断裂为阿尔金断裂带西南段，阿尔金断裂带是亚洲大陆内部一条长达2000 km三叠纪以来长期左旋走滑的岩石圈断裂. 地貌上表现为醒目的线状展布特征，呈东西至北东东向，是青藏高原北边界断裂. 但在较短的时间内在同一地区发生2次7级以上地震，引起了地震学者的广泛关注［15-17］.

2 数据及其处理

MODIS也被称之为中分辨率成像光谱仪（Moderate-resolution Imaging Spectroradiometer），具有高辐射分辨率（36通道）、中等空间分辨率及高饱和度，覆盖范围广、时间分辨率为每天1～2次，具有很高的地震热异常监测精度. 用于研究的LST 数据来自于美国航空航天局（NSIDC，http：//nsidc.org/NASA/MODIS）网站下载的MODIS/Aqua 8d合成的MYD21A2夜间LST产品，空间分辨率为1 km，时间范围为2003—2019年（782幅图像）.

MYD21A2的LST数据，首先利用MYD21A2数据投影及格式转换工具（MODIS reprojection tool，MRT）进行拼接和重投影，投影转换为ALBERS，椭球体选为WGS84，图像文件转换为GeoTIFF格式；其次对于地震热异常提取，云层是一种干扰. 有云时MODIS数据或为陆表温度，或云顶温度，或混合温度，需要去除云的干扰. 所以，本文中地表温度数据经过云识别并做去云处理. 对于云像元干扰而形成的空值区域，为了解决此问题，低质量和缺失的值通常被排除在分析之外，或者被各种预测方法的预测结果所取代. 许多预测方法利用了数据的时间相关性，但迄今为止，只有很少人尝试利用数据的时空性质. 为了解决云干扰而形成的缺失值区域，本文依赖R包中的gapfill缺失值预测方法来替代地表温度数据集中的缺失值. gapfill方法是利用空间连贯性的特征趋势以及遥感数据集的时间季节性动态来预测缺失值. 此方法优点包括：①能够预测具有很大比例的缺失值（高达50%）的数据集中的所有缺失值. 该方法结果均方根相对误差（RMSPE，Root Mean Square Percentage Error）和平均绝对百分比误差（MAPE，Mean Absolute Percentage Error）都较低［18］. ②填补缺失值完成之后，利用STL分解（stats和raster包）、gesd法异常检测（Anomalize包）等处理检测遥感数据集中的异常信息.

3 研究方法说明

3.1 STL分解法

STL是非参数的统计方法，此方法把时间序列数据分成三项，分别是低频率的趋势项，高频率的周期项（季节项），不规则变化的趋势项：

图1 基于STL的地震热异常分解流程Fig.1 Workflow of earthquake thermal anomaly decomposition by STL

3.2 GESD检验

大多数情况下数据集中异常值不止一个，因此Rosner提出了GESD（泛化版ESD）. GESD（Generlized Extreme Studentized Deviate Test）是Rosner 教授基于Grubb’s Test（或ESD）改进的识别离散值的方法，此方法中比较统计量Rj与临界值λj，若Rj＞λj，该点为异常点.

式中：λj为第i 大的数据临界值；n 为离散点数量；tn-i-1,p是(n-i-1)自由度的t 分布右侧尾概率为p 时的值.

3.3 基于STL分解和GESD检验的异常检测算法

LST 的变化过程不仅具有明显的季节性变化，还有年际趋势变化的特征，是由多种因素共同作用的结果，不能直接对其进行监测分析. 因此首先利用STL分解法，取出地表温度空间数据集中周期分量、趋势分量，对剩余的残余分量（没有季节及趋势特征的时间序列）采用GESD检验进行异常监测. 此数据处理过程借助R语言所实现. 基于STL分解和GESD检验的异常检测算法的流程图如图2所示.

图2 异常检测算法流程图Fig.2 Flow chart of anomaly detection algorithm

4 震例结果与分析

4.1 平均LST时间序列的STL分解结果

图3是于田县Ms7.3地震震中附近2003—2019年LST的8 d平均时间序列数据用STL分解后的结果图.分解结果包括3 个部分：数据（data）反映时间序列数据的统计信息，季节项（seasonal）反映该时间序列的周期信息，趋势项（trend）为长期趋势，即该时间序列的发展方向，残余项（remainder）为去除周期信息以及长期趋势后的剩余部分. 从LST 夜间数据的STL 分解结果（图3）可以看出，LST 数据7 月或8 月份达到最大，最大值约-0.24 ℃. 每年的12月和次年的1月份达到最低值，最低值约为-31.59 ℃. 季节性变化是LST年内变化的重要特征，其变化范围为-11～12 ℃之间. 趋势分量的变化范围为-13～-15 ℃左右，研究区地表温度趋势的总体变化幅度为2 ℃左右，一般与全球气候变化有关. LST残余项的波动范围-6～5 ℃之间，可以看出，2003—2019年间LST残余项一直处于不规则的波动，而且波动幅度较大.

图3 震中附近LST时序STL分解结果图Fig.3 STL decomposition of time series data of LST around the earthquake center

通过对地表温度趋势项的分析可以看出，2003—2019年地表温度具有明显的变化. 为了明确地表温度趋势项的变化点，对研究区地表温度数据进行差分，文中第一步对地表温度趋势序列做一阶差分，初步得到其趋势变化点. 然后对一阶差分序列再做一次一阶差分，使数据保持在0值附近. 本文界限为二次差分序列的2倍标准差. 超过上界，下界就归为趋势发生变化. 看图4可知，2008年1月前后至2008年3月21日地震发生前地表温度具有上升趋势，2014年2月前后研究区地表温度偏增，呈上升趋势. 这两次地表温度上升趋势时间点临近地震发震期. 初步认为可能与研究区2008年与2014年7.3Ms地震有一定的关系.

图4 地表温度趋势项的二次差分结果Fig.4 The result of the quadratic differences of the surface temperature trend terms

4.2 LST异常面积的时间变化特征

图5为2003年1月至2019年12月发生异常像素数量的时间变化图. 从图5、6、7可以看出，2008年3月21日地震前，1月17日约2个月前出现异常的面积覆盖率达到19.5%，异常最明显. 之后1月25日的异常面积覆盖率为9.1%，2月2日的为16.6%. 2月2日之后异常面积覆盖率依次减少.

图5 震中附近LST时序热异常提取结果Fig.5 Results of thermal infrared anomalies extracted of LST time series around the earthquake center

同样，对于2014年2月12日地震，分析震前2个月，2013年12月份至2014年3月份研究区异常结果. 结果可知，异常较明显的是2013年12月27日研究区1.4%范围出现了异常，之后1月17日，1月25日异常覆盖率为1.8%，1.0%. 发震前4天（2月10日）异常面积覆盖率达到最大为2.1%.

4.3 LST异常的空间变化特征

以图6，2008年3月21日地震前1月至3月LST异常时空演化过程可见，空间异常范围主要出现在震中西北部. 随着发震期将近，空间异常范围扩展到震中东北部，异常区域面积1月25日达到最大，出现异常的面积覆盖率19.5%，2月18日起异常覆盖率逐渐减少.

图6 2008年于田Ms7.3地震热异常时空演化Fig.6 Space-time evolution processes of thermal anomaly before and after the Yutian Ms7.3 earthquake in 2008

图7 2014年于田Ms7.3地震热异常时空演化Fig.7 Space-time evolution processes of thermal anomaly before and after the Yutian Ms7.3 earthquake in 2014

2014年2月12日相隔6年于田县发生Ms7.3地震. 2013年12月至2014年2月地震热异常时空演化图如图7，此次地震发生在冬季，比较2008年地震出现异常面积不明显，异常范围比较零散，震中西北以及东北部出现小范围异常. 2月10日震前2天较小的异常区域出现在研究区东南部，且异常最明显，异常面积覆盖率达到2.1%. 发震之后，异常范围随着时间逐渐减少. 图8表示研究区2003—2019年地震热红外异常的空间分布，2003年6月10日最为明显，异常区域达到最大，热异常方向主要表现为东南部分.

图8 2003年6月10日的LST与热异常空间分布图Fig.8 Spatial distributions of LST and thermal anomaly on June 10，2003

5 讨论与结论

本文以2003—2019 年连续17 年的MODIS LST 数据，采用基于STL 分解和GESD 检验的异常检测法，对于田县两次Ms7.3地震，地震热异常演化过程进行了分析，得到如下结论：

1）以于田县地震为例，本论文中采用基于STL和GESD的异常检测. 首先利用STL分解法将时间序列进行分解，除了年变趋势及其季节因素周期性的影响后，用GESD 检验对STL 分解后的残余项进行异常检测.研究结果可见，震前研究区LST有明显的增温现象. 此方法对地震热异常的监测提供了参考依据.

2）对于2008年3月21日地震，异常区域面积1月25日达到最大，出现异常的面积覆盖率19.5%，2月18日起异常覆盖率逐渐减少. 对于2014年2月12日地震2月10日异常最明显异常面积覆盖率达到2.1%，发震之后，异常范围随着时间逐渐减少.

3）两次地震异常区域主要靠近震中西北及东北部，但2008年地震异常比2014年异常较明显，此研究区内异常范围较大.

总之，本文基于STL分解和GESD检验的异常检测算法提取了研究区地震热红外异常，探讨了此方法对地震热异常的监测有一定的参考价值，并总结了地震演化过程. 此外，由图可知，2003年6月10日异常达到最大，但此日前后没发生过地震. 热红外异常与地震的关联性可提高热异常对于地震的检测能力，这仍需要进一步研究.