基于D-InSAR技术的伊朗巴姆地震地表形变监测震灾防御技术

基于D-InSAR技术的伊朗巴姆地震地表形变监测1寇程柯长青（南京大学地理信息科学系，南京210046）摘要以伊朗巴姆地区为例，对伊朗巴姆地震造成的地表形变进行了差分干涉测量，得到了垂直向的同震三维形变场，并运用GIS三维分析技术对形变场进行了分析。实验结果表明，地震在巴姆城市的东侧造成了较大形变，在西侧也产生了微量形变。巴姆城市北部地块沉降，南部地块隆起。同时在巴姆城市南部可明显看到地震造成的断层。实验结果验证了基于C波段的SAR数据的D-InSAR技术在干燥地区监测地表形变方面的可行性。本文对产生去相关效应的原因进行了解释，认为对于干燥少植被的地区干涉效果较好。并指出，如果能够通过技术进步提高雷达干涉测量的精度并降低观测成本，同时将该技术与GPS、GIS等技术相结合，从而更好地研究形变机理，这将对地质灾害的研究产生重大意义。关键词：差分合成孔径雷达干涉测量（D-InSAR）同震形变场去相关效应ENVISATASAR数据伊朗巴姆地震引言传统的地表形变监测方法一般有水准测量、三角测量、三边测距（查显杰，2007）、GPS测量（刘大杰等，1999）等，但是这些测量方法都具有范围小、难度大、周期长、成本高的缺陷。合成孔径雷达干涉测量（InterferometricSyntheticApertureRadar,InSAR）技术是以合成孔径雷达（SyntheticApertureRadar,SAR）复数据提取的相位信息为数据源获取地表三维信息和变化信息的一项技术（王超等，2002），差分合成孔径雷达干涉测量（D-InSAR）技术是InSAR技术的延伸，差分是指在研究区形变前后的SAR影像产生的干涉图中去除地球曲率、地形因素等产生的相位，从而得到地表形变信息。它是迄今为止独一无二的基于面观测的形变遥感监测手段（刘国祥，2006）。与常规方法相比，D-InSAR技术监测地表形变具有大面积、连续、快速、准确的优势，同时观测无需布设地面控制点，从而降低了对危险地区1基金项目国家自然科学基金（40971044）和教育部新世纪人才项目（NCET-08-0276）联合资助[收稿日期]2012-09-16[作者简介]寇程，男，生于1987年。南京大学地理信息科学系硕士研究生。研究方向：遥感地学应用。E-mail:[e-mail]通讯作者柯长青，男，生于1969年。南京大学地理信息科学系教授。主要研究方向：遥感，GIS，GPS及其应用。E-mail:[e-mail]第8卷第1期2013年3月震灾防御技术TechnologyforEarthquakeDisasterPreventionVol.8,No.1Mar.,2013寇程，柯长青，2013．基于D-InSAR技术的伊朗巴姆地震地表形变检测．震灾防御技术，8（1）：72—80.1期寇程等：基于D-InSAR技术的伊朗巴姆地震地表形变监测73（比如火山地区）观测的危险性（马超等，2004）。国内外学者对D-InSAR技术监测地表形变监测进行了广泛研究。D-InSAR技术被用于台湾集集地区地震（刘国祥等，2002）、昆仑山口西地震（单新建等，2004）、张北地震（张红等，2002）、汶川地震（孙建宝等，2008）、Landers地震（Zebker等，1994）和Hyogoken-Nanbu地震（Ito等，2003）同震地表位移监测和同震形变场参数提取。这些研究证明了D-InSAR技术在地表形变监测方面的可行性，其监测精度可达到厘米级。2003年12月26日，伊朗巴姆地震发生后，许多研究者对巴姆地震进行了研究（季灵运等，2009；凌勇等，2006；罗扬等，2006；王志勇等，2008），得到了巴姆地震LOS方向上的同震形变场。本文采用二轨法利用ENVISARscape模块，基于伊朗巴姆地震前后的两景SAR影像和SRTMDEM数据进行差分干涉处理，并利用雷达波束的局部入射角得到垂直方向上的同震形变场。最后结合GIS三维空间分析技术分析了地震造成的形变特征，并探讨了去相关效应产生的原因。1研究区地理背景及数据源巴姆市位于伊朗克尔曼省首府克尔曼市东南195km处，地处卢特沙漠的边缘，海拔1000m左右（张保淑，2009），地理位置见图1。这里干燥少雨，多风（Bam(Iran),2003）。2003年12月26日，世界时01:56:52，巴姆地区发生了里氏6.7级地震，震中位于29.01°N，58.26°E（李毅等，2008），距巴姆城市约15km。巴姆古城在这次地震中几乎毁于一旦（舒宁，2003）。本文采用ENVISAT卫星C波段的两景ASARImageSLC数据，数据的部分参数见表1。利用重采样后的SRTMDEM数据进行地形相位的去除，格网大小为25m×25m，UTM投影40分带，参考椭球体采用WGS84。研究区SAR影像见图1。从图1可以看出，数据覆盖范围大部分地方地形平坦，西南角有部分覆盖了山区。图1巴姆地区地理位置及研究区范围Fig.1ThelocationofBamandthestudyarea震灾防御技术8卷74表1ENVISATASAR数据部分参数Table1SomeparametersofENVISATASARdata影像轨道号成像时间极化方式轨道类型时间间隔(天)主影像91922003-12-03VV降轨70副影像101942004-02-112D-InSAR技术原理Zebker等（1994）、舒宁（2003）、王超等（2002）对D-InSAR技术测量形变的基本原理已经有了非常详细论述，本文不再赘述。能够进行干涉处理的SAR影像必须处理成SLC（SingleLookComplex）格式，SLC影像上的每个象元可表示为：(,)(,)(,)eimnSmnAmnϕ=⋅（1）式中，A(m,n)是影像的强度信息即振幅；φ(m,n)是相位信息。干涉处理就是将两景配准后的影像逐象元进行复共轭相乘，记两景影像象元分别是（何秀凤等，2012)：12(,)11(,)22(,)(,)e(,)(,)eimnimnSmnAmnSmnAmnϕϕ=⋅=⋅（2）逐象元进行复共轭相乘后得到：12*12[(,)(,)]12(,)(,)(,)(,)(,)eimnmnSmnSmnSmnAmnAmnϕϕ−=⋅=⋅⋅（3）因此，干涉相位为：12(,)(,)(,)mnmnmnϕϕϕ=−（4）这也是最后干涉图上每个象元的干涉相位信息，φ(m,n)的取值范围为(-π,π)。但是，地球表面是一个具有一定曲率的球面，微波信号在大气中传播时还会受到大气的影响，而且雷达影像有一定的噪音，所以一般认为干涉相位由以下几部分组成（刘国祥等，2002）：flat_Earthtopodisppathnoise2nϕϕϕϕϕϕπ=+++++⋅（5）式中，φflat_Earth是地球曲率产生的相位信息，也就是高度不变的平地在干涉图中表示出来的干涉条纹随距离向和方位向的变化而呈周期性变化的现象（王超等，2002）；φtopo是地形产生的相位信息，是由于地形起伏造成的，可用于建立DEM；φdisp是由地表形变产生的相位信息，反映了地表形变信息；φpath是雷达信号在传播过程中受到大气影响产生的相位信息；φnoise是雷达影像噪声产生的相位信息；n·2π是相位的整周数，需要通过相位解缠获得。3实验验证3.1实验步骤采用D-InSAR技术进行地表形变监测主要包括以下几个步骤：基线估计、影像配准、主干涉图生成、平地相位及地形相位的消除和相位解缠形成形变图。数据处理流程见图2（邵芸等，2002）。1期寇程等：基于D-InSAR技术的伊朗巴姆地震地表形变监测75副影像主影像基线估计影像配准主干涉图的生成平地相位消除模拟干涉图生成DEM轨道参数差分干涉图生成相位解缠地理编码区域形变场图2二轨法差分干涉测量流程图（邵芸等，2002）Fig.2Thedataprocessingprocedurefortwo-passdifferentialradarinterferometry(fromShaoetal.,2002)①基线估计。干涉像对的基线大小直接影响着干涉图的质量，基线过大则无法产生干涉。因此在进行干涉处理前，首先要进行基线估计，以确定像对能否产生干涉。一般认为ENVISAT的C波段数据临界基线大约为1100m，最佳基线应小于200m，而用于地面变化监测的基线大小应小于70m（舒宁，2003），本文所用数据的垂直基线估计结果为9.1m，可以产生干涉。②影像配准。进行干涉的两幅SAR影像之间必须进行精确地配准，一般认为配准精度应达到亚像元级（Ito等，2003）。本文利用SAR数据自带的卫星轨道参数，进行自动精确配准（Klees等，1999）。配准精度可以通过相干系数图（图3）间接地进行检查。③主干涉图生成。配准后的主、副影像相应象元进行复共轭相乘，得到主干涉图（图4）。④平地相位及地形相位的去除。在主干涉图上去除平地相位和地形相位以得到差分干涉图（图5），其中，平地相位可通过对干涉条纹乘以复相位函数来去除（王超等，2002）。对于地形相位的消除，本文则是利用该地区的SRTMDEM数据和进行干涉的两景SAR影像的天线姿态参数生成地形的模拟干涉图，从而得到地形相位。⑤相位解缠形成形变图。最后对差分干涉图采用最小费用流（MinimumCostFlow,MCF）方法进行相位解缠得到形变图（图6）。3.2实验结果图4是两景SAR影像形成的主干涉图，干涉条纹从黑色到白色的一个颜色周期代表了相位差从-π到π的变化，可以看出在巴姆城市的东侧形成2个明显的花瓣状干涉条纹，周围的干涉条纹主要是由于地球曲率和地形造成的。图5是去除了平地相位和地形相位的差分干涉图，同样，从黑色到白色的一个颜色周期代表了相位差从-π到π的变化，可以看到在巴姆城市周围形成了4个花瓣状（或蝴蝶状）的干涉条纹，东侧的条纹多且密集，而西侧的条纹少且稀疏。北侧两组条纹的颜色变化从里到外是由白到黑，而南侧条纹的颜色变化从里到外是由黑到白，说明巴姆城市南北两侧的形变方向不同，这一点从最后得到的同震形变场也可以看到。在干涉图上可以看到，一些地区（巴姆市及其东南侧）干涉条纹十分杂乱，这种现象称为去相关现象，发生去相关现象的斑块与城市区域和植被区域吻合得很好，所以可以初步推断去相关现象的发生是由植被的变化和城市建筑的毁坏造成的。另外在巴姆市的西侧有条带状的去相关区域，通过光学遥感影像可以看出这里是河流，所以可以初步认定这里的去相关现象是由河流造成的。发生破坏的城市、植被和河流的散射特性极不稳定，在两景影像获取的70天之内发生了极大地变化，所以这些地方的干涉效果比较差，去相关现象十分严重。去震灾防御技术8卷76相关现象也可以从两景影像的相干系数图中看出，两景影像成像之间散射特性的改变会降低两景影像之间的相关性，去相关现象越严重的地方，相关系数越低，如图3中的黑色地方（相干系数接近于0）。相干系数图也可以用来评价干涉质量（舒宁，2003）。舒宁（2003）提出变形监测采用雷达干涉测量方法时必须满足2个条件：一是2次观测期间的变化梯度必须小于一定阈值；二是2次观测期间，一个象元内的所有地物散射体的特性应具有相似性，特别是表面位置变化的误差应只是雷达波长的10%—20%的量级。去相关效应与干涉的两景SAR影像之间发生变化的程度等有关，大气影响、人类活动等都有可能引起去相关效应（Klees等，1999)。因此，雷达干涉测量监测地表形变的方法适用于地表变化在一定范围内、干燥、少植被的地区。图3相干系数图Fig.3Coherencecoefficientimage图4主干涉图Fig.4Themajorinterferogram图5差分干涉图Fig.5Thedifferettntialinterferogram图6同震形变图（单位：cm）Fig.6Coseismicdeformationmap(unit:cm)3.3地表同震形变场分析从干涉图直接解缠得到的形变是沿着视线方向（LineofSight,