CN109960831A - 用于扭摆系统的微推力平滑降噪优化还原方法 - Google Patents

用于扭摆系统的微推力平滑降噪优化还原方法 Download PDF

Info

Publication number
CN109960831A
CN109960831A CN201711408273.XA CN201711408273A CN109960831A CN 109960831 A CN109960831 A CN 109960831A CN 201711408273 A CN201711408273 A CN 201711408273A CN 109960831 A CN109960831 A CN 109960831A
Authority
CN
China
Prior art keywords
thrust
system response
value
slip
micro
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN201711408273.XA
Other languages
English (en)
Other versions
CN109960831B (zh
Inventor
金星
常浩
周伟静
叶继飞
李南雷
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Peoples Liberation Army Strategic Support Force Aerospace Engineering University
Original Assignee
Peoples Liberation Army Strategic Support Force Aerospace Engineering University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Peoples Liberation Army Strategic Support Force Aerospace Engineering University filed Critical Peoples Liberation Army Strategic Support Force Aerospace Engineering University
Priority to CN201711408273.XA priority Critical patent/CN109960831B/zh
Publication of CN109960831A publication Critical patent/CN109960831A/zh
Application granted granted Critical
Publication of CN109960831B publication Critical patent/CN109960831B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G06F30/15—Vehicle, aircraft or watercraft design
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/06—Power analysis or power optimisation
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Data Mining & Analysis (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Evolutionary Computation (AREA)
  • Computer Hardware Design (AREA)
  • Operations Research (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Probability & Statistics with Applications (AREA)
  • Evolutionary Biology (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Automation & Control Theory (AREA)
  • Aviation & Aerospace Engineering (AREA)
  • Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)

Abstract

本发明公开了一种用于扭摆系统的微推力平滑降噪优化还原方法，实现微推力曲线的精确还原和误差评定。该方法综合运用正交多项式局部滑动拟合方法和推力反向计算的复化辛普森方法，从计算过程、计算方法和计算结果上减小推力误差，获取一系列推力初始估计值；再综合运用样条函数插值法和系统响应预测值的正向计算方法，通过系统响应预测值与实际系统响应测量值比较，确定推力最优估计值和上下包络线。该方法能够高精度获取微推力随时间变化曲线，并有效评估推力测量误差。

Description

用于扭摆系统的微推力平滑降噪优化还原方法
技术领域
本发明涉及航天器微推进领域的推力测量技术。
背景技术
推力是微推力器的重要推进性能指标，对工况参数选择、工况优选、推力器设计以及推力器应用都有很强的工程指导意义。目前，微推力器的推力一般采用扭摆二阶振动结构测量。如果推力的作用时间小于扭摆测量系统固有周期的0.25倍时，推力作用在扭摆系统上的效应相当于瞬间脉冲冲量作用的效应，只能根据测量系统的响应辨识所施加推力的冲量大小，冲量辨识技术已经较为成熟。如果推力的作用时间大于扭摆测量系统固有周期的0.25倍时，可以根据测量系统的响应辨识推力大小。通常，推力辨识方法是根据推力与扭摆系统稳态变化量之间的线性关系。该方法只能获取推力器试验条件下的单点推力值，无法还原推力随时间的变化情况，从而无法分析推力器开机阶段、平稳工作阶段和关机阶段的推进性能。同时，微推力器的推力测量也具有很大的挑战性，易受多种噪声影响，如测量环境的振动、测量环境的气流、位移传感器等。因此，如何在多噪声条件下辨识还原推力曲线是亟待解决的技术问题。
发明内容
本发明要解决的技术问题在于：针对多噪声影响下的微推力测量和辨识还原推力曲线的问题，提出一种平滑降噪优化还原方法。该方法综合运用正交多项式局部滑动拟合方法、推力反向计算的复化辛普森方法、系统响应预测值的正向计算方法三个核心方法，试图从计算方法、计算过程和计算结果上减小推力误差，获取推力估计值，并采用系统响应预测值与实际系统响应测量值比较方法，确定最佳的推力估计值和上下包络线。
具体技术方案为：
1、微推力平滑降噪优化还原方法的实施步骤
第一步，系统参数标定和校验：采用高精度阶跃标定力产生扭摆系统的阶跃响应；利用正交多项式局部滑动拟合方法确定系统响应的平均位置曲线；根据系统响应的平均位置曲线确定振动频率和阻尼比的初始估计值及置信区间；在振动频率和阻尼比置信区间内对系统参数进行微调和验校，确定最优振动频率、最优阻尼比和最后扭转刚度系数。
第二步，确定最小划分区间数目：采用复化辛普森积分方法，将推力和系统响应的积分方程离散化为线性方程组，反向计算推力；采用与待测推力变化范围和趋势相似的模拟推力，在测量噪声忽略不计的条件下，确定积分方程离散化的最小划分区间数目。
第三步，确定推力最优估计值：采用正交多项式局部滑动拟合方法，确定推力作用下实际系统响应的平均位置曲线和上下包络线；从最小划分区间数目开始，逐步增大划分区间数目以及局部滑动拟合函数，综合运用多次局部滑动拟合方法、积分方程离散化反向计算推力、样条函数插值获得推力连续函数和微分方程正向计算系统响应预测值等方法，产生一系列推力估计值与对应的一系列系统响应预测值；根据系统响应预测值是否在实际系统响应测量值的包络线以内，以及是否在实际系统响应测量值的平均位置曲线附近上下波动最小，确定推力最优估计值；对推力最优估计值，采用正交多项式局部滑动拟合方法，确定推力最优估计值的平均位置曲线和上下包络线。
2、系统参数标定和校验的步骤
第一步，给扭摆系统施加阶跃标定力f0，获得实际系统响应测量值Θ(tj)，j＝1,2,…,m，采用正交多项式局部滑动拟合方法对实际系统响应Θ(tj)平滑降噪，平均位置曲线估计值并判读出极值点对应时间tMi和扭转角θ(tMi)，i＝1,2,…,n。
第二步，确定振动频率的初始估计值和标准差
式中，
第三步，确定阻尼比的初始估计值和标准差
式中， 利用方程获得，i＝2,3,…,n。
第四步，系统参数微调和验校：首先，在振动频率的置信区间和阻尼比的置信区间内，采样小步长搜索方法，不断改变振动频率和阻尼比；其次，在相应的振动频率和阻尼比下，计算出一系列稳态扭转角估计值 以及一系列平均位置曲线最后，构造残差根据准则“残差为零均值、零均值附近对称分布、取正值和负值的次数基本相等”，寻找到最优振动频率最优阻尼比最优稳态扭转角估计值
第五步，确定扭转刚度系数最优估计值和误差
式中，f0和df0为标定力及其误差，Lf和dLf为标定力力臂及其误差，
3、确定最小划分区间数目的步骤
第一步，获取待测推力作用下的实际系统响应测量值ΘF(ti)(i＝1,2,…,m)，采用正交多项式局部滑动拟合方法获得实际系统响应的平均位置曲线测量值θF(ti)(i＝1,2,…,m)。
第二步，采用推力反向计算的复化辛普森方法计算推力值，通过改变划分区间数目M，分析待测推力变化范围和趋势。
第三步，根据被测推力变化范围和趋势，选择模拟推力f(τ)(连续函数表示)，在测量噪声忽略不计条件下，根据推力微分方程正向计算模拟推力作用下的系统响应θfi(i＝0,1,2,…,m)。
第四步，根据模拟推力下的系统响应θfi(i＝0,1,2,…,m)，采用推力反向计算的复化辛普森方法，计算推力值fi'(i＝0,1,…,m-1)。
第五步，根据推力相对截断误差Δfi/fi＝(fi'-fi)/fi(i＝0,1,…,m-1)，选择最小划分区间数目Mmin，将推力截断误差控制在要求的水平。
4、确定推力最优估计值的步骤
第一步，根据待测推力作用下的实际系统响应测量值 [ti,ΘF(ti)](i＝1,2,…,m)，采用正交多项式局部滑动拟合方法，确定实际系统响应的平均位置曲线[ti,θF(ti)](i＝1,2,…,m)，以及确定实际系统响应测量值的上包络线ΘFu(ti)＝ΘF(ti)[ΘF(ti)＞θF(ti)]和下包络线ΘFd(ti)＝ΘF(ti)[ΘF(ti)＜θF(ti)]。
第二步，从最小划分区间数目Mmin开始，逐步增大划分区间数目，以及采用不同局部数据窗长度(如3、5、7等)的二阶正交多项式滑动拟合方法，在任一划分区间数目和局部滑动拟合函数下，首先，根据实际系统响应的平均位置曲线数据[ti,θF(ti)](i＝1,2,…,m)，采用推力反向计算的复化辛普森方法，计算推力值其次，采用正交多项式局部滑动拟合方法，对计算得到的推力值平滑降噪滤波处理，确定推力估计值为最后，采用对推力估计值样条函数插值方法，获得推力的连续函数，采用系统响应预测值的正向计算方法，计算系统响应预测值[ti,θ'F(ti)](i＝1,2,…,m)；从而得到一系列推力估计值以及对应的一系列系统响应预测值[θ'F(ti)]j(i＝1,2,…,m；j＝Mmin,Mmin+1,…)。
第三步，采用系统响应预测值与实际系统响应测量值比较方法，选择推力最优估计值选择原则一为：首先，真实推力应使系统响应预测值包含在实际系统响应测量值的上下包络线之内，即ΘFd(ti)≤[θ'F(ti)]j≤ΘFu(ti)；选择原则二为：真实推力应使系统响应预测值在实际系统响应的平均位置曲线附近上下波动最小，即残差δi＝[θ'F(ti)]j-θF(ti)满足ΘFd(ti)-θF(ti)≤δi≤ΘFu(ti)-θF(ti)，且残差为零均值、在零均值附近波动最小及对称分布。
第四步，对推力最优估计值采用正交多项式局部滑动拟合方法，确定推力最优估计值的平均位置曲线以及上包络线和下包络线
5、正交多项式局部滑动拟合方法
正交多项式局部滑动拟合方法的基本思想为：1)以每个数据点为中心对称选择若干个数据点，作为数据窗，根据数据窗内数据点，利用正交多项式，采用最小二乘方法进行拟合；2)通过数据窗逐点局部滑动过程，实现整个时间区间上数据点的拟合，寻找整个时间区间上最佳拟合曲线。
每一个数据窗下的正交多项式局部拟合方法的实施步骤为：
第一步，对于试验数据点(ti,yi)(i＝0,1,2,…,m)，以拟合点ti为中心截取(ti-p,yi-p)，(ti-p+1,yi-p+1)，…，(ti,yi)，…，(ti+p-1,yi+p-1)，(ti+p,yi+p)等2p+1 个数据点(p＝1,2,…)，数据窗长度为L＝2p+1。
第二步，对数据窗内数据点构建正交多项式拟合曲线 P(t)＝a0p0(t)+a1p1(t)+…+anpn(t)(n为拟合曲线阶数)，正交多项式 pk(t)(k＝0,1,…,n)的递推公式为
其中， 系数aj(j＝0,1,…,n)根据线性方程组求解。
6、推力反向计算的复化辛普森方法
推力反向计算的复化辛普森方法的基本思想为：采用复化辛普森公式，将待测推力F(t)与系统响应Θ(t)的积分方程离散化为线性方程组，根据推力作用下的实际系统响应测量值为[ti,Θ(ti)](i＝0,1,2,…,m)递推计算出推力值[ti,F(ti)](i＝0,1,2,…,m)。
若实际系统响应测量值采样时间步长为h，设ti＝ih(i＝0,1,2,…,m)；设 t0＝0时，Θ(t0)＝0；令Θ(ti)＝Θi和F(ti)＝Fi；采用复化辛普森公式将推力积分方程离散化为线性方程组过程如下：
当i＝1时，a10F0＝CfΘ1，
当i＝2时，a20F0+a21F1＝CfΘ2，
当i＝3时，a30F0+a31F1+a32F2＝CfΘ3，
当i＝2k(k＝2,3,...)时，
当i＝2k+1(k＝2,3,...)时，
7、系统响应预测值的微分方程正向计算方法
系统响应预测值的微分方程正向计算方法的基本思想为：根据推力估计值，采用样条函数插值方法获得推力的连续函数S(t)，再根据振动微分方程计算系统响应预测值 [ti,θp(ti)](i＝1,2,…,m)，其中，t0＝0，θp(0)＝0。
已知推力估计值设将推力作用时间区间[0,T0]划分为n′＝m-1个小区间(节点数目为n′+1)，有 0＝t0＜t1＜t2＜…＜tn′-1＜tn′＝T0，小区间长度为hi-1＝ti-ti-1(1≤i≤n′)，在子区间[ti-1,ti](1≤i≤n′)上，三次样条函数为 S(t)＝Fi-1+(Fi-Fi-1)(t-ti-1)/(ti-ti-1)+[ai(t-ti-1)+bi(t-ti)](t-ti-1)(t-ti)，ai和bi为待定系数。
ai和bi的计算表达式为：
式中，Mi可由计算获得，其中，λi＝1-μi，
与现有技术相比，采用本发明可以达到以下技术效果：
1.本发明能够还原出待测推力的最佳推力曲线以及推力曲线的上下包络线，这是区别已有技术的最大特色；
2.本发明具有较好的通用性，适用于类似扭摆系统的二阶振动测量系统的推力曲线还原。
附图说明
图1为本发明的技术路线图。
具体实施方案
结合附图对用于扭摆系统的微推力平滑降噪优化还原方法做进一步详细描述。微推力平滑降噪优化方法的基本思路是：①采用实际系统响应测量值平滑降噪方法，从整个计算过程上减小测量噪声影响；采用高精度数值积分离散化方法，从计算方法上减小推力误差；采用推力计算值平滑降噪方法，从计算结果上减小推力噪声误差，最终获得推力估计值。②采用对推力估计值样条函数插值方法，获得推力的连续函数；采用振动微分方程正向计算方法，获得系统响应预测值；采用系统响应预测值与实际系统响应测量值比较方法，验证推力估计值的优化程度，在一系列推力估计值中确定最佳的推力估计值。
微推力平滑降噪优化还原方法的步骤概述为：
第一步，系统参数标定和校验：采用高精度阶跃标定力产生扭摆系统的阶跃响应；利用正交多项式局部滑动拟合方法确定系统响应的平均位置曲线；根据系统响应的平均位置曲线确定振动频率和阻尼比的初始估计值及置信区间；在振动频率和阻尼比置信区间内对系统参数进行微调和验校，确定最优振动频率、最优阻尼比和最后扭转刚度系数；
第二步，确定最小划分区间数目：采用复化辛普森积分方法，将推力和系统响应的积分方程离散化为线性方程组，反向计算推力；采用与待测推力变化范围和趋势相似的模拟推力，在测量噪声忽略不计的条件下，确定积分方程离散化的最小划分区间数目；
第三步，确定推力最优估计值：采用正交多项式局部滑动拟合方法，确定推力作用下实际系统响应的平均位置曲线和上下包络线；从最小划分区间数目开始，逐步增大划分区间数目以及局部滑动拟合函数，综合运用多次局部滑动拟合方法、积分方程离散化反向计算推力、样条函数插值获得推力连续函数和微分方程正向计算系统响应预测值等方法，产生一系列推力估计值与对应的一系列系统响应预测值；根据系统响应预测值是否在实际系统响应测量值的包络线以内，以及是否在实际系统响应测量值的平均位置曲线附近上下波动最小，确定推力最优估计值；对推力最优估计值，采用正交多项式局部滑动拟合方法，确定推力最优估计值的平均位置曲线和上下包络线。
图1为本发明的技术路线图，结合附图具体阐述微推力平滑降噪优化还原方法的推导原理和实施步骤。
1、系统参数标定和校验
(1)采样阶跃响应方法，利用高精度标定力，标定系统参数
在已知阶跃标定力f0作用下，扭摆测量系统的实际系统响应为
式中，稳态扭转角θ(∞)是未知量；振动频率ωd、阻尼比ζ和扭转刚度系数 k等是未知量；Δθ(t)～N(0,σ2)为标定过程中环境、位移传感器和标定力等产生的测量噪声，测量噪声方差σ2是未知量。
(2)标定振动频率和阻尼比
由于测量噪声Δθ(t)～N(0,σ2)的影响，实际系统响应Θ(t)在平均位置曲线θ(t)附近上下起伏波动，给极值点对应时间和扭转角判读带来困难，因此，采样正交多项式局部滑动拟合方法，对实际系统响应Θ(t)平滑降噪滤波处理后，再判读极值点对应时间和tMi扭转角θ(tMi)。
振动频率为
振动频率的估计值和标准差为
极值点对应扭转角为
可得
令
式中，0≤xζ≤1。得到xζ的方程
可解得阻尼比为
阻尼比的估计值和标准差为
(3)标定转刚度系数、稳态扭转角、测量噪声方差
在已知振动频率估计值阻尼比估计值条件下，设系统响应的平均位置曲线为
式中，稳态扭转角θ(∞)是待估计量。
根据实际系统响应测量值Θ(ti)(i＝1,2,…,m)和平均位置曲线估计值构造残差为采用最小二乘方法，使得残差的平方和最小，确定稳态扭转角θ(∞)的估计值
残差为
令其平方和为最小，即
可得稳态扭转角θ(∞)的估计值为
其方差为
通过稳态扭转角θ(∞)的估计值计算扭转刚度系数的估计值，为
其绝对误差和相对误差，分别为
式中，df0为标定力的误差，dLf为力臂的误差。
测量噪声Δθ(t)～N(0,σ2)方差σ2的估计值，为
从而得到平均位置曲线为
确定平均位置曲线的过程为：首先确定振动频率和阻尼比的估计值，其次采用最小二乘方法确定稳态扭转角。
(4)系统参数的验校和微调
振动频率估计值阻尼比估计值和稳态扭转角估计值等，是否为最佳估计值需要检验和校正。验校的准则为：最佳的振动频率估计值阻尼比估计值和稳态扭转角估计值等，应使得残差为零均值、零均值附近对称分布、取正值和负值的次数基本相等。
首先，在振动频率的置信区间和阻尼比的置信区间内，采样小步长搜索方法，不断改变振动频率和阻尼比。
其次，采用步骤③的方法，计算稳态扭转角估计值并确定平均位置曲线
最后，构造残差根据准则“残差为零均值、零均值附近对称分布、取正值和负值的次数基本相等”，寻找到最佳振动频率、阻尼比、稳态扭转角、扭转刚度系数。
2、确定最小划分区间数目
(1)根据被测推力特点，选择模拟推力
首先，根据实际系统响应测量值Θ(ti)(i＝1,2,…,m)，采用正交多项式局部滑动拟合方法，计算实际系统响应的平均位置曲线测量值θ(ti)(i＝1,2,…,m)。
其次，采用推力反向计算的复化辛普森方法，计算推力值，通过改变划分区间数目m，考察被测推力变化范围和趋势。
最后，根据被测推力变化范围和趋势，采用相似的模拟推力表示被测推力变化范围和趋势，模拟推力采用连续函数表示。
(2)根据模拟推力，选择最小划分区间数目
首先，采用模拟推力f(τ)(连续函数表示)，根据推力微分方程正向计算系统响应θi(i＝0,1,2,…,m)；
其次，在测量噪声忽略不计条件下，根据模拟推力下的系统响应θi(i＝0,1,2,…,m)，采用推力反向计算的复化辛普森方法，计算推力值 Fi(i＝0,1,…,m-1)。
最后，根据推力相对截断误差Δfi/fi＝(Fi-fi)/fi(i＝0,1,…,m-1)，选择最小划分区间数目mmin，将推力截断误差控制在要求的水平。
在推力作用时间区间[0,T0]内，合理选择最小划分区间数目mmin，将推力截断误差控制在要求范围内，同时由于划分区间数目最小，达到降低推力噪声误差的目的。例如，将推力的截断误差控制在0.1％～0.5％范围内，选择最小划分区间数目。
3、推力反算和响应预测
(1)根据实际系统响应测量值[ti,Θ(ti)]，采用正交多项式局部滑动拟合方法，确定系统响应的平均位置曲线达到从整个计算过程上减小测量噪声影响的目的。
(2)根据实际系统响应的平均位置曲线数据采用推力反向计算的复化辛普森方法，计算推力值
式中，达到采用高精度3点和4点交替的复化辛普森离散化方法，减小划分区间数目，降低推力误差的目的。
(3)采用正交多项式局部滑动拟合方法，对计算得到的推力值平滑降噪滤波处理，确定推力估计值为
从而达到从计算结果上再次减小推力噪声误差的目的，最终获得推力估计值。
(4)采用对推力估计值样条函数插值方法，获得推力的连续函数，采用系统响应预测值的正向计算方法，计算系统响应预测值。
首先，已知推力的估计值设将推力作用时间区间[0,T0]划分为n′＝m-1个小区间(节点数目为n′+1)，有
0＝t0＜t1＜t2＜…＜tn′-1＜tn′＝T0
小区间长度为hi-1＝ti-ti-1(1≤i≤n′)，在子区间[ti-1,ti](1≤i≤n′)上，推力估计值的三次样条函数为
S(t)＝Fi-1+F(ti-1,ti)(t-ti-1) +[ai(t-ti-1)+bi(t-ti)](t-ti-1)(t-ti)
其次，利用微分方程数值计算可方便地获得高精度数值解特点，采用振动微分方程正向计算方法，计算系统响应预测值。扭摆振动微分方程表示为
式中，S(t)为推力估计值的样条插值连续函数。
已知推力估计值采用系统响应预测值的正向计算方法，计算的系统响应预测值为[ti,θp(ti)](i＝1,2,…,m)。
因此，针对振动微分方程正向计算高精度特点、但要求推力是连续函数特点，实现了由推力估计值计算系统响应预测值。
4、确定推力最优估计值
(1)采用正交多项式局部滑动拟合方法，确定实际系统响应的平均位置曲线、测量噪声方差、实际系统响应的上下包络线
首先，根据实际系统响应测量值[ti,Θ(ti)]，采用正交多项式局部滑动拟合方法，确定系统响应的平均位置曲线
其次，通过残差分析，确定测量噪声的方差。测量噪声Δθ(t)～N(0,σ2)方差σ2的估计值，为
最后，确定实际系统响应测量值的上下包络线，上包络线和下包络线，分别为
(2)从最小划分区间数目开始，逐步增大划分区间数目，以及选择不同平滑降噪滤波函数，综合运用上述推力反算和响应预测的方法，产生一系列推力估计值和对应的一系列系统响应预测值。
首先，选择正交多项式局部滑动拟合方法，选择局部拟合能力突出的二阶正交多项式(抛物线正交多项式)，选择局部数据窗长度为 2p+1＝3,5,7,…。
其次，根据实际系统响应测量值[ti,Θ(ti)]，从最小划分区间数目mmin开始，逐步增大划分区间数目，综合运用上述推力反算和响应预测的方法，得到一系列推力估计值和对应的一系列系统响应预测值[θp(ti)]j，为
[θp(ti)]j(i＝1,2,…,m；j＝mmin,mmin+1,…)
(3)采用系统响应预测值与实际系统响应测量值比较方法，确定最佳的推力估计值。
首先，针对真实推力应使系统响应预测值包含在实际系统响应测量值的上下包络线之内的特点，使得Θd(ti)≤[θp(ti)]j≤Θu(ti)。
其次，针对真实推力应使系统响应预测值在实际系统响应的平均位置曲线附近上下波动最小的特点，考察残差随着时间变化，当残差满足且残差为零均值、在零均值附近波动最小及对称分布时，寻找到最佳的推力估计值。
最后，对最佳的推力估计值，采用正交多项式局部滑动拟合方法，平滑降噪滤波处理，确定推力估计值和推力误差。
因此，通过系统响应预测值最佳逼近实际系统响应的平均位置曲线，达到了寻找最佳推力估计值的目的、以及推力测量和误差分析的目的。
从上述平滑降噪优化还原方法的具体应用过程可以看出，涉及正交多项式局部滑动拟合方法、推力反向计算的复化辛普森方法和系统响应预测值的正向计算方法三个核心方法，下面分别阐述。
1、正交多项式局部滑动拟合方法
(1)实际系统响应测量值与平均位置曲线
实际系统响应测量值总是包含测量噪声，实际系统响应测量值为
Θ(t)＝θ(t)+Δθ(t)
式中，θ(t)为理想系统响应，Δθ(t)为测量噪声。
根据实际系统响应测量值[ti,Θ(ti)](i＝0,1,2,…,m)，采用正交多项式局部滑动拟合方法，经过平滑降噪滤波处理，确定平均位置曲线作为理想系统响应θ(t)的估计值。
(2)正交多项式局部滑动拟合方法
①正交多项式拟合方法
实际系统响应测量值为(ti,Θi)(i＝0,1,2,…,m)(Θi＝Θ(ti))，由已知线性无关函数构造多项式P(t)进行拟合，为
根据最小二乘法，可得线性方程组
通过该线性方程组可解得系数aj(j＝0,1,…,n)。
对于已知线性无关函数上述方程组存在唯一解，但是该线性方程组往往具有病态特性，导致系数的计算误差较大，因此，选取正交多项式进行拟合。
②正交多项式局部滑动拟合方法
根据实际系统响应测量值(ti,Θi)(i＝0,1,2,…,m)，利用正交多项式P(t) 拟合实际系统响应曲线，避免了求解病态线性方程组，提高了系数 aj(j＝0,1,…,n)的求解精度。但是，由于实际系统响应是多极值点、多拐点、曲率变化大的复杂曲线，在整个时间t0≤t≤tm上全程拟合，不能反映每个局部数据取值特点，出现拟合曲线振荡和拟合误差大现象。因此，提出正交多项式局部滑动拟合方法。
正交多项式局部滑动拟合方法是：1)以每个数据点为中心对称选择若干个数据点，作为数据窗，根据数据窗内数据点，利用正交多项式，采用最小二乘方法进行拟合；2)通过数据窗逐点局部滑动过程，实现整个时间区间上数据点的拟合，寻找最佳拟合曲线，避免了全程数据点拟合的误差。
对于试验数据点(ti,Θi)(i＝0,1,2,…,m)中，以拟合点ti为中心，取 (ti-p,Θi-p)，(ti-p+1,Θi-p+1)，…，(ti,Θi)，…，(ti+p-1,Θi+p-1)，(ti+p,Θi+p) 等2p+1点(p＝1,2,…)，数据窗长度为l+1＝2p+1，进行正交多项式局部滑动拟合。
具体构造正交多项式pk(t)(k＝0,1,…,n)的递推公式为
其中
最后，采用最小二乘方法，利用正交多项式拟合曲线为
P(t)＝a0p0(t)+a1p1(t)+…+anpn(t)
其中，n为拟合曲线阶数。
2、推力反向计算的复化辛普森方法
(1)推力微分方程的正向计算和推力积分方程的反向计算
对于扭摆测量系统，推力与系统响应之间关系，采用微分方程表示为
式中，ζ为阻尼比，ωn为固有频率，J为转动惯量。
上述推力与系统响应的微分方程，也可表示为推力与系统响应的积分方程形式，为
在已知推力f(t)条件下，由推力微分方程，可方便地高精度计算系统响应，这是推力微分方程的正向计算问题；在已知系统响应θ(t)条件下，由推力积分方程，也可计算推力，这是推力积分方程的反向计算问题。
(2)推力积分方程离散化为线性方程组
采用复化辛普森公式，将以下推力积分方程离散化为线性方程组
式中，Cf＝Jωd/Lf。
设已知实际系统响应测量值为[ti,Θ(ti)](i＝0,1,2,…,m)，采样时间步长为h，ti＝ih，当t0＝0时Θ(t0)＝0，令Θ(ti)＝Θi和F(ti)＝Fi。当ti＝ih时，可知
采用复化辛普森公式，将推力积分方程离散化为线性方程组，需要交替使用3点和4点辛普森公式，分以下5种情况。
①当i＝1时，根据梯形公式，可知
可得
②当i＝2时，根据3点辛普森公式，可知
可得
a20F0+a21F1＝CfΘ2，
③当i＝3时，根据4点辛普森公式，可知
可得
a30F0+a31F1+a32F2＝CfΘ3
④当i＝2k(k＝2,3,…)时，根据3点辛普森公式，可知
可得
⑤当i＝2k+1(k＝2,3,…)时，根据4点辛普森公式，可知
式中，在点i＝2k+1＝2(k-1)+3(k＝2,3,…)中，最后的2k-2,2k-1,2k,2k+1 等4个点采用4点辛普森公式，之前的0,1,…,2k-3,2k-2等采用3点辛普森公式。可得
即推力积分方程离散化线性方程组，为
式中，aij＝0(j≥i)。
因此，根据实际系统响应测量值为(ti,Θi)(i＝0,1,2,…,m)，可反向计算推力值Fi(i＝0,1,2,…,m-1)。这种根据实际系统响应测量值，反向计算推力的方法，称为推力反向计算的复化辛普森方法。推力反向计算的复化辛普森方法，截断误差与步长h4成正比，可高精度计算推力值。
(3)推力误差包括推力截断误差和推力噪声误差
实际系统响应测量值总是包含测量噪声，实际系统响应测量值为
Θ(t)＝θ(t)+Δθ(t)
式中，θ(t)为理想系统响应，Δθ(t)为测量噪声。
在没有测量噪声Δθ(t)的理想情况下，理想系统响应为θ(t)，推力积分方程为
式中，Cf＝Jωd/Lf，f(τ)为推力真值。
通过积分方程离散化为线性方程组，将积分表示为有限和形式，为
积分计算的截断误差，造成推力计算误差Δf(τi)＝F(τi)-f(τi)，这种推力积分方程离散化时截断误差引起的推力误差，称为推力截断误差。由于推力积分方程的被积分是振荡函数，划分区间数目足够多时，才能控制推力截断误差。
实际系统响应Θ(t)＝θ(t)+Δθ(t)包含测量噪声，推力积分方程离散化为线性方程组，为
测量噪声Δθ(t)造成推力计算误差Δf(τi)＝F(τi)-f(τi)，这种测量噪声引起的推力误差，称为推力噪声误差。由于离散化线性方程组的严重病态特性和系统响应测量值包含测量噪声，划分区间数目多时，造成推力噪声误差急剧增大。
推力反向计算的难点在于推力截断误差和推力噪声误差相互制约限制。划分区间数目少时，造成推力截断误差增大；划分区间数目多时，造成推力噪声误差急剧增大。
推力反向计算的复化辛普森方法，针对积分方程离散化的划分区间数目造成推力截断误差和推力噪声误差相互制约限制的难点，采用高精度的复化辛普森(3点和4点交替运用)数值积分方法，将推力积分方程离散化为线性方程组，可在划分区间数目较少前提下，实现尽量减小推力截断误差和推力噪声误差，达到减小推力误差的目的。
3、系统响应预测值的正向计算方法
(1)根据推力估计值，采用样条函数插值方法，获得推力的连续函数
针对振动微分方程正向计算高精度特点、但要求推力是连续函数特点，由推力估计值，采用样条函数插值方法，获得推力的连续函数。
利用样条函数具有函数光滑和导数连续的特点，采用样条函数插值方法，获得推力的连续函数。已知推力估计值设将推力作用时间区间[0,T0]划分为n′＝m-1个小区间(节点数目为n′+1)，有
0＝t0＜t1＜t2＜…＜tn′-1＜tn′＝T0
小区间长度为hi-1＝ti-ti-1(1≤i≤n′)，在子区间[ti-1,ti](1≤i≤n′)上，三次样条函数为
S(t)＝Fi-1+F(ti-1,ti)(t-ti-1) +[ai(t-ti-1)+bi(t-ti)](t-ti-1)(t-ti)
式中，F(ti-1,ti)＝(Fi-Fi-1)/(ti-ti-1)为一阶均差，ai和bi为待定系数。
当i＝1,2,…,n′-1时，可得方程
μiMi-1+2Mi+λiMi+1＝di
其中
显然，插值区间等份时，μi＝λi＝1/2。
需要补充2个边界条件方程，推力的左边和右边边界看作第一类边界问题，已知F′0和F′n。补充得到两个方程
从而得到三对角方程
解得Mi(i＝0,1,…,n′)，可求得待定系数ai(i＝1,2,…,n′)和bi(i＝1,2,…,n′)，为
式中，小区间长度为hi-1＝ti-ti-1。其中F′0和F′n，可采用数值微分公式计算，例如
(2)根据推力估计值的连续函数，采用微分方程正向计算方法，获得系统响应预测值
利用微分方程正向计算可方便地获得高精度数值解特点，采用振动微分方程正向计算方法，计算系统响应预测值。
扭摆振动微分方程表示为
式中，
S(t)为推力样条插值连续函数，为
S(t)＝Fi-1+F(ti-1,ti)(t-ti-1) +[ai(t-ti-1)+bi(t-ti)](t-ti-1)(t-ti)
已知推力估计值采用振动微分方程正向计算方法，计算系统响应预测值[ti,θp(ti)](i＝1,2,…,m)，其中，t0＝0和θp(0)＝0。
在振动微分方程正向计算中，推力估计值的插值点数目过少(或积分方程离散化的划分区间数目过少)，样条插值函数不能准确表示推力随着时间变化，造成系统响应预测值的误差，因此，离散化的划分区间数目不但影响推力截断误差和推力噪声误差，还影响系统响应预测值的误差。
本申请不局限于说明书和权利要求文字部分所限定的内容，任何本领域范围内公知的修改和变化都属于本申请的范围，说明书具体实施例部分仅是对本发明示例性的说明，不是对本发明的具体限定。

Claims (7)

1.用于扭摆系统的微推力平滑降噪优化还原方法，其特征在于：
用于扭摆系统的微推力平滑降噪优化还原方法的实施步骤为：
第一步，系统参数标定和校验：采用高精度阶跃标定力产生扭摆系统的阶跃响应；利用正交多项式局部滑动拟合方法确定系统响应的平均位置曲线；根据系统响应的平均位置曲线确定振动频率和阻尼比的初始估计值及置信区间；在振动频率和阻尼比置信区间内对系统参数进行微调和验校，确定最优振动频率、最优阻尼比和最后扭转刚度系数；
第二步，确定最小划分区间数目：采用复化辛普森积分方法，将推力和系统响应的积分方程离散化为线性方程组，反向计算推力；采用与待测推力变化范围和趋势相似的模拟推力，在测量噪声忽略不计的条件下，确定积分方程离散化的最小划分区间数目；
第三步，确定推力最优估计值：采用正交多项式局部滑动拟合方法，确定推力作用下实际系统响应的平均位置曲线和上下包络线；从最小划分区间数目开始，逐步增大划分区间数目以及局部滑动拟合函数，综合运用多次局部滑动拟合方法、积分方程离散化反向计算推力、样条函数插值获得推力连续函数和微分方程正向计算系统响应预测值等方法，产生一系列推力估计值与对应的一系列系统响应预测值；根据系统响应预测值是否在实际系统响应测量值的包络线以内，以及是否在实际系统响应测量值的平均位置曲线附近上下波动最小，确定推力最优估计值；对推力最优估计值，采用正交多项式局部滑动拟合方法，确定推力最优估计值的平均位置曲线和上下包络线。
2.如权利要求1所述的用于扭摆系统的微推力平滑降噪优化还原方法，其特征在于：
系统参数标定和校验的步骤为：
第一步，给扭摆系统施加阶跃标定力f0，获得实际系统响应测量值Θ(tj)，j＝1,2,…,m，采用正交多项式局部滑动拟合方法对实际系统响应Θ(tj)平滑降噪，平均位置曲线估计值并判读出极值点对应时间tMi和扭转角θ(tMi)，i＝1,2,…,n；
第二步，确定振动频率的初始估计值和标准差
式中，
第三步，确定阻尼比的初始估计值和标准差
式中， 利用方程获得，i＝2,3,…,n；
第四步，系统参数微调和验校：首先，在振动频率的置信区间和阻尼比的置信区间内，采样小步长搜索方法，不断改变振动频率和阻尼比；其次，在相应的振动频率和阻尼比下，计算出一系列稳态扭转角估计值 以及一系列平均位置曲线最后，构造残差根据准则“残差为零均值、零均值附近对称分布、取正值和负值的次数基本相等”，寻找到最优振动频率最优阻尼比最优稳态扭转角估计值
第五步，确定扭转刚度系数最优估计值和误差
式中，f0和df0为标定力及其误差，Lf和dLf为标定力力臂及其误差，
3.如权利要求1所述的用于扭摆系统的微推力平滑降噪优化还原方法，其特征在于：
确定最小划分区间数目的步骤为：
第一步，获取待测推力作用下的实际系统响应测量值ΘF(ti)(i＝1,2,…,m)，采用正交多项式局部滑动拟合方法获得实际系统响应的平均位置曲线测量值θF(ti)(i＝1,2,…,m)；
第二步，采用推力反向计算的复化辛普森方法计算推力值，通过改变划分区间数目M，分析待测推力变化范围和趋势；
第三步，根据被测推力变化范围和趋势，选择模拟推力f(τ)(连续函数表示)，在测量噪声忽略不计条件下，根据推力微分方程正向计算模拟推力作用下的系统响应θfi(i＝0,1,2,…,m)；
第四步，根据模拟推力下的系统响应θfi(i＝0,1,2,…,m)，采用推力反向计算的复化辛普森方法，计算推力值fi'(i＝0,1,…,m-1)；
第五步，根据推力相对截断误差Δfi/fi＝(fi'-fi)/fi(i＝0,1,…,m-1)，选择最小划分区间数目Mmin，将推力截断误差控制在要求的水平。
4.如权利要求1所述的用于扭摆系统的微推力平滑降噪优化还原方法，其特征在于：
确定推力最优估计值的步骤为：
第一步，根据待测推力作用下的实际系统响应测量值[ti,ΘF(ti)](i＝1,2,…,m)，采用正交多项式局部滑动拟合方法，确定实际系统响应的平均位置曲线[ti,θF(ti)](i＝1,2,…,m)，以及确定实际系统响应测量值的上包络线ΘFu(ti)＝ΘF(ti)[ΘF(ti)＞θF(ti)]和下包络线ΘFd(ti)＝ΘF(ti)[ΘF(ti)＜θF(ti)]；
第二步，从最小划分区间数目Mmin开始，逐步增大划分区间数目，以及采用不同局部数据窗长度(如3、5、7等)的二阶正交多项式滑动拟合方法，在任一划分区间数目和局部滑动拟合函数下，首先，根据实际系统响应的平均位置曲线数据[ti,θF(ti)](i＝1,2,…,m)，采用推力反向计算的复化辛普森方法，计算推力值其次，采用正交多项式局部滑动拟合方法，对计算得到的推力值平滑降噪滤波处理，确定推力估计值为最后，采用对推力估计值样条函数插值方法，获得推力的连续函数，采用系统响应预测值的正向计算方法，计算系统响应预测值[ti,θ'F(ti)](i＝1,2,…,m)；从而得到一系列推力估计值以及对应的一系列系统响应预测值[θ'F(ti)]j(i＝1,2,…,m；j＝Mmin,Mmin+1,…)；
第三步，采用系统响应预测值与实际系统响应测量值比较方法，选择推力最优估计值选择原则一为：首先，真实推力应使系统响应预测值包含在实际系统响应测量值的上下包络线之内，即ΘFd(ti)≤[θ'F(ti)]j≤ΘFu(ti)；选择原则二为：真实推力应使系统响应预测值在实际系统响应的平均位置曲线附近上下波动最小，即残差δi＝[θ'F(ti)]j-θF(ti)满足ΘFd(ti)-θF(ti)≤δi≤ΘFu(ti)-θF(ti)，且残差为零均值、在零均值附近波动最小及对称分布；
第四步，对推力最优估计值采用正交多项式局部滑动拟合方法，确定推力最优估计值的平均位置曲线以及上包络线和下包络线
5.如权利要求1、2、3、4所述的用于扭摆系统的微推力平滑降噪优化还原方法，其特征在于：
正交多项式局部滑动拟合方法的基本思想为：1)以每个数据点为中心对称选择若干个数据点，作为数据窗，根据数据窗内数据点，利用正交多项式，采用最小二乘方法进行拟合；2)通过数据窗逐点局部滑动过程，实现整个时间区间上数据点的拟合，寻找整个时间区间上最佳拟合曲线；
每一个数据窗下的正交多项式局部拟合方法的实施步骤为：
第一步，对于试验数据点(ti,yi)(i＝0,1,2,…,m)，以拟合点ti为中心截取(ti-p,yi-p)，(ti-p+1,yi-p+1)，…，(ti,yi)，…，(ti+p-1,yi+p-1)，(ti+p,yi+p)等2p+1个数据点(p＝1,2,…)，数据窗长度为L＝2p+1；
第二步，对数据窗内数据点构建正交多项式拟合曲线P(t)＝a0p0(t)+a1p1(t)+…+anpn(t)(n为拟合曲线阶数)，正交多项式pk(t)(k＝0,1,…,n)的递推公式为
其中， 系数aj(j＝0,1,…,n)根据线性方程组求解。
6.如权利要求1、3、4所述的用于扭摆系统的微推力平滑降噪优化还原方法，其特征在于：
推力反向计算的复化辛普森方法的基本思想为：采用复化辛普森公式，将待测推力F(t)与系统响应Θ(t)的积分方程离散化为线性方程组，根据推力作用下的实际系统响应测量值为[ti,Θ(ti)](i＝0,1,2,…,m)递推计算出推力值[ti,F(ti)](i＝0,1,2,…,m)；
若实际系统响应测量值采样时间步长为h，设ti＝ih(i＝0,1,2,…,m)；设t0＝0时，Θ(t0)＝0；令Θ(ti)＝Θi和F(ti)＝Fi；采用复化辛普森公式将推力积分方程离散化为线性方程组过程如下：
当i＝1时，a10F0＝CfΘ1，
当i＝2时，a20F0+a21F1＝CfΘ2，
当i＝3时，a30F0+a31F1+a32F2＝CfΘ3，
当i＝2k(k＝2,3,...)时，
当i＝2k+1(k＝2,3,...)时，
7.如权利要求1、3、4所述的用于扭摆系统的微推力平滑降噪优化还原方法，其特征在于：
系统响应预测值的微分方程正向计算方法的基本思想为：根据推力估计值，采用样条函数插值方法获得推力的连续函数S(t)，再根据振动微分方程计算系统响应预测值[ti,θp(ti)](i＝1,2,…,m)，其中，t0＝0，θp(0)＝0；
已知推力估计值设将推力作用时间区间[0,T0]划分为n′＝m-1个小区间(节点数目为n′+1)，有0＝t0＜t1＜t2＜…＜tn′-1＜tn′＝T0，小区间长度为hi-1＝ti-ti-1(1≤i≤n′)，在子区间[ti-1,ti](1≤i≤n′)上，三次样条函数为S(t)＝Fi-1+(Fi-Fi-1)(t-ti-1)/(ti-ti-1)+[ai(t-ti-1)+bi(t-ti)](t-ti-1)(t-ti)，ai和bi为待定系数；
ai和bi的计算表达式为：
式中，Mi可由计算获得，其中，λi＝1-μi，
CN201711408273.XA 2017-12-22 2017-12-22 用于扭摆系统的微推力平滑降噪优化还原方法 Active CN109960831B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201711408273.XA CN109960831B (zh) 2017-12-22 2017-12-22 用于扭摆系统的微推力平滑降噪优化还原方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201711408273.XA CN109960831B (zh) 2017-12-22 2017-12-22 用于扭摆系统的微推力平滑降噪优化还原方法

Publications (2)

Publication Number Publication Date
CN109960831A true CN109960831A (zh) 2019-07-02
CN109960831B CN109960831B (zh) 2023-02-28

Family

ID=67019691

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201711408273.XA Active CN109960831B (zh) 2017-12-22 2017-12-22 用于扭摆系统的微推力平滑降噪优化还原方法

Country Status (1)

Country Link
CN (1) CN109960831B (zh)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112231890A (zh) * 2020-09-03 2021-01-15 兰州空间技术物理研究所 基于扭摆测量系统的高平稳电推力器的推力测评方法
CN112800126A (zh) * 2021-01-13 2021-05-14 海南微氪生物科技股份有限公司 荧光光电检测仪器关于预测检测时间的处理方法及系统
CN114545863A (zh) * 2022-03-07 2022-05-27 中南大学 一种基于b样条曲线拟合的数控加工的轨迹平滑方法

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101514927A (zh) * 2009-03-20 2009-08-26 北京航空航天大学 弹性微牛级小推力测量系统
CN103335769A (zh) * 2013-07-03 2013-10-02 中国科学院力学研究所 一种电推进器弱力测量装置
US20130321818A1 (en) * 2012-05-29 2013-12-05 General Photonics Corporation Measuring polarization crosstalk in optical birefringent materials and devices based on reduction of line broadening caused by birefringent dispersion
CN104280168A (zh) * 2014-04-24 2015-01-14 北京航空航天大学 一种基于双光束干涉原理的高精度光学微小推力测量系统
CN107091705A (zh) * 2017-05-22 2017-08-25 河南理工大学 一种微推力测量方法及装置

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101514927A (zh) * 2009-03-20 2009-08-26 北京航空航天大学 弹性微牛级小推力测量系统
US20130321818A1 (en) * 2012-05-29 2013-12-05 General Photonics Corporation Measuring polarization crosstalk in optical birefringent materials and devices based on reduction of line broadening caused by birefringent dispersion
CN103335769A (zh) * 2013-07-03 2013-10-02 中国科学院力学研究所 一种电推进器弱力测量装置
CN104280168A (zh) * 2014-04-24 2015-01-14 北京航空航天大学 一种基于双光束干涉原理的高精度光学微小推力测量系统
CN107091705A (zh) * 2017-05-22 2017-08-25 河南理工大学 一种微推力测量方法及装置

Non-Patent Citations (6)

* Cited by examiner, † Cited by third party
Title
常浩等: "脉冲微推力作用下的扭摆系统响应特性仿真分析", 《红外与激光工程》 *
洪延姬等: "微推力测量方法及其关键问题分析", 《航空学报》 *
潘小青等: "扭摆法测量转动惯量实验数据处理及误差讨论", 《物理与工程》 *
罗乐乐等: "基于扭摆测量系统的微冲量测量方法", 《机电产品开发与创新》 *
金星等: "一种用于微小推力冲量测量的扭摆系统参数标定方法", 《推进技术》 *
金星等: "扭摆推力测量系统的量程设计方法", 《测试技术学报》 *

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112231890A (zh) * 2020-09-03 2021-01-15 兰州空间技术物理研究所 基于扭摆测量系统的高平稳电推力器的推力测评方法
CN112800126A (zh) * 2021-01-13 2021-05-14 海南微氪生物科技股份有限公司 荧光光电检测仪器关于预测检测时间的处理方法及系统
CN112800126B (zh) * 2021-01-13 2022-11-15 海南微氪生物科技股份有限公司 荧光光电检测仪器关于预测检测时间的处理方法及系统
CN114545863A (zh) * 2022-03-07 2022-05-27 中南大学 一种基于b样条曲线拟合的数控加工的轨迹平滑方法
CN114545863B (zh) * 2022-03-07 2024-02-13 中南大学 一种基于b样条曲线拟合的数控加工的轨迹平滑方法

Also Published As

Publication number Publication date
CN109960831B (zh) 2023-02-28

Similar Documents

Publication Publication Date Title
CN109960831A (zh) 用于扭摆系统的微推力平滑降噪优化还原方法
Djurovic et al. A hybrid CPF-HAF estimation of polynomial-phase signals: Detailed statistical analysis
KR102206324B1 (ko) 기술적 시스템의 스타팅 변수의 모델을 확인하기 위한 방법
CN106385211B (zh) 一种步进电机负载转矩估计方法
CN108470089B (zh) 一种基于最小二乘样本拟合的复信号时延估计方法
CN101262199A (zh) 电动机控制装置以及电动机控制系统
CN109977613B (zh) 一种可预先设定调整时间的自适应滑模末制导律设计方法
CN104112079A (zh) 一种模糊自适应变分贝叶斯无迹卡尔曼滤波方法
CN101860344A (zh) 一种选频滤波器的构建方法及采用该方法实现fir型、iir型滤波器的构建方法
CN109298351A (zh) 一种基于模型学习的新能源车载电池剩余寿命估计方法
CN102779238A (zh) 一种基于自适应卡尔曼滤波的无刷直流电机系统辨识方法
US20230150696A1 (en) Digital filter based method for measuring thrust response time of satellite-borne micro-thruster
CN103593519B (zh) 一种基于试验设计的运载火箭总体参数优化方法
Huang et al. Attitude determination method integrating square-root cubature Kalman filter with expectation-maximization for inertial navigation system applied to underwater glider
Li et al. Improved wavelet threshold denoising method for MEMS gyroscope
CN107121641B (zh) 一种基于粒子群优化的电池状态估计方法
CN114235072B (zh) 一种基于过零检测的科氏流量计相位差计算方法
CN109409194B (zh) 多模态时域信号模态分离、阻尼参数辨识方法及存储介质
CN110118979B (zh) 基于广义互熵的改进差分进化算法估计多径参数的方法
CN103793614B (zh) 一种突变滤波方法
CN110361967A (zh) 滑模观测器的构建方法
CN104407366B (zh) 一种对伪距进行平滑处理的方法
CN111409865B (zh) 基于交会概率的深空探测器接近段制导方法
CN106053936B (zh) 一种获取电学信号瞬时频率的方法及系统
CN106153046B (zh) 一种基于自适应Kalman滤波的陀螺随机噪声AR建模方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant