RTKLIB学习总结（一）VS配置RTKLIB、manual、矩阵、最小二乘和Kalman滤波基本函数

谢弘阔

2023-12-01

一、前言

本人是导航工程大二的学生，打算暑假自己给我的项目写RTK实时定位算法。大一就听老师说学习RTKLIB的源码，先看懂别人的代码，才能自己写代码，但一直没有系统的开始学。现在掌握了一些理论知识，可以开始学了。我计划以博客的形式总结学习的过程，既是为了对学习的巩固，也为了可以分享学习到的知识。

二、2.4.2manual目录翻译介绍

rtklib的manual解读：有对重要内容更详细的翻译，可以直接Ctrt+F搜索。

Overview：概述、

User Requirements：用户要求、

System Requirements：系统要求、

License：许可、

Instructions ：说明书、

Installation and Uninstallation：安装和卸载

Real‐Time Positioning with RTKNAVI：用RTKNAVI实时定位

Configure Input, Output and Log Streams for RTKNAVI：为RTKNAVI配置输入、输出和日志流

Post‐Processing Analysis with RTKPOS：用RTKPOST后处理分析

Configure Positioning Options for RTKNAVI and RTKPOST：配置RTKNAVI和RTKPOST的定位选项

Convert Receiver Raw Data to RINEX with RTKCONV：用RTKCONV将接收器原始数据转换为RINEX

View and Plot Solutions with RTKPLOT：用RTKPLOT查看和绘制结果

View and Plot Observation Data with RTKPLOT：用RTKPLOT查看和绘制观数据

Download GNSS Products and Data with RTKGET：用RTKGET下载GNSS数据

NTRIP Browser ：NTRIP客户端

Use CUI APs of RTKLIB：用RTKLIB的命令行程序

Build APs or Develop User APs with RTKLIB：用RTKLIB构建程序

Rebuild GUI and CUI APs on Windows：在Windows上重建GUI和CUI程序

Build CUI APs：构建CUI程序

Develop and Link User APs with RTKLIB：开发用户程序

附录A：CUI程序命令参考：RTKRCV、RNX2RTKP、POS2KML、CONVBIN、STR2STR

附录B：文件格式：Positioning Solution File、SBAS Log File、Solution Status File、Configuration File、URL List File for GNSS Data

附录C：API参考，RTKLIB函数的一句话简介

附录D：文件和信息：Supported RINEX Files、Supported Receiver Messages、Supported Signal IDs/Observation Types、Default Priorities of Multiple Signals、Receiver Dependent Input Options

附录E：模型和算法：

Time System：时间系统

Coordinates System：坐标系统

GNSS Signal Measurement Models：信号测量模型

GNSS Satellite Ephemerides and Clocks：卫星星历表和时钟

Troposphere and Ionosphere Models：对流层、电离层模型

Single Point Positioning：SPP

Kinematic, Static and Moving‐Baseline：差分定位

PPP (Precise Point Positioning)：PPP

附录F：GNSS信号详细说明

参考文献

三、将 RTKLIB2.4.3 b34 配置到VS2022中的流程

参考”珞珈山圈哥“的博客: rtklib学习笔记1：在Visual Studio 2019中调试rtklib 2.4.3程序，博客是用VS2019调试，我用VS2022，按他的步骤调试也成功了。
rnx2rtkp的功能是读取rinex文件计算接收机位置

1、下载RTKLIB

在RTKLIB官网选最新版 2.4.3 b34，点Source Programs and Data下面的GitHub进入GitHub页面，点开绿色的Code下拉菜单，再点Download ZIP，下载解压即可。

2、在VS2022中创建空C++项目、导入源码文件

创建C++空项目，可以勾选“解决方案和项目放在统一目录中”，记住创建的项目目录。

把RTKLIB源码文件中整个src文件夹复制到创建的项目文件目录中。

把RTKLIB源码文件中\app\consapp中的rnx2rtkp.c放到刚刚复制过去的src文件夹。

在解决源文件中添加名为“src”的筛选器，再在src筛选器下面添加名为“rcv”的筛选器。

右键添加现有项目把src/rcv文件夹中的所有文件加到src/rcv筛选器中，src中所有代码文件加到src筛选器中。

把主函数rnx2rtkp.c文件中的#include "rtklib.h"修改为#include "./rtklib.h“

把在src/rcv文件夹几个的.c文件中的#include "rtklib.h"修改为#include "../rtklib.h”

3、项目属性设置

打开项目属性，在链接器—输入—附加依赖项中添加依赖库winmm.lib和ws2_32.lib。

配置属性—高级—字符集中设置为用使用多字节字符集。

C/C++中的SDL检查设置为否，附加包含目录添加.\src 、预编译头为不使用预编译头。

预处理器中添加如下内容：

_LIB
_CRT_SECURE_NO_WARNINGS
_WINSOCK_DEPRECATED_NO_WARNINGS             
ENAGLO
DLL
WIN32

尤其主要加WIN32，好多博客都没加这一项，加了这一项后RTKLIB就不会用Linux下的<pthread.h>和<sys/select.h>，咱们项目要在Windows下编译运行的，不加会报”找不到<pthread.h>和<sys/select.h>“的错。

将常规中的目标文件名改为rnx2rtkp 。

改不改都行，默认目标文件名是项目名。

4、改代码的BUG

可能会报“使用了可能未初始化的本地指针变量“sbs”的错误，解决方式是对指针变量进行初始化，将ephemeris.c文件中的第579行改为“const sbssatp_t *sbs=NULL;”。

5、启用所有卫星系统

#define ENAGLO
#define ENAGAL
#define ENACMP
#define ENAQZS
#define ENAIRN

四、前置知识

1、stdarg库

stdarg.h是C语言中C标准函数库的头文件，stdarg是由standard（标准）和 arguments（参数）简化而来，主要目的为让函数能够接收可变参数。用法如下：

首先在函数里定义一具VA_LIST型的变量，这个变量是指向参数的指针。

然后用VA_START宏初始化变量刚定义的VA_LIST变量。

然后用VA_ARG返回可变的参数，VA_ARG的第二个参数是你要返回的参数的类型（如果函数有多个可变参数的，依次调用VA_ARG获取各个参数）

最后用VA_END宏结束可变参数的获取

#include <stdio.h>
#include <stdarg.h>
voidprintargs(intarg1, ...) /* 输出所有int类型的参数，直到-1结束 */
{
    va_listap; //定义一个va_list变量ap
    inti;
    va_start(ap, arg1); //执行ap = (va_list)&v + _INTSIZEOF(v)，ap指向参数v之后的那个参数的地址，即 ap指向第一个可变参数在堆栈的地址。
    for (i=arg1; i!=-1; i=va_arg(ap, int))    //va_arg(ap,t) ， ( *(t *)((ap += _INTSIZEOF(t)) - _INTSIZEOF(t)) )取出当前ap指针所指的值，并使ap指向下一个参数。 ap＋= sizeof(t类型)，让ap指向下一个参数的地址。然后返回ap-sizeof(t类型)的t类型*指针，这正是第一个可变参数在堆栈里的地址。然后 用*取得这个地址的内容。
    printf("%d ", i);
    va_end(ap); //清空va_list ap
    putchar('\n');
}
intmain(void)
{
    printargs(5, 2, 14, 84, 97, 15, 24, 48, -1);
    printargs(84, 51, -1);
    printargs(-1);
    printargs(1, -1);
    while(1);
    return0;
}

2、fatalerr 函数

fatalerr 函数：致命错误的处理，接收可变参数。如果没提供fatalfunc回调函数，则stderr错误消息，并以-9形式退出程序；如果提供了fatalfunc回调函数，则把msg交给fatalfunc回调函数处理，并以-9形式退出程序。

static void fatalerr(const char *format, ...)
{
    char msg[1024];
    va_list ap;
    va_start(ap,format); vsprintf(msg,format,ap); va_end(ap);
    if (fatalfunc) fatalfunc(msg);
    else fprintf(stderr,"%s",msg);
    exit(-9);
}

五、矩阵计算相关函数

GNSS处理的数据都是矩阵数据，RTKLIB的 rtkcmn.c中写了一些矩阵运算的函数。
也可看 manual的 Appendix C API References，里面有RTKLIB函数功能的简单描述

mat() ：New matrix、创建新矩阵，传入行数n，列数m，malloc开辟n*m个double空间，返回此空间的首地址。

extern double *mat(int n, int m)
{
    double *p;
    
    if (n<=0||m<=0) return NULL;
    if (!(p=(double *)malloc(sizeof(double)*n*m))) {
        fatalerr("matrix memory allocation error: n=%d,m=%d\n",n,m);
    }
    return p;
}

imat()、New integer matrix、与mat类似，区别在于iamt创建的矩阵是int类型的。（RTKLIB里的矩阵一般都是double类型的二维数组）

extern int *imat(int n, int m)
{
    int *p;
    
    if (n<=0||m<=0) return NULL;
    if (!(p=(int *)malloc(sizeof(int)*n*m))) {
        fatalerr("integer matrix memory allocation error: n=%d,m=%d\n",n,m);
    }
    return p;
}

zeros():New zero matrix、创建0矩阵、double类型。它比mat多了\#if NOCALLOC if ((p=mat(n,m))) for (n=n*m-1;n>=0;n--) p[n]=0.0;，如果预编译指令有NOCALLOC，会把所有位赋值0.0

extern double *zeros(int n, int m)
{
    double *p;
    
#if NOCALLOC
    if ((p=mat(n,m))) for (n=n*m-1;n>=0;n--) p[n]=0.0;
#else
    if (n<=0||m<=0) return NULL;
    if (!(p=(double *)calloc(sizeof(double),n*m))) {
        fatalerr("matrix memory allocation error: n=%d,m=%d\n",n,m);
    }
#endif
    return p;
}

eye()：New identity matrix、创建单位矩阵，对角线元素全赋值1.0 。

extern double *eye(int n)
{
    double *p;
    int i;
    
    if ((p=zeros(n,n))) for (i=0;i<n;i++) p[i+i*n]=1.0;
    return p;
}

dot()：Inner Product、点乘、求两个向量的内积（对应位元素相乘再相加）

extern double dot(const double *a, const double *b, int n)
{
    double c=0.0;
    
    while (--n>=0) c+=a[n]*b[n];
    return c;
}

norm()：Euclid norm、求向量的I2范数（向量自己的内积再开方）

extern double norm(const double *a, int n)
{
    return sqrt(dot(a,a,n));
}

cross3() 、Outer product of 3D vectors、三维向量的外积

extern void cross3(const double *a, const double *b, double *c)
{
    c[0]=a[1]*b[2]-a[2]*b[1];
    c[1]=a[2]*b[0]-a[0]*b[2];
    c[2]=a[0]*b[1]-a[1]*b[0];
}

normv3() Normalize 3D vector、规格化三维向量（各元素除以向量的I2范数），a为传入向量，b为传出向量。

extern int normv3(const double *a, double *b)
{
    double r;
    if ((r=norm(a,3))<=0.0) return 0;
    b[0]=a[0]/r;
    b[1]=a[1]/r;
    b[2]=a[2]/r;
    return 1;
}

matcpy()：Copy matrix、将B矩阵的所有元素赋值给A矩阵

extern void matcpy(double *A, const double *B, int n, int m)
{
    memcpy(A,B,sizeof(double)*n*m);
}

matmul()：Multiply matrix、矩阵A、B相乘返回C

tr：两个元素的字符数组，Tr、Tr[1]分别标识A、B的状态，取值N正常、T转置，可直接写TN、NN。

n、k、m：A为nm矩阵， B为mk矩阵、A的列数与B的行数相等，所有只用三个参数。作者之所以把n和k放在m前面，看起来顺序反逻辑，其实是矩阵相乘一共要循环n*k次，从编程的顺序来排的。

总体流程为： C=alphaAB+beta*C

extern void matmul(const char *tr, int n, int k, int m, double alpha,
                   const double *A, const double *B, double beta, double *C)
{
    double d;
    int i,j,x,f=tr[0]=='N'?(tr[1]=='N'?1:2):(tr[1]=='N'?3:4);
    
    for (i=0;i<n;i++) for (j=0;j<k;j++) {
        d=0.0;
        switch (f) {
            case 1: for (x=0;x<m;x++) d+=A[i+x*n]*B[x+j*m]; break;
            case 2: for (x=0;x<m;x++) d+=A[i+x*n]*B[j+x*k]; break;
            case 3: for (x=0;x<m;x++) d+=A[x+i*m]*B[x+j*m]; break;
            case 4: for (x=0;x<m;x++) d+=A[x+i*m]*B[j+x*k]; break;
        }
        if (beta==0.0) C[i+j*n]=alpha*d; else C[i+j*n]=alpha*d+beta*C[i+j*n];
    }
}

matinv()：Inverse of matrix、矩阵求逆

extern int matinv(double *A, int n)
{
    double d,*B;
    int i,j,*indx;
    
    indx=imat(n,1); B=mat(n,n); matcpy(B,A,n,n);
    if (ludcmp(B,n,indx,&d)) {free(indx); free(B); return -1;}
    for (j=0;j<n;j++) {
        for (i=0;i<n;i++) A[i+j*n]=0.0;
        A[j+j*n]=1.0;
        lubksb(B,n,indx,A+j*n);
    }
    free(indx); free(B);
    return 0;
}

solve()：Solve linear equation、求方程线性解、AX=Y，求X，A为nn，Y为nm

extern int solve(const char *tr, const double *A, const double *Y, int n,
                 int m, double *X)
{
    double *B=mat(n,n);
    int info;
    
    matcpy(B,A,n,n);
    if (!(info=matinv(B,n))) matmul(tr[0]=='N'?"NN":"TN",n,m,n,1.0,B,Y,0.0,X);
    free(B);
    return info;
}

ludcmp ()：LU分解，即把矩阵A分解LU乘积的形式，U为单位上三角矩阵和L单位为下三角矩阵两部分。

static int ludcmp(double *A, int n, int *indx, double *d)
{
    double big,s,tmp,*vv=mat(n,1);
    int i,imax=0,j,k;
    
    *d=1.0;
    for (i=0;i<n;i++) {
        big=0.0; for (j=0;j<n;j++) if ((tmp=fabs(A[i+j*n]))>big) big=tmp;
        if (big>0.0) vv[i]=1.0/big; else {free(vv); return -1;}
    }
    for (j=0;j<n;j++) {
        for (i=0;i<j;i++) {
            s=A[i+j*n]; for (k=0;k<i;k++) s-=A[i+k*n]*A[k+j*n]; A[i+j*n]=s;
        }
        big=0.0;
        for (i=j;i<n;i++) {
            s=A[i+j*n]; for (k=0;k<j;k++) s-=A[i+k*n]*A[k+j*n]; A[i+j*n]=s;
            if ((tmp=vv[i]*fabs(s))>=big) {big=tmp; imax=i;}
        }
        if (j!=imax) {
            for (k=0;k<n;k++) {
                tmp=A[imax+k*n]; A[imax+k*n]=A[j+k*n]; A[j+k*n]=tmp;
            }
            *d=-(*d); vv[imax]=vv[j];
        }
        indx[j]=imax;
        if (A[j+j*n]==0.0) {free(vv); return -1;}
        if (j!=n-1) {
            tmp=1.0/A[j+j*n]; for (i=j+1;i<n;i++) A[i+j*n]*=tmp;
        }
    }
    free(vv);
    return 0;
}

lubksb ()：LU回代,即把单位上三角矩阵U和单位下三角矩阵L矩阵回代为一个整体矩阵

static void lubksb(const double *A, int n, const int *indx, double *b)
{
    double s;
    int i,ii=-1,ip,j;
    
    for (i=0;i<n;i++) {
        ip=indx[i]; s=b[ip]; b[ip]=b[i];
        if (ii>=0) for (j=ii;j<i;j++) s-=A[i+j*n]*b[j]; else if (s) ii=i;
        b[i]=s;
    }
    for (i=n-1;i>=0;i--) {
        s=b[i]; for (j=i+1;j<n;j++) s-=A[i+j*n]*b[j]; b[i]=s/A[i+i*n];
    }
}

matprint()：Print matrix、打印矩阵到stdout

extern void matprint(const double A[], int n, int m, int p, int q)
{
    matfprint(A,n,m,p,q,stdout);
}

matfprint()：Print matrix to file、打印矩阵到文件中

extern void matfprint(const double A[], int n, int m, int p, int q, FILE *fp)
{
    int i,j;
    
    for (i=0;i<n;i++) {
        for (j=0;j<m;j++) fprintf(fp," %*.*f",p,q,A[i+j*n]);
        fprintf(fp,"\n");
    }
}

六、最小二乘与Kalman滤波

lsq()：Least square estimation、最小二乘估计

A：nm阶设计矩阵的转置，m<n则无法计算。

y：m阶观测残差，y=v=l-HX 。

X：传出参数、待估计的n阶参数向量的增量。

Q：传出参数、nn协方差阵。

extern int lsq(const double *A, const double *y, int n, int m, double *x,
               double *Q)
{
    double *Ay;
    int info;
    
    if (m<n) return -1;
    Ay=mat(n,1);
    matmul("NN",n,1,m,1.0,A,y,0.0,Ay); /* Ay=A*y */
    matmul("NT",n,n,m,1.0,A,A,0.0,Q);  /* Q=A*A' */
    if (!(info=matinv(Q,n))) matmul("NN",n,1,n,1.0,Q,Ay,0.0,x); /* x=Q^-1*Ay */
    free(Ay);
    return info;
}

filter_()、filter():卡尔曼滤波器

X：n阶参数向量，就是参数，不再是最小二乘中的参数的增量

P：nn阶协方差阵。X、P在filter()函数中既是传入参数，也是传出参数，迭代。

H：nm阶设计矩阵的转置

V：m阶新息向量

R：mm阶量测噪声协方差阵

Xp、Pp：预测参数和预测协方差

extern  int filter_(const double *x, const double *P, const double *H,
                   const double *v, const double *R, int n, int m,
                   double *xp, double *Pp)
{
    double *F=mat(n,m),*Q=mat(m,m),*K=mat(n,m),*I=eye(n);
    int info;
    
    matcpy(Q,R,m,m);
    matcpy(xp,x,n,1);
    matmul("NN",n,m,n,1.0,P,H,0.0,F);       /* Q=H'*P*H+R */
    matmul("TN",m,m,n,1.0,H,F,1.0,Q);
    if (!(info=matinv(Q,m))) {
        matmul("NN",n,m,m,1.0,F,Q,0.0,K);   /* K=P*H*Q^-1 */
        matmul("NN",n,1,m,1.0,K,v,1.0,xp);  /* xp=x+K*v */
        matmul("NT",n,n,m,-1.0,K,H,1.0,I);  /* Pp=(I-K*H')*P */
        matmul("NN",n,n,n,1.0,I,P,0.0,Pp);
    }
    free(F); free(Q); free(K); free(I);
    return info;
}

extern int filter(double *x, double *P, const double *H, const double *v,
                  const double *R, int n, int m)
{
    double *x_,*xp_,*P_,*Pp_,*H_;
    int i,j,k,info,*ix;
    
    ix=imat(n,1); for (i=k=0;i<n;i++) if (x[i]!=0.0&&P[i+i*n]>0.0) ix[k++]=i;
    x_=mat(k,1); xp_=mat(k,1); P_=mat(k,k); Pp_=mat(k,k); H_=mat(k,m);
    for (i=0;i<k;i++) {
        x_[i]=x[ix[i]];
        for (j=0;j<k;j++) P_[i+j*k]=P[ix[i]+ix[j]*n];
        for (j=0;j<m;j++) H_[i+j*k]=H[ix[i]+j*n];
    }
    info=filter_(x_,P_,H_,v,R,k,m,xp_,Pp_);
    for (i=0;i<k;i++) {
        x[ix[i]]=xp_[i];
        for (j=0;j<k;j++) P[ix[i]+ix[j]*n]=Pp_[i+j*k];
    }
    free(ix); free(x_); free(xp_); free(P_); free(Pp_); free(H_);
    return info;
}

smoother()：Kalman smoother 、结合前向滤波和后向滤波的卡尔曼滤波平滑

xf、Qf：前向滤波n阶结果和nn阶协方差

xb、Qb：后向滤波n阶结果和nn阶协方差

xs、Qs：平滑n阶结果和nn阶协方差

extern int smoother(const double *xf, const double *Qf, const double *xb,
                    const double *Qb, int n, double *xs, double *Qs)
{
    double *invQf=mat(n,n),*invQb=mat(n,n),*xx=mat(n,1);
    int i,info=-1;
    
    matcpy(invQf,Qf,n,n);
    matcpy(invQb,Qb,n,n);
    if (!matinv(invQf,n)&&!matinv(invQb,n)) {
        for (i=0;i<n*n;i++) Qs[i]=invQf[i]+invQb[i];
        if (!(info=matinv(Qs,n))) {
            matmul("NN",n,1,n,1.0,invQf,xf,0.0,xx);
            matmul("NN",n,1,n,1.0,invQb,xb,1.0,xx);
            matmul("NN",n,1,n,1.0,Qs,xx,0.0,xs);
        }
    }
    free(invQf); free(invQb); free(xx);
    return info;
}