寫在前面:研究方向從精密單點定位改成了載波相位差分定位和基線解算,因此不再使用GAMP,但是GAMP是基於rtklib二次開發的,因此代碼結構以及函數沒有很大的差別。
首先給出我覺得寫的比較好的博客:
一、主要類(結構體)的解釋
首先解釋一下各種主要的類,有助於理解接下來的函數作用以及調用關系
1.1 處理選項類prcopt_t
- 類名:
prcopt_t - 所在文件:rtklib.h
- 功能說明:記錄各種處理選項,按照用戶需求進行初始化定義
mode:定位模式選項,單點定位,差分定位,PPP等soltype:輸出結果的形式,有三種模式(0:forward,1:backward,2:combined),具體有什么區別暫時不清楚nf:頻數選項,(1:L1,2:L1+L2,3:L1+L2+L5)navsys:處理的導航系統elmin:截止高度角snrmask:信噪比標記sateph:星歷選項,精密星歷,廣播星歷等- 還有很多,不一一列出來,源代碼中有注釋,有的變量我也不是很清楚在哪使用
1.2 衛星文件類filopt_t
- 類名:
filopt_t - 所在文件:rtklib.h
- 功能說明:儲存各種衛星文件
satantp:衛星天線參數文件rcvantp:接收機天線參數文件stapos:站點位置文件geoid:大地水准面數據iono:電離層數據dcb:DCBeop:地球定向參數blq:潮汐tempdir:ftp/http 網址臨時目錄geexe:用於運行谷歌地球?搞不懂solstat:結果目錄trace:正在調試的文件,處理數據的過程中經常會調用trace函數,用來標記程序當前的運行位置
1.3 輸出結果選項類solopt_t
- 類名:
solopt_t - 所在文件:rtklib.h
- 功能說明:定義輸出的成果文件的格式
1.4 衛星狀態類ssat_t
- 類名:
ssat_t - 所在文件:rtklib.h
- 功能說明:衛星當前狀態參數,該衛星所屬的導航系統,坐標,高度角等等
sys:導航系統vs:有效衛星單一標志,不知道什么東西azel:方位角,高度角resp:偽距殘差resc:載波相位殘差vsat:有效衛星標志,不知道什么東西snr:信噪比,描述信號強度fix:模糊度的狀態,浮點解、固定解slip:周跳- 其他的以后再補上,這個類儲存的信息還是挺重要的
1.5 解算結果類sol_t
- 類名:
sol_t - 所在文件:rtklib.h
- 功能說明:儲存計算成果,定位坐標就存在這個類里面,還有接收機鍾差,解決狀態(未解決,浮點解,PPP解,差分解,單點解等)
1.6 RTK進程類rtk_t
- 類名:
rtk_t - 所在文件:rtklib.h
- 功能說明:包含上述的解算結果類
sol_t以及衛星狀態類ssat_t,在postpos.c文件中的propos函數中產生,由rtkinit函數初始化,再調入rtkpos函數中參與坐標解算,根據解算進程改變儲存的衛星狀態,並儲存最終的解算坐標
二、解算流程
根據上述對主要的幾個類的解釋,不難想到,處理選項類prcopt_t、輸出結果選項類solopt_t、衛星文件類filopt_t需要用戶在main函數中先根據用戶的需求進行定義,接下來就可以根據這幾個類進行解算了。
rtklib源碼中給出的示例的main函數是用戶在VS的命令行中寫入配置文件,然后用一段程序去讀取這個命令行,對上述的三個類進行賦值。有需求的可以自己更改,主要就是定義好這三個類,就可以進行解算,筆者會在下一次日志中專門給出一個示例
2.1 main函數
- 函數名:
main - 所在文件:main.c(rnx2rtkp.c)
- 功能說明:三個基本類的定義
prcopt_t、solopt_t、filopt_t,並將其傳入定位解算進程函數postpos
ret=postpos(ts,te,tint,0.0,&prcopt,&solopt,&filopt,infile,n,outfile,"","");
2.2 postpos函數
- 函數名:
postpos - 所在文件:postpos.c
- 功能說明:函數openses讀取天線參數文件以及大地水准面數據,接着根據處理時間進行判斷,具體為什么要做這個判斷以及不同判斷下對文件做了什么樣的操作暫時不清楚,有待后續解決,猜想可能是因為一般衛星文件名都會帶有一些信息,比如歷元,需要從這些文件名中獲取一些信息,然后再進行操作。最終都是要調用函數execses_b
/* execute processing session */
stat=execses_b(ts,te,ti,popt,sopt,fopt,1,ifile,index,n,ofile,rov,base);
ts:進程開始時間te:進程結束時間ti:處理間隔tu:處理單位時間popt:不解釋sopt:不解釋fopt:不解釋infile:輸入文件n:輸入文件數量outfile:輸出文件rov:流動站IDbase:基准站ID
2.3 execses_b函數
- 函數名:
execses_b - 所在文件:postpos.c
- 功能說明:通過函數readpreceph讀取精密星歷,然后為每個基准站執行操作,這個函數筆者暫時也不是了解的非常清楚,但是大致可以推斷出的是,用戶傳入的基准站可能不止一個,所以在執行處理流動站的函數execses_r之前進行了一些字符串操作,應該是在分割傳入的基准站的字符,然后循環處理每個基准站。最后都是要進入函數execses_r,處理每個流動站
stat=execses_r(ts,te,ti,popt,sopt,fopt,flag,infile,index,n,outfile,rov);
2.4 execses_r函數
- 函數名:
execses_r - 所在文件:postpos.c
- 功能說明:同樣的對於流動站的操作首先也需要分割字符串,識別每個流動站,這部分代碼與函數execses_b基本一致,最后進入函數execses,處理每一個流動站
/* execute processing session */
stat=execses(ts,te,ti,popt,sopt,fopt,flag,ifile,index,n,ofile);
2.5 execses函數
- 函數名:
execses - 所在文件:postpos.c
- 功能說明:讀取電離層參數
readtec、erp文件readerp(地球章動、極移參數)、讀取觀測值文件以及導航電文readobsnav、差分碼文件readdcb、設置天線參數setpcv(天線參數在函數postpos中已經讀取過)、地球潮汐改正文件readotl;之后開始根據開始定義的prcopt_t變量中儲存的定位模式獲取天線相位中心位置antpos;接着寫輸出的結果文件的頭文件outhead;然后根據prcopt_t變量中儲存的soltype變量即上文提到的輸出結果的形式,對一些索引變量進行賦值操作iobsu=iobsr=isbs=revs=aborts=0;,revs=1; iobsu=iobsr=obss.n-1; isbs=sbss.n-1;等,這里不同的soltype會導致輸出有什么不同暫時不清楚,待后續查明,最后進入函數procpos進行下一步處理
//選擇不同的soltype,索引參數會變化
iobsu=iobsr=isbs=revs=aborts=0;
/* forward */
procpos(fp,&popt_,sopt,0);
/* backward */
revs=1; iobsu=iobsr=obss.n-1; isbs=sbss.n-1;
procpos(fp,&popt_,sopt,0);
/* combined */
isolf=isolb=0;
procpos(NULL,&popt_,sopt,1);
revs=1; iobsu=iobsr=obss.n-1; isbs=sbss.n-1;
procpos(NULL,&popt_,sopt,1);
2.6 procpos函數
- 函數名:
procpos - 所在文件:postpos.c
- 功能說明:這個函數開始產生變量
rtk_t,obsd_t,sol_t,開始正式進入rtk解算流程,首先通過函數rtkinit初始化變量rtk_t;然后進入while循環,調用函數inputobs,這個函數的作用應該是將之前讀入的觀測值文件、導航電文等的信息傳入rtk_t、obsd_t這兩個變量之中,obsd_t記錄了衛星號接收機號、偽距、載波、信噪比等等觀測值信息,這些都是解算坐標最重要的信息;固定好導航系統后satsys,下一步進行了載波相位偏差糾正corr_phase_bias_ssr,使用的SSR方法,具體這個方法不清楚,但是這里存在一個疑惑,筆者整體看了代碼,發現解算差分定位之前並沒有進行周跳修復檢測,對於載波觀測值只有這里的偏差糾正,而PPP中是有的,可能是沒有找到,有待后續解決;然后進入函數rtkpos開始解算坐標;本函數后面的代碼用於輸出結果文件,同樣對於不同的soltype有不同的輸出方式
/* carrier-phase bias correction */
//載波相位偏差校正
if (!strstr(popt->pppopt,"-ENA_FCB")) {
corr_phase_bias_ssr(obs,n,&navs);
}
//定位
if (!rtkpos(&rtk,obs,n,&navs)) continue;
//輸出結果文件
if (mode==0) { /* forward/backward */
if (!solstatic) {
outsol(fp,&rtk.sol,rtk.rb,sopt);
}
else if (time.time==0||pri[rtk.sol.stat]<=pri[sol.stat]) {
sol=rtk.sol;
for (i=0;i<3;i++) rb[i]=rtk.rb[i];
if (time.time==0||timediff(rtk.sol.time,time)<0.0) {
time=rtk.sol.time;
}
}
}
else if (!revs) { /* combined-forward */
if (isolf>=nepoch) return;
solf[isolf]=rtk.sol;
for (i=0;i<3;i++) rbf[i+isolf*3]=rtk.rb[i];
isolf++;
}
else { /* combined-backward */
if (isolb>=nepoch) return;
solb[isolb]=rtk.sol;
for (i=0;i<3;i++) rbb[i+isolb*3]=rtk.rb[i];
isolb++;
}
2.7 rtkpos函數
- 函數名:
rtkpos - 所在文件:postpos.c
- 功能說明:結合代碼說明,見如下代碼
//1、這里定義了一個prcopt_t用來儲存傳入的rtk_t中的prcopt_t
prcopt_t *opt=&rtk->opt;
//2、設置基准站位置,滿足如下的if條件后(應該是單點定位不需要基准站,動基線不太懂),坐標不變
//,每個方向的變化速率變成0
/* set base staion position */
if (opt->refpos<=POSOPT_RINEX&&opt->mode!=PMODE_SINGLE&&
opt->mode!=PMODE_MOVEB) {
for (i=0;i<6;i++) rtk->rb[i]=i<3?opt->rb[i]:0.0;
}
//3、計算流動站/基准站觀測值數量,可用於后面判斷是否滿足差分條件
/* count rover/base station observations */
for (nu=0;nu <n&&obs[nu ].rcv==1;nu++) ;
for (nr=0;nu+nr<n&&obs[nu+nr].rcv==2;nr++) ;
//4、單點定位解算流動站坐標pntpos,單點定位解算的坐標可以作為初始值參與其他精密定位方法
/* rover position by single point positioning */
if (!pntpos(obs,nu,nav,&rtk->opt,&rtk->sol,NULL,rtk->ssat,msg)) {
errmsg(rtk,"point pos error (%s)\n",msg);
if (!rtk->opt.dynamics) {
outsolstat(rtk);
return 0;
}
}
//5、如果定位模式是單點定位,那么就可以直接輸出結果結束了
/* single point positioning */
if (opt->mode==PMODE_SINGLE) {
outsolstat(rtk);
return 1;
}
//6、精密單點定位
/* precise point positioning */
if (opt->mode>=PMODE_PPP_KINEMA) {
pppos(rtk,obs,nu,nav);
outsolstat(rtk);
return 1;
}
//7、接下來是需要基准站的差分定位和基線解算,因此首先需要檢查基站數據數量和差分時間
/* check number of data of base station and age of differential */
if (nr==0) {
errmsg(rtk,"no base station observation data for rtk\n");
outsolstat(rtk);
return 1;
}
//8、動基線與其他差分定位方式,動基線的基站坐標需要隨時間同步變化,所以需要計算出變化速率
//,解釋了為什么第二步除了單點定位,動基線也不參與基站解算,動基線在這里單獨解算
if (opt->mode==PMODE_MOVEB) { /* moving baseline */
/* estimate position/velocity of base station */
if (!pntpos(obs+nu,nr,nav,&rtk->opt,&solb,NULL,NULL,msg)) {
errmsg(rtk,"base station position error (%s)\n",msg);
return 0;
}
rtk->sol.age=(float)timediff(rtk->sol.time,solb.time);
if (fabs(rtk->sol.age)>TTOL_MOVEB) {
errmsg(rtk,"time sync error for moving-base (age=%.1f)\n",rtk->sol.age);
return 0;
}
for (i=0;i<6;i++) rtk->rb[i]=solb.rr[i];
/* time-synchronized position of base station */
for (i=0;i<3;i++) rtk->rb[i]+=rtk->rb[i+3]*rtk->sol.age;
}
else {
rtk->sol.age=(float)timediff(obs[0].time,obs[nu].time);
if (fabs(rtk->sol.age)>opt->maxtdiff) {
errmsg(rtk,"age of differential error (age=%.1f)\n",rtk->sol.age);
outsolstat(rtk);
return 1;
}
}
//9、上面的步驟只算了相對定位的差分時間和動基線坐標,並沒有進行坐標解算
//,這里進行相位定位,並輸出最終結果,到這里定位步驟全部完成
/* relative potitioning */
relpos(rtk,obs,nu,nr,nav);
outsolstat(rtk);
