0.   引言

地面沉降是指由于自然因素或人类工程活动引发的松散层固结压缩并导致一定区域范围内地面高程降低的地质现象，是一种不可补偿的永久性环境和资源损失^[1]。我国沿海的珠三角地区不仅是经济发达地区，也是地面沉降灾害的典型高易发区^[2]。东莞市处于粤港澳大湾区的核心区域，地质条件复杂，地下水位高，在西部水乡片区和滨海片区广泛分布有10～30 m的深厚海陆交互相的淤泥类软弱土层。近年来随着城市建设规模和资源开发力度的逐渐增大，引发了大量地面沉降问题，产生了建筑物及道路的不均匀沉降、地下管线破坏等一系列难题，对城市环境地质及国土空间安全造成严重不利影响。

目前，对地面沉降的观测研究常用的手段有水准量测、GPS量测、三角高程量测等，但存在耗时费力、成本高等方面的缺点。合成孔径雷达干涉技术（InSAR）使用两个或多个雷达图像，利用返回卫星的波的相位差异来计算目标地区地表三维位置及变化信息，具有非接触、高精度、范围广、全天候、全天时、成本经济等优点，已广泛应用于地表变形量测、灾害调查及监测预警中^[3-4]。但传统差分干涉合成孔径雷达（D-InSAR）容易受时间、空间失相干影响^[5]。在此基础上，广大学者又陆续提出了永久散射体技术（permanent scatters InSAR, PS-InSAR）^[6-7]、小基线子集技术（small baseline subset InSAR, SBAS-InSAR）^[8]、斯坦福大学PS法（stanford method for persistent scatterers, StaMPS）^[9]等多时相InSAR技术，充分有效地利用多景雷达影像，利用时序分析技术解决传统InSAR技术中的难题，使监测精度得以提升，在北京平原区及通州地区^[10-11]、黄河三角洲地区^[12]、延安新区^[13-14]、西安市^[15]、昆明市^[16]等多个城市的地面沉降监测中得到了广泛的应用，为灾害防治提供了有力支撑。但目前的研究主要针对永久散射体（PS-InSAR）或小基线集技术（SBAS-InSAR）等单一技术方法在某地区展开应用，缺少对不同技术方法的融合使用及交叉验证；尤其是针对地表环境复杂的广大城镇化地区，需要进一步探讨提高监测精度的组合方法及应用实践。

在此背景下，针对东莞市地理环境复杂，植被覆盖茂密，永久散射点目标稀疏，无法选取足够密度PS点的难题，本论文提出了一种改进时序InSAR技术以增强PS点的数量及密度，基于差分InSAR技术获得多时相干涉对，并将SBAS-InSAR技术和PS-InSAR技术结合起来对时序形变进行解算，实现提高地表形变监测精度的目的。以东莞市全市域范围为研究对象，对覆盖东莞市的137景单视复SLC数据（single looking complex，SLC）进行分析，对2015年6月至2020年6月的地表形变展开动态监测，以研究东莞市地面沉降的时空演变特征及影响规律，为地区地面沉降灾害防治提供重要的技术参考。

1.   研究区概况

东莞市位于广东省中南部，珠江口东岸，全市陆地面积2460.1 km²，研究区范围如图1所示。东莞市位于北东东向罗浮山断裂带南部边缘的北东向博罗大断裂南西部、东莞断凹盆地中，构造形迹较复杂。地势东南高、西北低，地貌以丘陵台地和冲积平原为主，其中三角洲平原位于东莞市西北部，发育有海陆交互相的强透水砂层及软弱淤泥层，厚度为10～30 m，对城乡规划及建设影响较大^[17-18]。随着城市的发展，大量道路交通、房屋建筑、基础设施工程逐步推进，地表地质环境平衡遭到破坏，城区多地出现了明显的沉降现象。为保证城市健康可持续发展，需进一步探明沉降灾害发育分布和动态形变特征。

图  1  研究区位置及范围

Figure  1.  Location of the study area

下载: 全尺寸图片 幻灯片

2.   数据准备及方法原理

2.1   数据来源

本文选用覆盖东莞地区的Sentinel-1雷达影像数据对地面沉降开展精细探测。Sentinel-1是一个全天时、全天候雷达成像系统，搭载C波段的合成孔径雷达，是欧洲委员会（EC）和欧洲航天局（ESA）针对哥白尼全球对地观测项目研制的首颗卫星，具有多种成像方式^[19]。本文使用2015年6月至2020年6月采集的137景干涉宽幅（IW）模式的SLC数据，具体参数见表1。影像覆盖区域如图2所示。

表  1  Sentinel-1卫星数据参数表

Table  1.  Parameters of Sentinel-1 satellite

Sentinel-1	参数值
轨道高度/km	700
重返周期/d	12
入射角/（°）	29～46
分辨率	5 m×20 m
幅宽/km	250
极化方式	HH+HV、VH+VV、HH、VV
影像时间	2015年6月15日—2020年6月12日
影像数量/景	137

下载: 导出CSV 

| 显示表格

图  2  研究区及影像覆盖范围

Figure  2.  Study area and image coverage

下载: 全尺寸图片 幻灯片

2.2   改进的时序InSAR方法原理

目前在城市地面形变观测广泛应用的永久散射体雷达干涉量测技术（PS-InSAR）基于高相干性和强稳定性的永久性散射体（即PS点，如人工建/构筑物、基岩等），采用以公共主影像为基准的点目标干涉相位的相位建模和模型参数迭代回归分析，求解平均形变速率和时序形变，因此，PS点的数量和准确选取成为影响其监测效果的因素之一。

在东莞市植被覆盖较广泛、地形起伏较大、地物及地理环境复杂的城镇区域展开地表沉降监测工作，采用传统的PS-InSAR技术通常无法获得数量足够的目标点，从而降低现场监测的效果。为了提高本区域地表形变监测精度，引入了改进的时序InSAR技术，其主要工作原理是在PS选点阶段增加选取的点目标数量，将SLC影像进行差分干涉分析，基于平均相干性指标优选PS点，从而获得数量足够的高相干点目标集^[20]。

同时，为进一步提升图像解译的效果，本研究采用解空间搜索法求解形变量测参数^[21]，根据大气相位和噪声的时空特征，利用时空滤波分理处非线性变形相位，获得目标点的大气相位差、形变相位差和高程差相位等特征数值，解算得到目标区地表形变速率和时间序列，为沉降变形演化特征研究提供数据依据。

2.3   数据处理及分析流程

2.3.1   图像预处理及差分干涉

对前文收集的137景影像数据采用D-InSAR“二轨法”处理生成差分干涉图。影像配准的精度优于0.1个像素，以美国宇航局（NASA）SRTM-30 m DEM数据模拟地形相位。按常规的方法选取高质量PS点后，将SAR影像按距离向和方位向5×1做多视处理^[20]，选择任意一景影像将其与时间上相连至少两景匹配生成干涉对，计算其平均相干性；设定相干性阈值为0.38，大于其设定值的像素作为补充PS点^[22]，从而实现获取的PS点数目更多，点目标密度分布更均匀的效果。覆盖东莞市的InSAR干涉对组成的小基线集网络连接图如图3所示。之后对干涉图进行滤波和相位解缠，得到差分干涉解缠图。

图  3  干涉对小基线集网络连接图

Figure  3.  Small baseline set network connection diagram of interference pair

下载: 全尺寸图片 幻灯片

干涉图中各像元的相位可以表示为：

式中：v——线性形变速率；

$\varDelta \textit{z} $——高程残差；

r——卫星到像元的距离；

θ——卫星入射角；

$\varDelta {\phi }_{{\rm{atm}}} $——大气延迟；

${\varDelta }_{{n}_{j}} $——噪声和非线性形变。

2.3.2   去除地形残差和大气误差

去除轨道误差后对干涉对进行滤波，然后利用线性回归模型估计高程残差。求得符合阈值条件的高程残差后，计算不能被模型表示的残余相位，其中包括大气延迟相位、非线性形变和噪声相位。将高程残差与原始DEM相加以更新高程数据。为了获取准确的大气延迟误差，在去除残余相位中与地形相关的大气误差后，对残余相位进行空间滤波去除其中的噪声，最后得到的即为大气相位（非与地形相关）。

为了获得更精准的结果，本文共进行三次线性回归去除误差，以保证干涉图中的地形残差和大气相位（包括与地形相关的和与地形非相关的大气相位）已经基本去除。在进行第三次线性回归时可求解线性形变速率和残余地形以及大气相位。图4表示的三次线性回归后的高程残差示意图。

图  4  三次线性回归后的高程残差示意图

Figure  4.  Schematic diagram of elevation residual after three linear regression

下载: 全尺寸图片 幻灯片

2.3.3   形变序列求解

此时分离出的大气延迟相位并不十分准确，需要将大气相位加回已去除地形残差的干涉图中，然后根据不同的解算模型，利用奇异值分解法求解每个点的形变序列和形变速率。

数据处理分析流程见图5所示。可以看到，本文的改进时序InSAR技术结合了PS-InSAR和SBAS-InSAR两种技术的优势，充分利用收集的多景影像数据和处理获得的均布的高相干目标点，通过干涉相位和垂直基线的线性回归求取高程残差，通过干涉相位和时间的回归求取线性形变速率，然后通过滤波去除噪声获取大气延迟相位^[23-24]。最终将大气延迟相位、形变相位和高程残差进行分离，利用模型解算时序形变，实现对地面形变更为精细的观测。

图  5  改进的时序InSAR技术流程图

Figure  5.  Flow chart of improved time-series InSAR

下载: 全尺寸图片 幻灯片

3.   结果与讨论

3.1   沉降空间分布

通过改进时序InSAR技术得到东莞市2015—2020历年的平均形变速率图（图6）。从沉降发育规模来看，大部分区域的年平均沉降速率在20 mm/a以内，小部分区域的年平均沉降速率超过20 mm/a。

图  6  东莞市2015—2020年间各年年平均形变速率

Figure  6.  Annual average subsiding rate between 2015 and 2020

下载: 全尺寸图片 幻灯片

对东莞市沉降变形规模及面积进行统计分析如表2所示。可以看出，监测范围内地表沉降总体较为稳定，其中沉降速率低于10 mm/a的弱发育及不发育区面积为2432.3 km²，占监测总面积2460.1 km²的98.87%；变形速率大于10 mm/a的区域约27.8 km²，占监测总面积的1.13%；近5年时间内，平均沉降速率大于30 mm/a的强发育区面积为0.95 km²，约占总监测面积的0.04%。

表  2  东莞市沉降速率统计表

Table  2.  Statistics data of annual average subsiding rate

发育程度	年均形变速率/（mm·a−1）	分布面积/km2
		2015年	2016年	2017年	2018年	2019年	2020年	平均值
弱发育	0～−10	2444.27	2453.04	2444.70	2409.44	2413.01	2429.43	2432.31
中等发育	−20～−10	13.70	6.41	11.76	45.35	42.83	26.59	24.44
	−30～−20	1.46	0.53	2.47	4.17	3.08	2.66	2.40
强发育	−40～−30	0.51	0.11	0.80	0.80	0.79	0.82	0.64
	−40及以上	0.16	0.01	0.37	0.34	0.39	0.60	0.31

下载: 导出CSV 

| 显示表格

进一步得到2020年1月至2020年6月期间全域典型平均沉降速率图（图7）。由图可以看出，市域内地面沉降灾害发育较为严重的区域主要集中在东莞市长安镇地区（图7中的A区）、淡水河沿岸地区（图7中的B区）、万江街道地区（图7中的C区），以及麻涌水道（图7中的D区）。主要沉降区的分布与图6呈现出相同的分布规律。

图  7  东莞市主要沉降变形区分布

Figure  7.  The distribution of serious subsidence area

下载: 全尺寸图片 幻灯片

3.2   沉降时序演化分析

为进一步研究东莞市地面沉降长时间序列演化特征，以虎门水道入海口的长安镇地区（图7中的A区）为主要特征点展开分析。该地区位于东莞市西南部，地层为冲积平原并且存在大面积滩涂。图8（a）—（f）显示了该地区2015—2020年地表年平均形变速率变化趋势；图8（g）为（e）图中白色监测点E所示位置处的2019年累计形变时序。该地区自2015年起就存在持续性的地面沉降，且变形区域范围逐渐增大，沉降较大区域年平均形变速率已达到了60 mm/a，所选时序点在2019年的累计形变量达到了60 mm。该沉降区域集中分布在近海沿岸地区，由影像分析及野外调研可知该区域为东莞市主要填海造地区，推测随着地基逐渐往外扩建，填海区逐年扩张，预测在后续城市建设阶段也存在着较为显著的变形增大趋势，需重点关注和防范。

图  8  长安镇典型监测区的时间序列形变特征

Figure  8.  Subsiding characteristics of time-series in typical monitoring area of Chang’an Town

下载: 全尺寸图片 幻灯片

3.3   地面沉降调查验证

通过野外实地调查进一步验证InSAR监测结果的准确性和可靠性。现场沉降行迹调查点位于长安镇内有色金属材料交易中心前，佳能二手汽车店附近。该区域建筑下部地基为珠江口软土地层。现场观测发现房屋地基与道路交叉处存在明显分裂，路面下沉较严重，最大差异沉降超过200 mm，地基差异沉降导致房屋开裂，裂缝宽度达到10 mm，如图9所示。形变监测结果与野外调查情况吻合。

图  9  长安镇调查点地面沉降变形特征

Figure  9.  Characteristics of land subsidence deformation in Chang’an Town

下载: 全尺寸图片 幻灯片

3.4   地面沉降成因分析

3.4.1   地质因素

东莞市软土地面沉降主要受软弱土层分布控制。软土主要为淤泥、淤泥质土，粒度成分以粉粒、细粒为主，矿物成分中除长石、石英、云母外，含有大量的黏土矿物；多呈疏松多孔的蜂窝状结构。具有含水率高（w=42.4%～91.8%）、孔隙比高、压缩性高（α₁₋₂=0.51～3.25 MPa⁻¹）、灵敏度高、密度低、渗透性低（k=7.67×10⁻⁷～5.32×10⁻⁵ cm/s）、强度低（c<20 kPa）、欠固结等特点。部分地区软弱土层还含有砂砾夹层，排水能力较强，在建筑及车辆等荷载的作用下，孔隙水被不断挤出，体积压缩，地层发生固结压缩形变。现场调查，软土的分布主要在水乡片区及滨海片区，软土空间分布也与地面沉降的主要发育区高度重合。图10展示的是滨海片区（图7中A区域）累计地面沉降与软土厚度等值线的对比图，可以看到该区域软土厚度均大于10 m，表现为较大的沉降变形规模。

图  10  重点区软土厚度与InSAR解译累计沉降量对比图

Figure  10.  Comparison of soft soil thickness and InSAR interpretation accumulated settlement in main investigation areas

下载: 全尺寸图片 幻灯片

3.4.2   人类工程活动

附加荷载和地下资源开采是东莞市地面沉降形成的重要外在因素。附加荷载主要表现为城市建设带来的建筑物的静荷载和车辆的动荷载两类。在沉降量大于100 mm的调查点中，140个为静荷载作用，12个为动荷载，荷载对地面沉降的影响明显。

此外，东莞市河网密布，地下水位呈季节性周期变化，居民用水、水产养殖用水、工厂抽水、地下矿产开采也导致地下水位变化。在枯水季地下水反向补给地表水，含水层水量减小，水位降低，土层发生固结，加之黏性土在释水过程中具有单向不可逆性，继而引起较显著的地面沉降。

4.   结论

本文采用改进时序InSAR技术对东莞市2015年6月—2020年6月的137景Sentinel-1影像数据进行处理，对东莞市地面沉降时空演化特征进行了监测分析，得到如下结论：

（1）东莞市地面沉降发育面积达503.85 km²，占全市域面积的34.6%，是研究区典型的地质灾害类型。从发育规模上，累计沉降量≤100 mm的范围面积382.84 km²，位于100 mm和300 mm的范围面积135.81 km²，≥300 mm的范围面积1.35 km²。软土地面沉降主要集中发育于麻涌镇、道滘镇、洪梅镇、中堂镇、沙田镇及滨海湾新区。

（2）在研究的6年内，全区地面变形整体较稳定，大部分沉降点处于缓慢发展阶段，部分处于强烈发展阶段。沉降速率低于10 mm/a的弱发育区面积2432.3 km²占监测总面积的98.87%；强发育区，沉降速率超过30 mm/a的面积为0.95km²，约占总监测面积的0.04%，需进一步严密监测其变化动态。

（3）通过地表形变监测和野外调查发现，地面沉降与当地的地质条件（如淤泥质软土的压实固结）和人类活动（如建构筑物施工、农业灌溉、鱼虾养殖等）的耦合作用有很强的时空相关性。经现场验证，采用改进时序InSAR技术的形变监测结果与野外调查情况吻合，可为后续灾害防治提供技术参考。