张小咏 李庆亭 高 娜
1) 中国北京100049中国地震应急搜救中心 2) 中国北京100094中国科学院遥感与数字地球研究所
基于中分辨率遥感影像的居民区震害信息提取*
1) 中国北京100049中国地震应急搜救中心 2) 中国北京100094中国科学院遥感与数字地球研究所
针对中分辨率遥感影像建筑物震害信息弱以及变化检测法受非震害信息影响大等弱点, 本文建立了一种基于变化检测的居民区震害信息快速提取方法. 该方法利用主成分变换增强震害信息, 采用监督分类法提取似居民区, 并用灯光影像数据进一步对似居民区提取结果进行优化, 从而很好地消除了变化检测方法中非震害因素的影响. 在此基础上, 以2001年印度MW7.6地震的极重灾区为研究区域, 利用震前、 震后Landsat卫星TM图像和震区灯光影像数据, 对本文算法进行了验证和分析. 结果表明, 在30—50 m中分辨率遥感影像上, 以建筑物为主的居民区震后图像变化最为显著的震害特征是反射率变大, 本文所建立的居民区震害信息提取方法在解决中分辨率遥感影像震害目标信息弱、 背景复杂等方面效果明显.
震害信息 居民区 变化检测 中分辨率遥感影像 灯光遥感数据
当前遥感技术已成为快速获取地震灾情信息及进行震害评估的有效手段, 为震后应急救援提供了及时的信息保障(Hayashietal, 2000; Yusufetal, 2001; 柳家航, 2003; Matsuoka, Yamazaki, 2004; 王晓青等, 2008; 雷莉萍等, 2010). 随着高分辨率商业卫星的迅猛发展以及米级和亚米级遥感影像的大量出现, 能够利用遥感影像直接获取单体建筑物的破损程度, 为灾害评估和应急救援提供了不可或缺的信息支持(张景发等, 2002; 曾涛等, 2010; 张小咏等, 2013; 董燕生等, 2014). 然而, 在追求获取高空间分辨率遥感图像的同时, 也降低了获取灾区数据的时间效率, 不能完成大震发生后覆盖震区的宏观遥感灾情的快速获取, 从而也不能提供快速有效的应急救援决策.
相比米级-亚米级的高分辨率遥感影像, 30—50 m中分辨率遥感影像一方面具有重访周期短、 数据获取效率高和光谱信息多等优点, 另一方面也能准确反映城市或居民区的整体震害情况, 因此, 中分辨率遥感影像也被广泛应用于地震灾情的快速评估中. 目前, 利用中分辨率遥感影像提取震害信息的方法主要有主成分分析法、 光谱特征匹配法和变化监测法等. Estrada和Yamazaki(2000)利用Landsat影像数据对1999年土耳其以兹米特MW7.4地震采用主成分分析法提取灾情信息, 并与实际调查数据结果进行了对比分析; Matsuoka和Yamazaki(1998)利用Landsat和SPOT数据分析1995年日本阪神MW6.9地震灾情信息, 采用光谱特征匹配法从影像中提取了严重倒塌建筑物信息; Yusuf等(2001, 2002)基于震前和震后的Landsat-7数据, 采用变化检测法提取2001年印度MW7.6地震的震害信息, 并与高精度的SAR图像进行对比, 验证了震害信息提取的精度; 窦爱霞(2003)基于IRS-1C, ERS-2SAR和ETM+震前震后图像, 采用变化检测法对1999年台湾集集MW7.7地震和1998年张北ML6.2地震进行震害评估. 然而上述3种方法中, 主成分分析法和光谱特征匹配法容易受到同谱异物的影响, 变化检测法容易受到其它非震害因素的影响, 这使得目前基于中分辨率遥感影像的震害信息提取精度相对较低. 因此, 充分发挥中分辨率遥感影像在地震灾情监测中费用低、 数据处理效率高及灾情监测范围大等优点, 探求具有更高精度的震害信息提取算法, 为震后第一时间应急救援提供宏观震害信息是一项非常必要的工作.
本文针对中分辨率遥感影像建筑物震害信息弱以及变化检测法受非震害信息影响大等弱点, 拟对中分辨率遥感影像的震害特征予以分析, 在此基础上建立基于中分辨率遥感影像的居民区震害信息快速提取的流程和方法, 并以2001年印度MW7.6 地震的极重灾区为研究区域, 综合利用震前、 震后中分辨率遥感影像和灯光影像, 对该算法进行验证和分析.
1.1 震害特征分析
30—50 m中分辨率遥感影像无法对单个建筑物进行分析, 但可以从居民区整体表现的辐射特征、 光谱特征和空间特征来进行分析. 辐射特征主要为居民区在单个波段上的灰度或辐射特征; 光谱特征主要考察在多个波段上的整体反射特性; 空间特征则主要表现为空间纹理特征.
图1给出了2001年印度古吉拉特邦MW7.6地震安贾尔震区震前和震后Landsat卫星TM遥感器第一波段的灰度值图像. 图2给出了该震区震前和震后图像的灰度值和相对反射率曲线. 从图1可以看出, 地震受灾严重地区震后图像的灰度值较震前明显变亮, 而图2中的灰度值曲线也证明了这一点, 其原因主要是地震引起大量建筑物倒塌, 而倒塌的建筑物出现很多新的断面, 从而导致辐射亮度变大. 另外, 对比图2中的灰度值曲线与相对反
图1 印度安贾尔震区震前(a)和震后(b)TM影像(波段1)
图2 印度安贾尔震区震前、 震后的灰度值(a)和相对反射率(b)变化
射率曲线可以看出: 震后图像波段1, 2, 3, 4的灰度值变大, 波段5和7的灰度值变小; 而震后图像相对反射率在所有波段均变大. 之所以出现该差异, 是因为图像灰度值受到诸如太阳光照条件、 遥感器参数设置、 成像时间及成像时大气状况等多种因素的影响; 而相对反射率在一定程度上消除了这些因素的影响, 能更真实地反映地物的辐射和光谱特性.
图3给出了2010年海地MW7.0地震太子港地区Landsat卫星TM遥感影像震前和震后的灰度值和相对反射率曲线. 可以看出, 灰度值和相对反射率曲线在地震前后具有相似的特征, 而2001年印度MW7.6地震反射率相对灰度值更能反映倒塌建筑物反射率增加的客观事实. 从上述两个震区的影像特征分析来看, 震后居民区图像的反射率相对灰度值明显变大.
图3 海地震区震前、 震后的灰度值(a)和相对反射率(b)变化
对于居民区在中分辨率遥感影像上的空间纹理特征, 本文采用常用的灰度共生矩阵二阶概率统计的特征参数来分析, 包括均值、 方差、 协同性、 对比度、 相异性、 信息熵、 二阶矩和相关性. 图4分别给出了印度安贾尔震区和海地震区震前和震后遥感影像纹理参数的变化情况. 可以看出: 对于印度地震, 居民区震后均值、 方差和对比度较震前变化较大, 其它几个参数则变化不大; 对于海地地震, 居民区震后较震前仅有对比度变化较大, 方差则变化不大. 这说明除了对比度之外, 其它纹理参数均不适合作为地震灾害的表征参数. 纹理特征之所以不能很好地反映震区的破坏信息, 是因为30—50 m的空间分辨率(单个像元面积为900—2500 m2)相对建筑物个体来说太低, 建筑物几乎是与背景混合在一起, 因此建筑物的破坏对混合像元的纹理特征影响不太明显.
图4 印度安贾尔震区(a)和海地震区(b)震前、 震后遥感影像的纹理参数变化
1.2 震害信息提取方法与流程
基于上述对居民区震害在遥感图像上的特征分析得到以下结论: ① 地震使居民区建筑物产生破坏, 破坏后的建筑物新鲜断面引起震后居民区影像亮度增加; ② 反射率图像由于消除了太阳光照条件、 遥感器参数设置和大气状况等多种因素的影响, 所以其相对图像灰度值更能体现建筑物的破坏特征; ③ 30—50 m中分辨率图像相对建筑物分辨率太低, 其空间纹理特征不能很好地表征震害特征. 总之, 对于中低分辨率遥感影像, 震后居民区反射率值增大, 因此可以利用震后与震前反射率差值来提取居民区震害信息. 然而, 反射率差值法提取的震害信息也会包含反射率变化的其它错误信息, 如云、 植被和水体等. 因此, 必须对背景信息进行剔除, 得到最终的有效震害信息. 同时, 居民区震后图像的各波段反射率值相对震前均增大, 且增加的幅度也基本一致(图2, 3), 即各波段间存在很好的相关性, 因此可以通过融合各个波段的光谱变化特征, 把相关波段的信息集中于第一主成分, 从而能更加有效地提取震害信息.
基于居民区震害特征分析和去伪存真的思路, 本文建立了以反射率差值为基础的中分辨率遥感影像居民区震害信息的提取方法, 主要包括以下5个步骤: ① 获取地震区震前及震后的遥感数据, 并利用监督分类法提取居民区, 考虑到分类误差导致提取的结果中除了居民区之外还有与居民区相似的其它地物, 将利用监督分类法提取的居民区结果称之为似居民区; ② 对震前和震后遥感影像进行空间几何配准和反射率反演等预处理, 其中反射率反演可以采用内部平均法(童庆禧等, 2006)计算图像的相对反射率, 相对反射率值即为图像每个像素的灰度值与整景图像灰度值均值的比值; ③ 计算地震后与地震前反射率差值图像, 并对反射率差值图像进行主成分变换; ④ 提取主成分变换后的第一主成分, 利用阈值分割初步提取震害区; ⑤ 利用第一步提取的居民区和灯光影像数据对震区结果进行两次优化, 得到最终的居民区震害信息. 该方法提取流程如图5所示.
图5 居民区震害信息提取流程图
2.1 数据
北京时间2001年1月26日11时16分印度古吉拉特邦发生MW7.6地震, 震中位于(23.2°N, 70°E). 该地震造成2万人死亡, 6万多人受伤, 20万人无家可归, 给古吉拉特邦地区的人民带来灾难性后果, 受此次地震影响最大的城市有帕焦、 普杰、 安贾尔和甘地特姆等(Narayanetal, 2002). 本文以安贾尔市附近的Landsat卫星TM图像为例, 提取该城市的震害信息. TM数据来源于美国USGS网站, 图像空间分辨率为30 m, 已经辐射校正和正射纠正. 印度MW7.6地震震前TM影像(图6a)于2001年1月8日获取, 图像整体质量较好, 左上部有少量薄云和阴影; 震后TM影像(图6b)于2001年2月9日获取, 图像质量好于震前. 另外, 从美国NOAA网站获取了该地区2001年常年稳定灯光的平均数据, 其空间分辨率为1 km.
图6 印度古吉拉特邦震区震前(a)和震后(b)TM影像
2.2 震害信息提取结果
由于TM数据和灯光数据均经过辐射和几何预处理, 且图像之间的几何配准精度高, 无需再进行校正和配准处理. 根据本文算法, 首先, 对TM影像采用内部平均法进行相对反射率反演; 其次, 以人工选择的居民区为参考样本, 利用光谱夹角匹配法(童庆禧等, 2006)提取似居民区专题图; 然后, 根据是否有灯光数值对图像进行二值化, 制作居民区掩膜, 用于后期对居民区的优化处理; 最后, 根据本文计算流程依次进行反射率差值、 主成分变换、 震害信息初步提取和震害信息优化等处理后, 得到最终的居民区震害信息.
图7a为基于震前与震后图像反射率差值的主成分图像初步提取的震害结果. 由于震后图像受薄云阴影(图中红色箭头所指区域)、 季节的差异及土壤水分含量等因素的影响, 薄云阴影、 山区大片植被和海边沙滩被错误地提取出来. 图7b为利用监督分类法提取的似居民区对图7a进行优化后的结果, 可以看出, 大部分非居民区的错误目标基本都被去除掉, 但是左上角有阴影的无植被覆盖的裸土地仍然存在. 图7c为利用灯光遥感数据对图7b进行优化的最终结果, 可以看出, 仅剩下几个居民区内有建筑物震害信息(图中红色圆圈内区域). 图8为图7c中红色圆圈内的放大图, 从左至右分别是安贾尔市震前、 震后以及震害信息叠加在震后影像上的效果. 对比图8震前与震后的图像可以看出, 本文方法提
图7 印度古吉拉特邦震区地震前后反射率图像第一主成分差值图像(a)、 利用似居民区数据
图8 安贾尔市震前(a)、 震后(b)TM图像及震害信息提取结果(c), 白色斑点区域即为震害区域
取的结果与安贾尔市反射率变大部分非常一致.
基于地震前后图像变化检测法提取震害信息不可避免地会受到多种因素的影响, 本文首先利用主成分变换对原始图像的震害信息进行了增强, 然后利用居民区数据和灯光数据等对非震害信息进行了剔除, 从而较好地解决了利用中分辨率遥感影像提取居民区震害信息时目标信息弱、 背景信息复杂等问题.
图9a, b给出了印度安贾尔市震前与震后相对反射率图像直接差值运算后的差值图像直方图, 图9c给出了对差值图像进行主成分变换后的主成分图像直方图. 可以看出, 主成分变换前所有波段的直方图具有相似的正态分布, 且直方图峰值均大于0, 这说明震后所有波段的反射率均大于震前, 该现象与图2, 3中相对反射率差值曲线一致. 但是经主成分变换后的直方图则发生了明显变化: 除第一主成分之外的其它主成分分量直方图均具有相似的正态分布, 但峰值均降低为0, 表明这几个主成分分量基本不包含震害信息; 第一主成分明显与其它分量分离开, 且第一主成分的峰值相对主成分变换前明显变大, 说明主成分变换能将平均分配在各个波段的震害信息集中到第一主成分量, 因此利用第一主成分分量更容易提取震害信息. 图10给出了基于相对反射率差值图像提取的震害信息和基于差值图像主成分变换后的第一主成分提取的震害信息. 可以看出, 经过主成分变换后的结果明显优于变换之前的结果.
图9 印度安贾尔市震前震后相对反射率差值图波段1(a)和波段2—4的直方图分布(b)
图10 基于相对反射率差值(a)和基于该差值主成分变换后(b)的印度安贾尔震区震害信息提取结果
许多非地震因素也会在地震前后差值图像上产生反射率差异, 使得利用差值算法提取的震害信息出现错误结果, 如图7a中云阴影、 植被和沙滩等都是容易产生误差的地物. 解决该问题的一种思路是采用间接法, 即将研究区内各种背景对象依次剔除; 另一种思路是采用直接法, 直接从图像中提取居民区进而提取震害区域. 理论分析和试验结果表明, 间接法需要针对不同的背景对象设计相应的算法, 由于背景对象复杂导致整体算法复杂, 且在每种背景对象消除过程中所产生的误差均会叠加到最终结果中. 相比而言, 直接法的算法和处理过程则相对简单, 只有一次误差累积过程; 但直接法会受到“同谱异物”的影响, 影像中与居民区具有相似特征的背景地物会被当作居民区而被提取(图7a). 为解决该问题, 本文引入灯光遥感数据, 采用直接法提取居民区. 众所周知, 居民区与裸土、 沙滩、 河床等地物最大的差别在于居民区有人居住, 而灯光是人类活动的典型标志之一, 因此利用灯光遥感影像可以将居民区与裸土等背景地物区分开. 具体作法分两步: ① 基于震前遥感数据, 以居民区为参照, 采用监督分类法提取似居民区, 由于植被、 水体、 云等与居民区在统计特征上差异较大, 所以似居民区内几乎不包括这类背景, 仅包含与居民区相似的裸土区或其它类似区域; ② 利用灯光图像对似居民区进行二次优化, 去除无人居住的裸土地等背景目标. 从图7的优化结果看, 利用似居民区能有效去除植被、 水体、 云等背景误差, 然后经过灯光遥感数据的再次优化, 居民区之外的误差几乎全部被消除.
本文利用中分辨率遥感影像进行地震灾情信息监测, 具有费用低、 数据获取机会多、 数据处理效率高及灾情监测范围大等优点, 对震后应急救援阶段快速提供宏观震害信息具有非常重要的意义. 通过对中分辨率遥感影像震害特征的分析表明, 在30—50 m中分辨率遥感影像上, 震后图像反射率变大是以建筑物为主的居民区最显著的震害特征, 因此本文建立了一种基于反射率变化检测的方法来提取居民区震害信息. 鉴于震后反射率相对震前反射率在遥感影像所有波段均具有相似的特征, 对震后-震前反射率差值采用主成分变换, 将均匀分布在各个波段的震害信息集中在第一主成分, 从而达到增强震害信息的目的. 另外, 采用监督分类法从震前影像中提取似居民区, 并利用灯光遥感数据提取居民区范围, 通过这两种方法有效剔除了居民区的非震害信息, 从而显著提高了直接基于反射率变化提取震害信息的精度.
董燕生, 陈洪萍, 潘耀忠, 方伟华. 2014. 基于震后高分辨率卫星遥感影像的建筑物瓦砾快速提取方法[J]. 应用基础与工程科学学报, 22(6): 1079--1088.
Dong Y S, Chen H P, Pan Y Z, Fang W H. 2014. A rapid method for extracting building rubble based on post-earthquake high-resolution satellite imagery[J].JournalofBasicScienceandEngineering, 22(6): 1079--1088 (in Chinese).
窦爱霞. 2003. 震害遥感图像变化检测技术研究[D]. 济南: 山东科技大学地球科学与工程学院: 1--21.
Dou A X. 2003.StudyonSeismicDamageChangeDetectingTechniqueBasedonRemoteSensingImages[D]. Ji’nan: College of Earth Science and Engineering, Shandong University of Science and Technology: 1--21 (in Chinese).
雷莉萍, 刘良云, 张丽, 毕建涛, 吴艳红, 焦全军, 张文娟. 2010. 汶川地震房屋倒塌的遥感监测与分析[J]. 遥感学报, 14(2): 333--344.
Lei L P, Liu L Y, Zhang L, Bi J T, Wu Y H, Jiao Q J, Zhang W J. 2010. Assessment and analysis of collapsing houses by aerial images in the Wenchuan earthquake[J].JournalofRemoteSensing, 14(2): 333--344 (in Chinese).
柳稼航. 2003. 利用遥感技术进行城市建筑震害的自动识别与分类方法研究[D]. 北京: 中国地震局地质研究所: 1--10.
Liu J H. 2003.AMethodStudyonAutomaticRecognitionandClassificationofEarthquake-CausedBuildingDamageinCitiesUsingRemoteSensing[D]. Beijing: Institute of Geology, China Earthquake Administration: 1--10 (in Chinese).
童庆禧, 张兵, 郑兰芬. 2006. 高光谱遥感: 原理、 技术与应用[M]. 北京: 高等教育出版社: 166--217.
Tong Q X, Zhang B, Zheng L F. 2006.HyperspectralRemoteSensing[M]. Beijing: Higher Education Press: 166--217 (in Chinese).
王晓青, 王龙, 王岩, 丁香, 窦爱霞, 张飞宇. 2008. 汶川8.0级大地震应急遥感震害评估研究[J]. 震灾防御技术, 3(3): 251--258.
Wang X Q, Wang L, Wang Y, Ding X, Dou A X, Zhang F Y. 2008. Study on the damage assessment of WenchuanMS8.0 earthquake based on remote sensing imagery[J].TechnologyforEarthquakeDisasterPrevention, 3(3): 251--258 (in Chinese).
曾涛, 杨武年, 黎小东, 刘汉湖, 彭立, 简季. 2010. 面向对象的高空间分辨率遥感影像信息提取: 汶川地震城市震害房屋案例研究[J]. 自然灾害学报, 19(5): 81--87.
Zeng T, Yang W N, Li X D, Liu H H, Peng L, Jian J. 2010. Information extraction of object-oriented high resolution remote sensing image: A case study on urban damaged buildings in Wenchuan earthquake[J].JournalofNaturalDisasters, 19(5): 81--87 (in Chinese).
张景发, 谢礼立, 陶夏新. 2002. 建筑物震害遥感图像的变化检测与震害评估[J]. 自然灾害报, 11(2): 59--64.
Zhang J F, Xie L L, Tao X X. 2002. Change detection of remote sensing image for earthquake-damaged buildings and its application in seismic disaster assessment[J].JournalofNaturalDisasters, 11(2): 59--64 (in Chinese).
张小咏, 买莹, 张凌. 2013. 基于数学形态学的高分遥感影像地震重点目标信息提取[J]. 地震, 33(2): 115--122.
Zhang X Y, Mai Y, Zhang L. 2013. Extraction of earthquake key objects from high spatial resolution images based on mathematical morphology[J].Earthquake, 33(2): 115--122 (in Chinese).
Estrada M, Yamazaki F. 2000. Use of Landsat images for the identification of damage due to the 1999 Kocaeli, Turkey earthquake[C]∥Proceedingsofthe21stAsianConferenceonRemoteSensing. Taipei: ACRS: 1--11.
Hayashi H, Hashitera S, Kohiyama M, Matsuoka M, Maki N, Fujita H, Elvidge C D. 2000. International collaboration for the early damaged area estimation system using DMSP/OLS nighttime images[C]∥ProceedingsofIEEE2000InternationalGeoscienceandRemoteSensingSymposium(CD-ROM). Honolulu: IEEE: 2697--2699.
Matsuoka M, Yamazaki F. 1998. Identification of damaged areas due to the 1995 Hyogoken-Nanbu earthquake using satellite optical images[C]∥Proceedingsofthe19thAsianConferenceonRemoteSensing. Traders Hotel, Manila: 1--6.
Matsuoka M, Yamazaki F. 2004. Use of SAR intensity imagery for earthquake damage detection[C]∥Proceedingsofthe2ndWorkshoponApplicationofRemoteSensingTechnologiesforDisasterResponse. Newport Beach, California: IEEE: 1--5.
Narayan J P, Sharma M L, Kumar A. 2002. A seismological report on the 26 January 2001 Bhuj, India earthquake[J].SeismolResLett, 73(3): 343--355.
Yusuf Y, Matsuoka M, Yamazaki F. 2001. Damage assessment after 2001 Gujarat earthquake using Landsat-7 satellite images[J].JIndianSocRemoteSens, 29(1/2): 17--22.
Yusuf Y, Matsuoka M, Yamazaki F. 2002. Detection of building damage due to the 2001 Gujurat, India earthquake, using satellite remote sensing[C]∥Proceedingsofthe7thNationalUSConferenceonEarthquakeEngineering. Oakland, California: Earthquake Engineering Research Institute: 21--25.
Extraction of seismic damage information of the residential area based on medium-resolution remote sensing image
1)NationalEarthquakeResponseSupportService,Beijing100049,China2)InstituteofRemoteSensingandDigitalEarth,ChineseAcademyofSciences,Beijing100094,China
There are some weak spots in medium-resolution remote sensing image, such as weak information in seismic disaster of buildings and change detection method prone to be affected by non-seismic disaster information. According to those problems, this paper established a method for quickly extracting residential seismic damage information based on the change detection method, which uses principal component transform to enhance damaged information, and adopts the supervised classification method to extract the similar residential area, which were further optimized by using the nighttime lights data, eliminating the influence of non-seismic disaster factors in the change detection method. And then the meizoseismal area of IndiaMW7.6 earthquake in 2001 is chosen as the studied area by using Landsat thematic mapper (TM) imagery pre- and post-earthquake and nighttime lights data of earthquake-stricken area so as to verify and analyze the method proposed in this paper. The results show that the reflectivity becomes larger is the the most significant damaged characteristic for the building residential area in the resolution of 30--50 m medium-resolution remote sensing image. The method for extracting residential earthquake damage information established in this paper has good effects in solving weakly damaged information and background noise based on medium-resolution images.
seismic damage information; residential area; change detection; medium-resolution remote sensing image; nighttime lights data
地震星火项目(XH14059)和高分辨率对地观测系统重大专项(31-Y30B09-9001-13/15)共同资助.
2015-11-09收到初稿, 2015-12-29决定采用修改稿.
e-mail: xyzhang2005@hotmail.com
10.11939/jass.2016.03.016
P315.9
A
张小咏, 李庆亭, 高娜. 2016. 基于中分辨率遥感影像的居民区震害信息提取. 地震学报, 38(3): 486--495. doi:10.11939/jass.2016.03.016.
Zhang X Y, Li Q T, Gao N. 2016. Extraction of seismic damage information of the residential area based on medium-resolution satellite image.ActaSeismologicaSinica, 38(3): 486--495. doi:10.11939/jass.2016.03.016.