运行和调试
main–>postpos–>execses_r–>execses_b–>execses
execses函数进入rtklib的实际定位流程，前面的函数多是读入文件和设置的一些有效性判断？？。
execses
-
读入观测数据和导航电文：readobsnav
-
读入天线相位中心误差：setpcv–>searchpcv（通过卫星号和时间（接收机只通过时间）在pcvs_t结构体pcvss（卫星）、pcvsr（接收机）中查找对应的pcv（天线参数））
-
读入海水负荷潮汐误差模型参数：readotl
-
计算基站位置：antpos（可选postype：1-使用单点定位平均值作为基站坐标，2- 从基站pos文件里读，3-从rinex文件头读取，似乎与rtklib软件界面有所不同，软件界面：平均值、在界面中输入给定的基站坐标），一般情况下有坐标信息就直接给定基站坐标，没有则使用平均值。
-
open solution statistics：没进去
-
写输出文件头outhead：（一般是.pos格式）的文件头（定位模式、误差模型选择，输入文件，基站坐标等信息）
-
判断滤波是forward、backward、combined后进入procpos函数–>每个历元调用rtkpos函数
单点定位
相对定位
rtkpos
nr基站卫星数，nu流动站卫星数
- 单点定位计算流动站位置
- 判断：单点定位（return）、精密单点定位（进入pppos解算）、相对定位模式都需要进入relpos
- 相对定位relpos
relpos相对定位算法在此实现
zdres非差残差
输出y矩阵：先存载波后存伪距【卫星1：载波（nf个频率）；伪距nf个频率
卫星2：……】
时间预测ustate
udstate函数解析，站间单差
udpos
与模糊度（即相位偏移量）无关的待估计状态量的时间预测：
-
若dynamics模式开启，则状态量中有速度和加速度，需要根据时间间隔和转移矩阵进行时间预测
-
dynamics模式关闭，位置直接等于单点定位结果，没有速度和加速度时，x(k+1)=x(k)，相当于不需要时间预测，每个历元的位置初值等于单点定位结果
udbias
模糊度状态量（站间单差模糊度）的时间预测
整周模糊度解算resamb_LAMBDA
固定模糊度为什么要用双差：单差模糊度里面还有很多误差没有消除，不一定具有整数特性，尤其是卫星钟差的影响，双差消除了大部分误差，整数特性更好。
ddmat
ddmat求得D矩阵与X相乘得到的双差模糊度与双差残差不对应，只是第一颗无失锁&无周跳&高度角>最小高度角的卫星和其他卫星相减，并不是参考星的模糊度和非参考星相减。
/* single to double-difference transformation matrix (D') */
/*卡尔曼滤波解算时，x状态量中是单差模糊度，y=D'*x，通过D矩阵转换为双差模糊度*/
D=zeros(nx,nx);
if ((nb=ddmat(rtk,D))<=0) {
errmsg(rtk,"no valid double-difference\n");
free(D);
return 0;
}
ny=na+nb; y=mat(ny,1); Qy=mat(ny,ny); DP=mat(ny,nx);
b=mat(nb,2); db=mat(nb,1); Qb=mat(nb,nb); Qab=mat(na,nb); QQ=mat(na,nb);
/* transform single to double-differenced phase-bias (y=D'*x, Qy=D'*P*D) 模糊度;单差->双差*/
matmul("TN",ny, 1,nx,1.0,D ,rtk->x,0.0,y );//y=D'*x
matmul("TN",ny,nx,nx,1.0,D ,rtk->P,0.0,DP);//DP=D'*P
matmul("NN",ny,ny,nx,1.0,DP,D ,0.0,Qy);//Qy=DP*D
/* phase-bias covariance (Qb) and real-parameters to bias covariance (Qab) */
for (i=0;i<nb;i++) for (j=0;j<nb;j++) Qb [i+j*nb]=Qy[na+i+(na+j)*ny];
for (i=0;i<na;i++) for (j=0;j<nb;j++) Qab[i+j*na]=Qy[ i+(na+j)*ny];
lambda
LD：对Q进行Choleshy分解，Q=L’DL
reduction：整数变换降相关（gauss）&D矩阵降序排列（perm）
reference :
[1] P.J.G.Teunissen, The least-square ambiguity decorrelation adjustment:
a method for fast GPS ambiguity estimation, J.Geodesy, Vol.70, 65-82,
1995
[2] X.-W.Chang, X.Yang, T.Zhou, MLAMBDA: A modified LAMBDA method for
integer least-squares estimation, J.Geodesy, Vol.79, 552-565, 2005
integer least-square estimation. reduction is performed by lambda (ref.[1]),
and search by mlambda (ref.[2])
/* lambda reduction (z=Z'*a, Qz=Z'*Q*Z=L'*diag(D)*L) (ref.[1]) -整数变换降相关&D矩阵降序排列--------------*/
static void reduction(int n, double *L, double *D, double *Z)
{
int i,j,k;
double del;
j=n-2; k=n-2;//n-2，n是模糊度数量，换到矩阵索引要-1，第j列的del要和j+1列比（j+1<=n-1）,所以n-2
while (j>=0) {
if (j<=k) for (i=j+1;i<n;i++) gauss(n,L,Z,i,j);//整数变换降相关
//D矩阵降序排列,调序变换在prem中实现
del=D[j]+L[j+1+j*n]*L[j+1+j*n]*D[j+1];/*《GPS测量与数据处理》 （5-72）*/
if (del+1E-6<D[j+1]) { /* compared considering numerical error （5-72）满足此条件需要调序变换*/
perm(n,L,D,j,del,Z);/* 调序变换 */
k=j; j=n-2;
}
else j--;
}
}
guass：整数变换降相关
/* integer gauss transformation -----整数高斯变换---------------------------------*/
static void gauss(int n, double *L, double *Z, int i, int j)
{
int k,mu;
if ((mu=(int)ROUND(L[i+j*n]))!=0) {
for (k=i;k<n;k++) L[k+n*j]-=(double)mu*L[k+i*n];
for (k=0;k<n;k++) Z[k+n*j]-=(double)mu*Z[k+i*n];
}
}
perm：D矩阵降序排列
/* permutations ---排序-------------------------------------------------------*/
static void perm(int n, double *L, double *D, int j, double del, double *Z)
{
int k;
double eta,lam,a0,a1;
eta=D[j]/del;
lam=D[j+1]*L[j+1+j*n]/del;
D[j]=eta*D[j+1]; D[j+1]=del;
for (k=0;k<=j-1;k++) {
a0=L[j+k*n]; a1=L[j+1+k*n];
L[j+k*n]=-L[j+1+j*n]*a0+a1;
L[j+1+k*n]=eta*a0+lam*a1;
}
L[j+1+j*n]=lam;
for (k=j+2;k<n;k++) SWAP(L[k+j*n],L[k+(j+1)*n]);
for (k=0;k<n;k++) SWAP(Z[k+j*n],Z[k+(j+1)*n]);
}
search：在阈值范围内找整周模糊度
:z变换后的整周模糊度浮点解
1.第一次搜索，N_n=int(N_n浮点解)从最后一个模糊度往前推前面的模糊度的范围，（5-80）
2.第一个搜索序列，一律取范围内的中间值，即int(Z_i)
如果后面的搜索N_k情况下有小于maxdist的，保存到s和zn数组，maxdist=该搜索结果下的newdist，继续搜索k+1。
似乎除了N_0搜索了int（Z_0）附近的三个整数之外，其他模糊度只搜索了Z_i左右的两个整数，f是一个凹函数吗，会不会陷入局部最优。
更新rtk结果
if (!(info=lambda(nb,2,y+na,Qb,b,s))) {
trace(4,"N(1)="); tracemat(4,b ,1,nb,10,3);
trace(4,"N(2)="); tracemat(4,b+nb,1,nb,10,3);
rtk->sol.ratio=s[0]>0?(float)(s[1]/s[0]):0.0f;//ratio赋值
if (rtk->sol.ratio>999.9) rtk->sol.ratio=999.9f;
/* validation by popular ratio-test */
/*
如果最优解和次优解的比值大于阈值，则利用公式（E.7.19），根据lambda算法得到的双差整周模糊度，
计算除了单差相位偏移之外的所有状态量（例如：位置+速度+加速度+电离层+对流层……），存入rtk->xa中。
这一步实际就是利用lambda算法得到的整数的整周模糊度对其他状态量进行修正。
*/
if (s[0]<=0.0||s[1]/s[0]>=opt->thresar[0]) {
/* transform float to fixed solution (xa=xa-Qab*Qb\(b0-b)) 根据固定解更新相位偏差无关的状态量，存入rtk->xa*/
for (i=0;i<na;i++) {
rtk->xa[i]=rtk->x[i];
for (j=0;j<na;j++) rtk->Pa[i+j*na]=rtk->P[i+j*nx];
}
for (i=0;i<nb;i++) {
bias[i]=b[i];
y[na+i]-=b[i];/* y矩阵na索引后的双差浮点解替换成(b0-b)，即减去固定解*/
}
if (!matinv(Qb,nb)) {/* Qb=Qb^(-1) */
matmul("NN",nb,1,nb, 1.0,Qb ,y+na,0.0,db);/* db=Qb^(-1)*(b0-b) */
matmul("NN",na,1,nb,-1.0,Qab,db ,1.0,rtk->xa); /*xa = xa -Qab*db= xa - Qab * Qb\(b0 - b)*/
/* covariance of fixed solution (Qa=Qa-Qab*Qb^-1*Qab')根据固定解更新方差阵，存入rtk->Pa */
matmul("NN",na,nb,nb, 1.0,Qab,Qb ,0.0,QQ);
matmul("NT",na,na,nb,-1.0,QQ ,Qab,1.0,rtk->Pa);
trace(3,"resamb : validation ok (nb=%d ratio=%.2f s=%.2f/%.2f)\n",
nb,s[0]==0.0?0.0:s[1]/s[0],s[0],s[1]);
/* restore single-differenced ambiguity 调用restamb函数，
利用lambda算法计算得到的双差整周模糊度，重新计算单差相位偏移，
并存入xa中，同时将得到的其他状态量也存入xa中。*/
restamb(rtk,bias,nb,xa);
}
else nb=0;
}
else { /* validation failed */
errmsg(rtk,"ambiguity validation failed (nb=%d ratio=%.2f s=%.2f/%.2f)\n",
nb,s[1]/s[0],s[0],s[1]);
nb=0;
}
}
}
else nb=0;
}
else { /* validation failed */
errmsg(rtk,"ambiguity validation failed (nb=%d ratio=%.2f s=%.2f/%.2f)\n",
nb,s[1]/s[0],s[0],s[1]);
nb=0;
}
}
扫一扫
2448
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
