導彈系統振動頻率自適應估計方法與仿真

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

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

導彈系統振動頻率自適應估計方法與仿真

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

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

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

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

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

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

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

1 導彈彈性振動模型

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

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

(1)

(2)

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

(3)

這是一個常系數的齊次線性偏微分方程。令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為狀態矢量，a=[A,B,C,D]T

(7)

(8)

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

(9)

由各段梁的傳遞方程可以得到系統的總傳遞方程為

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

(10)

系統總傳遞矩陣為

(11)

系統的邊界條件為

(12)

代入式可以得到

(13)

刪去與零元素有關的行和列可以得到

(14)

式中

上式存在非零解，則

(15)

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

導彈的廣義坐標為

(16)

導彈的動力響應為

(17)

2 導彈振動頻率估計方法

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

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

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

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

(18)

式中，z是復數變量，n為陷波個數，當考慮單個陷波時n=1，上式變為

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

(19)

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

如圖4為級聯型的二階陷波濾波器的計算流程圖。

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

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

(20)

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

(21)

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

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

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

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

(22)

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

(23)

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

(24)

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

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

(25)

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

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

(26)

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

(27)

如式(28)～(30)顯示的遞推最小二乘[7]法是用來估計未知參數向量θ=[a1…ap]T。時間平均自相關φk(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階段的估計參數由下式給出

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

(30)

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

(31)

3 仿真結果及實驗驗證

3.1 仿真結果

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

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

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

圖6 分別用多體系統傳遞矩陣法計算得到的與用自適應二階陷波濾波器級聯的方法估計得到的導彈主動段振動頻率曲線的對比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可清楚的看到，用本文提出的用自適應二階陷波濾波器級聯的方法估計得到的導彈振動頻率結果與多體系統傳遞矩陣法計算得到的導彈振動頻率結果非常接近。仿真結果表明，該方法對于導彈的振動頻率估計效果較好，符合預期效果。

3.2 實驗驗證

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

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

圖7 捷聯慣組的振動臺實驗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 結 論

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

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

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

[1] 劉洋.飛機發動機振動信號數字處理與分析技術研究[D].哈爾濱：哈爾濱工業大學，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] 芮筱亭，贠來峰，陸毓琪，等.多體系統傳遞矩陣法及其應用[M].北京：科學出版社，2008.

[4] 張曉菲，劉振興，陳棟.頻率估計綜述[J].信息技術,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] 單東升,張培強,吳耀武,等.基于遞推最小二乘法的炮控系統參數辨識仿真[J].系統仿真學報, 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參數調整中的應用[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

高校博士點基金博導類：復雜多體系統變參數自適應控制方法(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