RTKLIB源碼解析（一）——單點定位(pntpos.c)

RTKLIB源碼解析（一）——單點定位(pntpos.c)

標籤： GNSS RTKLIB 單點定位

前段時間一直忙着寫畢業論文，因此也沒有太多時間來閱讀 RTKLIB源碼，最近好歹把 pntpos中的相關代碼看了一遍，知道了 RTKLIB是如何實現單點僞距定位的。這裏把每個函數都作成了小卡片的形式，每一個函數大都包含函數簽名、所在文件、功能說明、參數說明、處理過程、注意事項和個人疑惑這幾個部分，介紹了閱讀代碼時我本身的見解和疑惑。因此但願諸位看官能幫忙解答個人疑惑，與我交流，也但願能幫助後來也有須要閱讀 RTKLIB源碼的人，給他們多提供一份思路。總而言之，既爲人，也爲己。
這份文檔是使用 Cmd Markdown完成的，在做業部落上其格式顯式的很是完整，可是在博客園中目錄、代碼塊和流程圖彷佛都沒有顯示出來，因此這裏也貼上本文在做業部落上的連接 RTKLIB源碼解析（一）——單點定位(pntpos.c)，對格式「零容忍」的同窗請移步那裏。

目錄

• RTKLIB源碼解析（一）——單點定位(pntpos.c)
  • pntpos
  • satposs
  • estpos
  • raim_fde
  • estvel
  • ephclk
  • satpos
  • satsys
  • seleph
  • eph2clk
  • ephpos
  • eph2pos
  • rescode
  • lsq
  • valsol
  • matmul
  • dops
  • ecef2enu
  • xyz2enu
  • ecef2pos
  • geodist
  • satazel
  • prange
  • satexclude
  • ionocorr
  • tropcorr
  • varerr
  • testsnr
  • gettgd
  • ionmodel
  • iontec
  • iondelay
  • ionppp
  • interptec
  • tropmodel
  • resdop

---

pntpos

int pntpos (const obsd_t *obs, int n, const nav_t *nav, const prcopt_t *opt, sol_t *sol,
            double *azel, ssat_t *ssat, char *msg)

• 所在文件：pntpos.c
• 功能說明：依靠多普勒頻移測量值和僞距來進行單點定位，給出接收機的位置、速度和鐘差
• 參數說明：

函數參數，8個：
obsd_t     *obs      I     observation data
int        n         I     number of observation data
nav_t      *nav      I     navigation data
prcopt_t   *opt      I     processing options
sol_t      *sol      IO    solution
double     *azel     IO    azimuth/elevation angle (rad) (NULL: no output)
ssat_t     *ssat     IO    satellite status              (NULL: no output)
char       *msg      O     error message for error exit
返回類型:
int                  O     (1:ok,0:error)

• 調用關係：
  如無特別說明，本文所出現的流程圖中，縱向表明時間上的前後調用順序，橫向表明層次上的逐級調用順序。

  st=>start: pntpos
  satposs_=>operation: satposs
  estpos_=>operation: estpos
  raim_fde_=>operation: raim_fde
  estvel_=>operation: estvel
  e=>end: end
  
  st->satposs_->estpos_->raim_fde_->estvel_->e

  
  

pntpos satposs estpos raim_fde estvel

• 處理過程：

1. 檢查衛星個數是否大於 0
2. 當處理選項 opt中的模式不是單點模式時，電離層校訂採用 Klobuchar模型，對流層校訂則採用 Saastamoinen模型；相反，當其爲單點模式時，對輸入參數 opt不作修改。
3. 調用 satposs函數，按照所觀測到的衛星順序計算出每顆衛星的位置、速度、{鐘差、頻漂}。
4. 調用 estpos函數，經過僞距實現絕對定位，計算出接收機的位置和鐘差，順帶返回實現定位後每顆衛星的{方位角、仰角}、定位時有效性、定位後僞距殘差。
5. 調用 raim_fde函數，對上一步獲得的定位結果進行接收機自主正直性檢測（RAIM）。經過再次使用 vsat數組，這裏只會在對定位結果有貢獻的衛星數據進行檢測。
6. 調用 estvel函數，依靠多普勒頻移測量值計算接收機的速度。這裏只使用經過了上一步 RAIM_FDE操做的衛星數據，因此對於計算出的速度就沒有再次進行 RAIM了。
7. 首先將 ssat_t結構體數組的 vs(定位時有效性)、azel（方位角、仰角）、resp(僞距殘餘)、resc(載波相位殘餘)和 snr(信號強度)都置爲 0，而後再將實現定位後的 azel、snr賦予 ssat_t結構體數組，而 vs、resp則只賦值給那些對定位有貢獻的衛星，沒有參與定位的衛星，這兩個屬性值爲 0。

• 注意事項：

1. 這裏只計算了接收機的鐘差，而沒有計算接收機的頻漂，緣由在於 estvel函數中雖然計算獲得了接收機頻漂，但並無將其輸出到 sol_t:dtr中。
2. C語言中用 malloc申請的內存須要本身調用 free來予以回收，源碼中的 mat、imat、zeros等函數都只是申請了內存，並無進行內存的回收，在使用這些函數時，用戶必須本身調用 free來回收內存！源碼中將使用這些函數的代碼放置在同一行，在調用函數結尾處也統一進行內存回收，位置較爲明顯，不致於輕易忘記。

• 個人疑惑：

1. 源碼中將 obs[0].time做爲星曆選擇時間傳遞給 satposs函數，這樣對於每一顆觀測衛星，都要使用第一顆觀測衛星的數據接收時間做爲選擇星曆的時間標準。是否應該每顆衛星都使用本身的觀測時間？或者應該使用每顆衛星本身的信號發射時間？，仍是說這點差異對選擇合適的星曆其實沒有關係？
2. 這裏規定可以執行 raim_fde函數的前提是數目大於等於 6，感受不是隻要大於等於 5就能夠了嗎？

---

satposs

void satposs(gtime_t teph, const obsd_t *obs, int n, const nav_t *nav,
             int ephopt, double *rs, double *dts, double *var, int *svh)

• 所在文件：ephemeris.c
• 功能說明：按照所觀測到的衛星順序計算出每顆衛星的位置、速度、{鐘差、頻漂}
• 參數說明：

函數參數，9個：
gtime_t  teph      I   time to select ephemeris (gpst)
obsd_t   *obs      I   observation data
int      n         I   number of observation data
nav_t    *nav      I   navigation data
int      ephopt    I   ephemeris option (EPHOPT_???)
double   *rs       O   satellite positions and velocities，長度爲6*n，{x,y,z,vx,vy,vz}(ecef)(m,m/s)
double   *dts      O   satellite clocks，長度爲2*n， {bias,drift} (s|s/s)
double   *var      O   sat position and clock error variances (m^2)
int      *svh      O   sat health flag (-1:correction not available)
返回類型:
none

• 調用關係：

st_=>start: satposs
ephclk1=>operation: ephclk
satpos_=>operation: satpos
cond=>condition: 鐘差爲0？
ephclk2=>operation: ephclk
e=>end: end

st_->ephclk1->satpos_->cond
cond(no)->e
cond(yes,rigth)->ephclk2->e

satposs ephclk satpos

• 處理過程：

1. 按照觀測數據的順序，首先將將當前觀測衛星的 rs、dts、var和svh數組的元素置 0。
2. 經過判斷某一頻率下信號的僞距是否爲 0，來獲得此時所用的頻率個數。注意，頻率個數不能大於 NFREQ（默認爲 3）。
3. 用數據接收時間減去僞距信號傳播時間，獲得衛星信號的發射時間。
4. 調用 ephclk函數，由廣播星曆計算出當前觀測衛星的鐘差。注意，此時的鐘差是沒有考慮相對論效應和 TGD的。
5. 用 3中的信號發射時間減去 4中的鐘偏，獲得 GPS時間下的衛星信號發射時間。
6. 調用 satpos函數，計算信號發射時刻衛星的 P(ecef,m)、V(ecef,m/s)、C((s|s/s))。注意，這裏計算出的鐘差是考慮了相對論效應的了，只是尚未考慮 TGD。
7. 若是由 6中計算出的鐘偏爲 0，就再次調用 ephclk函數，將其計算出的衛星鐘偏做爲最終的結果。

• 注意事項：

1. ephclk函數計算的鐘偏是沒有考慮相對論和 TGD的，而 satpos函數考慮了相對論，沒有考慮 TGD。

• 個人疑惑：

1. 對於處理過程當中的第3步，數據接收時間減去僞距信號傳播時間，這裏的數據接收時間應該是有接收機獲得的，自己應該是包含接收機鐘差的，因此最終獲得的信號發射時間應該也是不許確的。難道說接收機鐘差較小，在此時能夠忽略不計？
2. ephclk函數最終是經過調用 eph2clk來計算衛星鐘偏的，但對於後者，我有問題。見 eph2clk處個人疑惑
3. 爲何要進行 7中操做？

---

estpos

int estpos(const obsd_t *obs, int n, const double *rs, const double *dts,
           const double *vare, const int *svh, const nav_t *nav,
           const prcopt_t *opt, sol_t *sol, double *azel, int *vsat,
           double *resp, char *msg)

• 所在文件：pntpos.c
• 功能說明：經過僞距實現絕對定位，計算出接收機的位置和鐘差，順帶返回實現定位後每顆衛星的{方位角、仰角}、定位時有效性、定位後僞距殘差。
• 參數說明：

函數參數，13個：
obsd_t   *obs      I   observation data
int      n         I   number of observation data
double   *rs       I   satellite positions and velocities，長度爲6*n，{x,y,z,vx,vy,vz}(ecef)(m,m/s)
double   *dts      I   satellite clocks，長度爲2*n， {bias,drift} (s|s/s)
double   *vare     I   sat position and clock error variances (m^2)
int      *svh      I   sat health flag (-1:correction not available)
nav_t    *nav      I   navigation data
prcopt_t *opt      I   processing options
sol_t    *sol      IO  solution
double   *azel     IO  azimuth/elevation angle (rad)
int      *vsat     IO  表徵衛星在定位時是否有效
double   *resp     IO  定位後僞距殘差 (P-(r+c*dtr-c*dts+I+T))
char     *msg      O   error message for error exit
返回類型:
int                O     (1:ok,0:error)

• 調用關係：

st=>start: estpos
rescode_=>operation: rescode
lsq_=>operation: lsq
valsol_=>operation: valsol
e=>end: end

st->rescode_->lsq_->valsol_->e

estpos rescode lsq valsol

• 處理過程：

1. 將 sol->rr的前 3項賦值給 x數組。
2. 開始迭代定位計算，首先調用 rescode函數，計算在當前接收機位置和鐘差值的狀況下，定位方程的右端部分 v(nv\*1)、幾何矩陣 H(NX*nv)、此時所得的僞距殘餘的方差 var、全部觀測衛星的 azel{方位角、仰角}、定位時有效性 vsat、定位後僞距殘差 resp、參與定位的衛星個數 ns和方程個數 nv。
3. 肯定方程組中方程的個數要大於未知數的個數。
4. 以僞距殘餘的標準差的倒數做爲權重，對 H和 v分別左乘權重對角陣，獲得加權以後的 H和 v。
5. 調用 lsq函數，根據 $\Delta x=(HH^T){-1}Hv$和 $Q=(HH^T){-1}$，獲得當前 x的修改量和定位偏差協方差矩陣中的權係數陣。
6. 將 5中求得的 x加入到當前 x值中，獲得更新以後的 x值。
7. 若是 5中求得的修改量小於截斷因子(目前是1e-4)，則將 6中獲得的 x值做爲最終的定位結果，對 sol的相應參數賦值，以後再調用 valsol函數確認當前解是否符合要求（僞距殘餘小於某個 $\chi^2$值和 GDOP小於某個門限值）。不然，進行下一次循環。
8. 若是超過了規定的循環次數，則輸出發散信息後，返回 0。

• 注意事項：

1. 關於第 1步，若是是第一次定位，即輸入的 sol爲空，則 x初值爲 0；若是以前有過定位，則經過 1中操做能夠將上一曆元的定位值做爲該曆元定位的初始值。
2. 關於定位方程，書中的寫法以下：
   $$
   G \begin{bmatrix} \Delta x \ \Delta y \ \Delta z \ \delta t_u \end{bmatrix} = b
   $$
   其中，
   $$
   G = \begin{bmatrix}
   -e^1_x & -e^1_y & -e^1_z & 1 \
   -e^2_x & -e^2_y & -e^2_z & 1 \
   \cdots & \cdots & \cdots & \cdots \
   -e^N_x & -e^N_y & -e^N_z & 1
   \end{bmatrix}
   \qquad
   b = \begin{bmatrix}
   \rho_r^1-(r_r1+c·dt_r-c·dt_s^1+I1+T^1) \
   \rho_r^2-(r_r2+c·dt_r-c·dt_s^2+I2+T^2) \
   \cdots \
   \rho_r^N-(r_rN+c·dt_r-c·dt_s^N+IN+T^N)
   \end{bmatrix}
   $$
   而 RTKLIB中採用的方程是下面這樣的（這裏先假設未知數的個數爲 4）：
   $$
   H^T \begin{bmatrix} \Delta x \ \Delta y \ \Delta z \ \delta t_u \end{bmatrix} = b
   $$
   其中，$H=G^T$。
   
   
   
   
   
   
   
   
   
   
   
   
   
   
   
   
   
   
   
   
   
   
   
   
3. 關於加權最小二乘，這裏的權重值是對角陣，這是創建在假設不一樣測量值的偏差之間是彼此獨立的；另外，這個權重值並不單是距測量偏差的，而是方程右端 b總體的測量偏差。最後，大部分資料上這裏都是把權重矩陣 W保留到方程的解的表達式當中，而這裏是直接對 H和 v分別左乘權重對角陣，獲得加權以後的 H和 v，其表示形式像是沒有加權同樣。
4. 若是某次迭代過程當中步長小於門限值(1e-4)，但經 valsol函數檢驗後該解無效，則會直接返回 0，並不會再進行下一次迭代計算。
5. 由於該函數有兩個返回路徑，並且又調用了 mat函數來構建矩陣，因此在兩個返回路徑處都須要調用 free函數來回收內存。
6. 源碼中的 dtr實際上單位是 m，是乘以了光速以後的。
7. 在對 sol結構體賦值時，並無直接將接收機鐘差 dtr賦值給 sol_dtr，而是經過在 sol中存儲的是減去接收機鐘差後的信號觀測時間，來說該信息包括到 sol中了。
8. 源碼中定位方程的個數 nv要大於有效觀測衛星的個數 ns，這裏爲了防止虧秩，彷佛又加了 3個未知數和觀測方程。
9. 在每一次從新調用 rescode函數時，其內部並無對 v、H和 var清零處理，因此當方程數變少時，可能會存在尾部仍保留上一次數據的狀況，可是由於數組相乘時都包含所需計算的長度，因此這種狀況並不會對計算結果形成影響。

• 個人疑惑：

$\color{red}{1. 這裏的 NX=7不明白，應該只有 4個未知數的啊！}$

---

raim_fde

int raim_fde(const obsd_t *obs, int n, const double *rs,
             const double *dts, const double *vare, const int *svh,
             const nav_t *nav, const prcopt_t *opt, sol_t *sol,
             double *azel, int *vsat, double *resp, char *msg)

• 所在文件：pntpos.c
• 功能說明：對計算獲得的定位結果進行接收機自主正直性檢測（RAIM）。經過再次使用 vsat數組，這裏只會在對定位結果有貢獻的衛星數據進行檢測。
• 參數說明：

函數參數，13個：
obsd_t   *obs      I   observation data
int      n         I   number of observation data
double   *rs       I   satellite positions and velocities，長度爲6*n，{x,y,z,vx,vy,vz}(ecef)(m,m/s)
double   *dts      I   satellite clocks，長度爲2*n， {bias,drift} (s|s/s)
double   *vare     I   sat position and clock error variances (m^2)
int      *svh      I   sat health flag (-1:correction not available)
nav_t    *nav      I   navigation data
prcopt_t *opt      I   processing options
sol_t    *sol      IO  solution
double   *azel     IO  azimuth/elevation angle (rad)
int      *vsat     IO  表徵衛星在定位時是否有效
double   *resp     IO  定位後僞距殘差 (P-(r+c*dtr-c*dts+I+T))
char     *msg      O   error message for error exit
返回類型:
int                O     (1:ok,0:error)

• 調用關係：

st=>start: raim_fde
estpos_=>operation: estpos
e=>end: end

st->estpos_->e

raim_fde estpos

• 處理過程：

1. 關於觀測衛星數目的循環，每次捨棄一顆衛星，計算使用餘下衛星進行定位的定位值。
2. 捨棄一顆衛星後，將剩下衛星的數據複製到一塊兒，調用 estpos函數，計算使用餘下衛星進行定位的定位值。
3. 累加使用當前衛星實現定位後的僞距殘差平方和與可用微信數目，若是 nvsat<5，則說明當前衛星數目過少，沒法進行 RAIM_FDE操做。
4. 計算僞距殘差平方和的標準平均值，若是大於 rms，則說明當前定位結果更合理，將 stat置爲 1，從新更新 sol、azel、vsat(當前被捨棄的衛星，此值置爲0)、resp等值，並將當前的 rms_e更新到 `rms'中。
5. 繼續棄用下一顆衛星，重複 2-4操做。總而言之，將一樣是棄用一顆衛星條件下，僞距殘差標準平均值最小的組合所得的結果做爲最終的結果輸出。
6. 若是 stat不爲 0，則說明在棄用衛星的前提下有更好的解出現，輸出信息，指出棄用了哪科衛星。
7. 使用 free函數回收相應內存。

• 注意事項：

1. 使用了 mat和 malloc函數後要使用 free函數來回收內存。
2. 源碼中有不少關於 i、j、k的循環。其中，i表示最外面的大循環，每次將將第 i顆衛星捨棄不用，這是經過 if (j==i) continue實現的；j表示剩餘使用的衛星的循環，每次進行相應數據的賦值；k表示參與定位的衛星的循環，與 j一塊兒使用。

• 個人疑惑：

1. 這裏執行 RAIM至少要有 6顆可用衛星，但是我感受 5顆就夠了呀！

---

estvel

void estvel(const obsd_t *obs, int n, const double *rs, const double *dts,
            const nav_t *nav, const prcopt_t *opt, sol_t *sol,
            const double *azel, const int *vsat)

• 所在文件：pntpos.c
• 功能說明：依靠多普勒頻移測量值計算接收機的速度。
• 參數說明：

函數參數，9個：
obsd_t   *obs      I   observation data
int      n         I   number of observation data
double   *rs       I   satellite positions and velocities，長度爲6*n，{x,y,z,vx,vy,vz}(ecef)(m,m/s)
double   *dts      I   satellite clocks，長度爲2*n， {bias,drift} (s|s/s)
nav_t    *nav      I   navigation data
prcopt_t *opt      I   processing options
sol_t    *sol      IO  solution
double   *azel     IO  azimuth/elevation angle (rad)
int      *vsat     IO  表徵衛星在定位時是否有效
返回類型:
int                O     (1:ok,0:error)

• 調用關係：

st=>start: estvel
resdop_=>operation: resdop
lsq_=>operation: lsq
e=>end: end

st->resdop_->lsq_->e

estvel resdop lsq

• 處理過程：

1. 在最大迭代次數限制內，調用 resdop，計算定速方程組左邊的幾何矩陣和右端的速度殘餘，返回定速時所使用的衛星數目。
2. 調用 lsq函數，解出 {速度、頻漂}的步長，累加到 x中。
3. 檢查當前計算出的步長的絕對值是否小於 1E-6。是，則說明當前解已經很接近真實值了，將接收機三個方向上的速度存入到 sol_t.rr中；否，則進行下一次循環。
4. 執行完全部迭代次數，使用 free回收內存。

• 注意事項：

1. 最終向 sol_t類型存儲定速解時，並無存儲所計算出的接收器時鐘頻漂。
2. 這裏不像定位時，初始值可能爲上一曆元的位置(從 sol中讀取初始值)，定速的初始值直接給定爲 0.

---

ephclk

int ephclk(gtime_t time, gtime_t teph, int sat, const nav_t *nav, double *dts)

• 所在文件：ephemeris.c
• 功能說明：經過廣播星從來肯定衛星鐘偏
• 參數說明：

函數參數，5個：
gtime_t  time      I   transmission time by satellite clock
gtime_t  teph      I   time to select ephemeris (gpst)
int      sat       I   satellite number (1-MAXSAT)
nav_t    *nav      I   navigation data
double   *dts      O   satellite clocks，長度爲2*n， {bias,drift} (s|s/s)
返回類型:
int                O     (1:ok,0:error)

• 調用關係：

st=>start: ephclk
satsys_=>operation: satsys
seleph_=>operation: seleph
eph2clk_=>operation: eph2clk
e=>end: end

st->satsys_->seleph_->eph2clk_->e

ephclk satsys seleph eph2clk

• 處理過程：

1. 首先調用 satsys函數，根據衛星編號肯定該衛星所屬的導航系統和該衛星在該系統中的 PRN編號。
2. 對於 GPS導航系統，調用 seleph函數來選擇 toe值與星曆選擇時間標準 teph最近的那個星曆。
3. 調用 eph2clk函數，經過廣播星曆和信號發射時間計算出衛星鐘差。

• 注意事項：

1. 此時計算出的衛星鐘偏是沒有考慮相對論效應和 TGD的。
2. 目前我還只關心 RTKLIB對於 GPS導航所進行的數據處理，因此這裏在選擇星曆和計算鐘差時都只介紹與 GPS系統有關的函數。

• 個人疑惑：

1. 見 eph2clk處

---

satpos

int satpos(gtime_t time, gtime_t teph, int sat, int ephopt, const nav_t *nav, 
           double *rs, double *dts, double *var, int *svh)

• 所在文件：ephemeris.c
• 功能說明：計算信號發射時刻衛星的 P(ecef,m)、V(ecef,m/s)、C((s|s/s))
• 參數說明：

函數參數，9個：
gtime_t time      I   time (gpst)
gtime_t teph      I   time to select ephemeris (gpst)
int     sat       I   satellite number
nav_t   *nav      I   navigation data
int     ephopt    I   ephemeris option (EPHOPT_???)
double  *rs       O   sat position and velocity {x,y,z,vx,vy,vz} (ecef)(m|m/s)
double  *dts      O   sat clock {bias,drift} (s|s/s)
double  *var      O   sat position and clock error variance (m^2)
int     *svh      O   sat health flag (-1:correction not available)
返回類型:
int               O     (1:ok,0:error)

• 調用關係：

st=>start: satpos
ephpos_=>operation: ephpos
e=>end: end

st->ephpos_->e

satpos ephpos

• 處理過程：

1. 判斷星曆選項的值，若是是 EPHOPT_BRDC，調用 ephpos函數，根據廣播星曆計算出算信號發射時刻衛星的 P、V、C

• 注意事項：

1. 此時計算出的衛星鐘差考慮了相對論，尚未考慮 TGD。
2. 目前還只閱讀了如何從廣播星曆中計算衛星 P、V、C的代碼，關於如何從精密星曆中計算，等對精密星曆的理論背景有了更多瞭解以後再予以添加。

---

satsys

int satsys(int sat, int *prn)

• 所在文件：rtkcmn.c
• 功能說明：根據衛星編號肯定該衛星所屬的導航系統和該衛星在該系統中的 PRN編號
• 參數說明：

函數參數，2個：
int    sat       I   satellite number (1-MAXSAT)
int    *prn      IO  satellite prn/slot number (NULL: no output)
返回類型:
int         satellite system (SYS_GPS,SYS_GLO,...)

• 處理過程：

1. 爲處理意外狀況（衛星號不在 1-MAXSAT以內），先令衛星系統爲 SYS_NONE。
2. 按照 TRKLIB中定義相應宏的順序來判斷是不是 GPS、GLO、GAL系統等，判斷標準是將衛星號減去前面的導航系統所擁有的衛星個數，來判斷剩餘衛星個數是否小於等於本系統的衛星個數。
3. 肯定所屬的系統後，經過加上最小衛星編號的 PRN再減去 1，獲得該衛星在該系統中的 PRN編號。

• 注意事項：

1. 這裏的衛星號是從 1開始排序的，這也是不少函數中與之有關的數組在調用時形式寫爲 A[B.sat-1]

---

seleph

eph_t *seleph(gtime_t time, int sat, int iode, const nav_t *nav)

• 所在文件：ephemeris.c
• 功能說明：從導航數據中選擇星曆。當 iode>=0時，選擇與輸入期號相同的星曆；不然，選擇 toe值與星曆選擇時間標準 time最近的那個星曆。
• 參數說明：

函數參數，4個：
gtime_t  time      I   time to select ephemeris (gpst)
int      sat       I   satellite number (1-MAXSAT)
int      iode      I   星曆數據期號
nav_t    *nav      I   navigation data
返回類型:
eph_t *         星曆數據

• 處理過程：

1. 根據該衛星所屬的導航系統給與星曆參考時間的最大時間間隔 tmax賦予相應的值。
2. 遍歷導航數據，首先確保所查找星曆的衛星號是否相同，接着確保星曆期號是否相同。
3. 接着確保星曆選擇時間與代查找星曆的參考時間的間隔是否小於 tmax。
4. 對於經過了 2-3驗證的星曆，若是此時對於輸入的期號，有 iode>=0，則該星曆就是所要尋找的星曆；不然，驗證 3中的 t是否知足 t<=tmin。知足的話，就令 tmin=t，該星曆目前是距離參考時間最近的。
5. 循環 2-4步操做，遍歷完導航數據中的全部星曆。若是都沒有符合條件的，就輸出信息並返回 NULL；不然，返回所查找到的星曆。

• 注意事項：

1. 對於 GPS系統，星曆數據期號每 2h更新一次，因此 RTKLIB中對 MAXDTOE的定義爲 7200。
2. 若是存在兩個相鄰時間段的星曆，經過 4中令 tmin=t的操做能夠最終查找出參考時間距 time最近的那個星曆。

• 個人疑惑：

1. 爲何 tmax和tmin都要加 1？
2. IODE正常狀況下應該都是 >=0的，爲何還要考慮 <0的狀況？
3. 考慮到星曆中衛星的健康情況，這裏在選擇星曆時是否也應該驗證 eph.svh==0呢？

---

eph2clk

int eph2clk (gtime_t time, const eph_t *eph)

• 所在文件：ephemeris.c
• 功能說明：根據信號發射時間和廣播星曆，計算衛星鐘差
• 參數說明：

函數參數，2個
gtime_t time     I   time by satellite clock (gpst)
eph_t    *eph    I   broadcast ephemeris
返回類型：
double    satellite clock bias (s) without relativeity correction

• 處理過程：

1. 計算與星曆參考時間的誤差 dt = t-toc。
2. 利用二項式校訂計算出衛星鐘差，從 dt中減去這部分，而後再進行一次上述操做，獲得最終的 dt。
3. 使用二項式校訂獲得最終的鐘差。

- **個人疑惑**： > 1. **看不懂上述處理過程，跟以往資料上都不同。咋回事？**

---

ephpos

int ephpos(gtime_t time, gtime_t teph, int sat, const nav_t *nav,
           int iode, double *rs, double *dts, double *var, int *svh)

• 所在文件：ephemeris.c
• 功能說明：根據廣播星曆計算出算信號發射時刻衛星的 P、V、C
• 參數說明：

函數參數，9個
gtime_t  time      I   transmission time by satellite clock
gtime_t  teph      I   time to select ephemeris (gpst)
int      sat       I   satellite number (1-MAXSAT)
nav_t    *nav      I   navigation data
int      iode      I   星曆數據期號
double   *rs       O   satellite positions and velocities，長度爲6*n，{x,y,z,vx,vy,vz}(ecef)(m,m/s)
double   *dts      O   satellite clocks，長度爲2*n， {bias,drift} (s|s/s)
double   *var      O   sat position and clock error variances (m^2)
int      *svh      O   sat health flag (-1:correction not available)
返回類型：
int                O    (1:ok,0:error)

• 調用關係：

st=>start: ephpos
satsys_=>operation: satsys
seleph_=>operation: seleph
eph2pos_=>operation: eph2pos
e=>end: end

st->satsys_->seleph_->eph2pos_->e

ephpos satsys seleph eph2pos

• 處理過程：

1. 調用 satsys函數，肯定該衛星所屬的導航系統。
2. 若是該衛星屬於 GPS系統，則調用 seleph函數來選擇廣播星曆。
3. 根據選中的廣播星曆，調用 eph2pos函數來計算信號發射時刻衛星的 位置、鐘差和相應結果的偏差。
4. 在信號發射時刻的基礎上給定一個微小的時間間隔，再次計算新時刻的 P、V、C。與 3結合，經過擾動法計算出衛星的速度和頻漂。

• 注意事項：

1. 這裏是使用擾動法計算衛星的速度和頻漂，並無使用那些位置和鐘差公式對時間求導的結果。
2. 因爲是調用的 eph2pos函數，計算獲得的鐘差考慮了相對論效應，尚未考慮 TGD

---

eph2pos

void eph2pos(gtime_t time, const eph_t *eph, double *rs, double *dts, double *var)

• 所在文件：ephemeris.c
• 功能說明：根據廣播星曆計算出算信號發射時刻衛星的位置和鐘差
• 參數說明：

函數參數，5個
gtime_t  time      I   transmission time by satellite clock
eph_t    *eph      I   broadcast ephemeris
double   *rs       O   satellite positions and velocities，長度爲6*n，{x,y,z,vx,vy,vz}(ecef)(m,m/s)
double   *dts      O   satellite clocks，長度爲2*n， {bias,drift} (s|s/s)
double   *var      O   sat position and clock error variances (m^2)
返回類型：
none

• 處理過程：

1. 與大部分資料上計算衛星位置和鐘差的過程是同樣的，只是這裏在計算偏近點角 E時採用的是牛頓法來進行迭代求解。
2. 計算偏差直接採用 URA值來標定，具體對應關係可在 ICD-GPS-200C P83中找到。

• 注意事項：

1. 這裏在計算鐘差時，就與大部分資料上介紹的同樣了，只進行了一次二項式校訂。另外，這裏的鐘差考慮了相對論效應，並無考慮 TGD。
2. 在計算偏近點角 E時，這裏採用的是牛頓法來進行迭代求解，這裏與不少資料上可能不同，具體內容可在 RTKLIB manual P143中找到。
3. 補充幾個程序中的參數說明：

mu，     地球引力常數
eph->A， 軌道長半徑
omgea,   地球自轉角速度

rescode

int rescode(int iter, const obsd_t *obs, int n, const double *rs,
                   const double *dts, const double *vare, const int *svh,
                   const nav_t *nav, const double *x, const prcopt_t *opt,
                   double *v, double *H, double *var, double *azel, int *vsat,
                   double *resp, int *ns)

• 所在文件：pntpos.c
• 功能說明：計算在當前接收機位置和鐘差值的狀況下，定位方程的右端部分 v(nv\*1)、幾何矩陣 H(NX*nv)、此時所得的僞距殘餘的方差 var、全部觀測衛星的 azel{方位角、仰角}、定位時有效性 vsat、定位後僞距殘差 resp、參與定位的衛星個數 ns和方程個數 nv。
• 參數說明：

函數參數，17個
int      iter      I    迭代次數
obsd_t   *obs      I    observation data
int      n         I    number of observation data
double   *rs       I   satellite positions and velocities，長度爲6*n，{x,y,z,vx,vy,vz}(ecef)(m,m/s)
double   *dts      I   satellite clocks，長度爲2*n， {bias,drift} (s|s/s)
double   *vare     I   sat position and clock error variances (m^2)
int      *svh      I   sat health flag (-1:correction not available)
nav_t    *nav      I   navigation data
double   *x        I   本次迭代開始以前的定位值
prcopt_t *opt      I   processing options
double   *v        O   定位方程的右端部分，僞距殘餘
double   *H        O   定位方程中的幾何矩陣
double   *var      O   參與定位的僞距殘餘方差
double   *azel     O   對於當前定位值，每一顆觀測衛星的 {方位角、高度角}
int      *vsat     O   每一顆觀測衛星在當前定位時是否有效
double   *resp     O   每一顆觀測衛星的僞距殘餘， (P-(r+c*dtr-c*dts+I+T))
int      *ns       O   參與定位的衛星的個數
返回類型：
int                O   定位方程組的方程個數

• 調用關係：

st=>start: rescode
ecef2pos_=>operation: ecef2pos
satsys_=>operation: satsys
geodist_=>operation: geodist
satazel_=>operation: satazel
prange_=>operation: prange
satexclude_=>operation: satexclude
ionocorr_=>operation: ionocorr
tropcorr_=>operation: tropcorr
varerr_=>operation: varerr
e=>end: end

st->ecef2pos_->satsys_->geodist_->satazel_->prange_->satexclude_->ionocorr_->tropcorr_->varerr_->e

rescode ecef2pos satsys geodist satazel prange satexclude ionocorr tropcorr varerr

• 處理過程：

1. 將以前獲得的定位解信息賦值給 rr和dtr數組，以進行關於當前解的僞距殘餘的相關計算。
2. 調用 ecef2pos函數，將上一步中獲得的位置信息由 ECEF轉化爲大地座標系。
3. 將 vsat、azel和 resp數組置 0，由於在先後兩次定位結果中，每顆衛星的上述信息都會發生變化。
4. 調用 satsys函數，驗證衛星編號是否合理及其所屬的導航系統。
5. 檢測當前觀測衛星是否和下一個相鄰數據重複。是，則 i++後繼續下一次循環；否，則進入下一步。
6. 調用 geodist函數，計算衛星和當前接收機位置之間的幾何距離和 receiver-to-satellite方向的單位向量。而後檢驗幾何距離是否 >0。
7. 調用 satazel函數，計算在接收機位置處的站心座標系中衛星的方位角和仰角，檢驗仰角是否 $\geq$截斷值。
8. 調用 prange函數，獲得僞距值。
9. 能夠在處理選項中事先指定只選用哪些導航系統或衛星來進行定位，這是經過調用 satexclude函數完成的。
10. 調用 ionocorr函數，計算電離層延時(m)。
11. 上一步中所得的電離層延時是創建在 L1信號上的，當使用其它頻率信號時，依據所用信號頻組中第一個頻率的波長與 L1波長的關係，對上一步獲得的電離層延時進行修正。
12. 調用 tropcorr函數，計算對流層延時(m)。
13. 由 $\rho_r^i-(r_ri+dt_r-c·dt_s^i+Ii+T^i)$，計算出此時的僞距殘餘。
14. 組裝幾何矩陣，前 3行爲 6中計算獲得的視線單位向量的反向，第 4行爲 1，其它行爲 0。
15. 將參與定位的衛星的定位有效性標誌設爲 1，給當前衛星的僞距殘餘賦值，參與定位的衛星個數 ns加 1.
16. 調用 varerr函數，計算此時的導航系統偏差（可能會包括 IFLC選項時的電離層延時），而後累加計算用戶測距偏差(URE)。
17. 爲了防止虧秩，人爲的添加了幾組觀測方程。

• 注意事項：

1. 返回值 v和 resp的主要區別在於長度不一致， v是須要參與定位方程組的解算的，維度爲 nv1；而 resp僅表示全部觀測衛星的僞距殘餘，維度爲 n1，對於沒有參與定位的衛星，該值爲 0。
2. 源碼中 dtr的單位是 m。
3. 16中的 URE值包括 ①衛星星曆和鐘差的偏差 ②大氣延時偏差 ③僞距測量的碼偏移偏差 ④導航系統的偏差

• 個人疑惑：

1. 關於 5中的去除重複數據的過程，有如下幾個見解：
   ① 這樣作的前提是相同衛星的重複數據必須相鄰出現！
   ② 爲何在這裏要進行重複數據檢測，在構建 obsd_t結構體時就能夠進行這項工做呀？
   ③ 5中當數據重複時，i++後繼續下一次循環，這樣的話會直接略去後面所重複的數據，這樣作就會將兩個相鄰重複數據都忽略掉，就至關於重複數據都不使用。感受能夠用其中一個的啊！
   
   
   
2. 11步中，爲何要選擇所用信號頻組中第一個頻率的波長來進行電離層延時修正呢？還有，電離層延時的值發生了改變，那這裏的方差是否也須要跟着一塊兒改變呢？
3. 在計算電離/對流層延時時，均傳入了 iter>0?opt->ionoopt:IONOOPT_BRDC或 iter>0?opt->tropopt:TROPOPT_SAAS參數，都強調了當 iter==0時，會強制使用 Klobuchar或Saastamoinen模型。這會不會是由於這兩種模型都是屬於對接收機位置不是很是敏感的類型？
4. 這裏虧秩應該就是欠定方程的意思吧？這裏 17中的操做沒有看懂，也沒能找到相關理論依據

---

lsq

int lsq(const double *A, const double *y, int n, int m, double *x, double *Q)

• 所在文件：rtkcmn.c
• 功能說明：計算獲得方程 $x=(AA^T){-1}Ay$ 左邊的值和該值的協方差矩陣 $Q=(AA^T){-1}$。
• 參數說明：

函數參數，6個
double *A        I   transpose of (weighted) design matrix (n x m)
double *y        I   (weighted) measurements (m x 1)
int    n,m       I   number of parameters and measurements (n<=m)
double *x        O   estmated parameters (n x 1)
double *Q        O   esimated parameters covariance matrix (n x n)
返回類型：
int              O   (0:ok,0>:error)

• 調用關係：

st=>start: lsq
matmul_=>operation: matmul
matinv_=>operation: matinv
e=>end: end

st->matmul_->matinv_->e

lsq matmul

• 處理過程：

1. 首先計算右半部分 $A_y=Ay$。
2. 再計算左半部分括號裏面的值 $Q=AA^T$。
3. 計算 Q矩陣的逆 $Q^{-1}$，但仍存儲在 Q中，最後再右乘 $A_y$，獲得 x的值。

• 注意事項：

1. 對於加權最小二乘，能夠直接在調用該函數以前直接將 A、y進行加權處理，以後在調用該函數，這樣獲得的就是加權最小二乘的解。
2. 全部的矩陣都是列優先存儲的，對於整個源代碼來講，矩陣都是這樣存儲的。因此對於代碼中出現的一維矩陣，基本都應該是列向量。在閱讀數組下標時，記住這一點是很是重要的。
3. 矩陣求逆並不簡單，尤爲是對於接近奇異的矩陣。可是因爲這是個基本功能，並不打算繼續深刻下去。

---

valsol

int valsol(const double *azel, const int *vsat, int n,
           const prcopt_t *opt, const double *v, int nv, int nx, char *msg)

• 所在文件：pntpos.c
• 功能說明：確認當前解是否符合要求，即僞距殘餘小於某個 $\chi^2$值和 GDOP小於某個門限值）
• 參數說明：

函數參數，8個
double   *azel     I   azimuth/elevation angle (rad)
int      *vsat     I   表徵衛星在定位時是否有效
int      n         I   number of observation data
prcopt_t *opt      I   processing options
double   *v        I   定位後僞距殘差 (P-(r+c*dtr-c*dts+I+T))
int      nv        I   定位方程的方程個數
int      nx        I   未知數的個數
char     *msg      O   error message for error exit
返回類型：
int                O    (1:ok,0:error)

• 調用關係：

st=>start: valsol
dops_=>operation: dops
e=>end: end

st->dops_->e

valsol dops

• 處理過程：

1. 先計算定位後僞距殘餘平方加權和 vv。
2. 檢查是否知足 $vv>\chi^2(nv-nx-1)$。是，則說明此時的定位解偏差過大，返回 0；不然，轉到下一步。
3. 複製 azel，這裏只複製那些對於定位結果有貢獻的衛星的 zael值，而且統計實現定位時所用衛星的數目。
4. 調用 dops函數，計算各類精度因子(DOP)，檢驗是否有 0<GDOP<max。否，則說明該定位解的精度不符合要求，返回 0；是，則返回 1。

• 注意事項：

1. 關於 2中的 $\chi^2$檢驗，這裏這裏與 RTKLIB manual P160中不一致（但與教材一致），文檔中知足 $\frac{v^{Tv}{nv-nx-1}>\chi}2(nv-nx-1)$時就會不合格。與文檔中相比，這裏的寫法將會放寬對於位置解的檢驗。

---

matmul

源碼中定義了兩個 matmul函數，一個是在包含了 LAPACK/BLAS/MKL庫使用，調用其中的 degmn函數來完成矩陣相乘操做。這裏主要說明在沒有包含上述庫時自定義的矩陣相乘函數。

void matmul(const char *tr, int n, int k, int m, double alpha,
            const double *A, const double *B, double beta, double *C)

• 所在文件：rtkcmn.c
• 功能說明：可進行以下矩陣運算 C = alphaAB + beta*C，而且能經過 tr標誌來選擇是否對 A、B進行轉置
• 參數說明：

函數參數，6個
char     *tr       I   transpose flags ("N":normal,"T":transpose)
int      n,k,m     I   size of (transposed) matrix A,B
double   alpha     I   alpha
double   *A,*B     I   (transposed) matrix A (n x m), B (m x k)
double   beta      I   beta
double   *C        IO  matrix C (n x k)
返回類型：
none

• 處理過程：

1. 根據 tr的值肯定矩陣相乘標誌，即
   $$
   \begin{matrix}
   AB \rightarrow NN \rightarrow 1 \
   A
   B^T \rightarrow NT \rightarrow 2 \
   A^T*B \rightarrow TN \rightarrow 3 \
   A^T*BT \rightarrow TT \rightarrow 4
   \end{matrix}
   $$
   
   
   
   
   
   
   
2. 按照 f的值，分別執行相應的元素相乘並累加操做。
3. 對於 2中獲得的乘積 tmp，執行 C = alphatmp + betaC操做。

• 注意事項：

1. 在調用該函數時，要確保矩陣的維度是否和上述參數說明一致。
2. 全部的矩陣都是列優先存儲的，記住這一點，對於看明白 2中不一樣相乘方式時，元素如何相乘累加是相當重要的。
3. 矩陣求逆並不簡單，尤爲是對於接近奇異的矩陣。可是因爲這是個基本功能，並不打算繼續深刻下去。

---

dops

void dops(int ns, const double *azel, double elmin, double *dop)

• 所在文件：rtkcmn.c
• 功能說明：由站心座標系時的權係數矩陣表達式$\tilde{H} = (\tilde{G}\tilde{G}^T){-1}$，計算各類精度因子。
• 參數說明：

函數參數，4個
int      ns        I   number of satellites
double   *azel     I   satellite azimuth/elevation angle (rad)
double   elmin     I   elevation cutoff angle (rad)
double   *dop      O   DOPs {GDOP,PDOP,HDOP,VDOP}
返回類型：
none

• 調用關係：

st=>start: dops
matmul_=>operation: matmul
matinv_=>operation: matinv
e=>end: end

st->matmul_->matinv_->e

dops matmul matinv

• 處理過程：

1. 先按照以下站心座標系時的幾何矩陣 $\tilde{G}$的表達式求出其值。
   $$
   \tilde{G} = \begin{bmatrix}
   cos\theta^1sin\alpha1 & cos\theta^2sin\alpha2 & \vdots & cos\theta^Nsin\alphaN \
   cos\theta^1cos\alpha1 & cos\theta^2cos\alpha2 & \vdots & cos\theta^Ncos\alphaN \
   sin\theta^1 & sin\theta^2 & \vdots & sin\theta^N \
   1 & 1 & \vdots & 1
   \end{bmatrix}
   $$
   
   
   
   
   
   
   
   
2. 檢驗上述矩陣的列數是否$\geq$ 4。
3. 調用 matmul和matinv函數，計算出 $\tilde{H} = (\tilde{G}\tilde{G}^T){-1}$的值。
4. 若是能正確計算出逆矩陣，就按照順序計算出 GDOP、PDOP、HDOP和 VDOP的值。

• 個人疑惑：

1. 1中幾何矩陣 $\tilde{G}$與書中的不一致，前 3行均少了一個負號，我本身推導以後也以爲應該有負號，即爲
   $$
   \tilde{G} = \begin{bmatrix}
   -cos\theta^1sin\alpha1 & -cos\theta^2sin\alpha2 & \vdots & -cos\theta^Nsin\alphaN \
   -cos\theta^1cos\alpha1 & -cos\theta^2cos\alpha2 & \vdots & -cos\theta^Ncos\alphaN \
   -sin\theta^1 & -sin\theta^2 & \vdots & -sin\theta^N \
   1 & 1 & \vdots & 1
   \end{bmatrix}
   $$
   不過，因爲 $\tilde{H} = (\tilde{G}\tilde{G}^T){-1}$，在括號裏面的矩陣相乘時，是否有負號只對底邊靠左 3個元素和右邊靠上 3個元素有影響（多了個負號），而後再進行求逆以後，前 3行是否有負號就對對角線上的元素彷佛沒有影響了。
   
   
   
   
   
   
   
   
   
2. 感受在 rescode函數中，在檢驗一個觀測衛星的僞距信息是否可用時，已經進行過是否大於截斷高度角的檢測了。這裏所用的衛星仰角又都是屬於參與了定位的衛星，因此感受這裏應該不須要再進行一次高度角檢測了吧？

---

ecef2enu

void ecef2enu(const double *pos, const double *r, double *e)

• 所在文件：rtkcmn.c
• 功能說明：將 ECEF座標系(X、Y、Z)中的向量轉換成站心座標系。
• 參數說明：

函數參數，3個
double *pos      I   geodetic position {lat,lon} (rad)
double *r        I   vector in ecef coordinate {x,y,z}
double *e        O   vector in local tangental coordinate {e,n,u}
返回類型：
none

• 調用關係：

st=>start: ecef2enu
xyz2enu_=>operation: xyz2enu
matmul_=>operation: matmul
e=>end: end

st->xyz2enu_->matmul_->e

ecef2enu xyz2enu matmul

• 處理過程：

1. 先調用 xyz2enu函數計算此時的座標變換矩陣 E。
2. 調用 matmul計算 $E*r$的值，即爲目標值。

---

xyz2enu

void xyz2enu(const double *pos, double *E)

• 所在文件：rtkcmn.c
• 功能說明：計算將ECEF中的向量轉換到站心座標系中的轉換矩陣。
• 參數說明：

函數參數，2個
double *pos      I   geodetic position {lat,lon} (rad)
double *E        O   vector in local tangental coordinate {e,n,u}
返回類型：
none

• 處理過程：

1. 按照大部分資料上的寫法計算 3*3的矩陣，優先按列存儲。

---

ecef2pos

void ecef2pos(const double *r, double *pos)

• 所在文件：rtkcmn.c
• 功能說明：將 ECEF座標系(X、Y、Z)轉換成大地座標系($\lambda、\phi、h$)。
• 參數說明：

函數參數，2個
double  *r        I   ecef position {x,y,z} (m)
double  *pos      O   geodetic position {lat,lon,h} (rad,m)
返回類型：
none

• 處理過程：

這裏採用的方法與不少資料上的並不一致，而關於源碼中方法的具體理論推導和計算過程，見 ECEF和大地座標系的相互轉化一文。

---

geodist

double geodist(const double *rs, const double *rr, double *e)

• 所在文件：rtkcmn.c
• 功能說明：計算衛星和當前接收機位置之間的幾何距離和 receiver-to-satellite方向的單位向量。
• 參數說明：

函數參數，3個
double *rs       I   satellilte position (ecef at transmission) (m)
double *rr       I   receiver position (ecef at reception) (m)
double *e        O   line-of-sight unit vector (ecef)
返回類型：
double           O   geometric distance (m) (0>:error/no satellite position)

• 處理過程：

1. 檢查衛星到 WGS84座標系原點的距離是否大於基準橢球體的長半徑。
2. ps-pr，計算由接收機指向衛星方向的觀測矢量，以後在計算出相應的單位矢量。
3. 考慮到地球自轉，即信號發射時刻衛星的位置與信號接收時刻衛星的位置在 WGS84座標系中並不一致，進行關於 Sagnac效應的校訂。

• 個人疑惑：

1. 3中使用關於 Sagnac效應的校訂來考慮地球自轉對衛星位置的影響，與教材中的地球自轉校訂並不同，兩者是否描述的是同一個事情？

---

satazel

double satazel(const double *pos, const double *e, double *azel)

• 所在文件：rtkcmn.c
• 功能說明：計算在接收機位置處的站心座標系中衛星的方位角和仰角
• 參數說明：

函數參數，3個
double *pos      I   geodetic position {lat,lon,h} (rad,)
double *e        I   receiver-to-satellilte unit vevtor (ecef)
double *azel     IO  azimuth/elevation {az,el} (rad) (NULL: no output)  (0.0<=azel[0]<2*pi,-pi/2<=azel[1]<=pi/2)
返回類型：
double           O   elevation angle (rad)

• 調用關係：

st=>start: satazel
ecef2enu_=>operation: ecef2enu
e=>end: end

st->ecef2enu_->e

satazel ecef2enu

• 處理過程：

1. 調用 ecef2enu函數，將 pos處的向量轉換到以改點爲原點的站心座標系中。
2. 調用 atan2函數計算出方位角，而後再算出仰角。

• 注意事項：

1. 這裏在計算方位角時，要使用 atan2函數，而不能是 atan函數，詳細緣由見 C語言中的atan和atan2。

---

prange

double prange(const obsd_t *obs, const nav_t *nav, const double *azel,
              int iter, const prcopt_t *opt, double *var)

• 所在文件：pntpos.c
• 功能說明：
• 參數說明：

函數參數，6個
obsd_t   *obs      I   observation data
nav_t    *nav      I   navigation data
double   *azel     I   對於當前定位值，每一顆觀測衛星的 {方位角、高度角}
int      iter      I   迭代次數
prcopt_t *opt      I   processing options
double   *vare     O   僞距測量的碼偏移偏差
返回類型：
double             O   最終能參與定位解算的僞距值

• 調用關係：

st=>start: prange
satsys_=>operation: satsys
testsnr_=>operation: testsnr
gettgd_=>operation: gettgd
e=>end: end

st->satsys_->testsnr_->gettgd_->e

prange satsys testsnr gettgd

• 處理過程：

1. 首先調用 satsys肯定該衛星屬於 RTKLIB設計時給定的幾個導航系統之中。
2. 若是 NFREQ>=3且該衛星屬於 GAL/SBS系統，則 j=2。而若是出現 NFREQ<2||lam[i]==0.0||lam[j]==0.0中的其中一個，直接返回 0.
3. 當 iter>0時，調用 testsnr函數，測試第i、j(IFLC)個頻率信號的載噪比是否符合要求。
4. 計算出 $\gamma$值(f1^2/f22，見ICD-GPS-200C P90)，從 obs和 nav數據中讀出測量僞距值和 碼偏移值(?)。
5. 從數據中讀出的P1_P2==0，則使用 TGD代替，TGD值由 gettgd函數計算獲得。
6. 若是 ionoopt==IONOOPT_IFLC，根據 obs->code的值來決定是否對 P一、P2進行修正，以後再組合出 IFLC時的僞距值(ICD-GPS-200C P91)。不然，則是針對單頻接收即進行的數據處理。先對 P1進行修正，而後再計算出僞距值（PC）
7. 若是 sateph==EPHOPT_SBAS，則還要對 PC進行處理。以後給該函數計算出的僞距值的方差賦值。

• 個人疑惑：

1. 該函數到底在對僞距進行哪部分的計算？計算進行 C/A碼修正後的僞距值？
2. 在調用 testsnr函數時，爲何要有 iter>0的限制？爲何第一次迭代就不能調用這些函數呢？
3. 2中操做的含義不明白，還有爲何出現 3個條件中的一個，就要返回 0呢？
4. 5中關於 IFLC模型僞距的從新計算是看明白了，可是 P1_P2、P1_C1、P1_C2這些變量具體表明什麼含義，以及 P1_P2==0.0時使用 TGD代替和最後關於 sbas clock based C1的操做看不懂。。。

---

satexclude

int satexclude(int sat, int svh, const prcopt_t *opt)

• 所在文件：rtkcmn.c
• 功能說明：檢測某顆衛星在定位時是否須要將其排除，排除標準爲該衛星是否處理選項預先規定的導航系統或排除標誌。
• 參數說明：

函數參數，3個
int       sat     I   satellite number，從 1開始
int       svh     I   sv health flag(0:ok)
prcopt_t  *opt    I   processing options (NULL: not used)
返回類型：
int               O   1:excluded,0:not excluded

• 處理過程：

1. 首先調用 satsys函數獲得該衛星所屬的導航系統。
2. 接着檢驗 svh<0。是，則說明在 ephpos函數中調用 seleph爲該衛星選擇星曆時，並沒有合適的星曆可用，返回 1；否，則進入下一步。
3. 若是處理選項不爲空，則檢測處理選項中對該衛星的排除標誌的值(1:excluded,2:included)，以後再檢測該衛星所屬的導航系統與處理選項中預先設定的是否一致。否，返回 1；是，進入下一步。
4. 若是此時 svh>0，說明此時衛星健康情況出現問題，此衛星不可用，返回 1。

• 注意事項：

1. 3中再在比較該衛星與預先規定的導航系統是否一致時，使用了 sys&opt->navsys來進行比較。這樣作的好處是當 opt->navsys=sys或 opt->navsys=SYS_ALL時，結果都會爲真。之因此會這樣，是由於在 rtklib.h文件中定義這些導航系統變量的時候，所賦的值在二進制形式下都是隻有一位爲 1的數。

• 個人疑惑：

1. 對於 3中檢測，先驗證狀態排除標誌，後驗證導航系統，這樣就可能出現排除標誌符合要求而所屬系統不符合要求的情況，而 3中作法會將上述情況設爲 included!
2. 另外，注意到在 3中檢測以後仍驗證了 svh>0，那若是出現 svh不合乎要求而排除標誌符合要求的情況，3中作法卻會將上述情況設爲 included!

---

ionocorr

int ionocorr(gtime_t time, const nav_t *nav, int sat, const double *pos,
             const double *azel, int ionoopt, double *ion, double *var)

• 所在文件：pntpos.c
• 功能說明：計算給定電離層選項時的電離層延時(m)。
• 參數說明：

函數參數，8個
gtime_t  time      I   time
nav_t    *nav      I   navigation data
int      sat       I   satellite number
double   *pos      I   receiver position {lat,lon,h} (rad|m)
double   *azel     I   azimuth/elevation angle {az,el} (rad)
int      ionoopt   I   ionospheric correction option (IONOOPT_???)
double   *ion      O   ionospheric delay (L1) (m)
double   *var      O   ionospheric delay (L1) variance (m^2)
返回類型：
int                O   (1:ok,0:error)

• 調用關係：

st=>start: ionocorr
c1_=>condition: IONOOPT_BRDC?
ionmodel_=>operation: ionmodel
c2_=>condition: IONOOPT_TEC?
iontec_=>operation: iontec
e=>end: end

st->c1_(yes)->ionmodel_->e
c1_(no)->c2_(yes)->iontec_->e

ionocorr ionmodel iontec

• 處理過程：

1. 根據 opt的值，選用不一樣的電離層模型計算方法。當 ionoopt==IONOOPT_BRDC時，調用 ionmodel，計算 Klobuchar模型時的電離層延時 (L1，m)；當 ionoopt==IONOOPT_TEC時，調用 iontec，計算 TEC網格模型時的電離層延時 (L1，m)。

• 注意事項：

1. 當 ionoopt==IONOOPT_IFLC時，此時經過此函數計算獲得的延時和方差都爲 0。其實，對於 IFLC模型，其延時值在 prange函數中計算僞距時已經包括在裏面了，而方差是在 varerr函數中計算的，而且會做爲導航系統偏差的一部分給出。

---

tropcorr

int int tropcorr(gtime_t time, const nav_t *nav, const double *pos,
                 const double *azel, int tropopt, double *trp, double *var)

• 所在文件：pntpos.c
• 功能說明：計算對流層延時(m)。
• 參數說明：

函數參數，7個
gtime_t  time      I   time
nav_t    *nav      I   navigation data
double   *pos      I   receiver position {lat,lon,h} (rad|m)
double   *azel     I   azimuth/elevation angle {az,el} (rad)
int      tropopt   I   tropospheric correction option (TROPOPT_???)
double   *trp      O   tropospheric delay (m)
double   *var      O   tropospheric delay variance (m^2)
返回類型：
int                O   (1:ok,0:error)

• 調用關係：

st=>start: tropcorr
tropmodel_=>operation: tropmodel
e=>end: end

st->tropmodel_->e

dops tropmodel

• 處理過程：

1. 當 tropopt==TROPOPT_SAAS或一些其它狀況時，調用 tropmodel函數，計算 Saastamoinen模型下的對流層延時。

• 注意事項：

1. 貌似對流層延時與信號頻率無關，因此這裏計算獲得的值並非只針對於 L1信號！

---

varerr

double varerr(const prcopt_t *opt, double el, int sys)

• 所在文件：pntpos.c
• 功能說明：計算導航系統僞距測量值的偏差
• 參數說明：

函數參數，3個
prcopt_t  *opt    I   processing options
double    el      I   elevation angle (rad)
int       sys     I   所屬的導航系統
返回類型：
double            O   導航系統僞距測量值的偏差

• 處理過程：

1. 肯定 sys系統的偏差因子。
2. 計算由導航系統自己所帶來的偏差的方差。
3. 若是 ionoopt==IONOOPT_IFLC時，IFLC模型的方差也會添加到最終計算獲得的方差中。

• 個人疑惑：

1. 本函數總體究竟是爲了計算哪一部分的偏差，仍是沒搞清楚。
2. IFLC模型的方差爲何能夠用 varr*=SQR(3.0)計算？

---

testsnr

int testsnr(int base, int freq, double el, double snr,
            const snrmask_t *mask)

• 所在文件：rtkcmn.c
• 功能說明：檢測接收機所屬類型和頻率信號的有效性
• 參數說明：

函數參數，5個
int    base      I   rover or base-station (0:rover,1:base station)
int    freq      I   frequency (0:L1,1:L2,2:L3,...)
double el        I   elevation angle (rad)
double snr       I   C/N0 (dBHz)
snrmask_t *mask  I   SNR mask
返回類型：
int               O   (1:masked,0:unmasked)

• 調用關係：

st=>start: dops
matmul_=>operation: matmul
matinv_=>operation: matinv
e=>end: end

st->matmul_->matinv_->e

dops matmul matinv

• 處理過程：

1. 知足下列狀況之一 !mask->ena[base]||freq<0||freq>=NFREQ，返回 0.
2. 對 el處理變換，根據後面 i的值獲得不一樣的閾值 minsnr，而對 1<=i<=8的狀況，則使用線性插值計算出 minsnr的值。

• 個人疑惑：

1. 這個函數貌似是根據接收機高度角和信號頻率來檢測該信號是否可用，但 mask在這裏應該翻譯成什麼？看了下調用該函數的地方，返回 0(unmasked)彷佛是合理的、但願看到的狀況，即 snr>=minsnr。
2. 知足 1中條件的狀況，感受應該是不合理的情形，爲何反而返回 0呢？

---

gettgd

double gettgd(int sat, const nav_t *nav)

• 所在文件：pntpos.c
• 功能說明：檢測某顆衛星在定位時是否須要將其排除，排除標準爲該衛星是否處理選項預先規定的導航系統或排除標誌。
• 參數說明：

函數參數，2個
int       sat     I   satellite number，從 1開始
nav_t    *nav     I   navigation data
返回類型：
double            O   tgd parameter (m)

• 處理過程：

1. 從導航數據的星曆中選擇衛星號與 sat相同的那個星曆，讀取 tgd[0]參數後乘上光速。

---

ionmodel

double ionmodel(gtime_t t, const double *ion, const double *pos,
                const double *azel)

• 所在文件：rtkcmn.c
• 功能說明：計算採用 Klobuchar模型時的電離層延時 (L1，m)。
• 參數說明：

函數參數，4個
gtime_t t        I   time (gpst)
double *ion      I   iono model parameters {a0,a1,a2,a3,b0,b1,b2,b3}
double *pos      I   receiver position {lat,lon,h} (rad,m)
double *azel     I   azimuth/elevation angle {az,el} (rad)
返回類型：
double           O   ionospheric delay (L1) (m)

• 處理過程：

1. 主要都是數學計算，其過程能夠在 ICD-GPS-200C P148中找到。

• 注意事項：

1. 這裏計算的電離層延時是相對於 GPS-L1信號而言的，其它頻率信號須要進行一次轉換。
2. 計算過程當中不少角度的單位是半圓，即 $\pi$弧度。在閱讀代碼時，記住這一點很是重要！好比，雖然上述過程與 ICD-GPS-200C P148中一致，但可能與大部分資料上的過程仍是會有所區別。尤爲是下面這個公式。
   $ \Psi = \frac{0.0137}{E+0.01}-0.022 \qquad (ICD-GPS-200C) $
   $ EA = (\frac{445°}{el+20°})-4° \qquad (王虎，GPS精密單點定位中電離層延遲改正模型的研究與分析)$
   可是將下面公式的角度轉化成半圓，即左右兩邊都除以 180，就能夠獲得上面的公式了！
   
   
   

---

iontec

int iontec(gtime_t time, const nav_t *nav, const double *pos,
           const double *azel, int opt, double *delay, double *var)

• 所在文件：ionex.c
• 功能說明：由所屬時間段兩端端點的 TEC網格數據插值計算出電離層延時 (L1) (m)
• 參數說明：

函數參數，3個
gtime_t time     I   time (gpst)
nav_t  *nav      I   navigation data
double *pos      I   receiver position {lat,lon,h} (rad,m)
double *azel     I   azimuth/elevation angle {az,el} (rad)
int    opt       I   model option
                        bit0: 0:earth-fixed,1:sun-fixed
                        bit1: 0:single-layer,1:modified single-layer
double *delay    O   ionospheric delay (L1) (m)
double *var      O   ionospheric dealy (L1) variance (m^2)
返回類型：
int              O   (1:ok,0:error)

• 調用關係：

st=>start: iontec
iondelay_=>operation: iondelay
e=>end: end

st->iondelay_->e

iontec iondelay

• 處理過程：

1. 檢測高度角和接收機高度是否大於閾值。否，則延遲爲 0，方差爲 VAR_NOTEC，返回 1；是，則進入下一步。
2. 從 nav_t.tec中找出第一個 tec[i].time > time（輸入參數，信號接收時間）的 tec數據。而後經過 i==0||i>=nav->nt，確保 time是在所給出的 nav_t.tec包含的時間段之中！
3. 經過確認所找到的時間段的右端點減去左端點，來確保時間間隔 != 0。
4. 調用 iondelay來計算所屬時間段兩端端點的電離層延時。
5. 由兩端的延時，插值計算出觀測時間點處的值。而對於兩端延時的組合，有 3種狀況。
   ① 兩個端點都計算出錯，輸出錯誤信息，返回 0.
   ② 兩個端點都有值，線性插值出觀測時間點的值，返回 1.
   ③ 只有一個端點有值，將其結果做爲觀測時間處的值，返回 1.
   
   
   

• 注意事項：

1. 因爲是經過調用 iondelay來計算所屬時間段端點的電離層延時，因此這裏求出來的值是以 L1信號爲前提的。
2. 關於 5中的第 ②種狀況，RTKLIB manual P152 (E.5.20)式是錯誤的，左端點TEC值得時間權重值應該是 $(t_{i+1}-t)/(t_{i+1}-t_i)$。manual中多是搞反了，源碼中是正確的，與個人見解相同。

• 個人疑惑：

1. 1中當高度角和接收機高度較小時，爲何延遲要爲 0呢？
2. 多是個對最終結果沒有什麼影響的小細節， 雖然時間間隔 tt後面用獲得，可是因爲 2中的操做，其實3中的時間間隔確定是 >0的！
3. 目前關於 tec model，我尚未找到很好的相關方面的文章！

---

iondelay

int iondelay(gtime_t time, const tec_t *tec, const double *pos,
             const double *azel, int opt, double *delay, double *var)

• 所在文件：ionex.c
• 功能說明：根據當前電離層網格模型，計算電離層延時 (L1) (m)。
• 參數說明：

函數參數，3個
gtime_t time      I   time (gpst)
tec_t   *tec      I   tec grid data
double  *pos      I   receiver position {lat,lon,h} (rad,m)
double  *azel     I   azimuth/elevation angle {az,el} (rad)
int     opt       I   model option
                        bit0: 0:earth-fixed,1:sun-fixed
                        bit1: 0:single-layer,1:modified single-layer
double *delay    O    ionospheric delay (L1) (m)
double *var      O    ionospheric dealy (L1) variance (m^2)
返回類型：
int              O    (1:ok,0:error)

• 調用關係：

st=>start: iondelay
ionppp_=>operation: ionppp
interptec_=>operation: interptec
e=>end: end

st->ionppp_->interptec_->e

iondelay ionppp interptec

• 處理過程：

1. 按照 RTKLIB manual P152中的公式(E.5.19)，先計算出與頻率有關的中間項。
2. 總體過程是按照電離層的高度，從起始高度開始，逐層計算每一層的延時和方差，以後累加到一塊兒。下面再具體闡述。
3. 首先調用 ionppp函數，計算出在當前電離層高度時，電離層穿刺點的位置 {lat,lon,h} (rad,m)和傾斜率( $1/cos \zeta$)。
4. 按照 opt的值可能再次進行修正。opt&1，則按照 M-SLM映射函數從新計算傾斜率；opt&2，則在日固座標系中考慮地球自轉，從新計算穿刺點經度；
5. 由 TEC網格計算穿刺點的電子數總量，而後按照 (E.5.19)累加電離層延時(m)和方差。

• 注意事項：

1. 這裏在計算電離層延時(m)時，是假設信號爲 L1的！

• 個人疑惑：

1. 4中操做看不懂，貌似跟 地/日座標系和 SL/M-SL有關？
2. IONEX文件中的電離層模型，高度真的有不少層嗎？

---

ionppp

double ionppp(const double *pos, const double *azel, double re,
              double hion, double *posp)

• 所在文件：rtkcmn.c
• 功能說明：計算電離層穿刺點的位置 {lat,lon,h} (rad,m)和傾斜率( $1/cos \zeta$)。
• 參數說明：

函數參數，5個
double *pos      I   receiver position {lat,lon,h} (rad,m)
double *azel     I   azimuth/elevation angle {az,el} (rad)
double re        I   earth radius (km)
double hion      I   altitude of ionosphere (km)
double *posp     O   pierce point position {lat,lon,h} (rad,m)
返回類型：
double           O   傾斜率

• 處理過程：

1. 與處理過程相對應的公式，請見 RTKLIB manual P151

• 注意事項：

1. 說明文檔中的 z並非仰角azel[1]，而是仰角關於$\pi/2$的補角，因此程序中在計算 rp是採用的是 cos(azel[1])的寫法。
2. 可能由於後面再從 TEC網格數據中插值時，並不須要高度信息，因此這裏穿刺點位置 posp中的第三項高度，其實並無進行賦值，

---

interptec

int interptec(const tec_t *tec, int k, const double *posp, double *value,
              double *rms)

• 所在文件：ionex.c
• 功能說明：經過在經緯度網格點上進行雙線性插值，計算第 k個高度層時穿刺點處的電子數總量 TEC
• 參數說明：

函數參數，5個
tec_t   *tec       I   tec grid data
int     k          I   高度方向上的序號，能夠理解爲電離層序號
double  *posp      I   pierce point position {lat,lon,h} (rad,m)
double  *value     O   計算獲得的刺穿點處的電子數總量(tecu) 
double  *rms       O   所計算的電子數總量的偏差的標準差(tecu) 
返回類型：
int                O    (1:ok,0:error)

• 處理過程：

1. 將 value和 rms所指向的值置爲 0。
2. 檢驗 tec的緯度和經度間隔是否爲 0。是，則直接返回 0。
3. 將穿刺點的經緯度分別減去網格點的起始經緯度，再除以網格點間距，對結果進行取整，獲得穿刺點所在網格的序號和穿刺點所在網格的位置(比例)。
4. 按照下圖的順序，調用 dataindex函數分別計算這些網格點的 tec數據在 tec.data中的下標，從而獲得這些網格點處的 TEC值和相應偏差的標準差。
   
   
5. 若是四個網格點的 TEC值都 >0，則說明穿刺點位於網格內，使用雙線性插值計算出穿刺點的 TEC值；不然使用最鄰近的網格點值做爲穿刺點的 TEC值，不過前提是網格點的 TEC>0；不然，選用四個網格點中 >0的值的平均值做爲穿刺點的 TEC值。

• 注意事項：

1. 對於lats[3]，其含義分別爲起始緯度、終止緯度和間隔，對 lons[3]、hgts[3]，其含義也是相似的。
2. 對於 dataindex函數， i、j、k都是從 0開始的，意味着分別表明各自方向上第 i+一、j+一、k+1層，而且是按照緯度、經度、高度的優先順序來存儲網格點數據的。

• 個人疑惑：

1. 關於輸出參數 rms，按照其名稱，應該是 均方根值，可是在調用了該函數的 iondelay中，確是把它的平方當作方差的一部分進行累加。因此我估計 tec.rms值得應該是相應網格點數據值的方差？
2. 老實說，源碼中關於 tec->lons[2]大於或者小於 0所作的處理，並無看得太明白。另外，我的感受 3中減去網格點起始經度後的差值應該也不會超過 360°吧？
3. 總體由網格點數據插值穿刺點值得過程能夠明白，但 tec.data會 <0嗎？還有在網格外是指某一個網格的外面，仍是總體 TEC大網格的四個角外面？

---

tropmodel

double tropmodel(gtime_t time, const double *pos, const double *azel,
                 double humi)

• 所在文件：rtkcmn.c
• 功能說明：由標準大氣和 Saastamoinen模型，計算電離層延時(m)
• 參數說明：

函數參數，4個
gtime_t time     I   time
double *pos      I   receiver position {lat,lon,h} (rad,m)
double *azel     I   azimuth/elevation angle {az,el} (rad)
double humi      I   relative humidity
返回類型：
double           O   tropospheric delay (m)

• 處理過程：

1. 與處理過程相對應的公式，請見 RTKLIB manual P149

• 個人疑惑：

1. 源碼中關於 trph的計算，與大多數文獻和 RTKLIB manual P149 (E.5.4)有所不一樣，咋回事兒呢？

---

resdop

int resdop(const obsd_t *obs, int n, const double *rs, const double *dts,
           const nav_t *nav, const double *rr, const double *x,
           const double *azel, const int *vsat, double *v, double *H)

• 所在文件：pntpos.c
• 功能說明：計算定速方程組左邊的幾何矩陣和右端的速度殘餘，返回定速時所使用的衛星數目
• 參數說明：

函數參數，11個：
obsd_t   *obs      I   observation data
int      n         I   number of observation data
double   *rs       I   satellite positions and velocities，長度爲6*n，{x,y,z,vx,vy,vz}(ecef)(m,m/s)
double   *dts      I   satellite clocks，長度爲2*n， {bias,drift} (s|s/s)
nav_t    *nav      I   navigation data
double   *rr       I   receiver positions and velocities，長度爲6，{x,y,z,vx,vy,vz}(ecef)(m,m/s)
double   *x        I   本次迭代開始以前的定速值
double   *azel     I   azimuth/elevation angle (rad)
int      *vsat     I   表徵衛星在定速時是否有效
double   *v        O   定速方程的右端部分，速度殘餘
double   *H        O   定速方程中的幾何矩陣
返回類型:
int                O    定速時所使用的衛星數目

• 調用關係：

st=>start: resdop
ecef2pos_=>operation: ecef2pos
xyz2enu_=>operation: xyz2enu
matmul_=>operation: matmul
e=>end: end

st->ecef2pos_->xyz2enu_->matmul_->e

resdop ecef2pos xyz2enu matmul

• 處理過程：

1. 調用 ecef2pos函數，將接收機位置由 ECEF轉換爲大地座標系。
2. 調用 xyz2enu函數，計算此時的座標轉換矩陣。
3. 知足下列條件 obs[i].D[0]==0.0||lam==0.0||!vsat[i]||norm(rs+3+i*6,3)<=0.0之一，則當前衛星在定速時不可用，直接進行下一次循環。
4. 計算當前接收機位置下 ENU中的視向量，而後轉換獲得 ECEF中視向量的值。
5. 計算 ECEF中衛星相對於接收機的速度，而後再計算出考慮了地球自轉的用戶和衛星之間的幾何距離變化率，校訂公式見 RTKLIB manual P159 (E.6.29)。
6. 根據公式 計算出方程右端項的多普勒殘餘，而後再構建左端項的幾何矩陣。最後再將觀測方程數增 1.

• 注意事項：

1. 這裏與定位不一樣，構建幾何矩陣時，就只有 4個未知數，而定位時是有 NX個。而且沒有像定位那樣爲了防止虧秩而進行約束處理。
2. 多普勒定速方程中幾何矩陣 G與定位方程中的同樣，前三行都是 ECEF座標系中由接收機指向衛星的單位觀測矢量的反向。而因爲轉換矩陣 S自己是一個正交單位矩陣($S^T=S{-1}$)，因此這裏在計算 ECEF中的視向量時，對 E進行了轉置處理。
3. 這裏構建的左端幾何矩陣 H，也與定位方程中的同樣，是大部分資料上的幾何矩陣的轉置。