CN105116935A - 一种基于频域变参数lms算法和陷波器的磁悬浮转子谐波振动抑制方法 - Google Patents

Abstract

本发明公开了一种基于频域变参数LMS(Least？Mean？Square,LMS)算法和陷波器的磁悬浮转子谐波振动抑制方法，针对磁悬浮控制力矩陀螺，首先建立含质量不平衡和传感器谐波的磁悬浮转子动力学模型，然后设计了一种基于频域变参数LMS算法和陷波器的谐波振动控制方法，利用频域变参数LMS算法抑制倍频振动，并利用陷波器抑制同频振动。本发明能对磁悬浮转子位移和磁轴承线圈电流中的谐波分量进行抑制，适用于存在质量不平衡和传感器谐波的磁悬浮转子系统谐波振动的主动控制。

Description

一种基于频域变参数LMS算法和陷波器的磁悬浮转子谐波振动抑制方法

技术领域

本发明涉及磁悬浮转子谐波抑制的技术领域，具体涉及一种频域变参数LMS算法和陷波器的磁悬浮转子谐波抑制方法，用于对磁悬浮控制力矩陀螺谐波振动进行抑制，从而使以磁悬浮控制力矩陀螺为惯性执行机构的卫星平台达到“超静”卫星平台的要求。

背景技术

无论是高分辨率对地观测的成像质量还是卫星间激光通信的可靠性都需要卫星提供可靠、微振动的“超静”卫星平台，这使得对超静卫星平台的研究具有非常重要的现实意义。磁悬浮控制力矩陀螺采用磁轴承支承，具有无摩擦、长寿命、主动振动可控等优点，是实现“超静”卫星平台的理想执行机构之一。

根据磁悬浮控制力矩陀螺转子的主动控制自由度的多少，可分为全主动磁悬浮控制力矩陀螺和主被动磁悬浮控制力矩陀螺。主被动磁悬浮控制力矩陀螺的部分自由度由主动磁轴承控制，较全主动磁悬浮控制力矩陀螺来说，具有低功耗、结构紧凑、体积小等优点。虽然磁悬浮控制力矩陀螺无摩擦，但是仍然存在一些高频振动，影响卫星平台的指向精度和稳定度。磁悬浮CMG的高频振动主要由转子质量不平衡和传感器谐波引起，包括同频振动和倍频振动，合称为谐波振动。其中，转子质量不平衡是主要振动源，由转子质量的不平衡引起。由于传感器检测面的圆度误差、材质不理想、电磁特性不均匀等原因，位移传感器信号中存在同频和倍频成分，即传感器谐波。随着转速的上升，当倍频成分接近转子的框架或壳体的固有模态时，会引起共振，振动幅值会急剧增加，严重影响超静卫星平台的高精度性能，因此需要对谐波振动加以抑制。

现有技术主要针对单一频率的振动进行抑制，对于谐波振动抑制研究相对较少，主要有陷波器、LMS算法、频域LMS算法等。其中，陷波器和LMS算法都是通过并联多个子系统对不同频率的谐波进行抑制，计算量随着频率成分的增多而增加，且需要考虑不同滤波器间的收敛速度问题，设计起来比较复杂。频域LMS算法，是在频域完成权值向量的自适应更新，无需并联多个滤波器变可对不同频率的扰动进行抑制，相对传统的LMS算法来说计算量减小。此外，为了抑制谐波振动，不仅需要抑制谐波电流刚度力，还要抑制谐波位移刚度力，目前主要是采用前馈补偿的方法对位移刚度力进行抑制。现有技术有以下问题：(1)对位移刚度力前馈补偿需要考虑功放对模型的影响。(2)采用频域LMS算法进行谐波振动抑制的收敛速度和稳态误差有待提高。

发明内容

本发明的目的为：克服现有技术的不足，发明一种频域变参数LMS算法和陷波器的磁悬浮转子谐波振动控制方法。采用构造轴承力的方法抑制谐波振动，该方法无需考虑功放的影响。此外，通过改变步长和块长提高振动抑制效果。

本发明采用的技术方案为：一种基于频域变参数LMS算法和陷波器的磁悬浮转子谐波振动抑制方法，包括以下步骤：

步骤(1)建立含质量不平衡和传感器谐波的磁悬浮转子动力学模型

两自由度主被动磁悬浮CMG转子径向磁轴承控制转子径向两自由度实现悬浮，其余三个自由度由安装在转子和定子上的永磁环实现无源稳定悬浮。设转子质心所在的中心平面为Ι，径向磁轴承的定子中心线与面Ι交于点N。转子几何轴和惯性轴分别交面Ι于O、C两点。在Ι内，以N为原点建立惯性坐标系NXY，以O为原点建立旋转坐标系Oεη。设OC长度为l，OC与Oε坐标轴的夹角为θ，O、C在坐标系NXY中的坐标分别为(x,y)、(X,Y)。

主被动磁轴承轴承力由主动磁轴承电磁力和被动磁轴承磁力两部分组成，X通道中轴承力f_x可写为：

f_x＝f_ex+f_px

其中，f_ex为X通道主动磁轴承电磁力，f_px为X通道被动磁轴承磁力。被动磁轴承产生的力与位移呈线性关系，表示为：

f_px＝K_prx

其中，K_pr是被动磁轴承位移刚度。

当转子悬浮在磁中心附近时，主动磁轴承电磁力可近似线性化为：

f_ex≈K_erx+K_ii_x

其中，K_er、K_i分别为主动磁轴承位移负刚度、电流刚度，i_x为功放输出电流。

对于含有质量不平衡的转子系统，有：

X(t)＝x(t)+Θ_x(t)

其中，X(t)为转子质心位移，x(t)为转子几何中心位移，Θ_x(t)为质量不平衡引起的位移扰动，记为：

Θ_x(t)＝lcos(Ωt+θ)

其中，l为质量不平衡的幅值，θ为相位，Ω为转子转速。

在实际转子系统中，由于机械加工精度和材料的不均匀因素的影响，传感器谐波不可避免，传感器实际测得的位移可表示为：

x_s(t)＝x(t)+x_d(t)

其中，x_d(t)为传感器谐波，可写为：

$x_d\left(t\right)=\underset{i=1}{\overset{n}{\Σ}}{c}_i\sin \left(i\mathrm{\Ω}t+{\mathrm{\θ}}_i\right)$

其中，c_i是传感器谐波系数的幅值，θ_i是传感器谐波系数的相位，n为传感器谐波的最高次数。

将i_x、X、Θ_x、x_d依次进行拉普拉斯变换得i_x(s)、X(s)、Θ_x(s)、x_d(s)，写出转子动力学方程有：

ms²X(s)＝(K_er+K_pr)(X(s)-Θ_x(s))+K_ii_x(s)

其中，

i_x(s)＝-K_sK_iG_c(s)G_w(s)(X(s)-Θ_x(s)+x_d(s))

其中，K_s为位移传感器环节、G_c(s)为控制器环节和G_w(s)为功放环节。则，有：

f_x(s)＝ms²X(s)＝(K_er+K_pr)(X(s)-Θ_x(s))-K_sK_iG_c(s)G_w(s)(X(s)-Θ_x(s)+x_d(s))

从上式可以看出，转子质量不平衡Θ_x(s)不仅会通过控制器产生电磁力，还会通过磁轴承本身产生永磁力，而传感器谐波仅通过磁轴承本身产生电磁力，即电磁力中不仅包含同频振动，还包含倍频振动，而永磁力中只包含同频振动，因此谐波振动抑制时需要加以区分。

步骤(2)传感器谐波中同频信号的辨识和补偿

让磁悬浮转子落在保护轴承上，控制磁悬浮转子以1Hz的转速低速旋转，由于此时转速很低，因此同频输出信号中由质量不平衡引起的同频信号很少，此时可认为传感器的同频输出不包含由质量不平衡引起的同频信号，均为传感器谐波，从而辨识出传感器谐波的同频信号。转子高速旋转时，加入辨识出的传感器谐波同频信号即完成传感器谐波同频信号的补偿。

步骤(3)基于频域变参数LMS算法的倍频振动抑制

对传感器谐波中的同频信号进行补偿后，传感器谐波只包含倍频信号。传感器谐波通过控制器产生倍频电磁力而不包含永磁力。因此，以倍频电流为控制对象即可对倍频振动进行抑制。将传感器谐波引起的倍频振动通过频域变参数LMS算法进行滤除，谐波振动抑制模块的计算过程如下：

为了抑制倍频振动，采用频域变参数LMS算法进行倍频振动抑制，以与倍频振动相关的正弦信号作为参考输入，即由与传感器倍频频率成分相同的单位幅值的正弦信号相加组成，误差作为基本输入。为了权衡收敛速度和稳态误差这一矛盾，实时地改变块长和步长，以达到更快的收敛速度和更小的稳态误差，根据误差变化情况设计块长和步长的更新算法，从而改善频域LMS算法的收敛性能。

步骤(4)基于陷波器的不平衡振动抑制

完成频域LMS算法进行倍频电流抑制和传感器谐波同频信号的补偿后，只剩下由质量不平衡引起的同频电磁力和永磁力，即不平衡振动力，可利用陷波器抑制不平衡振动力。直接以不平衡振动力为控制目标，将构造出的不平衡振动力作为陷波器的输入，输出反馈至控制器的输入。

本发明基本原理：对磁悬浮转子系统来讲，质量不平衡和传感器谐波引起的谐波振动会降低卫星平台的指向精度和稳定度，必须加以抑制。首先建立含质量不平衡和传感器谐波的磁悬浮转子动力学模型，通过分析轴承力的表达式分析谐波振动的主要原因。其次，采用基于频域变参数LMS算法和陷波器分别对传感器谐波引起的倍频振动和质量不平衡引起的同频振动进行抑制。此外，对传感器谐波引起的同频信号进行辨识和补偿。

本发明与现有技术相比的优点在于：为了有效抑制磁悬浮转子系统中主要频率成分的谐波振动，本发明提出一种基于频域变参数LMS算法和陷波器的谐波振动抑制方法，频域变参数LMS算法能提高收敛速度，减小稳态误差。此外，对于轴承力的抑制无需再对功放进行补偿，适用于存在质量不平衡和传感器谐波的磁悬浮转子谐波振动抑制。

附图说明

图1为本发明的流程图；

图2为X通道谐波振动抑制原理框图；

图3为Y通道谐波振动抑制原理框图；

图4为主被动磁悬浮转子系统结构示意图；

图5为传感器谐波示意图；

图6为X通道磁轴承控制系统框图；

图7为Y通道磁轴承控制系统框图；

图8为X通道基于频域变参数LMS算法的倍频振动抑制模块；

图9为X通道基于陷波器的不平衡振动抑制模块。

具体实施方式

下面结合附图以及具体实施例进一步说明本发明。

如图1所示，一种基于频域变参数LMS算法和陷波器的磁悬浮转子谐波振动抑制方法的实施过程是：首先建立含质量不平衡和传感器谐波的磁悬浮转子动力学模型，然后设计一种基于快速块变换的谐波振动抑制方法。即通过建立含质量不平衡和传感器谐波的磁悬浮转子动力学模型，分析谐波振动的主要原因。为了抑制谐波振动，首先对传感器谐波引起的同频信号进行辨识和补偿，此时传感器谐波只会引起倍频振动，而同频振动完全由质量不平衡导致。采用基于频域变参数LMS算法对传感器谐波引起的倍频振动进行抑制，最后，用陷波器对质量不平衡引起的不平衡振动进行抑制。图2和图3为基于以上流程的谐波振动原理框图，图2为本发明所采用的磁悬浮转子系统X通道谐波振动抑制方法原理框图。图3为本发明所采用的磁悬浮转子系统Y通道谐波振动抑制方法原理框图。谐波振动抑制主要包括传感器同频信号补偿1、倍频振动抑制模块2和不平衡振动抑制模块3。对传感器输出的同频信号进行辨识，将辨识出的传感器同频信号补偿，在此基础上，采用频域LMS算法对倍频振动进行抑制，最后，构造出轴承力，直接以轴承力为控制对象用陷波器抑制不平衡振动。

(1)建立含质量不平衡和传感器谐波的磁悬浮转子动力学模型

两自由度主被动磁悬浮CMG转子径向磁轴承控制转子径向两自由度实现悬浮，其余三个自由度由安装在转子和定子上的永磁环实现无源稳定悬浮，如图4所示。磁悬浮控制力矩陀螺依据支承转子的磁轴承主动控制自由度个数，可以分为全主动磁悬浮控制力矩陀螺和主被动磁悬浮控制力矩陀螺。设转子质心所在的中心平面为Ι，径向磁轴承的定子中心线与面Ι交于点N。转子几何轴和惯性轴分别交面Ι于O、C两点。在Ι内，以N为原点建立惯性坐标系NXY，以O为原点建立旋转坐标系Oεη。设OC长度为l，OC与Oε坐标轴的夹角为θ，O、C在坐标系NXY中的坐标分别为(x,y)、(X,Y)。

主被动磁轴承轴承力由主动磁轴承电磁力和被动磁轴承磁力两部分组成。以X通道为例，轴承力f_x可写为：

f_x＝f_ex+f_px

其中，f_ex为X通道主动磁轴承电磁力，f_px为X通道被动磁轴承磁力。被动磁轴承产生的力与位移呈线性关系，表示为：

f_px＝K_prx

其中，K_pr是被动磁轴承位移刚度。

当转子悬浮在磁中心附近时，主动磁轴承电磁力可近似线性化为：

f_ex≈K_erx+K_ii_x

其中，K_er、K_i分别为主动磁轴承位移负刚度、电流刚度，i_x为功放输出电流。

对于含有质量不平衡的转子系统，有：

X(t)＝x(t)+Θ_x(t)

其中，X(t)为转子质心位移，x(t)为转子几何中心位移，Θ_x(t)为质量不平衡引起的位移扰动，记为：

Θ_x(t)＝lcos(Ωt+θ)

其中，l为质量不平衡的幅值，θ为相位，Ω为转子转速。

在实际转子系统中，由于机械加工精度和材料的不均匀等因素的影响，传感器谐波不可避免，如图5所示。传感器实际测得的位移可表示为：

x_s(t)＝x(t)+x_d(t)

其中，x_d(t)为传感器谐波，可写为：

$x_d\left(t\right)=\underset{i=1}{\overset{n}{\Σ}}{c}_i\sin \left(i\mathrm{\Ω}t+{\mathrm{\θ}}_i\right)$

其中，c_i是传感器谐波系数的幅值，θ_i是传感器谐波系数的相位，n为传感器谐波的最高次数。

将i_x、X(t)、Θ_x(t)、x_d(t)依次进行拉普拉斯变换得i_x(s)、X(s)、Θ_x(s)、x_d(s)，图6为X通道磁轴承控制系统框图，磁轴承控制系统由控制器G_c(s)、功率放大器G_w(s)、转子系统组成P(s)，写出转子动力学方程有：

ms²X(s)＝(K_er+K_pr)(X(s)-Θ_x(s))+K_ii_x(s)

其中，

i_x(s)＝-K_sK_iG_c(s)G_w(s)(X(s)-Θ_x(s)+x_d(s))

其中，K_s为位移传感器环节、G_c(s)为控制器环节和G_w(s)为功放环节。可写出f_x(s)的表达式，有：

f_x(s)＝ms²X(s)＝(K_er+K_pr)(X(s)-Θ_x(s))-K_sK_iG_c(s)G_w(s)(X(s)-Θ_x(s)+x_d(s))

相应的，Y通道磁轴承控制系统框图如图7所示，Y通道的轴承力f_y(s)可写为：

f_y(s)＝ms²Y(s)＝(K_er+K_pr)(Y(s)-Θ_y(s))-K_sK_iG_c(s)G_w(s)(Y(s)-Θ_y(s)+y_d(s))

其中，Θ_y(s)和y_d(s)分别为Θ_y(t)和y_d(t)的拉普拉斯变换，Θ_y(t)＝lsin(Ωt+θ)，可以看出，转子质量不平衡不仅会通过控制器产生电磁力，还会通过磁轴承本身产生永磁力，而传感器谐波仅通过磁轴承本身产生电磁力，即电磁力中不仅不含同频振动，还包含倍频振动，而永磁力中只包含同频振动，因此谐波振动抑制时需要加以区分。

在建立含质量不平衡和传感器谐波的磁悬浮转子动力学模型的基础上，通过分析轴承力的表达式，得到谐波振动的主要原因。为了抑制谐波振动，需要将电磁力和永磁力区分开，即将传感器谐波和质量不平衡产生的同频振动区分开，因此需要对传感器谐波引起的同频信号进行辨识，并进行补偿。

(2)传感器谐波中同频信号的辨识和补偿

质量不平衡引起的扰动与转速有关，而传感器谐波引起的同频干扰与转速无关，因此可以利用转速对传感器谐波同频信号辨识并进行补偿。

让磁悬浮转子落在保护轴承上，控制磁悬浮转子以1Hz的转速低速旋转，由于此时转速很低，因此同频输出信号中由质量不平衡引起的同频信号很少，此时可认为传感器的同频输出不包含由质量不平衡引起的同频信号，均为传感器谐波，从而辨识出传感器谐波的同频信号。转子高速旋转时，加入辨识出的传感器谐波同频信号即完成传感器谐波同频信号的补偿。补偿后的传感器输出x_s'(t)、y_s'(t)为：

$\left\{\begin{array}{c}{x_s}^{\′}\left(t\right)=x\left(t\right)+x_d\left(t\right)-{\hat{x}}_d\left(t\right)\\ {y_s}^{\′}\left(t\right)=y\left(t\right)+y_d\left(t\right)-{\hat{y}}_d\left(t\right)\end{array}\right.$

其中，y(t)为转子Y通道几何中心位移，为X通道的传感器的同频补偿信号，为Y通道传感器的同频补偿信号， ${\hat{x}}_d\left(t\right)={\hat{c}}_{1}sin(\mathrm{\Ω}t+\widehat{\mathrm{\θ}}),{\hat{y}}_d\left(t\right)={\hat{c}}_{1}cos(\mathrm{\Ω}t+\widehat{\mathrm{\θ}}),$ 为传感器谐波同频信号幅值和相位的补偿值。可进一步写出补偿后的传感器输出为：

$\left\{\begin{array}{c}{x_s}^{\′}\left(t\right)=\underset{i=2}{\overset{n}{\Σ}}{c}_i\sin \left(i\mathrm{\Ω}t+{\mathrm{\θ}}_i\right)+\[c_1\sin \left(\mathrm{\Ω}t+\mathrm{\θ}\right)-{\hat{c}}_1\sin \left(\mathrm{\Ω}t+\widehat{\mathrm{\θ}}\right)\]=\underset{i=2}{\overset{n}{\Σ}}{c}_i\sin \left(i\mathrm{\Ω}t+{\mathrm{\θ}}_i\right)\\ {y_s}^{\′}\left(t\right)=\underset{i=2}{\overset{n}{\Σ}}{c}_i\cos \left(i\mathrm{\Ω}t+{\mathrm{\θ}}_i\right)+\[c_1\cos \left(\mathrm{\Ω}t+\mathrm{\θ}\right)-{\hat{c}}_1\cos \left(\mathrm{\Ω}t+\widehat{\mathrm{\θ}}\right)\]=\underset{i=2}{\overset{n}{\Σ}}{c}_i\cos \left(i\mathrm{\Ω}t+{\mathrm{\θ}}_i\right)\end{array}\right.$

从上式看出，对传感器进行同频信号补偿后，传感器的输出只有倍频信号。

(3)基于频域变参数LMS算法的倍频振动抑制

对传感器谐波中的同频信号进行补偿后，传感器谐波只包含倍频信号。传感器谐波通过控制器产生倍频电磁力而不包含永磁力。因此，以倍频电流为控制对象即可对倍频振动进行抑制。本发明采用一种基于频域变参数LMS算法对倍频振动进行抑制。以X通道为例，倍频振动中只包含由传感器谐波引起的电磁力，因此本发明以电流为控制目标抑制倍频振动，采用频域LMS算法，以电流为基本输入，以与倍频电流相关的正弦信号作为参考输入，倍频振动抑制模块的输出反馈至原控制系统的功放输入端，其倍频振动抑制模块如图8所示，以电流为基本输入，以与倍频电流相关的正弦信号作为参考输入，倍频振动抑制模块的输出反馈至原控制系统的功放输入端，计算过程如下：

为了便于描述，在本文中，时域变量用小写字母表示，频域变量用大写字母表示，向量和矩阵用黑体表示，标量用斜体表示，F(·)表示傅里叶变换，F^-1(·)表示傅里叶逆变换。

设滤波器长度和块长均为N_x。即滤波器长度和块长大小一样。输入信号通过串并转换变为大小为N_x的块。定义第k块参考输入u_x(k)、基本输入e_x(k)、滤波器输出out_x(k)、N_x阶滤波器抽头权向量w_x(k)依次表示为：

u_x(k)＝[u_x(kN_x-N_x),u_x(kN_x-N_x+1),...,u_x(kN_x+N_x-1)]^T

e_x(k)＝[e_x(kN_x),e_x(kN_x+1),...,e_x(kN_x+N_x-1)]^T

out_x(k)＝[out_x(kN_x-N_x),out_x(kN_x-N_x+1),...,out_x(kN_x+N_x-1)]^T

w_x(k)＝[w_0,x(k),w_1,x(k),...,w_N-1,x(k)]^T

滤波器输出out_x(k)为：

out_x(k)＝k_xF^-1OUT_x(k)＝k_xF^-1[U_x(k)W_x(k)]

其中，OUT_x(k)为第k块滤波器的频域输出，U_x(k)＝diag{F[u_x(k)]}，W_x(k)为根据权值更新公式计算出的频域权值向量，k_x是N_x×2N_x阶约束矩阵：

$k_x=\[O_{N_x},I_{N_x}\]$

其中，是N_x×N_x阶零阵，是N_x×N_x阶单位阵。

基于随机梯度下降原则，滤波器抽头向量权值更新写为：

W_x(k+1)＝W_x(k)+μ_x(k)F{g_xF^-1[U_x ^H(k)E_x(k)]}

其中，W_x(k+1)为计算出的第k+1时刻的频域权值向量，U_x ^H(k)为参考输入频域向量的共轭，E_x(k)为第k时刻的频域误差向量，μ_x(k)为步长，g_x为梯度约束矩阵：

$g_x=\left[\begin{array}{cc}{I}_{N_x}& O_{N_x}\\ O_{N_x}& O_{N_x}\end{array}\right]$

步长μ_x(k)控制滤波器抽头权向量从算法的当前迭代到下一次迭代的增量变化，为了保持系统稳定性，需满足：

$0<{\mathrm{\&mu;}}_x\left(k\right)<\frac{1}{N_x{\mathrm{\&lambda;}}_{max}}$

其中，λ_max是参考输入信号相关矩阵R_x＝E[u_x(n)u_x(n)^T]的最大特征值。步长μ_x(k)的不同取值会影响收敛速度和精度，常采用平均时间常数τ_mse,av和失调系数M这两个量作为衡量频域LMS算法性能的指标。平均时间常数τ_mse,av写为：

${\mathrm{\&tau;}}_{mse,av}=\frac{N_x}{4{\mathrm{\&mu;}}_x\left(k\right){\mathrm{\&lambda;}}_{av}}$

式中，λ_av是矩阵R_x的特征值的平均值，平均时间常数能反映系统的收敛速度。失调系数M写为：

$M=\frac{{\mathrm{\&mu;}}_x\left(k\right)}{N_x}tr\&lsqb;R_x\&rsqb;$

其中，tr[R_x]是矩阵R_x的迹。失调系数可以反映稳态误差。

结合平均时间常数τ_mse,av和失调系数M的表达式可以看出，块长N_x和步长μ_x(k)均会影响收敛速度和稳态误差。因此，为了权衡收敛速度和稳态误差这一矛盾，可以实时地改变块长和步长，以达到更快的收敛速度和更小的稳态误差。在初始时刻，可以采用小块长、大步长提高收敛速度，随着误差不断减小，增大块长减小步长来减小稳态误差。具体方法如下：

设当前块长为N_x，下一块块长为N_x'定义估计误差e_x,ave(k)记为：

$e_{x,ave}\left(k\right)={\mathrm{\&alpha;e}}_{x,ave}(k-1)+\frac{(1-\mathrm{\&alpha;})}{N_x}\underset{i=0}{\overset{N_x-1}{\&Sigma;}}{e}_x({\mathrm{kN}}_x+i)$

其中，e_x,ave(k-1)是第k-1块的估计误差，α为一个常数且0<α<1，e_x(j)表示第j时刻X通道的基本输入，j＝kN_x,kN_x+1,...,kN_x+N_x-1。块长更新的主要思想是：如果当前块的估计误差小于上一块的估计误差，块长增大；如果当前块的估计误差大于上一块的估计误差，块长减小。然而，由于前后两块的估计误差肯定不完全相等，导致每块的块长都会改变，这样会增加计算量。因此，前后两块进行对比时留有一定变化范围，块长的具体更新算法如下：

若β₁e_x,ave(k-1)<e_x,ave(k)<β₂e_x,ave(k-1)，则视为估计误差基本没变化，N_x'＝N_x。其中，β₁、β₂为常数，且0<β₁<1，β₂>1，用来保证前后两块比较时留有一定的余量，即在当前块与上一块的估计误差相差在一定范围内时认为误差基本不变，不用改变块长。

若e_x,ave(k)≥β₂e_x,ave(k-1)，则视为当前块的估计误差大于上一块的估计误差，此时，N_x'＝N_x/2。

若e_x,ave(k)≤β₁e_x,ave(k-1)，则视为当前块的估计误差小于上一块的估计误差，此时，N_x'＝2N_x。

块长变化时，权值向量也需要相应地改变。权值向量的更新算法如下：

若N_x'＝2N_x，需要增加N_x个权值向量。由于在时域中补零等效于在频域中插值，因此可以通过对原系数时域补零再频域变换得到新的权系数，即

w_x'(k+1)＝[w_x(k+1),0,...,0]^T

W_x'(k+1)＝F[w_x'(k+1)]

其中，w_x'(k+1)为根据块长更新情况得到的第k+1时刻的时域权值向量，w_x(k+1)为计算出的第k+1时刻的时域权值向量，W_x'(k+1)为根据块长更新情况得到的第k+1时刻的频域权值向量。

若N_x'＝N_x/2，需要减少N_x/2个权值向量，有

${W_x}^{\&prime;}(k+1)=\&lsqb;W_{x,0}(k+1),W_{x,2}(k+1),\mathrm{...},W_{x,2N_x}(k+1)\&rsqb;$

即每隔一个舍弃一个权值向量，W_x,i(k+1),i＝0,1,...,2N_x为计算出的第k+1时刻的频域权值向量。

若N_x'＝N_x，权值向量无需改变。

对于每一块采用不同步长以更好地改善收敛性能。特别地，频域LMS的收敛速度可以在不影响最小均方误差的基础上，通过对每个可调权值赋予不同的步长获得改善。步长更新算法为：

${\mathrm{\&mu;}}_{x,i}\left(k\right)=\frac{{\mathrm{\&mu;}}_0\frac{e_{x,ave}\left(k\right)}{e_{x,ave}\left(k-1\right)}}{\left|\right|U_{x,i}\left(k\right)|{|}^2+\mathrm{\&gamma;}}$

其中，μ_x,i(k)为第k块中第i个步长，U_x,i(k)为第k块中第i个参考输入的频域信号，||·||表示模值，μ₀起到控制步长大小的作用，γ是一个常数，保证步长不会因输入信号能量过低而过大。

同理可得Y通道的更新算法。滤波器抽头向量权值更新表示为：

W_y(k+1)＝W_y(k)+μ_y(k)F{g_yF^-1[U_y ^H(k)E_y(k)]}

其中，W_y(k)为Y通道第k块频域权值向量，W_y(k+1)为计算出的Y通道第k+1块频域权值向量，U_y ^H(k)为参考输入频域向量的共轭，U_y(k)＝diag{F[u_y(k)]}，u_y(k)为参考输入的时域向量，E_y(k)为第k块频域误差向量，μ_y(k)为步长，g_y为梯度约束矩阵。滤波器输出out_y(k)为：

out_y(k)＝k_yF^-1OUT_y(k)＝k_yF^-1[U_y(k)W_y(k)]

其中，OUT_y(k)为第k块滤波器的频域输出。k_y是N_y×2N_y阶约束矩阵：

$k_y=\&lsqb;O_{N_y},I_{N_y}\&rsqb;$

其中，是N_y×N_y阶零阵，是N_y×N_y阶单位阵。设当前块长为N_y，下一块块长为N_y'，定义估计误差e_y,ave(k)记为：

$e_{y,ave}\left(k\right)={\mathrm{\&alpha;e}}_{y,ave}(k-1)+\frac{(1-\mathrm{\&alpha;})}{N_y}\underset{i=0}{\overset{N_y-1}{\&Sigma;}}{e}_y({\mathrm{kN}}_y+i)$

其中，e_y,ave(k-1)是第k-1块的估计误差，e_y(j)表示第j时刻Y通道的基本输入，j＝kN_y,kN_y+1,...,kN_y+N_y-1。块长的具体更新算法如下：

若β₁e_y,ave(k-1)<e_y,ave(k)<β₂e_y,ave(k-1)，则视为估计误差基本没变化，N_y'＝N_y。

若e_y,ave(k)≥β₂e_y,ave(k-1)，则视为当前块的估计误差大于上一块的估计误差，此时，N_y'＝N_y/2。

若e_y,ave(k)≤β₁e_y,ave(k-1)，则视为当前块的估计误差小于上一块的估计误差，此时，N_y'＝2N_y。

权值向量的更新算法如下：

若N_y'＝2N_y，需要增加N_y个权值向量。由于在时域中补零等效于在频域中插值，因此可以通过对原系数时域补零再频域变换得到新的权系数，即：

w_y'(k+1)＝[w_y(k+1),0,...,0]^T

W_y'(k+1)＝F[w_y'(k+1)]

其中，w_y'(k+1)为根据块长更新情况得到的第k+1时刻的时域权值向量，w_y(k+1)为计算出的第k+1时刻的时域权值向量，W_y'(k+1)为根据块长更新情况得到的第k+1时刻的频域权值向量。

若N_y'＝N_y/2，需要减少N_y/2个权值向量，有：

${W_y}^{\&prime;}(k+1)=\&lsqb;W_{y,0}(k+1),W_{y,2}(k+1),\mathrm{...},W_{y,2N_y}(k+1)\&rsqb;$

即每隔一个舍弃一个权值向量，W_y,i(k+1),i＝0,1,...,2N_y为计算出的第k+1时刻的频域权值向量。

若N_y'＝N_y，权值向量无需改变。

步长更新算法为：

${\mathrm{\&mu;}}_y\left(k\right)=\frac{{\mathrm{\&mu;}}_0\frac{e_{y,ave}\left(k\right)}{e_{y,ave}(k-1)}}{\left|\right|U_{y,i}\left(k\right)|{|}^2+\mathrm{\&gamma;}}$

其中，U_y,i(k)为第k块中第i个参考输入的频域信号。

(4)基于陷波器的不平衡振动抑制

完成频域LMS算法进行倍频振动抑制和传感器谐波同频信号的补偿后，只剩下由质量不平衡引起的同频电磁力和永磁力，即不平衡振动力，可利用陷波器抑制不平衡振动力。本发明直接以不平衡振动力为控制目标，以X通道为例，将构造出的不平衡振动力F_x(s)作为陷波器的输入，输出反馈至控制器的输入。构造出的不平衡振动力F_x(s)可写为：

F_x(s)＝K_ii_x(s)+(K_er+K_pr)x_s'(s)

如图9所示，不平衡振动抑制模块主要由跟踪滤波器N_f(s)构成，将构造出的不平衡振动力作为跟踪滤波器的输入，跟踪滤波器的输出只含与转速有关的信号，跟踪滤波器N_f(s)的传递函数为：

$N_f\left(s\right)=\frac{s}{s^2+{\mathrm{\&Omega;}}^2}$

以质量不平衡Θ_x(s)为输入，轴承力F_x(s)为输出，对应的传递函数为：

$\frac{F_x\left(s\right)}{{\mathrm{\&Theta;}}_x\left(s\right)}=\frac{K_{er}+K_{pr}-K_i{K}_s{G}_c\left(s\right)G_w\left(s\right)}{1+K_i{K}_s{G}_c\left(s\right)G_w\left(s\right)P\left(s\right)-(K_{er}+K_{pr})P\left(s\right)+{\mathrm{\&epsiv;K}}_i{N}_f\left(s\right)G_w\left(s\right)}$

其中，P(s)为转子系统传递函数。则，有：

$\underset{s\&RightArrow;j\mathrm{\&Omega;}}{\lim }{F}_x\left(j\mathrm{\&Omega;}\right)=0$

由上式可得，该陷波器可以对不平衡振动进行抑制。

同理可得Y通道不平衡振动抑制方法。将构造出的不平衡振动力F_y(s)作为陷波器的输入，输出反馈至控制器的输入。构造出的不平衡振动力F_y(s)可写为：

F_y(s)＝K_ii_y(s)+(K_er+K_pr)y_s'(s)

以质量不平衡Θ_y(s)为输入，轴承力F_y(s)为输出，对应的传递函数为：

$\frac{F_y\left(s\right)}{{\mathrm{\&Theta;}}_y\left(s\right)}=\frac{K_{er}+K_{pr}-K_i{K}_s{G}_c\left(s\right)G_w\left(s\right)}{1+K_i{K}_s{G}_c\left(s\right)G_w\left(s\right)P\left(s\right)-(K_{er}+K_{pr})P\left(s\right)+{\mathrm{\&epsiv;K}}_i{N}_f\left(s\right)G_w\left(s\right)}$

则，有：

$\underset{s\&RightArrow;j\mathrm{\&Omega;}}{\lim }{F}_y\left(j\mathrm{\&Omega;}\right)=0$

由上式可得，该陷波器可以对不平衡振动进行抑制。

本发明说明书中未作详细描述的内容属于本领域专业技术人员公知的现有技术。

Claims (2)

1.一种基于频域变参数LMS算法和陷波器的磁悬浮转子谐波振动抑制方法，其特征在于：包括以下步骤：

步骤(1)建立含质量不平衡和传感器谐波的磁悬浮转子动力学模型

两自由度主被动磁悬浮CMG(ControlMomentGyro,CMG)转子径向磁轴承控制转子径向两自由度实现悬浮，其余三个自由度由安装在转子和定子上的永磁环实现无源稳定悬浮，设转子质心所在的中心平面为Ι，径向磁轴承的定子中心线与面Ι交于点N，转子几何轴和惯性轴分别交面Ι于O、C两点，在Ι内，以N为原点建立惯性坐标系NXY，以O为原点建立旋转坐标系Oεη，设OC长度为l，OC与Oε坐标轴的夹角为θ，O、C在坐标系NXY中的坐标分别为(x,y)、(X,Y)；

主被动磁轴承轴承力由主动磁轴承电磁力和被动磁轴承磁力两部分组成，X通道中轴承力f_x可写为：

f_x＝f_ex+f_px

其中，f_ex为X通道主动磁轴承电磁力，f_px为X通道被动磁轴承磁力，被动磁轴承产生的力与位移呈线性关系，表示为：

f_px＝K_prx

其中，K_pr是被动磁轴承位移刚度；

当转子悬浮在磁中心附近时，主动磁轴承电磁力可近似线性化为：

f_ex≈K_erx+K_ii_x

其中，K_er、K_i分别为主动磁轴承位移负刚度、电流刚度，i_x为功放输出电流；

对于含有质量不平衡的转子系统，有：

X(t)＝x(t)+Θ_x(t)

其中，X(t)为转子质心位移，x(t)为转子几何中心位移，Θ_x(t)为质量不平衡引起的位移扰动，记为：

Θ_x(t)＝lcos(Ωt+θ)

其中，l为质量不平衡的幅值，θ为相位，Ω为转子转速；

在实际转子系统中，由于机械加工精度和材料的不均匀的影响，传感器谐波不可避免，传感器实际测得的位移可表示为：

x_s(t)＝x(t)+x_d(t)

其中，x_d(t)为传感器谐波，可写为：

$x_d\left(t\right)=\underset{i=1}{\overset{n}{\&Sigma;}}{c}_i\sin \left(i\mathrm{\&Omega;}t+{\mathrm{\&theta;}}_i\right)$

其中，c_i是传感器谐波系数的幅值，θ_i是传感器谐波系数的相位，n为传感器谐波的最高次数；

将i_x、X、Θ_x、x_d依次进行拉普拉斯变换得i_x(s)、X(s)、Θ_x(s)、x_d(s)，写出转子动力学方程有：

ms²X(s)＝(K_er+K_pr)(X(s)-Θ_x(s))+K_ii_x(s)

其中，

i_x(s)＝-K_sK_iG_c(s)G_w(s)(X(s)-Θ_x(s)+x_d(s))

其中，K_s为位移传感器环节、G_c(s)为控制器环节和G_w(s)为功放环节，则，有：

f_x(s)＝ms²X(s)＝(K_er+K_pr)(X(s)-Θ_x(s))-K_sK_iG_c(s)G_w(s)(X(s)-Θ_x(s)+x_d(s))

从上式可以看出，转子质量不平衡Θ_x(s)不仅会通过控制器产生电磁力，还会通过磁轴承本身产生永磁力，而传感器谐波仅通过控制器产生电磁力，即电磁力中不仅包含同频振动，还包含倍频振动，而永磁力中只包含同频振动，因此谐波振动抑制时需要加以区分；

步骤(2)传感器谐波中同频信号的辨识和补偿

让磁悬浮转子落在保护轴承上，控制磁悬浮转子以1Hz的转速低速旋转，由于此时转速很低，因此同频输出信号中由质量不平衡引起的同频信号很少，此时可认为传感器的同频输出不包含质量不平衡引起的同频信号，均为传感器谐波，从而辨识出传感器谐波的同频信号，转子高速旋转时，加入辨识出的传感器谐波同频信号即完成传感器谐波同频信号的补偿；

步骤(3)基于频域变参数LMS算法的倍频振动抑制

对传感器谐波中的同频信号进行补偿后，传感器谐波只包含倍频信号，传感器谐波通过控制器产生倍频电磁力而不包含永磁力，因此，以倍频电流为控制对象即可对倍频振动进行抑制，将传感器谐波引起的倍频振动通过频域变参数LMS(LeastMeanSquare,LMS)算法进行滤除，谐波振动抑制模块的计算过程如下：

为了抑制倍频振动，采用频域变参数LMS算法进行倍频振动抑制，倍频振动由传感器谐波引起，以与倍频振动相关的正弦信号作为参考输入，即由与传感器倍频频率成分相同的单位幅值的正弦信号相加组成，误差作为基本输入，为了权衡收敛速度和稳态误差这一矛盾，实时地改变块长和步长，以达到更快的收敛速度和更小的稳态误差，根据误差变化情况设计块长和步长的更新算法，从而改善频域LMS算法的收敛性能；

步骤(4)基于陷波器的不平衡振动抑制

完成频域LMS算法进行倍频电流抑制和传感器谐波同频信号的补偿后，只剩下由质量不平衡引起的同频电磁力和永磁力，即不平衡振动力，可利用陷波器抑制不平衡振动力，直接以不平衡振动力为控制目标，将构造出的不平衡振动力作为陷波器的输入，输出反馈至控制器的输入。

2.根据权利要求1所述的一种基于频域变参数LMS算法和陷波器的磁悬浮转子谐波振动抑制方法，其特征在于：所述的步骤(3)块长更新算法为：

设当前块长为N_x，下一块块长为N_x'，定义估计误差e_x,ave(k)记为：

$e_{x,ave}\left(k\right)={\mathrm{\&alpha;e}}_{x,ave}(k-1)+\frac{(1-\mathrm{\&alpha;})}{N_x}\underset{i=0}{\overset{N_x-1}{\&Sigma;}}{e}_x({\mathrm{kN}}_x+i)$

其中，e_x,ave(k-1)是第k-1块的估计误差，α为一个常数且0<α<1，e_x(j)表示第j时刻X通道的基本输入，j＝kN_x,kN_x+1,…,kN_x+N_x-1，块长更新的主要思想是：如果当前块的估计误差小于上一块的估计误差，块长增大；如果当前块的估计误差大于上一块的估计误差，块长减小，然而，由于前后两块的估计误差肯定不完全相等，导致每块的块长都会改变，这样会增加计算量，因此，前后两块进行对比时留有一定变化范围，块长的具体更新算法如下：

若β₁e_x,ave(k-1)<e_x,ave(k)<β₂e_x,ave(k-1)，则视为估计误差基本没变化，N_x'＝N_x；其中，β₁、β₂为常数，且0<β₁<1，β₂>1，用来保证前后两块比较时留有一定的余量，即在当前块与上一块的估计误差相差在一定范围内时认为误差基本不变，不用改变块长；

若e_x,ave(k)≥β₂e_x,ave(k-1)，则视为当前块的估计误差大于上一块的估计误差，此时，N_x'＝N_x/2；

若e_x,ave(k)≤β₁e_x,ave(k-1)，则视为当前块的估计误差小于上一块的估计误差，此时，N_x'＝2N_x；

所述的步骤(3)步长更新算法为：

对于每一块采用不同步长以更好地改善收敛性能，特别地，频域LMS的收敛速度可以在不影响最小均方误差的基础上，通过对每个可调权值赋予不同的步长获得改善，步长更新算法为：

${\mathrm{\&mu;}}_{x,i}\left(k\right)=\frac{{\mathrm{\&mu;}}_0\frac{e_{x,ave}\left(k\right)}{e_{x,ave}\left(k-1\right)}}{\left|\right|U_{x,i}\left(k\right)|{|}^2+\mathrm{\&gamma;}}$

其中，μ_x,i(k)为第k块中第i个步长，U_x,i(k)为第k块中第i个参考输入的频域信号，||·||表示模值，μ₀起到控制步长大小的作用，γ是一个常数，保证步长不会因输入信号能量过低而过大。