吴玮莹 单新建 屈春燕 李新艳
1)中国地震局地质研究所,地震动力学国家重点实验室,北京 100029 2)宁夏回族自治区地震局,银川 750001
地震是造成财产损失和人员伤亡最为严重的自然灾害之一。预测未来地震发生的时间、 空间、 强度信息是预防地震灾害最为有效的手段(Mukhopadhyayetal.,2019)。随着遥感技术的发展,大量利用热红外卫星数据的地震异常研究表明: 地震前存在地表温度(文翔等,2015; Ahamdetal.,2019; Khalilietal.,2019; Chenetal.,2020)、 长波辐射(Xiongetal.,2010; Kongetal.,2018; Zhongetal.,2020)等多种热参数的异常变化现象,热参数的时空分布特征对未来地震的震中、 发震断层具有一定的指示作用。
由于热参数受其他因素(如局部水热条件、 季节、 气候等)影响,其本身的变化具有周期性、 偶发性的特点。因此,采用合适的地震热参数异常提取方法抑制其他因素对热参数的影响、 突出地震因素的作用,对于地震热参数异常变化研究是十分重要的。目前,ZS法(Z-score,也称偏移指数法)和RST法(Robust Satellite Technology)是应用较为广泛的地震热参数异常提取方法(Filizzolaetal.,2004; Tramutolietal.,2005; Ouzounovetal.,2011; Genzanoetal.,2015; Khalilietal.,2019)。ZS法和RST法均利用卫星多年观测的热参数值的平均值作为热参数正常变化的期望值,以多年观测值的标准差描述观测值的可变性以构建背景场,以当前观测值和背景场的差作为分母部分,多年的标准差作为分子部分,利用该比值定义热参数变化。不同的是,RST法在构建背景场时,通过划分地物覆盖类型并求取不同地物类型的平均值,以降低地物覆盖类型对最终结果的影响。ZS法和RST法都是基于统计分析的方法,原理清晰明确,操作较为简便。
尽管大量研究成果表明地震前存在热参数时空分布特征的异常变化,但由于不同地震的热参数异常提取方法对同一参数提取的结果存在差异,使得热参数异常的时空特征及其与地震活动的关联受到质疑。因此,需要对不同地震热参数异常提取方法的适用性和优缺点进行比较分析,但目前尚缺乏相关研究(Jiaoetal.,2018)。 本研究以2014年MW6.9 于田地震为典型震例,以2008年汶川MW7.9 地震为验证震例,同时利用ZS法和RST法提取典型震例震前3个月的地表温度(Surface temperature,ST)和晴空条件下向外出射长波辐射(Outgoing longwave radiation,OLR)数据异常变化的时空特征,定性分析提取结果之间的差异。随后,利用地震平静年(无震级>5.0级地震发生的年份)数据对2种提取方法在异常变化敏感性、 对背景信息的抑制能力和对地震信息的指示性3个方面进行了定量比较,并利用验证震例验证定量化比较的结果,详细探讨了造成差异的原因,总结了各方法的优缺点和适用性,为后续进一步探究地震震前热参数异常的时空分布及其可能的物理机制奠定基础。
本文主要对ZS法和RST法对同一震例、 不同热参数在震前时空变化特征提取方面的适用性和可靠性进行对比研究,选取地表温度(Surface temperature,ST)和晴空条件下向外长波辐射(Clear-sky Outgoing longwave radiation,OLR)作为研究参数。地表温度是地表和大气边界层之间能量交换的重要影响因素,能够更直观地反映地震孕育过程中的热辐射变化。晴空条件下向外长波辐射是地球表面(陆地和海洋)、 大气反射、 吸收和发射的长波辐射最终离开地球进入太空的辐射能量,是整个地面-大气系统热辐射的综合反映。
热参数数据来自欧洲中期天气预测中心(ECMWF,European Centre for Medium-Range Weather Forecasts)第5代再分析数据集ERA5(1)https: ∥cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-single-levels?tab=overview。。该数据集是通过数据同化和再分析技术将先进的预报模型得到的预报结果与多种观测手段得到的大量观测数据结合为精度更高、 覆盖更广的数据集。与传统的卫星遥感数据(如MODIS(Moderate-resolution Imaging Spectroradiometer)数据)相比,再分析数据具有以下优势: 1)时间序列完整,目前能够提供从1979年至今的完整数据集; 2)空间上无空值,能够完整反映热参数的空间分布特征。与其他再分析数据集(如NCEP(National Centers for Environmental Prediction)数据)相比,ERA5数据集的空间分辨率(0.25°×0.25°)更高。
由于地震震前热异常往往发生在震中和断裂带附近的大尺度区域内,因此,我们采用Dobrovolsky等(1979)提出的地震影响区范围估计公式(式(1))划定研究区域的范围:
r=100.43M
(1)
其中,r为地震影响区域的半径(单位为km),M为地震震级。根据计算结果,2014年于田地震研究区域的范围为30.91°~40.91°N,77.59°~87.59°E,2008年汶川地震的研究区范围为26°~36°N,98.3°~108.3°E。 2次地震的研究时段为地震前3个月—地震当月,共获取13a的观测数据。除地震年外,其余年份内均无震级>5.0级的地震发生。热参数数据为夜间平均值(北京时间00:00—03:00,UTC: 16:00—19:00),以避免太阳辐射和人为活动的影响。
ZS法和RST法是目前地震热参数异常时空变化特征研究中较为常用的2种地震热参数异常提取方法。ZS法在图像差分的基础上引入了信噪比的概念,将地震活动区当前观测值与多年观测数据平均值之间的差值作为信号部分,并将多年观测数据的标准差作为噪声部分。该方法可表示为(Ouzounovetal.,2011)
(2)
式中,v(x,y,t)是地表温度的当前观测值,μ(x,y)为背景场,即同一位置多年同期观测数据的平均值,δ(x,y)为相应的标准差。ZS指数的大小不仅反映了当前观测值偏离多年背景值的程度,同时还反映了这种偏离程度相对于多年噪声水平的强弱。
RST法利用在类似观测条件(如同月、 同时、 同传感器)下获得的多年观测数据集的预期值(即多年平均值)及可变性(多年标准差)定义被测信号的异常变化。RST法将邻域差分与ZS法中的背景场归一化相结合,明确定义了统计学中地震前异常的数学表达式,热红外异常变化通过特定的指数: REIRA(Robust Estimator of TIR Anomalies)定义(Genzanoetal.,2015; Ouzounovetal.,2018; Khalilietal.,2019)。RST法可表示为
(3)
本文同时采用ZS法和RST法对2014年于田地震和2008年汶川地震震前的地表温度和长波辐射异常时空变化特征进行提取分析。具体步骤包括(图 1): 1)数据下载和预处理(包括数据格式转换、 重投影等); 2)对用于构建背景场的数据进行滑动平均以剔除偶发因素对背景场的影响; 3)针对研究区的地物类型,在RST法中,大致将于田地震和汶川地震研究区分为山区和盆地2种地物类型,并分别计算这2种地物类型各自的平均值; 4)按照ZS法和RST法的算法原理分别计算多年数据的平均值和标准差作为背景场; 5)按照式(2)和式(3)计算得到ZS指数和REIRA指数; 6)为完整展示研究时段内热参数异常变化的过程,我们参照MODIS等数据的形式(MODIS产品提供8d平均值数据),以8d为一个阶段绘制三维ZS指数和REIRA指数的平均值图像。
图 1 地震异常提取的流程图Fig. 1 Flow chart of seismic anomalies extraction.
在ZS法和RST法的实际应用中,通常假设数据集具有近似高斯分布的概率分布。在这个假设下,68.3%、 95.4%和99.7%的值分布在1、 2、 3个标准差的范围内(Jiaoetal.,2018)。为在统一异常阈值下得到可对比的异常提取结果并讨论不同异常强度阈值对异常提取结果的影响,我们分别将异常强度阈值设定为2和3,此时的置信区间分别为95%和99%。
本文以2014年于田地震为典型震例,因此本节将着重介绍2014年于田地震的概况及区域构造背景。据中国地震台网测定,北京时间2014年2月12日17时19分49秒,中国新疆维吾尔自治区和田地区于田县发生MW6.9 地震,震中位于无人区内,坐标为(36.1°N,82.5°E),震源深度12km,震源机制为左旋走滑。截至北京时间2014年2月13日14:00:00,共记录到余震1097个,其中MS5.0 ~5.9地震1个,MS4.0 ~4.9地震12个,MS3.0 ~3.9地震22个。
2014年于田地震发生在阿尔金断裂带的西南段(图 2)。阿尔金断裂带位于青藏高原西北边缘,是亚洲大陆内最长的走滑型活动断裂带,呈NEE走向,长约2000km。阿尔金断裂带分割了青藏高原和塔里木盆地,并划分了青藏高原北部不同的构造单元,对吸收印度-欧亚大陆碰撞而产生的地壳缩短具有重要作用。根据断裂带属性,阿尔金断裂带可分为3个部分: NNE—E向逆冲的喀拉喀什断裂、 E—NEE走向的阿什库勒断裂和NE—NEE走向的硝尔库勒断裂(Yinetal.,2002; 徐锡伟等,2011; Lietal.,2015; 宋春燕等,2015; Lietal.,2016)。该地区历史上共发生过2次震级>7级的地震,分别为1924年阿尔金断裂带西段7.2级地震,震中位于2014年于田地震震中东北约250km处; 2008年于田MW7.2 地震,震中位于2014年于田地震震中西南150km处。
图 2 研究区的区域构造背景图Fig. 2 Tectonic map of the study area.
图 3 2014于田地震震前3个月的地表温度和长波辐射ZS 8d平均值的三维时空特征分布图Fig. 3 Surface temperature(ST) and outgoing longwave radiation(OLR) three-dimensional maps based on 8-day averaged Z-score index values for 3 months before the 2014 Yutian earthquake.图m为包含发震时段的结果
图 4 2014年于田地震震前3个月的地表温度和长波辐射REIRA 8d平均值的三维时空特征分布图Fig. 4 Surface temperature(ST) and outgoing longwave radiation(OLR) three-dimensional maps based on 8-day averaged REIRA index values for 3 months before the 2014 Yutian earthquake.图m为包含发震时段的结果
图 5 地表温度(a)和长波辐射(b)敏感性测试数据时间序列图Fig. 5 Time series plot of sensitivity test data for ST(a) and OLR(b).
图 3 展示了研究时段内地表温度(ST)和晴空向外长波辐射(OLR)ZS指数8d平均值三维图像,其中红色部分表示ZS指数>2的区域,红色矩形框代表地震的发震时段。从图 3 中可以看出,在震前3个月,ST和OLR都出现多次高值异常。从时间序列上来看,随着地震的临近,ST和OLR异常变化的频次逐渐增加; 高值异常最早出现在地震前2个月(图3f),而地震前1个月出现持续的高值异常(图3j,3k,3l)。从空间分布上来看,随着地震的临近,ST和OLR的高值异常区域由研究区的北部和西南部(图3f,3j)逐渐向地震震中和阿尔金断裂带周围集中(图3k,3l),在包含发震时刻的时段内,ST无明显高值异常,而OLR的高值异常区则明显分布于震中周围。
图 4 展示了研究时段内ST和OLR的REIRA指数8d平均值的三维图像,其中红色部分表示REIRA指数>2的区域,红色矩形框代表地震的发震时段。对比ZS法和RST法的提取结果可发现: 1)对比图3f、 4f,3k、 4k,3l、 4l可知,RST法与ZS法的提取结果在地理空间分布上具有相似性,基本分布于研究区的同一位置; 2)在地震震前3个月,2种方法的提取结果差异较大(图3a—e,图4a—e),ZS法得到的结果无明显异常变化,而RST法得到的结果则有异常变化且分散分布于整个研究区内; 此外,2种方法的提取结果在异常出现的频率和幅度上明显不同,RST法得到的异常频次和幅度都高于ZS法的提取结果。为何使用同样的数据会得到不尽相同的时空分布特征?哪种方法得到的时空分布特征更为可靠?为了回答以上问题,我们在下节中针对2种方法的异常变化敏感性、 对背景信息抑制能力和对地震信息指示性3个方面进行探究。
热参数时空特征的异常变化受到地形、 气候变化、 人为活动等多种因素的影响。相较于其他影响因素,地震孕震过程对热参数异常变化的影响十分微弱。地震热参数异常提取方法的敏感性代表了其识别地震孕震过程中引起的热参数微弱异常变化的能力。本文选取2014年于田地震前(2009年)和地震后(2015年)2个无震级>5级地震发生年份的震中所在像元的时间序列数据作为测试数据(图5a 为ST数据,图5b 为OLR数据),在时间序列中随机选取20个数据值,人为设定1~15(单位为K或W/m2)的异常变化,随后利用2种方法提取异常,并统计在不同异常阈值和异常变化幅度下检测出的设定异常个数及其准确率(即检测出的设定异常个数与总异常个数之比)。
图 6 为利用2009年地表温度(a、 c)和长波辐射(b、 d)数据通过ZS法和RST法提取得到设定异常提取天数(a、 b)和准确率(c、 d)。图 7 为利用2015年地表温度(a、 c)和长波辐射(b、 d)数据通过ZS法和RST法提取得到的设定异常提取天数(a、 b)和准确率(c、 d)。异常阈值的设定对异常变化检测具有较大影响,当设定的异常阈值较高时,虽然不能检测出全部人为设定的异常变化,但异常检测的准确率较高。对比2种方法得到的结果可知,ZS法和RST法都对小幅度的异常变化具有一定的敏感性(异常变化幅度为1~5K(W/m2)时),但RST法的异常检测准确率优于ZS法。当异常变化幅度较大时(12~15K(W/m2)),利用2009年的数据,RST法的提取效果更好; 利用2015年的数据,当异常阈值为2.5时,ZS法的提取效果也较好,这可能是由于提取参数不同所致。综合上述研究成果,我们认为ZS法和RST法都对于小幅度的异常变化具有一定的敏感性,但RST法的提取效果优于ZS法。
图 6 利用2009年地表温度(a、 c)和长波辐射(b、 d)数据通过ZS法、 RST法提取得到设定异常提取天数(a、 b)和准确率(c、 d)Fig. 6 Artificial anomaly extraction days(a and b)and accuracy rate(c and d)obtained from 2009 surface temperature(a and c)and outgoing longwave radiation(b and d)data extracted by ZS and RST methods.
图 7 利用2015年地表温度(a、 c)和长波辐射(b、 d)数据通过ZS法、 RST法提取得到设定异常提取天数(a、 b)和准确率(c、 d)Fig. 7 Artificial anomaly extraction days(a and b)and accuracy rate(c and d)obtained from 2015 surface temperature(a and c)and outgoing longwave radiation(b and d)data extracted by ZS and RST methods.
地震热参数异常提取方法对背景信息的抑制能力决定了提取结果的准确性。本文提取了2个无震年(2009年和2015年)的地表温度和长波辐射异常时空变化特征。无震年的数据没有受到地震活动的影响,因此理想的异常提取结果应无明显异常变化。然而,由于异常提取方法无法完全剔除其他因素造成的影响,故提取结果中可能有少量由其他因素引起的热异常。本研究以无震年的异常像元密度,即异常像元个数与总像元个数的比值定量描述异常提取方法对背景信息的抑制能力。我们统计了利用ZS法提取无震年异常像元密度大于RST法提取的天数和利用RST法提取无震年异常像元密度大于ZS法提取的天数,分别记为ρZR和ρRZ。无震年异常像元密度越低,则该方法对其他因素引起的热异常的抑制能力越强,相应地,其在地震年提取出的热异常时空变化特征可能越可靠。
图 8 为2009年(a、 c)和2015年(b、 d)地表温度(a、 b)和长波辐射(c、 d)异常像元密度在不同异常标准下的统计结果。无论异常标准设定为多少,ρZR的值远小于ρRZ,说明ZS法在无震年得到的异常像元明显少于RST法,基于2009年和2015年数据得到的结果呈现类似的特征。由此可知,ZS法相较于RST法对其他因素引起的热参数异常变化的抑制作用更强。
图 8 2009年(a、 c)和2015年(b、 d)地表温度(a、 b)和长波辐射(c、 d)异常像元密度在不同异常标准下的统计结果Fig. 8 Statistical results of surface temperature(a and b)and outgoing longwave radiation(c and d)anomaly pixel densities under different anomaly standards in 2009(a and c)and 2015(b and d).
前人研究表明,地震震前热异常的空间分布与震中位置具有一定的相关性(Jiaoetal.,2018)。因此,地震震前热异常的空间分布是预测地震信息(发震时段和震中位置)的重要指标。本文在第3节中说明了2种方法提取的异常结果存在时空分布差异。那么,哪种方法提取得到的热异常空间分布能够更好地指示地震信息?为进一步探究该问题,我们根据地理学第一定律,即“所有事物都与其他事物相关,但近处的事物比远处的事物更相关”(Tobler,1970),通过归一化距离指数(Qietal.,2020)定量评估热参数异常时空分布对地震信息(发震时段和震中位置)的指示性作用(式4):
(4)
式中,(xi,yi)为异常像元坐标,(xepi,yepi)为震中所在的像元坐标,n为整幅图像异常像元总数。NWI值越大,说明得到的空间分布越集中于发震震中,指示性越好。
在对地震信息指示性进行对比时,需要同时关注2个方面: 1)在都具有异常像元分布时,哪种方法计算得到的NWI值更大; 2)哪种方法计算得到NWI最大值的时间离发震时段更近。因此,我们首先统计了当ZS法和RST法都具有异常像元分布时,其中一种方法计算得到的NWI值大于另一种方法的天数,然后对地震震前2种方法计算出的NWI最大值距地震发震时段的天数进行了统计。图 9 为基于地表温度(图9a)和长波辐射(图9b)数据的ZS法和RST法对地震震中指示性的比较结果图。从图中可以看出,ZS法和RST法基于不同数据得到的结果存在不同。对于ST数据,当同时具有异常像元分布时,ZS法计算得到的NWI值略高于RST法计算得到的NWI值的天数更多,特别是当异常标准设定为3时。然而,对于OLR数据,RST法计算得到的NWI值高于ZS法计算出的NWI值的天数更多。这一现象说明对于ST数据,ZS法对于地震震中分布的指示性略好于RST法; 而对于OLR数据,RST法得到的结果对地震震中空间分布的指示性更好。
表1为在震前时段内使用ZS法、 RST法计算得到的NWI最大值出现时间距地震发生时的天数。从表中可以看出,相较于RST法,ZS法得到的NWI最大值出现时间更靠近地震的发震时段,对地震发震时段的指示性更强。综合上述研究成果表明: 1)针对ST数据,ZS法对于地震震中位置的指示性略好于RST法,对于OLR数据则是RST法更优; 2)ZS法得到的NWI最大值出现时段距离发震时刻更近。
图 9 基于地表温度(a)和长波辐射(b)数据的ZS法和RST法对地震震中指示性的比较结果图Fig. 9 Comparison of epicentre indication ability of surface temperature(a)and long wave radiation(b)obtained by ZS and RST methods.
表 1 ZS法、 RST法计算得到最大NWI值时段距离地震发震的天数(数字表示震前天数)Table1 The days from the maximum NWI value period to the occurrence of earthquake calculated by ZS and RST methods (The number represents the days before the earthquake)
表 2 ZS法、 RST法针对验证震例计算得到的最大NWI值时段距离地震发震天数(数字表示震前天数)Table2 The days between the maximum NWI value period and the earthquake occurrence calculated by ZS and RST methods from the verification earthquake case(The number represents the days before the earthquake)
据中国地震台网测定,2008年5月12日14点28分(北京时间)四川省阿坝藏族自治州汶川县发生MW7.9 地震,震中坐标为(31.0°N,103.4°E),震源深度14km,是21世纪以来发生在中国大陆境内破坏性最强的地震。我们以此次地震作为验证地震,利用汶川地震震前平静年的数据(2007年)作为测试数据,采用与于田地震相同的方法,对ZS法和RST法在异常变化敏感性、 对背景信息抑制能力和对地震信息的指示性3个方面进行对比,验证以于田地震作为典型震例的对比结果的可靠性。ZS法和RST法的对比策略如4.1、 4.2和4.3节所述。
图 9 为验证震例异常变化敏感性(a、 b)、 对背景信息抑制能力(c、 d)和对地震信息指示性(e、 f)的测试结果图。表2 为NWI最大值出现时段距地震发震的天数统计。利用汶川地震数据得到的异常变化敏感性结果与于田地震略有不同,当异常阈值为2时ZS法对于微弱异常变化敏感性最好。2种方法对背景信息抑制能力的测试结果与于田地震类似,ZS法对于背景信息的抑制能力优于RST法。在对地震信息指示性的测试中,总体而言,RST法提取得到的异常空间分布对地震震中的指示性更好,而ZS法对于地震发震时段的指示性更好,NWI最大值出现的时间更靠近地震发震时段。
图 10 针对验证震例的异常变化敏感性(a、 b)、 对背景信息抑制能力(c、 d)和对地震信息指示性(e、 f)的测试结果图Fig. 10 The test results of sensitivity to anomaly change(a and b),ability to suppress background information(c and d)and indicate seismic information(e and f)of the verification earthquake case.
本研究以2008年于田地震为典型震例,以2008年汶川地震为验证震例,对ZS法和RST法在对异常变化的敏感性、 对背景信息的抑制能力和对地震信息的指示性3个方面的提取效果进行了定性和定量化比较。除异常变化敏感性外,在对背景信息的抑制能力和对地震信息的指示性上,2个震例得到的对比结果类似; ZS法对背景信息的抑制能力更好; RST法对地震震中的指示性更好; ZS法对地震发震时段的指示性更好。
图 11 ZS法和RST法信号部分和噪声部分在整个研究时段内每天最大值与最小值的差异时间序列曲线Fig. 11 The time series curve of difference values between the maximum and minimum values of the signal part and the noise part of the ZS method and the RST method for each day during the whole study period.a、 b 信号部分的时间序列曲线; c、 d 噪声部分的时间序列曲线; a、 c 利用ST数据得到的结果; b、 d 利用OLR数据得到的结果
为何 2 种方法在这3个方面具有显著的差异?如何说明在第3节中2种方法得到的异常变化时空特征的不同?这些问题仍需要进一步分析和讨论。ZS法和RST法都是以指数形式表示的异常提取方法。在计算过程中,2种方法的不同之处在于RST法针对区域进行了地物划分,求取了不同地物类型的平均值。在构建背景场时,首先计算多年数据中每个像元与其所在的地物类型平均值的差值,之后再计算多年平均值和标准差以构建背景场。ZS法则直接计算多年平均值和标准差作为背景场。为确定求取每个像元与其所在地物类型平均值的差值对最终信噪比结果的影响,我们统计了ZS法和RST法信号部分(即地震年与多年平均值的差值,如图11a 和11b所示)和噪声部分(即多年标准差,如图11c 和11d所示)在整个研究时段内每天的最大值与最小值的差异。从图 11 中可以看出,求取每个像元与其所在地物类型平均值的差值对信号部分影响不大,对噪声部分的影响较大。整体上,ZS法噪声部分的最大值、 最小值的差异大于RST法。因此,在异常变化幅度较小时,由于RST法计算得到的信噪比大于ZS法,使得前者对于异常变化更为敏感,在相同的异常标准下能够提取出更多的异常变化,且在地震平静年也能提取出大量异常。当ZS法提取结果无明显变化时,RST法也能提取出分布于研究区内的异常; 当2种方法提取结果同时存在异常变化时,RST法提取结果的空间分布面积更大。ZS法和RST法的异常变化敏感性在不同震例中得到了不同的比较结果,是由研究区地物分布类型不同所造成的。对于田地震而言,区域内两大地物类型——山区和沙漠的分布较为均一,不同地物类型的平均值差异较大。而汶川地震研究区内的地物类型则更为复杂,不同地物类型的平均值差异较小,影响了RST法对于微弱变化的敏感性。由于热参数影像的空间分辨率较大,一个像元内包含了大量地物,无法利用光学影像进行精确的地物类型分类,这也降低了RST法提取结果的可靠性。此外,异常标准的设定对于提取结果的影响也很大,但目前并没有有效定义异常标准的方法,大多数研究都依据经验定义异常标准,因此需要进一步研究定义异常标准的有效方法,以提高热异常提取的可靠性。
综上所述,我们认为ZS法和RST法对于异常变化都具有一定的提取能力,但相比较而言,由于RST法对于背景信息抑制能力弱,且在对地物类型进行分类过程中存在较大的人为因素影响,因此,ZS法是较为合适且简便的地震异常变化提取方法。