利用Sentinel-1 SAR数据及SBAS技术的大区域地表形变监测

张金盈，崔 靓，刘增珉，王新田，林 琳，徐凤玲

(1.山东省国土测绘院，山东 济南 250013; 2.山东省水利科学研究院，山东 济南 250014; 3.山东华峰地理信息科技有限公司，山东 济南 250101)

星载SAR系统对全球地表形变监测具有重要作用。近几年广泛应用于城市、矿区、滑坡监测的SAR数据有如下卫星数据：2014年日本发射的ALOS卫星，幅宽最大为70 km；2007年德国发射的TerraSAR卫星，提高了干涉数据的相干性，幅宽最大为30 km；2007年加拿大发射的radarsat卫星，提供高分辨率全极化影像，幅宽最大为150 km。这几种数据的共性为覆盖范围较小，在大范围监测中面临数据拼接等复杂问题。2014年欧空局发射Sentinel-1卫星，该卫星具有4种成像模式，最大幅宽达400 km，满足大范围的数据处理。同时欧空局承诺为用户免费提供哨兵家族数据以供科学研究，大大促进了时序雷达干涉测量(SAR interferometry,InSAR)技术的发展和在大区域的应用[1]。

相比永久散射体合成孔径雷达干涉测量(PS-InSAR技术)，小基线集(small baseline subsets InSAR，SBAS InSAR)技术是在限制时空基线阈值的前提下，对获取的SAR影像集进行自由干涉组合，从而得到一系列短基线差分干涉相位图集。这些差分干涉图能够较好地克服时空失相关的影响[2]，可以显著增加相干影像对的数目，一定程度上解决了大范围时间、空间失相关和大气效应对InSAR的影响。

黄河三角洲区域为黄河入海口，地势低平，地质成因晚，拥有油田、天然气等丰富矿产资源，由于地质原因和人类生产活动，存在地表沉降现象[3]。文献[4]利用SBAS技术，采用39景ERS1/2数据对面积为5400 km2的东营地区1992年至2000年间的地表情况进行了监测，监测出胜利油田中心开采区最大沉降速率为-33.2 mm/a。文献[5]利用28景ALOS数据，监测东营市胜利油田2007—2011年间的地表沉降，分区域分析了地表沉降及成因。大面积区域性地面沉降调查存在数据相干性问题、低相干区沉降信息难提取、大气波动影响等问题[6]。本文利用Sentinel-1数据，采用SBAS方法对覆盖山东近1/4范围的6.3×104km2的黄河三角洲区域进行沉降监测研究，并进行了精度评价和成因分析，验证了利用Sentinel-1数据的SABS方法在大区域地表沉降监测的适用性，总结了黄三角区域在两年半时间内的地表沉降时空变化的主要特点。

1 Sentinel-1数据及SBAS技术介绍

1.1 Sentinel-1数据

Sentinel-1卫星是欧洲航天局哥白尼计划(GMES)中的地球观测卫星，由两颗卫星组成，载有C波段合成孔径雷达，可提供连续图像(白天、夜晚和各种天气)。该卫星具有超高的辐射分辨率，能有效提升雷达图像参数反演精度，能实现全球80%的覆盖，单个卫星在中国境内12 d的重访周期且实现80%覆盖。两个卫星星座可缩短至6 d。其干涉宽测绘带模式(interferometric wide swath，IW)带宽达250 km，单景影像可覆盖面积为4.5×104m2，空间分辨率为5 m×20 m，两景影像即可实现黄河三角洲区域全覆盖。

1.2 SBAS技术介绍

SBAS方法通过增加干涉数据的采样，获得较高精度的沉降信息和高程[7]，即在一定范围内设置空间基线和时间基线，将SAR影像组合成若干个小基线集，分别选择一个主影像进行配准，利用差分干涉测量的方法获得每个集合中的差分干涉图[9]，根据相关系数提取相干性较高的点，应用奇异值分解法(singular value decomposition，SVD)将多个小基线集联合起来求解，最终获得时间序列沉降图及沉降速率[8]。

假设不考虑由于地形、大气延迟、噪声等造成的相位误差，可建立如下模型

Aφ=δφ

(1)

式中，A为M×N稀疏系数矩阵，每一行对应一幅差分干涉图，每一列对应不同采集日期的同一幅影像，主影像的系数为1，辅影像的系数为-1，其他为0；φ为N+1个成像时刻的雷达图像上未知沉降相位点组成的矩阵；δφ为M幅差分干涉图中观测相位值矩阵；由于SAR影像由多个小基线集组成，M大于N，系数矩阵A秩亏。对矩阵A进行奇异值分解(SVD)

A=USVT

(2)

式中，U为M×N正交矩阵，由AAT的特征向量ui组成；V为N×M的正交矩阵，由ATA的特征向量vi组成；S为一个M×N的对角矩阵，对角线元素是AAT的前R个非零特征值，后面的M-R个特征值为0，A+为A的伪逆矩阵

(3)

经过奇异值分解，可以得到沉降相位φ的最小范数意义上的最小二乘解

(4)

解算出线性沉降速率和高程误差，去除线性部分，得到残余相位。对残余相位在空间上进行低通滤波，在时间上进行高通滤波，分离出大气相位，就可以得到非线性沉降。将非线性沉降叠加到线性沉降中，就可以求出总的地表沉降信息[10-11]。

2 试验处理流程及结果

2.1 区域概况

研究区域是黄河携带泥沙在入海口淤积形成的冲积平原，地势低平。近年来，由于人类日益频繁的经济活动和黄河带来的巨量沉积物的自重固结影响，黄河三角洲区域出现明显的地表沉降。根据已有研究，1992—2000年该地区最大沉降速率达-33.2 mm/a[4]；2015年5月—2016年4月期间，该地区最大年平均沉降速率为224.6 mm/a[12]。

2.2 SBAS技术处理流程

本文研究使用Sentinel-1A IW模式升轨影像进行监测，对相邻2景影像进行拼接，共采用32期64景Sentinel-1A影像，时间跨度为2015年7月—2017年8月，极化方式为VV极化，C波段。处理范围如图1所示，该测区覆盖东营、潍坊、淄博、滨州等地，陆地面积约为6.3×104km2。轨道数据使用与Sentinel-1A数据间隔21 d分发的精密轨道数据，为了精细化模拟地形相位，采用30 m分辨率的SRTM DEM数据。

验证数据采用均匀覆盖山东全省的SDCORS基准网和2015年12月—2016年12月逐年观测的GPS数据点。

根据SBAS网络构建原则，通过设定不同的时间阈值及空间阈值，优化影像组合方式，以提高干涉图的相干性。根据干涉的效果和垂直基线大小，选取基线小于100 m时间基线小于120 d的194对干涉图，如图2所示。采用分辨率为30 m的SRTM DEM生成模拟地形相位如图3所示，然后对干涉图进行差分干涉处理[13-18]。

图2 SAR影像组合连接

图3 SRTM DEM

首先提取相干系数大于0.4的点为永久散射体点，然后利用Delaunay不规则三角网建立候选点之间的连接关系，对永久散射体点进行解缠和地理编码，最后利用矩阵奇异值分解的方法估算沉降速率和高程系数，对结果进行时间域的高通滤波，对每幅干涉图进行空间域的低通滤波，从而获取大气延迟相位，将其从原始差分干涉图中去除，得到去除了大气影响的沉降结果。黄三角区域监测时段内形变速率如图4所示。

图4 SBAS方法监测黄河三角洲地区形变速率

3 可靠性验证与分析

3.1 利用CORS站点的精度验证

为验证该方法反演地表沉降的精度，利用收集到的监测区域2016年3月和2017年3月两期CORS数据点沉降观测数据进行验证。

根据验证点的地理坐标，采用最近点法搜索站点5个像元(150 m)内的相干目标点，共搜索到符合条件的17组CORS验证数据如图5所示。17组验证数据见表1，分别对每组的相干点沉降观测值和验证点沉降测量值进行时间序列的曲线拟合，部分验证点时间序列沉降图如图6所示。SBAS技术观测的相干点变化趋势与CORS站及GPS点解算的沉降变化趋势基本一致，解算中误差为7.8 mm。

图5 17组CORS站位置

图6 部分SBAS-InSAR观测值和CORS站观测值曲线拟合

表1 17组CORS站验证数据组中误差 mm

3.2 形变特征分析

图7为研究区域提取出的沉降漏斗，地面沉降严重区域主要分布在黄河下游(A、B区域)，平均沉降速率在50～152 mm/a。A区域主要呈现6个沉降明显的漏斗，面积为1500 km2，最大沉降速率达152 mm/a。实地调查发现该区域分布石油化工、钢板厂等多家高耗水企业。同时通过文献收集[19]，发现该沉降区域与地下水降落漏斗区域基本重合。

图7 黄三角地区平均沉降速率点分布

B区域分布在黄三角沿海区域，存在3个典型沉降漏斗，面积为500 km2，最大沉降速率达111.72 mm/a。通过现场调查B地区，主要分布盐场和溴素等工厂，开采地下卤水造成一定程度沉降，同时该区域地层沉积较晚，地层岩性以粉土为主，因此土的先期固结较低，在自重压力下长期固结压密，在一定程度上也会造成地面沉降。

4 结 语

本文重点介绍了采用SBAS技术，利用Sentinel-1A影像获取6.3×104km2黄河三角洲区域地表沉降的方法。与CORS站点测量结果相比，中误差为6.9 mm，表明利用Sentinel-1数据结合SBAS技术可以进行大范围地表沉降监测。研究结果发现，黄河入海口的平原区域存在多个降落漏斗，面积约2000 km2。最大沉降达152 mm/a。根据文献调查及地质资料分析，分析地表下沉主要是地下水及沿海卤水长期超量开采且得不到补偿导致。