9.7 detslp*(): 周跳探测
9.7.1 周跳探测与修复概念
周跳的探测与修复就是运用一定的方法探测出载波相位是否发生了整周跳变,并求出丢失的整周数,然后将中断后的整周计数恢复为正确的计数,使这部分观测值能正常使用。在工程实践中,通常周跳探测结束后并不会进行修复操作,而是将该卫星从固定的卫星组中排除出去。
1. 周跳产生原因(信号测量层面)
接收机开机后,先测出载波相位中的“小数部分”,并把整数部分设为 0,随后开始计数。只要信号一直锁得住,每当相位走过整整一圈(2π→0),计数器就加 1。因此,在任何时刻,测得的累积相位 Δφ 就是“整数计数 N”加上“小数部分”。这里的 N 是个未知常数;只要不失锁,N 保持不变。一旦信号中断(比如被高楼、树叶或人体遮挡),计数器就停走。信号重新捕获后,小数部分仍然正确,但整数部分从零重新开始,于是 N 突然变了——这就是“周跳” [23]。
造成信号中断、进而导致周跳的常见原因有:
- 视线被遮挡(动态测量中最常见);
- 信噪比太低,如电离层扰动、多路径、接收机高速运动或卫星高度角过低;
- 接收机软件缺陷或卫星钟异常。
2. 周跳对定位的影响
超过 10 周的周跳在预处理阶段很容易被发现并剔除;真正棘手的是 1–5 周的小周跳,以及半周、1/4 周这类“微周跳”。它们常被当成随机噪声,结果把系统性偏差伪装成偶然误差,直接拉低定位精度。
哪怕只有一颗卫星出现 1 周跳,也可能让点位偏移几厘米。要知道,单点定位至少需要 4 颗卫星,因此最终误差取决于卫星数量、几何构型、周跳在各方向上的分量大小及其出现次数。
3. 周跳探测方法
表9.7-1 单频周跳探测方法汇总
方法 | 简述 | 优点 | 局限 |
---|---|---|---|
相位伪距组合法 | 相位与伪距差分,扣除电离层变化后,剩余模糊度作为检验量。 | 简单适用单频数据。 | 伪距噪声大,易受多路径影响。 |
多项式拟合法 | 用前m个历元拟合模型,预测下一历元值与实际值比较检测周跳。 | 适用于连续数据。 | 阶数和窗口大小难确定,对小周跳敏感性低。 |
多普勒积分法 | 相位差减多普勒积分值作为检验量,差值在噪声范围内则无周跳。 | 稳定,独立于相位,适合动态数据。 | 采样间隔影响结果。 |
小波变换法 | 分解信号为不同频率,检测阶跃点识别周跳。 | 对跳变敏感。 | 计算慢,函数选择难,需权衡去噪与速度。 |
表9.7-2 多频周跳探测方法汇总
方法 | 简述 | 优点 | 局限 |
---|---|---|---|
TurboEdit 方法 | MW组合+GF组合,伪距GF拟合消噪后与相位GF相减,检测周跳。 | 综合MW和GF优势,覆盖性强。 | 伪距噪声大,需拟合处理。 |
9.7.2 detslp_ll():根据 LLI 检测周跳
detslp_ll()
函数根据 LLI 失锁标识符进行周跳探测(来自接收机芯片)。
1. 基本理论
LLI 范围为 0-7,对应二进制数 000~111,其中包含有 3 个 bit 位(bit0、bit1 仅用于相位):
位 | 名称 | 掩码 | 含义(仅对载波相位有效) |
---|---|---|---|
bit0 | slip flag | 0×01 | =1 表示“前一历元到当前历元间发生过信号失锁”,提示可能发生周跳。 |
bit1 | half-cycle flag | 0×02 | =1 表示接收机对该历元做了半周模糊度解算,或软件无法处理半周数据而跳记;若此位前后历元不一致,也视为周跳。 |
bit2 | AS flag | 0×04 | =1 表示该观测值是在反欺骗(Anti-Spoofing)状态下获得,噪声可能增大,一些软件也会把它当作可疑标识。 |
一些性能比较好的接收机,其自身的 LLI 标识就能把载波数据的周跳、半周跳已足够有效,就不需要其他探测方法了。LLI 的是从信号层面进行周跳探测的,效果理应比软件层面的方法更为可靠。
2. 参数列表
/* args */
rtk_t *rtk IO rtk 控制结构体
const obsd_t *obs I 观测数据
int i I 当前观测的索引
int rcv I 接收机类型(1 为 rover,2 为 base)
/* return */
none
3. 执行流程
a. 数据有效性检查:
- 若载波数据为空或数据重复(
obs[i].time
与rtk->ssat[sat-1].pt[rcv-1][f]
接近):- 跳过当前频率,避免重复操作。
- 获取前次 LLI:
- 流动站 (
rcv = 1
):getbitu()
提取 slip[f] 的 bit[0:2]。 - 基站 (
rcv = 2
):getbitu()
提取 slip[f] 的 bit[2:4]。
- 流动站 (
b. 周跳检测:
- 正向处理 (
rtk->tt>=0.0
):- 使用当前历元的LLI:
slip = obs[i].LLI[f]
。
- 使用当前历元的LLI:
- 逆向处理 (
rtk->tt<0.0
):- 使用前一历元的LLI:
slip = LLI
。
- 使用前一历元的LLI:
c. 半周跳检测:
- 若 bit1 状态变化:
(LLI & 2) != (obs[i].LLI[f] & 2)
:- 标记周跳:
slip |= 1
。
- 标记周跳:
d. 状态更新:
- 保存当前 LLI:
- 流动站:使用
setbitu()
更新 slip[f] 的 bit[0:2])。 - 基站:使用
setbitu()
更新 slip[f] 的 bit[2:4])。
- 流动站:使用
- 更新周跳和半周标识:
- 保存周跳标识:
rtk->ssat[sat-1].slip[f] |= (uint8_t)slip
。 - 保存半周跳标识:
rtk->ssat[sat-1].half[f] = (obs[i].LLI[f] \& 2) ? 0 : 1
。
- 保存周跳标识:
4. 源码注释
点击查看代码
static void detslp_ll(rtk_t *rtk, const obsd_t *obs, int i, int rcv)
{
uint32_t slip,LLI;
int f,sat=obs[i].sat;
trace(3,"detslp_ll: i=%d rcv=%d\n",i,rcv);
// nf 是频点个数 (1:L1,2:L1+L2,3:L1+L2+L5)
for (f=0;f<rtk->opt.nf;f++) {
if (obs[i].L[f]==0.0||
fabs(timediff(obs[i].time,rtk->ssat[sat-1].pt[rcv-1][f]))<DTTOL) {
continue;
}
/* restore previous LLI */
if (rcv==1) LLI=getbitu(&rtk->ssat[sat-1].slip[f],0,2); /* rover */
else LLI=getbitu(&rtk->ssat[sat-1].slip[f],2,2); /* base */
/* detect slip by cycle slip flag in LLI */
if (rtk->tt>=0.0) { /* forward */
// bit0 为 1,表明当前历元可能发生了周跳
if (obs[i].LLI[f]&1) {
errmsg(rtk,"slip detected forward (sat=%2d rcv=%d F=%d LLI=%x)\n",
sat,rcv,f+1,obs[i].LLI[f]);
}
slip=obs[i].LLI[f];
}
else { /* backward */
// 前一历元 bit0 位为 1,表明前一历元发生可能发生了周跳
if (LLI&1) {
errmsg(rtk,"slip detected backward (sat=%2d rcv=%d F=%d LLI=%x)\n",
sat,rcv,f+1,LLI);
}
slip=LLI;
}
/* detect slip by parity unknown flag transition in LLI */
// (LLI&2)&&!(obs[i].LLI[f]&2):前一历元 bit1 位为 1 并且当前历元 bit1 位不为 1
// (!(LLI&2)&&(obs[i].LLI[f]&2):前一历元 bit1 位为 0 且当前历元 bit1 位为 1
// 即前后历元的 bit1 位不相同,表明可能发生了半周跳。
if (((LLI&2)&&!(obs[i].LLI[f]&2))||(!(LLI&2)&&(obs[i].LLI[f]&2))) {
errmsg(rtk,"slip detected half-cyc (sat=%2d rcv=%d F=%d LLI=%x->%x)\n",
sat,rcv,f+1,LLI,obs[i].LLI[f]);
slip|=1;
}
/* save current LLI */
if (rcv==1) setbitu(&rtk->ssat[sat-1].slip[f],0,2,obs[i].LLI[f]);
else setbitu(&rtk->ssat[sat-1].slip[f],2,2,obs[i].LLI[f]);
/* save slip and half-cycle valid flag */
rtk->ssat[sat-1].slip[f]|=(uint8_t)slip;
rtk->ssat[sat-1].half[f]=(obs[i].LLI[f]&2)?0:1;
}
}
9.7.3 detslp_gf():利用几何无关组合检测周跳
detslp_gf()
函数是 RTKLIB 中通过几何自由(geometry-free, GF)相位跳跃检测载波相位周跳(cycle slip)的静态函数。该函数利用多频率观测数据之间的几何自由线性组合(GF LC),通过相位跳跃幅度与用户定义的阈值比较,识别周跳并更新卫星状态。
1. 基本原理
a. 载波相位站间单差:
b. 载波测不间频差作差(几何无关组合):
c. 几何无关组合历元间作差:
其中:
- :L1 和 L2 载波的波长。
- :接收机 和 之间的 L1 和 L2 载波相位观测值。
- :接收机 和 之间的几何距离(包括传播延迟)。
- :L1 和 L2 载波的整数周数模糊度。
- :两个连续的观测历元。
上式理论等于0,若不为0(大于阈值0.05m),则认为检测到了周跳
2. 参数列表
/* args */
rtk_t *rtk IO rtk 控制结构体
const obsd_t *obs I 观测数据
int i I 流动站观测索引
int j I 基站观测索引
const nav_t *nav I 导航数据
/* return */
none
3. 执行流程
a. GF 组合计算与检测:
- 遍历频率
k
从 1 到rtk->opt.nf-1
(跳过 L1,处理 L2/L5 等):- 调用
gfobs(obs, i, j, k, nav)
计算单差 GF 相位组合gf1
:- (其中 为单差相位, 为光速, 为 L1 和 L2/L5 频率)。
- 若
gf1 = 0.0
(无效数据,如缺失观测),跳过。
- 获取前次 GF 值:
gf0 = rtk->ssat[sat-1].gf[k-1]
(存储在前次历元)。 - 更新当前 GF 值:
rtk->ssat[sat-1].gf[k-1] = gf1
(为下次比较准备)。
- 调用
b. 周跳检测:
- 若
gf0 != 0.0
且|gf1 - gf0| > rtk->opt.thresslip
:- 标记周跳:
rtk->ssat[sat-1].slip[0] |= 1
(L1 频率周跳)。rtk->ssat[sat-1].slip[k] |= 1
(当前频率 $ k $ 周跳)。
- 标记周跳:
4. 源码注释
点击查看代码
static void detslp_gf(rtk_t *rtk, const obsd_t *obs, int i, int j,
const nav_t *nav)
{
int k, sat = obs[i].sat; /* 频率索引、卫星编号 (从 obs[i].sat 获取) */
double gf0, gf1; /* 前次和当前几何自由相位组合值 (单位:米) */
trace(4, "detslp_gf: i=%d j=%d\n", i, j); /* 调试:输出流动站和基站索引 */
/* 跳过检查:若阈值关闭或已检测周跳 */
if (rtk->opt.thresslip == 0) return; /* 阈值关闭,跳过检测 */
for (k = 0; k < rtk->opt.nf; k++)
if (rtk->ssat[sat - 1].slip[k] & 1) return; /* 已检测周跳,跳过 */
/* 遍历频率,检测 GF 跳跃 (从 L2/L5 开始) */
for (k = 1; k < rtk->opt.nf; k++) {
/* 计算单差几何自由相位组合 */
if ((gf1 = gfobs(obs, i, j, k, nav)) == 0.0) continue; /* 无效数据跳过 */
gf0 = rtk->ssat[sat - 1].gf[k - 1]; /* 获取前次 GF 值 */
rtk->ssat[sat - 1].gf[k - 1] = gf1; /* 保存当前 GF 值供下次使用 */
/* 检测相位跳跃 */
if (gf0 != 0.0 && fabs(gf1 - gf0) > rtk->opt.thresslip) {
rtk->ssat[sat - 1].slip[0] |= 1; /* 标记 L1 周跳 */
rtk->ssat[sat - 1].slip[k] |= 1; /* 标记当前频率周跳 */
errmsg(rtk, "slip detected GF jump (sat=%2d L1-L%d dGF=%.3f)\n",
sat, k + 1, gf0 - gf1); /* 记录错误,输出跳跃差值 */
}
}
}
9.7.4 detslp_dop():利用多普勒探测周跳
detslp_dop()
函数是 RTKLIB 中通过多普勒和相位差检测载波相位周跳(cycle slip)的静态函数。该函数通过比较历元间载波相位差与多普勒积分的差异,结合平均多普勒差值,识别可能的周跳。然而,由于 RTKLIB 2.4.3 中因 clock-jump(时间跳变)问题,该函数已不再使用。
1. 参数列表
/* args */
rtk_t *rtk IO rtk 控制结构体
const obsd_t *obs I 观测数据
int ix I 卫星索引数组
int ns I 卫星数目
const nav_t *nav I 导航数据
/* return */
none
2. 执行流程
a. 多普勒差值计算:
- 遍历有效卫星 (索引为
ix[i]
,总数ns
):- 提取卫星编号:
sat = obs[ii].sat
(ii = ix[i])。 - 遍历频率 ():
- 初始化
dopdif[i][f] = 0
,tt[i][f] = 0.00
。 - 若
obs[ii].L[f] == 0.0
(相位无效)、obs[ii].D[f] == 0.0
(多普勒无效)或rtk->ssat[sat-1].ph[rcv-1][f] == 0.0
(前次相位无效),跳过。 - 计算时间间隔:。
- 若 (时间差过小),跳过。
- 计算相位差:(单位:周期/秒)。
- 计算多普勒积分:(单位:周期,近似 )。
- 计算差异:(单位:周期/秒)。
- 若 (非异常值),累加到
mean_dop
,计数ndop++
。
- 初始化
- 提取卫星编号:
b. 平均多普勒差值:
- 若
ndop
== 0(无有效数据,通常因钟差过大),返回。 - 计算平均差值:
mean_dop = mean_dop / ndop
(单位:周期/秒,近似钟差影响)。
c. 周跳检测:
- 再次遍历卫星 和频率 :
- 若
dopdif[i][f] == 0.00
(无有效数据),跳过。 - 若
fabs(dopdif[i][f] - mean_dop) > rtk->opt.thresdop
:- 标记周跳:
rtk->ssat[sat-1].slip[f] |= 1
。
- 标记周跳:
- 若
3. 源码注释
点击查看代码
/* detect cycle slip by doppler and phase difference -------------------------*/
/*
* 通过多普勒和相位差检测载波相位周跳。
* 参数:rtk - RTK 控制结构,obs - 观测数据,ix - 卫星索引数组,ns - 有效卫星数,rcv - 接收机编号,nav - 导航数据。
* 输出:更新 `rtk->ssat[sat-1].slip` (周跳标志,bit0 = 1 表示周跳)。
* 备注:因 clock-jump 问题,已在 RTKLIB 2.4.3 中废弃。
*/
static void detslp_dop(rtk_t *rtk, const obsd_t *obs, const int *ix, int ns,
int rcv, const nav_t *nav)
{
int i, ii, f, sat, ndop = 0, nf = rtk->opt.nf; /* 索引、卫星数、有效计数、频率数 */
double dph, dpt, mean_dop = 0; /* 相位差、多普勒积分、平均差值 */
double dopdif[MAXSAT][NFREQ], tt[MAXSAT][NFREQ]; /* 差值数组、时间间隔数组 */
trace(4, "detslp_dop: rcv=%d\n", rcv); /* 调试:输出接收机编号 */
if (rtk->opt.thresdop <= 0) return; /* 阈值关闭,跳过检测 */
/* 计算所有卫星和频率的多普勒差值 */
for (i = 0; i < ns; i++) {
ii = ix[i];
sat = obs[ii].sat;
for (f = 0; f < nf; f++) {
dopdif[i][f] = 0; tt[i][f] = 0.00;
if (obs[ii].L[f] == 0.0 || obs[ii].D[f] == 0.0 || rtk->ssat[sat - 1].ph[rcv - 1][f] == 0.0) continue;
if (fabs(tt[i][f] = timediff(obs[ii].time, rtk->ssat[sat - 1].pt[rcv - 1][f])) < DTTOL) continue;
/* 计算相位差和多普勒积分 (单位:周期/秒) */
dph = (obs[ii].L[f] - rtk->ssat[sat - 1].ph[rcv - 1][f]) / tt[i][f];
dpt = -obs[ii].D[f]; /* 矩形积分近似 */
dopdif[i][f] = dph - dpt;
/* 若非异常值,计入平均值 */
if (fabs(dopdif[i][f]) < 3 * rtk->opt.thresdop) {
mean_dop += dopdif[i][f];
ndop++;
}
}
}
/* 计算平均多普勒差值 (近似钟差影响) */
if (ndop == 0) return; /* 无有效数据,跳过 */
mean_dop = mean_dop / ndop;
/* 检测周跳:移除平均值后超出阈值 */
for (i = 0; i < ns; i++) {
sat = obs[ix[i]].sat;
for (f = 0; f < nf; f++) {
if (dopdif[i][f] == 0.00) continue;
if (fabs(dopdif[i][f] - mean_dop) > rtk->opt.thresdop) {
rtk->ssat[sat - 1].slip[f] |= 1; /* 标记周跳 */
errmsg(rtk, "slip detected doppler (sat=%2d rcv=%d dL%d=%.3f off=%.3f tt=%.2f)\n",
sat, rcv, f + 1, dopdif[i][f] - mean_dop, mean_dop, tt[i][f]);
}
}
}
}
9.7.5 detslp_mw():利用 MW 组合检测周跳
detslp_mw()
函数是 RTKLIB 中通过墨尔本-乌贝纳(Melbourne-Wubbena, MW)线性组合跳跃检测载波相位周跳(cycle slip)的静态函数。该方法利用载波相位和伪距的线性组合,生成对几何距离和离子层延迟不敏感的 MW 组合,通过比较历元间的跳跃幅度检测周跳。
1. 参数列表
/* args */
rtk_t *rtk IO RTK 控制结构体
const obsd_t *obs I 观测数据
int n I 有效观测数
const nav_t *nav I 导航数据
/* return */
none
2. 执行流程
a. MW 组合计算:
- 调用
mwmeas(obs+i, nav)
计算当前 MW 值 :- 若 (无效数据),跳过。
- 获取前次 MW 值:
w0 = rtk->ssat[obs[i].sat-1].mw[0]
。 - 更新当前 MW 值:
rtk->ssat[obs[i].sat-1].mw[0] = w1
。
b. 周跳检测:
- 若 且 :
- 遍历所有频率 (),标记周跳:rtk->ssat[obs[i].sat-1].slip[j] |= 1。
3. 源码解析
/* detect slip by Melbourne-Wubbena linear combination jump ------------------*/
/*
* 通过墨尔本-乌贝纳线性组合跳跃检测载波相位周跳。
* 参数:rtk - RTK 控制结构,obs - 观测数据,n - 有效观测数,nav - 导航数据。
* 输出:更新 `rtk->ssat[obs[i].sat-1].slip` (周跳标志,bit0 = 1 表示周跳)。
*/
static void detslp_mw(rtk_t *rtk, const obsd_t *obs, int n, const nav_t *nav)
{
double w0, w1; /* 前次和当前 MW 组合值 (单位:米) */
int i, j; /* 观测索引、频率索引 */
trace(4, "detslp_mw: n=%d\n", n); /* 调试:输出有效观测数 */
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[0]; /* 获取前次 MW 值 */
rtk->ssat[obs[i].sat - 1].mw[0] = w1; /* 保存当前 MW 值 */
trace(4, "detslip_mw: sat=%2d mw0=%8.3f mw1=%8.3f\n", obs[i].sat, w0, w1); /* 调试:输出 MW 值 */
if (w0 != 0.0 && fabs(w1 - w0) > THRES_MW_JUMP) {
trace(3, "detslip_mw: slip detected sat=%2d mw=%8.3f->%8.3f\n",
obs[i].sat, w0, w1); /* 记录周跳 */
for (j = 0; j < rtk->opt.nf; j++) rtk->ssat[obs[i].sat - 1].slip[j] |= 1; /* 标记所有频率周跳 */
}
}
}
9.7.6 gfobs():求几何无关组合观测值
1. 参数列表
/* args */
const obsd_t *obs I 观测数据
int i I 流动站观测索引
int j I 基站观测索引
int k I 频率索引
const nav_t *nav I 导航数据
/* return */
double gf O 几何无关组合观测值
2. 源码注释
/* single-differenced geometry-free linear combination of phase --------------*/
/*
* 计算单差几何自由相位组合,用于检测载波相位周跳。
* 参数:obs - 观测数据,i - 流动站观测索引,j - 基站观测索引,k - 目标频率索引,nav - 导航数据。
* 输出:返回单差 GF 组合值 (单位:米),无效时返回 0.0。
*/
static double gfobs(const obsd_t *obs, int i, int j, int k, const nav_t *nav)
{
double freq1, freq2, L1, L2; /* 频率 (Hz)、单差相位 (周期) */
/* 获取频率 */
freq1 = sat2freq(obs[i].sat, obs[i].code[0], nav); /* L1 频率 */
freq2 = sat2freq(obs[i].sat, obs[i].code[k], nav); /* 目标频率 (L2/L5) */
/* 计算单差相位 */
L1 = sdobs(obs, i, j, 0); /* L1 单差相位 */
L2 = sdobs(obs, i, j, k); /* 目标频率单差相位 */
/* 有效性检查 */
if (freq1 == 0.0 || freq2 == 0.0 || L1 == 0.0 || L2 == 0.0) return 0.0;
/* 计算 GF 组合 (单位:米) */
return L1 * CLIGHT / freq1 - L2 * CLIGHT / freq2;
}
9.7.7 mwmeas():求 MW 线性组合观测值
1. 基本原理
MW 线性组合(由 mwmeas()
计算):
其中:
- :L1 和 L2/L5 的载波相位(单位:周期)。
- :对应伪距(单位:米)。
- :L1 和 L2/L5 的频率(单位:Hz)。
- (光速,约 299792458 m/s)。
跳跃检测:
其中 是预定义阈值(单位:米,典型值 0.15-0.3 米)。
3. 参数列表
/* args */
const obsd_t *obs I 观测数据
const nav_t *nav I 导航数据
/* return */
double mw O MW 组合观测值
4. 源码解析
/* Melbourne-Wubbena linear combination --------------------------------------*/
/*
* 计算单个观测的 MW 线性组合。
* 参数:obs - 观测数据,nav - 导航数据。
* 输出:返回 MW 组合值 (单位:米),无效时返回 0.0。
*/
static double mwmeas(const obsd_t *obs, const nav_t *nav)
{
double freq1, freq2; /* L1 和 L2/L5 频率 (单位:Hz) */
freq1 = sat2freq(obs->sat, obs->code[0], nav); /* 获取 L1 频率 */
freq2 = sat2freq(obs->sat, obs->code[1], nav); /* 获取 L2/L5 频率 */
if (freq1 == 0.0 || freq2 == 0.0 || obs->L[0] == 0.0 || obs->L[1] == 0.0 ||
obs->P[0] == 0.0 || obs->P[1] == 0.0) return 0.0; /* 无效数据检查 */
trace(3, "mwmeas: %12.1f %12.1f %15.3f %15.3f %15.3f %15.3f %d %d\n",
freq1, freq2, obs->L[0], obs->L[1], obs->P[0], obs->P[1], obs->code[0], obs->code[1]); /* 调试:输出数据 */
/* 计算 MW 组合 */
return (obs->L[0] - obs->L[1]) * CLIGHT / (freq1 - freq2) -
(freq1 * obs->P[0] + freq2 * obs->P[1]) / (freq1 + freq2);
}