GPS从入门到放弃（七）、GPS卫星位置解算

A = A_sqrt**2 # 卫星轨道半长轴
print("A={}".format(A))

上一篇讲了开普勒轨道参数，根据这些参数就可以确定卫星的位置，这一篇我们来实际计算一下。

WGS-84基本参数

首先给出几个WGS-84坐标系中的基本参数：

Python代码如下：

# WGS-84基本参数
a = 6378137 # 基准椭球体长半径（m）
f = 1/298.257223563 # 基准椭球体扁率
Omega_e_Dot = 7.2921151467e-5 # 地球自转角速度（rad/s）
mu = 3.986005e14 # 地球引力常数GM（m^3/s^2）
c = 2.99792458e8 # 真空中的光速（m/s）

星历参数

从网站 ftp://cddis.nasa.gov/gnss/data 下载2019年国庆（10月1日）当天的星历数据，可以得到星历参数的值，Python代码如下：

# 星历参数
t_oe = 2.016000000000E+05
A_sqrt = 5.153681812286E+03
e = 1.475233526435E-02
i_0 = 9.590228562257E-01
Omega_0 = -1.091936976129E-01
omega = 6.837269280624E-01
M_0 = 1.699075304872E+00
Delta_n = 4.599120143243E-09
i_Dot = -3.957307694893E-10
Omega_Dot = -8.244629136182E-09
Cuc = -5.902722477913E-06
Cus = 9.264796972275E-06
Crc = 2.046875000000E+02
Crs = -1.155625000000E+02
Cic = -3.259629011154E-07
Cis = 5.774199962616E-08

下面开始计算。

计算卫星轨道半长轴

A = A_sqrt**2 # 卫星轨道半长轴
print("A={}".format(A))

可得 A=26560436.222287513（米）

计算规化时间

A = A_sqrt**2 # 卫星轨道半长轴
t = 1.993680000000E+05
t_k = t - t_oe
if t_k > 302400:
    t_k -= 604800
if t_k < -302400:
    t_k += 604800
if -7200 <= t_k and t_k <= 7200:
    print("t_k={}".format(t_k))
else:
    print("time t={} is not valid".format(t))

代码如下：

# 平近点角
M_k = M_0 + n*t_k
if M_k < 0:
    M_k += 2*math.pi
if M_k > 2*math.pi:
    M_k -= 2*math.pi
print("M_k={}".format(M_k))

# 偏近点角
E_old = M_k
E_new = M_k + e*math.sin(E_old)
i = 1
while abs(E_new - E_old)>1e-8:
    print("i={},E={}".format(i, E_new))
    E_old = E_new
    E_new = M_k + e*math.sin(E_old)
    i += 1
    if (i>10):
        break

E_k = E_new
print("E_k={}".format(E_k))

# 真近点角
cosNu_k = (math.cos(E_k) - e) / (1 - e*math.cos(E_k))
sinNu_k = (math.sqrt(1-e**2)*math.sin(E_k)) / (1 - e*math.cos(E_k))
print("cosNu_k={}".format(cosNu_k))
print("sinNu_k={}".format(sinNu_k))

if cosNu_k == 0:
    if sinNu_k > 0:
        Nu_k = math.pi/2
    else:
        Nu_k = -math.pi/2
else:
    Nu_k = math.atan(sinNu_k/cosNu_k)

if cosNu_k < 0:
    if sinNu_k >= 0:
        Nu_k += math.pi
    else:
        Nu_k -= math.pi

print("Nu_k={}".format(Nu_k))

根据以下关系计算

代码如下：

delta_u_k = Cus*math.sin(2*Phi_k) + Cuc*math.cos(2*Phi_k)
delta_r_k = Crs*math.sin(2*Phi_k) + Crc*math.cos(2*Phi_k)
delta_i_k = Cis*math.sin(2*Phi_k) + Cic*math.cos(2*Phi_k)
print("delta_u_k={}".format(delta_u_k))
print("delta_r_k={}".format(delta_r_k))
print("delta_i_k={}".format(delta_i_k))

u_k = Phi_k + delta_u_k
r_k = A*(1-e*math.cos(E_k)) + delta_r_k
i_k = i_0 + i_Dot*t_k + delta_i_k
print("u_k={}".format(u_k))
print("r_k={}".format(r_k))
print("i_k={}".format(i_k))

代码如下：

x_p_k = r_k * math.cos(u_k)
y_p_k = r_k * math.sin(u_k)
print("x_p_k={}".format(x_p_k))
print("y_p_k={}".format(y_p_k))

Omega_k = Omega_0 + (Omega_Dot - Omega_e_Dot)*t_k - Omega_e_Dot*t_oe
print("Omega_k={}".format(Omega_k))

x_k = x_p_k*math.cos(Omega_k) - y_p_k*math.cos(i_k)*math.sin(Omega_k)
y_k = x_p_k*math.sin(Omega_k) + y_p_k*math.cos(i_k)*math.cos(Omega_k)
z_k = y_p_k*math.sin(i_k)
print("x_k={}".format(x_k))
print("y_k={}".format(y_k))
print("z_k={}".format(z_k))