利用COST-G GRACE時變模型反演黃土高原陸地水儲量變化

禤鍵豪 張興福 陳鑑華 馬資穎

1 廣東工業大學土木交通工程學院，廣州市外環西路100號，510006

GRACE衛星可為陸地水儲量變化監測、海平面變化監測等提供重要數據支撐。CSR、JPL、GFZ、TU Graz等科研機構陸續發布了多個版本的GRACE月時變重力場模型，但由于不同機構發布的模型所采用的反演方法和數據處理策略不同，導致各模型精度略有差異，特別是在GRACE任務末期，不同模型間差異較大。Meyer等[1]采用方差分量估計方法融合AIUB-RL02、CSR-RL06、GFZ-RL06、GRGS-RL04和ITSG-Grace2018共5個月時變模型得到COST-G GRACE模型，理論上該模型噪聲水平較小，可為利用GRACE時變模型反演區域陸地水儲量變化提供更好的數據支持。

黃土高原位于我國中北部，是我國四大高原之一。隨著黃土高原植被種植面積的增大，植被蒸騰作用會使該區域的陸地水儲量減少[2]。此外，降水、氣溫、農業用水、工業用水等因素也會影響區域陸地水儲量的變化[3]。本文綜合利用COST-G GRACE月時變重力場模型、實測降水、實測氣溫、GLDAS(global land data assimilation system)模型、實測淺層地下水、原煤開采量、歸一化植被指數(NDVI)等多源數據對黃土高原陸地水儲量變化進行分析，最后利用偏最小二乘回歸法分析該區域陸地水儲量變化的驅動因素。

1 原理與方法

1.1 GRACE模型反演陸地水儲量變化

利用時變重力場球諧模型計算區域質量變化對應的等效水柱高，公式為[4]：

(1)

式中，a為地球赤道處半徑，ρave為地球平均密度，ρ為水密度，lm為歸一化締合勒讓德多項式，kl是l階勒夫數，ωlm為平滑核函數，本文采用DDK3濾波[5]，Δlm與Δlm為球諧系數變化量。

通常陸地水儲量變化可近似表達為地表水儲量變化與地下水儲量變化之和[6]：

ΔTWS=ΔSW+ΔGW

(2)

式中，ΔTWS為陸地水儲量變化，ΔSW為地表水儲量變化，ΔGW為地下水儲量變化。利用GRACE時變重力場模型時間序列可以估計ΔTWS，利用GLDAS模型可以計算ΔSW，再利用式(2)即可得到ΔGW。

1.2 實測淺層地下水變化

將實測淺層地下水井監測數據轉換為水儲量變化，計算公式為[7]：

(3)

式中，ΔGWwell為實測淺層地下水儲量變化，Si為子區平均給水度，本文取黃土高原平均給水度為0.03[7]，Wi為子區面積，Δi為相應子區監測井的平均地下水位變化。

1.3 采煤耗水量

煤炭開采會影響區域陸地水儲量變化，利用式(4)可將原煤開采量轉換為采煤耗水量[3]：

TWScoal=μMc

(4)

式中，TWScoal為采煤用水量，μ為耗水系數，取為0.87 m3/t，Mc為原煤開采量。黃土高原內含陜西、山西、青海、寧夏、內蒙古、河南、甘肅等7個省或自治區，由于部分省或自治區并不完全位于黃土高原內，本文使用面積加權法獲取該省或自治區位于黃土高原區域內的原煤開采量。

1.4 偏最小二乘回歸

在偏最小二乘回歸中，自變量對因變量影響的重要程度可由變量投影重要性(VIP)測定。一般來說，當VIP大于0.8時，自變量對因變量的變化具有重要解釋意義，VIP計算公式為[3]：

(5)

式中，p為自變量個數，m為提取主成分個數，th為第h個主成分，Y為因變量，R(Y,th)為Y與th的相關系數，Whi為自變量在主成分th中的權重。本文取降水、氣溫、NDVI、采煤用水量、農業用水、工業用水、生活用水作為自變量，陸地水儲量作為因變量計算VIP值。

2 陸地水儲量變化時空分析

2.1 研究數據及處理流程

表1為用于研究黃土高原陸地水儲量時空變化特征及其驅動因素的有關數據，共收集到74口淺層地下水監測井的水位數據，其中，45口井收錄于《陜西省地下水監測年鑒》，其余收錄于《中國地質環境監測地下水位年鑒》，其點位見圖1。

表1 數據來源Tab.1 Data sources

黃土高原邊界矢量數據來源于文獻[8]圖1 黃土高原淺層地下水監測井位置分布Fig.1 Distribution of shallow groundwater monitoring wells in the Loess plateau

將COST-G GRACE月時變重力場模型(下文簡稱為GRACE模型)截斷至60階次，采用圖2中數據處理流程反演區域陸地水儲量變化。

圖2 GRACE模型數據處理流程Fig.2 Data processing flowchart of GRACE model

2.2 陸地水儲量時間變化分析

圖3為GRACE模型反演的陸地水儲量變化與氣溫異常量(當月值減去平均值)、降水量對比。結果顯示：1)氣溫、降水存在季節性變化特征，但陸地水儲量季節性變化特征不明顯；2)GRACE模型反演的陸地水儲量變化在2003年末、2013年夏季、2016年秋季出現明顯峰值，而在2003年全年山西、陜西降水較常年多20%以上[13]，2013年夏季降水較多，導致汾川河發生超歷史洪水、渭河發生超警洪水[14]，2016-07中下旬黃土高原內多地降水較常年偏多，導致山西等地發生暴雨洪澇災害[15]；3)GRACE模型反演的陸地水儲量變化在2007年、2009年、2015年和2016年夏季出現較為顯著的谷值，2007年夏季黃土高原部分地區由于降水較少使得旱情偏重且持續時間長[15]，2009年冬小麥主產區發生冬春連旱、西北部分地區發生伏旱和秋旱等[15]，2015年全年降水分布南多北少，特別是西北東部和黃淮等區域偏少10%～30%[14]，2016-01～05北方冬麥區發生冬旱，其中山西降水較常年同期偏少60%[15]。由此可見，黃土高原區域陸地水儲量變化與氣溫、降水變化密切相關，GRACE模型反演的陸地水儲量變化時間序列中較為明顯的峰值與谷值往往與當時或之前的洪旱災害相關。

圖3 GRACE模型反演的黃土高原ΔTWS與降水、氣溫變化對比Fig.3 Comparison of ΔTWS derived by GRACE and precipitation and temperature in the Loess plateau

圖4為黃土高原GRACE模型反演的陸地水儲量變化、GLDAS模型計算的地表水儲量變化、實測淺層地下水儲量變化以及GRACE模型估計的地下水儲量變化。由圖4可知，4條曲線變化較為一致。2002-04～2016-12黃土高原陸地水儲量變化與地表水變化、淺層地下水變化以及GRACE估計的地下水變化的相關系數分別為0.53、0.68、0.89，表明該區域陸地水儲量變化與地下水變化相關性最高。2002-04～2016-12黃土高原陸地水儲量變化大致可分為3個階段：2002-04～2003-12以3.64±3.51 cm/a的速率上升(由于該時間段較短，下文不作定量分析)；2004-01～2009-12以1.64±0.25 cm/a的速率快速下降；2010-01～2016-12以0.57±0.29 cm/a的速率緩慢下降。在陸地水儲量上升階段，地表水、實測淺層地下水和GRACE估計的地下水均具有明顯的上升趨勢；在陸地水儲量虧損的2個階段，地表水分別以-0.57±0.18 cm/a和0.08±0.17 cm/a的速率變化，實測淺層地下水下降速率分別為0.35±0.08 cm/a和0.46±0.06 cm/a，GRACE估計的地下水下降速率分別為1.07±0.20 cm/a和0.65±0.24 cm/a。2002-04～2016-12 GRACE估計的地下水與實測淺層地下水的相關系數達到0.71，而RMSE為2.86 cm，兩者具有較好的符合度，2種方法得到的地下水儲量均在2004年年初達到峰值。結合圖3和圖4可知，地表水、降水也在2003年中后期出現極大值，而氣溫較低，由此推測地下水儲量上升的主要原因為降水增加且蒸散量減少，經過一段時間的滲透作用使地下水出現滯后的上升信號，因此陸地水儲量也出現相似的上升趨勢。結合上述實測淺層地下水下降速率與GRACE估計的地下水虧損速率推斷，2004-01～2009-12黃土高原深層地下水也可能存在明顯的虧損。

圖4 GRACE模型反演的黃土高原ΔTWS、GLDAS模型計算的ΔSW、實測淺層ΔGWwell和GRACE模型估計的ΔGW變化對比Fig.4 Comparison of ΔTWS derived by GRACE model, ΔSW calculated by GLDAS model, measured shallow ΔGWwell and ΔGW estimated by GRACE model in the Loess plateau

2.3 陸地水儲量空間變化分析

圖5為GRACE反演的陸地水儲量、GRACE估計的地下水、GLDAS模型計算的地表水、NDVI、降水、氣溫等數據變化的空間趨勢圖?？梢钥闯?，GRACE模型反演的2002-04～2016-12陸地水儲量在山西存在較大虧損，虧損中心位于黃土高原中東部，速率超過1 cm/a。GLDAS模型計算的地表水在研究時間內并無明顯的空間變化趨勢。GRACE探測到的2002-04～2016-12陸地水儲量虧損信號與該區域地下水下降關聯較大，且地下水虧損的空間位置與陸地水儲量下降的位置相似。NDVI在黃土高原區域內具有較大的上升趨勢，在上升最大的黃土高原中東部區域，降水也具有較大的上升趨勢，但由于植被增多之后，地下水與地表水的補給會由于植被蒸騰耗水與冠層截留而降低，且該區域氣溫具有較大的上升趨勢，結合地下水開采等因素，該區域陸地水儲量顯示為下降趨勢。氣溫在研究區大部分區域均呈現上升趨勢，這會增大蒸散作用，對陸地水儲量虧損存在一定影響。

圖5 黃土高原2002-04～2016-12陸地水儲量、地下水、地表水、NDVI、降水、氣溫變化趨勢Fig.5 Variation trend of TWS, GW, SW, NDVI, precipitation and temperature in the Loess plateau from April 2002 to December 2016

3 陸地水儲量變化驅動因素分析

GRACE反演黃土高原陸地水儲量的虧損大致分為2個階段：2004-01～2009-12(時段1)呈現較大的虧損趨勢，2010-01～2016-12(時段2)表現為較為緩慢的虧損趨勢。因此，基于這2個階段分析陸地水儲量變化的驅動因素。其中，時段2原煤開采量數據大部分缺失，故暫不考慮該時段該因素的影響。

降水、氣溫、NDVI、采煤用水、農業用水、工業用水、生活用水等可直接或間接影響黃土高原陸地水儲量變化，以上述7個因素作為自變量，以陸地水儲量變化作為因變量，進行偏最小二乘回歸分析，表2為各驅動因子的VIP值。在時段1中，采煤用水、生活用水、NDVI的VIP值大于0.8，均為陸地水儲量下降的重要驅動因素，其中，采煤用水與生活用水的VIP值為1.60，這2個因素均與人類活動有直接關系。由于退耕還林工程的推進，人工林面積急劇增加[2]，NDVI在時段1中以每年0.003±0.002的速率上升，導致該區域被植被吸收的地表水和地下水增加。因此，在時段1，人類活動與植被作用對陸地水儲量的下降起到推進作用。在時段2中，VIP值大于0.8的因素為工業用水、生活用水、氣溫，該時段中人類活動仍對陸地水儲量的下降起到重要作用。但值得注意的是，氣溫上升也為黃土高原陸地水儲量虧損的驅動因素：時段2中黃土高原地區的氣溫以0.16±0.14 ℃/a的速率上升，而NDVI以每年0.004±0.002的速率上升，黃土高原氣溫顯著上升使得蒸發作用增大，加上植被增多會促進蒸騰作用，最終導致陸地水儲量減少。

表2 偏最小二乘回歸計算的驅動因子的VIP值Tab.2 VIP values of driving factors obtained by partial least squares regression

4 結 語

本文利用COST-G GRACE月時變重力場模型反演黃土高原陸地水儲量變化，并與實測降水、實測氣溫、GLDAS模型計算的地表水變化、實測淺層地下水變化等進行對比，分析陸地水儲量的時空變化特征，最后使用偏最小二乘回歸分析陸地水儲量變化的驅動因素，得到以下結論：

1)2002-04～2016-12黃土高原陸地水儲量具有上升-下降-平緩下降的變化特征。其中，山西陸地水儲量虧損最為嚴重，最大虧損速率超過1 cm/a。對比實測淺層地下水與GRACE估計的地下水可知，兩者相關系數達到0.71，符合度較高，通過分析認為黃土高原陸地水儲量變化與地下水變化密切相關。

2)經偏最小二乘回歸分析可知，2004～2009年采煤用水和生活用水等人類活動以及植被作用是導致黃土高原陸地水儲量虧損的重要驅動因素；2010～2016年工業用水、生活用水以及氣溫變化是黃土高原陸地水儲量虧損的重要驅動因素。

致謝：感謝陳劍利教授、馮偉教授和ICGEM提供模型與軟件。