海潮負荷地心運動對GPS坐標時間序列 周期信號的影響分析

郭仕偉 施 闖 魏 娜 張 濤

1 北京航空航天大學電子信息工程學院，北京市學院路37號，100191 2 武漢大學衛星導航定位技術研究中心，武漢市珞喻路129號，430079

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

1 海潮負荷改正

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

(1)

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

2 數據處理策略

2.1 PPP坐標解算

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

圖1 132個測站分布

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

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

2.2 坐標時間序列分析

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

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

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

3 結果與分析

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

3.1 對測站坐標的影響

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

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

圖2 CM-OTL與CE-OTL之間的坐標差異

圖3 CM-OTL與CE-OTL之間的坐標差異

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

3.2 對周期信號的影響

為分析坐標時間序列中的周期信號，采用Lomb-Scargle周期圖計算所有測站的堆積功率譜。圖5和圖6分別為不同頻率下NO-OTL、CM-OTL和CE-OTL坐標殘差時間序列的堆積功率譜。圖5中垂直方向的實線表示周年和半周年信號，垂直方向的虛線表示GPS交點年信號的第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 不同海潮負荷改正策略解算的坐標時間序列堆積功率譜(0.7～10 cpy)

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

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

由于CM-OTL和CE-OTL采用的海潮負荷改正框架不同，因此二者計算的坐標時間序列功率譜存在顯著差異，主要表現為CM-OTL坐標時間序列的14 d和9.1 d信號更加明顯。表1(單位mm2)為采用CM-OTL和CE-OTL解算的14.8 d、14.2 d、13.6 d、9.1 d信號功率。統計結果表明，與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%。可以看出，海潮負荷地心改正對13.6 d、9.1 d信號的影響十分顯著。

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

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

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

3.3 討論

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

4 結 語

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

2)當海潮負荷改正所屬框架與GPS軌道所屬框架不一致時，海潮負荷效應可混疊產生周期為GPS交點年信號、14 d、9 d的異常信號，與海潮負荷周期接近的亞周日誤差可能是引入GPS交點年信號的原因之一。本文使用的GPS軌道考慮了海潮負荷效應引起的地心運動，與CM-OTL相比，CE-OTL中海潮負荷效應引入的周期信號顯著減弱，其中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軌道從慣性坐標系轉為地固坐標系時需要進行地心改正。對采用SP3軌道解算的PPP坐標而言，其框架原點近似為固體地球質心CE。采用IGS提供的SP3軌道計算PPP時，應采用CE框架進行海潮負荷改正。