基于随机子空间法的海洋平台模态特征实时提取方法研究

黄 焱，陈 涛，朱本瑞

(1.天津大学 水利工程仿真与安全国家重点实验室，天津 300350；2.天津大学 建筑工程学院，天津 300350)

海洋平台作为海洋油气资源开发的基础性设施，是海上生产作业和生活的基地。为了保障海洋平台结构的安全服役，避免重大的恶性事故发生，必须对海洋平台结构在服役期内进行定期或不定期的检测和安全评估。海洋平台结构模态特征实时提取就是以此目标出发而衍生出的一项重要手段。目前，常用的模态特征提取方法包括峰值拾取法、自然激励技术(Natural Excitation Technique,NET)、ITD(Ibrahim Time Domain,ITD)法、随机子空间方法(Stochastic Subspace Identification,SSI)等。其中，SSI因其计算模型参数少，无迭代运算，收敛性好，运算快捷，鲁棒性好等特点[1-3]，受到广泛关注。

随机子空间方法根据识别过程中相关性计算方式的不同，可以分为数据驱动随机子空间方法和协方差驱动随机子空间方法[4-5]。在同样有限数据的情况下，两种方法识别精度差别不大，但协方差驱动随机子空间方法表现出更高的计算效率[6]。对于长期处于复杂激励下的大型结构物，协方差驱动随机子空间方法在结构模态实时提取中更具优势。Peeters等[7]分析了随机子空间方法的性质，并利用该方法识别受风载作用下的钢杆的振动特性，取得了良好效果。樊可清等[8]推导了随机子空间辨识算法中的协方差矩阵构造方式，通过减少辨识过程中的重复运算量达到缩短运算时间的目的，从而提高了方法的实用性，并通过对香港汀九大桥状态监测系统的数据处理结果加以验证。Reyndersat等[9]研究了随机子空间识别模态参数的不确定性，提出用稳定图方法剔除虚假模态，并通过一个桥梁仿真模型验证了该方法在辨识密集模态时的有效性。

目前，随机子空间方法多运用于陆上钢结构、桥梁、建筑物等工程结构，而对于难以监测水下结构的海洋平台涉猎甚少。辛俊峰等[10]针对一个四桩腿单层甲板海洋平台，研究了冰激振动下的结构模态提取工作，并取得了较理想的效果。但其研究计算的结构比较简单，环境工况的设置单一，忽略了环境条件的多变性及随机性，缺乏对复杂工况下所带来的虚假模态的考虑。针对复杂环境下的海上平台而言，随机子空间方法计算效率不高，有限测点下无法反映结构整体特征等问题尚未有很好的解决办法，因此，如何能够有效获得海洋平台等大型结构模态参数的研究仍有待深入。为此，本文基于海洋平台的结构特点，提出了一种改进的协方差驱动随机子空间模态实时提取方法，为结构健康评估提供参数化指标。

1 理论基础

1.1 离散时间状态空间模型

对于N个自由度的系统，在外界环境激励下，其振动微分方程可用下述方程表示：

(1)

式中:M表示N×N维质量矩阵；C表示N×N维阻尼矩阵；K表示N×N维刚度矩阵；q(t)表示t时刻下的N维位移列向量；f(t)表示N维外部激励向量。

实际测试中，考虑到信号采集的时间离散性和测量时存在的噪声干扰现象，对式(1)进行采样离散化处理以及微分求解后，得到系统离散时间状态空间模型如下：

(2)

式中:yk和xk分别表示第k个采样时刻系统的输出量和状态向量；wk和vk分别表示环境干扰噪声序列和传感器自身测量误差所带来的噪声序列；A∈R2N×2N，表示系统状态矩阵；C∈Rm×2N，表示系统输出矩阵；N为系统阶数，m是所采用的测点个数。

此空间模型是由状态方程和输出方程组合，构成一个系统完整的动态描述。其基本假设前提为：

(1)噪声项为均值为零的平稳随机过程，即

E[wk]=0E[vk]=0

(3)

(2)噪声信号是均值为零的白噪声序列，且与结构真实状态无关，即

(4)

(3)系统为线性时不变系统(既满足叠加原理又具有时不变特性)，状态序列为平稳随机过程，即

(5)

式中：E表示数学期望；状态协方差矩阵Σ与时间k无关。为此，定义输出协方差矩阵Ri如下：

(6)

式中：G表示状态—输出协方差矩阵，其表达式为：

(7)

在实际测试中，通常只选取一段数据量(S个有限数据点)参与计算，因此，输出协方差矩阵可估计为：

(8)

1.2 协方差驱动随机子空间方法计算流程

协方差驱动随机子空间法的计算主要分为五个步骤：首先由响应数据y按反向对角线元素相等规则构成Hankel矩阵，并按矩阵的行列数将Hankel矩阵分为过去部分Yp与将来部分Yf；其次计算Yp与Yf的相关性；然后对计算后的矩阵进行奇异值分解(SVD分解)；再依据分解后结果获得系统状态矩阵A和输出矩阵C；最后通过矩阵A和C求得结构模态参数。

(1)构造Hankel矩阵，并分为过去部分Yp与将来部分Yf[11]：

(9)

式中:yi∈Rm×1，表示第i时刻下所有测点(m个)的响应数据。

(2)构造T型矩阵(Toeplitz矩阵)，计算Yp与Yf的相关性，即

(10)

式中:Ri为输出协方差矩阵。

将式(8)代入式(10)，并对T型矩阵分解可得到：

(11)

式中：表示协方差矩阵T1|i可表示为扩展可观测矩阵Oi∈Rmi×2N与扩展可控矩阵Γi∈R2N×mi的乘积。

(3)T型矩阵奇异值分解(SVD)

在对T矩阵进行分解之前，通常先对其进行加权处理[12]，再对加权过后T型矩阵进行奇异值分解，可得：

(12)

式中：S1∈RN×N为对角矩阵，U、V∈Rmi×2N，是正交矩阵，即满足UTU=VTV=I

(4)系统状态矩阵A和输出矩阵C

比较式(9)和(10)可知：

(13)

根据T型矩阵T1|i构造方式，构造下一时间延迟协方差矩阵T2|i+1

(14)

式中：矩阵T1|i和T2|i+1具有相同的结构，只是包含的协方差矩阵Rk时延变为从2到i+1。

观察矩阵T1|i和T2|i+1可知，输出矩阵C为拓展可观测矩阵Oi的前m行，系统状态矩阵A可用矩阵T2|i+1矩阵表示(*+表示广义逆)：

(15)

(5)系统模态识别

对矩阵A进行特征值分解：

A=ψΛψ-1

(16)

式中:Λ=diag(λ)∈RN×N，是一个对角矩阵，该矩阵中的元素为系统的离散时间复特征值；ψ为离散时间特征向量组成的矩阵。

通过系统矩阵的特征值及特征向量可以得到系统模态参数为

(17)

式中:ξ为阻尼比，ω为圆频率，Φ为振型。

随机子空间方法在识别结构模态信息时，因输入信号不满足白噪声的假定前提或输出信号受到环境干扰，计算时容易出现虚假模态现象。且由于其算法采用状态空间模型，需要预先人为确定系统阶次n，而为了保证计算时避免出现模态遗漏问题，通常需要对系统阶次进行过估计，因而，亦不可避免的产生虚假模态[13]。

对于虚假模态的消除，目前广泛采用稳定图[14]来实现。理由是虚假模态信息不是反映系统本身特征，是由于各种干扰产生，必然不如结构真实模态的计算结果稳定。稳定图的做法是先选取一定常数据段，求不同系统阶数下识别的模态结果反映到一个二维坐标系中(坐标图的横坐标为频率值，纵坐标为阶次)，最后呈现稳定的频率信号值便被认定为结构的真实模态信息。目前，随机子空间辨识法与稳定图相结合的做法在工程应用上取得了较好的效果。

2 随机子空间方法在海洋平台结构的应用

根据随机子空间方法基本原理可知，在识别过程中需要的人工干预少，算法本身抗噪能力强，识别结果的准确度高。但目前随机子空间方法的研究多限于测点能够全覆盖的陆上桥梁等工程结构，对于海洋工程结构，该方法的应用存在诸多阻碍。主要原因在于海洋工程结构既无法做到测点全覆盖布置，又时刻受到各种复杂海况的干扰，导致采集到的信号干扰成分过多并且很难保留结构完整信息，极大地影响了算法的识别效果。因此，随机子空间方法在海洋工程结构上的应用仍有待进一步研究。

本节以渤海某导管架平台为研究对象，来分析随机子空间方法在海洋平台模态参数的识别效果。该平台由六桩腿的导管架钢结构作为支撑，上部组块分为3层，作业水深12.6 m，平台总质量为9 610.4 t，见图1。采用有限元软件ANSYS建立该平台结构有限元模型，见图2。建模时对井口区隔水套管进行模拟，并充分考虑平台历次改造的结构与设备布置变化，保证了各层甲板重心以及模型的总重量重心与平台设计及改造资料一致，使模型在各方面尽可能接近平台结构的真实状态。采用随机波浪理论模拟现实海况，对模型平台施加随机波浪激励，模拟各种海况下平台所受的波浪冲击。根据该平台所处海域水文资料统计结果，90°波浪入射方向为该平台的主浪向，其它出现概率较大的波浪方向为0°、30°，波浪入射方向见图3。为尽量覆盖平台可能遭遇的波浪环境，针对这3种波浪入射方向均设置了3种出现概率最大的波高条件，有效波高分别为0.5 m，1.0 m，1.5 m，不同工况对应编号,如表1所示。

图1 渤海某导管架平台示意图

图2 平台有限元模型及测点位置示意图

图3 平台水线面处波浪入射方向示意图

表1 有限元模型计算工况及其编号

现实监测中，考虑到传感器的适用性与施工可行性，测点通常布置在平台上部甲板上。本文选取的目标平台共有三层主甲板六根桩腿，按照一般的监测原则，每层甲板可布置6个测点，共计18个测点。对应于现实结构，数模计算时，选取平台结构上部三层甲板与六根桩腿相交的18个节点作为加速度信号采集点，即测点数目m=18，采样频率为10 Hz，采样时间为9 000 s。

通过ANSYS模态分析得出平台前三阶频率信息，其中包括x方向、y方向以及扭转方向(z)的模态。如表2所示。

表2 导管架平台有限元模型模态信息

面对目标平台如此复杂的结构，经过前期的估算，发现至少要形成18个测点的布置方案，才有可能对系统形成较好的识别。针对该布置方案，采用传统的随机子空间方法的处理方式是：选用监测结构响应的所有18个测点按顺序组成一列向量构建Hankel矩阵，再利用随机子空间算法进行模态识别，最后结合稳定图提取结构稳定模态。九种工况下系统的识别结果如图4所示。需要说明的是，由于海洋工程结构的复杂性与测点布置的局限性，可以预判此时所形成的模态识别结果中必然会含有一定数量的虚假模态，因此，为与上面给出的结构真实模态计算结果相区分，将这些识别结果称作准模态。考虑到均值可以表征准确度，而方差则反映了稳定性水平，因此识别结果中分别给出了每组数据的均值μ与方差σ，如表3所示。

图4 传统方法不同工况稳定图

如表3数据所示，其中准五阶，准六阶，准七阶频率与结构前三阶真实频率信息有很好的对标性，但其余的准一至四阶与准八阶模态信息并不是结构本身所拥有的，即传统方法下出现了无法消除的虚假模态信息。

由表3中数据和图4可知，在各种工况下，传统随机子空间法不仅在真实模态频率上出现了稳定点，在许多非系统模态频率上也表现出稳定点集中现象。且稳定图中数据点波动较大，即各种工况下的识别结果存在较大波动。这主要是由于测点无法实现结构全覆盖，导致算法对结构信息识别的“发散式”推测，以及对结构系统阶数的过估计所致。

表3 不同工况频率识别结果

由此可见，传统方法往往会导致虚假的模态提取。因此，如何基于随机子空间方法基本原理，发展一套能在有限测点下增强反映结构整体特征，并且有效消除虚假模态的改进方法仍需进一步讨论。

3 改进随机子空间方法

海洋平台工作环境的限制，决定了其不可能和陆地结构物一样尽可能采用测点全覆盖的方式来进行监测信号获取。因此，如何根据平台水上结构有限测点数据信息推测整体模态特征是随机子空间法运用的关键技术难点。由此目标出发，本文提出了两点改进：Hankel矩阵元素重构，即每个元素重新嵌套子矩阵的方法和基于信号累积特征的虚假模态辨识方法与剔除准则。

3.1 Hankel矩阵元素重构

随机子空间方法的核心是Hankel矩阵的建立，其描述的是动态响应的时间连续继承关系，在数学上体现为过往响应与未来响应间的相关性。因此，当Hankel矩阵中描述的相关性仅是有限的结构片段时，就会在预先设定的模态阶数所控制的满秩矩阵运算中，产生对结构剩余部分的“发散式”推测。由随机子空间方法的基础理论可知，数字矩阵间的协方差关系是形成结构信息描述的基本手段。那么可以推知，如果将局部位置(测点)的数据作为独立的Hankel矩阵元素进行运算，实际上是抹除了局部位置间本来可能具备的关联性而进行全新的关联性重构，这也是造成“发散式”推测的根本原因。因此，对于现实中监测数据仅为结构的片段信息时，就必须在Hankel矩阵的构建中尽可能保留甚至增强反映结构整体特征的数字信息量，即以某种测点数据排列方式构建一个关联子矩阵作为Hankel矩阵运算的基本元素。

在结构有限测点的情况下，为了在Hankel矩阵中尽可能保留结构整体特征的数字信息量，可以在选择测点数据排列方式时尽可能的对标测点在实际结构中的空间排布，这里主要考虑两方面对标性：

(1)子矩阵元素排布对标实际测点空间排布：Hankel矩阵中每个元素yi表示在第i个时刻下各个测点响应信号值，那么只需要把各个测点在某个平面内的空间排布位置与子矩阵中各元素位置对应，这些测点便可以按一定的排列方式组成一个子矩阵，用来取代Hankel矩阵中对应元素。为了能够充分收集到结构各个方向的模态信息，可以选取结构各个横剖面和纵剖面内的测点组成多种子矩阵交替计算形式。为了能更清楚表明结构测点位置信息，将上述导管架平台测点位置表示成如图5的形式。

图5 测点位置示意图

工程实际中，不同的测点位置对结构整体信息敏感程度也不一样，如果仅选取一个面内的测点响应数据不足以表达结构的整体信息。为了尽可能方便提取结构所有方向模态信息，把能表达结构不同方向信息的横剖面和纵剖面测点数据分别构成多个子矩阵形式，具体如图6所示。图中，yki表示第k个传感器在第i时刻下的响应数据。

图6 子矩阵形式

(2)子矩阵元素之间关联性的保持对标响应数据获取方向：按照上述子矩阵构建方式，已经完成了测点位置信息的表达，但子矩阵内部仍然缺少元素方向关联性(运算时各子矩阵内沿同一方向振动的测点数据不被打散)的表达。为了补全这种方向关联性，子矩阵内部元素还需进行进一步的向量重构，即把子矩阵元素沿振动方法合并成多个行向量。所以子矩阵构建可进一步表达为：

(18)

其中:

(19)

(20)

或

(21)

可以看出，相比较传统随机子空间方法把所有m个测点信号作为一个列向量输入，改进后的Hankel矩阵有明显的两个优势：子矩阵中尽可能的包含了结构整体特征的数字信息，使得计算结果更能反映结构的真实模态信息。另外，新的输出协方差矩阵由原来的m×m维变为现在的m/t×(m/t)维(t为Hankel矩阵子矩阵列数)，极大地缩减了后续矩阵运算维度，明显的提高了随机子空间方法的计算效率。

为验证元素重构后的识别效果，将上述五种子矩阵组合的构造方式运用于模型数据计算，分别采用随机子空间算法，对各自识别结果进行稳定图处理，得到五组稳定图数据。将这些稳定图数据进行均化合并，得到一组综合稳定图，九种工况下识别结果,如表4所示。

表4 不同工况频率识别结果

对比表3和表4可知，相比较传统方法识别结果，Hankel矩阵元素重构后，表现出多方面的优点。①改进后方法中每一阶次数据的方差下降十倍以上，识别结果数据稳定程度大幅提升。②改进后方法识别的准五阶，准六阶，准七阶频率信息，其结果也更加接近结构的前三阶模态真实值，增加了识别结果的准确度。③经过Hankel矩阵元素重构，新的Toeplitz矩阵由原来的mi×mi维变为现在的mi/t×mi/t维(t为Hankel矩阵子矩阵列数)，即新的协方差矩阵维度缩减了3×3倍，极大地缩减了后续矩阵运算维度。在相同计算条件下(戴尔DELL T630塔式服务器，32 GB内存，双英特尔处理器)对两种算法分别完成一组模态识别所需的计算时间进行比较，结果表明：传统随机子空间算法用时34.8 min，而改进算法仅用10.5 min，计算效率提高了三倍。综合而言，通过Hankel矩阵元素重构后，改进算法既提高了模态识别的准确度与稳定性，又极大地的提升了计算效率，具有良好的优越性。

应注意的是，不同响应数据经过计算处理后，很多非系统真实特征的信息在稳定图上亦呈现出很稳定的结果，如表4中的准一阶，准二阶，准三阶，准四阶，准八阶模态。此外，在识别过程中，也很容易使得结构真实模态信息在不同的系统阶数下表现为不同的阶次，即出现模态越阶现象。那么这时候，传统的稳定图方法会把结构真实模态信息消除，从而出现误判。为此，必须引入一种新的能够有效消除虚假模态的准则，来辅助提取结构真实特征信息。

3.2 基于信号累积特征的虚假模态辨识方法与剔除准则

由于虚假模态的出现，是由于计算过程中存在各种干扰，从而引起的对结构模态信息的发散性推测。它反映的不是系统本身特征，且不同时段下结构所受到的干扰不同，因此，随着时间的推移，虚假模态的不稳定现象会更加突出，而真实模态则具有模态鲁棒性的时间累积增强效应。因此，可以通过采用不同时间段的数据计算结果来统计结构信息的稳定程度，从而甄选出结构的真实模态信息。

基于此原理，本文提出一种基于信号累积特征的虚假模态辨识方法与剔除准则，主要是通过不同时段数据结合不同子矩阵的形式依次做随机子空间方法计算，最后把识别结果在各个小区间下做分段统计。那么，结构的真实模态必然会在一个区间范围内保持较好稳定性，而虚假模态则会相对杂乱的排布在其他小区间，其表现出来的稳定性自然也会越差。其基本步骤如下：

步骤1：选取一较长时间段数据S并按数据长度均分成多段数据Si。

步骤2：选取结构不同截面测点组成j种子矩阵构建形式，并把把每段数据Si有时间交叉的分割成j小段数据Sij(数据Sij间有80%重叠率，使得结构真实模态信号更稳定而虚假模态产生较大变化，数据Sij表示在数据Si下的第j小段数据)。

步骤3：每一段数据Sij对应一种子矩阵规则进行计算，得到j种Hankel矩阵，依次做子空间运算识别结构模态，并把整体数据段S下的频率计算结果在各个区间下做分段统计，即每组识别结果按数值大小划分在不同区间(为使各区间足够把结构各阶次模态区分开，定义每个区间的长度为0.02 Hz)。

步骤4：把包含数据量少的区间作为不稳定区间(定义每个区间内包含P种子矩阵的识别结果以及每个区间内包含Q个数据，取P、Q阈值为p、q，则各区间数据结果小于p、q的定义为不稳定区间)，抛除不稳定区间，留下稳定区间。

步骤5：统计稳定区间内前几组连续区间，其中每组连续区间内数据期望值便作为结构的真实模态信息模态。

针对文中算例，采用上述步骤进行虚假模态提取，其流程图如图7所示。取数据S为9 000 s，均分成五组S1，S2，S3，S4，S5，并把每组数据Si在80%数据重叠率下分成五段数据Sij，每一段数据Sij各选取一种子矩阵规则进行运算，交替做随机子空间运算，识别出多组结构模态。最后把整体数据S下的频率计算结果在各个区间下做分段统计，把包含数据量少的区间作为不稳定区间。定义不稳定分段时，经过多次验证，P、Q阈值分别取4和6。抛除不稳定区间，留下稳定区间，其中每组连续稳定区间内数据期望值便作为结构的真实模态信息。以各连续稳定区间内数据数学期望值为横坐标，各连续稳定区间内数据总量为纵坐标，统计各工况模态计算结果，如图8所示。

图7 模态提取流程图

图8 不同工况频率识别结果

依据上述方法，对结构9种工况下响应数据识别结果进行分析。得出目标平台的前三阶频率识别结果。如表5和图9所示。

观察表5和图9可以发现：经过Hankel矩阵重构与虚假模态消除处理后，最终识别出的准一阶，准二阶，准三阶模态与结构前三阶真实模态有很好的对标性，有效剔除了传统方法下的众多虚假模态，识别结果准确度也保持在95%以上。对于有限测点下的海洋平台结构，改进后的随机子空间方法有很好的优越性。

图9 不同工况频率识别结果

表5 不同工况频率识别结果

(1)改进后的随机子空间算法能有效识别结构各阶次模态结果，并把识别结果集中在一定动态范围之内，有效排除了众多虚假模态的干扰，准确提取出系统模态所在区间。

(2)相对于Hankel矩阵元素重构后的识别结果而言(不考虑虚假模态的情况下)，稳定段内识别出的模态信息与系统真实模态非常接近，而且识别结果的波动也比较平稳。也就是说，在不损失计算结果稳定程度的基础上，识别结果准确度有了明显的提升。

总而言之，改进后的随机子空间方法识别结果与理论分析保持很好的一致性。采用子矩阵与信号累计的形式，使计算结果尽可能的包含了结构整体特征的数字信息，极大地缩减了矩阵运算维度，明显的提高了随机子空间方法的识别准确度和计算效率。相比原重复性计算同一组数据下做稳定图，能更加明显的排除虚假模态干扰以及有效提取结构真实模态信息。

4 结 论

针对海洋平台模态识别问题，本文在传统随机子空间方法理论基础上，提出了新的Hankel矩阵元素重构方式以及虚假模态辨识方法与剔除准则。为了更好反映结构各个方向的整体数字特征信息，在原有Hankel矩阵元素上重新嵌套新的子矩阵。对标测点实际空间位置建立多种子矩阵组合方式，并结合有时间重叠的响应数据做交替计算。为消除识别结果中的众多虚假模态以及修正各组次下相近模态计算的越阶现象，把所有识别结果采用各区间分段统计的方式来处理。通过对目标平台有限元模型响应信号进行计算，在有限测点下情况下，改进的方法能更加明显的排除虚假模态干扰并有效提取结构前三阶真实模态信息，也极大地提高了原有算法计算效率，具有良好的适用性和鲁棒性。