楊帥， 鄭偉,3*， 尹文杰， 劉杰

1 河南理工大學(xué)測(cè)繪與國(guó)土信息工程學(xué)院， 焦作河南 454003 2 中國(guó)空間技術(shù)研究院錢學(xué)森空間技術(shù)實(shí)驗(yàn)室， 北京 100094 3 遼寧工程技術(shù)大學(xué)測(cè)繪與地理科學(xué)學(xué)院， 阜新遼寧 123000

0 引言

自2002年重力恢復(fù)與氣候?qū)嶒?yàn)(Gravity Recovery and Climate Experiment，GRACE)衛(wèi)星發(fā)射以來(lái)，GRACE重力衛(wèi)星提供了高精度觀測(cè)地球時(shí)變重力場(chǎng)信號(hào)的新手段(Tapley et al., 2004).截止到目前，GRACE重力衛(wèi)星已被廣泛應(yīng)用于陸地水儲(chǔ)量(Terrestrial Water Storage，TWS)變化(Rodell et al., 2009; Seyoum and Milewski, 2016)、冰川變化及其質(zhì)量平衡(Velicogna and Wahr, 2006a,b; Mohajerani et al., 2018)、全球平均海洋質(zhì)量變化(Chen et al., 2019)、地震(Han et al., 2006; Chen et al., 2007; Zhou et al., 2018)等多領(lǐng)域研究.GRACE重力衛(wèi)星在2017年10月停止運(yùn)行，隨著GRACE Follow-on重力衛(wèi)星的成功發(fā)射，該衛(wèi)星反演的地球重力場(chǎng)模型繼續(xù)在地球科學(xué)領(lǐng)域發(fā)揮著極其重要的作用(鄭偉等, 2014).

基于Wahr等(1998)的理論，地球表面的質(zhì)量變化可以利用GRACE球諧系數(shù)解反演得到.然而，球諧系數(shù)中既含有隨系數(shù)階數(shù)而增加的隨機(jī)誤差，又包含與系統(tǒng)相關(guān)的系統(tǒng)誤差.因此，許多研究者相繼提出了空間平滑濾波(如高斯平滑濾波(Wahr et al., 1998)、維納濾波(Sasgen et al., 2007)等)和去條帶濾波(Swenson和Wahr(2006)的去條帶方法、Chen等(2009)的P4M6方法等).這些濾波方法在消除誤差影響的同時(shí)也會(huì)改變真實(shí)的物理信號(hào).Chen等(2015)通過(guò)實(shí)驗(yàn)?zāi)M得出截?cái)嗉捌交?00 km產(chǎn)生的泄露誤差在南極洲可達(dá)80.95 Gt·a-1(實(shí)測(cè)值-150.00 Gt·a-1，模擬值-69.05 Gt·a-1)；Jin和Zou(2015)研究表明格陵蘭島冰蓋質(zhì)量損失的信號(hào)中，42.4%是由泄漏誤差造成，對(duì)泄漏誤差進(jìn)行校正后，格陵蘭島的冰川融化速度從-15.58±0.93 Gt·a-2增至-26.19±1.67 Gt·a-2.因此，泄露誤差是GRACE反演質(zhì)量變化中的主要限制因素，校正泄露誤差對(duì)提高反演結(jié)果的準(zhǔn)確性至關(guān)重要.

目前，GRACE數(shù)據(jù)泄露誤差常用的恢復(fù)方法包括：(1)尺度因子法(Scaling Factor，SF).Long等(2015a)利用全球水文模型PCR-GLOBWB得到的尺度因子用于恢復(fù)長(zhǎng)江流域水儲(chǔ)量變化的信號(hào)損失，該方法對(duì)于提高研究區(qū)水儲(chǔ)量變化的準(zhǔn)確性具有重要價(jià)值；Huang等(2019)基于全球不同水文模型對(duì)尺度因子泄露誤差校正法進(jìn)行敏感性分析，對(duì)中國(guó)西南部水儲(chǔ)量變化進(jìn)行了評(píng)估.尺度因子法依賴先驗(yàn)水文模型，雖然水文模型已被廣泛應(yīng)用于GRACE數(shù)據(jù)泄露誤差校正的先驗(yàn)信息(Long et al., 2014; Chen et al., 2017; Pan et al., 2017)，但是水文模型存在一些缺點(diǎn)，例如，地表水和地下水表示成分不完整以及沒有考慮人類活動(dòng)影響等.(2)迭代恢復(fù)法(Forward Modeling，F(xiàn)M).該方法不需要借助先驗(yàn)?zāi)Ｐ停饕蹽RACE原始觀測(cè)值的約束.Chen等(2014)利用迭代恢復(fù)法恢復(fù)了印度西北長(zhǎng)期地下水變化，新估算的水儲(chǔ)量變化趨勢(shì)反映迭代恢復(fù)法可以有效減少GRACE估計(jì)中的泄露誤差；吳云龍等(2015)基于迭代恢復(fù)法可有效恢復(fù)黑河流域水儲(chǔ)量變化初始信號(hào)，泄露誤差減少，提高水儲(chǔ)量變化的空間分辨效果.迭代恢復(fù)法計(jì)算過(guò)程相對(duì)復(fù)雜，在部分區(qū)域可能存在信號(hào)恢復(fù)不全或過(guò)量恢復(fù)現(xiàn)象.與尺度因子法相比，迭代恢復(fù)法是一種完全不同于尺度因子法的泄露誤差恢復(fù)方法，恢復(fù)結(jié)果可與尺度因子法校正結(jié)果進(jìn)行驗(yàn)證對(duì)比分析.

不同于前人的已用研究，本文首次利用模型與GRACE數(shù)據(jù)的均方根誤差和相關(guān)性加權(quán)尺度因子，構(gòu)建新型尺度因子校正法(New Scale Factor Correction Method，NSFCM).該方法充分考慮多種水文模型與GRACE之間的一致性和準(zhǔn)確性信息，相比于傳統(tǒng)尺度因子法和迭代恢復(fù)法，可以提高GRACE陸地水儲(chǔ)量變化的準(zhǔn)確性.

1 新型尺度因子校正法的構(gòu)建

傳統(tǒng)尺度因子法使用單一水文模型計(jì)算的尺度因子依賴水文模型結(jié)果.然而，不同的水文模型由于模型的結(jié)果以及驅(qū)動(dòng)數(shù)據(jù)的差異，導(dǎo)致模型輸出結(jié)果間也存在一定的差異(Long et al., 2013).因此，本文根據(jù)多種水文模型計(jì)算的陸地水儲(chǔ)量與CSR-SHc數(shù)據(jù)之間的一致性和準(zhǔn)確性信息賦予權(quán)重，綜合多種尺度因子模型提高GRACE陸地水儲(chǔ)量變化的準(zhǔn)確性.

新型尺度因子校正法主要分為以下幾個(gè)步驟：(1)水文格網(wǎng)數(shù)據(jù)轉(zhuǎn)化為球諧數(shù)據(jù)；(2)高斯平滑濾波和去條帶濾波；(3)球諧系數(shù)轉(zhuǎn)化為格網(wǎng)數(shù)據(jù)；(4)基于最小二乘原理計(jì)算尺度因子；(5)計(jì)算水文模型的陸地水儲(chǔ)量與CSR-SHc數(shù)據(jù)均方根誤差和相關(guān)性；(6)基于模型權(quán)重加權(quán)尺度因子.具體方法如下：

(1)水文格網(wǎng)數(shù)據(jù)轉(zhuǎn)化為球諧數(shù)據(jù)

根據(jù)Wahr等(1998)的理論，可將水文模型陸地水儲(chǔ)量轉(zhuǎn)為與之對(duì)應(yīng)的球諧系數(shù)：

×sinθdθdλ，

(1)

(2)高斯平滑濾波和去條帶濾波

為了保持與GRACE數(shù)據(jù)處理相同的濾波，本文使用Swenson去條帶濾波和300 km的高斯平滑濾波.去條帶濾波通常采用以l階為中心，寬為w的二次多項(xiàng)式濾波窗口，Swenson和Wahr(2006)采用的平滑窗口寬度w和球諧系數(shù)次數(shù)m的關(guān)系為

(2)

其中，max表示取兩者的最大值，A=30和K=10表示選擇的經(jīng)驗(yàn)參數(shù).

高斯濾波是指對(duì)不同階次的球諧系數(shù)賦予不同權(quán)重，從而減小GRACE數(shù)據(jù)高階項(xiàng)誤差的方法.Jekeli(1981)提出用遞歸關(guān)系計(jì)算權(quán)重：

(3)

其中，r表示高斯平滑半徑，a表示地球平均半徑，e表示自然常數(shù)，Wl表示權(quán)重.

(3)球諧系數(shù)轉(zhuǎn)化為格網(wǎng)數(shù)據(jù)

根據(jù)Wahr等(1998)的理論，可將濾波后的球諧系數(shù)轉(zhuǎn)為格網(wǎng)數(shù)據(jù)：

(4)

(4)基于最小二乘原理計(jì)算尺度因子

尺度因子是基于水文模型陸地水儲(chǔ)量模擬濾波過(guò)程造成的信號(hào)損失，并以最小二乘減小原始模型時(shí)間序列與濾波后模型時(shí)間序列之間差異的方法.對(duì)于任意格網(wǎng)點(diǎn)，在t時(shí)間序列對(duì)應(yīng)的未濾波的陸地水儲(chǔ)量為Δhi(t)，濾波后的陸地水儲(chǔ)量為Δηi(t)，通過(guò)最小二乘原則(Landerer and Swenson, 2012)：

(5)

其中，N表示時(shí)間序列個(gè)數(shù)，SF表示所求尺度因子.

(5)計(jì)算水文模型陸地水儲(chǔ)量與CSR-SHc數(shù)據(jù)的均方根誤差和相關(guān)性

皮爾遜相關(guān)系數(shù)(PR，Pearson correlation coefficient)通常用來(lái)描述兩個(gè)時(shí)間序列之間在一致性，均方根誤差(RMSE，Root Mean Square Error)反映了模擬值與“真實(shí)”值之間的偏差程度.對(duì)于任意格網(wǎng)點(diǎn)相關(guān)系數(shù)計(jì)算如下(Yin et al., 2020)：

(6)

對(duì)于任意格網(wǎng)點(diǎn)均方根誤差計(jì)算如下(Yin et al., 2020)：

(7)

則該格網(wǎng)點(diǎn)分配的權(quán)重為

W=PR+RMSE-1,

(8)

其中，W表示權(quán)重，N表示時(shí)間序列個(gè)數(shù)，Δh(t)表示t時(shí)間序列水文模型的陸地水儲(chǔ)量，ΔSH(t)表示t時(shí)間序列CSR-SHc的陸地水儲(chǔ)量.

(6)基于模型權(quán)重加權(quán)尺度因子

基于尺度因子模型的權(quán)重，可得新型尺度因子校正法的公式如下：

(9)

其中，NSFCM-SF表示新型尺度因子校正法計(jì)算的尺度因子，SF表示模型尺度因子，n表示尺度因子模型數(shù)目(CLM，MOS，NOAH，VIC，WGHM).

2 研究區(qū)域和數(shù)據(jù)

2.1 研究區(qū)域

長(zhǎng)江流域地處中國(guó)南部，經(jīng)緯度范圍為85°E—122°E和24°N—36°N，橫跨中國(guó)西部、中部和東部三大區(qū)域共計(jì)19個(gè)省、市和自治區(qū)，長(zhǎng)江流域位置及其海拔高度如圖1所示.長(zhǎng)江流域年平均降水量為1067 mm且時(shí)空分布不均勻，約70～90%集中在5—10月.長(zhǎng)江干流宜昌以上為長(zhǎng)江上游，流域面積為100×104km2，宜昌至湖口為長(zhǎng)江中游，流域面積為68×104km2，湖口至長(zhǎng)江入?？跒殚L(zhǎng)江下游，流域面積為12×104km2.鑒于長(zhǎng)江下游流域面積較小，將GRACE數(shù)據(jù)應(yīng)用于小于GRACE分辨率的區(qū)域，反演結(jié)果可能存在較大不確定性(Long et al., 2015b).因此，本研究將長(zhǎng)江中游和下游統(tǒng)稱為長(zhǎng)江中下游，如圖1虛線以東所示，虛線以西為長(zhǎng)江上游.

圖1 長(zhǎng)江流域位置及其海拔高度Fig.1 Location and elevation of the Yangtze River basin

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

本研究使用的數(shù)據(jù)包括兩種網(wǎng)格數(shù)據(jù)產(chǎn)品：CSR-SHc和CSR-M(CSR Mascon).CSR-SHc數(shù)據(jù)是基于美國(guó)德克薩斯大學(xué)空間研究中心(Center for Space Research，CSR)、美國(guó)噴氣推進(jìn)實(shí)驗(yàn)室(Jet Propulsion Laboratory，JPL)和德國(guó)波茲坦地球科學(xué)研究中心(GeoForschungsZentrum，GFZ)的RL05球諧系數(shù)得到，數(shù)據(jù)最大階次為60階，時(shí)間跨度范圍為2002年4月—2017年1月，利用三次樣條插值得到月份缺失的數(shù)據(jù)，每月網(wǎng)格數(shù)據(jù)表示相對(duì)于2004年1月—2009年12月平均值的質(zhì)量變化.GRACE的C20值比衛(wèi)星激光測(cè)距(Satellite Laser Ranging，SLR)的C20值具有更大的不確定性，因此，用SLR的C20系數(shù)解代替GRACE的C20系數(shù)(Cheng et al., 2011).使用Swenson等(2008)的方法估計(jì)一階項(xiàng)，應(yīng)用于地心校正.基于Geruo等(2012)的模型進(jìn)行冰川均衡調(diào)整(Glacial Isostatic Adjustment，GIA)校正.為了減少GRACE月度圖中南北條帶的相關(guān)誤差，對(duì)數(shù)據(jù)進(jìn)行了去條帶處理.此外，數(shù)據(jù)還應(yīng)用了300 km高斯平滑濾波以抑制球諧系數(shù)中的高階隨機(jī)誤差.

CSR-M數(shù)據(jù)是基于Level-1B數(shù)據(jù)解算的，沒有應(yīng)用任何平滑或去條帶濾波.因此，CSR-M數(shù)據(jù)無(wú)需進(jìn)一步處理即可使用.CSR-M數(shù)據(jù)空間分辨率為0.5°× 0.5°，該數(shù)據(jù)集并非針對(duì)特定應(yīng)用定制，適用于所有感興趣的科學(xué)領(lǐng)域，例如海洋學(xué)、陸地表面水文學(xué)、冰凍圈等研究.研究表明，CSR-M數(shù)據(jù)解適用于恢復(fù)全球和大流域尺度(大于約3°× 3°)中不同質(zhì)量源的質(zhì)量再分配，而球諧系數(shù)解在有效處理時(shí)，在較小或局部尺度上優(yōu)于CSR-M數(shù)據(jù)解(Zhang et al., 2019).

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

全球陸地同化系統(tǒng)水文模型(Global Land Data Assimilation System，GLDAS)，是由美國(guó)國(guó)家航空局和美國(guó)國(guó)家海洋與大氣管理局共同開發(fā)(Rodell et al., 2004).該模型通過(guò)陸地表面建模和數(shù)據(jù)同化技術(shù)，輸出陸地表面各項(xiàng)參數(shù)(土壤溫度和土壤水分、蒸發(fā)量和地表徑流等)和水文變量.本文使用的GLDAS數(shù)據(jù)包括四個(gè)陸地表面模型NOAH、CLM、MOS和VIC：NOAH是一種土地表面過(guò)程模型，既將土壤表面和植被視為一個(gè)整體的模型；CLM包含三個(gè)重要的子模型：美國(guó)國(guó)家大氣研究中心的地表模型，生物圈—大氣轉(zhuǎn)移計(jì)劃和中國(guó)科學(xué)院大氣物理研究所的陸地表面模型；MOS是根據(jù)單元內(nèi)植被類型的分布，將每個(gè)模型網(wǎng)格單元?jiǎng)澐值哪Ｐ?；VIC模型是基于熱平衡和物理動(dòng)力機(jī)制的大規(guī)模半分布式水文模型.本研究使用的數(shù)據(jù)空間分辨率為1°×1°，時(shí)間分辨率為1個(gè)月，選取的數(shù)據(jù)時(shí)間范圍與GRACE數(shù)據(jù)一致.GLDAS模型具有不同深度的土壤水分，VIC模擬3層，最大深度1.90 m；NOAH模擬4層，最大深度2.00 m；CLM模擬10層，最大深度3.43 m；MOS模擬3層，最大深度3.50 m.本研究模型在計(jì)算TWS時(shí)均考慮了土壤水分、雪水當(dāng)量和總冠層蓄水量，并且每月數(shù)據(jù)均扣除了2004年1月—2009年12月數(shù)據(jù)平均值獲得TWS相對(duì)變化量.

WGHM(Water GAP Global Hydrology Model)水文模型由德國(guó)法拉克福大學(xué)自然地理研究所研制(D?ll et al., 2014)，不同于GLDAS水文模型，該模型包含地下水成分，時(shí)間跨度范圍為2000年1月—2016年12月，空間分辨率為0.5°×0.5°，每月一值.本研究選取了2002年4月—2016年12月的數(shù)據(jù)，為與CSR-SHc數(shù)據(jù)空間分辨率保持一致，對(duì)WGHM水文模型數(shù)據(jù)進(jìn)行了插值.

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

中國(guó)地面降水月值數(shù)據(jù)集是基于國(guó)家氣象中心基礎(chǔ)資料專項(xiàng)最新整編中國(guó)地面臺(tái)站降水資料(趙煜飛等, 2014).該數(shù)據(jù)利用薄盤樣條法進(jìn)行空間插值，從而生成1961年至今的分辨率0.5°×0.5°的中國(guó)降水月值格點(diǎn)數(shù)據(jù)，數(shù)據(jù)范圍為72°E—136°E和18°N—54°N.本文選取的數(shù)據(jù)時(shí)間范圍為2002年4月—2017年1月，與CSR-SHc數(shù)據(jù)時(shí)間范圍一致.

3 新型尺度因子校正法的驗(yàn)證

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

3.1.1 截?cái)嗉盀V波引起的信號(hào)衰減

GRACE球諧系數(shù)截?cái)嘁约叭l帶和高斯平滑濾波處理會(huì)導(dǎo)致GRACE反演結(jié)果存在泄露誤差(馮偉等, 2017)，主要分為兩種情況：(1)當(dāng)研究區(qū)內(nèi)部信號(hào)較強(qiáng)，而周圍信號(hào)較弱時(shí)，球諧系數(shù)截?cái)嘁约盀V波后，研究區(qū)內(nèi)部的信號(hào)會(huì)泄露到研究區(qū)外，則表明該區(qū)域由外泄露誤差控制，稱為外泄露誤差；(2)當(dāng)研究區(qū)內(nèi)部信號(hào)較弱，而周圍信號(hào)較強(qiáng)時(shí)，球諧系數(shù)截?cái)嘁约盀V波后，周圍的信號(hào)會(huì)泄露到研究區(qū)內(nèi)部，則表明該區(qū)域由內(nèi)泄露誤差控制，稱為內(nèi)泄露誤差.

本研究利用GLDAS-NOAH水文模型計(jì)算出TWS空間分布趨勢(shì)變化，扣除研究區(qū)域長(zhǎng)江流域外的數(shù)據(jù)，模擬長(zhǎng)江流域球諧系數(shù)截?cái)嗉盀V波引起的外泄露誤差.圖2a是由2002年4月—2017年1月的GLDAS-NOAH數(shù)據(jù)計(jì)算的長(zhǎng)江流域TWS空間分布趨勢(shì)變化.將圖2a中的數(shù)據(jù)經(jīng)過(guò)球諧展開并截?cái)嘀?0階，結(jié)果如圖2b所示.圖2c和圖2d分別采用300 km和500 km高斯平滑濾波處理所得.可以看出圖2a中A1區(qū)域的取值范圍為0.50～1.00 cm·a-1，圖2b中A2區(qū)域的取值范圍為0.25～0.50 cm·a-1，圖2c中A3區(qū)域的取值范圍為0.25～0.30 cm·a-1，圖2d中A4區(qū)域的取值范圍為0.10～0.15 cm·a-1，從A1→A2→A3→A4標(biāo)注的紅色區(qū)域不斷向周圍擴(kuò)散信號(hào)，而且該區(qū)域的TWS空間分布趨勢(shì)變化不斷變小.可以得出經(jīng)球諧系數(shù)展開及截?cái)嗪笮盘?hào)向周圍擴(kuò)散，再經(jīng)過(guò)高斯平滑濾波后信號(hào)進(jìn)一步衰減，并且高斯平滑半徑越大，信號(hào)損失越強(qiáng).

圖2 GLDAS-NOAH模型計(jì)算的長(zhǎng)江流域TWS空間分布趨勢(shì)變化Fig.2 The spatial distribution trends of TWS in the Yangtze River basin calculated by the GLDAS-NOAH model

為了評(píng)估截?cái)嗉盀V波的外泄露誤差效應(yīng)，分別計(jì)算了長(zhǎng)江流域內(nèi)的區(qū)域平均趨勢(shì)變化進(jìn)行對(duì)比分析，定義趨勢(shì)變化損失率(Rate)為

Rate=(NOAH-Model)/NOAH,

(10)

其中，NOAH為GLDAS-NOAH水文模型的區(qū)域平均趨勢(shì)，Model為水文模型濾波后的區(qū)域平均趨勢(shì).

據(jù)表1可知，初始數(shù)據(jù)計(jì)算長(zhǎng)江流域的區(qū)域平均趨勢(shì)為0.08 cm·a-1；經(jīng)過(guò)球諧展開并截?cái)嘀?0階，不進(jìn)行高斯平滑的區(qū)域平均趨勢(shì)為0.07 cm·a-1，相對(duì)于初始數(shù)據(jù)信號(hào)損失12%；經(jīng)過(guò)球諧展開并截?cái)嘀?0階，進(jìn)行300 km的高斯平滑的區(qū)域平均趨勢(shì)為0.06 cm·a-1，相對(duì)于初始數(shù)據(jù)信號(hào)損失25%；經(jīng)過(guò)球諧展開并截?cái)嘀?0階，進(jìn)行500 km的高斯平滑的區(qū)域平均趨勢(shì)為0.05 cm·a-1，相對(duì)于初始數(shù)據(jù)信號(hào)損失37%.這些數(shù)據(jù)也驗(yàn)證了截?cái)嘁约盀V波造成的信號(hào)衰減.因此，恢復(fù)GRACE數(shù)據(jù)由于截?cái)嘁约盀V波造成的信號(hào)損失，對(duì)于利用GRACE球諧數(shù)據(jù)反演研究區(qū)域質(zhì)量變化具有重要意義.

表1 不同實(shí)驗(yàn)?zāi)Ｐ陀?jì)算的區(qū)域平均趨勢(shì)變化Table 1 Regional average trend changes calculated by different experimental models

3.1.2 尺度因子模型

本文使用了GLDAS四種水文模型和WGHM計(jì)算的尺度因子模型，基于新型尺度因子校正法得出的模型如圖3a所示.全球水文模型CLM、MOS、NOAH、VIC和WGHM數(shù)據(jù)基于公式(5)應(yīng)用于1°×1°的所有網(wǎng)格點(diǎn)，扣除邊界外數(shù)據(jù)得到長(zhǎng)江流域網(wǎng)格化的尺度因子，如圖3b、3d、3e、3f和3g所示.為比較尺度因子結(jié)果，本文使用了由官方提供的CLM4.0尺度因子(Landerer and Swenson, 2012)，將CLM4.0尺度因子保留長(zhǎng)江流域區(qū)域內(nèi)數(shù)據(jù)如圖3c所示.尺度因子CLM-SF(圖3b)、MOS-SF(圖3d)、NOAH-SF(圖3e))和VIC-SF(圖3f)空間分布類似，NSFCM-SF(圖3a)、CLM4.0-SF(圖3c)和WGHM-SF(圖3c)相比于圖3b、3d、3e、3f在長(zhǎng)江下游的空間分布明顯不同.這些差異也反映了不同模型TWS數(shù)據(jù)空間分布的差異，這是由于不同水文模型數(shù)據(jù)考慮的組成成分不同造成.

圖3a—3f模型的尺度因子網(wǎng)格值的空間分布基本都大于1，而圖3g模型的尺度因子網(wǎng)格值的空間分布基本都小于1，這是由于WGHM水文模型與GLDAS水文模型相比包含了地下水成分，不同模型的尺度因子顯示出明顯的差異和空間異質(zhì)性.本文利用不同模型的尺度因子計(jì)算了長(zhǎng)江流域內(nèi)空間平均尺度因子如表2所示，長(zhǎng)江流域的GLDAS模型空間平均尺度因子都大于1，其中CLM4.0-SF空間平均尺度因子最大(1.42)，WGHM-SF空間平均尺度因子最小(0.88)，空間平均尺度因子從一定程度反映了對(duì)校正結(jié)果周年振幅和長(zhǎng)期趨勢(shì)的影響.

表2 不同模型計(jì)算的長(zhǎng)江流域空間平均尺度因子Table 2 Spatial mean scale factor of The Yangtze River Basin calculated by different models

圖3 不同模型計(jì)算的長(zhǎng)江流域尺度因子Fig.3 Scale factors of the Yangtze River Basin calculated by different models

3.1.3 迭代恢復(fù)法

迭代恢復(fù)法由Chen等(2006)首次提出，經(jīng)過(guò)改進(jìn)后，廣泛應(yīng)用于GRACE數(shù)據(jù)泄露誤差校正.迭代恢復(fù)法的具體步驟如下：(1)在1°×1°網(wǎng)格的每個(gè)網(wǎng)格點(diǎn)上(圖4b)，令初始質(zhì)量變化率為GRACE的質(zhì)量變化率(圖4a).(2)正演模型質(zhì)量變化率(圖4c)由步驟(1)的1°×1°網(wǎng)格化模型質(zhì)量變化率(圖4b)，進(jìn)行球諧展開，并將球諧系數(shù)截?cái)嘀?0階，0階和1階的系數(shù)設(shè)為0，然后應(yīng)用300 km的高斯平滑濾波得到.(3)在每個(gè)格網(wǎng)點(diǎn)上，將初始的GRACE的質(zhì)量變化率(圖4a)和正演模型質(zhì)量變化率(圖4c)之間的差值乘以1.20，回加到1°×1°網(wǎng)格化模型質(zhì)量變化率(圖4b)中，作為新的模型質(zhì)量變化率進(jìn)行步驟(2)，重復(fù)該過(guò)程.連續(xù)迭代后產(chǎn)生正演模型質(zhì)量變化率和初始GRACE質(zhì)量變化率之間的一致性增加(圖4a與圖4c).(4)當(dāng)初始的GRACE的質(zhì)量變化率(圖4a)和正演模型質(zhì)量變化率(圖4c)之間的差值小于特定的閾值或者迭代特定次數(shù)后，停止迭代.

經(jīng)過(guò)迭代60次后，將得到球諧信號(hào)轉(zhuǎn)化為格網(wǎng)信號(hào)即為迭代恢復(fù)法恢復(fù)的信號(hào)結(jié)果(圖4b)，比較圖4a和圖4c發(fā)現(xiàn)具有較強(qiáng)的一致性，他們之間的差值(圖4d)較小(色棒值-0.10～0.10 cm·a-1).因此，認(rèn)為迭代恢復(fù)法恢復(fù)的結(jié)果信號(hào)是可信的.

3.2 新型尺度因子校正方法對(duì)比分析

3.2.1 長(zhǎng)江流域水儲(chǔ)量變化空間分布趨勢(shì)對(duì)比分析

根據(jù)GLDAS水文模型、WGHM水文模型和新型尺度因子校正法計(jì)算的尺度因子和官方提供的CLM4.0尺度因子獲得校正后的GRACE數(shù)據(jù)時(shí)變信號(hào)如圖5所示.可以看出，尺度因子法(WGHM水文模型除外)和迭代恢復(fù)法(圖4b)校正結(jié)果相比于CSR-SHc反演結(jié)果(圖5e)空間分布趨勢(shì)信號(hào)明顯得到增強(qiáng)，尤其是在長(zhǎng)江中下游.

為進(jìn)一步對(duì)比不同泄露誤差校正方法恢復(fù)結(jié)果，本文以CSR-M數(shù)據(jù)校正結(jié)果作為CSR-SHc數(shù)據(jù)的真實(shí)信號(hào).不同泄露誤差校正方法恢復(fù)結(jié)果空間分布趨勢(shì)與CSR-M數(shù)據(jù)空間分布趨勢(shì)相比，NSFCM-Trend(圖5a)、CLM-Trend(圖5b)、MOS-Trend(圖5g)、NOAH-Trend(圖5h)、VIC-Trend(圖5i)與CSR-M-Trend(圖5d)一致，而CLM4.0-Trend(圖5c)標(biāo)注區(qū)域C與CSR-M-Trend相比(圖5d)明顯存在過(guò)量恢復(fù)現(xiàn)象.WGHM-Trend(圖5f)標(biāo)注區(qū)域D和FM-Trend(圖4b)標(biāo)注區(qū)域B與CSR-M-Trend相比(圖5d)明顯存在信號(hào)恢復(fù)不全現(xiàn)象.對(duì)于CSR-SHc-Trend(圖5e)標(biāo)注E區(qū)域，GLDAS模型增強(qiáng)了該區(qū)域的陸地水儲(chǔ)量虧損信號(hào)，CLM4.0-Trend(圖5c)和WGHM-Trend(圖5f)在該區(qū)域的信號(hào)與CSR-M-Trend接近，NSFCM-Trend(圖5a)與GLDAS模型相比在該區(qū)域的信號(hào)更接近CSR-M-Trend.本文也計(jì)算了長(zhǎng)江流域不同泄露誤差校正方法恢復(fù)結(jié)果的區(qū)域平均趨勢(shì)如表3所示.由區(qū)域平均趨勢(shì)可以得出CSR-M-Trend最大，WGHM-Trend最小，NSFCM-Trend區(qū)域平均趨勢(shì)大于GLDAS模型計(jì)算的尺度因子校正結(jié)果的區(qū)域平均趨勢(shì)，但是小于CLM4.0-Trend和FM-Trend區(qū)域平均趨勢(shì).綜上所述，NSFCM-Trend避免了FM-Trend、CLM4.0-Trend和WGHM-Trend存在信號(hào)恢復(fù)不全和過(guò)量恢復(fù)現(xiàn)象，并且NSFCM-Trend區(qū)域平均趨勢(shì)比GLDAS水文模型尺度因子校正結(jié)果的區(qū)域平均趨勢(shì)較好.

圖4 迭代恢復(fù)法恢復(fù)GRACE質(zhì)量變化信號(hào)Fig.4 Forward modeling method to restore GRACE quality change signal

圖5 長(zhǎng)江流域不同恢復(fù)方法的空間分布趨勢(shì)Fig.5 Spatial distribution trends of different restoration methods in the Yangtze River basin

表3 長(zhǎng)江流域不同恢復(fù)方法的區(qū)域平均趨勢(shì)Table 3 Regional average trend of different restoration methods in the Yangtze River basin

3.2.2 長(zhǎng)江流域水儲(chǔ)量變化時(shí)間序列對(duì)比分析

圖6顯示了長(zhǎng)江流域不同泄露誤差校正方法恢復(fù)結(jié)果得到的水儲(chǔ)量變化時(shí)間序列.不同泄露誤差校正方法恢復(fù)結(jié)果與CSR-M數(shù)據(jù)相比周年相位一致性較好，但周年振幅具有明顯差異，表明尺度因子法和迭代恢復(fù)法主要改變水儲(chǔ)量變化時(shí)間序列的周年振幅，而對(duì)水儲(chǔ)量時(shí)間序列的周年相位影響較小.為了進(jìn)一步比較不同泄露誤差校正方法恢復(fù)結(jié)果的差異，本文計(jì)算了長(zhǎng)江流域水儲(chǔ)量變化時(shí)間序列的周年振幅和長(zhǎng)期趨勢(shì)如圖7所示.

圖6 長(zhǎng)江流域不同泄露誤差校正方法恢復(fù)的TWS時(shí)間序列Fig.6 TWS time series recovered by different leakage error correction methods in the Yangtze River Basin

圖7 長(zhǎng)江流域不同泄露誤差恢復(fù)方法恢復(fù)的TWS周年振幅和長(zhǎng)期趨勢(shì) (a) 周年振幅； (b) 長(zhǎng)期趨勢(shì).Fig.7 Annual amplitude and long-term trend of TWS recovered by different leakage error recovery methods in the Yangtze River Basin (a) Annual amplitude; (b) Long-term trends.

圖7顯示了長(zhǎng)江流域不同泄露誤差校正方法恢復(fù)的水儲(chǔ)量變化長(zhǎng)期趨勢(shì)和周年振幅，可以得出CSR-SHc數(shù)據(jù)的周年振幅和長(zhǎng)期趨勢(shì)處于較低水平，存在信號(hào)泄露現(xiàn)象.泄露誤差校正后，周年振幅和長(zhǎng)期趨勢(shì)均得到增強(qiáng)，信號(hào)泄露得到一定程度上的恢復(fù).在周年振幅上，TWS-WGHM、CSR-M數(shù)據(jù)校正結(jié)果小于CSR-SHc數(shù)據(jù)，GLDAS水文模型和CLM4.0尺度因子校正結(jié)果均明顯大于CSR-SHc數(shù)據(jù).因?yàn)镚RACE原始觀測(cè)數(shù)據(jù)模擬的TWS在精細(xì)的空間比例上不敏感(Landerer and Swenson, 2012)，而CSR-M和CSR-SHc主要受GRACE衛(wèi)星原始觀測(cè)數(shù)據(jù)約束，因此CSR-M、CSR-SHc和基于FM法恢復(fù)的周年振幅較小；GLDAS水文模型模擬的TWS大多數(shù)變化發(fā)生在次季節(jié)與季節(jié)的時(shí)間尺度上，因此由模型模擬的TWS計(jì)算的尺度因子優(yōu)化了周年振幅.在長(zhǎng)期趨勢(shì)上，WGHM模型計(jì)算的尺度因子降低了CSR-SHc數(shù)據(jù)的長(zhǎng)期趨勢(shì)；FM法、CLM4.0尺度因子、NSFCM尺度因子和GLDAS水文模型計(jì)算的尺度因子優(yōu)化了CSR-SHc數(shù)據(jù)的長(zhǎng)期趨勢(shì).其中，TWS-FM(0.34±0.08 cm·a-1)、TWS-CLM4.0數(shù)據(jù)(0.34±0.11 cm·a-1)和TWS-NSFCM(0.29±0.09 cm·a-1)較好，更接近于CSR-M數(shù)據(jù)(0.39±0.07 cm·a-1)的長(zhǎng)期趨勢(shì).

4 新型尺度因子校正法的應(yīng)用

4.1 長(zhǎng)江子流域水儲(chǔ)量變化對(duì)比分析

NSFCM-SF尺度因子綜合GLDAS水文模型計(jì)算的尺度因子校正結(jié)果空間分布趨勢(shì)的優(yōu)點(diǎn)，避免了CLM4.0-SF、FM和WGHM校正結(jié)果空間分布趨勢(shì)不均勻的現(xiàn)象.在長(zhǎng)期趨勢(shì)上，NSFCM-SF尺度因子校正結(jié)果雖然低于CLM4.0-SF和FM，但優(yōu)于GLDAS水文模型校正結(jié)果計(jì)算的長(zhǎng)期趨勢(shì)，在周年振幅上，NSFCM-SF尺度因子校正結(jié)果明顯優(yōu)于FM和CSR-M的結(jié)果.綜合考慮，本文以下主要討論NSFCM-SF尺度因子校正結(jié)果.選取CSR-M、CSR-SHc、CLM4.0和TWS-NSFCM模型，進(jìn)一步分析NSFCM-SF在長(zhǎng)江子流域(長(zhǎng)江上游和長(zhǎng)江中下游)的校正結(jié)果.

圖8顯示了長(zhǎng)江子流域CSR-M、CSR-SHc、TWS-CLM4.0和TWS-NSFCM模型的水儲(chǔ)量時(shí)間序列.據(jù)圖8可知，長(zhǎng)江子流域的TWS-NSFCM和TWS-CLM4.0周年振幅明顯優(yōu)于CSR-SHc的結(jié)果，長(zhǎng)江上游TWS-NSFCM周年振幅優(yōu)于TWS-CLM4.0，長(zhǎng)江中下游TWS-CLM4.0周年振幅優(yōu)于TWS-NSFCM.長(zhǎng)江中下游的時(shí)間序列相比于長(zhǎng)江上游更復(fù)雜，季節(jié)性的TWS變化更大，這可能與長(zhǎng)江中下游季節(jié)性地表水存儲(chǔ)變化有關(guān)(Long et al., 2015a).本文計(jì)算的長(zhǎng)江子流域水儲(chǔ)量變化長(zhǎng)期趨勢(shì)和周年振幅如表4所示.長(zhǎng)江子流域的長(zhǎng)期趨勢(shì)，TWS-NSFCM與TWS-CLM4.0相比更接近于CSR-M數(shù)據(jù).長(zhǎng)江上游CSR-SHc長(zhǎng)期趨勢(shì)為0.07 cm·a-1，經(jīng)NSFCM-SF尺度因子校正后TWS-NSFCM長(zhǎng)期趨勢(shì)為0.14 cm·a-1，長(zhǎng)期趨勢(shì)提高了100%，長(zhǎng)江中下游CSR-SHc長(zhǎng)期趨勢(shì)為0.40 cm·a-1，經(jīng)NSFCM-SF尺度因子校正后TWS-NSFCM長(zhǎng)期趨勢(shì)為0.49 cm·a-1，長(zhǎng)期趨勢(shì)提高了23%，可以得出長(zhǎng)江流域水儲(chǔ)量上升趨勢(shì)主要集中在長(zhǎng)江中下游.

表4 長(zhǎng)江子流域的TWS長(zhǎng)期趨勢(shì)和周年振幅Table 4 Long-term trends and annual amplitudes of TWS in the Yangtze River sub-basin

圖8 長(zhǎng)江子流域的TWS時(shí)間序列 (a) 長(zhǎng)江上游； (b) 長(zhǎng)江中下游.Fig.8 TWS time series of the Yangtze River sub-basin (a) The upper Yangtze River; (b) The middle and lower Yangtze River.

4.2 長(zhǎng)江流域降水與水儲(chǔ)量變化的比較

本文基于NSFCM-SF尺度因子校正結(jié)果時(shí)間序列與國(guó)家氣象中心的降水月值數(shù)據(jù)進(jìn)行了比較，如圖9所示.可以看出長(zhǎng)江流域以及子流域降水主要集中在5—10月，降水明顯增多時(shí)，TWS-NSFCM水儲(chǔ)量變化明顯增加，并且滯后降水?dāng)?shù)據(jù)1—2月.降水?dāng)?shù)據(jù)變化與TWS-NSFCM水儲(chǔ)量變化一致，這表明降水是長(zhǎng)江流域水儲(chǔ)量變化的主要原因.在2002年4月—2007年1月長(zhǎng)江流域和子流域的水儲(chǔ)量和降水均呈現(xiàn)上升趨勢(shì)，長(zhǎng)江流域水儲(chǔ)量上升趨勢(shì)為0.29 cm·a-1，降水上升趨勢(shì)為0.06 cm·a-1；上游水儲(chǔ)量上升趨勢(shì)為0.14 cm·a-1，降水上升趨勢(shì)為0.03 cm·a-1；中下游水儲(chǔ)量上升趨勢(shì)為0.49 cm·a-1，降水上升趨勢(shì)為0.10 cm·a-1.從而可以得出降水的上升趨勢(shì)主要集中在中下游，這與長(zhǎng)江流域水儲(chǔ)量上升趨勢(shì)主要集中在中下游的特征一致，進(jìn)一步驗(yàn)證了降水是長(zhǎng)江流域水儲(chǔ)量變化的主要原因.本文也計(jì)算了長(zhǎng)江流域年降水異常如圖10所示，由于2010年、2011年和2015年降水?dāng)?shù)據(jù)不完整，文中不討論降水?dāng)?shù)據(jù)不完整的年份.據(jù)圖10可知，多數(shù)年份降水量接近或低于多年平均值，但2012年、2014年和2016年降水量明顯高于多年平均值，這可能與近年來(lái)全球變暖、氣候極端事件頻發(fā)有關(guān)(例如厄爾尼諾等).結(jié)合圖9，降水較多的年份，TWS-NSFCM水儲(chǔ)量年際變化明顯增多，因此，長(zhǎng)江流域TWS-NSFCM水儲(chǔ)量年際變化特征與降水的年際變化規(guī)律一致.

圖10 長(zhǎng)江流域年降水異常Fig.10 Annual precipitation anomalies in the Yangtze River Basin

5 結(jié)論

GRACE衛(wèi)星的出現(xiàn)為研究大尺度區(qū)域的TWS變化提供了有效的手段.然而，GRACE數(shù)據(jù)處理中的泄露誤差是反演質(zhì)量變化的主要限制因素，傳統(tǒng)尺度因子法校正泄露誤差依賴水文模型結(jié)果，并且不同的水文模型由于模型的結(jié)果以及驅(qū)動(dòng)數(shù)據(jù)的差異，導(dǎo)致模型輸出結(jié)果間也存在一定的差異.因此，本文提出新型尺度因子校正法，旨在綜合多種尺度因子模型提高GRACE陸地水儲(chǔ)量變化的準(zhǔn)確性.

(1)針對(duì)GRACE數(shù)據(jù)在信號(hào)處理過(guò)程中存在的泄露誤差，本文在傳統(tǒng)尺度因子法的基礎(chǔ)上，根據(jù)水文模型的陸地水儲(chǔ)量與GRACE-SHc數(shù)據(jù)的均方根誤差和相關(guān)性賦予權(quán)重，構(gòu)建了新型尺度因子校正法.相比于傳統(tǒng)尺度因子法使用單一水文模型計(jì)算的尺度因子，新型尺度因子校正法充分考慮多種模型與GRACE-SHc數(shù)據(jù)一致性和準(zhǔn)確性信息，校正數(shù)據(jù)處理過(guò)程中的泄露誤差，恢復(fù)真實(shí)的水文信號(hào)，從而提高GRACE反演區(qū)域陸地水儲(chǔ)量變化的準(zhǔn)確性.這對(duì)于更好地理解區(qū)域陸地水儲(chǔ)量變化，進(jìn)一步認(rèn)識(shí)洪澇災(zāi)害以及對(duì)區(qū)域水資源管理具有重要意義.

(2)對(duì)新型尺度因子校正法在長(zhǎng)江流域進(jìn)行驗(yàn)證，空間分布趨勢(shì)上，GLDAS-SF和NSFCM-SF校正結(jié)果與CSR-M相近，而WGHM-SF、FM和CLM4.0-SF校正結(jié)果與CSR-M相比存在信號(hào)恢復(fù)不全和過(guò)量恢復(fù)現(xiàn)象.GLDAS-SF、CLM4.0-SF和NSFCM-SF校正結(jié)果較好地恢復(fù)了區(qū)域水儲(chǔ)量的周年振幅，但長(zhǎng)期趨勢(shì)恢復(fù)較弱，F(xiàn)M方法主要恢復(fù)了區(qū)域水儲(chǔ)量的長(zhǎng)期趨勢(shì)，相比于尺度因子周年振幅恢復(fù)結(jié)果較弱，而NSFCM-SF與GLDAS-SF相比，長(zhǎng)期趨勢(shì)優(yōu)于GLDAS-SF校正結(jié)果.

(3)基于新型尺度因子校正法校正結(jié)果顯示，長(zhǎng)江流域、上游和中下游在2002年4月—2017年1月水儲(chǔ)量呈現(xiàn)上升趨勢(shì)，分別為0.29 cm·a-1、0.14 cm·a-1和0.49 cm·a-1，相比于長(zhǎng)江流域、上游和中下游CSR-SHc上升趨勢(shì)0.21 cm·a-1、0.07 cm·a-1和0.40 cm·a-1，分別提高了38%、100%和23%.長(zhǎng)江流域水儲(chǔ)量上升趨勢(shì)主要集中在中下游.長(zhǎng)江流域降水上升趨勢(shì)為0.06 cm·a-1，上游降水上升趨勢(shì)為0.03 cm·a-1，中下游降水上升趨勢(shì)為0.10 cm·a-1，降水的上升趨勢(shì)主要集中在中下游，這與長(zhǎng)江流域水儲(chǔ)量上升趨勢(shì)主要集中在中下游的特征一致，進(jìn)而表明降水是長(zhǎng)江流域水儲(chǔ)量變化的主要原因.本研究方法對(duì)青藏高原冰川變化虧損信號(hào)泄露入長(zhǎng)江流域校正結(jié)果還存在不足之處，后期會(huì)根據(jù)一些實(shí)測(cè)數(shù)據(jù)進(jìn)一步改進(jìn)本研究方法.

致謝感謝噴氣推進(jìn)實(shí)驗(yàn)室提供的GRACE每月質(zhì)量網(wǎng)格陸地產(chǎn)品數(shù)據(jù)(https:∥grace.jpl.nasa.gov/data/get-data/monthly-mass-grids-land/)；感謝德克薩斯大學(xué)空間研究中心提供的GRACE網(wǎng)格數(shù)據(jù)產(chǎn)品CSR Mascon數(shù)據(jù)；感謝戈達(dá)德地球科學(xué)數(shù)據(jù)和信息服務(wù)中心提供的GLDAS水文模型數(shù)據(jù)(https:∥disc.gsfc.nasa.gov/)；感謝德國(guó)法拉克福大學(xué)自然地理研究所提供的WGHM水文模型數(shù)據(jù)；感謝國(guó)家氣象科學(xué)數(shù)據(jù)中心提供的降水?dāng)?shù)據(jù)(http:∥data.cma.cn/)；感謝GMT開源畫圖軟件和馮偉博士的開源代碼.感謝兩位匿名審稿人的意見，這些意見有助于提高文章質(zhì)量.楊帥、鄭偉和尹文杰為共同第一作者，鄭偉和尹文杰為共同通訊作者.