3.1伪距观测值与伪距单点定位原理
总共有三个时钟,第一个是卫星时钟,第二个接收机时钟,第三个是用户时钟。卫星和接收机在各自时钟的控制下,同步产生相同的测距码。
假设卫星时钟快于真实时钟,则钟差为
δ
t
p
\delta t^p
δtp,假设接收机时钟也快于真实时间,钟差为
δ
t
k
\delta t_k
δtk
中间的那一条是接收机接收到的码元,和自身码元作对比(利用平移求自相关函数)得到码元差
M
k
p
M_k^p
Mkp,码元差转换成传播时间
τ
k
p
\tau_k^p
τkp
测得的传输时间与钟差和理论时间有如下关系:
τ
k
p
=
t
k
p
+
δ
t
k
−
δ
t
p
\tau_k^p=t_k^p+\delta t_k-\delta t^p
τkp=tkp+δtk−δtp
左右同乘光速可得到距离等式
ρ
^
k
p
=
c
⋅
t
k
p
+
c
⋅
δ
t
k
−
c
⋅
δ
t
p
\hat{\rho}_k^p=c\cdot t_k^p+c\cdot \delta t_k-c\cdot \delta t^p
ρ^kp=c⋅tkp+c⋅δtk−c⋅δtp
对于传输的误差主要可以分为对流层误差和电离层的误差,所以距离公式可以写为:
ρ
^
k
p
=
ρ
k
p
+
δ
I
k
p
+
δ
T
k
p
+
c
⋅
δ
t
k
−
c
⋅
δ
t
p
\hat{\rho}_k^p=\rho _k^p+\delta I_k^p+\delta T^p_k+c\cdot \delta t_k-c\cdot \delta t^p
ρ^kp=ρkp+δIkp+δTkp+c⋅δtk−c⋅δtp
对星地之间的几何距离,在测站的近似坐标
(
x
k
0
,
y
k
0
,
z
k
0
)
(x_{k0},y_{k0},z_{k0})
(xk0,yk0,zk0)处用泰勒级数展开并且去掉高阶项可得:
ρ
^
k
p
=
ρ
k
0
p
+
l
k
0
p
δ
x
k
+
m
k
0
p
δ
y
k
+
n
k
0
p
δ
z
k
+
δ
I
k
p
+
δ
T
k
p
+
c
⋅
δ
t
k
−
c
⋅
δ
t
p
\hat{\rho}_k^p=\rho _{k0}^p+l_{k0}^p\delta x_k+m_{k0}^p\delta y_k+n_{k0}^p\delta z_k+\delta I_k^p+\delta T^p_k+c\cdot \delta t_k-c\cdot \delta t^p
ρ^kp=ρk0p+lk0pδxk+mk0pδyk+nk0pδzk+δIkp+δTkp+c⋅δtk−c⋅δtp
考虑到伪距观测量由观测值和随机误差两部分组成:可得伪距观测方程:
ρ
^
k
p
+
ε
k
p
=
ρ
k
0
p
+
l
k
0
p
δ
x
k
+
m
k
0
p
δ
y
k
+
n
k
0
p
δ
z
k
+
δ
I
k
p
+
δ
T
k
p
+
c
⋅
δ
t
k
−
c
⋅
δ
t
p
\hat{\rho}_k^p +\varepsilon_k^p =\rho _{k0}^p+l_{k0}^p\delta x_k+m_{k0}^p\delta y_k+n_{k0}^p\delta z_k+\delta I_k^p+\delta T^p_k+c\cdot \delta t_k-c\cdot \delta t^p
ρ^kp+εkp=ρk0p+lk0pδxk+mk0pδyk+nk0pδzk+δIkp+δTkp+c⋅δtk−c⋅δtp
将电离层误差,对流层误差,卫星钟差合并成常数项,则可以得到一个仅含四个待求未知数的的伪距误差方程:
v
k
p
=
l
k
0
p
δ
x
k
+
m
k
0
p
δ
y
k
+
n
k
0
p
δ
z
k
+
c
⋅
δ
t
k
−
L
k
p
L
k
p
=
ρ
^
k
p
−
ρ
k
0
p
−
δ
I
k
p
−
δ
T
k
p
+
c
⋅
δ
t
p
v_k^p =l_{k0}^p\delta x_k+m_{k0}^p\delta y_k+n_{k0}^p\delta z_k+c\cdot \delta t_k-L^p_k\\ L^p_k=\hat{\rho}_k^p-\rho _{k0}^p-\delta I_k^p-\delta T^p_k+c\cdot \delta t^p
vkp=lk0pδxk+mk0pδyk+nk0pδzk+c⋅δtk−LkpLkp=ρ^kp−ρk0p−δIkp−δTkp+c⋅δtp
之后可以把每一个观测到的卫星列出来,进行计算。
观测到四颗卫星以上时,可以利用最小二乘法,通过解算可以得到各个测站的估值,再加上测站的近似值,最后可以得到测站的坐标。
为了减小前面线性化过程中,舍去高阶项误差的影响,需在最新的测站坐标处,重新进行泰勒级数展开,迭代计算未知参数,直至坐标值变化小于一定的阈值。
在实际工作中,通常习惯将某些估值所对应的协因数矩阵对角线相应元素加和的开方,作为评价测量精度的指标,这个指标称作精度因子,也称做DOP值。根据评价对象的不同,有不同的DOP值,如三维位置精度因子PDOP 、时钟精度因子TDOP 、几何精度因子GDOP 等。
3.2载波观测值和载波观测方程
伪距测量精度一般为码元宽度的1/100,对于P为0.29m,对于C/A码为2.9m。而对于载波来说,波长大约在20个厘米左右。目前,对载波相位的测量精度,换算成距离,大约为1/100个波长,也即载波测量精度可达约2个毫米。
载波观测值
与测距码类似,卫星和接收机在各自时钟的控制下,同步产生相同的载波。
由于卫星时钟存在误差,使得载波的产生存在时间偏差。假设卫星时钟快于真实时间,则卫星端载波相位A产生的时间将会提前。同理,假设接收机时钟快于真实时间,则接收机端载波相位A产生的时间也会提前。由于两个时钟误差不同,相位A在卫星和接收机上产生的时间也不相同,存在着时间差。
然而,两相位相减,只能得到不到一周的小数部分,中间有多少个整周数并不知道。
这时, 通过相位比对,只能得到不到一整周的小数部分F(t1),整周数是个未知数N,这个未知整周数,我们称作模糊度。我们从第一个历元开始,不断跟踪卫星,由于模糊度在各个历元都一样,从而可以求出模糊度。
卫星载波相位可以通过频率和时延来相乘求得:
Φ
k
p
=
f
(
t
k
p
+
δ
t
k
−
δ
t
p
)
\Phi_k^p=f(t_k^p+\delta t_k-\delta t^p)
Φkp=f(tkp+δtk−δtp)
所以载波观测值和载波真实传播时间的关系为:
φ
k
p
=
f
t
k
p
+
f
δ
t
k
−
f
δ
t
p
−
N
k
p
\varphi_k^p=ft_k^p+f\delta t_k-f\delta t^p-N_k^p
φkp=ftkp+fδtk−fδtp−Nkp
带入星地几何距离,并且用C替代载波速度,之后加入电离层对流层误差,最终得到载波观测:
φ
k
p
+
ε
k
p
=
f
C
(
ρ
k
p
−
δ
I
k
p
+
δ
T
k
p
)
+
f
δ
t
k
−
f
δ
t
p
−
N
k
p
\varphi_k^p + \varepsilon_k^p =\frac fC(\rho_k^p-\delta I_k^p+\delta T^p_k)+f\delta t_k-f\delta t^p-N_k^p
φkp+εkp=Cf(ρkp−δIkp+δTkp)+fδtk−fδtp−Nkp
3.3载波相对定位
载波相对定位的核心思想就是利用两站对各卫星的观测值中,含有相近或相同的系统误差,通过对观测方程相减,削弱或消去这些系统误差项,从而实现高精度相对定位。
载波观测方程的线性组合
假设安置在k站和s站的两GNSS接收机,分别对卫星p和卫星q,于t1历元和t2历元进行同步观测,则可得到8个独立的载波相位观测值。
φ
k
p
(
t
1
)
,
φ
k
q
(
t
1
)
,
φ
s
p
(
t
1
)
,
φ
s
p
(
t
1
)
,
φ
k
p
(
t
2
)
,
φ
k
q
(
t
2
)
,
φ
s
p
(
t
2
)
,
φ
s
p
(
t
2
)
,
\varphi_k^p (t_1),\varphi_k^q (t_1),\varphi_s^p (t_1),\varphi_s^p (t_1),\\ \varphi_k^p (t_2),\varphi_k^q (t_2),\varphi_s^p (t_2),\varphi_s^p (t_2),\\
φkp(t1),φkq(t1),φsp(t1),φsp(t1),φkp(t2),φkq(t2),φsp(t2),φsp(t2),
对于t1历元的四个载波观测值来说,可相应列出四个载波观测方程,如式(1)至式(4)所示。下面,我们要对这四个方程进行线性组合。
φ
k
p
+
ε
k
p
=
f
c
(
ρ
k
p
−
δ
I
k
p
+
δ
T
k
p
+
c
δ
t
k
−
c
δ
t
p
)
−
N
k
p
(
1
)
φ
s
p
+
ε
s
p
=
f
c
(
ρ
s
p
−
δ
I
s
p
+
δ
T
s
p
+
c
δ
t
s
−
c
δ
t
p
)
−
N
s
p
(
2
)
φ
k
q
+
ε
k
q
=
f
c
(
ρ
k
q
−
δ
I
k
q
+
δ
T
k
q
+
c
δ
t
k
−
c
δ
t
q
)
−
N
k
q
(
3
)
φ
s
q
+
ε
s
q
=
f
c
(
ρ
s
q
−
δ
I
s
q
+
δ
T
s
q
+
c
δ
t
s
−
c
δ
t
q
)
−
N
s
q
(
4
)
\varphi_k^p + \varepsilon_k^p = \frac fc(\rho_k^p-\delta I_k^p+\delta T^p_k+c\delta t_k-c\delta t^p)-N_k^p (1)\\ \varphi_s^p + \varepsilon_s^p = \frac fc(\rho_s^p-\delta I_s^p+\delta T^p_s+c\delta t_s-c\delta t^p)-N_s^p (2)\\ \varphi_k^q + \varepsilon_k^q = \frac fc(\rho_k^q-\delta I_k^q+\delta T^q_k+c\delta t_k-c\delta t^q)-N_k^q (3)\\ \varphi_s^q + \varepsilon_s^q = \frac fc(\rho_s^q-\delta I_s^q+\delta T^q_s+c\delta t_s-c\delta t^q)-N_s^q (4)\\
φkp+εkp=cf(ρkp−δIkp+δTkp+cδtk−cδtp)−Nkp(1)φsp+εsp=cf(ρsp−δIsp+δTsp+cδts−cδtp)−Nsp(2)φkq+εkq=cf(ρkq−δIkq+δTkq+cδtk−cδtq)−Nkq(3)φsq+εsq=cf(ρsq−δIsq+δTsq+cδts−cδtq)−Nsq(4)
单差
单差,也就是两个观测方程之间求差。在具体做法上,可对同一测站进行星间求差,也可以对同一卫星进行测站间求差。
式(1)和式(2)是两测站对于卫星p的载波观测方程。两式相减可得式(5)
φ
k
p
+
ε
k
p
=
f
c
(
ρ
k
p
−
δ
I
k
p
+
δ
T
k
p
+
c
δ
t
k
−
c
δ
t
p
)
−
N
k
p
(
1
)
φ
s
p
+
ε
s
p
=
f
c
(
ρ
s
p
−
δ
I
s
p
+
δ
T
s
p
+
c
δ
t
s
−
c
δ
t
p
)
−
N
s
p
(
2
)
(
φ
k
p
+
ε
k
p
)
−
(
φ
s
p
+
ε
s
p
)
=
f
c
(
ρ
k
p
−
δ
I
k
p
+
δ
T
k
p
+
c
δ
t
k
−
c
δ
t
p
)
−
N
k
p
−
[
f
c
(
ρ
s
p
−
δ
I
s
p
+
δ
T
s
p
+
c
δ
t
s
−
c
δ
t
p
)
−
N
s
p
]
(
5
)
\varphi_k^p + \varepsilon_k^p = \frac fc(\rho_k^p-\delta I_k^p+\delta T^p_k+c\delta t_k-c\delta t^p)-N_k^p (1)\\ \varphi_s^p + \varepsilon_s^p = \frac fc(\rho_s^p-\delta I_s^p+\delta T^p_s+c\delta t_s-c\delta t^p)-N_s^p (2)\\ (\varphi_k^p + \varepsilon_k^p) -(\varphi_s^p + \varepsilon_s^p) = \frac fc(\rho_k^p-\delta I_k^p+\delta T^p_k+c\delta t_k-c\delta t^p)-N_k^p \\-[\frac fc(\rho_s^p-\delta I_s^p+\delta T^p_s+c\delta t_s-c\delta t^p)-N_s^p](5)\\
φkp+εkp=cf(ρkp−δIkp+δTkp+cδtk−cδtp)−Nkp(1)φsp+εsp=cf(ρsp−δIsp+δTsp+cδts−cδtp)−Nsp(2)(φkp+εkp)−(φsp+εsp)=cf(ρkp−δIkp+δTkp+cδtk−cδtp)−Nkp−[cf(ρsp−δIsp+δTsp+cδts−cδtp)−Nsp](5)
电离层误差、对流层误差认为近似相等,卫星时钟误差相同,于是可以
消去这些误差项,进而得到式(6)。式(6)即为t1历元,卫星p的站间单差载波观测方程,同理可以通过(3)和(4)得到卫星q的站间单差(7):
Δ
φ
k
s
p
+
Δ
ε
k
s
p
=
f
c
(
ρ
k
p
−
ρ
s
p
+
c
δ
t
k
s
)
−
Δ
N
k
s
p
(
6
)
Δ
φ
k
s
q
+
Δ
ε
k
s
q
=
f
c
(
ρ
k
q
−
ρ
s
q
+
c
δ
t
k
s
)
−
Δ
N
k
s
q
(
7
)
\Delta\varphi_{ks}^p + \Delta\varepsilon_{ks}^p = \frac fc(\rho_k^p-\rho_s^p+c\delta t_{ks})-\Delta N_{ks}^p (6)\\ \Delta\varphi_{ks}^q + \Delta\varepsilon_{ks}^q = \frac fc(\rho_k^q-\rho_s^q+c\delta t_{ks})-\Delta N_{ks}^q (7)\\
Δφksp+Δεksp=cf(ρkp−ρsp+cδtks)−ΔNksp(6)Δφksq+Δεksq=cf(ρkq−ρsq+cδtks)−ΔNksq(7)
双差
在单差方程,式(6)和式(7)中,含有相同的接收机钟差残差项。于是可对单差方程,式(6)和式(7)进行求差,进而消去接收机钟差残差项,得式(8)。这个过程称作双差,式(8)称作双差载波观测方程。
∇
Δ
φ
k
s
p
q
+
∇
Δ
ε
k
s
p
q
=
f
c
(
ρ
k
p
−
ρ
s
p
−
ρ
k
q
+
ρ
s
q
)
−
Δ
∇
N
k
s
p
q
(
8
)
\nabla\Delta\varphi_{ks}^{pq} +\nabla \Delta\varepsilon_{ks}^{pq} = \frac fc(\rho_k^p-\rho_s^p-\rho_k^q+\rho_s^q)-\Delta\nabla N_{ks}^{pq} (8)\\
∇Δφkspq+∇Δεkspq=cf(ρkp−ρsp−ρkq+ρsq)−Δ∇Nkspq(8)
三差
式(8)是t1历元的双差载波观测方程。 同理,可得t2历元的双差载波观测方程,如式(9)所示。
∇
Δ
φ
k
s
p
q
(
t
1
)
+
∇
Δ
ε
k
s
p
q
(
t
1
)
=
f
c
(
ρ
k
p
(
t
1
)
−
ρ
s
p
(
t
1
)
−
ρ
k
q
(
t
1
)
+
ρ
s
q
(
t
1
)
)
−
∇
N
k
s
p
q
(
8
)
∇
Δ
φ
k
s
p
q
(
t
2
)
+
∇
Δ
ε
k
s
p
q
(
t
2
)
=
f
c
(
ρ
k
p
(
t
2
)
−
ρ
s
p
(
t
2
)
−
ρ
k
q
(
t
2
)
+
ρ
s
q
(
t
2
)
)
−
∇
N
k
s
p
q
(
9
)
\nabla\Delta\varphi_{ks}^{pq}(t_1) +\nabla \Delta\varepsilon_{ks}^{pq} (t_1)= \frac fc(\rho_k^p(t_1)-\rho_s^p(t_1)-\rho_k^q(t_1)+\rho_s^q(t_1))-\nabla N_{ks}^{pq}(8)\\ \nabla\Delta\varphi_{ks}^{pq}(t_2) +\nabla \Delta\varepsilon_{ks}^{pq} (t_2)= \frac fc(\rho_k^p(t_2)-\rho_s^p(t_2)-\rho_k^q(t_2)+\rho_s^q(t_2))-\nabla N_{ks}^{pq}(9)\\
∇Δφkspq(t1)+∇Δεkspq(t1)=cf(ρkp(t1)−ρsp(t1)−ρkq(t1)+ρsq(t1))−∇Nkspq(8)∇Δφkspq(t2)+∇Δεkspq(t2)=cf(ρkp(t2)−ρsp(t2)−ρkq(t2)+ρsq(t2))−∇Nkspq(9)
和式(9)中,含有相同的双差模糊度项,两式相减,也即历元间求差,消去双差模糊度项,可得式(10),这个过程称作三差。
δ
∇
Δ
φ
k
s
p
q
(
t
12
)
+
δ
∇
Δ
ε
k
s
p
q
(
t
12
)
=
f
c
(
ρ
k
p
(
t
1
)
−
ρ
s
p
(
t
1
)
−
ρ
k
q
(
t
1
)
+
ρ
s
q
(
t
1
)
)
−
f
c
(
ρ
k
p
(
t
2
)
−
ρ
s
p
(
t
2
)
−
ρ
k
q
(
t
2
)
+
ρ
s
q
(
t
2
)
)
(
10
)
\delta\nabla\Delta\varphi_{ks}^{pq}(t_{12}) +\delta\nabla \Delta\varepsilon_{ks}^{pq} (t_{12})= \frac fc(\rho_k^p(t_1)-\rho_s^p(t_1)-\rho_k^q(t_1)+\rho_s^q(t_1))-\\ \frac fc(\rho_k^p(t_2)-\rho_s^p(t_2)-\rho_k^q(t_2)+\rho_s^q(t_2))(10)\\
δ∇Δφkspq(t12)+δ∇Δεkspq(t12)=cf(ρkp(t1)−ρsp(t1)−ρkq(t1)+ρsq(t1))−cf(ρkp(t2)−ρsp(t2)−ρkq(t2)+ρsq(t2))(10)
小结
从表中可见,双差已消除了除模糊度外的各系统误差项。三差消除了包括模糊度等各系统误差项。由于模糊度存在整数特性,这个特性对提高定位精度具有非常重要的作用,应保留在观测方程中。因此,实际数据处理中,普遍采用双差的形式。
3.4基线的运算
01基线方程推导
我们先以测站k和测站s与卫星P的几何关系为例,建立卫星P至两测站的距离与基线向量的数学关系。
如图,设测站k至卫星P的单位矢量为
u
k
p
u_k^p
ukp,测站s至卫星p的单位矢量为
u
s
p
u_s^p
usp。并得到两单位矢量的差,如式(1)所示;以及两单位矢量的平均值,如式(2)所示。
Δ
u
p
=
u
k
p
−
u
s
p
(
1
)
u
m
p
=
1
2
(
u
k
p
+
u
s
p
)
(
2
)
\Delta u^p=u_k^p-u_s^p(1)\\ u_m^p=\frac12(u_k^p+u_s^p)(2)
Δup=ukp−usp(1)ump=21(ukp+usp)(2)
然后在图上作辅助线形成等腰三角形,辅助线向量为
b
1
p
b_1^p
b1p,可表示为:
b
1
p
=
ρ
s
p
u
s
p
−
ρ
s
p
u
k
p
=
ρ
s
p
(
u
s
p
−
u
k
p
)
=
−
ρ
s
p
Δ
u
p
(
3
)
b_1^p=\rho_s^pu_s^p-\rho _s^pu_k^p\\ =\rho _s^p(u_s^p-u_k^p)=-\rho _s^p\Delta u^p(3)
b1p=ρspusp−ρspukp=ρsp(usp−ukp)=−ρspΔup(3)
两地的星地距离差可用(4)表示,向量可用(5)表示:
Δ
ρ
p
=
ρ
k
p
−
ρ
s
p
(
4
)
Δ
ρ
p
u
k
p
=
b
+
b
1
p
(
5
)
\Delta \rho^p=\rho_k^p-\rho_s^p(4)\\ \Delta \rho^pu_k^p=b+b_1^p(5)
Δρp=ρkp−ρsp(4)Δρpukp=b+b1p(5)
将式(3)代入式(5),星地距离差的矢量可以进一步表示成式(6)形式。顾及到式(2),对式(6)两边点乘两单位矢量的平均值,可得式(7)。如式(8)所示,两单位矢量的平均值与其差值点乘的结果等于零。所以式(7)可以简化成式(9)的形式。
Δ
ρ
p
u
k
p
=
b
−
ρ
s
p
Δ
u
p
(
6
)
Δ
ρ
p
u
m
p
u
k
p
=
u
m
p
b
−
ρ
s
p
(
u
m
p
Δ
u
p
)
(
7
)
\Delta \rho^pu_k^p=b-\rho _s^p\Delta u^p(6)\\ \Delta \rho^pu_m^pu_k^p=u_m^pb-\rho _s^p(u_m^p\Delta u^p)(7)
Δρpukp=b−ρspΔup(6)Δρpumpukp=umpb−ρsp(umpΔup)(7)
u
m
p
Δ
u
p
=
1
2
(
u
k
p
+
u
s
p
)
(
u
k
p
−
u
s
p
)
=
1
2
(
u
k
p
u
k
p
+
u
s
p
u
k
p
−
u
k
p
u
s
p
−
u
s
p
u
s
p
)
(
8
)
=
1
2
(
1
−
1
)
=
0
u_m^p\Delta u^p=\frac12(u_k^p+u_s^p)(u_k^p-u_s^p)\\= \frac12(u_k^pu_k^p+u_s^pu_k^p-u_k^pu_s^p-u_s^pu_s^p)(8)\\= \frac12(1-1)=0\\
umpΔup=21(ukp+usp)(ukp−usp)=21(ukpukp+uspukp−ukpusp−uspusp)(8)=21(1−1)=0
Δ
ρ
p
u
m
p
u
k
p
=
u
m
p
b
(
9
)
\Delta \rho^pu_m^pu_k^p=u_m^pb(9)
Δρpumpukp=umpb(9)
将式(2)代入式(9), 可得式(10)。
u
m
p
b
=
Δ
ρ
p
u
m
p
u
k
p
=
Δ
ρ
p
1
2
(
u
k
p
+
u
s
p
)
u
k
p
=
Δ
ρ
p
1
2
(
1
+
u
s
p
u
k
p
)
(
10
)
u_m^pb=\Delta\rho^pu_m^pu_k^p=\Delta\rho^p\frac12(u_k^p+u_s^p)u_k^p\\=\Delta\rho^p\frac12(1+u_s^pu_k^p)(10)
umpb=Δρpumpukp=Δρp21(ukp+usp)ukp=Δρp21(1+uspukp)(10)
设星端顶角为欧米伽p,则两星地方向单位矢量,点乘,可以表示成,星端顶角欧米伽p的函数,如式(11)所示.
u
s
p
u
k
p
=
∣
u
s
p
∣
∣
u
k
p
∣
c
o
s
ω
p
=
c
o
s
ω
p
(
11
)
u_s^pu_k^p=\vert u_s^p\vert\vert u_k^p\vert cos\omega^p=cos\omega^p(11)
uspukp=∣usp∣∣ukp∣cosωp=cosωp(11)
将式(11)代入式(10),可得式(12)。对式(12)进行整理得到式(13),从而实现,将卫星p与两测站的星地距离差表示成基线向量的函数。同理,对于卫星q,也可以得到类似的函数形式,如式(14)所示。
u
m
p
b
=
Δ
ρ
p
1
2
(
1
+
c
o
s
ω
p
)
=
Δ
ρ
p
⋅
c
o
s
2
w
p
2
(
12
)
Δ
ρ
p
=
s
e
c
2
w
p
2
⋅
u
m
p
b
(
13
)
Δ
ρ
q
=
s
e
c
2
w
q
2
⋅
u
m
q
b
(
14
)
u_m^pb=\Delta\rho^p\frac12(1+cos\omega^p)=\Delta\rho^p\cdot cos^2\frac{w^p}2(12)\\ \Delta\rho^p=sec^2\frac{w^p}2\cdot u_m^pb(13)\\ \Delta\rho^q=sec^2\frac{w^q}2\cdot u_m^qb(14)\\
umpb=Δρp21(1+cosωp)=Δρp⋅cos22wp(12)Δρp=sec22wp⋅umpb(13)Δρq=sec22wq⋅umqb(14)
现在回到双差载波观测方程,为便于讨论,我们将双差载波观测方程重新写在这里,如式(15)所示。我们将式(13)、式(14)带入式(15),可得式(16)。当基线长度小于40km时,可将星端顶角欧米伽p近似看做零,也即式(17)近似成立。
∇
Δ
φ
k
s
p
q
+
∇
Δ
ε
k
s
p
q
=
f
c
(
ρ
k
p
−
ρ
s
p
−
ρ
k
q
+
ρ
s
q
)
−
∇
Δ
N
k
s
p
q
=
f
c
(
Δ
ρ
p
−
Δ
ρ
q
)
−
∇
Δ
N
k
s
p
q
(
15
)
∇
Δ
φ
k
s
p
q
+
∇
Δ
ε
k
s
p
q
=
f
c
(
s
e
c
2
w
p
2
⋅
u
m
p
b
−
s
e
c
2
w
q
2
⋅
u
m
q
b
)
−
∇
Δ
N
k
s
p
q
(
16
)
s
e
c
2
w
p
2
≈
s
e
c
2
w
q
2
≈
1
(
17
)
\nabla\Delta\varphi_{ks}^{pq} +\nabla \Delta\varepsilon_{ks}^{pq} = \frac fc(\rho_k^p-\rho_s^p-\rho_k^q+\rho_s^q)-\nabla\Delta N_{ks}^{pq} \\=\frac fc(\Delta\rho^p-\Delta\rho^q)-\nabla\Delta N_{ks}^{pq} (15)\\ \nabla\Delta\varphi_{ks}^{pq} +\nabla \Delta\varepsilon_{ks}^{pq}=\frac fc(sec^2\frac{w^p}2\cdot u_m^pb-sec^2\frac{w^q}2\cdot u_m^qb)-\nabla\Delta N_{ks}^{pq} (16)\\ sec^2\frac{w^p}2\approx sec^2\frac{w^q}2\approx 1(17)
∇Δφkspq+∇Δεkspq=cf(ρkp−ρsp−ρkq+ρsq)−∇ΔNkspq=cf(Δρp−Δρq)−∇ΔNkspq(15)∇Δφkspq+∇Δεkspq=cf(sec22wp⋅umpb−sec22wq⋅umqb)−∇ΔNkspq(16)sec22wp≈sec22wq≈1(17)
将式(17)代入式(16),并顾及到式(2),可得式(18)。
∇
Δ
φ
k
s
p
q
+
∇
Δ
ε
k
s
p
q
=
f
c
[
u
m
p
−
u
m
q
]
⋅
b
−
∇
Δ
N
k
s
p
q
=
1
2
f
c
[
(
u
k
p
+
u
s
p
)
−
(
u
k
q
+
u
s
q
)
]
⋅
b
−
∇
Δ
N
k
s
p
q
(
18
)
\begin{aligned} &\nabla\Delta\varphi_{ks}^{pq} +\nabla \Delta\varepsilon_{ks}^{pq}\\ &=\frac fc[u_m^p-u_m^q]\cdot b-\nabla\Delta N_{ks}^{pq} \\ &=\frac 12\frac fc[(u_k^p+u_s^p)-(u_k^q+u_s^q)]\cdot b-\nabla\Delta N_{ks}^{pq} \end{aligned}(18)
∇Δφkspq+∇Δεkspq=cf[ump−umq]⋅b−∇ΔNkspq=21cf[(ukp+usp)−(ukq+usq)]⋅b−∇ΔNkspq(18)
在式(18)中,含有各星地距离单位矢量。这些单位矢量可以根据式(19),用卫星坐标和测站近似坐标,计算得到:
u
k
p
=
1
ρ
k
p
(
Δ
X
k
p
Δ
Y
k
p
Δ
Z
k
p
)
u
k
q
=
1
ρ
k
q
(
Δ
X
k
q
Δ
Y
k
q
Δ
Z
k
q
)
u
s
p
=
1
ρ
s
p
(
Δ
X
s
p
Δ
Y
s
p
Δ
Z
s
p
)
u
s
q
=
1
ρ
s
q
(
Δ
X
s
q
Δ
Y
s
q
Δ
Z
s
q
)
(
19
)
\begin{aligned} u_k^p=\frac1{\rho_k^p} \begin{pmatrix} \Delta X_k^p \\ \Delta Y_k^p \\ \Delta Z_k^p \end{pmatrix}\\ \end{aligned} \begin{aligned} u_k^q=\frac1{\rho_k^q} \begin{pmatrix} \Delta X_k^q \\ \Delta Y_k^q \\ \Delta Z_k^q \end{pmatrix}\\ \end{aligned}\\ \begin{aligned} u_s^p=\frac1{\rho_s^p} \begin{pmatrix} \Delta X_s^p \\ \Delta Y_s^p \\ \Delta Z_s^p \end{pmatrix}\\ \end{aligned} \begin{aligned} u_s^q=\frac1{\rho_s^q} \begin{pmatrix} \Delta X_s^q \\ \Delta Y_s^q \\ \Delta Z_s^q \end{pmatrix}\\ \end{aligned}(19)
ukp=ρkp1⎝⎛ΔXkpΔYkpΔZkp⎠⎞ukq=ρkq1⎝⎛ΔXkqΔYkqΔZkq⎠⎞usp=ρsp1⎝⎛ΔXspΔYspΔZsp⎠⎞usq=ρsq1⎝⎛ΔXsqΔYsqΔZsq⎠⎞(19)
并将基线向量符号b,表示成如式(20)所示的坐标增量形式。
b
=
(
Δ
X
k
s
Δ
Y
k
s
Δ
Z
k
s
)
(
20
)
\begin{aligned} b= \begin{pmatrix} \Delta X_{ks} \\ \Delta Y_{ks} \\ \Delta Z_{ks} \end{pmatrix}\\ \end{aligned}(20)
b=⎝⎛ΔXksΔYksΔZks⎠⎞(20)
将式(19)和式(20)代入,则式(18)可以写成式(21)形式。式(21)即为基线方程。
∇
Δ
φ
k
s
p
q
+
∇
Δ
ε
k
s
p
q
=
f
2
c
⋅
(
Δ
X
k
p
ρ
k
p
+
Δ
X
s
p
ρ
s
p
−
Δ
X
k
q
ρ
k
q
−
Δ
X
s
q
ρ
s
q
Δ
Y
k
p
ρ
k
p
+
Δ
Y
s
p
ρ
s
p
−
Δ
Y
k
q
ρ
k
q
−
Δ
Y
s
q
ρ
s
q
Δ
Z
k
p
ρ
k
p
+
Δ
Z
s
p
ρ
s
p
−
Δ
Z
k
q
ρ
k
q
−
Δ
Z
s
q
ρ
s
q
)
(
Δ
X
k
s
Δ
Y
k
s
Δ
Z
k
s
)
−
∇
Δ
N
k
s
p
q
=
f
2
c
⋅
(
Δ
X
k
p
ρ
k
p
+
Δ
X
s
p
ρ
s
p
−
Δ
X
k
q
ρ
k
q
−
Δ
X
s
q
ρ
s
q
)
Δ
X
k
s
+
f
2
c
⋅
(
Δ
Y
k
p
ρ
k
p
+
Δ
Y
s
p
ρ
s
p
−
Δ
Y
k
q
ρ
k
q
−
Δ
Y
s
q
ρ
s
q
)
Δ
Y
k
s
+
f
2
c
⋅
(
Δ
Z
k
p
ρ
k
p
+
Δ
Z
s
p
ρ
s
p
−
Δ
Z
k
q
ρ
k
q
−
Δ
Z
s
q
ρ
s
q
)
Δ
Z
k
s
−
∇
Δ
N
k
s
p
q
=
a
k
s
p
q
Δ
X
k
s
+
b
k
s
p
q
Δ
Y
k
s
+
c
k
s
p
q
Δ
Z
k
s
−
∇
Δ
N
k
s
p
q
(
21
)
\begin{aligned} &\nabla\Delta\varphi_{ks}^{pq} +\nabla \Delta\varepsilon_{ks}^{pq}\\ &=\frac f{2c}\cdot \begin{pmatrix} \frac{\Delta X_k^p}{\rho_k^p}+ \frac{\Delta X_s^p}{\rho_s^p}- \frac{\Delta X_k^q}{\rho_k^q}- \frac{\Delta X_s^q}{\rho_s^q} \\ \frac{\Delta Y_k^p}{\rho_k^p}+ \frac{\Delta Y_s^p}{\rho_s^p}- \frac{\Delta Y_k^q}{\rho_k^q}- \frac{\Delta Y_s^q}{\rho_s^q} \\ \frac{\Delta Z_k^p}{\rho_k^p}+ \frac{\Delta Z_s^p}{\rho_s^p}- \frac{\Delta Z_k^q}{\rho_k^q}- \frac{\Delta Z_s^q}{\rho_s^q} \\ \end{pmatrix} \begin{pmatrix} \Delta X_{ks} \\ \Delta Y_{ks} \\ \Delta Z_{ks} \end{pmatrix} -\nabla\Delta N_{ks}^{pq} \\ &=\frac f{2c}\cdot \begin{pmatrix} \frac{\Delta X_k^p}{\rho_k^p}+ \frac{\Delta X_s^p}{\rho_s^p}- \frac{\Delta X_k^q}{\rho_k^q}- \frac{\Delta X_s^q}{\rho_s^q} \end{pmatrix}\Delta X_{ks}\\ &+\frac f{2c}\cdot \begin{pmatrix} \frac{\Delta Y_k^p}{\rho_k^p}+ \frac{\Delta Y_s^p}{\rho_s^p}- \frac{\Delta Y_k^q}{\rho_k^q}- \frac{\Delta Y_s^q}{\rho_s^q} \end{pmatrix}\Delta Y_{ks}\\ &+\frac f{2c}\cdot \begin{pmatrix} \frac{\Delta Z_k^p}{\rho_k^p}+ \frac{\Delta Z_s^p}{\rho_s^p}- \frac{\Delta Z_k^q}{\rho_k^q}- \frac{\Delta Z_s^q}{\rho_s^q} \end{pmatrix}\Delta Z_{ks}-\nabla\Delta N_{ks}^{pq} \\ &=a_{ks}^{pq}\Delta X_{ks}+b_{ks}^{pq}\Delta Y_{ks}+c_{ks}^{pq}\Delta Z_{ks}-\nabla\Delta N_{ks}^{pq}\\ \end{aligned}\quad \quad(21)
∇Δφkspq+∇Δεkspq=2cf⋅⎝⎜⎜⎛ρkpΔXkp+ρspΔXsp−ρkqΔXkq−ρsqΔXsqρkpΔYkp+ρspΔYsp−ρkqΔYkq−ρsqΔYsqρkpΔZkp+ρspΔZsp−ρkqΔZkq−ρsqΔZsq⎠⎟⎟⎞⎝⎛ΔXksΔYksΔZks⎠⎞−∇ΔNkspq=2cf⋅(ρkpΔXkp+ρspΔXsp−ρkqΔXkq−ρsqΔXsq)ΔXks+2cf⋅(ρkpΔYkp+ρspΔYsp−ρkqΔYkq−ρsqΔYsq)ΔYks+2cf⋅(ρkpΔZkp+ρspΔZsp−ρkqΔZkq−ρsqΔZsq)ΔZks−∇ΔNkspq=akspqΔXks+bkspqΔYks+ckspqΔZks−∇ΔNkspq(21)
02基线方程的解算
接下来,我们讲解利用基线方程,解算基线向量的方法。
将式(21)的基线方程写成误差方程的形式,如式(22)所示。
v
k
s
p
q
=
a
k
s
p
q
Δ
X
^
k
s
+
b
k
s
p
q
Δ
Y
^
k
s
+
c
k
s
p
q
Δ
Z
^
k
s
−
∇
Δ
N
^
k
s
p
q
−
∇
Δ
φ
k
s
p
q
=
a
k
s
p
q
δ
X
^
k
s
+
b
k
s
p
q
δ
Y
^
k
s
+
c
k
s
p
q
δ
Z
^
k
s
−
∇
Δ
N
^
k
s
p
q
−
w
k
s
p
q
(
22
)
v_{ks}^{pq}=a_{ks}^{pq}\Delta\hat{X}_{ks}+b_{ks}^{pq}\Delta\hat{Y}_{ks}+c_{ks}^{pq}\Delta\hat{Z}_{ks}-\nabla\Delta \hat N_{ks}^{pq}-\nabla\Delta \varphi_{ks}^{pq}\\= a_{ks}^{pq}\delta\hat{X}_{ks}+b_{ks}^{pq}\delta\hat{Y}_{ks}+c_{ks}^{pq}\delta\hat{Z}_{ks}-\nabla\Delta \hat N_{ks}^{pq}-w_{ks}^{pq}(22)
vkspq=akspqΔX^ks+bkspqΔY^ks+ckspqΔZ^ks−∇ΔN^kspq−∇Δφkspq=akspqδX^ks+bkspqδY^ks+ckspqδZ^ks−∇ΔN^kspq−wkspq(22)
式中:
Δ
X
^
k
s
=
Δ
X
^
k
s
0
+
δ
X
^
k
s
Δ
X
^
k
s
0
=
X
^
s
0
−
X
^
k
0
Δ
Y
^
k
s
=
Δ
Y
^
k
s
0
+
δ
Y
^
k
s
Δ
Y
^
k
s
0
=
Y
^
s
0
−
Y
^
k
0
Δ
Z
^
k
s
=
Δ
Z
^
k
s
0
+
δ
Z
^
k
s
Δ
Z
^
k
s
0
=
Z
^
s
0
−
Z
^
k
0
w
k
s
p
q
=
a
k
s
p
q
Δ
X
^
k
s
0
+
b
k
s
p
q
Δ
Y
^
k
s
0
+
c
k
s
p
q
Δ
Z
^
k
s
0
−
∇
Δ
φ
k
s
p
q
\Delta\hat{X}_{ks}=\Delta\hat{X}_{ks0}+\delta\hat{X}_{ks}\quad \Delta\hat{X}_{ks0}=\hat{X}_{s0}-\hat{X}_{k0}\\ \Delta\hat{Y}_{ks}=\Delta\hat{Y}_{ks0}+\delta\hat{Y}_{ks}\quad \Delta\hat{Y}_{ks0}=\hat{Y}_{s0}-\hat{Y}_{k0}\\ \Delta\hat{Z}_{ks}=\Delta\hat{Z}_{ks0}+\delta\hat{Z}_{ks}\quad \Delta\hat{Z}_{ks0}=\hat{Z}_{s0}-\hat{Z}_{k0}\\ w_{ks}^{pq}=a_{ks}^{pq}\Delta\hat{X}_{ks0}+b_{ks}^{pq}\Delta\hat{Y}_{ks0}+c_{ks}^{pq}\Delta\hat{Z}_{ks0}-\nabla\Delta \varphi_{ks}^{pq}
ΔX^ks=ΔX^ks0+δX^ksΔX^ks0=X^s0−X^k0ΔY^ks=ΔY^ks0+δY^ksΔY^ks0=Y^s0−Y^k0ΔZ^ks=ΔZ^ks0+δZ^ksΔZ^ks0=Z^s0−Z^k0wkspq=akspqΔX^ks0+bkspqΔY^ks0+ckspqΔZ^ks0−∇Δφkspq
如果以某一个卫星作为参考卫星,某一历元在测站i和测站j上,同时观测了S个卫星,则可列出(S-1)个形如式(22)的误差方程,相应地要引入(S-1)个双差模糊度未知数,再考虑到有3个基线向量未知数,则该历元共有s+2个未知数,未知数个数大于方程数,此时,误差方程组不可解。若测站i和测站j对所有S个卫星进行连续观测,随着历元数的增加,误差方程数也在增加,但未知数没有变。假设观测了n个历元,则总共有n*(s-1)个误差方程,方程数大于未知数个数,进而实现了误差方程组可解。
将n个历元的误差方程写成如式(23)的矩阵形式。已知观测值的方差矩阵为D。
V
=
A
X
^
−
N
^
−
W
=
(
A
−
E
)
(
X
^
N
^
)
−
W
(
23
)
V=A\hat X-\hat N-W=(A-E) \begin{pmatrix} \hat X \\ \hat N \end{pmatrix}-W(23)
V=AX^−N^−W=(A−E)(X^N^)−W(23)
式中:
{
V
=
(
v
1
,
v
2
,
…
,
v
s
)
T
X
^
=
(
δ
X
,
δ
Y
,
δ
Z
)
T
N
^
=
(
∇
Δ
N
^
1
,
∇
Δ
N
^
2
,
⋯
,
∇
Δ
N
^
s
−
1
)
T
W
=
(
w
1
,
w
2
,
…
,
w
o
)
T
\begin{cases} V=(v_1,v_2,\dots,v_s)^T \\ \hat X = (\delta X,\delta Y,\delta Z)^T \\ \hat N=(\nabla\Delta\hat N_1,\nabla\Delta\hat N_2,\cdots,\nabla\Delta\hat N_{s-1})^T\\ W=(w_1,w_2,\dots,w_o)^T \end{cases}
⎩⎪⎪⎪⎨⎪⎪⎪⎧V=(v1,v2,…,vs)TX^=(δX,δY,δZ)TN^=(∇ΔN^1,∇ΔN^2,⋯,∇ΔN^s−1)TW=(w1,w2,…,wo)T
A
=
(
a
1
b
1
c
1
a
2
b
2
c
2
⋮
⋮
⋮
a
o
b
o
c
o
)
A=\begin{pmatrix} a_1&b_1&c_1 \\ a_2&b_2&c_2 \\ \vdots&\vdots&\vdots\\ a_o&b_o&c_o \end{pmatrix}
A=⎝⎜⎜⎜⎛a1a2⋮aob1b2⋮boc1c2⋮co⎠⎟⎟⎟⎞
则利用最小二乘,通过式(24)可得基线向量与双差模糊度在实数域的最优估值,通过式(25)可得相应的方差矩阵。
(
X
^
N
^
)
=
(
(
A
T
−
E
)
D
−
1
(
A
T
−
E
)
)
−
1
(
A
T
E
)
D
−
1
W
(
24
)
\begin{pmatrix} \hat X\\\hat N\end{pmatrix}=\left( \begin{pmatrix} A^T\\-E\end{pmatrix}D^{-1} \begin{pmatrix} A^T&-E\end{pmatrix}\right)^{-1} \begin{pmatrix} A^T\\E\end{pmatrix}D^{-1}W(24)
(X^N^)=((AT−E)D−1(AT−E))−1(ATE)D−1W(24)
(
D
X
^
X
^
D
X
^
N
^
D
N
^
X
^
D
N
^
N
^
)
=
(
(
A
T
−
E
)
D
−
1
(
A
T
−
E
)
)
−
1
(
25
)
\begin{pmatrix} D_{\hat X\hat X}&D_{\hat X\hat N}\\D_{\hat N\hat X}&D_{\hat N\hat N}\end{pmatrix}= \left(\begin{pmatrix} A^T\\-E\end{pmatrix}D^{-1} \begin{pmatrix} A^T&-E\end{pmatrix}\right)^{-1}(25)
(DX^X^DN^X^DX^N^DN^N^)=((AT−E)D−1(AT−E))−1(25)
对于模糊度,其真值应该为整数,这是一个重要的约束信息,但在前述解算过程中,并没有顾及到这一约束信息。因此,下一步,还需要利用模糊度的整数特性,进一步提高基线的解算精度。具体方法,我们将在后续内容中进行介绍。