Method for determining modal parameters

Field of the invention

The present invention is related to the field of modal analysis of mechanical structures.

Background of the invention

In the last two decades modal analysis has become a key technology in structural dynamics analysis. Experimental modal testing has become a commonly used technique for studying the behaviour of mechanical and civil structures such as cars, aircrafts, bridges, offshore platforms and industrial machinery. Experimental modal analysis (EMA) identifies a modal model from the measured forces (or any other input) applied to the test system and the measured vibration responses. However it is not always possible to measure these inputs or apply these inputs yourself, especially when the system of interest is in operation. These operating systems are in general excited by external forces such as the wind, waves, tire-road interaction or any other form of input. Therefore system identification techniques have been developped to identify the modal model from the structure under its operational conditions using output-only data. These techniques, referred to as operational modal analysis (OMA) or output-only modal analysis, take advantage of the ambient excitation due to e.g. traffic and wind. The OMA approach uses the system responses to the unknown inputs to retrieve the modal parameters (i.e. the system poles and mode shapes). An important advantage of OMA is that the system can remain in its normal operating condition. Further it is possible to combine the fields of EMA and OMA in what is called OMAX.

Popular OMA techniques are known e.g. from patent EP1250579 and from the paper 'Reference-based stochastic subspace identification for output-only modal analysis' (B. Peeters et al, Mechanical Systems and Signal Processing, 13(6):855 - 878, 1999). However, these methods are in general derived from existing experimental modal analysis techniques and therefore require some sort of knowledge of what is exciting the structure. The general assumption is that the inputs are broad spectrum white noise. The white noise assumption is in general acceptable for loads like wind and traffic (in a small frequency range). However, in the presence of rotating components this white noise assumption is severely violated. Classic examples, where rotating components occur, are wind turbines, hydro-electrical dams, rotating machinery, vehicles,... There are some remedies available to alleviate the issue of the colored (i.e. non-white) input. They in general identify the coloring by measuring additional data (for instance rotational speed) or by applying statisticals tools. These techniques then try to remove the "spoiled" data by averaging or cutting out these data points, possibly losing valuable data. Moreover it is not always possible to do this, especially when the data to-be-removed is near the system poles. Other techniques try to incorporate the harmonic (or coloring) into the used system model, be it starting from a known rotational speed or by modelling a varying harmonic. However, although these methods show some potential, they are never independent from the input's spectral content. Any deviation of the input to the made assumptions introduces modelling errors.

In order to have operational modal analysis without the influence of the input's spectral content the use of transmissibility functions are interesting, as these functions can become independent of the input's spectral content. Transmissibility functions are defined as the relation between output variables. However it was long thought that it was impossible to retrieve modal parameters from transmissibility functions.

Recent research has proven the possibility of determining modal parameters using transmissibility functions. This started with the observation that multiple instances of a scalar (a.k.a. local) transmissibility function (i.e. a function obtained by taking the ratio between two single signals corresponding to different output degrees of freedom), measured under different loading conditions, intersect at the poles of a system. This knowledge allows determining modal parameters in a way that benefits from the transmissibility function properties. The localisation of the intersections indeed brings forth the system poles and the associated mode shapes. More precisely, it becomes feasible to perform near deterministic operational modal analysis, without any prior knowledge nor assumptions considering the input spectrum and in which a method is considered deterministic when it becomes input-spectrum independent. It is this property that turned the spotlight on transmissibility based modal analysis (TOMA) techniques.

However, it is known that these scalar transmissibility functions, although still applicable, lose their input spectrum independent properties when multiple (not fully correlated) forces excite the structure. To resolve this problem multivariable transmissibility functions need to be introduced, i.e. a matrix expressing the relation between various sets (e.g. vectors) of output variables. While these functions are deterministic even with multiple uncorrelated forces, i.e. different input sources, they did not seem to intersect at the system poles. Yet it has been shown in 'An operational modal analysis approach based on parametrically identified multivariable transmissibilities' (P. Guillaume et al., Proceedings of IOMAC, Copenhagen, 2007) that nonetheless modal parameters can be deduced using these multivariable transmissibility functions, but unfortunately in a rather complicated manner and at the cost of an increased number of required loading conditions. With different loading conditions is meant a change in conditions which result in significant changes in the (relative) participations of the investigated modes. Different loading conditions can for example be a shift in the point of application of the forces or a change in the distribution of correlated forces over the structure. This situation occurs when the wind direction changes. An increased measurement time is required to obtain the increased number of loading conditions, for example a different wind direction can serve as a new loading condition. Another example of different loading conditions can be by using a signal starting from different moments in time after an impact on the structure. While it was possible for scalar transmissibilities to determine the modal parameters based on data obtained during two different wind directions, one now has to wait for a third loading condition before the modal parameters can be obtained. This also limits the possibility to use the mulitreference method described in the above-mentioned paper when there was no access to a third loading condition.

A more mathematical background of transmissibility function based techniques is now provided. As already mentioned, until recently, it did not seem possible to use transmissibility functions for modal parameter determination. Consider for instance scalar transmissibilitiy functions, obtained by taking the ratio of two response spectra (X_l(s) and X_r(s)). Assuming a single input F located at, say, the input degree of freedom (DOF) k, it is readily verified that the transmissibility function reduces to $$T_{\mathit{lr}}^k\left(s\right)=\frac{X_l\left(s\right)}{X_r\left(s\right)}=\frac{H_{\mathit{lk}}\left(s\right)F_k\left(s\right)}{H_{\mathit{rk}}\left(s\right)F_k\left(s\right)}=\frac{N_{\mathit{lk}}\left(s\right)}{N_{\mathit{rk}}\left(s\right)}$$

with N_lk and N_rk the numerator polynomials occurring in the transfer function models $H_{\mathit{lk}}=\frac{N_{\mathit{lk}}\left(s\right)}{D\left(s\right)}$ and $H_{\mathit{rk}}=\frac{N_{\mathit{rk}}\left(s\right)}{D\left(s\right)}$ Note that both fractions make use of the common denominator D(s), containing the system poles λ_m, which are global parameters. In equation (1) these denominators disappear by taking the ratio of the two response spectra. So, it is perfectly understandable why transmissibility functions were thought to be useless in the determination of system poles. After all the denominators, containing all information of the poles, were lost in the calculation of these functions.

While lost in a direct manner, poles can still be determined in an implicit manner. One does observe that the scalar transmissibility function (1) is dependent on the input location k. This dependency can be used to find the system's poles. To prove this one can make use of the modal model between input DOF k and output DOF l $$H_{\mathit{lk}}\left(s\right)={\displaystyle \sum _{m=1}^{N_m}}\frac{{\varphi }_{\mathit{lm}}{L}_{\mathit{km}}}{s-{\lambda }_m}+\frac{{\varphi }_{\mathit{lm}}^{*}{L}_{\mathit{km}}^{*}}{s-{\lambda }_m^{*}}$$

whereby N_m denotes the number of system poles λ_m. Using the modal model it is derived that the limit value of equation (1) for the Laplace variable s converging to the system's poles, λ_m is $$\underset{s\to {\lambda }_m}{\lim }{T}_{\mathit{lr}}^k\left(s\right)=\frac{{\varphi }_{\mathit{lm}}{L}_{\mathit{km}}}{{\varphi }_{\mathit{rm}}{L}_{\mathit{km}}}=\frac{{\varphi }_{\mathit{lm}}}{{\varphi }_{\mathit{rm}}}$$

and this is again independent of the (unknown) input and its input location k. This illustrates the fact that all scalar transmissibility functions become independent of the location of an applied input only in the system poles. As a consequence scalar transmissibility functions, obtained under different loading conditions, intersect at the system's poles. Hence, localization of these intersections yields the system's poles, in which the system poles contain the natural frequencies and damping ratios.

Unfortunately one of the advantages of using transmissibility functions is lost when scalar transmissibilities are used in the presence of multiple dominant not-fully correlated inputs. Equation (1) adapted for multiple (N_i) inputs interacting with the system is given by: $$T_{\mathit{lr}}^k\left(s\right)=\frac{X_l\left(s\right)}{X_r\left(s\right)}=\frac{{\displaystyle \sum _{k=1}^{N_i}{H}_{\mathit{lk}}\left(s\right)F_k\left(s\right)}}{{\displaystyle \sum _{k=1}^{N_i}{H}_{\mathit{rk}}\left(s\right)F_k\left(s\right)}}$$

This equation clearly shows that the scalar transmissibility functions are no longer completely independent of the input's spectral content (due to the presence of F_k(s)), and hence, in general, no longer deterministic.

The significance of the loss of deterministic behaviour is visualized in Fig.1a. Nonetheless the scalar transmissibilities still intersect in the system's poles and it is therefore still possible to use so called base functions, i.e. user built functions using different instances of transmissibility functions as building blocks designed to have poles in the system poles, like frequency response functions, e.g. Δ^-1T_lr(s)=(T_lr(s)-T_lr(s))^-1. Frequency domain estimators in order to obtain the system's poles, but at decreased quality and in need of a larger amount of data to allow averaging.

The only way to circumvent this spectral dependency and to regain the initial high quality, is by expanding the transmissibility function to its multivariable form. The multivariable transmissibility function can be defined as follows: $$\left\{X_l\left(s\right)\right\}={\left[T_{\mathit{lr}}\left(s\right)\right]}_k\left\{X_r\left(s\right)\right\}$$

whereby {X_r(s)} is the complex reference response vector of size N_rx1 (with N_r denoting the number of reference responses), and {X_l(s)} the complex non-reference response vector of size (N_o-N_r)x1, with N_o the number of obtained output signals. The decision on which responses to use as reference lies completely with the analyst and can influence the quality of the results significantly.

It can be shown that the matrix [T(s)] is given as $${\left[T_{\mathit{lr}}\left(s\right)\right]}_k=\left[H_{\mathit{lk}}\left(s\right)\right]{\left[H_{\mathit{rk}}\left(s\right)\right]}^{-1}$$

with [H_lk(s)] the complex transfer matrix of size (N_o-N_r)xN_i describing the dynamic behaviour between the N_i input locations and the non-reference outputs, while [H_rk(s)] is the complex transfer matrix of size N_rxN_i between the input locations and the reference outputs. In order to ensure the existence of [H_rk(s)]^-1 it is necessary that N_r = N_i.

One can clearly recognize that the multivariable transmissibility function is independent of the spectral content of the inputs, yet still dependent on its input locations k. Nevertheless, unlike scalar transmissibilities, multivariable transmissibility functions do no longer intersect in the system's poles (see for instance Fig. 1b) and, thus, a different approach is needed to deduce the system's poles. In the above-mentioned paper by Guillaume one way to resolve this is proposed, yet this method requires an increased number of loading conditions to retrieve the system's poles.

As already mentioned, when sufficient references are selected the multi-reference transmissibility functions again become deterministic, or input-spectrum independent. In the paper by Guillaume on multivariable transmissibility based modal analysis, [T(s)] is estimated in a parametric way. The obtained [T(s)] is then used to create pseudo-transmissibility functions, which are put in a base function. This base function, designed to have similar properties as Frequency Response Functions (FRFs), is fed into a parametric estimator for frequency response functions and the modal parameters are obtained which describe the inherent dynamic properties of the system or structure under investigation. While effective, one should note that this technique requires at least three different loading conditions and the selection of reference signals.

More recently a new approach has been introduced in 'Operational Modal Parameter Identification from Power Spectrum Density Transmissibility' (W.J. Yan et al., Computer-Aided Civil and Infrastructure Engineering, 27(3), pp.202-217, March 2012) that operates the use of so called power spectral density transmissibility functions in a single loading condition. Although possible, this method yields limited advantages compared to classical OMA techniques as it requires conditions similar to the white noise assumption and yields erroneous modal parameters in the presence of dominant harmonic sources of excitation. Moreover the technique presented in the paper by Yan et al. is based upon scalar transmissibility functions and can therefore not benefit from the advantages of the multivariable transmissibility functions.

Consequently, there is a need for a method for determining modal parameters of a system under test from response measurements, that is less cumbersome and more generally applicable than the prior art solutions described above.

Summary of the invention

It is an object of embodiments of the present invention to provide for a method for determining modal parameters of a system under test that is based on response measurements, little user interaction, allows for a reduced testing time and does not require any knowledge of the input's spectral content.

The above objective is accomplished by a method according to the present invention.

In a first aspect the invention relates to a method for determining modal parameters of a system under test. The method comprises the steps of

• taking a first set of values of output signals obtained under a first loading condition of the system and a second set of values of the output signals obtained under a second loading condition of the system,
• determining for each loading condition a parametric matrix model for the system, said parametric matrix model representing a space orthogonal to the space spanned by the output signals,
• determining for each loading condition parameters of the parametric matrix model,
• expressing a polynomial eigenvalue problem based on combining the parametric matrix models and their corresponding parameters for each loading condition
• obtaining the modal parameters by solving the polynomial eigenvalue problem.

The claimed method indeed allows achieving the above-mentioned aims. The proposed approach for modal analysis is matrix-based, in other words it is a multivariable solution. Basically, only two loading conditions are required for determining the modal parameters, i.e. the system poles and modal shapes. The fact that less loading conditions are required allows an increased applicability of the method in a wider range of fields. Having less loading conditions also implies the testing time is reduced, as there is no need to wait for more loading conditions. Further, less user interaction is needed. It should be noted that two (or more) loading conditions can be present in one measurement record.

The proposed method starts from system response data, e.g. obtained by measurement. The data comprises a first set of values of output signals obtained under a first loading condition and a second set of values of the same output signals obtained under a second loading condition. For each loading condition the corresponding set of values of the output signals spans a subspace. A parametric matrix model of a subspace orthogonal to that subspace of the output signals is determined and the parameters of the matrix model are estimated. The matrix models with parameters corresponding to the first and second loading condition are then combined to derive a polynomial eigenvalue problem. Solving this eigenvalue problem finally allows determining the desired modal parameters.

In a preferred embodiment the parametric matrix model has a number of rows not exceeding the difference between the number of the output signals and the number of uncorrelated inputs acting on the system during loading. The inputs be forces, pressures, moments, but also electrical quantities like voltages or currents. In general, any actuation capable of exciting a system and thus generating output signals can be employed.

In one embodiment the parametric matrix model is obtained from a matrix polynomial representation of a multivariable transmissibility function, which in essence describes the relation between output signals. The step of determining said parameters then advantageously comprises estimating coefficients of the matrix polynomials. In this embodiment the coefficients are estimated before the polynomial eigenvalue problem is derived. In such embodiment a step may be performed of selecting among the output signals a number of reference signals for the multivariable transmissibility function, said number being equal to or exceeding the number of uncorrelated inputs acting on the system during loading.

In an alternative embodiment the parametric matrix model is represented by a single matrix polynomial and the step of determining parameters comprises estimating coefficients of the single matrix polynomial. In this embodiment a direct estimation of the unknown polynomial coefficients is performed.

In a preferred embodiment the step of obtaining the modal parameters comprises discriminating system poles from non-physical poles. That step of discriminating advantageously comprises construction of a stabilization diagram to evaluate the numerical stability of the obtained modal parameters.

In another preferred embodiment in the step of determining (estimating) parameters a transient term is taken into account, said transient term reflecting a transient excitation and/or the influence of the system's initial and final conditions. This is especially advantageous when short datasets are used. In an advantageous embodiment this transient is used for reducing the effects of a leakage error, i.e. an error caused by the ill-modeling of the structure's initial and final conditions, and the influence of said transient excitation.

In an advantageous embodiment more than two sets of values of measured output signals are taken, obtained under more than two different loading conditions. This might reduce the error upon the results obtained when very similar loading conditions were used.

In one embodiment the method starts with actually measuring the output signals with vibrational sensors. If necessary, the data is processed (for example : averaged, filtered, ...) as is common practice in the fields of application of the proposed method.

The invention relates to a program, executable on a programmable device containing instructions, which when executed, perform the method as previously described. The invention also relates to a data carrier containing such program.

For purposes of summarizing the invention and the advantages achieved over the prior art, certain objects and advantages of the invention have been described herein above. Of course, it is to be understood that not necessarily all such objects or advantages may be achieved in accordance with any particular embodiment of the invention. Thus, for example, those skilled in the art will recognize that the invention may be embodied or carried out in a manner that achieves or optimizes one advantage or group of advantages as taught herein without necessarily achieving other objects or advantages as may be taught or suggested herein.

The above and other aspects of the invention will be apparent from and elucidated with reference to the embodiment(s) described hereinafter.

Brief description of the drawings

The invention will now be described further, by way of example, with reference to the accompanying drawings, wherein like reference numerals refer to like elements in the various figures, and in which:

Fig.1 illustrates the behaviour of scalar transmissibility functions (Fig.1a) and of multivariable transmissibility functions (Fig.1b). The system's poles are indicated with dotted lines.

Fig.2 illustrates the properties of [T̃(Ω)] for a (No= 2, N_r = 1) system with two different loading conditions.

Fig.3 represents the 4-DOF model used in the simulations.

Fig.4 illustrates loading conditions as simulated.

Fig.5 represents a stabilisation diagram obtained for a least-squares complex frequency (LSCF) estimator on cross power data (s: numerically stable pole, f: numerically stable in frequency, d : stable in damping, o: unstable).

Fig.6 represents a stabilisation diagram obtained for the LSCF estimation on the pseudo-inverse of the scalar transmissibility matrix (s: numerically stable pole, f: numerically stable in frequency, d : stable in damping, o: unstable).

Fig.7 represents a stabilisation diagram obtained using the multireference transmissibility approach according to an embodiment of the invention (s: numerically stable pole, f: numerically stable in frequency, d : stable in damping, o: unstable).

Fig.8 represents a stabilisation diagram obtained using the reference-free approach according to an embodiment of the invention (s: numerically stable pole, f: numerically stable in frequency, d : stable in damping, o: unstable).

Fig.9 represents loading conditions as simulated in an example.

Fig.10 represents stabilisation diagrams to show the influence of leakage.

Detailed description of illustrative embodiments

The present invention will be described with respect to particular embodiments and with reference to certain drawings but the invention is not limited thereto but only by the claims.

Furthermore, the terms first, second and the like in the description and in the claims, are used for distinguishing between similar elements and not necessarily for describing a sequence, either temporally, spatially, in ranking or in any other manner. It is to be understood that the terms so used are interchangeable under appropriate circumstances and that the embodiments of the invention described herein are capable of operation in other sequences than described or illustrated herein.

It is to be noticed that the term "comprising", used in the claims, should not be interpreted as being restricted to the means listed thereafter; it does not exclude other elements or steps. It is thus to be interpreted as specifying the presence of the stated features, integers, steps or components as referred to, but does not preclude the presence or addition of one or more other features, integers, steps or components, or groups thereof. Thus, the scope of the expression "a device comprising means A and B" should not be limited to devices consisting only of components A and B. It means that with respect to the present invention, the only relevant components of the device are A and B.

Reference throughout this specification to "one embodiment" or "an embodiment" means that a particular feature, structure or characteristic described in connection with the embodiment is included in at least one embodiment of the present invention. Thus, appearances of the phrases "in one embodiment" or "in an embodiment" in various places throughout this specification are not necessarily all referring to the same embodiment, but may. Furthermore, the particular features, structures or characteristics may be combined in any suitable manner, as would be apparent to one of ordinary skill in the art from this disclosure, in one or more embodiments.

Similarly it should be appreciated that in the description of exemplary embodiments of the invention, various features of the invention are sometimes grouped together in a single embodiment, figure, or description thereof for the purpose of streamlining the disclosure and aiding in the understanding of one or more of the various inventive aspects. This method of disclosure, however, is not to be interpreted as reflecting an intention that the claimed invention requires more features than are expressly recited in each claim. Rather, as the following claims reflect, inventive aspects lie in less than all features of a single foregoing disclosed embodiment. Thus, the claims following the detailed description are hereby expressly incorporated into this detailed description, with each claim standing on its own as a separate embodiment of this invention.

Furthermore, while some embodiments described herein include some but not other features included in other embodiments, combinations of features of different embodiments are meant to be within the scope of the invention, and form different embodiments, as would be understood by those in the art. For example, in the following claims, any of the claimed embodiments can be used in any combination.

It should be noted that the use of particular terminology when describing certain features or aspects of the invention should not be taken to imply that the terminology is being redefined herein to be restricted to include any specific characteristics of the features or aspects of the invention with which that terminology is associated.

In the description provided herein, numerous specific details are set forth. However, it is understood that embodiments of the invention may be practiced without these specific details. In other instances, well-known methods, structures and techniques have not been shown in detail in order not to obscure an understanding of this description.

The methodology proposed in this invention employs a parametric approach to directly estimate the model [T(Ω)] and obtain the system's poles directly from this parametric model. The minimal number of required loading conditions can therefore be reduced to two.

To correctly model the multivariable transmissibility functions it is preferable to use parametric models, like e.g. numerator-denominator models, such as for instance the left matrix fraction description (LMFD) : $$\left[T\left(\mathrm{\Omega }\right)\right]={\left[D_L\left(\mathrm{\Omega }\right)\right]}^{-1}\left[N_L\left(\mathrm{\Omega }\right)\right]$$

with [D_L(Ω)] the square (N₀ - N_r)-dimensional denominator matrix polynomial and [N_L(Ω)] E ${\mathbb{C}}^{\left(N_o-N_r\right)\times N_r}$ the complex numerator matrix polynomial of size (N_o-N_r)xN_r,, which can be described as follows : $$\left[D_L\left({\mathrm{\Omega }}_k\right)\right]={\displaystyle \sum _{i=0}^{n_D}}\left[A_i\right]{\mathrm{\Omega }}_k^i$$ $$\left[N_L\left({\mathrm{\Omega }}_k\right)\right]={\displaystyle \sum _{i=0}^{n_N}}\left[B_i\right]{\mathrm{\Omega }}_k^i$$

in which [A_i] and [B_i] denote the matrices to be estimated and n_D, n_N are the model orders defined by the user.

In the frequency domain one can choose for Ω_k, the base function, e.g. polynomial base function, between a continuous-time and a discrete-time model to approximate the discrete (e.g. sampled) time signal. The models are formulated as a Laplace domain (continuous-time) model if Ω_k = jω_k, and as a Z-domain (discrete-time) model for ${\mathrm{\Omega }}_k=e^{-j{\omega }_k{T}_s},$ with T_s the sampling period, but other choices are possible.

The proposed approach is not limited to the frequency domain and can also be implemented in the (discrete) time domain by using the time domain equivalences of the numerator and denominator matrix polynomials (8) and (9) : $$\left[D_L\left(q\right)\right]={\displaystyle \sum _{i=0}^{n_D}}\left[a_i\right]q^i$$ $$\left[N_L\left(q\right)\right]={\displaystyle \sum _{i=0}^{n_N}}\left[b_i\right]q^i$$

in which g (i) is the shift operator for the discrete time domain: qⁱ x(uT_s) = x((u - i)T_s) or, for the continuous time domain, the derivative: $q^{i}x\left(t\right)=\frac{d^{i}x\left(t\right)}{\mathit{dt}}\mathrm{.}$ Again [a_i] and [b_i] are the matrices to be estimated and n_D, n_N the model orders defined by the user.

The estimation process to obtain [A_i] and [B_i] or [a_i] and [b_i] can be performed in a variety of ways. Various types of input-output time or frequency domain, deterministic or stochastic estimators can be applied, as known in the art. The choice of the estimation method can however influence the performance of the proposed methodology. For instance, the use of the linear least squares cost function allows for a one step determination of the matrix coefficients and a fast working implementation of the approach of the current invention. A good model choice, as is necessary in the proposed numerator-denominator model, can influence the properties of the non-physical poles and therefore ease the distinction between physical and non-physical poles. When only noisy data is available, more advanced estimators (such as Maximum Likelihood Estimator (MLE)) are to be preferred.

Once the matrix polynomial coefficients are known, the numerator and denominator matrices are combined : $$\left\{\stackrel{}{X}\left(\omega \right)\right\}=\left\{\begin{array}{c}{X}_l\left(\omega \right)\\ X_r\left(\omega \right)\end{array}\right\},\left[\stackrel{}{T}\left(\mathrm{\Omega }\right)\right]=\left[D_L\left(\mathrm{\Omega }\right)-N_L\left(\mathrm{\Omega }\right)\right]$$

This yields the following : $$\left[\stackrel{}{T}\left(\mathrm{\Omega }\right)\right]\left\{\stackrel{}{X}\left(\omega \right)\right\}=0$$

Algebraically speaking [T̃(Ω)] is in fact a set of (row) vectors orthogonal to the measured response and therefore a subset of what is known as the nullspace of {X̃(ω)}. Yet, while the nullspace of {X̃(ω)} is stochastic (because orthogonal to the stochastic variable {X̃(ω)}), the subset provided in [T̃(Ω)] is not, due to its origin as a transmissibility function. This behaviour can be readily shown starting from the combination of (6) and (7) : $${\left[D_L\left(\mathrm{\Omega }\right)\right]}^{-1}\left[N_L\left(\mathrm{\Omega }\right)\right]=\left[H_{\mathit{lk}}\left(\mathrm{\Omega }\right)\right]{\left[H_{\mathit{rk}}\left(\mathrm{\Omega }\right)\right]}^{-1}$$ $$\left[N_L\left(\mathrm{\Omega }\right)\right]\left[H_{\mathit{rk}}\left(\mathrm{\Omega }\right)\right]=\left[D_L\left(\mathrm{\Omega }\right)\right]\left[H_{\mathit{lk}}\left(\mathrm{\Omega }\right)\right]$$

Further rearranging the terms in this equation yields $$\left[D_L\left(\mathrm{\Omega }\right)-N_L\left(\mathrm{\Omega }\right)\right]\left[\begin{array}{c}{H}_{\mathit{lk}}\left(\mathrm{\Omega }\right)\\ H_{\mathit{rk}}\left(\mathrm{\Omega }\right)\end{array}\right]=0$$

Using (13) this results in one of the most important properties of [T̃(Ω)] $$\left[\stackrel{}{T}\left(\mathrm{\Omega }\right)\right]\left[\begin{array}{c}{H}_{\mathit{lk}}\left(\mathrm{\Omega }\right)\\ H_{\mathit{rk}}\left(\mathrm{\Omega }\right)\end{array}\right]=0$$

This means that the subset of (row) vectors in [T̃(Ω)] spans in fact the (N₀-N_r)-dimensional subspace orthogonal to the columns k of the deterministic system's transfer matrix. Consequently, this means that [T̃(Ω)] is unaffected by the inputs' spectral content yet dependent of the input locations k. This is a property shared with all (multivariable) transmissibility functions thus far used in modal analysis.

To retrieve modal parameters using [T̃(Ω)] one can still rely on the properties of transmissibility functions. As set out earlier, scalar transmissibility functions depend on the loading condition's input location k for all values of Ω, except for the system's poles. Being scalar this means they always take the same value in the poles, and curves would intersect. To retrieve the modal parameters using [T̃(Ω)] a similar line of thought can be followed. In equation (17) it can be seen that in general [T̃(Ω)] depends on the input locations k. On the other hand, it is commonly known that [H(λ_m)] is rank deficient, its column space is solely spanned by the mode shape vectors {φ_m}. This means that in a pole the space spanned by the row vectors of [T̃(λ_m)] falls inside the nullspace of the mode shape vector {φ_m} which is independent of the loading conditions.

In order to understand this property of [T̃(Ω)], Fig.2 provides an illustration for a (N₀=2, N_r=1) system. There are two different loading conditions (k = 1 and k = 2). In Fig.2a h₁, h₂ (columns of [H(s)]) and consequently T̃₁ and T̃₂ are still dependent of the input location. In Fig.2b this dependency is lost and T̃₁ and T̃₂ always coincide.

Accordingly, N_i sets of signals of the N_o output variables of interest, obtained under N_i different loading conditions (typically by measurement), can be used to perform an estimation of [T̃(s)]_l for every loading condition l. All these transmissibility functions will be different, describing different nullspaces. In the poles all these functions are orthogonal to the corresponding mode shape vector, which can be formulated as follows: $${\left[\stackrel{}{T}\left({\lambda }_m\right)\right]}_l\left\{{\varphi }_m\right\}=0\forall m=1,2,\dots ,N_m\mathit{and}\forall l=1,2,\dots ,N_l$$

So, when combining all N_i different [T̃(s)]_k into one complex large matrix of size N_l(N_o-N_r)xN_o, i.e. $\left[\mathrm{\Upsilon }\left(\mathrm{\Omega }\right)\right]\in {\mathbb{C}}^{N_l\left(N_o-N_r\right)\times N_o}$ $$\left[\mathrm{\Upsilon }\left({\lambda }_m\right)\right]\left\{{\varphi }_m\right\}=\left[\begin{array}{c}\stackrel{}{T}{\left({\lambda }_m\right)}_1\\ \stackrel{}{T}{\left({\lambda }_m\right)}_2\\ ⋮\\ \stackrel{}{T}{\left({\lambda }_m\right)}_{N_l}\end{array}\right]\left\{{\varphi }_m\right\}=\left({\displaystyle \sum _{i=0}^{n_{\mathit{max}}}}\left[\begin{array}{c}{A}_{i1},-B_{i1}\\ A_{i2},-B_{i2}\\ ⋮\\ A_{{\mathit{iN}}_l,}-B_{{\mathit{iN}}_l}\end{array}\right]{\lambda }_m^i\right)\left\{{\varphi }_m\right\}=0$$

One so obtains an expression close to a polynomial eigenvalue problem, in which N_max = max{n_D,n_N) and A_i = 0 for i > n_D and B_i = 0 for i > n_N. However, in this family of problems one has to start from square matrix coefficients. This means the matrix coefficients (20) $$\left[\begin{array}{c}{A}_{i1},-B_{i1}\\ A_{i2},-B_{i2}\\ ⋮\\ A_{{\mathit{iN}}_l,}-B_{{\mathit{iN}}_l}\end{array}\right]$$

should be transformed into square matrices before a solution can be found.

While these matrices might be square (when N_l(N_o - N_r) = N_o), they in general are not. As the skilled person will readily appreciate, multiple methods exists to square the matrix coefficients, for instance by rejecting all rows past row number N_o, or by addition of fictional columns in order to have N_l(N_o - N_r) columns, and these choices might influence the performance of the method.

Once squared the polynomial eigenvalue problem can be solved (e.g. by means of optimization techniques, the use of companion matrices or any other method known by people skilled in the art) and within the set of eigenvalues and eigenvectors, respectively, the system's poles ( λ_m) and its mode shapes {φ_m} are contained.

To assure a unique solution the rank of [Υ(Ω)] must be N₀ for all values of Ω ≠ λ_m. When [Υ(Ω)] is not of full (row) rank, there is always a vector orthogonal to [Υ(Ω)] even when Ω ≠ λ_m; in this case the nullspace of [Υ(Ω)] indeed contains more than just the mode shapes. These vectors are also solutions to the eigenvalue problem posed in equation (19) and therefore will incorrectly be interpreted as if they were mode shapes.

In order to avoid this, one must have a complete set of loading conditions that assures the full (column) rank of [Υ(Ω)]. When the set is incomplete, the technique produces false results, yet this error can easily be prevented by detecting ill-conditioning of [Υ(Ω)] and rejecting all results when [Υ(Ω)] is rank insufficient.

An example is provided to illustrate the completeness of a set. Consider following two loading conditions: $${\left\{\stackrel{}{X}\right\}}_1={\stackrel{ʹ}{h}}_1\left(\omega \right)F_{1,1}\left(\omega \right)+{\stackrel{ʹ}{h}}_2\left(\omega \right)F_{2,1}\left(\omega \right)$$ $${\left\{\stackrel{}{X}\right\}}_2={\stackrel{ʹ}{h}}_1\left(\omega \right)F_{1,2}\left(\omega \right)+{\stackrel{ʹ}{h}}_3\left(\omega \right)F_{3,2}\left(\omega \right)$$

In both (21) and (22) a force is applied on DOF 1. A second force shifts from DOF 2 towards DOF 3, so there are two different loading conditions. However, as a set of loading conditions it is incomplete. The reason for this can be found in the recurrence of a force at DOF 1. Indeed, [T̃(Ω)] is orthogonal to h₁ and h₂, while [T̃(Ω)] is orthogonal to h₁ and h₃. Since both matrices are orthogonal to h₁, their combination remains orthogonal to h₁. Therefore the condition that the nullspace of [Υ(Ω)] consist solely out of the mode shapes is not met. To complete the set of loading conditions a new loading condition should be added. For instance : $${\left\{\stackrel{}{X}\right\}}_3={\stackrel{ʹ}{h}}_1\left(\omega \right)F_{2,3}\left(\omega \right)+{\stackrel{ʹ}{h}}_3\left(\omega \right)F_{3,3}\left(\omega \right)$$

Addition of this loading condition to the already existing set would resolve the issue. [T̃(Ω)]₃ is no longer orthogonal to h₁ and the new [Υ(Ω)] is of full (column) rank. Obviously, the validation of the full (column) rank of [Υ(Ω)] is essential in the demand for a reliable estimator.

As mentioned before the polynomial eigenvalue problem results in N_l(N_o - N_r) × n_max eigenvalues and eigenvectors, a number that may exceed the number of physical poles, N_m. To detect and reject these additional mathematical poles one can for example look at criteria such as excessive damping ratios (e.g. ξ > 10%), unstable poles (ξ ≤ 0) and pure mathematical poles, or use the properties of the stabilization chart to discriminate between numerical and mathematical poles. Further alternatives can be envisaged as well.

Once these operations are performed and a complete set of loading conditions is used, one has effectively obtained the complete set of system poles and mode shapes in a (near) deterministic manner, while only operational data was used.

The various steps of the method for determining modal parameters according to one embodiment of the invention are now described :

• obtain the system response as in N_o output signals, in which N_o exceeds the number of uncorrelated inputs exciting the structure, for instance by using a set of N_o vibrational sensors, of the system excited by an unknown input
• (optional) process the data in such a manner that the influence of measurement noise and leakage can be reduced, using commonly known methodologies
• select the number of references (N_r) in such a manner that the multi-reference transmissibility functions become (near) deterministic
• select the reference signals,
• repeat the previous steps for N_l loading conditions, until a complete set of loading conditions is obtained.
• estimate for every loading condition the matrix coefficients [A_i],[B_i], for the user-defined model orders n_D, n_N, in the preferred domain, and construct the matrix [T̃(Ω)]_k for the data linked to one loading condition.
• construct the combined matrix [Υ(Ω)] and solve the polynomial eigenvalue problem : $$\left({\displaystyle \sum _{i=0}^{n_{\mathit{max}}}}\left[\begin{array}{c}{A}_{i1},-B_{i1}\\ A_{i2},-B_{i2}\\ ⋮\\ A_{{\mathit{iN}}_l,}-B_{{\mathit{iN}}_l}\end{array}\right]{\lambda }_m^i\right)\left\{{\varphi }_m\right\}=0$$
• the system poles and mode shapes are now found within the set of solutions of the polynomial eigenvalue problem (24)
• to discriminate between physical and non-physical poles within the set of solutions one can use commonly used criteria such as for example : excessive damping ratios (e.g. ξ > 10%), unstable poles (ξ ≤ 0) and pure mathematical poles.

In addition there also exists the possibility to construct what is called a stabilization diagram to further distinguish between physical and non-physical poles. This can be done by repeating the steps of matrix coefficient estimation, matrix construction, eigenvalue problem solving and discrimination between found poles for different model orders n_D and n_N and evaluating the numerical stability of the obtained results over these different model orders, as is common practice within the modal analysis community.

The method as described above yields promising results. However, it suffers from one major drawback common to most existing techniques. There exists no clear criterium of what is a good choice for the reference signals. To alleviate the user from this decision a reference-less approach could not only save the user a significant amount of time, but also allow for an easier automatization of the technique.

Whereas above the multivariable transmissibility function was estimated as defined by formula (5), it is also possible to directly estimate the coefficients of [T̃(Ω)] as described in equation (13). This can be done by using a single matrix polynomial. An example is given for the frequency domain : $$\left[\stackrel{}{T}\left({\mathrm{\Omega }}_k\right)\right]={\displaystyle \sum _{i=0}^{n_C}}\left[C_i\right]{\mathrm{\Omega }}_k^i$$

wherein the [C_i], the complex matrix coefficients of size N_rows x N_o, are the parameters to be estimated, n_C the used model order and N_rows the preferred number of rows of matrix coefficiënts. Similar to the previous, N_rows should be smaller than or equal to N_o minus the number of uncorrelated forces acting on the system of interest during this loading condition. To estimate the coefficients [C_i] one can again rely on a variety of common estimators in the Laplace, Z- , continuous time, discrete time or any other preferred domain. Again the properties of these estimation techniques differ and may result in different performance properties of the proposed method.

While at a first glance there is only a minor difference between this adapted form and the previous, the impact of this change is rather fundamental. By estimating the nullspace directly rather than constructing it from the parametric model of the mulitreference transmissibility function, according to (12), there no longer is a need to define the references. Indeed, a closer look at (13) reveals that there is no longer a discrimination between {X_r(ω)} and {X_l(ω)} and the complete response vector can be used as a whole, relieving the user of the need to define references.

To retrieve the system poles, the method as previously described can be reused with the minor alterations as illustrated below. The adapted algorithm can be described in the following steps :

• obtain the system response as in N_o output signals, in which N_o exceeds the number of uncorrelated inputs exciting the structure, for instance by using a set of N_o vibrational sensors, of the system excited by an unknown input
• (optional) process the data in such a manner that the influence of measurement noise and leakage can be reduced, using the commonly known methodologies.
• select the number of rows of the matrix coefficients [C_i] and estimate the matrix coefficients [C_i], for the user-defined model orders n_C, in the preferred domain,
• repeat the previous steps for N_l loading conditions, until a complete set of loading conditions is obtained
• construct the combined matrix [Υ(Ω)] and solve the polynomial eigenvalue problem : $$\left({\displaystyle \sum _{i=0}^{n_c}}\left[\begin{array}{c}{C}_{i1}\\ C_{i2}\\ ⋮\\ C_{{\mathit{iN}}_l}\end{array}\right]{\lambda }_m^i\right)\left\{{\varphi }_m\right\}=0$$
• the system poles and mode shapes are now found within the set of solutions of the polynomial eigenvalue problem
• to discriminate between physical and non-physical poles within the set of solutions one can use commonly used criteria such as for example : excessive damping ratios (e.g. ξ > 10%), unstable poles (ξ ≤ 0) and pure mathematical poles.

Again it is possible to use different model orders in a stabilization diagram as mentioned previously.

In situations wherein only short datasets are available, for instance flight flutter, road-testing, or fast-varying boundary conditions, so called leakage errors can cause an issue. In general the presented methods already have a relative good robustness against leakage errors. While preprocessing in order to reduce the leakage error is possible and common practice (for instance using the Hanning window) it would be beneficial if this preprocessing could be avoided. In this way the required amount of data can be reduced and also the error on the obtained results, since some windows introduce an error upon the estimated system poles.

For classic modal analysis there is a way to incorporate these boundary conditions by addition of a transient term to the model. This concept can be applied to extend the model given in (13) as follows : $$\left[\stackrel{}{T}\left(\mathrm{\Omega }\right)\right]\left\{X\left(\omega \right)\right\}+\tau \left(\mathrm{\Omega }\right)=0$$

wherein τ(Ω), the transient polynomial, can be incorporated within T̃(Ω) as follows : $$\left[\stackrel{}{T}\left(\mathrm{\Omega }\right),\tau \left(\mathrm{\Omega }\right)\right]\left\{\begin{array}{c}X\left(\omega \right)\\ 1\end{array}\right\}=\left[\stackrel{}{T}\left(\mathrm{\Omega }\right),\tau \left(\mathrm{\Omega }\right)\right]\left\{\stackrel{}{X}\left(\omega \right)\right\}=0$$

A polynomial model can be estimated from the extended system response vector {X̃(ω) using the following : $$\left[\stackrel{}{T}\left({\mathrm{\Omega }}_k\right),\tau \left(\mathrm{\Omega }\right)\right]={\displaystyle \sum _{i=0}^{n_{\stackrel{}{C}}}}\left[{\stackrel{}{C}}_i\right]{\mathrm{\Omega }}_k^i$$

in which the [C̃_i], the complex matrix coefficients of size N_rows x (N_o+1), are the coefficients to be estimated. Since the transient effects, i.e. the cause of a leakage error, are accounted for within the last columns of the matrix coefficients [C̃_i] one can reduce the influence of the transients by removing these columns, and hence reduce the leakage error. This allows for system identification on shorter data sets.

The adapted algorithm can then be described using the following steps :

• obtain the system response as in N_o output signals, in which N_o exceeds the number of uncorrelated inputs exciting the structure, for instance by using a set of N₀ vibrational sensors, of the system excited by an unknown input,
• (optional) process the data in such a manner that the influence of measurement noise can be reduced, using the commonly known methodologies, and extend the response vectors as described in (28)
• select the number of rows of the matrix coefficiënts [C̃_i] and
• estimate the matrix coefficients [C̃_i], for the user-defined model order n_C, in the preferred domain and repeat this for every loading condition
• remove the columns associated with the transient terms from all instances of [C̃_i], in order to obtain the "leakage free" matrix coefficients [C̃_i]
• repeat the previous steps for N_l loading conditions until a complete set of loading conditions is obtained
• construct the combined matrix [Υ(Ω)] and solve the polynomial eigenvalue problem : $$\left({\displaystyle \sum _{i=0}^{n_{\mathit{max}}}}\left[\begin{array}{c}{C}_{i1}\\ C_{i2}\\ ⋮\\ C_{{\mathit{iN}}_l}\end{array}\right]{\lambda }_m^i\right)\left\{{\varphi }_m\right\}=0$$
• the system poles and mode shapes are now found within the set of solutions of the polynomial eigenvalue problem
• to discriminate between physical and non-physical poles within the set of solutions one can use commonly used criteria such as for example : excessive damping ratios (e.g. ξ > 10%), unstable poles (ξ ≤ 0) and pure mathematical poles.

Again there exists the possibility to use different model orders in a stabilization diagram as was mentioned previously.

To illustrate the performance of the method according to this invention a series of simulations using a simple mass dampener spring model, as depicted in Fig.3 is performed. The resonant frequencies and damping ratios are shown in the following table :

No.

ƒ_res (Hz)

ξ(%)

1

3.1105

0.4886

2

5.9163

0.9594

3

8.1428

1.2792

4

9.5721

1.5038

To illustrate the potential of the different techniques these series of discrete time domain simulations consist out of the measurement of the system responses to a leakage free input signal (i.e. a multisine), with similar properties to white noise. However, to illustrate the improved applicability of the method of this invention a colouring is added to the signal between the third and fourth mode, simulating the presence of a harmonic vibration source. Finally a small amount of additive measurement noise was added to increase the realistic character of the simulated results.

As a benchmark for comparison with existing techniques a cross power based approach using the Least Squares Complex Frequency (LSCF) estimator is used. As a benchmark for comparing with the existing scalar transmissibility function based approaches the pseudo-inverse approach is used.

As shown in Fig.3, a 4-DOF system was used to simulate a linear system and its responses to a series of input spectra F(ω). All degrees of freedom are considered as outputs, as if four vibrational sensors were placed on the four masses of the model. The number of outputs (N₀) is therefore equal to 4.

The simulation is performed as follows. Over two loading conditions the forces shift input location. During the first loading condition forces are applied on degrees of freedom 2 and 4. The force spectrum at DOF 4 (F₄(ω)) is a full-bandwidth multisine, the input at DOF 2 (F₂(ω)) is contaminated with a coloring at 8.6Hz, in between the 3^rd and 4^th mode, the remaining force spectra (F₁(ω) and F₃(ω)) are considered negligible compared to the two dominant input spectra. In a second loading condition the forces change input locations and now a broadband spectrum is applied on DOF 3, and the colored input is applied at DOF 1. Again the remaining force inputs are considered negligible compared to the two dominant inputs. The input spectra for these two loading conditions are also illustrated in Fig.4.

The system response spectra can be processed into auto and cross power spectra. These power spectra are averaged over the two loading conditions, this to provide the same amount of information within the cross power approach as is available for the newly developed method. In Fig.4 the stabilization diagram is given for the OMA approach using the power spectra and applying the Least Squares Complex Frequency (LSCF) estimator upon this data.

A stabilization diagram (Fig.5) is a method for discriminating mathematical poles (numerically unstable : o,d,f) from physical poles (numerically stable : s). A user is hence looking for's' indicators, preferably a line of subsequent instances of 's' indicating the presence of the pole in a wide span of model orders. In Fig.5 it is clear where the poles are located, since there are clear lines of pole indicators linked to the peaks of the background image, this is the absolute value of a cross power. However, these lines are not completely stable and are dominated by instances of 'f'. However since in every line there is at least one 's' present a relatively good estimation is possible. However the presence of the input coloring in the vicinity of the third and fourth mode influences the damping estimation of these modes in a negative manner. This can be confirmed by looking at the table below.

No.

ƒ_res (Hz)

ξ(%)

1

3.1116

0.3423

2

5.9084

0.5138

3

8.1698

0.1410

4

9.5824

0.1150

Moreover the coloring is incorrectly interpreted as a possible pole, as there is a line of pole indicators above the peak caused by the coloring of the input spectrum.

In the past the scalar transmissibility function was often suggested to resolve the influence of the input coloring. However, since in this case there are multiple uncorrelated forces exciting the structure at any given instance, scalar transmissibility functions cannot be used without further preprocessing. This of course requires an increased amount of required data. To illustrate this, a stabilization diagram (Fig.6) was constructed using an advanced pseudo-inverse approach. It is clear that without further data the scalar transmissibility functions are of limited use to solve the problem at hand.

An illustration of the multi-reference transmissibility approach is now provided. As a first illustration the algorithm for the multi-reference transmissibility approach using a linearized Least Squares estimator, implemented in the Z-domain will be given. For example the following base function can be used: $$z\left(\omega \right)=\mathit{exp}\left(\frac{j*\left(\omega -{\omega }_1\right)}{{\omega }_{N_f}-2{\omega }_1+{\omega }_2}\right)$$

with $\lfloor {\omega }_1,{\omega }_{N_f}\rfloor$ the (pulsation) band of interest and Δω the resolution within the band of interest. Of course several other possibilities exist and this is merely one illustration. After retrieving and preprocessing the response data the obtained response spectra {X(ω)} can be used in the observation matrix: $$\left[K_d{K}_n\right]\left\{\mathrm{\Theta }\right\}=0$$

with : $$K_d=\left[\begin{array}{ccccc}{X}_l^T\left({\omega }_1\right)& z{\left({\omega }_1\right)}^{-1}{X}_l^T\left({\omega }_1\right)& z{\left({\omega }_1\right)}^{-2}{X}_l^T\left({\omega }_1\right)& \cdots & z{\left({\omega }_1\right)}^{-n_D}{X}_l^T\left({\omega }_1\right)\\ X_l^T\left({\omega }_2\right)& z{\left({\omega }_2\right)}^{-1}{X}_l^T\left({\omega }_2\right)& z{\left({\omega }_2\right)}^{-2}{X}_l^T\left({\omega }_2\right)& \cdots & z{\left({\omega }_2\right)}^{-n_D}{X}_l^T\left({\omega }_2\right)\\ ⋮& ⋮& ⋮& & ⋮\\ X_l^T\left({\omega }_{N_f}\right)& z{\left({\omega }_{N_f}\right)}^{-1}{X}_l^T\left({\omega }_{N_f}\right)& z{\left({\omega }_{N_f}\right)}^{-2}{X}_l^T\left({\omega }_{N_f}\right)& \cdots & z{\left({\omega }_{N_f}\right)}^{-n_D}{X}_l^T\left({\omega }_{N_f}\right)\end{array}\right]$$ $$K_n=-\left[\begin{array}{ccccc}{X}_r^T\left({\omega }_1\right)& z{\left({\omega }_1\right)}^{-1}{X}_r^T\left({\omega }_1\right)& z{\left({\omega }_1\right)}^{-2}{X}_r^T\left({\omega }_1\right)& \cdots & z{\left({\omega }_1\right)}^{-n_N}{X}_l^T\left({\omega }_1\right)\\ X_r^T\left({\omega }_2\right)& z{\left({\omega }_2\right)}^{-1}{X}_r^T\left({\omega }_2\right)& z{\left({\omega }_2\right)}^{-2}{X}_r^T\left({\omega }_2\right)& \cdots & z{\left({\omega }_2\right)}^{-n_N}{X}_r^T\left({\omega }_2\right)\\ ⋮& ⋮& ⋮& & ⋮\\ X_r^T\left({\omega }_{N_f}\right)& z{\left({\omega }_{N_f}\right)}^{-1}{X}_r^T\left({\omega }_{N_f}\right)& z{\left({\omega }_{N_f}\right)}^{-2}{X}_r^T\left({\omega }_{N_f}\right)& \cdots & z{\left({\omega }_{N_f}\right)}^{-n_N}{X}_r^T\left({\omega }_{N_f}\right)\end{array}\right]$$

wherein X_r(ω) are the N_r reference responses. Considering the injected input in the system (two uncorrelated input spectra) N_r was set equal to two. The third and last responses in the response vector were chosen to serve as references : $$\mathrm{\Theta }=\left\{\begin{array}{c}{A}_0^T\\ A_1^T\\ ⋮\\ A_{n_D}^T\\ B_0^T\\ B_1^T\\ ⋮\\ B_{n_N}^T\end{array}\right\}$$

In order to be able to estimate the coefficients [A_i],[B_i] it is necessary to apply a constraint on the coefficients to be estimated. In this example $A_{n_D}$ is set to be a (N_o - N_r) × (N_o - N_r) identity matrix. The matrix coefficients are then found by : $${\mathrm{\Theta }}_{\mathit{LS}}={\left[K\right]}^{+}b$$

with b the negative of the (n_D+1)^th column of [K_d, K_n] and [K] is [K_d,K_n] in which the (n_D+1)^th column has been removed and wherein the operator [K]⁺ denotes the Moore Penrose pseudoinverse of [K]. From Θ_LS the matrix coefficients are found for the first loading condition. This is repeated for the second loading condition and from the combination in [Υ(Ω)] the modal parameters can be obtained from the polynomial eigenvalue problem as described in (24). One should note that the obtained system poles are positioned within the Z-domain.

In order to discriminate between physical and non-physical poles, one can rely on the aforementioned criteria and/or on the construction of a stabilization diagram. This is done by obtaining the matrix coefficients for different model orders in this case from 3 to 32. Physical poles will manifest as stable over the different model orders while mathematical poles will in general differ for each model order. Stabilization charts are commonly used to visualize the stability over different model orders and are user friendly tools to select physical poles.

Using the same data as in the benchmark tests the stabilization chart in Fig.7 was obtained using the inventive method herein introduced. It is immediately clear that now long stable lines are present within the stabilization chart, making pole selection straightforward. However, there are also stable lines present at the borders of the frequency band of interest. This is typical for the proposed method, these lines are however easily recognized as mathematical by comparing with the background image (pseudo-inverse of the transmissibility matrix) that clearly doesn't yield peaks at the border lines. The poles obtained using this method are given in table below and confirm the small error of the resulting system poles compared to the real system poles.

No.

ƒ_res(Hz)

ξ(%)

1

3.1104

0.4853

2

5.9165

0.9258

3

8.1429

1.2789

4

9.5737

1.5364

To asses whether the new method also yields correct mode shapes, the Modal Assurance Criteria (MAC) can be used to compare the estimated mode shapes with the exact system modes. In essence, the co-linearity of the estimated mode shapes with the exact mode shapes is quantified, ranging from 1, an exact correlation, to 0, perfectly orthogonal. An ideal result would therefore be an identity matrix with ones on the main diagonal and zeros elsewhere. The results for the newly developed method are given in the following table and approximate this ideal situation very well.

To illustrate the reference-less approach a Total Least Squares approach is used based on the discrete time domain data. Using the discrete time domain data makes it unnecessary to transform the measurement data in the frequency domain by means of a Fourier Transform, typically conducted using the Fast Fourier transform (FFT) algorithm.

Once the discrete time domain response data {x(u)} ∈ ${\Re }^{N_o\times 1}$ for N_p time instances, is preprocessed to the liking of the users, it can directly be fed into the following observation matrix without the need to define any references. $$K=\left[\begin{array}{ccccc}{x}^T\left(N_c\right)& x^T\left(N_c-1\right)& x^T\left(N_c-2\right)& \cdots & x^T\left(0\right)\\ x^T\left(N_c+1\right)& x^T\left(N_c\right)& x^T\left(N_c-1\right)& \cdots & x^T\left(1\right)\\ ⋮& ⋮& ⋮& & ⋮\\ x^T\left(N_p\right)& x^T\left(N_p-1\right)& x^T\left(N_p-2\right)& \cdots & x^T\left(N_p-N_c\right)\end{array}\right]$$

According to the TLS approach for obtaining the matrix coefficients it is sufficient to calculate the singular value decomposition of the observation matrix and retaining V₂, i.e. the last 2 (= N_o-N_forces) columns of the matrix V = [V₁,V₂] and with N_forces the number of uncorrelated forces acting on the system. $$\mathit{svd}\left(K\right)=U\mathrm{\Sigma }{\left[V_1{V}_2\right]}^H$$

With $${\mathrm{\Theta }}_{\mathit{TLS}}=V_2=\left\{\begin{array}{c}{C}_0^T\\ C_1^T\\ ⋮\\ C_{N_c}^T\end{array}\right\}$$

The matrix coefficients have to be estimated for the two different loading conditions and the two sets of coefficients can then be combined into one larger polynomial expression. The solution of the polynomial eigenvalue problem results, as described above, in a set of eigenvalues and eigenvectors containing amongst others the system poles and mode shape vectors. To further discriminate between mathematical poles and real physical poles the aforementioned criteria and /or the construction of a stabilization chart can be used.

One of the major advantages of using (discrete) time domain data lies in the easy updating of the observation matrix, and possibilities of fast online preprocessing by, for instance, using different filters. The use of the (discrete) time domain is to be expected within the online monitoring of modal parameters to asses structural integrity and expected lifetime.

Using the same data as in the benchmark tests the stabilization chart in Fig.8 was obtained using the inventive reference-less method herein introduced. One can clearly see that while the stabilization diagram contains more non-stable poles the real physical poles clearly manifest as stable poles from a very low model order. These long lines of stable poles indicate high quality results which is confirmed by looking at the obtained damping and resonance frequencies :

No.

ƒ_res (Hz)

ξ(%)

1

3.1105

0.4886

2

5.9163

0.9293

3

8.1428

1.2792

4

9.5721

1.5041

The results are almost identical to the true system values. Also the obtained mode shapes correlate very well to the true mode shapes :

The leakage free input is now replaced by broad-spectrum white noise, so leakage errors are present. A comparison will be made between the multireference transmissibility approach and the leakage compensated reference less approach.

To illustrate the leakage compensated approach again a Z-domain implementation of the Total Least Squares estimator is used. After preprocessing of the response data, without a need for windowing, the Fast Fourier Transform of the system response data can be fed directly into the following observation matrix: $$K=\left[\begin{array}{ccccc}{\stackrel{}{X}}^T\left({\omega }_1\right)& z{\left({\omega }_1\right)}^{-1}{\stackrel{}{X}}^T\left({\omega }_1\right)& z{\left({\omega }_1\right)}^{-2}{\stackrel{}{X}}^T\left({\omega }_1\right)& \cdots & z{\left({\omega }_1\right)}^{-n_C}{\stackrel{}{X}}^T\left({\omega }_1\right)\\ {\stackrel{}{X}}^T\left({\omega }_2\right)& z{\left({\omega }_2\right)}^{-1}{\stackrel{}{X}}^T\left({\omega }_2\right)& z{\left({\omega }_2\right)}^{-2}{\stackrel{}{X}}^T\left({\omega }_2\right)& \cdots & z{\left({\omega }_2\right)}^{-n_C}{\stackrel{}{X}}^T\left({\omega }_2\right)\\ ⋮& ⋮& ⋮& & ⋮\\ {\stackrel{}{X}}^T\left({\omega }_{N_f}\right)& z{\left({\omega }_{N_f}\right)}^{-1}{\stackrel{}{X}}^T\left({\omega }_{N_f}\right)& z{\left({\omega }_{N_f}\right)}^{-2}{\stackrel{}{X}}^T\left({\omega }_{N_f}\right)& \cdots & z{\left({\omega }_{N_f}\right)}^{-n_C}{\stackrel{}{X}}^T\left({\omega }_{N_f}\right)\end{array}\right]$$

in which $$\stackrel{}{X}\left(\omega \right)=\left\{\begin{array}{c}X\left(\omega \right)\\ 1\end{array}\right\}$$

Again by taking the singular value decomposition of the observation matrix K and retaining the 2 (= N_o-N_forces) last columns of the matrix V , the matrix coefficients are obtained: $${\mathrm{\Theta }}_{\mathit{TLS}}=\left\{\begin{array}{c}{\stackrel{}{C}}_0^T\\ {\stackrel{}{C}}_1^T\\ ⋮\\ {\stackrel{}{C}}_{n_C}^T\end{array}\right\}$$

Since the 1 was added as the (N_o+1)=5^th element of the response vector it is necessary to remove the fifth column of the matrix coefficients [C̃_i] in order to obtain the leakage compensated matrix coefficients [C̃_i]. These steps have to be performed for the two available loading conditions, after which the two sets of matrix coefficients can be combined and the polynomial eigenvalue problem can be solved to retrieve a set of eigenvalues and eigenvectors which contains the system poles and mode shapes. To discriminate in this set between mathematical and physical poles one can again rely on the discrimination criteria and stabilization charts.

As already said before the leakage free multisine signals have been replaced by broadspectrum noise. To further illustrate the power of the methods of this invention again some amount of coloring was added around 8.6Hz, in a similar way as was injected in the benchmark tests. This results in the input spectra given in Fig. 9.

To illustrate the influence of leakage, consider the stabilization diagram in Fig.10a obtained by using the multireference approach and compare it to the stabilization chart obtained using the leakage compensated approach, Fig. 10b. While it is obvious that the stabilization diagram of the multireference approach is much clearer, there is no detection of the fourth mode, moreover as one can clearly see in the table below the errors on the estimates obtained with the multi reference approach are considerably larger then these obtained with the leakage compensated technique.

No.

ƒ_res(Hz)

ξ(%)

ƒ_res(Hz)

ξ(%)

1

3.1129

0.3855

3.1104

0.4762

2

5.9782

0.9580

5.9163

0.9259

3

8.1913

11.5372

8.1419

1.2783

4

9.6712

0.1628

9.5715

1.5407

While the invention has been illustrated and described in detail in the drawings and foregoing description, such illustration and description are to be considered illustrative or exemplary and not restrictive. The foregoing description details certain embodiments of the invention. It will be appreciated, however, that no matter how detailed the foregoing appears in text, the invention may be practiced in many ways. The invention is not limited to the disclosed embodiments.

Other variations to the disclosed embodiments can be understood and effected by those skilled in the art in practicing the claimed invention, from a study of the drawings, the disclosure and the appended claims. In the claims, the word "comprising" does not exclude other elements or steps, and the indefinite article "a" or "an" does not exclude a plurality. A single processor or other unit may fulfil the functions of several items recited in the claims. The mere fact that certain measures are recited in mutually different dependent claims does not indicate that a combination of these measures cannot be used to advantage. A computer program may be stored/distributed on a suitable medium, such as an optical storage medium or a solid-state medium supplied together with or as part of other hardware, but may also be distributed in other forms, such as via the Internet or other wired or wireless telecommunication systems. Any reference signs in the claims should not be construed as limiting the scope.