CN114218750A

CN114218750A - 基于数字滤波器的星载微推力器推力响应时间测量方法

Info

Publication number: CN114218750A
Application number: CN202111352455.6A
Authority: CN
Inventors: 洪延姬; 冯孝辉; 金星; 叶继飞; 崔海超; 冯高平
Original assignee: Peoples Liberation Army Strategic Support Force Aerospace Engineering University
Current assignee: Peoples Liberation Army Strategic Support Force Aerospace Engineering University
Priority date: 2021-11-16
Filing date: 2021-11-16
Publication date: 2022-03-22
Anticipated expiration: 2041-11-16
Also published as: US20230150696A1; CN114218750B

Abstract

本发明属于航天卫星推进技术领域，具体涉及一种基于数字滤波器的星载微推力器推力响应时间测量方法。本发明的推力响应时间测量方法包括步骤，S1：对扭摆推力测量系统的非零初始条件进行归零处理,得到变量代换后的推力测量系统的振动微分方程；S2：在变量代换后的推力测量系统后串联一个数字滤波器,得到等效灵敏度高频推力测量系统；S3：确定等效灵敏度高频推力测量系统的系统响应；S4：从系统响应判读星载微推力器的推力响应时间；S5：利用系统响应反演计算待测推力，进一步确认推力响应时间。本发明的等效灵敏度高频推力测量系统的固有频率大幅提高并极大缩短了进入稳态时间，实现了快变推力响应时间的测量。

Description

基于数字滤波器的星载微推力器推力响应时间测量方法

技术领域

本发明涉及航天卫星推进技术领域，特别是一种基于数字滤波器的星载微推力器推力响应时间测量方法。

背景技术

在卫星精确控制场合，对微推力器提出了快速调节推力要求，需要将推力响应时间控制在指定范围内。推力响应时间可用推力上升时间和推力下降时间表示，当推力从较小稳定推力变化到较大稳定推力时，推力响应时间采用推力上升时间表示；当推力从较大稳定推力变化到较小稳定推力时，推力响应时间采用推力下降时间表示，此时，推力响应时间是微推力器的考核指标之一。

推力响应时间的测量难点在于：推力变化频率远大于推力测量系统的固有频率，或推力测量系统的固有频率太低，不满足快变推力响应时间测量要求。首先，由于用于微推力测量的推力测量系统的高灵敏度要求，推力测量系统的刚度系数很小，使得固有频率降低；其次，由于推力器搭载在推力测量系统上，推力测量系统的转动部件(推力器、配重、结构梁等)的转动惯量很大，也使得固有频率降低；最后，如果测量10ms级推力响应时间，则要求推力测量系统的固有线频率为100Hz以上，这是目前任何(扭摆、吊摆、倒摆、弯振)推力测量系统，在物理结构上都难以实现的。

现有的星载微推力器推力响应时间测量方法主要缺陷如下：

(1)不了解推力响应时间与推力测量系统的系统响应之间内在联系。实际上，推力测量系统的系统响应，包括暂态系统响应和稳态系统响应，只有推力响应时间大于推力测量系统的暂态持续时间(或进入稳态时间)，才能从推力测量系统的系统响应反映出推力响应时间。

(2)针对快变推力响应时间的动态推力计算方法不成熟，难以反映准确的推力上升过程或推力下降过程。在动态推力计算中，以差商作为导数的近似值反向求解振动微分方程(振动方程反问题求解)，一是误差大：差商作为导数的近似值，一阶导数误差很大，二阶导数误差更大；二是收敛性差：步长较大时截断误差增大，步长较小时舍入误差增大。

发明内容

有鉴于此，本发明的目的在于解决快变推力响应时间的测量问题。

为解决上述技术问题，本发明提供了一种基于数字滤波器的星载微推力器推力响应时间测量方法，包括以下步骤：

S1：对扭摆推力测量系统的非零初始条件进行归零处理,得到变量代换后的推力测量系统的振动微分方程；

S2：在变量代换后的推力测量系统后串联一个数字滤波器,得到等效灵敏度高频推力测量系统；

S3：确定等效灵敏度高频推力测量系统的系统响应；

S4：从系统响应判读星载微推力器的推力响应时间；

S5：利用等效灵敏度高频推力测量系统反演计算待测推力，进一步确认推力响应时间。

步骤S1中，在待测推力f′(t)作用下系统响应θ′(t)由位移传感器测量得到，根据非零初始条件

进行变量代换，令

得到变量代换后的推力测量系统的振动微分方程为：

f(t)＝f′(t)+f₀(t)

其中，f₀(t)是与非零初始条件相关的等效推力，f′(t)为待测推力，θ′(t)为该推力作用下由位移传感器测量得到的系统响应，ζ为阻尼比，ω_n为固有振动频率，J为转动惯量，L_f为力臂；θ₀表示初始扭转角，θ₀和

为常数，θ(t)为变量代换后系统响应，满足初始条件θ(0)＝0和

步骤S2中，数字滤波器的传递函数为：

数字滤波器的单位脉冲响应函数为

其中，ω_n为固有振动频率，ζ为阻尼比，ω_n1为等效灵敏度高频推力测量系统的固有频率，ω_n1＝C_ωω_n，C_ω为频率比，ζ₁为等效灵敏度高频推力测量系统的阻尼比，

为改进后振动频率。其中，C_ω＞1，ζ₁＝0.7～0.9。

步骤S3中，数字滤波器的输入为θ(τ)，通过位移传感器测量值θ′(τ)计算得到θ(τ)，数字滤波器输出的系统响应θ₁(t)为：

其中,τ为积分变量。

步骤S4中，当系统响应与稳态系统响应比较相差1％时所对应的时间为推力响应时间。

推力响应时间包括推力上升时间和推力下降时间。

步骤S5中，等效灵敏度高频推力测量系统的单位脉冲响应函数h₁(t)为

已知数字滤波器的输出为θ₁(t)，利用等效灵敏度高频推力测量系统的推力积分方程：

计算名义推力f(τ)。

待测推力f′(t)的估计值F′(t)为：F′(t)＝F(t)-f₀(t),其中，

F(t)是名义推力f(t)的估计值，ζ为阻尼比，ω_n为固有振动频率，J为转动惯量，L_f为力臂，θ₀表示初始扭转角，θ₀和

为常数；当系统响应与稳态系统响应比较相差1％时所对应的时间为推力响应时间。

有益效果

本发明通过在固有频率低的推力测量系统上串联数字滤波器，构建了等效灵敏度高频推力测量系统，采用数字滤波器的设计方法，使得等效灵敏度高频推力测量系统的固有频率大幅提高并极大缩短了进入稳态时间，实现了快变推力响应时间的测量，并且等效灵敏度高频推力测量系统与原推力测量系统的灵敏度相等，实现了等灵敏度测量。通过数字滤波器的设计，可调节振动频率和阻尼比，根据推力上升时间和推力下降时间的测量需求，通过调节振动频率和阻尼比，实现快变推力上升时间和推力下降时间的测量。

附图说明

图1是本发明推力响应时间测量方法的流程图。

图2是本发明等效灵敏度高频推力测量系统示意图。

图3是推力上升时间的变化图。

图4是原推力测量系统的系统响应测量值图。

图5是变量代换后系统响应的变化图。

图6是数字滤波器的系统响应变化图。

图7是推力上升时间25ms下反演计算的推力和误差图。

图8是推力上升时间50ms下反演计算的推力和误差图。

图9是推力上升时间100ms下反演计算的推力和误差图。

图10是推力逐渐下降时数字滤波器的系统响应变化图。

图11是推力逐渐下降时反演计算的推力图。

图12是推力逐渐下降时反演计算的推力误差图。

具体实施方式

下面将参考附图对本发明的具体实施方式进行详细说明。

如图1所示，是本发明推力响应时间测量方法的流程图，下面对实施过程进行详细说明。

1.建立推力响应时间与推力测量系统的进入稳态时间关系

建立推力响应时间与推力测量系统的进入稳态时间关系，是推力测量系统能够测量给定推力响应时间的必要条件。

1)推力响应时间与振动微分方程初始条件

微推力器的推力响应时间不但需要考虑推力上升时间还要考虑推力下降时间。推力由较小稳定推力上升到较大稳定推力时，推力响应时间可用推力上升时间表示；推力由较大稳定推力下降到较小稳定推力时，推力响应时间可用推力下降时间表示。研究推力响应时间必须考虑初始条件的影响。

扭摆推力测量系统的振动微分方程为

式中，f′(t)为待测推力，θ′(t)为该推力作用下系统响应(由位移传感器测量得到)，ζ为阻尼比，ω_n为固有振动频率，J为转动惯量，L_f为力臂。

当推力从一个稳定推力变化到另一个稳定推力时，初始条件不为零，为

式中，若是从一个非零稳定推力开始变化，可近似认为

若是从零推力开始变化，可近似认为θ′(0)≈0和

显然，推力测量系统的振动微分方程的初始条件，反映了推力是从怎样初始推力开始变化的，例如，从非零稳定推力开始变化，还是从零推力开始变化。因此，测量推力响应时间必须明确推力测量系统的振动微分方程的初始条件。

2)建立推力响应时间与进入稳态时间关系

扭摆推力测量系统的振动微分方程为

如果初始条件不为零，为θ′(0)＝θ₀和

设推力的拉普拉斯变换为L[f′(t)]＝F′(s)，系统响应及导数的拉普拉斯变换为

L[θ′(t)]＝Θ′(s) (4)

振动方程两边取拉普拉斯变换，整理得

两边取拉普拉斯反变换，可得

令

θ′(t)＝θ′₁(t)+θ′₂(t) (9)

在θ₁′(t)中包括暂态系统响应和稳态系统响应。例如，在θ₁′(t)中代入周期推力f′(t)＝A_f cosωt(A_f为常数表示周期推力的振幅，ω为周期推力的频率)，推力测量系统的系统响应，为

其中

式中，S＝L_f/k为推力测量系统的灵敏度，

为扭转刚度系数，R_ω＝ω/ω_n为频率比。显然，在θ₁′(t)中系统响应的第一项为暂态系统响应，代表暂态过程随着时间变化，暂态过程的振动频率为推力测量系统的振动频率；第二项为稳态系统响应，代表稳态过程随着时间变化，稳态过程的振动频率为周期推力的振动频率。θ₂′(t)也为暂态系统响应，其振动频率为推力测量系统的振动频率。

总之，暂态系统响应的振幅是逐渐衰减的，衰减项为

设衰减项

时所对应时间(暂态振幅衰减到1％所对应时间)，定义为暂态持续时间或进入稳态时间t_s，有

式中，T_n＝2π/ω_n为推力测量系统的固有周期。

推力测量系统的进入稳态时间或暂态持续时间，是推力测量系统的固有属性，具有以下结论：

(1)推力测量系统的固有频率越大则周期越小，进入稳态时间越短；

(2)推力测量系统的阻尼比越小，进入稳态时间越大；

(3)推力响应时间只有大于推力测量系统的进入稳态时间条件下，才能从系统响应曲线表现出推力响应时间的变化特点，或才能从系统响应曲线判读推力响应时间。

例如，在常用阻尼比0.1≤ζ≤0.3范围内，推力测量系统的进入稳态时间为固有周期的3～7倍，设推力上升时间为t_s，推力测量系统的固有周期为T_n，则有

t_s≥7T_n (16)

如果要求推力上升时间t_s≤50ms，则要求推力测量系统的固有周期为T_n≤7.1ms，也就是线频率为f_n≥140Hz，这是任何结构的推力测量系统都难以做到的。说明测量快变推力响应时间，仅考虑从物理结构设计上提高推力测量系统的固有频率是远远不够的，需要另辟蹊径。

2.非零初始条件的归零处理方法

在推力响应时间测量中初始条件不为零，但是，在后续数字滤波器设计和构建等效灵敏度高频推力测量系统中，由于利用传递函数建立模型，必须将非零初始条件进行归零处理。

扭摆推力测量系统的振动微分方程为

初始条件不为零，为

式中，θ₀和

为常数。

为了保证初始条件为零，进行变量代换，令

代入振动方程，整理得

代入初始条件，可得

将变量代换后的推力测量系统的振动微分方程，改写为

f(t)＝f′(t)+f₀(t) (25)

式中，f₀(t)是与非零初始条件相关的等效推力。

当推力从一个稳定推力变化到另一个稳定推力时，有θ′(0)＝θ₀和

可得

θ(t)＝θ(t)-θ₀ (27)

说明：代换后系统响应θ(t)与代换前系统响应θ′(t)，暂态系统响应和稳态系统响应都相同只是系统响应曲线上下移动一定距离。

非零初始条件的归零处理方法为：

(1)在待测推力f′(t)作用下系统响应θ′(t)是由位移传感器测量得到，根据非零初始条件

进行变量代换，为

式中，θ(t)为变量代换后系统响应，满足初始条件为零θ(0)＝0和

条件。

(2)当推力从一个稳定推力变化到另一个稳定推力时(θ′(0)＝θ₀和

)，变量代换后系统响应θ(t)与代换前系统响应θ′(t)的进入稳态时间相同，因此，从代换后系统响应θ(t)曲线可判读推力响应时间。

(3)根据变量代换后系统响应θ(t)，利用变量代换后振动微分方程，可计算名义推力f(t)，从而求得待测推力f′(t)为

f′(t)＝f(t)-f₀(t) (29)

式中，f₀(t)是与非零初始条件相关的等效推力。

3.利用数字滤波器构建等效灵敏度高频推力测量系统

为了解决推力测量系统固有频率过低，不能满足快变推力响应时间测量的难题，采用数字滤波器构建等效灵敏度的高频推力测量系统，大幅提高固有频率，满足快变推力响应测量的需求。

如图2所示，是本发明等效灵敏度高频推力测量系统示意图。采用数字滤波器构建等效灵敏度高频推力测量系统的方法为：

(1)在变量代换后的推力测量系统的传递函数上，串联一个数字滤波器，构造等效灵敏度高频推力测量系统；

(2)所构造的等效灵敏度高频推力测量系统，灵敏度与变量代换前的推力测量系统相同，但是，通过数字滤波器的设计，可调节振动频率和阻尼比；

(3)根据推力上升时间和推力下降时间的测量需求，通过调节振动频率和阻尼比，实现推力上升时间和推力下降时间的测量。

4.数字滤波器设计方法与系统响应计算方法

1)数字滤波器设计方法

设变量代换后的推力测量系统的输入为名义推力f(t)(拉普拉斯变化为F(s))，其系统响应为θ(t)(由位移传感器测量值θ′(t)计算得到，拉普拉斯变化为Θ(s))；数字滤波器的传递函数为D(s)，数字滤波器的输入为θ(t)，数字滤波器的输出为θ₁(t)(由θ(t)计算得到，拉普拉斯变化为Θ₁(s))；所构造的等效灵敏度高频推力测量系统的传递函数为H₁(s)＝H(s)D(s)，等效灵敏度高频推力测量系统的输入为名义推力f(t)，输出为θ₁(t)。

扭摆推力测量系统的传递函数为

式中，L_f为力臂，J为转动惯量，ζ为阻尼比，ω_n为固有振动频率，

为振动频率，

为扭转刚度系数。

数字滤波器设计方法是采用数字滤波器与变量代换后的推力测量系统的输出端串联，构成改进的扭摆推力测量系统。改进后扭摆推力测量系统的传递函数为

式中，改进后扭摆推力测量系统的固有振动频率为ω_n1，阻尼比为ζ₁，扭转刚度系数为k₁。

设数字滤波器的传递函数为D(s)，则有

为了减小进入稳态时间提高推力响应时间的测量能力，需要增大固有振动频率，令ω_n1＝C_ωω_n(C_ω＞1)和k₁＝k(灵敏度相等)，有

数字滤波器的传递函数为

式中，

为改进后振动频率。进一步简化，可得数字滤波器的传递函数为

设数字滤波器的单位脉冲响应函数为d(t)，则有

可得数字滤波器的单位脉冲响应函数为

数字滤波器的输入为θ(τ)(通过位移传感器测量值θ′(τ)计算得到)，设数字滤波器的输出为θ₁(t)，则有

从而，获得等效灵敏度高频推力测量系统的系统响应。

数字滤波器的设计方法为：

(1)以灵敏度相等和增大固有振动频率为目标，设计数字滤波器的传递函数；

(2)根据数字滤波器的传递函数，通过拉普拉斯反变换获得数字滤波器的单位脉冲响应函数；

(3)利用变量代换后的推力测量系统的系统响应，计算数字滤波器的系统响应。

2)数字滤波器的系统响应计算方法

数字滤波器的输出量需要采用其单位脉冲响应函数，利用数值积分方法进行计算。

设数字滤波器的单位脉冲响应函数为d(t)，有

令

d(t)＝d₁(t)+d₂(t) (42)

代入以下方程

可得

式中，

可根据系统响应θ(t)直接计算，

的计算需要采用数值积分方法。

为了计算

的积分方程，设时间采样步长为h，t_i＝ih(i＝0,1,2,…)，有

由于[t₀＝0,t_i＝ih]划分为i等分，对于τ_j＝jh(j＝0,1,2,…,i)，令函数

g(t_i,τ_j)＝θ(τ_j)d₂(t_i-τ_j) (50)

的计算方法为采用梯形积分公式起步计算，再交替采用3点和4点辛普森积分公式计算，具体为：

(1)当子区间数目n＝1时，采用梯形积分公式起步计算。

(2)当子区间数目n＝2时，采用3点(2个区间)辛普森积分公式计算。

(3)当子区间数目n＝3时，采用4点(3个区间)辛普森积分公式计算。

(4)当子区间数目n＝2k(k＝2,3,4,…)时，采用3点辛普森积分公式计算。

用2k+1个节点t_i＝ih(i＝0,1,2,…,2k)，将积分区间[t₀,t_2k]划分为2k等分，子区间长度为h＝t_2k/(2k)，每个子区间上反复使用3点辛普森积分公式，可得复化辛普森公式为

(5)当子区间数目n＝2k+1(k＝2,3,4,…)时，交替采用3点和4点辛普森积分公式计算。

用2k+2个节点t_i＝ih(i＝0,1,2,…,2k+1)，将积分区间[t₀,t_2k+1]划分为2k+1等分，子区间长度为h＝t_2k+1/(2k+1)，在前面的2k-2(k＝2,3,4,…)个子区间反复使用3点辛普森积分公式，在后面2k-2、2k-1、2k、2k+1等4个点使用4点辛普森积分公式，可得复化辛普森公式为

在数字滤波器的系统响应计算中，变量代换前的推力测量系统的固有频率ω_n和阻尼比ζ是已知的，根据推力响应时间测量需要，选取等效灵敏度高频推力测量系统的固有频率ω_n1＝C_ωω_n(通过选取频率比C_ω＞1实现)，选取阻尼比ζ₁＝0.7～0.9进一步达到减小进入稳态时间的目的。

5.利用等效灵敏度高频推力测量系统反演计算推力方法

由于等效灵敏度高频推力测量系统的固有频率大幅提高，可实现快变推力响应时间的测量。

一方面，从数字滤波器的系统响应，也就是等效灵敏度高频推力测量系统的系统响应，可间接判读推力上升时间和推力下降时间；另一方面，根据数字滤波器的系统响应(也就是等效灵敏度高频推力测量系统的系统响应)，利用等效灵敏度高频推力测量系统的传递函数和单位脉冲响应函数，反演计算推力，可以进一步确认推力上升时间和推力下降时间。此时，推力响应时间的判读更为合理和可靠。

1)建立等效灵敏度高频推力测量系统的推力积分方程

等效灵敏度高频推力测量系统的传递函数为

等效灵敏度高频推力测量系统的单位脉冲响应函数h₁(t)为

已知数字滤波器的输出为θ₁(t)，设名义推力f(τ)，等效灵敏度高频推力测量系统的推力积分方程为

上述方程称为推力积分方程，根据推力积分方程，利用θ₁(t)可反演计算名义推力f(τ)。

2)推力积分方程离散化计算推力

将等效灵敏度高频推力测量系统的单位脉冲响应函数h₁(t)代入推力积分方程，可得

已知等效灵敏度高频推力测量系统的系统响应为[t_i,θ₁(t_i)](i＝0,1,2,3,…)，t₀＝0时θ₁(t₀)＝0，采样时间步长为h，t_i＝ih，令θ₁(t_i)＝θ_1i和F(t_i)＝F_i。对于任意给定t_i＝ih(i＝1，2,3…)，推力积分方程组为

式中，用F(τ)代替f(τ)是表示F(τ)为推力f(τ)的估计值。

将区间[0,t_i]划分为i等分(i＝1,2,3,…)，节点为τ_j＝jh(j＝0,1,…,i)，被积分函数为

满足

采用复化辛普森离散化方法，将推力积分方程组离散化为推力线性方程组，需要交替使用3点和4点辛普森积分公式，分以下5种情况讨论。

(1)当i＝1时，根据梯形积分公式起步计算，可得

可得

(2)当i＝2时，根据3点辛普森积分公式，可知

可得

a₂₀F₀+a₂₁F₁＝C_fθ₁₂ (66)

(3)当i＝3时，根据4点辛普森积分公式，可知

可得

a₃₀F₀+a₃₁F₁+a₃₂F₂＝C_fθ₁₃ (68)

(4)当i＝2k(k＝2,3,…)时，根据3点辛普森积分公式，可得

可得

(5)当i＝2k+1(k＝2,3,…)时，在2k+2个点中，前面的i＝0,1,…,2k-3,2k-2等点采用3点辛普森积分公式；最后的4个点i＝2k-2,2k-1,2k,2k+1等采用4点辛普森积分公式，可得

即有

(6)计算得到名义推力f(τ)的估计值F(τ)，设待测推力f′(τ)的估计值为F′(τ)，可得待测推力估计值为

F′(t)＝F(t)-f₀(t) (80)

如果在步骤(1)至步骤(5)中，线性方程组表示为矩阵形式，则有

推力的复化辛普森计算方法，需要交替使用3点和4点辛普森公式，计算方法构造复杂，但是由于复化辛普森公式的截断误差与时间步长h⁴成正比，计算精度大幅提高。

如图1所示，本发明基于数字滤波器的星载微推力器推力响应时间测量方法，包括以下步骤：

S3：确定等效灵敏度高频推力测量系统的系统响应；

S4：从系统响应判读星载微推力器的推力响应时间。

S5：利用等效灵敏度高频推力测量系统反演计算待测推力，可以进一步确认推力响应时间。

下面结合一个具体的实施例进行详细说明。

实施例

已知扭摆推力测量系统的扭转刚度系数为k＝0.075N·m/rad，振动频率为ω_d＝0.25rad/s，阻尼比为ζ＝0.2，固有频率为ω_n＝0.255155rad/s(固有线频率f_n＝0.040609Hz，固有周期为24.6s)，转动惯量为J＝1.152002kg·m²，力臂为L_f＝0.4m。

在扭摆转动部件中，搭载3kg推力头和3kg配重其转动惯量为2×3×0.4²＝0.96kg·m²，再加上阻尼器、标定力装置、横梁等，转动部件的转动惯量为J＝1.152002kg·m²。

采用该推力测量系统，利用所提出的推力响应时间的数字滤波器测量方法，测量推力响应时间分别为25ms、50ms、100ms的情况。

初始条件为：θ₀＝50μrad，

由于灵敏度为S＝L_f/k，初始条件为θ₀＝50μrad和

对应的推力为：

这表示推力从9.375μN开始变化。

1.建立推力响应时间与推力测量系统的进入稳态时间关系

该推力测量系统进入稳态时间为

因此，该推力测量系统只能测量推力响应时间大于89.8s的情况，远不能满足推力响应时间为25ms、50ms、100ms的测量要求。

2.非零初始条件的归零处理方法

在待测推力f′(t)作用下系统响应θ′(t)是由位移传感器测量得到。此处，为了检验推力响应时间的数字滤波器测量方法，取推力函数形式为：

f′(t)＝20[1-exp(-t/τ)](单位：μN)

式中，τ为时间常数表征推力上升时间大小，推力上升时间为t_f≈5τ。由f′(t)根据变量代换前振动微分方程，计算出θ′(t)作为位移传感器的测量值。

为了研究推力响应时间分别为25ms、50ms、100ms的情况，推力的时间常数分别取τ＝5ms、10ms、20ms。如图3所示，①、②、③分别为推力响应时间为25ms、50ms、100ms时，推力上升时间的变化(推力从9.375μN开始变化到29.375μN，推力增量为f′(t))。

如图4所示，①、②、③分别为推力响应时间为25ms、50ms、100ms时，原推力测量系统的系统响应测量值。

根据初始条件θ₀＝50μrad和

进行变量代换，变量代换后系统响应为：

如图5所示，①、②、③分别为推力响应时间为25ms、50ms、100ms时，变量代换后系统响应的变化。

变量代换后，推力测量系统的振动微分方程为

其中，f(t)＝f′(t)+f₀(t)

式中，f₀(t)是与非零初始条件相关的等效推力。

3.利用数字滤波器构建等效灵敏度高频推力测量系统

在变量代换后，推力测量系统的传递函数上，串联一个数字滤波器，构造等效灵敏度高频推力测量系统。

如图2所示，在变量代换后，推力测量系统的传递函数为H(s)，其输入为待求名义推力f(t)，其输出为计算得到的θ(t)；所串联的数字滤波器的传递函数为D(s)，数字滤波器的输入为θ(t)，数字滤波器的输出为θ₁(t)(由θ(t)计算得到)；所构造的等效灵敏度高频推力测量系统的传递函数为H₁(s)＝H(s)D(s)，等效灵敏度高频推力测量系统的输入为名义推力f(t)，输出为θ₁(t)。

4.数字滤波器设计方法与系统响应计算方法

为了减小进入稳态时间提高推力响应时间的测量能力，需要增大固有振动频率，令ω_n1＝C_ωω_n(C_ω＞1)和k₁＝k(灵敏度相等)，数字滤波器的输出为θ₁(t)，则有：

式中，d(t)为数字滤波器的单位脉冲响应函数，由数字滤波器设计得到。

采用数字滤波器的系统响应计算方法，计算得到等效灵敏度高频推力测量系统的系统响应θ₁(t)。

选取等效灵敏度高频推力测量系统的频率比为C_ω＝1000，则固有频率ω_n1＝255.155rad/s(固有频率由变量代换前的0.255155rad/s，提升为变量代换后的255.155rad/s，即固有周期变为24.6ms)，并且选取阻尼比ζ₁＝0.7，变量代换后推力测量系统进入稳态时间为：

如图6所示，①、②、③分别为数字滤波器的系统响应(由图5数据计算得到)。由图6可判读所施加推力的推力响应时间分别为25ms、50ms、100ms，并且推力响应时间25ms是刚好能够判读的情况。

5.利用等效灵敏度高频推力测量系统反演计算推力方法

根据推力积分方程，利用θ₁(t)可反演计算名义推力f(τ)。

根据所构造的推力积分方程，采用数值积分离散化方法建立推力线性方程组，求解名义推力f(τ)的估计值F(τ)，设待测推力f′(τ)的估计值为F′(τ)，可得待测推力估计值为

F′(t)＝F(t)-f₀(t)

待测推力f′(τ)估计值F′(τ)的误差为

式中，推力函数形式为

f′(t)＝20[1-exp(-t/τ)](单位：μN)

式中，推力的时间常数分布取τ＝5ms、10ms、20ms。

如图7所示，为推力上升时间25ms下反演计算推力和误差。从图7可清晰确认推力时间响应为25ms，并且20ms以后推力误差小于0.5％。

如图8所示，为推力上升时间50ms下反演计算推力和误差。从图8可清晰确认推力时间响应为50ms，并且20ms以后推力误差小于0.5％。

如图9所示，为推力上升时间100ms下反演计算推力和误差。从图9可清晰确认推力时间响应为100ms，并且20ms以后推力误差小于0.5％。

6.推力下降时间测量的情况

下面举例说明在推力逐渐下降时，推力下降时间测量情况。

取推力函数形式为：f′(t)＝-20[1-exp(-t/τ)](单位：μN)

式中，当推力的时间常数取τ＝5ms、10ms、20ms时，推力下降时间分别为τ_f＝25ms、50ms、100ms。

将上述推力代入以下扭摆推力测量系统的振动微分方程

初始条件为：θ₀＝150μrad，

式中，阻尼比为ζ＝0.2，固有频率为ω_n＝0.255155rad/s，转动惯量为J＝1.152002kg·m²，力臂为L_f＝0.4m。采用常微分方程数值计算方法，求得θ′(t)作为系统响应的测量值。

初始条件θ₀＝150μrad和

表示推力从28.125μN开始逐渐下降，减小的增量为20μN(推力从28.125μN开始减小到8.125μN)。

如图10所示，为数字滤波器的系统响应变化(由θ′(t)计算得到)。从图10可判读推力下降时间分别为25ms、50ms、100ms，并且推力下降时间25ms是刚好能够判读的情况。

如图11所示，为反演计算的推力(由θ₁(t)计算得到)。从图11可清晰确认推力下降时间分别为25ms、50ms、100ms。

如图12所示，为反演计算的推力误差。在上述三种情况下，在20ms以后，推力误差都小于1％。

综上，本发明的主要技术特征和优点为：

(1)针对目前不了解推力响应时间与推力测量系统的系统响应之间内在联系的现状，本发明建立了推力响应时间与推力测量系统的进入稳态时间关系。首先，指明了推力响应时间测量中，推力从一个稳定推力开始变化，意味着推力测量系统的初始条件不为零；其次，推力测量系统的进入稳态时间是推力测量系统的固有属性，推力测量系统的固有频率越大则进入稳态时间越短，推力响应时间只有大于推力测量系统的进入稳态时间条件下，才能从系统响应曲线判读推力响应时间；最后，指明了10ms级推力响应时间测量，要求推力测量系统的固有频率为100Hz以上，这是目前任何(扭摆、吊摆、倒摆、弯振)推力测量系统在物理结构上都难以实现的，因此，需要另辟蹊径。

(2)针对目前快变推力响应时间的动态推力计算方法不成熟和不准确的现状，本发明通过固有频率低的推力测量系统上串联数字滤波器，构建等效灵敏度高频推力测量系统，采用数字滤波器设计方法，使得等效灵敏度高频推力测量系统的固有频率大幅提高并极大缩短了进入稳态时间，实现了快变推力响应时间的测量。首先，根据等效灵敏度高频推力测量系统的系统响应曲线判读推力响应时间；其次，通过推力反演计算由推力曲线确认推力响应时间，使得推力响应时间评估结论更加准确和可靠。

(3)星载微推力器推力响应时间的数字滤波器测量方法，所构建的等效灵敏度高频推力测量系统的固有频率提升500-1500倍，阻尼比调节范围由0.1-0.3变为0.7-0.9，大幅减小了进入稳态时间，实现了10ms级推力响应时间的测量。

以上所述仅为本发明的具体实施方式而已，并不用来限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。