何國(guó)慶

（1. 上海市測(cè)繪院，上海 200063）

GRACE 衛(wèi)星由美國(guó)國(guó)家航空航天局（NASA）和德國(guó)航空中心（DLR）共同研制、部署，其科學(xué)意義是探測(cè)地球時(shí)變重力場(chǎng)信號(hào)分析地球質(zhì)量密度變化與地球重力場(chǎng)的時(shí)變特征[1-9]。本文以華北平原為研究對(duì)象，利用最新的GRACE RL06重力衛(wèi)星數(shù)據(jù)反演華北平原近年來(lái)水儲(chǔ)量變化，分析其時(shí)空變化特征，旨在實(shí)現(xiàn)華北平原水儲(chǔ)量時(shí)空連續(xù)監(jiān)測(cè)，為華北平原合理開發(fā)利用水資源提供依據(jù)。

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

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

GRACE之所以能反演水儲(chǔ)量變化，是基于John[10]的理論依據(jù)，利用GRACE 時(shí)變重力場(chǎng)模型數(shù)據(jù)反演水儲(chǔ)量變化時(shí)，通常以等效水高的形式表示：

式中，θ為地心余緯；λ為地心經(jīng)度；a為地球半徑；ρa(bǔ)ve為地球平均密度（ρa(bǔ)ve=5 517 kg/m3）；kl為l階負(fù)荷勒夫數(shù)；ρw為水的密度（ρw=1 000 kg/m3）；l、m為球諧系數(shù)的階數(shù)和次數(shù)；Pˉlm(cosθ)為完全規(guī)格化的締合勒讓德函數(shù)；ΔClm、ΔSlm為完全規(guī)格化的球諧系數(shù)變化量；Wlm為平滑函數(shù)。

本文采用GSR公布的自2002-08～2016-08共154個(gè)月的RL06 版本GSM 數(shù)據(jù)，其中存在部分月份數(shù)據(jù)丟失的情況。RL06 數(shù)據(jù)相比于RL05 采用了新的背景模型，改善了數(shù)據(jù)處理方法，使得數(shù)據(jù)精度明顯提高。在利用式（1）進(jìn)行計(jì)算時(shí)需要注意：受衛(wèi)星軌道及星載儀器等影響，GRACE球諧系數(shù)中低階項(xiàng)（一階項(xiàng)和二階項(xiàng)）觀測(cè)精度較差，需要對(duì)其進(jìn)行替換；GRACE時(shí)變重力場(chǎng)數(shù)據(jù)中高階項(xiàng)存在著較大的誤差，并且誤差隨著階數(shù)的升高而增大，對(duì)于該項(xiàng)誤差常用高斯濾波或各向異性濾波來(lái)處理；GRACE球諧系數(shù)奇（偶）數(shù)階項(xiàng)數(shù)之間存在一定的相關(guān)性，這種相關(guān)性導(dǎo)致了反演結(jié)果在圖上表現(xiàn)為“南北條帶”誤差，針對(duì)該誤差處理方法，通常采用去相關(guān)濾波進(jìn)行處理。

1.2 Mascon數(shù)據(jù)

Mascon 數(shù)據(jù)由CSR、GFZ、JPL 三家機(jī)構(gòu)利用GRACE Level-1 數(shù)據(jù)各自解算的地球重力場(chǎng)模型數(shù)據(jù)。Watkins[11]等國(guó)外學(xué)者通過(guò)實(shí)驗(yàn)發(fā)現(xiàn)Mascon 數(shù)據(jù)較球諧系數(shù)反演的重力場(chǎng)模型結(jié)果有更高的時(shí)空分辨率和信噪比。Mascon數(shù)據(jù)是重力場(chǎng)模型的另一種表達(dá)方式，數(shù)據(jù)已經(jīng)替換過(guò)了低階項(xiàng)、改正了GIA（冰川均衡調(diào)整）等步驟，不用進(jìn)行后處理和空間濾波平滑。與球諧系數(shù)不同的是，球諧系數(shù)沒有特定的位置信息，而Mascon數(shù)據(jù)中每個(gè)質(zhì)量集中塊都具有特定的地球物理位置。

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

收集了研究期間內(nèi)的水文模型數(shù)據(jù)和降水?dāng)?shù)據(jù)。其中水文模型采用GLDAS 數(shù)據(jù)，GLDAS 由美國(guó)國(guó)家航空航天局NASA 和美國(guó)國(guó)家環(huán)境預(yù)報(bào)中心NCEP 共同建立，作為高精度的水文模型，能通過(guò)衛(wèi)星與地面觀測(cè)數(shù)據(jù)約束模型，提供最優(yōu)化的近實(shí)時(shí)的地球表面質(zhì)量變化信息[12]。降水?dāng)?shù)據(jù)來(lái)自GPCC，它是由德國(guó)氣象局DWD 成立的，根據(jù)地面觀測(cè)數(shù)據(jù)整合了全球陸地表面降水?dāng)?shù)據(jù)集，并且還接受了世界氣象組織和全球電信傳輸系統(tǒng)的每天地面天氣和月氣候數(shù)據(jù)[13]。GLDAS 提供的水文模型數(shù)據(jù)有Noah、VIC、Mosaic和CLM；其時(shí)間分辨率有3 h和1個(gè)月，空間分辨率提供0.25°×0.25°和1°×1°；GPCC數(shù)據(jù)提供了4種空間分辨率：0.25°×0.25°、0.5°×0.5°、1°×1°、2.5°×2.5°。本文均采用1°×1°、時(shí)間分辨率為1個(gè)月的水文數(shù)據(jù)，時(shí)間跨度選取與GRACE數(shù)據(jù)相同。

1.4 基于尺度因子的恢復(fù)

由于GRACE 重力模型數(shù)據(jù)提供的球諧系數(shù)是有限的，無(wú)法完全展開，這就造成了部分信號(hào)無(wú)法反映的截?cái)嗾`差；另一方面，使用濾波方法對(duì)球諧系數(shù)進(jìn)行處理時(shí)，也會(huì)產(chǎn)生誤差，使得GRACE 估計(jì)的結(jié)果是有偏差的。本文將采用尺度因子方法來(lái)恢復(fù)上述內(nèi)容中“泄露”的信號(hào)，即利用水文模型數(shù)據(jù)來(lái)模擬泄露誤差的影響，再利用尺度因子恢復(fù)GRACE 反演信號(hào)。具體操作步驟是：①對(duì)水文模型數(shù)據(jù)進(jìn)行球諧展開，并展開到和GRACE 數(shù)據(jù)一樣的階數(shù)；②對(duì)展開后的水文模型球諧系數(shù)進(jìn)行和GRACE數(shù)據(jù)同樣的后處理操作，得到研究區(qū)域的質(zhì)量變化時(shí)間序列ΔS′；③計(jì)算原始水文數(shù)據(jù)得到的研究區(qū)域質(zhì)量變化的時(shí)間序列ΔS；④將處理后的時(shí)間序列ΔS′和原始數(shù)據(jù)的時(shí)間序列ΔS進(jìn)行最小二乘擬合，計(jì)算出尺度因子k；⑤對(duì)GRACE反演的研究區(qū)域質(zhì)量變化時(shí)間序列信號(hào)乘上尺度因子k。其中尺度因子k滿足處理后的時(shí)間序列ΔS′要和原始數(shù)據(jù)得到的時(shí)間序列ΔS的殘差平方和最小。即：

2 結(jié)果分析

2.1 數(shù)據(jù)處理

本文采用CSR提供的GRACE RL06 GSM重力模型數(shù)據(jù)產(chǎn)品，截取時(shí)間跨度為2002-08～2016-08 共計(jì)154 個(gè)月份。對(duì)球諧系數(shù)中一階項(xiàng)利用Chambers 解算出的一階項(xiàng)數(shù)據(jù)進(jìn)行替代，對(duì)球諧系數(shù)中二階項(xiàng)數(shù)據(jù)利用SLR 提供的數(shù)據(jù)進(jìn)行替代。對(duì)于南北條帶誤差，采用Duan 滑動(dòng)窗口去相關(guān)方法進(jìn)行去條帶處理。采用平滑半徑為300 km 的扇形濾波方法抑制高階項(xiàng)噪聲。最后，對(duì)得到的全球水儲(chǔ)量變化格網(wǎng)數(shù)據(jù)進(jìn)行按緯度的余弦值定權(quán)，并根據(jù)研究區(qū)域的選擇，對(duì)研究區(qū)域的格網(wǎng)數(shù)據(jù)進(jìn)行加權(quán)平均計(jì)算出該研究區(qū)域的平均水儲(chǔ)量變化的時(shí)間序列。

利用式（2）計(jì)算出在Duan 滑動(dòng)窗口去相關(guān)濾波和300 km扇形濾波方法處理下的尺度因子k，這里計(jì)算出的尺度因子k=1.573，本文利用此值乘以GRACE計(jì)算出的結(jié)果，作為GRACE反演水儲(chǔ)量的最終估計(jì)值。

為了分析華北平原水儲(chǔ)量變化在時(shí)間尺度上的特征，通過(guò)最小二乘擬合對(duì)華北平原水儲(chǔ)量時(shí)間序列進(jìn)行以下形式擬合：

式中，b為趨勢(shì)項(xiàng)；Ai為振幅；φi為相位；Ti為周期（T1為周年，T2為半周年）。通過(guò)以上擬合公式可以求出華北平原水儲(chǔ)量變化的年際變化趨勢(shì)、周年和半周年項(xiàng)等變化特征。

2.2 GRACE反演結(jié)果

華北平原是我國(guó)第二大平原，坐落于我國(guó)東部，西到太行山，北達(dá)燕山，東鄰渤海，南至黃河，是我國(guó)主要的經(jīng)濟(jì)、文化、政治中心。利用GRACE時(shí)變重力模型數(shù)據(jù)計(jì)算出的研究期間2002-08～2016-08 全球及華北平原水儲(chǔ)量變化速率，從圖1可以看出華北平原在研究期間內(nèi)水儲(chǔ)量整體上呈遞減趨勢(shì)。在空間上，華北平原西南部遞減的速率最大，東北部遞減速率相對(duì)較小，遞減速率呈北緩南快的空間分布特征。

圖1 2002年8月至2016年8月華北平原水儲(chǔ)量變化速率分布

圖2 給出了由GRACE 反演的華北平原2002-08～2016-08 水儲(chǔ)量變化的時(shí)間序列圖，以便對(duì)華北平原年際變化特征進(jìn)行具體分析。如圖所示，在2002-08～2016-08 期間，華北平原水儲(chǔ)量變化總體趨勢(shì)上呈現(xiàn)下降現(xiàn)象，下降速率為1.38 cm/a。2002年，我國(guó)北方遭遇旱災(zāi)，黃河遭遇特枯、京杭大運(yùn)河山東濟(jì)寧段斷航50 d、河北遼河干流斷流158 d[14]，受2002年極端干旱氣候的影響,華北平原水儲(chǔ)量在2002—2003 年處于較低值。2004年我國(guó)雨水豐富，部分中小河流發(fā)生較大洪水，所以在2002—2004年期間，華北水儲(chǔ)量變?yōu)樯仙厔?shì)，上升速率達(dá)到11.13 cm/a，2004 年華北平原水儲(chǔ)量到達(dá)研究期間內(nèi)的峰值；2004—2008 年間，華北平原水儲(chǔ)量以2.87 cm/a的速率下降；2009—2013年，華北平原水儲(chǔ)量表現(xiàn)為緩慢增長(zhǎng)的趨勢(shì)，增長(zhǎng)速率為0.28 cm/a；2013—2016年，華北平原水儲(chǔ)量減少速率最快為4.29 cm/a；2016 年7 月我國(guó)發(fā)生大面積的洪澇災(zāi)害，華北平原水儲(chǔ)量在該期間內(nèi)增幅較大。

圖2 華北平原2002-08～2016-08水儲(chǔ)量變化

由圖2 可以看出，華北平原水儲(chǔ)量變化存在明顯的季節(jié)性變化，分別在冬季達(dá)到最大值，在夏季達(dá)到最小值。了解周年和半周年變化振幅和相位有利于了解華北平原水儲(chǔ)量季節(jié)性變化特征，圖3 計(jì)算出了由GRACE 解算出的2002 年8 月至2016 年8 月期間華北平原水儲(chǔ)量變化的周年、半周年項(xiàng)。

圖3 華北平原2002-08～2016-08水儲(chǔ)量變化周年與半周年項(xiàng)

圖3 表明華北平原南部屬于黃河下游區(qū)域，雨水量較充分，因此華北平原南部周年和半周年振幅大于北方區(qū)域；周年振幅由華北平原西南向東北逐漸減小，最大值出現(xiàn)在黃河下游河南境內(nèi)達(dá)到31.5 mm；華北平原周年相位表現(xiàn)為東部略大于西部，整個(gè)華北平原平均周年相位大約為325 d，說(shuō)明華北平原水儲(chǔ)量峰值一般出現(xiàn)在10～11 月附近，這一結(jié)論與圖2 華北平原水儲(chǔ)量變化時(shí)間序列所給出的信息相一致；與華北平原水儲(chǔ)量變化的周年振幅相比，半周年振幅相對(duì)較小，其空間分布特征與周年振幅分布特征相似，也是由西南向東北方向遞減，半周年振幅最大值達(dá)16.1 mm；華北平原半周年相位在空間上呈現(xiàn)西部大于東部的現(xiàn)象，說(shuō)明華北平原半周年變化振幅東部快于西部，整個(gè)華北平原半周年相位平均約為130 d。

為了更好地揭示華北平原年內(nèi)水儲(chǔ)量變化的時(shí)空特征，本文利用GRACE 數(shù)據(jù)計(jì)算了華北平原2002-08～2016-08期間內(nèi)各個(gè)月份的平均水儲(chǔ)量，如圖4所示，給出了華北平原研究期間內(nèi)多年各月平均水儲(chǔ)量的空間分布。

從月尺度上分析，有利于剔除時(shí)間上的復(fù)雜性和空間上的異質(zhì)性，由圖4可知從1月份至6月份華北平原水儲(chǔ)量不斷減少，6 月份至9 月份水儲(chǔ)量不斷增加，即華北平原水儲(chǔ)量冬春之際處于不斷虧損狀態(tài)，夏季屬于盈余補(bǔ)充狀態(tài)。其中，1、2月份水儲(chǔ)量空間分布差異不大，都是東部水儲(chǔ)量高于西部，這是因?yàn)槿A北平原東鄰渤海，東部區(qū)域氣候相對(duì)濕潤(rùn)的原因。整個(gè)華北平原1 月份水儲(chǔ)量略高于2 月水儲(chǔ)量，水儲(chǔ)量變化范圍在9～17 mm 之間，整個(gè)華北平原3 月份水儲(chǔ)量總值大概為0 左右，表明了華北平原3 月份水儲(chǔ)量達(dá)到了1 a內(nèi)的均值。4～7月華北平原水儲(chǔ)量處于負(fù)值，是1 a 中最缺水的幾個(gè)月份，且河南境內(nèi)缺水表現(xiàn)最為嚴(yán)重，北京天津地區(qū)缺水情況相對(duì)較輕，其中4～6 月華北平原水儲(chǔ)量空間分布呈現(xiàn)為西南地區(qū)低于東北地區(qū)，六月份整個(gè)華北平原水儲(chǔ)量達(dá)到1 a當(dāng)中最小值，其均值約為-50 mm；7～9月份華北平原水儲(chǔ)量不斷增加，由虧損狀態(tài)轉(zhuǎn)為盈余狀態(tài)，9月份水儲(chǔ)量達(dá)到了1 a之中最大值，且空間分布上表現(xiàn)為華北平原西南地域盈余最大達(dá)到34 mm，東北部相對(duì)較小，最小地區(qū)只有16 mm；10～11月份華北平原水儲(chǔ)量空間分布情況與8月份恰好相反，西南地區(qū)水儲(chǔ)量高于東北部?？傮w來(lái)說(shuō)，華北平原春夏水儲(chǔ)量較少，秋冬水儲(chǔ)量較多；且冬春水儲(chǔ)量整體不斷減少，在6 月份達(dá)到谷值，夏季水儲(chǔ)量不斷補(bǔ)充在9月份達(dá)到峰值。

圖4 華北平原水儲(chǔ)量各月空間分布

2.3 GRACE與GLDAS和Mascon結(jié)果對(duì)比

本文采用GLDAS 水文模型中的Noah 數(shù)據(jù)和CSR提供的Mascon 數(shù)據(jù)與GRACE 觀測(cè)結(jié)果進(jìn)行比較分析。前文已經(jīng)介紹了GLDAS可提供降雨、降雪、蒸發(fā)量、大氣壓、風(fēng)速等多種數(shù)據(jù)，本文將選取其提供的地表積雪數(shù)據(jù)、土壤水含量及植被冠層水含量數(shù)據(jù)來(lái)模擬研究區(qū)域的水儲(chǔ)量變化，其中土壤水含量分為4 層，包含了地下0～10 cm、10～40 cm、40～100 cm 及100～200 cm的土壤水。

圖5給出了由GRACE、GLDAS Noah和Mascon 3種數(shù)據(jù)得到的華北平原水儲(chǔ)量變化時(shí)間序列。由圖可知，通過(guò)3 種數(shù)據(jù)計(jì)算得到的華北平原水儲(chǔ)量長(zhǎng)期變化趨勢(shì)基本一致。經(jīng)計(jì)算，GRACE 與GLDAS 反演結(jié)果相關(guān)系數(shù)為0.76；Mascon與GLDAS結(jié)果相關(guān)系數(shù)為0.71；GRACE 與Mascon 相關(guān)性最高為0.94，3 種數(shù)據(jù)反演結(jié)果呈強(qiáng)相關(guān)。從反演結(jié)果振幅上來(lái)看，GRACE反演的結(jié)果振幅大于Mascon 估算的結(jié)果，而Mascon估算的結(jié)果振幅又大于GLDAS水文模型的結(jié)果，這是因?yàn)镚RACE 反演的結(jié)果相較于Mascon 進(jìn)行了信號(hào)恢復(fù)，GRACE和Mascon反演結(jié)果又相較于GLDAS水文模型不僅包含了土壤表層水、地表積雪和植被冠層水，同時(shí)還包括了深層地下水和各種形式存在的水元素。與GRACE 和Mascon 結(jié)果相比，GLDAS 水文模型估算的結(jié)果在2011 年以后高于GRACE 和Mascon 結(jié)果，原因是近年來(lái)華北平原地下水遭到了嚴(yán)重開采，農(nóng)業(yè)灌溉和生活用水主要來(lái)源于地下水，資料顯示[15]2012年華北平原地下水開采量已達(dá)當(dāng)年水資源利用量的6成以上，地下水又是GRACE反演結(jié)果中主要的組成部分，所以地下水的超采驟減會(huì)相應(yīng)導(dǎo)致GRACE反演結(jié)果相較于平均水平驟縮，對(duì)GLDAS估算結(jié)果沒有直接影響。GRACE、Mascon 和GLDAS 反演結(jié)果都顯示出2004 年和2016 年水儲(chǔ)量增幅最大，2002 年和2009年水儲(chǔ)量變化偏小，這說(shuō)明華北平原在2002年、2004 年、2009 年及2016 年發(fā)生了水文事件。查閱中國(guó)水旱災(zāi)害公報(bào)發(fā)現(xiàn)華北平原在2002 年和2009 年出現(xiàn)嚴(yán)重的旱災(zāi)；2004年和2016年7月全國(guó)出現(xiàn)大面積的洪澇災(zāi)害，華北平原受災(zāi)嚴(yán)重。

圖5 GRACE、GLDAS和Mascon估算的華北平原水儲(chǔ)量變化

2.4 與GPCC降水?dāng)?shù)據(jù)對(duì)比

區(qū)域水儲(chǔ)量的變化是由區(qū)域內(nèi)降水、蒸發(fā)、徑流共同作用的，其水平衡方程表達(dá)式可寫作為：

式中，ΔS為區(qū)域水儲(chǔ)量變化量；P為降水量；E為蒸散發(fā)；R為地表、地下徑流。我們可以看出降水量是影響區(qū)域水儲(chǔ)量變化中的唯一輸入量，因此降水量對(duì)區(qū)域水儲(chǔ)量的影響至關(guān)重要。

因?yàn)镚RACE 觀測(cè)的某月區(qū)域水儲(chǔ)量也包括了研究區(qū)域內(nèi)該月份之前的降水部分，因此GRACE 觀測(cè)結(jié)果會(huì)較降水結(jié)果出現(xiàn)一定的滯后現(xiàn)象，為了解決這一問(wèn)題，使GRACE 觀測(cè)結(jié)果與降水結(jié)果更具有直接相關(guān)性，本文對(duì)GRACE 得到的華北平原水儲(chǔ)量變化時(shí)間序列做進(jìn)一步處理：首先針對(duì)GRACE 中缺失的月份數(shù)據(jù)利用三次樣條插值進(jìn)行補(bǔ)充，然后將2002-08～2016-08 水儲(chǔ)量變化數(shù)據(jù)用該月份數(shù)據(jù)減去前一個(gè)月份的數(shù)據(jù)，這樣就得到了該月的水儲(chǔ)量變化量。通過(guò)上述步驟計(jì)算得到的華北平原各個(gè)月份水儲(chǔ)量變化量去和降水?dāng)?shù)據(jù)進(jìn)行比較發(fā)現(xiàn)，由GRACE 計(jì)算出的月水儲(chǔ)量變化量與GPCC月降水量峰值和谷值對(duì)應(yīng)較好，峰值出現(xiàn)在1 a 的夏秋季；谷值出現(xiàn)在1 a的冬春季，說(shuō)明華北平原水儲(chǔ)量流量的季節(jié)性變化與降水的季節(jié)性規(guī)律有較好的一致性；水儲(chǔ)量流量的變化走勢(shì)與降水量的走勢(shì)較為吻合，說(shuō)明降水量是影響華北平原水儲(chǔ)量變化的重要因素，但水儲(chǔ)量的變化還受蒸散發(fā)和徑流的影響，所以兩者走勢(shì)并不會(huì)完全相同。

3 結(jié) 語(yǔ)

本文利用2002-08～2016-08 的GRACE RL06 版本GSM 數(shù)據(jù)反演了華北平原水儲(chǔ)量變化，同時(shí)結(jié)合了GLDAS水文模型數(shù)據(jù)、Mascon數(shù)據(jù)及GPCC降水?dāng)?shù)據(jù)對(duì)華北平原區(qū)域水儲(chǔ)量時(shí)空變化進(jìn)行了詳細(xì)的分析，得到如下結(jié)論。

1）華北平原水儲(chǔ)量變化呈現(xiàn)出明顯的季節(jié)性和周年性變化，春夏水少、秋冬水多，水儲(chǔ)量在1 a中6月份達(dá)到最小值，9月份達(dá)到最大值。研究期間內(nèi)華北平原水儲(chǔ)量以1.38 cm/a 的速度減少，其中2004—2008 年間，以2.87 cm/a 的速率減少；2009—2013 年，以0.28 cm/a的速率緩慢增長(zhǎng)；2013—2016年，下降速率達(dá)4.29 cm/a。空間分布整體表現(xiàn)為由東北部向西南部減少速率不斷加快。

2）GRACE 反演結(jié)果與GLDAS 和Mascon 計(jì)算結(jié)果具有較強(qiáng)的一致性，GRACE 觀測(cè)的季節(jié)性變化振幅高于水文模型模擬的結(jié)果，與Mascon計(jì)算結(jié)果相關(guān)度最高。三者都能探測(cè)到華北平原在研究期間內(nèi)發(fā)生了干旱和洪澇災(zāi)害。

3）GRACE 反演水儲(chǔ)量流量變化與GPCC 降水?dāng)?shù)據(jù)具有較好的相關(guān)性，2016年華北平原水儲(chǔ)量的增加和2002年華北平原的旱情，都在同時(shí)期降水?dāng)?shù)據(jù)上有著良好的反映，說(shuō)明降雨是華北平原水儲(chǔ)量變化的重要因素之一。