自适应弱敏秩卡尔曼滤波方法及其应用与流程

[0001]
本发明非线性系统状态估计技术领域，具体涉及一种自适应弱敏秩卡尔曼滤波方法及其应用。


背景技术：

[0002]
非线性状态估计问题广泛应用于航空航天器导航、运动目标跟踪、电力系统等众多领域。学者们提出了多种滤波方法来解决非线性系统的状态估计问题。秩卡尔曼滤波(rkf)适用于高斯分布、多元t分布、多元极值分布等非线性滤波，是建立在秩采样方法上的一种滤波方法。但是，rkf方法对系统模型的参数高度敏感，只有当系统模型的参数是精确已知时才能得到状态最优估计，当系统模型参数不确定时，状态估计精度下降。
[0003]
弱敏秩卡尔曼滤波(全称，drkf)是在rkf的基础上，将基于状态估计误差敏感性和敏感性权重加权的惩罚函数引入rkf的代价函数中，建立弱敏代价函数，并通过将该函数最小化获得弱敏最优增益，从而在一定程度上解决了不确定参数带来的状态估计误差敏感性问题。


技术实现要素：

[0004]
本发明要解决的技术问题是提供一种自适应弱敏秩卡尔曼滤波(adrkf)方法，并将其应用于非线性系统的状态估计过程，旨在解决非线性系统状态估计中的滤波精度不够所导致的准确性差的技术问题。
[0005]
由于弱敏秩卡尔曼滤波方法中没有给出合适的敏感性权重矩阵的取值方式，本发明基于该领域长期的实践研究，并结合量测残差的正交原理，设计了敏感性权重的自适应因子，以解决敏感性权重的自适应问题，进而降低非线性系统模型中不确定参数的影响，提高滤波精度。
[0006]
为解决上述技术问题，本发明采用如下技术方案：
[0007]
设计一种基于自适应弱敏秩卡尔曼滤波的非线性系统状态估计方法，主要包括以下步骤：
[0008]
(一)建立非线性系统的状态方程及量测方程
[0009]
x
k
＝f(x
k-1
,c)+w
k-1
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(1)
[0010]
z
k
＝h(x
k
,c)+v
k
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(2)
[0011]
其中，x
k
和z
k
分别是系统的状态向量和量测向量，k指第k步，代表t
k
时刻，f(
·
)和h(
·
)为非线性函数向量，c为不确定参数向量；w
k
和v
k
是零均值高斯白噪声，方差分别为q
k
和r
k
，满足：
[0012]
[0013]
其中，δ
kj
为kroneckerδ函数，当k＝j时，δ
kj
＝1；当k≠j时，δ
kj
＝0；为w
k
的转置矩阵为v
k
的转置矩阵。
[0014]
(二)对步骤(一)中的状态方程和量测方程进行自适应弱敏秩卡尔曼滤波估计
[0015]
(a)初始化非线性系统状态方程、量测方程的状态和状态误差方差阵
[0016][0017][0018]
其中，p0和q
k
、r
k
不相关。
[0019]
(b)计算状态和量测的秩采样点以及协方差和量测方差
[0020]
设第k-1步的状态估计值为误差方差阵为计算第k步的秩采样点集：
[0021][0022]
式中，u
p1
与u
p2
为正态偏量，用中位秩计算p
i
＝(i+2.7)/5.4,i＝1,2，p1＝0.6852,p2＝0.8704,u
p1
＝0.4823,u
p2
＝1.1281；n表示状态向量x的维数，为平方根的第j列向量；
[0023]
时间更新：计算状态一步预测
[0024][0025]
其中，
[0026][0027]
计算一步预测误差的方差阵：
[0028][0029]
其中，上标
“-”
表示变量的先验估计；为协方差权重系数。
[0030]
(c)量测更新
[0031]
重新秩采样，得到采样点集：
[0032][0033]
量测采样点集：
[0034][0035]
状态估计：
[0036][0037]
估计误差的方差阵：
[0038][0039]
式(13)中，
[0040][0041][0042]
其中，p
xz,k
为状态和量测的协方差，p
zz,k
为量测方差。
[0043]
(d)秩采样点的敏感性传播
[0044]
第一步，计算k-1步秩采样点的敏感性：
[0045][0046]
式中，为第k-1步的敏感性，上标“+”表示变量的后验估计。
[0047]
更新秩采样点集：
[0048][0049]
第二步，计算先验状态估计和先验协方差矩阵的敏感性
[0050][0051][0052]
式中为第k步的敏感性的先验估计；
[0053]
第三步，计算重新秩采样点集的敏感性和量测秩采样点的敏感性
[0054][0055][0056]
计算量测的敏感性：
[0057][0058]
第四步，计算状态和量测协方差的敏感性和量测方差的敏感性：
[0059][0060][0061]
第五步，计算状态估计的敏感性和状态误差方差阵的敏感性：
[0062][0063][0064]
其中，
[0065][0066][0067]
其中，是一个斜对称矩阵，满足γ
t
＝-γ，ψ和θ均为非奇异矩阵，且满足
[0068]
(e)计算滤波增益矩阵
[0069][0070]
其中，l为不确定参数个数，λ
k
为渐消因子，w
i,k
对应第i个不确定参数c
i
的敏感性权重；w
i,k
通过下述方法求得：
[0071][0072]
其中，bc
i
和ac
i
分别对应c
i
的上下限，为c
i
的均值；λ
k
通过下述方法求得：
[0073]
令：
[0074][0075]
式中，
[0076][0077][0078]
式(32)中v
k
为残差矩阵，通过下式估计：
[0079][0080]
渐消因子取为：
[0081][0082]
敏感性代价函数：
[0083][0084]
(f)计算敏感性矩阵
[0085][0086]
(g)计算状态更新
[0087][0088]
以上(a)至(g)步循环迭代，得到系统的实时状态监测结果。
[0089]
所述非线性系统包括弹道再入目标跟踪系统、电机转速估计、电力系统控制等。
[0090]
(ⅰ)建立弹道再入目标跟踪系统的模型：
[0091]
[0092][0093]
其中，x(t)＝[x1(t) x2(t) x3(t)]
t
为状态变量，分别为位置x1(t)(m)、速度x2(t)(m/s)和弹道常数x3(t)，c为具有不确定性的弹道系数；m和h为观测雷达的位置坐标；v为量测的零均值gauss白噪声。
[0094]
(ⅱ)建立电机转速状态估计的模型：
[0095]
异步电机离散状态方程：
[0096][0097]
其中，dt对应于构建量测方程步骤的采样时间间隔，是t
k-1
时刻的第一定子电压控制输入，为t
k-1
时刻的第二定子电压控制输入，x＝[x1,x2,x3,x4,x5]
t
为状态向量，x1和x2是定子电流，x3和x4是转子磁链，x5是角速度；x
k
表示t
k
时刻的状态向量；j是转子惯性；p
n
是极对数；u1和u2是定子电压控制输入；c＝[c
1 c2]，是不确定参数向量，c1和c2分别是定子电阻和转子电阻；w是零均值高斯白噪声；l
s
、l
r
和l
m
分别为转子电感、定子电感和互感。其它模型参数为：
[0098][0099][0100]
式中，k对应于t
k
时刻的步数；u
n
是三相对称电源的额定电压；f是供电频率；dt对应于构建量测方程步骤的采样时间间隔；u
n
＝[u
n1
,u
n2
,u
n3
]
t
。
[0101]
将测得的定子电流和转子磁链角速度作为量测值，以此建立相应的量测模型，则相应的量测方程为：
[0102]
z
k
＝h(x
k
,c)+v
k
＝[x1,x3]
t
+v
k
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(44)
[0103]
其中，w
k
和v
k
是相互独立的零均值高斯白噪声序列，且w
k
和v
k
的方差分别为q
k
和r
k
，且满足
[0104][0105]
其中，δ
kj
为kronecker δ函数，当k＝j时，δ
kj
＝1；当k≠j时，δ
kj
＝0；
[0106]
与现有技术相比，本发明的主要有益技术效果在于：
[0107]
1.本发明基于量测残差正交原理，结合工程应用实践经验，创造性的设计了敏感性权重自适应因子，提供了弱敏秩卡尔曼滤波(drkf)中敏感性权重的自适应取值方法，降低了因不确定参数而造成的状态估计误差，从而提高滤波精度，改善滤波性能。
[0108]
2.本发明滤波方法可广泛应用于航空航天器导航、运动目标跟踪、电机转速状态估计、电力系统控制等非线性系统或工程中，通过提高或改善其系统状态估计中的滤波精度和滤波性能，进而改善非线性系统监测的准确性。
附图说明
[0109]
图1为本发明自适应弱敏秩卡尔曼滤波的原理示意图。
[0110]
图2为本发明试验例与drkf和rkf在弹道目标再入过程中对于目标的状态监测结果的均方根误差对比图。
[0111]
图3为本发明试验例与drkf和rkf在感应电机空载启动过程中对于感应电机的状态监测结果的均方根误差对比图。
具体实施方式
[0112]
下面结合实施例来说明本发明的具体实施方式，但以下实施例只是用来详细说明本发明，并不以任何方式限制本发明的范围。
[0113]
在以下实施例中所涉及的仪器设备如无特别说明，均为常规仪器设备；所涉及的试验方法或计算方法，如无特别说明，均为常规方法；各公式或方程中所涉及的字母、符合，如无特别说明，其代表本领域常规的物理或数学含义。
[0114]
实施例一、基于自适应弱敏秩卡尔曼滤波的非线性系统状态估计方法，主要步骤如下：
[0115]
步骤(一)，建立非线性系统的状态方程及量测方程：
[0116]
x
k
＝f(x
k-1
,c)+w
k-1
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(1)
[0117]
z
k
＝h(x
k
,c)+v
k
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(2)
[0118]
其中，x
k
和z
k
分别是系统的状态向量和量测向量，k指第k步，代表t
k
时刻，f(
·
)和h(
·
)为非线性函数向量，c为不确定参数向量；w
k
和v
k
是零均值高斯白噪声，方差分别为q
k
和r
k
，满足：
[0119][0120]
其中，δ
kj
为kronecker δ函数，当k＝j时，δ
kj
＝1；当k≠j时，δ
kj
＝0；为w
k
的转置
矩阵为v
k
的转置矩阵。
[0121]
步骤(二)，对步骤(一)中的状态方程和量测方程进行自适应弱敏秩卡尔曼滤波估计(参见图1)：
[0122]
1.对离散的状态方程、量测方程的状态和状态误差方差阵分别进行初始化
[0123][0124][0125]
其中，p0、q
k
和r
k
均不相关。
[0126]
2.计算状态和量测的秩采样点以及协方差和量测方差
[0127]
设第k-1步的状态估计值和误差方差阵分别为和则第k步的秩采样点集为：
[0128][0129]
式中，u
p1
与u
p2
为正态偏量，用中位秩计算p
i
＝(i+2.7)/5.4,i＝1,2，p1＝0.6852,p2＝0.8704,u
p1
＝0.4823,u
p2
＝1.1281；n表示状态向量x的维数，为平方根的第j列向量。
[0130]
时间更新：状态一步预测为：
[0131][0132]
式(7)中，
[0133][0134]
一步预测误差的方差阵：
[0135][0136]
其中，上标
“-”
表示变量的先验估计；为协方差权重系数。
[0137]
量测更新：重新秩采样，得到采样点集：
[0138][0139]
量测均值：
[0140][0141]
状态估计：
[0142][0143]
估计误差的方差阵：
[0144][0145]
式(13)中，
[0146][0147][0148]
其中，p
xz,k
为状态和量测的协方差，p
zz,k
为量测方差。
[0149]
3.秩采样点的敏感性传播
[0150]
第一步，计算k-1步秩采样点的敏感性：
[0151][0152]
式中，为第k-1步的敏感性的后验估计，上标“+”表示变量的后验估计；更新秩采样点集：
[0153][0154]
第二步，计算先验状态估计和先验协方差矩阵的敏感性
[0155][0156][0157]
式中为第k步的敏感性的先验估计；
[0158]
第三步，计算重新秩采样点集和预测量测秩采样点的敏感性
[0159][0160][0161]
计算量测均值的敏感性：
[0162][0163]
第四步，计算状态和量测协方差以及量测方差的敏感性：
[0164][0165][0166]
最后，计算状态估计和状态误差方差阵的敏感性：
[0167][0168][0169]
式中，
[0170][0171]
式(27)中，
[0172][0173]
其中，是一个斜对称矩阵，满足γ
t
＝-γ，ψ和θ均为非奇异矩阵，且满足
[0174]
4.计算自适应弱敏秩卡尔曼滤波的卡尔曼增益k
k
[0175][0176]
其中，l为不确定参数个数，λ
k
为渐消因子，w
i,k
对应第i个不确定参数c
i
的敏感性权重；w
i,k
通过下述方法求得：
[0177][0178]
其中，bc
i
和ac
i
分别对应c
i
的上、下限，为c
i
的均值；；λ
k
通过下述方法求得：
[0179]
令：
[0180][0181]
式(31)中，
[0182][0183][0184]
以上式中v
k
为残差矩阵，通过下式估计：
[0185][0186]
渐消因子可取为：
[0187][0188]
敏感性代价函数：
[0189][0190]
5.计算敏感性矩阵
[0191][0192]
6.第k步的状态估计
[0193][0194]
以上6步循环迭代，得到系统的实时状态监测结果。
[0195]
试验例1
[0196]
以一个弹道目标的再入问题为例，通过对目标的状态检测，对本发明的方法进行仿真验证。首先建立如下状态方程和量测方程：
[0197]
状态方程：
[0198][0199]
量测方程：
[0200][0201]
其中，x(t)＝[x1(t) x2(t) x3(t)]
t
为状态变量，分别为位置x1(t)(m)、速度x2(t)(m/s)和弹道常数x3(t)，c为具有不确定性的弹道系数，经验值为m＝105m和h＝105m为观测雷达的位置坐标；v为量测的零均值gauss白噪声，其方差为r＝104m2。
[0202]
假设参数c服从均匀分布状态的初始真值和初始估计值分别为
[0203]
[x1(0),x2(0),x3(0)]＝[3
×
105m,-2
×
104m/s,1
×
10-3
]
ꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(41)
[0204][0205]
方差阵的初始值为
[0206]
p0＝diag{1
×
106m2,4
×
106m2/s2,1
×
10-4
}
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(43)
[0207]
整个monte carlo仿真时间为60s，采样时间为0.1s，仿真次数为1000次；连续状态方程式采用四阶runge-kutta方法进行离散化；仿真中，imp.rkf算法和adrkf算法的敏感性权重均为w＝104。
[0208]
使用本发明方法与本领域常规方法drkf和rkf对此弹道目标再入问题中目标的位置，速度以及弹道参数进行实时状态监测结果对比如下：
[0209]
在弹道目标再入过程中对于目标的状态监测，本发明方法与drkf、rkf监测结果的均方根误差对比如图2所示。
[0210]
从图2中可以看出，per.rkf算法中系统模型的参数是精确已知的，因此取得了最好的估计效果；imp.rkf，drkf与adrkf算法中没有给出具体的不确定参数。可以看出，与imp.rkf和drkf算法相比，本发明方法的均方根误差值最小，具有更好的精确度。
[0211]
特别地，per.rkf算法为已知参数真实值的rkf(perfect rkf)，imp.rkf为仅知道参数经验值而不知真实值的rkf(imperfect rkf)，drkf为具有解析增益的弱敏rkf(desensitized rank kalman filter)。
[0212]
试验例2
[0213]
再以异步电机空载启动的状态监测为例，对本发明的方法进行仿真验证。
[0214]
建立异步电机的离散状态方程：
[0215]
[0216]
其中，dt对应于构建量测方程步骤的采样时间间隔，是t
k-1
时刻的第一定子电压控制输入，为t
k-1
时刻的第二定子电压控制输入，x＝[x1,x2,x3,x4,x5]
t
为状态向量，x1和x2是定子电流，x3和x4是转子磁链，x5是角速度；x
k
表示t
k
时刻的状态向量；j是转子惯性；p
n
是极对数；u1和u2是定子电压控制输入；c＝[c
1 c2]，是不确定参数向量，c1和c2分别是定子电阻和转子电阻；w是零均值高斯白噪声；其它模型参数为：
[0217][0218][0219]
其中，转子电感l
s
＝0.265[h]，定子电感l
r
＝0.265[h]，互感l
m
＝0.253[h]，转子惯性j＝0.02[kg
·
m2]，极对数p
n
＝2，k对应于t
k
时刻的步数；u
n
是三相对称电源的额定电压；f是供电频率；dt对应于构建量测方程步骤的采样时间间隔；u
n
＝[u
n1
,u
n2
,u
n3
]
t
。
[0220]
将测得的定子电流和转子磁链角速度作为量测值，以此建立相应的量测模型，则相应的量测方程为：
[0221]
z
k
＝h(x
k
,c)+v
k
＝[x1,x3]
t
+v
k
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(47)
[0222]
其中，w
k
和v
k
是相互独立的零均值高斯白噪声序列，且w
k
和v
k
的方差分别为q
k
和r
k
，且满足
[0223][0224]
其中，δ
kj
为kroneckerδ函数，当k＝j时，δ
kj
＝1；当k≠j时，δ
kj
＝0；
[0225]
经保密实验，发明人获取的系统噪声方差阵q
k
和量测噪声方差阵r
k
矩阵如下：
[0226]
其中，观测次数为n＝150，总采样时间为t＝0.15s
[0227]
初始状态初始误差方差矩阵p0、q
k
和r
k
均不相关；
[0228]
使用本发明方法与本领域常规方法drkf和rkf对感应电机的定子电流、转子磁链、角速度参数进行实时状态监测结果对比如下：
[0229]
在感应电机空载启动过程中对于感应电机的状态监测，本发明方法与drkf、rkf监测结果的均方根误差对比图如图3所示。
[0230]
从图3中可以看出，本发明的均方根误差值较小，具有更好的精确度。
[0231]
以上两个仿真试验例均是在matlab(r2016b)软件上进行的建模和仿真，并在cpu为i5-7400、内存为8g的电脑上进行运行。仿真过程中，在matlab(r2016b)软件上通过程序的编写搭建了仿真模型，并进行了初始数据的输入(上面的具体实施过程中已给出)，然后再通过运行matlab(r2016b)软件进行计算。图2与图3均是通过matlab(r2016b)仿真计算获取。
[0232]
上面结合实施例对本发明作了详细的说明；但是，所属技术领域的技术人员能够理解，在不脱离本发明构思的前提下，还可以对上述实施例中的各个具体参数进行变更，或者对相关技术手段进行等同替代，从而形成多个具体的实施例，均为本发明的常见变化范围，在此不再逐一详述。