基于夏式最小二乘的轨道控制推力拟合系数计算方法技术

本发明专利技术提供了一种基于夏式最小二乘的轨道控制推力拟合系数计算方法，综合分析历次轨道控制理论速度增量、实际速度增量、轨道控制推力系数的相关性，采用夏式最小二乘估计表征推力的拟合系数，建立一个能模仿真实系统行为的模型，用当前可测量的系统输入和输出预测系统输出的未来演变。本发明专利技术利用历史轨道控制数据作为输入，运用夏氏法广义最小二乘计算实际推力拟合系数，具有计算效率高、不需反复数据过滤等特点；通过迭代可以消去最小二乘轨迹中的偏差，因此该方法的计算相对简单，预测精度高。

Calculation method of thrust fitting coefficient for trajectory control based on Xia type least squares

The invention provides a method for calculating the fitting coefficient of orbit control thrust based on the Xia-Square Least Square, comprehensively analyzes the correlation between the theoretical velocity increment, the actual velocity increment and the orbit control thrust coefficient of previous orbit control, and adopts the Xia-Square Least Square estimation to characterize the fitting coefficient of thrust, so as to establish a simulation system. The model of unified behavior predicts the future evolution of the system output with the currently measurable system input and output. The invention uses the historical orbit control data as input, calculates the actual thrust fitting coefficient by using the generalized least squares of the Charles method, and has the characteristics of high calculation efficiency and no need for repeated data filtering; the deviation in the least squares trajectory can be eliminated by iteration, so the calculation of the method is relatively simple and the prediction accuracy is high.

【技术实现步骤摘要】
基于夏式最小二乘的轨道控制推力拟合系数计算方法
本专利技术涉及在轨航天器测控管理领域，适用于近地卫星轨道控制期间的推力系数拟合计算。
技术介绍
近地卫星由于大气阻力导致轨道衰减，其地面轨迹会逐渐偏离设定的运行轨迹，需定期通过轨道控制推力器进行轨道维持控制，使卫星地面轨迹始终满足轨迹网设计要求。卫星轨道控制推力器推力是以推进剂贮箱压力为自变量的多项式函数，其系数为装订常数，在卫星的整个寿命器保持不变。随着卫星轨控次数的增加，推进剂贮箱压力不断降低，以理论推力为基础计算的卫星速度增量与实际速度增量出现较大误差，导致轨道控制精度降低。综上所述，需要一种轨道控制推力系数拟合计算方法。传统的轨道控制速度增量计算仅仅以前一次推力器标定系数作为计算依据，并且推力系数维持不变，没有有效利用历次轨道控制数据进行优化计算，导致理论速度增量与实际速度增量误差不可测，给轨道控制策略制定带来了一定困难。
技术实现思路
为了克服现有技术的不足，本专利技术提供一种基于夏式最小二乘的轨道控制推力拟合系数计算方法，综合分析历次轨道控制理论速度增量、实际速度增量、轨道控制推力系数的相关性，采用夏式最小二乘估计表征推力的拟合系数，建立一个能模仿真实系统行为的模型，用当前可测量的系统输入和输出预测系统输出的未来演变，具有预测精度高的特点。本专利技术解决其技术问题所采用的技术方案包括以下步骤：1)计算两个贮箱在第i次轨控前的推进剂密度ρix＝1025.5-0.875×(Tix-273.15)，然后计算两个贮箱在第i次轨控前的燃料剩余质量其中，ρ0x、P0x、T0x、VTx和m0x分别表示两个贮箱加注时的气体推进剂密度、气体压强、气体绝对温度、贮箱容积和贮箱初始燃料质量，Tix和Pix分别为第i次轨控前的贮箱绝对温度和贮箱压力，下标x表示燃料储箱编号，x＝1,2；2)计算第i次轨控前的卫星质量msati＝msat0-m01+mi1-m02+mi2，其中msat0为两个贮箱加注时包含燃料在内的卫星质量；然后计算轨控推力器所需提供压力所对应的总冲IiEx＝msatiΔVx*，ΔVx*是理论计算的速度变化量；并根据比冲与压力的关系计算每次轨控的比冲IiSx＝G0+G1Pix+G2Pix2+G3Pix3，G0、G1、G2、G3为设定的比冲压力拟合系数；最后计算轨控后的贮箱压力3)根据夏氏最小二乘估计法，定义Pix、Δtix、Tix和ΔVx*作为输入量，其中Δtix是推力作用时间，ΔVix作为输出量，拟合多项式系数D0、D1、D2、D3作为待求系数；对于i次轨控的i组测量值，计算测量矩阵计算输出向量4)计算待求推力系数向量θ＝[D0,D1,D2,D3]T的最小二乘估值计算残差在第1次计算残差e时取然后利用残差构造计算噪声拟合系数的一致无偏估计计算推力系数误差Γ＝(ΦTΦ)-1ΦT；5)如果与上一次计算结果的误差小于1e-5，则进入步骤6)；否则返回步骤4)；6)计算θ的夏氏最小二乘估值得到推力系数估计值7)计算第k次轨控的轨控前推力和轨控后推力计算平均推力最后，计算得到该次轨控的预测速度变化量其中，Δtk为轨控推力持续时间，msatk为卫星当前质量。所述的比冲压力拟合系数G0、G1、G2、G3分别为1888.8、635.1、-403.6和89.61。本专利技术的有益效果是：1)本专利技术综合分析历次轨道控制理论速度增量、实际速度增量、轨道控制推力系数的相关性，采用夏式最小二乘估计表征推力的拟合系数，建立一个能模仿真实系统行为的模型，用当前可测量的系统输入和输出预测系统输出的未来演变，预测精度高，平均误差在1.5％以内。2)本专利技术有效利用历次轨道控制数据进行优化计算，改变传统的轨道控制速度增量计算仅仅以前一次推力器标定系数作为计算依据，并且推力系数维持不变的缺点。3)本专利技术对在轨卫星轨道控制策略制定具有指导意义。综上，利用历史轨道控制数据作为输入，运用夏氏法广义最小二乘计算实际推力拟合系数，具有计算效率高、不需反复数据过滤等特点；通过迭代可以消去最小二乘轨迹中的偏差，因此该方法的计算相对简单，预测精度高。附图说明图1是本专利技术的方法流程图。具体实施方式下面结合附图和实施例对本专利技术进一步说明，本专利技术包括但不仅限于下述实施例。利用卫星轨道控制期间动量与冲量之间的关系，定义输入量：轨道控制前贮箱压力P(单位：Mpa)、推力作用时间Δt(单位：s)、轨道控制前贮箱绝对温度Tx(单位：K)、理论计算的速度变化量ΔVx*(脚标“x”代表燃料贮箱编号，x＝1，2)，定义输出量：实际速度变化量ΔV(单位：m/s)，建立推力拟合系数辨识模型。通过n组测量数据，通过逐次改善偏差项精度的交替求解算法，逐步消去最小二乘轨迹中的偏差，从而得到估计系数。具体过程如下：根据动量与冲量之间的关系可得：式中是推力平均值(单位：N)、Δt是推力作用时间(单位：s)、msat是卫星当前质量(单位：kg)、ΔV是标定后的速度变化量(单位：m/s)。F是轨控前推力，Ff是轨控后推力。Di(i＝1,2,3,4)是拟合多项式系数，P是轨控前贮箱压力(单位：Mpa)，Pf是轨控后贮箱压力，脚标“x”代表燃料贮箱编号(x＝1，2)，推力器工作时贮箱1和贮箱2根据压强变化互相切换。其中，贮箱1、2当前推进剂密度(单位：kg/cm3)：ρx＝1025.5-0.875×(Tx-273.15)(5)燃料剩余质量(单位：kg)：比冲与压力的关系(单位：N·s/kg)：IS＝G0+G1Px+G2Px2+G3Px3(7)本次机动推力器所需提供压力所对应的总冲(单位：N·s)：IEx＝msatΔVx*(8)卫星当前质量(单位：kg)：msat＝msat0-m01+m1-m02+m2(9)其中，VTx为贮箱容积(单位：m3)，Tx为轨控前贮箱绝对温度(单位：K)，T0x、P0x为加注时气体绝对温度和压力(单位：K、Mpa)，m01、m02为贮箱1、2初始燃料质量(单位：kg)，ρ0x初始燃料密度(单位：kg/cm3)，msat0初始卫星质量(单位：kg)，ΔVx*是理论计算的速度变化量，脚标“x”代表燃料贮箱编号(x＝1，2)。综上，根据夏氏最小二乘估计法，定义Px、Δtx、Tx、ΔVx*是输入量，ΔVx是输出量，D0、D1、D2、D3是待求系数。假设有n组测量值，则有：式(10)中i代表第i组测量值。将式(10)展开后可得到：根据式(11)，可得到测量矩阵Φ、待求参数向量θ、输出向量y：当考虑系统噪声ξ影响时，定义噪声拟合系数为f＝[f1f2f3f4]T，系统噪声ξ的拟合公式为当考虑系统噪声影响时，系统的差分方程可表示为：其中，f为噪声拟合系数向量，ε为均值为0的白噪声向量，为不相关的随机序列，故可用最小二乘法得到系数f的一致无偏估计。若有n组测量数据，考虑到计算效率，采用夏氏改良法，则有其中，e为残差，Ω的具体表达式如下综上，参照迭代计算步骤即可求得估计系数。选取某低轨卫星共n＝45次轨控进行实例计算，该星有两个推进剂贮箱，每次轨控根据贮箱压力差进行切换使用。1)数据初始化。数据输入，包括两个贮箱加注时气体推进剂密度ρ0x(单位：kg/cm3)、气体压强P0x(单位：Mpa)、气体绝对温度T0x(单位：K)，贮箱容积VTx(单位：本文档来自技高网...

【技术保护点】
1.一种基于夏式最小二乘的轨道控制推力拟合系数计算方法，其特征在于包括下述步骤：1)计算两个贮箱在第i次轨控前的推进剂密度ρix＝1025.5‑0.875×(Tix‑273.15)，然后计算两个贮箱在第i次轨控前的燃料剩余质量

【技术特征摘要】
1.一种基于夏式最小二乘的轨道控制推力拟合系数计算方法，其特征在于包括下述步骤：1)计算两个贮箱在第i次轨控前的推进剂密度ρix＝1025.5-0.875×(Tix-273.15)，然后计算两个贮箱在第i次轨控前的燃料剩余质量其中，ρ0x、P0x、T0x、VTx和m0x分别表示两个贮箱加注时的气体推进剂密度、气体压强、气体绝对温度、贮箱容积和贮箱初始燃料质量，Tix和Pix分别为第i次轨控前的贮箱绝对温度和贮箱压力，下标x表示燃料储箱编号，x＝1,2；2)计算第i次轨控前的卫星质量msati＝msat0-m01+mi1-m02+mi2，其中msat0为两个贮箱加注时包含燃料在内的卫星质量；然后计算轨控推力器所需提供压力所对应的总冲IiEx＝msatiΔVx*，ΔVx*是理论计算的速度变化量；并根据比冲与压力的关系计算每次轨控的比冲IiSx＝G0+G1Pix+G2Pi...

【专利技术属性】
技术研发人员：卞燕山，王西京，陈军，崔鹏，张莹，李伟，张雷，
申请(专利权)人：中国人民解放军六三七八九部队，
类型：发明
国别省市：陕西,61