黄龙霄，张旭晴，赵强，高明久，程微，安继魁，吴迪

1.吉林大学 地球探测科学与技术学院,长春 130026；2.辽宁工程勘察设计院，辽宁 锦州 121000

0 引言

地面沉降又称地面下陷，是由于地下松散地层固结压缩，导致地壳表面标高降低的一种局部的工程地质现象[1]。合成孔径雷达干涉测量技术(Interferometric Synthetic Aperture Radar，InSAR)是近20年发展起来的一种先进的遥感技术，补充了已有的地面沉降监测方法:精密水准测量和GPS监测[2]。Gabriel等人将D--InSAR技术应用在地面沉降监测领域[3]，D--InSAR技术虽然可以获取密度较大的形变值[4]，但它着重研究短时间间隔内的单一形变情况。为了弥补D--InSAR技术的不足，学者Ferretti提出PS--InSAR技术[5]。PS--InSAR 技术只关注稳定性较高点，可以进行长时间序列的地面沉降监测，但在地表情况变化较大的一些非城市区域使用PS--InSAR技术时很难获取高密度的沉降监测结果。因此，Berardino P et al.提出 SBAS--InSAR 技术[6]，该技术利用时空基线较短的SAR数据集形成干涉对[7]，充分利用短时空基线影像的相干信息，有效地抑制相位噪声对地形相位的影响，从而获取长时间缓慢地表形变的演变规律[8]。SBAS技术不仅可以减轻时空失相关、大气延迟的影响，还可以应用有限数量的影像得到毫米级的时序沉降量[9]，空间分辨率的降低也相应的减少了数据处理的运算量。

由于辽宁义县中东部地区地形较复杂，为了得到好的监测结果，综上考虑，采用SBAS--InSAR技术对该区域进行地面沉降监测研究。

1 SBAS--InSAR技术原理

SBAS技术的主要思想为通过设置时间和空间基线阈值得到L个小基线集合，每个集合之内，干涉对基线较小，而集合之间，干涉对基数较大，共包括M幅差分干涉图，筛选出基线较短的干涉对组成不同的集合。集合与集合之间的基线较长，可通过使用最小二乘法得到每个集合的形变序列，再通过奇异值分解法将每个集合联合求解，这种方法可以在有限的数据基础上得到高度精度反演，同时不要求主影像必须相同[10--13](图1)。

图1 技术流程Fig.1 Technical process

选取同区域的N+1幅SAR影像，获取时间依次为t0,t1,t2,…,tn,根据设置的干涉条件可以组合得到M幅干涉图，其中M满足下式：

(1)

现选取tA和tB两时刻获得的SAR影像，然后进行差分干涉处理，形成第i个差分干涉对，则差分干涉相位δφi可由下式表示：

(2)

式中：d(tA)和d(tB)表示tA和tB时刻相对于参考时刻t0(d(t0)≡0)的累积形变信息，φ(tA)和φ(tB)表示相应相位值。

现用向量φT=[φ(t1),…,φ(tN)]表示N个相位图，向量δφT=[δφ1,…,δφM]x表示M个干涉相位图，如果主图像先于辅图像的获得时间，则有：

δφj=φtSj-φtMmj,j=1,…,M

(3)

式中：Mm、S分别表示干涉像对中主、辅图像。所有参与生成干涉图的SAR影像可以在不同的短基线集中。现定义一个含有M个方程，N个未知参数的方程组，其表达式为：

δφ=Aφ

(4)

式中：矩阵A[M×N]中每一行对应一个干涉对，每一列对应一副SAR影像，A[j,Sj]=1和A[j,Mmj]=-1，其余元素为0。矩阵A是一个由干涉图组合方式决定的近似关联矩阵。当所有干涉对属于同一个基线集时，A是一个列满秩矩阵。当M=N时，方程组(8)有唯一解；当M>N时，方程组是超定方程，此时可使用最小二乘法求此方程组的唯一解：

(5)

当所有干涉对属于不同基线集时，矩阵A的秩等于N-L+1，出现秩亏，可采用奇异值分解(SVD)的方法计算矩阵A的广义逆矩阵A+，进而求解方程组(4)的最小范数解，其过程如下：

A+=US+VT

(6)

式中：正交矩阵U[M×N]，其前N列是ATA特征向量,称为矩阵A的左奇异矩阵，正交矩阵V[M×N]，其所有列ATA称为矩阵A的右奇异矩阵。S为M×N矩阵。则：

A+VS+UT

(7)

(8)

式中：ui和vi分别表示正交矩阵U和V的列向量。

将方程组(4)中对相位的求解转化为对平均相位变化速率的求解问题，则待求参数向量为：

(9)

则有：

δφ=Bv

(10)

式中：B[M×N]，其中,B[i,j]=tj+1-tj,(Sj+1≤j≤Mmi,i=1,…,M)其他元素值为0。对B进行奇异值分解，即可解出各时间段平均速度v。

2 地面沉降监测

2.1 研究区概括

辽宁义县东部为剥蚀构造中低山区，由变质岩和花岗岩侵入体组成，山势较陡峭，海拔在400～800 m之间，山脉走向由东北向西南延展，山体多为直坡，局部为凸坡，坡度角为30°～50°，山间谷地，多呈树枝状分布。义县的中部为冲积平原区，分布于大凌河、细河东岸；上部为亚砂土覆盖；下部为砂砾石层；底部混土，地面较平，微倾向河床。中部区域部分岩层性质为白垩系：灰白色砂砾岩夹砂岩和页岩及上部夹薄层煤，中部赋存煤层(图2)。

图2 研究区域Fig.2 Research area

2.2 监测数据

Sentinel--1数据具有较好的空间基线控制，并且重返周期较短，适合地表形变监测[13]。由于辽宁地处中国东北部，冬季多雪，而C波段影像受降雪影响会出现严重的失相干现象，导致成果获取的相干点数量降低，误差增大，所以在本次监测中未采用冬季影像。最终选取由欧空局提供的IW宽条干涉的降轨SAR影像，C波段，分辨率为5 m×20 m，极化方式为VV，2017年—2018年2月到11月的18景影像数据，具体时间见表1。同时使用Sentinel--1卫星精密轨道数据(POD精密定轨星历数据)和由地理空间数据云将ASTER GDEM(V1)数据进行加工得来的GDEMDEM 30M 分辨率数字高程数据作为辅助数据。

表1 Sentiel--1数据

2.3 处理流程

利用SAR Scape 软件，通过设置时间基线、空间基阈值控制生成干涉对的数量。本次监测设置270 d的时间基线和最大45%的临界空间基线阈值，生成88对像对，生成的时空连接如图3所示。

图3 时空连接图Fig.3 Space-time connection map

然后基于干涉对进行SLC影像配准，干涉图的生成、去平，经过自适应滤波、相干性生成及相位解缠获取一系列相位图。解缠相干系数阈值为0.2，解缠方法为Delaunay MCF，滤波方法为Goldstein。然后再选择其中1景干涉效果中等的干涉图，根据相干性图选取相干性高、相位好的GCP点，为轨道精炼和重去平做准备。选取近100个控制点，通过多次筛选将各个点的误差降到了1.5以下。然后经过两次SBAS反演，精确估计且去除地形残余相位、大气效应相位，最后结合DEM数据进行地理编码后获取WGS--84坐标系下的各期累积形变图和平均形变速率图[1,10]。

2.4 监测结果与分析

本次监测使用SBAS--InSAR技术通过对Sentinel--1数据的处理，结合ArcGIS软件，得到2017—2018年成果沉降速率渲染图(图4)、沉降速率统计图(图5)和2017—2018年成果沉降速率等值线图(图6)。

图4 2017年—2018年成果沉降速率渲染图Fig.4 Rending map of results settlement rate from 2017 to 2018

图5 沉降速率统计图Fig.5 Statistical map of settlement rate

由图4可知，在研究区的左下角部分即义县县城周边存在较明显的沉降，义县中心地区沉降最为集中，由图5可知，2017 年 2 月 21日—2018年 10月26日期间，研究区范围内大部分PS点的沉降速率为5 mm/a±，平均沉降速率为3.5 mm/a，最大沉降速率为 98 mm/a，说明研究区地面沉降空间分布有较大差异，地面沉降不均衡。

使用克里金法进行插值，绘制出该区域的形变量等值线图(图6)，叠加到谷歌地图上，便于沉降区定位；发现在A、B、C 3个小区域存在沉降漏斗且沉降量较大，对这3个小区域进行重点分析，从这3个区域中分别选择一个能够代表本区域形变特征的、形变量较大的 PS 点进行时间序列分析。 时间序列分析是处理变形观测数据的一种有效方法，通过时间序列分析可以发现各形变点的时间变化规律，可以对沉降区域进行短期的预测与分析[14]。由图7～9各曲线可得，沉降区A的沉降中心沉降量>-180 mm；沉降区B最大沉降量大约为259 mm；沉降区C有3个较明显的沉降中心，沉降量>-275 mm。

图6 2017年—2018年成果沉降速率等值线图Fig.6 Contour map of results settlement rate from 2017 to 2018

图7 A区PS点形变时间序列Fig.7 Deformation time series of PS point in region A

图8 B区PS点形变时间序列Fig.8 Deformation time series of PS point in region B

图9 C区PS点形变时间序列Fig.9 Deformation time series of PS point in region C

图10 义县地质灾害详细调查图Fig.10 Detailed survey map of geological hazards in Yixian

为了检验结果的可靠性，将图10与监测结果对比，发现义县地面塌陷情况与所得结果基本吻合。地面沉降主要是由地壳构造活动、矿产资源开采、地下水过量采集与城市建筑荷载增加等因素造成的。为分析义县地面沉降成因，通过查阅相关资料发现：研究区内矿产资源较丰富，有煤矿21个、黏土矿10个、采石厂11个、铁矿2个和硅石矿1个(表2)。义县矿业开采强烈，尤其是对地质环境影响较大且易发生灾害的煤矿、花岗岩矿的开采，其中对煤矿、黏土矿的开采会造成地面塌陷、地裂缝灾害。据统计，截至2004年义县中部地区的聚粮屯煤矿和九道岭煤矿采空区，使2 km2多耕地成为沉陷区。一些个体企业无规则开采对山体破坏严重，极容易发生崩塌、滑坡等地质灾害。沉降区的分布与义县的矿区分布较为一致，且处于采空塌陷易发区，表明研究区地面沉降主要是由人类工程经济活动导致的，其中矿产资源开采是最大的影响因素。

表2 矿产资源汇总

3 讨论

本次监测改进之处：①数据处理过程中为了能够选出形变量更小、稳定性更好、误差更小的 GCP 点，通过控制每个GCP点的误差值，使结果更加的精确，减弱了人为因素对轨道精炼误差的影响。②监测结果分析采用了等值线分析法、时间序列分析法和统计学法，从不同的角度研究了义县地面沉降情况，成功地获取研究区的平均沉降速率和形变序列并结合野外实地验证，确定了主要沉降区的位置、最大沉降量以及沉降趋势。

4 结论

(1)监测发现义县地表沉降分布不均匀，在义县中部地区存在多处较明显的沉降区，如聚粮屯乡和前杨乡附近。在监测时段内研究区沉降量在-275～5 mm范围内。

(2)义县的沉降区与矿产资源的分布有着密切的关系。在矿石资源开采中，土体开挖卸载时开挖面土体向开采面内移动、地下结构整体下沉均可造成地表沉降。

(3)义县地表沉降趋势显示下降，时间序列与季节变化存在一定的相关性，雨水充沛期和冬季降雪期，沉降速率都会出现减缓或与总趋势相反的小幅度变化。