多參考最小二乘復(fù)頻域法的數(shù)值問題分析及優(yōu)化

章國穩(wěn) 湯寶平 陳卓

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

關(guān)鍵詞： 模態(tài)分析; 多參考最小二乘復(fù)頻域法; 特征值分解; 奇異值分解; 穩(wěn)定圖

引 ?言

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

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

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

1 多參考最小二乘復(fù)頻域法

2 存在的問題分析

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

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

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

3 改進(jìn)多參考最小二乘復(fù)頻域法

3.1 M矩陣求解

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

4 實(shí)例分析

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

對于環(huán)境激勵(lì)響應(yīng)分析，采用系統(tǒng)正功率譜密度[9]替代頻響函數(shù)。為了充分識別所有模態(tài)，以所有測點(diǎn)為參考，即Ni=12，No=12，一共可以得到144個(gè)正功率譜密度函數(shù)，每通道頻譜線數(shù)為256。假設(shè)最高階為50，采用原始p_LSCF和本文提出的改進(jìn)p_LSCF分別進(jìn)行分析并構(gòu)造穩(wěn)定圖。其中原始p_LSCF經(jīng)過多次嘗試，本文列出三個(gè)效果較好的穩(wěn)定圖，如圖5所示?？梢钥闯鲈诓煌膮?shù)下，識別結(jié)果差距很大，其中當(dāng)TR=1-5，TM=1-25時(shí)效果較好，本文采用該結(jié)果用以對比分析。圖6為本文提出方法的穩(wěn)定圖，結(jié)果較為清晰。為了便于參考，本文也采用了COV?SSI方法[10]對數(shù)據(jù)進(jìn)行分析，所得結(jié)果對比如表3所示?？梢钥闯?，三者識別結(jié)果總體一致，在阻尼比上差異相對較大（由于數(shù)據(jù)較?。渲性紁_LSCF與改進(jìn)p_LSCF結(jié)果更接近。從實(shí)例分析可知，本文方法在不失精度的前提下可自動(dòng)得到清晰穩(wěn)定圖，更有利于參數(shù)跟蹤。在同一臺(tái)計(jì)算機(jī)基于Matlab實(shí)現(xiàn)算法分析上述數(shù)據(jù)，原始和改進(jìn)p_LSCF分別需要27.11 s和31.14 s，可見原始和改進(jìn)p_LSCF具有相似計(jì)算量。

5 結(jié) ?論

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

參考文獻(xiàn)：

[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] 石海榮，趙海峰，周國強(qiáng). PolyMAX法在自升式平臺(tái)模型損傷檢測中的應(yīng)用[J].機(jī)械設(shè)計(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.