基于SBAS-InSAR的淮安市地面沉降监测研究

夏斯雨，陆藩藩，晏王波

(1.江苏省天地图地理信息工程技术公司，江苏 南京 210000; 2.江苏省测绘研究所，江苏 南京 210000)

0 引 言

地面沉降又被称为地面下沉或地陷，是在自然条件和人为因素作用下，因地壳表层土压缩而导致区域性地面标高降低的一种环境地质现象[1-2]。

为了防治地面沉降带来的灾害，多种技术被用于监测地表形变，如传统的大地测量(GPS和精密水准)，合成孔径雷达干涉测量(InSAR)技术[3-4]等。常规的差分干涉测量(D-InSAR)技术在实际运用中容易受空间、时间失相干及大气延迟的制约，该方法主要用于大尺度地表突变的监测，而对于长时间的缓慢的地表形变过程监测受到很大的限制。近年来发展的永久散射体干涉测量(PS-InSAR)技术和短基线(SBAS-InSAR)[6]技术，在城市地面沉降监测领域逐步取代了传统的D-InSAR技术，成为利用SAR数据获取地表形变的主要手段，能有效提高干涉图的相干性，获取地表形变的沉降规律，大幅提高监测的应用场景和效益[10]。

淮安市地处江苏省北部中心地域，北接连云港市，东毗盐城市，南连扬州市和安徽省滁州市，西邻宿迁市。东西最大直线距离132 公里，南北最大直线距离150 公里，面积10 072平方公里。项目获取了研究区淮安市的数据，选用加拿大约25米分辨率的2012-2015年25景RADARSAT-2雷达影像作为数据源。选取2012年2月获取的影像为主影像，其余24景影像为辅影像进行精密配准，计算所有像对的垂直基线。DEM采用的是公开的30 m分辨率DEM数据。

1 SBAS-InSAR法求解沉降量原理

Berardino和Lanari等在2002年提出SBAS方法，即通过设定一定的时间、空间阈值自由组合干涉对，削弱时间、空间失相关的影响，该方法是InSAR监测中常见的方法之一。在SAR影像数量一定的情况下，可以形成更多的干涉像对，提高相位观测量，增加多余观测量，有利于提高最终解算结果的精度。

假设存在覆盖同一地区的N+1幅SAR影像图，通过设定空间基线与时间基线的阈值，组合干涉对，从干涉图中去除参考椭球面和地形的相位贡献，得到M幅差分干涉图， M满足下列条件：

(N+1)/2≤M≤(N*(N+1))/2

(1)

对于第j幅干涉图上的一点，其干涉相位为：

δΦj=Φ(tB)-Φ(tA)≈4π[d(tB)-d(tA)]/λ

(2)

式(2)中：λ为雷达波波长；Φ(tA)和Φ(tB)分别为tA和tB时刻对应的相位；d(tA)和d(tB)是相对于基准时刻t0的雷达视线向的地表形变量。

假设IE=[IE1,……,EM]和IS=[IS1,……,ISM]分别为干涉数据处理时根据时间顺序排列的主影像时间序列和从影像时间序列，且满足IEj>ISj。则所有的差分干涉相位可以组成如下观测方程：

δΦj=Φ(tIEj)-Φ(tISj)，j=1,……K

(3)

上式是由N个未知数的M个方程式组成的方程式，即：

δΦj=AΦ

(4)

式(4)中，A是一个M×N的系数矩阵，如果δΦj=Φ(t4)-Φ(t2)，δΦ2=Φ(t3)-Φ(t2)，则A表示如下：

当影像分多组时，上式是秩亏的，方程组没有唯一解。此时采用奇异值分解法对矩阵A进行分解，求解出方程的最小范数意义的最小二乘解。

采用奇异值分解方法直接对相位进行求解，得到的形变结果在时间上出现跳跃，即不连续，此时没有物理意义。鉴于此，将未知数转变为相邻影像获取时间内像元点沿LOS向的平均速率，即：

(5)

带入式δΦj=AΦ得：

(6)

Bv=δΦ

(7)

B是一个M×N的矩阵，对于第j行，位于主从影像之间的列B(j,i)=(ti-ti-1)，其他情况下B(j,i)=0。此时，对B进行奇异值分解，求取最小范数意义上的LOS方向的相位平均速率，当知道地表形变的先验信息时，可以利用形变相位在时间域上满足某些约束条件(即形变模型)来简化相位平均速率计算的复杂度。

考虑到高程误差对相位的影响，可以使用沉降模型来联合估算DEM误差，再根据各时段的沉降速率，对它们进行时间域上的积分，就可以得到各时间段的形变量。

2 淮安地区沉降速率计算与验证

2.1 淮安地区InSAR数据处理

SBAS方法的处理步骤主要包括(图1)：

①首先将监测区的SAR影像配准到同一的坐标空间下，并生成干涉图集；②通过滤波的方式去除由于失相干引起的噪声等；③利用辅助的DEM数据模拟相位去除地形和平地效应；④利用干涉图进行相位解缠并进行大气相位的去除工作；⑤得出淮安监测区沉降的平均位移速率和位置时间序列成果。

图1 SBAS方法处理步骤

为了保证解缠的效果，确定设定研究区内相干性阈值为0.35。在选取点的基础上，根据相干点构建Delaunay三角网，以三角网的边为单位进行形变建模，采用最小费用流量法进行差分干涉的解缠，将解缠相位校正到相干性较高的参考点上，求解得出2012—2015年内的灌河口周边地面沉降的平均形变速率(图2)。

图2 淮安市年平均沉降速率图(2012年-2015年)

通过在时间域上对各时间段的形变速率进行积分即得到各个时间点的累积形变量，并最终得到时间序列形变。可以看出：淮安市涟水县存在多个沉降漏斗且已经连接成片，最为严重的地区位于朱码镇，最大沉降速率达-42 mm/a，涟城镇、陈师镇也有较大沉降。

2.2 淮安地区CORS数据验证

为验证SBAS方法得到的结果，搜集了12个灌河口周边2012-2015年逐年观测的CORS数据，比较CORS实测点与最邻近点InSAR点目标测量得到的形变值，根据两者的差异评估InSAR地表沉降监测结果的精度，具体的CORS点分布图如图3所示。

由于获取InSAR监测点覆盖淮安市，CORS站点也分布在研究区范围内，所以在验证时，认为CORS站点与最邻近InSAR点的地质条件及受人类活动影响相同，故不做额外的误差分析。

将实测结果与SBAS方法所得监测结果进行比较(表1)：

图3 CORS点分布图

表1 CORS点与InSAR监测结果对比mm/a

二者插值最大值为1.85 mm/a、最小值为-1.34 mm/a。通过SBAS-InSAR技术获取的地表形变量与已有的CORS点检测的沉降值进行比较，利用均方根误差(Root-Mean-Square Error，RMSE)作为评定精度的指标，其公式如下：

(8)

其中，N为站点的个数；ΔHi为CORS点监测到第i个点的沉降结果；Δhi为InSAR检测到第i个CORS点的沉降结果。计算得到RMSE为3.04 mm/a。

3 结 语

本文利用短基线方法(SBAS)监测了淮安市2012-2015年的地表形变，获取了时间形变的特征、累积沉降量和年平均沉降速率。研究区内年平均沉降速率最大达-42 mm/a。

(1)从时间角度来看，在监测周期内，2012-2015年，沉降漏斗逐渐形成且有成区连片趋势，主要原因是各化工园区企业规模数量不断增加，生产效能提高，导致地面沉降，故地面沉降逐年累加。

(2)从空间角度来看，沉降区域相对较为集中，主要分布在涟水县朱码镇、涟城镇、陈师镇，大规模的人类生活生产活动导致了地面的沉降塌陷。

(3)从实验方法来看，在以平原地貌为主的灌河口周边利用短基线方法(SBAS)进行地面沉降的监测依然有效，这也为大规模进行全省地面沉降监测奠定了基础，有益于推进全省减灾防灾工作。