基于載波相位差分的形變監測高精度定位算法

陳凱 孫希延 紀元法 王守華 陳紫強

摘 要：傳統載波相位差分算法在形變監測領域適用性不足，實時動態定位（RTK）精度難以滿足要求，而載波雙差靜態相對定位連續解算時形變跟蹤性能較低等。針對這些問題，在對動靜態算法深入研究的基礎上，提出一種基于載波相位差分的動靜態自適應融合算法。通過方差變化法實時判斷定位結果是否收斂，自適應調節擴展卡爾曼濾波（EKF）狀態先驗估計過程。在收斂時刻增大位置參數的先驗估計誤差的協方差值，使EKF后驗估計過程傾向于信賴測量值;未收斂時刻通過EKF迭代，使EKF后驗估計過程傾向于信賴狀態預測值。實驗結果表明：相比傳統RTK新算法精度有明顯提高，水平達±2mm內，高程達到±4mm內;相比靜態定位則縮減了觀測周期，提高了微小形變跟蹤性能。

關鍵詞：形變監測;實時動態定位;靜態定位;方差變化;定位精度

中圖分類號：TN967.1

文獻標志碼：A

文章編號：1001-9081（2019）04-1234-06

Abstract： Traditional carrier phase difference algorithm is not suitable for deformation monitoring， the accuracy of Real-Time Kinematic （RTK） positioning cannot meet requirements and static relative positioning based on carrier double-differential phase has poor deformation tracking performance with continuous calculations. To solve these problems， based on the deep research of dynamic and static algorithms， a dynamic and static adaptive fusion algorithm based on carrier phase difference was proposed. The convergence of positioning results was judged by variance-change method in real time， then the state priori estimation process of Extended Kalman Filter （EKF） was adaptively adjusted. In the process， the covariance value of priori estimation error of position parameters was increased at the convergence time， so that the posteriori process of EKF tended to trust measured value. EKF iteration was used at the non-convergence time， so that the posteriori process of EKF tended to trust state predicted value. The experimental results show that compared with traditional RTK， the accuracy of the new algorithm is improved， with horizontal accuracy of ±2mm， and altitudinal accuracy of ±4mm. Compared with static positioning， the observation period is reduced， and the tracking performance of micro-deformation is improved.

Key words： deformation monitoring; Real-Time Kinematic （RTK） positioning; static positioning; variance-change; positioning accuracy

0?引言

全球衛星導航系統（Global Navigation Satellite System， GNSS）定位技術在形變監測領域具備自動化程度高、可同時測定三維位移及全天候等優勢，已經廣泛應用在如建筑物、橋梁、大壩和滑坡等監測領域。在實際工程應用中，當形變監測精度要求較高時（例如大壩、船閘監測的精度要求毫米級），或者在沉降監測領域，傳統定位算法很難滿足監測需要[1]。

目前，形變監測常用的GNSS定位算法有實時動態定位（Real-Time Kinematic positioning， RTK）[2]和載波雙差靜態相對定位（簡稱靜態定位）[3]兩種。傳統RTK算法實時性和形變跟蹤性能較高，但是受限于厘米級精度水平[2，4];靜態定位精度可達毫米級，但實際工程應用中需周期性反復采集觀測數據進行處理分析，且連續解算缺乏較好的形變跟蹤能力，

文獻[5]給出了相應的分析，并采用周期為10h的靜態算法監測地面沉降。

針對傳統算法的不足，張小紅等[1]認為微小形變對雙差觀測值的影響遠小于一周，在數據處理時可繞開周跳探測及整周模糊度確定問題;但是該方法依然會受到其他觀測誤差的影響，精度提高有限。文獻[6]提出單雙頻接收機混合監測模式，通過利用雙頻接收機數據建立區域電離層模型，提高單頻機的在沉降領域的監測精度;但是該方法因使用雙頻接收機導致布設成本高，且不適用于小區域監測，應用領域受限。文獻[7]探討了精密單點定位（Precise Point Positioning， PPP）算法在地面沉降、緩變性地質災害監測中適用性問題，雖然采取精細誤差修正模型改正，但精度依然難以滿足毫米級要求。文獻[8]構建的觀測值域內的多路徑誤差模型實時改正觀測值中的多路徑誤差后，三維位置精度可提高約30%。

針對傳統算法在工程應用中的不足，本文提出了基于載波相位差分的動靜態自適應融合算法。該算法采用方差變化法實時檢測解算結果是否達到收斂狀態，在非收斂時刻根據靜態算法的擴展卡爾曼濾波（Extended Kalman Filter， EKF）策略不斷迭代處理，提升其定位精度水平，而在收斂時刻根據實時動態算法的EKF策略自適應調節位置狀態參數及其協方差值，改善其形變跟蹤性能。

1?載波相位差分定位數學模型

監測站接收機r和基準站接收機b基于共視衛星i的載波相位觀測值φ（i）r和φ（i）b分別為：

2?方差變化法

檢測出靜態定位結果的收斂時刻能有效提高新算法的形變跟蹤性能。靜態定位結果時間序列具有由非平穩逐漸收斂到平穩的特點，確定收斂時間，只需要判斷該時間序列中非平穩部分和平穩部分的過渡點。若判斷出時間序列收斂點，根據該收斂點的觀測歷元就可以得到快速靜態定位收斂時間[10]。

非平穩序列的方差值較大，波動也較大;而平穩序列的方差值接近零，波動比較小。因此，本文采用方差變化法，從方差的角度判斷出序列收斂點。

方差變化法的流程如圖1所示。

為了具體說明，需要定義經驗閾值d1、d2。當滿足穩定條件則可認為相鄰方差所對應的定位結果時間序列坐標段已經達到平穩，其收斂點認為是下一個段中觀測歷元最小的點xjk，收斂時間為T=[jk]t（t表示采樣時間）。

經實驗，當監測區域觀測環境較好時，靜態定位能較快達到收斂要求，且發散概率很小。但是惡劣環境下，受到多路徑和周跳等觀測誤差影響，會在某時段造成靜態定位收斂時間較長，甚至出現發散的情況，從而導致方差變化法在經驗閾值下不能得到收斂時間。考慮到現場環境且在滿足GNSS技術監測要求條件下，方差變化法必須滿足以下條件：

1）定義數據段迭代時間限值Tmax，取較好觀測環境下數據段迭代收斂時間的最大值。如果迭代時間超過該值還未收斂，則停止迭代并進入下一數據段處理。

2）閾值d1、d2根據監測環境選取經驗值，但要滿足在Tmax下有85%以上的數據段可收斂。

3?動靜態自適應融合載波相位差分算法

3.1?載波相位差分定位算法

在定位算法中，采用最小二乘法（Least Squares， LS）的處理結果通常會受到誤差與噪聲的影響，顯得既粗糙又雜亂，而EKF能有效降低、分離信號中所含的噪聲量，得到平滑、準確的定位結果[11]。

其中：D為整周模糊度單差轉雙差矩陣;X為狀態向量，包括監測站位置三維坐標、三維速度、三維加速度參數和單差整周模糊度參數，但是在形變精密監測領域中，載體不會發生高動態運動，因此狀態矩陣中通常會舍棄三維速度和三維加速度參數;xr、 yr、zr為監測站三維位置參數;Nirb是基于共視衛星i的站間單差整周模糊度參數，為了避免參考星的變換對雙差整周模糊度固定的影響，在對雙差整周模糊度搜索處理前才將Nirb轉換為雙差形式。

在載波相位差分定位中，當前位置參數利用前一步的預測值確定，如果不發生周跳，認為雙差整周模糊度參數是連續的[9]，因此給出狀態方程及預測值的方差陣：

其中：X^k-1表示k-1歷元的狀態后驗估計向量;A表示狀態轉移矩陣;X-k表示k歷元的狀態先驗估計向量（狀態預測向量）;P^k-1表示k-1歷元的后驗估計誤差協方差矩陣;P-k表示先驗估計誤差的協方差矩陣;Q表示過程噪聲向量的協方差矩陣。

其中：R表示觀測誤差協方差矩陣;yk表示觀測值向量。

實現載波差分定位的核心問題之一是解算整周模糊度的整數解[2]。雙差載波相位方程組（5）中包含監測站位置和雙差整周模糊度未知參量，對于實時動態和靜態定位系統，希望盡可能地快速、簡單而又可靠地完成整周模糊度的求解。

是常用的雙差整周模糊快速求解算法，可在少數歷元下確定雙差整周模糊度的最優值。由式（24）可得到三維坐標固定解：

其中：為三維坐標浮點解;為三維坐標固定解;為整周模糊度固定解向量;為整周模糊度浮點解向量;Q是坐標參數和整周模糊度參數的協方差陣;Q-1是整周模糊度協方差的逆矩陣。

3.2?RTK算法中的EKF先驗估計過程

在RTK定位算法中，如果不發生周跳，認為雙差整周模糊度參數是連續的，在狀態矩陣先驗估計過程中不需要對其進行處理。又因為狀態矩陣中不包含三維速度和三維加速度參數，根據式（18），在EKF的狀態先驗估計過程中，基于當前歷元k的狀態轉移矩陣可設計為：

其中：I為單位矩陣;n為兩站間共視的衛星顆數。為了保證RTK算法的實時動態性，對基于式（18）的EKF狀態矩陣先驗估計過程進行了調整。分為兩種情況：

可以看出，在RTK算法的EKF先驗估計過程中，在每一歷元均基于當前歷元的單點坐標更新狀態矩陣的位置參數部分，并將相應的先驗估計誤差的方差數值設置為900。增大先驗估計誤差的方差，說明當前預測狀態矩陣中三維位置相對不可靠或者測量值相對準確，會使式（20）得到的增益值傾向于信任測量值yk而減少對預測值的依賴[11]。該過程能有效保證算法的實時動態性，在實際監測中，可及時反映出監測體的形變。

3.3?靜態定位中的EKF先驗估計過程

靜態定位中采用了和RTK同樣的狀態轉移矩陣A，為了提高定位的精度，在EKF先驗估計過程中，采用了不同于RTK的處理方式。靜態定位過程中先驗估計狀態矩陣及其協方差矩陣更新為：

靜態定位不同于動態定位，它在非首歷元沒有對狀態矩陣及其協方差矩陣進行更新和初始化，而是通過不斷迭代收斂過程，達到高精度。隨著不斷迭代，其先驗估計誤差的協方差會逐漸變小，說明當前狀態先驗矩陣中三維位置相對可靠或者測量值相對不準確，會使式（20）得到的增益值傾向于信任預測值而減少對測量值的依賴，導致定位結果不斷趨于穩定，精度不斷提高。但是，該算法收斂到穩定狀態時，由于先驗方差值較小，將無法及時反映出監測體形變[14]。

3.4?動靜態自適應融合算法

實際工程應用中，定位精度和形變跟蹤能力是算法性能的兩個重要指標。新算法在基于動靜態算法研究的基礎上，通過方差變化法對收斂時間的判斷，自適應融合動靜態算法EKF的狀態先驗估計過程，使其在滿足精度的同時提高形變跟蹤性能。動靜態自適應融合算法的具體步驟如下：第1步?建立載波相位雙差方程組，如式（5）。針對方程組的非線性，進行線性化得到式（14）。

第2步?對式（14）采用EKF的先驗估計過程進行處理，該過程分為3種情況：1）若當前歷元為首歷元，則先驗估計狀態矩陣及其協方差矩陣根據式（27）～（28）進行更新。

2）若當前歷元為非首歷元，且基于第5步，根據方差變化法判斷當前歷元未滿足收斂要求，則將先驗估計狀態矩陣及其協方差矩陣按照式（31）～（32）進行更新。

3）若當前歷元為非首歷元，且基于第5步，根據方差變化法判斷當前歷元滿足收斂要求，則將先驗估計狀態矩陣及其協方差矩陣按照式（29）～（30）進行更新。

第3步?對式（14）進行EKF的狀態更新過程（狀態矩陣的后驗估計過程），如式（20）～（22）。

第4步?采用LAMBDA/MLAMBDA搜索整周模糊度，通過式（24）得到當前歷元的三維位置固定解。

第5步?保存解算結果，對解算結果形成的時間序列通過方差變化法進行收斂性判斷。若當前歷元判斷收斂，則輸出收斂結果;否則，進入下一歷元解算。

4?算例與分析

為了驗證動靜態自適應融合算法的有效性并比較RTK、靜態定位和本文算法在精密形變監測領域中的實際效果，做了如下3個實驗：1）檢驗方差變化法的有效性;2）基于靜態數據處理實驗，比較動靜自適應融合算法和傳統算法的精度;3）模擬微小形變實驗，檢驗動靜自適應融合算法的形變跟蹤性能。

4.1?方差變化法有效性檢驗

實驗1的是用來檢驗方差變化法的有效性，先構建原始觀測時間序列如下：

式（33）是一個Matlab函數，生成滿足C～N（0，0.005），時間序列長度為1000的隨機數，并采用卡爾曼濾波（Kalman Filter， KF）算法對該序列進行濾波處理，處理結果作為方差變化法輸入的時間序列，最后通過方差變化法對該時間序列進行收斂性分析。實驗中，方差變化法的搜索窗口、收斂閾值和迭代時間限值分別設置為：k=10，d1=5×10-4、d2=1×10-4，Tmax=400s。如圖2所示，實驗分別得到隨機時間序列曲線（原始觀測數據曲線）、KF濾波結果曲線、序列段方差曲線和方差變化曲線。

根據圖2，在歷元達到200s后，方差曲線和方差變化曲線趨于平穩，其中，方差曲線還呈收斂趨勢。通過數據分析，在滿足序列段方差值小于d12，且方差變化值小于d22時，序列收斂時間在210s處。基于KF濾波曲線可以看出，達到收斂時間后，濾波曲線波動明顯變小，已經達到平穩狀態。

基于式（39），隨機生成8組時間序列，分別通過方差變化法處理后得到收斂時間，結果如表1所示。

根據表1，不同隨機時間序列得到的收斂時間之間差異很大，說明不同序列收斂快慢可能不一致。在序列收斂較快時能及時檢測出收斂時間，可有效提高新算法的形變跟蹤能力。

4.2?解算精度分析

實驗2在某大學圖書館樓頂進行，基準站和監測站均采用ublox M8T型兼容美國全球衛星定位系統（Global Positioning System， GPS）和中國北斗衛星系統（BeiDou System， BDS）的雙模單頻接收機采集數據。兩站均保持固定不動，距離約為64m，且已知的精確的坐標采用精密星歷并基于GAMIT軟件處理得到。

數據采集時間為2018年6月20日至21日，觀測86400個歷元，數據采樣間隔為30s，衛星高度截止角為15°，接收到的BDS衛星有9～11顆，GPS衛星有6～10顆，如圖3所示。

實時數據分別采用RTK算法、靜態算法和動靜態自適應融合算法處理。如圖4～6，分別為基于不同算法處理后的坐標與已知精確坐標的差值三維曲線圖。其中，動靜態自適應融合算法中方差變化法的搜索窗口長度、收斂閾值和最長迭代時間分別取經驗值：k=60s，d1=d2=0.001m，Tmax=3600s，實驗過程未出現數據段迭代時間超過Tmax值的情況。

通過將處理得到的定位結果和已知的移動站精確坐標進行比較，可求得各個算法定位結果的外符合精度。表2是基于不同算法下定位結果和已知精確坐標絕對差值的2Sigma統計。

根據圖4～6和表2中的結果可以看出，RTK、靜態定位（24h）和動靜態自適應融合算法的南北向外符合精度MN分別為：±4.2mm、±0.14mm、±1.5mm，東西向外符合精度ME分別為：±4.6mm、±0.26mm、±1.2mm，高程MU外符合精度分別為：±10.2mm、±1.4mm、±3.6mm，其中靜態定位外符合精度與迭代時間成正比。動靜態自適應融合算法水平外符合精度達到±2mm內，高程外符合精度達到±4mm內，較RTK算法有了顯著提高。

4.3?形變跟蹤性能分析

為模擬實際監測環境，實驗3將監測站天線放置到觀測環境較差的教學樓頂，基準站位置不變，兩站之間距離約521m。其中，監測天線北側受到移動公司基站機房遮擋，南側受到樹木影響。實驗中將監測站天線放置到模擬形變的可移動平臺上，該平臺可在相互垂直的三個方向移動，測微器可精確測定移動量，精確至0.01mm，故可視為已知值，未移動前的精確坐標可基于實驗2通過GAMIT軟件得到。

實驗在2018年6月25日至26日進行。開始時間為25日10時，從實驗開始至26日10時保持監測天線固定不動，以后每隔1h在南北向移動0.5mm，在高程方向移動1mm，共移動10h。觀測172800個歷元，數據采樣間隔為30s，衛星高度截止角為25°，以減小低仰角衛星的影響，接收到的BDS衛星有7～10顆，GPS衛星有5～8顆。

實時數據分別由RTK、靜態定位、動靜態自適應融合算法處理，其中動靜態自適應融合算法中方差變化法的搜索窗口長度k、收斂限值和迭代時間限值分別設置為：k=60s，d1=d2=0.0025m，Tmax=3600s，相比實驗2，由于為保證惡劣環境下收斂，增大收斂閾值，處理過程中未出現收斂時間超出Tmax值的情況。

如圖7～8所示，為不同算法處理結果偏移量和實際偏移量在南北和高程方向的比較。根據表3和圖7～8，RTK、靜態定位（24～48h）和動靜態自適應融合算法的南北向外符合精度MN分別為：±5.1mm、-3.3mm、±1.8mm，高程外符合精度MU分別為：±18.3mm、-7.8mm、±3.9mm。其中，RTK算法外符合精度大于微小形變量，處理結果無法直觀反映出微小形變;監測天線移動后，靜態定位算法外符合精度明顯下降，說明其連續解算時形變跟蹤性能較差;動靜態自適應融合算法不管在水平方向還是高程方向形變跟蹤性能均優于其他兩種算法。

5?結語

本文介紹了RTK和靜態算法基于EKF濾波過程的不同對精度和形變跟蹤性能的影響，提出了一種適用于精密形變監測領域的動靜態自適應融合新算法。從實驗結果可看出，與RTK算法相比，新算法精度顯著提高，符合精密形變領域的監測精度要求;相比靜態定位和RTK算法，其微小形變跟蹤性能較好。這表明新算法融合了高精度和微小形變跟蹤性能的優勢，有利于GNSS技術在精密形變監測領域的應用。

參考文獻（References）

[1] 張小紅， 李征航， 徐紹銓.高精度GPS形變監測的新方法及模型研究[J]. 武漢大學學報（信息科學版）， 2001， 26（5）： 451-454. （ZHANG X H， LI Z H， XU S Q. A new method for high accuracy deformation monitor with GPS[J]. Geomatics and Information Science of Wuhan University， 2001， 26（5）： 451-454.）

[2] 王世進， 秘金鐘， 李得海， 等.GPS/BDS的RTK定位算法研究[J]. 武漢大學學報（信息科學版）， 2014， 39（5）： 621-625. （WANG S J， BEI J Z， LI D H， et al. Real-time kinematic positioning algorithm of GPS/BDS[J]. Geomatics and Information Science of Wuhan University， 2014， 39（5）： 621-625.）

[3] 袁本銀， 鮑志雄， 潘國富.GPS/BD靜態相對定位解算研究與實現[J]. 測繪， 2012， 35（6）： 247-250. （YUAN B Y， BAO Z X， PAN G F. Research and realization of GPS/BD static relative positioning solution[J]. Surveying Mapping， 2012， 35（6）： 247-250.）

[4] HE H B， LI J L， YANG Y X， et al. Performance assessment of single-?and dual-frequency BeiDou/GPS single-epoch kinematic positioning[J]. GPS Solutions， 2014， 18（3）： 393-403.

[5] 王利， 張勤， 范麗紅， 等. 北斗/GPS融合靜態相對定位用于高精度地面沉降監測的試驗與結果分析[J]. 工程地質學報， 2015， 23（1）： 119-125. （WANG L， ZHANG Q， FAN L H， et al. Experiment and results of high precision land subsidence monitoring using fused BDS/GPS data and static relative positioning[J]. Journal of Engineering Geology， 2015， 23（1）： 119-125.）

[6] 張超， 戴吾蛟， 石強， 等. 電離層延遲對單頻GPS點的影響及改正方法研究[J]. 武漢大學學報（信息科學版）， 2018， 43（3）： 471-477. （ZHANG C， DAI W J， SHI Q， et， al. Influence of ionosphere delay on single frequency GPS point and its correction method[J]. Geomatics and Information Science of Wuhan University， 2018， 43（3）： 471-477.）

[7] 王利. 地質災害高精度GPS監測關鍵技術研究[J]. 測繪學報， 2015， 44（7）： 826-826. （WANG L. A study on key technology of high precision GPS monitoring for geological hazard[J]. Acta Geodaetica et Cartographica Sinica， 2015， 44（7）： 826-826.）

[8] 易清根， 劉心龍， 劉萬科. 基于EMD和小波的GPS/BDS變形監測中的多路徑誤差削弱方法[J]. 大地測量與地球動力學， 2017， 37（5）： 462-466. （YI Q G， LIU X L， LIU W K. Multipath mitigation method by combined EMD and wavelet in GPS/BDS deformation monitoring[J]. Journal of Geodesy and Geodynamics， 2017， 37（5）： 462-466.）

[9] 王趁香， 葛茂榮， 祝會忠， 等.BDS/GPS組合靜態相對定位算法探討[J]. 導航定位學報， 2017， 5（2）： 93-97， 102. （WANG C X， GE M R， ZHU H Z， et al. Discussion on static relation positioning algorithm of combined BDS/GPS[J]. Journal of Navigation and Positioning， 2017， 5（2）： 93-97， 102.）

[10] 葛春蕾.時間序列的變點估計的相合性和收斂速度[D]. 合肥： 合肥工業大學， 2007： 25-29. （GE C L. The consistency and convergence rate of the change-point estimation of rime series[D]. Hefei： Heifei Institute of Technology， 2007： 25-29.）

[11] 謝鋼.GPS原理與接收機設計[M]. 北京： 電子工業出版社， 2009： 126-131. （XIE G. Principles of GPS and Receiver Design[M]. Beijing： Publishing House of Electronic Industry， 2009： 126-131.）

[12] TEUNISSEN P J G. The least-squares ambiguity decorrelation adjustment： a method for fast GPS integer ambiguity estimation[J]. Journal of Geodesy，1995，70（2）：65-82.

[13] CHANG X W， YANG X， ZHOU T. MLAMBDA： a modified LAMBDA method for integer least-squares estimation[J]. Journal of Geodesy， 2005， 79（9）： 552-565.

[14] 陳鵬.基于單歷元解算的實時變形監測技術的研究與算法實現[D].沈陽：東北大學，2014：1-24.（CHEN P. Research and algorithm development of real-time deformation monitoring technology based on GPS single epoch calculation. Shenyang： Northeastern University， 2014： 1-24.