一种FDI攻击下非线性信息物理系统滤波方法及系统

一种fdi攻击下非线性信息物理系统滤波方法及系统
技术领域
1.本发明涉及信息物理系统安全技术领域，具体为一种fdi攻击下非线性信息物理系统滤波方法及系统。


背景技术：

2.信息物理系统(cyber physical systems,简称cps)作为计算进程和物理进程的统一体，是集成计算、通信与控制于一体的下一代智能系统，信息物理系统通过人机交互接口实现和物理进程的交互，使用网络化空间以远程的、可靠的、实时的、安全的、协作的方式操控一个物理实体，随着网络技术和计算机技术的发展，信息物理系统在智慧城市、智能电网、和工业生产等领域均得到了广泛的应用，其中工业信息物理系统的安全运行稳定事关国家经济安全和人身安全。采用滤波方法对系统进行状态估计是工业信息物理系统的热点研究问题之一。卡尔曼滤波是一种常见的工业信息物理滤波方法，它是线性系统和高斯分布假设下的最优滤波器。然而，实际系统很难同时满足线性和高斯的假设。扩展卡尔曼滤波器通过对非线性系统的直接线性化操作，可实现对非线性高斯系统的状态估计。但扩展卡尔曼滤波不适用于高度非线性化的系统，而改进的高阶扩展卡尔曼滤波器会使计算复杂性增加。卡尔曼滤波方法假设状态的先验密度服从高斯分布，且系统过程噪声和测量噪声的协方差矩阵均为已知且服从高斯分布。因此，基于卡尔曼滤波的方法适用于高斯分布假设下的工业信息物理系统。
3.随着信息化和工业化深度融合进程的推进，工业信息物理系统大多采用通用协议与公共网络连接。由于协议漏洞等的存在，系统可能遭受网络攻击的威胁。虚假数据注入(false data injected,fdi)攻击是工业信息物理系统中的常见攻击之一，攻击者可通过向系统的执行器信号中注入虚假数据实施攻击，扰乱系统的正常运行且难以被发现。fdi攻击主要影响系统的状态估计。当工业信息物理系统遭受fdi攻击时，系统可能变成非高斯且不再稳定，将对工业生产运行造成重大隐患。此时采用卡尔曼滤波方法可能得到发散的滤波结果，因此无法直接用卡尔曼滤波方法很好实现fdi攻击下非线性非高斯工业信息物理系统的状态估计。
4.由于遭受fdi攻击的工业信息物理系统大多表现为非线性非高斯形式，目前针对工业信息物理系统的状态估计问题的大多数方法适用于线性高斯系统，且存在精度较低、计算量较高等问题，这些方法难以精确和高效地实现fdi攻击下非线性非高斯工业信息物理系统的状态估计。


技术实现要素：

5.(一)解决的技术问题
6.针对现有技术的不足，本发明提供了一种fdi攻击下非线性信息物理系统方法及系统，解决了上述背景技术中提出的目前针对信息物理系统的状态估计问题滤波的大多数方法适用于线性高斯系统，且存在精度较低、计算量较高等问题，这些方法难以精确和高效
地实现fdi攻击下非线性非高斯信息物理系统的状态估计的问题。
7.(二)技术方案
8.为实现以上目的，本发明通过以下技术方案予以实现：
9.一种fdi攻击下非线性信息物理系统状态估计滤波方法，包括：
10.步骤1.建立fdi攻击下非线性信息物理系统的状态空间模型；
11.步骤2.设置初始时刻k＝1，并进行模型参数初始化；
12.步骤3.在k+1时刻输入k时刻状态滤波值和滤波误差协方差矩阵，运行无迹高斯和滤波方法，实现状态滤波值的预测；
13.步骤4.计算未知参数的预测值
14.步骤5.k+1时刻的滤波值和滤波误差协方差矩阵更新；
15.步骤6.更新k+1时刻的未知参数估计值
16.步骤7.判断时刻k是否小于最后时刻n，若k＜n，则设置k＝k+1，输出k+1时刻的状态估计值返回步骤3，否则结束计算。
17.优选地，所述建立fdi攻击下非线性信息物理系统的状态空间模型，包括：
18.1.1)一个离散时间非线性信息物理系统模型的表达形式为：
[0019][0020]
其中xk∈rm表示系统在k时刻的状态向量，r表示实数域，m表示系统状态的维数，zk∈rn表示系统在k时刻的测量向量，n表示测量向量的维数，uk∈rm表示系统在k时刻的控制输入，f(xk)表示非线性状态转移函数，h(xk)表示非线性测量函数，nk表示未知的非高斯过程噪声，vk表示k时刻的测量噪声，是协方差矩阵为lk的零均值高斯白噪声，协方差矩阵lk是n
×
n维对称方阵；
[0021]
1.2)考虑系统(1)会遭受执行器fdi攻击，引入执行器fdi攻击信号α
ak
，将系统模型(1)写成以下模型形式：
[0022][0023]
其中α
ak
表示执行器fdi攻击信号，服从非高斯分布；
[0024]
1.3)重新定义系统的过程噪声，将式(2)写成以下模型形式：
[0025][0026]
其中表示未知的总和非高斯过程噪声；
[0027]
1.4)将总和非高斯噪声的概率密度函数近似表示为以下高斯混合分布：
[0028][0029]
其中n(.|μ,σ)表示均值为μ、方差为σ的高斯概率密度函数；nn表示高斯混合的
数量，表示每个高斯分布占全部分布的权重。
[0030]
优选地，所述设置初始时刻k＝1，并进行模型参数初始化，包括：
[0031][0032][0033]
其中，表示初始时刻状态滤波期望的均值，x0表示初始时刻系统状态；p
0|0
表示初始时刻滤波误差期望协方差矩阵；e(.)表示数学期望函数；()
t
表示转置运算。
[0034]
优选地，所述在k+1时刻输入k时刻状态滤波值和滤波误差协方差矩阵，运行无迹高斯和滤波方法，实现状态滤波值的预测；包括：
[0035]
在k+1时刻输入k时刻状态滤波值和滤波误差协方差矩阵p
k|k
，运行无迹高斯和滤波方法，实现状态滤波值的预测；
[0036]
3.1)生成2m+1个状态转移函数的sigma点集包括：
[0037]
令i＝0，重复执行2m+1次以下操作：
[0038]
若i＝0，得到若1≤i≤m，得到
[0039]
若m+1≤i≤2m，得到
[0040]
其中表示k时刻预测k+1时刻的状态向量，表示k时刻的状态滤波值，m表示系统状态的维数；参数λ＝α2(m+κ)-m；参数α表示缩放因子；参数κ表示缩放参数；
[0041]
3.2)生成2m+1个状态转移函数sigma点集对应权重包括：令i＝0，重复执行2m+1次以下操作：
[0042]
若i＝0，得到
[0043]
若1≤i≤2m，得到
[0044]
i＝i+1；
[0045]
3.3)重定义状态转移函数的sigma点如下：
[0046][0047][0048]
其中，表示第i个状态转移函数的sigma点；表示所有的线性加权和；
[0049]
3.4)将状态转移函数的sigma点代入式zk＝h(xk)+vk，得到：
[0050][0051]
其中，表示第i个测量函数的sigma点，h(xk)表示非线性测量函数；
[0052]
3.5)通过式(10)计算的均值
[0053][0054]
其中，表示的均值；
[0055]
3.6)通过式(11)计算互协方差矩阵
[0056][0057]
其中表示互协方差矩阵，用于公共卡尔曼增益矩阵的计算；
[0058]
3.7)通过式(12)计算一步预测协方差矩阵
[0059][0060]
其中表示一步预测协方差矩阵，用于公共卡尔曼增益矩阵的计算，lk为系统测量噪声的协方差矩阵；
[0061]
3.8)令系统(1)在k时刻状态的一步预测分布p(x
k+1
|z
1:k
)为如下高斯混合模型：
[0062][0063]
其中为未知参数，m为系统状态的维数。
[0064]
优选地，所述计算未知参数的预测值包括：
[0065]
4.1)设计未知参数p
k+1|k
的条件先验分布p(p
k+1|k
|z
1:k
)为如下形式：
[0066][0067]
其中，未知参数用于公共增益矩阵的计算；为n
×
n维矩阵a的逆威沙特分布概率密度函数；e(a)＝t(t-n-1)，e(.)表示期望函数；exp{}表示指数函数；t表示自由度；γn(.)表示伽马函数；tr(.)表示迹函数；
[0068]
4.2)令k+1时刻第i个逆威沙特分布中未知参数的预测值和如下：
[0069][0070]
其中，表示式(14)中逆威沙特分布的参数之一，ρ表示遗忘因子；
[0071][0072]
其中，表示式(14)中逆威沙特分布的参数之一；
[0073]
4.3)计算未知参数的估计值
[0074][0075]
其中可用于状态滤波值更新步骤的计算。
[0076]
优选地，所述k+1时刻的滤波值和滤波误差协方差矩阵更新，包括：
[0077]
5.1)通过公式(18)计算公共卡尔曼增益矩阵，通过公式(19)计算k+1时刻第i个滤波器的滤波值，通过公式(20)计算k+1时刻第i个滤波器的滤波误差协方差矩阵，通过公式(21)计算k+1时刻第i个滤波器的权重，通过公式(22)计算k+1时刻经过归一化处理后的第i个滤波器权重；
[0078][0079][0080][0081][0082][0083]
其中，表示k+1时刻公共卡尔曼增益矩阵；表示k+1时刻第i个滤波器的滤波值；表示k+1时刻第i个滤波器的滤波误差协方差矩阵；表示k+1时刻第i个滤波器的权重；表示k+1时刻经过归一化处理后的第i个滤波器权重；
[0084]
5.2)通过公式(23)计算k+1时刻滤波值通过公式(24)计算p
k+1|k+1
滤波误差协方差矩阵；
[0085][0086][0087]
其中，表示k+1时刻滤波值，是2m+1个滤波器滤波值的线性加权和；p
k+1|k+1
表示k+1时刻的滤波误差协方差矩阵，是2m+1个滤波器的滤波误差协方差矩阵线性加权和。
[0088]
优选地，所述更新k+1时刻的未知参数估计值包括：
[0089]
6.1)令系统总和过程噪声的后验分布由如下混合高斯分布近似表示：
[0090][0091]
其中，m为系统状态的维数，参数的后验概率密度函数为的后验概率密度函数为表示参数的估计值；
[0092]
6.2)采用变分贝叶斯方法近似后验联合概率分布即：
[0093][0094][0095]
其中，
[0096][0097]
6.3)将系统状态x
k+1
和未知参数估计值的联合后验概率密度函数用以下式子近似表示：
[0098]
p(x
k+1
,p
k+1|k
|z
1:k+1
)≈q(x
k+1
)q(p
k+1|k
)，
ꢀꢀꢀ
(27)
[0099]
其中z
1:k+1
为系统从1时刻到k+1时刻的测量输出；
[0100]
6.4)定义真实后验分布p(x
k+1
,p
k+1|k
|z
1:k+1
)和近似分布q(x
k+1
)q(p
k+1|k
)之间的kl散度如下：
[0101][0102]
其中∫表示积分操作，log()表示取对数；
[0103]
6.5)通过最小化式(28)，可得如下式子：
[0104]
logq(θ)＝e
ξ-θ
[logp(ξ,z
1:k+1
)]+c
θ
，
ꢀꢀꢀ
(29)
[0105][0106]
其中θ表示集合ξ中的任意元素；c
θ
表示与θ相关的常数；
[0107]
6.6)由式(29)计算式(26)的对数期望值，可得：
[0108][0109]
其中，表示与参数p
k+1|k
相关的常数；tr(.)表示迹函数；式(31)中如下：
[0110][0111]
6.7)将式(32)表示为如下形式：
[0112][0113][0114]
由式(31)可知，未知参数估计值p
k+1|k
概率密度函数的近似值q(p
k+1|k
)为2m+1个逆
威沙特分布的概率密度函数的乘积，即
[0115][0116][0117][0118]
其中，q(p
k+1|k
)表示未知参数估计值p
k+1|k
概率密度函数的近似值，q(p
k+1|k
)是2m+1个逆威沙特分布的乘积；表示第i个逆威沙特分布的后验参数；表示第i个逆威沙特分布的后验参数；
[0119]
6.8)由式(35)、式(36)和逆威沙特分布的性质得到：
[0120][0121]
其中表示k+1时刻更新后的参数。
[0122]
本发明还提供一种fdi攻击下非线性信息物理系统状态估计系统，包括：
[0123]
状态空间建模模块：用于建立fdi攻击下非线性信息物理系统的状态空间模型；
[0124]
初始化设置模块：用于设置初始时刻k＝1，并进行模型参数初始化；
[0125]
状态滤波值预测模块：用于在k+1时刻输入k时刻状态滤波值和滤波误差协方差矩阵，运行无迹高斯和滤波方法，实现状态滤波值的预测；
[0126]
状态估计模块：用于计算未知参数的预测值
[0127]
k+1时刻的滤波值和滤波误差协方差矩阵更新；
[0128]
更新k+1时刻的未知参数估计值
[0129]
判断时刻k是否小于最后时刻n，若k＜n，则设置k＝k+1，输出k+1时刻的状态估计值，否则结束计算。
[0130]
本发明还提供一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，所述计算机程序包括程序指令，所述程序指令当被处理器所执行时执行如前任一所述的一种fdi攻击下非线性信息物理系统滤波方法。
[0131]
本发明还提供一种fdi攻击下非线性信息物理系统滤波终端，包括：输入设备、输出设备、存储器、处理器；所述输入设备、所述输出设备、所述存储器和所述处理器相互连接，其中，所述存储器用于存储计算机程序，所述计算机程序包括程序指令，所护处理器被配置调用所述程序指令，执行如前任一所述的一种fdi攻击下非线性信息物理系统滤波方法。
[0132]
(三)有益效果
[0133]
本发明提供了一种fdi攻击下非线性信息物理系统滤波方法及系统。具备以下有益效果：
[0134]
本发明技术方案为实现fdi攻击下非线性非高斯信息物理系统的实时状态估计，引入无迹高斯和滤波方法和偏变分贝叶斯方法，保证状态估计精度的同时且有效减少方法的计算复杂度。
[0135]
本发明将fdi攻击建模为非高斯信号，无需假定fdi攻击服从某种特殊的分布(如伯努利分布、正态分布等)，将fdi攻击与系统的未知过程噪声描述为系统的总和非高斯过程噪声，进而将fdi攻击下非线性信息物理系统建模为具有未知非高斯过程噪声的非线性系统，可以将fdi攻击下非线性信息物理系统的滤波问题转化为具有未知过程噪声的非线性非高斯系统的滤波问题。
[0136]
为减少计算复杂度，采用无迹高斯和滤波方法进行滤波，可有效减少方法的计算复杂度，并将无迹高斯和滤波方法推广至具有未知非高斯过程噪声的非线性系统，采用偏变分贝叶斯方法，可实现未知参数的推断。
附图说明
[0137]
图1为本发明提供的一种fdi攻击下非线性信息物理系统滤波方法流程图；
[0138]
图2为本发明提供的一种fdi攻击下非线性信息物理系统滤波系统结构图。
具体实施方式
[0139]
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述。
[0140]
本发明实施例提供一种fdi攻击下非线性信息物理系统滤波方法，如图1所示，包括：
[0141]
s1.建立fdi攻击下非线性信息物理系统的状态空间模型；
[0142]
s2.设置初始时刻k＝1，并进行模型参数初始化；
[0143]
s3.在k+1时刻输入k时刻状态滤波值和滤波误差协方差矩阵，运行无迹高斯和滤波方法，实现状态滤波值的预测；
[0144]
s4.计算未知参数的预测值
[0145]
s5.k+1时刻的滤波值和滤波误差协方差矩阵更新；
[0146]
s6.更新k+1时刻的未知参数估计值
[0147]
s7.判断时刻k是否小于最后时刻n，若k＜n，则设置k＝k+1，输出k+1时刻的状态估计值返回步骤3，否则结束计算。
[0148]
在一个实施例中，所述建立fdi攻击下非线性信息物理系统的状态空间模型，包括：
[0149]
1.1)一个离散时间非线性信息物理系统模型的表达形式为：
[0150][0151]
其中xk∈rm表示系统在k时刻的状态向量，r表示实数域，m表示系统状态的维数，zk∈rn表示系统在k时刻的测量向量，n表示测量向量的维数，uk∈rm表示系统在k时刻的控制输入，f(xk)表示非线性状态转移函数，h(xk)表示非线性测量函数，nk表示未知的非高斯过程噪声，vk表示k时刻的测量噪声，是协方差矩阵为lk的零均值高斯白噪声，协方差矩阵lk是n
×
n维对称方阵；
[0152]
1.2)考虑系统(1)会遭受执行器fdi攻击，引入执行器fdi攻击信号α
ak
，将系统模型(1)写成以下模型形式：
[0153][0154]
其中α
ak
表示执行器fdi攻击信号，服从非高斯分布；
[0155]
1.3)重新定义系统的过程噪声，将式(2)写成以下模型形式：
[0156][0157]
其中表示未知的总和非高斯过程噪声；
[0158]
1.4)将总和非高斯噪声的概率密度函数近似表示为以下高斯混合分布：
[0159][0160]
其中n(.|μ,σ)表示均值为μ、方差为σ的高斯概率密度函数；nn表示高斯混合的数量，表示每个高斯分布占全部分布的权重。
[0161]
在一个实施例中，所述设置初始时刻k＝1，并进行模型参数初始化，包括：
[0162][0163][0164]
其中，表示初始时刻状态滤波期望的均值，x0表示初始时刻系统状态；p
0|0
表示初始时刻滤波误差期望协方差矩阵；e(.)表示数学期望函数；()
t
表示转置运算。
[0165]
在一个实施例中，所述在k+1时刻输入k时刻状态滤波值和滤波误差协方差矩阵，运行无迹高斯和滤波方法，实现状态滤波值的预测；包括：
[0166]
在k+1时刻输入k时刻状态滤波值和滤波误差协方差矩阵p
k|k
，运行无迹高斯和滤波方法，实现状态滤波值的预测；
[0167]
3.1)生成2m+1个状态转移函数的sigma点集包括：
[0168]
令i＝0，重复执行2m+1次以下操作：
[0169]
若i＝0，得到若1≤i≤m，得到
[0170]
若m+1≤i≤2m，得到
[0171]
其中表示k时刻预测k+1时刻的状态向量，表示k时刻的状态滤波值，m表示系统状态的维数；参数λ＝α2(m+κ)-m；参数α表示缩放因子；参数κ表示缩放参数；
[0172]
3.2)生成2m+1个状态转移函数sigma点集对应权重包括：令i＝0，重复执行2m+1次以下操作：
[0173]
若i＝0，得到
[0174]
若1≤i≤2m，得到
[0175]
i＝i+1；
[0176]
3.3)重定义状态转移函数的sigma点如下：
[0177][0178][0179]
其中，表示第i个状态转移函数的sigma点；表示所有的线性加权和；
[0180]
3.4)将状态转移函数的sigma点代入式zk＝h(xk)+vk，得到：
[0181][0182]
其中，表示第i个测量函数的sigma点，h(xk)表示非线性测量函数；
[0183]
3.5)通过式(10)计算的均值
[0184][0185]
其中，表示的均值；
[0186]
3.6)通过式(11)计算互协方差矩阵
[0187][0188]
其中表示互协方差矩阵，用于公共卡尔曼增益矩阵的计算；
[0189]
3.7)通过式(12)计算一步预测协方差矩阵
[0190][0191]
其中表示一步预测协方差矩阵，用于公共卡尔曼增益矩阵的计算，lk为系统测量噪声的协方差矩阵；
[0192]
3.8)令系统(1)在k时刻状态的一步预测分布p(x
k+1
|z
1:k
)为如下高斯混合模型：
[0193][0194]
其中为未知参数，m为系统状态的维数。
[0195]
在一个实施例中，所述计算未知参数的预测值包括：
[0196]
4.1)设计未知参数p
k+1|k
的条件先验分布p(p
k+1|k
|z
1:k
)为如下形式：
[0197][0198]
其中，未知参数用于公共增益矩阵的计算；为n
×
n维矩阵a的逆威沙特分布概率密度函数；e(a)＝t(t-n-1)，e(.)表示期望函数；exp{}表示指数函数；t表示自由度；γn(.)表示伽马函数；tr(.)表示迹函数；
[0199]
4.2)令k+1时刻第i个逆威沙特分布中未知参数的预测值和如下：
[0200][0201]
其中，表示式(14)中逆威沙特分布的参数之一，ρ表示遗忘因子；
[0202][0203]
其中，表示式(14)中逆威沙特分布的参数之一；
[0204]
4.3)计算未知参数的估计值
[0205][0206]
其中可用于状态滤波值更新步骤的计算。
[0207]
在一个实施例中，所述k+1时刻的滤波值和滤波误差协方差矩阵更新，包括：
[0208]
5.1)通过公式(18)计算公共卡尔曼增益矩阵，通过公式(19)计算k+1时刻第i个滤波器的滤波值，通过公式(20)计算k+1时刻第i个滤波器的滤波误差协方差矩阵，通过公式(21)计算k+1时刻第i个滤波器的权重，通过公式(22)计算k+1时刻经过归一化处理后的第i个滤波器权重；
[0209][0210][0211][0212][0213][0214]
其中，表示k+1时刻公共卡尔曼增益矩阵；表示k+1时刻第i个滤波器的滤波值；表示k+1时刻第i个滤波器的滤波误差协方差矩阵；表示k+1时刻第i个滤波器的权重；表示k+1时刻经过归一化处理后的第i个滤波器权重；
[0215]
5.2)通过公式(23)计算k+1时刻滤波值通过公式(24)计算p
k+1|k+1
滤波误差协方差矩阵；
[0216][0217][0218]
其中，表示k+1时刻滤波值，是2m+1个滤波器滤波值的线性加权和；p
k+1|k+1
表示k+1时刻的滤波误差协方差矩阵，是2m+1个滤波器的滤波误差协方差矩阵线性加权和。
[0219]
在一个实施例中，所述更新k+1时刻的未知参数估计值包括：
[0220]
6.1)令系统总和过程噪声的后验分布由如下混合高斯分布近似表示：
[0221][0222]
其中，m为系统状态的维数，参数的后验概率密度函数为的后验概率密度函数为表示参数的估计值；
[0223]
6.2)采用变分贝叶斯方法近似后验联合概率分布即：
[0224][0225]
其中，
[0226][0227]
6.3)将系统状态x
k+1
和未知参数估计值的联合后验概率密度函数用以下式子近似表示：
[0228]
p(x
k+1
,p
k+1|k
|z
1:k+1
)≈q(x
k+1
)q(p
k+1|k
)，
ꢀꢀꢀ
(27)
[0229]
其中z
1:k+1
为系统从1时刻到k+1时刻的测量输出；
[0230]
6.4)定义真实后验分布p(x
k+1
,p
k+1|k
|z
1:k+1
)和近似分布q(x
k+1
)q(p
k+1|k
)之间的kl散度如下：
[0231][0232]
其中∫表示积分操作，log()表示取对数；
[0233]
6.5)通过最小化式(28)，可得如下式子：
[0234]
logq(θ)＝e
ξ-θ
[logp(ξ,z
1:k+1
)]+c
θ
，
ꢀꢀꢀ
(29)
[0235]
[0236]
其中θ表示集合ξ中的任意元素；c
θ
表示与θ相关的常数；
[0237]
6.6)由式(29)计算式(26)的对数期望值，可得：
[0238][0239]
其中，c
pk+1|k
表示与参数p
k+1|k
相关的常数；tr(.)表示迹函数；式(31)中如下：
[0240][0241]
6.7)将式(32)表示为如下形式：
[0242][0243]
由式(31)可知，未知参数估计值p
k+1|k
概率密度函数的近似值q(p
k+1|k
)为2m+1个逆威沙特分布的概率密度函数的乘积，即
[0244][0245][0246][0247]
其中，q(p
k+1|k
)表示未知参数估计值p
k+1|k
概率密度函数的近似值，q(p
k+1|k
)是2m+1个逆威沙特分布的乘积；表示第i个逆威沙特分布的后验参数；表示第i个逆威沙特分布的后验参数；
[0248]
6.8)由式(35)、式(36)和逆威沙特分布的性质得到：
[0249][0250]
其中表示k+1时刻更新后的参数。
[0251]
如图2所示，本发明实施例还提供一种fdi攻击下非线性信息物理系统滤波系统，包括：
[0252]
状态空间建模模块：用于建立fdi攻击下非线性信息物理系统的状态空间模型；
[0253]
初始化设置模块：用于设置初始时刻k＝1，并进行模型参数初始化；
[0254]
状态滤波值预测模块：用于在k+1时刻输入k时刻状态滤波值和滤波误差协方差矩阵，运行无迹高斯和滤波方法，实现状态滤波值的预测；
[0255]
状态估计模块：用于计算未知参数的预测值
[0256]
k+1时刻的滤波值和滤波误差协方差矩阵更新；
[0257]
更新k+1时刻的未知参数估计值
[0258]
判断时刻k是否小于最后时刻n，若k＜n，则设置k＝k+1，输出k+1时刻的状态估计值，否则结束计算。
[0259]
本发明实施例还提供一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，所述计算机程序包括程序指令，所述程序指令当被处理器所执行时执行如前任一所述的一种fdi攻击下非线性信息物理系统滤波方法。
[0260]
本发明实施例还提供一种fdi攻击下非线性信息物理系统滤波终端，包括：输入设备、输出设备、存储器、处理器；所述输入设备、所述输出设备、所述存储器和所述处理器相互连接，其中，所述存储器用于存储计算机程序，所述计算机程序包括程序指令，所护处理器被配置调用所述程序指令，执行如前任一所述的一种fdi攻击下非线性信息物理系统滤波方法。
[0261]
尽管已经示出和描述了本发明的实施例，对于本领域的普通技术人员而言，可以理解在不脱离本发明的原理和精神的情况下可以对这些实施例进行多种变化、修改、替换和变型，本发明的范围由所附权利要求及其等同物限定。