relpos函数解读

relpos函数解读

在配合网络的一些解读和一些自己的理解，将RTKLIB2.3.2的relpos函数进行解读
流程：
zdres—ddres—ekf求出浮点解----lambda固定解-----固定解带回ddres

1.源程序

略

2.大概框架

1. 计算卫星位置、速度
2. 计算基站非差残差项
3. 选择共视卫星
4. 计算非差残差项和双差
5. 使用kalman滤波计算浮点解
6. 重新进行双差计算
7. 计算固定解

3.各个步骤的关键算法如下：

1.计算卫星位置速度

  对应代码： relpos--satposs
  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)
 input:观测值（obs）,导航数据（nav）
 output:none
 function：利用obs和nav计算卫星位置，速度和时钟并保存在rs，dts,var等参数中
note：
rs [（0：2）+ i * 6] = obs [i]卫星位置{x，y，z}（m）
* rs [（3：5）+ i * 6] = obs [i]卫星速度{vx，vy，vz}（m / s）
* dts [（0：1）+ i * 2] = obs [i]卫星时钟{bias，drift}（s | s / s）
* var [i] = obs [i卫星位置和时钟误差方差（m ^ 2）
* svh [i] = obs [i]卫星健康状态标记
 *如果没有导航数据，则将0设置为rs []，dts []，var []和svh []
*卫星位置和时钟是信号传输时的值
*卫星位置参考天线相位中心
*卫星时钟不包括代码偏差校正（tgd或bgd）
*始终需要任何伪距和广播星历来获取信号传输时间

2. 计算基站非差残差项

 对应代码： relpos--zdres
/* undifferenced phase/code residuals ----------------------------------------*/
static int zdres(int base, const obsd_t *obs, int n, const double *rs, const double *dts, const int *svh, const nav_t *nav, const double *rr, const prcopt_t *opt, int index, double *y, double *e, double *azel)
输入：obs,n
输出：y(+nu*nf*2):基准站的无差观测向量(存了无差相位和码的残差)
		e(+nu*3):基准站的视线向量（卫星坐标-基准站坐标）
		azel(+nu*2)：基准站的方位角仰角向量
功能:计算基站的无差残差（预测值和真实值之差）
for（遍历基站卫星）
{
geodist：计算几何距离
satazel：计算高度角和仰角
satexclude：查看卫星是否被排除
r+=-CLIGHT*dts[i*2]：计算卫星钟差
tropmodel：通过模型计算大气层延迟
tropmapf：通过NMF计算对流层映射函数
antmodel：接收机天线相位中心校正
zdres_sat：计算卫星非差相位/码残差（观察值-伪距）
}

3. 选择基准站与移动站共视卫星

   对应代码： relpos--selsat
static void udstate(rtk_t *rtk, const obsd_t *obs, const int *sat,
                    const int *iu, const int *ir, int ns, const nav_t *nav)
将基准站卫星，将通过“高度角阈值”的卫星，作为共视卫星；

找到共同观察到的卫星，然后顺便把他们在数组中的索引 index（iu,ir）记录下来

4.kalman 时间/状态更新 udstate

static void udstate(rtk_t *rtk, const obsd_t *obs, const int *sat,
                    const int *iu, const int *ir, int ns, const nav_t *nav)
主要更新了：
udpos ：kalman中状态更新方程、状态协方差方程更新；
udion ： 电离层状态、协方差更新；
udtrop ： 对流层状态、协方差更新；
udrcvbias：接收器的时间更新
udbias ：更新单差模糊度状态、协方差；
存放到了rtk参数中

5. 1计算移动站非差残差项

if (!zdres(0,obs,nu,rs,dts,svh,nav,xp,opt,0,y,e,azel))

此函数与上边计算基站的无差观测为同一个函数，只不过总体说来有两点不同：
站点坐标使用的是$xp$,这个是kalman滤波时间更新中得到的新坐标，
所有的输出，如$y,e,azel$，起始存储位置为0，而上边的基准站无差观测时为$nu\ast size$

5.2 计算双差

 if ((nv=ddres(rtk,nav,dt,xp,Pp,sat,y,e,azel,iu,ir,ns,v,H,R,vflg))<1)

上边函数最重要的三个输出$v,H,R$都是卡尔曼滤波量测更新所必须的。
ekf的量测向量是双差载波和双差伪距，这里为什么观测变成双差残差了呢？

残差/v=【参考星（移动站非差残差 - 基准站非差残差）- 非参考星（移动站非差残差 - 基准站非差残差）】-【双差模糊度】
H = -参考星站星单位矢量 + 非参考星站星单位矢量；移动站下
H 参考星对应的列为1 ，非参考星对应的列为-1
保存伪距、载波残差；对新息阈值检测，卫星相关情况、相关值保存到
ssat_t中 resc、resp、rejc、vsat

5.3 kalman滤波状态量更新

if ((info=filter(xp,Pp,H,v,R,rtk->nx,nv)))

在这里，其实浮点解已经出来了：
$$\begin{gathered} 状态更新公式: \\ Q=H^{\prime }\ast P\ast H+R \\ K=P\ast H\ast Q^{-} \\ xp=x+K\ast v \\ Pp=(I-K\ast H^{\prime })\ast P \end{gathered}$$
kalman滤波状态更新主要内容有3个，即a 状态转移矩阵的确定 b 初值确定 c 状态向量确定

rtklib中状态向量选择移动站的位置速度 以及单差整周模糊度（x,v,N1,N2,N3…Nm）
初值选择：坐标速度为第一步中计算的概略坐标，整周模糊度初值则使用伪距载波相位组合确定uppos upbias

状态转移矩阵的选择：对于位置速度可分为静态模型和动态模型 高动态模型，具体模型与单点定位kalman滤波模型一致uppos
对于整周模糊度，未发生周跳则整周模糊度不变，发生周跳时，则重置滤波器，重新使用伪距载波组合计算概略整周模糊度为初值进行滤波

计算量测值与预测值差v和观测矩阵H(ddres)，并进行卡尔曼滤波的量测更更新filter,此步得到浮点解

周跳的判断：rtklib使用电离层残差法进行周跳的探测
对应代码udbias-detslp_gf_L1L2/L1L5

5. 重新进行双差计算

重新计算双差ddres

6. 计算固定解

使用LAMBDA算法计算固定解
对应代码： resamb_LAMBDA

下边才是重头戏，单差整周模糊度的固定，只有固定了模糊度，才能得到固定解，也就是真正意义上的高精度解。

rtklib中的整周固定解通过lambda算法求解，整周固定解的求解过程其数学本质是一个混合整数最小二乘问题，即待求解的变量中有一部分是浮点数，一部分是整数（单差模糊度）。如果不考虑计算效率，那么我们通过估计误差设置一个搜索区间，带入所有模糊度的组合，取残差最小的那一组就是我们要求解的整周模糊度。但是，这个计算量无疑是一个天文数字！想象一下，我们有14*3=42个单差模糊度要求解，每个浮点周围我们取10个待算整数，那么组合数量就是$10^{42}$，即便是后处理应用，这个量级也挺吓人的。
所以，lambda算法本质上是解决计算速度的问题，relpos的核心算法或者之一，
这里我们需要知道，通过这个函数之后固定解rtk->xa得到了

参考文献

链接：
https://blog.csdn.net/jack909633117/article/details/87997123
https://blog.csdn.net/iceboy314159/article/details/105422468
https://blog.csdn.net/wuwuku123/article/details/107410950
https://blog.csdn.net/weixin_42918498/article/details/106747386
https://www.cnblogs.com/taqikema/p/8819798.html