利用GRACE衛(wèi)星Level-1B數(shù)據(jù)反演陸地水儲(chǔ)量變化的方法研究

郭飛霄1,2，肖云1,2，汪菲菲3，苗岳旺3

(1.西安測繪研究所，陜西西安710054;2.地理信息工程國家重點(diǎn)實(shí)驗(yàn)室，陜西西安710054;3.西安測繪信息技術(shù)總站，陜西西安710054)

摘要：衛(wèi)星重力探測技術(shù)為監(jiān)測陸地水儲(chǔ)量變化提供了新的技術(shù)手段。對(duì)利用GRACE衛(wèi)星Level-1B數(shù)據(jù)反演地球陸地水儲(chǔ)量變化的重力位差法和Mascon方法的數(shù)學(xué)模型作了詳細(xì)推導(dǎo)分析，總結(jié)兩種方法的特點(diǎn)和解算處理步驟。推導(dǎo)過程表明：重力位差法和Mascon方法在反演時(shí)只采用衛(wèi)星飛臨研究區(qū)域上空時(shí)的觀測數(shù)據(jù)，能夠提高反演結(jié)果的空間分辨率，比傳統(tǒng)的球諧系數(shù)法更具優(yōu)勢；Mascon方法在解算時(shí)還引入了時(shí)空約束方程，進(jìn)一步提高了反演結(jié)果的時(shí)空分辨率。

關(guān)鍵詞：衛(wèi)星重力探測；陸地水儲(chǔ)量；GRACE衛(wèi)星；Mascon方法；重力位差

中圖分類號(hào)：TP72;P228

收稿日期：2014-05-21

基金項(xiàng)目：國家自然科學(xué)基金資助項(xiàng)目(41374083)；國家重點(diǎn)基礎(chǔ)研究發(fā)展計(jì)劃 (2013CB733303-2)；大地測量與地球動(dòng)力學(xué)國家重點(diǎn)實(shí)驗(yàn)室開放基金資助項(xiàng)目(SKLGED2013-3-3-E)

作者簡介：郭飛霄(1988-)，男，碩士.

A study of the methods of recovering continental water storage variation using GRACE Level-1B data

GUO Fei-xiao1,2,XIAO Yun1,2,WANG Fei-fei3,MIAO Yue-wang3

(1.Xi’an Research Institute of Surveying and Mapping,Xi’an 710054,China;2.State Key Laboratory of Geo-Information Engineering,Xi’an 710054,China;3.Technical Division of Surveying and Mapping,Xi’an 710054,China.)

Abstract：Technology of satellite gravimetry provides a new technical means for monitoring continental water storage variation.Mathematical models of Mascon and situ geopotential using GRACE Level-1B data for recovering continental water storage variation are deduced detailedly.The characteristic and processing steps of the two methods are summarized.It is showed that methods of Mascon and situ geopotential can promote the spatial resolution of result and have more advantages than the method of spheric harmonic coefficient because of using range rate observations only when satellites flying over study area.Constraint equation is introduced into Mascon method processing,and further promotes the spatial and temporal resolution.

Key words：satellite gravimetry；continental water storage；GRACE satellite；mascon method；situ geopotential

地球重力場及其時(shí)變效應(yīng)反映了地球表層及內(nèi)部物質(zhì)的重新分布和質(zhì)量遷移變化。當(dāng)前，全球性環(huán)境問題如海平面上升、極地冰川融化等與地球表層的物質(zhì)遷移緊密相關(guān)。因此，研究地球系統(tǒng)的質(zhì)量遷移和重新分布對(duì)監(jiān)測全球環(huán)境和氣候變化具有重要的意義[1]。水資源變化是每個(gè)國家經(jīng)濟(jì)和社會(huì)發(fā)展最關(guān)心的問題之一，陸地水儲(chǔ)量變化的準(zhǔn)確估計(jì)對(duì)氣候變化的研究和預(yù)報(bào)、農(nóng)業(yè)生產(chǎn)和防災(zāi)減災(zāi)等方面具有重要的現(xiàn)實(shí)意義。我國是水資源匱乏國家，近年來隨著人口的持續(xù)增長和經(jīng)濟(jì)的高速發(fā)展，水資源匱乏問題逐步凸顯。因此，監(jiān)測陸地水儲(chǔ)量變化趨勢將有著巨大的經(jīng)濟(jì)效益和社會(huì)意義。但是，地基臺(tái)站觀測覆蓋范圍小且空間分布不均勻，傳統(tǒng)的衛(wèi)星遙感技術(shù)也只能觀測到地表十幾厘米厚度的土壤濕度變化[2]。2002年3月，GRACE衛(wèi)星的成功發(fā)射開創(chuàng)了高時(shí)空分辨率監(jiān)測全球重力場的新紀(jì)元，也為監(jiān)測陸地水儲(chǔ)量變化提供了新的技術(shù)手段。衛(wèi)星重力探測技術(shù)作為遙感的一種新技術(shù)[3]，適用于中長空間尺度陸地水儲(chǔ)量變化的監(jiān)測，且觀測尺度統(tǒng)一、全球分布均勻，在中長空間尺度上可有效彌補(bǔ)遙感衛(wèi)星只能觀測地表土壤濕度以及地基觀測臺(tái)站空間分布不均勻等不足，目前已廣泛應(yīng)用于全球及區(qū)域陸地水儲(chǔ)量變化研究[4-7]。

當(dāng)前，利用GRACE衛(wèi)星數(shù)據(jù)反演陸地水儲(chǔ)量變化大多采用球諧系數(shù)法[8]，該方法采用GRACE衛(wèi)星Level-2數(shù)據(jù)的時(shí)變重力場模型，反演結(jié)果的時(shí)間分辨率為30 d，空間分辨率約為800 km。為進(jìn)一步提高陸地水儲(chǔ)量變化反演結(jié)果的時(shí)空分辨率，國外學(xué)者利用星間距離變率觀測值對(duì)地球表層質(zhì)量變化異常敏感的特性，用Level-1B級(jí)數(shù)據(jù)的KBRR(K-Band Range Rate)、GPS和加速度計(jì)觀測值直接推求陸地水儲(chǔ)量變化。根據(jù)反演原理的不同有兩種方法：重力位差法[9-11]和Mascon(Mass concentration)方法[12-14]。研究結(jié)果表明：利用Level-1B數(shù)據(jù)可以有效提高反演結(jié)果的時(shí)空分辨率，反演結(jié)果時(shí)間分辨率甚至可達(dá)10 d左右，空間分辨率可提高到200 km左右[11-12]，比球諧系數(shù)法更具優(yōu)勢。國內(nèi)針對(duì)利用Level-1B數(shù)據(jù)反演陸地水儲(chǔ)量變化的研究起步較晚。本文針對(duì)利用Level-1B數(shù)據(jù)反演陸地水儲(chǔ)量變化的方法進(jìn)行了研究，詳細(xì)推導(dǎo)分析了重力位差法和Mascon方法的數(shù)學(xué)模型。

1重力位差法

兩顆衛(wèi)星間的瞬時(shí)重力位差是恢復(fù)地球重力場的重要觀測量，可以建立起與地球重力場的直接顯式關(guān)系式，獲得較高精度的地球重力場模型。并且，利用重力位差觀測量還可以高頻率直接感知獲取區(qū)域重力場信號(hào)隨時(shí)間的變化，進(jìn)而反演地球表層區(qū)域的質(zhì)量變化[15]。

1.1基本原理

根據(jù)能量守恒原理，推導(dǎo)兩顆衛(wèi)星的瞬時(shí)重力位差觀測方程如下[16]:

(1)

(2)

地球表層的質(zhì)量遷移變化都會(huì)引起兩顆衛(wèi)星重力位差的變化，兩顆衛(wèi)星的重力位差可表示為如下形式:

(3)

行星引力、潮汐、大氣和海洋等效應(yīng)當(dāng)前已有較好的模型去去除，因此在去除這些效應(yīng)引起的重力位差后，由式(3)可以得到

(4)

1.2觀測方程

地球表層區(qū)域的物質(zhì)質(zhì)量遷移變化會(huì)引起地球重力場的局部變化，對(duì)于陸地區(qū)域，地球質(zhì)量變化主要是陸地水儲(chǔ)量的變化。因此，根據(jù)牛頓萬有引力定律，采用點(diǎn)質(zhì)量模型，如圖1所示將研究區(qū)域劃分成一定大小的網(wǎng)格，t時(shí)刻某區(qū)域陸地水儲(chǔ)量變化引起的兩顆衛(wèi)星重力位差[17]

(5)

其中:

(6)

式中:lj為網(wǎng)格j中心點(diǎn)到衛(wèi)星的距離;網(wǎng)格(r1,θ1,λ1)為前星的地心向徑、余緯和經(jīng)度;(r2,θ2,λ2)為后星的地心向徑、余緯和經(jīng)度;Δm(θj,λj,t)為第j個(gè)網(wǎng)格區(qū)域t時(shí)刻的水儲(chǔ)量變化;(θj,λj)為第j個(gè)區(qū)域網(wǎng)格中心點(diǎn)的余緯和經(jīng)度;R為地球半徑;G為萬有引力常數(shù);n為網(wǎng)格總數(shù)。

進(jìn)一步，Δm(θj,λj,t)可表示為

(7)

式中:ρw為水的密度，(RΔθ)(RΔλsinθj)為網(wǎng)格區(qū)域j的水平面積，Δhj(t)為網(wǎng)格j內(nèi)質(zhì)量變化的等效水柱高度。將式(7)代入式(5)可得

(8)

式(4)和式(8)表明，可以通過星載GPS觀測值、KBRR觀測值和加速度計(jì)數(shù)據(jù)求得兩顆衛(wèi)星的重力位差，進(jìn)而反演陸地水儲(chǔ)量變化。

圖1　重力位差法反演陸地水儲(chǔ)量變化原理

1.3參數(shù)估計(jì)

將式(8)可寫成誤差方程組形式為

(9)

采用衛(wèi)星測量數(shù)據(jù)確定地球表面物理量向下延拓的過程是個(gè)病態(tài)問題。因此，采用最小二乘平差解算上述誤差方程組，解算時(shí)法方程組是病態(tài)的。引入?yún)?shù)的先驗(yàn)信息有助于解決該病態(tài)問題[9]。設(shè)待估參數(shù)向量x的先驗(yàn)期望值為x0，Cx為先驗(yàn)協(xié)方差陣，于是可得

(10)

(11)

2Mascon方法

Mascon方法最初被應(yīng)用于研究月球表面質(zhì)量異常導(dǎo)致的引力擾動(dòng)[18]，Rowland將該方法應(yīng)用于GRACE衛(wèi)星KBRR觀測數(shù)據(jù)恢復(fù)區(qū)域時(shí)變重力場的研究[12]。目前，Mascon方法已在多個(gè)領(lǐng)域取得了豐富的研究成果，是當(dāng)前GRACE時(shí)變重力場研究的熱點(diǎn)。

2.1基本原理

(12)

σj(t)可寫成σj(t)=10×Hj(t)。比例因子10 的含義是1 m2區(qū)域內(nèi)質(zhì)量增加或減少10 kg等效為區(qū)域內(nèi)水柱升高或減少1 cm，Hj(t)為t時(shí)刻網(wǎng)格j內(nèi)質(zhì)量變化的等效水柱高度，單位為cm，Hj(t)即Mascon參數(shù)。由于重力衛(wèi)星對(duì)地球表層區(qū)域質(zhì)量變化的時(shí)間分辨率是有限的，故認(rèn)為在一定時(shí)間段內(nèi)同一網(wǎng)格區(qū)域j的Mascon參數(shù)Hj(t)是不變的。因此式(12)又可改寫為如下形式：

(13)

(14)

對(duì)研究區(qū)域內(nèi)所有網(wǎng)格采取以上相同處理，可得到式(15)，式中N1為Mascon參數(shù)個(gè)數(shù)。

(15)

2.2觀測方程

根據(jù)基于星間距離變率數(shù)據(jù)確定地球重力場模型的數(shù)學(xué)推導(dǎo)[19]，可建立星間距離變率觀測值與Mascon參數(shù)間的關(guān)系式如下：

(16)

(17)

研究結(jié)果表明，地球表層區(qū)域的質(zhì)量變化對(duì)GRACE衛(wèi)星KBRR觀測數(shù)據(jù)的影響在空間范圍是有限的[20]，因此Mascon方法在反演陸地水儲(chǔ)量變化時(shí)僅采用衛(wèi)星飛臨研究區(qū)域上空時(shí)的KBRR觀測數(shù)據(jù)。式(16)將Mascon參數(shù)與KBRR觀測值聯(lián)系起來，是Mascon反演陸地水儲(chǔ)量變化的觀測方程。

2.3時(shí)間-空間約束方程

在空間上鄰近的網(wǎng)格區(qū)域，受相近的環(huán)境影響，水儲(chǔ)量變化應(yīng)該在數(shù)值上是連續(xù)的，而不應(yīng)存在跳躍變化。并且對(duì)同一個(gè)網(wǎng)格區(qū)域，在鄰近時(shí)間段內(nèi)該網(wǎng)格區(qū)域內(nèi)的質(zhì)量變化數(shù)值上也是趨于一致的。因此，可引入時(shí)間-空間約束方程如下[13]：

(18)

對(duì)任意兩個(gè)Mascon參數(shù)Hn和Hk，其中n≠k，滿足上述約束方程，可以認(rèn)為兩者在數(shù)值上趨于一致的“可能性”滿足一個(gè)與它們之間的空間和時(shí)間“距離”相關(guān)的指數(shù)函數(shù)[14]。

(19)

式中:tnk是兩個(gè)Mascon參數(shù)在時(shí)間上的“距離”，dnk是兩Mascon參數(shù)所在網(wǎng)格的空間距離。T和D是表征了Mascon參數(shù)解的約束程度，與研究區(qū)域的網(wǎng)格劃分有關(guān)。從式(19)可以看出，兩個(gè)Mascon參數(shù)在時(shí)空上的“距離”越小，它們在數(shù)值上接近的可能性越大；反之，二者差異越大。

2.4參數(shù)估計(jì)

將時(shí)間-空間約束方程(18)看作虛擬觀測方程。在對(duì)Mascon參數(shù)進(jìn)行平差解算時(shí)，存在兩類觀測方程：KBRR觀測方程和虛擬觀測方程。因此，需要對(duì)虛擬觀測方程賦予合適的權(quán)因子Pnk，分析知，可設(shè)置為[21]

(20)

式中，S為虛擬觀測方程相對(duì)于KBRR觀測方程取單位權(quán)時(shí)的比例因子，取值一般是經(jīng)驗(yàn)性的，也可以采用赫爾默特方差分量估計(jì)。

聯(lián)立式(16)和式(18)，寫成誤差方程組的形式

(21)

其中:δX為待估參數(shù)，包括Mascon參數(shù)H和基線參數(shù)y0，A和B分別為觀測方程(16)和(18)線性化后的系數(shù)矩陣，L為觀測值。設(shè)P1和P2分別是KBRR觀測方程組和虛擬觀測方程組的權(quán)陣，對(duì)式(21)進(jìn)行平差解算可得

[ATP1A+BTP2B]-1[ATP1L1+BTP2L2]=δX.

(22)

解該線性方程組，即可求得Mascon參數(shù)解。

3結(jié)束語

由重力位差法和Mascon方法的推導(dǎo)過程可知，兩種方法在利用GRACE衛(wèi)星Level-1B數(shù)據(jù)反演陸地水儲(chǔ)量變化時(shí)都只采用了衛(wèi)星飛臨研究區(qū)域上空時(shí)的觀測數(shù)據(jù)，理論上能夠提高反演結(jié)果的空間分辨率。并且，Mascon方法在解算時(shí)還引入了時(shí)空約束方程，能夠進(jìn)一步提高反演結(jié)果的時(shí)間分辨率，并使得解算結(jié)果更加符合實(shí)際。與球諧系數(shù)法相比，利用GRACE衛(wèi)星Level-1B數(shù)據(jù)的反演結(jié)果時(shí)間和空間分辨率更高，具有更大的優(yōu)勢。

參考文獻(xiàn)：

[1]朱廣彬.利用GRACE位模型研究陸地水儲(chǔ)量的時(shí)變特征[D].北京：中國測繪科學(xué)研究院，2007.

[2]國家測繪地理信息局,遙感衛(wèi)星[EB/OL]，http://sasmac.sbsm.gov.cn/article/kpzs/200912/20091200059269.shtml.

[3]鐘敏，段建賓，許厚澤，等.利用衛(wèi)星重力觀測研究近5年中國陸地水量中長空間尺度的變化趨勢[J].科學(xué)通報(bào)，2009，54(9)：1290-1294.

[4]羅志才，李瓊，鐘波.利用GRACE時(shí)變重力場反演黑河流域水儲(chǔ)量變化[J].測繪學(xué)報(bào)，2012,41(5):676-682.

[5]李瓊，羅志才，鐘波．利用GRACE時(shí)變重力場探測2010年中國西南干旱陸地水儲(chǔ)量變化[J].地球物理學(xué)報(bào)，2012，41(5):676-673.

[6]刑樂林，李輝，劉冬至,等.利用GRACE時(shí)變重力場監(jiān)測中國及其周邊地區(qū)水儲(chǔ)量月變化[J].大地測量與地球動(dòng)力學(xué)，2007,27(4)：35-39.

[7]程朋根，鄭守住，陳曉勇，等.GNSS-R技術(shù)在海冰反演中的應(yīng)用[J].測繪工程，2014,23(4)：1-3.

[8]WAHRJ,MOLEVARARM,BRYANF.TimevariabilityoftheEarth’sgravityfield：hydrologicalandoceaniceffectsandtheirpossibledetectionusingGRACE[J].J.Geophys.Res,1998，103,30205-30229.

[9]HANSC,SHUMCK,BRAUNK.High-resolutioncontinentalwaterstoragerecoveryfromlow-lowsatellite-to-satellitetracking[J].JournalofGeodynamics,2005,39:11-28.

[10]HANSC,SHUMCK,JEKELIC,etal.ImprovedestimationofterrestrialwaterstoragechangesfromGRACE[J].GeophysicalResearchLetters,2005,Vol.32,L07302.

[11]王雷.GRACE時(shí)變重力場[D].武漢:武漢大學(xué)，2006.

[12]ROWLANDSDD,LUTHCKESB,KLOSKOSM,etal.ResolvingmassfluxathighspatialandtemporalresolutionusingGRACEintersatellitemeasurements[J].GeophysicalResearchLetters，2005,Vol.32,L04310.

[13]ROWLANDSDD,LUTHCKESB,MCCARTHYJJ,etal.GlobalmassfluxsolutionsfromGRACE：Acomparisonofparameterestimationstrategies-MassconcentrationsversusStokescoefficients[J].JournalofGeophysicalResearch,2010,Vol.115,B01403.

[14]張坤.基于Mascon方法確定時(shí)變地球重力場的研究[D].武漢:武漢大學(xué)，2011.

[15]王正濤，姜衛(wèi)平，晁定波.衛(wèi)星跟蹤衛(wèi)星測量確定地球重力場的理論和方法[M].武漢:武漢大學(xué)出版社，2011.

[16]JEKELIC.Thedeterminationofgravitationalpotentialdifferencesfromsatellite-to-satellitetracking[J].Celest.Mech.Dyn.Astron,1999,75:85-100.

[17]PAILR.Syntheticglobalgravitymodelforplanetarybodiesandapplicationsinsatellitegravitygradiometry[D].Ph.D.dissertation,TechnicalUniversityGraz,1999.

[18]邢樂林，李建成，沈飛,等.點(diǎn)質(zhì)模型求定月球重力場及其特征分析[J].測繪科學(xué)，2008,33(1)：36-40.

[19]肖云.基于衛(wèi)星跟蹤衛(wèi)星數(shù)據(jù)恢復(fù)地球重力場的研究[D].鄭州：信息工程大學(xué)，2006.

[20]LEMOINEFG,LUTHCKESB,ROWLANDSDD,etal.Theuseofmasconstoresolvetime-variablegravityfromGRACE[J].InternationalAssociationofGeodesySymposia,2007,130:231-236.

[21]KLOSKOS,ROWLANDSDD,LUTHCKESB,etal.EvaluationandvalidationofmasconrecoveryusingGRACEKBRRdatawithindependentmassfluxestimateintheMississippiBasin[J].JGeod,2009,83:817-827.廣彬.利用GRACE位模型研究陸地水儲(chǔ)量的時(shí)變特征[D].北京：中國測繪科學(xué)研究院，2007.

[2]國家測繪地理信息局,遙感衛(wèi)星[EB/OL]，http://sasmac.sbsm.gov.cn/article/kpzs/200912/20091200059269.shtml.

[3]鐘敏，段建賓，許厚澤，等.利用衛(wèi)星重力觀測研究近5年中國陸地水量中長空間尺度的變化趨勢[J].科學(xué)通報(bào)，2009，54(9)：1290-1294.

[4]羅志才，李瓊，鐘波.利用GRACE時(shí)變重力場反演黑河流域水儲(chǔ)量變化[J].測繪學(xué)報(bào)，2012,41(5):676-682.

[5]李瓊，羅志才，鐘波．利用GRACE時(shí)變重力場探測2010年中國西南干旱陸地水儲(chǔ)量變化[J].地球物理學(xué)報(bào)，2012，41(5):676-673.

[6]刑樂林，李輝，劉冬至,等.利用GRACE時(shí)變重力場監(jiān)測中國及其周邊地區(qū)水儲(chǔ)量月變化[J].大地測量與地球動(dòng)力學(xué)，2007,27(4)：35-39.

[7]程朋根，鄭守住，陳曉勇，等.GNSS-R技術(shù)在海冰反演中的應(yīng)用[J].測繪工程，2014,23(4)：1-3.

[8]WAHRJ,MOLEVARARM,BRYANF.TimevariabilityoftheEarth’sgravityfield：hydrologicalandoceaniceffectsandtheirpossibledetectionusingGRACE[J].J.Geophys.Res,1998，103,30205-30229.

[9]HANSC,SHUMCK,BRAUNK.High-resolutioncontinentalwaterstoragerecoveryfromlow-lowsatellite-to-satellitetracking[J].JournalofGeodynamics,2005,39:11-28.

[10]HANSC,SHUMCK,JEKELIC,etal.ImprovedestimationofterrestrialwaterstoragechangesfromGRACE[J].GeophysicalResearchLetters,2005,Vol.32,L07302.

[11]王雷.GRACE時(shí)變重力場[D].武漢:武漢大學(xué)，2006.

[12]ROWLANDSDD,LUTHCKESB,KLOSKOSM,etal.ResolvingmassfluxathighspatialandtemporalresolutionusingGRACEintersatellitemeasurements[J].GeophysicalResearchLetters，2005,Vol.32,L04310.

[13]ROWLANDSDD,LUTHCKESB,MCCARTHYJJ,etal.GlobalmassfluxsolutionsfromGRACE：Acomparisonofparameterestimationstrategies-MassconcentrationsversusStokescoefficients[J].JournalofGeophysicalResearch,2010,Vol.115,B01403.

[14]張坤.基于Mascon方法確定時(shí)變地球重力場的研究[D].武漢:武漢大學(xué)，2011.

[15]王正濤，姜衛(wèi)平，晁定波.衛(wèi)星跟蹤衛(wèi)星測量確定地球重力場的理論和方法[M].武漢:武漢大學(xué)出版社，2011.

[16]JEKELIC.Thedeterminationofgravitationalpotentialdifferencesfromsatellite-to-satellitetracking[J].Celest.Mech.Dyn.Astron,1999,75:85-100.

[17]PAILR.Syntheticglobalgravitymodelforplanetarybodiesandapplicationsinsatellitegravitygradiometry[D].Ph.D.dissertation,TechnicalUniversityGraz,1999.

[18]邢樂林，李建成，沈飛,等.點(diǎn)質(zhì)模型求定月球重力場及其特征分析[J].測繪科學(xué)，2008,33(1)：36-40.

[19]肖云.基于衛(wèi)星跟蹤衛(wèi)星數(shù)據(jù)恢復(fù)地球重力場的研究[D].鄭州：信息工程大學(xué)，2006.

[20]LEMOINEFG,LUTHCKESB,ROWLANDSDD,etal.Theuseofmasconstoresolvetime-variablegravityfromGRACE[J].InternationalAssociationofGeodesySymposia,2007,130:231-236.

[21]KLOSKOS,ROWLANDSDD,LUTHCKESB,etal.EvaluationandvalidationofmasconrecoveryusingGRACEKBRRdatawithindependentmassfluxestimateintheMississippiBasin[J].JGeod,2009,83:817-827.