基于GRACE 的黃河流域陸地水儲量時空變化研究

楊 飛，張成業(yè)，李 軍，宋子恒，郭添玉

(1.中國礦業(yè)大學(xué)(北京) 煤炭資源與安全開采國家重點實驗室，北京 100083；2.中國礦業(yè)大學(xué)(北京) 地球科學(xué)與測繪工程學(xué)院，北京 100083；3.煤炭開采水資源保護與利用國家重點實驗室，北京 102211)

黃河流域是我國重要的經(jīng)濟地帶，構(gòu)成了我國重要的生態(tài)屏障。推進黃河流域水資源節(jié)約集約利用是黃河流域生態(tài)保護和高質(zhì)量發(fā)展的主要目標(biāo)任務(wù)之一。彭蘇萍等[1]指出，水資源是支撐黃河流域煤炭開發(fā)、區(qū)域經(jīng)濟發(fā)展和生態(tài)環(huán)境協(xié)調(diào)的基礎(chǔ)資源，而水資源短缺也是黃河流域中上游生態(tài)環(huán)境治理的瓶頸。因此，了解并掌握黃河流域水資源及其變化顯得尤為重要，是當(dāng)前黃河流域及其礦區(qū)生態(tài)環(huán)境研究的熱點方向。傳統(tǒng)監(jiān)測手段大多利用地面水資源動態(tài)監(jiān)測網(wǎng)點及相應(yīng)的監(jiān)測系統(tǒng)來實現(xiàn)，會受到監(jiān)測站點分布和數(shù)目的影響。

GRACE(Gravity Recovery and Climate Experiment)重力衛(wèi)星計劃由美國國家宇航局和德國空間飛行中心聯(lián)合實施，衛(wèi)星于2002 年發(fā)射并于2004 年公開其重力數(shù)據(jù)[2]。GRACE 衛(wèi)星采用極地近圓軌道低衛(wèi)星跟蹤技術(shù)，通過獲取星間距離的變化及其殘差來恢復(fù)時變重力場信息[3-4]。在消除天體引潮力對海洋、大氣和相關(guān)地球動力過程引起的時變重力場影響后，GRACE 時變重力場反映的主要是兩極冰蓋、山岳冰川、陸地水儲量以及海平面變化信息[5-6]。因此，利用GRACE 時變重力場可以反演陸地水儲量變化信息，成為當(dāng)前進行中長尺度陸地水資源時空變化分析的重要手段[7-9]。J.Wahr 等[10-11]推導(dǎo)了GRACE 數(shù)據(jù)計算陸地水儲量變化的基礎(chǔ)理論公式，首次利用GRACE衛(wèi)星11 個月的數(shù)據(jù)反演密西西比河流域、亞馬孫河流域和孟加拉灣地區(qū)的儲水量變化，證明GRACE 具有在月尺度上探測到這些流域儲水量變化的潛力。隨后，眾多學(xué)者提出多種基于GRACE 的陸地水儲量反演的改進方法[12-17]，并成功應(yīng)用于世界范圍內(nèi)不同區(qū)域陸地水儲量的變化研究中[18-21]。

當(dāng)前，利用GRACE 衛(wèi)星數(shù)據(jù)計算陸地水儲量變化的方法主要包括球諧系數(shù)方法和Mascons 方法[22-23]。鄒賢才等[24]詳細(xì)對比了上述2 種方法并發(fā)現(xiàn)Mascons 方法存在巨大優(yōu)勢；其他學(xué)者的研究結(jié)果同樣證明Mascons 方法在一些流域的表現(xiàn)要優(yōu)于球諧系數(shù)方法[25-26]。2015 年以來，美國得克薩斯州州立大學(xué)空間研究中心(CSR)和美國噴氣推進實驗室(JPL)發(fā)布用Mascons 方法計算陸地水儲量的產(chǎn)品，分別稱為CSRM 和JPL-M。2 種產(chǎn)品在不同區(qū)域陸地水儲量變化的比較和研究中表明，CSR-M 產(chǎn)品可以更好地反映整個區(qū)域地表陸地水遷移的實際情況[27-28]。

基于此，利用CSR 提供的RL06 Mascons 水儲量變化產(chǎn)品，開展黃河流域上、中、下游近20 a 陸地水儲量變化的研究。通過緯圈長度加權(quán)平均計算每個子流域的水儲量變化均值，采用時間序列分解方法得到黃河流域水儲量變化的趨勢、年周期等整體規(guī)律，利用空間分析比較不同區(qū)域水儲量變化的特征，以期為更好地了解和掌握黃河流域水資源及其變化規(guī)律提供依據(jù)，為流域礦區(qū)的生態(tài)保護與可持續(xù)發(fā)展提供參考。

1 黃河流域概況

黃河發(fā)源于青海省巴顏喀拉山脈，橫跨青藏高原、內(nèi)蒙古高原、黃土高原和黃淮海平原，先后流經(jīng)青海、四川、甘肅、寧夏、內(nèi)蒙古、陜西、山西、河南和山東9 個省/自治區(qū)。其東西、南北跨度分別長約1 900 km和1 100 km，經(jīng)緯度介于96°E-119°E 和32°N-42°N，流域總面積79.5 萬km2。受大氣環(huán)流和季風(fēng)環(huán)流的影響作用比較復(fù)雜，黃河流域內(nèi)不同地區(qū)的氣候差異顯著，主要屬于南溫帶、中溫帶和高原氣候區(qū)。根據(jù)黃河水利委員會的劃分方案[29]，黃河上游指的是從源頭到內(nèi)蒙古河口鎮(zhèn)，河長3 472 km，流域面積38.5 萬km2；黃河中游指的是從內(nèi)蒙古河口鎮(zhèn)到河南鄭州桃花峪，河長1 206 km，流域面積34.3 萬km2；黃河下游指的是從鄭州桃花峪到渤海，河長786 km，流域面積2.3 萬km2。黃河流域中上游分布有晉北、晉東、晉中、黃隴、陜北、神東和寧東7 個國家大型煤炭基地，是我國目前主要煤炭生產(chǎn)區(qū)和調(diào)出區(qū)；而黃河下流分布有河南、魯西等煤炭基地。黃河流域及其分界線如圖1 所示。

圖1 黃流流域分布Fig.1 Map of the Yellow River Basin

2 數(shù)據(jù)及分析方法

本文采用美國德克薩斯大學(xué)空間研究中心CSR提供的2002 年4 月到2017 年6 月的最新版本RL06 Mascons 產(chǎn)品。較RL05 版本，RL06 使用了改進的背景模型、新的地球系統(tǒng)參考框架和新的數(shù)據(jù)處理策略，顯著提升了數(shù)據(jù)質(zhì)量。CSR RL06 Mascons 產(chǎn)品是對GRACE 衛(wèi)星的高精度星間觀測數(shù)據(jù)Level1B 進行直接解算，C20 項采用衛(wèi)星激光測月(SLR)數(shù)據(jù)的C20項替換，同時考慮地心一階項和冰后回彈改正，采用Tikhonov 正則化進行方程組的求解。最終，CSR RL06 Mascons 數(shù)據(jù)采用NetCDF 格式，以0.25°×0.25°分辨率格網(wǎng)的形式發(fā)布。

為分析2002-2017 年黃河流域上中下游各個區(qū)域的水儲量變化，需要計算每個區(qū)域的平均水儲量變化值。采用對每個格網(wǎng)點的水儲量變化數(shù)據(jù)直接取平均值的傳統(tǒng)方法，相當(dāng)于認(rèn)為每個格網(wǎng)點等權(quán)。然而，CSR RL06 Mascons 產(chǎn)品為0.25°×0.25°格網(wǎng)點，數(shù)據(jù)在地球不同緯度之間代表的實際面積差異較大，需要考慮地球曲面的影響。本文采用緯圈長度加權(quán)方法，根據(jù)地理緯度為每個格網(wǎng)點賦予權(quán)重來獲取不同區(qū)域的水儲量變化均值。計算公式如下：

式中：Wr為某個區(qū)域的水儲量變化均值；Δh為某一個格網(wǎng)點的水儲量變化值；λi和φi分別為格網(wǎng)點經(jīng)度和緯度；n為該區(qū)域總格網(wǎng)點數(shù)量。

本文采用時間序列分解方法來分析黃河流域上中下游不同區(qū)域水儲量變化的時間序列特征。將各個區(qū)域水儲量變化均值時間序列W(t)分解為趨勢項、年周期項和半年周期項。利用最小二乘法進行參數(shù)擬合，擬合形式如下：

式中：t為時間；a為常數(shù)項；b為趨勢項；c1、c2和d1、d2分別為年周期變化和半年周期變化。

3 結(jié)果分析

通過緯圈長度加權(quán)方法得到每個月黃河流域上中下游的水儲量變化均值，采用等效水高表示[3-4,6-7]，其時間序列及線性擬合趨勢如圖2 左側(cè)所示，可以看到，2002-2017 年，黃河上游水儲量變化均值的波動范圍較小，最大和最小值分別為32.8、-56.9 mm，線性擬合趨勢項為-2.3 mm/a，水儲量變化遞減的趨勢不明顯。而黃河流域中游和下游的水儲量變化均值表現(xiàn)出明顯的遞減趨勢，其趨勢項分別達到-10.3、-17.5 mm/a。在黃河流域中游，水儲量變化均值的最大和最小值分別為108.9、-168.6 mm；在黃河流域下游，它們分別達到了143.8、-301.9 mm。上述結(jié)果表明近20 a 黃河流域3 個流域段水儲量變化的差異性，上游水儲量變化穩(wěn)定，月均值差異較小，趨勢性弱；中游的水儲量變化呈遞減趨勢；下游的水儲量變化范圍最大且遞減趨勢尤為明顯。

去趨勢后，黃河流域上中下游水儲量變化均值的時間序列如圖2 右側(cè)所示，3 個流段都表現(xiàn)出周年性及季節(jié)性變化。值得注意的是，黃河中下游水儲量變化在2003 年秋季有顯著的提升。這是因為在2003 年上半年，黃河流域來水量較常年減少40%～70%，為實測記錄以來最小值；進入汛期以后，黃河流域大部分地區(qū)降水量增加，較常年增多10%～70%，洪水時水庫得以蓄水；所以在2003 年秋季，黃河流域水儲量變化迅速從負(fù)值變化到正值。

圖2 黃河流域上中下游水儲量變化均值的時間序列及趨勢和去趨勢后的時間序列Fig.2 Time series and trend of the mean water storage anomalies(left),and the time series after detrending(right) in the upper,middle and lower Yellow River Basin

為進一步分析黃河流域不同流段水儲量變化的月度規(guī)律以及數(shù)據(jù)的變化幅度，本文統(tǒng)計了2002 年4 月到2017 年6 月3 段流域各個月份的水儲量變化均值并繪制箱形圖，如圖3 所示。可以發(fā)現(xiàn)，黃河上中下游的水儲量變化月度范圍有顯著的區(qū)域差異，上游整體處于-40～40 mm，而中游和下游分別處于-100～100 mm和-200～200 mm。黃河上中下游每個月的數(shù)據(jù)都在一個相對穩(wěn)定的范圍內(nèi)波動，春冬季節(jié)的波動明顯大于夏秋季節(jié)，這與該流域“冬干春旱，夏秋多雨”的氣候特點密切相關(guān)。另外，有少數(shù)波動表現(xiàn)為異常值，這些異常值集中于夏季，受到夏季降水和高溫蒸發(fā)等因素影響。3 個流域段多年平均水儲量變化在9 月份達到最大值，且在整個夏季和秋季維持較大值。

圖3 黃河流域上中下游水儲量變化月數(shù)據(jù)箱形圖Fig.3 The boxplots of monthly mean data of the terrestrial water storage anomalies in the upper,middle and lower Yellow River Basin

通過對黃河流域上中下游水儲量變化均值的時間序列進行分解，得到其趨勢、年周期和半年周期等信息，見表1。可以看到，時間序列分解得到的趨勢項與線性擬合的趨勢差別很小。黃河上游的趨勢項只有-2.17 mm/a，年振幅和半年振幅都較小，說明該區(qū)域水儲量變化在近20 a 一直處于比較穩(wěn)定的狀態(tài)。相較于黃河中游，黃河下游的趨勢項、年振幅及半年振幅的增幅分別達到73%、38%和98%。上述結(jié)果表明，黃河流域上中下游水儲量變化存在不同的趨勢項和周期性。確定系數(shù)的范圍為0～1，表征最小二乘法參數(shù)擬合的可靠程度的統(tǒng)計指標(biāo)。上中下游的確定系數(shù)均大于0.5，表明時間序列分解得到較好的模型擬合參數(shù)。其中，中游確定系數(shù)為0.70，參數(shù)擬合結(jié)果優(yōu)于上游和下游。

表1 黃河流域上中下游水儲量變化均值的時間序列分解統(tǒng)計Table 1 Statistics of the water storage anomalies in the upper,middle and lower Yellow River Basin through time series decomposition

為表現(xiàn)黃河流域水儲量變化的空間分布情況，本文對該區(qū)域0.25°×0.25°格網(wǎng)點數(shù)據(jù)進行時間序列分解，獲取每個格網(wǎng)點2002 年6 月到2017 年6 月的時間序列擬合參數(shù)。基于此，開展黃河流域水儲量變化空間分析。圖4 表示水儲量變化趨勢(圖4a)和確定系數(shù)(圖4b)的空間分布。可以看到，黃河流域水儲量變化為增加趨勢的地區(qū)主要集中在黃河上游源頭附近區(qū)域，其年增加量在5 mm 左右。黃河下游及中游靠近黃河的區(qū)域表現(xiàn)出明顯的遞減趨勢，年減少量基本達到15 mm 以上。其他一些區(qū)域趨勢項大多處于-10 mm/a 以下。黃河流域水儲量變化存在隨經(jīng)度變化的規(guī)律，由西向東表現(xiàn)出越來越明顯的遞減趨勢。從確定系數(shù)的空間分布可以看到黃河流域擬合結(jié)果的南北差異，流域北側(cè)擬合可靠度整體優(yōu)于流域南側(cè)。整個流域確定系數(shù)大于0.5 的格網(wǎng)點占64%，大于0.3的格網(wǎng)點達到85%。

圖4 黃河流域水儲量參數(shù)的空間分布Fig.4 Spatial distribution of water storage factors in the Yellow River Basin

圖5 表示黃河流域水儲量變化的年周期振幅和半年周期振幅的空間分布。可以看到，整個流域年周期振幅都明顯大于其半年周期振幅。黃河上游源頭附近區(qū)域、中游桃花峪及黃河下游區(qū)域年周期振幅較大，達到30 mm 及以上。在上述區(qū)域，高山融雪、降水量等季節(jié)性差異是導(dǎo)致其年周期振幅相對較大的原因。黃河流域中上游其他區(qū)域的年周期振幅較小，處于10～20 mm。黃河流域水儲量變化的半年周期振幅表現(xiàn)出由西北向東南增大的趨勢，在黃河上游小于10 mm，黃河中游處于10～20 mm，黃河下游達到了30 mm。這與不同區(qū)域所處的氣候環(huán)境等因素有直接關(guān)系。

圖5 黃河流域水儲量不同變化周期的空間分布Fig.5 Spatial distribution of water storage anomalies in different period in the Yellow River Basin

4 結(jié) 論

a.黃河上游水儲量變化的趨勢性不明顯，只有-2.3 mm/a；中下游表現(xiàn)出顯著的遞減趨勢，分別達到-10.3 mm/a 和-17.5 mm/a。GRACE 水儲量變化數(shù)據(jù)很好地反映出2003 年黃河流域來水量、降水量的特殊性。

b.GRACE 水儲量變化數(shù)據(jù)表明黃河流域水儲量變化在9 月達到最大值，而春冬季節(jié)波動較大，與該流域“冬干春旱，夏秋多雨”的氣候特點密切相關(guān)。水儲量變化的異常值集中于夏季，與夏季降水和高溫蒸發(fā)等因素相關(guān)。

c.黃河上游源頭附近區(qū)域的水儲量變化呈增長趨勢，增長率為5 mm/a 左右。黃河流域水儲量變化存在隨經(jīng)度變化的趨勢，由西向東表現(xiàn)出越來越明顯的遞減趨勢，最大遞減率超過20 mm/a。

d.黃河流域年周期振幅明顯大于其半年周期振幅。受高山融雪、降水量等季節(jié)性差異的影響，年周期振幅在上游源頭附近區(qū)域、中游桃花峪及下游區(qū)域達到較大值。半年周期振幅表現(xiàn)出由西北向東南增大的趨勢，可能與區(qū)域的氣候環(huán)境等因素相關(guān)。