中國東北三省地下水儲(chǔ)量時(shí)空變化特征及其影響因素分析

王子龍 孫昌鴻 姜秋香 劉傳興 單家珣

摘要：為完善區(qū)域地下水開發(fā)利用措施、規(guī)劃區(qū)域地下水資源管理，利用GRACE衛(wèi)星評估2002—2017年中國東北三省地下水儲(chǔ)量變化規(guī)律。結(jié)合GRACE和GLDAS估算地下水儲(chǔ)量變化，與實(shí)測地下水儲(chǔ)量變化對比驗(yàn)證，并探究其影響因素。結(jié)果表明：GRACE模擬地下水儲(chǔ)量變化與實(shí)測地下水儲(chǔ)量變化相關(guān)性較強(qiáng)，為0.72;地下水儲(chǔ)量在2013年盈余最大，2008年虧損最大，平均增長率為2.23 mm/a，秋冬兩季有明顯虧損，夏季發(fā)生盈余;地下水儲(chǔ)量空間分布有明顯差異性，2013年前東北少西南多，2013年后東北多西南少，黑龍江省變化較為明顯，遼寧省和吉林省受旱災(zāi)影響虧損過多;降水量和農(nóng)業(yè)用水量變化與地下水儲(chǔ)量變化極顯著相關(guān)，冬季地下水儲(chǔ)量變化與降雪顯著相關(guān)。研究東北三省地下水儲(chǔ)量時(shí)空變化對中國乃至全球水資源優(yōu)化配置和生態(tài)環(huán)境可持續(xù)發(fā)展具有參考價(jià)值。

關(guān)鍵詞：地下水儲(chǔ)量變化;GRACE;GLDAS;降水量;農(nóng)業(yè)用水

中圖分類號：P228;S127

文獻(xiàn)標(biāo)志碼：A

文章編號：1001-6791（2023）03-0360-14

收稿日期：2023-01-04;

網(wǎng)絡(luò)出版日期：2023-06-01

網(wǎng)絡(luò)出版地址：https：∥kns.cnki.net/kcms2/detail/32.1309.P.20230531.1450.004.html

基金項(xiàng)目：國家自然科學(xué)基金資助項(xiàng)目（52179035）;黑龍江省自然科學(xué)基金資助項(xiàng)目（YQ2019E004）

作者簡介：王子龍（1982—），男，山東膠州人，教授，博士，主要從事寒區(qū)水土資源高效利用方面研究。

E-mail：wangzilong2017@126.com

通信作者：姜秋香，E-mail：jiangqiuxiang2017@163.com

地下水作為水資源的重要組成部分，影響著作物生長、土壤質(zhì)地以及生態(tài)環(huán)境平衡［1］。地下水儲(chǔ)量作為衡量水資源量的重要評價(jià)指標(biāo)之一，探究其時(shí)空變化情況及影響因素已成為亟須解決的科學(xué)問題［2］。此外，地下水儲(chǔ)量的觀測可用于地下水文研究和水文災(zāi)害預(yù)警［3］。因此，研究地下水儲(chǔ)量變化特征對探究東北三省水儲(chǔ)量時(shí)空分布及理解跨流域大尺度水資源優(yōu)化配置有重要意義。

目前，相比于傳統(tǒng)的地下水儲(chǔ)量監(jiān)測方法，重力反演與氣候?qū)嶒?yàn)（Gravity Recovery and Climate Experiment，GRACE）衛(wèi)星技術(shù)打破了傳統(tǒng)地基觀測成本高、氣象站點(diǎn)監(jiān)測結(jié)果分布不均等局限性，解決了大尺度水文數(shù)據(jù)獲取難的問題，數(shù)據(jù)較為精確且監(jiān)測尺度統(tǒng)一。國內(nèi)外學(xué)者將GRACE衛(wèi)星廣泛應(yīng)用于全球、區(qū)域和流域陸地水儲(chǔ)量監(jiān)測，進(jìn)而計(jì)算地下水儲(chǔ)量［4-6］。其中，國外學(xué)者針對流域范圍的陸地水儲(chǔ)量研究集中于亞馬孫河流域和密西西比河流域，國內(nèi)學(xué)者則關(guān)注于長江、黃河、雅魯藏布江、黑河以及珠江等流域［7-10］。部分學(xué)者也通過監(jiān)測歐洲和北極地區(qū)的陸地水儲(chǔ)量來分析當(dāng)?shù)氐乃臍夂蛱匦裕?1-12］，開展對寒旱地區(qū)的地下水儲(chǔ)量估算研究，如印度、中國北方地區(qū)及華北平原地區(qū)等［13-15］。

中國東北三省包括黑龍江省、遼寧省和吉林省，總面積約78萬km2，分布范圍為38°36′N—53°36′N，118°36′E—135°06′E，屬溫帶大陸性季風(fēng)氣候，是中國重要的商品糧基地。地下水資源總量相對豐富，但地下水供水量基本消耗地下水資源一半左右，年際變化不穩(wěn)定且空間分布不均，部分地區(qū)中度缺水［16］。因冬季寒冷，監(jiān)測地下水較為困難，長期存在數(shù)據(jù)稀缺問題，對評估地下水資源和識別地下水影響因素有一定的限制。對于東北三省水儲(chǔ)量的研究，大部分學(xué)者停留于利用GRACE衛(wèi)星監(jiān)測陸地水儲(chǔ)量，分析地下水儲(chǔ)量時(shí)空變化特征及其影響因素的研究仍然有限。

本文基于GRACE新一代數(shù)據(jù)產(chǎn)品RL06結(jié)合GLDAS水文模型反演地下水儲(chǔ)量，并與基于實(shí)測地下水?dāng)?shù)據(jù)所估算的儲(chǔ)量變化對比驗(yàn)證，通過離散小波變換法研究東北三省地下水儲(chǔ)量的變化趨勢，運(yùn)用Mann-Kendall突變檢驗(yàn)法分析地下水儲(chǔ)量的季節(jié)變化，借助相關(guān)分析法和頻譜圖研究自然因素和人類活動(dòng)對其變化的影響。

1?數(shù)據(jù)來源與處理方法

1.1?GRACE數(shù)據(jù)

GRACE衛(wèi)星由美國國家航空航天局（NASA）和德國宇航中心（German Aerospace Center，DLR）合作研發(fā)，是監(jiān)測并記錄重力場變化的重力衛(wèi)星［17］。主要為地球系統(tǒng)科學(xué)提供高分辨率平均時(shí)變地球重力場模型;利用GPS無線電掩星（Radio Occultation）技術(shù)獲取全球溫濕度剖面圖;為研究海洋、冰川冰蓋質(zhì)量變化和水儲(chǔ)量變化提供高精度的重力場模型［18］。

GRACE數(shù)據(jù)產(chǎn)品主要由美國得克薩斯大學(xué)空間研究中心（Center for Space Research，CSR）、德國地學(xué)研究中心（Geo Forschungs Zentrum，GFZ）和美國噴氣動(dòng)力實(shí)驗(yàn)室（Jet Propulsion Laboratory，JPL）3家機(jī)構(gòu)對外開放，可在ICGEM（International Centre for Global Earth Models）中心下載（http：∥icgem.gfz-potsdam.de/series）。本文采用GFZ中心提供的RL06月重力場數(shù)據(jù)，扣除高頻非潮汐大氣和海洋的質(zhì)量變化影響［19］，時(shí)間范圍為2002年4月至2017年6月，共163個(gè)月（部分月份缺失），空間分辨率為1.00°×1.00°。

對GRACE月重力場數(shù)據(jù)進(jìn)行預(yù)處理，運(yùn)用海洋模型和GRACE Stokes系數(shù)所估算的1階重力系數(shù)替換原有系數(shù);采用高精度衛(wèi)星激光測距（Satellite Laser Ranging，SLR）得到的月估計(jì)值替換原有的C20球諧系數(shù)?；赟wenson提出的去相關(guān)濾波器法去除GRACE數(shù)據(jù)南北條紋誤差，以半徑為300 km的高斯濾波平滑信號噪聲，濾波后的數(shù)據(jù)截取為60階次，由此計(jì)算出球諧系數(shù)，并去除沿海區(qū)域海平面信號泄露的影響。將球諧系數(shù)轉(zhuǎn)換為空間分辨率更高的0.25°×0.25°柵格數(shù)據(jù)，同時(shí)利用gmt_cs2grid函數(shù)估算數(shù)據(jù)處理過程中的偏差和泄漏從而降低誤差。上述處理過程基于馮偉設(shè)計(jì)的開源程序GRACE_Matlab_Toolbox（GRAMAT）［20］，完善部分代碼后計(jì)算陸地水儲(chǔ)量為等效水高度，單位為mm。

1.2?GLDAS數(shù)據(jù)

GLDAS（Global Land Data Assimilation Systems）是由美國國家航空航天局戈達(dá)德空間飛行中心（GSFC）與美國海洋和大氣局國家環(huán)境預(yù)報(bào)中心（NCEP）聯(lián)合研發(fā)提供［21］，數(shù)據(jù)可在美國宇航局戈達(dá)德地球科學(xué)數(shù)據(jù)和信息服務(wù)中心GES DISC下載（https：∥disc.gsfc.nasa.gov）。GLDAS數(shù)據(jù)集中包含CLSM和NOAH 2個(gè)陸面模型和VIC水文模型。本研究從GLDAS-2 NOAH模型數(shù)據(jù)中提取出逐月土壤水分?jǐn)?shù)據(jù)，空間分辨率為0.25°×0.25°，時(shí)間序列為2002—2017年。

為削減GRACE后處理過程中信號和精度改變引起的誤差，使其與真實(shí)陸地水儲(chǔ)量近乎相似。對數(shù)據(jù)進(jìn)行2種處理，一是將GLDAS格網(wǎng)數(shù)據(jù)進(jìn)行同GRACE處理方法相同的濾波平滑處理，二是對0～2 m土層的土壤水分求和并距平。比較2種處理后的GLDAS模型數(shù)據(jù)，算得尺度因子可以反映出GRACE信號衰減幅度，其中距平處理得到的土壤水分變化即為土壤水儲(chǔ)量變化。本文采用單一尺度因子恢復(fù)法，2種處理方法得到的GLDAS數(shù)據(jù)之間殘差平方和最小時(shí)，計(jì)算得尺度因子為1.94。將其與GRACE估算的陸地水儲(chǔ)量相乘，得到最終的陸地水儲(chǔ)量變化估計(jì)值［22］。

1.3?實(shí)測地下水和影響因素?cái)?shù)據(jù)

本文根據(jù)《中國地質(zhì)環(huán)境監(jiān)測地下水位年鑒》獲取東北三省的地下水位實(shí)測數(shù)據(jù)，時(shí)間范圍從2005—2015年共計(jì)119個(gè)月（部分月份缺失）。研究區(qū)實(shí)測地下水類型選取為潛水和承壓水，選取的地下水位數(shù)據(jù)為地下水水位變幅。地下水實(shí)際監(jiān)測點(diǎn)位于東北三省36個(gè)城市，共計(jì)56個(gè)監(jiān)測點(diǎn)，基本覆蓋整個(gè)研究區(qū)，部分監(jiān)測點(diǎn)地下水位實(shí)測值如表1所示，研究區(qū)地下水位監(jiān)測點(diǎn)位置分布如圖1所示。

降水?dāng)?shù)據(jù)來源于中國氣象數(shù)據(jù)網(wǎng)（http：∥data.cma.cn/）和國家青藏高原科學(xué)數(shù)據(jù)中心（http：∥data.tpdc.ac.cn/zh-hans/），2種數(shù)據(jù)的精度同為0.1 mm。本研究選取了位于黑龍江省的氣象站降水量逐日變化資料，將其整理匯總為月度數(shù)據(jù)?；谳^高分辨率的《中國1 km分辨率逐月降水量數(shù)據(jù)集（1901—2020）》獲取吉林省和遼寧省的降水?dāng)?shù)據(jù)［23-24］。根據(jù)降水?dāng)?shù)據(jù)集計(jì)算得黑龍江省的月值降水量與氣象站點(diǎn)的降水量較為相符，可用于計(jì)算東北三省的月度降水?dāng)?shù)據(jù)。

潛在蒸散發(fā)數(shù)據(jù)來源于國家青藏高原科學(xué)數(shù)據(jù)中心，本研究根據(jù)《中國1 km逐月潛在蒸散發(fā)數(shù)據(jù)集（1990—2020）》獲取黑龍江省、吉林省和遼寧省的潛在蒸散發(fā)數(shù)據(jù)［25-26］，精度為0.1 mm。

地表徑流數(shù)據(jù)來自GLDAS的CLSM陸面模型，時(shí)間分辨率為1 d，空間分辨率為0.25°×0.25°，時(shí)間范圍為2003—2017年。將其整合處理為月值數(shù)據(jù)后，計(jì)算得地表徑流量。

人口數(shù)據(jù)來源于國家青藏高原科學(xué)數(shù)據(jù)中心的《全國各地區(qū)常住人口規(guī)模（2007—2018）》［27］和《全國各地區(qū)人口出生率、死亡率、自然增長率（2001—2008）》［28］。農(nóng)業(yè)用水量數(shù)據(jù)來源于黑龍江省、吉林省和遼寧省的水資源公報(bào)。

1.4?數(shù)據(jù)處理方法

1.4.1?水量平衡原理

陸地水儲(chǔ)量包含有地表徑流量、地下水量、土壤水量、雪水儲(chǔ)量及植被冠層儲(chǔ)水量［29］。由于GLDAS監(jiān)測植被含水量和冰川積雪儲(chǔ)量變化量相對于陸地總水儲(chǔ)量變化量微小，參考已有研究在計(jì)算水儲(chǔ)量變化時(shí)將其忽略。同時(shí)地表水受水循環(huán)的作用處于相對穩(wěn)定狀態(tài)，探究水儲(chǔ)量變化時(shí)也基本可以忽略［30-31］。由此，構(gòu)建水量平衡方程如下：

式中：ΔSGW為地下水儲(chǔ)量變化，mm;ΔSTW為陸地水儲(chǔ)量變化，mm;ΔSSM為土壤水儲(chǔ)量變化，mm。

為便于驗(yàn)證，需要給實(shí)測地下水?dāng)?shù)據(jù)乘以給水度，給水度反映的是含水層釋水能力的強(qiáng)弱［32］。地下水實(shí)際監(jiān)測點(diǎn)大部分位于松嫩平原、遼河平原和白城扇形地，含水層類型為松散巖類孔隙含水層，結(jié)合中國地質(zhì)調(diào)查局所提供的含水層巖性給水度經(jīng)驗(yàn)值，計(jì)算得區(qū)域平均給水度為0.26［33-34］。將觀測的地下水位變幅插值平均，得到東北三省地下水位的月度數(shù)據(jù)。計(jì)算實(shí)測地下水儲(chǔ)量變化量公式如下：

式中：ΔSGW，me為實(shí)測地下水儲(chǔ)量變化，mm;H為地下水水位變幅月值，m。

1.4.2?評價(jià)方法

為評估模型的準(zhǔn)確性，運(yùn)用了皮爾遜相關(guān)系數(shù)（r）、均方根誤差（ERMS）和標(biāo)準(zhǔn)化均方根誤差（ENRMS）3種評價(jià)指標(biāo)。r的絕對值越接近于1說明相關(guān)性越好;ERMS越小表示數(shù)據(jù)一致性越好;ENRMS有助于比較不同比例的數(shù)據(jù)集或模型，ENRMS越接近或小于10%時(shí)，模擬值和實(shí)測值的一致性較好，但當(dāng)其大于30%時(shí)，模擬值和實(shí)測值一致性較差［35］。計(jì)算公式如下：

式中：X為實(shí)測值;Y為模擬值;Ymax、Ymin分別為模擬值的最大值和最小值；σX、σY分別為樣本X和Y的標(biāo)準(zhǔn)差;n為樣本個(gè)數(shù)。

1.4.3?Mann-Kendall突變檢驗(yàn)法

Mann-Kendall（MK）檢驗(yàn)是一種不受少數(shù)異常值影響的非參數(shù)統(tǒng)計(jì)檢驗(yàn)，可用于判斷水文氣象數(shù)據(jù)的突變年份和變化趨勢。通常給定顯著性水平α=0.05，相應(yīng)臨界值U0.05=±1.96，UF和UB分別是按照時(shí)間序列正順序和逆順序計(jì)算出的統(tǒng)計(jì)序列。UF和UB在臨界區(qū)內(nèi)交點(diǎn)所對應(yīng)的時(shí)間即為原時(shí)間序列突變開始的時(shí)間，當(dāng)UF>0時(shí)序列呈上升趨勢，反之呈下降趨勢，位于臨界值外時(shí)表示變化趨勢顯著［36］。

1.4.4?離散小波變換法

離散小波變換法（Discrete wavelet transform，DWT）是對連續(xù)小波變換的尺度、位移按2次方進(jìn)行離散化得到的，也稱為二進(jìn)制小波變換，其中高頻對應(yīng)高時(shí)間分辨率，低頻對應(yīng)高頻率精度，可在時(shí)間和頻率上取得一定平衡。將時(shí)間序列分解為低頻和高頻部分，能有效判斷時(shí)間序列中的潛在趨勢，其中最大分解水平上的低頻序列可用來表示時(shí)間序列的變化趨勢［37］。

2?結(jié)果分析與討論

2.1?地下水儲(chǔ)量空間分布

對于東北三省地區(qū)有部分國內(nèi)學(xué)者使用CSR和JPL機(jī)構(gòu)提供的GRACE數(shù)據(jù)模擬陸地水儲(chǔ)量變化，為保證準(zhǔn)確性和更高的精度，本文獲取CSR和JPL機(jī)構(gòu)提供的數(shù)據(jù)，采取相同的處理方式估算地下水儲(chǔ)量變化。分析其與實(shí)測地下水量的關(guān)系，實(shí)測地下水量與GFZ、CSR和JPL的ERMS分別為23.7 mm/月、24.8 mm/月和24.0 mm/月，其在研究區(qū)內(nèi)的精度較好，且GFZ的信噪比值大于其他2種產(chǎn)品，由此可知GFZ產(chǎn)品更適用于東北三省的陸地水儲(chǔ)量監(jiān)測。

補(bǔ)充GRAMAT已有的代碼程序，計(jì)算當(dāng)年12個(gè)月數(shù)據(jù)的平均值，以此來代表該年水儲(chǔ)量，再將16 a的年水儲(chǔ)量距平，得到東北三省地下水儲(chǔ)量年變化的空間分布，并運(yùn)用克里金法進(jìn)行空間插值如圖2所示。地下水儲(chǔ)量變化呈現(xiàn)出明顯波動(dòng)趨勢，2002—2008年有逐漸下降趨勢，2009年短暫上升后出現(xiàn)虧損狀態(tài)，于2013年增加到最大后緩慢下降。整體空間上地下水儲(chǔ)量在2013年前表現(xiàn)為南多北少，2013年后表現(xiàn)為南少北多。

2002—2012年東北部黑龍江省一帶地下水儲(chǔ)量減少，虧損最大區(qū)域位于小興安嶺地區(qū)、松嫩平原、三江平原以及東部山地一帶；大興安嶺地區(qū)地下水儲(chǔ)量相對穩(wěn)定，對比多年平均長期處于弱增長狀態(tài)；遼寧省、吉林省大部分地區(qū)和長白山地一帶地下水儲(chǔ)量變化相對較小，處于長期穩(wěn)定狀態(tài)。2013年研究區(qū)地下水儲(chǔ)量達(dá)到峰值，在黑龍江省西北部增長最大，與該年東北三省洪災(zāi)受損嚴(yán)重地區(qū)較為吻合。2013年后遼寧省和吉林省部分地區(qū)受旱災(zāi)影響地下水儲(chǔ)量開始處于虧損狀態(tài)且長期難以恢復(fù)至正常水平，虧損趨勢逐漸向黑龍江省擴(kuò)散；長白山地一帶受洪災(zāi)影響較小，地下水儲(chǔ)量在2013年短暫增長后重新歸于穩(wěn)定，僅在西南段有一定程度的虧損；黑龍江省大興安嶺地區(qū)地下水儲(chǔ)量在2013—2017年間雖然存在相對緩慢的減少趨勢，但較多年平均仍有一定盈余。

2.2?地下水儲(chǔ)量時(shí)間變化分析

分析地下水儲(chǔ)量變化量與實(shí)測地下水儲(chǔ)量變化量的相關(guān)性（圖3（a））。2002—2017年地下水儲(chǔ)量與實(shí)測地下水儲(chǔ)量變化量具有極顯著相關(guān)性且相關(guān)系數(shù)為0.72，表明GRACE模擬的地下水儲(chǔ)量變化可以很好地解釋東北三省地下水變化。GRACE模擬地下水儲(chǔ)量平均增長率為2.23 mm/a，實(shí)測地下水位變化也在逐年增長，平均增長率為3.19 mm/a。地下水儲(chǔ)量變化具有時(shí)間差異性，在2013年底達(dá)到最大，為72.03 mm;在2008年初最小，為-79.36 mm;2012—2016年整體變化幅度相對多年平均有所增加。實(shí)測地下水儲(chǔ)量變化幅度在-36～56 mm之間。根據(jù)圖3（b）發(fā)現(xiàn)，實(shí)測地下水儲(chǔ)量變化同ΔSSM具有一定程度的負(fù)相關(guān)，滿足水量平衡原理，當(dāng)某一時(shí)段土壤水儲(chǔ)量減少時(shí)，實(shí)測地下水儲(chǔ)量變化量處于增長狀態(tài)。實(shí)測地下水儲(chǔ)量變化同陸地水儲(chǔ)量變化的趨勢和波動(dòng)更為一致，分析兩者相關(guān)性，其相關(guān)系數(shù)為0.64，且呈極顯著相關(guān)關(guān)系，說明根據(jù)GRACE模擬的陸地水儲(chǔ)量中地下水儲(chǔ)量占據(jù)其重要組成部分。

對地下水儲(chǔ)量變化的時(shí)間序列進(jìn)行MK突變分析（圖3（c））和DWT趨勢分析（圖3（d），D5、D6分別表示分解層數(shù)為5和6的趨勢）。由MK突變檢驗(yàn)分析表明，地下水儲(chǔ)量大部分時(shí)間呈下降趨勢，在2011年底發(fā)生突變后呈現(xiàn)上升趨勢并在2014年初顯著上升。為更好地觀察地下水儲(chǔ)量變化的趨勢，借助DWT將地下水儲(chǔ)量序列分解為6個(gè)子序列，其中分解層數(shù)為6的低頻重構(gòu)序列反映的變化趨勢同MK結(jié)果更為相符，趨勢更加顯著，由此可確定東北三省2002—2017年地下水儲(chǔ)量總體呈波動(dòng)趨勢，先下降后上升最終趨于平穩(wěn)。

分離出2002—2017年間東北三省季節(jié)性地下水儲(chǔ)量變化，并通過MK法分析其趨勢和突變時(shí)段（圖4）。地下水儲(chǔ)量呈現(xiàn)明顯的季節(jié)性變化，振幅基本為±25 mm，與實(shí)測地下水儲(chǔ)量有較好的一致性。春季ΔSGW在2002—2012年間波動(dòng)較為明顯并存在一定虧損，于2012年發(fā)生突變后呈現(xiàn)上升趨勢;夏季ΔSGW趨勢明晰，在2002—2010年間持續(xù)下降，于2010年發(fā)生突變后逐漸呈上升趨勢，并在8月產(chǎn)生最大盈余，盈余24.83 mm;秋冬兩季地下水儲(chǔ)量在2002—2012年呈現(xiàn)持續(xù)下降趨勢，其中秋季ΔSGW存在短期增加，冬季2月份發(fā)生最大虧損，虧損-19.04 mm，兩季于2012年發(fā)生突變后逐漸恢復(fù)上升趨勢，推測地下水儲(chǔ)量減少的主要原因?yàn)樵摷竟?jié)降水較少，而生產(chǎn)生活供水較大，需要地下水供給。

由圖3（c）和圖4可知，地下水儲(chǔ)量在年際、春季、秋季和冬季的突變點(diǎn)均在2012年左右。根據(jù)中國氣象局發(fā)布的氣候事件，春季發(fā)生突變是由于黑龍江省該年春季土壤干旱程度最為嚴(yán)重，土壤濕度存在明顯下降趨勢；秋冬兩季發(fā)生突變是由于秋季多個(gè)臺(tái)風(fēng)登陸引起的暴雨天氣和冬季大范圍暴雪天氣，陸地水儲(chǔ)量呈持續(xù)上升趨勢。

2.3?地下水儲(chǔ)量變化的影響因素探究

地下水參與水循環(huán)路徑會(huì)受到降水入滲、地表徑流入滲補(bǔ)給、河道排泄和潛水層蒸發(fā)等自然因素的影響，同時(shí)灌溉補(bǔ)給和地下水開采等人類活動(dòng)也會(huì)對其造成影響。東三省被認(rèn)為是農(nóng)業(yè)重點(diǎn)發(fā)展區(qū)，其地下水作為農(nóng)業(yè)用水中的重要部分會(huì)受到一定程度影響［38-39］。

運(yùn)用皮爾遜相關(guān)分析法，分析自然因素和人為因素對地下水儲(chǔ)量的影響。自然因素包括降水、潛在蒸散發(fā)、土壤濕度和地表徑流，人為因素有人口數(shù)量和農(nóng)業(yè)用水。由于人為因素缺乏相關(guān)月度數(shù)據(jù)，因此只比對自然因素與ΔSGW的月相關(guān)系數(shù)。由圖5可知，自然因素與ΔSGW均存在極顯著相關(guān)，其中，呈明顯負(fù)相關(guān)的是土壤濕度，相關(guān)系數(shù)為-0.55，地表徑流的相關(guān)性較小為-0.34；呈較明顯正相關(guān)的是降水，相關(guān)系數(shù)為0.33。土壤濕度、降水和地表徑流成為對ΔSGW影響較大的因素。在年相關(guān)關(guān)系中，人口數(shù)量與地下水儲(chǔ)量變化量相關(guān)性不顯著，但農(nóng)業(yè)用水與ΔSGW顯著相關(guān)，相關(guān)系數(shù)為0.58。

本文在探究地下水開采對地下水儲(chǔ)量變化量的響應(yīng)關(guān)系時(shí)，以2002—2011年松遼流域水資源公報(bào)所提供數(shù)據(jù)為例（其大部分區(qū)域均與研究區(qū)重合），分析發(fā)現(xiàn)松遼流域地下水開采量以0.58 mm/a的趨勢增加，東北三省地下水儲(chǔ)量變化量在同期雖處于下降趨勢，但二者p值遠(yuǎn)大于0.05，無顯著性差異。而且由于地下水系統(tǒng)的復(fù)雜性和非線性特征也會(huì)影響地下水開采對地下水儲(chǔ)量變化量的響應(yīng)關(guān)系，因此暫未考慮地下水開采對地下水儲(chǔ)量變化量的影響。相對于地下水開采，農(nóng)業(yè)用水與地下水儲(chǔ)量變化量相關(guān)性最為明顯，根據(jù)黑龍江省水資源公報(bào)可知，當(dāng)?shù)叵滤畠?chǔ)量盈余時(shí)大部分被用于灌溉農(nóng)田。

基于偏相關(guān)分析法研究冬季地下水儲(chǔ)量變化與自然因素的相關(guān)關(guān)系，降水與其存在顯著相關(guān)性，成為氣象因子中的主導(dǎo)因素。探究土壤濕度、潛在蒸散發(fā)、地表徑流、綜合變量（包括前3種變量）和無變量情況下降水和ΔSGW的相關(guān)關(guān)系（表2）。5種情況下降水和ΔSGW均表現(xiàn)出極顯著相關(guān)，且在排除地表徑流因素的干擾時(shí)，其相關(guān)系數(shù)達(dá)到最高。土壤濕度對其相關(guān)性沒有明顯影響，而排除潛在蒸散發(fā)和綜合因素的干擾后相關(guān)性明顯下降，就對ΔSGW的影響而言，冬季潛在蒸散發(fā)與降水存在拮抗關(guān)系。由于東北地區(qū)冬季降水基本以降雪形態(tài)出現(xiàn)，所以冬季地下水儲(chǔ)量變化在一定程度上受降雪變化影響。

如圖6所示，借助降水?dāng)?shù)據(jù)發(fā)現(xiàn)降雪和地下水儲(chǔ)量的變化趨勢較為一致，有明顯的周期性且年、半年周期一致，地下水儲(chǔ)量變化量的振幅與降雪振幅較為相符。變化特征均呈先下降后上升的趨勢，降雪在2012年底達(dá)到峰值，后續(xù)2013年底也有較大幅度增加，地下水儲(chǔ)量隨之在2013年底產(chǎn)生最大盈余。其中，2007年1月降雪量驟減現(xiàn)象并沒有引起當(dāng)月地下水儲(chǔ)量的改變，但其1—2月明顯地減少趨勢對整體地下水儲(chǔ)量變化產(chǎn)生了一定影響，即在2008年2月達(dá)最大虧損。

由圖7可知，地下水儲(chǔ)量、降水、地表徑流、潛在蒸散發(fā)和土壤濕度均呈現(xiàn)出很明顯的周期性變化，且各年、半年周期變化規(guī)律基本一致，地表徑流因缺少部分?jǐn)?shù)據(jù)在周期性上有些微差異。對地表徑流影響較大的降水量與潛在蒸散發(fā)的年振幅分別為259.6 mm和500.4 mm，相較于其他氣象因素的振幅差異明顯;地下水儲(chǔ)量年振幅與土壤濕度最為接近，在13～18 mm之間;地表徑流年振幅最小，為4 mm。

地表徑流在很大程度上受降水和潛在蒸散發(fā)的影響，尤其受潛在蒸散發(fā)影響產(chǎn)生較大損失。地下水儲(chǔ)量的降水補(bǔ)給部分，在很大程度上被地表徑流和潛在蒸散發(fā)消耗。ΔSGW、降水量、地表徑流、潛在蒸散發(fā)和土壤濕度年振幅的ERMS分別為2.60、35.66、0.61、63.92和1.84 mm/月，由此分析ΔSGW和氣象因子的ENRMS均在15.45%以下且差異性較小，表明ΔSGW和氣象因子的頻譜圖較為可信。

綜合相關(guān)系數(shù)和頻譜圖2種分析方法對比可知，自然因素中土壤濕度和降水量對ΔSGW的影響最明顯，又因本研究ΔSGW是依據(jù)土壤濕度計(jì)算而來，所以降水成為地下水儲(chǔ)量變化的主要影響因素。

2.4?降水量與地下水儲(chǔ)量的時(shí)空分布關(guān)系分析

計(jì)算得到的東北三省2002—2017年月度降水?dāng)?shù)據(jù)如圖8（a）所示，規(guī)律與趨勢同降水量的變化基本一致，2013年降水量增加，同期ΔSGW也有顯著上升。地下水儲(chǔ)量對降水的響應(yīng)存在一定延遲，不能及時(shí)反映季節(jié)性變化。在6—8月降水增加的背景下，地下水儲(chǔ)量則延遲至7—10月開始改變，相較于每年降水量的最大值和最小值ΔSGW都有少許延遲。

分析降水量年際變化對地下水儲(chǔ)量變化的影響。由圖8（b）可知，降水量均集中在年平均值425.6 mm以下，說明該時(shí)段產(chǎn)生明顯虧損，與同年ΔSGW有較好的一致性。降水累計(jì)3 a低于年均值，影響其對地下水的補(bǔ)給，從而在2008年初引起最大虧損。由圖8（c）可知，降水量于2011年集中在年均值偏下處，存在一定的虧損，但在2012—2013年降水量基本集中在年均值附近，尤其是2012年下半年產(chǎn)生明顯盈余，2012年和2013年上半年受異常值影響增加比較明顯，盈余量一定程度上彌補(bǔ)2011年的虧損量，同期ΔSGW也有較強(qiáng)程度的增幅。降水量突然增加給東北三省帶來了嚴(yán)重的洪水災(zāi)害，松花江流域發(fā)生了自1998年以來最大洪水。即使2014年降水量大幅減少，地下水儲(chǔ)量變化量也在多年平均值之上。

利用克里金插值方法繪制GRACE反演地下水儲(chǔ)量多年平均變化速率如圖9（a）所示，2002—2017年地下水儲(chǔ)量年均變化速率在空間上有明顯差異性，大部分地區(qū)基本穩(wěn)定，部分地區(qū)存在明顯的增減趨勢。在黑龍江省大興安嶺一帶增長速率最大，最高達(dá)68.48 mm/a；在小興安嶺南部、三江平原西北部、松嫩平原東北部和吉林省西北部等地降低速率較大，最低達(dá)-65.39 mm/a；遼寧省、吉林省中部地區(qū)以及黑龍江省部分地區(qū)年均盈余和虧損達(dá)到平衡，地下水儲(chǔ)量變化保持相對穩(wěn)定。

根據(jù)降水?dāng)?shù)據(jù)集整理出2002—2017年降水量年均值，采用克里金插值法繪制其空間分布如圖9（b）所示，研究時(shí)段內(nèi)東北三省年均降水量在220～821 mm之間，自沿海地區(qū)向內(nèi)陸地區(qū)逐漸減少。黑龍江省西部和吉林省西北部地區(qū)降水量最少，年均降水量在333 mm左右浮動(dòng)；遼河平原和長白山脈西南段降水量較多，年均降水量在670～820 mm之間。對比年均地下水儲(chǔ)量變化速率和降水量的空間分布，黑龍江西北部降水量相對較小，而地下水儲(chǔ)量的年均變化速率有增有減，推測是由于松嫩平原蒸發(fā)量大于降水量引起的干旱導(dǎo)致；黑龍江省大興安嶺山地即使降水較少地下水儲(chǔ)量依然存在增加趨勢，推測其土壤涵養(yǎng)水源能力良好且蒸發(fā)量少，可在一定程度上降低水土流失；遼河平原和長白山脈西南段，降水量最大，地下水儲(chǔ)量變化速率最穩(wěn)定，推測其受暴雨影響，侵蝕程度和地表徑流的提升速率相對東北三省其他地區(qū)更高。

3?結(jié)?論

本文基于GRACE衛(wèi)星數(shù)據(jù)估算東北三省2002—2017年地下水儲(chǔ)量變化，并與實(shí)測地下水儲(chǔ)量變化進(jìn)行對比驗(yàn)證，探究地下水儲(chǔ)量的時(shí)空變化規(guī)律，并分析影響地下水變化的主要因素。主要結(jié)論如下：

（1） 地下水儲(chǔ)量表現(xiàn)出多年波動(dòng)變化趨勢，其空間分布表現(xiàn)為2002—2012年東北少西南多，2013—2017年東北多西南少;東北三省地下水儲(chǔ)量變化在2011年底顯著增加，平均增長率為2.23 mm/a，在2013年底達(dá)最大盈余，并與實(shí)測值有極顯著相關(guān)性;地下水儲(chǔ)量秋冬兩季有明顯虧損，春季逐漸回升，夏季發(fā)生盈余。

（2） 降水和農(nóng)業(yè)用水是影響地下水儲(chǔ)量變化的主導(dǎo)因素。冬季降雪量變化密切影響地下水儲(chǔ)量變化，在時(shí)間分布上有明顯一致性?；谄嚓P(guān)分析法能有效探究冬季降雪對地下水儲(chǔ)量變化量的影響，但對于冬季潛在蒸散發(fā)對降雪和地下水儲(chǔ)量變化的間接影響，需要更加詳細(xì)的水文資料去分析。

（3） 雖然本文基于GRACE反演方法能有效估算東北三省地下水儲(chǔ)量變化，但相對于地下水監(jiān)測站仍無法提供更高分辨的監(jiān)測數(shù)據(jù)，隨著下一代衛(wèi)星GRACE Follow-on的發(fā)射有望改善其空間分辨率以拓寬重力衛(wèi)星在水文方面的應(yīng)用。

參考文獻(xiàn)：

［1］吳彬，杜明亮，穆振俠，等.1956—2016年新疆平原區(qū)地下水資源量變化及其影響因素分析［J］.水科學(xué)進(jìn)展，2021，32（5）：659-669.（WU B，DU M L，MU Z X，et al.Analysis on the variation of groundwater resources and influencing factors in Xinjiang plain area from 1956 to 2016［J］.Advances in Water Science，2021，32（5）：659-669.（in Chinese））

［2］劉國緯.水文科學(xué)的基本問題及當(dāng)代前沿［J］.水科學(xué)進(jìn)展，2020，31（5）：685-689.（LIU G W.Basic problems and frontier of hydrology［J］.Advances in Water Science，2020，31（5）：685-689.（in Chinese））

［3］尹文杰，胡立堂，王景瑞.基于GRACE重力衛(wèi)星的甘肅北山地區(qū)地下水儲(chǔ)量變化規(guī)律研究［J］.水文地質(zhì)工程地質(zhì)，2015，42（4）：29-34.（YIN W J，HU L T，WANG J R.Changes of groundwater storage variation based on GRACE data at the Beishan area，Gansu Province［J］.Hydrogeology & Engineering Geology，2015，42（4）：29-34.（in Chinese））

［4］SINGH A，REAGER J T，BEHRANGI A.Estimation of hydrological drought recovery based on precipitation and Gravity Recovery and Climate Experiment（GRACE） water storage deficit［J］.Hydrology and Earth System Sciences，2021，25（2）：511-526.

［5］羅志才，李瓊，鐘波.利用GRACE時(shí)變重力場反演黑河流域水儲(chǔ)量變化［J］.測繪學(xué)報(bào)，2012，41（5）：676-681.（LUO Z C，LI Q，ZHONG B.Water storage variations in Heihe River basin recovered from GRACE temporal gravity field［J］.Acta Geodaetica et Cartographica Sinica，2012，41（5）：676-681.（in Chinese））

［6］汪漢勝，王志勇，袁旭東，等.基于GRACE時(shí)變重力場的三峽水庫補(bǔ)給水系水儲(chǔ)量變化［J］.地球物理學(xué)報(bào)，2007，50（3）：730-736.（WANG H S，WANG Z Y，YUAN X D，et al.Water storage changes in Three Gorges water systems area inferred from GRACE time-variable gravity data［J］.Chinese Journal of Geophysics，2007，50（3）：730-736.（in Chinese））

［7］TAPLEY B D，BETTADPUR S，RIES J C，et al.GRACE measurements of mass variability in the earth system［J］.Science，2004，305（5683）：503-505.

［8］SWENSON S，WAHR J，MILLY P C D.Estimated accuracies of regional water storage variations inferred from the Gravity Recovery and Climate Experiment （GRACE）［J］.Water Resources Research，2003，39（8）：375-384.

［9］聶寧，張萬昌，陳豪，等.長江流域水儲(chǔ)量變化的時(shí)間變化特征及歸因分析［J］.水科學(xué)進(jìn)展，2021，32（3）：396-407.（NIE N，ZHANG W C，CHEN H，et al.Temporal variation characteristics and attribution analysis of terrestrial water storage change in the Yangtze River basin［J］.Advances in Water Science，2021，32（3）：396-407.（in Chinese））

［10］鄧梓鋒，吳旭樹，王兆禮，等.基于GRACE重力衛(wèi)星數(shù)據(jù)的珠江流域干旱監(jiān)測［J］.農(nóng)業(yè)工程學(xué)報(bào)，2020，36（20）：179-187.（DENG Z F，WU X S，WANG Z L，et al.Drought monitoring based on GRACE data in the Pearl River basin，China［J］.Transactions of the Chinese Society of Agricultural Engineering，2020，36（20）：179-187.（in Chinese））

［11］NIU G Y，SEO K W，YANG Z L，et al.Retrieving snow mass from GRACE terrestrial water storage change with a land surface model［J］.Geophysical Research Letters，2007，34（15）：L15704.

［12］ANDERSEN O B，SENEVIRATNE S I，HINDERER J，et al.GRACE-derived terrestrial water storage depletion associated with the 2003 European heat wave［J］.Geophysical Research Letters，2005，32（18）：109-127.

［13］RODELL M，VELICOGNA I，F(xiàn)AMIGLIETTI J S.Satellite-based estimates of groundwater depletion in India［J］.Nature，2009，460（7258）：999-1002.

［14］YIN W J，HU L T，JIAO J J.Evaluation of groundwater storage variations in Northern China using GRACE data［J］.Geofluids，2017，2017：1-13.

［15］HUANG Z Y，PAN Y，GONG H L，et al.Subregional-scale groundwater depletion detected by GRACE for both shallow and deep aquifers in North China Plain［J］.Geophysical Research Letters，2015，42（6）：1791-1799.

［16］張琦，李松森，夏慧琳.基于模糊綜合評判模型的東北三省水資源承載力研究［J］.水土保持通報(bào)，2019，39（5）：179-188，193.（ZHANG Q，LI S S，XIA H L.A study on water resource carrying capacity of three provinces in Northeast China based on fuzzy comprehensive evaluation model［J］.Bulletin of Soil and Water Conservation，2019，39（5）：179-188，193.（in Chinese））

［17］ZHENG W，HSU H，ZHONG M，et al.Efficient accuracy improvement of GRACE global gravitational field recovery using a new Inter-satellite Range Interpolation Method［J］.Journal of Geodynamics，2012，53：1-7.

［18］吳庭濤，鄭偉，尹文杰，等.地球衛(wèi)星重力場模型及其應(yīng)用研究進(jìn)展［J］.科學(xué)技術(shù)與工程，2020，20（25）：10117-10132.（WU T T，ZHENG W，YIN W J，et al.Review of earth′s gravity field model and application from satellite［J］.Science Technology and Engineering，2020，20（25）：10117-10132.（in Chinese））

［19］徐永明.利用GRACE RL06時(shí)變重力場的水儲(chǔ)量變化反演研究與應(yīng)用［D］.徐州：中國礦業(yè)大學(xué)，2020.（XU Y M.Research on the recovery method of water storage variations and applications using GRACE RL06 temporal gravity filed［D］.Xuzhou：China University of Mining and Technology，2020.（in Chinese））

［20］FENG W.GRAMAT：a comprehensive Matlab toolbox for estimating global mass variations from GRACE satellite data［J］.Earth Science Informatics，2019，12（3）：389-404.

［21］RODELL M，HOUSER P R，JAMBOR U，et al.The global land data assimilation system［J］.Bulletin of the American Meteorological Society，2004，85（3）：381-394.

［22］李婉秋，王偉，章傳銀，等.利用GRACE衛(wèi)星重力數(shù)據(jù)監(jiān)測關(guān)中地區(qū)地下水儲(chǔ)量變化［J］.地球物理學(xué)報(bào)，2018，61（6）：2237-2245.（LI W Q，WANG W，ZHANG C Y，et al.Monitoring groundwater storage variations in the Guanzhong area using GRACE satellite gravity data［J］.Chinese Journal of Geophysics，2018，61（6）：2237-2245.（in Chinese））

［23］彭守璋.中國1km分辨率逐月降水量數(shù)據(jù)集（1901—2020）［DB］.北京：國家青藏高原科學(xué)數(shù)據(jù)中心，2020.（PENG S Z.1-km monthly precipitation dataset for China （1901—2020）［DB］.Beijing：National Tibetan Plateau Data Center，2020.（in Chinese））

［24］PENG S Z，GANG C C，CAO Y，et al.Assessment of climate change trends over the Loess Plateau in China from 1901 to 2100［J］.International Journal of Climatology，2018，38（5）：2250-2264.

［25］彭守璋.中國1 km逐月潛在蒸散發(fā)數(shù)據(jù)集（1990—2021）［DB］.北京：國家青藏高原科學(xué)數(shù)據(jù)中心，2022.（PENG S Z.1 km monthly potential evapotranspiration dataset in China （1990—2021）［DB］.Beijing：National Tibetan Plateau Data Center，2022.（in Chinese））

［26］DING Y X，PENG S Z.Spatiotemporal trends and attribution of drought across China from 1901—2100［J］.Sustainability，2020，12（2）：477.

［27］青海省統(tǒng)計(jì)局.全國各地區(qū)常住人口規(guī)模（2007—2018）［DB］.北京：國家青藏高原科學(xué)數(shù)據(jù)中心，2021.（Qinghai Provincial Bureau of Statistics.Size of permanent resident population in different regions of China （2007—2018）［DB］.Beijing：National Tibetan Plateau Data Center，2021.（in Chinese））

［28］青海省統(tǒng)計(jì)局.全國各地區(qū)人口出生率、死亡率、自然增長率（2001—2008）［DB］.北京：國家青藏高原科學(xué)數(shù)據(jù)中心，2021.（Qinghai Provincial Bureau of Statistics.Birth rate，death rate and natural growth rate of population in different regions of China （2001—2008）［DB］.Beijing：National Tibetan Plateau Data Center，2021.（in Chinese））

［29］龍笛，李雪瑩，李興東，等.遙感反演2000—2020年青藏高原水儲(chǔ)量變化及其驅(qū)動(dòng)機(jī)制［J］.水科學(xué)進(jìn)展，2022，33（3）：375-389.（LONG D，LI X Y，LI X D，et al.Remote sensing retrieval of water storage changes and underlying climatic mechanisms over the Tibetan Plateau during 2000—2020［J］.Advances in Water Science，2022，33（3）：375-389.（in Chinese））

［30］孫倩.基于GRACE和GLDAS的新疆水資源時(shí)空變化研究［D］.烏魯木齊：新疆大學(xué)，2015.（SUN Q.GRACE and GLDAS data-based estimation of spatial variations in terrestrial water variations over Xinjiang［D］.Urumqi：Xinjiang University，2015.（in Chinese））

［31］STRASSBERG G，SCANLON B R，RODELL M.Comparison of seasonal terrestrial water storage variations from GRACE with groundwater-level measurements from the High Plains Aquifer（USA）［J］.Geophysical Research Letters，2007，34（14）：L14402.

［32］冉全，潘云，王一如，等.GRACE衛(wèi)星數(shù)據(jù)在海河流域地下水年開采量估算中的應(yīng)用［J］.水利水電科技進(jìn)展，2013，33（2）：42-46，67.（RAN Q，PAN Y，WANG Y R，et al.Estimation of annual groundwater exploitation in Haihe River basin by use of GRACE satellite data［J］.Advances in Science and Technology of Water Resources，2013，33（2）：42-46，67.（in Chinese））

［33］康健.基于GRACE衛(wèi)星數(shù)據(jù)反演貴州巖溶水儲(chǔ)量變化［D］.北京：中國地質(zhì)大學(xué)（北京），2016.（KANG J.Deriving Karst groundwater storage changes of Guizhou Province based on GRACE data［D］.Beijing：China University of Geosciences，2016.（in Chinese））

［34］王一平，王本明.東北平原區(qū)地下水資源的推求［J］.東北水利水電，2006，24（8）：35-36，46.（WANG Y P，WANG B M.Calculation of groundwater resources in Northeast Plain area［J］.Water Resources & Hydropower of Northeast China，2006，24（8）：35-36，46.（in Chinese））

［35］陶惠林，徐良驥，馮海寬，等.基于無人機(jī)數(shù)碼影像的冬小麥株高和生物量估算［J］.農(nóng)業(yè)工程學(xué)報(bào)，2019，35（19）：107-116.（TAO H L，XU L J，F(xiàn)ENG H K，et al.Estimation of plant height and biomass of winter wheat based on UAV digital image［J］.Transactions of the Chinese Society of Agricultural Engineering，2019，35（19）：107-116.（in Chinese））

［36］SANG Y F，WANG Z G，LIU C M.Comparison of the MK test and EMD method for trend identification in hydrological time series［J］.Journal of Hydrology，2014，510：293-298.

［37］SANG Y F，SUN F B，SINGH V P，et al.A discrete wavelet spectrum approach for identifying non-monotonic trends in hydroclimate data［J］.Hydrology and Earth System Sciences，2018，22（1）：757-766.

［38］陳飛，徐翔宇，羊艷，等.中國地下水資源演變趨勢及影響因素分析［J］.水科學(xué)進(jìn)展，2020，31（6）：811-819.（CHEN F，XU X Y，YANG Y，et al.Investigation on the evolution trends and influencing factors of groundwater resources in China［J］.Advances in Water Science，2020，31（6）：811-819.（in Chinese））

［39］郭林茂，王根緒，宋春林，等.多年凍土區(qū)下墊面條件對坡面關(guān)鍵水循環(huán)過程的影響分析［J］.水科學(xué)進(jìn)展，2022，33（3）：401-415.（GUO L M，WANG G X，SONG C L，et al.Analysis on the influence of underlying surface conditions in permafrost regions on the key water cycle processes at slope scale［J］.Advances in Water Science，2022，33（3）：401-415.（in Chinese））

Abstract：In order to improve regional groundwater development and utilization measures and to aid regional groundwater resource management，groundwater storage variability characteristics were assessed using GRACE data in three provinces of Northeast China from 2002 to 2017.GRACE and GLDAS data were combined to estimate groundwater storage variations.The estimated groundwater storage variations were verified using measured groundwater storage variations and the factors influencing groundwater storage variations were explored.The results showed that the correlation between groundwater storage variations estimated by GRACE and the measured groundwater storage variations was strong，with a correlation coefficient of 0.72.Groundwater storage variations increased significantly at the end of 2011，and the average growth rate was 2.23 mm/a.The largest surplus in groundwater storage was at the end of 2013，and the largest deficit was at the beginning of 2008.There was an evident loss in groundwater storage in the winter and autumn，and an evident surplus in the summer.Furthermore，there were differences in the spatial distribution of groundwater storage.In 2002—2012，there was less groundwater storage in the northeast than in the southwest，and in 2013—2017 there was more groundwater storage in the northeast than in southwest.There were evident groundwater storage variations in Heilongjiang Province.Furthermore，due to drought，there were significant losses in groundwater storage in Liaoning Province and Jilin Province.Groundwater storage variations were strongly correlated with changes in precipitation and had the strongest correlation with agricultural irrigation，with a correlation coefficient of 0.58.In the winter，groundwater storage variations were the most significantly correlated with snowfall，with a correlation coefficient of up to 0.426.By studying the spatiotemporal variability characteristics of groundwater storage in three provinces of Northeast China，this study aims to provide a reference for the optimal allocation of water resources and the sustainable development of the environment in China and globally.

Key words：groundwater storage;GRACE;GLDAS;precipitation;agricultural water