砲撃体位置明確化のためのシステムおよび方法
【課題】超音速飛翔体の砲撃体を位置決めするための方法を提供する。
【解決手段】方法はアンテナを形成する５以上の間隔をあけた音響センサにより、衝撃波単独信号を計測し、かつタイミングエラー分布を推定し、単独信号から、センサ対に到着する時間差を、信頼度に基づいて決定することにより、飛翔体軌跡を得る。
【解決手段】方法はアンテナを形成する５以上の間隔をあけた音響センサにより、衝撃波単独信号を計測し、かつタイミングエラー分布を推定し、単独信号から、センサ対に到着する時間差を、信頼度に基づいて決定することにより、飛翔体軌跡を得る。
【発明の詳細な説明】
【技術分野】
【０００１】
政府契約
合衆国政府は、本発明において、支払い済みのライセンスと、限定された状況において、ＤＡＲＰＡＡＴＯにより裁定された契約番号第ＨＲ００１１−０４−Ｃ−００３５の条項によって規定される適切な条件で特許権者に他にライセンスすることを要求する権利とを有する。
発明の分野
本発明は、法的措置技術および安全に関し、より詳細には超音速飛翔体の発生源と移動の方向を決定するための方法とシステムに関する。本方向およびシステムは、砲撃体（ｓｈｏｏｔｅｒ）とセンサ間の距離が大きいときでも、砲口音からの信号が無いか微弱なときでも、砲撃体位置を決定し明確化することができる方法とシステムである。
【背景技術】
【０００２】
飛翔体によって生成された衝撃波に伴うパラメータを計測することで、弾丸や砲弾等の超音速飛翔体の概略の方向と軌跡を決定することができるシステムと方法が公知である。このようなシステムの１つは、米国特許第５，２４１，５１８号の記載によると、少なくとも３つの間隔を空けたセンサを含んでおり、各センサには１つの平面に配置された３つの音響トランスデューサが組み込まれている。センサは、その発生源に対する方位と仰角に関連する衝撃波に反応する。衝撃波単独計測では、センサと衝撃波の発生源間の距離を決定することができない。通常、距離情報は、火光や砲口爆風から得られる。
【特許文献１】米国特許第５，２４１，５１８号公報
【０００３】
センサ位置に対する砲撃体の方位と仰角は、通常は角センサにおける砲口信号と衝撃波信号の「ＴＯＡ」すなわち到着時間（Ｔｉｍｅ ｏｆ Ａｒｒｉｖａｌ）情報を計測することで決定される。各センサは、異なる時間において信号と遭遇して砲口および衝撃波圧力に応答した信号を生成する。各種のセンサからの信号が処理され、砲口および衝撃波の発生源に対するセンサからの方向（方位および仰度）とともに飛翔体の軌跡を決定することができる。
【０００４】
従来のシステムでは、マイクロフォンが採用され、これは、比較的近接したり（例えば１ｍ間隔等）、広く分散（例えば、車輌に搭載されたり、戦場の兵士が持ち歩いたり等）しており、それぞれの位置で砲口および衝撃波圧力を全方位的に計測する。しかし、センサが比較的広く間隔が空けられ、および／または軌跡がアンテナ内にあると、衝撃波単独解を正確に得るのに必要なタイミング精度が非常に高度となり、特別な技術が要求される。
【０００５】
大きなアンテナサイズは、例えば車輌搭載システムにおいては大変不利である。さらに、小さい時間分解能しか有さないシステでは、不明確な解を出し得、これでは、所与の一式のセンサにおける、２つの鏡面対称の砲撃体位置に対する衝撃波のＴＯＡ時間情報がほぼ同一となる。
【０００６】
従来のアルゴリズムは、少なくとも４つの衝撃波および砲口検出を要求するので、衝撃波ＴＯＡに関して、４×４行列が反転されて平面波がマッピング可能とする。衝撃および砲口ＴＯＡ決定における微小な誤差が、距離（ｒａｎｇｅ）推定における大きな誤差を生み出し得る。さらに、従来のアルゴリズムは、弾丸軌跡に沿って一定の弾丸速度を仮定しており、これでは、約３００ｍを超える距離から発射されている長距離弾に対しては、距離推定が不正確となる。
【０００７】
従って、長い砲撃体距離を正確に推測することを可能とする迅速に変換するアルゴリズムが必要とされている。また、砲撃体方向に対する衝撃波単独解を明確化する必要もあった。さらに、砲口爆風に関係のない音響的サインに隠され得る砲口信号を抽出する必要がある。
【発明の開示】
【０００８】
発明の要旨
本発明は、様々な実施の形態により、従来技術の欠点に対処したものであり、特に、砲口信号が弱いか検出チャンネルにおいて検出された数が不十分であるときの、長距離に対する砲撃体距離（ｓｈｏｏｔｅｒ ｒａｎｇｅ）の推定のための方法およびシステムを提供する。開示の方法およびシステムは、最適化処理中で、微弱なおよび／または信頼度の低い検出砲口信号を含むことにより、さらに改善を達成して衝撃波単独砲撃体軌跡解の明確化も改善する。
【０００９】
本発明の１つの局面によると、衝撃波単独信号から飛翔体軌跡を明確化するための方法であって、アンテナを形成する５以上の間隔を空けた音響センサにおいて、衝撃波単独信号の少なくとも初期部分を計測する工程と、前記音響センサに対するタイミングエラー分布を推定する工程と、前記衝撃波単独信号の計測された初期部分から、センサ対に対する到着の時間差（Ｔｉｍｅ−Ｄｉｆｆｅｒｅｎｃｅ−ｏｆ−Ａｒｒｉｖａｌ：ＴＤＯＡ）を、前記推定されたタイミングエラー分布より大きい時間分解能で決定する工程と、明確化に対する規定された信頼度レベルに基づいて、および前記音響センサのＴＤＯＡに対する残差の値に基づいて、前記明確化した飛翔体軌跡を選択する工程と、を含む方法である。
【００１０】
本発明の別の局面によると、衝撃波単独信号から飛翔体軌跡を明確化するための方法であって、アンテナを形成する５以上の間隔を空けた音響センサにおいて、衝撃波単独信号の少なくとも初期部分を計測する工程と、各センサ対に対して、計測された衝撃波単独信号初期部分から到着の時間差（ＴＤＯＡ）を決定する工程と、遺伝的アルゴリズムを、飛翔体軌跡想定を備える初期染色体に、所定の世代数適用する工程と、前記遺伝的アルゴリズムから前記染色体とともに得られた解に対する残差を算出する工程と、最小の残差を有する解およびその不明確代替解について傾斜サーチを実行する工程と、前記最小の残差を有する解およびその不明確代替解の比率が、所定値より大きい場合は、前記最小の残差を有する解を前記明確化した飛翔体軌跡として指定する工程とを含む。
【００１１】
さらに本発明の別の局面によると、衝撃波および砲口爆風を検出することにより、砲撃体距離を推定するための方法において、衝撃波単独信号および砲口爆風信号が計測される。初期の砲撃体距離が、初期の弾丸速度および弾丸抗力係数を想定して、計測された衝撃波および砲口爆風信号から推定され、弾丸軌跡に沿った瞬時の弾丸速度が、反復的に算出されて更新された砲撃体距離を得る。砲口爆風検出チャンネルの数は、通常は衝撃波検出チャンネルの数よりも少ない。
【００１２】
有利な実施の形態は、１以上の次の特徴を含むことができる。衝撃波単独信号および砲口爆風信号の間の到着の時間差（ＴＤＯＡ）ならびに到着角度が計算されて、初期の砲撃体距離を決定する。特定数の反復が実行され得、すなわち、引き続いて決定された更新された砲撃体距離の間の関係が収束基準を満足すると、更新された砲撃体距離が、最終の砲撃体距離であるとみなされる。例えば、収束基準は、引き続いて決定された更新された砲撃体距離の間の差が、または引き続いて決定された更新された砲撃体距離の間のパーセント変化が、所定値よりも小さいように選択されることができる。真の解を得るためには、算出された弾丸速度は、常に少なくとも音速になるように設定される。解は、整合性がチェックされる。例えば、弾丸軌跡角度および到着角度が所定値よりも大きいと決定されると、更新された砲撃体距離は無効であるとみなされる。
【００１３】
算出された砲撃体距離が、無効であると決定されても、解は、なお遺伝的アルゴリズム（ｇｅｎｅｔｉｃ ａｌｇｏｒｉｔｈｍ：ＧＡ）を適用することで得ることができる。例えば、所定数の個体を有するＧＡの初期の個体群が規定されることができ、各個体は、弾丸軌跡の推定した砲撃体距離、外れた方位（ＭＡ）、および外れた仰度（ＭＥ）を含む３タプルにより表現される。ＧＡは、所定数の世代に対して実行され、各世代における個体に対する残差が算出される。各世代において、最小残差を有する解が変異なく生存する個体として選択される。最小残差を有する解が、更新された砲撃体距離として選択される。解は、所定数の反復において各３タプルに対して実行されることで、洗練され得て、訂正した砲撃体距離を算出することができる。ここで、各世代における個体に対する残差が、訂正した砲撃体距離とともに算出される。
【００１４】
ＧＡは、クロスオーバおよび変異演算子を含む。クロスオーバ演算子は、世代の個体群からの２つの個体の間で、外れた方位および外れた仰度の少なくとも１つを交換する。一方、変異演算子は、（３タプルの値をランダムに選択した値に置換する）フィールド変異（ｆｉｅｌｄ ｍｕｔａｔｉｏｎ）、増加的変異（（ｉｎｃｒｅｍｅｎｔａｌ ｍｕｔａｔｉｏｎ）３タプルの全フィールドにおける微小変異を含む）、および（世代における個体を変更しないままとする）変異なしを備える。
【００１５】
本発明のさらに別の局面によると、衝撃波信号から、および限定された数の砲口爆風信号から、飛翔体軌跡を明確化するための方法であり、５以上の間隔を空けた音響センサにおける衝撃波単独信号を計測する工程と、多くても４つのセンサ上で砲口爆風信号を計測する工程と、前記衝撃波単独信号から各センサ対に対する到着の時間差（ＴＤＯＡ）情報を決定する工程とを含む。方法は、さらに、各個体が砲撃体方位、砲撃体仰度、外れた方位、および外れた仰度を含む４タプルで表現されている所定数の個体を含む初期個体群を用いて、所定数の世代に対して遺伝的アルゴリズムを実行する工程と、各世代における個体に対する残差を算出し、その残差がＴＤＯＡ衝撃波および砲口爆風信号の組み合わせの最小自乗適合を含んでいる工程と、を含む。最小残差を有する解およびその不明確代替解の比率が、所定値（例えば少なくとも２の値）よりも大きい場合は、最小の算出された残差を有する解は、明確化した飛翔体軌跡として指定される。
【００１６】
本発明のさらに別の局面によると、衝撃波信号の存在下で砲口波からの信号を抽出するための方法は、砲口波がセンサアレイを横断するのに要する時間に対応する時間ウィンドウを規定する工程と、前記衝撃波信号を検出する工程とを含む。衝撃波信号の検出に続いて、ウィンドウが適時に進行され、前記ウィンドウで受容した全エネルギーが進行時間（ａｄｖａｎｃｅ ｔｉｍｅ）の関数として計測される。計測された全エネルギーが砲口信号に関連づけられる。
【００１７】
スプリアス信号が衝撃波波形として解釈されることを防止するため、計測された衝撃波波形の音響エネルギーが、所定の周波数帯（例えば、約７００Ｈｚと１０ｋＨｚとの間の周波数）に亘って所定の閾値を下回ると、飛翔体軌跡は、誤りであるとして排除され得る。あるいは、またはそれに加えて、計測された衝撃波波形が正値を有する時間間隔が、最小時間を下回るか、最大時間を超える場合は、例えば、約７０μｓ未満か約３００μｓより大きい場合は、飛翔体軌跡は、誤りであるとして排除される。
【００１８】
有利な実施の形態において、全エネルギーは、好ましくは、衝撃波エコーにより生じた検出信号の部分を無視することにより、ウィンドウに亘って計測されたエネルギーを積分することにより決定されることができる。有利には、ピーク信号値は、最大の全エネルギーを生成するウィンドウにおいて決定することができる。ピーク信号値が、ウィンドウにおいて計測された全エネルギーよりも所定の率だけ大きい場合は、そのピーク信号は砲口爆風信号に関連していると認識され得る。
【００１９】
本発明の実施の形態は、１以上の次の特徴を含むことができる。アンテナおよび／または音響センサのタイミング誤差分布が、ゲイン変化、サンプリング変化、およびアンテナセンサのセンサ位置変化に関連し得る。明確化のための信頼度レベルは、アンテナのサイズに依存し、従って、比較的小さいアンテナには比較的高い計測精度が要求される。２つの不明確解が存在すると、明確化した飛翔体軌跡は、２つの不明確解に対する残差の比率に基づいて選択される。
【００２０】
別の有利な実施の形態において、各センサ対に対する到着の時間差（ＴＤＯＡ）が、最初に衝撃波に遭遇するセンサを基準センサとして指定する工程と、例えば、前記基準センサにおける衝撃波単独信号の初期部分の振幅が閾値を交差するとき、タイミング回路の第１のラッチを設定する工程と、により決定され得る。第１のラッチは、その他のセンサのそれぞれに対する開始カウンタを作動させる。このカウンタは、その他のセンサのそれぞれが衝撃波に遭遇するまで対応するセンサにおいて動作する。その他のセンサの１つが、例えば、衝撃波単独信号の初期部分に遭遇すると、これは、そのセンサに対する開始カウンタをそのセンサが停止する第２のラッチを設定する。基準センサに対するその他のセンサに対するＴＤＯＡ値が次いで記録される。
【００２１】
本発明の他の特徴と有利性が、次の好適な実施の形態の説明と特許請求の範囲から明確となるであろう。
【図面の簡単な説明】
【００２２】
本発明のこのような特徴および他の特徴および有利性は、添付の図面を参照して次の図示によりさらに充分に理解される。ここで、要素は、同様の参照符号で符号が付され、寸法が表わされているものではない。
【図１】図１は、アンテナとのマッハ円錐交差の断面図を概略的に示す。
【図２】図２は、７方向音響センサを有する代表的なセンサアレイの概略を示す。
【図３】図３は、衝撃波単独軌跡決定における固有の不明確さを概略的に示す。
【図４】図４は、マッハ円錐の曲率を決定するための、到着計測の時間差に対する確率密度を概略的に示す。
【図５】図５は、砲撃体軌跡間で正確に明確化する確率を概略的に示す。
【図６】図６は、相関過程の概略図を示す。
【図７】図７は、砲撃体軌跡間で正確に明確化するために使用される遺伝的アルゴリズムの処理フローを示す。
【図８】図８は、非衝撃波信号に対して区別するための処理フローである。
【図９】図９は、衝撃波ＴＯＡ（到着の時間）モデルの概略図である。
【図１０】図１０は、距離推定のための概略処理フロー図を示す。
【図１１】図１１は、距離推定のための遺伝的アルゴリズムの概略処理フロー図を示す。
【発明を実施するための形態】
【００２３】
発明の詳細な説明
上述の要旨に記載したように、本発明は、砲撃体距離推定および飛翔体軌跡の明確化のためのさまざまな実施の形態で、方法およびシステムを提供する。これらのシステムと方法は、正確な解に要求されるには不十分な数のパラメータが検出されたときや、そのようなパラメータを信頼して検出できないときに特に有用であり有利である。
【００２４】
超音速飛翔体軌跡は、近接する間隔の数個のセンサによって計測される飛翔体衝撃波到着時間からもっぱら推定される。このようなセンサは、アンテナと呼ばれる「微小な」計測容積の至る所に分布している。計測容積は、センサ間隔が２ｍ以下であれば微小とみなされる。一旦、飛翔体の軌跡が認識されると、砲撃体の位置は、軌跡に沿った距離（ｒａｎｇｅ）をさかのぼることなく知られる。アンテナが砲口爆風音の到着時間も得れば、この距離を見いだすことができる。しかし、砲口爆風は、いつも検出可能というわけではないので、正確な衝撃波単独解が軌跡の決定に重要である。
【００２５】
図１を参照して、衝撃波面が拡張してゆく円錐面を考慮する。この円錐面は、弾丸軌跡に一致してその軸を有する。衝撃波面は、マッハ円錐とも呼ばれる。衝撃波単独解を得るために、拡張する円錐面の到着角、曲率の半径、曲率半径の空間的勾配の３つの特性が、５個以上のアンテナセンサから計測された到着時間から決定されるべきである。
【００２６】
アンテナに最初に到着する円錐面起点（ｃｏｎｉｃａｌ ｓｕｒｆａｃｅ−ｇｅｎｅｒａｔｏｒ）の到着角度は、アンテナでの到着角度について弾丸軌跡の２つの可能性のある相対角度（しばしば「不明確」角度と呼ばれる）を決定する。この「不明確」角度については、図３を参照してさらに詳細に説明する。アンテナでの円錐面の曲率半径は、軌跡の距離と方向の両方を決定する。表面生成起点の経路に沿った曲率半径の勾配は、弾丸がどの方向に移動しているかを決定し、それにより、２つの可能性のある方向間で「不明確さ」を除去している。この２つの可能性のある「不明確な」軌跡角度の間で、これらの３つの衝撃波特性を正確に適切に決定するためには、非常に精緻な計測が要求される。例えば、２つの選択的な砲撃体アスペクト角間の確率的誤差は、約１μｓ以下とすべきである。
【００２７】
要求される精度は、図１に示す衝撃波の伝搬特性を考慮することで推定することができる。図２も参照して、アンテナ２０は、Ｎ個のセンサ（Ｎ＝７）を含み、進行する円錐衝撃波の到着時間を決定することが可能である。入ってくる弾丸軌跡は、どこからでも発生することが本質的に期待できるため、アンテナエレメント２３−２８は、有利には、（Ｃｘ０，Ｃｙ０，Ｃｚ０）の中心に位置する１エレメント２２を有して、球面を覆う位置Ｃ（Ｃｘｊ，Ｃｙｊ，Ｃｚｊ）に一様に分布することができるため、均一なセンサ開口が到着角度に無関係に提供されている。基準センサと呼ぶ第1のセンサが進行する円錐面を検出する時刻をｔ０とする。その他のセンサは、進行する円錐面をそれに続くｔｉとする時刻で検出する。進行する円錐面の方向での音の伝搬距離は、各時間差を局所の音速ｃで乗算することにより、すなわちｄｉ＝ｃ・（ｔｉ−ｔｉ）で得られる。もし、計測誤差が無ければ、Ｎ個の点の３次元座標が衝撃波円錐の全パラメータを理想的に決定して、その他の（Ｎ−１）個のセンサで基準センサを通過する円錐面を決定することもできる。しかし、上述のように、到着時間計測とセンサ座標における誤差が、衝撃波円錐に対する誤ったパラメータと、それによる飛翔体の軌跡の誤りをも生じ得る。次に、２つの別の不明確軌跡角度についての正確な決定をするために必要な到着精度の時間差について記載する。
【００２８】
このシステムは、有利には、車輌ノイズ、振動、風ノイズ、およびＥＭＩ等の砲撃体についての非弾道性信号を確実に間違えないようにする特徴を包含する。例えば、センサマストは、ガタガタ音を立てないよう嵌め合い接続の弾性スリーブで車輌（図示せず）に搭載することができる。センサは、背柱の端部に弾性カップリングで装着することができ、背柱の震動から隔絶するため約１Ｈｚの低周波共振を有する。センサの背柱は、アナログ電子回路を有する共通ハブに取り付けることができ、この共通ハブはまた、センサマストに弾性緩衝取り付けで装着してその振動から隔絶することができる。
【００２９】
さらに、衝撃波由来の信号に通常見られる特徴を欠く信号を除外するため、次の決定アルゴリズムを採用することができる。全ての値はパラメータ化され、すなわち、相対的であり、外部から調節することができる。リストに挙げた値は例示目的のみの提供である。
【００３０】
図８を参照して、プロセス８００は、検出された信号が衝撃波に由来するかを決定する。プロセス８００は、工程８０２で開始し、信号が衝撃としてカウントするのに充分な音量イベントか、例えば、ピーク信号値が所与のパラメータ化された、閾値（例えば５００）を超えるか工程８０４でチェックする。これが当てはまれば、工程８００は工程８０６へと続き、ゼロからピーク信号値への急峻な移行があるかをチェックし、このピーク値への移行が、重要な大きさ（例えばピーク信号値の１／１６）を有する別の信号によって先行したものではないことを確証する。
【００３１】
これが当てはまる場合は、工程８００は、工程８０８へと続き、衝撃波の最小と最大の間の時間が、充分に大きい値（例えば２００−４００μｓ）であるかをチェックする。これが当てはまれば、工程８００は、工程８１０へと続き、最小と最大ピーク信号振幅の大きさが、例えば３５％以内に互いに近接するかをチェックする。これが当てはまれば、工程８００は、工程８１２へと続き、最小ピーク信号からゼロへの圧力ピーク移行が急峻であるかを、実質的に工程８０６における同じ基準を使用してチェックする。これが当てはまれば、工程８００は、工程８１４へと続き、最大信号値とゼロ交差との間の時間、およびゼロ交差と最小信号値との間の時間が、例えば約１８０μｓ以内で比較可能であるかをチェックする。全ての工程で肯定的な応答を生じている場合は、工程８００は、信号が衝撃波である可能性があり、その信号は処理される（工程８１６）。逆に、もし６個の決定工程の１つが否定的な答えである場合は、検出された信号は、衝撃波から由来するものではない（工程８１８）。
【００３２】
図１に戻って、飛翔体軌跡がｘ軸と一致すると仮定する。マッハ角はΘ＝ａｒｃｓｉｎ（１／Ｍ）、ここでＭは、発射速度Ｖを音速ｃで除算して定義されるマッハ数である。Ｌは、アンテナの特性長を意味する。アンテナ２０の２つの端部における円錐の曲率半径は、ｒ１、ｒ２である。図の左半分の端面図は、曲率ｒ１の計測法を示す。距離ｄは、ｄ＝ｒ１・ｃｏｓ（Φ）に等しい。角Φは、ｓｉｎ（Φ）＝Ｌ／２ｒ１で規定されるので、微小角Φに対してΦ〜Ｌ／２ｒ１を得る。円錐面を半径ｒ１で二分するアンテナ面上の点間の曲率の時間差計測は、ｄｔ１＝Δｄ／ｃ＝（ｒ１−ｄ）／ｃ〜ｒ１Φ２／２ｃ＝Ｌ２／（８・ｒ１・ｃ）に等しい。ｒ２＝ｒ１−Ｌ・ｓｉｎ（Θ）における曲率の時間差計測がｒ１をｒ２に置換して同じ表現で与えられる。従って、ｄｔ２＝ｄｔ１＋Ｌ３ｓｉｎ（Θ）／８ｒ１２ｃである。
【００３３】
計測誤差にバイアスが無いとすると、すなわち、計測時間差ｄｔ１とｄｔ２が、異なる平均値ｄｔ１とｄｔ２であるが同じ統計的に決定された偏差σを有するランダムに分散した値であるとすると、アレイの２つの端部における平均計測値は、そこでの局所の曲率を正確に決定する。時間差ｄｔ１とｄｔ２に対する計測値の代表的な分布を図４に示す。
【００３４】
端部２でなされたサンプル計測をＸとして示す。端部２における曲率半径（半径ｒ２）は、端部１におけるもの（半径ｒ１）よりも小さい。従って、Ｘより大きい値を有する１の端部でなされた全計測は、端部１における曲率が端部２におけるよりも大きいという正確な決定が得られる。端部２における計測がＸに等しいとき、正確な決定がなされる確率は、
【数１】
で与えられる。
ｘに関して積分し、変数を代入すると、次の正確な決定を生成する確率が得られる。
【数２】
【００３５】
今度は図５を参照して、正確な決定の確率、すなわち明確化（ｄｉｓａｍｂｉｇｕａｔｉｏｎ）に対する信頼レベルを、２つの代表的なアンテナサイズＬ＝１ｍ、およびＬ＝２ｍについて、飛翔体の軌跡とアンテナ２０間の最接近点（ＣＰＡ）ｒに対してプロットしている。音速は、ｃ＝３４０ｍ／ｓと仮定している。大きなアンテナが不明確な衝撃波単独解に対する極めて拡大した範囲を有していることが明白である。大きなＣＰＡ値に対して、アンテナの２つの端部（ｒ１とｒ２）における曲率の差は、区別するには小さすぎるので、正確に決定する確率は、５０％に、すなわち完全に不明確に近づていく。従って、信頼レベルは、アンテナのサイズに、すなわちその径すなわち空間的広がりに依存する。
【００３６】
上述のように、誤差は、タイミングエラーとセンサ座標の不確定性から生じる。センサ座標の不確定性は、衝撃波到着角度の非常に可変的な関数であるバイアス誤差に寄与している。しかし、ランダム到着角度に対しては、センサ座標誤差は、ランダムな時間差誤差に見える。
【００３７】
タイミング誤差は、チャンネルからチャンネルへのゲインと信号強度の変化の両方からも生じる。到着の時間（Ｔｉｍｅｓ）は、センサ出力が所定の閾値Ｖ０に上昇するときに得られる。ゲイン変化ｄｇにより生じるタイミング誤差ｄｔは、そのチャンネルに対する時間当たりの電圧増加に依存する。ここで
【数３】
である。
【００３８】
タイミング誤差は、信号強度が開口に関して変化するときにも生じる。長さＬの開口で距離ｒにおける円柱形音源に対して、開口に亘る最大信号レベル変化は、ｐ０（Ｌ／２ｒ）に等しい。ここで、ｐ０は、開口中心における音圧である。上記のタイミング誤差式は、このタイプの誤差に対しても適用され、相対速度変分ｄｇ／ｇが式Ｌ／（２ｒ）に置換される。振幅誤差は、センサ間でランダムではなく、開口の全体に亘って最大から中心のゼロへと均一に変化する。１０ｍを超える距離では、１ｍ開口に対して、最大振幅係数は、０．０５未満であり、これは０．２のチャンネルゲイン変分パラメータよりも小さいので、振幅誤差による影響は無視することができる。逆に、上述のように、約１０ｍ未満の距離においては、マッハ円錐半径は、１ｍの長さの開口に対して充分に小さいので、計測誤差はあまり重要ではない。
【００３９】
センサの不確定性により生じるタイミング誤差の現実的な推定については、誤差ベクトルの大きさが統計的に独立しており、０と１ｍｍの間に均一に分散しており、誤差角度が統計的に独立であると仮定すると、均一に分散したランダム時間差誤差方程式の標準偏差は、
【数４】
に等しい。１ＭＨｚでサンプリングするシステムに対する二項分布したランダム時間サンプリング誤差の標準偏差は、０．２５μｓに等しい。ゲイン変化によるタイミング誤差は、約０．０２Ｖ／μｓの電圧変化に対応して、約１８ｋＨｚのチャンネル帯域幅を有する代表的なシステムについて約０．７５μｓと推定される。採用された各アレイに対する音響センサは、±１，５ｄＢ以内の感度を有するように選択された。従って、チャンネル相対ゲイン変化は、０．８４と１．１９の間にほぼ均一に分布しているため、相対ゲインの標準偏差は、ほぼ
【数５】
に等しい。閾値電圧は、Ｖ０＝０．１５Ｖであり、タイミング誤差の標準偏差は約０．７５μｓという結果となる。
【００４０】
全体の計測タイミング誤差は、チャンネルゲイン変化分、サンプリング変化分、センサ位置変化分が全て統計的に独立であるという仮定により推定される。次いで、タイミング誤差標準偏差は、
【数６】
として推定される。
【００４１】
このような精度をアナログからデジタルへの変換を伴って達成することは、高いサンプリング率とその後の補完が必要とされるため困難である。到着の時間差（ＴＤＯＡ）を正確に計測するための２つの異なる回路が開示のシステムに採用されている。
【００４２】
１つの実施の形態において、代表的なシステムは、各チャンネルにおいて１ＭＨｚクロックを使用するアナログの到着時間差（ＴＤＯＡ）回路を使用する。このクロックは、センサ信号が基準センサにおける閾値レベルを超えるとトリガされる。基準センサは、上記のように、センサが最初に衝撃波に遭遇するときに規定されたものである。上述のように、１ＭＨｚクロック速度は、実際の時間−サンプル誤差の重さを除去するために充分である。このシステムは、閾値の検出に依存してアナログモードで操作され、次の機能を実行するデジタル論理を有する。
１．衝撃波に最初に遭遇する基準センサにおけるチャンネル信号振幅が閾値を交差するとき第1のラッチが設定される。
２．第1のラッチが各チャンネルに対するスタートカウンタを設定する。これは各クロックサイクルにおいて１カウント増加する。プロセッサが警告を受ける。
３．各チャンネルのカウンタが、対応するセンサが衝撃波に遭遇するまで作動する。これは、そのチャンネルにおけるカウントを停止する第２のラッチを設定する。第２のラッチが設定されない場合は、対応するカウントは上限値まで作動する。
４．各カウンタの最後のカウントの数がデジタルＴＤＯＡレジスタに記録される。
５．プロセッサがＴＤＯＡレジスタを読み取る。
６．ブロセッサは、次の衝撃波を受容するためにカウンタをリセットする。
【００４３】
別の実施の形態において、ハードウェアＴＤＯＡ検出の時間を中央とした時間セグメントに対して、各チャンネルの１つおきのチャンネルとの相関が算出される。Ｃｏｒｒ（ｇ，ｈ）で示す２つの関数の相関は、
【数７】
で定義される。
【００４４】
この相関は、ｔの関数であり、「遅れ（ｌａｇ）」と呼ばれる。従って、時間領域にあり、次の特性を有する。
【数８】
ここで、ｇとｈは、時間の実関数である。Ｇ（ｆ）は、ｇ（ｔ）のフーリエ変換であり、Ｈ（ｆ）は、ｈ（ｔ）のフーリエ変換である。
【００４５】
信号の全パワーは、
【数９】
【００４６】
到着時間信号は、有限の長さを有するので、積分（または離散データに対する総和）は、到着時間のあたりを中心とした有限の時間間隔について実行されればよく、公知のように、１つまたは両方のチャンネルのデータの長さは、ゼロ挿入（ｚｅｒｏ ｐａｄｄｉｎｇ）により拡張され得るので、２つの信号の期間が整合する。
【００４７】
次の議論において、連続関数の積分が簡単のために使用されるが、実際のデータはデジタル化された離散値である。当業者は、容易に積分を総和に置き換えるであろう。
【００４８】
図６を参照して、工程６０において、衝撃波信号時間データｇｉ（ｔ），ｇｊ（ｔ）が、各チャンネルｉ，ｊ、工程６０１，６０２において取得され、時間の関数として記録される。工程６０３、６０４において、チャンネルｉにおける全信号パワーは、次の相関の正規化に対して算出される。
【数１０】
【００４９】
チャンネルｉに対する衝撃波信号時間データｇｉ（ｔ）のフーリエ変換Ｇｉ（ｆ）が計算され、その共役Ｇｉ（−ｆ）が形成される（工程６０５）。同様に、衝撃波信号時間データｇｊ（ｔ）のフーリエ変換Ｇｊ（ｆ）が、その他の全てのチャンネルｊに対して算出される（工程６０６）。その後、相互相関Ｇｉ（−ｆ）・Ｇｊ（ｆ）が、各チャンネルペア（ｉ，ｊ）に対して形成される（工程６０８）。これは、「遅れ」ｔの関数ｆｉ，ｊ（ｔ）である。各チャンネルペアに対するＴＤＯＡは、ｆ（ｔ）が最大値を有する時間ｔｍａｘである（工程６１０）。チャンネルｉとｊの間の相関は、
【数１１】
として定義することができる。
【００５０】
チャンネルｉの残差は、工程６１２に示すように、全てのセンサｊに亘ってセンサｉに対する平均値を計算することで算出される。
【数１２】
最良（すなわち最小）の全残差を有するそのチャンネルに対するＴＤＯＡおよび相関は、次いで「最善」の解として選択される（工程６１４）。
【００５１】
上述のように、チャンネルデータは、通常、例えば４１，６６６．６６サンプル／秒の所定のサンプリングレートの離散時間間隔でサンプリングされる。これは、２４μｓのビン幅に対応し、受信信号に対する時間分解能を反映している。相関処理は、３３３３３３サンプル／秒を取ることで３μｓに対して８の係数を取ることで改善された時間分解能でなされる。
【００５２】
一旦センサ間の様々な到着時間差（ＴＤＯＡ）が衝撃波単独の信号から決定されると、砲撃体方位と仰度、および弾丸軌跡を決定することができる。砲撃体位置、すなわちセンサアレイから砲撃体の距離は、砲口爆風信号が既知であればさらに決定することができる。
【００５３】
アレイの中央を中心とする直交座標系において、すなわち｛（Ｃｘ０，Ｃｙ０，Ｃｚ０）＝（０，０，０）｝、所与のセンサ（Ｃｘｊ，Ｃｙｊ，Ｃｚｊ）（図２参照）における衝撃波の到着時間（ＴＯＡ）は、
【数１３】
で与えられる。
ここで、ｃは音速であり、Ｍはマッハ数である。βは、砲撃体位置と弾丸軌跡の「外れ（ミス）角度」を表わし、これには方位角と仰角の両方が含まれる。直接のヒットであれば、β＝０に対応する。マッハ角度Θは、
【数１４】
で定義される。
【００５４】
上述および図３に示すように、所与の砲撃体位置と弾丸軌跡に対して、所与の一式のセンサにおける衝撃波のＴＯＡがほぼ等しい別の砲撃体位置と弾丸軌跡が存在する。２つの不明確な解は、簡単化したモデルにおいては事実上同一であり、衝撃波は、平面波としてセンサアレイを横断して伝搬するものと思われる。ＴＤＯＡの分解能が、衝撃波の曲率を分解するのに充分高ければ、２つのほぼ同一の解は明確化することができる。本質的に不明確な衝撃波単独ＴＤＯＡ解を図３に示す。
【００５５】
充分に正確なＴＯＡ計測を仮定すると、砲撃体位置と弾丸軌跡の真の解は、計測および算出された衝撃波ＴＤＯＡの残差の二乗平均平方根（ＲＭＳ）を最小化する砲撃体／軌道の組み合わせを算出することで得ることができる。
【数１５】
ここで総計は全センサに関してなされる。
【００５６】
この問題を解く１つのアプローチは、米国特許第５，９３０，２０２号にその詳細が記載のＬ１Ｌｅｖｅｎｂｅｒｇ−Ｍａｒｑｕａｒｄｔアルゴリズムである。通常の古典的逐一アルゴリズムは、最適解にアプローチするために決定論的処理を使用し、当て推量解から開始して事前に規定された移行則に基づいてサーチ方向を指定する。このようなものには、目的関数と制約値を使用した直接法、ならびに１階および２階微分を使用した勾配法がある。しかし、このような方法は、例えば、最適解が選択された初期解に依存し、アルゴリズムが、局所最小値や費用関数面が平坦な谷を有する所等の準最適解で「スタック」するという欠点を有する。このためさらに繰り返しても結果が改善しない。
【００５７】
進化型（ｅｖｏｌｕｔｉｏｎａy)遺伝的アルゴリズム（ＧＡ）を使用することで、砲撃体方向と飛翔体軌跡の全体の最小値を、より速くより信頼度を高く明確化して計算することができることがわかった。ＧＡは、自然の進化原理を模倣し、サーチや最適化処理に適用される。
【００５８】
ＧＡの概略フロー図を図７に示す。解に対して単一の推量から開始する代わりに、ＧＡ工程７０は、解のランダムな個体群を初期化すること（工程７１）でそのサーチを開始し、初期解集合を示すゼロに世代カウンタを設定する（工程７２）。一旦、解のランダム個体群が生成されると、それぞれは非線形計画法問題のコンテキストにおいて評価され（工程７３）、適性（優劣）が各解に割り当てられる（工程７４）。この適性は、計算解と計測解の間のユークリッド距離Δτｍｉｎによって表わすことができる。
【数１６】
【００５９】
直感的には、小さい値のΔτｍｉｎを有するアルゴリズムが良い。
【００６０】
例えば、ＧＡを砲撃体方向と飛翔体軌跡に対する解を明確化するために適用するときは、代表的なＧＡは、２００ ４−ｓの初期個体群を染色体として使用する。ここで、各４−は、次の値を有する。
［方位砲撃体、仰度砲撃体、方位外れ、仰度外れ］，［方位砲撃体，仰度砲撃体］は、角度（Θ＋β）により規定され、ここで、［方位外れ、仰角外れ］は、角度β（図３参照）で規定される。砲口爆風は、先述の衝撃波単独アプローチを使用しないので、センサアレイと砲撃体間の１００ｍの基準距離を仮定する。
【００６１】
初期個体群は、重要で合理的な範囲に亘る値（全ての値は角度である。）である４−ｓのランダムな選択により生成され、
方位砲撃体＝｛０，．．．．，３６０｝，
仰角砲撃体｛−１０，．．．．，３０｝，
方位外れ｛−２０，．．．．，２０｝，および
仰角外れ｛−２０，．．．．，２０｝である。
【００６２】
工程７５において、例えば２５に設定可能であるＧＡの反復の最大数に達したかがチェックされる。反復の最大数に達していた場合は、工程７０は、工程８０で停止し、結果が許容できるか別に評価されるかのいずれかであり得る。そうでない場合は、工程７６は、事前に設定した適性基準が満足されたかどうかをチェックする。
【００６３】
適性基準は、例えば、計算された＜１５°の外れた方位および／または２つの不明確解の残差の割合であり得る。適性基準が満足されている場合は、処理７０は工程８０で停止し、そうでなければ、新たな個体群が、クロスオーバ（工程７７）、変異（工程７８）を通じて作成され、生成カウンタが1刻みで増加する（工程７９）。
【００６４】
各世代において、「最良の」個体は、変異することなく生存を許され、それに対して、トップ１００の個体も、その適性により判断されて生存するが、表１に列挙したクロスオーバ／変異演算子を伴って、これらの生存体のペアから次の１００の個体を作成するために使用される。
【００６５】
次の代表的なクロスオーバおよび変異演算子が、工程７０を実行するために使用された。
【表１】
【００６６】
明確化は、最良解と対応する代替解の勾配サーチ（ｇｒａｄｉｅｎｔ ｓｅａｒｃｈ）を実行することで達成および／または改善される。両方の不明確な解に対して、残差と残差の比率が算出される。算出された外れた方位が＜１５°の場合、「近い」砲撃であり、残差の比率が＞２の場合は、低い方の残差を有する解が選択される。さもなくば、実際の選択がなされず、低い方の残差を有する解が「主要」解と標示され、他方の解は「代替」解と標示される。
【００６７】
衝撃波単独検出に際して、ＧＡアルゴリズムは、Ｌｉｎｕｘオペレーティグシステムを作動する１ＧＨｚコンピュータで０．１５秒で広範なシミュレートされた砲撃において１つの解を生成した。９７％のシミュレートされた砲撃が外れた方位の１５°以内で、シミュレートされた砲撃の８６％が外れた方位の５°以内であった。先述の明確化アルゴリズムを使用して、近い砲撃、すなわち１５°未満の外れた方位を有する砲撃は、その時の９５％が明確化された。明確化アルゴリズムは、その時の７０％の遠い砲撃に対して正確な結果を生成した。明確化の正確性さは、センサアレイ配置と推定された砲撃の分布に基づいて変化すると予測され、低い仰度を有する砲撃は明確化が容易である。
【００６８】
飛翔体軌跡に対する先述の解は、砲口爆風を検出することなく得られた。しかし、弱い砲口信号または限定された数のチャンネル上のみにおいて受信された砲口信号でさえ、距離決定と明確化の改善のために有利に使用することができる。
【００６９】
図９は、到着時間（ＴＯＡ）モデルを概略的に示す図であり、米国特許第６，１７８，１４１号にさらに詳細が記載されている。ＴＯＡモデルは、センサ位置に対する飛翔体の軌跡と砲撃体方向を推定するために使用することができる。ＴＯＡモデルは、飛翔体の飛程に関する特定の物理特性を考慮した弾道モデルに基づいている。このような物理特性には、空気密度（温度に関連）、砲撃体の位置
【数１７】
ライフル砲口の方位角および仰角、飛翔体の砲口速度（または相当のマッハ数）、および音速（温度／空気密度により変化）等がある。この弾道モデルを使用して、衝撃波と砲口爆風が空間の特定の点に到達する時間を正確に計算することができる。
【００７０】
図９に示すように
砲撃体は、原点（０，０，０）に対して点
【数１８】
に位置し、多様なセンサが点
【数１９】
に配置され、弾丸軌跡は、
【数２０】
の方向に砲撃体から生じているとして示される。砲撃体からＪ番目のセンサへのベクトルは、
【数２１】
であり、ｊ番目のセンサへの弾丸の最接近点（ＣＰＡ）は、
【数２２】
であり、衝撃波が軌跡からｊ番目のセンサへ生じる点からの経路は、
【数２３】
（センサの添字ｊは省略）である。
弾丸のマッハ角は、θ＝ｓｉｎ−１（１／Ｍ），Ｍ＝Ｖ／ｃ０である。Ｍは飛翔体のマッハ数であり、Ｖは飛翔体の超音速であり、ｃ０は、音速（圧力と温度に依存）である。軌跡とｊ番目のセンサの間の「外れ角度」はβである。軌跡は、ｘ−ｙ平面においてＸ軸から反時計回りに計測された方位角によって、およびｘ−ｙ平面から上方に計測された仰角によって特徴づけられる。衝撃波到着時間ｔｊとｊ番目のセンサにおける単位ベクトルを規定する方程式は、これらの幾何学的量という形で記述される。
【００７１】
到着の時間は、飛翔体が、ｊ番目のセンサに向けて音が放射される点までの距離
【数２４】
を進むのに要する時間
【数２５】
に、衝撃波が放射点からｊ番目のセンサまでの距離
【数２６】
を進むのに要する時間
【数２７】
を加えたものに等しい。
【数２８】
ここで、ｔ０は、時間基準（発射時間）であり、ｃ０は、音速である。マッハ角Θも、図９に示す。
【００７２】
飛翔体の速度Ｖは、センサ間隔に対応する距離に亘って一定のままであると安全に仮定することができるので、飛翔体が異なるセンサへと放射する時間の間では大きな速度の損失は存在しない。しかし、比較的長距離に亘っては、飛翔体は空気抵抗により減速することが知られている。空気抵抗は、抗力係数Ｃｂで表現することができ、これは弾丸形状と弾丸口径に依存する。物理的原理から得られる数学的弾道モデルは、空中における任意の点での衝撃波の到着時間を、あらかじめ既知の、飛翔体を記述する全パラメータ一式（例えばその抗力係数Ｃｂにより）、その初期速度、および周囲空気の密度の関数として予測することができる。
【００７３】
正確な計算に必要なパラメータは、戦場などの実際的な設定においては通常は既知ではない。しかし、距離推量は、図１０の処理フロー図２００の形態に示す処理を繰り返すことにより、かなり改善することができる。この処理においては、軌跡に沿った飛翔体の速度の減速を考慮している。処理２００は、次の仮定、
ｃ０＝外気温度／空気圧に対して修正された音速（〜３４０ｍ／ｓ）
Ｃｂ＝予測の武器に関して平均した公称抗力係数
Ｖ０＝予測の武器に関して平均した、飛翔体の発射時の初期速度
Ｍ０＝Ｖ０／ｃ＝飛翔体の初期マッハ数
【００７４】
砲撃体距離Ｄ０の第1の推定は、方程式
【数２９】
に従って、計測された到着の時間差（ＴＤＯＡ）τｍｓ、およびセンサアレイにおける衝撃と砲口音の間の到着角度αを使用して、ならびに初期、定常速度Ｖ０、およびマッハ数Ｍ０を仮定することにより、工程２０４において計算される。
【００７５】
これらの仮定とともに、砲撃体位置
【数３０】
からの距離における飛翔体の速度が、工程２０６において方程式
【数３１】
から、計算され得るので、飛翔体が軌跡に沿った距離ａを移動する時間（工程２０８）は、
【数３２】
となる。
【００７６】
角度Θは、方程式
【数３３】
により、マッハ数Ｍａに関連する。ここで、マッハ数Ｍａは、初期的にはＭ０に設定される。瞬時の弾丸速度は、算出された弾丸速度が音速よりも小さくなる場合は音速（すなわちＭａ＝１）に設定されることに注意すべきである。工程２１０における訂正した距離
【数３４】
は、
【数３５】
となる。
【００７７】
角度α、β、およびΘは、方程式（α＋β＋Θ）＝９０°という関係にある。処理２００は、次いで上記のＭａとＴａについての方程式において算出した距離ａを挿入することにより工程２０６へループバックし、移動距離ａに対して更新したマッハ数Ｍａと更新した弾丸移動時間Ｔａをそれぞれ得る。計測されたＴＤＯＡτｍｓとＴａとａに対する算出した更新値が次いで砲撃体距離に対する値Ｄを連続的に更新するために使用される。
【数３６】
【００７８】
この処理は、工程２１２で決定するように、反復回数の最大数が達したか距離値Ｄが収束するかのいずれかまで繰り返される。
【００７９】
処理２００は、砲撃体とセンサアレイ間の距離に対する訂正した距離値
【数３７】
が「合理的な」値であるかも工程２１４でチェックする。該当する場合は、処理２００は工程２１６で終了する。例えば、飛翔体が移動する距離ａとセンサと音波が飛翔体からセンサへ放射される点までの間の距離
【数３８】
とが有効な数値である場合、すなわちＮＡＮ（非数）でない場合は、距離Ｄに対する値は有効であると考えることができる。ＮＡＮは、数値演算の結果を有効な数値に戻すことができないことを示す特別な浮動小数点値であり、通常は、計算中のエラーの伝搬を阻止するために使用される。加えて、αとβは両方とも所定の閾値未満である必要があり、これは、飛翔体が確かにセンサアレイに向けて発射されたことを示すものである。
【００８０】
上述のように、数値対（τｍｓ、α）が零次近似の砲撃体距離Ｄ０の算出に初期的に使用される。これにおいては、軌跡に沿った飛翔体の速度の変化が無視される。上述の反復処理２００が数値対（τｍｓ、α）を支持する一致構造を返さない場合は、その解は廃棄される。
【００８１】
正確な解が得られない場合でも、計測衝撃ＴＤＯＡと計測砲口ＴＤＯＡに最も近く整合する、砲撃体距離Ｄ、外れた方位角と仰角（これらはβに関連する。）に対する値を見つけることが目標である。既述のように、多様なセンサ間の衝撃単独ＴＤＯＡが、通常の状況において容易に計測することができる。砲撃体方位と砲撃体仰度は、既知のセンサアレイ座標（Ｓｘｊ，Ｓｙｊ，Ｓｚｊ）を使用して衝撃単独ＴＤＯＡから決定することができるが、砲撃体距離についてはできない。検出された衝撃波と砲口音の間のＴＤＯＡτｍｓも計測され得るが、砲口音はそれによっては全てのセンサにより検出することができないことが考えられる。
【００８２】
反復処理２００が有効な結果を返さないと工程２１４において決定されると、次いで処理２００は、進化型遺伝的アルゴリズム（ＧＡ）３００を呼び出すことによって砲撃体距離を算出する。ＧＡは、自然の進化原理を模倣し、サーチや最適化手順に適用される。ＧＡは、そのサーチを、１つのみの解ではなく、ランダムな解の集合から開始する。一旦、ランダムな解の個体群が生成されると、各々は、非線形計画問題のコンテキストにおいて評価され、適性（優劣）が各解に対して割り当てられる。１つの実施の形態において、この適性は、計算解と計測解間のユークリッド距離で、例えば、
【数３９】
で表現することができる。
【００８３】
直感的には、比較的小さい値τｍｉｎを生成するアルゴリズムが良い。
【００８４】
ＧＡ処理３００の概略フロー図を図１１に示す。処理３００は、工程２００に対して先に計測された到着の時間差（ＴＤＯＡ）τｍｓと角度αを使用する（工程３０２）。値｛ＲＡＮＧＥ（距離）、ＭＡ、ＭＥ｝を有する３タプル（３−ｔｕｐｅｓ）の典型的な数が初期個体群として規定される（工程３０４）。ここでＲＡＮＧＥは、図９に示す砲撃体距離
【数４０】
であり、ＭＡは、外れた方位であり、ＭＥは外れた仰度である。ＭＡとＭＥ値は、方位と仰度空間において弾丸が標的を外れた程度を示す。図示の例の標的は、センサアレイと仮定している。工程３０４における初期個体群は、ランダムな選択の３タプル（３−ｔｕｐｅｌｓ）により生成され、この中には重要で合理的な値の範囲
距離砲撃体（Ｒａｎｇｅｓｈｏｏｔｅｒ）＝｛１０００，．．．，３０００｝［メートル］
方位外れた（ＡｚｉｍｕｔｈＭｉｓｓｅｄ）＝｛−２０，．．．，２０｝［度］、および
仰度外れた（ＥｌｅｖａｔｉｏｎＭｉｓｓｅｄ）＝｛−２０，．．．，２０｝［度］
が含まれている。
【００８５】
算出は、衝撃波単独解に対する上記概説したものと同様の処理に従う。初期的には、世代Ｇｅｎ＝０に対して（工程３０６）、砲撃体位置ベクトル
【数４１】
が、先に決定された砲撃体方位と仰度、および特定の３タプルに対する仮定した距離（ＲＡＮＧＥ）を使用して、各３タプルに対して算出される。初期のマッハ数をＭ０とすると、ベクトル
【数４２】
すなわち、衝撃音が放射する位置が、各３タプルに対するＭＡとＭＥ値を用いて計算される（工程３０８）。砲撃体と衝撃波を検出した各センサｊ間の距離
【数４３】
も算出される。
【００８６】
各３タプルに対して、角度βが方程式
【数４４】
から算出される。ここでシンボル「・」は、２つのベクトル間のスカラー積を示す。βに対して算出した値、およびＭａ＝Ｍ０に対して初期的に仮定した値、およびａを、Ｍａ、Ｔａ、ａおよびＤに対する上記の方程式に挿入することにより、距離ａに対して更新された値、飛翔体の距離ａに亘る移動時間Ｔａ、およびマッハ数が算出される（工程３１２）。この処理は、各３タプルに対して数回、例えば工程３１２で決めているように３回繰り返され、その後、砲口信号を含む上記で規定された残差Δτｍｉｎが各３タプルに対して算出される（工程３１４）。
【００８７】
工程３１６において、ＧＡに対する反復の最大数（例えば２５反復）が達成されたかがチェックされる。反復の最大数に達した場合は、処理３００は工程３２０で停止し、最小残差を有する３タプルを返す。達していない場合は、処理３００は、新たな個体群をクロスオーバと変異演算から作成し（工程３１８）、世代カウンタが１増加する（工程３２２）。
【００８８】
各世代において、「最良の」個体は、変異することなく生存を許され、それに対して、トップ１００の個体も、その適性により判断されて生存するが、以下の表２に列挙したクロスオーバ／変異演算子を使って、これらの生存体のペアから次の１００の個体を作成するために使用される。
【００８９】
次の例示のクロスオーバと変異演算子が処理３００を実行するために使用された。
【表２】
【００９０】
ＧＡ処理３００が、２００の異なる３タプルの初期個体群で、５０の補充率（ｒｅｆｉｌｌ ｒａｔｅ）で、合計２５世代に対して実行される。ＧＡは、異なる初期の３タプルのセットと並行して５回作動し、最小の残差を有する解が、砲撃体の距離（ＲＡＮＧＥ）、外れた方位、および外れた仰度に対する最終的解として選択され、これがベクトル
【数４５】
の算出を可能とする。
【００９１】
最近の実験的施行によると、５以上の衝撃波チャンネルに加えて少なくとも１つの砲口信号チャンネルを使用することにより、同じデータセット上で、不明確な砲撃が９５％から８％へと減少を示した。これは、衝撃波単独解に関して重要な改善である。
【００９２】
空気抵抗による飛翔体の経路に沿った減速を考慮しない計算は、過大評価した距離となる傾向がある。特定の幾何的配置と充分に遠い砲撃に対して、この過大評価は２０％を超えることもある。先述の処理は、長距離砲撃検出に対する距離推定からこのバイアスを除去する。
【００９３】
上述のように、２つの異なる軌跡から残差を比較し、小さい方の残差を有する軌跡を選択することで、不明確な衝撃波単独解は、明確化され得ることがしばしばである。
【００９４】
砲口爆風信号は、４以上のセンサチャンネルで検出されると、衝撃チャンネルの数に関係なく、先述の衝撃砲口アルゴリズムを使用して、明確に砲撃体位置を決定することができる。砲口爆風信号が４未満のセンサで検出されるが、衝撃波信号が５以上の衝撃波チャンネルで検出されると、修正した費用関数または残差を有する先述のＧＡを使用して、どのような砲口信号が得られても最適化関数に「混合」されて衝撃波単独解を明確化し、および／または砲撃体距離の推定精度を上げる。しかし、３未満の砲口チャンネルおよび５未満の衝撃波チャンネルが検出されると、砲撃体位置を特定しようとすることなく警告が起動する。
【００９５】
砲口信号は、全てのチャンネルで確実に検出されるわけではない。その理由は、
１．１以上のチャンネルでの検出レベルでは、確信的に検出するには低すぎる。
２．砲口エネルギーが生信号中では区別できず、システムに関連するノイズを生じさせ、ＴＤＯＡ推量の信頼度を低くする。
３．衝撃波からのエコーが砲口爆風よりも強い場合があり、砲口爆風よりも早く到着することがあり、システムが衝撃を砲口として誤検出することがある。
【００９６】
数チャンネルで検出されたのみの砲口爆風信号に関して、この状況での残差は、
【数４６】
として定義することができる。ここで、砲口爆風の最初の項は、砲口爆風を検出する減少した数のセンサ（＜４）について総計されており、ｊは、衝撃波を検出するセンサ（通常は全センサ）に関して総計されている。
【００９７】
上述の例に示すように、砲口爆風信号は、衝撃波単独解に比べて、砲撃体方位とそれによる飛翔体の軌跡についての重要な情報を提供するので、算出された軌跡の解は、不明確な解の１つとより近くに整合する。すなわち、このようにして解を明確化している。
【００９８】
少なくともいくつかの信頼できる砲口信号が無ければ、特に長い砲撃体距離においてかなりの数の不明確な衝撃波単独解が生成され得るがあまり正確な解ではなく、これは比較的少ない数の明確なものほど望ましいものではない。
【００９９】
信頼度が無いかもしれない砲口検出の場合は、砲口爆風信号を検出する試みが最初はなされ得る。例えば、ノイズの多い信号の中の砲口爆風痕跡を見いだして、得られるＴＤＯＡを算出することである。チャンネル間で充分な相互相関を有して充分な数のセンサで砲口信号が見い出される場合であって、砲口信号と各チャンネルの対応する未処理のバンド（フィルタ遅延を補償するために数ビンだけオフセットされている。）との間に充分に強い相関がある場合は、砲口検出は信頼度があるものと考えられる。
【０１００】
そうでない場合は、不十分な相関を示す砲口信号が少なくとも消去され、次の「粗い砲口検出」ロジックが呼び出される。
−衝撃に続く衝撃エネルギーのピークを探す。これらのピークを「衝撃エコー」の可能性ありとしてフラグし、これによって砲口爆風としては排除する。
−砲口波がセンサアレイを横断するのに要するであろう時間を決定し、対応する所用時間を有する「ウィンドウ」を定義する。検出された衝撃波に続いて、衝撃エコーとして認識された検出信号の部分を無視しつつ、実質的に全部の検出器チャンネルに亘ってこのウィンドウを動かすことにより砲口エネルギーピークを探す。ウィンドウに関してエネルギーを積分する。すなわち、
【数４７】
の最大値を探す。
ここで、
【数４８】
の２乗は、ｎ番目センサによって計測された砲口爆風のエネルギー等のエネルギーの測度を示し、ｉ＋ｊは、検出チャンネルを示し、ｉは、衝撃波が検出された時間とウィンドウの開始との間の離散時間間隔を示し、ｊは、ウィンドウの開始からの計測された時間間隔を示す。
【０１０１】
ノイズに対して区別するため、関数ｆｍａｘ（ｉ）の最大値を生成するウィンドのピークエネルギーがチェックされ、その最大値におけるエネルギーピークが、ウィンドウに亘るエネルギーよりも所与の比率係数だけ大きいかを決定する。該当する場合は、ウィンドウの信号が砲口検出として認識され、砲口爆風バンドｍｂの全チャンネルで相互相関が実行されて砲口ＴＤＯＡを決定する。
【０１０２】
検出された砲口爆風信号は、次いで上述のように、砲撃体距離の決定および／または衝撃波信号の明確化のために使用される。
【０１０３】
要するに、記載のシステムは、衝撃波単独計測に基づいて、正確に、迅速に、そしてしばしば明確に、砲撃体方向と弾丸軌跡を提供する。明確化は、改善することができ、砲撃体距離は、微弱な砲口爆風波形が検出される場合でも推定することができる。システムは、車輌震動やノイズ、風音にはそうではないが、花火やシステムから遠ざかる方向の近くの発砲、に応答する誤った砲撃体指示に比較的鈍感である。
【０１０４】
衝撃波信号を検出するシステムは、信号が本当に衝撃波に由来するかを決定するため、初期波形上で２つのテストを実施する点が記載されるべきである。第１に、約７００Ｈｚと１０ｋＨＺの間の周波数帯における計測された全エネルギーを経験的閾値と比較する。この閾値を超える場合に限り、信号形態は衝撃波から生じているものと考えることができる。第２に、検出された初期正圧ピークの時間間隔が、約７０μｓより大きく約３００μｓより小さくなければならない。これらの基準は、花火や威嚇性のない発砲などのインパルス的ノイズからシステムの免疫を提供する。これらのテストを通過しない場合は、検出された波形は、衝撃波ではないと考えられ、砲撃体解は試行されない。
【０１０５】
本発明は、図示および詳述な記載の好適な実施の形態に関連して開示してきたが、これに関する多様な改変および改善が本発明の精神と範囲から逸脱することなくなされ得る。例えば、例示の実施の形態においては、マイクロフォン等の音響センサを有するとして図示しているが、これはそうである必要はない。代わりに、別のタイプの圧力感応性機構や電気的センサが使用され得る。さらに、表１と２の多様な演算子のための値は、例示のみのための意図であり、実際の現場の状況に応じて別の値が選択されることがあり得る。従って、本発明の精神と範囲は、添付の特許請求の範囲のみによって限定されるべきものである。
【技術分野】
【０００１】
政府契約
合衆国政府は、本発明において、支払い済みのライセンスと、限定された状況において、ＤＡＲＰＡＡＴＯにより裁定された契約番号第ＨＲ００１１−０４−Ｃ−００３５の条項によって規定される適切な条件で特許権者に他にライセンスすることを要求する権利とを有する。
発明の分野
本発明は、法的措置技術および安全に関し、より詳細には超音速飛翔体の発生源と移動の方向を決定するための方法とシステムに関する。本方向およびシステムは、砲撃体（ｓｈｏｏｔｅｒ）とセンサ間の距離が大きいときでも、砲口音からの信号が無いか微弱なときでも、砲撃体位置を決定し明確化することができる方法とシステムである。
【背景技術】
【０００２】
飛翔体によって生成された衝撃波に伴うパラメータを計測することで、弾丸や砲弾等の超音速飛翔体の概略の方向と軌跡を決定することができるシステムと方法が公知である。このようなシステムの１つは、米国特許第５，２４１，５１８号の記載によると、少なくとも３つの間隔を空けたセンサを含んでおり、各センサには１つの平面に配置された３つの音響トランスデューサが組み込まれている。センサは、その発生源に対する方位と仰角に関連する衝撃波に反応する。衝撃波単独計測では、センサと衝撃波の発生源間の距離を決定することができない。通常、距離情報は、火光や砲口爆風から得られる。
【特許文献１】米国特許第５，２４１，５１８号公報
【０００３】
センサ位置に対する砲撃体の方位と仰角は、通常は角センサにおける砲口信号と衝撃波信号の「ＴＯＡ」すなわち到着時間（Ｔｉｍｅ ｏｆ Ａｒｒｉｖａｌ）情報を計測することで決定される。各センサは、異なる時間において信号と遭遇して砲口および衝撃波圧力に応答した信号を生成する。各種のセンサからの信号が処理され、砲口および衝撃波の発生源に対するセンサからの方向（方位および仰度）とともに飛翔体の軌跡を決定することができる。
【０００４】
従来のシステムでは、マイクロフォンが採用され、これは、比較的近接したり（例えば１ｍ間隔等）、広く分散（例えば、車輌に搭載されたり、戦場の兵士が持ち歩いたり等）しており、それぞれの位置で砲口および衝撃波圧力を全方位的に計測する。しかし、センサが比較的広く間隔が空けられ、および／または軌跡がアンテナ内にあると、衝撃波単独解を正確に得るのに必要なタイミング精度が非常に高度となり、特別な技術が要求される。
【０００５】
大きなアンテナサイズは、例えば車輌搭載システムにおいては大変不利である。さらに、小さい時間分解能しか有さないシステでは、不明確な解を出し得、これでは、所与の一式のセンサにおける、２つの鏡面対称の砲撃体位置に対する衝撃波のＴＯＡ時間情報がほぼ同一となる。
【０００６】
従来のアルゴリズムは、少なくとも４つの衝撃波および砲口検出を要求するので、衝撃波ＴＯＡに関して、４×４行列が反転されて平面波がマッピング可能とする。衝撃および砲口ＴＯＡ決定における微小な誤差が、距離（ｒａｎｇｅ）推定における大きな誤差を生み出し得る。さらに、従来のアルゴリズムは、弾丸軌跡に沿って一定の弾丸速度を仮定しており、これでは、約３００ｍを超える距離から発射されている長距離弾に対しては、距離推定が不正確となる。
【０００７】
従って、長い砲撃体距離を正確に推測することを可能とする迅速に変換するアルゴリズムが必要とされている。また、砲撃体方向に対する衝撃波単独解を明確化する必要もあった。さらに、砲口爆風に関係のない音響的サインに隠され得る砲口信号を抽出する必要がある。
【発明の開示】
【０００８】
発明の要旨
本発明は、様々な実施の形態により、従来技術の欠点に対処したものであり、特に、砲口信号が弱いか検出チャンネルにおいて検出された数が不十分であるときの、長距離に対する砲撃体距離（ｓｈｏｏｔｅｒ ｒａｎｇｅ）の推定のための方法およびシステムを提供する。開示の方法およびシステムは、最適化処理中で、微弱なおよび／または信頼度の低い検出砲口信号を含むことにより、さらに改善を達成して衝撃波単独砲撃体軌跡解の明確化も改善する。
【０００９】
本発明の１つの局面によると、衝撃波単独信号から飛翔体軌跡を明確化するための方法であって、アンテナを形成する５以上の間隔を空けた音響センサにおいて、衝撃波単独信号の少なくとも初期部分を計測する工程と、前記音響センサに対するタイミングエラー分布を推定する工程と、前記衝撃波単独信号の計測された初期部分から、センサ対に対する到着の時間差（Ｔｉｍｅ−Ｄｉｆｆｅｒｅｎｃｅ−ｏｆ−Ａｒｒｉｖａｌ：ＴＤＯＡ）を、前記推定されたタイミングエラー分布より大きい時間分解能で決定する工程と、明確化に対する規定された信頼度レベルに基づいて、および前記音響センサのＴＤＯＡに対する残差の値に基づいて、前記明確化した飛翔体軌跡を選択する工程と、を含む方法である。
【００１０】
本発明の別の局面によると、衝撃波単独信号から飛翔体軌跡を明確化するための方法であって、アンテナを形成する５以上の間隔を空けた音響センサにおいて、衝撃波単独信号の少なくとも初期部分を計測する工程と、各センサ対に対して、計測された衝撃波単独信号初期部分から到着の時間差（ＴＤＯＡ）を決定する工程と、遺伝的アルゴリズムを、飛翔体軌跡想定を備える初期染色体に、所定の世代数適用する工程と、前記遺伝的アルゴリズムから前記染色体とともに得られた解に対する残差を算出する工程と、最小の残差を有する解およびその不明確代替解について傾斜サーチを実行する工程と、前記最小の残差を有する解およびその不明確代替解の比率が、所定値より大きい場合は、前記最小の残差を有する解を前記明確化した飛翔体軌跡として指定する工程とを含む。
【００１１】
さらに本発明の別の局面によると、衝撃波および砲口爆風を検出することにより、砲撃体距離を推定するための方法において、衝撃波単独信号および砲口爆風信号が計測される。初期の砲撃体距離が、初期の弾丸速度および弾丸抗力係数を想定して、計測された衝撃波および砲口爆風信号から推定され、弾丸軌跡に沿った瞬時の弾丸速度が、反復的に算出されて更新された砲撃体距離を得る。砲口爆風検出チャンネルの数は、通常は衝撃波検出チャンネルの数よりも少ない。
【００１２】
有利な実施の形態は、１以上の次の特徴を含むことができる。衝撃波単独信号および砲口爆風信号の間の到着の時間差（ＴＤＯＡ）ならびに到着角度が計算されて、初期の砲撃体距離を決定する。特定数の反復が実行され得、すなわち、引き続いて決定された更新された砲撃体距離の間の関係が収束基準を満足すると、更新された砲撃体距離が、最終の砲撃体距離であるとみなされる。例えば、収束基準は、引き続いて決定された更新された砲撃体距離の間の差が、または引き続いて決定された更新された砲撃体距離の間のパーセント変化が、所定値よりも小さいように選択されることができる。真の解を得るためには、算出された弾丸速度は、常に少なくとも音速になるように設定される。解は、整合性がチェックされる。例えば、弾丸軌跡角度および到着角度が所定値よりも大きいと決定されると、更新された砲撃体距離は無効であるとみなされる。
【００１３】
算出された砲撃体距離が、無効であると決定されても、解は、なお遺伝的アルゴリズム（ｇｅｎｅｔｉｃ ａｌｇｏｒｉｔｈｍ：ＧＡ）を適用することで得ることができる。例えば、所定数の個体を有するＧＡの初期の個体群が規定されることができ、各個体は、弾丸軌跡の推定した砲撃体距離、外れた方位（ＭＡ）、および外れた仰度（ＭＥ）を含む３タプルにより表現される。ＧＡは、所定数の世代に対して実行され、各世代における個体に対する残差が算出される。各世代において、最小残差を有する解が変異なく生存する個体として選択される。最小残差を有する解が、更新された砲撃体距離として選択される。解は、所定数の反復において各３タプルに対して実行されることで、洗練され得て、訂正した砲撃体距離を算出することができる。ここで、各世代における個体に対する残差が、訂正した砲撃体距離とともに算出される。
【００１４】
ＧＡは、クロスオーバおよび変異演算子を含む。クロスオーバ演算子は、世代の個体群からの２つの個体の間で、外れた方位および外れた仰度の少なくとも１つを交換する。一方、変異演算子は、（３タプルの値をランダムに選択した値に置換する）フィールド変異（ｆｉｅｌｄ ｍｕｔａｔｉｏｎ）、増加的変異（（ｉｎｃｒｅｍｅｎｔａｌ ｍｕｔａｔｉｏｎ）３タプルの全フィールドにおける微小変異を含む）、および（世代における個体を変更しないままとする）変異なしを備える。
【００１５】
本発明のさらに別の局面によると、衝撃波信号から、および限定された数の砲口爆風信号から、飛翔体軌跡を明確化するための方法であり、５以上の間隔を空けた音響センサにおける衝撃波単独信号を計測する工程と、多くても４つのセンサ上で砲口爆風信号を計測する工程と、前記衝撃波単独信号から各センサ対に対する到着の時間差（ＴＤＯＡ）情報を決定する工程とを含む。方法は、さらに、各個体が砲撃体方位、砲撃体仰度、外れた方位、および外れた仰度を含む４タプルで表現されている所定数の個体を含む初期個体群を用いて、所定数の世代に対して遺伝的アルゴリズムを実行する工程と、各世代における個体に対する残差を算出し、その残差がＴＤＯＡ衝撃波および砲口爆風信号の組み合わせの最小自乗適合を含んでいる工程と、を含む。最小残差を有する解およびその不明確代替解の比率が、所定値（例えば少なくとも２の値）よりも大きい場合は、最小の算出された残差を有する解は、明確化した飛翔体軌跡として指定される。
【００１６】
本発明のさらに別の局面によると、衝撃波信号の存在下で砲口波からの信号を抽出するための方法は、砲口波がセンサアレイを横断するのに要する時間に対応する時間ウィンドウを規定する工程と、前記衝撃波信号を検出する工程とを含む。衝撃波信号の検出に続いて、ウィンドウが適時に進行され、前記ウィンドウで受容した全エネルギーが進行時間（ａｄｖａｎｃｅ ｔｉｍｅ）の関数として計測される。計測された全エネルギーが砲口信号に関連づけられる。
【００１７】
スプリアス信号が衝撃波波形として解釈されることを防止するため、計測された衝撃波波形の音響エネルギーが、所定の周波数帯（例えば、約７００Ｈｚと１０ｋＨｚとの間の周波数）に亘って所定の閾値を下回ると、飛翔体軌跡は、誤りであるとして排除され得る。あるいは、またはそれに加えて、計測された衝撃波波形が正値を有する時間間隔が、最小時間を下回るか、最大時間を超える場合は、例えば、約７０μｓ未満か約３００μｓより大きい場合は、飛翔体軌跡は、誤りであるとして排除される。
【００１８】
有利な実施の形態において、全エネルギーは、好ましくは、衝撃波エコーにより生じた検出信号の部分を無視することにより、ウィンドウに亘って計測されたエネルギーを積分することにより決定されることができる。有利には、ピーク信号値は、最大の全エネルギーを生成するウィンドウにおいて決定することができる。ピーク信号値が、ウィンドウにおいて計測された全エネルギーよりも所定の率だけ大きい場合は、そのピーク信号は砲口爆風信号に関連していると認識され得る。
【００１９】
本発明の実施の形態は、１以上の次の特徴を含むことができる。アンテナおよび／または音響センサのタイミング誤差分布が、ゲイン変化、サンプリング変化、およびアンテナセンサのセンサ位置変化に関連し得る。明確化のための信頼度レベルは、アンテナのサイズに依存し、従って、比較的小さいアンテナには比較的高い計測精度が要求される。２つの不明確解が存在すると、明確化した飛翔体軌跡は、２つの不明確解に対する残差の比率に基づいて選択される。
【００２０】
別の有利な実施の形態において、各センサ対に対する到着の時間差（ＴＤＯＡ）が、最初に衝撃波に遭遇するセンサを基準センサとして指定する工程と、例えば、前記基準センサにおける衝撃波単独信号の初期部分の振幅が閾値を交差するとき、タイミング回路の第１のラッチを設定する工程と、により決定され得る。第１のラッチは、その他のセンサのそれぞれに対する開始カウンタを作動させる。このカウンタは、その他のセンサのそれぞれが衝撃波に遭遇するまで対応するセンサにおいて動作する。その他のセンサの１つが、例えば、衝撃波単独信号の初期部分に遭遇すると、これは、そのセンサに対する開始カウンタをそのセンサが停止する第２のラッチを設定する。基準センサに対するその他のセンサに対するＴＤＯＡ値が次いで記録される。
【００２１】
本発明の他の特徴と有利性が、次の好適な実施の形態の説明と特許請求の範囲から明確となるであろう。
【図面の簡単な説明】
【００２２】
本発明のこのような特徴および他の特徴および有利性は、添付の図面を参照して次の図示によりさらに充分に理解される。ここで、要素は、同様の参照符号で符号が付され、寸法が表わされているものではない。
【図１】図１は、アンテナとのマッハ円錐交差の断面図を概略的に示す。
【図２】図２は、７方向音響センサを有する代表的なセンサアレイの概略を示す。
【図３】図３は、衝撃波単独軌跡決定における固有の不明確さを概略的に示す。
【図４】図４は、マッハ円錐の曲率を決定するための、到着計測の時間差に対する確率密度を概略的に示す。
【図５】図５は、砲撃体軌跡間で正確に明確化する確率を概略的に示す。
【図６】図６は、相関過程の概略図を示す。
【図７】図７は、砲撃体軌跡間で正確に明確化するために使用される遺伝的アルゴリズムの処理フローを示す。
【図８】図８は、非衝撃波信号に対して区別するための処理フローである。
【図９】図９は、衝撃波ＴＯＡ（到着の時間）モデルの概略図である。
【図１０】図１０は、距離推定のための概略処理フロー図を示す。
【図１１】図１１は、距離推定のための遺伝的アルゴリズムの概略処理フロー図を示す。
【発明を実施するための形態】
【００２３】
発明の詳細な説明
上述の要旨に記載したように、本発明は、砲撃体距離推定および飛翔体軌跡の明確化のためのさまざまな実施の形態で、方法およびシステムを提供する。これらのシステムと方法は、正確な解に要求されるには不十分な数のパラメータが検出されたときや、そのようなパラメータを信頼して検出できないときに特に有用であり有利である。
【００２４】
超音速飛翔体軌跡は、近接する間隔の数個のセンサによって計測される飛翔体衝撃波到着時間からもっぱら推定される。このようなセンサは、アンテナと呼ばれる「微小な」計測容積の至る所に分布している。計測容積は、センサ間隔が２ｍ以下であれば微小とみなされる。一旦、飛翔体の軌跡が認識されると、砲撃体の位置は、軌跡に沿った距離（ｒａｎｇｅ）をさかのぼることなく知られる。アンテナが砲口爆風音の到着時間も得れば、この距離を見いだすことができる。しかし、砲口爆風は、いつも検出可能というわけではないので、正確な衝撃波単独解が軌跡の決定に重要である。
【００２５】
図１を参照して、衝撃波面が拡張してゆく円錐面を考慮する。この円錐面は、弾丸軌跡に一致してその軸を有する。衝撃波面は、マッハ円錐とも呼ばれる。衝撃波単独解を得るために、拡張する円錐面の到着角、曲率の半径、曲率半径の空間的勾配の３つの特性が、５個以上のアンテナセンサから計測された到着時間から決定されるべきである。
【００２６】
アンテナに最初に到着する円錐面起点（ｃｏｎｉｃａｌ ｓｕｒｆａｃｅ−ｇｅｎｅｒａｔｏｒ）の到着角度は、アンテナでの到着角度について弾丸軌跡の２つの可能性のある相対角度（しばしば「不明確」角度と呼ばれる）を決定する。この「不明確」角度については、図３を参照してさらに詳細に説明する。アンテナでの円錐面の曲率半径は、軌跡の距離と方向の両方を決定する。表面生成起点の経路に沿った曲率半径の勾配は、弾丸がどの方向に移動しているかを決定し、それにより、２つの可能性のある方向間で「不明確さ」を除去している。この２つの可能性のある「不明確な」軌跡角度の間で、これらの３つの衝撃波特性を正確に適切に決定するためには、非常に精緻な計測が要求される。例えば、２つの選択的な砲撃体アスペクト角間の確率的誤差は、約１μｓ以下とすべきである。
【００２７】
要求される精度は、図１に示す衝撃波の伝搬特性を考慮することで推定することができる。図２も参照して、アンテナ２０は、Ｎ個のセンサ（Ｎ＝７）を含み、進行する円錐衝撃波の到着時間を決定することが可能である。入ってくる弾丸軌跡は、どこからでも発生することが本質的に期待できるため、アンテナエレメント２３−２８は、有利には、（Ｃｘ０，Ｃｙ０，Ｃｚ０）の中心に位置する１エレメント２２を有して、球面を覆う位置Ｃ（Ｃｘｊ，Ｃｙｊ，Ｃｚｊ）に一様に分布することができるため、均一なセンサ開口が到着角度に無関係に提供されている。基準センサと呼ぶ第1のセンサが進行する円錐面を検出する時刻をｔ０とする。その他のセンサは、進行する円錐面をそれに続くｔｉとする時刻で検出する。進行する円錐面の方向での音の伝搬距離は、各時間差を局所の音速ｃで乗算することにより、すなわちｄｉ＝ｃ・（ｔｉ−ｔｉ）で得られる。もし、計測誤差が無ければ、Ｎ個の点の３次元座標が衝撃波円錐の全パラメータを理想的に決定して、その他の（Ｎ−１）個のセンサで基準センサを通過する円錐面を決定することもできる。しかし、上述のように、到着時間計測とセンサ座標における誤差が、衝撃波円錐に対する誤ったパラメータと、それによる飛翔体の軌跡の誤りをも生じ得る。次に、２つの別の不明確軌跡角度についての正確な決定をするために必要な到着精度の時間差について記載する。
【００２８】
このシステムは、有利には、車輌ノイズ、振動、風ノイズ、およびＥＭＩ等の砲撃体についての非弾道性信号を確実に間違えないようにする特徴を包含する。例えば、センサマストは、ガタガタ音を立てないよう嵌め合い接続の弾性スリーブで車輌（図示せず）に搭載することができる。センサは、背柱の端部に弾性カップリングで装着することができ、背柱の震動から隔絶するため約１Ｈｚの低周波共振を有する。センサの背柱は、アナログ電子回路を有する共通ハブに取り付けることができ、この共通ハブはまた、センサマストに弾性緩衝取り付けで装着してその振動から隔絶することができる。
【００２９】
さらに、衝撃波由来の信号に通常見られる特徴を欠く信号を除外するため、次の決定アルゴリズムを採用することができる。全ての値はパラメータ化され、すなわち、相対的であり、外部から調節することができる。リストに挙げた値は例示目的のみの提供である。
【００３０】
図８を参照して、プロセス８００は、検出された信号が衝撃波に由来するかを決定する。プロセス８００は、工程８０２で開始し、信号が衝撃としてカウントするのに充分な音量イベントか、例えば、ピーク信号値が所与のパラメータ化された、閾値（例えば５００）を超えるか工程８０４でチェックする。これが当てはまれば、工程８００は工程８０６へと続き、ゼロからピーク信号値への急峻な移行があるかをチェックし、このピーク値への移行が、重要な大きさ（例えばピーク信号値の１／１６）を有する別の信号によって先行したものではないことを確証する。
【００３１】
これが当てはまる場合は、工程８００は、工程８０８へと続き、衝撃波の最小と最大の間の時間が、充分に大きい値（例えば２００−４００μｓ）であるかをチェックする。これが当てはまれば、工程８００は、工程８１０へと続き、最小と最大ピーク信号振幅の大きさが、例えば３５％以内に互いに近接するかをチェックする。これが当てはまれば、工程８００は、工程８１２へと続き、最小ピーク信号からゼロへの圧力ピーク移行が急峻であるかを、実質的に工程８０６における同じ基準を使用してチェックする。これが当てはまれば、工程８００は、工程８１４へと続き、最大信号値とゼロ交差との間の時間、およびゼロ交差と最小信号値との間の時間が、例えば約１８０μｓ以内で比較可能であるかをチェックする。全ての工程で肯定的な応答を生じている場合は、工程８００は、信号が衝撃波である可能性があり、その信号は処理される（工程８１６）。逆に、もし６個の決定工程の１つが否定的な答えである場合は、検出された信号は、衝撃波から由来するものではない（工程８１８）。
【００３２】
図１に戻って、飛翔体軌跡がｘ軸と一致すると仮定する。マッハ角はΘ＝ａｒｃｓｉｎ（１／Ｍ）、ここでＭは、発射速度Ｖを音速ｃで除算して定義されるマッハ数である。Ｌは、アンテナの特性長を意味する。アンテナ２０の２つの端部における円錐の曲率半径は、ｒ１、ｒ２である。図の左半分の端面図は、曲率ｒ１の計測法を示す。距離ｄは、ｄ＝ｒ１・ｃｏｓ（Φ）に等しい。角Φは、ｓｉｎ（Φ）＝Ｌ／２ｒ１で規定されるので、微小角Φに対してΦ〜Ｌ／２ｒ１を得る。円錐面を半径ｒ１で二分するアンテナ面上の点間の曲率の時間差計測は、ｄｔ１＝Δｄ／ｃ＝（ｒ１−ｄ）／ｃ〜ｒ１Φ２／２ｃ＝Ｌ２／（８・ｒ１・ｃ）に等しい。ｒ２＝ｒ１−Ｌ・ｓｉｎ（Θ）における曲率の時間差計測がｒ１をｒ２に置換して同じ表現で与えられる。従って、ｄｔ２＝ｄｔ１＋Ｌ３ｓｉｎ（Θ）／８ｒ１２ｃである。
【００３３】
計測誤差にバイアスが無いとすると、すなわち、計測時間差ｄｔ１とｄｔ２が、異なる平均値ｄｔ１とｄｔ２であるが同じ統計的に決定された偏差σを有するランダムに分散した値であるとすると、アレイの２つの端部における平均計測値は、そこでの局所の曲率を正確に決定する。時間差ｄｔ１とｄｔ２に対する計測値の代表的な分布を図４に示す。
【００３４】
端部２でなされたサンプル計測をＸとして示す。端部２における曲率半径（半径ｒ２）は、端部１におけるもの（半径ｒ１）よりも小さい。従って、Ｘより大きい値を有する１の端部でなされた全計測は、端部１における曲率が端部２におけるよりも大きいという正確な決定が得られる。端部２における計測がＸに等しいとき、正確な決定がなされる確率は、
【数１】
で与えられる。
ｘに関して積分し、変数を代入すると、次の正確な決定を生成する確率が得られる。
【数２】
【００３５】
今度は図５を参照して、正確な決定の確率、すなわち明確化（ｄｉｓａｍｂｉｇｕａｔｉｏｎ）に対する信頼レベルを、２つの代表的なアンテナサイズＬ＝１ｍ、およびＬ＝２ｍについて、飛翔体の軌跡とアンテナ２０間の最接近点（ＣＰＡ）ｒに対してプロットしている。音速は、ｃ＝３４０ｍ／ｓと仮定している。大きなアンテナが不明確な衝撃波単独解に対する極めて拡大した範囲を有していることが明白である。大きなＣＰＡ値に対して、アンテナの２つの端部（ｒ１とｒ２）における曲率の差は、区別するには小さすぎるので、正確に決定する確率は、５０％に、すなわち完全に不明確に近づていく。従って、信頼レベルは、アンテナのサイズに、すなわちその径すなわち空間的広がりに依存する。
【００３６】
上述のように、誤差は、タイミングエラーとセンサ座標の不確定性から生じる。センサ座標の不確定性は、衝撃波到着角度の非常に可変的な関数であるバイアス誤差に寄与している。しかし、ランダム到着角度に対しては、センサ座標誤差は、ランダムな時間差誤差に見える。
【００３７】
タイミング誤差は、チャンネルからチャンネルへのゲインと信号強度の変化の両方からも生じる。到着の時間（Ｔｉｍｅｓ）は、センサ出力が所定の閾値Ｖ０に上昇するときに得られる。ゲイン変化ｄｇにより生じるタイミング誤差ｄｔは、そのチャンネルに対する時間当たりの電圧増加に依存する。ここで
【数３】
である。
【００３８】
タイミング誤差は、信号強度が開口に関して変化するときにも生じる。長さＬの開口で距離ｒにおける円柱形音源に対して、開口に亘る最大信号レベル変化は、ｐ０（Ｌ／２ｒ）に等しい。ここで、ｐ０は、開口中心における音圧である。上記のタイミング誤差式は、このタイプの誤差に対しても適用され、相対速度変分ｄｇ／ｇが式Ｌ／（２ｒ）に置換される。振幅誤差は、センサ間でランダムではなく、開口の全体に亘って最大から中心のゼロへと均一に変化する。１０ｍを超える距離では、１ｍ開口に対して、最大振幅係数は、０．０５未満であり、これは０．２のチャンネルゲイン変分パラメータよりも小さいので、振幅誤差による影響は無視することができる。逆に、上述のように、約１０ｍ未満の距離においては、マッハ円錐半径は、１ｍの長さの開口に対して充分に小さいので、計測誤差はあまり重要ではない。
【００３９】
センサの不確定性により生じるタイミング誤差の現実的な推定については、誤差ベクトルの大きさが統計的に独立しており、０と１ｍｍの間に均一に分散しており、誤差角度が統計的に独立であると仮定すると、均一に分散したランダム時間差誤差方程式の標準偏差は、
【数４】
に等しい。１ＭＨｚでサンプリングするシステムに対する二項分布したランダム時間サンプリング誤差の標準偏差は、０．２５μｓに等しい。ゲイン変化によるタイミング誤差は、約０．０２Ｖ／μｓの電圧変化に対応して、約１８ｋＨｚのチャンネル帯域幅を有する代表的なシステムについて約０．７５μｓと推定される。採用された各アレイに対する音響センサは、±１，５ｄＢ以内の感度を有するように選択された。従って、チャンネル相対ゲイン変化は、０．８４と１．１９の間にほぼ均一に分布しているため、相対ゲインの標準偏差は、ほぼ
【数５】
に等しい。閾値電圧は、Ｖ０＝０．１５Ｖであり、タイミング誤差の標準偏差は約０．７５μｓという結果となる。
【００４０】
全体の計測タイミング誤差は、チャンネルゲイン変化分、サンプリング変化分、センサ位置変化分が全て統計的に独立であるという仮定により推定される。次いで、タイミング誤差標準偏差は、
【数６】
として推定される。
【００４１】
このような精度をアナログからデジタルへの変換を伴って達成することは、高いサンプリング率とその後の補完が必要とされるため困難である。到着の時間差（ＴＤＯＡ）を正確に計測するための２つの異なる回路が開示のシステムに採用されている。
【００４２】
１つの実施の形態において、代表的なシステムは、各チャンネルにおいて１ＭＨｚクロックを使用するアナログの到着時間差（ＴＤＯＡ）回路を使用する。このクロックは、センサ信号が基準センサにおける閾値レベルを超えるとトリガされる。基準センサは、上記のように、センサが最初に衝撃波に遭遇するときに規定されたものである。上述のように、１ＭＨｚクロック速度は、実際の時間−サンプル誤差の重さを除去するために充分である。このシステムは、閾値の検出に依存してアナログモードで操作され、次の機能を実行するデジタル論理を有する。
１．衝撃波に最初に遭遇する基準センサにおけるチャンネル信号振幅が閾値を交差するとき第1のラッチが設定される。
２．第1のラッチが各チャンネルに対するスタートカウンタを設定する。これは各クロックサイクルにおいて１カウント増加する。プロセッサが警告を受ける。
３．各チャンネルのカウンタが、対応するセンサが衝撃波に遭遇するまで作動する。これは、そのチャンネルにおけるカウントを停止する第２のラッチを設定する。第２のラッチが設定されない場合は、対応するカウントは上限値まで作動する。
４．各カウンタの最後のカウントの数がデジタルＴＤＯＡレジスタに記録される。
５．プロセッサがＴＤＯＡレジスタを読み取る。
６．ブロセッサは、次の衝撃波を受容するためにカウンタをリセットする。
【００４３】
別の実施の形態において、ハードウェアＴＤＯＡ検出の時間を中央とした時間セグメントに対して、各チャンネルの１つおきのチャンネルとの相関が算出される。Ｃｏｒｒ（ｇ，ｈ）で示す２つの関数の相関は、
【数７】
で定義される。
【００４４】
この相関は、ｔの関数であり、「遅れ（ｌａｇ）」と呼ばれる。従って、時間領域にあり、次の特性を有する。
【数８】
ここで、ｇとｈは、時間の実関数である。Ｇ（ｆ）は、ｇ（ｔ）のフーリエ変換であり、Ｈ（ｆ）は、ｈ（ｔ）のフーリエ変換である。
【００４５】
信号の全パワーは、
【数９】
【００４６】
到着時間信号は、有限の長さを有するので、積分（または離散データに対する総和）は、到着時間のあたりを中心とした有限の時間間隔について実行されればよく、公知のように、１つまたは両方のチャンネルのデータの長さは、ゼロ挿入（ｚｅｒｏ ｐａｄｄｉｎｇ）により拡張され得るので、２つの信号の期間が整合する。
【００４７】
次の議論において、連続関数の積分が簡単のために使用されるが、実際のデータはデジタル化された離散値である。当業者は、容易に積分を総和に置き換えるであろう。
【００４８】
図６を参照して、工程６０において、衝撃波信号時間データｇｉ（ｔ），ｇｊ（ｔ）が、各チャンネルｉ，ｊ、工程６０１，６０２において取得され、時間の関数として記録される。工程６０３、６０４において、チャンネルｉにおける全信号パワーは、次の相関の正規化に対して算出される。
【数１０】
【００４９】
チャンネルｉに対する衝撃波信号時間データｇｉ（ｔ）のフーリエ変換Ｇｉ（ｆ）が計算され、その共役Ｇｉ（−ｆ）が形成される（工程６０５）。同様に、衝撃波信号時間データｇｊ（ｔ）のフーリエ変換Ｇｊ（ｆ）が、その他の全てのチャンネルｊに対して算出される（工程６０６）。その後、相互相関Ｇｉ（−ｆ）・Ｇｊ（ｆ）が、各チャンネルペア（ｉ，ｊ）に対して形成される（工程６０８）。これは、「遅れ」ｔの関数ｆｉ，ｊ（ｔ）である。各チャンネルペアに対するＴＤＯＡは、ｆ（ｔ）が最大値を有する時間ｔｍａｘである（工程６１０）。チャンネルｉとｊの間の相関は、
【数１１】
として定義することができる。
【００５０】
チャンネルｉの残差は、工程６１２に示すように、全てのセンサｊに亘ってセンサｉに対する平均値を計算することで算出される。
【数１２】
最良（すなわち最小）の全残差を有するそのチャンネルに対するＴＤＯＡおよび相関は、次いで「最善」の解として選択される（工程６１４）。
【００５１】
上述のように、チャンネルデータは、通常、例えば４１，６６６．６６サンプル／秒の所定のサンプリングレートの離散時間間隔でサンプリングされる。これは、２４μｓのビン幅に対応し、受信信号に対する時間分解能を反映している。相関処理は、３３３３３３サンプル／秒を取ることで３μｓに対して８の係数を取ることで改善された時間分解能でなされる。
【００５２】
一旦センサ間の様々な到着時間差（ＴＤＯＡ）が衝撃波単独の信号から決定されると、砲撃体方位と仰度、および弾丸軌跡を決定することができる。砲撃体位置、すなわちセンサアレイから砲撃体の距離は、砲口爆風信号が既知であればさらに決定することができる。
【００５３】
アレイの中央を中心とする直交座標系において、すなわち｛（Ｃｘ０，Ｃｙ０，Ｃｚ０）＝（０，０，０）｝、所与のセンサ（Ｃｘｊ，Ｃｙｊ，Ｃｚｊ）（図２参照）における衝撃波の到着時間（ＴＯＡ）は、
【数１３】
で与えられる。
ここで、ｃは音速であり、Ｍはマッハ数である。βは、砲撃体位置と弾丸軌跡の「外れ（ミス）角度」を表わし、これには方位角と仰角の両方が含まれる。直接のヒットであれば、β＝０に対応する。マッハ角度Θは、
【数１４】
で定義される。
【００５４】
上述および図３に示すように、所与の砲撃体位置と弾丸軌跡に対して、所与の一式のセンサにおける衝撃波のＴＯＡがほぼ等しい別の砲撃体位置と弾丸軌跡が存在する。２つの不明確な解は、簡単化したモデルにおいては事実上同一であり、衝撃波は、平面波としてセンサアレイを横断して伝搬するものと思われる。ＴＤＯＡの分解能が、衝撃波の曲率を分解するのに充分高ければ、２つのほぼ同一の解は明確化することができる。本質的に不明確な衝撃波単独ＴＤＯＡ解を図３に示す。
【００５５】
充分に正確なＴＯＡ計測を仮定すると、砲撃体位置と弾丸軌跡の真の解は、計測および算出された衝撃波ＴＤＯＡの残差の二乗平均平方根（ＲＭＳ）を最小化する砲撃体／軌道の組み合わせを算出することで得ることができる。
【数１５】
ここで総計は全センサに関してなされる。
【００５６】
この問題を解く１つのアプローチは、米国特許第５，９３０，２０２号にその詳細が記載のＬ１Ｌｅｖｅｎｂｅｒｇ−Ｍａｒｑｕａｒｄｔアルゴリズムである。通常の古典的逐一アルゴリズムは、最適解にアプローチするために決定論的処理を使用し、当て推量解から開始して事前に規定された移行則に基づいてサーチ方向を指定する。このようなものには、目的関数と制約値を使用した直接法、ならびに１階および２階微分を使用した勾配法がある。しかし、このような方法は、例えば、最適解が選択された初期解に依存し、アルゴリズムが、局所最小値や費用関数面が平坦な谷を有する所等の準最適解で「スタック」するという欠点を有する。このためさらに繰り返しても結果が改善しない。
【００５７】
進化型（ｅｖｏｌｕｔｉｏｎａy)遺伝的アルゴリズム（ＧＡ）を使用することで、砲撃体方向と飛翔体軌跡の全体の最小値を、より速くより信頼度を高く明確化して計算することができることがわかった。ＧＡは、自然の進化原理を模倣し、サーチや最適化処理に適用される。
【００５８】
ＧＡの概略フロー図を図７に示す。解に対して単一の推量から開始する代わりに、ＧＡ工程７０は、解のランダムな個体群を初期化すること（工程７１）でそのサーチを開始し、初期解集合を示すゼロに世代カウンタを設定する（工程７２）。一旦、解のランダム個体群が生成されると、それぞれは非線形計画法問題のコンテキストにおいて評価され（工程７３）、適性（優劣）が各解に割り当てられる（工程７４）。この適性は、計算解と計測解の間のユークリッド距離Δτｍｉｎによって表わすことができる。
【数１６】
【００５９】
直感的には、小さい値のΔτｍｉｎを有するアルゴリズムが良い。
【００６０】
例えば、ＧＡを砲撃体方向と飛翔体軌跡に対する解を明確化するために適用するときは、代表的なＧＡは、２００ ４−ｓの初期個体群を染色体として使用する。ここで、各４−は、次の値を有する。
［方位砲撃体、仰度砲撃体、方位外れ、仰度外れ］，［方位砲撃体，仰度砲撃体］は、角度（Θ＋β）により規定され、ここで、［方位外れ、仰角外れ］は、角度β（図３参照）で規定される。砲口爆風は、先述の衝撃波単独アプローチを使用しないので、センサアレイと砲撃体間の１００ｍの基準距離を仮定する。
【００６１】
初期個体群は、重要で合理的な範囲に亘る値（全ての値は角度である。）である４−ｓのランダムな選択により生成され、
方位砲撃体＝｛０，．．．．，３６０｝，
仰角砲撃体｛−１０，．．．．，３０｝，
方位外れ｛−２０，．．．．，２０｝，および
仰角外れ｛−２０，．．．．，２０｝である。
【００６２】
工程７５において、例えば２５に設定可能であるＧＡの反復の最大数に達したかがチェックされる。反復の最大数に達していた場合は、工程７０は、工程８０で停止し、結果が許容できるか別に評価されるかのいずれかであり得る。そうでない場合は、工程７６は、事前に設定した適性基準が満足されたかどうかをチェックする。
【００６３】
適性基準は、例えば、計算された＜１５°の外れた方位および／または２つの不明確解の残差の割合であり得る。適性基準が満足されている場合は、処理７０は工程８０で停止し、そうでなければ、新たな個体群が、クロスオーバ（工程７７）、変異（工程７８）を通じて作成され、生成カウンタが1刻みで増加する（工程７９）。
【００６４】
各世代において、「最良の」個体は、変異することなく生存を許され、それに対して、トップ１００の個体も、その適性により判断されて生存するが、表１に列挙したクロスオーバ／変異演算子を伴って、これらの生存体のペアから次の１００の個体を作成するために使用される。
【００６５】
次の代表的なクロスオーバおよび変異演算子が、工程７０を実行するために使用された。
【表１】
【００６６】
明確化は、最良解と対応する代替解の勾配サーチ（ｇｒａｄｉｅｎｔ ｓｅａｒｃｈ）を実行することで達成および／または改善される。両方の不明確な解に対して、残差と残差の比率が算出される。算出された外れた方位が＜１５°の場合、「近い」砲撃であり、残差の比率が＞２の場合は、低い方の残差を有する解が選択される。さもなくば、実際の選択がなされず、低い方の残差を有する解が「主要」解と標示され、他方の解は「代替」解と標示される。
【００６７】
衝撃波単独検出に際して、ＧＡアルゴリズムは、Ｌｉｎｕｘオペレーティグシステムを作動する１ＧＨｚコンピュータで０．１５秒で広範なシミュレートされた砲撃において１つの解を生成した。９７％のシミュレートされた砲撃が外れた方位の１５°以内で、シミュレートされた砲撃の８６％が外れた方位の５°以内であった。先述の明確化アルゴリズムを使用して、近い砲撃、すなわち１５°未満の外れた方位を有する砲撃は、その時の９５％が明確化された。明確化アルゴリズムは、その時の７０％の遠い砲撃に対して正確な結果を生成した。明確化の正確性さは、センサアレイ配置と推定された砲撃の分布に基づいて変化すると予測され、低い仰度を有する砲撃は明確化が容易である。
【００６８】
飛翔体軌跡に対する先述の解は、砲口爆風を検出することなく得られた。しかし、弱い砲口信号または限定された数のチャンネル上のみにおいて受信された砲口信号でさえ、距離決定と明確化の改善のために有利に使用することができる。
【００６９】
図９は、到着時間（ＴＯＡ）モデルを概略的に示す図であり、米国特許第６，１７８，１４１号にさらに詳細が記載されている。ＴＯＡモデルは、センサ位置に対する飛翔体の軌跡と砲撃体方向を推定するために使用することができる。ＴＯＡモデルは、飛翔体の飛程に関する特定の物理特性を考慮した弾道モデルに基づいている。このような物理特性には、空気密度（温度に関連）、砲撃体の位置
【数１７】
ライフル砲口の方位角および仰角、飛翔体の砲口速度（または相当のマッハ数）、および音速（温度／空気密度により変化）等がある。この弾道モデルを使用して、衝撃波と砲口爆風が空間の特定の点に到達する時間を正確に計算することができる。
【００７０】
図９に示すように
砲撃体は、原点（０，０，０）に対して点
【数１８】
に位置し、多様なセンサが点
【数１９】
に配置され、弾丸軌跡は、
【数２０】
の方向に砲撃体から生じているとして示される。砲撃体からＪ番目のセンサへのベクトルは、
【数２１】
であり、ｊ番目のセンサへの弾丸の最接近点（ＣＰＡ）は、
【数２２】
であり、衝撃波が軌跡からｊ番目のセンサへ生じる点からの経路は、
【数２３】
（センサの添字ｊは省略）である。
弾丸のマッハ角は、θ＝ｓｉｎ−１（１／Ｍ），Ｍ＝Ｖ／ｃ０である。Ｍは飛翔体のマッハ数であり、Ｖは飛翔体の超音速であり、ｃ０は、音速（圧力と温度に依存）である。軌跡とｊ番目のセンサの間の「外れ角度」はβである。軌跡は、ｘ−ｙ平面においてＸ軸から反時計回りに計測された方位角によって、およびｘ−ｙ平面から上方に計測された仰角によって特徴づけられる。衝撃波到着時間ｔｊとｊ番目のセンサにおける単位ベクトルを規定する方程式は、これらの幾何学的量という形で記述される。
【００７１】
到着の時間は、飛翔体が、ｊ番目のセンサに向けて音が放射される点までの距離
【数２４】
を進むのに要する時間
【数２５】
に、衝撃波が放射点からｊ番目のセンサまでの距離
【数２６】
を進むのに要する時間
【数２７】
を加えたものに等しい。
【数２８】
ここで、ｔ０は、時間基準（発射時間）であり、ｃ０は、音速である。マッハ角Θも、図９に示す。
【００７２】
飛翔体の速度Ｖは、センサ間隔に対応する距離に亘って一定のままであると安全に仮定することができるので、飛翔体が異なるセンサへと放射する時間の間では大きな速度の損失は存在しない。しかし、比較的長距離に亘っては、飛翔体は空気抵抗により減速することが知られている。空気抵抗は、抗力係数Ｃｂで表現することができ、これは弾丸形状と弾丸口径に依存する。物理的原理から得られる数学的弾道モデルは、空中における任意の点での衝撃波の到着時間を、あらかじめ既知の、飛翔体を記述する全パラメータ一式（例えばその抗力係数Ｃｂにより）、その初期速度、および周囲空気の密度の関数として予測することができる。
【００７３】
正確な計算に必要なパラメータは、戦場などの実際的な設定においては通常は既知ではない。しかし、距離推量は、図１０の処理フロー図２００の形態に示す処理を繰り返すことにより、かなり改善することができる。この処理においては、軌跡に沿った飛翔体の速度の減速を考慮している。処理２００は、次の仮定、
ｃ０＝外気温度／空気圧に対して修正された音速（〜３４０ｍ／ｓ）
Ｃｂ＝予測の武器に関して平均した公称抗力係数
Ｖ０＝予測の武器に関して平均した、飛翔体の発射時の初期速度
Ｍ０＝Ｖ０／ｃ＝飛翔体の初期マッハ数
【００７４】
砲撃体距離Ｄ０の第1の推定は、方程式
【数２９】
に従って、計測された到着の時間差（ＴＤＯＡ）τｍｓ、およびセンサアレイにおける衝撃と砲口音の間の到着角度αを使用して、ならびに初期、定常速度Ｖ０、およびマッハ数Ｍ０を仮定することにより、工程２０４において計算される。
【００７５】
これらの仮定とともに、砲撃体位置
【数３０】
からの距離における飛翔体の速度が、工程２０６において方程式
【数３１】
から、計算され得るので、飛翔体が軌跡に沿った距離ａを移動する時間（工程２０８）は、
【数３２】
となる。
【００７６】
角度Θは、方程式
【数３３】
により、マッハ数Ｍａに関連する。ここで、マッハ数Ｍａは、初期的にはＭ０に設定される。瞬時の弾丸速度は、算出された弾丸速度が音速よりも小さくなる場合は音速（すなわちＭａ＝１）に設定されることに注意すべきである。工程２１０における訂正した距離
【数３４】
は、
【数３５】
となる。
【００７７】
角度α、β、およびΘは、方程式（α＋β＋Θ）＝９０°という関係にある。処理２００は、次いで上記のＭａとＴａについての方程式において算出した距離ａを挿入することにより工程２０６へループバックし、移動距離ａに対して更新したマッハ数Ｍａと更新した弾丸移動時間Ｔａをそれぞれ得る。計測されたＴＤＯＡτｍｓとＴａとａに対する算出した更新値が次いで砲撃体距離に対する値Ｄを連続的に更新するために使用される。
【数３６】
【００７８】
この処理は、工程２１２で決定するように、反復回数の最大数が達したか距離値Ｄが収束するかのいずれかまで繰り返される。
【００７９】
処理２００は、砲撃体とセンサアレイ間の距離に対する訂正した距離値
【数３７】
が「合理的な」値であるかも工程２１４でチェックする。該当する場合は、処理２００は工程２１６で終了する。例えば、飛翔体が移動する距離ａとセンサと音波が飛翔体からセンサへ放射される点までの間の距離
【数３８】
とが有効な数値である場合、すなわちＮＡＮ（非数）でない場合は、距離Ｄに対する値は有効であると考えることができる。ＮＡＮは、数値演算の結果を有効な数値に戻すことができないことを示す特別な浮動小数点値であり、通常は、計算中のエラーの伝搬を阻止するために使用される。加えて、αとβは両方とも所定の閾値未満である必要があり、これは、飛翔体が確かにセンサアレイに向けて発射されたことを示すものである。
【００８０】
上述のように、数値対（τｍｓ、α）が零次近似の砲撃体距離Ｄ０の算出に初期的に使用される。これにおいては、軌跡に沿った飛翔体の速度の変化が無視される。上述の反復処理２００が数値対（τｍｓ、α）を支持する一致構造を返さない場合は、その解は廃棄される。
【００８１】
正確な解が得られない場合でも、計測衝撃ＴＤＯＡと計測砲口ＴＤＯＡに最も近く整合する、砲撃体距離Ｄ、外れた方位角と仰角（これらはβに関連する。）に対する値を見つけることが目標である。既述のように、多様なセンサ間の衝撃単独ＴＤＯＡが、通常の状況において容易に計測することができる。砲撃体方位と砲撃体仰度は、既知のセンサアレイ座標（Ｓｘｊ，Ｓｙｊ，Ｓｚｊ）を使用して衝撃単独ＴＤＯＡから決定することができるが、砲撃体距離についてはできない。検出された衝撃波と砲口音の間のＴＤＯＡτｍｓも計測され得るが、砲口音はそれによっては全てのセンサにより検出することができないことが考えられる。
【００８２】
反復処理２００が有効な結果を返さないと工程２１４において決定されると、次いで処理２００は、進化型遺伝的アルゴリズム（ＧＡ）３００を呼び出すことによって砲撃体距離を算出する。ＧＡは、自然の進化原理を模倣し、サーチや最適化手順に適用される。ＧＡは、そのサーチを、１つのみの解ではなく、ランダムな解の集合から開始する。一旦、ランダムな解の個体群が生成されると、各々は、非線形計画問題のコンテキストにおいて評価され、適性（優劣）が各解に対して割り当てられる。１つの実施の形態において、この適性は、計算解と計測解間のユークリッド距離で、例えば、
【数３９】
で表現することができる。
【００８３】
直感的には、比較的小さい値τｍｉｎを生成するアルゴリズムが良い。
【００８４】
ＧＡ処理３００の概略フロー図を図１１に示す。処理３００は、工程２００に対して先に計測された到着の時間差（ＴＤＯＡ）τｍｓと角度αを使用する（工程３０２）。値｛ＲＡＮＧＥ（距離）、ＭＡ、ＭＥ｝を有する３タプル（３−ｔｕｐｅｓ）の典型的な数が初期個体群として規定される（工程３０４）。ここでＲＡＮＧＥは、図９に示す砲撃体距離
【数４０】
であり、ＭＡは、外れた方位であり、ＭＥは外れた仰度である。ＭＡとＭＥ値は、方位と仰度空間において弾丸が標的を外れた程度を示す。図示の例の標的は、センサアレイと仮定している。工程３０４における初期個体群は、ランダムな選択の３タプル（３−ｔｕｐｅｌｓ）により生成され、この中には重要で合理的な値の範囲
距離砲撃体（Ｒａｎｇｅｓｈｏｏｔｅｒ）＝｛１０００，．．．，３０００｝［メートル］
方位外れた（ＡｚｉｍｕｔｈＭｉｓｓｅｄ）＝｛−２０，．．．，２０｝［度］、および
仰度外れた（ＥｌｅｖａｔｉｏｎＭｉｓｓｅｄ）＝｛−２０，．．．，２０｝［度］
が含まれている。
【００８５】
算出は、衝撃波単独解に対する上記概説したものと同様の処理に従う。初期的には、世代Ｇｅｎ＝０に対して（工程３０６）、砲撃体位置ベクトル
【数４１】
が、先に決定された砲撃体方位と仰度、および特定の３タプルに対する仮定した距離（ＲＡＮＧＥ）を使用して、各３タプルに対して算出される。初期のマッハ数をＭ０とすると、ベクトル
【数４２】
すなわち、衝撃音が放射する位置が、各３タプルに対するＭＡとＭＥ値を用いて計算される（工程３０８）。砲撃体と衝撃波を検出した各センサｊ間の距離
【数４３】
も算出される。
【００８６】
各３タプルに対して、角度βが方程式
【数４４】
から算出される。ここでシンボル「・」は、２つのベクトル間のスカラー積を示す。βに対して算出した値、およびＭａ＝Ｍ０に対して初期的に仮定した値、およびａを、Ｍａ、Ｔａ、ａおよびＤに対する上記の方程式に挿入することにより、距離ａに対して更新された値、飛翔体の距離ａに亘る移動時間Ｔａ、およびマッハ数が算出される（工程３１２）。この処理は、各３タプルに対して数回、例えば工程３１２で決めているように３回繰り返され、その後、砲口信号を含む上記で規定された残差Δτｍｉｎが各３タプルに対して算出される（工程３１４）。
【００８７】
工程３１６において、ＧＡに対する反復の最大数（例えば２５反復）が達成されたかがチェックされる。反復の最大数に達した場合は、処理３００は工程３２０で停止し、最小残差を有する３タプルを返す。達していない場合は、処理３００は、新たな個体群をクロスオーバと変異演算から作成し（工程３１８）、世代カウンタが１増加する（工程３２２）。
【００８８】
各世代において、「最良の」個体は、変異することなく生存を許され、それに対して、トップ１００の個体も、その適性により判断されて生存するが、以下の表２に列挙したクロスオーバ／変異演算子を使って、これらの生存体のペアから次の１００の個体を作成するために使用される。
【００８９】
次の例示のクロスオーバと変異演算子が処理３００を実行するために使用された。
【表２】
【００９０】
ＧＡ処理３００が、２００の異なる３タプルの初期個体群で、５０の補充率（ｒｅｆｉｌｌ ｒａｔｅ）で、合計２５世代に対して実行される。ＧＡは、異なる初期の３タプルのセットと並行して５回作動し、最小の残差を有する解が、砲撃体の距離（ＲＡＮＧＥ）、外れた方位、および外れた仰度に対する最終的解として選択され、これがベクトル
【数４５】
の算出を可能とする。
【００９１】
最近の実験的施行によると、５以上の衝撃波チャンネルに加えて少なくとも１つの砲口信号チャンネルを使用することにより、同じデータセット上で、不明確な砲撃が９５％から８％へと減少を示した。これは、衝撃波単独解に関して重要な改善である。
【００９２】
空気抵抗による飛翔体の経路に沿った減速を考慮しない計算は、過大評価した距離となる傾向がある。特定の幾何的配置と充分に遠い砲撃に対して、この過大評価は２０％を超えることもある。先述の処理は、長距離砲撃検出に対する距離推定からこのバイアスを除去する。
【００９３】
上述のように、２つの異なる軌跡から残差を比較し、小さい方の残差を有する軌跡を選択することで、不明確な衝撃波単独解は、明確化され得ることがしばしばである。
【００９４】
砲口爆風信号は、４以上のセンサチャンネルで検出されると、衝撃チャンネルの数に関係なく、先述の衝撃砲口アルゴリズムを使用して、明確に砲撃体位置を決定することができる。砲口爆風信号が４未満のセンサで検出されるが、衝撃波信号が５以上の衝撃波チャンネルで検出されると、修正した費用関数または残差を有する先述のＧＡを使用して、どのような砲口信号が得られても最適化関数に「混合」されて衝撃波単独解を明確化し、および／または砲撃体距離の推定精度を上げる。しかし、３未満の砲口チャンネルおよび５未満の衝撃波チャンネルが検出されると、砲撃体位置を特定しようとすることなく警告が起動する。
【００９５】
砲口信号は、全てのチャンネルで確実に検出されるわけではない。その理由は、
１．１以上のチャンネルでの検出レベルでは、確信的に検出するには低すぎる。
２．砲口エネルギーが生信号中では区別できず、システムに関連するノイズを生じさせ、ＴＤＯＡ推量の信頼度を低くする。
３．衝撃波からのエコーが砲口爆風よりも強い場合があり、砲口爆風よりも早く到着することがあり、システムが衝撃を砲口として誤検出することがある。
【００９６】
数チャンネルで検出されたのみの砲口爆風信号に関して、この状況での残差は、
【数４６】
として定義することができる。ここで、砲口爆風の最初の項は、砲口爆風を検出する減少した数のセンサ（＜４）について総計されており、ｊは、衝撃波を検出するセンサ（通常は全センサ）に関して総計されている。
【００９７】
上述の例に示すように、砲口爆風信号は、衝撃波単独解に比べて、砲撃体方位とそれによる飛翔体の軌跡についての重要な情報を提供するので、算出された軌跡の解は、不明確な解の１つとより近くに整合する。すなわち、このようにして解を明確化している。
【００９８】
少なくともいくつかの信頼できる砲口信号が無ければ、特に長い砲撃体距離においてかなりの数の不明確な衝撃波単独解が生成され得るがあまり正確な解ではなく、これは比較的少ない数の明確なものほど望ましいものではない。
【００９９】
信頼度が無いかもしれない砲口検出の場合は、砲口爆風信号を検出する試みが最初はなされ得る。例えば、ノイズの多い信号の中の砲口爆風痕跡を見いだして、得られるＴＤＯＡを算出することである。チャンネル間で充分な相互相関を有して充分な数のセンサで砲口信号が見い出される場合であって、砲口信号と各チャンネルの対応する未処理のバンド（フィルタ遅延を補償するために数ビンだけオフセットされている。）との間に充分に強い相関がある場合は、砲口検出は信頼度があるものと考えられる。
【０１００】
そうでない場合は、不十分な相関を示す砲口信号が少なくとも消去され、次の「粗い砲口検出」ロジックが呼び出される。
−衝撃に続く衝撃エネルギーのピークを探す。これらのピークを「衝撃エコー」の可能性ありとしてフラグし、これによって砲口爆風としては排除する。
−砲口波がセンサアレイを横断するのに要するであろう時間を決定し、対応する所用時間を有する「ウィンドウ」を定義する。検出された衝撃波に続いて、衝撃エコーとして認識された検出信号の部分を無視しつつ、実質的に全部の検出器チャンネルに亘ってこのウィンドウを動かすことにより砲口エネルギーピークを探す。ウィンドウに関してエネルギーを積分する。すなわち、
【数４７】
の最大値を探す。
ここで、
【数４８】
の２乗は、ｎ番目センサによって計測された砲口爆風のエネルギー等のエネルギーの測度を示し、ｉ＋ｊは、検出チャンネルを示し、ｉは、衝撃波が検出された時間とウィンドウの開始との間の離散時間間隔を示し、ｊは、ウィンドウの開始からの計測された時間間隔を示す。
【０１０１】
ノイズに対して区別するため、関数ｆｍａｘ（ｉ）の最大値を生成するウィンドのピークエネルギーがチェックされ、その最大値におけるエネルギーピークが、ウィンドウに亘るエネルギーよりも所与の比率係数だけ大きいかを決定する。該当する場合は、ウィンドウの信号が砲口検出として認識され、砲口爆風バンドｍｂの全チャンネルで相互相関が実行されて砲口ＴＤＯＡを決定する。
【０１０２】
検出された砲口爆風信号は、次いで上述のように、砲撃体距離の決定および／または衝撃波信号の明確化のために使用される。
【０１０３】
要するに、記載のシステムは、衝撃波単独計測に基づいて、正確に、迅速に、そしてしばしば明確に、砲撃体方向と弾丸軌跡を提供する。明確化は、改善することができ、砲撃体距離は、微弱な砲口爆風波形が検出される場合でも推定することができる。システムは、車輌震動やノイズ、風音にはそうではないが、花火やシステムから遠ざかる方向の近くの発砲、に応答する誤った砲撃体指示に比較的鈍感である。
【０１０４】
衝撃波信号を検出するシステムは、信号が本当に衝撃波に由来するかを決定するため、初期波形上で２つのテストを実施する点が記載されるべきである。第１に、約７００Ｈｚと１０ｋＨＺの間の周波数帯における計測された全エネルギーを経験的閾値と比較する。この閾値を超える場合に限り、信号形態は衝撃波から生じているものと考えることができる。第２に、検出された初期正圧ピークの時間間隔が、約７０μｓより大きく約３００μｓより小さくなければならない。これらの基準は、花火や威嚇性のない発砲などのインパルス的ノイズからシステムの免疫を提供する。これらのテストを通過しない場合は、検出された波形は、衝撃波ではないと考えられ、砲撃体解は試行されない。
【０１０５】
本発明は、図示および詳述な記載の好適な実施の形態に関連して開示してきたが、これに関する多様な改変および改善が本発明の精神と範囲から逸脱することなくなされ得る。例えば、例示の実施の形態においては、マイクロフォン等の音響センサを有するとして図示しているが、これはそうである必要はない。代わりに、別のタイプの圧力感応性機構や電気的センサが使用され得る。さらに、表１と２の多様な演算子のための値は、例示のみのための意図であり、実際の現場の状況に応じて別の値が選択されることがあり得る。従って、本発明の精神と範囲は、添付の特許請求の範囲のみによって限定されるべきものである。
【特許請求の範囲】
【請求項１】
衝撃波単独信号から飛翔体軌跡を明確化するための方法であって、以下、
アンテナを形成する５以上の間隔を空けた音響センサにおいて、衝撃波単独信号の少なくとも初期部分を計測する工程と、
前記音響センサに対するタイミングエラー分布を推定する工程と、
前記衝撃波単独信号の計測された初期部分から、センサ対に対する到着の時間差（ＴＤＯＡ）を、前記推定されたタイミングエラー分布より大きい時間分解能で決定する工程と、
明確化に対する規定の信頼度レベルに基づいて、および前記音響センサのＴＤＯＡに対する残差の値に基づいて、前記明確化した飛翔体軌跡を選択する工程と、
を備える、方法。
【請求項２】
前記タイミングエラー分布が、ゲイン変動、サンプリング変動および／またはアンテナセンサのセンサ位置変動に関連している、請求項１に記載の方法。
【請求項３】
前記明確化に対する信頼度レベルが前記アンテナのサイズに依存する、請求項１または２に記載の方法。
【請求項４】
前記明確化した飛翔体軌跡が、２つの不明確解に対する比率に基づいて選択される、請求項１〜３のいずれかに記載の方法。
【請求項５】
センサ対に対する前記到着の時間差（ＴＤＯＡ）が、
前記衝撃波に最初に遭遇するセンサを基準センサとして指定する工程と、
前記基準センサにおける衝撃波単独信号の初期部分の振幅が閾値を交差するとき、タイミング回路の第1のラッチを設定する工程であって、前記第1のラッチが、その他の各センサに対する開始カウンタを始動させ、前記対応するセンサが前記衝撃波に遭遇するまで前記その他の各センサのカウンタが動作する工程と、
前記その他のセンサの１つが、前記衝撃波単独信号の初期部分に遭遇するとき、前記その他のセンサに対する開始カウンタを停止する前記対応する他のセンサに対する第２のラッチを設定する工程と、
前記基準センサに対する前記その他のセンサに対するＴＤＯＡ値を記録する工程と、
により決定される、請求項１〜４のいずれかに記載の方法。
【請求項６】
前記計測された衝撃波形の音響エネルギーが、所定の周波数帯に亘って所定の閾値を下回るとき、飛翔体軌跡が誤りであるとして除去される、請求項１〜５のいずれかに記載の方法。
【請求項７】
前記所定の周波数帯が、約７００Ｈｚと１０ｋＨｚとの間の周波数を備える、請求項６に記載の方法。
【請求項８】
計測された衝撃波波形が正値を有する時間間隔が、最小時間を下回るか、最大時間より大きいとき、飛翔体軌跡が誤りであるとして除去される。請求項１〜７のいずれかに記載の方法。
【請求項９】
前記最小時間が、約７０μｓであり、前記最大時間が約３００μｓである、請求項８に記載の方法。
【請求項１０】
前記明確化飛翔体軌跡が、他の算出されたいずれの飛翔体軌跡よりも小さい値を有する、請求項１〜９のいずれかに記載の方法。
【請求項１】
衝撃波単独信号から飛翔体軌跡を明確化するための方法であって、以下、
アンテナを形成する５以上の間隔を空けた音響センサにおいて、衝撃波単独信号の少なくとも初期部分を計測する工程と、
前記音響センサに対するタイミングエラー分布を推定する工程と、
前記衝撃波単独信号の計測された初期部分から、センサ対に対する到着の時間差（ＴＤＯＡ）を、前記推定されたタイミングエラー分布より大きい時間分解能で決定する工程と、
明確化に対する規定の信頼度レベルに基づいて、および前記音響センサのＴＤＯＡに対する残差の値に基づいて、前記明確化した飛翔体軌跡を選択する工程と、
を備える、方法。
【請求項２】
前記タイミングエラー分布が、ゲイン変動、サンプリング変動および／またはアンテナセンサのセンサ位置変動に関連している、請求項１に記載の方法。
【請求項３】
前記明確化に対する信頼度レベルが前記アンテナのサイズに依存する、請求項１または２に記載の方法。
【請求項４】
前記明確化した飛翔体軌跡が、２つの不明確解に対する比率に基づいて選択される、請求項１〜３のいずれかに記載の方法。
【請求項５】
センサ対に対する前記到着の時間差（ＴＤＯＡ）が、
前記衝撃波に最初に遭遇するセンサを基準センサとして指定する工程と、
前記基準センサにおける衝撃波単独信号の初期部分の振幅が閾値を交差するとき、タイミング回路の第1のラッチを設定する工程であって、前記第1のラッチが、その他の各センサに対する開始カウンタを始動させ、前記対応するセンサが前記衝撃波に遭遇するまで前記その他の各センサのカウンタが動作する工程と、
前記その他のセンサの１つが、前記衝撃波単独信号の初期部分に遭遇するとき、前記その他のセンサに対する開始カウンタを停止する前記対応する他のセンサに対する第２のラッチを設定する工程と、
前記基準センサに対する前記その他のセンサに対するＴＤＯＡ値を記録する工程と、
により決定される、請求項１〜４のいずれかに記載の方法。
【請求項６】
前記計測された衝撃波形の音響エネルギーが、所定の周波数帯に亘って所定の閾値を下回るとき、飛翔体軌跡が誤りであるとして除去される、請求項１〜５のいずれかに記載の方法。
【請求項７】
前記所定の周波数帯が、約７００Ｈｚと１０ｋＨｚとの間の周波数を備える、請求項６に記載の方法。
【請求項８】
計測された衝撃波波形が正値を有する時間間隔が、最小時間を下回るか、最大時間より大きいとき、飛翔体軌跡が誤りであるとして除去される。請求項１〜７のいずれかに記載の方法。
【請求項９】
前記最小時間が、約７０μｓであり、前記最大時間が約３００μｓである、請求項８に記載の方法。
【請求項１０】
前記明確化飛翔体軌跡が、他の算出されたいずれの飛翔体軌跡よりも小さい値を有する、請求項１〜９のいずれかに記載の方法。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【公開番号】特開２０１１−５９１２８（Ｐ２０１１−５９１２８Ａ）
【公開日】平成２３年３月２４日（２０１１．３．２４）
【国際特許分類】
【出願番号】特願２０１０−２４８３６６（Ｐ２０１０−２４８３６６）
【出願日】平成２２年１１月５日（２０１０．１１．５）
【分割の表示】特願２００７−５３００９０（Ｐ２００７−５３００９０）の分割
【原出願日】平成１７年８月２４日（２００５．８．２４）
【公序良俗違反の表示】
（特許庁注：以下のものは登録商標）
１．Ｌｉｎｕｘ
【出願人】（５０７３４８０１２）レイセオン ビービーエヌ テクノロジーズ コーポレーション (6)
【Ｆターム（参考）】
【公開日】平成２３年３月２４日（２０１１．３．２４）
【国際特許分類】
【出願日】平成２２年１１月５日（２０１０．１１．５）
【分割の表示】特願２００７−５３００９０（Ｐ２００７−５３００９０）の分割
【原出願日】平成１７年８月２４日（２００５．８．２４）
【公序良俗違反の表示】
（特許庁注：以下のものは登録商標）
１．Ｌｉｎｕｘ
【出願人】（５０７３４８０１２）レイセオン ビービーエヌ テクノロジーズ コーポレーション (6)
【Ｆターム（参考）】
[ Back to top ]