本科生如何入门GNSS算法(四)- rtklib单点定位解算源码说明

目录

单点定位

LSQ

设计矩阵

残差阵

后记

欢迎关注个人公众号


单点定位

单点定位的过程就是最小二乘平差的过程,所涉及到的卫星位置计算/大气误差计算等都是为了计算残差。关于最小二乘,知乎有大量优秀的回答,比如 最小二乘法的本质是什么? - 马同学的回答 - 知乎 。 使用 - 知乎 最小二乘法的本质是什么? - Mr.Bo的回答 回答中的公式来说明整个单点定位的流程。最小二乘步骤:

​对应的在blog中的 定位方程解算和定位精度 相同,即我们要构造一个方程组: A*x=b或者G*x=b

一般我们将A或者G称为设计矩阵(design matrix),或者blog中的雅可比矩阵。而b可以称为先验残差阵,定位的过程就是为了计算这两个矩阵。

​先让我们打开rtklib中最小二乘的函数,在rtkcnm.c的1247行。即如下图,其中的A就是设计矩阵,y就是残差,即观测值残差。x就是最小二乘估计得到的状态量 x=(A*A')^-1*A*y,也就是我们的定位结果。在这里因为,我们依赖概略位置进行了残差求解,所以估计得到的状态量是相对概略位置的修正量。


LSQ

lsq则在pntpos.c 388行estpos函数中进行调用;estpos函数则被pntpos.c中的pntpos函数调用。

pntpos函数是整个单点定位的入口,来仔细看一下函数说明和函数变量。

如下图中,obs/nav/opt分别为观测值/星历/配置,而且都是不可修改的const变量,都是从文件读取后存储到指定结构体中,这些都是函数的输入量。具体的如何把文件存储为结构体数据,可自行学习;而sol/azel/ssat等都是可修改的变量,传入的都是指针,表明这些是输出量。

即这个函数实现了输入原始数据,并输出定位结果的功能。如果对一些变量名字有疑问,可以参考该 网页中的介绍 estpos ,有详细的调用关系。


设计矩阵

初步的调用关系已经明白了,那我们要看一下是如何计算得到的设计矩阵和残差阵。

从各个blog介绍以及学习中,已知道设计矩阵计算需要概略位置和每个卫星的位置。对于初次定位我们一般将概略位置设置为地球质心,也就是我们参考坐标系的零点(0,0,0);而卫星位置计算,则在satposs函数中。

对于我们此次的目标是单点定位,你只需要搞清楚这个函数的输入以及输出即可。这也是模块化编程或者面向对象编程的目的之一,不需要了解具体实现,只需要你输入什么数据,输出的数据是如何存储的。具体函数解释可自行理解,或者参考我上面贴的网址介绍。

​设计矩阵与残差阵的计算都在pntpos.c 268行rescode()函数中,可以看到该函数前面都是const或者值传递,都是输入量,而H就是输出的设计矩阵,v就是输出的残差阵,var则是观测值的权重。具体的输出量是如何计算的,自行对照公式理解。

​需要注意的是,设计矩阵的具体计算其实是在geodist函数中,该函数被rescode调用。该函数的中的输出 e (line-of-sight vector(ecef))就是各个卫星的设计矩阵,也就是视线方向ecef框架下(坐标系)三方向的系数。


残差阵

以下语句就是计算伪距残差,其中r是卫地距,由接收机概略位置和卫星位置计算得到;dtr是接收机钟差,可给0或者上一次迭代结果;dts是卫星钟差,由广播星历计算得到;dion是电离层误差,由广播星历中的电离层模型或者经验模型得到;dtrp是对流层误差,由经验模型计算得到;

​前面文章已经介绍,如果电离层或者对流层都是已有模型计算,可完全不关注计算过程。

至此,设计矩阵和残差阵均已得到,输入到lsq就可以估计得到结果。中间有观测值的方差计算我没有讲,这一部分可以自行研究,简单一点就是依据高度角或者信噪比计算每个观测值的方差信息,可简单看一些教材或者论文,自行研究。

另外需要注意的是,计算过程都需要迭代,比如初值我们一般给地球质心,经过多次迭代即可收敛到正确结果。此处rtklib直接给了一个10次的迭代循环,这样设置不太科学;我们完全可以从状态量是否收敛来终止循环,具体方法可自行研究。


后记

至此,单点定位的入门已初步完结。写的过程中我发现网上其实有很多的rtklib相关的资料,比我读研究生多的多的多。当然大家的侧重点都不一致,对于初学者可能会有摸不到门路的感觉。

还是要多看多研究,功夫到了,知识点就可以串起来了。

下一节我会写一些单点的优化方向以及下一步的计划,其实单点定位在工程应用中还是非常重要,尤其是在复杂场景下,单点做不好,没有好的初始位置,设计矩阵都存在不可忽略的误差,后续的rtk更没有办法进行;甚至有些场景只能输出单点解,单点定位结果尤其重要。

欢迎关注个人公众号

个人公众号 GNSS和自动驾驶,会持续更新GNSS的基础教程/进阶教程/GNSS在自动驾驶中的应用/自动驾驶技术进展等。

其他相关文章

测绘工程本科生如何入门GNSS算法 - 引言_梧桐Fighting的博客-CSDN博客 

本科生如何入门GNSS算法 (一)- 在Visual Studio 2017中跑通RTKLIB开源代码_梧桐Fighting的博客-CSDN博客  
本科生如何入门GNSS算法(二)- rtklib定位解算过程中的GNSS数据格式以及基本概念_梧桐Fighting的博客-CSDN博客
本科生如何入门GNSS算法(三)- rtklib单点定位解算基础知识储备.docx_梧桐Fighting的博客-CSDN博客
本科生如何入门GNSS算法(四)- rtklib单点定位解算源码说明_梧桐Fighting的博客-CSDN博客
 

  • 1
    点赞
  • 24
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
在Matlab中,可以使用最小二乘法进行GNSS单点定位解算。以下是一个简单的示例代码: 假设我们有颗卫星,它们的位置和伪距信息如下: ``` % 卫星位置 x1 = 15600; y1 = 7540; z1 = 20140; x2 = 18760; y2 = 2750; z2 = 18610; x3 = 17610; y3 = 14630; z3 = 13480; x4 = 19170; y4 = 610; z4 = 18390; % 伪距信息 p1 = 20700.86; p2 = 25549.05; p3 = 21688.08; p4 = 23672.17; ``` 接下来,我们需要构造误差方程,并使用最小二乘法求解。 ``` % 构造误差方程 A = [-2*x1+2*x2, -2*y1+2*y2, -2*z1+2*z2, 1; -2*x1+2*x3, -2*y1+2*y3, -2*z1+2*z3, 1; -2*x1+2*x4, -2*y1+2*y4, -2*z1+2*z4, 1; -2*x2+2*x3, -2*y2+2*y3, -2*z2+2*z3, 1; -2*x2+2*x4, -2*y2+2*y4, -2*z2+2*z4, 1; -2*x3+2*x4, -2*y3+2*y4, -2*z3+2*z4, 1]; b = [p1^2-x1^2-y1^2-z1^2-p2^2+x2^2+y2^2+z2^2; p1^2-x1^2-y1^2-z1^2-p3^2+x3^2+y3^2+z3^2; p1^2-x1^2-y1^2-z1^2-p4^2+x4^2+y4^2+z4^2; p2^2-x2^2-y2^2-z2^2-p3^2+x3^2+y3^2+z3^2; p2^2-x2^2-y2^2-z2^2-p4^2+x4^2+y4^2+z4^2; p3^2-x3^2-y3^2-z3^2-p4^2+x4^2+y4^2+z4^2]; % 使用最小二乘法求解 result = pinv(A)*b; % 输出结果 x = result(1) y = result(2) z = result(3) t = result(4) / 299792458 % 时间偏差 ``` 在这个示例中,我们使用了矩阵求逆的方法来求解最小二乘问题,其中 `pinv` 函数用于求解矩阵的伪逆。最终输出的 `x`、`y`、`z` 分别表示接收机的位置,`t` 表示接收机时间偏差。需要注意的是,这个代码只是一个简单的示例,实际应用中可能需要进行更多的误差分析和校正。

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值