基于GRACE RL06 數據探測華北平原水儲量變化

何國慶

（1. 上海市測繪院，上海 200063）

GRACE 衛星由美國國家航空航天局（NASA）和德國航空中心（DLR）共同研制、部署，其科學意義是探測地球時變重力場信號分析地球質量密度變化與地球重力場的時變特征[1-9]。本文以華北平原為研究對象，利用最新的GRACE RL06重力衛星數據反演華北平原近年來水儲量變化，分析其時空變化特征，旨在實現華北平原水儲量時空連續監測，為華北平原合理開發利用水資源提供依據。

1 數據與方法

1.1 GRACE重力數據

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

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

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

1.2 Mascon數據

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

1.3 水文數據

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

1.4 基于尺度因子的恢復

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

2 結果分析

2.1 數據處理

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

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

為了分析華北平原水儲量變化在時間尺度上的特征，通過最小二乘擬合對華北平原水儲量時間序列進行以下形式擬合：

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

2.2 GRACE反演結果

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

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

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

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

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

圖3 華北平原2002-08～2016-08水儲量變化周年與半周年項

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

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

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

圖4 華北平原水儲量各月空間分布

2.3 GRACE與GLDAS和Mascon結果對比

本文采用GLDAS 水文模型中的Noah 數據和CSR提供的Mascon 數據與GRACE 觀測結果進行比較分析。前文已經介紹了GLDAS可提供降雨、降雪、蒸發量、大氣壓、風速等多種數據，本文將選取其提供的地表積雪數據、土壤水含量及植被冠層水含量數據來模擬研究區域的水儲量變化，其中土壤水含量分為4 層，包含了地下0～10 cm、10～40 cm、40～100 cm 及100～200 cm的土壤水。

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

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

2.4 與GPCC降水數據對比

區域水儲量的變化是由區域內降水、蒸發、徑流共同作用的，其水平衡方程表達式可寫作為：

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

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

3 結 語

本文利用2002-08～2016-08 的GRACE RL06 版本GSM 數據反演了華北平原水儲量變化，同時結合了GLDAS水文模型數據、Mascon數據及GPCC降水數據對華北平原區域水儲量時空變化進行了詳細的分析，得到如下結論。

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

2）GRACE 反演結果與GLDAS 和Mascon 計算結果具有較強的一致性，GRACE 觀測的季節性變化振幅高于水文模型模擬的結果，與Mascon計算結果相關度最高。三者都能探測到華北平原在研究期間內發生了干旱和洪澇災害。

3）GRACE 反演水儲量流量變化與GPCC 降水數據具有較好的相關性，2016年華北平原水儲量的增加和2002年華北平原的旱情，都在同時期降水數據上有著良好的反映，說明降雨是華北平原水儲量變化的重要因素之一。