# An Extended Time Series Algorithm for Modal Identification from Nonstationary Ambient Response Data Only

Authors: Chang-Sheng Lin; Dar-Yun Chiang; Tse-Chuan Tseng
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

Modal Identification is considered from response data of structural systems under nonstationary ambient vibration. The conventional autoregressive moving average (ARMA) algorithm is applicable to perform modal identification, however, only for stationary-process vibration. The ergodicity postulate which has been conventionally employed for stationary processes is no longer valid in the case of nonstationary analysis. The objective of this paper is therefore to develop modal-identification techniques based on the nonstationary time series for linear systems subjected to nonstationary ambient excitation. Nonstationary ARMA model with time-varying parameters is considered because of its capability of resolving general nonstationary problems. The parameters of moving averaging (MA) model in the nonstationary time-series algorithm are treated as functions of time and may be represented by a linear combination of base functions and therefore can be used to solve the identification problem of time-varying parameters. Numerical simulations confirm the validity of the proposed modal-identification method from nonstationary ambient response data.


---

## Body

## 1. Introduction

Experimental identification of modal parameters of a structure is usually carried out by measuring both its input and its corresponding output. Some modal testing techniques use free or impulse responses of structures so that the input excitation need not be measured. However, there are situations where controlled excitation or initial excitation cannot be employed, such as the case of in-operation testing or in-flight measurement. Consequently, it is desirable to develop techniques for modal-parameter identification without the need of input measurement [1, 2].Modal-parameter identification from ambient vibration data has gained considerable attention in recent years [3, 4]. The feature of experimental modal analysis lies in the emphasis on the importance of the modal-identification process of a structure from measured input/output data when compared with the classic vibration analysis. Normally, the analysis of system identification acquires the estimation of important parameters from the data measured in the experiments. Based on the historical development of time series, Yule was the first to purpose the autoregressive (AR) model of a single variable in 1927. He developed the application of AR model for investigating periodicities in disturbed series, with special reference to Wolfer’s sunspot numbers. In 1970, Box et al. [5] proposed a modeling cycle for the AR model, which assumes that future values of a time series can be expressed as a linear combination of its past values. Applications of the Box et al. methodology spread in the following decades, covering a wide range of scientific areas such as Biology, Astronomy, or Econometrics. The time-series algorithm is also extensively used in the modal identification [6, 7]. In 1980, Pandit and Wu [8] proposed a time-series method based on the mathematical theory of statistics and further developed the so-called data dependent system (DDS) using ARMA models for modal-parameter identification of structural systems [9, 10]. The main advantages of DDS are that they are accurate for extracting information from noisy signals, that estimation of the parameter covariance matrix is natural part of the estimation, and that they can be formulated as exact covariance equivalent discrete time domain models of a set of second order differential equations.The conventional ARMA algorithm is applicable to perform modal identification, however, only for stationary-process vibration. In this paper, we will develop modal-identification techniques based on the nonstationary time series for linear systems subjected to nonstationary ambient excitation. Nonstationary ARMA model with time-varying parameters is considered because of its capability of resolving general nonstationary problems. The parameters of MA model in the nonstationary time-series algorithm are treated as functions of time and may be represented by a linear combination of base functions and therefore can be used to solve the identification problem of time-varying parameters. Numerical simulations confirm the validity of the proposed modal-identification method from nonstationary ambient response data.


## 2. Time-Series Model of Structural Vibration Responses

The standard matrix equations of motion of a discrete linear system can be expressed as follows:(1)Mx¨(t)+Dx˙(t)+Kx(t)=f(t).
In (1), M, D, and K are,  respectively,  the system mass, damping, and stiffness matrices, while f(t) and x(t) are the excitation and displacement vectors, respectively. Applying Laplace transform to (1) yields
(2)X(s)=Z(s)-1F(s)=H(s)F(s),
where Z(s)=Ms2+Ds+K. Z(s) is called the system impedance matrix, which is an inverse matrix of the transfer function matrix H(s). According to adjoint matrix method, the transfer function H(s) can be expressed as
(3)H(s)=Z(s)-1=adjZ(s)det⁡|Z(s)|,
where adjZ(s) is the adjoint matrix of Z(s). It follows from (3) that the transfer function Hlp(s) at the lth degree of freedom (DOF) due to the input at the pth DOF is
(4)Hlp(s)=b0+b1s+b2s2+⋯+bmsma0+a1s+a2s2+⋯+ansn,
where m=2N-2, n=2N. For a real structure vibration system with damping, the transfer function, (4), can be decomposed into
(5)Hlp(s)=∑r=1N(Alprs-sr+Alpr*s-sr*),
where sr is a pole corresponding to the rth mode of a system and Alpr is the residue in the form of a complex number. To obtain the time-dependent impulse response hlp(t), it is, therefore, necessary to perform an inverse Laplace transformation of Hlp(s) as follows:
(6)hlp(t)=∑r=1N(Alpresrt+Alpr*esr*t).
Define t=kΔ as the discrete time step, and (6) can be written as
(7)hlp(k)=∑r=1N(AlpresrkΔ+Alpr*esr*kΔ),
where k is the time step and Δ is the sampling interval. Let Zr=esrΔ be a Z-transform factor; (7) can be rewritten as
(8)hlp(k)=∑r=1N(AlprZrk+Alpr*Zr*k).
The transfer function Hlp(Z) of the discrete system results from the Z-transform of impulse response function hlp(k) as follows:
(9)Hlp(Z)=∑r=1N(Alpr1-ZrZ-1+Alpr*1-Zr*Z-1).
Comparing the relationship between input and output in original transfer function, (9) can be further expressed as follows:
(10)Hlp(Z)=θ0+θ1Z-1+θ2Z-2+⋯+θ2N-1Z-2N+11+φ1Z-1+φ2Z-2+⋯+φ2NZ-2N=X(Z)a(Z).
One can, therefore, derive
(11)(1+φ1Z-1+φ2Z-2+⋯+φ2NZ-2N)X(Z)=(θ0+θ1Z-1+θ2Z-2+⋯+θ2N-1Z-2N+1)a(Z).
By inverse transforming both sides of the Z-domain (11) and rearranging, we obtain the following equation in the time domain:
(12)Xt+φ1Xt-1+φ2Xt-2+⋯+φ2NXt-2N=at+θ1at-1+θ2at-2+⋯+θ2N-1at-2N+1.
Equation (12) is a normal ARMA (2N,2N-1) model. The poles of transfer function can then be obtained by solving the polynomial equation a(Z) in (10) as follows:
(13)∑k=02NφkZ-k=1+φ1Z-1+φ2Z-2+⋯+φ2NZ-2N=0.
The relationship between the poles, Zr and Zr*, of transfer function as well as natural frequencies ωr and damping ratios ζr of a structural system can be expressed as follows:
(14)Zr=esrΔ=e(-ζrωr+jωr1-ζr2)Δ,Zr*=esr*Δ=e(-ζrωr-jωr1-ζr2)Δ
from which the natural frequencies ωr and damping ratios ζr of a structural system can then be obtained to be
(15)ωr=1Δln⁡Zrln⁡Zr*,ζr=-ln⁡(ZrZr*)2ln⁡Zrln⁡Zr*.The preceding results show that the ARMA(2N,2N-1) model is applicable to identify the modal parameters of an N-DOF structural system. However, the original ARMA model was effective only for stationary process. Remarkably, the relationship between AR and MA parts in ARMA model is approximated as a relationship between input and output. In the light of this concept, we come up with a new idea that if the nonstationary ambient vibration data can also be approximated directly as the aforementioned relationship of input and output, then the ARMA model can be extended for nonstationary process. In the following, we will propose a modification of ARMA model using the aforementioned concept to construct a nonstationary ARMA model, which can be applied to modal identification from nonstationary ambient vibration data.


## 3. Estimation of ARMA Parameters

The ARMA algorithm is based on the assumption of a stationary input while the others assume that the input is a Gaussian white noise process. However, estimating the parameters is more involved for the ARMA model than it is for the AR model because the white noise signal in the ARMA model is also unknown, so a three-stage least-square approach or nonlinear approach has to be proposed. As mentioned in the previous section, it is important to estimate the model parameter for reconstructing the response data of a system through the ARMA algorithm. In this paper, we estimate the parameters of ARMA by using the nonlinear least-square method. However, it requires a good initial estimation as the initial parameters of the ARMA model when performing the nonlinear least-square method. We therefore propose a three-stage least-square method, which is applicable for nonstationary process, to obtain the initial estimated parameters of ARMA model.


### 3.1. Three-Stage Least-Square Method for the Estimation of Initial Parameters of ARMA Model

Consider an ARMA model as follows:(16)Wt+φ1Wt-1+φ2Wt-2+⋯+φpWt-p=at+θ1at-1+θ2at-2+⋯+θqat-q,
where at is white noise, p and q are the orders of an ARMA model, and φ1,φ2,…, and φp as well as θ1,θ2,…, and θq are the parameters of ARMA model. By treating the linear combination of at-q in the right-hand side of (16) as new error term at′, (16) can be rewritten as
(17)Wt+φ1Wt-1+φ2Wt-2+⋯+φpWt-p=at′.
Equation (17) can be treated as a p-order AR model, whose parameters can be determined through the least-square method as follows:
(18)y_=[X]α_,
where
(19)y_=[WkWk+1⋮Wn],[X]=[-Wk-1-Wk-2⋯-Wk-p-Wk-Wk-1⋯-Wk+1-p⋮⋮⋯⋮-Wn-1-Wn-2⋯-Wn-p],α_=[φ1φ2⋮φp],kkkkkkkkkkkkkkkkkkkkkkkkk=p+1.
The estimations α^_ can then be obtained as follows:
(20)α^_=[[X]T[X]]-1[X]Ty_.
A new set of estimations Wt′ can be composed of the obtained AR parameters, φ1,φ2,…, and φp, as follows:
(21)Wt′=Wt+φ1Wt-1+φ2Wt-2+⋯+φpWt-p.
Substituting the new observability data into the AR part in the original ARMA model, the following equation can be derived:
(22)Wt′=at+θ1at-1+θ2at-2+⋯+θqat-qorWt′=Θ(B)at;Θ(B)=1+θ1B+θ2B2+⋯+θqBq.
Equation (22) can be treated as q-order MA model. Due to the fact that the finite-order MA model can be approximated as the high-order AR model [9], (22) can be rewritten as
(23)at=Θ-1(B)Wt′=(∑j=0∞φjBj)Wt′≈Φ(B)Wt′,
where Φ(B)=1+φ1′B+φ2′B2+⋯+φnp′Bnp and np is the order of high-order AR model. The parameters of high-order AR model can also be obtained through least-square method. Comparing (22) with (23), one can derive the following:
(24)Φ(B)Wt′=at,Wt′=Θ(B)at,
from which the relationship between Θ(B) and Φ(B) can be obtained to be
(25)Θ(B)Φ(B)=1.
The above equation can be rewritten as the standard form of least-square method, and the parameters of MA model are further determined as follows:
(26)Y_=[X]μ_,
where
(27)Y_=[-φ1′-φ2′⋮-φnp′0⋮0]np+q,[X]=[10⋯⋯0φ1′10⋯0⋮⋮⋮⋮⋮φq-1′⋮⋮⋮1⋮⋮⋮⋮⋮φnp′⋮⋮⋮⋮0⋮⋮⋮⋮⋮⋮⋮⋮⋮0⋯⋯0φnp′](np+q)×q,μ_=[θ1θ2⋮θq]q.
The estimations μ^_ can then be obtained as follows:
(28)μ^_=[[X]T[X]]-1[X]TY_.
In the preceding, it has been shown that the initial estimations of all parameters of ARMA model can be obtained through the three-stage least-square method without the assumption of stationary process. To obtain the relatively precise estimation of ARMA model, the nonlinear least-square method in conjunction with the Gauss-Newton algorithm will be employed, as described in detail next.


### 3.2. Nonlinear Least-Square Method for Estimation of Parameters of ARMA Model

Consider an ARMA model, which can be expressed as(29)Wt=-φ1Wt-1-φ2Wt-2-⋯-φpWt-p+θ1at-1+⋯+θqat-q+at,
where p and q are the orders of an ARMA model. φi, i=1,2,…,p, and θi, i=1,2,…,q, are constants, which are parameters of an ARMA model, as well as at being white noise. It follows from (29) that at can be derived as follows:
(30)at=Wt+φ1Wt-1+φ2Wt-2+⋯+φpWt-p-θ1at-1-⋯-θqat-q.
By integrating the terms of parameters of ARMA, (29) can be rewritten as
(31)Wt=fi(φi,θi,Wi,ai)+at.
Let unknown at be a set of the initial estimation of model parameters in conjunction with (30), whose vector form is as follows:
(32)Y=F(W,β)+e.
In (32), F(W,β) is a function of W:wk,wk+1,…,wN and β={φi,i=1,2,…,p;θj,j=1,2,…,q}, where Y=[Wk,Wk+1,…,Wn]T, F(W,β)=[fk,fk+1,…,fn]T, e=[ak,ak+1,⋯,an]T, k=p+1, and N is the number of samples in the input (or output) data. Define the prediction errors Q as follows:
(33)Q=[Y-F(W,β)]T[Y-F(W,β)].
Equation (30) is necessarily to be employed to construct F(W,β) due to the fact that at is unknown. F(W,β) will in general be nonlinear; that is, it is not the linear function of β, and the minimization of Q cannot be performed. Thus, to obtain a practically applicable approach, through the Gauss-Newton algorithm, which is a method usually applied to solve nonlinear least squares problems, F(W,β) can be linearized using a first-order Taylor expansion at the operating point β=β0 and then extract the linear part of F(W,β) as follows:
(34)F(W,β)=[fp+1(β1,β2,…,βm)fp+2(β1,β2,…,βm)⋮fN(β1,β2,…,βm)]=[fp+1(β10,β20,…,βm0)fp+2(β10,β20,…,βm0)⋮fN(β10,β20,…,βm0)]+[∂fp+1∂β1∂fp+1∂β2⋯∂fp+1∂βm∂fp+2∂β1∂fp+2∂β2⋯∂fp+2∂βm⋮∂fN∂β1∂fN∂β2⋯∂fN∂βm][β1-β10β2-β20⋮βm-βm0],
where β0 is an initial estimation set of ARMA parameters and can be expressed as
(35)β0=[β10,β20,…,βm0],m=p+q.
Due to the fact that W is a known observability data matrix when performing the Gauss-Newton algorithm and to emphasize that F(W,β) is function of β, F(W,β) can thus be written as F(β). Equation (34) can be therefore expressed as
(36)F(β)=F(β0)+X(β0)(β-β0).
Substituting (36) into (34), the following equation can be derived:
(37)Q=[Y-F(β0)-X(β0)(β-β0)]T×[Y-F(β0)-X(β0)(β-β0)].
According to the theory of least-square method, through the minimization of Q, the estimated parameter β^ can be expressed as
(38)β^=β0+[XT(β0)X(β0)]-1XT(β0)[Y-F(β0)].
In practice, the evaluation result of β^ is an iterative form of running to a set of initial estimation β0 for the beginning. Through the iteration evaluation of a set of β^(k), β^(k+1) will be yielded better and then continue iterations until the final convergence of the results. The iteration form can be represented by
(39)β^(k+1)=β^(k)+[X(k)TX(k)]-1X(k)T(Y-F(k)),
and the criterion of the convergence of calculation results to stop iterations can be defined as follows:
(40)|∥β^(k+1)∥-∥β^(k)∥∥β^(k)∥|<δ,
where ∥*∥ is norm of parametric vector and δ is small value to be chosen and may be chosen as 10-4 or 10-5.In the preceding, it has been shown that the initial estimations of all parameters of ARMA model can be obtained through the three-stage least-square method without the assumption of stationary process. By using the nonlinear least-square method in conjunction with the Gauss-Newton algorithm, the relatively precise estimation of ARMA model can then be obtained. In the following, we will reconstruct a nonstationary ARMA algorithm and then further extend the time-series algorithm to general nonstationary ambient vibration problems.


## 3.1. Three-Stage Least-Square Method for the Estimation of Initial Parameters of ARMA Model

Consider an ARMA model as follows:(16)Wt+φ1Wt-1+φ2Wt-2+⋯+φpWt-p=at+θ1at-1+θ2at-2+⋯+θqat-q,
where at is white noise, p and q are the orders of an ARMA model, and φ1,φ2,…, and φp as well as θ1,θ2,…, and θq are the parameters of ARMA model. By treating the linear combination of at-q in the right-hand side of (16) as new error term at′, (16) can be rewritten as
(17)Wt+φ1Wt-1+φ2Wt-2+⋯+φpWt-p=at′.
Equation (17) can be treated as a p-order AR model, whose parameters can be determined through the least-square method as follows:
(18)y_=[X]α_,
where
(19)y_=[WkWk+1⋮Wn],[X]=[-Wk-1-Wk-2⋯-Wk-p-Wk-Wk-1⋯-Wk+1-p⋮⋮⋯⋮-Wn-1-Wn-2⋯-Wn-p],α_=[φ1φ2⋮φp],kkkkkkkkkkkkkkkkkkkkkkkkk=p+1.
The estimations α^_ can then be obtained as follows:
(20)α^_=[[X]T[X]]-1[X]Ty_.
A new set of estimations Wt′ can be composed of the obtained AR parameters, φ1,φ2,…, and φp, as follows:
(21)Wt′=Wt+φ1Wt-1+φ2Wt-2+⋯+φpWt-p.
Substituting the new observability data into the AR part in the original ARMA model, the following equation can be derived:
(22)Wt′=at+θ1at-1+θ2at-2+⋯+θqat-qorWt′=Θ(B)at;Θ(B)=1+θ1B+θ2B2+⋯+θqBq.
Equation (22) can be treated as q-order MA model. Due to the fact that the finite-order MA model can be approximated as the high-order AR model [9], (22) can be rewritten as
(23)at=Θ-1(B)Wt′=(∑j=0∞φjBj)Wt′≈Φ(B)Wt′,
where Φ(B)=1+φ1′B+φ2′B2+⋯+φnp′Bnp and np is the order of high-order AR model. The parameters of high-order AR model can also be obtained through least-square method. Comparing (22) with (23), one can derive the following:
(24)Φ(B)Wt′=at,Wt′=Θ(B)at,
from which the relationship between Θ(B) and Φ(B) can be obtained to be
(25)Θ(B)Φ(B)=1.
The above equation can be rewritten as the standard form of least-square method, and the parameters of MA model are further determined as follows:
(26)Y_=[X]μ_,
where
(27)Y_=[-φ1′-φ2′⋮-φnp′0⋮0]np+q,[X]=[10⋯⋯0φ1′10⋯0⋮⋮⋮⋮⋮φq-1′⋮⋮⋮1⋮⋮⋮⋮⋮φnp′⋮⋮⋮⋮0⋮⋮⋮⋮⋮⋮⋮⋮⋮0⋯⋯0φnp′](np+q)×q,μ_=[θ1θ2⋮θq]q.
The estimations μ^_ can then be obtained as follows:
(28)μ^_=[[X]T[X]]-1[X]TY_.
In the preceding, it has been shown that the initial estimations of all parameters of ARMA model can be obtained through the three-stage least-square method without the assumption of stationary process. To obtain the relatively precise estimation of ARMA model, the nonlinear least-square method in conjunction with the Gauss-Newton algorithm will be employed, as described in detail next.


## 3.2. Nonlinear Least-Square Method for Estimation of Parameters of ARMA Model

Consider an ARMA model, which can be expressed as(29)Wt=-φ1Wt-1-φ2Wt-2-⋯-φpWt-p+θ1at-1+⋯+θqat-q+at,
where p and q are the orders of an ARMA model. φi, i=1,2,…,p, and θi, i=1,2,…,q, are constants, which are parameters of an ARMA model, as well as at being white noise. It follows from (29) that at can be derived as follows:
(30)at=Wt+φ1Wt-1+φ2Wt-2+⋯+φpWt-p-θ1at-1-⋯-θqat-q.
By integrating the terms of parameters of ARMA, (29) can be rewritten as
(31)Wt=fi(φi,θi,Wi,ai)+at.
Let unknown at be a set of the initial estimation of model parameters in conjunction with (30), whose vector form is as follows:
(32)Y=F(W,β)+e.
In (32), F(W,β) is a function of W:wk,wk+1,…,wN and β={φi,i=1,2,…,p;θj,j=1,2,…,q}, where Y=[Wk,Wk+1,…,Wn]T, F(W,β)=[fk,fk+1,…,fn]T, e=[ak,ak+1,⋯,an]T, k=p+1, and N is the number of samples in the input (or output) data. Define the prediction errors Q as follows:
(33)Q=[Y-F(W,β)]T[Y-F(W,β)].
Equation (30) is necessarily to be employed to construct F(W,β) due to the fact that at is unknown. F(W,β) will in general be nonlinear; that is, it is not the linear function of β, and the minimization of Q cannot be performed. Thus, to obtain a practically applicable approach, through the Gauss-Newton algorithm, which is a method usually applied to solve nonlinear least squares problems, F(W,β) can be linearized using a first-order Taylor expansion at the operating point β=β0 and then extract the linear part of F(W,β) as follows:
(34)F(W,β)=[fp+1(β1,β2,…,βm)fp+2(β1,β2,…,βm)⋮fN(β1,β2,…,βm)]=[fp+1(β10,β20,…,βm0)fp+2(β10,β20,…,βm0)⋮fN(β10,β20,…,βm0)]+[∂fp+1∂β1∂fp+1∂β2⋯∂fp+1∂βm∂fp+2∂β1∂fp+2∂β2⋯∂fp+2∂βm⋮∂fN∂β1∂fN∂β2⋯∂fN∂βm][β1-β10β2-β20⋮βm-βm0],
where β0 is an initial estimation set of ARMA parameters and can be expressed as
(35)β0=[β10,β20,…,βm0],m=p+q.
Due to the fact that W is a known observability data matrix when performing the Gauss-Newton algorithm and to emphasize that F(W,β) is function of β, F(W,β) can thus be written as F(β). Equation (34) can be therefore expressed as
(36)F(β)=F(β0)+X(β0)(β-β0).
Substituting (36) into (34), the following equation can be derived:
(37)Q=[Y-F(β0)-X(β0)(β-β0)]T×[Y-F(β0)-X(β0)(β-β0)].
According to the theory of least-square method, through the minimization of Q, the estimated parameter β^ can be expressed as
(38)β^=β0+[XT(β0)X(β0)]-1XT(β0)[Y-F(β0)].
In practice, the evaluation result of β^ is an iterative form of running to a set of initial estimation β0 for the beginning. Through the iteration evaluation of a set of β^(k), β^(k+1) will be yielded better and then continue iterations until the final convergence of the results. The iteration form can be represented by
(39)β^(k+1)=β^(k)+[X(k)TX(k)]-1X(k)T(Y-F(k)),
and the criterion of the convergence of calculation results to stop iterations can be defined as follows:
(40)|∥β^(k+1)∥-∥β^(k)∥∥β^(k)∥|<δ,
where ∥*∥ is norm of parametric vector and δ is small value to be chosen and may be chosen as 10-4 or 10-5.In the preceding, it has been shown that the initial estimations of all parameters of ARMA model can be obtained through the three-stage least-square method without the assumption of stationary process. By using the nonlinear least-square method in conjunction with the Gauss-Newton algorithm, the relatively precise estimation of ARMA model can then be obtained. In the following, we will reconstruct a nonstationary ARMA algorithm and then further extend the time-series algorithm to general nonstationary ambient vibration problems.


## 4. Nonstationary Time-Series Model for Modal Identification

The original time-series model is the well-approximately equivalent mathematical model for stationary process and is further available in practice. For vibration problems in engineering, under the assumption of stationary process, we can construct a mathematical model for the response of a system through the ARMA model. By extracting the AR parameters, the parameters of a system can then be obtained. We also treated the MA part as an equivalent excitation force for the vibration process. The aforementioned time-series algorithm was effective only for stationary process. Most ambient excitations encountered in engineering problems are, however, nonstationary in nature; that is, the characteristics of nonstationary stochastic process may change with time. The original time-series algorithm is therefore not applicable to identify the modal parameters of a structural system subjected to nonstationary ambient excitation.In the following, we will improve the time-series algorithm to get rid of the previous restriction of stationary assumption. Without changing the AR parameters, we construct a nonstationary time-series model to effectively describe a wide variety of nonstationary behavior and then extract the modal parameters of a structural system. Furthermore, by introducing the base functions to construct a nonstationary time-series model, we extend the time-series algorithm to identify the modal parameters of a structural system subjected to nonstationary white noise.


### 4.1. Nonstationary ARMA Algorithm Using the Temporal Root-Mean-Square Function from Structural Response Data

The conventional stationary ARMA model, whose parameters are all constants, can be expressed as(41)Wt+φ1Wt-1+φ2Wt-2+⋯+φpWt-p=at+θ1at-1+θ2at-2+⋯+θqat-q.
In a previous study of the authors [4], it has been shown that if the excitation can be modeled as nonstationary white noise with a slowly time-varying envelope function, then the nonstationary responses of the system can also be modeled approximately as a product model with the same envelope function. In addition, if the original nonstationary response data could be represented by the product model with a slowly varying envelope function, the temporal root-mean-square functions of the data also have the same nonstationary trend as that of the original response data. Based on the above mentioned concept and assumption, in this paper, we first determine the temporal root-mean-square function σt and so the envelope function by using the interval average and then applying curve-fitting technique [4]. We can then acquire the approximate stationary responses by dividing the nonstationary responses of each DOF with the same σt, and a new set of observability data as quasistationary response data can then be obtained. The parameters of ARMA model can then be determined through the proposed algorithm in this paper. Define the nonstationary white noise in the product model by using σt as follows:
(42)at′=σt·at.
By substituting at′ into at in (41), the nonstationary time-series model can be constructed as follows:
(43)Wt+φ1Wt-1+φ2Wt-2+⋯+φpWt-p=at′+θ1at-1′+θ2at-2′+⋯+θqat-q′.
It indicates that the right-hand side of (43) and so the MA part in the ARMA model describes a nonstationary behavior of time-varying amplitude. The parameters of AR part (in the left-hand side of (43)) are all constants and can then be employed to extract the modal parameters of a structural system using (15).


### 4.2. Nonstationary ARMA Algorithm Using the Base Functions

The original ARMA model can be expressed as follows:(44)Wt+φ1Wt-1+φ2Wt-2+⋯+φpWt-p=at+θ1at-1+θ2at-2+⋯+θqat-q,
where the parameters of the original ARMA model are constants due to the stationary assumption. To effectively describe a time-varying stochastic characteristic of nonstationary behavior, we propose the general nonstationary time-series model, in which the parameters of model are described as function of time as follows:
(45)Wt+φ1(t)Wt-1+φ2(t)Wt-2+⋯+φp(t)Wt-p=at+θ1(t)at-1+θ2(t)at-2+⋯+θq(t)at-q.Consider a time-invariant linear dynamic system subjected to the nonstationary excitation; the MA part can be treated as an equivalent excitation force for the vibration process, and (45) can be simplified as
(46)Wt+φ1Wt-1+φ2Wt-2+⋯+φpWt-p=at+θ1(t)at-1+θ2(t)at-2+⋯+θq(t)at-q.
The concept of this nonstationary time-series model describes a wide variety of nonstationary excitation by using the parametric function of MA part only, and the parameters of AR model are constants to determine the modal parameters of time-invariant system. Due to the fact that the parametric functions θq(t)s in this model are all known ones, the general estimation method of model parameters cannot be employed. In this paper, we further come up with a new idea by using the known base function gj(t) in conjunction with the combination of unknown coefficient kji instead of the original unknown parametric function θq(t) in (46) as
(47)θi(t)=∑j=0mkji·gj(t);i=1,2,…,q.
The original ARMA model can then be rewritten as
(48)Wt+φ1Wt-1+φ2Wt-2+⋯+φpWt-p=at+[∑j=0m(gj(t)kj1)]at-1+[∑j=0m(gj(t)kj2)]at-2+⋯+[∑j=0m(gj(t)kjq)]at-q.
Due to the fact that the base function gj(t) is known, we can transfer the original problem to solve parametric function into the problem to solve the linear combination of at-q. The choice of the base function depends on the adequate expression to nonstationary response data. In this paper, we choose the polynomial function as base function as follows:
(49)gi(t)=(tN)i,i=0,1,…,m,
where N is the number of samples in the input (or output) data and i is the order of base function to be chosen. When introducing the base function to construct a nonstationary time-series model, the reconstructed ARMA model form is similar to that of the original time-series model. The major difference is that the MA parameters in the nonstationary base-function ARMA model are treated as functions of time and may be represented by a linear combination of suitable base functions. As long as by properly expanding the data matrix in the nonlinear and three-stage least-square method, respectively, the AR parameters and the coefficients of MA part corresponding to base function can be obtained. Then we construct nonstationary base-function ARMA model to describe the behavior of the nonstationary response data and further determine the modal parameters of a structural system. Note that although the structural modal parameters can be identified from the AR parameters only, the estimation of AR parameters and MA parametric function will affect each other in practice; estimation for AR parameters and MA parametric function cannot therefore be separately and solely performed.


## 4.1. Nonstationary ARMA Algorithm Using the Temporal Root-Mean-Square Function from Structural Response Data

The conventional stationary ARMA model, whose parameters are all constants, can be expressed as(41)Wt+φ1Wt-1+φ2Wt-2+⋯+φpWt-p=at+θ1at-1+θ2at-2+⋯+θqat-q.
In a previous study of the authors [4], it has been shown that if the excitation can be modeled as nonstationary white noise with a slowly time-varying envelope function, then the nonstationary responses of the system can also be modeled approximately as a product model with the same envelope function. In addition, if the original nonstationary response data could be represented by the product model with a slowly varying envelope function, the temporal root-mean-square functions of the data also have the same nonstationary trend as that of the original response data. Based on the above mentioned concept and assumption, in this paper, we first determine the temporal root-mean-square function σt and so the envelope function by using the interval average and then applying curve-fitting technique [4]. We can then acquire the approximate stationary responses by dividing the nonstationary responses of each DOF with the same σt, and a new set of observability data as quasistationary response data can then be obtained. The parameters of ARMA model can then be determined through the proposed algorithm in this paper. Define the nonstationary white noise in the product model by using σt as follows:
(42)at′=σt·at.
By substituting at′ into at in (41), the nonstationary time-series model can be constructed as follows:
(43)Wt+φ1Wt-1+φ2Wt-2+⋯+φpWt-p=at′+θ1at-1′+θ2at-2′+⋯+θqat-q′.
It indicates that the right-hand side of (43) and so the MA part in the ARMA model describes a nonstationary behavior of time-varying amplitude. The parameters of AR part (in the left-hand side of (43)) are all constants and can then be employed to extract the modal parameters of a structural system using (15).


## 4.2. Nonstationary ARMA Algorithm Using the Base Functions

The original ARMA model can be expressed as follows:(44)Wt+φ1Wt-1+φ2Wt-2+⋯+φpWt-p=at+θ1at-1+θ2at-2+⋯+θqat-q,
where the parameters of the original ARMA model are constants due to the stationary assumption. To effectively describe a time-varying stochastic characteristic of nonstationary behavior, we propose the general nonstationary time-series model, in which the parameters of model are described as function of time as follows:
(45)Wt+φ1(t)Wt-1+φ2(t)Wt-2+⋯+φp(t)Wt-p=at+θ1(t)at-1+θ2(t)at-2+⋯+θq(t)at-q.Consider a time-invariant linear dynamic system subjected to the nonstationary excitation; the MA part can be treated as an equivalent excitation force for the vibration process, and (45) can be simplified as
(46)Wt+φ1Wt-1+φ2Wt-2+⋯+φpWt-p=at+θ1(t)at-1+θ2(t)at-2+⋯+θq(t)at-q.
The concept of this nonstationary time-series model describes a wide variety of nonstationary excitation by using the parametric function of MA part only, and the parameters of AR model are constants to determine the modal parameters of time-invariant system. Due to the fact that the parametric functions θq(t)s in this model are all known ones, the general estimation method of model parameters cannot be employed. In this paper, we further come up with a new idea by using the known base function gj(t) in conjunction with the combination of unknown coefficient kji instead of the original unknown parametric function θq(t) in (46) as
(47)θi(t)=∑j=0mkji·gj(t);i=1,2,…,q.
The original ARMA model can then be rewritten as
(48)Wt+φ1Wt-1+φ2Wt-2+⋯+φpWt-p=at+[∑j=0m(gj(t)kj1)]at-1+[∑j=0m(gj(t)kj2)]at-2+⋯+[∑j=0m(gj(t)kjq)]at-q.
Due to the fact that the base function gj(t) is known, we can transfer the original problem to solve parametric function into the problem to solve the linear combination of at-q. The choice of the base function depends on the adequate expression to nonstationary response data. In this paper, we choose the polynomial function as base function as follows:
(49)gi(t)=(tN)i,i=0,1,…,m,
where N is the number of samples in the input (or output) data and i is the order of base function to be chosen. When introducing the base function to construct a nonstationary time-series model, the reconstructed ARMA model form is similar to that of the original time-series model. The major difference is that the MA parameters in the nonstationary base-function ARMA model are treated as functions of time and may be represented by a linear combination of suitable base functions. As long as by properly expanding the data matrix in the nonlinear and three-stage least-square method, respectively, the AR parameters and the coefficients of MA part corresponding to base function can be obtained. Then we construct nonstationary base-function ARMA model to describe the behavior of the nonstationary response data and further determine the modal parameters of a structural system. Note that although the structural modal parameters can be identified from the AR parameters only, the estimation of AR parameters and MA parametric function will affect each other in practice; estimation for AR parameters and MA parametric function cannot therefore be separately and solely performed.


## 5. Numerical Simulation

To demonstrate the feasibility and effectiveness of the proposed method, we consider a linear 3-DOF chain model with viscous damping. A schematic representation of this model is shown in Figure1. Assume that the three blocks of mass are 3, 1.5, and 1 Kg, respectively, as well as all spring constants being 4, 1, 1, and 2 N/m, respectively. The mass matrix M, stiffness matrix K, and the proportional damping matrix D of the system are given as follows:
(50)M=[30001.50001][kg],K=100[5-10-12-10-13][N/m],D=0.36M+0.002K[N·sec/m].Figure 1
Schematic plot of the 3-DOF chain model.Note that the system has proportional damping, since the damping matrixD can be expressed as a linear combination of M and K. Consider that the ambient vibration input can be modeled as nonstationary white noise as represented by the product model. The stationary white noise is generated using the spectrum approximation method [11] as a zero-mean band-pass noise, whose standard deviation is 0.02 N2·sec/rad with a frequency range from 0 to 50 Hz. The sampling interval is chosen as Δt=0.01 sec, and the sampling period is T=655.36 sec. The stationary white noise simulated is then multiplied by an amplitude-modulating function, as shown in Figure 2, to obtain the nonstationary white noise, which serves as the excitation input acting on the 3rd mass point of the system. The time signal of a simulated sample of the nonstationary white noise and the power spectrum of the corresponding stationary part are shown in Figures 3 and 4, respectively. The displacement responses of the system obtained through Newmark’s method [12] are shown in Figure 5. By examining the Fourier spectra associated with each of the response channel, also as shown in Figure 5, we chose the response of the 3rd channel, X3(t), which contains rich overall frequency information, as the reference response to construct a nonstationary ARMA model of the system.Figure 2
A typical plot of the amplitude-modulating function.Figure 3
Sample function of nonstationary white noise in time domain.Figure 4
Power spectrum associated with the stationary part of the nonstationary white noise.Figure 5
Typical displacement responses and the corresponding Fourier spectra of the 3-DOF chain system subject to nonstationary white input.According to the theory presented in the previous sections, by evaluating temporal root-mean-square function in conjunction with curve-fitting technique, as shown in Figure6, nonstationary problem may reduce to a stationary problem if we can extract the amplitude-modulating function from the nonstationary data, as shown in Figure 7. Therefore, we can follow the same procedures as those for stationary problems, and the original time-series algorithm can thus be performed. We then consider ARMA (6,5) model to perform the modal-parameter identification of the 3-DOF structural system. The initial parameters of the ARMA model can then be estimated through the three-stage least-square method, as shown in Table 1, and the relatively precise ARMA parameters can further be determined by using the nonlinear least-square methods, as listed in Table 2. We construct a nonstationary time-series model by introducing the aforementioned temporal root-mean-square function and then extracting the AR parameters from the above constructed nonstationary time-series model. The modal parameters of a structural system can then be obtained from (15), as summarized in Table 3, which show that the errors in natural frequencies are less than 1% and the maximum error in damping ratios is about 12%.Table 1
The results of initial estimation of ARMA model for quasistationary response data through the three-stage least-square method.
AR parameter
MA parameter
φ
1
−5.958
θ
1
0.404
φ
2
14.851
θ
2
0.940
φ
3
−19.820
θ
3
0.365
φ
4
14.937
θ
4
0.273
φ
5
−6.027
θ
5
0.137
φ
6
1.017Table 2
The results of relatively precise estimation of ARMA model for quasistationary response data through the nonlinear least-square method.
AR parameter
MA parameter
φ
1
−5.916
θ
1
0.312
φ
2
14.641
θ
2
0.771
φ
3
−19.403
θ
3
0.127
φ
4
14.521
θ
4
0.103
φ
5
−5.819
θ
5
0.001
φ
6
0.976Table 3
Results of natural frequency and damping ratio identification of the 3-DOF system subjected to nonstationary white noise through the nonstationary ARMA algorithm with curve-fitting technique.
Mode
Natural frequency (rad/sec)
Damping ratio (%)
Exact
ARMA
Error (%)
Exact
ARMA
Error (%)
1
8.84
8.84
0.01
2.95
3.30
11.86
2
13.66
13.63
0.22
2.74
2.91
6.20
3
18.31
18.29
0.11
2.90
2.97
2.41Figure 6
Schematic plot of the temporal root-mean-square function through curve-fitting technique.Figure 7
Quasistationary displacement response transformed from nonstationary one.In addition, we also consider ARMA(6,5) model in conjunction with a linear combination of suitable base functions to perform the modal identification of the previous 3-DOF chain-model system. The 5-order polynomial functions of time are chosen instead of the unknown MA parameter function in the original ARMA algorithm. The initial AR parameters and the coefficients of known base functions can then be estimated through the three-stage least-square method, as shown in Table 4, and the aforementioned AR parameters and coefficients of base functions can further be relatively precisely determined by using the nonlinear least-square methods, as listed in Table 5. We extract the AR parameters from the above constructed nonstationary time-series model, and the modal parameters of a structural system can then be obtained from (15), as summarized in Table 6, which show that the errors in natural frequencies are less than 1% and the maximum error in damping ratios is about 2%.Table 4
The results of initial estimation for the AR parameters and the coefficients of polynomial base function from the nonstationary response data through the three-stage least-square method (*θi(t)=∑j=0mkji·gj(t), m=, i=1,2,…,5).
AR parameter
The coefficient of polynomial base function
i
=
1
i
=
2
i
=
3
i
=
4
i
=
5
φ
1
−5.946
k
0
i
0.653
1.122
0.712
0.355
0.173
φ
2
14.790
k
1
i
2.705
−0.751
−0.847
0.492
−0.271
φ
3
−19.697
k
2
i
−10.108
3.982
3.532
−1.762
0.981
φ
4
14.815
k
3
i
10.482
−6.499
−4.345
2.678
−2.100
φ
5
−5.966
k
4
i
−2.871
3.336
1.576
−1.392
1.373
φ
6
1.005Table 5
The results of relatively precise estimation for the AR parameters and the coefficients of polynomial base function from the nonstationary response data through the nonlinear least-square method (*θi(t)=∑j=0mkji·gj(t), m=, i=1,2,…,5).
AR parameter
The coefficient of polynomial base function
i
=
1
i
=
2
i
=
3
i
=
4
i
=
5
φ
1
−5.917
k
0
i
0.591
0.859
0.259
0.049
−0.030
φ
2
14.648
k
1
i
2.493
1.075
1.155
1.086
0.714
φ
3
−19.417
k
2
i
−9.857
−4.254
−4.470
−3.609
−2.805
φ
4
14.535
k
3
i
10.635
4.915
4.673
4.482
4.141
φ
5
−5.826
k
4
i
−3.112
−1.687
−1.333
−1.905
−2.071
φ
6
0.977Table 6
Results of natural frequency and damping ratio identification of the 3-DOF system subjected to nonstationary white noise through the nonstationary ARMA algorithm with polynomial base function.
Mode
Natural frequency (rad/sec)
Damping ratio (%)
Exact
ARMA
Error (%)
Exact
ARMA
Error (%)
1
8.84
8.82
0.23
2.95
3.00
1.69
2
13.66
13.65
0.07
2.74
2.69
1.82
3
18.31
18.29
0.14
2.90
2.92
0.69From Tables3 and 6, it is indicated that the errors in results of natural frequency and damping ratio identification through the nonstationary ARMA algorithm with polynomial base function (as summarized in Table 3) are relatively smaller than those with curve-fitting technique (as summarized in Table 6), especially for the identified damping ratios. It should be mentioned that the errors of identified damping ratios are somewhat larger than those of identified natural frequencies. This is due to the fact that the system response generally has lower sensitivity to damping ratios than to the natural frequencies. In general, the amplitude curve-fitting technique has relatively high computational efficiency when compared with base-function method; however, it is only available under the assumption of product model for nonstationary input process. The base-function ARMA algorithm is well-performed without the assumption of product model for nonstationary input process and treatment of curve-fitting technique for extracting envelope function. In addition, the frame of mathematical model of base-function ARMA algorithm is relatively adjustable and is therefore applicable to general nonstationary ambient vibration data.To further verify the effectiveness of the proposed method for the higher DOF structural systems, as well as for the substitution of the exact input with the data disturbed by some noise corresponding to the experience, we consider a 20-DOF chain model with nonproportional damping excited by a vibration practically recorded at Sun Moon Lake on the 21st of September, 1999, when the Chi-Chi Earthquake of a moment magnitude 7.6 occurred in central Taiwan. The sampling interval and period of this seismic record areΔt=0.005 sec and T=59.995 sec, respectively. A sample of the seismic record, which serves as the excitation input acting on the 20th mass of the model, is shown in Figure 8. Assume that each mass is 1 kg and all spring constants are 600 N/m. The damping matrix of the system is assumed to be
(51)C=0.05M+0.001K+0.02[1…1⋮⋱⋮1…1]20×20[N·s/m].
In this example, we use the practical nonstationary input acting on the twentieth mass of the chain model, and then the corresponding displacement responses obtained by Newmark’s method [12] are used for modal identification. The results of identification are summarized in Table 7, which shows that the modal identification in this case is satisfactory. It is also noted that the higher modes generally are not identified so accurately as the lower ones, because the contribution of higher modes to the system response is somewhat less than that of the lower modes.Table 7
Results of natural frequency and damping ratio identification of the 20-DOF chain system subjected to a recorded sample of the Chi-Chi Earthquake through the nonstationary base-function ARMA algorithm in conjunction with the proposed modified least-square method.
Mode
Natural frequency (rad/sec)
Damping ratio (%)
Exact
ARMA
Error (%)
Exact
ARMA
Error (%)
1
3.66
3.63
0.93
5.49
5.49
0.07
2
7.30
7.30
0.03
0.71
0.68
4.79
3
10.90
10.89
0.11
0.94
1.00
6.17
4
14.44
14.41
0.19
0.89
0.82
8.31
5
17.90
17.85
0.27
1.06
1.04
2.08
6
21.26
21.18
0.37
1.17
1.10
5.81
7
24.49
24.37
0.49
1.33
1.30
2.11
8
27.60
27.42
0.67
1.45
1.39
3.86
9
30.54
30.32
0.72
1.59
1.78
11.70
10
33.32
32.97
1.04
1.71
1.68
1.64
11
35.91
35.43
1.33
1.83
1.81
0.87
12
38.30
37.88
1.11
1.94
1.87
3.40
13
40.48
39.90
1.43
2.03
2.14
5.42
14
42.43
41.85
1.38
2.12
1.77
16.32
15
44.14
43.37
1.74
2.20
2.43
10.64
16
45.60
44.78
1.79
2.26
2.01
10.97
17
46.81
45.46
2.88
2.26
2.53
12.04
18
47.76
46.48
2.69
2.20
1.90
13.55
19
48.44
47.14
2.68
2.51
2.67
6.29
20
48.85
48.03
1.67
2.41
2.64
9.54Figure 8
A recorded sample of the Chi-Chi Earthquake.Recall that the conventional ARMA algorithm is effective only for stationary process [6, 7]. Therefore, in order to convert the original nonstationary responses into quasistationary data, we need additional treatments, such as employing the curve-fitting technique, to extract the amplitude-modulating function from the original nonstationary data in the form of a product model [4]. In contrast, in the present method, by introducing the base functions to construct a nonstationary time-series model, the MA parameters in the nonstationary time series are treated as functions of time and may be represented by a linear combination of suitable base functions. Without using input data or any additional treatment of converting the original nonstationary data into the form of quasistationary vibration, we directly perform the modal identification from nonstationary response data through the nonstationary base-function ARMA algorithm in conjunction with the proposed modified least-square method. Therefore, we avoid a distortion in the modal parameters of identification induced by the error involved in the approximate quasistationary response obtained through curve-fitting technique.In this paper, we developed the nonstationary base-function ARMA algorithm in conjunction with the proposed modified least-square method for modal identification from nonstationary ambient response data only. We also demonstrated the validity of these methods through numerical simulations without using the practical response data. From the vibration behavior of realistic ambient excitation and the theory of the proposed method, we know that the assumptions and approximations are consistent with the time-varying nature of ambient excitation in practice. Thus, the proposed methods are generally applicable in identifying the modal parameters of a structure from the identification results obtained through numerical simulations.


## 6. Conclusions

An extended time-series algorithm is presented in this paper for modal identification from nonstationary ambient vibration data only. Nonstationary ARMA model with time-varying parameters is considered because of its capability in resolving general nonstationary problems. Without changing the AR parameters, we construct a nonstationary time-series model to effectively describe a wide variety of nonstationary behavior and then extract the modal parameters of a structural system. Furthermore, by introducing the base functions to construct a nonstationary time-series model, the MA parameters in the nonstationary time series are treated as functions of time and may be represented by a linear combination of suitable base functions. We can then solve the identification problem of time-varying parameters and extend the time-series algorithm to identify the modal parameters of a structural system subjected to nonstationary white noise. However, identification of the mode shapes is still a challenging problem to be resolved in the proposed method from only one single sample of response signal. Through numerical simulation, applicability and effectiveness of the proposed method of modal parameter identification from nonstationary ambient vibration data are demonstrated.



---

*Source: 391815-2014-07-17.xml*