基于GRACE 衛(wèi)星的黃淮海地區(qū)地下水儲量變化與歸因分析

王 坤，李玉龍，李琛曦，賴燕玲，楊傳國,3

(1.河海大學 水文水資源學院, 江蘇 南京 210098； 2.河北省地質(zhì)資源環(huán)境監(jiān)測與保護重點實驗室, 河北 石家莊050021； 3.河海大學 水災(zāi)害防御全國重點實驗室, 江蘇 南京 210098)

地下水是水資源的重要組成部分，由于具有分布廣、水質(zhì)好、可持續(xù)開發(fā)利用等許多優(yōu)點，已成為最重要的和最佳的供水選擇[1]，特別是在水資源相對匱乏的北方地區(qū)。黃淮海地區(qū)是中國重要的糧食生產(chǎn)基地[2]，是人均水資源量較少的地區(qū)之一，也是人類活動對水資源影響最嚴重的區(qū)域[3]。大規(guī)模開發(fā)地下水資源導致地下水位下降嚴重，產(chǎn)生了水質(zhì)惡化、地面沉降、海水入侵等一系列水安全問題[4]，已成為制約地區(qū)發(fā)展的主要因素。全面監(jiān)測地下水并掌握其變化規(guī)律是當務(wù)之急。傳統(tǒng)人工監(jiān)測地下水的方法耗費人力物力，且由于監(jiān)測站點布設(shè)的困難，難以實現(xiàn)地下水的大面積整體監(jiān)測。隨著重力衛(wèi)星探測技術(shù)的發(fā)展，通過測量地球重力場的變化反演全球陸地水儲量變化，為實現(xiàn)地下水的區(qū)域整體監(jiān)測提供了新途徑[5]。國內(nèi)外諸多學者利用GRACE （Gravity Recovery and Climate ）重力衛(wèi)星在不同區(qū)域開展了地下水研究，證實了在大尺度范圍內(nèi)反演地下水資源量變化的可靠性[6]。Famiglietti 等[7]和Scanlon 等[8]利用GRACE 數(shù)據(jù)反演加利福尼亞州中央山谷地下水儲量變化，與研究區(qū)實測地下水位驗證結(jié)果大體一致。國內(nèi)學者驗證了GRACE 與GLDAS （Global Land Data Assimilation System）數(shù)據(jù)反演黃河流域[9]、海河流域[10]、三江源[11]等流域地下水儲量變化具有較好適用性。

黃淮海地區(qū)地下水儲量變化引起越來越多的研究者關(guān)注，對其影響的主導因素開展了相關(guān)研究。Guo 等[12]研究發(fā)現(xiàn)人類活動在海河流域地下水儲量變化中占主導作用，農(nóng)業(yè)灌溉是地下水儲量消耗的主要因素，而降水和蒸發(fā)量相對穩(wěn)定；林雪敏等[13]研究表明農(nóng)業(yè)用水量對京津冀地區(qū)地下水儲量變化的影響程度明顯大于降水量與蒸散量；秦歡歡等[14]利用模型設(shè)定不同情景模擬華北平原水資源利用狀況，結(jié)果表明農(nóng)業(yè)節(jié)水措施和南水北調(diào)工程可以有效恢復地下水儲量。南水北調(diào)通水后黃淮海地區(qū)供水結(jié)構(gòu)發(fā)生了改變，相關(guān)研究較少，亟需采用多源數(shù)據(jù)深入分析區(qū)域地下水儲量變化的特征與成因。因此，本文利用最新的GRACE 數(shù)據(jù)和GLDAS 陸面同化數(shù)據(jù)，定量計算地下水儲量變化，評估黃淮海地區(qū)地下水儲量變化的時空特征，結(jié)合降水、蒸散發(fā)、地下水供水、跨流域調(diào)水等自然和人為因素，探討影響地下水儲量變化的原因，從而為科學評估與合理利用地下水提供依據(jù)。

1 研究區(qū)概況

本文選擇黃淮海地區(qū)，包括海河流域、淮河流域及黃河流域花園口以下地區(qū)為研究區(qū)域（圖1）。地理位置位于北緯31°～43°，東經(jīng)112°～122°，北起燕山南麓，南抵桐柏山、大別山，西起太行山，東臨渤海、黃海，總面積為67.8 萬km2。黃淮海地區(qū)屬溫帶季風氣候，季節(jié)變化明顯，冬季干燥寒冷，夏季高溫多雨[15]；年平均氣溫為13.1 ℃；年平均降水量為500～800 mm；人均水資源量為462 m3，為全國平均水平的五分之一，水資源承載能力與經(jīng)濟社會發(fā)展矛盾突出，是南水北調(diào)東線和中線工程的受水區(qū)[16]。黃淮海地區(qū)土層深厚，土質(zhì)肥沃，農(nóng)業(yè)規(guī)模大，主要種植小麥和玉米。地下水開采量的持續(xù)增加，導致一些地區(qū)的地下水位下降嚴重。

圖1 研究區(qū)位置Fig.1 The location of the study area

2 數(shù)據(jù)與方法

2.1 數(shù)據(jù)來源

（1）GRACE 重力衛(wèi)星數(shù)據(jù)。GRACE 計劃是美國國家航空航天局和德國航空中心的合作項目，通過衛(wèi)星監(jiān)測地球重力場，估算陸地水儲量的變化。本文選用美國德克薩斯大學研究中心（CSR）最新發(fā)布的RL06 Mascon 版本（https://www2.csr.utexas.edu/grace/）提供的重力衛(wèi)星陸地水儲量變化數(shù)據(jù)。相比于傳統(tǒng)的球諧系數(shù)，Mascon 產(chǎn)品的最大特點是無需做任何后處理，可以直接提供格網(wǎng)產(chǎn)品。該數(shù)據(jù)以2004—2009 年的陸地水儲量等效水深為基準，空間分辨率為0.25°×0.25°，時間分辨率為1 個月，時間序列從2003 年2 月至2021 年12 月，個別缺測時段采用奇異譜分析法[17]插補缺失值。

（2）GLDAS 陸面同化數(shù)據(jù)。全球陸面數(shù)據(jù)同化系統(tǒng)（GLDAS）包含4 個陸面模型（Noah、CLM、Mosaic、VIC），本文采用的是GLDAS-2.1 Noah 模型陸面同化數(shù)據(jù)，可解決內(nèi)陸水域無效數(shù)據(jù)的問題。從中選取總冠層水量（Canoplnt）、雪水當量（SWE）、土壤含水量（SMS）和蒸散發(fā)（Evap）數(shù)據(jù)，其中SMS 為4 層土壤含水量數(shù)據(jù)之和（0～10 cm、10～40 cm、40～100 cm、100～200 cm）。時間長度及空間分辨率與GRACE 數(shù)據(jù)一致。

（2）降水和供水量數(shù)據(jù)。數(shù)據(jù)來源于水利部和流域部委水資源公報，時間序列為2003—2021 年。選取各流域降水、地表水供水量、地下水供水量、跨流域調(diào)水和其他水源等數(shù)據(jù)，其他水源包括污水處理回水、集雨工程供水、海水淡化供水等。

2.2 研究方法

2.2.1 地下水儲量變化估算 陸地水儲量包括地表水、地下水、土壤水、生物水、冰雪、植被冠層水等，其中地表水與生物水相對于其他成分，其多年變化量可以忽略不計[18-19]。因此，陸地水儲量變化可表示為：

由此可得到地下水儲量變化量的計算：

式中：ΔGWS為地下水儲量變化（cm）；ΔTWS為陸地水儲量變化（cm）；ΔSMS為土壤含水量變化（cm）；ΔSWE為雪水當量變化（cm）；ΔCanoplnt為冠層水量變化（cm）。陸地水儲量變化由GRACE 重力衛(wèi)星反演得到，土壤含水量、雪水當量和植被冠層水的變化由GLDAS 數(shù)據(jù)計算。

2.2.2 Theil-Sen Median 和Mann-Kendall 法 Theil-Sen Median 方法又稱Sen 斜率估計，是一種穩(wěn)健的非參數(shù)統(tǒng)計的趨勢計算方法[20]，用于估計地下水儲量的變化速率 β，計算公式為：

式中：x為 各時段的ΔGWS；i和j為時段月份數(shù)；Median 表示取中值。若β大于0，則表示地下水儲量為增長趨勢，反之則為下降趨勢。

Mann-Kendall 檢驗是一種非參數(shù)的時間序列趨勢性檢驗方法，廣泛應(yīng)用于水文氣象數(shù)據(jù)分析，定義檢驗統(tǒng)計量S：

使用檢驗統(tǒng)計量Z進行趨勢檢驗，Z值計算方法為：

式中：Svar為統(tǒng)計量S的方差。取顯著性水平α=0.05，當|Z|≥1.96 時，表示該序列通過95%置信度檢驗，地下水儲量變化趨勢顯著，反之變化趨勢不顯著。

2.2.3 隨機森林模型 隨機森林（Random Forest）是一種基于統(tǒng)計學習理論的組合分類智能算法?；締卧菦Q策樹，使用自抽樣（Bootstrap）法從原始數(shù)據(jù)集中隨機重復采樣，得到多個決策樹模型，并將所有結(jié)果的平均值作為最終回歸結(jié)果。根據(jù)基尼指數(shù)判斷多個自變量對因變量的所做貢獻，可以有效進行特征重要性排序，隨機森林算法具有訓練速度快、精度高、不易過擬合等優(yōu)點，在水文地質(zhì)學領(lǐng)域應(yīng)用越來越廣泛[21]。本文采用python 語言中sklearn 庫建立隨機森林模型，以地下水儲量變化作為因變量，將影響地下水儲量變化的因素作為自變量，計算影響因素的貢獻率。

基尼指數(shù)代表了模型的不純度，基尼指數(shù)越小，則不純度越低，特征越好。其計算公式為：

式中：K為數(shù)據(jù)集分類總數(shù)；pk為某個子集屬于第k類的概率。

特征A在內(nèi)部某個節(jié)點m的重要性Vim,m為分枝前后基尼指數(shù)的變化量：

式中：Gini,m為節(jié)點m處的基尼指數(shù)，Gini,l、Gini,r為分枝后出現(xiàn)新節(jié)點處的基尼指數(shù)；l、r為節(jié)點m處分枝后的新節(jié)點。

則特征A在隨機森林模型中基尼指數(shù)評分為：

式中：Vim,i為第i棵決策樹各節(jié)點處分支前后基尼指數(shù)變化量，n為節(jié)點數(shù)。

3 結(jié)果與分析

3.1 時間變化特征

根據(jù)GRACE 和GLDAS 數(shù)據(jù)，采用式（2）計算研究區(qū)各網(wǎng)格地下水儲量的逐月變化。整個黃淮海地區(qū)相對于基準期，地下水儲量在2003—2021 年總體波動范圍為?24.94～5.99 cm（見圖2），整體呈減少趨勢，近20 年內(nèi)共減少23.18 cm，平均每年減少1.22 cm。其中在2003—2010 年平均減少0.45 cm/a，在2011—2019 年間減少速度加快為1.68 cm/a，其原因可能與經(jīng)濟快速發(fā)展及人口增長有關(guān)，在2020 年以后以1.05 cm/a 的速度開始增加。

圖2 黃淮海地區(qū)地下水儲量時間變化趨勢Fig.2 Temporal fluctuations in groundwater storage within the Huang-Huai-Hai region

逐一分析比較黃淮海地區(qū)內(nèi)各二級子流域地下水儲量變化（見圖3），具體變化速度見表1。從表1 和圖3 可見研究期內(nèi)各子流域的地下水儲量均處于持續(xù)減少狀態(tài)，海河流域平均下降速率為1.75 cm/a，淮河流域平均下降速率為0.57 cm/a。在2003—2010 年間，地下水儲量下降最快的是海河南系，達1.94 cm/a，而淮河流域整體處于緩慢上升階段。2011—2019 年研究區(qū)地下水儲量虧損程度擴大，虧損地區(qū)增多，除淮河下游變化不明顯，各子流域呈現(xiàn)不同程度的下降趨勢，花園口以下地區(qū)以2.92 cm/a 的速率下降，海河南系下降速率為2.86 cm/a。2020 年后，部分子流域地下水儲量開始上升，花園口以下地區(qū)上升速率為7.22 cm/a。

表1 黃淮海地區(qū)子流域地下水儲量變化速度Tab.1 The rates of change in groundwater storage within the Huang-Huai-Hai region

圖3 黃淮海地區(qū)各子流域地下水儲量變化趨勢Fig.3 Temporal fluctuations in groundwater storage within the sub-basins of the Huang-Huai-Hai region

各子流域地下水儲量變化在年內(nèi)變化表現(xiàn)出差異性（圖4）。灤河及冀東沿海、海河北系地區(qū)地下水儲量逐月變化較為穩(wěn)定，山東半島沿海諸島、淮河下游地區(qū)地下水儲量在4—8 月份處于高位。海河南系、徒駭馬頰河、花園口以下、淮河上游等主要農(nóng)作物耕種地區(qū)的地下水儲量下降明顯，且在6—7 月份明顯偏少。黃淮海地區(qū)水儲量持續(xù)減少主要是因農(nóng)業(yè)灌溉進行的大規(guī)模開采地下水[22]，用于農(nóng)作物灌溉的地下水供水量占地下水開采總量的80%。黃淮海地區(qū)的農(nóng)作物以冬小麥和夏玉米為主，5—8 月為作物需水關(guān)鍵期，需水量占整個生長發(fā)育過程的79%[23]，地下水下降最快，地下水儲量處于最低狀態(tài)。冬小麥開始大量播種時間為10 月中旬，越冬期不需要灌溉，降水基本能滿足需要，此時地下水儲量得到恢復。提高農(nóng)業(yè)用水效率和減少地下水開采量是緩解區(qū)域地下水資源短缺的重要措施。

圖4 黃淮海地區(qū)地下水儲量逐月變化Fig.4 The monthly average of groundwater storage in the Huang-Huai-Hai region

3.2 空間變化趨勢

從空間分布看，研究區(qū)地下水儲量變化存在明顯的空間差異性，整體上由東至西呈現(xiàn)下降幅度增大趨勢，變化范圍為?21.9～2.27 cm，負值為相對于基準期處于減少狀態(tài)。圖5 表明地下水儲量減幅最大的是海河南系，主要集中在邯鄲、安陽、鶴壁，淮河流域上游較為盈余；最低值出現(xiàn)在鶴壁市，最高值出現(xiàn)在淮安市。位于海河南系的白洋淀流域是華北最大的淡水湖，已成為雄安新區(qū)發(fā)展的重要生態(tài)水體，該地區(qū)地下水儲量最低值為?10.75 cm，最高值為?4.07 cm，虧損較為嚴重。

圖5 黃淮海地區(qū)多年平均地下水儲量變化分布Fig.5 The spatial pattern of groundwater storage within the Huang-Huai-Hai region

基于Sen 斜率估計和Mann-Kendall 法分析黃淮海地區(qū)2003—2021 年地下水儲量變化趨勢，計算出每個格網(wǎng)的β斜率和顯著性Z值，研究區(qū)88.87%區(qū)域通過置信水平0.05 的檢驗，其中94.2%的區(qū)域呈現(xiàn)出顯著的負向變化，即地下水儲量減少。為了進一步比較各個地區(qū)地下水儲量的變化情況，將研究區(qū)域劃分為極顯著增加、顯著增加、微顯著增加、不顯著增加、不顯著減少、微顯著減少、顯著減少、極顯著減少等8 個等級（圖6）。從圖6 可以看出：黃淮海地區(qū)地下水儲量減少速率整體表現(xiàn)出自南向北逐步增大的特點，海河南系大部分區(qū)域處于極顯著減少狀態(tài)，水儲量增加的地區(qū)位于淮河流域上游及中游部分區(qū)域。

圖6 Sen+M-K 計算地下水儲量變化趨勢分布Fig.6 The trend of variations in groundwater storage calculated using Sen+M-K methods

3.3 地下水儲量變化成因

3.3.1 降水、蒸散發(fā)變化影響 降水入滲補給是影響地下水儲量變化最直接的因子?；春恿饔蚪邓S沛，年均降水量為873.3 mm，充足的降水可促進地下水補給，有利于恢復地下水儲量。海河流域年均降水量為532.7 mm，降雨徑流經(jīng)水庫截留等損失后流域下滲較少，對地下水儲量變化影響的并不明顯（圖7）。同時，由于海河流域農(nóng)業(yè)部門耗水最多，降水量變化會影響灌溉量，枯水年灌溉量增多，地下水開采量會隨之增多。蒸散發(fā)過程會消耗土壤水分，減少地下水的補給。黃淮海地區(qū)多年平均蒸散量表現(xiàn)為增加趨勢，淮河流域的蒸散量由2003—2015 年的722.4 mm 增加至2016—2021 年的782.7 mm，海河流域的蒸散量由533.9 mm 增加至597.3 mm，花園口以下的蒸散量由629.7 mm 增加至678.0 mm。

圖7 降水、蒸散發(fā)與地下水儲量變化關(guān)系Fig.7 The interrelationship among precipitation, evapotranspiration, and groundwater storage

降水在黃淮海地區(qū)時間空間上的分布差異明顯?；春恿饔蚪邓觌H變化更大，直接影響地下水儲量變化速度，地下水儲量變化與降水變化趨勢具有較明顯的一致性。而海河流域降水對地下水儲量變化的影響不顯著，整個海河流域年降水量與地下水儲量之間的相關(guān)系數(shù)僅為0.02，灤河及冀東沿海、海河北系甚至出現(xiàn)負相關(guān)，這表明影響海河流域地下水儲量變化的主要因素不是降水量，而是人類開采的地下水量，這與之前學者研究結(jié)果一致[24]。

3.3.2 供水量變化影響 流域供水來源包括當?shù)氐乇硭⒌叵滤?、跨流域調(diào)水等，其中淮河流域外調(diào)水為引長江水與黃河水，海河流域外調(diào)水在2015 年之前為引黃河水，南水北調(diào)引長江水主要從2015 年開始?；春恿饔虻叵滤┧壤鹉昃徛郎p少，從2003 年占總供水量的32.3%降至2021 年的22%，主要原因是跨流域調(diào)水量呈增加趨勢；地表水多年平均供水比例為55.8%，是其最主要的供水來源，地下水次之，平均占比26.9%（圖8）。黃河花園口以下地區(qū)地下水供水比例呈減小趨勢，但總供水量呈增加趨勢。

圖8 黃淮海地區(qū)年供水量變化Fig.8 The annual fluctuations in water supply within the Huang-Huai-Hai region

海河流域地下水是主要的供水來源，2015 年以前地下水供水占比64.2%，2015 年南水北調(diào)工程通水后開始顯著減少，2019—2021 年已降到39.1%，但仍是海河流域最主要的供水來源。南水北調(diào)工程通水后，海河流域供水結(jié)構(gòu)發(fā)生了改變，地下水供水量明顯減少，但由于地下水多年持續(xù)超采嚴重，目前仍維持在近150 億m3/a，地下水儲量減少速率有所減緩，但下降趨勢未得到根本遏制[10]。

3.3.3 成因分析 以降水量、蒸散量、地下水供水量和跨流域調(diào)水量為自變量建立隨機森林模型，將計算所得影響因素重要性評分進行歸一化處理，定量分析比較不同因子對黃淮海地區(qū)地下水儲量變化的相對影響程度（圖9）。結(jié)果表明，淮河流域地下水儲量變化的影響因子相對貢獻率為降水＞蒸散發(fā)＞地下水供水＞跨流域調(diào)水，而海河流域影響因子的相對貢獻率為地下水供水＞跨流域調(diào)水＞降水＞蒸散發(fā)，黃河流域花園口以下的相對貢獻率為地下水供水＞蒸散發(fā)＞降水。各個流域地下水儲量變化的主導因素不同，降水是淮河流域的主導因素，貢獻率28.3%；地下水供水是海河流域和黃河流域花園口以下地區(qū)的主導因素，貢獻率分別達52.6%和58.0%。

圖9 黃淮海地區(qū)地下水儲量變化影響因素的相對貢獻度Fig.9 The relative contributions of factors influencing variations in groundwater storage within the Huang-Huai-Hai region

降水和蒸發(fā)等自然因素對淮河流域地下水儲量變化的影響程度顯著大于海河流域，而地下水供水是影響海河流域和黃河流域花園口以下地區(qū)地下水儲量下降的主導因素。研究表明人為因素對黃淮海地區(qū)地下水儲量變化影響逐漸增大[25-26]，2003—2021 年間，地下水消耗遠遠超出了地下水系統(tǒng)的自然補給能力，導致整個海河流域和淮河流域北部地下水儲量大幅下降。海河流域跨流域調(diào)水的相對貢獻率大于淮河流域，隨著地下水超采治理及南水北調(diào)供水工程持續(xù)推進，地下水儲量下降速率近年來有所減緩，北京、天津等受水區(qū)城市地下水位上升明顯[27-28]，但農(nóng)業(yè)區(qū)地下水位仍在下降。

4 結(jié) 語

（1）2003—2021 年間，黃淮海地區(qū)地下水儲量呈波動減少趨勢，平均變化速率為?1.22 cm/a，其中海河流域為?1.75 cm/a，淮河流域為?0.57 cm/a；海河南系虧損程度最嚴重，淮河上游和下游變化不顯著。2010 年代的下降速率高于2000 年代，隨著南水北調(diào)工程調(diào)水量的增加，地下水儲量下降速率有所減緩，沿線部分城市地下水儲量有所回升。

（2）地下水儲量變化存在明顯的空間差異性，海河流域虧損程度大于淮河流域。黃淮海地區(qū)94.2%的區(qū)域地下水儲量呈下降狀態(tài)，下降速率整體表現(xiàn)出自南向北逐步增加的特征；海河南系為極顯著減少趨勢，淮河流域上游及中游部分地區(qū)呈增加趨勢。

（3）黃淮海地區(qū)地下水儲量變化的主導因素不同，降水對淮河流域地下水儲量變化具有主導作用，而影響海河流域和黃河流域花園口以下地區(qū)最主要的因素是地下水供水量，多年持續(xù)超采造成地下水儲量大幅減少。南水北調(diào)工程通水后海河流域地下水開采量減少，改善了區(qū)域供水結(jié)構(gòu)，地下水儲量下降得到緩解，但還需進一步提高農(nóng)業(yè)用水效率，加強地下水資源管理。