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

Abstract

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

Description

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

技术领域

本发明涉及航天器微推进领域的推力测量技术。

背景技术

推力是微推力器的重要推进性能指标，对工况参数选择、工况优选、推力器设计以及推力器应用都有很强的工程指导意义。目前，微推力器的推力一般采用扭摆二阶振动结构测量。如果推力的作用时间小于扭摆测量系统固有周期的0.25倍时，推力作用在扭摆系统上的效应相当于瞬间脉冲冲量作用的效应，只能根据测量系统的响应辨识所施加推力的冲量大小，冲量辨识技术已经较为成熟。如果推力的作用时间大于扭摆测量系统固有周期的0.25倍时，可以根据测量系统的响应辨识推力大小。通常，推力辨识方法是根据推力与扭摆系统稳态变化量之间的线性关系。该方法只能获取推力器试验条件下的单点推力值，无法还原推力随时间的变化情况，从而无法分析推力器开机阶段、平稳工作阶段和关机阶段的推进性能。同时，微推力器的推力测量也具有很大的挑战性，易受多种噪声影响，如测量环境的振动、测量环境的气流、位移传感器等。因此，如何在多噪声条件下辨识还原推力曲线是亟待解决的技术问题。

发明内容

本发明要解决的技术问题在于：针对多噪声影响下的微推力测量和辨识还原推力曲线的问题，提出一种平滑降噪优化还原方法。该方法综合运用正交多项式局部滑动拟合方法、推力反向计算的复化辛普森方法、系统响应预测值的正向计算方法三个核心方法，试图从计算方法、计算过程和计算结果上减小推力误差，获取推力估计值，并采用系统响应预测值与实际系统响应测量值比较方法，确定最佳的推力估计值和上下包络线。

具体技术方案为：

1、微推力平滑降噪优化还原方法的实施步骤

第一步，系统参数标定和校验：采用高精度阶跃标定力产生扭摆系统的阶跃响应；利用正交多项式局部滑动拟合方法确定系统响应的平均位置曲线；根据系统响应的平均位置曲线确定振动频率和阻尼比的初始估计值及置信区间；在振动频率和阻尼比置信区间内对系统参数进行微调和验校，确定最优振动频率、最优阻尼比和最后扭转刚度系数。

第二步，确定最小划分区间数目：采用复化辛普森积分方法，将推力和系统响应的积分方程离散化为线性方程组，反向计算推力；采用与待测推力变化范围和趋势相似的模拟推力，在测量噪声忽略不计的条件下，确定积分方程离散化的最小划分区间数目。

第三步，确定推力最优估计值：采用正交多项式局部滑动拟合方法，确定推力作用下实际系统响应的平均位置曲线和上下包络线；从最小划分区间数目开始，逐步增大划分区间数目以及局部滑动拟合函数，综合运用多次局部滑动拟合方法、积分方程离散化反向计算推力、样条函数插值获得推力连续函数和微分方程正向计算系统响应预测值等方法，产生一系列推力估计值与对应的一系列系统响应预测值；根据系统响应预测值是否在实际系统响应测量值的包络线以内，以及是否在实际系统响应测量值的平均位置曲线附近上下波动最小，确定推力最优估计值；对推力最优估计值，采用正交多项式局部滑动拟合方法，确定推力最优估计值的平均位置曲线和上下包络线。

2、系统参数标定和校验的步骤

第一步，给扭摆系统施加阶跃标定力f₀，获得实际系统响应测量值Θ(t_j)，j＝1,2,…,m，采用正交多项式局部滑动拟合方法对实际系统响应Θ(t_j)平滑降噪，平均位置曲线估计值

并判读出极值点对应时间t_Mi和扭转角θ(t_Mi)，i＝1,2,…,n。

第二步，确定振动频率的初始估计值

和标准差

式中，

第三步，确定阻尼比的初始估计值

和标准差

式中，

利用方程

获得，i＝2,3,…,n。

第四步，系统参数微调和验校：首先，在振动频率的置信区间

和阻尼比的置信区间

内，采样小步长搜索方法，不断改变振动频率和阻尼比；其次，在相应的振动频率和阻尼比下，计算出一系列稳态扭转角估计值

以及一系列平均位置曲线

最后，构造残差

根据准则“残差为零均值、零均值附近对称分布、取正值和负值的次数基本相等”，寻找到最优振动频率

最优阻尼比

最优稳态扭转角估计值

第五步，确定扭转刚度系数最优估计值

和误差

式中，f₀和df₀为标定力及其误差，L_f和dL_f为标定力力臂及其误差，

3、确定最小划分区间数目的步骤

第一步，获取待测推力作用下的实际系统响应测量值Θ_F(t_i)(i＝1,2,…,m)，采用正交多项式局部滑动拟合方法获得实际系统响应的平均位置曲线测量值θ_F(t_i)(i＝1,2,…,m)。

第二步，采用推力反向计算的复化辛普森方法计算推力值，通过改变划分区间数目M，分析待测推力变化范围和趋势。

第三步，根据被测推力变化范围和趋势，选择模拟推力f(τ)(连续函数表示)，在测量噪声忽略不计条件下，根据推力微分方程正向计算模拟推力作用下的系统响应θ_fi(i＝0,1,2,…,m)。

第四步，根据模拟推力下的系统响应θ_fi(i＝0,1,2,…,m)，采用推力反向计算的复化辛普森方法，计算推力值f_i'(i＝0,1,…,m-1)。

第五步，根据推力相对截断误差Δf_i/f_i＝(f_i'-f_i)/f_i(i＝0,1,…,m-1)，选择最小划分区间数目M_min，将推力截断误差控制在要求的水平。

4、确定推力最优估计值的步骤

第一步，根据待测推力作用下的实际系统响应测量值 [t_i,Θ_F(t_i)](i＝1,2,…,m)，采用正交多项式局部滑动拟合方法，确定实际系统响应的平均位置曲线[t_i,θ_F(t_i)](i＝1,2,…,m)，以及确定实际系统响应测量值的上包络线Θ_Fu(t_i)＝Θ_F(t_i)[Θ_F(t_i)＞θ_F(t_i)]和下包络线Θ_Fd(t_i)＝Θ_F(t_i)[Θ_F(t_i)＜θ_F(t_i)]。

第二步，从最小划分区间数目M_min开始，逐步增大划分区间数目，以及采用不同局部数据窗长度(如3、5、7等)的二阶正交多项式滑动拟合方法，在任一划分区间数目和局部滑动拟合函数下，首先，根据实际系统响应的平均位置曲线数据[t_i,θ_F(t_i)](i＝1,2,…,m)，采用推力反向计算的复化辛普森方法，计算推力值

其次，采用正交多项式局部滑动拟合方法，对计算得到的推力值平滑降噪滤波处理，确定推力估计值为

最后，采用对推力估计值样条函数插值方法，获得推力的连续函数，采用系统响应预测值的正向计算方法，计算系统响应预测值[t_i,θ'_F(t_i)](i＝1,2,…,m)；从而得到一系列推力估计值

以及对应的一系列系统响应预测值[θ'_F(t_i)]^j(i＝1,2,…,m；j＝M_min,M_min+1,…)。

第三步，采用系统响应预测值与实际系统响应测量值比较方法，选择推力最优估计值

选择原则一为：首先，真实推力应使系统响应预测值包含在实际系统响应测量值的上下包络线之内，即Θ_Fd(t_i)≤[θ'_F(t_i)]^j≤Θ_Fu(t_i)；选择原则二为：真实推力应使系统响应预测值在实际系统响应的平均位置曲线附近上下波动最小，即残差δ_i＝[θ'_F(t_i)]^j-θ_F(t_i)满足Θ_Fd(t_i)-θ_F(t_i)≤δ_i≤Θ_Fu(t_i)-θ_F(t_i)，且残差为零均值、在零均值附近波动最小及对称分布。

第四步，对推力最优估计值

采用正交多项式局部滑动拟合方法，确定推力最优估计值的平均位置曲线

以及上包络线

和下包络线

5、正交多项式局部滑动拟合方法

正交多项式局部滑动拟合方法的基本思想为：1)以每个数据点为中心对称选择若干个数据点，作为数据窗，根据数据窗内数据点，利用正交多项式，采用最小二乘方法进行拟合；2)通过数据窗逐点局部滑动过程，实现整个时间区间上数据点的拟合，寻找整个时间区间上最佳拟合曲线。

每一个数据窗下的正交多项式局部拟合方法的实施步骤为：

第一步，对于试验数据点(t_i,y_i)(i＝0,1,2,…,m)，以拟合点t_i为中心截取(t_i-p,y_i-p)，(t_i-p+1,y_i-p+1)，…，(t_i,y_i)，…，(t_i+p-1,y_i+p-1)，(t_i+p,y_i+p)等2p+1 个数据点(p＝1,2,…)，数据窗长度为L＝2p+1。

第二步，对数据窗内数据点构建正交多项式拟合曲线 P(t)＝a₀p₀(t)+a₁p₁(t)+…+a_np_n(t)(n为拟合曲线阶数)，正交多项式 p_k(t)(k＝0,1,…,n)的递推公式为

其中，

系数a_j(j＝0,1,…,n)根据线性方程组

求解。

6、推力反向计算的复化辛普森方法

推力反向计算的复化辛普森方法的基本思想为：采用复化辛普森公式，将待测推力F(t)与系统响应Θ(t)的积分方程

离散化为线性方程组，根据推力作用下的实际系统响应测量值为[t_i,Θ(t_i)](i＝0,1,2,…,m)递推计算出推力值[t_i,F(t_i)](i＝0,1,2,…,m)。

若实际系统响应测量值采样时间步长为h，设t_i＝ih(i＝0,1,2,…,m)；设 t₀＝0时，Θ(t₀)＝0；令Θ(t_i)＝Θ_i和F(t_i)＝F_i；采用复化辛普森公式将推力积分方程离散化为线性方程组过程如下：

当i＝1时，a₁₀F₀＝C_fΘ₁，

当i＝2时，a₂₀F₀+a₂₁F₁＝C_fΘ₂，

当i＝3时，a₃₀F₀+a₃₁F₁+a₃₂F₂＝C_fΘ₃，

当i＝2k(k＝2,3,...)时，

当i＝2k+1(k＝2,3,...)时，

7、系统响应预测值的微分方程正向计算方法

系统响应预测值的微分方程正向计算方法的基本思想为：根据推力估计值，采用样条函数插值方法获得推力的连续函数S(t)，再根据振动微分方程

计算系统响应预测值 [t_i,θ_p(t_i)](i＝1,2,…,m)，其中，

t₀＝0，θ_p(0)＝0。

已知推力估计值

设

将推力作用时间区间[0,T₀]划分为n′＝m-1个小区间(节点数目为n′+1)，有 0＝t₀＜t₁＜t₂＜…＜t_n′-1＜t_n′＝T₀，小区间长度为h_i-1＝t_i-t_i-1(1≤i≤n′)，在子区间[t_i-1,t_i](1≤i≤n′)上，三次样条函数为 S(t)＝F_i-1+(F_i-F_i-1)(t-t_i-1)/(t_i-t_i-1)+[a_i(t-t_i-1)+b_i(t-t_i)](t-t_i-1)(t-t_i)，a_i和b_i为待定系数。

a_i和b_i的计算表达式为：

式中，M_i可由

计算获得，其中，

λ_i＝1-μ_i，

与现有技术相比，采用本发明可以达到以下技术效果：

1.本发明能够还原出待测推力的最佳推力曲线以及推力曲线的上下包络线，这是区别已有技术的最大特色；

2.本发明具有较好的通用性，适用于类似扭摆系统的二阶振动测量系统的推力曲线还原。

附图说明

图1为本发明的技术路线图。

具体实施方案

结合附图对用于扭摆系统的微推力平滑降噪优化还原方法做进一步详细描述。微推力平滑降噪优化方法的基本思路是：①采用实际系统响应测量值平滑降噪方法，从整个计算过程上减小测量噪声影响；采用高精度数值积分离散化方法，从计算方法上减小推力误差；采用推力计算值平滑降噪方法，从计算结果上减小推力噪声误差，最终获得推力估计值。②采用对推力估计值样条函数插值方法，获得推力的连续函数；采用振动微分方程正向计算方法，获得系统响应预测值；采用系统响应预测值与实际系统响应测量值比较方法，验证推力估计值的优化程度，在一系列推力估计值中确定最佳的推力估计值。

微推力平滑降噪优化还原方法的步骤概述为：

第一步，系统参数标定和校验：采用高精度阶跃标定力产生扭摆系统的阶跃响应；利用正交多项式局部滑动拟合方法确定系统响应的平均位置曲线；根据系统响应的平均位置曲线确定振动频率和阻尼比的初始估计值及置信区间；在振动频率和阻尼比置信区间内对系统参数进行微调和验校，确定最优振动频率、最优阻尼比和最后扭转刚度系数；

第二步，确定最小划分区间数目：采用复化辛普森积分方法，将推力和系统响应的积分方程离散化为线性方程组，反向计算推力；采用与待测推力变化范围和趋势相似的模拟推力，在测量噪声忽略不计的条件下，确定积分方程离散化的最小划分区间数目；

第三步，确定推力最优估计值：采用正交多项式局部滑动拟合方法，确定推力作用下实际系统响应的平均位置曲线和上下包络线；从最小划分区间数目开始，逐步增大划分区间数目以及局部滑动拟合函数，综合运用多次局部滑动拟合方法、积分方程离散化反向计算推力、样条函数插值获得推力连续函数和微分方程正向计算系统响应预测值等方法，产生一系列推力估计值与对应的一系列系统响应预测值；根据系统响应预测值是否在实际系统响应测量值的包络线以内，以及是否在实际系统响应测量值的平均位置曲线附近上下波动最小，确定推力最优估计值；对推力最优估计值，采用正交多项式局部滑动拟合方法，确定推力最优估计值的平均位置曲线和上下包络线。

图1为本发明的技术路线图，结合附图具体阐述微推力平滑降噪优化还原方法的推导原理和实施步骤。

1、系统参数标定和校验

(1)采样阶跃响应方法，利用高精度标定力，标定系统参数

在已知阶跃标定力f₀作用下，扭摆测量系统的实际系统响应为

式中，稳态扭转角θ(∞)是未知量；振动频率ω_d、阻尼比ζ和扭转刚度系数 k等是未知量；Δθ(t)～N(0,σ²)为标定过程中环境、位移传感器和标定力等产生的测量噪声，测量噪声方差σ²是未知量。

(2)标定振动频率和阻尼比

由于测量噪声Δθ(t)～N(0,σ²)的影响，实际系统响应Θ(t)在平均位置曲线θ(t)附近上下起伏波动，给极值点对应时间和扭转角判读带来困难，因此，采样正交多项式局部滑动拟合方法，对实际系统响应Θ(t)平滑降噪滤波处理后，再判读极值点对应时间和t_Mi扭转角θ(t_Mi)。

振动频率为

振动频率的估计值和标准差为

极值点对应扭转角为

可得

令

式中，0≤x_ζ≤1。得到x_ζ的方程

可解得阻尼比为

阻尼比的估计值和标准差为

(3)标定转刚度系数、稳态扭转角、测量噪声方差

在已知振动频率估计值

阻尼比估计值

条件下，设系统响应的平均位置曲线为

式中，稳态扭转角θ(∞)是待估计量。

根据实际系统响应测量值Θ(t_i)(i＝1,2,…,m)和平均位置曲线估计值

构造残差为

采用最小二乘方法，使得残差的平方和最小，确定稳态扭转角θ(∞)的估计值

残差为

令其平方和为最小，即

可得稳态扭转角θ(∞)的估计值

为

其方差为

通过稳态扭转角θ(∞)的估计值

计算扭转刚度系数的估计值，为

其绝对误差和相对误差，分别为

式中，df₀为标定力的误差，dL_f为力臂的误差。

测量噪声Δθ(t)～N(0,σ²)方差σ²的估计值，为

从而得到平均位置曲线为

确定平均位置曲线的过程为：首先确定振动频率和阻尼比的估计值，其次采用最小二乘方法确定稳态扭转角。

(4)系统参数的验校和微调

振动频率估计值

阻尼比估计值

和稳态扭转角估计值

等，是否为最佳估计值需要检验和校正。验校的准则为：最佳的振动频率估计值

阻尼比估计值

和稳态扭转角估计值

等，应使得残差

为零均值、零均值附近对称分布、取正值和负值的次数基本相等。

首先，在振动频率的置信区间

和阻尼比的置信区间

内，采样小步长搜索方法，不断改变振动频率和阻尼比。

其次，采用步骤③的方法，计算稳态扭转角估计值

并确定平均位置曲线

最后，构造残差

根据准则“残差为零均值、零均值附近对称分布、取正值和负值的次数基本相等”，寻找到最佳振动频率、阻尼比、稳态扭转角、扭转刚度系数。

2、确定最小划分区间数目

(1)根据被测推力特点，选择模拟推力

首先，根据实际系统响应测量值Θ(t_i)(i＝1,2,…,m)，采用正交多项式局部滑动拟合方法，计算实际系统响应的平均位置曲线测量值θ(t_i)(i＝1,2,…,m)。

其次，采用推力反向计算的复化辛普森方法，计算推力值，通过改变划分区间数目m，考察被测推力变化范围和趋势。

最后，根据被测推力变化范围和趋势，采用相似的模拟推力表示被测推力变化范围和趋势，模拟推力采用连续函数表示。

(2)根据模拟推力，选择最小划分区间数目

首先，采用模拟推力f(τ)(连续函数表示)，根据推力微分方程正向计算系统响应θ_i(i＝0,1,2,…,m)；

其次，在测量噪声忽略不计条件下，根据模拟推力下的系统响应θ_i(i＝0,1,2,…,m)，采用推力反向计算的复化辛普森方法，计算推力值 F_i(i＝0,1,…,m-1)。

最后，根据推力相对截断误差Δf_i/f_i＝(F_i-f_i)/f_i(i＝0,1,…,m-1)，选择最小划分区间数目m_min，将推力截断误差控制在要求的水平。

在推力作用时间区间[0,T₀]内，合理选择最小划分区间数目m_min，将推力截断误差控制在要求范围内，同时由于划分区间数目最小，达到降低推力噪声误差的目的。例如，将推力的截断误差控制在0.1％～0.5％范围内，选择最小划分区间数目。

3、推力反算和响应预测

(1)根据实际系统响应测量值[t_i,Θ(t_i)]，采用正交多项式局部滑动拟合方法，确定系统响应的平均位置曲线

达到从整个计算过程上减小测量噪声影响的目的。

(2)根据实际系统响应的平均位置曲线数据

采用推力反向计算的复化辛普森方法，计算推力值

式中，

达到采用高精度3点和4点交替的复化辛普森离散化方法，减小划分区间数目，降低推力误差的目的。

(3)采用正交多项式局部滑动拟合方法，对计算得到的推力值平滑降噪滤波处理，确定推力估计值为

从而达到从计算结果上再次减小推力噪声误差的目的，最终获得推力估计值。

(4)采用对推力估计值样条函数插值方法，获得推力的连续函数，采用系统响应预测值的正向计算方法，计算系统响应预测值。

首先，已知推力的估计值

设

将推力作用时间区间[0,T₀]划分为n′＝m-1个小区间(节点数目为n′+1)，有

0＝t₀＜t₁＜t₂＜…＜t_n′-1＜t_n′＝T₀

小区间长度为h_i-1＝t_i-t_i-1(1≤i≤n′)，在子区间[t_i-1,t_i](1≤i≤n′)上，推力估计值的三次样条函数为

S(t)＝F_i-1+F(t_i-1,t_i)(t-t_i-1) +[a_i(t-t_i-1)+b_i(t-t_i)](t-t_i-1)(t-t_i)

其次，利用微分方程数值计算可方便地获得高精度数值解特点，采用振动微分方程正向计算方法，计算系统响应预测值。扭摆振动微分方程表示为

式中，

S(t)为推力估计值的样条插值连续函数。

已知推力估计值

采用系统响应预测值的正向计算方法，计算的系统响应预测值为[t_i,θ_p(t_i)](i＝1,2,…,m)。

因此，针对振动微分方程正向计算高精度特点、但要求推力是连续函数特点，实现了由推力估计值计算系统响应预测值。

4、确定推力最优估计值

(1)采用正交多项式局部滑动拟合方法，确定实际系统响应的平均位置曲线、测量噪声方差、实际系统响应的上下包络线

首先，根据实际系统响应测量值[t_i,Θ(t_i)]，采用正交多项式局部滑动拟合方法，确定系统响应的平均位置曲线

其次，通过残差

分析，确定测量噪声的方差。测量噪声Δθ(t)～N(0,σ²)方差σ²的估计值，为

最后，确定实际系统响应测量值的上下包络线，上包络线和下包络线，分别为

(2)从最小划分区间数目开始，逐步增大划分区间数目，以及选择不同平滑降噪滤波函数，综合运用上述推力反算和响应预测的方法，产生一系列推力估计值和对应的一系列系统响应预测值。

首先，选择正交多项式局部滑动拟合方法，选择局部拟合能力突出的二阶正交多项式(抛物线正交多项式)，选择局部数据窗长度为 2p+1＝3,5,7,…。

其次，根据实际系统响应测量值[t_i,Θ(t_i)]，从最小划分区间数目m_min开始，逐步增大划分区间数目，综合运用上述推力反算和响应预测的方法，得到一系列推力估计值

和对应的一系列系统响应预测值[θ_p(t_i)]^j，为

[θ_p(t_i)]^j(i＝1,2,…,m；j＝m_min,m_min+1,…)

(3)采用系统响应预测值与实际系统响应测量值比较方法，确定最佳的推力估计值。

首先，针对真实推力应使系统响应预测值包含在实际系统响应测量值的上下包络线之内的特点，使得Θ_d(t_i)≤[θ_p(t_i)]^j≤Θ_u(t_i)。

其次，针对真实推力应使系统响应预测值在实际系统响应的平均位置曲线附近上下波动最小的特点，考察残差

随着时间变化，当残差满足

且残差为零均值、在零均值附近波动最小及对称分布时，寻找到最佳的推力估计值。

最后，对最佳的推力估计值，采用正交多项式局部滑动拟合方法，平滑降噪滤波处理，确定推力估计值和推力误差。

因此，通过系统响应预测值最佳逼近实际系统响应的平均位置曲线，达到了寻找最佳推力估计值的目的、以及推力测量和误差分析的目的。

从上述平滑降噪优化还原方法的具体应用过程可以看出，涉及正交多项式局部滑动拟合方法、推力反向计算的复化辛普森方法和系统响应预测值的正向计算方法三个核心方法，下面分别阐述。

1、正交多项式局部滑动拟合方法

(1)实际系统响应测量值与平均位置曲线

实际系统响应测量值总是包含测量噪声，实际系统响应测量值为

Θ(t)＝θ(t)+Δθ(t)

式中，θ(t)为理想系统响应，Δθ(t)为测量噪声。

根据实际系统响应测量值[t_i,Θ(t_i)](i＝0,1,2,…,m)，采用正交多项式局部滑动拟合方法，经过平滑降噪滤波处理，确定平均位置曲线

作为理想系统响应θ(t)的估计值。

(2)正交多项式局部滑动拟合方法

①正交多项式拟合方法

实际系统响应测量值为(t_i,Θ_i)(i＝0,1,2,…,m)(Θ_i＝Θ(t_i))，由已知线性无关函数

构造多项式P(t)进行拟合，为

根据最小二乘法，可得线性方程组

通过该线性方程组可解得系数a_j(j＝0,1,…,n)。

对于已知线性无关函数

上述方程组存在唯一解，但是该线性方程组往往具有病态特性，导致系数的计算误差较大，因此，选取正交多项式进行拟合。

②正交多项式局部滑动拟合方法

根据实际系统响应测量值(t_i,Θ_i)(i＝0,1,2,…,m)，利用正交多项式P(t) 拟合实际系统响应曲线，避免了求解病态线性方程组，提高了系数 a_j(j＝0,1,…,n)的求解精度。但是，由于实际系统响应是多极值点、多拐点、曲率变化大的复杂曲线，在整个时间t₀≤t≤t_m上全程拟合，不能反映每个局部数据取值特点，出现拟合曲线振荡和拟合误差大现象。因此，提出正交多项式局部滑动拟合方法。

正交多项式局部滑动拟合方法是：1)以每个数据点为中心对称选择若干个数据点，作为数据窗，根据数据窗内数据点，利用正交多项式，采用最小二乘方法进行拟合；2)通过数据窗逐点局部滑动过程，实现整个时间区间上数据点的拟合，寻找最佳拟合曲线，避免了全程数据点拟合的误差。

对于试验数据点(t_i,Θ_i)(i＝0,1,2,…,m)中，以拟合点t_i为中心，取 (t_i-p,Θ_i-p)，(t_i-p+1,Θ_i-p+1)，…，(t_i,Θ_i)，…，(t_i+p-1,Θ_i+p-1)，(t_i+p,Θ_i+p) 等2p+1点(p＝1,2,…)，数据窗长度为l+1＝2p+1，进行正交多项式局部滑动拟合。

具体构造正交多项式p_k(t)(k＝0,1,…,n)的递推公式为

其中

最后，采用最小二乘方法，利用正交多项式拟合曲线为

P(t)＝a₀p₀(t)+a₁p₁(t)+…+a_np_n(t)

其中，n为拟合曲线阶数。

2、推力反向计算的复化辛普森方法

(1)推力微分方程的正向计算和推力积分方程的反向计算

对于扭摆测量系统，推力与系统响应之间关系，采用微分方程表示为

式中，ζ为阻尼比，ω_n为固有频率，J为转动惯量。

上述推力与系统响应的微分方程，也可表示为推力与系统响应的积分方程形式，为

在已知推力f(t)条件下，由推力微分方程，可方便地高精度计算系统响应，这是推力微分方程的正向计算问题；在已知系统响应θ(t)条件下，由推力积分方程，也可计算推力，这是推力积分方程的反向计算问题。

(2)推力积分方程离散化为线性方程组

采用复化辛普森公式，将以下推力积分方程离散化为线性方程组

式中，C_f＝Jω_d/L_f。

设已知实际系统响应测量值为[t_i,Θ(t_i)](i＝0,1,2,…,m)，采样时间步长为h，t_i＝ih，当t₀＝0时Θ(t₀)＝0，令Θ(t_i)＝Θ_i和F(t_i)＝F_i。当t_i＝ih时，可知

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

①当i＝1时，根据梯形公式，可知

可得

②当i＝2时，根据3点辛普森公式，可知

可得

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

③当i＝3时，根据4点辛普森公式，可知

可得

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

④当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点辛普森公式。可得

即推力积分方程离散化线性方程组，为

式中，a_ij＝0(j≥i)。

因此，根据实际系统响应测量值为(t_i,Θ_i)(i＝0,1,2,…,m)，可反向计算推力值F_i(i＝0,1,2,…,m-1)。这种根据实际系统响应测量值，反向计算推力的方法，称为推力反向计算的复化辛普森方法。推力反向计算的复化辛普森方法，截断误差与步长h⁴成正比，可高精度计算推力值。

(3)推力误差包括推力截断误差和推力噪声误差

实际系统响应测量值总是包含测量噪声，实际系统响应测量值为

Θ(t)＝θ(t)+Δθ(t)

式中，θ(t)为理想系统响应，Δθ(t)为测量噪声。

在没有测量噪声Δθ(t)的理想情况下，理想系统响应为θ(t)，推力积分方程为

式中，C_f＝Jω_d/L_f，f(τ)为推力真值。

通过积分方程离散化为线性方程组，将积分表示为有限和形式，为

积分计算的截断误差，造成推力计算误差Δf(τ_i)＝F(τ_i)-f(τ_i)，这种推力积分方程离散化时截断误差引起的推力误差，称为推力截断误差。由于推力积分方程的被积分是振荡函数，划分区间数目足够多时，才能控制推力截断误差。

实际系统响应Θ(t)＝θ(t)+Δθ(t)包含测量噪声，推力积分方程离散化为线性方程组，为

测量噪声Δθ(t)造成推力计算误差Δf(τ_i)＝F(τ_i)-f(τ_i)，这种测量噪声引起的推力误差，称为推力噪声误差。由于离散化线性方程组的严重病态特性和系统响应测量值包含测量噪声，划分区间数目多时，造成推力噪声误差急剧增大。

推力反向计算的难点在于推力截断误差和推力噪声误差相互制约限制。划分区间数目少时，造成推力截断误差增大；划分区间数目多时，造成推力噪声误差急剧增大。

推力反向计算的复化辛普森方法，针对积分方程离散化的划分区间数目造成推力截断误差和推力噪声误差相互制约限制的难点，采用高精度的复化辛普森(3点和4点交替运用)数值积分方法，将推力积分方程离散化为线性方程组，可在划分区间数目较少前提下，实现尽量减小推力截断误差和推力噪声误差，达到减小推力误差的目的。

3、系统响应预测值的正向计算方法

(1)根据推力估计值，采用样条函数插值方法，获得推力的连续函数

针对振动微分方程正向计算高精度特点、但要求推力是连续函数特点，由推力估计值，采用样条函数插值方法，获得推力的连续函数。

利用样条函数具有函数光滑和导数连续的特点，采用样条函数插值方法，获得推力的连续函数。已知推力估计值

设

将推力作用时间区间[0,T₀]划分为n′＝m-1个小区间(节点数目为n′+1)，有

0＝t₀＜t₁＜t₂＜…＜t_n′-1＜t_n′＝T₀

小区间长度为h_i-1＝t_i-t_i-1(1≤i≤n′)，在子区间[t_i-1,t_i](1≤i≤n′)上，三次样条函数为

S(t)＝F_i-1+F(t_i-1,t_i)(t-t_i-1) +[a_i(t-t_i-1)+b_i(t-t_i)](t-t_i-1)(t-t_i)

式中，F(t_i-1,t_i)＝(F_i-F_i-1)/(t_i-t_i-1)为一阶均差，a_i和b_i为待定系数。

当i＝1,2,…,n′-1时，可得方程

μ_iM_i-1+2M_i+λ_iM_i+1＝d_i

其中

显然，插值区间等份时，μ_i＝λ_i＝1/2。

需要补充2个边界条件方程，推力的左边和右边边界看作第一类边界问题，已知F′₀和F′_n。补充得到两个方程

从而得到三对角方程

解得M_i(i＝0,1,…,n′)，可求得待定系数a_i(i＝1,2,…,n′)和b_i(i＝1,2,…,n′)，为

式中，小区间长度为h_i-1＝t_i-t_i-1。其中F′₀和F′_n，可采用数值微分公式计算，例如

(2)根据推力估计值的连续函数，采用微分方程正向计算方法，获得系统响应预测值

利用微分方程正向计算可方便地获得高精度数值解特点，采用振动微分方程正向计算方法，计算系统响应预测值。

扭摆振动微分方程表示为

式中，

S(t)为推力样条插值连续函数，为

S(t)＝F_i-1+F(t_i-1,t_i)(t-t_i-1) +[a_i(t-t_i-1)+b_i(t-t_i)](t-t_i-1)(t-t_i)

已知推力估计值

采用振动微分方程正向计算方法，计算系统响应预测值[t_i,θ_p(t_i)](i＝1,2,…,m)，其中，t₀＝0和θ_p(0)＝0。

在振动微分方程正向计算中，推力估计值的插值点数目过少(或积分方程离散化的划分区间数目过少)，样条插值函数不能准确表示推力随着时间变化，造成系统响应预测值的误差，因此，离散化的划分区间数目不但影响推力截断误差和推力噪声误差，还影响系统响应预测值的误差。

本申请不局限于说明书和权利要求文字部分所限定的内容，任何本领域范围内公知的修改和变化都属于本申请的范围，说明书具体实施例部分仅是对本发明示例性的说明，不是对本发明的具体限定。

Claims (7)

1.用于扭摆系统的微推力平滑降噪优化还原方法，其特征在于：

用于扭摆系统的微推力平滑降噪优化还原方法的实施步骤为：

第一步，系统参数标定和校验：采用高精度阶跃标定力产生扭摆系统的阶跃响应；利用正交多项式局部滑动拟合方法确定系统响应的平均位置曲线；根据系统响应的平均位置曲线确定振动频率和阻尼比的初始估计值及置信区间；在振动频率和阻尼比置信区间内对系统参数进行微调和验校，确定最优振动频率、最优阻尼比和最优扭转刚度系数；

第二步，确定最小划分区间数目：采用复化辛普森积分方法，将推力和系统响应的积分方程离散化为线性方程组，反向计算推力；采用与待测推力变化范围和趋势相似的模拟推力，在测量噪声忽略不计的条件下，确定积分方程离散化的最小划分区间数目；

第三步，确定推力最优估计值：采用正交多项式局部滑动拟合方法，确定推力作用下实际系统响应的平均位置曲线和上下包络线；从最小划分区间数目开始，逐步增大划分区间数目以及局部滑动拟合函数，综合运用多次局部滑动拟合方法、积分方程离散化反向计算推力、样条函数插值获得推力连续函数和微分方程正向计算系统响应预测值方法，产生一系列推力估计值与对应的一系列系统响应预测值；根据系统响应预测值是否在实际系统响应测量值的包络线以内，以及是否在实际系统响应测量值的平均位置曲线附近上下波动最小，确定推力最优估计值；对推力最优估计值，采用正交多项式局部滑动拟合方法，确定推力最优估计值的平均位置曲线和上下包络线。

2.如权利要求1所述的用于扭摆系统的微推力平滑降噪优化还原方法，其特征在于：

系统参数标定和校验的步骤为：

第一步，给扭摆系统施加阶跃标定力f₀，获得实际系统响应测量值Θ(t_j)，j＝1,2,…,m，采用正交多项式局部滑动拟合方法对实际系统响应Θ(t_j)平滑降噪，平均位置曲线估计值

并判读出极值点对应时间t_Mi和扭转角θ(t_Mi)，其中M表示极值点，i＝1,2,…,n；

第二步，确定振动频率的初始估计值

和标准差

式中，

i＝1,2,…,n；

第三步，确定阻尼比的初始估计值

和标准差

式中，

表示求解阻尼比的衰减量，利用方程

获得，i＝2,3,…,n得到；

第四步，系统参数微调和验校：首先，在振动频率的置信区间

和阻尼比的置信区间

内，采样小步长搜索方法，不断改变振动频率和阻尼比；其次，在相应的振动频率和阻尼比下，计算出一系列稳态扭转角估计值

以及一系列平均位置曲线

最后，构造残差

根据准则“残差为零均值、零均值附近对称分布、取正值和负值的次数相等”，寻找到最优振动频率

最优阻尼比

最优稳态扭转角估计值

第五步，确定扭转刚度系数最优估计值

和误差

式中，f₀和df₀为标定力及其误差，L_f和dL_f为标定力力臂及其误差，

3.如权利要求2所述的用于扭摆系统的微推力平滑降噪优化还原方法，其特征在于：

确定最小划分区间数目的步骤为：

第一步，获取待测推力作用下的实际系统响应测量值

Θ_F(t_i)(i＝1，2，…，m)，

采用正交多项式局部滑动拟合方法获得实际系统响应的平均位置曲线测量值θ_F(t_i)(i＝1，2，…，m)；

第二步，采用推力反向计算的复化辛普森方法计算推力值，通过改变划分区间数目M，分析待测推力变化范围和趋势；

第三步，根据被测推力变化范围和趋势，选择连续函数表示的模拟推力f(τ)，在测量噪声忽略不计条件下，根据推力微分方程正向计算模拟推力作用下的系统响应θ_fi(i＝0，1，2，…，m)；

第四步，根据模拟推力下的系统响应θ_fi(i＝0，1，2，…，m)，采用推力反向计算的复化辛普森方法，计算推力值f_i′(i＝0，1，…，m-1)；

第五步，根据推力相对截断误差Δf_i/f_i＝(f_i′-f_i)/f_i(i＝0，1，…，m-1)，选择最小划分区间数目M_min，将推力截断误差最小化。

4.如权利要求3所述的用于扭摆系统的微推力平滑降噪优化还原方法，其特征在于：

确定推力最优估计值的步骤为：

第一步，根据待测推力作用下的实际系统响应测量值

[t_i，Θ_F(t_i)](i＝1，2，…，m)，采用正交多项式局部滑动拟合方法，确定实际系统响应的平均位置曲线[t_i，θ_F(t_i)](i＝1，2，…，m)，以及确定实际系统响应测量值的上包络线Θ_Fu(t_i)＝Θ_F(t_i)[Θ_F(t_i)＞θ_F(t_i)]和下包络线Θ_Fd(t_i)＝Θ_F(t_i)[Θ_F(t_i)＜θ_F(t_i)]；

第二步，从最小划分区间数目M_min开始，逐步增大划分区间数目，以及采用不同局部数据窗长度的二阶正交多项式滑动拟合方法，在任一划分区间数目和局部滑动拟合函数下，首先，根据实际系统响应的平均位置曲线数据[t_i，θ_F(t_i)](i＝1，2，…，m)，采用推力反向计算的复化辛普森方法，计算推力值

为利用复化辛普森方法离散化推力方程获得的下三角系数矩阵的逆，

推力导致的扭摆转角变化；其次，采用正交多项式局部滑动拟合方法，对计算得到的推力值平滑降噪滤波处理，确定推力估计值为

最后，采用对推力估计值样条函数插值方法，获得推力的连续函数，采用系统响应预测值的正向计算方法，计算系统响应预测值[t_i，θ′_F(t_i)](i＝1，2，…，m)；从而得到一系列推力估计值

以及对应的一系列系统响应预测值[θ′_F(t_i)]^j(i＝1，2，…，m；j＝M_min，M_min+1，…)；

第三步，采用系统响应预测值与实际系统响应测量值比较方法，选择推力最优估计值

选择原则一为：首先，真实推力应使系统响应预测值包含在实际系统响应测量值的上下包络线之内，即

Θ_Fd(t_i)≤[θ′_F(t_i)]^j≤Θ_Fu(t_i)；

选择原则二为：真实推力应使系统响应预测值在实际系统响应的平均位置曲线附近上下波动最小，即残差δ_i＝[θ′_F(t_i)]^j-θ_F(t_i)满足

Θ_Fd(t_i)-θ_F(t_i)≤δ_i≤Θ_Fu(t_i)-θ_F(t_i)，且残差为零均值、在零均值附近波动最小及对称分布；

第四步，对推力最优估计值

采用正交多项式局部滑动拟合方法，确定推力最优估计值的平均位置曲线

以及上包络线

和下包络线

5.如权利要求4所述的用于扭摆系统的微推力平滑降噪优化还原方法，其特征在于：

正交多项式局部滑动拟合方法的基本思想为：1)以每个数据点为中心对称点选择若干个数据点，作为数据窗，根据数据窗内数据点，利用正交多项式，采用最小二乘法进行拟合；2)通过数据窗逐点局部滑动过程，实现整个时间区间上数据点的拟合，寻找整个时间区间上最佳拟合曲线；

每一个数据窗下的正交多项式局部拟合方法的实施步骤为：

第一步，对于试验数据点(t_i，y_i)(i＝0，1，2，…，m)，以拟合点t_i为中心截取(t_i-p，y_i-p)，(t_i-p+1，y_i-p+1)，...，(t_i，y_i)，...，(t_i+p-1，y_i+p-1)，(t_i+p，y_i+p)，2p+1个数据点(p＝1，2，…)，数据窗长度为L＝2p+1；

第二步，对数据窗内数据点构建正交多项式拟合曲线

P(t)＝a₀p₀(t)+a₁p₁(t)+…+a_np_n(t)，其中n为拟合曲线阶数，正交多项式p_k(t)(k＝0，1，…，n)的递推公式为

其中，

系数a_j(j＝0，1，…，n)根据线性方程组

求解。

6.如权利要求5所述的用于扭摆系统的微推力平滑降噪优化还原方法，其特征在于：

推力反向计算的复化辛普森方法的基本思想为：采用复化辛普森公式，将待测推力F(t)与系统响应Θ(t)的积分方程

离散化为线性方程组，其中，F(τ)为对应于时间微元τ的待测微元推力，根据推力作用下的实际系统响应测量值为[t_i，Θ(t_i)](i＝0，1，2，…，m)递推计算出推力值

若实际系统响应测量值采样时间步长为h，设t_i＝ih(i＝0，1，2，…，m)；设t₀＝0时，Θ(t₀)＝0；令Θ(t_i)＝Θ_i和F(t_i)＝F_i；采用复化辛普森公式将推力积分方程离散化为线性方程组过程如下：

当i＝1时，a₁₀F₀＝C_fΘ₁，

当i＝2时，a₂₀F₀+a₂₁F₁＝C_fΘ₂，

当i＝3时，a₃₀F₀+a₃₁F₁+a₃₂F₂＝C_fΘ₃，

当i＝2k(k＝2，3，...)时，

其中，l为与k相关的序数；

当i＝2k+1(k＝2，3，...)时，

7.如权利要求6所述的用于扭摆系统的微推力平滑降噪优化还原方法，其特征在于：

系统响应预测值的微分方程正向计算方法的基本思想为：根据推力估计值，采用样条函数插值方法获得推力的连续函数S(t)，再根据振动微分方程

计算系统响应预测值

[t_i，θ_p(t_i)](i＝1，2，…，m)，

其中，

t₀＝0，θ_p(0)＝0；θ表示扭转角，

角加速度，

表示角速度；

已知推力估计值

设

将推力作用时间区间[0，T₀]划分为n′＝m-1个小区间，其中节点数目为n′+1，有0＝t₀＜t₁＜t₂＜…＜t_n′-1＜t_n′＝T₀，小区间长度为h_i-1＝t_i-t_i-1(1≤i≤n′)，在子区间[t_i-1，t_i](1≤i≤n′)上，三次样条函数为S(t)＝F_i-1+(F_i-F_i-1)(t-t_i-1)/(t_i-t_i-1)+[a_i(t-t_i-1)+b_i(t-t_i)](t-t_i-1)(t-t_i)，a_i和b_i为待定系数；

a_i和b_i的计算表达式为：

式中，

其中，

中间参数

λ_i＝1-μ_i，

Images