寻求模型火箭发动机推力数据插值与多项式拟合方案

架构老林

2026-6-30

问题描述

我正在开发一个处理模型火箭发动机静态测试TXT数据的程序，已经实现了给数据添加时间戳、生成推力曲线的功能，但在处理大量数据时，插值和提取最优拟合多项式遇到了麻烦：用SciPy文档里的方法会因为数据量过大报错，不知道怎么用NumPy实现相关功能。

已完成代码

# 初始目标是给每一行数据添加对应时间戳，用于后续绘图
input_file_path = 'data.txt'
output_file_path = 'output_file.txt'

with open(input_file_path, 'r') as input_file:
    lines = input_file.readlines()

# 处理每一行并添加时间戳
for i, line in enumerate(lines):
    time_value = i
    lines[i] = f'{time_value} {line}'

# 将新内容写入另一个TXT文件，用于后续计算
with open(output_file_path, 'w') as output_file:
    output_file.writelines(lines)

# 用Pandas读取TXT文件并提取数据
import pandas as pd
import matplotlib.pyplot as plt

data = pd.read_csv('output_file.txt', sep=' ', header=None)
data = pd.DataFrame(data)

# 用Matplotlib绘制推力曲线
x = data[0]
y = data[1]
fig = plt.figure(figsize=(15,10))
ax1 = fig.add_subplot(111)
ax1.set_title("Thrust Curve")
ax1.set_xlabel('Time (seconds)')
ax1.set_ylabel('Thrust (Newtons)')
plt.plot(x, y, color='blue', linestyle='-')
plt.show()

测试TXT数据示例

解决方案

针对大数据量的插值和多项式拟合，优先使用NumPy的矢量化运算（避免循环，效率远高于SciPy部分方法），以下是具体实现：

1. 优化数据读取与时间戳添加

原代码通过读写中间文件处理时间戳，大数据量下效率低，直接用Pandas/NumPy一步完成：

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# 直接读取数据并生成时间戳（无需中间文件）
input_file_path = 'data.txt'
# 读取数据，假设每行是单个推力值
thrust_data = pd.read_csv(input_file_path, header=None, names=['thrust'])
# 生成时间戳（从0开始，步长1秒）
thrust_data['time'] = np.arange(len(thrust_data))
# 提取NumPy数组用于后续运算
x = thrust_data['time'].values
y = thrust_data['thrust'].values

2. NumPy实现高效插值

针对大数据量，用np.interp实现线性插值（矢量化操作，无性能瓶颈），如果需要高阶插值，用numpy.polynomial.Polynomial的拟合插值：

线性插值（适合平滑补全数据）

# 生成需要插值的目标时间点（比如原数据是1秒步长，现在要0.1秒步高）
target_x = np.linspace(x.min(), x.max(), num=len(x)*10)
# 执行线性插值
interpolated_y = np.interp(target_x, x, y)

# 绘制插值后的曲线
plt.figure(figsize=(15,10))
plt.plot(x, y, 'b-', label='原始数据')
plt.plot(target_x, interpolated_y, 'r--', label='插值后数据')
plt.title("推力曲线（含插值）")
plt.xlabel('Time (seconds)')
plt.ylabel('Thrust (Newtons)')
plt.legend()
plt.show()

高阶多项式插值（适合拟合趋势）

如果需要高阶插值，用numpy.polynomial.Polynomial.fit，比SciPy的方法更适合大数据量：

# 拟合5阶多项式作为插值模型（可根据需求调整阶数）
poly_model = np.polynomial.Polynomial.fit(x, y, deg=5)
# 生成插值结果
interpolated_y_high = poly_model(target_x)

# 绘制对比曲线
plt.figure(figsize=(15,10))
plt.plot(x, y, 'b-', label='原始数据')
plt.plot(target_x, interpolated_y_high, 'g--', label='5阶多项式插值')
plt.title("推力曲线（高阶多项式插值）")
plt.xlabel('Time (seconds)')
plt.ylabel('Thrust (Newtons)')
plt.legend()
plt.show()

3. 自动提取最优拟合多项式

要找到最优拟合阶数，通过计算不同阶数下的均方误差（MSE），选择MSE最小且阶数合理的模型：

def find_optimal_poly(x, y, max_deg=10):
    """遍历阶数，找到最优拟合多项式"""
    mse_list = []
    models = []
    for deg in range(1, max_deg+1):
        model = np.polynomial.Polynomial.fit(x, y, deg=deg)
        y_pred = model(x)
        mse = np.mean((y - y_pred)**2)
        mse_list.append(mse)
        models.append(model)
    
    # 找到MSE最小的阶数
    optimal_idx = np.argmin(mse_list)
    optimal_deg = optimal_idx + 1
    optimal_model = models[optimal_idx]
    
    print(f"最优拟合阶数：{optimal_deg}，对应MSE：{mse_list[optimal_idx]:.4f}")
    return optimal_model, optimal_deg

# 获取最优模型
optimal_poly, optimal_deg = find_optimal_poly(x, y, max_deg=10)
# 生成拟合曲线数据
fit_y = optimal_poly(x)

# 绘制原始数据与最优拟合曲线
plt.figure(figsize=(15,10))
plt.plot(x, y, 'b-', label='原始数据')
plt.plot(x, fit_y, 'm--', label=f'最优拟合（{optimal_deg}阶多项式）')
plt.title("推力曲线（最优多项式拟合）")
plt.xlabel('Time (seconds)')
plt.ylabel('Thrust (Newtons)')
plt.legend()
plt.show()

关键说明

NumPy的矢量化运算完全基于C实现，处理百万级数据也不会出现性能问题，避免了SciPy部分方法在大数据量下的内存/计算瓶颈。
最优阶数选择时，建议限制max_deg在10以内，避免过拟合（火箭发动机推力曲线通常是低阶趋势）。

内容的提问来源于stack exchange，提问作者Ilzy Lima Vieira

寻求模型火箭发动机推力数据插值与多项式拟合方案

已完成代码 #

测试TXT数据示例 #

1. 优化数据读取与时间戳添加 #

2. NumPy实现高效插值 #

线性插值（适合平滑补全数据） #

高阶多项式插值（适合拟合趋势） #

3. 自动提取最优拟合多项式 #

关键说明 #