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

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

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

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

1 原理與方法

1.1 GRACE模型反演陸地水儲(chǔ)量變化

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

(1)

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

通常陸地水儲(chǔ)量變化可近似表達(dá)為地表水儲(chǔ)量變化與地下水儲(chǔ)量變化之和[6]：

ΔTWS=ΔSW+ΔGW

(2)

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

1.2 實(shí)測(cè)淺層地下水變化

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

(3)

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

1.3 采煤耗水量

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

TWScoal=μMc

(4)

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

1.4 偏最小二乘回歸

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

(5)

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

2 陸地水儲(chǔ)量變化時(shí)空分析

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

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

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

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

將COST-G GRACE月時(shí)變重力場(chǎng)模型(下文簡(jiǎn)稱為GRACE模型)截?cái)嘀?0階次，采用圖2中數(shù)據(jù)處理流程反演區(qū)域陸地水儲(chǔ)量變化。

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

2.2 陸地水儲(chǔ)量時(shí)間變化分析

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

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

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

圖4 GRACE模型反演的黃土高原ΔTWS、GLDAS模型計(jì)算的ΔSW、實(shí)測(cè)淺層ΔGWwell和GRACE模型估計(jì)的ΔGW變化對(duì)比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 陸地水儲(chǔ)量空間變化分析

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

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

3 陸地水儲(chǔ)量變化驅(qū)動(dòng)因素分析

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

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

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

4 結(jié) 語(yǔ)

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

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

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

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