范雪婷李明巨张大骞赵朝贺

(1.江苏省基础地理信息中心,江苏 南京210013;2.江苏地质调查研究院,江苏 南京210018)

0 引 言

地表沉降是一项重要的地理国情信息。持续的地表沉降会导致建(构)筑物开裂、路面塌陷等灾害发生,破坏市政基础设施和生产生活设施,严重时易引发城市内涝,给城市建设和人民群众生产生活造成较大影响。目前,传统的地表沉降监测技术手段主要是GNSS测量及水准测量技术,这些方法虽然监测精度高,但是存在监测点空间密度低、投入成本大,无法提供监测全区域的细节沉降信息等问题。近几年来发展起来的时间序列InSAR(Interferometry Synthetic Aperture Radar)技术可以获取大范围、高密度、连续时间序列的地表沉降数据,逐步被应用于各大城市的区域地面沉降监测。本文主要采用多主影像相干目标小基线InSAR方法(MCTSB-InSAR)[1],该方法的主体思路是将干涉相位分析与处理焦点完全集中在那些具备稳定散射特性的像元点集上,基于这些稳定点目标进行差分干涉相位建模,继而从干涉相位中分离出单一信号分量并获得所需的形变信息[2]。本文以无锡市为例,探索MCTSB-InSAR方法在大区域地表沉降监测中的应用效果,获取无锡市2012年2月-2016年1月间的地表形变时空分布,并对监测结果进行精度评定。

1 监测方法

MCTSB-InSAR技术是在永久散射体(Persistent Scatterer Interferometry,PS-InSAR)和小基线集(Small Baseline Subset,SBAS)方法基础上进行改进的一种时间序列InSAR技术,该方法综合了PS-InSAR和SBAS两种方法的优势,克服了相应缺点。它不像PS-InSAR方法选取单一SAR影像为主影像,而是考虑满足小基线要求的SAR影像任意组合干涉像对,从而降低了对影像数量的要求。MCTSB-InSAR也无需像SBAS方法要首先完成每个差分干涉图的相位解缠,再进行时序分析,而是通过考察空间邻近的两个高相干稳定点的二次差分,同时完成形变参数解算和相位解缠[2]。另外MCTSB-InSAR方法可以实现多幅影像,特别是不同轨道沉降提取结果的无缝拼接,在跨轨道大区域沉降监测中具有较大优势。MCTSB-InSAR技术监测地面沉降算法原理主要包括组合干涉像对、稳定点目标提取、建立形变模型、非线性形变量估计、多幅InSAR地表沉降结果拼接5个方面。

1.1 基于特征区域相干性组合干涉像对

传统小基线干涉像对组合方法主要为时空基线阈值法,该方法生成差分干涉图后需要人工筛选,对干涉相位质量较差的干涉对进行人工剔除,导致工作量非常繁重。本文基于区域相干特性组合干涉像对,首先选取特征区域,通常选择包含建筑物较多的高相干城区,特征区域范围应较小,以利于快速处理。经反复试验,一般为500×500像元较为适宜。然后对所有干涉图的相干性进行定量评估,将平均相干系数由高到低排序,按影像数2-3倍自动筛选出排序在前的干涉像对,最后生成整个研究区域的干涉图。与传统通过设置时空基线硬阈值方法相比,该方法可实现优质干涉像对的自动、完备提取。

1.2 三级阈值稳定点目标提取方法

单一阈值法用于稳定点目标的识别具有局限性。本文利用振幅、相干系数、振幅离差指数信息三级阈值探测法进行稳定点目标提取[3]。首先设置绝对平均相干系数阈值,选取高相干点HCP1;设置相对平均相干系数阈值和平均幅度阈值T A,选取那些相干性相对较高,并且平均幅度又很大的高相干点HCP2;然后根据相干目标的稳定性,设置幅度离差指数阈值D T,选取稳定点目标HCP3;依据公式HCP=(HCP1∪HCP2)∩HCP3,将获取的点目标作为最终稳定点目标[4]。

1.3 基于区域Delaunay三角网的线性形变速率与高程误差估计

1.3.1 Delaunay三角网建立

构建合理的稳定点目标网络是探测线性形变的关键步骤之一,MCTSB-InSAR技术利用局部Delaunay三角网连接方法[1],将整景影像分割成相互重叠的网格,在各网格中使用Delaunay三角网进行连接。三角网中一条边就代表相邻点目标构成的一个点目标对,由此构建的三角网既简洁,又增加了连接的冗余度,可有效避免整体网络断裂成不连通的子网络。

1.3.2 线性形变速率与高程误差参数计算

通过局部Delaunay三角网连接所有的点目标,考察第k个差分干涉图上i、j两顶点之间的差分干涉相位:

式(1)中,wrap{}表示相位缠绕算子;系数a1由该点的入射角、斜距以及雷达波长λ决定,a2=是垂直基线;Δh i,j为高程误差之差;t k为干涉图的时间基线;Δv i,j为两点间的线性形变速率之差;残余相位增量,为两点间非线性相位、大气相位和噪声相位之差的和。

考虑三角网两点间非线性形变相位差决定性较小,大气影响为低频信号,在空间上存在一个相关距离,因此设置三角网两点间距在大气影响一定距离内(一般1～2km),大气相位差可以忽略,以及各高相干点的单点噪声很小,满足在该条件下,利用整体相位的相关系数最大值进行优化处理,建立如下目标函数:

式(2)中,J表示单位复数;γ为整体相位相干系数,其最大值在[0,1]之间,γ值越大表示模型估计值对差分相位的拟合程度越好,Δv i,j和Δh i,j的估值越可靠。

待网络上的边全部实现最大化求解后,对模型相干系数进行最大值检验计算,判断哪些连接边计算出来的结果是可靠的[5]。一般设定适当阈值(本文设为0.7),将大于该阈值的模型相干系数最大值所对应的边作为可靠的边。

1.3.3 三角网集成

以某一具有已知形变量和DEM误差的高相干点为参考点,对筛选出的优质三角网连接边,建立联合大型方程组。对连接边上解算出的速率和高程误差增量进行积分,获得各相关点的形变速率和高程误差的绝对值[6]。

1.4 非线性形变量估计

从原始差分干涉图相位减去高程误差和线性形变速率分量可获得残余相位,主要包含非线性形变,大气影响及噪声等相位。非线性形变速率由于空间上不相关性,因此为高频信号,而在时间上为低频信号;大气相位在时间上具有不相关性,为高频信号;噪声相位体现为在时间和空间上都不相关的随机高频信号[7]。因此,可以利用残余相位三个分量的不同频率特征,采取适当的频域滤波方法,将三者分离出来。将得到的非线性形变结果与线性形变结果叠加,从而得到时间序列的累计形变量[1]。

1.5 多幅大区域InSAR地表沉降结果拼接

大范围监测区域需要多幅影像才能实现全覆盖,但是各景影像获取时间、影像数目都不一致,得到的地表沉降结果具有差异性,因此实现多景影像沉降结果的无缝拼接是一大挑战。MCTSB-InSAR方法从测量平差角度出发,通过提取影像间重叠区域的同名点结果,采用最小二乘平差方法算出各影像结果的改正数,实现多幅点目标结果的一次性无缝拼接。

2 实验及结果分析

2.1 实验数据介绍

以无锡市作为实验区。选取加拿大C波段的Radarsat-2宽幅雷达影像作为沉降监测的数据源,其单景影像覆盖面积为150 km×150 km,影像分辨率约为30 m,共需两幅影像实现无锡市全覆盖(图1)。每幅各有25期影像,数据获取时段为2012年2月—2016年1月。采用30 m分辨率SRTM DEM数据去除地形贡献项。

图1 无锡市Radarsat-2影像覆盖范围

2.2 数据处理

借助中国测绘科学研究院张永红课题组研发的InSAR地表形变监测系统(GDEMSI)进行数据预处理和时序SAR影像的地表形变反演工作。利用特征区域相干性法组合干涉像对,选取包含建(构)筑物较多的500×500像元大小的城镇地区作为特征区,共生成63个最佳干涉像对(图2)。

图2 小基线对时空分布图

设置绝对平均相干系数阈值0.72,相对平均相干系数阈值0.65和平均幅度阈值1.2,幅度离差指数阈值0.3,根据三级阈值法探测高精度稳定点目标。在此基础上,创建Delaunay三角网差分相位函数模型,求解获得线性形变信息和DEM误差。最后依据残余相位分量的不同频率特征分离出非线性形变相位,与线性形变叠加获取该幅影像完整形变值。在获取两幅监测点结果后,通过提取影像间重叠区域的同名点结果,采用最小二乘平差方法对两幅点目标结果进行拼接,最终形成区域协调一致的沉降结果。

2.3 结果分析

经处理,最终获取无锡市489 781个点目标,平均每平方千米监测点数达106个,可见监测点分布密度之大。从无锡市2012年2月—2016年1月沉降速率分布,可以看出,在监测时段内,无锡市的沉降相对较为集中,主要分布在江阴南部和锡山北部地区,涉及江阴市云亭街道、周庄镇、徐霞客镇、长泾镇、华士镇以及锡山区东港镇等。全市最大沉降速率为28 mm/a,四年最大累积沉降量达94.6 mm,位于江阴市周庄镇;沉降速率超过10 mm/a的面积达97.3 km2。笔者收集到的无锡市2007—2012年的InSAR监测结果,对比发现无锡市沉降范围略有扩大,沉降速率有所趋缓。

无锡市是中国重要的纺织品制造基地和出口基地,至2016年末,全市共有规模以上纺织服装企业772家,每年生成大量的纺织品。通过现场调研发现,无锡沉降区存在大量的纺织企业,由于纺织行业是用水较多的行业之一,据估算,平均工业产值每100元就消耗3 m3的水,需水压力较大[9]。因此,长期过量开采地下水是无锡市地面沉降的主要原因。2000年,江苏省通了《关于在苏锡常地区限期禁止开采地下水的决定》,2000-2005年无锡市共封井1 100眼,区域地质环境发生了显著改善。但根据InSAR监测结果,江阴南部和锡山北部地区地面沉降仍相对较为严重,推测深部地下水禁采以后,江阴的工厂纷纷转为开采浅部地下水,导致沉降范围扩大。

图3 无锡市年沉降速率分布图(2012年2月—2016年1月)

3 精度评定

为了验证InSAR监测结果的精度,搜集了无锡市30个水准点数据,2012—2015年每年12月份施测一次。因InSAR监测点和水准点位置的不完全一致性,根据约定临近点原则对初始获取的水准点进行了筛选:即以实测点位为中心,如果在一定距离范围内(本文设为80 m)存在至少一个InSAR点目标,则选择该实测点参与精度评定。实验中利用19个水准点数据对InSAR监测得到的无锡市地表平均沉降速率进行比对验证分析[8](表1)。

表1 无锡市水准与InSAR沉降速率对比/mm·a-1

从表1可以看出,二者差值绝对值最大值为4.7 mm/a,误差均方差为1.9 mm/a。该误差实质上还包括两种数据时间不一致性和空间不一致性导致的误差,可见MCTSB-InSAR方法的理论测量精度应高于该值。

为了进一步分析InSAR时序沉降监测结果的可靠性,选取4个典型水准监测点(CJDL33、CJ-123、CJ078、I-145),利用2012-2015年的时序实测成果与InSAR监测点的时序沉降量进行对比分析(图4)。水准点CJDL33位于东港镇沉降中心,CJ-123位于祝塘镇沉降中心,两个水准点表现为不断下沉;CJ078位于惠山区堰桥街道的京沪高速铁路沿线附近,呈轻微抬升趋势;I-145位于江阴青阳镇,水准点和InSAR监测点表现都较为稳定,形变量在+5～-5 mm范围内波动。通过对比发现,无论是位于沉降中心还是抬升区亦或稳定区,二者的监测结果具有较高的一致性。

综上分析,利用时间序列InSAR技术监测得到的年平均形变速率和时序累积形变量具都有较高的监测精度,实验结果能够反映了区域的基本沉降特征。

图4 无锡市时序InSAR沉降量与水准沉降量曲线对比图

4 结 语

(1)基于区域相干特性组合干涉像对,实现了优质干涉像对的自动、完备提取,极大提高了干涉图质量与生成效率。

(2)本文既顾及相干点回波信号的高信噪特征,又兼顾相干点的稳定性,采用振幅、相干系数、振幅离差指数信息三级阈值探测法实现了稳定点目标的高精度、大密度提取。

(3)本文以 Radarsat-2为数据源,利用MCTSB-InSAR方法实现了无锡市2012年2月-2016年1月的地表形变的沉降信息反演。结果表明,在该阶段内,无锡市沉降主要发生在江阴南部和锡山北部地区,涉及江阴市云亭街道、周庄镇、徐霞客镇、长泾镇、华士镇以及锡山区东港镇等,最大沉降速率为28 mm/a,发生在江阴市周庄镇。利用水准实测数据开展精度评定,二者具有较高一致性,验证了MCTSB-InSAR方法利用Radarsat-2宽幅雷达影像获取大区域地表形变信息的可靠性。