顧及GRACE季節(jié)影響的華北平原水儲(chǔ)量變化反演

李 圳，章傳銀，柯寶貴，劉 陽，李婉秋，尹 財(cái)

1. 山東科技大學(xué)測(cè)繪科學(xué)與工程學(xué)院，山東 青島 266510； 2. 中國(guó)測(cè)繪科學(xué)研究院，北京 100830

近年來，GRACE(gravity recovery and climate experiment)在監(jiān)測(cè)陸地水儲(chǔ)量時(shí)變、海平面變化和冰蓋質(zhì)量變化等領(lǐng)域展示出了巨大的應(yīng)用潛力。在監(jiān)測(cè)陸地水儲(chǔ)量時(shí)變方面，文獻(xiàn)[1]利用GRACE RL04數(shù)據(jù)成功監(jiān)測(cè)到了2005年亞馬遜流域的干旱。文獻(xiàn)[2—3]利用GRACE和水文模式數(shù)據(jù)發(fā)現(xiàn)了印度北部地下水超采引起的陸地水變化。文獻(xiàn)[4]利用GRACE數(shù)據(jù)估計(jì)2006—2009年間美國(guó)加利福尼亞州地下水含量總共減少了(31±3)km3。此外，我國(guó)學(xué)者對(duì)三峽、華北平原、長(zhǎng)江黃河流域等地區(qū)陸地水儲(chǔ)量變化的GRACE監(jiān)測(cè)開展了廣泛研究[5-10]。但這些研究對(duì)于影響陸地水儲(chǔ)量變化的因素均缺少系統(tǒng)的分析。

GRACE月重力場(chǎng)產(chǎn)品不僅本身包含有誤差，而且在對(duì)其進(jìn)行后處理時(shí)也會(huì)產(chǎn)生誤差，通常稱之為后處理誤差[11]。例如對(duì)球諧系數(shù)的去相關(guān)處理、高斯濾波平滑及研究區(qū)域問題時(shí)的信號(hào)“泄漏”等。目前，主要使用尺度因子來減少其后處理誤差的影響，對(duì)于尺度因子的確定方法，不同學(xué)者提出了不同的意見。文獻(xiàn)[12]提出了格網(wǎng)尺度因子的方法，試驗(yàn)表明這種方法與傳統(tǒng)的區(qū)域尺度因子方法結(jié)果的精度高低取決于是否位于各自的適用區(qū)域。文獻(xiàn)[13]提出一種基于不同時(shí)刻水文數(shù)據(jù)作為先驗(yàn)信息的尺度因子確定方法，但該方法過于依賴水文模型。文獻(xiàn)[14—15]使用單一尺度因子結(jié)合區(qū)域核函數(shù)來減少數(shù)據(jù)后處理影響，在亞馬遜流域和華北地區(qū)也取得了不錯(cuò)的效果。水儲(chǔ)量的變化通常包含長(zhǎng)周期和季節(jié)性變化的趨勢(shì)項(xiàng)，文獻(xiàn)[12]指出，如果能將這兩種趨勢(shì)項(xiàng)分開，使用不同的尺度因子處理，可能會(huì)得到更高精度的處理結(jié)果。

本文基于移去-恢復(fù)方法，利用GRACE Level-2數(shù)據(jù)研究了華北平原2003年1月—2014年6月的陸地水變化，提出一種新的顧及季節(jié)影響尺度因子的計(jì)算方法，用于減小GRACE后處理誤差。該方法對(duì)水文模型的依賴性不強(qiáng)，而且還保留了尺度因子的季節(jié)性變化特征。將本文方法所得結(jié)果與水文模式和降水模型結(jié)合進(jìn)行比較分析，并使用交叉小波譜分析了陸地水儲(chǔ)量變化較降水?dāng)?shù)據(jù)的滯后現(xiàn)象。

1 數(shù)據(jù)與處理方法

1.1 GRACE數(shù)據(jù)

本文使用美國(guó)得克薩斯大學(xué)空間研究中心(Center for Space Research，University of Texas at Austin，CSR)提供的2003年1月—2014年6月共138個(gè)月的GRACE Level-2(Release-5)數(shù)據(jù)。期間缺失數(shù)據(jù)的月份采用多年平均值補(bǔ)齊。該數(shù)據(jù)產(chǎn)品已將非潮汐的海洋、大氣信號(hào)，以及各種潮汐影響扣除[16]。本文利用ECCO地心運(yùn)動(dòng)模型估計(jì)GRACE數(shù)據(jù)中的1階球諧系數(shù)[17]；利用SLR提供的C20項(xiàng)替代基于GRACE軌道解算的C20項(xiàng)[18]；采用改進(jìn)的P3M9去相關(guān)濾波去除“南-北條帶”誤差[19-20]；利用250 km扇形濾波來平滑高階球諧系數(shù)的噪聲[21]。為消除地球長(zhǎng)周期的影響和平均地球重力場(chǎng)的影響，對(duì)這138個(gè)月的地球重力場(chǎng)模型的球諧系數(shù)求平均，再用各個(gè)月的球諧系數(shù)減去這一平均值，得到球諧系數(shù)變化(ΔC、ΔS)。陸地水質(zhì)量的變化是引起重力場(chǎng)變化的主要原因[22]，陸地水儲(chǔ)量的變化Δh可用式(1)計(jì)算[23]

ΔSlmsin(mλ))Plm(cosθ)

(1)

式中，θ和λ分別表示計(jì)算點(diǎn)的余緯和經(jīng)度；ρe表示地球的平均密度；ρw為水的平均密度；a為赤道半徑；kl為負(fù)荷Love數(shù)；Plm(cosθ)表示l階m次的完全規(guī)格化勒讓德函數(shù)。高斯平滑函數(shù)(Wl、Wm)可以通過遞推公式(2)獲取[23]

(2)

1.2 水文數(shù)據(jù)

本文采用的水文模型是美國(guó)國(guó)家海洋和大氣局氣象預(yù)報(bào)中心的CPC(Climate Prediction Center)模型和美國(guó)宇航局哥達(dá)德飛行中心提供的全球陸地資料同化系統(tǒng)(global land data assimilation system，GLDAS)。CPC水文模型是根據(jù)全球觀測(cè)到的降水分布而建立的，主要反映表層土壤水含量和積雪變化。其時(shí)間間隔為1個(gè)月，經(jīng)緯度格網(wǎng)為0.5°×0.5°。GLDAS水文模型采用了NOAH、VIC和CLM等多種地表模型。本文采用基于NOAH地面模型的GLDAS，主要反映的是4層土壤水和冰雪變化，空間分辨率為1°×1°。本文只選用與GRACE數(shù)據(jù)相同時(shí)間范圍的數(shù)據(jù)。使用FFT(fast Fourier transformation)算法將格網(wǎng)形式的水文數(shù)據(jù)做球諧展開[24]，截取至60階，并同樣采用改進(jìn)P3M9去相關(guān)濾波和250 km扇形濾波的組合方法對(duì)其濾波。

1.3 降水?dāng)?shù)據(jù)

為了研究華北地區(qū)水儲(chǔ)量與降水的關(guān)系，使用了全球降水氣候?qū)W中心(Global Precipitation Climatology Center，GPCC)提供的降水?dāng)?shù)據(jù)，這一數(shù)據(jù)集是基于全球分布的約67 200個(gè)監(jiān)測(cè)站10年以上的觀測(cè)數(shù)據(jù)建立的。其時(shí)間間隔為一個(gè)月，格網(wǎng)分辨率為1°×1°。由于數(shù)據(jù)獲取有限，本文所分析使用的是2003年1月—2013年12月之間的數(shù)據(jù)。

2 顧及季節(jié)影響的尺度因子確定

為了更精確地利用GRACE研究陸地水儲(chǔ)量變化，通常使用水文模型數(shù)據(jù)來估計(jì)GRACE數(shù)據(jù)后處理導(dǎo)致的信號(hào)變化。水儲(chǔ)量的變化除了長(zhǎng)期的變化趨勢(shì)外，還包含了季節(jié)性變化，這兩種趨勢(shì)可能會(huì)有不同的尺度因子，在這種情況下使用單一的尺度因子就難以較好地恢復(fù)數(shù)據(jù)后處理導(dǎo)致的信號(hào)變化。

以華北平原為例，圖1給出了利用CPC水文模型展開的球諧系數(shù)計(jì)算的表層水儲(chǔ)量變化在濾波前后的時(shí)間序列?？梢詮膱D1中明顯地發(fā)現(xiàn)，華北平原的表層水儲(chǔ)量具有明顯的季節(jié)性變化特征，大約每年的春、秋季節(jié)，也就是圖中時(shí)間序列的谷值和峰值附近，CPC數(shù)據(jù)濾波前后的差異較其他月份大。使用了固定的尺度因子k后可以一定程度地減小這一差異，但在某些月份例如2007年秋季、2011年春季效果仍然不佳(圖1中黃色矩形框)。

該結(jié)果一定程度上說明了長(zhǎng)期變化趨勢(shì)和季節(jié)變化趨勢(shì)的尺度因子存在差異性，若只使用單一的尺度因子，那么季節(jié)變化對(duì)于結(jié)果的影響會(huì)比較大。文獻(xiàn)[25]利用以下方法來減小季節(jié)變化對(duì)GRACE后處理結(jié)果的影響：首先利用正向建模方法(forward modeling)確定水儲(chǔ)量變化的長(zhǎng)期趨勢(shì)，然后將這一長(zhǎng)期趨勢(shì)移去，并利用水文模型確定尺度因子，從而減少GRACE后處理過程中的振幅衰減和泄漏誤差。這一過程可以用式(3)表示

h(t)=khdetrend(t)+at

(3)

式中，h(t)是t時(shí)刻的水儲(chǔ)量；k是利用水文模型求出的尺度因子；hdetrend(t)是t時(shí)刻去除長(zhǎng)周期項(xiàng)之后的h(t)；a即為從正演模型中求出的長(zhǎng)周期趨勢(shì)項(xiàng)。在分析區(qū)域GRACE后處理影響時(shí)，也可以通過使用不同的尺度因子將長(zhǎng)期變化趨勢(shì)和季節(jié)變化趨勢(shì)分開處理，即

Δh=klongΔhlong+kseasonΔhseason

(4)

式中，Δhlong和Δhseason分別表示水儲(chǔ)量的長(zhǎng)期變化趨勢(shì)量和季節(jié)性變化量；klong和kseason分別為長(zhǎng)周期尺度因子和季節(jié)性尺度因子。為了分別求取長(zhǎng)周期尺度因子和季節(jié)性尺度因子，減少季節(jié)變化對(duì)水儲(chǔ)量反演的影響，本文使用一種新的基于移去-恢復(fù)方法的季節(jié)尺度因子計(jì)算方法，可以通過以下步驟實(shí)現(xiàn)：

(1) 時(shí)間序列分析確定長(zhǎng)期趨勢(shì)項(xiàng)與周期項(xiàng)。

(2) 分別移去濾波前后數(shù)據(jù)的長(zhǎng)期趨勢(shì)項(xiàng)與周期項(xiàng)，計(jì)算長(zhǎng)周期尺度因子和季節(jié)性尺度因子。

(3) 利用上一步計(jì)算得到的尺度因子恢復(fù)數(shù)據(jù)后處理過程對(duì)信號(hào)的影響和移去過程中的剩余殘差。

具體的計(jì)算方法是：首先利用CPC模型分析了華北平原表層水儲(chǔ)量變化的長(zhǎng)期趨勢(shì)，基于均方根值(root mean square，RMS)最小原則進(jìn)行線性擬合。圖1中給出了CPC數(shù)據(jù)未濾波結(jié)果的分段一次線性擬合結(jié)果，可以發(fā)現(xiàn)該地區(qū)的表層水儲(chǔ)量在2005年之前呈上升趨勢(shì)，在2005—2007年間呈下降趨勢(shì)，而2007年之后表層水儲(chǔ)量又重新變?yōu)樯仙厔?shì)，這一變化的趨勢(shì)與降水有著密切聯(lián)系。由于CPC的趨勢(shì)項(xiàng)在2005年初和2008年初發(fā)生了改變，將這11年水儲(chǔ)量變化的長(zhǎng)期趨勢(shì)項(xiàng)看作是相同的并不妥當(dāng)，所以本文采用分段處理的方法，即2003年1月—2005年6月為第1段，2005年7月—2007年12月為第2段，2008年1月—2014年12月為第3段。

分別對(duì)CPC未濾波和濾波之后結(jié)果進(jìn)行連續(xù)小波變換，分析其季節(jié)性變化的周期。圖2、圖3分別給出了CPC未濾波和濾波之后的連續(xù)小波變換能量譜。黑色粗實(shí)線圈內(nèi)表示通過了95%置信水平的紅噪聲檢驗(yàn)，黑色細(xì)實(shí)線下方是數(shù)據(jù)邊緣效應(yīng)影響較大的小波影響錐區(qū)域[26]。結(jié)合圖2、圖3可以發(fā)現(xiàn)，濾波并沒有改變數(shù)據(jù)的能量分布，CPC數(shù)據(jù)在濾波前后均有一個(gè)十分明顯的12個(gè)月的周期T1，對(duì)應(yīng)季節(jié)的變化，此外還有一個(gè)24個(gè)月左右的能量相對(duì)較弱的周期項(xiàng)T2。

然后，本文使用了一種移去-恢復(fù)技術(shù)，即分別移去濾波前后的水儲(chǔ)量變化長(zhǎng)周期項(xiàng)和季節(jié)性周期項(xiàng)，分析計(jì)算尺度因子，然后利用尺度因子恢復(fù)數(shù)據(jù)后處理過程對(duì)信號(hào)的影響并恢復(fù)移去過程中剩余的殘差項(xiàng)。那么，可以使用下式來表示CPC的表層水儲(chǔ)量變化

(5)

式中，A表示參考時(shí)刻t0的Δh；B表示長(zhǎng)期變化的線型趨勢(shì)項(xiàng)；Δt表示第i個(gè)月ti減去t0；C和S分別表示周期項(xiàng)的振幅，其下標(biāo)a和b分別對(duì)應(yīng)于周期T1和T2；ε表示移去過程中的剩余殘差。

使用式(5)擬合CPC濾波前的水儲(chǔ)量變化，時(shí)間序列擬合參數(shù)記為a1、b1、C1a、S1a、C1b和S1b。同理，用式(5)擬合CPC濾波后的水儲(chǔ)量變化時(shí)間序列，擬合參數(shù)記為a2、b2、C2a、S2a、C2b和S2b。將擬合的時(shí)間序列從原序列中“移出”，并利用b1=b2k1、C1a=C2ak2、S1a=S2ak3、C1b=C2bk4和S1b=S2bk5求得長(zhǎng)周期尺度因子k1與季節(jié)性尺度因子k2、k3、k4和k5，這就是移去過程。這一過程移去的序列與原序列之間的差值即為剩余殘差ε。

最后，利用移去過程求得的各尺度因子，通過下式恢復(fù)數(shù)據(jù)后處理過程對(duì)信號(hào)的影響和移去過程中的剩余殘差

(6)

使用單一尺度因子(filter+k)、顧及季節(jié)性影響的多尺度因子恢復(fù)濾波后的CPC數(shù)據(jù)(filter+kn)，與未濾波的CPC數(shù)據(jù)的時(shí)間序列(original)對(duì)比如圖4所示。未使用尺度因子(filter)、filter+k和filter+kn與original之間的差值如圖5所示。綜合兩圖可以看出，顧及季節(jié)性影響的多尺度因子可以進(jìn)一步恢復(fù)后處理影響(圖5中藍(lán)色矩形框)，較單一尺度因子法有明顯提高。綜合來看2003—2014這11年的數(shù)據(jù)，可以發(fā)現(xiàn)本文所提出的方法恢復(fù)的表層水儲(chǔ)量變化與濾波之前序列的偏差基本是在±1 cm等效水高這一區(qū)間內(nèi)浮動(dòng)，結(jié)果更加穩(wěn)定。計(jì)算求得filter、filter+k、filter+kn方法的RMS分別為13.3、9.9和6.9 mm。本文所提出的方法分別比前兩種方法的RMS降低了48%和30%左右。需要說明的是，由于是利用CPC進(jìn)行模擬，當(dāng)使用Klees方法對(duì)CPC水文模式進(jìn)行改正時(shí)，RMS為0。本文將利用CPC數(shù)據(jù)求得的尺度因子應(yīng)用于GLDAS水文模式數(shù)據(jù)，求得filter+k和filter+kn方法的RMS分別為7.8 mm和7.1 mm，而Klees方法的RMS為25.1 mm。這一結(jié)論說明了filter+k和filter+kn方法對(duì)水文模型的依賴性均相對(duì)較小，而filter+kn方法RMS更小，說明可以更好地降低GRACE數(shù)據(jù)的后處理誤差。Klees方法在沒有可靠水文模式數(shù)據(jù)的基礎(chǔ)上，結(jié)果并不可靠，這與Klees等的結(jié)論一致[13]。

為驗(yàn)證使用了移去-恢復(fù)技術(shù)的顧及季節(jié)性變化的多尺度因子并沒有改變?cè)夹盘?hào)的頻譜特征，及引入新的誤差，圖6給出了這一方法恢復(fù)的等效水高變化的連續(xù)小波變換能量譜?？梢园l(fā)現(xiàn)，使用該方法恢復(fù)的序列頻譜能量分布與未濾波的CPC數(shù)據(jù)基本相同，均表現(xiàn)出具有一個(gè)12個(gè)月的明顯的周期和24個(gè)月左右的能量相對(duì)較弱的周期。頻譜能量并沒有增加，只在高頻上有幅度較弱的減小趨勢(shì)，說明了該方法的RMS主要來源于高頻部分。這在一定程度上驗(yàn)證了本文方法的可靠性。

3 GRACE反演華北平原陸地水量變化

使用上節(jié)所述的基于移去-恢復(fù)技術(shù)的顧及季節(jié)性變化的多尺度因子恢復(fù)GRACE數(shù)據(jù)后處理影響，利用2003年1月—2014年6月的GRACE數(shù)據(jù)和CPC數(shù)據(jù)，估計(jì)了華北平原的陸地水的月變化，并結(jié)合GPCC月降水量數(shù)據(jù)做了進(jìn)一步的分析?？紤]GRACE的分辨率大致為300 km，本文的研究區(qū)域?yàn)?5.5°N—42.5°N，114.5°E—120.5°E范圍內(nèi)的陸地區(qū)域，主要由華北平原、燕山山脈和山東丘陵的部分地區(qū)組成。溫帶地區(qū)植物根區(qū)表層水儲(chǔ)量對(duì)陸地水量的變化貢獻(xiàn)顯著，本文利用GRACE估計(jì)的陸地水量扣除CPC估計(jì)的地表水變化，就可以得到華北平原的地下水變化情況。由于軌道調(diào)整等原因，2003年前幾個(gè)月的GRACE精度較差[27]，為了得到更精確的水儲(chǔ)量變化值，下文在采用線性擬合時(shí)未采用上節(jié)的分段方法，只以2008年作為分界線將數(shù)據(jù)分為了兩段擬合?？紤]到華北平原的水儲(chǔ)量季節(jié)性變化顯著，2014年7—12 月的數(shù)據(jù)未參與解算，所以在擬合第2段時(shí)使用的數(shù)據(jù)截至2013年12月。取擬合結(jié)果2倍的中誤差作為變化速率的不確定度[5]。

圖7和圖8分別為華北平原陸地水量變化趨勢(shì)和地下水儲(chǔ)量變化。從圖7中可以看出，GRACE反演的陸地水量與CPC模型反演的地表水均表現(xiàn)出季節(jié)性特征，主要表現(xiàn)為春季和秋季達(dá)到最小和最大值，但地表水的季節(jié)性特征更加明顯，而陸地水儲(chǔ)量變化波動(dòng)較大，這主要與GRACE的觀測(cè)誤差和后處理過程中引入的誤差有關(guān)。GRACE資料與CPC資料得到的結(jié)果在振幅上存在差異，主要是因?yàn)镚RACE得到的是綜合各種因素后的水儲(chǔ)量變化，而CPC則主要反映表層土壤水和積雪的變化，并沒有包含地下水等的變化。二者在2008年之前符合較好，這一時(shí)期發(fā)生了多次比較嚴(yán)重的干旱災(zāi)害，所以陸地水和地表水分別以(-7.9±2.4)mm/a和(-7.3±2.8)mm/a的速度下降。如圖8所示，在這一期間地下水水量基本保持不變，僅以(-0.6±1.4)mm/a的速度減少。而在2008年以后，GRACE與CPC的差異逐步增加，陸地水和地表水分別以(4.4±1.3)mm/a和(10.9±2.1)mm/a的速度上升。這一差異產(chǎn)生的原因主要是因?yàn)槿A北平原植物稀少且根系不發(fā)達(dá)，地表蒸發(fā)量大，難以保留水分[27]。這一期間超采的地下水幾乎全部流失，導(dǎo)致了地下水下降，但地下水下降速度有減緩趨勢(shì)。通過對(duì)2008—2014年間的數(shù)據(jù)進(jìn)行二次多項(xiàng)式擬合發(fā)現(xiàn)，地下水平均以(-6.5±1.6)mm/a的速度下降，下降速度以0.9 mm/a2的趨勢(shì)減小。在整個(gè)研究期間內(nèi)，華北平原的陸地水和地表水整體分別呈現(xiàn)(-2.0±0.6)mm/a和(2.9±0.7)mm/a的變化趨勢(shì)，而地下水整體呈現(xiàn)(-4.8±0.7)mm/a的下降趨勢(shì)。

結(jié)合圖9中GPCC提供的降水?dāng)?shù)據(jù)，可以發(fā)現(xiàn)華北平原的降水量在2008年之前整體呈現(xiàn)減小趨勢(shì)，但在2008年之后降水量開始增加，尤其是2010年之后增加顯著，再加上南水北調(diào)工程的逐步推進(jìn)，一定程度上解釋了2008年前后華北平原水儲(chǔ)量變化產(chǎn)生趨勢(shì)改變的原因，也說明了降水是影響華北平原水儲(chǔ)量變化的重要因素。華北平原地區(qū)在每年的上半年大量依靠地下水來灌溉農(nóng)作物[28]，特別是每年的4、5月份，正值小麥生長(zhǎng)的灌漿期，需水量大，作物灌溉用水占到了總用水量的57%[29]，因此，圖8中可以看到地下水在每年的前半段呈現(xiàn)下降趨勢(shì)。而由于夏季降水增多，下半年的農(nóng)業(yè)用水更多地依賴于降雨，地下水的開采量相應(yīng)減小[30]，因此下半年地下水呈現(xiàn)一定程度的回升。雖然2008年之后夏季降水量增多，但冬春季節(jié)降水仍未有明顯變化，春旱災(zāi)害仍然頻發(fā)，期間地下水的開采量劇增，這一點(diǎn)可以通過圖8中冬春季節(jié)地下水變化的振幅上看出。因?yàn)槎航邓^少，可以認(rèn)為地下水的變化主要是由人類開采引起的?？梢园l(fā)現(xiàn)2008年之后地下水的振幅以(12.9±1.2)mm/a的速度增大，這也說明了地下水的季節(jié)性變化較之前更為顯著，冬春季節(jié)的地下水消耗量更大。夏季降水是地下水的主要補(bǔ)充來源，上述研究表明夏季地下水的補(bǔ)充量尚不足以彌補(bǔ)冬秋季節(jié)的開采量，且補(bǔ)充量與開采量的差異在2008年后增加顯著，但這種差異隨著降水的持續(xù)增多也會(huì)減小。綜合上述因素，也就解釋了2008年之后地下水加速下降，但下降速度以0.9 mm/a2的趨勢(shì)減小的現(xiàn)象。

通過對(duì)比圖9中的GRACE陸地水儲(chǔ)量和GPCC降水量，發(fā)現(xiàn)降水量大于水儲(chǔ)量的變化，這是因?yàn)榻邓坎粫?huì)完全轉(zhuǎn)化為區(qū)域的陸地水，還有一部分通過地表徑流，蒸發(fā)作用流失。還可以發(fā)現(xiàn)，降水量和陸地水的變化有一個(gè)明顯的相位差。產(chǎn)生這一相位差的主要原因是：

(1) GRACE觀測(cè)到的陸地水包含了保留在本地區(qū)該月之前的降水，而降水又需要一定的時(shí)間通過滲透等作用轉(zhuǎn)化成華北區(qū)域的陸地水儲(chǔ)量。

(2) 根據(jù)水量平衡方程：GRACE前后月份之差=降雨-蒸發(fā)-徑流，華北平原的降水主要集中在夏季，夏季氣溫高，蒸發(fā)作用十分劇烈，同時(shí)降水增多，地表徑流也會(huì)相應(yīng)增大[31]。

為了進(jìn)一步分析這一相位變化，本文對(duì)陸地水儲(chǔ)量變化和降水量的時(shí)間序列進(jìn)行了交叉小波分析，結(jié)果如圖10所示。圖中箭頭方向反映GRACE與GPCC的相位關(guān)系，箭頭向右表示相位相同，向左表示相位反向，箭頭向上表示GRACE比GPCC滯后1/4周期，向下則表示GRACE比GPCC提前1/4周期[26]。

圖10表明，GRACE與GPCC有一個(gè)約12個(gè)月的共振周期。2008年之前，在這一周期上GRACE的相位比GPCC平均滯后77°，這也就是說GRACE需要大概2.6個(gè)月的時(shí)間來響應(yīng)降水的變化。但在2008—2011年期間，GRACE與GPCC在12個(gè)月的周期上的共振能量減小，相位差也發(fā)生了比較大的改變。2003—2008年期間，地下水的變化并不顯著，GRACE主要反映的是地表水的變化，因而GRACE與降水量在12個(gè)月的周期上具有較大的共振能量，表現(xiàn)出比較明顯的季節(jié)性特征。2008年是華北平原的地下水開采量開始顯著增加的年份，尤其是2009—2010年，由圖7可以發(fā)現(xiàn)，GRACE的谷值對(duì)應(yīng)了CPC的峰值，說明GRACE主要反映的是地下水儲(chǔ)量的變化。這一期間人類生產(chǎn)活動(dòng)加劇了地下水的降低，一定程度上破壞了華北平原的陸地水量變化的季節(jié)性特征。2010年之后，地下水的開采量沒有減小趨勢(shì)，而降水量的增加對(duì)地下水的補(bǔ)充作用更為明顯，在這一時(shí)期，地下水的下降趨勢(shì)趨于穩(wěn)定，GRACE的季節(jié)性變化特征也就逐漸恢復(fù)。但GRACE的相位滯后現(xiàn)象更為顯著，平均達(dá)到了136°，即4.6個(gè)月，這種滯后現(xiàn)象的產(chǎn)生與地下水的超采、生產(chǎn)活動(dòng)取用水量增加以及地面建筑增多等有一定的關(guān)系，說明了人類活動(dòng)是華北平原陸地水儲(chǔ)量變化的另一重要決定因素。

圖1 CPC未濾波與濾波之后時(shí)間序列對(duì)比Fig.1 Time series comparison of the CPC before and after filtered

圖2 CPC濾波前時(shí)間序列能量譜Fig.2 Energy spectrum of CPC time series before filter

圖3 CPC濾波后時(shí)間序列能量譜Fig.3 Energy spectrum of CPC time series after filter

圖4 兩種尺度因子法恢復(fù)后處理影響效果對(duì)比Fig.4 Result comparison of two scale factor methods on restore post-processing effects

圖5 與未濾波表層水儲(chǔ)量之間的差值Fig.5 Difference from the unfiltered surface water storage

圖6 顧及季節(jié)影響尺度因子恢復(fù)濾波后時(shí)間序列的能量譜Fig.6 Energy spectrum of the recovered filtered time series using the scale factors considered seasonal influence

圖7 華北平原陸地水量變化趨勢(shì)Fig.7 Trends of land water storage change in North China Plain

圖8 華北平原地下水儲(chǔ)量變化Fig.8 Variation of groundwater storage in North China Plain

圖9 華北平原水儲(chǔ)量與降水量變化趨勢(shì)Fig.9 Trends of water storage and precipitation in North China Plain

圖10 GRACE水儲(chǔ)量與降水量的交叉小波譜分析Fig.10 Cross-spectral analysis of GRACE and precipitation

4 討 論

目前利用GRACE衛(wèi)星解算的時(shí)變重力場(chǎng)在高階項(xiàng)上信噪比較低。為了提高信噪比，減小誤差，通常會(huì)引入濾波技術(shù)，但在濾波過程中難免會(huì)削弱信號(hào)振幅。為了恢復(fù)被削弱的振幅，本文基于移去-恢復(fù)方法，提出了一種新的顧及季節(jié)影響尺度因子計(jì)算方法。該方法比單一尺度因子法[13]和Klees方法[14]表現(xiàn)更好，且比文獻(xiàn)[25]中的計(jì)算方法簡(jiǎn)便。然而，雖然文中將CPC計(jì)算求得的尺度因子應(yīng)用于GLDAS時(shí)仍可以得到比較好的結(jié)果，但是實(shí)際上這種尺度因子確定方法同單一尺度因子法等其他3種方法一樣，仍然依賴于水文模型的精度，只是這種依賴性相對(duì)于Klees方法要小。由于不同的水文模型所使用的數(shù)據(jù)源不同，其中包含的信息、精度等也存在差異，同時(shí)測(cè)站在全球范圍內(nèi)的分布并不均勻，測(cè)站分布較少的區(qū)域其水文模型精度也會(huì)受到一定影響。在實(shí)際應(yīng)用時(shí)，應(yīng)當(dāng)同時(shí)兼顧研究區(qū)域內(nèi)水文模型的精度和陸地水的時(shí)間變化規(guī)律等特征，來確定尺度因子的計(jì)算方法。例如文獻(xiàn)[14]在研究亞馬遜流域水儲(chǔ)量變化時(shí)，由于區(qū)域內(nèi)水文模型精度較低，而Klees方法又強(qiáng)烈依賴水文模型精度[13]，所以使用單一尺度因子法即可以取得比較理想的結(jié)果。文獻(xiàn)[12]在研究印度北部的水儲(chǔ)量變化時(shí)，發(fā)現(xiàn)研究區(qū)域內(nèi)水儲(chǔ)量變化的季節(jié)變化影響比較大，若使用單一尺度因子則反演結(jié)果不能體現(xiàn)這種季節(jié)變化趨勢(shì)，所以將長(zhǎng)周期項(xiàng)和季節(jié)性周期項(xiàng)分開，分別計(jì)算尺度因子。另外，由于水文模型通常在兩極地區(qū)存在數(shù)據(jù)缺失，在大洋上也沒有數(shù)據(jù)，所以在利用GRACE研究?jī)蓸O冰蓋融化、海平面變化時(shí)，這類利用水文模式數(shù)據(jù)求取尺度因子的方法并不適用。

近來，也有學(xué)者結(jié)合水井實(shí)測(cè)數(shù)據(jù)研究了華北地區(qū)的地下水儲(chǔ)量變化情況，文獻(xiàn)[15]利用京津冀平原地區(qū)的水井觀測(cè)數(shù)據(jù)得出地下水大致以-6 mm/a 的速度下降。水井觀測(cè)數(shù)據(jù)主要反映的是淺層地下水的變化，對(duì)于深層承壓水的變化并不能很好地探測(cè)，華北地區(qū)深層承壓水的超采現(xiàn)象十分嚴(yán)重。而GRACE反映的是各含水層總體的變化，所以水井觀測(cè)的地下水變化量遠(yuǎn)小于GRACE探測(cè)的地下水變化[15]。由于文獻(xiàn)[15]所使用的水井?dāng)?shù)據(jù)分布并不均勻，而華北地區(qū)地下水儲(chǔ)量減少主要發(fā)生在太行山山前平原地區(qū)，本文研究區(qū)域內(nèi)的燕山山脈和山東丘陵均缺少測(cè)站分布，這些地區(qū)的地下水開采量較小，因此在分析區(qū)域內(nèi)總體的地下水儲(chǔ)量變化時(shí)，其結(jié)果較利用京津冀平原地區(qū)水井觀測(cè)數(shù)據(jù)得到的地下水變化結(jié)果要小。如果可以使用在研究區(qū)域均勻分布的水井觀測(cè)數(shù)據(jù)，則可以進(jìn)一步驗(yàn)證GRACE反演結(jié)果。但本文地下水的反演結(jié)果與文獻(xiàn)[5](-4.1 mm/a，2003—2012年)和文獻(xiàn)[7](-5 mm/a，2002—2010年)的結(jié)果大致相同。存在差異的原因一方面是本文使用了顧及季節(jié)影響的尺度因子方法來減小GRACE數(shù)據(jù)后處理誤差，結(jié)果更具有可靠性，另外一方面是因?yàn)檠芯康膮^(qū)域和時(shí)段存在一定的差異。

5 結(jié) 論

本文提出了一種顧及季節(jié)的影響尺度因子計(jì)算方法，用于降低GRACE后處理影響，并利用GRACE RL05數(shù)據(jù)、CPC水文數(shù)據(jù)和GPCC降水?dāng)?shù)據(jù)，研究了華北平原的陸地水變化，得到如下結(jié)論：

(1) 提出的基于移去-恢復(fù)方法計(jì)算求得顧及季節(jié)影響的尺度因子對(duì)比不使用尺度因子和使用單一的尺度因子恢復(fù)GRACE后處理影響時(shí)，效果分別提升了48%和30%，并通過試驗(yàn)證明了這種方法對(duì)水文模型的依賴性不強(qiáng)。對(duì)多尺度因子恢復(fù)濾波后時(shí)間序列進(jìn)行連續(xù)小波變換，發(fā)現(xiàn)該方法并不改變信號(hào)的能量分布和大小，只在高頻部分有微弱減小。

(2) GRACE和CPC均表現(xiàn)出明顯的季節(jié)性特征，但GRACE的波動(dòng)更大，二者在振幅上也存在差異。2008年之前GRACE與CPC符合較好，華北平原的陸地水和地表水分別以(-7.9±2.4)mm/a、(-7.3±2.8)mm/a的速度下降，但2008年之后，GRACE與CPC的差異逐步增加，陸地水和地表水分別以(4.4±1.3)mm/a、(10.9±2.1)mm/a的速度上升。在整個(gè)研究期間內(nèi)，華北平原的陸地水和地表水整體分別呈現(xiàn)(-2.0±0.6)mm/a、(2.9±0.7)mm/a的變化趨勢(shì)。

(3) 華北平原地下水在2008年之前基本保持不變，僅以(-0.6±1.4)mm/a的速度減少，2008年之后地下水平均以(-6.5±1.6)mm/a的速度下降，下降速度以0.9 mm/a2的趨勢(shì)減小，地下水的振幅以(12.9±1.2)mm/a的速度增大，說明超采現(xiàn)象有加重趨勢(shì)，但降水的補(bǔ)充作用提升。研究期間內(nèi)地下水整體呈現(xiàn)(-4.8±0.7)mm/a的下降趨勢(shì)。

(4) GPCC降水量變化趨勢(shì)與GRACE和CPC基本一致，但其變化幅度略大于GRACE和CPC，說明降水是華北平原陸地水量變化的重要決定因素之一。利用交叉小波譜分析了GRACE和GPCC時(shí)間序列，發(fā)現(xiàn)二者均表現(xiàn)周年際的季節(jié)性變化特征，2008年之前GRACE約滯后GPCC序列2.6個(gè)月。季節(jié)性變化特征在2009年前后減弱，主要是因?yàn)榈叵滤_采量驟增，破壞了該地區(qū)的水平衡，2011年之后這種特征恢復(fù)，但GRACE滯后時(shí)間拉長(zhǎng)為4.6個(gè)月，這也說明了人類活動(dòng)是陸地水量變化的另一重要決定因素。