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

1. 廣東工業(yè)大學(xué)測(cè)繪工程系，廣東 廣州 510006； 2. 同濟(jì)大學(xué)測(cè)繪與地理信息學(xué)院，上海 200092； 3. 波恩大學(xué)大地測(cè)量與地理信息學(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é)、大地測(cè)量學(xué)和冰川學(xué)等科學(xué)研究提供高精度的中高階靜態(tài)地球重力場(chǎng)模型[1]。GOCE衛(wèi)星在任務(wù)期間共獲取了約42個(gè)月的有效重力梯度觀測(cè)值，由于重力梯度觀測(cè)值包含了時(shí)變重力場(chǎng)變化信號(hào)，在高階靜態(tài)重力場(chǎng)反演中需考慮時(shí)變重力場(chǎng)變化改正的影響。ESA在Level 2級(jí)別的GOCE衛(wèi)星重力梯度數(shù)據(jù)中給出了ESA計(jì)算的時(shí)變重力場(chǎng)變化改正值[2]，目前國(guó)內(nèi)外絕大多數(shù)科研機(jī)構(gòu)均直接利用GOCE衛(wèi)星Level 2重力梯度數(shù)據(jù)進(jìn)行重力場(chǎng)反演[3-5]，但時(shí)變重力場(chǎng)變化改正是實(shí)現(xiàn)由GOCE衛(wèi)星Level 1b級(jí)別數(shù)據(jù)到Level 2級(jí)別數(shù)據(jù)轉(zhuǎn)換的重要過(guò)程，也是完整GOCE衛(wèi)星數(shù)據(jù)處理過(guò)程中的一個(gè)重要環(huán)節(jié)，而ESA有關(guān)文檔中并未給出詳細(xì)的數(shù)據(jù)處理過(guò)程，自主實(shí)現(xiàn)由GOCE衛(wèi)星Level 1b重力梯度數(shù)據(jù)直接進(jìn)行重力場(chǎng)反演具有重要的科學(xué)意義，可為我國(guó)獨(dú)立開(kāi)展相關(guān)GOCE類衛(wèi)星任務(wù)的數(shù)據(jù)處理研究提供重要參考。

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

1 原理與方法

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

1.1 直接潮汐改正

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

(1)

1.2 固體潮汐改正

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

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

(2)

(3)

第2步：計(jì)算2階項(xiàng)的日潮、半日潮及長(zhǎng)周期潮改正，其中計(jì)算日潮改正的公式為

(4)

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

1.3 海洋潮汐改正

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

(5)

1.4 極潮改正

(6)

式中，m1、m2參數(shù)含義及其具體計(jì)算過(guò)程可參考文獻(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)模型計(jì)算，其海洋非潮汐質(zhì)量變化改正采用OMCT(ocean model for circulation and tides)模型數(shù)據(jù)計(jì)算。由于GFZ研發(fā)的AOD1B(Atmosphere and Ocean De-aliasing level-1B)模型給出了大氣和海洋非潮汐變化模型，因此，本文直接采用AOD1B模型中提供的有關(guān)重力位系數(shù)的改正數(shù)進(jìn)行計(jì)算。

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

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

2 計(jì)算結(jié)果驗(yàn)證與更新比較

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

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

2.1 時(shí)變重力場(chǎng)變化改正結(jié)果驗(yàn)證

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

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

圖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的時(shí)變重力場(chǎng)變化改正值比較Fig.6 Comparisons of temporal gravity field variations corrections between Tongji and ESA

3 時(shí)變重力場(chǎng)變化改正對(duì)高階靜態(tài)重力場(chǎng)反演的影響分析

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

考慮到低軌道GOCE衛(wèi)星對(duì)重力場(chǎng)信號(hào)更敏感，本文采用GOCE衛(wèi)星2012年10月—2013年4月期間新版的Level 1b重力梯度觀測(cè)數(shù)據(jù)[23]，按照上述3種方案進(jìn)行高階靜態(tài)重力場(chǎng)反演。首先選擇GIF48模型截?cái)嘀?00階次作為本次反演的參考模型，利用頻帶范圍為5～100 mHz的IIR帶通濾波器對(duì)擾動(dòng)重力梯度值進(jìn)行帶通濾波，采用直接法構(gòu)建反演220階次重力場(chǎng)模型的法方程；然后與180階次的Tongji-Grace02s模型[24]法方程進(jìn)行融合，反演220階次的無(wú)約束重力場(chǎng)模型。在數(shù)據(jù)融合過(guò)程中單位權(quán)中誤差取GRACE衛(wèi)星的軌道精度，GRACE KBR觀測(cè)值精度取驗(yàn)后統(tǒng)計(jì)精度，而重力梯度分量的精度根據(jù)濾波之后的擾動(dòng)重力梯度值進(jìn)行統(tǒng)計(jì)確定。3種方案解算模型和Tongji-Grace02s模型的大地水準(zhǔn)面階誤差如圖7所示，比對(duì)模型為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)面階誤差幾乎一致，不能明顯地分析時(shí)變重力場(chǎng)變化改正對(duì)靜態(tài)重力場(chǎng)反演的影響。融合了GOCE衛(wèi)星重力梯度數(shù)據(jù)后，模型大地水準(zhǔn)面階誤差相對(duì)于Tongji-Grace02s模型在72階次后逐漸有了改善，在108階次后其改善效果更為明顯，說(shuō)明相比GRACE衛(wèi)星，GOCE衛(wèi)星重力梯度數(shù)據(jù)對(duì)中高階重力位系數(shù)的貢獻(xiàn)更大，因此在中高階的位系數(shù)精度改善效果會(huì)比較明顯。為了更清晰地分析時(shí)變重力場(chǎng)變化對(duì)靜態(tài)重力場(chǎng)反演的影響，本文計(jì)算了全球1°×1°的大地水準(zhǔn)面(220階次)，其比較結(jié)果見(jiàn)圖8(a)、圖8(b)和圖8(c)，同時(shí)圖8(d)給出了100～220階次不同方案反演模型所對(duì)應(yīng)大地水準(zhǔn)面差異的統(tǒng)計(jì)結(jié)果。

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

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

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

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

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

4 結(jié)束語(yǔ)

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

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

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