On the use of weighted least squares for time domain modal parameter identification

An approach for the use of the Weighted Least Squares method for time domain modal parameter identification is introduced. The approach enables standard identification methods to be employed in the usual manner, except that sections of the time records with high signal-noise ratios are given a greater emphasis in the identification process. It is shown how it is possible to apply the method to a wide range of applications, including the analysis of impulse response, input-output and time varying data. A number of different weighting schemes are explored, based upon exponential weighting and the envelope of the measured data. A series of simulated tests demonstrate how the performance of conventional time domain algorithms can be improved significantly with little extra computational effort though the use of a Weighted Least Squares scheme.


Introduction
Modal testing is performed to determine experimentally the structural modal properties, often with the goal of finite element model validation and updating in mind.A number of recent applications have used system identification techniques to monitor the condition of aircraft, cars, railways, etc.The science of linear parameter identification is relatively mature [1,10] and reasonably widespread in industry and there is a wide range of different algorithms that can be used to estimate the modal parameters.These methods can be categorised as to whether the identification is performed using frequency or time domain data.It is the latter type of method that will be considered in this work.
Time domain modal parameter identification methods were originally developed as time series analysis techniques for the analysis of such diverse fields as finance, geophysics and astronomy.Their use by the modal community has been widespread since the 1980s, however, a number of difficulties in implementation (determination of model order, length of data set to use, selection of poles) necessitate a fair degree of expertise in order to obtain good parameter estimates.
Time domain methods can be used in a variety of ways, depending primarily upon what type of experiment has been performed and what data is available.By far the most common approach is to curve-fit impulse responses obtained either by: -measuring the response to some impact to the structure, -via the inverse FFT of the Frequency Response Function (FRF), or -from the auto-correlation of the response to an unmeasured random excitation ISSN 1070-9622/04/$17.00 2004 -IOS Press and the authors.All rights reserved One difficulty in the use of an impulse response is that the noise-signal ratio gets progressively worse as the decay progresses.Conventional identification methods, such as Least Squares Complex Exponential (LSCE) and Polyreference [8], ERA [9], ITD [12], etc. are based upon the Least Squares method [2] and do not take this effect into account although it can be severe when noisy data (such as flight flutter test data [5] is analysed.The engineer simply has to define how much of the impulse response is used for the analysis, and this is set by using as much data as possible whilst not including data in the noise floor.In a lot of commercial software the user does not have any idea how much of the impulse response function is being used in the curve-fit, let alone being able to define it.All of the above approaches place an equal weight on all parts of the decay, and consequently the performance of the identification will be degraded the greater the amount of decay that is used.
The key tool in the modal analyst's armoury to overcome biased estimates is that of model order overspecification [2].Rather than trying to define beforehand how many modes there are in the data, the user simply estimates the modal parameters with an increasing number of modes in the model.Through the use of stability plots [5], the engineer is able to not only determine how many modes that there are in the data, but also to eliminate any bias that will occur in the modal parameter estimates when there is noise present on the data.When considering noisy data, a significant amount of overspecification is required to eliminate any bias, therefore increasing the complexity of distinguishing system modes from spurious modes.The system identification literature is full of methods that are superior statistically to the Least Squares approach, however, the use of such methods has not become widespread for Modal Analysis, due to the high order, lightly damped systems that are usually considered.
The same argument can be applied to the analysis of input -output data sets.A number of approaches have considered the analysis of such data sets to identify modal parameters, though far fewer than the impulse response methods described above.It is usual to employ some type of autoregressive moving average (ARMA) or ARMAX (with exogeneous inputs) model and the parameters can be estimated using Least Squares or some other more appropriate method to eliminate the bias (e.g.Instrumental Variables, Maximum Likelihood, etc.) If the excitation sequence is random, and particularly if a random sequence possesses a good crest factor, then the signal-noise ratio can be thought of as remaining more or less constant throughout the data acquisition.However, there are many cases when a "chirp" signal is used, and this will result in sections of the data away from resonance where the signal noise ratio will be poor, particularly if a linear or log sweep is used.
Other time domain approaches where the same arguments apply are the analysis of stochastic data sets (unmeasured input) and on-line identification of time varying systems.A relatively recent advance in system identification has been the advent of subspace methods [13] which can be used either for the analysis of stochastic data sets (the input is unmeasured and assumed to be random), or for the identification of the full state-space model given measured input and output data.
In this paper, the Weighted Least Squares method is applied to a number of time domain modal parameter estimation methods in order to give a greater emphasis to those data points with a superior signal noise ratio.The use of exponential weighting and also the Hilbert Transform to define the weighting schemes is investigated.Simulated studies show how traditional modal parameter estimation approaches can be improved with little increase in computational complexity.

Impulse response based methods
In order to demonstrate the approach proposed here, the response of a single-input single-output (SISO) two Degree of Freedom (DOF) system is identified using the Least Squares Complex Exponential (LSCE) method [8].A similar approach could be used for a wide range of time domain methods.The key step of the LSCE method is to find the coefficients of the auto-regressive difference equation from measured data y j using the well known Least Squares (LS) method to minimise the squares of the residual ε.
All of the terms in Eq. ( 1) are scalar as a SISO system is being considered, however, the equation can be reformulated using matrices in order to model both the single-input multi-output (SIMO) and multi-input multi-output (MIMO) cases [3].Although LS is the most widely used approach, other more sophisticated system identification techniques These more advanced methods can be shown to give better estimates with less model order overspecification being required on simple simulated data sets.However, when applied to large order real structures, the advantages of using more sophisticated algorithms are less clear cut and they can be difficult to implement.Consequently, most work in the modal community has used the LS based methods.Other techniques such as ERA and Subspace methods have been developed to produce a minimum order realisation, through use of the Singular Value Decomposition (SVD), which reduces significantly the number of estimated modes that have to be considered following an analysis.However, LS is still the fundamental element of such techniques and a large amount of model order overspecification is employed in the solution to reduce the bias.
The difference equation coefficients are found by writing Eq. ( 1) for different time instants and then solving      y 5 y 6 . . .
where the a i terms are the difference equation coefficients that have to be estimated based upon measured time data y j in order to find the natural frequencies and damping ratios.This is done in a Least Squares fashion, minimising the sum of the squares of the residuals ε j such that θ = (φ T φ) −1 φ T Y (3) Such an approach is the fundamental element in algorithms such as Least Squares Complex Exponential, Polyreference, ERA, ITD, some stochastic subspace methods, etc.Other methods (e.g.ERA/DC [11], subspace methods, etc.) make use of the instrumental variables technique in some form to enable the bias that arises to be eliminated.Should the MIMO case be considered, then the terms in Eq. ( 2) simply become matrices whose dimensions depend upon the number of measurement stations and input cases.
The key point to note is that all of the above methods apply an equal weighting to the points in the decay response.Therefore if a large amount of the decay is chosen for the analysis, a lot of the data considered will have a poor signal noise ratio.Should the data be truncated, so that only part of the decay is considered, then it is likely that the estimates for the lower frequency modes will be biased as relatively few cycles for these modes are measured.Figure 1 shows a typical noisy impulse response and it can be seen that the noise floor is not negligible.

A weighted least squares complex exponential algorithm
The Weighted Least Squares (WLS) approach is a method that can be applied when greater or lesser emphasis needs to be applied to points in a curve-fit without actually removing the points from the fit altogether.It is often applied in the fitting of experimental data where there is less confidence in some of the points.Now, this is exactly the situation that we are considering here, and we have an idea as to how much confidence we have in the data by simply examining the decay visually and noting how quickly it gets absorbed into the noise floor.
For the WLS version of LSCE, we simply take the set of equations above and before applying Least Squares, multiply each row of the data matrices by a constant α i , such that where 0 α i 1.If we have a lot of confidence in a measured row of data (i.e.near the start of the decay) we set α i at, or close to, unity.When there is less confidence in the data (i.e. when the decay has reached the noise floor, or the response is small) α i is set close to zero.The equations are then solved in the same Least Squares manner as before.Note that there is no alteration to the solution of the equations, and no damping is added to the solution.
Even when parts of the data are given a smaller weighting, they still have some influence on the curve-fit which will benefit the estimation of the low frequency modes.
In order to make this process robust, some form of automation needs to be introduced for the definition of the weights.

Exponential weighting approach
Due to the exponential damping characteristic of an impulse response, it makes sense to set where ∆t is the sampling interval.Figure 2 shows a typical set of these curves for 0 β 1.It can be seen how an increase in the value of β quickly puts a significant decrease in the exponential function.Having applied the weighting, the solution is then weighted in a much greater manner towards the early part of the decay.Factor β needs to be set so that virtually zero weighting is given to the data where the decay meets the noise floor.
One problem with the exponential weighting is that the best value of β is data dependent.Consequently, some approach that can be used in some automated sense can be thought of as being superior.

Envelope function approach
An approach to automate the weighting process is to determine the "envelope" of the response via the Hilbert Transform.As long as the excitation is such that leakage is not a problem, then the envelope can be calculated by: 1. transforming time signal y into the frequency domain 2. introducing a 90 • phase shift through multiplication by j 3. transformation of this new signal back into time to give y H 4. the envelope is found as y 2 + y 2 H Such an approach also enables that straightforward extension of the WLS approach to all mathematical model types rather than just impulse responses.Figure 3 shows how the envelope can be calculated for an impulse response.

Application to other mathematical models
Although derived above for the LSCE method, the WLS approach can be applied to virtually any other time domain system identification algorithm and model.

Hankel matrix based impulse response methods
The fundamental element of methods such as ERA, ERA/DC and the Subspace techniques are based upon the relation between two different data matrices.For instance, in the ERA method, two data matrices are defined such that and the frequencies and damping ratios can determined by finding the Singular Value Decomposition of H 0 such that and then solving the eigensolution of the system matrix This solution can be weighted in a WLS sense simply by multiplying each of the columns of data matrices in Eq. ( 6) by the appropriate weighting factor.The same approach can be used for the ERA/DC method.

Input/Output data
Such data is usually modelled by some form of ARMAX model such that where the unknowns (a i and b i ) are found in a similar matrix approach to that of Eqs (1)(2)(3)(4).Such a formulation naturally lends itself to the WLS approach, with the weighting depending upon the output measurements (assuming that a constant level input force has been used).

On-line analysis
It is sometimes necessary to employ a recursive formulation of Eq. ( 9), particularly if the system is time varying when some form of forgetting factor must be used to emphasis the more recent data [4,7].Writing Eq. ( 9) in the form y = Φθ + ε (10) then the recursive solution for the nth time step is well known as where with The WLS weighting can be applied in the same manner as above, with Eqs (11-13) being altered such the data sequences for each new time step (β n+1 , χ n+1 , y n+1 ) are multiplied by the required weighting parameter.

Implementation on simulated data
Consider the impulse response of a 2 DOF system in Fig. 1 corrupted with 10% measurement noise.The effect of varying the weighting factor β was investigated via a Monte Carlo approach considering 30 different noise samples with the same signal noise.95% confidence bands of the frequency and damping ratio estimates without using model order overspecification are shown in Fig. 4. The envelope approach was also implemented on the same data sets.
For this simple 2 DOF system (10 & 15 Hz, 1% and 2% damping) corrupted with 10% noise, increasing the value of β improves the identified frequency and damping ratio estimates from the curve-fit.Although not eliminating the bias altogether, it is reduced by a significant degree.Note that if β is set too high, this can reduce the effectiveness of the exponential weighting (effectively truncating the data).The results were achieved without increasing the model order.The equivalent envelope weighted damping results were unbiased, although with a bigger scatter on the error bounds (0.26% to1.17% and 0.8% to 2.41% respectively) being found.
The WLS approach was now applied using the conventional method of increasing the model order and considering stability plots [6] in order to determine the system modes.A single data sequence corrupted with 10% noise was considered.The frequencies and damping ratios were estimated using the LSCE and WLSCE methods for model orders between 1 and 20.Results from consecutive model orders are compared on a stability plot, as shown in Fig. 5  for the LSCE approach.If the change in frequency of a particular mode between different orders is less than 1%, then a '+' is shown on the plot.An 'o' indicates that the change in frequency is less than 1% and that of the damping ratio is less than 5%.Such modes are said to have "stabilised".Estimates at each model order where there is no stability are represented by an 'x'.
Although the over-specification procedure works and stabilised results are found for the two modes in question, the interpretation is not very clear.In Fig. 6, the effect of applying exponential weighting of β = 0.3 to the identification process is shown.The latter plot is much easier to interpret with the stabilisation occurring much sooner and with less scatter on the stabilised columns.Similar findings were obtained using the envelope weighting approach.
These results show that the Weighted Least Squares methodology can be used to enhance current methods without requiring a vast change to current analysis software.The basic tools and algorithms remain the same, and all that is required is the extra weighting step.

Conclusions
A Weighted Least Squares approach for time domain modal parameter identification has been introduced.It has been shown that it is possible to obtain improved modal parameter estimates for little extra computational effort.The approach can be applied to a wide range of different time domain system identification algorithms and models.Both exponential and envelope weighting approaches were considered, with the latter being the approach of choice as no weighting constant needs to be chosen.The envelope approach can be used for all types of time domain model.Further work is progressing to evaluate the methodology on experimental data.

Fig. 1 .
Fig. 1.Typical Impulse Response with 10% Noise/Signal Ratio.havebeen used to replace the fundamental LS component, such as Total Least Squares, Instrumental Variables, Correlation Fit and Maximum Likelihood.These more advanced methods can be shown to give better estimates with less model order overspecification being required on simple simulated data sets.However, when applied to large order real structures, the advantages of using more sophisticated algorithms are less clear cut and they can be difficult to implement.Consequently, most work in the modal community has used the LS based methods.Other techniques such as ERA and Subspace methods have been developed to produce a minimum order realisation, through use of the Singular Value Decomposition (SVD), which reduces significantly the number of estimated modes that have to be considered following an analysis.However, LS is still the fundamental element of such techniques and a large amount of model order overspecification is employed in the solution to reduce the bias.The difference equation coefficients are found by writing Eq. (1) for different time instants and then solving 