动力系统故障下运载火箭推力与质量辨识方法

2023-07-24于海森谭述君毛玉明

宇航学报 2023年6期

于海森，谭述君，2，刘浩，毛玉明

(1. 大连理工大学工业装备结构分析优化与CAE软件全国重点实验室，大连 116024；2. 大连理工大学辽宁省空天飞行器前沿技术重点实验室，大连 116024；3. 上海宇航系统工程研究所，上海 201109)

0 引言

运载火箭是航天运输系统的主要组成部分,是目前人类进入空间的主要工具,是发展空间技术、确保空间安全的基石[1]。运载火箭系统结构复杂,发射成本高,确保其安全可靠非常重要。在国内外失败的发射案例中,相当一部分由动力系统故障[2]引起。动力系统故障一般表现为推力下降或关机[3]。运载火箭某台发动机的推力下降会产生推力不平衡干扰力矩[4],而推进剂消耗速度发生变化则导致质量、质心等模型参数与预设值偏差变大,同时推力下降也导致发动机摆角效率的下降,即引起姿控系统伺服机构故障,从而引发严重的控制问题。事实上,动力系统发生故障有时并不是灾难性的,即发生故障后,火箭仍具有继续飞行并完成或部分完成发射任务的能力[5-6]。因此,及时地诊断出火箭推力大小和火箭当前质量有助于后续的控制重构和飞行任务调整[7]。

故障诊断常用的方法有基于信号处理的诊断方法[8]、基于知识的诊断方法[9]以及基于模型的故障诊断方法[10]。火箭的自主诊断技术主要采用基于模型的分析方法[11]。基于模型的方法已经得到了广泛的研究,可细分为参数估计方法[12]、状态估计方法[13]、等价空间法[14]等。在参数估计领域中,最小二乘法是一种得到最广泛应用的参数估计方法,可用于动态、静态、线性和非线性系统等[15]。例如文献[16]采用渐消记忆的递推最小二乘法对火箭推力参数进行了辨识。卡尔曼滤波是处理随机估计应用最广泛的方法[17],为扩展适用范围,卡尔曼滤波出现了多种改进形式,包括扩展卡尔曼滤波、无迹卡尔曼滤波[18]、UD分解滤波等[19]。例如文献[20]基于卡尔曼滤波和故障因子设计了一种针对液体火箭发动机稳态故障的检测与诊断方法。文献[21]提出了线性-二次后退时域算法进行故障检测,该算法利用扩展卡尔曼滤波获取残差来估计推力损失程度,同时利用过载故障方向定位方法对故障发动机定位。

针对运载火箭推力参数辨识的方向主要有两个方面:一种是固体火箭利用遥测的发动机燃烧室压强和预先建立的内弹道程序估计性能参数[22];另一种是利用视加速度等信息对推力参数进行辨识,从而为后续的控制重构和飞行任务调整提供信息[23]。目前利用火箭视加速度等信息对火箭推力进行辨识的研究中,均假设火箭的质量已知,在这种情况下采用成熟的辨识方法显然可以估计出推力大小。然而火箭实际飞行过程中质量与预设质量是有偏差的,尤其是当火箭发生推力下降的情况下,燃料的消耗速度跟预先设定速度会偏差很大,导致火箭实际质量与预设质量偏差很大,此时质量不能再看作是已知的。而且测量的视加速度信息中运载火箭的推力和质量耦合在一起,直接采用传统的最小二乘或卡尔曼滤波方法都无法区分出推力参数和质量参数,给辨识工作带来了很大的困难。同时推进系统故障将导致推力参数的突变,提高辨识方法的跟踪能力也是一个难点。

精确辨识火箭动力系统故障下的质量和推力,可以更好地进行动力学模型修正、控制重构以及飞行任务调整等。因此,本文针对运载火箭推力与质量辨识相互耦合的问题,利用火箭惯性敏感器件测量的视加速度信息,提出了一种最小二乘法与卡尔曼滤波相结合的联合校正辨识方法,可以精确有效地辨识火箭的推力参数与质量参数。进一步,提出了一种基于视加速度向量2范数差分的推力参数故障检测方法,提高了对推力参数突变的跟踪能力。

1 推力参数与质量参数的辨识问题

本文以两级运载火箭二级飞行段为研究对象[24],假设发动机比冲已知。定义地心惯性坐标系OX1Y1Z1：原点在地心,X1轴在赤道平面内指向发射时刻本初子午线方向,Z1轴垂直赤道平面指向北极,Y1轴满足右手定则。在地心惯性坐标系中建立火箭的二级飞行段动力学方程以及质量消耗方程如下:

(1)

(2)

(3)

本文采用地心惯性坐标系下的视加速度信息作为辨识的观测量,则辨识的观测方程可以表示为:

y=Fu/m+d

(4)

其中,d表示视加速度的测量噪声。

为了后续描述上的方便,记观测矩阵为

h(m,u)=u/m

(5)

则观测方程也可以表示成:

y=h(m,u)F+d

(6)

当发动机出现动力系统故障后,运载火箭的推进剂质量消耗率会随之发生变化,不能再采用质量的预设值;而式(4)或式(6)则表明,测量的视加速度信息中推力与质量耦合在一起的,因此需要将质量和推力同时作为未知量进行辨识,这也是本文研究的辨识问题。

2 联合校正辨识方法

2.1 基本思路

视加速度测量方程式(6)可以看作是对F的观测,然而,由于运载火箭的推力和质量均是未知的、待估计的量,它们是耦合在一起的,因此,无法直接采用最小二乘方法进行推力参数辨识。而且,当运载火箭推力发生故障后,推力参数会发生突变,因此需要采用跟踪性强的辨识方法。

图1 运载火箭推力参数与质量联合校正辨识算法框图

上文的联合校正递推辨识算法的关键在于可以在每个递推步中实现对质量和推力参数的有效更新:在状态(含质量参数)估计模块中,采用的卡尔曼滤波算法包含基于模型的预估和基于测量信息的校正两个环节,因此即使推力参数存在一定偏差,也可以综合测量信息实现对状态(含质量参数)的校正,这正是卡尔曼滤波算法的优点;在推力参数辨识模块,采用最小二乘算法具有很强的抗干扰能力,可以有效降低质量参数偏差带来的影响,同时采用渐消记忆的递推最小二乘算法,通过降低历史数据的权重来实现对时变推力参数的跟踪。

进一步,由于当火箭发生动力系统故障时,推力参数发生突变是重要的特征,以前的观测数据(视加速度)不再反映突变后的推力信息。因此为了提高对推力参数的跟踪性能,在上述联合校正递推辨识算法框架的基础上,又提出了推力突变的检测方法和相应的修正措施。

2.2 联合校正辨识算法的实现

运载火箭推力与质量联合校正递推辨识算法将卡尔曼滤波与渐消记忆的递推最小二乘结合,本节具体介绍该方法的具体实现过程。

选取状态变量x=[rxryrzvxvyvzm]T,则有火箭的质量m=x7。根据式(1)、式(2)以及式(6)将运载火箭动力学模型转化为如下状态空间形式:

(7)

其中:f(x,u,F)和g(x,u,F)表示动力学方程和观测方程中的相应函数,具体表达形式如下:

(9)

卡尔曼滤波算法中对状态的预估方程形式如下:

(10)

状态预估值的协方差矩阵方程为:

P(k+1|k)=Φ(k+1|k)P(k)ΦT(k+1|k)+

(11)

式中:P(k+1|k)表示k+1时刻的状态预估协方差矩阵;Φ(k+1|k)表示状态转移矩阵;P(k)表示k时刻的状态估计协方差矩阵;Γ(k+1|k)表示离散噪声分布矩阵;Qw(k)表示离散过程噪声方差。其中,离散噪声分布矩阵Γ(k+1|k)的计算公式为:

Γ(k+1|k)=Φ(k+1|k)C(k+1)Δt

(12)

式(11)中状态转移矩阵Φ(k+1|k)的计算公式为:

Φ(k+1|k)=eA(k+1)Δt

(13)

式中:A表示函数f(x,u,F)的雅可比矩阵,其具体形式为:

(14)

其中,元素的具体计算公式为:

(15)

卡尔曼滤波增益矩阵K的计算公式为:

K(k+1)=P(k+1|k)HT(k+1)[H(k+1)P(k+1|k)HT(k+1)+Rd(k+1)]-1

(16)

式中:Rd(k+1)表示离散观测噪声方差;H表示函数g(x,u,F)的雅可比矩阵,其具体形式表示为:

(17)

状态估计协方差矩阵的计算公式如下:

P(k+1)=[Is-K(k+1)H(k+1)]P(k+1|k)

(18)

式中:Is表示对应状态维度单位阵。

卡尔曼滤波的状态估计值为:

(19)

(20)

根据k+1时刻运载火箭质量的辨识结果,通过递推最小二乘估计k+1时刻的推力参数。由于推力是变化的,且当发生故障时,推力会产生突变,为了提高参数突变情况下的辨识方法跟踪性能,采用渐消记忆的递推最小二乘算法。引入渐消因子ρ,ρ∈(0, 1],通过降低历史数据的权重进而降低历史数据对当前时刻参数估计的影响。虽然渐消因子越小,辨识的跟踪性越好,但是由于旧数据的减少,往往会降低对噪声的抗干扰能力,因此渐消因子一般选取0.95≤ρ<1。

渐消记忆的递推最小二乘增益矩阵为:

(21)

式中:Ib表示对应观测维度单位阵;PLs表示信息矩阵,且有:

(22)

从而推力参数的估计值为:

(23)

运载火箭推力与质量的联合校正递推辨识算法流程如图2所示。可以看出,每个循环周期内包含质量参数估计模块和推力参数估计模块。两个模块依次交替进行,质量参数估计模块为推力参数估计的最小二乘算法提供观测矩阵,而推力参数估计模块为质量参数估计的卡尔曼滤波算法提供模型参数,二者互相校正,最终实现质量和推力参数的联合辨识。

图2 联合校正辨识方法流程图

2.3 动力系统故障的检测

在联合校正辨识算法中,推力参数辨识采用渐消记忆的递推最小二乘算法,通过不断相乘渐消因子降低对历史数据的权重来实现对推力参数变化的跟踪。因此当发生推力突变后,渐消因子需要一定的时间才能比较彻底的消除突变前测量数据的影响,这就导致也需要一定的时间才能跟踪上突变后的推力。当然,通过降低渐消因子也可以提高跟踪能力,但这同时会降低对噪声的抗干扰能力。因此,如果能够检测出动力系统故障导致的推力突变的时刻,那么就可以直接删除掉突变前的测量数据,从突变时刻重新开始最小二乘辨识,从而大大提高对推力参数突变的跟踪能力。

本文的观测信息是地心惯性坐标系下的视加速度,根据式(4),视加速度向量的2范数与推力成正比,当运载火箭推力发生突变时,由于运载火箭质量不会发生突变,因此推力的突变必然导致视加速度的突变。

理论上,发生突变时视加速度的变化率趋近无穷大。然而工程上只能利用采样数据的差分来近似变化率,所以提出通过计算视加速度向量2范数的差分来判断推力参数突变时刻。当视加速度向量2范数的差分在某一时刻的绝对值超过设定阈值时,认为在该时刻运载火箭的推力发生了突变。因此,得到判断推力参数突变的条件如下:

(24)

不等式左边表示视加速度向量2范数差分的绝对值,不等式右边中η为给定变化率的阈值。当不等式成立时,认为推力参数在k+1时刻发生突变。在推力没有突变的时间内,视加速度向量2范数差分的绝对值非常小,约10-3量级;考虑到实际测量的信号中含有噪声,因此η的设置需要比10-3量级大一些,不妨设置为η=10-1。

实际上,视加速度变化率阈值η的设定可以根据发动机正常工作下推力的变化率和视加速度采样周期进行预估。同时,在使用式(24)进行判断时,如果仅采用某一个时刻的计算值进行判别,可靠性过低,容易引起误判。因此一方面可将测量信号先通过低通滤波器、滤除测量噪声的影响;另一方面,利用卡尔曼滤波算法本身具有的滤噪特性,也可以将视加速度测量值替换为对应的视加速度预估值进行判断,也能在一定程度上提高检测可靠性。

检测出突变时刻后,在推力参数辨识的最小二乘递推算法中就可以删除突变前的测量数据,从而消除突变前测量数据对突变后推力参数辨识的影响。

有故障检测的运载火箭质量和推力参数联合校正递推辨识算法如框图3所示。与前面无故障检测的算法框图1相比,主要增加了故障检测模块和对推力参数辨识模块的修正。具体修正如下。

图3 有故障检测的运载火箭推力参数与质量联合校正辨识算法框图

u(k+1)))-1

(25)

并利用k+1时刻测量信息直接进行k+1时刻推力参数的最小二乘估计,如下:

u(k+1))y(k+1)

(26)

3 仿真校验

本节以火箭二级飞行段出现推力下降故障为研究对象,验证本文所提出运载火箭推力参数和质量参数联合校正递推辨识算法的有效性。

火箭二级飞行段仿真时间选取239～249 s,采样周期为Δt=0.01 s,假设火箭在243.9 s发生推力突变下降故障。动力学方程式(1)和(2)中的引力常数为μ=3.986 0×1014,火箭比冲为Isp=342.786 4 s,海平面处重力加速度为g0=9.806 7 m/s2。过程噪声方差选取Qw=diag(0.01,0.01,0.01,0.01,0.01,0.01,0.01)。

通过在视加速度仿真值上施加高斯白噪声来模拟含观测噪声的视加速度测量值[17]。为了验证本文方法对噪声不确定性的适应性,适当增大噪声的方差,噪声的方差取Rd=diag(0.05,0.05,0.05)。

图4(a)、图4(b)和图4(c)给出了推力下降20%故障下的视加速度仿真模拟数据曲线,可以看出,当火箭推力发生突变后,视加速度也随之突变。

图4 推力下降20%故障下视加速度模拟值

-1.6×103,3.4×103,-1.7×103,3.8×104]T

(27)

且有

(28)

3.1 联合校正辨识算法中的参数影响

本节验证联合校正辨识算法中渐消因子的大小对辨识结果的影响,以及验证有故障检测的效果。在无故障检测的情况下,渐消因子分别选择0.99,0.97与0.95;在有故障检测情况下,渐消因子选择0.97。则采用本文联合校正辨识算法的运载火箭推力与推力误差辨识结果如图5、图6所示。

图5 不同渐消因子与有无故障检测的推力参数辨识结果

图6 不同渐消因子与有无故障检测的推力参数辨识误差

从图5、图6中的不同渐消因子-无故障检测曲线可以看出,渐消因子选取的越小,在推力发生故障突变后,辨识的跟踪性越好;在辨识结果跟踪上理论推力的这些时间内,推力辨识误差小于1%,渐消因子选取越小,推力辨识误差越大。因此,渐消因子选取越小虽然可以有效提高辨识的跟踪性,但是会增加推力辨识误差,渐消因子选取不宜过大或过小。从图中渐消因子均取0.97时,有故障检测与无故障检测的两条曲线可以看出,在推力发生故障突变后,有故障检测的推力辨识结果可以立刻跟踪上理论的推力;在243.9 s至246.4 s时间内,有故障检测的推力辨识误差与无故障检测的推力辨识误差几乎一致,且有故障检测的结果在推力辨识结果跟踪突变后推力的这段时间里推力误差也只有2.5%。

质量与质量误差的辨识结果如图7、图8所示。从图7、图8可以看出,渐消因子越小,在故障下推力突变后,质量跟踪越好;在渐消因子均取0.97,有故障检测和无故障检测的结果可以看出,在故障下推力突变后,有故障检测的质量误差更低,小于0.002 5%。这是因为在推力突变后,减小渐消因子以及有故障检测的推力参数估计值较快得接近理论值,进而在联合校正辨识算法中通过推力估计结果来辨识的质量结果更精确。

图7 不同渐消因子与有无故障检测的质量辨识结果

图8 不同渐消因子与有无故障检测的质量辨识误差

3.2 联合校正辨识方法的有效性

本节验证联合校正辨识方法的有效性,包括不同推力初值对辨识结果的影响以及不同推力下降故障下的辨识结果。本节采用渐消因子取0.97、有故障检测的方法,推力初值在理论值的基础上分别取0.9倍理论值和1.1倍理论值,则有不同推力初值下的运载火箭推力与质量辨识结果如图9、图10所示。

图9 不同推力初值下的推力辨识结果

图10 不同推力初值下的质量辨识结果

从图9中可以看出,当推力初值有偏差时,联合校正辨识方法可以较快地跟踪上理论的推力;从图10中可以看出,当推力初值有偏差时,联合校正辨识方法对质量的辨识效果依然理想。

当火箭发生推力突变下降故障后,往往在故障发生后的一段时间并不是保持不变,而是变化的,即二次故障。因此,本节验证本文的有故障检测的联合校正辨识方法对不同故障大小以及二次故障下的适应性。通过假设在推力发生突变下降后经过一段正弦波动变化后保持在某个数值来模拟二次故障,推力初值取理论值。

则有在不同故障大小以及二次故障情况下的辨识结果如图11至图14所示。

图11 不同故障下的推力辨识结果

从图11、图12中可以看出,当发生不同大小故障或二次故障时,有故障检测的联合校正辨识仍然可以快速跟踪上变化的推力,在验证不同大小故障或二次故障时,需获得不同的观测值,由于噪声是随机的高斯白噪声,因此即便在发生故障时间前,推力辨识误差仍不相同。在不同推力大小故障下,除了在故障时刻附近推力辨识误差有时会增加达到2%左右,其它时间内均小于0.6%。