利用PPP反演海潮負荷位移參數

鄭 凱 馬 蘭 季 偉 郭博峰

1 武漢大學測繪學院，武漢市珞喻路129號，430079

海潮負荷效應會引起測站的位移［1－3］。目前，對海潮負荷引起的地表位移改正的主要方法是利用全球海潮模型計算，但由于利用衛星測高等資料建立的海潮模型在一些地區（特別是高緯度地區）模型改正精度不高，一些學者提出利用GPS技術測定海潮負荷位移參數的研究思路。Zhang等［4］用PPP技術獲得了南極Amery冰架地區的海潮波形，采用5d連續觀測的GPS數據反演得到8個分潮在垂直分量上的振幅和相位；Thomas等［5］比較了GPS 和VLBI技術測定海潮負荷位移參數的差異，表明利用GPS和VLBI技術測定海潮負荷位移參數是可行的，且PPP固定解計算得到的參數并沒有較浮點解得到的參數在精度上有明顯改善；Yuan等［6］利用香港12個測站3～7 a的GPS觀測資料，采用與Thomas等類似的海潮負荷位移參數估算策略，將分潮參數與PPP模型其余參數一起估計。該方法雖然理論上比較嚴密，但增加了PPP 數據處理的復雜度，過度的參數化還可能造成分潮正余弦系數和模糊度參數及鐘差參數之間的耦合相關性，也在一定程度上影響PPP定位的收斂速度。本文提出直接利用未顧及海潮負荷位移改正的動態PPP 解算得到的長時間（1a）測站坐標序列，采用傅里葉變換方法提取海潮周期信息，進而反演8個分潮的海潮負荷位移參數的方法，將反演結果與最新的5個全球海潮模型比較，分析評價該方法反演海潮負荷位移參數的可行性。

1 海潮參數反演方法

測站的海潮負荷位移可由11個主分潮的海潮負荷位移疊加獲得，其中3 個為長周期分潮。由于反演長周期分潮需要長時間的觀測數據，且長周期潮對測站位移的影響較其他分潮要小［7－8］，故在此僅研究8 個主分潮，即4 個半日潮（M2、S2、N2、K2）和4個周日潮（K1、O1、P1、Q1）。

測站的海潮負荷位移可表示為：

式中，Hk，i和φk，i分別為分潮i在k方向上的振幅和格林尼治相位；t為格林尼治時間；Vi為天文幅角初相（本文天文變量參考J2000.0系統）；ωi為分潮角速度，可查表獲得；fi和μi分別為節點因數和天文相角，可由天文變量計算得到；ξk為海潮負荷引起的測站位移偏差［9］。為便于求解Hk，i和φk，i，將式（1）改為：

式中，A＝fiHk，icos（φk，i－μi），B＝fiHk，isin（φk，i－μi）。由式（2）得：

式中，A、B為待估參數，ξk為PPP計算得到的E、N、U3個方向的測站坐標序列。將式（2）以法方程形式表達，用經典最小二乘準則求解待估參數。考慮到φ變化范圍從－PI到PI，實際解算結果需要加減180°。

2 PPP處理策略

海潮負荷對單點定位的影響在沿海地區可達mm級甚至cm級［1，4］，因此在對其他誤差精確改正之后，得到的測站位移即由海潮負荷引起。

本文采用無電離層組合的PPP模型［10－11］：其中，P為無電離層的偽距觀測值，Φ為無電離層載波相位觀測值，ρ為真實幾何距離，dt為衛星鐘差，dT為接收機鐘差，Tr為對流層延遲，N為無電離層組合模糊度，λ為組合后的載波波長，ξ為系統誤差（相對論效應、固體潮改正、地球自轉改正、天線相位中心偏差及變化），ξoc為海潮負荷引起的位移偏差，εp和εΦ分別為偽距觀測噪聲和相位觀測噪聲。解算過程中，將測站坐標、接收機鐘差、無電離層組合模糊度及對流層天頂延遲視為未知數，采用Kalman濾波作為參數估計器，衛星高度角設為7°，用TurboEidt方法探測周跳（不作修復）。同時，對接收機鐘跳進行探測修復［12－13］，采用絕對天線相位中心模型（IGS＿08），用選權迭代的驗后殘差分析方法進行質量控制，并將模糊度當作實數處理。為了得到包含海潮信息的測站位移序列，不對海潮負荷引起的位移進行改正。衛星鐘差、軌道以及硬件延遲選用IGS的精密產品。

由于相鄰天之間解算數據的不連續以及精密產品不一致等原因，可能導致天與天的結果文件出現跳變。為解決此問題，將全年觀測文件以7d為1組進行合并，如，1～7d，6～12d，…，共73個組合觀測文件。對每個組合文件而言，因上述原因導致數據跳變的點只可能出現在文件的首尾時段。因此，對組合觀測文件進行PPP動態解算之后，截取每組2～6d的解算坐標，然后融合成全年的測站坐標系列。

3 實例分析

3.1 數據

選取高緯度地區的3 個CORS 站AC67（57.798°N，152.336°W）、AC43（59.521°N，149.709°W）、AC39（58.610°N，152.394°W）以及中緯度地區的SHAO 測站（31.099°N，121，2°E）2012年全年的觀測數據進行分析。由于缺乏測站坐標真值，故將各測站解算結果與首歷元計算得到的測站坐標求差，并轉換為E、N、U3個方向的坐標序列。實驗表明，此方法會產生因地殼板塊運動導致的坐標漂移。由于海潮負荷引起的位移僅是余弦波的疊加，故可通過一階線性擬合的方式予以消除。考慮到數據量過大對計算機的性能要求也越高，將數據的采樣率由30s稀疏為5min。

3.2 海潮負荷位移參數反演

為驗證動態PPP 解算結果是否包含分潮信息，利用傅里葉變換對結果進行頻譜分析。以測站AC43為例，將測站單天每個方向的坐標序列進行傅里葉變換后所得的振幅在頻域內進行疊乘［14］，可以分離出分潮信號與噪聲。

圖1給出了測站AC43的分潮周期圖。可以很明顯地看出周日潮和半日潮的周期信息，其次是1／3周日潮、1／4周日潮。可以看到，對測站而言，主要受周日潮和半日潮的影響，而1／3潮、1／4潮的影響僅為兩主潮的1／6～1／7，就目前PPP精度而言可以不予考慮。圖中長周期的分潮信息不明顯，這是因為本文選用的數據時間跨度還不夠長，不足以將更長周期的分潮信息提取出來。

圖1 測站AC43海潮頻譜分析周期Fig.1 Ocean tidal periodogram spectrum analysis about station AC43

進一步從量化的角度分析。節點因數和天文相角近似線性變化如圖2、3所示，相較于其他分潮，K2分潮的節點因數和天文相角的變動略為明顯，分別為8%和1%左右。盡管如此，兩者對海潮負荷位移參數的影響也在亞mm級［7］，故在單歷元參數處理中不予考慮，僅對最終結果取節點因數和天文相角的年均值進行改正。將測站位移代入式（2）、（3）式進行參數解算，解算結果與TPXO.7.2海潮模型計算結果進行比較，如表1、2所示。

圖2 全年節點因數變化趨勢Fig.2 Yearly change trend about nodal corrections which account for the modulating effect of the lunar node

圖3 全年天文相角變化趨勢Fig.3 Yearly change trend about astronomical phase angle

3.3 結果評價

選用5個最新的全球海潮模型作為參考（表3），并采用均方根（RMS）［7］來評定PPP反演的參數和海潮模型計算值的差異。計算方法如下：

其中，A、φ分別代表振幅和相位，n（1，2，…，N）為測站號，j為分潮號，k為坐標分量，i為虛數符號。

圖4給出了4個測站反演結果與TPXO.7.2模型計算結果的均方根。造成測站間差異的主要原因與測站所處的地理環境有關。例如，不同地區海潮負荷對地殼形變的影響程度不同，其中SHAO 測站在U方向上RMS 較小，尤其是K2分潮RMS約為其他3個站的1／3，是因為利用海潮模型計算的K2分潮在SHAO 測站的振幅僅為亞mm級，而其他幾個測站K2分潮的振幅達到了幾個mm，說明SHAO 測站受K2分潮的影響（尤其在垂直方向）較其他測站小很多。此外，在淺海地區由于復雜的地勢，海潮模型并不能十分準確地反映分潮信息［15］。4 個測站在K1、K2分潮的均方根最為顯著，主要是受軌道誤差和多路徑效應的影響［14］。

為進一步說明分潮的反演情況以及各模型間的差異，圖5給出了PPP反演參數與5個海潮模型計算值之間的均方根誤差統計結果。可以看到，除U方向Q1分潮各模型計算值有較大差異外，所選的5個全球海潮模型差異不大，說明在研究區域5個模型具有一致的精度。圖中還發現，利用PPP反演出與太陽有關的分潮參數（S2、K2、K1、P1）與模型間的偏差較大，其中在水平方向上K2、K1略為明顯，而在垂直方向上P1分潮的偏差則要高于K1。N方向的均方根誤差整體優于2.5mm，E方向次之優于3.5mm，U方向相對偏差最大但仍優于6 mm，主要是垂直方向的定位精度低于水平精度所致。而影響均方根誤差的原因除了和PPP 反演結果與模型間的系統偏差之外，數據長度也是一個重要原因。

表1 測站半日潮負荷位移參數反演結果與TPXO.7.2模型半日潮計算結果比較Tab.1 Compared PPP－derived results with TPXO.7.2about semi－diurnal constituents

表2 測站周日潮負荷位移參數反演結果與TPXO.7.2模型周日潮計算結果比較Tab.2 Compared PPP－derived results with TPXO.7.2about ter－diurnal constituents

表3 全球海潮模型Tab.3 Global ocean tide model

圖4 PPP反演結果與TPXO.7.2模型之間的比較Fig.4 RMS of the PPP－derived OTL displacement estimates to the TPX.O.7.2model estimates

圖5 PPP反演結果與模型的均方根誤差Fig.5 RMS of the PPP－derived OTL displacement estimates to all ocean tide model estimates

在不考慮系統差的前提下，利用1a的GPS觀測數據反演得到的參數與模型間的均方根誤差與Yuan等［6］計算結果基本相當，這也說明利用本文方法反演海潮負荷位移參數可以較快地得到收斂。

4 結 語

1）本文方法計算的海潮負荷位移參數與模型間的均方根誤差為mm級，與Yuan等［6］方法解算結果基本相當，說明利用本文方法計算海潮負荷位移參數是可行的。

2）除Q1分潮TPXO.7.2 和FES2004 與其他海潮模型有差異外，在中高緯地區本文所選5個海潮模型差異不明顯。

［1］周江存，孫和平.海潮負荷對GPS 基線的影響［J］.大地測量與地球動力學，2005，25（4）：27－32（Zhou Jiangcun，Sun Heping.Effect of Ocean Tide Loading on GPS Baseline Measurement［J］.Journal of Geodesy and Geodynamics，2005，25（4）：27－32）

［2］鄭祎，伍吉倉，王解先，等.GPS精密定位中的海潮位移改正［J］.武 漢 大 學 學 報：信 息 科 學 版，2003，28（4）：405－408（Zheng Yi，Wu Jicang，Wang Jiexian，et al.Ocean Tidal Displacement Corrections in GPS Precision Positioning［J］.Geomatics and Information Science of Wuhan University，2003，28（4）：405－408）

［3］Dragert H，James T S.Ocean Loading Corrections for Continuous GPS：A Case Study at the Canadian Coastal Site Holberg［J］.Geophysical Research Letters，2005，27（14）：2 045－2 048

［4］Zhang X H，Andersen O B.Surface Ice Flow Velocity and Tide Retrieval of the Amery Ice Shelf Using Precise Point Positioning［J］.J Geod，2006，80：171－176

［5］Thomas I D，Peter M A K，Clarke J.A Comparison of GPS，VLBI and Model Estimates of Ocean Tide Loading Displacement［J］.J Geod，2007，81（5）：359－368

［6］Yuan L G，Ding X L，Zhong P.Estimates of Ocean Tide Loading Displacement and Its Impact on Position Time Series in Hong Kong Using a Dense Continuous GPS Network［J］.J Geod，2009，83：999－1 015

［7］郗欽文.精密引潮位展開的精度評定［J］.地球物理學報，1992，35（2）：150－153（Xi Qinwen.The Evaluation on the Precision of the Development of the Tidal Generating Potential［J］.Acta Geophysica Sinica，1992，35（2）：150－153）

［8］周旭華，吳斌，李軍.高精度大地測量中的海潮位移改正［J］.測 繪 學 報，2001，30（4）：327－330（Zhou Xuhua，Wu Bin，Li Jun.Ocean Tide Displacement Corrections in High Precision Geodesy［J］.Acta Geodaetica et Cartographica Sinica，2001，30（4）：327－330）

［9］黃祖珂，黃磊.潮汐原理與計算［M］.青島：中國海洋大學出版社，2005（Huang Zuke，Huang Lei.Tidal Theory and Calculation［M］.Qindao：Publishing House of Ocean University of China，2005）

［10］Héroux P，Kouba J.Precise Point Positioning Using IGS Orbit and Clock Products［J］.GPS Solution，2001，5（2）：12－28

［11］Kouba J.A Guide to Using International GNSS Service（IGS）Products，Natural Resources Canada［EB／OL］.http：／／igscb.jpl.nasa.gov／components／usage.html，2012

［12］Zhang X H，Guo B F，Guo F.Influence of Clock Jump on the Velocity and Acceleration Estimation with a Single GPS Receiver Based on Carrier－Phase－Derived Doppler［J］.GPS Solutions，2013，17（4）：549－559

［13］Ruymbeke M V，Zhu P，CadicheanuNat N.Very Weak Signals（VWS）Detected by Stacking Method According to Different Astronomical Periodicities（HiCum）［J］.Hazards Earth Syst Sci，2007（7）：651－656

［14］King M.Kinematic and Static GPS Techniques for Estimating Tidal Displacements with Application to Antarctica［J］.J Geodyn，2006，41：77－86

［15］李大煒，李建成，金濤勇.利用驗潮站資料評估全球海潮模型的精度［J］.大地測量與地球動力學，2012，32（4）：106－110（Li Dawei，Li Jiancheng，Jin Taoyong.Accuracy Estimation of Recent Global Ocean Tide Models Using Tide Gauge Data［J］.Journal of Geodesy and Geodynamics，2012，32（4）：106－110）