Magnetic levitation bearing gap millimeter wave measuring device and compensation method

Abstract

The present application relates to the technical field of magnetic suspension bearing, and discloses a magnetic suspension bearing gap millimeter wave measuring device and a compensation method, which comprises the following steps: deploying a millimeter wave radar array on the inner wall of a magnetic suspension bearing stator to perform multi-point gap measurement and speed extraction, and obtaining millimeter wave data; collecting multi-source data, calculating confidence weight, and obtaining fusion data; performing space-time consistency verification and anomaly detection, and adjusting the confidence weight; performing adaptive harmonic extraction and rotation frequency synchronous wave trapping, and extracting a synchronous disturbance component; generating a compensation force and generating a total control instruction; the present application improves the measurement reliability and environmental adaptability of the magnetic suspension bearing in a harsh industrial environment by constructing a multiple-input multiple-output measuring array based on a millimeter wave radar sensor, and can fully capture the circumferential gap distribution characteristics of the rotor by combining a multiple signal classification spatial spectrum estimation algorithm to reconstruct the gap spatial distribution field, thereby improving the dynamic stiffness and suspension stability of the magnetic suspension bearing in the full rotation speed range.

Description

Technical Field

[0001] This invention relates to the field of magnetic levitation bearing technology, and more specifically, to a millimeter-wave measurement device and compensation method for the gap of magnetic levitation bearings. Background Technology

[0002] Magnetic levitation bearings achieve contactless support through electromagnetic force and are widely used in high-speed, high-precision applications such as high-speed electric spindles, vacuum pumps, and energy storage flywheels.

[0003] Chinese Patent CN113659911B discloses a rotor vibration suppression system and method in a magnetic levitation bearing system. When the rotor is subjected to a disturbance force, the rotor disturbance compensation observer receives the rotor vibration displacement signal and the rotor vibration acceleration signal, and then sends a rotor vibration suppression compensation control signal. The control current is then adjusted by a power amplifier to provide electromagnetic force. When the base is subjected to a disturbance force, the base acceleration feedforward compensator receives the base vibration acceleration signal and sends a base acceleration feedforward compensation signal. The position controller then provides a reference displacement. When the magnetic levitation bearing system becomes uncontrollable, the position error control channel is cut off, and a rotor vibration suppression control signal after instability is added to suppress the full clearance eddy current generated by the rotor in the protective bearing. At the same time, the rotor disturbance compensation observer sends a rotor vibration suppression compensation control signal, and the control current of the magnetic levitation bearing is adjusted by a power amplifier to provide electromagnetic force.

[0004] However, in real production environments, especially in high-speed machining centers, oil mist and cutting fluid vapor are commonly present inside the machine tool housing, along with flying metal chips and particles. Temperature and humidity are high, and condensation may occur. Commonly used optical or eddy current displacement sensors are prone to bias and drift under the cover of oil mist, water film, and contaminant layers. These factors significantly reduce the measurement reliability of existing technologies. The double differentiation will highly amplify noise and quantization errors, causing acceleration estimation to tend to be distorted in noisy environments. Notch filters are mostly designed for steady-state single-frequency operation and are difficult to maintain effective extraction under rapid speed changes and multi-harmonic interference. This leads to instability and decreased reliability of gap measurements under critical operating conditions, making it difficult to support high-precision, long-term, and continuous operation without speed reduction. Summary of the Invention

[0005] This invention provides a millimeter-wave measurement device and compensation method for the clearance of magnetic levitation bearings to solve the above problems.

[0006] The present invention provides a millimeter-wave measurement device and compensation method for the clearance of magnetic levitation bearings, comprising:

[0007] Millimeter-wave data is obtained by multi-point gap measurement and velocity extraction using a millimeter-wave radar array deployed on the inner wall of the stator of the magnetic levitation bearing.

[0008] Collect multi-source data, calculate confidence weights based on the multi-source data and millimeter-wave data, and perform confidence assessment and fusion based on the confidence weights to obtain fused data;

[0009] Based on the millimeter-wave data and fused data, spatiotemporal consistency verification and anomaly detection are performed, and the confidence weights are adjusted according to the detection results.

[0010] Adaptive harmonic extraction and frequency conversion synchronization notch filtering are performed on the fused data to extract synchronization disturbance components;

[0011] Based on the fused data and synchronous disturbance components, predictive state observation and look-ahead compensation are performed to generate compensation force and generate overall control command.

[0012] Furthermore, the multi-point gap measurement and velocity extraction include:

[0013] Millimeter-wave radar sensors are deployed at a predetermined number of equally spaced circumferential positions on the inner wall of the magnetic levitation bearing stator to form a multi-input multi-output measurement array;

[0014] Each radar sensor transmits a linear frequency modulated continuous wave signal, and the received echo signal is mixed with the local oscillator signal to obtain the intermediate frequency signal.

[0015] Millimeter-wave data includes millimeter-wave measurement gap and micro-Doppler velocity. Based on the intermediate frequency signal, the beat frequency is extracted through Fourier transform, the millimeter-wave measurement gap is calculated based on the beat frequency, and the micro-Doppler velocity is extracted through phase change. The millimeter-wave measurement gap is calculated by multiplying the preset speed of light by the preset sweep period by the beat frequency and then dividing by twice the preset sweep bandwidth. The micro-Doppler velocity is calculated by multiplying the preset millimeter-wave wavelength by the phase difference between adjacent sweep periods, dividing by four times pi and then dividing by the preset sweep period.

[0016] Furthermore, the confidence assessment and fusion include:

[0017] The multi-source data includes signal power, noise power, and eddy current displacement sensor output; the confidence weights include millimeter-wave measurement confidence weights and sensor confidence weights; and the fused data includes gap estimates and fusion velocity estimates.

[0018] The signal-to-noise ratio and multipath interference index are calculated based on the signal power and noise power, and the confidence weight of millimeter-wave measurement is calculated through the signal-to-noise ratio and multipath interference index.

[0019] Based on the output and saturation range of the eddy current displacement sensor, the confidence weight of the sensor is calculated.

[0020] Based on the millimeter-wave measurement confidence weights and sensor confidence weights, an extended Kalman filter state vector and measurement vector are constructed, and a weighted fusion update is performed to obtain the fusion gap estimate and fusion speed estimate.

[0021] Furthermore, the weighted fusion update includes:

[0022] The measurement vector includes a first element, a second element, and a third element. The first element is the millimeter-wave measurement confidence weight multiplied by the millimeter-wave average gap. The second element is the sensor confidence weight multiplied by the eddy current displacement sensor output. The third element is the average micro-Doppler velocity.

[0023] The state vector includes a first state element, a second state element, a third state element, and a fourth state element. The first state element is the fusion gap, the second state element is the fusion velocity, the third state element is the fusion acceleration, and the fourth state element is the eccentricity angle. The fusion velocity represents the first derivative of the fusion gap with respect to time, the fusion acceleration represents the second derivative of the fusion gap with respect to time, and the eccentricity angle represents the eccentricity angle of the rotor attitude.

[0024] The state vector is updated by the state transition equation. The state vector at the next time step is equal to the state transition matrix multiplied by the state vector at the current time step, plus the process noise vector.

[0025] The measurement vector is updated by the measurement equation. The measurement vector at time t is equal to the measurement matrix multiplied by the state vector at time t, plus the measurement noise vector.

[0026] The state estimate is updated by the Kalman gain matrix to obtain the fusion gap estimate and the fusion velocity estimate, where the fusion gap estimate is equal to the fusion gap and the fusion velocity estimate is equal to the fusion velocity.

[0027] Furthermore, the signal-to-noise ratio (SNR) is 10 multiplied by a logarithmic function to the base 10, where the independent variable of the logarithmic function is the ratio of signal power to noise power; the multipath interference index is the sum of the absolute amplitude values ​​of the multipath components from the second to the Mth multipath components, divided by the absolute amplitude value of the first multipath component, where M is the total number of multipath components; the millimeter-wave measurement confidence weight is the sigmoid function applied to the comprehensive index, where the comprehensive index is the SNR term minus the multipath interference term, the SNR term is a preset SNR adjustment coefficient multiplied by the millimeter-wave SNR, and the multipath interference term is a preset multipath interference adjustment coefficient multiplied by the multipath interference index;

[0028] The saturation range is a closed interval from the preset lower saturation limit to the preset upper saturation limit. When the absolute value of the difference between the output of the eddy current displacement sensor at the current moment and the output at the previous moment is less than the preset normal change threshold, and the output of the eddy current displacement sensor at the current moment is not within the saturation range, the sensor confidence weight is set to 1, where the previous moment is the current moment minus the preset time interval; otherwise, the sensor confidence weight is a natural exponential function with the product of the negative value of the preset attenuation coefficient and the abnormal change amount as the exponent, where the abnormal change amount is the absolute value of the difference between the output of the eddy current displacement sensor at the current moment and the output at the previous moment minus the preset normal change threshold.

[0029] Furthermore, the spatiotemporal consistency verification and anomaly detection include:

[0030] The standard deviation of the circumferential clearance is calculated. When the standard deviation of the circumferential clearance is greater than a preset spatial consistency threshold, or when the smallest millimeter-wave measurement clearance among a preset number of circumferential positions is less than a preset clearance warning threshold, a spatial anomaly flag is triggered. The standard deviation of the circumferential clearance is the square root of the clearance square term divided by a preset number. The clearance square term is the sum of the clearance difference terms at the preset number of circumferential positions. The clearance difference term is the difference between the millimeter-wave measurement clearance and the average millimeter-wave clearance.

[0031] Calculate the velocity residual. When the velocity residual is greater than the preset velocity consistency threshold, trigger the velocity inconsistency flag. The velocity residual is the absolute value of the fused velocity estimate minus the average microDoppler velocity.

[0032] When a spatial anomaly flag or a velocity inconsistency flag is triggered, the confidence weights are adjusted and normalization is performed.

[0033] Furthermore, adjusting the confidence weights includes:

[0034] When only the spatial anomaly flag is triggered and the velocity inconsistency flag is not triggered, the sensor confidence weight is updated to the original sensor confidence weight multiplied by the preset sensor weight attenuation factor, and the millimeter wave measurement confidence weight is updated to the original millimeter wave measurement confidence weight multiplied by the preset millimeter wave weight enhancement factor.

[0035] When the speed inconsistency flag is triggered, regardless of whether the spatial anomaly flag is triggered at the same time, the millimeter-wave measurement confidence weight is updated to the original millimeter-wave measurement confidence weight multiplied by the preset millimeter-wave weight attenuation factor, and the sensor confidence weight is updated to the original sensor confidence weight multiplied by the preset sensor weight enhancement factor.

[0036] Furthermore, the adaptive harmonic extraction and frequency-switching synchronous notch filtering includes:

[0037] A notch filter bank is constructed based on real-time frequency conversion, and the notch filter bank includes multiple notch filters.

[0038] A notch filter bank is applied to the estimated fusion gap to extract the synchronization disturbance component. The synchronization disturbance component is composed of the superposition of harmonic components from the fundamental frequency to the preset harmonic order. Each harmonic component is the amplitude of the kth harmonic multiplied by a sine function. The independent variable of the sine function is the product of k times the real-time frequency and time t plus the phase of the kth harmonic. The amplitude and phase of the kth harmonic are estimated in real time using the least mean square adaptive algorithm.

[0039] Furthermore, the predictive state observation and look-ahead compensation include:

[0040] A rotor dynamics state-space model is established. The dynamics state-space model is as follows: the derivative of the rotor dynamics state vector with respect to time is equal to the sum of the product of the rotor dynamics system matrix and the rotor dynamics state vector, the product of the rotor dynamics input matrix and the electromagnetic force command, and the product of the rotor dynamics disturbance matrix and the external disturbance.

[0041] Using the Luneburger observer, the current rotor dynamic state estimate and external disturbance estimate are estimated based on the fusion gap estimate and fusion velocity estimate;

[0042] Perform state prediction at a preset prediction time step to obtain the predicted rotor dynamics state vector;

[0043] The compensation force is calculated based on the predicted gap and the preset safety threshold. The predicted gap is the first element of the predicted rotor dynamic state vector.

[0044] Generate a general control command, which is the sum of the output of the traditional PID controller, the compensation force, the notch compensation term, and the disturbance feedforward compensation term. The notch compensation term is the negative value of the preset notch gain multiplied by the synchronous disturbance component, and the disturbance feedforward compensation term is the negative value of the external disturbance estimate.

[0045] Furthermore, the predicted rotor dynamics state vector is the product of the first matrix exponential function and the current rotor dynamics state estimate, plus a first integral term, plus a second integral term. The first integral term is the integral from 0 to the preset prediction time step, and its integrand is the second matrix exponential function multiplied by the rotor dynamics input matrix and then by the electromagnetic force command. The second integral term is the integral from 0 to the preset prediction time step, and its integrand is the third matrix exponential function multiplied by the rotor dynamics disturbance matrix and then by the external disturbance estimate. The base of the first matrix exponential function is a natural constant, and its exponent is the rotor dynamics system matrix multiplied by the preset prediction time step. The base of the second matrix exponential function is a natural constant, and its exponent is the integral variable of the rotor dynamics system matrix multiplied by the first integral term. The base of the third matrix exponential function is a natural constant, and its exponent is the integral variable of the rotor dynamics system matrix multiplied by the second integral term.

[0046] The compensation force is the difference between the preset safety threshold and the predicted gap multiplied by the preset position gain, plus the velocity term, which is the preset velocity gain multiplied by the predicted velocity, where the predicted velocity is the second element of the predicted rotor dynamic state vector.

[0047] The beneficial effects of this invention are as follows: By constructing a multi-input multi-output measurement array based on an industrial-grade frequency-modulated continuous wave millimeter-wave radar sensor, this invention significantly improves the measurement reliability and environmental adaptability of magnetic levitation bearings in harsh industrial environments. Utilizing the high spatial resolution and excellent penetration characteristics of the millimeter-wave band, this device can effectively penetrate oil mist and condensation generated during high-speed operation and resist the influence of strong electromagnetic interference, solving the problem that traditional optical or capacitive displacement sensors are prone to contamination failure or signal distortion under complex working conditions. Simultaneously, by combining a multi-signal classification spatial spectrum estimation algorithm to reconstruct the gap spatial distribution field, it can capture the circumferential gap distribution characteristics of the rotor from all directions, achieving high-precision, full-circumferential monitoring of the rotor's attitude eccentricity angle and eccentricity, overcoming the limitations of single-point measurement methods in spatial coverage.

[0048] This invention employs multi-source data confidence assessment and fusion technology, significantly enhancing the robustness and data accuracy of the measurement system. By real-time acquisition and analysis of millimeter-wave measurement data and eddy current displacement sensor data, and leveraging the complementary advantages of their physical characteristics and error distributions, a comprehensive confidence assessment model is constructed, incorporating signal-to-noise ratio, multipath interference indicators, and sensor saturation state. Based on this model, an adaptive weighted fusion mechanism automatically reduces the weight of a particular type of sensor and increases the weight of high-confidence sensors when a sensor experiences specific interference, signal attenuation, or enters the nonlinear saturation region. This dynamic adjustment strategy effectively avoids measurement deviations caused by single sensor failures or sudden environmental changes, ensuring that the fusion gap and fusion speed estimates remain at the optimal confidence level.

[0049] This invention significantly optimizes the high-speed dynamic performance and control accuracy of a magnetic levitation bearing system through an integrated adaptive harmonic extraction and predictive state observation compensation method. Utilizing an adaptive notch filter bank based on real-time rotational frequency, it accurately extracts and filters out periodic disturbance components synchronized with the rotational speed, eliminating harmonic noise introduced by rotor mass imbalance and sensor self-jumping. More importantly, by employing a rotor dynamic state observer based on a Luneburger observer and a look-ahead prediction algorithm, this method can estimate the rotor's dynamic state in real time and effectively compensate for computational delays, sampling delays, and actuator response delays in the control loop. This look-ahead compensation mechanism solves the problem of decreased system stability margin caused by phase lag during high-speed operation, greatly improving the dynamic stiffness and levitation stability of the magnetic levitation bearing across the entire speed range. Attached Figure Description

[0050] Figure 1 This is a flowchart illustrating the magnetic levitation bearing clearance millimeter-wave measurement device and compensation method of the present invention.

[0051] Figure 2 This is an example diagram illustrating the fusion data obtained by the magnetic levitation bearing gap millimeter-wave measurement device and compensation method of the present invention;

[0052] Figure 3 This is an example diagram illustrating the extraction of synchronous disturbance components in the magnetic levitation bearing clearance millimeter-wave measurement device and compensation method of the present invention;

[0053] Figure 4 This is an example diagram of the overall control command for generating the magnetic levitation bearing gap millimeter-wave measurement device and compensation method of the present invention. Detailed Implementation

[0054] The subject matter described herein will now be discussed with reference to exemplary embodiments. It should be understood that these embodiments are discussed only to enable those skilled in the art to better understand and implement the subject matter described herein, and changes may be made to the function and arrangement of the elements discussed without departing from the scope of this specification. Various processes or components may be omitted, substituted, or added as needed in the examples. Furthermore, features described in some examples may be combined in other examples.

[0055] Millimeter-wave measurement equipment and compensation methods for the clearance of magnetic levitation bearings, such as Figure 1 As shown, it includes:

[0056] Step 100: Multi-point gap measurement and velocity extraction are performed by a millimeter-wave radar array deployed on the inner wall of the stator of the magnetic levitation bearing to obtain millimeter-wave data.

[0057] Industrial-grade frequency-modulated continuous wave millimeter-wave radar sensors are deployed at a predetermined number of equally spaced circumferential positions on the inner wall of the magnetic levitation bearing stator, forming a multi-input multi-output measurement array. The default number of sensors is 8, determined based on the circumferential dimensions of the magnetic levitation bearing stator and the required spatial resolution, ensuring that the circumferential clearance distribution characteristics of the rotor can be fully captured. Because the magnetic levitation bearing rotor generates complex dynamic eccentricity and vibration modes during high-speed rotation, single-point measurement cannot fully reflect the true clearance distribution between the rotor and stator. Multi-point measurement can reconstruct the entire circumferential clearance spatial distribution field in real time, thereby accurately identifying the rotor's eccentricity direction, eccentricity amount, and possible local anomaly proximity areas. The angular positions of these predetermined number of circumferential positions are: the angular position of the i-th position is equal to i multiplied by the predetermined angular interval, where i ranges from 1 to the predetermined number, and the predetermined angular interval is equal to 360 degrees divided by the predetermined number. When the predetermined number is 8, the predetermined angular interval is 45 degrees, meaning that a radar sensor is deployed every 45 degrees circumferentially. Equally spaced deployment ensures uniformity in spatial sampling, avoiding spatial aliasing caused by uneven distribution of measurement points. It also facilitates subsequent use of spatial spectrum estimation algorithms such as multiple signal classification for high-precision angle resolution and gap field reconstruction. Each radar sensor employs a preset carrier frequency, with a default value of 77 GHz. This default value selects the millimeter-wave band to achieve high spatial resolution and good penetration. The default preset sweep bandwidth is 4 GHz, determined based on the required range resolution. The default preset sweep period is 1 millisecond, determined based on the rotor's maximum speed and the required measurement update rate. The 77 GHz millimeter-wave band has a wavelength of only 3.9 mm, enabling micrometer-level range resolution, meeting the high-precision requirements for magnetic levitation bearing gap measurement. Furthermore, millimeter waves have excellent penetration capabilities through aerosol media such as oil mist and condensation, making them less susceptible to the effects of oil mist obstruction in the working environment compared to optical measurement methods such as lasers and infrared. Additionally, millimeter-wave measurement is non-contact, applying no additional mechanical or thermal load to the rotor surface, making it suitable for continuous online monitoring in high-speed rotation scenarios.

[0058] During millimeter-wave signal transmission and reception, each radar sensor synchronously transmits a linear frequency modulated continuous wave signal. The time-domain expression of the transmitted signal is a preset transmission amplitude multiplied by a cosine function. The phase angle of the cosine function consists of two parts: the first part is twice pi multiplied by the preset carrier frequency and then multiplied by time t; the second part is pi multiplied by the preset sweep bandwidth divided by the preset sweep period and then multiplied by the square of time t. This signal is valid within the range of time t from 0 to the preset sweep period. The default value of the preset carrier frequency is 77 GHz, and the preset transmission amplitude is the amplitude value corresponding to the preset transmission power. The default value of the preset transmission power is determined according to the measurement distance range and receiving sensitivity requirements, and is usually set to 10 dB / mW to ensure sufficient echo signal strength within the gap range of the magnetic levitation bearing. Linear frequency modulated (LFM) continuous wave radar can simultaneously achieve high-precision distance and velocity measurements. By transmitting a continuous wave signal with a frequency that varies linearly with time, the beat frequency generated after the received echo is mixed with the local oscillator signal directly reflects the target distance, while the phase change between adjacent sweep cycles reflects the target's radial velocity. This dual measurement capability allows for the simultaneous acquisition of rotor position and motion state information during a single sweep. Compared to pulse radar, LFM continuous wave radar offers advantages such as lower transmission power, simpler circuit implementation, and lower cost. Furthermore, due to its continuous wave operation, it achieves higher average transmission power and a better signal-to-noise ratio, making it particularly suitable for short-range, high-precision measurement scenarios. The necessity of synchronous transmission among radar sensors ensures phase consistency between multiple measurement channels, a prerequisite for subsequent virtual array construction and spatial spectrum estimation. If the sensors are not synchronized in their transmission times, additional phase errors will be introduced, leading to decreased accuracy in spatial distribution field reconstruction or even the appearance of false peaks.

[0059] During the intermediate frequency (IF) signal acquisition process, echo signals from the rotor surface are received. The received echo signals are mixed with the local oscillator signal to obtain the IF signal. The time-domain expression of the IF signal is the IF amplitude multiplied by a cosine function. The phase angle of the cosine function consists of two parts: the first part is twice pi multiplied by twice the preset sweep bandwidth divided by the preset speed of light, then divided by the preset sweep period, multiplied by the distance from the rotor surface to the stator inner wall, and then multiplied by time t. The second part is the preset initial phase. The default value of the preset speed of light is 3 x 10⁸ meters per second, which is the speed of electromagnetic waves in a vacuum. The preset initial phase is a fixed phase offset introduced during the mixing process. The default value of the preset initial phase is determined through the system calibration process, typically obtained after system installation using a standard reflector.

[0060] The use of mixing to acquire intermediate frequency (IF) signals allows for down-conversion of high-frequency millimeter-wave echo signals to a lower IF range, significantly reducing the sampling rate requirements and computational complexity of subsequent analog-to-digital conversion and digital signal processing. Simultaneously, the mixing process preserves the distance and velocity information contained in the echo signal, enabling the system to achieve high-precision measurements at a lower hardware cost. The frequency of the IF signal is directly proportional to the distance from the rotor surface to the stator inner wall. This linear relationship simplifies distance calculation; a fast Fourier transform of the IF signal directly retrieves the distance information from the spectral peak position. Calibration of the preset initial phase is crucial for ensuring measurement accuracy. This phase shift is introduced by various factors, including local oscillator leakage from the mixer, circuit delay, and antenna phase center deviation. By using a standard reflector at a known distance for calibration after system installation, the value of the preset initial phase can be accurately determined. This fixed phase shift is then compensated for in subsequent measurements, eliminating systematic errors and achieving absolute rather than relative gap measurement.

[0061] In the distance and velocity calculation process, a fast Fourier transform (FFT) is performed on the intermediate frequency signal in the fast time dimension to extract the beat frequency corresponding to the peak value of the distance dimension spectrum. The millimeter-wave measurement gap is calculated by multiplying the preset speed of light by the preset sweep period by the beat frequency, and then dividing the result by twice the preset sweep bandwidth. Micro-Doppler velocity is extracted through phase change in the slow time dimension. The micro-Doppler velocity is calculated by multiplying the preset millimeter-wave wavelength by the phase difference between adjacent sweep periods, dividing the result by four times pi, and then dividing by the preset sweep period. The preset millimeter-wave wavelength is equal to the preset speed of light divided by the preset carrier frequency. When the preset carrier frequency is 77 GHz, the preset millimeter-wave wavelength is 3.9 mm. The FFT for distance calculation has extremely high computational efficiency, completing spectrum analysis in milliseconds, meeting the stringent real-time requirements of the magnetic levitation bearing control system. Simultaneously, the FFT can obtain the spectrum information of all distance units, facilitating the identification of multipath reflections and clutter interference. The linear relationship between beat frequency and gap distance eliminates the need for complex nonlinear corrections in distance calculation. Measurement accuracy primarily depends on the estimation accuracy of spectral peaks and the stability of the preset sweep bandwidth. Micro-Doppler velocity extraction leverages the continuity of echo signal phase between adjacent sweep cycles. By calculating the phase difference, the radial velocity of the rotor surface relative to the radar sensor can be obtained. This velocity measurement method offers advantages such as high accuracy, large dynamic range, and immunity to distance measurement errors. The necessity of micro-Doppler velocity information is twofold: firstly, it can be cross-validated with the fused velocity estimated by the Kalman filter to detect temporal consistency of measurement data and identify anomalous jumps; secondly, micro-Doppler velocity reflects the changing trend of rotor gap and can be used to predict the gap state at future moments, providing crucial input information for look-ahead compensation control.

[0062] In the spatial distribution field reconstruction process, a virtual array is constructed based on a multi-input multi-output (MIMO) measurement array formed by a preset number of radar sensors. This virtual array is equivalent to a preset number of virtual array elements, where the preset number is equal to the square of the preset number. When the preset number is 8, the preset number of virtual array elements is 64. A multi-signal classification spatial spectrum estimation algorithm is used to reconstruct the gap spatial distribution field. The gap spatial distribution field is defined as the set of millimeter-wave measurement gap values ​​at a preset number of circumferential positions. The rotor attitude eccentricity angle and eccentricity are obtained through spatial spectrum peak search. The eccentricity represents the offset distance of the rotor's geometric center relative to the stator's geometric center. By combining the transmission and reception of the MIMO array, the equivalent aperture and spatial resolution of the array can be significantly improved without increasing the number of physical array elements. When the preset number is 8, only 8 physical radar sensors are needed to obtain the spatial sampling capability equivalent to 64 virtual array elements. This virtual array technology significantly reduces hardware costs and system complexity. The advantage of employing the multi-signal classification algorithm lies in its super-resolution spatial spectrum estimation method. Its angle resolution capability is not constrained by the Rayleigh limit, enabling it to achieve angle resolution accuracy far exceeding traditional beamforming methods with a limited number of array elements. This makes it particularly suitable for applications like magnetic levitation bearings, which require precise identification of rotor eccentricity direction and magnitude. The multi-signal classification algorithm separates the signal and noise subspaces by performing eigenvalue decomposition on the array covariance matrix. It constructs a spatial spectrum function using the orthogonality between the signal steering vector and the noise subspace, and obtains high-precision angle estimation by searching for the peak position of the spatial spectrum. The reconstruction of the gap spatial distribution field allows the system to comprehensively understand the three-dimensional spatial relationship between the rotor and stator. It can not only identify the overall eccentricity of the rotor but also detect local gap anomalies, such as bulges on the rotor surface, deformation of the stator inner wall, or abnormal proximity at specific angles. This information is of significant value for fault diagnosis and predictive maintenance of magnetic levitation bearings.

[0063] Step 200: Collect multi-source data, calculate confidence weights based on the multi-source data and millimeter-wave data, perform confidence assessment and fusion based on the confidence weights, and obtain fused data, as detailed below. Figure 2 As shown.

[0064] During the multi-source data acquisition process, measurement information from different sensing principles is collected and integrated in real time. The multi-source data includes signal power, noise power, and eddy current displacement sensor output; millimeter-wave measurement data acquired by a millimeter-wave radar sensor array; and traditional displacement monitoring data acquired by an eddy current displacement sensor. The millimeter-wave measurement data is based on the industrial-grade frequency-modulated continuous wave millimeter-wave radar sensor array described in step 100. It is obtained by transmitting a linear frequency-modulated continuous wave and receiving the rotor surface echo, followed by mixing, fast Fourier transform, and multiple signal classification algorithms. This data includes not only the millimeter-wave measurement gap value and micro-Doppler velocity at each preset circumferential position, but also the received signal power, received channel noise power, and multipath component amplitude spectrum information used for confidence assessment. The traditional displacement monitoring data is acquired by an eddy current displacement sensor pre-installed at a specific position on the stator of the magnetic levitation bearing. This sensor senses changes in rotor surface position based on the principle of electromagnetic induction, outputs an analog voltage signal, and converts it into a digital displacement reading via an analog-to-digital converter. These two heterogeneous data sources are complementary in physical characteristics and error distribution, providing rich information redundancy for subsequent high-precision fusion. Because a single sensor cannot guarantee continuous measurement reliability under complex operating conditions, millimeter-wave radar sensors, while offering advantages such as non-contact operation, resistance to oil mist, and high precision, suffer from decreased measurement accuracy under strong multipath interference or deteriorating signal-to-noise ratio. Eddy current displacement sensors, while technologically mature and with fast response, have limitations including limited measurement range, susceptibility to electromagnetic interference, and sensitivity to rotor surface materials and temperature. By fusing measurement data from these two sensors, their respective advantages can be fully utilized and their shortcomings compensated for. When millimeter-wave measurement quality is good, the eddy current sensor is used as the primary source to obtain high precision and multi-point spatial information; when millimeter-wave measurement is affected by interference, the stable output of the eddy current sensor becomes more crucial, thus achieving highly reliable measurement under all operating conditions. The necessity of received signal power and received channel noise power lies in the fact that these two parameters directly reflect the quality status of millimeter-wave measurements. The magnitude of signal power characterizes the sufficiency of the echo intensity, while the magnitude of noise power characterizes the level of electromagnetic interference in the measurement environment. The ratio of these two, i.e., the signal-to-noise ratio, is a core indicator for evaluating measurement confidence. The necessity of multipath component amplitude spectrum information lies in the fact that there are multiple reflecting surfaces inside the magnetic levitation bearing, such as the stator inner wall, rotor surface, and end cover. Millimeter wave signals will produce multipath propagation phenomena. The main path echo corresponds to the actual rotor surface distance, while the multipath echo will generate false peaks in the spectrum. By analyzing the amplitude distribution of multipath components, the severity of multipath interference can be assessed, thereby dynamically adjusting the confidence weight of millimeter wave measurement data.

[0065] In the confidence assessment of millimeter-wave measurements, the signal-to-noise ratio (SNR) at each measurement point is first calculated. The SNR of millimeter-wave measurements is calculated by multiplying 10 by a logarithmic function to base 10, where the independent variable is the ratio of signal power to noise power. Signal power represents the power of the received effective echo signal, and noise power represents the noise power in the receiving channel. The SNR directly reflects the quality status of millimeter-wave measurements. A high SNR indicates that the echo signal strength is much greater than the noise level, resulting in smaller random errors and higher reliability in the measurement results. Conversely, a low SNR indicates that the echo signal is overwhelmed by noise, increasing the uncertainty of the measurement results. In this case, the weight of this measurement data in the fusion process should be reduced. The necessity of using signal power and noise power as inputs lies in the fact that these two parameters can quantify the strength of the effective information and the level of interference background, respectively. The ratio between them determines the SNR characteristics of the measurement data.

[0066] Next, the multipath interference index is evaluated. The multipath interference index is calculated by summing the absolute amplitude values ​​of the multipath components from the second to the Mth multipath component, and dividing by the absolute amplitude value of the main path (i.e., the first multipath component), where M is the total number of multipath components. The multipath interference index is introduced because the internal space of a magnetic levitation bearing contains multiple metal reflective surfaces such as the stator inner wall, rotor surface, end cover, and coil frame. Millimeter-wave signals undergo multiple reflections and scattering between these reflective surfaces, forming multipath propagation. The main path echo corresponds to the direct reflection from the rotor surface, carrying the true gap information, while the multipath echo corresponds to the signal reaching the receiving antenna after multiple reflections. These multipath echoes generate false peaks in the distance-dimensional spectrum, interfering with the identification of the real target. The multipath interference index quantifies the relative intensity of multipath interference by calculating the ratio of the amplitudes of non-main path components to the main path components. A larger index indicates stronger multipath echoes, which can easily cause measurement errors. In this case, the confidence weight of the millimeter-wave measurement data should be reduced.

[0067] Then, the confidence weights for millimeter-wave measurements are defined. The calculation method for the millimeter-wave measurement confidence weights is to apply a sigmoid function to a comprehensive index, which equals the signal-to-noise ratio (SNR) term minus the multipath interference term. The SNR term is calculated by multiplying a preset SNR adjustment coefficient by the millimeter-wave SNR, and the multipath interference term is calculated by multiplying a preset multipath interference adjustment coefficient by the multipath interference index. The default value for the preset SNR adjustment coefficient is 0.1, determined based on the typical variation range of the millimeter-wave SNR to ensure that the millimeter-wave measurement confidence weights vary within a reasonable range. The default value for the preset multipath interference adjustment coefficient is 2.0, determined based on the degree of impact of multipath interference on measurement accuracy to appropriately reduce the weight when multipath interference is severe.

[0068] In the traditional sensor reliability assessment process, the detection saturation range of an eddy current displacement sensor is a closed interval from a preset lower saturation limit to a preset upper saturation limit. The default value of the preset lower saturation limit is -2 mm, and the default value of the preset upper saturation limit is +2 mm. These two default values ​​are determined based on the linear measurement range of the selected eddy current displacement sensor, typically taking the boundary value of the sensor's nominal linear range. The default value of the preset normal change threshold is 50 micrometers. This default value is determined based on the maximum rate of change of rotor displacement and the sampling period under normal operating conditions of the magnetic levitation bearing, ensuring the ability to distinguish between normal vibration and abnormal jumps. Because eddy current sensors operate based on the principle of electromagnetic induction, there is a non-linear relationship between their output voltage and the distance from the coil to the measured metal surface. They maintain good linearity only within a certain distance range. When the rotor displacement exceeds this linear range, the sensor output enters the saturation region. At this point, the output value can no longer accurately reflect the actual displacement. If the measurement data under saturation is used for fusion, it will introduce a large systematic error. Therefore, it is necessary to determine the validity of the measurement data by detecting whether the sensor output is within the saturation range. The necessity of detecting the output change rate is that when the eddy current displacement sensor is working normally, its output should change smoothly and continuously with the rotor displacement, and the output difference between adjacent sampling times should be within a reasonable range. If a sudden change occurs that exceeds the preset normal change threshold, it usually means that the sensor has been affected by abnormal factors such as electromagnetic interference, signal line failure, or uneven rotor surface material. At this time, the reliability of the measurement data decreases, and its weight in the fusion process needs to be reduced.

[0069] The sensor confidence weight is calculated in a piecewise manner. When the absolute value of the difference between the current output and the previous output of the eddy current displacement sensor is less than a preset normal change threshold, and the current output of the eddy current displacement sensor is not within the saturation range, the sensor confidence weight is set to 1. The previous time interval is the current time minus a preset time interval, which is equal to the preset sampling period. Otherwise, the sensor confidence weight is a natural exponential function with the product of the negative of a preset attenuation coefficient and the abnormal change amount as the exponent. The default value of the preset attenuation coefficient is 0.05. This default value is determined based on the desired weight attenuation rate, ensuring that the weight decreases rapidly when the abnormal change amount is large. The abnormal change amount is equal to the absolute value of the difference between the current output and the previous output of the eddy current displacement sensor minus the preset normal change threshold. This abnormal change amount represents the magnitude of change exceeding the preset normal change threshold. The reason for using a segmented approach to calculate sensor confidence weights is to distinguish between normal and abnormal sensor states. Under normal conditions, the sensor output is stable and reliable, and should be assigned the maximum weight of 1 to fully utilize its measurement information. However, under abnormal conditions, the sensor output is unreliable, and the weight should be dynamically reduced according to the degree of abnormality. The advantage of using a natural exponential function to describe weight decay under abnormal conditions is that this function has a smooth decay characteristic; the greater the abnormal change, the faster the weight decays. This effectively suppresses the impact of severely abnormal data on the fusion result. At the same time, the continuity of the exponential function ensures a smooth transition in weight adjustment, avoiding abrupt changes in the fusion result that may be caused by hard threshold decisions.

[0070] In the design of the Extended Kalman Filter (EKF), an EKF state vector is constructed. This state vector is a column vector with a preset state dimension, the default value of which is 4. This default value is determined by the number of rotor dynamic state variables to be estimated. The state vector includes a first state element, a second state element, a third state element, and a fourth state element. The first state element is the fusion gap, the second state element is the fusion velocity, the third state element is the fusion acceleration, and the fourth state element is the eccentricity angle. The fusion velocity represents the first derivative of the fusion gap with respect to time, the fusion acceleration represents the second derivative of the fusion gap with respect to time, and the eccentricity angle represents the eccentricity angle of the rotor attitude. The reason for using the EKF for multi-source data fusion is that this algorithm can estimate the system state in real time through a recursive method, even in the presence of measurement noise and process noise. Compared with simple weighted averaging or median filtering methods, the Kalman filter fully utilizes the dynamic model information and the statistical characteristics of the measurement data, and can obtain the optimal state estimate with the minimum mean square error. The advantage of the Extended Kalman Filter (EKF) over the Standard Kalman Filter (SKF) lies in its ability to handle nonlinear systems. Although the state transition equation in this embodiment is linear, the measurement equation involves the weighted fusion of multiple heterogeneous sensors, exhibiting certain nonlinear characteristics. The EKF can effectively handle this nonlinear relationship through local linearization. The necessity of using fusion gap, fusion velocity, fusion acceleration, and eccentricity angle as state vector elements lies in the fact that these four state variables comprehensively describe the rotor's position, motion, and attitude information. The fusion gap is the core variable of greatest concern to the control system. The fusion velocity and fusion acceleration reflect the changing trend and dynamic characteristics of the gap, which are crucial for predicting future states and achieving look-ahead compensation. The eccentricity angle describes the rotor's attitude deviation relative to the stator and is an important basis for identifying rotor imbalance and implementing active avoidance.

[0071] The state vector is updated through a state transition equation, which describes the evolution of the state vector from the current time step to the next. The state vector at the next time step is equal to the state transition matrix multiplied by the state vector at the current time step, plus the process noise vector. The state transition matrix is ​​a matrix with a preset state dimension of rows and columns. When the preset state dimension is 4, the elements in the first row are 1, the preset sampling period, the square of the preset sampling period divided by 2, and 0. The elements in the second row are 0, 1, the preset sampling period, and 0. The elements in the third row are 0, 0, 1, and 0. The elements in the fourth row are 0, 0, 0, and 1. The default value for the preset sampling period is 0.1 milliseconds. This default value is determined based on the required response speed and data processing capability of the control system, ensuring that the dynamic characteristics of the rotor are captured while meeting real-time requirements. The process noise vector represents the uncertainty of the system dynamics model. The first row of the state transition matrix describes the evolution of the fusion gap, indicating that the fusion gap at the next time step is equal to the fusion gap at the current time step plus the product of the fusion velocity and the preset sampling period, plus the product of the fusion acceleration and half the square of the preset sampling period. This is the expression of the displacement formula for uniformly accelerated motion in a discrete-time system. The second row describes the evolution of the fusion velocity, indicating that the fusion velocity at the next time step is equal to the fusion velocity at the current time step plus the product of the fusion acceleration and the preset sampling period. This is the velocity formula for uniformly accelerated motion. The third row indicates that the fusion acceleration remains constant between adjacent time steps. This is a simplification assumption applicable to scenarios where acceleration changes relatively slowly. The fourth row indicates that the eccentricity angle remains constant between adjacent time steps because the eccentricity angle is mainly determined by the unbalanced mass distribution of the rotor and will not change significantly in a short time. The introduction of the process noise vector is necessary because real-world systems contain modeling errors, unknown disturbances, and random fluctuations during the state transition process. The process noise vector is used to describe the impact of these uncertainties on the state evolution. The Kalman filter quantifies the statistical characteristics of this uncertainty through the process noise covariance matrix and compensates for it during the state estimation process.

[0072] In the weighted fusion update process, a measurement vector is first constructed. The measurement vector is a column vector with a preset measurement dimension, the default value of which is 3. This default value is determined based on the number of measurement sources participating in the fusion. The measurement vector includes a first element, a second element, and a third element. The first element is the millimeter-wave measurement confidence weight multiplied by the millimeter-wave average gap; the second element is the sensor confidence weight multiplied by the eddy current displacement sensor output; and the third element is the average micro-Doppler velocity. The millimeter-wave average gap is calculated by summing the millimeter-wave measurement gap values ​​at a preset number of circumferential positions and dividing by the preset number. Specifically, it is calculated by summing all millimeter-wave measurement gap values ​​at circumferential position i from 1 to the preset number and then dividing by the preset number. The average micro-Doppler velocity is calculated by summing the micro-Doppler velocity values ​​at a preset number of circumferential positions and then dividing by the preset number. The first and second elements of the measurement vector correspond to the gap information measured by the millimeter-wave sensor and the eddy current sensor, respectively. The advantage of weighting by confidence level is that it dynamically adjusts the contribution of each sensor to the fusion result based on the real-time measurement quality. When a sensor's measurement quality is good, its weight is automatically increased; when its measurement quality is poor, its weight is automatically decreased, thus achieving adaptive data fusion. Calculating the average millimeter-wave gap is necessary because the millimeter-wave gap values ​​at a preset number of circumferential positions reflect the circumferential gap distribution of the rotor. Averaging these measurements yields the overall average gap of the rotor. This average value can suppress random errors and local anomalies from single-point measurements, providing a more stable and reliable gap estimate. The necessity of the average micro-Doppler velocity as the third element of the measurement vector lies in the fact that this velocity information provides a motion state observation independent of position measurements. It can be cross-validated with the fused velocity estimated by the Kalman filter, improving the robustness of the state estimation.

[0073] The measurement vector is updated through a measurement equation, which describes the relationship between the measurement vector and the state vector. The measurement vector at time t equals the measurement matrix multiplied by the state vector at time t, plus the measurement noise vector. The measurement matrix is ​​a matrix with preset measurement dimensions and preset state dimensions. When the preset measurement dimension is 3 and the preset state dimension is 4, the elements in the first row are 1, 0, 0, 0; the elements in the second row are 1, 0, 0, 0; and the elements in the third row are 0, 1, 0, 0. The measurement noise vector represents the uncertainty in the measurement process. The first and second rows of the measurement matrix are both 1, 0, 0, 0, indicating that the first and second elements of the measurement vector are observations of the first element of the state vector, i.e., the fusion gap. This reflects that both millimeter-wave measurement and eddy current sensor measurement are for rotor gap measurement, providing redundant observations of the same physical quantity. The third row of the measurement matrix is ​​0, 1, 0, 0, indicating that the third element of the measurement vector, i.e., the average micro-Doppler velocity, is an observation of the second element of the state vector, i.e., the fusion velocity, providing independent velocity measurement information. The introduction of measurement noise vectors is necessary because all sensor measurement processes involve random errors, including electronic noise, quantization errors, environmental interference, and other factors. Measurement noise vectors are used to describe the statistical characteristics of these random errors. Kalman filters quantify the noise level of each measurement channel by measuring the noise covariance matrix and rationally allocate the weights of each measurement data during the state update process.

[0074] The state estimate is updated using the Kalman gain matrix. After the update, high-confidence fusion gap and fusion velocity estimates are obtained, with the fusion gap estimate equal to the fusion gap and the fusion velocity estimate equal to the fusion velocity. The Kalman gain matrix is ​​a core parameter of the Kalman filter. This matrix is ​​calculated in real-time based on the process noise covariance matrix, measurement noise covariance matrix, and state estimation error covariance matrix. Its function is to determine the degree of correction of the state estimate by the measurement data. When the measurement noise is low, the Kalman gain matrix is ​​large, and the state estimate relies more on the measurement data for correction. When the measurement noise is high, the Kalman gain matrix is ​​small, and the state estimate relies more on the model prediction. Through this adaptive gain adjustment mechanism, the Kalman filter can achieve an optimal balance between model prediction and measurement update, obtaining the optimal state estimate in the sense of minimum mean square error. The fusion gap estimate and fusion velocity estimate, as outputs of step 200, provide high-quality input data for subsequent spatiotemporal consistency verification, harmonic extraction, predictive compensation, and critical approach warning.

[0075] Step 300: Based on the millimeter-wave data and fused data, perform spatiotemporal consistency verification and anomaly detection, and adjust the confidence weight according to the detection results.

[0076] In the spatial distribution consistency analysis, the standard deviation of the circumferential clearance is first calculated. The method for calculating the standard deviation of the circumferential clearance is to divide the squared clearance term by a preset number and then take the square root. The squared clearance term is the sum of the clearance difference terms at a preset number of circumferential positions, where the clearance difference term represents the difference between the millimeter-wave measured clearance and the average millimeter-wave clearance. The specific calculation steps are as follows: For each circumferential position i from 1 to a preset number, calculate the difference between the millimeter-wave measured clearance value at that position and the average millimeter-wave clearance, square this difference, sum the squared values ​​for all preset number of circumferential positions, divide the sum by the preset number, and finally take the square root of the quotient.

[0077] Next, the minimum gap position is determined. The minimum gap position angle is the circumferential angle position that minimizes the millimeter-wave measurement gap value. Specifically, it is determined by finding the angle corresponding to the position that minimizes the millimeter-wave measurement gap value among all circumferential positions from 1 to a preset number.

[0078] During the spatiotemporal consistency detection process, preset spatial consistency thresholds and preset gap warning thresholds are set. The default value of the preset spatial consistency threshold is 25 micrometers. This default value is determined based on the measurement accuracy of the millimeter-wave radar sensor and the circumferential gap uniformity during normal rotor operation, and is used to identify spatial distribution anomalies such as rotor eccentricity or stator deformation. The default value of the preset gap warning threshold is 150 micrometers. This default value is determined based on the safety clearance margin of the magnetic levitation bearing and the response capability of the control system, allowing sufficient warning time. When the standard deviation of the circumferential gap is greater than the preset spatial consistency threshold, or when the minimum millimeter-wave measured gap value among a preset number of circumferential positions is less than the preset gap warning threshold, a spatial anomaly flag is triggered. The spatial anomaly flag is set to 1 when triggered. The reason for setting dual criteria to trigger the spatial anomaly flag is that it is necessary to simultaneously monitor the uniformity of the gap distribution and the safety of the absolute value of the gap. The standard deviation of the circumferential gap reflects the degree of eccentricity of the rotor attitude, while the minimum millimeter-wave measured gap value is directly related to the risk of rubbing. The two criteria complement each other and can more comprehensively identify spatial anomalies.

[0079] During velocity consistency cross-validation, the velocity residual between the fused velocity and the micro-Doppler velocity is calculated. The velocity residual is calculated by subtracting the average micro-Doppler velocity from the estimated fused velocity and taking the absolute value. A preset velocity consistency threshold is set, with a default value of 0.3 mm / s. This default value is determined based on the accuracy of millimeter-wave Doppler velocimetry and the velocity estimation accuracy of the Kalman filter, and is used to detect the temporal consistency of the measurement data. When the velocity residual exceeds the preset velocity consistency threshold, a velocity inconsistency flag is triggered. The velocity inconsistency flag is set to 1 when triggered.

[0080] When a spatial anomaly flag or a velocity inconsistency flag is triggered, adaptive weight adjustment is performed. During this process, the source of the anomaly is first determined, and then the sensor weights are adjusted according to the anomaly type. If only the spatial anomaly flag is triggered without the velocity inconsistency flag, it indicates an abnormal dispersion in the circumferential gap distribution or an excessively small absolute gap value, but the velocity measurement remains consistent. In this case, the anomaly mainly stems from electromagnetic interference affecting the eddy current displacement sensor or its entry into the nonlinear saturation region. However, the multi-point spatial resolution and velocity measurement capabilities of millimeter-wave measurements remain reliable. Therefore, the sensor confidence weight is updated to the original sensor confidence weight multiplied by a preset sensor weight attenuation factor, and the millimeter-wave measurement confidence weight is updated to the original millimeter-wave measurement confidence weight multiplied by a preset millimeter-wave weight enhancement factor. The default value of the preset sensor weight attenuation factor is 0.1, determined based on the degree of reliability reduction of the eddy current sensor under electromagnetic interference or saturation conditions. The default value of the preset millimeter-wave weight enhancement factor is 1.2, determined based on the relative advantages of millimeter-wave measurements under such abnormal conditions.

[0081] When the velocity inconsistency flag is triggered, regardless of whether the spatial anomaly flag is triggered simultaneously, it indicates a significant deviation between the fused velocity estimate and the average micro-Doppler velocity. This anomaly primarily stems from multipath interference or signal-to-noise ratio degradation in millimeter-wave measurements, leading to micro-Doppler velocity measurement distortion, or insufficient accuracy of the millimeter-wave array measurement resulting in errors in the reconstruction of the gap spatial distribution field. Since single-point measurements by eddy current displacement sensors are often more stable and reliable under such conditions, the millimeter-wave measurement confidence weight is updated by multiplying the original millimeter-wave measurement confidence weight by a preset millimeter-wave weight attenuation factor. The sensor confidence weight is updated by multiplying the original sensor confidence weight by a preset sensor weight enhancement factor. The default value of the preset millimeter-wave weight attenuation factor is 0.5, determined based on the degree of reliability reduction in millimeter-wave measurements under multipath interference or signal-to-noise ratio degradation. The default value of the preset sensor weight enhancement factor is 1.5, determined based on the relative advantage of eddy current sensors under such abnormal conditions. After adjustment, all millimeter-wave measurement confidence weights and sensor confidence weights are normalized to ensure the sum of the weights is 1. The reason for performing adaptive weight adjustment is that different types of anomaly flags indicate different measurement failure modes. Spatial anomaly flags mainly reflect measurement problems of eddy current sensors, while velocity inconsistency flags mainly reflect a degradation in the quality of millimeter-wave measurements. By adjusting the weights according to the anomaly type, more reliable measurement sources under the current operating conditions can be accurately identified and automatically switched, thereby ensuring the accuracy and robustness of the fusion results under various abnormal operating conditions. The necessity of normalization is to ensure that the sum of the weights of each measurement source is always 1. This is a basic requirement of the weighted fusion algorithm and also guarantees the physical meaning and numerical stability of the fusion gap estimate.

[0082] Step 400: Perform adaptive harmonic extraction and frequency conversion synchronization notch filtering on the fused data to extract the synchronization disturbance component, specifically as follows: Figure 3 As shown.

[0083] In the design of the adaptive notch filter bank, a notch filter bank is constructed based on the real-time rotation frequency, comprising multiple notch filters. The real-time rotation frequency is obtained through encoder measurement and represents the actual rotation frequency of the rotor at the current moment, ranging from 0 to a preset maximum rotation frequency. The default value of the preset maximum rotation frequency is 1000 Hz, determined based on the maximum operating speed of the rotor supported by the magnetic levitation bearing. The k-th notch frequency is calculated by multiplying k by the real-time rotation frequency and dividing by twice pi, where k ranges from 1 to a preset harmonic order. A number of valley notch filters with a preset harmonic order are designed, with a default value of 5. This default value is determined based on the spectral distribution characteristics of the rotor unbalance force and electromagnetic force harmonic components, covering the fundamental frequency to the 5th harmonic. The reason for using an adaptive notch filter bank is that the magnetic levitation bearing rotor generates periodic disturbances synchronized with the rotation frequency during rotation. These disturbances mainly originate from rotor mass imbalance, geometric eccentricity, and electromagnetic force harmonic components. These synchronous disturbances cause periodic fluctuations in the clearance measurement value, interfering with the control system's judgment of the true clearance state. By designing a notch filter bank, harmonic components synchronized with the rotor frequency can be selectively suppressed, while asynchronous random disturbances and true gap variation information can be extracted. The adaptive characteristic is reflected in the dynamic adjustment of the notch frequency according to the real-time rotor frequency. When the rotor speed changes, the notch frequency automatically tracks the fundamental frequency and all harmonic frequencies of the rotor, ensuring that the notch filter always operates at the correct frequency. This adaptive capability enables the system to maintain good disturbance suppression performance over a wide speed range. The necessity of using the real-time rotor frequency as input lies in the fact that this parameter determines the notch frequency setting. Only by accurately acquiring the real-time rotor frequency can the notch filter be precisely aligned with the harmonic components that need to be suppressed.

[0084] An infinite impulse response (IIR) notch filter is used. The numerator of the transfer function of the k-th notch filter is 1 minus the cosine term plus z to the power of -2. The cosine term is twice the cosine function multiplied by z to the power of -1. The independent variable of the cosine function is twice pi multiplied by the k-th notch frequency and then multiplied by the preset sampling period. The transfer function is given in the Z-domain, where z is a complex variable of the Z-transform, and z to the power of -1 represents the unit delay of one sampling period. The denominator of the transfer function is 1 minus the notch cosine term plus the notch term. The notch cosine term is twice the preset notch coefficient multiplied by the cosine function multiplied by z to the power of -1. The notch term is the square of the preset notch coefficient multiplied by z to the power of -2. The independent variable of the cosine function is the same as in the numerator. The default value of the preset notch coefficient is 0.95. This default value is determined based on the trade-off between the required notch depth and bandwidth, and is used to control the notch bandwidth. A larger preset notch coefficient produces a narrower notch bandwidth and a deeper notch depth. The advantages of using an infinite impulse response (IRR) notch filter lie in its simple structure, low computational cost, and good real-time performance. It requires only a small amount of historical data and coefficient storage to achieve efficient frequency-selective filtering, making it particularly suitable for magnetic levitation bearing control systems that require real-time processing. The transfer function design of the IRR notch filter results in extremely deep attenuation at the notch frequency, effectively suppressing this frequency component while having minimal impact on other frequency components. This high selectivity ensures that while suppressing synchronization disturbances, it does not over-filter out useful signal components. The selection of the preset notch coefficient reflects the trade-off between notch depth and bandwidth. A larger preset notch coefficient, such as 0.95, produces a very narrow notch bandwidth and a very deep notch depth, accurately suppressing the target frequency while minimizing the impact on adjacent frequencies. However, this also requires higher precision in setting the notch frequency. If the preset notch coefficient is too small, the notch bandwidth will be too wide, potentially filtering out too many useful signal components.

[0085] During the extraction of synchronization disturbance components, a notch filter bank is applied to the fused gap estimate to extract periodic disturbances synchronized with the rotation frequency. The synchronization disturbance component is composed of the superposition of harmonic components from the fundamental frequency to a preset harmonic order. The fundamental frequency is the real-time rotor rotation frequency, with k equal to 1. Each harmonic component is equal to the amplitude of the k-th harmonic multiplied by a sine function, where the independent variable of the sine function is the product of k times the real-time rotation frequency and time t, plus the phase of the k-th harmonic. A least mean square adaptive algorithm is used to estimate the amplitude and phase of the k-th harmonic in real time. The default value of the preset step size parameter of this algorithm is determined based on a trade-off between convergence speed and steady-state error. The necessity of extracting the synchronization disturbance component lies in the fact that although these periodic disturbances synchronized with the rotation frequency cause fluctuations in the gap measurement, they are essentially predictable deterministic signals. By accurately extracting these synchronization disturbance components, corresponding compensating forces can be applied in the control commands for active cancellation, thereby significantly reducing the synchronous vibration amplitude of the rotor and improving the gap control accuracy. The advantage of using the least mean square adaptive algorithm lies in its ability to estimate the amplitude and phase of each harmonic online in real time, without the need for offline identification or table lookup. It possesses strong adaptive capability; when the rotor's imbalance state or electromagnetic force characteristics change, the algorithm can automatically track these changes and update the harmonic parameter estimates. The least mean square algorithm has low computational complexity, involving only simple multiplication and addition operations, making it easy to implement in real-time control systems.

[0086] In the process of identifying asynchronous disturbances, the asynchronous disturbance component is calculated. The asynchronous disturbance component equals the estimated fusion gap value minus the synchronous disturbance component, and is used to identify random shocks and aperiodic interference. In addition to periodic disturbances synchronized with the rotation frequency, magnetic levitation bearing systems are also subject to various asynchronous random disturbances, such as airflow pulsations, mechanical shocks, and sudden load changes. These asynchronous disturbances are unpredictable and cannot be suppressed by simple harmonic compensation. However, by separating them from the estimated fusion gap value, important information can be provided for fault diagnosis and anomaly detection. When the amplitude or frequency characteristics of the asynchronous disturbance component show abnormal changes, it may indicate a potential fault or an increase in external interference, requiring corresponding countermeasures.

[0087] Step 500: Based on the fused data and synchronization disturbance components, perform predictive state observation and look-ahead compensation, generate compensation force, and generate overall control commands, specifically as follows: Figure 4 As shown.

[0088] In the design process of the rotor dynamics state observer, a rotor dynamics state space model is established. The dynamics state space model is as follows: the derivative of the rotor dynamics state vector with respect to time is equal to the sum of the product of the rotor dynamics system matrix and the rotor dynamics state vector, the product of the rotor dynamics input matrix and the electromagnetic force command, and the product of the rotor dynamics disturbance matrix and the external disturbance. The rotor dynamics state vector is a preset dynamics state dimension column vector, with a default value of 3, determined according to the order of the rotor dynamics model. The three elements are, in order, the radial clearance of the rotor relative to the stator, the first derivative of this clearance with respect to time, and the second derivative of this clearance with respect to time. The radial clearance of the rotor relative to the stator refers to the real-time radial distance between the outer surface of the rotor and the inner surface of the stator. This clearance value is calculated in real-time using a Luneburger observer. The calculation process uses the fused clearance estimate and fused velocity estimate output in step 200 as observation values. The preset observer gain matrix is ​​used to correct the state predicted by the model, thereby obtaining a smooth clearance state quantity with predictive characteristics. External disturbances refer to unknown external forces acting on the rotor, including rotor imbalance forces, airflow disturbance forces, mechanical impact forces, and other disturbances that cannot be directly measured. In the dynamic state-space model, these external disturbances act on the time derivative of the rotor dynamic state vector through the rotor dynamic disturbance matrix. The rotor dynamic state-space model can describe the rotor's motion under electromagnetic forces and external disturbances, providing a theoretical basis for state observation and predictive control. The state-space model uses a system of first-order differential equations, facilitating the application of observer design and optimal control methods from modern control theory. The advantage of explicitly incorporating external disturbances into the state-space model lies in the ability to estimate external disturbances in real time through an extended observer, thereby achieving feedforward compensation based on disturbance estimation. This active disturbance rejection strategy, compared to traditional feedback control, can more quickly suppress the impact of disturbances on the system, improving control performance.

[0089] The rotor dynamics system matrix is ​​a matrix with preset dynamic state dimensions in rows and columns. When the preset dynamic state dimension is 3, the elements in the first row are 0, 1, 0, the elements in the second row are 0, 0, 1, and the elements in the third row are negative preset stiffness coefficient divided by preset rotor mass, negative preset damping coefficient divided by preset rotor mass, and negative preset quadratic damping coefficient divided by preset rotor mass. The default value for the preset rotor mass is 5 kg, determined based on the actual measured rotor mass. The default value for the preset stiffness coefficient is 2 x 10^6 Newtons per meter, determined based on the magnetic circuit design and electromagnetic stiffness characteristics of the magnetic levitation bearing. The default value for the preset damping coefficient is 1000 N / s / m, determined based on the damping characteristics obtained from system identification experiments. The default value for the preset quadratic damping coefficient is 50 N / s² / m, obtained based on the nonlinear damping effect identification during high-speed operation. The rotor dynamics input matrix is ​​a column vector with one row and one column, representing the preset dynamic state dimension. Its three elements are 0, 0, and 1 divided by the preset rotor mass, respectively. This matrix describes the effect of electromagnetic force commands on the time derivative of the rotor dynamics state vector. The rotor dynamics disturbance matrix is ​​also a column vector with one row and one column, representing the preset dynamic state dimension. Its three elements are 0, 0, and 1 divided by the preset rotor mass, respectively. This matrix describes the effect of external disturbances on the time derivative of the rotor dynamics state vector. Its structure is the same as the rotor dynamics input matrix, indicating that both external disturbances and electromagnetic force commands act on the rotor's acceleration term. The reason for introducing a preset quadratic damping coefficient in the third row of the rotor dynamics system matrix is ​​that when a magnetically levitated bearing rotor rotates at high speed, in addition to the linear damping effect, there is also a nonlinear damping effect proportional to the square of the velocity. This nonlinear damping mainly originates from airflow resistance and eddy current losses. By introducing a preset quadratic damping coefficient into the dynamic model, the actual motion characteristics of the rotor can be described more accurately, improving the accuracy of state observation and prediction. The necessity of using preset rotor mass, preset stiffness coefficient, preset damping coefficient, and preset secondary damping coefficient as inputs lies in the fact that these parameters fully describe the dynamic characteristics of the rotor and are the basis for establishing an accurate dynamic model. The values ​​of these parameters need to be obtained through actual measurement or system identification experiments to ensure the consistency between the model and the actual system.

[0090] During the state estimation process, a Luneburger observer is used to estimate the current rotor dynamics state estimate and external disturbance estimate based on the fusion gap estimate and fusion velocity estimate output in step 200. The state estimation equation of the Luneburger observer is as follows: the derivative of the current rotor dynamics state estimate with respect to time is equal to the sum of the rotor dynamics system matrix multiplied by the current rotor dynamics state estimate, the rotor dynamics input matrix multiplied by the electromagnetic force command, the rotor dynamics disturbance matrix multiplied by the external disturbance estimate, and the preset observer gain matrix multiplied by the observation error. The observation error is equal to the fusion gap estimate minus the first element of the current rotor dynamics state estimate. This observation error reflects the deviation between the actual measured value and the observer estimate. This deviation is fed back into the state estimation equation through the preset observer gain matrix to achieve real-time correction of the current rotor dynamics state estimate. The external disturbance estimate is a real-time estimate of the external disturbance, obtained through an extended Luneburger observer. The extended observer treats the external disturbance as an augmented state variable, assuming that the derivative of the external disturbance with respect to time is zero or changes slowly. The external disturbance is tracked and estimated through dynamic adjustment of the observer. The time derivative of the external disturbance estimate is equal to the preset disturbance observer gain multiplied by the observation error. The default value of the preset disturbance observer gain is determined using the extended observer pole placement method. This default value is determined based on the desired convergence speed and noise immunity of the disturbance estimate, and is typically chosen to be 0.1 to 0.5 times the largest element in the preset observer gain matrix to ensure the smoothness and stability of the disturbance estimate. The advantage of using a Luneburger observer lies in its simple structure and high computational efficiency. By introducing an observation error feedback mechanism, it can use actual measurements to correct the model's predicted state in real time, thereby overcoming the influence of model errors and process noise, and obtaining a more accurate state estimate than simple model prediction. The convergence characteristics of the Luneburger observer can be adjusted by the preset observer gain matrix. A well-designed gain matrix can ensure that the observation error converges quickly to zero, allowing the state estimate to rapidly track the true state. The advantage of extending the Luneburger observer to treat external disturbances as augmented states is that it can infer the unknown disturbance force acting on the rotor by observing measurable states such as gaps and speeds without directly measuring the disturbance. This disturbance estimation capability provides key information for realizing active disturbance rejection control. By introducing negative external disturbance estimates into the control commands, the influence of disturbances on the rotor can be actively offset, significantly improving the system's disturbance rejection performance.

[0091] The preset observer gain matrix is ​​designed using the pole placement method to ensure dynamic stability of the estimation error. The preset observer gain matrix is ​​a column vector with one row and one column representing a preset dynamic state dimension. The default value of this column vector is obtained by placing the closed-loop poles of the observer at preset pole locations. The default values ​​of the preset pole locations are determined based on the desired observer response speed and noise immunity, and are typically chosen to be 2 to 5 times the real part of the open-loop poles of the rotor dynamics system. When the preset dynamic state dimension is three, the three elements of the preset observer gain matrix correspond to the gap observation gain, velocity observation gain, and acceleration observation gain, respectively. These three gain values ​​are obtained by solving the pole placement equation, ensuring that the characteristic polynomial roots of the observer are located at the preset pole locations. The advantage of using the pole placement method to design the preset observer gain matrix is ​​that this method can directly specify the dynamic characteristics of the observer, and by selecting appropriate pole locations, an optimal balance can be achieved between the observer response speed and noise immunity. The reason for placing the closed-loop poles of the observer at positions 2 to 5 times the real part of the open-loop poles of the rotor dynamics system is that such a configuration ensures that the convergence speed of the observer is faster than the dynamic response speed of the system itself, so that the state estimation can track the changes in the system state in a timely manner. At the same time, the pole position should not be too far to the left, otherwise the observer will be too sensitive to measurement noise, causing drastic fluctuations in the state estimation.

[0092] During the predictive calculation process, state prediction is performed using a preset prediction time step. The default value for the preset prediction time step is 1 millisecond. This default value is determined based on the execution delay of the control system and the desired look-ahead compensation effect, and is usually set to the sum of the controller calculation delay, communication delay, and actuator response delay. The predicted rotor dynamics state vector is equal to the product of the first matrix exponential function and the current rotor dynamics state estimate, plus the first integral term, plus the second integral term. The first integral term is the integral from 0 to the preset prediction time step, and its integrand is the second matrix exponential function multiplied by the rotor dynamics input matrix and then multiplied by the electromagnetic force command. The second integral term is the integral from 0 to the preset prediction time step, and its integrand is the third matrix exponential function multiplied by the rotor dynamics disturbance matrix and then multiplied by the external disturbance estimate. The base of the first matrix exponential function is a natural constant, and the exponent is the rotor dynamics system matrix multiplied by the preset prediction time step. The base of the second matrix exponential function is a natural constant, and the exponent is the rotor dynamics system matrix multiplied by the integral variable of the first integral term. The base of the third matrix exponential function is a natural constant, and the exponent is the rotor dynamics system matrix multiplied by the integral variable of the second integral term. The predicted rotor dynamic state vector is a function of time t plus a preset prediction time step. This prediction process not only considers the natural evolution of the current rotor dynamic state estimate and the effect of electromagnetic force commands, but also incorporates the influence of external disturbance estimates into the prediction calculation through a second integral term, thereby improving prediction accuracy and enabling the predicted rotor dynamic state vector to more accurately reflect the actual motion state of the rotor at future moments. The necessity of state prediction lies in the inherent execution delay in the magnetic levitation bearing control system, including the controller's calculation time, communication transmission time, and the response time of the electromagnetic actuator. These delays mean that by the time the control command is applied to the rotor, the actual state of the rotor has already changed. If the control command is calculated based on the current state, control lag will occur, affecting the control effect and even leading to system instability. By predicting the rotor state after a preset prediction time step and calculating the control command based on the predicted state, the impact of execution delay can be effectively compensated, achieving look-ahead control. The advantage of using matrix exponential functions and integral forms for state prediction lies in the fact that this method provides an exact solution for the state transition of a linear time-invariant system. Compared to simple numerical integration methods such as the Euler method or Runge-Kutta method, the matrix exponential method achieves higher prediction accuracy, especially when the prediction time step is large or the system dynamics are fast. The necessity of incorporating external disturbance estimates into the prediction calculation is that external disturbances continuously act on the rotor. If the influence of disturbances is ignored in the prediction process, the predicted state will deviate from the true state. By considering external disturbance estimates in the prediction calculation, the future trajectory of the rotor under the influence of disturbances can be predicted more accurately.

[0093] In the forward-looking compensation force calculation process, the compensation force is calculated based on the predicted clearance and the preset safety threshold. The predicted clearance is the first element of the predicted rotor dynamics state vector. The default value of the preset safety threshold is 200 micrometers, which is determined based on the protective clearance of the magnetic levitation bearing and the adjustment margin of the control system to ensure sufficient safe distance between the rotor and stator under various operating conditions. The compensation force is calculated by multiplying the difference between the preset safety threshold and the predicted clearance by the preset position gain, and then adding a velocity term, which is the preset velocity gain multiplied by the predicted velocity. The default value of the preset position gain is 5 x 10⁴ Newtons per meter, which is determined based on the desired position adjustment stiffness and system stability margin. The default value of the preset velocity gain, i.e., the preset damping coefficient, is 100 Newton-seconds per meter, which is determined based on the desired damping ratio and vibration suppression effect. The predicted velocity is the second element of the predicted rotor dynamics state vector. The advantage of calculating the compensation force based on predicted gap is that this method achieves true look-ahead control. The compensation force is not for the current gap state, but for the predicted gap state after a preset prediction time step. This look-ahead characteristic can effectively compensate for the execution delay of the control system, ensuring that the rotor state matches the design state of the control command when it is actually applied to the rotor, thereby improving the timeliness and accuracy of control. The compensation force adopts a proportional-derivative form, where the proportional term generates the restoring force based on the deviation between the predicted gap and the preset safety threshold, and the derivative term generates the damping force based on the predicted velocity. This proportional-derivative control structure can achieve rapid position adjustment and effective vibration suppression while ensuring system stability.

[0094] In the control command synthesis process, the total control command equals the output of the traditional PID controller plus the compensation force, the notch compensation term, and the disturbance feedforward compensation term. The notch compensation term equals the negative value of the preset notch gain multiplied by the synchronous disturbance component, and the disturbance feedforward compensation term equals the negative value of the external disturbance estimate. The default value of the preset notch gain is 1 x 10⁵ Newtons per meter. This default value is determined based on the amplitude range of the synchronous disturbance component and the desired disturbance suppression effect, and the optimal value is obtained through frequency domain analysis and experimental debugging. The disturbance feedforward compensation term actively cancels out external disturbances using the external disturbance estimate. Since the external disturbance estimate is a real-time estimate of the unknown external force acting on the rotor, introducing a negative external disturbance estimate into the total control command generates a compensation force equal in magnitude and opposite in direction to the external disturbance. This achieves feedforward compensation before or during the disturbance's action on the rotor, significantly reducing the impact of external disturbances on rotor clearance and improving the disturbance resistance and control accuracy of the magnetic levitation bearing system. This combination of feedforward compensation based on external disturbance estimates and feedback control of a traditional PID controller forms a composite control strategy of feedback plus feedforward. Feedback control is responsible for eliminating steady-state errors and tracking the setpoint, while feedforward compensation actively suppresses observable disturbances. The two work together to achieve high-precision control of the rotor clearance. The advantage of using multiple compensation terms to synthesize the overall control command lies in the fact that each compensation term targets different error sources and disturbance types. The traditional PID controller output provides basic feedback adjustment, the compensation force achieves look-ahead control based on predicted states, the notch compensation term specifically suppresses periodic disturbances synchronized with the rotational frequency, and the disturbance feedforward compensation term actively cancels unknown external disturbances. This multi-level compensation mechanism can comprehensively address various control challenges faced by magnetic levitation bearing systems, achieving superior overall performance compared to single control methods. The notch compensation term uses negative feedback because the synchronization disturbance component represents a periodic fluctuation in the clearance measurement synchronized with the rotational frequency. By applying the compensation force corresponding to the negative synchronization disturbance component, a control action opposite in phase to the synchronization disturbance can be generated, thereby canceling the influence of the synchronization disturbance on the rotor position and reducing the synchronous vibration amplitude of the rotor.

[0095] Step 600: During the clearance monitoring process, the system monitors the minimum clearance and the predicted minimum clearance in real time. The minimum clearance is calculated by taking the minimum value of the clearance measured by the millimeter wave from 1 to a preset number of circumferential positions i. The predicted minimum clearance is calculated by combining the eccentricity and eccentricity angle obtained in step 100 with the predicted rotor dynamic state vector to calculate the predicted clearance value at each circumferential position at the predicted time, and then taking the minimum predicted clearance value among all circumferential positions. The minimum clearance position angle is the circumferential angle position that makes the millimeter wave measured clearance value reach the minimum value; this angle position corresponds to the orientation where the rotor and stator are closest. The necessity of simultaneously monitoring the minimum clearance and the predicted minimum clearance lies in the fact that the minimum clearance reflects the minimum safety margin between the rotor and stator at the current moment and is a direct indicator for judging the risk of collision and rubbing, while the predicted minimum clearance reflects the minimum safety margin at future moments and can provide early warning of potential collision and rubbing hazards. By combining the dual monitoring of the current state and the predicted state, the system can achieve more timely and reliable safety protection.

[0096] In the design of the multi-level early warning mechanism, the system sets three levels of early warning thresholds to achieve hierarchical safety protection. The default value for the first-level early warning threshold is 250 micrometers, determined based on the gap fluctuation range during normal operation, and is used for early warning and status monitoring. The default value for the second-level early warning threshold is 180 micrometers, determined based on the critical condition for initiating proactive avoidance measures, and is used to activate proactive avoidance control. The default value for the third-level early warning threshold is 120 micrometers, determined based on the minimum safe gap required for emergency protection activation, and is used to trigger the emergency protection mechanism. These three thresholds satisfy the relationship that the third-level early warning threshold is less than the second-level early warning threshold, which is less than the first-level early warning threshold, which is less than the preset safety threshold, forming a progressive safety protection system. The advantage of using a multi-level early warning mechanism is that it can take corresponding measures of varying strength according to the degree of gap proximity, avoiding overreaction or underreaction. The first-level early warning provides an early warning signal, allowing sufficient response time for operators and the control system; the second-level early warning activates proactive avoidance control, actively increasing the gap by applying avoidance force; and the third-level early warning triggers the emergency protection mechanism, using the strongest control force to quickly stop the collision tendency. This tiered response strategy ensures both security and avoids unnecessary interference with the normal operation of the system due to strong controls.

[0097] During the Level 1 warning response, when the minimum gap is less than the preset Level 1 warning threshold but greater than or equal to the preset Level 2 warning threshold, the system enters the Level 1 warning state. At this time, the system increases the gap monitoring frequency, shortening the sweep period of the millimeter-wave radar sensor from the preset sweep period to a preset accelerated sweep period. The default value of the preset accelerated sweep period is 0.5 milliseconds, which is half the default value of the preset sweep period, used to improve the measurement update rate to more densely track gap changes. Simultaneously, the system activates the trend analysis function to calculate the gap change rate. The gap change rate is calculated by subtracting the minimum gap from the minimum gap at the previous moment from the current minimum gap, and then dividing the difference by the preset time interval. The system records the Level 1 warning event and sends a warning signal to the control system. Based on the warning signal, the control system appropriately increases the preset position gain and preset velocity gain by multiplying the original gain value by the preset Level 1 gain adjustment coefficient. The default value of the preset Level 1 gain adjustment coefficient is 1.1, determined according to the required control response enhancement level under Level 1 warning conditions, ensuring increased control stiffness without causing system oscillations. Furthermore, the system continuously monitors and predicts the minimum clearance during the Level 1 warning state. If the predicted minimum clearance shows a continued decreasing trend and the absolute value of the clearance change rate exceeds the preset clearance change rate threshold, it prepares to enter the Level 2 warning state in advance. The default value of the preset clearance change rate threshold is 0.5 micrometers per millisecond, which is determined based on the rotor's dynamic characteristics and the warning response time. Increasing the clearance monitoring frequency is necessary because the Level 1 warning state indicates that the system has deviated from its normal operating state, and the clearance is approaching the safety boundary. At this point, more intensive measurement data is needed to accurately track the dynamic changes in the clearance and promptly detect the trend of clearance deterioration. By shortening the frequency sweep cycle, the measurement update rate can be doubled, enabling the control system to acquire clearance information and respond more promptly. The advantage of activating the trend analysis function is that the clearance change rate reflects the trend and speed of clearance change. When the clearance change rate is negative and the absolute value is large, it indicates that the clearance is decreasing rapidly, posing a high risk of rubbing, requiring stronger countermeasures in advance. By monitoring the clearance change rate, the system can achieve a smooth transition from Level 1 to Level 2 warnings, avoiding lag in warning level switching. The reason for appropriately increasing the control gain is that the adjustment capability of the control system needs to be enhanced under the first-level warning state. By increasing the preset position gain, the influence of position error on the control force can be amplified and the position adjustment speed can be accelerated. By increasing the preset speed gain, the damping effect can be enhanced and the oscillation and overshoot of the gap can be suppressed. However, the gain adjustment range should not be too large, otherwise it may cause the system to become oversensitive and generate oscillation. The preset first-level gain adjustment coefficient is set to 1.1, which reflects the balance between enhancing control performance and system stability.

[0098] During the Level 2 warning response, when the minimum clearance is less than the preset Level 2 warning threshold but greater than or equal to the preset Level 3 warning threshold, the system enters the Level 2 warning state and initiates an active avoidance strategy. At this time, the system calculates the minimum clearance direction angle, which indicates the closest orientation between the rotor and stator. The system applies a radial offset avoidance force to actively increase the minimum clearance. The avoidance force is calculated by multiplying the negative preset avoidance force gain by the natural exponential function and then by the radial unit vector, where the exponent of the natural exponential function is the negative clearance term, and the clearance term is the difference between the minimum clearance and the preset Level 3 warning threshold divided by the preset attenuation length. The default value for the preset avoidance force gain is 500 Newtons, determined based on the rotor mass and desired avoidance acceleration to ensure the avoidance force effectively changes the rotor position without causing excessive dynamic impact. The default value for the preset attenuation length is 50 micrometers, determined based on the desired attenuation rate of the avoidance force as the clearance changes, causing the avoidance force to increase rapidly when approaching the preset Level 3 warning threshold and gradually decrease when moving away from it. The radial unit vector is the unit vector at the angle of the minimum clearance position, pointing radially from the stator geometric center to the rotor geometric center. The avoidance force is applied in the opposite direction, i.e., pushing the rotor away from the position closest to the stator inner wall. In the secondary warning state, the system simultaneously enhances the compensation force in step 500, further increasing the preset position gain and preset speed gain by multiplying the original gain value by the preset secondary gain adjustment coefficient. The default value of the preset secondary gain adjustment coefficient is 1.3, which is determined based on the stronger control response required in the secondary warning state. The system continuously monitors the avoidance effect. If the minimum clearance fails to rise above the preset secondary warning threshold within the preset avoidance response time after applying the avoidance force, or if the predicted minimum clearance shows that the clearance will continue to decrease to below the preset tertiary warning threshold, the system prepares to enter the tertiary warning state. The default value of the preset avoidance response time is 10 milliseconds, which is determined based on the rotor dynamic response time and the effect of the avoidance force. The necessity of activating the active avoidance strategy lies in the fact that when the clearance decreases below the preset secondary warning threshold, conventional feedback control and look-ahead compensation alone may not be able to prevent the clearance from decreasing further in time. Additional avoidance force is required to actively push the rotor away from the danger zone. The advantage of using an exponentially decaying avoidance force is that this form can automatically adjust the magnitude of the avoidance force according to how close the minimum clearance is to the preset tertiary warning threshold. When the minimum clearance is close to the preset tertiary warning threshold, the exponential function value is large, the avoidance force is strong, and it can generate a large thrust. When the minimum clearance is far from the preset tertiary warning threshold, the exponential function value is small, the avoidance force is weak, and excessive intervention in normal control is avoided. This adaptive avoidance force adjustment mechanism ensures both safety and maintains control smoothness.

[0099] During the Level 3 early warning response, when the minimum clearance is less than the preset Level 3 early warning threshold, or when the predicted minimum clearance is less than the preset Level 3 early warning threshold, the system enters the Level 3 early warning state and activates the emergency protection mechanism. At this time, the system determines that there is a risk of contact between the rotor and stator, requiring the highest level of protection measures. The system cuts off the conventional position error channel, i.e., suspends the output of the traditional PID controller in step 500, to avoid the delayed response of the conventional controller from delaying the emergency protection opportunity. The system injects an emergency suppression signal based on millimeter-wave real-time measurement, which directly acts on the electromagnetic force actuator to achieve the fastest response. The calculation method for the emergency suppression signal, i.e., the emergency control force, is: the preset emergency position gain multiplied by the difference between the preset safety threshold and the millimeter-wave measured clearance value at the minimum clearance position, minus the product of the preset emergency damping gain and the micro-Doppler velocity at the minimum clearance position. The millimeter-wave measured clearance value at the minimum clearance position refers to the millimeter-wave measured clearance value measured at the circumferential position corresponding to the angle of the minimum clearance position, and the micro-Doppler velocity at the minimum clearance position refers to the micro-Doppler velocity value measured at that circumferential position. The default value for the preset emergency position gain is 2 x 10⁵ Newtons per meter. This default value is determined based on the maximum control force and rapid response requirements under emergency conditions, and is typically set to 2 to 5 times the normal preset position gain to generate sufficient braking force. The default value for the preset emergency damping gain is 500 Newton-seconds per meter. This default value is determined based on the damping requirements and system stability during emergency braking, providing strong damping to quickly suppress rotor movement towards the stator. The direction of the emergency control force is from the minimum clearance position to the opposite direction of the stator's geometric center, i.e., pushing the rotor rapidly away from the rubbing danger zone. In the third-level warning state, the system simultaneously triggers an alarm signal, notifying the operator to take manual intervention measures, such as reducing the speed or activating the backup protection system. The system continuously monitors the minimum clearance and the predicted minimum clearance. If the minimum clearance recovers to above the preset level 3 warning threshold within the preset emergency response time after applying emergency control force, and the predicted minimum clearance indicates that the clearance will continue to increase, the system will gradually exit the level 3 warning state, restore the normal control channel, and determine whether to maintain the level 2 or level 1 warning state based on the value of the minimum clearance. The default value of the preset emergency response time is 5 milliseconds, which is determined based on the timeliness requirements of emergency protection and the rotor dynamic response characteristics. If the minimum clearance fails to recover effectively within the preset emergency response time, the system will maintain the emergency protection mechanism until the clearance returns to a safe level or the system shutdown protection procedure is initiated. The necessity of initiating the emergency protection mechanism lies in the fact that when the clearance decreases below the preset level 3 warning threshold, the safety margin between the rotor and stator is already very small. If no emergency measures are taken, the rotor may rub against the stator in a very short time, causing serious mechanical damage or even system failure.The reason for cutting off the conventional position error channel is that traditional PID controllers rely on position error feedback for adjustment, which introduces a certain response delay. In emergencies, this delay can lead to untimely protection measures. By cutting off the conventional channel and directly injecting emergency control force based on real-time millimeter-wave measurement, the fastest response can be achieved, minimizing the time required for protection action. The advantage of calculating the emergency control force directly based on the millimeter-wave measured gap value and micro-Doppler velocity at the minimum gap position is that these two measurements are real-time status information at that dangerous position. Without complex filtering and fusion processing, they can reflect the rotor's true state at that position with minimal delay, thus achieving the most timely emergency control. The reason for setting the preset emergency position gain and preset emergency damping gain to several times the normal gain is that in emergency conditions, a braking force much greater than the normal control force is needed to effectively change the rotor's motion trend in a very short time and prevent rubbing.

[0100] The embodiments of the present invention have been described above. However, the embodiments are not limited to the specific implementation methods described above. The specific implementation methods described above are merely illustrative and not restrictive. Those skilled in the art can make more equivalent embodiments under the guidance of the present embodiments, and all of them are within the protection scope of the present embodiments.

Claims

1. A millimeter-wave measurement device and compensation method for the clearance of magnetic levitation bearings, characterized in that, include: Millimeter-wave data is obtained by multi-point gap measurement and velocity extraction using a millimeter-wave radar array deployed on the inner wall of the stator of the magnetic levitation bearing. Collect multi-source data, calculate confidence weights based on the multi-source data and millimeter-wave data, and perform confidence assessment and fusion based on the confidence weights to obtain fused data; Based on the millimeter-wave data and fused data, spatiotemporal consistency verification and anomaly detection are performed, and the confidence weights are adjusted according to the detection results. Adaptive harmonic extraction and frequency conversion synchronization notch filtering are performed on the fused data to extract synchronization disturbance components; Based on the fused data and synchronous disturbance components, predictive state observation and look-ahead compensation are performed to generate compensation force and generate overall control command.

2. The millimeter-wave measurement device and compensation method for magnetic levitation bearing clearance according to claim 1, characterized in that, The multi-point gap measurement and velocity extraction include: Millimeter-wave radar sensors are deployed at a predetermined number of equally spaced circumferential positions on the inner wall of the magnetic levitation bearing stator to form a multi-input multi-output measurement array; Each radar sensor transmits a linear frequency modulated continuous wave signal, and the received echo signal is mixed with the local oscillator signal to obtain the intermediate frequency signal. Millimeter-wave data includes millimeter-wave measurement gap and micro-Doppler velocity. Based on the intermediate frequency signal, the beat frequency is extracted through Fourier transform, the millimeter-wave measurement gap is calculated based on the beat frequency, and the micro-Doppler velocity is extracted through phase change. The millimeter-wave measurement gap is calculated by multiplying the preset speed of light by the preset sweep period by the beat frequency and then dividing by twice the preset sweep bandwidth. The micro-Doppler velocity is calculated by multiplying the preset millimeter-wave wavelength by the phase difference between adjacent sweep periods, dividing by four times pi and then dividing by the preset sweep period.

3. The millimeter-wave measurement device and compensation method for magnetic levitation bearing clearance according to claim 1, characterized in that, The confidence assessment and fusion include: The multi-source data includes signal power, noise power, and eddy current displacement sensor output; the confidence weights include millimeter-wave measurement confidence weights and sensor confidence weights; and the fused data includes gap estimates and fusion velocity estimates. The signal-to-noise ratio and multipath interference index are calculated based on the signal power and noise power, and the confidence weight of millimeter-wave measurement is calculated through the signal-to-noise ratio and multipath interference index. Based on the output and saturation range of the eddy current displacement sensor, the confidence weight of the sensor is calculated. Based on the millimeter-wave measurement confidence weights and sensor confidence weights, an extended Kalman filter state vector and measurement vector are constructed, and a weighted fusion update is performed to obtain the fusion gap estimate and fusion speed estimate.

4. The millimeter-wave measurement device and compensation method for magnetic levitation bearing clearance according to claim 3, characterized in that, Weighted fusion updates include: The measurement vector includes a first element, a second element, and a third element. The first element is the millimeter-wave measurement confidence weight multiplied by the millimeter-wave average gap. The second element is the sensor confidence weight multiplied by the eddy current displacement sensor output. The third element is the average micro-Doppler velocity. The state vector includes a first state element, a second state element, a third state element, and a fourth state element. The first state element is the fusion gap, the second state element is the fusion velocity, the third state element is the fusion acceleration, and the fourth state element is the eccentricity angle. The fusion velocity represents the first derivative of the fusion gap with respect to time, the fusion acceleration represents the second derivative of the fusion gap with respect to time, and the eccentricity angle represents the eccentricity angle of the rotor attitude. The state vector is updated by the state transition equation. The state vector at the next time step is equal to the state transition matrix multiplied by the state vector at the current time step, plus the process noise vector. The measurement vector is updated by the measurement equation. The measurement vector at time t is equal to the measurement matrix multiplied by the state vector at time t, plus the measurement noise vector. The state estimate is updated by the Kalman gain matrix to obtain the fusion gap estimate and the fusion velocity estimate, where the fusion gap estimate is equal to the fusion gap and the fusion velocity estimate is equal to the fusion velocity.

5. The millimeter-wave measurement device and compensation method for magnetic levitation bearing clearance according to claim 3, characterized in that, The signal-to-noise ratio (SNR) is calculated by multiplying 10 by a logarithmic function to the base 10, where the independent variable of the logarithmic function is the ratio of signal power to noise power. The multipath interference index is calculated by summing the absolute amplitude values ​​of the multipath components from the second to the Mth multipath component and dividing by the absolute amplitude value of the first multipath component, where M is the total number of multipath components. The millimeter-wave measurement confidence weight is calculated by applying a Sigmoid function to the comprehensive index, where the comprehensive index is the SNR term minus the multipath interference term. The SNR term is calculated by multiplying a preset SNR adjustment coefficient by the millimeter-wave SNR, and the multipath interference term is calculated by multiplying a preset multipath interference adjustment coefficient by the multipath interference index. The saturation range is a closed interval from the preset lower saturation limit to the preset upper saturation limit. When the absolute value of the difference between the output of the eddy current displacement sensor at the current moment and the output at the previous moment is less than the preset normal change threshold, and the output of the eddy current displacement sensor at the current moment is not within the saturation range, the sensor confidence weight is set to 1, where the previous moment is the current moment minus the preset time interval; otherwise, the sensor confidence weight is a natural exponential function with the product of the negative value of the preset attenuation coefficient and the abnormal change amount as the exponent, where the abnormal change amount is the absolute value of the difference between the output of the eddy current displacement sensor at the current moment and the output at the previous moment minus the preset normal change threshold.

6. The millimeter-wave measurement device and compensation method for magnetic levitation bearing clearance according to claim 1, characterized in that, The spatiotemporal consistency verification and anomaly detection include: The standard deviation of the circumferential clearance is calculated. When the standard deviation of the circumferential clearance is greater than a preset spatial consistency threshold, or when the smallest millimeter-wave measurement clearance among a preset number of circumferential positions is less than a preset clearance warning threshold, a spatial anomaly flag is triggered. The standard deviation of the circumferential clearance is the square root of the clearance square term divided by a preset number. The clearance square term is the sum of the clearance difference terms at the preset number of circumferential positions. The clearance difference term is the difference between the millimeter-wave measurement clearance and the average millimeter-wave clearance. Calculate the velocity residual. When the velocity residual is greater than the preset velocity consistency threshold, trigger the velocity inconsistency flag. The velocity residual is the absolute value of the fused velocity estimate minus the average microDoppler velocity. When a spatial anomaly flag or a velocity inconsistency flag is triggered, the confidence weights are adjusted and normalization is performed.

7. The millimeter-wave measurement device and compensation method for magnetic levitation bearing clearance according to claim 1, characterized in that, Adjusting the confidence weights includes: When only the spatial anomaly flag is triggered and the velocity inconsistency flag is not triggered, the sensor confidence weight is updated to the original sensor confidence weight multiplied by the preset sensor weight attenuation factor, and the millimeter wave measurement confidence weight is updated to the original millimeter wave measurement confidence weight multiplied by the preset millimeter wave weight enhancement factor. When the speed inconsistency flag is triggered, regardless of whether the spatial anomaly flag is triggered at the same time, the millimeter-wave measurement confidence weight is updated to the original millimeter-wave measurement confidence weight multiplied by the preset millimeter-wave weight attenuation factor, and the sensor confidence weight is updated to the original sensor confidence weight multiplied by the preset sensor weight enhancement factor.

8. The millimeter-wave measurement device and compensation method for magnetic levitation bearing clearance according to claim 1, characterized in that, The adaptive harmonic extraction and frequency-switching synchronous notch filter includes: A notch filter bank is constructed based on real-time frequency conversion, and the notch filter bank includes multiple notch filters. A notch filter bank is applied to the estimated fusion gap to extract the synchronization disturbance component. The synchronization disturbance component is composed of the superposition of harmonic components from the fundamental frequency to the preset harmonic order. Each harmonic component is the amplitude of the kth harmonic multiplied by a sine function. The independent variable of the sine function is the product of k times the real-time frequency and time t plus the phase of the kth harmonic. The amplitude and phase of the kth harmonic are estimated in real time using the least mean square adaptive algorithm.

9. The millimeter-wave measurement device and compensation method for magnetic levitation bearing clearance according to claim 1, characterized in that, The predictive state observation and look-ahead compensation include: A rotor dynamics state-space model is established. The dynamics state-space model is as follows: the derivative of the rotor dynamics state vector with respect to time is equal to the sum of the product of the rotor dynamics system matrix and the rotor dynamics state vector, the product of the rotor dynamics input matrix and the electromagnetic force command, and the product of the rotor dynamics disturbance matrix and the external disturbance. Using the Luneburger observer, the current rotor dynamic state estimate and external disturbance estimate are estimated based on the fusion gap estimate and fusion velocity estimate; Perform state prediction at a preset prediction time step to obtain the predicted rotor dynamics state vector; The compensation force is calculated based on the predicted gap and the preset safety threshold. The predicted gap is the first element of the predicted rotor dynamic state vector. Generate a general control command, which is the sum of the output of the traditional PID controller, the compensation force, the notch compensation term, and the disturbance feedforward compensation term. The notch compensation term is the negative value of the preset notch gain multiplied by the synchronous disturbance component, and the disturbance feedforward compensation term is the negative value of the external disturbance estimate.

10. The millimeter-wave measurement device and compensation method for magnetic levitation bearing clearance according to claim 9, characterized in that, The predicted rotor dynamics state vector is the product of the first matrix exponential function and the current rotor dynamics state estimate, plus the first integral term, plus the second integral term. The first integral term is the integral from 0 to the preset prediction time step, and its integrand is the second matrix exponential function multiplied by the rotor dynamics input matrix and then multiplied by the electromagnetic force command. The second integral term is the integral from 0 to the preset prediction time step, and its integrand is the third matrix exponential function multiplied by the rotor dynamics disturbance matrix and then multiplied by the external disturbance estimate. The base of the first matrix exponential function is a natural constant, and the exponent is the rotor dynamics system matrix multiplied by the preset prediction time step. The base of the second matrix exponential function is a natural constant, and the exponent is the integral variable of the rotor dynamics system matrix multiplied by the first integral term. The base of the third matrix exponential function is a natural constant, and the exponent is the integral variable of the rotor dynamics system matrix multiplied by the second integral term. The compensation force is the difference between the preset safety threshold and the predicted gap multiplied by the preset position gain, plus the velocity term, which is the preset velocity gain multiplied by the predicted velocity, where the predicted velocity is the second element of the predicted rotor dynamic state vector.

Images