CN112231890A - 基于扭摆测量系统的高平稳电推力器的推力测评方法 - Google Patents

Abstract

本发明公开了基于扭摆测量系统的高平稳电推力器的推力测评方法，涉及航天推进技术及测试技术领域，将电推力器搭载在扭摆测量系统上，根据待测推力作用下系统响应测量值，测量和评估推力，为高平稳电推力器的推力测评提供了有效方法。采用推力积分方程辛普森离散化反演测量推力的方法，获得推力测量值；用推力噪声作用下系统响应的方差分析方法，获得推力噪声方差预估计值。本发明采用推力蒙特卡洛模拟与系统响应对比寻优评估推力方法，搜索到优化的平均推力和推力噪声方差，实现了平均推力和推力噪声方差的高精度搜索优化求解，为高平稳推力器提供了有效推力测评方法。

Description

基于扭摆测量系统的高平稳电推力器的推力测评方法

技术领域

本发明涉及航天推进技术及测试技术领域，具体涉及基于扭摆测量系统的高平稳电推力器的推力测评方法。

背景技术

扭摆测量系统具有严格的二阶振动微分方程、环境噪声干扰小、承重能力强等特点，适合各种推力器的搭载测量推力，推力测试条件与推力器在空间工作条件基本一致，因此推力测评结果可信度高。

高平稳电推力器要求推力输出很平稳且起伏波动幅度很小，一般采用平均推力和零均值推力噪声来描述，对推力测量和评估方法提出了新要求。例如，空间引力波探测中对微牛级精密微推力器，在5～50μN推力范围内，要求推力噪声小于0.1μN/(Hz)^1/2，推力分辨率小于或等于0.1μN，即为高平稳电推力器。

现有的星载推力器推力测量和评估方法主要缺陷为：

(1)高平稳推力器用于高精度姿轨控情况，其推力采用平均推力与推力噪声之和描述，对平均推力和推力噪声方差的测评精度要求很高，已有推力测评方法亟待改进。

(2)以往根据稳态系统响应测量稳态推力作为平均推力估计值，测量得到的是足够长时间后的平均推力测量值，没有考虑进入稳态前平均推力的测量问题。

(3)不了解推力噪声的功率谱密度、自相关函数、方差，对其系统响应的功率谱密度、自相关函数、方差的影响规律，没有解决推力噪声方差的测量问题。

发明内容

有鉴于此，本发明提供了基于扭摆测量系统的高平稳电推力器的推力测评方法，将电推力器搭载在扭摆测量系统上，根据待测推力作用下系统响应测量值，测量和评估推力，为高平稳电推力器的推力测评提供了有效方法。

为达到上述目的，本发明的技术方案为：采用扭摆测量系统搭载高平稳电推力器进行推力测评，包括如下步骤：

步骤1、由扭摆振动微分方程构建推力积分方程，并将推力积分方程进行辛普森离散化为推力线性方程，对推力线性方程进行求解得到推力测量值。

步骤2、针对推力测量值进行线性拟合，寻找推力测量值的平均位置直线，给出平均推力测量值。

步骤3、进行推力作用下系统的响应分析，包括两部分：一部分是在平均推力作用下系统响应，另一部分是在推力噪声作用下系统响应。

分析获得推力噪声作用下系统响应的功率谱密度、自相关函数、方差特性；

建立推力噪声方差与系统响应方差之间定量关系，根据推力噪声作用下系统响应的方差特性，得到推力噪声方差预估值。

步骤4、以平均推力测量值与推力噪声方差预估值为基准，通过平均推力、推力噪声方差的增大或减小搜索方法，结合蒙特卡洛模拟抽样方法，获得一系列不同平均推力、不同推力噪声方差的模拟推力。

根据一系列的模拟推力，采用辛普森数值离散化方法，利用推力积分方程，计算一系列的模拟推力作用下系统响应。

将一系列的模拟推力作用下系统响应，与待测推力作用下系统响应测量值进行逐一对比，若模拟推力作用下系统响应落入待测推力作用下系统响应包络线内，以当前模拟推力对应的平均推力和推力噪声方差作为优化的结果。

进一步地，步骤1具体包括如下步骤：

步骤101、由扭摆振动微分方程构建推力积分方程，具体为：

扭摆振动微分方程为：

式中，f(t)为t时刻的施加在扭摆测量系统上的推力，t为时间变量；θ(t)为t时刻扭摆测量系统的系统响应；

为扭摆测量系统的系统响应的二阶导数；ζ为阻尼比，ω_n为扭摆测量系统的固有振动频率；J为转动惯量；L_f为力臂；t为时间变量；

由此构建推力积分方程，为

其中扭摆测量系统的固有振动频率为

k为扭摆测量系统的扭转刚度系数；ω_n(t-τ)为扭摆测量系统在t-τ时刻的固有振动频率；扭摆测量系统的阻尼比为

c为扭摆测量系统的阻尼系数；扭摆测量系统的振动频率为

ω_d(t-τ)为扭摆测量系统在t-τ时刻的振动频率；

步骤102、对推力积分方程进行辛普森离散化为推力线性方程组，具体包括如下步骤：

步骤1020、在无测量噪声干扰的理想情况下，推力f(t)作用下，扭摆的系统响应为θ(t)；由于环境噪声干扰及位移传感器噪声干扰，实际系统响应测量值为Θ(t)＝θ(t)+Δθ(t)，Δθ(t)为系统响应噪声干扰；

已知系统响应测量值为[t_i,Θ(t_i)]，采样时间步长为h，t_i为第i个采样时间取值为：i＝0,1,2,…,n；t_i＝ih，当t₀＝0时Θ(t₀)＝0，令Θ_i＝Θ(t_i)；

当t_i＝ih时，推力积分方程为

式中，C_f为指代值，用于指代Jω_d/L_f；ω_d(ih-τ)为扭摆测量系统在第ih-τ时刻的振动频率；ω_n(ih-τ)为扭摆测量系统在第ih-τ时刻的固有振动频率；被积分函数F(τ)表示τ时刻测量得到的推力，与待测推力f(t)有区别；

将时间区间[t₀＝0,t_i＝ih](i＝1,2,…)划分为i等分，第j个节点为τ_j＝jh，j取值为j＝0,1,2,…,i；令F_j＝F(τ_j)表示τ_j时刻测量得到的推力，推力积分方程辛普森离散化为推力线性方程组过程，具体为：

步骤1021：当i＝1时，采用梯形公式起步计算，为

可得10号系数a₁₀：a₁₀F₀＝C_fΘ₁，

步骤1022：当i＝2时，采用3点辛普森公式计算，为

可得20号系数

21号系数

a₂₀F₀+a₂₁F₁＝C_fΘ₂。

步骤1023：当i＝3时，采用4点辛普森公式计算，为

可得30号系数

31号系数

以及32号系数

a₃₀F₀+a₃₁F₁+a₃₂F₂＝C_fΘ₃

步骤1024：当i＝2k；k＝2,3,…时，采用3点辛普森公式计算，为

可得i0号系数a_i0以及il号系数a_il：

步骤1025：当i＝2k+1(k＝2,3,…)时，在i＝2k+1个点中，前面的i＝0,1,…,2k-3,2k-2个点采用3点辛普森公式计算，最后的i＝2k-2,2k-1,2k,2k+1等4个点采用4点辛普森公式计算，为

得到：

从而得到下三角推力线性方程组，为

表示为矩阵形式，为

求解下三角推力线性方程组，推力测量值F_i，为

式中，矩阵

为下三角矩阵[a_ij]的逆阵，

也是下三角矩阵，采用下三角矩阵快速求逆方法计算。

进一步地，步骤2中，针对推力测量值进行线性拟合，寻找推力测量值的平均位置直线，给出平均推力测量值，具体为：

已知推力测量值为(t_i,F_i)；i＝0,1,2,…,n-1，采用线性拟合方法求解平均推力，设推力测量值满足

F_i＝α+βt_i+ε_i,(i＝0,1,2,…,n-1) (24)

式中，α为第一待定系数；β为第二待定系数，ε_i为独立等方差的正态随机变量。

根据最小二乘方法令残差的平方和J最小

可得第一待定系数α和第二待定系数β的估计值

和

分别为

式中，符号

表示样本均值，

为时间的样本均值；

为推力测量的样本均值，为

从而，得到平均推力与时间线性关系，即平均推力的估计值，为

进一步地，步骤3中，分析获得推力噪声作用下系统响应的功率谱密度、自相关函数、方差特性；建立推力噪声方差与系统响应方差之间定量关系，根据推力噪声作用下系统响应的方差特性，得到推力噪声方差预估值，具体包括如下步骤：

推力噪声Δf(t)的系统响应满足

设推力噪声Δf(t)的拉普拉斯变换为F(s)＝L[Δf(t)]，系统响应θ_Δf(t)的拉普拉斯变换Θ(s)＝L[θ_Δf(t)]，由扭摆振动方程可得传递函数为

两边取拉普拉斯反变换，可得系统单位脉冲响应函数，为

并且有

其中iω为复变量；ζ为阻尼比；

令

可得

ω_1,2＝±ω_d-ζω_ni

ω_1,2为系统稳定时单位脉冲响应函数振动频率的第一解；

令

可得

ω_3,4＝±ω_d+ζω_ni

ω_3,4为系统稳定时单位脉冲响应函数振动频率的第二解；

由于Δf(t)为零均值白噪声，功率谱密度为S_Δf(ω)＝σ²，σ²为常数，则自相关函数为R_Δf(s)＝σ²δ(s)，此时系统响应也是零均值平稳过程θ_Δf(t)，系统响应的功率谱密度为

S_Δθ(ω)＝H(-iω)H(iω)S_Δf(ω) (33)

系统响应的自相关函数为

系统响应的自相关函数R_Δθ(τ)转换为

式中，振动固有频率

自相关函数R_θ(τ)是偶函数，上述公式是τ＞0情况下得到的；

当τ→0时，系统响应的自相关函数为

式中，

推力噪声Δf(t)的方差为σ²，其系统响应方差为

则有

在高平稳电推力器的推力f(t)＝f_a+Δf(t)作用下，系统响应测量值为[t_i,Θ(t_i)](i＝0,1,2,…,n)，系统响应进入稳态后，系统响应测量值为[t_i,Θ(t_i)]；i＝m,m+1,…,n；m＜n，此时，推力噪声作用下系统响应方差的估计值，为

从而，得到推力噪声方差的预估计值

进一步地，步骤4包括如下具体步骤：

步骤401、以平均推力测量值与推力噪声方差预估值为基准，通过平均推力、推力噪声方差的增大或减小搜索方法，结合蒙特卡洛模拟抽样方法，获得一系列不同平均推力、不同推力噪声方差的模拟推力；包括如下具体步骤：

步骤4011、采用推力测量值线性拟合测量平均推力方法，平均推力估计值为

对于高平稳推力器的平均推力估计值为：

步骤4012：推力噪声ΔF(τ)服从零均值、方差为σ²的正态分布，第j个抽样值ΔF_j为

式中，r_j1和r_j2为(0,1)区间均匀分布随机数；j＝0,1,2,…；

步骤4013：推力蒙特卡洛模拟抽样值为

步骤4014：平均推力、推力噪声方差的增大或减小搜索方法为：

S1，平均推力增大搜索时

Δα＞0为平均推力步进值，平均推力减小搜索时

S2，推力噪声方差增大搜索时有σ＝σ+Δσ，Δσ＞0为推力噪声方差步进值，推力噪声方差减小搜索时有σ＝σ-Δσ；

S3，重复S1和S2，确定一系列不同平均推力、不同推力噪声方差的模拟推力；

步骤402、模拟推力作用下系统响应辛普森计算；时间采样步长为h，采样时刻t_i＝ih；i＝0,1,2,…，根据推力积分方程，可得采样时刻t_i对应的系统响应θ(t_i)：

当i＝0时θ(t₀)＝θ(0)＝0，将[t₀＝0,t_i＝ih]划分为i等分，i为得到的子区间数目，令中间量g(t_i,τ_j)为：

式中，

对应的系统响应[t_i,θ(t_i)](i＝1,2,…)计算方法，为

步骤4021：当子区间数目i＝1时，采用梯形公式起步计算。

步骤4022：当子区间数目i＝2时，采用3点辛普森公式计算，为

步骤4023：当子区间数目i＝3时，采用4点辛普森公式计算，为

步骤4024：当子区间数目i＝2k(k＝2,3,4,…)时，采用3点辛普森公式计算；

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

步骤4025：当子区间数目i＝2k+1；k＝2,3,4,…时，交替采用3点和4点辛普森公式计算。

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

步骤403、系统响应对比寻优评估推力，包括如下步骤：

步骤4031：在高平稳电推力器的待测推力f(t)＝f_a+Δf(t)作用下，获得系统响应测量值为[t_i,Θ(t_i)]；i＝0,1,2,…,n。

步骤4032：采用推力蒙特卡洛模拟方法，获得不同平均推力、不同推力噪声方差的模拟推力抽样值

步骤4033：采用模拟推力作用下系统响应辛普森计算方法，获得一系列不同平均推力、不同推力噪声方差的模拟推力作用下系统响应[t_i,θ(t_i)](i＝1,2,…)。

步骤4034：将一系列模拟推力作用下系统响应[t_i,θ(t_i)]，与待测推力作用下系统响应测量值[t_i,Θ(t_i)]，逐一对比，当[t_i,θ(t_i)]曲线第一次覆盖[t_i,Θ(t_i)]曲线时，或第一次落入[t_i,Θ(t_i)]曲线的上下包络线内，便搜索到优化的平均推力和推力噪声方差。

有益效果：

(1)本发明采用推力积分方程辛普森离散化反演测量推力的方法，获得推力测量值；利用平均推力不受系统响应噪声干扰的特点，采用推力测量值线性拟合测量平均推力的方法，根据推力测量值计算平均推力测量值，解决了高平稳推力器的平均推力测量问题。

(2)本发明采用推力噪声作用下系统响应的方差分析方法，明确了高平稳推力器系统响应包括两部分即平均推力作用下系统响应、推力噪声作用下系统响应；建立了推力噪声方差与系统响应方差之间定量关系；由系统响应测量值计算系统响应方差，再由系统响应方差获得推力噪声方差预估计值，解决了高平稳推力器的推力噪声方差测量问题。

(3)本发明采用推力蒙特卡洛模拟与系统响应对比寻优评估推力方法，一是，以所得到的平均推力测量值、推力噪声方差预估计值为基准，进行搜索，获得一系列不同平均推力、不同推力噪声方差的模拟推力；二是，采用辛普森数值离散化方法计算一系列模拟推力作用下系统响应；三是，将一系列模拟推力作用下系统响应，与待测推力作用下系统响应测量值，逐一对比，当模拟推力作用下系统响应刚好落入待测推力作用下系统响应包络线内，便搜索到优化的平均推力和推力噪声方差，实现了平均推力和推力噪声方差的高精度搜索优化求解，为高平稳推力器提供了有效推力测评方法。

(4)本发明所提供的推力测评方法，既适合系统响应的稳态阶段，也适合系统响应的非稳态初始阶段。

附图说明

图1为本发明技术方案的实施流程示意图；

图2为待测推力作用下系统响应的变化示意图；

图3为待测推力作用下系统响应及上下包络线(局部放大)示意图；

图4为推力计算值示意图；

图5系统响应对比寻优评估推力示意图。

具体实施方式

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

附图1给出了本发明技术方案的实施流程。本发明的目的是基于扭摆测量系统，提供一种星载高平稳电推力器推力测评方法。

本发明通过以下技术方案实现。

步骤1.推力积分方程辛普森离散化反演测量推力

(1)由扭摆振动微分方程构建推力积分方程

由扭摆振动微分方程构建推力积分方程，建立系统响应与推力之间积分方程关系式。

(2)推力积分方程辛普森离散化为推力线性方程组并求解推力

一是采用辛普森数值积分方法，将推力积分方程离散化为推力线性方程组；二是利用推力线性方程组的系数矩阵为下三角矩阵特点，采用下三角矩阵快速求逆方法，求解推力线性方程组，实现由系统响应测量值反演计算推力的目的。

步骤2.推力测量值线性拟合测量平均推力

针对很小的系统响应噪声干扰将造成推力计算值很大波动的难点，利用平均推力不受系统响应噪声干扰的特点，通过推力积分方程辛普森离散化反演测量推力的方法，计算推力测量值，再采用一元线性拟合方法，寻找推力测量值的平均位置直线，给出平均推力测量值。

步骤3.推力噪声作用下系统响应的方差分析

(1)平均推力作用下系统响应特性分析

通过扭摆测量系统的振动微分方程分析，表明平均推力作用下系统响应为开始起伏波动、后期进入水平稳态的确定曲线。

(2)推力噪声作用下系统响应功率谱密度、自相关函数、方差特性分析

通过扭摆测量系统的振动微分方程分析，表明在零均值平稳随机过程的推力噪声作用下，扭摆系统响应也是零均值的平稳随机过程，并且建立了两个平稳随机过程的功率谱密度、自相关函数、方差的定量关系。

(3)推力噪声方差的预估计

利用推力作用下系统响应方差在稳态和非稳态过程都不变的特性，根据所建立的推力噪声方差与系统响应方差的定量关系，由系统响应测量值计算系统响应方差，再由系统响应方差给出推力噪声方差的预估计值，为推力噪声方差的优化评估提供了基础。

步骤4.推力蒙特卡洛模拟与系统响应对比寻优评估推力

(1)推力蒙特卡洛模拟方法

以所得到的平均推力测量值、推力噪声方差的预估计值为基准，通过平均推力、推力噪声方差的增大或减小搜索方法，结合蒙特卡洛模拟抽样方法，获得一系列不同平均推力、不同推力噪声方差的模拟推力。

(2)模拟推力作用下系统响应辛普森计算

根据所获得的不同平均推力、不同推力噪声方差的模拟推力，采用辛普森数值离散化方法，利用推力积分方程，计算一系列模拟推力作用下系统响应。

(3)系统响应对比寻优评估推力

将不同平均推力、不同推力噪声方差的模拟推力作用下系统响应，与待测推力作用下系统响应测量值，逐一对比，当模拟推力作用下系统响应刚好落入待测推力作用下系统响应包络线内，便搜索到优化的平均推力和推力噪声方差。

实例：采用扭摆测量系统搭载某型高平稳电推力器，测评其推力，根据以往经验该推力器推力小于50μN，要求推力分辨率≤0.1μN、推力噪声<0.1μN/(Hz)^1/2，附图2为待测推力作用下系统响应的变化。

所采用扭摆测量系统的系统参数估计值，扭转刚度系数、振动频率、阻尼比分别为

k＝0.232N·m/rad，ω_d＝0.43rad/s，ζ＝0.21

采用变极距式电容传感器测量位移，位移分辨率为30nm，力臂为L_f＝0.5m，测量臂为L_s＝0.5m。

高平稳电推力器的推力表示为f(t)＝f_a+Δf(t)，f_a为平均推力，Δf(t)为推力噪声(推力的起伏波动部分)，推力噪声采用零均值平稳随机过程描述。推力噪声可看作为零均值白噪声，功率谱密度为S_Δf(ω)＝(0.1)²(单位：(μN)²/(Hz))，其自相关函数为R_Δf(τ)＝(0.1)²δ(τ)(单位：(μN)²)，也就是推力噪声在任意两个不同时刻独立同分布。

与待测推力f(t)＝f_a+Δf(t)相对应，测量得到的推力为F(t)＝F_a+ΔF(t)，F_a为平均推力f_a的测量评估值，ΔF(t)为推力噪声Δf(t)的测量评估值。

1.推力积分方程辛普森离散化反演测量推力

推力积分方程辛普森离散化反演测量推力的方法：首先，由扭摆振动微分方程构建推力积分方程，建立系统响应与推力之间积分方程关系式；其次，采用辛普森数值积分方法，将推力积分方程离散化为推力线性方程组，通过下三角推力线性方程组求解，实现由系统响应测量值反演计算推力的目的。

(1)由扭摆振动微分方程构建推力积分方程

扭摆振动微分方程为

式中，f(t)为t时刻的施加在扭摆测量系统上的推力，t为时间变量；θ(t)为t时刻扭摆测量系统的系统响应；

为扭摆测量系统的系统响应的二阶导数；ζ为阻尼比，ω_n为扭摆测量系统的固有振动频率；J为转动惯量；L_f为力臂；t为时间变量。

由此构建推力积分方程，为

其中扭摆测量系统的固有振动频率为

k为扭摆测量系统的扭转刚度系数；ω_n(t-τ)为扭摆测量系统在t-τ时刻的固有振动频率；扭摆测量系统的阻尼比为

c为扭摆测量系统的阻尼系数；扭摆测量系统的振动频率为

ω_d(t-τ)为扭摆测量系统在t-τ时刻的振动频率。

从而建立了推力与系统响应的积分方程关系式。

扭摆测量系统的固有振动频率为

式中，k为扭摆测量系统的扭转刚度系数。扭摆测量系统的阻尼比为

式中，c为扭摆测量系统的阻尼系数。扭摆测量系统的振动频率为

(2)推力积分方程辛普森离散化为推力线性方程组并求解推力

以梯形公式起步，交替使用3点辛普森和4点辛普森公式，将推力积分方程辛普森离散化为推力线性方程组，通过下三角推力线性方程组求解计算推力。

在无测量噪声干扰的理想情况下，推力f(t)作用下，扭摆的系统响应为θ(t)，由于环境噪声干扰(环境位移激励、环境外力激励、推力加载冲击等干扰)及位移传感器噪声干扰，实际系统响应测量值为Θ(t)＝θ(t)+Δθ(t)，Δθ(t)为系统响应噪声干扰。

已知系统响应测量值为[t_i,Θ(t_i)](i＝0,1,2,…,n)，采样时间步长为h，t_i＝ih，当t₀＝0时Θ(t₀)＝0，令Θ_i＝Θ(t_i)。当t_i＝ih时，推力积分方程为

式中，C_f＝Jω_d/L_f，被积分函数F(τ)表示测量得到的推力，与待测推力f(t)有区别。

将时间区间[t₀＝0,t_i＝ih](i＝1,2,…划分为i等分，节点为τ_j＝jh(j＝0,1,2,…,i)，令F_j＝F(τ_j)，推力积分方程辛普森离散化为推力线性方程组过程，具体为：

步骤1：当i＝1时，采用梯形公式起步计算，为

可得

步骤2：当i＝2时，采用3点辛普森公式计算，为

可得

a₂₀F₀+a₂₁F₁＝C_fΘ₂ (8)

步骤3：当i＝3时，采用4点辛普森公式计算，为

可得

a₃₀F₀+a₃₁F₁+a₃₂F₂＝C_fΘ₃

步骤4：当i＝2k(k＝2,3,…)时，采用3点辛普森公式计算，为

可得

步骤5：当i＝2k+1(k＝2,3,…)时，在i＝2k+1个点中，前面的i＝0,1,…,2k-3,2k-2个点采用3点辛普森公式计算，最后的i＝2k-2,2k-1,2k,2k+1等4个点采用4点辛普森公式计算，为

可得

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

从而得到下三角推力线性方程组，为

表示为矩阵形式，为

求解下三角推力线性方程组，推力测量值，为

式中，矩阵

为下三角矩阵[a_ij]的逆阵，也是下三角矩阵，可采用下三角矩阵快速求逆方法计算，其中

附图2所示为待测推力作用下系统响应的变化，系统响应最大值为81.5μrad，稳态系统响应为53.9μrad。附图3为待测推力作用下系统响应及上下包络线(局部放大)，是附图2的局部放大图，约60s后系统响应进入稳态，并且在系统响应噪声干扰下，稳态系统响应的变化范围为53.25～54.75μrad。

根据附图2所示的待测推力作用下系统响应，采用推力积分方程辛普森离散化反演测量推力的方法，计算推力，附图4为推力计算值。从附图4可看出，由于系统响应噪声干扰，推力计算值波动范围很大为-400～400μN。说明很小的系统响应噪声干扰，将造成很大的推力计算值的波动。

2.推力测量值线性拟合测量平均推力

推力测量值线性拟合测量平均推力的方法：针对很小的系统响应噪声干扰将造成推力计算值很大波动的难点，利用平均推力不受系统响应噪声干扰特点，线性拟合寻找推力测量值的平均位置直线，给出平均推力测量值。

由于实际系统响应为Θ_k＝θ_k+Δθ_k，包含系统响应噪声干扰Δθ_k，将造成推力测量值波动，为

当时间步长h→0时

所以很小的系统响应噪声干扰Δθ_k，将造成很大的推力计算值的波动。

利用零均值噪声特点E(Δθ_k)＝0，方程两边取均值，可得平均推力为

式中，E(·)为求均值运算。显然，平均推力不受系统响应噪声干扰的影响。

已知推力计算值为(t_i,F_i)(i＝0,1,2,…,n-1)，采用线性拟合方法求解平均推力，设推力测量值满足

F_i＝α+βt_i+ε_i,(i＝0,1,2,…,n-1) (24)

式中，α和β为待定系数，ε_i为独立等方差的正态随机变量

根据最小二乘方法令残差的平方和最小

可得系数α和β的估计值

和

分别为

式中，符号

表示样本均值，为

从而，得到平均推力与时间线性关系，即平均推力的估计值，为

根据附图4所示推力计算值，采用推力测量值线性拟合测量平均推力的方法，可得平均推力估计值为

(单位：μN)

显然，平均推力随着时间变化很小，可取平均推力估计值为

F_a≈25.10732μN

3.推力噪声作用下系统响应的方差分析

推力噪声作用下系统响应的方差分析方法：首先，通过推力作用下系统响应分析，明确了推力作用下系统响应由两部分组成，一部分是在平均推力作用下系统响应(确定性部分)，一部分是在推力噪声作用下系统响应(平稳随机性部分)；其次，通过推力噪声作用下系统响应功率谱密度、自相关函数、方差分析，建立了推力噪声方差与系统响应方差之间定量关系；最后，根据系统响应测量值计算系统响应方差，由系统响应方差给出了推力噪声方差的预估计值，为推力噪声方差优化估计提供了基准值。

(1)平均推力作用下系统响应特性分析

平均推力作用下系统响应为开始起伏波动、后期进入水平稳态的确定曲线。

在高平稳电推力器的推力f(t)＝f_a+Δf(t)作用下，扭摆测量系统的系统响应为

即系统响应包括两部分，一是平均推力作用下系统响应；二是推力噪声作用下系统响应。

根据扭摆振动微分方程

可得平均推力作用下系统响应，为

当时间足够长时，稳态系统响应为

(2)推力噪声作用下系统响应功率谱密度、自相关函数、方差特性分析

首先，在零均值平稳随机过程的推力噪声作用下，扭摆系统响应也是零均值的平稳随机过程；其次，建立了推力噪声的功率谱密度、自相关函数、方差等，与系统响应的功率谱密度、自相关函数、方差等的定量关系。

根据扭摆振动微分方程

可得推力噪声作用下系统响应功率谱密度、自相关函数、方差特性。

推力噪声Δf(t)的系统响应满足

设推力噪声Δf(t)的拉普拉斯变换为F(s)＝L[Δf(t)]，系统响应θ_Δf(t)的拉普拉斯变换Θ(s)＝L[θ_Δf(t)]，由扭摆振动方程可得传递函数为

两边取拉普拉斯反变换，可得系统单位脉冲响应函数，为

并且有

其中iω为复变量；ζ为阻尼比；

令

可得

ω_1,2＝±ω_d-ζω_ni

ω_1,2为系统稳定时单位脉冲响应函数振动频率的第一解；

令

可得

ω_3,4＝±ω_d+ζω_ni

ω_3,4为系统稳定时单位脉冲响应函数振动频率的第二解；

由于Δf(t)为零均值白噪声，功率谱密度为S_Δf(ω)＝σ²(常数)，则自相关函数为R_Δf(s)＝σ²δ(s)，此时系统响应也是零均值平稳过程θ_Δf(t)，系统响应的功率谱密度为

S_Δθ(ω)＝H(-iω)H(iω)S_Δf(ω) (33)

系统响应的自相关函数为

由于被积分函数是有理函数，分母阶数高于分子，在实轴上没有极点，采用留数定理表示为

式中，τ＞0情况取上半复平面的极点，即ω_3,4＝±ω_d+ζω_ni，可得

其中

由于ω_3,4＝±ω_d+ζω_ni时，分母为零但分母导数不为零，可简化为

并且有

为了简化方便，令

z＝(ζω_n)²+ζω_nω_di，

可得

代入可得

因此，进一步简化为

可得，系统响应的自相关函数为

式中，振动固有频率

自相关函数R_θ(τ)是偶函数，上述公式是τ＞0情况下得到的。

当τ→0时，系统响应的自相关函数为

式中，

推力噪声Δf(t)的方差为σ²，其系统响应方差为

则有

(3)推力噪声方差的预估计

利用推力作用下系统响应方差在稳态和非稳态过程都不变的特性，根据所建立的推力噪声方差与系统响应方差的定量关系，由系统响应方差给出了推力噪声方差的预估计值，为推力噪声方差的优化评估提供了基础。

在高平稳电推力器的推力f(t)＝f_a+Δf(t)作用下，系统响应测量值为[t_i,Θ(t_i)](i＝0,1,2,…,n)，系统响应进入稳态后，系统响应测量值为[t_i,Θ(t_i)](i＝m,m+1,…,n)(m＜n)，此时，推力噪声作用下系统响应方差的估计值，为

从而，得到推力噪声方差的预估计值

根据附图2和附图3可知，60s后系统响应进入稳态，利用稳态系统响应测量值，计算稳态系统响应的均值和标准差为

并进一步，计算推力噪声标准差的预估计值

σ＝1.917078×10^-1μN

4.推力蒙特卡洛模拟与系统响应对比寻优评估推力

推力蒙特卡洛模拟与系统响应对比寻优评估推力的方法：首先，采用推力蒙特卡洛模拟方法，以平均推力估计值、推力噪声方差预估计值为基准，通过增大或减小平均推力、推力噪声方差进行上下搜索，获得一系列不同平均推力、不同方差的推力平稳随机过程；其次，采用模拟推力作用下系统响应辛普森计算方法，获得一系列不同平均推力、不同推力噪声方差作用下系统响应计算值；最后，采用系统响应对比寻优评估推力的方法，将所获得的一系列不同平均推力、不同推力噪声方差的模拟推力作用下系统响应计算值，与待测推力作用下系统响应测量值，逐一对比，当模拟推力作用下系统响应刚好落入待测推力作用下系统响应测量值时(或刚好落入其上下包络线内)，便得到待测推力的优化估计值。

(1)推力蒙特卡洛模拟方法

以平均推力估计值、推力噪声方差的预估计值为基准，通过平均推力、推力噪声方差的增大或减小搜索方法，结合蒙特卡洛模拟抽样方法，获得一系列不同平均推力、不同推力噪声方差的推力平稳随机过程。

推力蒙特卡洛模拟方法为：

步骤1：采用推力测量值线性拟合测量平均推力方法，平均推力估计值为

对于高平稳推力器的平均推力，可令

步骤2：设推力噪声ΔF(τ)服从零均值、方差为σ²的正态分布，抽样值为

式中，r_j1和r_j2(j＝0,1,2,…)为(0,1)区间均匀分布随机数。

步骤3：推力蒙特卡洛模拟抽样值为

步骤4：平均推力、推力噪声方差搜索方法为，一是，平均推力增大搜索时

平均推力减小搜索时

二是，推力噪声方差增大搜索时有σ＝σ+Δσ(Δσ＞0)，推力噪声方差减小搜索时有σ＝σ-Δσ(Δσ＞0)；三是，重复步骤1和步骤2，可得到一系列不同平均推力、不同推力噪声方差的模拟推力抽样值。

(2)模拟推力作用下系统响应辛普森计算

设时间采样步长为h，t_i＝ih(i＝0,1,2,…)，根据推力积分方程，可得

当i＝0时θ(t₀)＝θ(0)＝0，将[t₀＝0,t_i＝ih]划分为i等分，令

式中，

对应的系统响应[t_i,θ(t_i)](i＝1,2,…)计算方法，为

步骤1：当子区间数目i＝1时，采用梯形公式起步计算。

步骤2：当子区间数目i＝2时，采用3点辛普森公式计算，为

步骤3：当子区间数目i＝3时，采用4点辛普森公式计算，为

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

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

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

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

(3)系统响应对比寻优评估推力

步骤1：在高平稳电推力器的待测推力f(t)＝f_a+Δf(t)作用下，获得系统响应测量值为[t_i,Θ(t_i)](i＝0,1,2,…,n)。

步骤2：采用推力蒙特卡洛模拟方法，获得不同平均推力、不同推力噪声方差的模拟推力抽样值

步骤3：采用模拟推力作用下系统响应辛普森计算方法，获得一系列不同平均推力、不同推力噪声方差的模拟推力作用下系统响应[t_i,θ(t_i)](i＝1,2,…)。

步骤4：将一系列模拟推力作用下系统响应[t_i,θ(t_i)]，与待测推力作用下系统响应测量值[t_i,Θ(t_i)]，逐一对比，当[t_i,θ(t_i)]曲线刚好覆盖[t_i,Θ(t_i)]曲线时(或刚好落入其上下包络线内)，便搜索到优化的平均推力和推力噪声方差。

采用推力测量值线性拟合测量平均推力的方法，可取平均推力估计值为

F_a≈25.10732μN

采用推力噪声作用下系统响应的方差分析方法，计算推力噪声标准差的预估计值为

σ＝1.917078×10^-1μN

平均推力搜索步长为Δα＝0.01，推力噪声标准差搜索步长为Δσ＝0.01，当平均推力为F_a＝25.05μN和推力噪声标准差为σ＝0.09μN时，该模拟推力作用下系统响应刚好落入待测推力作用下系统响应测量值的上下包络线内，如附图5所示为系统响应对比寻优评估推力。因此，该推力器的推力表示为

F(t)＝25.05+ΔF(单位：μN)

式中，ΔF～N[0,(0.09)²]为正态分布随机变量，推力噪声标准差为0.09μN，满足推力噪声小于0.1μN要求。

扭摆测量系统，采用变极距式电容传感器测量位移，位移分辨率为30nm，力臂为L_f＝0.5m，测量臂为L_s＝0.5m，扭转刚度系数为k＝0.232N·m/rad，稳态扭转角为

式中，Δh为位移测量值。令位移分辨率为

可得，推力分辨率为

因此，推力分辨率为0.02784μN，满足小于或等于0.1μN要求。

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

Claims (5)

1.基于扭摆测量系统的高平稳电推力器的推力测评方法，其特征在于，采用扭摆测量系统搭载所述高平稳电推力器进行推力测评，包括如下步骤：

步骤1、由扭摆振动微分方程构建推力积分方程，并将所述推力积分方程进行辛普森离散化为推力线性方程，对所述推力线性方程进行求解得到推力测量值；

步骤2、针对所述推力测量值进行线性拟合，寻找推力测量值的平均位置直线，给出平均推力测量值；

步骤3、进行推力作用下系统的响应分析，包括两部分：一部分是在平均推力作用下系统响应，另一部分是在推力噪声作用下系统响应；

分析获得推力噪声作用下系统响应的功率谱密度、自相关函数、方差特性；

建立推力噪声方差与系统响应方差之间定量关系，根据推力噪声作用下系统响应的方差特性，得到推力噪声方差预估值；

步骤4、以所述平均推力测量值与所述推力噪声方差预估值为基准，通过平均推力、推力噪声方差的增大或减小搜索方法，结合蒙特卡洛模拟抽样方法，获得一系列不同平均推力、不同推力噪声方差的模拟推力；

根据一系列的模拟推力，采用辛普森数值离散化方法，利用推力积分方程，计算一系列的模拟推力作用下系统响应；

将一系列的模拟推力作用下系统响应，与待测推力作用下系统响应测量值进行逐一对比，若模拟推力作用下系统响应落入待测推力作用下系统响应包络线内，以当前模拟推力对应的平均推力和推力噪声方差作为优化的结果。

2.如权利要求1所述的方法，其特征在于，所述步骤1具体包括如下步骤：

步骤101、由扭摆振动微分方程构建推力积分方程，具体为：

所述扭摆振动微分方程为：

式中，f(t)为t时刻的施加在扭摆测量系统上的推力，t为时间变量；θ(t)为t时刻扭摆测量系统的系统响应；

为扭摆测量系统的系统响应的二阶导数；ζ为阻尼比，ω_n为扭摆测量系统的固有振动频率；J为转动惯量；L_f为力臂；t为时间变量；

由此构建推力积分方程，为

其中扭摆测量系统的固有振动频率为

k为扭摆测量系统的扭转刚度系数；ω_n(t-τ)为扭摆测量系统在t-τ时刻的固有振动频率；扭摆测量系统的阻尼比为

c为扭摆测量系统的阻尼系数；扭摆测量系统的振动频率为

ω_d(t-τ)为扭摆测量系统在t-τ时刻的振动频率；

步骤102、对所述推力积分方程进行辛普森离散化为推力线性方程组，具体包括如下步骤：

步骤1020、在无测量噪声干扰的理想情况下，推力f(t)作用下，扭摆的系统响应为θ(t)；由于环境噪声干扰及位移传感器噪声干扰，实际系统响应测量值为Θ(t)＝θ(t)+Δθ(t)，Δθ(t)为系统响应噪声干扰；

已知系统响应测量值为[t_i,Θ(t_i)]，采样时间步长为h，t_i为第i个采样时间取值为：i＝0,1,2,…,n；t_i＝ih，当t₀＝0时Θ(t₀)＝0，令Θ_i＝Θ(t_i)；

当t_i＝ih时，推力积分方程为

式中，C_f为指代值，用于指代Jω_d/L_f；ω_d(ih-τ)为扭摆测量系统在第ih-τ时刻的振动频率；ω_n(ih-τ)为扭摆测量系统在第ih-τ时刻的固有振动频率；被积分函数F(τ)表示τ时刻测量得到的推力，与待测推力f(t)有区别；

将时间区间[t₀＝0,t_i＝ih](i＝1,2,…)划分为i等分，第j个节点为τ_j＝jh，j取值为j＝0,1,2,…,i；令F_j＝F(τ_j)表示τ_j时刻测量得到的推力，推力积分方程辛普森离散化为推力线性方程组过程，具体为：

步骤1021：当i＝1时，采用梯形公式起步计算，为

可得10号系数a₁₀：a₁₀F₀＝C_fΘ₁，

步骤1022：当i＝2时，采用3点辛普森公式计算，为

可得20号系数

21号系数

a₂₀F₀+a₂₁F₁＝C_fΘ₂；

步骤1023：当i＝3时，采用4点辛普森公式计算，为

可得30号系数

31号系数

以及32号系数

a₃₀F₀+a₃₁F₁+a₃₂F₂＝C_fΘ₃

步骤1024：当i＝2k；k＝2,3,…时，采用3点辛普森公式计算，为

可得i0号系数a_i0以及il号系数a_il：

步骤1025：当i＝2k+1(k＝2,3,··)时，在i＝2k+1个点中，前面的i＝0,1,…,2k-3,2k-2个点采用3点辛普森公式计算，最后的i＝2k-2,2k-1,2k,2k+1等4个点采用4点辛普森公式计算，为

得到：

从而得到下三角推力线性方程组，为

表示为矩阵形式，为

求解下三角推力线性方程组，推力测量值F_i，为

式中，矩阵

为下三角矩阵[a_ij]的逆阵，

也是下三角矩阵，采用下三角矩阵快速求逆方法计算。

3.如权利要求1或2所述的方法，其特征在于，所述步骤2中，针对所述推力测量值进行线性拟合，寻找推力测量值的平均位置直线，给出平均推力测量值，具体为：

已知推力测量值为(t_i,F_i)；i＝0,1,2,…,n-1，采用线性拟合方法求解平均推力，设推力测量值满足

F_i＝α+βt_i+ε_i,(i＝0,1,2,…,n-1) (24)

式中，α为第一待定系数；β为第二待定系数，ε_i为独立等方差的正态随机变量；

根据最小二乘方法令残差的平方和J最小

可得第一待定系数α和第二待定系数β的估计值

和

分别为

式中，符号“—”表示样本均值，

为时间的样本均值；

为推力测量的样本均值，为

从而，得到平均推力与时间线性关系，即平均推力的估计值，为

4.如权利要求3所述的方法，其特征在于，所述步骤3中，分析获得推力噪声作用下系统响应的功率谱密度、自相关函数、方差特性；建立推力噪声方差与系统响应方差之间定量关系，根据推力噪声作用下系统响应的方差特性，得到推力噪声方差预估值，具体包括如下步骤：

推力噪声Δf(t)的系统响应满足

设推力噪声Δf(t)的拉普拉斯变换为F(s)＝L[Δf(t)]，系统响应θ_Δf(t)的拉普拉斯变换Θ(s)＝L[θ_Δf(t)]，由扭摆振动方程可得传递函数为

两边取拉普拉斯反变换，可得系统单位脉冲响应函数，为

并且有

其中iω为复变量；ζ为阻尼比；

令

可得

ω_1,2＝±ω_d-ζω_ni

ω_1,2为系统稳定时单位脉冲响应函数振动频率的第一解；

令

可得

ω_3,4＝±ω_d+ζω_ni

ω_3,4为系统稳定时单位脉冲响应函数振动频率的第二解；

由于Δf(t)为零均值白噪声，功率谱密度为S_Δf(ω)＝σ²，σ²为常数，则自相关函数为R_Δf(s)＝σ²δ(s)，此时系统响应也是零均值平稳过程θ_Δf(t)，系统响应的功率谱密度为

S_Δθ(ω)＝H(-iω)H(iω)S_Δf(ω) (33)

系统响应的自相关函数为

系统响应的自相关函数R_Δθ(τ)转换为

式中，振动固有频率

自相关函数R_θ(τ)是偶函数，上述公式是τ＞0情况下得到的；

当τ→0时，系统响应的自相关函数为

式中，

推力噪声Δf(t)的方差为σ²，其系统响应方差为

则有

在高平稳电推力器的推力f(t)＝f_a+Δf(t)作用下，系统响应测量值为[t_i,Θ(t_i)](i＝0,1,2,…,n)，系统响应进入稳态后，系统响应测量值为[t_i,Θ(t_i)]；i＝m,m+1,0，…,n；m＜n，此时，推力噪声作用下系统响应方差的估计值，为

从而，得到推力噪声方差的预估计值

5.如权利要求1、2、3、4任一所述的方法，其特征在于，所述步骤4包括如下具体步骤：

步骤401、以所述平均推力测量值与所述推力噪声方差预估值为基准，通过平均推力、推力噪声方差的增大或减小搜索方法，结合蒙特卡洛模拟抽样方法，获得一系列不同平均推力、不同推力噪声方差的模拟推力；包括如下具体步骤：

步骤4011、采用所述的推力测量值线性拟合测量平均推力方法，平均推力估计值为

对于高平稳推力器的平均推力估计值为：

步骤4012：推力噪声ΔF(τ)服从零均值、方差为σ²的正态分布，第j个抽样值ΔF_j为

式中，r_j1和r_j2为(0,1)区间均匀分布随机数；j＝0,1,2,…；

步骤4013：推力蒙特卡洛模拟抽样值为

步骤4014：平均推力、推力噪声方差的增大或减小搜索方法为：

S1，平均推力增大搜索时

Δα＞0为平均推力步进值，平均推力减小搜索时

S2，推力噪声方差增大搜索时有σ＝σ+Δσ，Δσ＞0为推力噪声方差步进值，推力噪声方差减小搜索时有σ＝σ-Δσ；

S3，重复S1和S2，确定一系列不同平均推力、不同推力噪声方差的模拟推力；

步骤402、模拟推力作用下系统响应辛普森计算；时间采样步长为h，采样时刻t_i＝ih；i＝0,1,2,…，根据推力积分方程，可得采样时刻t_i对应的系统响应θ(t_i)：

当i＝0时θ(t₀)＝θ(0)＝0，将[t₀＝0,t_i＝ih]划分为i等分，i为得到的子区间数目，令中间量g(t_i,τ_j)为：

式中，

对应的系统响应[t_i,θ(t_i)](i＝1,2,…)计算方法，为

步骤4021：当子区间数目i＝1时，采用梯形公式起步计算：

步骤4022：当子区间数目i＝2时，采用3点辛普森公式计算，为

步骤4023：当子区间数目i＝3时，采用4点辛普森公式计算，为

步骤4024：当子区间数目i＝2k(k＝2,3,4,…)时，采用3点辛普森公式计算；

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

步骤4025：当子区间数目i＝2k+1；k＝2,3,4,…时，交替采用3点和4点辛普森公式计算；

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

步骤403、系统响应对比寻优评估推力，包括如下步骤：

步骤4031：在高平稳电推力器的待测推力f(t)＝f_a+Δf(t)作用下，获得系统响应测量值为[t_i,Θ(t_i)]；i＝0,1,2,…,n；

步骤4032：采用所述的推力蒙特卡洛模拟方法，获得不同平均推力、不同推力噪声方差的模拟推力抽样值

步骤4033：采用所述的模拟推力作用下系统响应辛普森计算方法，获得一系列不同平均推力、不同推力噪声方差的模拟推力作用下系统响应[t_i,θ(t_i)](i＝1,2,…)；

步骤4034：将一系列模拟推力作用下系统响应[t_i,θ(t_i)]，与待测推力作用下系统响应测量值[t_i,Θ(t_i)]，逐一对比，当[t_i,θ(t_i)]曲线第一次覆盖[t_i,Θ(t_i)]曲线时，或第一次落入[t_i,Θ(t_i)]曲线的上下包络线内，便搜索到优化的平均推力和推力噪声方差。

Images