光刻机是集成电路生产线中的核心装备，包括光源、光学透镜系统、掩模台、硅片台和对准测量系统等^[1-2]。其中搭载硅片的硅片台采用粗微叠层结构，现有的粗动台采用磁悬浮平面电机驱动，能够有效地解决原有的“直线电机+气浮导轨”驱动方式加速度受限的缺点，可实现高速、高加减速的大行程运动。微动台由音圈电机驱动，在粗动台运动的基础上进行微调，以实现纳米级精度的小行程运动，二者结合可使得硅片台在大行程范围内实现纳米级的高精度运动^[3-4]。

为提升光刻机产率，需要提高磁悬浮平面电机的速度和加速度，这要求其提供较大的电机推力并实现结构的轻量化设计，但会导致磁悬浮平面电机的结构固有频率降低。同时，为提高光刻机的套刻精度，需要磁悬浮平面电机实现更高的运动精度，因此需提高系统的控制带宽。结构固有频率的降低和控制带宽的提升，使得磁悬浮平面电机的柔性结构动力学特性成为制约硅片台系统性能提高的瓶颈。

准确估计磁悬浮平面电机的模态参数并采取相应的控制策略可以抑制其振动。本文中的模态参数估计特指柔性模态参数估计，共包含系统辨识和模态分析2部分。系统辨识是指通过参数优化算法得到系统参数化频率响应函数(frequency response function，FRF)的过程；而模态分析是指根据系统的参数化频率响应函数得到系统模态参数的过程，其中系统待估计的模态参数包括固有频率、模态阻尼系数、特征振型和各线圈的柔性模态力系数。

磁悬浮平面电机的系统辨识问题为非凸优化问题，难以求得全局最优解，Guillaume等^[5]将原优化问题进行线性化近似，并提出非迭代的正规方程法和降维正规方程法求解待优化参数。为了减小线性化近似所引入的误差，Sanathanan等^[6]提出在每次迭代中将权重系数设为上次迭代所得到的分母，然而该方法仍无法消除目标函数线性化带来的误差。另有学者利用Gauss-Newton法对原非线性问题进行直接迭代求解，并将权重系数设为相对误差权重^[7-8]，然而该方法会影响参数化频率响应函数在低频共振点处的估计精度，进而影响后续的模态参数估计精度。

系统辨识得到的系统参数化模型中，既包含系统真实存在的物理模态，也包含由于系统测量噪声或者模型参数偏差等所引入的数学模态。模态分析的目标是从模态集合中对物理模态进行提取和分类，并进行物理模态参数估计。现有的模态分析方法中最简单直接的是稳态图(stabilization diagram)法，它基于如下的经验观测：物理模态的各模态参数随模型阶数增加而保持不变，而数学模态的各模态参数随模型阶数增加而发生变化，进而通过人为设定各模态参数的误差阈值确定系统的物理模态^[9-11]。稳态图法简单直观，但需要人工设定阈值，比较依赖经验。为了减少需要人工设定的参数数量，实现自动化的模态分析，多阶段聚类方法^[12-16]被大量学者研究，并得到广泛应用，然而该类方法过程复杂，且现有的模态验证标准(modal assurance criterion，MAC)假设系统质量分布均匀且测量点数量较多^[17]，不适用于质量分布不均匀且测量点少的磁悬浮平面电机。

为了解决磁悬浮平面电机模态参数的准确估计问题，本文提出了一种基于密度聚类的模态参数估计方法。该方法首先通过一种两步迭代的系统辨识算法，得到系统的参数化频率响应函数。然后，利用DBSCAN聚类算法对系统进行模态分析，去除不稳定的数学模态点，并对物理模态点进行基于概率密度的离群点剔除，最终获得系统稳定的物理模态参数。通过仿真对所提方法的模态参数估计精度进行了分析，并通过实验验证了所提方法的有效性。

1 系统辨识算法 1.1 问题描述

磁悬浮平面电机为包含刚体模态和柔性模态的多输入多输出(multi-input multi-output，MIMO)系统，且z方向的柔性特性最为突出。系统输入为各线圈电流，对于模态参数估计而言，系统输出为动子各点z方向位移。单线圈电流输入、单点z方向位移输出所构成的单输入单输出(single-input single-output，SISO)系统传递函数为

$ g(s)=\frac{k_z}{m s^2}+\sum\limits_{n=1}^N \frac{\phi_n k_n}{s^2+2 \xi_n \omega_n s+\omega_n^2} . $

(1)

其中：m为动子质量，k_z为单线圈z向推力系数，ω_n为第n阶固有频率，ξ_n为第n阶模态阻尼系数，ϕ_n为第n阶特征振型，k_n为单线圈第n阶柔性模态力系数，N为建模所考虑的系统柔性模态阶次。由于平面电机系统具有统一的固有频率和模态阻尼系数，因此其MIMO系统频率响应函数矩阵中各项应具有相同的分母多项式^[18]，其参数化结构可表示为

$ \begin{gathered} \hat{\boldsymbol{G}}(\omega, \boldsymbol{\theta})=\frac{\hat{\boldsymbol{N}}(\omega, \boldsymbol{\beta})}{\hat{D}(\omega, \boldsymbol{\alpha})}= \\ \frac{1}{\sum\limits_{q=0}^{n_{\mathrm{a}}} \alpha_q \mathit{\Omega}_q(\omega)}\left(\sum\limits_{q=0}^{n_{\mathrm{b}}} \beta_{i j}^q \Omega_q(\omega)\right)_{N_{\mathrm{o}} \times N_{\mathrm{i}}} . \end{gathered} $

(2)

其中：N_o为输出维数，N_i为输入维数，待优化参数θ=[β₁^T, β₂^T, …, β_(i-1)Ni+j^T, …, β_NoNi^T, α^T]^T，β_(i-1)Ni+j=[β_ij⁰, β_ij¹, …, β_ij^n_b]^T，α=[α₀, α₁, …, α_na]^T，i=1, 2, …, N_o, j=1, 2, …, N_i, n_a为分母多项式最高阶次，n_b为分子多项式最高阶次，α_q为频率响应函数矩阵分母多项式的第q阶系数，β_ij^q为频率响应函数矩阵第i行第j列元素分子多项式的第q阶系数，Ω_q(ω)=(jω)^q为多项式基函数。

系统辨识问题的目标函数为

$ L_{\mathrm{NLS}}(\boldsymbol{\theta})=\sum\limits_{r=1}^{N_{\mathrm{o}} N_{\mathrm{i}}} \sum\limits_{f=1}^{N_{\mathrm{f}}}\left[W_r\left(\omega_f\right)\left(\hat{G}_r\left(\omega_f, \boldsymbol{\theta}\right)-G_r\left(\omega_f\right)\right)\right]^2 . $

(3)

其中：G_r(ω_f)为实验辨识得到的系统非参数化频率响应函数，下标r表示频率响应函数矩阵中第r个元素，ω_f为第f个频率采样点的频率，N_f为频率采样点数目，W_r(ω_f)为权重系数。

1.2 两步迭代系统辨识

目标函数式(3)与待优化参数θ之间呈非线性关系，且待优化参数α的存在使得原优化问题为非凸优化问题，因此基于梯度的迭代求解算法对参数初始化状态敏感度高，即不同的参数初始化状态经迭代算法后将收敛到不同的局部最优值，难以求得全局最优参数。为此，本文采取一种两步迭代式的系统辨识算法，来获取系统参数化频率响应函数的最优参数。

首先将式(3)线性化：

$ \begin{array}{*{20}{c}} & L_{\mathrm{LS}}(\boldsymbol{\theta})= \sum\limits_{r=1}^{N_{\mathrm{o}} N_{\mathrm{i}}} \sum\limits_{f=1}^{N_{\mathrm{f}}}\left[W _ { r } ( \omega _ { f } ) \left(\hat{N}_r\left(\omega_f, \boldsymbol{\beta}_r\right)-\right.\right. \\ &\left.\left.\hat{D}\left(\omega_f, \boldsymbol{\alpha}\right) G_r\left(\omega_f\right)\right)\right]^2 . \end{array} $

(4)

线性化误差可以表示为

$ \begin{gathered} {e_{\mathrm{LS}}(\boldsymbol{\theta})=\boldsymbol{J} \boldsymbol{\theta}=} \\ {\left[\begin{array}{ccccccc} \boldsymbol{P}_1 & 0 & \cdots & 0 & \cdots & 0 & \boldsymbol{Q}_1 \\ 0 & \boldsymbol{P}_2 & \cdots & 0 & \cdots & 0 & \boldsymbol{Q}_2 \\ \vdots & \vdots & & \vdots & & \vdots & \vdots \\ 0 & 0 & \cdots & \boldsymbol{P}_r & \cdots & 0 & \boldsymbol{Q}_r \\ \vdots & \vdots & & \vdots & & \vdots & \vdots \\ 0 & 0 & \cdots & 0 & \cdots & \boldsymbol{P}_{N_{\mathrm{o}} N_{\mathrm{i}}} & \boldsymbol{Q}_{N_{\mathrm{o}} N_{\mathrm{i}}} \end{array}\right]\left[\begin{array}{c} \boldsymbol{\beta}_1 \\ \boldsymbol{\beta}_2 \\ \vdots \\ \boldsymbol{\beta}_r \\ \vdots \\ \boldsymbol{\beta}_{N_{\mathrm{o}} N_{\mathrm{i}}} \\ \boldsymbol{\alpha} \end{array}\right] .} \end{gathered} $

(5)

其中J为误差对于待优化参数的Jacobi矩阵，且

$ \boldsymbol{P}_r=\left[\begin{array}{c} W_r\left(\omega_1\right)\left[\mathit{\Omega}_0\left(\omega_1\right), \mathit{\Omega}_1\left(\omega_1\right), \cdots, \mathit{\Omega}_{n_{\mathrm{b}}}\left(\omega_1\right)\right] \\ W_r\left(\omega_2\right)\left[\mathit{\Omega}_0\left(\omega_2\right), \mathit{\Omega}_1\left(\omega_2\right), \cdots, \mathit{\Omega}_{n_{\mathrm{b}}}\left(\omega_2\right)\right] \\ \vdots \\ W_r\left(\omega_{N_{\mathrm{f}}}\right)\left[\mathit{\Omega}_0\left(\omega_{N_{\mathrm{f}}}\right), \mathit{\Omega}_1\left(\omega_{N_{\mathrm{f}}}\right), \cdots, \mathit{\Omega}_{n_{\mathrm{b}}}\left(\omega_{N_{\mathrm{f}}}\right)\right] \end{array}\right] . $

(6)

$ \boldsymbol{Q}_r=\left[\begin{array}{c} -W_r\left(\omega_1\right) G_r\left(\omega_1\right)\left[\mathit{\Omega}_0\left(\omega_1\right), \mathit{\Omega}_1\left(\omega_1\right), \cdots, \mathit{\Omega}_{n_{\mathrm{a}}}\left(\omega_1\right)\right] \\ -W_r\left(\omega_2\right) G_r\left(\omega_2\right)\left[\mathit{\Omega}_0\left(\omega_2\right), \mathit{\Omega}_1\left(\omega_2\right), \cdots, \mathit{\Omega}_{n_{\mathrm{a}}}\left(\omega_2\right)\right] \\ \vdots \\ -W_r\left(\omega_{N_{\mathrm{f}}}\right) G_r\left(\omega_{N_{\mathrm{f}}}\right)\left[\mathit{\Omega}_0\left(\omega_{N_{\mathrm{f}}}\right), \mathit{\Omega}_1\left(\omega_{N_{\mathrm{f}}}\right), \cdots, \mathit{\Omega}_{n_{\mathrm{a}}}\left(\omega_{N_{\mathrm{f}}}\right)\right] \end{array}\right] . $

(7)

对待优化参数进行随机初始化，并应用SK算法^[6]选取每次迭代的权重系数：

$ W_r^{\langle p\rangle}\left(\omega_f\right)= \begin{cases}1, & p=1 ; \\ \frac{1}{\hat{D}\left(\omega_f, \boldsymbol{\alpha}^{\langle p-1\rangle}\right)}, & p \neq 1 .\end{cases} $

(8)

其中p为迭代次数，并应用Gauss-Newton法对待定参数进行优化，其参数更新律为

$ \boldsymbol{\theta}^{\langle p+1\rangle}=\boldsymbol{\theta}^{\langle p\rangle}-\lambda\left(\operatorname{Re}\left(\boldsymbol{J}^{\mathrm{H}} \boldsymbol{J}\right)\right)^{-1} \operatorname{Re}\left(\boldsymbol{J}^{\mathrm{H}}\right) \operatorname{Re}\left(\boldsymbol{e}_{\mathrm{LS}}\left(\boldsymbol{\theta}^{\langle p\rangle}\right)\right) . $

(9)

其中：λ为迭代步长，H表示Jacobi矩阵的共轭转置。由于SK算法对原问题的非线性目标函数进行了线性化处理，其所得的参数并非最优。但SK算法所得的参数可以为原非凸优化问题提供一个理想的初始化状态，因此可用第1步迭代中所得到的最优参数进行参数初始化：

$ \boldsymbol{\theta}^{\langle 0\rangle}=\underset{\boldsymbol{\theta}}{\operatorname{argmin}} L_{\mathrm{NLS}}(\boldsymbol{\theta}) . $

(10)

然后对原非凸优化问题进行第2步迭代求解，以获取原问题的最优参数。原问题的非线性误差关于待优化参数θ的梯度可以表示为

$ \begin{gathered} \nabla{\boldsymbol{e}_{\mathrm{NLS}}(\boldsymbol{\theta})=\boldsymbol{J}_{\mathrm{NL}}=} \\ {\left[\begin{array}{ccccccc} \boldsymbol{U}_1 & 0 & \cdots & 0 & \cdots & 0 & \boldsymbol{V}_1 \\ 0 & \boldsymbol{U}_2 & \cdots & 0 & \cdots & 0 & \boldsymbol{V}_2 \\ \vdots & \vdots & & \vdots & & \vdots & \vdots \\ 0 & 0 & \cdots & \boldsymbol{U}_r & \cdots & 0 & \boldsymbol{V}_r \\ \vdots & \vdots & & \vdots & & \vdots & \vdots \\ 0 & 0 & \cdots & 0 & \cdots & \boldsymbol{U}_{N_{\mathrm{o}} N_{\mathrm{i}}} & \boldsymbol{V}_{N_{\mathrm{o}} N_{\mathrm{i}}} \end{array}\right] .} \end{gathered} $

(11)

其中J_NL为原问题非线性误差对于待优化参数的Jacobi矩阵，且

$ \boldsymbol{U}_r=\left[\begin{array}{c} W_r\left(\omega_1\right) / \hat{D}\left(\omega_1, \boldsymbol{\alpha}\right)\left[\mathit{\Omega}_0\left(\omega_1\right), \mathit{\Omega}_1\left(\omega_1\right), \cdots, \mathit{\Omega}_{n_{\mathrm{b}}}\left(\omega_1\right)\right] \\ W_r\left(\omega_2\right) / \hat{D}\left(\omega_2, \boldsymbol{\alpha}\right)\left[\mathit{\Omega}_0\left(\omega_2\right), \mathit{\Omega}_1\left(\omega_2\right), \cdots, \mathit{\Omega}_{n_{\mathrm{b}}}\left(\omega_2\right)\right] \\ \vdots \\ W_r\left(\omega_{N_{\mathrm{f}}}\right) / \hat{D}\left(\omega_{N_{\mathrm{f}}}, \boldsymbol{\alpha}\right)\left[\mathit{\Omega}_0\left(\omega_{N_{\mathrm{f}}}\right), \mathit{\Omega}_1\left(\omega_{N_{\mathrm{f}}}\right), \cdots, \mathit{\Omega}_{n_{\mathrm{b}}}\left(\omega_{N_{\mathrm{f}}}\right)\right] \end{array}\right], $

(12)

$ \boldsymbol{V}_r=\left[\begin{array}{c} -W_r\left(\omega_1\right) \hat{N}_r\left(\omega_1, \boldsymbol{\beta}_r\right) /\left(\hat{D}\left(\omega_1, \boldsymbol{\alpha}\right)\right)^2\left[\mathit{\Omega}_0\left(\omega_1\right), \mathit{\Omega}_1\left(\omega_1\right), \cdots, \mathit{\Omega}_{n_{\mathrm{a}}}\left(\omega_1\right)\right] \\ -W_r\left(\omega_2\right) \hat{N}_r\left(\omega_2, \boldsymbol{\beta}_r\right) /\left(\hat{D}\left(\omega_2, \boldsymbol{\alpha}\right)\right)^2\left[\mathit{\Omega}_0\left(\omega_2\right), \mathit{\Omega}_1\left(\omega_2\right), \cdots, \mathit{\Omega}_{n_{\mathrm{a}}}\left(\omega_2\right)\right] \\ \vdots \\ -W_r\left(\omega_{N_f}\right) \hat{N}_r\left(\omega_{N_{\mathrm{f}}}, \boldsymbol{\beta}_r\right) /\left(\hat{D}\left(\omega_{N_{\mathrm{f}}}, \boldsymbol{\alpha}\right)\right)^2\left[\mathit{\Omega}_0\left(\omega_{N_{\mathrm{f}}}\right), \mathit{\Omega}_1\left(\omega_{N_{\mathrm{f}}}\right), \cdots, \mathit{\Omega}_{n_{\mathrm{a}}}\left(\omega_{N_{\mathrm{f}}}\right)\right] \end{array}\right] . $

(13)

此外，由于第1步迭代所得参数已接近全局最优参数，为了避免参数化频率响应函数在低频共振段估计误差增加而影响后续的模态参数估计精度，第2步迭代中不同频率点的权重系数均取为W_r(ω_f)=1。并应用Gauss-Newton法对待定参数进行优化，其参数更新律为

$ \begin{array}{*{20}{c}} \boldsymbol{\theta}^{\langle p+1\rangle}=\boldsymbol{\theta}^{\langle p\rangle}-\lambda\left(\operatorname{Re}\left(\boldsymbol{J}_{\mathrm{NL}}^{\mathrm{H}} \boldsymbol{J}_{\mathrm{NL}}\right)\right)^{-1} \cdot\\ \operatorname{Re}\left(\boldsymbol{J}_{\mathrm{NL}}^{\mathrm{H}}\right) \operatorname{Re}\left(\boldsymbol{e}_{\mathrm{NLS}}\left(\boldsymbol{\theta}^{\langle p\rangle}\right)\right). \end{array} $

(14)

2 基于密度聚类的模态分析 2.1 密度聚类算法

1.2节的系统辨识算法可以得到系统不同阶次的参数化模型：

$ \hat{G}_r(s)=\frac{\sum\limits_{q=0}^{n_{\mathrm{b}}} \beta_r^q s^q}{\sum\limits_{q=0}^{n_{\mathrm{a}}} \alpha_q s^q} . $

(15)

对其进行部分分式展开即可计算待估计模态参数的原始估计结果：

$ \hat{G}_r(s)=\frac{\hat{R}_r^0}{s^2}+\sum\limits_{n=1}^{n_a / 2-1} \frac{\hat{R}_r^n}{s^2+2 \hat{\xi}_n \hat{\omega}_n s+\hat{\omega}_n^2} . $

(16)

其中：$ \hat{\omega}_n$为第n阶固有频率估计值，$ \hat{\xi}_n$为第n阶模态阻尼系数估计值，$ \hat{R}_r^n$为第n阶特征振型与第n阶柔性模态力系数的乘积估计值。原始估计结果中既包含稳定的物理模态参数，也包含由于系统测量噪声或者模型参数偏差等所引入的数学模态参数，2类模态极点在极点复平面的密度分布不同，因此可利用密度聚类算法进行物理模态提取和分类。具体执行步骤如下。

首先对原始估计结果中的模态极点进行硬性准则筛选，即物理模态极点需满足实部小于0的条件，不满足该条件的模态极点均为数学模态极点，在该步骤中被剔除。此外由于物理模态极点均以共轭形式存在，因此仅保留第二象限极点进行后续分析。

然后利用DBSCAN聚类算法^[19]对剩余的模态极点进行模态分析。DBSCAN聚类算法是经典的密度聚类算法，其待定超参数包括邻域半径ε和数量阈值M。该算法首先将样本集中ε-邻域内点的数量大于或者等于M的点定义为核心点，然后从核心点出发，不断向密度可达的区域扩张，最终得到一个包含核心点和边界点的最大化密度相连集合，该集合记为样本集中的一个簇，然后重复上述过程直到样本集中的所有点被划分完毕。该算法的优点在于可以对任意形状的稠密数据集进行聚类，对算法初始值和数据集中的异常点不敏感，且适用于模态极点这类低维数据。

针对模态分析问题，本文给出ε和M的确定方法如下：

1) 首先计算每个模态极点与其距离最近的第k个点的距离，记为每个极点的k-dist参数，极点之间的距离度量采用Euclid距离。

2) 将所有极点的k-dist参数从高到低排序，画出极点编号与k-dist参数的关系图，即k-dist图。

3) 根据系统参数化模型阶次与所要估计的柔性模态数量，分别确定ε和M：

$ \begin{aligned} \varepsilon & =k-\operatorname{dist}\left(N_{\mathrm{p}}-\left(\max \left(n_{\mathrm{a}}\right)-\right.\right. \\ & \left.\left.\min \left(n_{\mathrm{a}}\right)+1\right) \times N+1\right), \end{aligned} $

(17)

$ M=k+1. $

(18)

其中：N_p为硬性准则筛选后的模态极点总数，所要估计的柔性模态数量为N，max(n_a)和min(n_a)分别表示参数化模型中分母多项式最高阶次的最大值和最小值。k的选取虽然会影响超参数ε和M的取值，但由于不同k所得的k-dist图趋势基本一致，因此对聚类结果影响不大。

2.2 基于概率分布的离群点剔除

DBSCAN聚类算法可以得到属于各阶柔性模态的物理模态点集合，还需要在每个集合中选择合适的代表元素以确定最终的模态参数估计结果。现有文献大多直接取每个集合中的参数均值作为最终的模态参数。本文研究发现，不同阶模型得到的同阶模态参数符合正态分布的概率密度曲线，且参数集中在区间[μ-σ, μ+σ](该结论将在3.3节中验证)，因此该区间之外的点可视为离群点剔除，并对各集合中的剩余点取均值作为最终的模态参数估计值。

2.3 分子乘积矩阵奇异值分解

经过模态分析得到系统的传递函数矩阵可表示为

$ \hat{\boldsymbol{G}}(s)=\frac{\hat{\boldsymbol{R}}^0}{s^2}+\sum\limits_{n=1}^N \frac{\hat{\boldsymbol{R}}^n}{s^2+2 \hat{\xi}_n \hat{\omega}_n s+\hat{\omega}_n^2} . $

(19)

其中$ \hat{\boldsymbol{R}}^n \in \mathbb{R}^{N_{\rm{o}} \times N_{\rm{i}}}$为第n阶特征振型与第n阶柔性模态力系数的乘积矩阵，对其进行主奇异值分解即可得到各阶特征振型与各阶柔性模态力系数：

$ \left\{\begin{array}{l} \hat{\boldsymbol{R}}^n=\boldsymbol{U}_n \boldsymbol{S}_n \boldsymbol{V}_n^{\mathrm{T}}, \\ \hat{\boldsymbol{\phi}}_n=\left[\boldsymbol{U}_n\right]^1, \\ \hat{\boldsymbol{k}}_n=\left[\boldsymbol{S}_n\right]_1^1\left[\boldsymbol{V}_n^{\mathrm{T}}\right]_1 . \end{array}\right. $

(20)

其中：[A]_i和[A]^j分别表示矩阵A的第i行和第j列向量，[A]_i^j表示矩阵A的第i行第j列元素。

3 仿真验证 3.1 两步迭代系统辨识算法验证

仿真采用的多自由度振动系统模型为：

$ G_1(s)=\frac{1}{s^2}+\sum\limits_{n=1}^6 \frac{R_1^n}{s^2+2 \xi_n \omega_n s+\omega_n^2}, $

(21)

$ G_2(s)=\frac{1}{s^2}+\sum\limits_{n=1}^6 \frac{R_2^n}{s^2+2 \xi_n \omega_n s+\omega_n^2} . $

(22)

其中模态参数如表 1所示，该系统包含刚体模态与6阶柔性模态，并在系统输出中加入均值为0，标准差为1 μm的Gauss噪声以模拟测量噪声，系统辨识算法与模态分析方法仅对模型中的前4阶柔性模态参数进行估计。

估计出前4阶柔性模态参数所需的最低模型阶次为10，因此n_a的最小值取为10。n_a的最大值取为40，n_b=n_a-2。频率采样点范围为10~350 Hz，包含了前4阶柔性模态固有频率。

以10阶参数化模型为例，图 1给出了两步迭代系统辨识算法的误差收敛曲线，从图 1a中可以看到，在第1步SK算法迭代过程中，式(3)的原问题非线性误差不会单调收敛，而在第28次迭代时达到最优值。将该次迭代后的参数作为第2步迭代的初始参数，利用式(14)进行第2步迭代，从图 1b中可以看到，原问题非线性误差可单调收敛。最终两步迭代辨识算法所得到的10阶参数化模型如图 2所示。

利用两步迭代辨识算法对10~40阶参数化模型分别进行参数优化，并根据式(16)对其结果进行部分分式展开，即可得到待估计模态参数的原始估计结果如图 3和4所示。可以看到，前4阶柔性模态参数的原始估计结果均在真值附近，验证了所提两步迭代系统辨识算法的有效性。

3.2 基于密度聚类的模态分析方法验证

本节利用所提出的基于密度聚类的模态分析方法对原始估计结果进行物理模态的提取和分类。首先对原始估计结果中的模态极点进行硬性准则筛选，以剔除部分数学模态极点，剔除前后的模态极点分别如图 5和6所示，其中真实极点位于半径为固有频率的圆上。

从剔除后的剩余极点分布可以发现，稳定的物理模态点集中在真实极点所在圆附近的高密度区域，而不稳定的数学模态点则分布在复平面的低密度区域，因此可利用密度聚类算法进行物理模态提取和分类。

然后计算筛选后每个模态极点的k-dist参数，其中，并画出不同k值的k-dist图如图 7所示，可以看到不同k值的k-dist曲线趋势一致。根据式(17)和(18)分别确定ε=0.42和M=11，利用DBSCAN聚类算法对筛选后的模态极点进行物理模态极点的提取和分类，图 8a给出了前4阶柔性模态极点的聚类结果，其中η₁~η₄表示前4阶柔性模态。特别地，图 8b给出了第1阶柔性模态极点的聚类结果，以验证密度聚类算法可以有效地实现物理模态极点的提取和分类。

3.3 模态参数估计误差

图 9给出了前4阶模态阻尼系数的概率分布，可以发现，不同阶模型得到的同阶模态阻尼系数符合正态分布的概率密度曲线且参数集中在[μ-σ, μ+σ]区间，因此该区间之外的点可视为离群点剔除。其余模态参数同样符合正态分布的概率密度曲线且集中在[μ-σ, μ+σ]区间，限于篇幅不再一一列出。

离群点剔除前后的模态阻尼系数分布对比如图 10所示，离群点剔除后的模态阻尼系数更集中地分布在真值附近。对各集合中的剩余点取均值作为最终的模态参数估计结果，并计算离群点剔除前后各模态参数的相对估计误差，分别如表 2和3所示，离群点剔除后大部分模态参数(表 3中加粗项)的估计误差相比剔除前有所下降，且各模态参数的相对估计误差均小于3%，证明所提方法可以实现模态参数的准确估计。

4 实验验证 4.1 单线圈激励闭环辨识

实验在图 11的磁悬浮平面电机上进行，定子为二维Halbach永磁阵列，动子为线圈阵列，结构原理如图 12所示。在模态参数估计实验中系统输入为12个线圈电流，系统输出为8个位于线圈阵列角部的电涡流传感器信号，用于测量各点的z向位移。由于磁悬浮平面电机为开环不稳定系统，因此辨识实验在闭环条件下进行，利用PID反馈控制器实现线圈阵列的闭环稳定悬浮，每次实验时在单线圈中加入均值为0、标准差为3 A的Gauss噪声激励电流i_e，并记录该线圈的总输入电流i_a和系统8个测量点的输出z，通过频谱估计计算2个中间模型：

$ g_{i_{\mathrm{a} i}{ }_{\mathrm{e}}}(\omega)=\frac{\mathit{\Phi}_{i_{\mathrm{a}} i_{\mathrm{e}}}}{\mathit{\Phi}_{i_{\mathrm{e}}}}, $

(23)

$ \boldsymbol{G}_{z i_{\mathrm{e}}}(\omega)=\frac{\mathit{\pmb{\Phi}}_{z i_{\mathrm{e}}}}{\mathit{\Phi}_{i_{\mathrm{e}}}} . $

(24)

其中：Φ_ie为i_e的频谱估计，Φ_iaie为i_a与i_e的互谱估计，Φ_zie为z与i_e的互谱估计。从而可得到系统的频域辨识模型为

$ \boldsymbol{G}^1(\omega)=\frac{\boldsymbol{G}_{z i_{\mathrm{e}}}(\omega)}{g_{i_\mathrm{a} i_{\mathrm{e}}}(\omega)} . $

(25)

其中G¹(ω)表示MIMO系统频率响应函数矩阵的第1列。对12个线圈分别进行上述闭环辨识实验，即可得到MIMO系统频率响应函数矩阵。

4.2 模态参数估计

同样只对系统前4阶柔性模态参数进行估计，n_a的最小值取为10，n_a的最大值取为40，n_b=n_a-2。根据频谱中模态共振峰对应的频率位置，选取频率采样点范围为100~500 Hz，包含了前4阶柔性模态固有频率。图 13给出了两步迭代辨识算法对系统频率响应函数的估计结果，其中8个子图分别对应8个电涡流传感器的输出通道。从估计结果可以看出，所提出的两步迭代系统辨识算法可以有效地实现系统频率响应函数的参数化。

然后对原始估计结果中的模态极点进行硬性准则筛选，并计算筛选后每个模态极点的k-dist参数，根据式(17)和(18)确定ε=49和M=11，利用DBSCAN聚类算法对筛选后的模态极点进行物理模态极点的提取和分类，图 14给出了前4阶柔性模态极点的聚类结果，进一步验证了密度聚类算法可以有效地实现对物理模态极点的提取和分类。所得到的前4阶柔性模态固有频率分别为151.64、241.02、324.76和492.05 Hz。

图 15给出了前4阶柔性模态的模态阻尼系数的概率分布，可以看到，不同阶模型得到的同阶模态阻尼系数符合正态分布的概率密度曲线且参数集中在[μ-σ, μ+σ]区间，该区间之外的点可视为离群点剔除。其余模态参数同样符合正态分布的概率密度曲线且集中在[μ-σ, μ+σ]区间，限于篇幅不再一一列出。对各集合中的剩余点取均值作为最终的模态参数估计结果。

上述模态参数估计过程可得到系统8×12维的传递函数矩阵的第一列元素，根据式(19)所示的传递函数矩阵同分母特性，在进行其余线圈激励的闭环辨识实验时，可复用此次估计的分母参数(包括固有频率和模态阻尼系数)，而只估计频率响应函数的分子参数，从而将原模态参数估计问题转化为标准最小二乘问题，以提高计算效率。

对系统8×12维的传递函数矩阵进行式(20)的主奇异值分解，即可得到系统在8个测量点的特征振型和12个线圈的柔性模态力系数。保持线圈阵列在全局坐标系的y方向坐标p_y不变，x方向坐标p_x从36.5 mm至121.3 mm每隔$ \sqrt{2} \tau / 4$进行一次特征振型和柔性模态力系数估计，其中磁钢极距τ=60 mm。图 16和17分别给出了p_x为57.7和100.1 mm处8个测量点的特征振型估计值与自然邻点差值^[20]结果，其中x和y为线圈阵列坐标系的坐标，曲面为自然邻点插值结果。需要说明的是，由于测量点均位于线圈阵列的角部，因此自然邻点差值所得到的各阶特征振型在线圈阵列中心会与实际振型有所偏差，该插值结果主要为了验证线圈阵列的特征振型在不同位置保持不变。图 18给出了单线圈柔性模态力系数在5个位置的估计值，其变化规律符合周期为$ \sqrt{2} \tau$的正弦曲线。

更进一步地，图 19给出了12个线圈前四阶柔性模态力系数随p_x变化情况，各线圈的柔性模态力系数均会随线圈阵列在全局坐标系的位置p_x呈周期性变化，且周期均为$ \sqrt{2} \tau$。

5 结论

针对磁悬浮平面电机的模态参数估计问题，本文提出了一种基于密度聚类的模态参数估计方法。首先通过两步迭代的系统辨识算法，获取了系统不同阶的参数化频率响应函数。在此基础上利用DBSCAN聚类算法对系统进行模态分析，并对提取出的物理模态点进行基于正态分布的离群点剔除，获取了系统最终的模态参数。仿真和实验结果验证了所提方法的有效性。模态参数的估计结果为后续设计可抑制系统振动的电流分配算法奠定了基础。