摘要

河西地区地下水利用比重不断上升导致地下水位显著下降,引起了局部地区地面沉降。研究河西地区地下水变化与地面沉降滞后性对当地水资源管理、土地利用规划和农业发展具有重要意义。利用GRACE与GLDAS数据得到研究区2010—2017年地下水变化速率,结合监测井实测数据验证了反演地下水变化数据的可靠性,利用小基线集合成孔径雷达干涉测量(small baseline subset interferometry synthetic aperture Radar,SBAS-InSAR)技术得到局部沉降区2014年10月—2017年6月的地表形变速率,并用永久散射体合成孔径雷达干涉测量(persistent scatters interferometry synthetic aperture Radar,PS-InSAR)技术对结果进行对比验证,运用快速傅里叶变换和时滞相关性分析对地下水变化与地表沉降数据解算分析。结果表明,临泽、甘州、凉州、金川沉降区地面沉降较地下水变化滞后时间分别为74~86 d,61~80 d,80~99 d,74~99 d; 相关系数分别在0.541~0.593,0.589~0.689,0.600~0.750,0.543~0.630之间。研究结果可为河西地区水资源管理、土地利用规划和农业发展提供科学依据。

引用

魏小强, 杨国林, 刘涛, 邵明, 马志刚. 基于GRACE与InSAR数据地下水变化与地面沉降滞后性研究[J]. 自然资源遥感, 2025, 37(1): 122-130.

0 引言

河西走廊位于内陆干旱地区,水资源短缺,地下水作为地表水的补充,是河西地区农业灌溉、工矿和城市的重要水源,近年来,地下水利用比重不断上升导致地下水位显著下降,引起了局部地区地面沉降,直接影响了当地水资源管理、农业生产、生态保护和城市建设的可持续发展。因此,揭示地下水变化与地面沉降滞后性对该地区发展具有广泛而深远的影响。

当前监测地下水变化主要是利用地下水井观测,该方法观测范围小,难以满足大尺度地下水变化,而由美国航天局与德国航天局联合研制的GRACE卫星因其不受地形、气候等因素影响为反演地下水动态变化提供了新的途径。早期用于地表形变监测方法包括水准测量、全球定位系统等,然而这些方法获取数据所需成本高、覆盖范围小,耗时费力。随着合成孔径雷达干涉测量技术(interferometry synthetic aperture Radar,InSAR)发展,差分合成孔径雷达干涉测量(differential interferometric synthetic aperture Radar,D-InSAR)技术得到应用,与传统测量相比具有精度高、效率快、覆盖范围广等优点,已广泛运用于形变监测、滑坡灾害预测及地震形变,但D-InSAR技术易受时空失相干与大气延迟等因素制约。为克服这些因素影响,研究人员提出了永久散射体InSAR技术(persistent scatters -InSAR,PS-InSAR)和小基线集InSAR技术(small baseline subset -InSAR,SBAS-InSAR)。由此,Guo等利用GRACE数据和SBAS-InSAR监测技术得到北京地区地下水变化和地表沉降数据以及两者之间的相关性; Vasco等利用GRACE和合成孔径雷达(synthetic aperture Radar,SAR)数据得到图莱里盆地地表沉降与地下水储量变化的关系; Massoud等通过GRACE数据反演得到黎巴嫩平原地下水变化趋势并结合SAR数据和监测井数据证实了研究结果。基于上述研究可知,将GRACE与InSAR数据结合能够更好地揭示地下水变化对地表沉降的影响,但针对河西地区相关研究较少,未能有效、明确地阐述当地地面沉降对地表水变化的滞后响应。

本文利用GRACE和GLDAS数据解算得到研究区地下水年内、年际变化数据(缺失月份数据采用奇异谱插值法将其补齐),并用监测井实测数据验证反演地下水变化数据的可靠性。采用SBAS-InSAR技术得到地表形变数据,用PS-InSAR技术对形变结果进行对比验证,在此基础上运用快速傅里叶变换(fast Fourier transform,FFT)将地下水变化和地表沉降数据进行解算得到地面沉降较地下水变化滞后天数,通过时滞互相关分析两者的相关性,为今后河西地区在水资源管理、土地利用规划和农业发展方面提供了科学依据。

1、研究区概况及数据源1.1 研究区概况

河西走廊中东部是甘肃省主要的农业地区,包括张掖、金昌、武威3市,地处我国西北内陆干旱地区,区域内属于典型的山地、戈壁、绿洲复合生态系统,自然条件严酷,生态条件脆弱。开发和利用地下水资源在其经济社会发展中占有十分重要的地位,然而由于过度开发地下水资源,含水层孔隙水位下降,引起区域性地表沉降。由甘肃省公布的地下水超开采区域中,浅层地下水和深层承压水共划分出47处,本文选取河西走廊的张掖市临泽县、甘州区,武威市凉州区,金昌市金川区这4处超采区来探究地表沉降对地下水变化滞后的响应。

1.2 数据源

SAR数据采用Sentinel-1A单式复数影像及同期卫星精密轨道数据; 数字高程模型(digital elevation model,DEM)为GDEMV2 30m空间分辨率数据; GRACE数据采用德克萨斯大学空间研究中心(Center for Space Research at the University of Teaxs,CSR)提供的Level-2 RL06版本的GSM数据模型,时间范围为2003年3月—2017年6月,阶次为96; GLDAS数据为GLDAS_NOAH025_M_2.1版本的GLDAS模型数据[17],时间尺度为1个月,空间分辨率为0.25°××0.25°,时间范围为2003年3月—2017年6月; 实测地下水数据来源于《中国地质环境监测地下水位年鉴》,时间范围为2003—2016年。

2 研究方法2.1 PS-InSAR技术原理

PS-InSAR技术主要通过SAR数据定标和配准、差分干涉SAR数据处理、永久性散射体点(persistent scatters,PS)选择、多影像稀疏格网相位解缠、大气相位屏的估计和去除、像素点的时序分析和PS点形变值估计等步骤处理得到地面形变时间序列,处理流程如图1所示。


2.2 SBAS-InSAR技术原理

SBAS-InSAR技术是利用多幅主影像与基于短时空基线准则的InSAR时间序列分析方法,该方法首先将若干基线进行组合,产生若干对高相干的干涉图,然后基于最小二乘原理,通过相位解缠、空间滤波等方法对差分干涉图进行地表形变估计运算,再利用矩阵奇异值分解(singular value decomposition,SVD)法解算小基线集合获得地面形变时间序列,处理流程如图2所示。


2.3 GRACE时变重力场反演地下水储量

本文将GRACE数据中的低阶项C20用卫星激光测距所测精度更高的C20项替换,加入GRACE数据计算的地心改正项修正一阶系数

2.4 地面沉降与地下水变化滞后性

地下水变化与地面沉降是相关过程,通过FFT处理地下水变化与地面沉降数据,可得到二者在不同频率下的幅度和相位信息,若2个信号存在滞后关系,则它们在某些频率下的相位信息会发生变化,通过分析相位变化,可得到二者滞后时间。本文运用FFT将二者数据的时域信号转化为频域信号,然后共轭相乘,再进行FFT逆变换得到交叉相关性延迟数组,找到数组中最大值所对应的延迟位置,即为滞后时间。

3 结果与分析3.1 地面沉降3.1.1 形变监测结果

利用SBAS-InSAR技术对研究区Sentinel-1A数据进行地表形变特征信息提取,获得沉降区地表形变速率,如图3所示。由图可知,临泽沉降区地表形变速率在-45~40 mm/a之间,在平川镇、新华镇沉降达到了-36~-20 mm/a; 甘州沉降区地表形变速率在-30~35 mm/a之间,在上秦镇、长安镇地表形变速率达到了-25~-10 mm/a; 凉州沉降区地表形变速率在-45~30 mm/a之间, 在吴家井镇沉降达到了-30~-15 mm/a; 金川沉降区地表形变速率在-65~45 mm/a之间,在康盛村沉降达到了-45~-15 mm/a。实验结果与已有研究成果一致。从图3可知,研究区局部区域有抬升现象,依据现有的文献资料,原因有: ①随着时间推移,局部地区地下水位会通过灌溉水、降雨等方式得到恢复,导致局部区域出现隆起现象; ②人为土地利用也会使地表产生不同程度的抬升与沉降; ③季节性或年际地下水位波动可能会导致局部区域抬升。


3.1.2 PS-InSAR结果差异性对比

由于缺少同期地面观测数据,本文对各个沉降区利用PS-InSAR技术得到的PS点和SBAS-InSAR技术得到的散射体位移场点(scatterer displacement field point, SDFP)进行统计分析,结果如图4所示。



图4 年平均沉降速率分布

由图4可知,沉降区由PS-InSAR与SBAS-InSAR所获得的地表沉降点值在各个数值上分布趋势较为相似,说明2种方法得到的年平均沉降速率基本一致,但在选择备选点时由于2种方法局限性等因素,得到的结果仍有细微差别。

为进一步探究2种监测结果之间的差异,利用ArcGIS软件在沉降区域选取同名点(具有相同地理坐标)建立回归方程,分析2种结果之间的相关性,如图5所示。由图可知,PS点和SDFP点形变速率具有较高相关性,4个沉降区R2分别为0.837 23,0.781 44,0.813 82,0.796 93,说明2种方法得到的年平均沉降速率具有一致性,也证明了数据的可靠性。


3.2 地下水储量变化3.2.1 地下水储量变化速率空间分布

用反演的地下水变化数据可得到2010年1月—2017年6月年际地下水变化速率,如图6所示,临泽县、甘州区地下水变化速率为0~4 mm/a,水储量呈现缓慢增长趋势;金川区、凉州区地下水变化速率为-10~0 mm/a,水储量变化较临泽县、甘州区亏损严重。


3.2.2 地下水变化验证

为验证反演得到的地下水变化数据,本文采用研究区监测井实测数据与其进行对比,由于与GRACE结果相当的等效水位只能通过实测地下水与给水度之积得到,因此只对二者趋势进行定性分析,如图7所示。


由图7(a)可知,反演地下水与实测地下水在整体上波动较为一致,但反演地下水波动趋势比实测地下水更大。为了更加直观地表达二者趋势变化,将研究时间段多年数据取月平均处理,如图7(b)所示,在1—2月,二者均表现为下降趋势,3—5月缓慢上升,6—8月份水量最为充盈,反演地下水在7月份水量达到最高,较实测地下水提前一个月,从9月份开始,二者均呈下降趋势,对两者进行年内相关性分析,其相关系数为0.867,为极强相关,说明二者具有较高的季节性相关,也验证了反演地下水数据的可靠性。

3.3 地面沉降与地下水变化滞后性分析3.3.1 滞后时间计算

采用FFT将不同沉降区地下水变化与地面沉降数据进行解算,结果如表1所示。

表1地面沉降与地下水变化滞后时间


由表1可知,甘州沉降区沉降速率小于临泽、凉州、金川沉降区,通过FFT计算地表沉降较地下水变化滞后时间,临泽、甘州沉降区滞后时间少于凉州、金川沉降区。由图6地下水变化速率可知,临泽、甘州地区水储量较金川、凉州地区丰富,故地表水补充到地下含水层所需时间较短,与本文计算得出的滞后时间差具有一致性。

3.3.2 时滞互相关分析

依据式(3)时滞互相关将不同沉降区地下水变化与地面沉降数据进行解算,结果如表2所示。由表2可知,地面沉降与地下水变化之间的相关系数越大,地下水变化对地面沉降影响越明显,所占比例越高。

表2地面沉降与地下水变化相关系数


3.3.3 地面沉降与地下水变化关系

通过对每个沉降区若干点进行多项式拟合,得到各个沉降区地下水变化与地面沉降时间变化曲线,如图8所示。


由图8可知,地面沉降与地下水变化趋势较为一致且两者之间表现出一定的滞后性,说明当地下水位上升或下降时,地表沉降速率会相应的减慢或加快,但由于滞后性,二者之间有一定的时间差,与表1—2得出的数据具有一致性。

引起地表沉降的因素很多,且不同地区所属地理环境与地质土层结构不一样,为了进一步分析滞后时间与相关系数在空间上的表现,采用反距离平均插值法得到滞后时间及相关系数的空间变化,如图9所示。


由表1和图9可知,临泽、甘州地区地面沉降较地下水变化的滞后时间为61~86 d,凉州、金川地区为74~99 d; 临泽、甘州地区地下水变化与地面沉降相关系数也小于凉州、金川地区,临泽、甘州地区二者相关系数在0.541~0.689之间,凉州、金川地区相关系数在0.543~0.750之间。由此可得,地面沉降与地下水变化时滞相关性变化空间分布规律为: 由河西走廊东部向中部先逐渐减弱,再增大; 滞后时间也符合此规律,二者具有一致性。

4 结论

利用SBAS-InSAR技术和GRACE时变重力场反演数据得到了研究区地下水变化与地面沉降数据,并采用PS-InSAR技术和实测地下水数据验证了实验数据的可靠性。运用FFT和时滞互相关分析解算了地面沉降较地下水变化滞后时间及相关系数,主要结论如下:

1)2010年1月—2017年6月,凉州、金川地区地下水变化速率在-10~0 mm/a,水储量呈现亏损状态; 临泽、甘州地区地下水变化速率在0~4 mm/a。通过与实测地下水数据进行年内相关性分析(R为0.867),为极强相关,验证了地下水数据的可靠性。

2)利用SBAS-InSAR技术,得到了沉降区地表形变数据,采用PS-InSAR技术对形变数据进行对比验证,得出2种方法得到的年平均沉降速率具有一致性,证明了数据的可靠性。

3)采用FFT计算地面沉降较地下水变化滞后时间,得到临泽、甘州、凉州、金川沉降区地面沉降较地下水变化滞后时间分别为74~86 d,61~80 d,80~99 d,74~99 d,通过时滞互相关计算,得到相关系数区间分别为0.541~0.593,0.589~0.689,0.600~0.750,0.543~0.630。

4)通过对地面沉降较地下水变化滞后时间及二者相关系数空间分析可知,当地下水变化时,地表会发生相应的变化,二者相关性与滞后时间变化规律具有一致性。

本文通过FFT和时滞互相关分析得出了河西地区地下水变化滞后时间及相关系数,为河西地区水资源管理、土地利用规划和农业发展等提供科学依据,但本研究只是单一分析了地下水变化与地面沉降滞后性,对研究区部分区域抬升现象和不同沉降等级下二者滞后性问题触及较少,下一步需结合水准、水位和降雨等数据进一步分析滞后性规律和抬升区域驱动机制。

(原文有删减)

来源:测绘学术资讯

ad1 webp
ad2 webp
ad1 webp
ad2 webp