D-InSAR技术在库区滑坡监测上的实例分析

王桂杰，谢谟文，柴小庆，王立伟，董晨曦

(1.北京科技大学土木与环境工程学院，北京 100083;2.中国地质环境监测院，北京 100081)

滑坡是一种常见的自然地质灾害，每年都会导致数以千计的受害者和数百亿元的经济损失。对滑坡的监测、预报和早期防治，一直是世界各国广泛关注和研究的问题。然而，传统的滑坡监测方法，已经不能满足广域大范围滑坡监测[1]。随着遥感技术的发展，合成孔径雷达差分干涉D-InSAR技术，正好能够解决连续大面积上非常小的地面活动监测，而且具有高精度、高分辨率、大面积、全天候、监测成本低，以及能够监测人员无法进入的区域等优点。因此，已经被广泛应用在滑坡灾害调查的科学研究和风险评估上[2-10]。

1 D-InSAR技术基本原理和处理流程

1.1 D-InSAR技术的基本原理

合成孔径雷达差分干涉测量D-InSAR技术，是通过对SAR数据的处理，将同一场景与地形和形变有关的两幅干涉SAR图像的相位差，转变为地表微小形变值的技术。其中，两幅干涉SAR图像中的一幅，是由形变前获得的两帧单视复数(SLC)SAR数据产生的，仅仅包含地形相位项；另一幅是由形变前获取的一帧SLC和形变后获取的另一帧SLC产生的，既包含地形相位项，又包含形变相位项。利用已知的卫星轨道数据和地形数据(如已知的数字高程DEM)，能够很容易去除地形相位项，差分后余下的相位项即与地面形变(沿卫星视线方向)成正比。其测量精度为其所用电磁波波长的量级，因此在微波频率范围内，可达到厘米级的精度[6]。

通常，D-InSAR技术有三种方法，即，已知DEM的两轨法、三轨法和四轨法。它们的几何原理基本相同。与三轨法比较，两轨法中其地形数据对是从已知的DEM获得的，而四轨法的地形数据对是从两个不同的位置分别获得的。三轨法的基本原理如图1所示[11,12]。

图1 三轨法差分干涉测量成像几何示意图

A1、A2、A3是卫星三次对同一地区成像的位置。由图1中几何关系及B1<

(1)

式(1)表明，用干涉测量法所得到的相位差与视线方向的基线分量成正比。这里设定在A1处获得第一幅影像且为主影像，假设在地表未发生形变前在A2处获取第二幅影像，所以第二幅影像与A1处的主影像形成的干涉SAR图像，其干涉相位仅包含地形信息，即相位差可表示为：

(2)

式中：星载重复轨道ρ=2，λ为微波波长。假设发生形变后，在A3处获取了第三幅影像，所以第三幅影像与主影像形成的干涉SAR图像的相位，既包含地形信息，又包含地表形变信息。且由于获得的影像间要求基线足够小，所以可近似看作θ不变，即相位差可表示为：

(3)

由式(2)和式(3)得：

(4)

假设ΔRd为视线向形变量，则由式(2)、式(3)、式(4)推得由视线向形变量引起的干涉条纹图相位差Δφd可表示为：

(5)

式(5)中左边的各量，可由干涉条纹图的相位和轨道参数计算得到，进而可确定影像每点的视线向形变量ΔRd，分解后得到水平形变量和垂直形变量等。

1.2 D-InSAR技术获得形变的数据处理流程

D-InSAR技术三轨法和两轨法的数据处理流程，分别见图2、图3。在图2中，三帧单视复数数据(SLC) 通过精确地配准和多视产生两幅干涉图。由地形数据对产生的干涉图，仅仅包含地形相位项；由地形/形变数据对产生的干涉图，既包含地形相位项，又包含形变相位项。地形数据对和地形/形变数据对必须被正确的选取。对于地形数据对，要求有较短的时间基线；而地形/形变数据对，要求具有较短的空间基线。然后，将生成的两幅干涉图进行相位差分，消除地形相位项，得到仅含形变相位项的差分干涉条纹图。对差分干涉条纹图进行滤波处理，消除相位噪声，同时衡量干涉质量的相干系数图也被产生。差分干涉图中的差分相位仅仅是0～2π主值范围内的，即以2π为模数的，因此相位解缠是D-InSAR技术中较重要的步骤，通过相位解缠得到绝对的形变相位值。最后，将解缠后的相位值转换为形变值，直接进行地理编码得到地表的微小形变值。

在图3已知DEM的两轨法中，选取两帧单视复数数据(SLC)作为地形/形变数据对，产生地形/形变干涉图，已知DEM被合成地形相位项。然后，通过相位差分、干涉滤波、相干系数计算、相位解缠和地理编码等，最后产生地表微小形变值。处理过程与三轨法是基本相同的。

图2 D-InSAR技术三轨法数据处理流程

图3 D-InSAR两轨法数据处理流程

2 D-InSA技术在滑坡辨识上的试验研究

2.1 研究区域和应用的数据

研究区域位于云南省禄劝县和四川省会东县交界的乌东德水电站库区。图4为研究区域背景图，区域内标记出正处活动状态的滑坡，黑色线圈标记的区域是本研究的20km×20km的主研究区域。图5 是主研究区域详图，其为Sport5与1∶2000的高精度航片融合的立体影像。图6为主研究区域上一正处活动状态的L1R-6号滑坡详图，并有一些GPS监测点位于该滑坡上。

图4 研究区域背景图

图5 主研究区域详图

图6 L1R-6号滑坡详图

在我们的研究中，有5景ALOS卫星PALSAR传感器的SAR数据被获取，波长是23cm(L波段)，最长时间跨度是368d，ALOS卫星 PALSAR的数据参数见表1。

在三轨法中，将2008年1月12日获取的影像作为主影像，将2008年1月12日与2008年2月13日获取的两幅影像组成的数据对作为地形数据对AT，将2008年1月12日与2008年4月13日获取的两幅影像组成的数据对作为地形/形变数据对AD1。

在两轨法中，将2007年7月12日获取的影像作为主影像，将2007年7月12日与2008年7月14日获取的两幅影像组成的数据对作为地形/形变数据对AD2。将2.5m分辨率与ASTER30m分辨率镶嵌后的DEM，作为已知DEM合成地形相位项。

表1 所用ALOS PALSAR卫星的数据参数

2.2 微小形变位移图的产生

通过对获取的5景单视复数SAR数据的处理和分析，获得了研究区域内滑坡的形变位移图。首先，三组单视复数图像对(AT1、 AD1、AD2)经精确地配准和干涉处理，产生三幅干涉图(AT1、 AD1、AD2)。如图7所示，我们能够看到，时间基线和空间基线的影响是非常重要的，AT1的干涉条纹最清晰。干涉图AT1是由三轨法中的地形数据对产生的，AD1、AD2分别是由三轨法及两轨法中的地形/形变数据对产生的。另外，在两轨法中，由于地形数据是由已知镶嵌的DEM提供的，因此，可不用直接形成干涉图。然后，将干涉图AD1与干涉图AT1、干涉图AD2与镶嵌DEM合成的相位项分别进行差分处理，得到AD1与AT1以及AD2与合成地形相位项差分后仅含地形相位的差分干涉图，并进行滤波去噪处理，得到滤波后的差分干涉图，如图8所示。相干系数图也被产生，如图9所示。

最后，通过相位解缠、地理编码等，获得研究区域上两个研究期间内详细的滑动位移图，将其转换为沿垂直方向和水平方向的位移图来表示，如图10、图11所示。在整景数据区域中，四边形区域是主研究区域，不规则多边形区域是2.5m分辨率DEM镶嵌的区域，图中沉降区域和抬升区域都采用逐级的灰度图来表示，并用正、负号区分，负值区域表示的是沉降区域及其沉降值，正值区域表示的是抬升区域及其抬升值。由于滑坡的滑动主要伴随着沉降，因此在主研究区域中，我们主要研究沉降区域。

在主研究区域的垂直滑动位移图中，抬升区域用小于0.02m和大于等于0.02m的两个区域来表示；在图10中，沉降区域用四个不同的灰级来表示四个不同的位移值区间：0～0.05m、0.05～0.1m, 0.1～0.2m、>0.2m。在图11中，沉降区域用五个不同的灰级来表示五个不同的位移值区间：0～0.05m、 0.05～0.1m、 0.1～0.2m、 0.2～0.5m、 >0.5m。

图7 初始干涉条纹图

图8 滤波后的差分干涉图

图9 相干系数图

在主研究区域的水平滑动位移图中，按水平位移与垂直位移所成比例的倍数，将抬升区域的水平位移用小于0.774m和大于等于0.744m的两个区域来表示。在图10中，沉降区域的水平位移值用四个不同的灰级来表示四个不同的位移值区间：0～0.194m、 0.194～0.387m、 0.388～0.774m、≥0.775m。在图11中，沉降区域的水平位移值，用五个不同的灰级来表示五个不同的位移值区间：0～0.194m、 0.194～0.387m、 0.388～0.774m、 0.775～1.938m、≥1.938m。其中，水平位移是指北南向的水平位移，图中负号表示的是沉降区域的水平位移，并不代表方向。从图10和图11中，我们能够清晰的辨识出研究区域的位移状态，将垂直和水平滑动位移较大的区域判定为潜在的滑坡区域或滑坡的危险区域，并在图10和图11中分别用灰色线圈标记出。

2.3 D-InSAR监测结果及分析

由图10、图11可知，抬升区域垂直位移值绝大部分都处在小于20cm的区间上。对于沉降区域，在92d的形变位移图10中，其垂直位移最大为25.3cm；在368d的形变位移图11中，其垂直位移值最大为96.6cm。因此，由于PALSAR数据所用电磁波长为L波段23cm，再考虑到各种可能误差, 能够得出两个时间期间抬升区域的位移值几乎都处在误差影响的范围内。因此可知，92d (2008-01-12～2008-04-13)研究时间内，滑动位移是很小的。我们将垂直滑动位移大于20cm(水平位移77.4cm)的区域用灰色线圈标记出(图中除黑色线圈LIR-6号滑坡标记外的其他两个线圈)，作为潜在的滑坡或滑坡灾害区域。在368d (2007-01-12～2008-07-14)的形变位移图中，我们用灰色线圈标记出垂直位移大于50cm(水平位移193.8cm)的区域(图中最大的不规则线圈标记的那个区域)，作为潜在的和较危险的滑坡区域，同时垂直位移值大于20cm(水平位移77.4cm) 也能被较好的辨识出。比较两个时间段的形变位移图10和图11，我们发现，在92d的时间期内发生较大位移的区域，在368d的形变图中没有发生较大的改变。这表明在92d的形变位移图中，发生位移较大的两个区域，可能是物质的季节性流动而不是潜在的滑坡。最后，经光学影像确认，这两个区域位于金沙江右岸一个较大泥石流沟的沟口和中部汇水区。因此，在92d位移图中，这两个区域位移的改变，最可能的原因就是雨水导致沟中物质迁移所引起的。在368d(2007/07/12～2008/07/14)的形变位移图中，发生较大位移的不规则线圈的区域，被确认为是几个较大泥石流沟和三个正处活动状态滑坡( L6R-5, L9R-4, L10R-3)的汇集区。因此，这个不规则线圈圈定的区域应该被作为重点研究和控制区域。

图10 研究区域92d的形变状态图(单位：m)

图11 研究区域368d的形变状态(单位：m)

另外，在368d的形变位移图中，我们也用灰色线圈标记出已发生过滑坡的区域(并在其上标有滑坡编号)，用较细的灰色线圈标记出现场监测正发生滑动的区域。从图11中我们看到，这些区域都位于沉降区域上垂直位移值大于20cm的区域上。因此，能够得出D-InSAR监测结果与现场监测结果十分一致的结论。

3 L1R-6号滑坡分析及与GPS监测结果的比较

3.1 L1R-6号滑坡的详细分析

为了对研究区域上正处活动状态并对水库正常运营有重要影响的L1R-6号滑坡进行进一步研究。两个时间期间内，L1R-6号滑坡的D-InSAR监测结果如图12和图13所示，位移图中每一点的位移值都能够被读出。但由于位移值数据庞大，在图12和图13中用灰色分级图来表示不同的位移值区间。通过位移图，不同区域的滑动速率能够清晰地被辨识出。在位移图12 和图13中，可以看到No.L1R-6号滑坡被划分为四个区域：Ⅰ、Ⅱ、Ⅲ、Ⅳ[13,14]。在92d的位移图12中，我们能够看到Ⅰ区和Ⅳ区的位移值是较小的(垂直小于0.1m,水平小于0.387m),在Ⅱ区和Ⅲ区的部分区域位移值是较大的(垂直是0.1～0.133m，水平是0.387～0.514m)。在368d的位移图13中，Ⅱ区部分区域的位移值和Ⅳ区一小部分的位移值是较大的(垂直>0.2m,水平>0.775m)，另外的区域，垂直位移值都小于0.2m(水平0.775m)。L1R-6号滑坡上白色字体标记的白色点为GPS监测点。

图12 L1R-6号滑坡92 d的形变图(单位：m)

图13 L1R-6号滑坡368 d的形变图(单位：m)

3.2 与GPS监测结果的比较

为了验证D-InSAR技术的实用性和有效性，以及检测单点监测精度, 对D-InSAR在GPS监测点上获得的92d和368d的形变位移值与GPS监测的位移值进行比较。 在与92dD-InSAR监测结果的比较中，GPS的监测时间在Ⅱ区上是从2008年1月11日到2008年4月19日，在Ⅰ、Ⅲ、Ⅳ区上的监测时间是从2008年1月12日到2008年4月18日。在与368dD-InSAR监测结果的比较中，GPS的监测时间在Ⅱ区上是从2007年6月23日到2008年7月15日，在Ⅰ、Ⅲ、Ⅳ区上的监测时间是从2007年6月28日到2008年7月14日。其中，点名为AL01D、AL02D、AL03D、TP44、TP45等5个监测点位于Ⅰ区，TP01- TP15以及AL02C和AL03C等17个监测点位于Ⅱ区， AL01A、AL02A、AL03A、TN10、TN11、AL02B、AL03B等7个监测点位于Ⅲ区, AL01B与AL01C两个监测点位于Ⅳ区。垂直位移值是以向下为正，水平位移值是以向南为正，相反方向为负号。比较结果如表2所示。

从表2我们可以得出，尽管D-InSAR监测结果的单点精度并不完美，D-InSAR处在厘米的量级上，而GPS监测结果为毫米的量级，但其整体形变和运动趋势是基本一致的。两个时期的两种调查结果都表明，Ⅱ区的形变是最大的，必须被作为重点灾害区进行进一步的调查, 其他区域的形变是相对较小的。由于误差的影响，尽管在92d形变位移图的Ⅲ区是略有不同的,但在其他区域以及368d的形变图中，两种技术的监测结果与实际情况是完全相符的。

表2 D-InSAR与GPS的结果比较数据表

4 结 论

本文通过D-InSAR技术，利用PALSAR传感器获取的5景SAR数据，对金沙江下游乌东德水电站库区内大面积区域上的滑坡进行辨识和监测，获得了不同区域的滑动状态，以及研究时间期间内研究区域上详细的位移图，潜在的滑坡滑动区域和滑坡危险区域被确定。同时，L1R-6号滑坡被详细的研究，并与GPS监测结果相比较，得到如下结论：通过PALSAR数据的D-InSAR技术，能够辨识出潜在的滑坡区域；D-InSAR与GPS监测结果的位移趋势是基本一致的，证明了D-InSAR技术在滑坡监测上的有效性及可行性。

因此，随着遥感技术、先进地球探测技术以及数字处理技术的发展，D-InSAR技术必将成为滑坡早期辨识和监测的重要技术之一。

[1] Yao G, and Mu J. D-InSAR Technique for land subsidence monitoring. Earth Science Frontiers, 2008,15(4): 239-243.

[2] Fruneau B, Achace J, Delacourt C. Observation and modeling of the Sant-Etienne-de Tine′e landslide using SAR interferometry. Tectonophysics , 1996, 265:181-190.

[3] Achache J, Fruneau B, Delacourt C. Applicability of SAR interferometry for operational monitoring of landslides. In: Proceedings of the Second International Workshop on ERS Applications. London, 6-8 December 1995, 165-168.

[4] Vietmeier J, Wagner W, Dikau R.. Monitoring moderate slope movements (landslides) in the southern French Alps using differential SAR interferometry. In: Proceedings of 2nd International Workshop on ERS SAR Interferometry.Liège, Belgium, 10-12 November. 1999.

[5] Rizzo V., Tesauro M. 2000.SAR interferometry and field data of Randazzo landslide (Eastern Sicily, Italy). Physics and Chemistry of the Earth , 25(9): 771-780.

[6] Paolo B, Mario C, Giorgio F, et al. Use of differential SAR interferometry in monitoring and modelling large slope instability at Maratea (Basilicata, Italy). Engineering Geology, 2003, 68: 31-51.

[7] Squarzoni C, Delacourt C, Allemand P. Nine years of spatial and temporal evolution of the La Valette landslide observed by SAR interferometry. Engineering Geology, 2003, 68(1-2): 53-66

[8] Colesanti C, Wasowski J. Investigating landslides with space-borne Synthetic Aperture Radar (SAR) Interferometry. Engineering Geology ,2006, 88: 173-199.

[9] Singhroy V, Couture R, Molch K, et al. InSAR monitoring of post-landslide activity. International Geoscience and Remote Sensing Symposium (IGARSS). Denver, Colorado: July 31-Aug.4, 2006, 1635-1638.

[10] Singhroy V, Alasset P-J, Couture R, et al. InSAR monitoring of landslides on permafrost terrain in Canada. International Geoscience and Remote Sensing Symposium (IGARSS) , Boston, Massachusetts, USA: July 6-11, 2008,2451-2454.

[11] 廖明生，林珲.合成孔径雷达干涉原理和信号处理[M]. 北京：测绘出版社.

[12] 王桂杰，谢谟文，邱骋，江崎哲郎. D-INSAR 技术在大范围滑坡监测中的应用[J].岩土力学, 2010,31(4):1337-1344.

[13] 李会中，王团乐，段伟峰，等.金坪子滑坡形成机制分析与河段河谷地貌演化地质研究[J].长江科学院院报，2006，23(4)：17-22.

[14] 谢谟文，江崎哲郎，邱骋，等.空间三维滑坡敏感性分区工具及其应用[J].地学前缘, 2007,14(6):73-84.

[15] Ferretti A, Prati C, Rocca F. Permanent Scatterers InSAR interferometry. IEEE Transactions on Geoscience and Remote Sensing, 2001,39(1): 8-20.

[16] Gili J- A, Corominas J, Rius J. Using Global Positioning System techniques in landslide monitoring.Engineerin. Geology, 2000,55:167-192.