引 言

当前遥感技术已成为快速获取地震灾情信息及进行震害评估的有效手段,为震后应急救援提供了及时的信息保障(Hayashi et al,2000；Yusuf et al,2001；柳家航,2003；Matsuoka,Yamazaki,2004；王晓青等,2008；雷莉萍等,2010)．随着高分辨率商业卫星的迅猛发展以及米级和亚米级遥感影像的大量出现,能够利用遥感影像直接获取单体建筑物的破损程度,为灾害评估和应急救援提供了不可或缺的信息支持(张景发等,2002；曾涛等,2010；张小咏等,2013;董燕生等,2014)．然而,在追求获取高空间分辨率遥感图像的同时,也降低了获取灾区数据的时间效率,不能完成大震发生后覆盖震区的宏观遥感灾情的快速获取,从而也不能提供快速有效的应急救援决策．

相比米级-亚米级的高分辨率遥感影像,30—50 m中分辨率遥感影像一方面具有重访周期短、数据获取效率高和光谱信息多等优点,另一方面也能准确反映城市或居民区的整体震害情况,因此,中分辨率遥感影像也被广泛应用于地震灾情的快速评估中．目前,利用中分辨率遥感影像提取震害信息的方法主要有主成分分析法、光谱特征匹配法和变化监测法等． Estrada和Yamazaki(2000)利用Landsat影像数据对1999年土耳其以兹米特M_W7.4地震采用主成分分析法提取灾情信息,并与实际调查数据结果进行了对比分析；Matsuoka和Yamazaki(1998)利用Landsat和SPOT数据分析1995年日本阪神M_W6.9地震灾情信息,采用光谱特征匹配法从影像中提取了严重倒塌建筑物信息；Yusuf等(2001,2002)基于震前和震后的Landsat-7数据,采用变化检测法提取2001年印度M_W7.6地震的震害信息,并与高精度的SAR图像进行对比,验证了震害信息提取的精度；窦爱霞(2003)基于IRS-1C,ERS-2SAR和ETM+震前震后图像,采用变化检测法对1999年台湾集集M_W7.7地震和1998年张北M_L6.2地震进行震害评估．然而上述3种方法中,主成分分析法和光谱特征匹配法容易受到同谱异物的影响,变化检测法容易受到其它非震害因素的影响,这使得目前基于中分辨率遥感影像的震害信息提取精度相对较低．因此,充分发挥中分辨率遥感影像在地震灾情监测中费用低、数据处理效率高及灾情监测范围大等优点,探求具有更高精度的震害信息提取算法,为震后第一时间应急救援提供宏观震害信息是一项非常必要的工作．

本文针对中分辨率遥感影像建筑物震害信息弱以及变化检测法受非震害信息影响大等弱点,拟对中分辨率遥感影像的震害特征予以分析,在此基础上建立基于中分辨率遥感影像的居民区震害信息快速提取的流程和方法,并以2001年印度M_W7.6地震的极重灾区为研究区域,综合利用震前、震后中分辨率遥感影像和灯光影像,对该算法进行验证和分析．

1.   研究方法

1.1   震害特征分析

30—50 m中分辨率遥感影像无法对单个建筑物进行分析,但可以从居民区整体表现的辐射特征、光谱特征和空间特征来进行分析．辐射特征主要为居民区在单个波段上的灰度或辐射特征；光谱特征主要考察在多个波段上的整体反射特性；空间特征则主要表现为空间纹理特征．

图 1给出了2001年印度古吉拉特邦M_W7.6地震安贾尔震区震前和震后Landsat卫星TM遥感器第一波段的灰度值图像．图 2给出了该震区震前和震后图像的灰度值和相对反射率曲线．从图 1可以看出,地震受灾严重地区震后图像的灰度值较震前明显变亮,而图 2中的灰度值曲线也证明了这一点,其原因主要是地震引起大量建筑物倒塌,而倒塌的建筑物出现很多新的断面,从而导致辐射亮度变大．另外,对比图 2中的灰度值曲线与相对反射率曲线可以看出：震后图像波段1,2,3,4的灰度值变大,波段5和7的灰度值变小；而震后图像相对反射率在所有波段均变大．之所以出现该差异,是因为图像灰度值受到诸如太阳光照条件、遥感器参数设置、成像时间及成像时大气状况等多种因素的影响；而相对反射率在一定程度上消除了这些因素的影响,能更真实地反映地物的辐射和光谱特性．

图  1  印度安贾尔震区震前(a)和震后(b)TM影像(波段1)

Figure  1.  TM images (band 1) pre-earthquake (a) and post-earthquake (b) in Anjar seismic zone, India

下载: 全尺寸图片 幻灯片

图  2  印度安贾尔震区震前、 震后的灰度值(a)和相对反射率(b)变化

Figure  2.  Variation of gray value (a) and relative reflectance (b) pre- and post-earthquake in Anjar seismic zone, India

下载: 全尺寸图片 幻灯片

图 3给出了2010年海地M_W7.0地震太子港地区Landsat卫星TM遥感影像震前和震后的灰度值和相对反射率曲线．可以看出,灰度值和相对反射率曲线在地震前后具有相似的特征,而2001年印度M_W7.6地震反射率相对灰度值更能反映倒塌建筑物反射率增加的客观事实．从上述两个震区的影像特征分析来看,震后居民区图像的反射率相对灰度值明显变大．

图  3  海地震区震前、 震后的灰度值(a)和相对反射率(b)变化

Figure  3.  Variation of gray value (a) and relative reflectance (b) pre- and post-earthquake in Haiti seismic zone

下载: 全尺寸图片 幻灯片

对于居民区在中分辨率遥感影像上的空间纹理特征,本文采用常用的灰度共生矩阵二阶概率统计的特征参数来分析,包括均值、方差、协同性、对比度、相异性、信息熵、二阶矩和相关性．图 4分别给出了印度安贾尔震区和海地震区震前和震后遥感影像纹理参数的变化情况．可以看出：对于印度地震,居民区震后均值、方差和对比度较震前变化较大,其它几个参数则变化不大；对于海地地震,居民区震后较震前仅有对比度变化较大,方差则变化不大．这说明除了对比度之外,其它纹理参数均不适合作为地震灾害的表征参数．纹理特征之所以不能很好地反映震区的破坏信息,是因为30—50 m的空间分辨率(单个像元面积为900—2500 m²)相对建筑物个体来说太低,建筑物几乎是与背景混合在一起,因此建筑物的破坏对混合像元的纹理特征影响不太明显．

图  4  印度安贾尔震区(a)和海地震区(b)震前、 震后遥感影像的纹理参数变化

Figure  4.  Texture parameter variation pre- and post-earthquake in Anjar, India (a) and Haiti seismic zones (b)

下载: 全尺寸图片 幻灯片

1.2   震害信息提取方法与流程

基于上述对居民区震害在遥感图像上的特征分析得到以下结论： ①地震使居民区建筑物产生破坏,破坏后的建筑物新鲜断面引起震后居民区影像亮度增加；②反射率图像由于消除了太阳光照条件、遥感器参数设置和大气状况等多种因素的影响,所以其相对图像灰度值更能体现建筑物的破坏特征；③ 30—50 m中分辨率图像相对建筑物分辨率太低,其空间纹理特征不能很好地表征震害特征．总之,对于中低分辨率遥感影像,震后居民区反射率值增大,因此可以利用震后与震前反射率差值来提取居民区震害信息．然而,反射率差值法提取的震害信息也会包含反射率变化的其它错误信息,如云、植被和水体等．因此,必须对背景信息进行剔除,得到最终的有效震害信息．同时,居民区震后图像的各波段反射率值相对震前均增大,且增加的幅度也基本一致(图 2,3),即各波段间存在很好的相关性,因此可以通过融合各个波段的光谱变化特征,把相关波段的信息集中于第一主成分,从而能更加有效地提取震害信息．

基于居民区震害特征分析和去伪存真的思路,本文建立了以反射率差值为基础的中分辨率遥感影像居民区震害信息的提取方法,主要包括以下5个步骤： ①获取地震区震前及震后的遥感数据,并利用监督分类法提取居民区,考虑到分类误差导致提取的结果中除了居民区之外还有与居民区相似的其它地物,将利用监督分类法提取的居民区结果称之为似居民区；②对震前和震后遥感影像进行空间几何配准和反射率反演等预处理,其中反射率反演可以采用内部平均法(童庆禧等,2006)计算图像的相对反射率,相对反射率值即为图像每个像素的灰度值与整景图像灰度值均值的比值；③计算地震后与地震前反射率差值图像,并对反射率差值图像进行主成分变换；④提取主成分变换后的第一主成分,利用阈值分割初步提取震害区；⑤利用第一步提取的居民区和灯光影像数据对震区结果进行两次优化,得到最终的居民区震害信息.该方法提取流程如图 5所示．

图  5  居民区震害信息提取流程图

Figure  5.  Flow chart for extracting residential seismic damage information

下载: 全尺寸图片 幻灯片

2.   印度古吉拉特邦地震试验与结果

2.1   数据

北京时间2001年1月26日11时16分印度古吉拉特邦发生M_W7.6地震,震中位于(23.2°N,70°E)．该地震造成2万人死亡,6万多人受伤,20万人无家可归,给古吉拉特邦地区的人民带来灾难性后果,受此次地震影响最大的城市有帕焦、普杰、安贾尔和甘地特姆等(Narayan et al,2002)．本文以安贾尔市附近的Landsat卫星TM图像为例,提取该城市的震害信息． TM数据来源于美国USGS网站,图像空间分辨率为30 m,已经辐射校正和正射纠正．印度M_W7.6地震震前TM影像(图 6a)于2001年1月8日获取,图像整体质量较好,左上部有少量薄云和阴影；震后TM影像(图 6b)于2001年2月9日获取,图像质量好于震前．另外,从美国NOAA网站获取了该地区2001年常年稳定灯光的平均数据,其空间分辨率为1 km．

图  6  印度古吉拉特邦震区震前(a)和震后(b)TM影像

Figure  6.  TM images pre-earthquake (a) and post-earthquake (b) in Gujarat seismic zone, India

下载: 全尺寸图片 幻灯片

2.2   震害信息提取结果

由于TM数据和灯光数据均经过辐射和几何预处理,且图像之间的几何配准精度高,无需再进行校正和配准处理．根据本文算法,首先,对TM影像采用内部平均法进行相对反射率反演；其次,以人工选择的居民区为参考样本,利用光谱夹角匹配法(童庆禧等,2006)提取似居民区专题图；然后,根据是否有灯光数值对图像进行二值化,制作居民区掩膜,用于后期对居民区的优化处理；最后,根据本文计算流程依次进行反射率差值、主成分变换、震害信息初步提取和震害信息优化等处理后,得到最终的居民区震害信息．

图 7a为基于震前与震后图像反射率差值的主成分图像初步提取的震害结果．由于震后图像受薄云阴影(图中红色箭头所指区域)、季节的差异及土壤水分含量等因素的影响,薄云阴影、山区大片植被和海边沙滩被错误地提取出来．图 7b为利用监督分类法提取的似居民区对图 7a进行优化后的结果,可以看出,大部分非居民区的错误目标基本都被去除掉,但是左上角有阴影的无植被覆盖的裸土地仍然存在．图 7c为利用灯光遥感数据对图 7b进行优化的最终结果,可以看出,仅剩下几个居民区内有建筑物震害信息(图中红色圆圈内区域)．图 8为图 7c中红色圆圈内的放大图,从左至右分别是安贾尔市震前、震后以及震害信息叠加在震后影像上的效果．对比图 8震前与震后的图像可以看出,本文方法提取的结果与安贾尔市反射率变大部分非常一致．

图  7  印度古吉拉特邦震区地震前后反射率图像第一主成分差值图像(a)、 利用似居民区数据 对图7a第一次优化后结果(b)以及利用灯光数据第二次优化后结果(c)

红色箭头所指区域为受薄云影响区域， 红色圆圈内为提取的安贾尔市震害信息

Figure  7.  First principal component D-value image about reflectance in Gujarat seismic zone, India (a), the first optimization results using similar residential area extraction algorithm (b), and the second optimization results using nighttime lights data (c). Red arrow refers to the area affected by the cloud area, red circle delineates the earthquake damaged information in Anjar city

下载: 全尺寸图片 幻灯片

图  8  安贾尔市震前(a)、 震后(b)TM图像及震害信息提取结果(c)， 白色斑点区域即为震害区域

Figure  8.  TM images pre-earthquake (a) and post-earthquake (b) in Anjar city, and the extraction results of seismic damage information (c), where the white spots are the seismic damage zone

下载: 全尺寸图片 幻灯片

3.   讨论与结论

基于地震前后图像变化检测法提取震害信息不可避免地会受到多种因素的影响,本文首先利用主成分变换对原始图像的震害信息进行了增强,然后利用居民区数据和灯光数据等对非震害信息进行了剔除,从而较好地解决了利用中分辨率遥感影像提取居民区震害信息时目标信息弱、背景信息复杂等问题．

图 9a,b给出了印度安贾尔市震前与震后相对反射率图像直接差值运算后的差值图像直方图,图 9c给出了对差值图像进行主成分变换后的主成分图像直方图．可以看出,主成分变换前所有波段的直方图具有相似的正态分布,且直方图峰值均大于0,这说明震后所有波段的反射率均大于震前,该现象与图 2,3中相对反射率差值曲线一致．但是经主成分变换后的直方图则发生了明显变化：除第一主成分之外的其它主成分分量直方图均具有相似的正态分布,但峰值均降低为0,表明这几个主成分分量基本不包含震害信息；第一主成分明显与其它分量分离开,且第一主成分的峰值相对主成分变换前明显变大,说明主成分变换能将平均分配在各个波段的震害信息集中到第一主成分量,因此利用第一主成分分量更容易提取震害信息．图 10给出了基于相对反射率差值图像提取的震害信息和基于差值图像主成分变换后的第一主成分提取的震害信息.可以看出,经过主成分变换后的结果明显优于变换之前的结果．

图  9  印度安贾尔市震前震后相对反射率差值图波段1(a)和波段2—4的直方图分布(b) 以及差值图像经过主成分变换后的直方图分布(c)

Figure  9.  Histogram distribution of band 1 (a) and band 2--4 (b) of relative reflectance D-value pre-earthquake and post-earthquake and that under principal component transformation (c) in Anjar city, India

下载: 全尺寸图片 幻灯片

图  10  基于相对反射率差值(a)和基于该差值主成分变换后(b)的印度安贾尔震区震害信息提取结果

Figure  10.  Extraction result of seismic damage information in Anjar seismic zone, India based on relative reflectance D-value (a) and that under principal component transformation (b)

下载: 全尺寸图片 幻灯片

许多非地震因素也会在地震前后差值图像上产生反射率差异,使得利用差值算法提取的震害信息出现错误结果,如图 7a中云阴影、植被和沙滩等都是容易产生误差的地物．解决该问题的一种思路是采用间接法,即将研究区内各种背景对象依次剔除；另一种思路是采用直接法,直接从图像中提取居民区进而提取震害区域．理论分析和试验结果表明,间接法需要针对不同的背景对象设计相应的算法,由于背景对象复杂导致整体算法复杂,且在每种背景对象消除过程中所产生的误差均会叠加到最终结果中．相比而言,直接法的算法和处理过程则相对简单,只有一次误差累积过程；但直接法会受到“同谱异物”的影响,影像中与居民区具有相似特征的背景地物会被当作居民区而被提取(图 7a)．为解决该问题,本文引入灯光遥感数据,采用直接法提取居民区．众所周知,居民区与裸土、沙滩、河床等地物最大的差别在于居民区有人居住,而灯光是人类活动的典型标志之一,因此利用灯光遥感影像可以将居民区与裸土等背景地物区分开．具体作法分两步： ①基于震前遥感数据,以居民区为参照,采用监督分类法提取似居民区,由于植被、水体、云等与居民区在统计特征上差异较大,所以似居民区内几乎不包括这类背景,仅包含与居民区相似的裸土区或其它类似区域；②利用灯光图像对似居民区进行二次优化,去除无人居住的裸土地等背景目标．从图 7的优化结果看,利用似居民区能有效去除植被、水体、云等背景误差,然后经过灯光遥感数据的再次优化,居民区之外的误差几乎全部被消除．

本文利用中分辨率遥感影像进行地震灾情信息监测,具有费用低、数据获取机会多、数据处理效率高及灾情监测范围大等优点,对震后应急救援阶段快速提供宏观震害信息具有非常重要的意义．通过对中分辨率遥感影像震害特征的分析表明,在30—50 m中分辨率遥感影像上,震后图像反射率变大是以建筑物为主的居民区最显著的震害特征,因此本文建立了一种基于反射率变化检测的方法来提取居民区震害信息．鉴于震后反射率相对震前反射率在遥感影像所有波段均具有相似的特征,对震后-震前反射率差值采用主成分变换,将均匀分布在各个波段的震害信息集中在第一主成分,从而达到增强震害信息的目的．另外,采用监督分类法从震前影像中提取似居民区,并利用灯光遥感数据提取居民区范围,通过这两种方法有效剔除了居民区的非震害信息,从而显著提高了直接基于反射率变化提取震害信息的精度．