核心问题:TDCP 为什么能利用载波相位直接估计短时位移?它与测站间单差、星间单差有何不同?怎样测速、怎样防止周跳和模型错误?
符号约定:使用 ECEF 坐标;历元差一律定义为“当前历元 kkk 减前一历元 k−1k-1k−1”;时间间隔 Δt=tk−tk−1>0\Delta t=t_k-t_{k-1}>0Δt=tk−tk−1>0。所有载波相位公式以米为单位,除非明确注明 cycle。
1. TDCP 是什么?为什么需要它?
TDCP(Time-Differenced Carrier Phase)即历元间载波相位差分:对同一接收机、同一卫星、同一信号在两个历元的载波相位观测值作时间差。
ΔtLks=Lks−Lk−1s=λs(Φks−Φk−1s).(1) \Delta_t L_k^s=L_k^s-L_{k-1}^s =\lambda^s\left(\Phi_k^s-\Phi_{k-1}^s\right). \tag{1} ΔtLks=Lks−Lk−1s=λs(Φks−Φk−1s).(1)
- sss:卫星标识;k−1,kk-1,kk−1,k:前后两个历元。
- Φks\Phi_k^sΦks:接收机输出的载波相位观测值,单位 cycle(周)。
- λs=c/fs\lambda^s=c/f^sλs=c/fs:该信号载波波长,单位 m;ccc 为光速(m/s),fsf^sfs 为该信号频率(Hz)。
- Lks=λsΦksL_k^s=\lambda^s\Phi_k^sLks=λsΦks:载波相位距离观测量,单位 m。
- Δt\Delta_tΔt:时间差分算子,不是微分求导;整个时间间隔内必须使用同一频率和一致的相位定义。
核心痛点:载波相位观测精度高,但单历元绝对测距含有未知整周模糊度;如果只想知道短时间内移动了多少,不一定要先固定这个整数。在同一接收机、同一卫星、同一信号的连续跟踪弧内,模糊度不变,前后两历元作差就能消去它。[1][2]
这里的“连续”不只是文件中有两条记录:还要求接收机相位基准未发生不可处理的跳变,且时间标记、频率、波长和观测定义一致。失锁后重新捕获、周跳或设备重置都可能使这一前提失效。
因此,TDCP 主要用于短时位移与区间平均速度估计。它不是一种独立的绝对定位技术,也不能不加处理就把原始相位差除以时间当作接收机速度。
1.1 一句话区分常见方法
| 方法 | 直接利用的主要观测 | 主要结果 | 主要特点 |
|---|---|---|---|
| SPP | 伪距 | 绝对坐标、接收机钟差 | 易获得起始坐标,但位置噪声相对大 |
| Doppler | 多普勒频移 | 近似瞬时速度、钟漂 | 单历元可测速,较少依赖跨历元相位连续性 |
| TDCP | 两历元的连续载波相位 | 区间位移、平均速度、钟差变化 | 不需要先固定整数模糊度;极其敏感于周跳 |
| RTK / PPP-AR | 伪距、相位及必要的改正 | 高精度位置 | 需要相应的轨钟、大气、模糊度和/或基准站处理 |
2. 核心原理:TDCP 究竟消掉了什么?
以同一卫星 sss、同一频率的未组合载波相位为例,采用一致的符号约定:
Lks=ρks+c(δtr,k−δtks)+Tks−Iks+λsNks+Bks+εks.(2) L_k^s=\rho_k^s+c(\delta t_{r,k}-\delta t_k^s) +T_k^s-I_k^s+\lambda^s N_k^s+B_k^s+\varepsilon_k^s. \tag{2} Lks=ρks+c(δtr,k−δtks)+Tks−Iks+λsNks+Bks+εks.(2)
| 符号 | 含义 | 单位 |
|---|---|---|
| ρks\rho_k^sρks | 接收机与卫星的几何传播距离,计算时应考虑卫星发射时刻、地球自转等 | m |
| δtr,k\delta t_{r,k}δtr,k、δtks\delta t_k^sδtks | 接收机钟差、卫星钟差;须明确所采用的时间系统和偏差基准 | s |
| TksT_k^sTks | 对流层传播延迟 | m |
| IksI_k^sIks | 一阶电离层相位超前的距离绝对值;相位方程前为负号 | m |
| NksN_k^sNks | 整周模糊度,在相位连续跟踪弧内为常数 | cycle |
| BksB_k^sBks | 其余相位相关项,如相位缠绕、天线相位中心改正和硬件相位偏差 | m |
| εks\varepsilon_k^sεks | 多路径、热噪声与未建模误差 | m |
按“当前减前一历元”相减:
ΔtLks=Δtρks+cΔtδtr,k−cΔtδtks+ΔtTks−ΔtIks+λsΔtNks+ΔtBks+Δtεks.(3) \begin{aligned} \Delta_t L_k^s={}&\Delta_t\rho_k^s +c\Delta_t\delta t_{r,k}-c\Delta_t\delta t_k^s\\ &+\Delta_t T_k^s-\Delta_t I_k^s +\lambda^s\Delta_t N_k^s+\Delta_t B_k^s+\Delta_t\varepsilon_k^s. \end{aligned} \tag{3} ΔtLks=Δtρks+cΔtδtr,k−cΔtδtks+ΔtTks−ΔtIks+λsΔtNks+ΔtBks+Δtεks.(3)
只要该卫星、该信号在两历元之间保持相位连续且没有周跳,就有
ΔtNks=Nks−Nk−1s=0.(4) \boxed{\Delta_t N_k^s=N_k^s-N_{k-1}^s=0.} \tag{4} ΔtNks=Nks−Nk−1s=0.(4)
这一步严格消去的是保持不变的模糊度项;固定不变的硬件相位偏差也可在时间差中抵消,但随时间变化的偏差不会。卫星移动、卫星钟差变化、接收机钟差变化及大气变化不会自动消失。[1][3]
2.1 为什么“差分后精度更高”不是无条件成立?
TDCP 的高精度潜力来自载波相位本身的低噪声以及一些慢变化误差在短间隔内的相关性,不是因为做差能让所有噪声变小。若相邻两历元的相位白噪声独立且标准差均为 σL\sigma_LσL(m),则差分后为 2σL\sqrt2\sigma_L2σL,反而略大。第 8 章给出完整关系。
3. 时间单差、测站间单差、星间单差有何区别?
“单差”只是减法次数;关键在于沿哪个维度做差。 令 Lr,ksL_{r,k}^sLr,ks 表示测站 rrr、历元 kkk、卫星 sss 的同一信号载波相位距离观测值。
| 对比项目 | 时间单差(TDCP) | 测站间单差 | 星间单差 |
|---|---|---|---|
| 差分维度 | 时间 | 接收机/测站 | 卫星 |
| 差分对象 | 同站、同星、同信号,前后两个历元 | 同星、同历元、两个站 | 同站、同历元、两颗星 |
| 基本形式 | Lr,ks−Lr,k−1sL_{r,k}^{s}-L_{r,k-1}^{s}Lr,ks−Lr,k−1s | LB,ks−LA,ksL_{B,k}^{s}-L_{A,k}^{s}LB,ks−LA,ks | Lr,ks−Lr,kqL_{r,k}^{s}-L_{r,k}^{q}Lr,ks−Lr,kq |
| 主要消去 | 连续无周跳时不变的整周模糊度及静态偏差 | 在同步、同星及一致建模条件下抵消公共卫星钟差;近站还削弱共同传播误差 | 同一系统、同一时间基准下的公共接收机钟差 |
| 仍保留的钟差 | 接收机钟差变化、卫星钟差变化 | 两接收机钟差之差 | 不同卫星钟差之差 |
| 仍存在的模糊度 | 无周跳时消去 | 两接收机的模糊度之差 | 两卫星的模糊度之差 |
| 是否要求第二个站 | 不需要 | 需要 | 不需要 |
| 常见用途 | 短时位移、平均测速 | 相对定位、基线模型 | 消除公共钟差、构造双差 |
| 主要风险 | 周跳、参考位置与历元间模型不连续 | 同步误差、基线相关性 | 参考星变化、系统/频率偏差与差分相关性 |
其中 A,BA,BA,B 为两台接收机,s,qs,qs,q 为两颗卫星。差分不保证所有公共项严格为零:测站间单差通常假定两站时间同步并使用一致的卫星钟差模型;实际信号发射时刻不同,长基线下仍可能出现残差。星间单差要求所消除的接收机钟差确实为同一个参数。不同频率、跨星座或跨时间基准时,还须考虑硬件偏差和系统间偏差。[5]
3.1 双差与三差不等于 TDCP
站星双差:
∇sΔrLk=(LB,ks−LA,ks)−(LB,kq−LA,kq).(5) \nabla_s\Delta_r L_k =(L_{B,k}^{s}-L_{A,k}^{s})-(L_{B,k}^{q}-L_{A,k}^{q}). \tag{5} ∇sΔrLk=(LB,ks−LA,ks)−(LB,kq−LA,kq).(5)
站星时间三差:
Δt∇sΔrLk=∇sΔrLk−∇sΔrLk−1.(6) \Delta_t\nabla_s\Delta_r L_k =\nabla_s\Delta_r L_k-\nabla_s\Delta_r L_{k-1}. \tag{6} Δt∇sΔrLk=∇sΔrLk−∇sΔrLk−1.(6)
- Δr\Delta_rΔr:测站间差分算子(B−AB-AB−A)。
- ∇s\nabla_s∇s:星间差分算子(s−qs-qs−q)。
- Δt\Delta_tΔt:时间差分算子(k−(k−1)k-(k-1)k−(k−1))。
三差还包含测站差和卫星差;而本文的 TDCP 是单台接收机对同星同信号的时间单差。二者虽然都利用时间差消去稳定的模糊度,但不是同一个观测组合。
4. 为什么一颗卫星的 TDCP 还不能测速?
先修正已知的卫星钟差、大气延迟和其他相位效应:
L~ks=Lks+cδt^ks−T^ks+I^ks−B^ks.(7) \widetilde L_k^s =L_k^s+c\widehat{\delta t}_k^s-\widehat T_k^s +\widehat I_k^s-\widehat B_k^s. \tag{7} Lks=Lks+cδtks−Tks+Iks−Bks.(7)
式中上标帽号表示用于改正的模型值;L~\widetilde LL 为改正后相位距离观测值(m)。未完全改正的误差归入残差。于是两历元时间差满足
ΔtL~ks≃ρs(rk,k)−ρs(rk−1,k−1)+Δbk+ηks,(8) \Delta_t\widetilde L_k^s \simeq \rho^s(\mathbf r_k,k)-\rho^s(\mathbf r_{k-1},k-1) +\Delta b_k+\eta_k^s, \tag{8} ΔtLks≃ρs(rk,k)−ρs(rk−1,k−1)+Δbk+ηks,(8)
其中:
- rk−1\mathbf r_{k-1}rk−1、rk\mathbf r_krk:接收机前后两历元 ECEF 坐标(m)。
- ρs(r,k)\rho^s(\mathbf r,k)ρs(r,k):以接收机坐标 r\mathbf rr 和历元 kkk 为输入的完整几何距离模型(m)。
- bk=cδtr,kb_k=c\delta t_{r,k}bk=cδtr,k:接收机钟差的距离形式(m);Δbk=bk−bk−1\Delta b_k=b_k-b_{k-1}Δbk=bk−bk−1。
- ηks\eta_k^sηks:时间差中尚未建模的误差(m),假定未发生未检测周跳。
即使测站静止,卫星仍然在运动,所以式(8)中的几何距离变化不是接收机位移本身;且 Δbk\Delta b_kΔbk 对同一时间系统下各卫星共同存在。
4.1 先扣掉“测站不动时卫星造成的距离变化”
设前历元参考接收机坐标为 r0\mathbf r_0r0(m),未知位移为 Δr=rk−rk−1\Delta\mathbf r=\mathbf r_k-\mathbf r_{k-1}Δr=rk−rk−1,定义
gks=ρs(r0,k)−ρs(r0,k−1),yks=ΔtL~ks−gks.(9) \begin{aligned} g_k^s&=\rho^s(\mathbf r_0,k)-\rho^s(\mathbf r_0,k-1),\\ y_k^s&=\Delta_t\widetilde L_k^s-g_k^s. \end{aligned} \tag{9} gksyks=ρs(r0,k)−ρs(r0,k−1),=ΔtLks−gks.(9)
- gksg_k^sgks:接收机固定于 r0\mathbf r_0r0 时,由卫星运动等几何变化引起的理论距离差(m)。
- yksy_k^syks:扣除该理论几何差后的 TDCP 观测残差(m)。
若前历元真值等于参考坐标 rk−1=r0\mathbf r_{k-1}=\mathbf r_0rk−1=r0,在小位移条件下线性化:
yks≈−uksTΔr+Δbk+ηks.(10) \boxed{y_k^s\approx-\mathbf u_k^{sT}\Delta\mathbf r+\Delta b_k+\eta_k^s.} \tag{10} yks≈−uksTΔr+Δbk+ηks.(10)
其中 uks\mathbf u_k^suks 是当前历元从接收机指向卫星的视线单位向量(无量纲):
uks≈sks−r0∥sks−r0∥.(11) \mathbf u_k^s\approx\frac{\mathbf s_k^s-\mathbf r_0} {\|\mathbf s_k^s-\mathbf r_0\|}. \tag{11} uks≈∥sks−r0∥sks−r0.(11)
sks\mathbf s_k^ssks 为与距离模型一致、经过信号传播和地球自转处理的卫星位置(m),∥⋅∥\|\cdot\|∥⋅∥ 表示向量模长;−uksT-\mathbf u_k^{sT}−uksT 为几何距离对接收机坐标的偏导数(行向量)。式(11)是便于理解的近似几何表达。实际高精度实现中,距离计算要与卫星发射时刻、地球自转、天线相位中心等模型一致;雅可比也应是所采用距离模型对接收机坐标的导数。
一颗卫星只提供一个标量方程,但未知有 Δx,Δy,Δz,Δbk\Delta x,\Delta y,\Delta z,\Delta b_kΔx,Δy,Δz,Δbk 四项。因此需要多颗具有足够几何差异的卫星联合求解;“观测数足够”本身不代表几何满秩。
4.2 初始参考坐标不准会怎样?
如果前历元真实坐标是 rk−1=r0+e0\mathbf r_{k-1}=\mathbf r_0+\mathbf e_0rk−1=r0+e0,其中 e0\mathbf e_0e0 是真实位置减去参考位置(m),则一阶近似下额外出现
yks≈−uksTΔr+Δbk−(uks−uk−1s)Te0+ηks.(12) \boxed{y_k^s\approx-\mathbf u_k^{sT}\Delta\mathbf r+\Delta b_k -(\mathbf u_k^s-\mathbf u_{k-1}^s)^T\mathbf e_0+\eta_k^s.} \tag{12} yks≈−uksTΔr+Δbk−(uks−uk−1s)Te0+ηks.(12)
uk−1s\mathbf u_{k-1}^suk−1s 是前历元视线单位向量。可见,参考坐标误差通过两历元视线方向的差影响 TDCP;当时间间隔较短时,此项往往较小,但并不是严格为零。[7]
这解释了 TDCP 为什么通常能从 SPP 粗略坐标启动,又为什么在高精度、较长历元间隔或几何变化较大时,需要重视起始坐标质量。式(12)是一阶误差近似,并不意味着任何参考位置误差都可按这一式子精确扣除。[7]
5. 多星 TDCP:四未知数最小二乘如何求位移?
对同一时间系统下的 mmm 颗有效卫星建立观测方程:
y=Hx+η,x=[ΔxΔyΔzΔbk]T,(13) \mathbf y=\mathbf H\mathbf x+\boldsymbol\eta, \quad \mathbf x=\begin{bmatrix}\Delta x&\Delta y&\Delta z&\Delta b_k\end{bmatrix}^{T}, \tag{13} y=Hx+η,x=[ΔxΔyΔzΔbk]T,(13)
H=[−ux,k1−uy,k1−uz,k11−ux,k2−uy,k2−uz,k21⋮⋮⋮⋮−ux,km−uy,km−uz,km1].(14) \mathbf H=\begin{bmatrix} -u_{x,k}^{1}&-u_{y,k}^{1}&-u_{z,k}^{1}&1\\ -u_{x,k}^{2}&-u_{y,k}^{2}&-u_{z,k}^{2}&1\\ \vdots&\vdots&\vdots&\vdots\\ -u_{x,k}^{m}&-u_{y,k}^{m}&-u_{z,k}^{m}&1 \end{bmatrix}. \tag{14} H=−ux,k1−ux,k2⋮−ux,km−uy,k1−uy,k2⋮−uy,km−uz,k1−uz,k2⋮−uz,km11⋮1.(14)
- y\mathbf yy:由 yksy_k^syks 组成的 m×1m\times1m×1 向量(m)。
- x\mathbf xx:三个坐标位移分量(m)及接收机钟差距离变化(m)的 4×14\times14×1 向量。
- ux,ks,uy,ks,uz,ksu_{x,k}^s,u_{y,k}^s,u_{z,k}^sux,ks,uy,ks,uz,ks:视线单位向量在 ECEF 三轴的无量纲分量。
- H\mathbf HH:m×4m\times4m×4 设计矩阵;η\boldsymbol\etaη:m×1m\times1m×1 差分误差向量(m)。
加权最小二乘(WLS):
x^=argminx(y−Hx)TR−1(y−Hx).(15) \widehat{\mathbf x}= \arg\min_{\mathbf x}(\mathbf y-\mathbf H\mathbf x)^T \mathbf R^{-1}(\mathbf y-\mathbf H\mathbf x). \tag{15} x=argxmin(y−Hx)TR−1(y−Hx).(15)
R=Cov(η)\mathbf R=\operatorname{Cov}(\boldsymbol\eta)R=Cov(η) 为差分观测误差协方差矩阵(m²),R−1\mathbf R^{-1}R−1 为权阵(m⁻²)。若各卫星的差分观测误差可近似独立,可用对角阵 R=diag(σ12,…,σm2)\mathbf R=\operatorname{diag}(\sigma_1^2,\ldots,\sigma_m^2)R=diag(σ12,…,σm2);其中 σs\sigma_sσs 应当是 TDCP 差分后(不是原始单历元相位)的标准差。若使用的卫星钟差、轨道或其他改正引入跨星相关误差,则需评估 R\mathbf RR 的非对角项。
对于满秩、线性化合理、且 R\mathbf RR 正确描述观测误差的情况,可用下面的公式近似评估仅由观测噪声引起的位移与测速协方差:
Px≃(HTR−1H)−1,Pv≃[Px]1:3,1:3Δt2. \mathbf P_x\simeq(\mathbf H^T\mathbf R^{-1}\mathbf H)^{-1}, \qquad \mathbf P_v\simeq\frac{[\mathbf P_x]_{1:3,1:3}}{\Delta t^2}. Px≃(HTR−1H)−1,Pv≃Δt2[Px]1:3,1:3.
Px\mathbf P_xPx 为四未知数估计协方差矩阵(m²),[Px]1:3,1:3[\mathbf P_x]_{1:3,1:3}[Px]1:3,1:3 为其前三维位移的协方差子阵(m²),Pv\mathbf P_vPv 为区间平均速度协方差(m²/s²)。这只是条件协方差,未自动包含参考坐标误差、未建模轨钟/大气误差和周跳风险;不能把其对角元开方直接当作完整测速精度保证。
必要条件:至少四条几何独立的有效观测,且 H\mathbf HH 满列秩。要注意:
- 四条信号 ≠ 四个独立卫星方向;同一卫星的多频信号不能凭空增加三维位置的几何秩。
- 恰好四颗卫星、四个未知数时通常没有残差冗余,难以通过单次拟合残差发现周跳。
- 卫星几何不良时即使满秩,也会严重放大测速噪声;须检查条件数、残差和观测数量。
- GPS、Galileo、BDS 等混合时,若系统间偏差随时间变化,应增加对应的 ΔISB\Delta ISBΔISB 参数,或使用经过验证的偏差校准;不能一概使用同一个四参数模型。GLONASS FDMA 与不同信号的硬件偏差也需要单独处理。
6. TDCP 如何测速?与 Doppler 有何不同?
TDCP 直接估计的是位移 Δr^\widehat{\Delta\mathbf r}Δr(m),除以历元间隔得到
v‾k−1:k=Δr^Δt,Δt=tk−tk−1>0.(16) \boxed{\overline{\mathbf v}_{k-1:k} =\frac{\widehat{\Delta\mathbf r}}{\Delta t}}, \qquad \Delta t=t_k-t_{k-1}>0. \tag{16} vk−1:k=ΔtΔr,Δt=tk−tk−1>0.(16)
v‾k−1:k\overline{\mathbf v}_{k-1:k}vk−1:k 是两历元之间的 ECEF 平均速度(m/s)。它不能在存在明显加速度或转弯时直接等同于末端时刻 tkt_ktk 的瞬时速度。若需要比较 TDCP 和 Doppler 的精度,须对齐两者的速度定义和参考时刻;在速度随时间近似线性变化的条件下,TDCP 平均速度更接近时间区间中点的瞬时速度,而不是终点速度。
Doppler 则利用载波相位的变化率。在接收机采用 L˙≈−λD\dot L\approx-\lambda DL˙≈−λD 的多普勒符号约定时,有近似关系
Lks−Lk−1s≈−∫tk−1tkλsDs(t) dt.(17) L_k^s-L_{k-1}^s\approx -\int_{t_{k-1}}^{t_k}\lambda^s D^s(t)\,dt. \tag{17} Lks−Lk−1s≈−∫tk−1tkλsDs(t)dt.(17)
Ds(t)D^s(t)Ds(t) 为多普勒观测(Hz,即 cycle/s),ttt 为积分时间变量,λs\lambda^sλs 为波长。不同接收机的多普勒符号、平滑与采样定义可能不同;式(17)说明物理联系,不表示设备输出的多普勒积分一定与相位差逐位相等。[4]
| 比较维度 | TDCP 测速 | Doppler 测速 | 两次 SPP 位置相减测速 |
|---|---|---|---|
| 核心观测 | 两历元载波相位差 | 单历元多普勒 | 两历元伪距位置 |
| 输出含义 | 区间平均速度 | 近似瞬时速度 | 区间平均速度 |
| 无需模糊度固定 | 是(无周跳前提下) | 是 | 是 |
| 依赖跨历元相位连续 | 强依赖 | 相对较少 | 不依赖相位 |
| 精度潜力 | 在优良环境中高,可达 mm/s 量级 | 常见为 cm/s 量级,视设备和环境而定 | 短时间间隔时往往较差 |
| 高动态响应 | 区间平均会掩盖短时加速度 | 通常更适合快速变化 | 受粗位置噪声强烈影响 |
| 主要风险 | 周跳、星历/钟差切换、初值与差分误差 | Doppler 噪声、跟踪环路与钟漂模型 | SPP 噪声被时间差放大 |
| 合适角色 | 高精度短时位移与测速 | 稳健测速与 TDCP 回退 | 粗略估计或比较基线 |
输出导航坐标系速度时,可将 ECEF 平均速度旋转到当地 ENU:
v‾ENU=RENU←ECEF(φ,ℓ)v‾ECEF. \overline{\mathbf v}_{\mathrm{ENU}}= \mathbf R_{\mathrm{ENU}\leftarrow\mathrm{ECEF}}(\varphi,\ell) \overline{\mathbf v}_{\mathrm{ECEF}}. vENU=RENU←ECEF(φ,ℓ)vECEF.
RENU←ECEF\mathbf R_{\mathrm{ENU}\leftarrow\mathrm{ECEF}}RENU←ECEF 为由测站参考大地纬度 φ\varphiφ、经度 ℓ\ellℓ 确定的 ECEF 到东-北-天(ENU)坐标旋转矩阵;其仅改变表达坐标系,不改变速度是区间平均量的含义。
精度量级只是特定环境下的文献结果,不能保证 TDCP 在所有车载、手机或遮挡场景都优于 Doppler。 例如 Ji 等(2022)在动态加减速场景下观察到 Doppler 更适合跟踪速度快速变化。[2][4][6]
7. 一个直观数字例子
设相邻历元间隔 Δt=1 s\Delta t=1\,\mathrm sΔt=1s。对一颗卫星:
- 接收机不动时,卫星自身运动使理论几何距离增加 600.000 m600.000\,\mathrm m600.000m;
- 完成卫星钟差、大气和其他相位改正后的相位距离增加 600.040 m600.040\,\mathrm m600.040m。
则扣除卫星运动后的观测为
yks=600.040−600.000=0.040 m.(18) y_k^s=600.040-600.000=0.040\,\mathrm m. \tag{18} yks=600.040−600.000=0.040m.(18)
假设该卫星当前视线方向 u=(0.6,0.8,0)T\mathbf u=(0.6,0.8,0)^Tu=(0.6,0.8,0)T,真实位移 Δr=(0.05,0,0)T m\Delta\mathbf r=(0.05,0,0)^T\,\mathrm mΔr=(0.05,0,0)Tm,而接收机钟差距离变化 Δb=0.070 m\Delta b=0.070\,\mathrm mΔb=0.070m,则
−uTΔr+Δb=−0.6×0.05+0.070=0.040 m.(19) -\mathbf u^T\Delta\mathbf r+\Delta b =-0.6\times0.05+0.070=0.040\,\mathrm m. \tag{19} −uTΔr+Δb=−0.6×0.05+0.070=0.040m.(19)
一颗卫星只能观测到“位移沿视线的投影 + 公共钟差变化”,所以必须使用多个不同方向的卫星;单星 0.0400.0400.040 m 并不等于设备运动了 0.0400.0400.040 m。
8. 误差、周跳与协方差:TDCP 的真正工程难点
8.1 时间差分不会自动消除噪声
设某条卫星相位噪声 εk\varepsilon_kεk 在相邻历元的方差分别为 σk2,σk−12\sigma_k^2,\sigma_{k-1}^2σk2,σk−12,则
Var(εk−εk−1)=σk2+σk−12−2Cov(εk,εk−1).(20) \operatorname{Var}(\varepsilon_k-\varepsilon_{k-1}) =\sigma_k^2+\sigma_{k-1}^2 -2\operatorname{Cov}(\varepsilon_k,\varepsilon_{k-1}). \tag{20} Var(εk−εk−1)=σk2+σk−12−2Cov(εk,εk−1).(20)
若独立且均为 σL\sigma_LσL:
σΔL=2σL,σLOS-speed≃2σLΔt.(21) \sigma_{\Delta L}=\sqrt{2}\sigma_L, \qquad \sigma_{\text{LOS-speed}}\simeq\frac{\sqrt2\sigma_L}{\Delta t}. \tag{21} σΔL=2σL,σLOS-speed≃Δt2σL.(21)
这里 σLOS-speed\sigma_{\text{LOS-speed}}σLOS-speed 是单条观测的速度噪声尺度,不是三维速度解的标准差。若 σL=3\sigma_L=3σL=3 mm、Δt=1\Delta t=1Δt=1 s,则为约 4.24 mm/s;真正的三维速度精度还取决于卫星几何、钟差参数和其他误差。
采样间隔并非越小越好:对于相邻历元独立的相位白噪声,间隔缩短会增大折算成速度后的噪声;间隔过长又容易增大大气变化、参考坐标误差、周跳和模型不一致的影响。若噪声具有时间相关性,需使用式(20)而不能只套用 2σL/Δt\sqrt 2\sigma_L/\Delta t2σL/Δt。
8.2 连续 TDCP 之间也可能相关
若三个历元原始相位白噪声相互独立且方差为 σL2\sigma_L^2σL2,定义
ηk=εk−εk−1,ηk+1=εk+1−εk, \eta_k=\varepsilon_k-\varepsilon_{k-1}, \qquad \eta_{k+1}=\varepsilon_{k+1}-\varepsilon_k, ηk=εk−εk−1,ηk+1=εk+1−εk,
则
Cov(ηk,ηk+1)=−σL2.(22) \operatorname{Cov}(\eta_k,\eta_{k+1})=-\sigma_L^2. \tag{22} Cov(ηk,ηk+1)=−σL2.(22)
原因是两次 TDCP 共用了中间历元的相位噪声。把连续 TDCP 放进滑窗平差、序贯估计或与 IMU 融合时,不应不加判断地将其当成完全独立的观测。这一结论直接来自上面的协方差恒等式。
8.3 周跳为何致命?
发生一周周跳时,差分式中出现 λsΔNs\lambda^s\Delta N^sλsΔNs:
δLslips=λsΔNs.(23) \delta L_{\text{slip}}^s=\lambda^s\Delta N^s. \tag{23} δLslips=λsΔNs.(23)
- ΔNs\Delta N^sΔNs:前后历元整数模糊度之差(cycle)。
- δLslips\delta L_{\text{slip}}^sδLslips:周跳引入的额外相位距离差(m)。
GPS L1 的波长约为 0.19030.19030.1903 m;在 1 s 间隔内,单周周跳对应某条 TDCP 观测约 0.19030.19030.1903 m/s 的伪速度尺度,但不是最终三维速度必然偏差 0.19030.19030.1903 m/s。故必须检查 lock time、LLI、半周标志、相位创新、多频几何无关组合及 Doppler 一致性等;单一阈值不能保证检测出所有小周跳。[3]
8.4 还有哪些误差不能忽略?
| 误差或异常 | 为什么影响 TDCP | 工程措施 |
|---|---|---|
| 卫星运动与地球自转 | 两历元几何距离本来就变 | 每历元按卫星发射时刻、地球自转一致计算 |
| 卫星钟差与星历切换 | 轨钟模型切换可形成不连续台阶 | 检查 IOD、跨历元轨钟连续性;必要时采用一致星历段或平滑处理 |
| 接收机钟差跳变 | 公共钟差变化很大;设备时间/相位重置还可能破坏相位连续性 | 联合估计 Δb\Delta bΔb,并检查接收机钟差重置、时间标签和相位基准一致性 |
| 电离层和对流层变化 | 时间差并非严格为零 | 按频率、采样间隔、精度需求建模或合理定权 |
| 多路径、低高度角与遮挡 | 差分后仍有显著变化 | 仰角/SNR 定权、抗差检验、剔除失锁信号 |
| 天线相位缠绕/安装姿态变化 | 相位模型随姿态、卫星方向变化 | 连续一致地处理相位缠绕、天线 PCO/PCV |
| 多系统、频率硬件偏差 | 不能保证共用同一钟差变化 | 增加 ΔISB\Delta ISBΔISB 等参数,或经验证后合并 |
| 起始坐标误差 | 两历元视线方向差会放大参考误差影响 | 提供可靠的 SPP/PPP 起点,并按精度需求评价误差传递 |
Freda 等(2015)特别指出:用不同广播星历参数集计算相邻历元轨道/钟差,可能引起 TDCP 伪不连续,应作为工程重点检查。[2]
9. TDCP 的优势与局限:选型比较
| 维度 | 优势 | 局限与代价 | 实际建议 |
|---|---|---|---|
| 模糊度 | 无周跳时直接消去不变的整周数 | 只要发生未检测周跳,估计就可能严重失真 | 周跳检测是首要条件 |
| 位移精度 | 可利用毫米级载波相位观测变化 | 仍受轨钟、几何、初值和大气变化影响 | 以真实数据验证有效精度 |
| 测速 | 短时平均速度精度潜力高 | 不是末端瞬时速度;小间隔白噪声会放大 | 需要高动态时结合 Doppler/IMU |
| 单机独立性 | 不要求本地基准站,也无需先固定模糊度 | 仍需卫星星历、钟差和必要改正模型 | 适合独立接收机辅助运动测量 |
| 启动 | 不必等待模糊度固定 | 首个历元无法形成时间差 | 先用 SPP 提供参考坐标 |
| 遮挡环境 | 良好跟踪时仍可输出细致位移信息 | 共同卫星减少、失锁和周跳频繁 | 不足时回退 Doppler,不强行输出 |
| 长时间导航 | 每一步都能提供短时运动约束 | 连续积分会积累漂移 | 定期用绝对位置约束校正 |
| 软件实现 | 最小二乘核心并不复杂 | 数据预处理与模型一致性要求高 | 优先做好观测匹配与质量控制 |
最重要的权衡:TDCP 省去了整数模糊度固定,但要求相位连续、前后历元模型一致,并需要控制短时误差与周跳。
10. 从观测到测速:推荐工程流程
最小工程清单:
- 按
constellation + satellite + observation signal + time tag匹配同一跟踪弧;相位须统一为米,时间差须来自一致的时间系统。即使两历元 RINEX 数据都存在,也应检查锁定与相位连续性。 - 周跳预检与异常剔除先于位移求解;失锁、半周不确定和频率切换不能被当成正常连续相位。
- 卫星发射时刻、地球自转、星历与钟差应逐历元计算,同时保证模型和偏差基准一致;遇到广播星历段、精密产品或钟差基准切换时,检查是否引入伪不连续。
- 使用相同符号约定处理电离层、对流层、相位缠绕、相位中心与时间系统。
- 检查有效独立卫星数、设计矩阵秩和条件数;有冗余观测时再做残差/粗差检验。仅四颗独立卫星的四参数解通常没有残差冗余,不能靠拟合残差证明“无周跳”。
- 输出明确标注为“k−1k-1k−1 到 kkk 的区间平均速度”;若要当前瞬时速度,需要额外模型和时间标记约定。
11. C++11 + Eigen 核心示例:迭代求解 TDCP 位移
以下代码与第 4 章的模型一致,采用非线性几何距离 + 迭代最小二乘,而不是仅做一次线性化。它只承担 TDCP 估计核心,不负责 RINEX/RTCM 解码、同星同信号匹配、周跳检测和完整 GNSS 改正。
使用前提:上游已提供两历元同一跟踪弧的改正后相位距离观测值(m),以及按一致时标和轨钟模型计算的卫星位置。代码仅用欧氏距离教学展示;真正高精度实现须将几何距离及其偏导一起替换为考虑信号发射时刻、地球自转等因素的一致模型。这里的 sigma_tdcp_m 是差分后的标准差,不是单历元相位标准差。
#include <Eigen/Dense>
#include <vector>
#include <stdexcept>
#include <cmath>
struct TdcpObs {
Eigen::Vector3d sat_prev_m; // 前历元卫星位置,m
Eigen::Vector3d sat_curr_m; // 当前历元卫星位置,m
double phase_prev_m; // 前历元已改正的相位距离,m
double phase_curr_m; // 当前历元已改正的相位距离,m
double sigma_tdcp_m; // 差分观测标准差,m
};
struct TdcpResult {
Eigen::Vector3d dr_m; // 前历元到当前历元的 ECEF 位移
Eigen::Vector3d vel_mps; // 区间 ECEF 平均速度
double db_m; // 接收机钟差的距离变化
double weighted_res_norm; // 无量纲加权残差二范数
};
TdcpResult solveTdcp(const std::vector<TdcpObs>& obs,
const Eigen::Vector3d& receiver_prev_m,
double dt_s)
{
if (obs.size() < 4 || !(dt_s > 0.0) || !std::isfinite(dt_s) ||
!receiver_prev_m.allFinite()) {
throw std::invalid_argument("invalid TDCP input");
}
Eigen::Vector4d x = Eigen::Vector4d::Zero(); // [dx,dy,dz,db]
bool converged = false;
Eigen::MatrixXd A(obs.size(), 4);
Eigen::VectorXd y(obs.size());
// 迭代:接收机前历元坐标视为已给定参考值。
for (int iter = 0; iter < 6; ++iter) {
const Eigen::Vector3d receiver_curr_m = receiver_prev_m + x.head<3>();
for (std::size_t i = 0; i < obs.size(); ++i) {
const TdcpObs& o = obs[i];
if (!(o.sigma_tdcp_m > 0.0) || !std::isfinite(o.sigma_tdcp_m) ||
!o.sat_prev_m.allFinite() || !o.sat_curr_m.allFinite() ||
!std::isfinite(o.phase_prev_m) || !std::isfinite(o.phase_curr_m)) {
throw std::invalid_argument("invalid TDCP observation");
}
const Eigen::Vector3d d_prev = o.sat_prev_m - receiver_prev_m;
const Eigen::Vector3d d_curr = o.sat_curr_m - receiver_curr_m;
const double rho_prev = d_prev.norm();
const double rho_curr = d_curr.norm();
if (!(rho_prev > 0.0 && rho_curr > 0.0)) {
throw std::runtime_error("invalid satellite range");
}
const Eigen::Vector3d los = d_curr / rho_curr;
const double tdcp_m = o.phase_curr_m - o.phase_prev_m;
const double pred_m = rho_curr - rho_prev + x(3);
const double w = 1.0 / o.sigma_tdcp_m;
A.row(i).head<3>() = (-w * los).transpose();
A(i, 3) = w;
y(i) = w * (tdcp_m - pred_m);
}
// SVD 既用于检查几何条件,也避免形成 H^T W H 的平方条件数。
Eigen::JacobiSVD<Eigen::MatrixXd> svd(
A, Eigen::ComputeThinU | Eigen::ComputeThinV);
const Eigen::VectorXd singular = svd.singularValues();
if (!(singular(3) > 0.0) ||
singular(0) / singular(3) > 1e6) { // 演示阈值,应按业务验证
throw std::runtime_error("poor TDCP geometry");
}
const Eigen::Vector4d step = svd.solve(y);
if (!step.allFinite()) {
throw std::runtime_error("TDCP solve failed");
}
x += step;
if (step.norm() < 1e-5) { // 教学阈值,单位 m
converged = true;
break;
}
}
if (!converged) {
throw std::runtime_error("TDCP iteration did not converge");
}
// 以最终解重新计算残差,避免使用上一次线性化的残差。
double sum_sq = 0.0;
const Eigen::Vector3d r_curr = receiver_prev_m + x.head<3>();
for (const TdcpObs& o : obs) {
const double model = (o.sat_curr_m - r_curr).norm()
- (o.sat_prev_m - receiver_prev_m).norm() + x(3);
const double residual = (o.phase_curr_m - o.phase_prev_m) - model;
const double scaled = residual / o.sigma_tdcp_m;
sum_sq += scaled * scaled;
}
TdcpResult out;
out.dr_m = x.head<3>();
out.vel_mps = out.dr_m / dt_s;
out.db_m = x(3);
out.weighted_res_norm = std::sqrt(sum_sq);
return out;
}
11.1 代码与数学公式对照
这段实现实际上求解的是下面的非线性加权最小二乘,而不是简单“原始相位差除以时间”:
minΔr,Δb∑s=1m[ΔtL~ks−{ρs(r0+Δr,k)−ρs(r0,k−1)+Δb}σs]2.(24) \boxed{ \min_{\Delta\mathbf r,\Delta b} \sum_{s=1}^{m} \left[\frac{ \Delta_t\widetilde L_k^s -\left\{\rho^s(\mathbf r_0+\Delta\mathbf r,k) -\rho^s(\mathbf r_0,k-1)+\Delta b\right\} }{\sigma_s}\right]^2.} \tag{24} Δr,Δbmins=1∑m[σsΔtLks−{ρs(r0+Δr,k)−ρs(r0,k−1)+Δb}]2.(24)
- r0\mathbf r_0r0:前历元接收机参考坐标(m)。
- Δr\Delta\mathbf rΔr:待求的前后历元位移(m)。
- Δb\Delta bΔb:待求的接收机钟差距离变化(m)。
- ρs(⋅,k)\rho^s(\cdot,k)ρs(⋅,k):当前历元卫星 sss 的几何距离模型(m)。
- σs\sigma_sσs:第 sss 条 TDCP 观测的标准差(m)。
迭代时定义预测值、观测减计算残差和雅可比行:
fs(x)=ρs(r0+Δr,k)−ρs(r0,k−1)+Δb,vs=ΔtL~ks−fs(x),hs=[−uksT, 1],x=[Δx,Δy,Δz,Δb]T.(25) \begin{aligned} f_s(\mathbf x)&=\rho^s(\mathbf r_0+\Delta\mathbf r,k) -\rho^s(\mathbf r_0,k-1)+\Delta b,\\ v_s&=\Delta_t\widetilde L_k^s-f_s(\mathbf x),\\ \mathbf h_s&=\left[-\mathbf u_k^{sT},\ 1\right],\\ \mathbf x&=[\Delta x,\Delta y,\Delta z,\Delta b]^T. \end{aligned} \tag{25} fs(x)vshsx=ρs(r0+Δr,k)−ρs(r0,k−1)+Δb,=ΔtLks−fs(x),=[−uksT, 1],=[Δx,Δy,Δz,Δb]T.(25)
其中 fsf_sfs 为模型预测(m);vsv_svs 为残差(m);hs\mathbf h_shs 是 1×41\times41×4 雅可比行;x\mathbf xx 为当前迭代的未知参数向量(m)。将 vsv_svs 和 hs\mathbf h_shs 同乘以 1/σs1/\sigma_s1/σs,再用 SVD 求参数增量 δx\delta\mathbf xδx:
minδx∥Aδx−y∥22,x←x+δx.(26) \min_{\delta\mathbf x} \left\|\mathbf A\delta\mathbf x-\mathbf y\right\|_2^2, \qquad \mathbf x\leftarrow\mathbf x+\delta\mathbf x. \tag{26} δxmin∥Aδx−y∥22,x←x+δx.(26)
A\mathbf AA、y\mathbf yy 分别为按 1/σs1/\sigma_s1/σs 加权后堆叠得到的设计矩阵和残差向量;∥⋅∥2\|\cdot\|_2∥⋅∥2 表示欧氏二范数。
必须理解代码边界:
-
示例的迭代未知数只有当前位移和接收机钟差增量;将前历元参考坐标当作已知。因此它没有联合估计前历元位置不确定度,代码计算得到的残差也不能反映由参考坐标误差引起的全部偏差。
-
示例中的
sat_prev_m、sat_curr_m和欧氏距离只是一致几何的教学近似,不是完整的信号传播模型;真实工程需替换为同一约定下的range()/雅可比实现。 -
phase_prev_m、phase_curr_m必须按同一卫星钟差及硬件偏差基准和同一信号处理,单位均为米。不能把 RINEX 相位周数直接作为米输入。未改正的原始相位并非在数学上完全不可用,但此代码没有为它们建立相应的卫星钟差、大气及其他误差项。 -
代码默认一个时间系统、一个公共接收机钟差变化参数;未处理多星座 ISB、GLONASS FDMA、周跳、观测间相关和抗差粗差识别。
-
weighted_res_norm只是残差指标,不代表测速精度,也不是自动周跳判定;当只有 4 条独立观测时,残差本身缺少冗余检核能力。 -
1e6条件数阈值、1e-5m 停止阈值是教学示例数值,不能直接作为任何接收机产品的通用门限;工程上还应配置最大迭代、收敛质量标记和异常数据回退。
12. TDCP 在 GNSS 工程中怎样使用?
| 场景 | 推荐用法 | 关键注意 |
|---|---|---|
| 车载 GNSS/机器人 | 提供短时位移和平均速度增量 | 遮挡、周跳、高动态时间差 |
| 静态基准站 | 检测异常位移、短时间稳定性 | 仍需消除卫星运动及钟差变化 |
| 结构形变监测 | 辅助高频微位移、振动分析 | 多路径和高采样率速度噪声 |
| GNSS/IMU 紧组合或松组合 | 提供时间区间位移约束 | 时间同步、天线杠杆臂、跨历元噪声相关 |
| SPP 启动后的测速 | 用 SPP 作参考坐标,TDCP 补充高精度运动 | 需评价起始坐标误差传播 |
| PPP/RTK 暂时不稳定时 | 维持短时运动估计辅助信息 | 不替代绝对定位;漂移逐步累积 |
12.1 TDCP 是否应完全替代 Doppler?
通常不建议。 连续跟踪、周跳检核通过、卫星数量足够且模型一致时,TDCP 可以提供很精细的短时位移;首次启动、严重遮挡、共同卫星不足、星历不连续或动态剧变时,则应保留 Doppler 结果作为可用性与异常检测的补充。[4][6]
12.2 TDCP 能否提供绝对坐标?
TDCP 自身只约束位置变化。若上一历元坐标有参考值,可以做
r^k=r^k−1+Δr^k.(27) \widehat{\mathbf r}_k=\widehat{\mathbf r}_{k-1} +\widehat{\Delta\mathbf r}_k. \tag{27} rk=rk−1+Δrk.(27)
r^k\widehat{\mathbf r}_krk 是当前 ECEF 位置估计(m),r^k−1\widehat{\mathbf r}_{k-1}rk−1 是前历元位置估计,Δr^k\widehat{\Delta\mathbf r}_kΔrk 是本次 TDCP 位移估计。这种逐历元累加不是无漂移的绝对定位;应定期使用 SPP、PPP 或 RTK 的绝对位置约束修正。
13. 初学者最容易混淆的概念
| 常见说法 | 准确理解 |
|---|---|
| TDCP 是“测站间单差” | 不是;它是时间单差,不需要第二台接收机 |
| TDCP 等于 RTK 三差 | 不是;三差还包含测站差与卫星差 |
| TDCP 必须固定整周模糊度 | 连续无周跳时不必固定;模糊度差为零 |
| 差分可以消掉所有误差 | 只能严格消去相同的公共常数;时变误差仍存在 |
| 原始相位差除以时间就是速度 | 错;卫星运动和两类钟差必须正确处理 |
| TDCP 测得时刻 kkk 的瞬时速度 | 直接得到的是**k−1k-1k−1 至 kkk 的平均速度** |
| 四个信号就够了 | 必须有四个几何独立的观测方向,且质量合格 |
| 参考坐标精度无所谓 | 会通过两历元视线方向差影响位移估计 |
| TDCP 一定优于 Doppler | 无此保证;TDCP 受周跳、遮挡、高动态和模型不连续限制 |
| TDCP 可长期积分得到无漂移坐标 | 不行;累积误差必须由绝对定位约束修正 |
14. 参考文献
以下资料依照出版卷期年份标注(在线优先发表时间可能更早)。[1] 支撑载波相位观测模型;[3] 支撑周跳检测;[5] 用于介绍 RTK 测站—卫星双差;[2][4][6] 为 TDCP 测速及方法比较;[7] 直接研究首历元坐标误差的影响。本文的简化推导、教学数值算例和程序门限,不应误写为论文原文结论。
[1] Sanz Subirana, J., Juan Zornoza, J. M., Hernández-Pajares, M. (2011). GNSS Basic Observables. ESA Navipedia. https://gssc.esa.int/navipedia/index.php/GNSS_Basic_Observables
[2] Freda, P., Angrisano, A., Gaglione, S., Troisi, S. (2015). Time-differenced carrier phases technique for precise GNSS velocity estimation. GPS Solutions, 19, 335–341. https://doi.org/10.1007/s10291-014-0425-1 (2014 年在线发表,2015 年卷期出版)
[3] Sanz Subirana, J., Juan Zornoza, J. M., Hernández-Pajares, M. (2011). Carrier Phase Cycle-Slip Detection. ESA Navipedia. https://gssc.esa.int/navipedia/index.php/Carrier_Phase_Cycle-Slip_Detection
[4] Angrisano, A., Cappello, G., Gaglione, S., Gioia, C. (2024). Velocity Estimation Using Time-Differenced Carrier Phase and Doppler Shift with Different Grades of Devices: From Smartphones to Professional Receivers. Algorithms, 17(1), 2. https://doi.org/10.3390/a17010002 (2023 年 12 月在线发表,归属 2024 年卷期)
[5] GMV (ed.). (2011). RTK Fundamentals. ESA Navipedia. https://gssc.esa.int/navipedia/index.php/RTK_Fundamentals
[6] Ji, L., Sun, R., Cheng, Q., Wang, J. (2022). Evaluation of the performance of GNSS-based velocity estimation algorithms. Satellite Navigation, 3, Article 18. https://doi.org/10.1186/s43020-022-00080-4
[7] Cherfi, H., Lesouple, J., Solà, J., Thevenon, P. (2025). Uncertainty Quantification of First Fix in a Time-Differenced Carrier Phase Observation Model. Sensors, 25(11), 3480. https://doi.org/10.3390/s25113480
结论:TDCP 在连续跟踪、无周跳条件下,用时间差消去不变的整周模糊度;正确处理卫星运动、卫星钟差与其他时变误差,并估计接收机钟差变化后,可以求短时三维位移,进而得到区间平均速度。TDCP 的优势来自载波相位的低噪声,而非“差分后噪声必然更小”。



被折叠的 条评论
为什么被折叠?



