遥感技术因具有覆盖范围广、 图像获取方便等特点,能够客观、 全面地反映地震后灾区的景观,能为震害调查、 损失快速评估提供科学依据,对地震经济损失评估也具有重要意义. 我国曾对1966年邢台地震、 1975年海城地震、 1976年唐山地震等大震进行了震后灾区航空摄影,并积累了丰富的震害遥感影像判读经验(陈鑫连,1995; 魏成阶等,1993). 在2003年的伽师6.8级地震中,王晓青等根据以往震害影像统计并结合本次地震震害遥感特征,提出了遥感震害分级分类标准和地震烈度划分标准,进而得到基于震害遥感影像的伽师地震等震线图(王晓青等,2003). 然而光学影像由于受天气状况的制约,如果震后出现云雨等恶劣天气,则利用光学影像进行震害调查就会受到极大限制. 而雷达因为具有全天候、 全天时的特点,就可以克服上述弱点. 张景发等(2002)曾利用合成孔径雷达(synthetic aperture radar,简写为SAR)图像对张北地震进行了震害评估研究; 另外,Yonezawa和Takeuchi(2001)对1995年日本兵库县南部(Hyogoken-nanbu)地震前后的ERS-1/SAR图像的相干性进行研究,发现地震前后的两幅SAR图像有明显的失相干性,这种失相干性的大小、 分布与地震中地面破坏程度、 房屋倒塌情况及其分布状况有明显的关系. 需要注意的是,这种方法受到数据参数(如基线长、 间隔时长、 波长、 地表变化及季节植被变化等因素的影响,需要对图像的失相干性原因进行认真分析,判断由哪种因素引起的,具有一定的局限性(Zebker,Villasenor,1992). 而SAR幅度图像则较少受到这些因素的制约,因此可以利用后向散射系数来评估地震、 洪水等自然灾害导致的受灾地区. Giovanna和Paolo(2008)利用SAR幅度图像研究2003年Algeria地震及2007年Peru 地震,结合GIS手段经过分类得到了震害分布图,准确率较高. Matsuoka和Yamazaki(2004)利用ERS的SAR数据,根据后向散射系数差异与图像相干系数之间的统计关系,给出了一个区分指数来对建筑物在地震中的受损程度进行分类,在神户地震的震例研究中取得了很好的效果.
汶川地震后,科学家们运用各种手段对此次地震的成因进行了研究,各个研究组的研究成果相继在学术期刊和Internet上公布(张培震等,2008; 王卫民等,2008; 滕吉文等,2008; 何宏林等,2008). 但利用遥感技术进行震害评估的并不多,这是因为此次地震破坏面积广,且大部分在山区,震后出现了多天的阴雨天气,给高分辨率光学影像的解译带来了很大困难. 基于SAR图像具有全天时、 全天候的特性,本文拟利用地震前后SAR幅度图像变化检测、 干涉相干图像的相干系数相结合,对地震灾害地区进行快速识别与评估.
1 研究区域及实验数据选取汶川县位于阿坝州境东南部的岷江两岸,县城威州镇位于县北部杂谷脑河与岷江交汇地,映秀镇地处阿坝州南大门,距成都88 km,是进出九寨沟、 卧龙、 四姑娘山的必经之地,是重要的交通要道. 2008年5月12日下午2时28分,汶川发生MS8.0强烈地震,汶川全县遭受严重损失,主要表现为: ① 受灾面积大,全县13个乡镇全面受灾. 截至6月27日全县死亡15 941人,受伤34 584人; ② 受灾程度深. 全县耕地面积7 100公顷(106 500亩),受灾6 000公顷(90 000亩),房屋倒塌21万间,损毁房屋30.35万间; ③ 基础设施毁损严重. 道路垮塌,水、 电、 气、 通讯全部中断,工矿商贸企业80%被损毁; ④ 地震次生灾害威胁特别严重. 汶川县原有地质灾害点160处,地震后新增3 590处,县城周边新增的79处地质灾害点随时威胁县城的安全,已形成“晴天沙尘暴,雨天泥石流”的现实状况,县城已被地震次生灾害威胁所包围. 随后几天多次发生强烈余震,映秀镇是震中和重灾区,全镇大部分房屋倒塌,到处山体滑坡,造成停水、 停电,通讯、 交通中断. 本文以映秀镇及周边地区为研究范围,图 1为研究区ETM合成图.
![]() |
图 1 (a)研究区震前(2007年9月18日)ETM 7,4,2波段合成图;(b)交通位置图 Fig. 1 (a)L and sat 7 pre-event image(ETM b and 7,4,2)of the test area;(b)location of Yingxiu |
地震发生后第二天早晨起,灾区大部分地区出现了10 mm以上的中雨,在震后的半个多月中,天气也一直以阴雨为主,这就给光学成像带来了很大的困难,合成孔径雷达(SAR)成为这次抗震救灾前期遥感信息保障的十分重要的数据源. 震后欧空局提供了大量的SAR数据,若要将覆盖灾区的影像全部都进行处理,数据量是相当庞大的,故只挑选覆盖映秀和都江堰地区的ASAR数据进行计算,使用到的数据如表 1所示.
![]() |
表 1 震前/震后的ASAR数据 Table 1 ASAR data and baseline distances in the test area |
实验用到的数据中B1和B2幅为震前获取的,A1幅为震后获取的,数据格式均为单视复数据(SLC),结合ENVISAT的精密轨道数据将其处理为多视幅度图像,并采用经典的Lee算法进行滤波处理. 一般而言,较大的估算窗口具有较高的估算性能,但会降低空间分辨率,考虑两者间的平衡,本文估算窗口的大小为21×21. 最后将所有图像都与2007 年8月6日的震前图像进行配准. 震前图像B1如图 2a所示,震后图像A1如图 2b所示. 图中白色方框为研究区范围; 图 2c,d分别为从图 2a,b中裁剪的研究区图像.
![]() |
图 2 幅度图像.(a)震前;(b)震后;(c)映秀地区震前;(d)映秀地区震后 Fig. 2 Amplitude images before(a) and after(b)the earthquake.(c) and (d)are amplified images for the box region in(a) and (b),respectively |
直接对比震前震后的SAR图像,我们很难得到一些变化信息,这是因为当我们面对SAR图像时,我们看到的是完全不同于人们肉眼所见的特征. 人的肉眼是一种被动式“传感器”,接收的是地物反射太阳光的能量. 而雷达则是一种工作在微波波段的主动式传感器,它发射某一特定波长的微波波段的电磁波,接收来自地面的后向散射电磁波能量,两者相干而成像. 上面所看到的幅度图像是一种灰度图像,图像上色调的变化取决于目标物的后向散射截面. 每一个接收到的回波被转换成电信号以某一具有特定值的、 用于表示亮度的数字化象元. 尽管坡度的变化、 含水量和表面粗糙程度是影响雷达图像色调的3个主要作用,但是表面粗糙度在决定雷达图像的灰度方面起着决定性的作用. 而地震带来的巨大破坏则会造成震前震后同一物体的粗糙程度发生很大的变化,从而使其回波强度发生变化,Matsuoka和Yamazaki(2000a)曾用一幅房屋倒塌后的废墟造成漫反射这一生动形象的图画来反映该现象. Aoki等在研究1995年神户地震中发现,人工建筑由于在结构体与地面之间形成一种“角反射”效应而具有较高的反射强度,并且反射具有较好的方向性; 而开阔的场地或者受损的建筑物由于其表面发生漫反射而回波强度相对较低(Aoki et al,1998; Matsuoka,Yamazaki,2000b). 因此,与震前的图像相比较,由于建筑物的倒塌或者废墟被清理而露出地表,都会造成震后的后向散射系数变小. 从图 2幅度图像上可以看出地震前后一些地区的灰度强度发生了变化,这其中有些是由于两幅图像接收季节不同,时间间隔较长一些地物发生的自然变化,然而相当一部分应该是由地震带来的破坏而造成的. 尽管如此,仅凭肉眼我们还是很难看出这些发生变化的地区的分布情况,也就是地震影响较大地区的分布. 下面通过进一步的数据处理来进行分析.
2 利用SAR幅度图像做比值变化检测要想对两幅图像进行变化检测,必须要对不同期次的图像进行精确配准. 航空图像由于不同成像期次成像轨道的差异给像素级匹配造成困难,而对于SAR图像来说,因为采用同一卫星系统,这并不造成太大的困难. 利用SAR图像做变化检测有两种情况(Rignot,van Zyl,1993). 其中一种情况是检测目标形态的变化,这种变化检测要求不同时相的图像要精确配准,但并不需要绝对定标数据. 对于精确配准的两个时相图像来检测变化,最常规的方法是图像相减或比值处理. 其它的处理方法,如多时相数据分类及主成分变化等,已被光学遥感数据处理的经验证明效果不如图像相减或比值方法. 而比值的分布只依赖于相对变化,因此从统计模型来看,比值方法比图像相减方法更适用于变化检测(郭华东,2000). 选用比值方法还有一个理由就是比值方法对辐射误差有更强的适应性,因为很多辐射误差是乘性的. 图 3a为2007年8月6日与2008年3月3日幅度图像做比值检测的结果,两幅均为震前图像. 从图 3a可以看出,图右下角都江堰市及周边乡镇颜色比其它山区的颜色深一些,出现一条明显的山区与平原地区的分界线. 这是由于平原地区人工建筑多于山区,从而地物特征随季节变化而出现这样的差异; 平原地区修建的永久性建筑较多,随时间推移变化不大; 而其它山区植被覆盖广泛,由于做比值检测的两幅图像一幅在夏季8月,一幅在春季3月,植被差异较大,造成颜色较浅.
![]() |
图 3 震前及同震比值图像.(a)B1与B2的比值图像;(b)B2与A1的比值图像 Fig. 3 Ratio of pre-earthquake and coseismic amplitude images.(a)Ratio of B1 over B2;(b)B2 over A1 |
为观察细节变化,从整幅图像中截取出研究区域. 图 4a,b是分别从图 3a,b中截取的映秀镇地区的图像,图中白色椭圆处为映秀镇所在地. 由于地形影响,处于山体阴影中的象元同样具有极低的反射系数值,与水体一样呈现暗色调,因此在图 2的SAR图像中并不能辨别出紫坪铺水库的轮廓. 经过对幅度图像进行比值处理,在图 4中紫坪铺水库已经清晰可见,并且在震前两景幅度图像间的比值图像与同震的比值图像中表现出一定的差异,这种差异可能是由于水质的变化造成后向反射系数变化造成的. 在图 4a的两幅震前SAR比值图像中,映秀镇的整体色调呈灰色,表明在接收这两幅图像的期间没有发生什么变化; 在图 4b中,图像中映秀镇及周边的色调已经不再是单一的灰色调,出现了从黑到白的变化. 黑色区域分布在沿江的建筑区,表明建筑物在地震中受到了巨大破坏,与3个月前的图像相比发生了巨大变化; 白色区域分布在山坡上,显然是由于滑坡所致. 图中方框圈住的区域距离映秀镇大约3 km,在图 4a中岷江河道清晰可见,而在图 4b中已模糊不清,这段正是在地震中严重破坏的都江堰至汶川公路的一段. 当时进入映秀救援的冲锋舟也只能到达这里,然后徒步翻越这几公里的山体滑坡路段进入映秀. 其实这段已不能称之为“路段”,因为整片的山坡因地震而坍塌下来,几米高、 甚至十几米高的巨大岩石从山顶遍布到山脚下的岷江中,将原来的公路全部掩埋.
![]() |
图 4 映秀、 漩口地区SAR比值图像 (a)震前B1与B2的比值图像;(b)同震B2与A1的比值图像 Fig. 4 SAR ratio image in Yingxiu and Xuankou (a)Pre-earthquake,B1 over B2;(b)coseismic,B2 over A1 |
干涉产品中的相干系数通常用于评价干涉相位的质量. 以前对相干系数的应用研究比较少,相干系数对SAR图像散射体特性的变化十分敏感,因此它能用于地表变化的探测、 土地利用的分类和绘制地震破坏程度评估等. 重复观测获得的雷达图像的失相干可能由以下3个因素引起: ① 系统热噪声引起的失相干; ② 空间基线失相干; ③ 时相失相干. 总的相干系数可表示为
式中
SNR是信噪比,B是两次观测卫星天线之间的距离,Ry是方位向的分辨率,θ是视角,λ是波长,R是雷达斜距. ρa是通过复数图像来计算的,即将信号的幅度和相位用复数表示; ρn是热噪声失相干系数; ρB是空间基线失相干系数. 由此可以看出,已知SNR和B即可由ρa推算出变化检测所需要的ρs. 因为对一景SAR图像来说,在幅宽范围内,R和cosθ的变化不到10%和5%,而Ry与λ一般是常数,B一般很小,所以在有大的SNR情况下,可以简单地通过ρa来检测变化. 因此相干系数可表示为
式中,S1和S2为经过预滤波的单视SAR图像,w为进行相干性估计的矩形窗口.
时间失相干是在重复轨道干涉中,雷达在不同时间飞临同一目标区域,但目标已发生某种变化所引起的失相干效应. 而地震发生造成建筑物破坏则正好符合这一条件. 因此在震害识别中关心的正是这种失相干. 下面利用SAR相干系数图像的失相干性来探讨震害程度. 数据处理步骤如图 5所示.
![]() |
图 5 震害估计图像处理流程图 Fig. 5 Flow chart of change detection using remote sensing before and after the earthquake |
基本工作包括以下诸方面: ① 准备3幅研究区域的ASAR幅度图像(其中震前两幅,震后一幅),以其中一幅震前数据为基准,另外两幅图像与它配准; ② 对配准后的图像进行滤波处理,采用经典的散斑噪声滤除算法 Lee 滤波算法,窗口大小为21×21; ③ 分别求取震前两幅图像的相干系数图rbb及震前震后两幅图像的相干系数图rab,采用的窗口大小为5×25(距离向和方位向),从而求得两幅相干系数图的比值图像rnd=rab/rbb.
图 6为计算结果. 其中图 6a为用B1和B2两幅图像干涉处理得到的相干图像,将其称之为rbb; 图 6b为用B1和A1干涉处理得到的相干图像,将其称之为rab; 图 6c为用rbb和rab得到的比值图rnd. 因为采用的ENVISAT为C波段,波长只有5.6 cm,在山区由于雷达图像叠掩、 前缩、 阴影,以及致密的植被覆盖等因素的影响,雷达回波中噪声加大,很难得到干涉信息,因此在图 6a中大部分区域为黑色,表现为失相干; 在都江堰平原地区为一片亮色区域,表明变化不大; 在映秀镇所在地也可以观察到一块白色区域(图中圆圈处),说明这些地方都保持了较高的相干性. 而在图 6b中这些相干性较高的区域颜色也已经变得很暗淡. 而用来干涉的两幅图像B1和A1的基线比图 6a中用来干涉的两幅图像B1和B2的基线还要短,这说明由于地震的原因,使两幅图像发生了时间失相干. 但是从图 6b中还难以圈定地震破坏区域,为此将图 6a,b两幅图像做比值处理,得到图 6c,从而使映秀镇等在地震中受损严重的区域在图中以白色显现了出来.
![]() |
图 6 映秀及周边地区失相干图像 (a)震前B1和B2相干图像rbb;(b)同震B2和A1相干图像rab;(c)rab与rbb比值图像rnd Fig. 6 Decorrelation image in and around Yingxiu (a)Pre-earthquake coherence image rbb;(b)coseismic coherence image rab; (c)ratio coherence image rnd of rab and rbb |
从图 6中可以看出,由灾害引起的失相干区域的分布范围与图 4中比值图像的结果非常相似. 以上述及的变化可通过计算地震前后图像之间的相干性等参数来定量地表示. 通常,房屋倒塌率越高,相干性越小,平均灰度差异性越大,平均灰度方差差异性越高. 因此,可以通过比较相干性系数、 平均灰度差异和平均灰度方差差异等来评估震害.
文中挑选了几个典型地区做相关统计(表 2). 统计结果表明,由于地震的破坏,造成了图像的相干系数大大下降. 当然这并不排除别的因素的影响,最主要的有基线长度的变化、 时间失相干等. 用来计算的震前与同震图像之间的基线长度均在300 m左右,时间间隔为3个月,因此由于建筑物的破坏造成的失相干占主要因素. 当然由于各种因素的影响,可能还存在异常,比如汶川县城震前、 震后相干系数相差不大,可能是由于建筑物虽然受到破坏,但倒塌率并不高的缘故所致.
![]() |
表 2 汶川地震前后部分乡镇图像相干系数统计表 Table 2 Statistics of correlation of SAR interferograms for the village houses before and after the Wenchuan earthquake |
我们也试图想得到建筑物损毁程度与失相干性之间的定量关系,然而这需要大量的统计样本及实地震害调查的数据,经过统计分析确定一定的阈值范围,将一定阈值范围内的相干系数归入相应的破坏程度范围内. 但判断破坏严重程度的阈值的确定是很困难的,不同地区可能有较大的差异.
4 讨论与结论SAR是一种工作在微波波段的相干成像雷达,其主动发射电磁波能量的特点,对地表物理特性和几何形状比较敏感的特性,多波段、 多极化散射特征、 极化测量及干涉成像等方式,正成为对地观测中最重要的前沿领域之一. 由于具有全天候、 全天时数据获取能力,SAR已成为遥感变化检测数据获取的重要技术手段,尤其对传统的光学传感器成像困难地区有着特别重要的意义,如多云多雨地区、 南北极地区等. 因此,其在震后灾害评估领域的应用前景十分广阔,国内外也已经开展了一系列的研究工作,尤其在利用相干图像识别地震破裂带方面取得了较好的效果(Fielding et al,2005; Liu et al,2001).
本文利用SAR幅度图像,通过比值法变化检测及干涉相关图像失相干分析方法对汶川地震中受损区域进行了识别,结果显示两种方法均能识别出映秀镇在地震前后发生的巨大变化. 直接利用震前、 震后的SAR幅度图像做比值处理来做变化检测在紫坪铺水库周边效果比较突出,能识别水体范围及震前震后水位的变化; 而利用SAR干涉相干图像的失相干分析在都江堰城区等平原地区效果较明显,并且建筑物的破坏程度与相干系数变化指数的大小高度相关,这为SAR图像在震害识别方面开辟了一条新途径. 另外,本文用到的ENVISAT数据的波长为5.6 cm,该波段数据在山区及植被茂密地区严重失相干,如果采用波长较长的ALOS数据,效果会好些. 考虑到这一因素,本文用震前的干涉相干图像与同震的干涉相干图像做了比值处理,这样可以去除一些空间及系统热噪声因素引起的图像失相干,而只保留时间失相干对图像的影响,更有利于识别震害地区的分布.
在震后的SAR幅度图像上,可以发现一些有用的震害信息,然而用人工辨别这些震害信息需要花费大量的时间,还可能造成遗漏. 本文应用SAR图像做变化检测,能够快速圈定一些受损严重区域,从实际的分析结果发现能定性或定量地检测出震害的一些信息. 然而由于SAR数据源及质量的限制,本文没有定量地给出地震前后图像的相干性大小与震害程度之间的关系,这也是今后要进一步开展的工作. 相信在不久的将来,随着更多高分辨率SAR卫星的上天,以及极化干涉SAR信息提取技术的发展,利用相干系数及后向散射强度的变化信息来提取震害信息将会发挥更大的作用.
本文得到欧空局提供的ASAR数据,荷兰Delft大学提供的Doris处理软件,CREASO 提供的SARscape软件及技术支持,NASA提供的DEM数据. 在此一并表示感谢.
[1] |
Matsuoka M, Yamazaki F. 2000b. Interferometric characterization of areas damaged by the 1995 Kobe earthquake using satellite SAR images[C/CD]∥Proceedings of the12th World Conference on Earthquake Engineering, 2, CD-ROM, ID2141.(![]() |
[2] |
Rignot E J M, van Zyl J J. 1993. Change detection techniques for ERS-1 SAR data[J]. IEEE Transactions on, Geoscience and Remote Sensing, 31(4): 896-906.(![]() |
[3] |
Yonezawa C, Takeuchi S. 2001. Decorrelation of SAR data uban damages caused by the 1995 Hyogoken-nanbu earthquake[J]. International Journal of Remote Sensing, 22(8): 1585-1600.(![]() |
[4] |
Zebker H A, Villasenor J. 1992. Decorrelation in interferometric radar echoes[J]. IEEETransactions on Geoscience and Remote Sensing, 30(5): 950-959.(![]() |
[5] |
陈鑫连. 1995. 地震灾害的航空遥感信息快速评估与救灾决策[M]. 北京:地震出版社: 1-123.(![]() |
[6] |
郭华东. 2000. 雷达对地观测理论与应用[M] . 北京:科学出版社: 1-517.(![]() |
[7] |
何宏林, 孙昭民, 王世元, 董绍鹏. 2008. 汶川MS8. 0地震地表破裂带[J]. 地震地质, 30(2): 359-362.(![]() |
[8] |
滕吉文, 白登海, 杨辉, 闫雅芬, 张洪双, 张永谦, 阮小敏. 2008. 2008汶川MS8. 0地震发生的深层过程和动力学响应 [J]. 地球物理学报, 51(5): 1385-1402.(![]() |
[9] |
王卫民, 赵连锋, 李娟, 姚振兴. 2008. 四川汶川8. 0级地震震源过程[J]. 地球物理学报, 51(5): 1403-1410.(![]() |
[10] |
王晓青, 魏成阶, 苗崇刚, 张景发, 单新建, 马庆尊. 2003. 震害遥感快速提取研究:以2003年2月24日巴楚-伽师 6. 8级地震为例[J]. 地学前沿, 10(s1): 285-291.(![]() |
[11] |
魏成阶, 朱博勤, 张宗科. 1993. 地震灾害航空遥感快速调查技术研究:以唐山地震区作模拟试验场[C]∥何建邦, 田国良, 王劲峰编. 重大自然灾害遥感监测与评估研究进展. 北京:科学出版社: 12-23.(![]() |
[12] |
张景发, 谢礼立, 陶夏新. 2002. 建筑物震害遥感图像的变化检测与震害评估[J]. 自然灾害学报, 11(2): 59-64.(![]() |
[13] |
张培震, 徐锡伟, 闻学泽, 冉勇康. 2008. 2008年汶川8. 0级地震发震断裂的滑动速率、复发周期和构造成因[J]. 地球物理学报, 51(4): 1066-1073.(![]() |
[14] |
Aoki H, Matsuoka M, Yamazaki F. 1998. Characteristics of satellite SAR images in the damaged areas due to the Hyogoken-Nanbu earthquake[C]∥Proceedings of the19th Asian Conference of Remote Sensing, C7. Manila: Asian Association on Remote Sensing: 1-6.(![]() |
[15] |
Fielding E J, Talebian M, Rosen P A, Nazari H, Jackson J A, Ghorashi M, Walker R. 2005. Surface ruptures and building damage of the 2003 Bam, Iran, earthquake mapped by satellite synthetic aperture radar interferometric correlation[J]. J Geophy Res, 110(B03302): 1-15.(![]() |
[16] |
Giovanna T, Paolo G. 2008. Damage Detection from SAR Imagery: Application to the 2003 Algeria and 2007 Peru earthquakes[J]. International Journal of Navigation and Observation. Article ID 762378, 8 pages.(![]() |
[17] |
Liu, J G, Black A, Lee H, Hanaizumi. 2001. Land surface change detection in a desert area in Algeria using multi-temporal ERS SAR coherence images[J].Int J Remote Sensing, 22(13): 2463-2477.(![]() |
[18] |
Matsuoka M, Yamazaki F. 2004. Use of satellite SAR intensity imagery for detecting building areas damaged due to earthquakes[J]. Earthquake Spectra, 20(3): 975-994.(![]() |
[19] |
Matsuoka M, Yamazaki F. 2000a. Use of interferometric satellite SAR for earthquake damage detection[C]∥Proceedings of the6th International Conference on Seismic Zonation. Earthquake Engineering Research Institute: 103-108.(![]() |