CN116954075A - 非合作连续推力航天器推力参数辨识方法、系统和设备 - Google Patents

非合作连续推力航天器推力参数辨识方法、系统和设备 Download PDF

Info

Publication number
CN116954075A
CN116954075A CN202310897508.5A CN202310897508A CN116954075A CN 116954075 A CN116954075 A CN 116954075A CN 202310897508 A CN202310897508 A CN 202310897508A CN 116954075 A CN116954075 A CN 116954075A
Authority
CN
China
Prior art keywords
thrust
inclination angle
track
spacecraft
theoretical
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202310897508.5A
Other languages
English (en)
Other versions
CN116954075B (zh
Inventor
赵帅龙
李智
陶雪峰
张雅声
程文华
刁华飞
王训
王雅如
汪夏
薛锦妍
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Peoples Liberation Army Strategic Support Force Aerospace Engineering University
Original Assignee
Peoples Liberation Army Strategic Support Force Aerospace Engineering University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Peoples Liberation Army Strategic Support Force Aerospace Engineering University filed Critical Peoples Liberation Army Strategic Support Force Aerospace Engineering University
Priority to CN202310897508.5A priority Critical patent/CN116954075B/zh
Publication of CN116954075A publication Critical patent/CN116954075A/zh
Application granted granted Critical
Publication of CN116954075B publication Critical patent/CN116954075B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • B—PERFORMING OPERATIONS; TRANSPORTING
    • B64—AIRCRAFT; AVIATION; COSMONAUTICS
    • B64G—COSMONAUTICS; VEHICLES OR EQUIPMENT THEREFOR
    • B64G1/00—Cosmonautic vehicles
    • B64G1/22—Parts of, or equipment specially adapted for fitting in or to, cosmonautic vehicles
    • B64G1/24—Guiding or controlling apparatus, e.g. for attitude control
    • B64G1/242—Orbits and trajectories
    • G—PHYSICS
    • G05—CONTROLLING; REGULATING
    • G05B—CONTROL OR REGULATING SYSTEMS IN GENERAL; FUNCTIONAL ELEMENTS OF SUCH SYSTEMS; MONITORING OR TESTING ARRANGEMENTS FOR SUCH SYSTEMS OR ELEMENTS
    • G05B13/00—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion
    • G05B13/02—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric
    • G05B13/04—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric involving the use of models or simulators
    • G05B13/042—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric involving the use of models or simulators in which a parameter or coefficient is automatically adjusted to optimise the performance
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G06F17/13—Differential equations
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Computational Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Theoretical Computer Science (AREA)
  • Software Systems (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Remote Sensing (AREA)
  • General Engineering & Computer Science (AREA)
  • Artificial Intelligence (AREA)
  • Automation & Control Theory (AREA)
  • Medical Informatics (AREA)
  • Evolutionary Computation (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Operations Research (AREA)
  • Chemical & Material Sciences (AREA)
  • Health & Medical Sciences (AREA)
  • Aviation & Aerospace Engineering (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Combustion & Propulsion (AREA)
  • Computing Systems (AREA)
  • Navigation (AREA)

Abstract

本发明公开了一种非合作连续推力航天器推力参数辨识方法、系统和设备，解决了无法对非合作连续推力航天器推力参数辨识的技术问题；属于卫星轨道控制领域；包括：将接收到的非合作连续推力航天器的雷达数据转换为非初始时刻的轨道面数据，得到真实轨道倾角模型，并输出真实轨道倾角集；将理论轨道倾角解析表达，得到理论轨道倾角模型，输出理论轨道倾角集；采用最小二乘的思想，设置轨道倾角误差集的欧几里得范数作为评价指标；对欧几里得范数进行极值寻优，得到推力参数。本发明求解速度快，精度高，通过对非合作连续推力航天器的推力参数进行求解，能够研判机动策略，减少其碰撞风险。

Description

非合作连续推力航天器推力参数辨识方法、系统和设备
技术领域
本发明属于高端装备制造领域中的卫星轨道控制技术领域，涉及一种非合作连续推力航天器推力参数辨识方法、系统、设备和存储介质。
背景技术
连续推力航天器是指航天器的推进方式为连续推力，一般情况下认为各方向上的推力加速度相对连续推力航天器保持不变，推进器持续向连续推力航天器施加推力。
对于传统只搭载脉冲推力载荷的航天器轨道机动问题，可通过分析其轨道参数改变量进行数值计算，从而得知其脉冲推力加速度大小和方向，对其机动策略和目的进行研判。但对于搭载连续推力载荷的航天器而言，研究简单的轨道参数变化量不足以解析出其推力参数。
非合作连续推力航天器泛指不能提供有效合作信息的连续推力航天器，一般是指没有装备通信应答机或者其他传感器的航天器，其他航天器无法采用电子讯问及发射信号等手段实现对所述非合作连续推力航天器空间目标进行判别或定位。
非合作连续推力航天器具有以下特点：没有安装特征块和合作标志器；没有安装特殊设计的对接接口；不能主动传送其姿态信息。非合作连续推力航天器通常包括己方未配置合作接口的卫星、安装合作接口但发生故障或燃料耗尽的卫星、失效卫星的空间碎片或其他无合作信息的航天器等。
随着星链卫星等非合作连续推力航天器在近地空间中的快速高密度部署，非合作可机动目标编目、航天器间碰撞预警等太空态势感知任务的难度日益增加。
发明人在研究的过程中发现，传统的计算方法求解只搭载脉冲推力载荷的航天器的轨道机动问题是在摄动条件下求解圆锥曲线问题，即航天器在空间中的轨道可近似为圆锥曲线；但用此方法对非合作连续推力航天器施加连续推力，相当于给非合作连续推力航天器在空间中增加二或三个方向上的加速度参数，导致非合作连续推力航天器在空间中的轨道无法用普通的圆锥曲线描述，增加的参数以及无法用圆锥曲线描述轨道给问题的求解增加了难度，因此，传统的计算方法和航天器动力学方程不再适用于非合作连续推力航天器的推力参数求解。
此类航天器处于连续轨道机动状态中，由于该类航天器在未知其推进器各方向推力参数的前提下，无法对其推力参数辨识，难以评估其轨道误差，也无法研判其机动策略，故与其他航天器存在碰撞风险。
发明内容
针对非合作连续推力航天器无法预知推力参数导致无法对其进行轨道预推以及无法分析其机动策略的技术问题，本发明提供了一种非合作连续推力航天器推力参数辨识方法、系统和设备，本发明通过将接收到的非合作连续推力航天器的雷达数据转换得到轨道倾角模型，输出真实轨道倾角集以及对理论轨道倾角积分解析表达得到理论轨道倾角模型并输出理论轨道倾角集；采用最小二乘的思想，设置理论轨道倾角集与真实轨道倾角集之间的轨道倾角误差集的欧几里得范数作为评价指标；对轨道倾角误差集的欧几里得范数进行极值寻优，得到非合作连续推力航天器的推力参数；本发明针对传统计算脉冲推力加速度大小和方向的方法不再适用星链卫星等非合作连续推力航天器，寻求新的方法对非合作连续推力航天器空间目标进行轨道确定，根据非合作连续推力航天器极短弧信息中轨道倾角的变化情况，通过多次数学变换，采用最小二乘的数学思想，将非合作连续推力航天器的推力反演问题转化为最小二乘优化的数学问题，反演出本段机动过程中航天器所受各方向的推力参数，研判其机动策略，减小非合作连续推力航天器与其他航天器存在的碰撞风险。
本发明的目的通过以下技术方案来具体实现：
本发明公开了一种非合作连续推力航天器推力参数辨识方法，该方法包括：
步骤一、将接收到的非合作连续推力航天器的雷达数据转换为非初始时刻的轨道面数据，基于非初始时刻轨道面与赤道平面之间夹角得到真实轨道倾角模型；将与非初始时刻时间集对应的位置坐标集输入真实轨道倾角模型中，得到真实轨道倾角集；
步骤二、将非合作连续推力航天器的理论轨道倾角解析表达，得到理论轨道倾角模型，将非初始时刻时间集输入理论轨道倾角模型中，得到理论轨道倾角集；
步骤三、采用最小二乘的思想，设置理论轨道倾角集与真实轨道倾角集之间的轨道倾角误差集的欧几里得范数作为评价指标；
步骤四、对作为评价指标的轨道倾角误差集的欧几里得范数进行极值寻优得到最小值，得到非合作连续推力航天器各方向上的推力参数。
步骤一中，真实轨道倾角模型为：
其中，ii为非初始时刻的真实轨道倾角，为非初始时刻轨道面与赤道平面之间夹角；e表示赤道平面的法向量；hi表示非初始时刻轨道面的法向量；|e|为赤道平面法向量的模；|hi|为非初始时刻轨道面法向量的模。
步骤二中，理论轨道倾角模型为：
其中，i为非初始时刻的理论轨道倾角；i0为初始时刻已知的轨道倾角；和/>为推力参数，/>为单位质量所受的切向力，即切向推力加速度；/>为单位质量所受的法向力，即法向推力加速度；μ为地球引力常数；r0为轨道半径；t为推力施加时间。
步骤三中，采用最小二乘的思想，设置的理论轨道倾角集与真实轨道倾角集之间的轨道倾角误差集的欧几里得范数为：
其中，Fs为欧几里得范数；t为推力施加时间，t＝(1,2,3,…T)，T为推力施加总时长，T>0；errt为推力施加时间对应的轨道倾角误差；errT为推力施加总时长对应的轨道倾角误差。
步骤三中，理论轨道倾角集与真实轨道倾角集之间的轨道倾角误差集为：
ERR＝Is-I＝(err1,err2,err3,…errT)；
其中，ERR为轨道倾角误差集；Is为真实轨道倾角集；I为理论轨道倾角集；errT为推力施加总时长对应的轨道倾角误差；T为推力施加总时长，T>0。
步骤四中，对作为评价指标的轨道倾角误差集的欧几里得范数进行极值寻优的方法，包括：通过LM算法、遗传算法、粒子群算法或蚁群算法对欧几里得范数进行极值寻优，得到欧几里得范数的最小值。
优选的，本发明利用遗传算法对轨道倾角误差集的欧几里得范数进行极值寻优，步骤包括：
步骤1：设置所求轨道倾角误差集的欧几里得范数的待求参数切向推力加速度和法向推力加速度的取值边界；
步骤2：设置进化迭代计数器的初始值及进化代数的最大值，并随机生成待求参数切向推力加速度和法向推力加速度的多个个体组成初始群体；
步骤3：计算初始群体中每个个体的适应度值，并根据适应度值对所有个体进行排序；
步骤4：排序后的初始群体经过遗传操作中的选择运算、交叉运算和变异运算之后得到下一代种群群体；进化迭代计数器计时，初始值增加一次；
步骤5：若进化迭代计数器的数值不大于最大进化代数，则进化迭代计数器继续计时，循环执行步骤3和步骤4；若进化迭代计数器的数值大于最大进化代数，则进化迭代计数器停止计时，此进化过程中所得到的具有最大适应度值的个体作为最优解输出，输出的最优解为非合作连续推力航天器的切向推力加速度和法向推力加速度。
本发明还提供了一种非合作连续推力航天器推力参数辨识系统，包括：
真实轨道倾角集模块，用于将接收到的非合作连续推力航天器的雷达数据转换为非初始时刻的轨道面数据，基于非初始时刻轨道面与赤道平面之间夹角得到真实轨道倾角模型；将与非初始时刻时间集对应的位置坐标集输入真实轨道倾角模型中，得到真实轨道倾角集；
理论轨道倾角集模块，用于对非合作连续推力航天器的理论轨道倾角积分解析表达，得到理论轨道倾角模型，将非初始时刻时间集输入理论轨道倾角模型中，得到理论轨道倾角集；
评价指标模块，用于采用最小二乘的思想，设置理论轨道倾角集与真实轨道倾角集之间的轨道倾角误差集的欧几里得范数作为评价指标；
推力参数求解模块，用于对作为评价指标的轨道倾角误差集的欧几里得范数进行极值寻优得到最小值，得到非合作连续推力航天器各方向上的推力参数。
本发明还提供了一种非合作连续推力航天器推力参数辨识设备，包括存储器和处理器，所述存储器存储有计算机程序，所述处理器执行所述计算机程序时实现一种非合作连续推力航天器推力参数辨识方法的步骤。
本发明还提供了一种计算机可读存储介质，其存储有计算机程序，所述计算机程序被处理器执行时实现一种非合作连续推力航天器推力参数辨识方法的步骤。
本发明的有益效果是：
本发明通过将接收到的非合作连续推力航天器的雷达数据转换为非初始时刻的轨道面数据，基于非初始时刻轨道面与赤道平面之间夹角得到轨道倾角模型；将与非初始时刻时间集对应的位置坐标集输入轨道倾角模型中，得到真实轨道倾角集；利用连续推力作用下，航天器的轨道升交点赤经受J2项非球形摄动发生进动这一特点，空间转化后将雷达观测数据集转换为轨道倾角模型；
对轨道倾角模型积分解析表达，得到轨道倾角变化模型，将非初始时刻时间集输入轨道倾角变化模型中，得到理论轨道倾角集；解决了连续推力圆轨道的轨道倾角无法解析表达的问题；
采用最小二乘的思想，设置理论轨道倾角集与真实轨道倾角集之间的轨道倾角误差集的欧几里得范数作为评价指标；将非合作连续推力航天器的推力反演问题转化为最小二乘优化的数学问题。
利用遗传算法对轨道倾角误差集的欧几里得范数进行极值寻优，得到非合作连续推力航天器的切向推力加速度和法向推力加速度。本发明采用遗传算法的优势：
由于本发明针对不能提供有效合作信息的非合作连续推力航天器的推力参数辨识，发明人经过试验，遗传算法基本上不用搜索空间的知识或其它辅助信息，而仅用适应度函数值来评估个体，在此基础上进行遗传操作。适应度函数不仅不受连续可微的约束，而且其定义域可以结合本发明针对的非合作连续推力航天器进行设定。遗传算法的搜索从轨道倾角误差集的欧几里得范数群体出发，利用内在并行性进行分布式计算，同时处理群体中的多个个体，即对搜索空间中的多个解进行评价，减少了陷入局部最优解的风险，也加快了求解速度，提高了求解精度。
附图说明
下面根据附图和实施例对本发明作进一步详细说明。
图1是天东北坐标系下的空间几何关系示意图。
图2是天东北坐标系与地心坐标系的空间几何关系示意图。
图3是地心坐标系与地心惯性坐标系的空间几何关系示意图。
图4是相邻两时刻轨道变化示意图。
图5是地基雷达可见性几何示意图。
图6是地基雷达可见性时间窗口示意图。
图7是切向推力加速度求解误差情况示意图。
图8是法向推力加速度求解误差情况示意图。
具体实施方式
实施例一
本发明实施例一提供了一种非合作连续推力航天器推力参数辨识方法，该方法包括：
步骤一、将接收到的非合作连续推力航天器的雷达数据转换为非初始时刻的轨道面数据，基于非初始时刻轨道面与赤道平面之间夹角得到真实轨道倾角模型；将与非初始时刻时间集对应的位置坐标集输入真实轨道倾角模型中，得到真实轨道倾角集。
步骤具体包括：
(一)地球观测站将接收到的非合作连续推力航天器的雷达数据坐标，由天东北坐标系下的坐标转换成地心惯性坐标系下的坐标：
地球观测站接收到的非合作连续推力航天器的雷达观测数据为(ρ,Az,h)，其中，ρ为观测距离，Az为方位角，h为俯仰角；现有地基观测手段以雷达和光学望远镜为主，地球观测站在地球上的位置固定不变；如图1所示，P为航天器位置，P’为其在测站平面Ls上的投影位置，天东北坐标系原点S位于地球观测站，各坐标轴的定义如下：
(1)SZ(XS)轴：沿过观测站的铅垂线方向，指向天顶。
(2)SE(YS)轴：位于过观测站的地平面内，指向正东。
(3)SN(ZS)轴：位于过观测站的地平面内，指向正北。
由于该坐标系是固连在地球上的，所以会随地球一起自转，是一个动坐标系，用字符S表示。利用该坐标系可以计算出卫星相对于观测站的位置。将雷达数据转化为天东北坐标系S下的坐标Ps：
Ps＝(ρsinh,ρcoshsinAz,ρcoshcosAz) (1)
由图2所示，地心坐标系原点位于地心OE,各坐标轴的定义如下：
(1)OEXE轴：位于在赤道平面内，由地心指向某时刻t0的起始子午线(通常取格林尼治天文台所在的子午线)，显然该坐标系是随着地球的自转而转动的。
(2)OEZE轴：垂直于赤道平面，与地球自转轴重合，指向北极。
(3)OEYE轴：位于赤道平面内，其方向满足右手直角坐标系准则。
由坐标系的定义可知，该坐标系的OEXE轴和OEYE轴随着地球的自转而转动，是一个动坐标系，用字符E表示,其中，为测站纬度，λs为测站经度。
天东北坐标系S与地心坐标系E之间的转换关系为：
如图3所示，地心惯性坐标系原点位于地心OE，各坐标轴的定义如下：
(1)OEXI轴：位于赤道平面内，由地心指向平春分点。由于春分点是随着时间而变化的，所以此处的平春分点规定为2000年1月1日12：00的平春分点。
(2)OEYI轴：垂直于赤道平面，与地球自转轴重合，指向北极。
(3)OEZI轴：位于赤道平面内，其方向满足右手直角坐标系准则。
由坐标系的定义可知，该坐标系的各坐标轴在惯性空间保持方向不变，是一个惯性坐标系，用字符I表示。
地心坐标系E与地心惯性坐标系I的坐标原点和OEZI均重合，而差别在于OEXI指向平春分点，而OEXE指向所讨论时刻格林尼治天文台所在子午线与赤道的交点。地心坐标系E与地心惯性坐标系I的夹角可以通过天文年历表查算得到，记为ΩG。由于OEXI轴是固定的，而OEXE轴是随着地球转动的，所以角ΩG随所讨论的时刻不同而不同。因此，不难解出地心坐标系E与地心惯性坐标系I之间转换矩阵关系为：
EI＝M3[-ΩG] (3)
综上所述，根据雷达观测数据，可得非合作连续推力航天器在地心惯性坐标系中的坐标为：
其中，三个坐标系之间的初等转换矩阵为第一初等转换矩阵M1、第二初等转换矩阵M2和第三初等转换矩阵M3：
式中，θ为转动的角度。
(二)利用非合作连续推力航天器的雷达数据在地心惯性坐标系下的坐标得到非合作连续推力航天器在非初始时刻的位置坐标及非初始轨道面参数；基于非初始时刻轨道面与赤道平面之间夹角得到真实轨道倾角模型的步骤包括：
如图4所示，在空间中建立空间直角坐标系O-xyz，将地球看作质点O，赤道平面E在xOy平面上，赤道平面可表示为z＝0。
近地空间中，航天器沿其轨道从南向北运动时与赤道平面的交点称为升交点，从北向南运动时与赤道平面的交点称为降交点，轨道面与赤道平面的交线(即升交点与降交点之间连线所在的直线)称为“交点线”。航天器连续长期机动中升交点赤经不受推力影响，但在J2项摄动影响下，交点线会发生进动。
其中，为t时刻的升交点赤经，J2为摄动系数，aE为地球半径，/>为切向推力加速度，v0为初始速度，t为时刻，Ω0为初始升交点赤经。意味着初始轨道面L0和平面Li与赤道平面E分别相交形成的交点线l0与li存在角度差，为ΔΩi。交点线l0可表示为
ti时刻的升交点赤经为Ωi＝Ω0+ΔΩi，交点线li可表示为
点P0为连续推力航天器在初始时刻t0的位置，点Pi表示航天器在ti时刻的位置；L0为初始轨道面，Li为ti时刻的轨道面，点P0:(x0,y0,z0)和点Pi:(x0,y0,z0)的表达式如下
直线OP0表达式为
OP0:(x-x0)＝(y-y0)＝(z-z0) (10)
直线OPi表达式为
OPi:(x-xi)＝(y-yi)＝(z-zi) (11)
交点线l0和直线OP0都在初始轨道面L0上，且在升交点角距≠0时不平行，面L0法向量h0为
h0＝l0×OP0(12)
e＝(0,0,1),为赤道平面E法向量。面L0与赤道平面E所成夹角为初始轨道倾角i0，
同上述推导，ti时刻的轨道倾角ii为
式中，hi为面Li法向量
hi＝li×OPi(15)
(三)将与非初始时刻时间集对应的位置坐标集输入真实轨道倾角模型中，得到真实轨道倾角集Is。
步骤二、对非合作连续推力航天器的理论轨道倾角积分解析表达，得到理论轨道倾角模型，将非初始时刻时间集输入理论轨道倾角模型中，得到理论轨道倾角集。
具体步骤包括：
(一)对非合作连续推力航天器的理论轨道倾角积分解析表达，得到理论轨道倾角模型：
连续推力推进器比冲高，短时间内其质量消耗可忽略不计，故认为其质量短时间内不发生改变。在实际工程任务中，航天器始终保持近圆轨道，故推力加速度矢量径向分量通常为0，此时控制方程有：
式中，r为圆轨道半径，i为轨道倾角，Ω为轨道升交点赤经，u为升交点角距，为切向推力加速度，/>为法向推力加速度，；v为航天器速度，n为航天器角速度。上式摄动运动方程在计算过程中不存在奇异点。式(14)可简化为：
圆轨道的活力公式可简化为
对上式微分，有
则有
下面对轨道半径进行变形，
对上式进行积分，
式中，Cr为积分常数，当t＝0时，故，
Cr＝0(21)
有，
对轨道倾角进行推导
β＝±π/2时变号，倾角在一个周期内的改变量为
一个轨道周期Δt＝2πa/v，有
对上式进行积分
式中，Ci是积分常数，当t＝0时，故
有，
(二)将非初始时刻时间集输入理论轨道倾角模型中，得到理论轨道倾角集；
将t＝(t1,t2,t3,…tT)代入式(28)可得理论轨道倾角集I。
步骤三、采用最小二乘的思想，设置理论轨道倾角集与真实轨道倾角集之间的轨道倾角误差集的欧几里得范数作为评价指标；
具体步骤包括：
定义轨道倾角误差集ERR如下：
ERR＝Is-I＝(err1,err2,err3,…errT) (29)
采用最小二乘的思想，设置轨道倾角误差集ER的欧几里得范数Fs作为评价指标，用来描述理论轨道倾角集I与真实轨道倾角集Is之间的误差，从而评价所求切向推力加速度和法向推力加速度的求解精度。
观察式(28)、(29)和(30)可知，理论轨道倾角集I与参数r0,i0,Ω0,t,J2,aE,μ有关，其中，
r0,i0,Ω0,J2,aE,μ作为已知参数，故代入t＝(t1,t2,t3,…tT)后，理论轨道倾角集I是关于切向推力加速度和法向推力加速度/>的函数。所以轨道倾角误差集的欧几里得范数Fs是关于切向推力加速度/>和法向推力加速度/>的复杂二元函数。
步骤三通过多次数学变换，采用最小二乘的数学思想，将非合作连续推力航天器的推力反演问题转化为二元函数求极值的数学问题。
步骤四、对作为评价指标的轨道倾角误差集的欧几里得范数进行极值寻优得到最小值，将最小值作为评价指标得到非合作连续推力航天器各方向上的推力参数。
Fs的表达式复杂、数据量大，难以通过传统求偏导的方式对函数求极值，可通过对二元函数求极值的方法，包括但不限于：LM算法、遗传算法、粒子群算法或蚁群算法对函数求极值。
本发明实施例使用遗传算法进行极值寻优，利用遗传算法对轨道倾角误差集的欧几里得范数进行极值寻优的步骤包括：
步骤1：设置所求轨道倾角误差集的欧几里得范数的待求参数切向推力加速度和法向推力加速度的取值边界；
步骤2：设置进化迭代计数器的初始值及进化代数的最大值，并随机生成待求参数切向推力加速度和法向推力加速度的多个个体组成初始群体；
设置进化迭代计数器k＝0，设置最大进化代数K，随机生成N个体作为初始群体P(0)。
步骤3：计算初始群体中每个个体的适应度值，并根据适应度值对所有个体进行排序；
例如：将目标函数的倒数作为适应度函数，每个染色体的适应度值就是它的目标函数值的倒数：
目标函数值越小，其适应度函数值越大，表示其更优，被选择作为下一步中的父代的概率更大。
还可以通过对适应度函数进行尺度变换的方法进行计算，包括但不限于：线性变换、幂函数变换或指数变换的方法。
步骤4：初始群体经过遗传操作中的选择运算、交叉运算和变异运算之后得到下一代种群群体；进化迭代计数器计时，初始值增加一次；
将选择算子作用于群体，由个体的适应度大小来决定选择概率，选择一些优良个体遗传到下一代群体，例如：
式中,pj表示个体被选择遗传的概率；
将交叉算子作用于群体，对选中的成对个体，交换它们之间的部分染色体，产生新的染色体，例如：采用最简单的单点交叉，交叉点随机产生，交叉率设为0.6；
将变异算子作用于群体，对选中的个体，改变某一个或一些基因值为其他的等位基因；例如：采取基因位反转变异法，即0变为1，1变为0。要进行变异的基因位的选取也是随机的。变异率为0.2。
群体经过选择、交叉和变异运算之后得到下一代群体。
步骤5：若进化迭代计数器的数值不大于最大进化代数，则进化迭代计数器继续计时，循环执行步骤3和步骤4；若进化迭代计数器的数值大于最大进化代数，则进化迭代计数器停止计时，此进化过程中所得到的具有最大适应度值的个体作为最优解输出，输出的最优解为非合作连续推力航天器的切向推力加速度和法向推力加速度。
终止条件判断：若k≤K，则k＝k+1，转到步骤3；若k＞K，则此进化过程中所得到的具有最大适应度的个体作为最优解输出，终止计算。
下面提供一具体例子对本发明提供的方法详细说明：
为了对连续推力作用下的非合作连续推力航天器的推力参数辨识进行分析，需要在实际场景中进行全模型仿真，从而验证本文参数求解方法的精度，仿真条件设置如表1所示。
表1仿真参数设置
在各种摄动力作用下，非合作连续推力航天器沿其轨道进行连续机动，考虑到雷达观测实际，限定非合作连续推力航天器出测站地平5°之后与入测站地平线5°之前为雷达可观测范围，如图5所示，为地基雷达可见性几何示意图，灰色区域为可观测区域；如图6所示，为地基雷达可见性时间窗口示意图，灰色区段为可观测窗口期，低轨卫星轨道周期短，可观测窗口较多，但每次窗口的观测时间短，一般为10min左右。仿真中设置采样步长为60s，每次采样时间根据观测窗口时长而定，对2020.04.02 03:22:22 UTC至2020.04.0403:22:22 UTC时间段内所有可观测窗口进行观测采样，时间跨度48h，共得到307组观测数据，记为(ti,ρi,Azi,Eli)，i＝1,2,…,307。
本文使用LM算法进行参数估计。此方法综合最速下降法和线性化方法，能够较快地找到轨道参数和推力参数最优值。采用新建的连续推力轨道参数集描述求解结果如下。
表2求解结果(UTC:2020.04.02 -2020.04.04)
轨道参数 实际值 求解结果 误差
切向推力加速度(m/s2) 0.00433 0.004326 0.000004
法向推力加速度(m/s2) 0.00250 0.00245 0.00005
利用上述方法进行推力参数求解和轨道确定，共进行了1000次蒙特卡洛仿真，航天器切向推力加速度求解结果的相对误差在4×10-6m/s2左右，法向推力加速度求解结果的相对误差在5×10-5m/s2左右。图7和图8显示了切向推力加速度和法向推力加速度求解误差情况。
从上述结果可以看出航天器切向和法向推力加速度参数的估计精度很高，能够应用于工程实际，切向推力加速度的求解精度更高是因为其只涉及轨道半径的变化，求解过程简单，存在更少的数据观测误差和系统误差。整个求解过程能够解析表达，不涉及复杂的数学计算。
综上，本发明实施例在求解非合作连续推力航天器推力参数时误差小，精度高，求解速度快，不依赖大量观测数据，利用短弧观测数据即可进行求解推力参数，在航天器推力参数改变的情况下能够实时更新其参数信息，提高航天器运动控制精度，更好地进行轨道预报。
实施例二
为了执行上述实施例一对应的方法，以实现相应的功能和技术效果，本发明实施例二提供了一种非合作连续推力航天器推力参数辨识系统，包括：
真实轨道倾角集模块，用于将接收到的非合作连续推力航天器的雷达数据转换为非初始时刻的轨道面数据，基于非初始时刻轨道面与赤道平面之间夹角得到真实轨道倾角模型；将与非初始时刻时间集对应的位置坐标集输入真实轨道倾角模型中，得到真实轨道倾角集；
理论轨道倾角集模块，用于对非合作连续推力航天器的理论轨道倾角积分解析表达，得到理论轨道倾角模型，将非初始时刻时间集输入理论轨道倾角模型中，得到理论轨道倾角集；
评价指标模块，用于采用最小二乘的思想，设置理论轨道倾角集与真实轨道倾角集之间的轨道倾角误差集的欧几里得范数作为评价指标；
推力参数辨识模块，用于对作为评价指标的轨道倾角误差集的欧几里得范数进行极值寻优得到最小值，得到非合作连续推力航天器各方向上的推力参数。。
实施例三
本发明还提供了一种非合作连续推力航天器推力参数辨识设备，包括存储器和处理器，所述存储器存储有计算机程序，所述处理器执行所述计算机程序时实现一种非合作连续推力航天器推力参数辨识方法的步骤。
另外，本发明还提供了一种计算机可读存储介质，其存储有计算机程序，所述计算机程序被处理器执行时实现一种非合作连续推力航天器推力参数辨识方法的步骤。
本发明实施例的有益效果是：
本发明的实施例通过将接收到的非合作连续推力航天器的雷达数据转换为非初始时刻的轨道面数据，基于非初始时刻轨道面与赤道平面之间夹角得到轨道倾角模型；将与非初始时刻时间集对应的位置坐标集输入轨道倾角模型中，得到真实轨道倾角集；利用连续推力作用下，航天器的轨道升交点赤经不会发生改变这一特点，将雷达观测数据集转换为轨道倾角模型；
对轨道倾角模型积分解析表达，得到轨道倾角变化模型，将非初始时刻时间集输入轨道倾角变化模型中，得到理论轨道倾角集；解决了连续推力圆轨道的轨道倾角无法解析表达的问题；
采用最小二乘的思想，设置理论轨道倾角集与真实轨道倾角集之间的轨道倾角误差集的欧几里得范数作为评价指标；将非合作连续推力航天器的推力反演问题转化为二元函数求极值的数学问题。
利用遗传算法对轨道倾角误差集的欧几里得范数进行极值寻优，得到非合作连续推力航天器的推力参数。本发明采用遗传算法的优势：
由于本发明针对不能提供有效合作信息的非合作连续推力航天器的推力参数辨识，发明人经过试验，一般的数值迭代方法容易陷入局部极小的陷阱而出现"死循环"现象，使迭代无法进行。遗传算法基本上不用搜索空间的知识或其它辅助信息，而仅用适应度函数值来评估个体，在此基础上进行遗传操作。适应度函数不仅不受连续可微的约束，而且其定义域可以结合本发明针对的非合作连续推力航天器进行设定。遗传算法的搜索从轨道倾角误差集的欧几里得范数群体出发，利用内在并行性进行分布式计算，同时处理群体中的多个个体，即对搜索空间中的多个解进行评价，减少了陷入局部最优解的风险，也加快了求解速度，提高了求解精度。
以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以所述权利要求的保护范围为准。

Claims (9)

1.一种非合作连续推力航天器推力参数辨识方法，其特征在于，该方法包括：
步骤一、将接收到的非合作连续推力航天器的雷达数据转换为非初始时刻的轨道面数据，基于非初始时刻轨道面与赤道平面之间夹角得到真实轨道倾角模型；将与非初始时刻时间集对应的位置坐标集输入真实轨道倾角模型中，得到真实轨道倾角集；
步骤二、将非合作连续推力航天器的理论轨道倾角解析表达，得到理论轨道倾角模型，将非初始时刻时间集输入理论轨道倾角模型中，得到理论轨道倾角集；
步骤三、采用最小二乘的思想，设置理论轨道倾角集与真实轨道倾角集之间的轨道倾角误差集的欧几里得范数作为评价指标；
步骤四、对作为评价指标的轨道倾角误差集的欧几里得范数进行极值寻优得到最小值，得到非合作连续推力航天器各方向上的推力参数。
2.如权利要求1所述的方法，其特征在于，步骤一中，真实轨道倾角模型为：
其中，ii为非初始时刻的真实轨道倾角，为非初始时刻轨道面与赤道平面之间夹角；e表示赤道平面的法向量；hi表示非初始时刻轨道面的法向量；|e|为赤道平面法向量的模；|hi|为非初始时刻轨道面法向量的模。
3.如权利要求1所述的方法，其特征在于，步骤二中，理论轨道倾角模型为：
其中，i为非初始时刻的理论轨道倾角；i0为初始时刻已知的轨道倾角；和/>为推力参数，/>为单位质量所受的切向力，即切向推力加速度；/>为单位质量所受的法向力，即法向推力加速度；μ为地球引力常数；r0为轨道半径；t为推力施加时间。
4.如权利要求1所述的方法，其特征在于，步骤三中，采用最小二乘的思想，设置的理论轨道倾角集与真实轨道倾角集之间的轨道倾角误差集的欧几里得范数为：
其中，Fs为欧几里得范数；t为推力施加时间，t＝(1,2,3,…T)，T为推力施加总时长，T>0；errt为推力施加时间对应的轨道倾角误差；errT为推力施加总时长对应的轨道倾角误差。
5.如权利要求1或4所述的方法，其特征在于，步骤三中，理论轨道倾角集与真实轨道倾角集之间的轨道倾角误差集为：
ERR＝Is-I＝(err1,err2,err3,…errT)；
其中，ERR为轨道倾角误差集；Is为真实轨道倾角集；I为理论轨道倾角集；errT为推力施加总时长对应的轨道倾角误差；T为推力施加总时长，T>0。
6.如权利要求1所述的方法，其特征在于，步骤四中，对作为评价指标的轨道倾角误差集的欧几里得范数进行极值寻优的方法，包括：通过LM算法、遗传算法、粒子群算法或蚁群算法对欧几里得范数进行极值寻优，得到欧几里得范数的最小值。
7.一种非合作连续推力航天器推力参数辨识系统，其特征在于，包括：
真实轨道倾角集模块，用于将接收到的非合作连续推力航天器的雷达数据转换为非初始时刻的轨道面数据，基于非初始时刻轨道面与赤道平面之间夹角得到真实轨道倾角模型；将与非初始时刻时间集对应的位置坐标集输入真实轨道倾角模型中，得到真实轨道倾角集；
理论轨道倾角集模块，用于对非合作连续推力航天器的理论轨道倾角积分解析表达，得到理论轨道倾角模型，将非初始时刻时间集输入理论轨道倾角模型中，得到理论轨道倾角集；
评价指标模块，用于采用最小二乘的思想，设置理论轨道倾角集与真实轨道倾角集之间的轨道倾角误差集的欧几里得范数作为评价指标；
推力参数求解模块，用于对作为评价指标的轨道倾角误差集的欧几里得范数进行极值寻优得到最小值，得到非合作连续推力航天器各方向上的推力参数。
8.一种非合作连续推力航天器推力参数辨识设备，包括存储器和处理器，所述存储器存储有计算机程序，所述处理器执行所述计算机程序时实现如权利要求1至6中任意一项所述的方法的步骤。
9.一种计算机可读存储介质，其特征在于，其存储有计算机程序，所述计算机程序被处理器执行时实现如权利要求1至6中任意一项所述的方法的步骤。
CN202310897508.5A 2023-07-20 2023-07-20 非合作连续推力航天器推力参数辨识方法、系统和设备 Active CN116954075B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202310897508.5A CN116954075B (zh) 2023-07-20 2023-07-20 非合作连续推力航天器推力参数辨识方法、系统和设备

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202310897508.5A CN116954075B (zh) 2023-07-20 2023-07-20 非合作连续推力航天器推力参数辨识方法、系统和设备

Publications (2)

Publication Number Publication Date
CN116954075A true CN116954075A (zh) 2023-10-27
CN116954075B CN116954075B (zh) 2024-04-19

Family

ID=88447155

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202310897508.5A Active CN116954075B (zh) 2023-07-20 2023-07-20 非合作连续推力航天器推力参数辨识方法、系统和设备

Country Status (1)

Country Link
CN (1) CN116954075B (zh)

Citations (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2016063923A1 (ja) * 2014-10-24 2016-04-28 株式会社アストロスケール 非協力接近に関する誘導方法
CN105737847A (zh) * 2014-12-09 2016-07-06 上海新跃仪表厂 非合作目标条件下闭环自主导航的试验系统
CN109165415A (zh) * 2018-07-28 2019-01-08 西北工业大学 一种基于人工合成引力势场的连续推力轨道设计方法及其应用
CN109190158A (zh) * 2018-07-26 2019-01-11 西北工业大学 一种考虑非合作目标禁飞区约束的最优轨道设计方法
CN110850719A (zh) * 2019-11-26 2020-02-28 北京航空航天大学 一种基于强化学习的空间非合作目标参数自整定追踪方法
RU2775092C1 (ru) * 2021-12-27 2022-06-28 Федеральное государственное автономное образовательное учреждение высшего образования "Омский государственный технический университет" Способ увода объектов крупногабаритного космического мусора и устройство для его реализации
CN114815869A (zh) * 2022-04-29 2022-07-29 西北工业大学 一种基于星群指纹库的非合作航天器运动参数识别方法
CN115935519A (zh) * 2022-12-28 2023-04-07 中国科学院数学与系统科学研究院 飞行器纵向气动参数的分组加权在线最小二乘辨识方法
CN116050133A (zh) * 2023-01-10 2023-05-02 西北工业大学 抵近空间站的非合作目标协同探测及参数识别方法及系统
FR3129237A1 (fr) * 2021-11-17 2023-05-19 Airbus Defence And Space Sas Procédé d'acquisition d'images d'un objet spatial en orbite terrestre par un engin spatial en orbite terrestre

Patent Citations (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2016063923A1 (ja) * 2014-10-24 2016-04-28 株式会社アストロスケール 非協力接近に関する誘導方法
CN105737847A (zh) * 2014-12-09 2016-07-06 上海新跃仪表厂 非合作目标条件下闭环自主导航的试验系统
CN109190158A (zh) * 2018-07-26 2019-01-11 西北工业大学 一种考虑非合作目标禁飞区约束的最优轨道设计方法
CN109165415A (zh) * 2018-07-28 2019-01-08 西北工业大学 一种基于人工合成引力势场的连续推力轨道设计方法及其应用
CN110850719A (zh) * 2019-11-26 2020-02-28 北京航空航天大学 一种基于强化学习的空间非合作目标参数自整定追踪方法
FR3129237A1 (fr) * 2021-11-17 2023-05-19 Airbus Defence And Space Sas Procédé d'acquisition d'images d'un objet spatial en orbite terrestre par un engin spatial en orbite terrestre
RU2775092C1 (ru) * 2021-12-27 2022-06-28 Федеральное государственное автономное образовательное учреждение высшего образования "Омский государственный технический университет" Способ увода объектов крупногабаритного космического мусора и устройство для его реализации
CN114815869A (zh) * 2022-04-29 2022-07-29 西北工业大学 一种基于星群指纹库的非合作航天器运动参数识别方法
CN115935519A (zh) * 2022-12-28 2023-04-07 中国科学院数学与系统科学研究院 飞行器纵向气动参数的分组加权在线最小二乘辨识方法
CN116050133A (zh) * 2023-01-10 2023-05-02 西北工业大学 抵近空间站的非合作目标协同探测及参数识别方法及系统

Non-Patent Citations (5)

* Cited by examiner, † Cited by third party
Title
卫晓娜;董云峰;郝朝;: "基于遗传编程的航天器有限推力逼近轨迹规划", 北京航空航天大学学报, no. 05 *
崔红正;刘文玲;唐歌实;宋柏延;葛茂荣;: "不同推力下的非合作空间目标轨道机动检测", 宇航学报, no. 03, 30 March 2016 (2016-03-30) *
张莹;王西京;袁博;孔大林;卞燕山;: "基于夏氏最小二乘的轨道控制力系数辨识", 中国空间科学技术, no. 02, 27 March 2019 (2019-03-27) *
潘迅;泮斌峰;唐硕;: "考虑J2项摄动的小推力燃料最优转移轨道设计", 哈尔滨工业大学学报, no. 10 *
郭延宁;韩旭;郭增千;张瑶;马广富;: "基于双目视觉的非合作目标逼近控制系统设计与仿真", 空间控制技术与应用, no. 05, 15 October 2015 (2015-10-15) *

Also Published As

Publication number Publication date
CN116954075B (zh) 2024-04-19

Similar Documents

Publication Publication Date Title
CN103528587B (zh) 自主组合导航系统
CN102156478B (zh) 一种基于蚁群Unscented粒子滤波算法的组合定姿方法
CN111102981B (zh) 一种基于ukf的高精度卫星相对导航方法
Ahn et al. Fast alignment using rotation vector and adaptive Kalman filter
CN109901598A (zh) 基于随机模型预测控制技术的自主水下机器人路径跟踪方法
CN101059349A (zh) 微型组合导航系统及自适应滤波方法
CN109708663B (zh) 基于空天飞机sins辅助的星敏感器在线标定方法
CN109656133B (zh) 一种针对空间走廊跟踪观测的分布式卫星群优化设计方法
CN106772524A (zh) 一种基于秩滤波的农业机器人组合导航信息融合方法
CN110262241B (zh) 基于高斯过程预测控制的航天器轨道控制方法
CN115523927B (zh) 基于光学传感器观测的geo航天器机动检测方法
CN107300700B (zh) 敏捷合成孔径雷达卫星聚束模式姿态机动需求计算方法
CN114580224A (zh) 一种分布式气动融合轨道耦合姿态摄动分析方法
CN112731281B (zh) 一种空间碎片测角数据仿真方法
US12078716B2 (en) System and method of hypersonic object tracking
CN112161632A (zh) 一种基于相对位置矢量测量的卫星编队初始定位算法
CN104729510A (zh) 一种空间目标相对伴飞轨道确定方法
CN110146092B (zh) 基于导航信息评价的双体小行星探测轨迹优化方法
CN115314101A (zh) 一种基于并行计算的低轨通信卫星星座快速建模方法
CN111637895B (zh) 一种基于q学习的导航观测目标选取方法
CN109781374A (zh) 一种实时在线快速估计飞行器推力的方法
CN116954075B (zh) 非合作连续推力航天器推力参数辨识方法、系统和设备
CN110231619B (zh) 基于恩克法的雷达交接时刻预报方法及装置
CN115421483B (zh) 一种无人船操纵运动预报方法
CN115270643A (zh) 基于观测能效最优的光测设备空间目标测量站址选取系统及方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant