近Nyquist频率情况下的数字化陷波器设计方法研究

赵艳辉，张拥军，李海峰，赵 霞

(中国空空导弹研究院，河南 洛阳 471009)

0 引 言

受气动力和惯性力的作用，导弹在飞行过程中会产生弹性形变和模态振动。弹体的弹性振动信号通过弹上传感器进入控制系统，导致控制面的附加输出，使弹性振动加剧，严重时使控制回路失去稳定性。因此,在设计飞行控制系统时,需对弹性振动进行抑制，使控制系统的频率特性在弹体模态频率处获得适当的增益衰减，保证控制系统的伺服弹性稳定性。文献[1-2]指出，高超声速飞行器在气动/结构/推进系统间存在严重的耦合效应，弹性振动是控制系统必须解决的问题之一。文献[3-5]对弹箭和飞机控制系统设计中弹性振动抑制问题的解决方案进行了全面的综述,并指出合理设计陷波器是实现振动模态抑制的基本方法之一。文献[6] 利用递归最小二乘算法在线估计弹体一阶弯曲频率，并以估计值作为结构滤波器的中心频率，实现了对弹体弯曲模态振动的抑制。文献[7]以气动伺服弹性系统的弹性模态频响峰值最小作为优化目标，以刚体模态频响特性作为设计约束，建立了一种基于多目标遗传算法的结构陷幅滤波器优化设计方法。文献[8]提出一种基于修正反正切函数的陷波器自适应算法，克服了梯度自适应算法收敛速度慢、对噪声敏感等缺点。文献[9]提出一种并行主动噪声滤波器自适应算法，可以获得多个频率的精确估计值，有效实现了主动噪声控制。文献[10] 通过粒子群多参数寻优，以时域峭度最大原则对陷波器中心频率及带宽进行自适应选取，提高了陷波效率及准确性。文献[11] 提出一种基于小波变换与奇异值分解的弹性自适应陷波滤波方法，在线估计弹性振型频率，实现自适应陷波。上述文献提供了多种结构滤波器的设计方法，部分方法在连续域设计结构滤波器，并分析其性能；在线辨识弹性模态频率的方法对振荡信号的信噪比有特定要求，限制了其在工程实际中的应用。

在工程实际中，结构滤波器的离散化过程是实现弹性模态振动抑制的关键环节，需要在嵌入飞控计算机的飞控软件中实现与连续域结构滤波器性能一致的数字化滤波器。目前，结构滤波器的数字化基于控制系统采样频率远高于滤波中心频率这一假设展开。在飞控计算机采样频率一定的情况下，这种假设有时不能成立，甚至出现陷波中心频率接近Nyquist频率(0.5倍采样频率)的情况。采样定理和数字控制系统设计理论表明：随着滤波中心频率接近Nyquist频率，离散化后的数字滤波器会发生频率特性畸变。从公开发表的文献看，中心陷波频率接近Nyquist频率情况下的结构滤波器设计和离散化问题，目前仍是一项技术空白，如何解决这一问题是本文研究的重点。

1 常用结构滤波器离散化方法及其缺点

通常，数字化结构滤波器的设计步骤是：按照陷波中心频率处的增益衰减要求，在连续域设计结构滤波器；通过合适的离散化方法，将连续域的滤波器传递函数离散化，并应用于数字化飞行控制系统。

复变量z和复变量s之间的转换关系为

z=eTss

(1)

式中：Ts为控制系统的采样频率。

常用的离散化方法有双线性变换法、预畸校正法和零极点匹配法[12]。

1.1 双线性变换法(Tustin Method)

利用式(1)，在s=0的领域内进行1阶Taylor展开，得到z的近似表达式：

(2)

利用式(2)反解连续域复变量s，得到双线性变换公式：

(3)

双线性变换方法简单，应用广泛，离散化后的滤波器一定是稳定的。但由于连续域的频率ωs和离散域的频率ωd的非线性关系，造成ωd的扭曲现象：

在离散域看到的频率ωd存在如下关系：

(4)

通过双线性变换，将ωs映射到离散域后，得到的ωd小于ωs。随着ωs接近于Nyquist频率(0.5倍采样频率)，扭曲现象愈发严重。ωs和ωd之间的非线性关系如图1所示。

图1 ωs和ωd之间的非线性关系Fig.1 Nonlinear relation between ωs and ωd

这种频率扭曲现象随着频率逐渐接近Nyquist频率越发严重，甚至会使结构滤波器失效。

1.2 预畸校正法(Prewaping Method)

为改善双线性变换离散化方法的缺陷，需要对结构滤波器的陷波中心频率进行预畸校正。根据离散域的中心频率的期望值，对连续域滤波器的中心频率进行预畸校正，使ωd与ωs一致，预畸校正过程如下：

简化得

(5)

通过预畸校正可以使离散化的滤波器和连续域滤波器的中心频率和衰减深度一致。但随着滤波中心频率接近Nyquist频率，结构滤波器在中心频率附近的陷波宽度逐渐减小。这是由于预畸校正方法仅能对某一个特定频率进行预畸校正所造成的。

1.3 零极点匹配法(Matched Pole-Zero Method)

利用式(1)可以将连续域滤波器的零极点映射为离散域滤波器的零极点。按照如下步骤进行连续域滤波器的零极点匹配离散化过程：

(1) 将连续域传递函数Gc(s)的有限极点在z域中用z=eTss表示。

(2) 将连续域传递函数Gc(s)的有限零点在z域中用z=eTss表示。

(6)

(7)

式中：m≤n；zi=-e-μiTs，i=1,…,m；pk=-e-λkTs，k=1,…,n。

(4) 调整Gd(z)的增益Kz。

对于结构滤波器设计，应满足约束条件：

Gd(z)|z=1=Gc(s)|s=0

(8)

零极点匹配方法避免了预畸校正方法带来的陷波宽度变窄的问题，但是随着滤波中心频率接近Nyquist频率，滤波器的相位滞后和幅值衰减程度远超过连续域传递函数的设计结果，需要通过合理的设计方法加以克服。

2 数字化结构滤波器设计

2.1 数字化结构滤波器解析表达式

结构滤波器连续域传递函数：

(9)

双线性变换后得到的数字滤波器：

(10)

预畸校正后，再经双线性变换得到数字滤波器：

(11)

零极点匹配离散化后得到数字滤波器：

(12)

利用欧拉公式得到化简后的离散化传递函数：

(13)

考虑到控制系统的有限带宽特性，利用低频匹配关系式(8)确定结构滤波器归一化增益Kz：

(14)

式(10)～(11)和式(13)为不同离散化方法对应的数字化结构滤波器解析表达式，可应用于在线实时更新滤波器参数的情形，也可以作为众多结构滤波器设计方法离散化过程的有益补充。

2.2 近Nyquist频率情况下的滤波性能改进方法

当结构滤波器中心频率接近Nyquist频率时，现有离散化方法存在缺点，本文给出2种通过修正连续域结构滤波器分母多项式阻尼比ξ2至一个期望值ξ2c来改进零极点匹配离散化设计结果的优化方法。

2.2.1 方法1：衰减增益约束条件下的修正方法

假设修正后的滤波器离散化结果如下：

(15)

若Fnf(z)在指定的关键频率ωc处的增益和式(9)所示的连续域滤波器的增益相同，即满足约束关系式：

|Fnf(z)|z=ejωcTs|=|Fnf(s)|s=jωc|

(16)

式中：ωc=2π(fm-Δf)，fm为滤波器中心频率，Δf为弹性模态频率的最大摄动量。

式(16)为超越方程，难以获得解析解，可以通过数值解法获得ξ2c的数值解。式(16)右端的模为常值，左端的模为随ξ2c变化的函数，二者交点对应的横坐标即为ξ2c的数值解。

2.2.2 方法2：相位滞后约束条件下的修正方法

对于式(15)，若Fnf(z)在指定的角频率ωl处的相位和式(9)所示的连续域滤波器的相位相同，即满足约束关系式：

∠Fnf(z)|z=ejωlTs=∠Fnf(s)|s=jωl

(17)

式中：ωl=2πfl，fl为控制系统设计者所关心的与控制系统回路带宽有关的低频段频率值。

式(17)为超越方程，难以获得解析解，可以通过数值解法获得ξ2c的数值解。式(17)右端的相位为常值，左端的相位为随ξ2c变化的函数，二者交点对应的横坐标即为ξ2c的数值解。

3 设计实例与仿真验证

通过设计实例演示现有离散化方法的缺点，并证明近Nyquist频率情况下滤波性能改进方法的有效性。设控制系统的离散化周期Ts为0.004 s，分别在陷波中心频率远小于及接近Nyquist频率的情况下进行仿真验证。

陷波器中心频率40 Hz及110 Hz情况下3种离散化方法的频率特性比较结果如图2所示。

图2 陷波器中心频率40 Hz,110 Hz情况下不同离散化方法的频率特性比较结果Fig.2 Frequency characteristic comparison of different discretization methods when center frequency is 40 Hz and 110 Hz

从图2 (a) 可知,滤波器中心频率远小于Nyquist频率的情况下，双线性变换形成的离散化滤波器的中心频率为37 Hz,偏离了设计值40 Hz；通过预畸校正可以使离散化后的滤波器中心频率与连续域设计结果一致；零极点匹配离散化结果的幅频特性更接近连续域的设计结果。

从图2 (b) 可知,滤波器中心频率接近Nyquist频率的情况下，通过预畸校正仍可以使离散化后的滤波器中心频率与连续域设计结果一致，但滤波器的陷波宽度变窄；弹体弹性模态频率的摄动超过一定边界的情况下，经预畸校正获得的离散化滤波器将失去弹性滤波作用；采用零极点匹配方法获得的数字化结构滤波器最为有效，但需要克服低频相位滞后偏大的缺点。

通过方法1得到修正阻尼比ξ2c的数值解如图3所示，数字化结构滤波器的性能如图4所示。

图3 采用方法1获得的ξ2c的数值解Fig.3 Numerical solution of ξ2c obtained by method 1

图4 采用方法1获得的数字化结构滤波器设计结果Fig.4 Design results of digital notch filter obtained by method 1

从图3～4可以看出，方法1通过将分母多项式的阻尼比ξ2c的数值修正为0.4，可以使得近Nyquist频率情况下的离散域滤波器的陷波特性更接近连续域设计结果，同时相位滞后也相应减小。

通过方法2得到修正阻尼比ξ2c的数值解如图5所示，数字化结构滤波器的性能如图6所示。

图6 采用方法2获得的数字化结构滤波器设计结果Fig.6 Design results of digital notch filter obtained by method 2

从图5～6可以看出,方法2通过将分母多项式的阻尼比ξ2c的数值修正为0.35，可以使得近Nyquist频率情况下的离散域滤波器低频相位特性几乎与连续域设计结果重合，陷波特性也更接近连续域设计结果。

图5 采用方法2获得的ξ2c的数值解Fig.5 Numerical solution of ξ2c obtained by method 2

假设滤波器输入信号为幅值1°、频率110 Hz的正弦信号，对上述滤波器的陷波效果进行检验。弹体弹性模态频率为名义值时的仿真结果如图7所示，弹体弹性模态频率摄动情况下的仿真结果如图8～9所示。

图7 弹性模态频率取名义值110Hz时的仿真结果Fig.7 Simulationresultswhenelasticmodalfrequencyis110Hz图8 弹性模态频率取边界值105Hz时的仿真结果Fig.8 Simulationresultswhenelasticmodalfrequencyis105Hz图9 弹性模态频率取边界值115Hz时的仿真结果Fig.9 Simulationresultswhenelasticmodalfrequencyis115Hz

从图8～9可以看出,在弹体弹性模态频率摄动的情况下，采用预畸校正离散化方法形成的数字化结构滤波器已经不能满足弹性滤波需求；而本文给出的2种修正零极点匹配方法的滤波性能依旧能够和连续域设计结果保持一致，证明了本文提出的近Nyquist频率情况下修正零极点匹配离散方法的有效性。

4 结 论

本文回顾了常用的结构滤波器离散化方法，给出了3种数字化结构滤波器的解析表达式，并分析了陷波中心频率接近Nyquist频率的情况下存在的缺点。通过分析离散域滤波器的频率特性，得出了陷波中心频率接近Nyquist频率的情况下适合采用零极点匹配方法离散化方法的结论，并给出了2种修正分母多项式阻尼比的数值方法，使修正零极点匹配方法获得的数字化结构滤波器性能和连续域设计结果一致；最后，通过输入信号频率存在摄动的滤波仿真试验验证了本文提出的近Nyquist频率情况下修正零极点匹配离散化方法的有效性。