卫星高度角和方位角的计算
文章目录
Part.I IntroductionPart.II 计算方法Chap.I 概念Chap.II 计算流程
Part.III 关于是否要考虑坐标误差的探讨Chap.I 基础数据Chap.II 计算流程Chap.III 考虑坐标误差的情况Chap.IV 部分代码
Reference
Part.I Introduction
本文将介绍如何计算卫星的高度角和方位角。
在了解下面的内容之前，首先需要了解地心地固坐标系（XYZ），站心坐标系（ENU），大地坐标系（BLH），以及高度角、方位角的一些概念。下面的链接（Reference）中已经有详细的介绍，这里就不再赘述。
Part.II 计算方法
Chap.I 概念
如下图所示：卫星在太空中飞，地球上的测站位于下图左边的空心小圆处，地球半径为
r
e
r_e
re，卫星到地心的距离为
r
s
r_s
rs，卫星到测站间的距离为
d
d
d。从测站出发作地球的一个切面（在二维图形上看就是上面的那条垂线）那么卫星的高度角就是上图中的
E
I
EI
EI 。
把上面的“那条线”（也就是切面）拿出来，如上面右图所示（随手画的，很抽象）。测站正北方向为
N
N
N 方向，正东方向为
E
E
E 方向，测站正上方向为
U
U
U 方向。图中的橙线为站星连线在测站坐标系
N
O
E
NOE
NOE 平面的投影，投影线与
N
N
N 方向的夹角（图中的
A
Z
AZ
AZ）即为方位角。
Chap.II 计算流程
假设卫星在
t
0
t_0
t0 时刻播发的信号在
t
1
t_1
t1 时刻被卫星接收到，卫星在
t
0
t_0
t0 时刻的坐标为
(
x
0
s
,
y
0
s
,
z
0
s
)
(x^s_0,y^s_0,z^s_0)
(x0s,y0s,z0s)；测站在
t
1
t_1
t1 时刻的坐标为
(
x
r
1
,
y
r
1
,
z
r
1
)
(x_{r1},y_{r1},z_{r1})
(xr1,yr1,zr1)。
值得注意的是，上面的坐标都是地心地固系的坐标（在信号传播的过程中，地球一直在旋转，旋转的角度为
(
ρ
/
c
)
∗
ω
(\rho/c)*\omega
(ρ/c)∗ω『 传播距离 /光速 * 地球自转角速度』），换言之，坐标系在变化。
想要计算卫星的高度角，首先需要进行坐标系的统一（地球自转改正）。一般情况下，都是在信号接收时刻的地心地固系进行数据处理的，
下面简要写一下计算过程。
§1: 根据测站坐标与卫星坐标计算信号传播距离
ρ
\rho
ρ，进而计算出传播时间
t
=
(
ρ
/
c
)
t=(\rho/c)
t=(ρ/c)，需要改正的角度为
ϕ
=
(
ρ
/
c
)
∗
ω
\phi=(\rho/c)*\omega
ϕ=(ρ/c)∗ω。
§2: 将卫星坐标左乘一个旋转矩阵（因为相当于是绕
z
z
z 轴旋转）
R
z
(
θ
)
=
[
c
o
s
θ
−
s
i
n
θ
0
s
i
n
θ
c
o
s
θ
0
0
0
1
]
R_z(\theta)=\left[ \begin{array}{ccc} cos\theta & -sin\theta & 0 \\ sin\theta & cos\theta & 0 \\ 0 & 0 & 1 \\ \end{array}\right]
Rz(θ)=
cosθsinθ0−sinθcosθ0001
得到卫星在
t
1
t_1
t1 时刻的地心地固系下的坐标
(
x
1
s
,
y
1
s
,
z
1
s
)
(x_{1}^s,y_{1}^s,z_{1}^s)
(x1s,y1s,z1s)。
§3: 然后以
(
x
r
1
,
y
r
1
,
z
r
1
)
(x_{r1},y_{r1},z_{r1})
(xr1,yr1,zr1) 为站心坐标，计算出卫星在站心坐标系中的坐标
(
n
1
s
,
e
1
s
,
u
1
s
)
(n_1^s,e_1^s,u_1^s)
(n1s,e1s,u1s)。
§4: 那么卫星高度角
E
E
E 与
(
n
1
s
,
e
1
s
,
u
1
s
)
(n_1^s,e_1^s,u_1^s)
(n1s,e1s,u1s) 存在如下关系(将
(
n
1
s
,
e
1
s
,
u
1
s
)
(n_1^s,e_1^s,u_1^s)
(n1s,e1s,u1s) 简单记作
(
n
,
e
,
u
)
(n,e,u)
(n,e,u))：
c
o
s
E
=
n
2
+
e
2
n
2
+
e
2
+
u
2
cosE=\sqrt{\frac{n^{2}+e^2}{n^{2}+e^2+u^2}}
cosE=n2+e2+u2n2+e2
根据这个公式就可以计算出高度角。
§5: 卫星的方位角
A
A
A 与
(
n
,
e
,
u
)
(n,e,u)
(n,e,u) 存在如下关系：
t
a
n
A
=
e
n
tanA=\frac{e}{n}
tanA=ne根据这个公式可以计算出方位角。
上面的计算流程大致一看没什么问题，但是仔细考虑的话会存在如下几个问题： 1、一般测站的坐标是通过 SPP 简单计算出来，误差很大，卫星的坐标可能是通过广播星历得到的，误差也很大，要不要考虑坐标误差呢？ 2、如果要考虑坐标误差的话，第一次算出来的站星距离根本不准，进而得到的地球自转角度也不准，这里面是否还要迭代？
Part.III 关于是否要考虑坐标误差的探讨
直观感受影响不是很大，下面拿个具体算例计算一下。
Chap.I 基础数据
导航卫星到地面的距离大约为 20,200 km，现假设一颗星的坐标为
G03 12712.882254 23247.798196 -2637.709427
接收机坐标
-2267752.0605993434, 5009151.1456511570, 3221301.4797024932
常量
OMEGA = 7292115.1467e-11 # rad/s
CLIGHT = 2.99792458e+8 # m/s
Chap.II 计算流程
卫星至接收机的距离（两万四千多千米）
distance = np.linalg.norm(sat_crd - rec_crd)
# 卫星至接收机的距离: 24318627.829295974 m
信号传播时间
tau = distance / CLIGHT # [s]
# 信号传播时间: 0.08111821088339712 s
地球自转角度
ang = OMEGA * tau # [rad]
ang_deg = ang * 180 / G_PI # [deg]
# 地球旋转角度: 0.00033891790536376496 °
卫星坐标改正量
delta_X_sat = OMEGA * ang * sat_crd[1]
delta_Y_sat = -OMEGA * ang * sat_crd[0]
# 卫星坐标改正量 (delta_X_sat,delta_Y_sat): (0.010027836078424127, -0.005483646160923473)
m
经地球自转改正之后的卫星坐标
sat_crd_new = sat_crd + np.array([delta_X_sat, delta_Y_sat, 0])
# 经地球自转改正后的卫星坐标: [12712882.26402784 23247798.19051635 -2637709.427] m
卫星坐标转为站心坐标
sat_enu = bc.XYZ2ENU(sat_crd_new.tolist(),rec_crd.tolist()) # e,n,u
# 以接收机为站心的卫星站心坐标: [-21169312.671131082, -10348729.992622295, 6013289.297207948] m
计算卫星高度角和方位角
e,n,u = sat_enu[0], sat_enu[1], sat_enu[2]
cos_ELE = math.sqrt((n**2+e**2)/(n**2+e**2+u**2))
tan_AZI = e/n
ELE = math.acos(cos_ELE)
AZI = math.atan(tan_AZI)
ELE_deg = ELE * 180 / G_PI
AZI_deg = AZI * 180 / G_PI
# 卫星高度角为: 14.31607742027866 °
# 卫星方位角为: 63.948059502576534 °
Chap.III 考虑坐标误差的情况
如果把初始接收机的坐标加
(
100
,
100
,
100
)
(100,100,100)
(100,100,100)，按照同样的流程，得到的卫星的高度角和方位角为
卫星高度角为: 14.316743335284418 °
卫星方位角为: 63.94682589822006 °
由此可见，坐标误差对于高度角、方位角的计算完全可以忽略不计（小数点后第三位）。因此，就不需要迭代，用接收机的概略坐标即可！
Chap.IV 部分代码
下面展示了测试所用的部分代码
import math
import numpy as np
def compute():
# Initail value
sat_crd = np.array([12712882.254, 23247798.196, -2637709.427])
rec_crd = np.array([-2267752.0605993434, 5009151.1456511570, 3221301.4797024932])
rec_crd += np.array([100,100,100])
OMEGA = 7292115.1467e-11 # [rad/s]
CLIGHT = 2.99792458e+8 # [m/s]
G_PI = 3.14159265358979311599796346854419e0
# Calculate
distance = np.linalg.norm(sat_crd - rec_crd)
tau = distance / CLIGHT # [s]
ang = OMEGA * tau # [rad]
ang_deg = ang * 180 / G_PI # [deg]
delta_X_sat = OMEGA * ang * sat_crd[1]
delta_Y_sat = -OMEGA * ang * sat_crd[0]
sat_crd_new = sat_crd + np.array([delta_X_sat, delta_Y_sat, 0])
sat_enu = bc.XYZ2ENU(sat_crd_new.tolist(),rec_crd.tolist()) # e,n,u
e,n,u = sat_enu[0], sat_enu[1], sat_enu[2]
cos_ELE = math.sqrt((n**2+e**2)/(n**2+e**2+u**2))
tan_AZI = e/n
ELE = math.acos(cos_ELE)
AZI = math.atan(tan_AZI)
ELE_deg = ELE * 180 / G_PI
AZI_deg = AZI * 180 / G_PI
# Output
print("卫星至接收机的距离: ", distance, "m")
print("信号传播时间: ", tau, "s")
print("地球旋转角度: ", ang_deg, "°")
print("经地球自转改正后的卫星坐标: ", sat_crd_new, "m")
print("以接收机为站心的卫星站心坐标: ", sat_enu, "m")
print("卫星高度角为: ", ELE_deg, "°")
print("卫星方位角为: ", AZI_deg, "°")
Reference
关于大地测量领域常用的角度知识汇总计算卫星高度角和方位角大地测量（含导航定位）中常用的坐标系统概念简介【GNSS】地球自转改正