RTKLIB单点定位处理流程之一（postpos/后处理）

绝地灬酷狼 2023-03-12 10:53 21阅读 0赞

> ### 版本： ###
> 
> ### \#define VER\_RTKLIB "2.4.3" /\* library version \*/ ###
> 
> ### \#define PATCH\_LEVEL "b29" /\* patch level \*/ ###

# **处理参数设置：****主要是三个结构体：** #

> prcopt\_t opt = prcopt\_default; //processing options
> 
> solopt\_t sopt = solopt\_default; //solution options
> 
> filopt\_t fopt = \{ 0 \}; //file options

#include "rtklib.h"
    int main() {
    	char* infile[] = { {"移动站观测数据"},
    			    {"移动站星历文件"} };
    	char* ofile = { "输出结果文件/myposbyVS.pos" };
    	gtime_t ts = { 0 }, te = { 0 };
    	double st[] = { 2019,12,5,0,0,0 }; //开始处理时间
    	double ed[] = { 2019,12,5,1,0,0 }; //结束处理时间
    	ts = epoch2time(st);
    	te = epoch2time(ed);
    	double ti = 0.0;      //  processing interval(s) (0:all) 时间间隔
    	double tu = 0.0;    //   processing unit time(s) (0:all)
    	prcopt_t opt = prcopt_default;  //processing options
    	solopt_t sopt = solopt_default; //solution options
    	filopt_t fopt = { 0 }; //file options
    	int n = 2;  //number of input files
            //GPS单点定位参数设置
    	opt.navsys = SYS_GPS;
    	opt.mode = PMODE_SINGLE; // PMODE_SINGLE; // PMODE_MOVEB; // ; PMODE_FIXED           
    	opt.nf = 1;             /* number of frequencies (1:L1,2:L1+L2,3:L1+L2+L5) */
    	opt.elmin = 10 * D2R;
    	opt.ionoopt = IONOOPT_BRDC;       /* ionosphere option (IONOOPT_???) */
    	opt.tropopt = TROPOPT_SAAS;       /* troposphere option (TROPOPT_???) */
    	opt.dynamics = 0;       /* dynamics model (0:none,1:velociy,2:accel) */
    	opt.tidecorr = 0;       /* earth tide correction (0:off,1:solid,2:solid+otl+pole) */
    	opt.niter = 1;          /* number of filter iteration */
    	opt.codesmooth = 0;     /* code smoothing window size (0:none) */
    	opt.rovpos = 0;         /* rover position for fixed mode */
    	opt.refpos = 0;         /* base position for relative mode */
    						/* (0:pos in prcopt,  1:average of single pos, */
    						/*  2:read from file, 3:rinex header, 4:rtcm pos) */
    	//opt.eratio[] = {100.0,100.0}; /* code/phase error ratio */
    	//opt.err[] = { 100.0,0.003,0.003,0.0,1.0 }; /* err[] */
    	/* measurement error factor */
    	/* [0]:reserved */
    	/* [1-3]:error factor a/b/c of phase (m) */
    	/* [4]:doppler frequency (hz) */
    	//opt.std[] = { 30.0,0.03,0.3 };            /* std[] */
    	/* initial-state std [0]bias,[1]iono [2]trop */
    	//opt.prn[] = { 1E-4,1E-3,1E-4,1E-1,1E-2 ,0.0 };
    	/* process-noise std [0]bias,[1]iono [2]trop [3]acch [4]accv [5] pos */
    
    //求解选项定义
    	sopt.posf = SOLF_LLH; // SOLF_XYZ; //  SOLF_ENU;
    	sopt.times = TIMES_GPST; /* time system: gps time */
    	sopt.timef = 1; /* time format (0:sssss.s,1:yyyy/mm/dd hh:mm:ss.s) */
    	sopt.timeu = 3; /* time digits under decimal point */
    	sopt.degf = 0; /* latitude/longitude format (0:ddd.ddd,1:ddd mm ss) */
    	sopt.outhead = 1;        /* output header (0:no,1:yes) */
    	sopt.outopt = 1;         /* output processing options (0:no,1:yes) */
    	sopt.outvel = 1;         /* output velocity options (0:no,1:yes) */
    	sopt.datum = 0;          /* datum (0:WGS84,1:Tokyo) */
    	sopt.height = 0;         /* height (0:ellipsoidal,1:geodetic) */
    	sopt.geoid = 0;          /* geoid model (0:EGM96,1:JGD2000) */
    	sopt.solstatic = 0;      /* solution of static mode (0:all,1:single) */
    	sopt.sstat = 1;          /* solution statistics level (0:off,1:states,2:residuals) */
    	sopt.trace = 3;          /* debug trace level (0:off,1-5:debug) */
    	//sopt.nmeaintv[] = { 0.0,0.0 }; /* nmea output interval (s) (<0:no,0:all) */
    					/* nmeaintv[0]:gprmc,gpgga,nmeaintv[1]:gpgsv */
    	//sopt.sep[] = " ";       /* field separator */
    	//sopt.prog[] = "";      /* program name */
    	//sopt.maxsolstd = 0;   /* max std-dev for solution output (m) (0:all) */
    
    //文件设置选项 /* file options type */
    	// #define MAXSTRPATH  1024        /* max length of stream path */
    	//fopt.satantp[MAXSTRPATH] = { 0 }; /* satellite antenna parameters file */
    	//fopt.rcvantp[MAXSTRPATH] = { 0 }; /* receiver antenna parameters file */
    	//fopt.stapos[MAXSTRPATH] = { 0 }; /* station positions file */
    	//fopt.geoid[MAXSTRPATH] = { 0 }; /* external geoid data file */
    	//fopt.iono[MAXSTRPATH] = { 0 }; /* ionosphere data file */
    	//fopt.dcb[MAXSTRPATH] = { 0 }; /* dcb data file */
    	//fopt.eop[MAXSTRPATH] = { 0 }; /* eop data file */
    	//fopt.blq[MAXSTRPATH] = { 0 }; /* ocean tide loading blq file */
    	//fopt.tempdir[MAXSTRPATH] = { 0 }; /* ftp/http temporaly directory */
    	//fopt.geexe[MAXSTRPATH] = { 0 }; /* google earth exec file */
    	//fopt.solstat[MAXSTRPATH] = { 0 }; /* solution statistics file */
    	//fopt.trace[MAXSTRPATH] = { 0 }; /* debug trace file */
    
    
    	postpos(ts,te,ti,tu,&opt,&sopt,&fopt,infile,n,ofile,"","");
    	printf("*********postpos over*********\n");
    	getchar();
    	return 0;
    }

函数调用关系：

首先通过函数：postpos(ts,te,ti,tu,&opt,&sopt,&fopt,infile,n,ofile,"","");进入子函数；

*  # openses函数：预处理 #

目的：通过读取filopt\_t fopt中的参数，其中filopt\_t结构体定义如下，实现对‘文件的读取’，读取子函数为openses！

typedef struct {        /* file options type */
        char satantp[MAXSTRPATH]; /* satellite antenna parameters file */
        char rcvantp[MAXSTRPATH]; /* receiver antenna parameters file */
        char stapos [MAXSTRPATH]; /* station positions file */
        char geoid  [MAXSTRPATH]; /* external geoid data file */
        char iono   [MAXSTRPATH]; /* ionosphere data file */
        char dcb    [MAXSTRPATH]; /* dcb data file */
        char eop    [MAXSTRPATH]; /* eop data file */
        char blq    [MAXSTRPATH]; /* ocean tide loading blq file */
        char tempdir[MAXSTRPATH]; /* ftp/http temporaly directory */
        char geexe  [MAXSTRPATH]; /* google earth exec file */
        char solstat[MAXSTRPATH]; /* solution statistics file */
        char trace  [MAXSTRPATH]; /* debug trace file */
    } filopt_t;

/* open procssing session ----------------------------------------------------*/
    static int openses(const prcopt_t *popt, const solopt_t *sopt,
                       const filopt_t *fopt, nav_t *nav, pcvs_t *pcvs, pcvs_t *pcvr)

然后为：判断处理时间是否正确，否则程序终止；

然后为：通过malloc函数为ifile\[\]，动态分配内存，否则程序终止；

然后为：进行文件的读取，对于sp3/SP3、.eph/EPH文件的时间，起使和终止相差1小时；brdc文件的开始和终止时间差2小时。

然后为：路径替换，实现语句：

nf+=reppaths(infile[j],ifile+nf,MAXINFILE-nf,tts,ttte,"","");
    子函数内实现语句：
    reppath(path,rpath[n],time,rov,base);

经过两次读取循环，已经将观测、导航文件读取到ifie数组中：

![20200512094440399.png][]

将输出结果文件路径保存到ofile数组中：

if (!reppath(outfile,ofile,tts,"","")&&i>0) flag=0;
    子函数实现代码：
    strcpy(rpath,path);

*  # execses\_b #

> 在main中，输入文件中：第一个为观测文件，第二个为导航文件；

## 1、readpreceph ##

/\* 读取文件为：read prec ephemeris, sbas data, lex data, tec grid and open rtcm \----------\*/  
static void readpreceph(char \*\*infile, int n, const prcopt\_t \*prcopt,nav\_t \*nav, sbs\_t \*sbs, lex\_t \*lex)

> 读取文件的函数依次为：
> 
> /\* read precise ephemeris files \*/
> 
> /\* read precise clock files \*/
> 
> /\* read satellite fcb files \*/
> 
> /\* read solution status files for ppp correction \*/
> 
> /\* read sbas message files \*/
> 
> /\* read lex message files \*/
> 
> /\* allocate sbas ephemeris \*/
> 
> /\* set rtcm file and initialize rtcm struct \*/

以函数：readsp3为例

/* read precise ephemeris files */
        for (i=0;i<n;i++) {
            if (strstr(infile[i],"%r")||strstr(infile[i],"%b")) continue;
            readsp3(infile[i],nav,0);
        }

第一次循环读入的是观测文件，因此进入函数之后，动态分配内存、释放内存，结束！

第二次循环读入的是导航文件，因此进入函数之后，动态分配内存、释放内存，结束！

在函数readrnxc中：读取钟差文件时，读取导航文件时，将导航文件，文件头的相关参数保存到nav中：

/* read precise clock files */
        for (i=0;i<n;i++) {
            if (strstr(infile[i],"%r")||strstr(infile[i],"%b")) continue;
            readrnxc(infile[i],nav);
        }

导航文件，文件头为：

![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3d1d3VrdTEyMw_size_16_color_FFFFFF_t_70][]

nav:

![20200512105333577.png][]

![20200512105402228.png][]

# 2、execses\_r #

stat=execses\_r(ts,te,ti,popt,sopt,fopt,flag,infile,index,n,outfile,rov);

## 2.1 execses ##

stat=execses(ts,te,ti,popt,sopt,fopt,flag,infile,index,n,outfile);

如果main函数中，设置了调试等级：

sopt.sstat = 1;   /* solution statistics level (0:off,1:states,2:residuals) */
    	sopt.trace = 3;   /* debug trace level (0:off,1-5:debug) */

则在此处新建生成的调试数据，保存的路径、文件名

/* open debug trace */
        if (flag&&sopt->trace>0) {
            if (*outfile) {
                strcpy(tracefile,outfile);
                strcat(tracefile,".trace"); //输出文件名，后加.trace
            }
            else {
                strcpy(tracefile,fopt->trace);
            }
            traceclose(); //关闭记录文件
            traceopen(tracefile); //开始对 trace文件进行写处理
            tracelevel(sopt->trace); //读取调试等级
        }

顺序执行，子函数如下：

1、/\* read ionosphere data file \*/ 读取数据放到如下结构体中：

typedef struct {        /* TEC grid type */
        gtime_t time;       /* epoch time (GPST) */
        int ndata[3];       /* TEC grid data size {nlat,nlon,nhgt} */
        double rb;          /* earth radius (km) */
        double lats[3];     /* latitude start/interval (deg) */
        double lons[3];     /* longitude start/interval (deg) */
        double hgts[3];     /* heights start/interval (km) */
        double *data;       /* TEC grid data (tecu) */
        float *rms;         /* RMS values (tecu) */
    } tec_t;

2、/\* read erp data \*/ 读取数据放到如下结构体中：

typedef struct {        /* earth rotation parameter type */
        int n,nmax;         /* number and max number of data */
        erpd_t *data;       /* earth rotation parameter data */
    } erp_t;

3、/\* read obs and nav data \*/

其中观测数据结构体的使用：

typedef struct {        /* observation data record */
        gtime_t time;       /* receiver sampling time (GPST) */
        unsigned char sat,rcv; /* satellite/receiver number */
        unsigned char SNR [NFREQ+NEXOBS]; /* signal strength (0.25 dBHz) */
        unsigned char LLI [NFREQ+NEXOBS]; /* loss of lock indicator */
        unsigned char code[NFREQ+NEXOBS]; /* code indicator (CODE_???) */
        double L[NFREQ+NEXOBS]; /* observation data carrier-phase (cycle) */
        double P[NFREQ+NEXOBS]; /* observation data pseudorange (m) */
        float  D[NFREQ+NEXOBS]; /* observation data doppler frequency (Hz) */
    } obsd_t;
    
    typedef struct {        /* observation data */
        int n,nmax;         /* number of obervation data/allocated */
        obsd_t *data;       /* observation data records */
    } obs_t;

### 2.1.1 readobsnav ###

函数： static int readobsnav(gtime\_t ts, gtime\_t te, double ti, char \*\*infile,const int \*index, int n, const prcopt\_t \*prcopt,obs\_t \*obs, nav\_t \*nav, sta\_t \*sta)

程序中，每次都是将相关的变量置0、指针变量赋值为null

obs->data=NULL; obs->n =obs->nmax =0;
        nav->eph =NULL; nav->n =nav->nmax =0;
        nav->geph=NULL; nav->ng=nav->ngmax=0;
        nav->seph=NULL; nav->ns=nav->nsmax=0;
        nepoch=0;

> 2.1.1.1 readrnxt

/\* read rinex obs and nav files ------------------------------------------------

extern int readrnxt(const char \*file, int rcv, gtime\_t ts, gtime\_t te,double tint, const char \*opt, obs\_t \*obs, nav\_t \*nav,sta\_t \*sta)

> 2.1.1.1.1 readrnxfile

调用子函数：

stat=readrnxfile(files\[i\],ts,te,tint,opt,0,rcv,&type,obs,nav,sta);

调用子函数：

if (sta) init\_sta(sta);

> 2.1.1.1.1.1 readrnxfp

stat=readrnxfp(fp,ts,te,tint,opt,flag,index,type,obs,nav,sta);

调用子函数：先读头文件，判断.o或.n；然后进行body读取

> 2.1.1.1.1.1.1 readrnxh

\[.1\] static int readrnxh(FILE \*fp, double \*ver, char \*type, int \*sys, int \*tsys,

char tobs\[\]\[MAXOBSTYPE\]\[4\], nav\_t \*nav, sta\_t \*sta)

先读取文件头1部分，再通过调用子函数读取2部分（见下图）

读取Rinex文件头的代码实现：

while (fgets(buff,MAXRNXLEN,fp)) {
            
            if (strlen(buff)<=60) continue;
            
            else if (strstr(label,"RINEX VERSION / TYPE")) {
                *ver=str2num(buff,0,9);
                *type=*(buff+20);
                
                /* satellite system */
                switch (*(buff+40)) {
                    case ' ':
                    case 'G': *sys=SYS_GPS;  *tsys=TSYS_GPS; break;
                    case 'R': *sys=SYS_GLO;  *tsys=TSYS_UTC; break;
                    case 'E': *sys=SYS_GAL;  *tsys=TSYS_GAL; break; /* v.2.12 */
                    case 'S': *sys=SYS_SBS;  *tsys=TSYS_GPS; break;
                    case 'J': *sys=SYS_QZS;  *tsys=TSYS_QZS; break; /* v.3.02 */
                    case 'C': *sys=SYS_CMP;  *tsys=TSYS_CMP; break; /* v.2.12 */
                    case 'I': *sys=SYS_IRN;  *tsys=TSYS_IRN; break; /* v.3.03 */
                    case 'M': *sys=SYS_NONE; *tsys=TSYS_GPS; break; /* mixed */
                    default :
                        trace(2,"not supported satellite system: %c\n",*(buff+40));
                        break;
                }
                continue;
            }
            else if (strstr(label,"PGM / RUN BY / DATE")) continue;
            else if (strstr(label,"COMMENT")) { /* opt */
                
                /* read cnes wl satellite fractional bias */
                if (strstr(buff,"WIDELANE SATELLITE FRACTIONAL BIASES")||
                    strstr(buff,"WIDELANE SATELLITE FRACTIONNAL BIASES")) {
                    block=1;
                }
                else if (block) {
                    /* cnes/cls grg clock */
                    if (!strncmp(buff,"WL",2)&&(sat=satid2no(buff+3))&&
                        sscanf(buff+40,"%lf",&bias)==1) {
                        nav->wlbias[sat-1]=bias;
                    }
                    /* cnes ppp-wizard clock */
                    else if ((sat=satid2no(buff+1))&&sscanf(buff+6,"%lf",&bias)==1) {
                        nav->wlbias[sat-1]=bias;
                    }
                }
                continue; 
            }
            /* file type */
            switch (*type) {
                case 'O': decode_obsh(fp,buff,*ver,tsys,tobs,nav,sta); break;
                case 'N': decode_navh (buff,nav); break;
                case 'G': decode_gnavh(buff,nav); break;
                case 'H': decode_hnavh(buff,nav); break;
                case 'J': decode_navh (buff,nav); break; /* extension */
                case 'L': decode_navh (buff,nav); break; /* extension */
            }
            if (strstr(label,"END OF HEADER")) return 1;
            
            if (++i>=MAXPOSHEAD&&*type==' ') break; /* no rinex file */
        }

如下示例一个观测文件的头文件：

![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3d1d3VrdTEyMw_size_16_color_FFFFFF_t_70 1][]

其中，数据类型为：C1 L1 P1 C2 L2 P2 C5 L5；rinex版本为：2.11；在程序中数据类型保存为：

![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3d1d3VrdTEyMw_size_16_color_FFFFFF_t_70 2][]

调用子函数：

/\* decode obs header ---------------------------------------------------------\*/  
static void decode\_obsh(FILE \*fp, char \*buff, double ver, int \*tsys,  
char tobs\[\]\[MAXOBSTYPE\]\[4\], nav\_t \*nav, sta\_t \*sta)

> 2.1.1.1.1.1.2 readrnxobs

\[.2\]读取Rinex Body /\* read rinex body \*/

具体实现代码如下：

switch (*type) {
            case 'O': return readrnxobs(fp,ts,te,tint,opt,index,ver,&tsys,tobs,obs,
                                        sta);
            case 'N': return readrnxnav(fp,opt,ver,sys    ,nav);
            case 'G': return readrnxnav(fp,opt,ver,SYS_GLO,nav);
            case 'H': return readrnxnav(fp,opt,ver,SYS_SBS,nav);
            case 'J': return readrnxnav(fp,opt,ver,SYS_QZS,nav); /* extension */
            case 'L': return readrnxnav(fp,opt,ver,SYS_GAL,nav); /* extension */
            case 'C': return readrnxclk(fp,opt,index,nav);
        }

读取观测文件，首先为： obsd\_t \*data; 分配内存：

if (!(data=(obsd\_t \*)malloc(sizeof(obsd\_t)\*MAXOBS))) return 0;

其中; \#define MAXOBS 64 /\* max number of obs in an epoch \*/

> 2.1.1.1.1.1.2.1 readrnxobsb 循环读取观测数据：

static int readrnxobsb(FILE \*fp, const char \*opt, double ver, int \*tsys,char tobs\[\]\[MAXOBSTYPE\]\[4\], int \*flag, obsd\_t \*data,  
sta\_t \*sta)

/* read rinex obs data body */
        while ((n=readrnxobsb(fp,opt,ver,tsys,tobs,&flag,data,sta))>=0&&stat>=0) {
            
            for (i=0;i<n;i++) {
                
                /* utc -> gpst */
                if (*tsys==TSYS_UTC) data[i].time=utc2gpst(data[i].time);
                
                /* save cycle-slip */
                saveslips(slips,data+i);
            }
            /* screen data by time */
            if (n>0&&!screent(data[0].time,ts,te,tint)) continue;
            
            for (i=0;i<n;i++) {
                
                /* restore cycle-slip */
                restslips(slips,data+i);
                
                data[i].rcv=(unsigned char)rcv;
                
                /* save obs data */
                if ((stat=addobsdata(obs,data+i))<0) break;
            }
        }

读入Rinex观测文件每一个历元的时间后，进入子函数（decode\_obsepoch）：读取每一个历元的卫星个数，即下图中的第一个红框

/\* decode obs epoch ----------------------------------------------------------\*/  
static int decode\_obsepoch(FILE \*fp, char \*buff, double ver, gtime\_t \*time,int \*flag, int \*sats)

![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3d1d3VrdTEyMw_size_16_color_FFFFFF_t_70 3][]

读取上图中的第二个红框：

子函数： /\* decode obs data \*/  
if (decode\_obsdata(fp,buff,ver,mask,index,data+n)&&n<MAXOBS) n++;

每一历元的观测数据读取完成后，对LLI标识的周跳进行处理：

![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3d1d3VrdTEyMw_size_16_color_FFFFFF_t_70 4][]

**输入文件中，第一个为移动站观测文件；当读取完第一个观测文件后，又开始如下图的循环：**

![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3d1d3VrdTEyMw_size_16_color_FFFFFF_t_70 5][]

本例程中，存在两个输入文件分别为：移动站观测文件、移动站导航文件；

注意点：当读取到导航文件头文件中的GPS电离层α、β参数时，

![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3d1d3VrdTEyMw_size_16_color_FFFFFF_t_70 6][]

程序跳转到此处：

![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3d1d3VrdTEyMw_size_16_color_FFFFFF_t_70 7][]

将观测、导航文件的数据读取出来，并保存到相应的结构体中：obs\_t \*obs, nav\_t \*nav

下一步就是：

/\* sort observation data \*/  
nepoch=sortobs(obs);

4、/\* read dcb parameters \*/

保存在：nav\_t \*nav 里面的cbias、rbias中：

> double cbias\[MAXSAT\]\[3\]; /\* satellite dcb (0:p1-p2,1:p1-c1,2:p2-c2) (m) \*/  
> double rbias\[MAXRCV\]\[2\]\[3\]; /\* receiver dcb (0:p1-p2,1:p1-c1,2:p2-c2) (m) \*/

5、/\* set antenna paramters 设置天线参数\*/

6、/\* read ocean tide loading parameters 设置潮汐参数\*/

7、/\* rover/reference fixed position 获取基准站坐标\*/

#define PMODE_SINGLE 0                  /* positioning mode: single */
    #define PMODE_DGPS   1                  /* positioning mode: DGPS/DGNSS */
    #define PMODE_KINEMA 2                  /* positioning mode: kinematic */
    #define PMODE_STATIC 3                  /* positioning mode: static */
    #define PMODE_MOVEB  4                  /* positioning mode: moving-base */
    #define PMODE_FIXED  5                  /* positioning mode: fixed */
    #define PMODE_PPP_KINEMA 6              /* positioning mode: PPP-kinemaric */
    #define PMODE_PPP_STATIC 7              /* positioning mode: PPP-static */
    #define PMODE_PPP_FIXED 8               /* positioning mode: PPP-fixed */

#define POSOPT_POS   0                  /* pos option: LLH/XYZ */
    #define POSOPT_SINGLE 1                 /* pos option: average of single pos */
    #define POSOPT_FILE  2                  /* pos option: read from pos file */
    #define POSOPT_RINEX 3                  /* pos option: rinex header pos */
    #define POSOPT_RTCM  4                  /* pos option: rtcm station pos */
    #define POSOPT_RAW   5                  /* pos option: raw station pos */

8、/\* open solution statistics \*/输出统计信息：solution statistics level (0:off,1:states,2:residuals)

建立路径和保存文件名！

if (flag&&sopt->sstat>0) {
            strcpy(statfile,outfile);
            strcat(statfile,".stat");
            rtkclosestat();
            rtkopenstat(statfile,sopt->sstat);
        }

\* $POS,week,tow,stat,posx,posy,posz,posxf,posyf,poszf  
\* week/tow : gps week no/time of week (s)  
\* stat : solution status  
\* posx/posy/posz : position x/y/z ecef (m) float  
\* posxf/posyf/poszf : position x/y/z ecef (m) fixed  
\*  
\* $VELACC,week,tow,stat,vele,veln,velu,acce,accn,accu,velef,velnf,veluf,accef,accnf,accuf  
\* week/tow : gps week no/time of week (s)  
\* stat : solution status  
\* vele/veln/velu : velocity e/n/u (m/s) float  
\* acce/accn/accu : acceleration e/n/u (m/s^2) float  
\* velef/velnf/veluf : velocity e/n/u (m/s) fixed  
\* accef/accnf/accuf : acceleration e/n/u (m/s^2) fixed  
\*  
\* $CLK,week,tow,stat,clk1,clk2,clk3,clk4  
\* week/tow : gps week no/time of week (s)  
\* stat : solution status  
\* clk1 : receiver clock bias GPS (ns)  
\* clk2 : receiver clock bias GLO-GPS (ns)  
\* clk3 : receiver clock bias GAL-GPS (ns)  
\* clk4 : receiver clock bias BDS-GPS (ns)

\*  
\* $ION,week,tow,stat,sat,az,el,ion,ion-fixed  
\* week/tow : gps week no/time of week (s)  
\* stat : solution status  
\* sat : satellite id  
\* az/el : azimuth/elevation angle(deg)  
\* ion : vertical ionospheric delay L1 (m) float  
\* ion-fixed: vertical ionospheric delay L1 (m) fixed  
\*  
\* $TROP,week,tow,stat,rcv,ztd,ztdf  
\* week/tow : gps week no/time of week (s)  
\* stat : solution status  
\* rcv : receiver (1:rover,2:base station)  
\* ztd : zenith total delay (m) float  
\* ztdf : zenith total delay (m) fixed  
\*

\* $HWBIAS,week,tow,stat,frq,bias,biasf  
\* week/tow : gps week no/time of week (s)  
\* stat : solution status  
\* frq : frequency (1:L1,2:L2,...)  
\* bias : h/w bias coefficient (m/MHz) float  
\* biasf : h/w bias coefficient (m/MHz) fixed  
\*  
\* $SAT,week,tow,sat,frq,az,el,resp,resc,vsat,snr,fix,slip,lock,outc,slipc,rejc  
\* week/tow : gps week no/time of week (s)  
\* sat/frq : satellite id/frequency (1:L1,2:L2,...)  
\* az/el : azimuth/elevation angle (deg)  
\* resp : pseudorange residual (m)  
\* resc : carrier-phase residual (m)  
\* vsat : valid data flag (0:invalid,1:valid)  
\* snr : signal strength (dbHz)  
\* fix : ambiguity flag (0:no data,1:float,2:fixed,3:hold,4:ppp)  
\* slip : cycle-slip flag (bit1:slip,bit2:parity unknown)  
\* lock : carrier-lock count  
\* outc : data outage count  
\* slipc : cycle-slip count  
\* rejc : data reject (outlier) count  
\*  
\*-----------------------------------------------------------------------------\*/

9、 /\* write header to output file \*/

输出结果数据格式如下：

#define SOLF_LLH    0                   /* solution format: lat/lon/height */
    #define SOLF_XYZ    1                   /* solution format: x/y/z-ecef */
    #define SOLF_ENU    2                   /* solution format: e/n/u-baseline */
    #define SOLF_NMEA   3                   /* solution format: NMEA-183 */

![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3d1d3VrdTEyMw_size_16_color_FFFFFF_t_70 8][]

输出 结果文件头：

![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3d1d3VrdTEyMw_size_16_color_FFFFFF_t_70 9][]

输出结果类型结构体：sol\_t /\* solution type \*/

typedef struct {        /* solution type */
        gtime_t time;       /* time (GPST) */
        double rr[6];       /* position/velocity (m|m/s) */
                            /* {x,y,z,vx,vy,vz} or {e,n,u,ve,vn,vu} */
        float  qr[6];       /* position variance/covariance (m^2) */
                            /* {c_xx,c_yy,c_zz,c_xy,c_yz,c_zx} or */
                            /* {c_ee,c_nn,c_uu,c_en,c_nu,c_ue} */
        float  qv[6];       /* velocity variance/covariance (m^2/s^2) */
        double dtr[6];      /* receiver clock bias to time systems (s) */
        unsigned char type; /* type (0:xyz-ecef,1:enu-baseline) */
        unsigned char stat; /* solution status (SOLQ_???) */
        unsigned char ns;   /* number of valid satellites */
        float age;          /* age of differential (s) */
        float ratio;        /* AR ratio factor for valiation */
        float thres;        /* AR ratio threshold for valiation */
    } sol_t;

RTK配置信息结构体：rtk\_t; /\* RTK control/result type \*/

typedef struct {        /* RTK control/result type */
        sol_t  sol;         /* RTK solution */
        double rb[6];       /* base position/velocity (ecef) (m|m/s) */
        int nx,na;          /* number of float states/fixed states */
        double tt;          /* time difference between current and previous (s) */
        double *x, *P;      /* float states and their covariance */
        double *xa,*Pa;     /* fixed states and their covariance */
        int nfix;           /* number of continuous fixes of ambiguity */
        ambc_t ambc[MAXSAT]; /* ambibuity control */
        ssat_t ssat[MAXSAT]; /* satellite status */
        int neb;            /* bytes in error message buffer */
        char errbuf[MAXERRMSG]; /* error message buffer */
        prcopt_t opt;       /* processing options */
    } rtk_t;

ambc\_t /\* ambiguity control type \*/

typedef struct {        /* ambiguity control type */
        gtime_t epoch[4];   /* last epoch */
        int n[4];           /* number of epochs */
        double LC [4];      /* linear combination average */
        double LCv[4];      /* linear combination variance */
        int fixcnt;         /* fix count */
        char flags[MAXSAT]; /* fix flags */
    } ambc_t;

卫星状态信息： ssat\_t; /\* satellite status type \*/

typedef struct {        /* satellite status type */
        unsigned char sys;  /* navigation system */
        unsigned char vs;   /* valid satellite flag single */
        double azel[2];     /* azimuth/elevation angles {az,el} (rad) */
        double resp[NFREQ]; /* residuals of pseudorange (m) */
        double resc[NFREQ]; /* residuals of carrier-phase (m) */
        unsigned char vsat[NFREQ]; /* valid satellite flag */
        unsigned char snr [NFREQ]; /* signal strength (0.25 dBHz) */
        unsigned char fix [NFREQ]; /* ambiguity fix flag (1:fix,2:float,3:hold) */
        unsigned char slip[NFREQ]; /* cycle-slip flag */
        unsigned char half[NFREQ]; /* half-cycle valid flag */
        int lock [NFREQ];   /* lock counter of phase */
        unsigned int outc [NFREQ]; /* obs outage counter of phase */
        unsigned int slipc[NFREQ]; /* cycle-slip counter */
        unsigned int rejc [NFREQ]; /* reject counter */
        double  gf;         /* geometry-free phase L1-L2 (m) */
        double  gf2;        /* geometry-free phase L1-L5 (m) */
        double  mw;         /* MW-LC (m) */
        double  phw;        /* phase windup (cycle) */
        gtime_t pt[2][NFREQ]; /* previous carrier-phase time */
        double  ph[2][NFREQ]; /* previous carrier-phase observable (cycle) */
    } ssat_t;

11、单点定位&前向处理

if (popt\_.mode==PMODE\_SINGLE||popt\_.soltype==0)

执行子函数：procpos

/\* process positioning -------------------------------------------------------\*/  
static void procpos(FILE \*fp, const prcopt\_t \*popt, const solopt\_t \*sopt,int mode)

*  # procpos #

rtkinit(&rtk,popt);

## 1、RTKINIT 初始化 ##

/\* initialize rtk control ------------------------------------------------------  
\* initialize rtk control struct  
\* args : rtk\_t \*rtk IO rtk control/result struct  
\* prcopt\_t \*opt I positioning options (see rtklib.h)  
\* return : none  
\*-----------------------------------------------------------------------------\*/  
extern void rtkinit(rtk\_t \*rtk, const prcopt\_t \*opt)

## 2、当前历元观测卫星多少、观测值； ##

/\* input obs data, navigation messages and sbas correction -------------------\*/  
static int inputobs(obsd\_t \*obs, int solq, const prcopt\_t \*popt)

循环执行：每一个历元、每一个历元进行:

while ((nobs=inputobs(obs,rtk.sol.stat,popt))>=0) {
            
            /* exclude satellites */
            for (i=n=0;i<nobs;i++) {
                if ((satsys(obs[i].sat,NULL)&popt->navsys)&&
                    popt->exsats[obs[i].sat-1]!=1) obs[n++]=obs[i]; //去除配置中不参与的卫星的观测数据
            }
            if (n<=0) continue;
            
            /* carrier-phase bias correction */
            if (navs.nf>0) {    //nf 为 number of satellite fcb data 
                corr_phase_bias_fcb(obs,n,&navs);
            }
            else if (!strstr(popt->pppopt,"-DIS_FCB")) {
                corr_phase_bias_ssr(obs,n,&navs);
            }
            /* disable L2 */
    #if 0
            if (popt->freqopt==1) {
                for (i=0;i<n;i++) obs[i].L[1]=obs[i].P[1]=0.0;
            }
    #endif
            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.1 rtkpos #

extern int rtkpos(rtk\_t \*rtk, const obsd\_t \*obs, int n, const nav\_t \*nav)

> ## 2.1.1 pntpos 单点定位 ##

/\* single-point positioning ----------------------------------------------------  
\* compute receiver position, velocity, clock bias by single-point positioning  
\* with pseudorange and doppler observables

extern 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)

> ### 2.1.1.1 satposs 计算卫星位置和钟差 ###

/\* satellite positions and clocks ----------------------------------------------  
\* compute satellite positions, velocities and clocks  
extern 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)

> 2.1.1.1.1 ephclk 利用广播星历计算卫星钟差

/\* satellite clock with broadcast ephemeris ----------------------------------\*/  
static int ephclk(gtime\_t time, gtime\_t teph, int sat, const nav\_t \*nav,double \*dts)

调用子函数：

1、首先选择与当前历元观测量对应的卫星星历：

/\* select ephememeris --------------------------------------------------------\*/  
static eph\_t \*seleph(gtime\_t time, int sat, int iode, const nav\_t \*nav)

2、计算观测量对应卫星的钟差

/\* broadcast ephemeris to satellite clock bias ---------------------------------

extern double eph2clk(gtime\_t time, const eph\_t \*eph)

> 2.1.1.1.2 **satpos**/\* satellite position and clock at transmission time \*/

/\* satellite position and clock -----------------\* compute satellite position, velocity and clock

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)

注意：satpos中，根据系统对星历的配置，进行不同的卫星位置、速度求解！

#define EPHOPT_BRDC 0                   /* ephemeris option: broadcast ephemeris */
    #define EPHOPT_PREC 1                   /* ephemeris option: precise ephemeris */
    #define EPHOPT_SBAS 2                   /* ephemeris option: broadcast + SBAS */
    #define EPHOPT_SSRAPC 3                 /* ephemeris option: broadcast + SSR_APC */
    #define EPHOPT_SSRCOM 4                 /* ephemeris option: broadcast + SSR_COM */
    #define EPHOPT_LEX  5                   /* ephemeris option: QZSS LEX ephemeris */

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);
            case 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;
            case EPHOPT_LEX   :
                if (!lexeph2pos(time,sat,nav,rs,dts,var)) break; else return 1;
        }

调用子函数：

1、根据设置的星历类型，进行对应的卫星位置、钟差计算：

/\* satellite position and clock by broadcast ephemeris -----------------------\*/  
static int ephpos(gtime\_t time, gtime\_t teph, int sat, const nav\_t \*nav,int iode, double \*rs, double \*dts, double \*var, int \*svh)

1.1 选择观测量对应的卫星星历：

/\* select ephememeris --------------------------------------------------------\*/  
static eph\_t \*seleph(gtime\_t time, int sat, int iode, const nav\_t \*nav)

2、计算卫星位置、速度、钟差

/\* broadcast ephemeris to satellite position and clock bias --------------------

extern void eph2pos(gtime\_t time, const eph\_t \*eph, double \*rs, double \*dts,double \*var)

**2.1.1.1 节对应的函数调用顺序如下：**

![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3d1d3VrdTEyMw_size_16_color_FFFFFF_t_70 10][]

> ### 2.1.1.2 estpos 利用伪距计算接收机位置 ###

/\* estimate receiver position ------------------------------------------------\*/  
static 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)

估计参数7个：3个位置、1个钟差、3个系统间钟差参数；

x\[NX\]=\{0\},dx\[NX\],Q\[NX\*NX\]

v=mat(n+4,1); H=mat(NX,n+4); var=mat(n+4,1);

**下面进行循环迭代：求解接收机位置、速度；迭代终止条件为：(norm(dx,NX)<1E-4)**

> 2.1.1.2.1 伪距残差

/\* pseudorange residuals -----------------------------------------------------\*/  
static 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)

调用如下子函数：

1、计算卫星-接收机之间几何距离：

主要去除地球自转影响！

/\* geometric distance ----------------------------------------------------------  
\* compute geometric distance and receiver-to-satellite unit vector

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

2、对伪距进行码偏差改正: 单频点定位考虑TGD

/\* psendorange with code bias correction -------------------------------------\*/  
static double prange(const obsd\_t \*obs, const nav\_t \*nav, const double \*azel,  
int iter, const prcopt\_t \*opt, double \*var)

3、剔除预先设置中不参与定位的卫星

/\* test excluded satellite -----------------------------------------------------

extern int satexclude(int sat, int svh, const prcopt\_t \*opt)

4、电离层改正：

/\* ionospheric correction ------------------------------------------------------  
\* compute ionospheric correction

extern int ionocorr(gtime\_t time, const nav\_t \*nav, int sat, const double \*pos,const double \*azel, int ionoopt, double \*ion, double \*var)

调用子函数：广播星历中电离层α、β模型：

5、对流层改正：

/\* tropospheric correction -----------------------------------------------------  
\* compute tropospheric correction

extern int tropcorr(gtime\_t time, const nav\_t \*nav, const double \*pos, const double \*azel, int tropopt, double \*trp, double \*var)

调用子函数：compute tropospheric delay by standard atmosphere and saastamoinen model

6、伪距残差：

/\* pseudorange residual \*/ //nv表示方程的个数  
v\[nv\]=P-(r+dtr-CLIGHT\*dts\[i\*2\]+dion+dtrp);

7、设计/构造矩阵 H

/\* design matrix \*/ //  
for (j=0;j<NX;j++) H\[j+nv\*NX\]=j<3?-e\[j\]:(j==3?1.0:0.0);

8、如果是多系统参与定位，需要估计系统间时间偏差

/\* time system and receiver bias offset correction \*/  
if (sys==SYS\_GLO) \{v\[nv\]-=x\[4\]; H\[4+nv\*NX\]=1.0; mask\[1\]=1;\}  
else if (sys==SYS\_GAL) \{v\[nv\]-=x\[5\]; H\[5+nv\*NX\]=1.0; mask\[2\]=1;\}  
else if (sys==SYS\_CMP) \{v\[nv\]-=x\[6\]; H\[6+nv\*NX\]=1.0; mask\[3\]=1;\}  
else mask\[0\]=1;

其具体参数为：

> sol->dtr\[0\]=x\[3\]/CLIGHT; /\* receiver clock bias (s) \*/  
> sol->dtr\[1\]=x\[4\]/CLIGHT; /\* glo-gps time offset (s) \*/  
> sol->dtr\[2\]=x\[5\]/CLIGHT; /\* gal-gps time offset (s) \*/  
> sol->dtr\[3\]=x\[6\]/CLIGHT; /\* bds-gps time offset (s) \*/

9、伪距测量中各部分误差项中方差和

/\* pseudorange measurement error variance ------------------------------------\*/  
static double varerr(const prcopt\_t \*opt, double el, int sys)

函数代码：

/* pseudorange measurement error variance ------------------------------------*/
    static double varerr(const prcopt_t *opt, double el, int sys)
    {
        double fact,varr;
        fact=sys==SYS_GLO?EFACT_GLO:(sys==SYS_SBS?EFACT_SBS:EFACT_GPS);
        varr=SQR(opt->err[0])*(SQR(opt->err[1])+SQR(opt->err[2])/sin(el));
        if (opt->ionoopt==IONOOPT_IFLC) varr*=SQR(3.0); /* iono-free */
        return SQR(fact)*varr;
    }

对应的项为：

结构体：prcopt\_t

double err\[5\]; /\* measurement error factor \*/  
/\* \[0\]:reserved \*/  
/\* \[1-3\]:error factor a/b/c of phase (m) \*/  
/\* \[4\]:doppler frequency (hz) \*/

默认参数配置为：\{100.0,0.003,0.003,0.0,1.0\}, /\* err\[\] \*/

为了保持“卫星个数”大于参数个数(7个)，保持为超定方程：设置H矩阵及对应的方差，代码如下所示：

/* constraint to avoid rank-deficient */
        for (i=0;i<4;i++) {
            if (mask[i]) continue;
            v[nv]=0.0;
            for (j = 0; j < NX; j++)
            {
                H[j + nv * NX] = j == i + 3 ? 1.0 : 0.0;
                printf("%f ", H[j + nv * NX]);
                if (j == NX - 1) printf("\n");
            }
            var[nv++]=0.01;
            //printf("%f ",H[j + nv * NX]);
        }

某一个历元观测到卫星个数为5，故需要增加3个“伪卫星”，构成超定方程；代码运行如下：

![20200514141342621.png][]

> 2.1.1.2.2 对伪距残差进行定权，高度角定权：

/* weight by variance */
    		printf("伪距误差标准差：\n");
            for (j=0;j<nv;j++) {
                sig=sqrt(var[j]);
                v[j]/=sig;
    			printf("sig[%d]=%f%t ",j,sig);
                for (k=0;k<NX;k++) H[k+j*NX]/=sig;
            }

> 2.1.1.2.3 加权最小二乘：

/\* least square estimation 用最小二乘法，解Ax=b-----------------------------------------------------  
\* least square estimation by solving normal equation (x=(A\*A')^-1\*A\*y)

\* least square estimation by solving normal equation (x=(A\*A')^-1\*A\*y)  
\* args : 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) m为方程个数，测量个数  
\* double \*x O estmated parameters (n x 1)  
\* double \*Q O esimated parameters covariance matrix (n x n)  
\* return : status (0:ok,0>:error)  
\* notes : for weighted least square, replace A and y by A\*w and w\*y (w=W^(1/2))  
\* matirix stored by column-major order (fortran convention)  
\*-----------------------------------------------------------------------------\*/

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

> 2.1.1.2.4 保存结果：接收机位置、钟差、方差、协方差、参与定位卫星个数：

if (norm(dx,NX)<1E-4) {
                sol->type=0;
                sol->time=timeadd(obs[0].time,-x[3]/CLIGHT); //求解位置的时刻为：扣除接收机钟差影响
                sol->dtr[0]=x[3]/CLIGHT; /* receiver clock bias (s) */
                sol->dtr[1]=x[4]/CLIGHT; /* glo-gps time offset (s) */
                sol->dtr[2]=x[5]/CLIGHT; /* gal-gps time offset (s) */
                sol->dtr[3]=x[6]/CLIGHT; /* bds-gps time offset (s) */
                for (j=0;j<6;j++) sol->rr[j]=j<3?x[j]:0.0;
                for (j=0;j<3;j++) sol->qr[j]=(float)Q[j+j*NX];
                sol->qr[3]=(float)Q[1];    /* cov xy */
                sol->qr[4]=(float)Q[2+NX]; /* cov yz */
                sol->qr[5]=(float)Q[2];    /* cov zx */
                sol->ns=(unsigned char)ns;
                sol->age=sol->ratio=0.0;
                
                /* validate solution */
                if ((stat=valsol(azel,vsat,n,opt,v,nv,NX,msg))) {
                    sol->stat=opt->sateph==EPHOPT_SBAS?SOLQ_SBAS:SOLQ_SINGLE;
                }

同时还包括：有效性判断：

1、chisqr

2、最大DOP值

/\* validate solution ---------------------------------------------------------\*/  
static int valsol(const double \*azel, const int \*vsat, int n,  
const prcopt\_t \*opt, const double \*v, int nv, int nx,  
char \*msg)

**2.1.1.2 estpos 对应的函数调用层级如下：**

![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3d1d3VrdTEyMw_size_16_color_FFFFFF_t_70 11][]

> ### 2.1.1.3 raim\_fde 有效性判断： ###

进入raim\_fde前进行的判断：

/\* raim fde \*/  
if (!stat&&n>=6&&opt->posopt\[4\]) \{ stat=raim\_fde(obs,n,rs,dts,var,svh,nav,&opt\_,sol,azel\_,vsat,resp,msg);\}

1、参与定位卫星个数6个及以上；

2、opt->posopt\[4\] 为 1 ？表示啥意思？

/\* raim fde (failure detection and exclution) -------------------------------\*/  
static 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)

> ### 2.1.1.4 estvel 计算接收机速度： ###

/\* estimate receiver velocity ------------------------------------------------\*/  
static 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)

调用子函数：

1、多普勒残差：

/\* doppler residuals ---------------------------------------------------------\*/  
static 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)

2、最小二乘计算速度：

> 2.1.1.4 ssat 保存相关信息：

ssat\_t ssat;

> ## 2.1.2 outsolstat ##

/\* output solution status ----------------------------------------------------\*/  
static void outsolstat(rtk\_t \*rtk)

调用子函数：

> ### 2.1.2.1 rtkoutstat ###

/\* write solution status to buffer -------------------------------------------\*/  
extern int rtkoutstat(rtk\_t \*rtk, char \*buff)

保存的为如下信息到文件：myposbyVS.pos.stat

/* write ppp solution status to buffer */
        if (rtk->opt.mode>=PMODE_PPP_KINEMA) {
            return pppoutstat(rtk,buff);
        }
        est=rtk->opt.mode>=PMODE_DGPS;
        nfreq=est?nf:1;
        tow=time2gpst(rtk->sol.time,&week);
        
        /* receiver position */
        if (est) {
            for (i=0;i<3;i++) xa[i]=i<rtk->na?rtk->xa[i]:0.0;
            p+=sprintf(p,"$POS,%d,%.3f,%d,%.4f,%.4f,%.4f,%.4f,%.4f,%.4f\n",week,tow,
                       rtk->sol.stat,rtk->x[0],rtk->x[1],rtk->x[2],xa[0],xa[1],
                       xa[2]);
        }
        else {
            p+=sprintf(p,"$POS,%d,%.3f,%d,%.4f,%.4f,%.4f,%.4f,%.4f,%.4f\n",week,tow,
                       rtk->sol.stat,rtk->sol.rr[0],rtk->sol.rr[1],rtk->sol.rr[2],
                       0.0,0.0,0.0);
        }
        /* receiver velocity and acceleration */
        if (est&&rtk->opt.dynamics) {
            ecef2pos(rtk->sol.rr,pos);
            ecef2enu(pos,rtk->x+3,vel);
            ecef2enu(pos,rtk->x+6,acc);
            if (rtk->na>=6) ecef2enu(pos,rtk->xa+3,vela);
            if (rtk->na>=9) ecef2enu(pos,rtk->xa+6,acca);
            p+=sprintf(p,"$VELACC,%d,%.3f,%d,%.4f,%.4f,%.4f,%.5f,%.5f,%.5f,%.4f,%.4f,%.4f,%.5f,%.5f,%.5f\n",
                       week,tow,rtk->sol.stat,vel[0],vel[1],vel[2],acc[0],acc[1],
                       acc[2],vela[0],vela[1],vela[2],acca[0],acca[1],acca[2]);
        }
        else {
            ecef2pos(rtk->sol.rr,pos);
            ecef2enu(pos,rtk->sol.rr+3,vel);
            p+=sprintf(p,"$VELACC,%d,%.3f,%d,%.4f,%.4f,%.4f,%.5f,%.5f,%.5f,%.4f,%.4f,%.4f,%.5f,%.5f,%.5f\n",
                       week,tow,rtk->sol.stat,vel[0],vel[1],vel[2],
                       0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0);
        }
        /* receiver clocks */
        p+=sprintf(p,"$CLK,%d,%.3f,%d,%d,%.3f,%.3f,%.3f,%.3f\n",
                   week,tow,rtk->sol.stat,1,rtk->sol.dtr[0]*1E9,rtk->sol.dtr[1]*1E9,
                   rtk->sol.dtr[2]*1E9,rtk->sol.dtr[3]*1E9);
        
        /* ionospheric parameters */
        if (est&&rtk->opt.ionoopt==IONOOPT_EST) {
            for (i=0;i<MAXSAT;i++) {
                ssat=rtk->ssat+i;
                if (!ssat->vs) continue;
                satno2id(i+1,id);
                j=II(i+1,&rtk->opt);
                xa[0]=j<rtk->na?rtk->xa[j]:0.0;
                p+=sprintf(p,"$ION,%d,%.3f,%d,%s,%.1f,%.1f,%.4f,%.4f\n",week,tow,
                           rtk->sol.stat,id,ssat->azel[0]*R2D,ssat->azel[1]*R2D,
                           rtk->x[j],xa[0]);
            }
        }
        /* tropospheric parameters */
        if (est&&(rtk->opt.tropopt==TROPOPT_EST||rtk->opt.tropopt==TROPOPT_ESTG)) {
            for (i=0;i<2;i++) {
                j=IT(i,&rtk->opt);
                xa[0]=j<rtk->na?rtk->xa[j]:0.0;
                p+=sprintf(p,"$TROP,%d,%.3f,%d,%d,%.4f,%.4f\n",week,tow,
                           rtk->sol.stat,i+1,rtk->x[j],xa[0]);
            }
        }
        /* receiver h/w bias */
        if (est&&rtk->opt.glomodear==2) {
            for (i=0;i<nfreq;i++) {
                j=IL(i,&rtk->opt);
                xa[0]=j<rtk->na?rtk->xa[j]:0.0;
                p+=sprintf(p,"$HWBIAS,%d,%.3f,%d,%d,%.4f,%.4f\n",week,tow,
                           rtk->sol.stat,i+1,rtk->x[j],xa[0]);
            }
        }

# 2.2 outsol #

/\* output solution body --------------------------------------------------------  
\* output solution body to file  
\* args : FILE \*fp I output file pointer  
\* sol\_t \*sol I solution  
\* double \*rb I base station position \{x,y,z\} (ecef) (m)  
\* solopt\_t \*opt I solution options  
\* return : none  
\*-----------------------------------------------------------------------------\*/  
extern void outsol(FILE \*fp, const sol\_t \*sol, const double \*rb, const solopt\_t \*opt)

调用子函数：

> ## 2.2.1 outsols ##

/\* output solution body --------------------------------------------------------  
\* output solution body to buffer  
\* args : unsigned char \*buff IO output buffer  
\* sol\_t \*sol I solution  
\* double \*rb I base station position \{x,y,z\} (ecef) (m)  
\* solopt\_t \*opt I solution options  
\* return : number of output bytes  
\*-----------------------------------------------------------------------------\*/  
extern int outsols(unsigned char \*buff, const sol\_t \*sol, const double \*rb,const solopt\_t \*opt)

> ### 2.2.1.1 outpos ###

判断语句：根据设置的输出位置形式，进行选择：

switch (opt->posf) {
            case SOLF_LLH:  p+=outpos (p,s,sol,opt);   break;
            case SOLF_XYZ:  p+=outecef(p,s,sol,opt);   break;
            case SOLF_ENU:  p+=outenu(p,s,sol,rb,opt); break;
            case SOLF_NMEA: p+=outnmea_rmc(p,sol);
                            p+=outnmea_gga(p,sol); break;
        }

/\* output solution as the form of lat/lon/height -----------------------------\*/  
static int outpos(unsigned char \*buff, const char \*s, const sol\_t \*sol,const solopt\_t \*opt)

**2.2节函数调用层级，如下：**

![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3d1d3VrdTEyMw_size_16_color_FFFFFF_t_70 12][]

**GPS伪距单频/L1 单点定位整个过程包括：**

*  1、读取移动站观测数据、导航数据；将整个文件的相关数据放入结构体中；进行排序和重复数据删除；
 *  2、进行单点定位，测速；保存数据.pos、.pos.trace、.pos.stat文件；
 *  3、最终freeobsnav(&obss,&navs);

# 参考： #

[RTKLIB 源码阅读笔记][RTKLIB]

[RTKLIB源码解析（一）——单点定位(pntpos.c)][RTKLIB_pntpos.c]

RTKLIB-MANUAL（手册）

[20200512094440399.png]: /images/20230312/202ea95c273f4257b2bfa2475e3325db.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3d1d3VrdTEyMw_size_16_color_FFFFFF_t_70]: /images/20230312/225f96c9437245e3b1ac0eb11933f24a.png
[20200512105333577.png]: /images/20230312/3e29108ba0524cbfbde580f79516f639.png
[20200512105402228.png]: /images/20230312/d9e97d6248e9446e8a825607a63552e6.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3d1d3VrdTEyMw_size_16_color_FFFFFF_t_70 1]: /images/20230312/cf32c230ad9e4851b23babd372a53d73.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3d1d3VrdTEyMw_size_16_color_FFFFFF_t_70 2]: /images/20230312/3ffe27a52e8049c793fd74d14ffa6626.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3d1d3VrdTEyMw_size_16_color_FFFFFF_t_70 3]: /images/20230312/151ba1a7ebd2409bb9d72361d10c7fd6.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3d1d3VrdTEyMw_size_16_color_FFFFFF_t_70 4]: /images/20230312/8c2c5fa2d213446dacab4b88735679a8.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3d1d3VrdTEyMw_size_16_color_FFFFFF_t_70 5]: /images/20230312/41ef062e9f9c4d868d434b6742a11588.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3d1d3VrdTEyMw_size_16_color_FFFFFF_t_70 6]: /images/20230312/9d5448ee1e964f2caa8f006d229cf716.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3d1d3VrdTEyMw_size_16_color_FFFFFF_t_70 7]: /images/20230312/7e84149d3ced45df93e319a0ea0a0d23.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3d1d3VrdTEyMw_size_16_color_FFFFFF_t_70 8]: /images/20230312/6f47f9fb5ddc4eba9933bb3611ae0a87.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3d1d3VrdTEyMw_size_16_color_FFFFFF_t_70 9]: /images/20230312/8eab49f3791640feb6da91f938b200f9.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3d1d3VrdTEyMw_size_16_color_FFFFFF_t_70 10]: /images/20230312/d6a648ecf92f4df79faacf4217464a4d.png
[20200514141342621.png]: /images/20230312/9d691b971a8d44ec877123a57aea1685.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3d1d3VrdTEyMw_size_16_color_FFFFFF_t_70 11]: /images/20230312/5eaaaad90a3c4efd86a3c1615e610590.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3d1d3VrdTEyMw_size_16_color_FFFFFF_t_70 12]: /images/20230312/1d92b269c1e341f887feb4ac8351b59b.png
[RTKLIB]: https://zhuanlan.zhihu.com/p/78359579
[RTKLIB_pntpos.c]: https://www.zybuluo.com/taqikema/note/1101465