发明名称 一种粘性阻尼振动信号中的模态参数提取方法
摘要 本发明公开了一种粘性阻尼振动信号中的模态参数提取方法,其通过选用较大的系统阶次的初始值,利用稀疏优化OMP算法计算状态矩阵,之后对可观矩阵中的每一行都运用OMP算法计算输出矩阵,并对所有的输出矩阵求取均值,较大程度上减少了噪声对结果的影响,从而提高了本发明方法的消噪能力和识别精度;然后计算固有频率、固有阻尼比和固有模态振型系数,并运用K‑means算法从众多模态参数中选出有效模态参数,来消除虚假模态,从而大大消弱了系统阶次对模态参数提取精度的影响。
申请公布号 CN106225914A 申请公布日期 2016.12.14
申请号 CN201610547068.0 申请日期 2016.07.11
申请人 宁波大学 发明人 李玉刚;叶庆卫;周宇;王晓东
分类号 G01H17/00(2006.01)I;G01M99/00(2011.01)I;G06F17/16(2006.01)I 主分类号 G01H17/00(2006.01)I
代理机构 宁波奥圣专利代理事务所(普通合伙) 33226 代理人 周珏
主权项 一种粘性阻尼振动信号中的模态参数提取方法,其特征在于包括以下步骤:①对连续的粘性阻尼振动信号进行奈奎斯特均匀采样,采样间隔为T<sub>S</sub>秒,采样点数为2N点,得到包含2N个采样值的采样信号,记为x,将x以列向量形式表示为x=(x<sub>1</sub>,x<sub>2</sub>,…,x<sub>2N‑1</sub>,x<sub>2N</sub>)<sup>T</sup>;其中,T<sub>S</sub>的取值满足奈奎斯特采样定律,N的取值要求远大于目标系统的系统阶次的预估值,(x<sub>1</sub>,x<sub>2</sub>,…,x<sub>2N‑1</sub>,x<sub>2N</sub>)<sup>T</sup>为(x<sub>1</sub>,x<sub>2</sub>,…,x<sub>2N‑1</sub>,x<sub>2N</sub>)的转置,x<sub>1</sub>,x<sub>2</sub>,…,x<sub>2N‑1</sub>,x<sub>2N</sub>对应表示x中的第1个采样值、第2个采样值、…、第2N‑1个采样值、第2N个采样值;②利用x中的所有采样值构建一个Hankel矩阵,记为H,H中的每一条副对角线上的元素都相等,<img file="FDA0001045539370000011.GIF" wi="701" he="599" />然后将H分成两个子矩阵,记为Y<sub>p</sub>和Y<sub>f</sub>,<img file="FDA0001045539370000012.GIF" wi="1388" he="306" />其中,H的维数为2i×j,Y<sub>p</sub>和Y<sub>f</sub>的维数均为i×j,i与j满足关系:2i+j‑1=2N,i的取值为人为确定的常数,且n&lt;i&lt;N,n表示目标系统的系统阶次的预估值,x<sub>3</sub>、x<sub>j</sub>、x<sub>j+1</sub>、x<sub>i</sub>、x<sub>i+1</sub>、x<sub>i+2</sub>、x<sub>i+3</sub>、x<sub>i+j‑1</sub>、x<sub>i+j</sub>、x<sub>i+j+1</sub>、x<sub>2i</sub>、x<sub>2i+1</sub>、x<sub>2i+j‑1</sub>对应表示x中的第3个采样值、第j个采样值、第j+1个采样值、第i个采样值、第i+1个采样值、第i+2个采样值、第i+3个采样值、第i+j‑1个采样值、第i+j个采样值、第i+j+1个采样值、第2i个采样值、第2i+1个采样值、第2i+j‑1个采样值;③利用Y<sub>p</sub>和Y<sub>f</sub>构建一个维数为i×i的Toeplitz矩阵,记为T,<img file="FDA0001045539370000021.GIF" wi="206" he="72" />然后选取T的第1行至第i‑1行构成维数为(i‑1)×i的第一子矩阵,记为T<sub>1</sub>;并选取T的第2行至第i行构成维数为(i‑1)×i的第二子矩阵,记为T<sub>2</sub>;其中,<img file="FDA0001045539370000022.GIF" wi="54" he="70" />为Y<sub>p</sub>的转置;④根据状态空间方程的定义对T<sub>1</sub>进行分解,得到T<sub>1</sub>的可观矩阵和可控反转矩阵,对应记为Γ<sub>1</sub>和Δ<sub>1</sub>,T<sub>1</sub>=Γ<sub>1</sub>Δ<sub>1</sub>,Γ<sub>1</sub>=(C<sub>,</sub>CA,…,CA<sup>i‑2</sup>)<sup>T</sup>;并对T<sub>1</sub>进行SVD分解,得到T<sub>1</sub>的第一正交矩阵、第二正交矩阵和正奇异矩阵组成的对角矩阵,对应记为U<sub>1</sub>、V<sub>1</sub>和S<sub>1</sub>,T<sub>1</sub>=U<sub>1</sub>S<sub>1</sub>V<sub>1</sub><sup>T</sup>;然后根据T<sub>1</sub>=Γ<sub>1</sub>Δ<sub>1</sub>和T<sub>1</sub>=U<sub>1</sub>S<sub>1</sub>V<sub>1</sub><sup>T</sup>,令<img file="FDA0001045539370000023.GIF" wi="244" he="63" />同样,根据状态空间方程的定义对T<sub>2</sub>进行分解,得到T<sub>2</sub>的可观矩阵和可控反转矩阵,对应记为Γ<sub>2</sub>和Δ<sub>2</sub>,T<sub>2</sub>=Γ<sub>2</sub>Δ<sub>2</sub>,Γ<sub>2</sub>=(CA,CA<sup>2</sup>,…,CA<sup>i‑1</sup>)<sup>T</sup>;并对T<sub>2</sub>进行SVD分解,得到T<sub>2</sub>的第一正交矩阵、第二正交矩阵和正奇异矩阵组成的对角矩阵,对应记为U<sub>2</sub>、V<sub>2</sub>和S<sub>2</sub>,<img file="FDA0001045539370000024.GIF" wi="271" he="71" />然后根据T<sub>2</sub>=Γ<sub>2</sub>Δ<sub>2</sub>和<img file="FDA0001045539370000025.GIF" wi="274" he="71" />令<img file="FDA0001045539370000026.GIF" wi="251" he="71" />其中,Γ<sub>1</sub>的维数为(i‑1)×1,Δ<sub>1</sub>的维数为1×i,C表示维数为1×(i‑1)的输出矩阵,A表示维数为(i‑1)×(i‑1)的状态矩阵,(C,CA,…,CA<sup>i‑2</sup>)<sup>T</sup>为(C,CA,…,CA<sup>i‑2</sup>)的转置,A<sup>i‑2</sup>为A的i‑2次方,V<sub>1</sub><sup>T</sup>为V<sub>1</sub>的转置,<img file="FDA0001045539370000027.GIF" wi="76" he="63" />为S<sub>1</sub>的<img file="FDA0001045539370000028.GIF" wi="47" he="110" />次方,Γ<sub>2</sub>的维数为(i‑1)×1,Δ<sub>2</sub>的维数为1×i,(CA,CA<sup>2</sup>,…,CA<sup>i‑1</sup>)<sup>T</sup>为(CA,CA<sup>2</sup>,…,CA<sup>i‑1</sup>)的转置,A<sup>2</sup>为A的2次方,A<sup>i‑1</sup>为A的i‑1次方,<img file="FDA0001045539370000029.GIF" wi="63" he="71" />为V<sub>2</sub>的转置,<img file="FDA00010455393700000210.GIF" wi="78" he="71" />为S<sub>2</sub>的<img file="FDA00010455393700000211.GIF" wi="42" he="111" />次方;⑤根据Γ<sub>1</sub>=(C<sub>,</sub>CA,…,CA<sup>i‑2</sup>)<sup>T</sup>和Γ<sub>2</sub>=(CA<sub>,</sub>CA<sup>2</sup>,…,CA<sup>i‑1</sup>)<sup>T</sup>,确定Γ<sub>1</sub>与Γ<sub>2</sub>之间的关系为:Γ<sub>2</sub><sup>T</sup>=Γ<sub>1</sub><sup>T</sup>Λ;然后将Γ<sub>1</sub><sup>T</sup>分解为<img file="FDA00010455393700000212.GIF" wi="218" he="63" />其中,Γ<sub>2</sub><sup>T</sup>为Γ<sub>2</sub>的转置,Γ<sub>1</sub><sup>T</sup>为Γ<sub>1</sub>的转置,Λ表示维数为(i‑1)×(i‑1)的对角矩阵,<img file="FDA00010455393700000213.GIF" wi="446" he="295" />P<sub>1</sub>表示维数为1×(i‑1)的行向量,P<sub>1</sub>=(C,C,…,C),Q<sub>1</sub>表示维数为(i‑1)×(i‑1)的对角矩阵,<img file="FDA0001045539370000031.GIF" wi="500" he="295" />I表示i‑1阶的单位矩阵;⑥根据稀疏优化的原理,确定稀疏度小于或等于i;并根据目标系统的系统阶次的预估值,确定稀疏度大于2n;然后确定稀疏度的取值范围为大于2n且小于或等于i;接着在稀疏度的取值范围内选择一个正整数k作为稀疏度的确定值,k∈(2n,i];再将大于2n的值i赋值给目标系统的系统阶次作为目标系统的系统阶次的初始值;⑦令q表示执行的次数,令Q表示重复执行的总次数;其中,q的初始值为1,Q≥2;⑧利用第q次执行时随机产生的第一高斯随机测量矩阵H<sub>1</sub>对Γ<sub>2</sub><sup>T</sup>=Γ<sub>1</sub><sup>T</sup>Λ进行观测,构建得到第q次执行时的第一稀疏优化模型,描述为:<img file="FDA0001045539370000032.GIF" wi="533" he="174" />并利用第q次执行时随机产生的第二高斯随机测量矩阵H<sub>2</sub>对<img file="FDA0001045539370000033.GIF" wi="194" he="63" />进行观测,构建得到第q次执行时的第二稀疏优化模型,描述为:<img file="FDA0001045539370000034.GIF" wi="581" he="182" />其中,H<sub>1</sub>的维数为M×1,H<sub>2</sub>的维数为M×(i‑1),M∈[k‑10,k+10],min为取最小值函数,符号“|| ||<sub>1</sub>”为求矩阵的1‑范数符号,s.t.表示“受约束于…”,符号“|| ||”为求欧氏距离符号,P<sub>1</sub><sup>T</sup>为P<sub>1</sub>的转置,Q<sub>1</sub><sup>T</sup>为Q<sub>1</sub>的转置,σ<sub>1</sub>和σ<sub>2</sub>均为值很小的常数;⑨根据目标系统的系统阶次的初始值和稀疏度的确定值,利用OMP方法对<img file="FDA0001045539370000035.GIF" wi="510" he="175" />进行稀疏求解,得到对角矩阵Λ;然后根据<img file="FDA0001045539370000036.GIF" wi="427" he="295" />和求解得到的Λ,计算得到A,并令A<sub>q</sub>=A;之后对A<sub>q</sub>进行特征值分解,得到A<sub>q</sub>的特征值向量,记为D<sub>q</sub>;接着根据<img file="FDA0001045539370000037.GIF" wi="475" he="299" />和得到的A,计算得到Q<sub>1</sub>;同样,根据目标系统的系统阶次的初始值和稀疏度的确定值,利用OMP方法对<img file="FDA0001045539370000041.GIF" wi="550" he="191" />进行稀疏求解,得到P<sub>1</sub><sup>T</sup>;最后根据P<sub>1</sub>=(C,C,…,C)得到C,并令C<sub>q</sub>=C;其中,A<sub>q</sub>和C<sub>q</sub>的初始值均为0;⑩根据A<sub>q</sub>和C<sub>q</sub>计算第q次执行得到的连续的粘性阻尼振动信号的固有频率向量、固有阻尼比向量和固有模态振型系数向量,对应记为<img file="FDA0001045539370000042.GIF" wi="147" he="71" />和<img file="FDA0001045539370000043.GIF" wi="82" he="71" />其中,<img file="FDA0001045539370000044.GIF" wi="150" he="71" />和<img file="FDA0001045539370000045.GIF" wi="64" he="71" />的维数均为1×(i‑1);<img file="FDA0001045539370000046.GIF" wi="54" he="55" />判断q≥Q是否成立,如果成立,则结束重复执行的过程,得到D<sub>1</sub>,D<sub>2</sub>,…,D<sub>q</sub>,…,D<sub>Q</sub>、C<sub>1</sub>,C<sub>2</sub>,…,C<sub>q</sub>,…,C<sub>Q</sub>、<img file="FDA0001045539370000047.GIF" wi="846" he="70" />和<img file="FDA0001045539370000048.GIF" wi="445" he="71" />然后执行步骤<img file="FDA0001045539370000049.GIF" wi="67" he="55" />否则,令q=q+1,然后返回步骤⑧继续执行;其中,q=q+1中的“=”为赋值符号;<img file="FDA00010455393700000410.GIF" wi="51" he="55" />对D<sub>1</sub>,D<sub>2</sub>,…,D<sub>q</sub>,…,D<sub>Q</sub>各自中的部分元素强制置零,再从D<sub>1</sub>,D<sub>2</sub>,…,D<sub>q</sub>,…,D<sub>Q</sub>中任意选择一个作为状态矩阵的最终特征值向量,记为D<sup>*</sup>;同样,对C<sub>1</sub>,C<sub>2</sub>,…,C<sub>q</sub>,…,C<sub>Q</sub>各自中的部分元素强制置零,再从C<sub>1</sub>,C<sub>2</sub>,…,C<sub>q</sub>,…,C<sub>Q</sub>中任意选择一个作为输出矩阵的最终值,记为C<sup>*</sup>;并计算Q次执行共得到的Q个固有频率向量<img file="FDA00010455393700000411.GIF" wi="395" he="70" />的均值向量<img file="FDA00010455393700000412.GIF" wi="83" he="55" />Q个固有阻尼比向量<img file="FDA00010455393700000413.GIF" wi="366" he="70" />的均值向量<img file="FDA00010455393700000414.GIF" wi="51" he="63" />和Q个固有模态振型系数向量<img file="FDA00010455393700000415.GIF" wi="422" he="69" />的均值向量<img file="FDA00010455393700000416.GIF" wi="83" he="62" /><img file="FDA00010455393700000417.GIF" wi="51" he="55" />根据D<sup>*</sup>中的每个零元素的索引,将<img file="FDA00010455393700000418.GIF" wi="54" he="54" />中对应索引的模态参数确定为虚假模态;而根据D<sup>*</sup>中的每个非零元素的索引,将<img file="FDA00010455393700000419.GIF" wi="53" he="55" />中对应索引的模态参数确定为有效模态参数,并提取出;根据D<sup>*</sup>中的每个零元素的索引,将<img file="FDA00010455393700000420.GIF" wi="55" he="68" />中对应索引的模态参数确定为虚假模态;而根据D<sup>*</sup>中的每个非零元素的索引,将<img file="FDA0001045539370000051.GIF" wi="51" he="63" />中对应索引的模态参数确定为有效模态参数,并提取出;根据C<sup>*</sup>中的每个零元素的索引,将<img file="FDA0001045539370000052.GIF" wi="57" he="53" />中对应索引的模态参数确定为虚假模态;而根据C<sup>*</sup>中的每个非零元素的索引,将<img file="FDA0001045539370000053.GIF" wi="59" he="47" />中对应索引的模态参数确定为有效模态参数,并提取出。
地址 315211 浙江省宁波市江北区风华路818号
您可能感兴趣的专利