利用时序InSAR技术反演邯郸平原区地表形变与含水层参数

白林，李振洪，宋莎，刘东

1 长安大学地质工程与测绘学院，西安 710054 2 西部矿产资源与地质工程教育部重点实验室，西安 710054 3 长安大学地学与卫星大数据研究中心，西安 710054 4 机械工业勘察设计研究院有限公司，西安 710043

0 引言

从承压含水层抽取地下水时，会引起水头下降，孔隙水压力减小，导致有效应力增大，含水系统骨架发生压缩(即地面沉降)，造成储水能力损失(Jiang et al., 2018).作为我国水资源最为短缺和环境问题最为突出的地区之一，华北平原地下水供水量达总供水量的近70%(邵景力等，2009).20世纪70年代以来，华北平原地下水经历了长期的超采过程，已形成二十多个大型复合水位降落漏斗(张兆吉等, 2009; 葛大庆等，2014)，导致了大面积的快速地面沉降，成为我国沉降面积最大、类型最复杂的地区之一.不均匀地面沉降不仅会对地下和地面建筑物、构筑物产生破坏，威胁高速、高铁等重大基础设施安全运营，而且会导致地面塌陷和地裂缝等地质灾害，给人们的生产、生活带来巨大的危害及损失(Bai et al., 2016).同时，地面沉降还会导致滨海平原地面高程资源损失，与全球海平面上升产生累加效应，加剧了风暴潮、海岸侵蚀、海水入侵等生态环境灾害的致灾程度(Xue et al., 2005).随着华北平原地面沉降灾害的不断发展，已极大地影响和制约着当地经济建设可持续发展及社会安定.因此，为了预防和治理华北平原地面沉降地质灾害，亟需查明地面沉降的时空演变特征，厘清承压地下水变化与地面沉降的耦合关系.

GPS、精密水准测量能够获取高精度的地面沉降监测结果，但因其需要现场布设且监测点稀疏，难以探测区域尺度地面沉降的时空分布特征.作为近年来发展迅速的空间大地测量新技术，时间序列InSAR(Multi Temporal Interferometric Synthetic Aperture Radar，MT-InSAR)技术具有全天时、全天候、精度高、覆盖范围大、时间/空间分辨率高等优点(Hooper et al., 2012; Bai et al., 2016)，已广泛应用于华北平原地面沉降监测(葛大庆, 2013; Yang et al., 2018).

目前已有相关学者成功将时序InSAR应用于含水层参数反演(Hoffmann et al., 2001; Chaussard et al., 2014).Hoffmann等(2001)首次提出了基于InSAR观测反演含水层参数的方法，估算了Las Vegas Valley地区测井附近的弹性骨架释水系数，结果与实测数据较为一致.当水头低于前期固结水头时，弱透水层往往存在延迟排水现象，若不考虑其影响，会导致非弹性骨架释水系数低估.因此，Hoffmann等(2003)提出了顾及延迟排水的含水层参数反演算法，估算了Antelope Valley地区的时间常数及非弹性骨架释水系数.随着时序InSAR技术的发展和SAR数据的积累，获取某一地区形变时间序列成为可能，大大拓展了InSAR在地下水相关研究的应用领域.许文斌等(2012)利用InSAR和地下水水位数据估算了洛杉矶地区的地表形变及含水层参数.Chaussard等(2014, 2017)分别根据1992—2011年多源SAR数据和2011—2017年COSMO-SkyMed数据估算了Santa Clara Valley地区的弹性骨架释水系数和弹性水储量变化.Reeves等(2011, 2014)反演了San Luis Valley农业区的弹性骨架释水系数，并用其预测了水头变化.Chen等(2016)使用ALOS数据估算了San Luis Valley地区的弹性骨架释水系数和迟滞时间.近年来，国内学者利用时序InSAR技术在华北平原开展了一些系列研究，但这些研究多集中于北京(Guo et al., 2020; Yu et al., 2020)、沧州(Jiang et al., 2018; Bai et al., 2022)等地区.华北平原水文地质条件复杂，地表形变对水头变化的响应差异显著(何庆成等，2006)，因此开展更大范围的地面沉降及含水层参数反演研究具有重要的意义.

本文将以华北平原承压地下水超采、地面沉降灾害严重的邯郸平原区为研究区，利用2015—2019年Sentinel-1A数据监测地表形变的时空分布特征，并结合水头数据分析含水系统对水头变化的不同响应；根据承压地下水变化与地表形变的耦合关系，估算邯郸平原区承压含水系统的弹性/非弹性骨架释水系数及时间常数.

1 研究区与数据

1.1 研究区概况

邯郸位于河北省南部，太行山脉东麓，地理坐标为东经113°28′—115°28′，北纬36°04′—37°01′，总面积约1.2万km2(如图1所示).邯郸市地势自西向东呈阶梯状下降，高差悬殊，地貌类型复杂多样.以京广铁路为界，西部为中、低山丘陵地貌，东部为华北平原.邯郸平原区约占全市面积的54%，按成因和地貌形态特征可分为山前冲积洪积平原和中东部冲积湖积平原(吴旭等，2012).平原区第四系含水层组厚度为350～500 m，可分为四个含水层组(第I～IV含水层组)，第Ⅰ、Ⅱ含水层组普遍分布有咸水、微咸水，淡水主要赋存于第II含水层组的下段及第III、IV含水层组，地下水主要开采层组为第II含水层组的下段及第III含水层组(石立明，2020).平原区地下水开采始于20世纪70年代末，随着区域经济迅猛发展，开采量不断扩大，水位逐年迅速下降，至2010年中心水位下降总量超70 m，形成了面积达3063 km2的深层承压地下水降落漏斗，并与衡水、沧州复合大型漏斗相连(付丹平，2013).承压地下水长期持续超采，导致了严重的地面沉降，至2015年邯郸平原区沉降中心累计沉降量已超1500 mm，沉降量大于500 mm的面积为3576 km2，占平原区总面积的47.5%(河北省国土资源厅, 2016).为了保护地下水资源，控制地面沉降灾害，邯郸市自2014年被列为全省地下水超采综合治理首批试点市以来，实施了一系列的地下水压采措施，并增加了地表水供给，使得区内地下水超采得到了有效遏制，截至2018年年底，累计压减地下水超采量4.02亿m3(石立明，2020).

图1 邯郸市地理位置白色方框表示水头观测井位置，蓝色矩形表示所用Sentinel-1A数据覆盖范围，绿色实线表示含砂率，棕色实线表示断裂.Fig.1 The location of HandanThe white boxes represent the locations of observation wells, the blue rectangle marks the coverage of Sentinel-1A data, the green lines represent the sand content, and the brown lines show the faults.

1.2 实验数据

本文所用SAR数据为2015年7月30日—2019年12月12日获取的117景Sentinel-1A升轨数据，所处波段为C波段，极化方式为VV极化，方位向和距离向分辨率分别为20 m和5 m.原始SAR数据幅宽为250 km，选取了其中170 km×130 km的感兴趣区，覆盖范围如图1所示.此外，本文使用覆盖实验区的30 m分辨率SRTM DEM数据去除干涉图中地形相位并进行地理编码.

本文收集了邯郸平原区6个承压地下水水头观测井2010年至2017年的观测数据(空间分布如图1所示)，用于分析地面沉降与水头变化的响应关系和计算含水层参数.本文未收集到井位的观测层位及深度信息，但观测井最初用于地下水开采，因此观测井记录的是第II含水层组的下段及第III含水层组的水头变化.

2 数据处理方法

2.1 时序InSAR数据处理

邯郸平原区以农田为主，人工目标点较为稀少，此外，华北平原地下水超采引起的地表形变往往伴随着明显的季节性波动(葛大庆等，2014)，因此本文选择适用于非城市区域，且具有反演地表周期性过程能力的StaMPS方法(Hooper et al., 2004, 2007a)进行邯郸地区地表形变监测.本节将对该方法进行简要介绍.

首先，为了抑制由于时间/空间基线过长引起的失相干，设置时间基线阈值为36 d，空间基线阈值为150 m，并根据生成的干涉图质量及相干性进行挑选，最终获取340个干涉对用于时间序列分析，具体干涉对组合如图2所示.为进一步提高相干性，分别在方位向去除不重叠的多普勒频率，在距离向进行频域滤波(Hooper, 2008).然后，利用轨道数据和SRTM DEM数据去除干涉图中地形相位，获得差分干涉图，并对其进行地理编码.

图2 Sentinel-1A干涉对时空基线图Fig.2 Temporal/spatial baselines of Sentinel-1A interferometric pairs

StaMPS方法是基于滤波相位失相干缓慢(Slowly-Decorrelating Filtered Phase，SDFP)的目标点进行时间序列分析的(Hooper, 2008).为提高计算效率，首先对所有像元进行幅度稳定性分析，设置幅度差离散度阈值为0.6选取SDFP候选点.然后，对获取的SDFP候选点进行相位稳定性分析，剔除不满足相位稳定性标准的候选点，获取SDFP点集(Hooper, 2008).利用三维解缠算法进行相位解缠(Hooper et al., 2007b)，获取SDFP点的解缠相位.为获取形变相位，需从解缠相位中消除DEM误差、轨道误差、大气延迟等干扰项的影响.本文根据DEM误差与垂直基线的线性关系反演估算并去除其影响(Hooper et al., 2007a).因为邯郸平原地区地势平坦(高程范围为20～60 m)，所以大气延迟中垂直分层部分对形变监测的影响较小(Liu et al., 2016)，通过估算并去除每个干涉图的相位趋势面以消除轨道误差及大气延迟长波段部分的影响(Hooper et al., 2007a).最后，根据形变相位与形变量的线性关系即可反演获取每个SDFP点视线向形变量时间序列及平均形变速率.本文未对大气延迟中湍流分量进行改正，通过计算山前平原稳定区形变时间序列的标准差来评估其对形变结果的影响，结果表明该区域形变时间序列标准差的平均值为1.5 cm，因湍流分量具有明显的随机性(杨成生等, 2011)，其对形变速率结果影响较小.邯郸地区沉降主要由承压地下水超采引起，与垂直形变相比，地下水超采引起的水平形变相对较小(Jiang et al., 2018; Yang et al., 2018)，因此本文假设水平形变可忽略，那么根据InSAR形变监测结果及SAR影像入射角即可估算邯郸地区垂直形变.

2.2 骨架释水系数估算方法

承压含水系统由含水层(砂、砾石)和弱透水层(黏土、淤泥)组成.骨架释水系数是表征含水系统释水能力的参数，指水头下降(升高)一个单位时，单位面积含水系统由于骨架压缩(膨胀)而释出(储存)的水的体积.承压含水系统骨架释水系数可根据含水系统厚度变化量与水头变化的线性关系估算(Riley, 1969)：

Sk=Sskb0=Δb/Δh,

(1)

式中，Sk为承压含水系统骨架释水系数，Ssk为承压含水系统骨架释水率，b0为承压含水系统原始厚度，Δb为地表形变量，假设其等于承压含水系统厚度变化，Δh为承压含水系统水头变化.

骨架释水系数在有效应力超过历史最大有效应力前后差别很大，这个历史最大有效应力称为前期固结应力，对应的水头为前期固结水头(Leake, 1990).有效应力小于前期固结应力时，若有效应力增加(或水头降低)，则导致含水层和弱透水层发生弹性压缩，当有效应力恢复到初值时，压缩可恢复.如果有效应力超过前期固结应力，弱透水层将发生非弹性压缩，这种非弹性压缩是沉积物颗粒的移动和重新排列导致的，很大程度上是永久性的.相同有效应力变化条件下，弱透水层的非弹性压缩远远大于弹性压缩，前者可以比后者大一到两个量级(Hoffmann et al., 2001).含水层中砂、砾石强度较高一般不发生非弹性压缩(Smith et al., 2017).

为了刻画有效应力超过前期固结应力前后骨架释水系数的显著变化，经常用到以下两个参数：

(2)

式中，Ske为弹性骨架释水系数，Skv为非弹性骨架释水系数，σe为有效应力，σe(max)为前期固结应力.非弹性骨架释水系数通常比弹性骨架释水系数大一到两个量级(Hoffmann et al., 2003).

前期固结水头难以直接测量，且会随着地下水开采而发生变化，在实际工作中往往难以确定地表形变是否为弹性，这给弹性/非弹性骨架释水系数计算带来了困难.Miller和Shirzaei(2015)提出了基于季节性形变及水头变化估算弹性骨架释水系数的方法.这一方法的基本假设是季节性形变是弹性形变的一部分，是由水头的季节性变化引起的.则根据形变与水头变化的线性关系(式(1))，弹性骨架释水系数可用(3)式计算：

Ske=Δbs/Δhs,

(3)

式中，Δbs为含水系统厚度季节性变化，即季节性形变，Δhs为水头季节性变化.由于含水层不发生非弹性压缩，因此无论水头是否高于前期固结水头，含水层骨架压缩都是弹性的.当水头高于前期固结水头时，含水层和弱透水层骨架压缩都是弹性的，因此式(3)计算的是含水层和弱透水层弹性骨架释水系数之和.当水头低于前期固结水头时，只有含水层骨架压缩为弹性，因此式(3)计算的仅为含水层弹性骨架释水系数.

由于弱透水层由细颗粒沉积物组成，垂向渗透系数较低，含水系统中弱透水层的水头变化通常滞后于周围含水层(Hoffmann et al., 2003).因此弱透水层压缩也滞后于水头变化，滞后的时间称为时间常数，水头下降引起的93%压缩量是在该时间段内发生的.非弹性骨架释水系数及时间常数可用(4)式计算：

(4)

式中，τ0为时间常数，Δbinelastic为非弹性形变，可用(5)式计算：

Δbinelastic=Δb-SkeΔh.

(5)

2.3 延迟沉降分析

当水头停止下降时，由于弱透水层的渗透系数较低，弱透水层延迟排水导致的延迟沉降仍将继续，可用(6)式表示：

Δb(t)=M(ekt-1),

(6)

式中，M为幅度；k为衰减系数，取值范围为-1～0，表征含水系统对水头变化的响应速度.k越小，沉降速率减缓越快；k越大，沉降速率减缓越慢；如果k趋近于0则表明沉降无减缓趋势，地面沉降不是由弱透水层延迟排水引起的.

2.4 季节性信号分离

为估算弹性骨架释水系数，本文使用谐波函数提取季节性形变及季节性水头变化(Neely et al., 2021)：

Y(t)=vt+Acos(2π(t-T))+B,

(7)

式中，Y(t)为形变/水头变化时间序列，v为形变/水头线性变化速率，A为季节性变化幅度，t为时间，T为季节性信号达到峰值的时刻，B为常数.

2.5 水头变化与形变的相关性分析

水头的季节性变化会引起明显的季节性形变，由于含水系统存在低渗透介质，形变往往滞后于水头变化.利用(8)式对时间序列进行相关性分析，可估算形变滞后水头变化的时间：

(8)

式中，τlag为形变滞后水头变化的时间，corr表示相关系数，Δhdetrend为去除线性趋势后的残余水头时间序列，Δbdetrend为去除线性趋势后的残余形变时间序列.在进行相关性分析时，需将残余形变时间序列线性插值从而获取与水头观测时间对应的形变结果.

3 结果与讨论

3.1 地表形变监测结果与解译

2015—2019年邯郸平原区地表形变监测结果如图3所示，地面沉降分布广泛，西部受承压含水层分布控制，东部受馆陶西断裂及元村断裂控制，北部与邢台沉降区相连.邯郸平原区可分为两个沉降区：北部沉降区(图3中区域S1，以曲周县、肥乡县、广平县等为中心的有咸水区)和南部沉降区(图3中区域S2，以临漳县为中心)，两个沉降区已连为一片，最大沉降速率可达14 cm·a-1.地面沉降主要是由承压地下水超采导致的含水系统骨架压缩引起的，为进一步分析含水系统对水头变化的不同响应，本文比较了2个典型观测井水头变化与其附近SDFP点的形变时间序列，结果如图4所示.D3观测井2015—2017年水头呈上升趋势，但地面沉降仍在继续(图4a)，表明该区域地面沉降主要是由于弱透水层延迟排水引起的.受逐年升高的水头影响，沉降速率快速减缓，衰减系数达-0.74，到2017年末沉降速率已趋近于0.D4观测井2015—2017年水头变化趋势和沉降趋势大体上是一致的，地面沉降速率无明显的减缓趋势(图4b)，表明该沉降主要是由2015—2017年水头下降引起.相较于2015—2017年，2个典型观测井附近区域2018—2019年沉降速率明显加快，可能的原因是该区域承压地下水发生了快速下降.两个观测井含水系统对水头变化的不同响应表明邯郸平原区含水层和弱透水层的空间分布及厚度差异显著，这一现象和沧州地区较为相似(Jiang et al., 2018).由于部分地区地面沉降发展趋势在2018年前后发生了明显变化，本文分别计算了2015—2017年和2018—2019年的邯郸平原区地面沉降的衰减系数.由图5a可知，2015—2017年西部山前平原、魏县东部及馆陶县东部衰减系数在-1～-0.4之间，尤其是邯郸市区及磁县、永年县、魏县、馆陶县城区衰减系数可达-1，地面沉降存在明显的减缓趋势；但在中部沉降区衰减系数均大于-0.1，地面沉降无明显的减缓趋势.2018—2019年邯郸平原区沉降衰减系数明显增大，除西部山前平原部分地区、魏县东部及馆陶县城区外，衰减系数均大于-0.1(图5b)，地面沉降无明显的减缓趋势，其空间分布在东部与沉降分布大致一致.造成这一现象可能的原因是2019年降水量仅为多年平均量的87%，导致承压地下水开采量增加，地面沉降加剧.以上结果表明地下水压采措施在部分城市地区取得了较好的效果，但在地面沉降严重的中部地区效果不显著，主要是因为该区域地表水资源匮乏，难以满足该地区的用水需求.

图3 邯郸平原区2015—2019年平均形变速率图白色十字表示参考点位置，白色实线表示断裂，白色方框表示2个典型观测井的位置，灰色实线表示承压含水层边界.Fig.3 The mean deformation rate map across Handan from 2015 to 2019The white cross represents the location of the reference point. The white lines represent the faults. The white boxes represent the 2 typical observation wells. The grey line marks the boundary of confined aquifer system.

3.2 季节性形变结果

承压水水头的季节性变化导致形变存在明显的季节性波动(图4c和4d)，且形变滞后于水头变化约45 d(图6)，这主要是由于该含水系统垂直渗透系数较低，被开采含水层水头需较长时间才能与相邻含水层或弱透水层达到平衡.基于谐波函数提取了典型井位的季节性形变(2015—2016年D3观测井无明显季节性水头变化及季节性形变，因此仅提取了2017年季节性形变)，其幅度为7～15 mm，形变峰值对应时间为每年的2—4月.从图4中可看出，2015—2017年及2018—2019年两个时段的季节性形变幅度有明显变化，可能的原因是降水变化导致的地下水开采强度变化引起的.为分析邯郸平原区季节性形变的时空分布特征，分别计算了2015—2017年及2018—2019年的季节性形变幅度及峰值时间(图7).结果表明，2015—2017年邯郸平原区北部沉降区及东南部区域存在明显的季节性形变，北部沉降区幅度为5～20 mm，东南部区域幅度为10～30 mm(图7a).受降雨减少的影响，2018—2019年沉降区域地下水超采加剧，引起的水头季节性变化幅度增加导致了沉降区季节性形变幅度增大，最大幅度可达25mm(图7b).东南部区域地面沉降呈现减缓趋势(图5b)，表明该区域深层承压地下水开采强度有所降低，导致季节性形变幅度减小至5～15 mm(图7b).自然状态下，地下水受降水影响，水头在每年的8—9月份达到最高值.但邯郸平原区地下水变化主要受农业灌溉控制，水头在每年的1—3月份达到最高值，6—7月份达到最低值(付丹平，2013).北部沉降区浅层地下水(潜水)为咸水，农业灌溉主要依赖承压地下水，因此2015—2017年该区域季节性形变峰值时间为1—3月(图7c).其他区域峰值时间主要为8—9月，与自然状态下承压地下水最高水头时间相符，表明其灌溉主要依赖浅层地下水，深层承压地下水受浅层地下水越流补给及侧向径流补给影响，在8—9月达到最高水头.受降雨减少的影响，2018—2019年承压地下水超采范围增大，季节性形变峰值时间为1—3月的范围扩展到整个沉降区域及山前平原北部(图7d)，表明这些区域的深层地下水开采量增加，引起了明显的季节性形变.

图4 (a)观测井D3，(b)观测井D4水头时间序列与其附近SDFP点形变时间序列对比图.(c)观测井D3，(d)观测井D4残余水头时间序列与残余形变时间序列对比结果Fig.4 Head time series vs. land deformation time series for nearby SDFP pixels at (a) D3 and (b) D4. Detrend head time series vs. detrend land deformation time series for nearby pixels at (c) D3 and (d) D4

图5 (a)2015—2017年和(b)2018—2019年邯郸平原区地面沉降衰减系数白色实线表示断层，紫色实线表示沉降速率60 mm·a-1等值线.Fig.5 The map of the decay coefficient across Handan during (a) 2015—2017 and (b) 2018—2019The white lines represent the faults, and the purple lines represent the contours of subsidence rate of 60 mm·a-1.

图6 去除线性趋势后的残余形变与残余水头时间序列相关性分析Fig.6 Correlation analysis between the detrend deformation and the detrend head time series

图7 (a)2015—2017年和(b)2018—2019年季节性形变幅度.(c)2015—2017年和(d)2018—2019年季节性形变峰值时间紫色实线表示沉降速率60 mm·a-1等值线.Fig.7 The amplitude of seasonal deformation during (a) 2015—2017 and (b) 2018—2019. The peak time of seasonal deformation during (c) 2015—2017 and (d) 2018—2019The purple lines represent the contours of subsidence rate of 60 mm·a-1.

3.3 骨架释水系数估算结果

利用季节性形变及季节性水头变化估算了邯郸平原区弹性骨架释水系数.由于观测井D2的水头没有明显的季节性变化，未用于弹性骨架释水系数估算，因此本文最终估算了5个观测井位的弹性骨架释水系数(表1).邯郸平原区弹性骨架释水系数估算结果介于1.51×10-3～4.05×10-3之间，与抽水试验结果1.0×10-3～8.0×10-3(张兆吉等，2009)较为相符.弹性骨架释水系数由北向南大体上呈增大的趋势(图8)，可能的原因是北部地区浅层地下水为咸水，为满足农业用水需求大量开采深层承压地下水，导致该区域开采厚度较大(式(1)中较大的b0).

当有效应力超过前期固结应力(水头低于前期固结水头)时，表征单位面积含水系统释水能力的参数称为非弹性骨架释水系数.非弹性骨架释水系数通常比弹性骨架释水系数大一到两个量级(Riley, 1969).本节采用2.2节所述方法估算非弹性骨架释水系数和时间常数，需挑选水头低于前期固结水头且仍在下降，并引起明显地面沉降的观测井位.结合2000—2017年水头观测资料及InSAR地面沉降监测结果，本文选取了2个满足该条件的观测井.非弹性骨架释水系数估算结果为3.62×10-2～4.57×10-2(表1).在观测井D4，非弹性骨架释水系数约为弹性骨架释水系数的11.75倍，表明当水头低于前期固结水头时，地表形变以非弹性形变(永久性沉降)为主，可恢复的弹性形变占比仅为8%左右.观测井D2和D4时间常数分别为0.47 a和0.77 a，而由图4a可知，观测井D3的时间常数应大于2 a(水头上升期：2016年1月—2017年12月)，但该观测井距离D2和D4观测井较远(图8)，水文地质条件可能存在较大差异，导致时间常数相差较大.此外，邯郸地区时间常数明显小于华北平原东部(6 a)(Bai et al., 2022)，表明该区域弱透水层渗透性较好，对水头变化响应速度较快.

表1 邯郸平原区骨架释水系数估算结果Table 1 The results of the estimated skeletal storativity across the Handan plain

图8 邯郸平原区弹性骨架释水系数插值结果观测井位置如图中白色方框.Fig.8 Interpolation map of elastic skeletal storativity across the Handan plainThe white squares represent the well locations.

4 结论

本文利用2015—2019年117景Sentinel-1A数据对邯郸平原区地表形变进行监测，获取了该区形变时间序列和年均形变速率，并结合水头数据分析了不同区域含水系统对水头变化的不同响应.基于谐波函数分离了地表形变及水头变化的季节性变化，并用其估算了邯郸平原区空间差异变化的弹性骨架释水系数；基于顾及弱透水层延迟排水的一维水头变化-形变模型，反演了非弹性骨架释水系数和时间常数.结果表明：

(1)邯郸平原区以沉降为主，区内共有2个沉降速率大于20 mm·a-1的沉降区，其中北部沉降区地面沉降最为严重，年均沉降速率达14 cm·a-1，沉降区地面沉降无减缓趋势.

(2)研究区形变对水头变化的响应差异明显，地面沉降主要是由承压含水层水头下降及弱透水层的延迟排水引起.

(3)水头的季节性变化引起了明显的季节性形变，沉降区季节性形变幅度可达25 mm，峰值时间为1—3月.

(4)研究区弹性骨架释水系数介于1.51×10-3～4.05×10-3之间，与抽水试验结果较为相符；非弹性骨架释水系数为3.62×10-2～4.57×10-2，时间常数为0.47～0.77 a；当水头低于前期固结水头时，地表形变以非弹性形变为主，弹性形变占比仅为8%左右.

本文中6个承压地下水水头观测井均未位于沉降中心，因此其水头变化难以正确反映邯郸平原区深层地下水的整体变化趋势.利用少数承压地下水头观测，分析地面沉降对水头变化的响应规律，在此基础上，根据InSAR地表形变监测结果研究邯郸平原区(尤其是沉降区)深层地下水的整体变化趋势，一定程度上克服了水井观测较为稀疏，难以准确精细刻画区域地下水变化的不足.表明利用时序InSAR技术监测地下水变化，估算承压含水层参数具一定的优势及良好的应用潜力.此外，近年来时序InSAR技术已逐步用于承压地下水储量变化估算(Smith et al., 2017; Jiang et al., 2018; Bai et al., 2022)，随着新一代SAR数据愈加丰富，时序InSAR技术将成为地下水变化研究的重要手段.

致谢感谢欧洲空间局提供的Sentinel-1A数据，英国利兹大学Hooper博士提供的StaMPS软件.