哨兵-1A合成孔径雷达的深圳围填海区域地面沉降监测

周 吕，施显健，任 超，黄远林，朱子林

(1.桂林理工大学测绘地理信息学院，桂林 541004；2.城市空间信息工程北京市重点实验室，北京 100038；3.广西空间信息与测绘重点实验室，桂林 541004；4.北部湾大学资源与环境学院，钦州 535011)

改革开放以来，深圳地区的经济得到了快速的发展，但随着城市的快速扩张和土地需求的迅速增长，深圳土地资源紧缺问题日益凸显，其中围填海成为深圳缓解土地供应的重要途径之一[1]。资料显示，仅在1990—2015年，深圳就通过围填海获得了超过80 km2的土地面积[2]。然而深圳围填海区域的地层结构性质复杂，随着围填海区域地面的逐渐巩固和地面建筑负荷的持续作用，极易引发地面沉降，进而威胁城市发展和人民生命财产安全[3]。因此，及时掌握深圳围填海区域在填海造地工程结束后的地面沉降情况具有重大现实意义。

传统的地面沉降监测方法一般是利用全球定位系统(global positioning system，GPS)测量方法或水准测量方法通过布设地面监测网的方式进行地面形变监测[4]。利用上述方法采集到的数据的精度较高，但是需要投入较大人力和时间成本。近年来快速兴起的合成孔径雷达干涉测量技术(interferometric synthetic aperture radar, InSAR)具有大范围[5]、高精度[6]、长时间监测[7]等优点，为大面积地面沉降监测提供新思路[8]。然而InSAR技术容易受到失相关、大气及轨道误差的影响[9]。为了克服这些问题，中外众多学者在合成孔径雷达差分干涉测量技术(differential interferometric synthetic aperture radar, D-InSAR)基础上提出了时间序列合成孔径雷达(time-series interferometric synthetic aperture radar, TS-InSAR)技术[10-11]，通过解析在时间序列中保持稳定散射特性的像素点来提取地面形变速率[12]，较好地消除了失相关、大气及轨道误差对地面形变监测的影响[13]，实现了较高精度的地面形变时序监测[14]。

深圳围填海区域的SAR数据相对紧缺和商用SAR卫星数据定制价格高昂的问题在某种程度上制约了InSAR技术用于围填海区域的地面形变监测。哨兵-1A(Sentinel-1A)卫星是欧洲航天局(ESA)哥白尼计划(global monitoring for environment and security, GMES)的新一代对地观测卫星，于2014年4月3日发射，具有影像分辨率高、全球大范围监测、重访周期短和数据公开等优势[15]，可为深圳围填海区域的地面形变时序监测提供较好的数据支持。综合上述分析，现利用TS-InSAR技术和深圳围填海区域2016年12月6日—2019年10月4日的18景Sentinel-1A数据提取深圳围填海区域的地面形变结果，以此为基础重点分析深圳围填海区域的地面沉降规律及沉降原因。

1 研究区和数据源

深圳位于中国华南地区、珠江口东岸，西濒伶仃洋和珠江口，东临大鹏湾和大亚湾。改革开放后凭借着濒临香港、背靠珠三角的天然地理优势和其独有的经济特区政策，深圳经济迅猛发展，现已成为中国超一线城市中最年轻的城市以及世界一线城市。然而，随着城市规模的急剧扩张，深圳土地资源紧缺问题日益凸显。为了缓解土地供应紧张问题，深圳开始了较大规模的围填海，图1为深圳围填海区域概况。

图1 深圳围填海区域概况

基于TS-InSAR技术和深圳填海区2016年12月6日—2019年10月4日的18景Sentinel-1A数据，利用ESA发放的精密定轨(precise orbit determination，POD)去除轨道误差，利用航天飞机雷达地形测绘使命(shuttle radar topography mission，SRTM)数字高程模型(digital elevation model，DEM)Version4 90 m去除地形相位，提取了深圳围填海区域的地面沉降结果。表1为Sentinel-1A SAR数据参数。

表1 Sentinel-1A数据参数

2 TS-InSAR原理

首先，基于研究区的18景Sentinel-1A SAR影像构造时空基线，图2为试验的垂直基线分布图。

图2 垂直基线分布

根据生成的时空基线，从研究区(t0,t1,…,tN)时刻采集的N+1景SAR影像中筛选出主影像，将主影像和剩余影像配准并生成M景干涉图满足

(1)

对tA时刻采集的SAR影像和tB时刻采集的SAR影像(tA

δφj(x,r)=φ(tB,x,r)-φ(tA,x,r)≈

φnoise,j(x,r)+φatm,j(x,r)+

φdef,j(x,r)+φtopo,j(x,r)

(2)

式(2)中：x、r分别为像素点的方位向坐标、距离向坐标；φ(tB,x,r)、φ(tA,x,r)分别为tB、tA时刻采集的SAR影像上某一像素点(x,r)的相位；φnoise,j(x,r)为噪声相位；φatm,j(x,r)为大气相位的误差；φdef,j(x,r)为tA—tB时刻的卫星视线向(line of sight, LOS)形变相位；φtopo,j(x,r)为外部DEM引入的地形相位误差。

大气相位的误差计算公式为

φatm,j(x,r)=φatm,j(tB,x,r)-φatm,j(tA,x,r)

(3)

式(3)中：φatm,j(tB,x,r)、φatm,j(tA,x,r)分别为tB、tA时刻采集的SAR影像上的大气相位分量。

卫星视线向形变相位计算公式为

(4)

式(4)中：λ为波长；d(tB,x,r)、d(tA,x,r)分别为tB、tA时刻到t0时刻间卫星LOS累积形变量。

外部DEM引入的地形相位误差的误差可表示为

(5)

式(5)中：B⊥j为干涉图j的垂直基线；R为卫星与目标点间的距离；Δz为DEM误差；θ为视角(SAR波束和地面垂直直线形成的夹角)。

此外，式(2)干涉相位δφj(x,r)中的噪声相位分量、大气延迟相位误差分量和地形相位误差分量需要精确估算出并排除，则有

Aφ=δφ

(6)

式(6)中：A为M×N系数矩阵(M为干涉图数量，N为SAR影像数量)，且∀j=1,2,…,M，向量φT=[φ(t1),φ(t2),…,φ(tN)]由每景SAR影像上高相干点目标的相位值组成，向量δφT=[δφ(t1),δφ(t2),…,δφ(tN)]由各个差分干涉图的解缠相位组成。

用两景SAR影像之间的平均相位速率替换相位值以接解算影像上各高相干点目标的位移形变速率，式(6)可表示为

Bv=δφ

(7)

若B秩亏(M

3 结果与分析

3.1 形变结果与分析

利用TS-InSAR技术和深圳围填海区域2016年12月6日—2019年10月4日的18景Sentinel-1A数据提取的地面形变结果如图3所示，图中正值表示地面在垂直方向上隆起，负值表示地面在垂直方向上下沉。由图3的地面形变结果可以发现，2016—2019年，深圳西海岸的地面形变速率为-19.85～15.60 mm/a，深圳东部BA地区(宝安区)和西部FT地区(福田区)的地面沉降相对稳定，于NS地区(南山区)监测到大面积的地面沉降。

图3 深圳围填海区域的地面形变速率图

为检测沿海岸线的形变情况，从图3中沿AA′提取了地面形变剖面图，如图4所示。由图4可知，AA′沿线地面自西向东呈现逐渐下沉的趋势，AA′中最大沉降点A1、A2(分别靠近南山区与宝安区、福田区的边界)的垂直沉降分别为-13.50、-15.03 mm。可以发现，沿海岸线的沉降漏斗主要分布在南山区，则重点对南山区进行沉降机理的分析。

图4 沿A-A′提取的地面形变剖面图

3.2 精度分析

由于没有获取到外部水准数据进行辅助验证，因此，为了更好地评估Sentinel-1A InSAR形变反演的内部精度，通过计算速度的线性拟合偏差得到形变标准差，并对形变标准差进行统计分析，图5给出本次试验平均形变速率的标准差分布。若点目标表现出较强的非线性运动(即表现为高标准差值)，则相对于线性模型会产生较大的残差。由图5可以发现Sentinel-1A InSAR形变反演得出的平均形变速率的标准差为1.39 mm/a，点目标点最大标准差达到9.26 mm/a，96.75%的形变点目标的标准差小于2.25 mm/a。上述分析表明Sentinel-1A InSAR形变反演具有较高的精度和可靠性。

图5 形变速率标准差的分布

3.3 南山区沉降机理分析

从地理上看，NS(南山区)濒临南海，在亚热带海洋性季风气候的作用下，四季温润、年均降雨量多，再加上该地区平地、山丘、台地相间，地势北高南低，地面容易产生向向南倾斜下沉的趋势。发展上看，南山区不仅是深圳的科研、娱乐、教育以及体育中心，还是深圳新的经济增长极，自国务院2016年5月将深圳市南山区设为中国首批双创示范基地，越来越多的企业开始入驻。随着围填海区域经济建设的快速推进，容易加速诱发南山区的不均匀沉降。一般的，形变点目标临近区域具有相同的形变趋势。为分析南山区地面沉降的时间演化及其与降水的关系，从图3中选取形变点目标a、b、c、d与CW(赤湾)点监测的月平均降雨量进行对比，图6为结果分析图。由图6可以发现，所选形变点目标的沉降时间序列随降雨量的季节变化而非线性递减：平均降水量增加的时间里，可能是由于降水有效地补充了地下水，所选点目标的地面沉降速率较为明显地减缓了。从空间位置上看，靠近海岸的形变点目标明显受降雨的影响明显大于其他点，这可能与深圳的沉积层厚度有关。

图6 地面沉降的时间演化及其与降水的关系

图7为深圳围填海区域的沉积层厚度分布图。从时空演变的角度来看，深圳是珠江三角洲流水的长期沉积作用而汇集成的冲积平原。深圳围填海区域广泛分布着第四纪的海洋沉积物，这些沉积物形成的沉积层具有孔隙率大、含水量高、压缩性高等特点，因而沉积层厚度大的地面在施工建设和后期土体自固结过程更易发生形变。通过对比分析图7和图3，可以发现沉积层厚度为10～20 m的地区的地面沉降相比沉积层厚度为3～10 m的地区要明显。

图7 深圳围填海区域沉积层厚度分布

4 结论

利用TS-InSAR技术和深圳填海区2016年12月6日—2019年10月4日的18景Sentinel-1A数据，提取了深圳围填海区域的地面沉降结果，在此基础上重点分析了围填海区域的地面沉降规律和地面沉降原因。结果显示：深圳西海岸的地面形变速率为-19.85～15.60 mm/a，深圳东部宝安区和西部福田区的地面沉降相对稳定，于南山区监测到大面积的地面沉降；与年均降雨量的对比分析显示了南山区的地面沉降存在明显的季节变化规律，表明该地区的地面沉降与降水有关；靠近海岸的形变点目标明显受降雨的影响明显大于其他点，这可能与深圳的沉积层厚度有关。为了更好地探究围填海区域地面沉降，下一步，将采集填海区地下水数据和沉积层钻孔测试数据进行更进一步的研究。