利用PPP反演海潮負(fù)荷位移參數(shù)

鄭 凱 馬 蘭 季 偉 郭博峰

1 武漢大學(xué)測(cè)繪學(xué)院，武漢市珞喻路129號(hào)，430079

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

1 海潮參數(shù)反演方法

測(cè)站的海潮負(fù)荷位移可由11個(gè)主分潮的海潮負(fù)荷位移疊加獲得，其中3 個(gè)為長(zhǎng)周期分潮。由于反演長(zhǎng)周期分潮需要長(zhǎng)時(shí)間的觀測(cè)數(shù)據(jù)，且長(zhǎng)周期潮對(duì)測(cè)站位移的影響較其他分潮要?。?－8］，故在此僅研究8 個(gè)主分潮，即4 個(gè)半日潮（M2、S2、N2、K2）和4個(gè)周日潮（K1、O1、P1、Q1）。

測(cè)站的海潮負(fù)荷位移可表示為：

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

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

式中，A、B為待估參數(shù)，ξk為PPP計(jì)算得到的E、N、U3個(gè)方向的測(cè)站坐標(biāo)序列。將式（2）以法方程形式表達(dá)，用經(jīng)典最小二乘準(zhǔn)則求解待估參數(shù)?？紤]到φ變化范圍從－PI到PI，實(shí)際解算結(jié)果需要加減180°。

2 PPP處理策略

海潮負(fù)荷對(duì)單點(diǎn)定位的影響在沿海地區(qū)可達(dá)mm級(jí)甚至cm級(jí)［1，4］，因此在對(duì)其他誤差精確改正之后，得到的測(cè)站位移即由海潮負(fù)荷引起。

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

由于相鄰天之間解算數(shù)據(jù)的不連續(xù)以及精密產(chǎn)品不一致等原因，可能導(dǎo)致天與天的結(jié)果文件出現(xiàn)跳變。為解決此問(wèn)題，將全年觀測(cè)文件以7d為1組進(jìn)行合并，如，1～7d，6～12d，…，共73個(gè)組合觀測(cè)文件。對(duì)每個(gè)組合文件而言，因上述原因?qū)е聰?shù)據(jù)跳變的點(diǎn)只可能出現(xiàn)在文件的首尾時(shí)段。因此，對(duì)組合觀測(cè)文件進(jìn)行PPP動(dòng)態(tài)解算之后，截取每組2～6d的解算坐標(biāo)，然后融合成全年的測(cè)站坐標(biāo)系列。

3 實(shí)例分析

3.1 數(shù)據(jù)

選取高緯度地區(qū)的3 個(gè)CORS 站AC67（57.798°N，152.336°W）、AC43（59.521°N，149.709°W）、AC39（58.610°N，152.394°W）以及中緯度地區(qū)的SHAO 測(cè)站（31.099°N，121，2°E）2012年全年的觀測(cè)數(shù)據(jù)進(jìn)行分析。由于缺乏測(cè)站坐標(biāo)真值，故將各測(cè)站解算結(jié)果與首歷元計(jì)算得到的測(cè)站坐標(biāo)求差，并轉(zhuǎn)換為E、N、U3個(gè)方向的坐標(biāo)序列。實(shí)驗(yàn)表明，此方法會(huì)產(chǎn)生因地殼板塊運(yùn)動(dòng)導(dǎo)致的坐標(biāo)漂移。由于海潮負(fù)荷引起的位移僅是余弦波的疊加，故可通過(guò)一階線性擬合的方式予以消除?？紤]到數(shù)據(jù)量過(guò)大對(duì)計(jì)算機(jī)的性能要求也越高，將數(shù)據(jù)的采樣率由30s稀疏為5min。

3.2 海潮負(fù)荷位移參數(shù)反演

為驗(yàn)證動(dòng)態(tài)PPP 解算結(jié)果是否包含分潮信息，利用傅里葉變換對(duì)結(jié)果進(jìn)行頻譜分析。以測(cè)站AC43為例，將測(cè)站單天每個(gè)方向的坐標(biāo)序列進(jìn)行傅里葉變換后所得的振幅在頻域內(nèi)進(jìn)行疊乘［14］，可以分離出分潮信號(hào)與噪聲。

圖1給出了測(cè)站AC43的分潮周期圖。可以很明顯地看出周日潮和半日潮的周期信息，其次是1／3周日潮、1／4周日潮。可以看到，對(duì)測(cè)站而言，主要受周日潮和半日潮的影響，而1／3潮、1／4潮的影響僅為兩主潮的1／6～1／7，就目前PPP精度而言可以不予考慮。圖中長(zhǎng)周期的分潮信息不明顯，這是因?yàn)楸疚倪x用的數(shù)據(jù)時(shí)間跨度還不夠長(zhǎng)，不足以將更長(zhǎng)周期的分潮信息提取出來(lái)。

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

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

圖2 全年節(jié)點(diǎn)因數(shù)變化趨勢(shì)Fig.2 Yearly change trend about nodal corrections which account for the modulating effect of the lunar node

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

3.3 結(jié)果評(píng)價(jià)

選用5個(gè)最新的全球海潮模型作為參考（表3），并采用均方根（RMS）［7］來(lái)評(píng)定PPP反演的參數(shù)和海潮模型計(jì)算值的差異。計(jì)算方法如下：

其中，A、φ分別代表振幅和相位，n（1，2，…，N）為測(cè)站號(hào)，j為分潮號(hào)，k為坐標(biāo)分量，i為虛數(shù)符號(hào)。

圖4給出了4個(gè)測(cè)站反演結(jié)果與TPXO.7.2模型計(jì)算結(jié)果的均方根。造成測(cè)站間差異的主要原因與測(cè)站所處的地理環(huán)境有關(guān)。例如，不同地區(qū)海潮負(fù)荷對(duì)地殼形變的影響程度不同，其中SHAO 測(cè)站在U方向上RMS 較小，尤其是K2分潮RMS約為其他3個(gè)站的1／3，是因?yàn)槔煤３蹦Ｐ陀?jì)算的K2分潮在SHAO 測(cè)站的振幅僅為亞mm級(jí)，而其他幾個(gè)測(cè)站K2分潮的振幅達(dá)到了幾個(gè)mm，說(shuō)明SHAO 測(cè)站受K2分潮的影響（尤其在垂直方向）較其他測(cè)站小很多。此外，在淺海地區(qū)由于復(fù)雜的地勢(shì)，海潮模型并不能十分準(zhǔn)確地反映分潮信息［15］。4 個(gè)測(cè)站在K1、K2分潮的均方根最為顯著，主要是受軌道誤差和多路徑效應(yīng)的影響［14］。

為進(jìn)一步說(shuō)明分潮的反演情況以及各模型間的差異，圖5給出了PPP反演參數(shù)與5個(gè)海潮模型計(jì)算值之間的均方根誤差統(tǒng)計(jì)結(jié)果?？梢钥吹?，除U方向Q1分潮各模型計(jì)算值有較大差異外，所選的5個(gè)全球海潮模型差異不大，說(shuō)明在研究區(qū)域5個(gè)模型具有一致的精度。圖中還發(fā)現(xiàn)，利用PPP反演出與太陽(yáng)有關(guān)的分潮參數(shù)（S2、K2、K1、P1）與模型間的偏差較大，其中在水平方向上K2、K1略為明顯，而在垂直方向上P1分潮的偏差則要高于K1。N方向的均方根誤差整體優(yōu)于2.5mm，E方向次之優(yōu)于3.5mm，U方向相對(duì)偏差最大但仍優(yōu)于6 mm，主要是垂直方向的定位精度低于水平精度所致。而影響均方根誤差的原因除了和PPP 反演結(jié)果與模型間的系統(tǒng)偏差之外，數(shù)據(jù)長(zhǎng)度也是一個(gè)重要原因。

表1 測(cè)站半日潮負(fù)荷位移參數(shù)反演結(jié)果與TPXO.7.2模型半日潮計(jì)算結(jié)果比較Tab.1 Compared PPP－derived results with TPXO.7.2about semi－diurnal constituents

表2 測(cè)站周日潮負(fù)荷位移參數(shù)反演結(jié)果與TPXO.7.2模型周日潮計(jì)算結(jié)果比較Tab.2 Compared PPP－derived results with TPXO.7.2about ter－diurnal constituents

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

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

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

在不考慮系統(tǒng)差的前提下，利用1a的GPS觀測(cè)數(shù)據(jù)反演得到的參數(shù)與模型間的均方根誤差與Yuan等［6］計(jì)算結(jié)果基本相當(dāng)，這也說(shuō)明利用本文方法反演海潮負(fù)荷位移參數(shù)可以較快地得到收斂。

4 結(jié) 語(yǔ)

1）本文方法計(jì)算的海潮負(fù)荷位移參數(shù)與模型間的均方根誤差為mm級(jí)，與Yuan等［6］方法解算結(jié)果基本相當(dāng)，說(shuō)明利用本文方法計(jì)算海潮負(fù)荷位移參數(shù)是可行的。

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

［1］周江存，孫和平.海潮負(fù)荷對(duì)GPS 基線的影響［J］.大地測(cè)量與地球動(dòng)力學(xué)，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］鄭祎，伍吉倉(cāng)，王解先，等.GPS精密定位中的海潮位移改正［J］.武 漢 大 學(xué) 學(xué) 報(bào)：信 息 科 學(xué) 版，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］郗欽文.精密引潮位展開的精度評(píng)定［J］.地球物理學(xué)報(bào)，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］周旭華，吳斌，李軍.高精度大地測(cè)量中的海潮位移改正［J］.測(cè) 繪 學(xué) 報(bào)，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］黃祖珂，黃磊.潮汐原理與計(jì)算［M］.青島：中國(guó)海洋大學(xué)出版社，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］李大煒，李建成，金濤勇.利用驗(yàn)潮站資料評(píng)估全球海潮模型的精度［J］.大地測(cè)量與地球動(dòng)力學(xué)，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）