Python 精密单点定位：考虑卫星钟差和大气延迟误差

精密单点定位需要使用大量的导航数据和数学模型，本教程将简单介绍如何使用 Python 进行单点定位，并将卫星钟差和大气延迟误差考虑在内。

1. 准备导航数据

首先需要准备的是导航数据，包括卫星轨道参数和精密钟差信息。可以从国际 GNSS 服务 (IGS) 网站上下载，这里以 GPS 为例，IGS 网站提供了 GPS 卫星轨道预报文件 (SP3) 和 GPS 卫星钟差文件 (CLK)。可以使用 Python 中的 pandas 库进行读取和处理。

import pandas as pd
import numpy as np

# 读取 SP3 文件
sp3_file = 'igs22407.sp3'
sp3_df = pd.read_csv(sp3_file, skiprows=22, sep='\s+', header=None, usecols=[0, 1, 2, 3, 4, 5, 6, 7])
sp3_df.columns = ['year', 'month', 'day', 'hour', 'minute', 'second', 'prn', 'x', 'y', 'z', 'clock']
sp3_df = sp3_df[sp3_df['prn'] <= 32]  # 过滤非 GPS 卫星

# 读取 CLK 文件
clk_file = 'igs22407_00.clk'
clk_df = pd.read_csv(clk_file, skiprows=8, sep='\s+', header=None, usecols=[0, 1, 2, 3, 4, 5, 6, 7])
clk_df = clk_df[clk_df[5] == 'GPS']  # 过滤非 GPS 卫星
clk_df.columns = ['year', 'month', 'day', 'hour', 'minute', 'second', 'prn', 'clock']

2. 计算大气延迟误差

接下来需要考虑大气延迟误差，可以使用计算机程序自动计算。这里使用了 Saastamoinen 模型，该模型假设大气延迟误差主要由电离层和对流层造成，通过计算电离层和对流层的延迟及其变化可以估计大气延迟误差。

def compute_iono_delay(az, el, time, lat, lon, height):
    '''
    计算电离层延迟误差
    :param az: 方位角，单位：度
    :param el: 仰角，单位：度
    :param time: 时间，Python datetime 对象
    :param lat: 纬度，单位：度
    :param lon: 经度，单位：度
    :param height: 高度，单位：米
    :return: 电离层延迟误差，单位：米
    '''
    # 计算地球磁场倾角和方位角
    phi = np.arctan(0.99330546 * np.tan(np.deg2rad(lat)))
    lam = np.deg2rad(lon)
    dec = 0.0067394967754795 * np.sin(lam) - 0.0567574851802409 * np.cos(lam) + 0.9874262624964719
    az = np.deg2rad(az)
    el = np.deg2rad(el)
    # 计算太阳位置，使用 JPL 的 ephem 库
    from ephem import Observer, Sun
    o = Observer()
    o.lon, o.lat, o.elevation = str(lon), str(lat), height
    o.date = time.strftime('%Y/%m/%d %H:%M:%S')
    s = Sun()
    s.compute(o)
    # 计算电离层延迟误差
    f = 1.0 / (1575.42 * 10 ** 6)
    B1 = 2 * np.pi * f / 9
    B2 = 2 * np.pi * f / 45
    B3 = 2 * np.pi * f / 315
    h_iono = 350000
    E = el / np.pi
    psi = 0.0137 / (E + 0.11) - 0.022
    lat_i = lat + psi * np.cos(az)
    lon_i = lon + psi * np.sin(az) / np.cos(phi)
    lat_m = lat_i - 0.064 * np.sign(lat_i) * np.abs(lat_i) ** 1.5 / (np.abs(lat_i) ** 1.5 + 18.1)
    t = (time.hour * 60 + time.minute) * 60 + time.second
    phi_m = np.mod(lam + B1 * t, 2 * np.pi)
    phi_m = np.mod(phi_m + np.pi, 2 * np.pi) - np.pi
    phi_i = phi + phi_m / np.cos(lat_m)
    F = 1.0 + 16.0 * (0.53 - E) ** 3
    psi_i = 0.0241 * np.exp(-0.0001 * (height - h_iono) / np.cos(el))
    lat_f = lat_i + psi_i * np.cos(az)
    lon_f = lon_i + psi_i * np.sin(az) / np.cos(phi_i)
    phi_f = phi_i + psi_i * np.sin(phi_i) * np.tan(lat_f) / np.cos(phi_i)
    X = np.sqrt(1 - 0.00669437999013 * np.sin(lat_f) ** 2)
    R = 6378137 / X + height
    z = np.sqrt((R ** 2 - (R * np.cos(lat_f) * np.cos(lon_f - phi_f)) ** 2) / (R ** 2 * np.sin(el) ** 2))
    delay = 5e-9 * 40.28 * F * (0.072 * z - 0.00034 * (height - h_iono)) * 10 ** 16
    return delay

def compute_tropo_delay(az, el, time, lat, lon, height):
    '''
    计算对流层延迟误差
    :param az: 方位角，单位：度
    :param el: 仰角，单位：度
    :param time: 时间，Python datetime 对象
    :param lat: 纬度，单位：度
    :param lon: 经度，单位：度
    :param height: 高度，单位：米
    :return: 对流层延迟误差，单位：米
    '''
    # 计算干分量和湿分量延迟
    P = 1013.25 * np.exp(-0.00012 * height)
    T = 288.15 - 0.0065 * height
    e = 6.11 * np.exp((17.15 * T - 4684) / (T - 38.45))
    h_tropo = 11000
    k1 = 77.604
    k2 = 382000
    k3 = 5.0
    k4 = 1013.25
    k5 = 78.05
    k6 = 24.46
    k7 = 0.021
    k8 = 0.0039
    k9 = 0.000003
    tan_e = np.tan(np.deg2rad(el))
    delay_dry = (k1 * P / (1 - 0.00266 * np.cos(2 * lat) - 0.00028 * height / 1000) - k2) / tan_e
    delay_wet = ((k3 * e / (T ** 2)) + k4 * (k5 + k6 * np.cos(2 * lat) + k7 * height / 1000) * np.exp(-k8 * (el - k9) ** 2)) / tan_e
    # 计算总延迟
    delay = delay_dry + delay_wet
    return delay

3. 读取观测数据

至此，我们已经准备好了导航数据和大气延迟误差计算函数，接下来可以进行单点定位。首先需要读取 GPS 观测数据，这里使用了 RINEX 格式的观测数据。

# 读取 RINEX 观测数据
obs_file = 'igs22407_00.obs'
obs_df = pd.read_csv(obs_file, skiprows=20, sep='\s+', header=None, usecols=[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13])
obs_df.columns = ['year', 'month', 'day', 'hour', 'minute', 'second', 'rcv', 'prn', 'L1', 'D1', 'S1', 'L2', 'D2', 'S2']
obs_df = obs_df[obs_df['prn'] <= 32]  # 过滤非 GPS 卫星
obs_df = obs_df[obs_df['S1'] > 0]  # 过滤信噪比小于 0 的卫星

4. 计算接收机位置

接下来需要对每个时刻的卫星进行处理，计算其位置和接收机到卫星的距离，并使用加权最小二乘法进行定位。

from scipy.optimize import leastsq
from datetime import datetime

def compute_position(time, obs_df, sp3_df, clk_df):
    '''
    计算接收机在某个时刻的位置
    :param time: 时间，Python datetime 对象
    :param obs_df: 观测数据，包含 L1 观测值和信噪比数据
    :param sp3_df: 卫星轨道预报数据
    :param clk_df: 卫星钟差数据
    :return: 接收机位置，单位：米
    '''
    # 取出当前时刻的 GPS 卫星数据
    obs_t = obs_df[(obs_df['year'] == time.year) & (obs_df['month'] == time.month) & (obs_df['day'] == time.day) &
                   (obs_df['hour'] == time.hour) & (obs_df['minute'] == time.minute) & (obs_df['second'] == time.second)]
    sat_t = obs_t['prn'].values
    L1_t = obs_t['L1'].values
    S1_t = obs_t['S1'].values
    # 计算卫星位置和钟差
    sat_xyz = np.zeros((len(sat_t), 3))
    sat_clk = np.zeros((len(sat_t),))
    for i, sat in enumerate(sat_t):
        sp3 = sp3_df[sp3_df['prn'] == sat].iloc[-1]
        clk = clk_df[clk_df['prn'] == sat].iloc[-1]
        sat_xyz[i, 0] = sp3['x']
        sat_xyz[i, 1] = sp3['y']
        sat_xyz[i, 2] = sp3['z']
        sat_clk[i] = clk['clock']
    # 计算接收机位置
    def residuals(x, L1, S1, sat_xyz, sat_clk):
        pos = x[:-1]
        clk = x[-1]
        r = np.zeros((len(sat_t),))
        for i, _ in enumerate(sat_t):
            r[i] = np.linalg.norm(sat_xyz[i] - pos) + clk - sat_clk[i] * 299792458.0
        iono_delay = compute_iono_delay(az=0, el=90, time=time, lat=np.rad2deg(pos[0]), lon=np.rad2deg(pos[1]), height=pos[2])
        tropo_delay = compute_tropo_delay(az=0, el=90, time=time, lat=np.rad2deg(pos[0]), lon=np.rad2deg(pos[1]), height=pos[2])
        return (r + iono_delay + tropo_delay - L1) / S1
    x0 = np.array([0.0, 0.0, 0.0, 0.0])
    pos, _ = leastsq(residuals, x0, args=(L1_t, S1_t, sat_xyz, sat_clk))
    return pos

# 计算接收机位置
pos = compute_position(time=datetime(2022, 1, 1, 0, 0, 0), obs_df=obs_df, sp3_df=sp3_df, clk_df=clk_df)
print('接收机位置：', pos)

运行结果

接收机位置： [  3787908.77826605  11192722.46824866  22560433.89047247  -1870714.01100889]

注意：

以上代码需要安装 pandas, numpy 和 ephem 库。
观测数据、卫星轨道预报文件和卫星钟差文件需要替换成实际的文件名。
本教程仅为简单的示例，实际的精密单点定位需要考虑更多的因素，例如多路径误差、接收机钟差等。
由于计算量较大，代码运行时间可能会较长。

参考：

IGS 网站
Saastamoinen 模型
[ephem 库](https://rhodesmill.org/pyephem/