RTKLIB与GAMP的周跳探测实现

阅读RTKLIB源码后发现，RTKLIB中采取观测值文件中LLI标志（卫星失锁）、GF组合和MW组合三种方式联合探测周跳，GF组合与MW组合在RTKLIB中采取经验阈值来判断是否发生周跳。

常用周跳探测方法

  载波相位周跳探测是高精度测量的关键技术，常用的周跳探测方法各具优势，但同时也存在一些探测失效的情况，如 MW 组合周跳探测精度高，但当两个频率同时发生同样大小周跳时该方法探测失效，无几何组合去除系统性影响造成的误差十分有效，但由于噪声较大探测精度不高；电离层残差组合利用双频数据，难于确定单个载波上的周跳，而且存在对一些周跳组合不敏感的问题，但如同RTKLIB联合使用上述方法时，则是一种周跳探测的理想方案。

双频组合周跳探测

基于MW组合的周跳探测方法

$$L_{mw}=\lambda_{mw}(\phi_2-\phi_2)=\rho+\frac{f_1f_2}{f^2_1-f^2_2}I+\lambda_{mw}N_{mw}$$ $$P_{mw}=\frac{f_1P_1+f_2P_2}{f_1+f_2}=\rho+\frac{f_1f_2}{f^2_1-f^2_2}I$$

MW组合为：

$$N_{mw}=\lambda_{mw}(\phi_1-\phi_2)-\frac{f_1P_1+f_2P_2}{f_1+f_2}$$

  可看出MW组合消除了电离层延迟、钟差（卫地）、站星几何距离，只包含宽巷模糊度和噪声，非常适合探测周跳。在RTKLIB中，MW组合的周跳探测是根据经验阈值(pppos.c)来确定的：

#define THRES_MW_JUMP 10.0
static void detslp_mw(rtk_t *rtk, const obsd_t *obs, int n, const nav_t *nav)

{
    double w0,w1;
    int i,j;

    for (i=0;i<n&&i<MAXOBS;i++) {
        if ((w1=mwmeas(obs+i,nav))==0.0) continue;
        w0=rtk->ssat[obs[i].sat-1].mw;
        rtk->ssat[obs[i].sat-1].mw=w1;            
        if (w0!=0.0&&fabs(w1-w0)>THRES_MW_JUMP) {
            for (j=0;j<rtk->opt.nf;j++) rtk->ssat[obs[i].sat-1].slip[j]|=1;
        }
    }
}

  代码中mwmeas为构造MW组合的函数，当相邻历元之间的MW构造量之差大于THRES_MW_JUMP（10）时，RTKLIB将其记为周跳。
  而在很多论文中，Blewitt方法中，是根据递推方法求MW构造量之差的平均值和均方根中误差来探测的：

$${N}(i )=\overline{N}(i-1)+\frac{1}{i}[N(i)-\overline{N}(i-1)]$$ $$\sigma^2(i)=\sigma^2(i-1)+\frac{1}{i}[(N(i)-\overline{N}(i-1))^2-\sigma^2(i-1)]$$

其中，$\overline{N}(i-1)$为前i个历元宽巷模糊度的平均值；$\sigma(i)$便是前i个历元的标准差，若满足如下两个条件，则认为当前历元存在周跳：

$$|N(i)-\overline{N}(i-1)|\geq4\sigma(i-1)$$ $$|N(i)-{N}(i+1)|\leq1$$

  由于RTKLIB中原函数并非在数据读取时进行预处理，它考虑了实时定位，因此MW组合探测周跳函数放在pppos函数初始化协方差中，每次都只传入单历元的数据计算mw值，联合上个历元计算得到的mw值，判断本历元是否发生了周跳。而要根据公式看出，递推方法不仅需要上个历元及本次历元，还需要下一个历元共三个历元，因此在只传入一组观测历元的情况下是无法完成周跳探测的，若做事后的精密定位，可考虑将周跳探测单独拿出来做，或者添加一个全局变量。
  在GAMP中（基于RTKLIB二次开发的多系统非差非组合精密单点定位开源程序）改进了RTKLIB的周跳探测代码，具体的做法是定义了一个全局结构体：

PPPGlobal_t PPP_Glo;

  此全局结构体中存储了采样时间和卫星高度角，考虑这两种因素来给周跳检测赋予不同的经验阈值。具体公式为：

  公式中$E$为高度角，$R$为采样时间；$R_{GF}$单位为米或m，$R_{MW}$单位为周或m，分别为GF组合与MW组合的阈值。GAMP中单点定位的周跳探测函数为detecs在myPpp.c中，
  先贴出gamp中实现MW组合的代码：

static void detslp_mw(rtk_t *rtk, const obsd_t *obs, int n, const nav_t *nav)
{
	int i,j,nd=0,sat;
	int bSlip[MAXSAT],bLowElev=0;
	double wl0,wl1,elev,el,thres,thres0,delta[50],dmw[MAXSAT],dtmp,fact;
	double std_ex,ave_ex;
	const double rad_20=20.0*D2R;

	for (i=0;i<MAXSAT;i++) {
		dmw[i]=0.0;
		bSlip[i]=0;
	}

	//the thresh values is suitable for 30s interval
	fact=1.0;
	if (PPP_Glo.sample>=29.5) {
		if (PPP_Glo.delEp<=2) fact=1.0;
		else if (PPP_Glo.delEp<=4) fact=1.25;
		else if (PPP_Glo.delEp<=6) fact=1.5;
		else	 fact=2.0;
	}

	for (i=0;i<n&&i<MAXOBS;i++) {
		sat=obs[i].sat;
		if (timediff(PPP_Glo.tNow,PPP_Glo.ssat_Ex[sat-1].tLast)>PPP_Glo.sample) {
			PPP_Glo.ssat_Ex[sat-1].mw[1]=0.0;
			PPP_Glo.ssat_Ex[sat-1].mwIndex=0;
		}

		wl1=wlAmbMeas(obs+i,nav);
		wl0=PPP_Glo.ssat_Ex[sat-1].mw[1];

		/*sprintf(PPP_Glo.chMsg,"sat=%s, mw0=%13.3f, mw1=%13.3f\n",PPP_Glo.sFlag[sat-1].id,wl1,
			PPP_Glo.ssat_Ex[sat-1].mw[1]);
		outDebug(1,1,1);*/

		if (wl1==0.0||wl0==0.0) continue;

		elev=rtk->ssat[sat-1].azel[1];
		el=elev;

		bLowElev=0;

		if (elev<rtk->opt.elmin) {
			el=rtk->opt.elmin;
			bLowElev=1;
		}

		dtmp=el*R2D;

		if (el>=rad_20) thres=PPP_Glo.prcOpt_Ex.csThresMW;
		else	 thres=-PPP_Glo.prcOpt_Ex.csThresMW*0.1*dtmp+3*PPP_Glo.prcOpt_Ex.csThresMW;

		dmw[sat-1]=wl1-wl0;

		thres0=MIN(thres*fact,6.0);
		if (fabs(wl1-wl0)>MIN(thres*fact,6.0)) {
			bSlip[sat-1]=1;
			delta[nd++]=wl1-wl0;

			if (bLowElev) continue;

			sprintf(PPP_Glo.chMsg, "%s MW CS mw_new=%13.3f, mw_old=%13.3f, diff_abs=%13.3f, thres=%13.3f, elev=%4.1f\n", 
				PPP_Glo.sFlag[sat-1].id,wl1,wl0,wl1-wl0,thres0,elev);
			outDebug(OUTWIN,OUTFIL,OUTTIM);
		}
	}

	if (2*nd+1<=n&&nd<n-3&&nd<=3) {

	}
	else if (nd>2) {
		i=findGross(0,0,delta,nd,0,&std_ex,&ave_ex,NULL,4.0,0.3,0.2);

		if (i<=1&&std_ex<=1.0&&fabs(ave_ex)<=10.0 ) {
			sprintf(PPP_Glo.chMsg,"*** WARNING: MW CS abnormally at this epoch %4.1f  %4.2f %d %d\n",ave_ex,std_ex,nd,n);
			outDebug(OUTWIN,OUTFIL,OUTTIM);

			for (i=0;i<n&&i<MAXOBS;i++) {
				sat=obs[i].sat;
				dmw[sat-1]=0.0;
				bSlip[sat-1]=0;
			}
		}
	}

	for (i=0;i<n&&i<MAXOBS;i++) {
		sat=obs[i].sat;

		if (!bSlip[sat-1]) {
			if (fabs(dmw[sat-1])>=1.0) 
				PPP_Glo.ssat_Ex[sat-1].obsStat.iCycle=1;
		}
		else {
			for (j=0;j<rtk->opt.nf;j++) 
				rtk->ssat[sat-1].slip[j]|=1;

			PPP_Glo.ssat_Ex[sat-1].obsStat.iCycle=2;
		}
	}
}

简单的解释一下代码：

1. dmw指的是两历元间mw组合差值，bSlip为标称其是否为周跳，将两者初始化

2. fact为间隔了多个历元的倍值，delEp为间隔历元（由于可能存在历元缺失或者启用某些历元，因此用delEp标记上一次历元和当前历元差了几个）

3. 按这一历元每颗卫星遍历进行周跳判断
    3.1 wlAmbMeas计算当前历元的mw组合，当前mw记为wl1，上一历元为wl0
    3.2 在detect这个函数中，GAMP没有存当下的mw值，而是在gamPos这个函数最后的keepEpInfo函数中存入PPP_GLO.ssat_Ex.mw中（这什么操作？？）

4. 前文提到，GAMP中是根据采样时间和高度角配置阈值的。因此，接下来就是读取该历元卫星的高度角，并检查高度角是否小于我们设置的最小值（默认10°）。
GMAP中是以20°分开计算两种不同的阈值配置方式。

5. PPP_Glo.prcOpt_Ex.csThresMW这个值是在execess函数中开头位置calCsThres函数计算而来，参考前面贴出的公式。

6. 其实阈值最大就为6:*thres0=MIN(thres*fact,6.0);*

7. nd是探测出的周跳数量，后这个处理基于一个什么判断暂时没弄懂

基于GF组合的周跳探测方法

  GF组合为：

$$L_{gf}=\lambda_1\phi_1-\lambda_2\phi_2=\lambda_1N_1-\lambda_2N_2-I_1+I_2 =\lambda_1N_1-\lambda_2N_2-\frac{A}{f^2_1}+\frac{A}{f^2_2}$$ $$P_{gf}=P_1-P_2=I_1-I_2=\frac{A}{f^2_1}-\frac{A}{f^2_2}$$

  GF组合[相位-相位]组合消除了除电离层和整周模糊度外的其他项，[伪距-伪距]则消除了除电离层延迟外的其他项，当电离层变化平缓或采样时间较短时（把其当做常数），GF组合能有效探测周跳；但当电离层变化活跃时，可能无法探测较小的周跳甚至发生误判，为解决此不足，一般引入伪距观测值进行组合进一步消除电离层的影响，引入伪距观测值的同时也引入了其观测噪声等误差，使得GF观测值误差变大，因此对$P_{GF}(i)$进行拟合生成$Q_i$来滤去伪距噪声的影响，多项式拟合阶数必须满足条件$m=min[(N/100+1),6]$，$N$为历元个数。得到GF组合周跳探测模型：

$$\delta L_{GF}=L_{GF}-Q=\lambda_1N_1-\lambda_2N_2$$

  若满足以下两个条件则认为当前历元存在周跳：

$$|\delta L_{GF}(i)-\delta L_{GF}(i-1)|>6(\lambda_2-\lambda_1)$$ $$|\delta L_{GF}(i+1)-\delta L_{GF}(i)|<(\lambda_2-\lambda_1)$$

  上述方法对$P_{GF}$进行了多项式拟合，而多项式拟合设计多项式类型、拟合阶数等人为因素，降低了数据处理的准确性、客观性，且该算法的探测性能依赖于电离层的变化，对采样率要求较高。基于此，对GF组合做出一些改进，较为常用的为相邻历元求差，其去除了噪声较大的伪距观测值，仅将相邻的两个历元的相位GF组合观测值求差来探测周跳：

$$\delta L_{GF}(i)=L_{GF}(i)-L_{GF}(i-1)=\delta I+\lambda_1\delta N_1-\lambda_2\delta N_2$$

  与MW组合相似的，RTKLIB考虑实时的因素，GF组合也是用固定阈值判断是否发生了周跳，下面贴出RTKLIB中GF组合探测周跳的代码：

static void detslp_gf(rtk_t *rtk, const obsd_t *obs, int n, const nav_t *nav)
{
    double g0,g1;
    int i,j;
    for (i=0;i<n&&i<MAXOBS;i++) {
        if ((g1=gfmeas(obs+i,nav))==0.0) continue;        
        g0=rtk->ssat[obs[i].sat-1].gf;
        rtk->ssat[obs[i].sat-1].gf=g1;
        if (g0!=0.0&&fabs(g1-g0)>rtk->opt.thresslip) {      
            for (j=0;j<rtk->opt.nf;j++) rtk->ssat[obs[i].sat-1].slip[j]|=1;
        }
    }
}

  其中gfmeas为构造GF组合的函数，若g1-g0的绝对值大于用户设定的阈值thresslip，则认为该历元发生了周跳（默认阈值为0.05）。
  再说GAMP中，GF组合探测周跳是如何实现的。探测方式是与RTKLIB一样的，利用GF组合当前观测值和上一历元组合的观测值做差来检验是否发生周跳，只是在GAMP中根据高度角和采样间隔进行了阈值的配置，下面看一下代码：

/* detect cycle slip by geometry free phase jump -----------------------------*/
static void detslp_gf(rtk_t *rtk, const obsd_t *obs, int n, const nav_t *nav)
{
	double g0,g1,elev,thres=0.0,deltaEpoch=1.0;
	double elev0,thres0;
	int i,j,sat;

	deltaEpoch=fabs(rtk->tt/PPP_Glo.sample);

	for (i=0;i<n&&i<MAXOBS;i++) 
	{
		sat=obs[i].sat;
		elev=rtk->ssat[sat-1].azel[1]*R2D;
		elev0=elev;

		g1=gfmeas(obs+i,nav);

		if (g1==0.0) continue;

		if (elev<rtk->opt.elmin*R2D) elev=rtk->opt.elmin*R2D;


		if (elev>=15.0) thres=PPP_Glo.prcOpt_Ex.csThresGF;
		else thres=-PPP_Glo.prcOpt_Ex.csThresGF/15.0*elev+2*PPP_Glo.prcOpt_Ex.csThresGF;

		g0=rtk->ssat[sat-1].gf;

		thres0=MIN(thres*deltaEpoch,1.5);
		if (g0!=0.0&&fabs(g1-g0)>MIN(thres*deltaEpoch,1.5)) 
		{
			for (j=0;j<rtk->opt.nf;j++) 
				rtk->ssat[sat-1].slip[j]|=1;

			sprintf(PPP_Glo.chMsg,"%s GF CS gf_new=%13.3f, gf_old=%13.3f, diff_abs=%13.3f, thres=%13.3f, elev=%4.1f\n", 
				PPP_Glo.sFlag[sat-1].id,g1,g0,g1-g0,thres0,elev0);
			outDebug(OUTWIN,OUTFIL,OUTTIM);
		}
	}
}

  简单浏览代码，我们可以发现其和MW组合的实现方式是类似的：

1. 计算相隔历元 deltaEpoch，为后续计算阈值使用

2. 取得当前卫星的高度角，计算当前历元GF组合观测值

3. 检查GF观测值和卫星高度角的值

4. 以15°的高度角为分界，两种不同权值计算阈值，取其和1.5中的最小值：*thres0=MIN(thres*deltaEpoch,1.5);*
5. 历元间GF做差与阈值比较

——————————————————————————————————————————————

  GAMP中的周跳探测较RTKLIB有较大改善，具体实验算例在下次[有空….时]做。另，许多论文都针对TurboEdit（MW+GF）法有许多新的改进方法，也会在下次中配合算例实现。

最后更新： 2019年07月01日 00:00

原始链接： https://wu941202.github.io/2019/07/14/2019-07-14/

WuFF