差分干涉合成孔径雷达技术在广域滑坡动态辨识上的实验研究

D010.13374斤.is0053x.20l.02.021 第33卷第2期 北京科技大学学报 Vo133 No 2 2011年2月 Journal ofUniversity of Science and Technobgy Bejjing Feb 2011 差分干涉合成孔径雷达技术在广域滑坡动态辨识上的 实验研究 王桂杰12》 谢谟文1) 邱骋”吴维伦D黄洁慧1》 1)北京科技大学土木与环境工程学院北京1000832)中国地质环境监测院北京100081 通信作者，Ema到movenx@used!m 摘要通过AIOS卫星PALSAR传感器获得的三景合成孔径雷达(SAR)数据，利用合成孔径雷达差分干涉(D-SAR技术 对金沙江下游乌东德水电站库区内的滑玻活动进行了详细的研究。获得研究区域内地表高精度形变位移值.将其按滑动速率 及位移大小进行分级显示，清晰地表明了研究区内不同区域的地面活动位移状态，辨识出研究区内可能发生滑动和发生滑动 较大的区域，从而确定了滑坡活动的风险区.对库区内正处活动状态的山R6号滑坡进行了详细的研究.结果表明：DS4R 技术与全球定位系统GS在监测整体形变和运动趋势上基本一致.对DSAR结果存在的单点误差进行了分析，提出了 DSAR与GPS湘融合的栅格函数误差插值消减法提高了D-mSAR技术的监测精度. 关键词合成孔径雷达(SAR,滑坡：遥感：监测：全球定位系统(GPS 分类号TP79 Experm ent research ofD-InSAR technique on identify ing landslidemoving in a w ide area WANGGui j 2)XE Mowen),QU Cheng).WUWei-un,HUANG Jiehup 1)SchoolofCivil and Enviramenl Engineerng University of Science and Technokgy Beijing Beijng 100083 China 2)China Institute ofGeenvirommentMonitorng Beijng 100081 China Correspand ng author Email movenxa@usth edu c ABSTRACT Based on hree frame syn hetic aperture radar(SAR)dana derived fron AlOS satellite sensor PALSAR a three pass d ifferentialSAR nterferomnetry (DInSAR)technique was used p anayze lndslide activities in the Wudongde Hyd opower Reservoir area at pwer Jnsha River and high accuracy gound disp acement values were acquired by this mehod The classificatpn of these ground displacen ent vales by sliding velocity and d spacement cleary iljustrated the groind activity defomation states of varjous zones in the sud ed area were ceary ob taned detem ined sone poren talmovng landslides and active lands lides and dentified he danger ous zaes of lndslides For an active landslide numbered No 1R-6 he depmation detected by D hSAR anaysis shows an accord ance tendency with the one by gpbal pos itponing system (GPS monitoring In the end single ponterrors which existed in he result ofD-InSAR technique were analyzed and a gril fmction error n terpoltion metol ofD hSAR and GPSwas prposed p mpove he ma itoring accuracy ofD InSAR technique KEY WORDS syn thetic aperture radar SAR)ndslides remote sensng mon ioring global positian ng system (GPS) 滑坡是一种常见的自然地质灾害，给国家建设 全球定位系统(GPS)、时域反射法(tme donain 和人民生命财产造成巨大损失，对滑坡的监测、预 e作ctmetry TDR技术和分布式光纤等电子远程 报和早期防治一直是世界各国广泛关注和研究的问 监测技术等.这些技术的优点是单点监测精度高， 题.传统上对滑坡的监测和预警主要采用钴孔倾斜 但无法满足人员难以进入的区域以及广域大面积区 仪、位移计和自动伸缩计等现场仪器观测，以及利用 域的地面活动监测需要.随着合成孔径雷达干涉 收稿日期：2010-04-08 基金项目：国家自然科学基金资助项目(NQ4972229)

第 33卷 第 2期 2011年 2月 北 京 科 技 大 学 学 报 JournalofUniversityofScienceandTechnologyBeijing Vol.33 No.2 Feb.2011 差分干涉合成孔径雷达技术在广域滑坡动态辨识上的 实验研究 王桂杰 1, 2) 谢谟文 1) 邱 骋 1) 吴维伦 1) 黄洁慧 1) 1)北京科技大学土木与环境工程学院, 北京 100083 2)中国地质环境监测院, 北京 100081 通信作者, E-mail:mowenxie@ustb.edu.cn 摘 要 通过 ALOS卫星 PALSAR传感器获得的三景合成孔径雷达(SAR)数据, 利用合成孔径雷达差分干涉(D-InSAR)技术 对金沙江下游乌东德水电站库区内的滑坡活动进行了详细的研究, 获得研究区域内地表高精度形变位移值.将其按滑动速率 及位移大小进行分级显示, 清晰地表明了研究区内不同区域的地面活动位移状态, 辨识出研究区内可能发生滑动和发生滑动 较大的区域 , 从而确定了滑坡活动的风险区.对库区内正处活动状态的 L1R--6号滑坡进行了详细的研究.结果表明:D-InSAR 技术与全球定位系统 GPS在监测整体形变和运动趋势上基本一致.对 D-InSAR结果存在的单点误差进行了分析, 提出了 D-InSAR与 GPS相融合的栅格函数误差插值消减法, 提高了 D-InSAR技术的监测精度. 关键词 合成孔径雷达(SAR);滑坡;遥感;监测;全球定位系统(GPS) 分类号 TP79 ExperimentresearchofD-InSARtechniqueonidentifyinglandslidemovingina widearea WANGGui-jie1, 2) , XIEMo-wen1) , QIUCheng1) , WUWei-lun1) , HUANGJie-hui1) 1)SchoolofCivilandEnvironmentalEngineering, UniversityofScienceandTechnologyBeijing, Beijing100083, China 2)ChinaInstituteofGeo-environmentMonitoring, Beijing100081, China Correspondingauthor, E-mail:mowenxie@ustb.edu.cn ABSTRACT Basedonthree-framesyntheticapertureradar(SAR)dataderivedfromALOSsatellitesensorPALSAR, athree-pass differentialSARinterferometry(D-InSAR)techniquewasusedtoanalyzelandslideactivitiesintheWudongdeHydropowerReservoir areaatlowerJinshaRiver, andhighaccuracygrounddisplacementvalueswereacquiredbythismethod.Theclassificationofthese grounddisplacementvaluesbyslidingvelocityanddisplacementclearlyillustratedthegroundactivitydeformationstatesofvariouszones inthestudiedareawereclearlyobtained, determinedsomepotentialmovinglandslidesandactivelandslidesandidentifiedthedanger￾ouszonesoflandslides.ForanactivelandslidenumberedNo.L1R-6, thedeformationdetectedbyD-InSARanalysisshowsanaccord￾ancetendencywiththeonebyglobalpositioningsystem(GPS)monitoring.Intheend, single-pointerrorswhichexistedintheresult ofD-InSARtechniquewereanalyzedandagridfunctionerrorinterpolationmethodofD-InSARandGPSwasproposedtoimprovethe monitoringaccuracyofD-InSARtechnique. KEYWORDS syntheticapertureradar(SAR);landslides;remotesensing;monitoring;globalpositioningsystem(GPS) 收稿日期:2010--04--08 基金项目:国家自然科学基金资助项目(No.4972229) 滑坡是一种常见的自然地质灾害, 给国家建设 和人民生命财产造成巨大损失 .对滑坡的监测 、预 报和早期防治一直是世界各国广泛关注和研究的问 题 .传统上对滑坡的监测和预警主要采用钻孔倾斜 仪 、位移计和自动伸缩计等现场仪器观测 ,以及利用 全球定位系统 (GPS)、时域 反射法 (timedomain reflectometry, TDR)技术和分布式光纤等电子远程 监测技术等.这些技术的优点是单点监测精度高, 但无法满足人员难以进入的区域以及广域大面积区 域的地面活动监测需要 .随着合成孔径雷达干涉 DOI :10 .13374 /j .issn1001 -053x .2011 .02 .021

。132 北京科技大学学报 第33卷 ynthetic aperture radar interferon etry InSAR) 相位项，形变后的单视复数数据与主影像形成的干 的发展，合成孔径雷达差分干涉(differen tial 涉SR图像既包含地形相位项又包含形变相位项， synthetic aperture radar n erfermetry D-InSAR) 利用己知的卫星轨道数据和地形数据（如己知的 术正好能够解决连续大面积上非常小的地面活动监 DEM能够很容易去除地形相位项，差分后仅剩下 测，而且具有高精度、高分辨率、全天候、监测成本低 与地表形变有关的相位项，然后将其转变为微小形 以及能够监测人员无法进入的区域等优点，是非常 变值.其测量精度为其所用电磁波波长的量级，因 具有潜力的地面活动辨识和监测技术，成为滑坡、沉 此在微波频率范围内可达到厘米级的精度匀，与永 降、地震和火山等地质灾害活动调查监测实践应用 久散射体PS技术及GPS技术等结合，精度可达到 上的研究热点.20世纪90年代以来，在滑坡监测领 毫米级0 域的研究和应用中已有许多报道1.1995年 D-ISAR技术根据其消除地形相位的方法不 Achache等山首先将D-nSAR技术应用在大规模小 同，主要分为DEM的两轨法、三轨法和四轨法它们 位移的滑坡监测上，获得了与地面监测相同精度的 的原理是基本相同的.本文采用三轨法，其基本原 结果.1996年Frnea等获得了La Clapierey滑坡 理如图1所示山. 位移活动的范围，揭示了整个滑坡的滑动趋势及地 面不连续监测不能监测到的滑坡上小规模的不稳定 QBwB,只4 性.1999年Vemeer等3也对La Valet附近的滑 B. 坡进行了位移的提取，获得了与地面监测相一致的 1.2md的位移精度.2000年Rz0等揭示了 Randazzo滑坡整个滑坡体的形变位移状态.2003年 Paolo等l对Maratea Valley内的滑坡进行了研究， 获得了D-ISAR GPS汲电子测距仪DM三种监测 基本相一致的结果，并在此基础上提出了调查期间 内滑坡的不稳定模型.2003年S9 uarzon等I9也对 La Vaetter滑坡进行了研究，从15幅差分干涉图中 图1三轨法差分干涉测量成像几何示意图 推演出滑坡的位移值。得到了与地面监测相一致的 Fg Geomn etrical skech of theeorbitD hSAR 结果.2006年Cojesant垮I对Liech enstein Ap滑 坡进行了研究，得出了D-S技术最具吸引力和 A、A和A是卫星三次对同一地区成像的位 最可靠的优势就是对广域区域上稳定与不稳定区域 置，下和为卫星三次成像位置到地面监测点P 的定性辨识，以及大地质灾害的危险区定性划分和 的距离，B和B分别为卫星在A、A及A、A位置 对缓慢滑坡上不同滑动速率区域的辨识.2006年与 成像时的空间基线，和？分别为B和B与水平 2008年，Singh roy等8-列利用D-nS4R技术分别对 向倾角，0为传感器电磁波的入射角，B、B1和 老滑坡和冻融性滑坡活动进行了研究，都取得了较 B从、B分别为基线B和B沿视线向的平行和垂 好的结果.但是，目前我国利用DSAR技术在滑 直分量.由图中几何关系及B《十可得 坡方面研究的应用实例还很少. 子=+居-2 r B cos-0+A 2 1D-InSAR技术基本原理和处理流程 I一Bs(0-月)会BM (1) 1.1D-nSAR技术的基本原理 式(1表明用干涉测量法所得到的相位差与视线方 合成孔径雷达差分干涉(D-SR)技术，是通 向的基线分量成正比.这里设定在A处获得第一 过对SAR数据的处理将同一场景与地形和形变有 幅影像且为主影像，假设在地表未发生形变前在A 关的两幅干涉SR图像的相位差转变为地表微小 处获取第二幅影像，所以第二幅影像与A处的主影 形变值的技术.两幅干涉SAR图像是由形变前的 像形成的干涉SR图像，其干涉相位差仅包含地形 一帧单视复数(SIC)SAR数据和形变后的一帧单视 信息，即两次对地面同一点P成像的干涉相位差 复数(SLC)SAR数据，分别与选定的SC主影像数 △中12可表示为 据进行精确的配准和干涉产生的.形变前的单视复 数数据与主影像形成的干涉SAR图像仅包含地形 △p:-,-4:-交p(-D=

北 京 科 技 大 学 学 报 第 33卷 (syntheticapertureradarinterferometry, InSAR)技术 的 发 展, 合 成 孔 径 雷 达 差 分 干 涉 (differential syntheticapertureradarinterferometry, D-InSAR)技 术正好能够解决连续大面积上非常小的地面活动监 测 ,而且具有高精度、高分辨率 、全天候、监测成本低 以及能够监测人员无法进入的区域等优点 , 是非常 具有潜力的地面活动辨识和监测技术 ,成为滑坡、沉 降 、地震和火山等地质灾害活动调查监测实践应用 上的研究热点.20世纪 90年代以来 ,在滑坡监测领 域的研究和应 用中已有许多报道 [ 1--9] .1995 年 Achache等 [ 1] 首先将 D-InSAR技术应用在大规模小 位移的滑坡监测上, 获得了与地面监测相同精度的 结果.1996年 Fruneau等 [ 2]获得了 LaClapiere滑坡 位移活动的范围 ,揭示了整个滑坡的滑动趋势及地 面不连续监测不能监测到的滑坡上小规模的不稳定 性 .1999年 Vietmeier等 [ 3] 也对 LaValette附近的滑 坡进行了位移的提取, 获得了与地面监测相一致的 1.2cm·d -1的位移精度 .2000 年 Rizzo等 [ 4] 揭示了 Randazzo滑坡整个滑坡体的形变位移状态.2003年 Paolo等 [ 5]对 MarateaValley内的滑坡进行了研究 , 获得了 D-InSAR、GPS及电子测距仪 EDM三种监测 基本相一致的结果, 并在此基础上提出了调查期间 内滑坡的不稳定模型 .2003年 Squarzoni等 [ 6] 也对 LaValette滑坡进行了研究, 从 15幅差分干涉图中 推演出滑坡的位移值, 得到了与地面监测相一致的 结果.2006年 Colesanti等 [ 7] 对 LiechtensteinAlps滑 坡进行了研究, 得出了 D-InSAR技术最具吸引力和 最可靠的优势就是对广域区域上稳定与不稳定区域 的定性辨识 ,以及大地质灾害的危险区定性划分和 对缓慢滑坡上不同滑动速率区域的辨识.2006年与 2008年 , Singhroy等 [ 8--9] 利用 D-InSAR技术分别对 老滑坡和冻融性滑坡活动进行了研究, 都取得了较 好的结果.但是, 目前我国利用 D-InSAR技术在滑 坡方面研究的应用实例还很少 . 1 D-InSAR技术基本原理和处理流程 1.1 D-InSAR技术的基本原理 合成孔径雷达差分干涉 (D-InSAR)技术, 是通 过对 SAR数据的处理将同一场景与地形和形变有 关的两幅干涉 SAR图像的相位差转变为地表微小 形变值的技术.两幅干涉 SAR图像是由形变前的 一帧单视复数(SLC)SAR数据和形变后的一帧单视 复数(SLC)SAR数据 ,分别与选定的 SLC主影像数 据进行精确的配准和干涉产生的.形变前的单视复 数数据与主影像形成的干涉 SAR图像仅包含地形 相位项,形变后的单视复数数据与主影像形成的干 涉 SAR图像既包含地形相位项又包含形变相位项, 利用已知的卫星轨道数据和地形数据 (如已知的 DEM)能够很容易去除地形相位项 , 差分后仅剩下 与地表形变有关的相位项, 然后将其转变为微小形 变值 .其测量精度为其所用电磁波波长的量级, 因 此在微波频率范围内可达到厘米级的精度 [ 5] , 与永 久散射体 PS技术及 GPS技术等结合 , 精度可达到 毫米级 [ 10] . D-InSAR技术根据其消除地形相位的方法不 同, 主要分为 DEM的两轨法 、三轨法和四轨法,它们 的原理是基本相同的 .本文采用三轨法, 其基本原 理如图 1所示 [ 11] . 图 1 三轨法差分干涉测量成像几何示意图 Fig.1 Geometricalsketchofthree-orbitD-InSAR A1 、A2 和 A3 是卫星三次对同一地区成像的位 置, r1 、r2 和 r3 为卫星三次成像位置到地面监测点 P 的距离, B1 和 B2 分别为卫星在 A1 、A2 及 A1 、A3 位置 成像时的空间基线 , 1 和 2 分别为 B1 和 B2 与水平 向倾角 , θ为传感器电磁波的入射角, B1∥ 、 B1⊥ 和 B2∥ 、B2⊥分别为基线 B1 和 B2 沿视线向的平行和垂 直分量.由图中几何关系及 B1 r1 +r2 可得 r 2 2 =r 2 1 +B 2 1 -2r1 B1cos π 2 -θ+ 1 r1 -r2 B1sin(θ- 1) B1∥ (1) 式(1)表明用干涉测量法所得到的相位差与视线方 向的基线分量成正比 .这里设定在 A1 处获得第一 幅影像且为主影像 ,假设在地表未发生形变前在 A2 处获取第二幅影像 ,所以第二幅影像与 A1 处的主影 像形成的干涉 SAR图像 ,其干涉相位差仅包含地形 信息 ,即两次对地面同一点 P成像的干涉相位差 Δ 12可表示为 Δ 12 = 1 - 2 = 2π λ ρ(r2 -r1)= · 132·

第2期 王桂杰等：差分干涉合成孔径雷达技术在广域滑坡动态辨识上的实验研究 ·133 B0-网)=资5 △中12B/ 入 (2) △13一马0 (4) 式中，中，和中2分别为卫星在A和A处对地面点P 假设△R为视线向形变量，则由式(2)~(4推 成像的相位值，入为微波波长，星载重复轨道P=2 得由视线向形变量引起的干涉条纹图相位差△中： 假设发生形变后在A处获取了第三幅影像，所以第 可表示为 三幅影像与主影像形成的干涉SAR图像的干涉相 B △中≈△中3-△中12=△中2 B/ 红△R 位差，既包含地形信息又包含地表形变信息.且由 入 于获得的影像间要求基线足够小，所以可近似看作 (5) 0不变，即此时A和A两次对地面同一点P成像 式(5左边的各量可由干涉条纹图的相位和轨道参 的干涉相位差位差△中：可表示为： 数计算得到，进而可确定影像每点的视线向形变量 △e=,-,-交p(-)= △R分解后得到水平形变量和垂直形变量. 1.2D-nSAR技术的数据处理流程 Bm(0-)=年B 三轨法是指采用卫星三次成像的三景数据进行 入 (3) 分析，其数据处理流程如图2所示.具体包括以下 由式(2和(3得： 几个过程 单视复数铺影液？ 单视复数主影像1 单视复数辅影像3 相关系数阳行准 相关系数租配准 相干系数情准 相干系数精配准 父共纯相彩 复共矩相乘 83时多混处理 83的多视处厘 生或地形对干沙条纹图1 生成地形/形变对干法条纹图。 粗精度去除平地效纹 相情度W去除半地效应 更形均值半清连波 形均值平滑滤波 区成州达电位解新+ 区城增长法相仪解缩 机通参数的提取 轨道参数的提取 次力分生成力分干涉图 差分干涉图矩形均值平滑滤波 区域线墙长法麦分干涉图相位解推 相 地 理 改小形变图 图2DSAR技术三轨法数据处理流程 Fg 2 Fl chart ofD-InSAR da n processing (1)选择合适的雷达卫星数据，并转换为干涉 (2)选择合适的SC地形数据对和地形形变 处理的单视复数(S○数据.处理过程中首先必须 数据对.为了保证数据的分析质量需要首先对单视 将获得的卫星原始数据经聚焦处理为SC数据，作 复数对进行基线的预估计，根据预估计的时间基线、 为后期处理的输入数据，本研究获取的是己处理好 空间基线和多普勒频移差来选取数据对及确定主辅 的单视复数AIOS卫星数据. 图像.对三轨差分干涉主要考虑时间基线和空间基

第 2期 王桂杰等:差分干涉合成孔径雷达技术在广域滑坡动态辨识上的实验研究 - 4π λ B1 sin(θ- 1 )=- 4π λ B1∥ (2) 式中, 1 和 2 分别为卫星在 A1 和 A2 处对地面点 P 成像的相位值, λ为微波波长, 星载重复轨道 ρ=2. 假设发生形变后在 A3 处获取了第三幅影像 ,所以第 三幅影像与主影像形成的干涉 SAR图像的干涉相 位差, 既包含地形信息又包含地表形变信息.且由 于获得的影像间要求基线足够小, 所以可近似看作 θ不变 , 即此时 A1 和 A3 两次对地面同一点 P成像 的干涉相位差位差 Δ 13可表示为: Δ 13 = 1 - 3 = 2π λ ρ(r3 -r1 )= - 4π λ B2 sin(θ- 2)=- 4π λ B2∥ (3) 由式(2)和(3)得: Δ 12 Δ 13 = B1∥ B2∥ (4) 假设 ΔRd为视线向形变量 ,则由式(2)～ (4)推 得由视线向形变量引起的干涉条纹图相位差 Δ d 可表示为 Δ d≈Δ 13 -Δ 12 =Δ 12 1 - B2∥ B1∥ =- 4π λ ΔRd (5) 式(5)左边的各量可由干涉条纹图的相位和轨道参 数计算得到,进而可确定影像每点的视线向形变量 ΔRd,分解后得到水平形变量和垂直形变量. 1.2 D-InSAR技术的数据处理流程 三轨法是指采用卫星三次成像的三景数据进行 分析 ,其数据处理流程如图 2所示 .具体包括以下 几个过程 . 图 2 D-InSAR技术三轨法数据处理流程 Fig.2 FlowchartofD-InSARdataprocessing (1)选择合适的雷达卫星数据, 并转换为干涉 处理的单视复数 (SLC)数据.处理过程中首先必须 将获得的卫星原始数据经聚焦处理为 SLC数据 ,作 为后期处理的输入数据 , 本研究获取的是已处理好 的单视复数 ALOS卫星数据 . (2)选择合适的 SLC地形数据对和地形 /形变 数据对.为了保证数据的分析质量需要首先对单视 复数对进行基线的预估计,根据预估计的时间基线、 空间基线和多普勒频移差来选取数据对及确定主辅 图像 .对三轨差分干涉主要考虑时间基线和空间基 · 133·

。134 北京科技大学学报 第33卷 线，原则上时间基线最好小于6个月，PALSAR传感 条纹现象.平地效应的存在，使干涉条纹过于密集， 器空间基线要小于1200m越小相关性越好，所得 加大了相位解缠的难度，因此必须采用精确的轨道 结果也更精确些).对于地形数据对原则上要求有 数据和已知的辅助DM予以去除.即去除后的地 较短的时间基线，而地形形变数据对要求具有较短 形干涉条纹图仅含有地形相位项，地形形变干涉条 的空间基线. 纹图只含有地形相位项和形变相位项.另外在干涉 (3)对选好的地形和地形形变数据对进行 图形成过程中会有许多原因产生干涉相位噪声，所 SC主辅图像的精确配准.SAR的SC数据是由离 以在平地效应去除后还要进行滤波处理 散的像素构成，这些像素是对目标地面特性的离散 (6)相位解缠.上面干涉条纹图的生成己经提 点取样.D-S技术要依据地面上同一目标点在 到，干涉图上的相位为[一π，π)主值范围内的相位 两次成像中的相位差，因此必须要掌握地面每一目 差值，要把它表现为实际真实的相位值必须进行相 标点在两次成像中的对应位置.只有使一次成像中 位解缠展开为每点的真实相位. 各像素对应在地面上的离散样点与另一次成像的各 (7)二次差分产生差分干涉图.将解缠后的地 像素对应在地面上的离散样点正好完全重合，才可 形形变干涉条纹图与地形干涉条纹图进行差分，消 以用两次成像中位置相同的各像素构成的像素对来 去地形相位项，得到只含地表形变项的差分干涉图， 求取地面上每一样点在两次成像间的相位差.因 并对差分后的差分干涉图也进行滤波、解缠等处理. 此，在干涉图生成之前必须进行主辅图像的像素 (8)提取轨道参数、相位信息转换为形变信息、 配准。 地理编码.最后利用提取的轨道参数将相位信息转 (4)地形和地形形变数据对干涉条纹图的产 化为形变信息并进行地理编码，得到研究区域地表 生.在完成数据配准之后，将主、辅图像对应像元复 的微小形变图. 数值进行共轭相乘： ua=4y=|y||y|e2= 2滑坡动态辨识实验研究 141g1e*2) (6) 2.1研究区域和数据 式中，=||，号=|y|马表示号的共轭 研究区域位于四川省金沙江流域下游乌东德水 复数. 电站库区，如图3所示.整幅数据研究区（底图为 Google Eart影像截图)中的紫色线框内的区域为获 把所得复数的相位信息单独显示出来称为 取的整景SAR数据的区域：绿色线框内为20m义 干涉相位图或干涉图.凸的相位信息是主、辅图像 每一同名点上的相位差，在干涉图上用相位主值 20的主研究区域，即主研究区域详图（底图为 [一，π)显示出来，这个相位称为缠绕相位. SPo5与12000的高精度航片融合的影像：山R- 由于原始数据在方位向上像元分辨率高于斜距 6号滑坡（底图为航测影像立体视图为主研究区域 方向的分辨率，在未进行多视处理的SC图像中其 中对水库正常运营有重要影响的正处活动状态的 方位向和距离向按其分辨率显示，则两个方向的比 滑坡. 例尺不一致，图像被拉伸变形，显示比例差别较大. 本文所用的数据为AIOS卫星PALSAR传感器 为了使显示出来的影像与地形图对照时同名地物形 获得的三帧升轨SR数据，其极化方式为水平极化 状基本一致及对干涉条纹进行显示和评估，在干涉 (HH)沮分辨率为6.259波长采用的是23a的L 条纹图生成之后需要进行多视处理.多视处理即按 波段，其数据参数如表1所示.按前面处理流程中 某一比例作压缩处理.例如按方位向与距离向之比 所叙述的原则将编号为AISP104820510(2008- 为8压缩，压缩后原始数据中方位向八行变为一 01-12)数据作为主影像.选取ALS04820510 行，距离向三列变为一列. (200801-12和APSn11530510(200802-27) (5)干涉平地效应的去除和滤波.干涉条纹图 作为地形数据对，ALPSRI04820510(2008-01-12) 上的相位由两部分组成，一是地形的相对高度及地 和AIPSRP18240510(2008-04-13)作为地形形 表形变引起的相位成分，另一部分是由平地效应引 变数据对.将30m分辨率的ASIER DEM与等高线 起的.平地效应是指平地相位在干涉条纹中所表现 生成的2.5m分辨率DM镶嵌后的数字高程数据 出来的随距离和方位的变化而呈周期性变化的密集 作为数据处理过程中的辅助DM

北 京 科 技 大 学 学 报 第 33卷 线 ,原则上时间基线最好小于 6个月, PALSAR传感 器空间基线要小于 1 200 m(越小相关性越好 ,所得 结果也更精确些).对于地形数据对原则上要求有 较短的时间基线 ,而地形 /形变数据对要求具有较短 的空间基线 . (3)对选好的地形和地形 /形变数据对进行 SLC主辅图像的精确配准 .SAR的 SLC数据是由离 散的像素构成,这些像素是对目标地面特性的离散 点取样 .D-InSAR技术要依据地面上同一目标点在 两次成像中的相位差, 因此必须要掌握地面每一目 标点在两次成像中的对应位置 .只有使一次成像中 各像素对应在地面上的离散样点与另一次成像的各 像素对应在地面上的离散样点正好完全重合 , 才可 以用两次成像中位置相同的各像素构成的像素对来 求取地面上每一样点在两次成像间的相位差 .因 此 ,在干涉图生成之前必须进行主辅图像的像素 配准. (4)地形和地形 /形变数据对干涉条纹图的产 生 .在完成数据配准之后 ,将主、辅图像对应像元复 数值进行共轭相乘: uint=u1u ＊ 2 = u1 e j 1 u2 e -j 2 = u1 u2 e j( 1 - 2 ) (6) 式中, u1 = u1 e j 1 , u2 = u2 e j 2 , u ＊ 2 表示 u2 的共轭 复数. 把所得复数 uint的相位信息单独显示出来称为 干涉相位图或干涉图 .uint的相位信息是主 、辅图像 每一同名点上的相位差 , 在干涉图上用相位主值 [ -π, π)显示出来 ,这个相位称为缠绕相位 . 由于原始数据在方位向上像元分辨率高于斜距 方向的分辨率,在未进行多视处理的 SLC图像中其 方位向和距离向按其分辨率显示, 则两个方向的比 例尺不一致 ,图像被拉伸变形, 显示比例差别较大 . 为了使显示出来的影像与地形图对照时同名地物形 状基本一致及对干涉条纹进行显示和评估 , 在干涉 条纹图生成之后需要进行多视处理.多视处理即按 某一比例作压缩处理 .例如按方位向与距离向之比 为 8∶3压缩, 压缩后原始数据中方位向八行变为一 行 ,距离向三列变为一列. (5)干涉平地效应的去除和滤波 .干涉条纹图 上的相位由两部分组成 , 一是地形的相对高度及地 表形变引起的相位成分 , 另一部分是由平地效应引 起的.平地效应是指平地相位在干涉条纹中所表现 出来的随距离和方位的变化而呈周期性变化的密集 条纹现象 .平地效应的存在 ,使干涉条纹过于密集, 加大了相位解缠的难度, 因此必须采用精确的轨道 数据和已知的辅助 DEM予以去除 .即去除后的地 形干涉条纹图仅含有地形相位项 ,地形 /形变干涉条 纹图只含有地形相位项和形变相位项.另外在干涉 图形成过程中会有许多原因产生干涉相位噪声, 所 以在平地效应去除后还要进行滤波处理 . (6)相位解缠 .上面干涉条纹图的生成已经提 到, 干涉图上的相位为 [ -π, π)主值范围内的相位 差值 ,要把它表现为实际真实的相位值必须进行相 位解缠,展开为每点的真实相位. (7)二次差分产生差分干涉图.将解缠后的地 形 /形变干涉条纹图与地形干涉条纹图进行差分 ,消 去地形相位项 ,得到只含地表形变项的差分干涉图, 并对差分后的差分干涉图也进行滤波、解缠等处理. (8)提取轨道参数 、相位信息转换为形变信息、 地理编码 .最后利用提取的轨道参数将相位信息转 化为形变信息并进行地理编码, 得到研究区域地表 的微小形变图 . 2 滑坡动态辨识实验研究 2.1 研究区域和数据 研究区域位于四川省金沙江流域下游乌东德水 电站库区 , 如图 3所示 .整幅数据研究区 (底图为 GoogleEarth影像截图 )中的紫色线框内的区域为获 取的整景 SAR数据的区域;绿色线框内为 20 km× 20 km的主研究区域 , 即主研究区域详图 (底图为 Sport5与 1∶2 000的高精度航片融合的影像 );L1R-- 6号滑坡(底图为航测影像立体视图 )为主研究区域 中对水库正常运营有重要影响的正处活动状态的 滑坡 . 本文所用的数据为 ALOS卫星 PALSAR传感器 获得的三帧升轨 SAR数据, 其极化方式为水平极化 (HH)且分辨率为 6.25 m,波长采用的是 23 cm的 L 波段 ,其数据参数如表 1所示 .按前面处理流程中 所叙述的原则将编号为 ALPSRP104820510(2008-- 01--12)数据作为主影像 .选取 ALPSRP104820510 (2008--01--12)和 ALPSRP111530510 (2008--02--27) 作为地形数据对 , ALPSRP104820510(2008--01--12) 和 ALPSRP118240510 (2008--04--13)作为地形 /形 变数据对 .将 30m分辨率的 ASTERDEM与等高线 生成的 2.5 m分辨率 DEM镶嵌后的数字高程数据 作为数据处理过程中的辅助 DEM. · 134·

第2期 王桂杰等：差分干涉合成孔径雷达技术在广域滑坡动态辨识上的实验研究 135 整望数据研究区城 丰研究区域详露 1.R-6滑坡 图3滑坡研究区域 Fg 3 Sudied area of landslides 表1所用ALOS PALSAR卫星的数据参数 Table 1 Parame ters ofAIOS PALSAR 像对编号 卫星传感器 景编号 获取日期 时间间隔/d垂直基线m多普勒频移差Hz2π模糊度高度m ALOS PALSAR A1P9Rh04820510 2008-01-12 Ar 46 438883 -14625 131.471 ALOS PALSAR AIPh11530510 2008-02-27 ALOS PALSAR AIP9Rh04820510 2008-01-12 92 813466 -33191 70.981 ALOS PALSAR A1Ph182405102008-04-13 2.2微小形变位移图的产生 过处理后的干涉条纹图则仅含地形和形变相位信 选定的地形数据对A及地形形变数据对A 息，如图5所示.从图5可看出去除平地效应和滤 经配准、复数共轭（干涉和多视处理后，得到如图4 波后的干涉图，不仅将均匀密集的小条纹去除而且 所示的初始干涉条纹图.由图可明显看到由于时间 干涉条纹较清晰. 引起的失相关及垂直基线距等原因，46的地形对 干涉条纹图A明显好于92d的形变对干涉条纹 图A 滤波后的整幅干涉条纹图 左上角小红线框内的放大图 L整幅初始干涉条纹图 左上角小红线相框内的放大图 1,整平和滤波后的整幅干浅条纹图 左上角小红线框内的放大 1整制初始干涉条纹图 左上角小红线框内的放大图 图5整平和滤波后的干涉图 Fg 5 hterfe ogem a fter flatten ing and filtering 图4初步形成的干涉图 F4 Inital interfeogram 另外，在干涉处理过程中主辅图像的相干性好 为了去除干涉图中由平地效应而引起的较密集 坏，是影响数据质量和结果精确与否的重要因素. 均匀的小干涉条纹以及干涉中产生的各种噪声影 相干系数是衡量相干性好坏的标准.图6为干涉数 响，对产生的干涉条纹图进行整平和滤波处理.经 据对A和A的相干系数图，其为[0]间的值.相

第 2期 王桂杰等:差分干涉合成孔径雷达技术在广域滑坡动态辨识上的实验研究 图 3 滑坡研究区域 Fig.3 Studiedareaoflandslides 表 1 所用 ALOSPALSAR卫星的数据参数 Table1 ParametersofALOSPALSAR 像对编号 卫星传感器 景编号 获取日期 时间间隔 /d 垂直基线 /m 多普勒频移差 /Hz 2π模糊度高度 /m AT ALOSPALSAR ALPSRP104820510 2008--01--12 46 438.883 -14.625 131.471 ALOSPALSAR ALPSRP111530510 2008--02--27 AD ALOSPALSAR ALPSRP104820510 2008--01--12 92 813.466 -33.191 70.931 ALOSPALSAR ALPSRP118240510 2008--04--13 2.2 微小形变位移图的产生 选定的地形数据对 AT及地形 /形变数据对 AD 经配准 、复数共轭 (干涉)和多视处理后,得到如图 4 所示的初始干涉条纹图.由图可明显看到由于时间 引起的失相关及垂直基线距等原因, 46 d的地形对 干涉条纹图 AT明显好于 92 d的形变对干涉条纹 图 AD. 图 4 初步形成的干涉图 Fig.4 Initialinterferogram 为了去除干涉图中由平地效应而引起的较密集 均匀的小干涉条纹以及干涉中产生的各种噪声影 响 ,对产生的干涉条纹图进行整平和滤波处理.经 过处理后的干涉条纹图则仅含地形和形变相位信 息, 如图 5所示 .从图 5可看出去除平地效应和滤 波后的干涉图 ,不仅将均匀密集的小条纹去除而且 干涉条纹较清晰. 图 5 整平和滤波后的干涉图 Fig.5 Interferogramafterflatteningandfiltering 另外 , 在干涉处理过程中主辅图像的相干性好 坏, 是影响数据质量和结果精确与否的重要因素. 相干系数是衡量相干性好坏的标准.图 6为干涉数 据对 AT和 AD的相干系数图 ,其为[ 0, 1]间的值 .相 · 135·

。136 北京科技大学学报 第33卷 将去除平地效应、滤波后的干涉相位图进行相 位解缠、二次差分和轨道精化等处理.最后，将相位 信息转换为形变信息并地理编码得到沿卫星视线向 的微小形变图.将微小形变位移图转换为沿垂直方 向和水平方向的位移图，并根据滑动位移的移动速 率及位移的大小进行分级处理即得到各个区域的滑 L数据对相十系数图 A数据对相干系数图 动位移情况.如图7所示，在整个研究区域形变图 中的四边形区域为主研究区域不规则的多边形区 图6相干系数图 域为2.5分辨率DM镶嵌的区域，图中红色区域 Fg 6 Map ofccherence coefficient 为抬升区域，绿色区域为沉降区域.因为滑坡的滑 干系数值大的点，说明干涉数据对在此点的相干性 移主要伴随的是地表沉降，因此这里主要研究沉降 好即形变结果较精确，图中用灰度图来显示越亮的 区域. 地方相干系数越大 ■-0.403-0200 。020一0100 000050 ■-0403-020 02-0 0-0050 00050 ■0.050-0.100 ■00s0-0100 ■0100-0.200 ■0.100-0.200 ■02000306 ■0.200-0.25 整景研究区域上的垂直滑动状态图 主研究区域上的垂直滑动状态图 ■-1531-0774 ■-0774-0384 厘0384-0194 ■-0.0. 0194-0 ■0-0.194 -0.38 ■0.194-0387 ■0387-，774 -0 m07741.538 整景研究区域上的水平滑动状态图 b 主研究区域上的水平滑动状态图 图7研究区域滑动位移状态分级图（单位：m)()研究区域92的垂直滑动状态图：(b)研究区域92d的水平滑动状态图 Fg7 Depma tonmaps of the suded a宠a(unit m).(两g2 d vertical depmationm现(b)2 d horiz知alde fma t知m即 在主研究区域垂直滑动位移形变状态图中，将为滑动速率为0~0.0167m月-1（位移值为0~ 抬升区域由位移值大小按红色由浅到深划分为0~ 0.05m).滑动速率为0.0167~0.0334m月-（位 0.2m和大于0.2m的两个区域.将垂直沉降滑动 移值为005~01m)、滑动速率为0.0334~0.0667 位移值按移动速率及位移大小按绿色由浅到深划分 m月-（位移值为01~02m)以及滑动速率大于

北 京 科 技 大 学 学 报 第 33卷 图 6 相干系数图 Fig.6 Mapofcoherencecoefficient 干系数值大的点 ,说明干涉数据对在此点的相干性 好即形变结果较精确, 图中用灰度图来显示越亮的 地方相干系数越大. 将去除平地效应 、滤波后的干涉相位图进行相 位解缠、二次差分和轨道精化等处理 .最后, 将相位 信息转换为形变信息并地理编码得到沿卫星视线向 的微小形变图 .将微小形变位移图转换为沿垂直方 向和水平方向的位移图, 并根据滑动位移的移动速 率及位移的大小进行分级处理即得到各个区域的滑 动位移情况.如图 7所示, 在整个研究区域形变图 中的四边形区域为主研究区域, 不规则的多边形区 域为 2.5 m分辨率 DEM镶嵌的区域, 图中红色区域 为抬升区域, 绿色区域为沉降区域 .因为滑坡的滑 移主要伴随的是地表沉降, 因此这里主要研究沉降 区域 . 图 7 研究区域滑动位移状态分级图(单位:m).(a)研究区域 92d的垂直滑动状态图;(b)研究区域 92d的水平滑动状态图 Fig.7 Deformationmapsofthestudiedarea(unit:m):(a)92dverticaldeformationmap;(b)92dhorizontaldeformationmap 在主研究区域垂直滑动位移形变状态图中 ,将 抬升区域由位移值大小按红色由浅到深划分为 0 ～ 0.2m和大于 0.2 m的两个区域 .将垂直沉降滑动 位移值按移动速率及位移大小按绿色由浅到深划分 为滑动速率为 0 ～ 0.016 7 m·月 -1 (位移值为 0 ～ 0.05 m)、滑动速率为 0.016 7 ～ 0.033 4 m·月 -1 (位 移值为 0.05 ～ 0.1 m)、滑动速率为 0.033 4 ～ 0.0667 m·月 -1(位移值为 0.1 ～ 0.2 m)以及滑动速率大于 · 136·

第2期 王桂杰等：差分干涉合成孔径雷达技术在广域滑坡动态辨识上的实验研究 ·137° 0.0667m月-1（位移值大于0.2m)的四个区域来 同时进行GPS监测的山R6号滑坡为例，将其 显示.在水平滑动位移形变状态图中，按水平位移 D-ISA监测结果进行详细分析.图8为山R-6号 与垂直位移所成比例的倍数，将抬升区域由水平位 滑坡上D-SAR监测的详细滑动形变图.在垂直滑 移绝对值大小按红色由浅到深划分为0~0.774m 动形变位移图中按滑动速率小于0.0167m月-1 和大于0.774两个区域.将沉降区的水平位移按 (位移值小于0.05m,滑动速率为0.0167~ 移动速率绝对值及位移绝对值大小也进行相同的划 0.0334m月-1（位移值为0.05~01m)和滑动速 分，按绿色由浅到深划分为滑动速率为0~0.0647 率大于0.0334m月1（位移值大于0.1m从慢到 m月-1（位移值为0~0194m)、滑动速率为 快分为三个区域，在水平形变位移图中按绝对值滑 0.0647-0.129m月（位移值为0.194~0.387 动位移速率小于0.0647m月-1（位移值小于0194 m)滑动速率为0.129-0.258m月-（位移值为 m)、滑动速率为0.0647~0.129m月-1（位移值为 0.387-0774m以及滑动速率大于0.258m月-1 0.194~0.387m和滑动速率大于0.129m月（位 (位移值为大于0.774m的四个区域.其中水平位 移值大于0.387m)也从慢到快分为三个区域分别 移是指南北向的水平位移，图中向南符号为负，向北 用黄、蓝和红三种颜色来显示.从形变图中能清晰 符号为正.形变图中沉降区域的垂直和水平位移分 地分辨出不同区域的滑动速率和位移情况.图中将 级图之所以按着一定的所成比例来划分，是因为得 山R-6号滑坡划分为四个区域2-，其中I区的部 到的研究区域的滑动位移是按着沿卫星视线向的形 分区域处在垂直向滑动速率小于0.0167m月- 变，根据卫星侧视倾角和轨道参数分解而得到的，所 水平向小于00647m月一)的滑动区间上，Ⅱ区 以相同区域的水平位移值与垂直位移值具有基本相 的部分区间和区处在垂直向滑动速率大于 同的比例系数，因而图中水平形变分级图和垂直形 0.0334m月-1（冰平向大于0.129m月-）的滑 变分级图的变化区域是基本相一致的，水平形变分 动区间上，其余区域都处在垂直向滑动速率为 级图仅反映了其水平滑动的速率和位移大小趋势. 0.0164-00334m·月-1（水平向0.0647~ 从研究区域的滑动形变位移图可知，其沉降区 0.129m月-1)的滑动区间上.图中红色线体为长 域最大滑动位移值为25.3四由于所用A1OS卫星 江三峡院现场监测划定的重点灾害区，可知处在 PALSAR传感器的波长为I波段23四其形变监测 D-InSR监测结果滑动位移值为5~13.3m的区 灵敏度为12波长即为12.5m(沿卫星视线向 间上. 的，再考虑到各种可能误差（后面分析，其区域内 a 观测时段92内的变动位移很小.可看到其红色抬 升区域的抬升位移值绝大多数都在小于20m的范 00.018-0050 ■0.050-0100 围内，因此抬升区域的抬升位移值大部分都处在误 m0.100-0.133 差影响范围内.因此对研究区域垂直沉降滑动位移 速率大于667m月-(0.0667m月-1位移值大 于20四水平滑动位移速率25.8m月-1(0.258 m月1)位移值大于77.4的区域划定为可能发 b 生滑坡和滑坡活动的风险区域，即为图中用紫色线 ■-0.5140387 圈圈出的两个区域.经与研究区域地质资料及光学 ■-0387-0194 ■-0.194-0.070 影像确认，此两区域一个区域位于金沙江右岸一个 较大泥石流沟的沟口附近，另一区域位于该泥石流 沟沟中汇水区的中部.其形变位移值较大的原因除 可能因为崩滑原因外，也可能是由于沟中物质迁移 图8L1R-6号滑坡的滑动形变图（单位：m,()92d的垂直 所引起的. 滑动形变图：(b)92的水平滑动形变图 F琴8 Depmation map of山R-6and时de(ni诉m:(a)92d 3与GPS对比研究 venical defmation map (b)92 d horgontal depmation map 3.1与GP结果的对比 将DSR获得的山R-6号滑坡的形变位移 为证明DSAR技术在滑坡辨识监测上的可 值与己获得的地面GPS监测结果相比较，图8中绿 行性及正确性，本文以研究区域内正处活动状态并 色字体编号的紫色点为现场GPS监测点，其中点名

第 2期 王桂杰等:差分干涉合成孔径雷达技术在广域滑坡动态辨识上的实验研究 0.066 7 m·月 -1 (位移值大于 0.2 m)的四个区域来 显示.在水平滑动位移形变状态图中 , 按水平位移 与垂直位移所成比例的倍数 ,将抬升区域由水平位 移绝对值大小按红色由浅到深划分为 0 ～ 0.774 m 和大于 0.774 m两个区域.将沉降区的水平位移按 移动速率绝对值及位移绝对值大小也进行相同的划 分 ,按绿色由浅到深划分为滑动速率为 0 ～ 0.064 7 m·月 -1 (位移值为 0 ～ 0.194 m)、滑动 速率为 0.064 7 ～ 0.129 m·月 -1 (位移值为 0.194 ～ 0.387 m)、滑动速率为 0.129 ～ 0.258 m·月 -1 (位移值为 0.387 ～ 0.774 m)以及滑动速率大于 0.258 m·月 -1 (位移值为大于 0.774 m)的四个区域 .其中水平位 移是指南北向的水平位移 ,图中向南符号为负 ,向北 符号为正.形变图中沉降区域的垂直和水平位移分 级图之所以按着一定的所成比例来划分, 是因为得 到的研究区域的滑动位移是按着沿卫星视线向的形 变 ,根据卫星侧视倾角和轨道参数分解而得到的,所 以相同区域的水平位移值与垂直位移值具有基本相 同的比例系数,因而图中水平形变分级图和垂直形 变分级图的变化区域是基本相一致的, 水平形变分 级图仅反映了其水平滑动的速率和位移大小趋势. 从研究区域的滑动形变位移图可知, 其沉降区 域最大滑动位移值为 25.3cm,由于所用 ALOS卫星 PALSAR传感器的波长为 L波段 23 cm, 其形变监测 灵敏度为 1 /2 波长即为 12.5 cm(沿卫星视线向 的 ),再考虑到各种可能误差 (后面分析),其区域内 观测时段 92 d内的变动位移很小 .可看到其红色抬 升区域的抬升位移值绝大多数都在小于 20 cm的范 围内, 因此抬升区域的抬升位移值大部分都处在误 差影响范围内.因此对研究区域垂直沉降滑动位移 速率大于 6.67 cm·月 -1 (0.066 7 m·月 -1 )位移值大 于 20 cm、水平滑动位移速率 25.8 cm·月 -1 (0.258 m·月 -1 )位移值大于 77.4 cm的区域划定为可能发 生滑坡和滑坡活动的风险区域 ,即为图中用紫色线 圈圈出的两个区域.经与研究区域地质资料及光学 影像确认,此两区域一个区域位于金沙江右岸一个 较大泥石流沟的沟口附近 , 另一区域位于该泥石流 沟沟中汇水区的中部 .其形变位移值较大的原因除 可能因为崩滑原因外, 也可能是由于沟中物质迁移 所引起的. 3 与 GPS对比研究 3.1 与 GPS结果的对比 为证明 D-InSAR技术在滑坡辨识监测上的可 行性及正确性,本文以研究区域内正处活动状态并 同时进行 GPS监测的 L1R-6 号滑坡为例 , 将其 D-InSAR监测结果进行详细分析.图 8为 L1R--6号 滑坡上 D-InSAR监测的详细滑动形变图.在垂直滑 动形变位移图中按滑动速率小于 0.016 7 m·月 -1 (位 移 值 小 于 0.05 m)、滑 动 速 率 为 0.016 7 ～ 0.033 4 m·月 -1 (位移值为 0.05 ～ 0.1 m)和滑动速 率大于 0.033 4 m·月 -1 (位移值大于 0.1 m)从慢到 快分为三个区域, 在水平形变位移图中按绝对值滑 动位移速率小于 0.064 7 m·月 -1 (位移值小于 0.194 m)、滑动速率为 0.064 7 ～ 0.129 m·月 -1 (位移值为 0.194 ～ 0.387m)和滑动速率大于 0.129m·月 -1 (位 移值大于 0.387 m)也从慢到快分为三个区域, 分别 用黄 、蓝和红三种颜色来显示 .从形变图中能清晰 地分辨出不同区域的滑动速率和位移情况 .图中将 L1R--6号滑坡划分为四个区域 [ 12--13] , 其中 Ⅰ区的部 分区域处在垂直向滑动速率小于 0.016 7 m·月 -1 (水平向小于 0.064 7 m·月 -1 )的滑动区间上, Ⅱ区 的部分区 间和 Ⅲ区 处在垂直 向滑动速率 大于 0.033 4 m·月 -1 (水平向大于 0.129 m·月 -1 )的滑 动区间上, 其余区域都处在垂直向滑动速率为 0.016 4 ～ 0.033 4 m·月 -1 (水 平 向 0.064 7 ～ 0.129 m·月 -1 )的滑动区间上 .图中红色线体为长 江三峡院现场监测划定的重点灾害区, 可知处在 D-InSAR监测结果滑动位移值为 5 ～ 13.3 cm的区 间上 . 图 8 L1R--6号滑坡的滑动形变图(单位:m).(a)92d的垂直 滑动形变图;(b)92d的水平滑动形变图 Fig.8 DeformationmapofL1R-6 landslide(unite:m):(a)92d verticaldeformationmap;(b)92dhorizontaldeformationmap 将 D-InSAR获得的 L1R--6号滑坡的形变位移 值与已获得的地面 GPS监测结果相比较, 图 8中绿 色字体编号的紫色点为现场 GPS监测点, 其中点名 · 137·

。138 北京科技大学学报 第33卷 为AI01DAI02DA103DT44和TP45等5个监 的形变监测时间间隔是2008年1月12日一4月 测点位于I区，Tb1~T5、AI02C和A103C等17 13日，GPSI、和V区的监测时间间隔是2008年 个监测点位于Ⅱ区，AIg1AAIg2AAI03AN0 1月12日一4月18日，Ⅱ区的时间间隔是2008年1 TN1、AIg2B和AIg3B等7个监测点位于I区， 月11日一4月19日.图9为两种技术监测结果的 A101 B.与AIg1C两个监测点位于V区.D-SAR 比较. 140c 120- 100 +D-INSAR监测值 一G>监测值 监测点 600 500 400 200 ◆)-INSAR监测值 监测值 100 -100 监测点 图9DR与GS滥测结果的比较.(a)见的垂直滑动形变：(b)92d的水平滑动形变 F9 Comparison ofD NS\R and GPS monitorng results (a)92d vertical defmatiop (b)92 d horizonldefmation 从图9可以看出，尽管D-SAR与GPS弹点监 探讨. 测结果的量级不同，D-SAR技术由于所用数据微 3.2误差分析 波波长的限制其量级仅为厘米且部分点由于处在相 通过分析研究可以看出D-SAR技术获得的 干性较差的区域上无法获取其形变信息，地面GS 形变值，从整体上能够辨识出各个区域的滑动大小 的监测结果为毫米精度，但其整体形变和运动趋势 和趋势，识别出可能的滑坡和危险区域.只是单点 是基本相一致的，两者都表明山R-6滑坡的I区变 数据的精确性还不够理想，由于用ALOS数据的波 动最小，Ⅱ区由于所处位置和地质结构发生形变相 长的限制监测结果处于厘米的量级，且受数据对相 对较大，而V区形变情况处于I、Ⅱ区之间.对于Ⅲ 干性的影响较大，在形变数据对相干系数较小的区 区比较结果有所不同，DSAR结果表明该区域发 域无法获取其形变信息，而在相干性较好的条件下 生的形变较大，而GPS结果显示该区域发生形变较 各点都有微小形变值等.这些误差既来源于航天 小且很多点处于抬升状态.因此在小于20m的范 SAR传感器本身也来源于DSR技术.综合分析 围内原因不排除由于误差的影响，有待于进一步 本次研究所采用的技术和数据源，其形变值误差主

北 京 科 技 大 学 学 报 第 33卷 为 AL01D、AL02D、 AL03D、TP44和 TP45等 5 个监 测点位于Ⅰ区, TP01 ～ TP15、AL02C和 AL03C等 17 个监测点位于 Ⅱ区, AL01A、AL02A、AL03A、TN10、 TN11、AL02B和 AL03B等 7 个监测点位于 Ⅲ区 , AL01B与 AL01C两个监测点位于 Ⅳ区 .D-InSAR 的形变监测时间间隔是 2008年 1月 12日 — 4 月 13日, GPSⅠ 、Ⅲ和Ⅳ区的监测时间间隔是 2008年 1月 12日 — 4月 18日, Ⅱ区的时间间隔是 2008年 1 月 11日— 4月 19日.图 9为两种技术监测结果的 比较 . 图 9 D-InSAR与 GPS监测结果的比较.(a)92d的垂直滑动形变;(b)92d的水平滑动形变 Fig.9 ComparisonofD-INSARandGPSmonitoringresults:(a)92dverticaldeformation;(b)92dhorizontaldeformation 从图 9可以看出 ,尽管 D-InSAR与 GPS单点监 测结果的量级不同, D-InSAR技术由于所用数据微 波波长的限制其量级仅为厘米且部分点由于处在相 干性较差的区域上无法获取其形变信息 ,地面 GPS 的监测结果为毫米精度 , 但其整体形变和运动趋势 是基本相一致的 ,两者都表明 L1R--6滑坡的 Ⅰ区变 动最小 , Ⅱ区由于所处位置和地质结构发生形变相 对较大 ,而Ⅳ区形变情况处于 Ⅰ 、Ⅱ区之间.对于Ⅲ 区比较结果有所不同 , D-InSAR结果表明该区域发 生的形变较大,而 GPS结果显示该区域发生形变较 小且很多点处于抬升状态 .因此在小于 20 cm的范 围内原因不排除由于误差的影响 , 有待于进一步 探讨 . 3.2 误差分析 通过分析研究可以看出 D-InSAR技术获得的 形变值,从整体上能够辨识出各个区域的滑动大小 和趋势,识别出可能的滑坡和危险区域 .只是单点 数据的精确性还不够理想 , 由于用 ALOS数据的波 长的限制监测结果处于厘米的量级 , 且受数据对相 干性的影响较大, 在形变数据对相干系数较小的区 域无法获取其形变信息, 而在相干性较好的条件下 各点都有微小形变值等 .这些误差既来源于航天 SAR传感器本身也来源于 D-InSAR技术 .综合分析 本次研究所采用的技术和数据源 , 其形变值误差主 · 138·

第2期 王桂杰等：差分干涉合成孔径雷达技术在广域滑坡动态辨识上的实验研究 ·139° 要来源以下几个方面 与GP在每一点监测结果的位移差其计算式可表 (1)SR传感器自身的限制产生的误差.受 示如下式： SR传感器所发射电磁波波长的限制，目前 Z(X)=S()一品(X￥) DSA支术本身进行监测所得处理结果的形变位 =l2,N (7) 移值受所用波长的限制：另外，受卫星重复周期的限 式中，为研究区上GPS监测点的数量，⑧( 制，对同一地区进行相同的观测至少需要十几天甚 )为研究区上第i个GP9监测点上的D-InSAR:技 至几十天的时间间隔，这使得发生较快地面活动的 术监测位移值，⑧（￥￥）为研究区上第个GPS 区域发生严重失相关因此DSAR技术目前只适 监测点上的GPS滥测位移值，Z(￥Y为研究区 合监测每月毫米至厘米到每年毫米至厘米级的滑坡 上第个GPS滥测点上的D-nSAR与GPS滥测技 活动状况. 术的位移差.利用式(7得到有GPS监测点的位置 (2)失相关产生的误差.失相关是限制 处两种监测结果单点的误差函数值. DSA技术应用的一个较为严重的问题除了由 以这些单点误差函数值为标准，进行有GPS监 于干涉数据的时间去相关外，由于基线长度、轨道轻 测点的监测区域内栅格数据反距离加权插值 微不平行、变形运动过快、植被覆盖以及在连续获取 (DW),.其插值原理公式如下： 数据期间发生的滑坡形变过多等都会导致相位的失 n Z(ǚ名)=∑ 入Z(￥号) (8) 相关问题，严重影响数据处理结果的精度，且对于数 据对上相关性较小的区域无法获得其形变信息. 式中，Z(8?)为点（吝8石）处的预测误差 (3)电磁波在大气中传播延时产生的误差.这 值，为预算计算过程中要使用的预测点周围样点 个延时相当于在干涉图中附加一个相位，这个附加 的数量，入；为预测计算过程中使用的各样点的权 相位给最终的测量结果增加了不确定性.也就是说 重，等号)为在点（芩号处获得的误差值. 只要相干性满足一定条件即使实际没有形变的点， 权重入的确定和计算公式为 也有微小位移值结果.因延时引起的位移残差一般 在所用波长量级的范围内变动. A=1会：名=1 (9) (4)数据处理过程中处理方法的异同也会对结 式中：为指数值。本研究取为2为预测点 果产生误差.像去除平地效应的整平、滤除系统噪 (爷)与各已知样点（芩之间的距离.插 声的滤波以及相位解缠等处理的好坏都会直接影响 值后得到D-SAR研究区域上有GPS滥测点的监 测量结果误差的大小， 测区域上各点的误差函数值，假设其与位置坐标函 因此，要消去DSAR技术监测的单点误差， 数关系如式(10所示： 提高监测精度，可将其与其他技术结合处理.本研 Z(xy3=〔Xy3 (10) 究根据所拥有数据提出了D-InSAR与GPS湘融合 式中，Z(xy为对应坐标点(xY处的误差函 的栅格函数误差插值消减法，对山R-6号滑坡有 数值. GPS监测点的区域D-nSAR监测结果的误差进行 在获得插值后的误差函数值后，对DSAR技 了消减 术监测结果与误差函数值Z(Xy作差，如下式 所示： 4D-InSAR与GPS相融合的栅格函数误差 9(2=(Y2-ZY2 (11) 插值消减法 式中，S(xY为误差消减后点（￥y处的形 由于D-InSAR:技术为面状监测、单点位移精度 变值，S(xy马为D-SAR技术获得的(xy处 还不理想，而GPS技术为单点监测、位移精度高达 的形变值. 毫米量级：因此利用两种技术各自的优点进行融合 栅格函数误差插值消减后，实现了D-SAR技 处理，则可消去DSR技术监测的单点误差值， 术监测结果在有GP四监测点区域上的误差的消除. 获得连续面状监测毫米精度的滑动位移值.因此， 本研究将D-InSAR与GPS相融合的栅格函数 本研究提出D.FSAR&GPS相融合的栅格函数误 误差插值消减法应用在山R-6号滑坡上，消减掉该 差插值消减法，其基本思想是以反距离权重(nese 区域上设有GPS滥测点的区域DSAR监测结果 distance weighted Dw)插值为基础.首先，将GPS 的误差值，其结果如图10所示.可以看到误差消减 监测的单点位移值作为标准位移值，计算D-SAR 后，其垂直滑动位移值大部分都处在小于3的区

第 2期 王桂杰等:差分干涉合成孔径雷达技术在广域滑坡动态辨识上的实验研究 要来源以下几个方面 . (1)SAR传感器自身的限制产生的误差 .受 SAR传 感 器 所发 射 电磁 波 波 长的 限 制, 目前 D-InSAR技术本身进行监测所得处理结果的形变位 移值受所用波长的限制;另外, 受卫星重复周期的限 制 ,对同一地区进行相同的观测至少需要十几天甚 至几十天的时间间隔, 这使得发生较快地面活动的 区域发生严重失相关, 因此 D-InSAR技术目前只适 合监测每月毫米至厘米到每年毫米至厘米级的滑坡 活动状况. (2)失 相 关产 生 的误 差.失 相关 是 限制 D-InSAR技术应用的一个较为严重的问题, 除了由 于干涉数据的时间去相关外,由于基线长度、轨道轻 微不平行、变形运动过快、植被覆盖以及在连续获取 数据期间发生的滑坡形变过多等都会导致相位的失 相关问题,严重影响数据处理结果的精度 ,且对于数 据对上相关性较小的区域无法获得其形变信息 . (3)电磁波在大气中传播延时产生的误差.这 个延时相当于在干涉图中附加一个相位, 这个附加 相位给最终的测量结果增加了不确定性.也就是说 只要相干性满足一定条件即使实际没有形变的点 , 也有微小位移值结果 .因延时引起的位移残差一般 在所用波长量级的范围内变动 . (4)数据处理过程中处理方法的异同也会对结 果产生误差 .像去除平地效应的整平 、滤除系统噪 声的滤波以及相位解缠等处理的好坏都会直接影响 测量结果误差的大小 . 因此,要消去 D-InSAR技术监测的单点误差 , 提高监测精度, 可将其与其他技术结合处理.本研 究根据所拥有数据提出了 D-InSAR与 GPS相融合 的栅格函数误差插值消减法 , 对 L1R--6号滑坡有 GPS监测点的区域 D-InSAR监测结果的误差进行 了消减 . 4 D-InSAR与 GPS相融合的栅格函数误差 插值消减法 由于 D-InSAR技术为面状监测 、单点位移精度 还不理想,而 GPS技术为单点监测、位移精度高达 毫米量级;因此利用两种技术各自的优点进行融合 处理, 则可消去 D-InSAR技术监测的单点误差值 , 获得连续面状监测毫米精度的滑动位移值 .因此 , 本研究提出 D-InSAR＆ GPS相融合的栅格函数误 差插值消减法,其基本思想是以反距离权重(inverse distanceweighted, IDW)插值为基础 .首先 , 将 GPS 监测的单点位移值作为标准位移值, 计算 D-InSAR 与 GPS在每一点监测结果的位移差, 其计算式可表 示如下式 : Z(xi, yi, zi)=SD(xi, yi, zi)-SG(xi, yi, zi), i=1, 2, …, N (7) 式中 , N为研究区上 GPS监测点的数量 , SD (xi, yi, zi)为研究区上第 i个 GPS监测点上的 D-InSAR技 术监测位移值 , SG(xi, yi, zi)为研究区上第 i个 GPS 监测点上的 GPS监测位移值 , Z(xi, yi, zi)为研究区 上第 i个 GPS监测点上的 D-InSAR与 GPS监测技 术的位移差.利用式(7)得到有 GPS监测点的位置 处两种监测结果单点的误差函数值. 以这些单点误差函数值为标准, 进行有 GPS监 测点的监 测区域内 栅格数据 反距离加权 插值 (IDW).其插值原理公式如下 : Z(x0 , y0 , z0 )=∑ n j=1 λjZ(xj, yj, zj) (8) 式中 , Z(x0 , y0 , z0 )为点 (x0 , y0 , z0 )处的预测误差 值, n为预算计算过程中要使用的预测点周围样点 的数量, λj为预测计算过程中使用的各样点的权 重, Z(xj, yj, zj)为在点(xj, yj, zj)处获得的误差值. 权重 λj的确定和计算公式为 λj=d -p j0 ∑ n i=1 d -p i0 , ∑ n j=1 λj=1 (9) 式中 :p为指数值, 本研究取为 2;dj0为预测点 (x0 , y0 , z0 )与各已知样点(xj, yj, zj)之间的距离 .插 值后得到 D-InSAR研究区域上有 GPS监测点的监 测区域上各点的误差函数值, 假设其与位置坐标函 数关系如式(10)所示 : Z(x, y, z)=f(x, y, z) (10) 式中 , Z(x, y, z)为对应坐标点 (x, y, z)处的误差函 数值 . 在获得插值后的误差函数值后, 对 D-InSAR技 术监测结果与误差函数值 Z(x, y, z)作差, 如下式 所示 : ST(x, y, z)=SD(x, y, z)-Z(x, y, z) (11) 式中 , ST(x, y, z)为误差消减后点 (x, y, z)处的形 变值 , SD(x, y, z)为 D-InSAR技术获得的 (x, y, z)处 的形变值 . 栅格函数误差插值消减后 ,实现了 D-InSAR技 术监测结果在有 GPS监测点区域上的误差的消除. 本研究将 D-InSAR与 GPS相融合的栅格函数 误差插值消减法应用在 L1R--6号滑坡上, 消减掉该 区域上设有 GPS监测点的区域 D-InSAR监测结果 的误差值 ,其结果如图 10所示.可以看到误差消减 后,其垂直滑动位移值大部分都处在小于 3 cm的区 · 139·

。140 北京科技大学学报 第33卷 间上，水平滑动位移值大部分都处在小于8m的滑 动区间上 00046-0.050 ■0050-0.10 ■-0.503--0.400 ■-0.400--0200 ■0.100-0.130 -0.200--0.178 结 ●0.0002-0.050 -0493-0.400 ■0.050-0.100 ■0.100-0.127 -030 0.2000.001 插值格数据 0-0.011-0 ■-0.470--0.080 ■0-0.030 ■-0.080-0 00.068 ■0.30-0.122 92d的垂直位移插值消减结果图 图92d的水平位移插值消减结果图 图10D-S1R&GS融合消减掉误差后的形变图（单位：m) Fg 10 De fommation map ofD.hSAR&GPS eror reduction(unite m) 为了能够清晰地表明消减掉误差后的滑动形变 但D-ISAR差分位移单点监测精度低于GPS滥测 位移值状态，图11给出了研究区域上研究期间GS 精度.给出了D-SAR技术单点监测误差的可能 监测点处的D-SAR监测、GPS滥测和栅格函数误 原因. 差插值消减法消减掉误差后的三种情况下的滑动形 (3)提出了基于反距离权重差值的D-InSAR& 变位移值状态.图中清晰地表明了采用栅格函数误 GPS湘融合的栅格函数误差插值消减法，消除掉了 差插值消减法消减掉误差后有GPS滥测点处的滑 有GPS滥测点的监测区域上DSR技术监测结 动形变位移值与GPS滥测结果基本相同，达到了 果的误差值，提高了DSA时技术的监测精度. GPS监测量级的精度.因此，区域上其他各点的位 本实验研究表明，基于D-SR技术并使用 移值也都处在相同量级上，即消除误差后获得了 PALSAR数据可以辨识出研究区域内滑坡滑动的可 D-S4技术在有GPS监测点的监测区域上高精度 能区域确定滑动风险区域.由于对滑坡的单点位 的滑动位移值 移监测精度上还有一定的困难，对精度要求较高的 这种方法的优点是以GPS滥测值为标准，一次 区域可结合GPS等对地观测技术加以解决. 性消减掉D-InSAR技术由于各种原因产生的误差， 较适合大面积区域辨识的情况下对重点区域的 参考文献 监测 Achache J Fruneau B De lcaun C Applicability of SAR nter fe rmetry for operatonal monitoring of andslides//P rooeed ings of 5结论 he Seomnd ERS APP licationsWorkshop Londap 1995 165 2 Fruneau B Achace J DelcourtC Observation and malell ing of (1)利用D.FSAR技术获得了研究区域内详 he Sante tienne deTiee landslide usng SAR nterferometry 细的垂直和水平滑动形变位移值.通过利用滑动位 Tec0 nothysics1996265(3):181 移速率及形变大小对位移值的分级显示，表明了研 [31 Viemeer J Wagner W Dkau R Monitoring malerate slope 究区域上各部分区域的变动状态，辨识出研究区域 movements(landslides)in he southem French Alps using differ 上可能发生滑动的区域和发生滑动较大的区域，并 ential SAR interferometry/Poceedings of the 2nd hema timnal 通过与研究区域其他地质资料及光学影像对比，分 Workshop an ERS SAR In erferomnety FRNGE 99 Lige Bel gi1999 析了滑动位移较大区域的可能原因. [49 RizzoV TesauoM SAR interfermet and fel dam of Ran (2)对单体山R-6号滑坡进行了详细研究，得 dazzo landslide Eastem Sicily Italy).PhysChem Earh PartB 到滑坡上各个区域的滑动位移状态.与现场GS滥 200025(9)771 测进行了比较研究，两者监测结果的位移趋势一致， [5 Paob B MarpC GiogioF et al Use of differentialSR nter

北 京 科 技 大 学 学 报 第 33卷 间上,水平滑动位移值大部分都处在小于 8 cm的滑 动区间上 . 图 10 D-InSAR＆GPS融合消减掉误差后的形变图(单位:m) Fig.10 DeformationmapofD-InSAR＆ GPSerrorreduction(unite:m) 为了能够清晰地表明消减掉误差后的滑动形变 位移值状态 ,图 11给出了研究区域上研究期间 GPS 监测点处的 D-InSAR监测 、GPS监测和栅格函数误 差插值消减法消减掉误差后的三种情况下的滑动形 变位移值状态.图中清晰地表明了采用栅格函数误 差插值消减法消减掉误差后有 GPS监测点处的滑 动形变位移值与 GPS监测结果基本相同, 达到了 GPS监测量级的精度 .因此, 区域上其他各点的位 移值也都处在相同量级上, 即消除误差后获得了 D-InSAR技术在有 GPS监测点的监测区域上高精度 的滑动位移值. 这种方法的优点是以 GPS监测值为标准 , 一次 性消减掉 D-InSAR技术由于各种原因产生的误差 , 较适合大面积区域辨识的情况下对重点区域的 监测. 5 结论 (1)利用 D-InSAR技术获得了研究区域内详 细的垂直和水平滑动形变位移值, 通过利用滑动位 移速率及形变大小对位移值的分级显示, 表明了研 究区域上各部分区域的变动状态, 辨识出研究区域 上可能发生滑动的区域和发生滑动较大的区域 ,并 通过与研究区域其他地质资料及光学影像对比 ,分 析了滑动位移较大区域的可能原因. (2)对单体 L1R--6号滑坡进行了详细研究 ,得 到滑坡上各个区域的滑动位移状态.与现场 GPS监 测进行了比较研究,两者监测结果的位移趋势一致 , 但 D-InSAR差分位移单点监测精度低于 GPS监测 精度 .给出了 D-InSAR技术单点监测误差的可能 原因 . (3)提出了基于反距离权重差值的 D-InSAR＆ GPS相融合的栅格函数误差插值消减法, 消除掉了 有 GPS监测点的监测区域上 D-InSAR技术监测结 果的误差值,提高了 D-InSAR技术的监测精度 . 本实验研究表明, 基于 D-InSAR技术并使用 PALSAR数据可以辨识出研究区域内滑坡滑动的可 能区域,确定滑动风险区域 .由于对滑坡的单点位 移监测精度上还有一定的困难, 对精度要求较高的 区域可结合 GPS等对地观测技术加以解决 . 参 考 文 献 [ 1] AchacheJ, FruneauB, DelacourtC.ApplicabilityofSARinter￾ferometryforoperationalmonitoringoflandslides//Proceedingsof theSecondERSApplicationsWorkshop.London, 1995:165 [ 2] FruneauB, AchaceJ, DelacourtC.Observationandmodellingof theSaint- tienne-de-TinéelandslideusingSAR interferometry. Tectonophysics, 1996, 265(3):181 [ 3] VietmeierJ, WagnerW, DikauR.Monitoringmoderateslope movements(landslides)inthesouthernFrenchAlpsusingdiffer￾entialSAR interferometry//Proceedingsofthe2ndInternational WorkshoponERSSARInterferometry, FRINGE99, Liège.Bel￾gium, 1999 [ 4] RizzoV, TesauroM.SARinterferometryandfielddataofRan￾dazzolandslide(EasternSicily, Italy).PhysChemEarthPartB, 2000, 25(9):771 [ 5] PaoloB, MarioC, GiorgioF, etal.UseofdifferentialSARinter- · 140·