基于多源數(shù)據(jù)的科爾沁沙地陸地水及地下水儲量變化研究

趙珍珍， 馮建迪

(山東理工大學 建筑工程學院， 山東 淄博 255000)

水是生命的起源，更是人類生存和發(fā)展的基本條件之一。作為一種基礎(chǔ)資源，水被廣泛應(yīng)用于工業(yè)、農(nóng)業(yè)、生活、娛樂等人類活動。對于整個生態(tài)系統(tǒng)，水資源的地位舉足輕重，影響著全球氣候演變、全球水循環(huán)和生物多樣性等基本的生態(tài)因素[1-2]。根據(jù)不同的空間分布，水資源一般分為兩大類：地表水和地下水。地下水資源具有水質(zhì)好、蓄水穩(wěn)定、分布廣泛、易于開發(fā)等特點，是重要的淡水資源之一。在人們的生產(chǎn)生活中，地下水也擔任了重要的角色。很多地方，地下水可被直接抽取飲用，作為農(nóng)村和城市居民生活用水。在農(nóng)業(yè)生產(chǎn)中，地下水是優(yōu)質(zhì)的灌溉水源，被大量抽取使用。另外，由于其獲取途徑方便、成本低等因素，地下水也被廣泛應(yīng)用于工業(yè)和礦產(chǎn)開發(fā)生產(chǎn)中。

然而，作為一種極其寶貴的淡水資源，地下水正面臨著污染和被過量開采的現(xiàn)狀。目前，水資源匱乏已發(fā)展成為一個日益嚴重的全球性生態(tài)問題，特別是在干旱和半干旱地區(qū)尤為突出[3-5]。根據(jù)2012—2016年世界經(jīng)濟論壇年度全球風險調(diào)查，水資源危機是全球十大風險之一，排在前三位。

科爾沁沙地位于中國北方，屬于半干旱和半濕潤的過渡性氣候。年度降水量變化很大，近年來呈現(xiàn)波動下降的態(tài)勢。地表水資源日益減少，河流幾近干涸。農(nóng)業(yè)、工業(yè)、畜牧業(yè)及城市生活用水越來越依賴于地下水。因長期超量開采地下水，地下水位已嚴重下降。

水資源的實時動態(tài)監(jiān)測是實現(xiàn)水資源的保護和管理的重要參考之一。水資源的監(jiān)測手段有很多。傳統(tǒng)的方法有河流/水庫監(jiān)測站、地下水位監(jiān)測井等。傳統(tǒng)的監(jiān)測站和監(jiān)測井，易受地理環(huán)境影響，測站數(shù)量相對稀少、總體分布不均勻，監(jiān)測能力有限。另外，有些水文數(shù)據(jù)屬于內(nèi)部參考資料，出于信息安全等因素，不對外公布。以上因素很大程度上限制了對水資源的監(jiān)測、管理、保護和研究。

始于2002年的GRACE(gravity recovery and climate experiment，重力恢復(fù)和氣候探測任務(wù))衛(wèi)星重力任務(wù)為水資源的監(jiān)測研究開辟了一條嶄新的道路。GRACE提供的月重力場模型，可用來監(jiān)測地球表面淺層物質(zhì)質(zhì)量遷移。對陸地水儲量變化的監(jiān)測是GRACE的重要應(yīng)用之一。另外，結(jié)合陸地水文資料或水文模型，GRACE也可用于地下水儲量變化的監(jiān)測[6]。目前，利用GRACE重力數(shù)據(jù)監(jiān)測陸地水儲量變化(地下水儲量變化)的研究很多，研究區(qū)域大多是水量豐富的大河流域以及水資源匱乏的干旱和半干旱區(qū)。例如，密西西比河流域[7-8]、尼羅河流域[9]、亞馬遜河流域[10-12]和恒河流域[13-14]等。在中國區(qū)域，華北平原[15]、新疆天山山區(qū)[16]、海河流域[17]、黑河流域[18-19]和長江流域[19]等區(qū)域的水資源監(jiān)測一直是一個熱門的話題。

以上研究多集中在大河流域以及部分干旱和半干旱區(qū)，很少涉及科爾沁沙地水資源的研究。為了實現(xiàn)對科爾沁沙地的水資源監(jiān)測，本文擬聯(lián)合GRACE重力數(shù)據(jù)、水文模型、監(jiān)測點水井數(shù)據(jù)及給水度信息，建立科爾沁沙地水資源監(jiān)測體系，反演該地區(qū)陸地水儲量變化和地下水儲量變化，以期為當?shù)厮Y源的保護和可持續(xù)利用提供參考。

1 研究區(qū)概況

科爾沁沙地是中國最大的沙地(北緯42°30′—45°30′，東經(jīng)119°—124°)，平均海拔高度為178.5 m，呈現(xiàn)西高東低，行政區(qū)域上包括18個縣(市、旗)，總面積5.06×104km2。科爾沁沙地屬于半干旱農(nóng)牧交錯區(qū)，耕地和草地之間不斷發(fā)生互相轉(zhuǎn)換，土質(zhì)松散、黏連性差，為土地沙化提供了物質(zhì)基礎(chǔ)。同時，該地區(qū)降水量少，年平均降水量為371.3 mm，而且年際降雨量波動性大，1980年的年均降雨量只有198.7 mm，1991年年均降雨量達到了577.6 mm，經(jīng)常出現(xiàn)干旱年份，以及偶爾發(fā)生洪澇災(zāi)害；年內(nèi)降雨不均，降雨主要集中在夏季，占總降水量的69.6%，為豐水期，而春季和秋季總降水比例不足30%，冬季降水量只占百分之一左右。降雨量直接影響該地區(qū)的土壤含水量，以及地表水面積和地下水水位變化。

同時，科爾沁沙地人口數(shù)量穩(wěn)步上升，從1980年的2.37×106人，增長到2016年的3.13×106人，35 a的時間里人口增加了32.1%。伴隨著人口的快速增長，為滿足生活需求，耕地面積和牲畜存欄量也在不斷增加。科爾沁沙地經(jīng)濟水平不斷提高，GDP從1999年的1.35×1010元，經(jīng)過16年的發(fā)展，到2016年末達到了1.95×1011元，增長了14倍；農(nóng)牧民人均純收入從1999年的2 028元，增長到2014年的9 932元，增加到原來的近5倍；產(chǎn)業(yè)結(jié)構(gòu)由最初的以第一產(chǎn)業(yè)為主，逐漸轉(zhuǎn)變?yōu)榈诙a(chǎn)業(yè)占據(jù)主要地位，第三產(chǎn)業(yè)次之，第一產(chǎn)業(yè)的比重最小?？茽柷呱车靥厥獾牡乩砦恢谩夂蛱卣骱腿藶橐蛩?，共同決定了該地區(qū)地表水和地下水的不斷減少、沙地范圍大并呈現(xiàn)出波動性變化等生態(tài)環(huán)境問題。

2 反演方法

利用GRACE時變重力場反演地表面質(zhì)量變化的計算模型可以表示如下[20]：

〔ΔClmcos(mλ)+ΔSlmsin(mλ)〕

(1)

由于GRACE任務(wù)的飛行特點、共軌星間跟蹤觀測模式、軌道誤差、大氣和海洋模型誤差等因素，由原始觀測數(shù)據(jù)解算的地球重力場模型球諧系數(shù)存在條帶誤差和高頻誤差。為了削弱這些誤差的影響，本文采用去相關(guān)濾波與高斯濾波(平滑半徑為300 km)的組合濾波方法，經(jīng)過去相關(guān)濾波和高斯濾波以后，公式(1)變?yōu)椋?/p>

(2)

式中：Wl——l階高斯平滑核函數(shù)。該函數(shù)可以通過一組遞推公式表示：

W0=1

(3)

(4)

(5)

其中，

(6)

式中：r——“半寬”常數(shù)，即高斯平滑半徑；R——地球的平均半徑。

經(jīng)過去相關(guān)濾波和高斯平滑濾波數(shù)據(jù)處理，GRACE數(shù)據(jù)中的干擾因素已被剔除，剩余的信息為地球表面質(zhì)量遷移帶來的重力異常。在較短的時間尺度內(nèi)，地球表面質(zhì)量遷移一般是指陸地水儲量變化。陸地水儲量變化可進一步轉(zhuǎn)換為等效水高的形式，公式為

(7)

式中：ρ——水的密度(1 000 kg/m3)；H——等效水高。

本文采用Landerer和Swenson(2012)提出的尺度因子法，其基本原理是利用最小二乘技術(shù)，以濾波前GLDAS-NOAH得到的水儲量信號與濾波后的水儲量信號之差的最小為目標，求解得到尺度因子[21]為0.7。一個區(qū)域的陸地水儲量變化包括：土壤含水量、地表徑流、積雪和冰川以及地下水的變化，可以用公式(8)表示[22]。將土壤含水量、地表徑流、積雪和冰川的變化從GRACE計算的陸地水儲量變化中扣除，即可得到地下水儲量變化，從而實現(xiàn)利用GRACE重力數(shù)據(jù)監(jiān)測地下水儲量的目的[23]。

ΔTWS=ΔSMS+ΔSMES+ΔRESS+ΔGWS

(8)

其中：ΔTWS——陸地水儲量變化； ΔSMS——土壤含水量變化； ΔSWES——積雪和冰川變化； ΔRESS——地表徑流變化； ΔGWS——地下水儲量變化。

鑒于科爾沁沙地的實際情況和已有的資料，在計算科爾沁沙地的地下水儲量變化時，本文不考慮積雪和冰川以及地表徑流的變化，僅計算了占地表水主導(dǎo)地位的土壤含水量變化。那么，公式(8)簡化為：

ΔTWS=ΔSMS+ΔGWS

(9)

該簡化的水資源平衡公式已被應(yīng)用于多個地區(qū)的地下水儲量變化的研究中[18,24-25]。

3 結(jié)果與分析

3.1 基于GRACE數(shù)據(jù)反演陸地水儲量變化

本文采用德克薩斯大學CSR(center of space research)空間研究中心提供的200307—201012 GRACE Release-05 Level-2數(shù)據(jù)，對其進行去相關(guān)濾波與高斯濾波(平滑半徑為300 km)的組合濾波方法處理后，反演了科爾沁沙地陸地水儲量變化時間序列(圖1)。

圖1 基于GRACE數(shù)據(jù)和CPC水文模型的科爾沁沙地陸地水儲量變化(200307-201012)

由GRACE時變重力場數(shù)據(jù)反演得到科爾沁沙地200307—201012的陸地水儲量下降速率為-13.2±2.6 mm/a(下降速率即為等效水高的變化率，下同)。同時，本文采用全球水文模型CPC (climate prediction center)進行對比分析，該模型是根據(jù)全球觀測的降水分布而建立，采用的數(shù)據(jù)包括全球降水數(shù)據(jù)、溫度、潮濕度、水平風速、地表大氣壓和太陽輻射數(shù)據(jù)等。CPC提供的產(chǎn)品包括地表積雪分布及其厚度、地表以下4層的土壤含水量、土壤溫度等。作為驗證數(shù)據(jù)，CPC水文模型得到的陸地水儲量變化與GRACE的結(jié)果進行比較。從圖1可以看出：①GRACE數(shù)據(jù)和CPC水文模型反演的科爾沁沙地陸地水儲量結(jié)果整體上均呈現(xiàn)持續(xù)減小趨勢，具有較好的一致性；②二者具有相似的季節(jié)變化特性，相關(guān)系數(shù)為0.72。需要注意的是，CPC水文模型計算的水儲量變化振幅略小于GRACE的結(jié)果。這是因為GRACE的水儲量變化是所有陸地水變化的總和，相對于CPC的結(jié)果來說，還包含了其他的地表質(zhì)量變化(如地表生物總含水量變化和固體地球質(zhì)量變化)。

3.2 基于GRACE和GLDAS水文模型反演地下水儲量變化

本文構(gòu)建了兩種科爾沁沙地地下水儲量變化反演方法：一是聯(lián)合GRACE數(shù)據(jù)和水文模型反演地下水儲量變化；二是利用監(jiān)測點水井數(shù)據(jù)和給水度數(shù)據(jù)確定地下水儲量變化。最后，將這兩種方法計算的科爾沁沙地地下水儲量變化結(jié)果進行比較，相互驗證。

本文利用GLDAS水文模型計算了科爾沁沙地的土壤含水量變化(如圖2所示)。從圖2中可以看出，科爾沁沙地的土壤含水量存在明顯的季節(jié)變化特性，在夏季達到最大值，在春冬季達到最小值?？茽柷呱车匚挥谥袊狈剑瑢儆跍貛Т箨懶詺夂蛱匦?。夏季炎熱，受東南季風的影響，降水較多，占全年降水量的70%～80%；春冬季多偏北風和西北風，空氣干燥，降水稀少。

圖2 基于GLDAS水文模型的科爾沁沙地土壤含水量變化(200307-201012)

圖3記錄了2003年到2010年科爾沁沙地的降雨量月均值序列。通過對比圖2和圖3可知，科爾沁沙地的土壤含水量與降雨量是密切相關(guān)的，二者存在明顯的正相關(guān)關(guān)系。根據(jù)GRACE數(shù)據(jù)計算的陸地水儲量變化減去GLDAS水文模型計算的土壤含水量變化，得到了科爾沁沙地的地下水儲量變化(如圖4所示)。從圖4中可以看出，2003年7月到2010年12月期間，科爾沁沙地的地下水儲量總體上呈現(xiàn)下降趨勢，總體下降速率為-13.5±1.9 mm/a。其中，2003年7月到2005年7月，地下水儲量是增加的，速率為30.2±4.1 mm/a；從2005年1月開始，地下水儲量急劇下降，到2009年12月，地下水儲量下降速率為-26.7±2.7 mm/a；2009年7月到2010年12月，地下水儲量依然處于下降態(tài)勢，但下降速率明顯減緩，為-6.0±1.9 mm/a。

圖3 2003-2010年科爾沁沙地月平均降水量變化

圖4 基于GRACE數(shù)據(jù)和GLDAS水文模型的200307-201012科爾沁沙地地下水儲量變化

3.3 基于監(jiān)測點水井數(shù)據(jù)和給水度反演地下水儲量變化

3.3.1 基于監(jiān)測點水井數(shù)據(jù)反演地下水位變化 本文使用了《中國地質(zhì)環(huán)境監(jiān)測地下水位年鑒》記錄的地下水數(shù)據(jù)，研究了科爾沁沙地的地下水位變化。根據(jù)年鑒記錄，內(nèi)蒙古自治區(qū)通遼市科爾沁沙地區(qū)域共分布著6個國家級地下水監(jiān)測點。將每個監(jiān)測點的2005—2010年地下水月平均水位提取出來，計算得到6個監(jiān)測點的地下水位變化序列。再將6個監(jiān)測點的地下水位變化序列求平均值，得到了科爾沁沙地2005—2012年的月平均地下水變化序列(如圖5所示)。從圖5可以看出，科爾沁沙地地下水位一直處于下降的趨勢。

圖5 2005-2010年科爾沁沙地地下水位月平均變化量

3.3.2 基于監(jiān)測點水井數(shù)據(jù)和給水度反演地下水儲量變化 本文根據(jù)研究區(qū)域內(nèi)估算的潛水含水層的給水度結(jié)果和相關(guān)的給水度室內(nèi)試驗數(shù)據(jù)，利用Kringing模型算法，估計了科爾沁沙地的給水度數(shù)值，估值的均方根誤差優(yōu)于0.001?？紤]到區(qū)域均值的不確定性，科爾沁沙地給水度的取值區(qū)間為：(0.09～0.1)。將監(jiān)測點水井數(shù)據(jù)得到的地下水位變化乘以給水度，即可得到地下水儲量變化。圖6給出了200501—201012科爾沁沙地的地下水儲量變化，包含了給水度取0.09和0.10的計算結(jié)果。根據(jù)地下水消耗的速率，將科爾沁沙地的地下水儲量變化分為兩個階段：200501—200912，地下水下降速度較快，下降速率為-29.1 mm/a(給水度取0.1)和-26.04 mm/a(給水度取0.09)；200907—201012，地下水下降速度有所減緩，下降速率為-6.1 mm/a給水度取0.1)和-5.5 mm/a(給水度取0.09)。

3.4 科爾沁沙地地下水儲量變化結(jié)果比較分析

基于GRACE數(shù)據(jù)和GLDAS水文模型得到的地下水儲量變化，與監(jiān)測點水井數(shù)據(jù)得到的地下水儲量變化進行比較(如表1所示)。由于《中國地質(zhì)環(huán)境監(jiān)測地下水位年鑒》記錄的地下水信息始于2005年，所以對比時間段選取200501—201012。從表1可以看出，這兩種方法得到的地下水儲量變化在此期間是一致的。兩種方法相互驗證，增加了科爾沁沙地地下水儲量變化結(jié)果的可靠性。

表1 不同階段GRACE和水井數(shù)據(jù)監(jiān)測的地下水變化結(jié)果比較

綜合分析GRACE數(shù)據(jù)、GLDAS水文模型、監(jiān)測點水井數(shù)據(jù)和給水度信息得到的科爾沁沙地地下水儲量變化可知，2003—2010年科爾沁沙地的地下水儲量呈現(xiàn)消耗的趨勢，消耗速率為-13.5±1.9 mm/a，其中200501—200912月，地下水儲量消耗最快，2009年7月以后，地下水儲量有所恢復(fù)，下降速率減緩。

3.5 科爾沁沙地水儲量變化成因分析

結(jié)合氣象數(shù)據(jù)和當?shù)剞r(nóng)業(yè)模式分析科爾沁沙地水儲量變化的成因。首先，分析了科爾沁沙地降水量與陸地水儲量變化的關(guān)系。圖7給出了2003—2010年科爾沁沙地的年降水量數(shù)據(jù)，其中，年降水量的最大值出現(xiàn)在2010年，最小值出現(xiàn)在2009年。根據(jù)相關(guān)研究報道，2009年中國北方發(fā)生了嚴重的干旱[26]。同時結(jié)合圖2可以看出，由GRACE數(shù)據(jù)得到的2009年科爾沁沙地陸地水儲量下降最為明顯。隨著降水量的增加，從2010年開始陸地水儲量有所恢復(fù)。此外，2006—2009年的降水量均處于較低的水平，平均約為310 mm，是導(dǎo)致科爾沁沙地陸地水儲量多年下降態(tài)勢的主要原因之一。

然后，從科爾沁沙地的農(nóng)業(yè)生產(chǎn)模式方面分析了當?shù)氐叵滤畠α肯牡脑?。科爾沁沙地的農(nóng)作物種類繁多，包括玉米、向日葵、果樹、各種蔬菜等。其中，玉米為當?shù)氐闹饕Z食作物，被大面積廣泛種植。研究表明，當?shù)赜衩椎姆N植面積占農(nóng)作物總種植面積的70%[27-28]。玉米是一年生單子葉草本植物，秸稈很高，枝葉茂盛，生長過程中需水量很大。

圖7 2003-2010年科爾沁沙地年降水量變化趨勢

在科爾沁沙地，年平均降水量僅為350 mm，發(fā)生在玉米生長周期內(nèi)的降水占全年降水的70%左右。玉米生長需水量很大，僅依靠自然降雨是不能滿足玉米對水分需求的。本文引用張銅會等[27]的研究結(jié)果(如表2所示)，展示了科爾沁沙地玉米耗水量及用水來源[29]。玉米的生長周期分為4個階段：播種—拔節(jié)、拔節(jié)—抽穗、抽穗—灌漿、灌漿—成熟。從表2可以看出，在不同的生長階段，玉米對水的需求是不同的。此外，從該表中還可以看出，玉米總的耗水量已超過當?shù)亟涤炅?，因此，需要抽取地下水進行灌溉。灌溉水量占到了玉米總耗水量的37%左右。特別是在播種—拔節(jié)和拔節(jié)—抽穗階段，灌溉水量所占比例高達60%左右。

表2 科爾沁沙地玉米不同生長階段的需水特點[27]

結(jié)合圖7可知，2006—2009年間，科爾沁沙地的降水量嚴重不足，低于平均值。為了抵御干旱，確保農(nóng)作物的產(chǎn)量，更多的地下水被抽取到地面，用于灌溉。與此同時，地表蒸發(fā)量也隨之更加。水資源大量被消耗，地下水得不到足夠的補給，導(dǎo)致地下水位不斷下降。

圖8是利用2005—2010年《中國地質(zhì)環(huán)境監(jiān)測地下水位年鑒》地下水位數(shù)據(jù)計算的科爾沁沙地地下水位季節(jié)變化。從圖中可以看出，每年的地下水位從3月份開始下降，一直到8月份下降至最低點。8月份以后，地下水得以補給，水位上升。根據(jù)前文描述的當?shù)剞r(nóng)業(yè)模式和降水特點，地下水消耗和恢復(fù)的過程可以理解為：春季降雨稀少，當?shù)鼐用翊罅砍槿〉叵滤詫勾汉?，?dǎo)致地下水從3月份開始下降；夏季農(nóng)作物生長的過程中，水的需求量是巨大的，降雨量不能滿足農(nóng)作物的生長需求，因此大量地下水被不間斷的抽取到地面，導(dǎo)致地下水位持續(xù)性下降；秋季作物開始成熟，對水的需求也在逐漸減少，地下水得以恢復(fù)；秋季以后直到冬季結(jié)束，該地區(qū)無農(nóng)作物生長，地下水在冬末春初時，恢復(fù)到最大值。另外，從圖8中還可以發(fā)現(xiàn)，雖然地下水存在消耗和恢復(fù)的周期性，但是每年的地下水位都不能恢復(fù)到最初的水位，地下水位在不斷的下降。

圖8 基于監(jiān)測點水井數(shù)據(jù)的2005-2010年科爾沁沙地地下水位季節(jié)變化

4 討論與結(jié)論

本文采用200307—201012的GRACE時變重力場模型Release-05 Level-2，采用去相關(guān)濾波與高斯平滑濾波(平滑半徑300 km)結(jié)合的濾波方法，反演了科爾沁沙地陸地水儲量變化，將其與CPC水文模型反演結(jié)果進行了對比分析；結(jié)合GLDAS水文模型計算了土壤含水量變化，并分析了科爾沁沙地地下水處理變化的時空特征；利用監(jiān)測點水井數(shù)據(jù)和給水度信息確定地下水儲量變化，并與GRACE反演結(jié)果進行了初步驗證。研究結(jié)果表明：①科爾沁沙地陸地水儲量整體呈現(xiàn)減少趨勢，減少速率為-13.2±2.6 mm/a，與CPC水文模型反演結(jié)果具有較好的一致性；②科爾沁沙地地下水儲量整體上呈現(xiàn)下降趨勢，下降速率為-13.5±1.9 mm/a，與監(jiān)測點水井數(shù)據(jù)反演結(jié)果具有較好的一致性。

結(jié)合氣象數(shù)據(jù)和當?shù)剞r(nóng)業(yè)模式分析了科爾沁沙地水儲量變化的成因：首先，降水量變化與陸地水儲量變化有直接關(guān)系，2003-2010年期間，年降水量最小值出現(xiàn)在2009年，此時對應(yīng)的科爾沁沙地的陸地水儲量下降最為明顯，年降水量最大值出現(xiàn)在2010年，此時期的陸地水儲量也開始有所恢復(fù)；第二，根據(jù)科爾沁沙地玉米不同生長階段需水量以及對應(yīng)階段降水量的特點，發(fā)現(xiàn)干旱和農(nóng)業(yè)灌溉是導(dǎo)致地表水減少和地下水位降低的重要原因之一。

基于GRACE衛(wèi)星時變重力場數(shù)據(jù)反演長時間序列、大尺度的陸地水儲量變化和地下水儲量變化，該方法有效地彌補了傳統(tǒng)水文監(jiān)測中因觀測站點分布不均勻、觀測范圍有限等因素導(dǎo)致資料獲取不充分、數(shù)據(jù)不均勻的缺陷。隨著后續(xù)衛(wèi)星重力探測計劃的相繼實施、時變重力場反演方法的改進和水文模型的日趨完善，將進一步提高衛(wèi)星重力觀測結(jié)果中地球物理信號的時空分辨率，以及信號分離和提取的精度，為研究固體地球物質(zhì)遷移以及精化和驗證全球和局部水文模型提供可靠的地球物理信息。