开源GNSS软件RTKLIB学习总结（一）MarkDown、VS配置RTKLIB、矩阵、最小二乘和Kalman滤波基本函数

本文从昨晚写到今天上午，比较粗糙，且因为本人对GNSS和RTKLIB掌握的还很少，可能错误比较多，各位可以批判着看。

一、前言

本人是导航工程大二的学生，计划读研读博当老师，可能要和GNSS打几十年的交道。从大一暑假开始断断续续自学GNSS，看过B站上武大的平差课，哈工程的GNSS原理课，矿大王坚的《GNSS程序设计讲义》、简单看了一遍吴云主编的《最优估计基础》和安理工的《GNSS导航定位原理》、目前在看杨元喜的《自适应动态导航定位》，昨天刚买了本杨元喜主编的《GNSS网络RTK技术原理与工程应用》，打算暑假自己给我的项目写RTK实时定位算法。

大一就听老师说学习RTKLIB的源码，先看懂别人的代码，才能自己写代码，但一直没有系统的开始学。现在掌握了一些理论知识，可以开始学了。

我计划以博客的形式总结学习的过程，既是为了对学习的巩固，也为了可以分享学习到的知识。计划每周一篇，下周我打算写卫星导航时间系统、坐标系统、和Rinex文件，相关的知识、函数。

二、MarkDown的学习总结

本人第一次使用Markdown，先看的”我是管小亮“这篇博客:Markdown入门指南【我为什么要推荐你学习Markdown？】、使用Typora软件

1、MarkDown概述

Markdown 是一种可以使用普通文本编辑器编写的轻量级标记语言，通过简单的标记语法，它可以使普通文本内容具有一定的格式，可以导出 PPT、LaTex、HTML、Word、PDF、Epub、JPG 等多种格式的文档，文档后缀为 .md。

Markdown 的语法可以说是非常简单，简单到常用的语法操作不超过十个，相较于同为标记语言的 HTML 来说，Markdown 被称为小儿科也不为过，学习成本极低，上手率超高，并且一旦熟悉这种语法规则，稍加练习，就会有一劳永逸的效果，回报率极高。

听说 Markdown 现在是互联网上最流行的写作语言！！！无论是博客，笔记还是论文等等，都是基于它写的！！！

Markdown 再牛，也不过是个工具，最终还是要看你的文字和图片等内容，所以很多的 Markdown 工具集成了不同的功能（基础功能大致相同），例如流程图与时序图，复杂表格与复杂公式的呈现，虽然好像挺有用，但是还是与简单书写的初衷不符，用其他专业绘图工具不香嘛？不要舍本逐末！！！

2、MarkDown基础语法

标题：#文字，一个#代表一级标题，共六级

文本：*斜体*、_斜体_、**粗体**、__粗体__、 ***粗斜体***、___粗斜体___

列表：无序：*、+、-+空格+文字

有序：：：数字+.+空格

控制列表层级：在 -前Tab

下划线：<u>文字</u>

分割线：三个#或-，同时在上面空一行

删除线：文字前后各加两个波浪线~~

引用：>文字

表格：|分隔不同单元格，-分隔表头和其它行

:-：左对齐、 -:：右对齐、:-:：居中对齐

图片：！[描述](地址)

脚注：[文字](接受 "脚注名称")

链接：[链接名称](链接地址)

代码块：行内：前后加反引号`,（英文输入法ESC下面）

块：前后一行各三个反引号，第一行上写语言名

3.快捷键

撤销：Ctrl + Z

重做：Ctrl + Y

加粗：Ctrl + B

斜体：Ctrl + I

标题：Ctrl + Shift + H

无序列表：Ctrl + Shift + U

有序列表：Ctrl + Shift + O

检查列表：Ctrl + Shift + C

插入代码：Ctrl + Shift + K

插入链接：Ctrl + Shift + L

插入图片：Ctrl + Shift + G

查找：Ctrl + F

替换：Ctrl + G

三、将 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;”。

四、前置知识

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变量apinti;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 apputchar('\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 NOCALLOCif ((p=mat(n,m))) for (n=n*m-1;n>=0;n--) p[n]=0.0;
#elseif (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);}
#endifreturn 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;
}