CN105021311B - 一种在轨卫星推力器温度异常实时诊断方法 - Google Patents

Abstract

本发明公开了一种在轨卫星推力器温度异常实时诊断方法。本发明根据在轨卫星历史遥测数据，确定四种工况下推力器温度的物理模型中的参数：向深冷空间的热辐射系数、推力器单位质量因推力器的加热器加热产生的温度增量、推力器单位质量因推进剂燃烧产生的温度增量和推力器单位质量因太阳光照产生的温度增量；根据四种工况下推力器温度的物理模型，对当前监测时刻所处工况下的推力器温度进行实时预测，并根据卫星实际情况外扩一定范围后得到动态的报警门限，实现对轨控过程中推力器的工作状态的实时监测和异常温度报警，进行推力器温度异常实时诊断和预警研究。

Description

一种在轨卫星推力器温度异常实时诊断方法

技术领域

本发明涉及在轨卫星推进系统，具体涉及一种在轨卫星推力器温度异常实时诊断方法。

背景技术

目前我国中高轨道在轨卫星推进系统均采用双组元统一推进系统，10N推力器是推进分系统中提供推力的唯一执行部件，其工作正常与否直接决定了卫星轨道控制的成败。现今主要通过两种手段开展10N推力器的状态监测：(1)通过卫星姿态角变化进行监测；(2)通过控制结束后轨道测量结果进行评估。上述两种手段存在以下不足：

(1)一般地，卫星成对推力器的推力存在差异，正常情况下采用成对推力器点火也会引起姿态波动，并且在推力器点火过程中如果卫星动量装置参与姿态控制，则推力器喷气产生的干扰力矩会逐渐被动量装置吸收，推力器异常工况无法完全反映到姿态波动上，因此采用姿态角变化监测方法无法发现引起姿态波动与正常情况姿态波动差异不明显的推力器异常工况。

(2)采用控制结束后轨道测量结果评估虽然能够发现10N推力器所有异常，但缺乏实时性。例如某通信卫星某次轨控10N推力器发生故障，在轨道控制过程中，对卫星姿态角进行监测，未发现异常，但控制结束后通过轨道测量发现轨道未达到预期目标，进一步分析确定是因为10N推力器异常，后续虽然通过补充控制使卫星轨道达到预期目标，但是延长了卫星的轨控时间，影响了卫星点波束天线等有效载荷的正常使用。因此，需要寻找一种方法解决10N推力器异常发现不及时的难题。

10N推力器温度是唯一可以实时全面表征推力器工作正常与否的重要参数。无论是贮箱问题，还是10N推力器自身问题，甚至是电磁阀管路问题，其直接表现均为推力不稳定或推力减小，最终均能通过10N推力器温度表征。此外，10N推力器温度变化趋势也是反映推进剂是否排空或变为单组元的直接判据，是监测管理寿命末期卫星的重要依据。但现今仍依靠人工对10N推力器温度趋势进行判断分析，这种方式存在以下问题：首先，面对十二五末期上百颗在轨卫星的局面，依靠人工不可能完成实时判断分析；其次，这种方式得到的结论严重依赖在轨监测人员的经验，不利于早期发现10N推力器的故障或性能下降。

发明内容

有鉴于此，本发明提供了一种在轨卫星推力器温度异常实时报警方法，能够在任意时刻，对各种工况下的推力器温度进行预测，并根据卫星实际情况外扩一定范围后得到动态的报警门限，实现对轨控过程中推力器的工作状态的实时监测，进行推力器温度异常实时诊断和预警研究。

本发明解决其技术问题所采用的技术方案是：

一种在轨卫星推力器温度异常实时诊断方法，具体包括如下步骤：

步骤一、构建推力器温度的物理模型：

T(t+Δt)＝T(t)+ΔT_太阳(t,t+Δt)+ΔT_推进剂+ΔT_加热器-μη(T(t))⁴Δt (1)

其中，t表示当前监测时间，T(t)表示t时刻推力器的温度，T(t+Δt)表示下一个时刻推力器温度的实时预测值，Δt为监测推力器温度的时间步长；

η表示向深冷空间的热辐射系数，μ为修正系数；

ΔT_推进剂表示推进剂燃烧产生的温度增量；

ΔT_加热器表示推力器的加热器加热产生的温度增量；

ΔT_太阳(t,t+Δt)表示太阳光照产生的温度增量；

步骤二、根据公式(1)对不同工况下的推力器温度进行预测和参数拟合：

工况一：推力器未点火、推力器的加热器关闭，并且不受太阳光照，推力器温度的物理模型为：T(t+Δt)＝T(t)-η(T(t))⁴Δt (2)

选择与工况一最相近的推力器温度的在轨遥测数据，并代入公式(2)中确定向深冷空间的热辐射系数η；

工况二：推力器点火、推力器的加热器关闭，并且不受太阳光照，推力器温度的物理模型为：T(t+Δt)＝T(t)+ΔT_推进剂-μη(T(t))⁴Δt (3)

选择与工况二最相近的推力器温度的在轨遥测数据，并代入公式(3)中确定参数ΔT_推进剂和修正系数μ；

工况三：推力器未点火、推力器的加热器打开，并且不受太阳光照，推力器温度的物理模型为：T(t+Δt)＝T(t)+ΔT_加热器-η(T(t))⁴Δt (4)

选择与工况三最相近的推力器温度的在轨遥测数据，并代入公式(4)中确定参数ΔT_加热器；

工况四：推力器未点火、推力器的加热器关闭，并且受太阳光照，推力器温度的物理模型为：T(t+Δt)＝T(t)+ΔT_太阳(t,t+Δt)-η(T(t))⁴Δt (5)

选择与工况四最相近，并且太阳光照条件与在轨卫星当前位置相同的推力器温度的在轨遥测数据，并代入公式(5)中确定ΔT_太阳(t,t+Δt)；

步骤三、推力器温度实时预测及异常报警门限的设置：

判断当前监测时刻的工况，根据步骤二中该工况下推力器温度的物理模型得到推力器温度的实时预测值，并将该实时预测值外扩±(3～5)℃作为推力器温度异常的报警门限；

步骤四、推力器温度异常实时诊断：

当监测到的推力器温度在所述报警门限的范围内，则推力器工作正常；

当监测到的推力器温度高于所述报警门限的范围时，确定推进剂的燃烧效率高的原因，当温度持续过高时，停止推力器工作；

当监测到的推力器温度低于所述报警门限的范围时，确定推进剂的燃烧效率低的原因，及时进行推力器温度异常的处理。

此外，监测推力器温度的时间步长Δt满足：t₁＜Δt＜t₂；

其中，t₁表示推力器温度变化一个分层值所需的时间；t₂表示推力器温度变化二个分层值所需的时间。

在工况一时，选取≥5组的推力器温度的在轨遥测数据进行拟合，得到向深冷空间的热辐射系数η；在工况三时，选取≥5组的推力器温度的在轨遥测数据进行拟合，得到ΔT_加热器。

在工况二下，选取的推力器温度的在轨遥测数据，将推力器点火初期1min～1.5min的数据剔除后再利用公式(3)确定参数ΔT_推进剂和修正系数μ。

有益效果：

(1)本发明根据推力器的热平衡方程，考虑空间外热流等影响因素，应用推力器在轨实测温度数据，通过系数修正和实时解算的方法建立了精确度较高的在轨卫星推力器温度数学模型，进一步提供了一种在轨卫星推力器温度异常实时诊断方法；在卫星姿态稳定的前提下，本发明中得到的推力器温度拟合公式仅与时间相关，利用该公式预测推力器任意时刻的温度，能够对轨控过程中推力器的工作状态进行实时监测，填补了目前在轨卫星轨控过程中推力器异常无法实时自主发现的空白，简单高效、可行可靠。

(2)由本发明在推力器温度预测值的基础上，根据卫星实际情况经过外扩一定范围后得到动态的报警门限，适用于在轨监视工作、大幅优化了原有的报警上下限，可在卫星寿命末期进行推进剂夹气监视，以及作为卫星离轨推进剂排空的判据等。

(3)本发明采用集总思路，方法简单便于工程实现；本发明的拟合方法是基于卫星在轨真实的遥测数据进行的半物理方法，与物理仿真、测试等数据相比，卫星在轨遥测数据能更真实地反映卫星在轨性能的变化情况，而且经过在轨遥测数据对物理模型进行修正，大幅提高拟合精度，具有一定的普适性，可应用于在轨卫星1N、5N和10N等推力器温度异常实时诊断和预警研究，也可应用于性能分析。

附图说明

图1为在轨卫星推力器温度异常实时诊断方法的流程图。

图2为某在轨卫星10N推力器未点火、推力器的加热器关闭，并且不受太阳光照时的10N推力器温度变化图；

图3为某在轨卫星10N推力器未点火、推力器的加热器关闭，并且不受太阳光照时的10N推力器实际温度与拟合温度的关系图，其中散点为10N推力器实际温度，光滑曲线为10N推力器拟合温度；

图4为某在轨卫星10N推力器未点火、推力器的加热器关闭，并且不受太阳光照时的10N推力器拟合温度误差图；

图5为某在轨卫星10N推力器点火、推力器的加热器关闭，并且不受太阳光照时的10N推力器温度变化图；

图6为某在轨卫星10N推力器点火、推力器的加热器关闭，并且不受太阳光照时10N推力器实际温度与拟合点火温度的关系图，其中散点为10N推力器实际温度，光滑曲线为10N推力器拟合温度；

图7为某在轨卫星10N推力器点火、推力器的加热器关闭，并且不受太阳光照时的10N推力器拟合温度误差图；

图8为某在轨卫星10N推力器未点火、推力器的加热器打开、并且不受太阳光照时的10N推力器温度变化图；

图9为某在轨卫星10N推力器未点火、推力器的加热器打开，并且不受太阳光照时的10N推力器实际温度与拟合温度的关系图，其中散点为10N推力器实际温度，光滑曲线为10N推力器拟合温度；

图10为某在轨卫星10N推力器未点火、推力器的加热器打开、并且不受太阳光照时的10N推力器拟合温度误差图；

图11为某在轨卫星10N推力器未点火、推力器的加热器关闭，并且受太阳光照且时的10N推力器温度变化图；

图12为某在轨卫星10N推力器未点火、推力器的加热器关闭，并且受太阳光照时的10N推力器实际温度与拟合温度的关系图，其中散点为10N推力器实际温度，光滑曲线为10N推力器拟合温度；

图13为某在轨卫星10N推力器未点火、推力器的加热器关闭，并且受太阳光照时的10N推力器拟合温度误差图；

图14为某在轨卫星某时间段10N推力器点火、推力器的加热器关闭，并且受太阳光照时的10N推力器温度异常报警门限示意图；其中，直线a1和a2为设计门限，曲线b1和b2为本发明的动态门限，曲线c为设计门限无法识别的温度异常变化曲线，曲线d为设计门限不能及时识别的温度异常曲线，t_d为报警滞后时间。

具体实施方式

下面结合附图并举实施例，对本发明进行详细描述。

本实施例以10N推力器为例，10N推力器点火过程温度遥测的影响因素主要包括：空间外热流(包括帆板红外热流及卫星本体遮挡等)、10N推力器的加热器通断状态、推进剂燃烧所产生的热耗和测温热敏电阻特性(精度及遥测分层值误差)。若想通过在轨遥测数据将所有因素一一剥离，得到10N推力器温度的精确数学模型进行预测是极难实现的。

本发明为解决10N推力器温度精细化监测需求与精确数学模型难以构建之间的矛盾，在考虑了以上各影响因素的基础上，提出了一种经过简化的10N推力器温度的半物理模型，并在此基础上构建10N推力器温度的动态报警门限，可服务于轨控过程中的精细化监测。本发明提出的方法逻辑复杂度适中，计算量适中，便于工程实现。

本发明所提出的在轨卫星推力器温度异常实时诊断方法的流程如附图1所示，具体步骤如下：

步骤一、构建推力器温度的物理模型：

对于推进系统采用双组元统一推进系统，且推力器采用恒定脉宽或者连续方式点火的在轨卫星，10N推力器的温度的理想物理模型：采用集总思路，将10N推力器假设为均温物体，即10N推力器各部分的比热容C相同；空间外热流仅考虑太阳光照，同时考虑推进剂燃烧、推力器的加热器加热及10N推力器向深冷空间的热辐射，则10N推力器温度的理想物理模型满足：

其中：

M为10N推力器质量；T为集总温度，表示10N推力器的平均温度，单位开尔文(K)；t为时间，单位秒(s)；

q_太阳(t)是时间的函数，为太阳光照引起的10N推力器温度变化率；

q_推进剂为常数，表示推进剂燃烧引起的10N推力器温度变化率；

q_加热器为常数，表示推力器的加热器加热引起的10N推力器温度变化率；

σ为斯蒂芬-波尔兹曼常数，ε为10N推力器表面发射率，A为10N推力器的表面积，C为推力器的比热容，表示10N推力器向深冷空间的放热速率。

将公式(1)进行离散化积分得到：

其中：

Δt为相邻两个温度值之间的时间间隔，即监测推力器温度的时间步长，单位为秒，Δt是综合热敏电阻的灵敏度、卫星遥测参数下传的时间间隔及卫星在轨遥测数据设置的，要求Δt满足：t₁＜Δt＜t₂；

其中，t₁表示10N推力器温度变化一个分层值所需的时间；t₂表示10N推力器温度变化二个分层值所需的时间；即Δt的选择可以识别出10N推力器温度的变化；

T(t)为时间t的温度函数，表示t时刻10N推力器温度；T(t+Δt)表示t+Δt时刻10N推力器温度；

ΔT_太阳(t,t+Δt)表示t到t+Δt时间推力器单位质量因太阳光照产生的温度增量；

ΔT_推进剂表示t到Δt时间推力器单位质量因推进剂燃烧产生的温度增量；

ΔT_加热器表示t到Δt时间推力器单位质量因推力器的加热器加热产生的温度增量；

为Δt时间推力器单位质量因向深冷空间热辐射产生的温度增量，其中为向深冷空间的热辐射系数，用η代替。

然而，实际上，10N推力器是非均温体，太阳光照、推进剂燃烧、推力器的加热器加热的不均匀性必然导致10N推力器本体温度场存在梯度。尤其当推力器点火时，两种燃料在燃烧室混合并开始燃烧、通过拉伐尔喷管喷出，因此推力器末端温度很高。由于推力器热容及自身导热热阻的存在，使得点火过程中沿推力器轴向必然形成较大温度梯度。在轨卫星10N推力器热敏电阻安装于法兰一端，所以在温度梯度较大的情况下，温度测点所下传的温度便不能代表推力器的平均温度，即集总温度，因此需要对公式(2)中η(T(t))⁴Δt项进行修正。引入修正系数μ，则10N推力器温度物理模型为：

T(t+Δt)＝T(t)+ΔT_太阳(t,t+Δt)+ΔT_推进剂+ΔT_加热器-μη(T(t))⁴Δt (3)

步骤二、根据公式(3)和不同工况下的卫星在轨遥测数据对推力器温度进行预测和参数拟合，下面以某颗在轨卫星10N推力器4A为例：

根据Δt需要满足的要求，在本实施例中设置时间步长Δt＝24s，即下面步骤中10N推力器温度在轨遥测数据选择时，数据的采样间隔为24s。

工况一：确定常值系数η

工况一为推力器未点火、推力器的加热器关闭，并且不受太阳光照，即在公式(3)中，ΔT_太阳(t,t+Δt)、ΔT_推进剂和ΔT_加热器为0，由于推力器还未点火，所以不考虑修正系数μ的影响，则工况一的推力器温度的物理模型为：

T(t+Δt)＝T(t)-η(T(t))⁴Δt (4)

选择一个与工况一(推力器的工作状态和工作环境)最相近的轨道周期的数据，即选取其中推力器未点火、推力器的加热器关闭且不受照时的10N推力器温度的在轨遥测数据，利用公式(4)进行拟合，得到向深冷空间的热辐射系数η。由于不受照时10N推力器的加热器控温门限较窄，推力器的加热器关闭的时间较短，因此需要选取≥5组的推力器的加热器关闭时的推力器温度数据进行拟合，以便提高常值系数η的拟合精度。

选择多组数据拟合得到多个常值系数η，可以根据历史经验选择最接近10N推力器实际工况的拟合值；或对拟合得到的多个η求取均值，作为最终的拟合值。

本实施例中，连续选择7组，推力器未点火、推力器的加热器关闭且不受太阳光照时的10N推力器温度的在轨遥测数据，按时间先后顺序连接起来如图2所示。对所选7组数据，独立选择自变量T(t)和因变量T(t+Δt)，利用公式(4)拟合常值系数η。

拟合曲线和拟合误差分别如图3和图4所示，η的拟合值为2.858×10^-12，拟合曲线中的最大的温度误差为0.62℃。

工况二：拟合推进剂燃烧产生的温度增量ΔT_推进剂和修正系数μ

工况二为推力器点火、推力器的加热器关闭，且不受太阳光照，即在公式(3)中，ΔT_太阳(t,t+Δt)和ΔT_加热器为0，考虑修正系数μ的影响，使用工况一中拟合的常值系数η，则工况二下10N推力器温度的物理模型为：

T(t+Δt)＝T(t)+ΔT_推进剂-μη(T(t))⁴Δt (5)

选择一个与工况二(推力器的工作状态和工作环境)最相近的轨道周期的数据，即选取其中推力器点火、推力器的加热器关闭，且不受太阳光照时的10N推力器温度在轨遥测数据，利用公式(5)进行曲线拟合，得到ΔT_推进剂和修正系数μ。由于推力器自身热容的影响，点火开始后推进剂燃烧产生的热耗需经过一定时间才能传至推力器法兰端热敏电阻处，因此在进行曲线拟合时，10N推力器点火初期1～1.5min的数据需剔除。

本实施例中，选取不受太阳光照，推力器的加热器关闭，10N推力器点火时的推力器温度正常变化的在轨遥测数据，将点火初期1.5分钟的10N推力器温度数据去掉，如图5所示。应用公式(5)拟合推进剂燃烧产生的温度增量ΔT_推进剂和修正系数μ。

拟合曲线和拟合误差分别如图6和图7所示，推进剂燃烧产生的温度增量ΔT_推进剂＝4.427，修正系数μ＝4.003。拟合曲线的最大温度误差为-1.68℃。

工况三：拟合推力器的加热器加热产生的温度增量ΔT_加热器

工况三为推力器未点火、推力器的加热器打开，并且不受太阳光照，即在公式(3)中，ΔT_太阳(t,t+Δt)和ΔT_推进剂为0，不考虑修正系数μ的影响，使用工况一中拟合的常值系数η，则10N推力器温度的物理模型为：

T(t+Δt)＝T(t)+ΔT_加热器-η(T(t))⁴Δt (6)

选择一个与工况三(推力器的工作状态和工作环境)最相近的轨道周期的数据，选取其中推力器未点火，推力器的加热器打开，并且不受太阳光照时的10N推力器温度的在轨遥测数据，利用公式(6)进行曲线拟合，可得到ΔT_加热器。由于不受照时10N推力器的加热器控温门限较窄，10N推力器的加热器打开的时间较短，因此需要选取≥5组以上的10N推力器的加热器打开时的推力器温度数据，以便提高ΔT_加热器的拟合精度。

选择多组数据拟合得到多个ΔT_加热器，可以根据历史经验选择最接近10N推力器实际工况的拟合值；或对拟合得到的多个ΔT_加热器求取均值，作为最终的拟合值。

本实施例中，连续选择6组推力器未点火、推力器的加热器打开且不受太阳光照时的10N推力器温度在轨遥测数据，按时间先后顺序连接起来如图8所示。对所选6组数据，独立选择自变量T(t)和因变量T(t+Δt)，利用公式(6)拟合推力器的加热器加热产生的温度增量ΔT_加热器。

拟合曲线和拟合误差分别如图9和图10所示，推力器的加热器加热产生的温度增量ΔT_加热器＝1.198，拟合曲线的最大温度误差为-1.01℃。

工况四：计算太阳光照产生的温度增量ΔT_太阳(t,t+Δt)

工况四为推力器未点火、推力器的加热器关闭，并且受太阳光照，即在公式(3)中，ΔT_加热器和ΔT_推进剂为0，不考虑修正系数μ的影响，使用工况一中拟合的常值系数η，则10N推力器温度的物理模型为：

T(t+Δt)＝T(t)+ΔT_太阳(t,t+Δt)-η(T(t))⁴Δt (7)

为了减弱轨道漂移对10N推力器温度变化周期的影响，选择一个与工况四(推力器的工作状态和工作环境)最相近的轨道周期的数据，选取其中推力器未点火、推力器的加热器关闭，受太阳光照，并且太阳光照条件与在轨卫星当前位置相同的10N推力器温度在轨遥测数据，按照公式(7)计算ΔT_太阳(t,t+Δt)；

步骤三、推力器温度的实时预测以及异常报警门限的设置：

在工况一到工况四下对各个工况下物理模型中的参数进行拟合后，得到的10N推力器温度是时间的函数，根据监测时的工况以及时间，根据相应工况下的推力器温度的物理模型，可得到10N推力器温度的实时预测值。

实时预测：选取某次10N推力器工作时的在轨温度遥测数据，其工况为10N推力器点火、推力器的加热器关闭且受太阳光照，如图11所示；基于上述的拟合结果，计算10N推力器温度监测时间的预测值，如图12所示。10N推力器温度实际值与预测值的误差如图13所示，最大温度误差为-2.06℃。

异常报警门限选取：判断当前监测时刻的工况，根据步骤二中该工况下推力器温度的物理模型得到推力器温度的实时预测值，然后根据在轨遥测数据的波动情况，将实时预测值外扩±(3～5)℃作为推力器温度异常的报警门限，即选取T(t+Δt)±(3～5)℃作为异常报警门限，形成适用于该卫星10N推力器温度的正常值范围，进行在轨异常发现。

如图14所示，为某时间段10N推力器受太阳光照、推力器的加热器关闭时的报警门限情况。其中，二条水平直线a1和a2为设计门限，曲线b1和b2为本发明的动态门限，曲线c为设计门限无法识别的温度异常变化曲线，曲线d为设计门限不能及时识别的温度异常曲线，t_d为报警滞后时间。

从图14中，可以观察到，设计提供的正常值范围门限(简称设计门限)范围较宽；而本发明提供的动态报警门限阈值的跨度仅为5℃，大大缩小了报警门限范围。此外，本发明提供的动态报警门限能够发现设计门限无法识别的温度异常；并且能够相对于设计门限能够提前发现异常，为异常处理争取宝贵时间；因此本发明得出的动态报警门限能够更贴切地反映10N推力器温度的真实变化情况，从而更准确合理地进行在轨监视和异常发现。

步骤四、推力器温度异常实时诊断：

当监测到的10N推力器温度在所述报警门限的范围内，判断推力器工作正常；

当监测到的10N推力器温度高于所述报警门限的范围时，分析推进剂的燃烧效率高的原因，当温度持续过高时，停止推力器工作；

当监测到的10N推力器温度低于所述报警门限的范围时，判断推进剂的燃烧效率低的原因，及时进行异常分析处理，例如：排查推进剂是否用完，检查电磁阀工作是否异常，及时对推力器效率系数进行修正，进一步增加控制量进行轨道机动补充控制。

本发明的推力器温度实时预测方法对于在轨卫星普遍适用，该方法已利用多颗卫星的在轨遥测数据进行验证。现再举例给出对2颗在轨卫星的预测结果，分别称为卫星1、卫星2。

卫星1：向西位置保持过程中，10N推力器4A点火，推力器的加热器关闭，并且不受太阳光照的情况下，10N推力器4A温度拟合公式为：

T(t)＝T(t-24)+4.427-2.746×10^-10(T(t-24)+273)⁴

拟合绝对误差为[-1.68℃,1.33℃]。

卫星2：向东位置保持过程中，10N推力器2A点火，推力器的加热器关闭，并且不受太阳光照的情况下，10N推力器2A温度拟合公式为：

T(t)＝T(t-24)+5.398-3.409×10^-10(T(t-24)+273)⁴

拟合绝对误差为[-1.71℃,1.33℃]。

综上所述，以上仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (4)

1.一种在轨卫星推力器温度异常实时诊断方法，其特征在于，包括：

步骤一、构建推力器温度的物理模型：

T(t+Δt)＝T(t)+ΔT_太阳(t,t+Δt)+ΔT_推进剂+ΔT_加热器-μη(T(t))⁴Δt (1)

其中，t表示当前监测时间，T(t)表示t时刻推力器的温度，T(t+Δt)表示下一个时刻推力器温度的实时预测值，Δt为监测推力器温度的时间步长；

η表示向深冷空间的热辐射系数，μ为修正系数；

ΔT_推进剂表示推进剂燃烧产生的温度增量；

ΔT_加热器表示推力器的加热器加热产生的温度增量；

ΔT_太阳(t,t+Δt)表示太阳光照产生的温度增量；

步骤二、根据公式(1)对不同工况下的推力器温度进行预测和参数拟合：

工况一：推力器未点火、推力器的加热器关闭，并且不受太阳光照，推力器温度的物理模型为：T(t+Δt)＝T(t)-η(T(t))⁴Δt (2)

选择与工况一最相近的推力器温度的在轨遥测数据，并代入公式(2)中确定向深冷空间的热辐射系数η；

工况二：推力器点火、推力器的加热器关闭，并且不受太阳光照，推力器温度的物理模型为：T(t+Δt)＝T(t)+ΔT_推进剂-μη(T(t))⁴Δt (3)

选择与工况二最相近的推力器温度的在轨遥测数据，并代入公式(3)中确定参数ΔT_推进剂和修正系数μ；

工况三：推力器未点火、推力器的加热器打开，并且不受太阳光照，推力器温度的物理模型为：T(t+Δt)＝T(t)+ΔT_加热器-η(T(t))⁴Δt (4)

选择与工况三最相近的推力器温度的在轨遥测数据，并代入公式(4)中确定参数ΔT_加热器；

工况四：推力器未点火、推力器的加热器关闭，并且受太阳光照，推力器温度的物理模型为：T(t+Δt)＝T(t)+ΔT_太阳(t,t+Δt)-η(T(t))⁴Δt (5)

选择与工况四最相近，并且太阳光照条件与在轨卫星当前位置相同的推力器温度的在轨遥测数据，并代入公式(5)中确定ΔT_太阳(t,t+Δt)；

步骤三、推力器温度实时预测及异常报警门限的设置：

判断当前监测时刻的工况，根据步骤二中该工况下推力器温度的物理模型得到推力器温度的实时预测值，并将该实时预测值外扩±(3～5)℃作为推力器温度异常的报警门限；

步骤四、推力器温度异常实时诊断：

当监测到的推力器温度在所述报警门限的范围内，则推力器工作正常；

当监测到的推力器温度高于所述报警门限的范围时，确定推进剂的燃烧效率高的原因，当温度持续过高时，停止推力器工作；

当监测到的推力器温度低于所述报警门限的范围时，确定推进剂的燃烧效率低的原因，及时进行推力器温度异常的处理。

2.如权利要求1所述的一种在轨卫星推力器温度异常实时诊断方法，其特征在于，监测推力器温度的时间步长Δt满足：t₁＜Δt＜t₂；

其中，t₁表示推力器温度变化一个分层值所需的时间；t₂表示推力器温度变化二个分层值所需的时间。

3.如权利要求1所述的一种在轨卫星推力器温度异常实时诊断方法，其特征在于，在工况一时，选取≥5组的推力器温度的在轨遥测数据进行拟合，得到向深冷空间的热辐射系数η；在工况三时，选取≥5组的推力器温度的在轨遥测数据进行拟合，得到ΔT_加热器。

4.如权利要求1所述的一种在轨卫星推力器温度异常实时诊断方法，其特征在于，在工况二下，选取的推力器温度的在轨遥测数据，将推力器点火初期1min～1.5min的数据剔除后再利用公式(3)确定参数ΔT_推进剂和修正系数μ。