多參考最小二乘復頻域法的數值問題分析及優化

章國穩 湯寶平 陳卓

摘要： 分析多參考最小二乘復頻域法在實際應用中的數值問題，提出一種多參考最小二乘復頻域優化方法。由于實際系統階次未知，為了避免模態遺漏，需先過估計系統階次。通過理論與數值分析發現：由于系統階次過估計，識別過程需要對奇異矩陣進行偽逆計算，偽逆計算方法和參數的選擇對識別結果有很大影響。利用特征值分解計算奇異矩陣Ro的偽逆矩陣，通過奇異值分解計算奇異矩陣Msub的偽逆矩陣，以能量為指標自動確定階次，在不失精度的前提下可自動得到清晰穩定圖。利用優化多參考最小二乘復頻域法對一個7自由度線性時不變系統和葡萄牙Infante D. Henrique大橋辨識，實驗結果表明本文方法可以在保持精度的前提下容易地識別系統模態。

關鍵詞： 模態分析; 多參考最小二乘復頻域法; 特征值分解; 奇異值分解; 穩定圖

引 ?言

20世紀90年代，針對頻域識別中數值病態問題，Guillaume等提出了最小二乘復頻域法[1];后來通過使用右矩陣分式模型代替公分母模型改進了識別模型，提出了多參考最小二乘復頻域法（本文簡稱為p_LSCF）[2]，在LMS的Test. Lab系統[3]中被稱為PolyMAX。該方法具有較為清晰的穩定圖，易于實現參數自動識別，識別精度高，被廣泛應用于各種領域[4?5]。

由于實際系統階次未知，為了避免模態遺漏，通常先過估計系統階次，接著采用穩定圖等方式從計算結果中提取系統真實模態[6]。系統階次過估計導致了計算過程中的部分矩陣為奇異矩陣，而奇異矩陣逆矩陣只能以偽逆矩陣替代。目前偽逆矩陣的計算方式有多種，不同的計算方法以及誤差容許設置對識別結果有較大影響。Ro和Msub是p_LSCF中需要求逆的中間矩陣，在系統階次過估計時為奇異矩陣，而如何確定其偽逆矩陣的最佳估計是一個值得研究的問題。

本文針對p_LSCF計算過程中奇異矩陣的逆矩陣計算問題，提出了一種改進多參考最小二乘復頻域法。利用特征值分解計算奇異矩陣Ro的偽逆矩陣，通過奇異值分解計算奇異矩陣Msub的偽逆矩陣，以能量為指標自動確定階次，可得到清晰穩定圖。對一個7自由度線性時不變系統以及葡萄牙Infante D.Henrique大橋進行模態參數識別驗證了該方法的有效性。

1 多參考最小二乘復頻域法

2 存在的問題分析

由于實際系統的階次未知，為了不遺漏系統模態，一般先對系統階次過估計，然后通過穩定圖等提取系統真實模態。由于系統階次過估計，即n大于系統實際階次，將導致算法識別過程中部分矩陣為奇異矩陣。式（8）在求解矩陣M時需要n+1維方陣Ro的逆矩陣，式（9）在求解a時需要nNi階方陣Msub的逆矩陣，當n大于系統真實階次時，Ro和Msub為奇異矩陣，無法得到它們的逆矩陣。因此，需采用它們的偽逆矩陣替代逆矩陣。目前偽逆矩陣的數值計算方法有多種，不同的計算方法和誤差容許設置對結果有較大影響。

以一7自由度線性時不變系統模態參數識別為例，如圖1所示。其中每個質量塊的質量都是1 kg，K1=10 kN/m，K2=20 kN/m，阻尼矩陣為C=0.2M+0.0003K，其中M，K分別為系統的質量矩陣、剛度矩陣。依次對各質量塊施加脈沖激勵，采集每個質量塊的速度響應，采樣頻率為500 Hz，采集時間為50 s，并且在響應信號中含有10%的測量噪聲。基于激勵響應得到49個頻響函數，譜線數為512。本文通過Matlab平臺對原p_LSCF進行實現，假設最高階為50，在計算Ro和Msub的偽逆矩陣時采用Matlab的pinv函數，其中容差Tol分別采用1-1，1-10和1-25。構造穩定圖如圖2所示，其中，特征頻率、阻尼比、模態振型的容差分別為：0.01，0.1，0.02，‘s表示穩定點，‘v表示頻率和振型穩定的點，‘d表示頻率和阻尼比穩定的點，‘f表示頻率穩定的點。

從圖2可以看出，在不同容差Tol下，本文將Ro和Msub的容差分別記為TR和TM，原p_LSCF算法得到的結果具有很大差異。當TR和TM都為1-1時，穩定圖中只出現了前2階的極點，剩下5階模態未識別出;當TR和TM都為1-10時，能較好計算出所有模態;當TR和TM都為1-15時，計算出了大量虛假模態，從中難以提取系統真實模態。算法實現過程中容差的選取對結果影響很大，但由于Matlab的偽逆算法未對外公開，對于實際系統識別過程中的容差選取也沒有理論依據，因此難以確定最佳容差參數。針對以上問題，本文對識別過程進行優化。

3 改進多參考最小二乘復頻域法

3.1 M矩陣求解

采用本文方法對第2節中的7自由度系統進行分析。計算過程中能量閾值Ethd設為97%，可以看出，在本文方法所得穩定圖中系統固有頻率形成了穩定軸，并且虛假頻率點較少，從中將非常容易獲取系統極點。表1為原始p_LSCF（TR=TM=1-10）與本文提出的p_LSCF所得的結果對比。表2給出了不同噪聲水平下，兩種方法的最大識別誤差，可以看出，兩種方法識別精度接近。該系統3和4階以及5和6階的頻率非常接近，從識別結果看出本文方法對密集模態也具有高識別精度。與原始方法相比，本文方法在保持識別精度的前提下可自動得到清晰穩定圖，避免人為嘗試設置不同參數。在同一臺計算機基于Matlab實現兩種算法分析以上數據，原始和改進p_LSCF分別需要7.08 s和7.97 s，可見原始和改進p_LSCF具有相似計算量。

4 實例分析

將本文方法應用于葡萄牙InfanteD.Henrique大橋的模態參數識別中，具體描述參考文獻[8]。橋梁在環境激勵下振動，在橋梁的4個縱向截面上布置了12個加速度傳感器，每一縱向截面上布置3個傳感器，一個垂直于側面，另2個垂直于橋面，如圖4所示。采樣頻率為12.5 Hz，對半小時的數據進行分析，即每個通道22500個數據。

對于環境激勵響應分析，采用系統正功率譜密度[9]替代頻響函數。為了充分識別所有模態，以所有測點為參考，即Ni=12，No=12，一共可以得到144個正功率譜密度函數，每通道頻譜線數為256。假設最高階為50，采用原始p_LSCF和本文提出的改進p_LSCF分別進行分析并構造穩定圖。其中原始p_LSCF經過多次嘗試，本文列出三個效果較好的穩定圖，如圖5所示。可以看出在不同的參數下，識別結果差距很大，其中當TR=1-5，TM=1-25時效果較好，本文采用該結果用以對比分析。圖6為本文提出方法的穩定圖，結果較為清晰。為了便于參考，本文也采用了COV?SSI方法[10]對數據進行分析，所得結果對比如表3所示。可以看出，三者識別結果總體一致，在阻尼比上差異相對較大（由于數據較小），其中原始p_LSCF與改進p_LSCF結果更接近。從實例分析可知，本文方法在不失精度的前提下可自動得到清晰穩定圖，更有利于參數跟蹤。在同一臺計算機基于Matlab實現算法分析上述數據，原始和改進p_LSCF分別需要27.11 s和31.14 s，可見原始和改進p_LSCF具有相似計算量。

5 結 ?論

采用多參考最小二乘復頻域法對實際結構分析時，通常需要對奇異矩陣Ro和Msub進行偽逆計算，不同偽逆計算方法和參數設置會直接影響識別結果正確性。本文主要針對Ro和Msub偽逆矩陣計算過程進行了優化：1）基于Ro是對稱矩陣的特點，采用特征值分解求解Ro偽逆矩陣，根據矩陣特征值確定矩陣的秩;2）在計算Msub偽逆矩陣時，采用奇異值分解計算Msub的偽逆矩陣，根據矩陣奇異值確定矩陣的秩。數值仿真和實例分析表明，相比于原始方法，本文方法在保持識別精度的前提下可得到更清晰穩定圖，能更加容易和可靠地識別系統模態。

參考文獻：

[1] Guillaume P， Verboven P， Vanlanduit S. Frequency-domain maximum likelihood identification of modal parameters with confidence intervals[C]. Proceeding of ISMA 23， 1998.

[2] Guillaume P， Verboven P， Vanlanduit S， et al. A poly-reference implementation of the least-squares complex frequency domain estimator[C]. Proceeding of the 21st International Modal Analysis Conference， Kissimmee， USA， 2003.

[3] LMS International. The LMS Theory and Background Book[M]. Leuven， Belgium， 2000.

[4] Lin Chang-Sheng. Frequency-domain approach for the parametric identification of structures with modal interference[J]. Journal of Mechanical Science and Technology， 2019，33 （1）： 4081-4091.

[5] 石海榮，趙海峰，周國強. PolyMAX法在自升式平臺模型損傷檢測中的應用[J].機械設計與制造工程，2016，45（3）：93-97.

Shi Hairong， Zhao Haifeng， Zhou Guoqiang. Application of PolyMAX method in the damage detection of offshore platform[J]. Machine Design and Manufacturing Engineering，2016，45（3）：93-97.

[6] Zhang Guowen， Ma Jinghua， Chen Zhou， et al. Automated eigensystem realisation algorithm for operational modal analysis[J]. Journal of Sound and Vibration， 2014，333 （15）： 3550-3563.

[7] Rolain Y， Pintelon R， Xu K Q， et al. Best conditioned parametric identification of transfer function models in the frequency domain[J]. IEEE Transactions on Automatic Control， 1995， 40（11）： 1954-1960.

[8] Magalhaes Filipe， Cunha Alvaro， Caetano Elsa. Online automatic identification of the modal parameters of a long span arch bridge[J]. Mechanical Systems and Signal Processing， 2009，23（2）： 316-329.

[9] Peeters B， Van der Auweraer H， Vanhollebeke F， et al. Operational modal analysis for estimating the dynamic properties of a stadium structure during a football game[J]. Shock and Vibration， 2007，14（4）： 283-303.

[10] Reynders E， Pintelon R， De Roeck G. Uncertainty bounds on modal parameters obtained from stochastic subspace identification[J]. Mechanical Systems and Signal Processing， 2008， 22（4）： 948-969.