近Nyquist頻率情況下的數(shù)字化陷波器設(shè)計(jì)方法研究

趙艷輝，張擁軍，李海峰，趙 霞

(中國(guó)空空導(dǎo)彈研究院，河南 洛陽 471009)

0 引 言

受氣動(dòng)力和慣性力的作用，導(dǎo)彈在飛行過程中會(huì)產(chǎn)生彈性形變和模態(tài)振動(dòng)。彈體的彈性振動(dòng)信號(hào)通過彈上傳感器進(jìn)入控制系統(tǒng)，導(dǎo)致控制面的附加輸出，使彈性振動(dòng)加劇，嚴(yán)重時(shí)使控制回路失去穩(wěn)定性。因此,在設(shè)計(jì)飛行控制系統(tǒng)時(shí),需對(duì)彈性振動(dòng)進(jìn)行抑制，使控制系統(tǒng)的頻率特性在彈體模態(tài)頻率處獲得適當(dāng)?shù)脑鲆嫠p，保證控制系統(tǒng)的伺服彈性穩(wěn)定性。文獻(xiàn)[1-2]指出，高超聲速飛行器在氣動(dòng)/結(jié)構(gòu)/推進(jìn)系統(tǒng)間存在嚴(yán)重的耦合效應(yīng)，彈性振動(dòng)是控制系統(tǒng)必須解決的問題之一。文獻(xiàn)[3-5]對(duì)彈箭和飛機(jī)控制系統(tǒng)設(shè)計(jì)中彈性振動(dòng)抑制問題的解決方案進(jìn)行了全面的綜述,并指出合理設(shè)計(jì)陷波器是實(shí)現(xiàn)振動(dòng)模態(tài)抑制的基本方法之一。文獻(xiàn)[6] 利用遞歸最小二乘算法在線估計(jì)彈體一階彎曲頻率，并以估計(jì)值作為結(jié)構(gòu)濾波器的中心頻率，實(shí)現(xiàn)了對(duì)彈體彎曲模態(tài)振動(dòng)的抑制。文獻(xiàn)[7]以氣動(dòng)伺服彈性系統(tǒng)的彈性模態(tài)頻響峰值最小作為優(yōu)化目標(biāo)，以剛體模態(tài)頻響特性作為設(shè)計(jì)約束，建立了一種基于多目標(biāo)遺傳算法的結(jié)構(gòu)陷幅濾波器優(yōu)化設(shè)計(jì)方法。文獻(xiàn)[8]提出一種基于修正反正切函數(shù)的陷波器自適應(yīng)算法，克服了梯度自適應(yīng)算法收斂速度慢、對(duì)噪聲敏感等缺點(diǎn)。文獻(xiàn)[9]提出一種并行主動(dòng)噪聲濾波器自適應(yīng)算法，可以獲得多個(gè)頻率的精確估計(jì)值，有效實(shí)現(xiàn)了主動(dòng)噪聲控制。文獻(xiàn)[10] 通過粒子群多參數(shù)尋優(yōu)，以時(shí)域峭度最大原則對(duì)陷波器中心頻率及帶寬進(jìn)行自適應(yīng)選取，提高了陷波效率及準(zhǔn)確性。文獻(xiàn)[11] 提出一種基于小波變換與奇異值分解的彈性自適應(yīng)陷波濾波方法，在線估計(jì)彈性振型頻率，實(shí)現(xiàn)自適應(yīng)陷波。上述文獻(xiàn)提供了多種結(jié)構(gòu)濾波器的設(shè)計(jì)方法，部分方法在連續(xù)域設(shè)計(jì)結(jié)構(gòu)濾波器，并分析其性能；在線辨識(shí)彈性模態(tài)頻率的方法對(duì)振蕩信號(hào)的信噪比有特定要求，限制了其在工程實(shí)際中的應(yīng)用。

在工程實(shí)際中，結(jié)構(gòu)濾波器的離散化過程是實(shí)現(xiàn)彈性模態(tài)振動(dòng)抑制的關(guān)鍵環(huán)節(jié)，需要在嵌入飛控計(jì)算機(jī)的飛控軟件中實(shí)現(xiàn)與連續(xù)域結(jié)構(gòu)濾波器性能一致的數(shù)字化濾波器。目前，結(jié)構(gòu)濾波器的數(shù)字化基于控制系統(tǒng)采樣頻率遠(yuǎn)高于濾波中心頻率這一假設(shè)展開。在飛控計(jì)算機(jī)采樣頻率一定的情況下，這種假設(shè)有時(shí)不能成立，甚至出現(xiàn)陷波中心頻率接近Nyquist頻率(0.5倍采樣頻率)的情況。采樣定理和數(shù)字控制系統(tǒng)設(shè)計(jì)理論表明：隨著濾波中心頻率接近Nyquist頻率，離散化后的數(shù)字濾波器會(huì)發(fā)生頻率特性畸變。從公開發(fā)表的文獻(xiàn)看，中心陷波頻率接近Nyquist頻率情況下的結(jié)構(gòu)濾波器設(shè)計(jì)和離散化問題，目前仍是一項(xiàng)技術(shù)空白，如何解決這一問題是本文研究的重點(diǎn)。

1 常用結(jié)構(gòu)濾波器離散化方法及其缺點(diǎn)

通常，數(shù)字化結(jié)構(gòu)濾波器的設(shè)計(jì)步驟是：按照陷波中心頻率處的增益衰減要求，在連續(xù)域設(shè)計(jì)結(jié)構(gòu)濾波器；通過合適的離散化方法，將連續(xù)域的濾波器傳遞函數(shù)離散化，并應(yīng)用于數(shù)字化飛行控制系統(tǒng)。

復(fù)變量z和復(fù)變量s之間的轉(zhuǎn)換關(guān)系為

z=eTss

(1)

式中：Ts為控制系統(tǒng)的采樣頻率。

常用的離散化方法有雙線性變換法、預(yù)畸校正法和零極點(diǎn)匹配法[12]。

1.1 雙線性變換法(Tustin Method)

利用式(1)，在s=0的領(lǐng)域內(nèi)進(jìn)行1階Taylor展開，得到z的近似表達(dá)式：

(2)

利用式(2)反解連續(xù)域復(fù)變量s，得到雙線性變換公式：

(3)

雙線性變換方法簡(jiǎn)單，應(yīng)用廣泛，離散化后的濾波器一定是穩(wěn)定的。但由于連續(xù)域的頻率ωs和離散域的頻率ωd的非線性關(guān)系，造成ωd的扭曲現(xiàn)象：

在離散域看到的頻率ωd存在如下關(guān)系：

(4)

通過雙線性變換，將ωs映射到離散域后，得到的ωd小于ωs。隨著ωs接近于Nyquist頻率(0.5倍采樣頻率)，扭曲現(xiàn)象愈發(fā)嚴(yán)重。ωs和ωd之間的非線性關(guān)系如圖1所示。

圖1 ωs和ωd之間的非線性關(guān)系Fig.1 Nonlinear relation between ωs and ωd

這種頻率扭曲現(xiàn)象隨著頻率逐漸接近Nyquist頻率越發(fā)嚴(yán)重，甚至?xí)菇Y(jié)構(gòu)濾波器失效。

1.2 預(yù)畸校正法(Prewaping Method)

為改善雙線性變換離散化方法的缺陷，需要對(duì)結(jié)構(gòu)濾波器的陷波中心頻率進(jìn)行預(yù)畸校正。根據(jù)離散域的中心頻率的期望值，對(duì)連續(xù)域?yàn)V波器的中心頻率進(jìn)行預(yù)畸校正，使ωd與ωs一致，預(yù)畸校正過程如下：

簡(jiǎn)化得

(5)

通過預(yù)畸校正可以使離散化的濾波器和連續(xù)域?yàn)V波器的中心頻率和衰減深度一致。但隨著濾波中心頻率接近Nyquist頻率，結(jié)構(gòu)濾波器在中心頻率附近的陷波寬度逐漸減小。這是由于預(yù)畸校正方法僅能對(duì)某一個(gè)特定頻率進(jìn)行預(yù)畸校正所造成的。

1.3 零極點(diǎn)匹配法(Matched Pole-Zero Method)

利用式(1)可以將連續(xù)域?yàn)V波器的零極點(diǎn)映射為離散域?yàn)V波器的零極點(diǎn)。按照如下步驟進(jìn)行連續(xù)域?yàn)V波器的零極點(diǎn)匹配離散化過程：

(1) 將連續(xù)域傳遞函數(shù)Gc(s)的有限極點(diǎn)在z域中用z=eTss表示。

(2) 將連續(xù)域傳遞函數(shù)Gc(s)的有限零點(diǎn)在z域中用z=eTss表示。

(6)

(7)

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

(4) 調(diào)整Gd(z)的增益Kz。

對(duì)于結(jié)構(gòu)濾波器設(shè)計(jì)，應(yīng)滿足約束條件：

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

(8)

零極點(diǎn)匹配方法避免了預(yù)畸校正方法帶來的陷波寬度變窄的問題，但是隨著濾波中心頻率接近Nyquist頻率，濾波器的相位滯后和幅值衰減程度遠(yuǎn)超過連續(xù)域傳遞函數(shù)的設(shè)計(jì)結(jié)果，需要通過合理的設(shè)計(jì)方法加以克服。

2 數(shù)字化結(jié)構(gòu)濾波器設(shè)計(jì)

2.1 數(shù)字化結(jié)構(gòu)濾波器解析表達(dá)式

結(jié)構(gòu)濾波器連續(xù)域傳遞函數(shù)：

(9)

雙線性變換后得到的數(shù)字濾波器：

(10)

預(yù)畸校正后，再經(jīng)雙線性變換得到數(shù)字濾波器：

(11)

零極點(diǎn)匹配離散化后得到數(shù)字濾波器：

(12)

利用歐拉公式得到化簡(jiǎn)后的離散化傳遞函數(shù)：

(13)

考慮到控制系統(tǒng)的有限帶寬特性，利用低頻匹配關(guān)系式(8)確定結(jié)構(gòu)濾波器歸一化增益Kz：

(14)

式(10)～(11)和式(13)為不同離散化方法對(duì)應(yīng)的數(shù)字化結(jié)構(gòu)濾波器解析表達(dá)式，可應(yīng)用于在線實(shí)時(shí)更新濾波器參數(shù)的情形，也可以作為眾多結(jié)構(gòu)濾波器設(shè)計(jì)方法離散化過程的有益補(bǔ)充。

2.2 近Nyquist頻率情況下的濾波性能改進(jìn)方法

當(dāng)結(jié)構(gòu)濾波器中心頻率接近Nyquist頻率時(shí)，現(xiàn)有離散化方法存在缺點(diǎn)，本文給出2種通過修正連續(xù)域結(jié)構(gòu)濾波器分母多項(xiàng)式阻尼比ξ2至一個(gè)期望值ξ2c來改進(jìn)零極點(diǎn)匹配離散化設(shè)計(jì)結(jié)果的優(yōu)化方法。

2.2.1 方法1：衰減增益約束條件下的修正方法

假設(shè)修正后的濾波器離散化結(jié)果如下：

(15)

若Fnf(z)在指定的關(guān)鍵頻率ωc處的增益和式(9)所示的連續(xù)域?yàn)V波器的增益相同，即滿足約束關(guān)系式：

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

(16)

式中：ωc=2π(fm-Δf)，fm為濾波器中心頻率，Δf為彈性模態(tài)頻率的最大攝動(dòng)量。

式(16)為超越方程，難以獲得解析解，可以通過數(shù)值解法獲得ξ2c的數(shù)值解。式(16)右端的模為常值，左端的模為隨ξ2c變化的函數(shù)，二者交點(diǎn)對(duì)應(yīng)的橫坐標(biāo)即為ξ2c的數(shù)值解。

2.2.2 方法2：相位滯后約束條件下的修正方法

對(duì)于式(15)，若Fnf(z)在指定的角頻率ωl處的相位和式(9)所示的連續(xù)域?yàn)V波器的相位相同，即滿足約束關(guān)系式：

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

(17)

式中：ωl=2πfl，fl為控制系統(tǒng)設(shè)計(jì)者所關(guān)心的與控制系統(tǒng)回路帶寬有關(guān)的低頻段頻率值。

式(17)為超越方程，難以獲得解析解，可以通過數(shù)值解法獲得ξ2c的數(shù)值解。式(17)右端的相位為常值，左端的相位為隨ξ2c變化的函數(shù)，二者交點(diǎn)對(duì)應(yīng)的橫坐標(biāo)即為ξ2c的數(shù)值解。

3 設(shè)計(jì)實(shí)例與仿真驗(yàn)證

通過設(shè)計(jì)實(shí)例演示現(xiàn)有離散化方法的缺點(diǎn)，并證明近Nyquist頻率情況下濾波性能改進(jìn)方法的有效性。設(shè)控制系統(tǒng)的離散化周期Ts為0.004 s，分別在陷波中心頻率遠(yuǎn)小于及接近Nyquist頻率的情況下進(jìn)行仿真驗(yàn)證。

陷波器中心頻率40 Hz及110 Hz情況下3種離散化方法的頻率特性比較結(jié)果如圖2所示。

圖2 陷波器中心頻率40 Hz,110 Hz情況下不同離散化方法的頻率特性比較結(jié)果Fig.2 Frequency characteristic comparison of different discretization methods when center frequency is 40 Hz and 110 Hz

從圖2 (a) 可知,濾波器中心頻率遠(yuǎn)小于Nyquist頻率的情況下，雙線性變換形成的離散化濾波器的中心頻率為37 Hz,偏離了設(shè)計(jì)值40 Hz；通過預(yù)畸校正可以使離散化后的濾波器中心頻率與連續(xù)域設(shè)計(jì)結(jié)果一致；零極點(diǎn)匹配離散化結(jié)果的幅頻特性更接近連續(xù)域的設(shè)計(jì)結(jié)果。

從圖2 (b) 可知,濾波器中心頻率接近Nyquist頻率的情況下，通過預(yù)畸校正仍可以使離散化后的濾波器中心頻率與連續(xù)域設(shè)計(jì)結(jié)果一致，但濾波器的陷波寬度變窄；彈體彈性模態(tài)頻率的攝動(dòng)超過一定邊界的情況下，經(jīng)預(yù)畸校正獲得的離散化濾波器將失去彈性濾波作用；采用零極點(diǎn)匹配方法獲得的數(shù)字化結(jié)構(gòu)濾波器最為有效，但需要克服低頻相位滯后偏大的缺點(diǎn)。

通過方法1得到修正阻尼比ξ2c的數(shù)值解如圖3所示，數(shù)字化結(jié)構(gòu)濾波器的性能如圖4所示。

圖3 采用方法1獲得的ξ2c的數(shù)值解Fig.3 Numerical solution of ξ2c obtained by method 1

圖4 采用方法1獲得的數(shù)字化結(jié)構(gòu)濾波器設(shè)計(jì)結(jié)果Fig.4 Design results of digital notch filter obtained by method 1

從圖3～4可以看出，方法1通過將分母多項(xiàng)式的阻尼比ξ2c的數(shù)值修正為0.4，可以使得近Nyquist頻率情況下的離散域?yàn)V波器的陷波特性更接近連續(xù)域設(shè)計(jì)結(jié)果，同時(shí)相位滯后也相應(yīng)減小。

通過方法2得到修正阻尼比ξ2c的數(shù)值解如圖5所示，數(shù)字化結(jié)構(gòu)濾波器的性能如圖6所示。

圖6 采用方法2獲得的數(shù)字化結(jié)構(gòu)濾波器設(shè)計(jì)結(jié)果Fig.6 Design results of digital notch filter obtained by method 2

從圖5～6可以看出,方法2通過將分母多項(xiàng)式的阻尼比ξ2c的數(shù)值修正為0.35，可以使得近Nyquist頻率情況下的離散域?yàn)V波器低頻相位特性幾乎與連續(xù)域設(shè)計(jì)結(jié)果重合，陷波特性也更接近連續(xù)域設(shè)計(jì)結(jié)果。

圖5 采用方法2獲得的ξ2c的數(shù)值解Fig.5 Numerical solution of ξ2c obtained by method 2

假設(shè)濾波器輸入信號(hào)為幅值1°、頻率110 Hz的正弦信號(hào)，對(duì)上述濾波器的陷波效果進(jìn)行檢驗(yàn)。彈體彈性模態(tài)頻率為名義值時(shí)的仿真結(jié)果如圖7所示，彈體彈性模態(tài)頻率攝動(dòng)情況下的仿真結(jié)果如圖8～9所示。

圖7 彈性模態(tài)頻率取名義值110Hz時(shí)的仿真結(jié)果Fig.7 Simulationresultswhenelasticmodalfrequencyis110Hz圖8 彈性模態(tài)頻率取邊界值105Hz時(shí)的仿真結(jié)果Fig.8 Simulationresultswhenelasticmodalfrequencyis105Hz圖9 彈性模態(tài)頻率取邊界值115Hz時(shí)的仿真結(jié)果Fig.9 Simulationresultswhenelasticmodalfrequencyis115Hz

從圖8～9可以看出,在彈體彈性模態(tài)頻率攝動(dòng)的情況下，采用預(yù)畸校正離散化方法形成的數(shù)字化結(jié)構(gòu)濾波器已經(jīng)不能滿足彈性濾波需求；而本文給出的2種修正零極點(diǎn)匹配方法的濾波性能依舊能夠和連續(xù)域設(shè)計(jì)結(jié)果保持一致，證明了本文提出的近Nyquist頻率情況下修正零極點(diǎn)匹配離散方法的有效性。

4 結(jié) 論

本文回顧了常用的結(jié)構(gòu)濾波器離散化方法，給出了3種數(shù)字化結(jié)構(gòu)濾波器的解析表達(dá)式，并分析了陷波中心頻率接近Nyquist頻率的情況下存在的缺點(diǎn)。通過分析離散域?yàn)V波器的頻率特性，得出了陷波中心頻率接近Nyquist頻率的情況下適合采用零極點(diǎn)匹配方法離散化方法的結(jié)論，并給出了2種修正分母多項(xiàng)式阻尼比的數(shù)值方法，使修正零極點(diǎn)匹配方法獲得的數(shù)字化結(jié)構(gòu)濾波器性能和連續(xù)域設(shè)計(jì)結(jié)果一致；最后，通過輸入信號(hào)頻率存在攝動(dòng)的濾波仿真試驗(yàn)驗(yàn)證了本文提出的近Nyquist頻率情況下修正零極點(diǎn)匹配離散化方法的有效性。