海潮負(fù)荷地心運(yùn)動對GPS坐標(biāo)時間序列 周期信號的影響分析

郭仕偉 施 闖 魏 娜 張 濤

1 北京航空航天大學(xué)電子信息工程學(xué)院，北京市學(xué)院路37號，100191 2 武漢大學(xué)衛(wèi)星導(dǎo)航定位技術(shù)研究中心，武漢市珞喻路129號，430079

海潮負(fù)荷效應(yīng)在垂直方向上引起的測站位移可達(dá)數(shù)cm[1]，此誤差在高精度GPS數(shù)據(jù)處理中不可忽略。諸多學(xué)者分析不同海潮模型引起的測站海潮負(fù)荷位移差異[2-3]以及不同模型的海潮負(fù)荷改正對GPS精密定位[4-5]和坐標(biāo)時間序列的影響[6]。海潮負(fù)荷改正一般通過全球海潮模型和格林函數(shù)褶積積分計(jì)算得出，格林函數(shù)采用負(fù)荷勒夫數(shù)導(dǎo)出，其中一階負(fù)荷勒夫數(shù)與選用的地球參考框架原點(diǎn)有關(guān)[7]。地球參考框架的原點(diǎn)可以定義為包括大氣和海洋在內(nèi)的整個地球系統(tǒng)的質(zhì)心CM或者固體地球的質(zhì)心CE[8]。海潮負(fù)荷改正可以基于CM框架或者CE框架進(jìn)行計(jì)算，CM和CE框架下不同的一階負(fù)荷勒夫數(shù)會得到不同的一階形變，兩者間的差異即為海潮負(fù)荷引起的地心運(yùn)動。本文分別采用CM框架和CE框架的海潮負(fù)荷改正計(jì)算132個全球站2002～2019年長達(dá)18 a的GPS坐標(biāo)時間序列，重點(diǎn)分析海潮負(fù)荷引起的地心運(yùn)動對GPS坐標(biāo)時間序列的影響，尤其是對異常周期信號的影響。

1 海潮負(fù)荷改正

對于指定測站，海潮負(fù)荷位移可通過11個主潮波引起的位移進(jìn)行疊加計(jì)算[9]，即

(1)

式中，Δc為t時刻測站某一方向上(N、E、U)的海潮負(fù)荷位移；N為疊加的潮波數(shù)；ωk和χk(t0)分別為分潮波k的角頻率和天文幅角初相；fk和μk分別為分潮波k的交點(diǎn)因子和交點(diǎn)訂正角，與月球軌道升交點(diǎn)的經(jīng)度有關(guān)；Ack和φck分別為分潮波k在測站相應(yīng)方向上的振幅和相位，可根據(jù)具體的海潮模型計(jì)算得到。11個主潮波包括4個半日潮波M2、S2、N2、K2，4個全日潮波K1、O1、P1、Q1，以及3個長周期潮波Mf、Mm、Ssa。IERS采用一種更加嚴(yán)密的方法計(jì)算海潮負(fù)荷位移。首先利用Hans-Georg Scherneck網(wǎng)頁工具計(jì)算得到11個主潮波在CM和CE框架下的振幅和相位；然后利用樣條插值法將11個主潮波分量進(jìn)行內(nèi)插，得到共計(jì)342個分潮波的振幅和相位；最終計(jì)算得到總的海潮負(fù)荷位移。

2 數(shù)據(jù)處理策略

2.1 PPP坐標(biāo)解算

本文在IGS2014框架下選取132個全球站(如圖1所示，其中紅色測站用于Helmert轉(zhuǎn)換)，利用武漢大學(xué)研發(fā)的PANDA軟件，采用靜態(tài)精密單點(diǎn)定位PPP模式解算測站坐標(biāo)。測站總體觀測時段為2002-01-01～2019-12-31。精密衛(wèi)星軌道和鐘差采用德國地學(xué)研究中心GFZ提供的第二次重處理產(chǎn)品和事后產(chǎn)品(均為SP3格式)。需要注意的是，自2017-01-29(GPS 1 934周)起，GFZ軌道采用IGS2014框架，而之前的軌道產(chǎn)品均采用IGb2008框架。

圖1 132個測站分布

采用無電離層組合偽距/相位非差觀測值進(jìn)行數(shù)據(jù)處理，數(shù)據(jù)采樣率為300 s，衛(wèi)星截止高度角為7°。根據(jù)不同時期軌道產(chǎn)品的參考框架選用igs08.atx或igs14.atx絕對天線相位中心改正模型。利用雙頻無電離層組合觀測值消除電離層延遲一階項(xiàng)，忽略高階項(xiàng)。采用GMF/GPT2模型改正對流層延遲[10-11]，天頂對流層延遲每小時估計(jì)一次，將接收機(jī)鐘差作為白噪聲過程進(jìn)行估計(jì)。采用文獻(xiàn)[12]的方法，將模糊度參數(shù)固定為整數(shù)，固體潮、極潮引起的測站位移根據(jù)IERS 2010協(xié)議[9]進(jìn)行改正。

為分析海潮負(fù)荷引起的地心運(yùn)動對GPS精度定位和坐標(biāo)時間序列的影響，采用3種處理策略對海潮負(fù)荷位移進(jìn)行改正：1)不改正(NO-OTL)；2)CM框架的海潮負(fù)荷改正(CM-OTL)，在計(jì)算各潮波的振幅和相位時，需要考慮海潮負(fù)荷引起的地心運(yùn)動；3)CE框架的海潮負(fù)荷改正(CE-OTL)，在計(jì)算各潮波的振幅和相位時，不考慮海潮負(fù)荷引起的地心運(yùn)動。針對CM框架和CE框架，采用FES2004海潮模型(GFZ軌道計(jì)算時亦采用)計(jì)算11個主潮波的振幅和相位[13]。在此基礎(chǔ)上采用IERS 2010協(xié)議推薦的方法[9]，將11個主潮汐分量應(yīng)用樣條插值法擴(kuò)展為342個小潮波，以此計(jì)算總的海潮負(fù)荷位移改正。除海潮負(fù)荷位移改正的處理策略不同外，本文其他數(shù)據(jù)處理策略保持不變。

2.2 坐標(biāo)時間序列分析

為消除GFZ軌道產(chǎn)品解算PPP測站坐標(biāo)時存在的框架差異，采用Helmert轉(zhuǎn)換將52個核心站(圖1中紅色測站)的PPP天解坐標(biāo)對準(zhǔn)到IGS2014框架中，最終獲得框架一致的坐標(biāo)時間序列，用于后續(xù)的時間序列分析。

本文采用Hector軟件包進(jìn)行時間序列分析[14]。首先對原始坐標(biāo)時間序列進(jìn)行預(yù)處理，在考慮跳變、周年項(xiàng)和半周年項(xiàng)的基礎(chǔ)上，采用最小二乘擬合法得到時間序列的線性趨勢；然后采用四分位距法探測去除趨勢項(xiàng)后的坐標(biāo)殘差時間序列異常值[15]，從而去除原始坐標(biāo)時間序列中的粗差。在此基礎(chǔ)上，基于閃爍噪聲(FN)+白噪聲(WN)模型，估計(jì)時間序列中的跳變、線性趨勢、周年項(xiàng)和半周年項(xiàng)、GPS交點(diǎn)年信號[16]的第1～9個諧波。由地震、儀器變化等因素引起的跳變信息參考國際GNSS服務(wù) IGS提供的soln_IGS14.snx文件。

PPP解算的測站坐標(biāo)時間序列大多是含噪聲的非均勻采樣序列，而快速傅立葉變換等常用的頻譜分析方法要求序列均勻采樣，因此常用的頻譜分析方法在處理坐標(biāo)時間序列時存在明顯不足。Lomb-Scargle周期圖[17]是一種適用于非均勻采樣序列的頻譜分析方法，常用于GNSS 時間序列頻譜分析領(lǐng)域[8, 16,18-19]?；诖?，本文采用Lomb-Scargle周期圖對去除跳變、線性趨勢的坐標(biāo)殘差序列進(jìn)行頻譜分析，提取坐標(biāo)時間序列中的周期信號。由于坐標(biāo)時間序列中的周期信號(特別是海潮負(fù)荷效應(yīng)引起的周期信號)比較微弱，單個測站坐標(biāo)時間序列的功率譜容易受噪聲影響，因此本文首先計(jì)算單個測站坐標(biāo)殘差時間序列的功率譜，然后將所有測站的功率譜進(jìn)行堆積處理[16]，即按照N、E和U坐標(biāo)分量將不同測站的功率譜聚合成單個功率譜，按頻率排序后采用0.05 cpy的平滑窗口進(jìn)行平滑處理，最終得到堆積功率譜。

3 結(jié)果與分析

根據(jù)上述數(shù)據(jù)處理方法，采用NO-OTL、CM-OTL和CE-OTL 三種處理策略計(jì)算得到3組測站坐標(biāo)時間序列。

3.1 對測站坐標(biāo)的影響

為分析海潮負(fù)荷引起的地心運(yùn)動對測站坐標(biāo)的影響，對同一測站的CM-OTL和CE-OTL坐標(biāo)時間序列作差，以IGS2014框架下2002.0歷元的參考坐標(biāo)為基準(zhǔn)，提取N、E和U方向上的坐標(biāo)差異時間序列。圖2為波蘭的BOR1(17.07°E，52.38°N)和瑞典的ONSA(11.93°E，57.40°N)2個測站 CM-OTL與CE-OTL之間的坐標(biāo)差異時間序列，BOR1和ONSA測站分別距離海岸線約200 km和0.1 km。圖3為加拿大東南部地區(qū)NRC1(75.62°W，45.45°N)和HLFX(63.61°W，44.68°N)2個測站CM-OTL與CE-OTL之間的坐標(biāo)差異時間序列，NRC1和HLFX測站分別距離海岸線約400 km和0.3 km。圖中水平方向紅色虛線表示坐標(biāo)差異時間序列的平均值。由圖2和圖3可以看出，CM-OTL和CE-OTL之間的坐標(biāo)差異沒有明顯的系統(tǒng)性偏差，但會隨時間呈周期性變化。圖2的N方向和圖3的N、E方向均存在顯著的周年信號，圖3的N、E方向還包含明顯的長周期信號。對比圖2和圖3可以發(fā)現(xiàn)，不同地區(qū)測站坐標(biāo)差異時間序列的變化趨勢差別較大，但相同區(qū)域測站坐標(biāo)差異時間序列的變化趨勢較為一致。例如圖3中NRC1和HLFX測站在N、E方向上具有幾乎相同的周期信號。由此可知，對同一地區(qū)的測站而言，海潮負(fù)荷地心運(yùn)動引起的坐標(biāo)差異具有相似的趨勢特征，與測站至海岸線的距離無顯著相關(guān)性。

圖4為全球132個測站CM-OTL與CE-OTL之間坐標(biāo)差異時間序列的RMS。由圖可見，海潮負(fù)荷相關(guān)的地心運(yùn)動引起的水平方向上(N、E)的坐標(biāo)差異為0.4～1.4 mm，平均值為0.7 mm；垂直方向(U)上坐標(biāo)差異為0.6～2.8 mm，平均值為1.3 mm。總體來看，海潮負(fù)荷地心運(yùn)動引起的垂直方向坐標(biāo)差異約為水平方向的2倍。結(jié)合圖2、圖3、圖4可以發(fā)現(xiàn)，海潮負(fù)荷地心運(yùn)動引起的坐標(biāo)差異量級與測站至海岸線的距離無顯著相關(guān)性，但具有一定的區(qū)域分布特征。在北美洲北部、歐洲西北部和南極洲地區(qū)，測站坐標(biāo)差異略低于平均水平；而在赤道附近，測站坐標(biāo)差異略高于平均水平。

圖2 CM-OTL與CE-OTL之間的坐標(biāo)差異

圖3 CM-OTL與CE-OTL之間的坐標(biāo)差異

圖4 132個測站CM-OTL與CE-OTL之間坐標(biāo)差異時間序列的RMS

3.2 對周期信號的影響

為分析坐標(biāo)時間序列中的周期信號，采用Lomb-Scargle周期圖計(jì)算所有測站的堆積功率譜。圖5和圖6分別為不同頻率下NO-OTL、CM-OTL和CE-OTL坐標(biāo)殘差時間序列的堆積功率譜。圖5中垂直方向的實(shí)線表示周年和半周年信號，垂直方向的虛線表示GPS交點(diǎn)年信號的第1～9個諧波；圖6中垂直方向的虛線依次表示27.6 d、14.8 d、14.2 d、13.6 d、9.6 d、9.4 d、9.1 d信號。

圖5 不同海潮負(fù)荷改正策略解算的坐標(biāo)時間序列堆積功率譜(0.7～10 cpy)

圖6 不同海潮負(fù)荷改正策略解算的坐標(biāo)時間序列堆積功率譜(10～50 cpy)

由圖5和圖6可知， NO-OTL的功率譜在半周年、第2個GPS交點(diǎn)年諧波、27.6 d、14 d(14.8 d、14.2 d、13.6 d)、9 d(9.6 d、9.4 d、9.1 d)信號附近存在顯著尖峰，說明若不改正海潮負(fù)荷位移，海潮負(fù)荷效應(yīng)會引入周期信號。上述周期信號可以分為2類：1)周期信號來自海潮信號的直接傳遞，例如半周年、27.6 d、13.6 d的周期信號分別對應(yīng)長周期潮波Ssa、Mm、Mf；2)周期信號與亞周日(接近24 h或12 h)信號的混疊有關(guān)[18]，比如在24 h數(shù)據(jù)處理間隔下，M2、O1、N2、Q1的混疊周期分別為14.8 d、14.2 d、9.6 d、9.4 d；在衛(wèi)星軌道重復(fù)周期(23.93 h)采樣條件下，M2和O1的混疊周期接近13.6 d，N2和Q1的混疊周期接近9.1 d。與NO-OTL相比，CM-OTL和CE-OTL中海潮負(fù)荷效應(yīng)引起的周期信號顯著減小，其中半周年信號和第2個GPS交點(diǎn)年諧波信號輕微減弱，14 d和9 d的周期信號顯著減弱，27.6 d信號消失。

由于CM-OTL和CE-OTL采用的海潮負(fù)荷改正框架不同，因此二者計(jì)算的坐標(biāo)時間序列功率譜存在顯著差異，主要表現(xiàn)為CM-OTL坐標(biāo)時間序列的14 d和9.1 d信號更加明顯。表1(單位mm2)為采用CM-OTL和CE-OTL解算的14.8 d、14.2 d、13.6 d、9.1 d信號功率。統(tǒng)計(jì)結(jié)果表明，與CM-OTL相比，CE-OTL中14.8 d信號在N、E、U方向上分別減弱17.0%、17.1%、11.6%，14.2 d信號分別減弱47.9%、34.7%、19.7%，13.6 d信號分別減弱67.7%、83.5%、88.0%，9.1 d信號分別減弱24.0%、29.6%、55.4%?？梢钥闯觯３必?fù)荷地心改正對13.6 d、9.1 d信號的影響十分顯著。

表1 采用CM-OTL和CE-OTL解算的信號功率

為比較海潮負(fù)荷引起的地心運(yùn)動對坐標(biāo)時間序列周期信號的影響，對各測站CM-OTL與CE-OTL坐標(biāo)時間序列作差，可消除公共的地球物理信號和其他模型誤差，得到的差分時間序列將包含海潮負(fù)荷地心運(yùn)動引起的坐標(biāo)差異以及相應(yīng)的周期信號。圖7為CM-OTL與CE-OTL坐標(biāo)時間序列之差的堆積功率譜。由圖可知，坐標(biāo)差異時間序列存在明顯的周年、GPS交點(diǎn)年、27.6 d、14 d以及9 d信號，說明海潮負(fù)荷地心運(yùn)動可引入GPS交點(diǎn)年信號。由此推測，未完全模型化的海潮負(fù)荷效應(yīng)或與海潮負(fù)荷效應(yīng)周期類似的亞周日信號，可能是引入GPS交點(diǎn)年信號的原因之一。

圖7 CM-OTL與CE-OTL坐標(biāo)時間序列之差的堆積功率譜

3.3 討論

由于受到軌道動力學(xué)的約束，GPS衛(wèi)星軌道的框架原點(diǎn)指向地心CM。根據(jù)IGS規(guī)定，衛(wèi)星軌道從慣性坐標(biāo)系轉(zhuǎn)為地固坐標(biāo)系時需要進(jìn)行地心改正(center of mass correction, CMC)，此時的框架原點(diǎn)指向測站網(wǎng)的中心(center of network, CN)。由于IGS及其分析中心提供的SP3軌道已經(jīng)考慮了海潮負(fù)荷引起的地心改正[20]，因此就海潮負(fù)荷而言，其框架原點(diǎn)指向CN。同樣，采用SP3軌道解算的PPP坐標(biāo)屬于CN框架，CN可近似為CE，但是CN與CM之間的差異為CMC[21]。因此，CN框架的SP3軌道與CE框架的海潮負(fù)荷改正具有更好的框架一致性。本文研究結(jié)果表明，采用CN框架(針對海潮負(fù)荷)的SP3軌道和采用CM框架的海潮負(fù)荷改正將引入GPS交點(diǎn)年信號和14 d信號。利用IGS及其分析中心提供的精密軌道進(jìn)行PPP解算時，推薦采用CE框架進(jìn)行海潮負(fù)荷改正。

4 結(jié) 語

1)在全球范圍內(nèi)，海潮負(fù)荷地心運(yùn)動引起的坐標(biāo)差異在N、E、U方向上的平均值分別為0.7 mm、0.7 mm和1.3 mm，垂直方向的坐標(biāo)差異約為水平方向的2倍。海潮負(fù)荷地心運(yùn)動引起的坐標(biāo)差異具有顯著的區(qū)域分布特征，鄰近測站的坐標(biāo)差異序列具有相似的變化趨勢。

2)當(dāng)海潮負(fù)荷改正所屬框架與GPS軌道所屬框架不一致時，海潮負(fù)荷效應(yīng)可混疊產(chǎn)生周期為GPS交點(diǎn)年信號、14 d、9 d的異常信號，與海潮負(fù)荷周期接近的亞周日誤差可能是引入GPS交點(diǎn)年信號的原因之一。本文使用的GPS軌道考慮了海潮負(fù)荷效應(yīng)引起的地心運(yùn)動，與CM-OTL相比，CE-OTL中海潮負(fù)荷效應(yīng)引入的周期信號顯著減弱，其中14.8 d信號減弱11.6%～17.1%、14.2 d信號減弱19.7%～47.9%、13.6 d信號減弱67.7%～88.0%、9.1 d信號幾乎消失。

3)IGS及其分析中心提供的SP3軌道從慣性坐標(biāo)系轉(zhuǎn)為地固坐標(biāo)系時需要進(jìn)行地心改正。對采用SP3軌道解算的PPP坐標(biāo)而言，其框架原點(diǎn)近似為固體地球質(zhì)心CE。采用IGS提供的SP3軌道計(jì)算PPP時，應(yīng)采用CE框架進(jìn)行海潮負(fù)荷改正。