基于InSAR时序分析的北京地区地表沉降监测研究

李红梅,陶秋香

(1.青岛滨海学院,青岛 266510;2.山东科技大学,山东省基础地理信息与数字化技术重点实验室,青岛 266590)



基于InSAR时序分析的北京地区地表沉降监测研究

李红梅1,陶秋香2

(1.青岛滨海学院,青岛 266510;2.山东科技大学,山东省基础地理信息与数字化技术重点实验室,青岛 266590)

针对常规大地测量技术在城市地表沉降监测过程中时空分辨率较低的的局限性,本文利用PS和SBAS两种InSAR时序分析技术对城市地表形变进行监测。论文首先对SBAS-InSAR技术的原理进行详细推导和分析,在此基础上,对覆盖北京地区的29景ENVISAT-ASAR数据采用PS和SBAS技术进行处理,分别提取了该区域地表年平均形变速率图。通过与已发表成果对比发现,PS和SBAS在大范围、持续缓慢城市地表形变监测中能够探测一致的地表形变运动趋势;并且在数据量有限的情况下SBAS技术能够获得更为可靠的监测结果。

PS;SBAS;InSAR时序分析;城市地表形变监测

0　引　言

随着国民经济的迅速发展,城市地下水开采量也急剧增加,由此引发了一系列的城市地表沉降问题,严重制约着经济与环境的协调和可持续发展。北京地区拥有2000万人口的特大型都市,城区人口密集,有关资料显示该地区城市地表沉降越来越严重。因此,北京地区开展城市地面沉降监测与防范相关的技术研究具有重要意义。

合成孔径雷达差分干涉测量技术(D-InSAR)[1]以其监测周期短、覆盖范围广以及成本低的优势被广泛应用于较大区域的城市地表沉降监测中。常规D-InSAR技术侧重研究短时间间隔的单次形变,无法获取研究区域地面沉降随时间的演化情况,而且其测量精度易受时空失相关和大气延迟等因素的影响[2]。为了克服这些困难,米兰工业大学的Ferretti与Prati等人基于干涉点目标分析的思想,提出了永久散射体(PS)法[3],并成功应用在Ancona地区的滑坡监测中。国内对PS-InSAR技术应用在城市地表形变监测中的研宄工作起步较晚,最早由西南交通大学陈强、刘国祥等人于2006年在上海地区成功的利用PS技术提取了上海地面沉降场。随着研究的不断深入,学者们发现PS-InSAR技术监测精度比较依赖于主影像的质量,并且基线过长容易导致几何失相干,因此该方法往往受到数据量较少的限制。为此Berardino等人于2002年提出了短基线集(SBAS)方法[4],该方法选择基线较短的干涉对组成不同集合,集合之间基线较长,采用最小二乘法得到每个基线集合的地表形变时间序列,再利用奇异值分解法将多个集合联合起来求解。不要求具有同一主影像,因此可以充分利用SAR数据生成较多的干涉图,实现有限数据量条件下的高精度反演。

本文在对PS和SBAS技术进行分析的基础上,以北京地区为试验区,对覆盖该试验区2007～2010年的29景ENVISAT-ASAR数据进行处理,获取该区域地表年平均形变速率图,分析了引起北京沉降的主要原因。通过对比PS和SBAS所得的结果进行表明, 两种时序分析方法在大范围、持续缓慢城市地表形变的监测中可以获取一致的地表形变运动趋势;通过与已发表的成果对比可以发现,在数据量有限的情况下SBAS技术监测的形变结果更可靠。

1　InSAR时序分析技术原理

1.1PS-InSAR技术原理

2000年,Ferretti等[5]为克服常规D-InSAR技术受大气延迟、时空去相关和DEM误差等因素的限制,提出了PS-InSAR技术。该技术利用二维线性(或非线性) 相位回归分析模型,对长时间间隔内依然保持高相干的永久散射体(PS)的差分干涉相位进行回归分析来获取形变信息。

现有N+1景SAR影像。其中一景影像为公共主影像,另外N景影像为辅影像,将N景辅影像分别与公共主影像进行配准并做干涉处理。利用外部DEM做干涉处理后,得到N干涉图,其干涉相位为

φdiff=φdef+φtopo+φatm+φflat+φnoise,

(1)

式中: φdef为两次成像期间内的地表形变产生的相位; φtopo表示由外部DEM误差引起的地形残差相位; φflat表示由卫星轨道误差导致的残余平地相位; φatm为两次成像时雷达信号穿过不均匀大气层时引起的延迟相位; φnoise为噪声相位。式中前两项是主要相位贡献值,后面三项为残余相位,表示为δ φres,则:

φdiff=φdef+φtopo+φres

(2)

假设Δd=v·T,其中v为形变速率,则有

φdef=k1v·T+k2·Δherror+φres,

(3)

其中:T为两次成像时间间隔;φres为残余相位; Δherror为高程模糊数;k1,k2为系数。

从上述公式可以看出,相位差分模型为二维线性相位模型,通过多次回归分析可以得到地表形变的线性速率和残余DEM误差及其各类噪声相位。

1.2SBAS技术原理

2001年,Berardino和Lanari等[6]提出了小基线集(SBAS)InSAR技术。该技术主要解决PS技术中时间上采样过于稀疏的问题[7]。

设研究区有N+1幅SAR影像,获取时间依次为t0,t1,…tn.根据时空基线条件选择其中一幅影像作为公共主影像,其他N幅影像作为辅影像,得到M个干涉像对。假设在tA和tB时刻获得的两景SAR影像进行差分干涉处理后生成第j个差分干涉对,则差分干涉相位δφj可以表示为

δφj(r,x)=φ(tB)-φ(tA)

(4)

式中:d(tA)和d(tB)表示tA和tB时刻相对于参考时刻t0(假设d(t0)≡0)的累积形变信息;φ(tA)和φ(tB)表示相应相位值。

现用向量φT=[φ(t1),…,φ(tN)]表示N个相位图,向量δφT=[δφ1,…,δφM]表示M个干涉相位图,假设主图像获取时间先于辅图像获取时间,则:

δφj=φtISj-φtIMj,j=1,…,M,

(5)

式中:IM表示干涉像对中主图像;IS表示干涉像对中辅图像。

定义一个含有M个方程,N个未知参数的方程组,其表达式表示为

δφ=Aφ,

(6)

式中: 矩阵A[M×N]每一行对应一个干涉对,A[j,ISj]=1和A[j,IMj]=-1,其余元素为零。矩阵A是取决于干涉图组合方式的一个近似关联矩阵。当所有干涉对属于同一个基线集时,A是一个列满秩矩阵。当M=N时,方程组(9)有唯一解; 当M>N时,方程组是一个超定方程,此时在最小二乘约束下可以求得其唯一解

(7)

当所有干涉对属于不同的基线集时(设有L(L>1)个基线集)，A的秩为N-L+1,秩亏,为此,采用奇异值分解的方法计算矩阵A的广义逆矩阵A+,进而求解方程组(6)的最小二乘最小范数解,其过程为

对矩阵A进行奇异值分解

A+=US+VT,

(8)

式中,U为M×N的正交矩阵,其前N列是ATA特征向量,称为矩阵A的左奇异矩阵,V为N×N的正交矩阵,它的所有列ATA是的特征向量,称为矩阵A的右奇异矩阵,S为M×N矩阵。则:

A+=VS+UT,

(9)

则方程组(6)的最小二乘最小范数解为

(10)

式中: ui和vi分别是各自矩阵U和V的列向量。

将方程组(6)中对相位的求解转化为对平均相位变化速率的求解问题,则待求参数向量为

(11)

则有:

δφ=Bv.

(12)

式中: B为M×N的矩阵,其中B[i,j]=tj+1-tj, (ISi+1≤j≤IMi,i=1,…,M),其它元素值为零。对B进行奇异值分解,解出各时间段平均速度v.

2　实验分析

2.1实验数据

本文研究区如图1所示,覆盖整个北京地区,中心地理坐标为北纬40.02°、东经116.25°,覆盖面积约为4 200km2.数据选用29景Track490的ENVISATASAR降轨数据,时间跨度为2007年1月至2010年10月。实验过程中还使用了荷兰Delft大学发布的DORIS精轨数据和SRTM的 90m分辨率的外部DEM.

图1　研究区域及周边环境示意图

2.2实验处理与结果

1) 数据预处理

试验收集的SAR影像是原始雷达数据,进行时序分析之前需要对数据进行预处理,生成时序分析可以处理的干涉图。主要包括原始数据成像、图像配准、生成干涉图主要步骤。干涉图是数据处理的基础,生成干涉图的前提是两景影像的精确配准,数据处理中采用相干系数法[8]逐一像素进行配准,保证配准精度。

2)PS时序分析

PS-处理过程中主要包括以下几个关键环节,分别是:PS点选择、相位解缠、轨道误差去除、地形相位去除,DEM误差相位的估算与去除、大气误差相位的估算与去除等。

其中,永久散射特性的像素点(PS点)是PS时序分析法的处理对象。实验采用Hooper博士提出的方法对PS点进行选取。研究区覆盖城区,建筑物较为密集,干涉效果比较好。试验共选取了302 387个点,其密度大约为34个/km2.

确定PS点后,采用基于Delaunay不规则三角网法对缠绕的相位进行解缠[9]。解缠的方式是在选取的PS点上构建三角格网,对三角格网按照时间域上的相位差和空间域上的相位差的顺序逐步进行解缠。根据公式(1)显示,每一个PS点的相位值包含地表形变相位、DEM误差相位、大气延迟和轨道误差相位。为了提取地表形变相位值,需对其他相位成分加以估计并剔除。

在选取PS点时,计算出DEM误差与垂直基线成比例关系,通过计算比例值大小判断影像中对应的像素是否足够稳定,进而考虑其是否可以确定为PS点。在选取PS点的同时也估算出了DEM误差相位。试验中使用分辨率为90m、高程精度为16m的SRTM-DEM数据。为了消除大气效应对数据处理的影响,2001年Ferretti假设大气相位和轨道误差相位是空间相关的,并且在时间域上和空间域上对解缠相位分别进行高通滤波和低通滤波来估算个影像的大气效应误差。图2示出了PS处理过程中估计的各类误差相位示例。

图2　DEM误差相位、大气延迟和轨道误差相位图示例：(a) DEM误差相位；(b) 主影像大气延迟和轨道误差相位

3)SBAS时序分析

与PS时序分析数据处理不同,SBAS时序分析的首要步骤是选取合适的差分干涉像对。实验中选取了满足空间基线小于800m、时间基线小于692天、相关系数大于0.65三个条件的110个干涉像对,组合方式如图3所示。将110景干涉像对使用相干系数法进行精确配准,并做复共轭相乘得到原始干涉图。

图3　小基线集干涉图组合方式

SBAS时序分析技术中选取的相位稳定的点称为缓慢失相关滤波相位像元点(SDFP),SDFP像元点是基于相位特征来定义的。试验中共选取了290 364个SDFP点,其密度大约为33个/km2.

确定出SDFP点后,实验采用与PS技术相位解缠相同的方法进行相位解缠,获取解缠后差分相位。建立如公式(9)的矩阵方程,并求解得到地表形变量和形变速率。图4示出了SBAS-InSAR处理过程中估计的各类误差相位的示例。

经过以上一系列的时间序列分析处理,最终得到的研究区年平均雷达视线向(LOS) 形变速率场,如图 5所示。其中得到的形变量都是相对于主图像20 080 908里面的PS/SDFP点的形变量进行分析得到的。形变速率图中,红色点位(负值)表示背离卫星的方向运动,即该点的运动趋势是沿雷达视线向方向下降;蓝色点位(正值)表示朝着卫星的方向运动,即该点的运动趋势是沿雷达视线向方向上升。为了更好的显示研究区的沉降情况,分别将PS法与SBAS法得到的高相干点的年平均沉降速率文件导入ArcGIS,利用Kriging插值方法得到整个研究区的沉降图如图5所示。

图4　SBAS-InSAR技术估算的DEM误差相位、大气延迟和轨道误差相位图示例(a) DEM误差相位；(b) 主影像大气延迟和轨道误差相位

图5　研究区2007～2010年期间的沉降图(a) PS法; (b) SBAS法

2.3结果分析

如图5(a)表示PS技术获取的研究区地表年平均沉降速率图,该区域年平均沉降速率为21.2～-51.6mm/yr.北京城区范围内呈现红色,表现为地表下沉的运动趋势;周边区为蓝色,表现为地表抬升的运动趋势。图5(b)表示SBAS技术获取的研究区地表年平均沉降速率图,该区域年平均沉降速率为31.7～-109.9mm/yr.同样的该监测结果显示城区范围内呈现红色,表现为地表下沉的运动趋势;周边区为蓝色,表现为地表抬升的运动趋势。两种技术监测出的该区域地表沉降保持一致。

仔细分析后发现,该区域沉降由中部到东部逐步加重的原因有两点:1) 该区域主要是冲积扇的中下部地区,土质较细,过度抽取地下水会导致地面沉降; 2) 随着一些重工业的东移,地表构筑物的建设以及地下工程建设加深了该区域的地表沉降。

西部以及周边区域地表略有抬升的原因是该区域相干性低,受噪声污染严重。仔细分析发现:1) 该区域植被覆盖较密集,这一点造成了数据处理过程中选取PS点和SDFP点过程产生误差,增加了后续处理过程中噪声的产生。2) 该区域地处山区,地形起伏比较大,地形相位估计需要高分辨率的DEM数据,但本次试验中由于各种原因仅收集到了90m分辨率的数据,精度较低,这一点也增加了数据处理的噪声。

2.4精度分析

由于种种原因,未能获取研究区外部的水准数据或GPS数据,因此本文将PS和SBAS技术所得结果与已发表的文献所述结果进行对比,表1给出了两种方法所得结果与文献[10]的对比情况。

从表1可知,本文利用SBAS技术获取的北京地区地表形变监测结果相对于PS技术的结果在数值上更接近已发表文献的监测结果,主要有以下两个原因:1)SAR数据量较少,PS-InSAR技术在识别PS点时对数据量有一定的要求,实验中由于资料有限仅收集了29景数据,影响了PS-InSAR的监测结果;2)SBAS-InSAR技术将29景SAR数据组合成若干个集合,形成了109景时间和空间基线均较小的差分干涉图,弥补了PS技术依赖于单一公共主影像质量的缺陷,从而获取了相对更可靠的监测结果。

表1　两种方法所得的北京地区年平均沉降速

3　结束语

本文利用PS和SBAS技术分别对覆盖北京地区的29景ENVISAT ASAR数据进行了处理,获取了该地区2007年至2010年期间的沉降时间序列图,发现研究区的东部沉降比较严重,部分地区年平均沉降速率达到90 mm/a以上。同时分析了整个研究区的沉降情况及形成原因,得到造成严重沉降的主要原因为地下水的过量开采。PS和SBAS法所得年平均速率标准差以及与已发表文献的对比均表明, SBAS技术相对于PS技术监测地面沉降能够得到更为可靠的结果,与实际情况比较相符。

[1] GABRIEL K, GOLDSTEIN M, ZEBKER H A. Mapping small elevation changes over large areas:differential radar interferometer[J].Journal of Geophysical Research, 1989, 94(B7): 9183-9191.

[2] 许文斌,李志伟,丁晓丽,等. 利用MERIS 水汽数据改正ASAR干涉图中的大气影响[J]. 地球物理学报,2010, 53(5): 1073-1084.

[3] FERRETTI A, ROCCA F, PRATI C. Permanent scatterers in SAR interferometry [C]//Proceedings of IGARSS 1999, Hamburg, Germany, 28 June-2 July, 1999.

[4] BERARDINO P, FORNARO G, LANARI R,etal. A new algorithm for surface deformation monitoring based on small baseline differential SAR interferograms[J]. IEEE Transactions on Geosciences and Remote Sensing, 2002, 40(11), 2375-2383.

[5] FERRETTI A, PRATI C, ROCCA F. Permanent scatterers in SAR interferometry[J].Geoscience and Remote Sensing, 2001, 39(1):8-20.

[6] BERARDINO P, FORNARO G, LANARI R,etal. A new algorithm for surface deformation monitoring based on small baseline differential SAR interferograms[J].Geoscience and Remote Sensing, 2002, 40(11):2375-2383.

[7] LANARI R, MORA O, MANUNTA M,etal. A small-baseline approach for investigating deformations on full-resolution differential SAR interferograms[J].Geoscience and Remote Sensing, 2004, 42(7):1377-1386

[8] 陈富龙,张红,王超.高分辨率SAR影像同名点自动匹配技术[J].中国图象图形学报,2006,11(9):1276-1281.

[9] 吴涛,张红,王超,等.多基线距DInSAR技术反演城市地表缓慢形变[J].科学通报,2008,53(15):1849-1857.

[10]李永生,张景发,李振洪等.利用短基线集干涉测量时序分析方法监测北京市地面沉降[J].武汉大学学报(信息科学版),2013,38(11):1374-1377.

Research of the land Deformation Monitoring in Beijing Based on the InSAR Tme Series Analysis Method

LI Hongmei1,TAO Qiuxiang2

(1.QingdaoBinhaiUniversity,Qingdao266510,China; 2.GeomaticsCollege,ShandongUniversityofScienceandTechnology,Qingdao266590,China)

Aiming at the limitations of low spatial and temporal resolution of the traditional geodetic measurements in the process of urban surface subsidence monitoring. Two InSAR time series analysis methods were introduced into the surface subsidence monitoring. The paper firstly derived and analysed the basic principles of SBAS-InSAR technique detailedly. On this basis, we processed and analysed 29 scenes of ENVISAT-ASAR images covering the area using the PS-InSAR and SBAS-InSAR technique, getting the annual average rate of deformation of the area. Compared with published resultsshows that PS and SBAS are of consistency in deformation tendency of the wide range and slowly urban surface subsidence and in the case of a limited quantity of the data SBAS technique can obtain a more reliable monitoring result.

PS; SBAS; InSAR time series analysis; urban surface subsidence monitoring

2015-12-11

山东省优秀中青年科学家科研奖励基金项目(编号:BS2013SF013)

P642.26

A

1008-9268(2016)03-0078-06

李红梅(1978-),女,硕士,讲师,研究方向为短时交通流预测及深基础变形监测预测。

陶秋香(1977-),女,博士,副教授,研究方向为微波遥感技术与应用、InSAR关键技术、数据处理等。

doi:10.13442/j.gnss.1008-9268.2016.03.016

联系人： 陶秋香 E-mail: qiuxiangtao@163.com