RTKLIB源码解析（一）——单点定位(pntpos.c)
RTKLIB源碼解析（一）——單點定位(pntpos.c)
標簽： GNSS RTKLIB 單點定位
前段時間一直忙著寫畢業論文，所以也沒有太多時間來閱讀 RTKLIB源碼，最近好歹把 pntpos中的相關代碼看了一遍，知道了 RTKLIB是如何實現單點偽距定位的。這里把每一個函數都做成了小卡片的形式，每個函數大都包含函數簽名、所在文件、功能說明、參數說明、處理過程、注意事項和我的疑惑這幾個部分，介紹了閱讀代碼時我自己的看法和疑惑。所以希望諸位看官能幫忙解答我的疑惑，與我交流，也希望能幫助后來也有需要閱讀 RTKLIB源碼的人，給他們多提供一份思路。總而言之，既為人，也為己。
這份文檔是使用 Cmd Markdown完成的，在作業部落上其格式顯式的非常完整，但是在博客園中目錄、代碼塊和流程圖似乎都沒有顯示出來，所以這里也貼上本文在作業部落上的鏈接 RTKLIB源碼解析（一）——單點定位(pntpos.c)，對格式“零容忍”的同學請移步那里。
目錄RTKLIB源碼解析（一）——單點定位(pntpos.c)pntpossatpossestposraim_fdeestvelephclksatpossatsysselepheph2clkephposeph2posrescodelsqvalsolmatmuldopsecef2enuxyz2enuecef2posgeodistsatazelprangesatexcludeionocorrtropcorrvarerrtestsnrgettgdionmodelionteciondelayionpppinterptectropmodelresdop
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
處理過程：
檢查衛星個數是否大于 0
當處理選項opt中的模式不是單點模式時，電離層校正采用Klobuchar模型，對流層校正則采用Saastamoinen模型；相反，當其為單點模式時，對輸入參數opt不做修改。
調用 satposs函數，按照所觀測到的衛星順序計算出每顆衛星的位置、速度、{鐘差、頻漂}。
調用 estpos函數，通過偽距實現絕對定位，計算出接收機的位置和鐘差，順帶返回實現定位后每顆衛星的{方位角、仰角}、定位時有效性、定位后偽距殘差。
調用 raim_fde函數，對上一步得到的定位結果進行接收機自主正直性檢測（RAIM）。通過再次使用vsat數組，這里只會在對定位結果有貢獻的衛星數據進行檢測。
調用 estvel函數，依靠多普勒頻移測量值計算接收機的速度。這里只使用通過了上一步RAIM_FDE操作的衛星數據，所以對于計算出的速度就沒有再次進行RAIM了。
首先將ssat_t結構體數組的 vs(定位時有效性)、azel（方位角、仰角）、resp(偽距殘余)、resc(載波相位殘余)和 snr(信號強度)都置為 0，然后再將實現定位后的 azel、snr賦予ssat_t結構體數組，而 vs、resp則只賦值給那些對定位有貢獻的衛星，沒有參與定位的衛星，這兩個屬性值為 0。
注意事項：
這里只計算了接收機的鐘差，而沒有計算接收機的頻漂，原因在于 estvel函數中雖然計算得到了接收機頻漂，但并沒有將其輸出到
sol_t:dtr中。
C語言中用malloc申請的內存需要自己調用free來予以回收，源碼中的mat、imat、zeros等函數都只是申請了內存，并沒有進行內存的回收，在使用這些函數時，用戶必須自己調用free來回收內存！源碼中將使用這些函數的代碼放置在同一行，在調用函數結尾處也統一進行內存回收，位置較為明顯，不致于輕易忘記。
我的疑惑：
源碼中將
obs[0].time作為星歷選擇時間傳遞給 satposs函數，這樣對于每一顆觀測衛星，都要使用第一顆觀測衛星的數據接收時間作為選擇星歷的時間標準。是否應該每顆衛星都使用自己的觀測時間？或者應該使用每顆衛星自己的信號發射時間？，還是說這點差別對選擇合適的星歷其實沒有關系？
這里規定能夠執行 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
處理過程：
按照觀測數據的順序，首先將將當前觀測衛星的 rs、dts、var和svh數組的元素置 0。
通過判斷某一頻率下信號的偽距是否為 0，來得到此時所用的頻率個數。注意，頻率個數不能大于NFREQ（默認為 3）。
用數據接收時間減去偽距信號傳播時間，得到衛星信號的發射時間。
調用 ephclk函數，由廣播星歷計算出當前觀測衛星的鐘差。注意，此時的鐘差是沒有考慮相對論效應和 TGD的。
用 3中的信號發射時間減去 4中的鐘偏，得到 GPS時間下的衛星信號發射時間。
調用 satpos函數，計算信號發射時刻衛星的 P(ecef,m)、V(ecef,m/s)、C((s|s/s))。注意，這里計算出的鐘差是考慮了相對論效應的了，只是還沒有考慮 TGD。
如果由 6中計算出的鐘偏為 0，就再次調用 ephclk函數，將其計算出的衛星鐘偏作為最終的結果。
注意事項：
ephclk函數計算的鐘偏是沒有考慮相對論和 TGD的，而 satpos函數考慮了相對論，沒有考慮 TGD。
我的疑惑：
對于處理過程中的第3步，數據接收時間減去偽距信號傳播時間，這里的數據接收時間應該是有接收機得到的，本身應該是包含接收機鐘差的，所以最終得到的信號發射時間應該也是不準確的。難道說接收機鐘差較小，在此時可以忽略不計？
ephclk函數最終是通過調用 eph2clk來計算衛星鐘偏的，但對于后者，我有問題。見 eph2clk處我的疑惑
為什么要進行 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
處理過程：
將
sol->rr的前 3項賦值給 x數組。
開始迭代定位計算，首先調用 rescode函數，計算在當前接收機位置和鐘差值的情況下，定位方程的右端部分v(nv*1)、幾何矩陣H(NX*nv)、此時所得的偽距殘余的方差var、所有觀測衛星的azel{方位角、仰角}、定位時有效性vsat、定位后偽距殘差resp、參與定位的衛星個數ns和方程個數nv。
確定方程組中方程的個數要大于未知數的個數。
以偽距殘余的標準差的倒數作為權重，對 H和 v分別左乘權重對角陣，得到加權之后的 H和 v。
調用 lsq函數，根據 $Delta x=(HHT){-1}Hv$和 $Q=(HHT){-1}$，得到當前 x的修改量和定位誤差協方差矩陣中的權系數陣。
將 5中求得的 x加入到當前 x值中，得到更新之后的 x值。
如果 5中求得的修改量小于截斷因子(目前是1e-4)，則將 6中得到的 x值作為最終的定位結果，對sol的相應參數賦值，之后再調用 valsol函數確認當前解是否符合要求（偽距殘余小于某個 $chi^2$值和GDOP小于某個門限值）。否則，進行下一次循環。
如果超過了規定的循環次數，則輸出發散信息后，返回 0。
注意事項：
關于第 1步，如果是第一次定位，即輸入的
sol為空，則 x初值為 0；如果之前有過定位，則通過 1中操作可以將上一歷元的定位值作為該歷元定位的初始值。
關于定位方程，書中的寫法如下：
$$
G egin{bmatrix} Delta x Delta y Delta z delta t_u end{bmatrix} = b
$$
其中，
$$
G = egin{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 = egin{bmatrix}
ho_r1-(r_r1+c·dt_r-c·dt_s1+I1+T^1)
ho_r2-(r_r2+c·dt_r-c·dt_s2+I2+T^2)
cdots
ho_rN-(r_rN+c·dt_r-c·dt_sN+IN+T^N)
end{bmatrix}
$$
而 RTKLIB中采用的方程是下面這樣的（這里先假設未知數的個數為 4）：
$$
H^T egin{bmatrix} Delta x Delta y Delta z delta t_u end{bmatrix} = b
$$
其中，$H=G^T$。
關于加權最小二乘，這里的權重值是對角陣，這是建立在假設不同測量值的誤差之間是彼此獨立的；另外，這個權重值并不單是距測量誤差的，而是方程右端b整體的測量誤差。最后，大部分資料上這里都是把權重矩陣W保留到方程的解的表達式當中，而這里是直接對 H和 v分別左乘權重對角陣，得到加權之后的 H和 v，其表示形式像是沒有加權一樣。
如果某次迭代過程中步長小于門限值(1e-4)，但經 valsol函數檢驗后該解無效，則會直接返回 0，并不會再進行下一次迭代計算。
因為該函數有兩個返回路徑，而且又調用了mat函數來構建矩陣，所以在兩個返回路徑處都需要調用free函數來回收內存。
源碼中的 dtr實際上單位是 m，是乘以了光速之后的。
在對sol結構體賦值時，并沒有直接將接收機鐘差 dtr賦值給sol_dtr，而是通過在sol中存儲的是減去接收機鐘差后的信號觀測時間，來講該信息包括到sol中了。
源碼中定位方程的個數nv要大于有效觀測衛星的個數ns，這里為了防止虧秩，似乎又加了 3個未知數和觀測方程。
在每一次重新調用 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
處理過程：
關于觀測衛星數目的循環，每次舍棄一顆衛星，計算使用余下衛星進行定位的定位值。
舍棄一顆衛星后，將剩下衛星的數據復制到一起，調用 estpos函數，計算使用余下衛星進行定位的定位值。
累加使用當前衛星實現定位后的偽距殘差平方和與可用微信數目，如果nvsat<5，則說明當前衛星數目過少，無法進行RAIM_FDE操作。
計算偽距殘差平方和的標準平均值，如果大于rms，則說明當前定位結果更合理，將stat置為 1，重新更新sol、azel、vsat(當前被舍棄的衛星，此值置為0)、resp等值，并將當前的rms_e更新到 `rms'中。
繼續棄用下一顆衛星，重復 2-4操作。總而言之，將同樣是棄用一顆衛星條件下，偽距殘差標準平均值最小的組合所得的結果作為最終的結果輸出。
如果 stat不為 0，則說明在棄用衛星的前提下有更好的解出現，輸出信息，指出棄用了哪科衛星。
使用free函數回收相應內存。
注意事項：
使用了
mat和malloc函數后要使用free函數來回收內存。
源碼中有很多關于 i、j、k的循環。其中，i表示最外面的大循環，每次將將第 i顆衛星舍棄不用，這是通過if (j==i) continue實現的；j表示剩余使用的衛星的循環，每次進行相應數據的賦值；k表示參與定位的衛星的循環，與 j一起使用。
我的疑惑：
這里執行
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
處理過程：
在最大迭代次數限制內，調用 resdop，計算定速方程組左邊的幾何矩陣和右端的速度殘余，返回定速時所使用的衛星數目。
調用 lsq函數，解出 {速度、頻漂}的步長，累加到 x中。
檢查當前計算出的步長的絕對值是否小于1E-6。是，則說明當前解已經很接近真實值了，將接收機三個方向上的速度存入到sol_t.rr中；否，則進行下一次循環。
執行完所有迭代次數，使用free回收內存。
注意事項：
最終向
sol_t類型存儲定速解時，并沒有存儲所計算出的接收器時鐘頻漂。
這里不像定位時，初始值可能為上一歷元的位置(從 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
處理過程：
首先調用 satsys函數，根據衛星編號確定該衛星所屬的導航系統和該衛星在該系統中的
PRN編號。
對于 GPS導航系統，調用 seleph函數來選擇toe值與星歷選擇時間標準teph最近的那個星歷。
調用 eph2clk函數，通過廣播星歷和信號發射時間計算出衛星鐘差。
注意事項：
此時計算出的衛星鐘偏是沒有考慮相對論效應和 TGD的。
目前我還只關心 RTKLIB對于 GPS導航所進行的數據處理，所以這里在選擇星歷和計算鐘差時都只介紹與 GPS系統有關的函數。
我的疑惑：
見 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
處理過程：
判斷星歷選項的值，如果是
EPHOPT_BRDC，調用 ephpos函數，根據廣播星歷計算出算信號發射時刻衛星的 P、V、C
注意事項：
此時計算出的衛星鐘差考慮了相對論，還沒有考慮 TGD。
目前還只閱讀了如何從廣播星歷中計算衛星 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-MAXSAT之內），先令衛星系統為SYS_NONE。
按照 TRKLIB中定義相應宏的順序來判斷是否是 GPS、GLO、GAL系統等，判斷標準是將衛星號減去前面的導航系統所擁有的衛星個數，來判斷剩余衛星個數是否小于等于本系統的衛星個數。
確定所屬的系統后，通過加上最小衛星編號的PRN再減去 1，得到該衛星在該系統中的PRN編號。
注意事項：
這里的衛星號是從 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 * 星歷數據
處理過程：
根據該衛星所屬的導航系統給與星歷參考時間的最大時間間隔
tmax賦予相應的值。
遍歷導航數據，首先確保所查找星歷的衛星號是否相同，接著確保星歷期號是否相同。
接著確保星歷選擇時間與代查找星歷的參考時間的間隔是否小于tmax。
對于通過了 2-3驗證的星歷，如果此時對于輸入的期號，有iode>=0，則該星歷就是所要尋找的星歷；否則，驗證 3中的t是否滿足t<=tmin。滿足的話，就令tmin=t，該星歷目前是距離參考時間最近的。
循環 2-4步操作，遍歷完導航數據中的所有星歷。如果都沒有符合條件的，就輸出信息并返回 NULL；否則，返回所查找到的星歷。
注意事項：
對于 GPS系統，星歷數據期號每 2h更新一次，所以 RTKLIB中對
MAXDTOE的定義為 7200。
如果存在兩個相鄰時間段的星歷，通過 4中令tmin=t的操作可以最終查找出參考時間距time最近的那個星歷。
我的疑惑：
為什么
tmax和tmin都要加 1？
IODE正常情況下應該都是 >=0的，為什么還要考慮 <0的情況？
考慮到星歷中衛星的健康狀況，這里在選擇星歷時是否也應該驗證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
處理過程：
計算與星歷參考時間的偏差 dt = t-toc。
利用二項式校正計算出衛星鐘差，從 dt中減去這部分，然后再進行一次上述操作，得到最終的 dt。
使用二項式校正得到最終的鐘差。
- **我的疑惑**：
> 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
處理過程：
調用 satsys函數，確定該衛星所屬的導航系統。
如果該衛星屬于 GPS系統，則調用 seleph函數來選擇廣播星歷。
根據選中的廣播星歷，調用 eph2pos函數來計算信號發射時刻衛星的 位置、鐘差和相應結果的誤差。
在信號發射時刻的基礎上給定一個微小的時間間隔，再次計算新時刻的 P、V、C。與 3結合，通過擾動法計算出衛星的速度和頻漂。
注意事項：
這里是使用擾動法計算衛星的速度和頻漂，并沒有使用那些位置和鐘差公式對時間求導的結果。
由于是調用的 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
處理過程：
與大部分資料上計算衛星位置和鐘差的過程是一樣的，只是這里在計算偏近點角 E時采用的是牛頓法來進行迭代求解。
計算誤差直接采用URA值來標定，具體對應關系可在ICD-GPS-200C P83中找到。
注意事項：
這里在計算鐘差時，就與大部分資料上介紹的一樣了，只進行了一次二項式校正。另外，這里的鐘差考慮了相對論效應，并沒有考慮 TGD。
在計算偏近點角 E時，這里采用的是牛頓法來進行迭代求解，這里與很多資料上可能不一樣，具體內容可在RTKLIB manual P143中找到。
補充幾個程序中的參數說明：
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
處理過程：
將之前得到的定位解信息賦值給
rr和dtr數組，以進行關于當前解的偽距殘余的相關計算。
調用 ecef2pos函數，將上一步中得到的位置信息由 ECEF轉化為大地坐標系。
將vsat、azel和resp數組置 0，因為在前后兩次定位結果中，每顆衛星的上述信息都會發生變化。
調用 satsys函數，驗證衛星編號是否合理及其所屬的導航系統。
檢測當前觀測衛星是否和下一個相鄰數據重復。是，則i++后繼續下一次循環；否，則進入下一步。
調用 geodist函數，計算衛星和當前接收機位置之間的幾何距離和receiver-to-satellite方向的單位向量。然后檢驗幾何距離是否 >0。
調用 satazel函數，計算在接收機位置處的站心坐標系中衛星的方位角和仰角，檢驗仰角是否 $geq$截斷值。
調用 prange函數，得到偽距值。
可以在處理選項中事先指定只選用哪些導航系統或衛星來進行定位，這是通過調用 satexclude函數完成的。
調用 ionocorr函數，計算電離層延時(m)。
上一步中所得的電離層延時是建立在 L1信號上的，當使用其它頻率信號時，依據所用信號頻組中第一個頻率的波長與 L1波長的關系，對上一步得到的電離層延時進行修正。
調用 tropcorr函數，計算對流層延時(m)。
由 $ho_ri-(r_ri+dt_r-c·dt_si+Ii+T^i)$，計算出此時的偽距殘余。
組裝幾何矩陣，前 3行為 6中計算得到的視線單位向量的反向，第 4行為 1，其它行為 0。
將參與定位的衛星的定位有效性標志設為 1，給當前衛星的偽距殘余賦值，參與定位的衛星個數ns加 1.
調用 varerr函數，計算此時的導航系統誤差（可能會包括IFLC選項時的電離層延時），然后累加計算用戶測距誤差(URE)。
為了防止虧秩，人為的添加了幾組觀測方程。
注意事項：
返回值
v和resp的主要區別在于長度不一致，v是需要參與定位方程組的解算的，維度為 nv1；而 resp僅表示所有觀測衛星的偽距殘余，維度為 n1，對于沒有參與定位的衛星，該值為 0。
源碼中dtr的單位是 m。
16中的URE值包括 ①衛星星歷和鐘差的誤差 ②大氣延時誤差 ③偽距測量的碼偏移誤差 ④導航系統的誤差
我的疑惑：
關于 5中的去除重復數據的過程，有以下幾個看法：
① 這樣做的前提是相同衛星的重復數據必須相鄰出現！
② 為什么在這里要進行重復數據檢測，在構建obsd_t結構體時就可以進行這項工作呀？
③ 5中當數據重復時，i++后繼續下一次循環，這樣的話會直接略去后面所重復的數據，這樣做就會將兩個相鄰重復數據都忽略掉，就相當于重復數據都不使用。感覺可以用其中一個的啊！
11步中，為什么要選擇所用信號頻組中第一個頻率的波長來進行電離層延時修正呢？還有，電離層延時的值發生了改變，那這里的方差是否也需要跟著一起改變呢？
在計算電離/對流層延時時，均傳入了iter>0?opt->ionoopt:IONOOPT_BRDC或iter>0?opt->tropopt:TROPOPT_SAAS參數，都強調了當iter==0時，會強制使用Klobuchar或Saastamoinen模型。這會不會是因為這兩種模型都是屬于對接收機位置不是非常敏感的類型？
這里虧秩應該就是欠定方程的意思吧？這里 17中的操作沒有看懂，也沒能找到相關理論依據
lsq
int lsq(const double *A, const double *y, int n, int m, double *x, double *Q)
所在文件：rtkcmn.c
功能說明：計算得到方程 $x=(AAT){-1}Ay$ 左邊的值和該值的協方差矩陣 $Q=(AAT){-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
處理過程：
首先計算右半部分 $A_y=Ay$。
再計算左半部分括號里面的值 $Q=AA^T$。
計算 Q矩陣的逆 $Q^{-1}$，但仍存儲在 Q中，最后再右乘 $A_y$，得到 x的值。
注意事項：
對于加權最小二乘，可以直接在調用該函數之前直接將 A、y進行加權處理，之后在調用該函數，這樣得到的就是加權最小二乘的解。
所有的矩陣都是列優先存儲的，對于整個源代碼來說，矩陣都是這樣存儲的。所以對于代碼中出現的一維矩陣，基本都應該是列向量。在閱讀數組下標時，記住這一點是非常重要的。
矩陣求逆并不簡單，尤其是對于接近奇異的矩陣。但是由于這是個基本功能，并不打算繼續深入下去。
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
處理過程：
先計算定位后偽距殘余平方加權和
vv。
檢查是否滿足 $vv>chi^2(nv-nx-1)$。是，則說明此時的定位解誤差過大，返回 0；否則，轉到下一步。
復制 azel，這里只復制那些對于定位結果有貢獻的衛星的 zael值，并且統計實現定位時所用衛星的數目。
調用 dops函數，計算各種精度因子(DOP)，檢驗是否有0<GDOP<max。否，則說明該定位解的精度不符合要求，返回 0；是，則返回 1。
注意事項：
關于 2中的 $chi^2$檢驗，這里這里與
RTKLIB manual P160中不一致（但與教材一致），文檔中滿足 $frac{vTv}{nv-nx-1}>chi2(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
處理過程：
根據
tr的值確定矩陣相乘標志，即
$$
egin{matrix}
AB ightarrow NN ightarrow 1
AB^T ightarrow NT ightarrow 2
A^T*B ightarrow TN ightarrow 3
AT*BT ightarrow TT ightarrow 4
end{matrix}
$$
按照f的值，分別執行相應的元素相乘并累加操作。
對于 2中得到的乘積 tmp，執行 C = alphatmp + betaC操作。
注意事項：
在調用該函數時，要確保矩陣的維度是否和上述參數說明一致。
所有的矩陣都是列優先存儲的，記住這一點，對于看明白 2中不同相乘方式時，元素如何相乘累加是至關重要的。
矩陣求逆并不簡單，尤其是對于接近奇異的矩陣。但是由于這是個基本功能，并不打算繼續深入下去。
dops
void dops(int ns, const double *azel, double elmin, double *dop)
所在文件：rtkcmn.c
功能說明：由站心坐標系時的權系數矩陣表達式$ ilde{H} = ( ilde{G} ilde{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
處理過程：
先按照如下站心坐標系時的幾何矩陣 $ ilde{G}$的表達式求出其值。
$$
ilde{G} = egin{bmatrix}
cos heta1sinalpha1 & cos heta2sinalpha2 & vdots & cos hetaNsinalphaN
cos heta1cosalpha1 & cos heta2cosalpha2 & vdots & cos hetaNcosalphaN
sin heta^1 & sin heta^2 & vdots & sin heta^N
1 & 1 & vdots & 1
end{bmatrix}
$$
檢驗上述矩陣的列數是否$geq$ 4。
調用 matmul和matinv函數，計算出 $ ilde{H} = ( ilde{G} ilde{G}T){-1}$的值。
如果能正確計算出逆矩陣，就按照順序計算出 GDOP、PDOP、HDOP和 VDOP的值。
我的疑惑：
1中幾何矩陣 $ ilde{G}$與書中的不一致，前 3行均少了一個負號，我自己推導之后也覺得應該有負號，即為
$$
ilde{G} = egin{bmatrix}
-cos heta1sinalpha1 & -cos heta2sinalpha2 & vdots & -cos hetaNsinalphaN
-cos heta1cosalpha1 & -cos heta2cosalpha2 & vdots & -cos hetaNcosalphaN
-sin heta^1 & -sin heta^2 & vdots & -sin heta^N
1 & 1 & vdots & 1
end{bmatrix}
$$
不過，由于 $ ilde{H} = ( ilde{G} ilde{G}T){-1}$，在括號里面的矩陣相乘時，是否有負號只對底邊靠左 3個元素和右邊靠上 3個元素有影響（多了個負號），然后再進行求逆之后，前 3行是否有負號就對對角線上的元素似乎沒有影響了。
感覺在 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
處理過程：
先調用 xyz2enu函數計算此時的坐標變換矩陣 E。
調用 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
處理過程：
按照大部分資料上的寫法計算 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)
處理過程：
檢查衛星到 WGS84坐標系原點的距離是否大于基準橢球體的長半徑。
ps-pr，計算由接收機指向衛星方向的觀測矢量，之后在計算出相應的單位矢量。
考慮到地球自轉，即信號發射時刻衛星的位置與信號接收時刻衛星的位置在 WGS84坐標系中并不一致，進行關于Sagnac效應的校正。
我的疑惑：
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
處理過程：
調用 ecef2enu函數，將
pos處的向量轉換到以改點為原點的站心坐標系中。
調用atan2函數計算出方位角，然后再算出仰角。
注意事項：
這里在計算方位角時，要使用
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
處理過程：
首先調用 satsys確定該衛星屬于 RTKLIB設計時給定的幾個導航系統之中。
如果NFREQ>=3且該衛星屬于GAL/SBS系統，則j=2。而如果出現NFREQ<2||lam[i]==0.0||lam[j]==0.0中的其中一個，直接返回 0.
當iter>0時，調用 testsnr函數，測試第i、j(IFLC)個頻率信號的載噪比是否符合要求。
計算出 $gamma$值(f12/f22，見ICD-GPS-200C P90)，從obs和nav數據中讀出測量偽距值和碼偏移值(?)。
從數據中讀出的P1_P2==0，則使用 TGD代替，TGD值由 gettgd函數計算得到。
如果ionoopt==IONOOPT_IFLC，根據obs->code的值來決定是否對 P1、P2進行修正，之后再組合出 IFLC時的偽距值(ICD-GPS-200C P91)。否則，則是針對單頻接收即進行的數據處理。先對 P1進行修正，然后再計算出偽距值（PC）
如果sateph==EPHOPT_SBAS，則還要對 PC進行處理。之后給該函數計算出的偽距值的方差賦值。
我的疑惑：
該函數到底在對偽距進行哪部分的計算？計算進行 C/A碼修正后的偽距值？
在調用 testsnr函數時，為什么要有iter>0的限制？為什么第一次迭代就不能調用這些函數呢？
2中操作的含義不明白，還有為什么出現 3個條件中的一個，就要返回 0呢？
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
處理過程：
首先調用 satsys函數得到該衛星所屬的導航系統。
接著檢驗svh<0。是，則說明在 ephpos函數中調用 seleph為該衛星選擇星歷時，并無合適的星歷可用，返回 1；否，則進入下一步。
如果處理選項不為空，則檢測處理選項中對該衛星的排除標志的值(1:excluded,2:included)，之后再檢測該衛星所屬的導航系統與處理選項中預先設定的是否一致。否，返回 1；是，進入下一步。
如果此時svh>0，說明此時衛星健康狀況出現問題，此衛星不可用，返回 1。
注意事項：
3中再在比較該衛星與預先規定的導航系統是否一致時，使用了
sys&opt->navsys來進行比較。這樣做的好處是當opt->navsys=sys或opt->navsys=SYS_ALL時，結果都會為真。之所以會這樣，是因為在rtklib.h文件中定義這些導航系統變量的時候，所賦的值在二進制形式下都是只有一位為1的數。
我的疑惑：
對于 3中檢測，先驗證狀態排除標志，后驗證導航系統，這樣就可能出現排除標志符合要求而所屬系統不符合要求的狀況，而 3中做法會將上述狀況設為
included!
另外，注意到在 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
處理過程：
根據
opt的值，選用不同的電離層模型計算方法。當ionoopt==IONOOPT_BRDC時，調用 ionmodel，計算Klobuchar模型時的電離層延時 (L1，m)；當ionoopt==IONOOPT_TEC時，調用 iontec，計算TEC網格模型時的電離層延時 (L1，m)。
注意事項：
當
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
處理過程：
當
tropopt==TROPOPT_SAAS或一些其它情況時，調用 tropmodel函數，計算Saastamoinen模型下的對流層延時。
注意事項：
貌似對流層延時與信號頻率無關，所以這里計算得到的值并不是只針對于
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 導航系統偽距測量值的誤差
處理過程：
確定 sys系統的誤差因子。
計算由導航系統本身所帶來的誤差的方差。
如果ionoopt==IONOOPT_IFLC時，IFLC模型的方差也會添加到最終計算得到的方差中。
我的疑惑：
本函數整體到底是為了計算哪一部分的誤差，還是沒搞清楚。
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
處理過程：
滿足下列情況之一
!mask->ena[base]||freq<0||freq>=NFREQ，返回 0.
對el處理變換，根據后面i的值得到不同的閾值minsnr，而對1<=i<=8的情況，則使用線性插值計算出minsnr的值。
我的疑惑：
這個函數貌似是根據接收機高度角和信號頻率來檢測該信號是否可用，但
mask在這里應該翻譯成什么？看了下調用該函數的地方，返回 0(unmasked)似乎是合理的、希望看到的情況，即snr>=minsnr。
滿足 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)
處理過程：
從導航數據的星歷中選擇衛星號與
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)
處理過程：
主要都是數學計算，其過程可以在
ICD-GPS-200C P148中找到。
注意事項：
這里計算的電離層延時是相對于 GPS-L1信號而言的，其它頻率信號需要進行一次轉換。
計算過程中很多角度的單位是半圓，即 $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
處理過程：
檢測高度角和接收機高度是否大于閾值。否，則延遲為 0，方差為
VAR_NOTEC，返回 1；是，則進入下一步。
從nav_t.tec中找出第一個tec[i].time>time（輸入參數，信號接收時間）的tec數據。然后通過i==0||i>=nav->nt，確保time是在所給出的nav_t.tec包含的時間段之中！
通過確認所找到的時間段的右端點減去左端點，來確保時間間隔!=0。
調用 iondelay來計算所屬時間段兩端端點的電離層延時。
由兩端的延時，插值計算出觀測時間點處的值。而對于兩端延時的組合，有 3種情況。
① 兩個端點都計算出錯，輸出錯誤信息，返回 0.
② 兩個端點都有值，線性插值出觀測時間點的值，返回 1.
③ 只有一個端點有值，將其結果作為觀測時間處的值，返回 1.
注意事項：
由于是通過調用 iondelay來計算所屬時間段端點的電離層延時，所以這里求出來的值是以
L1信號為前提的。
關于 5中的第 ②種情況，RTKLIB manual P152 (E.5.20)式是錯誤的，左端點TEC值得時間權重值應該是 $(t_{i+1}-t)/(t_{i+1}-t_i)$。manual中可能是搞反了，源碼中是正確的，與我的看法相同。
我的疑惑：
1中當高度角和接收機高度較小時，為什么延遲要為 0呢？
可能是個對最終結果沒有什么影響的小細節， 雖然時間間隔tt后面用得到，但是由于 2中的操作，其實3中的時間間隔肯定是>0的！
目前關于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
處理過程：
按照
RTKLIB manual P152中的公式(E.5.19)，先計算出與頻率有關的中間項。
整體過程是按照電離層的高度，從起始高度開始，逐層計算每一層的延時和方差，之后累加到一起。下面再具體闡述。
首先調用 ionppp函數，計算出在當前電離層高度時，電離層穿刺點的位置{lat,lon,h} (rad,m)和傾斜率( $1/cos zeta$)。
按照opt的值可能再次進行修正。opt&1，則按照M-SLM映射函數重新計算傾斜率；opt&2，則在日固坐標系中考慮地球自轉，重新計算穿刺點經度；
由 TEC網格計算穿刺點的電子數總量，然后按照 (E.5.19)累加電離層延時(m)和方差。
注意事項：
這里在計算電離層延時(m)時，是假設信號為 L1的！
我的疑惑：
4中操作看不懂，貌似跟 地/日坐標系和
SL/M-SL有關？
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 傾斜率
處理過程：
與處理過程相對應的公式，請見
RTKLIB manual P151
注意事項：
說明文檔中的
z并不是仰角azel[1]，而是仰角關于$pi/2$的補角，所以程序中在計算rp是采用的是cos(azel[1])的寫法。
可能因為后面再從 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)
處理過程：
將
value和rms所指向的值置為 0。
檢驗tec的緯度和經度間隔是否為 0。是，則直接返回 0。
將穿刺點的經緯度分別減去網格點的起始經緯度，再除以網格點間距，對結果進行取整，得到穿刺點所在網格的序號和穿刺點所在網格的位置(比例)。
按照下圖的順序，調用dataindex函數分別計算這些網格點的tec數據在tec.data中的下標，從而得到這些網格點處的TEC值和相應誤差的標準差。
如果四個網格點的 TEC值都 >0，則說明穿刺點位于網格內，使用雙線性插值計算出穿刺點的 TEC值；否則使用最鄰近的網格點值作為穿刺點的 TEC值，不過前提是網格點的 TEC>0；否則，選用四個網格點中 >0的值的平均值作為穿刺點的 TEC值。
注意事項：
對于lats[3]，其含義分別為起始緯度、終止緯度和間隔，對 lons[3]、hgts[3]，其含義也是類似的。
對于dataindex函數， i、j、k都是從 0開始的，意味著分別代表各自方向上第 i+1、j+1、k+1層，并且是按照緯度、經度、高度的優先順序來存儲網格點數據的。
我的疑惑：
關于輸出參數
rms，按照其名稱，應該是均方根值，但是在調用了該函數的 iondelay中，確是把它的平方當做方差的一部分進行累加。所以我估計tec.rms值得應該是相應網格點數據值的方差？
老實說，源碼中關于tec->lons[2]大于或者小于 0所做的處理，并沒有看得太明白。另外，個人感覺 3中減去網格點起始經度后的差值應該也不會超過 360°吧？
整體由網格點數據插值穿刺點值得過程可以明白，但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)
處理過程：
與處理過程相對應的公式，請見
RTKLIB manual P149
我的疑惑：
源碼中關于
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
處理過程：
調用 ecef2pos函數，將接收機位置由 ECEF轉換為大地坐標系。
調用 xyz2enu函數，計算此時的坐標轉換矩陣。
滿足下列條件obs[i].D[0]==0.0||lam==0.0||!vsat[i]||norm(rs+3+i*6,3)<=0.0之一，則當前衛星在定速時不可用，直接進行下一次循環。
計算當前接收機位置下 ENU中的視向量，然后轉換得到 ECEF中視向量的值。
計算 ECEF中衛星相對于接收機的速度，然后再計算出考慮了地球自轉的用戶和衛星之間的幾何距離變化率，校正公式見RTKLIB manual P159 (E.6.29)。
根據公式 計算出方程右端項的多普勒殘余，然后再構建左端項的幾何矩陣。最后再將觀測方程數增 1.
注意事項：
這里與定位不同，構建幾何矩陣時，就只有 4個未知數，而定位時是有
NX個。并且沒有像定位那樣為了防止虧秩而進行約束處理。
多普勒定速方程中幾何矩陣 G與定位方程中的一樣，前三行都是 ECEF坐標系中由接收機指向衛星的單位觀測矢量的反向。而由于轉換矩陣 S本身是一個正交單位矩陣($ST=S{-1}$)，所以這里在計算 ECEF中的視向量時，對E進行了轉置處理。
這里構建的左端幾何矩陣 H，也與定位方程中的一樣，是大部分資料上的幾何矩陣的轉置。
總結
以上是生活随笔為你收集整理的RTKLIB源码解析（一）——单点定位(pntpos.c)的全部內容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： Python学习16 正则表达式3 练习
- 下一篇： 17 款程序员必备 Chrome扩展插件