Abstract

Operational modal parameter identification is a tough problem in aerospace engineering due to the complex mechanics environment, various noises, and limited computational resources. In this paper, a novel, recursive, robust, and high-efficiency modal parameter identification approach is proposed for this issue. The kernelized time-dependent autoregressive moving average (TARMA) model is adopted to model the nonstationary responses, a recursive estimator is established based on the maximum correntropy criterion, and sliding-window technique is applied to fix the computational complexity, which ensures the approach its estimation accuracy, robustness, and high efficiency. Finally, steps of the identification procedure and model selection are presented. An experimental scheme is proposed for validation, and the proposed approach is comparatively assessed against the classical recursive pseudo-linear regression TARMA method via Monte Carole tests of a time-varying experimental system. The results of the comparative study demonstrate that the proposed method achieves similar estimation accuracy and higher computation efficiency under the Gaussian environment. Moreover, a superior estimation accuracy and enhanced robustness are rendered under additive non-Gaussian impulsive noise.

1. Introduction

Aerospace structures are typical time-varying structures, such as airplanes with varying additional aerodynamic effects in flight [1, 2], launch vehicle with varying fuel mass [35], and satellite with deployable antenna [6, 7]. The great increase in the complexity and size of aerospace structures has brought remarkable difficulties and challenges to the system level dynamics simulation and test on the ground. Therefore, online time-varying operational modal parameter identification is imperative for better flight performance, longer service lifetime, and more reliable control ability in aerospace engineering.

The problem of operational modal parameter identification of time-varying structures involving exclusively available response signals is also referred as the output-only time-varying identification problem [8]. Time-dependent autoregressive moving average (TARMA) model-based identification methods have gained much more attention due to their attractive features [9, 10].

Over the past decades, many efforts have been undertaken to develop the batch identification methods based on the TARMA model [1113]. Most batch identification methods have good tracking and estimation accuracy. However, great computational cost and postprocessing manner limit their application to online identification. Also, a lot of work has been done in recursive methods based on the TARMA model. Xie and Evans [14] first used exponentially weighted Legendre functions to model time-varying parameters of the TARMA model. Poulimenos and Fassois [15] further extended the recursive methods by selecting basis functions of orthogonal functions and proposed the recursive pseudo-linear regression TARMA (RPLR-TARMA) method. However, as running basis methods mentioned above, both of them encountered numerical problems when basis sequences are not bounded. In order to improve the numerical problems, Yang [16] proposed a kernelized TAR model, which represented the time-varying parameters by using a series combination of kernel functions. Ma et al. [17] extended the kernelized TAR model to the TARMA model and proposed a recursive estimator. However, the kernel-based TARMA method does not perform well when signal-to-noise ratio is low, or there is a significant non-Gaussian impulsive noise, which limit its application to aerospace structures with rugged mechanical environment. Yu et al. [18] proposed a maximum correntropy criterion based on the TARMA model; it shows charming ability to deal with significant non-Gaussian impulsive noise and provide a new manner to extend the recursive TARMA method to the complex mechanics environment.

In order to extend the recursive TARMA method for online operational modal parameter identification of aerospace structures, based on the conditions of the complex mechanical environment, various noises, and limited computational resources, a novel, robust, recursive estimator is established based on the kernelized TARMA model with the maximum correntropy criterion. The kernelized TARMA model represents time-varying parameters as a linear combination of kernel functions; by using “kernel trick” and sliding-window technique, the computation complexity is fixed and only related to the window length. Moreover, by considering the ridge regression criterion and the maximum correntropy criterion [1923], the proposed estimator is also capable of the ill-posed problem and the non-Gaussian noise problem. Finally, a recursive workflow of the sliding-window exponentially weighted recursive maximum correntropy ridge regression-kernelized TARMA (SWRMCRR-KTARMA) method is presented by taking the assumption that the time-varying system changes slowly in a certain time scale.

2. Method

2.1. Kernelized TARMA Model and Maximum Correntropy Criterion

The TARMA model with na and nc designating, respectively, its autoregressive (AR) and moving average (MA) orders, is of the general form [12]with designating the normalized discrete time, the discrete-time nonstationary vibration signal modeled, an unobservable uncorrelated innovation sequence with zero mean and time-varying nonsingular covariance matrix , and the model’s time-dependent AR and MA parameter matrices. stands for normally independently distributed random variables with indicated mean and covariance.

By using kernel transformation, the AR and MA parameter matrices can be expressed as

The kernelized TARMA model is given by [24]where is the Kronecker product, is the time-dependent parameter vector, is the regression matrix, and and are the corresponding kernelized counterparts. is the transformed feature vector lying in the feature space; by using a Gaussian kernel, the “kernel trick” can be expressed as follows:where designates the adjustable parameter.

The system’s “frozen-time” characteristic matrix equation is given by [25]with designating the backshift operator, . The roots of the system can be derived from a general eigenvalue problem [26]where and stand for the rth eigenvalue and eigenvector of . The matrix is constructed from the AR parameter matrices as

The “frozen-time” modal frequencies and damping ratios of the system can be computed by

In general, correntropy is a local similarity measure between two arbitrary random variables of , with joint distribution function , and is defined by [27]where is a kernel function that satisfies Mercer theory and denotes the expectation operator. In this note, the shift-invariant Mercer kernel is the Gaussian kernel given bywhere , and stands for the kernel bandwidth.

For a ridge regression problem, with parameter weight vector , the corresponding optimal problem based on the maximum correntropy criterion can be expressed as follows:

2.2. Recursive Maximum Correntropy Ridge Regression-Kernelized TARMA Method

The exponentially weighted mechanism is used to put more emphasis on the recent data. In order to fix the computational complexity, a sliding-window technique is adopted, allowing the proposed method can be operated in a recursive manner with low computational cost. Given the forgetting factor and the regularization parameter , the weighted regularized maximum correntropy cost function is defined asand the solution to this optimal problem is given bywithwhere .

Denote and ; the following equation can be received:where can be obtained by removing the first row and column of .

By using “kernel trick,” can be computed as follows:where .

By using the upsized matrix inversion formula, the update process of iswith and .

Define

By using the downsized matrix inversion formula, the following equation can be obtained:

Considering equations (19) and (13), the kernelized parameter matrix can be obtained:

Finally, time-dependent AR and MA parameter matrices can be computed by

The workflow of the SWRMCRR-KTARMA method is summarized in (Algorithm 1).

Initialization:
Computation :
Computation :
 Downsizing: obtain by removing the first row and column of .
 Compute according to the downsized matrix inversion formula.
TARMA model parameter estimation:
2.3. Model Structure Selection

The model structure parameters of the proposed SWRMCRR-KTARMA method include the kernel bandwidth and , the forgetting factor , the sliding-window length , the regularization parameter , and AR and MA orders na and nc.

Modal structure selection is generally based upon trail-and-error or integer optimization schemes, and the “fitness” one is finally selected. The “fitness” function can be established based on many rules; in this note, a normalized residual sum of squares (RSS) criterion is chosen as the “fitness” function, which can be seen as

Model structure selection of SWRMCRR-KTARMA can be tackled via the following three basic steps:Step 1: given initial values of , , , , and based on prior knowledge of model parameter variation.Step 2: select AR and MA orders based on the minimization of the normalized RSS:(a)Given the integer search space with na = nc, and the fitness AR order is selected when the decrease of RSS values is quite subtle(b)Given the integer search space , and the best fit MA order can be selected in the same way as that of the AR orderStep 3: optimizing initial parameters based on the model orders obtained in Step 2.

This process can be done while the fitness or suitable one is found. A detailed optimization scheme is beyond the scope of the current note, and the initial parameters will be selected based on prior knowledge or the normalized RSS criterion.

Remark 1. The kernel bandwidth reflects the nonlinear power of the proposed method; it should be selected based on different problems. The kernel bandwidth controls the robustness to impulsive noise of the proposed method; it should be selected properly to avoid the “over-robust” problem.

Remark 2. The selection of the forgetting factor and sliding-window length is a coupled problem. Both of them affect the tracking ability of the method, and the forgetting factor should be very close to one in case N is relatively small. In this note, the forgetting factor is initially recommended, and the selection of is based on the normalized RSS criterion.

Remark 3. The regularization parameter deals with the overfitting problem, which controls the variance of estimated regression parameters; however, it does this at the expense of adding bias to the estimate. The selection of is a compromise issue.

3. Experimental Verification

In order to effectively and reliably assess the SWRMCRR-KTARMA method, a time-varying barrel experimental structure is adopted to simulate the typical mass-varying feature of the launch vehicle by draining off contained water for validation, along with the classical exponentially weighted recursive pseudo-linear regression TARMA (EWRPLR-TARMA) method [15] for comparison. The experimental scheme is as follows.

The “frozen” modal experiment is first taken out, and the corresponding “frozen” modal parameters are referred as baseline value. Then, the time-varying experiment is carried out, and the nonstationary response signal is recorded by transducers. Various intensities of numerical non-Gaussian impulsive noise are added into the response signal directly for evaluating the robustness under non-Gaussian noise, as it is limited by the present experimental condition to simulate the non-Gaussian impulsive phenomena physically. Finally, the SWRMCRR-KTARMA method and the EWRPLR-TARMA method are validated and compared based on the experimental data.

3.1. Mass-Varying Test Structure

A free-free barrel mass-varying laboratory system is set up to simulate the mass-varying feature of the carrier rocket by draining off water. The test structure is hung from the aluminum alloy support frame by the nylon rope to simulate the free-free boundary conditions. The steel barrel is 2.5 m high, and the outer diameter is 0.05 m, and the thickness is 0.001 m. Total weight of the barrel is 2.4 kg, and two 15 kg added masses are attached to it in order to reduce the natural frequency of the system, as shown in Figure 1. The tankage of the barrel allows mass of water varying from 15.7 kg (full state) to 0 kg (empty state), and the diameter of the hole at the bottom controls the flow rate. The exciter system consists of an exciter (Modal Shop 2025E) and a power amplifier (SmartAmp 2100E21-400). An impedance head (PCB 288D01), a laser range finder (Y1TA100MHT88), and 9 accelerometers (PCB 333B30) are, respectively, used as the force, distance, and motion transducers. The acquisition module is a LMS SCADAS III system. Control systems consist of a steering engine (JX PDI6221MG) and a remote control (Futaba 14SG).

The experimental mass-varying system is characterized by the position of the water level, and the water level is time-dependent. Considering the frozen-time assumption, the time-varying experimental system can be treated as a series time-invariant system, and the time-invariant modal analysis methods can be applied. And the frozen modal experiments are conducted as follows.

During the experiment, the water level starts at 0.25 m (the bottom of the barrel is the zero point) and ends at 2 m, which means the measurements of the laser range finder start at 2.25 m and end at 0.5 m because the laser range finder is on the top of the barrel. We divided the level range, 0.25 m–2 m, into 35 equal segments of 0.05 m. The “frozen-time” experiment is carried out by adding about 395 ml water into the barrel in each run, and a random excitation is added at 0.5 m, and 9 channels of accelerometers are settled along the barrel from bottom to top, as shown in Figure 1. Figure 2 shows the first two modal frequency and damping ratio estimates. The horizontal axis is the measurement of the laser range finder; the vertical axis is the modal frequency in Figure 2(a) and damping ratio in Figure 2(b).

The time-varying experiment is operated by opening the hole at the bottom of the barrel, and a random excitation is added at 0.5 m at the same time. Total 9 channels are used to obtain the nonstationary response signals. The response signals are sampled at 1,024 Hz, the time duration is 16 s, and the total number of the samples is 16,384.

3.2. Non-Gaussian Noise Simulation

In order to model the outliers/impulsive noise that exist in real test cases, a class of alpha-stable distribution is discussed. There is no closed-form expression for the probability density function of alpha-stable distributions, but its characteristic function [28] is givenwhereand is the characteristic exponent satisfying , which controls the thickness of tails in the distribution. is the location parameter , which corresponds to the mean for and the median for . is the dispersion parameter (), which determines the spread of the density around its location parameter. is the symmetry parameter (); when , the distribution is symmetric around the location parameter and referred to as the symmetry alpha-stable distribution (). The overall shape and the tails of the probability density functions of alpha-stable distributions are presented in Figure 3. In order to compare the level of a desired signal to an impulsive noise, a mixture signal-to-noise ratio (MSNR) [29] is defined aswhere is the variance of the signal and is the dispersion parameter of the alpha-stable distribution. In this note, a class of symmetric alpha-stable distribution () is mainly considered.

3.3. Comparative Results and Discussion

Various capabilities of the previously mentioned methods are examined in this section. The main focus here is placed upon the assessment of the following characteristics: (a) achievable time-dependent natural frequency accuracy, (b) computational complexity, and (c) robustness to non-Gaussian impulsive noise.

The proposed SWRMCRR-KTARMA method is employed to identify the experimental time-varying structure. The initial values, , are selected via prior knowledge, and model order selection is achieved via the RSS criterion as shown in Figure 4, which suggests the fitness AR and MA order . The forgetting factor is finally selected as the sliding-window length is relatively small. Also, the suggested AR, MA orders and forgetting factor of the EWRPLR-TARMA method are , and , which can be achieved via the RSS criterion as shown in Figure 5. The finally selected characteristics of the SWRMCRR-KTARMA method and its LS counterpart EWRPLR-TARMA method are summarized in Table 1.

The modal frequency identification results of SWRMCRR-KTARMA and EWRPLR-TARMA can be seen in Figure 6. Evidently, the performance of the SWRMCRR-KTARMA method achieves similar estimation and tracking accuracy with the EWRPLR-TARMA method under none additive noise, and it shows superior accuracy and robustness under non-Gaussian additive noise.

In order to compare these methods by using natural frequency estimates in a quantitative manner, a mean absolute error (MAE) is introduced aswith designating the baseline value of the modal frequency, its estimated value based on the ith Monte Carlo run, and the effective number of estimations for the ith Monte Carlo run.

The MAE values and the average effective number obtained by each method are summarized in Figures 7 and 8, which illustrate the performance of the SWRMCRR-KTARMA method in estimation accuracy is a bit poor (quite similar) than the EWRPLR-TARMA method under none additive noise (relative MAE error is less than 0.5% for two modes), and the performance in robustness is not as good as the EWRPLR-TARMA method (the relative effective number is 3.9% and 10.1% for mode 1 and mode 2). Under additive non-Gaussian impulsive noise, with MSNR decreasing, the performance of the EWRPLR-TARMA method degenerates gradually, while the SWRMCRR-KTARMA method keeps an accurate and robust modal frequency estimation. The main reason to explain this comparative result is that the EWRPLR-TARMA method is based on the LS estimator, which is optimal under the Gaussian noise assumption. In another word, the maximum correntropy criterion-based SWRMCRR-KTARMA method loses its optimality to enhance its robustness under non-Gaussian noise.

The time and memory complexities of the SWRMCRR-KTARMA method are , which is mainly determined by the sliding-window length. However, the time and memory complexities of the EWRPLR-TARMA method are [17], which is related to the dimensionality of the regression matrix and covariance matrix. Figure 9 illustrates the normalized processing CPU times required for model parameter estimation, which indicate the computational complexity associated with each method. Evidently, the proposed method attains much better computational performance than the EWRPLR-TARMA method, especially when the number of signal channels (parameter ) is quite large which is common in the large aerospace structure test.

4. Conclusion

A novel kernel recursive method was proposed for operational modal parameter identification of aerospace structures. The SWRMCRR-KTARMA method was validated by an experimental time-varying system and compared with the conventional least-squares criterion-based EWRPLR-TARMA method. The results indicated that the proposed method was able to track the dynamics of the time-varying system and achieved similar accuracy, lower computational complexity, and enhanced online identification capability compared with the existing EWRPLR-TARMA method under Gaussian conditions. Moreover, it exhibited its advantages on superior estimation accuracy and robustness under additive non-Gaussian noise. The proposed method extends the application range of recursive TARMA methods in aerospace engineering to the conditions of the complex mechanics environment, various noises, and limited computational resources.

Nomenclature

:Time-dependent autoregressive parameter matrix
:Time-dependent moving average parameter matrix
:Characteristic matrix
:Uncorrelated innovations, m
:rth natural frequency
:Transformed feature vector
:rth eigenvector
:Time-dependent parameter vector
:Discrete-time nonstationary vibration signal, m
:Covariance matrix
:Regression matrix
:rth eigenvalue
:rth damping ratio
:Characteristic exponent
:Normalized discrete time, s
:Gaussian kernel bandwidth parameter
:Forgetting factor
:Regularization parameter
:Sliding-window length
na:Autoregressive order
nc:Moving average order
:Location parameter
:Dispersion parameter
:Symmetry parameter.

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that they have no conflicts of interest to this work.