A Novel Data-Driven Feature Extraction Strategy and Its Application in Looseness Detection of Rotor-Bearing System

: During the service of rotating machinery, the pedestal bolts are prone to looseness due to the vibration environment, which affects the performance of rotating machinery and generate potential safety hazard. To monitor the occurrence and deterioration of the pedestal looseness in time, this paper proposes a data-driven diagnosis strategy for the rotor-bearing system. Firstly, the ﬁnite element model of a rotor-bearing system is established, which considers the piecewise nonlinear force caused by pedestal looseness. Taking the vibration response as output and periodic unbalanced force as input, the system’s NARX (Nonlinear Auto-Regressive with exogenous inputs) model is identiﬁed. Then GALEs (Generalized Associated Linear Equations) are used to evaluate NOFRFs (Nonlinear Output Frequency Response Functions) of the NARX model. Based on the ﬁrst three-order NOFRFs, the analytical expression of the second-order optimal weighted contribution rate is derived and used as a new health indicator. The simulation results show that compared with the conventional NOFRFs-based health indicator, the new indicator is more sensitive to weak looseness. Finally, a rotor-bearing test rig was built, and the pedestal looseness was performed. The experiment results show that as the degree of looseness increases, the new indicator SRm shows a monotonous upward trend, increasing from 0.48 in no looseness to 0.90 in severe looseness, a rise of 89.7%. However, the traditional indicator Fe2 has no monotonicity, which further veriﬁes the sensitivity of the ﬁrst three-order NOFRFs-based second-order optimal weighted contribution rate and the effectiveness of the proposed data-driven feature extraction strategy.


Introduction
Affected by installation, vibration environment, and structural deformation, looseness is prone to occur between the bearing seat and foundation (or casing) in rotating machinery [1]. When the looseness fault becomes serious, the vibration of rotating parts will exceed the allowable value, which is likely to cause other coupling faults, such as rub impact [2,3]. Therefore, early detection of looseness is essential to the service life and safe operation of rotating machinery.
The study on the fault feature extraction of looseness in rotor-bearing systems can be divided into three categories. The first is the feature extraction based on the physical model and failure mechanism. By exploring the vibration characteristics of the rotor-bearing system caused by looseness, the corresponding fault diagnosis methods are proposed. Jiang [4] proposed to linearize the nonlinear model of a rotor-bearing-pedestal system and use the difference between the linearized model and the nonlinear model as an indicator to measure the degree of nonlinearity caused by pedestal looseness. Zhang [5] considered nonlinear support and loss foundation when modelling a rotor-bearing system and analyzed the influence of the two kinds of faults on the transient and steady-state response

Data-Driven Feature Extraction
This paper mainly uses the NARX model to conduct the data-driven modeling of a rotor-bearing system. In engineering practice, the NARX model can represent a large class of nonlinear systems. The structure of a discrete single-input and single-output NARX model is [21]: where M and K are integers, k represents discrete time, p + q = m and ∑ K k 1 ,k p+q =1 = ∑ K k 1 =1 · · · ∑ K k p+q =1 ; c p,q k 1 , · · · , k p+q are the coefficients of the model terms. y(k) and u(k) represent the discrete output and input, respectively. After the NARX model identification based on the real-time or historical vibration data, the Model Predicted Output (MPO) method is used to verify the model's effectiveness. Then the GALEs method mentioned in Ref. [20] is used to decompose the NARX model and obtain the analytical expression of the system's first N-order output y n (k). Aiming at the general NARX model structure represented by Equation (1), the specific decomposition of GALEs is c p,q k 1 , · · · , k p+q y K n−q,p (k) c p,0 k 1 , · · · , k p y K n,p (k) (2) where n = 1, . . . , N, K = k 1 , . . . , k p+q represents the integer set of delay corresponding to coefficients c p,q k 1 , · · · , k p+q .
y i k − k p y K n−i,p−1 (k) y K n,1 (k) = y n (k − k 1 ) The NARX model that satisfies the MPO criterion can be expanded to any order. At this time, the NOFRFs of the NARX model can be obtained from the GALEs as [22] G n (jω) = Y n (jω) U n (jω) (4) where the nth (n = 1, . . . , N) order output spectra Y n (jω) and the input spectra U n (jω) can be obtained from the Discrete Time Fourier Transform (DTFT) of y n (k) and u n (k), respectively. At this time, the output spectrum of the nonlinear system can be expressed as G n (jω)U n (jω) (5) Previously, identifying NOFRFs under harmonic excitation required two excitations with different excitation intensities but the same excitation frequency combining the least square method, which limits the application of NOFRFs in the fault diagnosis of rotating machinery. For the GALEs method, the NOFRFs can be obtained only by the NARX model obtained from one harmonic excitation, which significantly facilitates the application of NOFRFs in engineering. NOFRFs can directly characterize the nonlinearity of a system, but NOFRFs are multi-valued functions, and some orders of NOFRFs may be not sensitive and effective. Therefore, it is necessary to further extract sensitive features from the obtained NOFRFs.
For NOFRFs, the high-order G n (jω) (n > 1) can reflect the nonlinearity of a system properly, while G 1 (jω) can only reflect the linear characteristics. In addition, the higher the order, the smaller the G n (jω). However, the sensitive features of faults are sometimes included in high-order G n (jω). For this reason, this paper proposes a weighted contribution rate of NOFRFs to amplify the weak features in high-order G n (jω). To improve the contribution rate of the high-order G n (jω) (n > 1), the G n (jω) is subjected to feature weighting. The specific calculation method is as where, T n (jω) is the weighted G n (jω), n ρ is the weighting factor, and ρ is an indefinite constant, which can be selected based on different faults. From Equation (6), the weighted contribution rate proposed has the following characteristics.
The weighting method can amplify the high-order G n (jω) and increase its proportion to the contribution rate. Based on this, the weighted contribution rate of NOFRFs (Rn) has been proposed as [23,24] In combining contribution rate and feature weighting with NOFRFs, the order n is introduced into calculating the weighted contribution rate. With the increase of the order n, the higher-order G n (jω) (n > 1) take up more and more weight. It solves the problem that the magnitude of high-order G n (jω) is small by amplifying the contribution rate of high-order G n (jω).
When ρ is 0, the weighted contribution rate of NOFRFs (Rn) proposed is the NOFRFsbased indicator Fe in [25]. Therefore, the weight contribution rate of NOFRFs can be regarded as a generalized expression of the indicator Fe. Moreover, it can be seen from Equation (8), when ρ is greater than 0, the value of the weighting factor is greater than 1, and the higher the order, the larger the corresponding weighting factor. It means the contribution rate of the high-order G n (jω) (n > 1) reduces rather than increases. Therefore, ρ = 0 is not an optimal choice.
Only when ρ is less than 0, the weighting factor n ρ is less than 1. It means that G n (jω) (n > 1) is multiplied by a multiplication factor greater than 1, and the higher the order, the larger the multiplication factor. In other word, the contribution of higher-order G n (jω) (n > 1) is larger. Therefore, the contribution rate of the G n (jω) will be amplified when ρ is only less than 0. Next, it needs to find an optimal ρ value. This paper mainly determines the optimal ρ value based on the second-order weighted contribution rate. According to Equation (8), the expression of Rn2 is as Mathematics 2023, 11, 2769 5 of 18 where, Rn2(ρ) is only a function of ρ. To simplify the process of solving, denote ϕ i = +∞ −∞ |G i (jω)|dω,ϕ i is only related to the nonlinearity of the system. When only the first three orders NOFRFs are considered, the function Rn2(ρ) can be rewrite as Calculate the derivative of the function Rn2(ρ) about ρ, as When the function dRn2(ρ)/dρ is equal 0, there is only one critical point ρ m Thus, the analytical expression of the maximum of the function Rn2(ρ) is SRm is defined as the second-order optimal weighted contribution rate based on the first three orders of NOFRFs, and the derived NOFRFs-based feature will be used as a health indicator to characterize the weak faults. The data-driven sensitive feature extraction process of the rotor system proposed is shown in Figure 1.
When the function dRn2(ρ)/dρ is equal 0, there is only one critical point ρm Thus, the analytical expression of the maximum of the function Rn2(ρ) is SRm is defined as the second-order optimal weighted contribution rate based on the first three orders of NOFRFs, and the derived NOFRFs-based feature will be used as a health indicator to characterize the weak faults. The data-driven sensitive feature extraction process of the rotor system proposed is shown in Figure 1. The proposed algorithm was performed by using the software MATLAB R2022b on a laptop equipped with an Intel Core i7 processor. The implementation of the proposed datadriven feature extraction strategy in the bearing-rotor system in Figure 1 is as follows.
Step 1: Acquire the vibration displacement signal of the shaft under the health state and operating state of a rotor system form simulation data or experiment data. For the acquired experiment data, the TSA (Time synchronous Averaging) is used to improve signal smoothness and signal-to-noise ratio.
Step 2: Establish the NARX models characterising the rotor system of health state The proposed algorithm was performed by using the software MATLAB R2022b on a laptop equipped with an Intel Core i7 processor. The implementation of the proposed data-driven feature extraction strategy in the bearing-rotor system in Figure 1 is as follows. Step 1: Acquire the vibration displacement signal of the shaft under the health state and operating state of a rotor system form simulation data or experiment data. For the acquired experiment data, the TSA (Time synchronous Averaging) is used to improve signal smoothness and signal-to-noise ratio.
Step 2: Establish the NARX models characterising the rotor system of health state and operating state using the downsampled vibration signal. Use the FROLS (Forward Regression with Orthogonal Least Squares) algorithm to determine the model terms and their coefficients and apply the stability and accuracy criteria to validate the model.
Step 3: Extract the sensitive feature from the identified NARX models. Introduce the unbalanced excitation to calculate the models' NOFRFs (G n (jω), n = 1, . . . , N) using the GALEs in Equation (2). Then, obtain the derived NOFRFs-based feature SRm as a health indicator to characterize the health condition of the rotor system.

Dynamic Modelling of a Rotor-Bearing System with Pedestal Looseness
To verify the proposed fault feature extraction method, a rotor-bearing system with pedestal looseness is considered. First, based on a rotor-bearing test rig, considering the coupling structure, the finite element model of the rotor-bearing system is established, as shown in Figure 2. The specific parameters of the finite element model of the rotor-bearing system are shown in Table 1. The rotor structural parameters in Table 1 are consistent with the real dimensions of the test rig, while the supporting stiffness and damping are given by reference to Ref. [26] and combining the measured rotor vibration responses in the experiment. GALEs in Equation (2). Then, obtain the derived NOFRFs-based feature SRm as a health indicator to characterize the health condition of the rotor system.

Dynamic Modelling of a Rotor-Bearing System with Pedestal Looseness
To verify the proposed fault feature extraction method, a rotor-bearing system with pedestal looseness is considered. First, based on a rotor-bearing test rig, considering the coupling structure, the finite element model of the rotor-bearing system is established, as shown in Figure 2. The specific parameters of the finite element model of the rotor-bearing system are shown in Table 1. The rotor structural parameters in Table 1 are consistent with the real dimensions of the test rig, while the supporting stiffness and damping are given by reference to Ref. [26] and combining the measured rotor vibration responses in the experiment.  Assume that the bearing seat at the non-driving end has a looseness fault, as shown in Figure 3. The nonlinearity caused by the looseness is simulated by a piecewise linear spring element [26,27], and the equivalent looseness stiffness is as

Parameters Values
Shaft diameter at left bearing/mm 55 Shaft diameter at disc/mm 55 The span between two bearings l/mm 800 Coupling cantilever length l c/ mm 165 The span between the right bearing and disc l/mm 164. Assume that the bearing seat at the non-driving end has a looseness fault, as shown in Figure 3. The nonlinearity caused by the looseness is simulated by a piecewise linear spring element [26,27], and the equivalent looseness stiffness is as where, k f1 is the tensile stiffness of the bolt, k f2 is the stiffness of the loosened bolt, k f3 is the stiffness when the bolt is not loosened, δ is the loosening gap, and z 3 is the vibration displacement of the mass point m 3 in the z direction. where, kf1 is the tensile stiffness of the bolt, kf2 is the stiffness of the loosened bolt, kf3 is the stiffness when the bolt is not loosened, δ is the loosening gap, and z3 is the vibration displacement of the mass point m3 in the z direction. Assuming that the stiffness of the unloosed bolt is approximately equal to the tensile stiffness,

Right pedestal
At this time, the governing equation of the rotor-bearing system can be expressed as M is the mass matrix considering the mass of the shaft, coupling and disc, D is the damping matrix including bearing damping and rotor gyro moment, K is the stiffness matrix, Q is the external force vector, and q is the displacement vector. Assume that the looseness occurs between the right bearing seat and the foundation, the looseness gap δ = 0.6 mm, the looseness mass m3 = 12 kg, and the foundation stiffness before looseness is kf3 = 1 × 10 9 N/m. Assume that the eccentricity only exists at the disc, and the amount of unbalanced m·r = 880 g·mm. The rotor speed is ω = 2400 r/min, and the sampling frequency is fs = 10,240 Hz. This paper only considers that a single pedestal is loose, and the looseness gap is much larger than the vibration amplitude of the pedestal mass. Therefore, the change in stiffness only occurs when the pedestal is in contact with the foundation. As the degree of pedestal looseness increases, the corresponding looseness stiffness gradually decreases. The simulation mainly considers four different looseness stiffness. The stiffness of the unloose bolt is kf3 = 1 × 10 9 N/m, and the stiffness of the loose bolt kf2 is 1 × 10 9 N/m, 1 × 10 8 N/m, 5 × 10 7 N/m, and 1 × 10 7 N/m. When kf2 = 1 × 10 9 N/m, the rotor-bearing system has no pedestal looseness fault.

Influence of Looseness Stiffness on Vibration Response
After obtaining the differential motion equation of the rotor-bearing system, since the equation is a typical nonlinear system, the Newmark-β method is used to obtain the numerical solution. This paper mainly analyzes the influence of the looseness stiffness of the right pedestal on the vibration response, including the time domain, frequency domain, and orbit of the mbr mass point. Assuming that the stiffness of the unloosed bolt is approximately equal to the tensile stiffness, At this time, the governing equation of the rotor-bearing system can be expressed as M is the mass matrix considering the mass of the shaft, coupling and disc, D is the damping matrix including bearing damping and rotor gyro moment, K is the stiffness matrix, Q is the external force vector, and q is the displacement vector. Assume that the looseness occurs between the right bearing seat and the foundation, the looseness gap δ = 0.6 mm, the looseness mass m 3 = 12 kg, and the foundation stiffness before looseness is k f3 = 1 × 10 9 N/m. Assume that the eccentricity only exists at the disc, and the amount of unbalanced m·r = 880 g·mm. The rotor speed is ω = 2400 r/min, and the sampling frequency is f s = 10,240 Hz. This paper only considers that a single pedestal is loose, and the looseness gap is much larger than the vibration amplitude of the pedestal mass. Therefore, the change in stiffness only occurs when the pedestal is in contact with the foundation. As the degree of pedestal looseness increases, the corresponding looseness stiffness gradually decreases. The simulation mainly considers four different looseness stiffness. The stiffness of the unloose bolt is k f3 = 1 × 10 9 N/m, and the stiffness of the loose bolt k f2 is 1 × 10 9 N/m,

Influence of Looseness Stiffness on Vibration Response
After obtaining the differential motion equation of the rotor-bearing system, since the equation is a typical nonlinear system, the Newmark-β method is used to obtain the numerical solution. This paper mainly analyzes the influence of the looseness stiffness of the right pedestal on the vibration response, including the time domain, frequency domain, and orbit of the m br mass point.
As the looseness stiffness decreases, the vibration amplitude of the rotating shaft gradually increases, as shown in Figure 4. The increase from 0.5 µm, when there is no looseness to the final 5 µm, is nearly ten times larger, and the increasing direction is only along the positive z-axis. This proves that the vibration amplitude is less than the looseness gap, and the pedestal can vibrate freely along the positive z-axis without restriction, which aligns with the assumption.
Mathematics 2023, 11, 2769 8 along the positive z-axis. This proves that the vibration amplitude is less than looseness gap, and the pedestal can vibrate freely along the positive z-axis with restriction, which aligns with the assumption.  From the orbits in Figure 5, as the degree of looseness increases, the orbit gradu evolves from a regular ellipse when the base is not loose to an irregular ellipse wi "local cusp" and restricted on one side. It shows that a serious unilateral collision occu in the system at this time. At the same time, from the spectrum in Figure 6, as the loose stiffness decreases, the frequency components gradually increase. For the healthy r system that does not have looseness, only rotational frequency is included in spectrum. When looseness occurs and deteriorates, namely, the looseness stiffnes gradually decreases from 1 × 10 9 N/m to 1 × 10 7 N/m, the amplitude of the dou rotational frequency in the spectrum increases from 0 µm to 2.4125 µm. The amplitud the double rotational frequency changes obviously, which indicates the degree nonlinearity of the system is increasing.  From the orbits in Figure 5, as the degree of looseness increases, the orbit gradually evolves from a regular ellipse when the base is not loose to an irregular ellipse with a "local cusp" and restricted on one side. It shows that a serious unilateral collision occurred in the system at this time. At the same time, from the spectrum in Figure 6, as the looseness stiffness decreases, the frequency components gradually increase. For the healthy rotor system that does not have looseness, only rotational frequency is included in the spectrum. When looseness occurs and deteriorates, namely, the looseness stiffness k f2 gradually decreases from 1 × 10 9 N/m to 1 × 10 7 N/m, the amplitude of the double rotational frequency in the spectrum increases from 0 µm to 2.4125 µm. The amplitude of the double rotational frequency changes obviously, which indicates the degree of nonlinearity of the system is increasing. stiffness decreases, the frequency components gradually increase. For the healthy rotor system that does not have looseness, only rotational frequency is included in the spectrum. When looseness occurs and deteriorates, namely, the looseness stiffness kf2 gradually decreases from 1 × 10 9 N/m to 1 × 10 7 N/m, the amplitude of the double rotational frequency in the spectrum increases from 0 µm to 2.4125 µm. The amplitude of the double rotational frequency changes obviously, which indicates the degree of nonlinearity of the system is increasing.

Data-Driven Fault Feature Extraction
To extract the fault feature of the pedestal looseness in the rotor-bearing system, first take the responses of mbr in z direction obtained from the numerical analysis as output, and take the synchronous excitation generated by the unbalanced force of the disc as input. Based on the FROLS method, the NARX model under different looseness stiffness is identified. Before identifying the NARX model, the initial model parameters must be determined through the trial-and-error method to make the final NARX model satisfies the MPO criterion. The number of NARX model terms is 14, and the maximum time lag of input and output is 6, respectively. The specific NARX model terms under four different looseness stiffness are shown in Table 2. The * in [*] represents the corresponding coefficient of the terms, and the 9th line represents the constant terms.
To ensure the accuracy of the identified NARX model, the MPO prediction output of the NARX model and the actual output obtained from the numerical simulation are compared and analyzed in the time domain and frequency domain, as shown from Figures  7-10. The predicted output of the NARX model is the same as the original simulation waveform, and the frequency domain characteristics of the first seven times the rotational frequency are consistent. It shows that the NARX model identified under each looseness stiffness has high accuracy and reliability. The effective NARX model lays the foundation for the subsequent NOFRF evaluation and further sensitive feature extraction.

Data-Driven Fault Feature Extraction
To extract the fault feature of the pedestal looseness in the rotor-bearing system, first take the responses of m br in z direction obtained from the numerical analysis as output, and take the synchronous excitation generated by the unbalanced force of the disc as input. Based on the FROLS method, the NARX model under different looseness stiffness is identified. Before identifying the NARX model, the initial model parameters must be determined through the trial-and-error method to make the final NARX model satisfies the MPO criterion. The number of NARX model terms is 14, and the maximum time lag of input and output is 6, respectively. The specific NARX model terms under four different looseness stiffness are shown in Table 2. The * in [*] represents the corresponding coefficient of the terms, and the 9th line represents the constant terms.
To ensure the accuracy of the identified NARX model, the MPO prediction output of the NARX model and the actual output obtained from the numerical simulation are compared and analyzed in the time domain and frequency domain, as shown from Figures 7-10. The predicted output of the NARX model is the same as the original simulation waveform, and the frequency domain characteristics of the first seven times the rotational frequency are consistent. It shows that the NARX model identified under each looseness stiffness has high accuracy and reliability. The effective NARX model lays the foundation for the subsequent NOFRF evaluation and further sensitive feature extraction.    After the identified NARX model has been decomposed by GALEs, the first threeorder output of the system can be obtained. Combined with the first three-order input of the system, the NOFRFs can be evaluated with Equation (4). Finally, the sensitive feature indicator SRm can be obtained according to Equation (13). The comparison between SRm under different looseness stiffness and the conventional NOFRFs-based indicator Fe2 is shown in Figure 11. It can be seen that when the looseness stiffness is greater than 1 × 10 8 N/m, the changing trend of the two indicators is the same, with the looseness stiffness decreasing. When the looseness stiffness is less than 1 × 10 8 N/m, the indicator SRm is significantly larger than Fe2, indicating that the sensitivity of the new indicator SRm to looseness stiffness is better than that of the conventional NOFRFs-based indicator Fe2. After the identified NARX model has been decomposed by GALEs, the first threeorder output of the system can be obtained. Combined with the first three-order input of the system, the NOFRFs can be evaluated with Equation (4). Finally, the sensitive feature indicator SRm can be obtained according to Equation (13). The comparison between SRm under different looseness stiffness and the conventional NOFRFs-based indicator Fe2 is shown in Figure 11. It can be seen that when the looseness stiffness is greater than 1 × 10 8 N/m, the changing trend of the two indicators is the same, with the looseness stiffness decreasing. When the looseness stiffness is less than 1 × 10 8 N/m, the indicator SRm is significantly larger than Fe2, indicating that the sensitivity of the new indicator SRm to looseness stiffness is better than that of the conventional NOFRFs-based indicator Fe2.

Experimental Verification
To verify the proposed data-driven feature extraction method and simulation results, a rotor-bearing test rig was used to analyze the influence of the pedestal looseness on the vibration response. The specific structure of the test rig is shown in Figure 12. The rotorbearing system is a single-disc two-fulcrum structure. The test rig considers the structural characteristics of the high-pressure rotor of the aero-engine. For the high-pressure rotor, the structure has the characteristic of mass distribution bias. In addition, the support bearing of the high-pressure rotor is a cylindrical roller bearing at one end and an angular contact ball bearing at the other. The models of the two types of bearings selected for this test rig are N209E and QJ210M, respectively. The coupling is a rope coupling, which can better prevent the vibration of the motor from being transmitted to the shaft. The bearing lubrication

Experimental Verification
To verify the proposed data-driven feature extraction method and simulation results, a rotor-bearing test rig was used to analyze the influence of the pedestal looseness on the vibration response. The specific structure of the test rig is shown in Figure 12. The rotorbearing system is a single-disc two-fulcrum structure. The test rig considers the structural characteristics of the high-pressure rotor of the aero-engine. For the high-pressure rotor, the structure has the characteristic of mass distribution bias. In addition, the support bearing of the high-pressure rotor is a cylindrical roller bearing at one end and an angular contact ball bearing at the other. The models of the two types of bearings selected for this test rig are N209E and QJ210M, respectively. The coupling is a rope coupling, which can better prevent the vibration of the motor from being transmitted to the shaft. The bearing lubrication method adopts oil-air lubrication. The pedestal at the non-driving end is loose. There are three different degrees of looseness: no looseness, slight looseness, and severe looseness.

Experimental Verification
To verify the proposed data-driven feature extraction method and simulation results, a rotor-bearing test rig was used to analyze the influence of the pedestal looseness on the vibration response. The specific structure of the test rig is shown in Figure 12. The rotorbearing system is a single-disc two-fulcrum structure. The test rig considers the structural characteristics of the high-pressure rotor of the aero-engine. For the high-pressure rotor, the structure has the characteristic of mass distribution bias. In addition, the support bearing of the high-pressure rotor is a cylindrical roller bearing at one end and an angular contact ball bearing at the other. The models of the two types of bearings selected for this test rig are N209E and QJ210M, respectively. The coupling is a rope coupling, which can better prevent the vibration of the motor from being transmitted to the shaft. The bearing lubrication method adopts oil-air lubrication. The pedestal at the non-driving end is loose. There are three different degrees of looseness: no looseness, slight looseness, and severe looseness. Laser sensor: 1pulse/rev Bearing-rotor system Measurement system Figure 12. Test rig structure and sensors installation location.
The measurement system mainly acquires the vibration displacement signal of the shaft, the acceleration signal of the bearing seat, and the rotational speed signal. The layout of the sensors and their detailed information are shown in the dashed box in Figure  12. The eddy current sensors (model RP660677) were used to measure the vibration displacement in the horizontal and vertical directions of the shaft, respectively. The accelerometers (model BH5031EX-050-VL) were used to measure the vertical vibration acceleration of the two bearing housings at the drive and non-drive ends, respectively. A Figure 12. Test rig structure and sensors installation location.
The measurement system mainly acquires the vibration displacement signal of the shaft, the acceleration signal of the bearing seat, and the rotational speed signal. The layout of the sensors and their detailed information are shown in the dashed box in Figure 12. The eddy current sensors (model RP660677) were used to measure the vibration displacement in the horizontal and vertical directions of the shaft, respectively. The accelerometers (model BH5031EX-050-VL) were used to measure the vertical vibration acceleration of the two bearing housings at the drive and non-drive ends, respectively. A photoelectric sensor (model OGP700) was used to measure the rotational speed of the shaft through the reflective sheet pasted on the shaft. The NI-9234 card and Cdaq-9178 chassis are used to collect and store the acquired data by the self-made Labview program. The motor speed is set to 2400 rpm, and the sampling frequency f s of the displacement and acceleration signals are both 5120 Hz.
The original horizontal vibration signal of the shaft near the non-driving end under no looseness, slight looseness, and severe looseness is shown in Figure 13. As the degree of looseness increases, the vibration amplitude decreases. The rotor-bearing test rig structure has a lot of mating gaps, and many factors affect vibration, such as assembly. When the degree of looseness increases, the vibration amplitude decreases, which may be due to the initial misalignment of the rotor system. The misalignment weakens gradually when the looseness occurs and deteriorates. This also shows that because the rotor system sometimes couples a variety of faults, the diagnosis of looseness faults simply by using changes in vibration amplitude is unreliable. Therefore, it is necessary to extract sensitive fault features to characterize looseness from the vibration signal better. From the spectrum in Figure 13, the amplitude of rotational frequency decreases, and the amplitude of double frequency increases. As the degree of looseness increases, the lower right part of the orbit becomes more and more "sharp" in Figure 14. It shows that the bearing seat collides with the foundation, consistent with the simulation results.
sometimes couples a variety of faults, the diagnosis of looseness faults simply by using changes in vibration amplitude is unreliable. Therefore, it is necessary to extract sensitive fault features to characterize looseness from the vibration signal better. From the spectrum in Figure 13, the amplitude of rotational frequency decreases, and the amplitude of double frequency increases. As the degree of looseness increases, the lower right part of the orbit becomes more and more "sharp" in Figure 14. It shows that the bearing seat collides with the foundation, consistent with the simulation results.  To effectively identify the data-driven model of the rotor-bearing test rig system and improve the accuracy of identification, the horizontal vibration signal is first comb-filtered by the Time Synchronous Averaging method to filter out noise or other frequency components that are not related to the rotational frequency. Then the vibration signal is downsampled with the downsampling frequency of 1280 Hz. Finally, the NARX models of the system in different looseness states are identified for further feature extraction. The terms and their corresponding coefficients of the NARX model under three different degrees of looseness are shown in Table 3.  To effectively identify the data-driven model of the rotor-bearing test rig system and improve the accuracy of identification, the horizontal vibration signal is first combfiltered by the Time Synchronous Averaging method to filter out noise or other frequency components that are not related to the rotational frequency. Then the vibration signal is downsampled with the downsampling frequency of 1280 Hz. Finally, the NARX models of the system in different looseness states are identified for further feature extraction. The terms and their corresponding coefficients of the NARX model under three different degrees of looseness are shown in Table 3. Table 3. NARX model terms and coefficients under three different looseness stiffness in the experiment.

No. Normal Weak Serious
The predicted output of the identified NARX model is compared with the experimental data processed by time synchronous averaging, as shown in Figure 15. The predicted output is highly consistent with the experimental results, which verifies the effectiveness of the NARX model under three different degrees of looseness. The identified NARX model is used to obtain the first three-order NOFRFs using the GALEs method. The second-order optimal weighted contribution rate of NOFRFs indicator SRm proposed in this paper is used to extract sensitive features and is compared with the conventional NOFRFs-based indicator Fe2. The two indicators under different looseness states are shown in Figure 16. As the degree of looseness increases, Fe2 first decreases and then increases. It is not monotonic and cannot be used to detect the change in the looseness state. The new indicator SRm showed a monotonous upward trend, increasing from 0.48 when no looseness occurs to 0.90 when severe looseness occurs, an increase of 89.7%. Therefore, the data-driven sensitive feature indicator SRm proposed in this paper can better reflect the occurrence and deterioration of system nonlinear damage. then increases. It is not monotonic and cannot be used to detect the change in the looseness state. The new indicator SRm showed a monotonous upward trend, increasing from 0.48 when no looseness occurs to 0.90 when severe looseness occurs, an increase of 89.7%. Therefore, the data-driven sensitive feature indicator SRm proposed in this paper can better reflect the occurrence and deterioration of system nonlinear damage.

Conclusions
This paper proposes a data-driven sensitive feature extraction method. First, the FROLS (Forward Regression with Orthogonal Least Squares) method is used to identify the inspected system's NARX (Nonlinear autoregressive with exogenous inputs) model. Then the system's NOFRFs (Nonlinear Output Frequency Response Functions) are obtained according to the GALEs (Generalized Associated Linear Equations) method. Based on the first three-order NOFRFs, the concept of the weighted contribution rate of NOFRFs is proposed, and the analytical expression of the second-order optimal weighted contribution rate of NOFRFs is derived and used as a sensitive feature. Further, the proposed method is applied to the feature extraction of the pedestal looseness of the rotor- then increases. It is not monotonic and cannot be used to detect the change in the looseness state. The new indicator SRm showed a monotonous upward trend, increasing from 0.48 when no looseness occurs to 0.90 when severe looseness occurs, an increase of 89.7%. Therefore, the data-driven sensitive feature indicator SRm proposed in this paper can better reflect the occurrence and deterioration of system nonlinear damage.

Conclusions
This paper proposes a data-driven sensitive feature extraction method. First, the FROLS (Forward Regression with Orthogonal Least Squares) method is used to identify the inspected system's NARX (Nonlinear autoregressive with exogenous inputs) model. Then the system's NOFRFs (Nonlinear Output Frequency Response Functions) are obtained according to the GALEs (Generalized Associated Linear Equations) method. Based on the first three-order NOFRFs, the concept of the weighted contribution rate of NOFRFs is proposed, and the analytical expression of the second-order optimal weighted contribution rate of NOFRFs is derived and used as a sensitive feature. Further, the proposed method is applied to the feature extraction of the pedestal looseness of the rotor-

Conclusions
This paper proposes a data-driven sensitive feature extraction method. First, the FROLS (Forward Regression with Orthogonal Least Squares) method is used to identify the inspected system's NARX (Nonlinear autoregressive with exogenous inputs) model. Then the system's NOFRFs (Nonlinear Output Frequency Response Functions) are obtained according to the GALEs (Generalized Associated Linear Equations) method. Based on the first three-order NOFRFs, the concept of the weighted contribution rate of NOFRFs is proposed, and the analytical expression of the second-order optimal weighted contribution rate of NOFRFs is derived and used as a sensitive feature. Further, the proposed method is applied to the feature extraction of the pedestal looseness of the rotor-bearing system. The dynamic model of the bearing-rotor system is established, and the piecewise linear function is used to characterize the pedestal looseness. The Newmark-β is used to solve the governing equation. Taking the unbalanced excitation of the rotor system as the input and the node response as the output, the NARX models under different degrees of looseness are identified. The model's predicted output is compared with the numerical response, and the weighted contribution rate is used to extract the sensitive feature of looseness. Simulation results show that the sensitive feature indicator proposed in this paper is more effective than the conventional NOFRFs-based indicator. Finally, a singledisc two-fulcrum rotor-bearing test rig was established, and different degrees of pedestal looseness was simulated. The simulation analysis and the effectiveness of the proposed data-driven feature extraction method are verified. In future, the application of data-driven diagnostic methods in quantitatively analyzing faults and detecting the location of damage to engineering structures will be further explored.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.