基于GRACE數據的長江流域地表水儲量變化

葉勤玉，張 強，李小龍，金濤勇，何澤能，楊世琦

(1.重慶市氣象科學研究所，重慶 401147；2.重慶市氣象局，重慶 401147；3.武漢大學測繪學院，湖北 武漢 430072)

0 引言

陸地水儲量是水文循環的關鍵一環，是評估區域水量收支狀態的重要參數。目前對陸地水儲量變化的監測方法主要有：站點觀測法[1]、衛星遙感反演法[2-3]、物理模型模擬法[4]、GRACE重力衛星反演法[5]。站點觀測法受實際地形地貌的影響，很難對大尺度區域進行測量；衛星遙感反演方法只能反演表層土壤水分，不能對深層水進行探測；而物理模型模擬方法中觀測站點的分布不均影響了其精度和應用范圍。GRACE重力衛星具有全球觀測尺度統一、分布均勻等優勢，被廣泛的應用于全球以及區域尺度的陸地水儲量變化監測評估中[6]。

近年來，利用GRACE重力衛星數據研究陸地水儲量變化特征備受研究者的關注。馮貴平等[7]利用GRACE重力衛星數據反演了2002年4月—2011年2月的月間隔全球陸地水總儲存量，并分析了其季節性和長期變化特征。結果表明，GRACE反演的全球地下水儲量具有明顯的季節性變化和長期趨勢。尼勝楠等[8]利用GRACE數據反演長江、黃河流域的水儲量變化，并與水文模型進行對比，研究其季節性和年際變化特征。結果表明，GRACE反演水儲量季節變化振幅略大于水文模型，可揭示厘米級等效水高變化。劉任莉等[9]利用GRACE衛星數據反演了2003—2010年中國陸地水儲量變化，結果表明，GRACE反演的陸地水儲量變化結果與基于實測水文、氣象數據的結果符合較好，對2009—2010年西南五省特大旱災也有著較好的反映。徐永明等[10]利用GRACE RL06數據及GLDAS數據對云貴高原陸地、地表以及地下水儲量變化進行了反演，并對水儲量變化時序進行了特征分析。結果表明，2002—2010年，GRACE與GLDAS反演的陸地水儲量的相關系數達0.94，且兩者的變化趨勢一致。魏光輝等[11]對GRACE估算的陸地水儲量進行降尺度，分析了2002—2015年間塔里木河流域局部的干旱特征。結果表明，流域內GRACE陸地水儲量與GLDAS土壤水的相關系數達0.35。陳智偉等[12]反演了2005—2012年珠江流域的陸地水儲量變化，結果表明，ITSG-GRACE 2018模型反演的珠江流域陸地水儲量季節性變化特征與GLDAS水文模型、降水數據以及地下水測井監測數據具有較好的一致性。丁一航等[13]利用GRACE研究發現我國四大流域水儲量變化信號具有多種周期特性，研究結果與GLDAS結果具有很強的相關性。曹艷萍等[14]利用GRACE重力衛星反演河南省水儲量變化，發現河南省陸地水儲量變化呈現典型的“余弦曲線”分布特征。

長江流域是我國面積最大、水量最豐富的流域，同時也是我國人口密度最大的區域之一[15]。研究長江流域的水儲量變化對探究流域水儲量變化特征、分析不同水文地質條件對水儲量變化影響、指導自然資源的合理開發和生產活動順利進行、維護和保障國家安全有著重要的戰略意義。本研究利用2002年4月—2017年6月的GRACE數據，反演了長江流域的水儲量變化，以期為長江流域水資源合理利用與配置提供參考。

1 研究區與數據

1.1 研究區域

長江是我國最長的河流，發源于青藏高原，自西向東流經青海、西藏、四川、云南、重慶、湖北、湖南、江西、安徽、江蘇、上海等11個省市，并最終注入東海。長江流域，是指長江干流和支流流經的區域，是世界第三大流域，流域總面積約180萬 km2，維護和保障了我國的生態安全、供水安全、糧食安全、能源安全[16]。受西太平洋副熱帶高壓的強度以及位置的影響[17]，長江流域年降雨量與暴雨時空分布不均，年降雨量約1 100 mm，水儲量變化具有很強的周期性[18]。

圖1 研究區高程圖Fig.1 DEM of Yangtze River Basin

1.2 數據

GRACE(Gravity Recovery and Climate Experiment)衛星由美國國家航空航天局(NASA)和德國航空航天中心(DLR)聯合研制，于2002年3月發射成功，旨在提供高精度、高分辨率的地球重力場靜態分布數據，監測月尺度上地球重力場的時空變化特征。GRACE數據由德州大學空間研究中心(CSR)、德國地學中心(GFZ)和噴氣推進實驗室(JPL)3家機構處理并發布，數據形式為月重力場球諧系數，經過相應的反演程序即可獲取區域質量變化信息。除球諧系數外，3家機構還同時發布了Mascon產品，該產品為進行適當空間平滑后的全球格網點重力異常，可以用來進行后續驗證和進一步對比分析。

本文選擇的GRACE數據為CSR于2018年6月發布的球諧系數產品，版本為最新的RL06，數據類型為GSM。RL06版本可以提供最大階數為96階和60階的兩組月重力場模型，而GSM表示扣除了高頻海洋信號和非潮汐因素等的影響，產品時間跨度為2002年4月—2017年6月。

為了對GRACE反演結果進行驗證，我們選擇了3種氣象數據模型：降雨模型選擇TRMM 3B43 v7數據集。該數據空間分辨率0.25°×0.25°，時間分辨率為1月，觀測范圍為50°S～50°N，180°W～180°E，覆蓋包括長江流域在內的我國大部分區域，時間跨度為2002年1月—2018年12月；徑流數據選擇水文站實測數據，包括大通站、宜昌站等，時間跨度為2000年1月—2020年10月；蒸散模型選擇ESA于2019年5月發布的GLEAM v3.3a全球蒸散數據集，該模型的時間跨度為1980—2018年共39 a，模型時間分辨率為月，空間分辨率為0.5°×0.5°。所有要素的時間間隔統一為2002年4月—2017年6月。表1為數據信息匯總表。

表1 數據信息匯總表Tab.1 Data information table

2 地表水儲量反演與驗證方法

地球重力場通常可以用大地水準面的形狀來描述，地球上某一點的大地水準面N可以表示為一組球諧系數的總和[19]：

(1)

(2)

大地水準面的變化會導致地表密度的重新分布，記密度變化為△ρ(r，θ，λ)，球諧系數的變化可表示為：

(3)

其中，ρavs為地球平均密度，數值為5 517 kg·m-3。不妨假設△ρ為一厚度為H的薄層，并設地球表面面密度變化為△σ，則可以將其表示為△ρ的積分：

(4)

(5)

公式(5)反映了地表質量變化對大地水準面的影響，而在地表物質遷移的同時，作為滯彈性體的地球也會在表面負荷變化的作用下產生形變，并引起大地水準面的變化，即[20]：

(6)

其中kl為l階負荷勒夫數，表征滯彈性體對表面負荷的形變響應。綜合公式(5)、(6)，大地水準面的總變化為：

(7)

將面密度變化△σ展開：

(8)

其中ρw表示水的密度，數值為1 000 kg·m-3，△lm、△?lm為無量綱規格化球諧系數。根據公式(8)，有：

(9)

聯立公式(5)、(6)、(7)、(9)可得：

(10)

聯立公式(2)、(10)得：

(△Clmcos(mλ)+△Slmsin(mλ))

(11)

在陸地水儲量變化分析中，通常使用等效水高EWH(Equivalent Water Height)來表示地表質量的變化，△EWH和△σ有如下關系：

(12)

聯立公式(11)、(12)，有：

(13)

其中ρave為地球平均密度，其數值為5 517 kg·m-3，ρw表示水的密度，數值為1000 kg·m-3。

GRACE發布的數據反映了地球靜態重力場，動態變化信息會被大量靜態信號淹沒，因此需要將背景場去除掉，以放大時變信號，通常選擇2004—2009年的平均值作為背景場，并從每月球諧系數中扣除[21]。注意到GRACE衛星軌道中心為地球的瞬心，因此與地球質心有關的C10、C11、S113項無法反映真實的重力場信號；同時GRACE衛星軌道近圓，對反映地球扁率的C20項也不敏感。因此，在利用GRACE數據反演地表質量變化前，需要對以上數據進行替換，通常選擇基于海洋模型的一階項系數和基于衛星激光測距SLR(Satellite Laser Ranging)的C20項系數作為替代[22-23]。

同時，在反演過程中，需要將系數進行截斷以降低噪聲，但這種方法并不能將噪聲清除干凈，這是因為此時的球諧系數除仍存在著高頻噪聲和南北方向上的條帶誤差。這些誤差有兩種來源：一是GRACE重力場模型高階位系數中存在著較大的隨機誤差，需要用高斯濾波進行去除，通常選擇300 km高斯濾波[24]；二是相同次數、不同階數的系數之間存在著相關關系，在空間域上表現為南北方向出現的明顯條帶，需要通過于PxMy模型方法進行降噪，通常選擇P3M9法[25]。上述兩種濾波方法雖然能較大程度上去除噪聲，但同時也帶來了新的誤差：泄露誤差，這是一種伴隨著信號截斷和空間平滑出現的誤差，需要使用基于尺度因子法的泄露誤差改正方法進行改正[26]。

冰后回彈效應PGR(Post Glacial Rebound)，又稱冰川均衡調整GIA(Glacial Isostatic Adjustment)，揭示了始于7萬年前、終于1.15萬年前的末次冰期對黏彈性地球的影響。在末次冰期，地球上絕大部分區域都被大冰原覆蓋，地表在冰原巨大的重力作用下不斷凹陷；1.8萬年前冰川開始消退，大片冰原消融，但由于地幔的黏彈性，原來負載冰蓋的地區會在壓力減小后向上彈起，對地殼運動、重力場反演、地表質量遷移等都有著重要的影響。本文使用基于最新的GIA模型——A13模型的冰川均衡調整方法[27]。

在獲取指定區域地表水儲量變化TWSC(Terrestial Water Storage Change)時間序列后，需要對序列進行頻譜分析以得到該區域水儲量變化的頻譜域特征，包括趨勢項、季節項、異常項等信號。因此本文選用基于最小二乘的水儲量變化時間序列分析方法[28]，可由以下函數模型近似表達：

(14)

式中，a0為擬合殘差，a1為趨勢值，fi(i=1,2)為信號頻率，單位為1/a，φi為信號相位。當i=1時對應周年信號，此時fi=1，bi為周年信號振幅；當i=2時對應半周年信號，此時fi=2，bi、ci為半周年信號振幅。

基于GRACE衛星的水儲量變化反演及分析流程如圖2所示：①系數替換。將2004—2009年的系數平均值作為背景場，并從每月球諧系數中扣除；②使用海洋模型和SLR數據對一階項和C20系數進行替換；③對重力場模型進行300 km高斯濾波和P3M9去相關濾波，得到基本消除高階誤差和南北條帶的質量變化空間分布；④利用尺度因子法對泄露誤差進行改正，恢復被信號截斷和高斯濾波污染的信號；⑤利用A13模型消除冰后回彈效應的影響；⑥利用最小二乘法對獲得的地表水儲量變化時間序列進行分析和處理，從而得到序列的趨勢、周年信號、半周年信號等重要參數，并分析它們的時空變化特征。

為驗證水儲量變化結果的準確性，采用水文要素模型對GRACE結果進行評估。基于流域水量平衡方程，區域水儲量變化與降雨(P)、徑流(R)、蒸散(ET)之間存在著如下關系：

△TWSC=P-ET-R

(15)

因此，在獲取氣象要素時間序列和空間分布后，可以擬合水儲量變化的時空規律，擬合效果越好，就說明GRACE結果越精確。同時，由于長江流域位于東亞季風帶上，降雨變化和后續的徑流蒸散過程是導致其水儲量變化的主要影響因素，因此對氣象要素的研究有助于了解水儲量變化的主要驅動力和驅動過程，有助于評估流域尺度的水量收支狀態。

3 結果與分析

3.1 長江流域水儲量時間變化特征分析

時間變化特征即為時間序列的趨勢，周年、半周年項的振幅以及相位等信息。使用基于最小二乘的水儲量變化時間序列分析方法，對時間序列進行頻譜分析，獲取序列趨勢項、季節項等變化特征。

圖2 地表水儲量變化反演與分析流程Fig.2 Flow chart of terrestrial water storage changes inversion and analysis

長江流域2002年4月—2017年6月的地表水儲量變化時間序列及其最大影響源：降雨的時間序列對比見圖3，特征分析結果如表2所示[29]。

圖3 長江流域2002年4月—2017年6月地表水儲量變化與降雨時間序列對比Fig.3 The time series terrestrial water storage changes in the Yangtze River Basin from April 2002 to June 2017

結合表2，從趨勢項可以推斷2002年4月—2017年6月這16 a間，長江流域整體水儲量處在一個緩慢上升的狀態，長江流域水儲量變化速率為4.1±1.1 mm·ɑ-1，總體水量在此期間增加了大約1.1×1011T，說明長江流域在近16 a來處于豐水狀態，平均每年增加約1.6個太湖的水量。

表2 長江流域2002年4月—2017年6月地表水儲量變化與降雨時間序列特征分析表Tab.2 The characteristics of the time series terrestrial water storage changes and precipitation in the angtze River Basin from April 2002 to June 2017

考慮到長江流域復雜的地質構造、特殊地形以及人為因素的影響，長江流域水儲量變化的周期信號有很多，最主要的是周年信號和半周年信號。表2的結果表明，長江流域地表水儲量變化中有著56.0 mm等效水柱高的周年性波動，半周年波動的振幅為10.9 mm等效水柱高；而降雨的周年振幅是地表水儲量周年振幅的約1.5倍，半周年振幅則與地表水儲量相當，相關系數高達0.75，反映了降雨對長江流域地表水儲量變化的深刻影響。地表水儲量變化的周年相位約-6.6°，且振幅為負，揭示了長江流域水儲量在每年8月左右達到最大值，在2月左右達到最小值，這與實際情況相吻合。長江流域受亞熱帶季風氣候的影響，夏汛集中在5—8月，降雨的增加會導致水儲量的大幅度增加，于是在8月左右出現峰值；而汛期過后，水儲量逐漸平穩，隨著秋冬季節的到來，降水減少，加之大量蒸散發作用，水儲量逐漸減少，并在2月左右出現谷值，這一變化規律同樣體現在降雨序列中。降雨周年項的相位為-0.1°，表明降雨在每年的7月達到峰值，并在1月達到谷值，這一結果與水儲量之間存在著約1個月的間隔，這是因為降雨是一個動態過程，且降雨除了蓄積在流域內，還有部分以徑流、水蒸氣的形式流出，所以降雨過程結束后需要一段時間才能在水儲量中有所體現。

長江流域地表水儲量半周年相位為-6.6°，振幅為10.9 mm，說明長江流域水儲量除了“2月谷值—8月峰值—2月谷值”這種水儲量周期變化模式外，還存在著“2月峰值—5月谷值—8月峰值—11月谷值—2月峰值”的變化規律，且規模只有前者的20%左右；而降雨的半周年相位為0.2°，揭示了降雨“12月峰值—3月谷值—6月峰值—9月谷值—12月峰值”的變化規律，與水儲量變化之間存在2個月的時延，與周年項的1個月不同，推測可能是因為振幅較小，所以需要累積更長時間才能被GRACE衛星捕捉到，而降雨的半周年振幅為10.1 mm，幾乎與水儲量變化振幅相同。

3.2 長江流域水儲量空間變化特征分析

空間分布特征體現了區域內不同格網點、不同參數的分布特征和規律。圖4a、4b分別為長江流域2002年4月—2017年6月1°×1°格網點上地表水儲量變化趨勢項、周年項振幅絕對值的空間分布，可以分辨出長江流域上各格網點、各區域、各水系的水儲量變化分布特征。

圖4 長江流域2002年4月—2017年6月地表水儲量變化的趨勢(a)和周年項振幅絕對值(b)空間分布Fig.4 Trends and absolute values of annual signals’ amplitudes of water storage change in the Yangtze River Basin from April 2002 to June 2017

總的來看，長江流域大部分區域水儲量在2002年4月—2017年6月間保持著增加的趨勢，且在中下游區域增速較快，而金沙江流域中部大片區域以及岷江流域、嘉陵江流域、漢江流域北部等區域的水儲量處在持續減少的狀態。年際信號的振幅反映了水儲量周期變化的幅度，可以看到長江流域周年振幅較為平穩，統計后發現全流域平均振幅為64.9 mm，78.4%的格網振幅均小于100 mm，而在金沙江流域南部以及洞庭湖流域東部、鄱陽湖流域等中下游區域振幅較大，部分地區達到200 mm以上。

3.3 水儲量變化的水文氣象驗證

為了驗證GRACE水儲量變化結果的準確性，基于流域水量平衡方程，我們選擇長江流域上游2007—2016這10 a間的水儲量和水文氣象數據，探究GRACE衛星在流域尺度的觀測效果和適應性，同時分析水儲量變化的驅動因素。圖5為長江流域上游水儲量變化和降雨—蒸散—徑流的時間序列對比，發現二者的變化規律十分接近。2006—2016年，長江流域上游水量收支狀態基本平衡，即水儲量的變化量(-3.7 mm·a-1)≈輸入量(降雨0.4 mm·a-1)-輸出量(蒸散0.1 mm·a-1+徑流3.8 mm·a-1)，同時收支狀態基本平衡也反映出這些水文要素模型的準確性，它們在長江流域上游的這10 a間有著很好的擬合效果，總體上呈現出良好的評估效果。

圖5 長江流域上游2007年1月—2016年12月水量收支狀態Fig.5 The water storage changes of upper Yangtze River from January 2007 to December 2016

4 結論

本文利用GRACE重力衛星數據對長江流域地表水儲量變化進行了分析。結果表明：2002年4月—2017年6月近16 a間，長江流域整體水儲量處在一個緩慢上升的狀態，地表水儲量變化速率為4.1±1.1 mm·ɑ-1，總體水量在這一期間增加了大約1.1×1011T。在長江中下游區域，水儲量增速較快，而金沙江流域中部大片區域以及岷江流域、嘉陵江流域、漢江流域北部等區域的水儲量處在持續減少的狀態；流域大部分區域水儲量的年際波動較為平緩，78.4%的格網振幅均小于100 mm。降雨對長江流域水儲量變化周期項的影響程度較大。基于水文氣象數據的統計分析顯示這一結果能夠較為準確的體現長江流域近16 a的水量收支狀態。

目前GRACE衛星任務已結束，后續衛星GRACE Follow-On已發射，但數據尚未發布。未來將繼續關注GRACE Follow-On衛星數據，對長江流域地表水儲量變化進行更加詳細和深入的研究與分析，并研究與GRACE衛星空間時間段數據的填補問題，實現對長江流域陸地水儲量的長期監測。