Equipment state monitoring method based on vibration signals

Abstract

An equipment state monitoring method based on a vibration signal comprises an initial sample acquisition stage and a monitoring stage, and the initial sample acquisition stage comprises the following steps: placing a three-axis acceleration sensor and a tested instrument and setting three-axis directions; collecting an acceleration data sample; preprocessing and de-noising; establishing a naive Bayes classifier to calculate initial probability distribution; the monitoring stage comprises real-time data acquisition and acquired data analysis, and specifically comprises the following steps: preprocessing and denoising the acquired data; and analyzing the current state of the data output instrument. According to the invention, whether the instrument is actually operated on site manually or is in a standby running state can be distinguished accurately, objectivity of monitoring the use condition of the instrument is improved, the operation state of the instrument can be further distinguished when the instrument is operated manually, and accuracy of monitoring the state of the instrument is improved.

Description

Technical Field

[0001] The present invention belongs to the field of instrument and equipment monitoring, and relates to classifying and judging the use status of instruments and equipment, and is specifically a method for monitoring equipment status based on vibration signals. Background Art

[0002] Research institutes, schools, hospitals, and other institutions are equipped with numerous instruments and equipment. As a key component of fixed assets, these instruments and equipment typically require scientific planning and management to fully tap their potential and maximize their value. To achieve this goal, managers need to monitor instrument usage in real time to rationally allocate and optimize instrument resources. By promoting resource sharing of instruments and equipment, resource utilization can be further improved, duplicate purchases can be avoided, and stronger support can be provided for scientific research, teaching, and medical activities.

[0003] Traditional equipment management relies primarily on manual operations, a method characterized by cumbersome procedures, low efficiency, and prone to errors. Manually recording equipment usage, operation times, and maintenance status is not only time-consuming and labor-intensive, but also prone to inaccurate data due to oversight or recording errors, compromising the scientific nature of management decisions. Furthermore, manual management makes it difficult to monitor equipment status in real time and respond quickly to equipment idleness, resulting in low resource utilization and high management costs. In recent years, many new instrument management methods have emerged. For example, the patent "Instrument and Equipment Usage Recording Device" (Publication No.: CN204631962U) incorporates an external card swiping device that allows instrument use and records user information and usage time. However, the recorded time does not represent the actual time the instrument was used. The patent "An Intelligent Data Management System and Method for Laboratory Instruments Based on RFID Electronic Tags" (Publication No.: CN119398355A) uses an auxiliary reader to continuously scan and read RFID tags on laboratory instruments, including information such as status, location, and usage history. However, the information read requires the user to enter it manually, which is cumbersome and does not represent the actual usage of the instrument. The patent "A big data analysis method based on instrument current data" (publication number: CN110414768B) determines whether someone is using the instrument by monitoring the current status during operation, but the fact that the instrument is powered on does not completely mean that the instrument is being used manually.

[0004] In practical applications, existing methods often suffer from several issues, resulting in ineffective monitoring of instrument usage. Furthermore, the accuracy of instrument usage data is poor, making it difficult to accurately assess the actual status of instrument usage in real time. These issues hinder the precise assessment of instrument asset utilization, impacting the scientific and effective nature of resource management. Summary of the Invention

[0005] In response to the shortcomings of the above-mentioned prior art, the purpose of the present invention is to address the shortcomings of existing instrument and equipment management methods, especially the problems of cumbersome supervision procedures, poor real-time performance, and low accuracy. A vibration signal-based instrument usage classification method is proposed. Based on the instrument's usage status, it is distinguished whether the instrument is being used. If it is being used, it is further determined in which mode the instrument is being used, such as whether a button is pressed or a knob is rotated. The frequency of device usage and the number of times each button is used are counted at the device terminal.

[0006] The device status monitoring method based on vibration signals of the present invention includes an initial sample collection phase and a monitoring phase. The initial sample collection phase includes the following steps:

[0007] A-1. Place the triaxial accelerometer, device management terminal, and at least one device under test in the same positions as they will be during normal use, allowing the triaxial accelerometer to detect the motion of the device under test. Set the triaxial directions of the triaxial accelerometer.

[0008] A-2. Collect X, Y and Z axis acceleration data samples of the three-axis acceleration sensor at the same sampling frequency. The acceleration data samples refer to the acceleration data samples of the tested instrument in different working states C i The acceleration values ​​of the X, Y and Z axes under each working state are collected multiple times and the average is taken;

[0009] A-3. The device management terminal performs preprocessing and denoising on the acceleration data collected in step A-2 to obtain a denoised signal sample;

[0010] A-4. Build a naive Bayes classifier to calculate the initial probability distribution, specifically:

[0011] A-4-1. Determine the classification category C = (C1, C2…C H ), where C i Indicates the categories corresponding to different working states, and H is the total number of states;

[0012] A-4-2. Slice each denoised signal sample obtained in step A-3, where the slicing time is greater than the time required for the basic operation of the tested instrument;

[0013] A-4-3. Calculate the feature vectors in each slice time period, including: X i =(X1: X-axis peak value, X2: Y-axis peak value, X3: Z-axis peak value, X4: X-axis variance, X5: Y-axis variance, X6: Z-axis variance, ..., X 3m-2 : Number of X-axis peaks, X 3m-1 : Number of peaks on the Y axis, X 3m: number of Z-axis spectrum peaks), i = 1, 2…3m, m is the total number of eigenvector types;

[0014] A-4-4. Establish a discretization standard, divide different vibration signal features into different intervals, normalize the vibration signal features in the same interval to the same eigenvalue, and obtain the normalized eigenvector X i

[0015] A-4-5. Use the maximum likelihood estimation method to estimate the normalized eigenvector X i Calculate the conditional probability P(X i |C n ),Right now:

[0016] P(X1|C1),P(X2|C1),…P(X 3m |C1);P(X1|C2),P(X2|C2),…P(X 3m |C n );

[0017] ...; P(X1|C2), P(X2|C n ),…P(X 3m |C H );

[0018] A-4-6. Initially, the probability of each category appearing in each category is equal, that is, the probability of each category appearing

[0019] P(C1)=P(C2)=…=P(C H )=1 / H;

[0020] The monitoring phase includes:

[0021] D. Real-time data collection: Use the same placement method and sampling frequency as step A-1 and step A-2 to collect data and upload it in real time;

[0022] E. Data collection and analysis, including:

[0023] C-1. The collected data is preprocessed and denoised in the same manner as in step A-3;

[0024] C-3. Analyze the number of times each usage scenario occurs in the data, specifically:

[0025] C-3-1. Perform the same slicing process on the data set as in step A-4-2, calculate all eigenvalues ​​for each slice, and output the normalized eigenvector y according to the discretization standard established in step A-4-4;

[0026] C-3-2. Based on the normalized feature vector y calculated in step C-3-1, and the conditional probability established in step A-4-5 and the prior probability in step A-4-6, calculate the probability of the data set in different classification categories C n The first posterior probability in :

[0027]

[0028] C-3-3. Find the category with the highest probability among all categories, that is

[0029] P(C k |y)=max{P(C1|y), P(C2|y),…,P(C n |y)}, then this set of data is considered to belong to category C k ; That is, the output instrument is currently in state C k middle.

[0030] Preferably, after step C-1 and before step C-3, step C-2 is further included. The displacement MX of the processed signal is calculated. If it is found that the displacement MA is greater than the set movement threshold, the subsequent steps are terminated, otherwise the subsequent steps are performed.

[0031] Preferably, the A-3 step specifically includes:

[0032] A-3-1. De-compensate the acceleration data by first calculating the average acceleration of multiple samples taken along the X, Y, and Z axes when the instrument is stationary. Then, subtract the average of the corresponding axis from all measured acceleration data while stationary. This allows the corrected acceleration value to fluctuate around zero, eliminating the effects of inherent sensor bias.

[0033] A-3-2. Perform EMD decomposition on the acceleration data samples to obtain several intrinsic modal components imf i , i=1,2…n-1. i represents the serial number of different frequency components, n is the total number of frequency components, and the acceleration signal x(t) is expressed as:

[0034] x(t)=imf1(t)+imf2(t)+...+imf n-1 (t)+r n (t);

[0035] Among them, imf i (t) represents the intrinsic modal components with different frequency components, arranged from high frequency to low frequency, r n (t) is the residual component, t represents time;

[0036] A-3-3. Calculate adjacent IMF i (t) mutual information entropy I(imfi , imf i+1 ), first calculate each imf i Energy entropy of (t)

[0037] The total energy E i is the energy entropy of the i-th component, P i is the energy of the i-th component. Then the inherent energy entropy H(IMF) corresponding to each connotation modal component is i )for:

[0038]

[0039] The mutual information entropy of adjacent components is:

[0040] I(imf i , imf i+1 )=H(IMF i )+H(IMF i+1 )-H(IMF i , IMF i+1 )

[0041] in

[0042]

[0043] P(x,y) represents the joint probability distribution of IMFi and IMFi+1;

[0044] A-3-4. Find I(imf i , imf i+1 ) is the first minimum turning point, the minimum turning point is defined as: if I(imf i , imf i+1 ) decreases and I(imf i , imf i+2 ) increases, then the i value at this time is the minimum turning point, and the first minimum value that appears is the high-low frequency dividing point. Let the minimum turning point be k, and rewrite x(t) with k as the dividing point as follows:

[0045] x(t)=h i (t)+I j (t)+r n (t)

[0046] where h i (t) is the high frequency signal, i=1,2…k,

[0047] I j (t) is the low-frequency signal, j=k+1, k+2…n-1.

[0048] A-3-5. For high frequency signal h i (t) is subjected to adaptive threshold denoising to obtain the denoised high-frequency signal hp(t); the noise threshold is calculated as follows: Where N is the signal length, σ is the standard deviation of the noise component, and x i is the high frequency IMF component, that is, x i =h i (t), median represents the median function;

[0049] Calculate h separately i (t) The corresponding noise threshold T i , h i (t) is below the noise threshold T i The value of is set to 0 to obtain the high-frequency signal hp(t) after noise reduction;

[0050] A-3-6. Add the denoised high-frequency signal to the low-frequency signal and discard the residual component r n Get the denoised signal sample:

[0051] xp(t)=hp(t)+I j (t).

[0052] Preferably, in step A-3, before step A-3-2 begins, step A-3-1 is also included. The acceleration data samples are de-compensated. The specific method is: first calculate the average acceleration of m samples of the X, Y, and Z axes when the instrument is stationary, then subtract the average of the corresponding axis from all the measured acceleration data when it is stationary, and then proceed to step A-3-2.

[0053] Preferably, step A-3-5 is specifically as follows:

[0054] Calculating the noise threshold Where N is the signal length, σ is the standard deviation of the noise component, and x i is the high frequency IMF component, that is, x i =h i (t), median represents the median function;

[0055] Calculate h separately i (t) The corresponding noise threshold T i , h i (t) is below the noise threshold T i The value of is set to 0 to obtain the high-frequency signal hp(t) after noise reduction.

[0056] Preferably, step C further comprises:

[0057] C-4. Count the number of occurrences of each state and the total number of times the instrument is used (MA) to obtain the second posterior probability of different states;

[0058] C-5. Continuously update the posterior probability during the detection process, specifically:

[0059] The first posterior probability is used to replace the conditional probability obtained in step A-4-5, and the second posterior probability is used to replace the occurrence probability of each category in step A-4-6, and the updated data is used for detection in the subsequent detection process.

[0060] The present invention is based on the three-dimensional vibration signals of the instrument when it is in use. By establishing a feature set to classify the signals, it has the following beneficial effects:

[0061] 1. The present invention can accurately distinguish whether the instrument is currently being actually operated on-site or in a standby state, thereby increasing the objectivity of monitoring the usage status of the instrument.

[0062] 2. When the instrument is manually operated, it can further distinguish which operating state the instrument is in, such as key operation or knob operation, thereby improving the accuracy of instrument status monitoring.

[0063] 3. The method of the present invention is applicable to instruments and equipment with an operation panel, has strong universality, and can be widely used in monitoring the usage status of different instruments and equipment.

[0064] 4. The present invention can automatically realize equipment status monitoring with only one acceleration sensor module at the acquisition end, effectively reducing the cost of equipment status monitoring.

[0065] 5. The usage status of instruments and equipment obtained by the present invention is more real and objective, and can provide reliable data support for optimizing instrument resource allocation and improving asset utilization. BRIEF DESCRIPTION OF THE DRAWINGS

[0066] Figure 1 This is a schematic diagram of a specific arrangement of the three-axis acceleration sensor module, the device to be detected, and the device management terminal in step A of the present invention.

[0067] Figure 2 This is a flow chart of a specific implementation of step A of the present invention;

[0068] Figure 3 This is a schematic flow chart of a specific implementation of step C of the present invention. DETAILED DESCRIPTION

[0069] To make the objectives, technical solutions, and advantages of the present invention more clear, the technical solutions of the present invention will be clearly and completely explained below in conjunction with the specific embodiments of the present invention and the corresponding drawings. Obviously, the embodiments described are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts are within the scope of protection of the present invention.

[0070] The present invention will be further described below in conjunction with the accompanying drawings and specific preferred embodiments, but the scope of protection of the present invention is not limited thereby.

[0071] The present invention is implemented in an instrument management system based on two principles:

[0072] 1. Extract features from the collected acceleration data and establish a classifier.

[0073] 2. Match the actual collected data with the instrument usage.

[0074] Before the detection of the present invention, initial sample collection is required. The initial sample collection stage is completed by the following steps:

[0075] A-1. The layout diagram of each module of the entire system is as follows Figure 1 As shown, the placement position should be exactly the same as the position during subsequent normal use, so that the spatial relationship of the instruments during subsequent normal use is consistent with that at this time.

[0076] Each module includes a three-axis acceleration sensor S2, a device management terminal S4, and at least one tested instrument S1. The tested instrument may include an oscilloscope, a spectrum analyzer, and other commonly used laboratory instruments with button and knob functions. Figure 1 In the figure, S1 is the instrument to be tested, which includes a knob and a button S3. When in use, the three-axis acceleration sensor can be installed on the upper side of the casing of the instrument to be tested, usually with the positive direction of the X-axis parallel to the operation panel, for collecting acceleration information of the X, Y and Z axes. On the right is the equipment management terminal, which is specifically a fixed computer device, including components such as a host, a display, a keyboard and a mouse. The host is equipped with core hardware such as a central processing unit (CPU), memory, hard disk and power supply for data processing and storage. The display is used to display the operation interface and operation results. The desktop computer is wirelessly connected to the instrument management terminal for receiving, storing and processing the data collected by the sensor, and performing functions such as data analysis, status classification and resource management. The wireless communication methods include WIFI, Bluetooth, LoRa and other methods.

[0077] A-2. Collect the X, Y and Z axis acceleration data samples of the three-axis acceleration sensor respectively and upload them to the instrument management terminal. The acceleration data samples refer to the X, Y and Z axis acceleration values ​​of the instrument under different working states. The different working states here include the vibration state of the instrument when it is used manually during normal use and the vibration state when it is not used manually. The vibration state of the instrument when it is used manually can be further divided into categories such as key operation and knob operation according to the operation method, depending on the design and use of the instrument; the vibration state when it is not used manually is more complicated, and may include categories such as static, desktop tapping, desktop shaking, etc., depending on the environment in which the instrument is located.

[0078] When collecting acceleration sample data in step A-2, each signal type must be sampled at the same sampling frequency w. The selection of the sampling frequency must comprehensively consider the timing constraints of button presses and rebounds, as well as the vibration characteristics of rigid materials such as the paddle spring, to ensure accurate capture of the dynamic changes in the signal. Therefore, the sampling frequency should meet the minimum requirement of w > 33 Hz to ensure data accuracy and reliability.

[0079] When collecting data in step A-2, the single variable principle should be adhered to, meaning that a set of data contains only one operating mode. To more closely resemble actual conditions, the operating modes of the same operable device on the panel of the instrument being tested can be used once, and the resulting data set can then be uploaded to the backend management terminal in real time.

[0080] In step A-2, the average acceleration values ​​of the X, Y, and Z axes of all the measured instruments are usually collected first when they are stationary. Then, the measured instruments are grouped and measured according to their different working modes.

[0081] A-3. The device management terminal performs pre-processing and denoising on the collected acceleration data. The steps involved are as follows: Figure 2 As shown, the following is a detailed explanation of the process:

[0082] A-3-1. De-compensate the acceleration data samples by first calculating the average acceleration of 100 samples of the X, Y, and Z axes when the instrument is stationary. Then, subtract the average value of the corresponding axis from all measured acceleration data when the instrument is stationary. This allows the corrected acceleration value to fluctuate around zero, eliminating the influence of the sensor's inherent bias.

[0083] A-3-2. Perform EMD (Empirical Mode Decomposition) on the acceleration data sample to obtain several intrinsic modal components imf i , i=1,2…n-1. i represents the sequence number of different frequency components, n is the total number of frequency components, and the acceleration signal x(t) can be expressed as:

[0084] x(t)=imf1(t)+imf2(t)+...+imf n-1 (t)+r n (t)

[0085] Among them, imf i (t) represents the intrinsic modal components with different frequency components, arranged from high frequency to low frequency, r n (t) is the residual component, and t represents time.

[0086] A-3-3. Calculate adjacent IMF i (t) mutual information entropy I(imf i , imf i+1 ), first calculate each imf i Energy entropy of (t)

[0087] The total energy E i is the energy entropy of the i-th component, P i is the energy of the i-th component. Then the inherent energy entropy H(IMF) corresponding to each connotation modal component is i )for:

[0088]

[0089] At this time, the mutual information entropy of adjacent components can be calculated as:

[0090] I(imf i , imf i+1 )=H(IMF i )+H(IMF i+1 )-H(IMF i , IMF i+1 )

[0091] in

[0092]

[0093] P(x,y) represents the joint probability distribution of IMFi and IMFi+1

[0094] A-3-4. The mutual information entropy of adjacent components will show the characteristics of first decreasing and then increasing. Find I(imf i , imf i+1 ) is the first minimum turning point, the minimum turning point is defined as: if I(imf i , imf i+1 ) decreases and I(imf i , imf i+2) increases, then the i value at this time is the minimum turning point, and the first minimum value that appears is the high-frequency and low-frequency dividing point of the IMF. Let the minimum turning point be k, and rewrite x(t) with k as the dividing point as:

[0095] x(t)=h i (t)+I j (t)+r n (t)

[0096] where h i (t) is the high frequency signal, i=1,2…k,

[0097] I j (t) is the low-frequency signal, j=k+1, k+2…n-1.

[0098] A-3-5. For high frequency signal h i (t) is subjected to adaptive threshold noise reduction to obtain the high-frequency signal hp(t) after noise reduction. The noise threshold is calculated as follows: Where N is the signal length, σ is the standard deviation of the noise component, and x i is the high frequency IMF component, that is, x i =h i (t), median represents the median function;

[0099] Calculate h separately i (t) The corresponding noise threshold T i , h i (t) is below the noise threshold T i The value of is set to 0 to obtain the high-frequency signal hp(t) after noise reduction.

[0100] A-3-6. Add the denoised high-frequency signal to the low-frequency signal and discard the residual component r n Get the denoised signal sample:

[0101] xp(t)=hp(t)+I j (t)

[0102] A-4. Establish a naive Bayes classifier to calculate the initial probability distribution. The steps involved are as follows Figure 3 As shown, the following is a detailed explanation of the process:

[0103] A-4-1. Determine the classification category C = (C1, C2…C n ), such as C1: Instrument still, C2: Instrument pressed, C3: Instrument knob, C4: Desktop knock, C5: Desktop shaking, etc.

[0104] A-4-2. Slice each denoised signal sample obtained in step A-3-6. Set the slicing time t1 to a value generally ranging from 0.8 to 1.5 seconds, which is slightly longer than the time required for a basic operation of the instrument under test, such as normal key presses.

[0105] A-4-3. Determine the eigenvector. In the present invention, the types of eigenvectors include but are not limited to peak value, variance, kurtosis and skewness in the time domain, as well as the main peak frequency, main peak intensity and number of spectral peaks in the frequency domain after calculating its power spectrum. The above eigenvectors need to extract the features of the X, Y and Z axes respectively, that is: X i =(X1: X-axis peak value, X2: Y-axis peak value, X3: Z-axis peak value, X4: X-axis variance, X5: Y-axis variance, X6: Z-axis variance, ..., X 3m-2 : Number of X-axis peaks, X 3m-1 : Number of peaks on the Y axis, X 3m : number of Z-axis spectral peaks), i = 1, 2…3m. Calculate the values ​​of the above features in the slice time period to obtain the feature vector.

[0106] A-4-4. Establish a discretization standard. Since the collected vibration signal fluctuates within a range of values, it is necessary to discretize this range. Different vibration signal characteristics are assigned to different ranges, and vibration signal characteristics within the same range are normalized to the same characteristic value.

[0107] For example, set the interval division criteria: 1: 0-50, 2: 50-250, 3: 250-450, 4: 450-650. This means that when the vibration signal feature value is between 0 and 50, the feature value is normalized to 1, and when the vibration signal feature value is between 50 and 250, the feature value is normalized to 2. Setting the peak value criterion normalizes the discretized vibration signal feature value to a limited number of values, thereby reducing the subsequent calculation workload.

[0108] A-4-5. Use maximum likelihood estimation to calculate the conditional probability of each feature in each category

[0109] P(X i |C n ),Right now:

[0110] P(X1|C1),P(X2|C1),…P(X 3m |C1);P(X1|C2),P(X2|C2),…P(X 3m |C n );

[0111] ...; P(X1|C2), P(X2|C n ),…P(X 3m |C H );

[0112] A-4-6. Initial setting of each category C i The probability of each category appearing is equal, that is, the probability of each category appearing is equal.

[0113] P(C1)=P(C2)=…=P(C H )=1 / H.

[0114] The monitoring phase includes:

[0115] B. Real-time data collection steps.

[0116] B-1. Install the triaxial accelerometer in the same manner and position as in step A-1. Collect data at the same sampling frequency as in step A-2, and then upload the acceleration data collected by the triaxial accelerometer in real time.

[0117] B-2. The instrument management terminal receives and stores the received data.

[0118] C. Collect data and analyze it, the steps involved are as follows Figure 2 As shown, the following is a detailed explanation of the process:

[0119] C-1. The collected data is pre-processed in the same manner as in step A-3.

[0120] C-2. Calculate the displacement MX of the processed signal. If the displacement MX is greater than a set motion threshold, such as 1 cm, the instrument or sensor is being moved as a whole. Generally, the instrument is not functioning properly at this point, and subsequent testing is not performed. If no displacement is detected, proceed to the next step.

[0121] In this step, the displacement data can be directly calculated using the collected acceleration information. This is common physical knowledge in this field, and the specific calculation method will not be repeated here.

[0122] C-3. Analyze the number of times each usage scenario occurs in this set of data.

[0123] C-3-1. Slice the data in the same way as in step A-4-2, calculate all eigenvalues ​​for each slice, and output the normalized eigenvectors based on the feature criteria established in step A-4-4. For example, if the X-axis peak value in this slice is 423, the Y-axis peak value is 38, and the Z-axis peak value is 25, then according to the previously set interval division criteria: 1: 0-50, 2: 50-250, 3: 250-450, 4: 450-650, the feature data can be output as follows:

[0124] y=(4,1,1...y i ),i=1,2,3...3m.

[0125] C-3-2. Based on the characteristic data y calculated in step C-3-1, the conditional probability established in step A-4-5 and the prior probability in step A-4-6, calculate the probability of this group of data in different classification categories C n The first posterior probability in :

[0126]

[0127] C-3-3. Find the category with the highest probability among all categories, that is

[0128] P(C k |y)=max{P(C1|y), P(C2|y),…,P(C n |y)}, then this set of data is considered to belong to category C k ; That is, the output instrument is currently in state C k middle.

[0129] C-4. The backend counts the number of times each state occurs, recording the frequency of instrument usage in each applicable state and thus deriving a second-order posterior probability of each state occurring. For example, each time an instrument button is pressed, the cumulative number of times the button is pressed (MN) increases by 1. Simultaneously, the backend also counts the total number of times the instrument is used (MA).

[0130] Thus, the second posterior probability PN=MN / MA of the button pressed state can be obtained;

[0131] C-5. The posterior probability is continuously updated during the detection process, specifically

[0132] The first posterior probability is used to replace the conditional probability obtained in step A-4-5, and the second posterior probability is used to replace the probability of occurrence of each category in step A-4-6. In the subsequent detection process, the updated data is used for detection, so as to achieve the purpose of continuous updating and repeated adjustment, making the results more and more accurate.

[0133] Example 1

[0134] This embodiment provides a method for classifying instrument usage based on instrument vibration signals, including the following steps:

[0135] A. Preparation steps:

[0136] A-1 The layout diagram of each module of the entire system is as follows Figure 1 As shown, it includes a three-axis acceleration sensor, a device management terminal, and an instrument for testing. The white cube in the figure is the instrument to be tested. In this embodiment, the instrument is an oscilloscope. It should be noted that the oscilloscope includes knobs and buttons, such as Figure 1The instrument surface is shown as a rectangle and a cylinder. The dark gray rectangle above the instrument in the figure is a three-axis acceleration sensor module, which is used to collect acceleration information of the three axes of X, Y and Z, and send it to the instrument management terminal. In this embodiment, the X-axis of the three-axis acceleration sensor is parallel to the instrument operation panel. On the right is the device management terminal, which is specifically a fixed computer device, including components such as a host, a display, a keyboard and a mouse. The desktop computer is connected to the instrument management terminal via a wired or wireless network, and is used to receive, store and process data collected by the sensor, and perform functions such as data analysis, status classification and resource management.

[0137] A-2 collects acceleration data samples from the X, Y, and Z axes of the triaxial accelerometer during operation and uploads them to the instrument management terminal. These acceleration data samples refer to the X, Y, and Z axis acceleration values ​​of the instrument under different operating conditions. The instrument's vibration state when in use and when not in use can be further categorized as keystroke and knob operation based on the mode of operation. The vibration state when not in use can be categorized as static, tapped on the desktop, and shaken.

[0138] When collecting acceleration sample data in step A-2, each signal type must be sampled at the same sampling frequency, w. The selection of the sampling frequency must comprehensively consider the timing constraints of button presses and rebounds, as well as the vibration characteristics of rigid materials such as the paddle spring, to ensure accurate capture of the dynamic changes in the signal. Therefore, the sampling frequency should meet the minimum requirement of w > 33 Hz to ensure data accuracy and reliability. In this embodiment, the sampling rate is 35 Hz.

[0139] When collecting data in step A-2, the single variable principle must be adhered to, meaning that a set of data represents only one operating mode. To more closely resemble actual conditions, all identical operable components on the device panel can be used. This embodiment includes static, key, knob, desktop tap, and desktop shake operations. After determining the operating mode, the sampled data is uploaded to the desktop computer in real time.

[0140] The vibration signal obtained by the A-3 sensor contains external vibration data and the sensor's own noise, so it is necessary to select an appropriate algorithm to perform the same preprocessing on each set of data, such as decompensation value processing and filtering.

[0141] A-3-1 The decompensation process employed in this embodiment is as follows: First, the average acceleration values ​​of the X, Y, and Z axes, sampled m = 2100 times (i.e., the number of samples in a minute) when the instrument is stationary, are calculated. The average values ​​for the X axis are -592.96, the Y axis is -308, and the Z axis is 14268. The average values ​​for the corresponding axes are then subtracted from all measured acceleration data, calibrating the acceleration values ​​to fluctuate around zero, thereby eliminating the influence of the inherent bias of the triaxial acceleration sensor.

[0142] A-3-2 The specific filtering process is as follows: perform EMD decomposition on the acceleration data, calculate the mutual information entropy of adjacent IMF components, find the boundary between the high-frequency and low-frequency IMF components, and then perform threshold noise reduction on the low-frequency portion. Finally, the noise-reduced signal is added to the low-frequency portion to obtain the final noise-reduced signal. This achieves filtering of the sampled data.

[0143] x(t)=imf1(t)+imf2(t)+...+imf 11 (t)+r 12 (t)

[0144] A-3-3. Calculate adjacent IMF i (t) mutual information entropy I(imf i , imf i+1 ), first calculate each imf i Energy entropy of (t)

[0145] The total energy E i is the energy entropy of the i-th component, P i is the energy of the i-th component. Then the inherent energy entropy H(IMF) corresponding to each connotation modal component is i )for:

[0146]

[0147] At this time, the mutual information entropy of adjacent components can be calculated as:

[0148] I(imf i , imf i+1 )=H(IMF i )+H(IMF i+1 )-H(IMF i , IMF i+1 )

[0149] in

[0150]

[0151] P(x,y) represents the joint probability distribution of IMFi and IMFi+1

[0152] A-3-4. The mutual information entropy of adjacent components will first decrease and then increase. Considering that I does not necessarily have only one minimum value, the information entropy obtained by A-3-4 is:

[0153] I(imf1, imf2) = 8.35, I(imf2, imf3) = 8.151, I(imf3, imf4) = 8.73, so

[0154] It can be determined that the value of k is 3. At this time, x(t) is rewritten with k=3 as the dividing point as follows:

[0155] x(t)=h i (t)+I j (t)+r 12 (t), i=1,2,3; j=4,5,...11.

[0156] A-3-5 For high frequency signal h i (t) is subjected to adaptive threshold noise reduction to obtain the high-frequency signal hp(t) after noise reduction. The noise threshold is calculated as follows: Where N is the signal length, σ is the standard deviation of the noise component, and x i is the high frequency IMF component, that is, x i =h i (t), median represents the median function. The calculation result is

[0157] T1=147.29, T2=75.08, T3=52.03. i The values ​​below the noise threshold in (t) are set to 0, and the high-frequency signal hp(t) after noise reduction is obtained.

[0158] A-3-6 adds the high-frequency denoised signal to the low-frequency signal and discards the residual component to obtain the denoised signal.

[0159] A-4 extracts the time-frequency domain features of the X, Y, and Z axes based on the sample data, and establishes a discretization standard and a naive Bayes classifier.

[0160] A-4-1 Determine the classification category. The classification categories adopted in this embodiment are specifically: C1: Instrument is stationary, C2: Instrument is pressed, C3: Instrument knob is turned, C4: Desktop is knocked, and C5: Desktop is shaken.

[0161] A-4-2 slices each sample and sets the slicing time t1 to 1 second.

[0162] A-4-3. Determine the eigenvector. In the present invention, the types of eigenvectors include but are not limited to peak value, variance, kurtosis and skewness in the time domain, as well as the main peak frequency, main peak intensity and number of spectral peaks in the frequency domain after calculating its power spectrum. The above eigenvectors need to extract the features of the X, Y and Z axes respectively, that is: X i =(X1: X-axis peak value, X2: Y-axis peak value, X3: Z-axis peak value, X4: X-axis variance, X5: Y-axis variance, X6: Z-axis variance, ..., X 3m-2 : Number of X-axis peaks, X 3m-1 : Number of peaks on the Y axis, X 3m: number of Z-axis peaks), i = 1, 2…3m. Calculate the values ​​of the above features in the slice to obtain the feature vector.

[0163] A-4-4 establishes discretization criteria. Because the characteristic values ​​of the collected vibration signal fluctuate within a range of threshold values, the threshold intervals need to be discretized. Specifically, appropriate discretization algorithms are selected for different features, discretization criteria are established, and segmented thresholds are established for each feature. The established segmented threshold criteria are shown in Table 1. The threshold criteria here are valid only for this embodiment.

[0164] Table 1

[0165]

[0166] A-4-5. Use maximum likelihood estimation to calculate the conditional probability of each feature in each category. This yields the conditional probability for each feature, as shown in Tables 2 through 4. In Tables 2, 3, and 4, P represents the conditional probability, X represents the feature value, and C represents the different categories. The conditional probabilities here are valid only for this embodiment; the specific conditional probability values ​​depend on the preprocessing results of the collected data.

[0167] Table 2

[0168]

[0169] Table 3

[0170]

[0171] Table 4

[0172]

[0173] A-4-6. Assume that the probability of each category appearing in each category is equal, that is, the prior probability

[0174] P(C1)=P(C2)=…=P(C5)=1 / 5.

[0175] B. Real-time data collection:

[0176] The installation method and position of the B-1 triaxial acceleration sensor are the same as those in preparation step A-1, and then the data collected by the instrument are uploaded in real time.

[0177] The B-2 instrument management terminal receives and stores the received data.

[0178] C. Data collection and analysis:

[0179] C-1 pre-processes the collected data in the same manner as step A-3.

[0180] C-2 calculates the displacement of the processed signal and obtains a displacement of 0.003 cm. This indicates that the instrument has not moved, so the subsequent steps are continued.

[0181] C-3 analyzes the number of times each usage scenario occurs in the data set.

[0182] C-3-1 performs the same slicing process on the data as in A-4-2, calculating all eigenvalues ​​for each slice and outputting a feature vector X based on the feature criteria established in A-4-4, creating the data to be classified y. X = (453, 38, 18, 24588, 317, 275, 10, 5, 6, -2, -0.2, 0.1). The output partial feature data y = (4, 2, 1, 3, 1, 1, 3, 2, 2, 1, 3, 2).

[0183] C-3-2 calculates the posterior probability of the data set in different classification categories based on the calculated feature data y, the conditional probability established in step A-4-5, and the prior probability in step A-4-6.

[0184] The calculation results are:

[0185] P(C1|y)=1.3300*10 -9 , P(C2|y)=0.47235, P(C3|y)=0.27387,

[0186] P(C4|y)=0.12777, P(C5|y)=0.12601.

[0187] C-3-3 finds the category with the highest probability among all categories. Based on C-3-2, we can conclude that the maximum value is P(C2|y), so the output shows that the instrument is currently in state C2. This means that the instrument is in the key-pressing state, which is the same as actual operation.

[0188] The foregoing are the preferred embodiments of the present invention. Unless the preferred implementation modes in each preferred embodiment are obviously self-contradictory or based on a certain preferred implementation mode, each preferred implementation mode can be arbitrarily superimposed and used in combination. The embodiments and the specific parameters in the embodiments are only for the purpose of clearly describing the inventor's invention verification process, and are not intended to limit the patent protection scope of the present invention. The patent protection scope of the present invention shall still be based on its claims. Any equivalent structural changes made using the contents of the description and drawings of the present invention should also be included in the protection scope of the present invention.

Claims

1. A device status monitoring method based on vibration signals, characterized in that , including an initial sample collection phase and a monitoring phase, wherein the initial sample collection phase includes the following steps: A-1. Place the triaxial accelerometer, device management terminal, and at least one device under test in the same positions as they will be during normal use, allowing the triaxial accelerometer to detect the motion of the device under test. Set the triaxial directions of the triaxial accelerometer. A-2. Collect X, Y and Z axis acceleration data samples of the three-axis acceleration sensor at the same sampling frequency. The acceleration data samples refer to the acceleration data samples of the tested instrument in different working states C i The acceleration values ​​of the X, Y and Z axes under each working state are collected multiple times and the average is taken; A-3. The device management terminal performs preprocessing and denoising on the acceleration data collected in step A-2 to obtain a denoised signal sample; A-4. Build a naive Bayes classifier to calculate the initial probability distribution, specifically: A-4-1. Determine the classification category C = (C1, C2…C H ), where C i Indicates the categories corresponding to different working states, and H is the total number of states; A-4-2. Slice each denoised signal sample obtained in step A-3, where the slicing time is greater than the time required for the basic operation of the tested instrument; A-4-3. Calculate the feature vectors in each slice time period, including: X i =(X1: X-axis peak value, X2: Y-axis peak value, X3: Z-axis peak value, X4: X-axis variance, X5: Y-axis variance, X6: Z-axis variance, ..., X 3m-2 : Number of X-axis peaks, X 3m-1 : Number of peaks on the Y axis, X 3m : number of Z-axis spectrum peaks), i = 1, 2…3m, m is the total number of eigenvector types; A-4-4. Establish a discretization standard, divide different vibration signal features into different intervals, normalize the vibration signal features in the same interval to the same eigenvalue, and obtain the normalized eigenvector X i A-4-5. Use the maximum likelihood estimation method to estimate the normalized eigenvector X i Calculate the conditional probability P(X i |C n ),Right now: P(X1|C1),P(X2|C1),…P(X 3m |C1);P(X1|C2),P(X2|C2),…P(X 3m |C n ); …;P(X1|C2),P(X2|C n ),…P(X 3m |C H ); A-4-6. Initially, the probability of each category appearing in each category is equal, that is, the probability of each category appearing P(C1)=P(C2)=…=P(C H )=1 / H; The monitoring phase includes: B. Real-time data collection: Use the same placement method and sampling frequency as step A-1 and step A-2 to collect data and upload it in real time; C. Collect data and analyze it, including: C-1. The collected data is preprocessed and denoised in the same manner as in step A-3; C-3. Analyze the number of times each usage scenario occurs in the data, specifically: C-3-1. Perform the same slicing process on the data set as in step A-4-2, calculate all eigenvalues ​​for each slice, and output the normalized eigenvector y according to the discretization standard established in step A-4-4; C-3-2. Based on the normalized feature vector y calculated in step C-3-1, and the conditional probability established in step A-4-5 and the prior probability in step A-4-6, calculate the probability of the data set in different classification categories C n The first posterior probability in : C-3-3. Find the category with the highest probability among all categories, that is P(C k |y)=max{P(C1|y), P(C2|y),…,P(C n |y)}, then this set of data is considered to belong to category C k ; That is, the output instrument is currently in state C k middle.

2. The device status monitoring method according to claim 1, wherein: After step C-1 and before step C-3, step C-2 is also included. The displacement MX of the processed signal is calculated. If it is found that the displacement MA is greater than the set movement threshold, the subsequent steps are terminated, otherwise the subsequent steps are performed.

3. The device status monitoring method according to claim 1, wherein: The A-3 step specifically includes: A-3-1. De-compensate the acceleration data by first calculating the average acceleration of multiple samples taken along the X, Y, and Z axes when the instrument is stationary. Then, subtract the average of the corresponding axis from all measured acceleration data while stationary. This allows the corrected acceleration value to fluctuate around zero, eliminating the effects of inherent sensor bias. A-3-2. Perform EMD decomposition on the acceleration data samples to obtain several intrinsic modal components imf i , i=1,2…n-1. i represents the serial number of different frequency components, n is the total number of frequency components, and the acceleration signal x(t) is expressed as: x(t)=imf1(t)+imf2(t)+...+imf n-1 (t)+r n (t); Among them, imf i (t) represents the intrinsic modal components with different frequency components, arranged from high frequency to low frequency, r n (t) is the residual component, t represents time; A-3-3. Calculate adjacent IMF i (t) mutual information entropy I(imf i , imf i+1 ), first calculate each imf i Energy entropy of (t) The total energy E i is the energy entropy of the i-th component, P i is the energy of the i-th component. Then the inherent energy entropy H(IMF) corresponding to each connotation modal component is i )for: The mutual information entropy of adjacent components is: I(imf i ,imf i+1 )=H(IMF i )+H(IMF i+1 )-H(IMF i ,IMF i+1 ) in P(x,y) represents the joint probability distribution of IMFi and IMFi+1; A-3-4. Find I(imf i , imf i+1 ) is the first minimum turning point, the minimum turning point is defined as: if I(imf i , imf i+1 ) decreases and I(imf i , imf i+2 ) increases, then the i value at this time is the minimum turning point, and the first minimum value that appears is the high-low frequency dividing point. Let the minimum turning point be k, and rewrite x(t) with k as the dividing point as follows: x(t)=h i (t)+I j (t)+r n (t) where h i (t) is the high frequency signal, i=1,2…k, I j (t) is the low-frequency signal, j=k+1, k+2…n-1. A-3-5. For high frequency signal h i (t) is subjected to adaptive threshold denoising to obtain the denoised high-frequency signal hp(t); the noise threshold is calculated as follows: Where N is the signal length, σ is the standard deviation of the noise component, and x i is the high frequency IMF component, that is, x i =h i (t), median represents the median function; Calculate h separately i (t) The corresponding noise threshold T i , h i (t) is below the noise threshold T i The value of is set to 0 to obtain the high-frequency signal hp(t) after noise reduction; A-3-6. Add the denoised high-frequency signal to the low-frequency signal and discard the residual component r n Get the denoised signal sample: xp(t)=hp(t)+I j (t)。 4. The device status monitoring method according to claim 3, wherein: In step A-3, before step A-3-2 begins, step A-3-1 is also included. The acceleration data samples are processed to remove the compensation value. The specific method is: first calculate the average acceleration of m samples of the X, Y, and Z axes when the instrument is stationary, then subtract the average of the corresponding axis from all the measured acceleration data when it is stationary, and then proceed to step A-3-2.

5. The device status monitoring method according to claim 3, wherein: The specific steps of A-3-5 are: Calculating the noise threshold Where N is the signal length, σ is the standard deviation of the noise component, and x i is the high frequency IMF component, that is, x i =h i (t), median represents the median function; Calculate h separately i (t) The corresponding noise threshold T i , h i (t) is below the noise threshold T i The value of is set to 0 to obtain the high-frequency signal hp(t) after noise reduction.

6. The device status monitoring method according to claim 1, wherein: Step C also includes: C-4. Count the number of occurrences of each state and the total number of times the instrument is used (MA) to obtain the second posterior probability of different states; C-5. Continuously update the posterior probability during the detection process, specifically: The first posterior probability is used to replace the conditional probability obtained in step A-4-5, and the second posterior probability is used to replace the occurrence probability of each category in step A-4-6. The updated data is used for detection in the subsequent detection process.

Images