数据驱动随机子空间法矩阵维数选择与噪声问题研究

辛峻峰，盛进路，张永波

(1.中国海洋大学 工程学院，青岛 266100;2.大连海事大学 航海学院，大连 116026;3.青岛国家海洋科学研究中心，青岛 266071)

有效的模态识别方法能精确地得到结构的相关参数，从而可以准确地掌握大型结构的健康状况。国内外许多学者［1-3］提出了基于时域响应的模态参数识别方法，如时间序列法、随机减量法、自然激励技术、数据驱动随机子空间法等。其中，数据驱动随机子空间法是目前为止较为先进的环境激励下模态参数识别方法之一，它能够更精确地提取海洋平台等大型结构的模态参数。但是，确定Hankel矩阵的维数是该方法有效应用的关键［4］，不同的Hankel矩阵维数会导致数据驱动随机子空间法消噪能力的变化。如何确定Hankel矩阵的维数?目前对这个问题的研究还鲜见报道。针对于此，本文推导了Hankel矩阵的维数与数据驱动随机子空间法消噪能力的理论关系，继而提出了数据驱动随机子空间法Hankel矩阵维数的选择方法，探讨了不同的Hankel矩阵构建方式与数据驱动随机子空间法消噪能力的关系。

1 数学理论

1.1 数据驱动随机子空间识别法

在仅考虑随机噪声的前提下，振动系统的离散状态空间方程可表示为:

其中:xk∈Rn表示在离散时间k时，系统的状态向量，n为系统的阶数;A∈Rn×n表示离散状态矩阵;C∈Rl×n表示输出矩阵，描述内部状态怎样转化为外界测量值，l表示测点的数目;wk∈Rn表示由于干扰和模型误差造成的过程噪声;vk∈Rl表示由于传感器误差等造成的测量噪声。这两种噪声都是不可测量的，在推导过程中经常被假设为零均值，平稳的白噪声。两种噪声的协方差矩阵可以用下式表示:

其中:E表示数学期望算子;Q∈Rn×n，S∈Rn×l，X∈Rl×l;δpq为 Kronecker函数，δpq=1(p=q)，p和q表示不同的时间点，δpq=0(p≠q);E［wp］=0，E［vp］=0。

数据驱动随机子空间识别法详细讨论可见文献［5-7］。

1.2 Hankel矩阵维数选择与噪声关系

通常数据驱动随机子空间法的投影矩阵P可以进一步表示成如下:

其中:S代表真实信号;N代表需要剔除的噪声;i代表Hankel矩阵1/2行数;j代表矩阵列数;m代表数据总量。

继而在奇异值分解(SVD)后，设定阈值(模态阶次)n，式(3)变为:

其中:PS=US∑S，PN=UN∑N，U∈Ri×i，∑∈Ri×j，V∈Rj×j，∑代表分解后测量信号投影后包含的所有奇异值;∑S代表设定阈值n后，真实信号投影后包含的奇异值。∑N代表设定阈值n后，噪声信号投影后包含的奇异值。由式(4)可得，奇异值分解可以将含噪声的投影矩阵P划分为两个互不相关的空间，即真实信号投影PS和噪声信号投影PN。下面讨论式(4)中噪声投影PN与所有测量信号投影P的维数(行数i或者列数j)及其奇异值分解阈值n两个参数之间的关系。

根据式(4)有:

由于真实信号与噪声不相关，所以:

由式(5)及式(6)可得:

因为测量信号、真实信号以及噪声的投影的协方差分别为:

式(7)可改写为:

结合式(4)与式(8)，式(9)可以写为:

由式(8)到式(10)可得:

式(11)左右两边同时左乘iVT和右乘V得:

式(12)中，在数据总量m不变的前提下，V与∑S取决于Hankel矩阵维数(行数2i或者矩阵列数j，2i+j-1=m)和投影矩阵奇异值分解的阈值n，∑取决于Hankel矩阵维数(行数2i或者矩阵列数j，2i+j-1=m)。所以，如果投影奇异值分解的阈值n不变，Hankel矩阵的维数是影响噪声CN的主要因素。

1.3 稳定图

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

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

2 数据驱动随机子空间法－矩阵维数设定评估方法

归一化奇异值(SVD)、稳定图以及有限元模态识别结果(FE)三维一体评估方法。操作流程和每一步的功能性如图1所示。

3 Hankel矩阵构建方案

设数据总量为m(2i+j-1=m)，分别对四种Han-kel矩阵维数设定方案(2i/m=2/10，3/10，4/10，5/10)进行探讨，其中，方案一(2i/m=2/10)，表示 Hankel矩阵 H是一个较扁平的非方阵;而方案四(2i/m=5/10)，表示Hankel矩阵H是一个方阵，因此方案一到方案四的变化过程也是Hankel矩阵从一个扁平的非方阵逐渐变为方阵的过程。

图1 评估方法操作流Fig.1 Operation process of evaluation method

4 数值算例

本节将通过五自由度质量-弹簧-阻尼系统白噪声激励试验(模型如图2所示)，探讨不同Hankel矩阵的维数对数据驱动随机子空间法消噪能力的影响，同时验证SVD、稳定图、FE组合评估方法的有效性。

4.1 五自由度质量－弹簧－阻尼系统

每个单元有相同的质量、刚度和阻尼系数，分别被设定为:mn=50 kg，kn=2.9 ×107N/m，cn=1 000 Ns/m。xn代表位移，并且n=1，…，5。通过特征值分析，得到 5阶模态频率的理论值为:34.499、100.700、158.730、203.880、232.520 Hz;5 阶模态阻尼比的理论值为:0.003 737、0.010 909、0.017 197、0.022 092、0.025 198［8-9］。

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

4.2 白噪声激励加载及数据获取

采用Matlab中randn函数生成均值为零，标准差为1.006的高斯白噪声作为输入(图3)将其右向加载在靠固定端的第一个质量块上(图2)，从lsim函数生成振动响应，以500 Hz的采样频率提取第一个质量块的1 024个响应数据，如图4所示。

4.3 不同方案对比

首先使用归一化奇异值对比不同方案能识别出的模态的数目。

图3 输入的白噪声序列Fig.3 White noise loading

图4 靠固定端质量块响应Fig.4 Response of mass close to fixed point

图5 四种方案归一化奇异值曲线Fig.5 Normalized singular values in four cases

归一化奇异值的结果(图5)清晰地显示出:当Hankel矩阵从扁平(方案一，2i/m=2/10)趋于方阵(方案四，2i/m=5/10)的过程中，数据驱动随机子空间法能清晰识别的模态数目从4个(方案一，红线最后的较大落差在横轴为8的时候)减少到1个(方案四，紫色线最后的较大落差在横轴为2的时候)。这说明Hankel矩阵趋于方阵的过程中，数据驱动随机子空间法的消噪能力逐渐减弱，从而导致了其识别出的模态越来越少。

接下来，利用稳定图分析结果对比不同方案下识别出的模态的稳定性，如图6所示。

图6 SSI/data四种方案的稳定图分析结果(a)～(d)对应着方案一，2i/m=2/10到方案四2i/m=5/10，代表稳定，o代表不稳定)Fig.6 Stability Diagrams for SSI-data in four cases(a-d responding to case 1，2i/m=2/10 to case 4，2i/m=5/10，*stable，°not stable)

稳定图清晰地显示了在Hankel矩阵趋于方阵的过程中(图6，(a)-(d)):

(1)数据驱动随机子空间法逐渐地不能稳定识别出区间在200 Hz左右的第四阶模态和区间在230 Hz左右的第五阶模态(代表不稳定的圆圈逐渐增多;对应峰值的直线逐渐消失或者发生混乱)，这个现象说明数据驱动随机子空间法消噪能力在逐渐减弱，从而导致其对弱一些的模态越来越不敏感。

(2)数据驱动随机子空间法对最强的两个模态(区间分别在30 Hz和100 Hz左右的第一阶模态和第二阶模态峰值最大)的估计逐渐变的不稳定(代表不稳定的圆圈逐渐替代代表稳定的星号)，这证明了数据驱动随机子空间法的消噪能力在逐步减弱，从而导致了其对较强模态的识别也减弱了。

表1 SSI/data四种方案下的模态频率估计结果Tab.1 Frequencies estimated in four cases using

最终，通过参考有限元的模态分析结果(FE)，确定不同方案识别出的稳定模态的真假。表1和表2分别罗列了不同方案估计的模态频率和阻尼比。

表2 SSI/data四种方案下的模态阻尼比估计结果Tab.2 Damping ratios estimated in four cases using

表1和表2，不但确定了方案一识别出的5个模态的真实性，而且能够观察到:Hankel矩阵从方阵(方案四)趋于非方阵(方案一)的过程中，数据驱动随机子空间法对模态参数的估计逐渐精确，尤其是对模态阻尼比的估计。

此外，本文分别使用了含有不同噪声水平的响应数据，进行分析得到的结论与上文一致。

总之，数值算例证明了:① 数据驱动随机子空间法的Hankel矩阵应该设置成非方阵的形式;② SVD、稳定图以及FE相结合的方法能够有效评估不同Hankel矩阵维数对数据驱动随机子空间法模态识别能力的影响。

5 导管架平台振动台试验

本节将进行钢导管架平台物理模型白噪声激励试验进一步验证上文的研究结论。试验方案为:使用液压振动台产生白噪声激励，在钢导管架平台物理模型底部加载(与x轴和y轴成45度角)，获得结构的振动响应数据，继而对比不同方案下的数据驱动随机子空间法的消噪能力。

5.1 导管架平台物理模型及试验

本导管架物理模型采用钢管制作，在Y向设3根主梁，X向设2根，有3层水平横撑，分别位于-88 cm、142 cm和188 cm处设水平横撑，在88 cm至-142 cm处设竖向斜撑(图7)，主要构件的尺寸为:主腿25×2.5 mm，水平撑杆 16 ×1.5 mm，斜撑杆 16 ×1.5 mm［10］。

钢导管架平台有限元模型的前9阶模态如表3所示。

用液压振动台生成白噪声激励，如图8所示。以500 Hz的采样频率，在y方向的传感器A提取了1 024个数据点(图8)。

图7 导管架平台物理模型，白噪声激励加载位置以及传感器位置Fig7.Physical model of jacket platform，white noise ground motion loading position and the sensor A location

表3 有限元模型前9阶模态频率Tab.3 The first nine modal frequencies of FE

图8 白噪声激励传感器A的响应信号Fig.8 Response of sensor A from white noise ground motion

5.2 不同方案对比

跟数值算例中的分析过程一致，首先用归一化奇异值对比不同方案。

归一化奇异值(图9)清晰地显示出:当Hankel矩阵H趋于方阵的过程中，数据驱动随机子空间法能识别出的模态数目不断增加。从4个(方案一(2i/m=2/10)，红线最终较大落差在横轴为8的时候)逐渐增加到6个(方案四(2i/m=5/10)，紫色线最终较大落差在横轴为12的时候)。这结果似乎在表明方阵的Hankel矩阵使数据驱动随机子空间法识别更多的模态，为了证实这一点，下文结合稳定图继续分析。

图9 导管架试验，SSI/data四种方案归一化奇异值曲线Fig.9 Normalized singular values in four cases in jacket platform model test

图10 导管架试验，SSI/data四种方案的稳定图分析结果(图a-d对应着2i/s=2/10到2i/s=5/10，*代表稳定，o代表不稳定)Fig.10 Stability Diagrams for SSI-data in four cases in jacket platform model test.(a-d responding to case 1，2i/s=2/10 to case 4，2i/s=5/10，*stable，°not stable)

伴随着Hankel矩阵趋于方阵(图10，a到d)，从稳定图中可以观察到:

(1)在30 Hz、50 Hz、130 Hz三个最大的峰值处，数据驱动随机子空间法的识别结果越来越不稳定(代表不稳定的圆圈逐渐增多。

(2)方案一(2i/m=2/10，图10(a))与方案二(2i/m=3/10，图10(b))都能稳定地识别出4个模态，比其他方案识别的模态多，但是它们在第四阶模态的识别上有所不同:方案一识别出的第四阶模态在60 Hz左右，而方案二在130 Hz左右。

(3)只有方案一(2i/m=2/10，图10(a))与方案二(2i/m=3/10，图10(b))稳定图的分析结果与归一化奇异值分析一致:能识别出4个模态，这个现象表明本文提出的评估方法的第二步(稳定图分析)的必要性。

在稳定图的分析中，方案一与方案二的识别出的模态较多，但是在第四阶的模态识别上出现了不一致。由此，需要FE作为参考来辨别识别出的模态的真假。表4和5中罗列了不同方案的评估结果(模态阶数为20)。

表4 SSI/data四种方案下的模态频率估计结果Tab.4 Frequencies estimated in four cases using SSI/data

表5 SSI/data四种方案下的模态阻尼比估计结果Tab.5 Damping ratios estimated in four cases using SSI/data

在表4和5中，参照FE的结果，可以观察到:

(1)第四阶模态应该是在63 Hz左右的扭转模态，方案一(2i/m=2/10)的识别结果是正确的。

(2)只有方案一(2i/m=2/10)能识别出63 Hz左右的第四阶模态，证明了方案一的优势，这个结果表明了本文提出的评估方法第三步(FE验证)的必要性。

总之，导管架平台白噪声试验与数值试验的分析结果一致:本文提出的评估方法有效而且非方阵的Hankel矩阵使数据驱动随机子空间法具备更强的消噪能力和更高的模态识别精度。

6 结论

本文利用SVD、稳定图和FE相结合的方式对数据驱动随机子空间法的Hankel矩阵维数进行选择。其中，首先用SVD对比不同Hankel矩阵下数据驱动随机子空间法能识别出的模态数目，筛选出能识别较多模态的Hankel矩阵;继而使用稳定图对比模态的稳定性，进一步筛选出稳定模态最多的Hankel矩阵;最后参考FE的结果辨别稳定模态的真假，再进一步筛选出真模态最多的Hankel矩阵。基于此方法，通过数值分析和导管架平台振动台试验，探讨了数据驱动随机子空间法的Hankel矩阵维数与其消噪能力之间的关系。系统地证明了:为了最大程度地提高数据驱动随机子空间法的消噪能力和对弱模态的敏感度，Hankel矩阵应采用非方阵的形式(即行数大于列数)。本文的工作可为今后数据驱动随机子空间法的有效推广和应用提供参考。

［1］Wang S Q. Comparativestudyofoutput-based modal identification methods using measured signal from an offshore platform，Proceedings ofthe ASME 29th International Conference on Ocean［C］. Offshore and Arctic Engineering，2010.

［2］ Hu S L J，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:1943-1963.

［3］ Cunha A，Caetano E.Experimental modal analysis of civil engineering structures［J］.Sound and Vibration，2006:12-20.

［4］刘进明，应怀樵，沈 松，等.时域模态分析方法的研究及软件研发［J］.振动与冲击，2004，23(4):123-126.

LIU Jin-ming，YING Huai-jiao，SHEN Song，et al.Method research and software developing of time domain modal analysis［J］.Journal of Vibration and Shock，2004，23(4):123-126.

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

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

［7］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.

［8］ Hu S L J，Bao X，Li H.Model order determination and noise removal for modal parameter estimation［J］.Mechanical Systems and Signal Processing，2010，24:1605-1620.

［9］ The MathWorks，Control system toolbox user’s guide，The MathWorks Inc.，Natic，MA，2004.

［10］ ANSYS.Inc，ed.，2004.ANSYS Academic Research.Release 9.0.Help System.Coupled Field Analysis Guide10.