摘要： 貝葉斯FFT算法是運營模態(tài)分析的最新一代算法，以其準確性高、計算速度快、可有效進行不確定性度量等優(yōu)點受到廣泛關(guān)注。然而，現(xiàn)有貝葉斯FFT算法針對不同情況（稀疏模態(tài)、密集模態(tài)、多步測試等）需采用不同優(yōu)化算法，且編程實現(xiàn)極為復(fù)雜。為此，本文旨在提出針對不同情況的貝葉斯FFT算法的統(tǒng)一框架，并實現(xiàn)模態(tài)參數(shù)的高效求解；視結(jié)構(gòu)模態(tài)響應(yīng)為隱變量，建立貝葉斯模態(tài)識別單步測試和多步測試的隱變量模型框架；針對提出的隱變量模型運用期望最大化算法實現(xiàn)各種情況下模態(tài)參數(shù)的統(tǒng)一貝葉斯推斷，利用隱變量解耦模態(tài)參數(shù)優(yōu)化過程，采用Louis等式間接求取似然函數(shù)的Hessian矩陣。通過兩個實際工程測試案例，并與現(xiàn)有方法對比，驗證所提方法的準確性和高效性。分析結(jié)果表明，本文所提算法與現(xiàn)有方法結(jié)果相同，但其推導(dǎo)簡單、易編程，尤其對于密集模態(tài)識別問題具有明顯的計算優(yōu)勢。本文為貝葉斯模態(tài)識別建立起統(tǒng)一的隱變量模型框架，在很大程度上簡化原本繁瑣且冗長的推導(dǎo)過程，提高計算效率，同時也為應(yīng)用變分貝葉斯、吉布斯采樣等算法求解貝葉斯模態(tài)識別問題提供了可能。

關(guān)鍵詞： 運營模態(tài)分析； 參數(shù)識別； 隱變量模型； 期望最大化； 不確定性

中圖分類號： TU311.3；U441+.2 文獻標志碼： A 文章編號： 1004-4523（2024）09-1476-09

DOI：10.16385/j.cnki.issn.1004-4523.2024.09.004

引 言

由于數(shù)據(jù)采集相對簡便且不干擾正常運營，運營模態(tài)分析已逐步發(fā)展為土木工程結(jié)構(gòu)振動測試的重要手段。該方法利用結(jié)構(gòu)在環(huán)境振動下的響應(yīng)數(shù)據(jù)識別模態(tài)參數(shù)（自振頻率、阻尼比和振型等），識別結(jié)果反映結(jié)構(gòu)在正常運營狀態(tài)下的動力學(xué)特征，為結(jié)構(gòu)健康監(jiān)測［1］、振動控制［2］和設(shè)計驗證［3］等工作的開展提供了基礎(chǔ)工具。

基于貝葉斯統(tǒng)計的運營模態(tài)識別方法可在實現(xiàn)模態(tài)參數(shù)估計的同時量化識別不確定性。它將模態(tài)參數(shù)視為隨機變量，參數(shù)識別結(jié)果以后驗概率分布形式給出，該后驗分布正比于先驗分布和似然函數(shù)的乘積。目前已開發(fā)出時域算法［4?5］、功率譜密度算法［6］和貝葉斯FFT算法［7］等模態(tài)識別方法。貝葉斯FFT算法得到了最廣泛的關(guān)注，被應(yīng)用于稀疏模態(tài)［8］、密集模態(tài)［9］、多步測試［10］和時間異步數(shù)據(jù)［11］等各種實際工況。

貝葉斯FFT算法采用Laplace逼近，即用高斯分布逼近后驗分布，將均值選作模態(tài)參數(shù)的點估計，以協(xié)方差矩陣表征識別不確定性［12］。由于在運營模態(tài)分析中數(shù)據(jù)量一般較大，根據(jù)中心極限定理，Laplace逼近可以較準確地求解模態(tài)參數(shù)的后驗概率分布。貝葉斯FFT算法通過數(shù)值優(yōu)化算法最大化似然函數(shù)，然后采用解析求導(dǎo)獲得似然函數(shù)的二階導(dǎo)數(shù)矩陣（Hessian矩陣），進而求得模態(tài)參數(shù)的協(xié)方差矩陣［13］。該方法易于理解，但其編程實現(xiàn)極為復(fù)雜，尤其是密集模態(tài)條件下計算效率大打折扣。

運營模態(tài)分析只利用結(jié)構(gòu)振動響應(yīng)數(shù)據(jù)，實驗?zāi)B(tài)分析卻可同時利用輸入、輸出信息，從而更準確地識別結(jié)構(gòu)模態(tài)參數(shù)。受此啟發(fā)，如果能夠先識別輸入荷載，進而可把運營模態(tài)分析近似轉(zhuǎn)化為實驗?zāi)B(tài)分析。需指出地是，這種條件下獲得的輸入荷載只是估計結(jié)果且具有較大不確定性，并不能直接套用實驗?zāi)B(tài)分析方法。為此，本文將運營模態(tài)分析中未被測量的輸入荷載視作“隱變量”，并應(yīng)用隱變量模型理論［14?15］考慮輸入荷載估計的不確定性，從而簡化貝葉斯FFT算法中后驗均值和Hessian矩陣的求解。實際上隱變量模型在土木工程中早有應(yīng)用：例如，KULLAA［16］將溫、濕度等當作隱變量用以消除環(huán)境因素對橋梁模態(tài)頻率的影響，NAIR等［17］使用混合高斯模型對結(jié)構(gòu)損傷識別中的特征向量進行概率建模。上述工作均從數(shù)據(jù)分析角度直接應(yīng)用不同隱變量模型，并未涉及與物理模型的結(jié)合，無法解釋實際工程問題中構(gòu)建和求解隱變量模型的具體意義。

本文首先介紹統(tǒng)計學(xué)習(xí)中隱變量模型框架及常用的期望最大化算法；然后從結(jié)構(gòu)動力學(xué)的基本物理關(guān)系出發(fā)將貝葉斯FFT算法中的概率模型歸結(jié)為隱變量模型，針對單步測試和多步測試提出求解貝葉斯FFT后驗均值和協(xié)方差矩陣的方法；最后，本文以單步測試和多步測試的兩個工程實例驗證上述方法的正確性和有效性。

1 隱變量模型概述

為解釋觀察到的現(xiàn)象，常用方法是構(gòu)建參數(shù)模型，建立輸入向量和輸出向量的函數(shù)關(guān)系，其中為參數(shù)向量。由于建模誤差和測量誤差難以避免，實際觀測向量和模型輸出之間必然存在差異，即：

（1）

式中為誤差向量。

為建立與式（1）對應(yīng)的概率分布模型，隨機嵌入方法［18］被廣泛采用，即首先通過最大熵原理［19］給定誤差的概率分布，并據(jù)此構(gòu)建似然函數(shù)，然后給定參數(shù)的先驗分布 （反映研究人員對未知參數(shù)的知識和理解），進而獲得聯(lián)合概率分布模型：

（2）

給定測量數(shù)據(jù) ，貝葉斯定理提供了計算參數(shù)后驗分布的方法：

（3）

式中為歸一化常數(shù)，常被稱為“證據(jù)函數(shù)”，廣泛應(yīng)用于模型選擇問題［20］。

由于儀器數(shù)量有限、數(shù)據(jù)缺失或無法直接獲取感興趣的物理量，并非所有輸入向量或觀測向量中的變量都能被實際測量。此時，引入隱變量，式（1）可寫成：

（4）

式中 和表示未被觀測到的隱變量，記作。應(yīng)用隨機嵌入方式，可以獲得隱變量模型：

（5）

式中 概率分布一般被稱作完全似然函數(shù)，表征了隱變量和觀測值與輸入向量的相互關(guān)系。其與似然函數(shù)的關(guān)系為：

（6）

為了區(qū)分，在隱變量模型中似然函數(shù) 又被稱作非完全似然函數(shù)。

在隱變量模型（5）中，給定測量數(shù)據(jù)，參數(shù)的后驗概率分布仍需通過式（3）求解。但由于需首先利用式（6）求得非完全似然函數(shù)，其后驗推斷過程往往更為復(fù)雜。近年來，針對隱變量模型，新的推斷算法不斷被提出，如期望最大化［21］、變分貝葉斯［22］和吉布斯采樣［23］等。這些算法極大促進了隱變量模型在實際工程中的應(yīng)用。在貝葉斯參數(shù)辨識中，期望最大化算法適用于Laplace逼近，要求大樣本和參數(shù)具備全局可辨識性；變分貝葉斯可用于解決小樣本推斷問題，且可規(guī)避Hessian矩陣的計算難題；吉布斯采樣具有廣泛的適用性，但計算效率偏低。已有經(jīng)驗表明，在大樣本條件下三種算法的計算結(jié)果極為接近。本文將貝葉斯FFT算法中的概率模型轉(zhuǎn)化為隱變量模型，通過挖掘隱變量模型的優(yōu)勢，利用期望最大化算法實現(xiàn)模態(tài)參數(shù)后驗均值和協(xié)方差矩陣的求解。

2 期望最大化算法

根據(jù)中心極限定理，在數(shù)據(jù)量較大條件下，參數(shù)的后驗概率分布可用高斯分布近似，即Laplace逼近。對于該高斯分布，其均值為極大后驗估計量，協(xié)方差矩陣由負對數(shù)后驗函數(shù)在處Hessian矩陣的逆矩陣給出，即：

（7）

（8）

式中 表示對數(shù)后驗分布函數(shù)；“”表示最大值運算；“”表示矩陣的求逆運算。期望最大化算法天然適用于Laplace逼近，本節(jié)將對該算法進行詳細介紹。

在隱變量模型中，由于非完全似然函數(shù)通常是觀測值和未知參數(shù)的復(fù)雜函數(shù)，因此求解和十分困難。在這方面，期望最大化算法提供了一種系統(tǒng)且易于處理的方法：采用迭代計算估計。設(shè)定初始值，期望最大化算法迭代有以下兩個步驟：

（1） E步：給定輸入數(shù)據(jù)，觀測數(shù)據(jù)和當前參數(shù)估計值，計算完整數(shù)據(jù)對數(shù)后驗分布函數(shù)的條件期望：

（9）

式中 “”表示針對隱變量的期望運算；

（2） M步：最大化上述，并令：

（10）

重復(fù)執(zhí)行上述兩步直至收斂。

此外，對于后驗協(xié)方差矩陣的求取，LOUIS等式［24］給出了式（8）中函數(shù)的Hessian矩陣的求解公式：

（11）

通常在期望最大化算法M步最大化函數(shù)時已經(jīng)得到函數(shù)關(guān)于的一階導(dǎo)數(shù)，因此上式第二、三項可直接利用其結(jié)果，另外函數(shù)對參數(shù)的二階求導(dǎo)一般較對參數(shù)的二階求導(dǎo)更為簡便。因此，通過式（11）計算Hessian矩陣相比于直接法往往效率更高。

期望最大化算法充分利用“分而治之”的思想，將原問題分解為更易于求解的子問題。特別是當概率模型屬于指數(shù)分布族時，E步變?yōu)榍蠼獬浞纸y(tǒng)計量的期望，在M步中函數(shù)為關(guān)于未知參數(shù)二次函數(shù)，其最優(yōu)解可解析獲得［25］。貝葉斯FFT算法的概率模型屬于指數(shù)分布族，可利用期望最大化算法實現(xiàn)高效求解。在文獻［21］中已經(jīng)應(yīng)用期望最大化算法實現(xiàn)單步測試貝葉斯FFT算法的極大后驗估計，本文將進一步應(yīng)用公式（11）完成后驗協(xié)方差估計，并拓展其在多步測試中的應(yīng)用。

3 貝葉斯FFT模態(tài)識別

運營模態(tài)分析只利用結(jié)構(gòu)的動力響應(yīng)，而缺少相應(yīng)激勵信息，因此在運用隱變量模型進行概率建模時很自然地將激勵或與之相關(guān)的物理量視作隱變量。本節(jié)將詳細討論在頻域內(nèi)構(gòu)建運營模態(tài)分析中，兩類測試的隱變量模型和期望最大化算法求解的過程。

3.1 單步測試

單步測試是指用個拾振器同時測量結(jié)構(gòu)個測點的動力響應(yīng)，且其位置不隨時間改變，是結(jié)構(gòu)健康監(jiān)測、動力測試中最常用的方法。

3.1.1 隱變量模型

假定采樣頻率為（單位：Hz），采樣時長為（單位：s），測得的結(jié)構(gòu)動力響應(yīng)記為（），則其對應(yīng)頻域數(shù)據(jù)表達為，其中為頻率處的離散傅里葉變換值，為虛數(shù)單位。應(yīng)用模態(tài)疊加法［26］，對結(jié)構(gòu)自振頻率附近一個頻帶內(nèi)的頻域數(shù)據(jù)進行概率建模：

（12）

式中 為振型矩陣（為所選頻帶內(nèi)模態(tài)階數(shù)，一般不超過3個）；為噪聲項。模態(tài)響應(yīng)表達為，其中為頻響函數(shù)矩陣，為模態(tài)激勵。在假定線性時不變系統(tǒng)和經(jīng)典阻尼的情況下，頻響函數(shù)可表示為：

（13）

式中和分別表示結(jié)構(gòu)在所選頻帶內(nèi)的第階自振頻率和阻尼比；等于0，1和2分別對應(yīng)加速度、速度和位移數(shù)據(jù)。

貝葉斯FFT方法假定模態(tài)激勵和噪聲項分別服從相互獨立的復(fù)高斯分布，即和，和為相應(yīng)的功率譜密度，表示維度為n的單位矩陣。給定測量頻域數(shù)據(jù)（；為所選頻帶內(nèi)FFT點數(shù)），為識別未知參數(shù)（，），貝葉斯FFT方法假定先驗為均勻分布，然后基于Laplace逼近原理，通過直接優(yōu)化獲得極大后驗估計，對未知參數(shù)逐一進行二階求導(dǎo)獲得Hessian矩陣并以其偽逆作為協(xié)方差矩陣［13］（下文稱直接法）。

本文應(yīng)用隱變量模型，將模態(tài)響應(yīng)視作隱變量，可得其完整數(shù)據(jù)對數(shù)后驗分布函數(shù)：

（14）

式中 。另外，基于復(fù)高斯分布的性質(zhì)，條件分布仍為復(fù)高斯分布，其均值和協(xié)方差矩陣分別為和，其中。下面在利用期望最大化算法求解時應(yīng)用這一結(jié)論。另外需要指出的是，本文曾考慮將模態(tài)激勵視作隱變量，但發(fā)現(xiàn)求解過程更為復(fù)雜，且在計算效率上也無優(yōu)勢。因此，本文最終選擇模態(tài)響應(yīng)為隱變量，從而簡化推導(dǎo)過程。

3.1.2 期望最大化算法求解

在本小節(jié)中，將針對隱變量模型式（14）應(yīng)用期望最大化算法實現(xiàn)模態(tài)參數(shù)的高效求解。

根據(jù)本文第2節(jié)介紹，在E步中求取式（14）對條件分布的期望構(gòu)建函數(shù)，并在M步中通過優(yōu)化未知參數(shù)最大化函數(shù)，最后利用Louis等式計算Hessian矩陣。受篇幅所限，且公式較為冗長，本文省略推導(dǎo)過程，主要結(jié)果如表1所示。

輸出： 模態(tài)參數(shù)后驗估計及其不確定性

相比于直接法，由于隱變量的引入，除頻率和阻尼比外，其他未知參數(shù)在M步可以以解析解的形式依次被更新，避免了直接法中原有的大量迭代步，特別是對參數(shù)量巨大的振型（正比于數(shù)據(jù)通道總數(shù)）而言，有效地提高了優(yōu)化效率。在求解Hessian矩陣時，由于是復(fù)高斯分布且其均值與協(xié)方差矩陣在E步構(gòu)建函數(shù)已求得，因此對函數(shù)求二階導(dǎo)數(shù)較直接法對似然函數(shù)求二階導(dǎo)更為簡單。另外，需要說明的是，本文提出的隱變量模型計算框架理論上適用于多階模態(tài)參數(shù)的同時識別，矩陣求導(dǎo)規(guī)則的應(yīng)用更進一步統(tǒng)一了推導(dǎo)過程，因此無需如文獻［8?9］所示明確區(qū)分稀疏模態(tài)和密集模態(tài)。

3.2 多步測試

在實際工程中，考慮到經(jīng)濟成本及現(xiàn)場測試環(huán)境，多步模態(tài)測試方式也廣泛被采用，即通過在不同測點反復(fù)使用一定數(shù)目（少量）的傳感器來分布式測量具有較大自由度數(shù)目的結(jié)構(gòu)，從而獲得結(jié)構(gòu)的全局振型。

3.2.1 隱變量模型

多步測試中每個測步的數(shù)據(jù)建模過程與單步測試基本相同，但模態(tài)參數(shù)受溫濕度等因素影響往往體現(xiàn)出一定的時變性，因此本文建模時直接考慮這一影響，假定頻率、阻尼比等模態(tài)參數(shù)在不同測步內(nèi)取不同值，但在同一測步內(nèi)保持不變，同時振型不隨時間而變以獲得全局振型。具體而言，假定某次測試共有個測步，結(jié)構(gòu)自振頻率附近一個頻帶內(nèi)的第r測步頻域數(shù)據(jù)可建模為：

（15）

式中 為全局振型矩陣；為傳感器位置矩陣（表示測步的傳感器總數(shù)），如果第個傳感器測量全局振型的第個自由度，那么的元素為1，否則為0；其余參數(shù)與式（12）具有相同含義。模態(tài)響應(yīng)包含模態(tài)激勵和頻響函數(shù)矩陣信息，頻響函數(shù)的表達式與式（13）相同，反映頻率和阻尼比的影響。

與單步測試類似，假定模態(tài)激勵和噪聲項分別服從相互獨立的復(fù)高斯分布，即和，視模態(tài)響應(yīng)為隱變量并引入均勻分布先驗，可得完整數(shù)據(jù)對數(shù)后驗分布函數(shù)，其中：

（16）

式中 表示第r測步所選頻帶內(nèi)FFT數(shù)據(jù)點數(shù)，并引入局部振型矩陣以簡化表達式。

多步模態(tài)測試未知參數(shù)集包含，總數(shù)量為。當多步測試自由度數(shù)和測步數(shù)較大時，龐大的未知參數(shù)數(shù)量給貝葉斯模態(tài)參數(shù)識別造成巨大困難。對于密集模態(tài)問題（），由于模態(tài)之間的相互干擾，已有直接法［10］更是難以收斂。本文將基于上述多步模態(tài)測試的隱變量模型，應(yīng)用期望最大化算法實現(xiàn)模態(tài)參數(shù)的極大后驗估計和協(xié)方差矩陣快速求解。

3.2.2 期望最大化算法求解

基于3.1.2節(jié)單步測試期望最大化算法，本節(jié)將提出針對多步測試的模態(tài)參數(shù)識別算法。對比公式（16）和（14）可以看出，多步模態(tài)測試的對數(shù)后驗分布函數(shù)僅比單步測試情況多一層關(guān)于測步數(shù)的求和，所以其期望最大化算法E步和M步與單步測試基本相同。值得注意的是，多步測試模型估計的振型參數(shù)是全+6A74TNsQtk3Yfsc3KZCsw==局振型，但是函數(shù)卻是關(guān)于局部振型構(gòu)建的，因此在M步獲取全局振型的解析更新形式時需要采用復(fù)合函數(shù)求導(dǎo)。另外關(guān)于后驗協(xié)方差矩陣的求取，仍然可以使用公式（11）間接求取Hessian矩陣。多步測試模態(tài)參數(shù)識別期望最大化算法的偽代碼如表2所示。

4 實例分析

本節(jié)將以工程實例驗證本文提出的模態(tài)識別算法，并通過與現(xiàn)有貝葉斯FFT算法對比驗證其準確性和高效性。

4.1 單步測試

本小節(jié)以位于浙江省海寧市的某上承式鋼筋混凝土拱橋（如圖1所示）現(xiàn)場測試為例驗證本文單步測試模態(tài)識別方法。

該橋跨度27 m，寬5 m。如圖1（b）所示，沿橋兩側(cè)等間距布置4個三向伺服式加速度計，其內(nèi)置24位采集儀，外接大容量電池為其供電，以GPS形式進行同步。測試時間為15 min，采樣頻率為100 Hz。舍棄剛度過大的沿橋方向數(shù)據(jù)，共測量16個自由度。相應(yīng)的功率譜密度（PSD）和奇異值譜（SV）如圖2所示，PSD圖中，不同顏色的線代表不同自由度的功率譜密度；SV圖中，不同顏色的線代表功率譜密度矩陣的各個特征值，并以中括號形式給出了每個模態(tài)所選頻帶范圍，奇異值譜繪制了功率譜密度矩陣的特征值隨頻率的變化情況，最大特征值的峰值代表結(jié)構(gòu)模態(tài)的可能位置，在結(jié)構(gòu)頻率附近明顯大于其他特征值的特征值數(shù)量代表了頻帶內(nèi)包含的模態(tài)數(shù)（通常不超過3個）。頻帶的選取需綜合考慮統(tǒng)計誤差和模型誤差，當所選頻帶過窄時實際采用數(shù)據(jù)有限進而造成較高的統(tǒng)計誤差，但當頻帶過寬時則易引起模型誤差，因為本文模型假定模態(tài)激勵和測量噪聲在所選頻帶內(nèi)滿足白噪聲假定。實際操作時，可以最大特征值與其他特征值的交接處作為選取頻帶的邊界。

實際結(jié)構(gòu)各個方向的激勵大小可能存在顯著差異，加之測量噪聲對低頻段影響，可能導(dǎo)致某些方向的響應(yīng)被掩蓋。如圖2中0～5 Hz的兩階模態(tài)，前一階模態(tài)的峰值明顯被后一階模態(tài)掩蓋。為了解決此問題，本文僅采用一個方向的數(shù)據(jù)識別這些模態(tài)。圖2奇異值譜中以圓圈代表各階頻率的初始值，其中黑色代表使用全部測量通道數(shù)據(jù)進行模態(tài)分析，共5個；紅色代表只用z方向（豎橋向），共2個；藍色代表只用y方向（橫橋向），共2個。

因豎橋向和橫橋向的剛度相差較大，本文所選頻帶皆為稀疏模態(tài)。由于期望最大化算法結(jié)果與直接法結(jié)果在所考慮數(shù)值精度內(nèi)完全相同，而直接法結(jié)果的正確性已被反復(fù)驗證［7］，所以本文方法的合理性也可因此得到證明。本文只列出由期望最大化算法求得的9階模態(tài)，如表3所示。

變異系數(shù)為標準差除以均值，代表了模態(tài)參數(shù)的識別不確定性。表3給出的變異系數(shù)為采用Louis等式計算所得。耗時定義為針對某階模態(tài)同一電腦運行兩種方法10次的平均時間。對于單個稀疏模態(tài)，本文方法與直接法計算效率近似，但在總計算時間上比直接法節(jié)省約5.4 s。

模態(tài)振型的識別結(jié)果如圖3所示。各階模態(tài)以模態(tài)振型的形態(tài)命名，如TY1表示y方向平動一階，R3表示扭轉(zhuǎn)三階。括號中為振型識別的變異系數(shù)，定義為振型后驗協(xié)方差矩陣跡的平方根［12］。需要指出地是，振型TZ1， TZ3， TY2和TY3因激勵過小被其他模態(tài)掩蓋而取單個方向的數(shù)據(jù)分析。結(jié)果中缺少一階和二階扭轉(zhuǎn)振動可能是因為其在測試時間內(nèi)未被環(huán)境激勵充分激發(fā)。

4.2 多步測試

本小節(jié)以浙江大學(xué)海寧國際校區(qū)鐘樓（如圖4所示）為例驗證本文多步測試模態(tài)識別方法。該建筑呈立方體形，長寬均為4 m，共計8層。使用8個三向伺服式加速度計分5個測步對其進行測試，共計120個自由度。

各測步加速度計放置在螺旋向上的樓梯處，圖5展示了各個測步測點的布置位置，圖中數(shù)字代表高度，單位為m。測步1布置在最上面兩層的四個角落（即圖中F?7和F?8的4個方塊），其余四個測步沿旋轉(zhuǎn)樓梯平臺布置。測步1的8個位置分別作為其他四個測步的參考點位置，保證其與每個測步至少有一個振型值不為零的公共測點，以便集成全局振型。

每個測步測試15 min，加上在每兩個測步之間移動和調(diào)節(jié)儀器等工作，該測試大約花費3 h。測步1采集數(shù)據(jù)的功率譜密度和奇異值譜如圖6所示，所選頻帶和初始頻率值以中括號和圓圈表示。與單步測試算例相同，黑色代表使用全部測量通道數(shù)據(jù)進行模態(tài)分析，共5個；紅色代表僅采用z方向（豎向）數(shù)據(jù)，共1個。由于結(jié)構(gòu)在水平x和y方向的剛度和質(zhì)量分布相似，所以產(chǎn)生了明顯的密集模態(tài)，體現(xiàn)在奇異值譜中同一共振頻帶內(nèi)有超過一個特征值明顯大于其他特征值，如所選的第1，4和6個頻帶，均包含兩階模態(tài)。

本文方法與直接法識別的后驗均值在所考慮精度內(nèi)完全一樣，故此處只列出本文期望最大化算法識別到的9階模態(tài)，如表4和圖7所示，各階模態(tài)以模態(tài)振型的形態(tài)命名?？梢钥闯鲱l率和阻尼比均體現(xiàn)出一定的時變特性，由于識別誤差的影響，阻尼比的變化范圍相對更大。需要指出的是，模態(tài)參數(shù)的時變特性和識別變異系數(shù)實際上表征了兩種不同的不確定性，即偶然不確定性和認知不確定性。前者屬于結(jié)構(gòu)的自身屬性，難以有效降低，例如變化的溫度可以導(dǎo)致結(jié)構(gòu)彈性模量改變，進而引起結(jié)構(gòu)頻率的時變特性；后者表征在特定模型假設(shè)和給定數(shù)據(jù)下識別量的變異性，可通過改變模型、增加數(shù)據(jù)等方式有效降低其影響。

另外本文以第1，2階密集模態(tài)和第3階稀疏模態(tài)的阻尼比識別變異系數(shù)為例進行不確定性分析，結(jié)果如表5所示。可以看出，兩種方法計算結(jié)果極為接近（相差小于1%）。對于稀疏模態(tài)，變異系數(shù)計算時間相近，均為7 s左右；但對于密集模態(tài)，本文所提方法有明顯優(yōu)勢，計算時間僅為直接法的1/3左右。

5 結(jié) 論

本文構(gòu)建了貝葉斯運營模態(tài)識別的隱變量模型框架，應(yīng)用期望最大化算法求解了模態(tài)參數(shù)的極大后驗估計和后驗協(xié)方差矩陣，為模態(tài)識別單步測試和多步測試提供了統(tǒng)一的貝葉斯推斷方法。通過理論分析和實驗驗證，可得出以下結(jié)論：

（1）隱變量模型為貝葉斯運營模態(tài)分析提供了清晰的統(tǒng)一框架，將結(jié)構(gòu)模態(tài)響應(yīng)視作隱變量，經(jīng)合理處理可極大降低模態(tài)參數(shù)識別難度；

（2）期望最大化算法為已有貝葉斯FFT模態(tài)識別算法提供了替代解法，兩者對稀疏模態(tài)的識別效果相當，但前者在密集模態(tài)條件下表現(xiàn)更優(yōu)。

本文應(yīng)用Laplace逼近原理，采用期望最大化算法求解運營模態(tài)，識別隱變量模型。除此之外，變分貝葉斯和吉布斯采樣等也常用于隱變量模型推斷，這將在以后做進一步探究。

參考文獻：

［1］茅建校，桂桂，王浩，等.大跨度斜拉橋模態(tài)參數(shù)長期追蹤及其變異性分析［J］.振動工程學(xué)報， 2024，37（1）：52-59.

MAO Jianxiao， GUI Gui， WANG Hao， et al. Long-term tracking and variability analysis of modal parameters of a long-span cable-stayed bridge［J］. Journal of Vibration Engineering， 2024，37（1）：52-59.

［2］楊超，陳政清，華旭剛，等.超長拉索多模態(tài)控制的黏滯阻尼器參數(shù)優(yōu)化研究［J］.振動工程學(xué)報，2021，34（6）：1124-1132.

YANG Chao，CHEN Zhengqing，HUA Xugang， et al. Parameters optimization of viscous dampers for multi-modal control of ultra-long stay cables［J］. Journal of Vibration Engineering， 2021，34（6）：1124-1132.

［3］陳輝，張衡，李燁君，等.測量模態(tài)不確定的梁式結(jié)構(gòu)隨機有限元模型修正［J］.振動工程學(xué)報，2019，32（4）：653-659.

CHEN Hui，ZHANG Heng，LI Yejun，et al． Stochastic finite element model updating of beam structure based on uncertain measurement modes［J］． Journal of Vibration Engineering， 2019，32（4）： 653-659．

［4］YUEN K V， KATAFYGIOTIS L S. Bayesian time-domain approach for modal updating using ambient data［J］. Probabilistic Engineering Mechanics， 2001， 16（3）： 219-231.

［5］YUEN K V， KUOK S C. Bayesian methods for updating dynamic models［J］. Applied Mechanics Reviews， 2011， 64（1）： 010802.

［6］KATAFYGIOTIS L S， YUEN K V. Bayesian spectral density approach for modal updating using ambient data［J］. Earthquake Engineering & Structural Dynamics， 2010， 30（8）：1103-1123.

［7］AU S K. Operational Modal Analysis［M］. Singapore：Springer， 2017.

［8］AU S K. Fast Bayesian FFT method for ambient modal identification with separated modes［J］. Journal of Engineering Mechanics， ASCE， 2011， 137（3）：214-226.

［9］AU S K. Fast Bayesian ambient modal identification in the frequency domain， part Ⅰ： posterior most probable value［J］. Mechanical Systems and Signal Processing， 2012， 26（1）： 60-75.

［10］ZHANG F L， AU S K， LAM H F. Assessing uncertainty in operational modal analysis incorporating multiple setups using a Bayesian approach［J］. Structural Control and Health Monitoring， 2015， 22（3）：395-416.

［11］ZHU Y C， AU S K. Spectral characteristics of asynchronous data in operational modal analysis［J］. Structural Control and Health Monitoring， 2017， 24（11）： e1981.

［12］YUEN K V， KATAFYGLOTIS L S. Bayeslan fast Fourler transform approach for modal updating using ambient data［J］. Advances in Structural Engineering， 2023， 6（2）： 81-95.

［13］AU S K. Fast Bayesian ambient modal identification in the frequency domain， part Ⅱ： posterior uncertainty［J］. Mechanical Systems and Signal Processing，2012，26（1）：76-90.

［14］EVERIT B S. An Introduction to Latent Variable Models［M］. Dordrecht： Springer， 1984.

［15］BLEI D M. Build， compute， critique， repeat： data analysis with latent variable models［J］. Annual Review of Statistics and Its Application， 2014， 1（1）：203-232.

［16］KULLAA J. Eliminating environmental or operational influences in structural health monitoring using the missing data analysis［J］. Journal of Intelligent Material Systems and Structures， 2008， 20（11）：1351-1360.

［17］NAIR K K， KIREMIDJIAN A S. Time series based structural damage detection algorithm using Gaussian mixtures modeling［J］. Journal of Dynamic Systems， Measurement， and Control， 2007， 129（3）：285-293.

［18］BECK J L. Bayesian system identification based on probability logic［J］. Structural Control and Health Monitoring， 2010， 17（7）：825-847.

［19］JAYNES E T. Probability Theory： The Logic of Science［M］. Cambridge： Cambridge University Press， 2003.

［20］DEMPSTER A P， LAIRD N M， RUBIN D B. Maximum likelihood from incomplete data via the EM algorithm［J］. Journal of the Royal Statistical Society， Series B： Methodological， 1977， 39（1）：1-38.

［21］LI B B， AU S K. An expectation-maximization algorithm for Bayesian operational modal analysis with multiple （possibly close） modes［J］. Mechanical Systems and Signal Processing， 2019， 132： 490-511.

［22］LI B B， DER KIUREGHIAN A. Operational modal identification using variational Bayes［J］. Mechanical Systems and Signal Processing， 2017， 88： 377-398.

［23］LI B B， DER KIUREGHIAN A， AU S K. A Gibbs sampling algorithm for structural modal identification under seismic excitation［J］. Earthquake Engineering & Structural Dynamics， 2018， 47（14）： 2735-2755.

［24］LOUIS T A. Finding the observed information matrix when using the EM algorithm［J］. Journal of the Royal Statistical Society， Series B： Methodological， 1982， 44（2）： 226-233.

［25］BERNARDO J M， SMITH A F. Bayesian Theory［M］. Chichester： John Wiley & Sons， 2001.

［26］CHOPRA A K. Dynamics of Structures［M］. Pearson Education India， 2007.

Latent variable model and its application to Bayesian operational modal analysis

ZHU Wei1， LI Bin-bin1，2，3， XIE Yan-long4， CHEN Xiao-yu1

（1.ZJU-UIUC Institute， Zhejiang University， Haining 314400， China； 2.Center for Balance Architecture， Zhejiang University， Hangzhou 310058， China； 3.The Architectural Design and Research Institute of Zhejiang University Co.， Ltd.， Hangzhou 310058， China； 4.College of Civil and Transportation Engineering， Shenzhen University， Shenzhen 518060， China）

Abstract： As a method for operational modal analysis （OMA）， the Bayesian FFT algorithm has garnerd significant attention for its high accuracy and efficiency， as well as its ability of uncertainty quantification. However， different cases of OMA （e.g. well-separated mode， closely-spaced modes， and multi-setup OMA） require different optimization strategy， and it is tedious in computer coding. A new framework is proposed in this paper to unify the above-mentioned cases of OMA， and the implement is simplified as a consequence. Regarding the structural modal response as a latent variable， the single-setup and multi-setup Bayesian OMA is cast as latent variable models， which have been deeply investigated in statistics. An expectation-maximization （EM） algorithm is developed for both single-setup and multi-setup OMA. The introduction of latent variables decouples the parameter optimization in EM， and Louis identity is employed to calculate the Hessian matrix. Two field tests are applied to verify the performance of the proposed approach， with a comparison to the existing algorithm. Consistent results are obtained， and a great advantage in efficiency is observed in the case of closely-spaced modes. The proposed latent variable model unifies the cases of Bayesian OMA， with the advantage of simplified implementation and fast computation. It also paves a way for a further improvement of Bayesian OMA， e.g. with the approach of variational Bayes or Gibbs sampling.

Key words： operational modal analysis；parameter identification；latent variable model；expectation maximization；uncertainty

作者簡介： 朱 偉（1998―），男，碩士。E-mail：weiz.20@intl.zju.edu.cn。

通訊作者： 謝炎龍（1988―），男，博士，副研究員。E-mail：yanlongxie@szu.edu.cn。