基于SBAS-InSAR技术和Sentinel-1A的城市地表沉降监测

王 磊，连增增,刘 杰

（1.河南省有色金属地质矿产局第一地质大队，河南 郑州 450000；2.河南理工大学测绘与国土信息工程学院，河南 焦作 454000）

地面沉降是由多种因素引起的地面标高缓慢降低的地质现象，是一种主要由人类活动引起的环境地质问题。随着城市化进程的加快，地面沉降这一地质问题日益突出，城市基础设施遭到破坏，城市规划布局、土地有效利用及地下空间开发受到严重影响[1-2]。传统的地面沉降监测手段有GPS测量、水准测量等，虽然精度高，但具有观测周期长、人力成本高及空间分辨率低等缺点[3]。相对于传统大地测量手段，InSAR技术具有全天候、全天时、成本低、覆盖范围大、空间分辨率高等优势[4-5]。D-InSAR技术是基于InSAR干涉相位获取地表形变信息的一种新方法，可以探测亚厘米级形变，但D-InSAR技术易受时空失相关、大气延迟等因素影响[6]，进而导致形变结果精度的降低，甚至导致错误的形变提取结果，较难完成高精度、长时间跨度的地表形变监测[7]。

永久散射体合成孔径雷达干涉测量（persistent scatterer interferometric synthetic apertureradar，PS-In⁃SAR）技术[8]能够通过对地面高相干目标点进行时序差分干涉相位分析，并剔除低频大气的影响，较好地克服失相干与大气延迟影响问题[9]，但要得到比较可靠的监测结果需要较多的SAR影像数据，一般不少于25景[10-11]，一定程度上限制了PS技术的发展。小基线集合成孔径雷达干涉测量（small baseline subset inter⁃ferometric synthetic aperture radar，SBAS-InSAR）技术[12]将SAR影像数据组成若干个子集，采用最小二乘法求解子集的形变时间序列，同时利用奇异值分解法（singularvalue decomposition，SVD）将多个小基线集联合求解，获取整个时间跨度的形变序列[13]，相对于PS-InSAR技术，SBAS-InSAR需要的SAR影像数目较少且获取非线性形变信息的能力较强[10,14]。现有成果[15-19]已利用SBAS-InSAR技术对城市地表沉降进行了监测，均取得了理想成效。但是，由于Sentinel-1A影像数据的新颖性，基于SBAS-InSAR技术和Sen⁃tinel-1A影像的城市地表沉降监测方面的研究尚显不足。

本文在对SBAS-InSAR技术进行分析的基础上，以位于顺义-朝阳-通州沉降带的三河市燕郊镇为研究区，对覆盖该研究区2017—2018年的20景Senti⁃nel-1A数据进行处理，获取该区域地表年平均形变速率及时序形变信息。同时，本文还监测了PS-InSAR技术下的研究区地表沉降信息，并进一步对比分析了2种监测技术之间的差异性。最后，本文探讨分析了引起研究区地表沉降的主要成因。

1 SBAS-InSAR技术原理

设有按时间序列t0，…，ti，…，tN获取到N+1幅单视复数SAR影像，将它们以任意影像为主影像进行配准；设定垂直基线阈值，将垂直基线小于该阈值的SAR影像归为一组，设共分L组，对分组后的影像做差分干涉处理，可得到M幅差分干涉图，假设N为奇数，则差分干涉图的个数M可以表示为：

以t0为初始时刻，任意时刻ti(i=1，…，N)相对于t0时刻的差分相位φ(ti)为未知参数，差分干涉相位δφ(tk)(k=1，…，M)为观测量。对于第k(k=1，…，i，…，M)幅差分干涉图的任意像元(x，r)有

式中，λ为雷达波长；d(tA，x，r)和d(tB，x，r)分别为tA和tB时刻像元(x，r)相对于初始时刻t0的雷达LOS像地表形变，即有d(t0，x，r)=0。为获取研究区的时间形变序列，需要精确预估出地形相位误差分量、大气延迟相位误差分量以及噪声相位分量，并将这3个分量从干涉相位δφk(x，r)中去除。

根据式（2），在去除各项误差分量后，M幅差分干涉图可以得到M个方程，用矩阵表示为：

式中，A为M×N矩阵，由-1、0、1组成；φΤ=[φ(t1)，…，φ(tN)]为每一景SAR影像中高相干点对应的相位值所组成的向量；δφΤ=[δφ1，…，δφM]为各差分干涉图对应的解缠相位值所组成的向量。

为求解研究区域各高相干点的形变速率，可用两景SAR影像间的平均相位速率代替相位值，则式（3）变为：

式中，B为M×N的系数矩阵；vT可以表示为：

当系数矩阵B为满秩（即M≥N）时，可用最小二乘法则求解出形变速率；当M＜N时，矩阵B出现秩亏，可利用奇异值分解法（SVD）获取研究区域的形变速率[16]。获取形变速率后，依据研究区SAR影像时间区间求得相应时间段内的形变量[15]。

2 研究区概况与数据来源

2.1 研究区概况

燕郊镇位于河北省三河市，如图1所示，西临北京市通州区，东临天津市，中心地理坐标为39°56'N、116°48'E，覆盖面积约324 km2。区内地势平坦，位于潮白河畔，东南部为缓向渤海倾斜的平原。气候为暖温带半湿润半干旱季风气候，夏季高温多雨，冬季寒冷干燥，春、秋短促。现有资料显示，研究区属于京津冀地区沉降较为严重的区域，且与北京五环外的顺义-朝阳-通州沉降带相连接。目前，由于研究区人口密度大，导致本研究区相对于其他区域的沉降较为显著[20]。

图1 研究区概况图

2.2 数据来源

本文选取时间跨度为2017-03-21~2018-10-06的20景覆盖研究区域的Sentinel-1A1级影像数据作为研究数据，采用波长为5.6 cm的C波段，轨道方向为升轨，方位向与距离向分辨率分别为13.985 m与9.318 m，影像中心入射角约为43.8°，极化方式为VV极化。采用美国宇航局提供的30 m分辨率SRTM1 DEM数据去除地形相位影响，同时利用欧洲空间局发布的POD精密轨道数据去除轨道误差。

3 数据处理与结果分析

3.1 数据处理

本次研究采用SAR-scape系统进行数据处理，具体数据处理过程如下：

1）数据预处理：将Sentinel-1A原始影像裁剪至本次研究区域大小，如图1所示，对裁剪后的影像进行5:1（方位向×距离向）多视处理。

2）干涉处理：为避免空间失相干，设定时间基线阈值为180 d、空间基线阈值为150 m进行差分干涉处理，由系统自动选出日期为2017-11-16获取的影像为公共主影像，这样可以保证有足够多的像对，在对像对编辑时可以将相干性低的像对进行去除。将其余所有影像进行配准并重采样至公共主影像，并选用3D解缠，共生成92对小基线差分干涉像对，像对组合方式如图2所示，图2中黄颜色点代表公共主影像，线段长度代表干涉基线长度，统计组合后干涉相对最大时间基线为180 d，最大空间基线为110 m。

图2 时空基线组合

3）轨道精炼及重去平：对恒定相位进行估算去除，在此步骤进行之前需要对上一步生成的干涉结果进行质量检查，若在干涉图、解缠结果图以及相干性图中出现质量较差（失相干现象严重）的结果，则考虑将该像对做删除处理。利用GCP控制点对所有像对进行重去平，本次研究共选取了30个GCP点。

4）SBAS第一次反演：设置相干系数阈值为0.35，估算形变速率和残余地形，同时进行二次解缠对干涉图进行优化。

5）SBAS第二次反演：分别设置大气高通滤波、大气低通滤波2个参数进行大气滤波，估算并去除大气相位，计算时间序列位移信息。

6）地理编码：依据振幅与相位的稳定性筛选研究区内的高相干点，通过奇异值分解（SVD）法求解高相干目标点的沉降速率，将第二次反演结果从SAR坐标系转换到地理坐标系，最后得到研究区的LOS向形变速率及时序形变信息。

3.2 结果分析

1）研究区形变速率。由SBAS-InSAR技术获取到研究区2017-03-21~2018-10-06的地表年平均沉降速率信息，如图3所示，研究区下沉速率较大区域主要集中在西部与中部（图中所示多边形A区和B区）。整个研究区的年平均沉降速率范围达到-62~26 mm/a，不均匀沉降较为明显，平均沉降速率超过20 mm/a，同时存在的沉降漏斗平均沉降速率超过36 mm/a。

图3 SBAS技术研究区LOS向平均沉降速率

2）研究区时序形变。以2017-03-21为起始时间（参考形变量为零）获得研究区的时序形变信息，如图4所示，至2017-09-17研究区中部（图中深蓝色区域）已经发生了较大沉降，最大沉降达到34 mm；至2018-03-28研究区西部发生较大沉降，最大沉降达到45 mm，研究区沉降范围基本稳定；至2018-10-06，研究区中部与西部最大沉降均达到74 mm。

图4 SBAS技术研究区LOS向时序沉降变化

为验证SBAS-InSAR监测结果的准确性，同时采用PS-InSAR技术进行地表形变监测，监测结果如图5所示。由图5可知，沉降同样主要发生在研究区的中部与西部（图中所示多边形A区和B区），最大沉降速率达到24 mm/a。

图5 PS技术研究区LOS向平均沉降速率

为分析2种技术在地表形变监测结果的差异性，在研究区的集中沉降区域提取同名特征点，其中A5特征点位于桥面，其余所有特征点均位于建筑物表面，点位分布如图6所示。通过放大图6可以发现，SBAS技术在监测地表形变细节上更具优势，在A2、A3、B1点所在的特征点区域均形成明显的下沉盆地，且下沉范围连续集中；而在PS监测结果中，A2、A3、B1点所处区域没有监测到下沉盆地形态。此外，整体而言，相较于SBAS技术的监测结果，通过PS技术所监测到的高相干点出现了大面积缺失。

图6 同名特征点分布图

为定量分析SBAS与PS的监测结果，分别从研究区的A区和B区提取特征点的时序形变信息，提取结果如图7所示。从图7可以看出，SBAS获取的沉降值整体大于PS获取的沉降值。尽管两组曲线的走向是一致的，但是利用SBAS技术监测的特征点时序曲线较为平滑，而PS技术监测的特征点时序曲线出现较大幅度波动，且在A3、A6点无法获取有效监测值。

图7 SBAS与PS技术结果特征点时序累计沉降量

通过计算各特征点至2018-10-06的累计沉降值可知（表1），位于桥面的A5特征点累计沉降值为18 mm，两者沉降值差为零；A2与B1特征点沉降值差较大，原因是这2个特征点在SBAS结果中均出现在下沉较为集中区域，而在PS结果中此区域无高相干点被探测，且其周边也出现大面积缺失值；对其余各特征点沉降值差求取均方根误差为5.6 mm，说明SBAS-InSAR方法在监测地表形变是可靠的。

表1 各同名特征点最大沉降值差/mm

通过研究结果的交叉对比可知，SBAS技术和PS技术获取的沉降结果在监测细节上存在轻微差异，但是总体来讲，2种方法获取的形变区域具有较高的一致性，这也证明了该研究区域在2017年3月至2018年10月期间存在地面沉降问题，且平均沉降速率均达到20 mm/a。导致2种方法监测下沉信息有所差异的原因是由于SBAS是非线性沉降模型，PS是线性沉降模型，两者所描述的沉降过程不完全一样，存在局部一致或局部不一致的现象。本次研究使用20景Senti⁃nel-1A数据，数据量相对较少，在同等数据量的情况下，SBAS可以形成高密度的干涉像对组合，而PS方法过于依赖于单一公共主影像的质量，且在识别高相干点时对影像数据量有一定要求，因此在数据量有限时，SBAS-InSAR方法监测地表形变可以获得更为可靠的结果[21-22]。

3.3 形变成因分析

从SBAS和PS技术提取的结果可以看出，沉降区域位于研究区中部与西部，造成沉降的原因有：①该区域主要分布于潮白河、温榆河和泃河流域的冲积、洪积扇平原上，多条河流横跨研究区，属于冲积扇的中下部，土质主要由细沙黏土构成[22]，同时地下水过度开采会导致地面沉降；②近几年北京地区的城市化向东扩张依旧呈现逐年增加趋势，本研究区与北京市沉降区逐渐连成一片[20]；随着一些重工业的东移，地表构筑物的建设以及地下工程建设打破了原有地下岩土体应力平衡，且这些建设工程在空间分布不均匀，造成的区域性集中沉降分布不均匀，形变量级存在明显差异，造成该区域出现较为显著的地表沉降现象，形成明显的下沉盆地。

4 结语

本文利用SBAS-InSAR技术对20景Sentinel-1A数据进行处理，得到了研究区2017~2018年的沉降速率及时序形变信息，发现下沉速率较大区域主要集中在研究区中部与西部，平均沉降速率超过20 mm/a，部分集中沉降区域平均沉降速率超过36 mm/a，累计沉降达到74 mm。同时利用PS-InSAR做结果交叉验证，发现SBAS方法和PS方法获取的形变区域具有较高的一致性，但SBAS方法在监测地表细节更占优势，在数据量有限时，SBAS-InSAR监测地表形变可以获得更为可靠的结果。研究区沉降主要受到其地质构造、地下水开采及北京地区的城市化向东扩张有关。