CN113051092A - 一种基于优化核密度估计和js散度的故障诊断方法 - Google Patents

一种基于优化核密度估计和js散度的故障诊断方法 Download PDF

Info

Publication number
CN113051092A
CN113051092A CN202110158768.1A CN202110158768A CN113051092A CN 113051092 A CN113051092 A CN 113051092A CN 202110158768 A CN202110158768 A CN 202110158768A CN 113051092 A CN113051092 A CN 113051092A
Authority
CN
China
Prior art keywords
divergence
data
training data
kernel density
density estimation
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
CN202110158768.1A
Other languages
English (en)
Other versions
CN113051092B (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.)
National University of Defense Technology
Original Assignee
National University of Defense Technology
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 National University of Defense Technology filed Critical National University of Defense Technology
Priority to CN202110158768.1A priority Critical patent/CN113051092B/zh
Publication of CN113051092A publication Critical patent/CN113051092A/zh
Application granted granted Critical
Publication of CN113051092B publication Critical patent/CN113051092B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F11/00—Error detection; Error correction; Monitoring
    • G06F11/07—Responding to the occurrence of a fault, e.g. fault tolerance
    • G06F11/0703—Error or fault processing not based on redundancy, i.e. by taking additional measures to deal with the error or fault not making use of redundancy in operation, in hardware, or in data representation
    • G06F11/079—Root cause analysis, i.e. error or fault diagnosis
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/21—Design or setup of recognition systems or techniques; Extraction of features in feature space; Blind source separation
    • G06F18/214—Generating training patterns; Bootstrap methods, e.g. bagging or boosting

Landscapes

  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • General Physics & Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Evolutionary Biology (AREA)
  • Evolutionary Computation (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Artificial Intelligence (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Health & Medical Sciences (AREA)
  • Biomedical Technology (AREA)
  • Quality & Reliability (AREA)
  • Complex Calculations (AREA)

Abstract

本发明提供一种基于优化核密度估计和JS散度的故障诊断方法，用于诊断具有平稳分布特性的硬件设备的故障。首先，该方法先估算最优核密度带宽，从而得到最优核密度函数估计；其次，通过滑动采样窗口导出了样本密度分布和总体密度分布之间JS散度的分布特征；最后，给出基于JS散度的故障检测阈值与隔离阈值，进而构建了一种基于优化核密度估计和JS散度的故障诊断方法。针对轴承等具有平稳分布特性的硬件，该方法可以有效提高故障诊断的性能。

Description

一种基于优化核密度估计和JS散度的故障诊断方法
技术领域
本发明涉及故障检测领域，具体涉及一种基于优化核密度估计和JS散度的故障诊断方法。
背景技术
随着工业信息化的发展，各个领域都开始出现海量的数据，对于这些数据的处理成为了业内的难点问题，尤其是在故障诊断等领域。实际上，数据量爆炸式增长提供了更多的信息，在这种情况下，典型的数据分析理论存在应用弊端，其主要原因在于：典型的数据分析往往通过先验信息给出数据的分布类型，并在此假设基础上进行分析，但是一旦假设给出，后续工作只停留在参数的估计和分析上，而无法对假设本身进行修正。在故障诊断领域，本质问题是衡量样本之间的差异。通常使用频率直方图的方式表现两个样本的分布差异，但是该方法存在三个不足：一是大量离散操作较为浪费时间；二是离散间隔选取有较大主观性；三是没有直观指标反映差异大小。以滚动轴承为例，作为机械设备的关键部件，其发生故障会对设备的安全平稳运行造成严重影响，而对滚动轴承的早期故障检测可以避免设备带故障运行，避免造成严重的安全事故和经济损失，具有重要的现实意义和工程意义。与传统的故障诊断相比，滚动轴承的故障诊断更复杂，其主要表现在以下三点：第一、故障信号微弱。轴承数据通常是一种高频数据，而故障信号往往被这些高频信号所掩盖，导致传统的故障诊断方法失效。第二、数据高耦合。轴承数据通常以振动信号的形式反映出来，在不同维的信号中存在强耦合性，使得故障诊断存在较大的难点。第三、数据不均衡。滚动轴承多在正常状态下工作，能收集到的故障数据往往较少，使得数据不均衡，导致故障数据集不够完善，加大了故障检测的难度。
为解决这些问题，提出了基于趋势剔除和噪声消减的故障检测技术，其通过剔除趋势来增强信号趋势比，通过噪音消减增强信号噪声比，从而提高故障检测效果。但是该方法仍沿用传统的T2检测方法，并不能有效解决数据之间的耦合问题。还提出了基于PCA降维和模态分解特征提取的故障检测方法，对于高维数据先进行PCA降维处理，使得数据维数降低并消除不同维数之间的相关性，随后利用模态分解的方法提取各维度之间的特征进行故障检测。该方法有效解决了数据之间的强耦合性，但是在PCA降维处理过程中会损耗部分信息，导致故障检测效果降低。
发明内容
本发明实施例提供一种基于优化核密度估计和JS散度的故障诊断方法，通过基于最优带宽的核密度函数估计和JS散度构建了一种设备故障检测和辨识方法，并通过滑动采样窗口的方法，导出了样本密度分布和总体密度分布之间JS散度的分布特征，并基于此给出了故障检测的阈值，从而实现对不同故障的辨识,有效提高了设备故障诊断的性能。
为达到上述目的，一方面，本发明实施例提供了一种基于优化核密度估计和JS散度的故障诊断方法，包括：
通过传感器采集设备工作期间的运行数据，将采集到的运行数据作为待检测数据；
根据所述待检测数据和训练数据集，计算得到所述待检测数据与所述训练数据集中的每一训练数据之间的JS散度值，并且将获得的JS散度值构成第一JS散度集合；其中，所述训练数据集中的各训练数据是由传感器采集的设备工作期间的运行数据；所述训练数据集中的各训练数据与设备的已知的各状态模式标签相对应；所述状态模式标签用于标识设备的工作状态；
将所述第一JS散度集合中最小JS散度值对应的所述训练数据集中的训练数据，作为待选训练数据，并将所述最小JS散度值作为待选JS散度值；
利用滑窗原理对所述待选训练数据进行采样，得到滑窗数据集，计算所述滑窗数据集和所述待选训练数据之间的JS散度分布，并使用核密度估计的方法获得JS散度上界值；所述JS散度上界值用于作为设备故障诊断的检测阈值；以及，
根据所述待选JS散度值和所述JS散度上界值，确定所述待检测数据对应的设备的状态模式标签，具体包括：
若所述待选JS散度值小于或者等于所述JS散度上界值，则所述待检测数据对应的状态模式标签与所述待选训练数据对应的状态模式标签相同，或者，
若所述待选JS散度值状态大于所述JS散度上界值，则所述待检测数据对应设备的新工作状态。
进一步地，所述根据所述待检测数据和所述训练数据集，计算得到所述待检测数据与所述训练数据集的每一训练数据之间的JS散度值，并且将获得的JS散度值构成第一JS散度集合，包括：
将所述训练数据集的每一训练数据作为当前输入数据，将当前输入数据输入给基于最优带宽的核密度估计过程，获得最优核密度估计，并将当前获得的最优核密度估计作为所述训练数据集的当前训练数据相应的最优核密度估计；
将所述待检测数据作为当前输入数据，将当前输入数据输入给基于最优带宽的核密度估计过程，获得最优核密度估计，并将当前获得的最优核密度估计作为所述待检测数据的最优核密度估计；以及，
根据所述待检测数据的最优核密度估计和所述训练数据集的各训练数据各自相应的最优核密度估计，计算所述待检测数据与所述训练数据集中的每一训练数据之间的JS散度值，并且由计算得到的JS散度值构成第一JS散度集合。
进一步地，所述利用滑窗原理对所述待选训练数据进行采样，得到滑窗数据集，计算所述滑窗数据集和所述待选训练数据之间的JS散度分布，并使用核密度估计的方法获得JS散度上界值，包括：
根据指定宽度的滑动窗口在所述待选训练数据上滑动选取数据，得到至少一个滑窗训练数据，由所述至少一个滑窗训练数据构成所述滑窗数据集；
将所述滑窗数据集的每一滑窗训练数据作为当前输入数据，将当前输入数据输入给基于最优带宽的核密度估计过程，获得最优核密度估计，并将当前获得的最优核密度估计作为所述滑窗数据集的当前滑窗训练数据相应的最优核密度估计；
将所述待选训练数据作为当前输入数据，将当前输入数据输入给基于最优带宽的核密度估计过程，获得最优核密度估计，并将当前获得的最优核密度估计作为所述待选训练数据相应的最优核密度估计；
根据所述滑窗数据集的各滑窗训练数据各自相应的最优核密度估计和所述待选训练数据相应的最优核密度估计，计算所述滑窗数据集的每一滑窗训练数据与所述待选训练数据之间的JS散度值，由计算得到的JS散度值构成第二JS散度集合；
估计第二JS散度集合对应的JS散度密度函数；以及，
根据所述JS散度密度函数和指定的显著性水平值，获得JS散度上界值。
进一步地，所述将当前输入数据输入给基于最优带宽的核密度估计过程，获得最优核密度估计，包括：
设置带宽的初始值、给定的估计精度和最大迭代次数，并循环执行后续步骤，直到第一跳出条件或者第二跳出条件之一满足时，跳出循环；
根据带宽、核函数K(·)和核密度估计公式
Figure BDA0002934840890000021
计算核密度估计；
根据以下公式计算当前循环中的带宽：
Figure BDA0002934840890000022
判断若满足第一跳出条件，则已经得到最优带宽，跳出循环；
判断若满足第二跳出条件，则迭代次数已经超限，跳出循环；
保留当前循环计算得到的带宽，用于在下次循环中判断第一跳出条件；
其中：h表示带宽；K(·)表示核函数；rj表示当前输入数据中的第j个元素；第一跳出条件为当前循环计算得到的带宽减去上次循环计算得到的带宽，所得的差值的绝对值小于给定的估计精度；第二跳出条件为循环次数达到最大迭代次数。
进一步地，所述根据所述待检测数据的最优核密度估计和所述训练数据集的各训练数据各自相应的最优核密度估计，计算所述待检测数据与所述训练数据集的每一训练数据之间的JS散度值，并且由计算得到的JS散度值构成第一JS散度集合，包括：
依据以下公式，计算所述待检测数据与所述训练数据集的每一训练数据之间的JS散度值:
Figure BDA0002934840890000031
由计算得到的JS散度值构成第一JS散度集合：
{JS(Z,R1),JS(Z,R2),JS(Z,R3),…，JS(Z,Rq)}
其中：Z是待检测数据；Ri是训练数据集{R1,R2,R3,…，Rq}中的训练数据；
Figure BDA0002934840890000032
是训练数据的最优核密度估计；
Figure BDA0002934840890000033
是待检测数据的最优核密度估计。
进一步地，所述根据所述滑窗数据集的各滑窗训练数据各自相应的最优核密度估计和所述待选训练数据相应的最优核密度估计，计算所述滑窗数据集的每一滑窗训练数据与所述待选训练数据之间的JS散度值，由计算得到的JS散度值构成第二JS散度集合，包括:
根据以下公式，计算所述滑窗数据集的每一滑窗训练数据与所述待选训练数据之间的JS散度值:
Figure BDA0002934840890000034
由计算得到的JS散度值构成第二JS散度集合:
{JS1,JS2,JS3,…,JSm-p}
其中：R(j)表示滑窗数据集中的第j个滑窗训练数据；R表示待选训练数据；
Figure BDA0002934840890000035
是待选训练数据对应的最优核密度估计；
Figure BDA0002934840890000036
是第j个滑窗训练数据对应的最优核密度估计；H(·)表示求熵运算。
进一步地，所述估计第二JS散度集合对应的JS散度密度函数，包括：
根据如下公式估计第二JS散度集合对应的JS散度密度函数：
Figure BDA0002934840890000037
其中：JSj是第二JS散度集合中的第j个元素；K(·)是核函数；h是带宽；m是待选训练数据中的元素个数；p是以元素个数为单位的滑动窗口的宽度。
进一步地，所述根据所述JS散度密度函数和指定的显著性水平值，获得JS散度上界值，包括：
通过对如下公式做数值积分得到JS散度上界值：
Figure BDA0002934840890000038
JShigh＝h*i
其中：h是步长；i是指步长编号；
Figure BDA0002934840890000039
是JS散度密度函数；α是指定的显著性水平值；JShigh是JS散度上界值。
进一步地，在所述待选JS散度值大于所述JS散度上界值，则所述待检测数据对应于设备的新工作状态，之后还包括：
定义所述待检测数据对应的设备的新状态模式的状态模式标签，将所述待检测数据加入到所述训练数据集中。
区别于现有技术，上述技术方案具有如下有益效果：
根据本发明的上述技术方案，针对设备故障诊断的问题，提出了一种基于优化核密度估计和JS散度的故障诊断方法，通过基于最优带宽的核密度函数估计和JS散度构建了一种设备故障检测和辨识方法，将核密度估计方法拓展到高维数据上，避免了针对各维度单独进行核密度估计时造成的信息损失，从而更好地刻画数据的密度概率分布，同时，改进了传统方法中使用交叉熵函数作为密度分布差异的度量方法，采用了JS散度作为密度分布差异的度量，规避了采用交叉熵函数作为度量导致的相对性并通过滑动采样窗口的方法，导出了样本密度分布和总体密度分布之间JS散度的分布特征，并基于此给出了故障检测的阈值，从而实现对不同故障的辨识,有效提高了设备故障诊断的性能。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1为本发明实施例之一的基于优化核密度估计和JS散度的故障诊断方法的一种流程图；
图2为本发明实施例之一的基于最优带宽的核密度估计方法流程图；
图3为本发明实施例之一的基于最优带宽的故障诊断方法流程图；
图4为本发明实施例之一的驱动端加速度正常数据功率谱分布图；
图5为本发明实施例之一的风扇端加速度正常数据功率谱分布图；
图6为本发明实施例之一的驱动端加速度正常数据预处理后的幅值图；
图7为本发明实施例之一的风扇端加速度正常数据预处理后的幅值图；
图8为本发明实施例之一的训练集正常数据之一的驱动端幅值图；
图9为本发明实施例之一的训练集正常数据之一的风扇端幅值图；
图10为本发明实施例之一的训练集0.007英寸内滚道故障数据之一的驱动端幅值图；
图11为本发明实施例之一的训练集0.007英寸内滚道故障数据之一的风扇端幅值图；
图12为本发明实施例之一的训练集0.014英寸内滚道故障数据之一的驱动端幅值图；
图13为本发明实施例之一的训练集0.014英寸内滚道故障数据之一的风扇端幅值图；
图14为本发明实施例之一的训练集正常数据的二维频率直方图；
图15为本发明实施例之一的训练集0.007英寸内滚道故障数据二维频率直方图；
图16为本发明实施例之一的训练集0.014英寸内滚道故障数据二维频率直方图；
图17为本发明实施例之一的训练集正常数据的二维核密度估计示意图；
图18为本发明实施例之一的训练集0.007英寸内滚道故障数据二维核密度估计图；
图19为本发明实施例之一的训练集0.014英寸内滚道故障数据二维核密度估计图；
图20为本发明实施例之一的训练集正常数据的JS散度以及分布的核密度估计示意图；
图21为本发明实施例之一的训练集0.007英寸内滚道故障数据的JS散度以及分布的核密度估计图；
图22为本发明实施例之一的训练集0.014英寸内滚道故障数据的JS散度以及分布的核密度估计图；
图23为本发明实施例之一的测试集正常数据采用交叉熵函数的检测结果示意图；
图24为本发明实施例之一的测试集0.007英寸内滚道故障数据采用交叉熵函数的检测结果示意图；
图25为本发明实施例之一的测试集0.014英寸内滚道故障数据采用交叉熵函数的检测结果示意图；
图26为本发明实施例之一的测试集正常数据采用本发明方法的检测结果示意图；
图27为本发明实施例之一的测试集0.007英寸内滚道故障数据采用本发明方法的检测结果示意图；
图28为本发明实施例之一的测试集0.014英寸内滚道故障数据采用本发明方法的检测结果示意图；
图29为本发明实施例之一的0.014英寸外滚道故障数据驱动端加速度数据幅值图；
图30为本发明实施例之一的0.014英寸外滚道故障数据风扇端加速度数据幅值图；
图31为本发明实施例之一的不同窗宽下的检测效果对比图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
如图1所示，本发明提供的一种基于优化核密度估计和JS散度的故障诊断方法，包括：
步骤101，通过传感器采集设备工作期间的运行数据，将采集到的运行数据作为待检测数据；
步骤102，根据所述待检测数据和训练数据集，计算得到所述待检测数据与所述训练数据集中的每一训练数据之间的JS散度值，并且将获得的JS散度值构成第一JS散度集合；其中，所述训练数据集中的各训练数据是由传感器采集的设备工作期间的运行数据；所述训练数据集中的各训练数据与设备的已知的各状态模式标签相对应；所述状态模式标签用于标识设备的工作状态；
步骤103，将所述第一JS散度集合中最小JS散度值对应的所述训练数据集中的训练数据，作为待选训练数据，并将所述最小JS散度值作为待选JS散度值；
步骤104，利用滑窗原理对所述待选训练数据进行采样，得到滑窗数据集，计算所述滑窗数据集和所述待选训练数据之间的JS散度分布，并使用核密度估计的方法获得JS散度上界值；所述JS散度上界值用于作为设备故障诊断的检测阈值；以及，
步骤105，根据所述待选JS散度值和所述JS散度上界值，确定所述待检测数据对应的设备的状态模式标签，具体包括：
若所述待选JS散度值小于或者等于所述JS散度上界值，则所述待检测数据对应的状态模式标签与所述待选训练数据对应的状态模式标签相同，或者，
若所述待选JS散度值状态大于所述JS散度上界值，则所述待检测数据对应设备的新工作状态。
在步骤101中，通过传感器采集设备在工作期间的运行数据，设备包括运行数据具有平稳分布特性的硬件设备，传感器可以包括加速度计、摄像头、红外测距，激光测距，超声波测距、麦克风、红外测温、电流计、电压计等中的一种或多种，采集设备的数据可以包括震动、声音、发热、电流、电压等一种或多种数据；将采集到的运行数据作为待检测数据，通过后面的步骤，根据待检测数据识别设备当前的工作状态，设备的工作状态可以包括正常状态以及一种或多种的故障状态。在步骤102中，根据所述待检测数据和训练数据集，计算得到所述待检测数据与所述训练数据集中的每一训练数据之间的JS散度值，并且将获得的JS散度值构成第一JS散度集合；其中，所述训练数据集中的各训练数据是由传感器采集的设备工作期间的运行数据；所述训练数据集中的各训练数据与设备的已知的各状态模式标签相对应；所述状态模式标签用于标识设备的工作状态；训练数据集是在步骤102之前预先准备好的，训练数据集中包括一个或多个训练数据，每个训练数据也是通过传感器采集设备在工作期间的运行数据，每个训练数据可以包含一个数据元素或者指定时间范围内的连续或离散的设备运行数据；每个训练数据都与设备的一种状态模式标签对应；在训练数据集中的多个训练数据可以各自对应于不同的状态模式标签，也可以对应于同一个状态模式标签；一种状态模式标签对应的设备的一种工作状态；工作状态可以包括正常状态和一种或多种的故障状态；步骤103，将所述第一JS散度集合中的最小JS散度值对应的所述训练数据集中的训练数据，作为待选训练数据，并将所述最小JS散度值作为待选JS散度值；从步骤102中得到的第一JS散度集合中获得最小JS散度值，从而确定了在训练数据集中，与待检测数据的特征最接近的训练数据，并将这个训练数据作为待选训练数据，通过后面的步骤继续判断待检测数据是否与待选训练数据都对应于同一个状态模式标签；步骤104通过滑动采样窗口的方法，导出了样本密度分布和总体密度分布之间JS散度的分布特征，并基于此给出了JS散度上界值，即故障检测的阈值，步骤105，根据所述待选JS散度值和所述JS散度上界值，确定所述待检测数据对应的设备的状态模式标签，若所述待选JS散度值小于或者等于所述JS散度上界值，则所述待检测数据与所述待选训练数据对应相同的状态模式标签，从而可以确定待检测数据对应的设备工作状态；若所述待选JS散度值状态大于所述JS散度上界值，则所述待检测数据对应设备的新工作状态，此时在训练数据集中，没有找到与待检测数据特征近似的训练数据，当前的训练数据集没有包含待检测数据的特征与设备工作状态对应关系的信息；可以将待检测数据及待检测数据对应的设备当前的工作状态信息添加到训练数据集中，增加训练数据集对故障诊断的全面性。
进一步地，所述根据所述待检测数据和所述训练数据集，计算得到所述待检测数据与所述训练数据集的每一训练数据之间的JS散度值，并且将获得的JS散度值构成第一JS散度集合，包括：
将所述训练数据集的每一训练数据作为当前输入数据，将当前输入数据输入给基于最优带宽的核密度估计过程，获得最优核密度估计，并将当前获得的最优核密度估计作为所述训练数据集的当前训练数据相应的最优核密度估计；
将所述待检测数据作为当前输入数据，将当前输入数据输入给基于最优带宽的核密度估计过程，获得最优核密度估计，并将当前获得的最优核密度估计作为所述待检测数据的最优核密度估计；以及，
根据所述待检测数据的最优核密度估计和所述训练数据集的各训练数据各自相应的最优核密度估计，计算所述待检测数据与所述训练数据集中的每一训练数据之间的JS散度值，并且由计算得到的JS散度值构成第一JS散度集合。
进一步地，所述利用滑窗原理对所述待选训练数据进行采样，得到滑窗数据集，计算所述滑窗数据集和所述待选训练数据之间的JS散度分布，并使用核密度估计的方法获得JS散度上界值，包括：
根据指定宽度的滑动窗口在所述待选训练数据上滑动选取数据，得到至少一个滑窗训练数据，由所述至少一个滑窗训练数据构成所述滑窗数据集；
将所述滑窗数据集的每一滑窗训练数据作为当前输入数据，将当前输入数据输入给基于最优带宽的核密度估计过程，获得最优核密度估计，并将当前获得的最优核密度估计作为所述滑窗数据集的当前滑窗训练数据相应的最优核密度估计；
将所述待选训练数据作为当前输入数据，将当前输入数据输入给基于最优带宽的核密度估计过程，获得最优核密度估计，并将当前获得的最优核密度估计作为所述待选训练数据相应的最优核密度估计；
根据所述滑窗数据集的各滑窗训练数据各自相应的最优核密度估计和所述待选训练数据相应的最优核密度估计，计算所述滑窗数据集的每一滑窗训练数据与所述待选训练数据之间的JS散度值，由计算得到的JS散度值构成第二JS散度集合；
估计第二JS散度集合对应的JS散度密度函数；以及，
根据所述JS散度密度函数和指定的显著性水平值，获得JS散度上界值。
进一步地，所述将当前输入数据输入给基于最优带宽的核密度估计过程，获得最优核密度估计，包括：
设置带宽的初始值、给定的估计精度和最大迭代次数，并循环执行后续步骤，直到第一跳出条件或者第二跳出条件之一满足时，跳出循环；
根据带宽、核函数K(·)和核密度估计公式
Figure BDA0002934840890000061
计算核密度估计；
根据以下公式计算当前循环中的带宽：
Figure BDA0002934840890000062
判断若满足第一跳出条件，则已经得到最优带宽，跳出循环；
判断若满足第二跳出条件，则迭代次数已经超限，跳出循环；
保留当前循环计算得到的带宽，用于在下次循环中判断第一跳出条件；
其中：h表示带宽；K(·)表示核函数；rj表示当前输入数据中的第j个元素；第一跳出条件为当前循环计算得到的带宽减去上次循环计算得到的带宽，所得的差值的绝对值小于给定的估计精度；第二跳出条件为循环次数达到最大迭代次数。
进一步地，所述根据所述待检测数据的最优核密度估计和所述训练数据集的各训练数据各自相应的最优核密度估计，计算所述待检测数据与所述训练数据集的每一训练数据之间的JS散度值，并且由计算得到的JS散度值构成第一JS散度集合，包括：
依据以下公式，计算所述待检测数据与所述训练数据集的每一训练数据之间的JS散度值:
Figure BDA0002934840890000071
由计算得到的JS散度值构成第一JS散度集合：
{JS(Z,R1),JS(Z,R2),JS(Z,R3),…，JS(Z,Rq)}
其中：Z是待检测数据；Ri是训练数据集{R1,R2,R3,…，Rq}中的训练数据；
Figure BDA0002934840890000072
是训练数据的最优核密度估计；
Figure BDA0002934840890000073
是待检测数据的最优核密度估计。
进一步地，所述根据所述滑窗数据集的各滑窗训练数据各自相应的最优核密度估计和所述待选训练数据相应的最优核密度估计，计算所述滑窗数据集的每一滑窗训练数据与所述待选训练数据之间的JS散度值，由计算得到的JS散度值构成第二JS散度集合，包括:
根据以下公式，计算所述滑窗数据集的每一滑窗训练数据与所述待选训练数据之间的JS散度值:
Figure BDA0002934840890000074
由计算得到的JS散度值构成第二JS散度集合:
{JS1,JS2,JS3,…,JSm-p}
其中：R(j)表示滑窗数据集中的第j个滑窗训练数据；R表示待选训练数据；
Figure BDA0002934840890000075
是待选训练数据对应的最优核密度估计；
Figure BDA0002934840890000076
是第j个滑窗训练数据对应的最优核密度估计；H(·)表示求熵运算。
进一步地，所述估计第二JS散度集合对应的JS散度密度函数，包括：
根据如下公式估计第二JS散度集合对应的JS散度密度函数：
Figure BDA0002934840890000077
其中：JSj是第二JS散度集合中的第j个元素；K(·)是核函数；h是带宽；m是待选训练数据中的元素个数；p是以元素个数为单位的滑动窗口的宽度。
进一步地，所述根据所述JS散度密度函数和指定的显著性水平值，获得JS散度上界值，包括：
通过对如下公式做数值积分得到JS散度上界值：
Figure BDA0002934840890000078
JShigh＝h*i
其中：h是步长；i是指步长编号；
Figure BDA0002934840890000079
是JS散度密度函数；α是指定的显著性水平值；JShigh是JS散度上界值。
进一步地，在所述待选JS散度值大于所述JS散度上界值，则所述待检测数据对应于设备的新工作状态，之后还包括：
定义所述待检测数据对应的设备的新状态模式的状态模式标签，将所述待检测数据加入到所述训练数据集中。
下面结合具体的应用实例对本发明实施例上述技术方案进行详细说明，实施过程中没有介绍到的技术细节，可以参考前文的相关描述。
核密度估计作为一种非参数估计方法，其直接从样本数据出发，对数据的总体分布不做先验假定，更适用于当前的海量数据。随着数据的增长，会带来更多的信息，此时应当对数据的分布假设做成修正，但是随着维度的增长，多维的核密度估计变得复杂，其最优带宽公式并没有给出。一种方法通过对高维数据进行低维表征，而后在不同维度上分别进行核密度估计，从而从一定程度上对高维数据的分布进行了刻画。但该方法仍然没有解决高维数据的最优核密度估计问题。事实上，多维数据的最优核密度估计问题是一个值得深入研究的问题。事实上，核密度估计方法可以克服频率直方图主观离散操作的不足，而交叉熵等方法能对数据分布的差异进行定量度量，是解决该问题的思路之一。另一种方法通过核密度估计的方法对数据分布进行了重构，并构建了交叉熵函数对分布差异进行了度量，从而提高了故障检测结果。但是该方法针对各个维度分别进行核密度估计，无法反映出不同维数之间的关联性，同时交叉熵函数在密度分布的刻画上并不精细，导致故障检测效果降低，尤其是对未包含在故障集的非预期故障检测效果不佳。本发明将核密度估计方法拓展到高维数据上，避免了针对各维度单独进行核密度估计时造成的信息损失，从而更好地刻画数据的密度概率分布。同时，改进了传统方法中使用交叉熵函数作为密度分布差异的度量方法，采用了JS散度作为密度分布差异的度量，规避了采用交叉熵函数作为度量导致的相对性。大多数方法都仅基于距离度量进行故障辨识，但仅依靠距离度量无法有效检测出非预期故障。本发明在JS散度的基础上，利用滑窗原理，导出了样本密度分布和总体密度分布之间JS散度的分布特征，并基于此给出了故障辨识的检测阈值，从而实现对非预期故障的辨识。
观测数据通常可以分解为低频部分和高频部分两个部分，一般来说，低频部分主要表征系统的非平稳工作状态，具有一定的趋势性、单调性和周期性；而高频部分主要表征系统平稳工作状态，具有一定的零均值、高频振动和统计稳定性。对于低频部分而言，可以通过系统状态方程来刻画其变化规律，当低频部分发生故障时，故障引起的征兆变化是相对明显的，针对低频信号的故障检测方法较为完善。对于以设备为代表的高频振动系统，其微小故障往往容易被正常、大幅值、高频振动所掩盖。因此需要对观测数据进行深入的分析。观测数据通常可以分解为本征部分和非本征部分，一般来说，本征部分主要表征系统的主要工作状态，而非本征部分主要表征系统噪声等。对于本征部分，可以通过系统状态方程来刻画其变化规律，当本征部分发生故障时，引起的征兆相对显著，相应的故障检测方法较为成熟。但对于高频振动信号，其微小故障往往隐藏在非本征部分中，容易被噪声等掩盖。因此需要对观测数据进行深入的分析。
1.信号分解
在设备初始运行阶段，系统运行不平稳会造成数据波动较大，这部分数据不仅对系统趋势会产生较大的影响，还会影响数据的统计特性，所以需要对数据进行截断处理，去除不平稳信号。去除不平稳时段数据后的时间序列对应的时间为{t1,t2,…,tm}，得到如下m个观测数据
Y＝[y(t1),y(t2),…,y(tm)] (0.1)
每次采样y(ti)均包含n个特征，分量形式为
y(ti)＝[y1(ti),y2(ti),…,yn(ti)]T,i＝1,2,…,m (0.2)
数据Y可分解为
Figure BDA0002934840890000081
其中
Figure BDA0002934840890000082
表示本征部分，主要是由趋势项构成，R表示非本征部分，主要由噪音及故障数据构成。
本征部分一般由多种信号复合而成，选定合适的基函数f(t)＝[f1(t),f2(t),…,fs(t)]T可以对本征部分进行刻画，通过遍历m个数据对非线性数据Y建模，得到方程组：
Figure BDA0002934840890000083
记
Figure BDA0002934840890000084
则式(0.4)可表示为
Y＝βF (0.6)
从而β的有效估计为
Figure BDA0002934840890000085
利用式(0.3)，(0.7)可以得到信号分解为
Figure BDA0002934840890000091
2.传统检测统计量
简单起见，记ri＝r(ti),i＝1,2,…,m，则由公式(0.8)得到信号分解后的训练数据为R＝[r1,r2,…,rm]，通常认为它是期望为0的正态随机向量，即
ri～N(0,∑) (0.9)
其中Σ是总体协方差矩阵。在协方差矩阵Σ未知时，Σ的无偏估计由下式给出
Figure BDA0002934840890000092
设Z＝[z1,z2,…,zp]是待检测测试窗内数据，其样本均值为
Figure BDA0002934840890000093
若Z与训练数据R均来自相同的模式，则
Figure BDA00029348408900000913
仍然服从正态分布且
Figure BDA0002934840890000094
可以构造T2统计量
Figure BDA0002934840890000095
T2统计量的分布满足
Figure BDA0002934840890000096
故在给定显著性水平为α的情况下，若
Figure BDA0002934840890000097
则认为测试窗内数据Z与训练数据R均来自相同的模式，否则认为不同。这个判据的误判率为α。
3.最优核密度估计
传统的故障检测方法，主要包含了信号分解技术和基于T2统计量的模式判别方法。但基于T2统计量的检测方法假设的数据满足正态分布，而现实观测数据可能并不满足该假设，导致T2统计量的判别性能不能满足设计要求。另外，可以发现T2统计量主要从本征项
Figure BDA0002934840890000098
和协方差矩阵
Figure BDA0002934840890000099
的角度对数据进行了检验，这两个属性不足以刻画系统的所有统计特性，当微小故障被数据噪声淹没时，容易漏检。对此，本节构建了多维数据的核密度估计方法，使之能更精确地描述数据的概率统计特性。
3.1最优带宽定理
对于观测数据而言，通常可以用频率直方图直观表现其统计特征，但在实际应用过程中，由于频率直方图是离散的统计方法，直方图的区间数量不好划分，更重要的是离散化操作会给后续进一步数据处理带来不便。为了可以克服这些局限性，提出了核密度估计方法，该方法属于非参数估计方法，它通过采样数据直接对总体的概率密度分布进行估计。
对于任意一个点
Figure BDA00029348408900000910
假定某个模式下的概率密度为f(x)，依据预处理后的采样数据R＝[r1,r2,…,rm]对f(x)进行核密度估计，估计公式如下
Figure BDA00029348408900000911
其中，m是采样数据个数，n是采样数据维数，K(·)是核函数，hm是带宽。为后续讨论方便，在不产生疑义的情况下，记
Figure BDA00029348408900000912
一般而言，因为核函数满足
Figure BDA0002934840890000101
所以
Figure BDA0002934840890000102
从而有
Figure BDA0002934840890000103
这意味
Figure BDA0002934840890000104
同时满足正定性、有连续性和规范性，故用
Figure BDA0002934840890000105
作为f(x)的核密度估计是合理的。通常从下表选取核函数：
表0常用核函数表
Figure BDA0002934840890000106
本文用积分均方误差(MISE，mean integral square error)刻画核密度估计的性能，如下
Figure BDA0002934840890000107
可以认为MISE
Figure BDA0002934840890000108
越小，表明
Figure BDA0002934840890000109
对f(x)的估计越有效。
核函数的类型K(·)和带宽hm都对MISE
Figure BDA00029348408900001010
有影响。一方面，MISE
Figure BDA00029348408900001011
对核函数K(·)的选择不敏感，也就是选用不同核函数得到估计结果的积分均方误差几乎一致，这一点在后续的推导过程中也有所体现。另一方面，MISE
Figure BDA00029348408900001012
主要取决于带宽hm的选取，如果hm选的太小，密度估计值
Figure BDA00029348408900001013
会因为随机性增强而呈现不规则形状，而当hm选的太大，密度估计值
Figure BDA00029348408900001014
会过度平均化而无法展示足够多的细节。
本发明以定理的形式给出了最优带宽公式，这也是本发明主要的理论结果之一，如下.
定理0.1对于任意n维概率密度函数f(·)和表中任意一种核函数K(·)，若用(0.16)中的
Figure BDA00029348408900001015
估计f(·)，并且海塞矩阵的迹是可积的，即
Figure BDA00029348408900001016
存在，则当积分均方误差MISE
Figure BDA00029348408900001017
取最小值时，带宽hm满足
Figure BDA00029348408900001018
其中cK和dK是两个常值，如下
Figure BDA00029348408900001019
称公式(0.19)为最优带宽公式，对应地hm为最优带宽。
为了证明该定理成立，先给出两个等式
Figure BDA00029348408900001020
实际上，
Figure BDA00029348408900001021
而且
Figure BDA0002934840890000111
由公式(0.21)第一式得
Figure BDA0002934840890000112
由公式(0.21)两式得
Figure BDA0002934840890000113
通过公式(0.24)-(0.25)得
Figure BDA0002934840890000114
为后续推理方便，先给出如下定理。
定理0.2对于任意矩阵
Figure BDA0002934840890000115
K(·)是表中的任意一种核函数，则有
Figure BDA0002934840890000116
证明：若奇函数g(x)在
Figure BDA0002934840890000121
上可积，则必有
Figure BDA0002934840890000122
类似的可以验证，对于表中的任意一种核函数满足
Figure BDA0002934840890000123
故
Figure BDA0002934840890000124
从而定理2得证。
对于任意单位长向量
Figure BDA0002934840890000125
由Taylor展式得
Figure BDA0002934840890000126
如果带宽hm满足条件
Figure BDA0002934840890000127
则由公式(0.26)，(0.29)，(0.30)，(0.31)得：
Figure BDA0002934840890000128
实际上
Figure BDA0002934840890000129
若
Figure BDA0002934840890000131
可积，则有
Figure BDA0002934840890000132
当MISE
Figure BDA0002934840890000133
最小时，公式(0.34)关于hm求导为0，即
Figure BDA0002934840890000134
从而解得定理1中的最优带宽
Figure BDA0002934840890000135
综上。从公式(0.21)～(0.36)可知定理1得证。
备注给定样本序列R＝[r1,r2,…,rm]，通过公式(0.36)可以选取合适的带宽hm，从而在依据公式(0.16)给出核密度估计函数。影响核函数带宽hm选取的因素主要包括cK和dK，而这两者对核函数K(·)的选择不敏感，对最终的带宽hm选择几乎没有影响。
3.2最优带宽算法
公式(0.36)给出了最优带宽公式，但是公式(0.36)中f(x)是未知的，所以
Figure BDA0002934840890000136
也是未知的，此时可以用公式(0.16)中的
Figure BDA0002934840890000137
代替f(x)，可算得带宽参数hm的一个近似值。进一步可以用迭代算法算得一个更精确的带宽参数，下述定理表明该算法是收敛的，该定理是本文的另一个主要理论结果。
定理0.3对于任意n维概率密度函数f(·)和高斯核函数K(·)，用公式(0.16)中的
Figure BDA0002934840890000138
估计f(·)，用如下公式迭代计算hm
Figure BDA0002934840890000139
则该是收敛
Figure BDA00029348408900001310
是收敛的，记
Figure BDA00029348408900001311
证明：
对于高斯核函数K(u)，即
Figure BDA00029348408900001312
可知dK为自由度为n的卡方分布，期望等于自由度，故
dK＝∫uTuK(u)du＝n (0.39)
另外
Figure BDA00029348408900001313
将(0.38)(0.39)(0.40)代入(0.36)，并用(0.16)中的
Figure BDA00029348408900001314
代替f(x)得到计算hm的迭代形式
Figure BDA0002934840890000141
为了后续推理方便，先给出如下引理。
命题0.1对任意函数f1,f2,…,fn，下述不等式
Figure BDA0002934840890000142
成立，当且仅当f1(x)＝f2(x)＝…＝fn(x)几乎处处成立时等号成立。
实际上，对于任意的
Figure BDA0002934840890000143
有
0≤(f1(x)+f2(x)+…+fn(x))2≤n(f1(x)2+f2(x)2+…+fn(x)2) (0.43)
对(0.43)两边进行积分得到
Figure BDA0002934840890000144
显然，式(0.43)中等号成立的条件为：f1(x)＝f2(x)＝…＝fn(x)几乎处处成立。
进一步，公式(0.38)关于变量xi的二阶导为
Figure BDA0002934840890000145
所以
Figure BDA0002934840890000146
由引理和公式(0.46)得
Figure BDA0002934840890000147
值得注意的是，当hm,k足够大时，可以认为
Figure BDA0002934840890000148
几乎处处相同，此时可以认为公式(0.47)中的等号成立，此即
Figure BDA0002934840890000149
从而当hm,k较大时，该迭代过程是递减的。又由于hm,k是有下界的，从而该算法收敛。
综上基于本文主要理论结果，即定理1和定理3，给出如下基于最优带宽的核密度估计算法，对应的流程图如图2所示。
Figure BDA0002934840890000151
4.基于JS散度分布的故障检测方法
基于最优带宽的多维核密度估计方法，该方法能精确地描述多维数据的密度分布。并在此基础上利用JS散度度量分布差异，使其能更能凸显不同模式数据的统计特性差异。
4.1模式差异指标
利用核函数的方法得到了高维数据的概率密度估计，推导了最优的带宽公式，给出了最优算法的迭代计算方法，并证明了迭代的收敛性。当系统发生故障时，系统的状态必然发生改变，系统输出的统计特征也随之改变，导致观测数据的密度分布产生变化。对于两组样本窗口数据R,Z，可以用交叉熵H(R,Z)衡量R和Z的分布差异：
Figure BDA0002934840890000152
其中
Figure BDA0002934840890000153
分别表示由公式(0.16)得到的关于R和Z的最优核函数估计值。
值得注意的是，H(R,Z)并不满足距离的定义，因为H(R,Z)不一定满足正定性和对称性，也就是说可能H(R,Z)<0或者H(R,Z)≠H(Z,R)。
(1)R和Z的分布差异越小，H(R,Z)越小，这意味着即使H(R,Z)<0，用H(R,Z)衡量R和Z的分布差异也是合理的。
(2)但是，分布差异定量描述必须满足对称性，否则，交换位置，分布差异就不同，这些难以接受的。用JS散度JS(R,Z)作为R和Z的分布差异度量，如下
Figure BDA0002934840890000154
此时易得
Figure BDA0002934840890000155
本文利用公式(0.50)衡量测试数据Z和训练数据R的分布差异大小，从而实现故障检测与隔离。
4.2模式判别方法
若训练数据有q个模式{R1,R2,…,Rq}即训练数据集，可以算的测试数据Z(即待检测数据)与不同模式Ri(即训练数据集中的第i个训练数据)间的JS散度集合{JS(Z,R1),JS(Z,R2),…,JS(Z,Rq)}(即第一JS散度集合)，其中
Figure BDA0002934840890000156
若i0是最小JS散度(即第一JS散度集合中的最小JS散度值，即待选JS散度值)对应的模式标签(即状态模式标签)，即
i0＝argmin{JS(Z,R1),JS(Z,R2),…,JS(Z,Rq)} (0.53)
则有道理认为测试数据Z与训练数据
Figure BDA00029348408900001613
(即待选训练数据)属于同一个故障模式(即设备的一种工作状态)。
但是，应用中可能未知的新故障模式，此时公式(0.53)必然会把测试数据Z判断为第i0已知故障模式，这样显然是不合理的。
如果
Figure BDA00029348408900001614
过大，我们有理由认为测试数据Z来至未知的新故障模式，其标签记为q+1(即设备的新状态模式的状态模式标签)，但是如何给出
Figure BDA00029348408900001615
阈值JShigh呢？下面给出一种确定JShigh(即JS散度上界值)的方法。
对于第i0个模式的训练数据
Figure BDA00029348408900001616
利用公式(0.16)可以得到数据集的密度估计(即待选训练数据对应的最优核密度估计)如下
Figure BDA0002934840890000161
另外，固定采样窗口长度为p(p<m)，通过滑动采样窗口得到新的采样数据为
Figure BDA00029348408900001612
对于每一个R(j)，利用公式(0.16)可以得到数据集的密度估计(即滑窗数据集中的第j个滑窗训练数据对应的最优核密度估计)为
Figure BDA0002934840890000162
利用公式(0.50)可以得到样本数据R(j)和训练数据集R之间的散度(即滑窗数据集的每一滑窗训练数据与所述待选训练数据之间的JS散度值)为
Figure BDA0002934840890000163
利用公式(0.56)可以得到一系列的JS散度计算值集合JS＝{JS1,JS2,…,JSm-p}(即第二JS散度集合)，我们用该集合给出JS散度的密度函数fJS(x)的估计式
Figure BDA0002934840890000164
(即第二JS散度集合对应的JS散度密度函数)，如下
Figure BDA0002934840890000165
若取显著性水平为α，则
Figure BDA0002934840890000166
超过阈值JShigh的概率为
Figure BDA0002934840890000167
因为JS散度的分布类型不是常用随机分布类型的，所以分位数不能通过查表获得，只能通过数值积分获得，若h是步长，且
Figure BDA0002934840890000168
则有道理认为
JShigh＝h*i (0.60)
由公式(0.60)构建如下故障检测和隔离的准则。
准则1对于新的待检测数据Z＝[z1,z2,…,zp]，若i0是最小JS散度对应的模式标签，见公式(0.53)，第i0个模式的训练数据
Figure BDA0002934840890000169
JS散度上界为JShigh，见公式(0.60)，且满足
Figure BDA00029348408900001610
则人认为测试数据Z与训练数据
Figure BDA00029348408900001611
属于同一个故障模式(即待检测数据与待选训练数据对应相同的状态模式标签)，否则认为测试数据Z来自未知的新故障模式(即待检测数据对应设备的新工作状态)，其标签记为q+1(即定义待检测数据对应的设备的新状态模式的状态模式标签，将待检测数据加入到训练数据集中)。
综上，给出基于最优带宽的故障诊断方法，对应的故障诊断方法流程图如图3所示。
Figure BDA0002934840890000171
备注值得注意的是，公式(0.55)，(0.56)表明，JS散度的计算结果与采样数据长度是直接相关的，显而易见，随着采样数据长度的增加，利用公式(0.55)得到的密度估计也越能刻画样本的分布特性，从而能显著提高故障检测的精度。
5.下面通过将本方法具体应用于轴承数据，来说明本方法的应用过程和效果。
采用凯斯西储大学轴承数据中心的轴承数据作为故障诊断数据。轴承数据为电机负载为0马力时的运行数据，采样频率为12kHz。数据集包含四组样本数据：正常数据(f0)，0.007英寸内滚道故障数据(f1)，0.014英寸内滚道故障数据(f2)，0.014英寸外滚道故障数据(f3)。每组数据具有两个维度：驱动端加速度数据(fi-DE)，风扇端加速度数据(fi-FE)。
5.1数据预处理
轴承运行过程中的观测数据一般表现出存在较明显的周期性，为此需要对这种本征信号进行剔除。以正常数据f0为例，对观测信号进行快速傅立叶分析可以得到信号中的主要频率，f0的傅立叶频谱图如图4驱动端加速度正常数据功率谱分布图和图5风扇端加速度正常数据功率谱分布图所示。
从图4和图5中可以发现，其主要频率在Fs＝1036Hz，从而构造基函数为
f(t)＝[1 sin(1036×2πt) cos(1036×2πt)]T
利用公式(0.7)计算得到β的估计
Figure BDA0002934840890000172
为
Figure BDA0002934840890000173
从而得到f0预处理后数据如图6和图7所示的预处理前后数据对比，其中图6表示f0-DE预处理后的数据，图7表示f0-FE预处理后的数据。
在后文的故障检测过程中，所有的数据fi均类似于f0进行了如上的操作，其结果仍记为fi。
5.2预期故障检测效果
选取f0,f1,f2中前20480个样本点作为训练集，分别记为f0-train,f1-ttain,f2-train；后81920个样本点作为测试集，分别记为f0-test,f1-test,f2-test。每次检测使用128个样本点作为检测对象。训练集数据如图8到图13所示，其中图8，图9分别表示f0-train两个维度的数据f0-train-DE,f0-train-FE，图10和图11分别表示f1-train两个维度的数据f1-train-DE,f1-train-FE，图12和图13分别表示f2-train两个维度的数据f2-train-DE,f2-train-FE。
图8到图13表明：轴承数据大多为高频数据，同时轴承故障不改变观测的均值，但是改变了数据的散布特征或者数据之间的关联性，如前文所说，轴承数据的这些特点使得轴承故障检测极具挑战。
利用算法得到带宽为hm＝0.0445
通过公式(0.16)得到训练集的核密度估计结果如图14到图19所示，其中图14、图15和图16分别表示训练集f0-train,f1-train,f2-train的二维频率直方图，图17、图18和图19分别表示训练集f0-train,f1-train,f2-train的二维核密度估计结果。
图14到图19进一步表明，轴承故障主要是改变了数据的散布特征和数据之间的关联性。同时，图14到图19表明通过公式(0.16)得到的训练集核密度估计结果和训练集本身的数据分布吻合较好，从而此方法确实可以对高维数据的分布进行刻画。
利用公式(0.56)，(0.59)得出训练集数据中的f0,f1,f2的JS散度以及分布的核密度估计分别如图20到图22所示。利用公式(0.58)计算得到显著性水平为α＝0.05时，训练集的检测阈值分别为：
Figure BDA0002934840890000181
从而，在测试集上的检测结果如图23到图28所示，图中·代表检测正确，*代表检测错误；若检测点落在上下阈值之间，则待检测数据为正常；否则数据为故障。图23到图25分别是f0,f1,f2的采用交叉熵函数的检测结果；图26到图28分别是f0,f1,f2的本发明的方法的检测结果；
进一步，可以得到不同方法下的测试结果如下表所示：
表1不同方法下的检测准确率
方法\类型 T<sup>2</sup>检测量 交叉熵函数 JS散度
正常数据 95.80％ 96.95％ 97.03％
0.007英寸内滚道故障数据 83.47％ 94.41％ 95.81％
0.014英寸内滚道故障数据 78.11％ 94.19％ 95.36％
表2表明，本发明基于多维核密度估计和JS散度构建的轴承故障辨识在训练数据集上取得了比传统的T2检测方法更优的结果，对正常数据的检测率从95.08％提升到97.03％，对0.007英寸内滚道故障数据的检测率从81.33％提升到95.81％，对0.014英寸内滚道故障数据的检测率从70.69％提升到95.36％。同时，该方法相比较于交叉熵函数也有一定程度的提升，对正常数据的检测率从96.95％提升到97.03％，对0.007英寸内滚道故障数据的检测率从94.41％提升到95.81％，对0.014英寸内滚道故障数据的检测率从94.19％提升到95.36％。
5.3非预期故障检测效果
实际上，训练集并不一定包含所有的故障类型，而对非预期故障的检测一直是一个难题。对此，本节将f3作为非预期故障进行故障检测，值得注意的是训练集样本中并未包含任何f3的信息。非预期故障f3数据如图29和图30所示，其中图29表示其驱动端加速度数据(f3-DE)，图30表示其风扇端加速度数据(f3-FE)。图29和图30表明非预期故障f3的数据其他两类故障数据较为接近，如果故障检测方法不敏感则会导致检测率显著降低。利用不同方法对0.014英寸外滚道故障数据的检测结果如下表所示：
表2不同方法对非预期故障(0.014英寸外滚道故障数据)的检测结果
方法\类型 T<sup>2</sup>检测量 交叉熵函数 JS散度
0.014英寸外滚道故障数据 41.55％ 53.16％ 69.49％
表3表明，传统的T2检测方法对非预期故障的检测率较低，仅为41.55％，而利用交叉熵函数作为度量的方法仅能对非预期故障的检测率为53.16％，其效果并不明显。而本文构建的JS散度方法则较准确地辨识出了非预期故障，其检测率达到69.49％。这是由于JS散度在衡量分布间差异时更为精确。
5.4窗宽对故障诊断效果的影响
故障诊断效果与数据窗宽相关，故本节考察了在不同窗宽下的故障诊断效果。其结果如图31所示。
从图31可以发现，随着检测窗口的增加，本发明提出的方法对预期故障检测的检测率先上升随后表现趋于平稳，这是由于当检测窗口长度增加到一定程度后，待检测数据中就已经包含足够多的信息，此时检测窗口再继续增加对故障检测率的提升贡献率并不大。同时可以发现，对于非预期故障而言，其检测率随着检测窗口的变长而迅速增加，这是由于检测窗口越长，待检测数据中就包含更多的信息，也就更能表征出和已知故障之间的差异。
6.结论
本发明通过多维核密度函数估计和JS散度构建了一种设备故障检测和辨识方法，并应用本方法通过对轴承的运行数据和运行状态进行测量及分析，验证了本方法的有效性。本文通过滑动采样窗口的方法，导出了样本密度分布和总体密度分布之间JS散度的分布特征，并基于此给出了故障检测的阈值，从而实现对不同故障的辨识。理论表明，多维核密度估计方法能减少对各维度进行处理时造成的信息损失，而采用JS散度度量密度分布差异则比传统的交叉熵函数更为精确。实验验证了上述结论：其一，对于预期故障，该方法的检测效果明显优于传统方法，较交叉熵函数也有一定程度的提升。其二，对于非预期故障，由于传统的方法对分布差异的度量不够精细而无法有效检出，而本发明的方法则对非预期故障的检测效果有显著提升。
进一步，由于该方法依赖带宽，随着带宽的增加，检测效果先上升后降低。另外，随着检测窗口的增长，检测效果一直提升。故本文在给定窗宽的条件下，给出了多维核密度函数最优带宽的估计公式。实验表明，该公式对任意数据条件都适用，因而具有一定的普适性。
应该明白，公开的过程中的步骤的特定顺序或层次是示例性方法的实例。
基于设计偏好，应该理解，过程中的步骤的特定顺序或层次可以在不脱离本公开的保护范围的情况下得到重新安排。所附的方法权利要求以示例性的顺序给出了各种步骤的要素，并且不是要限于所述的特定顺序或层次。
在上述的详细描述中，各种特征一起组合在单个的实施方案中，以简化本公开。不应该将这种公开方法解释为反映了这样的意图，即，所要求保护的主题的实施方案需要比清楚地在每个权利要求中所陈述的特征更多的特征。相反，如所附的权利要求书所反映的那样，本发明处于比所公开的单个实施方案的全部特征少的状态。因此，所附的权利要求书特此清楚地被并入详细描述中，其中每项权利要求独自作为本发明单独的优选实施方案。
为使本领域内的任何技术人员能够实现或者使用本发明，上面对所公开实施例进行了描述。对于本领域技术人员来说；这些实施例的各种修改方式都是显而易见的，并且本文定义的一般原理也可以在不脱离本公开的精神和保护范围的基础上适用于其它实施例。因此，本公开并不限于本文给出的实施例，而是与本申请公开的原理和新颖性特征的最广范围相一致。
上文的描述包括一个或多个实施例的举例。当然，为了描述上述实施例而描述部件或方法的所有可能的结合是不可能的，但是本领域普通技术人员应该认识到，各个实施例可以做进一步的组合和排列。因此，本文中描述的实施例旨在涵盖落入所附权利要求书的保护范围内的所有这样的改变、修改和变型。此外，就说明书或权利要求书中使用的术语“包含”，该词的涵盖方式类似于术语“包括”，就如同“包括，”在权利要求中用作衔接词所解释的那样。此外，使用在权利要求书的说明书中的任何一个术语“或者”是要表示“非排它性的或者”。
本领域技术人员还可以了解到本发明实施例列出的各种说明性逻辑块(illustrative logical block)，单元，和步骤可以通过电子硬件、电脑软件，或两者的结合进行实现。为清楚展示硬件和软件的可替换性(interchangeability)，上述的各种说明性部件(illustrative components)，单元和步骤已经通用地描述了它们的功能。这样的功能是通过硬件还是软件来实现取决于特定的应用和整个系统的设计要求。本领域技术人员可以对于每种特定的应用，可以使用各种方法实现所述的功能，但这种实现不应被理解为超出本发明实施例保护的范围。
以上所述的具体实施方式，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施方式而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (9)

1.一种基于优化核密度估计和JS散度的故障诊断方法，其特征在于，包括：
通过传感器采集设备工作期间的运行数据，将采集到的运行数据作为待检测数据；
根据所述待检测数据和训练数据集，计算得到所述待检测数据与所述训练数据集中的每一训练数据之间的JS散度值，并且将获得的JS散度值构成第一JS散度集合；其中，所述训练数据集中的各训练数据是由传感器采集的设备工作期间的运行数据；所述训练数据集中的各训练数据与设备的已知的各状态模式标签相对应；所述状态模式标签用于标识设备的工作状态；
将所述第一JS散度集合中最小JS散度值对应的所述训练数据集中的训练数据，作为待选训练数据，并将所述最小JS散度值作为待选JS散度值；
利用滑窗原理对所述待选训练数据进行采样，得到滑窗数据集，计算所述滑窗数据集和所述待选训练数据之间的JS散度分布，并使用核密度估计的方法获得JS散度上界值；所述JS散度上界值用于作为设备故障诊断的检测阈值；以及，
根据所述待选JS散度值和所述JS散度上界值，确定所述待检测数据对应的设备的状态模式标签，具体包括：
若所述待选JS散度值小于或者等于所述JS散度上界值，则所述待检测数据对应的状态模式标签与所述待选训练数据对应的状态模式标签相同，或者，
若所述待选JS散度值状态大于所述JS散度上界值，则所述待检测数据对应设备的新工作状态。
2.如权利要求1所述的基于优化核密度估计和JS散度的故障诊断方法，其特征在于，所述根据所述待检测数据和所述训练数据集，计算得到所述待检测数据与所述训练数据集的每一训练数据之间的JS散度值，并且将获得的JS散度值构成第一JS散度集合，包括：
将所述训练数据集的每一训练数据作为当前输入数据，将当前输入数据输入给基于最优带宽的核密度估计过程，获得最优核密度估计，并将当前获得的最优核密度估计作为所述训练数据集的当前训练数据相应的最优核密度估计；
将所述待检测数据作为当前输入数据，将当前输入数据输入给基于最优带宽的核密度估计过程，获得最优核密度估计，并将当前获得的最优核密度估计作为所述待检测数据的最优核密度估计；以及，
根据所述待检测数据的最优核密度估计和所述训练数据集的各训练数据各自相应的最优核密度估计，计算所述待检测数据与所述训练数据集中的每一训练数据之间的JS散度值，并且由计算得到的JS散度值构成第一JS散度集合。
3.如权利要求1所述的基于优化核密度估计和JS散度的故障诊断方法，其特征在于，所述利用滑窗原理对所述待选训练数据进行采样，得到滑窗数据集，计算所述滑窗数据集和所述待选训练数据之间的JS散度分布，并使用核密度估计的方法获得JS散度上界值，包括：
根据指定宽度的滑动窗口在所述待选训练数据上滑动选取数据，得到至少一个滑窗训练数据，由所述至少一个滑窗训练数据构成所述滑窗数据集；
将所述滑窗数据集的每一滑窗训练数据作为当前输入数据，将当前输入数据输入给基于最优带宽的核密度估计过程，获得最优核密度估计，并将当前获得的最优核密度估计作为所述滑窗数据集的当前滑窗训练数据相应的最优核密度估计；
将所述待选训练数据作为当前输入数据，将当前输入数据输入给基于最优带宽的核密度估计过程，获得最优核密度估计，并将当前获得的最优核密度估计作为所述待选训练数据相应的最优核密度估计；
根据所述滑窗数据集的各滑窗训练数据各自相应的最优核密度估计和所述待选训练数据相应的最优核密度估计，计算所述滑窗数据集的每一滑窗训练数据与所述待选训练数据之间的JS散度值，由计算得到的JS散度值构成第二JS散度集合；
估计第二JS散度集合对应的JS散度密度函数；以及，
根据所述JS散度密度函数和指定的显著性水平值，获得JS散度上界值。
4.如权利要求2或3所述的基于优化核密度估计和JS散度的故障诊断方法，其特征在于，所述将当前输入数据输入给基于最优带宽的核密度估计过程，获得最优核密度估计，包括：
设置带宽的初始值、给定的估计精度和最大迭代次数，并循环执行后续步骤，直到第一跳出条件或者第二跳出条件之一满足时，跳出循环；
根据带宽、核函数K(·)和核密度估计公式
Figure FDA0002934840880000021
计算核密度估计；
根据以下公式计算当前循环中的带宽：
Figure FDA0002934840880000022
判断若满足第一跳出条件，则已经得到最优带宽，跳出循环；
判断若满足第二跳出条件，则迭代次数已经超限，跳出循环；
保留当前循环计算得到的带宽，用于在下次循环中判断第一跳出条件；
其中：h表示带宽；K(·)表示核函数；rj表示当前输入数据中的第j个元素；第一跳出条件为当前循环计算得到的带宽减去上次循环计算得到的带宽，所得的差值的绝对值小于给定的估计精度；第二跳出条件为循环次数达到最大迭代次数。
5.如权利要求2所述的基于优化核密度估计和JS散度的故障诊断方法，其特征在于，所述根据所述待检测数据的最优核密度估计和所述训练数据集的各训练数据各自相应的最优核密度估计，计算所述待检测数据与所述训练数据集的每一训练数据之间的JS散度值，并且由计算得到的JS散度值构成第一JS散度集合，包括：
依据以下公式，计算所述待检测数据与所述训练数据集的每一训练数据之间的JS散度值:
Figure FDA0002934840880000023
由计算得到的JS散度值构成第一JS散度集合：
{JS(Z,R1),JS(Z,R2),JS(Z,R3),…，JS(Z,Rq)}
其中：Z是待检测数据；Ri是训练数据集{R1,R2,R3,…，Rq}中的训练数据；
Figure FDA0002934840880000024
是训练数据的最优核密度估计；
Figure FDA0002934840880000025
是待检测数据的最优核密度估计。
6.如权利要求3所述的基于优化核密度估计和JS散度的故障诊断方法，其特征在于，所述根据所述滑窗数据集的各滑窗训练数据各自相应的最优核密度估计和所述待选训练数据相应的最优核密度估计，计算所述滑窗数据集的每一滑窗训练数据与所述待选训练数据之间的JS散度值，由计算得到的JS散度值构成第二JS散度集合，包括:
根据以下公式，计算所述滑窗数据集的每一滑窗训练数据与所述待选训练数据之间的JS散度值:
Figure FDA0002934840880000026
由计算得到的JS散度值构成第二JS散度集合:
{JS1,JS2,JS3,…,JSm-p}
其中：R(j)表示滑窗数据集中的第j个滑窗训练数据；R表示待选训练数据；
Figure FDA0002934840880000027
是待选训练数据对应的最优核密度估计；
Figure FDA0002934840880000028
是第j个滑窗训练数据对应的最优核密度估计；H(·)表示求熵运算。
7.如权利要求3所述的基于优化核密度估计和JS散度的故障诊断方法，其特征在于，所述估计第二JS散度集合对应的JS散度密度函数，包括：
根据如下公式估计第二JS散度集合对应的JS散度密度函数：
Figure FDA0002934840880000031
其中：JSj是第二JS散度集合中的第j个元素；K(·)是核函数；h是带宽；m是待选训练数据中的元素个数；p是以元素个数为单位的滑动窗口的宽度。
8.如权利要求3所述的基于优化核密度估计和JS散度的故障诊断方法，其特征在于，所述根据所述JS散度密度函数和指定的显著性水平值，获得JS散度上界值，包括：
通过对如下公式做数值积分得到JS散度上界值：
Figure FDA0002934840880000032
JShigh＝h*i
其中：h是步长；i是指步长编号；
Figure FDA0002934840880000033
是JS散度密度函数；α是指定的显著性水平值；JShigh是JS散度上界值。
9.如权利要求1所述的基于优化核密度估计和JS散度的故障诊断方法，其特征在于，在所述待选JS散度值大于所述JS散度上界值，则所述待检测数据对应于设备的新工作状态，之后还包括：
定义所述待检测数据对应的设备的新状态模式的状态模式标签，将所述待检测数据加入到所述训练数据集中。
CN202110158768.1A 2021-02-04 2021-02-04 一种基于优化核密度估计和js散度的故障诊断方法 Active CN113051092B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110158768.1A CN113051092B (zh) 2021-02-04 2021-02-04 一种基于优化核密度估计和js散度的故障诊断方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110158768.1A CN113051092B (zh) 2021-02-04 2021-02-04 一种基于优化核密度估计和js散度的故障诊断方法

Publications (2)

Publication Number Publication Date
CN113051092A true CN113051092A (zh) 2021-06-29
CN113051092B CN113051092B (zh) 2022-05-17

Family

ID=76508916

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110158768.1A Active CN113051092B (zh) 2021-02-04 2021-02-04 一种基于优化核密度估计和js散度的故障诊断方法

Country Status (1)

Country Link
CN (1) CN113051092B (zh)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114112390A (zh) * 2021-11-23 2022-03-01 哈尔滨工程大学 一种非线性复杂系统早期故障诊断方法
CN114995338A (zh) * 2022-05-30 2022-09-02 保控(南通)物联科技有限公司 一种基于规范变量分析与js散度融合的工业过程微小故障检测方法
CN115225528A (zh) * 2022-06-10 2022-10-21 中国科学院计算技术研究所 基于张量填充的网络流量数据分布式测量调度方法及系统

Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103868692A (zh) * 2014-03-18 2014-06-18 电子科技大学 基于核密度估计和k-l散度的旋转机械故障诊断方法
CN105241680A (zh) * 2015-08-26 2016-01-13 电子科技大学 一种基于概率密度函数的旋转机械健康状态评估方法
US20180257663A1 (en) * 2017-03-07 2018-09-13 GM Global Technology Operations LLC Method and apparatus for monitoring an on-vehicle controller
CN108896308A (zh) * 2018-07-02 2018-11-27 昆明理工大学 一种基于概率包络的轮对轴承故障诊断方法
CN109447187A (zh) * 2018-12-25 2019-03-08 中南大学 电机故障诊断方法及系统
CN109828168A (zh) * 2019-01-31 2019-05-31 福州大学 基于核密度估计的变换器故障诊断方法
CN111783286A (zh) * 2020-06-16 2020-10-16 中国人民解放军国防科技大学 基于故障趋势比及特征选择的微小故障诊断方法及系统
CN111812080A (zh) * 2020-06-24 2020-10-23 中国人民解放军国防科技大学 一种基于迭代最优幂的非线性定标精度分析方法
US20200371512A1 (en) * 2018-04-09 2020-11-26 Diveplane Corporation Anomalous Data Detection in Computer Based Reasoning and Artificial Intelligence Systems

Patent Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103868692A (zh) * 2014-03-18 2014-06-18 电子科技大学 基于核密度估计和k-l散度的旋转机械故障诊断方法
CN105241680A (zh) * 2015-08-26 2016-01-13 电子科技大学 一种基于概率密度函数的旋转机械健康状态评估方法
US20180257663A1 (en) * 2017-03-07 2018-09-13 GM Global Technology Operations LLC Method and apparatus for monitoring an on-vehicle controller
US20200371512A1 (en) * 2018-04-09 2020-11-26 Diveplane Corporation Anomalous Data Detection in Computer Based Reasoning and Artificial Intelligence Systems
CN108896308A (zh) * 2018-07-02 2018-11-27 昆明理工大学 一种基于概率包络的轮对轴承故障诊断方法
CN109447187A (zh) * 2018-12-25 2019-03-08 中南大学 电机故障诊断方法及系统
CN109828168A (zh) * 2019-01-31 2019-05-31 福州大学 基于核密度估计的变换器故障诊断方法
CN111783286A (zh) * 2020-06-16 2020-10-16 中国人民解放军国防科技大学 基于故障趋势比及特征选择的微小故障诊断方法及系统
CN111812080A (zh) * 2020-06-24 2020-10-23 中国人民解放军国防科技大学 一种基于迭代最优幂的非线性定标精度分析方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
WENDE TIAN等: "Fault monitoring based on mutual information feature engineering modeling in chemical process", 《CHINESE JOURNAL OF CHEMICAL ENGINEERING》 *
吴海燕等: "ELMD及核密度估计的滚动轴承故障诊断", 《机械强度》 *
王炯琦等: "基于最优估计的数据融合理论", 《应用数学》 *

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114112390A (zh) * 2021-11-23 2022-03-01 哈尔滨工程大学 一种非线性复杂系统早期故障诊断方法
CN114112390B (zh) * 2021-11-23 2024-03-22 哈尔滨工程大学 一种非线性复杂系统早期故障诊断方法
CN114995338A (zh) * 2022-05-30 2022-09-02 保控(南通)物联科技有限公司 一种基于规范变量分析与js散度融合的工业过程微小故障检测方法
CN115225528A (zh) * 2022-06-10 2022-10-21 中国科学院计算技术研究所 基于张量填充的网络流量数据分布式测量调度方法及系统
CN115225528B (zh) * 2022-06-10 2024-04-09 中国科学院计算技术研究所 网络流量数据分布式测量调度方法、系统和介质

Also Published As

Publication number Publication date
CN113051092B (zh) 2022-05-17

Similar Documents

Publication Publication Date Title
CN113051092B (zh) 一种基于优化核密度估计和js散度的故障诊断方法
Yu et al. A new statistical modeling and detection method for rolling element bearing faults based on alpha–stable distribution
Wen et al. Graph modeling of singular values for early fault detection and diagnosis of rolling element bearings
CN111562108A (zh) 一种基于cnn和fcmc的滚动轴承智能故障诊断方法
US20230316051A1 (en) Pre-alarming method for rotary stall of compressors based on temporal dilated convolutional neural network
CN111238843B (zh) 一种基于快速谱峭度分析的风机健康评价方法
CN112633098B (zh) 一种旋转机械故障诊断方法、系统及存储介质
CN105181336B (zh) 一种用于轴承故障诊断的特征选取方法
CN103868692A (zh) 基于核密度估计和k-l散度的旋转机械故障诊断方法
CN108956111B (zh) 一种机械部件的异常状态检测方法及检测系统
CN111881594B (zh) 一种核动力设备的非平稳信号状态监测方法及系统
CN112380992B (zh) 一种加工过程监控数据准确性评估与优化方法及装置
CN116756595A (zh) 一种导电滑环故障数据采集监测方法
CN109612726A (zh) 一种用于振动信号特征提取的多重超阶分析方法
CN114118219A (zh) 基于数据驱动的长期加电设备健康状态实时异常检测方法
CN113029242A (zh) 结构健康监测系统中光纤光栅传感器异常诊断方法
Lu et al. Early fault warning and identification in condition monitoring of bearing via wavelet packet decomposition coupled with graph
He et al. The diagnosis of satellite flywheel bearing cage fault based on two-step clustering of multiple acoustic parameters
CN116627116B (zh) 一种流程工业故障定位方法、系统及电子设备
CN109840386B (zh) 基于因子分析的损伤识别方法
CN116641941A (zh) 一种基于典型变量分析的液压系统早期故障动态检测方法
CN117370879A (zh) 一种风力机齿轮箱的实时在线故障诊断方法及系统
CN114674511B (zh) 一种用于剔除时变环境因素影响的桥梁模态异常预警方法
CN114090949A (zh) 一种鲁棒的振动信号特征值计算方法
CN110108489B (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