网站建设的系统设计,程序开发是什么专业,建设部作业证件查询官方网站,苏州市建设局网站首页文章目录 一、固定模糊度的前置工作1. 做好固定模糊度的准备2. 建立双差模糊度3. 问题与总结 版权声明#xff1a;本文为原创文章#xff0c;版权归 Winston Qu 所有#xff0c;转载请注明出处。 在上一篇文章中#xff0c;介绍了RTKLIB中manage_amb_LAMBDA()函数#xff… 文章目录 一、固定模糊度的前置工作1. 做好固定模糊度的准备2. 建立双差模糊度3. 问题与总结 版权声明本文为原创文章版权归 Winston Qu 所有转载请注明出处。 在上一篇文章中介绍了RTKLIB中manage_amb_LAMBDA()函数并详细介绍了其操作方式和工作方法。可见[RTKLIB]模糊度固定相关问题(一) 本篇文章中我们针对其模糊度固定函数resamb_LAMBDA()、ddidx()做详细的分析进一步的探索模糊度固定中有意思的部分。
一、固定模糊度的前置工作
1. 做好固定模糊度的准备 在前一篇文中我们说到为了提高浮点模糊度的固定率进行了部分模糊度*、延迟上星的操作。但究其根本还没有真正接触到模糊度固定的底层并没有创建模糊度双差数组也没有去建立模糊度求解的矩阵方程更还没尝试通过整数最小二乘方法去固定模糊度。 模糊度固定是一个比较复杂和繁琐的过程前期会有许多的准备过程。从状态量中提取我们需要的状态通过方差-协方差阵计算我们需要的Q矩阵。然后通过LAMBDA方法去求解固定后的状态向量并通过一定的方法转换到单差模糊度修正已有状态量。… 上述过程听着就挺复杂的我们可以通过函数仔细看看做了哪些工作从中理解设计和求解的思想。 附上源代码
/* resolve integer ambiguity by LAMBDA ---------------------------------------*/
static int resamb_LAMBDA(rtk_t *rtk, double *bias, double *xa,int gps,int glo,int sbs)
{prcopt_t *optrtk-opt;int i,j,nb,nb1,info,nxrtk-nx,nartk-na;double *DP,*y,*b,*db,*Qb,*Qab,*QQ,s[2];int *ix;double var0,coeff[3];double QQb[MAXSAT];trace(3,resamb_LAMBDA : nx%d\n,nx);rtk-sol.ratio0.0;rtk-nb_ar0;if (rtk-opt.modePMODE_DGPS||rtk-opt.modearARMODE_OFF||rtk-opt.thresar[0]1.0) {return 0;}/* skip AR if position variance too high to avoid false fix */for (i0;i3;i) varrtk-P[ii*rtk-nx];varvar/3.0; /* maintain compatibility with previous code */trace(3,posvar%.6f\n,var);if (varrtk-opt.thresar[1]) {errmsg(rtk,position variance too large: %.4f\n,var);return 0;}/* Create index of single to double-difference transformation matrix (D)used to translate phase biases to double difference */iximat(nx,2);if ((nbddidx(rtk,ix,gps,glo,sbs))(rtk-opt.minfixsats-1)) { /* nb is sat pairs */errmsg(rtk,not enough valid double-differences\n);free(ix);return -1; /* flag abort */}rtk-nb_arnb;/* nx# of float states, na# of fixed states, nb# of double-diff phase biases */ymat(nb,1); DPmat(nb,nx-na); bmat(nb,2); dbmat(nb,1); Qbmat(nb,nb);Qabmat(na,nb); QQmat(na,nb);/* phase-bias covariance (Qb) and real-parameters to bias covariance (Qab) *//* yD*xc, QbD*Qc*D, QabQac*D */for (i0;inb;i) {y[i]rtk-x[ix[i*2]]-rtk-x[ix[i*21]];}for (j0;jnx-na;j) for (i0;inb;i) {DP[ij*nb]rtk-P[ix[i*2](naj)*nx]-rtk-P[ix[i*21](naj)*nx];}for (j0;jnb;j) for (i0;inb;i) {Qb[ij*nb]DP[i(ix[j*2]-na)*nb]-DP[i(ix[j*21]-na)*nb];}for (j0;jnb;j) for (i0;ina;i) {Qab[ij*na]rtk-P[iix[j*2]*nx]-rtk-P[iix[j*21]*nx];}for (i0;inb;i) QQb[i]1000*Qb[ii*nb];trace(3,N(0) ); tracemat(3,y,1,nb,7,2);trace(3,Qb*1000 ); tracemat(3,QQb,1,nb,7,4);/* lambda/mlambda integer least-square estimation *//* return best integer solutions *//* b are best integer solutions, s are residuals */if (!(infolambda(nb,2,y,Qb,b,s))) {trace(3,N(1) ); tracemat(3,b ,1,nb,7,2);trace(3,N(2) ); tracemat(3,bnb,1,nb,7,2);rtk-sol.ratios[0]0?(float)(s[1]/s[0]):0.0f;if (rtk-sol.ratio999.9) rtk-sol.ratio999.9f;/* adjust AR ratio based on # of sats, unless minARmaxAR */if (opt-thresar[5]!opt-thresar[6]) {nb1nb50?nb:50; /* poly only fitted for upto 50 sat pairs *//* generate poly coeffs based on nominal AR ratio */for ((i0);i3;i) {coeff[i] ar_poly_coeffs[i][0];for ((j1);j5;j)coeff[i] coeff[i]*opt-thresar[0]ar_poly_coeffs[i][j];}/* generate adjusted AR ratio based on # of sat pairs */rtk-sol.thres coeff[0];for (i1;i3;i) {rtk-sol.thres rtk-sol.thres*1/(nb11)coeff[i];}rtk-sol.thres MIN(MAX(rtk-sol.thres,opt-thresar[5]),opt-thresar[6]);} elsertk-sol.thres(float)opt-thresar[0];/* validation by popular ratio-test of residuals*/if (s[0]0.0||s[1]/s[0]rtk-sol.thres) {/* init non phase-bias states and covariances with float solution values *//* transform float to fixed solution (xaxa-Qab*Qb\(b0-b)) */for (i0;ina;i) {rtk-xa[i]rtk-x[i];for (j0;jna;j) rtk-Pa[ij*na]rtk-P[ij*nx];}/* y differences between float and fixed dd phase-biasesbias fixed dd phase-biases */for (i0;inb;i) {bias[i]b[i];y[i]-b[i];}/* adjust non phase-bias states and covariances using fixed solution values */if (!matinv(Qb,nb)) { /* returns 0 if inverse successful *//* rtk-xa rtk-x-Qab*Qb^-1*(b0-b) */matmul(NN,nb,1,nb, 1.0,Qb ,y,0.0,db); /* db Qb^-1*(b0-b) */matmul(NN,na,1,nb,-1.0,Qab,db,1.0,rtk-xa); /* rtk-xa rtk-x-Qab*db *//* rtk-Partk-P-Qab*Qb^-1*Qab) *//* covariance of fixed solution (QaQa-Qab*Qb^-1*Qab) */matmul(NN,na,nb,nb, 1.0,Qab,Qb ,0.0,QQ); /* QQ Qab*Qb^-1 */matmul(NT,na,na,nb,-1.0,QQ ,Qab,1.0,rtk-Pa); /* rtk-Pa rtk-P-QQ*Qab */trace(3,resamb : validation ok (nb%d ratio%.2f thresh%.2f s%.2f/%.2f)\n,nb,s[0]0.0?0.0:s[1]/s[0],rtk-sol.thres,s[0],s[1]);/* translate double diff fixed phase-bias values to single diff fix phase-bias values, result in xa */restamb(rtk,bias,nb,xa);}else nb0;}else { /* validation failed */errmsg(rtk,ambiguity validation failed (nb%d ratio%.2f thresh%.2f s%.2f/%.2f)\n,nb,s[1]/s[0],rtk-sol.thres,s[0],s[1]);nb0;}}else {errmsg(rtk,lambda error (info%d)\n,info);nb0;}free(ix);free(y); free(DP); free(b); free(db); free(Qb); free(Qab); free(QQ);return nb; /* number of ambiguities */
}我们来拆解步骤并给出详细的解释。
准入条件判断根据定为模式、位置方差判断当前是否适合进入模糊度固定。ps:可以想想如果位置方差很大为什么不适合进入模糊度固定呢申请了一个im(nx,2)的空间用于保存构建双差浮点模糊度的索引值。之后通过ddidx()函数构建了双差浮点模糊度矩阵在此处做了一个模糊度个数的判断判断是否小于设定的最小固定卫星。ps:若出现三颗三频的卫星此处的nb(3-1)*36确实满足最小卫星4的条件但是此处的固定是否有意义呢我觉得这是一个版本迭代的BUG根据公式计算双差浮点模糊度y、双差协方差矩阵Qb、双差实参协方差矩阵Qab。ps:如果从理论角度来说在此处我们需要求四个矩阵分别为Qa、Qab、Qba、Qb但实质上我们只用到Qb、Qab就能完成整个计算流程。输出固定前的浮点双差模糊度N(0)和模糊度的方差Qb*1000在调试模糊度固定的时候这两个值是非常重要的两个值。使用lambda()函数进行浮点双差模糊度的固定如果解算失败就释放申请的内存若解算成功则进行ratio test。ps:这里的解算成功不是指成功求解整数模糊度而是指成功进行lambda()解算输出计算过后的整数双差模糊度的最优解和次优解如果设置了动态AR阈值还需要根据卫星数计算模糊度阈值。进行ratio test判断模糊度最优解与次优解的残差是否合规模糊度是否可以被采纳。将计算的整数双差模糊度通过计算更新到固定解位置和相应的方差-协方差解算矩阵上并通过restamb()函数将整数双差模糊度转换到对应的浮点单差模糊度。ps:保存在xa变量中了,此处的restamb()与ddidx()互为反映射关系返回模糊度个数nb。
2. 建立双差模糊度 如果要进行模糊度解算肯定需要在系统中简历双差模糊度方程如何做双差?选择怎样的卫星作为第一颗卫星?其中的奥秘都集中在ddidx()函数中。 附上源代码
/* index for single to double-difference transformation matrix (D) --------------------*/
static int ddidx(rtk_t *rtk, int *ix, int gps, int glo, int sbs)
{int i,j,k,m,f,n,nb0,nartk-na,nfNF(rtk-opt),nofix;double fix[MAXSAT],ref[MAXSAT];trace(3,ddmat: gps%d/%d glo%d/%d sbs%d\n,gps,rtk-opt.gpsmodear,glo,rtk-opt.glomodear,sbs);/* clear fix flag for all sats (1float, 2fix) */for (i0;iMAXSAT;i) for (j0;jNFREQ;j) {rtk-ssat[i].fix[j]0;}for (m0;m6;m) { /* m0:GPS/SBS,1:GLO,2:GAL,3:BDS,4:QZS,5:IRN *//* skip if ambiguity resolution turned off for this sys */nofix(m0gps0)||(m1glo0)||(m3rtk-opt.bdsmodear0); /* step through freqs */ for (f0,kna;fnf;f,kMAXSAT) {/* look for first valid sat (istate index, i-ksat index) */for (ik;ikMAXSAT;i) {/* skip if sat not active */if (rtk-x[i]0.0||!test_sys(rtk-ssat[i-k].sys,m)||!rtk-ssat[i-k].vsat[f]) {continue;}/* set sat to use for fixing ambiguity if meets criteria */if (rtk-ssat[i-k].lock[f]0!(rtk-ssat[i-k].slip[f]2)rtk-ssat[i-k].azel[1]rtk-opt.elmaskar!nofix) {rtk-ssat[i-k].fix[f]2; /* fix */break;/* break out of loop if find good sat */}/* else dont use this sat for fixing ambiguity */else rtk-ssat[i-k].fix[f]1;}if (rtk-ssat[i-k].fix[f]!2) continue; /* no good sat found *//* step through all sats (jstate index, j-ksat index, i-kfirst good sat) */for (n0,jk;jkMAXSAT;j) {if (ij||rtk-x[j]0.0||!test_sys(rtk-ssat[j-k].sys,m)||!rtk-ssat[j-k].vsat[f]) {continue;}if (sbs0 satsys(j-k1,NULL)SYS_SBS) continue; if (rtk-ssat[j-k].lock[f]0!(rtk-ssat[j-k].slip[f]2)rtk-ssat[j-k].vsat[f]rtk-ssat[j-k].azel[1]rtk-opt.elmaskar!nofix) {/* set D coeffs to subtract sat j from sat i */ix[nb*2 ]i; /* state index of ref bias */ix[nb*21]j; /* state index of target bias *//* inc # of sats used for fix */ref[nb]i-k1;fix[nb]j-k1;rtk-ssat[j-k].fix[f]2; /* fix */n; /* count # of sat pairs for this freq/constellation */}/* else dont use this sat for fixing ambiguity */else rtk-ssat[j-k].fix[f]1;}/* dont use ref sat if no sat pairs */if (n0) rtk-ssat[i-k].fix[f]1;}}if (nb0) {trace(3,refSats);tracemat(3,ref,1,nb,7,0);trace(3,fixSats);tracemat(3,fix,1,nb,7,0);}return nb;
}我们来拆解步骤并给出详细的解释。
对所有的卫星rtk-ssat[i].fix[j]结构体进行清空进行初始化。对每个星座的每个频率进行双重循环确定每个星座和每个频率中的参考卫星第一颗符合规定的卫星对符合条件的卫星设定rtk-ssat[i-k].fix[f]2; /* fix */对有载波观测值但不符合条件的卫星设定rtk-ssat[i-k].fix[f]1;完成卫星清洗工作选出参考卫星。ps:本段存在很大的问题在后面详细叙述。上一个循环已经确定了参考卫星在这个循环中确定同星座同频率的非参考卫星与参考卫星组pairs并记录对应的索引号、卫星号。返回双差卫星(模糊度)个数。
3. 问题与总结 通过上述的两个函数我们就完成了除整数最小二乘(LAMBDA)外的所有模糊度固定的流程。表面上似乎没问题但通过仔细的解析和细细回味可以发现很多不合理的地方。
选择参考星的时候并没有选择高度角最高的卫星ps:这个问题其实困扰我一段时间了似乎如果只更新xa、xp选择高度角最高的卫星和系统内的第一颗卫星并没有很大的差别双差模糊度nb和参与计算的卫星ns之争。在上述的代码中使用nb判断这次解算是否有效这似乎是一个随着版本迭代的bug。在单频状态下nb确实可以代表参与解算的卫星数但随着频率的增加多频模糊度固定时nb的数量可能会很大但是真实参与解算和固定的卫星ns可能不满足固定的最小卫星数此时固定的解算结果有意义吗ps:通过固定的模糊度修正的nx实际上是秩亏的个人认为是不正确的。restamb()与ddidx()互为映射关系但其中的判断条件太简单了且需要丰富判断条件。
版权声明本文为原创文章版权归 Winston Qu 所有转载请注明出处。
