舵机是气动伺服弹性系统中的执行机构，其将舵面的控制指令信号转化为舵面的运动，从而驱动舵面偏转产生控制力矩^[1]。目前，市面上许多公司的舵机伺服系统解决方案，其舵机本身频响特性较好，并给出完整的特性参数，提供配套的控制系统，可针对舵机实际工作状态调整参数，满足工作要求。但是，对于大部分小型民用无人机来说，由于受到设计空间、质量、设计成本等方面的限制，无法选用频响特性好的驱动系统，而是选用航模常用的普通小型舵机。
小型无人机选用舵机时首要考虑舵机扭力、转速等是否达到要求，舵机频响特性考虑较少。但是，一些大展弦比的无人机为避免飞机破坏，会加入阵风减缓、颤振主动抑制等控制环节。这些控制环节不仅对舵机的带宽和响应速度有较高的要求，还需要获得舵机准确的数学模型以进行主动控制。
一方面，普通小型舵机内部元件组成复杂，受元件实际特性的影响，单纯基于理论模型的数学推导无法获得精确的舵机传递函数^[2]。另一方面，普通小型舵机的特性参数不全面，缺少频响特性测试报告。为此，在选用舵机时，有必要对拟选舵机进行频响特性测试，获取舵机的频响特性，并辨识舵机的传递函数。
目前，舵机频响特性测试多针对单一舵机，对市面上常用的小型舵机缺少通用的舵机频响测试平台^[2-6]。舵机传递函数的辨识方法通常采用曲线拟合法，如Levy法^[4]，通过2个待定实系数的多项式的比值描述该系统，利用最小二乘拟合频响曲线得到多项式系数。曲线拟合的方法存在一定的局限性^[7]。一方面，由于曲线拟合方法利用数值迭代，最终结果可能收敛于局部最优解而非全局最优解；另一方面，当系统的阶数很高时，系统对多项式的零点极点更加敏感，鲁棒性变差。
子空间辨识方法是基于控制理论、线性代数和统计学发展而来的系统辨识方法，是一个“输入-状态-输出”的辨识过程^[8]。经过不断的发展，子空间辨识发展出了N4SID、MOESP、CVA^[9]等在内的诸多方法且应用广泛。与传统的曲线拟合方法相比，子空间辨识无需进行非线性和迭代优化^[10-11]，而且子空间辨识对于高阶系统具有很好的鲁棒性^[12]。
本文考虑舵机工作过程中受到的舵面惯性载荷和气动载荷，设计了一种舵机频响特性测试平台，并给出了舵机频响特性测试与系统辨识的方法。验证了测试平台的有效性，测试了舵机在有/无气动载荷下的频响特性差异。根据舵机实测的频响曲线，利用子空间辨识方法辨识并构建了舵机的数学模型，为舵机选型和控制设计提供参考。
1 舵机频响特性测试平台设计 舵机工作状态下的负载包括舵面惯性载荷和气动载荷，为模拟舵机实际工作中受到的载荷，设计了舵机频响特性测试系统。如图 1、图 2所示，该系统主要由舵机测试平台、测控平台两部分组成。

图 1 舵机频响特性测试系统框架 Fig. 1 Framework of actuator frequency response characteristic test system

图选项

图 2 舵机频响特性测试系统实物图 Fig. 2 Photo of actuator frequency response characteristic test system

图选项

舵机测试平台由待测舵机安装模块和舵面载荷模拟模块组成，分别紧固于直线滑轨上，如图 3所示。待测舵机安装模块包括待测舵机及其安装夹具(见图 4)，该舵机安装夹具可调整待测舵机转轴的位置，方便不同尺寸、不同类型舵机的安装；另外，通过调整待测舵机转轴高度，使舵机摇臂、舵机舵面连杆、舵面模拟摇臂构成一套平行四边形平面连杆机构。

图 3 舵机测试平台 Fig. 3 Actuator test platform

图选项

图 4 待测舵机安装模块放大图 Fig. 4 Enlarged view of installation module of actuator to be tested

图选项

舵面载荷模拟模块包括气动载荷模拟扭杆和惯性载荷模拟摇臂等主要部件，在舵面模拟摇臂的转轴处安装角位移传感器，测量摇臂转角，其具体结构如图 5所示。

图 5 舵面载荷模拟模块放大图 Fig. 5 Enlarged view of control surface load simulation module

图选项

1.1 舵面气动载荷模拟 舵面的铰链力矩，即舵面所受气动载荷对转轴处的力矩为

(1)

式中：ρ为大气密度；V为来流速度；S_e为舵面参考面积；C_{he, δ}为铰链力矩系数；δ为舵面偏角。
在本文所述的舵面偏角范围内，舵面气动力满足线性小扰动假设。定义舵面的铰链刚度K_{he, δ}为

(2)

则舵面的铰链力矩是关于舵面偏角δ的函数:H_e(δ)=K_{he, δ}δ。
本文所述舵机测试平台通过扭杆模拟舵面的铰链力矩，为便于安装，扭杆截面形状选为矩形。矩形截面有2个参数：矩形高(h)和矩形宽(b)，该矩形截面的极惯性矩为

(3)

式中：β为矩形截面系数，为简化设计参数，取扭杆截面形状为正方形，即h/b=1，此时β=0.141。
设矩形截面扭杆一端固支，工作长度为L_beam(见图 6)，扭杆自由端产生的扭矩为

图 6 矩形截面扭杆气动力模拟装置原理图 Fig. 6 Schematic diagram of aerodynamic simulation device by a rectangular-section torsion bar

图选项

(4)

式中：G=E/[2(1+μ)]为材料的剪切模量，E为材料的弹性模量，μ为泊松比。
由式(1)和式(4)有

(5)

式中:L_beam和h为正方形截面扭杆设计的待定参数，通过计算舵面偏转产生的铰链力矩，设计合适的扭杆尺寸，从而实现对舵面铰链力矩的模拟。
1.2 舵面惯性载荷模拟 如图 7所示，舵面的惯性载荷指舵面的惯性力在转轴处产生的力矩。为模拟舵面的转动惯量，设舵面的质量为m_g，舵面质心到转轴的距离为L_g。

图 7 舵面转动惯量示意图 Fig. 7 Schematic diagram of rotational interia of control surface

图选项

通过在舵面模拟摇臂上增加配重并调整位置(见图 8)，使得配重的质量m_w=m_g，配重的位置满足L_w=L_g时，该舵面模拟摇臂的转动惯量即等于舵面本身的转动惯量。

图 8 舵面惯性载荷模拟装置原理图 Fig. 8 Schematic diagram of simulation device of rotational inertial loads of control surface

图选项

2 频响函数测试方法 获得幅值和相位差的方法主要有傅里叶变换谱分析法、最小二乘法^[13]和相关分析法^[5]。对于步进正弦扫频测试，在单一频率点上有充足的数据，考虑到傅里叶变换谱分析法对噪声敏感，适合使用最小二乘法或者相关分析法。本文采用最小二乘法^[13]计算系统在每个频率点上的幅值和相位差。
从计算机输出的扫频信号中截取样本个数为N、频率为f_i下的信号片段，其离散表示为

(6)

式中：n=1, 2, …, N；U为扫频信号u(n)的幅值。
在实际工作过程中，一般将舵系统近似为线性系统。测量到的舵偏转信号形式为

(7)

式中：Y为舵偏转信号y(n)的绝对幅值; A为传递函数幅值; ф为传递函数相位; N_y(n)为测量信号中的噪声。
将式(7)改写为关于Acos ф和Asin ф的矩阵形式：

(8)

简写为

(9)

通过最小二乘法寻找x=[Acos ф ??Asin ф ]^T的最优解为

(10)

得到在频率f_i下的频响值为

(11)

对每个频率线都做以上的最小二乘运算，就能够得到系统整体的频响特性，并且最小二乘的计算过程本身具有很强的滤波功能，能减小测量噪声对计算结果的误差。
3 舵机动力学系统辨识方法 3.1 辨识问题描述 子空间辨识在实测频响函数曲线的基础上，辨识获得舵机的离散状态空间模型，再将离散域模型转换到连续域，获得连续状态空间模型，利用该模型重构得到舵机传递函数。
考虑n阶线性时不变连续状态空间模型，其表征为

(12)

式中：A_c、B_c、C_c和D_c分别为连续状态空间模型的系统矩阵、输入矩阵、输出矩阵和直接传递矩阵；x为系统的状态变量；u为系统的输入；y为系统的输出；下标c表示连续模型标识。离散后的状态空间模型为

(13)

式中：A、B、C、D分别为离散状态空间模型的系统矩阵、输入矩阵、输出矩阵和直接传递矩阵；y_r∈ R^p为系统的输出；u_r∈ R^m为系统的输入；下标r表示采样点，即满足x_r=x(rΔt)，Δt为采样周期。
状态空间模型(13)的脉冲响应g_k表示为

(14)

则状态空间模型(13)的频响函数为

(15)

即有

(16)

虽然，状态空间描述并不唯一，描述{ A, B, C, D }与之间可以通过转换矩阵T进行转换，即

(17)

但可以证明，式(16)与T无关。
于是，系统辨识的问题描述如下：
已知实测(包含噪声)M个样本点的频响函数：

(18)

式中：n_k为噪声。寻找有限维的频响函数使得

(19)

式中：

(20)

其中：σ₁定义为最大奇异值。当n_k= 0时，给定有限个样本点M，存在M₀ < ∞使得

(21)

一般情况下，连续状态空间模型与离散状态空间模型之间满足信号的零阶采样保持(ZOH)假设。但是，大部分航模舵机均采用50 Hz的PWM控制，对于这种低采样率的离散系统，需要采用双线性变换，减小由于假设引起的误差^[14]。
假设连续系统的传递函数为G_c(s)，离散系统的传递函数为G_d(z)，下标c表示连续模型标识，下标d表示离散模型标识。当获得了离散系统的传递函数时，即可通过逆变换获得连续系统的传递函数

(22)

连续和离散状态空间模型之间的双线性变换为^[14]

(23)

式中：ΔT为系统的采样周期，需要注意的是，2/ΔT不是连续系统的极点。
3.2 频域子空间辨识 本文采用的子空间辨识算法的主要步骤如下：
步骤1??扩展频响函数样本点

(24)

式中：(·)^*为(·)的共轭。
步骤2??定义为2M个样本点的离散傅里叶逆变换

(25)

步骤3??构建Hankel矩阵

(26)

式中：q>n，r≥n，q+r≤2M。
步骤4??对Hankel矩阵进行奇异值分解

(27)

式中：为矩阵的左奇异矩阵；为矩阵的右奇异矩阵；为矩阵的奇异值矩阵。
步骤5??确定系统阶数n，通过矩阵重排，使得包含前n个最大的奇异值，将式(27)按分块矩阵展开为

(28)

式中：为对应的奇异向量构成的左、右奇异矩阵；为剩余奇异值对角阵对应的奇异向量构成的左、右奇异矩阵。
步骤6??确定系统矩阵和

(29)

(30)

式中：[·]⁺表示[·]的Moore-Penrose逆。

(31)

步骤7??求解最小二乘问题确定和

(32)

式中：‖·‖_F为F范数。
步骤8??求解辨识得到的传递函数

(33)

步骤9??将离散域模型转换为连续域模型，获得系统传递函数。
4 舵机频响特性测试试验设计 针对文献[15]中选用的三款舵机分别进行了无气动载荷和有气动载荷的频响特性测试，并针对实际使用中的带气动载荷情况，进行了舵机传递函数参数辨识。
4.1 测试对象 1) 舵机A
JR SERVO DS8900G舵机(见图 9)为无芯数字高速伺服电机，其主要参数如下：伺服电压为4.8 V，转速为0.05 s/60°(4.8 V)，扭矩为0.336 N·m(4.8 V)，质量为59 g。该舵机采用50 Hz PWM控制。

图 9 JR SERVO DS8900G舵机 Fig. 9 JR SERVO DS8900G actuator

图选项

2) 舵机B
Hitec HS-7954SH舵机(见图 10)为空心杯大扭力数字舵机，其主要参数如下：伺服电压为6.0 V，转速为0.15 s/60°(6.0 V)，扭力为24 kg/cm(6.0 V)，质量为65 g。该舵机采用50 Hz PWM控制。

图 10 Hitec HS-7954SH舵机 Fig. 10 Hitec HS-7954SH actuator

图选项

3) 舵机C
MaxonMotor EC-4pole22-323218舵机(见图 11)为直流高功率无刷电机，与减速器、传感器和控制器组成舵机伺服系统。该舵机的主要参数如下：电机电阻为0.527 Ω，堵转转矩为0.639 N· m，额定转矩为0.045 1 N·m，最大允许转速为15 000 r/min，转矩常数为0.014 N·m/A，速度常数为680 r/(min·V)，机械时长为1.48 ms，转子惯量为5.54 g·cm²，减速器传动比i=50，质量为223 g。

图 11 MaxonMotor EC-4pole22-323218舵机 Fig. 11 MaxonMotor EC-4pole22-323218 actuator

图选项

4.2 测试布置 根据飞机模型的舵面参数^[15]确定舵机频响特性测试平台中配重的质量和位置，并且根据该舵面实际的飞行状态，计算确定测试平台中的扭杆尺寸，其结果如表 1所示。
表 1 舵机频响特性测试平台参数 Table 1 Parameters of actuator frequency response characteristic test platform

参数	数值
舵面转动惯量/(kg·mm)	2.40
舵面铰链刚度/(N·m·(°)-1)	0.002 2

表选项






利用舵机频响特性测试系统进行舵机的频响特性测试，通过给待测舵机输入正弦扫频信号，同时测量扫频信号与待测舵机的舵响应信号，绘制待测舵机的频响函数。舵机频响特性测试的输入与测量参数如表 2所示，在每个频率点的信号长度不少于10个周期。
表 2 舵机频响特性测试参数 Table 2 Parameters of actuator frequency response characteristic test

参数	航机A	航机B	航机C
采样率/Hz	50	50	1 000
扫频初始频率/Hz	0.2	0.1	0.1
扫频终止频率/Hz	8	10	15
步进步长/(Hz·次-1)	0.1	0.1	0.1

表选项






通过地面试验得到待测舵机的传递函数，利用子空间辨识获得舵机的传递函数。
5 结果与讨论 分别对无气动载荷情况、有气动载荷情况测试舵机A的频响函数。选用舵机时，通常认为舵机的幅频特性>-3 dB即满足使用条件，由此可确定舵机的正常使用带宽；另外，实际舵机受内部机械结构和舵机采样率的影响存在时滞环节，由相频特性反映的舵机时滞时长也是舵机选用时的重要参数。
5.1 舵机A
5.1.1 实测频响函数 无气动载荷情况的实测频响函数如图 12所示。图 12显示在0.1~10 Hz区间内，该舵机的幅频特性较好。可以看出，该舵机的正常使用带宽在8 Hz以上。但从相频特性可以看出，该舵机存在明显的时滞环节，不同扫频幅值下的时滞时长差异明显。当扫频幅值在2°时，该舵机的时滞时长为0.06 s；当扫频幅值为5°时，时滞时长变为0.24 s。当扫频幅值为2°与3°时，其相位特性一致。

图 12 无气动载荷时舵机A随扫频信号幅值变化的频响特性 Fig. 12 Frequency response characteristics of Actuator A with sweep signal amplitude without aerodynamic load

图选项

图 13显示了舵机A有气动载荷时的频响特性。与无气动载荷不同的是，当扫频幅值为3°时时滞最明显，时滞时长为0.22 s。当扫频幅值为5°时，时滞时长变为0.07 s。当扫频幅值为2°与4°时，其相位特性一致。

图 13 有气动载荷时舵机A随扫频信号幅值变化的频响特性 Fig. 13 Frequency response characteristics of Actuator A with sweep signal amplitude with aerodynamic loads

图选项

图 14对比了扫频幅值为5°时有/无气动载荷对舵机特性的影响。可以看出，气动载荷不仅降低了舵机A的幅频特性，同时也显著改变了该舵机的相频特性。

图 14 扫频幅值为5°时，有/无气动载荷下的舵机A频响特性 Fig. 14 Frequency response characteristics of Actuator A under 5° sweep signal amplitude with/without aerodynamic load

图选项

5.1.2 舵机传递函数辨识 在考虑舵机时滞的影响下，引入时滞环节(e^－τs)，并对实测的传递函数进行辨识。辨识得到状态空间模型之后，利用式(16)重构离散系统的传递函数，并通过双线性变换构建连续系统的传递函数。
该舵机辨识得到的频响函数曲线与实测曲线如图 15所示。可以看出，在0.1~10 Hz的频率区间内，其频响特性基本吻合。

图 15 扫频幅值为5°时，有气动载荷下的舵机A频响函数曲线对比 Fig. 15 Frequency response function curves identification of Actuator A under 5° sweep signal amplitude with aerodynamic loads

图选项

5.2 舵机B
5.2.1 实测频响函数 从实测频响函数曲线看出(见图 16和图 17)，在0.2~8 Hz区间内，该舵机无气动载荷时的幅频特性随着扫频频率的增加而下降。其整体的幅频特性随着扫频幅值的增加而降低，并且该舵机的正常使用带宽小于4 Hz。

图 16 无气动载荷时舵机B随扫频信号幅值变化的频响特性 Fig. 16 Frequency response characteristics of Actuator B with sweep signal amplitude without aerodynamic load

图选项

图 17 有气动载荷时舵机B随扫频信号幅值变化的频响特性 Fig. 17 Frequency response characteristics of Actuator B with sweep signal amplitude with aerodynamic loads

图选项

在该舵机中也存在明显的时滞，但时滞时长随扫频幅值的变化不大。从图 18看出，有/无气动载荷对该舵机的相频特性影响较小，施加气动载荷会小幅降低该舵机的幅频特性。

图 18 扫频幅值为5°时，有/无气动载荷下的舵机B频响特性 Fig. 18 Frequency response characteristics of Actuator B under 5° sweep signal amplitude with/without aerodynamic load

图选项

5.2.2 舵机传递函数辨识 考虑该舵机的时滞并对该舵机进行传递函数辨识，得到的频响函数曲线与实测曲线如图 19所示。可以看出，在0.2~8 Hz的频率区间内，其频响特性吻合。该舵机中存在有τ=0.1 s左右的时滞时长。

图 19 扫频幅值为5°时，有气动载荷下的舵机B频响函数曲线对比 Fig. 19 Frequency response function curves identification of Actuator B under 5° sweep signal amplitude with aerodynamic loads

图选项

5.3 舵机C
5.3.1 实测频响函数 如图 20和图 21所示，从实测频响曲线可以看出，该舵机的频响特性良好，不存在明显的时滞。在扫频幅值为7°时，该舵机正常使用带宽也达到了8 Hz。随着扫频幅值从2°到7°不断增大，舵机带宽也不断降低。

图 20 无气动载荷时舵机C随扫频信号幅值变化的频响特性 Fig. 20 Frequency response characteristics of Actuator C with sweep signal amplitude without aerodynamic load

图选项

图 21 有气动载荷时舵机C随扫频信号幅值变化的频响特性 Fig. 21 Frequency response characteristics of Actuator C with sweep signal amplitude with aerodynamic loads

图选项

在施加气动载荷之后，如图 22所示，该舵机的频响特性与无气动载荷时的特性一致。

图 22 扫频幅值为7°时，有/无气动载荷下的舵机C频响特性 Fig. 22 Frequency response characteristics of Actuator C under 7° sweep signal amplitude with/without aerodynamic load

图选项

5.3.2 舵机传递函数辨识 通过辨识发现，虽然该舵机存在有τ=14 ms的时滞时长，远大于舵机的机械时长1.48 ms，但时滞时长远小于正常使用带宽对应周期，对实际使用的影响较小。
对于该舵机，辨识得到的频响函数曲线与实测曲线如图 23所示，在0.1~15 Hz的频率区间内，辨识与实测结果一致。

图 23 扫频幅值为7°时，有气动载荷下的舵机C频响函数曲线对比 Fig. 23 Frequency response function curves identification of Actuator C under 7° sweep signal amplitude with aerodynamic loads

图选项

5.4 频响特性测试结果 通过测试结果对比，获得3款舵机的频响特性，舵机的主要频响特性参数如表 3所示。3款舵机均存在不同程度的时滞，一方面是由数字信号采样频率引起，另一方面是受舵机元件中的间隙和摩擦影响。
表 3 舵机频响特性测试结果 Table 3 Results of actuator frequency response characteristic test

舵机	扫频幅值/(°)	无气动载荷带宽/Hz	有气动载荷带宽/Hz	时滞时长/s
A	2	0.9	0.9	0.06~0.24
	3	>10	1.6
	4	>10	7.6
	5	>10	>10
B	2	3.6	3.6	0.1
	3	3.2	2.9
	4	3.0	3.6
	5	2.6	2.9
C	2	>15	>15	0.014
	3	12.5	12.5
	4	10.7	10.7
	5	9.6	9.6
	6	8.7	8.7
	7	8.1	8.1

表选项






6 结论 1) 针对航模舵机的频响特性测试设计了一种舵机频响特性测试平台，利用该平台可以模拟舵面惯性载荷，可以对大部分航模舵机进行有/无气动载荷的测试，提出了航模舵机频响特性测试方法，利用子空间辨识获得舵机的传递函数。
2) 在舵机的标称扭矩内，工作载荷的变化对普通小型舵机的幅频特性影响较大，载荷增加会使幅频曲线整体降低，扫频幅值的增加会降低航模舵机带宽。
3) 普通小型舵机大多采用50 Hz PWM信号控制，受采样频率影响，舵机天然存在一个0.02 s的时滞，但通过频响特性测试发现，舵机的真实时滞时长远超于采样时间步长，这可能由舵机元件中的间隙和摩擦造成的。在使用时，应考虑采样频率和时滞的影响。
4) 应用频域子空间辨识方法和双线性变换，辨识得到舵机的传递函数，通过与实测频响函数的对比，验证了本文辨识方法的有效性。

参考文献

[1]	杨超. 飞行器气动弹性原理[M]. 3版. 北京: 北京航空航天大学出版社, 2016: 203-208. YANG C. Princple of aircraft aeroelasticity[M]. 3rd ed. Beijing: Beihang University Press, 2016: 203-208. (in Chinese)
[2]	陈怀民, 庄皓玥, 马松辉, 等. 基于试验的电动舵机模型辨识[J]. 航空计算技术, 2014, 44(4): 6-9. CHEN H M, ZHUANG H Y, MA S H, et al. Model identification of electromechanical actuator based on tests[J]. Aeronautical Computing Technique, 2014, 44(4): 6-9. DOI:10.3969/j.issn.1671-654X.2014.04.002 (in Chinese)
[3]	马敏, 周盛春, 王伯波, 等. 关于飞行器的电动舵机伺服系统参数辨识[J]. 计算机仿真, 2016, 33(2): 88-92. MA M, ZHOU S C, WANG B B, et al. Parameter identification of electric actuator servo system on aircraft[J]. Computer Simulation, 2016, 33(2): 88-92. DOI:10.3969/j.issn.1006-9348.2016.02.019 (in Chinese)
[4]	朱萌, 曹国武, 张志伟, 等. 基于Levy法的气动舵机参数辨识[J]. 弹箭与制导学报, 2011, 31(6): 69-72. ZHU M, CAO G W, ZHANG Z W, et al. The system identification of pneumatic actuator based on Levy method[J]. Journal of Projectiles, Rockets, Missiles and Guidance, 2011, 31(6): 69-72. DOI:10.3969/j.issn.1673-9728.2011.06.020 (in Chinese)
[5]	章家保, 刘慧, 贾宏光, 等. 电动舵机伺服系统的模型辨识及其校正[J]. 光学精密工程, 2008, 16(10): 1971-1976. ZHANG J B, LIU H, JIA H G, et al. Model identification and corrector design for servo system of electromechanical actuator[J]. Optics and Precision Engineering, 2008, 16(10): 1971-1976. DOI:10.3321/j.issn:1004-924X.2008.10.030 (in Chinese)
[6]	高素军.基于实测数据的系统辨识及实验研究[D].哈尔滨: 哈尔滨工业大学, 2006: 25-38. GAO S J.Research and experiment of system identification based on real-gained data[D].Harbin: Harbin Institute of Technology, 2006: 25-38(in Chinese). http://cdmd.cnki.com.cn/Article/CDMD-10213-2006171627.htm
[7]	MACIEJOWSKI J M. Guaranteed stability with subspace methods[J]. System and Control Letters, 1995, 26(2): 153-156. DOI:10.1016/0167-6911(95)00010-7
[8]	OVERSCHEE P V, MOOR B D. Subspace identification for linear systems[M]. Dordrecht: Kluwer Academic Publishers, 1996: 8-11.
[9]	OVERSCHEE P V, MOOR B D. A unifying theorem for three subspace system identification algorithms[J]. Automatica, Special Issue on Trends in System Identification, 1995, 31(12): 1853-1864.
[10]	DE CALLAFON R, DE ROOVER D, VAN DEN HOF P.Multivariable least squares frequency domain identification using polynomial matrix fraction descriptions[C]//Proceedings of the 35th IEEE Conference on Decision and Control.Piscataway, NJ: IEEE Press, 1996, 2: 2030-2035. https://www.researchgate.net/publication/2432455_Multivariable_Least_Squares_Frequency_Domain_Identication_using_Polynomial_Matrix_Fraction_Descriptions
[11]	LJUNG L. System identification:Theory for the user[M]. 2nd ed. Upper Saddle River: Prentice Hall, 1999.
[12]	CLAES M, GRAHAM M R, CALLAFON R A. Frequency domain subspace identification of a tape servo system[J]. Microsystem Technologies, 2007, 13(8-10): 1439-1447. DOI:10.1007/s00542-007-0389-y
[13]	张仁嘉.飞行器气动伺服弹性若干关键问题研究[D].北京: 北京航空航天大学, 2015: 101-102. ZHANG R J.Extensional research on several critical aeroservoelastic problems of air vehicles[D].Beijing: Beihang University, 2015: 101-102(in Chinese).
[14]	MCKELVEY T, AKCAY H, LJUNG L. Subspace-based multivariable system identification from frequency response data[J]. IEEE Transactions on Automatic Control, 1996, 41(7): 960-979. DOI:10.1109/9.508900
[15]	黄超.柔性飞翼飞机颤振主动抑制系统建模、设计与验证[D].北京: 北京航空航天大学, 2018: 89-90. HUANG C.Modeling, design and verification of active flutter suppression system acting on flexible flying wing aircraft[D].Beijing: Beihang University, 2018: 89-90(in Chinese).