Wavelet spectrum analysis approach to model validation of dynamic systems

Mechanical Systems and Signal Processing 25 (2011) 575–590

Contents lists available at ScienceDirect

Mechanical Systems and Signal Processing journal homepage: www.elsevier.com/locate/jnlabr/ymssp

Wavelet spectrum analysis approach to model validation of dynamic systems Xiaomo Jiang a,, Sankaran Mahadevan b a b

4200 Wildwood Parkway, 42-3-12C-6, GE Energy, Atlanta, GA 30339, USA Department of Civil and Environmental Engineering, Box 1831-B, Vanderbilt University, Nashville, TN 37235, USA

a r t i c l e in fo

abstract

Article history: Received 19 October 2009 Received in revised form 8 April 2010 Accepted 17 May 2010 Available online 1 June 2010

Feature-based validation techniques for dynamic system models could be unreliable for nonlinear, stochastic, and transient dynamic behavior, where the time series is usually non-stationary. This paper presents a wavelet spectral analysis approach to validate a computational model for a dynamic system. Continuous wavelet transform is performed on the time series data for both model prediction and experimental observation using a Morlet wavelet function. The wavelet cross-spectrum is calculated for the two sets of data to construct a time–frequency phase difference map. The Box-plot, an exploratory data analysis technique, is applied to interpret the phase difference for validation purposes. In addition, wavelet time–frequency coherence is calculated using the locally and globally smoothed wavelet power spectra of the two data sets. Significance tests are performed to quantitatively verify whether the wavelet time-varying coherence is significant at a specific time and frequency point, considering uncertainties in both predicted and observed time series data. The proposed wavelet spectrum analysis approach is illustrated with a dynamics validation challenge problem developed at the Sandia National Laboratories. A comparison study is conducted to demonstrate the advantages of the proposed methodologies over classical frequencyindependent cross-correlation analysis and time-independent cross-coherence analysis for the validation of dynamic systems. & 2010 Elsevier Ltd. All rights reserved.

Keywords: Wavelet coherence analysis Wavelet cross-spectra Model validation Hypothesis test Dynamic system

1. Motivation Systematic model validation needs quantitative methods to compare model prediction with experimental observations, both under uncertainty, and to quantify the confidence in model predictive capability. During the past decade, fundamental concepts and methodologies for the validation of large-scale computational models and numerical simulations have been studied by organizations such as the United States Department of Defense [1], American Institute of Aeronautics and Astronautics [2], Advanced Simulation and Computing (ASC) program of the United Sates Department of Energy [3], and American Society of Mechanical Engineers (ASME PTC#60 [4]). Detailed discussion of model verification and validation concepts and methods can also be found in the literature (see e.g., [5–24]). Within the context of model validation, decisions need to be made to accept or reject the model with certain preferences, based upon available information and prior knowledge. In order to minimize the decision error and enhance the validation accuracy, all available data and information as well as the corresponding uncertainty should be taken into

 Corresponding author. Tel.: + 1 678 844 3690; fax: + 1 678 844 6748.

E-mail addresses: [email protected], [email protected] (X. Jiang). 0888-3270/$ - see front matter & 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.ymssp.2010.05.012

576

X. Jiang, S. Mahadevan / Mechanical Systems and Signal Processing 25 (2011) 575–590

account in the model validation. In addition, a proper decision threshold needs to be selected and the probability of decision errors may be assessed in the validation. Recently, Jiang and Mahadevan [21,23] developed Bayesian risk-based decision approaches to facilitate rational decisions in model validation under uncertainty, considering the risk of using the current model, data support for the current model, and cost of acquiring new information to improve the model. These approaches also address the conditional error probabilities of Type I error (reject a correct model) and Type II error (accept a wrong model) in the assessment of a computational model under uncertainty. Model validation of a dynamic system is an especially challenging problem where time series data are used. Two main elements of this exercise are feature extraction and quantitative assessment. The objectives of feature extraction are to reduce the dimensionality of the time series data and to improve the efficiency of model validation. Different features have been used in the validation, such as peak response, root mean square error of time series, principal components, and frequency response functions. Refer to Hemez and Doebling [25] for details of various features extracted from time history data for validation purposes. With the extracted feature, a proper validation metric is then needed as a quantitative measure of agreement between the two sets of feature data. Recently, statistical hypothesis testing-based techniques have been applied for the assessment of structural dynamics models to determine whether the prediction output matches with the real output of a dynamic system (e.g., [11,25,26]). Within the context of binary hypothesis testing for model validation, hypothesis tests are performed to test whether the difference of the features between model prediction and experimental observation is equal to zero (point hypothesis testing) or within an allowable limit (interval hypothesis testing). In both approaches, null hypothesis is defined to accept the model, while alternative hypothesis is defined to reject the model. These traditional techniques for model assessment of a dynamic system, however, are based on features extracted from either the frequency or the time domain. These methods provide only an overall perspective about the validity of the computational model by comparing two signals in the time and frequency domains separately. In particular, the time-independent features (e.g., Fourier power spectrum or energy) cannot identify how the discrepancy (e.g., amplitude or phase) varies over time. Validation techniques that use such features usually over-estimate the predictive capability of a computational model (i.e., more easily accept the model) because they do not consider the feature variation in the time domain. On the other hand, various frequency contents in a time series usually play different roles in the validation of the computational model. Validation techniques that use the frequency-invariant features (e.g., root mean square error) usually under-estimate the predictive ability of a computational model (i.e., more easily reject the model). For example, the features in dominant frequencies may play a critical role in model validation while those in other frequencies (e.g., noise) may not be important. Erroneously using noise-related features may result in wrongly rejecting a good model. In particular, the existing feature-based validation techniques for dynamic system models may be ineffective for nonlinear, stochastic, and transient dynamic behavior, where the time series representing the dynamic system is usually non-stationary. In this study, the time–frequency-based coherence (or amplitude) and phase relationships between the predicted and observed time series are used as quantitative measures to assess the validity of the computational model for a dynamic system. Coherence is the most common measure used to quantify the amplitude and phase synchrony of two processes. It is a direct measure of the correlation between the spectra of two random processes, which are represented by the observed and predicted time series data. The coherence function is derived from the normalization of the cross-spectrum by the individual spectra and measures the correlation between the signals as a function of frequency. If the processes are stationary, i.e., their spectra do not change in time, Fourier analysis can provide accurate estimates of these spectra and the corresponding coherence (e.g., [27,28]). The Fourier-based coherence analysis, however, has three main limitations: (1) it assumes stationarity of the underlying signals whereas signals in many dynamic systems may be non-stationary; (2) it does not separate the effects of amplitude and phase from each other because the signals are decomposed into a sum of sinusoidal waves with constant amplitude and phase, and (3) it cannot capture the temporal information such as transient covariance in the signals. Windowed Fourier analysis allows for non-stationarity by splitting the signal into segments [29]. Although this overcomes the assumption of global stationarity of the time series, it still requires stationarity within each segment. Therefore, a generalized time–frequency-based wavelet coherence measure is pursued in this paper to overcome the limitations of classical coherence analysis for the model validation of dynamic systems. A wavelet transform (WT) decomposes a temporal signal into a set of time-domain basis functions with various frequency resolutions. The WT is computationally similar to the windowed short-term Fourier transform (WFT). However, unlike the sine and cosine functions used in the WFT, the wavelet transform decomposes a signal into a set of orthogonal basis functions, called wavelets. The wavelets consist of a family of mathematical scaling functions used to represent a signal in both the time and frequency domains. Unlike Fourier coherence that employs a constant window width for all frequencies, the wavelet transform uses shorter windows for higher frequencies, which leads to a more ‘‘natural’’ localization (refer to Daubechies [29] for details on this topic). In general, the wavelet transform overcomes the aforementioned limitations of Fourier transform in three aspects: (1) it can capture the features of non-stationary signals due to the simultaneous time–frequency decomposition of a signal; (2) it provides both the time-varying power spectrum and the time-varying phase spectrum in the coherence analysis of non-stationary signals, and (3) it has a more accurate time–frequency resolution and more effective representation of discontinuities and transient covariance in signals due to the frequency-varying decomposition windows. In the past decade, there is a wealth of literature about the application of wavelet analysis for processing time series data. For example, continuous wavelet transform-based spectral analysis of a stationary discrete process has been described in detail by Chiann and Morettin [30]. In addition, continuous wavelet

X. Jiang, S. Mahadevan / Mechanical Systems and Signal Processing 25 (2011) 575–590

577

transformation (CWT)-based coherence analyses have been extensively applied to examine the relationships between two sets of time series in various disciplines, such as climatology [31], geophysics [32–35], meteorology [36], neuroscience [37–39], engineering [40], and physics [41,42]. In contrast to the Fourier transform-based coherence, wavelet coherence provides a quantitative measure of the correlation between two time series in time–frequency domain. It measures the linear relationship between two processes with respect to time and scale (or frequency), and therefore can be used to assess the agreement of the predicted process with the observed process. Wavelet coherence analysis is a very recent tool. To the best knowledge of the authors, the wavelet coherence analysis is first mentioned in Hudgins et al. [43] and has not been pursued for model validation of dynamic systems. In the following sections, the classical cross-correlation and cross-coherence analyses are first introduced. Then the wavelet spectrum analysis approach (including wavelet spectrum, cross-spectrum, and coherence analysis) is presented for the model validation of dynamic systems. The proposed quantitative method is illustrated with a validation challenge problem developed at Sandia National Laboratories [44]. A comparison study of the classical cross-coherence analysis with the wavelet time–frequency coherence approach is conducted on the illustrative example in the context of model validation. 2. Cross-correlation and cross-coherence The cross-correlation in time domain and cross-coherence in frequency domain have been widely used as two criteria to evaluate the similarity of two signals (e.g., [45–51]). The two criteria can be found in Jiang and Mahadevan [22] and therefore their formulations are skipped in this paper for the sake of brevity. If there is no cross-correlation or cross-coherence between the two signals (i.e., model prediction and experimental observation), then their crosscorrelation or cross-coherence is constantly equal to zero, which implies an unacceptable model with full confidence. On the other hand, if the two signals are equal, then their cross-correlation or cross-coherence is constantly equal to one, implying a perfect model. This case is nearly impossible due to the variability in both model prediction and experimental observation. The desirable result is that the cross-correlation or cross-coherence approaches unity in the first few time delays, which implies that the two signals are maximally correlated or have maximum similarity. It should be pointed out that the cross-correlation and cross-coherence analyses only provide an overall perspective about the validity of computational model through comparing two signals in the time and frequency domains separately. This approach cannot provide details about the discrepancy between the model predictions and experimental signals in both time and frequency domains simultaneously, which may easily accept (cross-correlation) or easily reject (crosscoherence) the computational model as discussed earlier. This issue can be addressed through wavelet spectrum analysis as described below. 3. Wavelet spectrum analysis method 3.1. Continuous wavelet transformation In mathematics, the continuous wavelet transform of a continuous, square-integrable function x(t) is usually defined as the integral of x(t) with a scaled and translated wavelet mother function over the time t. For the sake of efficiency, in this study the continuous wavelet transform of a given discrete time series s(tj), Ws(a,b), is calculated as the convolution of s(tj) with a scaled and translated version of the mother wavelet c(t) as follows [52]: Ws ða,bÞ ¼

N1 X

qffiffiffiffiffiffiffi  sðtj Þc ððtj bÞ=aÞ= 9a9,

a,b 2 R,

c 2 L2 ðRÞ,

ð1Þ

j¼0

where the parameters a 40 and b denote the frequency (or scale) and the time (or space) location, respectively, and R is the set of real numbers. The time series s(t) represents either the experimental observation y(tj) or the model prediction z(tj), where tj =jdt (j = 0, 1, 2, y, N  1), dt is a sampling interval and N denotes the number of data points. The superscript n denotes complex conjugate. The notation L2(R) represents the square summable space of the time series s(tj), where the superscript 2 denotes the square of the modulus of the time series. Note that by varying the parameters a and b one can construct a picture to show how the amplitude of any features varies with the scale and the time. Refer to Daubechies [29] and Holschneider [53] for a detailed discussion of the continuous wavelet transform. The wavelet function c(t) in continuous wavelet transform must be square integrable and have no energy at zero frequency. One of the most commonly used mother wavelets in practice is the Morlet wavelet, defined as [54]

cðtÞ ¼ p1=4 ei2pf0 t et

2

=2

,

ð2Þ

where the subscript j is dropped in the variable t of the wavelet function. The Morlet wavelet is a complex wave function, which is localized in time by a Gaussian envelope, i.e., exp( t2/2). Fig. 1 shows the real part of the Morlet wavelet function, a cosine wave truncated by the Gaussian envelope. Obviously, it is a locally periodic wave-train and conceptually related to windowed-Fourier analysis. Its Fourier transform is pffiffiffi 2 Fc ðf Þ ¼ p1=4 2eð2pf 2pf0 Þ =2 , ð3Þ

578

X. Jiang, S. Mahadevan / Mechanical Systems and Signal Processing 25 (2011) 575–590

1 0.8

Function value

0.6 0.4 0.2 0 -0.2 -0.4 -0.6 -0.8 -1 -4

-3

-2

Fig. 1. Morelet wavelet function.

-1

cos(2pt),

0 t

1

Exp(  t2/2),

2

3

4

Morlet wavelet.

where f0 is a non-dimensional characteristic (or center) frequency. Eq. (3) is a Gaussian function with the center of f0. The characteristic frequency f0 determines the wave numbers within the envelope. Note that the admissibility condition of R 2 the Morlet wavelet (i.e., 9cðf0 Þ9 =9f0 9df0 o þ1) is satisfied only if f0 40.8 [54]. Following Torrence and Compo [36], f0 = 0.955 (i.e., o0 = 2pf0 = 6) is used in this study. This selection has two merits for the continuous wavelet transform of a time series. First, it provides a good balance between time and frequency localizations. This value results in six waves, as is shown in Fig. 1. Second, the chosen value of f0 results in the Fourier period being approximately equivalent to the wavelet scale. According to the method of Meyers et al. [55], the relationship between the equivalent Fourier period and the wavelet scale qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi can be derived analytically for a particular wavelet function. For the Morlet wavelet, o ¼ 4pa=ðo0 þ 2þ ðo0 Þ2 Þ is obtained in which o = 2pf. Substituting o0 =6 yields the relationship o =1.033a, which makes it easy to relate the wavelet scale to the classical Fourier period. By the convolution theorem, the continuous wavelet transform of the time series s(tj) is obtained by the inverse Fourier transform as follows [36]: Ws ða,bÞ ¼

1 pffiffiffi NX a Fs ðfk ÞFc ðafk Þei2pfk b ,

ð4Þ

k¼0

where the subscript n in c indicates the complex conjugate. The function Fs(fk) represents the discrete Fourier transform (DFT) of the time series s(tj) at the frequency fk, and the parameter k is the frequency index of the Fourier transform, which corresponds to the wavelet scale a. Eq. (4) provides a convenient approach to compute the continuous wavelet transform, which only requires an inverse DFT of the product of Fs(fk) and Fc ðafk Þ. Note that the continuous wavelet transform is a complex function since the Morlet wavelet is complex. The time series can be reconstructed exactly using inverse wavelet analysis. Refer to wavelet literature (e.g., [28,35,56]) for details of inverse wavelet analysis. There are two important issues which need to be addressed regarding the continuous wavelet transform. First, a set of scales needs to be chosen for a specific wavelet function. The maximum number of scales J is empirically determined by the scale width da and the predefined smallest scale a0 using the following formula [57]: J ¼ log2 ½Ndt=ð3a0 Þ=da,

ð5Þ

aj ¼ a0 2jda ,

ð6Þ

j ¼ 0,1,. . .,J,

where N is the number of time series data points. The smaller value of da gives a finer resolution for the continuous wavelet transform of a time series. It is suggested by Torrence and Compo [36] that a value dar0.5 should be used for the Morlet function. The value da = 0.25 is taken for the example presented in this paper and found to achieve the enough fine resolution of the continuous wavelet transform. In addition, the smallest scale a0 =2dt is taken for the Morlet function [55]. In the following formulation, we drop the subscript j in the scale a for the sake of simplicity. Second, the Fourier transform of the wavelet function at each a [Fc ðafk Þin Eq. (4)] needs to be normalized to generate unit energy. Thus, the continuous wavelet transform of a time series at each scale is directly comparable to each other and to the transform of other time series. In this study, the following normalization is performed: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Fc ðafk Þ ¼ 2pa=dt Fc0 ðafk Þ, ð7Þ

X. Jiang, S. Mahadevan / Mechanical Systems and Signal Processing 25 (2011) 575–590

579

Thus, the total energy at each scale a is N 1 X

2

9Fc ðafk Þ9 ¼ N,

ð8Þ

k¼0

where 9.9 denotes the absolute value. The normalization results in the continuous wavelet transform in Eq. (4) being weighed only by the amplitude of the Fourier coefficients of the time series Fs(fk) and not by those of the wavelet function Fc ðafk Þ. Thus, the generalized CWT method is applicable for other wavelet functions as well. 3.2. Wavelet power spectrum The wavelet power spectrum, Ss(a,b), provides a measure of the time series variance at each time and each scale (period). It is defined as the square modulus of the wavelet transform 2

Ss ða,bÞ ¼ 9Ws ða,bÞ9 ,

ð9Þ

Eq. (9) is a time-varying spectral representation. It describes the energy distribution of the signal along the time and frequency range. Following the terminology of Fourier analysis, the wavelet spectrum of a time series is called either a wavelet scalogram [48] or wavelet periodogram [58]. It measures the contribution of the point at a specific time and frequency to the total energy, and quantifies the variance of the process at the certain time and frequency. The main difference between the wavelet power spectra and the Fourier power spectra for the time series is that the former identifies the frequency content along the operational time, which cannot be done with the latter. There are three critical issues regarding the wavelet spectrum analysis that need to be considered: (1) end effects of time series data, (2) normalization of the wavelet power spectrum, and (3) averaging the wavelet periodogram in both time and scale directions. First, the continuous wavelet transform at a point in time tj always contains information of neighboring data points. Thus, end effects (or errors) occur at the beginning and end of the time series. The edge effect problem has been well described and analyzed by several authors [34,36,59, 60]. The affected area is called the cone of influence (COI) [34,36]. In order to reduce the COI, Torrence and Compo [36] proposed zero padding method and Zhu [60] proposed an antisymmetric padding on both sides of the intrinsic mode functions derived from the empirical mode decomposition method. In this study, we follow Torrence and Compo’s method [36] due to its simplicity and pad the end of the time series with sufficient zeros to the next-higher power of two, and then remove them after continuous wavelet transform. Zero padding introduces discontinuities at the endpoints and decreases the amplitude near the edges at the larger scales. We choose the cone of influence as the area in which the wavelet power spectra has dropped by a factor e2 and assume that the edge effects beyond this point are negligible. Refer to Torrence and Compo [36] for more details about the cone of influence. Second, based on Eqs. (4) and (8), the expected value of Ss(a,b) is equal to N times the expected value of 9Fs(fk)92. For a white-noise time series, the expected value of its wavelet power spectrum is s2/N, where s2 is the variance of the white noise. Thus, the expected value of the CWT of a white-noise process with N data points is s2 at all a and b, and we can normalize the wavelet spectrum of any time series by 1/s2, which gives a measure of the power relative to the white noise. By doing this we can verify whether the wavelet spectrum at time tj and scale a is significant relative to the white noise through hypothesis testing. Third, averaging of the wavelet spectrum along the time–frequency domain is necessary for comparing two signals using wavelet spectrum analysis. The wavelet spectrum/coherency of the underlying processes is not calculated exactly but estimated approximately using time series data. The neighboring points in the time–frequency spectrum are correlated, resulting in the smooth wavelet periodogram and coherence map. In Fourier analysis, the transform for any time series is asymptotically Gaussian distributed and their neighboring frequencies are uncorrelated. Thus, the cross periodogram is w2 distributed and a test for significant coherency can be performed easily using an F distribution [36]. This is not the case for wavelet coherence with different scales and times, because for small scales the wavelet transform is not asymptotically Gaussian distributed and the neighboring wavelet times and scales are correlated. If a Chi-square distribution is used to model the wavelet periodogram, the distribution will have a very high variance. It was indicated in Torrence and Compo [36] that the corresponding coherence may be identical to one for any two time series even if the real coherence of the two processes is zero, for example, two completely independent Gaussian processes. To reduce the variance and produce a useful measure of coherence to compare two signals, averaging of the wavelet spectrum has been commonly used in wavelet spectrum analysis (e.g., [36,38,41]). Generally, there are two ways to smooth the wavelet spectrum along the time and frequency domains: local and global. The local smoothing allows the wavelet power spectrum (or wavelet coherence to be discussed later) to be estimated for a single trial, in particular when only a pair of predicted and observed signals is available for model validation. In such a situation, the local smoothing is carried out for the wavelet spectrum of each time series over a certain period na (timeaveraged) and with a boxcar filter with the smoothing scale width da0 over the scale range a1 to a2, (scale-averaged) as follows: 2 3 a2 na X X 2 i 2 4 ð10Þ Ss ða,bÞ ¼ ½da0 dt=ðna CÞ 9Ws ða,bÞ9 =ðas Þ5, i¼1

aj

580

X. Jiang, S. Mahadevan / Mechanical Systems and Signal Processing 25 (2011) 575–590

where the smoothing window na = fix(2a) is taken in this study, in which the function fix(.) is used to round the value toward zero. The parameter C is a constant for a chosen wavelet function. In order to avoid both over-averaging and under-averaging the wavelet spectrum, proper values should be selected for these parameters C, da0, a1, and a2. In this study, C =0.776 and da0 =0.6 are applied for the Morlet function [36], and the scale range between a1 = 0.7a and a2 = 1.3a is taken at each scale a. These values are selected because the resulting filters provide the minimal amount of smoothing necessary in both the time and scale directions for the Morlet-based continuous wavelet transform of a signal [31]. On the other hand, if both model prediction and experimental observations contain a series of repeated data sets, i.e., Y = {y1(t), y2(t), y, yK(t)} and Z= {z1(t), z2(t), y, zK(t)}, the wavelet power spectrum is analyzed by averaging across multiple data sets (global smoothing), either with or without local smoothing in a single signal. The smoothed wavelet power spectrum is used for wavelet cross-spectrum and coherence analysis as described below. 3.3. Wavelet cross-spectrum In the model validation of a dynamic system, the time–frequency representations Wy(a,b) and Wz(a,b) are obtained by using Eq. (4) for experimental observations y(tj) and prediction output z(tj), respectively. The time–frequency cross-spectrum is then obtained by Syz ða,bÞ ¼ Wy ða,bÞWz ða,bÞ:

ð11Þ

As mentioned earlier, the continuous wavelet transform of a time series using Morlet wavelet is a complex function, and therefore can be divided into the real part RfWs ða,bÞg and the imaginary part IfWs ða,bÞg. Accordingly, Eq. (11) is a complex function that can be decomposed into amplitude and phase, respectively, as follows: Ayz ða,bÞ ¼ 9Syz ða,bÞ9,

ð12Þ

Fyz ða,bÞ ¼ tan1 ½IfSyz ða,bÞg=RfSyz ða,bÞg,

ð13Þ

and Eq. (11) becomes Syz ða,bÞ ¼ Ayz ða,bÞeiFyz ða,bÞ :

ð14Þ

Model validation involves the quantitative evaluation of the relationship between the prediction output and the experimental data. The time–frequency cross-spectrum of Eq. (11) or (14) denotes the co-varying power of two processes, but it does not identify the relative dependence of two time series. Therefore, the cross-spectrum cannot be used directly to assess the validity of model prediction with regard to experimental observation. The wavelet phase spectrum represented by Eq. (13), however, provides an effective quantity to assess the model validity. It estimates the phase difference of two time series with respect to both time and scale (or Fourier period). The range of possible phase difference values is  p r f r p. If the phase difference approaches zero, then it should indicate that the model is perfectly acceptable. However, this perfect case is impossible due to the uncertainties in both the experimental observation and the model prediction. It should be noted that local smoothing can be also performed for the cross-spectrum between two signals defined in 2 Eqs. (11) or (14), with the 9Wsi ða,bÞ9 replaced by the Wy ða,bÞWz ða,bÞ in Eq. (10). Similarly, global smoothing is applied for the wavelet cross-spectrum of two signals with multiple tests. 3.4. Wavelet coherence The time–frequency coherence is defined as the squared absolute value of the smoothed cross-wavelet spectrum between y(tj) and z(tj), normalized by the squared absolute value of the smoothed wavelet power spectrum of y(tj) 2

2 ða,bÞ ¼ 9Syz ða,bÞ9 =½Sy ða,bÞ2 , Cyz

ð15Þ

where Syz ða,bÞ and Sy ða,bÞ represent the smoothed time–frequency cross-spectrum and the smoothed wavelet power spectrum, respectively. The proposed wavelet coherence is normalized with regard to only the wavelet power spectrum of experimental observations. Its purpose is to compare the predicted output with regard to the experimental observation. The coherence map provides a view of the localized correlation between two signals at each time and frequency. Similarly, the time–frequency-based coherence varies in the range [0,1]. If there is no cross-correlation between the two signals, then their coherence is constantly equal to zero, which implies that the two signals are completely unrelated or independent at time b on a scale a. On the other hand, if two signals are equal, then their cross-correlation is constantly equal to one, which implies a perfect linear relationship between y(tj) and z(tj) at the time and scale. 3.5. Significance test Significance tests are performed on the wavelet time–frequency coherence to validate the computational model for the dynamic system. The tests are used to infer whether two signals (model prediction and experimental observation) represent two coherent processes. The use of significance tests allows one to verify that the time–frequency coherences

X. Jiang, S. Mahadevan / Mechanical Systems and Signal Processing 25 (2011) 575–590

581

between two signals are statistically significant. The idea is to formulate a null hypothesis H0 ‘‘the two processes have coherence equal to one’’ (implying model is valid) and then to derive a significance test for H0. Therefore, the probability distribution of the coherence under H0 has to be estimated. If the coherence lies outside a chosen a-quantile (e.g. a = 0.05), the so-called critical value of this distribution, H0 is rejected and the coherence between two signals is significantly different from one at the (1 a) level (implying model is rejected). Now we need to derive the significance test for H0. As discussed earlier, the neighboring wavelet time–frequency power spectra are correlated, and the resulting wavelet cross-spectrum is not Chi-square distributed. Thus, the wavelet coherence does not follow an F distribution, and therefore it is difficult to perform an analytical test statistic on the wavelet coherence. Following the procedure presented in Ref. [34], Monte Carlo simulations are performed in this study to numerically estimate the distribution under H0 in five steps: (1) Generate m realizations (m= 10,000 in the example presented in this paper) of two independent red noise processes by using two univariate lage-1 autoregressive [AR(1) or Markov] models as follows:

(2) (3)

(4)

(5)

y^ n þ 1 ¼ ay y^ n þ ey ,

ð16Þ

z^ n þ 1 ¼ az z^ n þez ,

ð17Þ

where ay and az are the lag-1 autocorrelations of experimental observations y(tj) and prediction output z(tj), respectively, the initial values y^ 0 ¼ 0 and z^ 0 ¼ 0, and ey and ez are taken from Gaussian white noise with unit variance. Calculate the wavelet coherence between the two simulated random processes using Eq. (15) and the associated equations. Estimate the critical values at the (1  a) level (say 95%) from the obtained wavelet coherence for different time-scale smoothing windows and for various time and scale resolutions defined by o0. It should be pointed out that, the white noise is stationary in time, but the resulting critical values can be applied for significance tests on both stationary and non-stationary processes. This is because the normalized wavelet coherence is a local measure of two signals, which does not need to consider the stationarity of the investigated time series. Calculate the scale-independent critical value for each given time-scale smoothing window by smoothing a constant number o of neighboring scales. Note that, the critical value is time-independent automatically since the white noise is stationary. Normalize the obtained time- and scale-independent critical values by the wavelet coherence between y(tj) and z(tj) at each time-scale point. The normalized values are the critical ones for the investigated wavelet coherence. Thus, the confidence contour lines can be computed along the scale at the particular confidence level, as indicated in the example later.

In summary, a wavelet spectrum analysis approach is presented in this section for model validation of dynamic systems, using non-stationary time series data. The wavelet power spectrum is calculated for both model prediction and experimental observation using a Morlet wavelet function. Time–frequency phase difference map is constructed from the wavelet cross-spectrum of the two data sets. Further, wavelet time–frequency coherence map is created for the two data sets. Hypothesis tests at a given significance level (say 5%) are performed to test whether the wavelet coherence is one at a specific time and frequency point (i.e., to accept the model). The overall validity of the model is evaluated via the statistical exploratory analysis, as described in the illustrative example subsequently. It should be noted that Bayesian hypothesis testing may also be pursued here instead of the classical hypothesis testing. In the next section, the proposed method is illustrated for the validation of a structural dynamics model. 4. Numerical application 4.1. Problem description A structural dynamics example developed by Sandia National Laboratories [44] as a validation challenge problem is used in this study to investigate the effectiveness of the proposed wavelet spectrum analysis method for model validation. This validation problem consists of a three-mass subsystem and a practical application system (shown in Fig. 2). This subsystem is mounted on a simply supported beam by an additional weakly nonlinear connection and only vertical motions are permitted. Computational models are available to make simulation-based predictions of structural responses. Both subsystem-level and system-level data were virtually generated for illustrative purposes [44]. This problem has been previously investigated by Jiang and Mahadevan [22] using the wavelet energy-based Bayesian approach. Various solution techniques for this problem were also provided by researchers in a workshop organized by Sandia National Laboratories in 2006; these are compiled in a special issue of the journal Computer Methods in Applied Mechanics and Engineering (see e.g. [61,62]). Refer to Red-Horse and Paez [44] for details of this dynamic system problem. In this paper, the wavelet coherence analysis approach is validated with this challenge problem and demonstrated as an effective alternative solution to this problem.

582

X. Jiang, S. Mahadevan / Mechanical Systems and Signal Processing 25 (2011) 575–590

Subsystem

x12

m3 x11 xi response measurement locations

m2 x10

m1 x1

x2

x3

x4

x5

x6

35 in (0.889m) 37.5 in (0.9525m)

x7

x8

x9

Shock load f (t) in system validation

50 in (1.27m) Fig. 2. Structural dynamics challenge problem.

‘‘Weakly’’ nonlinear connection.

In the validation test of the three-mass subsystem, three shock excitations (i.e., low-, medium-, and high-level) with free decay are applied to m2 (Fig. 2) one at a time. The acceleration responses at m1, m2, and m3 in the vertical direction are obtained for every level of vibration excitation, resulting in three sets of time series data, each consisting of three acceleration responses for three masses separately. Each shock excitation input has three different amplitudes, thus resulting in an overall non-stationary dynamic system. The virtual test is repeated on 20 nominally identical subsystems (i.e., samples), resulting in 60 sets of response acceleration data, each containing 6144 data points with signal content up to 20,000 Hz (the data record time is about 0.3 s). Correspondingly, 20 sets of modal parameters for the subsystems are used as inputs to the computational model to create 60 sets of predicted response acceleration data. In this research, data collected from a high-level shock loading on m2 (Fig. 2) are used in the subsystem-level model validation for illustration purposes. In the context of system level validation, the simply supported uniform beam is subject to a single blast-type shock load as shown in Fig. 2. The vertical acceleration responses are measured at nine points along the beam as well as three masses of the associated subsystem, resulting in 12 sets of measured response acceleration time series, each having 12,288 data points with signal content up to 20,000 Hz (the data record time is about 0.6 s). Due to the high cost of building test specimens for full systems in practice, only three specimens test is performed in this study. Again, the data on m2 are used in the system-level model validation for illustration purposes. As an example, Fig. 3 shows the time series for both model prediction and experimental signal of the 5th subsystem sample and the 2nd system sample. Fig. 3a and b shows the high-level shock excitation from 0 to 0.3 s for the subsystem and the blast shock excitation from 0 to 0.6 s for the system, respectively. Fig. 3c shows the acceleration response time series plot of m2 for the subsystem sample under the input excitation shown in Fig. 3a, while Fig. 3d shows the time series plot for the system sample under the excitation shown in Fig. 3b. Fig. 3e and f shows the Fourier power spectra of the signals shown in Fig. 3c and d, respectively. The power spectrum of each time series is estimated using the fast Fourier transform with the transform length of 512. The solid and dashed lines in Fig. 3c–f represent the experimental data and the model prediction, respectively. An enlarged plot is embedded in Fig. 3c–f to visualize the local difference between two curves. Two observations are obtained from the plots in Fig. 3. First, the difference between the two sets of data at subsystem level cannot be identified graphically in both the time and frequency domains, implying that the model may be accepted in the subsystem level; while the difference at system level is significant visually, implying that the model may be rejected at the system level. Second, most of the frequency content (up to 99%) is concentrated in the range between 500 and 3000 Hz for the time series data of the subsystem and between 1000 and 3000 Hz for the system. Thus, the difference in the specific frequency ranges may play a key role in model validation. It should be noted that the time series at both the subsystem (Fig. 3c) and system (Fig. 3d) levels are non-stationary under shock (Fig. 3a) and blast (Fig. 3b) excitation loadings. For the sake of simplicity, two assumptions have been made in this challenge problem ([44]): (1) various configurations of the structural systems are assumed to be completely known and therefore all uncertainties come only from the subsystem; and (2) all experimental measurements are perfect and therefore do not contain any noise. 4.2. Classical cross-correlation and cross-coherence analysis The classical cross-correlation and the cross-coherence between the experimental signal and the predicted signal [22] are computed for each mass. As an example, Fig. 4 shows the cross-correlation and cross-coherence of m2 for the 5th

Force (unit)

X. Jiang, S. Mahadevan / Mechanical Systems and Signal Processing 25 (2011) 575–590

1000

1000

500

0

0

0

0.1

0.2

0.3

-1000

0

0.2

0.4 Time (s)

0.6

0.2

0.4 Time (s)

0.6

Acceleration (in/s2)

0

-5

3 Power spectrum

× 104

5

×104

0

-5 0

0.1

0.2 Time (s)

0

0.3

× 109

2

2 1 0

1000

2000

3000

Power spectrum

Acceleration (in/s2)

Time (s)

5

583

× 109

1 0

1000

2000

3000

0

0 0

2000 4000 6000 8000 10000 Frequency (Hz)

0

2000 4000 6000 8000 10000 Frequency (Hz)

Fig. 3. Time series of m2 for subsystem (left column) and system (right column): (a) high-level shock excitation from 0 to 0.3 s; (b) blast shock excitation from 0 to 0.6 s; (c) acceleration response data of 5th test for mass m2 under excitation (a); (d) acceleration response data of 2nd test for mass m2 under excitation (b); (e) power spectrum of signal (c) using FFT; and (f) power spectrum of signal (d) using FFT. Experiment, prediction.

subsystem test (left column) and the 2nd system test (right column). Fig. 4a and b shows the cross-correlation plots of time series shown in Fig. 3c and d, respectively, and Fig. 4c and d shows the cross-coherence of spectra in Fig. 3e and f, respectively. It is observed that the first value of both the cross-correlation (Fig. 4a) and the cross-coherence (Fig. 4c) approaches unity [r(0) = 0.977 and c(0) =0.9771] for the subsystem, implying that the experimental data support the computational model at the subsystem level; while the first cross-correlation value is smaller than one [r(0) =0.742] and the first cross-coherence is larger than one [c(0)= 1.654] for the system, implying that the experimental data do not support the model at the system level. The observation is in agreement with the visual result obtained from Fig. 3. It should be pointed out that both the frequency-independent cross-correlation analysis (Fig. 4a and b) and the time-independent cross-coherence analysis (Fig. 4c and d) cannot capture the non-stationary and transient features in the time series response data, which can be handled by wavelet time-frequency spectrum analysis effectively.

4.3. Wavelet spectrum analysis 4.3.1. Wavelet power spectrum For the subsystem time series, N = 6144, dt = 5  10  5, a0 = 2dt = 1  10  4, and da= 0.25. The largest value of scale J= 40 is obtained using Eq. (5), resulting in 41 scales ranging from 1  10  4 to 0.1024 based on Eq. (6). The equivalent Fourier period ranging from 1.033  10  4 to 0.1058 is obtained by 1.033 times the wavelet scales. Then, the continuous wavelet transform is performed using Eq. (4) where the Fourier transform of the wavelet function at each a is normalized by Eq. (7). Next, the wavelet power spectrum for every time series is calculated using Eq. (11). The wavelet spectrum for each time series is then locally smoothed and further globally smoothed for both the 20 sets of experimental data and the 20 sets of model prediction. Similarly, for the system time series, each having N =12,288 data points, the largest number of scale J =44 is obtained using Eq. (5), resulting in 45 scales ranging from 1  10  4 to 0.2048 based on Eq. (6). The equivalent Fourier period ranging from 1.033  10  4 to 0.2116 is obtained by 1.033 times the wavelet scales. The resulting wavelet spectrum for each time

X. Jiang, S. Mahadevan / Mechanical Systems and Signal Processing 25 (2011) 575–590

1

1

0.5

0.5 r (τ )

r (τ )

584

0 -0.5

0

-0.5

-1

-1 0

50

τ

100

0

50

100

τ 2

0.5

1

c (ω )

c (ω)

1

0

0

-0.5

-1 0

50

ω

100

0

50

100

ω

Fig. 4. Cross-correlation and coherence of m2 for 5th subsystem test (left column) and 2nd system test (right column): (a) cross-correlation of time series in Fig. 3c [r(0) = 0.977]; (b) cross-correlation of time series in Fig. 3d [r(0)= 0.742]; (c) cross-coherence of spectra in Fig. 3e, [c(0) = 0.9771]; and (d) crosscoherence of spectra in Fig. 3f [c(0) = 1.654].

series is locally smoothed and further globally smoothed for both the 3 sets of experimental data and the 3 sets of model prediction.

4.3.2. Wavelet cross-spectrum The wavelet cross-spectrum is calculated by substituting the obtained wavelet power spectra of y(ti) and z(ti) into Eq. (11). Their phase difference is obtained using Eq. (13). As an example, Fig. 6 shows the locally smoothed wavelet cross spectrum of the virtual experimental and model predicted acceleration response data of mass 2 for the 5th subsystem sample. As an example, Fig. 5a and b shows the cross-spectrum map and the phase difference map, respectively. The equivalent Fourier period ( =1.033 times the wavelet scale) is shown as the scale of the vertical axis in Fig. 5. It should be noted that a logarithmic scale with base 2 is used for the scale of the vertical axis for the sake of visual clarity of the illustration. The cone of influence is marked as cross-hatched regions. The cross-spectrum shown in Fig. 5a cannot identify the difference between two signals. However, the wavelet phase spectrum of the two signals shown in Fig. 5b provides useful information to validate the model prediction. It is observed that the phase difference at the period range between 1/3000 and 1/500 (corresponding to the frequency range 500–3000 Hz) falls in the range [  p/2, p/2]. In order to quantify the uncertainties of phase difference between two signals in the frequency domain, the Box-plots for time-averaged phase differences of Fig. 5b in the entire period range is shown in Fig. 6. The vertical axis represents the phase difference from  p to p, and the horizontal axis represents the period range in the logarithmic scale (corresponding to the vertical axis in Fig. 5). The Box-plot invented by Tukey [63] (also known as a box-and-whisker diagram or candlestick chart) is a convenient graphical approach to depict the five-number summary, including the minimum, lower quartile (25%), median, upper quartile (75%), and maximum value, and to indicate the outliers in the data. The outlier is defined as any data observation which lies more than 1.5nIQR (inter-quartile range) lower than the lower quartile (25%) or 1.5nIQR higher than the upper quartile (75%), in which the IQR is calculated by subtracting the lower quartile from the upper quartile. The Box-plot is an important exploratory data analysis technique that is able to visually show different types of populations without any assumptions of statistical distribution. Two observations are obtained from the Box-plots shown in Fig. 6. First, the median value of the phase difference is relatively small along the period range. In particular, the phase differences in the period range between 1/3000 and 1/500 are smaller than p/3, implying that the model is acceptable at the subsystem-level based on the phase difference analysis of a single trial data of m2. Second, a number of outliers are present in the region with high frequency values (larger than 5000 Hz), implying that the corresponding time series data may be affected by additional noise.

Period

X. Jiang, S. Mahadevan / Mechanical Systems and Signal Processing 25 (2011) 575–590

1/16 1/32 1/64 1/128 1/256 1/512 1/1024 1/2048 1/4096 1/8192

COI

Period

0

COI

0.05

0.1

0.15 Time (s)

0.2

0.25

585

64 16 4 1 1/4 1/16 1/64

0.3

π

1/16 1/32 1/64 1/128 1/256 1/512 1/1024 1/2048 1/4096 1/8192

COI

COI π/2 0 -π/2

0

0.05

0.1

0.15 Time (s)

0.2

0.25

0.3

-π

Fig. 5. Locally smoothed wavelet cross spectrum of virtual experiment and model predicted acceleration response data of mass 2 for 5th subsystem sample: (a) cross spectrum map (Eq. (11)) and (b) phase difference map (Eq. (13)). The cone of influence is marked as cross-hatched regions.

π

Phase difference

2π/3 π/3 0 -π/3 -2π/3 -π 1/10000

1/2500

1/893

1/316 1/94 Period (logarithm)

1/24

1/10

Fig. 6. Box-plot for time-averaged phase difference of Fig. 5b.

4.3.3. Wavelet coherence analysis 4.3.3.1. Subsystem level. Wavelet time–frequency-based coherence is calculated by substituting both the wavelet power spectrum and the wavelet cross-spectrum into Eq. (15). As an example, Fig. 7 shows the wavelet coherence map of virtual experiment and model predicted acceleration response data of mass 2 for the 5th subsystem sample with local smoothing. The critical value for each time–frequency coherence is estimated using the procedure described previously. The 95% confidence level for the wavelet coherence is shown as a thick contour in this figure. The areas of strong coherence in time and period scales between two data sets are shown inside the contour line. The cone of influence is marked as crosshatched regions. It is observed that the areas of strong coherence in time and period scales between the two data sets fall in the period range between 1/10,000 and 1/500 through the entire operation time, and the zero coherence is shown in other

586

X. Jiang, S. Mahadevan / Mechanical Systems and Signal Processing 25 (2011) 575–590

1

1/16 COI

COI

1/32

0.8

1/64 Period

1/128

95% Confidence

0.6

1/256 1/512

0.4

1/1024 1/2048

0.2

1/4096 1/8192

0 0

0.05

0.1

0.15 Time (s)

0.2

0.25

0.3

Period

Fig. 7. Wavelet coherence map of virtual experiment and model predicted acceleration response data of mass 2 for 5th subsystem sample with local smoothing. The 95% confidence level is shown as a thick contour. The cone of influence is marked as cross-hatched regions.

1/16 1/32 1/64 1/128 1/256 1/512 1/1024 1/2048 1/4096 1/8192

COI

0.8 0.6 95% Confidence 0.4 0.2 0 0

Period

1

COI

1/16 1/32 1/64 1/128 1/256 1/512 1/1024 1/2048 1/4096 1/8192

0.05

0.1

0.15 Time (s)

0.2

0.25

0.3

π

COI

COI

π/2 0 -π/2 -π 0

0.05

0.1

0.15 Time (s)

0.2

0.25

0.3

Fig. 8. Wavelet coherence of virtual experiment and model predicted acceleration response data of mass 2 for all 20 sets of subsystem samples with local and global smoothing: (a) coherence and (b) phase difference.

time–frequency domain. The result indicates that the model is acceptable at subsystem level, based on the wavelet time–frequency coherence analysis using a set of single trial data of m2. Fig. 8 shows the wavelet coherence and phase difference maps of the virtual experiment and model predicted acceleration response data of mass 2 for all 20 sets of subsystem samples with both local and global smoothing. In the coherence map (Fig. 8a), the 95% confidence level is shown as a thick contour. The areas of strong coherence in time and period scales between two data sets are shown inside the contour line. The cone of influence is marked as cross-hatched regions in both the wavelet coherence and phase difference maps. Comparing with the coherence map with only locally smoothed wavelet spectra (Fig. 7), it is observed that the areas of strong coherence in time and period scales between two data sets fall in the period range between 1/10,000 and 1/1000(1/500 in Fig. 7) through the entire operation time. The change of the period range is understandable because each experimental observation contains different frequency content due to uncertainty, and the global smoothing over all the 20 sets of time series includes the main frequency content in all scenarios. Therefore, the result indicates that the model is acceptable at subsystem level, based on the wavelet time-frequency coherence map using multiple data of m2. Fig. 9 shows the Box-plots for time-averaged phase differences of Fig. 9b in the entire period range. Similarly, the median value of the phase difference is relatively small along the period range. In particular, all the phase difference values

X. Jiang, S. Mahadevan / Mechanical Systems and Signal Processing 25 (2011) 575–590

587

π Maximum value

2π/3

Phase difference

Outlier

75% Quartile

π/3 Median 0 -π/3 -2π/3 25% Quartile

Minimum value

-π 1/10000

1/2500

1/893 1/316 Period (logarithm)

1/94

1/24

1/10

Period

Fig. 9. Box-plot for time-averaged phase difference of Fig. 8b.

1/16 1/32 1/64 1/128 1/256 1/512 1/1024 1/2048 1/4096 1/8192

COI 0.8 0.6 95% Confidence

0.4 0.2 0

0

Period

1

COI

1/16 1/32 1/64 1/128 1/256 1/512 1/1024 1/2048 1/4096 1/8192

0.05

0.1

0.15 Time (s)

0.2

0.25

0.3

π

COI

COI

π/2 0 -π/2

0

0.05

0.1

0.15 Time (s)

0.2

0.25

0.3

-π

Fig. 10. Wavelet coherence of virtual experiment and model predicted acceleration response data of mass 2 for all three system sample tests with local and global smoothing: (a) coherence and (b) phase difference.

in the period range between 1/3000 and 1/500 are smaller than p/4, implying that the model is acceptable at the subsystem level based on the phase difference analysis of multiple data of m2. Compared with Fig. 6 (only local smoothing for the single trial), Fig. 9 shows a smaller variation of the phase difference (narrow interval between 25% and 75% quartiles) in the principal period range (i.e., between 1/3000 and 1/500) due to the global smoothing.

4.3.3.2. System level. Fig. 10 shows the wavelet coherence and phase difference maps of the virtual experiment and model predicted acceleration response data of mass 2 for all three sets of system samples with both local and global smoothing. Similarly, the 95% confidence level is shown as a thick contour. The cone of influence is marked as cross-hatched regions in both the wavelet coherence and phase difference maps. It is observed that the areas of strong coherence in time and period

588

X. Jiang, S. Mahadevan / Mechanical Systems and Signal Processing 25 (2011) 575–590

π 75% Quartile

Phase difference

2π/3

Maximum value

π/3

Outlier Median

0 -π/3 -2π/3

Minimum value 25% Quartile

-π 1/10000

1/2500

1/893

1/316

1/94

1/24

1/10

1/4.9

Period (logarithm) Fig. 11. Box-plot for time-averaged phase difference of Fig. 10b.

scales between two data sets mostly fall in the period range between 1/500 and 1/16 through the entire operation time, and the zero coherence is observed in the principal frequency (period) range, i.e., 1000–3000 Hz shown in Fig. 3f. The result indicates that the model is unacceptable in the principal frequency range (zero coherence) at the system level, based on the wavelet time–frequency coherence map using multiple trials of m2. Fig. 11 shows the Box-plots for time-averaged phase differences of Fig. 11b in the entire period range. Unlike Fig. 10, the median value of the phase difference in Fig. 11 is relatively large in the period range between 1/100 and 1/15. The phase difference along the entire period has a very large variation. The results indicate that the model is unacceptable at the system level based on the phase difference analysis of multiple trials data of m2.

5. Concluding remarks To the best knowledge of the authors, no published article has been found to develop the wavelet spectrum analysis approach for model validation. This paper intends to be the first that systematically presents a generalized wavelet spectrum analysis approach to model validation of dynamic systems. Continuous wavelet transform is performed on both model prediction and experimental observation in the form of time series, using a Morlet wavelet function. The wavelet cross-spectrum is calculated for the two data sets to construct a time–frequency phase difference map. The Box-plot is applied to interpret the phase difference for validation purposes. Wavelet coherence is further calculated to assess the validity of the computational model at each time and frequency. Significant tests are performed to quantitatively verify whether the wavelet coherence is significant at a specific time and frequency point. The proposed wavelet spectrum analysis approach is illustrated with a dynamics validation challenge problem developed by Sandia. The traditional techniques for model assessment of a dynamic system are based on features extracted from either the frequency or time domain. These methods provide only an overall perspective about the validity of computational model through comparing two signals in the time and frequency domains separately. They may over-estimate (more easily accept) or under-estimate (more easily reject) a computational model because they do not consider the significant timevarying or frequency-dependent information in the signals. The proposed time–frequency spectrum analysis approach can more reliably identify how the discrepancy (e.g., amplitude or phase) varies over time and what kind of discrepancy plays a key role on the validation of the computational model. Therefore, the proposed approach provides more evidence about the model validity. In comparison to the classical cross-correlation analysis in the time domain and cross-coherence analysis in the frequency domain, the proposed method has demonstrated two unique features, which are advantageous for the validation of dynamic system. First, it can be used to effectively compare two sets of time series data in both time and frequency domains simultaneously, thus quantitatively identifying the discrepancy in the time-varying phase and amplitude of two time series data, which is needed for effective model development. Second, no assumption of stationarity is made for decomposing the time series using the wavelet spectral analysis approach; thus the method is able to compare two sets of non-stationary time series and capture the transient features, which cannot be handled by the classical Fourier-based coherence analysis. Therefore, the proposed time–frequency coherence analysis approach has a strong potential for

X. Jiang, S. Mahadevan / Mechanical Systems and Signal Processing 25 (2011) 575–590

589

effective application in model validation of complicated real dynamic systems that usually generate non-stationary time series response. In future, the proposed wavelet coherence analysis-based validation approach may be extended to more complicated real-world examples. A modified Morlet wavelet may improve the time–frequency resolution of the wavelet transform of a time series. Besides, an anti-symmetric padding method may be applied to reduce the edge effect in the wavelet coherence analysis of time series. A high-order autoregressive model may also be investigated in the hypothesis testing for quantitative model assessment. In addition, instead of classical hypothesis testing, the Bayesian inference approach and risk-based decision approach may be pursued to verify whether the wavelet coherence approaches unity at a specific time and frequency point. Furthermore, multiple frequency points may be compared simultaneously to evaluate the overall validity of a computational model by using either classical or Bayesian multivariate hypothesis testing techniques.

Acknowledgments This study was supported by funds from Sandia National Laboratories, Albuquerque, New Mexico (Contract no. BG-7732, project monitors: Dr. Thomas L. Paez, Dr. Laura P. Swiler, and Dr. Martin Pilch). The support is gratefully acknowledged. Parts of the MATLAB codes used in this study to generate wavelet time-frequency coherence map were provided by A. Grinsted (http://www.pol.ac.uk/home/research/waveletcoherence/), which is also highly appreciated. We also thank the anonymous reviewers for their thoughtful and constructive comments on this paper. References [1] DOD. Verification, Validation, and Accreditation (VV&A) Recommended Practices Guide. Department of Defense, Alexandria, VA 1996, /http://vva. dmso.mil/S. [2] AIAA, Guide for the Verification and Validation of Computational Fluid Dynamics Simulations, American Institute of Aeronautics and Astronautics AIAA-G-077-1998, Reston, VA, 1998. [3] DOE, Advanced Simulation and Computing (ASC) Program Plan, Office of Advanced Simulation & Computing, DOE/NNSA NA-114. /http://www. sandia.gov/NNSA/ASC/pdfs/ASC-PP-Final-Jan08.pdfS. [4] ASME. Guide for Verification and Validation in Computational Solid Mechanics, American Society of Mechanical Engineers, ASME V&V 10, 2006. [5] O. Balci, Verification, validation, and accreditation of simulation models, in: Proceedings of the 29th Conference on Winter Simulation, Atlanta, GA, December 1997. [6] P.J. Roache, Verification and validation in computational science and engineering, Science and engineering, Hermosa Publishers, Albuquerque, NM, 1998. [7] T.G. Trucano, R.G. Easterling, K.J. Dowding, T.L. Paez, A. Urbina, V.J. Romero, B.M. Rutherford, R.G. Hills. Description of the Sandia Validation Metrics Project, Technical Report SAND 2001-1339, Sandia National Laboratories, Albuquerque, NM, 2001. [8] R.G. Hills, T.G. Trucano. Statistical Validation of Engineering and Scientific Models: A Maximum Likelihood Based Metric, Technical Report SAND 2001-1783, Sandia National Laboratories, Albuquerque, NM, 2002. [9] W.L. Oberkampf, T.G. Trucano, Verification and validation in computational fluid dynamics, Progress of Aerospace Science 38 (2002) 209–272. [10] W.L. Oberkampf, T.G. Trucano. Design of and Comparison with Verification and Validation Benchmarks, Technical Report SAND 2006-5376C, Sandia National Laboratories, Albuquerque, NM, 2006. [11] T.L. Paez, A. Urbina, Validation of mathematical models of complex structural dynamic systems, in: Proceedings of the Ninth International Congress on Sound and Vibration, Orlando, FL, July 2002. [12] R.G. Hills, I. Leslie. Statistical Validation of Engineering and Scientific Models: Validation Experiments to Application, Technical Report SAND 20030706, Sandia National Laboratories, Albuquerque, NM, 2003. [13] W.L. Oberkampf, T.G. Trucano, C. Hirsch. Verification, Validation, and Predictive Capabilities in Computational Engineering and Physics, Technical Report SAND 2003-3769, Sandia National Laboratories, Albuquerque, NM, 2003. [14] I. Babuska, J.T. Oden, Verification and validation in computational engineering and science: basic concepts, Computer Methods in Applied Mechanics and Engineering 193 (2004) 4057–4066. [15] K.J. Dowding, R.G. Hills, I. Leslie, M. Pilch, B.M. Rutherford, M.L. Hobbs. Case Study for Model Validation: Assessing a Model for Thermal Decomposition of Polyurethane Foam, Technical Report SAND 2004-3632, Sandia National Laboratories, Albuquerque, NM, 2004. [16] K. Campbell, A brief survey of statistical model calibration ideas, in: Proceedings of International Conference on Sensitivity Analysis of Model Output, Santa Fe, NM, March 2004. [17] P. Sargent, Validation and verification of simulation models, in: Proceedings of the 2004 Winter Simulation Conference, Washington, D.C., December 2004. [18] S. Mahadevan, R. Rebba, Validation of reliability computational models using Bayes networks, Reliability Engineering and System Safety 87 (2005) 223–232. [19] R. Rebba, S. Mahadevan, Model predictive capability assessment under uncertainty, AIAA Journal 44 (2006) 2376–2384. [20] W.L. Oberkampf, M.F. Barone, Measures of agreement between computation and experiment: validation metrics, Journal of Computational Physics 217 (2006) 5–36. [21] X. Jiang, S. Mahadevan, Bayesian risk-based decision method for model validation under uncertainty, Reliability Engineering and System Safety 92 (2007) 707–718. [22] X. Jiang, S. Mahadevan, Bayesian wavelet method for multivariate model assessment of dynamic systems, Journal of Sound and Vibration 312 (2008) 694–712. [23] X. Jiang, S. Mahadevan, Bayesian validation assessment of multivariate computational models, Journal of Applied Statistics 35 (2008) 49–65. [24] X. Jiang, S. Mahadevan, Bayesian inference method for model validation and confidence extrapolation, Journal of Applied Statistics doi:10.1080/ 02664760802499295. [25] F.M. Hemez, S.W. Doebling, Validation of structural dynamics models at Los Alamos national laboratory, in: Proceedings of the 41st AIAA/ASME/ ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference, Atlanta, GA, April 2000. [26] S.W. Doebling, T.A. Butler, J.F. Schultze, F.M. Hemez, L.M. Moore, M.D. McKay, Validation of transient structural dynamics simulations: an example, in: Proceedings of SAMO 2001: Third International Conference on Sensitivity Analysis of Model Output, Madrid, Spain, June 2001. [27] T. Bullock, M. McClune, J. Achimowicz, V. Iragui-Madoz, R. Duckrow, S. Spencer, EEG coherence has structure in the millimeter domain: subdural and hippocampal recordings from epileptic patients, Electroencephalography and Clinical Neurophysiology 95 (1995) 161–177.

590

X. Jiang, S. Mahadevan / Mechanical Systems and Signal Processing 25 (2011) 575–590

[28] S. Weiss, P. Rappelsberger, EEG coherence within the 13–18 Hz band as a correlate of a distinct lexical organisation of concrete and abstract nouns in humans, Neuroscience Letters 209 (1996) 17–20. [29] I. Daubechies, Ten Lectures on Wavelets, Society for Industrial and Applied Mathematics, Pennsylvania, 1992. [30] C. Chiann, P.A. Morettin, A wavelet analysis for time series, Journal of Nonparametric Statistics 10 (1998) 1–46. [31] C. Torrence, P.J. Webster, Interdecadal changes in the ENSO-Monsoon system, Journal of Climate 12 (1999) 2679–2690. [32] W. Popin´ski, W. Kosek, H. Schuh, M. Schmidt, Comparison of two wavelet transform coherence and cross-covariance functions applied on polar motion and atmospheric excitation, Studia Geophysica Et Geodaetica 46 (2002) 455–468. [33] A. Grinsted, J.C. Moore, S. Jevrejeva, Application of the cross wavelet transform and wavelet coherence to geophysical time series, Nonlinear Processes in Geophysics 11 (2004) 561–566. [34] D. Maraun, J. Kurths, Cross wavelet analysis: significance testing and pitfalls, Nonlinear Processes in Geophysics 11 (2004) 505–514. [35] C.P. Stark, J. Stewart, C.J. Ebinger, Wavelet transform mapping of effective elastic thickness and plate loading: validation using synthetic data and application to the study of southern African tectonics, Journal of Geophysical Research 108 (2003) ETG3.1-ETG3.19. [36] C. Torrence, G.P. Compo, A practical guide to wavelet analysis, Bulletin of the American Meteorological Society 79 (1998) 61–78. [37] J.-P. Lachaux, A. Lutz, D. Rudrauf, D. Cosmelli, M. Le Van Quyen, J. Martinerie, F. Varela, Estimating the time-course of coherence between single-trial brain signals: an introduction to wavelet coherence, Neurophysiologie Clinique 32 (2002) 157–174. [38] Y. Zhan, D. Halliday, P. Jiang, X. Liu, J. Feng, Detecting time-dependent coherence between non-stationary electrophysiological signals—A combined statistical and time–frequency approach, Journal of Neuroscience Methods 156 (2006) 322–332. [39] K. Keissar, L.R. Davrath, S. Akselrod, Coherence analysis between respiration and heart rate variability using continuous wavelet transform, Philosophical Transactions of the Royal Society A 367 (2009) 1393–1406. [40] K. Gurley, T. Kijewski, A. Kareem, First- and higher-order correlation detection using wavelet transforms, ASCE Journal of Engineering Mechanics 129 (2003) 188–201. [41] D.S. Bloomfield, R.T. James McAteer, B.W. Lites, P.G. Judge, M. Mathioudakis, F.P. Keenan, Wavelet phase coherence analysis: application to a quietsun magnetic element, The Astrophysical Journal 617 (2004) 623–632. [42] D. Maraun, J. Kurths, M. Holschneider, Nonstationary Gaussian processes in wavelet domain: synthesis, estimation, and significance testing, Physical Review E 75 (2007) 016707.1-016707.14. [43] L. Hudgins, C.A. Friehe, M.E. Mayer, Wavelet transforms and atmopsheric turbulence, Physical Review Letter 71 (1993) 3279–3282. [44] J.R. Red-Horse, T.L. Paez, Model validation challenge problem: structural Dynamics Application, Computer Methods in Applied Mechanics and Engineering (special issue) 197 (2008) 2578–2584. [45] J. Kittler, Mathematics methods of feature selection in pattern recognition, International Journal of Man-Machine Studies 7 (1975) 609–637. [46] M.H. Chen, D. Lee, T. Pavlidis, Residual analysis for feature detection, IEEE Transactions on Pattern Analysis and Machine Intelligence 13 (1991) 30–40. [47] Y.S. Cho, S.B. Kim, E.L. Hixson, E.J. Powers, A digital technique to estimate second-order distortion using higher order coherence spectra, IEEE Transactions on Signal Processing 40 (1992) 1029–1040. [48] O. Rioul, P. Flandrin, Time-scale energy distributions: a general class extending wavelet transforms, IEEE Transaction on Signal Processing 40 (1992) 1746–1757. [49] D.S. Ornstein, B. Weiss, Entropy and data compression schemes, IEEE Transactions on Information Theory 39 (1993) 78–83. [50] Y. Wu, R. Du, Feature extraction and assessment using wavelet packets for monitoring of machining processes, Mechanical System and Signal Processing 10 (1996) 29–53. [51] R. Yan, R.X. Gao, An efficient approach to machine health diagnosis based on harmonic wavelet packet transform, Robotics and Computer-Integrated Manufacturing 21 (2005) 291–301. [52] C.K. Chui, An Introduction to Wavelets, Academic Press, Inc., San Diego, CA, 1992. [53] M. Holschneider, Wavelets: An Analysis tool, Oxford University Press, London, 1998. [54] P. Goupillaud, A. Grossmann, J. Morlet, Cycle-octave and related transforms in seismic signal analysis, Geoexploration 23 (1984) 85–102. [55] S.D. Meyers, B.G. Kelly, J.J. O’Brien, An introduction to wavelet analysis in oceanography and meteorology: with application to the dispersion of Yanai waves, Monthly Weather Review 121 (1993) 2179–2186. [56] D.B. Percieval, A.T. Walden, Wavelet Methods for Time Series Analysis, Cambridge University Press, 2000. [57] M. Farge, Wavelet transforms and their applications to turbulence, Annual Review of Fluid Mechanics 24 (1992) 395–458. [58] G.P. Nason, R. von Sachs, G. Kroisandt, Wavelet processes and adaptive estimation of the evolutionary wavelet spectrum, Journal of the Royal Statistical Society, Series B 62 (2000) 271–292. [59] D. Zheng, B.F. Chao, Y. Zhou, N. Yu, Improvement of edge effect of the wavelet time-frequency spectrum: application to the length-of-day series, Journal of Geodesy 74 (2000) 249–254. [60] X. Zhu, Reduction of edge effects in the continuous wavelet transform by the empirical mode decomposition method observation, theory and modeling of atmospheric variability, selected papers of Nanjing Institute of Meteorology Alumni in Commemoration of Professor Jijia Zhang, 2002, pp. 604–612. [61] R. Ghanem, A. Doostan, J. Red-Horse, A probabilistic construction of model validation, Computer Methods in Applied Mechanics and Engineering (special issue) 197 (2008) 2585–2595. [62] T. Hasselman, G. Lloyd, A top-down approach to calibration, validation, uncertainty quantification and predictive accuracy assessment, Computer Methods in Applied Mechanics and Engineering (special issue) 197 (2008) 2596–2606. [63] J.W. Tukey, Exploratory Data Analysis, Addison-Wesley, Reading, MA, 1977.