辛峻峰，王树青，刘福顺

(中国海洋大学 工程学院，山东 青岛 266100)

有效的模态识别方法能精确得到结构的相关参数，从而可准确掌握大型结构的健康状况。Wang等［1-3］提出基于时域响应的模态参数识别方法。如时间序列法、随机减量法、自然激励技术、随机子空间法等。其中随机子空间法是目前较先进的环境激励下模态参数识别方法之一，能精确提取海洋平台等大型结构的模态参数。

DeMoor等［4］证明了数据驱动随机子空间算法与协方差驱动随机子空间算法在理论上的完全一致。Peters等［5］认为在应用过程中，两种方法的模态识别结果也一致。但在对导管架平台模态分析过程中，两种方法的识别结果并不一致。对此，本文研究探讨两种方法产生差异的原因，继而进行对比分析。

1 随机子空间法

简单介绍随机子空间法要点，讨论详见文献［4-6］。在仅考虑随机噪声前提下，振动系统的离散状态空间方程可表示为:

式中:xk∈Rn为离散时间状态向量;A∈Rn×n为离散状态矩阵;C∈Rl×n为输出矩阵，描述内部状态转化为外界测量值过程;wk∈Rn为由干扰、模型误差造成的过程噪声;vk∈Rl为由传感器误差等造成的测量噪声。两种噪声均不可测量，在推导过程中常被假设为零均值，平稳的白噪声。两种噪声的协方差矩阵可表示为:

式中:E为数学期望算子;Q∈Rn×n，S∈Rn×l，X∈Rl×l;δpq表示 δpq=1(p=q)，δpq=0(p≠q)，且E［wp］=0，E［vp］=0。

两种随机子空间算法的本质区别在于“过去”与“将来”相关性的计算方式不同:数据驱动随机子空间法用投影计算:

式中:i为采样时间间隔;()+表示伪逆(pseudoinverse)。而协方差驱动随机子空间法用协方差矩阵计算:

由此数据驱动与协方差驱动随机子空间法的关系主要为投影与协方差矩阵之间的关系。

2 数据驱动随机子空间法与协方差驱动随机子空间法之关系

推导［4］，投影Pi与协方差矩阵Ci为相似矩阵。据相似矩阵性质知，二者奇异值相同，从而可得相同模态参数。与文献［5］结果相同。但本文研究发现，被忽略的2个细节能使两种方法在应用中产生差异。

2.1 数据有限性

理论上式(4)、(5)成立的前提为:数据是无限的。而实际使用中由于计算时间、计算精度等因素影响，所提数据必然有限，故式(4)、(5)的右半部分(Yf/Yp与YfYTp)分别变成对投影Pi与协方差矩阵Ci估计，因数值计算方式不同，必会出现不同的估计误差，故在应用中投影Pi与协方差Ci不再相似，因而两种方法的识别结果也会不一致。即有限的数据可能是造成差异的原因之一。

2.2 QR分解

除数据影响外，研究发现在数据驱动随机子空间法的计算过程中，QR算法的使用也可能是差异产生的另一主要因素。由于常用的数据驱动随机子空间法并非直接用式(4)计算投影，而在计算投影的过程中用QR分解算法:

因此数据驱动随机子空间法应准确为:基于QR分解的数据驱动随机子空间法。本文由此推测QR算法的使用也会影响识别结果，为造成差异的另一原因。

3 数值算例

用一系列数值试验证明差异的存在，继而验证差异原因，并对两种随机子空间法进行对比。建立五自由度质量-弹簧-阻尼系统的数值模型［7］，见图1。

图1 自由度弹簧-质量-阻尼系统Fig.1 A 5-DOF mass-spring-dashpot system

每个单元质量、刚度、阻尼系数分别为:mn=50 kg，kn=2.9 ×107N/m，cn=1 000 N·s/m。位移记为xn，(n=1，…，5)。通过特征值分析，获得五阶模态频率的理论值为:34.499 Hz，100.700 Hz，158.730 Hz，203.880 Hz，232.520 Hz;五阶模态阻尼比的理论值为:0.003 737，0.010 909，0.017 197，0.022 092，0.025 198。

3.1 白噪声激励实验

采用Matlab中lsim函数［7］，以均值为零、标准差为1.006的高斯白噪声为输入，将其右向加载在靠固定端的首个质量块m1上(图1)。采样频率500 Hz，从首个质量块采集9 000个数据点，约18 s，提取的响应数据见图2。

图2 白噪声激励下首个质量块输出信号Fig.2 Response signal of the first mass with white noise loading

本文构建Hankel矩阵(i=100)，作为所有计算过程使用数据。

3.2 差异现象验证

为验证两种方法差异的存在，分别用基于QR的数据驱动随机子空间法与协方差驱动随机子空间法分析测量数据。分析结果如图3所示。

图3 两种方法稳定图(*为稳定，○为不稳定)Fig.3 Two methods of stabilization diagram

图4 两种方法稳定图(*为稳定，○为不稳定)Fig.4 Two methods of stabilization diagram

图5 两种方法稳定图(*为稳定，○为不稳定)Fig.5 Two methods of stabilization diagram

本文稳定图用信号的傅里叶变换为背景，显示模态阶数在0～30之间模态识别结果情况。对同一模态，若前后两次识别结果同一模态频率误差为1%，同时响应阻尼识别结果误差为5%，即:

则此次识别结果标识为稳定，否则为不稳定［8］。

由图3看出，两种方法的识别结果并不一致，与基于QR分解的数据驱动随机子空间法(图3(b))相比，协方差驱动随机子空间法(图3(a))特点如下:

(1)在横轴100 Hz、150 Hz左右，纵轴10～25之间形成稳定的非常接近的虚假模态(放大图);

(2)需高的阶数才能识别较弱势的第4个模态(区间在200 Hz左右);

(3)未能识别出较弱的第5个模态(区间在230 Hz左右)。

以上分析表明:差异现象是存在的。以下将展开对造成差异原因的验证。

3.3 数据有限性影响验证

本节将验证数据有限性影响数据驱动随机子空间法的识别结果。

未用QR分解的数据驱动随机子空间法与协方差驱动随机子空间方法唯一区别在于:有限数据造成的投影与协方差矩阵的估计误差。若该两种方法的分析结果不一致，则表示数据有限性造成了差异。稳定图的分析结果见图4。

由图4看出，两种方法的模态识别结果不一致:图4(a)中仅识别出第三阶模态(区间在150左右)，而图4(b)中却识别出4个模态，由此表明:数据有限性导致了估计误差的产生，是造成差异的原因之一。

3.4 QR分解算法影响验证

本节将验证QR分解算法影响了数据驱动随机子空间法的识别结果。

基于QR分解与无QR分解的数据驱动随机子空间法之间唯一区别在于:是否使用QR分解技术。若该两种方法的模态识别结果不一致，即验证QR分解算法会造成差异。

由图5看出，两种方法的识别结果差别较大:图5(a)中仅识别出第三阶模态(区间150左右)，图5(b)中识别出全部五阶模态。由此证明，QR分解提高了数据驱动随机子空间法的计算精度与模态识别能力，是造成差异的另一个原因。

3.5 两种子空间法对比分析

针对基于QR的数据驱动随机子空间法与协方差驱动随机子空间法的识别结果的差异，将采用蒙特卡洛模拟对比两种方法的优劣。

用Matlab程序中randn函数生成白噪声激励，进行100次蒙特卡洛激励试验，并统计模态阶数为10的相关模态参数估计结果(模态频率、模态阻尼)。如图6、图7所示。图中行代表4个模态;列代表协方差驱动随机子空间法(SSI/cov)与基于QR的数据驱动随机子空间法(SSI/data)。横轴为激励次数，纵轴为用真实值标准化的估计值，且所有估计值均用点表示，红线代表每个模态的真实值，蓝线代表每个模态100次估计结果的平均值。

图6 100次仿真的模态频率估计结果Fig.6 Modal eigenfrequency estimation results from 100 simulations

由图6看出，基于QR分解的数据驱动随机子空间法(SSI/data)识别出4个模态(4行)，而协方差驱动随机子空间法(SSI/cov)仅识别出3个模态(3行)，说明基于QR分解的数据驱动随机子空间法能更敏感地识别弱势模态，估计结果更集中，计算精度更高。图7与图6类似。但模态阻尼估计结果计算精度稍差，相对较分散。

为进一步说明以上分析结果，分别统计模态频率和模态阻尼估计结果均值及变异系数见表1、表2。由二表可见，数据驱动随机子空间法的估计结果更接近真实值，同一模态变异系数更小。进一步表明基于QR分解的数据驱动随机子空间法的优势。

图7 100次仿真的模态阻尼估计结果Fig.7 Damping ratios estimation results from 100 simulations

表1 频率估计值均值与变异系数(定阶为10)Tab.1 Mean and variation coefficient of estimated frequencies

表2 阻尼估计值均值与变异系数(定阶为10)Tab.2 Mean and variation coefficient of estimated damping ratios

由蒙特卡洛试验结果知，无论在识别模态的数目上，或估计的精度上(尤其阻尼)，基于QR分解的数据驱动随机子空间法明显优于协方差驱动随机子空间法。

4 结论

本文分析并验证了造成数据驱动与协方差驱动随机子空间法差异的原因。

通过对数据有限性与QR分解算法的使用、研究发现，基于QR分解的数据驱动随机子空间法无论在识别参数的精度上(尤其阻尼)，或在对较弱势模态的识别能力上均明显优于协方差驱动随机子空间法。因此，本文推荐用基于QR分解的数据驱动随机子空间法作为海洋平台等大型结构的模态参数识别工具。

［1］ Wang S Q，Zhang Y T，Feng Y X.Comparative study of outputbased modal identification methods using measured signal from an offshore platform［C］.Proceedings of the ASME 29th International Conference on Ocean，Offshore and Arctic Engineering，2010，2:561-567.

［2］ Hu S，Li P，Vincent H，et al.Modal parameter estimation for jacket-type platforms using free-vibration data［J］.Journal of Waterway，Port，Coastal，and Ocean Engineering，2011，137(5):234-275.

［3］ Cunha A，Caetano E.Experimental modal analysis of civil engineering structures［M］.lvaro Cunha and Elsa Caetano:University of Porto，2006:12-20.

［4］ Overschee P，DeMoor B.Subspace identification for linear system:theory-implementation-applications［M］. Kluwer Academic Publishers Boston/London/Dordrecht，1996.

［5］ Peeters B，DeRoeck G.Stochastic system identification for operational modal analysis:a review［J］.Journal of Dynamic Systems，Measurement and Control，2001，123(4):659-667.

［6］Peeters B，DeRoeck G.Reference-based stochastic subspace identification for output-only modal analysis［J］.Mechanical Systems and Signal Processing，1999，13(6):855-878.

［7］ Hu S，Bao X，Li H.Model order determination and noise removal for modalparameterestimation［J］.Mechanical Systems and Signal Processing，2010，24(6):1605-1620.

［8］ Matlab.Control system toolbox user's guide［M］.Natic，Ma，2004.