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

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

【技术实现步骤摘要】
基于扭摆测量系统的高平稳电推力器的推力测评方法
本专利技术涉及航天推进技术及测试
，具体涉及基于扭摆测量系统的高平稳电推力器的推力测评方法。
技术介绍
扭摆测量系统具有严格的二阶振动微分方程、环境噪声干扰小、承重能力强等特点，适合各种推力器的搭载测量推力，推力测试条件与推力器在空间工作条件基本一致，因此推力测评结果可信度高。高平稳电推力器要求推力输出很平稳且起伏波动幅度很小，一般采用平均推力和零均值推力噪声来描述，对推力测量和评估方法提出了新要求。例如，空间引力波探测中对微牛级精密微推力器，在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为转动惯量；Lf为力臂；t为时间变量；由此构建推力积分方程，为其中扭摆测量系统的固有振动频率为k为扭摆测量系统的扭转刚度系数；ωn(t-τ)为扭摆测量系统在t-τ时刻的固有振动频率；扭摆测量系统的阻尼比为c为扭摆测量系统的阻尼系数；扭摆测量系统的振动频率为ωd(t-τ)为扭摆测量系统在t-τ时刻的振动频率；步骤102、对推力积分方程进行辛普森离散化为推力线性方程组，具体包括如下步骤：步骤1020、在无测量噪声干扰的理想情况下，推力f(t)作用下，扭摆的系统响应为θ(t)；由于环境噪声干扰及位移传感器噪声干扰，实际系统响应测量值为Θ(t)＝θ(t)+Δθ(t)，Δθ(t)为系统响应噪声干扰；已知系统响应测量值为[ti,Θ(ti)]，采样时间步长为h，ti为第i个采样时间取值为：i＝0,1,2,…,n；ti＝ih，当t0＝0时Θ(t0)＝0，令Θi＝Θ(ti)；当ti＝ih时，推力积分方程为式中，Cf为指代值，用于指代Jωd/Lf；ωd(ih-τ)为扭摆测量系统在第ih-τ时刻的振动频率；ωn(ih-τ)为扭摆测量系统在第ih-τ时刻的固有振动频率；被积分函数F(τ)表示τ时刻测量得到的推力，与待测推力f(t)有区别；将时间区间[t0＝0,ti＝ih](i＝1,2,…)划分为i等分，第j个节点为τj＝jh，j取值为j＝0,1,2,…,i；令Fj＝F(τj)表示τj时刻测量得到的推力，推力积分方程辛普森离散化为推力线性方程组过程，具体为：步骤1021：当i＝1时，采用梯形公式起步计算，为可得10号系数a10：a10F0＝CfΘ1，步骤1022：当i＝2时，采用3点辛普森公式计算，为可得20号系数21号系数a20F0+a21F1＝CfΘ2。步骤1023：当i＝3时，采用4点辛普森公式计算，为可得30号系数31号系数以及32号系数a30F0+a31F1+a32F2＝CfΘ3步骤1024：当i＝2k；k＝2,3,…时，采用3点辛普森公式计算，为可得i0号系数ai0以及il号系数ail：步骤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点辛普森公式计算，为得到：从而得到下三角推力线性方程组，为表示为矩阵形式，为求解下三角推力线性方程组，推力测量值Fi，为式中，矩阵为下三角矩阵[aij]的逆阵，也是下三角矩阵，采用下三角矩阵快速求逆方法计算。进一步地，步骤2中，针对推力测量值进行线性拟合，寻找推力测量值的平均位置直线，给出平均推力测量值，具体为：已知推力测量值为(ti,Fi)；i＝0,1,2,…,n-1，采用线性拟合方法求解平均推力，设推力测量值满足Fi＝α+βti+εi,(i＝0,1,2,…,n-1)(24)式中，α为第一待定系数；β为第二待定系数，εi为独立等方差的正态随机变量。根据最小二乘方法令残差的平方和J最小可得第一待定系数α和第二待定系数β的估计值和分别为式中，符号表示样本均值，为时间的样本均值；为推力测量的样本均值，为从而，得到平均推力本文档来自技高网...

【技术保护点】
1.基于扭摆测量系统的高平稳电推力器的推力测评方法，其特征在于，采用扭摆测量系统搭载所述高平稳电推力器进行推力测评，包括如下步骤：/n步骤1、由扭摆振动微分方程构建推力积分方程，并将所述推力积分方程进行辛普森离散化为推力线性方程，对所述推力线性方程进行求解得到推力测量值；/n步骤2、针对所述推力测量值进行线性拟合，寻找推力测量值的平均位置直线，给出平均推力测量值；/n步骤3、进行推力作用下系统的响应分析，包括两部分：一部分是在平均推力作用下系统响应，另一部分是在推力噪声作用下系统响应；/n分析获得推力噪声作用下系统响应的功率谱密度、自相关函数、方差特性；/n建立推力噪声方差与系统响应方差之间定量关系，根据推力噪声作用下系统响应的方差特性，得到推力噪声方差预估值；/n步骤4、以所述平均推力测量值与所述推力噪声方差预估值为基准，通过平均推力、推力噪声方差的增大或减小搜索方法，结合蒙特卡洛模拟抽样方法，获得一系列不同平均推力、不同推力噪声方差的模拟推力；/n根据一系列的模拟推力，采用辛普森数值离散化方法，利用推力积分方程，计算一系列的模拟推力作用下系统响应；/n将一系列的模拟推力作用下系统响应，与待测推力作用下系统响应测量值进行逐一对比，若模拟推力作用下系统响应落入待测推力作用下系统响应包络线内，以当前模拟推力对应的平均推力和推力噪声方差作为优化的结果。/n...

【技术特征摘要】
1.基于扭摆测量系统的高平稳电推力器的推力测评方法，其特征在于，采用扭摆测量系统搭载所述高平稳电推力器进行推力测评，包括如下步骤：
步骤1、由扭摆振动微分方程构建推力积分方程，并将所述推力积分方程进行辛普森离散化为推力线性方程，对所述推力线性方程进行求解得到推力测量值；
步骤2、针对所述推力测量值进行线性拟合，寻找推力测量值的平均位置直线，给出平均推力测量值；
步骤3、进行推力作用下系统的响应分析，包括两部分：一部分是在平均推力作用下系统响应，另一部分是在推力噪声作用下系统响应；
分析获得推力噪声作用下系统响应的功率谱密度、自相关函数、方差特性；
建立推力噪声方差与系统响应方差之间定量关系，根据推力噪声作用下系统响应的方差特性，得到推力噪声方差预估值；
步骤4、以所述平均推力测量值与所述推力噪声方差预估值为基准，通过平均推力、推力噪声方差的增大或减小搜索方法，结合蒙特卡洛模拟抽样方法，获得一系列不同平均推力、不同推力噪声方差的模拟推力；
根据一系列的模拟推力，采用辛普森数值离散化方法，利用推力积分方程，计算一系列的模拟推力作用下系统响应；
将一系列的模拟推力作用下系统响应，与待测推力作用下系统响应测量值进行逐一对比，若模拟推力作用下系统响应落入待测推力作用下系统响应包络线内，以当前模拟推力对应的平均推力和推力噪声方差作为优化的结果。


2.如权利要求1所述的方法，其特征在于，所述步骤1具体包括如下步骤：
步骤101、由扭摆振动微分方程构建推力积分方程，具体为：
所述扭摆振动微分方程为：



式中，f(t)为t时刻的施加在扭摆测量系统上的推力，t为时间变量；θ(t)为t时刻扭摆测量系统的系统响应；为扭摆测量系统的系统响应的二阶导数；ζ为阻尼比，ωn为扭摆测量系统的固有振动频率；J为转动惯量；Lf为力臂；t为时间变量；
由此构建推力积分方程，为



其中扭摆测量系统的固有振动频率为k为扭摆测量系统的扭转刚度系数；ωn(t-τ)为扭摆测量系统在t-τ时刻的固有振动频率；扭摆测量系统的阻尼比为c为扭摆测量系统的阻尼系数；扭摆测量系统的振动频率为ωd(t-τ)为扭摆测量系统在t-τ时刻的振动频率；
步骤102、对所述推力积分方程进行辛普森离散化为推力线性方程组，具体包括如下步骤：
步骤1020、在无测量噪声干扰的理想情况下，推力f(t)作用下，扭摆的系统响应为θ(t)；由于环境噪声干扰及位移传感器噪声干扰，实际系统响应测量值为Θ(t)＝θ(t)+Δθ(t)，Δθ(t)为系统响应噪声干扰；
已知系统响应测量值为[ti,Θ(ti)]，采样时间步长为h，ti为第i个采样时间取值为：i＝0,1,2,…,n；ti＝ih，当t0＝0时Θ(t0)＝0，令Θi＝Θ(ti)；
当ti＝ih时，推力积分方程为



式中，Cf为指代值，用于指代Jωd/Lf；ωd(ih-τ)为扭摆测量系统在第ih-τ时刻的振动频率；ωn(ih-τ)为扭摆测量系统在第ih-τ时刻的固有振动频率；被积分函数F(τ)表示τ时刻测量得到的推力，与待测推力f(t)有区别；
将时间区间[t0＝0,ti＝ih](i＝1,2,…)划分为i等分，第j个节点为τj＝jh，j取值为j＝0,1,2,…,i；令Fj＝F(τj)表示τj时刻测量得到的推力，推力积分方程辛普森离散化为推力线性方程组过程，具体为：
步骤1021：当i＝1时，采用梯形公式起步计算，为



可得10号系数a10：a10F0＝CfΘ1，
步骤1022：当i＝2时，采用3点辛普森公式计算，为



可得20号系数21号系数
a20F0+a21F1＝CfΘ2；
步骤1023：当i＝3时，采用4点辛普森公式计算，为



可得30号系数31号系数以及32号系数
a30F0+a31F1+a32F2＝CfΘ3
步骤1024：当i＝2k；k＝2,3,…时，采用3点辛普森公式计算，为



可得i0号系数ai0以及il号系数ail：












步骤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点辛普森公式计算，为



得到：


















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



表示为矩阵形式，为



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



式中，矩阵为下三角矩阵[aij]的逆阵，也是下三角矩阵，采用下三角矩阵快速求逆方法计算。


3.如权利要求1或2所述的方法，其特征在于，所述步骤2中，针对所述推力测量值进行线性拟合，寻找推力测量值的平均位置直线，给出平均推力测量值，具体为：
已知推力测量值为(ti,Fi)；i＝0,1,2,…,n-1，采用线性拟合方法求解平均推力，设推力测量值满足
Fi＝α+βti+εi,(i＝0,1,2,…,n-1)(24)
式中，α为第一待定系数；β为第二待定系数，εi为独立等方差的正态随机变量；
根据最小二乘方法令残差的平方和J最小



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



式中，符号“—”表示样本均值，为时间的样本均值；为推力测量的样本均值，为






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

...

【专利技术属性】
技术研发人员：李得天，洪延姬，金星，郭宁，陈新伟，
申请(专利权)人：兰州空间技术物理研究所，
类型：发明
国别省市：甘肃;62