SBAS-InSAR方法通过限定合适的时间、空间基线阈值, 将已有的SAR影像数据集划分为若干个小集合, 使得各个子集内时空基线较小, 从而克服时空失相干现象^[22], 并能充分利用SAR影像数据, 在降低解算复杂度的同时保持一定的观测精度。

假设在t₀…t_i…t_N时刻获取了N+1幅SLC SAR影像, 在时空基线阈值限定条件下将SAR影像分为L组, 之后将影像之间相互配对组成短基线子集, 每个短基线子集中的若干影像相互配对后差分干涉处理, L组影像共生成M幅差分干涉图, 规定N为奇数, 则M可以表示为

N+12≪M≪N+12N。 (1)

式中: N为SAR影像幅数, 幅; M为差分干涉图幅数, 幅。

以t₀为初始时刻, 某一像元(x, r)在 ti(i=1, …, N)相对于t₀时刻在视线向(Line Of Sight, LOS)的累积形变量d(t_i, x, r)(i=1, …, N)为待求量。假设t_A和t_B时刻(t_A早于t_B)的2幅SLC SAR影像配对生成第k (k=1, …, M)幅差分干涉图, 所得相位为 δφ(tk, x, r)(k=1, …, M)。

忽略大气效应和高程残余等噪声引起的相位误差, 则(x, r)处的相位为

δ(tk, x, r)=δφtB, x, r-δφtA, x, r≈4πλ[dtB, x, r-dtA, x, r]。 (2)

式中: δφtA, x, r为t_A时刻的影像相位; δφtB, x, r为t_B时刻影像相位; dtA, x, r为t_A时刻视线向形变量, mm; dtB, x, r为t_B时刻视线向形变量, mm。

这个含有M个等式的N元方程组可以用矩阵表示为

δφ=A[d]。 (3)

式中: A为邻接矩阵, 其类型由影像配对生成的干涉图形式决定。当小基线子集仅有一个时, L=1, 可知M≥ N, A是一个N阶矩阵, M=N时, 方程组精确可解; M> N时, 方程个数大于未知量个数, 为超定方程组, 用最小二乘法通常可以将φ 的估计值 φ^表示为

φ^=ATA-1ATδφ。 (4)

但是, 通常情况下影像数据并非仅有一个小基线子集(L> 1), 此时的 ATA就成为了一个奇异矩阵, A的秩为(N-L+1), 方程组有无限多解, 基于奇异值分解(singular value decomposition, SVD)方法即可求出方程组在最小范数意义上的最小二乘解。同时, 估计并去除残余地形相位和大气延迟相位, 即可得到地表形变的时间序列。

Sentinel-1 A是欧空局为“ 哥白尼计划” (前身是GMES计划)发射的首颗卫星, 由2颗 C 波段SAR构成。本次研究采用2017年5月20日至2021年12月1日间的Sentinel-1 A IW(interferometric wide swath)成像模式下的VV极化升轨数据(Ascending Path 142), 共计47景, 相邻两幅影像之间最短时间间隔均为36 d。此外, POD精密定轨星历数据(POD Precise Orbit Ephemerides)和SRTM1 DEM数字高程模型被用作计算轨道误差与模拟地形相位, 并从干涉相位中去除。影像参数信息如表1所示。

表1

Tab.1

表1(Tab.1)

表1 Sentinel-1 A雷达影像数据参数 Tab.1 Sentinenl-1 A data parameters 数据模式 时间范围 数据量 极化方式 入射角 升降轨 IW SLC 20170520— 20211201 47景 VV 38.97° 升轨

表1 Sentinel-1 A雷达影像数据参数 Tab.1 Sentinenl-1 A data parameters

首先对Sentinel-1 A数据进行裁剪, 后对覆盖研究区范围的图像进行SBAS处理(图2), 主要包括: 选取2017年11月16日的Sentinel-1A影像作为超级主影像, 其余所有从影像与主影像进行配准; 设定空间基线阈值130 m(2%), 时间基线阈值200 d, 数据多视比4:1(距离向:方位向), 共生成了220幅干涉图(图3); 基于SRTM1 DEM计算模拟地形相位对干涉图进行去平, 消除地形相位得到差分干涉图; 采用最小费用流(minimum cost flow, MCF)方法实现相位解缠, 去除相干性低和解缠结果不理想的相对; 创建地面控制点对数据对进行轨道精炼和重去平, 估算和去除残余的恒定相位和解缠后残存的相位斜坡; 基于SVD估算平均位移速率和残余地形相; 分别设定大气高通(365 d)和大气低通(1 200 m)进行时间域和空间域的滤波, 估算和去除大气延迟相位, 得到研究区在监测时间段内的地表时序形变; 对SBAS反演结果进行地理编码, 得到研究区形变监测结果。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 SBAS-InSAR技术流程Fig.2 Technical process of SBAS-InSAR

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 时空基线连接图Fig.3 Space-time Baseline Connection Diagram

通过对砀山县近5 a来的47期Sentinel-1 A IW SLC影像做SBAS反演处理, 获取了研究区2017— 2021年视线向的年平均地表形变速率。忽略水平位移的微弱影响, 将视线向监测结果除以cos β (β 为入射角)得到垂直方向形变速率(图4), 其中负值代表地面向远离卫星传感器方向运动(沉降), 而正值则代表地面向靠近卫星传感器方向移动(抬升)。

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 砀山县2017— 2021年平均垂直方向地表形变速率Fig.4 Average vertical surface deformation rate in Dangshan County from 2017 to 2021

砀山县地面在观测期间整体呈现出“ 北部轻微沉降, 南部稳定, 砀城镇局部沉降严重” 的特点。县城砀城镇南部以东城社区为中心, 半径5 km范围内均属严重沉降区域, 东城社区附近最大沉降速率达50 mm/a, 砀城火车站附近沉降速率约为30 mm/a; 高铁新区东部与经济开发区西部, 沉降速率约为40 mm/a, 这些区域的工程建设等人类活动比较活跃。官庄坝镇刘堤头村— 玄庙镇孙黑楼村— 周寨镇大孟庄村一线沉降速率也较高, 为25~30 mm/a, 可能与这一带存在第四系全新统淤泥质黏土、含淤泥质软土分布有关。综合来看, 以砀山县境内陇海铁路为界, 北部大片区域沉降速率基本为6~20 mm/a, 约占总形变面积的78.20%; 南部区域地表稳定, 零散分布微弱沉降。

以监测起始时间2017年11月16日的影像为参考基准计算其余时相的累计形变, 绘制得到砀山县2017— 2021年地表累计时序形变图(图5)。从图中可以看出, 监测期间内砀山县的沉降区域不断扩大, 累计沉降量逐年增大, 沉降区域主要集中在以砀城镇为中心的县城周围。砀城镇沉降中心从2018年初显形态, 此时沉降量为55 mm, 2019年之后逐渐明显, 至2019年下半年, 沉降范围基本确定, 即以县城为中心, 向外辐射半径5 km左右, 沉降面积约78.50 km²。随着时间的推移, 该区域沉降量日益增加, 到2021年12月形变中心沉降量达186 mm。此外还在县城与李庄镇之间监测到一条较明显的条带状沉降区, 其沉降进程与县城沉降中心基本一致, 陇海铁路、310国道正位于这条沉降带上。县城以北的官庄坝镇、玄庙镇, 在2017— 2021年期间, 累计沉降超过120 mm, 沉降速率约为30 mm/a。砀山县东北部分的葛集镇、唐寨镇沉降也较为明显, 葛集镇南部与唐寨镇北部相邻, 此处沉降量最大超过100 mm, 沉降速率约为20 mm/a。砀山县西南地区, 赵屯镇至关帝庙镇一带, 地表有轻微抬升, 推测可能与郑徐客运专线工程建设有关。

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 砀山县2017— 2021年地表累计时序形变量Fig.5 The cumulative displacement of Dangshan County from 2017 to 2021

为更直观地表达地表位移随时间的变化, 在严重沉降区、一般沉降区和稳定区分别选定Ⅰ 、Ⅱ 、Ⅲ 三点进行时序分析(三点位置见图4)。图6为Ⅰ 、Ⅱ 、Ⅲ 三点的时序累计形变曲线, 其中Ⅰ 点位于地面沉降漏斗中心, 沉降速率高达55.66 mm/a; Ⅱ 点靠近沉降中心, 沉降速率约为33.15 mm/a; Ⅲ 点处于整体稳定区。监测期间各点的沉降量与时间基本呈线性关系。

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图6 时间序列形变量Fig.6 Time series deformation

值得注意的是, Ⅰ 、Ⅱ 两点累计形变在2018年7— 8月出现小规模回弹, 可能与该时段的异常强降雨有关(图7)。受“ 温比亚” 台风的影响, 宿州市 8月16日开始降小到中雨, 17日普降大到暴雨, 18日普降大暴雨, 部分地区特大暴雨, 具有暴雨强度大、范围广等特点。18日, 暴雨中心区域位于砀山县全境, 砀山县1 d、2 d、3 d最大降雨量分别为216.80 mm、255 mm、264.40 mm, 均居有资料记录以来第一位, 重现期分别为 60 a、超百年、90 a。2018年砀山县地下水资源量比2017 年增长12.40%, 比2019年增长24%, 比多年平均值增长28.10%。

图7

Fig.7

	Figure Option ViewDownloadNew Window
	图7 2018年每月降水量统计(萧县气象站)Fig.7 Statistics of Monthly Precipitation in 2018 (Meteorological Station in Xiao County)

严重沉降区域分布在人口密集的城区, 且未采取相应防治措施, 沉降量级和范围仍在快速增大, 应尽快开展防治预警工作, 合理限制地下水开采, 开展第四系结构分布详细调查, 条件允许的情况下采取人工回灌措施。

前人大量研究证实, 过量开采地下水是导致地面沉降的外因, 而中等、高压缩性黏土层和承压含水层的压实则是地面沉降的内因^[23]。

砀山县城区有开采井20余眼, 井深多为150~300 m, 主要开采中、下更新统及新近系含水层。2016— 2020年, 砀山县深层地下水年平均开采量约2 268万 m³, 其中城区约占50%。据砀山县2008年水文调查资料, 县城区漏斗中心水位埋深约为44 m, 最大动水位埋深达61.20 m。与已发生地面沉降且地质背景相似的阜阳市类比^[24], 在降落漏斗区, 当区域地下水位埋深下降至20 m时, 即可引发地面沉降。

本文收集了砀山县12眼监测井2017— 2022年的年平均水位变化情况, 12眼监测井可按成井深度分为50 m、150 m和350 m(两眼200 m左右)3个层级, 结尾编号分别为A、B和C, 各监测井的历年水位变化情况如图8所示。DS01A、DS01B、DS01C这3眼井位于砀城镇李洼村, 紧邻300 mm厚度软土层, 2017— 2022年间, 3眼监测井地下水位不断下降, DS01A地下水位从42.06 m下降至45.96 m, 下降了3.9 m; DS01B地下水位从40.11 m下降至43.6 m, 下降了3.49 m; DS01C成井深度350.3 m, 6 a内地下水位下降了16.28 m, 周边村庄沉降速率为15~16 mm/a。210A、210B、210C位于玄庙镇花园村花园小学附近, 210A地下水位基本稳定, 210B和210C水位略有下降, 下降值3~4 m, 周边村庄基本稳定, 未见明显沉降。610A位于良梨镇派出所, 水位下降了3.04 m; 610B、610C位于良梨镇礼河集污水处理厂, 610B年均水位下降1.93 m, 610C缺失2017年监测数据, 后4 a年均水位抬升0.55 m, 周边地区基本稳定。614A、614B、614C位于关帝庙镇黄楼农场, 614A监测井缺失2018年和2019年数据, 其余4 a水位稳定; 614B水位有波动, 波动范围在2 m以内; 614C缺失2017、2018年监测数据, 后4 a年均水位下降了5.56 m。

图8

Fig.8

	Figure Option ViewDownloadNew Window
	图8 砀山县监测井年平均水位变化情况Fig.8 Variation of annual average water level of monitoring wells in Dangshan County

砀山县城区开采井过于集中, 而地表岩性为第四系全新统粉土, 地下松散沉积物(N+Q)厚度近500 m, 极易形成局部降落漏斗。砀山县的软土区主要分布在县城北部地区, 占全县总面积的40.71%, 地层岩性主要为第四系全新统淤泥质黏性土。县域北部2条软土埋深等值线显示软土埋深为10 m, 2条等值线之间正是官庄坝、玄庙镇和周寨镇(图9), 与黄河故道以北玄庙镇EW向沉降区空间位置基本吻合。可以推断, 地下水开采和软土变形是造成砀山县地面沉降的主要因素。

图9

Fig.9

	Figure Option ViewDownloadNew Window
	图9 砀山县软土埋深等值线Fig.9 Soft soil depth isoline in Dang Shan County