利用時(shí)序InSAR技術(shù)反演邯鄲平原區(qū)地表形變與含水層參數(shù)

白林，李振洪，宋莎，劉東

1 長安大學(xué)地質(zhì)工程與測繪學(xué)院，西安 710054 2 西部礦產(chǎn)資源與地質(zhì)工程教育部重點(diǎn)實(shí)驗(yàn)室，西安 710054 3 長安大學(xué)地學(xué)與衛(wèi)星大數(shù)據(jù)研究中心，西安 710054 4 機(jī)械工業(yè)勘察設(shè)計(jì)研究院有限公司，西安 710043

0 引言

從承壓含水層抽取地下水時(shí)，會(huì)引起水頭下降，孔隙水壓力減小，導(dǎo)致有效應(yīng)力增大，含水系統(tǒng)骨架發(fā)生壓縮(即地面沉降)，造成儲(chǔ)水能力損失(Jiang et al., 2018).作為我國水資源最為短缺和環(huán)境問題最為突出的地區(qū)之一，華北平原地下水供水量達(dá)總供水量的近70%(邵景力等，2009).20世紀(jì)70年代以來，華北平原地下水經(jīng)歷了長期的超采過程，已形成二十多個(gè)大型復(fù)合水位降落漏斗(張兆吉等, 2009; 葛大慶等，2014)，導(dǎo)致了大面積的快速地面沉降，成為我國沉降面積最大、類型最復(fù)雜的地區(qū)之一.不均勻地面沉降不僅會(huì)對地下和地面建筑物、構(gòu)筑物產(chǎn)生破壞，威脅高速、高鐵等重大基礎(chǔ)設(shè)施安全運(yùn)營，而且會(huì)導(dǎo)致地面塌陷和地裂縫等地質(zhì)災(zāi)害，給人們的生產(chǎn)、生活帶來巨大的危害及損失(Bai et al., 2016).同時(shí)，地面沉降還會(huì)導(dǎo)致濱海平原地面高程資源損失，與全球海平面上升產(chǎn)生累加效應(yīng)，加劇了風(fēng)暴潮、海岸侵蝕、海水入侵等生態(tài)環(huán)境災(zāi)害的致災(zāi)程度(Xue et al., 2005).隨著華北平原地面沉降災(zāi)害的不斷發(fā)展，已極大地影響和制約著當(dāng)?shù)亟?jīng)濟(jì)建設(shè)可持續(xù)發(fā)展及社會(huì)安定.因此，為了預(yù)防和治理華北平原地面沉降地質(zhì)災(zāi)害，亟需查明地面沉降的時(shí)空演變特征，厘清承壓地下水變化與地面沉降的耦合關(guān)系.

GPS、精密水準(zhǔn)測量能夠獲取高精度的地面沉降監(jiān)測結(jié)果，但因其需要現(xiàn)場布設(shè)且監(jiān)測點(diǎn)稀疏，難以探測區(qū)域尺度地面沉降的時(shí)空分布特征.作為近年來發(fā)展迅速的空間大地測量新技術(shù)，時(shí)間序列InSAR(Multi Temporal Interferometric Synthetic Aperture Radar，MT-InSAR)技術(shù)具有全天時(shí)、全天候、精度高、覆蓋范圍大、時(shí)間/空間分辨率高等優(yōu)點(diǎn)(Hooper et al., 2012; Bai et al., 2016)，已廣泛應(yīng)用于華北平原地面沉降監(jiān)測(葛大慶, 2013; Yang et al., 2018).

目前已有相關(guān)學(xué)者成功將時(shí)序InSAR應(yīng)用于含水層參數(shù)反演(Hoffmann et al., 2001; Chaussard et al., 2014).Hoffmann等(2001)首次提出了基于InSAR觀測反演含水層參數(shù)的方法，估算了Las Vegas Valley地區(qū)測井附近的彈性骨架釋水系數(shù)，結(jié)果與實(shí)測數(shù)據(jù)較為一致.當(dāng)水頭低于前期固結(jié)水頭時(shí)，弱透水層往往存在延遲排水現(xiàn)象，若不考慮其影響，會(huì)導(dǎo)致非彈性骨架釋水系數(shù)低估.因此，Hoffmann等(2003)提出了顧及延遲排水的含水層參數(shù)反演算法，估算了Antelope Valley地區(qū)的時(shí)間常數(shù)及非彈性骨架釋水系數(shù).隨著時(shí)序InSAR技術(shù)的發(fā)展和SAR數(shù)據(jù)的積累，獲取某一地區(qū)形變時(shí)間序列成為可能，大大拓展了InSAR在地下水相關(guān)研究的應(yīng)用領(lǐng)域.許文斌等(2012)利用InSAR和地下水水位數(shù)據(jù)估算了洛杉磯地區(qū)的地表形變及含水層參數(shù).Chaussard等(2014, 2017)分別根據(jù)1992—2011年多源SAR數(shù)據(jù)和2011—2017年COSMO-SkyMed數(shù)據(jù)估算了Santa Clara Valley地區(qū)的彈性骨架釋水系數(shù)和彈性水儲(chǔ)量變化.Reeves等(2011, 2014)反演了San Luis Valley農(nóng)業(yè)區(qū)的彈性骨架釋水系數(shù)，并用其預(yù)測了水頭變化.Chen等(2016)使用ALOS數(shù)據(jù)估算了San Luis Valley地區(qū)的彈性骨架釋水系數(shù)和遲滯時(shí)間.近年來，國內(nèi)學(xué)者利用時(shí)序InSAR技術(shù)在華北平原開展了一些系列研究，但這些研究多集中于北京(Guo et al., 2020; Yu et al., 2020)、滄州(Jiang et al., 2018; Bai et al., 2022)等地區(qū).華北平原水文地質(zhì)條件復(fù)雜，地表形變對水頭變化的響應(yīng)差異顯著(何慶成等，2006)，因此開展更大范圍的地面沉降及含水層參數(shù)反演研究具有重要的意義.

本文將以華北平原承壓地下水超采、地面沉降災(zāi)害嚴(yán)重的邯鄲平原區(qū)為研究區(qū)，利用2015—2019年Sentinel-1A數(shù)據(jù)監(jiān)測地表形變的時(shí)空分布特征，并結(jié)合水頭數(shù)據(jù)分析含水系統(tǒng)對水頭變化的不同響應(yīng)；根據(jù)承壓地下水變化與地表形變的耦合關(guān)系，估算邯鄲平原區(qū)承壓含水系統(tǒng)的彈性/非彈性骨架釋水系數(shù)及時(shí)間常數(shù).

1 研究區(qū)與數(shù)據(jù)

1.1 研究區(qū)概況

邯鄲位于河北省南部，太行山脈東麓，地理坐標(biāo)為東經(jīng)113°28′—115°28′，北緯36°04′—37°01′，總面積約1.2萬km2(如圖1所示).邯鄲市地勢自西向東呈階梯狀下降，高差懸殊，地貌類型復(fù)雜多樣.以京廣鐵路為界，西部為中、低山丘陵地貌，東部為華北平原.邯鄲平原區(qū)約占全市面積的54%，按成因和地貌形態(tài)特征可分為山前沖積洪積平原和中東部沖積湖積平原(吳旭等，2012).平原區(qū)第四系含水層組厚度為350～500 m，可分為四個(gè)含水層組(第I～I(xiàn)V含水層組)，第Ⅰ、Ⅱ含水層組普遍分布有咸水、微咸水，淡水主要賦存于第II含水層組的下段及第III、IV含水層組，地下水主要開采層組為第II含水層組的下段及第III含水層組(石立明，2020).平原區(qū)地下水開采始于20世紀(jì)70年代末，隨著區(qū)域經(jīng)濟(jì)迅猛發(fā)展，開采量不斷擴(kuò)大，水位逐年迅速下降，至2010年中心水位下降總量超70 m，形成了面積達(dá)3063 km2的深層承壓地下水降落漏斗，并與衡水、滄州復(fù)合大型漏斗相連(付丹平，2013).承壓地下水長期持續(xù)超采，導(dǎo)致了嚴(yán)重的地面沉降，至2015年邯鄲平原區(qū)沉降中心累計(jì)沉降量已超1500 mm，沉降量大于500 mm的面積為3576 km2，占平原區(qū)總面積的47.5%(河北省國土資源廳, 2016).為了保護(hù)地下水資源，控制地面沉降災(zāi)害，邯鄲市自2014年被列為全省地下水超采綜合治理首批試點(diǎn)市以來，實(shí)施了一系列的地下水壓采措施，并增加了地表水供給，使得區(qū)內(nèi)地下水超采得到了有效遏制，截至2018年年底，累計(jì)壓減地下水超采量4.02億m3(石立明，2020).

圖1 邯鄲市地理位置白色方框表示水頭觀測井位置，藍(lán)色矩形表示所用Sentinel-1A數(shù)據(jù)覆蓋范圍，綠色實(shí)線表示含砂率，棕色實(shí)線表示斷裂.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 實(shí)驗(yàn)數(shù)據(jù)

本文所用SAR數(shù)據(jù)為2015年7月30日—2019年12月12日獲取的117景Sentinel-1A升軌數(shù)據(jù)，所處波段為C波段，極化方式為VV極化，方位向和距離向分辨率分別為20 m和5 m.原始SAR數(shù)據(jù)幅寬為250 km，選取了其中170 km×130 km的感興趣區(qū)，覆蓋范圍如圖1所示.此外，本文使用覆蓋實(shí)驗(yàn)區(qū)的30 m分辨率SRTM DEM數(shù)據(jù)去除干涉圖中地形相位并進(jìn)行地理編碼.

本文收集了邯鄲平原區(qū)6個(gè)承壓地下水水頭觀測井2010年至2017年的觀測數(shù)據(jù)(空間分布如圖1所示)，用于分析地面沉降與水頭變化的響應(yīng)關(guān)系和計(jì)算含水層參數(shù).本文未收集到井位的觀測層位及深度信息，但觀測井最初用于地下水開采，因此觀測井記錄的是第II含水層組的下段及第III含水層組的水頭變化.

2 數(shù)據(jù)處理方法

2.1 時(shí)序InSAR數(shù)據(jù)處理

邯鄲平原區(qū)以農(nóng)田為主，人工目標(biāo)點(diǎn)較為稀少，此外，華北平原地下水超采引起的地表形變往往伴隨著明顯的季節(jié)性波動(dòng)(葛大慶等，2014)，因此本文選擇適用于非城市區(qū)域，且具有反演地表周期性過程能力的StaMPS方法(Hooper et al., 2004, 2007a)進(jìn)行邯鄲地區(qū)地表形變監(jiān)測.本節(jié)將對該方法進(jìn)行簡要介紹.

首先，為了抑制由于時(shí)間/空間基線過長引起的失相干，設(shè)置時(shí)間基線閾值為36 d，空間基線閾值為150 m，并根據(jù)生成的干涉圖質(zhì)量及相干性進(jìn)行挑選，最終獲取340個(gè)干涉對用于時(shí)間序列分析，具體干涉對組合如圖2所示.為進(jìn)一步提高相干性，分別在方位向去除不重疊的多普勒頻率，在距離向進(jìn)行頻域?yàn)V波(Hooper, 2008).然后，利用軌道數(shù)據(jù)和SRTM DEM數(shù)據(jù)去除干涉圖中地形相位，獲得差分干涉圖，并對其進(jìn)行地理編碼.

圖2 Sentinel-1A干涉對時(shí)空基線圖Fig.2 Temporal/spatial baselines of Sentinel-1A interferometric pairs

StaMPS方法是基于濾波相位失相干緩慢(Slowly-Decorrelating Filtered Phase，SDFP)的目標(biāo)點(diǎn)進(jìn)行時(shí)間序列分析的(Hooper, 2008).為提高計(jì)算效率，首先對所有像元進(jìn)行幅度穩(wěn)定性分析，設(shè)置幅度差離散度閾值為0.6選取SDFP候選點(diǎn).然后，對獲取的SDFP候選點(diǎn)進(jìn)行相位穩(wěn)定性分析，剔除不滿足相位穩(wěn)定性標(biāo)準(zhǔn)的候選點(diǎn)，獲取SDFP點(diǎn)集(Hooper, 2008).利用三維解纏算法進(jìn)行相位解纏(Hooper et al., 2007b)，獲取SDFP點(diǎn)的解纏相位.為獲取形變相位，需從解纏相位中消除DEM誤差、軌道誤差、大氣延遲等干擾項(xiàng)的影響.本文根據(jù)DEM誤差與垂直基線的線性關(guān)系反演估算并去除其影響(Hooper et al., 2007a).因?yàn)楹惼皆貐^(qū)地勢平坦(高程范圍為20～60 m)，所以大氣延遲中垂直分層部分對形變監(jiān)測的影響較小(Liu et al., 2016)，通過估算并去除每個(gè)干涉圖的相位趨勢面以消除軌道誤差及大氣延遲長波段部分的影響(Hooper et al., 2007a).最后，根據(jù)形變相位與形變量的線性關(guān)系即可反演獲取每個(gè)SDFP點(diǎn)視線向形變量時(shí)間序列及平均形變速率.本文未對大氣延遲中湍流分量進(jìn)行改正，通過計(jì)算山前平原穩(wěn)定區(qū)形變時(shí)間序列的標(biāo)準(zhǔn)差來評估其對形變結(jié)果的影響，結(jié)果表明該區(qū)域形變時(shí)間序列標(biāo)準(zhǔn)差的平均值為1.5 cm，因湍流分量具有明顯的隨機(jī)性(楊成生等, 2011)，其對形變速率結(jié)果影響較小.邯鄲地區(qū)沉降主要由承壓地下水超采引起，與垂直形變相比，地下水超采引起的水平形變相對較小(Jiang et al., 2018; Yang et al., 2018)，因此本文假設(shè)水平形變可忽略，那么根據(jù)InSAR形變監(jiān)測結(jié)果及SAR影像入射角即可估算邯鄲地區(qū)垂直形變.

2.2 骨架釋水系數(shù)估算方法

承壓含水系統(tǒng)由含水層(砂、礫石)和弱透水層(黏土、淤泥)組成.骨架釋水系數(shù)是表征含水系統(tǒng)釋水能力的參數(shù)，指水頭下降(升高)一個(gè)單位時(shí)，單位面積含水系統(tǒng)由于骨架壓縮(膨脹)而釋出(儲(chǔ)存)的水的體積.承壓含水系統(tǒng)骨架釋水系數(shù)可根據(jù)含水系統(tǒng)厚度變化量與水頭變化的線性關(guān)系估算(Riley, 1969)：

Sk=Sskb0=Δb/Δh,

(1)

式中，Sk為承壓含水系統(tǒng)骨架釋水系數(shù)，Ssk為承壓含水系統(tǒng)骨架釋水率，b0為承壓含水系統(tǒng)原始厚度，Δb為地表形變量，假設(shè)其等于承壓含水系統(tǒng)厚度變化，Δh為承壓含水系統(tǒng)水頭變化.

骨架釋水系數(shù)在有效應(yīng)力超過歷史最大有效應(yīng)力前后差別很大，這個(gè)歷史最大有效應(yīng)力稱為前期固結(jié)應(yīng)力，對應(yīng)的水頭為前期固結(jié)水頭(Leake, 1990).有效應(yīng)力小于前期固結(jié)應(yīng)力時(shí)，若有效應(yīng)力增加(或水頭降低)，則導(dǎo)致含水層和弱透水層發(fā)生彈性壓縮，當(dāng)有效應(yīng)力恢復(fù)到初值時(shí)，壓縮可恢復(fù).如果有效應(yīng)力超過前期固結(jié)應(yīng)力，弱透水層將發(fā)生非彈性壓縮，這種非彈性壓縮是沉積物顆粒的移動(dòng)和重新排列導(dǎo)致的，很大程度上是永久性的.相同有效應(yīng)力變化條件下，弱透水層的非彈性壓縮遠(yuǎn)遠(yuǎn)大于彈性壓縮，前者可以比后者大一到兩個(gè)量級(jí)(Hoffmann et al., 2001).含水層中砂、礫石強(qiáng)度較高一般不發(fā)生非彈性壓縮(Smith et al., 2017).

為了刻畫有效應(yīng)力超過前期固結(jié)應(yīng)力前后骨架釋水系數(shù)的顯著變化，經(jīng)常用到以下兩個(gè)參數(shù)：

(2)

式中，Ske為彈性骨架釋水系數(shù)，Skv為非彈性骨架釋水系數(shù)，σe為有效應(yīng)力，σe(max)為前期固結(jié)應(yīng)力.非彈性骨架釋水系數(shù)通常比彈性骨架釋水系數(shù)大一到兩個(gè)量級(jí)(Hoffmann et al., 2003).

前期固結(jié)水頭難以直接測量，且會(huì)隨著地下水開采而發(fā)生變化，在實(shí)際工作中往往難以確定地表形變是否為彈性，這給彈性/非彈性骨架釋水系數(shù)計(jì)算帶來了困難.Miller和Shirzaei(2015)提出了基于季節(jié)性形變及水頭變化估算彈性骨架釋水系數(shù)的方法.這一方法的基本假設(shè)是季節(jié)性形變是彈性形變的一部分，是由水頭的季節(jié)性變化引起的.則根據(jù)形變與水頭變化的線性關(guān)系(式(1))，彈性骨架釋水系數(shù)可用(3)式計(jì)算：

Ske=Δbs/Δhs,

(3)

式中，Δbs為含水系統(tǒng)厚度季節(jié)性變化，即季節(jié)性形變，Δhs為水頭季節(jié)性變化.由于含水層不發(fā)生非彈性壓縮，因此無論水頭是否高于前期固結(jié)水頭，含水層骨架壓縮都是彈性的.當(dāng)水頭高于前期固結(jié)水頭時(shí)，含水層和弱透水層骨架壓縮都是彈性的，因此式(3)計(jì)算的是含水層和弱透水層彈性骨架釋水系數(shù)之和.當(dāng)水頭低于前期固結(jié)水頭時(shí)，只有含水層骨架壓縮為彈性，因此式(3)計(jì)算的僅為含水層彈性骨架釋水系數(shù).

由于弱透水層由細(xì)顆粒沉積物組成，垂向滲透系數(shù)較低，含水系統(tǒng)中弱透水層的水頭變化通常滯后于周圍含水層(Hoffmann et al., 2003).因此弱透水層壓縮也滯后于水頭變化，滯后的時(shí)間稱為時(shí)間常數(shù)，水頭下降引起的93%壓縮量是在該時(shí)間段內(nèi)發(fā)生的.非彈性骨架釋水系數(shù)及時(shí)間常數(shù)可用(4)式計(jì)算：

(4)

式中，τ0為時(shí)間常數(shù)，Δbinelastic為非彈性形變，可用(5)式計(jì)算：

Δbinelastic=Δb-SkeΔh.

(5)

2.3 延遲沉降分析

當(dāng)水頭停止下降時(shí)，由于弱透水層的滲透系數(shù)較低，弱透水層延遲排水導(dǎo)致的延遲沉降仍將繼續(xù)，可用(6)式表示：

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

(6)

式中，M為幅度；k為衰減系數(shù)，取值范圍為-1～0，表征含水系統(tǒng)對水頭變化的響應(yīng)速度.k越小，沉降速率減緩越快；k越大，沉降速率減緩越慢；如果k趨近于0則表明沉降無減緩趨勢，地面沉降不是由弱透水層延遲排水引起的.

2.4 季節(jié)性信號(hào)分離

為估算彈性骨架釋水系數(shù)，本文使用諧波函數(shù)提取季節(jié)性形變及季節(jié)性水頭變化(Neely et al., 2021)：

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

(7)

式中，Y(t)為形變/水頭變化時(shí)間序列，v為形變/水頭線性變化速率，A為季節(jié)性變化幅度，t為時(shí)間，T為季節(jié)性信號(hào)達(dá)到峰值的時(shí)刻，B為常數(shù).

2.5 水頭變化與形變的相關(guān)性分析

水頭的季節(jié)性變化會(huì)引起明顯的季節(jié)性形變，由于含水系統(tǒng)存在低滲透介質(zhì)，形變往往滯后于水頭變化.利用(8)式對時(shí)間序列進(jìn)行相關(guān)性分析，可估算形變滯后水頭變化的時(shí)間：

(8)

式中，τlag為形變滯后水頭變化的時(shí)間，corr表示相關(guān)系數(shù)，Δhdetrend為去除線性趨勢后的殘余水頭時(shí)間序列，Δbdetrend為去除線性趨勢后的殘余形變時(shí)間序列.在進(jìn)行相關(guān)性分析時(shí)，需將殘余形變時(shí)間序列線性插值從而獲取與水頭觀測時(shí)間對應(yīng)的形變結(jié)果.

3 結(jié)果與討論

3.1 地表形變監(jiān)測結(jié)果與解譯

2015—2019年邯鄲平原區(qū)地表形變監(jiān)測結(jié)果如圖3所示，地面沉降分布廣泛，西部受承壓含水層分布控制，東部受館陶西斷裂及元村斷裂控制，北部與邢臺(tái)沉降區(qū)相連.邯鄲平原區(qū)可分為兩個(gè)沉降區(qū)：北部沉降區(qū)(圖3中區(qū)域S1，以曲周縣、肥鄉(xiāng)縣、廣平縣等為中心的有咸水區(qū))和南部沉降區(qū)(圖3中區(qū)域S2，以臨漳縣為中心)，兩個(gè)沉降區(qū)已連為一片，最大沉降速率可達(dá)14 cm·a-1.地面沉降主要是由承壓地下水超采導(dǎo)致的含水系統(tǒng)骨架壓縮引起的，為進(jìn)一步分析含水系統(tǒng)對水頭變化的不同響應(yīng)，本文比較了2個(gè)典型觀測井水頭變化與其附近SDFP點(diǎn)的形變時(shí)間序列，結(jié)果如圖4所示.D3觀測井2015—2017年水頭呈上升趨勢，但地面沉降仍在繼續(xù)(圖4a)，表明該區(qū)域地面沉降主要是由于弱透水層延遲排水引起的.受逐年升高的水頭影響，沉降速率快速減緩，衰減系數(shù)達(dá)-0.74，到2017年末沉降速率已趨近于0.D4觀測井2015—2017年水頭變化趨勢和沉降趨勢大體上是一致的，地面沉降速率無明顯的減緩趨勢(圖4b)，表明該沉降主要是由2015—2017年水頭下降引起.相較于2015—2017年，2個(gè)典型觀測井附近區(qū)域2018—2019年沉降速率明顯加快，可能的原因是該區(qū)域承壓地下水發(fā)生了快速下降.兩個(gè)觀測井含水系統(tǒng)對水頭變化的不同響應(yīng)表明邯鄲平原區(qū)含水層和弱透水層的空間分布及厚度差異顯著，這一現(xiàn)象和滄州地區(qū)較為相似(Jiang et al., 2018).由于部分地區(qū)地面沉降發(fā)展趨勢在2018年前后發(fā)生了明顯變化，本文分別計(jì)算了2015—2017年和2018—2019年的邯鄲平原區(qū)地面沉降的衰減系數(shù).由圖5a可知，2015—2017年西部山前平原、魏縣東部及館陶縣東部衰減系數(shù)在-1～-0.4之間，尤其是邯鄲市區(qū)及磁縣、永年縣、魏縣、館陶縣城區(qū)衰減系數(shù)可達(dá)-1，地面沉降存在明顯的減緩趨勢；但在中部沉降區(qū)衰減系數(shù)均大于-0.1，地面沉降無明顯的減緩趨勢.2018—2019年邯鄲平原區(qū)沉降衰減系數(shù)明顯增大，除西部山前平原部分地區(qū)、魏縣東部及館陶縣城區(qū)外，衰減系數(shù)均大于-0.1(圖5b)，地面沉降無明顯的減緩趨勢，其空間分布在東部與沉降分布大致一致.造成這一現(xiàn)象可能的原因是2019年降水量僅為多年平均量的87%，導(dǎo)致承壓地下水開采量增加，地面沉降加劇.以上結(jié)果表明地下水壓采措施在部分城市地區(qū)取得了較好的效果，但在地面沉降嚴(yán)重的中部地區(qū)效果不顯著，主要是因?yàn)樵搮^(qū)域地表水資源匱乏，難以滿足該地區(qū)的用水需求.

圖3 邯鄲平原區(qū)2015—2019年平均形變速率圖白色十字表示參考點(diǎn)位置，白色實(shí)線表示斷裂，白色方框表示2個(gè)典型觀測井的位置，灰色實(shí)線表示承壓含水層邊界.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 季節(jié)性形變結(jié)果

承壓水水頭的季節(jié)性變化導(dǎo)致形變存在明顯的季節(jié)性波動(dòng)(圖4c和4d)，且形變滯后于水頭變化約45 d(圖6)，這主要是由于該含水系統(tǒng)垂直滲透系數(shù)較低，被開采含水層水頭需較長時(shí)間才能與相鄰含水層或弱透水層達(dá)到平衡.基于諧波函數(shù)提取了典型井位的季節(jié)性形變(2015—2016年D3觀測井無明顯季節(jié)性水頭變化及季節(jié)性形變，因此僅提取了2017年季節(jié)性形變)，其幅度為7～15 mm，形變峰值對應(yīng)時(shí)間為每年的2—4月.從圖4中可看出，2015—2017年及2018—2019年兩個(gè)時(shí)段的季節(jié)性形變幅度有明顯變化，可能的原因是降水變化導(dǎo)致的地下水開采強(qiáng)度變化引起的.為分析邯鄲平原區(qū)季節(jié)性形變的時(shí)空分布特征，分別計(jì)算了2015—2017年及2018—2019年的季節(jié)性形變幅度及峰值時(shí)間(圖7).結(jié)果表明，2015—2017年邯鄲平原區(qū)北部沉降區(qū)及東南部區(qū)域存在明顯的季節(jié)性形變，北部沉降區(qū)幅度為5～20 mm，東南部區(qū)域幅度為10～30 mm(圖7a).受降雨減少的影響，2018—2019年沉降區(qū)域地下水超采加劇，引起的水頭季節(jié)性變化幅度增加導(dǎo)致了沉降區(qū)季節(jié)性形變幅度增大，最大幅度可達(dá)25mm(圖7b).東南部區(qū)域地面沉降呈現(xiàn)減緩趨勢(圖5b)，表明該區(qū)域深層承壓地下水開采強(qiáng)度有所降低，導(dǎo)致季節(jié)性形變幅度減小至5～15 mm(圖7b).自然狀態(tài)下，地下水受降水影響，水頭在每年的8—9月份達(dá)到最高值.但邯鄲平原區(qū)地下水變化主要受農(nóng)業(yè)灌溉控制，水頭在每年的1—3月份達(dá)到最高值，6—7月份達(dá)到最低值(付丹平，2013).北部沉降區(qū)淺層地下水(潛水)為咸水，農(nóng)業(yè)灌溉主要依賴承壓地下水，因此2015—2017年該區(qū)域季節(jié)性形變峰值時(shí)間為1—3月(圖7c).其他區(qū)域峰值時(shí)間主要為8—9月，與自然狀態(tài)下承壓地下水最高水頭時(shí)間相符，表明其灌溉主要依賴淺層地下水，深層承壓地下水受淺層地下水越流補(bǔ)給及側(cè)向徑流補(bǔ)給影響，在8—9月達(dá)到最高水頭.受降雨減少的影響，2018—2019年承壓地下水超采范圍增大，季節(jié)性形變峰值時(shí)間為1—3月的范圍擴(kuò)展到整個(gè)沉降區(qū)域及山前平原北部(圖7d)，表明這些區(qū)域的深層地下水開采量增加，引起了明顯的季節(jié)性形變.

圖4 (a)觀測井D3，(b)觀測井D4水頭時(shí)間序列與其附近SDFP點(diǎn)形變時(shí)間序列對比圖.(c)觀測井D3，(d)觀測井D4殘余水頭時(shí)間序列與殘余形變時(shí)間序列對比結(jié)果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年邯鄲平原區(qū)地面沉降衰減系數(shù)白色實(shí)線表示斷層，紫色實(shí)線表示沉降速率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 去除線性趨勢后的殘余形變與殘余水頭時(shí)間序列相關(guān)性分析Fig.6 Correlation analysis between the detrend deformation and the detrend head time series

圖7 (a)2015—2017年和(b)2018—2019年季節(jié)性形變幅度.(c)2015—2017年和(d)2018—2019年季節(jié)性形變峰值時(shí)間紫色實(shí)線表示沉降速率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 骨架釋水系數(shù)估算結(jié)果

利用季節(jié)性形變及季節(jié)性水頭變化估算了邯鄲平原區(qū)彈性骨架釋水系數(shù).由于觀測井D2的水頭沒有明顯的季節(jié)性變化，未用于彈性骨架釋水系數(shù)估算，因此本文最終估算了5個(gè)觀測井位的彈性骨架釋水系數(shù)(表1).邯鄲平原區(qū)彈性骨架釋水系數(shù)估算結(jié)果介于1.51×10-3～4.05×10-3之間，與抽水試驗(yàn)結(jié)果1.0×10-3～8.0×10-3(張兆吉等，2009)較為相符.彈性骨架釋水系數(shù)由北向南大體上呈增大的趨勢(圖8)，可能的原因是北部地區(qū)淺層地下水為咸水，為滿足農(nóng)業(yè)用水需求大量開采深層承壓地下水，導(dǎo)致該區(qū)域開采厚度較大(式(1)中較大的b0).

當(dāng)有效應(yīng)力超過前期固結(jié)應(yīng)力(水頭低于前期固結(jié)水頭)時(shí)，表征單位面積含水系統(tǒng)釋水能力的參數(shù)稱為非彈性骨架釋水系數(shù).非彈性骨架釋水系數(shù)通常比彈性骨架釋水系數(shù)大一到兩個(gè)量級(jí)(Riley, 1969).本節(jié)采用2.2節(jié)所述方法估算非彈性骨架釋水系數(shù)和時(shí)間常數(shù)，需挑選水頭低于前期固結(jié)水頭且仍在下降，并引起明顯地面沉降的觀測井位.結(jié)合2000—2017年水頭觀測資料及InSAR地面沉降監(jiān)測結(jié)果，本文選取了2個(gè)滿足該條件的觀測井.非彈性骨架釋水系數(shù)估算結(jié)果為3.62×10-2～4.57×10-2(表1).在觀測井D4，非彈性骨架釋水系數(shù)約為彈性骨架釋水系數(shù)的11.75倍，表明當(dāng)水頭低于前期固結(jié)水頭時(shí)，地表形變以非彈性形變(永久性沉降)為主，可恢復(fù)的彈性形變占比僅為8%左右.觀測井D2和D4時(shí)間常數(shù)分別為0.47 a和0.77 a，而由圖4a可知，觀測井D3的時(shí)間常數(shù)應(yīng)大于2 a(水頭上升期：2016年1月—2017年12月)，但該觀測井距離D2和D4觀測井較遠(yuǎn)(圖8)，水文地質(zhì)條件可能存在較大差異，導(dǎo)致時(shí)間常數(shù)相差較大.此外，邯鄲地區(qū)時(shí)間常數(shù)明顯小于華北平原東部(6 a)(Bai et al., 2022)，表明該區(qū)域弱透水層滲透性較好，對水頭變化響應(yīng)速度較快.

表1 邯鄲平原區(qū)骨架釋水系數(shù)估算結(jié)果Table 1 The results of the estimated skeletal storativity across the Handan plain

圖8 邯鄲平原區(qū)彈性骨架釋水系數(shù)插值結(jié)果觀測井位置如圖中白色方框.Fig.8 Interpolation map of elastic skeletal storativity across the Handan plainThe white squares represent the well locations.

4 結(jié)論

本文利用2015—2019年117景Sentinel-1A數(shù)據(jù)對邯鄲平原區(qū)地表形變進(jìn)行監(jiān)測，獲取了該區(qū)形變時(shí)間序列和年均形變速率，并結(jié)合水頭數(shù)據(jù)分析了不同區(qū)域含水系統(tǒng)對水頭變化的不同響應(yīng).基于諧波函數(shù)分離了地表形變及水頭變化的季節(jié)性變化，并用其估算了邯鄲平原區(qū)空間差異變化的彈性骨架釋水系數(shù)；基于顧及弱透水層延遲排水的一維水頭變化-形變模型，反演了非彈性骨架釋水系數(shù)和時(shí)間常數(shù).結(jié)果表明：

(1)邯鄲平原區(qū)以沉降為主，區(qū)內(nèi)共有2個(gè)沉降速率大于20 mm·a-1的沉降區(qū)，其中北部沉降區(qū)地面沉降最為嚴(yán)重，年均沉降速率達(dá)14 cm·a-1，沉降區(qū)地面沉降無減緩趨勢.

(2)研究區(qū)形變對水頭變化的響應(yīng)差異明顯，地面沉降主要是由承壓含水層水頭下降及弱透水層的延遲排水引起.

(3)水頭的季節(jié)性變化引起了明顯的季節(jié)性形變，沉降區(qū)季節(jié)性形變幅度可達(dá)25 mm，峰值時(shí)間為1—3月.

(4)研究區(qū)彈性骨架釋水系數(shù)介于1.51×10-3～4.05×10-3之間，與抽水試驗(yàn)結(jié)果較為相符；非彈性骨架釋水系數(shù)為3.62×10-2～4.57×10-2，時(shí)間常數(shù)為0.47～0.77 a；當(dāng)水頭低于前期固結(jié)水頭時(shí)，地表形變以非彈性形變?yōu)橹鳎瑥椥孕巫冋急葍H為8%左右.

本文中6個(gè)承壓地下水水頭觀測井均未位于沉降中心，因此其水頭變化難以正確反映邯鄲平原區(qū)深層地下水的整體變化趨勢.利用少數(shù)承壓地下水頭觀測，分析地面沉降對水頭變化的響應(yīng)規(guī)律，在此基礎(chǔ)上，根據(jù)InSAR地表形變監(jiān)測結(jié)果研究邯鄲平原區(qū)(尤其是沉降區(qū))深層地下水的整體變化趨勢，一定程度上克服了水井觀測較為稀疏，難以準(zhǔn)確精細(xì)刻畫區(qū)域地下水變化的不足.表明利用時(shí)序InSAR技術(shù)監(jiān)測地下水變化，估算承壓含水層參數(shù)具一定的優(yōu)勢及良好的應(yīng)用潛力.此外，近年來時(shí)序InSAR技術(shù)已逐步用于承壓地下水儲(chǔ)量變化估算(Smith et al., 2017; Jiang et al., 2018; Bai et al., 2022)，隨著新一代SAR數(shù)據(jù)愈加豐富，時(shí)序InSAR技術(shù)將成為地下水變化研究的重要手段.

致謝感謝歐洲空間局提供的Sentinel-1A數(shù)據(jù)，英國利茲大學(xué)Hooper博士提供的StaMPS軟件.