基于SBAS-InSAR技术监测廊坊市地面沉降

江培华,胡宝俊,唐 伟

(1. 华北科技学院建筑工程学院,北京 东燕郊 065201;2. 中国矿业大学(北京)地球科学与测绘工程学院,北京 100083)

0 引言

地面沉降是一种在自然或人为因素下导致地下松散岩层固结压缩,并引起地面一定范围内标高降低的地质现象[1,2]。对于这种前期难发现、影响范围广、治理难度大的地质灾害已成为影响我国经济社会可持续发展的重要因素。近年来,由于地下水过度开采而带来的地面沉降问题越来越引起人们的关注,华北平原更是此类地质灾害的受灾区,地表非均匀的沉降已经对人民的生产生活带来了严重的经济损失,特别是对城镇线状工程的破坏。为了避免此类地质灾害继续造成威胁,必须对地表沉降情况进行有效监测和防治。

目前,地面沉降的监测方法主要有分层标监测、高精度水准测量等。这些传统的监测手段不仅耗资巨大,而且效率不高,传统监测方法中水准测量复测时间跨度过大,而分层标和GPS高程测量空间分辨率低,均不利于局部地区形变场时空演化特征分析与研究。合成孔径雷达差分干涉测量(Differential Interferometric Synthetic Aperture Radar,D-InSAR)技术是近二十年发展起来的一种新型空间对地观测技术,可以利用重复观测的合成孔径雷达(Synthetic Aperture Radar,SAR)影像相位信息获取大范围、高精度、面状分布的地面沉降形变信息,且监测精度理论上可以达到mm级。同时,依据雷达卫星的过境频率,实现地面沉降的动态跟踪监测,对地质灾害的预防做出了突出的贡献。

然而,传统的D-InSAR技术由于其易受大气效应的影响,以及精度低的局限性,为此,基于传统D-InSAR技术的局限性研究分析,提出了永久散射体合成孔径雷达(Permanent Scatterer InSAR,PS-InSAR)[3]和小基线集合成孔径雷达干涉测量技术(Small Baseline Subset,SBAS-InSAR)[4],这两种时序InSAR技术通过提取散射信息稳定的点反演形变信息,克服时空失相关,并通过时空滤波方法改正大气延迟误差问题。

廊坊地区地下水的过度开采和城镇化建设发展迅速,蒋喆等利用2015—2017年Sentinel-1A数据基于SBAS-InSAR技术对廊坊市进行了监测分析[5],得出廊坊地区的广阳区南尖塔镇存在明显的地面沉降,最大沉降速率达60.4mm/a,累计最大沉降量达107.6mm。周吕等[6]获取2007—2010年覆盖廊坊地区18景Envisat ASAR影像数据,基于可同时获取点目标和分布式目标沉降信息的多时相InSAR(Multiple Temporal InSAR,MTInSAR)技术,并结合前人研究结果、地下水开采以及工业发展状况等研究分析了廊坊市的地面沉降时空特征。李海君等结合振幅离差指数与相干系数识别,采用2007—2010年覆盖廊坊北部地区的33景Envisat ASAR重轨影像数据使用PS-InSAR技术,选取PS点,借助GIS平台,对研究区进行沉降监测、沉降差异性特征与成因进行分析[7]。周旭等利用SABS算法发展而来的InSAR TS+AEM(InSAR time series with atmospheric model)时序分析模型,分析了24幅Envisat ASAR影像数据,研究表明利用时序InSAR技术进行城市地表沉降监测具有良好的精度和稳定性[8]。上述研究均只获得雷达视线方向(line-of-sight,LOS)的形变信息,并假设形变只发生在垂直向,将LOS形变简单投影为垂直方向,存在一定的误差。

廊坊市受城市工业、生活开采深层地下水的影响,形成了以城区为中心的地下水位漏斗,漏斗断面呈上宽下窄的“V”字形[9]。近年来,廊坊市陆续出台了地下水限采、禁采等管理办法,并关停了城区自备井,深层地下水位开采得到了有效控制,水位下降幅度减缓。

本文基于Sentinel-1卫星2017年1月至2020年12月升、降两轨数据,利用SBAS-InSAR方法获取了2个轨道的LOS形变,并将两者融合获取了垂直向和水平向形变。在此基础上,分析了廊坊市地面沉降的时空演变格局。

1 研究区和数据

1.1 研究区概况

廊坊市位于河北省中部偏东,由于1974年行政区划变动,北京和和天津地域相交,使廊坊形成北部三县(市)与境域主体不相连接的版图,其北起燕山南麓丘陵地区,和首都北京为邻,南接壤沧州,可抵黑龙港流域。西连保定市的涿州、雄县、高碑店,东交天津市武清区、宝坻区、蓟州区,地处京津冀城市群核心地带、环渤海湾腹地[10];廊坊市辖安次区、广阳区两个区、三河市、霸州市两个县级市以及香河县、永清县、固安县、文安县、大城县和大厂回族自治县六个县,总面积约6429km2,截至2021年末,常住人口为553.82万人。

廊坊市受地质构造的影响,大部处于凹陷地区,随着地壳下沉,地面逐渐被第四纪沉积物填平,致使新生界地层沉降厚度较大,全市地貌比较平缓单调,以平原为主,一般高程在2.5～30m之间,平均海拔13m左右。由于洪积、冲积作用和河流多次决口改道淤积,沉积物交错分布,加上风力及人为活动的影响,境内地貌差异性较大,缓岗、洼地、沙丘、小型冲积堆等遍布,全市地貌呈现大平小不平状态。

研究区属于华北平原中北部,廊坊市南部整个区域,介于东经116.06°～116.96°,北纬38.17°～39.65°之间,范围如图1中黑色粗线包含的区域。

图中黑色粗线为廊坊市行政区划。方框为Sentinel-1卫星的覆盖范围(红色方框为升轨,黄色方框为降轨)。

1.2 SAR数据

Sentinel-1卫星是欧洲航天局哥白尼计划(Global Monitoring for Environment and Security,GMES)中的地球观测卫星,由两颗载有C波段合成孔径雷达的卫星组成。本研究收集了2017年至2020年的升轨(P142F121、P142F126)和降轨(P47F461)两个轨道数据,其中升轨数据为2017年5月至2020年12月共102幅影像,覆盖范围如图1中红色框所示。降轨数据为2017年1月至2020年12月共97幅影像,覆盖范围如图1中黄色框所示。

本文在选取时空基线时,时间基线阈值选择50天,空间基线阈值选择100m,对于时空基线不相连的情况,特别是在降轨影像间,因此手动增加相关干涉对(图2(b)中红色线),干涉对组合时空基线如图2所示。最终,升轨数据共生成290个干涉对,降轨数据共产生248个干涉对。SBAS-InSAR算法重在应用短基线的差分干涉图集,集合内的影像对基线距小,而集合间的影像对基线距大,集合构成相位回归分析的序列,用均值相干系数作为相干目标识别的指标,再根据最小范数准则,利用奇异值分解算法逐个分离相位成分,实现每个相干目标形变序列的获取,减少空间和时间去相关影响,保证精度与可靠性[1]。

图2 升降轨数据干涉对组合时空基线

2 InSAR分析方法

2.1 SBAS-InSAR方法

SBAS方法为了在提高点目标密度的同时保持干涉图中的相干性,采用多主影像方式,通过组合时间和空间基线阈值内的干涉对,提取在一定时间内保持相干的分布式散射体目标。假设获取同一地区、时间分别为t0,t1,t2,…,tn的N+1幅SAR影像,这些影像可以生成的干涉图数量为M,且M满足:

(1)

其中第j幅干涉图为ta和tb两个时刻的SAR影像组合而成,第j幅干涉图在某一地面目标点处的干涉相位δφab表示为:

(2)

式中,λ为波长;da和db为ta和tb两个时刻相对起始时刻SAR卫星视线向LOS的累计形变量;φtopo为地形残差相位;φatm为大气延迟相位;φnoise为噪声相位。在估计并移除上述三个相位分量后,上式简化为:

(3)

为了获取时序形变信息,上式中的相位变化与时间成线性关系:

(4)

相应的,第j幅干涉图的相位值可以表示为:

(5)

进一步,上式用矩阵形式表示为:

Bv=δφ

(6)

式中,B矩阵大小为M×N,由于SBAS为多主影像策略,矩阵B容易存在秩亏,利用奇异值分解方法可以求出最小范数意义上的最小二乘解[4]。

2.2 升降轨融合方法

对升降轨数据进行融合分解可以获取地表水平东西向和垂直向形变分布信息[11]。如图3所示,展示了SAR侧视成像几何及三维坐标下的角度参数。

图3 SAR侧视成像几何及角度参数

InSAR干涉图中某个像素的视线向LOS的形变实际上是该像素点三维位移分量在视线向的投影之和,如下:

rLOS=uxcosφsinθ-uysinφsinθ-uzcosθ

(7)

式中,γLOS为视线向形变量;θ为雷达信号入射角,即视线向与垂直方向的夹角;φ为卫星航向角。地面的位移矢量u=[ux,uy,uz],ux为沿东西方向的形变分量,uy为沿南北方向的形变分量,uz为沿垂直方向的形变分量,可将视线向的沉降量r,由式(7)分解得

(8)

忽略南北方向分量,即令uy为零,且考虑升降轨两个视线向形变则列出以下两个方程,如下:

(9)

将式(9)列为矩阵形式:

r=Pu′

(10)

u′=[PTP]-1PTr

(11)

得到ux与uz,即东西方向与垂直方向的形变。

3 廊坊市地面沉降结果

3.1 升轨和降轨LOS向形变结果

经过SBAS解算后得到了研究区升、降两轨数据的LOS地表形变速率(见图4)。其中,形变速率的正值和负值分别表示朝向卫星方向移动和远离卫星方向移动。参考点位东经116.48°、北纬39.32°,远离地面沉降中心(见图4中黑色倒三角)。可以看出升轨和降轨得到的沉降空间分布基本一致,证明了InSAR方法提取形变速率的可行性。从图4中可以看出,研究区存在两个明显的沉降区域,如图4中椭圆实线包含区域所示,分别位于天津武清区-廊坊霸州市(图4右侧椭圆)和廊坊固安县-保定雄县(图4左侧椭圆),此外研究区中存在部分小局部区域的沉降,分别位于文安县中部和广阳区廊坊市中心。天津武清区-廊坊霸州市的沉降影响范围最大,沉降分布差异最明显,LOS向平均形变速率为-76.08～-1.1mm/a(表1)。廊坊固安县-保定雄县的沉降影响范围较小,LOS向平均形变速率范围在-59.6～2.6mm/a(表2)。相比于大范围沉降中心,研究区小局部区域沉降速率普遍较小,其中文安县内最大形变速率为-30mm/a。

表1 天津武清区-廊坊霸州市特征点LOS向速率

表2 廊坊固安县-保定雄县特征点LOS向速率

图4 LOS向平均形变速率

3.2 垂直和水平形变融合结果

InSAR获得的形变是地表三维形变的LOS向的一维投影,仅以LOS向作为观测结果难以反映真实的地面沉降漏斗。多角度InSAR观测数据融合通过组合来自两个或者多个独立InSAR观测几何的LOS测量值,可以分离垂直和水平分量。由于升、降轨数据提取的像素点不一致,分离前需对升降两轨数据进行重采样,以使二者具有相同的规则网格。LOS 测量值随入射角变化很大。详见2.2节中所述,当舍弃式 (8) 中的uy,即忽略LOS观测中的南北分量时,估计的东西方向和垂直方向分量会更加准确。这仅适用于南北分量与东西 分量具有相同数量级的地表变形现象。如果南北分量显著大于东西分量时(如地震和滑坡变形的情况下),则不应忽略南北分量。本文研究区属于典型地下水开采引起的形变,以垂直沉降为主,水平变形一般具有相似量级,且朝着沉降中心移动。因此,在融合升、降轨数据时,本文忽略南北方向形变,仅分析垂直向和东西向形变。图5展示了研究区垂直方向和东西方向的形变速率,垂直向沉降速率最大达到-100.2mm/a(图5(a)),东西向形变速率最大为-29.6mm/a(图5(b)),所以区域内的形变原因主要以垂直方向为主。廊坊固安县-保定雄县沉降区(图5左侧椭圆)没有明显的水平形变,而天津武清区-廊坊霸州市沉降区(图5右侧椭圆)水平变形较为显著,整体上向沉降中心移动。在图5中,对于东西向水平形变,负值表示向西移动,正值表示向东移动。

图5 研究区垂直方向和东西水平方向平均形变速率

4 分析与讨论

4.1 天津武清区—廊坊霸州市沉降

武清区和霸州市由于开采大量深层地下水用于工业和农业发展,早在1992—2000年就出现了沉降现象,沉降速率为30～60mm/a[12]。之后在2003—2004年两个区域的沉降速率进一步增大,其中武清区王庆坨镇和霸州市胜芳镇为两个区域的沉降中心,霸州市胜芳镇的沉降速率一度高达192mm/a。2003—2010年和2012—2014年武清区和霸州市的沉降影响范围分别以王庆坨镇和胜芳镇为中心进一步扩大,其中武清区王庆坨镇沉降速率从100mm/a增大为153mm/a[12]。如图6所示,2017—2020年,武清区和霸州市沉降区不仅连接在一起,沉降范围还进一步扩大至安次区南部、文安县东部、大成县北部、北辰区西部和静海区西北部,沉降影响面积达到约1470km2。天津武清区-廊坊霸州市主要包含两个显著的沉降中心,分别位于王庆坨镇在内的整个武清区南部和北辰区西部地区以及胜芳镇在内的霸州市中部地区。

为量化展示霸州市胜芳镇范围内的累积沉降量,选取横跨沉降中心的一条剖面线A-B(图6(a)中黑色直线),绘制了剖面线上地面沉降随时间的变化累积(图7)。A-B剖面上最大累积沉降量为187.2mm。

将2017—2020年以年为单位分别做累积沉降分析,如图8所示,发现年度累积沉降量在逐年减小,2017年累积最大沉降量89mm,2018年累积最大沉降量67mm,2019年累积最大沉降量64mm,2020年累积最大沉降量65mm,2018—2020年稳定控制在65mm/a。其中由于2020年受新冠疫情和降水偏丰影响,地下水供水量相比2019年减少0.9×108m3,降幅超23%,为2016年以来地下水供水量最大降幅[13],因此2020年累积沉降量略有所减缓。从侧面反映出2014年以来政府陆续出台的各项地下水限采、禁采政策以及南水北调工程的建成通水使研究区深层地下水明显回升,有效缓解了地面沉降。蒋喆等人通过对廊坊市沉降时空分布特征进行分析和提取监测,得到的平均沉降形变速率区间是-17.9～8.5mm/a,在沉降空间分布和幅度上与本文的结果基本一致[5]。

图8 天津武清区-廊坊霸州市沉降区2017-2020年逐年累积沉降量

此外,提取和了沉降区内6个点(图6中P1～P6点)的沉降时间序列,绘制于图9,可以看出这6个点均呈线性沉降趋势,最大沉降速率位于P6点,速率达到68.5mm/a。P1、P2点位于沉降区边缘,其累积沉降值最小,分别为195mm和200mm。P3、P4、P5、P6点累积沉降量分别为210mm、215mm、249mm、274mm。其中P2点在胜芳汽车站附近,P4点在胜芳火车站附近,P5点靠近胜芳镇政府,三者均在一定程度上受该沉降漏斗的影响。

4.2 廊坊固安县-保定雄县沉降

固安县是华北地热资源开发利用条件较好的地区。地热资源在上世纪七十年代石油地质勘查中被发现,有利的开采地段主要分布在牛驼镇地区,面积约120平方公里内。雄县地热资源丰富,境内60%以上的区域储藏着优质温泉,基岩热储面积320平方千米,占牛驼镇地热田总面积的50%。地热资源的开采成为两地交界处资源开采发展而引起的地面沉降的首要因素。此外,西侯留工业区、赵岗工业区、仇小王工业区、陈台工业区、米东工业区在内的17个工业园区在距离沉降区20km范围内,工业区的大量抽取水资源成为引起地面沉降的另一大因素。

如图10所示,该沉降区域东起固安县羊驼镇,贯穿固安县马庄镇直达保定雄县金三角经济开发区内的板西村小阳村,西至大营镇全长29.6km。水平方向变形不显著。

图10 廊坊固安-保定雄县沉降区垂直沉降和水平变形

本文获取了廊坊范围内的三个特征沉降点Q1、Q2、Q3(点的位置如图10(a)所示),绘制了它们的沉降时间序列,如图11所示。其中Q1和Q3在沉降中心边缘,累积沉降值Q1点最小(175.4mm),Q2点最大 (206.8mm)。沉降范围内有众多交通要道,荣乌高速、大广高速、霸州立交桥、京环线等交通线路都穿过沉降区。

图11 廊坊固安-保定雄县沉降区Q1～Q3点时序形变图

4.3 地表季节性形变分析

自然或人为地下水补给和排放造成的水头变化会引起季节性地表沉降和抬升,本文通过对研究区沉降产生原因分析主要以人为因素为主,忽略了自然因素。廊坊市地表水分布稀少,浅层地下水埋藏浅,主要接受大气降水补给,其次为侧向径流补给、河渠渗流补给、地表水灌溉和井灌回归补给[8]。以人工开采消耗为主,其次为蒸发及侧向排泄。深层地下水主要为侧向径流补给和少量越流补给,主要是人工开采产生的消耗。如图12所示,廊坊市降水季节分布不均,主要集中在夏季,6～8月三个月降水量一般可达全年总降水量的70%～80%。研究区地表水属海河流域系统的永定河水系,主要河流有永定河、龙河、凤河等,均为季节性河流。如图12所示,通过对廊坊固安-保定雄县区域特征点时序形变去线性化分析发现沉降存在周期性波动,周期1年。相对于降雨期延后三个月,周期性抬升在该年11月至次年1月,周期性下降在次年1月至4月,研究区内地下水漏斗区长系列动态变化为年际间有小的起伏,但总的趋势成波状形下降。

图12 2017—2020年廊坊固安-保定雄县形变去线性化和逐月降水量

本文主要获取2017—2020年覆盖研究区范围的Sentinel-1A数据,即使有升降轨数据互相验证,但缺少地面监测数据进行对比分析。后续还可获取GPS监测数据进而对比分析InSAR获取地表形变的可靠性。

由于廊坊市城区水源以深层地下水为主,本研究仅简单分析地表水流向,缺少地下水位相关数据,仅仅依赖地表水和大气降水无法完整的解释地表沉降产生原因,故将地下水概况加入沉降形成原因也必不可少。

地下水过度开采只是城市地面沉降的主要因素,该研究希望将地面荷载、城市地下建设等因素也考虑进去,以增加形变成因的分析研究。

5 结论

(1) 廊坊市地面不均匀沉降明显,存在天津武清-廊坊霸州市和廊坊固安县—保定雄县2 个明显沉降区。其中,天津武清—廊坊霸州市是最严重的沉降区,最大沉降速率达到了 100.2mm/a。此外,研究区内散布着许多小型沉降区。

(2) 基于升降轨融合技术,对LOS向平均沉降速率分别提取垂直向和东西向沉降速率,研究区内的形变主要来源于垂直向沉降,印证了地下水过度开采,造成深层水位下降,导致地面沉降增加。研究区内沉降漏斗的出现伴随着工业区而出现,沉降漏斗的空间分布与工业区和城市建设区存在较高的相关性,而沉降区周围交通发达,荣乌高速、大广高速、霸州立交桥以及胜芳火车站等交通相关部门应引起重视。

(3) 通过对研究区内两个明显沉降区以及年降水数据进行时序分析,廊坊固安县—保定雄县沉降区的沉降发育具有1年周期的季节性变动。