高分辨率TerraSAR-X时序差分干涉沉降监测及精度验证

于 冰， 谭青雪， 刘国祥， 刘福臻， 周志伟， 何智勇

(1.西南石油大学土木工程与测绘学院，成都 610500； 2.西南石油大学油气空间信息工程研究所，成都 610500； 3.西南石油大学油气藏地质及开发工程国家重点实验室，成都 610500； 4.中国科学院精密测量科学与技术创新研究院，大地测量与地球动力学国家重点实验室，武汉 430077; 5.西南交通大学地球科学与环境工程学院，成都 611756)

0 引言

地面沉降是指在一定的地面范围内发生的地面标高降低的现象，它通常是在自然和人为因素作用下引发的一种地质灾害[1-2]。地面沉降不仅会造成建筑物下沉、建筑基础变形、地面塌陷，还会造成水准点高程失真等，更严重的是地面沉降可能会间接地引发或加剧其他自然灾害的发生，如海平面上升、海水倒流、土地盐碱化、地面裂缝等[3-6]。天津市是环渤海地区的经济发展中心，是第一批沿海开放的城市，位于华北平原北部，地理位置较为特殊，地面沉降灾害尤为严重[6]，对天津市开展沉降监测和分析具有重要的现实意义。

传统的沉降监测方法，如精密水准测量、全球定位系统(global positioning system,GPS)测量等，人工成本较大，作业时易受天气影响，且监测速度慢，监测范围较小[7-8]。合成孔径雷达干涉(interferometric synthetic aperture Radar，InSAR)测量具有观测范围广，空间覆盖率和分辨率高、监测效率高等优势[9]。尤其是时序InSAR，如永久散射体干涉(permanent scatterer InSAR，PSInSAR)、小基线集(small baseline subsets，SBAS)干涉和相干点目标分析(interferometric point target analysis, IPTA)等，在地面沉降、滑坡形变监测等方面得到了较多的应用[10-14]。

IPTA是在PSInSAR技术基础上加以改进而来，相对于PSInSAR技术来说，IPTA技术提高了相干点的空间采样率和解算效率[10]。时序差分干涉测量(differential InSAR,DInSAR)虽有较多应用，但其精度验证仍是目前研究中的重要组成部分，现有文献大多会采用GPS或水准数据评估其精度。随着高时空分辨率SAR数据的投入使用，往往存在检核数据的时间采样率不足而无法对时序DInSAR所获取的时序沉降进行验证。比如，实际中GPS和水准观测周期比SAR卫星重访周期长，尤其是很多区域只能获取极少甚至只有一期观测数据，如何对时序沉降进行验证仍是一个需要解决的问题。

本文以天津市西青区及其周边区域为研究区域，以高分辨率的TerraSAR-X SAR影像数据为数据源，采用IPTA技术进行沉降信息的提取与分析，并提出一种最小二乘拟合的沉降时间序列验证方法，对监测结果进行精度验证，进一步检验高分辨率时序DInSAR应用于城市沉降监测的精度，并为缺少时序地面测量数据的情况下验证时序DInSAR结果提供一种有效的技术思路。

1 研究区概况及其数据源

1.1 研究区概况

天津市处于华北平原北部，向东邻近渤海，向北依靠燕山[15]。土地总面积约为11 966 km2，地质构造复杂。早期，天津市的地面沉降主要发生在城市地区，但是随着天津市经济与城市建设的外围发展，并加上主城区地下水开采控制，地面沉降逐渐从城市扩展到郊区[16]。研究区域主要包括西青区及北辰区部分区域，如图1所示，红色框为影像范围，底图为数字高程模型(digital elevation model,DEM)。

图1 研究区域示意图Fig.1 Study area

1.2 数据源

本文选取了2009年4月7日—2010年12月14日的34幅高分辨率TerraSAR-X降轨影像数据作为时序DInSAR处理的数据源，所有影像均为单视复数格式，波长为3.1 cm，图像成像空间分辨率约为3 m。在考虑影像获取时间(时间居中，时间最好在春冬季，且天气状况较好为最适宜)和基线长度(基线差较小为最好)的条件下，选取了2009年11月13日的影像为主影像，其他33幅影像为辅影像。确定一景影像为主影像后，要先对影像做预处理： 将其他33景辅影像配准重采样到主影像的几何空间。影像获取日期及辅影像相对于主影像的时空基线参数如表1所示。

表1 干涉对及时空基线Tab.1 Interferometric pairs and spatiotemporal baselines

2 IPTA基本原理及数据处理流程

IPTA的核心思路是从配准好的SAR影像集中选取具有稳定散射特性的相干点目标，提取其相位信息并以点为单位进行DInSAR处理，以线性形变模型和高程残差为基础建立二维线性回归模型，并迭代计算线性形变速率和高程残差，然后从扣除线性形变和高程残差的残余相位中分离非线性形变、大气延迟和轨道误差等，最后将线性形变和非线性形变相叠加获得时间序列[5,17]。

在本文中，基于所获取的34幅TerraSAR-X SAR影像，采用IPTA提供的时序振幅统计稳定性和光谱相关性2种方法选取相干点目标，并对二者进行合并，以此增加点目标的密度。然后，从34幅复数影像中提取相干点目标的原始相位进行干涉处理，本文采取公共主影像干涉组合模式，共获取33幅干涉相位图。此时，相干点上的干涉相位由参考椭球相位、地形相位、形变相位、大气延迟相位和失相干噪声等组成。接下来，基于所获取的外部DEM计算参考椭球相位和地形相位并从干涉相位中去除，得到33幅差分干涉图。此时，干涉图上第i个相干点的差分干涉相位可表达为：

φdiff,i=φdef,i+φtop_res,i+φatm,i+φnoi,i+φorb_res,i

(1)

式中：φdiff,i为相干点上的差分干涉相位；φdef,i为形变相位；φtop_res,i和φorb_res,i分别为残余地形及轨道误差相位；φatm,i和φnoi,i分别为大气延迟和噪声相位。

在得到相干点上的差分干涉相位后，IPTA基于邻域差分的思路建立相邻相干点之间的相位增量(差异)模型，该模型是关于相邻2点间线性形变速率增量和地形残差增量的线性模型，并包含其他噪声或者误差项，具体可表达为：

(2)

式中: Δφdiff,i-j为2相邻相干点间的二次差分相位； Δν为2相邻相干点间的线性形变速率增量；Βt为差分干涉图的时间基线；λ,R和θ分别为雷达波长、成像斜距和入射角； Δδh为2相干点间的高程改正增量；Bs为差分干涉图的空间基线；φorb_res,i-j,φatm,i-j和φnoi,i-j分别为相应的轨道误差增量、大气延迟增量和噪声增量。

在上面模型的基础上，IPTA方法数据处理的流程如图2所示。IPTA以相干点上所有的差分干涉相位为观测值，采用线性回归最优化求解线性形变速率增量和高程增量，并选定参考点求解所有相干点上相对于参考点的线性形变速率和高程改正。本文中IPTA解算的参考点经纬度为(E116.976 966 9°，N39.093 383 8°)。除此之外，IPTA提供迭代策略逐步精化这2个线性分量的解算结果。然后，从差分干涉相位中扣除线性形变和高程改正值，得到残余相位，主要由非线性形变、大气延迟、轨道误差和失相干噪声等分量组成。根据这些相位分量的时空频率特性不同，IPTA通过空间域和时间域滤波从残余相位中分离出非线性形变。将线性形变和非线性形变叠加即可得到形变时间序列，最后将形变时间序列转换到垂直向得到沉降时间序列。

图2 IPTA数据处理流程Fig.2 Flow chart of IPTA data processing

3 IPTA数据处理结果及精度验证

3.1 沉降速率精度验证

根据第2部分的数据处理方法和步骤，本文提取了1 231 705个相干点，获取了研究区域内的年沉降速率，如图3所示。

图3 相干点沉降速率及水准点分布Fig.3 Distribution of coherent point subsidence rate and leveling points

为了验证本文所提取沉降结果的可靠性与精度，根据水准实测沉降数据对所得结果进行检验(为保证后续结果分析的可靠性，先对监测结果进行验证)。现有研究也多是采用这种验证形式对雷达DInSAR结果进行检验[2,5,8,10,12]。在水准沉降监测中，采用二等精密水准测量，获取不同时期的高程数据，通过高程作差得到相应的沉降量。高程测量过程中每km高差中误差为1 mm，满足沉降监测规范要求，因此可以作为验证IPTA结果精度的基准数据。在此需要说明的是，由于时序DInSAR处理方法得到的沉降速率为相对沉降速率，而水准数据的沉降速率为绝对值，所以选用其中的一个水准点为校正点，来校正其他相干点，利用校正后的结果进行对比分析。水准点的空间分布如图3中黑色三角形所示，共计13个(1个基准点和12个检核点)。首先，从IPTA沉降结果中提取出与这13个水准点相近的相干点(如图3中红色圆点所示)，并记录沉降结果，用基准点对IPTA结果进行校正后，采用12个检核点进行验证。

在对二者进行对比分析时，采用差异均值和均方根误差[1,3](root mean square error,RMSE)2个指标，其公式分别为：

(3)

(4)

式中：k=12；Μv为水准点与相干点沉降速率差异的均值；vKi表示第i个水准点的沉降速率值；vLi表示第i个相干点的沉降速率值；N为点的数目；sv表示水准点与相干点沉降速率的RMSE。

本文选取12号点为校正点，其余12个水准点作为检核点。如图4所示为13个水准点和相干点沉降速率值的差异对比(12号点为校正点，差异为0)。图4中红色折线表示相干点的沉降速率，黑色折线表示水准点的沉降速率，蓝色折线表示水准点和相干点沉降速率的差异。从图4中可看出，这13个点的沉降速率值大致在-45 mm/a和-20 mm/a之间，12个检核点中，有8个点差值在3 mm/a内，差异最大值为-6.62 mm/a，最小值为0.06 mm/a，从图4中蓝色折线的走势也可看出相干点的沉降速率与水准点处的沉降速率偏差较小[2]。根据式(3)—(4)计算得到水准点和相干点沉降速率差异的均值和RMSE分别为： 1.02 mm/a和±3.15 mm/a。上述对比分析表明结合高分辨率TerraSAR-X数据的IPTA技术所提取的沉降速率与精密水准沉降监测结果具有较好的符合度。

图4 水准点与相干点沉降结果对比Fig.4 Comparison of subsidence between leveling points and cohevent points

3.2 基于最小二乘拟合的沉降时间序列精度验证

前已述及，在缺少时序水准数据(尤其是仅能获取1期水准数据)时，如何对时序DInSAR所获取的时序沉降进行验证是需要解决的问题。在本文中，为了进一步验证时序DInSAR处理方法的精度及可靠性，使用基于最小二乘的方法，根据沉降时间序列和观测时间段来拟合沉降速率，并分析拟合沉降速率与水准沉降速率的差异。对于每个相干点，可建立K个观测方程，即

(5)

式中：ΔtK表示第K景影像到第一景影像的时间间隔，d;v表示相干点上的沉降速率，每一个相干点对应于一个v;dK表示相干点在第K景影像所对应时间的沉降量;vdK为对应沉降量的改正数。将上述观测方程转换成矩阵形式为：

(6)

据此可解得每个相干点的沉降速率为：

v=(BΤPB)-1BTL

(7)

式中:B=[Δt1,Δt2,…,ΔtK]Τ;L=[d1,d2,…,dK]Τ;P为单位权阵。

在得到相干点上的拟合沉降速率后，采用和3.1节中相同的方案进行验证。相干点上拟合沉降速率和水准沉降速率的折线图如图5所示。

图5 拟合沉降速率验证结果Fig.5 Result of fitted subsidence rate verification

图5中蓝色折线表示水准沉降速率，绿色折线表示拟合的沉降速率，橙色折线表示水准沉降速率与拟合沉降速率差异值，差异绝对值的最大值为7.06 mm/a。根据式(3)—(4)计算得到水准沉降速率与拟合沉降速率差异的均值和RMSE分别为-3.4 mm/a和±3.25 mm/a，上述对比分析表明通过IPTA时序DInSAR处理方法所获取的时序沉降结果拟合所得的沉降速率与精密水准沉降监测结果具有较好的一致性。

4 研究区域沉降及成因分析

4.1 研究区域整体沉降分析

基于上述IPTA方法对2009年4月—2010年10月天津市西青区34景高分辨率TerraSAR-X SAR数据进行分析处理，获得了该段时间的年均沉降速率，结果如图6所示。在数据处理过程中共获取了1 231 705个相干点，这些相干点大多分布于建筑物、构筑物和道路等具有较高散射特性的地物上。从沉降速率结果中可看出，在监测期间研究区域地面沉降总体呈现不均匀性，其中与西青区邻近的北辰区沉降速率较大。

图6 沉降速率结果Fig.6 Result of subsidence rate

具体来说，由图6可看出研究区域不同地方的沉降情况，不同的颜色等级代表不同的沉降速率大小，从颜色分布情况来看，可将沉降速率分布大致分为3个级别。如图中暗红色、红色、橙色、黄色部分位于北辰区附近，沉降速率较大，在-128.41～-74.26 mm/a范围内，说明该地区沉降较为严重； 绿色、青色部分主要位于西青区附近，其沉降速率在-74.25～-33.65 mm/a范围内； 蓝色和深藏青色部分覆盖了北辰区部分地区、红桥区和南开区3个区域，沉降速率在-33.64～2.45 mm/a之间，该沉降区间覆盖的范围较大。由沉降结果及沉降分布情况可知，该研究区主要的沉降位于北辰区和西青区，结合相关资料可知该区域范围内存在较多的居民地、工业园区和高大建筑物聚集地等，工业生产用水、居民生活用水以及建筑物外部荷载等都会造成地面沉降。

为了对形变结果的空间分布规律进行进一步分析，分别在西青区和北辰区各提取了1个局部剖面AB和CD。剖面AB、剖面CD的沉降剖面图，如图7所示，图7中沉降速率是进行克里金插值后的值。如图7中所示，红色实线表示AB的沉降剖面，可看出在剖面上距离A点不同位置处的沉降速率，在距离A点约1.3 km处沉降速率最大。蓝色实线表示CD的沉降剖面，在距C点约2.1 km处沉降速率最大，是该区段的沉降漏斗中心。

图7 AB及CD的沉降剖面Fig.7 Subsidence profile of AB and CD

提取图6中点M，N，O，P，Q，R，S，T的沉降时间序列如图8所示，即根据累积沉降量与时间间隔的关系绘制的形变时间序列图。图8中所示，点M和点N处的累计沉降量最大，点M处的累积沉降量最大达到-170 mm，点O处最大沉降量达到-165 mm。点M和O所在区域为居民地和工业园区，年平均沉降速率在-100 mm/a左右。点N,P,Q,R,S,T沉降量在120 mm左右，所在区域为居民地、工业园区以及交通要道，年平均速率在-80 mm/a以上。

图8 沉降时间序列Fig.8 Subsidence history

4.2 典型沉降区及交通要道沉降分析

为了更好地对沉降区域的沉降情况及成因进行探讨和分析，本文从整个研究区域中提取了部分块状区域和线状地物的沉降速率，并获取了对应地物的可见光影像。图9(a)为北辰区的沉降速率，图9(b)为西青区的沉降速率。

(a) 北辰区沉降速率 (b) 西青区沉降速率

(c) 工业园区c (d) 居民地d (e) 居民地e (f) 工业园区f图9 2个地区的沉降速率及相应的用地类型Fig.9 Subsidence rate and land use category analysis of the two selected areas

在图9(a)中选取了c、d、e这3个沉降速率较大的区域，在图9(b)中选取了沉降速率较大的f区域，图9(c),(d),(e),(f)为其相应地区的用地类型。图9(d)和(e)所在土地用地为居民地，(c)和(f)用地类型为工业园区。结合图6和图9，可知图6中红色和黄色区域(图9中(c),(d),(e),(f))在实地上是工业园区、厂房以及成片的居民地等，图9(a)中c和d区域的沉降速率最大达到了-128.41 mm/a，说明该区域每年的沉降量最多。c和d分别对应于图9(c)和(d)，用地类型主要为工业区和居民地，该区域沉降如此明显，推测可能与工业生产和居民生活对地下水资源需求较大进而导致过度开采地下水相关。

研究区域内3条高速公路局部路段的沉降速率如图10所示，分别为津保高速、京台高速和荣乌高速公路。从图10中可知，津保高速公路的沉降速率最大，沉降速率大致在-92.3～-60.73 mm/a之间，高速公路附近有聚集的居民地、工厂和农田等，地下水资源需求较大，大量的地下水开采导致了严重的地面沉降。京台高速公路和荣乌高速公路的沉降速率比较接近，沉降速率大致在-51.69～-33.65 mm/a之间，高速公路附近也有农田和居民地，工厂较少，主要的地下水需求为灌溉和生活用水，开采量相对较小，所以沉降速率相对较小。

(a) 津保高速 (b) 京台高速 (c) 荣乌高速

图10 3条高速公路沉降速率Fig.10 Subsidence rate of three highways

5 结论

本文以天津市西青区及其周边地区为研究区域，以覆盖研究区的34幅X波段TerraSAR-X SAR影像为数据源，使用基于IPTA的时序差分干涉处理方法，提取了研究区域的地面沉降信息，并基于精密水准数据进行验证，给出一种基于最小二乘沉降速率拟合的沉降时间序列验证方法。最后基于验证后的沉降结果分析了研究区的区域沉降、主要道路沉降及其成因。研究表明：

1)经与水准数据对比，IPTA解算的沉降速率、基于最小二乘拟合的沉降速率与水准点的沉降速率均较为吻合，进一步验证了IPTA时序差分干涉处理方法的精度和可靠性。本文所提出的基于最小二乘拟合的时序差分干涉沉降时间序列验证方法可在缺少时序地面测量数据的情况下对监测所得时序结果进行验证，在一定程度上可以缓解缺少地面验证数据的问题，并完善时序差分干涉沉降监测结果验证的技术方案。

2)通过对整个研究区域沉降结果分析可知，在2009—2010年间，研究区域总体呈现出显著的不均匀性，区域内沉降最严重的区域在北辰区一带，最大年平均沉降速率达到了-128.41 mm/a，这与现有研究较为吻合。结合该研究区的人口分布特点、经济发展情况和地物分布情况，并对研究区域进行局部沉降分析可知，造成区域沉降不均匀的原因是不同地物类型对地下水的需求量不同，进而造成地下水开采程度有所差异，引发不均匀沉降。

综上所述，基于高分辨率TerraSAR-X数据的IPTA时序差分干涉处理方法在城市沉降速率及时间序列监测中具有较高的精度和可靠性，可以获取详细的地表沉降分布情况。同时，高分辨率SAR影像在监测道路等线状地物沉降方面表现出较好的效果。本文所提出的基于最小二乘拟合的时序沉降验证方法可有效缓解地面验证数据不足的问题。由于地面测量数据较难获取，本文未能采用更多的水准监测点进行结果验证，为了进一步完善时序沉降验证方案，后续研究将进一步对所提出的时序沉降验证方法进行更多的论证。此外，采用最新的SAR数据对该区域进行沉降监测并与历史沉降进行对比分析，探讨区域沉降的时空演变规律也是非常值得研究的问题。