1、1 基礎知識1.1 GPS精密單點定位的基本原理GPS精密單點定位一般采用單臺雙頻GPS接收機，利用IGS提供的精密星歷和衛星鐘差，基于載波相位觀測值進行的高精度定位。觀測值中的電離層延遲誤差通過雙頻信號組合消除，對流層延遲誤差通過引入未知參數進行估計。1.2 時間系統RTKLIB內部使用GPST(GPST時間)用于GNSS的數據處理和定位算法。數據在RTKLIB內部處理之前，需要轉換成GPST時間。使用GPST的原因是避免處理潤秒。RTKLIB使用以下結構體表示時間：typedef structtime_t time; /* time(s) expressed by standard tim

2、e_t */double sec; /* fraction of second under 1 s */ gtime_t;1.2.1 GPST和UTC(Universal Time Coordinated)關系參考【圖1】，參考【圖2】：圖1 轉換關系公式1 / 40圖 2通過使用GPS導航信息中的UTC參數，GPST到UTC或者UTC到GPST之前的轉換可以用更準確的表達方式，如【圖3】。圖 3這些參數是由GPS導航消息提供的。1.2.2 BDT（北斗導航衛星系統時間）BDT（北斗導航衛星系統時間）是一個連續的時間系統，沒有潤秒。開始歷元的時間是【UTC 2006年1月1號 00:00:00

3、】。北斗時間計算公式【圖4】：圖 4UTC和GPST時間轉換同上面的GPS一樣，只不過UTC參數來自與北斗導航信息中。1.3 坐標系統接收機和衛星的位置在RTKLIB中表示為在ECEF(地心地固坐標系)坐標系統中的X, Y, Z組件。1.3.1 大地坐標到ECEF坐標的轉換轉換公式如【圖5】。第三個公式最后一行有錯有錯，應該為: (v(1 e2)+h)sin圖 5參數說明： a : 地球參考橢球的長半徑 f : 地球參考橢球的扁平率h: 橢球高度 : 緯度 : 經度當前版本的RTKLIB使用的值為【圖6】：圖 6圖7 參考橢球體1.3.2 ECEF坐系到大地坐標的轉換轉換公式如【圖8】圖 81

4、.3.3 本地坐標到ECEF坐標的轉換在接收機位置的本地坐標，也被稱為ENU坐標，通常使用在GNSS導航處理。ECEF坐標到本地坐標轉換的旋轉矩陣表示為【圖9】。圖 9 Er旋轉矩陣參數說明：接收機位置的緯度：接收機位置的經度通過使用Er和接收機的坐標rr【ECEF】，坐標recef【ECEF】可以被轉換到本地坐標的坐標rlocal，公式如【圖10】。圖 102 RTKPPP定位算法2.1 單點定位（pntpos）1:satposs à 2: estpos à 3: estvel1. 計算計算衛星位置、速度和時鐘（satposs）a) 通過廣播星歷計算衛星鐘差（ephclk

5、）Ø 衛星編號到衛星系統的轉換（satsys）根據衛星的編號，獲取到對應的衛星導航系統。Ø 選擇星歷（seleph）1、傳入信號傳輸時間，衛星編號，導航數據等參數。2、遍歷導航數據，遍歷導航數據里面的星歷數據，判斷星歷數據的衛星編號是否和傳入的衛星編號相等。3、如果星歷數據的衛星編號和傳入的衛星編號相等，就計算星歷參考時間(toe)和信號傳輸時間的時間差。如果不相等，繼續處理下一條星歷數據。4、判斷計算出來的時間差，如果時間差大于了允許的最大時間差，繼續查找下一個星歷數據。否則，判斷時間差最小的星歷數據，記錄星歷數據的位置。5、返回之前記錄出來的星歷數據。Ø 使用

6、廣播星歷計算衛星時鐘偏差（eph2clk）1、 傳入信號發射時刻的時間和星歷數據。2、 計算信號發射時刻的時間和本時段鐘差參數參考時間(星歷參數toc)的時間差。3、 通過下式計算鐘差，這里還沒有處理相對論校正項和tgd：（代碼中有個迭代過程，資料上沒看到寫）衛星鐘差計算出來之后，信號發射時刻的時間還要加上這個鐘差。b) 計算衛星在信號發射時刻的位置、速度和時鐘（satpos）根據星歷表選項來選擇不同的處理，如下：Ø 廣播星歷（EPHOPT_BRDC）: ephpos （廣播星歷到衛星位置和鐘差）1、 根據公式計算出tk；2、 根據使用的衛星系統，選擇使用的地球引力常數(mu)和地球

7、的角速度(omge)3、 根據公式計算出平近點角M。4、 求解開普勒方程，按照以下公式迭代求解。5、 根據以下公式計算出u(改正后的緯度幅角), r(改正后的徑向), i(改正后的軌道傾角) 的值。6、 根據以下公式計算衛星在軌道平面內的坐標。7、 根據不同的衛星系統，做不同的計算。GPS計算方式：北斗計算方式：其中8、 按照公式計算出時間tc。9、 按照以下公式計算出鐘差和鐘漂。Ø 精密星歷（EPHOPT_PREC）: peph2posØ 廣播+SBAS（EPHOPT_SBAS）：satpos_sbasØ 廣播+SSR_APC（EPHOPT_SSRAPC）: s

8、atpos_ssrØ 廣播+SSR_COM（EPHOPT_SSRCOM）：satpos_ssrØ QZSS LEX星歷（EPHOPT_LEX）：lexeph2pos2. 使用偽距估算接收機的位置,返回估算狀態結果 （estpos）a) 偽距殘差（rescode）Ø 把ecef坐標系轉換成大地坐標系（ecef2pos）1、 按照以下公式做轉換，暫時還沒看懂。Ø 計算幾何距離和接收機到衛星的單位矢量（geodist）1、 用衛星的坐標向量做歐幾里德范數，返回值和地球長半軸(WGS84)比較。小于地球長半軸(WGS84)，返回-1；2、 計算衛星坐標和接收機坐

9、標的差值向量。3、 用差值向量做歐幾里德范數，再用上一步計算出來的差值向量和計算結果做除法，得到視線向量。按照如下公式：4、 使用以下公式計算幾何距離。Ø 計算衛星方位角/仰角（satazel）1、 把接收機ecef坐標轉換到大地坐標；2、 判斷高度是否大于地球半長軸（WGS84）的負數值；3、 如果高度小于等于地球半長軸（WGS84）的負數值，方位角為0，仰角為PI/2;4、 如果高度大于地球半長軸（WGS84）的負數值，把ECEF向量轉換到局部坐標。然后對轉換出來的坐標做內積。5、 使用如下公式計算衛星方位角和仰角。Ø 偽距使用編碼殘差改正（prange）暫時沒找到對應

10、的文檔對應。Ø 電離層改正（ionocorr）n 通過廣播電離層模型（klobuchar模型）計算出電離層延遲（ionmodel）1、 校驗傳入的電離層模型參數，校驗失敗，使用默認的電離層模型參數；2、 使用以下公式計算出地球為中心的角度（半圓）；3、 使用以下公式計算子的電離層的緯度/經度（半圓）；4、 使用以下公式計算地磁緯度。5、 計算本地時間，返回值按這個公式【tt-=floor(tt/86400.0)*86400.0;】處理，保證tt的范圍（0<=tt<86400）。6、 計算傾斜因子。7、 電離層延遲計算，公式如下。Ø 對流層改正（tropcorr）

11、n 通過標準大氣壓和saastamoinen模型計算對流層延遲 （tropmodel）1、 使用以下公式計算總氣壓。2、 使用以下公式計算絕對溫度。3、 使用以下公式計算水蒸汽的分壓。代碼中沒有使用。4、 使用以下公式計算【Saastamoinen模型】。Ø 偽距殘差殘差值=偽距-(物理距離+dtr-光速*時鐘偏差+電離層誤差+對流層誤差)Ø 時鐘系統和接收器的偏置補償Ø 偽距測量誤差方差（varerr）b) 方差權重值 (weight by variance)c) 最小二乘估計 （lsq）Ø 最小二乘估計通過求解正規方程（X =（A * A'）

12、 - 1* A * Y）1、 計算矩陣A*Y的結果，保存到矩陣Ay；2、 計算矩陣A*A的結果，保存到矩陣Q；3、 求的矩陣Q的逆矩陣,結果保存到矩陣Q；4、 最后x=Q*Ayd) 歐幾里德范數（norm）公式：1、 通過最小二乘法估算出來的參數dx【長度為4】，分別加到位置向量上x【長度為4】。2、 對這個速度向量做歐幾里德范數，返回值同1E-4比較，小于這個值時，就得到估算出的接收機的位置為向量x的值。e) 驗證求解 (valsol)1、 對偽距殘差值做內積，然后同卡方分布的自由度分布值(alpha=0.001)做比較，大于卡方分布的自由度分布值的數據，無效。2、 計算dops； 暫時還沒

13、沒找到資料。3、 用計算出來的dops和配置的最大dops閾值比較，大于配置的最大dops閾值的數據無效。3. 使用多普勒估算接收機速度（estvel）a) 多普勒殘差（resdop）Ø 把ecef坐標系轉換成大地坐標系（ecef2pos）Ø 把xyz坐標轉換成enu坐標（xyz2enu）Ø 計算ecef中的瞄準線向量Ø 計算相對于接收機在ECEF中的衛星速度Ø 多普勒殘差計算1、速率公式：b) 最小二乘估計 （lsq）Ø 最小二乘估計通過求解正規方程（X =（A * A'） - 1* A * Y）1、 計算矩陣A*Y的結果，

14、保存到矩陣Ay；2、 計算矩陣A*A的結果，保存到矩陣Q；3、 求的矩陣Q的逆矩陣,結果保存到矩陣Q；4、 最后x=Q*Ayc) 歐幾里德范數（norm）公式：3、 通過最小二乘法估算出來的參數dx【長度為4】，分別加到速度向量上x【長度為4】。4、 對這個速度向量做歐幾里德范數，返回值同1E-6比較，小于這個值時，就得到估算出的接收機速度為向量x的值。2.2 精確定位 （pppos）1:udstate_ppp à 2: satposs à 3: testeclipse à 4: res_ppp à 5: res_ppp à 6: filter

15、 à 7: res_ppp1. 暫時更新狀態（udstate_ppp）a) 位置更新b) 時鐘更新c) 對流層參數更新d) 相位偏差更新1、 通過LLI檢測周跳；2、 通過Geometry-Free相位跳變檢測周跳(如果雙頻測量值可用)；2. 計算衛星的位置和速度和時鐘（satposs）c) 通過廣播星歷計算衛星鐘差（ephclk）Ø 衛星編號到衛星系統的轉換（satsys）根據衛星的編號，獲取到對應的衛星導航系統。Ø 選擇星歷（seleph）1、傳入信號傳輸時間，衛星編號，導航數據等參數。2、遍歷導航數據，遍歷導航數據里面的星歷數據，判斷星歷數據的衛星編號是否和

16、傳入的衛星編號相等。3、如果星歷數據的衛星編號和傳入的衛星編號相等，就計算星歷參考時間(toe)和信號傳輸時間的時間差。如果不相等，繼續處理下一條星歷數據。4、判斷計算出來的時間差，如果時間差大于了允許的最大時間差，繼續查找下一個星歷數據。否則，判斷時間差最小的星歷數據，記錄星歷數據的位置。5、返回之前記錄出來的星歷數據。Ø 使用廣播星歷計算衛星時鐘偏差（eph2clk）1、 傳入信號傳輸時間和星歷數據。2、 傳入信號發射時刻的時間和星歷數據。3、 計算信號發射時刻的時間和本時段鐘差參數參考時間(星歷參數toc)的時間差。4、 通過下式計算鐘差，這里還沒有處理相對論校正項和tgd：（

17、代碼中有個迭代過程，資料上沒看到寫）d) 計算衛星位置、速度和時鐘（satpos）根據星歷表選項來選擇不同的處理，如下：Ø 廣播星歷（EPHOPT_BRDC）: ephpos （廣播星歷到衛星位置和鐘差）1、 根據公式計算出tk；2、 根據使用的衛星系統，選擇使用的地球引力常數(mu)和地球的角速度(omge)3、 根據公式計算出平近點角M。4、 求解開普勒方程，按照以下公式迭代求解。5、 根據以下公式計算出u(改正后的緯度幅角), r(改正后的徑向), i(改正后的軌道傾角) 的值。6、 根據以下公式計算衛星在軌道平面內的坐標。7、 根據不同的衛星系統，做不同的計算。GPS計算方式

18、：北斗計算方式：其中8、 按照公式計算出時間tc。9、 按照以下公式計算出鐘差和鐘漂。Ø 精密星歷（EPHOPT_PREC）: peph2posØ 廣播+SBAS（EPHOPT_SBAS）：satpos_sbasØ 廣播+SSR_APC（EPHOPT_SSRAPC）: satpos_ssrØ 廣播+SSR_COM（EPHOPT_SSRCOM）：satpos_ssrØ QZSS LEX星歷（EPHOPT_LEX）：lexeph2pos3. 排除對被遮蔽衛星的測量（testeclipse）a) 獲得太陽和月亮在ecef中的位置（sunmoonpos

19、）b) 歸一化三維矢量（normv3）c) 歐幾里德范數的向量（norm）d) 計算太陽-地球-衛星角（dot）4. 相位和編碼殘差（res_ppp）a) 地球潮汐改正（tidedisp）b) 計算幾何距離（geodist）1、 用衛星的坐標向量做歐幾里德范數，返回值和地球長半軸(WGS84)比較。小于地球長半軸(WGS84)，返回-1；2、 計算衛星坐標和接收機坐標的差值向量。3、 用差值向量做歐幾里德范數，再用上一步計算出來的差值向量和計算結果做除法，得到視線向量。4、 使用以下公式計算幾何距離。備注：這個公式是地球自轉改正公式。地球自轉改正產生地球自轉改正的機制在于地固系隨地球的自轉而旋

20、轉, 地固系是非慣性系統, 用地固系中的坐標計算衛星到接收機幾何距離時, 由于衛星位置和接收機位置是兩個不同時刻的位置矢量, 而且這兩個時刻的地固系相對于慣性系是變化的, 因此要考慮地球自轉引起的距離變化。其改正公式為: 其中, w 為地球自轉角速度; C 為光速; x si和x ri分別代表衛星位置矢量和測站位置矢量的分量, i= 1, 2 對應x , y 分量。地球自轉對緯度影響甚小, 對經度影響最大, 其次是高度。同時地球的自轉影響還跟測站的經緯度以及測站與衛星的幾何關系有關, 對兩極的測站, 影響為零; 對赤道上的測站影響最大; 當衛星在測站子午面內影響為零;衛星在測站東方時, 影響為

21、負, 衛星在測站西方時,影響為正。c) 衛星方位角/仰角計算（satazel）1、 把接收機ecef坐標轉換到大地坐標；2、 判斷高度是否大于地球半長軸（WGS84）的負數值；3、 如果高度小于等于地球半長軸（WGS84）的負數值，方位角為0，仰角為PI/2;4、 如果高度大于地球半長軸（WGS84）的負數值，把ECEF向量轉換到局部坐標。然后對轉換出來的坐標做內積。5、 使用如下公式計算衛星方位角和仰角。d) 對流層延遲修正：Ø 通過標準大氣壓和saastamoinen模型計算對流層延遲 （tropmodel）1、 使用以下公式計算總氣壓。2、 使用以下公式計算絕對溫度。3、 使用

22、以下公式計算水蒸汽的分壓。代碼中沒有使用。4、 使用以下公式計算【Saastamoinen模型】。e) 衛星天線模型（satantpcv）f) 接收機天線模型，通過天線相位中心參數計算天線偏移 （antmodel）g) 相位纏繞校正（windupcorr）h) 電離層和天線相位校正測量（corrmeas）i) 計算衛星時鐘和對流層延遲5. 卡爾曼濾波（filter）卡爾曼濾波按照以下的公式更新狀態：K = P * H * (H' * P * H + R)-1 ;xp = x + K * v ;Pp = (I K * H') * P ;參數解釋：x: 狀態向量 (n x 1)p:

23、 狀態的協方差矩陣 (n x n)H: 設計矩陣的轉置矩陣 (n x m)v: 創新(測量模型) (m x 1)R: 測量誤差的協方差矩陣 (m x m)n,m: 狀態和測量值的數量xp：狀態更新后的向量 (n x 1)Pp: 狀態更新后的協方差矩陣 (n x n)6. 組合殘差 （res_ppp）7. 解析PPP整周模糊度（pppamb）a) average LC（average_LC）b) 固定寬巷模糊（fix_amb_WL）Ø 計算LC波長度（m）（lam_LC）Ø 寬巷模糊Ø 寬巷模糊度的方差c) fix narrow-lane ambiguity 【AR

24、 mode: PPP-AR】d) fix narrow-lane ambiguity by ILS 【AR mode: PPP-AR ILS】3 以下是上面一些算法的參考資料，信息來自于Rtklib的手冊 3.1 接收機和衛星天線之間的幾何距離對應rtklib中的函數extern double geodist(const double *rs, const double *rr, double *e)3.2 衛星方向的方位角和仰角對應rtklib的函數extern double satazel(const double *pos, const double *e, double *azel)3

25、.3 GPS、Galileo和QNSS廣播星歷和時鐘對應rtklib的函數extern void eph2pos(gtime_t time, const eph_t *eph, double *rs, double *dts, double *var)3.4 北斗廣播星歷和時鐘對應rtklib的函數extern void eph2pos(gtime_t time, const eph_t *eph, double *rs, double *dts, double *var) 3.5 對流層和電流層模型3.5.1 對流層模型對應rtklib的函數extern double tropmodel(g

26、time_t time, const double *pos, const double *azel,double humi)3.5.2 電離層模型對應rtklib的函數extern double ionmodel(gtime_t t, const double *ion, const double *pos,const double *azel) 3.6 單點定位3.6.1 線性最小二乘估計對應rtklib的函數extern int lsq(const double *A, const double *y, int n, int m, double *x,double *Q)3.6.2 估算接收機的位置和鐘差stat