# Recursive Method in Modal Parameter Identification of Aerospace Structures under Non-Gaussian Noise

Authors: Lei Yu; Yong-li Zhang; Meng-di Yuan; Rui-qing Liu; Qi Zhang
Journal: Shock and Vibration
 (2020)

---

## 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.


---

## Body

## 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 [3–5], 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 [11–13]. 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 [19–23], 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 withna and nc designating, respectively, its autoregressive (AR) and moving average (MA) orders, is of the general form [12](1)xt+∑i=1naAitxt−i=et+∑i=1ncCitet−i,et∼NID0,Σt,with t designating the normalized discrete time, xt the discrete-time nonstationary vibration signal modeled, et an unobservable uncorrelated innovation sequence with zero mean and time-varying nonsingular covariance matrix Σt, and Ait,Cit the model’s time-dependent AR and MA parameter matrices. NID⋅ 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(2)Ait≈ω1,1AitTut…ω1,kAitTut⋮⋱⋮ωk,1AitTut…ωk,kAitTut,Cit≈ω1,1CitTut…ω1,kCitTut⋮⋱⋮ωk,1CitTut…ωk,kCitTut.The kernelized TARMA model is given by [24](3)xt=−∑i=1naAitxt−i+∑i=1ncCitet−i+et,⇔xt=wtTφt+et,⇔xt=w1,1A1tT…w1,kA1tT…w1,1CnctT…w1,kCnctT⋮⋱⋮⋱⋮⋱⋮wk,1A1tT…wk,kA1tT…wk,1CnctT…wk,kCnctT−xt−1⊗ut⋮et−nc⊗ut+et,⇔xt=wtφt⊗ut+et,⇔xt=wtTφt+et,where ⊗ is the Kronecker product, wt is the time-dependent parameter vector, φt is the regression matrix, and ωt and ϕt are the corresponding kernelized counterparts. ut is the transformed feature vector lying in the feature space; by using a Gaussian kernel, the “kernel trick” can be expressed as follows:(4)utiTutj=κkti,tj=exp−ti−tj22σk2,where σk designates the adjustable parameter.The system’s “frozen-time” characteristic matrix equation is given by [25](5)AB,t=I+∑i=1naAitBi=0,with B designating the backshift operator, Bixt=xt−i. The roots of the system can be derived from a general eigenvalue problem [26](6)Dt−λrtIVrt=0,where λrt and Vrt=λr−naLr,…,λr−1LrT stand for the rth eigenvalue and eigenvector of Dt. The matrix Dt is constructed from the AR parameter matrices as(7)Dt=0I…0⋮⋮⋱⋮00…I−Anat−Ana−1t…−A1t.The “frozen-time” modal frequencies and damping ratios of the system can be computed by(8)frt=lnλrt2πTs,ζrt=−Relnλrtlnλrt.In general, correntropy is a local similarity measure between two arbitrary random variables ofX,Y∈ℝ, with joint distribution function FXYx,y, and is defined by [27](9)VX,Y=EκcX,Y=∫κcx,ydFXYx,y,where κ⋅ is a kernel function that satisfies Mercer theory and E⋅ denotes the expectation operator. In this note, the shift-invariant Mercer kernel is the Gaussian kernel given by(10)κcx,y=Ge=exp−e22σc2,where e=x−y, and σc>0 stands for the kernel bandwidth.For a ridge regression problem, with parameter weight vectorw, the corresponding optimal problem based on the maximum correntropy criterion can be expressed as follows:(11)maxwJMCC=maxw∑i=1NEexp−e2i2σc2−γw2.


### 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 factor0<λ≤1 and the regularization parameter γ≥0, the weighted regularized maximum correntropy cost function is defined as(12)maxωt∑τ=t−N+1tλNκcxτ−ωtTϕτ−12λNγωt2,and the solution to this optimal problem is given by(13)ωt=Ψ⟶tΨ⟶tTΨ⟶t+λNγσcΓ⟶2t−1−1x⟶t,with(14)Ψ⟶t≜ϕt−N+1,ϕt−N+2,…,ϕt,x⟶t≜xt−N+1,xt−N+2,…,xtT,Γ⟶t≜diagλN−1Et−N+1,λN−2Et−N+2,…,Et,where Et≜exp−1/2σc2etTet.DenoteQ⟶t=Ψ⟶tTΨ⟶t+λNγσc2Γ⟶t−1 and k⟶t=ϕt−N+1,…,ϕt−1Tϕt; the following equation can be received:(15)Q⟶t=Q⟶∗t−1k⟶tk⟶tTϕtTϕt+λNγσc2Et−1,where Q⟶∗t−1 can be obtained by removing the first row and column of Q⟶t−1.By using “kernel trick,”ϕtTϕt can be computed as follows:(16)ϕtTϕt=ϕt⊗utTϕt⊗ut=ϕtTϕtκkt,t,where κkti,tj=exp−ti−tj2/2σk.By using the upsized matrix inversion formula, the update process ofQ⟶t−1 is(17)Q⟶t−1=Q⟶∗t−1−1+z⟶tz⟶tTr⟶t−1−z⟶trt−1−z⟶tTr⟶t−1r⟶t−1,with z⟶t=Q⟶∗t−1−1k⟶t and r⟶t=λNγσc2Et−1+ϕtTϕt−z⟶tTk⟶t.Define(18)Q⟶t−1−1=A⟶B⟶C⟶D⟶.By using the downsized matrix inversion formula, the following equation can be obtained:(19)Q⟶∗t−1−1=D⟶-C⟶A⟶−1B⟶.Considering equations (19) and (13), the kernelized parameter matrix ωt can be obtained:(20)ωt=Ψ⟶tQ⟶t−1x⟶t.Finally, time-dependent AR and MA parameter matrices can be computed by(21)wtT=A1t,…,Anat,C1t,…,Cnct=wtTIkna+nc⊗ut=Ψ⟶tQ⟶t−1x⟶tTIkna+nc⊗ut=Q⟶t−1x⟶tTφt−N+1T⊗κkt−N+1,t⋮φtT⊗κkt,t.The workflow of the SWRMCRR-KTARMA method is summarized in (Algorithm1).Algorithm 1: SWRMCRR-KTARMA method.
Initialization:w^0T=0,σc,σk,λ,γ,NComputation1≤t≤N:e^tt−1=xt−x^tt−1=xt−w^t−1TφtEt=exp−1/2σc2e^tt−1Te^tt−1kt=ψ˜t−1Tϕt=ϕ1Tϕt⋮ϕt−1Tϕt=φ1Tφtκk1,t⋮φt−1Tφtκkt−1,tzt=Qt−1−1ktrt=λtγσc2Et−1+φtTφtκkt,t−ztTktQt−1=Qt−1−1+ztztTrt−1−ztrt−1−ztTrt−1rt−1w^tT=Qt−1x˜tTφ1T⊗κk1,t⋮φtT⊗κkt,tComputationt>N:Downsizing: obtainQ⟶∗t−1 by removing the first row and column of Q⟶t−1.ComputeQ⟶∗t−1−1 according to the downsized matrix inversion formula.e^tt−1=xt−w^t−1TφtEt=exp−1/2σc2e^tt−1Te^tt−1k⟶t=Ψ⟶t−1Tϕt=ϕt−N+1Tϕt⋮ϕt−1Tϕt=ϕt−N+Tϕtκkt−N+1,t⋮ϕt−1Tϕtκkt−1,tz⟶t=Q⟶∗t−1−1ktr⟶t=λNγσc2Et−1+φtTφtκkt,t−z⟶tTk⟶tQ⟶t−1=Q⟶∗t−1−1+z⟶tz⟶tTr⟶t−1−z⟶tr⟶t−1−z⟶tTr⟶t−1r⟶t−1w^t=Q⟶t−1x⟶tTφt−N+1T⊗κkt−N+1,t⋮φtT⊗κkt,tTARMA model parameter estimation:Α1t,...,Αnat,C1t,...,Cnct=w^tT


### 2.3. Model Structure Selection

The model structure parameters of the proposed SWRMCRR-KTARMA method include the kernel bandwidthσc and σk, the forgetting factor λ, the sliding-window length N, 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(22)RSS≜∑t=1Ne^tTe^tN.Model structure selection of SWRMCRR-KTARMA can be tackled via the following three basic steps:Step 1: given initial values ofσc, σk, λ, N, 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 withna = nc, and the fitness AR order n^a is selected when the decrease of RSS values is quite subtle(b)
Given the integer search space0,n^a, and the best fit MA order n^c 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σk reflects the nonlinear power of the proposed method; it should be selected based on different problems. The kernel bandwidth σc 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 caseN is relatively small. In this note, the forgetting factor λ=1 is initially recommended, and the selection of N 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.


## 2.1. Kernelized TARMA Model and Maximum Correntropy Criterion

The TARMA model withna and nc designating, respectively, its autoregressive (AR) and moving average (MA) orders, is of the general form [12](1)xt+∑i=1naAitxt−i=et+∑i=1ncCitet−i,et∼NID0,Σt,with t designating the normalized discrete time, xt the discrete-time nonstationary vibration signal modeled, et an unobservable uncorrelated innovation sequence with zero mean and time-varying nonsingular covariance matrix Σt, and Ait,Cit the model’s time-dependent AR and MA parameter matrices. NID⋅ 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(2)Ait≈ω1,1AitTut…ω1,kAitTut⋮⋱⋮ωk,1AitTut…ωk,kAitTut,Cit≈ω1,1CitTut…ω1,kCitTut⋮⋱⋮ωk,1CitTut…ωk,kCitTut.The kernelized TARMA model is given by [24](3)xt=−∑i=1naAitxt−i+∑i=1ncCitet−i+et,⇔xt=wtTφt+et,⇔xt=w1,1A1tT…w1,kA1tT…w1,1CnctT…w1,kCnctT⋮⋱⋮⋱⋮⋱⋮wk,1A1tT…wk,kA1tT…wk,1CnctT…wk,kCnctT−xt−1⊗ut⋮et−nc⊗ut+et,⇔xt=wtφt⊗ut+et,⇔xt=wtTφt+et,where ⊗ is the Kronecker product, wt is the time-dependent parameter vector, φt is the regression matrix, and ωt and ϕt are the corresponding kernelized counterparts. ut is the transformed feature vector lying in the feature space; by using a Gaussian kernel, the “kernel trick” can be expressed as follows:(4)utiTutj=κkti,tj=exp−ti−tj22σk2,where σk designates the adjustable parameter.The system’s “frozen-time” characteristic matrix equation is given by [25](5)AB,t=I+∑i=1naAitBi=0,with B designating the backshift operator, Bixt=xt−i. The roots of the system can be derived from a general eigenvalue problem [26](6)Dt−λrtIVrt=0,where λrt and Vrt=λr−naLr,…,λr−1LrT stand for the rth eigenvalue and eigenvector of Dt. The matrix Dt is constructed from the AR parameter matrices as(7)Dt=0I…0⋮⋮⋱⋮00…I−Anat−Ana−1t…−A1t.The “frozen-time” modal frequencies and damping ratios of the system can be computed by(8)frt=lnλrt2πTs,ζrt=−Relnλrtlnλrt.In general, correntropy is a local similarity measure between two arbitrary random variables ofX,Y∈ℝ, with joint distribution function FXYx,y, and is defined by [27](9)VX,Y=EκcX,Y=∫κcx,ydFXYx,y,where κ⋅ is a kernel function that satisfies Mercer theory and E⋅ denotes the expectation operator. In this note, the shift-invariant Mercer kernel is the Gaussian kernel given by(10)κcx,y=Ge=exp−e22σc2,where e=x−y, and σc>0 stands for the kernel bandwidth.For a ridge regression problem, with parameter weight vectorw, the corresponding optimal problem based on the maximum correntropy criterion can be expressed as follows:(11)maxwJMCC=maxw∑i=1NEexp−e2i2σc2−γw2.


## 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 factor0<λ≤1 and the regularization parameter γ≥0, the weighted regularized maximum correntropy cost function is defined as(12)maxωt∑τ=t−N+1tλNκcxτ−ωtTϕτ−12λNγωt2,and the solution to this optimal problem is given by(13)ωt=Ψ⟶tΨ⟶tTΨ⟶t+λNγσcΓ⟶2t−1−1x⟶t,with(14)Ψ⟶t≜ϕt−N+1,ϕt−N+2,…,ϕt,x⟶t≜xt−N+1,xt−N+2,…,xtT,Γ⟶t≜diagλN−1Et−N+1,λN−2Et−N+2,…,Et,where Et≜exp−1/2σc2etTet.DenoteQ⟶t=Ψ⟶tTΨ⟶t+λNγσc2Γ⟶t−1 and k⟶t=ϕt−N+1,…,ϕt−1Tϕt; the following equation can be received:(15)Q⟶t=Q⟶∗t−1k⟶tk⟶tTϕtTϕt+λNγσc2Et−1,where Q⟶∗t−1 can be obtained by removing the first row and column of Q⟶t−1.By using “kernel trick,”ϕtTϕt can be computed as follows:(16)ϕtTϕt=ϕt⊗utTϕt⊗ut=ϕtTϕtκkt,t,where κkti,tj=exp−ti−tj2/2σk.By using the upsized matrix inversion formula, the update process ofQ⟶t−1 is(17)Q⟶t−1=Q⟶∗t−1−1+z⟶tz⟶tTr⟶t−1−z⟶trt−1−z⟶tTr⟶t−1r⟶t−1,with z⟶t=Q⟶∗t−1−1k⟶t and r⟶t=λNγσc2Et−1+ϕtTϕt−z⟶tTk⟶t.Define(18)Q⟶t−1−1=A⟶B⟶C⟶D⟶.By using the downsized matrix inversion formula, the following equation can be obtained:(19)Q⟶∗t−1−1=D⟶-C⟶A⟶−1B⟶.Considering equations (19) and (13), the kernelized parameter matrix ωt can be obtained:(20)ωt=Ψ⟶tQ⟶t−1x⟶t.Finally, time-dependent AR and MA parameter matrices can be computed by(21)wtT=A1t,…,Anat,C1t,…,Cnct=wtTIkna+nc⊗ut=Ψ⟶tQ⟶t−1x⟶tTIkna+nc⊗ut=Q⟶t−1x⟶tTφt−N+1T⊗κkt−N+1,t⋮φtT⊗κkt,t.The workflow of the SWRMCRR-KTARMA method is summarized in (Algorithm1).Algorithm 1: SWRMCRR-KTARMA method.
Initialization:w^0T=0,σc,σk,λ,γ,NComputation1≤t≤N:e^tt−1=xt−x^tt−1=xt−w^t−1TφtEt=exp−1/2σc2e^tt−1Te^tt−1kt=ψ˜t−1Tϕt=ϕ1Tϕt⋮ϕt−1Tϕt=φ1Tφtκk1,t⋮φt−1Tφtκkt−1,tzt=Qt−1−1ktrt=λtγσc2Et−1+φtTφtκkt,t−ztTktQt−1=Qt−1−1+ztztTrt−1−ztrt−1−ztTrt−1rt−1w^tT=Qt−1x˜tTφ1T⊗κk1,t⋮φtT⊗κkt,tComputationt>N:Downsizing: obtainQ⟶∗t−1 by removing the first row and column of Q⟶t−1.ComputeQ⟶∗t−1−1 according to the downsized matrix inversion formula.e^tt−1=xt−w^t−1TφtEt=exp−1/2σc2e^tt−1Te^tt−1k⟶t=Ψ⟶t−1Tϕt=ϕt−N+1Tϕt⋮ϕt−1Tϕt=ϕt−N+Tϕtκkt−N+1,t⋮ϕt−1Tϕtκkt−1,tz⟶t=Q⟶∗t−1−1ktr⟶t=λNγσc2Et−1+φtTφtκkt,t−z⟶tTk⟶tQ⟶t−1=Q⟶∗t−1−1+z⟶tz⟶tTr⟶t−1−z⟶tr⟶t−1−z⟶tTr⟶t−1r⟶t−1w^t=Q⟶t−1x⟶tTφt−N+1T⊗κkt−N+1,t⋮φtT⊗κkt,tTARMA model parameter estimation:Α1t,...,Αnat,C1t,...,Cnct=w^tT


## 2.3. Model Structure Selection

The model structure parameters of the proposed SWRMCRR-KTARMA method include the kernel bandwidthσc and σk, the forgetting factor λ, the sliding-window length N, 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(22)RSS≜∑t=1Ne^tTe^tN.Model structure selection of SWRMCRR-KTARMA can be tackled via the following three basic steps:Step 1: given initial values ofσc, σk, λ, N, 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 withna = nc, and the fitness AR order n^a is selected when the decrease of RSS values is quite subtle(b)
Given the integer search space0,n^a, and the best fit MA order n^c 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σk reflects the nonlinear power of the proposed method; it should be selected based on different problems. The kernel bandwidth σc 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 caseN is relatively small. In this note, the forgetting factor λ=1 is initially recommended, and the selection of N 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 Figure1. 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).Figure 1
Schematic diagram (a) and setup (b) of the experimental system.
(a)(b)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 Figure1. 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).Figure 2
Modal parameters varying with the water level position: (a) natural frequencies and (b) damping ratios.
(a)(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 given(23)Φt=expjμt−ηtα1+jβsigntωt,α,where(24)ωt,α=tanαπ2,α≠1,2πlogt,α=1,signt=1,t>0,0,t=0,−1,t<0,and α is the characteristic exponent satisfying 0<α≤2, which controls the thickness of tails in the distribution. μ is the location parameter −∞<μ<∞, which corresponds to the mean for 1<α≤2 and the median for 0<α≤1. η is the dispersion parameter (η>0), which determines the spread of the density around its location parameter. β is the symmetry parameter (−1≤β≤1); when β=0, the distribution is symmetric around the location parameter and referred to as the symmetry alpha-stable distribution (SαS). 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 as(25)MSNRdB=10⋅log10σsignal2η,where σsignal2 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 (α=0.9,β=0,andμ=0) is mainly considered.Figure 3
Probability density functions for alpha-stable distributions: (a) the overall shape and (b) the tails.
(a)(b)


### 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,σc=1,σk=105,N=175,λ=1,γ=3×10−3, 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 na=7,nc=1. The forgetting factor λ=1 is finally selected as the sliding-window length N=175 is relatively small. Also, the suggested AR, MA orders and forgetting factor of the EWRPLR-TARMA method are na=8, nc=0 and λ=0.990, 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.Figure 4
Model structure selection of the SWRMCRR-KTARMA method: (a) AR order and (b) MA order.
(a)(b)Figure 5
EWRPLR-TARMA model order selection: (a) AR order, (b) MA order, and (c) forgetting factor.
(a)(b)(c)Table 1
Identification methods, their characteristics, and the identification models.
Identification methodMethod characteristicsIdentified modelSWRMCRR-KTARMAσc=1,σk=105,N=175,λ=1,γ=3×10−3KTARMA (7, 1)EWRPLR-TARMAλ=0.990TARMA (8, 0)The modal frequency identification results of SWRMCRR-KTARMA and EWRPLR-TARMA can be seen in Figure6. 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.Figure 6
Modal frequencies estimated by SWRMCRR-KTARMA and EWRPLR-TARMA (green lines: baseline values; red dots: modal frequency estimates of all 10 Monte Carlo runs).In order to compare these methods by using natural frequency estimates in a quantitative manner, a mean absolute error (MAE) is introduced as(26)MAE=1R∑i=1R1Nei∑t=1Neift−f^it.with ft designating the baseline value of the modal frequency, f^it its estimated value based on the ith Monte Carlo run, and Nei 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 Figures7 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.Figure 7
The MAE values of modal frequency estimates based on the experimental response signal with various added numerical noises for all 10 Monte Carlo runs.Figure 8
The average effective number of modal frequency estimates based on the experimental response signal with added numerical noise for all 10 Monte Carlo runs.The time and memory complexities of the SWRMCRR-KTARMA method areON2, which is mainly determined by the sliding-window length. However, the time and memory complexities of the EWRPLR-TARMA method are Ok5na+nc2 [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 k) is quite large which is common in the large aerospace structure test.Figure 9
Monte Carlo comparison of different methods: normalized processing CPU times required for parameter estimation.


## 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 Figure1. 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).Figure 1
Schematic diagram (a) and setup (b) of the experimental system.
(a)(b)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 Figure1. 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).Figure 2
Modal parameters varying with the water level position: (a) natural frequencies and (b) damping ratios.
(a)(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 given(23)Φt=expjμt−ηtα1+jβsigntωt,α,where(24)ωt,α=tanαπ2,α≠1,2πlogt,α=1,signt=1,t>0,0,t=0,−1,t<0,and α is the characteristic exponent satisfying 0<α≤2, which controls the thickness of tails in the distribution. μ is the location parameter −∞<μ<∞, which corresponds to the mean for 1<α≤2 and the median for 0<α≤1. η is the dispersion parameter (η>0), which determines the spread of the density around its location parameter. β is the symmetry parameter (−1≤β≤1); when β=0, the distribution is symmetric around the location parameter and referred to as the symmetry alpha-stable distribution (SαS). 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 as(25)MSNRdB=10⋅log10σsignal2η,where σsignal2 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 (α=0.9,β=0,andμ=0) is mainly considered.Figure 3
Probability density functions for alpha-stable distributions: (a) the overall shape and (b) the tails.
(a)(b)


## 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,σc=1,σk=105,N=175,λ=1,γ=3×10−3, 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 na=7,nc=1. The forgetting factor λ=1 is finally selected as the sliding-window length N=175 is relatively small. Also, the suggested AR, MA orders and forgetting factor of the EWRPLR-TARMA method are na=8, nc=0 and λ=0.990, 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.Figure 4
Model structure selection of the SWRMCRR-KTARMA method: (a) AR order and (b) MA order.
(a)(b)Figure 5
EWRPLR-TARMA model order selection: (a) AR order, (b) MA order, and (c) forgetting factor.
(a)(b)(c)Table 1
Identification methods, their characteristics, and the identification models.
Identification methodMethod characteristicsIdentified modelSWRMCRR-KTARMAσc=1,σk=105,N=175,λ=1,γ=3×10−3KTARMA (7, 1)EWRPLR-TARMAλ=0.990TARMA (8, 0)The modal frequency identification results of SWRMCRR-KTARMA and EWRPLR-TARMA can be seen in Figure6. 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.Figure 6
Modal frequencies estimated by SWRMCRR-KTARMA and EWRPLR-TARMA (green lines: baseline values; red dots: modal frequency estimates of all 10 Monte Carlo runs).In order to compare these methods by using natural frequency estimates in a quantitative manner, a mean absolute error (MAE) is introduced as(26)MAE=1R∑i=1R1Nei∑t=1Neift−f^it.with ft designating the baseline value of the modal frequency, f^it its estimated value based on the ith Monte Carlo run, and Nei 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 Figures7 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.Figure 7
The MAE values of modal frequency estimates based on the experimental response signal with various added numerical noises for all 10 Monte Carlo runs.Figure 8
The average effective number of modal frequency estimates based on the experimental response signal with added numerical noise for all 10 Monte Carlo runs.The time and memory complexities of the SWRMCRR-KTARMA method areON2, which is mainly determined by the sliding-window length. However, the time and memory complexities of the EWRPLR-TARMA method are Ok5na+nc2 [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 k) is quite large which is common in the large aerospace structure test.Figure 9
Monte Carlo comparison of different methods: normalized processing CPU times required for parameter estimation.


## 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.



---

*Source: 2946709-2020-06-05.xml*