第一次看到“PPP”這三個字母十個人里有八個會先想到那個撥號上網時代的Point-to-Point Protocol。但在衛星導航這個圈子里PPP的另一個含義——Precise Point Positioning精密單點定位——才是真正值得耐心研究的硬骨頭。手里一臺雙頻GNSS接收機不架基準站不拉電臺和4G差分鏈路單靠接收機自身把定位精度做到厘米級背后就是一套完整的PPP處理流程和一套不能出錯的誤差修正公式。這篇內容適合正在做GNSS數據處理、RTK轉PPP、或者想搞懂“接收機輸出的坐標到底怎么算出來”的朋友。我會把PPP從觀測方程到誤差修正公式逐項拆開再把從RINEX文件到最終坐標的完整流程走一遍最后夾帶一些我實際處理數據時踩過的坑。1. 這個縮寫經常被認錯PPP到底是什么1.1 和撥號協議無關精密單點定位的前世今生精密單點定位的概念最早是1997年Zumberge等人在JPL提出來的核心思路非常簡單直白既然差分定位能消掉衛星鐘差、軌道誤差、電離層延遲這些共性誤差那我能不能不用基準站直接把誤差一個一個“算明白”最終靠單臺接收機實現精密絕對定位這個想法在當時非常超前因為要做到這一點需要一個前提——必須有高精度的衛星軌道和衛星鐘差產品。在IGS國際GNSS服務還沒成氣候的上世紀九十年代這幾乎是奢望。直到IGS的精密星歷和精密鐘差產品逐漸穩定發布PPP才從論文變成可用的工具。PPP的本質是絕對定位它輸出的坐標直接定義在ITRF這樣的全球參考框架下。和RTK依賴基準站的相對定位不同PPP沒有作用距離的概念只要你頭頂上有衛星、手上能拿到精密產品在全球任何地方都能得到一致的定位結果。1.2 在RTK壟斷的高精度賽道里PPP憑什么分一杯羹很多人會問RTK精度又高、收斂又快為什么還要用PPP這個問題的答案取決于你身處什么場景。RTK的精度確實好水平能到1-2厘米但它有個天生弱點必須依靠基準站。短基線RTK作用距離一般不超過20公里網絡RTK雖然覆蓋廣也需要運營商布設連續運行參考站網還要解決數據鏈路的問題。你在海上、沙漠、山區、邊境線上作業時很可能根本沒有CORS站可用這時候PPP就是唯一的答案。我給一個直觀的對比對比項RTKPPP基準站需求必須或依賴CORS網不需要定位類型相對定位絕對定位參考框架隨基準站而定ITRF等全球框架典型精度厘米級厘米級到分米級收斂時間秒級到分鐘級十到幾十分鐘作業半徑受基線長度限制全球覆蓋主要誤差處理差分消除逐項建模修正PPP的代價也很明確收斂慢。這是因為單臺接收機無法通過差分把模糊度快速固定收斂過程受觀測幾何、衛星數量和誤差模型質量影響很大。但如果你做的是長時段的靜態監測比如滑坡監測、海平面監測、形變測量PPP完全可以勝任而且省掉了一堆基準站維護的麻煩。2. 建模思路決定成敗觀測方程與兩種主流實現2.1 無電離層組合用加權差消除一階電離層PPP中最經典的觀測模型是無電離層組合Ionosphere-FreeIF它把雙頻觀測值按頻率平方的比例加權組合從而消掉一階電離層延遲。這個模型從上世紀九十年代沿用至今現在很多業務化PPP軟件仍然以它為主打。先說偽距和相位的原始觀測方程。對頻率$f_i$的偽距$P_i$忽略噪聲的情況下可以寫成$$ P_i \rho c(dt_r - dt^s) T \frac{I}{f_i^2} B_i $$載波相位觀測值$L_i$則多了一個模糊度項并且電離層延遲符號相反$$ L_i \rho c(dt_r - dt^s) T - \frac{I}{f_i^2} \lambda_i N_i b_i $$其中$\rho$是衛星到接收機的幾何距離$dt_r$是接收機鐘差$dt^s$是衛星鐘差$T$是對流層延遲$I/f_i^2$是電離層延遲$\lambda_i$和$N_i$分別是波長和整周模糊度。如果直接用$P_1$和$P_2$做線性組合取$$ P_{IF} \frac{f_1^2 P_1 - f_2^2 P_2}{f_1^2 - f_2^2} $$你會發現電離層延遲項被消掉了類似地載波相位也有$$ L_{IF} \frac{f_1^2 L_1 - f_2^2 L_2}{f_1^2 - f_2^2} $$組合后的PPP觀測方程就變成了$$ P_{IF} \rho c(dt_r - dt^s) T \varepsilon_P $$$$ L_{IF} \rho c(dt_r - dt^s) T \lambda_{IF} N_{IF} \varepsilon_L $$這就是經典PPP的核心。無電離層組合最大的優點是模型簡單不需要電離層參數直接用雙頻觀測值就能把一階電離層項“數學消除”。缺點也很明顯組合后的觀測噪聲大約是原始觀測值的兩到三倍而且模糊度$N_{IF}$不再具有整數特性導致無法直接固定整周模糊度。2.2 非差非組合模型直接估計電離層延遲最近十年越來越多的PPP軟件轉向非差非組合UncombinedUC模型。這個模型不強行消電離層而是把傾斜電離層延遲當作參數直接估計。它保留了原始偽距和相位觀測值模型寫起來更貼近物理本質$$ P_i \rho c(dt_r - dt^s) T \gamma_i I_r \varepsilon_{P_i} $$$$ L_i \rho c(dt_r - dt^s) T - \gamma_i I_r \lambda_i N_i \varepsilon_{L_i} $$這里$\gamma_i$是電離層系數$\gamma_i f_1^2 / f_i^2$$I_r$是接收機天頂方向的一階電離層延遲通常還要用投影函數把它投影到衛星視線方向。我個人的體會是非差非組合模型有四個實打實的好處觀測噪聲小因為不需要做組合放大三頻四頻衛導系統擴展起來非常自然直接多寫幾個方程就行除了定位坐標你還能順帶輸出電離層延遲和對流層延遲這對氣象和水汽反演很有用模糊度保持原始波長為后續固定整數模糊度提供了更好的基礎。代價是狀態向量里多了一堆電離層參數卡爾曼濾波的計算量變大而且電離層參數需要合理的隨機游走模型否則估計出來會震蕩得很厲害。2.3 參數估計里的隱藏變量無論你用IF模型還是UC模型最終都要落到參數估計上。PPP里最關鍵的部分不是那些顯眼的坐標參數而是一堆“隱藏變量”接收機鐘差、對流層天頂濕延遲、模糊度參數。以經典IF模型為例單系統PPP的狀態向量一般可以寫成$$ X \begin{bmatrix} x y z cdt_r ZWD N_{IF}^1 N_{IF}^2 \cdots N_{IF}^n \end{bmatrix}^T $$其中$(x,y,z)$是接收機坐標靜態模式里這三個參數過程噪聲設為零動態模式則要估計速度和加速度$dt_r$是接收機鐘差通常當作白噪聲或隨機游走過程$ZWD$是對流層天頂濕延遲一般用隨機游走建模后面的$N_{IF}$是各顆衛星的無電離層組合模糊度視為常數發生周跳時重置。估計方法主流是卡爾曼濾波。卡爾曼濾波的好處是實時性好適合逐歷元處理而且狀態預測和觀測更新天然支持PPP中“坐標常量鐘差白噪聲模糊度常量”這種混合模型。如果做后處理也可以用最小二乘批處理但批處理遇到幾百個小時的數據時矩陣規模很嚇人卡爾曼濾波還是更順手。3. 誤差修正公式逐項拆解從衛星端到地面PPP“精密”二字全靠誤差修正撐起來。少了任何一項坐標偏差就可能從厘米級變成分米級甚至米級。我把誤差來源分成三類逐個公式拆開講。3.1 衛星端軌道、鐘差、天線相位中心與相對論精密軌道和鐘差PPP的衛星位置不能靠廣播星歷必須用IGS的精密軌道產品。SP3格式的軌道文件通常每5分鐘或15分鐘一個節點你要用拉格朗日多項式或切比雪夫多項式插值到觀測時刻。我自己的經驗是用九階或十階拉格朗日插值用SP3節點以外的一點時間做收尾比如前后多外推一兩個歷元防止邊緣震蕩。精密鐘差文件CLK格式給出衛星鐘差通常30秒間隔。注意SP3和CLK必須來自同一產品系列而且時間基準要嚴格對齊否則會出現幾厘米到幾十厘米的系統差。衛星天線相位中心衛星天線相位中心改正要從igs14.atx天線文件中讀取。這個文件里有衛星天線的PCO相位中心偏移和PCV相位中心變化。衛星PCO是星固坐標系下的一個固定矢量但在衛星姿態旋轉時它對觀測距離的投影會變化。近似的改正公式為$$ \Delta \rho_{ant,sat} \mathbf{e}{rec}^{sat} \cdot \mathbf{R}{body} \cdot \mathbf{r}_{PCO,sat} $$其中$\mathbf{e}{rec}^{sat}$是接收機到衛星的單位矢量$\mathbf{R}{body}$是星固坐標系相對慣性系的旋轉矩陣$\mathbf{r}_{PCO,sat}$是衛星天線PCO在星固系下的坐標。PCV則通常按衛星天底角查表。相對論效應衛星鐘在軌的高速運動和地球重力場會導致星鐘頻率偏移。IGS精密鐘差產品已經吸收了大部分相對論效應但還遺留了一個周期性相對論項需要自行修正$$ \Delta \rho_{rel} - \frac{2}{c}(\mathbf{r}{sat} \cdot \mathbf{v}{sat}) $$這里的$\mathbf{r}{sat}$和$\mathbf{v}{sat}$是信號發射時刻衛星在地固系中的位置和速度矢量。這項改正雖然量級不大最大幾厘米但如果不修正會在坐標時間序列里引入明顯的長周期波動尤其在做靜態監測時特別刺眼。相位纏繞衛星天線在運行中會不斷旋轉以保持對日定向這會導致載波相位觀測值產生一個與幾何無關的額外變化叫相位纏繞Phase Wind-up。動態定位時必須修正靜態長時間觀測如果忽略它模糊度估計會受污染。相位纏繞的改正量可以寫成$$ \Delta \phi \delta \phi 2\pi N $$其中$\delta \phi$由天線體坐標系的等效偶極子方向和接收機天線方向決定$N$由連續跟蹤的相位纏繞整數計數決定。大部分軟件會自動處理但如果你自己寫代碼不要漏掉這一項。3.2 傳播路徑電離層高階項與對流層濕延遲電離層IF組合已經把一階電離層項消掉了但二階和三階電離層項仍然殘留。對絕大多數動態定位場景來說殘差只有毫米到厘米級可以不管。如果你做的是最高精度的靜態PPP比如地震形變研究就要引入高階電離層改正公式基于地磁場模型$$ \Delta I_2 \frac{7527c}{f^3} \int N_e B \cos\theta , ds $$二階項的修正量大約在0-2厘米范圍。大多數業務用戶忽略它問題不大但你要知道它存在。對流層對流層延遲分干濕兩項。干延遲比較穩定可以用Saastamoinen模型直接算$$ ZHD \frac{0.002277 \cdot P}{1 - 0.0026\cos(2\varphi) - 0.00028h} $$其中$P$是測站氣壓hPa$\varphi$是緯度$h$是海拔高度km。干延遲算出來以后不需要估計當作已知值改正掉就行。濕延遲$ZWD$比較難建模通常作為未知參數用卡爾曼濾波去估計。把天頂總延遲$ZTD$拆成干濕兩部分$$ ZTD ZHD ZWD $$無論是干延遲還是濕延遲從天頂方向投影到衛星視線方向都需要映射函數。現在常用的是GMFGlobal Mapping Function$$ mf(\varepsilon) \frac{1 \frac{a}{1\frac{b}{1c}}}{\sin\varepsilon \frac{a}{\sin\varepsilon \frac{b}{\sin\varepsilon c}}} $$其中$\varepsilon$是衛星仰角$a$、$b$、$c$是和測站位置、年積日相關的經驗系數。對流層延遲沿視線方向的總改正就是$$ T(\varepsilon) ZHD \cdot mf_{dry}(\varepsilon) ZWD \cdot mf_{wet}(\varepsilon) $$我剛開始做PPP時犯過一個錯干濕延遲用了同一個映射函數結果在低仰角衛星上總是出現異常殘差。后來換成分別用干濕映射函數情況立刻改善。3.3 地球物理效應固體潮、海潮與極潮很多人在做PPP時容易忽略地球物理效應因為它們不是“導航誤差”而是“地球形變”。但PPP是絕對定位坐標框架本身就在隨著固體潮、海潮變化所以必須修正。固體潮固體潮是日月引力導致的地殼彈性形變垂直方向最大可達30厘米水平方向也有幾厘米是所有地球物理效應里量級最大的絕對不能忽略。IERS2010公約給出了計算公式寫成簡化形式是$$ \Delta \mathbf{r} \sum_{j1}^{2} \frac{GM_j}{GM} \frac{r^4}{R_j^3} \left[ (3l_2(\hat{R}_j \cdot \hat{r}))\hat{R}_j \left(\frac{3h_2}{2} - \frac{3l_2}{2}\right)(\hat{R}j \cdot \hat{r})^2 \hat{r} \right] \Delta \mathbf{r}{perm} $$其中$GM_j$和$GM$分別是日月和地球的引力常數$r$是地心到測站的距離$R_j$是日月到地心的距離$\hat{R}j$和$\hat{r}$是對應的單位矢量$h_2 \approx 0.6078$和$l_2 \approx 0.0847$是Love數和Shida數。后面還要補一項永久潮汐改正$\Delta\mathbf{r}{perm}$。實際程序里這個過程并不復雜逐歷元算一下加到測站坐標上就行。但如果你不修正固體潮靜態PPP的垂直分量時間序列會呈現明顯的半日波和全日波解算精度會大打折扣。海洋負荷潮汐海洋負荷潮汐是潮汐水體對沿岸地殼的加載形變在離海近的測站上可達幾厘米。它需要用海潮模型比如FES2014b、TPXO9輸出負荷潮汐系數再結合潮汐諧波參數算改正。如果你的測站離海岸線幾十公里以內強烈建議加上海潮改正。否則你會看到坐標時間序列里有一個明顯的12小時左右周期分量而這個分量在對流層估計里是消不掉的。極潮極潮是由地球自轉軸相對地殼的微小擺動引起的形變最大量級約2.5厘米發生在45度緯度附近。IERS2010給出了簡便公式$$ \Delta r_{polar} -33\sin 2\varphi \cdot (m_1 \cos\lambda m_2 \sin\lambda) $$其中$m_1$、$m_2$是極移參數$\lambda$是經度。極潮改正量雖小但長期監測數據處理里加上后坐標殘差會明顯更平。3.4 與接收機相關的細節接收機天線相位中心同樣有PCO和PCV可以通過igs14.atx文件查表改正。如果你用的接收機天線在atx文件里查不到至少要把PCO大致值輸進去否則會造成厘米級系統偏差。還有一項容易忽略的是地球自轉改正Sagnac效應。信號從衛星傳到地面大約要0.07秒這個時間地球轉過了大約30米所以計算幾何距離時必須在衛星坐標上做地球自轉改正或者在方程里顯式加入Sagnac項$$ \Delta \rho_{sagnac} \frac{\omega_e}{c}(x^s y_r - y^s x_r) $$其中$\omega_e$是地球自轉角速度$(x^s,y^s,z^s)$是衛星坐標$(x_r,y_r,z_r)$是接收機近似坐標。這一項不改正定位結果會直接偏出去幾十米屬于那種“一步錯步步錯”的基礎操作。4. 從觀測值到坐標的完整處理鏈講了這么多公式現在把它們串成一條完整的處理鏈。一套標準的PPP處理流程大致分成四步。4.1 第一步數據采集與RINEX規范化接收機設置上PPP要求至少雙頻觀測值最好能用三頻或四頻。采樣率方面靜態測量建議1-30秒動態測量按需求來但要注意數據量和后處理時間之間的平衡。天線架設時天線高的量測要精確到毫米級因為PPP是絕對定位天線高誤差會直接進到坐標結果里。周圍不能有遮擋物截止仰角建議先放到5度后期處理再根據數據質量調整。采集完原始數據后你需要把接收機廠商的私有格式轉成標準RINEX格式。現在各家廠商的轉換軟件都做得比較友好但要注意RINEX版本。我強烈建議用RINEX 3.04以上版本因為它對多系統、多頻率的支持要好得多字段也更規范。4.2 第二步精密產品下載與對齊這一步決定了你能不能解算出高精度坐標。IGS精密產品按時間延遲分幾檔產品類型軌道精度鐘差精度時延IGS最終產品2.5 cm75 ps約14天IGR快速產品2.5 cm75 ps約1天IGU超快速產品3 cm預報5 cm約200 ps實時/預報實時產品RTS/SSR3 cm約100 ps實時廣播做高精度事后靜態處理首選IGS最終產品。做時效要求高的工程可以用IGR快速產品。實時PPP則要接收SSR改正數通過NTRIP協議從IGS RTS等實時服務商獲取。下載完產品之后做“對齊”檢查這一步經常被新手忽略檢查觀測數據的時標是GPST還是UTC和精密產品的時標是否一致檢查SP3的參考框架和atx天線文件版本是否匹配比如用igs14的SP3就要配igs14.atx別混用igs08的產品和igs14的天線文件檢查CLK文件的衛星鐘差是否包含所有你需要解算的衛星系統。4.3 第三步逐歷元解算和濾波調參完整流程可以描述成下面的循環讀入一個歷元的雙頻偽距和載波相位觀測值用歷元內所有衛星的觀測值做一次標準單點定位SPP給出接收機位置的初始值根據SP3插值出每顆衛星在信號發射時刻的位置根據CLK插值出衛星鐘差計算前面章節列出的所有誤差修正項相對論、天線相位中心、相位纏繞、固體潮、海潮、極潮、Sagnac效應等形成觀測方程代入濾波器的預測狀態計算新息執行新息檢驗如果某顆衛星的殘差超過預設閾值給它降權或剔除卡爾曼濾波更新得到當前歷元的位置、鐘差、ZWD和模糊度進入下一歷元重復上述過程。卡爾曼濾波的調參是最考經驗的部分。濾波器的過程噪聲設置如果太緊參數會被“鎖死”跟不上真實變化設置太松解算結果又會被觀測噪聲帶得亂跳。我平時常用的初始參數經驗值放在第6章的統一配置清單里。4.4 第四步質量評估與坐標輸出處理完成后軟件會輸出一個坐標時間序列。靜態PPP的最終坐標一般取收斂后一段時間比如最后30-60分鐘的平均值用這段時間的標準差來評價內符合精度。如果有已知坐標可以計算外符合精度。評價一個PPP結果質量時我會同時看三樣東西位置時間序列的收斂曲線是否平滑有沒有反復震蕩收斂后的標準差水平分量10毫米以內算不錯垂直分量20毫米以內算合格模糊度殘差和新息序列有沒有系統性偏移如果有說明某個誤差修正項沒過關。5. 收斂時間這個老大難決定因素與壓榨技巧5.1 收斂時間為什么跑不短PPP最讓用戶抓狂的問題就是收斂慢。傳統單系統GPS靜態PPP收斂到水平10厘米以內順利的話要30到60分鐘不順利的話兩三個小時也正常。收斂慢的根本原因是模糊度與接收機鐘差之間的強相關性。單臺接收機時接收機鐘差和模糊度參數幾乎線性相關卡爾曼濾波需要足夠多的歷元和衛星幾何變化才能把它們逐步分離。此外無電離層組合的噪聲被放大也拖慢了收斂速度。之前網上搞過不少“快速收斂”的噱頭實際上PPP收斂就是個信息量積累的過程。幾何分布差、衛星數少、觀測噪聲大的時候誰來了都沒辦法。5.2 我實測有效的幾招經過大量數據實測我把有效縮短收斂時間的策略按見效程度排個序多系統融合這是效果最明顯的一招。GPSBDSGalileoGLONASS四系統融合后衛星數輕松超過25顆幾何構型大幅改善收斂時間能比單GPS縮短一半以上。我做一個靜態測試時單GPS收斂用了大約42分鐘四系統融合大約14分鐘就進入了水平10厘米以內。多頻觀測BDS-3、Galileo都有三頻甚至四頻信號非差非組合模型可以直接利用多頻觀測值進一步降低噪聲。三頻觀測值在高仰角下對模糊度解算的幫助非常直觀。合理約束ZWD和電離層對流層濕延遲的隨機游走譜密度設置要匹配實際氣象變化。譜密度設太大ZWD參數會把位置誤差吸收掉設太小又跟不上真實的水汽變化。一般靜電力可以設在$3\times10^{-8}$到$3\times10^{-7} , \mathrm{m^2/s}$具體值要看數據處理時長。先驗坐標約束如果做的是長期監測測站坐標可以先做一個小時左右的快速解把結果作為先驗約束反饋進去能明顯加速重新收斂。合理設置收斂判定條件不要用“解算到第幾分鐘就穩定”這種模糊說法。科學的辦法是設定一個滑動窗口比如連續10個歷元內N、E、U三個方向的坐標變化量都小于某個閾值常見取2-5厘米再判定為收斂。這樣判定結果可重復也方便不同批次數據之間做對比。6. 實戰中繞不開的坑踩坑記錄與參數清單6.1 精密產品時延導致的時間錯位我曾經處理一批IGU超快速產品的數據結果N方向整體偏了大約10厘米。排查了好久最后發現是因為IGU產品的鐘差時間序列比SP3軌道時間序列晚了整整一個采樣間隔插值時沒有對齊產生了系統性偏差。這個坑的排查思路很簡單在解算前隨機選幾顆衛星把插值得到的衛星軌道、衛星鐘差和廣播星歷的對應值畫在同一個時間軸上看看有沒有時間偏移。如果發現錯位修正產品的時間標簽后再處理。6.2 周跳漏探引發的災難性模糊度重收斂PPP中周跳一旦漏探等于把一顆衛星的模糊度估計值錯誤地帶到后續所有歷元。最壞的情況是收斂完成后一個漏探周跳又把坐標拉偏到分米級需要重新收斂白干半小時。周跳探測我強烈建議用雙組合策略。GFGeometry-Free組合對周跳敏感$$ L_{GF} L_1 - L_2 $$它對電離層變化也敏感可能把電離層擾動誤判成周跳。MWMelbourne-Wübbena組合則不受幾何、鐘差和電離層影響$$ L_{MW} \frac{1}{f_1 - f_2}(f_1 L_1 - f_2 L_2) - \frac{1}{f_1 f_2}(f_1 P_1 f_2 P_2) $$MW組合的缺點是偽距噪聲大小周跳容易漏。把GF和MW聯合起來用再配合TurboEdit算法基本能覆蓋絕大部分周跳場景。而且要注意GF組合中的周跳$\Delta N_1$和$\Delta N_2$相等的時候GF組合檢測不出來但MW組合幾乎是這個特殊場景的唯一救星。6.3 一套可以照抄的PPP配置清單最后給出一套我常用的靜態PPP配置參數供你當作初始值參考。不同場景可以在這個基礎上微調參數項配置值觀測值類型雙頻或無電離層組合或多頻非差非組合截止仰角7度采樣間隔30秒后處理靜態衛星天線PCO/PCVigs14.atx接收機天線PCO/PCVigs14.atx對流層干延遲Saastamoinen模型 GMF干映射函數對流層濕延遲GMF濕映射函數 隨機游走估計固體潮IERS2010海潮FES2014b近海必用極潮IERS2010相位纏繞修正相對論修正周期項卡爾曼濾波過程噪聲ZWD$1\times10^{-7} , \mathrm{m^2/s}$接收機鐘差過程噪聲白噪聲或$10^4 , \mathrm{m^2/s}$模糊度初始方差$10^6 , \mathrm{cycle^2}$這套配置不是拍腦袋寫的都是實際處理數據過程中一步步調出來的經驗值。你在自己數據上跑出來的最優過程噪聲可能和我略有差異但至少能保證一個不錯的起點。做PPP數據處理三年多下來我最深的體會是這套技術不神秘但也不簡單。每一項誤差修正單獨看都不難難的是把它們全部串起來并在實際數據中排查出那一個讓你結果偏幾厘米的隱蔽錯誤。尤其是固體潮、海潮這類地球物理改正看起來和“定位”沒什么關系但在靜態高精度場景里它們恰恰是決定垂直分量能不能進2厘米的關鍵。所以如果你正在調PPP解算結果先別急著懷疑算法和軟件把你的誤差修正項逐項過一遍多數問題都在那里。