利用PALSAR数据提取地面沉降短基线集方法研究

王成

(重庆市勘测院，重庆 400020)

利用PALSAR数据提取地面沉降短基线集方法研究

王成∗

(重庆市勘测院，重庆 400020)

短基线集DInSAR方法可以很好地克服数据相干性差的困难，本文在分析短基线集方法的基础上，对大气相位和轨道误差相位的影响进行处理得到较为精确的形变相位。采用天津实验区从2007年1月17日～2009年10月25日的17景PALSAR单视幅图像提取地面相对沉降信息，将DInSAR与水准测量数据进行融合得到外部精度为3.8mm的沉降结果，同时确定了天津地面的沉降漏斗的分布和绝对形变量。

合成孔径雷达干涉测量；短基线集；地面沉降；轨道误差；大气相位误差

1 引 言

合成孔径雷达差分干涉测量技术(Differential SAR Interferometry，DInSAR)是在InSAR技术基础上进一步发展而来的，它通过比较地表形变前后干涉图之间的相位差来监测地表目标的位移[1]。具有全天时、全天候、广覆盖、高分辨率及高地表垂直形变灵敏度等优点。但在长时间微小地表形变监测中，时空去相干和大气相位延迟的影响制约了它的发展[2]。

Berardino[3]等提出的短基线集算法，可以很好地解决DInSAR长时间序列干涉纹图中缺乏大面积连续高相干区域的问题。它是在常规DInSAR基础上限制时空基线的大小后自由组合干涉对，利用短基线干涉图提取高相干区域，从而有效地进行形变反演。

本文收集了天津地区自2007年1月17日～2009年10月25日的16景PALSAR单视复影像，基于GAMMA软件平台进行常规二轨差分处理得到136组干涉像对。由于是对长时间地面形变进行监测，因此引入了短基线集方法，分别选取其中干涉质量好的垂直基线较短的78组干涉像对采用短基线集DInSAR方法进行形变反演，通过基线精化处理和拟合最佳轨道面的方法来削弱轨道误差影响，采用高斯低通滤波的方法来减弱大气相位的影响，通过阻尼最小二乘算法分离平均线性形变速率和DEM残余误差的改正量，得到最终的垂直形变量。

最后利用短基线集DInSAR方法所得到的垂直形变量与7个已知水准点同时段的沉降量进行融合，使得短基线DInSAR方法监测得到的相对形变量转化为绝对形变量。PALSAR数据监测结果的精度为3.8 mm，高相干点位密度为666个/km2。将监测结果绘制成沉降等值线并叠加到Google Earth地图上，由此分析确定天津地区的地面沉降漏斗中心位置分布及其绝对沉降量。

2 实验区概况

地面沉降是许多沿海城市面临的顽症，天津地区的地面沉降现象在中国具有代表性，其主要原因是地下水过量开采。由于地下流体资源过量开采和构造沉降活动，在宝坻断裂和蓟运河断裂以南约7 300 km2的广大地区，出现不同程度的地面沉降现象，南部和滨海地区尤为明显，并且与河北省沉降区连成一片，在这一大面积的沉降范围内，又形成了市区、塘沽、汉沽、大港及海河下游工业区形成了5个二级沉降中心。最近几年随着区县、乡镇经济的发展，出现了武清、津南、静海等新的沉降中心，其发展速度比当年老的沉降区域还要快，局部沉降速率超过80 mm/y[4]。

常规地面沉降监测技术主要有水准测量和GPS测量方法，但二者的空间点位稀疏、周期长的缺点大大限制其对区域地面沉降发展趋势的研究。DInSAR技术具有大范围、全天候、高精度、高分辨率地测量地表三维空间的微小变化，在地表变形的监测方面显示出前所未有的优越性。

3 短基线集DInSAR方法与原理

设有按时间序列t0，……，tN获取的N+1幅单视复数SAR影像，将它们以任意影像为主影像进行配准，设定垂直基线阈值(一般小于400 m)，垂直基线小于该阈值的影像作为一个基线集，并将其划分为L组，每组内的影像进行差分干涉处理，最终可得到M幅差分干涉图，且有:

现以t0时刻为参考时刻，任意时刻ti(i=1，…，N)相对于参考时刻的差分相位φ(ti)为未知数，数据处理获取的差分干涉相位δφ(tk)(i=1，…，M)为观测量，不考虑大气延迟相位和数字高程模型(DEM)误差的影响，对第k幅(k=1，……，M)差分干涉相位图中的任意像元(x，r)可以组成如下方程:

上式中λ为雷达波长，d(ti，x，r)和d(tj，x，r)分别为ti和tj时刻像元(x，r)相对于参考时刻t0的雷达视线方向(line of sight，LOS)地表形变。将式中的φ(ti，x，r)变为相邻影像像元沿LOS方向的平均相位变化速率v，有:

其中IE和IS分别为按时间顺序排列的主影像序列和从影像序列，满足IEj＞ISj，∀j=1，……，M化简即有:

其中B是M×N矩阵，B为奇异矩阵，可采用阻尼最小二乘算法进行求解(Levenberg，1944)，可得到线性形变速率为:

式中γ是阻尼因子，阻尼因子γ越大，分辨就越低；另一方面，当γ增大时分母就越大，从而减小了模型参数的误差[5]，阻尼因子γ对反演结果有一定影响，本文对阻尼因子进行实验分析选取了一个较为合适的数值，阻尼因子为0.25进行阻尼最小二乘的反演计算。

将式(3)带入式(2)得:

4 数据预处理与误差消除

4.1 实验区数据准备与预处理

针对天津实验区，选取了16景ALOSPALSAR从2007年1月17日～2009年10月25日高分辨率模式单视复影像，所获取为升轨时期的数据。将其中任意两幅图进行组合形成120个像对，使用GAMMA软件利用常规二轨差分干涉测量方法进行处理，对于PALSAR还需要进一步精化基线来削弱线性相位的影响。利用美国航天飞机雷达地形测绘任务(SRTM)3″分辨率(90 m)DEM去除地形相位和地平相位，利用最小费用流方法对干涉相位进行相位解缠，从而得到差分干涉图。

处理得到的差分相位由形变相位、残余轨道误差相位、DEM残余误差相位和大气相位等组成。由In-SAR原理可知垂直基线长度与高程模糊度成反比，基线越短，高程模糊度越大，地形对相位的影响也就越小，所以需要选择时空基线较短的干涉像对组合。通过实验分析，最终选取了4 000 m以下的干涉质量较好的78组干涉像对，采用短基线集方法进行形变反演分析，对以上误差相位进行分离并予以削弱，得到最终的形变相位。所选取的干涉对的时空基线分布如图1所示。

图1 PALSAR时空基线分布图

4.2 卫星轨道误差的消除

SAR卫星轨道误差是指用来计算的卫星星历中存在误差引起的基线计算不精确而导致后续的干涉相位误差。这里首先采取精化基线的方法对轨道误差进行一定的改正，如图2所示。为了进一步削弱残余轨道误差的影响，然后通过线性拟合的方法来对解缠后的干涉相位进行拟合，得到最佳拟合轨道面，即残余的轨道误差，并将其从干涉图中去除。按照如下方程来计算每个干涉图的最佳的轨道面[6]:

式中，[x，y]是每个像素点的在像素坐标系下的坐标，u，v是斜率参数，w是截距，如图3所示。

图2 PALSAR基线精化处理

图3 PALSAR轨道误差改正

4.3 消除DEM残余误差

此外，在实际处理过程中需要考虑DEM残余误差的影响，则式(5)表示为:

Dv+C·ε=△φ(8)

式中cT=[(4π/λ)(B⊥1/r sinθ)，…，(4π/λ)(B⊥N/ r sinθ)]，λ为雷达波长，B⊥为垂直基线，r是卫星到地面间的距离，θ是卫星视角。通过阻尼最小二乘算法计算出DEM残余误差的改正量ε，其大小在[-10，10]m范围内，结果如图4所示。同时可得到沿雷达视线方向的线性形变速率v乘以时间就得到雷达视线方向的形变量再除以cosθ，从而转化为垂直形变量。

图4 PALSAR DEM残余误差改正量

4.4 大气相位影响

短基线集方法是采用线性模型，为了进一步消去残余相位中的大气相位和非线性形变相位的影响，可以在线性模型的基础上进行适当的滤波处理[7]。根据大气相位空域的低频特性，本文采用在空间上实施高斯低通滤波处理，移除“形变信号”中上的高频部分。如图5所示，是对其中一组数据进行滤波的前后对比图。

图5 PALSAR滤波前后对比图

5 结果分析

5.1 获取天津地面相对形变量

此次形变反演过程中相干阈值设置为0.67，得到高相干点个数为599 728，影像面积约为900 km2，即点密度为666个/km2。DInSAR测量结果是相对于基准点而言的，该基准点是研究区内所有其他点形变量的参考基准，所以需要选择形变较为稳定的点作为基准点。同时选择干涉质量较好的2007年7月20日主影像为沉降基准时刻，将另外的16景影像与之差分，从而得到16幅天津实验区地面相对形变量变化图，如图6所示。图中“+”表示的参考基准点，是研究区内相对稳定的点，负值表示沉降。通过地理编码，将沉降结果由制图坐标系转换到WGS84坐标下，可以获取实验区内自2007年01月17日～2009年10月25日的累积相对形变量图，如图7所示。图中“★”表示参考基准点，“∗”表示累积沉降量最大的点，我们定义为A点，其相对形变量不超过15 cm，将其相对于参考基准点的沉降趋势表示出来，如图8所示。

图6 天津实验区地面相对形变量变化图

图7 2007年01月17日～2009年10月25日天津地面相对形变量累积图

图8 A点相对于基准点的沉降趋势图

5.2 DInSAR与水准数据的融合

通过短基线集DInSAR方法获得的是大范围相对形变量，而通过水准测量可以获得高精度的局部范围的绝对沉降信息。为了获取大范围的高精度地面沉降信息，则可以将二者进行数据融合。由于DInSAR测量的基准与水准测量基准不一致，所以需要进行基准的统一。通过数据分析，两者监测结果关系可采用下列函数表达[8]，令:

式中y为水准测量结果，x为雷达干涉测量结果。

基于2007.01～2009.10的水准测量结果和2007.01.17～2009.10.25的DInSAR结果(由于水准测量的周期较长，所以不能保证与DInSAR数据的时间段完全重合)，对水准与DInSAR结果进行了融合。首先按线性关系对两者进行了拟合分析，剔除偏差较大的点后，再用7个点分布均匀的数据对两者重新拟合，得到新的函数关系结果如图9所示。

图9 DInSAR数据与水准数据的线性关系图

用于拟合的点的校正后的DInSAR结果与水准数据的标准差为3.8 mm。则可利用水准与DInSAR之间的函数关系对DInSAR形变结果进行校正，可以得到整个区域的经过水准校正的地面沉降成果，从而获取天津实验区的从2007年01月17日～2009年10月25日的绝对形变量，结果如图10所示。5.3 实验区沉降结果分析

图10 绝对累积形变量

将经过水准数据修正之后的InSAR数据叠加到GoogleEarth地图上，各相干点目标沉降累积形变量的大小以不同的颜色分级表示，结果显示DInSAR数据点主要分布在天津市中心的范围，越靠近郊区的范围的相干点数据越为稀疏(主要为农田或鱼塘)，如图11所示。将沉降结果绘制成沉降等值线叠加到GoogleE-arth地图上，如图12所示。

从图5中沉降趋势分析，自2007年1月～2009年10月天津市主城区沉降速率较为稳定，累积沉降量约为4 cm左右。沉降速率从中心到郊区逐渐变大，说明主城区的地面沉降已经得到了有效的遏制，而郊区的沉降速率仍然偏大，且发现局部地区存在明显的沉降漏斗。

图11 天津实验区DInSAR相干点数据叠加到Google Earth图

图12 DInSAR沉降等值线叠加到Google Earth图

从图11和图12可以看出沉降漏斗主要分布在武清区、北辰区以及西青区。研究区域内最大(负值代表下沉)沉降漏斗中心位于西青区南河镇附近该沉降漏斗中心的累积沉降量约为26.6 cm和大寺镇附近该沉降漏斗中心的累积形变量达24.6 cm。武清区的沉降漏斗位于王庆坨镇附近，累积沉降量达18.5 cm。还有位于北辰区的西南方向的上河头镇附近的沉降漏斗，漏斗中心的沉降量约为21.6 cm。郊区沉降较大的原因是近几年天津工业快速向郊区转移和发展，郊区地下水抽取日益严重，导致地面沉降加速[9，10]。

6 结 语

本文针对天津地区地面沉降监测试验主要做了以下几项工作，首先是利用16景PALSAR单幅视影像进行常规二通差分干涉测量方法处理，采用精轨数据去除较大的轨道误差，利用SRTM DEM去除地形相位和地平相位的影响，利用最小费用流方进行相位解缠，从而得到差分干涉图。然后根据短基线集方法的原理，选取78个干涉对进行时间序列分析，采用拟合最佳轨道面的方法消除残余轨道误差的影响，并利用阻尼最小二乘算法求解DEM残余误差改正量以及所需要得到的平均线性形变速率。为了减弱大气相位的影响，这里做了空间低通滤波的处理，得到天津地区2007年1月～2009年10月的相对累积形变。最后为了获得绝对形变量，与7个水准点测量数据进行了融合，获得外部精度为3.8 mm的沉降结果，并确定天津地区沉降漏斗的分布及其绝对沉降量。

基于上述工作得到了以下几点认识:

(1)短基线集方法算法较PS方法简单，本文的特点在于在该算法基础上对轨道误差和大气相位的影响做了一定的分析和削弱；

(2)在利用阻尼最小二乘算法进行反演过程中，选择较为合适的阻尼因子是关键，一般选择在0.01～1之间；

(3)在水准数据与DInSAR数据融合过程中，由于水准点处没有高相干点数据，其水准点的DInSAR数据是通过插值得到，对最终结果有一定影响；

(4)利用水准数据可以将DInSAR的相对沉降信息转为绝对沉降信息，同时兼顾水准测量高精度的优点和DInSAR大面积、分辨率高的优点，能够在地面沉降监测的实际应用工作中发挥一定作用。

[1] 周春霞.星载SAR干涉测量技术及其在南极冰貌地形研究中的应用[D].武汉:武汉大学，2004.

[2] 王超，张红，刘智.星载合成孔径雷达干涉测量[M].北京:科学出版社，2002.

[3] Berardino，P.Fornaro，G.Lanari，R.et al.A new algorithm for surface deformationmonitoring based on smallbaseline differential interferograms[J].IEEE Transactions on Geoscience and Remote Sensing，Vol.40，No.11，2375～2383，2002.

[4] 董克刚，王威，于强等.天津地面沉降区土水比论述[J].水文地质工程地质，2008，19(3):54～59.

[5] 崔希璋，於宗俦，陶本藻等.广义测量平差(第二版)[M].武汉:武汉大学出版社，2009.

[6] 汤益先.基于相干目标的DInSAR方法及其地表沉降应用研究[D].北京:中国科学研究院，2006.

[7] 张红，王超，吴涛.基于相干目标的DInSAR方法研究[M].北京:科学出版社，2008.

[8] 张诗玉.干涉图大气效应校正方法研究及InSAR在天津地面沉降监测中的应用[D].武汉:武汉大学，2009.

[9] 范景辉，郭华东，郭小方等.基于相干目标的干涉图叠加方法监测天津地区地面沉降[J].遥感学报，2008，12 (1):111～118.

[10] 张诗玉，李淘，夏耶.基于InSAR技术的城市地面沉降灾害监测研究[J].武汉大学学报·信息科学版.2008，33 (8):850～853.

Study on Using PALSAR Data Extract Land Subsidence w ith Short Baseline Subset M ethod

Wang Cheng
(Chongqing Survey Insitute，Chongqing 400020，China)

The short baseline subset DInSAR technology can be very good to overcome the difficulties of poor coherence data，this paper is based on the analysis of shortbaseline subsetmethod toweak the atmospheric phase and orbital phase error，then gainmore accurate deformation phase.This paper collected 17 ALOSPALSAR single complex images since the January 17，2007 to October 2009，25which cover the Tianjin area，Intergration of DInSAR and leveling data to get external precision as the settlement results of 3.8mm.Finally，Tianjin ground subsidence funnel distribution and absolute variables are determined.

DInSAR；short baseline subset；groud subsidence；orbit error；atmospheric phase error

1672-8262(2013)05-81-06

P237

A

2013—03—18

王成(1988—)，女，工程师，主要从事航空摄影测量与遥感生产工作。

国家科技支撑计划(2011BAH12B07-03)