簡(jiǎn)介比例導(dǎo)引三維彈道仿真是空對(duì)空導(dǎo)彈制導(dǎo)研究中的經(jīng)典課題這份壓縮包提供了基于MATLAB與龍格庫(kù)塔算法的完整三維彈道仿真實(shí)現(xiàn)方案面向?qū)椫茖?dǎo)設(shè)計(jì)與網(wǎng)絡(luò)攻防技術(shù)研發(fā)人員也適合相關(guān)研究者系統(tǒng)學(xué)習(xí)比例導(dǎo)引律建模、微分方程數(shù)值求解與彈道軌跡分析方法。壓縮包共589個(gè)文件、約7.11MB以481個(gè)m格式源碼文件、26個(gè)mat數(shù)據(jù)文件和20個(gè)fig圖形文件為主體另有少量C/C輔助程序與PDF說(shuō)明文檔源碼覆蓋導(dǎo)引律核心算法、彈道解算與參數(shù)配置數(shù)據(jù)文件支持多工況仿真結(jié)果對(duì)比圖形文件直接呈現(xiàn)三維彈道可視化效果。該資源已有179人學(xué)習(xí)通過(guò)調(diào)整導(dǎo)彈速度、目標(biāo)加速度等參數(shù)可直觀復(fù)現(xiàn)不同攔截場(chǎng)景進(jìn)而評(píng)估比例導(dǎo)引法在復(fù)雜機(jī)動(dòng)目標(biāo)下的精度與穩(wěn)定性。資源還探索了比例導(dǎo)引思路向網(wǎng)絡(luò)安全防御的遷移為應(yīng)對(duì)快速變化的安全威脅提供了新穎的對(duì)抗策略。1. 攻擊水平機(jī)動(dòng)目標(biāo)的比例導(dǎo)引三維彈道仿真這套MATLAB工程解決什么問(wèn)題一個(gè)攔截彈在末端遭遇水平面內(nèi)連續(xù)轉(zhuǎn)彎的目標(biāo)時(shí)很多人在二維平面里調(diào)好的比例導(dǎo)引參數(shù)拿到三維空間會(huì)直接失效——這不是算法錯(cuò)了而是視線幾何在三維空間里發(fā)生了耦合。標(biāo)題里的“攻擊水平機(jī)動(dòng)目標(biāo)比例導(dǎo)引三維彈道仿真”要做的正是把比例導(dǎo)引律、龍格庫(kù)塔算法、三維彈道仿真三件事在MATLAB里串成一條可復(fù)現(xiàn)的鏈路建立導(dǎo)彈與目標(biāo)的相對(duì)運(yùn)動(dòng)方程用龍格庫(kù)塔推進(jìn)彈道微分方程最終量化評(píng)估脫靶量和需用過(guò)載。這套方案適合做制導(dǎo)控制課程設(shè)計(jì)、彈道方案預(yù)研或算法對(duì)比驗(yàn)證的讀者解決的是“公式看得懂、代碼寫(xiě)不出、結(jié)果不敢信”的斷層問(wèn)題。讀完這篇文章你可以直接照著架構(gòu)搭出自己的三維攔截仿真模型。2. 三維彈道運(yùn)動(dòng)學(xué)與比例導(dǎo)引律建模先把坐標(biāo)系和制導(dǎo)指令算對(duì)仿真里最容易出問(wèn)題的不是積分器而是狀態(tài)量怎么選、制導(dǎo)指令怎么算。這兩件事決定了后面所有代碼的形態(tài)值得先花整章講清楚。2.1 狀態(tài)量選取與坐標(biāo)系約定我一般使用慣性坐標(biāo)系下的位置和速度作為核心狀態(tài)量不用彈道傾角和彈道偏角。原因是歐拉角描述在鉛垂方向接近正負(fù)90度時(shí)會(huì)出現(xiàn)奇異點(diǎn)而三維攔截彈道中導(dǎo)彈和目標(biāo)的相對(duì)幾何經(jīng)常穿越這些角度區(qū)域一旦角度跳變光線角速率會(huì)被污染成脈沖尖刺。這里的狀態(tài)量是一個(gè)15維向量按順序排列為導(dǎo)彈位置3維、導(dǎo)彈速度3維、導(dǎo)彈加速度3維、目標(biāo)位置3維、目標(biāo)速度3維。導(dǎo)彈加速度作為狀態(tài)量而不是直接等于制導(dǎo)指令是為了表達(dá)駕駛儀的一階慣性延遲。實(shí)際飛行器中指令加速度不可能瞬時(shí)建立用一階慣性環(huán)節(jié)近似是工程常見(jiàn)做法時(shí)間常數(shù)tau通常在0.1到0.5秒之間這個(gè)參數(shù)對(duì)末端脫靶量影響非常顯著。狀態(tài)排列約定如下表所示后續(xù)所有代碼都按這個(gè)索引取值。狀態(tài)索引物理含義符號(hào)1~3導(dǎo)彈位置Rm4~6導(dǎo)彈速度Vm7~9導(dǎo)彈加速度Am10~12目標(biāo)位置Rt13~15目標(biāo)速度Vt坐標(biāo)系的約定是x為水平前向y為高度方向z為水平側(cè)向。目標(biāo)做水平機(jī)動(dòng)時(shí)y方向速度始終為0機(jī)動(dòng)發(fā)生在x-z平面內(nèi)這樣“水平機(jī)動(dòng)”在模型里就有了明確的數(shù)學(xué)約束也方便后面對(duì)仿真結(jié)果做降維驗(yàn)證。2.2 比例導(dǎo)引律的矢量形式視線角速率是怎么驅(qū)動(dòng)過(guò)載指令的比例導(dǎo)引的核心思想不是“朝目標(biāo)當(dāng)前位置飛”而是把視線角速率壓到零。目標(biāo)只要改變運(yùn)動(dòng)方向視線就會(huì)轉(zhuǎn)動(dòng)導(dǎo)引律隨即產(chǎn)生與視線角速率成正比的過(guò)載指令等效于一個(gè)比例控制器。這個(gè)邏輯對(duì)機(jī)動(dòng)目標(biāo)天然有自適應(yīng)能力是它成為戰(zhàn)術(shù)導(dǎo)彈最常見(jiàn)制導(dǎo)律的根本原因。在三維空間里比例導(dǎo)引的標(biāo)量形式不再適用需要采用矢量形式。視線向量由目標(biāo)位置減導(dǎo)彈位置得到視線角速率向量由相對(duì)位置叉乘相對(duì)速度得到制導(dǎo)指令進(jìn)一步由視線角速率叉乘相對(duì)速度生成。這個(gè)矢量形式避免了把制導(dǎo)問(wèn)題拆成兩個(gè)平面的做法保留了縱向和側(cè)向通道的幾何耦合。function a_cmd proportional_navigation(Rm, Vm, Rt, Vt, N) % 比例導(dǎo)引指令矢量形式單位 m/s^2 R_vec Rt - Rm; % 視線向量 R_norm norm(R_vec); V_rel Vt - Vm; % 相對(duì)速度 u_R R_vec / R_norm; % 視線單位向量 omega cross(R_vec, V_rel) / (R_norm^2); % 視線角速率向量 V_c -dot(V_rel, u_R); % 接近速度 if V_c 0 a_cmd zeros(3,1); % 目標(biāo)遠(yuǎn)離時(shí)停止制導(dǎo) return; end a_cmd N * V_c * cross(omega, u_R); % 指令加速度 end這里的關(guān)鍵參數(shù)是導(dǎo)航比N工程常見(jiàn)取值是3到5。N偏小時(shí)末端需用過(guò)載大彈道比較彎曲N偏大時(shí)初始段指令偏大容易觸發(fā)過(guò)載飽和。vector形式下不需要單獨(dú)處理歐拉角也就繞開(kāi)了角度跳變和象限判斷這些玄學(xué)問(wèn)題。V_c小于0的情況在迎頭攔截場(chǎng)景幾乎不會(huì)出現(xiàn)但防御性判斷必須寫(xiě)上否則目標(biāo)一旦掉頭遠(yuǎn)離仿真會(huì)進(jìn)入導(dǎo)彈反向飛行的錯(cuò)誤狀態(tài)。2.3 水平機(jī)動(dòng)目標(biāo)建模轉(zhuǎn)彎?rùn)C(jī)動(dòng)與蛇形機(jī)動(dòng)目標(biāo)水平機(jī)動(dòng)的建模方式直接決定仿真的說(shuō)服力。最簡(jiǎn)單的是勻速直線目標(biāo)用來(lái)做模型驗(yàn)證實(shí)用的是水平轉(zhuǎn)彎目標(biāo)用來(lái)評(píng)估比例導(dǎo)引對(duì)抗持續(xù)機(jī)動(dòng)的能力。水平轉(zhuǎn)彎的含義是目標(biāo)速度始終在x-z平面內(nèi)旋轉(zhuǎn)高度不變速度大小恒定。function a_t target_level_maneuver(Vt, t, mode) % 目標(biāo)水平機(jī)動(dòng)加速度作用在水平面內(nèi) psi_t atan2(Vt(3), Vt(1)); % 目標(biāo)速度方位角 if strcmp(mode, const) w_t 0.15; % 恒定轉(zhuǎn)彎角速率 rad/s else w_t 0.2 * sin(0.5 * t); % 蛇形機(jī)動(dòng)角速率 end a_t norm(Vt) * w_t * [-sin(psi_t); 0; cos(psi_t)]; end這個(gè)函數(shù)的本質(zhì)是給目標(biāo)施加一個(gè)水平面內(nèi)向心加速度。目標(biāo)速度方向由方位角psi_t確定加速度方向垂直于速度方向且保持在水平面內(nèi)這樣就實(shí)現(xiàn)了“水平轉(zhuǎn)彎”而不是“爬升轉(zhuǎn)彎”。蛇形機(jī)動(dòng)模式讓角速率隨時(shí)間正弦變化更貼近真實(shí)目標(biāo)規(guī)避時(shí)的連續(xù)變向。建議做兩組仿真對(duì)照一組恒定轉(zhuǎn)彎率一組蛇形機(jī)動(dòng)前者用來(lái)測(cè)需用過(guò)載的包線后者用來(lái)測(cè)制導(dǎo)律的魯棒性。初始場(chǎng)景參數(shù)按典型超聲速攔截彈設(shè)置如下表。參數(shù)值說(shuō)明導(dǎo)彈初始位置(0, 3000, 0) m與目標(biāo)存在高度差彈道三維化導(dǎo)彈初始速度(800, 0, 0) m/s超聲速攔截彈典型速度目標(biāo)初始位置(8000, 4000, 0) m前方高空目標(biāo)初始速度(-250, 0, 0) m/s亞聲速迎頭接近目標(biāo)轉(zhuǎn)彎角速率0.15 rad/s水平機(jī)動(dòng)強(qiáng)度導(dǎo)航比 N4比例導(dǎo)引常數(shù)駕駛儀時(shí)間常數(shù)0.2 s一階慣性延遲這套參數(shù)下導(dǎo)彈初始高度低于目標(biāo)且航向水平目標(biāo)前方迎頭飛行并持續(xù)水平轉(zhuǎn)彎導(dǎo)彈必須在爬升同時(shí)完成側(cè)向轉(zhuǎn)彎彈道自然呈現(xiàn)三維形態(tài)適合檢驗(yàn)三維比例導(dǎo)引的實(shí)際效果。3. 龍格庫(kù)塔算法求解彈道微分方程MATLAB里如何把連續(xù)模型推進(jìn)成軌跡比例導(dǎo)引給出的是加速度指令彈道需要從微分方程積分出來(lái)。對(duì)這套模型來(lái)說(shuō)積分器不是隨便選一個(gè)就能用的四階龍格庫(kù)塔算法是精度和實(shí)現(xiàn)復(fù)雜度之間的平衡點(diǎn)也是標(biāo)題明確指定的核心算法。3.1 合并后的彈道微分方程組把第2章的各個(gè)模塊合并后整條彈道的狀態(tài)方程是15個(gè)一階常微分方程組成的方程組。導(dǎo)彈位置導(dǎo)數(shù)等于導(dǎo)彈速度導(dǎo)彈速度導(dǎo)數(shù)等于當(dāng)前駕駛儀實(shí)際加速度導(dǎo)彈加速度導(dǎo)數(shù)由一階慣性延遲方程給出。目標(biāo)側(cè)的位置導(dǎo)數(shù)等于目標(biāo)速度目標(biāo)速度導(dǎo)數(shù)等于水平機(jī)動(dòng)加速度。這條方程組沒(méi)有解析解必須數(shù)值積分。龍格庫(kù)塔算法的思路是在一個(gè)積分步內(nèi)取多個(gè)中間點(diǎn)上的導(dǎo)數(shù)加權(quán)平均后推進(jìn)狀態(tài)四個(gè)階段分別對(duì)應(yīng)步長(zhǎng)起點(diǎn)的導(dǎo)數(shù)、兩個(gè)半步長(zhǎng)中間點(diǎn)的導(dǎo)數(shù)、以及步長(zhǎng)終點(diǎn)的導(dǎo)數(shù)。四階意味著局部截?cái)嗾`差是步長(zhǎng)的五次方對(duì)彈道仿真這種跨幾十秒的積分來(lái)說(shuō)精度充分。需要在微分方程內(nèi)部重新計(jì)算制導(dǎo)指令這是個(gè)容易忽略但很關(guān)鍵的細(xì)節(jié)。RK4的k2和k3階段會(huì)用到半步長(zhǎng)處的狀態(tài)值如果制導(dǎo)指令只在步長(zhǎng)起點(diǎn)算一次然后保持不變快速變化的視線角速率在末端會(huì)被嚴(yán)重低估。3.2 RK4步進(jìn)函數(shù)最簡(jiǎn)實(shí)現(xiàn)與調(diào)用約定function [t_next, X_next] rk4_step(t, X, h, f) % 四階龍格庫(kù)塔單步推進(jìn) k1 f(t, X); k2 f(t h/2, X h/2 * k1); k3 f(t h/2, X h/2 * k2); k4 f(t h, X h * k3); X_next X h/6 * (k1 2*k2 2*k3 k4); t_next t h; end這個(gè)函數(shù)是純數(shù)值方法不關(guān)心狀態(tài)量是什么物理含義每個(gè)階段都調(diào)用傳入的函數(shù)句柄f來(lái)計(jì)算狀態(tài)導(dǎo)數(shù)。k1到k4依次使用越來(lái)越靠后的時(shí)間點(diǎn)估計(jì)導(dǎo)數(shù)加權(quán)系數(shù)1/6、2/6、2/6、1/6滿足積分公式的精度條件。調(diào)用時(shí)f的寫(xiě)法是匿名的把導(dǎo)彈模型的參數(shù)像N、tau、a_max一次性捕獲進(jìn)去這樣積分器不需要知道模型內(nèi)部的細(xì)節(jié)。血淚經(jīng)驗(yàn)是別在彈道仿真里默認(rèn)用ode45當(dāng)黑匣子。ode45是變步長(zhǎng)算法在末端相對(duì)距離快速變化時(shí)它的誤差控制在突然變小的步長(zhǎng)上會(huì)導(dǎo)致步數(shù)激增而且默認(rèn)容差對(duì)攔截彈道不夠緊。固定步長(zhǎng)RK4的輸出時(shí)間點(diǎn)完全可控后面做步長(zhǎng)收斂性掃描和脫靶量拋物線插值時(shí)都要依賴等間距采樣這個(gè)優(yōu)勢(shì)在驗(yàn)證階段會(huì)體現(xiàn)得很明顯。3.3 仿真主循環(huán)狀態(tài)更新、終止條件與數(shù)據(jù)記錄% 參數(shù)與初始狀態(tài) N 4; tau 0.2; h 0.01; t_end 60; a_max 30 * 9.8; % 可用過(guò)載單位 m/s^2 X zeros(15,1); X(1:3) [0; 3000; 0]; % 導(dǎo)彈位置 X(4:6) [800; 0; 0]; % 導(dǎo)彈速度 X(7:9) [0; 0; 0]; % 導(dǎo)彈加速度初值 X(10:12) [8000; 4000; 0]; % 目標(biāo)位置 X(13:15) [-250; 0; 0]; % 目標(biāo)速度 % 預(yù)分配記錄數(shù)組 n_max ceil(t_end / h) 1; hist_t zeros(n_max,1); hist_X zeros(n_max,15); t 0; idx 0; while t t_end idx idx 1; hist_t(idx) t; hist_X(idx,:) X.; R_vec X(10:12) - X(1:3); R_norm norm(R_vec); if R_norm 10 || R_norm 30000 break; % 命中或飛散后終止 end X rk4_step(t, X, h, (tt,xx) eom_missile(tt, xx, N, tau, a_max, const)); t t h; end hist_t hist_t(1:idx); hist_X hist_X(1:idx,:);這個(gè)主循環(huán)用固定步長(zhǎng)0.01秒推進(jìn)對(duì)速度差約1050米每秒的迎頭場(chǎng)景每個(gè)積分步內(nèi)導(dǎo)彈和目標(biāo)相對(duì)位置變化約10米制導(dǎo)指令的刷新頻率足夠。終止條件有兩個(gè)相對(duì)距離小于10米視為命中相對(duì)距離超過(guò)30000米說(shuō)明彈道發(fā)散提前退出。預(yù)分配記錄數(shù)組是個(gè)好習(xí)慣MATLAB里循環(huán)內(nèi)動(dòng)態(tài)增長(zhǎng)數(shù)組會(huì)反復(fù)申請(qǐng)內(nèi)存幾千步仿真感知不明顯但要跑參數(shù)掃描時(shí)差距就出來(lái)了。微分方程函數(shù)eom_missile負(fù)責(zé)把制導(dǎo)指令、目標(biāo)機(jī)動(dòng)和狀態(tài)導(dǎo)數(shù)組合在一起完整實(shí)現(xiàn)如下function dX eom_missile(t, X, N, tau, a_max, mode) Rm X(1:3); Vm X(4:6); Am X(7:9); Rt X(10:12); Vt X(13:15); % 制導(dǎo)指令微分方程內(nèi)部重新計(jì)算 a_cmd proportional_navigation(Rm, Vm, Rt, Vt, N); if norm(a_cmd) a_max a_cmd a_cmd / norm(a_cmd) * a_max; % 過(guò)載限幅 end % 目標(biāo)水平機(jī)動(dòng) a_t target_level_maneuver(Vt, t, mode); dX zeros(15,1); dX(1:3) Vm; dX(4:6) Am; dX(7:9) (a_cmd - Am) / tau; % 駕駛儀一階延遲 dX(10:12) Vt; dX(13:15) a_t; end過(guò)載限幅放在制導(dǎo)指令之后體現(xiàn)了真實(shí)飛行器的物理約束。a_max取30g大約294米每平方秒這是中遠(yuǎn)程攔截彈的常見(jiàn)過(guò)載指標(biāo)。如果不限幅仿真會(huì)在目標(biāo)機(jī)動(dòng)較強(qiáng)的場(chǎng)景里給出一個(gè)脫靶量很小的漂亮結(jié)果但那個(gè)結(jié)果建立在導(dǎo)彈能輸出上百g過(guò)載的假設(shè)上實(shí)際不可實(shí)現(xiàn)。4. 三維彈道仿真的MATLAB工程結(jié)構(gòu)腳本組織、可視化和評(píng)估指標(biāo)模型能跑通之后接下來(lái)是工程化的問(wèn)題。這個(gè)標(biāo)題本質(zhì)上是仿真建模工作代碼組織是否清晰直接決定了參數(shù)掃描和算法對(duì)比階段的工作效率。4.1 工程文件劃分與數(shù)據(jù)流我一般按函數(shù)職責(zé)拆成五個(gè)文件不搞大而全的單一腳本。每個(gè)函數(shù)只做一件事錯(cuò)誤定位和參數(shù)修改都方便。文件職責(zé)關(guān)鍵函數(shù)簽名main_sim.m參數(shù)設(shè)置、主循環(huán)、數(shù)據(jù)記錄無(wú)腳本eom_missile.m彈道狀態(tài)方程dX eom_missile(t, X, N, tau, a_max, mode)proportional_navigation.m比例導(dǎo)引指令a_cmd proportional_navigation(Rm, Vm, Rt, Vt, N)target_level_maneuver.m目標(biāo)水平機(jī)動(dòng)模型a_t target_level_maneuver(Vt, t, mode)rk4_step.m四階龍格庫(kù)塔單步[t_next, X_next] rk4_step(t, X, h, f)plot_trajectory.m三維彈道可視化plot_trajectory(hist_t, hist_X)數(shù)據(jù)流是單向的主循環(huán)持有當(dāng)前狀態(tài)X調(diào)用rk4_steprk4_step內(nèi)部多次調(diào)用eom_missileeom_missile內(nèi)部調(diào)用比例導(dǎo)引函數(shù)和目標(biāo)機(jī)動(dòng)函數(shù)。這樣分層后替換目標(biāo)機(jī)動(dòng)模型或者修改駕駛儀模型都不會(huì)牽動(dòng)積分器和主循環(huán)。4.2 三維彈道與目標(biāo)軌跡的可視化實(shí)現(xiàn)function plot_trajectory(hist_t, hist_X) figure(Color,w); hold on; grid on; box on; plot3(hist_X(:,1), hist_X(:,2), hist_X(:,3), b-, LineWidth, 1.6); plot3(hist_X(:,10), hist_X(:,11), hist_X(:,12), r--, LineWidth, 1.4); % 每隔 200 步畫(huà)一次導(dǎo)彈速度箭頭 idx_vec 1:200:size(hist_X,1); quiver3(hist_X(idx_vec,1), hist_X(idx_vec,2), hist_X(idx_vec,3), ... hist_X(idx_vec,4), hist_X(idx_vec,5), hist_X(idx_vec,6), ... Color, [0 0.45 0.74]); xlabel(x (m)); ylabel(y (m)); zlabel(z (m)); legend(導(dǎo)彈彈道,目標(biāo)軌跡,導(dǎo)彈速度,Location,best); axis equal; view(3); endplot3畫(huà)三維軌跡線quiver3在彈道上按固定間隔疊加速度矢量箭頭能直觀看出導(dǎo)彈速度方向的變化速率。axis equal保證三個(gè)軸比例一致否則垂直方向被自動(dòng)拉伸后會(huì)嚴(yán)重誤導(dǎo)彈道曲率的判斷。view(3)給出默認(rèn)三維視角配合rotate3d可交互旋轉(zhuǎn)。對(duì)追逐場(chǎng)景來(lái)說(shuō)如果導(dǎo)彈速度箭頭始終指向目標(biāo)當(dāng)前位置說(shuō)明比例導(dǎo)引實(shí)際上被寫(xiě)成了追蹤法需要回查視線角速率計(jì)算是否正確。4.3 脫靶量與需用過(guò)載評(píng)估可視化只能定性判斷定量評(píng)估需要計(jì)算指標(biāo)。核心指標(biāo)是脫靶量、需用過(guò)載峰值、飛行時(shí)間和視線角速率峰值。指標(biāo)計(jì)算方式工程意義脫靶量相對(duì)距離序列最小值制導(dǎo)精度需用過(guò)載峰值指令加速度最大值 / g機(jī)動(dòng)需求是否超限飛行時(shí)間仿真的終止時(shí)刻攔截窗口視線角速率峰值視線角速率向量模最大值導(dǎo)引頭跟蹤能力約束脫靶量的基礎(chǔ)計(jì)算是取相對(duì)距離序列最小值一行代碼即可。但固定步長(zhǎng)下采樣點(diǎn)可能正好錯(cuò)過(guò)真實(shí)的最近距離點(diǎn)尤其是接近速度大于1000米每秒時(shí)0.01秒步長(zhǎng)意味著相鄰采樣點(diǎn)相差10米。要更精確地估計(jì)脫靶量需要做拋物線插值這個(gè)技巧放到最后一章專門(mén)展開(kāi)。5. 比例導(dǎo)引三維彈道仿真的避坑指南五個(gè)讓結(jié)果失真的常見(jiàn)問(wèn)題做這類仿真最大的問(wèn)題不是代碼跑不通而是跑通了但結(jié)果不可信。下面五個(gè)坑是我自己踩過(guò)、也在幫別人排查時(shí)反復(fù)見(jiàn)過(guò)的按現(xiàn)象、原因、解決的順序?qū)懬宄?.1 現(xiàn)象脫靶量隨步長(zhǎng)減小反而增大彈道出現(xiàn)鋸齒形抖動(dòng)原因如果用了歐拉法或者ode45默認(rèn)容差末端相對(duì)距離快速變化時(shí)數(shù)值誤差主導(dǎo)結(jié)果。更隱蔽的原因是制導(dǎo)指令在積分步外只算了一次RK4中間階段用的都是舊視線信息末端一個(gè)步長(zhǎng)內(nèi)視線角速率可能變化幾十個(gè)百分點(diǎn)舊指令自然產(chǎn)生系統(tǒng)性偏差。解決制導(dǎo)指令計(jì)算放在微分方程函數(shù)內(nèi)部讓RK4的每個(gè)階段都基于當(dāng)前狀態(tài)重新算指令。然后做步長(zhǎng)收斂性掃描分別跑h0.1、0.05、0.01秒對(duì)比脫靶量和需用過(guò)載峰值。如果結(jié)果隨步長(zhǎng)顯著變化繼續(xù)縮小步長(zhǎng)直到曲線重合。固定步長(zhǎng)RK4的好處在這里體現(xiàn)得最充分。5.2 現(xiàn)象彈道末端過(guò)載出現(xiàn)尖峰指令直接頂?shù)较薹翟蛳鄬?duì)距離趨近于零時(shí)視線角速率的計(jì)算公式里分母是R的平方R越小角速率增長(zhǎng)越快乘上不斷增大的接近速度后指令爆炸是數(shù)學(xué)上的必然。這不是比例導(dǎo)引的問(wèn)題而是末端幾何本身固有的奇異效應(yīng)。解決工程常見(jiàn)的處理是末端切換。相對(duì)距離小于某個(gè)閾值時(shí)凍結(jié)制導(dǎo)指令或者切換到比例導(dǎo)引的末端直線彈道。閾值一般取100到200米也可以用代碼里的過(guò)載限幅來(lái)兜底。我建議限幅和凍結(jié)同時(shí)做限幅防數(shù)值爆炸凍結(jié)防指令抖動(dòng)兩條合起來(lái)彈道才穩(wěn)定。5.3 現(xiàn)象視線角速率序列里出現(xiàn)孤立尖峰彈道沒(méi)有明顯擾動(dòng)原因這個(gè)現(xiàn)象幾乎都是因?yàn)橛脷W拉角表示視線方向再通過(guò)差分或者解析求導(dǎo)得到角速率。atan2在正負(fù)π交界處會(huì)跳變asin在正負(fù)90度附近存在多值問(wèn)題即使角度序列看起來(lái)連續(xù)差分后也會(huì)產(chǎn)生幅度極大的偽角速率。解決換用第2章的矢量叉積法用相對(duì)位置叉乘相對(duì)速度直接得到視線角速率向量全程不出現(xiàn)歐拉角尖峰問(wèn)題從根上消失。如果出于某種原因必須保留歐拉角輸出至少要對(duì)角度序列做unwrap處理再把角速率超過(guò)物理合理的值判定為野值剔除。5.4 現(xiàn)象目標(biāo)水平機(jī)動(dòng)時(shí)拆成兩個(gè)平面分別仿真的結(jié)果與三維仿真偏差很大原因把三維問(wèn)題拆成縱向平面和側(cè)向平面獨(dú)立做比例導(dǎo)引忽略了兩者通過(guò)視線幾何產(chǎn)生的耦合。目標(biāo)水平轉(zhuǎn)彎時(shí)視線在空間內(nèi)持續(xù)旋轉(zhuǎn)縱向平面的視線角速率實(shí)際上包含了側(cè)向運(yùn)動(dòng)的影響分開(kāi)算等于人為切斷了這個(gè)耦合通道。解決直接用矢量形式的三維比例導(dǎo)引不要自己拆平面。如果模型結(jié)構(gòu)上必須分開(kāi)至少要把側(cè)向平面的視線旋轉(zhuǎn)信息反饋到縱向平面的接近速度修正里。但這樣做代碼復(fù)雜度和出錯(cuò)概率遠(yuǎn)高于直接三維計(jì)算得不償失。5.5 現(xiàn)象MATLAB腳本里的中文注釋變成亂碼舊工程在新機(jī)器上打不開(kāi)原因MATLAB編輯器編碼經(jīng)歷了從本地編碼向UTF-8遷移的過(guò)程。舊版本的默認(rèn)編碼在中文Windows下通常是GBK而新版本默認(rèn)UTF-8用新版本打開(kāi)GBK編碼的.m文件時(shí)所有中文字符按UTF-8解析自然變成亂碼。這個(gè)問(wèn)題和彈道模型無(wú)關(guān)但遇到時(shí)非常影響效率。解決新工程統(tǒng)一使用UTF-8編碼并在保存時(shí)確認(rèn)編輯器編碼設(shè)置。老文件可以用fileread配合native2unicode按GBK讀出再轉(zhuǎn)成UTF-8重新保存。代碼里盡量不用中文變量名注釋亂碼不影響執(zhí)行但變量名亂碼會(huì)導(dǎo)致整個(gè)腳本無(wú)法運(yùn)行。6. 用降維驗(yàn)證與參數(shù)掃描收尾讓三維比例導(dǎo)引仿真可信模型做完不等于模型是對(duì)的。在做任何參數(shù)分析之前我習(xí)慣先跑兩個(gè)低成本的驗(yàn)證實(shí)驗(yàn)這兩個(gè)實(shí)驗(yàn)?zāi)芎Y掉大部分實(shí)現(xiàn)錯(cuò)誤。第一個(gè)是冒煙測(cè)試把導(dǎo)航比N設(shè)為0此時(shí)比例導(dǎo)引不產(chǎn)生任何過(guò)載指令導(dǎo)彈應(yīng)該沿初始速度方向直線飛行。如果彈道彎曲了說(shuō)明狀態(tài)量拼接順序錯(cuò)誤或者目標(biāo)加速度被錯(cuò)誤混入導(dǎo)彈狀態(tài)。第二個(gè)是降維驗(yàn)證把目標(biāo)機(jī)動(dòng)關(guān)掉初始條件限制在單一平面內(nèi)比如把y方向的初始位置和速度全部清零此時(shí)三維矢量比例導(dǎo)引應(yīng)該退化為經(jīng)典二維比例導(dǎo)引脫靶量隨N的變化趨勢(shì)與教科書(shū)一致。這兩個(gè)驗(yàn)證通過(guò)后再打開(kāi)目標(biāo)水平機(jī)動(dòng)你看到的三維效應(yīng)才是真實(shí)可信的。參數(shù)掃描是另一個(gè)值得做的檢驗(yàn)。對(duì)導(dǎo)航比N取3、4、5對(duì)步長(zhǎng)h取0.1、0.05、0.01分別記錄脫靶量和需用過(guò)載峰值。你會(huì)看到N增加時(shí)脫靶量先下降后上升N在某個(gè)中間值最優(yōu)同時(shí)需用過(guò)載峰值單調(diào)上升。如果N增大時(shí)脫靶量單調(diào)下降且過(guò)載峰值不變說(shuō)明某個(gè)環(huán)節(jié)丟失了物理約束。脫靶量的精確估計(jì)建議用三點(diǎn)拋物線插值補(bǔ)上固定步長(zhǎng)采樣可能錯(cuò)過(guò)真實(shí)最近點(diǎn)的問(wèn)題。代碼很短直接在記錄數(shù)組上操作。R_hist vecnorm(hist_X(:,10:12) - hist_X(:,1:3), 2, 2); [~, i_min] min(R_hist); if i_min 2 i_min length(R_hist) - 1 tt hist_t(i_min-1 : i_min1); rr R_hist(i_min-1 : i_min1); p polyfit(tt, rr, 2); % 二次多項(xiàng)式擬合 t_miss -p(2) / (2 * p(1)); % 拋物線頂點(diǎn)時(shí)間 miss polyval(p, t_miss); % 插值脫靶量 else miss R_hist(i_min); endvecnorm函數(shù)需要MATLAB R2017b以上版本舊版本可以用sqrt(sum(R_hist.^2, 2))替代。插值的前提是脫靶量附近相對(duì)距離隨時(shí)間的曲線近似拋物線這在制導(dǎo)末端是成立的。最后說(shuō)一個(gè)習(xí)慣我在主循環(huán)里保留一個(gè)model_verify開(kāi)關(guān)默認(rèn)關(guān)閉打開(kāi)時(shí)自動(dòng)執(zhí)行冒煙測(cè)試和降維驗(yàn)證。每次改完目標(biāo)機(jī)動(dòng)模型或者調(diào)整狀態(tài)量后先跑一遍開(kāi)關(guān)再開(kāi)始正式仿真。這個(gè)開(kāi)關(guān)相當(dāng)于給自己留了后悔藥不然參數(shù)調(diào)多了之后你很難判斷某個(gè)結(jié)果的異常是模型改錯(cuò)還是參數(shù)本身導(dǎo)致的。希望這套驗(yàn)證思路能幫你在自己的仿真建模里少走一段彎路。本文還有配套的精品資源點(diǎn)擊獲取