Skip to main content

Operational modal analysis

NoteIn this tutorial
  • You will learn how to identify the modal parameters of a structure from its responses only, when the excitation is neither controlled nor measured.
  • Prerequisites: the Experimental modal analysis tutorial, which introduces the stabilization diagram and the MAC.
  • Duration: about 20 minutes.

Large or in-service structures, such as bridges, buildings or machines in operation, cannot easily be excited by a hammer or a shaker. They are, however, permanently excited by their environment: wind, traffic, machinery, etc. Operational Modal Analysis (OMA) exploits these ambient vibrations to identify the modal parameters from the measured responses only. The price to pay is an additional assumption: the excitation must be broadband and random, so that it excites all the modes in the frequency band.

1 Step 1 – Simulating ambient vibrations

We consider the same simply supported beam as in the EMA tutorial, equipped with 19 accelerometers. The “environment” is modelled by three independent random forces, applied at points that are not known by the analyst.

using Random, Statistics

# Structure and sensors
L, b, h = 1., 0.03, 0.01
beam = Beam(L, b*h, b*h^3/12, 2.1e11, 7850.)
ξ = 0.01
x_acc = collect(0.05:0.05:0.95)

# Acquisition: 60 s at 2560 Hz → useful band up to 1 kHz
fs = 2560
T = 60.
t = 0.:1/fs:(T - 1/fs)

# Modal basis
ωn, kn = modefreq(beam, fs/2)
nmodes = length(ωn)
Φ_acc = modeshape(beam, kn, x_acc)

# Unknown ambient excitation: 3 independent white noises
x_env = [0.12, 0.43, 0.81]
Φ_env = modeshape(beam, kn, x_env)
Random.seed!(1)
F = reduce(vcat, [excitation(ColoredNoise(0., 0., T, 1.), t)' for _ in x_env])

# Measured accelerations with 30 dB SNR
prob_time = ForcedModalTimeProblem(ωn, Φ_acc, ξ, Φ_env'*F, (zeros(nmodes), zeros(nmodes)), t)
y = agwn(solve(prob_time).ddu, 30.)

# Ground truth in the analysis band
in_band = ωn .< 2π*1000
fn_exact = ωn[in_band]/2π
Φ_exact = Φ_acc[:, in_band]
sv_plot(t, y[4, :], xlabel = "Time (s)", ylabel = "Acceleration (m/s²)", title = "Accelerometer at x = 0.2 m")

Unlike an impact response, the measured signal looks like noise: the modal content is hidden in its statistics.

2 Step 2 – The OMA problem

The output-only methods do not work on FRFs but on the cross-spectral densities between the responses. They are estimated with the Welch method when the OMAProblem is built, from the time signals, the sampling rate and the block size.

bs = fs  # 1 s blocks → 1 Hz resolution
prob = OMAProblem(y, t, fs, bs, frange = [1., 1000.])
TipRule of thumb

The record must be long enough to average many blocks: OMA usually requires records of several hundred to a few thousand periods of the lowest mode of interest. Here, 60 s correspond to about 1400 periods of the first mode.

3 Step 3 – Stabilization diagram with Covariance-driven SSI

The Covariance-driven Stochastic Subspace Identification (CovSSI) fits a stochastic state-space model to the output correlations. As in EMA, the model order is unknown, and a stabilization diagram is used to separate physical and spurious poles.

stab = stabilization(prob, 20, CovSSI())
stabilization_plot(stab)

4 Step 4 – Selecting the physical modes

Let us extract the poles and the mode shapes at order 20.

poles, Ψ = modes_extraction(prob, 20, CovSSI())
fn, ξn = poles2modal(poles)

for (f, d) in zip(fn, ξn)
    println("f = ", lpad(round(f, digits = 1), 6), " Hz, ξ = ", lpad(round(100d, digits = 1), 5), " %")
end
f =   23.2 Hz, ξ =   3.9 %
f =   93.7 Hz, ξ =   1.1 %
f =  211.2 Hz, ξ =   0.9 %
f =  375.3 Hz, ξ =   1.0 %
f =  586.5 Hz, ξ =   1.0 %
f =  844.3 Hz, ξ =   1.0 %
f =  942.8 Hz, ξ =  40.9 %
f =  963.8 Hz, ξ =  39.5 %
f =  975.3 Hz, ξ =  41.5 %

Besides the six physical modes, the model contains poles with a damping ratio around 40 %. Such highly damped poles are typical noise modes: they model the measurement noise and the broadband content of the excitation rather than the structure. A simple physical criterion removes them.

physical = ξn .< 0.05 # Keep modes with less than 5 % damping
fn_id, ξn_id = fn[physical], ξn[physical]
Φ_id = real_normalization(Ψ[:, physical]) # Complex → real mode shapes
data = hcat(string.(eachindex(fn_id)), round.(fn_exact[eachindex(fn_id)], digits = 2),
            round.(fn_id, digits = 2), round.(100ξn_id, digits = 2))

pretty_table(
    HTML,
    data,
    column_labels = ["Mode", "f exact (Hz)", "f identified (Hz)", "ξ identified (%)"],
    alignment = :c
)
Mode f exact (Hz) f identified (Hz) ξ identified (%)
1 23.45 23.17 3.89
2 93.81 93.69 1.05
3 211.08 211.24 0.92
4 375.25 375.33 1.01
5 586.33 586.55 1.03
6 844.32 844.32 0.96

The natural frequencies are identified within about 1 %. The damping ratios are correct except for the first mode, which is overestimated: the damping of low-frequency modes is the most difficult quantity to estimate in OMA, because few periods are averaged and the frequency resolution is coarse compared with the half-power bandwidth of the mode (\(2\xi f_1 \approx\) 0.5 Hz here, for a resolution of 1 Hz).

5 Step 5 – Mode shapes

MAC = mac(Φ_id, Φ_exact)
println("Diagonal of the MAC matrix: ", round.(diag(MAC), digits = 3))
Diagonal of the MAC matrix: [1.0, 1.0, 1.0, 1.0, 1.0, 1.0]
fig = Figure(size = (700, 600))
for i in eachindex(fn_id)
    # OMA mode shapes are not scaled: normalize both shapes to a unit maximum for display
    ϕ_ref = Φ_exact[:, i]/maximum(abs, Φ_exact[:, i])
    ϕ_id = Φ_id[:, i]/maximum(abs, Φ_id[:, i])
    ϕ_id .*= sign(sum(ϕ_id .* ϕ_ref)) # Same sign convention

    axi = Axis(fig[(i - 1) ÷ 2 + 1, (i - 1) % 2 + 1], title = "Mode $i – $(round(fn_id[i], digits = 1)) Hz")
    lines!(axi, x_acc, ϕ_ref, label = "Exact")
    scatter!(axi, x_acc, ϕ_id, color = :orange, label = "OMA")
    xlims!(axi, 0., L)
    i == 1 && axislegend(axi, position = :cb)
end
fig

The mode shapes are perfectly identified (MAC = 1), but only up to a scaling factor.

ImportantEMA vs OMA

Because the excitation is not measured, OMA mode shapes cannot be mass-normalized. Consequently, an OMA model cannot be used to predict the response to a given force (FRF synthesis), unless the mode shapes are scaled by other means (mass-change method, finite element model,…). This is the main limitation of OMA compared with EMA.

6 Going further

  • DataSSI() works directly on the time data rather than on the correlations. It is more robust in some situations, but also much more expensive.
  • FSDD() (Frequency-Spatial Domain Decomposition) is a fast, peak-picking based method, well suited to a first analysis: modes_extraction(prob, FSDD(), pks_indices = ...).
  • The LSCF method used in EMA can also be applied to the half-spectra of an OMAProblem.

Try to reduce the record duration to T = 10. s, or the SNR to 10 dB, and observe the stabilization diagram and the damping estimates.

Reference pages: Operational modal analysis, Modal extraction.