Method and system for detecting the health of a packaging machine installation

Abstract

This application discloses a method and system for detecting the health status of packaging machinery, relating to the field of tobacco technology. The method constructs a complete equipment health assessment system using a support vector data description algorithm. First, the weighted aggregated temporal feature vectors are mapped to a high-dimensional feature space, and the Euclidean distance matrix of the samples is calculated using a Gaussian radial basis kernel function. Then, a quadratic programming optimization problem is solved to determine the center and radius of a hypersphere representing the normal equipment state. Subsequently, the distance from the sample to the center is calculated to quantify the degree of deviation from the normal state. Finally, based on hypersphere radius normalization, a multi-dimensional equipment health score set is obtained. This scheme can achieve quantitative assessment and anomaly detection of equipment health status, providing a scientific basis for equipment maintenance decisions.

Description

Technical Field

[0001] This application relates to the field of tobacco technology, and in particular to a method and system for detecting the health of packaging machinery. Background Technology

[0002] Packaging machines are key equipment in the tobacco production process, primarily used for packaging, sealing, and labeling cigarette products. Modern packaging machines typically employ complex mechanical transmission and electrical control systems, including multiple subsystems such as feeding mechanisms, forming mechanisms, sealing mechanisms, and detection mechanisms. During long-term high-speed operation, key components such as bearings, gears, motors, and transmission chains gradually experience wear, fatigue, and loosening, leading to potential malfunctions, decreased equipment performance, unstable product quality, and even machine shutdown. Real-time monitoring of the health of packaging machine equipment is crucial for ensuring product quality, reducing failure rates, and minimizing maintenance costs.

[0003] Currently, health monitoring of packaging machinery mainly relies on manual inspections and periodic maintenance. This approach suffers from problems such as long inspection cycles, strong subjectivity in human judgment, and inability to detect early faults in a timely manner. Traditional methods based on single parameter thresholds cannot comprehensively reflect the overall health status of the equipment and cannot accurately assess the trend of equipment performance degradation. Existing detection methods lack effective feature extraction and fusion techniques when processing multi-source heterogeneous signals, making it difficult to uncover deep-seated fault characteristics from multi-dimensional signals such as vibration, temperature, current, and speed. This results in low detection accuracy and a high false alarm rate, failing to meet the needs of modern production lines for real-time and accurate monitoring of equipment health status. Summary of the Invention

[0004] The main purpose of this application is to provide a health detection method and system for packaging machine equipment, in order to solve the problems that traditional soft measurement models in the prior art are difficult to effectively model the internal mechanism of the drying process, and that most existing deep learning methods lack interpretability and are difficult to provide a reliable basis for process analysis and operation adjustment.

[0005] To achieve the above objectives, this application provides the following technical solution: A method for detecting the health status of packaging machine equipment, the method comprising: Step S1: The vibration signal, temperature signal, current signal, and speed signal of the packaging machine are collected in real time through external sensors, and the collected raw signals are denoised and normalized. Step S2: Input the preprocessed vibration signal, preprocessed temperature signal, preprocessed current signal, and preprocessed rotation speed signal into the empirical mode decomposition algorithm to decompose them into several intrinsic mode function components and several residual components. Step S3: Perform Hilbert-Huang transform on each intrinsic mode function component and each residual component to calculate the time-frequency characteristics, so as to extract the instantaneous frequency, instantaneous amplitude, and marginal spectral energy density. Step S4: Input all instantaneous frequencies, all instantaneous amplitudes, and all marginal spectral energy densities into the Long Short-Term Memory network for temporal feature learning. Update the cell state and hidden state through the gating mechanism of forget gate, input gate, and output gate to obtain the temporal feature vector. Step S5: Input the temporal feature vector into the attention mechanism to calculate the weight coefficient of each time step, and calculate the weighted aggregated temporal feature vector by weighted summation and aggregating the feature information of key time steps; Step S6: Input the weighted aggregated temporal feature vector into the support vector data description algorithm to construct a hypersphere model of the device health status. Calculate the distance from the sample to the center of the hypersphere by mapping to a high-dimensional feature space through a kernel function to obtain the device health score.

[0006] Beneficial effects of steps S1 to S6: By organically combining technologies such as multi-sensor data acquisition, signal decomposition, time-frequency feature extraction, time-series learning, attention mechanism weighting, and anomaly detection, a complete health monitoring solution for packaging machines is formed. This solution can monitor the operating status of packaging machines in real time, automatically extract deep fault features from multi-source heterogeneous signals, effectively assess the health level of the equipment, and achieve early fault warning. It solves the problems of long detection cycles, strong subjectivity, and low accuracy associated with traditional methods that rely on manual inspection and single-parameter threshold judgment, significantly improving the real-time performance, accuracy, and reliability of equipment health status monitoring. Specifically, step S1 denoises and normalizes vibration, temperature, current, and speed signals to eliminate noise interference and dimensional differences, providing high-quality basic data for subsequent analysis; step S2 uses an empirical mode decomposition algorithm to adaptively decompose the signal into feature components of different scales, improving the separability of signal features; and step S3 extracts instantaneous frequency, instantaneous amplitude, and marginal spectral energy density through Hilbert-Huang transform, achieving accurate characterization of the signal's time-frequency features. Step S4 uses a long short-term memory network to capture the temporal dependencies in the feature sequence and obtain a feature vector containing deep temporal information; Step S5 calculates weight coefficients and aggregates key features through an attention mechanism to improve the effectiveness of feature representation; Step S6 constructs a hypersphere model of the device health status based on the support vector data description algorithm and calculates the health score to achieve quantitative assessment and anomaly detection of the device health status.

[0007] As a further improvement to this application, step S1 involves real-time acquisition of vibration signals, temperature signals, current signals, and rotation speed signals of the packaging machine via external sensors, and denoising and normalizing the acquired raw signals, including: Step S1.1: Acquire vibration acceleration signal through accelerometer, temperature signal through thermocouple sensor, current signal through Hall sensor, and rotation speed signal through encoder to obtain original vibration signal, original temperature signal, original current signal, and original rotation speed signal; Step S1.2: Input the original vibration signal, the original temperature signal, the original current signal and the original rotation speed signal into the wavelet threshold denoising algorithm respectively, and decompose and reconstruct them through the sym4 wavelet basis function to obtain the denoised vibration signal, the denoised temperature signal, the denoised current signal and the denoised rotation speed signal. Step S1.3: The denoised vibration signal, the denoised temperature signal, the denoised current signal, and the denoised rotation speed signal are normalized by Z-score standardization to obtain normalized vibration signal, normalized temperature signal, normalized current signal, and normalized rotation speed signal, respectively. Step S1.4: Align the normalized vibration signal, the normalized temperature signal, the normalized current signal, and the normalized rotational speed signal according to the timestamp, and fill in the missing data points using a linear interpolation method; Step S1.5: After filling in the missing data points, the signal is filtered by a moving average algorithm to eliminate high-frequency jitter interference, resulting in preprocessed vibration signal, preprocessed temperature signal, preprocessed current signal, and preprocessed rotation speed signal.

[0008] Beneficial effects of steps S1.1 to S1.5: By organically combining multi-sensor collaborative acquisition, wavelet threshold denoising, Z-score normalization, timestamp alignment and linear interpolation imputation, and moving average filtering, a complete signal preprocessing technology solution is formed. This solution can effectively eliminate noise interference, dimensional differences, time deviations, and data gaps in the original signal, significantly improving signal quality and consistency. It provides high-quality basic data for subsequent signal decomposition, feature extraction, and health assessment, solving problems such as residual noise, data mismatch, and poor feature extractability in traditional preprocessing methods when processing multi-source heterogeneous signals. The process involves several steps: Step S1.1 involves acquiring vibration acceleration, temperature, current, and rotational speed signals using an accelerometer, thermocouple, Hall effect sensor, and encoder, respectively, to comprehensively obtain key physical quantities reflecting the equipment's status during operation. Step S1.2 involves inputting the acquired raw signals into a wavelet threshold denoising algorithm, decomposing and reconstructing them using the sym4 wavelet basis function to effectively remove high-frequency noise and random interference, thus improving the signal-to-noise ratio. Step S1.3 involves normalizing the denoised signals using Z-score standardization to eliminate dimensional differences and inconsistent numerical ranges between signals from different sensors, ensuring that all signals are within a uniform numerical scale. Step S1.4 involves aligning the normalized signals according to timestamps and filling in missing data points using linear interpolation to ensure the synchronization and integrity of multiple signals in the time dimension. Step S1.5 involves applying a moving average filtering algorithm to the signals after filling in missing data points to eliminate high-frequency jitter interference and smooth instantaneous fluctuations in the signals, resulting in preprocessed vibration, temperature, current, and rotational speed signals.

[0009] As a further improvement of this application, step S2 involves inputting the preprocessed vibration signal, preprocessed temperature signal, preprocessed current signal, and preprocessed rotational speed signal into an empirical mode decomposition algorithm, which decomposes them into several intrinsic mode function components and several residual components, including: Step S2.1: Input the preprocessed vibration signal, preprocessed temperature signal, preprocessed current signal, and preprocessed rotation speed signal into the empirical mode decomposition algorithm to identify all local maxima and local minima in the signal. Step S2.2: Based on the same signal, the upper envelope is obtained by fitting the local maxima points through cubic spline interpolation, and the lower envelope is obtained by fitting the local minima points. Step S2.3: Calculate the mean envelope of the upper envelope and the lower envelope, and subtract the mean envelope from the original signal of the same signal to obtain the candidate component; Step S2.4: Determine whether the candidate component satisfies the intrinsic mode function condition. If not, use the candidate component as a new input and continue the screening iteration through steps S2.2 to S2.3 until the intrinsic mode function condition is satisfied. Step S2.5: Extract the candidate components that satisfy the intrinsic mode function conditions as intrinsic mode function components, and subtract the intrinsic mode function components from the original signal of the same signal to obtain the remaining signal; Step S2.6: The remaining signal is used as a new input and the decomposition process is repeated from step S2.2 to step S2.5 until the remaining signal changes monotonically and cannot be decomposed further, resulting in several intrinsic mode function components and several residual components.

[0010] Beneficial effects of steps S2.1 to S2.6: Through a progressive execution of steps including local extremum point identification, envelope fitting, mean envelope calculation, intrinsic mode function condition judgment, component extraction, and repeated decomposition, a complete empirical mode decomposition (EMD) technique is formed. This scheme can adaptively decompose complex nonlinear and nonstationary signals into several intrinsic mode function components and residual components with clear physical meanings, achieving multi-scale decomposition and feature separation of the signal. This provides an accurate and reliable hierarchical data foundation for subsequent time-frequency feature extraction, effectively improving the separability and interpretability of signal features. In step S2.1, the preprocessed vibration signal, preprocessed temperature signal, preprocessed current signal, and preprocessed rotation speed signal are input into the empirical mode decomposition algorithm to identify all local maxima and local minima in the signal, providing basic data for envelope fitting. In step S2.2, based on the same signal, the upper envelope is obtained by fitting the local maxima through cubic spline interpolation, and the lower envelope is obtained by fitting the local minima, accurately depicting the upper and lower boundaries of the signal. In step S2.3, the mean envelope of the upper and lower envelopes is calculated, and the mean envelope is subtracted from the original signal of the same signal to obtain candidate components, gradually extracting high-frequency components from the signal. Step S2.4: Determine whether the candidate component satisfies the intrinsic mode function (IMF) condition. If not, use the candidate component as a new input to continue the screening iteration until the IMF condition is met, ensuring that the extracted components have real physical meaning. Step S2.5: Extract the candidate components that satisfy the IMF condition as IMF components, and subtract the IMF components from the original signal of the same signal to obtain the remaining signal, realizing the step-by-step extraction of the signal layer. Step S2.6: Use the remaining signal as a new input to repeat the decomposition process until the remaining signal changes monotonically and cannot be decomposed further, obtaining several IMF components and several residual components, completing the full-scale decomposition of the signal.

[0011] As a further improvement to this application, step S3 involves performing Hilbert-Huang transform calculations on each intrinsic mode function component and each residual component to extract instantaneous frequency, instantaneous amplitude, and marginal spectral energy density, including: Step S3.1: Perform Hilbert transform on each intrinsic mode function component and each residual component, and calculate the real and imaginary parts of the analytic signal through convolution operation; Step S3.2: Calculate the amplitude of the analytic signal using the real part and the imaginary part, and obtain the instantaneous frequency by taking the time derivative of the amplitude of the analytic signal. Step S3.3: Calculate the magnitude of the analytic signal using the real part and the imaginary part to obtain the instantaneous amplitude; Step S3.4: Integrate the instantaneous frequency and the instantaneous amplitude over time to obtain the marginal spectral energy density.

[0012] Beneficial effects of steps S3.1 to S3.4: A complete Hilbert-Huang transform time-frequency feature extraction scheme is constructed through systematic implementation of steps including Hilbert transform, calculation of the real and imaginary parts of the analytic signal, calculation of the amplitude and time derivative of the analytic signal, calculation of the magnitude of the analytic signal, and integration of time-frequency features. This scheme can perform accurate time-frequency analysis on each intrinsic mode function component and each residual component, extracting multidimensional features including instantaneous frequency, instantaneous amplitude, and marginal spectral energy density, comprehensively characterizing the time-varying properties of the signal, fully mining equipment operating status information, and providing a rich and effective feature foundation for subsequent time-series feature learning and health assessment. Specifically, step S3.1 involves performing a Hilbert transform on each intrinsic mode function component and each residual component, and calculating the real and imaginary parts of the analytic signal through convolution operations to construct a complex-form analytic signal representation. Step S3.2 involves calculating the amplitude of the analytic signal using the real and imaginary parts, and obtaining the instantaneous frequency by taking the time derivative of the amplitude, accurately reflecting the instantaneous phase change rate of the signal and capturing the dynamic evolution of the signal's frequency components. Step S3.3 involves calculating the magnitude of the analytic signal using the real and imaginary parts to obtain the instantaneous amplitude, accurately describing the instantaneous intensity change of the signal and characterizing the instantaneous energy level of vibration, temperature, current, and rotational speed signals. Step S3.4 involves integrating the instantaneous frequency and instantaneous amplitude over time to obtain the marginal spectral energy density, comprehensively characterizing the signal's distribution characteristics across different frequency components from an energy perspective, and providing frequency domain energy features for equipment status identification.

[0013] As a further improvement to this application, in step S4, all instantaneous frequencies, all instantaneous amplitudes, and all marginal spectral energy densities are input into the Long Short-Term Memory network for temporal feature learning. The cell state and hidden state are updated through a gating mechanism involving a forget gate, an input gate, and an output gate, resulting in a temporal feature vector, including: Step S4.1: Concatenate all instantaneous frequencies, all instantaneous amplitudes, and all marginal spectral energy densities in chronological order to form a feature sequence; Step S4.2: Input the feature sequence into the Long Short-Term Memory (LSTM) network, and calculate the forgetting gate coefficient through the forgetting gate of the LSM network to control the degree of retention of the cell state at the previous time step; Step S4.3: Calculate the input gating coefficient and candidate cell state through the input gate of the Long Short-Term Memory network to determine the new information at the current moment; Step S4.4: Update the current cell state using the forgetting gating coefficient and the input gating coefficient; Step S4.5: Calculate the output gating coefficient through the output gate of the Long Short-Term Memory Network to determine the hidden state at the current moment based on the updated cell state; Step S4.6: Repeat steps S4.2 to S4.5 until all time steps have been processed to obtain a time-series feature vector with temporal dependencies.

[0014] Beneficial effects of steps S4.1 to S4.6: A complete temporal feature learning scheme is formed through cyclical execution of steps such as feature sequence concatenation, temporal learning via Long Short-Term Memory (LSTM) networks, calculation of forgetting gate input-output gate gating mechanisms, and updating of cell states and hidden states. This scheme can automatically capture temporal dependencies and long-term dependency patterns in feature sequences, extracting deep features with temporal contextual information from instantaneous frequency, instantaneous amplitude, and marginal spectral energy density sequences. This effectively improves the ability to identify equipment performance degradation trends and early fault symptoms, providing feature vectors containing dynamic temporal information for subsequent health assessments. In this process, step S4.1 involves concatenating all instantaneous frequencies, instantaneous amplitudes, and marginal spectral energy densities in chronological order to form a feature sequence, constructing a unified high-dimensional temporal input and providing a standard data format for Long Short-Term Memory (LSTM) network processing. Step S4.2 inputs the feature sequence into the LSM network, calculating the forgetting gate coefficient to control the retention of the cell state from the previous time step, thus achieving intelligent filtering of historical information and maintenance of long-term memory. Step S4.3 calculates the input gate coefficient and candidate cell states through the input gate to determine the newly added information at the current time step, realizing the contribution of the current input to updating the network state. Step S4.4 updates the current cell state using the forgetting gate coefficient and the input gate coefficient, balancing the influence ratio of historical information and current information to ensure dynamic optimization of the network state. Step S4.5 calculates the output gate coefficient through the output gate, determines the hidden state at the current time step based on the updated cell state, and outputs feature information useful for subsequent tasks. Step S4.6 repeats steps S4.2 to S4.5 until all time steps are processed, obtaining a temporal feature vector with temporal dependencies, and fully extracting long-term dependency information from the feature sequence.

[0015] As a further improvement to this application, step S5 involves inputting the temporal feature vector into an attention mechanism to calculate the weight coefficient for each time step, and then calculating the weighted aggregated temporal feature vector by weighted summation and aggregating the feature information of key time steps, including: Step S5.1: Input the temporal feature vector into the attention mechanism and calculate the relevance score between the query vector and the key vector through dot product operation; Step S5.2: Perform an exponential operation on the correlation score and divide it by the sum of the exponential operation results of all time steps to obtain the weight coefficient of each time step. Step S5.3: The feature information of each time step in the time series feature vector is weighted and summed using the weight coefficients to obtain the weighted aggregated time series feature vector.

[0016] Beneficial effects of steps S5.1 to S5.3: By leveraging the synergistic effects of attention mechanisms in calculating the correlation between query vectors and key vectors, performing correlation score exponent calculations and normalization, and weighted summation of weight coefficients, a complete temporal feature aggregation scheme is constructed. This scheme can automatically identify the importance of each time step in the temporal feature vector, highlighting key moment features that significantly contribute to equipment health status assessment, suppressing interference from irrelevant or redundant features, improving the focus and effectiveness of feature representation, and providing a compact feature representation containing key information for subsequent equipment health scoring. In step S5.1, the temporal feature vector is input into the attention mechanism. The relevance score between the query vector and the key vector is calculated through dot product operation to quantify the degree of correlation between features at different time steps and evaluate the relative importance of features at each time step to the current task. Step S5.2, the relevance score is exponentially calculated and divided by the sum of the exponential calculation results of all time steps to obtain the weight coefficient of each time step. The relevance score is converted into a probability distribution with values ​​ranging from zero to one to ensure that the sum of all weight coefficients is one, thereby achieving reasonable allocation and normalization of weights. Step S5.3, the feature information of each time step in the temporal feature vector is weighted and summed through the weight coefficient. Features at different times are aggregated according to their importance to form a weighted aggregated temporal feature vector that enhances the features of key time steps and suppresses the features of secondary time steps, thereby improving the discriminative ability and robustness of feature representation.

[0017] As a further improvement to this application, step S6 involves inputting the weighted aggregated temporal feature vector into a support vector data description algorithm to construct a hypersphere model of the device health status. The distance from the sample to the center of the hypersphere is calculated by mapping to a high-dimensional feature space using a kernel function, resulting in a device health score, including: Step S6.1: Map the weighted aggregated temporal feature vectors to a high-dimensional feature space using the support vector data description algorithm; Step S6.2: Calculate the Euclidean distance matrix between sample vectors in the high-dimensional feature space using the Gaussian radial basis kernel function; Step S6.3: Solve the quadratic programming optimization problem of the Euclidean distance matrix to determine the center and radius of the hypersphere in order to minimize the hypersphere volume containing all weighted aggregated temporal feature vectors; Step S6.4: Calculate the distance from each weighted aggregated temporal feature vector to the center of the hypersphere to obtain the sample distance value; Step S6.5: Normalize all sample distance values ​​according to the radius of the hypersphere to obtain a device health score set based on multiple dimensions.

[0018] Beneficial effects of steps S6.1 to S6.5: A complete equipment health assessment scheme is constructed through the systematic implementation of steps including high-dimensional feature space mapping using the support vector data description algorithm, calculation of the Euclidean distance matrix using the Gaussian radial basis function kernel, solving the quadratic programming optimization problem, and calculation and normalization of sample distance values. This scheme can construct a hypersphere model of equipment health status based on weighted aggregated temporal feature vectors. By calculating the distance from the sample to the center of the hypersphere, it achieves quantitative assessment and anomaly detection of equipment health status, accurately reflecting the current operating status and health level of the equipment, and providing a scientific basis for equipment maintenance decisions. Specifically, step S6.1 uses the support vector data description algorithm to map the weighted aggregated temporal feature vectors to a high-dimensional feature space, mapping low-dimensional linearly inseparable features to a high-dimensional linearly separable space, providing suitable feature representations for subsequent distance calculation and anomaly detection; step S6.2 calculates the Euclidean distance matrix between sample vectors in the high-dimensional feature space using the Gaussian radial basis function kernel, accurately measuring the similarity and difference between different samples in the feature space, providing basic distance information for the construction of the hypersphere model; step S6.3 solves the quadratic programming optimization problem of the Euclidean distance matrix to determine the hypersphere. The center and radius of the sphere are set to minimize the volume of the hypersphere and include all weighted aggregated temporal feature vectors to establish a decision boundary that can accurately characterize the distribution range of normal equipment status. Step S6.4 calculates the distance from each weighted aggregated temporal feature vector to the center of the hypersphere to obtain the sample distance value, quantifying the degree to which each sample deviates from the center of the normal status distribution. Step S6.5 normalizes all sample distance values ​​according to the radius of the hypersphere to obtain a set of equipment health scores based on multiple dimensions, unifying the numerical range of scores in different dimensions, and facilitating intuitive assessment and horizontal comparison of equipment health status.

[0019] To achieve the above objectives, this application also provides the following technical solutions: A health detection system for packaging machine equipment, wherein the health detection system is applied to the health detection method described above, and the health detection system includes: The raw signal processing module is used to collect vibration signals, temperature signals, current signals, and speed signals of the packaging machine in real time through external sensors, and to perform noise reduction and normalization processing on the collected raw signals. The original signal decomposition module is used to input the preprocessed vibration signal, preprocessed temperature signal, preprocessed current signal, and preprocessed rotation speed signal into the empirical mode decomposition algorithm to decompose them into several intrinsic mode function components and several residual components. The signal component transformation module is used to perform Hilbert-Huang transform on each intrinsic mode function component and each residual component to calculate the time-frequency characteristics, so as to extract the instantaneous frequency, instantaneous amplitude, and marginal spectral energy density. The signal component feature learning module is used to input all instantaneous frequencies, all instantaneous amplitudes, and all marginal spectral energy densities into the long short-term memory network for temporal feature learning. The cell state and hidden state are updated through the gating mechanism of forget gate, input gate, and output gate to obtain the temporal feature vector. The temporal feature vector weighted aggregation module is used to input the temporal feature vector into the attention mechanism to calculate the weight coefficient of each time step, and to calculate the weighted aggregated temporal feature vector by weighted summation and aggregation of the feature information of key time steps; The device health score calculation module is used to input the weighted aggregated time-series feature vector into the support vector data description algorithm to construct a hypersphere model of the device health status, and calculate the distance from the sample to the center of the hypersphere by mapping to a high-dimensional feature space through a kernel function to obtain the device health score.

[0020] To achieve the above objectives, this application also provides the following technical solutions: An electronic device includes a processor and a memory coupled to the processor, the memory storing program instructions executable by the processor; when the processor executes the program instructions stored in the memory, it implements the health detection method for packaging machine equipment as described above.

[0021] To achieve the above objectives, this application also provides the following technical solutions: A computer-readable storage medium storing program instructions that, when executed by a processor, can implement the health detection method for packaging machine equipment as described above. Attached Figure Description

[0022] Figure 1 This is a schematic flowchart illustrating the steps of an embodiment of a health detection method for packaging machine equipment according to this application; Figure 2 This is a functional module diagram of an embodiment of a health detection system for packaging machine equipment according to this application; Figure 3 This is a schematic diagram of the structure of an embodiment of the electronic device of this application; Figure 4 This is a schematic diagram of the structure of one embodiment of the storage medium of this application. Detailed Implementation

[0023] The technical solutions of the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only a part of the embodiments of this application, and not all of the embodiments. Based on the embodiments of this application, all other embodiments obtained by those of ordinary skill in the art without creative effort are within the scope of protection of this application.

[0024] The terms "first," "second," and "third" in this application are used for descriptive purposes only and should not be construed as indicating or implying relative importance or implicitly specifying the number of technical features indicated. Therefore, a feature defined as "first," "second," or "third" may explicitly or implicitly include at least one of that feature. In the description of this application, "multiple" means at least two, such as two, three, etc., unless otherwise explicitly specified. All directional indications (such as up, down, left, right, front, back, etc.) in the embodiments of this application are only used to explain the relative positional relationships and movements between components in a specific orientation (e.g., as shown in the figures). If the specific orientation changes, the directional indications also change accordingly. Furthermore, the terms "comprising" and "having," and any variations thereof, are intended to cover non-exclusive inclusion. For example, a process, method, system, product, or device that includes a series of steps or units is not limited to the listed steps or units, but may optionally include steps or units not listed, or may optionally include other steps or units inherent to these processes, methods, products, or devices.

[0025] In this document, the term "embodiment" means that a particular feature, structure, or characteristic described in connection with an embodiment may be included in at least one embodiment of this application. The appearance of this phrase in various places throughout the specification does not necessarily refer to the same embodiment, nor is it a mutually exclusive, independent, or alternative embodiment. It will be explicitly and implicitly understood by those skilled in the art that the embodiments described herein can be combined with other embodiments.

[0026] It should be noted that, due to the limited types and number of symbols or letters that can represent specific meanings, for embodiments with many formulas or codes, there may be situations where symbols or letters cannot meet the usage requirements. Therefore, the interpretation of formula symbols in the steps or sub-steps of the embodiments is only valid for the current step or sub-step.

[0027] If the same symbol has different interpretations in different steps or sub-steps, the interpretation in the current step or sub-step shall prevail; if the same symbol appears in different steps or sub-steps, but no interpretation is given in subsequent steps or sub-steps after its first appearance, the interpretation in the first step or sub-step shall be used.

[0028] like Figure 1 As shown, this embodiment provides an example of a health detection method for packaging machine equipment. In this embodiment, the health detection method includes the following steps: Step S1: The vibration signal, temperature signal, current signal and speed signal of the packaging machine are collected in real time by external sensors, and the collected raw signals are denoised and normalized.

[0029] Preferably, the core objective of this step is to obtain high-quality time-series signals of the core operating status of the packaging machine. Through synchronous acquisition by multiple physical quantity sensors and multi-stage signal cleaning and standardization processing, noise, dimensional differences, timing misalignment, and data loss in the original signal are eliminated, providing a stable and reliable input data source for subsequent time-frequency feature extraction and health assessment.

[0030] Furthermore, step S1 specifically includes the following steps: Step S1.1: Acquire vibration acceleration signal through accelerometer, acquire temperature signal through thermocouple sensor, acquire current signal through Hall sensor, and acquire rotation speed signal through encoder to obtain original vibration signal, original temperature signal, original current signal, and original rotation speed signal.

[0031] Preferably, this step involves hardware-level synchronous acquisition of multiple physical quantity signals for the core moving parts, electrical system, thermal system, and transmission system of the GDX500 soft-pack cigarette packaging machine.

[0032] For vibration signal acquisition, a piezoelectric accelerometer can be used. The sensor is installed at the bearing seats of the second, third, and fourth wheels of the packaging machine's main unit (X) and the bearing seat of the main motor transmission box of the main unit (CH). The acquisition dimensions are radial and axial vibration acceleration. The sampling frequency is calibrated to 20kHz, matching the meshing frequency corresponding to the highest speed of the packaging machine's main shaft and the Nyquist sampling requirement of bearing fault characteristic frequency more than 10 times. The range is calibrated to ±50g, and the frequency response range is 0.5Hz~10kHz. The output is the original vibration acceleration time-domain signal, i.e., the original vibration signal.

[0033] For temperature signal acquisition, a K-type thermocouple sensor can be used. The sensor is installed in the heater modules of the packaging machine, the lubrication circuit of the transmission box, and the stator winding of the main motor. The sampling frequency is calibrated to 10Hz to match the slow change characteristics of the temperature signal. The range is calibrated to 0℃ to 300℃, the measurement accuracy is ±0.2℃, and the original temperature signal is output.

[0034] For current signal acquisition, a closed-loop Hall current sensor can be used, connected in series in the power supply circuits of each main motor, servo motor, heater, and electromagnet of the packaging machine. The sampling frequency is calibrated to 1kHz, covering the dynamic response frequency band of motor start-stop and load fluctuation. The range is calibrated to 0A to 20A (motor circuit) and 0A to 5A (heater / electromagnet circuit) according to the load, with a measurement accuracy of 0.2 class, and outputting the raw current signal.

[0035] For speed signal acquisition, an incremental photoelectric encoder can be used, installed at the end of the X and CH main spindles of the packaging machine. The sampling frequency is calibrated to 100Hz, the encoder resolution is 1024 lines, and the output is the real-time speed signal of the spindle, i.e., the original speed signal.

[0036] In summary, all sensors are synchronized at the hardware level through a synchronous acquisition card. The trigger source is uniformly the zero-position pulse of the packaging machine spindle encoder. All raw signals are bound to millisecond-level hardware timestamps, and the timestamp reference is uniformly the system clock of the acquisition device, ensuring that the timing reference of multi-source signals is consistent. The final output is raw vibration signal, raw temperature signal, raw current signal, and raw speed signal with uniform timestamps.

[0037] Step S1.2: Input the original vibration signal, original temperature signal, original current signal and original rotation speed signal into the wavelet threshold denoising algorithm respectively, and decompose and reconstruct them through the sym4 wavelet basis function to obtain the denoised vibration signal, denoised temperature signal, denoised current signal and denoised rotation speed signal.

[0038] Preferably, the core of this step is to eliminate industrial electromagnetic interference, mechanical background noise, and white noise from the acquisition circuit in the original signal, retain the effective feature components in the signal that reflect the operating status of the equipment, and use a wavelet threshold denoising algorithm.

[0039] For the selection of wavelet basis functions, the sym4Symlets wavelet basis function can be used. This wavelet basis has tight support, linear phase characteristics and good time-domain localization ability, which matches the feature extraction requirements of non-stationary transient signals such as packaging machine vibration and current, and avoids signal distortion.

[0040] Specifically, for wavelet decomposition level calibration, five-level wavelet decomposition can be performed on the original vibration signal, original temperature signal, original current signal, and original rotational speed signal, respectively. This decomposes the original signal into one level of high-frequency detail components (D1), two levels of high-frequency detail components (D2), three levels of high-frequency detail components (D3), four levels of high-frequency detail components (D4), five levels of high-frequency detail components (D5), and five levels of low-frequency approximation components (A5). The high-frequency components correspond to noise and transient impact characteristics in the signal, while the low-frequency components correspond to the steady-state trend term in the signal.

[0041] For the threshold function and threshold calculation, a soft threshold function can be used to threshold the high-frequency detail components of each layer. The mathematical expression of the soft threshold function is: In the formula, w j,k Let k be the wavelet coefficient of the j-th level decomposition. Let λ represent the wavelet coefficients after thresholding, and λ be the threshold. The threshold λ is calculated using a general threshold (Donoho threshold), and the formula is as follows: In the formula, σ is the noise standard deviation, estimated using the median absolute deviation method, σ = median(|w j,k |) / 0.6745, where N is the total length of the wavelet coefficients corresponding to the high-frequency components.

[0042] For signal reconstruction, the high-frequency detail components of each layer after threshold processing can be subjected to inverse wavelet transform with the low-frequency approximation component A5 without threshold processing to complete signal reconstruction, and the denoised vibration signal, denoised temperature signal, denoised current signal, and denoised rotation speed signal are output respectively.

[0043] For example, this step can be implemented using the following pseudocode: import pywt import numpy as np def wavelet_denoising(original_signal, wavelet_base='sym4',decompose_level=5): # 1. Perform multi-level wavelet decomposition to separate low-frequency and high-frequency components. wavelet_coeffs = pywt.wavedec(original_signal, wavelet=wavelet_base, level=decompose_level) approx_coeff = wavelet_coeffs[0] detail_coeffs = wavelet_coeffs[1:] # 2. Perform soft thresholding on high-frequency detail components of each layer. denoised_detail_coeffs = [] for single_layer_coeff in detail_coeffs: # Estimate the standard deviation of noise in the current layer noise_sigma = np.median(np.abs(single_layer_coeff)) / 0.6745 # Calculate the common threshold for the current layer threshold = noise_sigma * np.sqrt(2 * np.log(len(single_layer_coeff))) # Soft threshold calculation processed_coeff = np.sign(single_layer_coeff) * np.maximum(np.abs(single_layer_coeff) - threshold, 0) denoised_detail_coeffs.append(processed_coeff) # 3. Reconstructing the denoised signal using inverse wavelet transform denoised_signal = pywt.waverec([approx_coeff] + denoised_detail_coeffs, wavelet=wavelet_base) return denoised_signal Step S1.3: The denoised vibration signal, denoised temperature signal, denoised current signal, and denoised speed signal are normalized by Z-score standardization to obtain normalized vibration signal, normalized temperature signal, normalized current signal, and normalized speed signal.

[0044] Preferably, the core of this step is to eliminate the differences in dimensions and numerical magnitudes between different physical quantity signals, so as to avoid feature learning bias caused by the imbalance of input data magnitudes in subsequent algorithms.

[0045] For the calibration of the standardized benchmark dataset, 72 hours of steady-state operating data of the GDX500 packaging machine under rated operating conditions and fault-free operation can be collected in advance to construct a standardized benchmark dataset. The global mean and global standard deviation of four types of signals, namely vibration, temperature, current and speed, are calculated in the benchmark dataset as fixed standardized parameters to avoid standardization benchmark drift caused by real-time calculation.

[0046] For Z-score standardization calculation, point-by-point standardization calculation can be performed on each data point in the denoised vibration signal, denoised temperature signal, denoised current signal, and denoised speed signal. Since Z-score standardization is a mature existing technology, the calculation process of Z-score standardization will not be described in detail.

[0047] Preferably, after performing full-sequence standardization calculations on each signal, normalized vibration signal, normalized temperature signal, normalized current signal, and normalized rotational speed signal are output respectively. The mean of the standardized signal sequence is 0, and the standard deviation is 1, completely eliminating the difference in dimensions and numerical magnitude.

[0048] Step S1.4: Align the normalized vibration signal, normalized temperature signal, normalized current signal, and normalized rotational speed signal according to the timestamp, and fill in the missing data points using a linear interpolation method.

[0049] Preferably, the core of this step is to solve the timing misalignment problem caused by different sampling frequencies and acquisition delays of multi-source signals, as well as the missing data problem caused by acquisition interference and signal interruption, so as to ensure strict synchronization of multi-source signals in the time dimension.

[0050] For a unified time reference and alignment resolution, the timestamps of all normalized signals can be uniformly converted into millisecond-level UTC timestamps. The time alignment resolution is set to 1ms, and a time axis reference sequence with equal time intervals is constructed. The start and end times of the time axis are completely matched with the acquisition time range of all signals.

[0051] For timestamp alignment, for each normalized signal, the signal data points can be mapped to a unified 1ms interval time axis based on its original timestamp; for signals with multiple data points in the same time interval, the mean method is used to merge them into a single data point; for time intervals without corresponding data points, they are marked as missing data points to be filled.

[0052] For the determination of missing data, the threshold for determining missing data can be set as follows: if there are no valid data points in two or more consecutive alignment time intervals, it is marked as a continuous missing segment; if there are no data points in a single interval, it is marked as a single missing point.

[0053] For linear interpolation imputation, a one-dimensional linear interpolation method is used to imput all marked missing data points. The formula for calculating linear interpolation is: x interp =x0+(x1-x0)×(t-t0) / (t1-t0). Where t0 and t1 are the timestamps corresponding to the nearest valid data points before and after the missing segment, x0 and x1 are the signal values ​​at the corresponding timestamps, t is the timestamp of the missing point to be filled, and x... interp The filler value is calculated using interpolation.

[0054] For boundary constraint processing, missing segments at the beginning and end of the signal can be filled by forward filling (missing segment at the beginning) and backward filling (missing segment at the end). The filling value is the value of the nearest valid steady-state data point to avoid interpolation distortion.

[0055] Preferably, after completing the alignment and filling of the entire time axis, the output includes aligned vibration signal, aligned temperature signal, aligned current signal, and aligned rotation speed signal with completely consistent time axis and no missing data points.

[0056] Step S1.5: After filling in the missing data points, the signal is filtered by a moving average algorithm to eliminate high-frequency jitter interference, resulting in preprocessed vibration signal, preprocessed temperature signal, preprocessed current signal, and preprocessed rotation speed signal.

[0057] Preferably, the core of this step is to further eliminate the high-frequency random jitter remaining in the signal after alignment and filling, smooth the signal curve, and at the same time retain the trend terms and effective transient characteristics in the signal that reflect changes in the device state.

[0058] Specifically, for the sliding window length calibration, an appropriate sliding window length can be set for the changing characteristics of different types of signals, achieving an optimal balance between noise suppression and feature preservation. ① Aligned vibration signal: A 5-point sliding window is used to match the high-frequency transient characteristics of the vibration signal, avoiding excessive smoothing that would cause the loss of fault impact characteristics.

[0059] ② Aligned temperature signal: A 10-point sliding window is used to match the slow-changing characteristics of the temperature signal and minimize random measurement jitter.

[0060] ③ Aligned current signal: A 7-point sliding window is used to balance the load dynamic response and noise suppression requirements of the current signal.

[0061] ④ Aligned speed signal: A 5-point sliding window is used to retain the transient fluctuation characteristics of the speed and eliminate encoder counting jitter.

[0062] ⑤ Moving average filtering calculation: For each data point in the signal sequence, take all data points within the corresponding window length with that point as the center, and calculate the arithmetic mean as the filtered data point.

[0063] Preferably, for data points at the beginning and end of the signal that cannot meet the length of the center window, the beginning and end of the signal are expanded by mirror filling method, with a filling length of (N-1) / 2, and then the moving average calculation is performed to avoid amplitude distortion at the beginning and end of the signal.

[0064] Preferably, after completing the full-sequence moving average filtering, the preprocessed vibration signal, preprocessed temperature signal, preprocessed current signal, and preprocessed rotation speed signal are finally output as the sole input data source for the subsequent step S2.

[0065] Beneficial effects of steps S1.1 to S1.5: By organically combining multi-sensor collaborative acquisition, wavelet threshold denoising, Z-score normalization, timestamp alignment and linear interpolation imputation, and moving average filtering, a complete signal preprocessing technology solution is formed. This solution can effectively eliminate noise interference, dimensional differences, time deviations, and data gaps in the original signal, significantly improving signal quality and consistency. It provides high-quality basic data for subsequent signal decomposition, feature extraction, and health assessment, solving problems such as residual noise, data mismatch, and poor feature extractability in traditional preprocessing methods when processing multi-source heterogeneous signals. The process involves several steps: Step S1.1 involves acquiring vibration acceleration, temperature, current, and rotational speed signals using an accelerometer, thermocouple, Hall effect sensor, and encoder, respectively, to comprehensively obtain key physical quantities reflecting the equipment's status during operation. Step S1.2 involves inputting the acquired raw signals into a wavelet threshold denoising algorithm, decomposing and reconstructing them using the sym4 wavelet basis function to effectively remove high-frequency noise and random interference, thus improving the signal-to-noise ratio. Step S1.3 involves normalizing the denoised signals using Z-score standardization to eliminate dimensional differences and inconsistent numerical ranges between signals from different sensors, ensuring that all signals are within a uniform numerical scale. Step S1.4 involves aligning the normalized signals according to timestamps and filling in missing data points using linear interpolation to ensure the synchronization and integrity of multiple signals in the time dimension. Step S1.5 involves applying a moving average filtering algorithm to the signals after filling in missing data points to eliminate high-frequency jitter interference and smooth instantaneous fluctuations in the signals, resulting in preprocessed vibration, temperature, current, and rotational speed signals.

[0066] Step S2: Input the preprocessed vibration signal, preprocessed temperature signal, preprocessed current signal, and preprocessed rotation speed signal into the empirical mode decomposition algorithm to decompose them into several intrinsic mode function components and several residual components.

[0067] Preferably, the core objective of this step is to perform adaptive multi-scale decomposition on the preprocessed vibration signal, preprocessed temperature signal, preprocessed current signal, and preprocessed rotation speed signal output from step S1. This decomposes the non-stationary time-series signal containing multi-band composite features into a series of intrinsic mode function (IMF) components with single-component characteristics and residual components that characterize the overall trend of the signal. This eliminates the interference of non-stationary characteristics on subsequent time-frequency feature extraction and provides an input data source that meets the analysis requirements for the Hilbert-Huang transform.

[0068] Furthermore, step S2 specifically includes the following steps: Step S2.1: Input the preprocessed vibration signal, preprocessed temperature signal, preprocessed current signal, and preprocessed rotation speed signal into the empirical mode decomposition algorithm to identify all local maxima and local minima in the signals.

[0069] Preferably, in this step, extreme point identification is performed independently for the preprocessed vibration signal, preprocessed temperature signal, preprocessed current signal, and preprocessed rotation speed signal.

[0070] The details of the entire process for identifying extreme points of a single signal are as follows: (1) Quantization judgment rule for discrete signal extreme points: For the discrete time sequence x(n) to be decomposed, n=0,1,2,...,N-1, where N is the total length of the signal sequence and n is the discrete time index, the full sequence traversal judgment is strictly performed according to the following rules: ① Local maximum point: if and only if x(n)>x(n-1) and x(n)>x(n+1), where the range of n is 1≤n≤N-2; for the first and last endpoints of the sequence n=0 and n=N-1, if x(0)>x(1) or x(N-1)>x(N-2), then it is determined to be a boundary maximum point.

[0071] ② Local minimum point: if and only if x(n)<x(n-1) and x(n)<x(n+1), where the range of n is 1≤n≤N-2; for the first and last endpoints of the sequence n=0 and n=N-1, if x(0)<x(1) or x(N-1)<x(N-2), then it is determined to be a boundary minimum point.

[0072] (2) Generation of extreme point sequence: After completing the full sequence traversal, the local maximum point sequence and the local minimum point sequence are generated in strict order of time. Each set of sequences stores the time index and signal amplitude corresponding to the extreme point synchronously.

[0073] Preferably, in response to the common peak and valley phenomena in industrial field signals, an extreme point deduplication rule is set: for three or more consecutive extreme points with the same amplitude, only the first and last extreme points are retained to avoid singular solutions in subsequent envelope fitting.

[0074] Preferably, for each preprocessed signal, the corresponding sequence of local maxima and local minima are output as the sole input for subsequent envelope fitting.

[0075] Step S2.2: Based on the same signal, the upper envelope is obtained by fitting local maxima points and the lower envelope is obtained by fitting local minima points through cubic spline interpolation.

[0076] Preferably, this step uses the local maximum and local minimum point sequences of the same signal to respectively fit and generate upper and lower envelopes covering the entire time series. The interpolation algorithm adopts the standard cubic spline interpolation algorithm, which can ensure that the fitted curve is second-order continuous and differentiable at all interpolation nodes, has optimal envelope smoothness, no Runge phenomenon, and fully adapts to the accuracy requirements of the envelope in the empirical mode decomposition algorithm, avoiding component screening deviation caused by envelope distortion.

[0077] Preferably, cubic spline interpolation, for a given sequence of interpolation nodes (t) at extreme points, i ,y i ), i=0,1,...,M-1, where t i The time index of the extreme point, y i For the corresponding amplitude, M is the number of extreme points, and the cubic spline interpolation function S(t) must simultaneously satisfy the following constraints: (1) Interpolation constraint: S(t) i )=y i That is, the fitted curve strictly passes through all extreme points.

[0078] (2) Smoothness constraint: In each interpolation node interval [t i ,t i+1 Within the interval, S(t) is a cubic polynomial, and its first derivative S'(t) and second derivative S''(t) are continuous throughout the entire interval.

[0079] (3) Boundary constraints: Natural boundary conditions are adopted, and the second derivative at the first and last interpolation nodes is set to 0, that is, S''(t0)=0, S''(t0)=0, and S''(t0)=0. M-1 =0, suppressing the endpoint effect at the beginning and end of the envelope.

[0080] Preferably, the upper envelope fitting uses the sequence of local maxima as interpolation nodes, and calculates an upper envelope sequence e of the same length as the original signal using a cubic spline interpolation algorithm. max (n), where each time index n corresponds to a unique upper envelope amplitude; the lower envelope fitting uses the sequence of local minima as interpolation nodes, and calculates the lower envelope sequence e of the same length as the original signal using a cubic spline interpolation algorithm. min (n), where each time index n corresponds to a unique lower envelope amplitude.

[0081] Step S2.3: Calculate the mean envelope of the upper and lower envelopes, and subtract the mean envelope from the original signal of the same signal to obtain the candidate components.

[0082] Preferably, this step calculates the mean envelope based on the upper and lower envelope sequences generated in the previous step, and extracts the candidate components for the screening iteration.

[0083] Specifically, for the point-by-point calculation of the mean envelope, the arithmetic mean of the upper and lower envelopes can be calculated for each time point n within the full-time index of the original signal, resulting in the mean envelope sequence m(n), which is exactly the same length as the original signal.

[0084] For candidate component extraction, the mean envelope sequence can be subtracted point by point from the input signal sequence of the current screening iteration to obtain the candidate component sequence h(n). The calculation formula is: h(n) = x input (n)-m(n), n=0,1,2,...,N-1 where, x input (n) is the input signal for the current sieving iteration. In the first iteration, it is the preprocessed original signal of the corresponding path, and in subsequent iterations, it is the candidate component sequence output from the previous iteration.

[0085] Step S2.4: Determine whether the candidate component satisfies the intrinsic mode function condition. If not, use the candidate component as a new input and continue the screening iteration through steps S2.2 to S2.3 until the intrinsic mode function condition is satisfied.

[0086] Preferably, the core of this step is to perform compliance quantification on the candidate component sequences, clarify the feasible judgment criteria for the intrinsic mode function, set iterative convergence rules and iteration upper limits, and ensure that the components obtained by decomposition meet the input requirements of the Hilbert-Huang transform.

[0087] Specifically, the quantization determination of the Intrinsic Mode Function (IMF) must simultaneously meet the following two core conditions: (1) Matching condition for the number of extreme points and zero crossings: The total number of local maxima and local minima in the entire signal sequence is no more than 1 different from the number of zero crossings of the signal.

[0088] (2) Local symmetry condition: At any point in time, the mean of the upper and lower envelopes of the signal approaches 0. In engineering, this condition is quantified by screening the stop threshold (SD) to avoid the impossibility of an absolute 0 value.

[0089] For the screening stop threshold (SD), the industry-standard convergence criterion SD can be used, and the calculation formula is as follows: Among them, h k-1 (n) represents the candidate component sequence from the previous iteration, h k (n) represents the candidate component sequence in the current iteration, and N is the signal sequence length. In this step, the SD threshold is calibrated to 0.25. This threshold can control the number of iterations while ensuring the physical meaning of the components, and avoid component distortion caused by over-screening.

[0090] Preferably, the maximum number of iterations can be set to 100 to prevent the algorithm from getting stuck in an infinite iteration loop due to signal anomalies.

[0091] Preferably, the iteration control rule is as follows: (1) If a candidate component satisfies both of the above IMF quantization conditions at the same time, or the SD value is less than or equal to 0.25, then the candidate component is determined to meet the intrinsic mode function requirements, and the screening iteration of the current round is stopped.

[0092] (2) If the candidate component does not meet the above conditions, the candidate component is used as a new input signal, and the process returns to steps S2.2 to S2.3 to perform a new round of screening iteration, while accumulating the number of iterations.

[0093] (3) If the convergence condition is not met even after the number of iterations reaches the upper limit of 100, the iteration is forcibly stopped and the current candidate component is output as an intrinsic mode function component.

[0094] Step S2.5: Extract the candidate components that satisfy the intrinsic mode function conditions as intrinsic mode function components, and subtract the intrinsic mode function components from the original signal of the same signal to obtain the remaining signal.

[0095] Preferably, this step is based on the compliant candidate components after the screening iteration convergence, to complete the formal extraction of the intrinsic mode function components, and to calculate the remaining signal after decomposition, so as to provide input for subsequent multi-round decomposition.

[0096] Specifically, for the extraction of intrinsic mode function components, the candidate component sequences that meet the IMF conditions as determined in step S2.4 are marked as the k-th order intrinsic mode function component (IMF). k(n), where k is the order of the currently extracted IMF component. When it is extracted for the first time, k=1. After each successful extraction of a component, k increases by 1. The length of the component sequence is completely consistent with the original preprocessed signal.

[0097] For the calculation of the residual signal, the input signal of the current decomposition round can be used to subtract the extracted intrinsic mode function components point by point to obtain the residual signal sequence r. k (n), the calculation formula is: r k (n)=x input (n)-imf k (n), n=0,1,2,...,N-1. Where, x input (n) is the input signal for the current decomposition round. In the first decomposition, it is the preprocessed original signal of the corresponding path. In subsequent decomposition rounds, it is the remaining signal sequence output from the previous decomposition round.

[0098] In step S2.6, the remaining signal is used as a new input and the decomposition process is repeated from step S2.2 to step S2.5 until the remaining signal changes monotonically and cannot be decomposed further, resulting in several intrinsic mode function components and several residual components.

[0099] Preferably, the core of this step is to determine the feasibility of decomposing the remaining signal, clarify the quantization conditions for terminating the decomposition, complete the full-sequence EMD decomposition of a single signal, and finally output all effective components.

[0100] Specifically, the EMD decomposition process for the current signal is terminated if any one of the following conditions is met: (1) The total number of local maxima and local minima of the remaining signal sequence is less than 2, so it is impossible to generate effective upper and lower envelopes by cubic spline interpolation fitting, and it is impossible to continue the sieving iteration.

[0101] (2) The remaining signal sequence is a monotonic function or a constant function. There is only one extreme point in the whole sequence. There are no fluctuation characteristics, and it is impossible to extract new intrinsic mode function components.

[0102] (3) The peak amplitude of the remaining signal is less than 0.01% of the maximum amplitude of the original signal, and has no practical physical meaning. Further decomposition has no engineering value.

[0103] Preferably, when the decomposition termination condition is met, the current remaining signal sequence is marked as the residual component res(n). This component characterizes the overall steady-state change trend of the original signal and has no high-frequency fluctuation characteristics. For a single preprocessed signal, after completing all decompositions, several intrinsic mode function components imf1(n), imf2(n), ..., imf are output in descending order of frequency. K(n) (K is the total number of IMF components obtained from the decomposition of the signal), and 1 residual component res(n). After performing independent EMD decomposition on the preprocessed vibration signal, preprocessed temperature signal, preprocessed current signal, and preprocessed rotation speed signal, the intrinsic mode function components and residual components of all channels are integrated to form a complete component set.

[0104] Beneficial effects of steps S2.1 to S2.6: Through a progressive execution of steps including local extremum point identification, envelope fitting, mean envelope calculation, intrinsic mode function condition judgment, component extraction, and repeated decomposition, a complete empirical mode decomposition (EMD) technique is formed. This scheme can adaptively decompose complex nonlinear and nonstationary signals into several intrinsic mode function components and residual components with clear physical meanings, achieving multi-scale decomposition and feature separation of the signal. This provides an accurate and reliable hierarchical data foundation for subsequent time-frequency feature extraction, effectively improving the separability and interpretability of signal features. In step S2.1, the preprocessed vibration signal, preprocessed temperature signal, preprocessed current signal, and preprocessed rotation speed signal are input into the empirical mode decomposition algorithm to identify all local maxima and local minima in the signal, providing basic data for envelope fitting. In step S2.2, based on the same signal, the upper envelope is obtained by fitting the local maxima through cubic spline interpolation, and the lower envelope is obtained by fitting the local minima, accurately depicting the upper and lower boundaries of the signal. In step S2.3, the mean envelope of the upper and lower envelopes is calculated, and the mean envelope is subtracted from the original signal of the same signal to obtain candidate components, gradually extracting high-frequency components from the signal. Step S2.4: Determine whether the candidate component satisfies the intrinsic mode function (IMF) condition. If not, use the candidate component as a new input to continue the screening iteration until the IMF condition is met, ensuring that the extracted components have real physical meaning. Step S2.5: Extract the candidate components that satisfy the IMF condition as IMF components, and subtract the IMF components from the original signal of the same signal to obtain the remaining signal, realizing the step-by-step extraction of the signal layer. Step S2.6: Use the remaining signal as a new input to repeat the decomposition process until the remaining signal changes monotonically and cannot be decomposed further, obtaining several IMF components and several residual components, completing the full-scale decomposition of the signal.

[0105] Step S3: Perform Hilbert-Huang transform on each intrinsic mode function component and each residual component to calculate the time-frequency characteristics, so as to extract the instantaneous frequency, instantaneous amplitude, and marginal spectral energy density.

[0106] Preferably, the core objective of this step is to perform component-wise Hilbert transform and time-frequency feature calculation on all intrinsic mode function components and all residual components output in step S2, converting the single-component signal in the time domain into a feature expression of the time-frequency joint distribution, extracting the instantaneous frequency and instantaneous amplitude that can accurately reflect the transient changes in the state of the packaging machine, as well as the marginal spectral energy density that characterizes the energy accumulation characteristics of the frequency dimension, so as to provide a high-dimensional feature set with strong fault sensitivity for subsequent time-series feature learning.

[0107] Furthermore, step S3 specifically includes the following steps: Step S3.1: Perform Hilbert transform on each intrinsic mode function component and each residual component, and calculate the real and imaginary parts of the analytic signal through convolution operation.

[0108] Preferably, the core of this step is to perform Hilbert transform on each intrinsic mode function component and each residual component, construct the corresponding analytic signal through convolution operation, separate the amplitude and phase information of the signal, and provide a mathematical basis for subsequent instantaneous feature calculation.

[0109] Specifically, the Hilbert transform is based on a continuous-time real-valued signal x(t) (i.e., the continuous form of a single intrinsic mode function component or residual component). The definition of is the convolution operation between the signal and 1 / (πt), and the formula is: In the formula, H[] is the Hilbert transform operator, ∗ is the convolution operator, and τ is the integration time variable.

[0110] This transformation is essentially a 90° phase-shifted all-pass filter that applies a -90° phase shift to the positive frequency components of the signal and a +90° phase shift to the negative frequency components, without changing the amplitude characteristics of the signal.

[0111] For the engineering implementation of discrete digital signals, the discrete time-series component sequence x(n), n=0,1,2,...,N-1 (N is the length of the component sequence, exactly the same as the length of the preprocessed original signal), output in step S2 can be implemented using the Fast Hilbert Transform based on the Fast Fourier Transform (FFT), avoiding the high computational cost of time-domain convolution. The specific execution flow is as follows: (1) Perform an N-point FFT transform on the discrete component sequence x(n) to obtain the frequency domain sequence X(k), k=0,1,2,...,N-1, where k is the frequency index.

[0112] (2) Construct the frequency domain response function H(k) of the Hilbert transform, according to the following rules: If N is odd, the frequency domain response rule is adjusted as follows: H(k)=1 when k=0, H(k)=2 when 1≤k≤(N-1) / 2, and H(k)=0 when (N+1) / 2≤k≤N-1.

[0113] (3) Multiply the frequency domain sequence X(k) and the frequency domain response function H(k) point by point to obtain the transformed frequency domain sequence X'(k)=X(k)·H(k).

[0114] Perform an N-point inverse fast Fourier transform (IFFT) on X'(k) to obtain the discrete sequence z(n) of the analytic signal.

[0115] For the separation of the real and imaginary parts of the analytic signal, the analytic signal z(n) obtained by the above transformation is a complex sequence that strictly satisfies the condition that the real part is the original input component and the imaginary part is the Hilbert transform result of that component. The separation formula is as follows: In the formula, j is the imaginary unit, and x(n) is the real part of the analytic signal. This is the imaginary part of the analytic signal (i.e., the result of the Hilbert transform of the original components).

[0116] Preferably, the Hilbert transform and analytic signal construction operation can be performed on all intrinsic mode function components and all residual components output in step S2 one by one to generate a corresponding real part sequence and imaginary part sequence for each component.

[0117] Step S3.2: Calculate the amplitude of the analytic signal using the real and imaginary parts, and obtain the instantaneous frequency by taking the time derivative of the amplitude of the analytic signal.

[0118] Preferably, the core of this step is to solve the instantaneous phase sequence of the signal based on the real and imaginary parts of the analytical signal, eliminate phase jump errors by phase unwrapping, and finally obtain the instantaneous frequency by taking the time derivative of the phase sequence, so as to accurately characterize the frequency change characteristics of the signal at each time point.

[0119] For the instantaneous phase calculation of the analytic signal, the real part sequence x(n) and the imaginary part sequence of the analytic signal corresponding to each component can be calculated. The instantaneous phase ϕ(n) of the signal is calculated point by point using the following formula: In the formula, arctan2() is the four-quadrant arctangent function, which can output phase values ​​covering the interval [-π,π], avoiding the phase ambiguity problem of the single-quadrant arctangent function.

[0120] In the phase unwrapping process, due to the output range limitation of the arctan2 function, the original phase sequence exhibits jumps every 2π. Direct differentiation leads to significant abnormal spikes in the instantaneous frequency. Therefore, phase unwrapping is necessary to eliminate these phase jumps and obtain a continuous, unwrapped phase sequence. This step employs a standard one-dimensional phase unwrapping algorithm. The core rule is: for the phase difference between two adjacent time points, if the absolute value of the difference is greater than π, then all subsequent phase values ​​are compensated by ±2π, with the compensation direction opposite to the jump direction. The specific quantization rules are as follows: In the formula, ϕ unwrapped (n) is the unwrapped continuous phase sequence, and round() is the rounding function.

[0121] For the instantaneous frequency calculation of discrete sequences, the physical definition of instantaneous frequency is the derivative of instantaneous phase with respect to time. For discrete time-series signals, the central difference method is used to realize the numerical differentiation of phase, balancing computational accuracy and noise resistance. At the same time, the sampling frequency is combined to complete the frequency dimension conversion. The specific calculation formula is as follows: (1) For the intermediate points of the sequence (1≤n≤N-2), the phase derivative is calculated using the central difference method: .

[0122] (2) For the first and last endpoints of the sequence (n=0 and n=N-1), forward difference and backward difference are used respectively for calculation: , .

[0123] (3) Combining the physical relationship between phase and frequency (angular frequency ω=2πf=dϕ / dt), the dimension conversion of instantaneous frequency is completed to obtain the instantaneous frequency sequence in Hz. In the formula, Δt is the sampling time interval of the signal, Δt = 1 / f s f s The sampling frequency of the corresponding signal (strictly consistent with the calibrated sampling frequency of each signal in step S1).

[0124] Preferably, based on the physical characteristics of different types of signals, an effective range threshold for the instantaneous frequency can be set to eliminate outliers caused by numerical calculations. (1) Vibration signal (sampling frequency 20kHz): The effective frequency range is calibrated to 0Hz to 10kHz (Nyquist frequency). Values ​​outside the range are replaced by linear interpolation of adjacent effective points.

[0125] (2) Current signal (sampling frequency 1kHz): The effective frequency range is calibrated to 0Hz to 500Hz.

[0126] (3) Temperature and speed signals (slow-changing signals): The effective frequency range is calibrated to 0Hz to 50Hz.

[0127] Preferably, the phase calculation, unwrapping, differentiation, and frequency transformation are performed on the real and imaginary part sequences of each component to generate the instantaneous frequency sequence corresponding to each component. The equal-length instantaneous frequency sequence corresponding to each component is output as one of the inputs for subsequent marginal spectral energy density calculations and is also included in the final feature set.

[0128] Step S3.3: Calculate the magnitude of the analytic signal using the real and imaginary parts to obtain the instantaneous amplitude.

[0129] Preferably, the core of this step is to solve the instantaneous amplitude sequence of the signal based on the real and imaginary parts of the signal, and to characterize the energy intensity and impact characteristics of the signal at each time point. This is the core sensitive feature for identifying mechanical faults and abnormal electrical loads of equipment.

[0130] In the mathematical definition and calculation of instantaneous amplitude, the magnitude of the analytic signal is the corresponding instantaneous amplitude. For discrete sequences, the square root of the sum of the squares of the real and imaginary parts of the analytic signal is calculated point by point. The calculation formula is as follows: In the formula, A(n) is the instantaneous amplitude at the nth time point, z(n) is the analytic signal sequence, and x(n) is the real part of the analytic signal. To analyze the imaginary part of the signal.

[0131] Preferably, the physical meaning of the instantaneous amplitude in this step is explained. The instantaneous amplitude sequence is the envelope of the original single-component signal, which can accurately capture the amplitude modulation characteristics of the signal. For example, the periodic impact signal generated by the failure of the outer ring of the packaging machine bearing will have its fault characteristics directly reflected in the instantaneous amplitude modulation law of the corresponding IMF component, which has a stronger fault identification degree compared with the original time domain signal.

[0132] Among them, for the non-negativity constraint, since the instantaneous amplitude is the result of the modulus calculation, it naturally satisfies non-negativity. For the extremely small negative values ​​(floating-point precision error) caused by numerical calculation, they are uniformly set to zero to ensure the physical rationality of the feature.

[0133] Preferably, the instantaneous amplitude calculation is performed on the real and imaginary part sequences of each component to generate an equal-length instantaneous amplitude sequence corresponding to each component. The instantaneous amplitude sequence corresponding to each component is output as one of the inputs for subsequent marginal spectral energy density calculation, and is also included in the final feature set.

[0134] Step S3.4: Integrate the instantaneous frequency and instantaneous amplitude over time to obtain the marginal spectral energy density.

[0135] The core of this step is to construct the Hilbert time spectrum based on the instantaneous frequency and instantaneous amplitude of each component, and then obtain the marginal spectral energy density, which characterizes the cumulative energy distribution in the frequency dimension, through integration in the time dimension, thus quantifying the energy contribution of different frequency components over the entire signal duration.

[0136] In the construction of the Hilbert time spectrum, the Hilbert time spectrum H(f,t) of a single component is constructed based on the instantaneous frequency sequence f(n) and the instantaneous amplitude sequence A(n). Its physical meaning is the instantaneous energy amplitude of the signal at time t and frequency f. In the discrete form, each time point n corresponds to a frequency point f(n) and an energy amplitude A(n), that is: H(f(n),n)=A(n). The Hilbert time spectra of all components are superimposed to obtain the joint Hilbert time spectrum of the entire signal, which completely covers the energy distribution of the signal across the entire time and frequency range.

[0137] Among them, the mathematical definition and discrete implementation of the marginal spectrum, the marginal spectrum h(f) is defined as the integral of the Hilbert time spectrum over the entire time dimension, representing the cumulative energy of a certain fixed frequency f over the entire time history. The continuous form of the calculation formula is as follows: In the formula, T represents the total duration of the signal. For discrete-time signals, the integration operation is converted into a summation operation with equal time intervals, while simultaneously discretizing the frequency axis and accumulating energy. The specific execution flow is as follows: (1) Frequency axis discretization calibration: Based on the sampling frequency and effective frequency range of each signal, set the resolution and coverage of the frequency axis. For example, the frequency axis range of the vibration signal is 0Hz to 10kHz, and the frequency resolution is calibrated to 10Hz. Generate an equally spaced discrete frequency interval sequence fb,b=0,1,2,...,B-1, where B is the total number of frequency intervals.

[0138] (2) Energy accumulation calculation: For each time point in the whole time series, the instantaneous amplitude A(n) of that point is accumulated into the total energy corresponding to the frequency interval to which its instantaneous frequency f(n) belongs.

[0139] (3) Marginal spectral energy density normalization: For the cumulative energy in each frequency interval, divide by the product of the total duration and the frequency resolution to obtain the energy density per unit frequency and per unit time, eliminating the influence of signal duration and frequency resolution on the numerical value, and finally obtain the marginal spectral energy density sequence E(f b ).

[0140] Preferably, the marginal spectral energy densities of all intrinsic mode function components and residual components are superimposed on a frequency-interval basis to obtain the comprehensive marginal spectral energy density of the entire signal, which fully characterizes the full-band energy distribution characteristics of the original multi-physical quantity signal. The instantaneous frequency, instantaneous amplitude, and marginal spectral energy density sequence of the entire frequency range corresponding to each time point are output. All features are integrated into a complete time-frequency feature set, which serves as the sole input data source for the subsequent step S4.

[0141] For example, this step can be implemented using the following pseudocode: # Pseudocode for Hilbert-Huang Transform Time-Frequency Feature Extraction import numpy as np from scipy.signal import hilbert def hht_feature_extraction(imf_list, fs): """ Input parameters: imf_list: The component list output in step S2, where each element is a discrete sequence of an IMF component or residual component. fs: The sampling frequency of the corresponding signal, consistent with the sampling frequency calibrated in step S1. Output parameters: inst_freq_all: The set of instantaneous frequency sequences of all components. inst_amp_all: The set of instantaneous amplitude sequences for all components. marginal_spectrum: The marginal spectral energy density sequence after full component integration freq_axis: The frequency axis sequence corresponding to the marginal spectrum """ # Initialize basic variables N = len(imf_list[0]) # Length of the signal sequence, all components are the same length as the original signal dt = 1 / fs total_time = N * dt inst_freq_all = [] inst_amp_all = [] # Initialize frequency axis and energy accumulation array freq_resolution = 10 # Frequency resolution is calibrated to 10Hz max_freq = fs / 2 # Upper limit of Nyquist frequency freq_axis = np.arange(0, max_freq + freq_resolution, freq_resolution) energy_accum = np.zeros_like(freq_axis) # Perform Hilbert transform and time-frequency characteristic calculation on a component-by-component basis for component in imf_list: # 3.1 Constructing analytic signals using Hilbert transform, separating the real and imaginary parts. analytic_signal = hilbert(component) x_real = np.real(analytic_signal) x_imag = np.imag(analytic_signal) # 3.3 Calculating the Instantaneous Amplitude Sequence inst_amp = np.abs(analytic_signal) inst_amp_all.append(inst_amp) # 3.2 Calculation of instantaneous phase, phase unwrapping, and instantaneous frequency inst_phase = np.arctan2(x_imag, x_real) unwrapped_phase = np.unwrap(inst_phase) # One-dimensional phase unwrapping # Numerical differentiation to calculate the phase derivative, adapting to the start and end points inst_freq = np.zeros_like(unwrapped_phase) inst_freq[1:-1] = (unwrapped_phase[2:]- unwrapped_phase[:-2]) / (2 * dt) inst_freq[0] = (unwrapped_phase[1]- unwrapped_phase[0]) / dt inst_freq[-1] = (unwrapped_phase[-1]- unwrapped_phase[-2]) / dt # Convert angular frequency to Hz unit, constrain the effective frequency range inst_freq = inst_freq / (2 * np.pi) inst_freq = np.clip(inst_freq, 0, max_freq) inst_freq_all.append(inst_freq) # 3.4 Marginal Spectral Energy Accumulation for n in range(N): current_freq = inst_freq[n] current_amp = inst_amp[n] # Match the frequency range to which the current frequency belongs bin_idx = np.argmin(np.abs(freq_axis - current_freq)) energy_accum[bin_idx] += current_amp * dt # Calculate the normalized marginal spectral energy density marginal_spectrum = energy_accum / (total_time * freq_resolution) return inst_freq_all, inst_amp_all, marginal_spectrum, freq_axis Beneficial effects of steps S3.1 to S3.4: A complete Hilbert-Huang transform time-frequency feature extraction scheme is constructed through systematic implementation of steps including Hilbert transform, calculation of the real and imaginary parts of the analytic signal, calculation of the amplitude and time derivative of the analytic signal, calculation of the magnitude of the analytic signal, and integration of time-frequency features. This scheme can perform accurate time-frequency analysis on each intrinsic mode function component and each residual component, extracting multidimensional features including instantaneous frequency, instantaneous amplitude, and marginal spectral energy density, comprehensively characterizing the time-varying properties of the signal, fully mining equipment operating status information, and providing a rich and effective feature foundation for subsequent time-series feature learning and health assessment. Specifically, step S3.1 involves performing a Hilbert transform on each intrinsic mode function component and each residual component, and calculating the real and imaginary parts of the analytic signal through convolution operations to construct a complex-form analytic signal representation. Step S3.2 involves calculating the amplitude of the analytic signal using the real and imaginary parts, and obtaining the instantaneous frequency by taking the time derivative of the amplitude, accurately reflecting the instantaneous phase change rate of the signal and capturing the dynamic evolution of the signal's frequency components. Step S3.3 involves calculating the magnitude of the analytic signal using the real and imaginary parts to obtain the instantaneous amplitude, accurately describing the instantaneous intensity change of the signal and characterizing the instantaneous energy level of vibration, temperature, current, and rotational speed signals. Step S3.4 involves integrating the instantaneous frequency and instantaneous amplitude over time to obtain the marginal spectral energy density, comprehensively characterizing the signal's distribution characteristics across different frequency components from an energy perspective, and providing frequency domain energy features for equipment status identification.

[0142] Step S4: Input all instantaneous frequencies, all instantaneous amplitudes, and all marginal spectral energy densities into the Long Short-Term Memory network for temporal feature learning. Update the cell state and hidden state through the gating mechanism of forget gate, input gate, and output gate to obtain the temporal feature vector.

[0143] Preferably, the core objective of this step is to perform deep temporal feature learning on all instantaneous frequencies, all instantaneous amplitudes, and all marginal spectral energy densities output in step S3. By using the gating mechanism of Long Short-Term Memory (LSTM) networks, the vanishing and exploding gradient problems in traditional recurrent neural network long sequence learning are solved. This accurately captures the gradual change pattern of the operating status features of the GDX500 packaging machine in long time series and the dependencies across time steps, filters out redundant temporal information unrelated to the equipment health status, and finally outputs a high-dimensional temporal feature vector containing the trend of equipment health degradation, providing input for the key feature weighting of the subsequent attention mechanism.

[0144] Furthermore, step S4 specifically includes the following steps: Step S4.1: Concatenate all instantaneous frequencies, all instantaneous amplitudes, and all marginal spectral energy densities in chronological order to form a feature sequence.

[0145] Preferably, the core of this step is to convert the multi-dimensional, heterogeneous time-frequency features output in step S3 into a standardized feature sequence that meets the input format requirements of the LSTM network and is strictly aligned in time, so as to ensure that the network can complete the time-dependent learning in chronological order.

[0146] Specifically, regarding the alignment of feature dimensions with the time series reference, all instantaneous frequency sequences and instantaneous amplitude sequences output in step S3 maintain a time step length and millisecond-level timestamp reference that are completely consistent with the original signal after preprocessing in step S1. All component feature sequences of a single signal are strictly of equal length, with no time misalignment. The feature sequences of vibration signals, current signals, temperature signals, and rotational speed signals are all unified to the same total number of time steps N through the timestamp alignment mechanism in step S1, ensuring that the features of different physical quantities strictly correspond at the same time step.

[0147] Specifically, for the construction of single-time-step feature vectors, for the t-th time step (t=0,1,2,...,N-1) in the time series, three types of features are integrated: the first type is instantaneous features, which extract the instantaneous frequency values ​​and instantaneous amplitude values ​​corresponding to all intrinsic mode function components and residual components at this time step to form a transient time-frequency feature group; the second type is global energy features, which broadcast the full-sequence marginal spectrum energy density sequence output in step S3 to each time step as a supplementary feature of global frequency domain energy distribution, ensuring that the single-time-step features simultaneously contain transient change information and full-cycle energy distribution information; the third type is operating condition auxiliary features, which synchronously incorporate the packaging machine spindle speed signal collected in step S1 as an auxiliary feature for operating condition normalization, eliminating the influence of different production speeds on feature amplitude.

[0148] In particular, for the secondary standardization of features, in order to eliminate the impact of the difference in numerical magnitude of different feature dimensions on network learning, the mean and standard deviation of the concatenated full feature dimensions are used to perform Z-score standardization again to ensure that the mean of all feature dimensions is 0 and the standard deviation is 1, so as to avoid the large-scale features from dominating the network gradient update.

[0149] In the temporal sequence concatenation and tensor transformation, the standardized D-dimensional feature vectors (D being the total feature dimension of a single time step) of each time step are arranged sequentially according to the time axis to form a two-dimensional temporal feature sequence of shape [N,D]. Then, according to the input requirements of the LSTM network, it is converted into a three-dimensional tensor format [batch_size, seq_len, feature_dim], where batch_size is the batch size, seq_len=N is the length of the time step sequence, and feature_dim=D is the feature dimension of a single time step, which serves as the sole input to the LSTM network.

[0150] Preferably, for long time-series data from continuous operation of the packaging machine, a non-overlapping sliding window strategy is adopted to divide the ultra-long sequence into fixed-length subsequences. The window length is set to 1000 time steps to ensure network computational efficiency while retaining a sufficiently long temporal dependency context information. The output standardized and temporally aligned LSTM input feature sequence serves as the sole input for subsequent gating computation and temporal feature learning.

[0151] Step S4.2: Input the feature sequence into the Long Short-Term Memory (LSTM) network, and calculate the forgetting gate coefficient through the forgetting gate of the LSM network to control the degree of retention of the cell state at the previous time step.

[0152] Preferably, the core of this step is to calculate the forgetting gate coefficient at the current moment through the forgetting gate structure of LSTM, quantify the degree of retention and discarding of historical information in the cell state at the previous moment, filter out redundant historical information that is irrelevant to the health status of the device, and retain long-term dependent features that are effective for health assessment.

[0153] Preferably, the input to the forget gate consists of two parts: one is the hidden state h output by the network at the previous time step. t-1 Secondly, the input feature vector x at the current moment. t (i.e., the feature vector of the feature sequence output in step S4.1 at time step t), the two are concatenated and then input into the fully connected layer of the forget gate.

[0154] Preferably, the core operation of the forget gate is a fully connected linear transformation + Sigmoid activation function. The Sigmoid function can compress the output value to the [0,1] interval, which naturally fits the physical meaning of the gating coefficient (0 represents completely discarding the corresponding information, and 1 represents completely retaining the corresponding information). The calculation formula is as follows: f t =σ(W f ·[h t-1 ,x t ]+b f In the formula, f t W is the forgetting gating coefficient at time step t, with the same dimension as the cell state dimension; σ() is the sigmoid activation function; W f Here is the weight matrix for the forget gate, with dimensions [hidden_dim, hidden_dim + feature_dim], where hidden_dim is the number of hidden units in the LSTM layer; b f Here is the bias vector for the forget gate, with dimension [hidden_dim]; [h t-1 ,x t This is an operation that concatenates the hidden state from the previous time step with the input features from the current time step, completing the merging along the feature dimension.

[0155] Preferably, for the health monitoring scenario of packaging machines, the weight matrix adopts the Xavier uniform initialization method, and the bias vector is initialized to 0.1 to avoid the gradient vanishing problem during network initialization. Simultaneously, a Dropout layer is added after the fully connected layer, with the Dropout rate calibrated to 0.2, to suppress network overfitting and improve the model's generalization ability to complex industrial conditions. The current forgetting gate coefficient f is output. t It serves as one of the core inputs for subsequent cell state updates.

[0156] Step S4.3: Calculate the input gating coefficient and candidate cell state through the input gate of the Long Short-Term Memory network to determine the new information at the current moment.

[0157] Preferably, the core of this step is to calculate the input gating coefficient and candidate cell state at the current time through the input gating structure, thereby quantifying the degree of retention of the new input information at the current time and the newly learned candidate state information, providing new information sources for updating the cell state. Specifically, for the calculation of the input gating coefficients, the input to the input gate is exactly the same as that of the forget gate, which is the hidden state h from the previous time step. t-1 Input feature x at the current time t The concatenated vector is processed by a Sigmoid activation function, which outputs a gating coefficient in the 0-1 interval to control the degree of new information being written at the current time. The calculation formula is: i t =σ(W i ·[h t-1 ,x t ]+b i In the formula: i t W represents the input gating coefficient at time step t; i Let b be the weight matrix of the input gate. i This is the bias vector of the input gate, with dimensions consistent with the weights and biases of the forget gate.

[0158] For candidate cell state calculation, the tanh activation function is used to generate the candidate cell state at the current time step, representing the new content in the input information that can be used to update the cell state at the current time step. The tanh function compresses the output value to the interval [-1, 1] to match the numerical range of the cell state. The calculation formula is as follows: In the formula: Let W be the candidate cell state at time step t; tanh() is the hyperbolic tangent activation function; W C Let b be the weight matrix of the candidate cell states. C This is the corresponding bias vector, whose dimensions match the aforementioned weights and biases.

[0159] Preferably, the scene adaptation description is as follows: the input gating coefficient is coordinated with the candidate cell state, it Decision made Which information will be written into the final cell state? For packaging machine health monitoring scenarios, this mechanism can effectively capture transient abnormal changes in equipment status and write fault-related new feature information into the cell state, achieving the memory and retention of early degradation features. Output: Output the input gating coefficient i at the current moment. t With candidate cell state It serves as the core input for subsequent cell state updates.

[0160] Step S4.4: Update the current cell state using the forgetting gating coefficient and the input gating coefficient.

[0161] Preferably, this step is the core memory update step of the LSTM network. Based on the forgetting gate coefficient, input gate coefficient, and candidate cell state output from the previous steps, the cell state at the current time step is iteratively updated. The cell state acts as the "long-term memory conveyor" of the LSTM, completing information updates only through linear interactions. This effectively preserves cross-time step dependency information in long sequences and avoids the gradient vanishing problem.

[0162] The cell state update consists of two parts: forgetting the historical cell state and writing the new information at the current moment. The two are linearly superimposed to obtain the latest cell state at the current moment. The calculation formula is as follows: In the formula, C t The cell state after the update at time step t is the carrier of the network's long-term memory; C t-1 is the cell state at time step t-1, i.e., the long-term memory information of the previous time step; ∗ is the Hadamard product (element-by-element multiplication), ensuring that the gating coefficients and the corresponding elements of the cell state are precisely matched.

[0163] Preferably, for the health monitoring needs of the GDX500 packaging machine, this update mechanism can achieve two core capabilities: one is through f t Preserve the baseline characteristics of long-term healthy operation of the equipment, and filter out irrelevant interference caused by random noise and fluctuations in operating conditions; secondly, through i t By writing the subtle abnormal features of early equipment degradation into long-term memory, we can continuously track the information of early signs of failure and prevent these subtle features from being overwhelmed by long sequences of normal data.

[0164] Preferably, a gradient clipping mechanism is incorporated during the cell state update process, with a gradient clipping threshold set to 1.0 to avoid gradient explosion caused by extreme outliers and ensure numerical stability during network training and inference. The updated cell state C at the current time step is output. t On the one hand, it serves as the historical information input for the calculation of the forget gate and input gate at the next moment; on the other hand, it serves as the core input for the calculation of the output gate at the current moment.

[0165] Step S4.5: Calculate the output gating coefficients through the output gate of the Long Short-Term Memory network to determine the hidden state at the current moment based on the updated cell state.

[0166] Preferably, the core of this step is to calculate the output gating coefficient at the current time through the output gate structure, and control how much information in the updated cell state will be output to the hidden state at the current time. The hidden state is the short-term temporal feature learned by the network at the current time step.

[0167] Specifically, for calculating the output gating coefficients, the input to the output gate is also the hidden state h from the previous time step. t-1 Input feature x at the current time t The concatenated vector is passed through a Sigmoid activation function to output the gating coefficients in the 0-1 interval. The calculation formula is: o t =σ(W o ·[h t-1 ,x t ]+b o In the formula, o t W represents the output gating coefficient at time step t. o Let b be the weight matrix of the output gate. o This is the bias vector of the output gate, and its dimension matches the weights and biases of the aforementioned gate.

[0168] Specifically, for the calculation of the hidden state at the current moment, the updated cell state C is first calculated. t Perform tanh activation to compress the value to the interval [-1, 1], then multiply it element-wise with the output gating coefficient to obtain the hidden state h at the current time step. t The calculation formula is: h t =o t ∗tanh(C t In the formula, h t The hidden state at time step t is a compressed feature representation of the network for the current moment and historical time sequence information, containing time-dependent information related to the health status of the device.

[0169] Preferably, the hidden state h t It also integrates long-term memory information (from cell state C). t The hidden state at each time step, along with the transient input information at the current moment, fully preserves all valid temporal information from the beginning of the sequence to the current moment. It is the core object for subsequent weighted processing by the attention mechanism. The output is the hidden state h at the current moment. t On the one hand, it serves as the input of the historical hidden state for the next time-gated computation; on the other hand, it incorporates the hidden state sequence of the entire time step, providing a foundation for the final construction of the temporal feature vector.

[0170] Step S4.6: Repeat steps S4.2 to S4.5 until all time steps have been processed to obtain a time-series feature vector with temporal dependencies.

[0171] Preferably, the core of this step is to complete the iterative calculation of all time steps of the entire time series, integrate the hidden states of all time steps, generate a time series feature vector with complete time series dependencies, and provide a unique input for the subsequent step S5.

[0172] For the full-time-step iteration rule, starting from the first time step t=0 of the time series, steps S4.2 to S4.5 of the gating calculation and state update process are executed sequentially until the last time step t=N-1 of the sequence is processed; during the iteration process, the cell state C output at each time step is... t Hidden state h t It is automatically used as the input parameter for the next time step t+1, completing the continuous transmission and accumulation of timing information.

[0173] Preferably, based on the operating characteristics and data features of the GDX500 packaging machine, the core hyperparameters of the LSTM network are fixedly calibrated to ensure a balance between feature learning effectiveness and computational efficiency: the number of hidden layer units (hidden_dim) is calibrated to 128 to balance feature expressive power and computational overhead; the number of stacked LSTM layers is calibrated to 2 layers, with the first LSTM being a bidirectional LSTM that simultaneously captures forward and backward temporal dependencies and outputs the hidden state at each time step, and the second LSTM being a unidirectional LSTM that further learns higher-order temporal features to improve feature robustness; the maximum sequence length is calibrated to 5000 time steps to cover the complete operating data of at least 10 spindle rotation cycles of the packaging machine, ensuring the capture of periodic fault features.

[0174] Preferably, after completing the full time-step iteration, the hidden states of all time steps are collected to form a temporal feature matrix of shape [N, hidden_dim]. This matrix is ​​a set of temporal feature vectors with temporal dependencies, where each row corresponds to the temporal feature of a time step, thus fully preserving the temporal sequence relationship and long-term dependency information of the features.

[0175] Beneficial effects of steps S4.1 to S4.6: A complete temporal feature learning scheme is formed through cyclical execution of steps such as feature sequence concatenation, temporal learning via Long Short-Term Memory (LSTM) networks, calculation of forgetting gate input-output gate gating mechanisms, and updating of cell states and hidden states. This scheme can automatically capture temporal dependencies and long-term dependency patterns in feature sequences, extracting deep features with temporal contextual information from instantaneous frequency, instantaneous amplitude, and marginal spectral energy density sequences. This effectively improves the ability to identify equipment performance degradation trends and early fault symptoms, providing feature vectors containing dynamic temporal information for subsequent health assessments. In this process, step S4.1 involves concatenating all instantaneous frequencies, instantaneous amplitudes, and marginal spectral energy densities in chronological order to form a feature sequence, constructing a unified high-dimensional temporal input and providing a standard data format for Long Short-Term Memory (LSTM) network processing. Step S4.2 inputs the feature sequence into the LSM network, calculating the forgetting gate coefficient to control the retention of the cell state from the previous time step, thus achieving intelligent filtering of historical information and maintenance of long-term memory. Step S4.3 calculates the input gate coefficient and candidate cell states through the input gate to determine the newly added information at the current time step, realizing the contribution of the current input to updating the network state. Step S4.4 updates the current cell state using the forgetting gate coefficient and the input gate coefficient, balancing the influence ratio of historical information and current information to ensure dynamic optimization of the network state. Step S4.5 calculates the output gate coefficient through the output gate, determines the hidden state at the current time step based on the updated cell state, and outputs feature information useful for subsequent tasks. Step S4.6 repeats steps S4.2 to S4.5 until all time steps are processed, obtaining a temporal feature vector with temporal dependencies, and fully extracting long-term dependency information from the feature sequence.

[0176] Step S5: Input the temporal feature vector into the attention mechanism to calculate the weight coefficient of each time step, and calculate the weighted aggregated temporal feature vector by weighted summation and aggregating the feature information of key time steps.

[0177] Preferably, the core objective of this step is to perform feature filtering and weighted aggregation of key time steps on the time-series feature vector output in step S4. By using an attention mechanism, the contribution weight of features at different time steps to the equipment health assessment is automatically learned, amplifying the feature influence of key time steps such as fault precursors and sudden state changes, while weakening meaningless features of stable operating condition time steps. Finally, a weighted aggregated time-series feature vector that integrates core time-series information is generated, providing high-value feature input for the subsequent health modeling of the support vector data description algorithm.

[0178] Furthermore, step S5 specifically includes the following steps: Step S5.1: Input the temporal feature vector into the attention mechanism and calculate the relevance score between the query vector and the key vector through dot product operation.

[0179] Preferably, the core of this step is to construct the query and key vectors of the attention mechanism, quantify the correlation between the temporal feature vector of each time step and the global features through dot product operation, provide a basic score for subsequent weight coefficient calculation, and accurately identify key time steps that are of great significance to the determination of device health status.

[0180] In terms of feature vector dimension and format adaptation, the temporal feature vector output in step S4 is a matrix of shape [seq_len, hidden_dim], where seq_len is the total number of time steps and hidden_dim is the dimension of the LSTM hidden layer (calibrated to 128). This matrix is ​​used as the key vector and value vector of the attention mechanism, and the dimensions of the two are completely consistent, both being [seq_len, hidden_dim].

[0181] Specifically, for the construction of the query vector, a globally shared learnable query vector with dimensions [1, hidden_dim] is constructed for the single-sample time-series analysis scenario of health monitoring of packaging machine equipment. This vector is automatically learned through network training to represent the global features of the equipment health status. The initial value is initialized with a normal distribution (mean 0, standard deviation 0.01) and updated iteratively with the model during training to ensure accurate matching of equipment health-related features.

[0182] For the dot product correlation calculation, a batch dot product operation can be performed between the query vector and the key vectors of all time steps to quantify the similarity between the key vector and the query vector at each time step, i.e., the correlation between the feature at that time step and the global features of device health. The calculation formula is as follows: In the formula, S t Let Q be the relevance score at time step t, Q be the query vector ([1, hidden_dim]), and K be the relevance score at time step t. t Let t be the key vector at time step t ([1, hidden_dim]), · is the matrix dot product operation, and T is the matrix transpose.

[0183] Preferably, to avoid excessively large dot product values ​​due to the high dimension of `hidden_dim`, which could lead to gradient saturation in subsequent exponential operations, all correlation scores are scaled by a factor equal to the square root of the hidden layer dimension. The formula for the scaled correlation score is: In this scenario, hidden_dim=128, so the scaling factor is approximately 11.31, ensuring that the relevance scores are distributed within a reasonable range.

[0184] Preferably, the above dot product and scaling operations are performed sequentially on the key vectors of all time steps to generate a correlation score sequence of shape [1, seq_len], which fully covers the feature correlation quantization results of all time steps. The scaled full-time-step correlation score sequence is then output. It serves as the sole input for subsequent weighting coefficient calculations.

[0185] Step S5.2: Perform an exponential operation on the correlation score and divide it by the sum of the exponential operation results of all time steps to obtain the weight coefficient for each time step.

[0186] Preferably, the core of this step is to index and normalize the correlation score, and use the Softmax function to convert the correlation score into a weight coefficient with a sum of 1, so as to realize the quantitative allocation of the importance of features at different time steps. The larger the weight coefficient, the higher the contribution of the feature at that time step to the assessment of equipment health.

[0187] Specifically, to enhance feature differentiation through exponential operations, a natural exponential operation is performed on the relevance score S~t at each scaled time step to amplify the difference between high and low relevance scores, thereby strengthening the feature discriminative power at key time steps. The calculation formula is as follows: In the formula, E t Let be the exponential score at time step t. The exponential operation ensures that the score is non-negative, which meets the physical requirements of the weighting coefficient.

[0188] Specifically, for the global normalization calculation of weight coefficients, the global sum of the exponential scores of all time steps is calculated. The exponential score of each time step is then divided by this sum to obtain the normalized weight coefficients. This ensures that the sum of the weight coefficients of all time steps is 1, thus achieving the global allocation of feature importance. The calculation formula is as follows: In the formula, αt is the weight coefficient at the t-th time step, and its value ranges from [0,1]. .

[0189] Preferably, to avoid overflow issues caused by excessively large values ​​during exponential operations, a numerical shift is performed on the correlation score sequence, i.e., the maximum value in the sequence is subtracted first, and then the exponential operation is performed. The optimized formula is as follows: This operation does not change the relative magnitude of the weighting coefficients, effectively ensuring the stability of numerical calculations and making it suitable for feature analysis under complex industrial conditions.

[0190] Preferably, the minimum threshold for the weight coefficients can be set to 1 / seq_len (i.e., the average weight) to avoid assigning zero weights to some time steps due to excessively low relevance scores, which would lead to the loss of effective feature information and ensure that feature information from all time steps is included in subsequent aggregation calculations. The output is the weight coefficient sequence α=[α1,α2,...,α] for all time steps. seq_lenThe shape is [1, seq_len], which strictly corresponds to the time step of the temporal feature vector.

[0191] Step S5.3: The feature information of each time step in the time series feature vector is weighted and summed using weight coefficients to obtain the weighted aggregated time series feature vector.

[0192] Preferably, the core of this step is to perform element-wise weighted summation of the weight coefficients and the temporal feature vectors of the corresponding time steps, thereby amplifying and fusing the features of key time steps, and finally generating a weighted aggregated temporal feature vector with fixed dimensions, thus eliminating the influence of the time step dimension on subsequent algorithms.

[0193] Specifically, for the dimensional expansion of the weight coefficients, the weight coefficient sequence of shape [1, seq_len] is expanded to a matrix of shape [seq_len, hidden_dim], ensuring that the weight coefficients at each time step are matched one-to-one with each dimension of the corresponding temporal feature vector, thus guaranteeing dimensional consistency in the weighting operation. The expansion rule is to increase the weight coefficient α of a single time step... t Copy to all hidden_dim dimensions of the feature vector at this time step.

[0194] Specifically, for the element-wise weighted operation, the expanded weight coefficient matrix and the time-series feature vector matrix output in step S4 are subjected to an element-wise Hadamard product operation. This weights the feature values ​​for each time step and each dimension, amplifying the feature values ​​of high-weight time steps and weakening the feature values ​​of low-weight time steps. The calculation formula is as follows: In the formula, V is the weighted feature vector ([1, hidden_dim]) at time step t, ∗ is the Hadamard product operation, and V t Let be the original time series feature vector at the t-th time step.

[0195] Specifically, for the global summation to generate aggregated features, a global dimensional summation operation is performed on the weighted feature vectors of all time steps. The weighted feature vectors of seq_len time steps are merged into a fixed-dimensional feature vector with dimensions [1, hidden_dim], which is the weighted aggregated time-series feature vector. The calculation formula is as follows: In the formula, V agg The weighted aggregated temporal feature vector has dimensions [1, 128]. This vector integrates the effective feature information of all time steps and highlights the feature contributions of key time steps.

[0196] Preferably, L2 normalization can be performed on the aggregated feature vectors to map the feature values ​​onto a unit hypersphere, eliminating the influence of feature magnitude differences on subsequent support vector data description algorithms. The normalization formula is as follows: In the formula, ||Vagg||2 is the L2 norm of the aggregated feature vector. This is the final weighted aggregated temporal feature vector. The output is a weighted aggregated temporal feature vector with dimensions [1, hidden_dim], which serves as the sole input data source for the subsequent step S6 device health modeling.

[0197] For example, this step can be implemented using the following pseudocode: # Pseudocode for feature weighted aggregation based on PyTorch's attention mechanism import torch import torch.nn as nn import torch.nn.functional as F class AttentionPooling(nn.Module): def __init__(self, hidden_dim=128): super(AttentionPooling, self).__init__() # Construct a learnable global query vector to adapt to global feature matching for packaging machine health monitoring. self.query = nn.Parameter(torch.randn(1, hidden_dim, dtype=torch.float32) * 0.01) self.hidden_dim = hidden_dim def forward(self, temporal_feature): """ Forward propagation function Input temporal_feature: The temporal feature vector output from step S4, with shape [seq_len, hidden_dim] Output: Weighted aggregated temporal feature vector, shape [1, hidden_dim] """ # 5.1 Calculate the relevance score using dot product and perform dimensional scaling. # query: [1, hidden_dim], temporal_feature(Key): [seq_len,hidden_dim] score = torch.matmul(self.query, temporal_feature.T) # Dot product operation, shape [1, seq_len] scale_score = score / torch.sqrt(torch.tensor(self.hidden_dim, dtype=torch.float32)) # Dimension scaling # 5.2 Calculate normalized weight coefficients using Softmax and perform numerical stability optimization. scale_score = scale_score - torch.max(scale_score) # Shift the value to prevent overflow exp_score = torch.exp(scale_score) # Exponentiation operation weight = exp_score / torch.sum(exp_score, dim=1, keepdim=True) # Normalization, shape [1, seq_len] # Minimum threshold constraint for weight coefficients to avoid zero weights min_weight = 1 / temporal_feature.shape[0] weight = torch.clamp(weight, min=min_weight) weight = weight / torch.sum(weight, dim=1, keepdim=True) # Renormalize # 5.3 Weighted summation for feature aggregation: dimensionality expansion + Hadamard product + global summation weight_expand = weight.unsqueeze(2).repeat(1, 1, self.hidden_dim) # Dimension expansion [1, seq_len, hidden_dim] temporal_feature_unsq = temporal_feature.unsqueeze(0) # shape[1,seq_len, hidden_dim] weighted_feature = weight_expand * temporal_feature_unsq # Element-wise weighting `agg_feature = torch.sum(weighted_feature, dim=1)` # Global summation, shape [1, hidden_dim # L2 normalization agg_feature = F.normalize(agg_feature, p=2, dim=1) return agg_feature # Example of Model Inference and Feature Aggregation if __name__ == "__main__": # The hyperparameters are strictly consistent with step S4 HIDDEN_DIM = 128 SEQ_LEN = 1000 # Total number of time steps, consistent with the length of the LSTM input sequence. # Initialize the attention aggregation model attention_model = AttentionPooling(hidden_dim=HIDDEN_DIM) attention_model.eval() # Simulated input: The temporal feature vector output from step S4, with shape [seq_len, hidden_dim] temporal_feature_vec = torch.randn(SEQ_LEN, HIDDEN_DIM, dtype=torch.float32) # Perform weighted feature aggregation to obtain the final aggregated features. with torch.no_grad(): weighted_agg_feature = attention_model(temporal_feature_vec) print("Weighted aggregation time series feature vector shape:", weighted_agg_feature.shape) # Output: torch.Size([1, 128]) Beneficial effects of steps S5.1 to S5.3: By leveraging the synergistic effects of attention mechanisms in calculating the correlation between query vectors and key vectors, performing correlation score exponent calculations and normalization, and weighted summation of weight coefficients, a complete temporal feature aggregation scheme is constructed. This scheme can automatically identify the importance of each time step in the temporal feature vector, highlighting key moment features that significantly contribute to equipment health status assessment, suppressing interference from irrelevant or redundant features, improving the focus and effectiveness of feature representation, and providing a compact feature representation containing key information for subsequent equipment health scoring. In step S5.1, the temporal feature vector is input into the attention mechanism. The relevance score between the query vector and the key vector is calculated through dot product operation to quantify the degree of correlation between features at different time steps and evaluate the relative importance of features at each time step to the current task. Step S5.2, the relevance score is exponentially calculated and divided by the sum of the exponential calculation results of all time steps to obtain the weight coefficient of each time step. The relevance score is converted into a probability distribution with values ​​ranging from zero to one to ensure that the sum of all weight coefficients is one, thereby achieving reasonable allocation and normalization of weights. Step S5.3, the feature information of each time step in the temporal feature vector is weighted and summed through the weight coefficient. Features at different times are aggregated according to their importance to form a weighted aggregated temporal feature vector that enhances the features of key time steps and suppresses the features of secondary time steps, thereby improving the discriminative ability and robustness of feature representation.

[0198] Step S6: Input the weighted aggregated temporal feature vector into the support vector data description algorithm to construct a hypersphere model of device health status. Calculate the distance from the sample to the center of the hypersphere by mapping to a high-dimensional feature space through a kernel function, and obtain the device health score.

[0199] Preferably, the core objective of this step is to perform unsupervised anomaly detection modeling on the weighted aggregated time-series feature vector output in step S5. A minimal hypersphere model containing only equipment health status samples is constructed using the Support Vector Data Description (SVDD) algorithm. The kernel function is used to map low-dimensional features to a high-dimensional nonlinear space to solve the nonlinear distribution problem of industrial equipment health features. Finally, the equipment health is quantified by the distance from the sample to the center of the hypersphere, generating an intuitive health score set.

[0200] Furthermore, step S6 specifically includes the following steps: Step S6.1: The weighted aggregated temporal feature vectors are mapped to a high-dimensional feature space using the support vector data description algorithm.

[0201] Preferably, the core of this step is to implicitly map the time-series feature vectors after low-dimensional weighted aggregation to a high-dimensional feature space through a kernel function, thereby avoiding the curse of dimensionality in explicit high-dimensional computation. At the same time, by utilizing the linear separability of the high-dimensional space, a more accurate hypersphere model of health status can be constructed.

[0202] Regarding the selection of the kernel function, considering the nonlinear and multimodal distribution characteristics of the health features of the GDX500 packaging machine, the Radial Basis Function (RBF) was chosen as the kernel function for SVDD. This kernel function has the following advantages: First, it can map low-dimensional nonlinear features to an infinite-dimensional Hilbert space, achieving linearization of arbitrarily complex distributions; second, only one kernel parameter needs to be adjusted, making the model complexity controllable; and third, it is sensitive to local features, making it suitable for capturing weak anomalies in the early stages of equipment degradation.

[0203] For the Gaussian RBF kernel function, for two sample vectors x in the low-dimensional space... i and x j (That is, the weighted aggregated temporal feature vector output in step S5), its inner product in the high-dimensional Hilbert space can be directly calculated using the RBF kernel function without explicit mapping, and the formula is: In the formula, ||xi-xj||2 is the L2 Euclidean distance between the two sample vectors, and σ is the bandwidth parameter of the RBF kernel, which controls the local scope of the kernel function.

[0204] For the calibration of kernel parameters, the optimal kernel parameters are determined by cross-validation. For the health monitoring scenario of packaging machines, the bandwidth parameter σ is calibrated to 0.5, which can achieve the optimal balance between local feature sensitivity and global model generalization. At the same time, in order to avoid overflow of kernel function numerical calculation, the L2 distance of all sample vectors is normalized, and the maximum distance is normalized to 1 to ensure that the kernel function value is distributed in a reasonable range.

[0205] Preferably, in the SVDD algorithm, there is no need to explicitly calculate the mapping vector ϕ(xi) in the high-dimensional space. All inner product and distance calculations in the high-dimensional space are indirectly completed through the RBF kernel function, which significantly reduces the computational complexity. After completing the implicit mapping preparation of the kernel function parameters and sample vectors, all subsequent high-dimensional space calculations are performed based on the RBF kernel function.

[0206] Step S6.2: Calculate the Euclidean distance matrix between sample vectors in the high-dimensional feature space using the Gaussian radial basis kernel function.

[0207] Preferably, the core of this step is to indirectly calculate the Euclidean distance matrix between all sample vectors in the high-dimensional Hilbert space based on the RBF kernel function, so as to provide basic data for solving the subsequent quadratic programming optimization problem.

[0208] In the derivation of the kernel function for the Euclidean distance in high-dimensional space, let ϕ(x) i ) is the sample x iIn a high-dimensional space, the squared Euclidean distance between two samples in that high-dimensional space can be derived using kernel function expansion: Since the RBF kernel function satisfies K(x) i ,x i )=exp(0)=1, that is Therefore, the above formula can be simplified to: This derivation avoids explicit high-dimensional computation; the high-dimensional distance can be obtained simply by using the kernel function values ​​of low-dimensional samples.

[0209] For the distance matrix construction, assuming the training set has M healthy state samples, construct an M×M Euclidean distance matrix D, where the matrix elements D i,j Let be the Euclidean distance between the i-th sample and the j-th sample in the high-dimensional space, calculated using the following formula: The distance matrix is ​​a symmetric matrix with all diagonal elements being 0 (the distance between a sample and itself is 0).

[0210] Preferably, when calculating the kernel function value, a very small numerical offset ϵ=10 can be added. -8 To avoid negative distance values ​​due to floating-point precision errors and ensure the positive definiteness of the distance matrix, the high-dimensional (M×M) sample Euclidean distance matrix D is output as the core input for solving the subsequent quadratic programming optimization problem.

[0211] Step S6.3 is a quadratic programming optimization problem to solve the Euclidean distance matrix, in order to determine the center and radius of the hypersphere and minimize the hypersphere volume containing all weighted aggregated temporal feature vectors.

[0212] Preferably, this step is the core modeling step of the SVDD algorithm. By solving a constrained quadratic programming optimization problem, the center and radius of the hypersphere in the high-dimensional space are determined, so as to minimize the volume of the hypersphere while including all healthy state samples.

[0213] In the original definition of the SVDD optimization problem, the goal of SVDD is to find a hypersphere with center a and radius R such that the mapping vector ϕ(x) of all healthy state samples is equal to the target vector ϕ(x). i All elements are located inside or on the surface of a hypersphere, with the hypersphere having the smallest possible volume. The original optimization problem is: .

[0214] The constraints are: In the formula, C is the regularization parameter, controlling the compactness and fault tolerance of the hypersphere, and ξ... i As a slack variable, it allows a small number of outlier samples (such as data acquisition noise) to exceed the hypersphere, thus avoiding model overfitting.

[0215] Among them, for the regularization parameter calibration, considering the low-noise characteristics of the health data of the GDX500 wrapping machine, the regularization parameter C is calibrated to 0.1. This value can ensure that the hypersphere tightly wraps the healthy state samples and has a certain tolerance for a small amount of acquisition noise. If the on-site data noise is large, the value of C can be appropriately increased.

[0216] Among them, for the conversion of the dual optimization problem, to simplify the solution, the original optimization problem is converted into a dual problem, and the Lagrange multipliers α i and β i are introduced through the Lagrange multiplier method, and finally the dual optimization problem is obtained: . The constraint conditions are: 0 ≤ α i ≤ C, . The dual problem only involves the calculation of the kernel function and does not require explicit high-dimensional mapping, which is more suitable for engineering implementation.

[0217] Among them, for the quadratic programming solution algorithm, the Sequential Minimal Optimization (SMO) algorithm is used to solve the dual problem. This algorithm decomposes the large-scale quadratic programming problem into multiple smallest sub-problems, and only optimizes two Lagrange multipliers each time, with high computational efficiency and being suitable for large-scale sample modeling in industrial sites.

[0218] Among them, for the determination of the hypersphere parameters: (1) Center a: The center in the high-dimensional space is the linear combination of the mapped vectors of the support vectors. The support vectors are the samples with Lagrange multiplier α i > 0. The formula for the center is: .

[0219] (2) Radius R: Select any support vector (usually select the sample with 0 <α i < C, that is, the sample on the hypersphere), and calculate the distance from its mapped vector to the center of the sphere, which is the radius of the hypersphere. The formula is: . In the formula, x s is the support vector sample.

[0220] Preferably, the set of Lagrange multipliers {α i}, the set of support vectors {x s }, and the radius R of the hypersphere calculated are used as the core parameters for subsequent sample distance calculation and health score.

[0221] Step S6.4, calculate the distance from each weighted aggregated time-series feature vector to the center of the hypersphere respectively to obtain the sample distance value.

[0222] Preferably, the core of this step is to calculate the distance from each weighted aggregated temporal feature vector (including training samples and new test samples) to the center of the high-dimensional hypersphere based on the determined hypersphere parameters, and obtain the sample distance value. The smaller the distance value, the better the health status of the device. For the kernel function calculation of sample distance, the squared distance from the mapping vector ϕ(x) of any sample x to be evaluated (training sample or newly acquired test sample) to the center a of the hypersphere can be indirectly calculated through the kernel function without explicit high-dimensional mapping. The formula is as follows: In the formula, x i For the training set samples, α i K represents the corresponding Lagrange multiplier, and K(·,·) is the RBF kernel function with calibrated parameters.

[0223] Specifically, for generating sample distance values, the square root of the squared distance value is taken to obtain the final sample distance value. This value characterizes the degree of deviation of the sample to be evaluated from the healthy state hypersphere.

[0224] Preferably, if d(x)≤R, it means that the sample is inside or on the surface of the healthy hypersphere, and the device is in a healthy state; if d(x)>R, it means that the sample is outside the healthy hypersphere, and the device may be abnormal or degraded.

[0225] Preferably, the above distance calculation can be performed on each of the weighted aggregated time-series feature vectors to be evaluated, generating a set of sample distance values ​​{d(x1), d(x2), ..., d(x...}}. N )}, where N is the total number of samples to be evaluated. Output the set of distance values ​​for all samples to be evaluated, which serves as the sole input for subsequent health score normalization.

[0226] Step S6.5: Normalize all sample distance values ​​according to the radius of the hypersphere to obtain a device health score set based on multiple dimensions.

[0227] The core of this step is to convert the sample distance values ​​into an intuitive 0-1 interval health score, normalize them using the hypersphere radius, and generate a multi-dimensional equipment health score set to facilitate health status monitoring and early warning in industrial settings.

[0228] For the health score, a linear normalization method can be used to map the sample distance values ​​to the 0-1 interval. A higher score indicates a better health status for the device. The formula is as follows: In the formula, Score(x) is the health score of sample x, d(x) is the sample distance value, and R is the radius of the hypersphere.

[0229] The scoring rules and thresholds are defined as follows: (1) Health status: When the score(x)≥0.8, the device is judged to be in a healthy state, the distance value is within 20% of the radius of the supersphere, and there are no obvious abnormal characteristics; (2) Warning state: When 0.5≤Score(x)<0.8, the equipment is judged to be in a warning state. The distance value is between 20%-50% of the radius of the supersphere, which shows early degradation characteristics and requires strengthened inspection. (3) Abnormal state: When Score(x) < 0.5, the equipment is judged to be in an abnormal state. The distance value exceeds 50% of the radius of the supersphere or exceeds the supersphere. The machine must be stopped immediately for inspection.

[0230] Specifically, for the multi-dimensional scoring set, independent SVDD hypersphere models are constructed for different components of the GDX500 packaging machine (such as the X main unit, CH auxiliary unit, and CT / CV auxiliary unit), generating corresponding component health scores, which are then integrated into a multi-dimensional equipment health scoring set to achieve precise health monitoring at the component level.

[0231] Preferably, for continuously running time-series data, a sliding window averaging method is used to smooth the health score. The window length is set to 5 time steps to avoid score jumps caused by fluctuations in single-step data, thereby improving the stability and reliability of the score. The output is a multi-dimensional equipment health score set, with each score corresponding to one sample to be evaluated, and the score range being 0-1, intuitively reflecting the equipment health status.

[0232] For example, this step can be implemented using the following pseudocode: # Pseudocode for SVDD device health modeling and scoring based on LIBSVM import numpy as np from sklearn.svm import OneClassSVM # OneClassSVM can implement SVDD by adjusting parameters from sklearn.preprocessing import normalize class SVDD_HealthAssessor: def __init__(self, kernel='rbf', sigma=0.5, C=0.1): """ Initialize the SVDD health assessment model, with parameters strictly consistent with the aforementioned implementation details. """ self.kernel = kernel self.gamma = 1 / (2 * sigma**2) # In sklearn, the RBF kernel parameter gamma = 1 / (2σ²) self.C = C self.model = None self.radius = None self.support_vectors = None def fit(self, healthy_features): """ Train the SVDD model to construct a healthy hypersphere Input healthy_features: Weighted aggregated features of health status output from step S5, with shape [M, hidden_dim] """ # Feature L2 normalization, consistent with the normalization in step S5.3. healthy_features = normalize(healthy_features, norm='l2', axis=1) # Initialize OneClassSVM, controlling regularization via the nu parameter (nu≈1 / C, adapting to the sklearn interface) nu = min(1 / self.C, 0.5) # nu must be in the interval (0,1] self.model = OneClassSVM( kernel=self.kernel, gamma = self.gamma, nu=nu, verbose=False ) self.model.fit(healthy_features) # Extracting Support Vectors self.support_vectors = self.model.support_vectors_ # Calculate the radius of the hypersphere: Select any support vector (located on the hypersphere) and calculate the distance. # The squared distance in higher-dimensional space = 2(1 - K(x,x_sv)), since K(x,x) = 1 sv = self.support_vectors[0:1] # Get the first support vector decision_func = self.model.decision_function(sv) # In OneClassSVM, decision_function = R² - d²(x), therefore R² = decision_func + d²(sv) # The support vector d²(sv) = R², therefore R² = decision_func + R² → This actually needs to be recalculated using a kernel function. # The radius is calculated directly using a kernel function here. K_sv_sv = np.exp(-self.gamma * np.sum((sv - sv)**2, axis=1)) K_sv_all = np.exp(-self.gamma * np.sum((self.support_vectors -sv)**2, axis=1)) alpha = self.model.dual_coef_[0] # Lagrange multiplier radius_sq = K_sv_sv - 2 * np.sum(alpha * K_sv_all) + np.sum(np.outer(alpha, alpha) * np.exp(-self.gamma * np.sum((self.support_vectors[:, None] - self.support_vectors[None, :])**2, axis=2))) self.radius = np.sqrt(np.maximum(radius_sq, 1e-8)) # Avoid negative values def predict_score(self, test_features): """ Predicted device health score Input test_features: The features to be evaluated output from step S5, with shape [N, hidden_dim] Output: Health score set, shape [N], range 0-1 """ # Feature L2 normalization test_features = normalize(test_features, norm='l2', axis=1) # Calculate the squared distance from each sample to the center of the hypersphere # The extended form of high-dimensional distance squared = 2(1 - K(x,x_sv)), calculated with the aid of decision_function. decision_func = self.model.decision_function(test_features) # OneClassSVM’s decision_function = R² - d²(x) → d²(x) = R² -decision_func dist_sq = self.radius**2 - decision_func dist_sq = np.maximum(dist_sq, 1e-8) # Avoid negative values dist = np.sqrt(dist_sq) # Normalized calculation of health score score = 1 - (dist / self.radius) score = np.maximum(score, 0) # Score is not lower than 0 # Smooth sliding window processing, window length 5 window_len = 5 if len(score) >= window_len: score = np.convolve(score, np.ones(window_len) / window_len, mode='same') return score # Model Training and Scoring Examples if __name__ == "__main__": # The hyperparameters are strictly consistent with the aforementioned implementation details. HIDDEN_DIM = 128 SIGMA = 0.5 C = 0.1 M = 1000 # Number of training samples for healthy states N = 100 # Number of test samples to be evaluated # Simulated Input: Health features and features to be evaluated output from step S5 healthy_feat = np.random.randn(M, HIDDEN_DIM) test_feat = np.random.randn(N, HIDDEN_DIM) # Initialize and train the SVDD model assessor = SVDD_HealthAssessor(sigma=SIGMA, C=C) assessor.fit(healthy_feat) print("Hypersphere radius:", assessor.radius) # Predicted Health Score health_scores = assessor.predict_score(test_feat) print("Health score set (first 5):", health_scores[:5]) Beneficial effects of steps S6.1 to S6.5: A complete equipment health assessment scheme is constructed through the systematic implementation of steps including high-dimensional feature space mapping using the support vector data description algorithm, calculation of the Euclidean distance matrix using the Gaussian radial basis function kernel, solving the quadratic programming optimization problem, and calculation and normalization of sample distance values. This scheme can construct a hypersphere model of equipment health status based on weighted aggregated temporal feature vectors. By calculating the distance from the sample to the center of the hypersphere, it achieves quantitative assessment and anomaly detection of equipment health status, accurately reflecting the current operating status and health level of the equipment, and providing a scientific basis for equipment maintenance decisions. Specifically, step S6.1 uses the support vector data description algorithm to map the weighted aggregated temporal feature vectors to a high-dimensional feature space, mapping low-dimensional linearly inseparable features to a high-dimensional linearly separable space, providing suitable feature representations for subsequent distance calculation and anomaly detection; step S6.2 calculates the Euclidean distance matrix between sample vectors in the high-dimensional feature space using the Gaussian radial basis function kernel, accurately measuring the similarity and difference between different samples in the feature space, providing basic distance information for the construction of the hypersphere model; step S6.3 solves the quadratic programming optimization problem of the Euclidean distance matrix to determine the hypersphere. The center and radius of the sphere are set to minimize the volume of the hypersphere and include all weighted aggregated temporal feature vectors to establish a decision boundary that can accurately characterize the distribution range of normal equipment status. Step S6.4 calculates the distance from each weighted aggregated temporal feature vector to the center of the hypersphere to obtain the sample distance value, quantifying the degree to which each sample deviates from the center of the normal status distribution. Step S6.5 normalizes all sample distance values ​​according to the radius of the hypersphere to obtain a set of equipment health scores based on multiple dimensions, unifying the numerical range of scores in different dimensions, and facilitating intuitive assessment and horizontal comparison of equipment health status.

[0233] In summary, the overall beneficial effects of steps S1 to S6 are as follows: By organically combining technologies such as multi-sensor data acquisition, signal decomposition, time-frequency feature extraction, time-series learning, attention mechanism weighting, and anomaly detection, a complete health monitoring solution for packaging machines is formed. This solution can monitor the operating status of packaging machines in real time, automatically extract deep fault features from multi-source heterogeneous signals, effectively assess the health level of the equipment, and achieve early fault warning. It solves the problems of long detection cycles, strong subjectivity, and low accuracy associated with traditional methods that rely on manual inspection and single-parameter threshold judgment, significantly improving the real-time performance, accuracy, and reliability of equipment health status monitoring. Specifically, step S1 denoises and normalizes vibration, temperature, current, and speed signals to eliminate noise interference and dimensional differences, providing high-quality basic data for subsequent analysis; step S2 uses an empirical mode decomposition algorithm to adaptively decompose the signal into feature components of different scales, improving the separability of signal features; and step S3 extracts instantaneous frequency, instantaneous amplitude, and marginal spectral energy density through Hilbert-Huang transform, achieving accurate characterization of the signal's time-frequency features. Step S4 uses a long short-term memory network to capture the temporal dependencies in the feature sequence and obtain a feature vector containing deep temporal information; Step S5 calculates weight coefficients and aggregates key features through an attention mechanism to improve the effectiveness of feature representation; Step S6 constructs a hypersphere model of the device health status based on the support vector data description algorithm and calculates the health score to achieve quantitative assessment and anomaly detection of the device health status.

[0234] like Figure 2 As shown, this embodiment provides an example of a health detection system for packaging machine equipment. In this embodiment, the health detection system is applied to the health detection method as described in the above embodiment.

[0235] Specifically, the health detection system includes a raw signal processing module 1, a raw signal decomposition module 2, a signal component transformation module 3, a signal component feature learning module 4, a time-series feature vector weighted aggregation module 5, and a device health score calculation module 6, which are connected electrically or through communication in sequence.

[0236] The system comprises the following modules: Raw signal processing module 1 collects vibration, temperature, current, and rotation speed signals from the packaging machine in real time via external sensors, and performs denoising and normalization on the collected raw signals; Raw signal decomposition module 2 inputs the preprocessed vibration, temperature, current, and rotation speed signals into an empirical mode decomposition algorithm to decompose them into several intrinsic mode function components and several residual components; Signal component transformation module 3 performs Hilbert-Huang transform on each intrinsic mode function component and each residual component to calculate time-frequency features, extracting instantaneous frequency, instantaneous amplitude, and marginal spectral energy density; Signal component feature learning module 4 learns all instantaneous frequencies and all... Instantaneous amplitude and energy density of all marginal spectra are input into a long short-term memory network for temporal feature learning. The cell state and hidden state are updated through a gating mechanism of forget gate, input gate, and output gate to obtain a temporal feature vector. The temporal feature vector weighted aggregation module 5 is used to input the temporal feature vector into an attention mechanism to calculate the weight coefficient of each time step, and calculate the weighted aggregated temporal feature vector by weighted summation and aggregation of feature information of key time steps. The device health score calculation module 6 is used to input the weighted aggregated temporal feature vector into a support vector data description algorithm to construct a hypersphere model of device health status. The distance from the sample to the center of the hypersphere is calculated by mapping to a high-dimensional feature space through a kernel function to obtain the device health score.

[0237] Figure 3 This is a schematic diagram of the structure of an electronic device according to an embodiment of this application. Figure 3 As shown, the electronic device 7 includes a processor 71 and a memory 72 coupled to the processor 71.

[0238] The memory 72 stores program instructions for implementing the health detection method of the packaging machine equipment in any of the above embodiments.

[0239] The processor 71 is used to execute program instructions stored in the memory 72 to perform health checks on the packaging machine equipment.

[0240] The processor 71 can also be referred to as a CPU (Central Processing Unit). The processor 71 may be an integrated circuit chip with signal processing capabilities. The processor 71 can also be a general-purpose processor, a digital signal processor (DSP), an application-specific integrated circuit (ASIC), a field-programmable gate array (FPGA), or other programmable logic devices, discrete gate or transistor logic devices, or discrete hardware components. A general-purpose processor can be a microprocessor or any conventional processor.

[0241] Furthermore, Figure 4 This is a schematic diagram of the structure of a storage medium according to an embodiment of this application. See also: Figure 4 The storage medium 8 in this embodiment stores program instructions 81 capable of implementing all the above methods. These program instructions 81 can be stored in the storage medium as a software product, including several instructions to cause a computer device (which may be a personal computer, server, or network device, etc.) or processor to execute all or part of the steps of the methods in each embodiment of this application. The aforementioned storage medium includes various media capable of storing program code, such as USB flash drives, portable hard drives, read-only memory (ROM), random access memory (RAM), magnetic disks, or optical disks, or terminal devices such as computers, servers, mobile phones, and tablets.

[0242] In the several embodiments provided in this application, it should be understood that the disclosed systems, methods, and approaches can be implemented in other ways. For example, the system embodiments described above are merely illustrative; for instance, the division of units is only a logical functional division, and in actual implementation, there may be other division methods. For example, multiple units or components may be combined or integrated into another system, or some features may be ignored or not executed. Furthermore, the mutual coupling or direct coupling or communication connection shown or discussed may be through some interfaces; the indirect coupling or communication connection between systems or units may be electrical, mechanical, signal, or other forms.

[0243] Furthermore, the functional units in the various embodiments of this application can be integrated into one processing unit, or each unit can exist physically separately, or two or more units can be integrated into one unit. The integrated units described above can be implemented in hardware or as software functional units. The above are merely embodiments of this application and do not limit the patent scope of this application. Any equivalent structural or procedural transformations made based on the description and drawings of this application, or direct or indirect applications in other related technical fields, are similarly included within the patent protection scope of this application.

Claims

1. A method of detecting the health of a packaging machine apparatus, characterized by, The health status detection method includes: Step S1: The vibration signal, temperature signal, current signal, and speed signal of the packaging machine are collected in real time through external sensors, and the collected raw signals are denoised and normalized. Step S2: Input the preprocessed vibration signal, preprocessed temperature signal, preprocessed current signal, and preprocessed rotation speed signal into the empirical mode decomposition algorithm to decompose them into several intrinsic mode function components and several residual components. Step S3: Perform Hilbert-Huang transform on each intrinsic mode function component and each residual component to calculate the time-frequency characteristics, so as to extract the instantaneous frequency, instantaneous amplitude, and marginal spectral energy density. Step S4: Input all instantaneous frequencies, all instantaneous amplitudes, and all marginal spectral energy densities into the Long Short-Term Memory network for temporal feature learning. Update the cell state and hidden state through the gating mechanism of forget gate, input gate, and output gate to obtain the temporal feature vector. Step S5: Input the temporal feature vector into the attention mechanism to calculate the weight coefficient of each time step, and calculate the weighted aggregated temporal feature vector by weighted summation and aggregating the feature information of key time steps; Step S6: Input the weighted aggregated temporal feature vector into the support vector data description algorithm to construct a hypersphere model of the device health status. Calculate the distance from the sample to the center of the hypersphere by mapping to a high-dimensional feature space through a kernel function to obtain the device health score.

2. The healthiness detection method according to claim 1, characterized in that, Step S1 involves real-time acquisition of vibration, temperature, current, and rotation speed signals from the packaging machine via external sensors, followed by denoising and normalization of the acquired raw signals, including: Step S1.1: Acquire vibration acceleration signal through accelerometer, temperature signal through thermocouple sensor, current signal through Hall sensor, and rotation speed signal through encoder to obtain original vibration signal, original temperature signal, original current signal, and original rotation speed signal; Step S1.2: Input the original vibration signal, the original temperature signal, the original current signal and the original rotation speed signal into the wavelet threshold denoising algorithm respectively, and decompose and reconstruct them through the sym4 wavelet basis function to obtain the denoised vibration signal, the denoised temperature signal, the denoised current signal and the denoised rotation speed signal. Step S1.3: The denoised vibration signal, the denoised temperature signal, the denoised current signal, and the denoised rotation speed signal are normalized by Z-score standardization to obtain normalized vibration signal, normalized temperature signal, normalized current signal, and normalized rotation speed signal, respectively. Step S1.4: Align the normalized vibration signal, the normalized temperature signal, the normalized current signal, and the normalized rotational speed signal according to the timestamp, and fill in the missing data points using a linear interpolation method; Step S1.5: After filling in the missing data points, the signal is filtered by a moving average algorithm to eliminate high-frequency jitter interference, resulting in preprocessed vibration signal, preprocessed temperature signal, preprocessed current signal, and preprocessed rotation speed signal.

3. The healthiness detection method according to claim 1, characterized in that, Step S2: Input the preprocessed vibration signal, preprocessed temperature signal, preprocessed current signal, and preprocessed rotational speed signal into the empirical mode decomposition algorithm to decompose them into several intrinsic mode function components and several residual components, including: Step S2.1: Input the preprocessed vibration signal, preprocessed temperature signal, preprocessed current signal, and preprocessed rotation speed signal into the empirical mode decomposition algorithm to identify all local maxima and local minima in the signal. Step S2.2: Based on the same signal, the upper envelope is obtained by fitting the local maxima points through cubic spline interpolation, and the lower envelope is obtained by fitting the local minima points. Step S2.3: Calculate the mean envelope of the upper envelope and the lower envelope, and subtract the mean envelope from the original signal of the same signal to obtain the candidate component; Step S2.4: Determine whether the candidate component satisfies the intrinsic mode function condition. If not, use the candidate component as a new input and continue the screening iteration through steps S2.2 to S2.3 until the intrinsic mode function condition is satisfied. Step S2.5: Extract the candidate components that satisfy the intrinsic mode function conditions as intrinsic mode function components, and subtract the intrinsic mode function components from the original signal of the same signal to obtain the remaining signal; Step S2.6: The remaining signal is used as a new input and the decomposition process is repeated from step S2.2 to step S2.5 until the remaining signal changes monotonically and cannot be decomposed further, resulting in several intrinsic mode function components and several residual components.

4. The healthiness detection method according to claim 1, characterized in that, Step S3 involves performing Hilbert-Huang transform on each intrinsic mode function component and each residual component to calculate time-frequency features, thereby extracting instantaneous frequency, instantaneous amplitude, and marginal spectral energy density, including: Step S3.1: Perform Hilbert transform on each intrinsic mode function component and each residual component, and calculate the real and imaginary parts of the analytic signal through convolution operation; Step S3.2: Calculate the amplitude of the analytic signal using the real part and the imaginary part, and obtain the instantaneous frequency by taking the time derivative of the amplitude of the analytic signal. Step S3.3: Calculate the magnitude of the analytic signal using the real part and the imaginary part to obtain the instantaneous amplitude; Step S3.4: Integrate the instantaneous frequency and the instantaneous amplitude over time to obtain the marginal spectral energy density.

5. The health status detection method according to claim 1, characterized in that, Step S4: Input all instantaneous frequencies, all instantaneous amplitudes, and all marginal spectral energy densities into the Long Short-Term Memory (LSTM) network for temporal feature learning. Update the cell state and hidden state through gating mechanisms of forget gate, input gate, and output gate to obtain the temporal feature vector, including: Step S4.1: Concatenate all instantaneous frequencies, all instantaneous amplitudes, and all marginal spectral energy densities in chronological order to form a feature sequence; Step S4.2: Input the feature sequence into the Long Short-Term Memory (LSTM) network, and calculate the forgetting gate coefficient through the forgetting gate of the LSM network to control the degree of retention of the cell state at the previous time step; Step S4.3: Calculate the input gating coefficient and candidate cell state through the input gate of the Long Short-Term Memory network to determine the new information at the current moment; Step S4.4: Update the current cell state using the forgetting gating coefficient and the input gating coefficient; Step S4.5: Calculate the output gating coefficient through the output gate of the Long Short-Term Memory Network to determine the hidden state at the current moment based on the updated cell state; Step S4.6: Repeat steps S4.2 to S4.5 until all time steps have been processed to obtain a time-series feature vector with temporal dependencies.

6. The health status detection method according to claim 1, characterized in that, Step S5 involves inputting the temporal feature vector into an attention mechanism to calculate the weight coefficient for each time step, and then calculating the weighted aggregated temporal feature vector by weighted summation and aggregating the feature information of key time steps, including: Step S5.1: Input the temporal feature vector into the attention mechanism and calculate the relevance score between the query vector and the key vector through dot product operation; Step S5.2: Perform an exponential operation on the correlation score and divide it by the sum of the exponential operation results of all time steps to obtain the weight coefficient of each time step. Step S5.3: The feature information of each time step in the time series feature vector is weighted and summed using the weight coefficients to obtain the weighted aggregated time series feature vector.

7. The health status detection method according to claim 1, characterized in that, Step S6: Input the weighted aggregated temporal feature vector into the support vector data description algorithm to construct a hypersphere model of the device health status. Calculate the distance from the sample to the center of the hypersphere by mapping to a high-dimensional feature space using a kernel function, and obtain the device health score, including: Step S6.1: Map the weighted aggregated temporal feature vectors to a high-dimensional feature space using the support vector data description algorithm; Step S6.2: Calculate the Euclidean distance matrix between sample vectors in the high-dimensional feature space using the Gaussian radial basis kernel function; Step S6.3: Solve the quadratic programming optimization problem of the Euclidean distance matrix to determine the center and radius of the hypersphere in order to minimize the hypersphere volume containing all weighted aggregated temporal feature vectors; Step S6.4: Calculate the distance from each weighted aggregated temporal feature vector to the center of the hypersphere to obtain the sample distance value; Step S6.5: Normalize all sample distance values ​​according to the radius of the hypersphere to obtain a device health score set based on multiple dimensions.

8. A health detection system for packaging machine equipment, wherein the health detection system is applied to the health detection method as described in any one of claims 1 to 7, characterized in that, The health detection system includes: The raw signal processing module is used to collect vibration signals, temperature signals, current signals, and speed signals of the packaging machine in real time through external sensors, and to perform noise reduction and normalization processing on the collected raw signals. The original signal decomposition module is used to input the preprocessed vibration signal, preprocessed temperature signal, preprocessed current signal, and preprocessed rotation speed signal into the empirical mode decomposition algorithm to decompose them into several intrinsic mode function components and several residual components. The signal component transformation module is used to perform Hilbert-Huang transform on each intrinsic mode function component and each residual component to calculate the time-frequency characteristics, so as to extract the instantaneous frequency, instantaneous amplitude, and marginal spectral energy density. The signal component feature learning module is used to input all instantaneous frequencies, all instantaneous amplitudes, and all marginal spectral energy densities into the Long Short-Term Memory network for temporal feature learning. The cell state and hidden state are updated through the gating mechanism of forget gate, input gate, and output gate to obtain the temporal feature vector. The temporal feature vector weighted aggregation module is used to input the temporal feature vector into the attention mechanism to calculate the weight coefficient of each time step, and to calculate the weighted aggregated temporal feature vector by weighted summation and aggregation of the feature information of key time steps; The device health score calculation module is used to input the weighted aggregated time-series feature vector into the support vector data description algorithm to construct a hypersphere model of the device health status, and calculate the distance from the sample to the center of the hypersphere by mapping to a high-dimensional feature space through a kernel function to obtain the device health score.

9. An electronic device, characterized in that, The method includes a processor and a memory coupled to the processor, the memory storing program instructions executable by the processor; when the processor executes the program instructions stored in the memory, it implements the health detection method as described in any one of claims 1 to 7.

10. A computer-readable storage medium, characterized in that, The computer-readable storage medium stores program instructions that, when executed by a processor, enable the health detection method as described in any one of claims 1 to 7.

Images