CN109765786A - 一种基于证据滤波的电推船电机转轴不平衡故障检测方法 - Google Patents

Abstract

本发明涉及一种基于证据滤波的电推船电机转轴不平衡故障检测方法，属于机电设备状态监测与故障诊断领域。从推进电机传动轴采集的振动数据进行预处理，利用连续的Sigmoid隶属度函数将处理后的振动数据转化成每时刻的报警证据；利用多阶证据滤波方法对每时刻报警证据进行融合，其中构建信度推理模型求解多阶证据滤波中的融合权重，并给出优化目标函数对推理模型中的参数进行优化。对于在线获取监测数据，预处理后转化为每时刻报警证据，利用优化后的多阶证据滤波方法得到全局报警证据，进行报警决策；本发明通过多阶证据滤波得到全局的报警证据，可以给出比单个时刻的报警证据更为准确的故障检测结果。

Description

一种基于证据滤波的电推船电机转轴不平衡故障检测方法

技术领域

本发明涉及一种基于证据滤波的中压电推船舶推进电机传动轴不平衡故障检测方法，属于机电设备状态监测与故障诊断领域。

背景技术

电力推进作为一种区别于传统推进的新型推进方式，目前在军用和工程类船舶中已经得到较为普遍的应用，其能否安全稳定的运行也越来越受到人们的关注。因此，准确及时地对可能发生的故障进行报警，为主动维修争取时间，这对保证推进系统以及船舶的安全运行至关重要。由于故障检测与报警是一种可以降低损失、增加经济效益的技术，所以它已经被广泛的应用在各个领域行业中。但是对于船舶电力推进系统的推进电机而言，它的机动性强，电机输出功率变化频繁，国内外对其故障检测技术的研究较少，所以对推进电机故障检测技术的研究具有十分重要的现实意义和实用价值。

发明内容

本发明的目的是提出一种基于证据滤波的中压电推船舶推进电机传动轴不平衡故障检测方法，对于每个时刻获得的振动数据，将其用Sigmoid隶属度阈值函数转换为相应的报警证据，通过多阶证据滤波对每一时刻报警证据进行融合，得到可靠的全局报警证据，利用全局的报警证据进行决策。由于全局报警证据包含了当前和历史时刻电机的运行振动信息，利用其做出的故障报警决策比任何单个时段振动信息给出的决策结果更为可靠和准确。

本发明包括以下各步骤：

(1)设定中压电推船舶推进电机传动轴不平衡故障报警器的辨识框架为Θ＝{NA,A}，其中NA表示设备处于正常运行状态，A表示设备处于不平衡故障的异常运行状态亦即报警状态。

(2)设x为安装在推进电机传动轴上的振动传感器所检测的振动加速度变量，令x(t)，t＝1,2,3,…,是x在t时刻的采样值，单位为m/s²，利用式(1)对x(t)进行预处理，得到处理后的变量y及其采样值为：

定义y_otp为最优阈值，记max(y)和min(y)分别是y的最大值和最小值，其中y_otp∈[min(y),max(y)]，采用直接门限法最小化误报率和漏报率的平方和求得最优阈值y_otp；利用Sigmoid隶属度函数构建正常状态NA、异常状态A的模糊隶属度函数μ_NA(y(t))、μ_A(y(t))，如式(2)-(3)所示：

其中a∈[0.1,100]是一个可调参数。

(3)将y(t)分别带入步骤(2)关于正常状态NA、异常状态A的模糊隶属度函数μ_NA(y(t))、μ_A(y(t))中，即可得到报警证据m_t(NA)、m_t(A)，分别如下(4)-(5)式所示：

m_t(NA)＝μ_NA(y(t))/(μ_A(y(t))+μ_NA(y(t))) (4)

m_t(A)＝μ_A(y(t))/(μ_A(y(t))+μ_NA(y(t))) (5)

那么，在每个采样时刻，能够获得报警证据m_t＝(m_t(NA),m_t(A))，其中m_t(A)和m_t(NA)表示x(t)支持“报警A”和“正常NA”这两个事件发生的程度。

(4)在步骤(3)的基础上得到了每一时刻的报警证据之后，利用证据的动态更新得到全局报警证据，记为m_0:t＝(m_0:t(A),m_0:t(NA))，具体步骤如下：

(4-1)当t＝1时，因在t＝1之前没有相关的证据信息,所以t＝1时全局的报警证据是该时刻的报警证据即m_0:1＝(m₁(A),m₁(NA))。

(4-2)当t＝2时，利用t＝1时刻的全局报警证据m_0:1和t＝2时刻的证据m₂进行加权融合，得到t＝2时刻的全局报警证据m_0:2，具体用于计算的证据滤波公式如下

其中τ₂＝0.75表示对m_0:1(A)的线性加权值，ξ₂＝0.25表示当前证据的加权值；若m₂(NA)>m₂(A)，则D＝NA，那么m₂(A|NA)＝0；若m₂(A)≥m₂(NA)，则D＝A，那么m₂(A|A)＝1；

(4-3)当t＝3时利用t＝2时刻的全局报警证据m_0:2和t＝3时刻的证据m₃进行加权融合，得到t＝3时刻的全局报警证据m_0:3＝(m_0:3(A),m_0:3(NA))，具体用于计算的证据滤波公式如下

若m₃(NA)>m₃(A)，则D＝NA，那么m₃(A|NA)＝0；若m₃(A)≥m₃(NA)，则D＝A，那么m₃(A|A)＝1；τ₃表示对m_0:2(A)的线性加权值，ξ₃表示当前证据的加权值，τ₃、ξ₃的计算步骤如下：

(4-3-1)定义m_0:1、m_0:2两条证据之间的距离如式(8)所示：

其中，上标T表示求向量(m_0:1-m_0:2)的转置，那么m_0:1和m_0:2两条证据之间的相似度如式(9)所示：

Sim(m_0:1,m_0:2)＝1-d_J(m_0:1,m_0:2) (9)

Sim(m_0:1,m_0:2)越大，两者越相似，同理可以计算出m₃、m_0:1和m_0:2之间两两相似度分别为Sim(m_0:1,m₃)以及Sim(m_0:2,m₃)。

(4-3-2)根据步骤(4-3-1)得到的相似度求出每个证据被其他两个证据支持的程度，分别为Sup(m_0:2)、Sup(m_0:1)以及Sup(m₃)，如下式(10)-(12)所示

Sup(m_0:2)＝Sim(m_0:1,m_0:2)+Sim(m_0:2,m₃) (10)

Sup(m_0:1)＝Sim(m_0:1,m_0:2)+Sim(m_0:1,m₃) (11)

Sup(m₃)＝Sim(m_0:1,m₃)+Sim(m_0:2,m₃) (12)

(4-3-3)在步骤(4-3-2)的基础上得到证据动态更新的权重因子τ₃和ξ₃，如下式(13)和(14)所示

τ₃＝1-ξ₃ (14)

(4-4)当t≥4时，求取m_0:t＝(m_0:t(A),m_0:t(NA))，具体用于计算的证据滤波公式如下

其中τ_t+ξ_t＝1，若m_t(NA)>m_t(A)，则D＝NA，那么m_t(A|NA)＝0；若m_t(A)≥m_t(NA)，则D＝A，那么m_t(A|A)＝1；同理可得m_t-1(A|D)，m_t-2(A|D)。

(5)式(15)中和的具体确定步骤如下：

(5-1)当t≥4时，依据步骤(4-3)中的子步骤(4-3-1)提供的方法得到m_t、m_t-1和m_t-2两两之间的相似度Sim(m_t,m_t-1)、Sim(m_t,m_t-2)和Sim(m_t-1,m_t-2)；

(5-2)依据步骤(4-3)中的子步骤(4-3-2)提供的方法分别得到m_t、m_t-1和m_t-2相对于其他两个证据的支持程度Sup(m_t)、Sup(m_t-1)和Sup(m_t-2)；

(5-3)在步骤(5-2)的基础上，和可以由式(16)-(18)所示

(6)对于步骤(4)中式(15)的参数ξ_t，计算过程如下：

(6-1)依据步骤(4-3)中的子步骤(4-3-1)分别得到m_t、m_t-1、m_t-2和m_0:t-3两两之间的相似度分别为Sim(m_t,m_t-1)、Sim(m_t,m_t-2)、Sim(m_t,m_0:t-3)、Sim(m_t-1,m_0:t-3)、Sim(m_t-1,m_t-2)以及Sim(m_t-2,m_0:t-3)。

(6-2)依据步骤(4-3)中的子步骤(4-3-2)分别得到m_t、m_t-1、m_t-2和m_0:t-3对于其他两个证据的支持程度，分别为Sup(m_0:t-3)、Sup(m_t)、Sup(m_t-1)以及Sup(m_t-2)。

(6-3)构造关于ξ_t的信度推理模型，用于描述输入量Sup(m_0:t-3)、Sup(m_t-1)、Sup(m_t-2)和Sup(m_t)与输出量线性加权值ξ_t之间存在的复杂非线性关系，具体步骤如下：

(6-3-1)将输入变量Sup(m_0:t-3)、Sup(m_t)、Sup(m_t-1)以及Sup(m_t-2)记为f_i(t)(i＝1,2,3,4)；设定ξ_t的参考值集合D＝{D_n|n＝1,...,N}，输入变量f_i(t)的参考值集1<N<10为ξ_t的结果参考值个数，1<J_i<10为输入变量的参考值个数。

(6-3-2)给定如表1所示的信度矩阵表来描述输入f_i(t)和输出ξ_t之间的关系，由表1中可知，当输入值f_i(t)取参考值时，结果值ξ_t为参考值D_n的信度为 并有定义对应于参考值的信度为同时给定输入f_i(t)的可靠度为r_i，满足0≤r_i≤1。

表1 输入f_i(t)的信度矩阵表

(6-3-3)当获得输入样本f₁(t)、f₂(t)、f₃(t)和f₄(t)之后，根据信度矩阵表1和可靠度r₁、r₂、r₃和r₄，可利用证据推理规则推理出线性加权值的估计值具体步骤如下：

a)对于输入值f_i(t)，其必然落入某两个参考值构成的区间此时这两个参考值对应的信度和被激活，α_i,j表示输入值f_i(t)匹配参考值的相似度，α_i,j+1表示输入值f_i(t)匹配参考值的相似度，由式(19)所示

则输入值f_i(t)的信度可由和以加权和的形式获得

e_i＝{(D_n,p_n,i),n＝1,...,N} (20)

b)利用式(20)和式(21)获得f₁(t)、f₂(t)、f₃(t)和f₄(t)的信度e₁、e₂、e₃和e₄，给定e_i的可靠度r_i为1，设e_i的证据权重初值w_i＝r_i＝1，利用证据推理规则对e₁、e₂、e₃和e₄进行融合，融合过程如下：

(a)首先对e₁和e₂进行融合，得到融合后的证据为

(b)将e₁和e₂融合后的概率证据设为e′₂＝{(D_n,p′_n,2),n＝1,...,N}，将其与e₃进行融合，其中令e′₂的证据权重w′₂为1，可靠度r′₂为1，得到融合后的概率证据为：

(c)将e′₂和e₃融合后的概率证据设为e″₂＝{(D_n,p″_n,2),n＝1,...,N}，将其与e₄进行融合，其中令e″₂的证据权重w″₂为1，可靠度r″₂为1，作为最终的融合结果，即有：

O(f(t))＝{(D_n,p″_n,2),n＝1,...,N} (25)

(d)根据步骤(c)得到融合结果O(f(t))，线性加权值的估计值可由下式推理得到

式(12)中线性加权值τ_t的估计值为

(7)基于Jousselme证据距离构建对信度矩阵表1中参数优化模型，具体步骤如下：

(7-1)确定优化参数集合

(7-2)将最小化Jousselme证据距离作为优化目标函数

s.t.0≤w_i≤1,i＝1,2,3,4 (28)

D₂<D₃<…<D_N-1 (30)

其中L为采集的训练样本数量总数，m表示为设备所处的真实状态下的理想向量(010)或异常状态该下的理想向量(100)，式(28)-(31)表示优化参数需满足的约束条件。

(8)利用基于梯度下降法方法或者非线性优化软件包，获得最优的参数集合P，从而得到优化后的输入f_i(t)的信度矩阵表；从在线运行的推进电机传动轴上采集振动加速度信号，将其利用步骤(2)中式(1)进行处理，按照一定间隔选取阈值，再重复步骤(3)-(6)，可得到全局报警证据m_0:t＝(m_0:t(A),m_0:t(NA))；根据得到的t时刻的全局报警证据m_0:t＝(m_0:t(A),m_0:t(NA))，将m_0:t转化为相应的pignistic概率如下

其中，所以基于pignistic概率报警决策准则为：若则报警器发出报警，反之，不发出报警；得到在不同阈值下误报率和漏报率的平方和，误报率和漏报率的平方和最小所对应的阈值即为最优阈值。

本发明提出的基于证据滤波的中压电推船舶推进电机传动轴不平衡故障检测方法。该方法先将传感器采集的数据进行预处理，利用连续的Sigmoid隶属度函数将处理后的振动数据的信息转化成每时刻的报警证据；利用多阶证据滤波方法对每时刻报警证据进行融合，其中构建信度推理模型求解多阶证据滤波中的融合权重，并给出优化目标函数对推理模型中的参数进行优化，得到最优参数集合；在线采集数据，预处理后转化为每时刻报警证据，使用优化后的信度推理模型参数集合，按一定间隔选取不同的阈值利用多阶证据滤波公式得到全局报警证据，进行报警决策，误报率和漏报率的平方和最小所对应的阈值为最优阈值。根据本发明方法编制的程序(编译环境Matlab)可以在计算机上运行，并联合传感器、数据采集器等硬件组成船舶推进电机故障检测系统，从而实现对电力推进船舶推进电机故障检测并报警。

附图说明

图1是本发明方法的程序流程框图；

图2本发明方法的实施例中用来训练的x(t)序列；

图3本发明方法的实施例中用来训练的x(t)处理后得到的y(t)序列；

图4本发明方法的实施例中用来测试的在线采集的X(t)序列；

图5本发明方法的实施例中用来在线采集的X(t)序列处理后得到的Y(t)序列.

具体实施方式

本发明提出的一种基于证据滤波的中压电推船舶推进电机传动轴不平衡故障检测方法，其流程框图如图1所示，包括以下各步骤：

(1)设定中压电推船舶推进电机传动轴不平衡故障报警器的辨识框架为Θ＝{NA,A}，其中NA表示设备处于正常运行状态，A表示设备处于不平衡故障的异常运行状态亦即报警状态。

(2)设x为安装在推进电机传动轴上的振动传感器所检测的振动加速度变量，令x(t)，t＝1,2,3,…,是x在t时刻的采样值，单位为m/s²，利用式(1)对x(t)进行预处理，得到处理后的变量y及其采样值为：

定义y_otp为最优阈值，记max(y)和min(y)分别是y的最大值和最小值，其中y_otp∈[min(y),max(y)]，采用直接门限法最小化误报率和漏报率的平方和求得最优阈值y_otp；利用Sigmoid隶属度函数构建正常状态NA、异常状态A的模糊隶属度函数μ_NA(y(t))、μ_A(y(t))，如式(2)-(3)所示：

其中a∈[0.1,100]是一个可调参数。

为了加深对步骤(2)的理解，这里举个例子解释正常状态NA、异常状态A的模糊隶属度函数μ_NA(y(t))、μ_A(y(t))计算过程，假设a＝100，y_otp＝1.1324，y(1)＝1.2，那么μ_NA(y(1))＝0.0012。

(3)将y(t)分别带入步骤(2)关于正常状态NA、异常状态A的模糊隶属度函数μ_NA(y(t))、μ_A(y(t))中，即可得到报警证据m_t(NA)、m_t(A)，分别如下(4)-(5)式所示：

m_t(NA)＝μ_NA(y(t))/(μ_A(y(t))+μ_NA(y(t))) (4)

m_t(A)＝μ_A(y(t))/(μ_A(y(t))+μ_NA(y(t))) (5)

那么，在每个采样时刻，能够获得报警证据m_t＝(m_t(NA),m_t(A))，其中m_t(A)和m_t(NA)表示x(t)支持“报警A”和“正常NA”这两个事件发生的程度。

(4)在步骤(3)的基础上得到了每一时刻的报警证据之后，利用证据的动态更新得到全局报警证据，记为m_0:t＝(m_0:t(A),m_0:t(NA))，具体步骤如下：

(4-1)当t＝1时，因在t＝1之前没有相关的证据信息,所以t＝1时全局的报警证据是该时刻的报警证据即m_0:1＝(m₁(A),m₁(NA))。

(4-2)当t＝2时，利用t＝1时刻的全局报警证据m_0:1和t＝2时刻的证据m₂进行加权融合，得到t＝2时刻的全局报警证据m_0:2，具体用于计算的证据滤波公式如下

其中τ₂＝0.75表示对m_0:1(A)的线性加权值，ξ₂＝0.25表示当前证据的加权值；若m₂(NA)>m₂(A)，则D＝NA，那么m₂(A|NA)＝0；若m₂(A)≥m₂(NA)，则D＝A，那么m₂(A|A)＝1。

(4-3)当t＝3时利用t＝2时刻的全局报警证据m_0:2和t＝3时刻的证据m₃进行加权融合，得到t＝3时刻的全局报警证据m_0:3＝(m_0:3(A),m_0:3(NA))，具体用于计算的证据滤波公式如下：

若m₃(NA)>m₃(A)，则D＝NA，那么m₃(A|NA)＝0；若m₃(A)≥m₃(NA)，则D＝A，那么m₃(A|A)＝1；τ₃表示对m_0:2(A)的线性加权值，ξ₃表示当前证据的加权值，τ₃、ξ₃的计算步骤如下：

(4-3-1)定义m_0:1、m_0:2两条证据之间的距离如式(8)所示：

其中，上标T表示求向量(m_0:1-m_0:2)的转置，那么m_0:1和m_0:2两条证据之间的相似度如式(9)所示：

Sim(m_0:1,m_0:2)＝1-d_J(m_0:1,m_0:2) (9)

Sim(m_0:1,m_0:2)越大，两者越相似，同理可以计算出m₃、m_0:1和m_0:2之间两两相似度分别为Sim(m_0:1,m₃)以及Sim(m_0:2,m₃)。

(4-3-2)根据步骤(4-3-1)得到的相似度求出每个证据被其他两个证据支持的程度，分别为Sup(m_0:2)、Sup(m_0:1)以及Sup(m₃)，如下式(10)-(12)所示

Sup(m_0:2)＝Sim(m_0:1,m_0:2)+Sim(m_0:2,m₃) (10)

Sup(m_0:1)＝Sim(m_0:1,m_0:2)+Sim(m_0:1,m₃) (11)

Sup(m₃)＝Sim(m_0:1,m₃)+Sim(m_0:2,m₃) (12)

(4-3-3)在步骤(4-3-2)的基础上得到证据动态更新的权重因子τ₃和ξ₃，如下式(13)和(14)所示：

τ₃＝1-ξ₃ (14)

(4-4)当t≥4时，求取m_0:t＝(m_0:t(A),m_0:t(NA))，具体用于计算的证据滤波公式如下：

其中τ_t+ξ_t＝1，若m_t(NA)>m_t(A)，则D＝NA，那么m_t(A|NA)＝0；若m_t(A)≥m_t(NA)，则D＝A，那么m_t(A|A)＝1；同理可得m_t-1(A|D)，m_t-2(A|D)。

为了加深步骤(4)中t＝1,2,3时刻证据动态更新的过程的理解，这里举个例子加以解释，假设t＝1,2,3时刻对应的每个时刻对应的报警证据如下表2所示。

表2 振动数据t时刻报警证据

t	m<sub>t</sub>＝(m<sub>t</sub>(A),m<sub>t</sub>(NA))
		1	(0.1,0.9)
2	(0.3,0.7)
		3	(0.57,0.43)

根据步骤(4)计算3个时刻的全局报警证据：

当t＝1时，根据步骤(4-1)可得m_0：1＝m₁＝(0.1,0.9)。

当t＝2时，根据步骤(4-2)得到m₂(NA)>m₂(A)，则D＝NA，那么m₂(A|NA)＝0，m_0:2(A)＝τ₂m_0:1(A)+ξ₂m₂(A|D)＝0.75×0.1+0.25×0＝0.075，m_0:2(NA)＝0.925，从而得到m_0：2＝(0.075,0.925)。

当t＝3时，根据步骤(4-3)得到m₃(A)>m₃(NA)，则D＝A，那么m₃(A|A)＝1,根据子步骤a)可以得到Sim(m_0:1,m_0:2)＝0.975，Sim(m_0:1,m₃)＝0.53，Sim(m_0:2,m₃)＝0.505。再根据子步骤b)可以得到Sup(m_0:2)＝1.48，Sup(m_0:1)＝1.505，Sup(m₃)＝1.035。再根据子步骤c)可以得到ξ₃＝0.257，τ₃＝0.743，因此可以得到全局报警证据m_0:3(A)＝0.313，m_0:3(NA)＝0.687，从而得到m_0：3＝(0.313,0.687)。

(5)式(15)中和的具体确定步骤如下：

(5-1)当t≥4时，依据步骤(4-3)中的子步骤(4-3-1)提供的方法得到m_t、m_t-1和m_t-2两两之间的相似度Sim(m_t,m_t-1)、Sim(m_t,m_t-2)和Sim(m_t-1,m_t-2)。

(5-2)依据步骤(4-3)中的子步骤(4-3-2)提供的方法分别得到m_t、m_t-1和m_t-2相对于其他两个证据的支持程度Sup(m_t)、Sup(m_t-1)和Sup(m_t-2)。

(5-3)在步骤(5-2)的基础上，和可以由式(16)-(18)所示

为了便于理解，这里举例说明和的确定：

假设t＝1,2,3,4时刻对应的每个时刻对应的报警证据如下表3所示。

表3 振动数据t时刻报警证据

t	m<sub>t</sub>＝(m<sub>t</sub>(A),m<sub>t</sub>(NA))
		1	(0.1,0.9)
2	(0.3,0.7)
		3	(0.57,0.43)
4	(0.7,0.3)

当t＝4时，根据子步骤(5-1)可以得到Sim(m₄,m₃)＝0.87，Sim(m₄,m₂)＝0.6，Sim(m₃,m₂)＝0.73，再根据子步骤(5-2)可以得到Sup(m₄)＝1.47，Sup(m₃)＝1.6，Sup(m₂)＝1.33，再根据子步骤(5-2)可以得到

(6)对于步骤(4)中式(15)的参数ξ_t，计算过程如下：

(6-1)依据步骤(4-3)中的子步骤(4-3-1)分别得到m_t、m_t-1、m_t-2和m_0:t-3两两之间的相似度分别为Sim(m_t,m_t-1)、Sim(m_t,m_t-2)、Sim(m_t,m_0:t-3)、Sim(m_t-1,m_0:t-3)、Sim(m_t-1,m_t-2)以及Sim(m_t-2,m_0:t-3)。

(6-2)依据步骤(4-3)中的子步骤(4-3-2)分别得到m_t、m_t-1、m_t-2和m_0:t-3对于其他两个证据的支持程度，分别为Sup(m_0:t-3)、Sup(m_t)、Sup(m_t-1)以及Sup(m_t-2)。

(6-3)构造关于ξ_t的信度推理模型，用于描述输入量Sup(m_0:t-3)、Sup(m_t-1)、Sup(m_t-2)和Sup(m_t)与输出量线性加权值ξ_t之间存在的复杂非线性关系，具体步骤如下：

(6-3-1)将输入变量Sup(m_0:t-3)、Sup(m_t)、Sup(m_t-1)以及Sup(m_t-2)记为f_i(t)(i＝1,2,3,4)；设定ξ_t的参考值集合D＝{D_n|n＝1,...,N}，输入变量f_i(t)的参考值集1<N<10为ξ_t的结果参考值个数，1<J_i<10为输入变量的参考值个数。

(6-3-2)给定如表4所示的信度矩阵表来描述输入f_i(t)和输出ξ_t之间的关系，由表1中可知，当输入值f_i(t)取参考值时，结果值ξ_t为参考值D_n的信度为 并有定义对应于参考值的信度为同时给定输入f_i(t)的可靠度为r_i，满足0≤r_i≤1。

表4 输入f_i(t)的信度矩阵表

(6-3-3)当获得输入样本f₁(t)、f₂(t)、f₃(t)和f₄(t)之后，根据信度矩阵表1和可靠度r₁、r₂、r₃和r₄，可利用证据推理规则推理出线性加权值的估计值具体步骤如下：

a)对于输入值f_i(t)，其必然落入某两个参考值构成的区间此时这两个参考值对应的信度和被激活，α_i,j表示输入值f_i(t)匹配参考值的相似度，α_i,j+1表示输入值f_i(t)匹配参考值的相似度，由式(19)所示：

则输入值f_i(t)的信度可由和以加权和的形式获得：

e_i＝{(D_n,p_n,i),n＝1,...,N} (20)

b)利用式(20)和式(21)获得f₁(t)、f₂(t)、f₃(t)和f₄(t)的信度e₁、e₂、e₃和e₄，给定e_i的可靠度r_i为1，设e_i的证据权重初值w_i＝r_i＝1，利用证据推理规则对e₁、e₂、e₃和e₄进行融合，融合过程如下：

(a)首先对e₁和e₂进行融合，得到融合后的证据为

(b)将e₁和e₂融合后的概率证据设为e′₂＝{(D_n,p′_n,2),n＝1,...,N}，将其与e₃进行融合，其中令e′₂的证据权重w′₂为1，可靠度r′₂为1，得到融合后的概率证据为

(c)将e′₂和e₃融合后的概率证据设为e″₂＝{(D_n,p″_n,2),n＝1,...,N}，将其与e₄进行融合，其中令e″₂的证据权重w″₂为1，可靠度r″₂为1，作为最终的融合结果，即有：

O(f(t))＝{(D_n,p″_n,2),n＝1,...,N} (25)

(d)根据步骤(c)得到融合结果O(f(t))，线性加权值的估计值可由下式推理得到：

式(12)中线性加权值τ_t的估计值为

为了便于理解ξ_t的信度推理模型构造过程，这里举例加以说明，步骤如下：

假设当t＝4时，f₁(4)＝0.6，f₂(4)＝1.6，f₃(4)＝2.6，f₄(4)＝2.9，ξ_t的参考值集合D＝{1，2，3，4}，给定输入f₁(t)、f₂(t)、f₃(t)和f₄(t)的可靠度r₁＝r₂＝r₃＝r₄＝1，同时给定输入f_i(t)的信度矩阵表如表5～表8所示：

表5 输入f₁(t)的信度矩阵表

表6 输入f₂(t)的信度矩阵表

表7 输入f₃(t)的信度矩阵表

表8 输入f₄(t)的信度矩阵表

根据步骤(6-3-3)子步骤a)可以得到f₁(t)在与之间，f₁(4)匹配参考值的相似度α_1,1＝0.4，f₁(4)匹配参考值的相似度α_1,2＝0.6，激活证据和输入值f₁(4)的信度e₁＝{(0.1,0.13),(0.2,0.266),(0.3，0.36),(0.4,0.244)}；f₂(4)在与之间，f₂(4)匹配参考值的相似度α_2,2＝0.4，f₂(4)匹配参考值的相似度α_2,3＝0.6，激活证据和输入值f₂(4)的信度e₂＝{(0.1,0.242),(0.2,0.2),(0.3，0.32),(0.4,0.238)}；f₃(4)在与之间，f₃(4)匹配参考值的相似度α_1,1＝0.4，f₃(4)匹配参考值的相似度α_1,2＝0.6，激活证据和输入值f₃(4)的信度e₃＝{(0.1,0.284),(0.2,0.148),(0.3，0.17),(0.4,0.398)}；f₄(t)在与之间，f₄(t)匹配参考值的相似度α_1,1＝0.1，f₄(t)匹配参考值的相似度α_1,2＝0.9，激活证据和输入值f₄(t)的信度e₄＝{(0.1,0.374),(0.2,0.479),(0.3，0.045),(0.4,0.102)}；令e_i的证据权重初值w_i＝r_i＝1，根据子步骤(a)对e₁和e₂进行融合，得到融合后的概率证据p′_n,2＝{0.122，0.206，0.447，0.225}，再根据子步骤(b)得到融合后的概率证据p″_n,2＝{0.150，0.132，0.330，0.388},再根据子步骤(d)得到最终融合后的概率证据p″′_n,2＝{0.323，0.364，0.085，0.228}，从而得到线性加权值的估计值线性加权值τ_t的估计值为

(7)基于Jousselme证据距离构建对信度矩阵表1中参数优化模型，具体步骤如下：

(7-1)确定优化参数集合

(7-2)将最小化Jousselme证据距离作为优化目标函数

s.t.0≤w_i≤1,i＝1,2,3,4 (28)

D₂<D₃<…<D_N-1 (30)

其中L为采集的训练样本数量总数，m表示为设备所处的真实状态下的理想向量(010)或异常状态该下的理想向量(100)，式(28)-(31)表示优化参数需满足的约束条件。

(8)利用基于梯度下降法方法或者非线性优化软件包，获得最优的参数集合P，从而得到优化后的输入f_i(t)的信度矩阵表；从在线运行的推进电机传动轴上采集振动加速度信号，将其利用步骤(2)中式(1)进行处理，按照一定间隔选取阈值，再重复步骤(3)-(6)，可得到全局报警证据m_0:t＝(m_0:t(A),m_0:t(NA))；根据得到的t时刻的全局报警证据m_0:t＝(m_0:t(A),m_0:t(NA))，将m_0:t转化为相应的pignistic概率如下：

其中，所以基于pignistic概率报警决策准则为：若则报警器发出报警，反之，不发出报警；得到在不同阈值下误报率和漏报率的平方和，误报率和漏报率的平方和最小所对应的阈值即为最优阈值。

为了便于理解，这里举例说明根据全局报警证据进行决策，报警结果如下表9。

表9 报警结果输出

以下结合附图，详细介绍本发明方法的实施例：

本发明方法的流程图如图1所示，核心部分是：从推进电机传动轴采集振动加速度信号；将传感器采集的数据进行特定的处理；采用直接门限法最小化误报率和漏报率的平方和求得最优阈值y_otp，利用Sigmoid隶属度函数获取每一时刻的报警证据；利用多阶证据滤波对每时刻报警证据进行融合，构建信度推理模型描述多阶证据滤波公式中权重的复杂关系，给定输入与输出的信度矩阵表；利用证据推理规则融合输入样本向量激活的证据并从融合结果推理出多阶证据滤波公式中权重估计值；优化信度推理模型的参数；构建参数优化模型训练参考值及证据权重构成的参数集合，优化信度推理模型的参数集合，得到最优参数集合；在线采集数据，预处理后转化为每时刻报警证据，选取不同的阈值利用多阶证据滤波公式得到全局报警证据，进行报警决策，误报率和漏报率的平方和最小所对应的阈值为最优阈值。

以下结合ABB公司生产的船首侧推电机AMI450L6L，详细介绍本发明方法的各个步骤，通过实验数据在最优阈值的条件下，得到的误报率与漏报率要比在传统的方法上更加低。

1、实验数据的采集

侧推电机转速为1000转/分钟，采样频率为1280Hz，使用振动传感器在侧推电机传动轴上测得正常运行状态下得振动加速度1000组，异常运行状态下得振动加速度1000组。将这2000个数据组成排列成一段序列x(t)，t＝1,2,3,…,2000，如附图2所示。

2、求得最优阈值y_otp

根据本发明方法步骤(2)得到处理后的变量序列y(t)，t＝1,2,3,…,2000，如附图3所示。采用直接门限法最小化误报率和漏报率的平方和求得最优阈值y_otp＝1.1324。利用Sigmoid隶属度函数构建正常状态NA、异常状态A的模糊隶属度函数μ_NA(y(t))、μ_A(y(t))，其中a＝15。

3、获取每一时刻的报警证据m_t＝(m_t(NA),m_t(A))

根据本发明方法步骤(3)，a＝15,可以得到每时刻的报警证据。

4、前三个时刻全局报警证据的获取

根据本发明方法步骤(4)可以得到t＝1,2,3时刻的全局报警证据。

5、和的确定

根据本发明方法步骤(5)可以得到t ≥4后每一步证据滤波的和值。

4、输入值f_i(t)、结果ξ_t参考值的选取和初始信度矩阵表的构建

根据本发明方法步骤(6)设定输入变量f₁的参考值集合A₁＝{0,0.5,1,1.5,2,2.5,3}，f₂的参考值集合A₂＝{0,0.5,1,1.5,2,2.5,3}，f₃的参考值集合A₃＝{0,0.5,1,1.5,2,2.5,3}，f₄的参考值集合A₄＝{0,0.5,1,1.5,2,2.5,3}，ξ_t的参考值集合D＝{0,0.1,0.15,0.2,0.25,0.3 0.35}，初始信度矩阵表如表10～表13所示，同时给定输入f_i(t)的可靠度为r_i＝1。

表10 输入f₁(t)的初始信度矩阵表

表11 输入f₂(t)的初始信度矩阵表

表12 输入f₃(t)的初始信度矩阵表

表13 输入f₄(t)的初始信度矩阵表

根据本发明方法步骤(6-3-3)利用证据推理规则推理出线性加权值的估计值例如输入f₁(2000)＝2.2946，f₂(2000)＝2.2945，f₃(2000)＝2.2946，f₄(2000)＝0.9875，输入f₁(t)以相似度α_1，5＝0.4109，α_1，6＝0.5891激活证据和输入f₂(t)以相似度α_2，5＝0.4109，α_2，6＝0.5891激活证据和输入f₃(t)以相似度α_3，5＝0.4109，α_3，6＝0.5891激活证据和输入f₄(t)以相似度α_4，2＝0.0251，α_4，3＝0.9749激活证据和得到e₁＝{(0,0.1372),(0.04,0.1382),(0.08,0.0899),(0.12,0.1125),(0.16,0.1790),(0.20,0.1740),(0.26,0.1692)}，e₂＝{(0,0.1793),(0.04,0.1130),(0.08,0.0945),(0.12,0.1166),(0.16,0.1841),(0.20,0.1594),(0.26,0.1531)}，e₃＝{(0,0.1891),(0.04,0.1185),(0.08,0.0993),(0.12,0.1091),(0.16,0.1807),(0.20,0.1554),(0.26,0.1479)}，e₄＝{(0,0.1784),(0.04,0.0866),(0.08,0.1937),(0.12,0.1347),(0.16,0.2171),(0.20,0.1136),(0.26,0.0759)}，根据式(22)-(25)，可以得到最终融合结果为：Of(t)＝{(0,2428),(0.04,0.0468),(0.08,0.0478),(0.12,0.0564),(0.16,0.3780),(0.20,0.1431),(0.26,0.0851)}。将上述融合结果代入(26)式，即可得线性加权值的估计值

5、根据本发明方法步骤(7)构建参数优化模型，可得优化后的信度矩阵表，分别如下表14～表17所示。

表14 输入f₁(t)优化后的信度矩阵表

表15 输入f₂(t)优化后的信度矩阵表

表16 输入f₃(t)优化后的信度矩阵表

表17 输入f₄(t)优化后的信度矩阵表

6、在线采集振动数据x(t)，在振动数据x′(t)中，x′(1)到x′(500)和x′(1001)到x′(1500)这1000个采样数据处于正常的工作状态，x′(501)到x′(1000)和x′(1501)到x′(2000)这1000个采样数据处于异常的工作状态，如附图所示。将这两组共2000个采样数据排列成一段序列X(t)，如附图4所示，前1000个是设备处于正常运行状态下的数据，后1000个是设备处于异常状态下的数据。根据本发明步骤(2)处理振动数据，得到Y(t)，如附图5所示，按0.1的间隔选取不同的阈值利用多阶证据滤波公式得到全局报警证据，进行报警决策，误报率和漏报率的平方和最小所对应的最优阈值为1.0549，此时对应的误报率与漏报率是FAR＝8.7％，MAR＝6.8％。采用直接门限法最小化误报率和漏报率的平方和求得的误报率与漏报率FAR＝33％，MAR＝38.1％，此时最优阈值为1.1549，从最后的报警结果可以看出，本发明在最优阈值的下产生的误报率和漏报率都比传统的方法更加精确。

Claims (1)

1.一种基于证据滤波的电推船电机转轴不平衡故障检测方法，其特征在于该方法包括以下步骤：

(1)设定中压电推船舶推进电机传动轴不平衡故障报警器的辨识框架为Θ＝{NA,A}，其中NA表示设备处于正常运行状态，A表示设备处于不平衡故障的异常运行状态亦即报警状态；

(2)设x为安装在推进电机传动轴上的振动传感器所检测的振动加速度变量，令x(t)，t＝1,2,3,…,是x在t时刻的采样值，利用式(1)对x(t)进行预处理，得到处理后的变量y及其采样值为：

定义y_otp为最优阈值，记max(y)和min(y)分别是y的最大值和最小值，其中y_otp∈[min(y),max(y)]，采用直接门限法最小化误报率和漏报率的平方和求得最优阈值y_otp；利用Sigmoid隶属度函数构建正常状态NA、异常状态A的模糊隶属度函数μ_NA(y(t))、μ_A(y(t))，如式(2)-(3)所示：

其中a∈[0.1,100]是一个可调参数；

(3)将y(t)分别带入步骤(2)关于正常状态NA、异常状态A的模糊隶属度函数μ_NA(y(t))、μ_A(y(t))中，即可得到报警证据m_t(NA)、m_t(A)，分别如下(4)-(5)式所示：

m_t(NA)＝μ_NA(y(t))/(μ_A(y(t))+μ_NA(y(t))) (4)

m_t(A)＝μ_A(y(t))/(μ_A(y(t))+μ_NA(y(t))) (5)

那么，在每个采样时刻，能够获得报警证据m_t＝(m_t(NA),m_t(A))，其中m_t(A)和m_t(NA)表示x(t)支持“报警A”和“正常NA”这两个事件发生的程度；

(4)在步骤(3)的基础上得到了每一时刻的报警证据之后，利用证据的动态更新得到全局报警证据，记为m_0:t＝(m_0:t(A),m_0:t(NA))，具体步骤如下：

(4-1)当t＝1时，因在t＝1之前没有相关的证据信息,所以t＝1时全局的报警证据是该时刻的报警证据即m_0:1＝(m₁(A),m₁(NA))；

(4-2)当t＝2时，利用t＝1时刻的全局报警证据m_0:1和t＝2时刻的证据m₂进行加权融合，得到t＝2时刻的全局报警证据m_0:2，具体用于计算的证据滤波公式如下

其中τ₂＝0.75表示对m_0:1(A)的线性加权值，ξ₂＝0.25表示当前证据的加权值；若m₂(NA)>m₂(A)，则D＝NA，那么m₂(A|NA)＝0；若m₂(A)≥m₂(NA)，则D＝A，那么m₂(A|A)＝1；

(4-3)当t＝3时利用t＝2时刻的全局报警证据m_0:2和t＝3时刻的证据m₃进行加权融合，得到t＝3时刻的全局报警证据m_0:3＝(m_0:3(A),m_0:3(NA))，具体用于计算的证据滤波公式如下

若m₃(NA)>m₃(A)，则D＝NA，那么m₃(A|NA)＝0；若m₃(A)≥m₃(NA)，则D＝A，那么m₃(A|A)＝1；τ₃表示对m_0:2(A)的线性加权值，ξ₃表示当前证据的加权值，τ₃、ξ₃的计算步骤如下：

(4-3-1)定义m_0:1、m_0:2两条证据之间的距离如式(8)所示：

其中，上标T表示求向量(m_0:1-m_0:2)的转置，那么m_0:1和m_0:2两条证据之间的相似度如式(9)所示：

Sim(m_0:1,m_0:2)＝1-d_J(m_0:1,m_0:2) (9)

Sim(m_0:1,m_0:2)越大，两者越相似，同理可以计算出m₃、m_0:1和m_0:2之间两两相似度分别为Sim(m_0:1,m₃)以及Sim(m_0:2,m₃)；

(4-3-2)根据步骤(4-3-1)得到的相似度求出每个证据被其他两个证据支持的程度，分别为Sup(m_0:2)、Sup(m_0:1)以及Sup(m₃)，如下式(10)-(12)所示

Sup(m_0:2)＝Sim(m_0:1,m_0:2)+Sim(m_0:2,m₃) (10)

Sup(m_0:1)＝Sim(m_0:1,m_0:2)+Sim(m_0:1,m₃) (11)

Sup(m₃)＝Sim(m_0:1,m₃)+Sim(m_0:2,m₃) (12)

(4-3-3)在步骤(4-3-2)的基础上得到证据动态更新的权重因子τ₃和ξ₃，如下式(13)和(14)所示

τ₃＝1-ξ₃ (14)

(4-4)当t≥4时，求取m_0:t＝(m_0:t(A),m_0:t(NA))，具体用于计算的证据滤波公式如下

其中τ_t+ξ_t＝1，若m_t(NA)>m_t(A)，则D＝NA，那么m_t(A|NA)＝0；若m_t(A)≥m_t(NA)，则D＝A，那么m_t(A|A)＝1；同理可得m_t-1(A|D)，m_t-2(A|D)；

(5)式(15)中和的具体确定步骤如下：

(5-1)当t≥4时，依据步骤(4-3)中的子步骤(4-3-1)提供的方法得到m_t、m_t-1和m_t-2两两之间的相似度Sim(m_t,m_t-1)、Sim(m_t,m_t-2)和Sim(m_t-1,m_t-2)；

(5-2)依据步骤(4-3)中的子步骤(4-3-2)提供的方法分别得到m_t、m_t-1和m_t-2相对于其他两个证据的支持程度Sup(m_t)、Sup(m_t-1)和Sup(m_t-2)；

(5-3)在步骤(5-2)的基础上，和可以由式(16)-(18)所示

(6)对于步骤(4)中式(15)的参数ξ_t，计算过程如下：

(6-1)依据步骤(4-3)中的子步骤(4-3-1)分别得到m_t、m_t-1、m_t-2和m_0:t-3两两之间的相似度分别为Sim(m_t,m_t-1)、Sim(m_t,m_t-2)、Sim(m_t,m_0:t-3)、Sim(m_t-1,m_0:t-3)、Sim(m_t-1,m_t-2)以及Sim(m_t-2,m_0:t-3)；

(6-2)依据步骤(4-3)中的子步骤(4-3-2)分别得到m_t、m_t-1、m_t-2和m_0:t-3对于其他两个证据的支持程度，分别为Sup(m_0:t-3)、Sup(m_t)、Sup(m_t-1)以及Sup(m_t-2)；

(6-3)构造关于ξ_t的信度推理模型，用于描述输入量Sup(m_0:t-3)、Sup(m_t-1)、Sup(m_t-2)和Sup(m_t)与输出量线性加权值ξ_t之间存在的复杂非线性关系，具体步骤如下：

(6-3-1)将输入变量Sup(m_0:t-3)、Sup(m_t)、Sup(m_t-1)以及Sup(m_t-2)记为f_i(t)(i＝1,2,3,4)；设定ξ_t的参考值集合D＝{D_n|n＝1,...,N}，输入变量f_i(t)的参考值集1<N<10为ξ_t的结果参考值个数，1<J_i<10为输入变量的参考值个数；

(6-3-2)给定如表1所示的信度矩阵表来描述输入f_i(t)和输出ξ_t之间的关系，由表1中可知，当输入值f_i(t)取参考值时，结果值ξ_t为参考值D_n的信度为 并有定义对应于参考值的信度为同时给定输入f_i(t)的可靠度为r_i，满足0≤r_i≤1；

表1 输入f_i(t)的信度矩阵表

(6-3-3)当获得输入样本f₁(t)、f₂(t)、f₃(t)和f₄(t)之后，根据信度矩阵表1和可靠度r₁、r₂、r₃和r₄，可利用证据推理规则推理出线性加权值的估计值具体步骤如下：

a)对于输入值f_i(t)，其必然落入某两个参考值构成的区间此时这两个参考值对应的信度和被激活，α_i,j表示输入值f_i(t)匹配参考值的相似度，α_i,j+1表示输入值f_i(t)匹配参考值的相似度，由式(19)所示

则输入值f_i(t)的信度可由和以加权和的形式获得

e_i＝{(D_n,p_n,i),n＝1,...,N} (20)

b)利用式(20)和式(21)获得f₁(t)、f₂(t)、f₃(t)和f₄(t)的信度e₁、e₂、e₃和e₄，给定e_i的可靠度r_i为1，设e_i的证据权重初值w_i＝r_i＝1，利用证据推理规则对e₁、e₂、e₃和e₄进行融合，融合过程如下：

(a)首先对e₁和e₂进行融合，得到融合后的证据为

(b)将e₁和e₂融合后的概率证据设为e′₂＝{(D_n,p′_n,2),n＝1,...,N}，将其与e₃进行融合，其中令e′₂的证据权重w′₂为1，可靠度r′₂为1，得到融合后的概率证据为

(c)将e′₂和e₃融合后的概率证据设为e″₂＝{(D_n,p″_n,2),n＝1,...,N}，将其与e₄进行融合，其中令e″₂的证据权重w″₂为1，可靠度r″₂为1，作为最终的融合结果，即有：

O(f(t))＝{(D_n,p″′_n,2),n＝1,...,N} (25)

(d)根据步骤(c)得到融合结果O(f(t))，线性加权值的估计值可由下式推理得到

式(12)中线性加权值τ_t的估计值为

(7)基于Jousselme证据距离构建对信度矩阵表1中参数优化模型，具体步骤如下：

(7-1)确定优化参数集合

(7-2)将最小化Jousselme证据距离作为优化目标函数

s.t. 0≤w_i≤1,i＝1,2,3,4 (28)

D₂<D₃<…<D_N-1 (30)

其中L为采集的训练样本数量总数，m表示为设备所处的真实状态下的理想向量(010)或异常状态该下的理想向量(100)，式(28)-(31)表示优化参数需满足的约束条件；

(8)利用基于梯度下降法方法或者非线性优化软件包，获得最优的参数集合P，从而得到优化后的输入f_i(t)的信度矩阵表；从在线运行的推进电机传动轴上采集振动加速度信号，将其利用步骤(2)中式(1)进行处理，按照一定间隔选取阈值，再重复步骤(3)-(6)，可得到全局报警证据m_0:t＝(m_0:t(A),m_0:t(NA))；根据得到的t时刻的全局报警证据m_0:t＝(m_0:t(A),m_0:t(NA))，将m_0:t转化为相应的pignistic概率如下

其中，所以基于pignistic概率报警决策准则为：若则报警器发出报警，反之，不发出报警；得到在不同阈值下误报率和漏报率的平方和，误报率和漏报率的平方和最小所对应的阈值即为最优阈值。