RTKLIB学习总结（六）导航电文、卫星位置计算
生活随笔
收集整理的這篇文章主要介紹了
RTKLIB学习总结（六）导航电文、卫星位置计算
小編覺得挺不錯的,現在分享給大家,幫大家做個參考.
文章目錄
- 一、導航電文
- 1、GNSS衛星信號的組成
- 2、導航電文的編排
- 3、遙測字（TLW）
- 4、交接字（HOW）
- 5、第一數據塊
- 6、第二數據塊
- 7、第三數據塊
- 二、衛星鐘差鐘漂改正
- 1、時鐘校正參數（ a f 0 、 a f 1 、 a f 2 a_{f0}、a_{f1}、a_{f2} af0?、af1?、af2?）改正
- 2、相對論效應校正 Δ t r \Delta t_r Δtr?
- 3、群波延遲校正 T G D T_{GD} TGD?
- 4、鐘漂校正
- 三、satposs()調用流程
- 1、調用流程圖
- 2、傳入參數
- 3、執行流程
- 4、RTKLIB中衛星和衛星系統的表示
- 1.衛星系統
- 2.衛星
- 5、調用函數
- 1.satpos()：單顆衛星計算入口函數
- 2.ephclk()：鐘差解算
- 3.ephpos()：位置解算
- 4.seleph()、selgeph()：選擇用于解算星歷數據
- 四、衛星位置、鐘差的具體計算函數
- 1、alm2pos()：由歷書計算衛星
- 2、eph2clk()：由廣播星歷計算衛星鐘差
- 3、geph2clk()：由廣播星歷計算GLONASS衛星鐘差鐘差
- 4、eph2pos()：由廣播星歷計算衛星位置鐘差
- 5、var_uraeph()：用URA用戶測距精度標定衛星位置方差。
- 6、GLONASS衛星位置計算： geph2pos() -> glorbit() -> deq()
- 1.glorbit()：龍格庫塔迭代
- 2.deq()：微分方程計算
- 五、精密星歷
- 1、精密星歷讀取流程
- 2、peph2pos()：精密星歷計算衛星位置、鐘差、速度、鐘漂
- 2、pephpos()：精密星歷計算衛星位置，鐘差
- 3、interppol()：Neville插值
- 4、pephclk()：精密鐘差計算衛星鐘差
一、導航電文
1、GNSS衛星信號的組成
GPS衛星信號由載波、偽碼、導航電文（數據碼）三個層次組成。數據碼首先與偽碼異或相加實現擴頻，然后二者的組合碼通過雙向移位鍵控（BPSK）對載波進行調制。用戶接收機首先對載波信號進行BPSK（調相調制）解調，使衛星信號的中心頻率從L1下變頻為0，之后再將載波解調后的衛星信號與接收機內部復制的C/A碼Gi做自相關運算，剝離衛星信號中的C/A碼，使信號頻寬變回到只含數據碼的基帶，以得到 50bps 數據碼，再按導航電文的格式最終將數據碼編譯成導航電文。
2、導航電文的編排
- 衛星將導航電文以幀與子幀的形式編排成數據流，每幀導航電文長1500比特，計30s，依次由5個子幀組成。每個子幀長300比特，計6s，依次由10個字組成，每字長30比特，最高為比特先發，每個子幀中的每個字均以6比特的奇偶校驗碼結束。每比特長20ms，C/A碼重復20個周期。
- 每一個子幀的前兩個字分別是遙測字（TLW）、交接字（HOW），后8個字組成數據塊。第1子幀中的數據塊稱為第一數據庫塊，2和3子幀數據塊合稱為第二數據塊，剩下的稱為第三數據塊。當衛星出現故障時，會在各大數據塊的8個字中交替發射1和0。
- GPS對第三數據塊采用了分頁結構，即一幀中的第4子幀和第5子幀為一頁，然后下一幀的第4和第5子幀繼續發送下一頁，第三數據塊內容共25頁，所以發送完一整套導航電文需要750s（12.5min） ，整個導航電文每12.5min重復一次。
- 在一個新GPS周開始時，導航電文從第一子幀開始播發，子幀4和5從第三數據塊的首頁開始播發。子幀123需要更新時，新的導航電文從幀的延邊處開始播發（對應GPS時間是30s的整數倍）。子幀45需要更新時，新導航電文可以在子幀45的任何一頁開始播發。
3、遙測字（TLW）
每一個子幀的第一個字均為遙測字，在導航電文中每6s出現一次，內部組成情況。1~8位是二進制固定值10001011的同步碼，9–22位提供特許用戶所需信息，23、24位保留，最后6位是奇偶校驗碼。
4、交接字（HOW）
1–17位是截短的周內時計數，18位是警告標志，為1時提醒非特許用戶自己承擔使用該衛星信號的風險，該衛星第一數據塊所提供的URA值有可能比其真實值大。19比特A–S標志，其值為1時表示對該衛星實施了反電子欺騙措施。20–22比特是子幀識別標志，共5個有效二進制：001表示第1子幀，010表示第2子幀，依次類推。后面是奇偶校驗碼。
5、第一數據塊
也稱時鐘數據塊，由第一子幀構成，包含如下內容：
-
星期數WN：GPS星期，來自Z計數值的高10比特，最大值1023約19年，若1023+1則返0，上次返0在2019年4月6號。
我國的北斗系統，也存在BD周數翻轉問題，但在設計時，其周計數用13bit表示，翻轉周期是8192周，大概是160年 。
-
用戶測距精度（URA）：對所有由GPS地面監控部分和空間星座部分引起的測距誤差大小的一個統計值，通過4比特表示用戶測距精度因子N，0<=N <=6，URA = 2^(1+N/2)，6<N<15，URA = 2的（N-2）次方， URA值越大，表示衛星信號得到的GPS距離測量值精度越低。
-
衛星健康狀況：6比特，最高位0表示正常1表示出錯，低5位表示具體出錯情況。
-
時鐘校正參數（ a f 0 、 a f 1 、 a f 2 a_{f0}、a_{f1}、a_{f2} af0?、af1?、af2?）：衛星時鐘模型校正方程的三個系數。
-
群波延時校正值（ T G D ? T_{GD}? TGD??）：針對單頻接收機。單頻接收機之所以有此項校正因為 a f 0 ? a_{f0}? af0??是針對雙頻測量值而言的。
-
時鐘數據期號（IODC）：10比特表示的時鐘數據塊期號，一個IODC對應一套時鐘校正參數，可用于快速檢測時鐘參數是否發生改變，IODC不變說明時鐘參數沒更新，就不用再重復讀取。
6、第二數據塊
提供衛星自身的星歷參數、由子幀2 3組成，內容如下：
衛星星歷是描述衛星運動軌道的信息。也可以說衛星星歷就是一組對應某一時刻的軌道參數及其變率。有了衛星星歷就可以計算出任意時刻的衛星位置及其速度。GPS廣播星歷參數共有16個，其中包括1個參考時刻，6個對應參考時刻的開普勒軌道參數和9個反映攝動力影響的參數 ：
7、第三數據塊
子幀4 5組成的第三數據塊提供所有衛星的歷書參數、電離層延時校正參數、GPS時與UTC間的關系及衛星健康狀況等信息。
歷書參數與星歷參數比較：
- 歷書與星歷都是表示衛星運行的參數。歷書包括全部衛星的大概位置，用于衛星預報；星歷只是當前接收機觀測到的衛星的精確位置，用于定位。
- 為了縮短衛星鎖定時間，GPS接收機需利用歷書、當地位置的時間來預報衛星運行狀態。
- 歷書是從導航電文中提取的，每12.5分鐘的導航電文才能得到一組完整的歷書。歷書的有效期為半年。
- 利用歷書和當地的位置， 我們可以計算出衛星的方位和高度角，由此可以計算出當地能觀測到的衛星和持續時間，即衛星高度角大于5°的出現時間。
- GPS衛星星歷參數包含在導航電文的第二和第三子幀中。從有效的星歷中，我們可解得衛星的較準確位置和速度，從而用于接收機定位和測速。GPS衛星星歷每30秒重復一次，有效期為以星歷參考時間為中心的4小時內 。
- 歷書信息存在alm_t 中，星歷存在eph_t 中
二、衛星鐘差鐘漂改正
1、時鐘校正參數（ a f 0 、 a f 1 、 a f 2 a_{f0}、a_{f1}、a_{f2} af0?、af1?、af2?）改正
? 相當于GPS時間，衛星上作為時間和頻率信號來源的原子鐘也存在時間偏差和頻率漂移。為確保各顆衛星的時鐘與GPS時間同步，GPS地面監控部分通過對衛星信號進行檢測，將衛星時鐘在GPS時間t的衛星鐘差 Δ t ( s ) \Delta t^{(s)} Δt(s)描述為如下二項式：
Δ t ( s ) = a f 0 + a f 1 ( t ? t o c ) + a f 2 ( t ? t o c ) 2 \Delta t^{(s)}=a_{f0}+a_{f1}(t-t_{oc})+a_{f2}(t-t_{oc})^2 Δt(s)=af0?+af1?(t?toc?)+af2?(t?toc?)2
2、相對論效應校正 Δ t r \Delta t_r Δtr?
? 綜合狹義相對論和廣義相對論，在高空中高速運行的衛星原子鐘比地面上一模一樣的原子鐘每天要快 38000ns ，每秒快 0.44ns 。如果不考慮相對論效應，GPS發上天兩分鐘內，衛星原子鐘就會失去定位作用。在地面上設計原子鐘時可以減小一點點它的頻率，上天以后其時鐘頻率在地面上看來正好等于設計值。同時因為GPS運行軌道是橢圓而不是圓，地面上計算機還有根據衛星當前位置做相對論效應的校如下：
Δ t r = F e s a s sin ? E k \Delta t_r=Fe_s\sqrt{a_s} \sin E_k Δtr?=Fes?as??sinEk?
3、群波延遲校正 T G D T_{GD} TGD?
? 由第一數據塊給出，只適用于單頻。這樣對于L1單頻接收機，衛星時鐘總鐘差值如下：
δ t ( s ) = Δ t ( s ) + Δ t r ? T G D \delta t^{(s)}=\Delta t^{(s)}+\Delta t_{r}-T_{G D} δt(s)=Δt(s)+Δtr??TGD?
4、鐘漂校正
? 對上面衛星時鐘總鐘差值求導得：
δ f ( s ) = a f 1 + 2 a f 2 ( t ? t o c ) + Δ t ˙ r \delta f^{(s)}=a_{f 1}+2 a_{f 2}\left(t-t_{o c}\right)+\Delta \dot{t}_r δf(s)=af1?+2af2?(t?toc?)+Δt˙r?
? 群波延遲校正 T G D T_{GD} TGD?的導數為0，相對論效應校正 Δ t r \Delta t_r Δtr?如下：
Δ t ˙ r = F e s a s E ˙ k cos ? E k \Delta \dot{t}_r=F e_s \sqrt{a_s} \dot{E}_k \cos E_k Δt˙r?=Fes?as??E˙k?cosEk?
三、satposs()調用流程
1、調用流程圖
2、傳入參數
gtime_t teph I (gpst) 用于選擇星歷的時刻 (gpst) obsd_t *obs I OBS觀測數據 int n I OBS數 nav_t *nav I NAV導航電文 int ephopt I 星歷選項 (EPHOPT_???) double *rs O 衛星位置和速度，長度為6*n，{x,y,z,vx,vy,vz}(ecef)(m,m/s) double *dts O 衛星鐘差，長度為2*n， {bias,drift} (s|s/s) double *var O 衛星位置和鐘差的協方差 (m^2) int *svh O 衛星健康標志 (-1:correction not available)3、執行流程
rs [(0:2)+i*6]= obs[i] sat position {x,y,z} (m) 衛星位置 rs [(3:5)+i*6]= obs[i] sat velocity {vx,vy,vz} (m/s) 衛星速度 dts[(0:1)+i*2]= obs[i] sat clock {bias,drift} (s|s/s) 衛星鐘差鐘漂- 遍歷每一個OBS觀測數據：
- 首先初始化，將對當前觀測數據的 rs、dts、var和svh數組的元素置 0
- 通過判斷某一頻率下信號的偽距是否為 0，來得到此時所用的頻率個數j
- 用數據接收時刻減去偽距信號傳播時間，得到衛星信號的發射時刻 time[i]
- 調用 ephclk() 函數，由廣播星歷計算出當前觀測衛星與 GPS 時間的鐘差 dt ,此時的鐘差是沒有考慮相對論效應和 TGD 的 ，dt僅作為satpos()的參數，不作為最終計算的鐘差。
- 信號發射時刻減去鐘差 dt，得到 GPS 時間下的衛星信號發射時刻
- 調用 satpos() 函數，計算信號發射時刻衛星的位置(ecef,m)、速度(ecef,m/s)、鐘差((s|s/s)) ,這里計算出的鐘差是考慮了相對論效應的了，只是還沒有考慮 TGD
- 如果沒有精密鐘差，則ephclk()用廣播星歷的鐘差替代，猜測可能是選了EPHOPT_PREC 精密星歷，但精密鐘差計算出問題，用廣播星歷重新計算鐘差。
4、RTKLIB中衛星和衛星系統的表示
1.衛星系統
-
表示衛星系統的字母：G：GPS、R：GLONASS、E：GALILEO、C：BDS、J：QZSS，I：IRNSS、S：SBAS
-
7位二進制碼表示，對應位寫1表示有對應的系統，做與運算可加系統。
static const int navsys[]={ /* satellite systems */SYS_GPS,SYS_GLO,SYS_GAL,SYS_QZS,SYS_SBS,SYS_CMP,SYS_IRN,0 }; #define SYS_NONE 0x00 /* navigation system: none */ #define SYS_GPS 0x01 /* navigation system: GPS */ #define SYS_SBS 0x02 /* navigation system: SBAS */ #define SYS_GLO 0x04 /* navigation system: GLONASS */ #define SYS_GAL 0x08 /* navigation system: Galileo */ #define SYS_QZS 0x10 /* navigation system: QZSS */ #define SYS_CMP 0x20 /* navigation system: BeiDou */ #define SYS_IRN 0x40 /* navigation system: IRNS */ #define SYS_LEO 0x80 /* navigation system: LEO */ #define SYS_ALL 0xFF /* navigation system: all */
2.衛星
可以表示為各系統的衛星ID（系統縮寫+PRN）：B02、C21，也可表示為連續的satellite number ，各種轉換函數如下：
- satno()：傳入衛星系統(SYS_GPS,SYS_GLO,…) ，和PRN碼，轉換為連續的satellite number。
- satsys()：傳入satellite number ，返回衛星系統(SYS_GPS,SYS_GLO,…) ，通過傳入的指針prn傳出PRN值。
- satid2no()：傳入衛星ID，返回satellite number。
- satno2id()：傳入衛星系統，和PRN，返回衛星ID(Gxx,Cxx)
- sat2code()：傳入satellite number，返回衛星ID(Gxx,Cxx)
- code2sys()：傳入衛星系統縮寫，返回系統二進制碼SYS_XXX。
- satexclude()：檢測某顆衛星在定位時是否需要將其排除
5、調用函數
1.satpos()：單顆衛星計算入口函數
針對單顆衛星，根據星歷選項調用對應的解算函數ephpos()、satpos_sbas()、satpos_ssr()、peph2pos()，解算信號發射時刻衛星的 P(ecef,m)、V(ecef,m/s)、C((s|s/s))
extern int satpos(gtime_t time, gtime_t teph, int sat, int ephopt,const nav_t *nav, double *rs, double *dts, double *var,int *svh) {trace(4,"satpos : time=%s sat=%2d ephopt=%d\n",time_str(time,3),sat,ephopt);*svh=0;switch (ephopt) { //根據星歷選項調用對應的解算函數case EPHOPT_BRDC : return ephpos (time,teph,sat,nav,-1,rs,dts,var,svh); //廣播星歷case EPHOPT_SBAS : return satpos_sbas(time,teph,sat,nav, rs,dts,var,svh); //sbascase EPHOPT_SSRAPC: return satpos_ssr (time,teph,sat,nav, 0,rs,dts,var,svh); //參考天線相位中心case EPHOPT_SSRCOM: return satpos_ssr (time,teph,sat,nav, 1,rs,dts,var,svh); //參考質心，還需要天線相位中心改正case EPHOPT_PREC : //精密星歷if (!peph2pos(time,sat,nav,1,rs,dts,var)) break; else return 1;}*svh=-1;return 0; }2.ephclk()：鐘差解算
根據衛星系統調用對應的函數eph2clk()、geph2clk()、seph2clk()，通過廣播星歷來確定衛星鐘差 、鐘漂
static int ephclk(gtime_t time, gtime_t teph, int sat, const nav_t *nav,double *dts) {eph_t *eph;geph_t *geph;seph_t *seph;int sys;trace(4,"ephclk : time=%s sat=%2d\n",time_str(time,3),sat);//調用 satsys 函數，根據衛星編號確定該衛星所屬的導航系統和該衛星在該系統中的 PRN編號sys=satsys(sat,NULL);if (sys==SYS_GPS||sys==SYS_GAL||sys==SYS_QZS||sys==SYS_CMP||sys==SYS_IRN) {if (!(eph=seleph(teph,sat,-1,nav))) return 0; //調用 seleph 函數來選擇最接近 teph 的那個星*dts=eph2clk(time,eph); //調用 eph2clk 函數，通過廣播星歷和信號發射時間計算出衛星鐘差}else if (sys==SYS_GLO) {if (!(geph=selgeph(teph,sat,-1,nav))) return 0;*dts=geph2clk(time,geph);}else if (sys==SYS_SBS) {if (!(seph=selseph(teph,sat,nav))) return 0;*dts=seph2clk(time,seph);}else return 0;return 1; }3.ephpos()：位置解算
由satpos()調用，執行流程如下：
- 根據衛星系統，先調用對應的星歷選擇函數，seleph()、selgeph()、seph2pos()。
- 調用對應的解算函數，eph2pos()、geph2pos()、seph2pos()，計算位置、鐘差。
- 增加一個極短的時間tt，再調用對應的解算函數計算位置、鐘差。
- 兩次的位置、鐘差相減再除以tt，得速度、鐘漂。
4.seleph()、selgeph()：選擇用于解算星歷數據
傳入sattle number，時間或IDOE，如果傳入IDOE>=0,按IDOE找星歷數據，否則取最接近時間的星歷
static eph_t *seleph(gtime_t time, int sat, int iode, const nav_t *nav) {double t,tmax,tmin;int i,j=-1,sys,sel;trace(4,"seleph : time=%s sat=%2d iode=%d\n",time_str(time,3),sat,iode);//根據傳入的sattle number，調用satsys()判斷衛星系統，賦值tmax，tmin，selsys=satsys(sat,NULL);switch (sys) {case SYS_GPS: tmax=MAXDTOE+1.0 ; sel=eph_sel[0]; break;case SYS_GAL: tmax=MAXDTOE_GAL ; sel=eph_sel[2]; break;case SYS_QZS: tmax=MAXDTOE_QZS+1.0; sel=eph_sel[3]; break;case SYS_CMP: tmax=MAXDTOE_CMP+1.0; sel=eph_sel[4]; break;case SYS_IRN: tmax=MAXDTOE_IRN+1.0; sel=eph_sel[5]; break;default: tmax=MAXDTOE+1.0; break;}tmin=tmax+1.0;//遍歷nav->eph[]for (i=0;i<nav->n;i++) {if (nav->eph[i].sat!=sat) continue; //eph[i]不是需要的衛星，進入下一次循環if (iode>=0&&nav->eph[i].iode!=iode) continue; //如果傳入了idoe時間不符合，也進行下一次循環if (sys==SYS_GAL) { //若是伽利略衛星還要判斷是I/NAV、F/NAVsel=getseleph(SYS_GAL);if (sel==0&&!(nav->eph[i].code&(1<<9))) continue; /* I/NAV */if (sel==1&&!(nav->eph[i].code&(1<<8))) continue; /* F/NAV */if (timediff(nav->eph[i].toe,time)>=0.0) continue; /* AOD<=0 */}if ((t=fabs(timediff(nav->eph[i].toe,time)))>tmax) continue; //時間差過大，也進行下一次循化if (iode>=0) return nav->eph+i; //傳入IDOE>=0，直接返回符合條件的星歷if (t<=tmin) {j=i; tmin=t;} /* toe closest to time */ //存下時間差最小的星歷下標和時間差}if (iode>=0||j<0) {trace(3,"no broadcast ephemeris: %s sat=%2d iode=%3d\n",time_str(time,0),sat,iode);return NULL;}return nav->eph+j; } static geph_t *selgeph(gtime_t time, int sat, int iode, const nav_t *nav) {double t,tmax=MAXDTOE_GLO,tmin=tmax+1.0;int i,j=-1;trace(4,"selgeph : time=%s sat=%2d iode=%2d\n",time_str(time,3),sat,iode);//變量nav->geph[]，for (i=0;i<nav->ng;i++) {if (nav->geph[i].sat!=sat) continue; //衛星不同，進行下一次循環if (iode>=0&&nav->geph[i].iode!=iode) continue; //傳入IDOE>0，IDOE不同，進行下一次循環if ((t=fabs(timediff(nav->geph[i].toe,time)))>tmax) continue;if (iode>=0) return nav->geph+i; //傳入IDOE>=0，直接返回符合條件的星歷if (t<=tmin) {j=i; tmin=t;} /* toe closest to time */ //存下時間差最小的星歷下標和時間差}if (iode>=0||j<0) {trace(3,"no glonass ephemeris : %s sat=%2d iode=%2d\n",time_str(time,0),sat,iode);return NULL;}return nav->geph+j; }四、衛星位置、鐘差的具體計算函數
公式模型見manual的142面，我在對應的函數語句后面標注了對應的公式編號。
1、alm2pos()：由歷書計算衛星
歷書信息計算衛星位置、鐘差。（這個函數在RTKLIB里好像沒怎么被調用過）
extern void alm2pos(gtime_t time, const alm_t *alm, double *rs, double *dts) {double tk,M,E,Ek,sinE,cosE,u,r,i,O,x,y,sinO,cosO,cosi,mu;int n;trace(4,"alm2pos : time=%s sat=%2d\n",time_str(time,3),alm->sat);tk=timediff(time,alm->toa);if (alm->A<=0.0) {rs[0]=rs[1]=rs[2]=*dts=0.0;return;}mu=satsys(alm->sat,NULL)==SYS_GAL?MU_GAL:MU_GPS;M=alm->M0+sqrt(mu/(alm->A*alm->A*alm->A))*tk;for (n=0,E=M,Ek=0.0;fabs(E-Ek)>RTOL_KEPLER&&n<MAX_ITER_KEPLER;n++) {Ek=E; E-=(E-alm->e*sin(E)-M)/(1.0-alm->e*cos(E));}if (n>=MAX_ITER_KEPLER) {trace(2,"alm2pos: kepler iteration overflow sat=%2d\n",alm->sat);return;}sinE=sin(E); cosE=cos(E);u=atan2(sqrt(1.0-alm->e*alm->e)*sinE,cosE-alm->e)+alm->omg;r=alm->A*(1.0-alm->e*cosE);i=alm->i0;O=alm->OMG0+(alm->OMGd-OMGE)*tk-OMGE*alm->toas;x=r*cos(u); y=r*sin(u); sinO=sin(O); cosO=cos(O); cosi=cos(i);rs[0]=x*cosO-y*cosi*sinO;rs[1]=x*sinO+y*cosi*cosO;rs[2]=y*sin(i);*dts=alm->f0+alm->f1*tk; }2、eph2clk()：由廣播星歷計算衛星鐘差
根據信號發射時間和廣播星歷，計算衛星鐘差,不考慮相對論效應和TGD ，二項式校正公式用了三遍。
3、geph2clk()：由廣播星歷計算GLONASS衛星鐘差鐘差
與eph2clk() 類似
extern double eph2clk(gtime_t time, const eph_t *eph) {double t,ts;int i;trace(4,"eph2clk : time=%s sat=%2d\n",time_str(time,3),eph->sat);t=ts=timediff(time,eph->toc); //計算與星歷參考時間的偏差 dt = t-toc//利用二項式校正計算出衛星鐘差，從 dt中減去這部分，然后再進行一次上述操作，得到最終的 dtfor (i=0;i<2;i++) {t=ts-(eph->f0+eph->f1*t+eph->f2*t*t); //(E.4.16)}//使用二項式校正得到最終的鐘差return eph->f0+eph->f1*t+eph->f2*t*t; } extern double geph2clk(gtime_t time, const geph_t *geph) {double t,ts;int i;trace(4,"geph2clk: time=%s sat=%2d\n",time_str(time,3),geph->sat);t=ts=timediff(time,geph->toe);for (i=0;i<2;i++) {t=ts-(-geph->taun+geph->gamn*t); //(E.4.26)}return -geph->taun+geph->gamn*t; }4、eph2pos()：由廣播星歷計算衛星位置鐘差
根據廣播星歷計算出算信號發射時刻衛星的位置和鐘差 ，gps, galileo, qzss, bds，由開普勒軌道參數和攝動改正參數計算，考慮的相對論效應，但沒考慮TGD。
extern void eph2pos(gtime_t time, const eph_t *eph, double *rs, double *dts,double *var) {double tk,M,E,Ek,sinE,cosE,u,r,i,O,sin2u,cos2u,x,y,sinO,cosO,cosi,mu,omge;double xg,yg,zg,sino,coso;int n,sys,prn;trace(4,"eph2pos : time=%s sat=%2d\n",time_str(time,3),eph->sat);if (eph->A<=0.0) { //通過衛星軌道半長軸 A 判斷星歷是否有效，無效則返回rs[0]=rs[1]=rs[2]=*dts=*var=0.0;return;}tk=timediff(time,eph->toe); //計算規化時間 tk (E.4.2)switch ((sys=satsys(eph->sat,&prn))) { //根據不同衛星系統設置相應的地球引力常數 mu 和 地球自轉角速度 omgecase SYS_GAL: mu=MU_GAL; omge=OMGE_GAL; break;case SYS_CMP: mu=MU_CMP; omge=OMGE_CMP; break;default: mu=MU_GPS; omge=OMGE; break;}M=eph->M0+(sqrt(mu/(eph->A*eph->A*eph->A))+eph->deln)*tk; //計算平近點角 M (E.4.3)//用牛頓迭代法來計算偏近點角 E。參考 RTKLIB manual P145 (E.4.19) (E.4.4)for (n=0,E=M,Ek=0.0;fabs(E-Ek)>RTOL_KEPLER&&n<MAX_ITER_KEPLER;n++) {Ek=E; E-=(E-eph->e*sin(E)-M)/(1.0-eph->e*cos(E));}if (n>=MAX_ITER_KEPLER) {trace(2,"eph2pos: kepler iteration overflow sat=%2d\n",eph->sat);return;}sinE=sin(E); cosE=cos(E);trace(4,"kepler: sat=%2d e=%8.5f n=%2d del=%10.3e\n",eph->sat,eph->e,n,E-Ek);//計算攝動改正后的 升交點角距u 衛星矢徑長度r 軌道傾角iu=atan2(sqrt(1.0-eph->e*eph->e)*sinE,cosE-eph->e)+eph->omg; //(E.4.5) (E.4.6) (E.4.10)r=eph->A*(1.0-eph->e*cosE); //(E.4.11)i=eph->i0+eph->idot*tk; //(E.4.12)sin2u=sin(2.0*u); cos2u=cos(2.0*u); u+=eph->cus*sin2u+eph->cuc*cos2u; //(E.4.7)r+=eph->crs*sin2u+eph->crc*cos2u; //(E.4.8)i+=eph->cis*sin2u+eph->cic*cos2u; //(E.4.9)x=r*cos(u); y=r*sin(u); cosi=cos(i);//北斗的MEO、IGSO衛星計算方法與GPS, Galileo and QZSS相同，只是一些參數不同//GEO衛星的 O 和最后位置的計算稍有不同 /* beidou geo satellite */if (sys==SYS_CMP&&(prn<=5||prn>=59)) { /* ref [9] table 4-1 */O=eph->OMG0+eph->OMGd*tk-omge*eph->toes; //(E.4.29)sinO=sin(O); cosO=cos(O);xg=x*cosO-y*cosi*sinO;yg=x*sinO+y*cosi*cosO;zg=y*sin(i);sino=sin(omge*tk); coso=cos(omge*tk);rs[0]= xg*coso + yg*sino*COS_5 + zg*sino*SIN_5; //ECEF位置(E.4.30)rs[1]=-xg*sino + yg*coso*COS_5 + zg*coso*SIN_5;rs[2]=-yg*SIN_5 + zg*COS_5;}else {O=eph->OMG0+(eph->OMGd-omge)*tk-omge*eph->toes; //計算升交點赤經O (E.4.13)sinO=sin(O); cosO=cos(O);rs[0]=x*cosO-y*cosi*sinO; //計算衛星ECEF位置存入 rs 中 (E.4.14)rs[1]=x*sinO+y*cosi*cosO;rs[2]=y*sin(i);}tk=timediff(time,eph->toc); //(E.4.15)*dts=eph->f0+eph->f1*tk+eph->f2*tk*tk; //利用三個二項式模型系數af0、af1、af2計算衛星鐘差/* relativity correction */ *dts-=2.0*sqrt(mu*eph->A)*eph->e*sinE/SQR(CLIGHT); //相對論效應改正衛星鐘差/* position and clock error variance */*var=var_uraeph(sys,eph->sva); //用 URA 值來標定方差 }5、var_uraeph()：用URA用戶測距精度標定衛星位置方差。
static double var_uraeph(int sys, int ura) {const double ura_value[]={ 2.4,3.4,4.85,6.85,9.65,13.65,24.0,48.0,96.0,192.0,384.0,768.0,1536.0,3072.0,6144.0};if (sys==SYS_GAL) { /* galileo sisa (ref [7] 5.1.11) */if (ura<= 49) return SQR(ura*0.01);if (ura<= 74) return SQR(0.5+(ura- 50)*0.02);if (ura<= 99) return SQR(1.0+(ura- 75)*0.04);if (ura<=125) return SQR(2.0+(ura-100)*0.16);return SQR(STD_GAL_NAPA);}else { /* gps ura (ref [1] 20.3.3.3.1.1) */return ura<0||14<ura?SQR(6144.0):SQR(ura_value[ura]);} }GLONASS URA設置為常數5m
extern void geph2pos(gtime_t time, const geph_t *geph, double *rs, double *dts,double *var) { ... *var=SQR(ERREPH_GLO); }6、GLONASS衛星位置計算： geph2pos() -> glorbit() -> deq()
由GLONASS星歷計算衛星坐標。GLONASS衛星播發的是PZ-90坐標系下參考時刻的衛星狀態向量，每半個小時廣播一次。如果需要得到某個時間的衛星位置必須通過運動模型積分得到。
1.glorbit()：龍格庫塔迭代
2.deq()：微分方程計算
```c extern void geph2pos(gtime_t time, const geph_t *geph, double *rs, double *dts,double *var) {double t,tt,x[6];int i;trace(4,"geph2pos: time=%s sat=%2d\n",time_str(time,3),geph->sat);t=timediff(time,geph->toe);*dts=-geph->taun+geph->gamn*t; //計算鐘差dts(E.4.26)for (i=0;i<3;i++) {x[i ]=geph->pos[i];x[i+3]=geph->vel[i];}//步長TSTEP:60sfor (tt=t<0.0?-TSTEP:TSTEP;fabs(t)>1E-9;t-=tt) {if (fabs(t)<TSTEP) tt=t;glorbit(tt,x,geph->acc);}for (i=0;i<3;i++) rs[i]=x[i];*var=SQR(ERREPH_GLO); //glonass衛星的方差直接定為 5*5 } static void glorbit(double t, double *x, const double *acc) {double k1[6],k2[6],k3[6],k4[6],w[6];int i;deq(x,k1,acc); for (i=0;i<6;i++) w[i]=x[i]+k1[i]*t/2.0;deq(w,k2,acc); for (i=0;i<6;i++) w[i]=x[i]+k2[i]*t/2.0;deq(w,k3,acc); for (i=0;i<6;i++) w[i]=x[i]+k3[i]*t;deq(w,k4,acc); for (i=0;i<6;i++) x[i]+=(k1[i]+2.0*k2[i]+2.0*k3[i]+k4[i])*t/6.0; } static void deq(const double *x, double *xdot, const double *acc) {double a,b,c,r2=dot(x,x,3), //r平方r3=r2*sqrt(r2), //r三次方omg2=SQR(OMGE_GLO); //omg平方if (r2<=0.0) { //計算出錯xdot[0]=xdot[1]=xdot[2]=xdot[3]=xdot[4]=xdot[5]=0.0;return;}/* ref [2] A.3.1.2 with bug fix for xdot[4],xdot[5] */a=1.5*J2_GLO*MU_GLO*SQR(RE_GLO)/r2/r3; /* 3/2*J2*mu*Ae^2/r^5 */b=5.0*x[2]*x[2]/r2; /* 5*z^2/r^2 */c=-MU_GLO/r3-a*(1.0-b); /* -mu/r^3-a(1-b) */xdot[0]=x[3]; xdot[1]=x[4]; xdot[2]=x[5]; //(E.4.22)xdot[3]=(c+omg2)*x[0]+2.0*OMGE_GLO*x[4]+acc[0]; //(E.4.23)xdot[4]=(c+omg2)*x[1]-2.0*OMGE_GLO*x[3]+acc[1]; //(E.4.24)xdot[5]=(c-2.0*a)*x[2]+acc[2]; //(E.4.25) } ```五、精密星歷
1、精密星歷讀取流程
nav->peph[]存精密星歷數據，nav->ne精密鐘差數量。
nav->pclk[]存精密鐘差數據，nav->nc精密鐘差數量。
- execses_b()中調用readpreceph()。
- readpreceph()中：readsp3()讀取精密星歷，readrnxc()讀取精密鐘差
- readsp3()中：readsp3h()讀文件頭，readsp3b()讀文件體，combpeph()對精密星歷按時間、index排序，再將相同星歷合并。
- readrnxc()中：readrnxfile()讀取精密星歷文件，combpclk()排序合并精密鐘差。
2、peph2pos()：精密星歷計算衛星位置、鐘差、速度、鐘漂
執行流程如下：
- 調用pephpos()、pephclk() ，計算衛星位置、鐘差。
- 增加一個極短的時間tt再計算衛星位置、鐘差。
- 調用satantoff()計算天線相位中心改正量dant[]。
- 衛星位置rs[i]：pephpos()計算值+satantoff()天線相位中心改正 。
- 衛星速度rs[i+3]：兩次的位置相減再除以tt。
- 鐘差dts[0]：相對論效應改正 。
- 鐘漂dts[1]：兩次的鐘差相減再除以tt
- 沒有精密鐘差，dts賦值0，satposs()中會ephclk()廣播星歷的鐘差替代
2、pephpos()：精密星歷計算衛星位置，鐘差
static int pephpos(gtime_t time, int sat, const nav_t *nav, double *rs,double *dts, double *vare, double *varc) {double t[NMAX+1],p[3][NMAX+1],c[2],*pos,std=0.0,s[3],sinl,cosl;int i,j,k,index;trace(4,"pephpos : time=%s sat=%2d\n",time_str(time,3),sat);rs[0]=rs[1]=rs[2]=dts[0]=0.0;if (nav->ne<NMAX+1|| //如果時間早于第一個精密星歷時間，或遲于最后一個超過15分鐘，return 0timediff(time,nav->peph[0].time)<-MAXDTE||timediff(time,nav->peph[nav->ne-1].time)>MAXDTE) {trace(3,"no prec ephem %s sat=%2d\n",time_str(time,0),sat);return 0;}/* binary search */ //二分查找nav->peph[]中時間差最接近的精密星歷的下標indexfor (i=0,j=nav->ne-1;i<j;) {k=(i+j)/2;if (timediff(nav->peph[k].time,time)<0.0) i=k+1; else j=k;}index=i<=0?0:i-1;/* polynomial interpolation for orbit */ //軌道多項式插值i=index-(NMAX+1)/2; if (i<0) i=0; else if (i+NMAX>=nav->ne) i=nav->ne-NMAX-1;for (j=0;j<=NMAX;j++) {t[j]=timediff(nav->peph[i+j].time,time);if (norm(nav->peph[i+j].pos[sat-1],3)<=0.0) {trace(3,"prec ephem outage %s sat=%2d\n",time_str(time,0),sat);return 0;}}for (j=0;j<=NMAX;j++) {pos=nav->peph[i+j].pos[sat-1];/* correciton for earh rotation ver.2.4.0 */ //地球自轉改正sinl=sin(OMGE*t[j]);cosl=cos(OMGE*t[j]);p[0][j]=cosl*pos[0]-sinl*pos[1];p[1][j]=sinl*pos[0]+cosl*pos[1];p[2][j]=pos[2];}for (i=0;i<3;i++) {rs[i]=interppol(t,p[i],NMAX+1); //內維爾多項式插值獲取衛星位置}if (vare) {for (i=0;i<3;i++) s[i]=nav->peph[index].std[sat-1][i];std=norm(s,3);/* extrapolation error for orbit */if (t[0 ]>0.0) std+=EXTERR_EPH*SQR(t[0 ])/2.0;else if (t[NMAX]<0.0) std+=EXTERR_EPH*SQR(t[NMAX])/2.0;*vare=SQR(std);}/* linear interpolation for clock */ //線性插值獲取鐘差t[0]=timediff(time,nav->peph[index ].time);t[1]=timediff(time,nav->peph[index+1].time);c[0]=nav->peph[index ].pos[sat-1][3];c[1]=nav->peph[index+1].pos[sat-1][3];//計算標準差if (t[0]<=0.0) {if ((dts[0]=c[0])!=0.0) {std=nav->peph[index].std[sat-1][3]*CLIGHT-EXTERR_CLK*t[0];}}else if (t[1]>=0.0) {if ((dts[0]=c[1])!=0.0) {std=nav->peph[index+1].std[sat-1][3]*CLIGHT+EXTERR_CLK*t[1];}}else if (c[0]!=0.0&&c[1]!=0.0) {dts[0]=(c[1]*t[0]-c[0]*t[1])/(t[0]-t[1]);i=t[0]<-t[1]?0:1;std=nav->peph[index+i].std[sat-1][3]+EXTERR_CLK*fabs(t[i]);}else {dts[0]=0.0;}if (varc) *varc=SQR(std);return 1; }3、interppol()：Neville插值
由兩個n-1次插值多項式構造一個n次多項式的線性逐次插值方法
static double interppol(const double *x, double *y, int n) {int i,j;for (j=1;j<n;j++) {for (i=0;i<n-j;i++) {y[i]=(x[i+j]*y[i]-x[i]*y[i+1])/(x[i+j]-x[i]);}}return y[0]; }4、pephclk()：精密鐘差計算衛星鐘差
static int pephclk(gtime_t time, int sat, const nav_t *nav, double *dts,double *varc) {double t[2],c[2],std;int i,j,k,index;trace(4,"pephclk : time=%s sat=%2d\n",time_str(time,3),sat);if (nav->nc<2||timediff(time,nav->pclk[0].time)<-MAXDTE||timediff(time,nav->pclk[nav->nc-1].time)>MAXDTE) {trace(3,"no prec clock %s sat=%2d\n",time_str(time,0),sat);return 1;}/* binary search */ //二分查找nav->peph[]中時間差最接近的精密星歷的下標indexfor (i=0,j=nav->nc-1;i<j;) {k=(i+j)/2;if (timediff(nav->pclk[k].time,time)<0.0) i=k+1; else j=k;}index=i<=0?0:i-1;/* linear interpolation for clock */ //鐘差線性插值t[0]=timediff(time,nav->pclk[index ].time);t[1]=timediff(time,nav->pclk[index+1].time);c[0]=nav->pclk[index ].clk[sat-1][0];c[1]=nav->pclk[index+1].clk[sat-1][0];if (t[0]<=0.0) {if ((dts[0]=c[0])==0.0) return 0;std=nav->pclk[index].std[sat-1][0]*CLIGHT-EXTERR_CLK*t[0];}else if (t[1]>=0.0) {if ((dts[0]=c[1])==0.0) return 0;std=nav->pclk[index+1].std[sat-1][0]*CLIGHT+EXTERR_CLK*t[1];}else if (c[0]!=0.0&&c[1]!=0.0) {dts[0]=(c[1]*t[0]-c[0]*t[1])/(t[0]-t[1]);i=t[0]<-t[1]?0:1;std=nav->pclk[index+i].std[sat-1][0]*CLIGHT+EXTERR_CLK*fabs(t[i]);}else {trace(3,"prec clock outage %s sat=%2d\n",time_str(time,0),sat);return 0;}if (varc) *varc=SQR(std);return 1; }[0]=timediff(time,nav->pclk[index ].time);t[1]=timediff(time,nav->pclk[index+1].time);c[0]=nav->pclk[index ].clk[sat-1][0];c[1]=nav->pclk[index+1].clk[sat-1][0];if (t[0]<=0.0) {if ((dts[0]=c[0])==0.0) return 0;std=nav->pclk[index].std[sat-1][0]*CLIGHT-EXTERR_CLK*t[0];}else if (t[1]>=0.0) {if ((dts[0]=c[1])==0.0) return 0;std=nav->pclk[index+1].std[sat-1][0]*CLIGHT+EXTERR_CLK*t[1];}else if (c[0]!=0.0&&c[1]!=0.0) {dts[0]=(c[1]*t[0]-c[0]*t[1])/(t[0]-t[1]);i=t[0]<-t[1]?0:1;std=nav->pclk[index+i].std[sat-1][0]*CLIGHT+EXTERR_CLK*fabs(t[i]);}else {trace(3,"prec clock outage %s sat=%2d\n",time_str(time,0),sat);return 0;}if (varc) *varc=SQR(std);return 1; }總結
以上是生活随笔為你收集整理的RTKLIB学习总结（六）导航电文、卫星位置计算的全部內容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： 送崔九
- 下一篇： 有些MP4只有音频没有视频的解决办法