陳鑑華，張興福，沈云中，陳秋杰，李偉超

1. 廣東工業(yè)大學(xué)測繪工程系，廣東 廣州 510006； 2. 同濟大學(xué)測繪與地理信息學(xué)院，上海 200092； 3. 波恩大學(xué)大地測量與地理信息學(xué)院，波恩 53121

2009年，歐空局(European Space Agency，ESA)發(fā)射了GOCE(gravity field and steady-state ocean circulation explorer)衛(wèi)星，其主要科學(xué)目標(biāo)是確定約100 km分辨率、1～2 cm精度的大地水準(zhǔn)面，可為固體地球物理學(xué)、海洋學(xué)、大地測量學(xué)和冰川學(xué)等科學(xué)研究提供高精度的中高階靜態(tài)地球重力場模型[1]。GOCE衛(wèi)星在任務(wù)期間共獲取了約42個月的有效重力梯度觀測值，由于重力梯度觀測值包含了時變重力場變化信號，在高階靜態(tài)重力場反演中需考慮時變重力場變化改正的影響。ESA在Level 2級別的GOCE衛(wèi)星重力梯度數(shù)據(jù)中給出了ESA計算的時變重力場變化改正值[2]，目前國內(nèi)外絕大多數(shù)科研機構(gòu)均直接利用GOCE衛(wèi)星Level 2重力梯度數(shù)據(jù)進(jìn)行重力場反演[3-5]，但時變重力場變化改正是實現(xiàn)由GOCE衛(wèi)星Level 1b級別數(shù)據(jù)到Level 2級別數(shù)據(jù)轉(zhuǎn)換的重要過程，也是完整GOCE衛(wèi)星數(shù)據(jù)處理過程中的一個重要環(huán)節(jié)，而ESA有關(guān)文檔中并未給出詳細(xì)的數(shù)據(jù)處理過程，自主實現(xiàn)由GOCE衛(wèi)星Level 1b重力梯度數(shù)據(jù)直接進(jìn)行重力場反演具有重要的科學(xué)意義，可為我國獨立開展相關(guān)GOCE類衛(wèi)星任務(wù)的數(shù)據(jù)處理研究提供重要參考。

GOCE衛(wèi)星重力梯度觀測值的時變重力場變化改正主要包括直接潮汐改正、固體潮汐改正、海洋潮汐改正、極潮改正和非潮汐質(zhì)量變化改正等[6-7]。鑒于GOCE衛(wèi)星重力梯度觀測值的時變重力場變化改正的重要性，國內(nèi)外已有學(xué)者對重力梯度的時變重力場變化改正展開研究[7-16]，分別從功率譜密度特性[13-15]以及時變重力場變化的高頻信號源[9-10]等方面分析了時變重力場變化改正對GOCE衛(wèi)星重力梯度觀測值的影響。但目前暫未見從重力場反演結(jié)果層面分析時變重力場變化改正及標(biāo)準(zhǔn)和背景模型更新對高階靜態(tài)重力場反演的影響等相關(guān)研究成果。基于上述背景，本文采用直接法進(jìn)行GOCE衛(wèi)星重力反演，并進(jìn)行有關(guān)分析，首先選用和ESA一致或同類的標(biāo)準(zhǔn)和背景模型計算出時變重力場變化改正值并與ESA結(jié)果進(jìn)行比較，驗證本文計算的正確性；然后采用新標(biāo)準(zhǔn)和背景模型對GOCE衛(wèi)星重力梯度觀測值進(jìn)行時變重力場變化改正，自主實現(xiàn)了由Level 1b重力梯度數(shù)據(jù)直接進(jìn)行重力反演；最后通過3種時變重力場變化改正方案分別反演220階次的重力場模型，從大地水準(zhǔn)面階誤差和不同方案模型間全球大地水準(zhǔn)面差異等方面來分析時變重力場變化改正對高階靜態(tài)重力場反演結(jié)果的影響，并利用實測大地水準(zhǔn)面模型對結(jié)果進(jìn)行檢核。

1 原理與方法

GOCE衛(wèi)星的主要任務(wù)是反演高精度的靜態(tài)地球重力場，但重力梯度觀測值包含了時變重力場變化信號，需要考慮對其進(jìn)行時變重力場變化改正，主要包括：直接潮汐改正、固體潮汐改正、海洋潮汐改正、極潮改正和非潮汐質(zhì)量變化改正等，其中，直接潮汐改正先計算引潮位，再將其轉(zhuǎn)換為重力梯度改正值，而固體潮汐改正、海洋潮汐改正、極潮改正、非潮汐質(zhì)量變化改正是先確定對應(yīng)的位系數(shù)改正數(shù)，然后將其轉(zhuǎn)換為重力梯度改正值。

1.1 直接潮汐改正

直接潮汐是由日月及其他天體對衛(wèi)星的引潮力引起的，對應(yīng)的引潮力位Vj可表示為[17]

(1)

1.2 固體潮汐改正

固體潮汐引起的重力位系數(shù)改正可分為兩步計算[18]。

第1步：計算固體潮汐引起的2階和3階重力位系數(shù)改正，以及2階對4階重力位系數(shù)的改正，公式為

(2)

(3)

第2步：計算2階項的日潮、半日潮及長周期潮改正，其中計算日潮改正的公式為

(4)

2階項的潮汐校正包含48個日潮、兩個半日潮以及21個長周期潮(其他具體計算公式可參考文獻(xiàn)[18])，f為各潮波分量,Ampf(ip/op)和其他相關(guān)參數(shù)含義可通過文獻(xiàn)[18]獲得。

1.3 海洋潮汐改正

由海洋潮汐引起的重力位系數(shù)改正數(shù)可用式(5)計算[18]

(5)

1.4 極潮改正

(6)

式中，m1、m2參數(shù)含義及其具體計算過程可參考文獻(xiàn)[18]。

1.5 非潮汐質(zhì)量變化改正

非潮汐質(zhì)量變化改正包括大氣與海洋非潮汐質(zhì)量變化改正、陸地水質(zhì)量變化改正和冰雪質(zhì)量變化改正等地表質(zhì)量季節(jié)性變化的改正。

1.5.1 大氣與海洋非潮汐質(zhì)量變化改正

ESA提供的大氣質(zhì)量變化改正采用ECMWF(European Centre for Medium-range Weather Forecast)模型計算，其海洋非潮汐質(zhì)量變化改正采用OMCT(ocean model for circulation and tides)模型數(shù)據(jù)計算。由于GFZ研發(fā)的AOD1B(Atmosphere and Ocean De-aliasing level-1B)模型給出了大氣和海洋非潮汐變化模型，因此，本文直接采用AOD1B模型中提供的有關(guān)重力位系數(shù)的改正數(shù)進(jìn)行計算。

1.5.2 地表質(zhì)量季節(jié)性變化改正

地表質(zhì)量季節(jié)性變化改正主要為陸地水質(zhì)量變化改正和冰雪質(zhì)量變化改正等，而GRACE月時變重力場模型已扣除了各種潮汐及大氣和海洋非潮汐質(zhì)量變化影響，剩余的時變重力場變化信號主要來自地表質(zhì)量季節(jié)性變化信號，因此，可以采用GRACE月時變重力場模型計算地表質(zhì)量季節(jié)性變化影響的改正值[6,10]。

2 計算結(jié)果驗證與更新比較

為驗證本文時變重力場變化改正的正確性，并分析標(biāo)準(zhǔn)和背景模型更新前后時變重力場變化改正結(jié)果的差異，本文選擇了2010年6月1日的GOCE衛(wèi)星實測數(shù)據(jù)，分別計算出各項潮汐改正值以及非潮汐質(zhì)量變化改正值，并進(jìn)行有關(guān)比較與分析，所采用的標(biāo)準(zhǔn)和背景模型信息見表1，其中“ESA”表示ESA標(biāo)準(zhǔn)，“Tongji”表示本文采用的新標(biāo)準(zhǔn)。后文計算結(jié)果中，“ESA”表示ESA發(fā)布結(jié)果；“Tongji”表示本文按照新標(biāo)準(zhǔn)計算結(jié)果；“Tongji(ESA)”表示本文按照ESA標(biāo)準(zhǔn)計算結(jié)果，其中計算大氣與海洋非潮汐質(zhì)量變化改正時，采用了AOD1B RL05模型來代替ECMWF和OMCT模型，其他與ESA標(biāo)準(zhǔn)一致。

表1 標(biāo)準(zhǔn)和背景模型信息

2.1 時變重力場變化改正結(jié)果驗證

ESA和Tongji(ESA)比較結(jié)果見圖1—圖5，兩者差值標(biāo)準(zhǔn)差的統(tǒng)計結(jié)果見表2。綜合圖1—圖4以及表2可得，當(dāng)采用的標(biāo)準(zhǔn)和背景模型與ESA一致時，本文的計算結(jié)果與ESA的結(jié)果是吻合的；本文計算值與ESA結(jié)果差值的標(biāo)準(zhǔn)差比對應(yīng)潮汐改正值的最大量級低1～3個量級，說明本文的潮汐改正結(jié)果是正確的。圖5結(jié)果顯示，本文非潮汐質(zhì)量變化計算值與ESA發(fā)布結(jié)果在數(shù)值上稍有差異，產(chǎn)生差異的主要原因是所使用的背景模型不完全一致(本文計算大氣與海洋非潮汐質(zhì)量變化采用的背景模型為AOD1B RL05模型)，但兩者在趨勢和量級上具有較好的一致性，說明本文非潮汐質(zhì)量變化改正值計算結(jié)果應(yīng)該是可靠的。綜合看，本文時變重力場變化改正的結(jié)果是正確的，在所有改正項中，海潮改正的量級最大，其次為非潮汐質(zhì)量變化改正，極潮改正的量級最小。

圖1 直接潮汐改正值與ESA結(jié)果比較Fig.1 Comparisons between direct tide corrections and ESA results

圖2 固體潮汐改正值與ESA結(jié)果比較Fig.2 Comparisons between solid tide corrections and ESA results

2.2 ESA標(biāo)準(zhǔn)和背景模型更新比較

由于ESA計算重力梯度觀測值的時變重力場變化改正所采用的標(biāo)準(zhǔn)和背景模型較老，理論上，新標(biāo)準(zhǔn)和背景模型能更好地模型化時變重力場變化信號，因此本節(jié)探討標(biāo)準(zhǔn)和背景模型更新前后時變重力場變化改正值的差異(新標(biāo)準(zhǔn)和背景模型采用表1中“Tongji”項)，計算結(jié)果比較見圖6。由圖6可知，更新的標(biāo)準(zhǔn)和背景模型會影響時變重力場變化改正的結(jié)果，但更新前后兩者的時變重力場變化改正值趨勢一致，大小稍有差異但差異不大，總量在同一量級，下文將詳細(xì)分析標(biāo)準(zhǔn)和背景模型更新前后對高階靜態(tài)重力場反演的影響。

圖3 海洋潮汐改正值與ESA結(jié)果比較Fig.3 Comparisons between ocean tide corrections and ESA results

圖4 極潮改正值與ESA結(jié)果比較Fig.4 Comparisons between pole tide corrections and ESA results

表2 Tongji(ESA)與ESA結(jié)果的差值標(biāo)準(zhǔn)差

圖5 非潮汐改正值與ESA值的比較Fig.5 Comparisons between non-tidal corrections and ESA results

圖6 Tongji和ESA的時變重力場變化改正值比較Fig.6 Comparisons of temporal gravity field variations corrections between Tongji and ESA

3 時變重力場變化改正對高階靜態(tài)重力場反演的影響分析

本節(jié)設(shè)計了3種方案來分析時變重力場變化改正對利用GOCE衛(wèi)星重力梯度數(shù)據(jù)反演高階靜態(tài)重力場的影響，分別為：①無改正方案(Case 1)，不進(jìn)行任何時變重力場變化改正；②ESA方案(Case 2)，采用ESA發(fā)布的時變重力場變化改正值；③Tongji方案(Case 3)，采用新標(biāo)準(zhǔn)和背景模型計算時變重力場變化改正值，3種方案差異僅在于時變重力場變化改正值，其他數(shù)據(jù)處理方法完全一致。

考慮到低軌道GOCE衛(wèi)星對重力場信號更敏感，本文采用GOCE衛(wèi)星2012年10月—2013年4月期間新版的Level 1b重力梯度觀測數(shù)據(jù)[23]，按照上述3種方案進(jìn)行高階靜態(tài)重力場反演。首先選擇GIF48模型截斷至300階次作為本次反演的參考模型，利用頻帶范圍為5～100 mHz的IIR帶通濾波器對擾動重力梯度值進(jìn)行帶通濾波，采用直接法構(gòu)建反演220階次重力場模型的法方程；然后與180階次的Tongji-Grace02s模型[24]法方程進(jìn)行融合，反演220階次的無約束重力場模型。在數(shù)據(jù)融合過程中單位權(quán)中誤差取GRACE衛(wèi)星的軌道精度，GRACE KBR觀測值精度取驗后統(tǒng)計精度，而重力梯度分量的精度根據(jù)濾波之后的擾動重力梯度值進(jìn)行統(tǒng)計確定。3種方案解算模型和Tongji-Grace02s模型的大地水準(zhǔn)面階誤差如圖7所示，比對模型為EIGEN-6C4模型。

圖7 各方案模型與EIGEN-6C4模型比較大地水準(zhǔn)面階誤差Fig.7 Geoid degree error of different models with respect to EIGEN-6C4

圖7結(jié)果顯示，3種方案解算模型的大地水準(zhǔn)面階誤差幾乎一致，不能明顯地分析時變重力場變化改正對靜態(tài)重力場反演的影響。融合了GOCE衛(wèi)星重力梯度數(shù)據(jù)后，模型大地水準(zhǔn)面階誤差相對于Tongji-Grace02s模型在72階次后逐漸有了改善，在108階次后其改善效果更為明顯，說明相比GRACE衛(wèi)星，GOCE衛(wèi)星重力梯度數(shù)據(jù)對中高階重力位系數(shù)的貢獻(xiàn)更大，因此在中高階的位系數(shù)精度改善效果會比較明顯。為了更清晰地分析時變重力場變化對靜態(tài)重力場反演的影響，本文計算了全球1°×1°的大地水準(zhǔn)面(220階次)，其比較結(jié)果見圖8(a)、圖8(b)和圖8(c)，同時圖8(d)給出了100～220階次不同方案反演模型所對應(yīng)大地水準(zhǔn)面差異的統(tǒng)計結(jié)果。

由圖8(a)和圖8(b)可知，時變重力場變化改正對重力場反演是有影響的，其在局部區(qū)域的大地水準(zhǔn)面差異的絕對值最大達(dá)到1.4 cm，差異主要體現(xiàn)在南北兩極、阿根廷東部沿海和歐洲等地區(qū)，這些地區(qū)也是該時間段時變重力場變化信號較為明顯的區(qū)域，如果采用的數(shù)據(jù)觀測時間段不同，時變重力場變化影響會有所差異。但由圖8(c)可知，采用不同的標(biāo)準(zhǔn)和背景模型同樣也會影響高階靜態(tài)重力場反演結(jié)果，兩者在南北兩極、亞馬孫河口附近和歐洲等區(qū)域中稍有差異。由圖8(d)可知，從全球大地水準(zhǔn)面差異看，在220階次，時變重力場變化改正對靜態(tài)重力場反演的整體影響超過1 mm，影響主要體現(xiàn)在中高階次位系數(shù)，低于120階次的影響很小。綜合圖8可知：GOCE衛(wèi)星重力梯度觀測值的時變重力場變化改正是需要的，且不同標(biāo)準(zhǔn)和背景模型對靜態(tài)重力場反演的影響稍有差異。

圖8 不同模型間的全球大地水準(zhǔn)面比較Fig.8 Comparison of global geoid height between different models

為進(jìn)一步檢核時變重力場變化改正對利用GOCE衛(wèi)星重力梯度數(shù)據(jù)反演高階靜態(tài)重力場的影響，同時考慮到南北兩極可能會受到極區(qū)空白的影響，而極區(qū)空白對遠(yuǎn)離兩極的區(qū)域大地水準(zhǔn)面的精度影響較小[25]，因此本文利用收集到的歐洲區(qū)域(圖8中時變重力場影響較大的區(qū)域)約1 cm精度的France水準(zhǔn)面模型RAF18(圖8(b)中的A區(qū)域)和Great Britain水準(zhǔn)面模型OSGM15(圖8(b)中的B區(qū)域)對本文3種反演方案獲得的重力場模型進(jìn)行精度評定，精度統(tǒng)計結(jié)果見表3(其中截斷階次后采用EIGEN-6C4模型填充至2160階次)。

表3 RAF18和OSGM15水準(zhǔn)面模型評估模型的精度統(tǒng)計

表3結(jié)果表明，雖然不同方案模型間差異不大，但綜合看，進(jìn)行時變重力場變化改正(Case 2和Case 3)方案比無改正方案(Case 1)解算的模型精度稍好，說明重力梯度觀測值的時變重力場變化改正是有效的，同時Tongji方案(Case 3)精度稍優(yōu)于ESA方案(Case 2)，說明背景模型和標(biāo)準(zhǔn)的更新更有利于扣除時變重力場變化信號，從而提高高階靜態(tài)重力場反演的精度。

4 結(jié)束語

本文研究了GOCE衛(wèi)星重力梯度觀測值的時變重力場變化改正方法，更新了ESA標(biāo)準(zhǔn)和背景模型，并從重力場反演角度分析了時變重力場變化改正對高階靜態(tài)重力場反演的影響，實現(xiàn)了從GOCE衛(wèi)星Level 1b重力梯度數(shù)據(jù)直接進(jìn)行重力場反演，主要結(jié)論如下。

(1) 在時變重力場變化改正中，海洋潮汐改正值的量級最大，其次為非潮汐質(zhì)量變化改正，極潮改正值的量級最小。利用3種時變重力場變化改正方案進(jìn)行高階靜態(tài)重力場模型反演，并進(jìn)行分析與比較，結(jié)果表明：從全球大地水準(zhǔn)面比較可看出時變重力場變化改正對重力場反演結(jié)果是有影響的，其影響在局部區(qū)域的大地水準(zhǔn)面差異絕對值最大達(dá)到1.4 cm；在全球范圍內(nèi)對重力場反演的整體影響超過1 mm，因此利用GOCE衛(wèi)星重力梯度觀測值反演高階靜態(tài)重力場時需要扣除時變重力場變化的影響。

(2) 利用France水準(zhǔn)面模型RAF18和Great Britain水準(zhǔn)面模型OSGM15檢核3種方案模型的精度，結(jié)果表明進(jìn)行時變重力場變化改正方案比無改正方案解算的模型精度稍好，且更新后的標(biāo)準(zhǔn)和背景模型更有利于扣除時變重力場變化的影響。