幫“不想偏航”轉(zhuǎn)發(fā)的專欄《RTKLIB代碼學習》
記錄一下這倆禮拜從零開始學習RTKLIB的一些經(jīng)驗,據(jù)說RTKLIB是功能好用,代碼難讀。如能幫助新入門的同道最好不過。如果哪里不對還請大神不吝在評論區(qū)指正。
調(diào)試代碼使用工具為VS2019,RTKLIB版本為2.4.3 b33,是2018.1.30最后更新的源代碼。
下載方面,國內(nèi)從GitHub上下載比較慢,碼云(Gitee)上有源碼和軟件包,下載比較快。
一、編譯運行
RTKLIB的編譯基本參照這個CSDN上的分享完成的https://blog.csdn.net/sd28you28/article/details/82911273
通常來說,大多數(shù)國內(nèi)學習者都是將rnx2rtkp.c作為主函數(shù),在此基礎(chǔ)上進行修改和運行的。路徑為..\RTKLIB\app\rnx2rtkp\rnx2rtkp.c? ? ? ? 在這個路徑里,msc文件夾里有一個msc.sln,這就是VS的解決方案,我們幾乎可以直接使用了。當然,穩(wěn)妥起見,最好把整個rtklib都備份一下,以免改著改著就崩了還忘了怎么改回去。打開后,如果你電腦的VS版本更高,通常會彈出升級工具包之類的窗口,按指示升級即可。以下走的是另外一條路,可供參考。
新建一個項目去做也可,我個人是新建了項目命名為PPPtest,然后拉進來了如下圖所示的文件:

新建項目PPPtest引入的文件
然后選擇Debug模式,在右邊的解決方案資源管理器中,源文件下找到rnx2rtkp.c,雙擊點開,我是按照https://blog.csdn.net/sd28you28/article/details/82911273的指示把程序調(diào)到可以運行的。感謝原作者,在此不再贅述。
rnx2rtkp.c可作修改,下圖是我試圖作精密單點定位(precise point position)時做的一些修改:新增了字符串指針optfile指向由rtkpost.exe生成的conf文件,argc設(shè)置為1跳過了下面很多的循環(huán),然后選擇了定位模式(prcopt.mode)、衛(wèi)星系統(tǒng)(prcopt.navsys),頻點選擇(prcopt.nf ,3為L1+L2+L5),導入了觀測數(shù)據(jù)(.obs)、廣播星歷(.nav)、精密星歷(.sp3),精密鐘差(.clk),注意要ret=postpos...前面的4行if判定,因為之前的一大串把n賦為0了,留著最后的if會直接return -2,后面就沒法運行了。

rnx2rtkp.c

一些關(guān)于PPP的修改
在運行中,由于我的星歷文件里只有北斗的數(shù)據(jù),運行時總是說no nav data就退出了,通過逐步調(diào)試發(fā)現(xiàn),問題在MINPRNCMP和MINPRNMAX都賦了0(rtklib.h),在【項目->屬性->C/C++->預(yù)處理器->預(yù)處理器定義->編輯】里,加上【ENACMP】(跟之前鏈接里處理格洛納斯的報錯的方式一樣),就啟用了北斗的相關(guān)數(shù)據(jù)。然后就可以運行了。
但是處理數(shù)據(jù)時發(fā)現(xiàn),同樣的觀測數(shù)據(jù)、星歷、設(shè)置,使用源代碼運行的結(jié)果和使用rtkpost.exe運行的結(jié)果有差異,可見下圖:

源代碼運行結(jié)果和rtkpost.exe結(jié)果有差異
查看修改時間可以知道,rtkpost.exe更新時間是2019.8.19,而源碼在2018年初就不再更新了,中間做了什么修改不得而知,亦或是更新的代碼放在別處了?
二、精密星歷讀取
精密星歷是區(qū)別于廣播星歷的星歷文件,許多研究中心都通過獲取世界各地的測站的數(shù)據(jù),分析計算推出了各自的精密星歷/鐘差產(chǎn)品。就我所下載少量個精密星歷而言,通常對于GPS的精密星歷,各類產(chǎn)品的差異都比較小,定位誤差大約在厘米級到分米級,鐘差差異在幾個納秒(ns)左右。而對于北斗,差異是比較大的,對于北斗的MEO衛(wèi)星,精密星歷的定位誤差通常都比較小,與GPS類似,但是IGSO就要差一些,GEO差異最大,在米級左右,相隔36000公里,一動不動的定位著實不易。鐘差的差異通常存在系統(tǒng)性誤差,來源大概是各機構(gòu)的時間基準等不一致。這個系統(tǒng)誤差精密單點定位中是會被接收機鐘差吸收掉的,因此影響有限。
精密星歷讀取的調(diào)用流程為:main里的postpos->(openses中讀取一些修正參數(shù),如天線相位中心修正等)->execses_b->readpreceph->readsp3,在readsp3中讀取精密星歷和鐘差,賦值到一個很大的結(jié)構(gòu)體nav中,這個nav指向的是一個靜態(tài)全局結(jié)構(gòu)體navs,閱讀它的定義會發(fā)現(xiàn)。。簡直包羅萬象,讓我感覺之前自己寫的程序都是幼兒園水平。
readsp3中,readsp3h函數(shù)讀取文件頭,第一行str2time把日期轉(zhuǎn)換為它的時間計數(shù)方式,rtklib的時間計數(shù)既不是GPST也不是utc,而是把1970年1月1日0時0分0秒作為起點,不閏秒不跳秒。

文件preceph.c中的readsp3h函數(shù)
82行處有所修改,原來是ns=(int)str2num(buff,4,2),這是因為近年導航衛(wèi)星發(fā)射頻繁,原來只要兩位數(shù)就夠用,現(xiàn)在已經(jīng)達到一百多顆衛(wèi)星,需要三位數(shù)存儲了。事實上,訪問北斗的一些官方網(wǎng)站,比如這個【http://www.csno-tarc.cn/system/constellation】。
84-88行是讀取精密星歷播發(fā)的衛(wèi)星編號,每3個字符一顆衛(wèi)星,94-95是讀取星歷中位置和鐘差的標準差參數(shù),進行誤差估計與預(yù)報。但是這是為只有GPS的sp3文件設(shè)計的,不同的sp3文件里,這個參數(shù)不一定在第15行,這時就讀不了了。
(這里本該有readsp3b函數(shù)截圖,但是太長我就不放了)
readsp3b函數(shù)讀取文件體,EOF是文件結(jié)尾字符,讀取即跳出。*開頭的一行都是時間,因此調(diào)用str2time轉(zhuǎn)換為軟件內(nèi)通用的時間制式。141-152行為每個時間的星歷生成存儲空間并初始化為0。每個時間的星歷前4為表示編號,第一位為P表示后面的是位置和鐘差,據(jù)說,若第一位為V則表示后面的是速度和鐘差變率,不過V開頭的精密星歷我還沒下載到過。157-162行是把帶字母的編號轉(zhuǎn)換成數(shù)字順序,畢竟數(shù)組調(diào)用時不認字母。接下來便是對位置和鐘差信息驗證和賦值,在199行調(diào)用addpeph函數(shù)把讀取完的一個時刻的精密星歷結(jié)構(gòu)體peph賦給nav結(jié)構(gòu)體。
在execses_b中,readpreceph接下來是讀取基站名稱之類的操作,個人感覺不是很重要,然后是讀取擴展路徑,進入execses_r函數(shù),r表示rover,移動站,與之前的b(base,基站)區(qū)分,然后也是讀取流動站名稱之類的操作,然后進入execses函數(shù),這個函數(shù)囊括了讀寫和運算的順序,本文在第四節(jié)將解釋其中的setpcv部分。在觀測數(shù)據(jù)、廣播星歷、天線相位修正、海洋潮汐負荷讀取結(jié)束,調(diào)試日志(.trace)和定位結(jié)果文件(.pos)文件頭寫入之后,進入處理模式選擇,包括前向、后向、聯(lián)合三種,和rtkpost.exe在options中第一頁的選擇一致。最終都會進入procpos函數(shù)進行處理。在這個函數(shù)中完成模式和衛(wèi)星系統(tǒng)的一些選擇,又指向rtkpos函數(shù)(文件rtkpos.c中)進行處理。在這里,根據(jù)之前設(shè)置的選擇配置,開始進入不同的定位函數(shù),相對定位的、單點定位的等等,這里假設(shè)進入pppos函數(shù)(ppp.c文件)。衛(wèi)星位置和鐘差計算在satposs函數(shù)中(ephemeris.c)。
satposs函數(shù)里,先讀取了當前接收機觀測時刻里所有的偽距信息,并由此除以光速,加上廣播星歷中的鐘差改正,計算出了衛(wèi)星信號播發(fā)的時刻。由于有精密星歷,進入satpos函數(shù)后,在swtich選擇中會進入peph2pos函數(shù),這就是計算給定時刻(即衛(wèi)星信號播發(fā)時刻)的由精密星歷計算的衛(wèi)星位置和鐘差的函數(shù)。
三、精密星歷下衛(wèi)星位置和鐘差的計算
在peph2pos函數(shù)(preceph.c文件)中,相關(guān)參數(shù)的計算流程如下:先計算衛(wèi)星信號播發(fā)時刻的位置和鐘差(調(diào)用了pephpos,如有clk文件還會調(diào)用pephclk文件),然后計算0.001秒后的位置和鐘差,兩個位置作差除以時間間隔即是速度(求導麻煩,編程中通常都靠自變量增加一個小量然后作差除以小量來近似計算導數(shù)),然后進行衛(wèi)星天線相位中心的計算(精密星歷確定的是衛(wèi)星幾何中心的位置,但是信號播發(fā)是天線完成的,這中間有一個桿臂矢量需要修正),如有精密鐘差產(chǎn)品還會進行相對論效應(yīng)的修正。
接下來學習pephpos函數(shù),satantoff將在第四節(jié)敘述。

pephpos函數(shù)的一部分
preceph.c文件405-410行,是對比了衛(wèi)星信號播發(fā)時刻是否在讀取的精密星歷起止時間范圍內(nèi)。因為rtklib里對星歷的插值是用的多項式插值法,數(shù)值分析中我們知道,高次插值會發(fā)生龍格現(xiàn)象,延伸段的結(jié)果可信度較低。因此如果給定的信號播發(fā)時刻距離精密星歷起止太遠,會顯示沒有精密星歷可用。
412-415行用的是二分法搜索,搜索距離衛(wèi)星信號播發(fā)時刻最近的星歷數(shù)組指針。三行代碼就凝練地完成了二分法搜索功能,再次告訴我我就是個辣雞,還是不可回收的那種。
419-428行算是在校驗參數(shù),因為讀取的是結(jié)構(gòu)體中的數(shù)據(jù),沒有讀取到的或者無效的部分都賦值為0了,,因此這一部分的觀測值時沒有精密星歷數(shù)據(jù)的,需要跳過。

pephpos函數(shù)的一部分
429-442行進行了地球自轉(zhuǎn)改正,衛(wèi)星位置是使用的ECEF系進行表示,因此z軸無需改正。
443-445行是調(diào)用了interppol函數(shù)進行多項式插值計算,使用了Neville算法,默認是11個節(jié)點進行10階多項式插值。其實我也不是很懂為什么要用這么高的次數(shù),但是插值精度確實很高,尤其是滿足前后超過5個節(jié)點時。
446-454行是對位置的計算結(jié)果進行誤差估計。
455-459行是對鐘差線性插值參數(shù)的計算,在472行對信號播發(fā)時刻的鐘差進行了線性插值。
pephclk函數(shù)對于鐘差的計算也是如此,但是clk文件鐘差給的比較頻繁,比如30s左右就給一個,比起sp3文件5分鐘、15分鐘給一個要頻繁得多,而鐘差通常變化比較慢,因此另外的clk文件確實能夠提供更高的鐘差插值精度。
四、衛(wèi)星天線相位中心修正

rtklib手冊對于衛(wèi)星天線相位修正和體坐標系的示意圖
衛(wèi)星天線相位修正通常是三步:一是讀取atx文件,二是計算天線相位中心偏移(pco)的修正量,三是根據(jù)地面接收機的位置計算相對于衛(wèi)星天線相位中心的方位,進行天線相位變化(pcv)的修正。由于atx文件中沒有北斗衛(wèi)星的pcv信息(都是0),因此本節(jié)只作pco修正的敘述,敬請諒解。
衛(wèi)星天線相位中心的讀取和設(shè)置的調(diào)用路徑為:main里的postpos->openses->readpcv,atx文件到現(xiàn)在十幾兆,很多并不是衛(wèi)星的相位中心修正,而是接收機的修正。接收機的修正也是調(diào)用了這個readpcv函數(shù)。對于atx文件,readpcv里通過readantex函數(shù)(rtkcmn.c文件中)進行讀取,存儲在pcvs結(jié)構(gòu)體中,這個結(jié)構(gòu)體最終指向的也是一個靜態(tài)全局結(jié)構(gòu)體pcvss,相應(yīng)的,接收機也有個pcvsr結(jié)構(gòu)體存儲信息。
在atx文件中,第60個字符后有該行的說明(如果該行是一大串數(shù)據(jù)就沒有)。START OF ANTENNA表示開始顯示一次新的天線相位中心修正信息,到END OF? ANTENNA。由于衛(wèi)星存在更替、相位測量結(jié)果存在變化,因此同一顆衛(wèi)星可能有好幾個相位中心修正信息,這個由有效期來控制,VALID FROM和VALID UNTIL分別給出了相位修正的有效起止時間。各個頻點的相位中心可能不一致,由START OF FREQUENCY和END OF FREQUENCY控制。這里也有個我比較奇怪的地方,就是下面截圖的C01星,是北斗二代的地球靜止衛(wèi)星,播發(fā)三個頻點的信號(B1I/B2I/B3I),怎么會有4個頻點信息呢?如有大神路過還請不吝指教。NORTH/EAST/UP是頻點的天線相位中心偏移,NOAZI是天線相位中心變化,隔一定角度給出一個值,這個角度由之前的ZEN1/ZEN2/DZEN計算給定。

igs14_2118.atx文件部分內(nèi)容截圖
readantex函數(shù)的源代碼就不全放了,主要是根據(jù)atx文件60位之后的說明信息不同,將各種信息存儲至pcv結(jié)構(gòu)體中,還是比較好懂的。比較搞不懂的一點在2089-2092行,見下圖。一顆衛(wèi)星的載波信號總數(shù)NFREQ為3,但是通過sscanf函數(shù)掃描近來的f只比較了前三個freqs(也即1,2,5號頻點),那么北斗的頻點在2,7,6號,就只能讀取一個2號頻點的(還有個不知怎么就有的1號頻點的),那么這樣的修正值顯然是不對的,這讓我百思不得其解,還請指教。

readantex函數(shù)部分截圖
另外,訪問【http://www.csno-tarc.cn/support/satelliteparameters】可以獲得更加準確的衛(wèi)星天線相位中心的偏移信息,不知道為什么atx文件迄今仍不采納。插句題外話,據(jù)說國外的北斗使用效果不如GPS的原因,一方面是衛(wèi)星更少,另一方面就國外設(shè)備不上心,從硬件制造時信號捕獲的能力不足,到軟件編寫的方面。
衛(wèi)星pcv信息在以年為單位時時頻繁變化的,但在實驗時就那么幾天甚至幾分鐘,在這段時間通常是沒有變化的。因此給定一個時間后就可以設(shè)置當前需要的pcv信息,并讀取到navs結(jié)構(gòu)體中。這里用到的函數(shù)就是setpcv,調(diào)用路徑為main里的postpos->execses_b->execses_r->execses->setpcv。這個函數(shù)在postpos.c文件中,其中第826-834行是衛(wèi)星天線pcv信息的設(shè)置。調(diào)用searchpcv函數(shù)讀取pcv,然后放置進nav結(jié)構(gòu)體中。
對于searchpcv函數(shù),也是比較容易讀懂的,不放源代碼截圖了。主要就是通過衛(wèi)星prn號和修正信息的起止時間,逐個比較,確定需要pcv的信息并返回。
最后是對于讀取的信息的使用,也即天線相位中心偏移的修正。讀取atx文件所得的數(shù)據(jù)是衛(wèi)星體坐標系下的表示,為了計算需要必須投影至ECEF系下表示,這就用到了satantoff函數(shù)(preceph.c文件中),它在獲得了精密星歷的位置和鐘差結(jié)果之后進行修正,調(diào)用順序在第三節(jié)中有。

satantoff函數(shù)截圖
初看這個函數(shù)我是懵逼的,滿臉這啥玩意兒啊咋回事啊,咋還和太陽位置扯上關(guān)系了呢。回到上面rtklib手冊對pco信息的示意圖,下面就是坐標系的解釋,衛(wèi)星發(fā)送信號的天線始終對地定向,據(jù)說通過調(diào)姿始終保持對地心±1°的偏差內(nèi),而同時,大功率信號發(fā)送需要足夠的電能,因此繞衛(wèi)星質(zhì)心到地心的連線旋轉(zhuǎn),需使太陽能帆板盡量正對太陽。衛(wèi)星對于活動部件其實是極為抗拒的,尤其是這種要服役十幾年的衛(wèi)星,太陽能帆板是不會活動的,只能依靠調(diào)姿進行,因此我們就可以假設(shè),衛(wèi)星體坐標系一軸指向地心,一軸指向太陽。這就是sunmoonpos函數(shù)的作用了,給出了ECEF系下衛(wèi)星到太陽的方向矢量,和衛(wèi)星到地心的矢量叉乘,即得到了體坐標系在ECEF下的表示,最后加上電離層修正(這一塊原理我不大懂),就得到了衛(wèi)星天線相位中心偏移量在ECEF下的表示,用dant表示。
接下來康康sunmoonpos函數(shù),時間要轉(zhuǎn)換為utc(感覺過不久就應(yīng)該調(diào)整差值為19秒了呢),再轉(zhuǎn)為ut1,然后計算太陽和月亮方向在地心慣性系下的表示,調(diào)用了sunmoonpos_eci函數(shù),這一步觸及了知識盲區(qū)。。各種計算我一個都看不懂。跳過這一步往下看看,是計算地心慣性系(ECI)到地心地固系(ECEF)的轉(zhuǎn)換函數(shù),然后坐標系轉(zhuǎn)換矩陣乘矢量,就得到了新坐標系下的矢量表示。

從頭懵到尾的eci2ecef函數(shù)
那么是如何轉(zhuǎn)換的呢?打開eci2ecef函數(shù)看一下。。先是世界時、GPS時間系統(tǒng)和格林尼治平恒星時的混亂關(guān)系。。然后計算了個不知道干啥用的天文參數(shù)(不負責任的大膽猜測可能是極移的參數(shù)),然后是使用了iau1976歲差模型,又使用了iau1980章動模型,計算了格林尼治真恒星時,最后就一通操作得到了該時刻ECI到ECEF系的轉(zhuǎn)換矩陣和格林尼治平恒星時(假裝看懂了)。至于中間的歲差、章動、極移(猜測)的函數(shù),有興趣可以繼續(xù)看看,我是已經(jīng)看不懂了。
好了,以上就是這篇學習經(jīng)歷的全部內(nèi)容了。不得不承認自己的水平實在有限,文中謬誤之處還請各位看官指正。rtklib作為開源軟件到現(xiàn)在已經(jīng)完善了十幾年了,感謝全世界參與其中工作的學者們,身為后輩希望能夠盡早追上前人的腳步,在學術(shù)上做一點微不足道的貢獻。