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