文章目录
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, "°")