WO2022156330A1 - 一种旋转设备故障诊断方法 - Google Patents

一种旋转设备故障诊断方法 Download PDF

Info

Publication number
WO2022156330A1
WO2022156330A1 PCT/CN2021/130848 CN2021130848W WO2022156330A1 WO 2022156330 A1 WO2022156330 A1 WO 2022156330A1 CN 2021130848 W CN2021130848 W CN 2021130848W WO 2022156330 A1 WO2022156330 A1 WO 2022156330A1
Authority
WO
WIPO (PCT)
Prior art keywords
fault
parameter
value
calculate
sample
Prior art date
Application number
PCT/CN2021/130848
Other languages
English (en)
French (fr)
Inventor
李倩
杨皓杰
孙丰诚
Original Assignee
杭州安脉盛智能技术有限公司
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 杭州安脉盛智能技术有限公司 filed Critical 杭州安脉盛智能技术有限公司
Publication of WO2022156330A1 publication Critical patent/WO2022156330A1/zh

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
    • G06F2218/08—Feature extraction
    • G06F2218/10—Feature extraction by analysing the shape of a waveform, e.g. extracting parameters relating to peaks
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01M—TESTING STATIC OR DYNAMIC BALANCE OF MACHINES OR STRUCTURES; TESTING OF STRUCTURES OR APPARATUS, NOT OTHERWISE PROVIDED FOR
    • G01M99/00—Subject matter not provided for in other groups of this subclass
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/14—Fourier, Walsh or analogous domain transformations, e.g. Laplace, Hilbert, Karhunen-Loeve, transforms
    • G06F17/148—Wavelet transforms
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/25—Fusion techniques
    • G06F18/253—Fusion techniques of extracted features
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis

Definitions

  • the invention relates to the technical field of fault diagnosis, in particular to a fault diagnosis method for rotating equipment.
  • the fault diagnosis of rotating equipment mainly relies on the frequency spectrum analysis and trend analysis of the vibration signals collected on site by technicians, and relies on the vibration signal processing theory and personal experience to realize fault judgment.
  • the requirements for technicians are too high, which not only requires rich on-site experience, but also professional theoretical knowledge; 2.
  • the signals collected on-site contain noise signals. It is very difficult to select parameters.
  • Classifier methods such as support vector machines, nearest neighbors, neural networks, decision trees, etc., implement data classification.
  • the extraction method of characteristic parameters is based on the extraction of time-domain characteristic parameters and frequency-domain characteristic parameters, and signal processing methods such as envelope spectrum analysis, wavelet analysis and empirical mode decomposition are used to process the vibration signal to obtain each component signal, and then Then extract the time domain and frequency domain characteristic parameters on different component signals.
  • This method does not extract the key fault signals in the equipment, resulting in a large amount of redundancy in the extracted feature parameters, and there are many invalid features.
  • different types of faults are processed with the same parameters, and the weights of characteristic parameters under different faults are not considered.
  • the main disadvantages of the prior art are: 1.
  • the extracted characteristic parameters are characteristic parameters of all frequency bands, which cannot adaptively extract characteristic parameters of the fault frequency band with the change of the fault mode, resulting in redundant characteristic parameters. 2.
  • the characteristic parameters are the same, and the weights between the characteristic parameters are not considered.
  • the invention solves the problems of difficulty in selecting fault frequency bands, consistent feature parameter weights, and easy influence on fault diagnosis results in the existing fault diagnosis methods for rotating equipment. , and can calculate the weight of each characteristic parameter under the model according to different fault models, realize the failure mode identification through the multi-parameter state estimation technology, improve the fault diagnosis accuracy, and realize the diagnosis efficiency. Model to realize failure mode identification of rotating equipment in industrial field.
  • a fault diagnosis method for rotating equipment comprising the steps of:
  • Step A Perform wavelet packet decomposition on the vibration signal, determine the number of layers k of wavelet packet decomposition according to the comprehensive evaluation index of kurtosis and energy entropy, and perform wavelet packet decomposition, extract the component signal with the largest comprehensive evaluation of kurtosis value and energy entropy, and analyze Extract the characteristic parameters of the selected component signal to obtain the characteristic parameter vector;
  • Step B Use reliefF algorithm to calculate the characteristic parameter weight, establish the sample parameter matrix under different fault types according to the characteristic parameter weight, and regenerate the fault model and the corresponding characteristic weight value;
  • Step C Calculate the estimated value of the characteristic parameters of the verification data under each model through the multi-parameter state estimation technology, and fuse the estimated values of all characteristic parameters and the residuals of the actual values into a corresponding distribution interval of a margin value set;
  • Step D Calculate the confidence probability of the data to be tested under different models, and obtain the fault category to which the data to be tested belongs according to the confidence rate.
  • the invention calculates the optimal decomposition layer number of wavelet packets through comprehensive evaluation of kurtosis value and information entropy, so as to realize adaptive extraction of vibration fault signals; through reliefF algorithm, characteristic parameters under fault types can be selected to obtain characteristic parameter weights and reduce characteristic parameters. Redundancy is established, and fault models under different failure modes are established; the margin value and confidence between the parameters to be measured and the fault model are calculated by the multi-parameter estimation algorithm, and the fault category can be judged by the confidence, which can effectively improve the efficiency of fault diagnosis.
  • the invention can adaptively determine the optimal decomposition level of the signal according to the characteristics of the vibration signal, decompose the signal to extract the characteristic signal of the fault and obtain the characteristic parameters, and establish different fault models according to the fault type. Not all the same. Finally, by calculating the confidence interval of the signal to be tested under different fault models, the fault pattern recognition result is obtained.
  • the step A specifically includes the following steps:
  • Step A1 Set the decomposition level k to increase from 2 to n, perform wavelet packet decomposition on the vibration signal in turn, obtain 2 k component signals, calculate the comprehensive evaluation index of the kurtosis and energy entropy of each component signal, and assume the length of the vibration signal is N, the component signal is expressed as x, and the formula of the evaluation index is:
  • pi is the energy ratio of each component signal
  • Step A2 record the maximum evaluation index under different decomposition levels, and determine the optimal decomposition level according to the maximum value in the evaluation index curve;
  • Step A3 After the number of decomposition layers is determined, wavelet packet decomposition is performed on the vibration signal, and feature parameter extraction is performed on the component signal with the largest evaluation index.
  • the feature parameter extraction includes effective value, peak value, kurtosis, kurtosis index, waveform index and frequency domain feature parameters.
  • the step B specifically includes the following steps:
  • Step B1 Construct a characteristic parameter data set Z, the characteristic parameter data set Z contains all the samples in the test data, the faults have a total of b groups of fault categories, b i represents the name of the ith fault category, and the reliefF algorithm is used to calculate the different faults of the equipment.
  • the feature parameter weight under the type, reliefF algorithm is a feature parameter selection method, which can effectively realize the feature parameter selection and feature weight confirmation.
  • Step B2 Set the weight value w of all feature parameters to 0, and set the feature parameter selection result to empty;
  • Step B4 Calculate the weight of each feature parameter, the calculation formula is:
  • P(b q ) represents the proportion of the number of samples of the b q fault category in the total number of samples
  • P(b i ) represents the proportion of the number of samples of the b i fault category in the total number of samples
  • diff(A, R 1 , R 2 ) represents the difference between samples
  • Step B5 Establish a sample parameter matrix under different fault types, the characteristic parameters under each fault type are screened according to the calculated weight value, D i is the sample parameter matrix under the i-th type of fault, and the rows represent the parameters in different samples value, the column represents all the feature parameters obtained under one sample, the number of feature parameters is n, and the number of samples is m, then the sample matrix D i is expressed as:
  • Step B6 Correct the corresponding parameter weights according to the characteristic parameters selected in the sample matrix.
  • the correction method is:
  • the step C specifically includes the following steps:
  • Step C1 Use the data of the validation set to calculate the estimated value of the characteristic parameters under each fault model. Assuming that the data of the validation set is x obs , the estimated value x est is calculated by the formula:
  • Step C2 Calculate the residual between the estimated value x est and the validation set data:
  • Step C3 Calculate the fused margin value y according to the sample weight:
  • Step C4 Calculate the margin value y of all the verification data under the corresponding fault category, calculate the mean value and variance, and obtain the Gaussian distribution of each fault model.
  • the step D specifically includes the following steps: using the methods from steps C1 to C3 to calculate the margin value of the sample to be tested under the fault model, and judging the probability probability of the margin value under several fault models, which has the largest The probability density is used as the actual fault diagnosis result.
  • the beneficial effects of the present invention are: through the evaluation indexes of energy kurtosis and energy entropy, the optimal decomposition layer number of wavelet packets can be determined adaptively, and fault characteristic signals and corresponding characteristic parameters can be extracted more accurately; multi-dimensional characteristic parameters are extracted, and The feature parameter weight calculation is realized by the reliefF algorithm, and the key feature parameters can be extracted according to the weight screening parameters to avoid parameter redundancy. Models with different weights are established according to the failure mode, and the failure mode identification is realized through the multi-parameter state estimation technology, which can improve the fault diagnosis accuracy and realize the diagnosis efficiency.
  • Figure 4 shows the peak parameter estimates and residuals of the inner-ring fault samples in the validation set under the inner-ring fault model
  • This embodiment proposes a fault diagnosis method for rotating equipment, with reference to FIG. 1 , including steps:
  • Step A Perform wavelet packet decomposition on the vibration signal, determine the number of layers k of wavelet packet decomposition according to the comprehensive evaluation index of kurtosis and energy entropy, and perform wavelet packet decomposition, extract the component signal with the largest comprehensive evaluation of kurtosis value and energy entropy, and analyze Extract the characteristic parameters of the selected component signal to obtain the characteristic parameter vector;
  • Step A specifically includes the following steps:
  • Step A1 Set the decomposition level k to increase from 2 to 6, perform wavelet packet decomposition on the vibration signal in turn, obtain 2 k component signals, calculate the comprehensive evaluation index of the kurtosis and energy entropy of each component signal, and assume the length of the vibration signal is N, the component signal is expressed as x, and the formula of the evaluation index is:
  • pi is the energy ratio of each component signal
  • Step A2 Record the maximum evaluation index under different decomposition levels, and determine the optimal decomposition level according to the maximum value in the evaluation index curve.
  • Figure 2 is the evaluation index curve of a certain group of data under the decomposition level of 2-6 layers.
  • the evaluation index value of the third group is the largest, so the number of decomposition levels is selected to be 3.
  • Step A3 Perform 3-layer wavelet packet decomposition on the vibration signal, and perform feature parameter extraction on the component signal with the largest evaluation index.
  • the characteristic parameters extracted here are: vibration RMS, kurtosis, kurtosis index, spectral center, square root amplitude, waveform index, skewness index, margin index and zero peak value, etc.
  • Step B Use reliefF algorithm to calculate the characteristic parameter weight, establish the sample parameter matrix under different fault types according to the characteristic parameter weight, and regenerate the fault model and the corresponding characteristic weight value;
  • Step B specifically includes the following steps:
  • Step B1 Construct the characteristic parameter data set Z, which contains all the samples in the test data, including 168 sets of inner ring faults, 168 sets of rolling element faults, 339 sets of outer ring faults and 254 sets of normal sample data, a total of 929 sets of samples , 4 types of faults, the reliefF algorithm is used to calculate the characteristic parameter weights under different types of equipment faults;
  • Step B2 Set the weight value w of all feature parameters to 0, and set the feature parameter selection result to empty;
  • Step B4 Calculate the weight of each feature parameter, the calculation formula is:
  • P(b q ) represents the proportion of the number of samples of the b q fault category in the total number of samples
  • P(b i ) represents the proportion of the number of samples of the b i fault category in the total number of samples
  • diff(A, R 1 , R 2 ) represents the difference between samples
  • Figure 3 shows the result of the average weight of feature parameters calculated for 266 times of data.
  • Step B5 Establish a sample parameter matrix under different fault types, the characteristic parameters under each fault type are screened according to the calculated weight value, D i is the sample parameter matrix under the i-th type of fault, and the rows represent the parameters in different samples value, the column represents all the feature parameters under a sample, the number of feature parameters is n, and the number of samples is m, then the sample matrix D i is expressed as:
  • Step B6 Correct the corresponding parameter weights according to the characteristic parameters selected in the sample matrix.
  • the correction method is:
  • Table 1 shows the weight results of the modified feature parameters.
  • Table 1 is the result of feature parameter weight
  • Characteristic Parameters Weights Characteristic Parameters Weights valid value 0.242 Waveform Indicator 0.163 kurtosis value 0.199 kurtosis index 0.093 peak indicator 0.047 peak 0.254
  • Step C Calculate the estimated values of the characteristic parameters of the verification data under each model through the multi-parameter state estimation technique, and fuse the estimated values of all characteristic parameters and the residuals of the actual values into a corresponding distribution interval of a margin value set; Step C is specific Include the following steps:
  • Step C1 Use the data of the validation set to calculate the estimated value of the characteristic parameters under each fault model. Assuming that the data of the validation set is x obs , the estimated value x est is calculated by the formula:
  • Step C2 Calculate the residual between the estimated value x est and the validation set data:
  • Figure 4 shows the real-time estimates and residuals of the rolling element fault data in the validation set under the rolling element fault model and the outer ring fault model.
  • Step C3 Calculate the fused margin value y according to the sample weight:
  • FIG. 5 is a result of fusion of margin values of peak parameters.
  • Step C4 Calculate the margin value y of all the verification data under the corresponding fault category, calculate the mean and variance, and obtain the Gaussian distribution of each fault model, as shown in Table 2.
  • Fault model name mean standard deviation Inner ring failure model 0.00175 0.00178 Outer ring failure model 0.005 0.0184 Rolling Element Failure Model 0.0025 0.0025 normal model 0.0022 0.0028
  • Step D Calculate the confidence probability of the data to be tested under different models, and obtain the fault category to which the data to be tested belongs according to the confidence rate.
  • Step D specifically includes the following steps: using the methods from steps C1 to C3 to calculate the margin value of the sample to be tested under the fault model, and judging the probability probability of the margin value under several fault models, the one with the largest probability density is used as the actual value. fault diagnosis results. Table 3 shows the final diagnosis results, and the total diagnosis accuracy rate is 98.53%.
  • the method proposed in the present invention can effectively identify the faults of rotating machinery equipment.
  • the invention calculates the optimal decomposition layer number of wavelet packets through comprehensive evaluation of kurtosis value and information entropy, so as to realize adaptive extraction of vibration fault signals; through reliefF algorithm, characteristic parameters under fault types can be selected to obtain characteristic parameter weights and reduce characteristic parameters. Redundancy is established, and fault models under different failure modes are established; the margin value and confidence between the parameters to be measured and the fault model are calculated by the multi-parameter estimation algorithm, and the fault category can be judged by the confidence, which can effectively improve the efficiency of fault diagnosis.
  • the invention can adaptively determine the optimal decomposition level of the signal according to the characteristics of the vibration signal, decompose the signal to extract the characteristic signal of the fault and obtain the characteristic parameters, and establish different fault models according to the fault type. Not all the same. Finally, by calculating the confidence interval of the signal to be tested under different fault models, the fault pattern recognition result is obtained.

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Data Mining & Analysis (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Mathematical Physics (AREA)
  • General Engineering & Computer Science (AREA)
  • Mathematical Optimization (AREA)
  • Computational Mathematics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Artificial Intelligence (AREA)
  • Mathematical Analysis (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Evolutionary Biology (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Evolutionary Computation (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • Signal Processing (AREA)
  • Testing Of Devices, Machine Parts, Or Other Structures Thereof (AREA)

Abstract

一种旋转设备故障诊断方法，包括：步骤A：对振动信号进行小波包分解，提取峭度值和能量熵综合评估最大的分量信号，并对选择出来的分量信号进行特征参数提取，得到特征参数向量；步骤B：使用reliefF算法计算特征参数权重，依照特征参数权重建立不同故障类型下的样本参数矩阵，并重新生成故障模型和对应特征权重值；步骤C：通过多参数状态估计技术计算验证数据在各模型下的特征参数估计值，根据所有特征参数的估计值与实际值的残差融合成一个裕度值集的对应分布区间；步骤D：计算待测数据在不同模型下的置信概率，并根据置信率得到待测数据所属故障类别。本方法可提高故障诊断精度，实现诊断效率。

Description

一种旋转设备故障诊断方法 技术领域
本发明涉及故障诊断技术领域，尤其是一种旋转设备故障诊断方法。
背景技术
在各种工业现场有大量的旋转设备，旋转设备结构复杂，受环境和噪声影响容易发生故障。任何这些设备出现故障都有可能会造成设备非计划停机，工作人员需要花大量时间来进行故障定位、维修、订购更换零件等，从而产生巨大的经济损失。因此及时发现设备运行异常，准确诊断出旋转设备故障类型，为维修人员检查和维修提供指导建议具有重大意义。
目前，旋转设备的故障诊断主要是通过技术人员对现场采集的振动信号进行频谱分析、趋势分析等，依靠振动信号处理理论和个人经验来实现故障判断。但这种方式存在两个问题：1、对技术人员的要求过高，不仅需要丰富的现场经验，也需要专业的理论知识；2、现场采集的信号中包含噪声信号，在对信号频带设置和参数选择方面有很大难度。近年来，随着机器学习研究的深入，数据模型也逐渐在旋转设备故障诊断中有一定的应用，常见方法是对采集的振动信号进行特征提取，对不同故障类型数据打标签，然后通过各种分类器方法，比如支持向量机、最近邻、神经网络、决策树等方式实现数据分类。特征参数的提取方法，由提取时域特征参数，频域特征参数基础上，引用信号处理方法如包络谱分析、小波分析和经验模态分解等，对振动信号进行处理得到各分量信号，然后再提取不同分量信号上的时域、频域特征参数。此种方法没有提取到设备中的关键故障信号，导致提取的特征参数冗余量较大，存在很多无效的特征。并且，不同类型故障均使用同样的参数做处理，没有考虑不同故障下的特征参数权重。
现有的技术最主要的缺点是：1.所提取的特征参数是所有频带的特征参数，其不能自适应的随着故障模式的变化而提取故障频带的特征参数，造成特征参数冗余。2.所有类型的故障模型中，特征参数相同，且没有考虑特征参数之间的权重。这两点缺点会导致有效特征参数和无效特征参数的权重一样，故障识别结果不够准确，诊断错误的概率较高。
发明内容
本发明解决了现有旋转设备故障诊断方法中，故障频带选择困难，特征参数权重一致，容易影响故障诊断结果的问题，提出一种考虑特征参数权重的，能自适应提取特征信号中的故障特征，并能够针对不同故障模型计算模型下各特征参数的权重，通过多参数状态估计技术实现故障模式识别，可提高故障诊断精度，实现诊断效率，通过relief方法和多参数状态估计算法建立旋转设备故障模型，实现工业现场旋转设备故障模式识别。
为实现上述目的，提出以下技术方案：
一种旋转设备故障诊断方法，包括步骤：
步骤A：对振动信号进行小波包分解，根据峭度和能量熵综合评价指标确定小波包分解层数k，并进行小波包分解，提取峭度值和能量熵综合评估最大的分量信号，并对选择出来的分量信号进行特征参数提取，得到特征参数向量；
步骤B：使用reliefF算法计算特征参数权重，依照特征参数权重建立不同故障类型下的样本参数矩阵，并重新生成故障模型和对应特征权重值；
步骤C：通过多参数状态估计技术计算验证数据在各模型下的特征参数估计值，根据所有特征参数的估计值与实际值的残差融合成一个裕度值集的对应分布区间；
步骤D：计算待测数据在不同模型下的置信概率，并根据置信率得到待测数据所属故障类别。
本发明通过峭度值和信息熵综合综合评价计算小波包最优分解层数，实现振动故障信号自适应抽取；通过reliefF算法可对故障类型下的特征参数进行选择，得到特征参数权重，降低特征冗余度，并建立不同故障模式下的故障模型；通过多参数估计算法计算待测参数与故障模型间的裕度值和置信度，通过置信度判断故障所属类别，可有效提高故障诊断效率。本发明能根据振动信号特征自适应确定信号的最优分解层数，对信号分解提取故障的特征信号并获取特征参数，并根据故障类型建立不同故障模型，各模型基于不同的特征参数，权重值也不尽相同。最后通过计算待测信号在不同故障模型下的置信区间，得到故障模式识别结果。
作为优选，所述步骤A具体包括以下步骤：
步骤A1：设置分解层数k由2递增至n，依次对振动信号进行小波包分解，得到2 k个分量信号，计算各分量信号的峭度和能量熵的综合评价指标，假设振动信号的长度为N，分量信号表示为x，评价指标的公式为：
Figure PCTCN2021130848-appb-000001
其中，p i为每个分量信号的能量占比；
步骤A2：记录在不同分解层数下最大的评价指标，根据评价指标曲线中的最大值确定最优分解层数；
步骤A3：分解层数确定后，对振动信号进行小波包分解，并对评价指标最大的分量信号进行特征参数提取。
作为优选，所述特征参数提取包括有效值、峰值、峭度、峭度指标、波形指标和频域特征参数。
作为优选，所述步骤B具体包括以下步骤：
步骤B1：构建特征参数数据集Z，所述特征参数数据集Z中包含测试数据中的所有样本，故障共有b组故障类别，b i代表第i个故障类别名称，使用reliefF算法计算设备故障不同类型下的特征参数权重，reliefF算法是一种特征参数选择方法，能够有效地实现特征参数选择及特征权重确认。
步骤B2：将所有特征参数的权重值w设置为0，特征参数选择结果设置为空；
步骤B3：从全量矩阵中有放回地选择出第j个样本，第j个样本对应的故障类别为b i，先在属于b i故障类别的矩阵中筛选出前s个最相邻的样本，构成H l(l＝1,2,...,s)，按照同样的方法筛选在其他故障类别的中的前s个最相邻的样本，构成M l(b q)(l＝1,2,...,s)，M l(b q)代表属于b q(q≠i)故障类别的下的第l个最相似的样本；
步骤B4：计算每个特征参数的权重，计算公式为：
Figure PCTCN2021130848-appb-000002
其中，P(b q)代表b q故障类别的样本个数在总样本个数中的占比，P(b i)代表b i故障类别的样本个数在总样本个数中的占比，diff(A,R 1,R 2)代表样本间的差异；
Figure PCTCN2021130848-appb-000003
步骤B5：建立不同故障类型下的样本参数矩阵，每个故障类型下的特征参数根据计算得到的权重值进行筛选，D i是第i类故障下的样本参数矩阵，行代表不同样本中的参数值，列代表一个样本下得到所有特征参数，特征参数个数为n，样本个数为m，则样本矩阵D i表示为：
Figure PCTCN2021130848-appb-000004
步骤B6：根据样本矩阵中选择的特征参数，修正对应参数权重，修正方法为：
Figure PCTCN2021130848-appb-000005
其中，
Figure PCTCN2021130848-appb-000006
代表选中的P个特征参数的权重和。
作为优选，所述步骤C具体包括以下步骤：
步骤C1：使用验证集的数据计算在各个故障模型下的特征参数估计值，假设验证集的数据为x obs,估计值x est通过公式计算得到：
Figure PCTCN2021130848-appb-000007
式中，
Figure PCTCN2021130848-appb-000008
指的是欧氏距离、皮尔逊相关系数、曼哈顿距离；
步骤C2：计算估计值x est与验证集数据的残差：
res＝x obs-x est
步骤C3：根据样本权重计算得到融合的裕度值y：
Figure PCTCN2021130848-appb-000009
步骤C4：计算所有验证数据在对应的故障类别下的裕度值y，计算均值和方差，得到各故障模型的高斯分布。
作为优选，所述步骤D具体包括以下步骤：采用步骤C1至步骤C3的方法计算待测样本在故障模型下的裕度值，并判断裕度值在几个故障模型下的概率概率，拥有最大概率密度的作为实际的故障诊断结果。
本发明的有益效果是：通过能峭度和能量熵的评价指标，能够自适应确定小波包最优分解层数，能够更准确地提取到故障特征信号及对应特征参数；提取多维特征参数，并通过reliefF算法实现特征参数权重计算，根据权重筛选参数能够提取出关键特征参数，避免参数冗余。根据故障模式建立具有不同权重的模型，通过多参数状态估计技术实现故障模式识别，可提高故障诊断精度，实现诊断效率。
附图说明
图1实施例的算法计算流程图；
图2实施例的小波包分解综合指标曲线图；
图3实施例的reliefF算法计算的特征权重结果图；
图4实施例的验证集内圈故障样本在内圈故障模型下的峰值参数估计值及残差；
图5实施例的峰值参数的裕度值融合结果。
具体实施方式
实施例：
结合附图和具体实施例对本发明做进一步的描述，一种旋转设备的故障诊断方法，以凯斯西储轴承模拟实验台上，转动速度为1797r/min状态下采集的驱动端加速度数据为例，采 样频率为12kHz。
本实施例提出一种旋转设备故障诊断方法，参考图1，包括步骤：
步骤A：对振动信号进行小波包分解，根据峭度和能量熵综合评价指标确定小波包分解层数k，并进行小波包分解，提取峭度值和能量熵综合评估最大的分量信号，并对选择出来的分量信号进行特征参数提取，得到特征参数向量；
步骤A具体包括以下步骤：
步骤A1：设置分解层数k由2递增至6，依次对振动信号进行小波包分解，得到2 k个分量信号，计算各分量信号的峭度和能量熵的综合评价指标，假设振动信号的长度为N，分量信号表示为x，评价指标的公式为：
Figure PCTCN2021130848-appb-000010
其中，p i为每个分量信号的能量占比；
步骤A2：记录在不同分解层数下最大的评价指标，根据评价指标曲线中的最大值确定最优分解层数。图2是某组数据在2-6层分解层数下的评价指标曲线，其中第3组的评价指标值最大，因此选择分解层数为3。
步骤A3：对振动信号进行3层小波包分解，并对评价指标最大的分量信号进行特征参数提取。此处提取的特征参数为：振动有效值、峭度、峭度指标、谱中心、方根幅值、波形指标、偏斜度指标、裕度指标和零峰值等。
步骤B：使用reliefF算法计算特征参数权重，依照特征参数权重建立不同故障类型下的样本参数矩阵，并重新生成故障模型和对应特征权重值；
步骤B具体包括以下步骤：
步骤B1：构建特征参数数据集Z，Z中包含测试数据中的所有样本，其中包含168组内圈故障、168组滚动体故障、339组外圈故障和254组正常样本数据，共计929组样本，4类故障，使用reliefF算法计算设备故障不同类型下的特征参数权重；
步骤B2：将所有特征参数的权重值w设置为0，特征参数选择结果设置为空；
步骤B3：从全量矩阵中有放回地选择出第j个样本，第j个样本对应的故障类别为b i，先在属于b i故障类别的矩阵中筛选出前6个最相邻的样本，构成H l(l＝1,2,...,s)，按照同样的方法筛选在其他故障类别的中的前6个最相邻的样本，构成M l(b q)(l＝1,2,...,s)，M l(b q)代表属于b q(q≠i)故障类别的下的第l个最相似的样本；其中设置s为10。
步骤B4：计算每个特征参数的权重，计算公式为：
Figure PCTCN2021130848-appb-000011
其中，P(b q)代表b q故障类别的样本个数在总样本个数中的占比，P(b i)代表b i故障类别的样本个数在总样本个数中的占比，diff(A,R 1,R 2)代表样本间的差异；
Figure PCTCN2021130848-appb-000012
图3为266次数据计算的特征参数平均权重结果。
步骤B5：建立不同故障类型下的样本参数矩阵，每个故障类型下的特征参数根据计算得到的权重值进行筛选，D i是第i类故障下的样本参数矩阵，行代表不同样本中的参数值，列代表一个样本下的所有特征参数，特征参数个数为n，样本个数为m，则样本矩阵D i表示为：
Figure PCTCN2021130848-appb-000013
步骤B6：根据样本矩阵中选择的特征参数，修正对应参数权重，修正方法为：
Figure PCTCN2021130848-appb-000014
其中，
Figure PCTCN2021130848-appb-000015
代表选中的P个特征参数的权重和。
选择权重值超过0.03的特征参数，并对权重值进行修正，表1为修正后的特征参数权重结果。
表1是特征参数权重结果
特征参数 权重值 特征参数 权重值
有效值 0.242 波形指标 0.163
峭度值 0.199 峭度指标 0.093
峰值指标 0.047 峰值 0.254
步骤C：通过多参数状态估计技术计算验证数据在各模型下的特征参数估计值，根据所有特征参数的估计值与实际值的残差融合成一个裕度值集的对应分布区间；步骤C具体包括以下步骤：
步骤C1：使用验证集的数据计算在各个故障模型下的特征参数估计值，假设验证集的数据为x obs,估计值x est通过公式计算得到：
Figure PCTCN2021130848-appb-000016
式中，
Figure PCTCN2021130848-appb-000017
指的是欧氏距离、皮尔逊相关系数、曼哈顿距离；此处使用欧式距离方法实现参数估计，欧式距离的计算公式是：
Figure PCTCN2021130848-appb-000018
步骤C2：计算估计值x est与验证集数据的残差：
res＝x obs-x est
图4是验证集中滚动体故障数据在滚动体故障模型和外圈故障模型下的实时估计值和残差。
步骤C3：根据样本权重计算得到融合的裕度值y：
Figure PCTCN2021130848-appb-000019
参考图5，图5是峰值参数的裕度值融合结果。
步骤C4：计算所有验证数据在对应的故障类别下的裕度值y，计算均值和方差，得到各故障模型的高斯分布,如表2所示。
表2故障模型裕度范围
故障模型名称 均值 标准差
内圈故障模型 0.00175 0.00178
外圈故障模型 0.005 0.0184
滚动体故障模型 0.0025 0.0025
正常模型 0.0022 0.0028
步骤D：计算待测数据在不同模型下的置信概率，并根据置信率得到待测数据所属故障类别。
步骤D具体包括以下步骤：采用步骤C1至步骤C3的方法计算待测样本在故障模型下的裕度值，并判断裕度值在几个故障模型下的概率概率，拥有最大概率密度的作为实际的故障诊断结果。表3为最终诊断结果，总诊断准确率为98.53％，本发明提出的方法可有效识别旋转机械设备故障。
表3最终诊断结果
Figure PCTCN2021130848-appb-000020
Figure PCTCN2021130848-appb-000021
本发明通过峭度值和信息熵综合综合评价计算小波包最优分解层数，实现振动故障信号自适应抽取；通过reliefF算法可对故障类型下的特征参数进行选择，得到特征参数权重，降低特征冗余度，并建立不同故障模式下的故障模型；通过多参数估计算法计算待测参数与故障模型间的裕度值和置信度，通过置信度判断故障所属类别，可有效提高故障诊断效率。本发明能根据振动信号特征自适应确定信号的最优分解层数，对信号分解提取故障的特征信号并获取特征参数，并根据故障类型建立不同故障模型，各模型基于不同的特征参数，权重值也不尽相同。最后通过计算待测信号在不同故障模型下的置信区间，得到故障模式识别结果。

Claims (6)

  1. 一种旋转设备故障诊断方法，其特征是，包括步骤：
    步骤A：对振动信号进行小波包分解，根据峭度和能量熵综合评价指标确定小波包分解层数k，并进行小波包分解，提取峭度值和能量熵综合评估最大的分量信号，并对选择出来的分量信号进行特征参数提取，得到特征参数向量；
    步骤B：使用reliefF算法计算特征参数权重，依照特征参数权重建立不同故障类型下的样本参数矩阵，并重新生成故障模型和对应特征权重值；
    步骤C：通过多参数状态估计技术计算验证数据在各模型下的特征参数估计值，根据所有特征参数的估计值与实际值的残差融合成一个裕度值集的对应分布区间；
    步骤D：计算待测数据在不同模型下的置信概率，并根据置信率得到待测数据所属故障类别。
  2. 根据权利要求1所述的一种旋转设备故障诊断方法，其特征是，所述步骤A具体包括以下步骤：
    步骤A1：设置分解层数k由2递增至n，依次对振动信号进行小波包分解，得到2 k个分量信号，计算各分量信号的峭度和能量熵的综合评价指标，假设振动信号的长度为N，分量信号表示为x，评价指标的公式为：
    Figure PCTCN2021130848-appb-100001
    其中，p i为每个分量信号的能量占比；
    步骤A2：记录在不同分解层数下最大的评价指标，根据评价指标曲线中的最大值确定最优分解层数；
    步骤A3：分解层数确定后，对振动信号进行小波包分解，并对评价指标最大的分量信号进行特征参数提取。
  3. 根据权利要求2所述的一种旋转设备故障诊断方法，其特征是，所述特征参数提取包括有效值、峰值、峭度、峭度指标、波形指标和频域特征参数。
  4. 根据权利要求1所述的一种旋转设备故障诊断方法，其特征是，所述步骤B具体包括以下步骤：
    步骤B1：构建特征参数数据集Z，所述特征参数数据集Z中包含测试数据中的所有样本，故障共有b组故障类别，b i代表第i个故障类别名称，使用reliefF算法计算设备故障不同 类型下的特征参数权重；
    步骤B2：将所有特征参数的权重值w设置为0，特征参数选择结果设置为空；
    步骤B3：从全量矩阵中有放回地选择出第j个样本，第j个样本对应的故障类别为b i，先在属于b i故障类别的矩阵中筛选出前s个最相邻的样本，构成H l(l＝1,2,...,s)，按照同样的方法筛选在其他故障类别的中的前s个最相邻的样本，构成M l(b q)(l＝1,2,...,s)，M l(b q)代表属于b q(q≠i)故障类别的下的第l个最相似的样本；
    步骤B4：计算每个特征参数的权重，计算公式为：
    Figure PCTCN2021130848-appb-100002
    其中，P(b q)代表b q故障类别的样本个数在总样本个数中的占比，P(b i)代表b i故障类别的样本个数在总样本个数中的占比，diff(A,R 1,R 2)代表样本间的差异；
    Figure PCTCN2021130848-appb-100003
    步骤B5：建立不同故障类型下的样本参数矩阵，每个故障类型下的特征参数根据计算得到的权重值进行筛选，D i是第i类故障下的样本参数矩阵，行代表不同样本中的参数值，列代表一个样本下得到所有特征参数，特征参数个数为n，样本个数为m，则样本矩阵D i表示为：
    Figure PCTCN2021130848-appb-100004
    步骤B6：根据样本矩阵中选择的特征参数，修正对应参数权重，修正方法为：
    Figure PCTCN2021130848-appb-100005
    其中，
    Figure PCTCN2021130848-appb-100006
    代表选中的P个特征参数的权重和。
  5. 根据权利要求1所述的一种旋转设备故障诊断方法，其特征是，所述步骤C具体包括以下步骤：
    步骤C1：使用验证集的数据计算在各个故障模型下的特征参数估计值，假设验证集的数据为x obs,估计值x est通过公式计算得到：
    Figure PCTCN2021130848-appb-100007
    式中，
    Figure PCTCN2021130848-appb-100008
    指的是欧氏距离、皮尔逊相关系数、曼哈顿距离；
    步骤C2：计算估计值x est与验证集数据的残差：
    res＝x obs-x est
    步骤C3：根据样本权重计算得到融合的裕度值y：
    Figure PCTCN2021130848-appb-100009
    步骤C4：计算所有验证数据在对应的故障类别下的裕度值y，计算均值和方差，得到各故障模型的高斯分布。
  6. 根据权利要求5所述的一种旋转设备故障诊断方法，其特征是，所述步骤D具体包括以下步骤：采用步骤C1至步骤C3的方法计算待测样本在故障模型下的裕度值，并判断裕度值在几个故障模型下的概率概率，拥有最大概率密度的作为实际的故障诊断结果。
PCT/CN2021/130848 2021-01-19 2021-11-16 一种旋转设备故障诊断方法 WO2022156330A1 (zh)

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
CN202110068527.8 2021-01-19
CN202110068527.8A CN112906473B (zh) 2021-01-19 2021-01-19 一种旋转设备故障诊断方法

Publications (1)

Publication Number Publication Date
WO2022156330A1 true WO2022156330A1 (zh) 2022-07-28

Family

ID=76115517

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/CN2021/130848 WO2022156330A1 (zh) 2021-01-19 2021-11-16 一种旋转设备故障诊断方法

Country Status (2)

Country Link
CN (1) CN112906473B (zh)
WO (1) WO2022156330A1 (zh)

Cited By (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115096634A (zh) * 2022-08-26 2022-09-23 启东市嘉信精密机械有限公司 一种机械设备运行故障检测方法及系统
CN115356598A (zh) * 2022-10-20 2022-11-18 国网四川省电力公司电力科学研究院 一种配电网短路故障区段定位方法及系统
CN115683687A (zh) * 2023-01-03 2023-02-03 成都大汇物联科技有限公司 一种水电机械设备动静碰磨故障诊断方法
CN116342108A (zh) * 2023-04-06 2023-06-27 安徽配隆天环保科技有限公司 一种环保监测设备运行状态分析系统及分析方法
CN116610941A (zh) * 2023-07-21 2023-08-18 山东科技大学 快速峭度图轴承复合故障诊断方法、系统、设备以及介质
CN116975541A (zh) * 2023-09-21 2023-10-31 深圳市盘古环保科技有限公司 一种垃圾填埋场存量垃圾自动筛分系统
CN117077505A (zh) * 2023-06-14 2023-11-17 中国人民解放军海军航空大学 一种基于测试性评价的装备机内故障诊断器设计方法
CN117189720A (zh) * 2023-09-14 2023-12-08 成都飞航智云科技有限公司 一种飞行器液压系统的故障诊断方法
CN117367570A (zh) * 2023-11-02 2024-01-09 中国人民解放军海军工程大学 一种基于单点声信号的空气压缩机智能故障诊断方法

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112906473B (zh) * 2021-01-19 2023-06-20 杭州安脉盛智能技术有限公司 一种旋转设备故障诊断方法
CN113487141B (zh) * 2021-06-11 2023-09-29 北京控制工程研究所 一种多源信息聚类融合的轴承状态评估方法
CN114414242A (zh) * 2021-12-17 2022-04-29 哈尔滨工业大学 一种基于多域分析及ReliefF的大型回转设备主轴状态信号特征提取方法

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102841131A (zh) * 2012-09-20 2012-12-26 西安科技大学 钢丝绳芯输送带缺陷智能识别方法及系统
CN106769052A (zh) * 2017-03-21 2017-05-31 桂林电子科技大学 一种基于聚类分析的机械系统滚动轴承智能故障诊断方法
CN107609574A (zh) * 2017-08-18 2018-01-19 上海电力学院 基于数据挖掘的风电机组故障预警方法
CN108844725A (zh) * 2018-04-24 2018-11-20 天津大学 一种汽车发动机轴瓦磨损故障诊断方法
CN112906473A (zh) * 2021-01-19 2021-06-04 杭州安脉盛智能技术有限公司 一种旋转设备故障诊断方法

Family Cites Families (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2009115075A (ja) * 2007-10-15 2009-05-28 Toyota Motor Corp エンジンの油圧制御装置
CN105701470B (zh) * 2016-01-13 2019-02-15 合肥工业大学 一种基于最优小波包分解的模拟电路故障特征提取方法
CN107870892A (zh) * 2017-07-01 2018-04-03 南京理工大学 一种城轨列车滚动轴承故障诊断方法
CN107701378B (zh) * 2017-09-29 2019-09-27 上海电力设计院有限公司 一种风力发电机故障预警方法
CN108227676A (zh) * 2017-12-28 2018-06-29 浙江工业大学 阀控缸电液伺服系统在线故障检测、估计及定位方法
CN108318249B (zh) * 2018-01-24 2020-04-17 广东石油化工学院 一种旋转机械轴承的故障诊断方法
CN109993105A (zh) * 2019-03-29 2019-07-09 北京化工大学 一种改进的自适应稀疏采样故障分类方法
CN110045716B (zh) * 2019-04-18 2020-08-11 中广核工程有限公司 一种闭环控制系统早期故障检测和诊断方法和系统
CN110907826B (zh) * 2019-11-14 2022-04-08 中车株洲电力机车研究所有限公司 一种基于卷积神经网络滤波的电机故障诊断方法及系统
CN111551872B (zh) * 2020-02-27 2021-10-22 西北工业大学 一种pmsm驱动系统逆变器开路故障在线诊断方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102841131A (zh) * 2012-09-20 2012-12-26 西安科技大学 钢丝绳芯输送带缺陷智能识别方法及系统
CN106769052A (zh) * 2017-03-21 2017-05-31 桂林电子科技大学 一种基于聚类分析的机械系统滚动轴承智能故障诊断方法
CN107609574A (zh) * 2017-08-18 2018-01-19 上海电力学院 基于数据挖掘的风电机组故障预警方法
CN108844725A (zh) * 2018-04-24 2018-11-20 天津大学 一种汽车发动机轴瓦磨损故障诊断方法
CN112906473A (zh) * 2021-01-19 2021-06-04 杭州安脉盛智能技术有限公司 一种旋转设备故障诊断方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
WANG YOURONG: "Applied Research on Relieff Weighted Feature Selection Algorithm in Fault Diagnosis of Rotating Machinery", CHINESE MASTER'S THESES FULL-TEXT DATABASE, 1 December 2014 (2014-12-01), pages 1 - 73, XP055952214 *

Cited By (15)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115096634A (zh) * 2022-08-26 2022-09-23 启东市嘉信精密机械有限公司 一种机械设备运行故障检测方法及系统
CN115356598A (zh) * 2022-10-20 2022-11-18 国网四川省电力公司电力科学研究院 一种配电网短路故障区段定位方法及系统
CN115356598B (zh) * 2022-10-20 2023-02-03 国网四川省电力公司电力科学研究院 一种配电网短路故障区段定位方法及系统
CN115683687A (zh) * 2023-01-03 2023-02-03 成都大汇物联科技有限公司 一种水电机械设备动静碰磨故障诊断方法
CN116342108B (zh) * 2023-04-06 2023-09-05 安徽配隆天环保科技有限公司 一种环保监测设备运行状态分析系统及分析方法
CN116342108A (zh) * 2023-04-06 2023-06-27 安徽配隆天环保科技有限公司 一种环保监测设备运行状态分析系统及分析方法
CN117077505A (zh) * 2023-06-14 2023-11-17 中国人民解放军海军航空大学 一种基于测试性评价的装备机内故障诊断器设计方法
CN117077505B (zh) * 2023-06-14 2024-05-28 中国人民解放军海军航空大学 一种基于测试性评价的装备机内故障诊断器设计方法
CN116610941A (zh) * 2023-07-21 2023-08-18 山东科技大学 快速峭度图轴承复合故障诊断方法、系统、设备以及介质
CN116610941B (zh) * 2023-07-21 2023-09-22 山东科技大学 快速峭度图轴承复合故障诊断方法、系统、设备以及介质
CN117189720A (zh) * 2023-09-14 2023-12-08 成都飞航智云科技有限公司 一种飞行器液压系统的故障诊断方法
CN116975541A (zh) * 2023-09-21 2023-10-31 深圳市盘古环保科技有限公司 一种垃圾填埋场存量垃圾自动筛分系统
CN116975541B (zh) * 2023-09-21 2024-01-09 深圳市盘古环保科技有限公司 一种垃圾填埋场存量垃圾自动筛分系统
CN117367570A (zh) * 2023-11-02 2024-01-09 中国人民解放军海军工程大学 一种基于单点声信号的空气压缩机智能故障诊断方法
CN117367570B (zh) * 2023-11-02 2024-04-12 中国人民解放军海军工程大学 一种基于单点声信号的空气压缩机智能故障诊断方法

Also Published As

Publication number Publication date
CN112906473B (zh) 2023-06-20
CN112906473A (zh) 2021-06-04

Similar Documents

Publication Publication Date Title
WO2022156330A1 (zh) 一种旋转设备故障诊断方法
CN110940539B (zh) 一种基于人工经验及声音识别的机器设备故障诊断方法
CN108375476B (zh) 一种水电机组健康评估方法
CN104729853B (zh) 一种滚动轴承性能退化评估装置及方法
CN110867196A (zh) 一种基于深度学习及声音识别的机器设备状态监测系统
CN112257530B (zh) 基于盲信号分离和支持向量机的滚动轴承故障诊断方法
CN109781411A (zh) 一种结合改进稀疏滤波器与kelm的轴承故障诊断方法
CN107145675A (zh) 基于bp神经网络算法的电力变压器故障诊断装置及方法
CN113188794B (zh) 一种基于改进pso-bp神经网络齿轮箱故障诊断方法及装置
CN112305388B (zh) 一种发电机定子绕组绝缘局部放电故障在线监测诊断方法
CN112599134A (zh) 一种基于声纹识别的变压器声音事件检测方法
CN108470570B (zh) 电机异音检测方法
CN105241665A (zh) 一种基于IRBFNN-AdaBoost分类器的滚动轴承故障诊断方法
CN114757269A (zh) 一种基于局部子空间-邻域保持嵌入的复杂过程精细化故障检测方法
CN115496108A (zh) 一种基于流形学习和大数据分析的故障监测方法及系统
CN109299201B (zh) 基于两阶段聚类的电厂生产子系统异常监测方法及装置
CN107463872A (zh) 一种旋转机械转轴裂纹故障诊断方法
CN115293188A (zh) 一种用于往复机械设备的故障诊断方法和装置
Li et al. Transformer-based meta learning method for bearing fault identification under multiple small sample conditions
CN111783941B (zh) 一种基于概率置信度卷积神经网络的机械设备诊断分类方法
CN113469252A (zh) 一种考虑不平衡样本的特高压换流阀运行状态评估方法
CN106250937B (zh) 一种基于非相似度指标的故障分类诊断方法
CN110046419B (zh) 一种基于ctma-dl算法的动设备故障类型在线诊断方法
CN109784777B (zh) 基于时序信息片段云相似度度量的电网设备状态评估方法
CN116578833A (zh) 基于优化随机森林模型的igbt模块老化故障诊断系统

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 21920706

Country of ref document: EP

Kind code of ref document: A1

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 21920706

Country of ref document: EP

Kind code of ref document: A1