马希彬， 陈章位， 赵玉刚， 王 伟， 栾强利

(1.浙江大学 流体动力与机电系统国家重点实验室，杭州 310027； 2.杭州亿恒科技有限公司，杭州 310011)

振动试验通过复现给定的功率谱或振动波形再现真实的振动环境，常用于考核设备在相应振动情况下的性能，广泛应用于航空航天、机械制造、车辆与船舶、桥梁与建筑等众多工程测试和结构减振领域。振动试验系统从激振方向上分为单轴和多轴振动试验系统，从激励数目上分为单点激励和多点激励振动试验系统[1]。多轴多激励振动试验系统可提供更大的推力、实现应力的非均匀分布、提高信号的信噪比、缩短试验时间、提高故障的发现率等。

Sloane等[2]提出了随机功率谱密度估计均衡处理再现试验方法，推动了多轴振动试验控制的研究和发展；Smallwood[3]分析了单轴多激励振动实验系统随机振动控制驱动信号闭环修正算法，解决了多输入多输出耦合控制问题。本文基于已有成果，利用HV频响函数估计法对三轴振动实验系统进行系统辨识；针对辨识函数矩阵出现奇异值情况，设置阈值进行奇异值截断保证算法的稳定性，并用迭代控制算法修正驱动谱以提高振动控制精度；最后用集成该算法的多输入多输出振动控制器与三轴振动台进行试验。实验结果表明，基于驱动谱修正迭代控制算法对功率谱的复现具有较好的精度和工程实用性。

1 三轴振动试验控制系统的辨识[4-5]

三轴振动试验控制系统可以表述为多输入多输出系统模型，实际多输入多输出系统模型如图1所示，系统有n个驱动信号，m个响应信号，每一个驱动信号d(t)将对所有的控制信号c(t)产生作用，该系统的频域数学模型可描述为

c(f)-n(f)=H(f)(d(f)-m(f))

(1)

式中:c(f)为m×1维响应信号；d(f)为n×1维输入信号；m(f)为n×1维输入噪声信号；n(f)为m×1维输出噪声信号；H(f)为m×n维系统的频响函数。

图1 多输入多输出系统模型

频响函数估计实际上是一种有误差的测试数据的估计，其误差一部分来自功率谱估计算法的固有误差，一部分来自输入输出信号存在干扰而造成频率函数估计误差。针对不同的误差情况，应建立相应的模型，确定误差最小估计准则，分析频率响应估计结果。

1.1 输出噪声估计模型——H1估计

一般认为噪声信号与输入信号、输出信号不相关，当假设系统输入不存在噪声，即m(f)=0,式(1)可写成

c(f)=H(f)d(f)+n(f)

(2)

对式(1)两端同时右乘d(f)*(“*”为共轭转置),并求数学期望得

Gcd(f)=H(f)Gdd(f)+Gnd(f)

(3)

式中：Gcd(f)为系统输入输出之间的互功率谱；Gdd(f)为系统输入的自功率谱；Gnd(f)为系统的输出噪声与输入之间的互功率谱。

当n(f)与d(f)不相关时，即Gnd(f)=0,代入式(3)可得

Gcd(f)=H(f)Gdd(f)

(4)

若Gdd(f)可逆，则系统频响函数估计式可简化为

H1(f)=Gcd(f)Gdd(f)-1

(5)

多输入多输出系统理论上(无噪声)的频响函数H0(f)表达式为

H0(f)=Gvu(f)Guu(f)-1=Gvv(f)Guv(f)-1

(6)

式中：Gvu(f)为系统真实输出与真实输入之间的互功率谱；Guu(f)为系统真实输入的自功率谱；Gvv(f)为系统真实输出的自功率谱；Guv(f)为系统真实输入与真实输出之间的互功率谱。

显然有

Gdd(f)=E[d(f)d(f)*]=E[(u(f)+

m(f))(u(f)+m(f))*]=Guu(f)+Gum(f)+

Gmu(f)+Gmm(f)

(7)

Gcd(f)=E[c(f)d(f)*]=Gvu(f)+

Gvm(f)+Gnu(f)+Gnm(f)

(8)

由于噪声信号与输入信号、输出信号不相关，即Gum(f)=0,Gmu(f)=0,Gvm(f)=0,Gnu(f)=0,Gnm(f)=0。

式(5)可转换为

H1(f)=Gvu(f)(Guu(f)+Gmm(f))-1=

(9)

由式(9)可以看出H1估计法是欠估计，并且输入噪声对精度有很大影响。

1.2 输入噪声估计模型——H2估计

当假设系统输出不存在噪声，即n(f)=0,式(1)可写成

c(f)=H(f)(d(f)+m(f))

(10)

对式(10)两端同时右乘c(f)*(“*”为共轭转置)，并求数学期望整理得

Gcc(f)=H(f)(Gdc(f)+Gmc(f))

(11)

式中：Gcc(f)为系统输出的自功率谱；Gdc(f)为系统输入输出之间的互功率谱；Gmc(f)为系统的输入噪声与输出之间的互功率谱。

当m(f)与c(f)不相关时，即Gmc(f)=0代入式(11)可得

Gcc(f)=H(f)Gdc(f)

(12)

若Gdc(f)可逆，则系统频响函数估计式化简为

H2(f)=Gcc(f)Gdc(f)-1

(13)

Gcc(f)=E[c(f)c(f)*]=

E[(v(f)+n(f))(v(f)+n(f))*]=

Gvv(f)+Gvn(f)+Gnv(f)+Gnn(f)

(14)

Gdc(f)=E[d(f)c(f)*]=

Guv(f)+Gun(f)+Gmv(f)+Gmn(f)

(15)

由于噪声信号与输入信号、输出信号不相关，即Gnv(f)=0,Gvn(f)=0,Gun(f)=0,Gmv(f)=0,Gmn(f)=0。

式(13)可转换为

H2(f)=(Gvv(f)+Gnn(f))Guv(f))-1=

(16)

由式(16)可以看出H2估计法是过估计，并且输出噪声对精度有较大影响。

1.3 输入输出噪声估计模型—HV估计[6-8]

为了综合考虑系统输入输出噪声，提高频响估计精度，可根据最小二乘原理，极小化误差矩阵的方法建立频响函数估计HV模型。将系统频域模型式移项变换为

c(f)-H(f)d(f)=n(f)-H(f)m(f)

(17)

将式(17)两端同乘各自的共轭转置得

(c(f)-H(f)d(f))(c(f)-H(f)d(f))*=

(n(f)-H(f)m(f))(n(f)-H(f)m(f))*

(18)

对式(18)两端求数学期望，整理得

(19)

对于H1估计法和H2估计法来说，都需要进行矩阵求逆运算，实际测量信号都存在误差，且Gdd(f),Gdc(f)矩阵不一定可逆，这些因素导致矩阵病态，矩阵求逆时产生很大误差，从而增大频响函数估计误差。对于系统共振频率附近点，系统的输入信号信噪比较低，输出信号信噪比较高，因此在共振点附近H2估计比H1估计的精度高，反共振点附近，由于信噪比情况相反，所以H1估计比H2估计精度高。HV估计法综合考虑了系统输入输出噪声，在估计过程中对每个频率点上的输入输出互谱矩阵进行特征值分解，避免了矩阵求逆运算。故本文三轴振动实验系统选用HV估计法进行系统频响函数的辨识。

2 三轴振动试验系统解耦研究

本小节将对控制系统展开研究，使三轴振动的响应与设定的参考谱的误差保持在一定范围内，辨识控制过程如图2所示。

图2 系统辨识过程

Fig.2 System identification process

从理论上来讲，为使响应谱和参考谱一致，需要引入一个控制矩阵(解耦补偿矩阵Z(ω))对系统进行补偿。解耦补偿矩阵本质上就是系统频响函数矩阵的逆。显然，当系统的频响函数矩阵为方阵并且非奇异，即激励点与响应测量点个数相同(频响函数为n×n阶矩阵)，同时频响函数矩阵的行列式在控制频带内不等于零时，解耦补偿矩阵是唯一[9]，即Z(ω)=H(ω)-1当实际振动试验中激励点与响应测量点个数不相等时，系统频响函数矩阵称为长方阵，例如m×n阶矩阵。三轴振动试验系统工作在共振点处，结构的传递函数矩阵产生奇异，导致补偿矩阵与频响函数矩阵的乘积不等于单位阵，并由此产生较大的误差，导致在共振点处的谱值难以控制[10]。在系统频响函数为长方阵或为奇异阵时，解耦补偿矩阵是不存在的，为了求出系统的频响函数的逆，引入系统频响函数矩阵的广义逆或是伪逆，求解系统解耦补偿矩阵。选用奇异值分解法来求解奇异阵的伪逆。

设Hm×n为系统估计频响函数，秩为r，U∈Cm×m、V∈Cn×n为两正交矩阵，存在对角矩阵Λ∈Rm×n，有m×n矩阵H的奇异值分解如下(上标“*”表示矩阵的共轭转置)

(20)

式中：U和V分别为矩阵H的左奇异矩阵和右奇异矩阵；Λ为矩阵H的奇异值，且有Λr=diag(σ1,σ2,…,σr),σ1≥σ2≥…≥σr>0是方阵H*H特征值的平方根。因此,H的Moore-Penrose逆矩阵Z为

Z=H-1=VΛ-1U*

(21)

当矩阵H接近奇异时，Λ中至少有一个值接近零，Λ-1存在但值不稳定，用奇异值截断法Moore-Penrose逆矩阵求解补偿矩阵。

奇异值分解方法是基于Householder变换，QR算法迭代进行求解[11]。奇异值分解的数值算法精度很高，误差略大于机器误差，只要大于算法的误差的奇异值都会被保留下来，得到的奇异值矩阵中会有很小的奇异值，而且系统辨识越接近奇异，最小的奇异值越小。

矩阵的奇异性可用矩阵的奇异值来描述，例如方阵A是奇异阵，则A必定存在为零的奇异值。矩阵的奇异程度可用条件数来表示，而矩阵的条件数可用矩阵的特征值来描述，设方阵A可逆，则A的条件数可表示为

(22)

式中：cond(A)为相对于从属范数的矩阵A的条件数；σmax为方阵A的最大奇异值；σmin为方阵A的最小奇异值。

若H(f)奇异或接近奇异，对H(f)进行奇异值分解后奇异值大小相差较大，较大的奇异值及其相应的特征向量表示辨识模型中可靠部分，而不可靠部分通常是由较小的奇异值及其相应的特征向量表示，由此设置一个正的参考阈值，当矩阵内的奇异值小于阈值时，此时的奇异值剔除，即被视为零。该方法称为奇异值截断法，其具体实现过程如下：

(2)根据实际情况设置矩阵条件数M，当cond(Λr)>M时，则H(f)判定为奇异阵，对其进行阈值剔除处理。

(3)设置阈值，剔除小于该值的奇异值，得到新奇异矩阵Λnew，同时剔除U,V中不可靠的小奇异值所对应的列向量，得到新的正交矩阵Unew,Vnew。

3 三轴振动试验系统功率谱复现控制算法

由于在实际试验中存在干扰和一些不确定因素的影响，图2的控制过程中各个响应点上的响应信号都存在一定程度上的波动，导致控制结果仍然存在较大误差。由于其为开环控制策略，系统易受外界因素影响，具有很大的不稳定性，易损坏试验系统部件，因此并不能将其直接用于实际的工程中。为此需要引入一个反馈修正方案，通过不断地对驱动信号进行迭代修正，使控制点的幅值和相位满足要求。引入反馈后的控制修正过程如图3所示。

图3 多输入多输出振动试验修正控制框图

对于上述修正迭代控制系统，系统的初次辨识矩阵的逆并未得到修正，外部干扰信号将导致系统阻抗函数的解偏离，无法满足控制精度要求，为了消除干扰噪声对系统的影响，提出奇异值截断控制算法，系统解耦控制原理如图4所示，其控制过程主要包括频响函数估计、阻抗计算和驱动信号控制三部分。

图4 解耦控制算法框图

3.1 三轴振动试验系统频响函数估计

由图4知，由于频响函数与系统阻抗是根据系统的特性得到的，所以频响函数的辨识和阻抗计算的精度对整个系统的控制效果至关重要。考虑耦合特性的影响，输入、输出信号中噪声影响，根据前面讨论频响函数估计方法特点，采用HV估计法辨识出系统频响函数矩阵HV(f)。

3.2 三轴振动试验系统阻抗计算

采用Moore-Penrose逆矩阵对系统阻抗函数Z(f)的辨识，将辨识得到的系统阻抗函数Z(f)作为控制器传递函数，应用到系统辨识模型H(f)中，实现对系统有效控制。其计算过程如下

(2)通过设定矩阵条件数的大小M,对奇异值采用截断法进行阀值剔除处理；

(3)根据奇异值截断法处理后得到的新正交阵Unew,Vnew计算系统阻抗函数为

(23)

3.3 三轴振动试验系统驱动信号控制

功率谱复现试验采用驱动功率谱迭代修正控制策略，驱动功率谱迭代修正控制具有算法简单、收敛速度快的特点。由于系统传递函数的辨识精度直接决定了驱动功率谱迭代修正算法的精度和收敛性，当系统中存在较强的非线性因素时，系统辨识精度下降，可能导致驱动谱迭代修正算法发散。因此需要对系统辨识函数进行迭代修正以提高其精度，进而保证驱动谱修正算法收敛。驱动谱迭代修正控制算法包括对振动台系统阻抗函数(逆传递函数)的修正及振动台驱动功率谱的实时迭代和修正。振动台驱动功率谱迭代修正方程为

基于驱动谱修正迭代控制算法的功率谱复现试验控制过程主要分为5个步骤：

步骤1采用HV估计法辨识试验系统初始传递函数辨识并求逆，辨识得到的传递函数及其逆(阻抗函数)为HV(f)0,Z(f)0,关系可以表示为

(24)

步骤2初始驱动功率谱计算如下式，并生成驱动信号激励系统

(25)

d(t)0=IFFT(Gdd(f)0)

(26)

步骤3测量响应功率谱与控制误差谱

Gcc(f)k=FFT(c(t)k)

(27)

Gee(f)k=Grr(f)-Gcc(f)k

(28)

式中：Gcc(f)k为第k次迭代的响应谱；Gee(f)k为第k次迭代的控制误差谱。

步骤4根据第k次辨识的系统阻抗函数及误差，运用迭代算法修正系统的驱动功率谱。

(29)

步骤5通过频域随机化和时域随机化将驱动谱转换成时域驱动信号驱动系统。

4 三轴振动试验控制系统试验研究

为了验证论文系统建模及控制算法的有效性及精度，搭建三轴振动试验系统试验平台进行试验。图5是由杭州亿恒科技有限公司提供的三轴振动实验台；图6是由杭州亿恒科技有限公司开发的Premax多轴振动控制器，具有110 dB动态范围、精度高，多通道DSP(Digital Signal Processor)并行高速实时处理，易扩展及高可靠性等特点，PREMAX独立于PC机完成信号的分析和控制，采用千兆以太网与PC机相连，通过PC机显示对信号进行监测。

图5 三轴向振动试验台

三轴振动实验控制系统频响函数HV辨识结果如图7所示:系统X、Y轴向传递函数基本相同，共振点在1 200 Hz和1 800 Hz左右，反共振点在750 Hz左右，Z轴向传递函数与X、Y向的有一定的差异，共振点在1 300 Hz左右，反共振点在1 000 Hz左右，三轴向辨识得到的传递函数虽然存在共振峰及反共振峰(奇异点)，但其峰值较小，均控制在试验规定的范围内，由此可知用HV频响函数估计法精度较高，满足实验需求。

图6 Premax多轴振动控制器

图7 系统HV频率响应函数估计结果

三轴振动台功率谱复现试验结果如图8所示，目标加速度有效值为4g，图8从上往下依次可得X、Y、Z轴测得控制有效值分别为3.968 3g、3.953 1g、4.100 8g，功率谱密度范围设置为20～2 000 Hz，整个频率范围内为平直段功率谱密度，幅值为0.01 g2/Hz。由三轴振动试验系统自闭环控制试验结果可知：设置的频率范围内的功率谱都达到很好的控制效果，各轴控制响应谱均在目标谱的±3 dB范围内，满足工程试验要求。

(a)X轴控制结果

(b)Y轴控制结果

(c)Z轴控制结果

图8 三轴振动试验台X、Y、Z轴功率谱复现结果

Fig.8X、Y、Zaxis power spectrum retrieval result on three-axis vibration experimental table

5 结 论

本文利用HV频响函数估计法对三轴振动实验系统进行了系统辨识；针对辨识函数出现奇异值情况，设置阈值进行奇异值截断保证算法的稳定性，并用迭代控制算法修正驱动谱以提高振动控制精度；最后运用集成该控制算法的多输入多输出振动控制器与三轴振动试验台搭建的平台进行试验，对整个实验过程与结果进行分析有如下结论：

(1)三轴振动控制系统HV频响函数辨识法所得系统X、Y轴向传递函数基本相同，共振点在1 200 Hz和1 800 Hz左右，反共振点在750 Hz左右，Z轴向传递函数与X、Y向的有一定的差异，共振点在1 300 Hz左右，反共振点在1 000 Hz左右，三轴向辨识得到的传递函数虽然存在共振峰及反共振峰(奇异点)，但其峰值较小，均控制在试验规定的范围内，系统频响函数估计精度较高。

(2)三轴振动台功率谱复现试验显示，目标有效值为4g，X、Y、Z轴测得控制有效值分别为3.968 3g、3.953 1g、4.100 8g，均在规定范围值内。

(3)三轴振动试验系统自闭环控制试验在设置的频率范围内的功率谱都达到很好的控制效果，各轴上控制响应谱均在目标谱的±3 dB范围内，满足工程试验要求。

[1] CHONG K S, GWEE B H. A 16-channel low-power nonuniform spaced filter bank core for digital hearing aids[J]. IEEE Trans. on Circuits and Systems II: Express Briefs, 2006, 53(9):853-857.

[2] SLOANE E A. Vibration testing environment or apparatus[J]. Sound and Vibration Bulletin, 1973, 2: 165-182.

[3] SMALLWOOD D O. A random vibration control system for testing a single test item with multiple inputs[C]//SAE Aerospace Meeting. Warrendale: SAE, 1982.

[4] TUSTIN W. Fifth edition of Harris’ shock and vibration handbook reviewed[J]. Test Engineering and Management, 2002, 64(1): 30.

[5] 栾强利. 液压振动实验控制系统关键技术研究[D]. 杭州：浙江大学，2015.

[6] 吴吉平. 多输入多输出频域正交多项式模态参数识别方法[D]. 长沙：中南大学, 2002.

[7] 刘良玉. 多输入多输出频域模型参数识别[D]. 西安：西安电子科技大学, 2013.

[8] 杜永昌, 管迪华. 多输入多输出频域模态识别算法的探讨[J]. 清华大学学报(自然科学版), 1997, 37(11): 63-66.

DU Yongchang, GUAN Dihua. Study of multiple input multiple output frequency domain parameter estimation method[J]. Journal of Tsinghua University (Science and Technology), 1997, 37(11): 63-66.

[9] 朱银龙. 多输入多输出正弦振动试验控制系统的研究[D]. 南京：南京航空航天大学, 2007.

[10] 韩军, 鲍明, 倪宏伟. 进化规划在多振动台随机振动控制中的应用[J]. 中国机械工程, 2003, 14(13): 26-29.

HAN Jun, BAO Ming, NI Hongwei. An improved method based on evolutionary programming for multi-shaker random vibration control[J]. China Mechanical Engineering, 2003, 14(13): 26-29.

[11] 范锐. 轮耦合道路模拟台波形再现控制算法的研究[D]. 哈尔滨:哈尔滨工业大学, 2013.