基于精细DEM的InSAR大气相位改正试验研究
2013-01-04葛大庆范景辉
李 曼,夏 耶,2,葛大庆,张 玲,范景辉,王 艳
(1.中国国土资源航空物探遥感中心,北京 100083;2.德国地球科学研究中心,波茨坦 14473)
0 引言
我国是世界上滑坡灾害最严重的国家之一。特别是库区滑坡,滑坡变形破坏直接危害坡体上人员的生命财产安全;滑坡体入江减少水库的有效库容,形成堰塞湖并堵塞河道,同时滑坡体高速失稳入江后所激发的次生涌浪对航运、沿岸居民或建筑物也造成威胁[1]。因此,库区滑坡的变形破坏及稳定性问题已成为制约水电开发的主要工程地质问题之一。开展库区滑坡长期监测逐渐成为防止和减轻库岸滑坡灾害的重要手段之一。传统上对滑坡的监测和预警主要是采用位移计、应力计、水位计等现场仪器观测,或利用GPS和TDR等远程电子设备监测,但这些技术方法很难应用到植被茂盛、人员难以进入地区的大面积滑坡的监测上[2],且因监测成本高、监测网络密度低,对变形缓慢的滑坡监测灵敏度不是很理想。
自20世纪90年代中期开始,合成孔径雷达干涉测量(InSAR)技术以其监测成本低、监测网络密度高且可以监测到整个滑坡体的形变全貌等独特优势,而被应用于对缓变形滑坡体形变的监测中。但由于早期的星载雷达数据分辨率多在20~30 m之间,重复周期一般为20~50 d,对于面积不大的滑坡体,其时间去相关现象非常严重;另一方面,大气波动、尤其是对流层湿度和温度的变化造成大气延迟在差分相位图上表现出非均一性,过大的大气延迟相位甚至会“淹没”研究区域的形变相位,使得雷达干涉技术在滑坡监测领域远远没有取得实际应用的价值。随着TerraSAR-X和Cosmo-SkyMed等星载高分辨率雷达数据的出现,以及相应高分辨率数字高程模型数据的成功获取,短时间跨距雷达数据的时间去相关现象已迎刃而解[3],而大气波动延迟成为雷达干涉技术监测缓变形滑坡短时间跨距形变场的主要制约因素。如何从干涉相位或差分相位中抑制或者分离出大气延迟相位以提高滑坡体待解算形变量的准确性,成为监测缓变形滑坡亟待解决的关键问题。
目前,在差分干涉合成孔径雷达(D-InSAR)技术的应用过程中,消除或抑制大气延迟效应的技术大致可以分为2类:①在数据处理过程中,引入外部大气数据,如实测地表气象参数、大气数值模型、地基GPS水汽反演结果以及MODIS水汽反演结果等,这些数据都能提供大气层中水汽含量信息,将水汽含量转换为大气相位,进而从差分干涉相位中消除[4-5];②充分考虑大气在时间域呈高频、在空间域呈相对低频的特点,对多时相SAR数据进行处理,采用SAR数据集进行自身大气校正,进而达到抑制或分离、消除大气干扰相位的目的。线性组合法、干涉图叠加法、随机滤波法以及长时序D-InSAR等方法在实际数据处理过程中经常被采用,其中永久散射体干涉测量技术(permanent scatterer interferometry,PSInSAR)[6-7]和短基线集技术(small baseline subset,SBAS)是长时序 D-InSAR技术的典型代表,也是目前最为先进的SAR数据集自身校正法。可见,现有消除SAR图像中大气延迟的技术方法都能在一定程度上削弱大气扰动引起的相位贡献,但它们又都有各自的适用条件和自身无法克服的缺点,在实际应用过程中都受到这样或那样不同程度的限制,如MODIS和MERIS的水汽反演数据不能在夜晚以及白天云量较多的情况下发挥作用;线性组合法和干涉图叠加法是假设地表形变符合线性模型;随机滤波法需要先确定相邻的稳定区域等[4,7]。
但在地形起伏较大、降水丰沛的高山峡谷区用InSAR技术进行地表形变监测时,如果无法获得外部大气数据,而雷达数据量又较少的情况下,采用大气湿延迟相位与地面高程之间的相关模型削弱大气扰动相位的方法逐渐得到了国内外学者的关注。Delacourt等[8]在应用InSAR技术监测意大利Etna火山时指出,大气对流层中湿延迟与大气层的温度和水蒸汽的分压之间满足一定的函数关系。Xia等[9]分析大气效应引起的延迟相位依赖于研究区的气候和高程,如果2次获取雷达图像时的气温基本相同,则由这2幅图像得到的差分干涉图中大气延迟效应产生的相位很小;反之,如果2次获取雷达图像时的研究区温度差别较大,则干涉图中大气相位与研究区的高程相关。蒋延臣[10]基于湿延迟与高程的线性回归模型校正于田地震ScanSAR差分干涉纹图的大气效应,在一定程度上改善了干涉纹图的质量,但在DEM梯度变化较大的地区,大气扰动的去除效果仍不是很好。基于此,本文以三峡库区秭归县树坪滑坡区为研究对象,根据大气湿延迟相位与地面高程的相关性,通过提取树坪滑坡体以外区域高相干点上的高程值及解缠后的差分干涉相位,从数学角度寻求差分干涉相位(大气湿延迟相位)与高程之间的最优函数模型。根据这一模型,由树坪滑坡体的DEM得出该区域的大气效应相位,将其从干涉图中去除,最后保留的相位主要是由地表形变引起。这种改正大气效应的方法,可有效地改善干涉图的质量,提高雷达干涉测量技术的形变监测精度。
1 数据源
1.1 高分SAR数据
本文采用德国DLR的Terra SAR-X聚束模式的高空间分辨率的雷达数据进行差分干涉处理。Terra-SAR卫星重访周期为11 d,载波使用X波段,波长3.1 cm。该卫星定位精度高,空间基线控制能力强,为利用传统InSAR方法测定滑坡体短期内的小尺度形变量和形变场的非均匀分布提供了重要保障[11]。
1.2 高程数据
高程数据通常用于补偿InSAR数据处理过程中的高程相位。本文所用DEM数据是美国莱卡公司生产的机载激光雷达系统ALS50-II获取的,经过了数据解算、系统误差检校、航带系统误差检查修正及激光点云分类、质量控制及DEM精度评价等处理,其高程绝对精度一般为10~20 cm,平面精度通常约0.5~1 m。由于客观条件限制,三峡库区大部分地形起伏很大,其平面精度会相应略有变化,但高程精度通常变化较小[12]。
图1 时间跨距11 d的差分干涉相位图Fig.1 Differential interferometric phase map with 11 d interval
2 数据处理及相位特性分析
2.1 InSAR数据处理
采用两轨法提取雷达目标的相位信息,获取地表形变场。将重轨获取的2景雷达图像精确配准并生成干涉图后,利用外部高精度DEM数据和雷达卫星成像的几何信息模拟生成高程相位,从干涉相位中减去平地相位、高程相位,此时得到的差分干涉纹图(图1)中任一像元x的相位φdiff(x)可视为由形变相位Δφd(x)、大气延迟相位Δφatm(x)、因DEM数据不准确引起的残余地形相位Δφε-h(x)、因空间基线不准确引起的残差相位Δφb(x)以及系统热噪声等引起的误差相位 Δw(x)组成[7,13],即
式中W{}表示缠绕相位。
2.2 差分干涉相位特性分析
差分干涉纹图中各类型相位的量级大小以及在干涉图中的分布特征与研究区的气候条件、选取数据的性能指标以及InSAR技术的处理方法等密切相关。
本研究选取的TerraSAR雷达卫星分辨率为1 m,其基线精度高(通常小于10 cm),影像覆盖范围为10 km×5 km。因此,由空间基线不准确产生的残差相位Δφb(x)会很小。
因DEM数据不准确引起的残余地形相位Δφε-h(x)也很小。这是因为,研究中采用了1 m分辨率的DEM数据对高程相位进行了补偿。假定DEM数据误差为1 m,根据雷达微波的基本参数以及SAR干涉测量成像的几何关系,通过
分别计算出图1两幅相位图中心点处由高程误差引起的干涉相位值分别为-0.03 rad(图1(a))和0.07 rad(图 1(b))。式(2)中:Δφε-h为残余地形相位;B⊥为2次获取数据时雷达卫星所处位置的垂直基线长度;RL为卫星与图像中心点之间的距离;H为卫星与星下点之间的距离;Δε为高程误差;θ0为卫星入射角;λ为雷达波的波长。
通过计算发现,高程误差相位在整幅干涉图中的值基本不变。说明差分干涉图中的相位主要并不是由高程误差引起的。这一点也可以从另一角度进一步印证:由式(2)可知,高程误差引起的相位与垂直基线B⊥长度呈正比,如果差分干涉图中的相位主要由高程误差引起,那么垂直基线B⊥=-84.96 m的差分干涉图(图1(b))中的相位值约为垂直基线B⊥=35.20 m的差分干涉图(图1(a))中相位值的2.41倍,但图1(a)(b)两幅差分干涉纹图中的相位值并不存在这种倍数关系。于是,可以肯定地说,在差分干涉相位图中,高程误差及空间基线不准确引起的相位值很小;同时,系统热噪声和数据处理误差产生的相位虽然无法避免,但也相对较小,基本可以忽略不计。这样,图1差分干涉图中的相位主要由地表形变和大气扰动引起,即
对树坪滑坡区地表形变的长期监测结果表明,只在树坪滑坡体的局部区域存在缓慢形变迹象,其他区域都相对比较稳定。因此,图1(a)中大范围明显的相位基本都是由大气延迟引起的,树坪滑坡区的形变相位完全被大气延迟相位所“淹没”。但由于雷达影像沿距离向覆盖范围仅10 km,雷达卫星入射角的变化也较小(约0.6°),入射角变化引起的大气扰动相位也不明显,干涉图中这一明显的大气相位基本是由大气层本身的局部相关波动引起。在距离雷达卫星较远区域(图1(a)右侧),相位值基本不变,这可能与相应区域的大气层厚度变化较小或者基本没有变化有关;在距雷达卫星较近区域(图1(a)左侧),大气延迟相位梯度变化明显,且与入射角无关,随着高程的增加,大气层厚度变化明显。
3 大气延迟相位的改正方法
以树坪滑坡体附近稳定点CR08(见第4节图2)为相位解缠基准点,采用最小费用流算法对差分干涉图进行相位解缠。解缠后的相位仍是由大气相位和形变相位组成,只有将大气相位去除后,才能得到研究区沿雷达视线向的地表形变信息。
由于大气厚度大、成分复杂且垂直成层分布,根据大气对重复轨道干涉测量的影响,可以将地表和遥感器之间的大气分为电离层、中间大气(平流层和中间层)和对流层3部分,其中,对流层中的空气对流很明显,云、雾、雨等现象都发生在这一层,水蒸气也几乎都在这一层内。对流层又可以细分为干延迟和湿延迟2部分,干延迟由非水蒸气部分的大气延迟产生,在时域内比较稳定,在空域内具有大尺度变化的特性;湿延迟是由云、雨等水蒸气引起,水蒸汽分布极不均匀,且与温度密切相关,其变化在时间域呈高频,在空间域呈相对低频[11]。
根据大气各组分在时域、空域中的变化特性不同,雷达影像经干涉、差分及相位解缠后,差分干涉相位中的大气延迟相位主要由对流层中的湿延迟引起,其他层对形变相位的影响都已经基本消除。对流层中的湿延迟取决于研究区的温度和水蒸气压力,而温度、水蒸气压力与所处区域的高程密切相关[8]。因此,湿延迟相位与相应区域高程呈间接相关,依此可以有效消除解缠相位中大气湿延迟引起的干扰相位。另外,研究区山体的坡度、坡向相对高程来说,对湿延迟的影响相对很小,且扰动程度和影响范围也不尽相同,因此本文不考虑坡度及坡向对湿延迟相位的影响。
大气湿延迟在空间上的相关距离相对较小,在整幅差分干涉图中(图1(a))的分布特征完全不同,左侧大气相位梯度大,且随高程升高而呈规律性变化;右侧相对较小,甚至与位置也存在一定的相关性。于是根据不同区域大气相位的变化特征不同,分别提取相应区域内的湿延迟相位、高程信息和沿距离向的坐标,基于最小二乘法,对湿延迟相位与高程、沿距离向的坐标值进行曲线或曲面拟合,在保证拟合曲线或曲面在高程(0,800)m范围内没有明显的转折,残差都处于[-2,2]rad之间,且在误差平方和最小的前提下,建立相应区域的大气效应校正模型。最后,将模拟的大气相位从解缠相位中去除,即可得出树坪滑坡沿雷达视线向的形变信息。
4 树坪滑坡形变信息提取
树坪滑坡的变形历史较长,1996年以前,坡体局部就出现过变形,从2003年6月三峡水库蓄水后,滑坡变形加剧,在滑坡体的上、中、下部都出现了明显的地面裂缝和房屋开裂现象。目前,树坪滑坡体局部区域仍有比较明显的变形,但变形速率相对较小,每年仅数10 cm,尤其在三峡水库蓄水及水位消落期间,再加上大气降水的影响,滑坡体局部区域变形速率较大[14]。本文采用高分辨率TerraSAR-X数据,通过传统的雷达干涉技术提取树坪滑坡形变信息。
图2是树坪滑坡在2012年1月份、时间跨距为11 d的差分干涉相位图。由于在2012年1月2日至1月13日期间,树坪滑坡区以小雨或阴天为主,且气温变化明显,因而图2显示大气湿延迟引起的扰动相位较大,甚至掩盖了树坪滑坡体的形变信息。在图2中Ⅰ区域,大气湿延迟相位仅与该区的高程密切相关,随着高程的升高,大气湿延迟相位变化明显;而在Ⅱ区域及树坪滑坡体区域,湿延迟相位不仅与研究区的高程有关,还与雷达与目标间的距离(即沿距离向的坐标值)有关,这就需要分别建立大气湿延迟效应校正模型。
图2 树坪滑坡时间跨距为11 d的差分干涉相位图(20120102—20120113)Fig.2 Differential interferometric phase map of Shuping landslide with 11 d interval
在图2Ⅰ区域内,由相干系数≥0.3点上的湿延迟相位和高程值得出该区域的最优函数校正模型,即根据Ⅰ区域的DEM数据,由式(4)即可模拟出该区域所有点的大气湿延迟相位,且大气残差相位都介于[-1.5,1.0]之间。
但在图2的II区域及树坪滑坡区,湿延迟相位不仅是高程的二次函数,同时也是图像距离向坐标X的线性函数,且湿延迟在这2个区域的分布特征基本一致,由Ⅱ区域高相干点(相干系数≥0.3)上的湿延迟相位、高程值HⅡ及距离向坐标X,通过曲面拟合,模拟出Ⅱ区域的最优函数校正模型为式(5),即
树坪滑坡区的湿延迟相位分布特征也遵循式(5)的表达规律。因此,基于Ⅱ区域及树坪滑坡区内所有点上的DEM值及沿距离向的坐标值X,由式(5)就可模拟出Ⅱ区域及树坪滑坡体上所有点上的大气湿延迟相位。湿延迟相位的残差相位Δφresi,Ⅱwet介于[-2.0,1.0]之间。在对图2差分干涉相位解缠时,解缠相位参考点CR08位于Ⅰ和Ⅱ区域的共同区域A中,解缠后A区域内湿延迟相位值接近于0。通过湿延迟相位校正模型(4)(5),分别得出A区域的模拟大气相位值差很小,取2次模拟相位的平均值作为A区域的湿延迟相位,即可得到图2整个区域平滑的湿延迟相位分布状态,将其从解缠相位中去除,进而恢复出了树坪滑坡体在2012年1月2日至1月13日之间11 d内的形变信息(图3)。
图3 树坪滑坡11 d内的形变信息(mm)Fig.3 Deformation of Shuping landslide with 11 d interval(mm)
从图3可以看出,尽管树坪滑坡周围某些局部区域大气湿延迟相位并没有完全去除,但图中红、黄区域仍十分清楚地框出了树坪滑坡体形变场在11 d中的位置、大小及其分布,且形变场是非均匀的。红色区域为最大下沉区,集中分布在树坪滑坡上部,沿雷达视线向最大沉降值达-4.4mm;黄色区域的形变值相对较小,在滑坡体中部及坡体前缘都有分布,树坪滑坡这种局部变形特性与所处地区的气候及三峡水库水位消落密切相关。2012年1月初,秭归小雨连绵不断,降雨对滑坡的诱发作用从滑坡体上部开始,由于坡体后缘存在较宽裂缝,雨水入渗使滑带岩土体发生软化、泥化作用,造成滑带的抗剪强度降低。另一方面,雨水在坡体内积聚形成较大的静水压力,削弱滑动面上的有效正应力,导致坡体抗滑力减小,这都会造成滑坡体上部局部区域出现较大变形。三峡水库水位消落直接诱发滑坡体前缘出现明显形变,三峡水库于2011年12月31日开始启动新一轮生态补水,2012年1—2月,三峡水库下泄量按6 000 m3/s左右控制。可见,在这一时期三峡水库水位一直在持续消落,而坡体组成物质的透水性弱,坡体内孔隙水压力消散速度滞后于库水位的降低速度,使得坡体前缘动水压力增大,进而导致滑坡体的下滑力增加。树坪滑坡体局部区域虽有宏观变形现象,但其内并没有形成贯通的滑动面,目前滑坡整体仍处于稳定状态。
5 结论
基于树坪滑坡区特殊的地理位置、复杂的地形地貌特征以及多雨潮湿的气候条件,本文通过分析大气湿延迟在空间的分布特征,分别建立空间局部区域大气效应的最优函数校正模型,将大气湿延迟相位从解缠相位中去除,最终恢复出树坪滑坡体的形变场。
需要说明的是,这种大气改正方法并没有考虑研究区山体的坡向和坡度角变化对湿延迟产生的影响;同时又假设了对流层中湿延迟相位在空间局部区域内具有相关性。任何函数模型都不能完全确切地描述湿延迟相位的空间分布规律,无法准确无误地将差分干涉图中的湿延迟相位完全去除,在树坪滑坡体的形变场中仍有湿延迟效应的影响,但由于树坪滑坡体局部区域变形较大,通过这种去除大气效应的方法,清楚地框出了树坪滑坡体形变场在11 d中的位置、大小及其分布,较为真实地恢复出了树坪滑坡体的形变信息,这有利于推动D-InSAR技术在三峡库区大面积缓变形滑坡形变监测及稳定性调查的广泛应用。
志谢:感谢中国国土资源航空物探遥感中心郑雄伟提供了高精度DEM数据精度评价资料。
[1] 李远耀.三峡库区渐进式库岸滑坡的预测预报研究[D].武汉:中国地质大学,2010:1-30.Li Y Y.Research on prediction and forecast of progressive bank landslide in the Three Gorges Reservoir[D].Wuhan:China University of Geosciences,2010:1-30.
[2] 王桂杰,谢谟文,邱 聘,等.D-InSAR技术在大范围滑坡监测中的应用[J].岩土力学,2010,31(4):1337-1344.Wang G J,Xie M W,Qiu P,et al.Application of D-InSAR technique to landslide monitoring[J].Rock and Soil Mechanics,2010,31(4):1337-1344.
[3] 夏 耶,范景辉,李 曼,等.高分辨率雷达数据三峡库区滑坡监测技术[C]//第十八届中国遥感大会论文集,2012:161-169.Xia Y,Fan JH,LiM,etal.Monitoring technique of Shuping landslide with high-resolution Radar images in Three Gorges Reservoir area,in China[C]//18th China Symposium on Remote Sensing,2012:161-169.
[4] 吴云孙,李振洪,刘经南,等.InSAR观测值大气改正方法的研究进展[J].武汉大学学报:信息科学版,2006,31(10):862-867.Wu Y S,Li ZH,Liu JN,et al.Atmospheric correction models for InSAR measurements[J].Geomatics and Information Science of Wuhan University,2006,31(10):862-867.
[5] 范青松,汤翠莲,陈 于,等.GPS与InSAR技术在滑坡监测中的应用研究[J].测绘科学,2006,31(5):60-62.Fan Q S,Tang C L,Chen Y,et al.Applications of GPS and InSAR in monitoring of landslide studies[J].Science of Surveying and Mapping,2006,31(5):60-62.
[6] Ferretti A,Prati C,Rocca F.Nonlinear subsidence rate estimation using permanent scatterers in differential SAR interferometry[J].IEEE Transactions on Geoscience and Remote Sensing,2000,38(5):2202-2212.
[7] 范景辉.基于相干目标的D-InSAR技术地表形变监测研究与应用[D].北京:中国科学院研究生院,2008,136-139.Fan JH.Research on coherent targets based on D-InSAR technique and its application to surface deformation monitoring[D].Beijing:Chinese Academy of Sciences,2008,136-139.
[8] Delacourt C,Briole P,Achache J.Tropospheric corrections of SAR interferograms with strong topography,application to Etna[J].Geophysical Research Letters,1998,25(15):2849-2852.
[9] Xia Y,Kaufmann K,Guo X F.Landslide monitoring in the Three Gorges area using D-InSAR and corner reflectors[J].Photogrammetric Engineering and Remote Sensing,2004,70(10):1167-1172.
[10] 蒋延臣.星载宽幅合成孔径雷达干涉测量形变监测理论与应用[D].武汉:武汉大学,2010:82-90.Jiang Y C.Study on scanning synthetic aperture radar interferometry theory and application for deformation monitoring[D].Wuhan:Wuhan University,2010:82-90.
[11] 李小凡,李 颖,曾琪明,等.应用与ASAR同步的MERIS对重复轨道InSAR进行大气校正[J].北京大学学报:自然科学版,2009,45(6):1012-1018.Li X F,Li Y,Zeng Q M,et al.Correction of atmospheric effects on repeat-pass interferometric SAR using MERIS and ASAR synchronous data[J].Acta Scientiarum Naturalium Universitatis Pekinensis,2009,45(6):1012-1018.
[12] 刘圣伟,郭大海,陈伟涛,等.机载激光雷达技术在长江三峡工程库区滑坡灾害调查和监测中的应用研究[J].中国地质,2012,39(2):507-517.Liu SW,Guo D H,Chen H T,et al.The application of airborne lidar technology in landslide investigation and monitoring of Three Gorges Reservoir area[J].Geology in China,2012,39(2):507-517.
[13] 葛大庆,王 艳,范景辉,等.地表形变D-InSAR监测方法及关键问题分析[J].国土资源遥感,2007,19(4):14-22.Ge D Q,Wang Y,Fan JH,et al.A study of surface deformation monitoring using differential SAR interferometry technique and an analysis of its key problems[J].Remote Sensing for Land and Resources,2007,19(4):14-22.
[14] 汪发武,张业明,王功辉,等.三峡库区树坪滑坡受库水位变化产生的变形特征[J].岩石力学与工程学报,2007,26(3):509-517.Wang FW,Zhang Y M,Wang G H,et al.Deformation features of Shuping landslide caused by water level changes in Three Gorges Reservoir area,China[J].Chinese Journal of Rock Mechanics and Engineering,2007,26(3):509-517.