基于InSAR的重庆东港地区地表形变监测

李华蓉,戴双璘

(重庆交通大学 智慧城市学院,重庆 400074)

大型机械长期在高填方区域作业导致了重庆东港集装箱码头区域地表形变较为严重,引发了东港区域支座墩柱出现裂缝以及轨道移位严重等现象,对工作人员的生命财产安全构成了极大威胁。传统地表形变监测方法(如水准测量监测、GNSS监测等)获取形变信息精度高,但作业周期长,劳动强度大,仅能获取少量离散测点的形变信息,且测点分布稀疏、易被破坏、难以实时更新,具有一定局限性。

合成孔径雷达干涉测量(interferometric synthetic aperture radar, InSAR)技术起源于20世纪60年代中后期,集成了合成孔径雷达技术和微波干涉技术[1]。20世纪90年代,D-InSAR技术在InSAR技术基础上发展而来,但D-InSAR技术极易受大气扰动和时空失相干影响,无法应用于高精度、变化缓慢的地表形变监测,只适用于地震等形变剧烈的地表形变监测[2]。为克服上述缺点,时序InSAR技术应运而生,时序InSAR技术主要以永久散射体干涉测量(PS-InSAR)技术和小基线集干涉测量(SBAS-InSAR)技术为代表[3-4],能实现低成本、高精度、高密度的大范围地表形变监测。

近年来,学者们对时序InSAR技术进行了深入研究。王义梅等[5]基于PS-InSAR和SBAS-InSAR技术,处理了2012—2013年的15景TerraSAR-X影像,提取了郑州市的沉降场,并结合地质条件及已有的研究成果进行分析,表明时序InSAR技术在大范围地表形变监测应用的有效性;杨国创等[6]联合时序InSAR和GNSS技术对深圳和香港海岸线的地表形变进行监测,将基于SBAS-InSAR获得的矢量反演结果与基于PS-InSAR技术获得的反演结果进行比较分析,研究表明:这二者结果相近,形变结果均表明位于填海区的地铁建筑及地基较浅的建筑物均存在明显沉降现象;卢旺达等[7]利用PS-InSAR技术处理了2015—2018年的24景Sentinel-1A数据,获取了天津地区地面沉降结果,并利用SBAS-InSAR技术进行验证,这两种方法得到的累积形变量一致性高,表明该结果可为天津市灾害防治提供数据支持。上述研究均表明:时序InSAR技术在地表形变监测领域均能取得较好效果。

基于此,笔者利用Sentinel-1数据和改进SBAS-InSAR 技术对重庆东港集装箱码头区域进行地表形变监测,将基于PS-InSAR技术获得的永久散射体(persistent scatter, PS)点作为地面控制点,并引入SBAS-InSAR技术处理流程,与人工选取地面控制点的SBAS-InSAR技术进行对比,并将两种方法得到的形变结果与二等水准测量数据进行精度评定,探究改进了SBAS-InSAR技术的可靠性。

1 数据集与研究方法

1.1 研究区概况与数据集

研究区域地处重庆市南岸区,该区域地跨东经106°44′—106°46′,北纬29°34′—29°36′,最低海拔为150 m,最高海拔为327 m,面积约为38.84 km2。研究区域的重庆东港集装箱码头位于长江江岸,地势平坦,植被覆盖较少。该区域为高填方区,长期有大型机械在地面作业,引发地面振动,导致填土不断被压实,地表发生形变。在自然、人为因素的双重作用下,该研究区在运营期间曾出现了缆车斜坡道盖梁支座墩柱脱空现象,岸侧盖梁与支座墩柱中心线偏差较大,部分墩柱出现裂缝,靠近江侧裂缝较多,轨道移位严重等。研究区地表形变较为严重,对缆车道运行造成较大安全隐患,需定期对该研究区进行地表形变监测。研究区域及东港集装箱码头区域如图1。

图1 研究区域Fig. 1 Research area

笔者选取31景Sentinel-1升轨、VV极化方式、单视复数、干涉宽幅模式的SAR卫星影像,时间跨度为2018年1月4日—2019年12月25日,平均影像入射角为33°18′18″,每景影像时间间隔为24 d。此外,笔者还利用定位精度优于5 cm的Sentine-1的卫星精密轨道(precise orbit ephemerides, POD)数据修正轨道信息,去除因轨道误差引起的系统性误差。利用30 m分辨率的SRTM-1数据作为外部DEM数据,去除平地相位和地形相位。

1.2 研究方法

在SBAS-InSAR处理流程中,须使用GCP文件修正干涉相位和解缠后的相位。GCP是被认为形变为零的地面控制点,GCP选取应该远离形变区域、解缠错误相位跃变区域,应考虑平地区域及相干性较高区域。传统SBAS-InSAR技术通常是操作人员根据GCP点的特性,手动选取GCP点,结果具有不稳定性且依赖操作人员经验判断。

考虑PS-InSAR技术可从影像中提取出在各个时间都能保持高相干性的PS点。PS点常用的提取方法有振幅离差指数法和相干系数阈值法[8],振幅离差指数DA的计算如式(1)。

(1)

式中:μ为振幅标准差;σ为振幅均值。

通过振幅离差指数来判断相位标准差大小,振幅离差指数越大,相位标准差越小[9]。相干系数阈值法是通过限制相干系数阈值对PS点进行筛选,相干系数位于0～1,越接近1,相干性越好。

为提高基于SBAS-InSAR技术获取的地表形变精度,减小人工选取GCP点带来的误差,笔者利用PS点选取特性与GCP选取特性一致,先对研究区基于PS-InSAR技术进行处理,获取稳定PS点,再通过限制振幅离差指数、相干性指数、年平均沉降速率这3个指标阈值,从PS点中筛选出满足条件的点作为GCP,引入SBAS-InSAR技术处理流程。改进的SBAS-InSAR技术处理流程如图2。

图2 改进的SBAS-InSAR技术处理流程Fig. 2 Process flow chart of the improved SBAS-InSAR technology

2 数据处理

1)数据预处理。使用精密轨道数据对SLC数据轨道信息进行修正;对原始SLC进行多视处理,得到强度信息图;基于地理坐标系对强度图进行裁剪,提取研究区域。

2)生成连接图。根据空间基线为临界基线的2%,时间基线为120 d的原则建立干涉对连接图。共生成117对干涉对,其中最大时间基线为120 d,最大的空间基线为95.95 m,干涉对基线的时间-空间分布如图3。图3中:空心点代表主影像,实心点代表辅影像,线条代表干涉对。

图3 干涉对基线的时间-空间分布Fig. 3 Time-spatial distribution of interference versus baseline

3)差分干涉处理。按时间基线和最小原则选取2018年5月28日的影像作为超级主影像[10],辅影像基于超级主影像进行配准,由每组像对的复数影像共轭相乘生成干涉相位图[11]。利用外部DEM去除干涉图中的平地相位和地形相位,采用最小费流法(minimum cost flow, MCF)进行相位解缠处理,解缠相干系数阈值设置为0.2,以去除相干性较低的像元。对该数据进行距离向为4方位向为1的多视,获得20 m制图分辨率的InSAR产品,图4为20180128、20180404像对的干涉结果。对生成的117对干涉结果进行人工检视,删除相干性较低、解缠效果较差的干涉对。

图4 20180128、 20180404像对的干涉结果Fig. 4 The result of interferometric graph between 20180128 and 20180404

4)进行轨道精炼和重去平。轨道精炼与重去平的目的是估算和去除解缠后的相位图中依然存在的残余相位[12]。笔者采用改进的SBAS-InSAR技术,基于PS-InSAR技术获取研究区内的PS点,通过三阈值判别法对PS点进行筛选,选取振幅离差指数大于等于3.2、相干性指数大于0.75、年平均沉降速率位于[-1,1]之间的PS点,共计28个。利用该GCP 点对基线参数进行重新定义,计算出相位偏移,修改相位解缠图像头文件中的轨道参数。基于PS-InSAR技术获取的GCP分布情况如图5。

图5 基于PS-InSAR技术获取的GCP分布Fig. 5 GCP distribution map obtained by PS-InSAR technology

5)进行SBAS反演估算与地理编码。SBAS反演估算分为两次。第一次反演是对轨道精炼和重去平后的干涉图进行残余地形和形变相位的估算。第二次反演是通过在时间域上进行高通滤波和在空间域上进行低通滤波计算,计算和去除大气相位,再使用上述GCP剔除残余相位。此时获取的形变结果在雷达坐标系下,需将反演估算形变结果转换到WGS-84坐标系下。基于改进SBAS-InSAR技术获取的地表形变如图6。

图6 改进SBAS-InSAR技术获取的形变Fig. 6 Deformation map obtained by the improved SBAS-InSAR technology

3 精度评定

笔者对研究区进行了同时段的二等水准监测,将获得的该区域形变数据作为真值,对上述基于SBAS-InSAR技术反演的形变数据进行精度评定。同时,与手动选取地面控制点反演出的变形数据进行对比研究,验证上述方法的可靠性和先进性。

3.1 二等水准监测结果

本研究对A区域(东港集装箱码头)进行二等水准测量,获取了2018年1月—2019年12月间的形变数据。共布设48个水准点,具体数据如表1。需注意的是:水准测量数据是基于重庆独立坐标系获得的,因此需要在ArcGIS中利用特征点对水准点进行空间校正将其转换到WGS 84坐标系下,随后才能与SBAS-InSAR技术反演的形变数据进行精度评定,空间校正后的水准点如图7。

表1 水准测量获取的形变数据Table 1 Deformation data obtained by leveling measurement

图7 空间校正后的水准点Fig. 7 The levelling point after spatial correction

3.2 手动选取GCP方法的反演结果

为探究基于PS-InSAR技术选取GCP方法的可靠性,笔者同时结合相位解缠图、相干系数图和Google earth历史影像,根据GCP选点原则手动选取28个GCP,分布情况如图8。与图5相比,手动选取的GCP位置更加集中,均位于公路和建筑物范围,而基于PS-InSAR技术获取的GCP分布大部分位于建筑物,少部分位于裸地,总体分布更加均匀。

图8 手动选取的GCP分布Fig. 8 Manually selected GCP distribution map

基于手动选取的GCP反演研究区地表形变的结果如图9。与图6相比,未改进SBAS-InSAR技术获取的形变结果最大形变量和最小形变量相差无几,但4个形变严重区域的平均沉降量较小。

图9 未改进SBAS-InSAR技术获取的形变Fig. 9 Deformation map obtained by the unimproved SBAS-InSAR technology

由于InSAR技术获取形变量为雷达视线方向(light of sight, LOS)的一维形变量,即地表各方向形变量在雷达视线方向上的投影。需先进行形变场解算,获取垂直方向形变量,再与水准测量获取的形变量进行精度评定[13]。根据文献[14],列出视线向形变量垂直向转换计算,如式(2)。

d=Δr/cosθ

(2)

式中:d为垂直向形变;Δr为视线向形变;θ为雷达入射角。

SBAS-InSAR技术获取的垂直形变量如表2。

表2 改进前后获取的垂直形变数据Table 2 Vertical deformation data obtained before and after improvement

3.3 对比分析

将基于改进SBAS-InSAR技术获取的垂直形变和基于未改进SBAS-InSAR技术获取的垂直形变分别与水准测量数据进行对比。采用平均绝对误差EMA和均方根误差ERMS作为精度评定指标[15],EMA和ERMS的计算如式(3)、式(4):

(3)

(4)

式中:d为SBAS-InSAR技术反演的垂直形变量;d1为水准测量获得的形变数据。

利用S-W test对绝对误差进行正态性检验。精度评定表、S-W检验结果和频率分布直方图如表3、表4和图10。

表3 精度评定Table 3 Accuracy evaluation

表4 S-W检验结果Table 4 S-W test results

图10 频率分布直方图Fig. 10 Frequency distribution histogram

4 地表形变结果与分析

笔者分别基于PS-InSAR技术筛选GCP和手动选取GCP的SBAS-InSAR技术获取研究区地表形变场。由结果可知:研究区形变量约为-60～30 mm,研究区共存在4个形变严重区域,其中:A区域为重庆东港集装箱码头,B区域为重庆铁路枢纽东环线拌和站,C区域为东港家园,D区域为重庆永翔现代物流产业园。A区域地表形变情况复杂,最大形变量达-52 mm,最大形变速率达-28 mm/a-1;C区域为环状形变,外环形变约为-15 mm,内部较为稳定,几乎没有形变;B、 D区域呈漏斗状形变场,最大形变值为-30 mm。结合Google Earth历史影像(图11)和相关资料可知:A区域为高填方区域,大型机械作业及货物运输使填土受压不均,导致地表沉降现象严重。2018—2019年,B、D区域正处于施工期,实际形变情况与漏斗状形变结果一致;C区域的东港家园周边处于施工期,与环状形变结果相符。

图11 Google Earth历史影像Fig. 11 Google Earth historical imagery

由精度评定结果可知:改进后SBAS-InSAR技术的EMA和ERMS分别为6.27、7.85 mm,未改进SBAS-InSAR技术的EMA和ERMS分别为6.67、 8.57 mm。这两种方法的ERMS都满足文献[14]要求,改进后SBAS-InSAR技术的EMA和ERM更小。利用S-W检验法对两种方法的EMA进行正态性检验,由检验结果可知:这两种方法的EMA均满足正态分布,其中,改进后SBAS-InSAR技术的EMA的S-W检验P值更大,更符合正态分布。

5 结 语

笔者分别将基于PS-InSAR技术筛选的GCP和手动选取的GCP引入SBAS-InSAR技术处理流程,对两种方法反演地表形变进行结果分析和精度评定。研究结果表明:SBAS-InSAR技术可作为一种有力的新型监测手段进行地表形变监测,利用PS-InSAR技术筛选GCP方法比手动筛选GCP方法更加自动化,获取的形变结果精度更高。

但本研究仍存在不足之处:如缺乏坐标转换参数,无法将水准点准确地转换至WGS-84坐标系,在空间上造成误差;SAR影像覆盖时间段与水准测量时间段不能完全一致,在时间上造成误差。