Kernel Regression Residual Decomposition-Based Polynomial Frequency Modulation Integral Algorithm to Identify Physical Parameters of Time-Varying Systems under Random Excitation

: The physical parameters (stiffness, damping) of time-varying (TV) systems under random excitation provide valuable information for their working condition but they are often overwhelmed by noise interference. To overcome this problem, this paper presents a novel multi-level kernel regression residual decomposition method, which can not only effectively separate each modal component from the raw vibration acceleration signal, but also eliminate noise interference. Additionally, the multiple degree-of-freedom (DOF) parameter identiﬁcation problem is transformed into a single DOF parameter identiﬁcation problem. Combined with the derived polynomial frequency modulation integral algorithm and the cross-correlation theory based on the fractional Fourier ambiguity function, a physical parameter identiﬁcation method is proposed. The method provides a new idea in modeling TV systems and identifying physical parameters under random excitation. To demonstrate the effectiveness of the proposed method, numerical simulations are conducted with three different cases of variation (variation, quadratic variation, and periodic variation) in time. Moreover, its robustness is evaluated by adding different signal-to-noise ratio levels of noise (20 dB, 50 dB, 100 dB) to the input vibration acceleration signal. The analysis results conﬁrm the performance of the proposed method for the parameter identiﬁcation of TV systems under random excitation.


Introduction
A large number of engineering structures in aerospace, civil engineering, and other fields are subjected to environmental random excitation and exhibit time-varying (TV) characteristics. However, the corresponding random excitation force is difficult to measure, which causes significant errors in parameter identification using traditional methods. Furthermore, even if the environmental random excitation force is stationary, the collected vibration signal of the TV structure is non-stationary, thus limiting the performance of parameter identification methods based on the short-term time-invariant assumption. In addition, the collected TV vibration signal lacks an exact mathematical analytical expression. These problems must be urgently addressed as they are crucial for structural design, vibration testing, and health monitoring. Therefore, there is a need to effectively carry out physical parameter identification for TV systems under random excitation [1][2][3].
Vibration analysis-based methods have been extensively studied and applied in the field of parameter identification [4,5], and have been proved as the most effective way. These methods can be classified into two types. The first type is based on the short-term time-invariant assumption [6][7][8], in which TV systems are analyzed. These methods utilize various signal processing tools to identify parameters within a short-term interval, achieving the purpose of identification through overall fitting. The second type is based on the accuracy of TV system parameter identification [30]. Polynomial chirplet transform (PCT) [31] is a recently developed TF method, which is more precise in TF characterization by introducing a frequency demodulated operator. The frequency demodulated operator can make the horizontal line fitting the TF characteristic rotate an angle adaptively, i.e., the parameter variation satisfies the short-term linear variation assumption. Therefore, the TF resolution is significantly improved, and more TV information [32] can be reflected. Compared with other TF methods based on the short-term time-invariant assumption, PCT based on the short-term linear variation assumption is more suitable for analyzing non-stationary signals. Peng et al. [33] investigated PCT and demonstrated that it can produce a good TF distribution for a wide range of signals with continuous instantaneous frequency (IF). Li et al. [34] optimized PCT and constructed a synthetic simulation signal to verify its effectiveness. Meanwhile, the optimized method was applied to field seismic data, and the tight sandstone gas reservoirs can be effectively identified. Deng et al. [35] compared PCT with two parameter identification methods based on the Hilbert transform and verified the parameter identification performance of PCT by experiments. However, PCT is not suitable for multi-modal compound signals. Therefore, the combination of PCT with multi-level KRRD can better track the TF variation.
The ambiguity function (AF) is widely used in the field of radar and sonar, which can reflect the resolution characteristics of two targets with different distances and speeds, i.e., the signal correlation in the time domain. Onur Celik et al. [36] proposed a generalized formulation based on AF, solving the problem of multiple-input multiple-output radar beampattern design successfully. Aiming at the uncertainty associated with random vibration signals, the fractional Fourier ambiguity function (FFAF) [37] was introduced, and the correlation function [38] was derived. Based on the cross-correlation theory, the relationship between the correlation function and TV physical parameters can be determined.
Based on the preceding analysis, the combination of the mentioned methods (multilevel KRRD, PCT, FFAF) shows great potential for application in identifying TV parameters under random excitation, and the remainder of the paper is organized as follows. The analysis process and specific theoretical derivations of multi-level KRRD, PCT and FFAF are described in Section 2. In Section 3, the random vibration model with TV parameters is presented. In Section 4, numerical simulations of a three DOF TV structure are given. In Section 5, the conclusions are presented to close the paper.

Overview of the Study
Parameter identification of TV systems is an important part of dynamic analysis. How to globally decompose the multi-component signal is the premise for addressing the parameter identification problem of a multiple DOF TV system under random excitation. The decomposed component signals contain valuable system information that can be utilized for further investigation. Furthermore, effectively modeling the decomposed component signals in the time domain, based on the short-term linear variation assumption, holds the key to enhancing the parameter identification accuracy. Consequently, this paper proposes a novel physical parameters identification method of TV systems under random excitation. To visually illustrate the proposed method, a detailed flowchart is presented in Figure 1, and the major study contents are as follows.
By employing multi-level KRRD, the collected raw acceleration signal can be decomposed into a set of acceleration components representing modal components. In this step, the multiple DOF parameter identification problem can thus be transformed into a single DOF parameter identification problem. At the same time, the SNR of the signal is improved by eliminating noise interference. Here, the single decomposed acceleration component is considered as the analysis object, which can be represented by a linear combination constructed by PCT basis functions. Firstly, the corresponding first and second derivative basis functions are derived using the basis function of PCT. Secondly, PCT is performed on the decomposed acceleration component based on the obtained derivative basis functions, and the corresponding velocity and displacement signals can be reconstructed, which exhibits significant potential in resolving TV differential equations. Then, the cross-correlation theory is employed, and a linear model is established between the TV physical parameters and the acceleration signal, thereby transforming the problem of solving the differential equation of the TV system into solving a linear equation. Finally, the TV physical parameters are identified by solving the simultaneous linear equations. this step, the multiple DOF parameter identification problem can thus be transformed into a single DOF parameter identification problem. At the same time, the SNR of the signal is improved by eliminating noise interference. Here, the single decomposed acceleration component is considered as the analysis object, which can be represented by a linear combination constructed by PCT basis functions. Firstly, the corresponding first and second derivative basis functions are derived using the basis function of PCT. Secondly, PCT is performed on the decomposed acceleration component based on the obtained derivative basis functions, and the corresponding velocity and displacement signals can be reconstructed, which exhibits significant potential in resolving TV differential equations. Then, the cross-correlation theory is employed, and a linear model is established between the TV physical parameters and the acceleration signal, thereby transforming the problem of solving the differential equation of the TV system into solving a linear equation. Finally, the TV physical parameters are identified by solving the simultaneous linear equations.

Multi-Level KRRD
The raw vibration acceleration signal s(t) can be expressed by multiple residual signals as follows:

Multi-Level KRRD
The raw vibration acceleration signal s(t) can be expressed by multiple residual signals as follows: in which n represents the number of the residual signals, r i (t) represents the i-th residual signal, and s ∑ 1 (t) represents the set of the decomposed residual signals when the decomposition level is [1, n].
When the decomposition level is [2, n], s ∑ 2 (t) is described as follows: Therefore, the 1-th decomposed residual signal can be obtained by subtraction, and the corresponding expression is shown as follows: Appl. Sci. 2023, 13, 8151 5 of 16 in which K 1 (·) represents the non-parametric NW kernel estimation with the Gaussian kernel function. λ represents the bandwidth which can control the estimated smoothness, and it is taken as the reciprocal of the sampling frequency. s(t) represents the Gaussian kernel function center, the same as the mean value of s(t).
Similarly, the 2-th residual signal can be given as follows: Based on the above analysis, the (n − 1)-th residual signal can be obtained by the following equation. Meanwhile, an interference signal r n (t) containing noise can be eliminated. To provide a more intuitive representation of the multi-level KRRD, the decomposition frame is shown in Figure 2.
Therefore, the 1-th decomposed residual signal can be obtained by subtraction, and the corresponding expression is shown as follows: in which ( ) 1 K ⋅ represents the non-parametric NW kernel estimation with the Gaussian kernel function. λ represents the bandwidth which can control the estimated smoothness, and it is taken as the reciprocal of the sampling frequency. ( ) s t represents the Gaussian kernel function center, the same as the mean value of s(t).
Similarly, the 2-th residual signal can be given as follows: Based on the above analysis, the (n − 1)-th residual signal can be obtained by the following equation. Meanwhile, an interference signal rn(t) containing noise can be eliminated. To provide a more intuitive representation of the multi-level KRRD, the decomposition frame is shown in Figure 2.

Polynomial Frequency Modulation Integral Algorithm
PCT is employed to better generate IF representation, and it can be defined as follows: in which t0 represents the initial moment, s(t) represents the analytic signal with the TV parameter, w(t − t0) represents the window function, = √−1, ω represents the IF, represents the adaptive frequency demodulated coefficient, and

Polynomial Frequency Modulation Integral Algorithm
PCT is employed to better generate IF representation, and it can be defined as follows: in which t 0 represents the initial moment, s(t) represents the analytic signal with the TV parameter, w(t − t 0 ) represents the window function, j = √ −1, ω represents the IF, c represents the adaptive frequency demodulated coefficient, and represents the TV demodulated operator.
For the analytic signal s(t), it can be defined by amplitude a(t) and frequency ϕ(t): In a short period of time, ϕ(t) shows a linear variation regularity, and the corresponding Taylor expansion is given as follows: in which ϕ (·) (t) represents the derivative of ϕ(t).
According to the Equations (6)- (8), the PCT amplitude of IF ϕ(t)(ω = ϕ(t)) can be obtained as follows: As shown in Equation (9), when , the amplitude takes the maximum value. Meanwhile, the TF energy concentration is the highest, i.e., the optimal TF representation is obtained.
In order to obtain the optimal demodulated operator, a novel parameter α is introduced here [39]. The fitting angle in the TF domain is defined as arctan(c), and the parameter α can be defined by sampling time t s and frequency f s as follows: The discrete form of parameter α and the update form of Equation (6) are shown as follows: The corresponding TF spectrum is given as follows: in which c m is the discrete value corresponding to the maximum TF amplitude [c i , ω i ].
[c i , ω i ] is given as follows: Based on Equation (14), the basis function h i t i ,ω i ,c i (t) of the i-th component signal can be described as follows: Then, the first and second derivatives of the basis function h i t i ,ω i ,c i (t) with respect to time are calculated as follows: .. In essence, PCT is still WT and supports signal reconstruction. In addition, the first derivative and the second derivative basis functions satisfy the admissible condition and decay condition [35], which are given as follows: whereĥ i (ω) is the Fourier transform of the basis function.
For the acceleration signal ..
x(t), the corresponding velocity .
x(t) and displacement signal x(t) are calculated through the once and twice integrations, which are defined as follows: .. ..
x i (t) (20) in which · and · represent the once and twice integrations, respectively.
PCT based on ..
x i (t), and the transform is carried out by the partial integration method. The corresponding computational process is shown as follows: Combined with Equation (19), the first term in the above equation is updated as follows: ..
Then, Equation (21) is updated as follows: Similarly, PCT based on ..
x i (t) and the corresponding equation is given as follows: The first two terms in Equation (24) can be written as follows: ..
Equation (24) is updated as follows: In order to reconstruct the velocity and displacement signals, Equations (23) and (26) can be further rewritten as follows: PCT..
Then, the corresponding signals are rewritten as follows: .. ..

Fractional Fourier Ambiguity Function
As a linear operator, fractional Fourier transform (FRFT) can realize the counterclockwise rotation at any angle of the signal around the TF axis, and the p-order FRFT is given by: in which K 2 (u,t) represents a kernel function, and it is given as follows: in which β = pπ/2. From Equation (32), it can be seen that FRFT degenerates into the traditional FT when p = 1. According to PCT, a polynomial frequency modulation signal is defined as s 1 (t) = A 1 · exp j ωt + at 2 , in which A 1 and a represent the amplitude and phase coefficient, and the corresponding AF is given as follows: (33) in which τ represents the time delay. The autocorrelation function of s 1 is shown as follows: Therefore, the FFAF is defined as follows: Based on the above analysis, the correlation function takes the maximum value when τ = 0.
in which E{·} represents the expectation operator. Meanwhile, a frequency modulation signal s 2 (t) = A 2 · exp j(ω + ω)t + at 2 , in which ω represents the variation of frequency, is constructed. Based on the orthogonality of s 1 (t) and s 2 (t), the corresponding cross-correlation function can be obtained as follows:

Physical Parameter Identification under Random Excitation
For an R DOF TV system, the vibration differential equation can be defined as follows: ..
x(t) and x(t) are the system acceleration, velocity, and displacement vectors; and f (t) are the excitation vectors exerted on the system.
On the basis of modal superposition theory, the system acceleration vector can be regarded as the sum of multiple modal components. The acceleration vector of r-th DOF can be expressed as ..
x ri (t)(r = 1, 2, . . . , R). Under the condition that only the acceleration signal is known, the PCT integral algorithm is employed, and Equation (38) is rewritten as follows: ..
x(t) = f (t) (39) Further, Equation (39) can be updated with amplitude A pq and phase as follows: The amplitude A pq is given as follows: Expectation operations are performed on both sides, and Equation (40) can be updated as follows: Here, random Gaussian white noise is set as the system excitation force, and its correlation function is defined as follows: Therefore, the right side of Equation (42) can be rewritten as follows: Based on Equations (37) and (43), the correlation function matrix can be obtained as follows: Equation (42) is rewritten as follows: Finally, the TV dynamic linear equation is updated, and the physical parameters are obtained by solving the simultaneous equations.

Numerical Simulation
In this section, a simulated three DOF model incorporating TV stiffness and damping is constructed to testify the effectiveness of the proposed method, which is illustrated in Figure 3. The initial structural parameters of the TV model are given as follows, and the excitation force was random Gaussian white noise. The sampling frequency was set to 1200 Hz, and the sampling time was set to 10 s.

Numerical Simulation
In this section, a simulated three DOF model incorporating TV stiffness and damping is constructed to testify the effectiveness of the proposed method, which is illustrated in Figure 3. The initial structural parameters of the TV model are given as follows, and the excitation force was random Gaussian white noise. The sampling frequency was set to 1200 Hz, and the sampling time was set to 10 s. This paper considers the identification of stiffness and damping, as these parameters are more likely to change compared to the structural mass in practical engineering structures. Moreover, assessing stiffness and damping is more meaningful in terms of structural damage and vibration control. In order to simulate the dynamic characteristics of practical structures more accurately, three different cases of variation in time of stiffness and damping parameters are investigated, namely, linear variation, quadratic variation, and periodic variation. The detailed cases are given as follows: ( ) The raw vibration acceleration signal under a random excitation force was computed through the Newmark algorithm, and the corresponding time-domain waveform is shown in Figure 4a. In order to extract the modal component corresponding to each order DOF, the multi-level KRRD was then applied to the raw acceleration signal. The decomposed components are shown in Figure 4b-d, respectively. Subsequently, the polynomial frequency modulation integral algorithm was employed to each modal component. Figure 5 illustrates the signal reconstruction result of velocity and displacement signals (blue) obtained from the integral algorithm. Additionally, the velocity and displacement signals This paper considers the identification of stiffness and damping, as these parameters are more likely to change compared to the structural mass in practical engineering structures. Moreover, assessing stiffness and damping is more meaningful in terms of structural damage and vibration control. In order to simulate the dynamic characteristics of practical structures more accurately, three different cases of variation in time of stiffness and damping parameters are investigated, namely, linear variation, quadratic variation, and periodic variation. The detailed cases are given as follows: The raw vibration acceleration signal under a random excitation force was computed through the Newmark algorithm, and the corresponding time-domain waveform is shown in Figure 4a. In order to extract the modal component corresponding to each order DOF, the multi-level KRRD was then applied to the raw acceleration signal. The decomposed components are shown in Figure 4b-d, respectively. Subsequently, the polynomial frequency modulation integral algorithm was employed to each modal component. Figure 5 illustrates the signal reconstruction result of velocity and displacement signals (blue) obtained from the integral algorithm. Additionally, the velocity and displacement signals (red) obtained from numerical simulation are also included in Figure 5 to further verify the reconstruction performance of the integral algorithm. The time-domain waveforms in the first line of Figure 5 represent the displacement signals, while the second line represents the velocity signals, with each column corresponding to a different DOF. By comparing the blue and red waveforms in Figure 5, it can be seen that the proposed integral algorithm exhibits good signal reconstruction performance.       Finally, the physical parameters (stiffness, damping) can be identified according to Equation (47), and the corresponding identification results of linear variation, quadratic variation, and periodic variation are given in Figures 6 and 7, respectively. As shown in Figures 6 and 7, the identification values fluctuate around the theoretical values for the three different cases. Consequently, the proposed method can effectively capture the variation of physical parameters with time and track the trend of various variations. However, it is worth noting that the identification of damping is more volatile than that of stiffness, as seen in Figure 7. This is because the damping coefficient, due to its nature and usually small values in structural contexts, is more prone to error than the stiffness. For the oscillation mode of the damping, one of the reasons is the oscillation of the frequency and amplitude of the single component signal when it is constructed by the integral algorithm under random excitation. Meanwhile, errors in different sampling points during the solution of the simultaneous equations can also affect the identification results. In addition to this, the interference of the random excitation force induces non-linear variation in the TV system, resulting in errors in the damping identification results.

Conclusions
Aiming at the difficulty of physical parameter identification under random excitation, a novel multi-level KRRD-based polynomial frequency modulation integral analysis method was proposed. The method utilized multi-level KRRD to decompose the collected vibration acceleration signal into modal component signals corresponding to each DOF. These modal component signals contain crucial TV system information, such as stiffness

Conclusions
Aiming at the difficulty of physical parameter identification under random excitation, a novel multi-level KRRD-based polynomial frequency modulation integral analysis method was proposed. The method utilized multi-level KRRD to decompose the collected vibration acceleration signal into modal component signals corresponding to each DOF. These modal component signals contain crucial TV system information, such as stiffness To further indicate the robustness of the proposed method, different SNR levels of noise (20 dB, 50 dB, 100 dB) were added to the raw vibration acceleration signal. Meanwhile, a statistical indicator called mean absolute percentage error (MAPE) was introduced here, which can evaluate the accuracy of the identified results. The same procedure shown in Section 3 was performed, and the detailed identification error is shown in Table 1.
in which N S represents the total sampling point number, ι i represents the theoretical value, and ι i represents the identification value. From Table 1, it can be seen that under different SNR levels, the error of stiffness can be maintained within 3%, and the damping can still ensure good identification accuracy. Through the aforementioned analysis, it can be deduced that the method proposed in this paper is effective in overcoming random noise and has good robustness and TV physical parameter traceability.

Conclusions
Aiming at the difficulty of physical parameter identification under random excitation, a novel multi-level KRRD-based polynomial frequency modulation integral analysis method was proposed. The method utilized multi-level KRRD to decompose the collected vibration acceleration signal into modal component signals corresponding to each DOF. These modal component signals contain crucial TV system information, such as stiffness and damping. In this step, the parameter identification problem of a multi-DOF system was transformed into the parameter identification problem of a single DOF system, which facilitated subsequent analysis. Moreover, the noise interference component was eliminated, and the anti-noise performance of the proposed method was improved. Additionally, the method employed the polynomial frequency modulation integral algorithm to reconstruct the velocity and displacement signals based on the single decomposed modal component signal. By combining cross-correlation theory, the dynamic differential equation of the TV system was modeled in the time domain, and the physical parameters were identified by solving the simultaneous linear equation.
Based on numerical simulation, the identification effect of physical parameters under three different cases of variation (linear variation, quadratic variation, periodic variation) was conducted on an established three DOF TV model. The numerical results demonstrated the capability of the proposed method to effectively track various variations of physical parameters over time. Furthermore, the robustness of the proposed method was verified by adding different SNR levels of noise (20 dB, 50 dB, 100 dB) to the raw vibration acceleration signal. The identification errors of stiffness and damping are guaranteed to remain within a stable range, indicating the high identification accuracy of the proposed method even under noise interference.
Only the vibration acceleration signal was used for parameter identification, which facilitates signal acquisition in practical engineering and improves the application po-tential of the proposed method. However, the proposed method is essentially based on the assumption of the short-term linear variation, and it has limitations in identification performance for the TV system with short-term interval non-linear variation characteristics. To enhance the adaptability and parameter identification accuracy of the proposed method, the authors will further study the non-linear fitting of the TV signal.