卫星高度角和方位角的计算

Part.I Introduction

本文将介绍如何计算卫星的高度角和方位角。

在了解下面的内容之前，首先需要了解地心地固坐标系（XYZ），站心坐标系（ENU），大地坐标系（BLH），以及高度角、方位角的一些概念。下面的链接（Reference）中已经有详细的介绍，这里就不再赘述。

Part.II 计算方法

Chap.I 概念

如下图所示：卫星在太空中飞，地球上的测站位于下图左边的空心小圆处，地球半径为 $r_e$，卫星到地心的距离为 $r_s$，卫星到测站间的距离为 $d$。从测站出发作地球的一个切面（在二维图形上看就是上面的那条垂线）那么卫星的高度角就是上图中的 $EI$ 。

把上面的“那条线”（也就是切面）拿出来，如上面右图所示（随手画的，很抽象）。测站正北方向为 $N$ 方向，正东方向为 $E$ 方向，测站正上方向为 $U$ 方向。图中的橙线为站星连线在测站坐标系 $NOE$ 平面的投影，投影线与 $N$ 方向的夹角（图中的 $AZ$）即为方位角。

Chap.II 计算流程

假设卫星在 $t_0$ 时刻播发的信号在 $t_1$ 时刻被卫星接收到，卫星在 $t_0$ 时刻的坐标为 $(x_0^s,y_0^s,z_0^s)$；测站在 $t_1$ 时刻的坐标为 $(x_{r1},y_{r1},z_{r1})$。

值得注意的是，上面的坐标都是地心地固系的坐标（在信号传播的过程中，地球一直在旋转，旋转的角度为 $(\rho /c)\ast \omega$『 传播距离 /光速 * 地球自转角速度』），换言之，坐标系在变化。

想要计算卫星的高度角，首先需要进行坐标系的统一（地球自转改正）。一般情况下，都是在信号接收时刻的地心地固系进行数据处理的，

下面简要写一下计算过程。

• §1: 根据测站坐标与卫星坐标计算信号传播距离 $\rho$，进而计算出传播时间 $t=(\rho /c)$，需要改正的角度为 $\varphi =(\rho /c)\ast \omega$。

• §2: 将卫星坐标左乘一个旋转矩阵（因为相当于是绕 $z$ 轴旋转）
  $R_z(\theta )=\left[\begin{array}{ccc}cos\theta & -sin\theta & 0\\ sin\theta & cos\theta & 0\\ 0& 0& 1\end{array}\right]$
  得到卫星在 $t_1$ 时刻的地心地固系下的坐标 $(x_1^s,y_1^s,z_1^s)$。

• §3: 然后以 $(x_{r1},y_{r1},z_{r1})$ 为站心坐标，计算出卫星在站心坐标系中的坐标 $(n_1^s,e_1^s,u_1^s)$。

• §4: 那么卫星高度角 $E$ 与 $(n_1^s,e_1^s,u_1^s)$ 存在如下关系(将$(n_1^s,e_1^s,u_1^s)$ 简单记作 $(n,e,u)$)：
  $$cosE=\sqrt{\frac{n^2+e^2}{n^2+e^2+u^2}}$$ 根据这个公式就可以计算出高度角。

• §5: 卫星的方位角 $A$ 与 $(n,e,u)$ 存在如下关系：
  $$tanA=\frac{e}{n}$$根据这个公式可以计算出方位角。

---

上面的计算流程大致一看没什么问题，但是仔细考虑的话会存在如下几个问题：
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)$，按照同样的流程，得到的卫星的高度角和方位角为

卫星高度角为:  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】地球自转改正