上海门户网站开发,美颜秘籍网站建设,wordpress 视频播放,网站上线 流程1 ddres( )参数介绍
rtklib中进行的单频解算 双差观测值#xff0c;单差的模糊度 单频点双差 DD (double-differenced) phase/code residuals ------------------------------ x 模糊度 P 方差-协方差阵 sat 共识卫星列表 ns 共识卫星数量 y…1 ddres( )参数介绍
rtklib中进行的单频解算 双差观测值单差的模糊度 单频点双差 DD (double-differenced) phase/code residuals ------------------------------ x 模糊度 P 方差-协方差阵 sat 共识卫星列表 ns 共识卫星数量 y 残差阵,数量ns*nf*2(每颗卫星的每个频率×2是伪距和载波) e dx dy dz设计矩阵,数量ns*3 freq 频率,数量ns*nf 每个频点的伪距和载波频率一样而且伪距也用不到频率是载波求波长时用的 v 完整残差阵(包括了模糊度) H dx dy dz N1 N2 N2等模糊度 设计矩阵 R 方差-协方差阵 static int ddres(rtk_t *rtk, const nav_t *nav, double dt, const double *x,const double *P, const int *sat, double *y, double *e,double *azel, double *freq, const int *iu, const int *ir,int ns, double *v, double *H, double *R, int *vflg)
2 代码解读
计算基线的长度计算基准站和流动站方差时会用到 blbaseline(x,rtk-rb,dr);ecef2pos(x,posu); ecef2pos(rtk-rb,posr);
Ri是参考卫星方差RJ是非参考卫星方差。用于后面计算双差后的观测值方差--协方差阵
tropu、tropr等是对流层、电离层参数在此我们不作讨论 Rimat(ns*nf*22,1); Rjmat(ns*nf*22,1); immat(ns,1);tropumat(ns,1); troprmat(ns,1); dtdxumat(ns,3); dtdxrmat(ns,3); 初始化伪距和载波相位残差为0每个卫星每个频点的用于以后记录双差残差 rtk-ssat[i]访问到了每颗卫星 double resp[NFREQ]; residuals of pseudorange (m) double resc[NFREQ]; residuals of carrier-phase (m) MAXSAT最大卫星数不管观没观测到都初始化。这也造成了空间上的浪费 NFREQ最大频率 for (i0;iMAXSAT;i) for (j0;jNFREQ;j) {rtk-ssat[i].resp[j]rtk-ssat[i].resc[j]0.0;}
对流层和电离层处理暂不考虑 for (i0;ins;i) {if (opt-ionooptIONOOPT_EST) {im[i](ionmapf(posu,azeliu[i]*2)ionmapf(posr,azelir[i]*2))/2.0;}if (opt-tropoptTROPOPT_EST) {tropu[i]prectrop(rtk-sol.time,posu,0,azeliu[i]*2,opt,x,dtdxui*3);tropr[i]prectrop(rtk-sol.time,posr,1,azelir[i]*2,opt,x,dtdxri*3);}}
之后就进入一个大循环最外层是先遍历每个卫星系统
for (m0;m6;m) /* m0:GPS/SBS,1:GLO,2:GAL,3:BDS,4:QZS,5:IRN */
二重循环根据所选择的定位模式进行相应的频点遍历。nf*2是因为每个系统的每个频点(*2 包括伪距和载波) 只利用伪距观测值nf---2*nf因为伪距残差是在y[]后半 是载波(动态定位模式及以后)那我们从0--2*nf因为同时使用伪距和载波残差 for (fopt-modePMODE_DGPS?0:nf;fnf*2;f) ns---共识卫星数量 这个就是找参考卫星的在这里是把高度角最高的那个作为参考卫星。我们也可以改写代码设置其他的规则。参考卫星下标 i for (i-1,j0;jns;j) {sysirtk-ssat[sat[j]-1].sys;if (!test_sys(sysi,m)) continue;//当前共识卫星的系统和当前遍历的系统是否一样if (!validobs(iu[j],ir[j],f,nf,y)) continue; //评价是否可以if (i0||azel[1iu[j]*2]azel[1iu[i]*2]) ij;//选择排序的那种,i最终是参考卫星的index} 又是一个大循环遍历所有的共识卫星这次是要计算双差残差和设计矩阵。组双差的话两台接收机是要都观测到相同的卫星的所以需要共识卫星处理。 以下代码均是在此循环中执行
for (j0;jns;j) {}
两个不同的卫星所以i≠j freq[f%nfiu[i]*nf]取频率 f区间为[0,2*nf]所以取余% freqifreqj会在电离层、对流层用到 if (ij) continue;sysirtk-ssat[sat[i]-1].sys;//参考卫星系统sysjrtk-ssat[sat[j]-1].sys;//非参考卫星系统freqifreq[f%nfiu[i]*nf];freqjfreq[f%nfiu[j]*nf];if (!test_sys(sysj,m)) continue;//判断和当前遍历系统一样否if (!validobs(iu[j],ir[j],f,nf,y)) continue;//检验数据有效yy[]
设计矩阵初始化nx为待估量的个数 if (H) {HiHnv*rtk-nx;for (k0;krtk-nx;k) Hi[k]0.0;} 求双差残差 (y[fiu[i]*nf*2]-y[fir[i]*nf*2]) 共识卫星站间单差是参考卫星i y[fiu[j]*nf*2]-y[fir[j]*nf*2] 共识卫星站间单差是参考卫星j 两个一减就是双差残差了 v[nv](y[fiu[i]*nf*2]-y[fir[i]*nf*2])-(y[fiu[j]*nf*2]-y[fir[j]*nf*2]); 设计矩阵e[ ]是dx dy dz的系数已在站间单差函数中求得我们只考虑单频e[ ]中只有单频的信息 e[]是dx dy dz设计矩阵大小流动站观测到的卫星基准站观测到的卫星。 因为在relpos()函数中站间单差[先是算基准站后是流动站]设计矩阵都填到e里面了。所以大小如此 设计矩阵非参考卫星-参考卫星 每次循环都会产生一组双差dx dy dz的系数。但最最外面还有一个频率循环那么不同频的dx dy dz的系数是一样的吗这我还不知道
这两个会是一样吗留到以后回答或者有谁知道
我的猜想是可能一样因为由计算系数的公式可知这个与卫星位置和接收机概略位置有关 if (H) {for (k0;k3;k) {Hi[k]-e[kiu[i]*3]e[kiu[j]*3];}} IB()不懂是什么。 哪个卫星哪个频点的模糊度位置索引
看上去像是索引但为什么x[ ] Hi[ ]都可以用。一个是设计矩阵位置一个是对应模糊度位置啥的 if (opt-ionoopt!IONOOPT_IFLC) {v[nv]-CLIGHT/freqi*x[IB(sat[i],f,opt)]-CLIGHT/freqj*x[IB(sat[j],f,opt)];if (H) {Hi[IB(sat[i],f,opt)] CLIGHT/freqi;Hi[IB(sat[j],f,opt)]-CLIGHT/freqj;}}else {//根据站星双差方差还要减去(模糊度的derT)见理论。而在y[]的计算中还没有考虑模糊度在这里给补上//x[] 模糊度v[nv]-x[IB(sat[i],f,opt)]-x[IB(sat[j],f,opt)];if (H) {Hi[IB(sat[i],f,opt)] 1.0; //IB() 哪个卫星哪个频点的模糊度位置索引不懂Hi[IB(sat[j],f,opt)]-1.0;}}} Hi[IB(sat[i],f,opt)] 1.0 这里是在补设计矩阵参考卫星系数为-1。而这里为1后续应该在某个计算中会改掉 v[ ]是计算好的双差残差将其存到ssat结构体中
之后检验残差是否在容许范围内 maxinno是设置的标准 if (fnf) rtk-ssat[sat[j]-1].resc[f ]v[nv];else rtk-ssat[sat[j]-1].resp[f-nf]v[nv];/* test innovation 检查残差是否符合标准 maxinno */if (opt-maxinno0.0fabs(v[nv])opt-maxinno) {if (fnf) {rtk-ssat[sat[i]-1].rejc[f];rtk-ssat[sat[j]-1].rejc[f];}errmsg(rtk,outlier rejected (sat%3d-%3d %s%d v%.3f)\n,sat[i],sat[j],fnf?L:P,f%nf1,v[nv]);continue;}
每颗卫星每个频率的观测值噪声(方差) Ri[nv]varerr(sat[i],sysi,azel[1iu[i]*2],bl,dt,f,opt);Rj[nv]varerr(sat[j],sysj,azel[1iu[j]*2],bl,dt,f,opt);
设置标志这颗卫星这个频点的伪距/载波观测值有效。标记一下 if (opt-modePMODE_DGPS) {if (fnf) rtk-ssat[sat[i]-1].vsat[f]rtk-ssat[sat[j]-1].vsat[f]1;}else {rtk-ssat[sat[i]-1].vsat[f-nf]rtk-ssat[sat[j]-1].vsat[f-nf]1;}
独特的编码由此码可知双差的具体信息。哪两个卫星作差载波还是伪距第几个频率
vflg[nv](sat[i]16)|(sat[j]8)|((fnf?0:1)4)|(f%nf); 至此共识卫星遍历结束 每个小频率有几颗卫星(残差有多少个)
nb[b]; 至此频点遍历结束
关于基线约束的 if (opt-modePMODE_MOVEBconstbl(rtk,x,P,v,H,Ri,Rj,nv)) {vflg[nv]34;nb[b];}
求双差后的协方差矩阵
ddcov(nb,b,Ri,Rj,nv,R);
3 ddcov( )
非对角线是Ri对角线是RiRj
static void ddcov(const int *nb, int n, const double *Ri, const double *Rj,int nv, double *R)
{int i,j,k0,b;trace(3,ddcov : n%d\n,n);//nv*nv 矩阵嘛for (i0;inv*nv;i) R[i]0.0;//初始化for (b0;bn;knb[b]) {for (i0;inb[b];i) for (j0;jnb[b];j) {R[ki(kj)*nv]Ri[ki](ij?Rj[ki]:0.0); //参考GNSS理论16.2讲。非对角线参考卫星方差}}trace(5,R\n); tracemat(5,R,nv,nv,8,6);
} 4 validobs()
评价也不太懂为什么可以不用管伪距 i--基准站共识卫星索引下标f--频点 0--nf 是载波 nf--2*nf 是伪距 当fnf载波频率残差 有载波就不管伪距了解这是f(0--2*nf)的情况随着循环进行伪距最终会遍历到(nf--2*nf)检查y[fi*nf*2]!0 检查基准站共识卫星的载波残差fnf||(y[f-nfi*nf*2]!0.0y[f-nfj*nf*2]!0.0) f-nf为真短路原则当fnf,伪距残差,载波也要有不然伪距不可用(注释)y[fi*nf*2] 每颗共识卫星伪距观测y[f-nfi*nf*2] 载波观测
*/
static int validobs(int i, int j, int f, int nf, double *y)
{/* if no phase observable, psudorange is also unusable */return y[fi*nf*2]!0.0y[fj*nf*2]!0.0(fnf||(y[f-nfi*nf*2]!0.0y[f-nfj*nf*2]!0.0));
}
5 全部代码
static int ddres(rtk_t *rtk, const nav_t *nav, double dt, const double *x,const double *P, const int *sat, double *y, double *e,double *azel, double *freq, const int *iu, const int *ir,int ns, double *v, double *H, double *R, int *vflg)
{prcopt_t *optrtk-opt;double bl,dr[3],posu[3],posr[3],didxi0.0,didxj0.0,*im;double *tropr,*tropu,*dtdxr,*dtdxu,*Ri,*Rj,freqi,freqj,*HiNULL;int i,j,k,m,f,nv0,nb[NFREQ*4*22]{0},b0,sysi,sysj,nfNF(opt);trace(3,ddres : dt%.1f nx%d ns%d\n,dt,rtk-nx,ns);//计算基线长度blbaseline(x,rtk-rb,dr);ecef2pos(x,posu); ecef2pos(rtk-rb,posr);/*Ri,Rj 基准站和流动站方差计算以后协方差阵对一些中间变量进行初始化将双差伪距残差和双差载波相位残差初始化为0*/Rimat(ns*nf*22,1); Rjmat(ns*nf*22,1); immat(ns,1);tropumat(ns,1); troprmat(ns,1); dtdxumat(ns,3); dtdxrmat(ns,3);/*double resp[NFREQ]; residuals of pseudorange (m) double resc[NFREQ]; residuals of carrier-phase (m) 将双差伪距残差和双差载波相位残差初始化为0每个卫星每个频点的为什么没有jNFREQ*2,因为伪距和载波分为resp和resc了*/for (i0;iMAXSAT;i) for (j0;jNFREQ;j) {rtk-ssat[i].resp[j]rtk-ssat[i].resc[j]0.0;}/* compute factors of ionospheric and tropospheric delay 如果卡尔曼滤波中包含对流层状态量调用prectrop函数计算基站和移动站的对流层延迟湿分量存在tropu[i]和tropur[i]中如果卡尔曼滤波器包含电离层状态量调用ionmapf函数分别计算基站和移动站处的投影函数*/for (i0;ins;i) {if (opt-ionooptIONOOPT_EST) {im[i](ionmapf(posu,azeliu[i]*2)ionmapf(posr,azelir[i]*2))/2.0;}if (opt-tropoptTROPOPT_EST) {tropu[i]prectrop(rtk-sol.time,posu,0,azeliu[i]*2,opt,x,dtdxui*3);tropr[i]prectrop(rtk-sol.time,posr,1,azelir[i]*2,opt,x,dtdxri*3);}}/* 对各个系统以及载波相位、伪距分别进行循环处理双重循环遍历每个系统的每个频点(*2 包括伪距和载波)哪种定位模式是伪距那我们从nf---2*nf因为伪距残差是在y[]后半是载波(动态定位模式及以后)那我们从0--2*nf因为同时使用伪距和载波*/for (m0;m6;m) /* m0:GPS/SBS,1:GLO,2:GAL,3:BDS,4:QZS,5:IRN */for (fopt-modePMODE_DGPS?0:nf;fnf*2;f) {/* search reference satellite with highest elevation选择最大的高度角卫星作为参考卫星当然我们也可以自己设置规则ssat_t ssat[MAXSAT]; satellite status 所有的共识卫星(全系统)每次都会被遍历一遍我们选出所需的系统然后再*/for (i-1,j0;jns;j) {sysirtk-ssat[sat[j]-1].sys;if (!test_sys(sysi,m)) continue;//当前共识卫星的系统和当前遍历的系统是否一样if (!validobs(iu[j],ir[j],f,nf,y)) continue; //评价不懂if (i0||azel[1iu[j]*2]azel[1iu[i]*2]) ij;//选择排序的那种,i最终是参考卫星的index}if (i0) continue;/* make DD (double difference) 做双差遍历共识卫星一个历元一个历元处理输入流动站和基准站的数据流动站卫星残差 基准站卫星残差 每颗卫星每个频点的载波和伪距y[............... ................]站间单差两个接收机(a,b)观测同一颗卫星(p)作差(M1)。卫星p肯定是它们的共识卫星星间双差参考卫星站间单差(M2),M1-M2*/for (j0;jns;j) {if (ij) continue;sysirtk-ssat[sat[i]-1].sys;//参考卫星系统sysjrtk-ssat[sat[j]-1].sys;//流动freqifreq[f%nfiu[i]*nf];freqjfreq[f%nfiu[j]*nf];if (!test_sys(sysj,m)) continue;//判断和当前遍历系统一样否if (!validobs(iu[j],ir[j],f,nf,y)) continue;//检验数据有效yy[]//设计矩阵初始化nx---待解量的个数if (H) {HiHnv*rtk-nx;for (k0;krtk-nx;k) Hi[k]0.0;}/* DD residual 作双差残差iu[i] 共识--参考卫星在基准站残差中的索引(y[fiu[i]*nf*2]-y[fir[i]*nf*2]) 共识参考卫星站间单差*/v[nv](y[fiu[i]*nf*2]-y[fir[i]*nf*2])-(y[fiu[j]*nf*2]-y[fir[j]*nf*2]);/* partial derivatives by rover position e[]是dx dy dz设计矩阵大小流动站观测到的卫星基准站观测到的卫星。 因为在relpos()函数中站间单差[先是算基准站后是流动站]设计矩阵都填到e里面了。所以大小如此设计矩阵非参考卫星-参考卫星因为我们求得是流动站位置所以只用iu[]即可至于不同频率,这里只有单频*/if (H) {for (k0;k3;k) {Hi[k]-e[kiu[i]*3]e[kiu[j]*3];}}/* DD ionospheric delay term */if (opt-ionooptIONOOPT_EST) {didxi(fnf?-1.0:1.0)*im[i]*SQR(FREQ1/freqi);didxj(fnf?-1.0:1.0)*im[j]*SQR(FREQ1/freqj);v[nv]-didxi*x[II(sat[i],opt)]-didxj*x[II(sat[j],opt)];if (H) {Hi[II(sat[i],opt)] didxi;Hi[II(sat[j],opt)]-didxj;}}/* DD tropospheric delay term */if (opt-tropoptTROPOPT_EST||opt-tropoptTROPOPT_ESTG) {v[nv]-(tropu[i]-tropu[j])-(tropr[i]-tropr[j]);for (k0;k(opt-tropoptTROPOPT_ESTG?1:3);k) {if (!H) continue;Hi[IT(0,opt)k] (dtdxu[ki*3]-dtdxu[kj*3]);Hi[IT(1,opt)k]-(dtdxr[ki*3]-dtdxr[kj*3]);}}/* DD phase-bias term 双差模糊度 */if (fnf) {if (opt-ionoopt!IONOOPT_IFLC) {v[nv]-CLIGHT/freqi*x[IB(sat[i],f,opt)]-CLIGHT/freqj*x[IB(sat[j],f,opt)];if (H) {Hi[IB(sat[i],f,opt)] CLIGHT/freqi;Hi[IB(sat[j],f,opt)]-CLIGHT/freqj;}}else {//根据站星双差方差还要减去(模糊度的derT)见理论。而在y[]的计算中还没有考虑模糊度在这里给补上//x[] 模糊度v[nv]-x[IB(sat[i],f,opt)]-x[IB(sat[j],f,opt)];if (H) {Hi[IB(sat[i],f,opt)] 1.0; //IB() 哪个卫星哪个频点的模糊度位置索引不懂Hi[IB(sat[j],f,opt)]-1.0;}}}//将计算好的残差存到ssat结构体if (fnf) rtk-ssat[sat[j]-1].resc[f ]v[nv];else rtk-ssat[sat[j]-1].resp[f-nf]v[nv];/* test innovation 检查残差是否符合标准 maxinno */if (opt-maxinno0.0fabs(v[nv])opt-maxinno) {if (fnf) {rtk-ssat[sat[i]-1].rejc[f];rtk-ssat[sat[j]-1].rejc[f];}errmsg(rtk,outlier rejected (sat%3d-%3d %s%d v%.3f)\n,sat[i],sat[j],fnf?L:P,f%nf1,v[nv]);continue;}/* SD (single-differenced) measurement error variances 每颗卫星每个频率的观测值噪声(方差) */Ri[nv]varerr(sat[i],sysi,azel[1iu[i]*2],bl,dt,f,opt);Rj[nv]varerr(sat[j],sysj,azel[1iu[j]*2],bl,dt,f,opt);/* set valid data flags 设置标志这个卫星这个频率有效 */if (opt-modePMODE_DGPS) {if (fnf) rtk-ssat[sat[i]-1].vsat[f]rtk-ssat[sat[j]-1].vsat[f]1;}else {rtk-ssat[sat[i]-1].vsat[f-nf]rtk-ssat[sat[j]-1].vsat[f-nf]1;}trace(4,sat%3d-%3d %s%d v%13.3f R%8.6f %8.6f\n,sat[i],sat[j],fnf?L:P,f%nf1,v[nv],Ri[nv],Rj[nv]);//标记一下双差的具体信息。哪两个卫星作差载波还是伪距第几个频率vflg[nv](sat[i]16)|(sat[j]8)|((fnf?0:1)4)|(f%nf);nb[b];//每个小频率有几颗卫星(残差有多少个)}b;//每个小频率}/* end of system loop *//* baseline length constraint for moving baseline 基线约束以后看 */if (opt-modePMODE_MOVEBconstbl(rtk,x,P,v,H,Ri,Rj,nv)) {vflg[nv]34;nb[b];}if (H) {trace(5,H\n); tracemat(5,H,rtk-nx,nv,7,4);}/* DD measurement error covariance 协方差矩阵 */ddcov(nb,b,Ri,Rj,nv,R);free(Ri); free(Rj); free(im);free(tropu); free(tropr); free(dtdxu); free(dtdxr);return nv;
}
