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

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

1 廣東工業(yè)大學(xué)土木交通工程學(xué)院，廣州市外環(huán)西路100號，510006

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

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

1 原理與方法

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

利用時變重力場球諧模型計算區(qū)域質(zhì)量變化對應(yīng)的等效水柱高，公式為[4]：

(1)

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

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

ΔTWS=ΔSW+ΔGW

(2)

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

1.2 實測淺層地下水變化

將實測淺層地下水井監(jiān)測數(shù)據(jù)轉(zhuǎn)換為水儲量變化，計算公式為[7]：

(3)

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

1.3 采煤耗水量

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

TWScoal=μMc

(4)

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

1.4 偏最小二乘回歸

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

(5)

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

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

2.1 研究數(shù)據(jù)及處理流程

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

表1 數(shù)據(jù)來源Tab.1 Data sources

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

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

圖2 GRACE模型數(shù)據(jù)處理流程Fig.2 Data processing flowchart of GRACE model

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

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

圖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估計的地下水變化的相關(guān)系數(shù)分別為0.53、0.68、0.89，表明該區(qū)域陸地水儲量變化與地下水變化相關(guān)性最高。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估計的地下水與實測淺層地下水的相關(guān)系數(shù)達到0.71，而RMSE為2.86 cm，兩者具有較好的符合度，2種方法得到的地下水儲量均在2004年年初達到峰值。結(jié)合圖3和圖4可知，地表水、降水也在2003年中后期出現(xiàn)極大值，而氣溫較低，由此推測地下水儲量上升的主要原因為降水增加且蒸散量減少，經(jīng)過一段時間的滲透作用使地下水出現(xiàn)滯后的上升信號，因此陸地水儲量也出現(xiàn)相似的上升趨勢。結(jié)合上述實測淺層地下水下降速率與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、降水、氣溫等數(shù)據(jù)變化的空間趨勢圖?？梢钥闯?，GRACE模型反演的2002-04～2016-12陸地水儲量在山西存在較大虧損，虧損中心位于黃土高原中東部，速率超過1 cm/a。GLDAS模型計算的地表水在研究時間內(nèi)并無明顯的空間變化趨勢。GRACE探測到的2002-04～2016-12陸地水儲量虧損信號與該區(qū)域地下水下降關(guān)聯(lián)較大，且地下水虧損的空間位置與陸地水儲量下降的位置相似。NDVI在黃土高原區(qū)域內(nèi)具有較大的上升趨勢，在上升最大的黃土高原中東部區(qū)域，降水也具有較大的上升趨勢，但由于植被增多之后，地下水與地表水的補給會由于植被蒸騰耗水與冠層截留而降低，且該區(qū)域氣溫具有較大的上升趨勢，結(jié)合地下水開采等因素，該區(qū)域陸地水儲量顯示為下降趨勢。氣溫在研究區(qū)大部分區(qū)域均呈現(xiàn)上升趨勢，這會增大蒸散作用，對陸地水儲量虧損存在一定影響。

圖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 陸地水儲量變化驅(qū)動因素分析

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

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

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

4 結(jié) 語

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

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

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

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