トップページへ 弾道ミサイル迎撃 ミサイルの弾道 弾道ミサイルの発射 ご意見ください
弾道ミサイルの迎撃段階には、ブースト・フェーズ、ミッド・コース及びターミナル・フェーズ（終末段階）がありますが、このホームページ「弾道ミサイルの特性」は、首都圏にペトリオット（ＰＡＣ－３）が配備されたことをきっかけに書き始めましたので、迎撃としては、ターミナル・フェーズにおける迎撃から説明します。
１９９１年の湾岸戦争におけるスカッド（Ｓｃｕｄ）改良型「アル・フセイン」弾道ミサルイル攻撃に対するペトリオット（ＰＡＣ－２）防空ミサイルによる迎撃戦闘は、当時、ＣＮＮテレビニュースをとおして報道され視聴者を釘づけにしました。しかし、攻・防ともに確たる技術をもって戦闘に臨んだわけでなく、心理的効果はともくとして、ともに実質的な戦果をあげたわけではありません。
ＰＡＣ－３システムの開発は、湾岸戦争の前（１９８３年）に着手されていましたが、湾岸戦争をとおして加速されたのは間違いないでしょう。この章では、ターミナル・フェーズにおいて、大気抗力が弾道に及ぼす影響、迎撃ミサイルの空力操舵に及ぼす影響等について説明し、迎撃追尾においてヘッドオン（Head-on）迎撃が最良なことを示します。
ミスデスタンスを少なくするためには、迎撃追尾の最終段階に迎撃ミサイルにコマンドとして 与える加速度をゼロにすることが必要です。そのためには、多くのミサイルに広く使われている「比例航法」（Proportional Navigation）に、さらに目標（Target）に関する新たな情報を追加する必要があります。これらの「Advance Guidance law」についても、説明を加えたいと思います。
ペトリオットの迎撃空間は、高々数１０Km の範囲ですから、第１章で説明した「Flat-Earth Model」が適用されまが、大気の存在しますから、それによる抗力を考慮しますと下図のようになります。
図３－１ ターミナル・フェーズにおける弾道ミサイルの座標空間
上図において弾道ミサルには、抗力 Fdrag と重力の加速度 ｇ のみが働きます。再突入（re-entory）時の目標の初期速度お及び角度をそれぞれ、Vt 及びγt とします。なお、抗力 Fdrag は速度 Vt に対して逆向きに働きます。抗力Fdrag が重力の加速度 ｇ より大きい場合は、目標は減速（decelerate）されます。
再突入角γt は、速度 Vt の（ｘ、ｙ）成分をそれぞれ、Vtxt 、Vty とすると次式で与えられます。
γt = tan-1（ -Vty / Vtx ） (3.1)
次に、目標に働く抗力を考えてみます。目標の重量をW、基準面積をS、抗力係数を CD0 とし、さらに、空気流の動圧（dynamic presure）をｑ とすると、目標には、次のような加速度が生じます。
dVtx/dt = -(Fdrag/m) cos(γt) = -(qSCD0g/W) cos(γt) = -(qg/β）cos(γt) (3.2)
dVty/dt = (Fdrag/m) sin(γt) - g = (qSCD0g/W) sin(γt) - g = (qg/β）sin(γt) - g (3.3)
q = ρ・Vｔ2/ 2 ; （ρ：大気密度） (3.4)
β = W / (SCD0) ; (弾道係数 ballistic coefficient) (3.5)
Vｔ = sqrt( Vtx2 + Vty2 ) (3.6)
上の式から、弾道係数が大きいほど大気抗力による減速が少なくなることがわかります。大気密度ρ は下図に示すように、高度が増大するにつれて減少し、対数グラフで表示すると、高度９千メートル付近を境にして２本の直線で近似できます。
図３－１ 高度対大気密度
ρ = 1.225・exp( -0.1051・h ) [kg/m3] ; 0 ≦ 高度 h [km] ＜10 [Km] (3.7)
ρ = 2.020・exp( -0.1558・h ) [kg/m3] ; 10 ≦ 高度 h [km] (3.8)
式（3.2）のところで、CD0を単に、抗力係数としましたが、正しくは「零揚力抗力係数」（又は「零離軸角抗力係数」）と呼ばれる無次元の係数です。弾頭の中心軸に空気の流れが平行にあたったときの抗力です。CD0は弾頭の形状に等に起因するもので、ノーズ（ヘッド）による造波（圧力）抗力分、弾底部抗力分及び表面摩擦抗力分などを加算したもです。超音速領域における砲弾などの例では、CD0の値は、概略、０.１～０.２です。ミサイルの弾頭は、砲弾ほど精密加工されていないと考えられるので、ここでは、０.１～０.4 としておきます。
弾道係数（β = W / (SCD0)）を「ノドン」を例にとって、見積もると概略次のようになります。
W = 1,000 [Kg] ; S = π・ｒ2 = π・（0.66)2 = 1.37 [m2] ; CD0 = ０.１～０.4
β ≒ 2,000 ～8,000 [Kg/m2] ; （端数切り上げ）
弾道係数は、弾頭（又は弾丸）の断面積当たりの重量（断面加重）ですから、この値が大きいということは、重量が大きく、断面積が小さいわけですから、大気抗力が小さくなります。下図に、弾道ミサイルの速度を一定に保って、弾道係数を変化させた場合の仮想的な減速度[G] 示します。実際には、ミサイルの速度は一定でなく、抗力にによって減少していきます。
図３－２－１ 大気抗力による減速度（パラメータ：弾道係数）
次の図は、弾道係数を一定（2,500[Kg/m2]）にして、弾道ミサイルの速度を 1,000[m/s]、2,000[m/s] 及び 3.000[m/s] に変化させた場合の仮想的な減速度[G] 示します。式（3.4） に示したように、動圧がミサイルの速度の２乗に比例しますので速度が大きいほど抗力も大きくなります。
図３－２－２ 大気抗力による減速度（パラメータ：速度）
前節で、大気密度による抗力を個々にみてきましたが、ここでは式（3.2）及び（3.3）に示した抗力についての程式を弾道ミサイルの弾道方程式に組み入れて総合的に観察することにします。使用したプログラムについては下記をクリックしてみてください。C言語によるプログラムですが拡張子は「.txt」にしてあります。
使用したプログラム： 「terminal_phase_simulation.cpp」
最初に、高度上どのあたりから大気の影響が現れるかを見るために、取敢えず高度１００ [Km] 以下を考えます。第２章の図２－１５に示したデータから高度１００ [Km] における弾道ミサイルの速度は 3,100 [m/s] ですので、再突入速度としては、この値を用い、再突入角度として、４５°、抗力の影響をみるため弾道係数を小さく見積もって、2,500 [Kg/m2] とします。結果を下図に示します。
図３－３ 大気抗力による減速度および速度
図に示すように、弾道ミサイルの速度は、高度５０ [Km] 付近までは、重力の加速度が優勢で、徐々に速度を増してきますが、５０ [Km] 以下に降下してくると、抗力の影響が大きくなって速度は急激に減少してきます。一方加速度の方は、大気密度の増加にともなって、急激に増大しますが、速度が減少するため、高度１１ [Km] 付近で、最大値（約２５[G]） を示し、以後急速に減少します。以上のことから、この節では、高度５０ [Km] 以下について各種のシミュレーション結果を示していくことにします。
シミュレーションのためのデータ設定： 再突入高度を５０ [Km] として、下表のように設定することにします。
|
項 目 |
記号 |
データ値 |
| 弾道ミサイルの速度 | Vt | 2,500 3,000 3,500 [m/s] |
| 同上 弾道係数 | β | 2,000 5,000 8,000 [Kg/m2] |
| 同上 再突入角 | γ | 30 45 60 [deg] |
表３－１ シミュレーションのためのデータ設定
（１） ターミナルフェーズにおける弾道の特性
＜大気による抗力を無視した場合＞
最初に、比較のため大気の影響を見無視した場合の「ノドン」クラスの弾道側視図及び加速度・速度図を示します。
図3-4-1
空気抗力を無視していますから、弾頭の運動エネルギーは保存され、Downrange に垂直に働く重力加速度によって、垂直方向にのみ落ち込み、弾着点で水平方向に約 ２．５[Km] （=gTf2/2=4.9×(22.42)2=2,463[m]）縮みます。弾道は、ほぼ直線とみなせます。つぎに、加速度・速度を示します。
図3-4-2
飛翔間、常に１ｇ の重力の加速度が働きますから、経過時間 Tf に比例して速度が増大します。速度の増分 ΔV は、近似的にに突入角が弾着まで維持さているとみなせば １５５[m] です。計算は次のとおりです。
ΔV = ｇ・Tf・cos(γ)=9.8×22.42×0.707=155 [m]
＜大気による抗力を考慮した場合＞
次に、再突入速度 Vt=３,０００[M/s]、再突入角度 γ=４５[deg]、弾道係数 β=２,０００[Kg/m2] の弾道ミサイルについて、減速度・速度図、弾道側視図及びDownrange に対する速度・飛翔時間図を示します。
図3-5-1
速度 ３,０００[m/s] で進入した目標（target）は、高度 ４５[Km] 辺りから速度を減じ、減速度がピークになる高度１３～１４[Km] で、１,８２０[m/s] になります。この値は、初期速度 ３,０００[m/s] に対して６０．６％になっています。最大減速度点において、速度が約 ６０% に減じるという関係は（正確には、exp(-0.5)=0.6065）は、理論的に導かれることでいつでも成り立ちます。次に、弾道側視図を示します。
図3-5-2
弾道側視図の著しい特徴は、弾道の直線性です。飛翔間、大きな減速度により急激な速度の低下を受けても弾道の直線性は維持されます。これは、空気抗力の向きが弾頭の中心線（弾軸）に平行であり、その抗力の突入角 γ の余弦分（cos(γ)）及び正弦分（sin(γ)）が水平成分及び垂直成分に正しく対応しているからです。弾頭の存速が低下して抗力が重力の加速度以下になると直線性は失われて下方に湾曲します。上の図でも、Downrange ４２[Km] 以降にこの現象が出ています。
大気を無視した場合の弾道側視図（図３－４－１）においても直線性は維持されていますが、この場合は、速度を一定の割合でわずかに上げながら、直線に沿って突入してきます。飛翔時間は ２２．４２ 秒と、非常に短いです。一方、大気中では、速度を急激に失いながら直線に沿って突入してきます。飛翔時間は ４７．１５ 秒です。この値は、大気を無視した場合の ２ 倍強です。突入時 ３,０００[m/s] あった速度は、最終的に ２２１[m/s] まで減少します。 Downrange に対する平均対地速度は、前者が ２,１８２[m/s]、後者が９６６[m/s] です。次に、Downrange に対する飛翔速度及び時間を示します。
図3-5-3
ターミナル・フェーズにおける弾道ミサイル迎撃においては、一般に、突入時の目標（target）速度が迎撃ミサイルの速度より大きいでしょう。そこで、追尾可能な空域に迎撃点（intercept point）を選定し、そこに至る直前までは、指令誘導等でミサイルを誘導し、迎撃可能な最終段階でミサイル自体にホーミング誘導させて迎撃することになります。このとき、飛翔時間 Tf のうちの残り時間が、初期誤差を詰めきるだけの余裕あることが必要です。図に示すように、飛翔時間の半分は Downrange 上 ３９ [km] 以降に残されています。
（２） パラメータ分析（速度、弾道係数、突入角が減速度に及ぼす影響）
ここでは、弾道ミサイルの速度 Vt、弾道係数 β 及びを突入角 γ を表３－１に示した値に設定して、弾道ミサイルに生ずる減速度を見ていくことにします。
＜速度 Vt：2,500、3,000、3,500 [m/s]＞
図3-6-1
図3-6-2
図3-6-3ｇ
弾道ミサイルの速度が受ける動圧は、速度 Vt の二乗に比例しますから、減速度も速度の二乗に比例して増大しますが、最大動圧が生起する高度は変わりません。この例では、高度１８[Km] で、最大減速度が生じています。
＜弾道係数 β：2,000、5,000、8,000 [Kg/m2]＞
図3-7-1
図3-7-2
図3-7-3
弾道係数βは、大気抗力を受ける断面当積（この面積は見かけ上、抗力係数に比例する）当たりの荷重ですから、この値が大きいほど、弾頭の存速は維持され、最大減速度が生ずる位置は低高度まで持ち越されます。ただし、最大減速度の値は変わりません。
＜突入角 γ：30、45、60 [deg]＞
図3-8-1
図3-8-2
図3-8-3
突入角 γ が大きいほど、大気を通過する経路長が短くなり（飛翔時間も短くなる）ので、最大減速度が生ずる高度は低くなります。ただし、最大減速度の値は、ｓｉｎ（γ） に比例しますので、γ が大きくなるにつれて増大します。
なお、最大減速度 Amax は、突入速度 Vt 及び突入角 γ にのみ依存して、次式で与えられます。
Amax = 3.30 × 10-6 × Vt2 × sin（γ） [G] (3.10)
また、最大減速度 Amax が生起する高度 AltAmax は弾道係数 β 及び突入角 γ にのみ依存して、次式で与えられます（ln ： 自然対数）。
AltAmax = 5.513 × ln ( 17450 / ( β × sin(γ) )) [km] (3.11)
式（3.10）及び(3.11)をそれぞれグラフ化して、次に示しします。
図3-9-1
図3-9-2
レーダ装置によって再突入する弾道ミサイルを捕捉・追尾することによって、突入角 γ と突入速度 Vt を知ることができるますから、式（3.10) を使って最大加速度を予測できます。しかし、弾道係数 β を知ることはできませんから、残念ながら、最大減速度が生起する高度を予測できません。予想邀撃点をできる限り正確に推定することが弾道ミサイルの迎撃においてキーポイントになりますので、弾道係数を推測する方法については後ほど説明を追加します。
ここでは、ペトリオットなど、ターミナル・フェーズにおける迎撃用ミサイルを単純化して円柱とみなして、大気中を超音速で飛翔する場合の発生揚力について説明します。円柱モデルの揚力係数 CL（Lift Coefficient）は、迎角（Atack Angle）をα[rad]、基準断面積を Sref[m2] 及び円柱の平面投射面積を Splan[m2] とすると次式で与えられます。
CL = 2α + 1.5α2 Splan / Sref (3.12)
また、円柱の長さを Ｌ[m] 、直径を Ｄ[m] とすると、Splan 及び Sref は次のようになります。
Splan = L・Ｄ [m2] (3.13)
Sref = π・Ｄ２ / 4 [m2] (3.14)
次に、空力操舵によりミサイルに発生する横加速度を nL [m/s2]、ミサイルの重量を W [kg] として、揚力係数を含めて、運動の法則を考えると字式が得られます。
F = m・a = W・nL / g = q・Sref・CL [N] （3.15）
上式で、動圧は式（3.4）で示したように、 q = ρVm2 / 2 ですから、ミサイルに発生する横加速度 nL は次のようになります。
nL / g = q・Sref・CL / w = ρVm2 ・Sref・CL / （2・w） [G] （3.16）
定式に、代表的なミサイルの寸法、重量として、Ｌ＝５ [m]、Ｄ＝０．２５ [m] 及び Ｗ＝３００ [Kg] を代入して、迎角α 及び ミサイル速度 Vm をパラメータとしてグラフ化した図を次に示します。
図3-10-1
図3-10-2
いま、弾道ミサイルの再突入角度及び速度を45 [deg] 及び 3,000 [m/s] （図3-8-2参照）とし、一方、迎撃ミサイルの速度を1,000[m/s]、迎角を２０[deg] （図3-10-2参照）とすると、高度約１0[Km] より高い空域では、空力操舵のみでは、迎撃ミサイルは目標（target）を追尾しきれないことがわかります。高空域で迎撃するには、サイド・スラスター等の空力操舵以外の手段必要です。この関係図を「図3-8-2」及び「図3-10-2」から必要な部分を抜粋してプロットしなおした図を次に示します。
図3-10-3
また、この図から迎撃ミサイルの速度が１５００[m/s]（弾道ミサイルの１／２）であれば、高度１5[Km] まで戦えることも分かります。
この節では、これまで説明してきたターミナル・フェーズにける弾道ミサイルの特性を踏まえて、「ノドン」クラスのミサイルを想定してコンピュータ上で迎撃シミュレーションを行い、ターミナル・フェーズの交戦のキーポイントを見ていくことにします。
① 交戦開始（迎撃ミサイル発射）時の目標高度 ： ６０ [Km]
② 交戦開始時の目標速度（目標初期速度） ： ３,０００ [m/s]
③ 目標の再突入角度 ； ４５ [deg]
④ 目標の弾道係数 ： ２,０００ [Kg/m2]
ターミナル・フェーズにおける目標の「弾道側視図」、「減速度・速度」及び「経過秒時」を以下に示します。
① 弾道側視図
図3-11-1
② 減速度及び速度
図3-11-2
③ 経過秒時
図3-11-3
ターミナル・フェーズにおける弾道ミサイルの軌道は図3-11-1 に示したように、再突入角を維持して、ほぼ直線ですから、この直線のどこで迎撃するかは、防御側が決定できます。迎撃高度を決定する要因は、時間的な要因と迎撃ミサイルの追尾能力（速度と加速度）によります。時間的要因は、迎撃点で目標と会合できるか否かです。追尾能力は、迎撃点に至るまでにミス・デスタンスをゼロにできるかです。これは、迎撃ミサイルの速度・加速度特性と目標の減速度・速度特性に依存します。
ここでは、迎撃点高度を次の２点を設定してシミュレーションすることにします。
| 迎撃点高度 [Km] | 経過秒時 [s] | 速度 [m/s] | 減速度 [G] |
| 15 | 21.69 | 2,065 | 20.5 |
| 10 | 25.59 | 1,285 | 17.0 |
① 迎撃ミサイル発射前に、所望迎撃高度及び前述の（１）に示した迎撃対象目標に関する初期条件に基づいて、迎撃点における目標位置座標及び経過秒時を算定します。C言語プログラム上では、関数「intercept_point( )」を使って計算しています。
② 迎撃ミサイルの目標追尾の方式（航法）としては、比例航法（proportional navigation）を用いていますが、参考のために、増強比例航法（augmented proportional navigation）を使用できるようになっています。増強比例航法は、通常の比例航法に目標に関する加速度情報を追加ことによって、迎撃ミサイルに要求させる横加速度を緩和させる効果があります。ただし、現実的には目標に関する新たな情報を獲得し、それを比例航法ループの中に組み入れなけばならないという問題があります。
③ 迎撃ミサイルの配備位置（発射機の座標位置）については、当初、目標の弾着予想位置（図3-11-1「弾道側側視図」上で、原点か ５４．５５ [Km] ）に配置し、その後、弾着予想点から ５ [Km] 単位で移動させてシミュレーションしてみます。
④ プログラムから出力されるデータは、経過秒時並びに経過秒時に応ずる目標座標位置、迎撃ミサイル座標位置・速度及び視準線（Line Of Sight; LOS, 迎撃ミサイルから目標を見通した線）に直交する目標加速度成分です。
⑤ 使用したプログラムについては、次をクリックして取得してください。（拡張子は「,txt」になっています）
プログラム「 engagement_simulation.cpp 」
＜所望迎撃高度：１５[km],、迎撃ミサイルを予想ｆ弾着点に配置さした場合＞
最初に、迎撃ミサイルを、予想弾着点という特別な位置に配置してシミュレーションを行ってみます。ターミナル・フェーズにおける弾道ミサイルの軌道はほぼ直線になり、航空機のように回避運動はありませんから、予想弾着点から発射された迎撃ミサイルは、特別に努力しなくても命中することが期待されます。
① 弾道ミサイルと迎撃ミサイルの弾道側視図
図3-12-1
上図において、赤色の実線で示したグラフが弾道ミサイルの弾道です。また、点線で示したグラフは迎撃ミサルの弾道です。予想弾着点（座標原点から約５５[Km]の位置）から発射され所望の迎撃点（高度１５ [Km] で目標に会合しています。
迎撃はほぼヘッドオン（head -on）に近い状態で生じています。点線で示した迎撃ミサイルの弾道は、弾道ミサイルの予想弾道を逆にたどった軌跡になっていることから「逆軌道（inverse trajectory）」と呼ばれています。
② 迎撃ミサイルに生起する横加速度コマンド（誘導コマンド）
図3-12-２
迎撃ミサイルが目標を追尾する基本則である比例航法の原理は、「目標の加速度→視準線（ＬＯＳ）の回転（角速度）→迎撃ミサイルの誘導コマンド（横加速度）」です。上図に示したように、目標加速度（減速度）は迎撃点付近で２０ [G] を越えていますが、弾道の直線性が維持されています。一方、視準線（ＬＯＳ）が弾道と一致しているため、ＬＯＳの回転は少なく、結果として、迎撃ミサイルは、ほぼ直線的に迎撃点に向かって命中しています。
＜迎撃ミサイルの配備位置がを目標（Target）の弾着点から離れている場合＞
迎撃ミサルの配備位置が目標の予想弾着点から離れるに従って視準線（ＬＯＳ）の回転速度が大きくなりますから、迎撃ミサイルに生起する誘導コマンド加速度は大きくなります。シミュレーション結果の一例として、迎撃ミサイルが目標の予想弾着点か２５[Km] 手前（座標原点から３０[Km]）に配備されている場合の弾道側視図及び迎撃ミサイルに生起する誘導コマンド加速度を次に示します。
図3-13-1
図3-13-2
迎撃機ミサイルが、木曜の予想弾着点に配備されている場合の迎撃ミサイル誘導コマンド度加速度の図3-12-２と図3-13-2を見比べてください。迎撃直前の前者の誘導コマンド度加速度は １ [G] ほどですが、後者では、１５ [G] 増加しています。シミュレーション結果から、配備位置が予想弾着点から、手前に１ [Km] ずれるとごとに約 ０.６ [G] 増加します。
このシミュレーションでは、迎撃ミサイルシステムを「ゼロ遅れ」システムとして扱っているため、機体は誘導コマンドに遅れなく応答していますが、実際には２次遅れシステムと考えられますから、迎撃の最終段階で誘導コマンド加速度が残ることは好ましくありません。
ミサイルの誘導則である比例航法に目標（Target）の加速度ｊ情報を追加することで、この問題の解決を図った「増強比例航法（augmented proportional navigation）があります。次に、そのシミュレーション結果を示します。
＜増強比例航法（augmented proportional navigation）を適用した場合＞
比較のため、目標（target）に関するデータは、前の図と同じです。
図3-14
図に示したように、迎撃の最終段階で誘導コマンド加速度がゼロに収れんしていくのがわかります。
＜所望迎撃高度：１０[Km],、迎撃ミサイル配備位置：座標原点か３０[Km]＞
これまでのシミュレーション例では、所望迎撃高度を １５ [Km] に設定していました。この高度は、図3-11-2 に示したように目標（Targrt）の減速度がピークにさしかかる高度で、速度は２，０００[m/s] に減速されています。これに対して、次に示すシミュレーションは、所望迎撃高度を １０ [Km] に設定した例です。この高度では、減速度がピークを過ぎて、速度がさらに １，３００ [m/s] に減速されています。
図3-15-1
図3-15-2
所望迎撃高度を １０ [Km] に設定した以外は、すべて、図3-13-1、-2 と同じ設定です。迎撃の最終段階で誘導コマンド加速度が減少に転じ、前者で１５[G] あった減速度は ８ [G] になっていることが分かります。
① ターミナル・フェーズにおける弾道ミサイル（目標）の軌道はほぼ、再突入角を維持した直線となる。この結果、防者の欲する空域（高度）で、迎撃できる。
② 目標は直線的な軌道を維持した状態で減速され、減速度の最大点を経て減速度は減少する。
・ 減速度の大きさは、目標の速度と再突入角度（γ）に依存し、それらの値がが大きいほど大きくなる。
・ 減速度が最大になる高度は、目標の弾道係数（β）と再突入角度（γ）に依存し、それぞれ、値が大きいきほど、低くなる。
・ 他の条件が許すなら、迎撃高度を最大減速度点以降に設定したほうが迎撃ミサイルの負担が少なくなる。
③ 迎撃ミサイルの配備位置が目標の予想弾着点に近いほどミス・デスタンスは小さくなる。
④ 迎撃ミサイルの配備位置が目標の予想弾着点から離れると、迎撃の最終段階に過度の誘導コマンド加速度が生ずるので、増強比例航法（augmented proportional navigation）等、目標運動に関する情報を取り込み、誘導コマンド加速度の緩和を図る必要がある。
⑤ 再突入目標に関する情報のうち、速度、再突入角度については、レーダ装置により容易に取得できが、迎撃高度を選定し、その迎撃点にミサイルを指向するためには、目標の弾道係数が必要である。
弾道係数は、超音速風洞試験等によって概略の値は推定できますが、ここでは、迎撃時にレーダ装置によって目標を補足追尾する過程で推定する拡張カルマンフィルタについて説明します。
PATRIOTの参考マニュアル：
「FM 3-01.85 PATRIOT BATTALION AND BATTERY OPERATION」
（１） フィルタ
フィルタ（filter）という語は日常の用語として例えば、クーラーの「フィルター」のように使われています。「濾（こ）す」という意味で、その器具は「ろ（濾）過器」ということになります。電気信号の場合は、特定の周波数の信号を通過させるという意味で「ろ（濾）波器」という言葉で使われます。
特定の周波数の電気信号とは、信号として意味のある信号成分で、この信号成分に雑音成分が加わっている場合、信号成分だけを取り出すことを「フィルタリング（filtering）」と呼びます。
「フィルタ」という用語は、いままでの説明で、フィルタリングする器具や装置（ハードウェア）のように取られますが、ハードウェアのみならず、アルゴリズム（ソフトウェア）をも含んでいます。これから説明する「カルマンフィルタ」もアルゴリズム（算法）です。
私たちが側測定（観測、受信）できる値（データ、情報など）は、何が信号成分（あるいは真の値）で何が雑音成分なのか分かりません。過去から現在までのデータを適当に処理して真の値として推定して用いていざるを得ません。
一般に、現時点に立って、過去のデータをいかなる視点から処理するかによって、次のように使い分けられています。
① 「最適予測」： 未来の最適な推定量を求めること。（prediction,）
② 「最適フィルタリング」： 現在の最適な推定値をもとめること。（filtering）
③ 「最適スムージング」： 過去の適当な推定量を求めること。（smoothing）
なお、①の問題は、取り扱う対象が線形の差分式で表される場合は、②の問題に帰着されます。
（２） カルマン・フィルタ
カルマン（Rudolf Emil Kalman, 米, 1930 年～、スタンフォード大学教授などを経てスイス連邦工科大学教授）は、１９６０年に発表した「On The General Theory of Control System」及び翌１９６１年、数学者の R.S. Bucy との連名で発表した「New Results on Linear Filtering and Prediction Theory」 により、制御工学に時代を画する根本的な変革をもたらしました。それまで、制御対象システムを伝達関数で表していた制御理論は、これを期にいわゆる「状態変数」で表現されるようになって、「現代制御理論」と呼ばれ、それ以前の制御理論は「古典制御理論」と呼ばれるようになりました。なお、カルマンは、第１回（1985年）京都賞を受賞しています。
＜問題の定式化＞
最適フィルタについての形式的な定式化は、カルマンフィルタ以前に、ノバート・ウィナー（N.Wiener,米,1894～1964）らによって、積分方程式（「Wiener-Hoph の方程式」）の形で与えられていた。しかし、残念ながら、この方程式を解く一般的な方法は見つかりませんでした。
ウィナーのフィルタリングの問題は、信号や雑音のスペクトル密度を既知のものとして、これらの統計量に基づいて、前述の積分方程式を解き、最適フィルタの伝達関数を求めるものでした。
これに対して、カルマンは、スペクトル密度関数の代わりに、信号や雑音の時系列を実時間で捉え、信号は雑音（白色雑音）によって駆動された「線形ダイナミカルシステム」の出力としました。
カルマンは、最初に「線形ダイナミカルシステム」から生成される信号を考え、一般的に、信号過程が次式で表されるものとして、フィルタリングの問題を定式化しました。
ｄｘ（ｔ）/ ｄｔ ＝ Ｆ（ｔ）ｘ（ｔ） ＋ Ｇ（ｔ）ｕ（ｔ） ＋ ｗ（ｔ） （3.17）
ここで、ｘ（ｔ）；状態変数列ベクトル、Ｆ（ｔ）；システムダイナミクス行列、ｕ（ｔ）：制御ベクトル、ｗ（ｔ）：白色雑音（平均値=0,分散=E[w wt])
また、状態変数の観測値は、雑音 ｖ（ｔ）（平均値=0,分散=E[v vt]) が存在する測定機構
ｚ（ｔ） ＝ Ｈｘ（ｔ） ＋ ｖ（ｔ） （3.18）
を通して得られるものと仮定して、過去のある時点から現時点までの測定データを処理して得られるｘ（ｔ）の推定値のうちで、誤差の二乗平均値が最小にする問題としました。
式（3.17）の Ｆ を ｎ×ｎ 実行列とするとき、すべての実数 ｔ に対して、べき級数
Ｉ ＋ Ｆｔ ＋ （１/２！）Ｆ2ｔ2 ＋ ・・・ （3.19）
を、Φ（ｔ） と書くことにします。そうすると、Ｆ を一定値とする微分方程式（3.17）のラプラス変換に基づく次の表現が得られます。Ｉ は単位行列、Ｌ-1 はラプラス逆変換、ｓ はラプラス演算子です。
Φ（ｔ） = Ｌ-1｛（ｓＩ - F）-1｝ （3.20）
Φ（ｔ）は Ｆ から導かれた状態遷移行列（state transition matrix）と呼ばれています。
離散時間シスでは、サンプリング時間を Ｔｓ とし、現時点（ｔ）を第 ｋ 番目とすると次のように表現できます。
Φk = Φ（Ｔs） （3.21）
このΦk を用いて、離散時間のカルマンフィルタの解（状態変数の推定値X^）次式で与えられます。
X^k = ΦkX^k-1 + Gkuk-1 + Kk(zk - HΦkX^k-1 - HGkuk-1) （3.22）
上式の Kk はカルマン・ゲインマトリクス、Gkは次式で与えられます。
Gk= integral(0,Ts)Φ（τ）Ｇ dτ ； integral( ) は∫記号 （3.23）
拡張カルマンフィルタ（Extended Kalman filter） では、式（3.17）及び式（3.18）は非線形な方程式ですので、システムダイナミック行列 Ｆ 及び測定（観測）行列 Ｈ については、リカッチ（Ricati）の方程式（注）の一次近似が使われます。 また、状態遷移行列いても、リカッチの方程式を満足することが要求されますが、これについても、通常、一次近似として、テイラー級数展開（Taylor series expansion）の第２項までとってj次式のように近似しています。
Φk ≒ Ｉ ＋ Ｆ・Ｔs （3.24）
（注） Ricati の微分方程式は、１８世紀イタリアの Ricati 伯爵（Jacopo Francesco, Conte Ricati, 1676-1754）による次の形をした微分方程式です。ｙ について、２次項を含んでいるのが特徴です。
ｄｙ/ｄｘ ＋ Ｐ（ｘ）ｙ ＋ Ｑ（ｘ）ｙ2 = Ｒ（ｘ）
式（3.22）に示したように、最適フィルタを決定するためには、カルマン・フィルタゲイン Ｋk の値を決める必要があります。そして、このゲインは、次に示すリカッチの行列方程式を再帰的に解くことよって得られます。行列の右上につけた添字（T）は転置行列です。
図3-16
上図で、① の信号処理雑音行列 Qk 次式で与えられます。
Ｑk = integral( 0, Ts ) Φ(τ)ＱΦTdτ （3.25）
カルマンゲイン Ｋk が決まると、状態変数の推定値 ｘ^k+1は、観測値 Zk をもとに次式から求められます。
ｘ^k+1 = ｘ＾k + Ｋk [ ｚk - ｈ（ｘ＾k）] （3.26）
（３） １次元に単純化した弾道ミサイルの微分方程式
ここでは、弾道ミサイルの弾道係数を推定するためのカルマンフィルタを検討することが目的です。弾道係数（β）は大気密度（ρ）に応じて、弾道ミサイルに対する抗力（Ｆdrag）として働きます。そして、大気密度（ρ） は高度に依存します。したがって、弾道係数を推定するための最適フィルタを検討するためには、弾道ミサイルの垂直方向(ｙ軸)の運動のみに着目すれば十分です。
弾道ミサイルが、地上に設置したレーダに向かって、垂直に再突入してきた場合の弾道ミサイルの運動方程式は次式で与えられます。
dVt / dt = （Fdrag / m） - ｇ = （ｑ・Ｓref・Ｃd0・ｇ/ Ｗ）- g = （ｑ・ｇ / β） - ｇ （3.27）
ｑ = ρ・Ｖt2 / 2 （3.28）
ρ = ２．２０２・exp（-0.155yｔ/1000） / g （3.29）
ｄｙ / ｄｔ = Ｖt （3.30）
ｄβ / dt = 0 （3.31）
（４） １次元弾道ミサイル軌道に対する拡張カルマンフィルタ
＜状態変数＞
状態変数は、目標（target 弾道ミサイル）の高度（Ｙt）、速度（ｙｔ）及び弾道係数（β）ですから、次式のような列（縦）ベクトルで表されます。
ｘ ＝ [ Ｙt Ｖt β ] ’ ； （[ ] ’は列ベクトルを示す。） （3.32）
ここで、目標位置（高度）（Ｙt）はレーダによって直接測定されます。測定雑音（ｖk）を含んだ測定値（Ｙt*）は状態方程式で、次のように表現されます。 average( ) は平均値を示します。
Ｙt*＝ Ｙt + ｖk = [ １ ０ ０ ] [Ｙt Ｖt β ] ’ + ｖk （3.32）
＜状態変数の推定式＞ 各文字の右肩の添字（＾; hat）は推定値を示す。
Ｙt^k ＝ Ｙt^k-1 + average（Ｙt）Ｔs + K1 * Residual （3.33）
Ｖt^k ＝ Ｖt^k-1 + average（Ｖt）Ｔs + Ｋ2 * Residual （3.34）
β^k ＝ β^k-1 + Ｋ3 * Residual （3.35）
ここで、
Residual ＝ Ｙt* - Ｙt^k-1 - （Ｖt^k-1）Ｔs （3.36）
average（Ｙt） ＝ Ｖt^k-1 （3.37）
average（Ｖt） ＝ ρｇＶt^k-12 / （2β^k-1） （3.37）
（５） 弾道係数推定プログラム
＜使用したプログラム＞
使用したプログラムについては、次をクリックして取得してください。（拡張子は「,txt」になっています）
プログラム「ballistic_coefficient_estimation.cpp」
＜初期設定値＞
| 項 目 | 真 値 | 真値 ＋ 誤差 | 誤差の分散 |
| 初期高度 | 30,000 [m] | 30,010 [m] | 10 × 10[m2] |
| 初期速度 | 3,000 [m/s] | 3,100 [m/s] | 100 × 100 [m2/s2] |
| 弾道係数 | 2,000 [kg/m2] | 3,000 [kg/m2] |
1,000 × 1,000 [kg2/m4] |
信号処理雑音（Φs）（process noise）については、当初、ゼロに設定しておき、シミュレーションの過程で適宜に加えることにします。
数値計算の積分時間刻み幅は、０．００１秒（１０[ms]）、サンプリング 時間（τ）は０．０５秒（５０[ms]）です。
データ出力間隔は、Ｅｘｃｅｌでグラフを描画する都合上、１２５[ms]（１秒間に８回, τ の２.5倍）に設定しています。
（６） シミュレーション結果
このプログラムをランさせた結果は、２つのテキストファイルに分けて次に示す内容が出力されます。
① 「DATFIL.TXT」：
経過秒時、高度、推定高度、速度、推定速度、弾道係数、推定弾道係数
② 「COVFIL.TXT」：
経過秒時、高度、高度推定誤差、同理論上の標準偏差（以下「±σ」）、速度推定誤差、同±σ、弾道係数推定誤差、同±σ
ここでは、出力結果のおうち、「弾道係数(真値）と推定弾道係数」及び「弾道係数推定誤差と同±σ」を以下にグラフ化して示すことにします。
＜弾道係数(真値）＝２,０００[kg/m2] に対して、推定弾道係数（β^）＝３０００[kg/m2]とした場合＞
図3-17-1
図3-17-2
上図から、当初１,０００[kg/m2] あった推定誤差が、高度が減少するにつれて真値（２,０００[kg/m2] に収斂して誤差がゼロになる様子が理解できるでしょう。
＜測定開始の初期高度（Alt）を ６０[km] に上げた場合＞
前の図では、高度３０ [km] から測定を開始しましたが、測定開始高度を高くした方がだ弾道係数推定上有利なのでしょうか。次の図では、高度を３０ [km] から２倍の６０ [km] にした場合の例を示します。
図3-18-1
図3-18-2
上図に示したように、高度６０ [km] から測定を開始した場合でも、弾道係数の推定値として使える程度に収束するには高度１５ [km] 以下になるまで計測を継続する必要があります。このことから、図3-17 と比べて、高度３０ [km] から測定を開始しても十分なことがわかります。
弾道係数（β）は、大気の密度を（ρ）、弾道ミサイルの速度を（V）とすると、「ρ・V2/（２β）」 という形（減速度）でミサイルの軌道に影響を及ぼしますが、この影響が現れ始めるのは高度 ３０[km] 付近からです。これより高い高度では、対気密度（ρ）が小さく、ミサイルの速度変化も少ないため、過去のデータを蓄積しても、そこから有効な推定値が得られないのです。
＜弾道係数の初期推定値を ３,０００ [kg/m2] から ６,０００ [kg/m2] に設定した場合＞
図3-19-1
図3-19-2
上図に示したように、弾道係数の初期推定値を真値（２,０００ [kg/m2] ）から３倍の６,０００ [kg/m2] に設定した場合でも、最終的に正しい値に収束しています。
＜弾道係数の初期推定値を ６,０００ [kg/m2] から ７,０００ [kg/m2] に設定した場合＞
図3-20-1
図3-20-2
上図に示したように、弾道係数の初期推定値を段々増加していき、７,０００ [kg/m2] に増加させると、最終的に誤差が 約 ３００ [kg/m2] 残ってしまいます。 すなわち、弾道係数の初期見積もりが過大になるとフィルタのロバスト性（robust；頑健性）が崩れてしまいます。
この設定で、プロセス雑音（Φs）（process noize）を追加させると、次に示すように、誤差をほぼゼロになるように収束させることができました。 プロセス雑音は式（3.17) のｗ（t）に対応する雑音（共分散）です。
＜弾道係数の初期推定値を７,０００ [kg/m2] にし、プロセス雑音（１００） を追加した場合＞
図3-21-1
図3-21-2
可能な初期推定値の変動内で最悪状態設計を行うことを、ロバスト制御の基本とするならば、ここで提示した拡張カルマンフィルタは、プロセス雑音を加えることによって実際の弾道係数（２,０００ [kg/m2]）に対して、最大７,３００ [kg/m2] の初期推定値まで、安定にフィルタ機能を果たします。
ターミナルフェーズにおける迎撃のポイントは、２つあります。その第１は、最適な迎撃高度の設定です。第２は、目標（target）に関する加速度情報を追尾航法に組み入れることです。
第１の最適な迎撃高度を設定すためには、大気抗力により目標（target）に最大の減速が生ずる高度を知る必要があります。その高度は式（3.11）で示されます。
AltAmax = 5.513 × ln ( 17450 / ( β × sin(γ) )) [km] （再掲） (3.11)
上式で、弾道係数（β）は、直接測定できませんから、カルマンフィルタのような予測プログラムをレーダ装置に組み込んで、個々の迎撃戦闘間に推定して特定の弾道ミサイルの弾頭についてデータを蓄積しておく必要があります。
第２の目標（target）に関する加速度（減速度）（At）は、大気密度を（ρ）、速度を（Vt）、重力の加速度を（ｇ）としますと次式で示されます。
Ａｔ ＝ρ・ｇ・Ｖｔ2 / （２β）
上式で、大気密度（ρ）、重力の加速度（ｇ）は既知であり、目標（target）の速度は直接レーダで算定可能ですから、弾道係数（β）の値が推定できれば、目標（target）に関する加速度情報も得られることになります。
（参考）「FM 3-01.85 PATRIOT BATTALION AND BATTERY OPERATION」
（以下第世４章へ）