基于HHT 的結構模態參數自動化識別方法和試驗驗證

何定橋，楊 軍

(清華大學土木工程安全與耐久教育部重點實驗室，北京 100084)

結構模態參數是結構的固有屬性，在地震荷載、風荷載、溫度荷載的作用下結構出現損傷時模態參數會發生變化，可以在一定程度上反映結構的健康狀態，因此，結構模態參數的識別是結構健康監測的重要內容[1-3]。

隨著MEMS 芯片、無線通訊、人工智能和云計算等技術快速發展，基于無線智能傳感器和云平臺的結構健康監測系統開始出現。這類系統解決了傳統建筑結構監測的傳感器布線難、系統可維護性差、數據時效性差、適用場景有限等諸多問題[4-5]。通過無線智能傳感器配合云平臺可以實現建筑結構的全天候不間斷監測。雖然結構模態參數識別方法已經非常成熟，種類繁多[6-7]，但與實驗室中的模態參數識別方法不同，利用結構健康監測系統對結構進行長期監測的過程中識別算法不能出現主觀的參數選擇過程。基于HHT 的模態參數識別中多個步驟中均需要研究人員進行主觀判斷與篩選[8]， HHT 的第一步EMD 會產生虛假的IMF 分量，對虛假分量的識別與剔除往往依賴研究人員的主觀判斷[9]。真實的IMF 分量也會存在模態混疊現象，即單個IMF 包含結構多階模態信息，如果直接對存在模態混疊現象的IMF 分量進行Hilbert 變換，會損失其中部分模態的有效信息。

本文基于HHT 提出了一種自動化結構模態參數識別方法，首先將K-L 散度與DNN 結合對 EMD產生的虛假IMF 分量進行了自動化識別與剔除，接著將SSA 與Butterworth 濾波器結合對產生模態混疊現象的IMF 進行了混疊模態分離，最后利用Hilbert 變換實現模態參數的識別。

1 理論基礎

HHT 算法由HUANG 等[10]在1998 年提出，是一種分析非線性系統非平穩信號的方法，由經驗模態分解EMD 和Hilbert 變換組成。

設x(t)為結構的監測數據，首先找到上包絡線fmax(t)與下包絡線fmin(t)，記m(t)為二者均值：

如果h1(t)滿足以下兩個條件，那么h1(t)為第1 個IMF：1)h1(t)極值點個數等于零點的個數或相差1；2)h1(t)上下包絡線的均值為0。若不滿足，則重復上述過程，假設經過k次重復得到的hk(t)符合以上兩個條件，則hk(t)為第1 個IMF，記作c1(t)，得到第1 個IMF 后，將其從原信號中除去，得到剩余部分：

令x1(t)為新的x(t)，重復上述步驟，以此得到c1(t),c2(t), ···,cn(t)，直到余量函數為單調函數或者小于某一閾值時停止循環，此時便將原信號分解為一組IMF 與一個余量函數。理論情況下，分解得到的每個IMF 均包含結構的一階模態信息，將IMF 進行Hilbert 變換構造復信號，再通過最小二乘法擬合復信號的幅頻函數與相頻函數即可得到結構的模態阻尼頻率(小阻尼比下近似為模態頻率)與模態阻尼比[11]。

傳統HHT 結構模態參數識別無法直接應用于結構健康監測的系統，主要包括兩個原因：1) HHT第1 步是將監測數據進行EMD 分解，得到多個IMF，部分IMF 包含結構的固有模態信息，為真實的EMD 分量，但部分IMF 不包含結構的固有模態信息，為虛假分量，對虛假分量的識別與剔除往往依賴研究人員的主觀判斷。2) 理想狀態下單個IMF 應包含結構的單階模態信息，但實際IMF可能會包含不同尺度分量，使后續的Hilbert 變換失去物理意義，因此，需要對存在模態混疊現象的IMF 進行分離。本文針對以上兩個問題對傳統HHT 結構模態參數識別方法進行了改進，提出了自動化的模態參數識別方法，使得基于HHT 的結構模態參數識別不包含人為干預的主觀過程，可以有效應用于建筑結構的全天候不間斷長期監測，整體流程可以用圖1 表示。

2 自動化模態參數識別

2.1 數值模型介紹

采用一6 層平面混凝土框架作為數值算例，模型使用SAP2000 商業有限元軟件建立，層高3 m，柱距6 m。梁均為高0.6 m，寬0.25 m 的矩形梁，柱為0.4 m×0.4 m 矩形截面柱，該模型為線彈性模型，樓面荷載3 kN/m2，混凝土采用C30，模型外觀如圖2 所示。模型前4 階頻率分別為6.27 Hz、18.38 Hz、29.27 Hz 和38.25 Hz。

2.2 EMD 虛假分量識別與剔除

經EMD 分解后得到的IMF 并不都反映結構的真實模態，由于采樣率不足以及樣條插值選取不適當會產生EMD 虛假分量，尤其是低頻分量，導致分解結果不滿足能量守恒。目前常見的EMD虛假分量識別方法包括相關系數法，然而IMF 分量間的相關系數較為接近，且判別過程較為主觀，容易造成誤判。PENG 等[12]提出通過WPT 將原信號分解為窄帶信號再進行EMD；祁泉泉等[13]提出可以采用濾波器對原信號進行濾波再進行EMD分解；黃迪山[14]提出考慮虛假模態導致的能量不守恒。以上方法需預估結構自振頻率的區間，若結構發生損傷自振頻率變化時，結構自振頻率會超出事先設定的截止頻率，不適于結構的長期監測。

2.2.1 K-L 散度的應用

K-L 散度是一種量化兩種概率分布之間差異的方式，又稱為相對熵，可以用來衡量IMF 分量與原信號的“關系遠近”。若EMD 的一個分量為真實分量，則其與原信號較為相似，兩者之間的K-L 散度較小，反之K-L 散度較大[15]。

對數值模型從底部輸入白噪聲激勵，時長6 s，PGA 為40 mm/s2，采樣頻率為500 Hz。數值模型在荷載激勵下第6 層加速度響應經過EMD 分解得到的IMF 與原信號的K-L 散度如圖3 所示。IMF1-3與原信號的K-L 散度分別為0.0274、0.0772、0.0823，均小于0.1，其余的IMF 為EMD 虛假分量，其與原信號的K-L 散度均大于0.5。

為了更好分離真實和虛假EMD 分量，真實分量的K-L 散度應該足夠小，而虛假EMD 分量的K-L 散度應該足夠大。在K-L 散度計算的核密度估計中，窗寬h的選擇對計算K-L 散度精度的影響很大，窗寬過大或過小都會導致真實和虛假IMF的K-L 散度接近，從而影響分離效果。宋娜等[16]利用遺傳算法尋找最優的h最大化分離度，但是文中預設了原信號的真實概率密度函數已知，這一點在自動化模態參數識別中難以符合要求，且K-L 散度閾值的仍然是主觀確定的。對于不同的信號，因為采樣頻率、信號質量、噪聲強度等影響，區分真實或虛假EMD 分量的閾值不同，在處理不同數據時需要對閾值進行調整，這一點在模態自動化識別中常常難以實現。

2.2.2 DNN 的應用

深度神經網絡可以有效地應用于EMD 虛假分量的識別。將IMF 在不同窗寬下的K-L 散度及相關數據作為輸入特征，將IMF 是否為虛假分量作為輸出結果對DNN 進行訓練。當獲得新的信號，對信號進行EMD 分解后，DNN 可以自動識別該信號EMD 分解得到的IMF 中虛假分量，而無需人工選擇K-L 散度的窗寬h和閾值[17]。

DNN 訓練及驗證樣本為人工模擬結構的加速度信號經過EMD 后得到的IMF。仿真信號的確定需要參考待識別結構的動力特性與傳感器特性，仿真信號的生成公式為：

該信號由n條頻率為f1,f2, ···,fn幅值為A1,A2, ···,An的正弦信號與白噪聲疊加而成。假設一條仿真信號經過EMD 后得到S個IMF。每個IMF 樣本共提取21 個特征，第k個IMF 的特征1-10 為該IMF 與原信號在不同窗寬下的K-L 散度值：

式中：N為信號的長度；特征1 為窗寬最小(N/200)時該IMF 與原信號的K-L 散度；特征10 為窗寬最大(N/20)時該IMF 與原信號的K-L 散度。第k個IMF 的特征11-20 為該IMF 與原仿真信號在K-L 散度值的歸一化排名：

第k個IMF 的特征21 為該IMF 與原信號在不同窗寬下的K-L 散度值的標準差(也即特征1-10的標準差)：

隨機生成500 條仿真信號，經過EMD 后得到4309 條IMF 信號。DNN 采用有監督式機器學習，對每一條IMF 標記是否為虛假EMD 分量。將IMF 進行傅里葉變換并提取主頻fIMF，與原信號包含的頻率成分f1,f2, ···,fn進行對比，若fIMF與原信號包含的任一頻率成分相差不超過5%，則將該IMF 標記為EMD 真實分量，否則標記為虛假分量。

經過試驗，確定DNN 網絡層數為6 層，包含1 個輸入層、4 個隱藏層、1 個輸出層，輸入層神經元數量與樣本的特征數量一致為21 個，隱藏層神經元數量根據經驗公式估計為10 個，判別是否是虛假EMD 分量屬于分類問題，因此，輸出層神經元數量為1 個。選擇Sigmoid 函數作為輸出層的激活函數，隱藏層的激活函數選擇ReLU。學習率和丟失率的設定通常需要經過多次試驗確定，模型學習率LR 設置為0.1。丟失率的設置是為了簡化網絡結構，防止模型過擬合，提高容錯率選擇為5%。損失函數選擇二元交叉熵函數，梯度下降算法選擇小批量梯度下降法。

將70%的數據劃分為訓練集，30%的數據劃分為驗證集對模型進行訓練，模型訓練用時23.89 s。如圖4 所示，模型在訓練集上損失函數由0.628 下降到0.304，準確率達87.20%；在驗證集上損失函數由0.576 下降到0.304，準確率達87.63%，在訓練集和驗證集上均表現出很高的準確度。

2.2.3 模型驗證

提取數值模型在白噪聲荷載激勵下的加速度數據作為響應信號。為模擬實際情況，在信號中部分時間片段中加入干擾噪聲。提取數值模型2、4、6、8、10 節點數據，添加5%、10%、20%、30%的噪聲干擾，共生成20 條信號。將原信號進行EMD分解得到161 條IMF，使用DNN 對IMF 識別結果如表1。

表1 EMD 虛假分量識別準確率Table 1 EMD false component identification accuracy

在5%、10%、20%、30%的噪聲干擾水平下，識別準確率達85.7%、81%、68.3%、60.5%，整體識別準確率達73.3%。在間段噪聲干擾下，DNN 仍可準確識別EMD 的虛假分量。

2.3 混疊模態分離

結構監測信號EMD 后得到的IMF 通常不僅僅包含單階模態信息，而有可能包含多階模態信息。模態混疊產生的原因主要包括[18]：

1) 輸入信號中混有間斷的高頻小幅值信號；

2) 信號中存在間段的噪聲；

3) 信號中存在頻率相近的分量。

模態混疊問題的解決方法包括集總經驗模態分解[19](EEMD)、掩膜信號法[20]、獨立分量分析方法[21](ICA)、解相關算法等。但是以上方法都具有一定不足之處，EEMD 需要預設信號的信噪比，這一點在模態自動化識別中難以實現；掩膜信號法中對不同的信號掩膜信號幅度和頻率難以確定；ICA 的幅度不確定性會限制對信號能量信息的判斷；解相關算法中的相關系數閾值需要主觀確定，不適用于模態自動化識別。

2.3.1 分離步驟

SSA 可以根據時間序列構造出軌跡矩陣，并對軌跡矩陣進行分解、重構，從而提取出代表原時間序列不同成分的信號，分離出包含多階模態信息的IMF 中的單階模態信息。再利用Butterworth阻帶濾波器將該階模態信息從原IMF 信號中過濾，從余量函數中繼續提取剩余的模態信息。直到所有的單階模態均被有效分離。

具體的分離步驟如圖5 所示：

1) 獲取結構在白噪聲荷載激勵下加速度響應，進行EMD 后去除虛假分量，得到反應結構真實模態信息的IMF，假設該IMF 為第i個EMD 分量，記為IMFi，計算IMF 自功率譜；

2) 根據自功率譜的峰值判斷IMFi包含的模態數量；

3) 若IMFi包含的模態信息多于1 階(假設包含n階)，則利用SSA 從IMFi中分離出自功率譜幅值最大的模態，記為IMFi.1；

4) 對IMFi使用Butterworth 阻帶濾波器，過濾IMFi.1，得到ri.1，利用SSA 從ri.1中分離出自功率譜幅值最大的模態，記為IMFi.2；

5) 對IMFi使用Butterworth 阻帶濾波器，過濾IMFi.2，得到ri.2；

6) 重復步驟4)～步驟5)，直至得到IMFi.1,IMFi.2, ···, IMFi.n。

2.3.2 模型驗證

對數值模型節點11 在白噪聲荷載激勵下的加速度時程進行EMD 分解，去除EMD 虛假分量，得到前3 階IMF 為EMD 的真實分量，原信號及前3 階IMF 的自相關函數如圖6 所示。

對上述自相關函數作FFT，得到它們的自功率譜如圖7 所示，原信號包含了結構的前4 階模態信息，理想狀態下IMF1～IMF4應該分別包含結構的第4 階～第1 階模態信息，實際IMF1包含了結構的第2 階、第3 階和第4 階模態信息，IMF2包含了結構的第1 階和第2 階模態信息，IMF3包含了結構的1 階模態信息。

IMF1中第4 階模態的自功率譜幅值最大。將IMF1進行SSA，IMF1經過SVD 分解后得到1200個特征值，將特征值排序后選擇前80%的特征值進行重構得到IMF1.1，IMF1.1反應了結構的第4 階模態信息，IMF1.1及其自功率譜如圖8 所示。

使用阻帶濾波器在IMF1中過濾IMF1.1的模態信息，可以得到r1.1。r1.1包含結構的第2 階和第3 階模態信息，按照相同的方法利用SSA 分離可以得到IMF1.2和IMF1.3。IMF1.2和IMF1.3的自相關函數和自功率譜如圖9和圖10 所示。

至此包含結構第2 階～第4 階模態信息的固有模態函數IMF1已經被分離為了IMF1.1、IMF1.2和IMF1.3，均只包含結構的單階模態信息。包含2 階模態信息的IMF2混疊模態分離方式與IMF1完全相同，在此不再贅述。

2.4 Hilbert 變換

最終得到的固有模態函數包括IMF3(第1 階模態)、IMF2.2(第1 階 模 態)、IMF2.1(第2 階 模 態)、IMF1.2(第2 階模態)、IMF1.3(第3 階模態)、IMF1.1(第4 階模態)。對IMF3、IMF1.2、IMF1.3、IMF1.1進行Hilbert 變換，并通過最小二乘法擬合幅頻函數與相頻函數，便可得到結構的模態頻率與模態阻尼比。模型對于頻率的識別精度較高，結構前4 階頻率識別最大誤差為1.8%，最小為0.3%。模型對高階模態阻尼比識別存在一定誤差。模態識別結果及與SAP2000 有限元模型計算結果對比如表2 所示。

表2 模態參數識別結果Table 2 Modal parameter identification results

3 試驗驗證

3.1 振動臺試驗介紹

振動臺試驗由清華大學土木系紀曉東課題組設計并在中國建筑科學研究院地震模擬振動臺上完成[22]。試驗結構為3 層鋼筋混凝土框架剪力墻模型，長邊框剪方向為2 跨，長4.7 m，短邊框架方向為1 跨，寬3 m，層高為2 m，總高6 m。模型總重47.4 t，材料為C30 混凝土與HRB400 鋼筋。模型每層均有附加配重，第1 層、第2 層、第3層頂板配重分別為8.94 t、8.94 t、9.36 t。模型通過基礎梁與振動臺連接，模型外觀如圖11 所示。

如表3 所示，加載分為7 個大組共25 個工況。第1 組為雙向白噪聲加載，第2 組～第7 組加載順序均為框剪方向的單向地震波加載、雙向白噪聲加載、框架方向的單向地震波加載、雙向白噪聲加載。第2 組～第5 組輸入地震波為JMA Kobe NS，加速度峰值分別為70 cm/s2、200 cm/s2、400 cm/s2、620 cm/s2，第6 組、第7 組輸入地震波為JR Takatori 00，加速度峰值分別為400 cm/s2、7620 cm/s2，白噪聲加載加速度峰值均為50 cm/s2。

表3 試驗工況加載次序Table 3 Test loading sequence

3.2 EMD 虛假分量識別

試驗所用加速度傳感器采樣頻率為200 Hz，根據Nyquist 定理，可以準確識別到的最高頻率為100 Hz，通過Opensees 有限元計算顯示該結構的第6 階頻率為12.82 Hz。目標是分別識別結構x和y方向的前3 階頻率，因為需要識別結構在地震荷載輸入發生損傷后的模態參數，結構損傷剛度下降后周期會延長，自振頻率會降低，模態將會在0 Hz～10 Hz 內更加密集，所以設定DNN 訓練所用的數值仿真信號包含6 階信號，最高頻率為30 Hz。每條仿真信號采樣時間為10 s，隨機生成1000 條仿真信號，經過EMD 后得到6652 條IMF 數據。

(3)相關手續、權責不清楚。在存貨模式中，物流等監管其中角色重要。但是法律上的職責界定存在模糊。具體職責確定，風險、收益等的情況。現實中存在監管企業承擔了過大的風險，但是收益過低的情況。造成一種不匹配情景。

將每條IMF 進行傅里葉變換，若IMF 的主頻與仿真信號包含的頻率成分相差5%以內，則將該IMF 標記為EMD 真實分量，否則標記為EMD 虛假分量。神經網絡特征的選取與本文第2 節相同：

1) 特征1-10 為該IMF 與原信號在不同窗寬下的K-L 散度值；

2) 特征11-20 為該IMF 與原仿真信號在K-L散度值的歸一化排名；

3) 特征21 為IMF 與原信號在不同窗寬下的K-L 散度值的標準差。

根據IMF 的特征與標簽建立訓練數據集。神經網絡超參數需要不斷調整優化，最終參數選擇如表4 所示。

表4 DNN 模型超參數選取Table 4 DNN model parameter selection

使用70%的數據作為訓練集，30%的數據作為驗證集對模型進行訓練，模型訓練批次為20 次(迭代次數20 次)，總用時27.48 s。如圖12，在訓練集上損失函數由0.631 下降到0.416，準確率達到82.42%；在驗證集上損失函數由0.592 下降到0.413，準確率達到81.96%，模型在訓練集和驗證集上均表現出很高的準確度。

提取表3 中工況1～工況7 中每個工況最后一次白噪聲激勵下結構第3 層加速度信號，信號時長取10 s，對信號進行EMD 分解后提取IMF 的主頻，使用DNN 進行分類，標記為“虛假”或“真實”模態，IMF 的主頻及DNN 識別結果如表5 所示。

3.3 混疊模態分離

對工況1-1 下結構第3 層x方向加速度響應進行EMD，可以得到6 條IMF，其中IMF4～IMF7為EMD 虛假分量， IMF1～IMF3為EMD 真實分量。由圖13(a)中原信號的自功率譜可以看出原信號包含了結構x方向的前4 階模態信息。

IMF1主要包含了結構的第2 階、第3 階和第4 階信息(圖13(b))，IMF2和IMF3包含了結構的第2 階和第1 階信息(圖13(c)、圖13(d))。

使用本文提出的方法對IMF1進行混疊模態分離。如圖14 所示，IMF1.1反應了結構的第3 階模態信息，IMF1.2反應了結構的第2 階模態信息，IMF1.3反應了結構的第4 階模態信息。

因篇幅限制，本節只給出工況1 的x方向數據混疊模態分離步驟，工況1 的y方向、工況2-7的x與y方向計算過程與工況1 的x方向數據相同，過程將不再贅述。

3.4 Hilbert 變換

對IMF3、IMF1.2、IMF1.1和IMF1.3進行Hilbert變換構造復信號，得到幅值函數與相位函數，通過最小二乘法擬合幅值函數與相位函數，如圖15～圖18 所示，便可得到結構x方向的模態頻率與模態阻尼比。

通過最小二乘法擬合幅值函數與相位函數，便可得到試驗結果x方向的模態頻率與模態阻尼比。第1 階～第4 階頻率識別結果為1.64 Hz、7.10 Hz、13.49 Hz 與16.21 Hz，第1 階～第4 階阻尼比識別結果為10.6%、5.2%、3.3%、5.8%。

為了評價模態參數識別的準確性，將自動化模態參數識別方法的結果分別與有限元模態分析結果和頻響函數法識別的模態參數進行對比。選用兩種方法作為對比的原因是：1)結構的有限元建模比較精確的前提下模態分析結果可以看作模態參數的理論值；2)頻響函數法利用了輸入激勵的數據，識別結果可以看作試驗結構真實的模態參數。

振動臺試驗可以提取結構基礎的加速度數據作為結構的輸入激勵并計算頻響函數，通過實頻曲線與虛頻曲線可以識別結構的頻率和阻尼比。圖19 顯示了頻響函數法識別的結構頻率和阻尼比。

表6 顯示了自動化模態參數識別方法與其他兩種方法得到模態參數的結果。自動化模態參數可以準確識別結構的前4 階模態，與有限元模型的結果相比，前4 階頻率識別的誤差分別為21.9%、0.1%、3.8%、12.8%；與頻響函數法對比，前4 階頻率識別的誤差分別為8.9%、0.0%、2.2%、7.4%，前4 階阻尼比識別的誤差分別為6.1%、0.9%、1.1%、4.1%，自動化模態參數識別方法在頻率和阻尼比兩方面識別均具有較高的精度。

表6 模態參數識別結果對比Table 6 Comparison of modal parameters identification results

3.5 自動化模態參數識別的應用

將經過驗證的自動化模態參數識別方法應用于試驗結構在工況1～工況7(共25 次加載)中地震激勵輸入作用下的結構，對發生損傷后的結構進行模態參數識別。圖20～圖22 中，加載次序為奇數時，輸入激勵為白噪聲激勵，白噪聲為雙向白噪聲，雙向PGA 均為50 cm/s2；加載次序為偶數時，輸入激勵為地震荷載，輸入激勵選用的地震波、PGA、加載方向均在圖20～圖22 中標出，也可見表3 中的試驗工況。

在結構輸入地震荷載后，前3 階頻率均出現下降。第1 次加載時(工況1-1)，結構完整，x方向前3 階頻率為1.64 Hz、7.10 Hz、13.49 Hz；y方向前3 階頻率為0.92 Hz、3.67 Hz、6.17 Hz。第25 次加載時(工況7-25)，結構x方向前3 階頻率為0.63 Hz、3.41 Hz、9.44 Hz，相較結構完整時分別下降62%、52%、30%；y方向前3 階頻率為0.55 Hz、1.49 Hz、3.34 Hz，相較結構完整時分別下降40%、59%、46%。結構構件在地震荷載下的損傷導致構件剛度以及結構的總體剛度發生下降，進而導致自振頻率下降。

4 結論

本文基于HHT 提出了一種自動化結構模態參數識別方法，該方法可以應用于結構健康監測系統對結構模態參數的長期監測中，有效地解決了傳統結構模態參數識別方法中包含主觀參數選擇過程的弊端。

(1) 將K-L 散度與DNN 結合對HHT 第1 步EMD產生的虛假分量進行了自動化識別與剔除；

(2) 將SSA 與Butterworth 濾波器結合對產生模態混疊現象的IMF 進行了混疊模態分離。

(3) 利用數值模型驗證了自動化結構模態參數識別的方法的有效性，將該方法應用于振動臺試驗實際監測數據，取得了較好的識別效果，方法具有一定工程價值。

本文不足之處在于結構模態參數自動化識別方法依賴于監測數據的質量，如果監測數據沒有包含結構的有效模態信息，或數據包含的噪聲強度太高，則后續的模態參數識別方法會失效。