FY-2G卫星红外遥感图像中的震前异常统计分析

乐应波， 陈福春， 陈桂林

1. 中国科学院上海技术物理研究所 红外探测与成像技术重点实验室, 上海 200083;

2. 中国科学院大学, 北京 100049

1 引 言

20 世纪80 年代，苏联学者分析热红外图像发现了震前热红外辐射异常（Gornyi 等，1988）。国内外学者对此进行了大量研究，包括机理研究、特征分析和算法等（申旭辉 等，2018）。强祖基等（1997）和Qiang 等（1998）提出了震前温度异常升高机理的初步模型：在高应力作用下，沿地震断裂带释放CH4、CO2和带电粒子进入大气，在太阳辐射和电场激发的作用下大气温度升高。这就是震前红外异常的理论模型。

随着卫星遥感技术和现代信息技术的发展，海量的热红外遥感数据和先进的信号处理方法被应用于提取震前红外异常信号（Bellaoui 等，2017；Jiao 等，2018；荆凤 等，2018；Ouzounov 等，2006；屈春燕 等，2006；张元生 等，2004）。目前，常用于地震红外遥感研究的数据主要来自美国发射的NOAA 系列卫星 （Ouzounov 等，2007）、EOS 系列卫星（Zoran 等，2014）和中国发射的风云系列气象卫星（Yao 和Qiang，2012）。静止轨道卫星可以大范围地连续观测地表和大气参数，在时间和空间上更具有可比性（Bormann 等，2003；Turk 等，2000）。Choudhury 等（2006）利用NOAA-AVHRR数据对伊朗的几次地震进行研究，发现该地区震前普遍存在大幅度的红外增温。

受到季节变化、天气条件、地质活动和人类活动的共同影响，无法直接观测到地质活动引起的红外增温，可靠的异常识别和分析算法至关重要（Shen 等，2013；Rivera 等，2012）。常见的震前红外异常提取算法包括RST （Robust Satellite Techniques） 算法（Tramutoli 等，2013）、四分位法、小波变换法、卡尔曼滤波法（Saradjian 和Akhoondzadeh，2011）、自回归整合移动平均法、人工神经网络和支持向量机（Akhoondzade，2013）和一些集成算法（Akhoondzade，2014）等。这些算法在个别典型震例中得到了验证，缺乏大量震例统计验证。在震前异常统计中，阳性预测值是异常现象中地震前兆的比例，用于表征地震预测的准确性；真阳性率是地震中有前兆的比例，用于表征震前异常的普遍性（Filizzola 等，2022；Zhang 等，2021；Fu 等，2020；Jiao和Shan，2021）。Lu 等（2016）利用FY-2E 卫星数据对西藏地区多年地震热红外异常进行了研究，结果发现5级以上地震震前普遍存在异常，但因未对所有异常信号进行统计，尚无法预测准确率；Zhang 和Meng（2019）对四川地区的热异常进行了统计分析，结果表明热红外异常与5 级以下的地震相关性不明显，与5级以上的地震有较低的阳性预测值和真阳性率，预测潜力有限，反映异常信号识别算法效果不太理想。

功率谱相对变化法是震例研究中常用的红外遥感地震信息提取算法（郭晓 等，2010；孟庆岩等，2017）。在典型震例的研究中，由于样本数量少，尚不能从统计学的角度验证算法的准确性和普适性；个别地区的统计分析中，存在研究区域小，无法观测到大面积异常现象的整体空间分布情况的问题（Zhang 等，2019；Fu 等，2020）。为此，本文提出一种基于连通域的方法对中国及周边地区的功率谱异常信号进行统计，从统计学的角度评估算法的准确性和普适性。该方法的核心思想就是将时间和空间上连续的异常体视为一个样本，计算不同参数条件下异常信号的阳性预测值和地震的真阳性率，并分析异常信号的幅度、空间范围特征和相关的地震信息。

2 数据和方法

2.1 功率谱相对变化法

本文所用的数据来自FY-2G 卫星搭载的多通道扫描辐射计，所用波段为长波红外波段（10.3—11.3 μm）。地表除了发出自身辐射外，还反射其他物体的辐射，为了避免太阳辐射的影响，选取午夜北京时间00：00—04：00 的亮度温度进行分析（张元生 等，2011）。原始数据经几何校正和辐射定标得到每个像素点对应的经纬度和亮度温度数据。热红外辐射难以穿透云层，并且云顶的温度远低于地面，在云层覆盖区域会形成孤立的低温区（Kato 等，2011）。本文中剔除低于平均值且差值超过1.5 倍标准差的数据，然后求每日的平均值作为当日值，时间序列中缺失的数据按近邻法进行插值，得到日亮度温度序列。

在日亮度温度数据中，地球温度场的年变化是低频信息，天气变化和人类活动引起的温度变化是高频信息（Ma 等， 2010）。小波分解可以分离信号中的高频信息和低频信息（Wang 等，2018；Wei 等，2013）。本文以7 阶小波分解的低频部分作为背景场，2 阶小波分解的高频部分可代表天气变化信息等，低频部分包含地震异常信息。用2 阶小波分解的低频信息减去7 阶小波分解的低频信息，可以在一定程度上消除季节和天气的影响，又包含有可能与地震相关的异常信号（张元生 等，2010）；最后，计算功率谱的相对变化值，从相对变化的时频数据中选取极大值作为特征频率，再寻找对应的时间序列。从而获得涉及到整个研究区域温度信息功率谱的时空分布数据。

2.2 异常信号识别

异常信号就是功率谱相对变化的时空分布数据中时间上和空间上连续的高幅值点集（Zhang等，2017）。每个点集被视为一个异常信号样本。

具体的算法流程如图1所示，包括二值化、形态学处理和连通域识别4个步骤。

（1）二值化：设定幅值阈值为A0，第d天第m行第n列的功率谱相对变化的幅值为P（d，m，n）。当P（d，m，n）≥A0，则认为该点为高幅值点，即δ（d，m，n）=1；当P（d，m，n）

（2）形态学处理：在数字图像处理中，闭运算是先膨胀后腐蚀的操作，用于填充空隙，连接相邻区域，在本文中用于消除高幅值区域内的小区域空隙，对比图1（b）和图1（c），闭运算增强了异常点之间的空间连续性，避免将相同因素引起的异常判断为多个异常。开运算是先腐蚀后膨胀的操作，可以消除小的离散点（Sundararajan，2017），在纤细部位分离两个区域，在本文中用于消除小区域的高幅值点，对比图1（c）和图1（d），开运算消除了部分小面积异常区域，避免将其识别为异常信号。

图1 功率谱异常信号识别算法Fig. 1 Recognition algorithm of abnormal signals in power spectrum

（3）连通域识别：对于一个区域D，若D中任意两个点A1（d1，m1，n1）和A2（d2，m2，n2）都可以用完全属于D的一条折线连接起来，则称区域D是连通区域。本文通过遍历和迭代的方法，识别所有高幅值点中的连通域。识别结果如图1（e）所示，图像中不同颜色表示不同的连通域，面积最大的连通域适合用于计算异常区域的位置和面积。

（4）样本选择：大气气团的覆盖面积为数十至数千公里，人类活动的影响范围约为数十公里（Song 等，2018）。Dobrovolsky 等（1979）依据理论模型得到孕震区域半径R与地震震级M之间满足R=100.43M。本文将覆盖面积较大的连通域标记为异常信号样本，剔除人为因素造成的小范围升温。如图1（f）所示，仅面积最大的连通域被标记为异常信号样本并用于下一步的统计分析。

2.3 统计分析

若在异常发生后t0天内，有5 级及以上的地震发生，且震中与异常区域的距离小于d0，则该异常信号可视为震前征兆，属于正样本，否则属于负样本。其中，t0为时间阈值，d0为空间阈值，反映地震预测的空间精度和时间精度。阳性预测值是被判断为正样本的样本中实际分类也为正样本的比例，反映地震预测的准确性；真阳性率指的是震前识别出异常的地震占所有地震的比例，反映算法预测地震的普适性（Zhang和Meng，2019）。阳性预测值（PPV）和真阳性率（TPR）如式（1）、（2）所示。其中，NA为异常信号样本总数，NP为正样本数量，NE为地震总数，N0为有震前异常的地震数量。

Molchan 图表法是经验性地震预测中常用的显著性检验方法（Zechar 和Jordan，2008）。图表的横坐标为异常时空占有率τ，即预测有震区域在研究时空范围内的占比；纵坐标为地震漏报率ν，即1 -TPR。其概率增益Gain如式（3）所示（蒋长胜 等，2011）。

为了研究震前异常的特征，本文对异常信号的空间分布参数进行了统计。异常信号样本中的峰值可以表征异常的幅度，幅度越大，则异常越明显（马未宇 等，2018）。由于异常区域并非规则的几何图形，本文用异常区域的最小外接矩形来描述异常的空间范围，如图2所示，最小外接矩形的4个顶角经几何校正后发生形变，用四边形的最长边长度来表示异常区域的长度。对正负样本的异常幅度和长度进行统计，可以直观地看出震前异常与其它异常的差异。震前异常与震级相关，也受到周边地形的影响（Jiao 和Shan，2022）。因此，本文对地震的时间、地点和震级进行了统计分析。

图2 异常区域长度Fig. 2 The length of abnormal region

3 研究区应用实例与结果分析

3.1 研究区选择与相关信息概述

选取2018 年功率谱相对变化的时空分布数据进行研究，该数据中共包含365幅图像，每幅图像包含1192819个像素点。整个数据集中的大部分幅值集中分布于0—5。

本文研究的区域为0°—60°N，70°E—140°E，包括了整个中国及其周边地区。在2018年2月1日至2019 年1 月1 日期间，该区域内发生5 级及以上地震共47 次，无7 级以上地震。震级最大的是2018 年12 月29 日11：39：08 在棉兰老岛附近海域发生的6.9级地震，震中位置为5.85°N，126.89°E，震源深度为50 km，地震前的相对功率谱图像如图3 所 示。地震前8 天，即2018 年12 月21 日，震中周边开始出现离散的高幅值点，范围逐渐扩大，强度逐渐增强；24 日到25 日期间，震中东北方向明显形成孤立的高温区域，然后逐渐消散；异常共持续了7天，在震前一天完全消失；震中距异常边缘约231 km，距异常峰值约683 km。根据Dobrovolsky（1979）的理论模型，5级以上的地震，受影响区域的半径应大于141.25 km，7 级以上的地震影响区域的半径应大于1023.29 km，此次异常出现的位置，在地震的影响范围内。本文中所研究的地震，其影响区域半径为数百公里。Tronin（2006）通过不同地区的卫星数据研究，发现地震前6—24天出现热异常，地震后持续约一周，且异常对大于4.5 级的地震敏感。因此，本文取幅值阈值A0分别为5、10、15、20、25，距离阈值d0分别为200 km、400 km、600 km、800 km、1000 km，时间阈值t0为30 d，阳性预测值和真阳性率如表1所示。A0越小，d0越大，则阳性预测值和真阳性率越大。约66%的地震在震前一个月内，出现幅值大于5的异常信号，且震中与异常区域的最短距离不超过400 km。阳性预测值较低，表明有大量的异常信号与短期内的地震无关。

图3 震前相对功率谱图像Fig. 3 Power spectrum image before the earthquake

本文通过Molchan 图表法 （蒋长胜 等，2011）进行统计检验，A0取15、20、25 的真阳性率偏低。A0取5 和10的统计检验结果如图4所示，概率增益最大为1.76，对应的阈值选择为A0=5，d0=400 km，对应的异常时空占有率为0.3753，漏报率为0.3403。通过表1 可知，该阈值条件对应的阳性预测值为0.2037（20.37%），真阳性率为0.6596（65.96%）。

表1 阳性预测值和真阳性率统计结果Table 1 The statistical results of positive predictive value and true positive rate

图4 Molchan图表法进行统计检验Fig. 4 Statistic test using Molchan diagram method

3.2 异常信号特征分析

为达到较高的空间精度和真阳性率，本文取A0=5，d0=400 km，t0=30 d。正负样本的峰值分布如图5 所示，长度分布如图6 所示，横坐标为异常信号的峰值和长度，纵坐标为峰值或长度在给定范围内的样本频率。正样本集中高幅度、大范围的异常信号所占的比例高于负样本集。异常样本的峰值和长度如图7 所示。峰值和长度小的样本中，绝大部分是负样本。峰值高于10，长度大于500 km 的异常信号样本，可达到33.71%的阳性预测值。峰值高于20，长度大于2000 km的异常信号样本，可达到80.00%的阳性预测值。

图5 异常信号峰值分布Fig. 5 The frequency distribution of abnormal peaks

图6 异常区域长度分布Fig. 6 The frequency distribution of the length of abnormal region

图7 异常的峰值和长度Fig. 7 The peak and length of anomalies

3.3 地震信息统计

在本文统计的47 次地震中，有25 次地震震级不超过5.4 级，不同震级的震前异常统计如图8 所示，横坐标为震级，纵坐标为对应震级中有震前异常和无震前异常的地震数量。5.4 级及以下的地震，有无异常的地震数量接近。22次大于5.4 级的地震中，仅有4 次地震在震前一个月内未发现异常，真阳性率可以达到81.82%。图9 为地震发生与异常消失的时间间隔分布图，横坐标为地震与异常结束的日期间隔，正数表示地震发生在异常消失后，负数表示地震发生在异常消失前，研究表明，大部分地震发生在异常消失后。图10 为震中的地理位置分布，日本、中国的台湾省、菲律宾和印度尼西亚位于环太平洋地震带上，这几个国家和地区在研究区域和时段内共发生25次地震，仅有4次地震在震前一个月内未发现异常。新疆维吾尔自治区、西藏自治区、青海、四川和云南位于地中海-喜马拉雅地震带上，该地震带上地震的真阳性率明显低于环太平洋地震带。西藏自治区和青海位于中国地壳厚度最厚的地区，在该地区的5次地震中仅有2次在震前一个月内发现了异常。

图8 不同震级的震前异常统计Fig. 8 Statistics of anomalies before earthquakes with different magnitude

图9 地震与异常消失的时间间隔分布Fig. 9 Time intervals distribution between the earthquake and the disappearance of the anomaly

图10 震中位置分布Fig. 10 Distribution of epicenters

4 结 论

本文通过功率谱相对变化方法提取红外遥感图像中的异常信号，并基于连通域方法对异常信号进行识别和统计分析。研究空间范围大，可以比较完整地观测大面积红外异常现象。本文在不同参数条件下进行统计分析和显著性检验，在A0=5，d0=400 km 的条件下在得到20.37%的阳性预测值和65.96%的真阳性率，其概率增益为1.76。本文的结果验证了功率谱相对变化法能在大部分地震前提取到红外遥感图像异常，但阳性预测值偏低，预测的准确率有待提升。本文对研究区域异常信号的特征参数和地震信息进行统计分析，得到以下结论：

（1）功率谱相对变化法提取出的高幅度、大范围的异常信号（峰值大于20，长度大于2000 km），可以达到80.00%的阳性预测值，说明剔除小幅度小范围的异常信号可以在一定程度上提高地震预测的阳性预测值；

（2）对于大于5.4 级的地震，真阳性率可以达到81.82%，说明本文的方法对大于5.4 级的地震敏感；

（3）真阳性率具有明显的地区差异，环太平洋地震带的地震真阳性率要高于地中海-喜马拉雅地震带。

本文提出的方法适用于大范围、长期的震前异常统计和特征分析，解决了样本数据少、研究区域小的问题。在此基础上，可以进行更深入的研究，包括但不限于以下几个方面：

（1）增加数据样本，分类统计，研究不同地震类型的震前异常情况。

（2）提取更多的异常信号特征，分析正负样本的差异，提高阳性预测值。

（3）结合多种数据源和数据处理方法，对比不同算法的准确性和普适性，探索更可靠的震前征兆。

志 谢本文所用部分静止气象卫星数据资料来自中国气象局国家卫星气象中心，感谢郭强研究员所提供的帮助和指导。