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

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

Info

Publication number
CN111609878A
CN111609878A CN202010524010.0A CN202010524010A CN111609878A CN 111609878 A CN111609878 A CN 111609878A CN 202010524010 A CN202010524010 A CN 202010524010A CN 111609878 A CN111609878 A CN 111609878A
Authority
CN
China
Prior art keywords
time
mode
sensor
state
value
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
CN202010524010.0A
Other languages
English (en)
Other versions
CN111609878B (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 BDA0002533164490000031
用两个独立的概率密度函数
Figure BDA0002533164490000032
Figure BDA0002533164490000033
来表示，即
Figure BDA0002533164490000034
Figure BDA0002533164490000035
表示系统k时刻的模态为s，y1:k为从时刻1到时刻k的测量值序列；
用一组加权粒子描述系统状态分布，用逆伽马分布描述测量噪声协方差；
(3)k时刻，先对每一模态下系统状态以及逆伽马分布参数进行预测，再根据k时刻的测量值对每一模态下系统状态和逆伽马分布参数进行迭代更新，输出k时刻的系统状态估计值和测量噪声协方差矩阵估计值。
本发明一个较佳实施例中，进一步包括，步骤(3)中，输出k时刻的系统状态估计值和测量噪声协方差矩阵估计值包括，
重采样获得每一模态下新的状态粒子及权重
Figure BDA0002533164490000036
并输出k时刻系统各模态下的状态估计值
Figure BDA0002533164490000037
Figure BDA0002533164490000038
其中，
Figure BDA0002533164490000039
为k时刻s模态下系统状态粒子，
Figure BDA00025331644900000310
为k时刻s模态下粒子权重，Np为粒子数，i为粒子数索引。
本发明一个较佳实施例中，进一步包括，步骤(3)中，输出k时刻的系统状态估计值和测量噪声协方差矩阵估计值还包括，
计算系统模态概率的更新值：
Figure BDA0002533164490000041
其中，
Figure BDA0002533164490000042
为系统在k时刻处于s模态下的概率，
Figure BDA0002533164490000043
m为系统模态总数；πns为模态n到模态s的转移概率；
Figure BDA0002533164490000044
为系统在k-1时刻处于n模态下的概率；
Figure BDA0002533164490000045
为系统状态粒子的预测值；
Figure BDA0002533164490000046
为k时刻s模态下粒子权重的预测值；δ为狄拉克δ函数；
Figure BDA0002533164490000047
为k时刻测量噪声协方差矩阵估计值；N(·)表示高斯分布；
根据系统模态概率的更新值，将各模态下的状态估计值进行融合得到k时刻的系统状态估计值和测量噪声协方差矩阵估计值。
本发明一个较佳实施例中，进一步包括，步骤(3)还包括，k-1时刻，使用交互式多模型算法将系统在每一模态下的状态粒子值进行融合，
Figure BDA0002533164490000048
融合后的状态粒子值作为k时刻每一模态下的状态粒子初始值；
其中，
Figure BDA0002533164490000049
为粒子在s模态下的初始权重，
Figure BDA00025331644900000410
为s模态下的融合初始粒子，
Figure BDA00025331644900000411
和
Figure BDA00025331644900000412
分别为k时刻逆伽马分布的形状参数和尺度参数的初始值。
本发明一个较佳实施例中，进一步包括，步骤(2)中，生成一组加权粒子来表示系统在k-1时刻处于n模态下的系统状态的概率分布，即
Figure BDA0002533164490000051
其中，
Figure BDA0002533164490000052
表示系统在k-1时刻处于n模态；y1:k-1为从时刻1到时刻k-1的测量值序列；Np为粒子数；
Figure BDA0002533164490000053
为粒子权重；δ为狄拉克δ函数；
Figure BDA0002533164490000054
为状态粒子，i为粒子数索引。
本发明一个较佳实施例中，进一步包括，步骤(2)中，设定系统测量噪声协方差与系统模态无关，且定义
Figure BDA0002533164490000055
其中，
Figure BDA0002533164490000056
为Rk的对角线元素，dy为yk的维度，使用逆伽马分布
Figure BDA00025331644900000514
来描述测量噪声协方差的概率密度，在k-1时刻有，
Figure BDA0002533164490000057
其中，αk-1,j和βk-1,j分别为逆伽马分布的形状参数和尺度参数，Rk-1为系统k-1时刻的测量噪声协方差矩阵，
Figure BDA00025331644900000515
为系统k-1时刻处于n模态，y1:k-1为从时刻1到时刻k-1的测量值序列，j为矩阵对角线元素索引，dy为yk的维度，
Figure BDA00025331644900000516
为k-1时刻测量噪声协方差矩阵第j个对角线元素，diag(·)为对角矩阵。
本发明一个较佳实施例中，进一步包括，步骤(3)中，使用状态预测模型对每一模态下系统状态进行预测，所述状态预测模型为，
Figure BDA0002533164490000058
Figure BDA0002533164490000059
为系统状态粒子预测值；
Figure BDA00025331644900000510
为根据过程噪声的分布产生的噪声采样粒子；
Figure BDA00025331644900000511
为s模态下的融合初始粒子；
Figure BDA00025331644900000512
表示系统k时刻的模态为s。
本发明一个较佳实施例中，进一步包括，步骤(3)中还包括计算每一模态下的粒子权重
Figure BDA00025331644900000513
计算公式为：
Figure BDA0002533164490000061
Figure BDA0002533164490000062
为k时刻s模态下粒子权重的预测值，exp(·)为以自然常数e为底的指数函数，(·)T为矩阵的转置，
Figure BDA0002533164490000063
为k时刻测量噪声协方差矩阵逆的期望。
本发明一个较佳实施例中，进一步包括，步骤(3)中还包括，判断是否满足l＝L，是则输出k时刻系统各模态下的粒子权重值和逆伽马分布的尺度参数更新值，否则l＝l+1，并计算测量噪声协方差矩阵逆的期望
Figure BDA0002533164490000064
Figure BDA0002533164490000065
Figure BDA0002533164490000066
Figure BDA0002533164490000067
其中，
Figure BDA0002533164490000068
为k时刻s模态下粒子权重，
Figure BDA0002533164490000069
为第L步迭代时粒子权重值，
Figure BDA00025331644900000610
为第L步迭代时逆伽马分布的尺度参数，βk,j为逆伽马分布的尺度参数；
Figure BDA00025331644900000611
为第l-1步迭代时逆伽马分布的形状参数，
Figure BDA00025331644900000612
为第l-1步迭代时逆伽马分布的尺度参数，且
Figure BDA00025331644900000613
本发明的有益效果：
本发明的三自由度直升机系统传感器运行状态监测方法，计算获得传感器的测量噪声协方差矩阵的估计值，若该测量噪声协方差矩阵估计值偏离所述传感器的标称值范围，判断出传感器的运行状态不正常；若该测量噪声协方差矩阵估计值包含在所述传感器的标称值范围内，判断出传感器的运行状态正常。该方法适用于三自由度直升机系统中传感器运行状态监测，通过对测量噪声协方差矩阵的准确估计，来实时监测传感器工作状态，以保证系统安全稳定运行。
附图说明
图1为三自由度直升机系统的结构图；
图2为本发明优选实施例中三自由度直升机系统传感器运行状态监测方法的流程图；
具体实施方式
下面结合附图和具体实施例对本发明作进一步说明，以使本领域的技术人员可以更好地理解本发明并能予以实施，但所举实施例不作为对本发明的限定。
实施例
本发明实施例公开一种三自由度直升机系统传感器运行状态监测方法，适用于三自由度直升机系统中传感器运行状态监测，参照图2所示，该方法包括以下步骤，
计算出三自由度直升机系统中传感器的测量噪声协方差矩阵的估计值；
分析所述传感器的测量噪声协方差矩阵估计值，判断出所述传感器的运行状态是否正常：若所述传感器的测量噪声协方差矩阵估计值偏离所述传感器的标称值范围，判断出所述传感器的运行状态不正常；若所述传感器的测量噪声协方差矩阵估计值包含在所述传感器的标称值范围内，判断出所述传感器的运行状态正常。
具体的，计算三自由度直升机系统中传感器的测量噪声协方差矩阵的估计值的方法为：
第一步、根据三自由度直升机系统的结构创建随机跳变系统模型：
根据图1所示的三自由度直升机系统的结构，定义系统状态为
Figure BDA0002533164490000081
和系统输入为
Figure BDA0002533164490000082
其中，x1为仰角，x2为俯仰角，x3为行程角，
Figure BDA0002533164490000083
为仰角角速度，
Figure BDA0002533164490000084
为俯仰角角速度，
Figure BDA0002533164490000085
为行程角角速度，Vf为前电机电压，Vb为后电机电压，那么，三自由度直升机系统可由以下离散状态方程描述：
Figure BDA0002533164490000086
其中，k为时间索引；
Figure BDA0002533164490000087
为系统状态变化量；xk-1为k-1时刻的系统状态；uk-1为k-1时刻的系统输入；ωk为过程噪声且ωk服从均值为零的高斯分布即ωk～N(0,Qk)，Qk为过程噪声协方差矩阵；系数矩阵
Figure BDA0002533164490000088
且有A1,4＝A2,5＝A3,6＝1，
Figure BDA0002533164490000089
其余元素均为零；系数矩阵
Figure BDA00025331644900000810
且有
Figure BDA00025331644900000811
B5.1＝-B5,2＝(0.5Kf)/mfLh，其余元素均为零；
Figure BDA00025331644900000812
为欧氏空间；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 BDA0002533164490000091
用两个独立的概率密度函数
Figure BDA0002533164490000092
Figure BDA0002533164490000093
来表示，即，
Figure BDA0002533164490000094
Figure BDA00025331644900000912
表示系统k时刻的模态为s，y1:k为从时刻1到时刻k的测量值序列。
第四步、生成一组加权粒子来表示系统在k-1时刻处于n模态下的系统状态的概率分布，即
Figure BDA0002533164490000095
其中，
Figure BDA0002533164490000096
表示系统在k-1时刻处于n模态，y1:k-1为从时刻1到时刻k-1的测量值序列，Np为粒子数，
Figure BDA0002533164490000097
为粒子权重，δ为狄拉克δ函数，
Figure BDA0002533164490000098
为状态粒子，i为粒子数索引。
第五步：设定系统测量噪声协方差矩阵与系统模态无关，且定义为
Figure BDA0002533164490000099
其中，diag(·)为对角矩阵，
Figure BDA00025331644900000910
为Rk的对角线元素，dy为yk的维度。使用逆伽马分布
Figure BDA00025331644900000911
来描述测量噪声协方差的概率分布，在k-1时刻有
Figure BDA0002533164490000101
其中，Rk-1为系统k-1时刻的测量噪声协方差矩阵，αk-1,j和βk-1,j分别为逆伽马分布的形状参数和尺度参数，j为矩阵对角线元素索引；
Figure BDA0002533164490000102
为k-1时刻测量噪声协方差矩阵第j个对角线元素。
第六步：设定系统参数初始值
Figure BDA0002533164490000103
μ0，α0，β0，fk，gk，yk，Qk，ρ，steps，Np，LL，
Figure BDA00025331644900001015
其中，
Figure BDA0002533164490000105
为粒子初始值，
Figure BDA0002533164490000106
为粒子权重初始值，μ0为系统初始模态概率，α0为逆伽马分布的形状参数初始值，β0为逆伽马分布的尺度参数初始值，fk为系统状态方程，gk为系统测量方程，yk为系统k时刻测量值，ρ为自启发式动态模型系数，steps为总的采样系数，L为每一时刻迭代的次数，
Figure BDA0002533164490000107
为系统模态取值的有限集合，m为系统模态总数，
Figure BDA0002533164490000108
为系统转移概率。
令k＝1；
第七步：使用交互式多模型算法，将k-1时刻系统在每一模态下的状态值进行融合，可得
Figure BDA0002533164490000109
其中，
Figure BDA00025331644900001010
为粒子在s模态下的初始权重，
Figure BDA00025331644900001011
为s模态下的融合初始粒子，
Figure BDA00025331644900001012
和
Figure BDA00025331644900001013
分别为k时刻逆伽马分布的形状参数和尺度参数的初始值。
第八步：对系统状态以及逆伽马分布参数进行预测，状态的预测计算公式为：
Figure BDA00025331644900001014
其中，
Figure BDA0002533164490000111
为系统状态粒子预测值，
Figure BDA0002533164490000112
为根据过程噪声的分布产生的噪声采样粒子，逆伽马分布参数的预测模型为自启发式动态模型，即
Figure BDA0002533164490000113
Figure BDA0002533164490000114
其中，
Figure BDA0002533164490000115
为逆伽马分布的形状参数的预测值，
Figure BDA0002533164490000116
为逆伽马分布的尺度参数的预测值，自启发式动态模型系数ρ∈(0,1]。
第九步：更新逆伽马分布的形状参数α，更新公式为：
Figure BDA0002533164490000117
其中，αk,j为逆伽马分布的形状参数的更新值；
令l＝1；
第十步：计算测量噪声协方差矩阵逆的期望
Figure BDA0002533164490000118
计算公式为：
Figure BDA0002533164490000119
其中，
Figure BDA00025331644900001110
为第l-1步迭代时逆伽马分布的形状参数，
Figure BDA00025331644900001111
为第l-1步迭代时逆伽马分布的尺度参数，且
Figure BDA00025331644900001112
第十一步：计算每一模态下的粒子权重
Figure BDA00025331644900001113
计算公式为：
Figure BDA00025331644900001114
其中，
Figure BDA0002533164490000121
为k时刻s模态下粒子权重的预测值，exp(·)为以自然常数e为底的指数函数，(·)T为矩阵的转置。
第十二步：更新逆伽马分布的尺度参数β，更新公式为：
Figure BDA0002533164490000122
其中，
Figure BDA0002533164490000123
为第l步迭代时逆伽马分布的尺度参数，
Figure BDA0002533164490000124
πns为模态n到模态s的转移概率，
Figure BDA0002533164490000125
表示系统在k-1时刻处于n模态下的概率，
Figure BDA0002533164490000126
表示矩阵对角线元素的平方。
判断是否满足l＝L，是则执行下一步，否则l＝l+1，并跳转到第六步；
第十三步：输出k时刻系统各模态下的粒子权重值以及逆伽马分布的尺度参数的更新值，即
Figure BDA0002533164490000127
Figure BDA0002533164490000128
其中，
Figure BDA0002533164490000129
为k时刻s模态下粒子权重，
Figure BDA00025331644900001210
为第L步迭代时粒子权重值，
Figure BDA00025331644900001211
为第L步迭代时逆伽马分布的尺度参数。
第十四步：进行重采样得到每一模态下新的状态粒子及粒子权重：
Figure BDA00025331644900001212
Figure BDA00025331644900001213
其中，
Figure BDA0002533164490000131
为k时刻s模态下系统状态粒子，
Figure BDA0002533164490000132
为k时刻s模态下粒子权重；
第十五步：输出k时刻系统各模态下的状态估计值：
Figure BDA0002533164490000133
其中，
Figure BDA0002533164490000134
为k时刻s模态下系统状态估计值；
第十六步：计算系统模态概率的更新值，计算公式为：
Figure BDA0002533164490000135
其中，
Figure BDA0002533164490000136
为系统在k时刻处于s模态下的概率，
似然
Figure BDA0002533164490000137
Figure BDA0002533164490000138
为k时刻测量噪声协方差矩阵估计值；
第十七步：将各模态下的系统状态值进行融合得到k时刻的系统状态估计值
Figure BDA0002533164490000139
以及测量噪声协方差矩阵估计值
Figure BDA00025331644900001310
即
Figure BDA00025331644900001311
Figure BDA00025331644900001312
第十八步：判断是否满足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×303×3]，逆伽马分布参数的初始值为α0＝[3,3,2]和β0＝[0.1,2,1]，自启发式动态模型的系数为ρ＝0.99，每一时刻迭代的次数为L＝3，状态初始值x0＝[000000]T，系统的转移概率矩阵为：
Figure BDA0002533164490000141
通过仿真可以发现，本发明所提出的三自由度直升机系统传感器运行状态监测方法，可以实时跟踪变化的测量噪声协方差矩阵，且该方法的均方根误差仅为0.75左右，并且通过分析测量噪声协方差矩阵的准确估计值可以实时监测系统传感器的工作状态，从而为系统决策制定提供有效信息。
以上所述实施例仅是为充分说明本发明而所举的较佳的实施例，本发明的保护范围不限于此。本技术领域的技术人员在本发明基础上所作的等同替代或变换，均在本发明的保护范围之内。本发明的保护范围以权利要求书为准。

Claims (10)

1.一种三自由度直升机系统传感器运行状态监测方法，其特征在于：包括以下步骤，
计算出传感器的测量噪声协方差矩阵的估计值；
分析所述传感器的测量噪声协方差矩阵估计值，判断出所述传感器的运行状态是否正常：若所述传感器的测量噪声协方差矩阵估计值偏离所述传感器的标称值范围，判断出所述传感器的运行状态不正常；若所述传感器的测量噪声协方差矩阵估计值包含在所述传感器的标称值范围内，判断出所述传感器的运行状态正常。
2.如权利要求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 FDA0002533164480000011
用两个独立的概率密度函数
Figure FDA0002533164480000012
Figure FDA0002533164480000021
来表示，即
Figure FDA0002533164480000022
Figure FDA0002533164480000023
表示系统k时刻的模态为s，y1:k为从时刻1到时刻k的测量值序列；
用一组加权粒子描述系统状态分布，用逆伽马分布描述测量噪声协方差；
(3)k时刻，先对每一模态下系统状态以及逆伽马分布参数进行预测，再根据k时刻的测量值对每一模态下系统状态和逆伽马分布参数进行迭代更新，输出k时刻的系统状态估计值和测量噪声协方差矩阵估计值。
3.如权利要求2所述的三自由度直升机系统传感器运行状态监测方法，其特征在于：步骤(3)中，输出k时刻的系统状态估计值和测量噪声协方差矩阵估计值包括，
重采样获得每一模态下新的状态粒子及权重
Figure FDA0002533164480000024
并输出k时刻系统各模态下的状态估计值
Figure FDA0002533164480000025
Figure FDA0002533164480000026
其中，
Figure FDA0002533164480000027
为k时刻s模态下系统状态粒子，
Figure FDA0002533164480000028
为k时刻s模态下粒子权重，Np为粒子数，i为粒子数索引。
4.如权利要求3所述的三自由度直升机系统传感器运行状态监测方法，其特征在于：步骤(3)中，输出k时刻的系统状态估计值和测量噪声协方差矩阵估计值还包括，
计算系统模态概率的更新值：
Figure FDA0002533164480000031
其中，
Figure FDA0002533164480000032
为系统在k时刻处于s模态下的概率，
Figure FDA0002533164480000033
m为系统模态总数；πns为模态n到模态s的转移概率；
Figure FDA0002533164480000034
为系统在k-1时刻处于n模态下的概率；
Figure FDA0002533164480000035
为系统状态粒子的预测值；
Figure FDA0002533164480000036
为k时刻s模态下粒子权重的预测值；δ为狄拉克δ函数；
Figure FDA0002533164480000037
为k时刻测量噪声协方差矩阵估计值；N(·)表示高斯分布；
根据系统模态概率的更新值，将各模态下的状态估计值进行融合得到k时刻的系统状态估计值和测量噪声协方差矩阵估计值。
5.如权利要求2所述的三自由度直升机系统传感器运行状态监测方法，其特征在于：步骤(3)还包括，k-1时刻，使用交互式多模型算法将系统在每一模态下的状态粒子值进行融合，
Figure FDA0002533164480000038
融合后的状态粒子值作为k时刻每一模态下的状态粒子初始值；
其中，
Figure FDA0002533164480000039
为粒子在s模态下的初始权重，
Figure FDA00025331644800000310
为s模态下的融合初始粒子，
Figure FDA00025331644800000311
和
Figure FDA00025331644800000312
分别为k时刻逆伽马分布的形状参数和尺度参数的初始值。
6.如权利要求2所述的三自由度直升机系统传感器运行状态监测方法，其特征在于：步骤(2)中，生成一组加权粒子来表示系统在k-1时刻处于n模态下的系统状态的概率分布，即
Figure FDA0002533164480000041
其中，
Figure FDA0002533164480000042
表示系统在k-1时刻处于n模态；y1:k-1为从时刻1到时刻k-1的测量值序列；Np为粒子数；
Figure FDA0002533164480000043
为粒子权重；δ为狄拉克δ函数；
Figure FDA0002533164480000044
为状态粒子，i为粒子数索引。
7.如权利要求2所述的三自由度直升机系统传感器运行状态监测方法，其特征在于：步骤(2)中，设定系统测量噪声协方差与系统模态无关，且定义
Figure FDA0002533164480000045
其中，
Figure FDA0002533164480000046
为Rk的对角线元素，dy为yk的维度，使用逆伽马分布
Figure FDA0002533164480000047
来描述测量噪声协方差的概率密度，在k-1时刻有，
Figure FDA0002533164480000048
其中，αk-1,j和βk-1,j分别为逆伽马分布的形状参数和尺度参数，Rk-1为系统k-1时刻的测量噪声协方差矩阵，
Figure FDA0002533164480000049
为系统k-1时刻处于n模态，y1:k-1为从时刻1到时刻k-1的测量值序列，j为矩阵对角线元素索引，dy为yk的维度，
Figure FDA00025331644800000410
为k-1时刻测量噪声协方差矩阵第j个对角线元素，diag(·)为对角矩阵。
8.如权利要求2所述的三自由度直升机系统传感器运行状态监测方法，其特征在于：步骤(3)中，使用状态预测模型对每一模态下系统状态进行预测，所述状态预测模型为，
Figure FDA00025331644800000411
Figure FDA00025331644800000412
为系统状态粒子预测值；
Figure FDA00025331644800000413
为根据过程噪声的分布产生的噪声采样粒子；
Figure FDA00025331644800000414
为s模态下的融合初始粒子；
Figure FDA00025331644800000415
表示系统k时刻的模态为s。
9.如权利要求2所述的三自由度直升机系统传感器运行状态监测方法，其特征在于：步骤(3)中还包括计算每一模态下的粒子权重
Figure FDA0002533164480000051
计算公式为：
Figure FDA0002533164480000052
Figure FDA0002533164480000053
为k时刻s模态下粒子权重的预测值，exp(·)为以自然常数e为底的指数函数，(·)T为矩阵的转置，
Figure FDA0002533164480000054
为k时刻测量噪声协方差矩阵逆的期望。
10.如权利要求9所述的三自由度直升机系统传感器运行状态监测方法，其特征在于：步骤(3)中还包括，判断是否满足l＝L，是则输出k时刻系统各模态下的粒子权重值和逆伽马分布的尺度参数更新值，否则l＝l+1，并计算测量噪声协方差矩阵逆的期望
Figure FDA0002533164480000055
Figure FDA0002533164480000056
Figure FDA0002533164480000057
Figure FDA0002533164480000058
其中，
Figure FDA0002533164480000059
为k时刻s模态下粒子权重，
Figure FDA00025331644800000510
为第L步迭代时粒子权重值，
Figure FDA00025331644800000511
为第L步迭代时逆伽马分布的尺度参数，βk,j为逆伽马分布的尺度参数；
Figure FDA00025331644800000512
为第l-1步迭代时逆伽马分布的形状参数，
Figure FDA00025331644800000513
为第l-1步迭代时逆伽马分布的尺度参数，且
Figure FDA00025331644800000514
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 true CN111609878A (zh) 2020-09-01
CN111609878B 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)

Cited By (1)

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

Citations (6)

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

Patent Citations (6)

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

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
李洪志: "《信息融合技术》", 30 June 1996, 国防工业出版社 *

Cited By (1)

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

Also Published As

Publication number Publication date
CN111609878B (zh) 2021-06-22

Similar Documents

Publication Publication Date Title
US20190272354A1 (en) Method for degradation modeling and lifetime prediction considering recoverable shock damages
JP5657209B2 (ja) ノベルティ検出
CN112247992B (zh) 一种机器人前馈力矩补偿方法
CN106125548A (zh) 工业机器人动力学模型参数辨识方法
Low et al. On the long-term fatigue assessment of mooring and riser systems
CN111143981B (zh) 虚拟试验模型验证系统及方法
CN112949026B (zh) 一种考虑年龄和状态依赖的退化设备剩余寿命预测方法
CN111783238A (zh) 涡轮轴结构可靠性分析方法、分析装置及可读存储介质
CN113221263A (zh) 一种考虑分布参数不确定性的机械产品结构失效优化方法
CN106909738B (zh) 一种模型参数辨识方法
CN111609878B (zh) 三自由度直升机系统传感器运行状态监测方法
CN115979310A (zh) 一种惯导系统性能退化评估方法、系统、电子设备及介质
CN115270239A (zh) 基于动力特性和智能算法响应面法的桥梁可靠性预测方法
US7194394B2 (en) Method and apparatus for detecting and correcting inaccuracies in curve-fitted models
CN111339718A (zh) 一种管件冲蚀速率计算方法、装置及服务器
JP2020086778A (ja) 機械学習モデル構築装置および機械学習モデル構築方法
CN111461329B (zh) 一种模型的训练方法、装置、设备及可读存储介质
Solari The role of analytical methods for evaluating the wind-induced response of structures
WO2009120190A1 (en) Stable equilibrium point (sep) calculation apparatus of power system
CN116720600A (zh) 落货量预测方法、装置、设备及存储介质
CN113408057B (zh) 一种基于em算法的飞机气流角高精度离线估计方法
Khodaparast et al. Fuzzy model updating and its application to the DLR AIRMOD test structure
CN111581905B (zh) 隧道二极管电路系统在未知测量噪声下的状态估计方法
CN110321650B (zh) 基于新型试验设计与权重响应面的结构可靠性分析方法
JPH05209805A (ja) ばね・質点系のパラメータ同定装置およびその方法

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