An unsupervised fault diagnosis method and system based on multi-scale feature enhancement and adaptive feature selection

Abstract

The application provides an unsupervised fault diagnosis method and system based on multi-scale feature enhancement and adaptive feature selection, and belongs to the field of industrial equipment fault diagnosis. In order to solve the problems that the existing feature extraction method has insufficient feature representation ability, single feature scale and serious redundant feature interference, the application extracts multi-dimensional statistical features such as quantile, quartile range and skewness, combines adaptive wavelet base selection of the information maximum redundancy minimum criterion, and significantly improves the representation ability of the multi-scale characteristics of the fault signal. The cosine similarity weighted Laplace score and variance double threshold mechanism are introduced, and the redundant features are effectively eliminated. Through the collaborative design of multi-scale interval enhancement and adaptive wavelet transform, the time-frequency domain features are effectively extracted, the time consumption of high-dimensional data processing is obviously reduced compared with the traditional elastic distance method, and the real-time demand of industrial online monitoring is met.

Description

Technical Field

[0001] This invention relates to the field of industrial equipment fault diagnosis technology, and more specifically, to an unsupervised fault diagnosis method and system based on multi-scale feature enhancement and adaptive feature selection. Background Technology

[0002] In recent years, the intelligentization of industrial equipment has accelerated, with sensor networks widely deployed in key components such as hydraulic systems and power units, generating multi-dimensional time-series data in real time, providing a data foundation for equipment condition monitoring and fault diagnosis. Existing unsupervised fault diagnosis methods face the following technical bottlenecks in industrial scenarios:

[0003] Traditional threshold detection methods identify anomalies by setting fixed thresholds for sensor signals (such as pressure or temperature exceeding limits). While simple and easy to implement, these methods struggle to adapt to dynamic changes under complex operating conditions, are insensitive to early, subtle faults, and are easily affected by environmental noise, resulting in a high false alarm rate. For example, in hydraulic pump pressure fluctuation scenarios, fixed thresholds cannot distinguish between normal operating condition fluctuations and anomalies caused by early wear.

[0004] Statistical models and shallow machine learning methods, such as Principal Component Analysis (PCA) and Hidden Markov Models (HMMs), achieve fault detection through dimensionality reduction or state transition modeling. However, these methods assume that the data follows a static distribution and have limited ability to represent nonlinear, multimodal time-series features. For example, when PCA fuses multi-sensor data in excavators, it may ignore the coupling characteristics between the hydraulic system and the power system, leading to missed detections of critical fault modes.

[0005] Supervised learning dependency issues: Fault diagnosis models based on Support Vector Machines (SVM) and deep neural networks rely on a large amount of labeled data, but fault samples are scarce and labeling costs are high in industrial scenarios, limiting the model's generalization ability. Furthermore, supervised models struggle to handle unknown fault types (such as leaks in novel hydraulic valves), limiting their practicality in online monitoring.

[0006] Existing unsupervised methods have the following limitations:

[0007] Traditional clustering methods (such as k-means) are sensitive to feature redundancy when directly clustering the original time series data, and lack the ability to extract features at multiple scales. They cannot effectively distinguish between patterns with similar waveforms but different fault types (such as vibration signals of pump cavitation and bearing wear).

[0008] Autoencoders (AE) and variational autoencoders (VAE): Although they can detect anomalies through reconstruction errors, the feature extraction process is susceptible to high-frequency noise interference, and they do not explicitly fuse statistical distribution characteristics (such as skewness and interquartile range), resulting in insufficient discriminative power for early progressive faults.

[0009] Similarity measures based on elastic distance (such as DTW) can alleviate the temporal alignment problem, but the computational complexity increases dramatically with the data dimension, making it difficult to meet the real-time requirements of industrial online diagnostics.

[0010] In summary, the aforementioned existing technologies generally suffer from the following drawbacks:

[0011] Insufficient feature representation capabilities: Existing feature extraction methods (such as TSFresh and Catch22) mostly rely on manual design, making it difficult to adaptively capture distribution characteristics (such as short-term fluctuations and long-term trends) and morphological asymmetry (such as skewed signal distribution caused by faults) within multi-scale intervals.

[0012] Feature representation and single feature scale: Existing feature extraction methods often extract features of fault signals in a single domain (e.g., time domain, frequency domain) within a single scale interval, making it difficult to obtain diverse multi-domain fusion features and reducing the feature's ability to represent faults.

[0013] Redundant features cause serious interference: Traditional unsupervised feature selection methods do not sufficiently suppress redundant features, reducing the separation degree of fault clusters. Summary of the Invention

[0014] The technical problem to be solved by this invention is:

[0015] To address the problems of insufficient feature representation capabilities, single feature scale, and severe interference from redundant features in existing feature extraction methods.

[0016] The technical solution adopted by the present invention to solve the above-mentioned technical problems is as follows:

[0017] This invention provides an unsupervised fault diagnosis method based on multi-scale feature enhancement and adaptive feature selection, comprising the following steps:

[0018] S100. Obtain the time series dataset of industrial equipment operating status, and divide the time series dataset into equal-length time series samples of length n.

[0019] S200. For the multiple time series samples of length n obtained in step S100, perform multi-scale interval division, and extract features on each interval to obtain the original time series feature set of industrial equipment operation status.

[0020] S300, Perform sequence feature enhancement on the industrial equipment operating status time series dataset obtained in step S100, including performing first-order difference transform and wavelet transform on the industrial equipment operating status time series obtained in step S100, extracting the transformed approximation coefficient sequence and detail coefficient sequence respectively for the wavelet transform, and performing multi-scale interval division and feature extraction on the transformed time series as described in step S200 to obtain the industrial equipment operating status first-order difference sequence feature set, approximation coefficient sequence feature set and detail coefficient sequence feature set;

[0021] S400: Perform feature fusion and feature selection on the original time series feature set of industrial equipment operating status obtained in step S200 and the first-order difference sequence feature set, approximation coefficient sequence feature set and detail coefficient sequence feature set of industrial equipment operating status obtained in step S300 based on feature enhancement, to obtain the filtered feature set.

[0022] S500 performs cluster analysis on the filtered feature set and outputs equipment status classification and fault warning signals.

[0023] Further, in step S200, the following are included:

[0024] S210, multi-scale interval division

[0025] For the input time series data S=(s1,s2,…,s…), n ), where n represents the length of the time series sample, s i Representing the i-th time series sample point, i = 1, 2, ..., n, the sequence S = (s1, s2, ..., sn) is divided into two parts by a fixed bisection interval. n ) are divided into There are intervals, each with a length of [missing information]. Actual partitioning depth From the formula Confirmed, where d is the preset depth parameter. This is the floor function;

[0026] An offset interval is generated as a supplement. The offset interval has the same length as the original interval, and its starting position is the original interval's starting position shifted backward by half the original interval's length. For a length of... interval S a:b =(s a ,s a+1 ,…,s b ), where s a With s b Representing interval S respectively a:b The starting and ending points, a and b, represent the interval S respectively. a:b The sample point indices of the start and end points, if Its supplementary offset interval is

[0027] All the fixed bisection intervals and the offset intervals together constitute a set of subintervals;

[0028] S220, Interval Feature Extraction

[0029] Quantile features are extracted from each sub-interval of the multi-scale intervals obtained in step S210. The quantiles are then centered using the interquartile range. Skewness features are extracted to characterize the skewness of the interval data distribution, including...

[0030] Given an interval of length n I = b-a+1 time series interval S a:b =(s a ,s a+1 ,…,s b ), in the time series interval S a:b Location Extract k from q One quantile feature, the number of interval quantiles k q Given the interval length n I And the dimension hyperparameter v is determined,

[0031] Extracting time series interval S a:b =(s a ,s a+1 ,…,s b The interquartile range characteristic of ) is the difference between the 75th percentile and the 25th percentile;

[0032] For the extracted k q Each quantile feature is centered by subtracting the interquartile range feature from the interval of each quantile feature.

[0033] Using formula The skewness characteristic of an interval is calculated to characterize the skewness of the data distribution within that interval, where μ represents the interval S. a:b =(s a ,s a+1 ,…,s b The mean of ) , σ represents the interval S a:b =(s a ,s a+1 ,…,s b The standard deviation of s i Let skew(S) represent the i-th sample point in the time series interval, i = a, a+1, ..., b. a:b ) represents the skewness characteristic of a time series interval;

[0034] The skewness features of all time series samples are normalized, and the normalized skewness features are then supplemented into the quantile feature sequence centered by the interquartile range to obtain the feature set of the time series interval.

[0035] Further, in step S300, the following are included:

[0036] S310. The time series of industrial equipment operating status obtained in step S100 is S = (s1, s2, ..., s...). n Perform the feature extraction described in step S200 to obtain the original time series feature set F. O ;

[0037] S320. The time series of industrial equipment operating status obtained in step S100 is S = (s1, s2, ..., s...). n Perform a first-order difference transformation to obtain the first-order difference sequence ΔS = (Δs1, Δs2, ..., Δs...). n-1 ), Δs i =s i+1 -s i For i = 1, 2, ..., n-1, perform feature extraction as described in step S200 on the first-order difference sequence to obtain the first-order difference sequence feature set F. D ;

[0038] S330, wavelet transform driven feature enhancement:

[0039] A predefined set of r wavelet bases b1, b2, ..., b1, b2, ..., b3, r3, r4, r5, r6, r7, r8, r9, r1, r1, r2, ..., r1, r2, r3, r4, r6, r7, r8, r9, r1, r9, r1, r2, ... ...1, r2, r3, r9, r1, r2 r The time series of industrial equipment operating status obtained in step S100 is S = (s1, s2, ..., s...). n Using wavelet bases b1, b2, ..., b r Discrete wavelet transform is performed to extract the approximate coefficient sequence and the detail coefficient sequence. Correlation analysis is conducted between the approximate coefficient sequence and the original time series S and the first-order difference sequence ΔS, respectively. The information entropy of the approximate coefficient sequence is calculated, and p wavelet bases are selected. p≤r; By performing correlation analysis on the detail coefficient sequence with the original time series S and the first-order difference sequence ΔS, and calculating the information entropy of the detail coefficient sequence, q wavelet bases are selected.

[0040] Using wavelet basis Perform discrete wavelet transform and extract their approximate coefficient sequences respectively. Perform the feature extraction described in step S200 to obtain the approximate coefficient sequence feature set. Using wavelet basis Perform discrete wavelet transform and extract the sequence of detail coefficients respectively. Perform the feature extraction described in step S200 to obtain a detailed coefficient sequence feature set

[0041] Further, in step S330, among the wavelet bases b1, b2, …, b in a predefined set containing r wavelet bases available for discrete wavelet transform r include Daubechies wavelet bases, Symlets wavelet bases, Coiflets wavelet bases, Biorthogonal wavelet bases, and ReconstructedBiorthogonal wavelet bases.

[0042] Further, in step S330, the wavelet bases and are screened as follows:

[0043] In a predefined set containing r wavelet bases b1, b2, …, b r available for discrete wavelet transform, use the i-th wavelet base b i to perform discrete wavelet transform on the industrial equipment operating status time series data set obtained in step S100 respectively, where i = 1, 2, …, r, to obtain an approximation coefficient sequence Perform the feature extraction described in step S200 to obtain an approximation coefficient feature set F i A , and use the Pearson correlation coefficient to perform the first correlation analysis on the approximation coefficient feature set F i A and the original time series feature set F O . After arranging them in ascending order of the correlation coefficient, eliminate the candidate wavelet bases with the lowest r percentage of rankings, where 0% < r < 100%, to obtain a candidate wavelet base set B A1 ;

[0044] For the remaining candidate wavelet bases in the candidate wavelet base set B A1 , calculate the Pearson correlation coefficient of the approximation coefficient feature set F i A corresponding to the remaining candidate wavelet bases and the first-order difference sequence feature set F D again, perform the second correlation analysis, and after arranging them in ascending order of the correlation coefficient, eliminate the candidate wavelet bases with the lowest r percentage of rankings, where 0% < r < 100%, to obtain a candidate wavelet base set B A2 ;

[0045] For the remaining candidate wavelet bases in the candidate wavelet base set B A2 , calculate their corresponding feature sets F [i] ∈{F i A |i ∈ B A2The information entropy of}, and the calculation formula is: where is the i-th dimension feature of the h-th time series sample in the feature matrix F [i] P(·) is the probability function. The candidate wavelet basis set B is selected by arranging in descending order of information entropy A2 The wavelet bases with the top p in the information entropy ranking of the corresponding approximation coefficient feature set in

[0046] In a predefined set of wavelet bases b1, b2,..., b r that can be used for discrete wavelet transform, the i-th wavelet basis b i is used to perform discrete wavelet transform on the industrial equipment operating state time series data set obtained in step S100 respectively, i = 1, 2,..., r, and the detail coefficient sequence is obtained Perform the feature extraction described in step S200 to obtain the detail coefficient feature set F i D , and use the Pearson correlation coefficient to perform the first correlation analysis on the detail coefficient feature set F i D and the original time series feature set F O After arranging in ascending order of the correlation coefficient, the candidate wavelet bases with the lowest r percentage of the ranking are removed, where 0% < r < 100%, and the candidate wavelet basis set B is obtained D1 ;

[0047] For the remaining candidate wavelet bases in the candidate wavelet basis set B D1 , calculate the Pearson correlation coefficient of the detail coefficient feature set F i D corresponding to the remaining candidate wavelet bases and the first-order difference sequence feature set F D again, perform the second correlation analysis, and after arranging in ascending order of the correlation coefficient, remove the candidate wavelet bases with the lowest r percentage of the ranking, where 0% < r < 100%, and the candidate wavelet basis set B is obtained D2 ;

[0048] For the remaining candidate wavelet bases in the candidate wavelet basis set B D2 , calculate the information entropy of their corresponding feature set F [i] ∈{F i D |i ∈ B D2}, and the calculation formula is: where is the i-th dimension feature of the h-th time series sample in the feature matrix F [i] P(·) is the probability function. The wavelet bases with the top q in the information entropy ranking of the corresponding detail coefficient feature set in the candidate wavelet basis set B are selected by arranging in descending order of information entropy D2 ​

[0049] Using wavelet basis Discrete wavelet transform is performed on the industrial equipment operating status time series dataset obtained in step S100 to extract p approximate coefficient sequences. Using wavelet basis Discrete wavelet transform is performed on the industrial equipment operating status time series dataset obtained in step S100 to extract q detail coefficient sequences. For the p approximation coefficient sequences With q detail coefficient sequences Perform the feature extraction described in step S200 to obtain the approximate coefficient sequence feature set. With detail coefficient sequence feature set And splice them together to form F W ',Right now:

[0050] Further, in step S400, the following is included:

[0051] S410, Feature set of approximation coefficient sequence With detail coefficient sequence feature set Concatenate into feature set F W ', Feature set F W The number of features is L;

[0052] S420, based on feature set F W A k-nearest neighbor graph is constructed using the cosine similarity between feature vectors, where k represents the number of nearest neighbors connected to each node in the graph. A weighted Laplacian score is calculated to quantify the discriminative power of the features. Features are ranked according to their weighted Laplacian scores, and the feature with the lowest weighted Laplacian score is selected and retained. The wavelet transform feature set F is obtained from these features. W ;

[0053] S430, Transfer the original time series feature set F O First-order difference sequence feature set F D With wavelet transform feature set F W The features are concatenated to form a feature set F', and the number of features in feature set F' is L'.

[0054] S440. Construct a k-nearest neighbor graph based on the cosine similarity between feature vectors in feature set F', calculate the weighted Laplace score to quantify the discriminative power of the features, rank the features according to the weighted Laplace score, and select the k features with the lowest weighted Laplace score. optimal The feature set F is obtained from the features.

[0055] Furthermore, the feature set F contains k features. optimal It is determined based on a variance distribution dual-threshold mechanism, specifically including:

[0056] Calculate the characteristic matrix The variance values ​​of each eigenvector are calculated and sorted in descending order of variance, where n F The number of features in a feature vector is represented by m, which represents the number of time series samples in the industrial equipment operating status time series dataset obtained in step S100. The dimension of each feature vector is m×1.

[0057] k is determined based on a relative threshold. optimal Value: k optimal Defined as T, where the variance after sorting first falls below the maximum variance. r The position index at %.

[0058] Boundary constraint verification: Determine k optimal / n F Does the ratio exceed the preset lower bound b? L or upper bound b U If the value exceeds the threshold, then k is redefined using an absolute threshold. optimal , making k optimal The variance at a given location must be greater than or equal to the absolute threshold, k. optimal The variance value at position +1 is less than the absolute threshold.

[0059] Final adjustment: Ensure the final selected k optimal Always located in [b L ·n F ,b U ·n F Within the interval, if k optimal L ·n F Then k optimal The value is b L ·n F If k optimal >b U ·n F Then k optimal The value is b U ·n F If k optimal Located in [b L ·n F ,b U ·n F Within the interval, then k optimal Retain the original value to avoid extreme cases.

[0060] Further, in steps S420 and S440, the calculation of the weighted Laplace fraction includes,​

[0061] Constructing a time series graph structure: This involves plotting each time series sample S in the time series collection. i =(s i,1 ,s i,2 ,…,s i,n ) as a node in graph G, where s i,p This represents the p-th data point of the i-th time series sample, where p = 1, 2, ..., n, based on the time series sample S. i =(s i,1 ,s i,2 ,…,s i,n ) and time series samples S j =(s j,1 ,s j,2 ,…,s j,n Cosine similarity of (i≠j) For each node S i Determine its k nearest neighbor time series nodes, and establish the connection edges between each node and its k nearest neighbor nodes in graph G;

[0062] Assign edge weights: If node S i With node S j If a connection is made, its edge is assigned a weight w. ij =w ji w ij Represented by node S i To node S j The direction is the weight of the edge, w ji Represented by node S j To node S i The direction is the weight of the edge, where the weight values ​​are uniformly distributed in the interval [w,1], and higher weights are assigned to node pairs with higher cosine similarity, where 0 < w < 1;

[0063] Construct the weighted matrix P: Define the matrix elements P based on the edge weights and Euclidean distance. ij , where P ij Let be the value of the element at the i-th row and j-th column of the weighted matrix P, and its expression is:

[0064] When S i With S j When connecting

[0065] Where T is a constant, and The normalized time series sample is n, where n is the length of the time series sample.

[0066] Calculate the Laplacian matrix: Define a diagonal matrix The Laplace matrix Q is obtained by Q=DP, where m is the number of time series samples in the time series dataset;

[0067] Eigenvector centering: for the r-th eigenvector Centralized processing is performed to obtain in, Let h = 1, 2, ..., m, representing the eigenvectors after centering.

[0068] Calculate the weighted Laplace score: using the formula Determine the eigenvector F [r] The discriminative ability, weighted Laplace score L r The smaller the value, the more representative the eigenvector F. [r] The stronger the ability to distinguish, the better.

[0069] Further, in step S500, the clustering analysis method is one or more combinations of K-means clustering, hierarchical clustering, density-based clustering, spectral clustering, fuzzy C-means clustering, or affinity propagation clustering.

[0070] This invention provides an unsupervised fault diagnosis system based on multi-scale feature enhancement and adaptive feature selection. The system has program modules corresponding to the above steps, and executes the steps in the above-mentioned unsupervised fault diagnosis method based on multi-scale feature enhancement and adaptive feature selection during runtime.

[0071] Compared with the prior art, the beneficial effects of the present invention are:

[0072] Multi-scale fault feature enhancement: This invention extracts multi-dimensional statistical features of quantiles, interquartile ranges, and skewness, and combines adaptive wavelet basis selection with multi-scale interval enhancement based on the maximum redundancy criterion. This significantly improves the ability to characterize the multi-scale characteristics of fault signals (such as high-frequency transient impacts and low-frequency trend drifts) and enhances the ability to capture complex fault modes of industrial equipment (such as intermittent anomalies and gradual degradation). Compared with traditional methods, the fault detection sensitivity is greatly improved.

[0073] Efficient unsupervised feature selection: This invention introduces a cosine similarity-weighted Laplace score and a variance dual threshold mechanism, which maximizes the discriminative power of features while effectively eliminating redundant features (such as environmental noise-related dimensions), thus solving the problem of balancing information preservation and noise suppression in existing technologies.

[0074] Optimized computational efficiency: The high-dimensional data processing time of this invention is significantly reduced compared to the traditional elastic distance method, meeting the real-time requirements of industrial online monitoring and ensuring the robustness and real-time performance of this invention in real industrial scenarios. Attached Figure Description

[0075] Figure 1 This is a flowchart of an unsupervised fault diagnosis method based on multi-scale feature enhancement and adaptive feature selection according to the present invention.

[0076] Figure 2 The following are scatter plots comparing the accuracy of the present invention with existing methods in the embodiments of the present invention: (a) is a scatter plot comparing the accuracy of the present invention with the ED method; (b) is a scatter plot comparing the accuracy of the present invention with the MSM method; (c) is a scatter plot comparing the accuracy of the present invention with the TSFresh method; (d) is a scatter plot comparing the accuracy of the present invention with the Summary method; (e) is a scatter plot comparing the accuracy of the present invention with the Catch22 method; and (d) is a scatter plot comparing the accuracy of the present invention with the QUANT method.

[0077] Figure 3 This is a graph showing the critical difference in accuracy between the present invention and the ED method, MSM method, QUANT method, TSFresh method, Catch22 method, and Summary method in the embodiments of the present invention;

[0078] Figure 4 This is a radar chart showing the performance of the present invention compared to the MSM method, QUANT method, and Summary method on different types of datasets in embodiments of the present invention;

[0079] Figure 5 This is a graph showing the critical differences between the present invention and the ED method, MSM method, QUANT method, TSFresh method, Catch22 method and Summary method in adjusting the Rand index (ARI) in embodiments of the present invention.

[0080] Figure 6 This is a graph showing the critical differences between the invention and the ED method, MSM method, QUANT method, TSFresh method, Catch22 method, and Summary method in terms of the Normalized Mutual Information (NMI) index in the embodiments of the present invention. Detailed Implementation

[0081] To make the above-mentioned objects, features and advantages of the present invention more apparent and understandable, specific embodiments of the present invention will be described in detail below with reference to the accompanying drawings.

[0082] Specific Implementation Scheme 1: This invention provides an unsupervised fault diagnosis method based on multi-scale feature enhancement and adaptive feature selection, the flowchart of which is as follows: Figure 1 As shown, it includes the following steps:

[0083] S100. Obtain the time series dataset of industrial equipment operating status, and divide the time series dataset into equal-length time series samples of length n.

[0084] S200. For the multiple time series samples of length n obtained in step S100, perform multi-scale interval division, and extract features from each interval to obtain the original time series feature set of industrial equipment operating status. The specific steps include:

[0085] S210, multi-scale interval division

[0086] For the input time series data S=(s1,s2,…,s…), n ), where n represents the length of the time series sample, s i Representing the i-th time series sample point, i = 1, 2, ..., n, the sequence S = (s1, s2, ..., sn) is divided into two parts by a fixed bisection interval. n ) are divided into There are intervals, each with a length of [missing information]. Actual partitioning depth From the formula Confirmed, where d is the preset depth parameter. This is the floor function;

[0087] An offset interval is generated as a supplement. The offset interval has the same length as the original interval, and its starting position is the original interval's starting position shifted backward by half the original interval's length. For a length of... interval S a:b =(s a ,s a+1 ,…,s b ), where s a With s b Representing interval S respectively a:b The starting and ending points, a and b, represent the interval S respectively. a:b The sample point indices of the start and end points, if Its supplementary offset interval is

[0088] All the fixed bisection intervals and the offset intervals together constitute a set of subintervals;

[0089] S220, Interval Feature Extraction

[0090] Quantile features are extracted from each sub-interval of the multi-scale intervals obtained in step S210. The quantiles are then centered using the interquartile range. Skewness features are extracted to characterize the skewness of the interval data distribution, including...

[0091] Given an interval of length n I = b-a+1 time series interval S a:b =(s a ,sa+1 ,…,s b ), in the time series interval S a:b Location Extract k from q One quantile feature, the number of interval quantiles k q Given the interval length n I And the dimension hyperparameter v is determined,

[0092] Extracting time series interval S a:b =(s a ,s a+1 ,…,s b The interquartile range characteristic of ) is the difference between the 75th percentile and the 25th percentile;

[0093] For the extracted k q Each quantile feature is centered by subtracting the interquartile range feature from the interval of each quantile feature.

[0094] Using formula The skewness characteristic of an interval is calculated to characterize the skewness of the data distribution within that interval, where μ represents the interval S. a:b =(s a ,s a+1 ,…,s b The mean of ) , σ represents the interval S a:b =(s a ,s a+1 ,…,s b The standard deviation of s i Let skew(S) represent the i-th sample point in the time series interval, i = a, a+1, ..., b. a:b ) represents the skewness characteristic of a time series interval;

[0095] The skewness features of all time series samples are normalized, and then supplemented with the quantile feature sequence centered using the interquartile range to obtain the feature set of the time series interval.

[0096] S300. Perform sequence feature enhancement on the industrial equipment operating status time series dataset obtained in step S100, including performing first-order difference transform and wavelet transform on the industrial equipment operating status time series obtained in step S100. For the wavelet transform, extract the transformed approximation coefficient sequence and detail coefficient sequence respectively. Perform multi-scale interval division and feature extraction on the transformed time series as described in step S200 to obtain the industrial equipment operating status first-order difference sequence feature set, approximation coefficient sequence feature set, and detail coefficient sequence feature set. The specific steps include:

[0097] S310. The time series of industrial equipment operating status obtained in step S100 is S = (s1, s2, ..., s...). n Perform the feature extraction described in step S200 to obtain the original time series feature set F. O ;

[0098] S320. The time series of industrial equipment operating status obtained in step S100 is S = (s1, s2, ..., s...). n Perform a first-order difference transformation to obtain the first-order difference sequence ΔS = (Δs1, Δs2, ..., Δs...). n-1 ), Δs i =s i+1 -s i For i = 1, 2, ..., n-1, perform feature extraction as described in step S200 on the first-order difference sequence to obtain the first-order difference sequence feature set F. D ;

[0099] S330, wavelet transform driven feature enhancement:

[0100] A predefined set of r wavelet bases b1, b2, ..., b1, b2, ..., b3, r3, r4, r5, r6, r7, r8, r9, r1, r1, r2, ..., r1, r2, r3, r4, r6, r7, r8, r9, r1, r9, r1, r2, ... ...1, r2, r3, r9, r1, r2 r The time series of industrial equipment operating status obtained in step S100 is S = (s1, s2, ..., s...). n Using wavelet bases b1, b2, ..., b r Discrete wavelet transform is performed to extract the approximate coefficient sequence and the detail coefficient sequence. Correlation analysis is conducted between the approximate coefficient sequence and the original time series S and the first-order difference sequence ΔS, respectively. The information entropy of the approximate coefficient sequence is calculated, and p wavelet bases are selected. p≤r; By performing correlation analysis on the detail coefficient sequence with the original time series S and the first-order difference sequence ΔS, and calculating the information entropy of the detail coefficient sequence, q wavelet bases are selected.

[0101] Using wavelet basis Perform discrete wavelet transform and extract their approximate coefficient sequences respectively. Perform the feature extraction described in step S200 to obtain the approximate coefficient sequence feature set. Using wavelet basis Perform discrete wavelet transform and extract the sequence of detail coefficients respectively. Perform the feature extraction described in step S200 to obtain the detail coefficient sequence feature set.

[0102] S400: Perform feature fusion and feature selection on the original time series feature set of industrial equipment operating status obtained in step S200 and the first-order difference sequence feature set, approximation coefficient sequence feature set, and detail coefficient sequence feature set of industrial equipment operating status obtained based on feature enhancement in step S300 to obtain the filtered feature set. The specific steps include:

[0103] S410, Feature set of approximation coefficient sequence With detail coefficient sequence feature set Concatenate into feature set F W ', Feature set F W The number of features is L;

[0104] S420, based on feature set F W A k-nearest neighbor graph is constructed using the cosine similarity between feature vectors, where k represents the number of nearest neighbors connected to each node in the graph. A weighted Laplacian score is calculated to quantify the discriminative power of the features. Features are ranked according to their weighted Laplacian scores, and the feature with the lowest weighted Laplacian score is selected and retained. The wavelet transform feature set F is obtained from these features. W ;

[0105] S430, Transfer the original time series feature set F O First-order difference sequence feature set F D With wavelet transform feature set F W The features are concatenated to form a feature set F', and the number of features in feature set F' is L'.

[0106] S440. Construct a k-nearest neighbor graph based on the cosine similarity between feature vectors in feature set F', calculate the weighted Laplace score to quantify the discriminative power of the features, rank the features according to the weighted Laplace score, and select the k features with the lowest weighted Laplace score. optimal The feature set F is obtained from the features;

[0107] The feature set F contains k features. optimal It is determined based on a variance distribution dual-threshold mechanism, specifically including:

[0108] Calculate the characteristic matrix The variance values ​​of each eigenvector are calculated and sorted in descending order of variance, where n F The number of features in a feature vector is represented by m, which represents the number of time series samples in the industrial equipment operating status time series dataset obtained in step S100. The dimension of each feature vector is m×1.

[0109] k is determined based on a relative threshold. optimal Value: k optimal Defined as T, where the variance after sorting first falls below the maximum variance.r The position index at %.

[0110] Boundary constraint verification: Determine k optimal / n F Does the ratio exceed the preset lower bound b? L or upper bound b U If the value exceeds the threshold, then k is redefined using an absolute threshold. optimal , making k optimal The variance at a given location must be greater than or equal to the absolute threshold, k. optimal The variance value at position +1 is less than the absolute threshold.

[0111] Final adjustment: Ensure the final selected k optimal Always located in [b L ·n F ,b U ·n F Within the interval, if k optimal L ·n F Then k optimal The value is b L ·n F If k optimal >b U ·n F Then k optimal The value is b U ·n F If k optimal Located in [b L ·n F ,b U ·n F Within the interval, then k optimal Retain the original value to avoid extreme cases;

[0112] S500 performs cluster analysis on the filtered feature set and outputs equipment status classification and fault warning signals;

[0113] The cluster analysis includes one or more combinations of the following methods:

[0114] K-means clustering: a distance-based iterative algorithm that partitions data by minimizing the squared error within clusters;

[0115] Hierarchical clustering: merging or splitting clusters layer by layer based on distance or similarity using a dendrogram (tree structure);

[0116] Density-based clustering (DBSCAN): Identifies high-density regions and separates low-density noise, adapting to clusters of arbitrary shapes;

[0117] ​Spectral clustering: After dimensionality reduction using the eigenvectors of the data similarity matrix, it is suitable for non-convex data sets;

[0118] Fuzzy C-Means Clustering (FCM): A soft partitioning method that allows data points to belong to multiple clusters with membership degrees;

[0119] Affinity Propagation Clustering: Automatically determines representative sample points through message passing.

[0120] Preferably, in step S330, among the predefined set of r wavelet bases available for discrete wavelet transform, the wavelet bases b1, b2, …, b r include Daubechies wavelet bases, Symlets wavelet bases, Coiflets wavelet bases, Biorthogonal wavelet bases, and Reconstructed Biorthogonal wavelet bases.

[0121] Preferably, in step S330, the wavelet bases and are screened as follows:

[0122] Among the predefined set of r wavelet bases b1, b2, …, b r available for discrete wavelet transform, using the i-th wavelet basis b i to perform discrete wavelet transform on the industrial equipment operating status time series data set obtained in step S100 respectively, where i = 1, 2, …, r, to obtain the approximation coefficient sequence perform the feature extraction described in step S200 to obtain the approximation coefficient feature set F i A , and use the Pearson correlation coefficient to perform the first correlation analysis on the approximation coefficient feature set F i A and the original time series feature set F O . After arranging in ascending order of the correlation coefficient, eliminate the candidate wavelet bases with the lowest r percentage of rankings, where 0% < r < 100%, to obtain the candidate wavelet basis set B A1 ;

[0123] For the remaining candidate wavelet bases in the candidate wavelet basis set B A1 , calculate the Pearson correlation coefficient of the corresponding approximation coefficient feature set F i A and the first-order difference sequence feature set F D again, perform the second correlation analysis, and after arranging in ascending order of the correlation coefficient, eliminate the candidate wavelet bases with the lowest r percentage of rankings, where 0% < r < 100%, to obtain the candidate wavelet basis set B A2 ;

[0124] For the remaining candidate wavelet bases in the candidate wavelet basis set B A2 calculate their corresponding feature sets F respectively [i] ∈{F i A | i ∈ B A2}, and the calculation formula for the information entropy is: where is the i-th dimensional feature of the h-th time series sample in the feature matrix F [i] , P(·) is the probability function. Arrange them in descending order of information entropy, and select the wavelet bases corresponding to the top p of the information entropy ranking of the approximate coefficient feature sets in the candidate wavelet basis set B A2

[0125] In a predefined set containing r wavelet bases b1, b2,..., b r that can be used for discrete wavelet transform, use the i-th wavelet basis b i to perform discrete wavelet transform on the industrial equipment operating state time series data set obtained in step S100 respectively, i = 1, 2,..., r, and obtain the detail coefficient sequence Perform the feature extraction described in step S200 to obtain the detail coefficient feature set F i D , and use the Pearson correlation coefficient to perform the first correlation analysis on the detail coefficient feature set F i D and the original time series feature set F O . After arranging them in ascending order of the correlation coefficient, eliminate the candidate wavelet bases with the lowest r percentage of ranking, where 0% < r < 100%, and obtain the candidate wavelet basis set B D1 ;

[0126] For the remaining candidate wavelet bases in the candidate wavelet basis set B D1 calculate the Pearson correlation coefficient of the remaining candidate wavelet basis corresponding detail coefficient feature set F i D and the first-order difference sequence feature set F D again, perform the second correlation analysis. After arranging them in ascending order of the correlation coefficient, eliminate the candidate wavelet bases with the lowest r percentage of ranking, where 0% < r < 100%, and obtain the candidate wavelet basis set B D2 D2 ;

[0127] For the remaining candidate wavelet bases in the candidate wavelet basis set B D2 calculate their corresponding feature sets F respectively [i] ∈{F i D | i ∈ B D2}, and the calculation formula for the information entropy is: in, The characteristic matrix F [i] The i-th dimension feature of the h-th time series sample, P(·) is the probability function, arranged in descending order of information entropy, and the candidate wavelet basis set B is selected. D2 The wavelet basis with the highest q-order information entropy of the corresponding detail coefficient feature set is as follows.

[0128] Using wavelet basis Discrete wavelet transform is performed on the industrial equipment operating status time series dataset obtained in step S100 to extract p approximate coefficient sequences. Using wavelet basis Discrete wavelet transform is performed on the industrial equipment operating status time series dataset obtained in step S100 to extract q detail coefficient sequences. For the p approximation coefficient sequences With q detail coefficient sequences Perform the feature extraction described in step S200 to obtain the approximate coefficient sequence feature set. With detail coefficient sequence feature set And splice them together to form F W ',Right now:

[0129] Preferably, in steps S420 and S440, the calculation step of the weighted Laplace fraction includes,

[0130] Constructing a time series graph structure: This involves plotting each time series sample S in the time series collection. i =(s i,1 ,s i,2 ,…,s i,n ) as a node in graph G, where s i,p This represents the p-th data point of the i-th time series sample, where p = 1, 2, ..., n, based on the time series sample S. i =(s i,1 ,s i,2 ,…,s i,n ) and time series samples S j =(s j,1 ,s j,2 ,…,s j,n Cosine similarity of (i≠j) For each node S i Determine its k nearest neighbor time series nodes, and establish the connection edges between each node and its k nearest neighbor nodes in graph G;

[0131] Assign edge weights: If node S i With node S j If a connection is made, its edge is assigned a weight w.ij =w ji w ij Represented by node S i To node S j The direction is the weight of the edge, w ji Represented by node S j To node S i The direction is the weight of the edge, where the weight values ​​are uniformly distributed in the interval [w,1], and higher weights are assigned to node pairs with higher cosine similarity, where 0 < w < 1;

[0132] Construct the weighted matrix P: Define the matrix elements P based on the edge weights and Euclidean distance. ij , where P ij Let be the value of the element at the i-th row and j-th column of the weighted matrix P, and its expression is:

[0133] When S i With S j When connecting

[0134] Where T is a constant, and The normalized time series sample is n, where n is the length of the time series sample.

[0135] Calculate the Laplacian matrix: Define a diagonal matrix The Laplace matrix Q is obtained by Q=DP, where m is the number of time series samples in the time series dataset;

[0136] Eigenvector centering: for the r-th eigenvector Centralized processing is performed to obtain in, Let h = 1, 2, ..., m, representing the eigenvectors after centering.

[0137] Calculate the weighted Laplace score: using the formula Determine the eigenvector F [r] The discriminative ability, weighted Laplace score L r The smaller the value, the more representative the eigenvector F. [r] The stronger the ability to distinguish, the better.

[0138] Preferably, the preset depth parameter d is set to 6, the dimensional hyperparameter v is set to 4, the number of wavelet bases p is set to 3, the number of wavelet bases q is set to 3, and the relative threshold T is set to... r The value of % is 1%, and the b L The value is 35%, and the b UThe value is 65%, the absolute threshold is 0.002, the elimination ratio r is 50%, and the number of k-nearest neighbors in the k-nearest neighbor graph constructed based on the cosine similarity between feature vectors in the feature set is 10.

[0139] Specific Implementation Scheme Two: This invention provides an unsupervised fault diagnosis method based on multi-scale feature enhancement and adaptive feature selection, comprising the following modules:

[0140] The data acquisition module is used to acquire multimodal sensor data of the industrial equipment to be diagnosed in real time and construct a time series dataset according to the acquisition time sequence.

[0141] The data processing module is used for multi-scale feature extraction and enhancement of time series data, including:

[0142] a. Multi-scale sub-intervals are generated by dividing the data using fixed bisection intervals and offset intervals;

[0143] b. Extract quantile, interquartile range (IQR) and normalized skewness features for each sub-interval, and fuse the time-frequency domain features through discrete wavelet transform;

[0144] The feature selection module, used to filter features based on weighted Laplace scores, includes:

[0145] a. Construct a k-nearest neighbor graph based on cosine similarity and calculate the weight matrix using Euclidean distance;

[0146] b. Quantify the distinguishing ability of features using weighted Laplace scores;

[0147] c. Adaptively determine the number of features to be retained based on a variance distribution dual threshold mechanism;

[0148] The unsupervised clustering module performs clustering analysis using k-means or other methods on the above feature set.

[0149] The fault diagnosis module is used to integrate clustering results, match them with a preset fault mode library, and output equipment status diagnosis reports and early warning signals.

[0150] The communication module enables real-time data transmission between modules through edge computing nodes, supporting cloud synchronization and remote monitoring.

[0151] The other combinations and connections in this implementation scheme are the same as in Specific Implementation Scheme 1.

[0152] When this invention is used in specialized equipment in the mechanical field, the operating status information includes vibration amplitude, sound spectrum characteristics, and rotational speed of mechanical transmission components; when this invention is used in special equipment in the chemical field, the operating status information includes three-dimensional temperature gradient distribution inside the reactor, fluid pressure fluctuation values ​​in pipelines, and material phase transformation characteristic parameters; when this invention is used in general equipment in the energy field, the operating status information includes turbine shaft vibration modes, generator winding temperature distribution matrix, and harmonic component characteristic values ​​of transmission lines.

[0153] Simulation Experiment 1

[0154] To verify the effectiveness of the method of this invention, a computer simulation classification comparison experiment was conducted on 112 datasets on the UCR open-source dataset. Each dataset used a default training / test set split, and each sample was standardized. Clustering accuracy (ACC) was used as the performance metric for evaluation. ACC is defined as the number of correctly predicted clusters divided by the total number of samples. A Wilcoxon signed-rank test with Holm correction and a p-value of 0.05 was used to measure the significance of the differences between the different algorithms. The results are shown in Table 1.

[0155] Existing methods include using the Move-Split-Merge (MSM) elastic distance metric and a k-means clusterer using the Euclidean distance metric (ED) with clusterer parameters set to random initialization, a maximum of 300 iterations, 10 restarts, and mean averaging. The results for both methods are derived from the reference C. Holder, M., Middlehurst, and A. Bagnall, “A review and evaluation of elastic distance functions for time series clustering,” Knowledge Information Systems, vol. 66, no. 2, pp. 765–809, Feb. 2024.

[0156] And existing feature-based methods:

[0157] a. A feature extraction algorithm based on 22 typical time series characteristics (22 Canonical Time-series Characteristics, abbreviated as catch22), the method of which comes from the reference CHLubba, SS Sethi, P. Knaute, et al., “catch22: CAnonical Time-series CHaracteristics Selected through highly comparative time-series analysis,” Data Mining Knowledge Discovery, vol.33, no.6, pp.1821–1852, Nov.2019;

[0158] b. A time series feature extraction method based on scalable hypothesis tests (TSFresh), which is derived from the reference M. Christ, N. Braun, J. Neuffer, and AWKempa-Liehr, “Time Series Feature Extraction on basis of Scalable Hypothesis tests (tsfresh – A Python package),” Neurocomputing, vol. 307, pp. 72–77, Sep. 2018;

[0159] c. A feature extraction method based on the minimalist interval method (QUANT), which is derived from the reference A. Dempster, DFSchmidt, and GIWebb, “QUANT: a minimalist interval method for time series classification,” Data Mining Knowledge Discovery, vol.38, pp.2377–2402, May 2024;

[0160] d. Summary Statistics-based Time Series Feature Extraction (Summary): This method extracts seven basic statistical features of the time series (including minimum, maximum, 0.1 quantile, 0.25 quantile, 0.5 quantile, 0.75 quantile, and 0.9 quantile).

[0161] All four existing feature-based methods mentioned above use a k-means clusterer after feature extraction. The clusterer parameters are set to random initialization, a maximum of 300 iterations, 10 restarts, mean averaging, and Euclidean distance metric.

[0162] Table 1

[0163]

[0164]

[0165]

[0166]

[0167] The scatter plot showing the accuracy of the method of this invention compared with existing methods on 112 UCR datasets is shown below. Figure 2 As shown, points falling above / above / below the diagonal represent the accuracy of the method of the present invention on the dataset that is inferior to / equal to / superior to the accuracy of the prior art method, respectively. Here, W represents the number of datasets where the classification accuracy of the method is superior to the comparison method, T represents the number of datasets where the classification accuracy of the method is equal to the comparison method, and L represents the number of datasets where the classification accuracy of the method is inferior to the comparison method.

[0168] From Table 1 and Figure 2 The results show that the method of the present invention not only outperforms existing methods in accuracy on most datasets, but also significantly outperforms them in accuracy on some datasets.

[0169] Figure 3 The graph shows the critical difference in accuracy between the method of this invention and existing methods. Methods with no statistically significant differences are connected by thick black lines. The numbers represent the average ranking of the method of this invention and existing methods on 112 UCR datasets. The method of this invention achieves the highest average ranking on the accuracy metric, demonstrating its superiority.

[0170] Simulation Experiment 2

[0171] To test the applicability of the method of this invention, it was compared with existing methods such as QUANT, Summary, and MSM on nine tasks of the UCR open-source dataset. The performance radar chart is attached. Figure 4 As shown, the method of the present invention exhibits the best performance in five categories of tasks: sensor, spectrum, device, human activity recognition, and other data; and ranks second in three categories of tasks: simulation, electrocardiogram, and image data. The method of the present invention is applicable to a wide range of problems.

[0172]

[0173] Simulation Experiment 3

[0174] To further verify the effectiveness of the method of this invention, a computer simulation classification comparison experiment was conducted on 112 datasets on the UCR open-source dataset, comparing the method of this invention with existing methods. Each dataset used the default training / test set split, and each sample was standardized. The Adjusted Rand Index (ARI) was used as the performance metric for evaluation. The Rand Index (RI) is calculated by dividing the number of logarithmic agreements between the true labels and the model-generated labels by the total number of logarithmic agreements. However, one limitation of the Rand Index is that the score is amplified as the number of clusters increases. The Adjusted Rand Index (ARI) compensates for this limitation by using the expected score of a purely random model. Furthermore, the Wilcoxon signed-rank test with Holm correction and a p-value of 0.05 was used to measure the significance of the differences between the different algorithms. The results are shown in Table 2.

[0175] Table 2

[0176]

[0177]

[0178]

[0179] Figure 5 The graph shows the critical difference between the method of this invention and existing methods in adjusting the Rand Index (ARI). Methods with no statistically significant differences are connected by thick black lines. The numbers represent the average rankings calculated by the method of this invention and existing methods on 112 UCR datasets. The method of this invention achieves the highest average ranking in adjusting the Rand Index (ARI), demonstrating its superiority.

[0180] Simulation Experiment 4

[0181] To further verify the effectiveness of the method of this invention, a computer simulation classification comparison experiment was conducted on 112 datasets on the UCR open-source dataset, comparing the method of this invention with existing methods. Each dataset used a default training / test set split, and each sample was standardized. Normalized Mutual Information (NMI) was used as the performance metric for evaluation. Mutual Information (MI) is a function based on entropy that measures the consistency between two clusters or between cluster labels and the true labels. Normalized Mutual Information (NMI) scales MI to the [0,1] interval. Furthermore, a Wilcoxon signed-rank test with Holm correction and a p-value of 0.05 was used to measure the significance of differences between different algorithms. The results are shown in Table 3.

[0182] Table 3

[0183]

[0184]

[0185]

[0186] Figure 6 This graph shows the critical difference between the method of this invention and existing methods on the Normalized Mutual Information (NMI) metric. Methods with no statistically significant differences are connected by thick black lines. The numbers represent the average rankings calculated by the method of this invention and existing methods on 112 UCR datasets. The method of this invention achieves the highest average ranking on the Normalized Mutual Information (NMI) metric, demonstrating its superiority.

[0187] While the present invention has been disclosed above, its scope of protection is not limited thereto. Those skilled in the art can make various changes and modifications without departing from the spirit and scope of the present invention, and all such changes and modifications will fall within the scope of protection of the present invention.

Claims

1. An unsupervised fault diagnosis method based on multi-scale feature enhancement and adaptive feature selection, characterized in that, Includes the following steps: S100. Obtain the time series dataset of industrial equipment operating status, and divide the time series dataset into equal-length time series samples of length n. S200. For the multiple time series samples of length n obtained in step S100, perform multi-scale interval division, and extract features on each interval to obtain the original time series feature set of industrial equipment operation status. This includes S210 multi-scale interval partitioning and S220 interval feature extraction. In step S220, interval feature extraction, Quantile features are extracted from each sub-interval of the multi-scale intervals obtained in step S210. The quantiles are then centered using the interquartile range. Skewness features are extracted to characterize the skewness of the interval data distribution, including... Given an interval of length n I =b-a+1 time series interval In the time series interval Location Extract k from q One quantile feature, the number of interval quantiles k q Given the interval length n I And the dimension hyperparameter v is determined, ; Extracting time series intervals The interquartile range is a characteristic of the interquartile range, which is the difference between the 75th percentile and the 25th percentile. For the extracted k q Each quantile feature is centered by subtracting the interquartile range feature from the interval of each quantile feature. Using formula The skewness characteristics of the intervals are calculated to characterize the skewness of the data distribution within the intervals, where, Representing an interval The mean, Representing an interval standard deviation, s i Let skew(S) represent the i-th sample point in the time series interval, i = a, a+1, ..., b. a:b ) represents the skewness characteristic of the time series interval; a and b represent the intervals respectively. The sample point indexes of the start and end points; The skewness features of all time series samples are normalized, and the normalized skewness features are supplemented into the quantile feature sequence centered by the interquartile range to obtain the feature set of the time series interval. S300, Perform sequence feature enhancement on the industrial equipment operating status time series dataset obtained in step S100, including performing first-order difference transform and wavelet transform on the industrial equipment operating status time series obtained in step S100, extracting the transformed approximation coefficient sequence and detail coefficient sequence respectively for the wavelet transform, and performing multi-scale interval division and feature extraction on the transformed time series as described in step S200 to obtain the industrial equipment operating status first-order difference sequence feature set, approximation coefficient sequence feature set and detail coefficient sequence feature set; S400: Perform feature fusion and feature selection on the original time series feature set of industrial equipment operating status obtained in step S200 and the first-order difference sequence feature set, approximation coefficient sequence feature set and detail coefficient sequence feature set of industrial equipment operating status obtained in step S300 based on feature enhancement, to obtain the filtered feature set. S500 performs cluster analysis on the filtered feature set and outputs equipment status classification and fault warning signals.

2. The unsupervised fault diagnosis method based on multi-scale feature enhancement and adaptive feature selection according to claim 1, characterized in that: Step S200 includes, S210, multi-scale interval division For the input time series data Where n represents the length of the time series sample, s i Representing the i-th time series sample point, i=1,2,…,n, the sequence is divided using a fixed bisection interval method. Divided into There are intervals, each with a length of [missing information]. Actual partitioning depth From the formula Confirmed, where d is the preset depth parameter. This is the floor function; An offset interval is generated as a supplement. The offset interval has the same length as the original interval, and its starting position is the original interval's starting position shifted backward by half the original interval's length. For a length of... interval , where s a With s b Representing intervals The starting and ending points, a and b, represent the intervals respectively. The sample point indices of the start and end points, if Its supplementary offset interval is , ; All the fixed bisection intervals and the offset intervals together constitute a set of subintervals.

3. The unsupervised fault diagnosis method based on multi-scale feature enhancement and adaptive feature selection according to claim 1, characterized in that: Step S300 includes, S310. The time series of industrial equipment operating status obtained in step S100. Perform the feature extraction described in step S200 to obtain the original time series feature set. ; S320. The time series of industrial equipment operating status obtained in step S100. Perform a first-order difference transform to obtain a first-order difference sequence. , For i=1,2,…,n-1, perform the feature extraction described in step S200 on the first-order difference sequence to obtain the first-order difference sequence feature set. ; S330, wavelet transform driven feature enhancement: A predefined set of r wavelet bases that can be used for discrete wavelet transform is defined. The time series of industrial equipment operating status obtained in step S100 Using wavelet basis respectively Perform discrete wavelet transform to extract the approximate coefficient sequence and the detail coefficient sequence. Then, compare the approximate coefficient sequence with the original time series. First-order difference sequence Correlation analysis was performed, and the information entropy of the approximate coefficient sequence was calculated. P wavelet bases were then selected. , By comparing the detail coefficient sequence with the original time series... First-order difference sequence Correlation analysis was performed, and the information entropy of the detail coefficient sequence was calculated. q wavelet bases were then selected. , ; Using wavelet basis Perform discrete wavelet transform and extract their approximate coefficient sequences respectively. Perform the feature extraction described in step S200 to obtain the approximate coefficient sequence feature set. Using wavelet bases Perform discrete wavelet transform and extract the sequence of detail coefficients respectively. Perform the feature extraction described in step S200 to obtain the detail coefficient sequence feature set. .

4. The unsupervised fault diagnosis method based on multi-scale feature enhancement and adaptive feature selection according to claim 3, characterized in that: In step S330, the predefined set of wavelet bases containing r wavelet bases that can be used for discrete wavelet transform is used. This includes the Daubechies wavelet basis, Symlets wavelet basis, Coiflets wavelet basis, Biorthogonal wavelet basis, and Reconstructed Biorthogonal wavelet basis.

5. The unsupervised fault diagnosis method based on multi-scale feature enhancement and adaptive feature selection according to claim 4, characterized in that: In step S330, wavelet basis and The screening steps are as follows: In a predefined set of r wavelet bases that can be used for discrete wavelet transform In the set, using the i-th wavelet basis Discrete wavelet transforms are performed on the industrial equipment operating status time series dataset obtained in step S100. The approximate coefficient sequence is obtained. Perform the feature extraction described in step S200 to obtain the approximate coefficient feature set. Using Pearson correlation coefficient to approximate coefficient feature set Compared with the original time series feature set After performing the first correlation analysis and sorting the candidates in ascending order of correlation coefficient, the candidate wavelet bases with the lowest rank (r) are removed, where 0% < r < 100%, resulting in the candidate wavelet base set B. A1 ; For the candidate wavelet basis set B A1 For the remaining candidate wavelet bases, calculate the approximate coefficient feature set corresponding to the remaining candidate wavelet bases again. Features of first-order difference sequences The Pearson correlation coefficient was used to perform a second correlation analysis. After sorting the candidates by correlation coefficient in ascending order, the candidate wavelet bases with the lowest r-proportion (0% < r < 100%) were removed, resulting in the candidate wavelet base set B. A2 ; For the candidate wavelet basis set B A2 For the remaining candidate wavelet bases, calculate their corresponding feature sets. The information entropy is calculated using the following formula: ,in, Characteristic matrix The i-th dimension feature of the h-th time series sample, P(⋅) is the probability function, arranged in descending order of information entropy, and the candidate wavelet basis set B is selected. A2 The wavelet basis with the highest approximation entropy of the feature set corresponding to the approximate coefficients in the middle is ranked p. m represents the number of time series samples in the industrial equipment operating status time series dataset obtained in step S100. In a predefined set of r wavelet bases that can be used for discrete wavelet transform In the set, using the i-th wavelet basis Discrete wavelet transforms are performed on the industrial equipment operating status time series dataset obtained in step S100. , thus obtaining the detail coefficient sequence Perform the feature extraction described in step S200 to obtain the detail coefficient feature set. Using Pearson correlation coefficient to analyze the feature set of detail coefficients Compared with the original time series feature set After performing the first correlation analysis and sorting the candidates in ascending order of correlation coefficient, the candidate wavelet bases with the lowest rank (r) are removed, where 0% < r < 100%, resulting in the candidate wavelet base set B. D1 ; For the candidate wavelet basis set B D1 For the remaining candidate wavelet bases, calculate the detail coefficient feature set corresponding to the remaining candidate wavelet bases again. Features of first-order difference sequences The Pearson correlation coefficient was used to perform a second correlation analysis. After sorting the candidates by correlation coefficient in ascending order, the candidate wavelet bases with the lowest r-proportion (0% < r < 100%) were removed, resulting in the candidate wavelet base set B. D2 ; For the candidate wavelet basis set B D2 For the remaining candidate wavelet bases, calculate their corresponding feature sets. The information entropy is calculated using the following formula: ,in, Characteristic matrix The i-th dimension feature of the h-th time series sample, P(⋅) is the probability function, arranged in descending order of information entropy, and the candidate wavelet basis set B is selected. D2 The wavelet basis with the highest q-order information entropy of the corresponding detail coefficient feature set is as follows. ; Using wavelet basis Discrete wavelet transform is performed on the industrial equipment operating status time series dataset obtained in step S100 to extract p approximate coefficient sequences. Using wavelet basis Discrete wavelet transform is performed on the industrial equipment operating status time series dataset obtained in step S100 to extract q detail coefficient sequences. For the p approximation coefficient sequences With q detail coefficient sequences Perform the feature extraction described in step S200 to obtain the approximate coefficient sequence feature set. With detail coefficient sequence feature set And splice them together ,Right now: .

6. The unsupervised fault diagnosis method based on multi-scale feature enhancement and adaptive feature selection according to claim 5, characterized in that: Step S400 includes, S410, Feature set of approximation coefficient sequence With detail coefficient sequence feature set splicing into a feature set Feature set The number of features is L; S420, based on feature sets A k-nearest neighbor graph is constructed using the cosine similarity between feature vectors, where k represents the number of nearest neighbors connected to each node in the graph. A weighted Laplacian score is calculated to quantify the discriminative power of the features. Features are ranked according to their weighted Laplacian scores, and the feature with the lowest weighted Laplacian score is selected and retained. The wavelet transform feature set is obtained from these features. ; S430, Convert the original time series feature set First-order difference sequence feature set With wavelet transform feature set splicing into a feature set Feature set The number of features is ; S440, based on feature sets A k-nearest neighbor graph is constructed using the cosine similarity between feature vectors. A weighted Laplacian score is calculated to quantify the discriminative power of the features. Features are ranked according to their weighted Laplacian scores, and the feature with the lowest weighted Laplacian score is selected and retained. The feature set is obtained from the features. .

7. The unsupervised fault diagnosis method based on multi-scale feature enhancement and adaptive feature selection according to claim 6, characterized in that: The feature set Number of features It is determined based on a variance distribution dual-threshold mechanism, specifically including: Calculate the characteristic matrix The variance values ​​of each eigenvector are calculated and sorted in descending order of variance. Let m represent the number of features of a feature vector, and m represent the number of time series samples in the industrial equipment operating status time series dataset obtained in step S100. The dimension of each feature vector is... ; k is determined based on a relative threshold. optimal Value: k optimal Defined as T, where the variance after sorting first falls below the maximum variance. r The position index at %. Boundary constraint verification: Determine k optimal / n F Does the ratio exceed the preset lower bound b? L or upper bound b U If the value exceeds the threshold, then k is redefined using an absolute threshold. optimal , making k optimal The variance at a given location must be greater than or equal to the absolute threshold, k. optimal The variance value at position +1 is less than the absolute threshold. Final adjustment: Ensure the final selected k optimal Always located in [b L ⋅n F ,b U ⋅n F Within the interval, if k optimal L ⋅n F Then k optimal The value is b L ⋅n F If k optimal >b U ⋅n F Then k optimal The value is b U ⋅n F If k optimal Located in [b L ⋅n F ,b U ⋅n F Within the interval, then k optimal Retain the original value to avoid extreme cases.​ 8. The unsupervised fault diagnosis method based on multi-scale feature enhancement and adaptive feature selection according to claim 6, characterized in that: In steps S420 and S440, the calculation of the weighted Laplace fraction includes, Constructing a time series graph structure: This involves plotting each time series sample in the time series collection. As a node in graph G, s i,p This represents the p-th data point of the i-th time series sample, where p = 1, 2, ..., n, based on the time series sample. With time series samples Cosine similarity of (i≠j) For each node Determine its k nearest neighbor time series nodes, and establish the connection edges between each node and its k nearest neighbor nodes in graph G; Assign edge weights: If node With nodes If a connection is made, its edges are assigned weights. w ij Represented by nodes To node The direction is the weight of the edge, w ji Represented by nodes To node The direction is the weight of the edge, where the weight value is in the interval [0, 1]. The nodes are uniformly distributed within the matrix, with higher weights assigned to nodes with higher cosine similarity. ; Construct the weighted matrix P: Define the matrix elements P based on the edge weights and Euclidean distance. ij , where P ij Let be the value of the element at the i-th row and j-th column of the weighted matrix P, and its expression is: Where T is a constant, and The normalized time series sample is n, where n is the length of the time series sample. Calculate the Laplacian matrix: Define a diagonal matrix And the Laplace matrix Q is obtained by Q=D−P, where m is the number of time series samples in the time series dataset; Eigenvector centering: for the r-th eigenvector Centralized processing is performed to obtain ,in, , Let h = 1, 2, ..., m, be the eigenvectors after centering. Calculate the weighted Laplace score: using the formula Determine the eigenvector The ability to discriminate, weighted Laplace score The smaller the value, the more representative the feature vector. The stronger the ability to distinguish, the better.

9. The unsupervised fault diagnosis method based on multi-scale feature enhancement and adaptive feature selection according to claim 1, characterized in that: In step S500, the clustering analysis method is one or more combinations of K-means clustering, hierarchical clustering, density-based clustering, spectral clustering, fuzzy C-means clustering, or affinity propagation clustering.

10. An unsupervised fault diagnosis system based on multi-scale feature enhancement and adaptive feature selection, characterized in that: The system has a program module corresponding to the steps of any one of the claims 1-9 above, and executes the steps in the above-described unsupervised fault diagnosis method based on multi-scale feature enhancement and adaptive feature selection when running.

Images