基于图像识别的震象云地震预测方法

谢 庭1a,1b,1c，陈 忠1a,1b,1c，李志平2，张宁新1a,1b,1c，郭莉莉1a,1b,1c

(1.华中科技大学 a.自动化学院；b.多谱信息处理技术国家级重点实验室；c.图像信息处理与智能控制教育部重点实验室，武汉 430074；2.大港油田集团通信公司，天津 300280)

1 概述

到目前为止，地震预报，尤其是短临地震预报，是一个世界性难题。地震预报必须同时包括时间、地点和强度。由于地震情况复杂，有些地震能预报，有些则无法预报，现在全球预报地震的准确率只有20%多。然而地震是极具破坏性的自然灾害，给人类带来了巨大的物质和精神损失，因此探索和发现预报地震的新方法、新途径，在防震减灾方面具有十分深远的意义。

随着高新技术的不断发展，卫星热红外遥感为地震预测预报提供了新的方法和手段。早在1971年，我国著名地震学家傅承义提出红肿学说，为后来利用卫星遥感、监测地表与低空大气的异常增温预报地震奠定了理论基础[1]。1988年，前苏联科学家Gorny、Salman和Tronin在研究中亚以及东地中海地区的地震时就发现在中强地震震前，震区及其周边区域对应的卫星红外图像上存在热红外异常[2]。意识到这一发现的重要性，我国一大批学者投入到了卫星热红外异常预测地震的研究中。例如，强祖基、赁常恭等利用卫星遥感热红外图像，经过多年对震前热红外异常进行分析和研究，取得了显著的地震预报成果[3-5]；徐秀登、徐向民通过分析反生在中国及相邻地区的40多次地震的热红外图像，总结了卫星红外临震异常的基本特征[6-7]；刘德富等研究了地球长波辐射时空变化，结果表明在地震发生前的震区呈现出比周围区域更显著的辐射增强变化特征[8]。还有其他很多学者也参与到了卫星热红外异常预测地震的研究中，并取得了一定的成果。在上述传统方法研究过程中，红外异常均由人眼识别或手工计算得出，而当今时代人工智能发展日趋成熟，利用图像识别技术提取红外异常是这一研究领域发展的必然趋势。

本文提出一种基于图像识别热红外遥感异常数据的自动地震预测方法。不仅为基于热红外遥感数据的地震预测研究领域打开了人工智能的大门，而且实现国产风云2号气象卫星数据在地震预测中的应用。传统方法使用的热红外数据取自国美NOAA卫星或日本MTSAT卫星，而本文方法使用国产风云2号气象卫星(FY_2C)1通道分区图。该图一般应用于天气预报，也是首次将其应用于地震预测预报中。

2 震象云地震预测原理及图像识别方法描述

2.1 震象云地震预测原理

通过对卫星云图的观察，发现地震前震中上方的卫星红外云图会出现增温型特殊云团和降温型特殊云团，一般接近发震时间时出现较频繁，称这种特殊云团为震象云。

震象云是地震前震区上空出现增温和降温异常的特殊云迹。图1所示是在汶川地震震前成都上空出现的震象云。

图1 汶川地震震前汶川附近出现的震象云

震象云与气象中的云团是有差别的，震象云一般位于对流层顶，与周围气象云团有温度差异，有一定的纹理特征，并且这些震象云移动速度要远小于周围气象云团的移动速度，且接近发震时间时出现较为频繁。因此，可以综合震象云的颜色、纹理以及浮现频率来预测地震。

2.2 图像识别方法描述

从红外云图中观察到增温型震象云浮现频率相对较为频繁，颜色、纹理特征较降温型震象云更为突出，因此，本文针对增温型震象云进行自动检测和识别。

红外云图中增温型的震象云表现为在高层云中出现中、低层云，灰度级集中出现在110～145之间，云团大小在一定的范围变化，形状不确定，有一定模式的纹理特征，而且震象云是在地震地点上空附近浮动，可以根据以上特征对其进行目标识别和跟踪。针对震象云的以上特点，对其采取了以下方案进行处理。首先对输入数据进行预处理，拉伸震象云所在目标灰度区段的灰度级，压缩非目标灰度区段的灰度级；利用灰度共生矩阵对热红外数据进行纹理特征提取，使用BP神经网络模型训练目标神经网络，将纹理特征输入目标神经网络进行识别，提取疑似目标；根据目标的大小范围滤掉非目标；对剩下的疑似目标进行跟踪，依据其浮现频率精确定位目标，提取得到震象云。主要分为4个阶段：(1)预处理阶段；(2)目标识别阶段；(3)目标过滤阶段；(4)目标跟踪阶段。处理流程如图2所示。

图2 震象云提取整体流程

2.2.1 预处理

根据观测，增温型震象云主要位于是中层云区段，在红外云图中灰度级集中于110～145之间。为了在后续流程中能够充分提取增温型震象云的纹理特征，必须扩大中层云灰度级的浮动范围。因此，要压缩低层云和高层云，拉伸中层云。本文采用分段线性拉伸方法增强震象云区域对比度，采用的公式如下所示：

其中，r是输入灰度级；s为输出灰度级。

2.2.2 目标识别

由于BP网络能学习和存储大量的输入-输出模式映射关系，而无需事前揭示描述这种映射关系的数学方程，并且具有相当长的研究历史，发展成熟。因此，本文利用BP神经网络模型训练目标网络，进行目标识别。

识别过程为：对某一像素点的n×n邻域进行纹理特征提取，将特征向量输入到目标神经网络中进行识别，输出该像素点的特性判别值，若特性判别值落在接受域(ac1,ac2)之中，则断定为疑似目标点(由于震象云的纹理特征不具备100%的识别力，因此用其特征进行区别时，会将一些非震象云区域的像素点也包含进来)，否则为非目标点。经实验表明，邻域n×n中取n=51时较好，因为取样时，发现在51× 51正方形区域内可以较好地包含震象云，这样可以尽量获取震象云的纹理特征。

此过程的算法流程如下：

Step1特征提取。扫描以某一像素点为中心51× 51邻域，统计目标灰度级区段(拉伸前是110～145，拉伸后是55～200)0°、45°、90°、135°等4 个方向上的灰度共生矩阵，然后分别计算每个矩阵的对比度、相关性、能量和同质性这4个特性的值。这样将获得16个纹理特征值组成的特征向量。主要计算公式如下：

上述公式中用到的频度 p(i,j/d,θ) 、期望(ux,uy)、方差(σx,σy)等计算公式及各元素含义可参考文献[9-10]。

Step2获取目标神经网络：应用BP神经网络学习模型训练目标网络。输入为震象云和非震象云的特征向量，特征向量按照Step1进行计算。期望输出为特性判别值，震象云的为1.0，非震象云的为–1.0。按照BP神经网络模型理论[11-12]，设置训练模型的相关参数，进行训练，最终输出目标神经网络。

Step3利用目标神经网络对输入图像进行目标识别和分类。扫描图片，对于灰度级落在目标灰度区域内的像素点统计其邻域纹理特征，输入到目标神经网络进行识别，并且将特性判别值落在接受域内的标记为疑似目标点，否则标记为非目标点；对于落在目标灰度区域之外的像素点，直接标记为非目标点，不需要进行识别。

同一片震象云中的像素点虽然纹理是相似的，但并非一模一样，因此，计算出的特征向量会有差异，通过神经网络映射后的特性判别值不一定会与正样本的特性判别值相等，而是会在其值附近波动。上文实验规定震象云样本的特性判别值为1.0，实验验证目标神经网络输出值接受域为(ac1,a c2)=(0.5,1.2)较好，这样可以完整包含震象云区域的像素点。

2.2.3 目标过滤

对整张图片的每个像素点进行处理之后，疑似目标点会聚积形成疑似目标区域。目标过滤就是过滤掉疑似目标区域中过大或过小的区域，其中过大的区域是背景的云层图像，过小的区域是图像中的噪声或云团运动过程中拉伸出来的小空洞。具体实现为：统计每个8连通疑似目标区域的像素点个数Area(即面积)，若Area的取值落在接受域(a c3,a c4)之外，则标记为非目标块，过滤掉；否则保留疑似目标。实验证明，疑似目标块面积Area接受域取值为(a c3,a c4)=(250,1 2 50)较为理想，既可以滤掉大的背景云团，又可滤掉小的噪声或空洞。

2.2.4 目标跟踪

目标跟踪是对多张图片的疑似目标块进行处理的操作，分为目标标准化和跟踪统计2个步骤：

(1)目标标准化：扫描每一个目标过滤剩下的疑似目标块，计算其几何中心(即点集在水平和竖直方向上的平均值)和半径radious(公式为radious=)，即用圆形拟合跟踪区域)。

(2)跟踪统计：这个步骤需要建立疑似目标集，对每个出现出新的疑似目标进行判断是否属于已经存在的目标集。判断依据：如果疑似目标中心点和集合中心点的距离不超过前集合半径和疑似目标半径中的最大值，则该疑似目标属于该集合。如果属于某一个目标集，则将该疑似目标作为新元素加入到该集合中，调整这个疑似目标集的中心(即计算集合中所有元素中心的平均值)和半径(取当前集合半径和加入元素半径中的最大值)；如果不属于任何一个已存在的集合，则新建一个集合，将该疑似目标加入进去，将集合的中心和半径初始化为疑似目标的中心和半径。

当集合中的元素(疑似目标)个数超过一定的阈值5时，可以报告这个范围出现的疑似目标为目标，根据目标出现的周期长短，预测地震的发生；否则，继续跟踪，若出现次数太少(如只有一两次)且长期未出现，则可以删除这个集合。这是因为目标(增温型震象云)区域所对应的地面断断续续释放能量，为其加热，故其位置可以保留；而非目标区域为高层云运动时云团内部速度不一致，导致撕裂偶然形成的，没有能源支撑，故其不能保留。

3 实验及结果分析

3.1 实验

本文以汶川地震为例进行实验。实验对象为风云2C卫星的长波红外1通道云图，云图格式为.AWX，波长范围为10.3～11.3。

以FY2C_ANI_IR1_R01_20071014_0800为例说明单张图片处理过程，该图像大小为1200×1200。首先根据云图格式[13]读入图像数据，进行预处理，提取特征输入到目标神经网络进行识别。

图3为长波红外云图的原始数据所呈现的图像。图4为灰度经过式(1)分段线性拉伸处理之后的图像，从这幅图像中可以明显地看到在图像中央的白色云团中有一个长条形的孤立灰黑色的云斑，这个云斑就是震象云，即图1中所标识区域。图5为经过目标神经网络识别后留下的疑似目标所形成的图像，从图像中可以看出图4中的震象云区域得到完整的保留，其中白色点集表示疑似目标。图6是对图5进行目标过滤之后剩下的疑似目标所形成的图像。从这一系列图中看出，经过一系列处理之后，疑似目标区域只剩下少数几个，为后续的跟踪奠定了基础。

图3 长波红外原始图像

图4 灰度分段线性拉伸

图5 目标神经网络识别

图6 目标过滤输出

目标跟踪：为了简化跟踪，先将目标标准化，再在标准化的目标区域进行跟踪，统计同一区域疑似目标出现的次数，最后输出结果。图7为图6目标标准化后的目标分布图，图8为从2007年10月到2008年4月底跟踪一系列晚间云图所得到的疑似云目标分布。

图7 单张图片的标准化目标

图8 多张图片的标准化目标

跟踪晚间云图的原因是晚间没有阳光辐射且人类活动也相对白天较少，能够较好地反映地壳释放能量对大气增温的影响，有利于震象云跟踪的准确性。图9(a)是进行疑似目标跟踪后，同一区域出现疑似目标的次数和分布的示意图，同一矩形区域内的点表示该区域出现过的疑似目标。当目标出现次数达到一定的值时，这些疑似目标就很有可能是同一个地震带活动释放能量所产生的震象云。图9(b)为出现2次或2次以上疑似目标区域的分布示意图，图9(c)为出现3次或3次以上疑似目标区域的分布示意图，图9(d)为出现8次或8次以上疑似目标区域的分布示意图。

图9 跟踪结果

在实验中进行目标跟踪后，得到同一区域疑似目标出现次数与区域个数的关系如下：只出现疑似目标1次的区域有84个，出现疑似目标2次的区域有9个，出现疑似目标8次的区域只有1个，没有区域出现3次～7次疑似目标。本实验最终结果是将出现8次的疑似目标断定为目标，即将在同一区域出现的8次疑似目标确定为震象云。这些震象云的几何中心均包含在图片中以(688,690)为中心，半径是19单位像素的圆形区域内。图10为本实验最终跟踪到的目标，按照时间序列由远及近8次出现时的区域形态图，即汶川地震震前的震象云所在区域图。

图10 汶川地震震前震象云出现区域

3.2 结果分析

根据各图震象云出现的位置，结合当时风级、风向对震象云带来的偏移效果，将图像坐标转化成地理坐标，预测震中坐标是北纬30.885°± 0.1°，东经103.632°± 0.1°，而汶川地震震中是北纬30.986°，东经103.364°，虽没有精准地与实际震中重合，但偏移范围较小。

从实验结果看，汶川地震发生前，其上空至少出现了8次的增温型震象云，第一次检测到该区域的震象云是2007年10月14号(云图名为FY2C_ANI_IR1_R01_20071014_0800)，最后一次是2008年4月13日(云图名为FY2C_ANI_IR2_R01_20080413_1000)，距离地震发生时间5月12日是29天，属于短期地震预测。

4 结束语

本文利用图像识别技术提取震象云，实现了地震预测的自动化，因为基于震象云的地震预测本身具有震中位置预测较准确的特点，这就弥补了传统热红外地震预测方法中的主要缺陷。

由于震象云形成过程非常复杂：图像上的震象云云团体现的是对应位置对流层云顶的温度分布；这个温度分布不仅与地壳活动强烈和释放的各物质成分有关，还和当时周围的大气坏境有关。下一步的研究方向是单纯从震象云的形态上准确预测地震震级。

[1]傅承义.地球十讲[M].北京:科学出版社,1976.

[2]Gorny V I,Salman A G,Tronin AA,et al.The Earth Outgoing IR Radiation as an Indicator of Seismic Activity[J].Proceedings of the USSR Academy of Sciences,1988,30(1):67-69.

[3]强祖基,孔令昌,王戈平,等.地球放气，热红外异常与地震活动[J].科学通报,1992,37(24):2259-2262.

[4]强祖基,孔令昌,郭满红,等.卫星热红外增温机制的实验研究[J].地震学报,1997,19(2):197-201.

[5]强祖基,赁常恭,李玲芝,等.卫星热红外图像亮温异常-短临震兆[J].中国科学:D辑,1998,28(6):564-573.

[6]徐秀登,张行才,李贵达.张北地震与大气增温异常[J].西北地震学报,2000,22(1):44-47.

[7]徐秀登,徐向民.地震前红外异常的基本特征与成因机制[J].西北地震学报,2001,23(3):310-312.

[8]刘德富,康春丽.地球长波辐射(OLR)遥感与重大自然灾害预测[J].地学前缘,2003,10(2):427-435.

[9]侯海苗,冀小平.基于灰度共生矩阵的纹理特征[J].长治学院学报,2008,25(5):31-32.

[10]宁顺刚,白万民,喻 钧.基于灰度共生矩阵的图像分割方法研究[J].电子科技,2009,22(11):69-71.

[11]樊振宇.BP神经网络模型与学习算法[J].软件导刊,2011,10(7):66-68.

[12]李秀珍,孔纪名,李朝凤.基于Matlab的BP神经网络在泥石流危险性评价中的应用[J].工程勘察,2010,38(1):47-50.

[13]许健明,张文建,杨 军,等.风云二号卫星业务产品与卫星数据格式实用手册[M].北京:气象出版社,2008:229-232.