導彈系統(tǒng)振動頻率自適應估計方法與仿真

張 衡， 芮筱亭， 楊富鋒， 陳剛利

(南京理工大學 發(fā)射動力學研究所，南京 210094)

導彈系統(tǒng)振動頻率自適應估計方法與仿真

張 衡， 芮筱亭， 楊富鋒， 陳剛利

(南京理工大學 發(fā)射動力學研究所，南京 210094)

為了提高處于導彈高過載、強沖擊、大振動復合動態(tài)環(huán)境中捷聯(lián)慣組的輸出精度，設(shè)計了級聯(lián)型二階自適應陷波濾波器來估計導彈的彈性振動頻率。該濾波器采用帶有遺忘因子的遞推最小二乘法實時優(yōu)化每個陷波濾波器的參數(shù),具有計算復雜度低，收斂速度快等特點。仿真和實驗結(jié)果表明，此方法對于處在復合動態(tài)環(huán)境中導彈振動信號的頻率估計效果較好，且具有較好的抗噪聲干擾能力，適合彈上計算機進行實時計算。

自適應陷波濾波器；頻率估計；遞推最小二乘法

導彈發(fā)射與飛行過程中，隨著燃料不斷消耗，質(zhì)量和剛度分布將不斷變化，引起導彈彈性振動頻率隨飛行時間的函數(shù)不斷變化。捷聯(lián)式慣組的陀螺儀和加速度計直接固聯(lián)在導彈上，振動將傳遞到慣組器件上，引起陀螺儀質(zhì)心偏移和結(jié)構(gòu)變形，使慣性器件產(chǎn)生動態(tài)測量誤差。當慣組工作頻率接近其共振頻率時，將造成慣組動態(tài)精度大幅下降或者永久性破壞甚至輸出完全錯誤的信息，從而導致導彈完全失控。因此，設(shè)計一個適用于彈上非平穩(wěn)信號的頻率的實時估計的算法非常重要。

頻譜估計是數(shù)字信號處理的主要內(nèi)容之一，傳統(tǒng)的方法是基于快速傅里葉變換(FFT)，但它的估計性能不高，因為FFT的頻率分辨率僅與數(shù)據(jù)序列的長度有關(guān)，要提高FFT算法的分辨率只能靠增加采樣時間來實現(xiàn)，這在許多場合下尤其對于導彈的飛行過程來說是不允許的。現(xiàn)代譜估計算法中的參量法具有估計精度高和分辨率高的特性，但參量法譜估計的高分辨率和高精度特性只能在高信噪比條件下才能實現(xiàn)[1]。導彈飛行的主動段、自由飛行段以及再入段的某些瞬間的沖擊信號或干擾可能會對慣組產(chǎn)生很大的影響，因此如何分析和捕獲這些瞬態(tài)分量，對其進行有效的分析對于彈用慣組的自適應減振系統(tǒng)的設(shè)計、制造和維護具有重要意義。

在許多應用中，正弦信號可能受到非線性的影響，其中可能會產(chǎn)生諧波頻率分量。在這樣的環(huán)境下，我們想要估計信號的基頻以及任何諧波頻率，采用單個二階陷波濾波器去估計基頻與諧波的頻率是不夠的，因為它只能容納一個頻率分量。另一方面，應用高階無限脈沖響應(IIR)濾波器可能效率不高，由于它采用多個自適應濾波器系數(shù)[2]。因此，本文采用將二階自適應陷波濾波器級聯(lián)的方法進行導彈頻率的實時估計，該濾波器采用遞推最小二乘法實時優(yōu)化每個陷波濾波器的參數(shù),具有計算復雜度低，收斂速度快等特點。結(jié)果表明，該方法對于處在復合動態(tài)環(huán)境中導彈振動信號的頻率估計效果較好，且具有較好的抗噪聲干擾能力。

1 導彈彈性振動模型

導彈主動段飛行是在推力T的作用下做加速運動，此時，導彈受到均勻的分布慣性力。導彈具有較大的長徑比，且各個部分的質(zhì)量分布及剛度分布不同，因此，將導彈看成由n段Euler-Bernoulli梁組成，每段梁都具有均勻的線密度和彎曲剛度。由于受到慣性力的作用，每段梁所受到的軸向力沿軸向呈不均勻分布。采用多體系統(tǒng)傳遞矩陣法建立導彈動力學模型，如圖1所示。梁微段受力分析，如圖2所示。

圖1 導彈動力學模型Fig.1 The missile dynamic model

圖2 梁微段受力分析Fig.2 The force analysis of the sub-beam

對圖2所示的梁微段進行受力分析，梁微段dx沿

橫向所受的外力有：剪力fs(x)和-fs(x+dx)，軸向力N(x)和-N(x+dx)在軸上的投影。

根據(jù)牛頓第二定律，忽略二階小量，梁微段橫向運動滿足

(1)

(2)

此方程為一個非線性偏微分方程，這里將梁分為多段，每一段的軸向力視為常力處理，則梁的彎曲自由振動方程為

(3)

這是一個常系數(shù)的齊次線性偏微分方程。令y(x,t)=Y(x)eiωt，代入上式得

(4)

Y(x)=Acoshαx+Bsinhαx+
Ccosβx+Dsinβx

(5)

對Euler-Bernoulli梁，將Θz,Mz,Qy的表達式寫成矩陣形式可以得到

Z(x)=B(x)a

(6)

式中，Z(x)=[YΘzMzQy]Tx為狀態(tài)矢量，a=[A,B,C,D]T

(7)

(8)

所以受軸向力Euler-Bernoulli梁的傳遞矩陣為

(9)

由各段梁的傳遞方程可以得到系統(tǒng)的總傳遞方程為

Zn,n+1=UnUn-1…U2U1Z1,0=UallZ1,0

(10)

系統(tǒng)總傳遞矩陣為

(11)

系統(tǒng)的邊界條件為

(12)

代入式可以得到

(13)

刪去與零元素有關(guān)的行和列可以得到

(14)

式中

上式存在非零解，則

(15)

由上式可求得α和β，也即得到γ1和γ2，由γ2得到對應的ωp(p=1,2,…,n)及其對應的增廣特征矢量Vp(p=1,2,…,n)。從而可求解出廣義坐標qp(t)(p=1,2,…,n)，以及導彈的動力響應ν。

導彈的廣義坐標為

(16)

導彈的動力響應為

(17)

2 導彈振動頻率估計方法

在傳統(tǒng)的信號處理中，分析處理信號最常用的方法是Fourier變換，但是它針對的是周期性平穩(wěn)信號，依賴信號的全局信息，并不能反映信號的局部特征，故對分析非平穩(wěn)信號不具有效性[4]。對于導彈的發(fā)射與飛行過程，振動信號是非平穩(wěn)信號，針對其特點，本文采用將二階自適應陷波濾波器級聯(lián)的方法來處理。

圖3給出了自適應濾波器的一般結(jié)構(gòu)，其中，x(k)表示輸入信號。y(k)為輸出信號，d(k)表示期望信號，誤差信號e(k)=d(k)-y(k)。

由數(shù)字濾波器設(shè)計的一般原則,陷波器的設(shè)計有兩個基本要求:

一是其傳遞函數(shù)的零點應在單位圓上，以使陷波的陷阱深度為無窮大。滿足這一要求的多項式的系數(shù)應具有鏡像對稱的形式，如下式所示A(z-1)=1+a1z-1+a2z-1+…+anz-1+…+a1z-1+z-2n

(18)

式中，z是復數(shù)變量，n為陷波個數(shù)，當考慮單個陷波時n=1，上式變?yōu)?/p>

A(z-1)=1+a1z-1+z-2

(19)

第二個基本要求是其傳遞函數(shù)的極零點必須匹配，這樣除了陷波的頻率之外，其余的頻率完全不受陷波的影響。綜合這兩個基本要求，若在單位圓上頻率ω0處設(shè)置一對半徑r=1的共扼零點，同時在單位圓內(nèi)同樣頻率ω0處設(shè)置一對半徑為ρ(0≤ρ<1)的共扼極點，將會產(chǎn)生一個幅頻特性在ω0處強烈衰減，而在其他頻率成分處信號基本不受影響的陷波器。

如圖4為級聯(lián)型的二階陷波濾波器的計算流程圖。

圖4 級聯(lián)型二階陷波濾波器計算流程圖Fig.4 The calculation flow chart of cascade of second-order notch filters

Z域上二階陷波濾波器形式[5]為

(20)

式(20)中的濾波器結(jié)構(gòu)難以進行自適應陷波濾波計算，為此選用基于Regalia[6]提出的自適應陷波濾波器，形式為

(21)

式(20)與式(21)形式保持一致可得到β0(1+β1)=ρa(1+b)，可化為

β0=ρa(1+b)/(1+ρ2b)

為得到自適應陷波濾波器，參數(shù)的數(shù)目應保持最小。如果我們假設(shè)ρ等于1，

β1=ρ2b?ρb
β0(1+β1)=ρa(1+b)?a(1+ρb)

(22)

因此，β0=a成立，Regalia提出的方程可改寫為

(23)

一般地，當零點位于單位圓上時，能保證陷波濾波器有更好性能。因此，在式(23)中，參數(shù)b等于1，它可以簡化為一個下列形式的關(guān)于ρ和a的函數(shù)：

(24)

其中，ρ位于0.9和0.99之間，如果它接近1，濾波器的帶寬變窄。

陷波濾波器的輸入信號，記為U(n)。然后每個部分的輸出信號可以表示為

(25)

對于級聯(lián)型陷波濾波器，式(25)被下式取代

yk(n)=xk(n)+xk(n-2)+2akxk(n-1)

(26)

式中，xk(n)表示通過濾波器的分母的信號。為了設(shè)計自適應算法，提出優(yōu)化的損失函數(shù)如下

(27)

如式(28)～(30)顯示的遞推最小二乘[7]法是用來估計未知參數(shù)向量θ=[a1…ap]T。時間平均自相關(guān)φk(n)和時間平均互相關(guān)αk(n)由下式給出[8]

φk(n)=λφk(n-1)+A2(n)

(28)

αk(n)=λαk(n-1)+A(n)·B(n)

(29)

其中A(n)=2xk(n-1),B(n)=xk(n)+xk(n-2)

k階段的估計參數(shù)由下式給出

ak(n)=-φk(n)-1αk(n)

(30)

由于陷波極點位于單位圓內(nèi)，所以參數(shù)ak也在范圍ak∈[-1,1]。因此，導彈的振動頻率可以使用采樣時間ΔT估計為如下[9]

(31)

3 仿真結(jié)果及實驗驗證

3.1 仿真結(jié)果

為驗證算法的頻率估計效果，特以某型導彈為例，對其主動段飛行過程進行仿真，用自適應二階陷波濾波器級聯(lián)的方法對多體系統(tǒng)傳遞矩陣法計算得到的響應數(shù)據(jù)進行頻率的估計。首先，用多體系統(tǒng)傳遞矩陣法計算得到導彈的振動頻率ωp(p=1,2,…,n)及其對應的增廣特征矢量Vp(p=1,2,…,n)，從而可求解出廣義坐標qp(t)(p=1,2,…,n)以及導彈的動力響應ν，其中，導彈的動力響應的計算如式(17)所示。然后，用本文提出的自適應二階陷波濾波器級聯(lián)的方法對多體系統(tǒng)傳遞矩陣法計算得到的導彈的動力響應進行頻率的估計。由于一階模態(tài)對系統(tǒng)動力響應的影響最大，本文僅以識別導彈一階頻率為例，驗證此算法的有效性。

通過多體系統(tǒng)傳遞矩陣法計算得到導彈的動力響應，如圖5所示。分別通過多體系統(tǒng)傳遞矩陣法計算得到的導彈主動段的振動頻率變化曲線和通過自適應二階陷波濾波器級聯(lián)的方法估計得到的導彈主動段的振動頻率變化曲線的對比，如圖6所示。

圖5 通過多體系統(tǒng)傳遞矩陣法計算得到的動力響應Fig.5 The dynamic response of the missile calculated by the transfer matrix method of multibody system

圖6 分別用多體系統(tǒng)傳遞矩陣法計算得到的與用自適應二階陷波濾波器級聯(lián)的方法估計得到的導彈主動段振動頻率曲線的對比Fig.6 The vibration frequency curve of the missile in the boost phase obtained respectively by using the transfer matrix method of multibody system compared with the method of cascade of second-order adaptive notch filters

由圖6可清楚的看到，用本文提出的用自適應二階陷波濾波器級聯(lián)的方法估計得到的導彈振動頻率結(jié)果與多體系統(tǒng)傳遞矩陣法計算得到的導彈振動頻率結(jié)果非常接近。仿真結(jié)果表明，該方法對于導彈的振動頻率估計效果較好，符合預期效果。

3.2 實驗驗證

導彈在主動段、再入段飛行中受到發(fā)動機推力、噴氣噪聲以及紊流邊界層壓力等綜合作用產(chǎn)生的寬頻帶、大幅值振動激勵。測試表明某導彈飛行的振動環(huán)境所對應的頻率范圍為20～2 000 Hz。

為進一步驗證該算法的頻率估計效果，用振動臺模擬導彈振動的環(huán)境，做了捷聯(lián)于它的慣組的振動實驗，實驗現(xiàn)場如圖7所示。圖8為振動臺掃頻幅值時間歷程，圖9為振動臺20～2 000 Hz掃頻頻率估計曲線。導彈發(fā)射、飛行的整個過程，頻率范圍基本在20～2 000 Hz，因此選定振動臺掃頻的范圍為20～2 000 Hz。從圖9可看出來，該掃頻實驗的頻率估計結(jié)果與實際情況一致，估計曲線比較平滑。

圖7 捷聯(lián)慣組的振動臺實驗Fig.7 The shaking table experiment of the SIMU

圖8 振動臺掃頻幅值時間歷程Fig.8 The time history of the sweep frequency amplitude of the shaking table

4 結(jié) 論

基于多體系統(tǒng)傳遞矩陣法，從系統(tǒng)動力學角度建立了導彈系統(tǒng)的彈性振動模型，并計算出導彈系統(tǒng)的振動頻率以及動力響應。

圖9 振動臺20～2 000 Hz掃頻頻率估計曲線Fig.9 The estimation curve of sweep frequency of shaking table whose range is 20～2 000 Hz

針對直接固聯(lián)在導彈上的捷聯(lián)慣組處于導彈高過載、強沖擊、大振動復合動態(tài)環(huán)境中，輸出精度差的問題，設(shè)計了將二階自適應陷波濾波器級聯(lián)的方法來估計導彈彈性振動頻率，仿真和實驗結(jié)果表明，該方法對于導彈振動信號的頻率識別效果較好，為下一步捷聯(lián)慣組的振動控制提供了依據(jù)。本文的研究為提高捷聯(lián)慣組的使用精度提供了理論和技術(shù)參考。

[1] 劉洋.飛機發(fā)動機振動信號數(shù)字處理與分析技術(shù)研究[D].哈爾濱：哈爾濱工業(yè)大學，2008.

[2] TAN L, JIANG J. Novel adaptive IIR filter for frequency estimation and tracking[J]. IEEE Signal Processing Magazine,2009(11):186-189.

[3] 芮筱亭，贠來峰，陸毓琪，等.多體系統(tǒng)傳遞矩陣法及其應用[M].北京：科學出版社，2008.

[4] 張曉菲，劉振興，陳棟.頻率估計綜述[J].信息技術(shù),2011(11):58-62.

ZHANG Xiaofei,LIU Zhenxing, CHEN Dong. Summary of Frequency Estimation[J].Information Technology, 2011(11):58-62.

[5] RAO D V, KUNG S Y. Adaptive notch filtering for the retrieval of sinusoids in noise[J].Signal Processing, 1984,32(4): 791-802.

[6] REGALIA P A. An improved lattice-based adaptive IIR notch filter[J].Signal Processing,1991, 39(9):2124-2128.

[7] 單東升,張培強,吳耀武,等.基于遞推最小二乘法的炮控系統(tǒng)參數(shù)辨識仿真[J].系統(tǒng)仿真學報, 2013,25(8):1726-1729.

SHAN Dongsheng, ZHANG Peiqiang, WU Yaowu, et al. Simulation of parameter identification for gun control system based on RLS [J].Journal of System Simulation, 2013,25(8):1726-1729.

[8] 李湘清,孫秀霞,王棟,等.遞推最小二乘法在LQR參數(shù)調(diào)整中的應用[J].彈箭與制導學報，2007,27(4):99-101.

LI Xiangqing, SUN Xiuxia, WANG Dong, et al. LQR parameter regulate using the recursive least square[J]. Journal of Projectiles, Rockets, Missiles and Guidance, 2007, 27(4):99-101.

[9] CHOI H D, BANG H C. An adaptive control approach to the attitude control of a flexible rocket[J]. Control Engineering Practice, 2000, 8(9):1003-1010.

Adaptive estimation and simulation on the vibration frequency of missile system

ZHANG Heng,RUI Xiaoting,YANG Fufeng,CHEN Gangli

(Institute of Launch Dynamics, Nanjing University of Science & Technology, Nanjing 210094, China)

To improve the output accuracy of inertial measurement unit working in the complex dynamic environment of missiles which is of high acceleration, strong impact and severe vibration, a cascade of second-order adaptive notch filters was designed to estimate the elastic vibration frequency of missiles. The recursive least square method was adopted in the filters and a forgetting factor was introduced to optimize the real-time parameters of each notch filter. The filters possess the characteristics of low computational complexity and fast convergence speed. The results of simulation and experiment show that the method has a good frequency estimation effect on the vibration signals processing of missiles in complex dynamic environment. It has a strong anti-noise ability, and is suitable for the real-time calculation of missile computers.

adaptive notch filter; frequency estimation; recursive least square method

高校博士點基金博導類：復雜多體系統(tǒng)變參數(shù)自適應控制方法(20133219110037)

2016-04-01 修改稿收到日期： 2016-09-09

張衡 男，博士，1985年5月生

芮筱亭 男，博士，教授，博士生導師,1956年8月生 E-mail:ruixt@163.net

TJ760.1

A

10.13465/j.cnki.jvs.2017.10.034