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

楊 飛，張成業，李 軍，宋子恒，郭添玉

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

黃河流域是我國重要的經濟地帶，構成了我國重要的生態屏障。推進黃河流域水資源節約集約利用是黃河流域生態保護和高質量發展的主要目標任務之一。彭蘇萍等[1]指出，水資源是支撐黃河流域煤炭開發、區域經濟發展和生態環境協調的基礎資源，而水資源短缺也是黃河流域中上游生態環境治理的瓶頸。因此，了解并掌握黃河流域水資源及其變化顯得尤為重要，是當前黃河流域及其礦區生態環境研究的熱點方向。傳統監測手段大多利用地面水資源動態監測網點及相應的監測系統來實現，會受到監測站點分布和數目的影響。

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

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

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

1 黃河流域概況

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

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

2 數據及分析方法

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

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

式中：Wr為某個區域的水儲量變化均值；Δh為某一個格網點的水儲量變化值；λi和φi分別為格網點經度和緯度；n為該區域總格網點數量。

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

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

3 結果分析

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

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

圖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

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

圖3 黃河流域上中下游水儲量變化月數據箱形圖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，年振幅和半年振幅都較小，說明該區域水儲量變化在近20 a 一直處于比較穩定的狀態。相較于黃河中游，黃河下游的趨勢項、年振幅及半年振幅的增幅分別達到73%、38%和98%。上述結果表明，黃河流域上中下游水儲量變化存在不同的趨勢項和周期性。確定系數的范圍為0～1，表征最小二乘法參數擬合的可靠程度的統計指標。上中下游的確定系數均大于0.5，表明時間序列分解得到較好的模型擬合參數。其中，中游確定系數為0.70，參數擬合結果優于上游和下游。

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

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

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

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

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

4 結 論

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

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

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

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