邓晓景，曲国庆，张建霞，席换，王晖

(山东理工大学 建筑工程学院, 山东 淄博 255049)

地面沉降是指由于自然因素和人类活动因素共同作用下引起地面土层压缩，从而发生的地质灾害[1]。传统地面沉降监测方法，如水准测量、GPS测量等，受设备、外部环境等因素影响，存在监测时间长、监测范围小、维护成本高等缺点，难以实现连续、准确的地面沉降监测。随着雷达卫星监测技术的不断发展，合成孔径雷达差分干涉测量技术(Differential Interferometric SAR, D-InSAR)已经成为一种全天时、全天候、高精度监测地面沉降的有效方法[2]。但D-InSAR技术易受空间失相干、时间失相干、大气延迟等因素影响，使其监测精度降低，无法进行长时间序列的地面沉降监测[3]。永久散射体干涉测量技术(Permanent Scatterer InSAR,PS-InSAR)通过获取高相干性目标点，采用时空滤波的方法，能有效的削弱基线误差、大气延迟误差等因素的影响，得到高精度时序地面沉降监测结果[4]。

在利用PS-InSAR技术进行大区域地面沉降监测时，由于单一轨道覆盖范围有限，无法满足研究区需要，往往需要将升轨和降轨沉降监测结果进行拼接。由于在同步观测时间下，相同研究区域内多种观测模式下获取的监测结果是对相同地面点的多角度采样，这为升降轨两组解算结果的精度互检提供了条件[5]。由于监测结果中垂直方向分量受位移影响较小，将雷达视线向形变量转换为垂直向形变量后再进行精度互检，这样可以提高监测结果内部互检的准确性。

升降轨PS-InSAR数据融合的基本条件是提取同名相干目标的形变量。然而，受雷达采样方式差异的影响，升降轨数据中PS点的分布位置不相同，导致在升降轨解算结果中，同名相干目标点数量极为有限，严重影响升降轨解算结果融合的可靠性[6]。针对上述问题，本文提出了利用抗差最小二乘拟合法对升降轨PS-InSAR解算结果进行精度内部互检及数据融合的方法。对升降轨数据进行融合，可以加密单轨条件下的观测序列，能够更加丰富地显示研究区地面沉降监测时序变化规律。

1 研究区概况及数据选取

研究区为山东省东营市，其大部分位于黄河三角洲冲积平原上，地势平坦，地形以平原为主[7]。该地区海洋资源丰富，沿海地区分布大量晒盐和养殖场，且石油、天然气等自然资源丰富，其地层覆盖主要为第四系，岩性由粉砂、细砂、黄色黏土和沙砾层组成[8]。东营市近年来由于矿产资源开采以及人类活动等因素影响，开始出现地面沉降现象。

图1 研究区概况及数据范围Fig.1 Overview of study area and data range

哨兵1号(Sentinel-1)卫星是欧空局哥白尼计划(GMES)中的地球观测卫星，发射于2014年，由两颗卫星组成，载有C波段合成孔径雷达，本研究选择覆盖东营市的Sentinel_1A雷达卫星数据，该数据采样模式为IW宽幅干涉模式，在该模式下采用TOPS (terrain observation with progressive scans)扫描方式，可以同时获取3个子条带，其幅宽为250 km,重访周期为12 d，空间分辨率为5 m×20 m[9]。实验数据选取分别为降轨数据(Track-76)左侧子条带和升轨数据(Track-69)中间子条带，数据参数见表1。参考DEM数据为美国宇航局提供的SRTM(shuttle radar topography mission)数据，其分辨率为90 m。

表1 Sentinel-1数据参数Tab.1 Sentinel-1 data parameters

2 研究方法

2.1 PS-InSAR技术

PS-InSAR技术将数据处理分为数据预处理、差分干涉计算和形变量计算三个步骤。首先，选取覆盖研究区的N+1幅SAR影像，根据基线参数和多普勒效应，选取一幅SAR影像作为公共主影像，其余影像为辅影像，分别与公共主影像配准进行干涉处理，生成N幅干涉图；通过精密轨道文件去除平地效应，同时利用参考DEM数据去除地形相位，对N幅干涉图进行差分处理，生成差分干涉图[10]。其中，任意差分干涉图像元的相位Δφ可以表示为：

Δφ=φdef+φtopo+φatm+φnoise

(1)

式中：φdef为由于地形位移而产生的形变相位,φtopo为外部DEM引起的残差地形相位,φatm为大气因素引起的延迟相位,φnoise为其他因素引起的噪声相位。

然后，通过设定合理振幅离差阈值，提取PS候选点，采用Delaunay三角网法构建PS网络；根据建立相邻PS点相位差分模型，获得地表形变信息；利用空间搜索法，将PS点中线性形变和高程误差去除[11]；通过时间域上高通滤波，空间域上低通滤波，获取并分离非线性形变相位和大气延迟相位，最终得到研究区时序地面形变量[12]。

2.2 升降轨数据融合简化模型

PS-InSAR技术获取的地表形变量可分为东西方向和垂直方向的地表形变量，对于同一研究区的升降轨数据，PS-InSAR技术获取的观测结果东西方向上水平形变量呈现相反的特点，而垂直方向上垂直形变量基本相同[13]。雷达视线向(line of sight,LOS)累计形变量dLOS为：

dLOS=dNsinφsinθ-dEcosφsinθ+dVcosθ

(2)

式中:dN为南北方向的形变分量,dE为东西方向的形变分量,dV为垂直方向的形变分量,φ为方位角,θ为入射角。由于南北方向的水平位移对垂直向和东西向的影响非常小[5]，所以雷达视线向累计形变量简化为：

dLOS≈dEsinθ+dVcosθ

(3)

地表形变量以垂直方向的垂直形变量为主，东西方向的水平形变很小，即dE≪dLOS,因此，垂直方向的垂直形变量可近似表示为：

dV≈dLOS/cosθ

(4)

针对升降轨数据融合中同名目标点数量稀少问题，传统方法采用选取一条轨道解算结果作为主轨道，另一条轨道的解算结果为辅轨道，采用插值方法提取升降轨道相同位置的形变量，根据公式(5)进行升降轨道解算结果的数据融合[14]。

(5)

在实际数据处理过程中，由于受到升降轨数据入射角以及采样方式的影响，即使相同位置的垂直形变量偏差也不尽相同，只用主副轨道垂直形变量的平均偏差进行整体修正，会严重影响修正后的辅轨道垂直形变量的可靠性。因此，通过对辅轨道进行插值处理，获得主副轨道相同位置的整体偏差，再利用抗差最小二乘方法对整体偏差进行曲面拟合，用偏差拟合曲面矫正辅轨道，从而进行整体形变场的解算。

2.3 抗差最小二乘曲面拟合

根据研究区内升降轨解算结果的覆盖情况，将降轨作为主轨道，升轨为辅轨道，以升降轨共同覆盖区域为研究对象，对升轨解算结果进行插值处理，获得降轨相干目标点同名位置的偏差，采用二次曲面方程进行偏差拟合,设二次曲面函数为：

Z=a0+a1x+a2y+a3xy+a4x2+a5y2

(6)

式中：Z为已知降轨相干目标点偏差值；x，y为已知降轨相干目标点的坐标值；a1，a2，a3，a4，a5，a6为待求系数。

其误差方程矩阵形式为:

V=KM-Z

(7)

式中:V为观测值的残差向量，K为待求系数矩阵；M为设计矩阵，Z为观测向量。按照最小二乘原理,VTPV=min的原则，求解得：

K=(MTPM)-1MTPZ

(8)

当偏差数据含有粗差时,如果使用初始权值,会使拟合曲面差生弯曲，影响拟合结果精度[15]。利用抗差估计理论，构造极值函数为：

(9)

式中:Pi为第i个观测值的权重，ρ为凸函数。利用等价权原理，由上式可得：

(10)

单位权中误差：

(11)

(12)

式中, |vi/σ|为第i个观测值经k次迭代后的标准化残差；根据经验c0一般为1～1.5,c1为3.0左右。通过抗差最小二乘曲面拟合获得偏差拟合曲面，将辅轨道整体矫正，使其与主轨道参考基准统一，最终实现升降轨PS-InSAR地面沉降监测结果数据融合。

3 升降轨解算结果互检与融合

3.1 升降轨垂直形变量解算结果

将38景降轨数据和46景升轨数据，利用PS-InSAR技术时序分析，获取了研究区内2018 年1月至2019年7月间地面沉降序列，同时利用式(4)将获取雷达视线向的形变量转换为垂直方向上的累计形变量，如图2所示。可以看出，研究区内分布着多个沉降中心，最大沉降量为-459 mm，除沉降中心外，其他区域的沉降量均小于50 mm。

3.2 升降轨解算结果内部互检

通过获得升降轨公共覆盖区域，提取该区域内升降轨模式下各自的PS点，共提取降轨模式下25 809个PS点，升轨模式下31 049个PS点。由于雷达入射角以及采样方式的差异，导致升降轨解算结果中同名相干点数极为有限，为实现升降轨解算结果互检，本研究通过分别获取抬升区和沉降区各30对同名相干点用于对比分析实验。如图3所示，沉降区内升降轨垂直形变量相关系数为0.866 5，线性拟合方程为y=0.941 4x-45.264 3； 抬升区内升降轨垂直形变量相关系数为0.802 6，线性拟合方程为y=1.056 7x-33.536 8，表明升降轨模式下同名点垂直形变量具有较高的相关性。

(a)降轨垂直向形变量

(b)升轨垂直向形变量图2 升降轨垂直向形变量图Fig.2 Vertical Deformation of Lift Rail

3.3 升降轨解算结果融合方法对比

通过对升轨解算结果进行插值处理，得到与降轨解算结果同名相干点的偏差数组，利用平均值法和抗差最小二乘拟合法分别对偏差数组进行处理，并对升轨解算结果进行校正，使升降轨解算结果参考基准统一，实现升降轨解算结果的数据融合。从图4中可明显看出，图4(a)中1、2区域存在多个与周围不同形变量级的PS点，而图4(b)中1、2区域几乎所有PS点形变量级保持一致，因此，抗差最小二乘拟合法升降轨融合效果明显优于平均值法。选取研究区域内升降轨均匀分布的30对同名相干点，将降轨相干点垂直形变量设为真值，分别验证平均值法和抗差最小二乘拟合法对升轨相干点的矫正程度，如图5所示。通过三组数据对比计算得到，平均值法条件下升降轨同名相干点均方根误差为6.8 mm，而抗差最小二乘拟合法条件下升降轨同名相干点均方根误差为3.4 mm，表明在校正升降轨参考基准中抗差最小二乘拟合法较平均值法精度提高了3.4 mm。

(a)沉降区线性拟合方程图(b)抬升区线性拟合方程图(c)沉降区垂直形变量相关性图(d)抬升区垂直形变量相关性图图3 升降轨垂直形变量相关性图Fig.3 Correlation diagram of vertical deformation of lift rail

(a)平均值法升降轨融合图(b)抗差最小二乘拟合法升降轨融合图图4 升降轨解算结果数据融合方法对比图Fig.4 Comparison chart of data fusion methods for lifting track solution results

图5 平均值法和抗差最小二乘拟合法对比图Fig.5 Comparison chart of mean method and robust least squares fitting method

4 地面沉降监测结果与成因分析

利用抗差最小二乘拟合法融合升降轨解算，对东营市进行地表沉降监测，图6为2018年1月至2019年7月东营市地表垂直形变量图。东营市地面沉降监测结果显示，东营市主要沉降区位于东部沿海区域，城镇沉降区主要位于广饶城区，主要存在以下4个较为明显的沉降区：

图6 融合升降轨数据东营市地表垂直形变量图Fig.6 Surface vertical deformation map of Dongying city from ascending and descending data fusion

1)A沉降区位于河口区仙河镇东部及西北部，其最大沉降量为-458 mm；B沉降区位于东营市胜利机场东部及东北部，其最大沉降量为-202.3 mm；C沉降区位于广南水库南部广饶盐场附近，其最大沉降量为-290.6 mm。A、B、C三个沉降区皆存在着大范围的盐场、盐田，沿海盐场的主要生产模式是通过抽取大量的地下卤水来晒盐。因此，地下卤水的大量开采是造成三个区域沉降的主要因素。

2)D沉降区位于广饶县城区北部及石村镇，其最大沉降量为-187.5 mm，沉降面积约为40 km2，为东营市面积最大的沉降区域，其呈现整体均匀沉降的趋势。该区域存在多个工业园区，对水资源需求量较大，导致该区域深层地下水大量开，。广饶县深层地下水过度开采情况一直存在，其年平均开采量在3000万m3以上 。因此，过度开采深层地下水是造成D区域沉降的主要因素。

根据A、B、C、D四个沉降区的分布特点，可以得知东营市地面沉降呈现东部沉降大于西部沉降，沿海地区沉降大于城区内陆沉降。除过度开采卤水和深层地下水之外，影响东营市地面沉降的因素还有很多，东营市在地质构造上位于郯庐断裂带的西侧，断裂带活动较为频繁，在断裂带的影响下，地壳发生持续缓慢的凹陷，导致东营市地面沉降现象不断发展。从地质土层方面来看，由于东营市地处黄河三角洲，受黄河冲积影响，沉积层多为黏性土层，其自然固结时间短，土质松软易压缩，容易发生地面沉降现象，随着东营市经济社会的不断发展，地面建筑的不断增多，加速了土层的压实作用，加剧了地面沉降的程度。

5 结论

采用PS-InSAR技术获取东营市地表形变量，并利用抗差最小二乘拟合法进行升降轨解算结果融合，实现PS-InSAR解算结果的内部互检及形变结果分析。

1)分别利用升降轨数据PS-InSAR技术进行地面沉降监测，实现解算结果的内部互检，升降轨解算结果中同名相干点垂直形变量相关系数在0.8以上，保证了地面沉降监测的可靠性。

2)利用平均值法和抗差最小二乘拟合法修正升轨解算结果参考基准，通过与降轨同名相干点位比较，抗差最小二乘拟合法均方根误差3.4 mm，优于平均值法的6.8 mm，提高了融合升降轨PS-InSAR解算结果精度。

3)通过融合升降轨PS-InSAR获取东营市地面沉降信息，加密单一轨道的观测序列，丰富了地面相干点的空间分布，监测结果显示：2018年1月至2019年7月期间东营市存在4个明显沉降区，最大沉降量为-458 mm，最大沉降面积为40 km2。东营市沉降整体呈现空间分布不均匀，且东部沿海大于西部城区，其沉降成因主要与过度开采地下卤水与深层地下水有关。