利用时间序列InSAR技术对珠江三角洲主要城市进行地面变形监测研究

戴宜韦，吴希文，闵新颖，王 华，朱焕廉，彭林才

（广东工业大学 土木与交通学院，广东 广州 510006）

地面沉降是全球性问题，由于人类活动的影响，加剧了这种地质灾害发生的程度与范围. 现如今，全球性气候变暖明显，导致海平面上升，地面沉降对沿海地区的影响尤为严重. 据国土资源部地质灾害应急办不完全统计，仅仅2017年上半年我国由于地质灾害造成的直接经济损失就达到2.6亿元.

近年来，随着合成孔径雷达技术(InSAR，Synthetic Aperture Radar Interferometry)的发展，InSAR是一种利用微波波段特性进行相干成像系统技术，在监测地面沉降方面已经成为一种可靠的新手段. 1999年，Ferretti等提出永久散射体干涉测量PSInSAR(Permanent ScattererInSAR)技术，利用多景SAR图像从中找出散射特性稳定的高相干点作研究对象，并成功获取长时间序列的大范围地面形变缓慢的演化过程[1]. 在先前研究中，2014年王华等使用美国喷气推进实验室(JPL)和加州理工学院联合开发的ROI-PAC软件对覆盖广佛地区的Envisat卫星影像进行D-InSAR干涉图生成，并利用π-rate软件对干涉图进行后续分析，获得了广佛地区的地面沉降速率与时间序列，验证了利用InSAR技术在广佛地区的沉降监测可以达到毫米级精度，广佛地区趋于平稳状态，局部地区发现有沉降[2]. 2018年Ng，A.H.M等利用时间序列InSAR原理对广佛地区2011年至2017年COSMO-SkyMed卫星数据进一步进行地面沉降分析，整体区域的沉降速率范围在-35～10 mm/a之间[3].由于原先的雷达数据很大部分都是需要购买的或需要向卫星所属机构申请获得，有一定的局限性. 随着2014年Sentinel-1A卫星的发射成功，2015年4月5日开始稳定运行，数据免费向公众开放，利用Sentinel-1A数据监测地面沉降可以进一步减少所需费用. 在过去的研究中并没有人利用过Sentinel-1A来研究珠江三角洲广佛及东莞地面变形情况，本文研究证明了Sentinel-1A可以应用于珠江三角洲地面变形监测，为以后Sentinel-1A应用于监测气候湿润的三角洲地区的地面变形提供了参考依据.

本文利用时间序列InSAR技术对珠江三角洲主要城市的主城区进行沉降观测，获取了广州、佛山、东莞在内的沉降演变特征，并分析各个城市的地面沉降速率图，预测出潜在的危险沉降点.

1 时序InSAR原理

Persistent Scatterer Interferometry(PSI)技术的基本原理是利用多景同一地区的SAR影像，通过统计分析时间序列上的幅度与相位信息的稳定性，探测不受时间、空间基线去相关影响的稳定点目标[4-5]. 对探测出的稳定目标点进行时间序列分析，从而反演出研究区地物的形变量. 其主要处理在于干涉图生成、利用DEM进行差分干涉图计算、候选点(PSC)选取以及对候选点进一步进行筛选得到最终的PS点.相位解缠所使用的函数模型如式(1)所示.

其中， φtopo表示由参考DEM误差引起的地形相位；φdefo表示点的形变相位；φatmo表示大气延迟相位；φnoise表示不相干的噪声相位. 在观测相位中含有φatmo大气延迟相位，为了去除减弱大气相位对最终结果的影响，通常通过相邻PS点的大气相位空间相关性，采用基于相邻PS点相位差分模型提取地表沉降信息[6]. 一般对PSI做的准备工作包括:基于振幅离差指数方法对PS候选点(PSC)选取，基于整体相干性最大法估计出线性形变速度增量及高程误差以及对大气相位的滤波去除[4]. 对于大气相位的消除，是基于大气相位在空间区域相关度很高，呈现出低频状态，在时间域上大气相位相关度不高，呈现出高频状态的特征[2，4，7]. 从而采用在时间上高通滤波，在空间上采用低通滤波来分离出非线性沉降形变与大气延迟相位;即整体相位减去模型相位得残余相位，残余相位中包含大气相位(APS)、非线性形变与其余随机误差(noise)[2，4，7]. 随机误差对最后的结果影响较少，可以忽略不计. 最终得到研究区的形变相位，进行时序分析.

2 研究区概况

珠江三角洲的植被覆盖率高，空气湿热，年均温度21～23 ℃，年降水量大约1 700 mm，以南亚热带季风气候为主. 多雨季节与高温季节同步，土壤肥沃，河道纵横. 珠江三角洲陆上部分由西江断裂、瘦狗岭断裂、南岗—太平断裂所围拢的断陷盆地组成. 根据广东地质调查局的报告，有三种地质构造，即该地区有软、疏松的沉积物、埋藏的岩溶和其他类型的基岩. 前两个地质单元均有沉降，对周围建筑物的安全构成威胁，基岩由分布在增城丘陵地区和广花盆地北部的不同类型的花岗岩组成[3，8]. 三角洲及滨海平原软土区则有软土地基变形、砂土液化、河港淤积，胀缩土等分布区，常见的地质灾害有地裂缝，部分地区出现地面沉降[9]. 本文研究的主要城市主城区如图1所示

图1 研究区域(广州、佛山和东莞)Fig.1 Research area (Guangzhou, Foshan and Dongguan)

3 Sentinel-1A雷达数据及处理方法

本文数据处理的范围覆盖广州、佛山和东莞等主城区. 本文处理数据采用2015年12月～2018年5月的Sentinel-1A卫星数据，总共66景IW模式SAR影像.具体卫星与数据参数见表1.

表1 Sentinel-1A卫星数据参数Tab.1 Sentinel-1A satellite data parameters

本次研究共处理了66景Sentinel-1A卫星影像. 卫星轨道数据采用欧空局(ESA)发布的精密轨道数据，本文采用ESA发布的精密轨道数据，可以认为轨道残余误差对最终结果的影响很小，进一步忽略轨道误差造成的影响. 采用美国JPL发布的SRTM的1’DEM数据去除地形相位. 本文采用的66景影像及生成65幅干涉图对应的时间基线与垂直基线如图2所示.

图2 Sentinel-1A空间基线与时间基线分布图Fig.2 Sentinel-1A spatial baseline and time baseline distribution map

首先对于时间序列InSAR，选取覆盖珠江三角洲的66幅单视复数图像(SLC)进行辐射矫正，根据时空基线，选取获取日期2017年3月12日作为主影像，其余影像作为辅影像，生成65幅干涉图. 利用通过设定幅度离差阈值方法进行PS候选点选取，选取幅度离差指数低于0.25的PS点构建PS网络. 通过导入外部DEM生成差分干涉图，进行模型的反演，估算出残余地形与形变速率. 基于振幅离差指数[1，6]自动选出一个或者多个像素(参考点)，参考点设置为25 km2一个参考点，计算出相位偏移[6]. 在干涉图生成后，就从所有的干涉图中去除相位偏移. 第一次反演后用线性模型从所有差分干涉图中估算的形变速率和残余高程信息，残差是大气信息(噪声和非线性目标). 对相干性高的像元进行计算，进行第二次反演，得到估算大气延迟相位，去除大气延迟相位得到最终的形变速率[6]；最后对沉降速率、高程残差以及时间沉降量进行地理编码，投影到当前所用的WGS84地理坐标系.

4 研究区沉降监测结果

由于本文所用卫星数据从2015年12月至2018年5月底，所得出的研究成果通过跟2015年之前的研究及本文研究期间他人利用其他卫星平台所发表的研究成果进行比较，从而验证本文研究成果的准确性.

4.1 广佛地区

本文最终计算珠江三角洲广佛主城区地面沉降速率精度误差及精度误差统计如图3所示.

图3 Sentinel-1A广佛雷达卫星影像沉降速率精度图及精度分布频率图Fig.3 Sentinel-1A radar satellite image subsidence rate accuracy map and accuracy distribution frequency chart in Guangzhou and Foshan

从广佛地区沉降速率精度图可以看出，研究区内处理的最终结果精度大部分分布在0.3～0.6 mm/a范围内，标准差(std)0.15 mm/a，平均精度(mean)0.51 mm/a，监测结果稳定可靠. 图4为处理最终的沉降速率分布图，红色代表抬升，蓝色代表沉降，从图中广佛整体的区域分布上研究，城区沉降量不大.

本次研究广佛地区共生成6 833 697个PS点，生成的PS点沉降速率的标准差(std)为2.22 mm/a，平均值(mean)为-0.06 mm/a，可以认为研究区域整体平稳，局部发生沉降，与原先王华等对广州佛山沉降研究结论整体相符[2].

4.2 东莞地区

东莞主城区地面沉降速率精度误差及精度误差统计如图5所示

图5 Sentinel-1A雷达卫星影像沉降速率精度图及精度分布频率图Fig.5 Sentinel-1A radar satellite image subsidence rate accuracy map and accuracy distribution frequency chart in Dongguan

从东莞地区沉降速率精度图可以看出，研究区内处理的最终结果精度大部分分布在0.3～0.7 mm/a范围内，标准差(std)为0.14 mm/a，精度平均值(mean)为0.50 mm/a，可以认为此次处理的结果精度稳定可靠. 此次实验研究东莞地区共生成4 213 293个PS点，沉降速率图如图6所示，标准(std)为1.56 mm/a，平均值(mean)为-0.05 mm/a，可以认为东莞研究区域整体平稳.

图6 东莞地区沉降速率分布图及速率频率统计图Fig.6 Subsidence rate distribution map and rate frequency chart of Dongguan

5 变形成因分析

本次研究重点分析了图4广佛地区沉降速率分布图与图6东莞地区沉降速率分布图中所标注的A、B、F、G、M、N、U及V点.

5.1 广佛地区变形成因分析

据广州日报报道，佛山五区范围内约有6 000～8 000口井抽取地下水，尤其在农村和城乡接合部地区地下水破坏尤为严重[10]. 此前有研究表明特别在工业密集区域及隧道施工区域，由于非法不合理的对地下水抽取，沉降现象更加明显. 例如:2007年因武广高铁施工而导致金沙洲发生数起地面塌陷事故[11-12]，再如2008年7月21日，佛山市顺德区北滘职业技术学校全校几乎所有的建筑都出现了不同程度的沉降现象，并导致许多墙体开裂[2，12]. 广州地下水抽取与北方地区例如北京、天津等城市因地下水开采造成大范围沉降有所不同[13-14]. 北方地区特别是华北地区地下水开采几乎都在深水层，而珠江三角洲地区地下水开采在潜水层. 潜水层可以通过降水或地表水下渗补给，由于地下水开采在广佛主城区地区大部分都是在潜水层，距离地面5 m左右[3]. 所以广佛地区包括周边相邻地区并未出现由于地下水开采而发生的大范围地面沉降，只是出现的局部地面沉降.

图4中A点标注的季华西路沿线一带在研究期间出现重大安全事故，2018年2月7日，该路段下方在建的佛山地铁二号线绿湖岛至湖涌盾构期间发生地面塌陷事故，造成11人死亡、1人失踪、8人受伤，直接经济损失达到5 323.8万元[15]. 此处地陷事故是由于该区间地质条件复杂，存在深厚富水粉砂层且临近强透水的中粗砂层，且地下水具有承压性. 据中央电视台央广网报道由于地铁右线施工突发透水，导致轨道管片变形受损，外侧土壤与管片应力失衡，引发地面季华西路30多米路段坍塌[16]. A点在2015年12月至2018年5月沉降的时间序列图，如图7所示可以明显看出季华西路一带在研究期间出现沉降，最高累计沉降量达到52 mm，2015年12月至事故发生时已出现明显沉降.

图7 A点时间序列图Fig.7 Point A time series chart

图4中B点位于佛山乐从镇大墩村，在建的佛山3号线经过此地，并且附近有明显关于地铁的施工建筑. 主要沉降区域发生在岭南大道南路与荷月路交口东约600 m处及岳步工业区附近一带，其中岭南大道南路与荷月路交口东约600 m处最大沉降速率为-23.17 mm/a沉降点附近交通繁杂，人口居住多，且沉降点西100 m处有工厂，此处安全风险较高. 岳步工业区东北500 m处有明显地铁开挖口，附近施工活动较多. 据南方都市报报道，大墩站于2016年8月31日先行施工. 从图8中可以看出在2016年9月(黑色五角星标注)左右出现沉降速率加快，此后的沉降趋势基本上为与施工周期保持一致. 2017年3月16日(黑色正方形标注)至2019年12月31日大墩村建设工程实施封闭施工，从图8可以看出从2017年3月中旬开始，沉降速率进一步加快. 可以认定B点沉降原因主要在于佛山3号线施工引起，由于地铁施工破坏土壤稳定，造成该处发生沉降. 图8为B点附近一沉降点的累计沉降量达到48 mm.

图8 B点时间序列图Fig.8 Point B time series chart

图4中F点位于佛山钢铁世界至黎湖工业区一带，沉降区域内钢铁厂密集分布，物流公司分布密集，还分布着众多家具厂，监测获得的最高沉降速率为-24.8 mm/a. F点附近两沉降点累计沉降量如图9所示，达到了46.36 mm. G点与F点类似，位于密集的工业区内，周围皮革厂众多，人口密集；G点最高沉降速率达到-25 mm/a，沉降量达到48.86 mm. 从G点与F点位置上分析，此处的沉降点与工业分布十分吻合，初步认定为工业不合理开采地下水及当地人员抽取地下水引起的地面沉降. NG等利用广佛地区2011年至2017年COSMO-SkyMed卫星数据做时序InSAR分析，也显示此处周围有沉降，但未具体分析[3].F点与G点周围部分沉降点累计沉降量如图9所示.

5.2 东莞地区变形成因分析

图9 F、G点时间序列图及工厂分布图Fig.9 Point F, G time series chart and factory distribution map

东莞主城区的沉降速率图如图6所示，整体平稳，局部发现地面沉降，范围较小，只是呈现出点状分布. 近年来东莞市经济飞速发展，有世界工厂之称.如图11所示，研究区内的点状分布的沉降点也主要集中在工厂密集区域内. M点、N点分别位于东莞市与深圳市、广州市接壤处，工业密集，有各种印刷厂、涂料厂、造纸厂、五金厂等等. 此前也有报道证明目前深圳关外地区(东莞)的工厂以前开挖了大量的非法地下水水井，目前很多水井还在运行，无序开采将使地下水环境进一步恶化，可能会导致地质灾害和海水倒灌[17]. 这两块区域中发现的沉降点初步归因于工厂及人口对地下水不合理开采以及人类活动所致. M点与N点区域的明显沉降点时间序列图如图10所示，点M与点N的累积沉降量分别为达到39 mm、38 mm.

图10 M、N时间序列图Fig.10 Point M and N time series chart

M、N、U与V点周围的工业分布如图11所示，图中红色数字代表的是工业聚集点，红色圆圈点代表大工业厂房分布地. 从图11中与图6东莞沉降图可以看出，工业聚集点M、N、U与V都有明显沉降，这同时也印证了文献[3， 7]中工厂及当地人员不合理抽取地下水导致地面沉降的观点.

图11 M、N、U、V点工厂分布图Fig.11 Industrial distribution around the M, N, U, and V points

U点沉降明显区域位于光正实验学校周围至茶山塘工业区一带，学校周围有明显的施工建筑，施工活动频繁. 工业区内密集分布着制衣作坊、纸厂及印刷厂等. V点位于东莞市田坑工业区，周围有村尾第一工业区、四海工业园等，工业区内密集分布着模具厂、塑胶厂、五金厂等各种耗水量巨大的工厂. 此外，此处人口密集，大量人口抽取地下水量也是巨大的，与M点、N点区域极为相似. 从图12可以看出，V点工业区沉降速率更快，超过-18 mm/a，累计沉降量也更为明显，U点与V点累积沉降量分别达到45 mm、52 mm.

图12 U、V时间序列图Fig.12 Point U and V time series chart

6 结论

本文利用Sentinel-1A数据研究了广州、佛山与东莞地区的地面变形情况，结果表明2015年12月至2018年6月期间，整体区域平稳，局部发现明显沉降.将此次研究的广佛地区与其他关于广佛地区的研究进行对比发现，此次研究结果可靠. 局部地面变形区域主要集中在工程施工区与工厂密集区. 本次研究的主要的几个危险性较高的沉降区域中，A点(季华西路)发生了重大事故，B点(大墩村附近)暂未发生事故，但不排除未来有事故发生的可能. 希望并建议有关部门做好风险防控措施，预防进一步的灾害发生.其中此次试验研究发现F点(佛山市六涌工业区与钢铁世界一带)及G点(佛山市亚洲国际家具材料交易中心周围一带)沉降明显，建议有关部门采取相关措施进一步查明沉降原因. 东莞市U点(光正实验学校至茶山塘边工业区一带)、V点(田坑工业区)、M点及N点周围沉降明显，希望有关部门进一步查明原因，给予防控措施.

总体来说，本次实验结果可靠，证明了利用Sentinel-1A也可以较好地监测出珠江三角洲所研究城市的地面变形情况，因本次研究缺少珠三角地区的水文资料、水准数据与后续地铁沿线施工情况数据，后续研究中将继续结合不同资料进一步分析珠江三角洲地区的沉降情况.

致谢：感谢欧洲空间局Sentinel-1A卫星数据的免费使用.