宋冬梅 时洪涛 单新建 刘雪梅 崔建勇 沈 晨屈春燕 邵红梅 王一博 臧 琳 陈伟民 孔 建
1)中国石油大学(华东),地球科学与技术学院,青岛 266580
2)中国地震局地质研究所,地震动力学国家重点实验室,北京 100029
3)中国石油大学(华东),理学院,青岛 266580
4)武汉大学,测绘学院,武汉 430079
21世纪以来,中国乃至世界各地地震频发,造成大量人员伤亡和财产损失(王贵武等,2006),地震预测工作对保障中国经济建设和社会繁荣稳定具有非常重要的意义。几十年来,人们一直尝试寻找地震前兆信息来预测地震,虽然取得了一些进展,但地震预测一直是困扰各国科学家的一个世界性难题(徐程等,2012;姜金征等,2013;Lazaridou-Varotsos,2013)。
目前的地震预测方法大致可以分为3类:地震地质方法、地震统计方法和地震前兆方法。其中前兆法是常用方法之一,涉及测震、地形变、重力、地电、地磁、地下水化学、地下水物理、气象和地温等多个指标(强祖基等,1998;马瑾等,2000;陈顺云等,2004;曾中超等,2009)。震前地表温度异常即是地震前兆异常之一。随着遥感技术的发展,其信息量大、范围广、可实时监测和客观准确等特点使得利用热红外遥感数据监测震前地表温度异常成为可能(屈春燕等,2007;张元生等,2010)。近年来,在利用热红外遥感数据提取地震热异常信息的算法研究及其与地震关系方面的探索取得了一些重要进展(张微等,2008;强祖基等,2010;张璇等,2013)。大量研究结果表明,地震发生前存在不同程度的“热震兆”现象(强祖基等,1998;Tronin et al.,2002;屈春燕等,2006a;Arun et al.,2008;Ma et al.,2011;Yao et al.,2012)。大量试验表明(强祖基,2000;吕月琳等,2009;闫丽莉,2012),遥感热红外信息在地震短临预测中有其独特的优势,是探索地震预报的可能途径之一。
众所周知,地震现象是一个非常复杂的过程,异常与地震之间有较强的不确定性。因震前各类异常与未来地震之间有较强的非线性关系,很难通过某种解析表达式来表达,即很难用某种线性关系进行拟合与描述(屈春燕等,2006b;林德明等,2009)。而容错能力强、自组织、自学习,且具有很强的非线性拟合能力的神经网络在地震预报方面有着独特的优势(杨居义等,2008;Dhawal,2011)。张金华等(2012)以BP神经网络为基础建立短期地震预测模型,较好地对未来几个月的最大震级进行了预测。李炜等(2011)对实验区地震数据时间序列进行分析并对大地震的发生时间进行预测。聂仙娥等(2011)利用BP神经网络对地震震级的预测也取得了较好的效果。然而,以往的预测研究仅是单一地对发震时间或震级进行预测,本研究尝试将热异常信息作为神经网络的输入,建立热异常信息-地震三要素的时空框架,以此指导BP神经网络结构的构建,并以中国境内及周边2002年以来100个5级以上震例和70个随机无震样本对神经网络进行训练与仿真,实现地震三要素的同步预测试验。
本文采用2002至2011年中国及周边发生的100个5级以上地震(主要分布在中国中西部地区和中俄边境等地),以及70个随机无震样本进行研究,包括2010年4月14日青海省玉树7.1级地震(96.6°E,33.2°N)、2008年 10月 5日新疆乌恰县 6.8级地震(73.9°E,39.5°N)、2008年5月12日四川汶川8级地震(103.4°E,31°N)、2008年3月21日新疆于田县7.3级地震(81.6°E,35.6°N)以及2011年3月24日缅甸7.2级地震(99.8°E,20.8°N)等震例。
在本研究中,用于热异常信息提取的遥感数据是8d合成的1km分辨率夜间MODIS地表温度数据(LST),即MODIS/Terra Land Surface Temprature/Emissivity 8-Day L3 Global 1km SIN Grid V005。数据来源网址:http:∥reverb.echo.nasa.gov/。数据下载后得到.hdf文件,再由MRT软件转换为.tif图像,并进行裁剪得到单个震例图像(图1)。
本文的实验流程主要包括数据预处理和BP神经网络的预测试验两大部分。
数据预处理包括MODIS 8d的LST数据的下载、转换、裁剪以及利用RST算法对热红外异常信息的提取;BP神经网络的预测试验包括热异常时空坐标系的构建、BP神经网络输入输出的确定以及BP神经网络的训练和仿真;最后对预测结果进行分析(图2)。
数据预处理包括MODIS数据下载、格式转换和图像裁剪等步骤。根据MODIS中国图像分区代码,选择中国范围内的分区,进行MODIS数据下载。下载完成后,用MRT软件将.hdf文件转换为.tif图像,并根据单个震例的数据范围,使用ENVI进行相应的裁剪。
图1 玉树2010年3月4号的MODIS图像Fig.1 The MODIS image of Yushu on March 4,2010.
图2 技术路线框架图Fig.2 Technical route frame.
经过多年研究,目前已有多种关于热异常提取的算法(杨杰,2013),它们对异常的提取各具一定的作用,但各种算法所考虑的因素不同。本文使用的RST算法不仅考虑到季节交替和气候变化对亮温的影响,而且考虑到正常年份云覆盖和地震等因素的影响,使得结果更加可靠(吴立新等,2009;温少妍,2011)。
RST(Robust Satellite Techniques)方法是一种基于多时相数据分析方法,考虑的是时空域上相对于平静状态下的异常。在同一条件下,使用长时间序列的卫星亮温数据来构建信号的正常背景并计算局部偏差指数Alice值(即热异常值)。当局部偏差指数Alice值大于某一阈值时,可用来识别干扰事件,即地震事件。Alice值是基于像素级水平定义的,其具体表达式为式(1)—(4)中,T(ri,t)是t时刻像元ri=(xi,yi)处的亮温值(单位,K);ΔT(ri,t)为ri处像元的亮温值与研究区平均亮温T(t)的差值(单位,K);μΔT(ri)与σΔT(ri)分别为研究区内历年同期有效(去除云干扰)的ΔT(ri,t)的均值和标准差;S(ri,t)=C(ri,t)·E(ri,t)表示在无云和平静期,与影像获取时刻t相一致的像元ri的标示值,无量纲;t为长时间序列影像每景的观测时间。
用RST算法提取的Alice指数值结果如图3所示。
图3 2010年4月14日玉树地震前2个月的遥感热异常提取图像Fig.3 The thermal anomaly image two months before Yushu earthquake on April 14,2010.
2.3.1 热异常时空坐标系的构建
地震事件的三要素与热异常首次出现的时间、异常面积及异常强度有关。进行基于热异常的神经网络地震综合预测,须首先建立起一个能够表达热异常时空特征的坐标系,确定空间原点和时间起点,以体现热异常信号的游走路径及时空演变规律,确保神经网络的训练结果具有明确的物理意义。
本文在预测试验中,以震中为中心选择10°×10°矩形框为研究区,矩形框的左下角为空间直角坐标系的坐标原点;以震前2个月为数据的时间搜索范围,以热异常首次出现的时间为时间起点(图4)。对于未来实际地震预测,震中位置未知,可根据活动断裂带信息和地震地质构造情况结合专家经验知识,选取矩形框为研究区,确定出空间原点;在预测试验的当天前2个月的时间范围内,找到热异常首次出现的日期作为时间起点。若在此2个月中无热异常现象出现,则不必采取本方法进行地震预测。
震例研究区内,存在热异常区域的面积之和在某时刻会出现最大值,即为热异常最大面积(单位,km2);最大异常强度是指某时刻热异常区域中像素的Alice值之和的最大值(无量纲)。以上对热异常信息的定量描述均作为BP神经网络预测模型的输入层参量。
为研究热红外异常的时空演变路径及发展特点,本文引入热异常区域的几何能量质心)这一概念,其计算公式为
2.3.2 以热异常为输入的地震预测的BP神经网络构建
神经网络的原理:BP神经网络是误差反传网络(Back-Propagation Network,BP),也是一个以误差反向传播算法为基础的多层前向神经网络。BP神经网络结构由输入层、隐藏层和输出层3部分组成(图5)。每一层都有若干个节点,前一层和后一层之间靠权值连接。信号的正向传播和误差的反向传播是BP神经网络的2个主要学习过程。BP网络通过不断学习调节,通常采用Sigmoid函数作为传递函数。
图4 热异常信号时空表达示意图Fig.4 The thermal anomaly signal space-time diagram.
图5 反向传播神经网络结构Fig.5 Back propagation neural network structure.
确定网络结构:网络为3层结构,输入层-隐层-输出层(图6)。确定各层神经元数目为d-N-m,其中d为输入神经元的特征数目,对应于表征热异常特征的9个指标,分别反映热异常的面积、强度、热异常起始位置、热异常起始时间和断裂带属性,共计12个输入神经元,即12维向量(每个位置信息均由2维向量刻画)。N为隐层神经元数目;m为输出神经元数目,对应于输出向量(以震中位置、震级、发震时间、发震标志为分量的5维向量)的维数。
本文采用100个5级以上震例和70个随机无震样本进行研究(其中80个震例和55个随机无震样本用于训练;20个震例和15个随机无震样本用于仿真)。训练调试表明,在隐层节点数为45时,误差达到0.001(图7),且收敛速度最快。
图6 BP神经网络结构示意图Fig.6 The schematic structure of BP neural network.
图7 误差下降曲线Fig.7 The error drop curve.
仿真结果的4种情形如下:1)对发生地震的样本预测出地震的发生,即预测准确的地震有16个,准确率为80%;2)对发生地震的样本未预测出发生地震的个数为4,漏警率是20%;3)对未发生地震的样本预测出地震的发生的个数为2,虚警率是13.33%;4)对未发生地震的样本未预测出发生地震的个数为13,准确率是86.67%。
在预测准确的16个地震中:一个震例的震级误差>3级,其余的都在3级以内;震级误差<2级的有11个,占68.75%;发震时间误差<1个月的有14个,占87.5%;震中位置误差在3°以内的有13个,占81.25%(表1)。
表1 震级、发震时间以及震中位置误差一览表Table 1 The error of magnitude,origin time and epicenter position of earthquake
本实验所用的震例中,5~6级震例较多,6~7级的震例较少,5~6级的地震和6~7级的地震预测试验的结果统计如表2。可以看出,对于不同震级的地震,5~6级预测精度在震级方面比6~7级地震的高,但在发震时间及震中位置的预测方面,其预测精度低于6~7级的地震。
表2 不同等级的地震预测的地震三要素误差对比表Table 2 The contrast table of errors of the three elements of earthquake prediction for different magnitude levels
隐层神经元的数量不同,神经网络结构不同,设置隐层神经元数目35到50进行了试验(表3),结果表明,当隐层神经元数目在45时,其预测的精度最高。因此,在未来实际地震预测中,隐层神经元数量的选择较为关键。
图8是实际震中与预测震中的空间位置对比一览图。图中,蓝色的点和红色的点分别表示实际的地震位置和预测的地震位置,分别对应蓝色和红色的数字注记。实际的地震位置和预测的地震位置在空间上是很接近的,总体上取得了较好的试验效果。
表3 不同隐层神经元的数量的神经网络预测效果Table 3 Neural network's predicting outcomes with different number of hidden layer neuron
本文结合神经网络的优点,提出将热异常信息结合断裂带信息作为地震预测的信息源,通过构建BP神经网络,进行地震三要素预测的新思路。本次地震预测试验是在以震中为中心的10°×10°矩形范围,以地震发生前2个月为时间范围来提取地震异常。对中国及周边100个5级以上震例,以及70个随机无震样本进行训练和仿真。试验结果表明,通过RST算法提取的震前热异常指数值,用于BP神经网络地震预测是可行的。但本试验是在已知震例情况下,试验区的中心位置即为震中,热异常研究的时间原点和时间范围均是在地震发生时间已知的前提下确定的,因此,试验结果(预测精度的高低)本质上是一种地震要素与热异常值间的非线性相关性的呈现,预测精度的高低反映了热异常与地震的相关性大小。由于在未来实际地震预测中,研究区位置并不是像本次试验这样,以已知震中位置去控制,而是需要专家知识结合活动断裂带来确定,因此未来实际的地震预测精度并不一定会达到本次试验的效果。试验中,研究区范围和神经网络结构(隐层神经元数量)对预测效果有较大的影响。
增加震例数量,强化神经网络训练。本次研究选用的是2001年以来的MODIS-LST遥感数据,样本数量对于神经网络预测来说还显不够,尤其是7级以上的大地震震例更少,今后的研究应将NOAA AVHRR遥感数据系列纳入进来,以增加震例数量,强化神经网络的训练。
关于热异常信息提取范围的选取问题。文中地震热异常提取与地震预测的研究范围是以震中为中心,按10°×10°的方形框,以地震发生前2个月为时间范围,使用已发生的震例对该方法进行预测及验证。对未来的地震预测,研究范围可参照构造断裂带及危险地震分区来进行框取,尝试地震的真正预测。
图8 实际震中和预测震中在断裂带中的位置Fig.8 The position of actual epicenter and predicted epicenter in the fault zone.
增加前兆输入信息和地震带属性信息。神经网络的输入仅考虑了断裂带等级和条数,将来的研究中,应适当增加某些重要的断裂带属性信息,如断裂带走向、断裂带倾向和力学性质等输入信息,探究其能否进一步提高预测精度。此外,热异常只是地震前兆之一,为加强地震预测的综合性,将来的研究应考虑辅之以更多的神经网络输入参数,如b值异常、地磁变化、电磁变化、地震活动周期以及累计频度等。