邵 璠124，王小亞1234

(1.中國科學(xué)院 上海天文臺，上海 200030；2.中國科學(xué)院大學(xué)，北京 100049；3.中國科學(xué)院 上海天文臺上海市空間導(dǎo)航與定位技術(shù)重點實驗室，上海 200030；4.宇航動力學(xué)國家重點實驗室，西安 710043)

1 引言

與地球相關(guān)的參考架原點主要有3個：包括海洋和大氣的整個地球的質(zhì)量中心(center of mass,CM)、不包括質(zhì)量負荷的固體地球的質(zhì)量中心(center of Earth,CE)和固體地球外表面的形狀中心(center of figure,CF)。CM通常被空間大地測量所采用，因為衛(wèi)星動力學(xué)是針對CM的，理想地球參考架定義要求地球參考架原點定義在CM。CE應(yīng)用于某些地球物理學(xué)理論研究，如勒夫數(shù)的計算等。CF通常應(yīng)用于與地面測量相關(guān)的學(xué)科，這里地面點間的幾何形狀是唯一可測量的量，所以一般應(yīng)用于與地面大地測量和形變測量等有關(guān)的研究[1,2]。國際地球自轉(zhuǎn)服務(wù)(International Earth Rotation and Reference System Service,IERS)發(fā)布的國際地球參考架(International Terrestrial Reference Frame,ITRF)是通過SLR、甚長基線干涉測量(Very Long Baseline Interferometry,VLBI)、全球定位系統(tǒng)(Global Positioning System,GPS)和星載多普勒定軌定位系統(tǒng)(Doppler orbitography and radio-positioning integrated by satellite,DORIS)四種空間大地測量技術(shù)建立和維持的，由各技術(shù)的全球觀測網(wǎng)數(shù)據(jù)得到的技術(shù)內(nèi)綜合解，經(jīng)由技術(shù)間再綜合得到臺站坐標(biāo)和速度場來實現(xiàn)，其最初目標(biāo)是要實現(xiàn)理想的地球參考架。由于四種空間大地測量技術(shù)臺站都固定在地殼上，由這些臺站坐標(biāo)和速度場確定的ITRF并不能直接反映地球整體質(zhì)量分布和地心運動，而是反映了由這些臺站構(gòu)成的多面體中心，即觀測站網(wǎng)的中心(center of network,CN)，其近似于地球形狀中心CF[1?4]。由于地球上海平面變化、冰川融化、大氣環(huán)流、地幔對流和液核振蕩等原因，地球質(zhì)量會發(fā)生遷移，導(dǎo)致地球質(zhì)心CM相對于地球形狀中心(即參考架原點)發(fā)生位移，這個位移被稱為地心運動[2]。地心運動必然會影響衛(wèi)星精密定軌、地面點的精密定位、地球定向參數(shù)的確定、地球低階引力場的確定以及所有以ITRF作為參考架的空間大地測量結(jié)果。因此，對地心運動的研究十分重要。

Cheng等人[5]利用5顆激光衛(wèi)星的SLR數(shù)據(jù)，按照IERS2003規(guī)范，通過估算地球引力位一階系數(shù)(采用的參考架為LPOD2005，在SLRF2005的基礎(chǔ)上更新了一些測站坐標(biāo)和速度)，研究了地心的周年運動。他們與GRACE得到的結(jié)果進行了比較，發(fā)現(xiàn)其X方向(2.5 mm)與ILRSX方向(2.7 mm)和GRACEX方向(3 mm)符合得很好，而Y和Z方向與ILRS和GRACE結(jié)果符合較差。Collilieux等人[6]以全球SLR數(shù)據(jù)解出的測站坐標(biāo)作為類觀測，進行參考架轉(zhuǎn)換，得到平移參數(shù)，但是，認為該平移參數(shù)因為測站分布的不均勻性和稀疏性而不能被嚴格地當(dāng)做地心運動，它們之間的差別定義為網(wǎng)形效應(yīng)。國內(nèi)也有一些學(xué)者對地心運動進行了研究，但地心運動的振幅與周期特性不完全一致。吳斌等人[7]認為地心變化幅度不超過5 cm，秦顯平和楊元喜[8]發(fā)現(xiàn)地心運動存在近周年和半年的周期分量，但地心運動的長期變化并不明顯；朱元蘭和馮初剛[9]研究表明地心X、Y方向變化幅度為1 cm，Z方向為3～4 cm；郭金運等人[10]發(fā)現(xiàn)地心運動存在長期和周期變化；李九龍[11]研究表明地心運動主要呈現(xiàn)出季節(jié)性周期變化，其中以周年項變化為主，半周年項的變化不夠明顯。究其原因，很大可能是基于老的IERS規(guī)范和針對單一衛(wèi)星進行解算得到的地心運動并不可靠。我們利用上海天文臺SHORD II軟件，處理了Lageos-1,Lageos-2,Etalon-1和Etalon-2共4顆衛(wèi)星從1993年1月―2017年12月的25年SLR數(shù)據(jù)，得到上海天文臺(Shanghai Astronomical Observatory,SHAO)周解；并與ILRS各分析中心ACs的周解進行技術(shù)內(nèi)綜合，獲得幾何法的地心運動序列；又利用SHORD II軟件，直接求解地心運動序列，并將這兩類解的結(jié)果與CSR動力法結(jié)果進行比較；最后分析和探討了不同地心運動確定方法結(jié)果的特征、可靠性和影響因素。

2 地心運動的計算方法

目前主要有四種空間大地測量技術(shù)，分別為VLBI,SLR,GPS和DORIS。其中，VLBI技術(shù)是一種相對測量，其觀測目標(biāo)是遙遠的河外射電源，它們并不圍繞地球轉(zhuǎn)動，且其數(shù)據(jù)處理采用的是幾何方法，所以VLBI技術(shù)自身無法確定地球參考架的原點。GPS技術(shù)由于受軌道模型誤差及其他因素如衛(wèi)星鐘差、整周模糊度、天線相位中心改正等相關(guān)性的影響，所以無法給出精確的地心運動變化，同時，GPS衛(wèi)星軌道較高且其數(shù)據(jù)處理經(jīng)常采用的差分法也降低了對地心變化的敏感性。DORIS技術(shù)由于受自身系統(tǒng)誤差較嚴重和在軌觀測衛(wèi)星數(shù)變化影響，不能精確地解算出地球質(zhì)心的位置[11]。SLR技術(shù)是目前最精密的地球參考架原點確定技術(shù)，其觀測的目標(biāo)大都是面質(zhì)比小、結(jié)構(gòu)簡單的地球動力學(xué)衛(wèi)星，且大部分激光衛(wèi)星離地球比較近，對地心運動很敏感，這為監(jiān)測地心運動提供了有利條件。

目前，采用SLR技術(shù)確定地心運動的方法有三種：幾何法、直接法和動力學(xué)法。

2.1 幾何法

幾何法就是先求解測站坐標(biāo)序列，然后再求解地心運動序列[8]。該方法先利用衛(wèi)星精密定軌得到測站在地心坐標(biāo)系下的坐標(biāo)序列，然后選擇國際地球參考框架ITRF(或其他較權(quán)威的地球參考框架)為參考基準，由HELMERT 7參數(shù)轉(zhuǎn)換到這個ITRF地球參考框架，得出測站坐標(biāo)相對于該ITRF地球參考架的平移量。該平移量即為地心運動時間序列。

其中，T1,T2,T3是3個平移參數(shù)，也是所要求取的地心運動3個分量；R1,R2,R3為3個旋轉(zhuǎn)參數(shù)；D為尺度因子。

幾何法受到網(wǎng)形的影響，因此，為了使計算的地心運動序列更加真實可靠，我們應(yīng)均勻選取測站。我們在Pavlis和Luceri[12]給出的SLR核心站基礎(chǔ)上，選擇性地剔除數(shù)據(jù)不好的站點(轉(zhuǎn)換到新的參考架后，與框架結(jié)果的差大于0.3 m)作為SLR核心站，見表1。通過這些核心站進行7參數(shù)轉(zhuǎn)換，得到的其中3個平移參數(shù)即是地心運動序列的3個分量。由于該方法與核心站的選取、數(shù)目和測站分布有關(guān)，所以其可靠性會受這些因素的影響。

表1 SLR核心站列表

我們采用了幾何法來確定地心運動序列。我們先通過對Lageos-1,Lageos-2,Etalon-1和Etalon-2這4顆衛(wèi)星進行弧長為7 d的松約束(站坐標(biāo)先驗標(biāo)準差為1 m，極移先驗標(biāo)準差為20 mas，LOD先驗標(biāo)準差為2 ms)定軌，得到SHAO周解，然后對SHAO周解與ILRS各分析中心的周解進行加權(quán)綜合，得到測站坐標(biāo)，具體綜合過程如圖1所示。對于SINEX解中存在的松約束不予消除，我們進行直接綜合，然后利用方差分量估算法對其進行加權(quán)綜合，這樣得到的綜合周解還是松約束解[13,14]。我們再對混合后的周解相對于SLRF2008和SLRF2014分別進行7參數(shù)轉(zhuǎn)換求平移參數(shù)，得到了兩個地心運動序列。為了在下文圖表中表述得清晰簡潔，我們將其分別命名為SLRF2008幾何法和SLRF2014幾何法。

圖1 SLR技術(shù)內(nèi)綜合流程圖

2.2 動力學(xué)法

動力學(xué)法是通過處理SLR數(shù)據(jù)，利用動力學(xué)法進行衛(wèi)星精密定軌來估算地球引力位一階系數(shù)，從而確定地心運動序列[9]。地球引力場是在地固坐標(biāo)系中以球諧函數(shù)展開的，如果坐標(biāo)原點就是地球質(zhì)心，那么引力位的一階項系數(shù)C11,S11,C10將全部為0。但實際上，地球質(zhì)心一直在運動，也就是說C11,S11,C10并不為0，動力學(xué)法就是在衛(wèi)星精密定軌過程中，將地球引力位一階系數(shù)作為待估參數(shù)，通過解算引力位一階系數(shù)來確定地球的質(zhì)心運動。地球引力位的球諧系數(shù)的級數(shù)展開形式一般為：

式中，G為地球引力常數(shù)，M為地球質(zhì)量，r為地心向徑，R為地球半徑，?為地心緯度，λ為地心經(jīng)度，Pnm為締合勒讓德函數(shù)，Cnm和Snm為位系數(shù)，n和m為位系數(shù)的階和次，Cnm和Snm的大小反映了地球內(nèi)部的物質(zhì)分布情況，其定義為：

式中，當(dāng)m=0時δnm=1，否則δnm=0；dM表示坐標(biāo)為(r,?,λ)的空間點處的質(zhì)量元。

在大多數(shù)情況下，大氣、海洋、陸地水等都近似處理為地球表面質(zhì)量負荷變化，這樣式(4)和式(5)的時變部分可以簡化為[11]：

式中，L(?,λ,t)為地球表面上隨時間t變化的一點，dS為該處的質(zhì)量元，式(3)中沒有一次項，這是假定CM在地球參考架原點，當(dāng)存在地心運動時，其定義為：

式中，x,y,z為地球表面一點的三維坐標(biāo)。如只考慮地球表面質(zhì)量負荷L(?,λ,t)變化引起的地球質(zhì)心變化，則上述方程可重寫為：

式中，?x,?y,?z為地球質(zhì)心運動的3個分量，C11,S11,C10為地球引力位一階系數(shù)。

2.3 直接法

直接法就是在衛(wèi)星精密定軌過程中將地心運動作為待估參數(shù)，在測量模型中加入地心運動改正項，通過循環(huán)迭代求得地心的位移量[9]。臺站坐標(biāo)定義在形心坐標(biāo)系下，衛(wèi)星軌道定義在質(zhì)心坐標(biāo)下，由于地球質(zhì)心相對于形狀中心存在地心運動，因此激光測距觀測值須增加這一項改正。地心運動的具體改正方程為：

該方法的一個本質(zhì)缺陷就是其地心運動估值與SLR測量誤差中的偏心改正強相關(guān)，SLR測站偏心改正公式為：

其中，?ρRO為測站偏心改正；ρ為衛(wèi)星到測站的距離，ρx,ρy,ρz分別為ρ在測站坐標(biāo)系中的3個分量；ECE為測量設(shè)備的測量中心相對于測站標(biāo)定坐標(biāo)在東向的偏心量；ECN為測量設(shè)備的測量中心相對于測站標(biāo)定坐標(biāo)在北向的偏心量；ECU為測量設(shè)備的測量中心相對于測站標(biāo)定坐標(biāo)在高程方向的偏心量。

在SLR數(shù)據(jù)處理中，通常不計算測站偏心改正，從式(11)和(12)可以看出，其公式非常相似，即地心改正實際估算了所有測站的共同偏心改正，測站的偏心改正誤差也在其中。

本文分別利用Lageos-1和Lageos-2的SLR數(shù)據(jù)，采用直接法進行了地心運動解算，為了表述清晰簡潔，其結(jié)果分別命名為Lageos-1直接法和Lageos-2直接法。直接法計算地心運動所采用的模型和解算策略如表2所示。

表2 SLR直接法計算地心運動所采用的模型和策略

3 地心運動特征分析

對獲得的地心序列首先進行了傅里葉頻譜分析，分析其中所含周期，然后假設(shè)地心運動同時具有長期項和周期項，則地心運動的時間序列可以表示為：

其中，c為長期變化中的常數(shù)項，d為長期變化中的系數(shù)項，k為周期函數(shù)的個數(shù)，Ah為對應(yīng)第h個周期項的振幅，fh為對應(yīng)的頻率，θh為對應(yīng)相位。將信號中的周期項展開可得：

將求解的地心運動?Rj當(dāng)作已知觀測量，代入式(13)中，則可得到誤差方程：

式中，V為地心運動?Rj的殘差向量，X為未知數(shù)向量(即所求的周期項和長期項系數(shù))，B為系數(shù)矩陣，L為由?Rj組成的向量。通過最小二乘法可得未知數(shù)向量的估值為：

其中，P為地心運動的權(quán)矩陣，此處取單位權(quán)。當(dāng)求得參數(shù)向量X后，可通過式(15)和式(16)得到周期項的振幅和相位，這樣就可以獲得地心運動的長期變化和周期變化特征。

4 結(jié)果分析

本文利用上海天文臺SHORD II軟件處理了Lageos-1,Lageos-2,Etalon-1和Etalon-2這4顆衛(wèi)星從1993年1月―2017年12月期間的SLR數(shù)據(jù)，將得到的SHAO周解與ILRS各分析中心(包括ASI,BKG,DGFI,ESA,GRGS,JCET和NSGF)提供的周解進行加權(quán)混合，由混合后的解進行7參數(shù)轉(zhuǎn)換，分別得到以SLRF2008和SLRF2014作為參考的兩個地心運動序列，并進行了傅里葉變換及長期和周期性特征分析，其結(jié)果見圖2―4。圖2顯示了SLRF2008幾何法確定的地心序列傅里葉頻譜分析，證實了地心在X,Y,Z三個方向上有明顯的周年項和半年項，周期分別為0.9996 a和0.4996 a，然后對地心序列中的周年項、半年項的振幅和相位以及長期項按照第3章方法進行最小二乘法擬合，這里只給出了基于SLRF2008的幾何法傅里葉頻譜分析結(jié)果，其他結(jié)果類似，并且周期也相同。圖3和圖4分別給出了SLRF2008幾何法和SLRF2014幾何法確定的地心運動及其特征擬合序列，從圖中可以看出周期擬合序列較好地模擬了原始的地心運動序列，但是在振幅上的擬合結(jié)果還不是很好，特別是在Z方向。

利用Lageos-1和Lageos-2的SLR數(shù)據(jù)并通過直接法確定了地心運動，結(jié)果如圖5和圖6所示。圖7給出了CSR利用動力學(xué)法確定的地心運動。圖中，黑色曲線代表最小二乘法擬合的線性項和周期項之和?？梢园l(fā)現(xiàn)，地心運動在X,Y方向變化幅度均比較小，而在Z方向變化幅度較大，且通過直接法求得的地心序列比幾何法和動力學(xué)法求得的變化幅度更大，更不規(guī)則，特別是在Z方向。這可能是因為直接法求解的質(zhì)心運動與SLR測站偏心改正強相關(guān)，其可靠性相對較低。

圖2 SLRF2008幾何法確定的地心運動序列傅里葉頻譜分析圖

圖3 SLRF2008幾何法確定的地心運動序列及其特征擬合序列

圖4 SLRF2014幾何法確定的地心運動序列及其特征擬合序列

圖5 Lageos-1直接法確定的地心運動序列

圖6 Lageos-2直接法確定的地心運動序列

圖7 CSR動力學(xué)法確定的地心運動序列

我們采用幾何法和直接法計算了地心運動與CSR提供的動力學(xué)法地心運動序列的平均統(tǒng)計結(jié)果(如表3所示)，其中CSR的結(jié)果來自于ftp.csr.utexas.edu/pub/slr/geocenter發(fā)布的地心運動序列，其時間跨度為2002―2017年，解算頻率為每30 d提供一組解，選擇的參考架原點為SLRF2014原點。從表3中可以看出，利用直接法求得的地心運動平均值比幾何法和動力學(xué)法的結(jié)果偏大，幾何法求得的地心運動與CSR提供的動力法結(jié)果較接近。

表3 SLR不同方法監(jiān)測地心運動的平均統(tǒng)計結(jié)果 mm

表4給出了本文和基于CSR提供的動力學(xué)地心運動序列分別擬合得到的地心運動長期變化率。從表中可以看出，利用參考SLRF2014的幾何法確定的地心長期運動速率比利用參考SLRF2008幾何法確定的結(jié)果明顯小，特別是Z分量，說明地球參考架SLRF2014相對于SLRF2008定義的原點更加穩(wěn)定。利用幾何法和直接法求得的地心運動速率Z分量與CSR結(jié)果存在差異，其原因可能是幾何法求得結(jié)果受網(wǎng)形影響，由于核心站南北分布不均勻，導(dǎo)致解算的地心運動Z分量可靠性降低。而直接法使用的是單星解，衛(wèi)星觀測幾何強度不夠，解算的地心運動Z分量存在誤差，但從整體結(jié)果來看，本文計算的地心運動速率與CSR提供的速率較一致，都比較小，說明SLR確定的地心比較穩(wěn)定，其長期變化并不明顯。其中，基于SLRF2014幾何法確定的地心穩(wěn)定性好于0.1 mm·a?1，其他方法確定的地心穩(wěn)定性差于 0.1 mm·a?1。

表 4 地心運動長期變化率 mm·a?1

本文對地心運動在X,Y,Z三個方向上的周年及半年項振幅和相位(相位的參考歷元為1月1日)進行了解算，并與CSR的結(jié)果進行對比(如表5所示)，CSR結(jié)果來自對其提供的地心運動序列進行類似的周年項及半年項最小二乘法擬合。從表5可以看出，使用Lageos-1直接法與Lageos-2直接法確定的地心運動序列，其半年項比較接近，但是周年項相位相差較大，這一方面可能與該方法參數(shù)之間的相關(guān)性較強有關(guān)，也可能與單顆衛(wèi)星解算地心運動時的觀測幾何強度不夠有關(guān)。通過比較還發(fā)現(xiàn)SLRF2014幾何法與SLRF2008幾何法得到的地心運動較接近，其中SLRF2014幾何法得到的周年項和半年項與CSR動力學(xué)法結(jié)果更接近，這說明通過對SHAO提供的SLR周解與ILRS各分析中心提供的周解進行技術(shù)內(nèi)加權(quán)綜合后，可以得到較可靠的地心運動序列。

表5 地心運動周年項的振幅和相位與半年項的比較

5 結(jié)論

本文處理了Lageos-1,Lageos-2,Etalon-1和Etalon-2這4顆衛(wèi)星在1993―2017年期間的全球SLR觀測數(shù)據(jù)，分別通過幾何法(參考SLRF2008和SLRF2014)和直接法(針對Lageos-1和Lageos-2)確定了SLR地心運動序列，利用傅里葉變換和最小二乘法分析了兩種方法解算的地心運動的長期變化率、周年項和半年項振幅和相位，并與CSR提供的動力法結(jié)果進行了比較，發(fā)現(xiàn)：(1)地心運動的長期項比周期項小1個量級；(2)與SLRF2008幾何法相比，SLRF2014幾何法得到的地心運動長期項明顯減小，特別是在Z方向，且其得到的地心運動周年項和半年項與CSR動力學(xué)法得到的結(jié)果很相似，說明綜合SHAO與ILRS ACs周解來確定地心運動是可行的；(3)將Lageos-1直接法與Lageos-2直接法確定的地心運動序列相比，其半年項有較好的一致性，但周年項相位相差較大，與幾何法及CSR提供的動力學(xué)結(jié)果在Z方向差異也很大，這可能與該方法解算的地心運動與測站偏心改正強相關(guān)有關(guān)，也可能與單顆衛(wèi)星解算地心運動幾何結(jié)構(gòu)不夠好有關(guān)，該方法的可靠性有待進一步提高。