CN111609878B - 三自由度直升机系统传感器运行状态监测方法 - Google Patents

三自由度直升机系统传感器运行状态监测方法 Download PDF

Info

Publication number
CN111609878B
CN111609878B CN202010524010.0A CN202010524010A CN111609878B CN 111609878 B CN111609878 B CN 111609878B CN 202010524010 A CN202010524010 A CN 202010524010A CN 111609878 B CN111609878 B CN 111609878B
Authority
CN
China
Prior art keywords
time
mode
value
sensor
state
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.)
Active
Application number
CN202010524010.0A
Other languages
English (en)
Other versions
CN111609878A (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.)
Jiangnan University
Original Assignee
Jiangnan 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 Jiangnan University filed Critical Jiangnan University
Priority to CN202010524010.0A priority Critical patent/CN111609878B/zh
Publication of CN111609878A publication Critical patent/CN111609878A/zh
Application granted granted Critical
Publication of CN111609878B publication Critical patent/CN111609878B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01D—MEASURING NOT SPECIALLY ADAPTED FOR A SPECIFIC VARIABLE; ARRANGEMENTS FOR MEASURING TWO OR MORE VARIABLES NOT COVERED IN A SINGLE OTHER SUBCLASS; TARIFF METERING APPARATUS; MEASURING OR TESTING NOT OTHERWISE PROVIDED FOR
    • G01D18/00—Testing or calibrating apparatus or arrangements provided for in groups G01D1/00 - G01D15/00

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Testing Or Calibration Of Command Recording Devices (AREA)

Abstract

本发明公开了一种三自由度直升机系统传感器运行状态监测方法，包括以下步骤，计算出传感器的测量噪声协方差矩阵的估计值；分析所述传感器的测量噪声协方差矩阵估计值，判断出所述传感器的运行状态是否正常：若所述传感器的测量噪声协方差矩阵估计值偏离所述传感器的标称值范围，判断出所述传感器的运行状态不正常；若所述传感器的测量噪声协方差矩阵估计值包含在所述传感器的标称值范围内，判断出所述传感器的运行状态正常。本发明的三自由度直升机系统传感器运行状态监测方法，通过对测量噪声协方差矩阵的准确估计，来实时监测传感器工作状态，以保证系统安全稳定运行。

Description

三自由度直升机系统传感器运行状态监测方法
技术领域
本发明涉及一种传感器运行状态监测方法，具体涉及一种三自由度直升机系统传感器运行状态监测方法。
背景技术
在充满不确定性的工业环境中，由于设备的老化以及外部随机干扰，传感器测量数据通常会带来一定的误差。为了解决这一问题，通常的做法是在系统测量方程中引入一个扰动变量ν，该扰动变量表示传感器所有可能遇到的测量噪声。而且，通常将扰动变量ν假定为服从均值为零方差为R的高斯分布，其中，方差R可以描述测量值的不确定性。此外，测量噪声ν的强度可能会由于许多外部和内部因素而随时间变化，特别是在恶劣的工业条件下，这就使得测量噪声协方差R随时间而变化。因此，为了解决噪声强度随时间变化的而带来的问题，应监测传感器是否处于正常运行状态，并在传感器受到噪声变化的影响时描绘出测量数据的不确定性，从而为系统决策制定提供可靠信息。
现有的关于传感器监测的方法都需要预先了解测量噪声协方差，但在大多数情况下，系统测量噪声协方差都是时变未知的，即使通过选择适当的测量噪声协方差来处理某些异常情况，仍然有一些潜在异常情况的重要信息会丢失。因此，实时了解测量噪声协方差对于传感器工作状态的准确监测具有重要价值。
三自由度直升机系统作为一种随机跳变系统，目前没有良好的方法对系统传感器工作状态进行实时准确监测。
发明内容
本发明要解决的技术问题是提供一种三自由度直升机系统传感器运行状态监测方法，适用于三自由度直升机系统中传感器运行状态监测，通过对测量噪声协方差的准确估计，来实时监测传感器工作状态，以保证系统安全稳定运行。
为了解决上述技术问题，本发明提供了一种三自由度直升机系统传感器运行状态监测方法，包括以下步骤，
计算出传感器的测量噪声协方差矩阵的估计值；
分析所述传感器的测量噪声协方差矩阵估计值，判断出所述传感器的运行状态是否正常：若所述传感器的测量噪声协方差矩阵估计值偏离所述传感器的标称值范围，判断出所述传感器的运行状态不正常；若所述传感器的测量噪声协方差矩阵估计值包含在所述传感器的标称值范围内，判断出所述传感器的运行状态正常。
本发明一个较佳实施例中，进一步包括，计算传感器的测量噪声协方差矩阵的估计值包括以下步骤，
(1)创建所述传感器所在的三自由度直升机系统的随机跳变系统模型，
xk＝fk(xk-1,rk)+ωk；
yk＝gk(xk,rk)+νk；
k为时间索引；xk为k时刻系统状态值；yk为k时刻系统测量值；fk(·)为系统状态方程；gk(·)为系统测量方程；rk为k时刻系统模态；ωk为过程噪声且ωk服从均值为零的高斯分布即ωk～N(0,Qk)，Qk为过程噪声协方差矩阵；νk为测量噪声且νk服从均值为零的高斯分布即νk～N(0,Rk)，Rk为k时刻系统测量噪声协方差矩阵；xk-1为k-1时刻系统状态值；
(2)基于变分贝叶斯理论将每一模态下系统状态和测量噪声协方差的联合后验概率密度函数
Figure GDA0003047311250000031
用两个独立的概率密度函数
Figure GDA0003047311250000032
来表示，即
Figure GDA0003047311250000033
Figure GDA0003047311250000034
表示系统k时刻的模态为s，y1:k为从时刻1到时刻k的测量值序列；
用一组加权粒子描述系统状态分布，用逆伽马分布描述测量噪声协方差；
(3)k时刻，先对每一模态下系统状态以及逆伽马分布参数进行预测，再根据k时刻的测量值对每一模态下系统状态和逆伽马分布参数进行迭代更新，输出k时刻的系统状态估计值和测量噪声协方差矩阵估计值。
本发明一个较佳实施例中，进一步包括，步骤(3)中，输出k时刻的系统状态估计值和测量噪声协方差矩阵估计值包括，
重采样获得每一模态下新的状态粒子及权重
Figure GDA0003047311250000041
并输出k时刻系统各模态下的状态估计值
Figure GDA0003047311250000042
Figure GDA0003047311250000043
其中，
Figure GDA0003047311250000044
为k时刻s模态下系统状态粒子，
Figure GDA0003047311250000045
为k时刻s模态下粒子权重，Np为粒子数，i为粒子数索引。
本发明一个较佳实施例中，进一步包括，步骤(3)中，输出k时刻的系统状态估计值和测量噪声协方差矩阵估计值还包括，
计算系统模态概率的更新值：
Figure GDA0003047311250000046
其中，
Figure GDA0003047311250000047
为系统在k时刻处于s模态下的概率，
Figure GDA0003047311250000048
m为系统模态总数；πns为模态n到模态s的转移概率；
Figure GDA0003047311250000049
为系统在k-1时刻处于n模态下的概率；
Figure GDA00030473112500000410
为系统状态粒子的预测值；
Figure GDA00030473112500000411
为k时刻s模态下粒子权重的预测值；δ为狄拉克δ函数；
Figure GDA00030473112500000412
为k时刻测量噪声协方差矩阵估计值；N(·)表示高斯分布；
根据系统模态概率的更新值，将各模态下的状态估计值进行融合得到k时刻的系统状态估计值和测量噪声协方差矩阵估计值。
本发明一个较佳实施例中，进一步包括，步骤(3)还包括，k-1时刻，使用交互式多模型算法将系统在每一模态下的状态粒子值进行融合，
Figure GDA0003047311250000051
融合后的状态粒子值作为k时刻每一模态下的状态粒子初始值；
其中，
Figure GDA0003047311250000052
为粒子在s模态下的初始权重，
Figure GDA0003047311250000053
为s模态下的融合初始粒子，
Figure GDA0003047311250000054
和
Figure GDA0003047311250000055
分别为k时刻逆伽马分布的形状参数和尺度参数的初始值。
本发明一个较佳实施例中，进一步包括，步骤(2)中，生成一组加权粒子来表示系统在k-1时刻处于n模态下的系统状态的概率分布，即
Figure GDA0003047311250000056
其中，
Figure GDA0003047311250000057
表示系统在k-1时刻处于n模态；y1:k-1为从时刻1到时刻k-1的测量值序列；Np为粒子数；
Figure GDA0003047311250000058
为粒子权重；δ为狄拉克δ函数；
Figure GDA0003047311250000059
为状态粒子，i为粒子数索引。
本发明一个较佳实施例中，进一步包括，步骤(2)中，设定系统测量噪声协方差与系统模态无关，且定义
Figure GDA0003047311250000061
其中，
Figure GDA0003047311250000062
为Rk的对角线元素，dy为yk的维度，使用逆伽马分布
Figure GDA00030473112500000613
来描述测量噪声协方差的概率密度，在k-1时刻有，
Figure GDA0003047311250000063
其中，αk-1,j和βk-1,j分别为逆伽马分布的形状参数和尺度参数，Rk-1为系统k-1时刻的测量噪声协方差矩阵，
Figure GDA0003047311250000064
为系统k-1时刻处于n模态，y1:k-1为从时刻1到时刻k-1的测量值序列，j为矩阵对角线元素索引，dy为yk的维度，
Figure GDA0003047311250000065
为k-1时刻测量噪声协方差矩阵第j个对角线元素，diag(·)为对角矩阵。
本发明一个较佳实施例中，进一步包括，步骤(3)中，使用状态预测模型对每一模态下系统状态进行预测，所述状态预测模型为，
Figure GDA0003047311250000066
Figure GDA0003047311250000067
为系统状态粒子预测值；
Figure GDA0003047311250000068
为根据过程噪声的分布产生的噪声采样粒子；
Figure GDA0003047311250000069
为s模态下的融合初始粒子；
Figure GDA00030473112500000610
表示系统k时刻的模态为s。
本发明一个较佳实施例中，进一步包括，步骤(3)中还包括计算每一模态下的粒子权重
Figure GDA00030473112500000611
计算公式为：
Figure GDA00030473112500000612
Figure GDA0003047311250000071
为k时刻s模态下粒子权重的预测值，exp(·)为以自然常数e为底的指数函数，(·)T为矩阵的转置，
Figure GDA0003047311250000072
为k时刻测量噪声协方差矩阵逆的期望。
本发明一个较佳实施例中，进一步包括，步骤(3)中还包括，判断是否满足l＝L，是则输出k时刻系统各模态下的粒子权重值和逆伽马分布的尺度参数更新值，否则l＝l+1，并计算测量噪声协方差矩阵逆的期望
Figure GDA0003047311250000073
Figure GDA0003047311250000074
Figure GDA0003047311250000075
Figure GDA0003047311250000076
其中，
Figure GDA0003047311250000077
为k时刻s模态下粒子权重，
Figure GDA0003047311250000078
为第L步迭代时粒子权重值，
Figure GDA0003047311250000079
为第L步迭代时逆伽马分布的尺度参数，βk,j为逆伽马分布的尺度参数；
Figure GDA00030473112500000710
为第l-1步迭代时逆伽马分布的形状参数，
Figure GDA00030473112500000711
为第l-1步迭代时逆伽马分布的尺度参数，且
Figure GDA00030473112500000712
本发明的有益效果：
本发明的三自由度直升机系统传感器运行状态监测方法，计算获得传感器的测量噪声协方差矩阵的估计值，若该测量噪声协方差矩阵估计值偏离所述传感器的标称值范围，判断出传感器的运行状态不正常；若该测量噪声协方差矩阵估计值包含在所述传感器的标称值范围内，判断出传感器的运行状态正常。该方法适用于三自由度直升机系统中传感器运行状态监测，通过对测量噪声协方差矩阵的准确估计，来实时监测传感器工作状态，以保证系统安全稳定运行。
附图说明
图1为三自由度直升机系统的结构图；
图2为本发明优选实施例中三自由度直升机系统传感器运行状态监测方法的流程图；
具体实施方式
下面结合附图和具体实施例对本发明作进一步说明，以使本领域的技术人员可以更好地理解本发明并能予以实施，但所举实施例不作为对本发明的限定。
实施例
本发明实施例公开一种三自由度直升机系统传感器运行状态监测方法，适用于三自由度直升机系统中传感器运行状态监测，参照图2所示，该方法包括以下步骤，
计算出三自由度直升机系统中传感器的测量噪声协方差矩阵的估计值；
分析所述传感器的测量噪声协方差矩阵估计值，判断出所述传感器的运行状态是否正常：若所述传感器的测量噪声协方差矩阵估计值偏离所述传感器的标称值范围，判断出所述传感器的运行状态不正常；若所述传感器的测量噪声协方差矩阵估计值包含在所述传感器的标称值范围内，判断出所述传感器的运行状态正常。
具体的，计算三自由度直升机系统中传感器的测量噪声协方差矩阵的估计值的方法为：
第一步、根据三自由度直升机系统的结构创建随机跳变系统模型：
根据图1所示的三自由度直升机系统的结构，定义系统状态为
Figure GDA0003047311250000091
和系统输入为
Figure GDA0003047311250000092
其中，x1为仰角，x2为俯仰角，x3为行程角，
Figure GDA0003047311250000093
为仰角角速度，
Figure GDA0003047311250000094
为俯仰角角速度，
Figure GDA0003047311250000095
为行程角角速度，Vf为前电机电压，Vb为后电机电压，那么，三自由度直升机系统可由以下离散状态方程描述：
Figure GDA0003047311250000096
其中，k为时间索引；
Figure GDA0003047311250000097
为系统状态变化量；xk-1为k-1时刻的系统状态；uk-1为k-1时刻的系统输入；ωk为过程噪声且ωk服从均值为零的高斯分布即ωk～N(0,Qk)，Qk为过程噪声协方差矩阵；系数矩阵
Figure GDA0003047311250000098
且有A1,4＝A2,5＝A3,6＝1，
Figure GDA0003047311250000099
其余元素均为零；系数矩阵
Figure GDA00030473112500000910
且有
Figure GDA00030473112500000911
B5.1＝-B5,2＝(0.5Kf)/mfLh，其余元素均为零；
Figure GDA00030473112500000912
为欧氏空间；Lw为行程轴到平衡物的距离；mw为平衡物的质量；La为行程轴到直升机机体的距离；mf为前螺旋桨的质量；Lh为螺距轴与每个电动机之间的距离；Kf为螺旋桨推力常数，且Kf的取值决定系统模态。
建立系统测量方程为：
yk＝Cxk+νk
其中，yk为k时刻的系统测量值，νk为k时刻的系统测量噪声且νk服从均值为零的高斯分布即νk～N(0,Rk)，Rk为测量噪声协方差矩阵；C为系统系数矩阵。
第二步、将上述(1)和(2)式所描述的三自由度直升机系统写为一般的随机跳变系统模型为：
xk＝fk(xk-1,rk)+ωk；
yk＝gk(xk,rk)+νk；
k为时间索引；xk为k时刻系统状态值；yk为k时刻系统测量值；fk(·)为系统状态方程；gk(·)为系统测量方程；rk为k时刻系统模态；ωk为过程噪声且ωk服从均值为零的高斯分布即ωk～N(0,Qk)，Qk为过程噪声协方差矩阵；νk为测量噪声且νk服从均值为零的高斯分布即νk～N(0,Rk)，Rk为k时刻系统测量噪声协方差矩阵；xk-1为k-1时刻系统状态值。
第三步、基于变分贝叶斯理论将每一模态下系统状态和测量噪声协方差的联合后验概率密度函数
Figure GDA0003047311250000101
用两个独立的概率密度函数
Figure GDA0003047311250000102
来表示，即，
Figure GDA0003047311250000103
Figure GDA0003047311250000111
表示系统k时刻的模态为s，y1:k为从时刻1到时刻k的测量值序列。
第四步、生成一组加权粒子来表示系统在k-1时刻处于n模态下的系统状态的概率分布，即
Figure GDA0003047311250000112
其中，
Figure GDA0003047311250000113
表示系统在k-1时刻处于n模态，y1:k-1为从时刻1到时刻k-1的测量值序列，Np为粒子数，
Figure GDA0003047311250000114
为粒子权重，δ为狄拉克δ函数，
Figure GDA0003047311250000115
为状态粒子，i为粒子数索引。
第五步：设定系统测量噪声协方差矩阵与系统模态无关，且定义为
Figure GDA0003047311250000116
其中，diag(·)为对角矩阵，
Figure GDA0003047311250000117
为Rk的对角线元素，dy为yk的维度。使用逆伽马分布
Figure GDA0003047311250000118
来描述测量噪声协方差的概率分布，在k-1时刻有
Figure GDA0003047311250000119
其中，Rk-1为系统k-1时刻的测量噪声协方差矩阵，αk-1,j和βk-1,j分别为逆伽马分布的形状参数和尺度参数，j为矩阵对角线元素索引；
Figure GDA00030473112500001110
为k-1时刻测量噪声协方差矩阵第j个对角线元素。
第六步：设定系统参数初始值
Figure GDA00030473112500001111
μ0，α0，β0，fk，gk，yk，Qk，ρ，steps，Bp，L，
Figure GDA00030473112500001112
其中，
Figure GDA00030473112500001113
为粒子初始值，
Figure GDA0003047311250000121
为粒子权重初始值，μ0为系统初始模态概率，α0为逆伽马分布的形状参数初始值，β0为逆伽马分布的尺度参数初始值，fk为系统状态方程，gk为系统测量方程，yk为系统k时刻测量值，ρ为自启发式动态模型系数，steps为总的采样系数，L为每一时刻迭代的次数，
Figure GDA0003047311250000122
为系统模态取值的有限集合，m为系统模态总数，
Figure GDA0003047311250000123
为系统转移概率。
令k＝1；
第七步：使用交互式多模型算法，将k-1时刻系统在每一模态下的状态值进行融合，可得
Figure GDA0003047311250000124
其中，
Figure GDA0003047311250000125
为粒子在s模态下的初始权重，
Figure GDA0003047311250000126
为s模态下的融合初始粒子，
Figure GDA0003047311250000127
和
Figure GDA0003047311250000128
分别为k时刻逆伽马分布的形状参数和尺度参数的初始值。
第八步：对系统状态以及逆伽马分布参数进行预测，状态的预测计算公式为：
Figure GDA0003047311250000129
其中，
Figure GDA00030473112500001210
为系统状态粒子预测值，
Figure GDA00030473112500001211
为根据过程噪声的分布产生的噪声采样粒子，逆伽马分布参数的预测模型为自启发式动态模型，即
Figure GDA0003047311250000131
Figure GDA0003047311250000132
其中，
Figure GDA0003047311250000133
为逆伽马分布的形状参数的预测值，
Figure GDA0003047311250000134
为逆伽马分布的尺度参数的预测值，自启发式动态模型系数ρ∈(0,1]。
第九步：更新逆伽马分布的形状参数α，更新公式为：
Figure GDA0003047311250000135
其中，αk,j为逆伽马分布的形状参数的更新值；
令l＝1；
第十步：计算测量噪声协方差矩阵逆的期望
Figure GDA0003047311250000136
计算公式为：
Figure GDA0003047311250000137
其中，
Figure GDA0003047311250000138
为第l-1步迭代时逆伽马分布的形状参数，
Figure GDA0003047311250000139
为第l-1步迭代时逆伽马分布的尺度参数，且
Figure GDA00030473112500001310
第十一步：计算每一模态下的粒子权重
Figure GDA00030473112500001311
计算公式为：
Figure GDA00030473112500001312
其中，
Figure GDA0003047311250000141
为k时刻s模态下粒子权重的预测值，exp(·)为以自然常数e为底的指数函数，(·)T为矩阵的转置。
第十二步：更新逆伽马分布的尺度参数β，更新公式为：
Figure GDA0003047311250000142
其中，
Figure GDA0003047311250000143
为第l步迭代时逆伽马分布的尺度参数，
Figure GDA0003047311250000144
πns为模态n到模态s的转移概率，
Figure GDA0003047311250000145
表示系统在k-1时刻处于n模态下的概率，
Figure GDA0003047311250000146
表示矩阵对角线元素的平方。
判断是否满足l＝L，是则执行下一步，否则l＝l+1，并跳转到第六步；
第十三步：输出k时刻系统各模态下的粒子权重值以及逆伽马分布的尺度参数的更新值，即
Figure GDA0003047311250000147
Figure GDA0003047311250000148
其中，
Figure GDA0003047311250000149
为k时刻s模态下粒子权重，
Figure GDA00030473112500001410
为第L步迭代时粒子权重值，
Figure GDA00030473112500001411
为第L步迭代时逆伽马分布的尺度参数。
第十四步：进行重采样得到每一模态下新的状态粒子及粒子权重：
Figure GDA0003047311250000151
Figure GDA0003047311250000152
其中，
Figure GDA0003047311250000153
为k时刻s模态下系统状态粒子，
Figure GDA0003047311250000154
为k时刻s模态下粒子权重；
第十五步：输出k时刻系统各模态下的状态估计值：
Figure GDA0003047311250000155
其中，
Figure GDA0003047311250000156
为k时刻s模态下系统状态估计值；
第十六步：计算系统模态概率的更新值，计算公式为：
Figure GDA0003047311250000157
其中，
Figure GDA0003047311250000158
为系统在k时刻处于s模态下的概率，
似然
Figure GDA0003047311250000159
Figure GDA00030473112500001510
为k时刻测量噪声协方差矩阵估计值；
第十七步：将各模态下的系统状态值进行融合得到k时刻的系统状态估计值
Figure GDA00030473112500001511
以及测量噪声协方差矩阵估计值
Figure GDA00030473112500001512
即
Figure GDA0003047311250000161
Figure GDA0003047311250000162
第十八步：判断是否满足k＝steps，是则结束；否则k＝k+1，并跳转至第三步。
参照图1所示的三自由度直升机系统，系统参数为Lw＝0.47m，mw＝1.87kg，La＝0.66m，mf＝0.713kg，Lh＝0.178m，Vf＝Vb＝6V，选择采样时间为Ts＝0.1s，系统有三种模态：Kf＝0.1188为模态1，Kf＝0.25为模态2，Kf＝0.36为模态3，C＝[I3×3 03×3]，逆伽马分布参数的初始值为α0＝[3,3,2]和β0＝[0.1,2,1]，自启发式动态模型的系数为ρ＝0.99，每一时刻迭代的次数为L＝3，状态初始值x0＝[0 0 0 0 0 0]T，系统的转移概率矩阵为：
Figure GDA0003047311250000163
通过仿真可以发现，本发明所提出的三自由度直升机系统传感器运行状态监测方法，可以实时跟踪变化的测量噪声协方差矩阵，且该方法的均方根误差仅为0.75左右，并且通过分析测量噪声协方差矩阵的准确估计值可以实时监测系统传感器的工作状态，从而为系统决策制定提供有效信息。
以上所述实施例仅是为充分说明本发明而所举的较佳的实施例，本发明的保护范围不限于此。本技术领域的技术人员在本发明基础上所作的等同替代或变换，均在本发明的保护范围之内。本发明的保护范围以权利要求书为准。

Claims (7)

1.一种三自由度直升机系统传感器运行状态监测方法，其特征在于：包括以下步骤，
计算出传感器的测量噪声协方差矩阵的估计值；
分析所述传感器的测量噪声协方差矩阵估计值，判断出所述传感器的运行状态是否正常：若所述传感器的测量噪声协方差矩阵估计值偏离所述传感器的标称值范围，判断出所述传感器的运行状态不正常；若所述传感器的测量噪声协方差矩阵估计值包含在所述传感器的标称值范围内，判断出所述传感器的运行状态正常；
计算传感器的测量噪声协方差矩阵的估计值包括以下步骤，
(1)创建所述传感器所在的三自由度直升机系统的随机跳变系统模型，
xk＝fk(xk-1,rk)+ωk；
yk＝gk(xk,rk)+νk；
k为时间索引；xk为k时刻系统状态值；yk为k时刻系统测量值；fk(·)为系统状态方程；gk(·)为系统测量方程；rk为k时刻系统模态；ωk为过程噪声且ωk服从均值为零的高斯分布即ωk～N(0,Qk)，Qk为过程噪声协方差矩阵；νk为测量噪声且νk服从均值为零的高斯分布即νk～N(0,Rk)，Rk为k时刻系统测量噪声协方差矩阵；xk-1为k-1时刻系统状态值；
(2)基于变分贝叶斯理论将每一模态下系统状态和测量噪声协方差的联合后验概率密度函数
Figure FDA0003047311240000021
用两个独立的概率密度函数
Figure FDA0003047311240000022
来表示，即
Figure FDA0003047311240000023
Figure FDA0003047311240000024
表示系统k时刻的模态为s，y1:k为从时刻1到时刻k的测量值序列；
用一组加权粒子描述系统状态分布，用逆伽马分布描述测量噪声协方差；
(3)k时刻，先对每一模态下系统状态以及逆伽马分布参数进行预测，再根据k时刻的测量值对每一模态下系统状态和逆伽马分布参数进行迭代更新，输出k时刻的系统状态估计值和测量噪声协方差矩阵估计值；
步骤(3)中，输出k时刻的系统状态估计值和测量噪声协方差矩阵估计值包括，
重采样获得每一模态下新的状态粒子及权重
Figure FDA0003047311240000025
并输出k时刻系统各模态下的状态估计值
Figure FDA0003047311240000026
Figure FDA0003047311240000027
其中，
Figure FDA0003047311240000028
为k时刻s模态下系统状态粒子，
Figure FDA0003047311240000029
为k时刻s模态下粒子权重，Np为粒子数，i为粒子数索引；
步骤(3)中，输出k时刻的系统状态估计值和测量噪声协方差矩阵估计值还包括，
计算系统模态概率的更新值：
Figure FDA0003047311240000031
其中，
Figure FDA0003047311240000032
为系统在k时刻处于s模态下的概率，
Figure FDA0003047311240000033
m为系统模态总数；πns为模态n到模态s的转移概率；
Figure FDA0003047311240000034
为系统在k-1时刻处于n模态下的概率；
Figure FDA0003047311240000035
为系统状态粒子的预测值；
Figure FDA0003047311240000036
为k时刻s模态下粒子权重的预测值；δ为狄拉克δ函数；
Figure FDA0003047311240000037
为k时刻测量噪声协方差矩阵估计值；N(·)表示高斯分布；
根据系统模态概率的更新值，将各模态下的状态估计值进行融合得到k时刻的系统状态估计值和测量噪声协方差矩阵估计值，能够确保状态估计值的准确度。
2.如权利要求1所述的三自由度直升机系统传感器运行状态监测方法，其特征在于：步骤(3)还包括，k-1时刻，使用交互式多模型算法将系统在每一模态下的状态粒子值进行融合，
Figure FDA0003047311240000038
融合后的状态粒子值作为k时刻每一模态下的状态粒子初始值；
其中，
Figure FDA0003047311240000041
为粒子在s模态下的初始权重，
Figure FDA0003047311240000042
为s模态下的融合初始粒子，
Figure FDA0003047311240000043
和
Figure FDA0003047311240000044
分别为k时刻逆伽马分布的形状参数和尺度参数的初始值；
Figure FDA0003047311240000045
为k-1时刻测量噪声协方差矩阵第j个对角线元素；j为矩阵对角线元素索引；Np为粒子数；dy为yk的维度。
3.如权利要求1所述的三自由度直升机系统传感器运行状态监测方法，其特征在于：步骤(2)中，生成一组加权粒子来表示系统在k-1时刻处于n模态下的系统状态的概率分布，即
Figure FDA0003047311240000046
其中，
Figure FDA0003047311240000047
表示系统在k-1时刻处于n模态；y1:k-1为从时刻1到时刻k-1的测量值序列；Np为粒子数；
Figure FDA0003047311240000048
为粒子权重；δ为狄拉克δ函数；
Figure FDA0003047311240000049
为状态粒子，i为粒子数索引。
4.如权利要求1所述的三自由度直升机系统传感器运行状态监测方法，其特征在于：步骤(2)中，设定系统测量噪声协方差与系统模态无关，且定义
Figure FDA00030473112400000410
其中，
Figure FDA00030473112400000411
为Rk的对角线元素，dy为yk的维度，使用逆伽马分布
Figure FDA00030473112400000412
来描述测量噪声协方差的概率密度，在k-1时刻有，
Figure FDA00030473112400000413
其中，αk-1,j和βk-1,j分别为逆伽马分布的形状参数和尺度参数，Rk-1为系统k-1时刻的测量噪声协方差矩阵，
Figure FDA0003047311240000051
为系统k-1时刻处于n模态，y1:k-1为从时刻1到时刻k-1的测量值序列，j为矩阵对角线元素索引，dy为yk的维度，
Figure FDA0003047311240000052
为k-1时刻测量噪声协方差矩阵第j个对角线元素，diag(·)为对角矩阵。
5.如权利要求1所述的三自由度直升机系统传感器运行状态监测方法，其特征在于：步骤(3)中，使用状态预测模型对每一模态下系统状态进行预测，所述状态预测模型为，
Figure FDA0003047311240000053
Figure FDA0003047311240000054
为系统状态粒子预测值；
Figure FDA0003047311240000055
为根据过程噪声的分布产生的噪声采样粒子；
Figure FDA0003047311240000056
为s模态下的融合初始粒子；
Figure FDA0003047311240000057
表示系统k时刻的模态为s。
6.如权利要求1所述的三自由度直升机系统传感器运行状态监测方法，其特征在于：步骤(3)中还包括计算每一模态下的粒子权重
Figure FDA0003047311240000058
计算公式为：
Figure FDA0003047311240000059
Figure FDA00030473112400000510
为k时刻s模态下粒子权重的预测值，exp(·)为以自然常数e为底的指数函数，(·)T为矩阵的转置，
Figure FDA00030473112400000511
为k时刻测量噪声协方差矩阵逆的期望；
Figure FDA00030473112400000512
为系统状态粒子预测值。
7.如权利要求6所述的三自由度直升机系统传感器运行状态监测方法，其特征在于：步骤(3)中还包括，判断是否满足l＝L，是则输出k时刻系统各模态下的粒子权重值和逆伽马分布的尺度参数更新值，否则l＝l+1，并计算测量噪声协方差矩阵逆的期望
Figure FDA0003047311240000061
Figure FDA0003047311240000062
Figure FDA0003047311240000063
Figure FDA0003047311240000064
其中，
Figure FDA0003047311240000065
为k时刻s模态下粒子权重，
Figure FDA0003047311240000066
为第L步迭代时粒子权重值，
Figure FDA0003047311240000067
为第L步迭代时逆伽马分布的尺度参数，βk,j为逆伽马分布的尺度参数；dy为yk的维度；
Figure FDA0003047311240000068
为第l-1步迭代时逆伽马分布的形状参数，
Figure FDA0003047311240000069
为第l-1步迭代时逆伽马分布的尺度参数，且
Figure FDA00030473112400000610
αk,j为逆伽马分布的形状参数的更新值；
Figure FDA00030473112400000611
为逆伽马分布的尺度参数的预测值。
CN202010524010.0A 2020-06-10 2020-06-10 三自由度直升机系统传感器运行状态监测方法 Active CN111609878B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010524010.0A CN111609878B (zh) 2020-06-10 2020-06-10 三自由度直升机系统传感器运行状态监测方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010524010.0A CN111609878B (zh) 2020-06-10 2020-06-10 三自由度直升机系统传感器运行状态监测方法

Publications (2)

Publication Number Publication Date
CN111609878A CN111609878A (zh) 2020-09-01
CN111609878B true CN111609878B (zh) 2021-06-22

Family

ID=72204784

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010524010.0A Active CN111609878B (zh) 2020-06-10 2020-06-10 三自由度直升机系统传感器运行状态监测方法

Country Status (1)

Country Link
CN (1) CN111609878B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115070765B (zh) * 2022-06-27 2023-06-13 江南大学 一种基于变分推断的机器人状态估计方法及系统

Family Cites Families (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103048991A (zh) * 2013-01-08 2013-04-17 南京航空航天大学 三自由度直升机系统故障估计观测器及容错控制方法
CN105182990B (zh) * 2015-09-30 2019-08-20 南京航空航天大学 具有输出受限的三自由度模型直升机的鲁棒控制方法
CN105716844B (zh) * 2016-01-30 2018-06-08 西北工业大学 建立机电作动器的卡尔曼滤波模型及故障诊断方法
CN107729585B (zh) * 2016-08-12 2020-08-28 贵州火星探索科技有限公司 一种对无人机的噪声协方差进行估算的方法
CN108036812A (zh) * 2017-11-13 2018-05-15 深圳市易成自动驾驶技术有限公司 传感器状态检测方法、装置及计算机可读存储介质
CN110957011B (zh) * 2019-11-25 2023-03-17 江南大学 连续搅拌反应器在未知时变测量噪声下的在线生产参数估计方法

Also Published As

Publication number Publication date
CN111609878A (zh) 2020-09-01

Similar Documents

Publication Publication Date Title
CN112247992B (zh) 一种机器人前馈力矩补偿方法
Si et al. Estimating remaining useful life with three-source variability in degradation modeling
Wen et al. Multiple-phase modeling of degradation signal for condition monitoring and remaining useful life prediction
Liu et al. Comparison of 4 numerical solvers for stiff and hybrid systems simulation
CN111143981B (zh) 虚拟试验模型验证系统及方法
CN113221263A (zh) 一种考虑分布参数不确定性的机械产品结构失效优化方法
CN111609878B (zh) 三自由度直升机系统传感器运行状态监测方法
CN106909738B (zh) 一种模型参数辨识方法
JP2019125102A (ja) 流体解析装置、流体解析方法、及び流体解析プログラム
CN115979310B (zh) 一种惯导系统性能退化评估方法、系统、电子设备及介质
US7194394B2 (en) Method and apparatus for detecting and correcting inaccuracies in curve-fitted models
CN114970341B (zh) 基于机器学习的低轨卫星轨道预报精度提升模型建立方法
CN115270239A (zh) 基于动力特性和智能算法响应面法的桥梁可靠性预测方法
CN113340324A (zh) 一种基于深度确定性策略梯度的视觉惯性自校准方法
CN113236506B (zh) 一种基于滤波的工业时滞系统故障检测方法
Khodaparast et al. Fuzzy model updating and its application to the DLR AIRMOD test structure
CN111581905B (zh) 隧道二极管电路系统在未知测量噪声下的状态估计方法
CN113408057A (zh) 一种基于em算法的飞机气流角高精度离线估计方法
CN110175740B (zh) 一种基于Kriging代理模型的场景分析实现方法
Kukreja et al. Nonlinear black-box modeling of aeroelastic systems using structure detection approach: application to F/A-18 aircraft data
Kaur et al. Comparative analysis of the software effort estimation models
JP2006195542A (ja) モデル同定装置およびモデル同定プログラム
CN114858200B (zh) 车辆传感器检测到的对象的质量评价方法及装置
CN113886947B (zh) 一种基于迭代策略的飞行器静气弹系统输出状态量区间确定方法
CN112416774B (zh) 一种添加权重的软件可靠性测试方法

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