王佳偉，祁克玉，楊愷華，梁 軻，閆 杰

（1. 西北工業(yè)大學 航天學院，西安 710072；2. 機電動態(tài)控制重點實驗室，西安 710065）

面對現(xiàn)代戰(zhàn)場武器裝備高精度和低效費比的雙重迫切需求，常規(guī)火炮不僅要延續(xù)其廣闊區(qū)域火力壓制的優(yōu)勢，還應具備在復雜環(huán)境中對小面積、低成本目標的精確打擊能力。二維彈道修正引信正是一種實現(xiàn)低成本精確打擊的前沿技術，通過更換彈道修正引信即可滿足常規(guī)無控彈藥智能化、靈巧化改造的迫切需求，是目前國際引信裝備領域最炙手可熱的發(fā)展方向。

基于中大口徑迫/榴彈應用平臺，采用慣性器件/衛(wèi)星定位組合測量彈體姿態(tài)及彈道參數(shù)，將修正引信與彈體通過解耦部件連接構成雙旋穩(wěn)定彈（Dual-Spin Stabilized Projectiles），其修正原理如圖1所示。

在氣動力的作用下，修正引信頭部的兩對翼面分別產(chǎn)生滾轉(zhuǎn)力矩和修正彈道所需控制力，這是目前國際上二維彈道修正引信技術中最具代表性的總體設計方案[1-2]。修正引信的實時滾轉(zhuǎn)角直接影響修正方向，是彈道控制模塊的重要輸入信息，其準確測量是二維彈道修正技術有待解決的關鍵技術。陀螺/衛(wèi)星定位組合測量方案作為滾轉(zhuǎn)角測量的一種主流方法，最早于1995年由David J. Lucia提出[3]。近年多有文獻進行進一步的闡述，其中HeeYoung Park等人繼承并改良了David J. Lucia的算法[4-9]，但這些方法都是在修正引信不大于3 r/s的微旋狀態(tài)且轉(zhuǎn)速恒定的假設條件下提出的。在實際彈道環(huán)境中，修正引信轉(zhuǎn)速從起控至修正末段的跨度可達數(shù)十轉(zhuǎn)/秒，遠遠大于文獻[4-9]所設定的理想情況，因此需要設計一種能夠適應較大轉(zhuǎn)速變化范圍的滾轉(zhuǎn)角測量方法以適應實際工程需求。

針對二維彈道修正引信全彈道轉(zhuǎn)速變化范圍較大的旋轉(zhuǎn)特點，本文建立了修正引信滾轉(zhuǎn)角與彈載陀螺/衛(wèi)星定位輸出的解析模型，將實測修正引信轉(zhuǎn)速作為補償引入系統(tǒng)量測方程，提出了一種基于擴展卡爾曼濾波（extended Kalman filter, EKF）的轉(zhuǎn)速補償滾轉(zhuǎn)角測量算法。

圖1 二維彈道修正原理圖Fig.1 Schematic of 2-D course correction

1 研究基礎

1.1 二維彈道修正原理

二維彈道修正引信通過解耦部件與彈丸主體連接，可以自由旋轉(zhuǎn)。引信頭部安裝四個斜置固定翼面，其中三個翼面帶有左向斜置角，另外一個翼面帶有右向斜置角。彈丸在飛行過程中，由于四個翼面均存在斜置角，都會產(chǎn)生一定的升力作用，其中：控制翼面1和2斜置角相反，其受到的氣動力方向相同，這兩個力對修正引信滾轉(zhuǎn)方向的力矩相互抵消，僅形成向上的氣動合力，即為彈道修正所需的控制力；導轉(zhuǎn)翼面1和2斜置角相同，受到的氣動力方向相反，力的作用相互抵消，卻形成促使修正引信旋轉(zhuǎn)的導轉(zhuǎn)力矩，其方向與彈丸自轉(zhuǎn)方向相反。當彈載控制器發(fā)出控制指令驅(qū)動滾轉(zhuǎn)執(zhí)行機構，使得修正翼面滾轉(zhuǎn)角度指向期望方向時，其產(chǎn)生的氣動升力即可改變彈丸原飛行軌跡，達到彈道修正的目的[10-11]。其原理如圖1所示。

1.2 擴展卡爾曼濾波

卡爾曼濾波器是針對觀測狀態(tài)變量的最優(yōu)線性濾波器，但對于復雜的非線性系統(tǒng)則需要用到一種非線性濾波方法，而EKF無疑是過去數(shù)十年應用最為廣泛的非線性狀態(tài)估計方法。

設k時刻帶有離散量測的連續(xù)時間系統(tǒng)如下所示：

其中，x為被估計的狀態(tài)變量，w為協(xié)方差為Q的系統(tǒng)噪聲，v為協(xié)方差為Rk的量測噪聲。如果記為狀態(tài)變量的估計值，則可得到系統(tǒng)偏微分矩陣在當前狀態(tài)估計的值為

此時狀態(tài)估計誤差的協(xié)方差Pk和傳播協(xié)方差表達式如下：

其中，A和L由式(2)給出。式(3)的積分方程將狀態(tài)變量的估計和估計誤差的協(xié)方差Pk由k-1時刻的后驗估計狀態(tài)傳播到k時刻的先驗估計狀態(tài)，隨后在每一個量測時刻按照離散卡爾曼濾波方程更新狀態(tài)估計和協(xié)方差，如式(4)所示：

至此得到k時刻的后驗狀態(tài)變量估計，濾波方程的第k步計算結(jié)束。

1.3 陀螺/衛(wèi)星定位組合算法

陀螺/衛(wèi)星定位組合解算滾轉(zhuǎn)角是利用衛(wèi)星定位接收機實時測量彈丸在發(fā)射坐標系下三個方向的速度分量，考慮到彈丸攻角通常很小可以省略，可以通過速度估算得到準彈體系下修正引信此刻的俯仰角速度˙和偏航角速度˙。結(jié)合彈體坐標系下固連于修正引信內(nèi)部的陀螺三個方向角速率輸出根據(jù)外彈道學相關理論，利用坐標變換建立彈體坐標系下陀螺角速率輸出與準彈體坐標系下俯仰/偏航角速率的物理解析關系，并從中解算得到滾轉(zhuǎn)角γ。算法詳見文獻[4]，原理框圖如圖2所示。

遺憾的是，該方法僅局限于不大于 3 r/s的微旋恒定轉(zhuǎn)速條件下的滾轉(zhuǎn)角測量，基于理論解析得到的滾轉(zhuǎn)角又極易受到量測噪聲的影響，因此需要設計一種濾波算法以實現(xiàn)轉(zhuǎn)速大范圍變化環(huán)境下的滾轉(zhuǎn)角準確測量。

圖2 陀螺/衛(wèi)星定位組合測量滾轉(zhuǎn)角算法原理Fig.2 Principle of roll estimation algorithm

2 濾波算法設計

以高速且時變轉(zhuǎn)速旋轉(zhuǎn)的修正引信會令傳統(tǒng)EKF方法中每一次濾波迭代周期內(nèi)滾轉(zhuǎn)角增量突變，導致狀態(tài)變量在濾波狀態(tài)更新過程中發(fā)散，最終使?jié)L轉(zhuǎn)角估計誤差逐漸增大。與傳統(tǒng)的EKF方法不同的是，轉(zhuǎn)速補償 EKF滾轉(zhuǎn)角測量方法在系統(tǒng)偏微分矩陣中增加了修正引信轉(zhuǎn)速信息，經(jīng)過積分得到單次迭代周期Δt內(nèi)的滾轉(zhuǎn)角增量Δγ，使其作為量測預估值kk H x的補償能夠最大限度地保證濾波算法迭代過程中滾轉(zhuǎn)角的預估值收斂。

待估狀態(tài)向量Xk包括轉(zhuǎn)速p、初始滾轉(zhuǎn)角γ0以及y、z軸陀螺測量誤差，則：

量測向量Zk為：

系統(tǒng)偏微分矩陣如式(7)～(10)所示：

3 仿真驗證

3.1 仿真條件

某型雙旋穩(wěn)定迫彈射角采用最大射程角，彈丸初速為名義初速，在炮兵標準氣象條件下全彈道彈丸穩(wěn)定飛行，外彈道曲線見圖3。

圖3 雙旋迫彈外彈道曲線Fig.3 Exterior trajectory curve of dual-spin mortar shell

圖 4是全彈道范圍彈丸和修正引信無控飛行的轉(zhuǎn)速變化曲線。彈丸發(fā)射后，迫彈尾翼產(chǎn)生導轉(zhuǎn)力矩使彈丸轉(zhuǎn)速由靜止狀態(tài)迅速提升至11 r/s左右，之后轉(zhuǎn)速隨動壓變化（由彈丸速度決定）呈倒“U”型變化趨勢且全彈道保持右旋（從彈尾向彈頭方向看去）。修正引信頭部安裝了與彈丸尾翼反向的導轉(zhuǎn)翼面，在反向?qū)мD(zhuǎn)力矩作用下轉(zhuǎn)速迅速提升至25 r/s左右，隨后在導轉(zhuǎn)力矩、摩擦力矩和滾轉(zhuǎn)阻尼力矩的合作用下自由旋轉(zhuǎn)，使其在全彈道內(nèi)與彈丸旋轉(zhuǎn)方向相反，即向左旋轉(zhuǎn)。

圖4 彈丸和引信全彈道轉(zhuǎn)速變化曲線Fig.4 Rotation speed curves of shell and fuze

假設北斗衛(wèi)星定位接收機在 10 s內(nèi)實現(xiàn)重捕定位，由圖4可以看到整個修正階段的引信轉(zhuǎn)速將在15 r/s至25 r/s之間波動，遠遠超過文獻[5-6]算法中微旋且恒定的轉(zhuǎn)速限制，因此需要一種既能滿足修正控制初期較高轉(zhuǎn)速條件，又能不失精度地滿足起控后微旋狀態(tài)下的修正引信滾轉(zhuǎn)角測量方法。表1給出了應用于滾轉(zhuǎn)角測量模型中的各測量參數(shù)。

表1 測量模型參數(shù)Tab.1 Parameters used in simulation

3.2 高轉(zhuǎn)速驗證

在3.1節(jié)所述仿真條件下，采用理論解析方法解算得到的滾轉(zhuǎn)角誤差，如圖5所示?？梢钥吹綕L轉(zhuǎn)角真值完全湮沒在噪聲中無法分辨，其誤差分布呈現(xiàn)出彈道中段較小而彈道初段及末段較大的趨勢，這是因為修正引信的轉(zhuǎn)速恰好在彈道中段最低，其帶來的一個仿真步長周期內(nèi)的滾轉(zhuǎn)角累積誤差也最小。這也側(cè)面印證了理論解析方法對于較高轉(zhuǎn)速條件下的滾轉(zhuǎn)角測量是完全失效的。

采用引入轉(zhuǎn)速補償?shù)?EKF算法對修正引信滾轉(zhuǎn)角進行優(yōu)化估計，圖6為使用濾波算法前后的滾轉(zhuǎn)角解算誤差對比，圖7為10 s后的解算誤差局部放大圖。

圖5 高轉(zhuǎn)速條件下理論解析法滾轉(zhuǎn)角解算誤差Fig.5 Roll estimation error by theoretical analytical method in the case of high rotation speed

可以看到采用引入轉(zhuǎn)速補償 EKF濾波算法的滾轉(zhuǎn)角估計值與真值之間的誤差迅速收斂，收斂后全彈道滾轉(zhuǎn)角解算絕對誤差不超過5°，經(jīng)統(tǒng)計誤差均值為-1.66°，方差為1.84°。宏觀來看，滾轉(zhuǎn)角解算誤差在彈道初段及末段有所增大，這是由于該時間段修正引信轉(zhuǎn)速較大而導致每一個仿真步長內(nèi)滾轉(zhuǎn)角增量過大，預測協(xié)方差和濾波器的增益不能隨新息同步改變從而失去對突變狀態(tài)變量的跟蹤能力，使估計誤差略有增加。以圖4為例，修正引信的轉(zhuǎn)速在彈道初段及末段超過了15 r/s，這意味著僅在一個步長周期內(nèi)預估狀態(tài)變量的增量就有 5.4°，因此采用轉(zhuǎn)速補償?shù)?EKF算法對修正引信滾轉(zhuǎn)角的估計效果是令人滿意的。

圖6 濾波前后的滾轉(zhuǎn)角解算誤差對比Fig.6 Comparison of roll estimation error with theoretical analytical method and EKF

圖7 誤差對比局部放大圖Fig.7 Magnification of roll estimation error

3.3 低轉(zhuǎn)速驗證

針對起控后彈道修正階段的低轉(zhuǎn)速環(huán)境，轉(zhuǎn)速補償?shù)腅KF算法同樣需要準確測量實時滾轉(zhuǎn)角。為了足夠真實地模擬修正引信較低轉(zhuǎn)速環(huán)境，而又不改變其在彈道飛行過程中所受升力和阻力等氣動條件，僅將彈道模型中修正引信的滾轉(zhuǎn)阻尼力矩系數(shù)增大，令其轉(zhuǎn)速降低至1～5 r/s。引入轉(zhuǎn)速補償?shù)腅KF滾轉(zhuǎn)角測量結(jié)果如圖8所示。

可以看到，引入轉(zhuǎn)速補償?shù)腅KF算法在低轉(zhuǎn)速條件下滾轉(zhuǎn)角測量精度依然表現(xiàn)出色，解算誤差迅速收斂且全彈道絕對誤差不大于 2.5°，經(jīng)統(tǒng)計解算誤差在收斂后的誤差均值為-0.34°，方差為0.35°。由此可見，轉(zhuǎn)速補償后的 EKF滾轉(zhuǎn)角測量算法有別于傳統(tǒng)理論解析方法，對于高速、低速且時變轉(zhuǎn)速條件下的滾轉(zhuǎn)角測量都具有很好的解算效果。

圖8 1～5 Hz轉(zhuǎn)速條件下滾轉(zhuǎn)角解算誤差Fig.8 Roll estimation error under low spinning rate

4 實驗驗證

以高精度MEMS傳感器三軸轉(zhuǎn)臺作為實驗平臺，驗證轉(zhuǎn)速補償EKF滾轉(zhuǎn)角測量算法，如圖9所示。

通過滾轉(zhuǎn)控制單元設置由30 r/s均勻減小至1 r/s的滾轉(zhuǎn)速率，固連在滾轉(zhuǎn)軸上的滾轉(zhuǎn)角測量硬件模塊會采集陀螺輸出信號，最后將采用轉(zhuǎn)速補償EKF滾轉(zhuǎn)角測量算法計算出的實時滾轉(zhuǎn)角與理論值進行對比，分析解算誤差。三軸轉(zhuǎn)臺滾轉(zhuǎn)角速率輸出以及對應角速率條件下滾轉(zhuǎn)角解算誤差曲線如圖10所示。

由圖10可以看到，在修正引信轉(zhuǎn)速自30 r/s至1 r/s范圍動態(tài)變化過程中，除了濾波迭代初期的振蕩之外，僅有的滾轉(zhuǎn)角解算失準皆是出現(xiàn)在轉(zhuǎn)速狀態(tài)突變所對應時刻（其收斂過程都不會超過 200ms），但轉(zhuǎn)速補償EKF滾轉(zhuǎn)角測量算法始終能夠保持較高的測量精度，且滾轉(zhuǎn)角解算絕對誤差不超過 4°。由此可見該算法完全能夠?qū)崿F(xiàn)大范圍轉(zhuǎn)速變化條件下的滾轉(zhuǎn)角準確測量。

圖9 MEMS三軸轉(zhuǎn)臺Fig.9 MEMS three-axis turntable

圖10 實驗室驗證滾轉(zhuǎn)角解算誤差曲線Fig.10 In-lab testing results with three-axis turntable

5 結(jié) 論

本文通過雙旋穩(wěn)定迫彈的外彈道模型分析了修正引信轉(zhuǎn)速在全彈道大范圍變化的轉(zhuǎn)動特性，提出了一種基于擴展卡爾曼濾波的轉(zhuǎn)速補償滾轉(zhuǎn)角測量方法。驗證結(jié)果為：1）在較高轉(zhuǎn)速仿真條件下，采用該估計方法500 ms內(nèi)就可以實現(xiàn)滾轉(zhuǎn)角解算快速收斂，全彈道滾轉(zhuǎn)角解算絕對誤差不大于5°，誤差均值為-1.66°，方差為1.84°；2）在低轉(zhuǎn)速仿真條件下，滾轉(zhuǎn)角解算誤差仍可迅速收斂且全彈道絕對誤差不大于 2.5°，收斂后的誤差均值為-0.34°，方差為0.35°；3）以三軸轉(zhuǎn)臺為實驗平臺進行驗證，該算法滾轉(zhuǎn)角解算絕對誤差不超過4°。以上驗證結(jié)果表明該算法完全可以滿足二維彈道修正引信系統(tǒng)中滾轉(zhuǎn)角測量精度的指標要求。

參考文獻（References）：

[1]Philippe W, Friedrich L, Lutz L, et al. Wind tunnel tests and open-loop trajectory simulations for a 155 mm canards guided spin stabilized projectile[C]//AIAA Atmospheric Flight Mechanics Conference and Exhibit. Honolulu,Hawaii, 2008: 881-898.

[2]Philippe W. Stability analysis for canard guided dual-spin stabilized projectiles[C]//AIAA Atmospheric Flight Mechanics Conference. Chicago, Illinois, 2009: 843-867.

[3]David J L. Estimation of the local vertical state for a guided munition shell with an embedded GPS/Micro-Mechanical inertial navigation system[D]. Cambridge,Massachusetts: MIT, 1995: 33-51.

[4]James M M. Efficient attitude estimation for a spin stabilized projectile[J]. Journal of Guidance, Control, and Dynamics, 2016, 39(2): 339-350.

[5]Park H Y, Kim K J, Lee J C, et al. Roll angle esti- mation for smart munition[C]//IFAC Symposium on Automatic Control in Aerospace. Toulouse, France: Pleiades Publishing. 2007: 6-12.

[6]Lee H S, Kim K J, Park H Y. Roll estimation of a smart munition using a magnetometer based on an unscented kalman filter[C]//AIAA Guidance Navigation and Control Conference and Exhibit. Honolulu, Hawaii: American Institute of Aeronautics and Astronautics. 2008:1-13.

[7]Adam R, Pecheur E, Bernard L, et al. In-flight roll angle estimation for a guided high spin projectile[C]//Sensors Applications Symposium. Catania, Italy, 2016: 308-313.

[8]王佳偉, 史凱, 徐國泰, 等. 基于擴展卡爾曼濾波的高轉(zhuǎn)速修正引信滾轉(zhuǎn)角測量方法[J]. 西北工業(yè)大學學報,2016, 34(6): 938-944.Wang J W, Shi K, Xu G T, et al. Roll estimation of high rotation speed correction fuze based on extended Kalman filter[J]. Journal of Northwestern Polytechnical University,2016, 34(6): 938-944.

[9]Thomas C, Stefer M, Bernard L. Roll angle estimation of a spinning object: experimental setup and preliminary results[C]// Proceedings of the International Conference on Microelectronics. Casablanca, Morocco, 2016: 1-4.

[10]Chang S J, Wang Z Y. Analysis of spin-rate property for dual-spin-stabilized projectiles with canards[J]. Journal of Spacecraft and Rockets, 2014, 51(3): 958-966.

[11]Julien S, Spilios T, Phillippe W. Flight control for a class of 155mm spin-stabilized projectile with reciprocating canards[C]//AIAA Guidance, Navigation, and Control Conference. Minneapolis, Minnesota, 2012: 685-694.