精密星历计算卫星位置（Lagrange）

2301_77671395

于 2025-08-05 17:29:27 发布

阅读量587

点赞数 14

CC 4.0 BY-SA版权

文章标签： python 开发语言 GNSS SP3 拉格朗日插值法

本文链接：https://blog.csdn.net/2301_77671395/article/details/149916074

小编最近在练习编程数据处理，最近在学习GNSS方面知识时，关注到可以用精密星历计算卫星位置，方法选用拉格朗日插值，看看内插法和真实值之间会有怎样的差别，于是想着用编程实现一下，说干就干！

1 数据来源

首先需要获取事后星历，毕竟我这是想对比一下预测值和真实值之间的差异，所以非快速 / 超快速星历数据就不太合适了。

https://igs.org/ : International GNSS Sevrice,里面有许多IGS产品；

GPS Time Calculator:GPS Time Calculator:GPS Time Calculator,方便获取GPS时间，查看精密星历数据的头部信息，获取文件类型；在文件命名时可以输入年月日等时间，从而获取GPS周，查看日历，获取数据时间的星期数据等等，简化精密星历的命名,头部信息如下：

#dP2023 1 30 0 0 0.00000000 289 d+D IGS20 FIT AIUB

IGS

文件可以按照上面信息，第2247 GPS 周，星期1，命名为： IGS22471.sp3；

2 Lagrange方法介绍

精密星历是按照一定时间间隔来给出卫星在空间中的三维坐标、卫星钟改正数等信息。观测瞬间的卫星位置及运动速度可采用内插法求得。其中拉格朗日内插法被广泛使用，因为这种内插法速度快且易编程。拉格朗日插值公式非常简单：已知函数 $y = f(x)$ 的n+1 个节点 $x_{0},x_{1},x_{2},\cdots ,x_{n}$ 及其对应的函数值 $y_{0},y_{1},y_{2},\cdots ,y_{n}$ ，对于插值时间内任意一点x ，可以用下面的拉格朗日插值多项式来计算函数值：

$f(x) = \sum_{k=0}^{n}\prod_{i=0,i\neq k}^{n}(\frac{x-x_{i}}{x_{k-x_{i}}})y_{k}$

方便大家理解，给大家讲一个好笑的事情，《GNSS原理与应用》这门课期中考试时，有一道大题是考这个插值方法的，给了一组数据，让我们预测下一个节点对应的值，小编一看数据，不就是线性拟合吗，给了个线性表达式，就直接写出结果了，根本没往这方面想！如题：

X	0	1	2	4
Y	1	3	5	？

按照线性方程来解的话，其实就是 $y = 2x+1$ ，所以，x=4,y=9;

但是按照拉格朗日插值方法来进行计算，过程应该是如下：

$f(4)=(\frac{4-1}{0-1})\times (\frac{4-2}{0-2})\times 1+(\frac{4-0}{1-0})\times (\frac{4-2}{1-2})\times 3+(\frac{4-0}{2-0})\times (\frac{4-1}{2-1})\times 5$

$=3-24+30$

$=9$

这个公式计算的时候容易失误，不过记住，分母左边 $x_{k}$ 在同一个乘式中不变，最后要乘以其对应的 $f(x_{k})$ ,分式右边 $x_{i}$ 上下一致,最后就是 $i\neq k$ ;

3 编程实现

其实我编程比较简单，用的都是比较基础的函数和语法，入门新手啦，边学边成长！之前基本都是自己写个大概，运行报错就是求助AI，不过这次是我慢慢调试-运行-调试出来的，真的发现编程不易，代码写到最后，都想着，反正出图了，就这样吧，哈哈哈，一边躺平一边哒哒哒地敲键盘，好在最后运行成功。

3.1 解读卫星星历

我使用的精密星历是5min的更新频率，还挺神奇的，之前没有了解到，竟然更新的这么快！首先，我们先来解读一下卫星星历存储格式，根据文件内容可知

+开头，代表可观测的卫星数及其编号，如这个文件，表明可以观测到117颗卫星，G开头代表GPS卫星；

P开头，代表卫星轨道位置信息；

了解到文件格式后，就可以根据每一行的标识符来进行文件解读啦。

def parse_obs_file(filename):
    data = []  # 接收数据
    current_time = None  # 获取历元信息
    prn = []
    with open(filename,'r') as f:
        for line in f:
            line = line.strip()
            if not line:
                continue
            # 读取卫星PRN,以 + 开头的，同时避免误读取 ++(方便查看)
            if line.startswith('+') and not line.startswith("++"):
                # 从9位开始就是prn,每三位为一个卫星标识编号
                for i in range(9,len(line),3):
                    prn_id = line[i:i+3]
                    # 过滤空值和填充值
                    if prn_id and not prn_id.startswith('0'):
                        prn.append(prn_id)

            # 当前卫星历元数据,年月日时分秒
            if line.startswith('*'):
                parts = line.split()
                try:
                    current_time = datetime(
                        int(parts[1]),int(parts[2]),int(parts[3]),
                        int(parts[4]),int(parts[5]),int(float(parts[6]))
                    )
                except:
                    continue
            # 获取卫星轨道位置数据，prn，X,Y,Z（有必要可以获得钟差数据）
            if line.startswith('P'):
                dataline = line.split()
                try:
                    data.append([current_time,dataline[0][1:],
                                 float(dataline[1]),float(dataline[2]),float(dataline[3])])
                except ValueError:
                    continue

    # 返回各卫星轨道数据
    return pd.DataFrame(data,columns=["Time","PRN","X","Y","Z"]),prn

我顺便读取了卫星prn编号等信息，可以帮助我们看到可观测卫星的数量及编号，虽然好像没啥用，不过，锻炼了我数据提取能力嘛。同时，我们将星历数据以【历元，卫星prn号，此刻历元的X，Y ，Z】的格式创建了DataFrame类型数据，并将结果返回。

为了查看是否读取成功，小编在函数内部编写了一段测试代码，输出部分信息，帮助我们了解数据格式和内容。

print("-"*100)
print("卫星PRN：")
for i in range(1,len(prn)+1):
    print(prn[i-1],end=' ')
    if i%17==0:
        print("\n")
#打印一部分轨道数据:前八行数据
print("-"*100,"查看前八行轨道数据",sep="\n")
for i in range(8):
    print(data[i])
print("-" * 100)

结果如下：

......

结果表明文件读取和输出正常，就可以开始下一步操作。本来输出时间应该是精确到秒，但原始数据秒数为 0，这里会省略显示（但实际仍存在）。因为本次数据的更新频率是5min，所以秒速基本是0，不用特别在意。

3.2 将时间转为数值

插值时，建立的是y = f(t)的函数，所以datetime类型的数据转换为数值型可以方便计算，在此统一转换为秒。由于数值太大，所以我以卫星轨道首历元的时间为参考基准，减小数值。

# 将时间转换成数值，以便插值
def time_to_numeric(times):
    """
    为了降低大数据的处理复杂性，将时间转化为秒，
    因为只是预测一天的情况，所以时间定格到第一个历元时间数据为参考点
    :param times: 历元列表
    :return: 相对首历元的时间差
    """
    ref_time = times[0]
    diff_time = [(t-ref_time).total_seconds() for t in times]
    return diff_time

3.3 Lagrange插值方法实现

刚开始，我是想着按照书上公式来进行编写代码，然后上网搜了一下，推荐用移动窗口，但是不知道为什么误差竟然在1km范围内浮动，不知道什么原因，可能我用了30min的时间间隔。预测没那么准了，所以最后还是采用了书上的方法。

Remondi的研究表明，对GPS卫星而言，如果要精确到 $10^{-8}$ ，用30min的历元间隔和9阶内插已足够保证精度。内插的阶数选取也很重要，由于我们数据的时间间隔比较短，所以可以使用5阶的进行内插。

代码实现如下：

# 拉格朗日内插法
def lagrange(data,time):
    # 一般五阶的内插也能获得较好的结果
    """
    :param data: 已知函数值（如X,Y,Z），
    :param time: 对应历元时间，包含待预测历元
    :return: 下一个历元的预测值
    """
    sum = 0
    n = len(data)
    for k in range(n):
        multi = 1
        for i in range(n):
            if i==k :
                continue
            else:
                multi *= float(time[-1]-time[i])/(time[k]-time[i])
        sum += multi*data[k]
    return sum

4 代码测试与结果

为了减少数据处理的复杂度，我只选取了一个卫星多历元数据进行数据处理与分析，先对原始数据的轨道路径进行成图，大致了解卫星轨迹；

其实，这个图像包含了卫星轨道真实值路径和预测数据，但是由于差距在万级的距离上来看就有点微小了，所以还需要对预测值与真实值进行残差分析，直观地感受一下使用拉格朗日内插方法，在X,Y,Z方向上的差异情况；

# 测试一下代码
if __name__ == "__main__":
    file_path = r"./igs22471.sp3"
    datas,prn = parse_obs_file(file_path)
    # 将时间值进行简化
    datas["simply_time"] = pd.Series(time_to_numeric(datas["Time"]))
    # 为了简化计算量，仅用一个卫星的轨道实测数据开展实验,历元间隔为15min；
    data_G06 = datas[datas["PRN"]=="G06"]
    predict_X_G06 = []
    predict_Y_G06 = []
    predict_Z_G06 = []
    # 进行插值
    for i in range(5,len(data_G06)):
        predict_X_G06.append(lagrange(data_G06["X"][i-5:i].tolist(),
                                      data_G06["simply_time"][i-5:i+1].tolist()))
        predict_Y_G06.append(lagrange(data_G06["Y"][i-5:i].tolist(),
                                      data_G06["simply_time"][i-5:i+1].tolist()))
        predict_Z_G06.append(lagrange(data_G06["Z"][i-5:i].tolist(),
                                      data_G06["simply_time"][i-5:i+1].tolist()))
    # 创建3D图形
    fig = plt.figure(figsize=(10, 8))
    # 3D图展示实测点和预测点整体差异
    ax1 = fig.add_subplot(121, projection='3d')
    ax1.scatter(data_G06["X"][5:], data_G06["Y"][5:], data_G06["Z"][5:], color="#a8e6cf")
    ax1.scatter(predict_X_G06, predict_Y_G06, predict_Z_G06, color="#008ba3")
    # 设置坐标标签
    ax1.set_xlabel('X', fontsize=12)
    ax1.set_ylabel('Y', fontsize=12)
    ax1.set_zlabel('Z', fontsize=12)
    # 显示X,Y,Z坐标值随历元变化情况
    t = data_G06["simply_time"][5:]/3600
    ax2 = fig.add_subplot(222)
    ax2.scatter(t,data_G06["X"][5:]-predict_X_G06,color="#a8e6cf",s=8,label="X_diff")
    ax2.scatter(t, data_G06["Y"][5:] - predict_Y_G06, color="#cef0d5", s=8,label="Y_diff")
    ax2.scatter(t, data_G06["Z"][5:] - predict_Z_G06, color="#4bb89f", s=8,label="Z_diff")
    # 在y=0处添加一条水平直线，设置为虚线以便区分
    ax2.axhline(y=0, color='#4dd0e1', linestyle='--', alpha=0.7)
    # 设置X轴刻度与数值一一对应
    ticks=data_G06["simply_time"][5:]/3600
    ax2.set_xticks(ticks[::10])
    plt.xticks(fontsize=6,rotation=90)
    # 设置轴标签
    ax2.set_xlabel('Time (hour)')
    ax2.set_ylabel('Residual (km)')
    plt.legend(["X_diff","Y_diff","Z_diff"],loc="upper right",fontsize=8)
    # 调整布局与显示
    plt.tight_layout()
    plt.show()

其实，我写的过于繁琐了一些，但是目前，小编的水平有限，本来想着让AI帮我简化，使得代码更加规范简洁一点，但想着这会覆盖我之前的思考过程，所以还是保留下来，将规范化的代码另存，以便后续对比学习。

最后，呈现的结果如下所示：（）