Modal Identification Using OMA Techniques : Nonlinearity Effect

This paper is focused on an assessment of the state of the art of operational modal analysis (OMA) methodologies in estimating modal parameters from output responses of nonlinear structures. By means of the Volterra series, the nonlinear structure excited by random excitation is modeled as best linear approximation plus a term representing nonlinear distortions. As the nonlinear distortions are of stochastic nature and thus indistinguishable from the measurement noise, a protocol based on the use of the random phase multisine is proposed to reveal the accuracy and robustness of the linear OMA technique in the presence of the system nonlinearity. Several frequencyand time-domain based OMA techniques are examined for the modal identification of simulated and real nonlinearmechanical systems.Theoretical analyses are also provided to understand how the systemnonlinearity degrades the performance of the OMA algorithms.


Introduction
Operational modal analysis (OMA) aims at identifying the modal properties of a structure based on response data of the structure excited by ambient sources.Unlike experimental modal analysis (EMA) [1], the OMA is performed in operational conditions of the vibrating structure, where the excitation is inaccessible to be measured or hard to be applied.The research activity around the theoretical basis of the OMA has been largely increased, and powerful OMA techniques have been developed in both the time-and frequency-domains, as collected in recent tutorial work [2,3].The time-domain methodologies estimate the modal model based on a statespace representation of the system obtained from the timedomain data [4], while the frequency-domain approaches identify the modal parameters from power spectral density functions of output responses [5,6].System identification methods for the OMA are extensively reviewed in [7,8], different OMA and EMA methods for the estimation of linear time-invariant (LTI) dynamical models are compared in an extensive Monte Carlo simulation study [8].The state of the art of OMA methodologies is also assessed for estimating modal parameters of a helicopter structure in a laboratory environment, where the helicopter structure behaves linearly by controlling the excitation level [9].
The OMA techniques have been utilized to estimate the modal properties of important structures for the purpose of health monitoring and maintenance, such as wind turbines [10], bridges [11], and historical aqueduct [12].Generally, with the real-life vibration data, the detection of structural damage by estimating the modal properties is complicated by the impact of changing environmental conditions and the nonwhiteness property of the ambient excitation.Both of them can induce significant changes of the monitored alert feature.The former problem is treated using factor analysis, and damage is detected using statistical process control [13].The latter problem can be ideally solved through the concept of the transmissibility, as the input excitation is canceled in computing transmissibility functions between output responses.The transmissibility functions are fully exploited for the OMA, as reported in [14,15].In addition to the above two factors, the output-only modal estimate is also perturbed by the system nonlinearity.Indeed, the OMA techniques are all derived based on LTI dynamical models, while all the actual structures are nonlinear to some extent.Moreover, the structures to which the OMA approaches are applied are assemblies, often utilize composite materials, and even have flexible components.Diverse types of nonlinearities can be present, such as geometric nonlinearity due to the large deformation experienced by a flexible structure, material nonlinearity owing to a nonlinear stressstrain constitutive law, contact nonlinearity resulting from boundary conditions.Very often, the real-life structures to be monitored are exposed to nonstationary and even severe ambient excitation, for instance, caused by earthquake or hurricane.Within such context, it is therefore of paramount importance to study the influence of the system nonlinearity on the output-only modal estimation for the purpose of the (fully automatic) structural damage detection.To the authors' knowledge, this fundamental issue still remains to be addressed.The way of dealing with the system nonlinearity can be split into two groups.The first group considers explicitly a nonlinear model such as nonlinear ARX (autoregressive with exogenous inputs) for parameter estimation [16,17].The second group still considers linear models while bounding the estimation variability due to the nonlinear distortions.The latter fits in the scope of this paper as the goal of the present work is to assess the performance of the LTI model based OMA technique with respect to the system nonlinearity.
The system nonlinearity strongly depends on the class of excitation signal.The ambient excitation is usually not of Gaussian type.However, approximate Gaussian distribution can occur in many situations, as explained by the Central Limit Theorem.As reported in literature, only the first-and second-order statistical properties of the ambient excitation are exploited in the OMA algorithms.So a normally distributed ambient excitation is considered herein.The present work is confined to nonlinear Wiener systems that can be approximated in the mean square sense by a Volterra system for Gaussian excitation.The major property of a Wiener system is that the steady state response to a periodic input is periodic with the same period as the input.The nonlinear structure of the Wiener system can be modeled as a related linear system plus a part representing the nonlinear distortions under random excitation [18,19].The latter is proved to be asymptotically normally distributed, mixing of order infinity over frequency under Gaussian excitation [18].Then, considering them in the weighting matrix of the cost function, applying LTI system identification method will return best (in the least square sense) linear approximation (BLA) of the nonlinear system.
Disturbed by the measurement noise in observation data, the stochastic nonlinear distortions are not recognizable, usually naively treated as independent noise.Therefore, the issues of separating them from noisy data and of generating their stochastic realizations should first be addressed in order to quantify the performance of the OMA algorithm.To this end, the random phase multisine (RPM) is advocated as an excitation signal for the OMA.The RPM is normally distributed by increasing the number of harmonic lines and advantageous over the Gaussian noise due to its flexible design and its periodicity.The amplitudes of harmonic components in the RPM are set constant to meet the whiteness assumption of the excitation, as required by most OMA algorithms.The random phases of the RPM then become the sole factors to determine the nonlinear distortions of the nonlinear system.Additionally, the periodic property of the RPM leaves Fourier transformation of input/output data free of the leakage error.
Despite the existence of a large number of alternative OMA techniques developed during the last decades, they are, however, just based on few basic principles.Therefore, the present work is mainly focused on three representative methods with completely distinct theoretical backgrounds: frequency-domain decomposition (FDD), stochastic subspace identification (SSI), and transmissibility based operational modal analysis (TOMA).The main contribution of the present work is to provide a protocol to assess the modal estimation of nonlinear structures by the OMA techniques, which is based on the use of the RPM.The protocol comprises two parts: one part is to evaluate the accuracy of the LTI model based OMA algorithm in comparison with the well established input-output approach; the other is to examine its robustness to the nonlinear distortions.Further, the protocol is carried out for a series of excitation levels in order to reveal extensively the performance of the OMA techniques in the presence of the system nonlinearity.Also, theoretical analyses are provided to understand how the system nonlinearity influences the accuracy of the OMA algorithms.
The rest of the paper is structured as follows.Section 2 is devoted to building a nonlinear framework for system identification.Section 3 recapitulates the theoretical background of the considered OMA techniques.Section 4 describes the protocol dedicated to investigating OMA techniques in the presence of the structural nonlinearity.Section 5 illustrates the output-only identification of simulated and real nonlinear structures.Concluding remarks are given in Section 6.

Nonlinear System Modeling
2.1.Random Phase Multisine.The RPM plays a crucial role in designing a protocol to assess the accuracy of estimating nonlinear structures by OMA techniques.With a uniform amplitude  0 , it takes the following form: where   is the clock frequency of the arbitrary waveform generator,  is the number of samples in one signal period, and the set of the phases {  } is a realization of an independent distributed random process on [0, 2) such that E[exp(  )] = 0 and  − =   .The amplitude of the RPM is used to set the excitation level, and its random phases determine the nonlinear distortions of the system.

Volterra Series of a Nonlinear System Output.
A description of a nonlinear system by means of the Volterra series Nonlinear dynamic system Best linear approximation Modeling of a nonlinear system under a normally distributed excitation. 0 () and  0 () are the input-output noise-free spectra,  BLA () is the best linear approximation of the nonlinear system,  Lin () is the linear part of  0 (), and   () represents the nonlinear distortions in  0 ().
is formally introduced and the stochastic property of the system nonlinearity is presented.The case of single output is considered for the ease of illustration.By the Volterra series, the output of a nonlinear system is decomposed into where  () () is the contribution of degree where () denotes the excitation (driven by the RPM) and ℎ () ( 1 , . . .,   ) ( > 1) is the generalization of the impulse function ℎ (1) () and is referred to as Volterra kernel, whose symmetrized frequency-domain representation is written as Applying the discrete Fourier transform to (2) along with (4), it follows that at the th frequency line () = ∑ +∞ =1  () (), with and

Best Linear Approximation and Nonlinear Distortions.
The frequency response function of a nonlinear system is by definition computed as where  0 () is the underling linear system and ∑ +∞ =2  () () represents the system nonlinearity.
Under the excitation of the RPM, ∑ +∞ =2  () () is further decomposed into a systematic part   () which depends only on the amplitude of the RPM and a stochastic part denoted by   () which depends on both the amplitude and the phase of the RPM.  () behaves as uncorrelated noise (over frequency line), and further   () is an asymptotically zeromean circular complex normal variable with respect to the random realization of the phases of the RPM [18]. 0 () along with   () constitutes the best linear approximation (BLA) of the nonlinear system in the least square sense, It is stressed that, through the term   (),  BLA () depends on the loading condition (level or location of the excitation) applied to the distributed nonlinear structure.By (7), the output of a nonlinear system is accordingly split into two parts: the first part that is related to the input and the second part   () that is uncorrelated with the input over the random phases of the RPM, as illustrated in Figure 1.The output nonlinear distortions   () are related to   () as and have similar stochastic properties as   () (see [18] for proof details). 0 () is the force applied by the exciter, which is controlled with the RPM defined in (1).

Theoretical Background of Operational Modal Analysis Techniques
In this section, the theoretical background of the OMA techniques considered in this paper is briefly recapitulated, Note that they have all been developed with the purpose of extracting the underlying system  0 () in (6).More details can be found in the cited references.

Frequency-Domain Decomposition (FDD).
The frequency response matrix G() is expressed as a sum of the contributions of system modes in a frequency band of interest, where R  =      with   the column vector of the th mode shape and   the column vector of the modal participation factor and   = (−  +√1 −  2  ) 0, with   the damping ratio and  0, the angular frequency (rad/s).
Theoretically, assuming that the ambient excitation has a constant spectrum S  , the noiseless output power spectral density matrix is written as Substituting G() by its modal decomposition in (9) and using the Heaviside partial fraction theorem for polynomial expansions, S YY () can be written as where In the vicinity of the th eigenfrequency,  ∈ sub( 0, ), it approximatively holds by exploiting the lightly damped property that Then for  ∈ sub( 0, ) it is derived that where Re() is the real part of the complex number .
Assuming the orthogonality of the mode shapes, ( 14) can be interpreted as a singular value decomposition (SVD) of S YY () where only one single component is dominant while  ∈ sub( 0, ).This suggests a simple procedure to estimate the modal parameters: decomposing the spectral response into a set of single degree of freedom systems using the SVD, each corresponding to an individual mode.The poles are estimated by the enhanced FDD [20].The unscaled mode shape is obtained by picking out the corresponding eigenvector of the estimated power spectral density matrix.
With the real-life data, the estimate of the output power spectral density matrix can be (heavily) biased due to the presence of the nonlinear distortions which are mutually correlated over space in the frequency-domain.The use of the SVD can alleviate this situation by discarding lower singular values which correspond to less significant components in the noisy data.

Stochastic Subspace Identification (SSI).
The SSI method estimates the modal parameters based on the stochastic discrete-time state-space model of a mechanical structure, where the subscript  denotes the time instant, x  is a vector of the system state, y  is an output vector, A is the discrete state matrix, C is the selection matrix, w  is a vector of noise due to the (white) random excitation, and k  is a vector of noise representing the sum of the random excitation and the measurement noise.Two versions of the SSI are typically used [4,21]: data driven SSI and covariance driven SSI (SSI-COV).The latter is taken for algorithmic explanation in what follows.The SSI-COV algorithm starts with the covariance matrix of the structural response with   reference outputs y ref  of them, where the expectation operation (E) is performed with the purpose of correlating out the unwanted dynamics to obtain the system description.However, the nonlinear distortions in the output noise term k are mixing of infinite order over time (self-correlated), whose components are also mutually correlated over space.This leads to a biased covariance matrix Λ ref  , and eventually decreases the estimation accuracy.
A block Toeplitz matrix is constructed by assembling the output covariance matrices Using where O   is called the extended observability matrix and C ref   is the reference-based stochastic controllability matrix.Theoretically, from the observability matrix O   , the state matrix A and the output selection matrix C can be obtained.The natural frequencies, damping ratios, and unscaled mode shapes are finally derived from the estimated matrices A and C.

Transmissibility Based Operational Modal Analysis (TOMA).
The case of the scalar transmissibility function is considered for the algorithmic illustration.The scalar transmissibility function is defined in the Laplace domain as the ratio of the output spectrum at the th degree of freedom (dof) and the one at the th dof under the excitation of an unknown force at the th dof, The kernel idea of the TOMA approach is that the scalar transmissibility functions of a linear structure, estimated with the response data from different loading conditions, cross each other at the poles of the system.Using ( 9), the limit value of the transmissibility function when  →   lim which is independent of the input.Combining the transmissibility functions of two distinct loading locations () and (), it follows that lim Different ways are reported in the literature to extract the modal parameters based on (21); see, for instance, [14,15].The approach employed here is based on the parametrically estimated scalar transmissibility functions.The transmissibility function admits a rational form, as indicated by (19).
Choosing the th output as reference, a common denominator rational model is used to parameterize the transmissibility functions, whose parameters are estimated using the sample maximum likelihood method [18].Then, using the transmissibility functions estimated from the loading conditions () and (), respectively, the function used for modal identification is defined as where D denotes the set of the dofs of the outputs except the th one.  can be easily derived from θ by means of [1,1] (t) [2,1] (t) [N  ,1] (t) [1,2] (t) [2,2] (t) [N  ,2] (t) [1,N  ] (t) [2,N  ] (t) [N  ,N  ] (t) symbolic computation.The poles are then obtained from   and the unscaled mode shapes are estimated by evaluating the identified transmissibility functions at the estimated poles.Note that the TOMA method is superior to other outputonly identification techniques in dealing with colored ambient excitation; however, the need of combining different loading conditions makes it vulnerable to the system nonlinearity.

Assessment Protocol
The output-only modal identification of nonlinear structures is conducted for a series of excitation levels.At each excitation level, a protocol is applied to assess the algorithm performance with respect to the nonlinearity, which is established in the following sections.

Uncertainty Bound Based on Multiple Experiments.
With the RPM defined in Section 2.1 as an excitation signal, the measurement noise accounts for the difference of the observed data over the periods while the nonlinear distortions are deterministic once the phases of the RPM are fixed in an experiment.Therefore, the measurement noise is kicked out by averaging vibration data over period and the denoised output data are mainly corrupted by (a realization of) the nonlinear distortions.The efficient way for examining the robustness of the algorithm is to provide an uncertainty bound induced by the system nonlinearity in the form of formula.However, unlike the independent disturbing noise, the nonlinear distortions are dependent on input signals.This will create higher order moments in the variance calculations for modal estimates, which cannot be captured in a linear system identification framework.As a result, a Monte Carlo analysis based on multiple experiments with different random realizations of the RPM is needed.
The strategy of multiple experiments aims at generating a set of realizations of the nonlinear distortions, over which the robustness of the OMA algorithm is completely assessed for a specified excitation level, as depicted in Figure 2.Each experiment is conducted with an independent phase realization of the RPM, and a number of consecutive periods of the steady state response are measured.The dataset, collecting output data from all the experiments, is used to investigate the uncertainties of the modal estimates in light of the following procedure: (ii) Repeat steps (a)-(c) for all the   experiments, a set of modal estimates being delivered, over which any order statistics can be empirically computed.

Identification of Best Linear Approximation.
The OMA approaches fail to extract the underlying linear system due to the presence of the bias term   in (7) induced by the nonlinear behavior of the system.Therefore, their accuracy is herein evaluated in comparison with the input-output EMA approach.In the presence of the nonlinear distortions, the EMA approach extracts the modal parameters based on the BLAs of the nonlinear structure.The identification of the BLA can be done, for instance, using the sample maximum likelihood method which is optimal in the sense that the variance of the estimate is close to the Cramér-Rao lower bound [18].

Applications
The OMA techniques are applied to identify nonlinear systems excited by a series of excitation levels.The excitation level is expressed in the form of signal-to-noise ratio defined as follows: where ĜBLA and  ĜBLA are the identified BLA and quantified standard deviation, respectively.

Simulated Case
5.1.1.Nonlinear System.The simulated example is a 4DOF mass-spring-damper model with local nonlinearity and the spring stiffness and damping at the ends are state-dependent, as described in Figure 3(a).

Parameter Setting.
The sampling frequency is 64 Hz,   = 4,   = 100.In order to illustrate the algorithmic behavior very clearly, we take the extreme point of view that no measurement noise is present.The simulated data are only corrupted by the nonlinear distortions after discarding the periods in response corrupted by the transient effect.
Three levels of excitation are considered, the corresponding SNR dB , 89.8, 46.7, and 37.6, are classified as weak, medium, and strong, respectively.The evolution of the BLAs and of the associated nonlinear distortions of the 3rd mass is shown in Figure 3(b) when the input is located at the 4th mass.The modal assurance criterion (MAC) index is set equal to 0.85 for the FDD method to choose the frequency vicinity around a pole.Consider   = 15 in (17) for the SSI-COV algorithm.The number of measured outputs (  ) directly influences the performance of the FDD, SSI, and TOMA through (10), (16), and (22), respectively.The more reliable results would be obtained with more output data.However,   is limited by the data acquisition system.Here,   is set as 4 for the output-only and input-output modal identification.Two loading conditions are created for the TOMA approach by shifting the excitation location from the 1st mass to the 4th one.

Results and Discussion
. The multiple simulated data are generated (as shown in Figure 1), based on which the previously described OMA and EMA approaches are applied to obtain the modal estimates.The estimated eigenfrequencies and damping ratios are normalized with respect to those identified using the input-output data at the lowest level of excitation.The mode shape estimated by the OMA is normalized with respect to the one obtained by the EMA through the use of the MAC.These normalized modal values can provide a clear view on the accuracy of the OMA algorithms, whose dispersion vividly demonstrates the algorithmic robustness with respect to the nonlinear distortions.
As expected, Figures 4-6 depict that in the case of the lowest excitation level the considered OMA techniques all identify well the system which behaves linearly.The modal estimates are shown to become more dispersed and even biased but still situated around those of the BLAs when severe nonlinear distortions are stimulated by increasing the excitation level.The performance of the FDD method is remarkable except for the mode shape estimates as they are obtained by picking the discrete values determined by the frequency resolution of the spectrum.As commented in Section 3.3, the eigenfrequency estimates by the TOMA approach generally deviate more from the EMA estimates than the other methods.However, it is worth pointing out that the hardening behavior of the nonlinear system is most captured by all the OMA techniques, as reflected by the estimated eigenfrequencies.

Real Case
5.2.1.Nonlinear System.The nonlinear system under investigation is a mechanical structure consisting of a beam, frame, and support attachment, as displayed in Figure 7(a).The main part of the system is formed by fastening two beams together, whose left side is clamped in a rigid support and right side is adhered to a steel frame put freely on the ground.The use of various assemblies and joins in this setup creates diverse underlying sources of nonlinearity, such as beam clamping which introduces nonlinear stiffness of (possibly high-order) polynomial form, frictional slips at loosened interfaces which introduce additional flexibility and hysteretic damping to the overall structural dynamics.Also, it commonly occurs that nonlinearity is unintentionally introduced in the measurement chain, for example, preloading the beam because of insufficient human checks, which generally creates a nonlinear stiffness of cubic form.All kinds of nonlinearity are treated in a unified way based on the established protocol, whose effects on structural dynamics are determined by identifying the BLA and bounding the stochastic nonlinear distortions.

Parameter Setting.
The sampling frequency is 8192 Hz,   = 256,   = 40.Following the same lines of the simulation, three levels of excitation are set by the SNR dB : 60.3 (weak), 48.7 (medium), and 33.5 (strong).Applying the shaker at the position (1), Figure 7(b) shows the apparent shift-down of eigenfrequencies of the considered system by raising the excitation level, which is in large part due to the use of various joints in the setup.In addition to the softening effect, two close but distinct poles are present around 300 Hz for low levels of excitation.These poles are transformed into a single pole when a larger amount of power is injected into the system.The presence of the system nonlinearity is classically indicated in terms of the coherence function in Figure 8(a), and Figure 8(b) displays several realizations of the nonlinear distortions.The MAC index is set equal to 0.9 for the FDD method.Consider   = 25 for the SSI-COV algorithm.The shaker is applied at locations (1) and (2) on the beam, respectively, creating distinct loading conditions for the TOMA.

Results and Discussion
. When an electrodynamic shaker is used to excite the structure in modal testing, it is not trivial to investigate the considered OMA techniques in a fair way.In fact, the force applied on the structure is the reaction force between the shaker and the beam when the armature mass and spider stiffness of the shaker are not negligible, whose magnitude and phase depend upon the characteristics of the structure and of the exciter.Thus the shaker driving force does not meet the fundamental assumption of white noise.Consequently, the FDD and SSI-COV algorithms actually identify the whole system including the nonlinear structure and the shaker part.Superior to both of OMA techniques, the TOMA approach can identify only the nonlinear structure of interest as the shaker part  (1) is simultaneously present in different outputs and further canceled by computing transmissibility functions.The output data are processed in line with the proposed protocol (see Section 4.1).As shown in Figures 9 and 10, the eigenfrequencies estimated by the SSI-COV and FDD algorithms agree well with those extracted from the BLAs for the three excitation levels.It is also found by comparing with EMA estimates that the eigenfrequency estimates are more accurate than those of damping ratios for the considered OMA techniques.As analyzed in Sections 3.1 and 3.2, the estimated S YY () used for the FDD and Λ ref  for the COV-SSI are biased in the presence of the correlated nonlinear distortions over space and time.The nonlinear perturbations in these covariance matrices are propagated to the identification of damping ratios, especially for the SSI-COV without applying any filter, as verified in Figure 10.
Shifting the shaker between the two locations induces a remarkable evolution of the BLA, as clearly demonstrated in Figure 11.The modal parameters are estimated by combining the output dataset from these two loading configurations, as seen in Figure 11.The eigenfrequency estimates still approach quite well those of the BLAs with a relative error up to only 2.5%.They also can reflect the shift-down property of the system, while the damping estimates are less meaningful.In fact, (21) does not hold any more in the presence of the system nonlinearity such that lightly damping properties are hardly extracted with reliable accuracy.
Although the approaches in time-and frequencydomains behave similarly, it is found from the simulated and real case studies that the FDD approach is characterized by the highest accuracy for the eigenfrequency estimates.

Conclusions
The OMA techniques have been applied to identify simulated and real nonlinear structures; a full view on the algorithmic performance is delivered by using a protocol that establishes the estimation accuracy and robustness with respect to the system nonlinearity.A theoretical motivation for the proposed protocol is also provided.The following conclusions and remarks can be drawn from our study.
(i) The output-only identification techniques are able to extract most linear dynamics of a nonlinear structure, but the estimates obtained by them are more biased and dispersed in the presence of increasing nonlinearity induced by the excitation change.The estimated eigenfrequency is a robust indicator of the system state induced by the nonlinearity.
(i) Consider ∀ = 1, . . .,   .(a) Feed the shaker with the RPM of the th phase realization and acquire the noisy output signals {y [,] ()}   =1 .(b) Average the output data over the periods, leaving them mainly corrupted by the nonlinear distortions, noise in ŷ[] () vanishes at the rate of 1/√  .(c)Apply each OMA algorithm to ŷ[] () to extract eigenfrequencies, damping ratios, and unscaled mode shapes.

Figure 4 :
Figure 4: The rows show four modes and the columns represent the probability distributions of eigenfrequencies obtained by the SSI-COV, FDD, and TOMA for three excitation levels (light gray: weak, dark gray: medium, and black: strong).Dashed lines: EMA estimates with the force applied at  4 .

Figure 5 :
Figure 5: The rows show four modes and the columns represent the probability distributions of damping ratios obtained by the SSI-COV, FDD, and TOMA for three excitation levels (light gray: weak, dark gray: medium, and black: strong).Dashed lines: EMA estimates with the force applied at  4 .

Figure 6 :
Figure 6: The rows show four modes and the columns represent the probability distributions of MAC indices obtained by the SSI-COV, FDD, and TOMA for three excitation levels (light gray: weak, dark gray: medium, and black: strong).

Figure 7 :
Figure 7: (a) Experiment setup of the real structure, (b) best linear approximations under three excitation levels (black: strong, dark gray: medium, and light gray: weak).

Figure 9 :
Figure9: The columns show four modes and the rows represent the probability distributions of eigenfrequencies obtained by the SSI-COV and FDD with the shaker at location (1) for three excitation levels (light gray: weak, dark gray: medium, and black: strong).Dashed lines: the EMA estimates.