(57)【要約】【課題】移動体の航跡に関するオペレータの先見情報の利用を可能とする移動体の航跡分析装置及び方法を提供する。【解決手段】移動体としての船舶の航跡の分析処理を実行する分析処理器４と、オペレータの操作を受け付ける端末部８とを備える。分析処理器４は、受動センサ１、観測装置２、観測データベース３を解して入力される観測結果から船舶航跡分析処理を実行する。分析処理器４は、この処理の実行中に、端末部８で受け付けたオペレータの船舶航跡の推測値や初期条件変更の判断に関する先見情報を分析処理に反映させる。これにより、より高精度な船舶航跡分析を短時間に実現することができる。【選択図】      図１
2004239678.tif 000002
【特許請求の範囲】
【請求項１】
移動体から放射される音波または電波を検出するセンサにより時系列に測定した観測データを用い、前記移動体の航跡を分析する移動体の航跡分析装置において、前記時系列に測定した観測データに基づいて前記移動体の航跡の分析処理を実行する分析処理部と、オペレータの操作を受け付ける端末部とを備え、前記端末部は、オペレータからの航跡に関する先見情報を取り入れ、前記分析処理部は、前記端末部で受け付けたオペレータからの先見情報を反映させて前記移動体の航跡の分析を行うことを特徴とする移動体の航跡分析装置。
【請求項２】
移動体から放射される音波または電波を検出するセンサにより時系列に測定した観測データを用い、前記移動体の航跡を分析する移動体の航跡分析装置において、前記観測データを格納する記憶部と、前記時系列に測定した観測データに基づいて前記移動体の航跡の分析処理を実行する分析処理部と、オペレータの操作を受け付ける端末部とを備え、前記端末部は、オペレータからの初期値変更の指示を取り入れ、前記分析処理部は、前記端末部で受け付けたオペレータからの初期値変更の指示に応じて、移動体の航跡の分析における初期値を変更し、前記記憶部に記憶されている観測データを用いて再度移動体の航跡の分析を行うことを特徴とする移動体の航跡分析装置。
【請求項３】
前記端末部は、前記先見情報の入力を受け付けることに加えて、当該先見情報入力後も継続して先見情報を入力するか、当該先見情報に基づき推定初期値を変更して再度分析を行うかを選択するための選択指示を受け付けることを特徴とする請求項１記載の移動体の航跡分析装置。
【請求項４】
前記移動体の航跡の分析結果を表示するための表示部と、前記表示部が表示している分析結果が、前記端末部が受け付けた前記先見情報を用いて得られたものである場合に、それをユーザに告知するための告知手段とをさらに備えることを特徴とする請求項１記載の移動体の航跡分析装置。
【請求項５】
前記表示部は、入力された先見情報の妥当性を評価するための時系列残差をグラフィカルに表示することを特徴とする請求項１に記載の移動体の航跡分析装置。
【請求項６】
移動体から放射される音波または電波を検出するセンサにより時系列に測定した観測データを用い、前記移動体の航跡を分析する移動体の航跡分析方法において、オペレータからの移動体の航跡または分析条件に関する先見情報を受け入れ、前記受け入れた先見情報を用いて前記移動体の航跡を分析することを特徴とする移動体の航跡分析方法。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は、船舶、航空機等の移動体の航跡分析装置及び分析方法に係り、移動体から放射される音波または電波を受動センサにより測定した観測データを用い、移動体の航跡（移動体の速力、針路、方位、距離等）を分析する移動体の航跡分析装置及び分析方法に関する。
【０００２】
【従来の技術】
一般に、移動体の航跡分析、例えば、船舶航跡分析は、移動している観測船が、一定時間間隔毎の各時刻に音響センサを使って被観測船舶（以下、単に船舶という）が放射する音波を受信し、到来音波の方位を観測し、ここで観測された方位、観測時刻及び観測船の位置の時系列のデータに基づき船舶の航跡を分析するという方法で行われている。
【０００３】
前述のような観測方位、観測時刻、観測船の位置の時系列データによる船舶航跡分析方法であってカルマンフィルタを用い、方位情報のみを入力データとする航跡分析方法に関する従来技術として、例えば、非特許文献１等に記載された技術が知られている。
【０００４】
【非特許文献１】
ＩＥＥＥ会報１９８５年１０月Ｖｏｌ．ＡＣ−３０，  Ｎｏ．１０，  ９４０〜９４９ページ「Ａ  Ｓｔｏｃｈａｓｔｉｃ  Ａｎａｌｙｓｉｓ  ｏｆ  ａ  Ｍｏｄｉｆｉｅｄ  Ｇａｉｎ  Ｅｘｔｅｎｄｅｄ  Ｋａｌｍａｎ  Ｆｉｌｔｅｒ  ｗｉｔｈ  Ａｐｐｌｉｃａｔｉｏｎ  ｔｏ  Ｅｓｔｉｍａｔｉｏｎ  ｗｉｔｈ  Ｂｅａｒｉｎｇｓ  Ｏｎｌｙ  Ｍｅａｓｕｒｅｍｅｎｔｓ」
【０００５】
【発明が解決しようとする課題】
前述した従来技術による移動体としての船舶の航跡分析方法により船舶の航跡分析を行っている場合、観測対象となる船舶の船舶航跡のいずれかの値がオペレータの先見情報等により推定されることがある。このような先見情報は、オペレータの観測経験に基づくものであり、各種の受動センサからの観測データに加えて有効に活用されることが望ましい。
【０００６】
しかし、前述した従来技術による船舶航跡分析方法は、このようなオペレータの先見情報を用いることについて全く考慮されておらず、観測結果の分析収束時間が多くなるという問題点を有している。
【０００７】
本発明の目的は、前述した従来技術による問題点を解決し、移動体のの航跡分析処理に関するオペレータの先見情報の利用を可能とし、これにより、より高精度な移動体の航跡分析を実現することができる移動体の航跡分析装置及び分析方法を提供することにある。
【０００８】
【課題を解決するための手段】
本発明によれば前記目的は、移動体から放射される音波または電波を検出するセンサにより時系列に測定した観測データを用い、前記移動体の航跡を分析する移動体の航跡分析装置において、前記時系列に測定した観測データに基づいて前記移動体の航跡の分析処理を実行する分析処理部と、オペレータの操作を受け付ける端末部とを備え、前記端末部が、オペレータからの航跡に関する先見情報を取り入れ、前記分析処理部が、前記端末部で受け付けたオペレータからの先見情報を反映させて前記移動体の航跡の分析を行うことにより達成される。
【０００９】
本発明による移動体の航跡分析装置及び分析方法は、移動体のの航跡に関するオペレータの先見情報を分析処理で利用可能とする手段、分析途中に得られる先見情報を各航跡の推定値または初期値に反映させる手段、先見情報の妥当性を評価する手段等を提供することにより、より高精度な移動体の航跡分析を実現することができる。
【００１０】
オペレータの先見情報の具体的な利用方法としては、例えば、分析途中において得られたオペレータの先見情報や各航跡分析値を、オペレータの判断によって初期値に反映させて再度初観測データから分析を始める、または、オペレータが入力した先見情報を反映した分析結果と観測データとの誤差を、時系列残差として表示する等がある。
【００１１】
また、本発明は、オペレータが入力した先見情報を反映した分析結果が真の航跡に収束しているか否かを評価するための時系列残差表示手段を有する。この時系列残差表示手段は、例えば、移動体の航跡分析処理の結果得られたデータと観測データとの時系列残差を示す表示か、あるいは、前記時系列残差に前記観測データの精度情報を付加した表示を行う。さらに、本発明は、これらの表示に加えて、時系列残差の２乗平均または時系列正規化残差２乗平均の数値情報を表示するようにすることもできる。
【００１２】
さらに、本発明は、前述の表示に加えて、表示している分析結果や実行されている分析処理に、オペレータによる先見情報が用いられている場合に、これをオペレータや他のユーザに音声や画像にて告知するようにすることができる。
【００１３】
本発明によれば、オペレータの移動体の航跡に関する先見情報を用いて分析を行うことにより、移動体が放射する音波または電波に関する観測データのみによって航跡分析を実行する従来方法と比べ、航跡分析値が移動体の真の航跡に収束するまでの時間を改善することが可能となる。
【００１４】
また、本発明によれば、分析途中に得られた先見情報ならびに各航跡分析値を初期値に反映させることにより、収束後の航跡分析値と移動体の真の航跡の誤差を改善することが可能となる。
【００１５】
さらに、本発明によれば、オペレータが入力した先見情報を反映した分析結果を時系列残差として表示することにより、先見情報の妥当性を評価することが可能となる。
【００１６】
【発明の実施の形態】
以下、本発明による本発明による移動体の航跡分析装置及び分析方法の実施形態を図面により詳細に説明する。なお、以下に説明する本発明の実施形態は、移動体を船舶として本発明を適用したものとして説明するが、本発明が適用可能な移動体としては、船舶のみならず航空機等の音波または電波を放射しながら移動する各種の移動体であってよい。
【００１７】
図１は本発明の一実施形態による船舶航跡分析装置の構成例を示すブロック図である。図１において、１は受動センサ、２は到来音波または電波観測装置、３は到来音波または電波観測データベース、４は分析処理器、５は船舶航跡分析値表示器、６は初期条件演算器、７は観測行列演算器、８は端末部、９は時系列残差演算器である。
【００１８】
一般に、受動センサには、音波または電波の到来方位のみを測定する装置、音波または電波の到来方位と周波数とを測定する装置があり、また、船舶航跡分析方法としては、重回帰分析法、最尤法、遺伝的アルゴリズム等を用いた分析法が知られている。ここで説明する本発明の実施形態は、受動センサから音波到来方位の観測データが得られる条件の下で、カルマンフィルタを応用した船舶航跡分析法の一例である。
【００１９】
本発明の実施形態による船舶航跡分析装置は、図１に示すように、到来音波または電波観測データベース３、分析処理器４、船舶航跡分析値表示器５、初期条件演算器６、観測行列演算器７、端末部８、及び、時系列残差演算器９を備えて構成されるものであって、例えば、船舶航跡分析方法を実施するためのプログラムを実行する中央処理部、メモリ、表示部等を具備するコンピュータシステムにより実現される。
【００２０】
本発明の実施形態による船舶航跡分析装置において、受動センサ１による受信信号１０は、到来音波または電波観測装置２により観測され、その観測データ１１は、観測した船舶の方位の情報として到来音波または電波観測データベース３に格納される。オペレータの先見情報及び各航跡分析値が収束したとの判断がない場合、到来音波または電波観測データベース３に格納されている観測方位の一番新しいものが、到来音波または電波観測データベース３→ライン１２→分析処理器４の経路で分析処理器４の観測データとして入力され、分析処理器４で更新処理が行われ、その分析処理結果が分析処理器４→ライン１３→船舶航跡分析値表示器５の経路で船舶航跡分析値表示器５に出力される。
【００２１】
前述した処理は、一定時間毎の観測により、到来音波または電波観測データベース３に新しいデータが貯えられる毎に行われ、逐次推定更新の処理が行われる。また、時系列残差演算器９は、到来音波または電波観測データベース３→ライン２２→時系列残差演算器９の経路で到来音波または電波観測データベース３から送られてくる観測データと、分析処理器４→ライン２３→時系列残差演算器９の経路で分析処理器４から送られてくる分析値とから、時系列残差を計算し、その結果である時系列残差を、時系列残差演算器９→ライン２４→船舶航跡分析値表示器５の経路で船舶航跡分析値表示器５に表示する。前述の時系列残差の算出は、分析処理器４から送られてくる現在の分析値に基づいて、１つ前の観測データを採取した時刻における予測観測データを演算し、この予測観測データとデータベース３からの１つ前の観測データとの差を演算することにより行われる。
【００２２】
端末部８は、入力された操作内容を分析し、該分析結果に応じた出力信号を出力する機能を備えているものとする。例えば、端末部８に、分析対象となる船舶の航跡のいずれかの推測値等がオペレータの先見情報として入力された場合、端末部８は、端末部８→ライン２１→分析処理器４の経路で、分析処理器４にオペレータの介入があることを知らせ、同様に、端末部８→ライン１８→観測行列演算器７→ライン１９→分析処理器４の経路で、端末部８に入力されたオペレータの先見情報に応じた観測行列を観測行列演算器７により計算させ、その計算結果を分析処理器４に通知する。また、同様に、端末部８→ライン２０→分析処理器４の経路でオペレータの各先見情報が分析処理器４に通知される。分析処理器４は、これらの値を用いた分析を行い、その分析結果を分析処理器４→ライン１３→船舶航跡分析値表示器５の経路で船舶航跡分析値表示器５に表示する。
【００２３】
また、端末部８に、初期条件の変更に関するオペレータの判断、例えば、各航跡が収束したとの判断が入力されると、端末部８は、端末部８→ライン１７→分析処理器４の経路で、オペレータによる各航跡の収束判断を示す信号を分析処理器４内の判断更新処理部４１に通知する。同様に、端末部８は、端末部８→ライン１４→初期条件演算器６の経路で、前述のオペレータによる各航跡の収束判断を初期値演算器６に通知する。一方、分析処理器４は、分析処理器４→ライン１６→初期値演算器６の経路で初期値演算器６に状態ベクトルを通知する。
【００２４】
初期値演算器６は、通知された前述の判断及び状態ベクトルにより初期値を計算し、その初期値を初期値演算器６→ライン１５→分析処理器４の経路で分析処理器４に通知する。また、前述と同様に、端末部８は、端末部８→ライン２５→到来音波または電波観測データベース３の経路で前述の収束判断を到来音波または電波観測データベース３に通知する。これにより、到来音波または電波観測データベース３は、到来音波または電波観測データベース３→ライン１２→分析処理４の経路で、格納している観測データをその最初のものから逐次分析処理器４に送信する。前述したようにして分析された船舶航跡分析諸元分析値が分析処理器４→ライン１３→船舶航跡分析値表示器５の経路で船舶航跡分析値表示器５に示される。
【００２５】
また、端末部８は、船舶航跡分析値表示器５の表示内容の選択に関する選択指示や、先見情報として入力した航跡に関する推測値を分析処理にどのように反映させるか等の分析処理に関する指示等の前述したオペレータの先見情報や判断の他に、オペレータからの操作指示をも受け付けるものとする。
【００２６】
また、端末部８は、先見情報の入力を受け付けることに加えて、先見情報入力後も継続して先見情報を入力するか、先見情報に基づき推定初期値を変更して再度分析を行うかを選択するための選択指示をも受け付けるものとする。
【００２７】
本発明の実施形態で分析方法の一例として挙げるカルマンフィルタについては後述において概説するが、このカルマンフィルタは、観測データと分析値との誤差を最小にするシステムの状態を、観測データが与えられる度に逐次推定更新するものである。従って、これを船舶航跡分析に応用することにより、誤差を含んだ入力データから船舶の航跡の分析が可能となる。すなわち、カルマンフィルタを要約すれば、
（ａ）  信号を生成するシステムの動特性
（ｂ）  雑音の統計的性質
（ｃ）  時々刻々与えられる観測データ
を用いてシステムの状態の最小２乗推定値を逐次的に与えるオンラインデータアルゴリズムである。
【００２８】
図３はカルマンフィルタの一般的な応用例について説明する図であり、次に、これについて説明する。図３において、システム３１は、既知の制御入力及び外乱（雑音）を受けており、システム３１の状態は観測装置３２を通して観測される。通常観測データは雑音に乱されていて、また、いくつかの状態は直接観測されないことも多い。このような状況の下で、観測データからシステム３１の状態に関する情報摘出を行うのがカルマンフィルタ３３である。
【００２９】
以上のことから、本発明の実施形態でオペレータの先見情報を入力としたカルマンフィルタを応用した分析を行う際、この先見情報の精度が高くなくても適切に利用すれば、船舶航跡分析の性能向上に役立つといえる。また、上に述べた、システムの状態を、観測データが与えられる度に逐次更新するということは、精度の高い初期値を与えれば、精度の高い分析値を得ることができることを表し、分析途中の各航跡分析値からカルマンフィルタ応用アルゴリズムの初期条件を算出し、この初期条件のもとで、再度観測データの初めから船舶航跡分析を始める初期条件変更リトライ処理を行うことにより、より精度のよい分析値を得ることが可能となる。
【００３０】
本発明の実施形態におけるカルマンフィルタの応用例について説明する前に、以下で、既存のカルマンフィルタを応用した船舶航跡分析のアルゴリズムと、海上における航行船舶の船舶航跡分析とを例にして説明をする。
【００３１】
図２は船舶と観測船との関係から水上を運動する船舶の航跡を分析する原理を説明する図であり、次に、これについて説明する。
【００３２】
船舶及び観測船が図２に示すように運動しているものとすると、観測船は、時刻ｔ、ｔ＋１、ｔ＋２、……に毎に、刻々と真の船舶方位Ｂｙを観測することができる。この観測船舶方位Ｂｙに正規性白色雑音が加わったものが観測方位として離散的に与えられる。このような離散的時間確率システムに対するカルマンフィルタは次式で与えられる。
【００３３】
●フィルタ方程式
【数１】
2004239678.tif 000003
【数２】
2004239678.tif 000004
●カルマンゲイン
【数３】
2004239678.tif 000005
●推定誤差共分散行列
【数４】
2004239678.tif 000006
【数５】
2004239678.tif 000007
●初期条件
【数６】
2004239678.tif 000008
【００３４】
数１〜数６において、
ｘ（ｔ＋１／ｔ）：  状態ベクトル予測値
ｘ（ｔ／ｔ）：  状態ベクトル濾波値
Δｘ（ｔ）：  状態ベクトル更新量
ｙ（ｔ）：  観測ベクトル
Ｆ（ｔ）：  状態遷移行列
Ｋ（ｔ）：  カルマンゲイン
Ｈ（ｔ）：  観測行列
Ｐ（ｔ）：  推定誤差共分散行列
ΔＰ（ｔ）：  推定誤差共分散行列更新量
Ｒ（ｔ）：  正定値対称行列
ｘ（０）：  予測初期値
Σ（０）：  共分散の初期値
なお、添字の（ｔ＋１／ｔ）、（ｔ／ｔ）等は次のような意味を持つ。すなわち、
（ｔ＋１／ｔ）：  時刻ｔまでのデータに基づく時刻ｔ＋１における予測値、
（ｔ／ｔ）：  時刻ｔまでのデータに基づく時刻ｔにおける濾波値である。
【００３５】
さらに、（ｔ）はΔｘ（ｔ），Ｆ（ｔ），Ｋ（ｔ），Ｈ（ｔ），Ｐ（ｔ），ΔＰ（ｔ），Ｒ（ｔ）に関しては時刻ｔまでのデータからの計算値、ｙ（ｔ）は時刻ｔでの観測データ、ｘ（ｔ）は状態ベクトルの真値であるとする。また、前述の各ｔは初観測時刻を０とした、サンプリング間隔ΔＴによって離散化された時刻を表すこととする。また、本来、カルマンフィルタアルゴリズムでは前述の数４に更に、Ｇ（ｔ）Ｑ（ｔ）Ｇ（ｔ）Ｔなる項が追加されるが、ここで説明する例では、システムノイズを考慮しなかったのでこの項は０となる。
【００３６】
前述において、まず、ｙ（０），ｘ（０），Σ（０）が与えられると、前記数６により初期値ｘ（０／−１），Ｐ（０／−１）が定まり、前記数３によりＫ（０）が計算される。更に、前記数２、数５によりｘ（０／０），Ｐ（０／０）が求まる。ついで、前記数１、数４によりｘ（１／０），Ｐ（１／０）、前記数３によりＫ（１）が計算される。ｔ＝１においてｙ（１）を観測し再び上記数２、数５よりｘ（１／１），Ｐ（１／１）が求められるので、以下同様の手順で推定値及び推定誤差共分散行列が計算される。このようにして新しい観測データが得られる度に古い推定値を修正して新しい推定値が計算される。
【００３７】
前述で説明した処理過程がいわゆる一般のカルマンフィルタである。そして、これはシステムが線形であるという仮定に基づいている。このため、船舶航跡分析は、本質的に非線形であるため構造式中に表れる非線形な関数を適当に変形して疑似線形化してカルマンフィルタを適用している。しかし、非線形な関数であってもこれが滑らかであるという仮定の下に推定値のまわりでテーラー展開することにより線形化することができる。このようにして得られるものが拡張カルマンフィルタのアルゴリズムであり、次のような式で表される。
【００３８】
●フィルタ方程式
【数７】
2004239678.tif 000009
【数８】
2004239678.tif 000010
●カルマンゲイン
【数９】
2004239678.tif 000011
●推定誤差共分散行列
【数１０】
2004239678.tif 000012
【数１１】
2004239678.tif 000013
●初期条件
【数１２】
2004239678.tif 000014
【００３９】
数７〜数１２において、ｘ（ｔ＋１／ｔ），ｘ（ｔ／ｔ），Δｘ（ｔ），ｙ（ｔ），Ｋ（ｔ），Ｐ（ｔ），ΔＰ（ｔ），Ｒ（ｔ），ｘ（０），Σ（０）は前述した一般のカルマンフィルタのものと同一であるが、ここで前述数７、数８で新しく現れたｆ，ｈは前記数１、数２のＦ（ｔ），Ｈ（ｔ）と同様の機能を持ち、それぞれ、状態ベクトルを状態遷移させ、状態ベクトルから観測データを推定する関数であり、前記数９、数１０、数１１のＦ（ｔ），Ｈ（ｔ）は、
【数１３】
2004239678.tif 000015
【数１４】
2004239678.tif 000016
で定義される行列である。
【００４０】
また、本来拡張カルマンフィルタアルゴリズムは、前記数１０にさらにＱ（ｔ）なる項が追加されるが、ここで説明している例では、システムノイズを考慮しなかったのでこの項は０となる。上記数７〜数１２は、前記数１〜数６にきわめて類似している。ｆ，ｈが線形であれば、このアルゴリズムは、カルマンフィルタに帰着する。カルマンフィルタにおいて、推定誤差共分散行列Ｐ（ｔ／ｔ−１），Ｐ（ｔ／ｔ）及びカルマンゲインＫ（ｔ）は、観測データとは無関係でありこれらはオフラインで計算することができる。しかし、前記数１３、数１４からＦ（ｔ），Ｈ（ｔ）は、ｘ（ｔ／ｔ），ｘ（ｔ／ｔ−１）の関数であり実際の観測データに依存する。従って、前記数９〜数１１も観測データに依存するのでＫ（ｔ）をオフラインで計算しておくことはできない。
【００４１】
図４は海上における航行船舶の船舶航跡分析のモデルを説明する図であり、前述した拡張カルマンフィルタを用いた船舶航跡分析のためのものである。まず、図４及び以後に説明する変数の定義を以下に示す。
【００４２】
座標原点：  初探知時刻における観測船の位置
ｒ：  座標原点からの船舶距離
θ：  座標原点からの船舶方位
ｘ０：  座標原点からの観測船ｘ軸方向距離
ｙ０：  座標原点からの観測船ｙ軸方向距離
Ｘｔ：  観測船から船舶のＸ軸方向距離
Ｙｔ：  観測船から船舶のＹ軸方向距離
Ｒ：  観測船からの船舶距離
Ｂｙ：  観測船からの観測方位
Ｍｔ：  速度
Ｃｔ：  針路
なお、ｒ，θ，Ｘｔ，Ｙｔ，Ｒ，Ｂｙ，Ｍｔ，Ｃｔに関し、添字（ｔ＋１／ｔ）がついた場合、時刻ｔまでのデータに基づく時刻ｔ＋１における予測値とし、これは以下に示す方法でｘ（ｔ＋１／ｔ）より求め、添字（ｔ／ｔ）がついた場合、時刻かまでのデータに基づく時刻かにおける濾波値とし、これは同様に以下に示す方法でｘ（ｔ＋１／ｔ）より求め、添字（ｔ）がついた場合には真値とする。ただし、ｘ０，ｙ０に関しては、ｘ０（ｔ），ｙ０（ｔ）しかなく、これは時刻ｔにおける観測データであるとする。また、図４において、上に・が付されたθは、θの微分値を示し、上に・が付されたｒは、ｒの微分値を示す。
【００４３】
そして、状態ベクトルを、
【数５７】
2004239678.tif 000017
と定義すると、
状態ベクトル更新量Δｘ（ｔ）は、
【数５８】
2004239678.tif 000018
となる。
【００４４】
次に、状態遷移行列の求め方について説明する。
【００４５】
図５は状態遷移行列を説明するための図であり、船舶が等速直進運動をしていると仮定して船舶の運動を示している。図５において、座標原点は観測時刻における観測船の位置、ｒは座標原点からの船舶距離、θは座標原点からの船舶方位である。
【００４６】
図５において、
【数１５】
2004239678.tif 000019
【数１６】
2004239678.tif 000020
が成り立つ。
【００４７】
この数１５、数１６により
【数１７】
2004239678.tif 000021
【数１８】
2004239678.tif 000022
が得られる。
【００４８】
数１７、数１８をそれぞれ２乗して加えると、
【数１９】
2004239678.tif 000023
を得ることができる。
【００４９】
また、図５に示すモデルにおいて、
【数２０】
2004239678.tif 000024
【数２１】
2004239678.tif 000025
が成り立ち、前記数１９、数２０より、
【数２２】
2004239678.tif 000026
を得ることができる。
【００５０】
そして、前記数１７を数１８で割ると、
【数２３】
2004239678.tif 000027
を得ることができる。
【００５１】
上記数２０は、
【数２４】
2004239678.tif 000028
のように書くこともでき、前記数２１、数２３、数２４から、
【数２５】
2004239678.tif 000029
を得ることができる。
【００５２】
以上により上記数７のｆが求められ、状態遷移行列は、
【数２６】
2004239678.tif 000030
と表すことができる。この偏微分を解いて状態遷移行列を得る。但し、前述の説明は、状態ベクトルの真値を例にしているが、実際の分析では濾波値を使用して状態遷移を推定する。
【００５３】
正定値行列は、このモデルにおいて観測雑音行列を意味し、分散値として音響センサより与えられるものとする。
【００５４】
ここで、Ｘｔ，Ｙｔ，Ｒ，Ｂｙ，Ｍｔ，Ｃｔの値は、状態ベクトルの要素を用いて以下に示す数式により求めることができる。
【００５５】
【数２７】
2004239678.tif 000031
【数２８】
2004239678.tif 000032
【数２９】
2004239678.tif 000033
【数３０】
2004239678.tif 000034
【数３１】
2004239678.tif 000035
【数３２】
2004239678.tif 000036
【数３３】
2004239678.tif 000037
【数３４】
2004239678.tif 000038
【００５６】
ここでは、前記の変換を便宜上以下のように、
【数３５】
2004239678.tif 000039
と表すものとする。
【００５７】
観測行列Ｈ（ｔ）は、船舶航跡の推定値を状態ベクトルｘ（ｔ／ｔ−１）で微分して、以下のように求めることができる。
【００５８】
【数３６】
2004239678.tif 000040
数３６において、α（ｔ）は船舶航跡の推定値であるとする。
【００５９】
ここで例示している船舶航跡分析で使用されている拡張カルマンフィルタは、音響センサで捕らえた船舶の方位情報をもとに分析を行うので、α（ｔ）は、観測方位情報となる。この場合の観測方位は、座標原点からの船舶の推定方位、推定距離、観測船位置から
【数３７】
2004239678.tif 000041
と表すことができるので、
【数３８】
2004239678.tif 000042
となる。
【００６０】
それぞれの成分は、微分により、
【数３９】
2004239678.tif 000043
として求められる。
【００６１】
図６は前述した拡張カルマンフィルタを用いた船舶航跡分析処理の一例を説明するフローチャートであり、次に、これについて説明する。
【００６２】
（１）処理が開始されて船舶が探知されると、まず、ｙ（０），ｘ（０），Σ（０）が与えられ、数１２に示す初期条件（ｘ（０／−１），Ｐ（０／−１））の設定を行い、分析時間の初期化（ｔ＝０）を行う（ステップ６０１、６０２）。
【００６３】
（２）次に、数３６〜数３９による観測行列Ｈ（１）の設定を行うと共に、正定値対称行列Ｒ（１）を設定する（ステップ６０３）。
【００６４】
（３）次に、前記数９により定数Ｋ（０）が計算される。さらに、前記数８、数１１によりｘ（０／０），Ｐ（０／０）を求める（ステップ６０４、６０５）。
【００６５】
（４）その後、数２６により、状態遷移行列Ｆ（１）を設定し、また、前記数７、数１０によりｘ（１／０），Ｐ（１／０）が計算され、状態ベクトル時間更新、推定誤差共分散行列更新の処理を行う。前述により得られた状態ベクトル時間更新のデータｘ（１／０）が現時点での船舶位置を示すことになる（ステップ６０６、６０７）。
【００６６】
分析時間ｔをｔ＋１に更新、この場合、ｔ＝１に更新し、ステップ６０３からの処理に戻り、観測行列Ｈ（１）の設定を行うと共に、前記数９によりカルマンゲインＫ（１）を計算し、ｔ＝１において、ｙ（１）を観測し再び前記数７、数１０によりｘ（１／１），Ｐ（１／１）が求められ、以下、同様の手順で推定値及び推定誤差共分散行列が計算される（ステップ６０８）。
【００６７】
前述したようにして、新しい観測データが得られる度に古い推定値が修正されて新しい推定値が計算される。
【００６８】
図７は拡張カルマンフィルタにオペレータ介入分析処理を加えた船舶航跡分析処理の一例を説明するフローチャートであり、次に、これについて説明する。ここでの処理は、オペレータが有する船舶の各航跡に対する先見情報を図６により説明した航跡分析処理にオペレータ介入分析処理のアルゴリズムを取り入れ、拡張カルマンフィルタを用いた場合の分析処理の例である。
【００６９】
図７におけるステップ７０１〜７０５、７１２〜７１４の処理は、それぞれ、図６のステップ６０１〜６０５、６０６〜６０８の処理と同一の処理であり、以下では、図６に説明した処理と異なる処理について詳細に説明する。以下では、オペレータの船舶速度に対する先見情報値をＭｔｏＰ  、船舶距離に対する先見情報値をＲｏＰ、船舶針路に対する先見情報値をＣｔｏＰと表す。
【００７０】
また、ここで説明する例は、オペレータの船舶の航跡に対する前述の先見情報値ＭｔｏＰ，ＲｏＰ，ＣｔｏＰが個別に入力されている場合、従来の観測方位情報をもとにした分析後に、Ｍｔ，Ｒ，Ｃｔをもとにした分析を追加して行う例である。なお、ここで説明する例において、ＭｔｏＰ，ＲｏＰ，ＣｔｏＰが同時に入力されている場合には、それぞれの入力に対する分析をシーケンシャルに処理する。
【００７１】
さらに、ここで説明するモデルは、正定値行列が入力の観測雑音行列を表しているが、オペレータの先見情報の入力に対しては、先見情報の誤差がどれぐらいなのか、オペレータが各先見情報に対する誤差の分散として入力することができる方法を備えているものとする。ここで、前述のこれらＭｔ，Ｒ，Ｃｔについての観測行列を導出する。
【００７２】
●Ｍｔに関する観測行列
Ｍｔを入力として分析を行う場合、
【数４０】
2004239678.tif 000044
となり、
【数４１】
2004239678.tif 000045
とすると、微分により各成分は、
【数４２】
2004239678.tif 000046
となる。
【００７３】
●Ｒに関する観測行列
Ｒに関しては、状態ベクトルとの関係より、その逆数をとり、
【数４３】
2004239678.tif 000047
となり、よって
【数４４】
2004239678.tif 000048
とすると、それぞれの成分は、微分により、
【数４５】
2004239678.tif 000049
となる。
【００７４】
●Ｃｔに関する観測行列
Ｃｔを入力として分析を行う場合、
【数４６】
2004239678.tif 000050
となり、
【数４７】
2004239678.tif 000051
とすると、微分により各成分は、
【数４８】
2004239678.tif 000052
となる。
【００７５】
ここで、前記数３８に示す船舶の方位情報から求めた観測行列を、他の観測行列の表記にあわせてＨＢｙ（ｔ）  と表すことにする。
【００７６】
図７のフローチャートに示すアルゴリズムにおいて、図６のフローの場合と同様に、まず、ｙ（０），ｘ（０），Σ（０）が与えられると、前記数１２により初期値ｘ（０／−１），Ｐ（０／−１）が定まり、前記数９により定数Ｋ（０）が計算される。さらに、前記数８、数１１によりｘ（０／０），Ｐ（０／０）を求める（ステップ７０１〜７０５）。
【００７７】
この結果を通常修正更新とし、次に、オペレータの先見情報値（ＭｔｏＰ，ＲｏＰ，ＣｔｏＰ）が入力されているか否かを判定し、入力されていた場合、この先見情報値を観測ベクトル（ｙ（ｔ））とし、観測行列（Ｈ（ｔ））をそれぞれの先見情報に応じた観測行列とし、ｘ（０／０），Ｐ（０／０）に対して修正更新をして分析を行う。この場合のｘ（ｔ／ｔ），Ｐ（ｔ／ｔ）に対する修正更新の式は、
【数４９】
2004239678.tif 000053
【数５０】
2004239678.tif 000054
と示すことができる。
【００７８】
例えば、ステップ７０５の処理の後、速度に対する先見情報値ＭｔｏＰが入力されているか否かを判定し、入力されていた場合、通常修正更新の後、Ｈ（ｔ）を前記数４１に示すＨ（ｔ）Ｍｔとし、ｙ（ｔ）をＭｔｏＰとして、前記数４９、数５０により再び修正更新を行う（ステップ７０６、７０７）。
【００７９】
次に、距離に対する先見情報値ＲｏＰが入力されているか否かを判定し、入力されていた場合、通常修正更新の後、Ｈ（ｔ）を前記数４４に示すＨ（ｔ）Ｒとし、ｙ（ｔ）をＲｏＰとして、前記数４９、数５０により再び修正更新を行う（ステップ７０８、７０９）。
【００８０】
次に、進路に対する先見情報値ＣｔｏＰが入力されているか否かを判定し、入力されていた場合、通常修正更新の後、Ｈ（ｔ）を前記数４７に示すＨ（ｔ）Ｃｔとし、ｙ（ｔ）をＣｔｏＰとして、前記数４９、数５０により再び修正更新を行う（ステップ７１０、７１１）。
【００８１】
その後、通常通り前記数７、数１０によりｘ（１／０），Ｐ（１／０）が計算される（ステップ７１２、７１３）。
【００８２】
分析時間ｔをｔ＋１に更新、この場合、ｔ＝１に更新し、ステップ７０３からの処理に戻り、前記数９によりカルマンゲインＫ（１）を計算し、ｔ＝１において、ｙ（１）を観測し再び前記数７、数１０によりｘ（１／１），Ｐ（１／１）が求められ、以下、同様の手順で推定値及び推定誤差共分散行列が計算される（ステップ７１４）。
【００８３】
前述したようにして、新しい観測データが得られる度に、オペレータの船舶の航跡に対する先見情報値を含めて、古い推定値が修正されて新しい推定値が計算される。
【００８４】
図８は拡張カルマンフィルタにオペレータ介入分析処理及び初期条件変更リトライ処理を加えた船舶航跡分析処理の一例を説明するフローチャートであり、次に、これについて説明する。ここでの処理は、オペレータの判断により、船舶航跡分析の途中の観測時刻ｔ＝ｍにおいて、各航跡分析値を航跡推定初値に反映させる方法（初期条件変更リトライ処理）である。以下では、図７により説明したフローに示す船舶航跡分析処理に、初期条件変更リトライ処理を付加したものを例として説明する。図８におけるステップ８０１〜８１１、８１４〜８１６の処理は、それぞれ、図７に示すフローのステップ７０１〜７１１、７１２〜７１４の処理と同一の処理であり、以下では、図７のフローの処理と異なる処理について説明する。
【００８５】
ｔ＝ｍのときの各航跡分析値は、ステップ８０１〜８１１までの処理により、船舶距離については  Ｒ（ｍ／ｍ）、観測方位についてはＢｙ（ｍ／ｍ）、速度についてはＭｔ（ｍ／ｍ）、的針についてはＣｔ（ｍ／ｍ）となっている。ここで、ｔ＝ｍにおける各航跡をもとに船舶初観測時（ｔ＝０）における各航跡を考える。説明している本発明実施形態は、船舶が基本的に等速直進運動をしていると仮定しているので、ｔ＝ｍにおける各航跡分析値から求めるｔ＝０における船舶の的速は、  Ｍｔ（０）（ｍ／ｍ）＝Ｍｔ（ｍ／ｍ）、進路は、Ｃｔ（０）（ｍ／ｍ）＝Ｃｔ（ｍ／ｍ）となる。なお、これらの式における左辺の添え字の表記はｘ（ｍ／ｍ）を使用して求めた観測時刻ｔ＝０の値を意味する。
【００８６】
また、距離については、ｔ＝ｍにおける各航跡から、
【数５１】
2004239678.tif 000055
となる。
【００８７】
以上のＭｔ（０）（ｍ／ｍ），Ｃｔ（０）（ｍ／ｍ），Ｒ（０）（ｍ／ｍ）及びｙ（０）を用いて、前記数１２におけるｘ（０）は、
【数５２】
2004239678.tif 000056
で与えられることとなる。
【００８８】
図８に示している分析処理では、ステップ８０１〜８１１の処理の終了後、端末部８を介してオペレータからの初期条件の変更に関する判断、例えば、航跡のいずれが収束したとの判断である初期条件変更リトライ処理要求が入力されたか否かを判定し、前述の要求が入力されていた場合、前記数５２で得られたｘ（０）を用い、到来音波方位観測データベースに格納してある到来音波方位のデータからの解析モードとして、ステップ８０３からの処理に戻って、再度観測データの初めから分析を行う（ステップ８１２、８１３）。
【００８９】
前述で説明した船舶航跡分析値を船舶航跡推定初期値に反映させる方法の他に、オペレータが分析途中に得た船舶の航跡の先見情報を船舶航跡推定初期値に反映させる方法もあるが、この方法は、この先見情報を状態ベクトルに直しその値をｘ（０）を用い、到来音波方位観測データベースに格納してある到来音波方位のデータを用いて、再度観測データの初めから分析を行えばよい。
【００９０】
次に、オペレータが入力する先見情報の妥当性を評価する方法について説明する。
【００９１】
いま、観測時刻ｔ＝ｋのときの各航跡分析値は、前述の処理により、船舶距離についてはＲ（ｋ／ｋ）、観測方位についてはＢｙ（ｋ／ｋ）、速度についてはＭｔ（ｋ／ｋ）、進路についてはＣｔ（ｋ／ｋ）として求められており、また、本発明の実施形態では船舶は、基本的に等速直進運動をしていると仮定しているので、上記Ｒ（ｋ／ｋ）、Ｂｙ（ｋ／ｋ）、Ｍｔ（ｋ／ｋ）、Ｃｔ（ｋ／ｋ）より求められるｔ＝ｎ（０＜ｎ＜ｋ）における観測方位推定値は、
【数５３】
2004239678.tif 000057
と表すことができる。
【００９２】
そして、時系列残差は、この推定方位と観測データとの差より
【数５４】
2004239678.tif 000058
と表すことができる。
【００９３】
前記数５４を用いれば、現在の分析値と過去の各時刻におけるデータとのずれを時系列的に表すことができ、これをグラフィカルに表示すれば解の妥当性を容易に評価することが可能となる。
【００９４】
図９は前述の時系列残差をグラフィック表示により表した例を説明する図、図１０は観測精度が付加された時系列残差をグラフィック表示により表した例を説明する図である。
【００９５】
図９に示す時系列残差のグラフィック表示は、図の上から下に時刻ｔ＝ｎ＝０（初観測時刻）から時刻ｔ＝ｎ＝ｋ（現在）までの数５４による残差値を、中心線を０として横軸方向に示したものである。また、個々で説明している例では、音響センサより観測雑音として、方位の誤差が分散値として与えられるので、図１０に示すようにデータのずれと共にこの分散値を示すことにより、オペレータに対して誤差の大きい観測データによる残差であることを示すことができる。
【００９６】
さらに、図９、図１０のグラフには、時系列残差あるいは時系列正規化残差の２乗平均を示すようにしており、これにより、解の妥当性をより明確に示すことができる。ここで、時系列残差２乗平均は、
【数５５】
2004239678.tif 000059
として表され、この表示により、解の妥当性が、現在の分析値と過去の観測データとのずれの２乗総和の平均値として確認できる。この表示は、分析値が真値に向かって収束しているのか否かの指標となり、オペレータの先見情報を入力しているときには、その先見情報の妥当性を表すものとなる。この数値の表示例は、図９に示されている。
【００９７】
また、先見情報が入力されている場合には、その状況として、先見情報使用中の表示９１をグラフと共に表示したり、あるいは、別途音や光で前記状況を示すための告知手段を備える構成としてもよい。
【００９８】
また、時系列正規化残差２乗平均は
【数５６】
2004239678.tif 000060
として表され、前記数５５に観測誤差の概念を取り入れたものである。この数値の表示例は、図１０に示している。
【００９９】
前述した本発明の実施形態での各処理は、処理プログラムとして構成することができ、この処理プログラムは、ＨＤ、ＤＡＴ、ＦＤ、ＭＯ、ＤＶＤ−ＲＯＭ、ＣＤ−ＲＯＭ等の記録媒体に格納して提供することができる。
【０１００】
前述した本発明の実施形態によれば、船舶航跡のいずれかの値が推定可能である場合、この推定値に関するオペレータからの入力を反映して船舶航跡分析値が船舶の真の航跡に収束するまでの時間を改善することができる。
【０１０１】
また、本発明によれば、分析途中で入力された推定値ならびに各航跡分析値を初期値に反映して再度分析するようにすることにより、収束後の航跡分析値の船舶の真の航跡に対する誤差を改善することができる。
【０１０２】
さらに、本発明によれば、オペレータが入力した推定値を反映した分析結果を時系列残差として表示することにより、オペレータの操作による分析収束時間を低減させることが可能となる。
【０１０３】
【発明の効果】
以上説明したように、本発明によれば、移動体の航跡分析処理に関するオペレータの先見情報の利用を可能とし、これにより、より高精度な移動体航跡分析を短時間に実現することができる。
【図面の簡単な説明】
【図１】本発明の一実施形態による船舶航跡分析装置の構成例を示すブロック図である。
【図２】船舶と観測船との関係から水上を運動する船舶の航跡を分析する原理を説明する図である。
【図３】カルマンフィルタの一般的な応用例について説明する図である。
【図４】海上における航行船舶の船舶航跡分析のモデルを説明する図である。
【図５】船舶が等速直進運動をしていると仮定して船舶の運動を示すと共に、状態遷移行列を説明するための図である。
【図６】拡張カルマンフィルタを用いた船舶航跡分析処理の一例を説明するフローチャートである。
【図７】拡張カルマンフィルタにオペレータ介入分析処理を加えた船舶航跡分析処理の一例を説明するフローチャートである。
【図８】拡張カルマンフィルタにオペレータ介入分析処理及び初期条件変更リトライ処理を加えた船舶航跡分析処理の一例を説明するフローチャートである。
【図９】時系列残差をグラフィック表示により表した例を説明する図である。
【図１０】観測精度が付加された時系列残差をグラフィック表示により表した例を説明する図である。
【符号の説明】
１  受動センサ
２  到来音波または電波観測装置
３  到来音波また  電波観測データベース
４  分析処理器
５  船舶航跡分析値表示器
６  初期条件演算器
７  観測行列演算器
８  端末部
９  時系列残差演算器
【図１】
2004239678.tif 000061
【図２】
2004239678.tif 000062
【図３】
2004239678.tif 000063
【図４】
2004239678.tif 000064
【図５】
2004239678.tif 000065
【図６】
2004239678.tif 000066
【図７】
2004239678.tif 000067
【図８】
2004239678.tif 000068
【図９】
2004239678.tif 000069
【図１０】
2004239678.tif 000070