CN106354695B - 一种仅输出线性时变结构模态参数辨识方法 - Google Patents
一种仅输出线性时变结构模态参数辨识方法 Download PDFInfo
- Publication number
- CN106354695B CN106354695B CN201610701225.9A CN201610701225A CN106354695B CN 106354695 B CN106354695 B CN 106354695B CN 201610701225 A CN201610701225 A CN 201610701225A CN 106354695 B CN106354695 B CN 106354695B
- Authority
- CN
- China
- Prior art keywords
- formula
- time
- function
- vtar
- varying
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Fee Related
Links
- 238000000034 method Methods 0.000 claims abstract description 102
- 239000011159 matrix material Substances 0.000 claims abstract description 31
- 238000012360 testing method Methods 0.000 claims abstract description 19
- 230000009467 reduction Effects 0.000 claims abstract description 17
- 230000008014 freezing Effects 0.000 claims abstract description 6
- 238000007710 freezing Methods 0.000 claims abstract description 6
- 238000012706 support-vector machine Methods 0.000 claims description 11
- 238000004458 analytical method Methods 0.000 claims description 7
- 238000013459 approach Methods 0.000 claims description 6
- 238000013461 design Methods 0.000 claims description 6
- 238000012544 monitoring process Methods 0.000 claims description 6
- 238000013016 damping Methods 0.000 claims description 5
- 230000001419 dependent effect Effects 0.000 claims description 5
- 230000036541 health Effects 0.000 claims description 5
- 238000003745 diagnosis Methods 0.000 claims description 4
- 238000002955 isolation Methods 0.000 claims description 4
- 238000006467 substitution reaction Methods 0.000 claims description 4
- 238000012549 training Methods 0.000 claims description 4
- 238000010276 construction Methods 0.000 claims description 3
- 238000013507 mapping Methods 0.000 claims description 2
- 230000017105 transposition Effects 0.000 claims description 2
- 238000009795 derivation Methods 0.000 claims 1
- 238000002790 cross-validation Methods 0.000 description 8
- 238000004364 calculation method Methods 0.000 description 7
- 230000008859 change Effects 0.000 description 7
- 230000004044 response Effects 0.000 description 6
- 238000005516 engineering process Methods 0.000 description 5
- 241001003127 Tarma Species 0.000 description 4
- 230000008901 benefit Effects 0.000 description 4
- 230000000694 effects Effects 0.000 description 4
- 230000005284 excitation Effects 0.000 description 4
- 241001123248 Arma Species 0.000 description 3
- 238000000342 Monte Carlo simulation Methods 0.000 description 2
- 238000012850 discrimination method Methods 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 230000009466 transformation Effects 0.000 description 2
- OYPRJOBELJOOCE-UHFFFAOYSA-N Calcium Chemical compound [Ca] OYPRJOBELJOOCE-UHFFFAOYSA-N 0.000 description 1
- 238000007476 Maximum Likelihood Methods 0.000 description 1
- 238000012952 Resampling Methods 0.000 description 1
- 230000003044 adaptive effect Effects 0.000 description 1
- 230000032683 aging Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000000052 comparative effect Effects 0.000 description 1
- 238000000354 decomposition reaction Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 230000018109 developmental process Effects 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 238000006073 displacement reaction Methods 0.000 description 1
- 238000009826 distribution Methods 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 239000000446 fuel Substances 0.000 description 1
- 238000009499 grossing Methods 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 238000009776 industrial production Methods 0.000 description 1
- 238000011835 investigation Methods 0.000 description 1
- 230000007257 malfunction Effects 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 230000007246 mechanism Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000002715 modification method Methods 0.000 description 1
- 230000035772 mutation Effects 0.000 description 1
- 230000000737 periodic effect Effects 0.000 description 1
- 230000008569 process Effects 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 238000012916 structural analysis Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Data Mining & Analysis (AREA)
- General Physics & Mathematics (AREA)
- Mathematical Optimization (AREA)
- Pure & Applied Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Mathematical Physics (AREA)
- Computational Mathematics (AREA)
- Mathematical Analysis (AREA)
- Life Sciences & Earth Sciences (AREA)
- Operations Research (AREA)
- Probability & Statistics with Applications (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Algebra (AREA)
- Evolutionary Biology (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- General Engineering & Computer Science (AREA)
- Bioinformatics & Computational Biology (AREA)
- Complex Calculations (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明公开的一种仅输出线性时变结构模态参数辨识方法，属于结构动力学技术领域。本发明首先推导出最小二乘支持向量机矢量时变自回归模型的代价函数；利用Wendland紧支径向基函数构造函数空间基于Gamma测试的非参数方法确定正则因子，基于实际经验给出基函数宽度减缩系数；依据贝叶斯信息量准则和赤池信息量准则确定时变自回归模型阶数；依据残差平方和与序列平方和的比值确定函数空间阶数；最后根据代价函数求最小二乘支持向量机矢量时变自回归模型系数矩阵表达式，并根据时间冻结方法求系统的模态频率，完成线性时变结构模态参数的辨识。本发明能够提高计算效率，增强系统鲁棒性，在结构动力学工程应用中广泛应用于线性时变结构的模态辨识。
Description
技术领域
本发明涉及一种仅输出线性时变结构模态参数辨识方法，尤其涉及一种基于矢量时变自回归模型和最小二乘支持向量机的仅输出线性时变结构模态参数辨识方法，属于结构动力学技术领域。
背景技术
在日常生活和工业生产中，由于工作条件的变化、结构的老化、正常磨损等不同原因，许多工程结构的模态参数都会随时间发生变化，表现出时变特征。例如行驶车辆激励下的车桥系统、燃料质量随时间变化的运载火箭、飞行过程中附加气动效应下的飞机，以及可展开的几何可变航天机构等均展现出时变特性。随着生活生产的需要和科学技术的发展，逐渐要求对时变工程结构进行设计、故障监测、振动控制等，这就需要对时变结构有清晰的认识，掌握其变化规律，发展对时变结构的分析方法。
因此，作为时变结构动力学特性分析的重要方法和途径，时变结构模态参数辨识研究将成为未来结构动力学研究的重点之一。时变结构模态参数辨识可以辨识时变结构的模态频率、模态振型和模态阻尼，这些参数具有重要的物理意义，可以为时变结构的结构设计、结构健康监测、结构故障诊断、结构振动控制等方面的应用提供有力的支持。
根据分析域的不同，目前针对时变结构的模态参数辨识方法主要分为两类：
第一类是时域方法，主要分为两类：基于状态空间模型的方法和基于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测试的非参数方法确定正则因子，基于实际经验给出基函数宽度减缩系数；依据贝叶斯信息量准则(Bayesian Information Criterion，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)所示：
式中，为系统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]的分量形式如式(5)所示：
式中，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)所示：
步骤2：利用Wendland紧支径向基函数构造步骤1中式(3)所示的函数空间使函数空间变得稀疏，提高计算效率。
采用的Wendland紧支径向基函数f[k]具体形式如式(7)所示：
式中，Nw是给定的Wendland紧支径向基函数序列的非零部分，且δ是径向基函数序列的无量纲宽度缩减因子。是φ1,λ[r]的离散形式，φ1λ[r]表达式如式(8)所示：
λ是多项式的次数，pλ,λ+1(r)是λ次的多项式，运算(a)+＝max{a,0}。当λ取值分别为0,1,2,3时，pλ,λ+1(r)的表达式分别如式(9)至式(12)所示。
φ1,λ(r)＝(1-r)+ (9)
式(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]所示：
Delta测试定义如式(16)所示：
式中，nn(i,k)为离数据点z[i]的第k个最近数据点，yl[nn(i,κ)]是数据点z[nn(i,k)] 对应的输出。
Gamma测试定义如式(17)所示：
式中，nn(i,k)为离数据点z[i]的第k个最近数据点，ρ为实数。
由Gamma测试能够估计输出的噪声如式(18)所示：
式中，
步骤3.1.3：利用步骤3.1.2式(18)中及输出变量求步骤1公式(6)中正则因子γ。
正则因子γ定义为输出变量与输出噪声变量的比，表达式如式(19)所示：
式中，
步骤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)和赤池信息量准则 (AkaikeInformation 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]求导，并整理得：
式中，γl是第l个输出通道对应的正则因子，I为N维单位矩阵，αl、xl、Ωl分别如式(22)、(23)所示：
其中，⊙表示对应元素相乘运算，F矩阵和Rm矩阵元素如式(24)所示：
式(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)所示：
三自由度系统的各个变量取值如表1所示。
表1三自由度时变系统参数
图2所示的三自由度系统的动力学控制方程如式(26)所示：
式中，M、C分别为质量矩阵和阻尼矩阵，K(t)为时变的刚度矩阵，x为系统的响应，f(t)为外界激励。M、C、K(t)如式(27)所示：
分别对三个质量块施加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)所示。
步骤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)和赤池信息量准则(AkaikeInformation 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]求导，并整理得：
式中，γl是第l个输出通道对应的正则因子，I为8000维单位矩阵，αl、xl、Ωl分别如式(32)、(33)所示：
其中，F矩阵和Rm矩阵元素如式(34)所示：
求解式(30)所示的线性方程，即能够得到αl的值。将αl的值代入式(31)中，即能够得到aij,lm的值。再将aij，lm代入式(5)中，能够得到式(2)中矢量时变自回归模型(VTAR)系数矩阵Ai[k]。采用时间冻结法，据系数矩阵Ai[k]求第k个时刻点系统的频率和阻尼，完成线性时变结构模态参数的辨识。本具体实施方式得到的三自由度弹簧-阻尼器-质量系统的频率随时间变化曲线如图8所示。
由图8能够看出，本实施例公开的一种仅输出线性时变结构模态参数辨识方法，相对于已有技术中的最小二乘方法(LS)，能够更好地辨识出线性时变结构的频率，精度高，耗时少，计算成本低，方便使用，在结构动力学领域具有良好的应用前景。
以上所述的具体描述，对发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例，用于解释本发明，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。
Claims (3)
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)所示：
式中，为系统N0个通道的输出向量，k为第k个时刻点，na为VTAR模型阶数，e[k]为k时刻的误差或不可观测的非平稳扰动，Ai[k]为第i阶与时刻相关的VTAR系数矩阵；
将矢量时变自回归模型(VTAR)的系数投影到如式(3)所示的函数空间中：
式中，pa为函数空间阶数，fj为第j阶函数向量，其中j的取值范围为j＝1,2,…,pa；
fj[1]到fj[N]表示fj中的N个元素；将Ai[k]展开，即得到函数空间表示的矢量时变自回归模型(VTAR)，如式(4)所示：
x[k]的分量形式如式(5)所示：
式中，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)所示：
步骤2：利用Wendland紧支径向基函数构造步骤1中式(3)所示的函数空间使函数空间变得稀疏；
采用的Wendland紧支径向基函数fj[k]具体形式如式(7)所示：
式中，Nw是给定的Wendland紧支径向基函数序列的非零部分，且其中表示对括号里的数向下取整；δ是基函数宽度减缩系数；是φ1,λ(r)的离散形式，φ1,λ(r)表达式如式(8)所示：
λ是多项式的次数，pλ,λ+1(r)是λ次的多项式，运算(a)+＝max{a,0}；当λ取值分别为0,1,2,3时，φ1,λ(r)的表达式分别如式(9)至式(12)所示；
φ1,λ(r)＝(1-r)+ (9)
步骤3：基于Gamma测试的非参数方法确定步骤1中公式(6)中正则因子γ，基于工程经验给出步骤1中公式(7)中基函数宽度减缩系数δ；
步骤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]求导，并整理得：
式中，γl是第l个输出通道对应的正则因子，I为N维单位矩阵，αl、Ωl分别如式(15)、(16)所示：
其中，⊙表示对应元素相乘运算，F矩阵和Rm矩阵元素如式(17)所示：
式(16)所示的矩阵Ωl称为Gram矩阵；
由于式(13)中除αl外，其他变量只与步骤1式(3)中基函数空间和系统输出有关，均为已知量，因此求解式(13)所示的线性方程，即能够得到αl的值；将αl的值代入式(14)中，即能够得到aij,lm的值；再将aij,lm代入式(5)中，能够得到步骤1式(2)中矢量时变自回归模型(VTAR)系数矩阵Ai[k]；采用时间冻结法，据系数矩阵Ai[k]求第k个时刻点系统的频率和阻尼，完成线性时变结构模态参数的辨识；
所述的步骤4具体实现方法包括如下步骤，
步骤4.1：确定步骤1式(6)中最小二乘支持向量机矢量时变自回归模型(LS-SVM-VTAR)阶数na；
最小二乘支持向量机矢量时变自回归模型(LS-SVM-VTAR)阶数na的确定依据贝叶斯信息量准则(Bayesian Information Criterion，BIC)和赤池信息量准则(AkaikeInformation 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值；
还包括应用步骤5辨识的结构模态参数指导结构动力学领域结构分析与设计的步骤6；
根据步骤5辨识的结构模态参数，能得到工程结构的振动频率范围，检测是否符合工程结构标准及隔振要求，指导工程结构的设计；另外，得到的结构模态参数还能为时变结构的健康监测、结构故障诊断、结构振动控制方面的应用提供支持。
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,k-1)||v[i]-v[j]|| (19)
即为移除前k-1个最近点之后的最近点，因此离数据点v[i]最近的第k个数据点为v[nn(i,k)]；
步骤3.1.2：求Delta测试和Gamma测试，用于步骤3.1.3中正则因子γ的确定；
定义式(4)所示VTAR模型中的输入输出分别如式(20)中z[i]和yl[i]所示：
Delta测试定义如式(21)所示：
式中，nn(i,k)为离数据点z[i]的第k个最近数据点，yl[nn(i,k)]是数据点z[nn(i,k)]对应的输出；
Gamma测试定义如式(22)所示：
式中，nn(i,k)为离数据点z[i]的第k个最近数据点，ρ为实数；
由Gamma测试能够估计输出的噪声如式(23)所示：
式中，
步骤3.1.3：利用步骤3.1.2式(23)中及输出变量求步骤1公式(6)中正则因子γ；
正则因子γ定义为输出变量与输出噪声变量的比，表达式如式(24)所示：
式中，
步骤3.2：基于实际经验给出步骤1式(7)中基函数宽度减缩系数δ；
基函数宽度减缩系数δ根据工程经验确定，取值为2到6。
3.一种仅输出线性时变结构模态参数辨识方法，其特征在于：首先结合支持向量机、最小二乘方法和矢量时变自回归模型，将矢量时变自回归模型(VTAR)的系数投影到由径向基函数序列表示的函数空间中，推导出最小二乘支持向量机矢量时变自回归(LS-SVM-VTAR)模型的代价函数；利用Wendland紧支径向基函数构造函数空间使函数空间变得稀疏；基于Gamma测试的非参数方法确定正则因子，基于实际经验给出基函数宽度减缩系数；依据贝叶斯信息量准则(Bayesian Information Criterion，BIC)和赤池信息量准则(Akaike Information Criterion，AIC)确定时变自回归模型(VTAR)阶数；依据残差平方和(residual sum of squares，RSS)与序列平方和(series sum of squares，SSS)的比值，即RSS/SSS确定函数空间阶数；最后根据代价函数求最小二乘支持向量机矢量时变自回归模型(LS-SVM-VTAR)系数矩阵表达式，并根据时间冻结方法求系统的模态频率，完成线性时变结构模态参数的辨识；
得到的结构模态参数能得到工程结构的振动频率范围，检测是否符合工程结构标准及隔振要求，指导工程结构的设计；另外，得到的结构模态参数还能为时变结构的健康监测、结构故障诊断、结构振动控制方面的应用提供支持。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201610701225.9A CN106354695B (zh) | 2016-08-22 | 2016-08-22 | 一种仅输出线性时变结构模态参数辨识方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201610701225.9A CN106354695B (zh) | 2016-08-22 | 2016-08-22 | 一种仅输出线性时变结构模态参数辨识方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN106354695A CN106354695A (zh) | 2017-01-25 |
| CN106354695B true CN106354695B (zh) | 2019-09-17 |
Family
ID=57844425
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201610701225.9A Expired - Fee Related CN106354695B (zh) | 2016-08-22 | 2016-08-22 | 一种仅输出线性时变结构模态参数辨识方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN106354695B (zh) |
Families Citing this family (13)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107561934B (zh) * | 2017-08-24 | 2020-04-21 | 北京理工大学 | 基于多参考时域传递率的仅输出结构模态参数辨识方法 |
| CN108416141A (zh) * | 2017-08-31 | 2018-08-17 | 北京理工大学 | 一种线性时变结构模态振型辨识方法 |
| CN109598175B (zh) * | 2017-09-30 | 2022-05-31 | 北京航空航天大学 | 一种基于多小波基函数和超正交前向回归的时频分析方法 |
| CN107943757B (zh) * | 2017-12-01 | 2020-10-20 | 大连理工大学 | 一种基于稀疏分量分析模态识别中的阶数确定方法 |
| CN107967395A (zh) * | 2017-12-11 | 2018-04-27 | 北京航空航天大学 | 一种基于beta小波基函数展开的时变非线性系统快速辨识方法 |
| CN109905187A (zh) * | 2017-12-11 | 2019-06-18 | 深圳先进技术研究院 | 一种非参数异常值检测方法、系统及电子设备 |
| CN108536971B (zh) * | 2018-04-13 | 2022-04-26 | 广州市建筑科学研究院有限公司 | 一种基于贝叶斯模型的结构损伤识别方法 |
| CN109521457B (zh) * | 2018-09-25 | 2022-10-21 | 中国辐射防护研究院 | 一种基于信息准则的γ辐射源项划分方法及系统 |
| CN109255333B (zh) * | 2018-09-25 | 2022-01-28 | 内蒙古工业大学 | 一种大型风电机组滚动轴承故障混合诊断方法 |
| CN111679651B (zh) * | 2020-06-08 | 2021-05-25 | 中国人民解放军火箭军工程大学 | 用于故障引起的变结构变参数系统的辨识方法及系统 |
| CN111985361A (zh) * | 2020-08-05 | 2020-11-24 | 武汉大学 | 小波降噪和emd-arima的电力系统负荷预测方法及系统 |
| CN112861074B (zh) * | 2021-03-09 | 2022-10-04 | 东北电力大学 | 基于Hankel-DMD的电力系统机电参数提取方法 |
| CN113808393A (zh) * | 2021-08-24 | 2021-12-17 | 东南大学 | 一种消除混杂控制对象影响的方法 |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102393645A (zh) * | 2011-11-07 | 2012-03-28 | 温州大学 | 一种高速电液比例调速系统的控制方法 |
| CN102982196A (zh) * | 2012-10-30 | 2013-03-20 | 北京理工大学 | 基于时变公分母模型的时频域时变结构模态参数辨识方法 |
| CN103001567A (zh) * | 2012-11-13 | 2013-03-27 | 江苏科技大学 | 六相永磁同步直线电机调速系统的多模逆模型辨识方法 |
| CN103728879A (zh) * | 2014-01-20 | 2014-04-16 | 华北电力大学 | 基于最小二乘支持向量机及在线更新的电站锅炉烟气软测量方法 |
-
2016
- 2016-08-22 CN CN201610701225.9A patent/CN106354695B/zh not_active Expired - Fee Related
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102393645A (zh) * | 2011-11-07 | 2012-03-28 | 温州大学 | 一种高速电液比例调速系统的控制方法 |
| CN102982196A (zh) * | 2012-10-30 | 2013-03-20 | 北京理工大学 | 基于时变公分母模型的时频域时变结构模态参数辨识方法 |
| CN103001567A (zh) * | 2012-11-13 | 2013-03-27 | 江苏科技大学 | 六相永磁同步直线电机调速系统的多模逆模型辨识方法 |
| CN103728879A (zh) * | 2014-01-20 | 2014-04-16 | 华北电力大学 | 基于最小二乘支持向量机及在线更新的电站锅炉烟气软测量方法 |
Non-Patent Citations (2)
| Title |
|---|
| Maximum likelihood estimator of operational modal analysis for linear time-varying structures in time-frequency domain;Si-Da Zhou et.;《Journal of Sound and Vibration》;20141231;第2339-2358页 |
| Parametric identification of a time-varying structure based on vector vibration response measurements;M.D. Spiridonakos et.;《Mechanical Systems and Signal Processing》;20091231;第2029-2048页 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN106354695A (zh) | 2017-01-25 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN106354695B (zh) | 一种仅输出线性时变结构模态参数辨识方法 | |
| CN107844627A (zh) | 一种仅输出时变结构模态参数贝叶斯估计方法 | |
| CN105843073B (zh) | 一种基于气动力不确定降阶的机翼结构气动弹性稳定性分析方法 | |
| Li et al. | Deep neural network for unsteady aerodynamic and aeroelastic modeling across multiple Mach numbers | |
| Zhang et al. | Intelligent fault diagnosis of rotating machinery using support vector machine with ant colony algorithm for synchronous feature selection and parameter optimization | |
| Batselier et al. | Tensor Network alternating linear scheme for MIMO Volterra system identification | |
| CN103048921B (zh) | 用于位置伺服系统的半周期重复控制器 | |
| KR20160021209A (ko) | 동적 모델 식별을 모니터링 및 가변 구조 또는 가변 동작 조건을 가진 동적 기기의 제어 방법 및 시스템 | |
| Zhou et al. | Output-only modal parameter estimator of linear time-varying structural systems based on vector TAR model and least squares support vector machine | |
| CN101916241B (zh) | 一种基于时频分布图的时变结构模态频率辨识方法 | |
| CN114065662B (zh) | 适用于网格拓扑可变的翼型流场快速预测方法 | |
| Cardot et al. | Confidence bands for Horvitz–Thompson estimators using sampled noisy functional data | |
| CN110826803A (zh) | 一种电力现货市场的电价预测方法及装置 | |
| CN104317195A (zh) | 一种基于改进极限学习机的非线性逆模型控制方法 | |
| CN115034337B (zh) | 一种轨道交通车辆中牵引电机状态辨识方法及装置、介质 | |
| CN109254321B (zh) | 一种地震激励下快速贝叶斯模态参数识别方法 | |
| CN110705041A (zh) | 一种基于easi的线性结构工作模态参数识别方法 | |
| Peterka et al. | Foundations of multivariate functional approximation for scientific data | |
| Gopinath | Efficient Fourier-based algorithms for time-periodic unsteady problems | |
| CN110210690B (zh) | 一种配电系统微型同步相量测量单元优化配置方法 | |
| CN107301266B (zh) | 一种磷酸铁锂电池loc估算方法和系统 | |
| Zhang et al. | A right-hand side function surrogate model-based method for the black-box dynamic optimization problem | |
| CN116885711A (zh) | 一种风电功率预测方法、装置、设备和可读存储介质 | |
| CN109903181B (zh) | 基于压缩感知的缺失数据集下的线损预测方法 | |
| CN103605880B (zh) | 一种精确诊断密集模态阻尼比的方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| CF01 | Termination of patent right due to non-payment of annual fee | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20190917 |