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

寇 程 柯长青

（南京大学地理信息科学系，南京 210046）

引言

传统的地表形变监测方法一般有水准测量、三角测量、三边测距（查显杰，2007）、GPS测量（刘大杰等，1999）等，但是这些测量方法都具有范围小、难度大、周期长、成本高的缺陷。合成孔径雷达干涉测量（Interferometric Synthetic Aperture Radar, InSAR）技术是以合成孔径雷达（Synthetic Aperture Radar, SAR）复数据提取的相位信息为数据源获取地表三维信息和变化信息的一项技术（王超等，2002），差分合成孔径雷达干涉测量（D-InSAR）技术是InSAR技术的延伸，差分是指在研究区形变前后的SAR影像产生的干涉图中去除地球曲率、地形因素等产生的相位，从而得到地表形变信息。它是迄今为止独一无二的基于面观测的形变遥感监测手段（刘国祥，2006）。与常规方法相比，D-InSAR技术监测地表形变具有大面积、连续、快速、准确的优势，同时观测无需布设地面控制点，从而降低了对危险地区（比如火山地区）观测的危险性（马超等，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方向上的同震形变场。本文采用二轨法利用 ENVI SARscape模块，基于伊朗巴姆地震前后的两景SAR影像和SRTM DEM数据进行差分干涉处理，并利用雷达波束的局部入射角得到垂直方向上的同震形变场。最后结合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波段的两景ASAR Image SLC数据，数据的部分参数见表1。利用重采样后的SRTM DEM数据进行地形相位的去除，格网大小为25m×25m，UTM投影40分带，参考椭球体采用WGS84。研究区SAR影像见图1。从图1可以看出，数据覆盖范围大部分地方地形平坦，西南角有部分覆盖了山区。

表1 ENVISAT ASAR数据部分参数Table 1 Some parameters of ENVISAT ASAR data

2 D-InSAR技术原理

Zebker等（1994）、舒宁（2003）、王超等（2002）对D-InSAR技术测量形变的基本原理已经有了非常详细论述，本文不再赘述。能够进行干涉处理的 SAR影像必须处理成 SLC（Single Look Complex）格式，SLC影像上的每个象元可表示为：

式中，A (m, n) 是影像的强度信息即振幅；φ (m, n) 是相位信息。

干涉处理就是将两景配准后的影像逐象元进行复共轭相乘，记两景影像象元分别是（何秀凤等，2012)：

逐象元进行复共轭相乘后得到：

因此，干涉相位为：

这也是最后干涉图上每个象元的干涉相位信息，φ (m, n) 的取值范围为 (-π, π)。但是，地球表面是一个具有一定曲率的球面，微波信号在大气中传播时还会受到大气的影响，而且雷达影像有一定的噪音，所以一般认为干涉相位由以下几部分组成（刘国祥等，2002）：

式中，φflat_Earth是地球曲率产生的相位信息，也就是高度不变的平地在干涉图中表示出来的干涉条纹随距离向和方位向的变化而呈周期性变化的现象（王超等，2002）；φtopo是地形产生的相位信息，是由于地形起伏造成的，可用于建立DEM；φdisp是由地表形变产生的相位信息，反映了地表形变信息；φpath是雷达信号在传播过程中受到大气影响产生的相位信息；φnoise是雷达影像噪声产生的相位信息；n·2π是相位的整周数，需要通过相位解缠获得。

3 实验验证

3.1 实验步骤

采用D-InSAR技术进行地表形变监测主要包括以下几个步骤：基线估计、影像配准、主干涉图生成、平地相位及地形相位的消除和相位解缠形成形变图。数据处理流程见图 2（邵芸等，2002）。

图2 二轨法差分干涉测量流程图（邵芸等，2002）Fig. 2 The data processing procedure for two-pass differential radar interferometry (from Shao et al., 2002)

①基线估计。干涉像对的基线大小直接影响着干涉图的质量，基线过大则无法产生干涉。因此在进行干涉处理前，首先要进行基线估计，以确定像对能否产生干涉。一般认为ENVISAT的C波段数据临界基线大约为1100m，最佳基线应小于200m，而用于地面变化监测的基线大小应小于 70m（舒宁，2003），本文所用数据的垂直基线估计结果为 9.1m，可以产生干涉。②影像配准。进行干涉的两幅SAR影像之间必须进行精确地配准，一般认为配准精度应达到亚像元级（Ito等，2003）。本文利用SAR数据自带的卫星轨道参数，进行自动精确配准（Klees等，1999）。配准精度可以通过相干系数图（图 3）间接地进行检查。③主干涉图生成。配准后的主、副影像相应象元进行复共轭相乘，得到主干涉图（图4）。④平地相位及地形相位的去除。在主干涉图上去除平地相位和地形相位以得到差分干涉图（图5），其中，平地相位可通过对干涉条纹乘以复相位函数来去除（王超等，2002）。对于地形相位的消除，本文则是利用该地区的SRTM DEM数据和进行干涉的两景SAR影像的天线姿态参数生成地形的模拟干涉图，从而得到地形相位。⑤相位解缠形成形变图。最后对差分干涉图采用最小费用流（Minimum Cost Flow, MCF）方法进行相位解缠得到形变图（图6）。

3.2 实验结果

图4是两景SAR影像形成的主干涉图，干涉条纹从黑色到白色的一个颜色周期代表了相位差从-π到π的变化，可以看出在巴姆城市的东侧形成2个明显的花瓣状干涉条纹，周围的干涉条纹主要是由于地球曲率和地形造成的。图5是去除了平地相位和地形相位的差分干涉图，同样，从黑色到白色的一个颜色周期代表了相位差从-π到π的变化，可以看到在巴姆城市周围形成了4个花瓣状（或蝴蝶状）的干涉条纹，东侧的条纹多且密集，而西侧的条纹少且稀疏。北侧两组条纹的颜色变化从里到外是由白到黑，而南侧条纹的颜色变化从里到外是由黑到白，说明巴姆城市南北两侧的形变方向不同，这一点从最后得到的同震形变场也可以看到。

在干涉图上可以看到，一些地区（巴姆市及其东南侧）干涉条纹十分杂乱，这种现象称为去相关现象，发生去相关现象的斑块与城市区域和植被区域吻合得很好，所以可以初步推断去相关现象的发生是由植被的变化和城市建筑的毁坏造成的。另外在巴姆市的西侧有条带状的去相关区域，通过光学遥感影像可以看出这里是河流，所以可以初步认定这里的去相关现象是由河流造成的。发生破坏的城市、植被和河流的散射特性极不稳定，在两景影像获取的70天之内发生了极大地变化，所以这些地方的干涉效果比较差，去相关现象十分严重。去相关现象也可以从两景影像的相干系数图中看出，两景影像成像之间散射特性的改变会降低两景影像之间的相关性，去相关现象越严重的地方，相关系数越低，如图3中的黑色地方（相干系数接近于0）。相干系数图也可以用来评价干涉质量（舒宁，2003）。

舒宁（2003）提出变形监测采用雷达干涉测量方法时必须满足2个条件：一是2次观测期间的变化梯度必须小于一定阈值；二是2次观测期间，一个象元内的所有地物散射体的特性应具有相似性，特别是表面位置变化的误差应只是雷达波长的10%—20%的量级。去相关效应与干涉的两景SAR影像之间发生变化的程度等有关，大气影响、人类活动等都有可能引起去相关效应（Klees等，1999)。因此，雷达干涉测量监测地表形变的方法适用于地表变化在一定范围内、干燥、少植被的地区。

图3 相干系数图Fig. 3 Coherence coefficient image

图4 主干涉图Fig. 4 The major interferogram

图5 差分干涉图Fig. 5 The differettntial interferogram

图6 同震形变图（单位：cm）Fig. 6 Coseismic deformation map (unit: cm)

3.3 地表同震形变场分析

从干涉图直接解缠得到的形变是沿着视线方向（Line of Sight, LOS）的，为了将形变转换到垂直方向上，假设LOS方向的形变量是垂直方向形变量的分量，根据公式（6）转换：

式中，D⊥是垂直形变量；DLOS是LOS方向的形变量；α是雷达波束的局部入射角。

差分干涉图经过相位解缠和形变值的转换得到垂直方向的同震形变场，如图6所示。图上的线条是形变等值线。震中的位置距离形变最大点距离分别是15km（最大沉降点）和12km（最大抬升点），而震中附近的形变只有1—2cm。

从图6所示的形变图可以看出，在远离巴姆城市的地区形变量基本在0cm左右，说明在长达70天的时间里，巴姆地区除了地震造成了严重形变外，几乎没有大的形变发生，所以得到的形变可以认为基本是由地震造成的。地震造成的地表形变主要发生在巴姆城市的东南侧和东北侧。其中最大的地形抬升发生在巴姆城市东南侧，如图中白点所在位置，形变量达到 33.0cm，这点的高程约为 1042m，最大的地形沉降发生在巴姆城市东北侧，如图中黑点所在位置，沉降量达到 18.3cm，这点的高程约为 1021m。在巴姆城市西南侧有少量的上升，上升量在2cm左右；巴姆城市西北侧有少量的下沉，沉降量为2—4cm左右。在图中黑色粗线所示位置可以明显地看到一条断层特征。断层呈现南北走向，长约10km。

在形变图上过最大形变点分别作2条东西走向和南北走向的剖面线（剖面线A、B、C、D），如图7所示。

图7 形变剖面图Fig. 7 The profiles (A, B, C and D) of the deformation map

图7 a和b是过沉降最大点的剖面线，图7c和d是过抬升最大点的剖面线，其中，图7a、c为东走向剖面图，图7b、d为南北走向剖面图。从这4个剖面图可以看出，巴姆地区的形变基本表现为巴姆城市北部地区的地块有所沉降，而巴姆城市南部地区的地块被抬升。

4 结论与讨论

本文利用D-InSAR技术对伊朗巴姆地区2003年12月26日6.7级地震造成的地表形变进行了研究，得到了该地区垂直向的三维同震形变场（图6）。

（1）巴姆地震造成的地表形变主要分布在巴姆市的东侧，垂直形变量在-18.3cm到33.0cm之间。这个结果与李毅等（2008）、王志勇等（2008）得到的结论接近，但是存在几毫米到几厘米的误差。主要是因为他们得到的是雷达视线方向LOS上的形变，而本文将这个形变量转换成垂直方向的形变。从三维同震形变场的图上可以看到巴姆地区的形变主要表现为：北部地块发生沉降，最大沉降地点位于巴姆城市东北侧，距震中约 15km；南部地块被抬升，最大抬升地点位于巴姆城市东南侧，距震中约12km。震中附近形变量只有1—2cm。另外，地震还在巴姆城市南侧造成了长约10km、南北走向的断层。

（2）在巴姆地区的某些区域发生了去相关效应，这些现象主要发生在巴姆城市、河流和植被茂盛的区域。这主要是河流和植被的散射特性不稳定，在相隔2个多月的时间里，植被和河流散射特性已经发生了较大的变化，从而产生了去相关现象。巴姆城市地区则由于建筑物被地震严重地损毁，散射特性发生了极大变化，所以产生了去相关效应。一般认为，沙漠地区比密林地区的干涉条件好，干燥环境比潮湿环境的干涉效果好，波长较长的雷达图像会得到较好效果的干涉图（舒宁，2003）。

（3）巴姆地区差分干涉图是利用精度较低的SRTM DEM数据去除地形相位信息的，这样造成的地形误差较大，在今后的研究中，在提高地形测量精度方面还需要进一步努力。另外，由于D-InSAR技术对地表散射体散射特性的稳定性要求较高，所以对于植被丰富，气候潮湿的地区，就很难产生较好的干涉效果，因此在去相关问题上还需进一步的研究。

利用雷达影像差分干涉测量得到的形变图仅仅包含了LOS向的形变信息（张红，2002），即使得到了垂直的形变量，也是基于LOS方向形变是垂直形变量的分量的假设。要得到地表形变的真实方向及大小，需要两个影像对进行D-InSAR形变分析，从而合成出真实的形变方向及大小。另外，为了对地表形变场认识得更加全面充分，还需要结合GPS、GIS等手段建立形变模型，从而确定研究区的形变机制。如果能够通过技术进步提高雷达干涉测量的精度，并降低观测成本，建立长期的观测机制，周期性的获取高质量的SAR影像用来观测地表形变量，则可以对地震、滑坡、火山爆发、地面沉降等灾害的预报提供有价值的观测数据，从而在灾难发生前后采取有效措施以降低损失，真正让这项技术发挥价值。

查显杰，2007. InSAR形变测量及其在地震学中应用的研究. 合肥：中国科学技术大学.

何秀凤，何敏，2012. InSAR对地观测数据处理方法. 北京：科学出版社.

季灵运，许建东，2009. 利用D-InSAR和AZO技术获取Bam地震同震三维形变场. 大地测量地球动力学，29（6）40—44.

李毅，柳林涛，薛怀平等，2008. 用D-InSAR研究巴姆地震形变场. 大地测量与地球动力学，28（1）：32—35.

刘大杰，胡丛玮，1999. 应用GPS监测城市地表形变的初步分析. 地壳形变与地震，19（1）：37—42.

刘国祥，2006. 利用雷达干涉技术监测区域地表形变. 北京：测绘出版社.

刘国祥，丁晓利，李志伟等，2002. ERS卫星雷达干涉测量：1999年台湾集集大地震震前和同震地表位移. 地球物理学报，45（增刊）：165—174.

凌勇，曾祺明，罗扬，赵永红，2006. 巴姆地震变形场和应力场：Ⅰ. 用差分干涉雷达和Okada方法求解. 岩石学报，22（9）：2367—2374.

罗扬，施旭，赵永红，2006. 巴姆地震变形场和应力场：Ⅱ. 用FEPG有限元方法求解. 岩石学报，22（9）：2375—2380.

马超，单新建，2004. 星载合成孔径雷达差分干涉测量（D-InSAR）技术在形变监测中的应用概述. 中国地震，20（4）：410—418.

单新建，柳稼航，马超，2004. 2001年昆仑山口西8.1级地震同震形变场特征的初步分析. 地震学报，26（5）：474—480.

邵芸，谭衢霖，刘浩等，2002. 雷达差分干涉测量技术及其应用研究——以1997年西藏玛尼7.9级地震区域形变测量为例. 地球物理学报，45（增刊）：205—213.

舒宁，2003. 雷达影像干涉测量原理. 武汉：武汉大学出版社.

孙建宝，梁芳，沈正康等，2008. 汶川 MS8.0地震 InSAR形变观测及初步分析. 地震地质，30（3）：789—795.

王超，张红，刘智，2002. 星载合成孔径雷达干涉测量. 北京：科学出版社.

王志勇，张继贤，张永红等，2008. 伊朗巴姆 6.5级地震同震形变场的获取与解译. 地震研究，31（1）：70—76.

张保淑，2009. 科技之谜：追寻巴姆古城堡. http://www.people.com.cn/GB/keji/1058 /2284182.html.

张红，王超，刘智，2002. 获取张北地震同震形变场的差分干涉测量技术. 中国图象图形学报，5A（6）：497—500.

Bam (Iran), December 2003 [2009-04-26]. http://earth.esa.int/ew/earthquakes/Bam_Iran_Earth- quake.

Ito Y., Hosokawa M., 2003. A degree-of-damage estimation model of earthquake damage using interferometric SAR data. Electrical Engineering in Japan, 143 (3): 49—57.

Klees R., Massonnet D., 1999. Deformation measurements using SAR interferometry: potenital and limitation.Geologie en Mijnbouw, 77 (2): 161—176.

Zebker H.A., Rousen P., 1994. On the derivation of coseismic displacement fields using differential radar interferometry: The Landers earthquake. Geoscience and Remote Sensing Symposium, 1994. IGARSS '94.Surface and Atmospheric Remote Sensing: Technologies, Data Analysis and Interpretation. 1: 286—288.