E.7 动态,静态和移动基线
主要是RTKLIB ver. 2.4.2 Manual E.7的翻译.
目录
RTKLIB利用EKF(扩展卡尔曼滤波),结合附录E.3中的GNSS信号测量模型和附录E.5中的对流层和电离层模型,获得了DGPS/DGNSS、静态、运动和移动基线模式下的最终解。
(1) EKF公式
利用EKF,可以用测量向量yk在历元tk处估计未知模型参数的状态向量x及其协方差矩阵P:
x
^
k
(
+
)
=
x
^
k
(
−
)
+
K
k
(
y
k
−
h
(
x
^
k
(
−
)
)
)
−
−
−
(
E
.
7.1
)
P
k
(
+
)
=
(
I
−
K
k
H
(
x
^
k
(
−
)
)
)
P
k
(
−
)
−
−
−
(
E
.
7.2
)
K
k
=
P
k
(
−
)
H
(
x
^
k
(
−
)
)
(
H
(
x
^
k
(
−
)
)
P
k
(
−
)
H
(
x
^
k
(
−
)
)
T
+
R
k
)
−
1
−
−
−
(
E
.
7.3
)
\begin{aligned} \hat x_k(+)&=\hat x_k(-)+K_k(y_k-h(\hat x_k(-)))---(E.7.1)\\ P_k(+)&=(I-K_kH(\hat x_k(-)))P_k(-)---(E.7.2)\\ K_k&=P_k(-)H(\hat x_k(-))(H(\hat x_k(-))P_k(-)H(\hat x_k(-))^T+R_k)^{-1}---(E.7.3) \end{aligned}
x^k(+)Pk(+)Kk=x^k(−)+Kk(yk−h(x^k(−)))−−−(E.7.1)=(I−KkH(x^k(−)))Pk(−)−−−(E.7.2)=Pk(−)H(x^k(−))(H(x^k(−))Pk(−)H(x^k(−))T+Rk)−1−−−(E.7.3)
其中
x
^
k
\hat x_k
x^k和Pk为在历元时间tk时的估计状态向量及其协方差矩阵。(-)和(+)表示EKF测量更新前后。
h
(
x
)
,
H
(
x
)
a
n
d
R
k
h(x),H(x)andR_k
h(x),H(x)andRk分别为测量模型向量、偏导数矩阵和测量误差协方差矩阵。假设系统模型是线性的,EKF的状态向量及其协方差矩阵的时间更新表示为:
x
^
k
+
1
(
−
)
=
F
k
k
+
1
x
^
k
(
+
)
−
−
−
(
E
.
7.4
)
P
k
+
1
(
−
)
=
F
k
k
+
1
P
k
(
+
)
F
k
k
+
1
T
+
Q
k
k
+
1
−
−
−
(
E
.
7.5
)
\begin{aligned} \hat x_{k+1}(-)&=F_k^{k+1}\hat x_k(+)---(E.7.4)\\ P_{k+1}(-)&=F_k^{k+1}P_k(+)F_k^{k+1T}+Q_k^{k+1}---(E.7.5) \end{aligned}
x^k+1(−)Pk+1(−)=Fkk+1x^k(+)−−−(E.7.4)=Fkk+1Pk(+)Fkk+1T+Qkk+1−−−(E.7.5)
其中,
F
k
k
+
1
F_k^{k+1}
Fkk+1和
Q
k
k
+
1
Q_k^{k+1}
Qkk+1是系统噪声从历时tk到tk+1的转移矩阵和协方差矩阵。
(2)短基线DD(双差)测量模型
对于探测器r和基站b之间的短基线(< 10 km)载波相对定位,通常使用以下DD测量方程来计算iL相位范围和伪距。在这些方程中,利用DD技术几乎消除了卫星和接收机的时钟偏差、电离层和对流层效应和其他小的修正项。
Φ
r
b
,
i
j
k
=
ρ
r
b
j
k
+
λ
i
(
B
r
b
,
i
j
−
B
r
b
,
i
k
)
+
d
Φ
r
,
i
s
+
ε
Φ
P
r
b
,
i
j
k
=
ρ
r
b
j
k
+
ε
P
−
−
−
(
E
.
7.6
)
\begin{aligned} \Phi_{rb,i}^{jk}&=\rho_{rb}^{jk}+\lambda_i(B_{rb,i}^{j}-B_{rb,i}^{k})+d\Phi_{r,i}^{s}+\varepsilon_\Phi\\ P_{rb,i}^{jk}&=\rho_{rb}^{jk}+\varepsilon_P---(E.7.6) \end{aligned}
Φrb,ijkPrb,ijk=ρrbjk+λi(Brb,ij−Brb,ik)+dΦr,is+εΦ=ρrbjk+εP−−−(E.7.6)
其中,
d
Φ
r
,
i
s
d\Phi_{r,i}^{s}
dΦr,is是载波相位校正项,在短基线情况下可以忽略,除了带有不同天线的接收机PCV项。为了得到方程中的几何范围
ρ
r
s
\rho_r^s
ρrs,除移动基线情况外,将基站位置
r
b
r_b
rb固定为预定值。
请注意,接收机之间的SD最好是在同一历元时间的测量之间进行。然而,由于接收机的时钟偏差不同,接收机并不是完全同步的。在某些典型情况下,移动站的采样间隔与基站不同,如10hz和1hz。为了控制SD, RTKLIB采用一个简单的标准来选择一个测量对。RTKLIB只是选择在探测器测量历元时间之前或等于这个时间的最后一次测量。移动站和基站之间的历时差有时被称为ʺ差异年龄ʺ。随着时间差的增大,由于卫星时钟漂移和电离层时延的变化,求解精度逐渐降低。为了仅补偿卫星时钟漂移,RTKLIB利用广播SV时钟参数对SD测量进行校正。最大ʺAge of Differentialʺ设置为处理选项ʺMAX Age of Diffʺ。
在卫星侧SD生成方面,RTKLIB逐次选择一颗具有最大俯仰角的参考卫星。请注意,不同导航系统的卫星之间不会产生卫星侧SD,比如GPS和GLONASS之间。这是因为对于不同导航系统的信号,即使它们具有相同的载波频率,接收机通常也有不同的群时延。接收机的组延迟差称为接收机ISB(系统间偏置)。
假设移动站和基站均采用三频GPS/GNSS接收机,则待估计的未知状态向量x可定义为:
x
=
(
r
r
T
,
v
r
T
,
B
1
T
,
B
2
T
,
B
5
T
)
T
−
−
−
(
E
.
7.7
)
x=(r_r^T,v_r^T,B_1^T,B_2^T,B_5^T)^T---(E.7.7)
x=(rrT,vrT,B1T,B2T,B5T)T−−−(E.7.7)
其中
B
i
=
(
B
r
b
,
i
1
,
B
r
b
,
i
2
,
.
.
.
,
B
r
b
,
i
m
)
B_i=(B_{rb,i}^1,B_{rb,i}^2,...,B_{rb,i}^m)
Bi=(Brb,i1,Brb,i2,...,Brb,im)为LiSD(单差)载波相位偏差(循环)。作为RTKLIB实现,它内部使用SD载波相位偏差而不是DD,以避免参考卫星的麻烦移交处理。SD偏差还有助于解决GLONASS FDMA信号中的整数歧义。
测量矢量y也被DD相位范围和伪距测量定义为:
y
=
(
Φ
1
T
,
Φ
2
T
,
Φ
5
T
,
P
1
T
,
P
2
T
,
P
5
T
)
T
−
−
−
(
E
.
7.8
)
y=(\Phi_1^T,\Phi_2^T,\Phi_5^T,P_1^T,P_2^T,P_5^T)^T---(E.7.8)
y=(Φ1T,Φ2T,Φ5T,P1T,P2T,P5T)T−−−(E.7.8)
其中:
Φ
i
=
(
Φ
r
b
,
i
12
,
Φ
r
b
,
i
13
,
Φ
r
b
,
i
14
,
.
.
.
,
Φ
r
b
,
i
1
m
)
T
P
i
=
(
P
r
b
,
i
12
,
P
r
b
,
i
13
,
P
r
b
,
i
14
,
.
.
.
,
P
r
b
,
i
1
m
)
T
\Phi_i=(\Phi_{rb,i}^{12},\Phi_{rb,i}^{13},\Phi_{rb,i}^{14},...,\Phi_{rb,i}^{1m})^T\\ P_i=(P_{rb,i}^{12},P_{rb,i}^{13},P_{rb,i}^{14},...,P_{rb,i}^{1m})^T
Φi=(Φrb,i12,Φrb,i13,Φrb,i14,...,Φrb,i1m)TPi=(Prb,i12,Prb,i13,Prb,i14,...,Prb,i1m)T
(3)短基线EKF的测量更新
由式(E.7.6),测量模型向量
h
(
x
)
h(x)
h(x),偏导数矩阵
H
(
x
)
H(x)
H(x),测量误差R的协方差矩阵可记为:
h
(
x
)
=
(
h
Φ
,
1
T
,
h
Φ
,
2
T
,
h
Φ
,
5
T
,
h
P
,
1
T
,
h
P
,
2
T
,
h
P
,
5
T
)
T
−
−
−
(
E
.
7.9
)
h(x)=({h_{\Phi,1}}^T,{h_{\Phi,2}}^T,{h_{\Phi,5}}^T,{h_{P,1}}^T,{h_{P,2}}^T,{h_{P,5}}^T)^T---(E.7.9)
h(x)=(hΦ,1T,hΦ,2T,hΦ,5T,hP,1T,hP,2T,hP,5T)T−−−(E.7.9)
H ( x ) = ∂ h ( x ) ∂ x ∣ x = x ^ = ( − D E 0 λ 1 D 0 0 − D E 0 0 λ 2 D 0 − D E 0 0 0 λ 5 D − D E 0 0 0 0 − D E 0 0 0 0 − D E 0 0 0 0 ) − − − ( E . 7.10 ) H(x)=\frac{\partial h(x)}{\partial x}\vert_{x=\hat x}= \begin{pmatrix} -DE& 0&\lambda_1D&0&0 \\ -DE& 0&0&\lambda_2D&0 \\ -DE& 0&0&0&\lambda_5D \\ -DE& 0&0&0&0 \\ -DE& 0&0&0&0 \\ -DE& 0&0&0&0 \\ \end{pmatrix}---(E.7.10) H(x)=∂x∂h(x)∣x=x^=⎝⎜⎜⎜⎜⎜⎜⎛−DE−DE−DE−DE−DE−DE000000λ1D000000λ2D000000λ5D000⎠⎟⎟⎟⎟⎟⎟⎞−−−(E.7.10)
R
=
(
D
R
Φ
,
1
D
T
D
R
Φ
,
2
D
T
D
R
Φ
,
5
D
T
D
R
P
,
1
D
T
D
R
P
,
2
D
T
D
R
P
,
5
D
T
)
−
−
−
(
E
.
7.11
)
R=\begin{pmatrix} DR_{\Phi,1}D^T&&&&& \\ &DR_{\Phi,2}D^T&&&& \\ &&DR_{\Phi,5}D^T&&& \\ &&&DR_{P,1}D^T&&\\ &&&&DR_{P,2}D^T&\\ &&&&&DR_{P,5}D^T\\ \end{pmatrix}---(E.7.11)
R=⎝⎜⎜⎜⎜⎜⎜⎛DRΦ,1DTDRΦ,2DTDRΦ,5DTDRP,1DTDRP,2DTDRP,5DT⎠⎟⎟⎟⎟⎟⎟⎞−−−(E.7.11)
其中:
h
(
x
)
=
(
ρ
r
b
12
+
λ
i
(
B
r
b
1
−
B
r
b
2
)
ρ
r
b
13
+
λ
i
(
B
r
b
1
−
B
r
b
3
)
⋮
ρ
r
b
1
m
+
λ
i
(
B
r
b
1
−
B
r
b
m
)
)
H
=
(
ρ
r
b
12
ρ
r
b
13
⋮
ρ
r
b
1
m
)
h(x)=\begin{pmatrix} \rho_{rb}^{12}+\lambda_i(B_{rb}^1-B_{rb}^2)\\ \rho_{rb}^{13}+\lambda_i(B_{rb}^1-B_{rb}^3)\\ \vdots\\ \rho_{rb}^{1m}+\lambda_i(B_{rb}^1-B_{rb}^m)\\ \end{pmatrix} H=\begin{pmatrix} \rho_{rb}^{12}\\ \rho_{rb}^{13}\\ \vdots\\ \rho_{rb}^{1m}\\ \end{pmatrix}
h(x)=⎝⎜⎜⎜⎛ρrb12+λi(Brb1−Brb2)ρrb13+λi(Brb1−Brb3)⋮ρrb1m+λi(Brb1−Brbm)⎠⎟⎟⎟⎞H=⎝⎜⎜⎜⎛ρrb12ρrb13⋮ρrb1m⎠⎟⎟⎟⎞
(
1
−
1
0
⋯
0
1
0
−
1
⋯
0
⋮
⋮
⋮
⋱
⋮
1
0
0
⋯
−
1
)
:
S
D
(
单
差
)
矩
阵
\begin{pmatrix} 1&-1&0&\cdots&0 \\ 1&0&-1&\cdots&0 \\ \vdots&\vdots&\vdots&\ddots&\vdots \\ 1&0&0&\cdots&-1 \\ \end{pmatrix}:SD(单差)矩阵
⎝⎜⎜⎜⎛11⋮1−10⋮00−1⋮0⋯⋯⋱⋯00⋮−1⎠⎟⎟⎟⎞:SD(单差)矩阵
E
:
(
e
r
1
,
e
r
2
,
.
.
.
,
e
r
m
)
T
R
Φ
,
i
:
d
i
a
g
(
2
σ
Φ
,
i
1
2
,
2
σ
Φ
,
i
2
2
,
.
.
.
,
2
σ
Φ
,
i
m
2
)
R
P
,
i
:
d
i
a
g
(
2
σ
P
,
i
1
2
,
2
σ
P
,
i
2
2
,
.
.
.
,
2
σ
P
,
i
m
2
)
σ
Φ
,
i
s
:
L
i
相
位
范
围
测
量
误
差
的
标
准
偏
差
(
m
)
σ
P
,
i
s
:
L
i
伪
距
测
量
误
差
的
标
准
偏
差
(
m
)
\begin{aligned} E&:(e_r^1,e_r^2,...,e_r^m)^T\\ R_{\Phi,i}&:diag({2\sigma_{\Phi,i}^1}^2,{2\sigma_{\Phi,i}^2}^2,...,{2\sigma_{\Phi,i}^m}^2 )\\ R_{P,i}&:diag({2\sigma_{P,i}^1}^2,{2\sigma_{P,i}^2}^2,...,{2\sigma_{P,i}^m}^2 )\\ \sigma_{\Phi,i}^s&:L_i相位范围测量误差的标准偏差(m)\\ \sigma_{P,i}^s&:L_i伪距测量误差的标准偏差(m)\\ \end{aligned}
ERΦ,iRP,iσΦ,isσP,is:(er1,er2,...,erm)T:diag(2σΦ,i12,2σΦ,i22,...,2σΦ,im2):diag(2σP,i12,2σP,i22,...,2σP,im2):Li相位范围测量误差的标准偏差(m):Li伪距测量误差的标准偏差(m)
利用这些方程求解EKF公式(E.7.1),得到了移动站天线的位置、速度和浮动SD载波相位偏差的历元时间tk。
(4)EKF时间更新
对于RTKLIB中带有接收机动力学的动态定位模式(positioning mode = kinematic and REC dynamics = ON), EKF (E.7.2)的时间更新表示为:
F
k
k
+
1
=
(
I
3
×
3
I
3
×
3
τ
r
I
3
×
3
I
(
3
m
−
3
)
×
(
3
m
−
3
)
)
,
Q
k
k
+
1
=
(
0
3
×
3
Q
v
0
(
3
m
−
3
)
×
(
3
m
−
3
)
)
−
−
−
(
E
.
7.12
)
F_k^{k+1}= \begin{pmatrix} I_{3\times3}& I_{3\times3}\tau_r& \\ & I_{3\times3}& \\ & &I_{(3m-3)\times(3m-3)} \\ \end{pmatrix}, Q_k^{k+1}= \begin{pmatrix} 0_{3\times3}&& \\ & Q_v& \\ & &0_{(3m-3)\times(3m-3)} \\ \end{pmatrix}---(E.7.12)
Fkk+1=⎝⎛I3×3I3×3τrI3×3I(3m−3)×(3m−3)⎠⎞,Qkk+1=⎝⎛03×3Qv0(3m−3)×(3m−3)⎠⎞−−−(E.7.12)
其中:
Q
v
=
E
r
T
d
i
a
g
(
σ
v
e
2
τ
r
,
σ
v
n
2
τ
r
,
σ
v
u
2
τ
r
)
E
r
Q_v={E_r}^Tdiag({\sigma_{ve}}^2\tau_r,{\sigma_{vn}}^2\tau_r,{\sigma_{vu}}^2\tau_r)E_r
Qv=ErTdiag(σve2τr,σvn2τr,σvu2τr)Er
而且
τ
r
=
t
k
−
1
−
t
k
\tau_r=t_{k-1}-t_k
τr=tk−1−tk为GPS/GNSS接收机采样间隔(s),
(
σ
v
e
,
σ
v
n
,
σ
v
u
)
(\sigma_{ve},\sigma_{vn},\sigma_{vu})
(σve,σvn,σvu)为探测器速度系统噪声东、北、上分量的标准差(
m
/
s
/
s
m/s /\sqrt s
m/s/s)。
对于没有接收机动力学的纯动态模式(Positioning mode = kinematic and REC dynamics = OFF),公式(E.7.9)代入:
F
k
k
+
1
=
(
I
3
×
3
I
3
×
3
I
(
3
m
−
3
)
×
(
3
m
−
3
)
)
,
Q
k
k
+
1
=
(
∞
3
×
3
0
3
×
3
0
(
3
m
−
3
)
×
(
3
m
−
3
)
)
−
−
−
(
E
.
7.13
)
F_k^{k+1}= \begin{pmatrix} I_{3\times3}& & \\ & I_{3\times3}& \\ & &I_{(3m-3)\times(3m-3)} \\ \end{pmatrix}, Q_k^{k+1}= \begin{pmatrix} \infty_{3\times3}&& \\ &0_{3\times3}& \\ & &0_{(3m-3)\times(3m-3)} \\ \end{pmatrix}---(E.7.13)
Fkk+1=⎝⎛I3×3I3×3I(3m−3)×(3m−3)⎠⎞,Qkk+1=⎝⎛∞3×303×30(3m−3)×(3m−3)⎠⎞−−−(E.7.13)
为了避免数字不稳定,请在接收器位置的方差中添加无限的处理噪声,接收器位置状态将重置为每个时期的初始猜测值,并且将足够大的过程噪声(
1
0
4
m
2
10^4m^2
104m2)添加到RTKLIB的方差中。 初始位置来自单点定位过程,该过程用于避免非线性信号测量模型的迭代。
在静态模式下(Positioning mode = static),公式(F.7.10)可简单替换为:
F
k
k
+
1
=
(
I
3
×
3
I
3
×
3
I
(
3
m
−
3
)
×
(
3
m
−
3
)
)
,
Q
k
k
+
1
=
(
0
3
×
3
0
3
×
3
0
(
3
m
−
3
)
×
(
3
m
−
3
)
)
−
−
−
(
E
.
7.14
)
F_k^{k+1}= \begin{pmatrix} I_{3\times3}& & \\ & I_{3\times3}& \\ & &I_{(3m-3)\times(3m-3)} \\ \end{pmatrix}, Q_k^{k+1}= \begin{pmatrix} 0_{3\times3}&& \\ &0_{3\times3}& \\ & &0_{(3m-3)\times(3m-3)} \\ \end{pmatrix}---(E.7.14)
Fkk+1=⎝⎛I3×3I3×3I(3m−3)×(3m−3)⎠⎞,Qkk+1=⎝⎛03×303×30(3m−3)×(3m−3)⎠⎞−−−(E.7.14)
在瞬时模糊度分辨模式(整数模糊度分辨=瞬时)中,SD载波相位偏差Bi的时间更新处理方式与上述略有不同。在这种模式下,载波相位偏置状态的值不会被EKF时间更新后续到下一个epoch。在每个时点,将偏差重置为初始的猜测值,并在方差中加入足够大的过程噪声(104 m2)。如果在测量数据中检测到一个周跳,相应的SD载波相位偏置状态也被重置为初始值。RTKLIB通过输入测量数据中的LLI(丢失锁定指示器)检测周期滑移,如果双频测量可用,则通过几何自由LC(线性组合)检测相位跳变。周期滑移阈值可以通过处理选项ʺslip Thresʺ来更改。
(5)整数模糊度求解
一旦EKF测量中得到的估计状态更新,浮动载波相位模糊可以被分解成整数值,以提高精度和收敛时间(整模糊度Res =连续、瞬时或固定保持)。首先,将估计状态及其协方差矩阵转化为DD形式:
x
^
k
′
=
G
x
^
k
(
+
)
=
(
r
^
r
T
,
v
^
r
T
,
N
^
T
)
T
−
−
−
(
E
.
7.15
)
P
k
′
=
G
P
k
(
+
)
G
T
=
(
Q
R
Q
N
R
Q
R
N
Q
N
)
−
−
−
(
E
.
7.16
)
\hat x'_k=G\hat x_k(+)=({\hat r_r}^T,{\hat v_r}^T,{\hat N}^T)^T---(E.7.15)\\ P'_k=GP_k(+)G^T= \begin{pmatrix} Q_R&Q_{NR}\\ Q_{RN}&Q_N\\ \end{pmatrix}---(E.7.16)
x^k′=Gx^k(+)=(r^rT,v^rT,N^T)T−−−(E.7.15)Pk′=GPk(+)GT=(QRQRNQNRQN)−−−(E.7.16)
其中:
G
=
(
I
6
×
6
D
D
D
)
:
S
D
到
D
D
的
变
换
矩
阵
G= \begin{pmatrix} I_{6\times 6}&&&\\ &D&&\\ &&D&\\ &&&D\\ \end{pmatrix}:SD到DD的变换矩阵
G=⎝⎜⎜⎛I6×6DDD⎠⎟⎟⎞:SD到DD的变换矩阵
在这种转换中,SD载波相位偏差被转移到DD载波相位形式,以消除接收机初始相位项,获得整数模糊
N
^
\hat N
N^及其协方差QN。在这些公式中,通过求解整数最小二乘问题得到最适合整数二义性的整数向量
N
˘
\breve N
N˘,表达式为:
N
˘
=
a
r
g
m
i
n
(
(
N
−
N
^
)
T
Q
N
−
1
(
N
−
N
^
)
)
(
N
∈
Z
)
−
−
−
(
E
.
7.17
)
\breve N=argmin((N-\hat N)^TQ_N^{-1}(N-\hat N))\quad (N\in Z)---(E.7.17)
N˘=argmin((N−N^)TQN−1(N−N^))(N∈Z)−−−(E.7.17)
为了解决盲降问题,RTKLIB采用了一种众所周知的高效搜索策略LAMBDA[66]及其扩展MLAMBDA[67]。LAMBDA和MLAMBDA提供了一个线性变换的组合,以缩小整数向量搜索空间和一个转换空间中熟练的树搜索程序。通过以下简单的ʺRatio‐Testʺ来验证这些程序的整数向量解。在ʺRatio‐Testʺ中,ratio‐factor R,定义为第二优解
N
˘
2
\breve N_2
N˘2的残差平方根的加权和与第一优解
N
˘
\breve N
N˘的加权和之比,用于检验解的可靠性。验证阈值Rthres可以通过处理选项ʺ最小修正歧义比率ʺ来设置。当前版本RTKLIB只支持固定的阈值。
R
=
(
N
˘
2
−
N
^
)
T
Q
N
−
1
(
N
˘
2
−
N
^
)
(
N
˘
−
N
^
)
T
Q
N
−
1
(
N
˘
−
N
^
)
>
R
t
h
r
e
s
−
−
−
(
E
.
7.18
)
R=\frac{(\breve N_2-\hat N)^T{Q_N}^{-1}(\breve N_2-\hat N)}{(\breve N-\hat N)^T{Q_N}^{-1}(\breve N-\hat N)}>R_{thres}---(E.7.18)
R=(N˘−N^)TQN−1(N˘−N^)(N˘2−N^)TQN−1(N˘2−N^)>Rthres−−−(E.7.18)
验证后,通过求解下式得到rover天线位置和速度
r
˘
r
\breve r_r
r˘r和
v
˘
r
\breve v_r
v˘r的ʺFIXEDʺ解。如果验证失败,RTKLIB输出ʺFLOATʺ解决方案
r
^
r
\hat r_r
r^r和
v
^
r
\hat v_r
v^r。
(
r
˘
r
v
˘
r
)
=
(
r
^
r
v
^
r
)
−
Q
R
N
Q
N
−
1
(
N
^
−
N
˘
)
−
−
−
(
E
.
7.19
)
\begin{pmatrix} \breve r_r\\ \breve v_r \end{pmatrix}= \begin{pmatrix} \hat r_r\\ \hat v_r \end{pmatrix}-Q_{RN}{Q_N}^{-1}(\hat N-\breve N)---(E.7.19)
(r˘rv˘r)=(r^rv^r)−QRNQN−1(N^−N˘)−−−(E.7.19)
如果处理选项设置为ʺFix and Holdʺ模式(整数模糊度解决方案=固定和保持),并且固定的解决方案经之前的测试正确验证,DD载波相位偏置参数被严格约束为固定的整数值。为此,RTKLIB向EKF输入以下ʺ伪ʺ测量值,并通过(F.7.1)更新EKF。
y
=
N
˘
h
(
x
)
=
G
x
H
(
x
)
=
G
R
=
d
i
a
g
(
σ
c
2
,
σ
c
2
,
σ
c
2
,
.
.
.
)
y= \breve N\\ h(x)=Gx\\ H(x)=G\\ R=diag({\sigma_c}^2,{\sigma_c}^2,{\sigma_c}^2,...)
y=N˘h(x)=GxH(x)=GR=diag(σc2,σc2,σc2,...)
其中:
G
=
(
0
D
0
D
0
D
)
:
S
D
到
D
D
的
变
换
矩
阵
σ
c
:
对
固
定
整
数
二
义
性
的
约
束
(
=
0.001
周
期
)
G= \begin{pmatrix} 0&D&&\\ 0&&D&\\ 0&&&D\\ \end{pmatrix}:SD到DD的变换矩阵\\ \sigma_c:对固定整数二义性的约束(= 0.001周期)
G=⎝⎛000DDD⎠⎞:SD到DD的变换矩阵σc:对固定整数二义性的约束(=0.001周期)
ʺFix and Holdʺ模式在RTKLIB版本中首次引入。2.4.0,以提高固定率,特别是在运动模式下跟踪移动接收机。
(6) Long baseline DD measurement model
对于移动站r与基站b之间的长基线处理,可以形成类似于短基线DD模型的DD测量方程:
Φ
r
b
,
i
j
k
=
ρ
r
b
j
k
−
I
r
b
,
k
j
k
+
T
r
b
j
k
+
λ
i
(
B
r
b
,
i
j
−
r
b
,
i
k
)
+
d
Φ
r
,
i
s
+
ε
Φ
P
r
b
,
i
j
k
=
ρ
r
b
j
k
+
I
r
b
,
i
j
k
+
T
r
b
j
k
+
ε
P
−
−
−
(
E
.
7.24
)
\begin{aligned} \Phi_{rb,i}^{jk}&=\rho_{rb}^{jk}-I_{rb,k}^{jk}+T_{rb}^{jk}+\lambda_i(B_{rb,i}^{j}-_{rb,i}^{k})+d\Phi_{r,i}^{s}+\varepsilon_\Phi\\ P_{rb,i}^{jk}&=\rho_{rb}^{jk}+I_{rb,i}^{jk}+T_{rb}^{jk}+\varepsilon_P---(E.7.24) \end{aligned}
Φrb,ijkPrb,ijk=ρrbjk−Irb,kjk+Trbjk+λi(Brb,ij−rb,ik)+dΦr,is+εΦ=ρrbjk+Irb,ijk+Trbjk+εP−−−(E.7.24)
其中, I r i s I_{ri}^{s} Iris和 T r s T_{r}^{s} Trs分别作为电离层延迟(m)和对流层延迟(m)添加到短基线DD模型中。应使用卫星位置的精确星历表来减轻100公里以上基线的广播星历表误差。在载波相位校正项 d Φ r , i s d\Phi_{r,i}^{s} dΦr,is中,对于基线大于500公里的情况,应考虑固体潮效应。为了消除电离层项,有时会形成无电离层LC(线性组合)。然而,RTKLIB没有使用这种明确的LC,而是使用EKF的双频或三频测量来直接估计电离层项,用于基线处理。
长基线情况下的未知状态向量x也可求解为:
假设移动站和基站均采用三频GPS/GNSS接收机,则待估计的未知状态向量x可定义为:
x
=
(
r
r
T
,
v
r
T
,
Z
r
,
G
N
,
r
,
G
E
,
r
,
Z
b
,
G
N
,
b
,
G
E
,
b
,
I
T
,
B
1
T
,
B
2
T
,
B
5
T
)
T
−
−
−
(
E
.
7.25
)
x=(r_r^T,v_r^T,Z_r,G_{N,r},G_{E,r},Z_b,G_{N,b},G_{E,b},I^T,B_1^T,B_2^T,B_5^T)^T---(E.7.25)
x=(rrT,vrT,Zr,GN,r,GE,r,Zb,GN,b,GE,b,IT,B1T,B2T,B5T)T−−−(E.7.25)
其中Zrand Zb为移动站和基站站点的ZTD(天顶总延迟),
G
N
,
r
,
G
E
,
r
,
G
N
,
b
,
G
E
,
b
G_{N,r},G_{E,r},G_{N,b},G_{E,b}
GN,r,GE,r,GN,b,GE,b是对流层梯度的北部和东部分量。
I
=
(
I
r
b
1
,
I
r
b
2
,
.
.
.
,
I
r
b
m
)
I=(I_{rb}^1,I_{rb}^2,...,I_{rb}^m)
I=(Irb1,Irb2,...,Irbm)也是在频率为L1 (m)的SD垂直电离层延迟。
测量模型向量h(x)和偏导数矩阵H(x)可表示为:
由式(E.7.6),测量模型向量
h
(
x
)
h(x)
h(x),偏导数矩阵
H
(
x
)
H(x)
H(x),测量误差R的协方差矩阵可记为:
h
(
x
)
=
(
h
Φ
,
1
T
,
h
Φ
,
2
T
,
h
Φ
,
5
T
,
h
P
,
1
T
,
h
P
,
2
T
,
h
P
,
5
T
)
T
−
−
−
(
E
.
7.26
)
h(x)=({h_{\Phi,1}}^T,{h_{\Phi,2}}^T,{h_{\Phi,5}}^T,{h_{P,1}}^T,{h_{P,2}}^T,{h_{P,5}}^T)^T---(E.7.26)
h(x)=(hΦ,1T,hΦ,2T,hΦ,5T,hP,1T,hP,2T,hP,5T)T−−−(E.7.26)
h Φ , i = ( ρ r b 12 + T r b 12 − γ k ( m I 1 I r b 1 − m I 2 I r b 2 ) + λ i ( B r b 1 − B r b 2 ) + d Φ r b , i 12 ρ r b 13 + T r b 13 − γ k ( m I 1 I r b 1 − m I 3 I r b 3 ) + λ i ( B r b 1 − B r b 3 ) + d Φ r b , i 13 ⋮ ρ r b 1 m + T r b 1 m − γ k ( m I 1 I r b 1 − m I m I r b m ) + λ i ( B r b 1 − B r b m ) + d Φ r b , i 1 m ) − − − ( E . 7.27 ) h_{\Phi,i}=\begin{pmatrix} \rho_{rb}^{12}+T_{rb}^{12}-\gamma_k(m_I^1I_{rb}^1-m_I^2I_{rb}^2)+\lambda_i(B_{rb}^1-B_{rb}^2)+d\Phi_{rb,i}^{12}\\ \rho_{rb}^{13}+T_{rb}^{13}-\gamma_k(m_I^1I_{rb}^1-m_I^3I_{rb}^3)+\lambda_i(B_{rb}^1-B_{rb}^3)+d\Phi_{rb,i}^{13}\\ \vdots\\ \rho_{rb}^{1m}+T_{rb}^{1m}-\gamma_k(m_I^1I_{rb}^1-m_I^mI_{rb}^m)+\lambda_i(B_{rb}^1-B_{rb}^m)+d\Phi_{rb,i}^{1m}\\ \end{pmatrix}---(E.7.27) hΦ,i=⎝⎜⎜⎜⎛ρrb12+Trb12−γk(mI1Irb1−mI2Irb2)+λi(Brb1−Brb2)+dΦrb,i12ρrb13+Trb13−γk(mI1Irb1−mI3Irb3)+λi(Brb1−Brb3)+dΦrb,i13⋮ρrb1m+Trb1m−γk(mI1Irb1−mImIrbm)+λi(Brb1−Brbm)+dΦrb,i1m⎠⎟⎟⎟⎞−−−(E.7.27)
h P , i = ( ρ r b 12 + T r b 12 + γ k ( m I 1 I r b 1 − m I 2 I r b 2 ) ρ r b 13 + T r b 13 + γ k ( m I 1 I r b 1 − m I 3 I r b 3 ) ⋮ ρ r b 1 m + T r b 1 m + γ k ( m I 1 I r b 1 − m I m I r b m ) ) − − − ( E . 7.27 ) h_{P,i}=\begin{pmatrix} \rho_{rb}^{12}+T_{rb}^{12}+\gamma_k(m_I^1I_{rb}^1-m_I^2I_{rb}^2)\\ \rho_{rb}^{13}+T_{rb}^{13}+\gamma_k(m_I^1I_{rb}^1-m_I^3I_{rb}^3)\\ \vdots\\ \rho_{rb}^{1m}+T_{rb}^{1m}+\gamma_k(m_I^1I_{rb}^1-m_I^mI_{rb}^m)\\ \end{pmatrix}---(E.7.27) hP,i=⎝⎜⎜⎜⎛ρrb12+Trb12+γk(mI1Irb1−mI2Irb2)ρrb13+Trb13+γk(mI1Irb1−mI3Irb3)⋮ρrb1m+Trb1m+γk(mI1Irb1−mImIrbm)⎠⎟⎟⎟⎞−−−(E.7.27)
H
(
x
)
=
(
−
D
E
0
D
M
T
,
r
D
M
T
,
b
−
γ
1
D
M
I
λ
1
D
−
D
E
0
D
M
T
,
r
D
M
T
,
b
−
γ
2
D
M
I
λ
2
D
−
D
E
0
D
M
T
,
r
D
M
T
,
b
−
γ
5
D
M
I
λ
5
D
−
D
E
0
D
M
T
,
r
D
M
T
,
b
−
γ
1
D
M
I
−
D
E
0
D
M
T
,
r
D
M
T
,
b
−
γ
2
D
M
I
−
D
E
0
D
M
T
,
r
D
M
T
,
b
−
γ
5
D
M
I
)
−
−
−
(
E
.
7.28
)
H(x)= \begin{pmatrix} -DE& 0&DM_{T,r}&DM_{T,b}&-\gamma_1DM_I&\lambda_1D&& \\ -DE& 0&DM_{T,r}&DM_{T,b}&-\gamma_2DM_I&&\lambda_2D& \\ -DE& 0&DM_{T,r}&DM_{T,b}&-\gamma_5DM_I&&&\lambda_5D \\ -DE& 0&DM_{T,r}&DM_{T,b}&-\gamma_1DM_I&&& \\ -DE& 0&DM_{T,r}&DM_{T,b}&-\gamma_2DM_I&&& \\ -DE& 0&DM_{T,r}&DM_{T,b}&-\gamma_5DM_I&&& \\ \end{pmatrix}---(E.7.28)
H(x)=⎝⎜⎜⎜⎜⎜⎜⎛−DE−DE−DE−DE−DE−DE000000DMT,rDMT,rDMT,rDMT,rDMT,rDMT,rDMT,bDMT,bDMT,bDMT,bDMT,bDMT,b−γ1DMI−γ2DMI−γ5DMI−γ1DMI−γ2DMI−γ5DMIλ1Dλ2Dλ5D⎠⎟⎟⎟⎟⎟⎟⎞−−−(E.7.28)
其中:
R
=
(
D
R
Φ
,
1
D
T
D
R
Φ
,
2
D
T
D
R
Φ
,
5
D
T
D
R
P
,
1
D
T
D
R
P
,
2
D
T
D
R
P
,
5
D
T
)
−
−
−
(
E
.
7.11
)
R=\begin{pmatrix} DR_{\Phi,1}D^T&&&&& \\ &DR_{\Phi,2}D^T&&&& \\ &&DR_{\Phi,5}D^T&&& \\ &&&DR_{P,1}D^T&&\\ &&&&DR_{P,2}D^T&\\ &&&&&DR_{P,5}D^T\\ \end{pmatrix}---(E.7.11)
R=⎝⎜⎜⎜⎜⎜⎜⎛DRΦ,1DTDRΦ,2DTDRΦ,5DTDRP,1DTDRP,2DTDRP,5DT⎠⎟⎟⎟⎟⎟⎟⎞−−−(E.7.11)
其中:
γ
k
=
λ
k
2
/
λ
1
2
M
T
,
r
=
(
m
W
G
,
r
1
(
E
l
r
1
)
m
W
,
r
1
(
E
l
r
1
)
c
o
t
E
l
r
1
c
o
s
A
z
r
1
m
W
,
r
1
(
E
l
r
1
)
c
o
t
E
l
r
1
s
i
n
A
z
r
1
m
W
G
,
r
2
(
E
l
r
2
)
m
W
,
r
2
(
E
l
r
2
)
c
o
t
E
l
r
2
c
o
s
A
z
r
2
m
W
,
r
2
(
E
l
r
2
)
c
o
t
E
l
r
2
s
i
n
A
z
r
2
⋮
⋮
⋮
m
W
G
,
r
m
(
E
l
r
m
)
m
W
,
r
m
(
E
l
r
m
)
c
o
t
E
l
r
m
c
o
s
A
z
r
m
m
W
,
r
m
(
E
l
r
m
)
c
o
t
E
l
r
m
s
i
n
A
z
r
m
)
M
I
=
(
m
I
1
,
m
I
2
,
.
.
.
,
m
I
m
)
T
\gamma_k=\lambda_k^2/\lambda_1^2\\ M_{T,r}=\begin{pmatrix} m_{WG,r}^1(El_r^1)&m_{W,r}^1(El_r^1)cotEl_r^1cosAz_r^1& m_{W,r}^1(El_r^1)cotEl_r^1sinAz_r^1\\ m_{WG,r}^2(El_r^2)&m_{W,r}^2(El_r^2)cotEl_r^2cosAz_r^2& m_{W,r}^2(El_r^2)cotEl_r^2sinAz_r^2\\ \vdots&\vdots&\vdots \\ m_{WG,r}^m(El_r^m)&m_{W,r}^m(El_r^m)cotEl_r^mcosAz_r^m& m_{W,r}^m(El_r^m)cotEl_r^msinAz_r^m\\ \end{pmatrix}\\ M_I=(m_I^1,m_I^2,...,m_I^m)^T
γk=λk2/λ12MT,r=⎝⎜⎜⎜⎛mWG,r1(Elr1)mWG,r2(Elr2)⋮mWG,rm(Elrm)mW,r1(Elr1)cotElr1cosAzr1mW,r2(Elr2)cotElr2cosAzr2⋮mW,rm(Elrm)cotElrmcosAzrmmW,r1(Elr1)cotElr1sinAzr1mW,r2(Elr2)cotElr2sinAzr2⋮mW,rm(Elrm)cotElrmsinAzrm⎠⎟⎟⎟⎞MI=(mI1,mI2,...,mIm)T
长基线情况下EKF的时间更新表示为:
F
k
k
+
1
=
(
I
3
×
3
I
3
×
3
τ
r
I
3
×
3
I
6
×
6
I
m
×
m
I
(
3
m
−
3
)
×
(
3
m
−
3
)
)
−
−
−
(
E
.
7.30
)
Q
k
k
+
1
=
(
0
3
×
3
Q
v
Q
T
Q
I
0
(
3
m
−
3
)
×
(
3
m
−
3
)
)
−
−
−
(
E
.
7.31
)
F_k^{k+1}= \begin{pmatrix} I_{3\times3}& I_{3\times3}\tau_r&&& \\ & I_{3\times3}&&& \\ && I_{6\times6}&& \\ &&& I_{m\times m}& \\ & &&&I_{(3m-3)\times(3m-3)} \\ \end{pmatrix}---(E.7.30)\\ Q_k^{k+1}= \begin{pmatrix} 0_{3\times3}&&&& \\ & Q_v&&& \\ & &Q_T&& \\ & &&Q_I& \\ & &&&0_{(3m-3)\times(3m-3)} \\ \end{pmatrix}---(E.7.31)
Fkk+1=⎝⎜⎜⎜⎜⎛I3×3I3×3τrI3×3I6×6Im×mI(3m−3)×(3m−3)⎠⎟⎟⎟⎟⎞−−−(E.7.30)Qkk+1=⎝⎜⎜⎜⎜⎛03×3QvQTQI0(3m−3)×(3m−3)⎠⎟⎟⎟⎟⎞−−−(E.7.31)
式中,
Q
T
Q_T
QT和
Q
I
Q_I
QI为电离层项和对流层项的过程噪声协方差矩阵。在方程中,每个卫星的探测器和基地站的ZTD和梯度参数以及SD垂直电离层延迟被简单地建模为随机游走。除了估算电离层和对流层术语外,在2.4.1版本中增加了ʺ局部固定ʺ特性,用于长基线处理。这意味着只有所有二义性的一部分被解析为整数值。除固定外的其他歧义仍然作为浮点值挂起。为了判断模糊度是否固定的简单判据在RTKLIB中实现了利用卫星仰角。如果卫星在高度的阈值以下,卫星的模糊度就不是固定的。只有超过阈值的卫星模糊度才被解决为整数。模糊度分辨率的高程阈值可以设置为处理选项ʺ固定支架最小高程ʺ和ʺ固定支架最小高程ʺ,以控制ʺFix and Holdʺ特征。
(7) 移动基站模型
如果移动站和基站接收器都在移动,并且需要漫游者相对于基站的唯一相对位置,则通常使用移动基线模式。移动基线模式可以通过在一个移动平台上安装两个天线来确定精确的姿态。在RTKLIB中,如果处理选项ʺ定位模式ʺ设置为ʺmoving‐Baseʺ,则应用移动基线模式。
在移动基线模式下,基站位置不是固定的,而是在逐时的单点定位过程中估计出来的。一旦获得了基站位置,则通过(1)‐(5)所述的短基线运动学模式来估算固定在估计位置和漫球车位置的基站位置。在这种情况下,只有相对位置是有意义的,这意味着漫游者和基站的绝对位置解只有与点指向模式解相同的精度。
除了移动基线模式的简单实现外,RTKLIB还可以校正漫游者和基站之间的时间差。漫游者接收机和基站接收机不同步。接收机的时钟差通常最大达到2毫秒。在移动非常快的平台中,时钟不同步会导致精度下降。为了校正时钟差,在基线处理之前,对基站位置
r
b
r_b
rb进行校正:
r
b
(
t
r
)
=
r
b
(
t
b
)
+
v
b
(
t
b
)
(
t
r
−
t
b
)
−
−
−
(
E
.
7.32
)
r_b(t_r)=r_b(t_b)+v_b(t_b)(t_r-t_b)---(E.7.32)
rb(tr)=rb(tb)+vb(tb)(tr−tb)−−−(E.7.32)
其中tr和tb分别为通过单点定位过程估计的漫游者和基站的信号接收时间。
v
b
(
t
b
)
v_b(t_b)
vb(tb)也是用多普勒测量估计的基站速度。对于由移动基线模式确定的姿态,如果处理选项ʺ基线长度约束ʺ启用,则可以应用基线长度约束。该约束在EKF测量更新中应用了以下伪测量。
y
=
(
r
b
a
s
e
l
i
n
e
)
−
−
−
(
E
.
7.33
)
h
(
x
)
=
(
∣
r
r
(
t
r
)
−
r
b
(
t
r
)
∣
)
−
−
−
(
E
.
7.34
)
H
(
x
)
=
(
(
r
r
(
t
r
)
−
r
b
(
t
r
)
)
T
r
r
(
t
r
)
−
r
b
(
t
r
)
)
−
−
−
(
E
.
7.35
)
R
=
(
σ
r
2
)
−
−
−
(
E
.
7.36
)
\begin{aligned} y&=(r_{baseline})&---(E.7.33)\\ h(x)&=(\vert r_r(t_r)-r_b(t_r)\vert)&---(E.7.34)\\ H(x)&=(\frac{(r_r(t_r)-r_b(t_r))^T}{r_r(t_r)-r_b(t_r)})&---(E.7.35)\\ R&=(\sigma_r^2)&---(E.7.36) \end{aligned}
yh(x)H(x)R=(rbaseline)=(∣rr(tr)−rb(tr)∣)=(rr(tr)−rb(tr)(rr(tr)−rb(tr))T)=(σr2)−−−(E.7.33)−−−(E.7.34)−−−(E.7.35)−−−(E.7.36)
在基线rbaseline是给定的提前确定基线长度(m)和r的约束是基线长度(米)。应对非量线性在短基线长度的情况下,卡尔曼滤波器的迭代测量更新支持通过设置处理选项ʺʺ滤波器迭代的数量多吗