基于SBAS-InSAR技术的安徽省砀山县地面沉降监测

2022-11-04牛地吴倩朱成林

中国地质调查 2022年5期

牛地，吴倩，朱成林

(1.安徽省地质矿产勘查局332地质队，安徽黄山 245000；2.安徽省地质实验研究所，安徽合肥 230041)

0 引言

地面沉降是一种主要由人类工程活动引起的地壳表面标高缓慢降低的环境地质现象，其变化过程缓慢，难以被及时察觉，且不易治理[1]。该现象自二十世纪20年代开始在我国东部沿海出现，现已广泛分布于100多个地级市[2]，严重威胁社会安全与经济健康发展，发展科学有效的监测手段是预防和减轻灾害损失的关键所在。

相较于传统的地面沉降监测方法，合成孔径雷达差分干涉测量(differential interferometric synthetic aperture radar，D-InSAR)技术具有全天候、全天时、高精度以及快速获取大范围地表高程和形变信息的优势[3]，但时空失相干和大气效应严重阻碍了其获取真实的地面形变信息。在此基础上发展的基于相干目标的时序InSAR处理技术，主要包括永久散射体(permanent scatterer-interferometric synthetic aperture radar，PS-InSAR)方法和小基线集(small baseline subsets-interferometric synthetic aperture radar，SBAS-InSAR)方法，能在很大程度上减少D-InSAR方法中失相干和大气异常引起的相位误差。SBAS-InSAR技术已广泛应用于地面沉降[4-5]、滑坡[6-8]、泥石流[9]等地质灾害的监测，其中于地面沉降方面的监测成果尤为丰富，致灾机理研究多指向地下水[10-15]和矿产开采[16-18]，并在地质构造活动[19]、工程建设[20-21]的形变监测中取得了良好的反演结果。上述研究都表明了小基线集方法在地表形变监测预警与防治方面有着显著的技术的优势。

本文基于SBAS-InSAR技术和Sentinel-1 A升轨数据，对宿州市砀山县2017—2021年间地表形变进行反演，得到了监测期间该区平均沉降速率和累积地面沉降量，并对地面沉降的时空分布和变化特征进行了细致的分析。在此基础上结合相关水文地质资料对该区地面沉降的主要诱发因素进行探究，以期通过对砀山县地面沉降的监测，为灾害预警提供决策参考依据。

1 研究区概况

砀山县地处淮北平原，位于安徽省最北部，地势开阔平坦，地貌特征以故黄河冲积平原为主，总面积1 193 km2，下辖13个镇。县内交通路网发达，陇海铁路和310国道贯穿全境，101省道和223省道及县乡道路构成以县城为中心的密集交通网络。县境内地层属华北地层大区晋冀鲁豫地层区徐淮地层分区淮北地层小区，据古近系、新近系和第四系厚度可知，古近纪区内南部地层相对上升，北部则相对下降接受沉积，沉积厚度达数百米；新近纪开始，全区主要表现为振荡下降，接受沉积，新近纪末期地壳上升，普遍形成沉积间断；第四纪以来以下降为主，接受沉积。砀山县地下水类型为松散岩类孔隙水，浅层孔隙水底板埋深50 m，深层孔隙水砂层顶板埋深165～195 m，城区地下水开采以深层为主。在深层地下水大量开采的漏斗区，人工开采成为深层地下水的主要排泄形式，同时地下水由漏斗外缘向漏斗中心径流。在砀山县北部及良梨镇东北区域，有第四系全新统淤泥质黏土、含淤泥质黏土(软土)分布，一般埋深4 ～10 m，最大埋深18.30 m，厚度小于15.56 m，对于上部建筑构成一定潜在威胁。

图1 研究区范围示意图(底图源自Google Earth)Fig.1 Schematic diagram of the study area(Base map from Google Earth)

2 研究方法与数据处理

2.1 短基线集(SBAS-InSAR)方法

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

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

(1)

式中：N为SAR影像幅数，幅；M为差分干涉图幅数，幅。

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

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

(2)

式中：δφ(tA,x,r)为tA时刻的影像相位；δφ(tB,x,r)为tB时刻影像相位；d(tA,x,r)为tA时刻视线向形变量，mm；d(tB,x,r)为tB时刻视线向形变量，mm。

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

δφ=A[d]。

(3)

(4)

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

2.2 数据简介

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 Sentinel-1 A雷达影像数据参数Tab.1 Sentinenl-1 A data parameters

2.3 SBAS-InSAR处理流程

首先对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 SBAS-InSAR技术流程Fig.2 Technical process of SBAS-InSAR

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

3 结果分析

3.1 年均形变速率

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

图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%；南部区域地表稳定，零散分布微弱沉降。

3.2 地面沉降场时空分布特征

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

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

3.3 时间序列分析

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

图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 2018年每月降水量统计(萧县气象站)Fig.7 Statistics of Monthly Precipitation in 2018(Meteorological Station in Xiao County)

3.4 沉降机理分析

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

砀山县城区有开采井20余眼，井深多为150～300 m，主要开采中、下更新统及新近系含水层。2016—2020年，砀山县深层地下水年平均开采量约2 268万 m3，其中城区约占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。

(a)DS01A 成井深度50.2m (b)210A 成井深度46.19m (c)614A 成井深度48.56m

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