一种仅输出线性时变结构模态参数辨识方法.pdf

本发明公开的一种仅输出线性时变结构模态参数辨识方法，属于结构动力学技术领域。本发明首先推导出最小二乘支持向量机矢量时变自回归模型的代价函数；利用Wendland紧支径向基函数构造函数空间基于Gamma测试的非参数方法确定正则因子，基于实际经验给出基函数宽度减缩系数；依据贝叶斯信息量准则和赤池信息量准则确定时变自回归模型阶数；依据残差平方和与序列平方和的比值确定函数空间阶数；最后根据代价函数求最小二乘支持向量机矢量时变自回归模型系数矩阵表达式，并根据时间冻结方法求系统的模态频率，完成线性时变结构模态参数的辨识。本发明能够提高计算效率，增强系统鲁棒性，在结构动力学工程应用中广泛应用于线性时变结构的模态辨识。

1.一种仅输出线性时变结构模态参数辨识方法，其特征在于：包括如下步骤，
步骤1：推导出最小二乘支持向量机矢量时变自回归模型(LS-SVM-VTAR)的代价函数Lp
(aij,lm,αl[k],el[k])，具体包括如下步骤；
步骤1.1：推导最小二乘支持向量机模型(LS-SVM)的代价函数L(w,b,e,α)如式(1)所
示：

式中，L(w,b,e,α)表示代价函数，w为参数矢量，b为实数，e＝[e1 e2 … eN]T为误差向
量，α＝[α1 α2 … αN]T为拉格朗日乘子向量，γ为正则因子，为支持向量机训练样
本数据，为N维至更高nh维的映射，上标T表示转置运算，N为训练样本点个
数；
步骤1.2：将矢量时变自回归模型(VTAR)的系数投影到由径向基函数序列表示的函数
空间中；
矢量时变自回归模型(VTAR)如式(2)所示：
$x\[k\]=-\underset{i=1}{\overset{n_a}{\Σ}}{A}_i\[k\]x\[k-i\]+e\[k\]---\left(2\right)$
式中，为系统N0个通道的输出向量，k为第k个时刻点，na为VTAR模型阶数，e
[k]为k时刻的误差或不可观测的非平稳扰动，Ai[k]为第i阶与时刻相关的VTAR系数矩阵；
将矢量时变自回归模型(VTAR)的系数投影到如式(3)所示的函数空间中：

式中，pa为函数空间阶数，fj(j＝1,2,…pa)为第j阶函数向量；将Ai[k]展开，即得到函数
空间表示的矢量时变自回归模型(VTAR)，如式(4)所示：
$x\[k\]=-\underset{i=1}{\overset{n_a}{\Σ}}\underset{j=1}{\overset{p_a}{\Σ}}{A}_{ij}{f}_j\[k\]x\[k-i\]+e\[k\]---\left(4\right)$
x[k]的分量形式如式(5)所示：
$x_l\[k\]=-\underset{i=1}{\overset{n_a}{\Σ}}\underset{j=1}{\overset{p_a}{\Σ}}{f}_j\[k\]\underset{m=1}{\overset{N_o}{\Σ}}{a}_{ij,lm}{x}_m\[k-i\]+e_l\[k\]---\left(5\right)$
式中，aij,lm为矩阵Aij的元素，l表示第l个输出通道；
步骤1.3：得出最小二乘支持向量机矢量时变自回归模型(LS-SVM-VTAR)的代价函数Lp
(aij,lm,αl[k],el[k])；
结合式(1)和式(5)得到所需的代价函数Lp(aij,lm,αl[k],el[k])如式(6)所示：
$\begin{array}{c}{L}_p(a_{ij,lm},{\mathrm{\α}}_l\[k\],e_l\[k\])=\frac{1}{2}\underset{i=1}{\overset{n_a}{\Σ}}\underset{j=1}{\overset{p_a}{\Σ}}\underset{l=1}{\overset{N_o}{\Σ}}\underset{m=1}{\overset{N_o}{\Σ}}{a}_{ij,lm}^2+\mathrm{\γ}\frac{1}{2}\underset{k=1}{\overset{N}{\Σ}}\underset{l=1}{\overset{N_o}{\Σ}}{e}_l^2\[k\]\\ -\underset{k=1}{\overset{N}{\Σ}}\underset{l=1}{\overset{N_o}{\Σ}}{\mathrm{\α}}_l\[k\](-\underset{i=1}{\overset{n_a}{\Σ}}\underset{j=1}{\overset{p_a}{\Σ}}{f}_j\[k\]\underset{m=1}{\overset{N_o}{\Σ}}{a}_{ij,lm}{x}_m\[k-i\]+e_l\[k\]-x_l\[k\])\end{array}---\left(6\right)$
步骤2：利用Wendland紧支径向基函数构造步骤1中式(3)所示的函数空间使函数空
间变得稀疏；
采用的Wendland紧支径向基函数fj[k]具体形式如式(7)所示：

式中，Nw是给定的Wendland紧支径向基函数序列的非零部分，且δ是
径向基函数序列的无量纲宽度缩减因子；是φ1,λ[r]的离散形式，φ1,λ[r]表达式如
式(8)所示：
${\mathrm{\φ}}_{1,\mathrm{\λ}}\left(r\right)={(1-r)}_{+}^{2\mathrm{\λ}+1}{p}_{\mathrm{\λ},\mathrm{\λ}+1}\left(r\right)---\left(8\right)$
λ是多项式的次数，pλ,λ+1(r)是λ次的多项式，运算(a)+＝max{a,0}；当λ取值分别为0,1,
2,3时，pλ,λ+1(r)的表达式分别如式(9)至式(12)所示；
φ1,λ(r)＝(1-r)+ (9)
${\mathrm{\φ}}_{1,\mathrm{\λ}}\left(r\right)={(1-r)}_{+}^3(1+3r)---\left(10\right)$
${\mathrm{\φ}}_{1,\mathrm{\λ}}\left(r\right)={(1-r)}_{+}^5(3+15r+24r^2)/3---\left(11\right)$
${\mathrm{\φ}}_{1,\mathrm{\λ}}\left(r\right)={(1-r)}_{+}^7(15+105r+285r^2+315r^3)/15---\left(12\right)$
步骤3：基于Gamma测试的非参数方法确定步骤1中公式(6)中正则因子γ，基于工程经
验给出步骤1中公式(6)中基函数宽度减缩系数δ；
步骤4：确定步骤1式(6)中最小二乘支持向量机矢量时变自回归模型(LS-SVM-VTAR)阶
数na和函数空间阶数pa；
步骤5：根据步骤1中式(6)所示的代价函数Lp(aij,lm,αl[k],el[k])求步骤1式(2)中矢量
时变自回归模型(VTAR)系数矩阵Ai[k]表达式，并根据时间冻结方法求系统的模态频率，完
成线性时变结构模态参数的辨识；
将式(6)中的代价函数Lp(aij,lm,αl[k],el[k])对各个自变量aij,lm,αl[k],el[k]求导，并
整理得：
$({\mathrm{\Ω}}_l+\frac{I}{{\mathrm{\γ}}_l}){\mathrm{\α}}_l={\tilde{x}}_l,(l=1,2,...N_0)---\left(13\right)$
$a_{ij,lm}=-\underset{k=1}{\overset{N}{\Σ}}{\mathrm{\α}}_l\[k\]f_j\[k\]x_m\[k-i\]---\left(14\right)$
式中，γl是第l个输出通道对应的正则因子，I为N维单位矩阵，αl、xl、Ωl分别如式(15)、
(16)所示：
${\mathrm{\α}}_l=\left[\begin{array}{c}{\mathrm{\α}}_l\left(1\right)\\ {\mathrm{\α}}_l\left(2\right)\\ \·\\ \·\\ \·\\ {\mathrm{\α}}_l\left(N\right)\end{array}\right],{\tilde{x}}_l=\left[\begin{array}{c}{x}_l\left(1\right)\\ x_l\left(2\right)\\ \·\\ \·\\ \·\\ x_l\left(N\right)\end{array}\right]---\left(15\right)$

其中，⊙表示对应元素相乘运算，F矩阵和Rm矩阵元素如式(17)所示：
$F_{ks}=\underset{j=1}{\overset{p_a}{\Σ}}{f}_{j,k}{f}_{j,s},R_{m,ks}=\underset{i=1}{\overset{n_a}{\Σ}}{x}_{m,k-i}{x}_{m,s-i}---\left(17\right)$
式(16)所示的矩阵Ωl一般称为Gram矩阵；
由于式(13)中除αl外，其他变量只与步骤1式(3)中基函数空间和系统输出有关，
均为已知量，因此求解式(13)所示的线性方程，即能够得到αl的值；将αl的值代入式(14)中，
即能够得到aij,lm的值；再将aij,lm代入式(5)中，能够得到步骤1式(2)中矢量时变自回归模
型(VTAR)系数矩阵Ai[k]；采用时间冻结法，据系数矩阵Ai[k]求第k个时刻点系统的频率和
阻尼，完成线性时变结构模态参数的辨识。
2.如权利要求1所述的一种仅输出线性时变结构模态参数辨识方法，其特征在于：所述
的步骤3具体实现方法包括如下步骤：
步骤3.1：基于Gamma测试的非参数方法确定正则因子γ；
步骤3.1.1：寻找最近的邻近点，所述的最近的邻近点指与特定点距离范数最小的数据
点；
离数据点v[i]最近的第1个数据点编号为：
nn(i,1)＝argmin1≤N,j≠i||v[i]-v[j]|| (18)
同理，定义离数据点v[i]最近的第k个数据点编号为：
nn(i,k)＝argmin1≤j≤N,j≠i,nn(i,1),...,nn(i,κ-1)||v[i]-v[j]|| (19)
即为移除前k-1个最近点之后的最近点，因此离数据点v[i]最近的第k个数据点为v[nn
(i,κ)]；
步骤3.1.2：求Delta测试和Gamma测试，用于步骤3.1.3中正则因子γ的确定；
定义式(4)所示VTAR模型中的输入输出分别如式(20)中z[i]和yl[i]所示：
$z\[i\]=\left[\begin{array}{c}x\[k-1\]\\ x\[k-2\]\\ \·\\ \·\\ \·\\ x\[k-n_a\]\end{array}\right],y_l\[i\]=x_l\[k\]---\left(20\right)$
Delta测试定义如式(21)所示：
${\mathrm{\Δ}}_{l,N,\mathrm{\κ}}=\frac{1}{2N}\underset{i=1}{\overset{N-n_a}{\Σ}}{(y_l\[i\]-y_l\[nn(i,k)\])}^2---\left(21\right)$
式中，nn(i,k)为离数据点z[i]的第k个最近数据点，yl[nn(i,κ)]是数据点z[nn(i,k)]
对应的输出；
Gamma测试定义如式(22)所示：
${\mathrm{\Γ}}_{N,\mathrm{\ρ},\mathrm{\κ}}=\frac{1}{N}\underset{i=1}{\overset{N-n_a}{\Σ}}\left|\right|z\[i\]-z\[nn(i,\mathrm{\κ})\]|{|}^{\mathrm{\ρ}}---\left(22\right)$
式中，nn(i,k)为离数据点z[i]的第k个最近数据点，ρ为实数；
由Gamma测试能够估计输出的噪声如式(23)所示：
${\mathrm{\σ}}_{l,e}^2\≈E\left({\mathrm{\Δ}}_{l,N,b}\right)-E\left({\mathrm{\Γ}}_{N,2,b}\right)\frac{\underset{b=1}{\overset{\mathrm{\κ}}{\Σ}}({\mathrm{\Γ}}_{N,2,b}-E({\mathrm{\Γ}}_{N,2,b})\left)({\mathrm{\Δ}}_{l,N,b}-E({\mathrm{\Δ}}_{l,N,b})\right)}{\underset{b=1}{\overset{\mathrm{\κ}}{\Σ}}{({\mathrm{\Γ}}_{N,2,b}-E({\mathrm{\Γ}}_{N,2,b}\left)\right)}^2}---\left(23\right)$
式中，
步骤3.1.3：利用步骤3.1.2式(23)中及输出变量求步骤1公式(6)中正则因子γ；
正则因子γ定义为输出变量与输出噪声变量的比，表达式如式(24)所示：
${\mathrm{\γ}}_l=\frac{{\mathrm{\σ}}_{l,y}^2}{{\mathrm{\σ}}_{l,e}^2}---\left(24\right)$
式中，
步骤3.2：基于实际经验给出步骤1式(6)中基函数宽度减缩系数δ；
基函数宽度缩减因子δ根据工程经验确定，取值优选2到6。
3.如权利要求1或2所述的一种仅输出线性时变结构模态参数辨识方法，其特征在于：
所述的步骤4具体实现方法包括如下步骤，
步骤4.1：确定步骤1式(6)中最小二乘支持向量机矢量时变自回归模型(LS-SVM-VTAR)
阶数na；
最小二乘支持向量机矢量时变自回归模型(LS-SVM-VTAR)阶数na的确定依据贝叶斯信
息量准则(BayesianInformationCriterion，BIC)和赤池信息量准则(Akaike Information
Criterion，AIC)；当最小二乘支持向量机矢量时变自回归模型(LS-SVM-VTAR)阶数达到最
优时，AIC和BIC的值取最小；因此，赤池信息量准则AIC和贝叶斯信息量准则BIC值取最小时
对应的最小二乘支持向量机矢量时变自回归模型(LS-SVM-VTAR)阶数即为na值；
步骤4.2：确定步骤1式(6)中函数空间阶数pa；
函数空间阶数pa的确定依赖于残差平方和(residual sum of squares，RSS)与序列平
方和(series sum of squares，SSS)的比值，即RSS/SSS；当函数空间阶数pa最优时，RSS/
SSS取最小值；因此，RSS/SSS取最小值时对应的函数空间阶数即为pa值。
4.如权利要求3所述的一种仅输出线性时变结构模态参数辨识方法，其特征在于：还包
括应用步骤5辨识的结构模态参数指导结构动力学领域结构分析与设计的步骤6。
根据步骤5辨识的结构模态参数，能得到工程结构的振动频率范围，检测是否符合工程
结构标准及隔振等要求，指导工程结构的设计；另外，得到的结构模态参数还能为时变结构
的健康监测、结构故障诊断、结构振动控制等方面的应用提供支持。
5.一种仅输出线性时变结构模态参数辨识方法，其特征在于：首先结合支持向量机、最
小二乘方法和矢量时变自回归模型，将矢量时变自回归模型(VTAR)的系数投影到由径向基
函数序列表示的函数空间中，推导出最小二乘支持向量机矢量时变自回归(LS-SVM-VTAR)
模型的代价函数；利用Wendland紧支径向基函数构造函数空间使函数空间变得稀
疏；基于Gamma测试的非参数方法确定正则因子，基于实际经验给出基函数宽度减缩系数；
依据贝叶斯信息量准则(BayesianInformationCriterion，BIC)和赤池信息量准则(Akaike
Information Criterion，AIC)确定时变自回归模型(VTAR)阶数；依据残差平方和
(residual sum of squares，RSS)与序列平方和(series sum of squares，SSS)的比值，即
RSS/SSS确定函数空间阶数；最后根据代价函数求最小二乘支持向量机矢量时变自回归模
型(LS-SVM-VTAR)系数矩阵表达式，并根据时间冻结方法求系统的模态频率，完成线性时变
结构模态参数的辨识。
6.如权利要求5所述的一种仅输出线性时变结构模态参数辨识方法，其特征在于：得到
的结构模态参数能得到工程结构的振动频率范围，检测是否符合工程结构标准及隔振等要
求，指导工程结构的设计；另外，得到的结构模态参数还能为时变结构的健康监测、结构故
障诊断、结构振动控制等方面的应用提供支持。

一种仅输出线性时变结构模态参数辨识方法

技术领域

本发明涉及一种仅输出线性时变结构模态参数辨识方法，尤其涉及一种基于矢量
时变自回归模型和最小二乘支持向量机的仅输出线性时变结构模态参数辨识方法，属于结
构动力学技术领域。

背景技术

在日常生活和工业生产中，由于工作条件的变化、结构的老化、正常磨损等不同原
因，许多工程结构的模态参数都会随时间发生变化，表现出时变特征。例如行驶车辆激励下
的车桥系统、燃料质量随时间变化的运载火箭、飞行过程中附加气动效应下的飞机，以及可
展开的几何可变航天机构等均展现出时变特性。随着生活生产的需要和科学技术的发展，
逐渐要求对时变工程结构进行设计、故障监测、振动控制等，这就需要对时变结构有清晰的
认识，掌握其变化规律，发展对时变结构的分析方法。

因此，作为时变结构动力学特性分析的重要方法和途径，时变结构模态参数辨识
研究将成为未来结构动力学研究的重点之一。时变结构模态参数辨识可以辨识时变结构的
模态频率、模态振型和模态阻尼，这些参数具有重要的物理意义，可以为时变结构的结构设
计、结构健康监测、结构故障诊断、结构振动控制等方面的应用提供有力的支持。

根据分析域的不同，目前针对时变结构的模态参数辨识方法主要分为两类：

第一类是时域方法，主要分为两类：基于状态空间模型的方法和基于ARMA模型的
方法。

在基于状态空间模型的方法中，Liu通过对一系列由输出响应和附加输入数据组
成的Hankel矩阵进行奇异值分解，得到线性时变系统的模态参数。Liu和Deng降低了状态空
间方法对噪声的敏感度，并将此方法应用到移动悬臂梁的实验中。

在基于ARMA模型的方法中，Petsounis和Fassois提出一种时变的ARMA模型，并用
于非平稳随机振动的建模。Poulimenos和Fassois调研并比较研究了多种基于TARMA模型的
非平稳随机振动的建模方法，包括非结构化参数演化、随机参数演化和确定性参数演化。
Spiridonakos和Fassois提出一种自适应的函数序列TARMA方法，采用B样条函数求解TARMA
模型的系数，并用于非平稳振动的建模中。Yang等提出一种基于移动Kriging型函数的矢量
TARMA模型，能够较好地处理突变问题。

第二类是时频域方法，主要分为两类：非参数化方法和参数化方法。

非参数化方法直接采用时域分析，例如小波变换、平滑伪Wigner-Ville分布、
Hilbert-Huang变换等，不依赖任何系统的参数化模型。Ghanem和Francesco提出一种基于
小波的辨识方法，通过解小波展开式方程辨识出模态参数。Roshan-Ghias等采用平滑伪
Wigner-Ville分布方法，辨识系统的自然频率和阻尼系数。Xu提出一种基于响应信号Gabor
展开的模态参数辨识方法。

参数化方法利用系统时频域的参数化模型表征系统并辨识出参数化模型的参数，
例如时变公分母模型。Zhou等提出一种两步最小二乘方法辨识时变结构的频率和模态振
型。Zhou等进一步基于时频域最大似然估计辨识时变结构的模态参数，该方法对模态阶数
不敏感，且能选择频率带宽。Louarroudi等基于噪声输入输出监督针对周期时变系统提出
一种模态参数辨识方法。将时变系统频域模态辨识方法引入气动弹性颤振的分析和预测
中。

目前，提高时变系统辨识效果和效率的方法主要针对两个方面：参数建模和估计
器。前者主要致力于系统本身和时变参数的建模，例如正交多项式、B样条、Kriging型函数
等。后者关注估计器的形式，例如最小二乘、最大似然、贝叶斯估计等，能提供估计参数模型
中参数的方法。

过去二十年，支持向量机技术(SVM)在分类和函数估计中得到了广泛的关注。然
而，支持向量机的解是一个约束凸二次规划问题，当数据规模很大时计算成本很高。
Suykens采用最小二乘修正最初的支持向量机，形成LS-SVM方法。LS-SVM方法考虑最初支持
向量机中的等式约束，并用误差平方和代价函数代替Vapnikε不灵敏损失函数，将低效率的
约束二次规划问题转化为线性问题。LS-SVM及其修正方法广泛应用于辨识、函数估计和预
测中，如气体预测、非线性系统辨识、人类活动识别及交通流预测等。

总之，由于支持向量机技术的解是一个约束凸二次规划问题，在实际求解过程中
计算量巨大，计算效率低，计算成本高昂。而且对模型阶数和样本敏感度高，而且对模型阶
数和正则因子等参数的选择均采用交叉验证的方法，没有一种简单易行的方法，使用者在
应用时很复杂。

发明内容

针对基于支持向量机的时变系统辨识方法存在的上述技术问题，本发明公开的一
种仅输出线性时变结构模态参数辨识方法要解决的技术问题是，提高基于支持向量机的线
性时变系统辨识方法的计算效率、降低计算成本。此外，降低对矢量时变自回归模型(VTAR)
阶数、紧支径向基函数空间阶数等参数的敏感度，使得该方法计算量小、鲁棒性强，方便使
用。本发明即使在缺乏专业知识背景的情况下也能进行操作，能够在结构动力学工程应用
中广泛应用于线性时变结构的模态辨识。

本发明的目的是通过下述技术方案实现的：

本发明公开的一种仅输出线性时变结构模态参数辨识方法，首先结合支持向量
机、最小二乘方法和矢量时变自回归模型，将矢量时变自回归模型(VTAR)的系数投影到由
径向基函数序列表示的函数空间中，推导出最小二乘支持向量机矢量时变自回归(LS-SVM-
VTAR)模型的代价函数；利用Wendland紧支径向基函数构造函数空间使函数空间变
得稀疏；基于Gamma测试的非参数方法确定正则因子，基于实际经验给出基函数宽度减缩系
数；依据贝叶斯信息量准则(BayesianInformationCriterion，BIC)和赤池信息量准则
(Akaike Information Criterion，AIC)确定时变自回归模型(VTAR)阶数；依据残差平方和
(residual sum of squares，RSS)与序列平方和(series sum of squares，SSS)的比值，即
RSS/SSS确定函数空间阶数；最后根据代价函数求最小二乘支持向量机矢量时变自回归模
型(LS-SVM-VTAR)系数矩阵表达式，并根据时间冻结方法求系统的模态频率，完成线性时变
结构模态参数的辨识。

根据所述的一种仅输出线性时变结构模态参数辨识方法得到的结构模态参数，能
得到工程结构的振动频率范围，检测是否符合工程结构标准及隔振等要求，能够指导工程
结构的设计。另外，得到的结构模态参数还能为时变结构的健康监测、结构故障诊断、结构
振动控制等方面的应用提供有力的支持，具有广泛的应用前景与效益。

本发明公开的一种仅输出线性时变结构模态参数辨识方法，包括以下步骤：

步骤1：推导出最小二乘支持向量机矢量时变自回归模型(LS-SVM-VTAR)的代价函
数Lp(aij,lm,αl[k],el[k])，具体包括如下步骤：

步骤1.1：推导最小二乘支持向量机模型(LS-SVM)的代价函数L(w,b,e,α)如式(1)
所示：

式中，L(w,b,e,α)表示代价函数，w为参数矢量，b为实数，e＝[e1 e2…eN]T为误差
向量，α＝[α1 α2…αN]T为拉格朗日乘子向量，γ为正则因子，为支持向量机训练样
本数据，为N维至更高nh维的映射，上标T表示转置运算，N为训练样本点个
数。

步骤1.2：将矢量时变自回归模型(VTAR)的系数投影到由径向基函数序列表示的
函数空间中。

矢量时变自回归模型(VTAR)如式(2)所示：

$x\[k\]=-\underset{i=1}{\overset{n_a}{\Σ}}{A}_i\[k\]x\[k-i\]+e\[k\]---\left(2\right)$

式中，为系统N0个通道的输出向量，k为第k个时刻点，na为VTAR模型阶
数，e[k]为k时刻的误差或不可观测的非平稳扰动，Ai[k]为第i阶与时刻相关的VTAR系数矩
阵。

将矢量时变自回归模型(VTAR)的系数投影到如式(3)所示的函数空间中：

式中，pa为函数空间阶数，fj(j＝1,2,…pa)为第j阶函数向量。将Ai[k]展开，即得
到函数空间表示的矢量时变自回归模型(VTAR)，如式(4)所示：

$x\[k\]=-\underset{i=1}{\overset{n_a}{\Σ}}\underset{j=1}{\overset{p_a}{\Σ}}{A}_{ij}{f}_j\[k\]x\[k-i\]+e\[k\]---\left(4\right)$

x[k]的分量形式如式(5)所示：

$x_l\[k\]=-\underset{i=1}{\overset{n_a}{\Σ}}\underset{j=1}{\overset{p_a}{\Σ}}{f}_j\[k\]\underset{m=1}{\overset{N_o}{\Σ}}{a}_{ij,lm}{x}_m\[k-i\]+e_l\[k\]---\left(5\right)$

式中，aij，lm为矩阵Aij的元素，l表示第l个输出通道。

步骤1.3：得出最小二乘支持向量机矢量时变自回归模型(LS-SVM-VTAR)的代价函
数Lp(aij,lm,αl[k],el[k])。

结合式(1)和式(5)得到所需的代价函数Lp(aij,lm,αl[k],el[k])如式(6)所示：

$\begin{array}{c}{L}_p(a_{ij,lm},{\mathrm{\α}}_l\[k\],e_l\[k\])=\frac{1}{2}\underset{i=1}{\overset{n_a}{\Σ}}\underset{j=1}{\overset{p_a}{\Σ}}\underset{l=1}{\overset{N_o}{\Σ}}\underset{m=1}{\overset{N_o}{\Σ}}{a}_{ij,lm}^2+\mathrm{\γ}\frac{1}{2}\underset{k=1}{\overset{N}{\Σ}}\underset{l=1}{\overset{N_o}{\Σ}}{e}_l^2\[k\]\\ -\underset{k=1}{\overset{N}{\Σ}}\underset{l=1}{\overset{N_o}{\Σ}}{\mathrm{\α}}_l\[k\]\left(-\underset{i=1}{\overset{n_a}{\Σ}}\underset{j=1}{\overset{p_a}{\Σ}}{f}_j\[k\]\underset{m=1}{\overset{N_o}{\Σ}}{a}_{ij,lm}{x}_m\[k-i\]+e_l\[k\]-x_l\[k\]\right)\end{array}---\left(6\right)$

步骤2：利用Wendland紧支径向基函数构造步骤1中式(3)所示的函数空间使
函数空间变得稀疏，提高计算效率。

采用的Wendland紧支径向基函数f[k]具体形式如式(7)所示：

式中，Nw是给定的Wendland紧支径向基函数序列的非零部分，且
δ是径向基函数序列的无量纲宽度缩减因子。是φ1,λ[r]的离散形式，φ1λ[r]表达式
如式(8)所示：

${\mathrm{\φ}}_{1,\mathrm{\λ}}\left(r\right)={(1-r)}_{+}^{2\mathrm{\λ}+1}{p}_{\mathrm{\λ},\mathrm{\λ}+1}\left(r\right)---\left(8\right)$

λ是多项式的次数，pλ,λ+1(r)是λ次的多项式，运算(a)+＝max{a,0}。当λ取值分别为
0,1,2,3时，pλ,λ+1(r)的表达式分别如式(9)至式(12)所示。

φ1,λ(r)＝(1-r)+ (9)

${\mathrm{\φ}}_{1,\mathrm{\λ}}\left(r\right)={(1-r)}_{+}^3(1+3r)---\left(10\right)$

${\mathrm{\φ}}_{1,\mathrm{\λ}}\left(r\right)={(1-r)}_{+}^5(3+15r+24r^2)/3---\left(11\right)$

${\mathrm{\φ}}_{1,\mathrm{\λ}}\left(r\right)={(1-r)}_{+}^7(15+105r+285r^2+315r^3)/15---\left(12\right)$

式(7)所示的fj[k]长度仅为2Nwδ，其余元素均为0。式(3)中的函数空间由大量
的fj[k]组成，因此，利用Wendland紧支径向基函数构造步骤1中式(3)所示的函数空间
能使函数空间变得稀疏，提高计算效率。

步骤3：基于Gamma测试的非参数方法确定步骤1中公式(6)中正则因子γ，基于工
程经验给出步骤1中公式(6)中基函数宽度减缩系数δ。具体实现步骤如下：

步骤3.1：基于Gamma测试的非参数方法确定正则因子γ。

步骤3.1.1：寻找最近的邻近点，所述的最近的邻近点指与特定点距离范数最小的
数据点。

离数据点v[i]最近的第1个数据点编号为：

nn(i,1)＝argmin1≤N,j≠i||v[i]-v[j]|| (13)

同理，定义离数据点v[i]最近的第k个数据点编号为：

nn(i,k)＝argmin1≤j≤N,j≠i,nn(i,1),...,nn(i,κ-1)||v[i]-v[j]|| (14)

即为移除前k-1个最近点之后的最近点，因此离数据点v[i]最近的第k个数据点为
v[nn(i,κ)]。

步骤3.1.2：求Delta测试和Gamma测试，用于步骤3.1.3中正则因子γ的确定。

定义式(4)所示VTAR模型中的输入输出分别如式(15)中z[i]和yl[i]所示：

$z\[i\]=\left[\begin{array}{c}x\[k-1\]\\ x\[k-2\]\\ .\\ .\\ .\\ x\[k-n_a\]\end{array}\right],y_l\[i\]=x_l\[k\]---\left(15\right)$

Delta测试定义如式(16)所示：

${\mathrm{\Δ}}_{l,N,\mathrm{\κ}}=\frac{1}{2N}\underset{i=1}{\overset{N-n_a}{\Σ}}{(y_l\[i\]-y_l\[nn(i,k)\])}^2---\left(16\right)$

式中，nn(i,k)为离数据点z[i]的第k个最近数据点，yl[nn(i,κ)]是数据点z[nn
(i,k)]对应的输出。

Gamma测试定义如式(17)所示：

${\mathrm{\Γ}}_{N,\mathrm{\ρ},\mathrm{\κ}}=\frac{1}{N}\underset{i=1}{\overset{N-n_a}{\Σ}}\left|\right|z\[i\]-z\[nn(i,\mathrm{\κ})\]|{|}^{\mathrm{\ρ}}---\left(17\right)$

式中，nn(i,k)为离数据点z[i]的第k个最近数据点，ρ为实数。

由Gamma测试能够估计输出的噪声如式(18)所示：

${\mathrm{\σ}}_{l,e}^2\≈E\left({\mathrm{\Δ}}_{l,N,b}\right)-E\left({\mathrm{\Γ}}_{N,2,b}\right)\frac{\underset{b=1}{\overset{\mathrm{\κ}}{\Σ}}({\mathrm{\Γ}}_{N,2,b}-E({\mathrm{\Γ}}_{N,2,b}\left)\right)({\mathrm{\Δ}}_{l,N,b}-E({\mathrm{\Δ}}_{l,N,b}\left)\right)}{\underset{b=1}{\overset{\mathrm{\κ}}{\Σ}}{({\mathrm{\Γ}}_{N,2,b}-E({\mathrm{\Γ}}_{N,2,b}\left)\right)}^2}---\left(18\right)$

式中，

步骤3.1.3：利用步骤3.1.2式(18)中及输出变量求步骤1公式(6)中正则因
子γ。

正则因子γ定义为输出变量与输出噪声变量的比，表达式如式(19)所示：

${\mathrm{\γ}}_l=\frac{{\mathrm{\σ}}_{l,y}^2}{{\mathrm{\σ}}_{l,e}^2}---\left(19\right)$

式中，

步骤3.2：基于实际经验给出步骤1式(6)中基函数宽度减缩系数δ。

基函数宽度缩减因子δ根据工程经验确定，对于大多数问题，取值在2到6时效果良
好。

已有技术中确定正则因子和基函数宽度缩减因子均根据交叉验证的方法确定，计
算量很大，计算成本高，在实际工程应用中十分不便。步骤3基于Gamma测试的非参数方法确
定正则因子γ，基于实际经验给出基函数宽度减缩系数δ，即直接给出正则因子和基函数宽
度减缩系数的确定方法，无需进行交叉验证，大大较少计算量，计算成本低，使用者在使用
过程中能够方便快速地确定正则因子和基函数宽度减缩系数，十分方便。

步骤4：确定步骤1式(6)中最小二乘支持向量机矢量时变自回归模型(LS-SVM-
VTAR)阶数na和函数空间阶数pa。包括如下步骤：

步骤4.1：确定步骤1式(6)中最小二乘支持向量机矢量时变自回归模型(LS-SVM-
VTAR)阶数na。

最小二乘支持向量机矢量时变自回归模型(LS-SVM-VTAR)阶数na的确定依据贝叶
斯信息量准则(BayesianInformationCriterion，BIC)和赤池信息量准则(Akaike
Information Criterion，AIC)。当最小二乘支持向量机矢量时变自回归模型(LS-SVM-
VTAR)阶数达到最优时，AIC和BIC的值取最小。因此，赤池信息量准则AIC和贝叶斯信息量准
则BIC值取最小时对应的最小二乘支持向量机矢量时变自回归模型(LS-SVM-VTAR)阶数即
为na值。

步骤4.2：确定步骤1式(6)中函数空间阶数pa。

函数空间阶数pa的确定依赖于残差平方和(residual sum of squares，RSS)与序
列平方和(series sum of squares，SSS)的比值，即RSS/SSS。当函数空间阶数pa最优时，
RSS/SSS取最小值。因此，RSS/SSS取最小值时对应的函数空间阶数即为pa值。

已有技术中，在最小二乘支持向量机矢量时变自回归模型(LS-SVM-VTAR)阶数na
达到最优时，AIC和BIC的值并不能取到最小值；函数空间阶数pa达到最优时，RSS/SSS的值
也不能取到最小值。因此，在已有技术中，根据AIC和BIC的值并不能取到最优的最小二乘支
持向量机矢量时变自回归模型(LS-SVM-VTAR)阶数na，根据RSS/SSS的值并不能取到最优的
函数空间阶数pa，使得AIC、BIC与RSS/SSS值失去指导参数选择的意义，同时使得应用变得
困难。

步骤5：根据步骤1中式(6)所示的代价函数Lp(aij,lm,αl[k],el[k])求步骤1式(2)中
矢量时变自回归模型(VTAR)系数矩阵Ai[k]表达式，并根据时间冻结方法求系统的模态频
率，完成线性时变结构模态参数的辨识。

将式(6)中的代价函数Lp(aij,lm,αl[k],el[k])对各个自变量aij,lm,αl[k],el[k]求
导，并整理得：

$({\mathrm{\Ω}}_l+\frac{I}{{\mathrm{\γ}}_l}){\mathrm{\α}}_l={\tilde{x}}_l,(l=1,2,...N_0)---\left(20\right)$

$a_{ij,lm}=-\underset{k=1}{\overset{N}{\Σ}}{\mathrm{\α}}_l\[k\]f_j\[k\]x_m\[k-i\]---\left(21\right)$

式中，γl是第l个输出通道对应的正则因子，I为N维单位矩阵，αl、xl、Ωl分别如式
(22)、(23)所示：

${\mathrm{\α}}_l=\left[\begin{array}{c}{\mathrm{\α}}_l\left(1\right)\\ {\mathrm{\α}}_l\left(2\right)\\ .\\ .\\ .\\ {\mathrm{\α}}_l\left(N\right)\end{array}\right],{\tilde{x}}_l=\left[\begin{array}{c}{x}_l\left(1\right)\\ x_l\left(2\right)\\ .\\ .\\ .\\ x_l\left(N\right)\end{array}\right]---\left(22\right)$

其中，⊙表示对应元素相乘运算，F矩阵和Rm矩阵元素如式(24)所示：

$F_{ks}=\underset{j=1}{\overset{p_a}{\Σ}}{f}_{j,k}{f}_{j,s},R_{m,ks}=\underset{i=1}{\overset{n_a}{\Σ}}{x}_{m,k-i}{x}_{m,s-i}---\left(24\right)$

式(23)所示的矩阵Ωl一般称为Gram矩阵。

由于式(20)中除αl外，其他变量只与步骤1式(3)中基函数空间和系统输出
有关，均为已知量，因此求解式(20)所示的线性方程，即能够得到αl的值。将αl的值代入式
(21)中，即能够得到aij，lm的值。再将aij，lm代入式(5)中，能够得到步骤1式(2)中矢量时变自
回归模型(VTAR)系数矩阵Ai[k]。采用时间冻结法，据系数矩阵Ai[k]求第k个时刻点系统的
频率和阻尼，完成线性时变结构模态参数的辨识。

步骤6：应用步骤5辨识的结构模态参数指导结构动力学领域的结构分析与设计。

根据步骤5辨识的结构模态参数，能得到工程结构的振动频率范围，检测是否符合
工程结构标准及隔振等要求，能够指导工程结构的设计。另外，得到的结构模态参数还能为
时变结构的健康监测、结构故障诊断、结构振动控制等方面的应用提供有力的支持，具有广
泛的应用前景与效益。

有益效果：

1、本发明公开的一种仅输出线性时变结构模态参数辨识方法，采用Wendland紧支
径向基函数序列空间表示矢量时变自回归模型(VTAR)的系数矩阵，使得Gram矩阵变得稀
疏，大大降低运算成本；

2、本发明公开的一种仅输出线性时变结构模态参数辨识方法，采用基于Gamma测
试的非参数方法确定正则因子γ，基于工程经验给出基函数宽度减缩系数δ，相对于已有技
术中交叉验证方法使用更方便、效率更高；

3、本发明公开的一种仅输出线性时变结构模态参数辨识方法，根据AIC、BIC最小
值确定VTAR模型的阶数，根据RSS/SSS最小值确定函数空间阶数，相对于已有技术中的方法
使用更方便；

4、本发明公开的一种仅输出线性时变结构模态参数辨识方法，在矢量时变自回归
模型和最小二乘等已有技术中，结合支持向量机方法。在支持向量机技术中，少量支持向量
决定了最终结果，能够抓住关键样本，剔除大量冗余样本，使得本发明对矢量时变自回归模
型(VTAR)阶数、函数空间阶数及样本大小不敏感，且对过估计问题具有很好的适应性。因
此，本发明的鲁棒性很强，在工程应用中，对参数的选择不敏感，即使在缺乏专业知识背景
的情况下也能进行操作。

附图说明：

图1为本发明一种仅输出线性时变结构模态参数辨识方法的流程图；

图2为具体实施方式中的三自由度弹簧-阻尼器-质量系统；

图3为具体实施方式中的三自由度时变系统的响应曲线。其中，图3(A)、图3(B)、图
3(C)分别表示图2中质量块m1、m2、m3的响应曲线；

图4为具体实施方式中采用Wendland紧支径向基函数构造的函数空间计算稀
疏Gram矩阵对应的计算机CPU计算时间与已有技术中计算全阶Gram矩阵方法CPU计算时间
比较图；

图5为具体实施方式中采用基于Gamma测试的非参数方法确定正则因子γ所需的
CPU计算时间与已有技术中的交叉验证方法确定正则因子γ所需的CPU计算时间对比图；

图6为具体实施方式中采用LS-SVM-VTAR模型得到的AIC与BIC随LS-SVM-VTAR模型
阶数的变化曲线。其中圆心实点表示AIC随LS-SVM-VTAR模型阶数na的变化曲线，三角形点
表示BIC随LS-SVM-VTAR模型阶数na的变化曲线；

图7为具体实施方式中采用LS-SVM-VTAR模型得到的RSS/SSS值随函数空间阶数pa
的变化曲线；

图8为具体实施方式中采用LS-SVM-VTAR方法与采用已有技术中的最小二乘方法
(LS)求得的时变系统模态频率随时间变化曲线，图8(a)表示采用LS-SVM-VTAR方法得到的
结果，图8(b)表示采用已有技术中的最小二乘方法(LS)得到的结果。其中，黑色空心点表示
100次蒙特卡洛模拟平均值，粉红色三角点表示一次蒙特卡洛模拟结果，红色圈点表示基准
值。

具体实施方式

为了更好地说明本发明的目的和优点，下面通过对一个随机激励下的三自由度时
变结构进行动力学分析，对本发明做出详细解释。

实施例1：

本实施例的三自由度弹簧-阻尼器-质量系统，如图2所示。本实施例的三自由度弹
簧-阻尼器-质量系统包含三个质量块m1、m2、m3，四个阻尼器c1、c2、c3、c4，四个弹簧k1(t)、k2
(t)、k3(t)、k4(t)，其中，三个质量块与四个阻尼器均是定常的，不随时间变化，而四个弹簧
的刚度随时间变化。

三自由度时变系统中弹簧刚度ki(t)随时间变化关系如式(25)所示：

$k_i\left(t\right)=k_{i,0}+k_{i,1}cos\left(\frac{2\mathrm{\π}t}{p_{i,1}}\right)+k_{i,2}sin\left(\frac{2\mathrm{\π}t}{p_{i,2}}\right),(i=1,2,3,4)---\left(25\right)$

三自由度系统的各个变量取值如表1所示。

表1三自由度时变系统参数

图2所示的三自由度系统的动力学控制方程如式(26)所示：

$M\stackrel{\·\·}{x}+C\stackrel{\·}{x}+K\left(t\right)x=f\left(t\right)---\left(26\right)$

式中，M、C分别为质量矩阵和阻尼矩阵，K(t)为时变的刚度矩阵，x为系统的响应，f
(t)为外界激励。M、C、K(t)如式(27)所示：

$\begin{array}{c}M=diag(10,15,25),C=\left[\begin{array}{ccc}6& -2& 0\\ -2& 4& -2\\ 0& -2& 5\end{array}\right],\\ K=\left[\begin{array}{ccc}{k}_1\left(t\right)+k_2\left(t\right)& -k_2\left(t\right)& 0\\ -k_2\left(t\right)& k_2\left(t\right)+k_3\left(t\right)& -k_3\left(t\right)\\ 0& -k_3\left(t\right)& k_3\left(t\right)+k_4\left(t\right)\end{array}\right]\end{array}---\left(27\right)$

分别对三个质量块施加Gauss白噪声激励，得到系统中三个质量块的响应。系统响
应采用变步长四阶龙格库塔法求解。采用三个自由度位移为辨识所用的响应信号样本，将
龙格库塔法的求解结果以f＝16Hz重新采样，记录时间为500s(t∈[0,500]),信号长度N＝
8000，三个自由度的位移响应曲线如图3所示。

本实施例公开的一种仅输出线性时变结构模态参数辨识方法，包括以下步骤：

步骤1：推导出最小二乘支持向量机矢量时变自回归模型(LS-SVM-VTAR)的代价函
数Lp(aij,lm，αl[k]，el[k])。

对此具体实施方式，在式(6)中，输出通道数N0＝3，样本长度N＝8000，因此得到具
体实施方式的最小二乘支持向量机矢量时变自回归模型(LS-SVM-VTAR)的代价函数Lp
(aij，lm，αl[k]，el[k])如式(28)所示。

$\begin{array}{c}{L}_p(a_{ij,lm},{\mathrm{\α}}_l\[k\],e_l\[k\])=\frac{1}{2}\underset{i=1}{\overset{n_a}{\Σ}}\underset{j=1}{\overset{p_a}{\Σ}}\underset{l=1}{\overset{3}{\Σ}}\underset{m=1}{\overset{3}{\Σ}}{a}_{ij,lm}^2+\mathrm{\γ}\frac{1}{2}\underset{k=1}{\overset{8000}{\Σ}}\underset{l=1}{\overset{3}{\Σ}}{e}_l^2\[k\]\\ -\underset{k=1}{\overset{8000}{\Σ}}\underset{l=1}{\overset{3}{\Σ}}{\mathrm{\α}}_l\[k\]\left(-\underset{i=1}{\overset{n_a}{\Σ}}\underset{j=1}{\overset{p_a}{\Σ}}{f}_j\[k\]\underset{m=1}{\overset{3}{\Σ}}{a}_{ij,lm}{x}_m\[k-i\]+e_l\[k\]-x_l\[k\]\right)\end{array}---\left(28\right)$

步骤2：利用Wendland紧支径向基函数构造式(3)所示的函数空间使Gram矩阵
变得稀疏，提高计算效率。

在本具体实施方式中，函数空间中元素fj[k]选取Wendland紧支径向基函数，
如式(29)所示。

式(29)所示的fj[k]长度仅为2Nwδ，其余元素均为0。在计算式(23)所示的Gram矩阵
Ωl时，可以大大减少计算量。例如，当取函数空间阶数pa＝30时，采用Wendland紧支径
向基函数构造的函数空间计算稀疏Gram矩阵对应的计算机CPU计算时间与已有技术中
计算全阶Gram矩阵方法CPU计算时间比较如图3所示。由图4可以看出，采用Wendland紧支径
向基函数构造式(3)所示的函数空间使得Gram矩阵变得稀疏，大大提高计算效率。

步骤3：基于Gamma测试的非参数方法确定式(28)中正则因子γ，基于工程经验给
出式(28)中基函数宽度减缩系数δ。

根据式(19)求得三自由度三个输出通道的正则因子γ值，基函数宽度减缩系数δ
根据工程经验取为3。

在已有技术中的方法中，正则因子γ的取值均基于交叉验证的方法，计算量巨大，
计算成本很高。采用基于Gamma测试的非参数方法确定正则因子γ所需的CPU计算时间与已
有技术中的交叉验证方法确定正则因子γ所需的CPU计算时间对比如图5所示。由图5可以
看出，采用基于Gamma测试的非参数方法确定正则因子γ与交叉验证方法相比，计算量大大
减少。

步骤4：确定式(28)中最小二乘支持向量机矢量时变自回归模型(LS-SVM-VTAR)阶
数na和函数空间阶数pa。

最小二乘支持向量机矢量时变自回归模型(LS-SVM-VTAR)阶数na的确定依据贝叶
斯信息量准则(BayesianInformationCriterion，BIC)和赤池信息量准则(Akaike
Information Criterion，AIC)。当最小二乘支持向量机矢量时变自回归模型(LS-SVM-
VTAR)阶数达到最优时，AIC和BIC的值取最小。因此，AIC和BIC值取最小时对应的最小二乘
支持向量机矢量时变自回归模型(LS-SVM-VTAR)阶数即为na值。

在本具体实施方式中，AIC与BIC随最小二乘支持向量机矢量时变自回归模型(LS-
SVM-VTAR)阶数na的变化曲线如图6所示，由图6得到na＝2。

函数空间阶数pa的确定依赖于残差平方和(residual sum of squares，RSS)与序
列平方和(series sum of squares，SSS)的比值，即RSS/SSS。当函数空间阶数pa最优时，
RSS/SSS取最小值。因此，RSS/SSS取最小值时对应的函数空间阶数即为pa值。

在本具体实施方式中，RSS/SSS值随函数空间阶数pa的变化曲线如图7所示，由
图7得到函数空间阶数pa可取15。

步骤5：根据式(28)所示的代价函数Lp(aij，lm，αl[k]，el[k])求式(2)中矢量时变自
回归模型(VTAR)系数矩阵Ai[k]表达式，并根据时间冻结方法求系统的模态频率，完成线性
时变结构模态参数的辨识。

将式(28)中的代价函数Lp(aij，lm，αl[k]，el[k])对各个自变量aij，lm，αl[k]，el[k]求
导，并整理得：

$({\mathrm{\Ω}}_l+\frac{I}{{\mathrm{\γ}}_l}){\mathrm{\α}}_l={\tilde{x}}_l,(l=1,2,3)---\left(30\right)$

$a_{ij,lm}=-\underset{k=1}{\overset{8000}{\Σ}}{\mathrm{\α}}_l\[k\]f_j\[k\]x_m\[k-i\]---\left(31\right)$

式中，γl是第l个输出通道对应的正则因子，I为8000维单位矩阵，αl、xl、Ωl分别
如式(32)、(33)所示：

${\mathrm{\α}}_l=\left[\begin{array}{c}{\mathrm{\α}}_l\left(1\right)\\ {\mathrm{\α}}_l\left(2\right)\\ .\\ .\\ .\\ {\mathrm{\α}}_l\left(8000\right)\end{array}\right],{\tilde{x}}_l=\left[\begin{array}{c}{x}_l\left(1\right)\\ x_l\left(2\right)\\ .\\ .\\ .\\ x_l\left(8000\right)\end{array}\right]---\left(32\right)$

其中，F矩阵和Rm矩阵元素如式(34)所示：

$F_{ks}=\underset{j=1}{\overset{15}{\Σ}}{f}_{j,k}{f}_{j,s},R_{m,ks}=\underset{i=1}{\overset{2}{\Σ}}{x}_{m,k-i}{x}_{m,s-i}---\left(34\right)$

求解式(30)所示的线性方程，即能够得到αl的值。将αl的值代入式(31)中，即能够
得到aij,lm的值。再将aij，lm代入式(5)中，能够得到式(2)中矢量时变自回归模型(VTAR)系数
矩阵Ai[k]。采用时间冻结法，据系数矩阵Ai[k]求第k个时刻点系统的频率和阻尼，完成线性
时变结构模态参数的辨识。本具体实施方式得到的三自由度弹簧-阻尼器-质量系统的频率
随时间变化曲线如图8所示。

由图8能够看出，本实施例公开的一种仅输出线性时变结构模态参数辨识方法，相
对于已有技术中的最小二乘方法(LS)，能够更好地辨识出线性时变结构的频率，精度高，耗
时少，计算成本低，方便使用，在结构动力学领域具有良好的应用前景。

以上所述的具体描述，对发明的目的、技术方案和有益效果进行了进一步详细说
明，所应理解的是，以上所述仅为本发明的具体实施例，用于解释本发明，并不用于限定本
发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应
包含在本发明的保护范围之内。