蘇 勇 魏 偉 李 瓊 谷延超

1 西南石油大學(xué)土木工程與測繪學(xué)院，成都市新都大道8號，610500 2 武漢大學(xué)測繪學(xué)院，武漢市珞喻路129號，430079

陸地水作為全球水循環(huán)的重要部分，對經(jīng)濟(jì)發(fā)展、全球氣候、生態(tài)環(huán)境都有重要影響，因此監(jiān)測全球范圍內(nèi)水儲量變化具有現(xiàn)實(shí)意義。大氣壓力、海底壓力、陸地水儲量變化等會引起地球重力場變化[1]，GRACE重力衛(wèi)星可捕捉該地球物理信號的波動，為反演陸地質(zhì)量遷移提供新途徑。

目前，采用時變重力場反演地球質(zhì)量變化主要有球諧系數(shù)法[2-3]、Mascon法[4]以及點(diǎn)質(zhì)量模型法[5]。球諧系數(shù)法相比Mascon法和點(diǎn)質(zhì)量模型法算法簡單、易編程實(shí)現(xiàn)、應(yīng)用廣泛，但球諧系數(shù)法的反演結(jié)果存在信號泄露等缺陷。此外，如要達(dá)到較好的去條帶噪聲效果，濾波方法也需要進(jìn)行詳細(xì)探究[6-8]。Mascon法可采用GRACE星間距離變率數(shù)據(jù)反演地表質(zhì)量遷移，反演結(jié)果能以更高的時空分辨率揭示地球物理變化[9]，但計(jì)算過程復(fù)雜。為進(jìn)一步提高時變重力場模型反演結(jié)果的時空分辨率，且減少計(jì)算復(fù)雜程度，國內(nèi)外學(xué)者提出了各種點(diǎn)質(zhì)量模型法[2,10-13]。采用點(diǎn)質(zhì)量模型法求解地面點(diǎn)質(zhì)量變化時，法方程會呈現(xiàn)病態(tài)，無法獲得穩(wěn)定解。該問題常用的解決方法是進(jìn)行正則化，構(gòu)建正則化矩陣的方法主要包括不考慮地表物理信息的純數(shù)學(xué)的正則化約束、基于時空相關(guān)性進(jìn)行約束、基于先驗(yàn)信息進(jìn)行約束等[14-16]。正則化約束是大地測量不適定問題中使用最為廣泛的一種方法，點(diǎn)質(zhì)量模型法中常采用零階Tikhonov正則化進(jìn)行約束，該方法的原理是在法方程系數(shù)矩陣對角線增加純數(shù)學(xué)的正則化參數(shù)，以抑制法方程病態(tài)，獲得穩(wěn)定解?；跁r空相關(guān)性進(jìn)行約束是基于地理相鄰點(diǎn)之間的時間、空間相關(guān)性，通過權(quán)值定義其關(guān)系，從而構(gòu)建約束矩陣。先驗(yàn)信息約束是通過已有的參考模型或先驗(yàn)誤差信息構(gòu)建約束矩陣。這3類方法都涉及正則化參數(shù)，而該參數(shù)的選取通常會影響反演結(jié)果，過大會造成過度平滑，丟失真實(shí)地面物理信號，過小則會保留較多噪聲，無法獲得準(zhǔn)確結(jié)果。正則化參數(shù)的選取一般可使用廣義交叉驗(yàn)證法、L曲線法。

本文基于附有空間約束的三維加速度點(diǎn)質(zhì)量模型法[15]，采用水文模型構(gòu)建地面點(diǎn)之間的空間相關(guān)性矩陣，進(jìn)而抑制法方程病態(tài)，并將該方法應(yīng)用于2009年秋至2010年春云南、貴州、四川地區(qū)特大干旱監(jiān)測，在剔除季節(jié)性信號后，采用主成分分析法對2003年至2012年西南地區(qū)水儲量進(jìn)行分析。

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

1.1 附有空間約束的三維加速度點(diǎn)質(zhì)量模型法基本原理

由于文獻(xiàn)[15]已對附有空間約束的三維加速度點(diǎn)質(zhì)量模型法的基本原理進(jìn)行詳細(xì)介紹，本文在此直接給出相關(guān)矩陣的構(gòu)成方法。在構(gòu)建相關(guān)矩陣時，將空間相關(guān)性進(jìn)行量化，采用水文模型提取各點(diǎn)之間的相關(guān)系數(shù)，從而計(jì)算相關(guān)矩陣。當(dāng)質(zhì)量點(diǎn)之間的球面距離不大于相關(guān)距離時(dij≤D)，由式(1)構(gòu)成相關(guān)矩陣:

Cij=|Rij|/|R|sum

(1)

式中，|Rij|為點(diǎn)i和點(diǎn)j相關(guān)系數(shù)的絕對值，|R|sum為在相關(guān)距離內(nèi)，各質(zhì)量點(diǎn)和點(diǎn)i相關(guān)系數(shù)的絕對值之和。當(dāng)質(zhì)量點(diǎn)之間的球面距離大于相關(guān)距離D時(dij≥D)，Cij=0。由此，可以生成約束矩陣。

1.2 主成分分析法

主成分分析法(principal component analysis, PCA)可對數(shù)據(jù)進(jìn)行壓縮提取，在地學(xué)領(lǐng)域通常采用該方法獲得空間模態(tài)及時間序列。在對利用GRACE重力衛(wèi)星反演的水儲量數(shù)據(jù)進(jìn)行主成分分析時，所計(jì)算的特征向量可稱為空間特征或空間模態(tài)，而提取的主成分則認(rèn)為是時間序列，因此主成分分析也稱為時空分析或正交經(jīng)驗(yàn)分解。主成分分析法計(jì)算過程參見文獻(xiàn)[17]。

1.3 數(shù)據(jù)處理及檢驗(yàn)方法

本文采用CSR發(fā)布的60階Level-2 RL06時變重力模型數(shù)據(jù)，時間范圍為2003～2012年。由于GRACE重力衛(wèi)星對球諧系數(shù)C20不敏感，在計(jì)算中采用SLR數(shù)據(jù)替換C20項(xiàng)[18]。為對比反演結(jié)果，分別采用300 km高斯濾波的球諧系數(shù)法、零階Tikhonov正則化約束的三維加速度點(diǎn)質(zhì)量模型法反演研究區(qū)質(zhì)量變化。在構(gòu)建水文模型約束矩陣時，為減小各種水文模型的不確定性，采用GLDAS(global land data assimilation system)中2種水文模式和CPC(climate prediction center)水文模型的均值作為最終估值，并采用球諧分析與球諧綜合轉(zhuǎn)化為1°×1°格網(wǎng)值。本文利用GLDAS中VIC、NOAH模式，提取0～2 m深度的土壤濕度變化。CPC水文模型數(shù)據(jù)時間分辨率為1個月，空間分辨率為0.5°，包含0～1.6 m深度的土壤濕度信息。

為檢驗(yàn)引入水文模型后的附有空間約束的三維加速度點(diǎn)質(zhì)量模型法的反演結(jié)果，本文采用信噪比(SNR)參數(shù)進(jìn)行評估，具體公式為[15]:

(2)

式中，ΔmM為各種方法反演的地面質(zhì)量變化，ΔmR為參考模型的地面質(zhì)量變化，當(dāng)SNR>0時，表明信號大于噪聲。

選用CSR發(fā)布的RL06 Mascon模型作為參考模型，同樣進(jìn)行格網(wǎng)一致化處理。

2 結(jié)果與分析

2.1 空間約束半徑選取及各類方法效果對比

構(gòu)建約束矩陣時，以研究點(diǎn)為圓心，約束半徑內(nèi)的點(diǎn)被認(rèn)為與研究點(diǎn)具有地理相關(guān)性，因此約束半徑的選擇至關(guān)重要。在反演東亞地區(qū)質(zhì)量變化時，考慮到在反演質(zhì)量變化時地理格網(wǎng)間隔為1°，本文設(shè)置間隔為100 km，約束半徑為200～1 100 km，采用SNR進(jìn)行結(jié)果評估。由于廣義交叉函數(shù)變化平緩，確定最小值存在一定困難，會影響最優(yōu)正則化參數(shù)的確定，因此本文采用L曲線法選取最優(yōu)正則化參數(shù)。

圖1為不同約束半徑下2003～2012年東亞區(qū)域信噪比比例折線圖，即信噪比大于0的點(diǎn)位所占比例。采用統(tǒng)計(jì)方法，選擇該時間段內(nèi)峰值最為頻繁的700 km約束半徑作為最優(yōu)約束距離。需要注意的是，本文采用信噪比比例作為約束半徑的選取準(zhǔn)則，而不同空間尺度的研究區(qū)所計(jì)算的信噪比比例會存在差異，進(jìn)而會影響最優(yōu)約束半徑的選擇，因此本文在計(jì)算西南地區(qū)水儲量變化時，約束半徑也進(jìn)行重新選取。

圖1 不同約束半徑信噪比比例

圖2為各種方法反演的2003-01水儲量變化空間分布，其中圖2(a)為引入水文模型的附有空間約束的三維加速度點(diǎn)質(zhì)量模型法的反演結(jié)果(WMC-PMA)，圖2(b)為CSR Mascon產(chǎn)品結(jié)果(Mascon)，圖2(c)為采用零階Tikhonov正則化進(jìn)行約束的三維加速度點(diǎn)質(zhì)量模型法的反演結(jié)果(ZOT)，圖2(d)為球諧系數(shù)法的反演結(jié)果(SH)。圖中(a)、(c)與(d)去條帶噪聲效果大致相同，但(a)、(c)中地面物理信號更強(qiáng)。(a)、(c)、(d)與(b)相比存在一定區(qū)別，但在我國黃河下游區(qū)域、印度北部、緬甸、珠穆朗瑪峰地區(qū)都有明顯的質(zhì)量變化。在信號方面，(b)在保留高空間分辨率情況下，去條帶噪聲效果更好。(a)與(c)相比在空域上區(qū)別較小，為進(jìn)一步對比2種約束方法的效果，利用信噪比對2種方法進(jìn)行評估，同時給出球諧系數(shù)法的相應(yīng)結(jié)果。

圖2 2003-01東亞區(qū)域水儲量變化

圖3中WMC-PMA方法(約束距離取700 km)相比ZOT方法信噪比大于0的比例更高，其在多數(shù)月份均優(yōu)于ZOT方法。SH方法和WMC-PMA方法信噪比大于0的比例在多數(shù)月份相差不大，但由圖2可知，在去條帶噪聲效果差異較小的情況下，后者在部分地區(qū)具有更明顯的地面質(zhì)量變化信號，如我國華北區(qū)域。因此，本文方法在保留較好信噪比的情況下，也可保留較強(qiáng)的地面物理信號。

圖3 各類方法信噪比比例

2.2 西南地區(qū)干旱監(jiān)測

2009年秋至2010年春云南、貴州、四川發(fā)生特大干旱。去除季節(jié)性信號后，采用本文方法計(jì)算2003～2012年中國西南地區(qū)水儲量變化(在后續(xù)計(jì)算分析中均采用已去除季節(jié)性信號的結(jié)果進(jìn)行分析)，水儲量變化趨勢空間分布如圖4所示，水儲量變化時間序列如圖5所示。

圖4 2003～2012年西南地區(qū)水儲量變化趨勢

圖5 2003～2012年西南地區(qū)水儲量變化時間序列

從圖4可以看出，昆明、貴陽、成都具有質(zhì)量增加趨勢，貴陽東北部增加趨勢尤為明顯，與已有研究結(jié)果大致相同[19]。但本文計(jì)算結(jié)果表明昆明西南方向質(zhì)量處于增加趨勢，文獻(xiàn)[19]利用GRACE重力衛(wèi)星數(shù)據(jù)反演的結(jié)果無該信號特征，而GLDAS反演結(jié)果具有質(zhì)量增加趨勢。原因可能為文獻(xiàn)[19]采用的時變重力場模型為CSR發(fā)布的RL04版本，而不同版本的數(shù)據(jù)精度和質(zhì)量均在逐步提高，因此反演結(jié)果與本文存在一定區(qū)別。

圖5為研究區(qū)水儲量變化時間序列，從圖中可以看出，2006-09存在明顯負(fù)異常，而2006-07～08云南地區(qū)為高溫?zé)o雨天氣，這是導(dǎo)致該異常的主要原因。2007年至2008年末，水儲量呈現(xiàn)上升趨勢，而2009年水儲量變化雖然存在波動，但整體呈現(xiàn)下降趨勢。2011年夏季同樣存在明顯的水儲量下降趨勢，出現(xiàn)夏季極端干旱事件，該次干旱事件的主要原因?yàn)?011年夏季降雨量相比同期明顯偏少且干旱區(qū)基本與2009年冬至2010年春一致。值得注意的是，在2009年冬，研究區(qū)水儲量處于上升趨勢，本文將在后續(xù)對該異常進(jìn)行分析說明。

為進(jìn)一步分析西南地區(qū)水儲量變化，采用主成分分析法(PCA)對反演結(jié)果進(jìn)行處理。經(jīng)過主成分分析法處理，一般可獲取部分貢獻(xiàn)率較大的成分，如這幾個成分的累積貢獻(xiàn)率超過閾值，則認(rèn)為可以反映當(dāng)前信息量。

從圖6可以看出，主成分分析法在累積貢獻(xiàn)率達(dá)到90.48%的情況下，共提取9個成分，其中前6個成分的累積貢獻(xiàn)率達(dá)到81.25%，可以解釋研究區(qū)水儲量變化。第1個成分的貢獻(xiàn)率為30.6%，遠(yuǎn)大于其他成分的貢獻(xiàn)率。西南地區(qū)水儲量主成分分析結(jié)果如圖7、8所示。圖8中所有主成分對應(yīng)的時間序列均已剔除季節(jié)性信號，因此不存在周期性特征。

圖6 PCA分解特征值及方差貢獻(xiàn)率

圖8中PCA1在2009～2010年相比同期具有明顯的下降趨勢，而對應(yīng)的空間模態(tài)則顯示貴陽、昆明具有比其他地區(qū)更明顯的正變化信號，因此該成分主要代表貴陽、昆明地區(qū)水儲量變化信號，且貴陽、昆明水儲量在該期間也處于下降趨勢。在2010年春末，貴陽地區(qū)水儲量開始恢復(fù)，這與圖9中水儲量變化的空間分布大致相同；但從圖9可以看出，昆明地區(qū)在2010年春末負(fù)異常信號增強(qiáng)，與主成分分析結(jié)果不一致。該現(xiàn)象可能是由于在計(jì)算時間序列主成分時，空間模態(tài)在貴陽地區(qū)正變化數(shù)值較大，與昆明地區(qū)信號相互作用導(dǎo)致PCA1時間序列在該時間處于上升趨勢。2003～2012年整個研究時間段內(nèi)PCA1處于增加趨勢，因此，在該期間內(nèi)貴陽、昆明地區(qū)水儲量處于增加態(tài)勢，與圖4所反演的水儲量變化趨勢相同。

圖8 西南地區(qū)主成分時間序列

PCA2空間模態(tài)在成都、長江中游地區(qū)為負(fù)變化信號，在我國沿海地區(qū)及昆明南部具有較弱的正變化信號，因此PCA2主要表示成都及長江中游地區(qū)特征。PCA2時間序列在2009年秋至2010年春處于下降趨勢，結(jié)合空間模態(tài)的正負(fù)變化信號可知，成都地區(qū)水儲量變化呈現(xiàn)增加趨勢，與PCA2時間序列的下降趨勢相反；昆明南部水儲量變化呈現(xiàn)減少趨勢，與PCA2時間序列趨勢相同。由圖9可知，在去除季節(jié)性特征后，成都地區(qū)在2009年秋至2010年受干旱影響較小，且多數(shù)月份表現(xiàn)為水儲量正異常特征，昆明南部則在干旱期內(nèi)表現(xiàn)為負(fù)異常。從研究時間段的整體趨勢來看，PCA2時間序列呈現(xiàn)減小趨勢，這表明2003～2012年成都及長江中游地區(qū)水儲量為增加趨勢，昆明南部水儲量具有減小趨勢，與圖4反演的變化趨勢相印證。值得注意的是，PCA2時間序列在2006年具有明顯的起伏特征，在迅速上升后快速下降，這與云南夏季高溫?zé)o雨具有一定關(guān)系。

圖9 西南地區(qū)2009-07～2010-06水儲量空間分布

由圖7可知，PCA3主要表示青藏高原及云南部分地區(qū)，且由該地區(qū)負(fù)變化信號可知，其與圖8中PCA3時間序列變化趨勢相反。2009年P(guān)CA3時間序列相比同期具有明顯的上升趨勢，因此青藏高原及云南部分地區(qū)在該時期表現(xiàn)為水儲量快速下降趨勢。PCA3時間序列在整個時間段內(nèi)呈現(xiàn)上升趨勢，表明青藏高原及云南部分地區(qū)在該時間段內(nèi)水儲量處于下降趨勢，與圖4中趨勢一致。

PCA4、PCA5、PCA6空間模態(tài)較為復(fù)雜，原始數(shù)據(jù)在空間模態(tài)上進(jìn)行投影后，無法結(jié)合空間模態(tài)與時間序列對水儲量信號變化進(jìn)行定量分析，但PCA1、PCA2、PCA3時間序列及空間模態(tài)已經(jīng)可以得出西南地區(qū)水儲量變化趨勢，且主成分分析法也可成功提取出2009年秋至2010年春西南地區(qū)水儲量虧損信號。

為更好地描述西南地區(qū)干旱情況，圖9為西南地區(qū)2009年秋至2010年春水儲量空間分布。從圖中可以看出，青藏高原在2009-07～2010-04水儲量一直處于負(fù)異常，在2010-02信號達(dá)到最大值，隨后減弱，并在2010-05恢復(fù)為正值。昆明地區(qū)則在秋季前期呈現(xiàn)水儲量正異常，在2009-09開始表現(xiàn)為負(fù)異常，且數(shù)值不斷增大，并在2009-11開始分別與青藏高原負(fù)異常、貴陽地區(qū)負(fù)異常融合；2010-03西南地區(qū)大部分區(qū)域呈現(xiàn)水儲量負(fù)異常，干旱情況加劇，與文獻(xiàn)[19]反演結(jié)果一致；此后，昆明地區(qū)水儲量負(fù)異常范圍開始縮小，但數(shù)值不斷增大，主要原因還是缺乏有效降雨。貴陽地區(qū)相比昆明地區(qū)受干旱影響較小，同樣在2010-03達(dá)到負(fù)異常最大值，此后負(fù)異常值逐漸減小，在2010-06水儲量迅速恢復(fù)。成都地區(qū)在此次干旱中受影響最小，多數(shù)月份表現(xiàn)為正異常。從圖5可以看出，2009年冬季研究區(qū)水儲量處于上升趨勢，結(jié)合圖9(d)～(f)可知，在該段時間內(nèi)青藏高原負(fù)異常信號有所減弱，且長江中游、秦嶺、成都西北部正異常信號增加，從而導(dǎo)致該時間段內(nèi)水儲量時間序列處于上升趨勢，但云南大部分區(qū)域?yàn)樗畠α控?fù)異常。

降雨量作為陸地水儲量的主要補(bǔ)給來源，可對干旱洪澇事件進(jìn)行輔助分析。圖10為已去除季節(jié)變化的西南地區(qū)2003～2012年熱帶降水測量計(jì)劃(tropical rainfall measuring missionn, TRMM)的降雨量距平時間序列。

圖10 西南地區(qū)降雨量距平時間序列

由圖10可知，大多年份的降雨量距平處于波動狀態(tài)，但可以明顯看出2009年降雨異常，該年降雨量距平負(fù)異常的月份集中，且數(shù)值較大，說明2009年降雨量相比同期明顯偏少，而雨水減少會削弱對陸地水儲量的補(bǔ)給，從而引發(fā)干旱，這與GRACE捕捉到的干旱事件相吻合。此外，圖10中2011年出現(xiàn)降雨量距平極端負(fù)異常的情況，說明西南地區(qū)在2011年也出現(xiàn)一定程度的旱情。

3 結(jié) 語

本文采用空間約束方法處理反演過程中的病態(tài)問題，引入水文模型將地理點(diǎn)之間的相關(guān)性進(jìn)行量化，并與球諧系數(shù)法、零階Tikhonov約束的三維加速度點(diǎn)質(zhì)量模型法進(jìn)行比較。此外，采用該方法反演2003～2012年我國西南地區(qū)水儲量變化，去除季節(jié)性變化特征后，利用主成分分析法對水儲量進(jìn)行分析。研究結(jié)果表明：

1)2種約束形式的三維加速度點(diǎn)質(zhì)量模型法均可處理病態(tài)問題，且在空間信號特征上差別較小，但本文方法在信噪比方面優(yōu)于零階Tikhonov約束的三維加速度點(diǎn)質(zhì)量模型法。

2)采用主成分分析法處理2003～2012年我國西南地區(qū)去除季節(jié)性信號后的水儲量變化數(shù)據(jù)，在累積方差貢獻(xiàn)率達(dá)到90.48%的情況下，可提取出9個主成分。結(jié)果表明，PCA1貢獻(xiàn)率最大，該成分主要代表西南地區(qū)大部分區(qū)域，其所對應(yīng)的時間序列顯示2009年秋至2010年春存在明顯的水儲量下降趨勢。

3)2009年秋至2010年春，云南地區(qū)水儲量負(fù)異常信號特征最明顯、范圍最廣、持續(xù)時間最長。貴陽地區(qū)水儲量信號在前期呈現(xiàn)正值，后期逐步呈現(xiàn)負(fù)異常，貴陽東部尤為明顯，但在2010-06，貴陽東部水儲量信號迅速恢復(fù)為正異常。成都地區(qū)在此次干旱中受影響較小。