説明

推定装置及びコンピュータプログラム及び推定方法

【課題】処理時間を低減しつつ目標航跡を精度良く抽出する。
【解決手段】航跡候補生成部140は、目標航跡を記述するパラメータのうちいくつかを固定した上でパラメータの候補を多数生成する。相関ゲート設定部150は、航跡候補生成部140が生成したパラメータ候補におけるいくつかの固定されたパラメータと目標の真のパラメータとの誤差を最大限考慮した相関ゲートを設定する。概探索部160は、航跡候補生成部140が生成したパラメータ候補と、観測データ管理部130が取り出した観測データとから、相関ゲート設定部150が設定した相関ゲートを用いてパラメータ候補の尤度を計算し、尤度が最大となるパラメータ候補を抽出する。詳細探索部170は、概探索部160が抽出したパラメータ候補を初期値として、尤度が大きくなる方向にパラメータを更新し、相関ゲートが小さくなるように更新して、尤度が最大となるパラメータを探索する。

【発明の詳細な説明】
【技術分野】
【0001】
この発明は、入力データの時系列から目標航跡を抽出することにより目標の航跡を推定し目標を追尾する推定装置(追尾装置)に関する。
【背景技術】
【0002】
センサによって取得された、多数の不要信号を含む入力データの時系列から、目標航跡を抽出する方式として、ML−PDA(Maximum Likelihood Probabilistic Data Association)がある。ML−PDAでは、目標の位置や速度などをあらかじめ仮定し、最尤法により目標航跡を抽出する。
ML−PDAにおいて、パラメータ探索を効率化する方式として、格子点探索方式による概探索と、概探索の結果を初期値とする準Newton法等の最適化アルゴリズムによる詳細探索とを組み合わせる方式がある。
更に、概探索を効率化する方式として、最初は粗く格子点を生成し、評価関数が予め設定した閾値以上となった格子点を含む領域において再び格子点を生成するということを繰り返す方式がある。
また、航跡の位置ベクトルと観測値ベクトルとの差が小さいほど評価関数値が増大することを利用して、格子点の位置成分を観測値ベクトルに置き換えることにより、速度や加速度など、直接観測できない高次のパラメータについてのみ格子点を生成する方式がある。
【先行技術文献】
【特許文献】
【0003】
【特許文献1】特開2008−96137号公報
【非特許文献】
【0004】
【非特許文献1】T.Kirubarajan他「Passive ranging of a low observable ballistic missile in a gravitational field」IEEE Transactions on Aerospace and Electronic Systems,vol.37,no.2,481〜494ページ、2001年。
【非特許文献2】W.R.Blanding他「Directed subspace search ML−PDA with application to active sonar tracking」IEEE Transactions on Aerospace and Electronic Systems,vol.44,no.1,201〜216ページ、2008年。
【発明の概要】
【発明が解決しようとする課題】
【0005】
ML−PDAを用いて目標を追尾する場合、例えば目標が平面内を等速直線運動しているとすると、推定すべきパラメータ(状態量)は、目標の初期位置ベクトルと速度ベクトルであり、各々2次元の計4次元空間から尤度が最大となるパラメータを探索する。一方、目標が平面内を等加速度運動しているとすると、推定すべきパラメータに加速度ベクトルが加わるため、6次元空間からパラメータを探索する必要がある。目標が空間内を運動しているとすると、等速直線運動の場合で6次元空間、等加速度運動の場合で9次元空間からパラメータを探索する必要がある。
したがって、推定すべきパラメータの数が増えると、追尾に必要な計算量が多くなり、処理に時間がかかる。
【0006】
格子点探索方式による概探索と、概探索の結果を初期値とする最適化アルゴリズムによる詳細探索とを組み合わせる方式において、概探索の処理時間は、生成する格子点数に依存する。また、詳細探索は、一般に概探索に比べて少ない処理で済む。探索精度を向上させるためには、概探索において格子点の数を十分多く生成する必要がある。格子点を増やすと、計算量が多くなるので、処理に時間がかかる。
【0007】
粗く格子点を生成し、段階的に格子点生成領域を絞り込む格子点探索方式を用いる場合、概探索における推定精度を見積もることはできるが、探索空間の次元が減るわけではないため、探索の処理に時間がかかる。
また、格子点の位置成分を観測値ベクトルに置き換える格子点探索方式を用いる場合、位置成分の候補数が観測値数に依存するため、観測値数が少ない場合は格子点数を減らすことができるが、高次のパラメータについての格子点数が減るわけではないため、高次のパラメータ数が多い場合は探索の処理に時間がかかる。
【0008】
この発明は、例えば上記のような課題を解決するためになされたものであり、処理時間を低減しつつ目標航跡を精度良く抽出することを目的とする。
【課題を解決するための手段】
【0009】
この発明にかかる推定装置は、データを処理する処理装置と、候補算出部と、相関部と、尤度算出部と、選択部と、推定値更新部と、推定値出力部とを有し、
上記候補算出部は、目標の複数の状態量のうち所定の第一の状態量を除く他の状態量について、上記処理装置を用いて、上記他の状態量の推定値の複数の候補値を算出し、
上記相関部は、上記候補算出部が算出した上記他の状態量の推定値の複数の候補値それぞれについて、上記処理装置を用いて、上記第一の状態量が所定の推定値であり、上記第一の状態量が上記推定値であると仮定することによって生じるモデル化誤差が所定のモデル化誤差上限値以下であり、上記他の状態量が上記候補値であるとの仮定に基づいて、上記目標を観測した観測値を含む複数の観測値のなかから、上記目標を観測した観測値を抽出し、
上記尤度算出部は、上記候補算出部が算出した上記他の状態量の推定値の複数の候補値それぞれについて、上記処理装置を用いて、上記相関部が仮定したモデル化誤差上限値と、上記相関部が抽出した観測値とに基づいて、上記他の状態量が上記候補値である尤度を算出し、
上記候補選択部は、上記処理装置を用いて、上記尤度算出部が算出した尤度に基づいて、上記候補算出部が算出した複数の候補値のなかから、上記他の状態量の推定値を選択し、
上記推定値更新部は、上記処理装置を用いて、上記相関部が仮定した第一の状態量の推定値と上記候補選択部が選択した他の状態量の推定値とを出発点として、上記目標の複数の状態量が推定値である尤度が高くなるよう、上記複数の状態量の推定値を更新し、
上記推定値出力部は、上記処理装置を用いて、上記推定値更新部が更新した上記複数の状態量の推定値を出力することを特徴とする。
【発明の効果】
【0010】
この発明にかかる推定装置によれば、処理時間を低減しつつ目標航跡を精度良く抽出することができる。
【図面の簡単な説明】
【0011】
【図1】実施の形態1における追尾システム800の全体構成の一例を示すシステム構成図。
【図2】実施の形態1における追尾装置100のハードウェア資源の一例を示すハードウェア構成図。
【図3】実施の形態1における追尾装置100の機能ブロックの構成の一例を示すブロック構成図。
【図4】実施の形態1における追尾処理S600の流れの一例を示すフローチャート図。
【図5】目標の航跡の一例を示す図。
【図6】実施の形態1における概探索ステップS620における相関ゲートと、得られる航跡との関係の一例を示す図。
【図7】実施の形態1における詳細探索ステップS630における更新前の相関ゲート及び航跡と更新後の相関ゲート及び航跡との関係の一例を示す図。
【図8】実施の形態1における詳細探索ステップS630により最終的に得られる相関ゲートと航跡との関係の一例を示す図。
【図9】実施の形態3における量子化誤差を説明するための図。
【図10】実施の形態6における追尾処理S600の流れの一例を示すフローチャート図。
【図11】実施の形態7における追尾装置100の機能ブロックの構成の一例を示すブロック構成図。
【図12】実施の形態7における追尾処理S600の流れの一例を示すフローチャート図。
【発明を実施するための形態】
【0012】
実施の形態1.
実施の形態1について、図1〜図8を用いて説明する。
【0013】
図1は、この実施の形態における追尾システム800の全体構成の一例を示すシステム構成図である。
追尾システム800(推定システム)は、目標700を観測し、観測した目標700の航跡を推定して、目標700を追尾する。追尾システム800は、センサ810と、追尾装置100とを有する。
センサ810(観測装置)は、目標を観測し、観測した観測結果を表わす信号を出力する。センサ810は、連続して、もしくは繰り返し、観測を行う。なお、センサ810は、追尾装置100の一部であってもよい。
追尾装置100(推定装置、航跡推定装置)は、センサ810が観測した観測結果に基づいて、目標700の航跡を推定する。
【0014】
図2は、この実施の形態における追尾装置100のハードウェア資源の一例を示すハードウェア構成図である。
追尾装置100は、例えば、コンピュータである。追尾装置100は、例えば、処理装置911と、入力装置912と、出力装置913と、記憶装置914とを有する。
処理装置911は、記憶装置914が記憶したプログラムを実行することにより、データを処理し、追尾装置100全体を制御する。
記憶装置914は、処理装置911が実行するプログラムや、処理装置911が処理するデータなどを記憶する。記憶装置914は、例えば、揮発性メモリ、不揮発性メモリ、半導体メモリ、磁気ディスク装置、光ディスク装置などである。
入力装置912は、追尾装置100の外部から信号を入力し、入力した信号を処理装置911が処理できるデータに変換する。入力装置912が変換したデータは、処理装置911が直接処理する構成でもよいし、一時的に記憶装置914が記憶する構成でもよい。入力装置912は、例えば、アナログ信号をデジタルデータに変換するアナログデジタル変換装置、他の装置が送信した信号を受信して復調する受信装置、操作者の操作を入力するキーボードやマウスなどの操作入力装置などである。
出力装置913は、処理装置911が処理したデータや記憶装置914が記憶したデータなどを変換して、追尾装置100の外部に出力する。出力装置913は、例えば、デジタルデータをアナログ信号に変換するデジタルアナログ変換装置、他の装置に対する信号を送信する送信装置、操作者に対して情報を提示する表示装置・音声出力装置などである。
以下に説明する追尾装置100の機能ブロックは、記憶装置914が記憶したプログラムを処理装置911が実行することにより実現される。なお、これは一例であり、追尾装置100の機能ブロックは、他の方式で実現したものであってもよい。追尾装置100の機能ブロックは、例えば、アナログ回路やデジタル回路などの電気的構成で実現したものであってもよいし、機械的構成など他の構成により実現したものであってもよい。また、追尾装置100は、物理的に1つの装置である必要はなく、物理的に複数の装置が、それぞれ追尾装置100の機能ブロックのうちのいくつかを実現し、全体として追尾装置100として機能する構成であってもよい。
【0015】
図3は、この実施の形態における追尾装置100の機能ブロックの構成の一例を示すブロック構成図である。
追尾装置100は、例えば、観測データ抽出部110と、観測データ蓄積部120と、観測データ管理部130と、航跡候補生成部140と、相関ゲート設定部150と、概探索部160と、詳細探索部170とを有する。
【0016】
観測データ抽出部110は、入力装置912を用いて、センサ810が出力した信号(入力データ)を入力する。観測データ抽出部110は、処理装置911を用いて、入力した信号から、目標を観測した観測値(観測データ)を抽出する。なお、センサ810が出力する信号は、ノイズなどの影響による不要信号を含むので、観測データ抽出部110が抽出する観測値は、必ずしも目標を観測した観測値ではない可能性がある。
観測データ蓄積部120は、記憶装置914を用いて、観測データ抽出部110が抽出した観測値を、センサ810が観測した観測時刻(フレーム番号)と対応づけて蓄積して記憶する。
観測データ管理部130は、処理装置911を用いて、目標の航跡の推定に使用する観測値の観測時刻の範囲を決定する。観測データ管理部130は、観測データ蓄積部120が記憶した観測値のなかから、観測時刻が決定した範囲内である観測値(観測データ時系列)を抽出する。
航跡候補生成部140は、処理装置911を用いて、目標の航跡の候補を生成する。航跡候補生成部140は、複数の候補(航跡候補群)を生成する。
相関ゲート設定部150は、処理装置911を用いて、相関ゲート行列を設定・更新する。
概探索部160は、処理装置911を用いて、観測データ管理部130が抽出した観測値と、相関ゲート設定部150が設定した相関ゲート行列とに基づいて、航跡候補生成部140が生成した候補のなかから、実際の目標の航跡に近い航跡(航跡候補)を抽出する。
詳細探索部170は、処理装置911を用いて、観測データ管理部130が抽出した観測値と、相関ゲート設定部150が設定した相関ゲート行列とに基づいて、概探索部160が抽出した航跡から出発して、更に実際の目標の航跡に近い航跡を算出する。詳細探索部170は、相関ゲート設定部150が更新した相関ゲート行列に基づいて、更新した航跡から更に実際の目標の航跡に近い航跡を算出する。詳細探索部170は、航跡が収束するまでこれを繰り返す。詳細探索部170は、処理装置911を用いて、算出した航跡(航跡推定値)を出力する。
【0017】
図4は、この実施の形態における追尾処理S600の流れの一例を示すフローチャート図である。
追尾処理S600(推定方法、航跡推定方法)において、追尾装置100は、センサ810が観測した観測結果に基づいて、目標を追尾し、目標の航跡を推定する。追尾処理S600は、例えば、データ入力ステップS611と、観測データ抽出・蓄積ステップS612と、使用フレーム選択ステップS613と、相関ゲート設定ステップS614と、航跡候補設定ステップS615と、概探索ステップS620と、詳細探索ステップS630と、相関ゲート更新ステップS643と、ゲーティングステップS644と、尤度計算ステップS645と、航跡推定値出力ステップS651とを有する。追尾装置100は、追尾処理S600をデータ入力ステップS611から開始する。
【0018】
まず、データ入力ステップS611において、センサ810は、データを取得する。センサ810は、例えば、レーダ、赤外線センサ、光学カメラなどである。
次に、観測データ抽出・蓄積ステップS612において、観測データ抽出部110は、入力データから観測データを抽出する。観測データ蓄積部120は、観測データを観測された時刻ごとに蓄積する。観測データは、例えば、2次元画像の位置座標である。観測データを抽出する方式としては、例えば、輝度値が予め設定した閾値を超えた画素を抽出する方式、輝度値が上位の画素を一定数個抽出する方式などがある。
【0019】
次に、使用フレーム選択ステップS613において、観測データ管理部130は、航跡の推定に使用する時刻フレーム(観測時刻)を任意に選択する。選択された時刻フレームをt、t、…、tとする。ただし、u<vなる任意のuとvについてt<tであるとする。
観測データ管理部130は、観測データ蓄積部120が蓄積した観測データ群から前記時刻フレームにおける式(1)で表わされる観測データ時系列Zを取り出す。
【数11】

ただし、Zは、時刻tにおける観測データである。時刻tにおける観測データZは、M個の観測値ベクトルから構成される。また、zk,jは、時刻tにおける観測値ベクトルである。xk,j,yk,jは、それぞれx座標、y座標を表わす。
【0020】
次に、相関ゲート設定ステップS614において、相関ゲート設定部150は、式(2)で表わされる相関ゲート行列Sを設定する。
【数12】

ただし、Rは、観測雑音共分散行列である。観測雑音共分散行列Rは、正定値対称行列であり、例えば、観測値ベクトルの観測精度の2乗値が対角成分に並んだ行列である。
は、モデル化誤差共分散行列である。モデル化誤差共分散行列Qは、半正定値対称行列であり、例えば、式(3)で表わされる、想定される最大加速度誤差ベクトルamax=(amax,amax)を加速度として目標が移動した場合の目標の変位量の2乗値が対角成分に並んだ行列である。
【数13】

ただし、tは、目標の位置、速度、加速度の初期値を定める基準時刻である。基準時刻tは、例えば、観測開始時刻tや、最新時刻tなどである。
また、最大加速度誤差ベクトルamax=(amax,amax)の値は事前に設定しておく。最大加速度誤差ベクトルamaxの値は、例えば、推定したい目標の種類に応じて設定した値を用いる。
【0021】
次に、航跡候補設定ステップS615において、航跡候補生成部140は、概探索のための航跡候補群を作成する。ここで、航跡は、基準時刻tにおける目標の位置ベクトル、速度ベクトル、加速度ベクトルを成分に持つ状態ベクトルxとして、式(4)で記述される。
【数14】

ただし、(x,y)は、基準時刻tにおける目標の位置ベクトルである。(xドット,yドット)は、基準時刻tにおける目標の速度ベクトルである。(xツードット,yツードット)は、基準時刻tにおける目標の加速度ベクトルである。目標の位置ベクトル、速度ベクトル、及び加速度ベクトルは、いずれも、xy座標で表わす。また、上付きのTは行列の転置を表わす。
【0022】
航跡候補生成部140は、加速度ベクトルの候補を、ある固定値afix=(afix,afix)とした上で、あらかじめ定められた範囲における航跡候補、すなわち状態ベクトル候補を複数個生成する。状態ベクトル候補を生成する範囲は、例えば、位置成分は、センサの探知可能領域の範囲、速度成分は、目標の想定される移動方向に沿った値の範囲である。ここで、航跡候補生成部140は、位置成分、速度成分が前記範囲内で等間隔になるように生成する構成であってもよいし、乱数を用いて不等間隔に生成する構成であってもよい。また、航跡候補生成部140は、位置成分の候補については、観測値ベクトルを使う構成であってもよい。
また、加速度ベクトル候補の固定値afixは、事前に設定しておく。加速度ベクトル候補の固定値afixは、例えば、0[m/s]である。あるいは、航跡候補生成部140は、前回の推定結果の加速度ベクトル推定値を、加速度ベクトル候補の固定値afixとして用いる構成であってもよい。
【0023】
次に、概探索ステップS620において、概探索部160は、概探索をする。概探索ステップS620は、状態遷移ステップS623と、ゲーティングステップS624と、尤度計算ステップS625と、最大尤度比航跡候補抽出ステップS627とを有する。概探索部160は、航跡候補生成部140が生成した航跡候補それぞれについて、状態遷移ステップS623、ゲーティングステップS624、尤度計算ステップS625を実行する。
【0024】
概探索ステップS620の状態遷移ステップS623において、概探索部160は、航跡候補ごとに、式(5)にしたがって、各時刻tにおける状態ベクトルxを計算する。
【数15】

ただし、Iは2行2列の単位行列である。0は、2行2列の零行列である。
【0025】
次に、概探索ステップS620のゲーティングステップS624において、概探索部160は、航跡候補ごとに、各時刻tにおける状態ベクトルxの位置成分と、相関ゲート設定部150が設定した相関ゲート行列Sに基づいて、楕円形状の相関ゲートGを、式(6)を使って設定する。
【数16】

ただし、Hは、状態ベクトルから位置成分を取り出す観測行列である。また、dは、事前に設定するゲートサイズパラメータである。ゲートサイズパラメータdの値は、自由度2のカイ2乗分布より定まる正定数である。上付きのTは、転置を表わす。上付きの−1は、逆行列を表わす。
【0026】
次に、概探索ステップS620の尤度計算ステップS625において、概探索部160は、航跡候補ごとに、各時刻tにおける状態ベクトルxと、相関ゲートGに含まれる観測値ベクトルzk,jとに基づいて、式(7)を使って尤度φを計算する。
【数17】

ただし、Pは、事前に設定する探知確率パラメータである。Pは、0<P<1なる定数である。また、βは、事前に設定する不要信号密度パラメータである。βは、正定数である。また、mは時刻tにおける観測データZを構成するM個の観測値ベクトルのうち、式(6)の相関ゲートGに含まれる観測値ベクトルの個数である。また、zk,jは、前記観測データZを構成する観測値ベクトルのうち、式(6)の相関ゲートGに含まれるものである。lnは、自然対数(ネイピア数を底とする対数関数)を表わす。detは、行列式を表わす。expは、ネイピア数を底とする指数関数を表わす。
【0027】
なお、式(6)及び式(7)は、目標の位置観測誤差の分布が正規分布で与えられる場合についての相関ゲート及び尤度の計算式の例である。観測誤差の分布が正規分布以外である場合などにおいては、相関ゲートは楕円形状以外の形状である構成であってもよい。また、尤度の計算式は任意に設定する構成であってもよい。例えば、尤度の計算式は、観測値ベクトルと状態ベクトルの位置成分との距離の逆2乗値を用いる構成であってもよい。また、尤度の計算式は、観測データの振幅情報や輝度情報に基づくものを用いる構成であってもよい。
【0028】
概探索部160は、状態遷移ステップS623、ゲーティングステップS624、尤度計算ステップS625を繰り返し、航跡候補生成部140が生成したすべての航跡候補を処理する。
【0029】
次に、概探索ステップS620の最大尤度比航跡候補抽出ステップS627において、概探索部160は、航跡候補の中から、尤度φが最大となるものを抽出する。そして、詳細探索部170は、最大尤度となる航跡候補(状態ベクトル候補)を、状態ベクトルの初期値として入力する。
【0030】
次に、詳細探索ステップS630において、詳細探索部170は、詳細探索をする。詳細探索ステップS630は、探索方向計算ステップS631と、航跡更新ステップS632と、状態遷移ステップS633と、ゲーティングステップS634と、尤度計算ステップS635と、探索収束条件分岐ステップS636とを有する。
【0031】
詳細探索ステップS630の探索方向計算ステップS631において、詳細探索部170は、前回の探索における状態ベクトルを起点として、より尤度φが大きくなる方向ベクトルを計算する。探索方向の計算方式は、例えば、最急降下法、Newton法、準Newton法、共役勾配法などがある。
【0032】
次に、詳細探索ステップS630の航跡更新ステップS632において、詳細探索部170は、前回の探索における状態ベクトルを、探索方向計算ステップS631で計算した方向ベクトルを用いて、前回の探索よりも尤度の大きい状態ベクトルに更新する。
【0033】
次に、詳細探索ステップS630の状態遷移ステップS633、ゲーティングステップS634及び尤度計算ステップS635において、詳細探索部170は、状態遷移、ゲーティング、尤度計算を行う。状態遷移ステップS633、ゲーティングステップS634及び尤度計算ステップS635における詳細探索部170の動作は、概探索の対応するステップ(状態遷移ステップS623、ゲーティングステップS624及び尤度計算ステップS625)における概探索部160の動作と同一である。
【0034】
そして、詳細探索ステップS630の探索収束条件分岐ステップS636において、詳細探索部170は、探索の収束を判定する。探索の収束判定方式には、例えば、尤度勾配ベクトルのノルムが十分小さい所定の閾値より小さくなった場合に収束したと判定する方式、前回の探索の尤度と今回の探索の尤度との差が十分小さい所定の閾値より小さくなった場合に収束したと判定する方式、探索収束条件分岐ステップS636を通過した回数があらかじめ定めた上限値を超えた場合に収束したと判定する方式などがある。
探索が収束したと判定しなかった場合、詳細探索部170は、相関ゲート更新ステップS643へ処理を進める。
探索が収束したと判定した場合、詳細探索部170は、航跡推定値出力ステップS651へ処理を進める。
【0035】
相関ゲート更新ステップS643において、相関ゲート設定部150は、相関ゲート行列Sを更新する。例えば、相関ゲート設定部150は、探索によって得られた加速度ベクトル推定値aest=(aest,aest)、加速度ベクトル固定値afix=(afix,afix)、及び想定される最大加速度誤差ベクトルamax=(amax,amax)を用いて、式(8)のように、相関ゲート行列Sを更新する。
【数18】

【0036】
最大加速度誤差ベクトルamaxは、加速度ベクトル固定値afixと実際の目標加速度ベクトルとの差の最大許容量を表すものである。また、加速度ベクトル固定値afixと加速度ベクトル推定値aestとの差は、加速度ベクトル固定値afixが実際の目標加速度ベクトルに近づくための補正量を表すものである。したがって、前記補正量の分だけ加速度誤差が減少したと考えられる。そこで、前記最大許容量から前記補正量を減算することにより、モデル化誤差共分散行列Qの大きさを小さくする。そして、更新後の加速度誤差最大許容量(aハットmax,aハットmax)を新たな最大加速度誤差ベクトルaハットmaxとして設定する。ただし、式(9)の右辺が負値となる場合、相関ゲート設定部150は、更新後の加速度誤差最大許容量aハットmaxあるいはaハットmaxを、0[m/s]とする。
【0037】
そして、ゲーティングステップS644及び尤度計算ステップS645において、詳細探索部170は、更新された相関ゲート行列Sを用いて、ゲーティング及び尤度計算を行う。ゲーティングステップS644及び尤度計算ステップS645における詳細探索部170の動作は、詳細探索ステップS630のゲーティングステップS634及び尤度計算ステップS635における動作と同様である。
そして、詳細探索部170は、再び、詳細探索ステップS630に処理を戻し、探索方向計算ステップS631以降の動作を行う。
【0038】
探索収束条件分岐ステップS636で探索が収束したと判定した場合、航跡推定値出力ステップS651において、詳細探索部170は、航跡推定値として探索した状態ベクトルを出力する。
【0039】
図5は、目標の航跡の一例を示す図である。
例えば、目標が一定の加速度で等加速度運動をしているものとする。航跡候補生成部140が生成した航跡の候補における加速度は、実際の目標の加速度と異なっている。
【0040】
図6は、この実施の形態における概探索ステップS620における相関ゲートと、得られる航跡との関係の一例を示す図である。
この例において、基準時刻tは、時刻tである。観測時刻と基準時刻tとの差が大きくなるにつれて、相関ゲートが大きくなる。これにより、加速度の固定値と実際の目標の加速度との間の差が大きい場合でも、観測値が相関ゲートに入るので、どの航跡の候補が実際の航跡に最も近いかを判定することができる。
【0041】
図7は、この実施の形態における詳細探索ステップS630における更新前の相関ゲート及び航跡と更新後の相関ゲート及び航跡との関係の一例を示す図である。
詳細探索部170は、基準時刻における目標の位置、速度及び加速度の推定値を更新することにより、推定航跡を実際の目標の航跡に近づける。これにより、相関ゲートを小さくしても、観測値が相関ゲートに入るようになる。
相関ゲート設定部150が小さくした相関ゲートに基づいて、詳細探索部170が航跡を更新することにより、更に、推定航跡を実際の目標の航跡に近づけることができる。
【0042】
図8は、この実施の形態における詳細探索ステップS630により最終的に得られる相関ゲートと航跡との関係の一例を示す図である。
航跡の更新と相関ゲートの更新とを、航跡が収束するまで繰り返すことにより、最終的に、実際の目標の航跡にほぼ一致する航跡を推定することができる。
【0043】
この実施の形態における追尾装置100は、入力データの時系列から目標航跡を抽出することによって目標を追尾する。
追尾装置100は、観測データ抽出部110と、観測データ蓄積部120と、観測データ管理部130と、航跡候補生成部140と、相関ゲート設定部150と、概探索部160と、詳細探索部170とを有する。
前記観測データ抽出部110は、時刻ごとに入力データから観測データを抽出する。
前記観測データ蓄積部120は、前記観測データを時刻ごとに蓄積する。
前記観測データ管理部130は、指定された複数の時刻の観測データを前記観測データ蓄積部120から取り出す。
前記航跡候補生成部140は、目標航跡を記述するパラメータのうちいくつかを固定した上でパラメータの候補を多数生成する。
前記相関ゲート設定部150は、前記航跡候補生成部140で生成されたパラメータ候補におけるいくつかの固定されたパラメータと目標の真のパラメータとの誤差を最大限考慮した相関ゲートを設定する。
前記概探索部160は、前記航跡候補生成部140で生成されたパラメータ候補と、前記観測データ管理部130で取り出された観測データとから、前記相関ゲート設定部150で設定された相関ゲートを用いて前記パラメータ候補の尤度を計算するとともに、前記尤度が最大となるパラメータ候補を抽出する。
前記詳細探索部170は、前記概探索部160で抽出された尤度が最大となるパラメータ候補を初期値として、前記尤度が大きくなる方向にパラメータを更新しながら、前記相関ゲートをより小さくなるように更新して、前記尤度が最大となるパラメータを探索することにより、目標航跡を記述するパラメータの推定値を出力する。
【0044】
前記詳細探索部170は、前記尤度が大きくなる方向にパラメータを更新するたびに、前記相関ゲートをより小さくなるように更新することを繰り返す。
【0045】
前記航跡候補生成部140は、航跡を記述するパラメータとして目標の位置、速度、加速度を設定するとともに、固定パラメータとして加速度を選択する。
前記相関ゲート設定部150は、前記航跡候補生成部140で固定された加速度と目標の真の加速度との誤差を最大限考慮した相関ゲートを設定する。
【0046】
前記相関ゲート設定部150は、目標の観測誤差共分散行列と、目標の真の加速度との最大誤差許容量を考慮した共分散行列を用いた式(2)のような相関ゲート行列を用いて相関ゲートを規定する。
【0047】
前記相関ゲート設定部150は、目標の真の加速度との最大誤差許容量を考慮した共分散行列として式(3)を用いる。
【0048】
前記概探索部160及び前記詳細探索部170は、前記相関ゲート行列を用いて相関ゲートを式(6)のように設定し、パラメータの尤度を式(7)のように計算する。
【0049】
前記詳細探索部170は、前記相関ゲートを式(8)のように更新する。
【0050】
この実施の形態における追尾装置100は、設定した加速度固定値と実際の目標加速度との差をモデル化誤差として表現し、モデル化誤差を考慮した相関ゲートを設定して概探索を行う。これにより、設定した加速度固定値が実際と異なっていても目標観測値を相関ゲート内に捕捉することができる。
さらに、この実施の形態における追尾装置100は、概探索においては位置ベクトルと速度ベクトルのみを探索し、加速度ベクトルのような高次のパラメータについては探索を行わない。これにより、加速度ベクトルも探索する場合に比べて航跡候補数を減らすことができ、処理時間を削減することができる。
さらに、この実施の形態における追尾装置100は、詳細探索においては加速度誤差を表すモデル化誤差共分散行列の大きさを、推定結果の加速度を基に計算した補正量を用いて段階的に小さくしながら、加速度ベクトルのような高次のパラメータも含めて探索を行う。これにより、相関ゲートを目標観測値の存在する近傍に絞り込みつつ、実際の目標の位置ベクトル、速度ベクトル、及び加速度ベクトルを精度良く推定することが可能となる。
【0051】
なお、この実施の形態では、平面内を等加速度運動する目標を追尾する構成について説明したが、追尾装置100は、他の運動モデルに基づいて目標を追尾する構成であってもよい。例えば、目標の位置、速度、加速度などを二次元ベクトルで表わすのではなく、三次元ベクトルで表わす構成であってもよい。あるいは、目標が、等加速度運動ではなく、例えば等速円運動など他の運動をしていると仮定して、目標を追尾する構成であってもよい。
【0052】
実施の形態2.
実施の形態2について、説明する。
なお、実施の形態1と共通する部分については、説明を省略する。
【0053】
実施の形態1における相関ゲート設定部150は、式(8)のように加速度ベクトル推定結果を用いて相関ゲート行列Sを更新する。これに対し、この実施の形態における相関ゲート設定部150は、加速度ベクトル推定結果を用いずに、相関ゲート行列Sを更新する。例えば、相関ゲート設定部150は、想定される最大加速度誤差ベクトルamax=(amax,amax)を用いて、式(10)のように、相関ゲート行列Sを更新する。
【数19】

ただし、γは、事前に設定する減衰係数パラメータである。γは、0<γ<1なる定数である。また、nは、詳細探索の反復回数(探索収束条件分岐ステップS636を通過した回数)である。
【0054】
この実施の形態における追尾装置100において、前記詳細探索部170は、前記相関ゲートを式(10)のように更新する。
【0055】
この実施の形態における追尾装置100は、式(10)による相関ゲート行列の更新をして、加速度ベクトルの推定結果に関わらず、一定率でモデル化誤差共分散行列Qの大きさを減少させる。これにより、詳細探索における相関ゲートの絞り込みを容易にすることができる。
【0056】
実施の形態3.
実施の形態3について、図9を用いて、説明する。
なお、実施の形態1及び実施の形態2と共通する部分については、説明を省略する。
【0057】
実施の形態1における概探索部160は、式(2)のような相関ゲート行列を設定する。これに対し、この実施の形態における概探索部160は、格子点探索方式によって概探索を行う場合、探索空間を量子化することと同等であることから、状態ベクトル候補には観測誤差の他に量子化誤差が含まれると考え、式(11)のような相関ゲート行列を設定する。
【数20】

ただし、Q’は量子化誤差共分散行列である。量子化誤差共分散行列Q’は、一定の半正定値対称行列である。
【0058】
図9は、この実施の形態における量子化誤差を説明するための図である。
例えば、航跡候補生成部140が生成する航跡の候補において、所定の間隔を有する複数の値を、位置や速度の推定値の候補とする場合、尤度が最も高い航跡の候補を選択しても、実際の目標の航跡の位置や速度との間にずれが生じる。これが量子化誤差である。
したがって、位置や速度の推定値の候補の間の間隔から、量子化誤差の絶対値の最大値を求めることができる。
概探索部160は、このようにして算出した量子化誤差の絶対値の最大値に基づいて、量子化誤差協分散行列Q’を算出する。
【0059】
詳細探索で探索を進めるごとに、量子化誤差の影響は小さくなる。そこで、相関ゲート設定部150は、例えば、式(10)のモデル化誤差共分散行列の更新式と同様に、量子化誤差共分散行列に詳細探索の反復回数に応じて一定の減衰係数を乗ずる。相関ゲート設定部150は、例えば、次の式により、相関ゲート行列Sを更新する。
【数21】

ただし、ζは、事前に設定する減衰係数パラメータである。ζは、0<ζ<1なる定数である。また、nは、詳細探索の反復回数(探索収束条件分岐ステップS636を通過した回数)である。
【0060】
あるいは、相関ゲート設定部150は、例えば、次の式により、相関ゲート行列Sを更新する。
【数22】

【0061】
この実施の形態における追尾装置100において、前記相関ゲート設定部150は、目標の観測誤差共分散行列と、目標の真の加速度との最大誤差許容量を考慮した共分散行列に加えて、さらに概探索におけるパラメータ候補の量子化誤差を考慮した誤差共分散行列を用いた式(11)のような相関ゲート行列を用いて相関ゲートを規定する。
【0062】
この実施の形態における追尾装置100は、概探索において量子化誤差の影響を考慮することにより、目標航跡を精度良く抽出することができる。
【0063】
実施の形態4.
実施の形態4について、説明する。
なお、実施の形態1〜実施の形態3と共通する部分については、説明を省略する。
【0064】
実施の形態1における相関ゲート設定部150は、最大加速度誤差ベクトルamax=(amax,amax)を事前に設定する。これに対し、この実施の形態における相関ゲート設定部150は、相関ゲートの大きさに上限値を事前に設定することにより、最大加速度誤差ベクトルamaxの値を計算する。例えば、式(6)のように相関ゲートGを設定する場合、相関ゲートGの大きさVは、式(12)で計算することができる。
【数23】

ただしdetは、行列式を表わす。相関ゲート設定部150は、このVが事前に設定した上限値を超えないような最大加速度誤差ベクトルamaxの値を求める。
【0065】
この実施の形態における追尾装置100において、前記相関ゲート設定部150は、前記相関ゲートの大きさに上限値を設定することにより、目標の真の加速度との最大誤差許容量を計算する。
【0066】
これにより、最大加速度誤差ベクトルの値を相関ゲートの大きさから直感的に設定することが可能となる。
【0067】
実施の形態5.
実施の形態5について、説明する。
なお、実施の形態1〜実施の形態4と共通する部分については、説明を省略する。
【0068】
実施の形態1における観測データ管理部130は、使用フレーム選択ステップS613において、航跡の推定に使用する時刻フレームを任意に選択する。これに対し、この実施の形態における観測データ管理部130は、使用する時刻フレームの長さに上限値を設定した上で、上限値を超えないように時刻フレームを選択する。
【0069】
あるいは、観測データ蓄積部120は、あらかじめ相関ゲートの大きさの上限を定めておき、相関ゲートの大きさが、あらかじめ定めた上限を超えないよう、時刻フレームの長さを設定する構成であってもよい。
実施の形態4で述べたように、相関ゲートの大きさは、式(12)を使って計算することができる。実施の形態4では、最大加速度誤差ベクトルを変えることにより相関ゲートの大きさを変える構成を説明したが、相関ゲートの大きさは、基準時刻からの経過時間によっても変化する。例えば、観測データ蓄積部120は、最大加速度誤差ベクトルをあらかじめ定めた値に設定したうえで、相関ゲートの大きさが上限を超えないよう、基準時刻からの経過時間の絶対値の最大値を算出する。
【0070】
あるいは、観測データ蓄積部120は、例えば、相関ゲートに含まれる観測値ベクトルの数に上限値を設定し、相関ゲートに含まれる観測値ベクトルの数が上限値を超えないよう、時刻フレームの長さを設定する構成であってもよい。
例えば、観測データ蓄積部120は、基準時刻tから最も離れた観測時刻t(例えば、t=tとするとt=tである。)において、推定結果に用いた相関ゲートに含まれる観測値ベクトルの数mが事前に設定した閾値Mを超えたときの時刻フレーム数Nを上限値として設定する。
【0071】
この実施の形態における追尾装置100において、前記観測データ管理部130は、前記相関ゲートの大きさを基に、使用する時刻の長さに上限値を指定する。
【0072】
前記観測データ管理部130は、前記詳細探索部170が出力した結果における最も相関ゲートの大きい時刻について、相関ゲートに含まれる観測データの数が事前に設定した閾値を超えたときに使用していた時刻の長さを上限値とする。
【0073】
これにより、使用する時刻フレームが長くなることで、基準時刻から最も離れた時刻における相関ゲートが大きくなりすぎて、目標以外の観測値ベクトルが推定結果に与える影響が大きくなることを緩和することができる。
【0074】
実施の形態6.
実施の形態6について、図10を用いて説明する。
なお、実施の形態1〜実施の形態5と共通する部分については、説明を省略する。
【0075】
図10は、この実施の形態における追尾処理S600の流れの一例を示すフローチャート図である。
この実施の形態における追尾処理S600は、実施の形態1における追尾処理S600のステップに加えて、相関ゲート収束条件分岐ステップS641を有する。相関ゲート更新ステップS643は、相関ゲート収束条件分岐ステップS641の次段にある。
【0076】
ある相関ゲート行列Sを用いたときの状態ベクトルxの尤度をφ(x;S)と記述する。詳細探索ステップS630の航跡更新ステップS632における更新前の状態ベクトルをxa、更新後の状態ベクトルをxb、相関ゲート更新ステップS643における更新前の相関ゲート行列をSa、更新後の相関ゲート行列をSbとする。
航跡更新ステップS632において詳細探索部170は尤度が大きくなるように状態ベクトルを更新するので、φ(xa;Sa)≦φ(xb;Sa)が成り立つ。
しかし、相関ゲート更新ステップS643において相関ゲート設定部150は尤度を考慮せずに相関ゲート行列を更新するので、φ(xb;Sa)>φ(xb;Sb)となる可能性がある。このため、φ(xa;Sa)≦φ(xb;Sb)となるとは限らない。
したがって、航跡を更新するたびに相関ゲートを更新すると、尤度が大きくならず、正しい航跡を推定できない可能性がある。詳細探索ステップS630の探索方向計算ステップS631における方向ベクトルが計算できなくなり、真の目標航跡に向けた状態ベクトルの更新が行われない可能性がある。
【0077】
これを回避するため、詳細探索ステップS630の探索収束条件分岐ステップS636において、探索が収束していないと判定した場合、相関ゲート行列を更新せず、一定の相関ゲート行列を保ったまま、詳細探索部170は、探索方向計算ステップS631に処理を戻す。
探索が収束した場合、詳細探索部170は、相関ゲート収束条件分岐ステップS641に処理を進める。
【0078】
相関ゲート収束条件分岐ステップS641において、相関ゲート設定部150は、相関ゲートが収束したか否かを判定する。相関ゲートの収束判定方式には、例えば、前回の相関ゲート行列の行列式と今回の相関ゲート行列の行列式との差が十分小さい所定の閾値よりも小さくなった場合に収束したと判定する方式、相関ゲート収束条件分岐ステップS641を通過した回数があらかじめ定めた上限値を超えた場合に収束したと判定する方式などがある。
相関ゲートが収束したと判定した場合、相関ゲート設定部150は、航跡推定値出力ステップS651へ処理を進める。
相関ゲートが収束したと判定しなかった場合、相関ゲート設定部150は、相関ゲート更新ステップS643へ処理を進める。
【0079】
相関ゲート更新ステップS643において、相関ゲート設定部150は、例えば実施の形態1と同様、式(8)及び式(9)により、相関ゲート行列を更新する。なお、相関ゲート設定部150は、実施の形態2や実施の形態3で説明した方式など、他の方式で相関ゲート行列を更新する構成でもよい。
【0080】
ゲーティングステップS644及び尤度計算ステップS645において、詳細探索部170は、更新された相関ゲート行列を用いて、ゲーティング及び尤度計算を行う。ゲーティングステップS644及び尤度計算ステップS645における詳細探索部170の動作は、詳細探索ステップS630のゲーティングステップS634及び尤度計算ステップS635における動作と同様である。
そして、詳細探索部170は、再び、詳細探索ステップS630に処理を戻し、探索方向計算ステップS631以降の動作を行う。
【0081】
相関ゲート収束条件分岐ステップS641で相関ゲートが収束したと判定した場合、航跡推定値出力ステップS651において、詳細探索部170は、航跡推定値として探索した状態ベクトルを出力する。
【0082】
この実施の形態における追尾装置100において、前記詳細探索部170は、前記尤度が大きくなる方向にパラメータを更新することを繰り返し、尤度が十分大きくなった段階で前記相関ゲートをより小さくなるように更新して、再び尤度が大きくなる方向にパラメータを更新することを繰り返す。
【0083】
この実施の形態における追尾装置100は、相関ゲートの更新により尤度が減少する可能性を回避することができ、詳細探索における相関ゲートの絞り込みを容易にすることができる。
【0084】
実施の形態7.
実施の形態7について、図11〜図12を用いて説明する。
なお、実施の形態1〜実施の形態6と共通する部分については、同一の符号を付し、説明を省略する。
【0085】
図11は、この実施の形態における追尾装置100の機能ブロックの構成の一例を示すブロック構成図である。
追尾装置100(推定装置、航跡推定装置)は、観測値取得部111と、観測値記憶部121と、候補算出部141と、候補記憶部142と、誤差上限初期化部151と、誤差上限値記憶部152と、ゲート行列算出部153と、相関部161と、尤度算出部162と、選択部163と、推定値記憶部171と、推定値更新部172と、収束判定部173と、誤差上限値更新部174と、終了判定部175と、推定値出力部176とを有する。
【0086】
観測値取得部111は、入力装置912を用いて、センサ810が出力した信号を入力する。観測値取得部111は、処理装置911を用いて、入力した信号に基づいて、センサ810が目標を観測した観測値を取得する。観測値取得部111は、処理装置911を用いて、取得した観測値を表わす観測値データを生成する。
【0087】
例えば、センサ810は、レーダである。センサ810は、電波を放射し、放射した電波が目標に当たって反射した反射波を受信する。センサ810は、受信した反射波の強さを表わす信号を出力する。センサ810は、電波を放射する方向を一次元的に走査する。これにより、センサ810が出力した信号から、目標が存在する方向を判別できる。また、電波を放射してから反射波を受信するまでの経過時間に基づいて、センサ810から目標までの距離を判別できる。
例えば、観測値取得部111は、処理装置911を用いて、センサ810が出力した信号が表わす反射波の強さが所定の閾値より大きくなったとき、センサ810が目標を観測したと判定する。観測値取得部111は、処理装置911を用いて、反射波の強さが閾値より大きくなったタイミングなどに基づいて、センサ810が観測した目標が存在する方向や、センサ810から目標までの距離などを算出する。観測値取得部111は、処理装置911を用いて、算出した目標の方向と、目標までの距離とに基づいて、センサ810の位置を原点とした直交座標系における目標の位置(例えばX座標及びY座標)を算出する。観測値取得部111は、処理装置911を用いて、算出した目標の位置を観測値とする観測値データを生成する。
【0088】
あるいは、例えば、センサ810は、光学カメラである。センサ810は、目標が放射する可視光線や赤外線などの光を受光する。センサ810は、受光した光の強さを画素の輝度とする2次元画像を生成する。センサ810は、生成した画像を表わす信号を出力する。
例えば、観測値取得部111は、処理装置911を用いて、センサ810が出力した信号が表わす画像において、輝度が高い順に所定の数の画素を抽出し、抽出した画素に目標が観測されたと判定する。なお、輝度が高い画素が連続している場合、観測値取得部111は、1つの目標とみなし、そのなかで一番輝度が高い画素だけを抽出する構成であってもよい。観測値取得部111は、処理装置911を用いて、抽出した画素の画像内における位置(例えばX座標及びY座標)を算出する。観測値取得部111は、処理装置911を用いて、算出した画素の位置を観測値とする観測値データを生成する。
【0089】
なお、センサ810は、目標の位置を二次元的に観測する構成に限らず、目標の位置を三次元的に観測する構成であってもよいし、ドップラー効果などに基づいて目標の速度を観測する構成であってもよい。
【0090】
センサ810が出力する信号には、ノイズなどの影響による不要信号が含まれている。このため、観測値取得部111が抽出する観測値は、必ずしも目標を観測したものではない可能性がある。逆に、目標が存在しても、観測された信号が弱い場合は、ノイズなどに紛れて、観測値取得部111が観測値として抽出しない可能性がある。
【0091】
観測値記憶部121は、記憶装置914を用いて、観測値取得部111が生成した観測値データを蓄積して記憶する。センサ810は、連続して、もしくは繰り返し、観測を行うので、時間の経過に伴い、観測値記憶部121が記憶する観測値データは増えていく。観測値記憶部121は、センサ810が目標を観測した時刻がわかるように、観測値データを記憶する。例えば、観測値記憶部121は、処理装置911を用いて、観測値データが表わす観測値の観測時刻やフレーム番号を表わすデータを生成して観測値データに付加し、記憶装置914を用いて、生成したデータ付加した観測値データを記憶する。あるいは、観測値記憶部121は、処理装置911を用いて、観測値データが表わす観測値の観測時刻やフレーム番号ごとに観測値データを分類し、記憶装置914を用いて、分類した観測値データを記憶する。なお、観測値記憶部121は、観測時刻から所定時間以上経過した観測値データを破棄し、観測時刻からの経過時間が所定時間内である観測値データのみを保持する構成であってもよい。
【0092】
追尾装置100は、所定の運動モデルに基づいて目標を追尾する。例えば、目標の状態量として、時刻tにおける三次元空間内における位置を表わす位置ベクトルp(t)、速度を表わす速度ベクトルpドット(t)、加速度を表わす加速度ベクトルpツードット(t)の三つの三次元ベクトルを考える場合、目標の状態量ベクトルx(t)は、この三つの三次元ベクトルの成分を成分とする九次元ベクトルで表わすことができる。
【数24】

ただし、p(t)は、時刻tにおける目標の位置の三次元直交座標系におけるX座標である。p(t)は、時刻tにおける目標の位置の三次元直交座標系におけるY座標である。p(t)は、時刻tにおける目標の位置の三次元直交座標系におけるZ座標である。
(t)は、時刻tにおける目標の速度の三次元直交座標系におけるX方向成分である。v(t)は、時刻tにおける目標の速度の三次元直交座標系におけるY方向成分である。v(t)は、時刻tにおける目標の速度の三次元直交座標系におけるZ方向成分である。
(t)は、時刻tにおける目標の加速度の三次元直交座標系におけるX方向成分である。a(t)は、時刻tにおける目標の加速度の三次元直交座標系におけるY方向成分である。a(t)は、時刻tにおける目標の加速度の三次元直交座標系におけるZ方向成分である。
【0093】
目標が等加速度運動をしているという運動モデルを用いる場合、目標の航跡は、ある基準時刻tにおける状態量ベクトルx(t)で表わすことができる。このとき、基準時刻tから時間Δtが経過した時刻t=t+Δtにおける目標の状態量ベクトルx(t+Δt)は、次の式で求めることができる。
【数25】

ただし、F(Δt)は、時間Δtの経過に伴う状態量ベクトルの変化を表わす遷移行列である。遷移行列F(Δt)は、例えば、次の式で表わされる。
【数26】

【0094】
例えば、センサ810が光学カメラである場合、時刻tにおける観測値を表わす観測値ベクトルz(t)と、その時刻tにおける目標の状態量ベクトルx(t)との間には、例えば、次のような関係が成り立つ。
【数27】

ただし、p及びpは、目標が観測されたと観測値取得部111が判定した画素の二次元座標である。画素の座標は、センサ810が撮影した画像の中心を原点とする。また、状態量ベクトルにおける三次元座標系は、センサ810の撮影方向をZ軸とし、センサ810が撮影した画像上の横方向をX軸、縦方向をY軸とする。また、vは、観測誤差を表わす二次元ベクトル(観測誤差ベクトル)である。hは、観測関数を表わす。観測関数hは、例えば、次の式で表わされる。
【数28】

【0095】
候補算出部141は、処理装置911を用いて、目標の航跡の候補を算出する。候補算出部141は、処理装置911を用いて、複数の候補を算出する。例えば、候補算出部141は、処理装置911を用いて、目標の航跡を表わす状態量ベクトルの九つの成分それぞれについて、1つ以上の候補を算出し、それを組み合わせることにより、複数の航跡の候補を算出する。
候補記憶部142は、記憶装置914を用いて、候補算出部141が算出した目標の航跡の候補を表わす航跡候補データを記憶する。
【0096】
例えば、候補算出部141は、次のようにして、複数の航跡の候補を算出する。
【0097】
候補算出部141は、観測値記憶部121が記憶した観測値データに基づいて、処理装置911を用いて、基準時刻tを決定する。例えば、候補算出部141は、観測値記憶部121が記憶した観測値データのなかから、観測時刻が最も新しい観測値データを抽出し、抽出した観測値データの観測時刻を基準時刻tとする。
【0098】
候補算出部141は、処理装置911を用いて、基準時刻tにおける目標の加速度ベクトルの推定値の候補として、一つの加速度ベクトルを設定する。目標の加速度ベクトルの推定値の候補を一つの加速度ベクトルに固定することにより、実際の目標の加速度ベクトルとの差がモデル化誤差となる。モデル化誤差の期待値は、設定した加速度ベクトルが、想定される実際の目標の加速度ベクトルの確率分布に基づく平均(期待値)に一致する場合に最小になる。想定される実際の目標の加速度ベクトルの確率分布について事前の情報がない場合、候補算出部141は、例えば、目標の加速度ベクトルの推定値の候補として、ゼロベクトルを設定する。あるいは、前回推定した航跡がある場合など、想定される実際の目標の加速度ベクトルの確率分布について事前の情報がある場合、候補算出部141は、目標の加速度ベクトルの推定値の候補として、例えば、事前の情報に基づいて加速度ベクトルの平均値(期待値)を算出する。例えば、候補算出部141は、前回推定した航跡を表わす状態量ベクトルから加速度ベクトルを取り出して、今回の目標の加速度ベクトルの推定値の候補とする。
【0099】
候補算出部141は、処理装置911を用いて、基準時刻tにおける目標の速度ベクトルの推定値の候補として、複数の速度ベクトルを設定する。例えば、候補算出部141は、想定される実際の目標の速度ベクトルの分布範囲のなかから、所定の数の速度ベクトルを、なるべく均等に分布するよう選択して、目標の速度ベクトルの推定値の候補とする。想定される実際の目標の速度ベクトルの分布範囲について事前の情報がない場合、候補算出部141は、目標の速度ベクトルの推定値の候補として、例えば、ゼロベクトルを中心とし、絶対値が所定の想定最大速度を超えない球状の範囲内から、所定の数の速度ベクトルを選択する。あるいは、前回推定した航跡がある場合など、想定される実際の目標の速度ベクトルの分布範囲について事前の情報がある場合、候補算出部141は、目標の速度ベクトルの推定値の候補として、例えば、事前の情報に基づいて目標の速度ベクトルの分布範囲を求め、その範囲内から、所定の数の速度ベクトルを選択する。
【0100】
候補算出部141は、処理装置911を用いて、基準時刻tにおける目標の位置ベクトルの推定値の候補として、複数の位置ベクトルを設定する。例えば、候補算出部141は、処理装置911を用いて、観測値記憶部121が記憶した観測値データのなかから、観測時刻と基準時刻tとの差が所定の閾値以内である観測値データを抽出する。閾値は、目標が存在しても観測されない欠測の可能性を考慮して、その時間の間に少なくとも1回は、目標が観測される時間に設定しておく。例えば、欠測の可能性が10%ある場合、センサ810が少なくとも3回は観測を行う時間を閾値とすれば、その間に目標が1回も観測されない可能性は、0.1%以下である。
候補算出部141は、処理装置911を用いて、抽出した観測値データに基づいて、観測時刻tにおける目標の位置ベクトルの推定値の候補を算出する。例えば、候補算出部141は、処理装置911を用いて、次の式の右辺を計算することにより、観測時刻tにおける目標の位置ベクトルの推定値の候補を算出する。
【数29】

ただし、p(t)は、候補算出部141が抽出した観測値データの観測時刻tにおける目標の位置ベクトルの推定値の候補である。p及びpは、候補算出部141が抽出した観測値データが表わす座標である。pzeは、目標の位置のZ座標(センサ810から目標までの距離)の推定値の候補である。例えば、センサ810の観測可能範囲に基づいて、想定される目標の位置のZ座標の範囲をあらかじめ定めておき、候補算出部141は、想定される範囲内から、複数のZ座標の推定値の候補を選択する。候補算出部141は、選択したZ座標の推定値の候補に基づいて、観測時刻tにおける目標の位置のX座標及びY座標の推定値の候補を算出する。候補算出部141は、目標の位置ベクトルの推定値の候補を、抽出した観測値データの数と選択したZ座標の推定値の候補の数との積に等しい数だけ算出する。
【0101】
候補算出部141は、処理装置911を用いて、算出した目標の速度ベクトルの推定値の候補に基づいて、目標の位置ベクトルの推定値の候補の時刻合わせを行い、時刻tにおける目標の位置ベクトルの推定値の候補を、基準時刻tにおける目標の位置ベクトルの推定値の候補へ変換する。例えば、候補算出部141は、処理装置911を用いて、次の式の右辺を計算することにより、基準時刻tにおける目標の位置ベクトルの推定値の候補を算出する。
【数30】

ただし、p(t)は、基準時刻tにおける目標の位置ベクトルの推定値の候補である。p(t)は、候補算出部141が抽出した観測値データの観測時刻tにおける目標の位置ベクトルの推定値の候補である。pドット(t)は、候補算出部141が算出した基準時刻tにおける目標の速度ベクトルの推定値の候補である。Δtは、観測時刻tから基準時刻tを差し引いた差である。
【0102】
候補算出部141は、処理装置911を用いて、算出した基準時刻tにおける目標の位置ベクトルの推定値の候補と、算出した基準時刻tにおける目標の速度ベクトルの候補と、算出した基準時刻tにおける目標の加速度ベクトルの候補とを組み合わせて、目標の航跡を表わす状態量ベクトルの推定値の候補を生成する。
【数31】

ただし、x(t)は、基準時刻tにおける目標の状態量ベクトルの推定値の候補である。pツードット(t)は、候補算出部141が算出した基準時刻tにおける目標の加速度ベクトルの推定値の候補である。
候補算出部141は、状態量ベクトルの推定値の候補を、算出した目標の位置ベクトルの推定値の候補の数と、算出した目標の速度ベクトルの推定値の候補の数との積に等しい数だけ算出する。
【0103】
候補算出部141は、処理装置911を用いて、算出した目標の状態量ベクトルの推定値の候補のなかに重複するものがあれば統合する。例えば、候補算出部141は、処理装置911を用いて、算出した状態量ベクトルの推定値の候補の間のユークリッド距離を算出し、算出したユークリッド距離が所定の閾値以内である場合に、重複と判断する。候補算出部141は、処理装置911を用いて、重複と判定した2つの位置ベクトルの推定値の候補を平均して、1つに統合する。3つ以上の候補が重複している可能性がある場合、候補算出部141は、例えば、距離が最も近い2つの候補を統合し、統合後の新しい候補と、他の候補との間の距離を算出して、重複しているか否かを判定する。
【0104】
なお、目標が観測された観測時刻tと基準時刻tとの間の経過時間Δtが短い場合、経過時間Δtの間における目標の移動を考慮せず、観測時刻tにおける目標の位置ベクトルの推定値の候補を、そのまま、基準時刻tにおける目標の位置ベクトルの推定値の候補とする構成であってもよい。
【0105】
誤差上限値記憶部152は、記憶装置914を用いて、候補算出部141が算出した目標の状態量ベクトルの推定値の候補において候補を複数算出しなかった状態量について、その状態量の推定値の候補を固定したことによるモデル化誤差の上限値を記憶する。例えば、候補算出部141が加速度ベクトルの推定値の候補を固定した場合、誤差上限値記憶部152は、実際の加速度ベクトルと、候補算出部141が固定した加速度ベクトルの推定値の候補との間の差をモデル化誤差とし、その上限値を記憶する。誤差上限値記憶部152が記憶するモデル化誤差の上限値は、誤差上限初期化部151や誤差上限値更新部174が算出する。例えば、誤差上限値記憶部152は、次の式で表わされる最大加速度誤差ベクトルを、モデル化誤差上限値ベクトルとして記憶する。
【数32】

ただし、Δpツードットmaxは、最大加速度誤差ベクトルである。Δamaxは、最大加速度誤差ベクトルのX方向成分である。Δamaxは、最大加速度誤差ベクトルのY方向成分である。Δamaxは、最大加速度誤差ベクトルのZ方向成分である。
【0106】
誤差上限初期化部151は、処理装置911を用いて、誤差上限値記憶部152が記憶するモデル化誤差の上限値を初期化する。例えば、想定される目標の加速度ベクトルの確率分布について事前の情報がなく、候補算出部141が目標の加速度ベクトルの推定値の候補としてゼロベクトルを設定した場合、誤差上限初期化部151は、目標の想定最大加速度を、モデル化誤差の上限値とする。あるいは、例えば、想定される目標の加速度ベクトルの確率分布について事前の情報があり、候補算出部141がその期待値を目標の加速度ベクトルの推定値の候補として設定した場合、誤差上限初期化部151は、想定される目標の加速度ベクトルの確率分布に基づいて、モデル化誤差の上限値を算出する。誤差上限値記憶部152は、記憶装置914を用いて、誤差上限初期化部151が算出したモデル化誤差の上限値を、モデル化誤差の上限値の初期値として記憶する。
【0107】
ゲート行列算出部153は、誤差上限値記憶部152が記憶したモデル化誤差の上限値に基づいて、ゲート行列を算出する。ゲート行列は、観測時刻ごとに異なるので、ゲート行列算出部153は、相関処理の対象となる観測値データの観測時刻ごとに、ゲート行列を算出する。例えば、ゲート行列算出部153は、処理装置911を用いて、観測値記憶部121が記憶した観測値データのなかから、相関処理の対象となる観測値データを抽出する。ゲート行列算出部153は、抽出した観測値データの観測時刻tに基づいて、ゲート行列を算出する。
例えば、ゲート行列算出部153は、処理装置911を用いて、次の式の右辺を計算することにより、モデル化誤差共分散行列を算出する。
【数33】

ただし、Q(t+Δt)は、基準時刻tから時間Δtが経過した時刻におけるモデル化誤差共分散行列である。diagは、n次元ベクトルのi番目の成分をi行i列成分とするn次の対角行列を表わす。
ゲート行列算出部153は、算出したモデル化誤差共分散行列に基づいて、処理装置911を用いて、次の式の右辺を計算することにより、時刻tにおけるゲート行列を算出する。
【0108】
S(t)=R(t)+Q(t)
ただし、S(t)は、ゲート行列算出部153が抽出した観測値データの観測時刻tにおけるゲート行列である。R(t)は、その観測時刻tにおける観測値データの観測誤差共分散行列である。Q(t)は、その観測時刻tについてゲート行列算出部153が算出したモデル化誤差共分散行列である。S(t)、R(t)、Q(t)は、いずれも、観測値ベクトルz(t)と同じ次数の正方行列である。
【0109】
相関部161は、候補記憶部142が記憶した目標の航跡の候補それぞれについて、処理装置911を用いて、ゲート行列算出部153が算出したゲート行列に基づいて、観測値記憶部121が記憶した観測値データのなかから、その航跡の目標を観測した観測値であると考えられる観測値を表わす観測値データを抽出する。
例えば、相関部161は、次のようにして、観測値データを抽出する。
【0110】
相関部161は、候補記憶部142が記憶した目標の航跡の候補それぞれについて、処理装置911を用いて、観測値記憶部121が記憶した観測値データのうち相関処理の対象となる観測値データの観測時刻tに基づいて、時刻tにおける目標の状態量ベクトルの推定値を算出する。相関部161は、例えば、次の式の右辺を計算することにより、時刻tにおける目標の状態量ベクトルの推定値を算出する。
【0111】
(t)=F(Δt)x(t
ただし、x(t)は、基準時刻tから時間Δtが経過した時刻t=t+Δtにおける目標の状態量ベクトルの推定値である。x(t)は、候補記憶部142が記憶した目標の航跡の候補である基準時刻tにおける目標の状態量ベクトルの推定値の候補である。
【0112】
相関部161は、算出した時刻tにおける目標の状態量ベクトルの推定値と、観測値データが表わす時刻tにおける観測値ベクトルとに基づいて、処理装置911を用いて、その観測値データの観測誤差の推定値を算出する。相関部161は、例えば、次の式の右辺を計算することにより、その観測値データの観測誤差の推定値を算出する。
【0113】
Δz=z−h(x(t))
ただし、Δzは、その観測値データの観測誤差の推定値である。zは、観測値データが表わす観測値ベクトルである。tは、その観測値データの観測時刻である。hは、観測関数である。
【0114】
相関部161は、算出した観測誤差の推定値に基づいて、処理装置911を用いて、その観測値データが、その航跡の目標を観測した観測値を表わす観測値データであるか否かを判定する。例えば、相関部161は、ゲート行列算出部153が算出したゲート行列に基づいて、次の式が表わす条件が満たされる場合に、その観測値データが、その航跡の目標を観測した観測値を表わす観測値データであると判断する。
【0115】
ΔzS(t)−1Δz≦d
ただし、上付きのTは、転置行列を表わす。上付きの−1は、逆行列を表わす。dは、所定の閾値である。
【0116】
尤度算出部162は、候補記憶部142が記憶した目標の軌跡の候補それぞれについて、処理装置911を用いて、ゲート行列算出部153が算出したゲート行列と、相関部161が抽出した観測値データとに基づいて、その目標の軌跡の候補の尤度を算出する。例えば、尤度算出部162は、次の式の右辺を計算することにより、その軌跡の候補の尤度を算出する。
【数34】

ただし、φは、その軌跡の候補の尤度である。Tは、相関処理の対象となる観測値データの観測時刻を要素とする集合である。logは、ネイピア数を底とする対数(自然対数)関数を表わす。Pは、探知確率パラメータである。探知確率パラメータPは、0より大きく1より小さい実数であり、センサ810の観測範囲内に目標が存在するとき、その目標を観測した観測値を観測値取得部111が抽出する確率を表わす。βは、不要信号密度パラメータである。不要信号密度パラメータβは、0より大きい実数であり、ノイズなどの影響による不要信号が観測値空間の単位体積あたりに存在する平均個数を表わす。探知確率パラメータPや不要信号密度パラメータβは、あらかじめ設定された値を尤度算出部162が記憶装置914を用いて記憶しておく。なお、尤度算出部162は、追尾結果などに基づいて、探知確率パラメータPや不要信号密度パラメータβを更新する構成であってもよい。
detは、行列式を表わす。πは、円周率である。S(t)は、ゲート行列算出部153が算出した時刻tにおけるゲート行列である。G(t)は、その航跡の候補について相関部161が抽出した観測値データのうち、観測時刻が時刻tである観測値データが表わす観測値ベクトルを要素とする集合である。expは、ネイピア数を底とする指数関数を表わす。Δzは、観測値ベクトルzについて相関部161が算出した観測誤差の推定値である。上付きのTは、転置行列を表わす。上付きの−1は、逆行列を表わす。
【0117】
選択部163は、尤度算出部162が算出した尤度に基づいて、処理装置911を用いて、候補記憶部142が記憶した航跡の候補のなかから、正しい航跡に近い航跡である可能性の高い候補を抽出する。選択部163は、航跡の候補を1つだけ抽出する構成であってもよいし、複数の候補を抽出する構成であってもよい。
例えば、選択部163は、処理装置911を用いて、候補記憶部142が記憶した航跡の候補のなかから、尤度算出部162が算出した尤度が大きい順に、所定の数の候補を抽出する。あるいは、例えば、選択部163は、処理装置911を用いて、尤度算出部162が算出した尤度のなかから、一番大きな尤度を抽出し、候補記憶部142が記憶した航跡の候補それぞれについて、尤度算出部162が算出した尤度を、抽出した最大尤度で割った商を算出する。選択部163は、処理装置911を用いて、算出した商が所定の閾値より大きい候補を抽出する。
【0118】
推定値記憶部171は、記憶装置914を用いて、選択部163が抽出した航跡の候補を表わす状態量ベクトルの推定値を記憶する。
【0119】
推定値更新部172は、処理装置911を用いて、推定値記憶部171が記憶した状態量ベクトルの推定値を出発点として、もっと尤度が大きい航跡を表わす状態量ベクトルの推定値を算出する。推定値更新部172は、算出した状態量ベクトルの推定値で、推定値記憶部171が記憶した状態量ベクトルの推定値を更新する。推定値記憶部171は、記憶装置914を用いて、更新された状態量ベクトルの推定値を記憶する。
【0120】
例えば、推定値更新部172は、推定値記憶部171が記憶した状態量ベクトルの推定値をわずかに変化させる。推定値更新部172が変化させた状態量ベクトルの推定値に基づいて、相関部161が相関処理をし、尤度算出部162が尤度を算出する。推定値更新部172は、尤度算出部162が算出した尤度に基づいて、推定値記憶部171が記憶した状態量ベクトルの推定値の近傍における尤度の変化傾向を把握する。推定値更新部172は、把握した変化傾向に基づいて、状態量ベクトルの推定値をどの方向にどのくらい変化させたら、もっと尤度が大きい航跡を表わす状態量ベクトルの推定値となるかを算出する。
具体的には、例えば次のようにして、もっと尤度が大きい航跡を表わす状態量ベクトルの推定値を算出する。
【0121】
推定値更新部172は、処理装置911を用いて、推定値記憶部171が記憶した状態量ベクトルの推定値それぞれについて、状態量ベクトルの各成分を微かに変化させた状態量ベクトルを算出する。推定値更新部172は、例えば、状態量ベクトルのある成分に、所定の微小増分を加えた状態量ベクトルと、その成分から微小増分を差し引いた状態量ベクトルとを算出する。例えば、状態量ベクトルが九次元ベクトルである場合、推定値更新部172は、各成分について2つの状態量ベクトルを算出するので、全部で18個の状態量ベクトルを算出する。
相関部161(詳細相関部)は、処理装置911を用いて、推定値更新部172が算出した18個の状態量ベクトルと、その元となった1つの状態量ベクトルとを合わせた合計19個の状態量ベクトルそれぞれについて、その状態量ベクトルが表わす航跡を、候補記憶部142が記憶した航跡の候補の代わりとして、観測値記憶部121が記憶した観測値データのなかから、観測値データを抽出する。尤度算出部162(詳細尤度算出部)は、処理装置911を用いて、相関部161と同様、19個の状態量ベクトルそれぞれについて、その状態量ベクトルが表わす航跡を、候補記憶部142が記憶した航跡の候補の代わりとして、その航跡の尤度を算出する。
推定値更新部172は、処理装置911を用いて、尤度算出部162が算出した尤度に基づいて、状態量ベクトルの各成分について、尤度が最大になると予想される値を算出する。
例えば、推定値更新部172は、状態量ベクトルの一つの成分に注目し、注目した成分以外の成分の値は変えずに、その成分の値だけを変えて算出した尤度を比較する。推定値更新部172は、その成分のもとの値uを使った状態量ベクトルについて算出した尤度φ(u)と、微小増分Δuを加えた値を使った状態量ベクトルについて算出した尤度φ(u+Δu)と、微小増分Δuを差し引いた値を使った状態量ベクトルについて算出した尤度φ(u−Δu)との間の大小関係を判定する。
φ(u−Δu)≧φ(u)≦φ(u+Δu)である場合、推定値更新部172は、もとの値uを、その成分の更新した値とする。
φ(u−Δu)<φ(u)>φ(u+Δu)である場合も、推定値更新部172は、もとの値uを、その成分の更新した値とする。
φ(u−Δu)<φ(u)≦φ(u+Δu)である場合、推定値更新部172は、尤度φ(u)から尤度φ(u−Δu)を差し引いた差Δφ=φ(u)−φ(u−Δu)と、尤度φ(u+Δu)から尤度φ(u)を差し引いた差Δφ=φ(u+Δu)−φ(u)とを算出する。推定値更新部172は、算出した2つの差の大小関係を比較する。
Δφ≧Δφである場合、推定値更新部172は、0超所定の値未満の乱数rを生成し、生成した乱数rをもとの値uに加えた和u+rを、その成分の更新した値とする。
Δφ<Δφである場合、推定値更新部172は、差Δφと差Δφとの平均値Δφ=(Δφ+Δφ)/2と、差Δφから差Δφを差し引いた差Δφ=Δφ−Δφとを算出する。推定値更新部172は、算出した平均値Δφを、算出した差Δφで割った商を微小増分Δuに乗じた積をもとの値uに加えた和u+Δφ/Δφ・Δuを、その成分の更新した値とする。
φ(u−Δu)≧φ(u)>φ(u+Δu)である場合、推定値更新部172は、尤度φ(u)から尤度φ(u+Δu)を差し引いた差Δφ=φ(u)−φ(u+Δu)と、尤度φ(u−Δu)から尤度φ(u)を差し引いた差Δφ=φ(u−Δu)−φ(u)とを算出する。推定値更新部172は、算出した2つの差の大小関係を比較する。
Δφ≦Δφである場合、推定値更新部172は、0超所定の値未満の乱数rを生成し、生成した乱数rをもとの値uから差し引いた差u−rを、その成分の更新した値とする。
Δφ>Δφである場合、推定値更新部172は、差Δφと差Δφとの平均値Δφ=(Δφ+Δφ)/2と、差Δφから差Δφを差し引いた差Δφ=Δφ−Δφとを算出する。推定値更新部172は、算出した平均値Δφを、算出した差Δφで割った商を微小増分Δuに乗じた積をもとの値uから差し引いた差u−Δφ/Δφ・Δuを、その成分の更新した値とする。
推定値更新部172は、例えばこのようにして状態量ベクトルのすべての成分について更新した値を算出することにより、更新した状態量ベクトルの推定値を算出する。
推定値記憶部171は、記憶装置914を用いて、推定値更新部172が更新した状態量ベクトルの推定値を記憶する。
【0122】
収束判定部173は、処理装置911を用いて、状態量ベクトルの推定値が収束したか否かを判定する。例えば、収束判定部173は、推定値更新部172が更新した更新後の状態量ベクトルの推定値と、推定値更新部172が更新する前に推定値記憶部171が記憶していた更新前の状態量ベクトルの推定値との差ベクトルを算出する。収束判定部173は、算出した差ベクトルの絶対値を所定の閾値と比較する。差ベクトルの絶対値が閾値より小さい場合、収束判定部173は、状態量ベクトルの推定値が収束したと判定する。
【0123】
誤差上限値更新部174は、状態量ベクトルの推定値が収束したと収束判定部173が判定した場合、処理装置911を用いて、誤差上限値記憶部152が記憶したモデル化誤差の上限値を更新する。例えば、誤差上限値更新部174は、誤差上限値記憶部152が記憶したモデル化誤差の上限値に、0超1未満の所定の係数を乗じた積を算出して、更新したモデル化誤差の上限値とする。
誤差上限値記憶部152は、記憶装置914を用いて、誤差上限値更新部174が更新したモデル化誤差の上限値を記憶する。
【0124】
終了判定部175は、状態量ベクトルの推定値が収束したと収束判定部173が判定した場合、処理装置911を用いて、状態量ベクトルの推定を終了するか否かを判定する。例えば、終了判定部175は、状態量ベクトルの推定値が収束したと収束判定部173が判定したとき推定値記憶部171が記憶している状態量ベクトルの推定値を記憶する。状態量ベクトルの推定値が収束したと収束判定部173が判定したのが一回目である場合、終了判定部175は、状態量ベクトルの推定を終了しないと判定する。状態量ベクトルの推定値が収束したと収束判定部173が判定したのが二回目もしくはそれ以降である場合、終了判定部175は、記憶している前回の状態量ベクトルの推定値と、そのとき推定値記憶部171が記憶している状態量ベクトルの推定値との差ベクトルを算出する。終了判定部175は、算出した差ベクトルの絶対値を所定の閾値と比較する。差ベクトルの絶対値が閾値より小さい場合、終了判定部175は、状態量ベクトルの推定を終了すると判定する。
【0125】
推定値出力部176は、状態量ベクトルの推定を終了すると終了判定部175が判定した場合、出力装置913を用いて、推定値記憶部171が記憶した状態量ベクトルの推定値を出力する。なお、推定値記憶部171が記憶した状態量ベクトルの推定値のなかに、ほぼ同一の推定値が複数ある場合、推定値出力部176は、そのなかで最も尤度が大きい推定値だけを出力する構成であってもよい。
【0126】
図12は、この実施の形態における追尾処理S600の流れの一例を示すフローチャート図である。
追尾処理S600(推定方法、航跡推定方法)は、例えば、観測値取得工程S661と、候補算出工程S662と、誤差上限初期化工程S663と、2つのゲート行列算出工程S664,S683と、候補選択工程S665と、2つの相関工程S666,S673と、2つの尤度算出工程S667,S674と、選択工程S668と、推定値選択工程S671と、参照生成工程S672と、推定値更新工程S675と、収束判定工程S676と、個別終了判定工程S677と、全体終了判定工程S681と、誤差上限値更新工程S682と、推定値出力工程S684とを有する。
追尾装置100は、観測値取得工程S661から追尾処理S600を開始する。
【0127】
観測値取得工程S661において、観測値取得部111は、入力装置912を用いて、センサ810が出力した信号を入力する。観測値取得部111は、処理装置911を用いて、入力した信号に基づいて、観測値を取得する。観測値記憶部121は、記憶装置914を用いて、観測値取得部111が取得した観測値を記憶する。
【0128】
候補算出工程S662において、候補算出部141は、処理装置911を用いて、基準時刻tを決定する。候補算出部141は、処理装置911を用いて、決定した基準時刻tにおける目標の状態量の推定値の候補を複数生成する。候補記憶部142は、記憶装置914を用いて、候補算出部141が生成した複数の候補を記憶する。
【0129】
誤差上限初期化工程S663において、誤差上限初期化部151は、処理装置911を用いて、モデル化誤差の上限値の初期値を設定する。誤差上限値記憶部152は、記憶装置914を用いて、誤差上限初期化部151が設定した初期値を、モデル化誤差の上限値として記憶する。
【0130】
ゲート行列算出工程S664において、ゲート行列算出部153は、処理装置911を用いて、誤差上限初期化工程S663で誤差上限値記憶部152が記憶したモデル化誤差の上限値と、候補算出工程S662で候補算出部141が決定した基準時刻tとに基づいて、複数の観測時刻tそれぞれにおけるゲート行列S(t)を算出する。
【0131】
候補選択工程S665において、選択部163は、処理装置911を用いて、候補記憶部142が記憶した複数の候補のなかから、まだ尤度を算出していない候補を一つ選択する。
候補記憶部142が記憶したすべての候補について尤度を算出済であり、選択すべき候補がない場合、選択部163は、処理装置911を用いて、選択工程S668へ処理を進める。
候補記憶部142が記憶した候補のなかに、まだ尤度を算出していない候補がある場合、選択部163は、処理装置911を用いて、まだ尤度を算出していない候補のなかから、候補を一つ選択する。選択部163は、処理装置911を用いて、相関工程S666へ処理を進める。
【0132】
相関工程S666において、相関部161は、処理装置911を用いて、ゲート行列算出工程S664でゲート行列算出部153が算出したゲート行列S(t)と、候補選択工程S665で選択部163が選択した候補とに基づいて、観測値記憶部121が記憶した観測値のなかから、その候補が表わす航跡の目標を観測したと考えられる観測値を抽出する。
【0133】
尤度算出工程S667において、尤度算出部162は、処理装置911を用いて、ゲート行列算出工程S664でゲート行列算出部153が算出したゲート行列S(t)と、候補選択工程S665で選択部163が選択した候補と、相関工程S666で相関部161が抽出した観測値とに基づいて、その候補の尤度を算出する。選択部163は、記憶装置914を用いて、尤度算出部162が算出した尤度を記憶する。選択部163は、処理装置911を用いて、候補選択工程S665に処理を戻し、次の候補を選択する。
【0134】
選択工程S668において、選択部163は、処理装置911を用いて、尤度算出工程S667で記憶した尤度に基づいて、候補記憶部142が記憶した候補のなかから、正しい航跡に近い航跡を表わす候補を抽出する。推定値記憶部171は、記憶装置914を用いて、選択部163が抽出した候補(基準時刻tにおける状態量の推定値)を記憶する。
【0135】
推定値選択工程S671において、推定値更新部172は、処理装置911を用いて、選択工程S668で推定値記憶部171が記憶した推定値のうち、まだ推定が終了していない推定値のなかから、まだ更新をしていない推定値を一つ選択する。
推定値記憶部171が記憶した推定値すべてについての更新が終わり、選択すべき推定値がない場合、推定値更新部172は、処理装置911を用いて、全体終了判定工程S681へ処理を進める。
推定値記憶部171が記憶した推定値のなかに、まだ更新をしていない推定値がある場合、推定値更新部172は、処理装置911を用いて、まだ更新をしていない推定値のなかから、推定値を一つ選択する。推定値更新部172は、処理装置911を用いて、参照生成工程S672へ処理を進める。
【0136】
参照生成工程S672において、推定値更新部172は、処理装置911を用いて、所定の複数の方向のなかから、推定値を変化させる方向を選択する。
すべての方向についての処理が終わり、選択すべき方向がない場合、推定値更新部172は、処理装置911を用いて、推定値更新工程S675へ処理を進める。
所定の複数の方向のなかに、まだ処理していない方向がある場合、推定値更新部172は、処理装置911を用いて、まだ処理していない方向のなかから、方向を一つ選択する。推定値更新部172は、処理装置911を用いて、選択した方向へ向けて、推定値選択工程S671で選択した推定値をわずかに変化させる。推定値更新部172は、処理装置911を用いて、相関工程S673へ処理を進める。
【0137】
相関工程S673において、相関部161(詳細相関部)は、処理装置911を用いて、ゲート行列算出工程S664またはS683でゲート行列算出部153が算出したゲート行列S(t)と、参照生成工程S672で推定値更新部172が変化させた推定値とに基づいて、観測値記憶部121が記憶した観測値のなかから、その推定値が表わす航跡の目標を観測したと考えられる観測値を抽出する。
【0138】
尤度算出工程S674において、尤度算出部162(詳細尤度算出部)は、処理装置911を用いて、ゲート行列算出工程S664またはS683でゲート行列算出部153が算出したゲート行列S(t)と、推定値更新部172で推定値更新部172が変化させた推定値と、相関工程S673で相関部161が抽出した観測値とに基づいて、その推定値の尤度を算出する。推定値更新部172は、記憶装置914を用いて、尤度算出部162が算出した尤度を記憶する。推定値更新部172は、処理装置911を用いて、参照生成工程S672に処理を戻し、次の方向を選択する。
【0139】
推定値更新工程S675において、推定値更新部172は、処理装置911を用いて、尤度算出工程S674で記憶した尤度に基づいて、更新前よりも尤度が大きくなる推定値を算出する。推定値記憶部171は、処理装置911を用いて、推定値選択工程S671で推定値更新部172が選択した推定値を、推定値更新部172が算出した推定値で更新し、記憶装置914を用いて、更新した推定値を記憶する。
【0140】
収束判定工程S676において、収束判定部173は、処理装置911を用いて、推定値選択工程S671で推定値更新部172が選択した推定値が収束したか否かを判定する。
収束したと判定した場合、収束判定部173は、処理装置911を用いて、個別終了判定工程S677へ処理を進める。
まだ収束していないと判定した場合、推定値更新部172は、処理装置911を用いて、参照生成工程S672に処理を戻し、更に尤度が大きくなる推定値を算出する。
【0141】
個別終了判定工程S677において、終了判定部175は、処理装置911を用いて、推定値選択工程S671で推定値更新部172が選択した推定値について、状態量の推定を終了するか否かを判定する。終了判定部175は、記憶装置914を用いて、判定した判定結果を記憶する。
推定値更新部172は、処理装置911を用いて、推定値選択工程S671に処理を戻し、次の推定値を選択する。
【0142】
全体終了判定工程S681において、終了判定部175は、処理装置911を用いて、個別終了判定工程S677で記憶した判定結果に基づいて、推定値記憶部171が記憶したすべての推定値について、状態量の推定を終了すると判定したか否かを判定する。
すべての推定値について状態量の推定を終了すると判定し、状態量の推定を終了しないと判定した推定値がない場合、終了判定部175は、処理装置911を用いて、推定値出力工程S684へ処理を進める。
状態量の推定を終了しないと判定した推定値がまだある場合、終了判定部175は、処理装置911を用いて、誤差上限値更新工程S682へ処理を進める。
【0143】
誤差上限値更新工程S682において、誤差上限値更新部174は、処理装置911を用いて、誤差上限値記憶部152が記憶したモデル化誤差の上限値を更新する。誤差上限値記憶部152は、記憶装置914を用いて、誤差上限値更新部174が更新したモデル化誤差の上限値を記憶する。
【0144】
ゲート行列算出工程S683において、ゲート行列算出部153は、処理装置911を用いて、誤差上限値更新工程S682で誤差上限値記憶部152が記憶した更新後のモデル化誤差の上限値に基づいて、観測時刻tそれぞれにおけるゲート行列S(t)を算出する。
推定値更新部172は、推定値選択工程S671に処理を戻す。推定値選択工程S671において、推定値更新部172は、推定値記憶部171が記憶した推定値のうち、個別終了判定工程S677で状態量の推定を終了しないと終了判定部175が判定した推定値について、また最初から推定値を一つずつ選択する。
【0145】
推定値出力工程S684において、推定値出力部176は、出力装置913を用いて、推定値記憶部171が記憶した推定値を出力する。
【0146】
以上、各実施の形態で説明した具体的な構成は一例であり、これに限定されるものではない。例えば、重要でない部分の構成を他の構成で置き換えた構成としてもよし、異なる実施の形態で説明した構成を組み合わせた構成としてもよい。
【0147】
推定装置(追尾装置100)は、データを処理する処理装置(911)と、候補算出部(141;航跡候補生成部140)と、相関部(161;概探索部160)と、尤度算出部(162;概探索部160)と、選択部(163;概探索部160)と、推定値更新部(172;詳細探索部170)と、推定値出力部(176;詳細探索部170)とを有する。
上記候補算出部は、目標の複数の状態量(例えば位置・速度・加速度)のうち所定の第一の状態量(例えば加速度)を除く他の状態量(例えば位置・速度)について、上記処理装置を用いて、上記他の状態量の推定値の複数の候補値を算出する。
上記相関部は、上記候補算出部が算出した上記他の状態量の推定値の複数の候補値それぞれについて、上記処理装置を用いて、上記第一の状態量が所定の推定値であり、上記第一の状態量が上記推定値であると仮定することによって生じるモデル化誤差が所定のモデル化誤差上限値以下であり、上記他の状態量が上記候補値であるとの仮定に基づいて、上記目標を観測した観測値を含む複数の観測値のなかから、上記目標を観測した観測値を抽出する。
上記尤度算出部は、上記候補算出部が算出した上記他の状態量の推定値の複数の候補値それぞれについて、上記処理装置を用いて、上記相関部が仮定したモデル化誤差上限値と、上記相関部が抽出した観測値とに基づいて、上記他の状態量が上記候補値である尤度を算出する。
上記候補選択部は、上記処理装置を用いて、上記尤度算出部が算出した尤度に基づいて、上記候補算出部が算出した複数の候補値のなかから、上記他の状態量の推定値を選択する。
上記推定値更新部は、上記処理装置を用いて、上記相関部が仮定した第一の状態量の推定値と上記候補選択部が選択した他の状態量の推定値とを出発点として、上記目標の複数の状態量が推定値である尤度が高くなるよう、上記複数の状態量の推定値を更新する。
上記推定値出力部は、上記処理装置を用いて、上記推定値更新部が更新した上記複数の状態量の推定値を出力する。
【0148】
第一の状態量が所定の推定値であると仮定して、状態量の推定値の候補値を算出するので、状態量の推定値の組み合わせの数が減り、相関部や尤度算出部の処理を減らすことができる。これにより、目標の航跡の推定に必要な計算量が減るので、処理能力の低い処理装置を用いてもリアルタイムでの処理が可能になるなど、追尾装置100のコストを下げることができる。
【0149】
上記推定装置(100)は、収束判定部(173;詳細探索部170)と、誤差上限値更新部(174;詳細探索部170)と、詳細相関部(相関部161;詳細探索部170)と、詳細尤度算出部(尤度算出部162;詳細探索部170)とを有する。
上記収束判定部は、上記処理装置を用いて、上記推定値更新部が更新した上記複数の状態量の推定値が収束したか否かを判定する。
上記誤差上限値更新部は、上記推定値更新部が上記複数の状態量の推定量を更新した場合に、上記処理装置を用いて、上記モデル化誤差上限値を小さくする。
上記詳細相関部は、上記処理装置を用いて、上記複数の状態量が、上記推定値更新部が更新した推定値であり、上記モデル化誤差が、上記誤差上限値更新部が小さくしたモデル化誤差上限値以下であるとの仮定に基づいて、上記複数の観測値のなかから上記目標を観測した観測値を抽出する。
上記詳細尤度算出部は、上記推定値更新部が更新した上記複数の状態量の推定値について、上記処理装置を用いて、上記誤差上限値更新部が小さくしたモデル化誤差上限値と、上記詳細相関部が抽出した観測値とに基づいて、上記複数の状態量が上記推定値である尤度を算出する。
上記推定値更新部は、上記推定値更新部が更新した上記複数の状態量の推定値が収束していないと上記収束判定部が判定した場合に、上記処理装置を用いて、更新した上記複数の状態量の推定値を出発点として、上記詳細尤度算出部が算出する尤度が更に高くなるよう、上記複数の状態量の推定値を更新する。
上記推定値出力部は、上記推定値更新部が更新した上記複数の状態量の推定値が収束したと上記収束判定部が判定した場合に、上記処理装置を用いて、上記推定値更新部が更新した上記複数の状態量の推定値を出力する。
【0150】
状態量の推定値が収束するまで更新処理を繰り返すので、航跡の推定精度を高くすることができる。
【0151】
上記推定装置(100)は、収束判定部(173;詳細探索部170)と、終了判定部(175;詳細探索部170)と、誤差上限値更新部(174;詳細探索部170)と、詳細相関部(相関部161;詳細探索部170)と、詳細尤度算出部(尤度算出部162;詳細探索部170)とを有する。
上記収束判定部は、上記処理装置を用いて、上記推定値更新部が更新した上記複数の状態量の推定値が収束したか否かを判定する。
上記終了判定部は、上記推定値更新部が更新した上記複数の状態量の推定値が収束したと上記収束判定部が判定した場合に、上記処理装置を用いて、上記複数の状態量の推定を終了するか否かを判定する。
上記誤差上限値更新部は、上記複数の状態量の推定を終了しないと上記終了判定部が判定した場合に、上記処理装置を用いて、上記モデル化誤差上限値を小さくする。
上記詳細相関部は、上記処理装置を用いて、上記推定値更新部が更新した上記複数の状態量の推定値と、上記相関部が仮定し若しくは上記誤差上限値更新部が小さくしたモデル化誤差上限値とに基づいて、上記複数の観測値のなかから上記目標を観測した観測値を抽出する。
上記詳細尤度算出部は、上記推定値更新部が更新した上記複数の状態量の推定値について、上記処理装置を用いて、上記相関部が仮定し若しくは上記誤差上限値更新部が小さくしたモデル化誤差上限値と、上記詳細相関部が抽出した観測値とに基づいて、上記複数の状態量が上記推定値である尤度を算出する。
上記推定値更新部は、上記推定値更新部が更新した上記複数の状態量の推定値が収束していないと上記収束判定部が判定した場合、および、上記複数の状態量の推定を終了しないと上記終了判定部が判定した場合に、上記処理装置を用いて、更新した上記複数の状態量の推定値を出発点として、上記詳細尤度算出部が算出する尤度が更に高くなるよう、上記複数の状態量の推定値を更新する。
上記推定値出力部は、上記複数の状態量の推定を終了すると上記終了判定部が判定した場合に、上記処理装置を用いて、上記推定値更新部が更新した上記複数の状態量の推定値を出力する。
【0152】
状態量の推定値が収束するまで、モデル化誤差上限値を変えずに更新処理を繰り返し、状態量の推定値が収束した場合に、モデル化誤差上限値を小さくして、更に、状態量の推定値が収束するまで更新処理を繰り返すので、状態量の推定精度を高くすることができる。
【0153】
上記詳細相関部は、上記処理装置を用いて、上記推定値更新部が更新した上記複数の状態量の推定値と、上記相関部が仮定し若しくは上記誤差上限値更新部が小さくしたモデル化誤差上限値とに基づいて、上記複数の観測値が上記目標を観測した観測値であるかと判定する範囲を算出するために使用するゲート行列を算出する。
上記終了判定部は、上記処理装置を用いて、上記詳細相関部が算出したゲート行列に基づいて、上記複数の状態量の推定を終了するか否かを判定する。
【0154】
ゲート行列に基づいて状態量の推定を終了するか否かを判定するので、精度よく終了判定をすることができる。
【0155】
上記終了判定部は、上記複数の状態量の推定を終了するか否かを判定した回数が所定の回数を超えた場合に、上記処理装置を用いて、上記複数の状態量の推定を終了すると判定する。
【0156】
判定回数に基づいて状態量の推定を終了するか否かを判定するので、終了判定に必要な計算量を少なくすることができる。
【0157】
上記誤差上限値更新部は、上記処理装置を用いて、次の式の右辺を計算することにより、上記モデル化誤差上限値を更新する。
max←√[(amax−Δa
ただし、amaxは、上記モデル化誤差上限値を表わすベクトルのi番目の成分。Δaは、上記推定値更新部が更新した後の第一の状態量の推定値を表わすベクトルのi番目の成分と、上記推定値更新部が更新する前の第一の状態量の推定値を表わすベクトルのi番目の成分との差。
【0158】
これにより、モデル化誤差上限値を最適に更新することができるので、全体としての繰り返し回数が減り、目標の航跡の推定に必要な計算量を少なくすることができる。
【0159】
上記誤差上限値更新部は、上記処理装置を用いて、次の式の右辺を計算することにより、上記モデル化誤差上限値を更新する。
max←λ・amax
ただし、λは、モデル化誤差上限値の逓減率であり、0超1未満の実数。
【0160】
これにより、モデル化誤差上限値の更新に必要な計算量を少なくすることができるので、目標の航跡の推定に必要な計算量を少なくすることができる。
【0161】
上記相関部及び上記詳細相関部は、上記処理装置を用いて、次の式の右辺を計算することにより、ゲート行列を算出し、算出したゲート行列を使用して、上記複数の観測値が上記目標を観測した観測値であると判定する範囲を算出し、上記複数の観測値のうち、算出した範囲内に入る観測値を、上記目標を観測した観測値であると判定する。
S←R+Q
ただし、Sは、ゲート行列であり、上記観測値を表わすベクトルの成分の数を次数とする正方行列。Rは、観測雑音共分散行列であり、上記観測値を表わすベクトルの成分の数を次数とし、上記観測値を表わすベクトルの各成分の観測精度の二乗を対角成分とする対角行列。Qは、モデル化誤差共分散行列であり、上記観測値を表わすベクトルの成分の数を次数とし、上記モデル化誤差上限値を表わすベクトルの各成分から算出される値を対角成分とする対角行列。
【0162】
観測誤差とモデル化誤差とを考慮したゲート行列を使って、目標を観測した観測値を判定するので、モデル化誤差が大きい場合でも、目標を観測した観測値を精度よく判定することができる。
【0163】
上記相関部及び上記詳細相関部は、上記処理装置を用いて、次の式の右辺を計算することにより、ゲート行列を算出し、算出したゲート行列を使用して、上記複数の観測値が上記目標を観測した観測値であると判定する範囲を算出し、上記複数の観測値のうち、算出した範囲内に入る観測値を、上記目標を観測した観測値であると判定する。
S←R+Q+Q’
ただし、Q’は、量子化誤差共分散行列であり、上記観測値を表わすベクトルの成分の数を次数とする対称行列。
【0164】
観測誤差とモデル化誤差と量子化誤差とを考慮したゲート行列を使って、目標を観測した観測値を判定するので、量子化誤差が大きい場合でも、目標を観測した観測値を精度よく判定することができる。
【0165】
上記推定装置は、複数の時刻における上記目標の位置を、上記目標を観測した観測値として入力し、所定の基準時刻における上記目標の位置及び速度及び加速度を、上記複数の状態量として推定する。
上記相関部は、上記処理装置を用いて、上記目標の加速度を上記第一の状態量として、上記加速度の推定値及びモデル化誤差上限値を仮定する。
上記相関部及び上記詳細相関部は、上記処理装置を用いて、次の式の右辺を計算することにより、上記モデル化誤差共分散行列の対角成分を算出する。
(t)←(amax・Δt/2)
ただし、q(t)は、時刻tにおけるモデル化誤差共分散行列Qのi行i列成分。amaxは、上記加速度のモデル化誤差上限値を表わすベクトルのi番目の成分。Δtは、上記基準時刻から上記時刻tまでの経過時間。
【0166】
目標の加速度のモデル化誤差上限値に基づいて目標の位置のモデル化誤差共分散行列を算出し、算出した目標の位置のモデル化誤差共分散行列に基づいて、ゲート行列を算出するので、目標を観測した観測値を精度よく判定することができる。
【0167】
上記推定装置は、複数の時刻における観測値を入力し、所定の基準時刻における複数の状態量を推定する。
上記相関部及び上記詳細相関部は、上記処理装置を用いて、上記基準時刻における上記目標の複数の状態量の推定値または推定値の候補に基づいて、上記複数の時刻それぞれにおける上記目標の複数の状態量を予測し、上記複数の観測値のうち次の不等式を満たす観測値を、上記目標を観測した観測値であると判定する。
[z−h(x)]・S−1・[z−h(x)]≦d
ただし、zは、上記観測値を表わす縦ベクトル。hは、観測関数を表わす。xは、時刻tにおける上記目標の複数の状態量の予測値を表わす縦ベクトル。上付きのTは、行列の転置を表わす。Sは、上記時刻tにおけるゲート行列。上付きの−1は、逆行列を表わす。dは、ゲートの大きさを表わすパラメータであり、0超の実数。
【0168】
観測時刻によって異なるゲート行列に基づいて目標を観測した観測値を判定するので、目標を観測した観測値を精度よく判定することができる。
【0169】
上記推定装置は、複数の時刻における観測値を入力し、所定の基準時刻における複数の状態量を推定する。
上記尤度算出部及び詳細尤度算出部は、上記処理装置を用いて、次の式の右辺を計算することにより、尤度を算出する。
【数35】

ただし、φ(Z|x)は、上記複数の状態量が上記推定値である尤度。lnは、自然対数を表わす。Pは、目標の探知確率を表わすパラメータであり、0超1未満の実数。βは、不要信号密度を表わすパラメータであり、0超の実数。detは、行列式を表わす。πは、円周率。Sは、時刻tにおけるゲート行列。mは、上記複数の観測値のうち、上記時刻tにおいて上記目標を観測した観測値であると、上記相関部または上記詳細相関部が判定した観測値の数。zk,jは、上記時刻tにおいて上記目標を観測した観測値であると、上記相関部または上記詳細相関部が判定した観測値を表わす縦ベクトル。hは、観測関数を表わす。xは、時刻tにおける上記目標の複数の状態量の予測値を表わす縦ベクトル。上付きのTは、行列の転置を表わす。上付きの−1は、逆行列を表わす。
【0170】
上の式を使って算出した尤度に基づいて、目標の航跡の候補を選択し、更新するので、航跡の推定精度を高くすることができる。
【0171】
上記収束判定部は、上記処理装置を用いて、上記推定値更新部が上記複数の状態量の推定値を更新したことによる上記尤度の変化に基づいて、上記複数の状態量の推定値が収束したか否かを判定する。
【0172】
状態量の推定値の更新による尤度の変化に基づいて収束判定をするので、収束判定を効率よくすることができ、目標の航跡の推定に必要な計算量を減らすことができる。
【0173】
上記収束判定部は、上記複数の状態量の推定値が収束したか否かを判定した回数が所定の回数を超えた場合に、上記処理装置を用いて、上記推定値更新部が上記複数の状態量の推定値が収束したと判定する。
【0174】
収束判定をした回数に基づいて収束判定をするので、収束判定に必要な計算量を減らすことができ、全体として、目標の航跡の推定に必要な計算量を減らすことができる。
【0175】
上記推定値更新部は、上記処理装置を用いて、上記尤度が高くなる上記複数の状態量の推定値の差分を算出し、上記複数の状態量それぞれについて、算出した差分を上記状態量の推定値に加えることにより、上記複数の状態量の推定値を更新する。
【0176】
上記推定値更新部は、上記処理装置を用いて、最急降下法およびニュートン法および準ニュートン法および共役分配法のいずれかにより、上記尤度が高くなる上記複数の状態量の推定値の差分を算出する。
【0177】
上記のような最適化方式を使って、状態量の推定値を更新するので、効率よく状態量の推定値を更新することができる。
【0178】
以上説明した推定装置(追尾装置100)は、データを処理する処理装置を有するコンピュータを上記推定装置として機能させるコンピュータプログラムを、コンピュータが実行することにより実現することができる。
【符号の説明】
【0179】
100 追尾装置、110 観測データ抽出部、111 観測値取得部、120 観測データ蓄積部、121 観測値記憶部、130 観測データ管理部、140 航跡候補生成部、141 候補算出部、142 候補記憶部、150 相関ゲート設定部、151 誤差上限初期化部、152 誤差上限値記憶部、153 ゲート行列算出部、160 概探索部、161 相関部、162 尤度算出部、163 選択部、170 詳細探索部、171 推定値記憶部、172 推定値更新部、173 収束判定部、174 誤差上限値更新部、175 終了判定部、176 推定値出力部、700 目標、800 追尾システム、810 センサ、911 処理装置、912 入力装置、913 出力装置、914 記憶装置。

【特許請求の範囲】
【請求項1】
データを処理する処理装置と、候補算出部と、相関部と、尤度算出部と、選択部と、推定値更新部と、推定値出力部とを有し、
上記候補算出部は、目標の複数の状態量のうち所定の第一の状態量を除く他の状態量について、上記処理装置を用いて、上記他の状態量の推定値の複数の候補値を算出し、
上記相関部は、上記候補算出部が算出した上記他の状態量の推定値の複数の候補値それぞれについて、上記処理装置を用いて、上記第一の状態量が所定の推定値であり、上記第一の状態量が上記推定値であると仮定することによって生じるモデル化誤差が所定のモデル化誤差上限値以下であり、上記他の状態量が上記候補値であるとの仮定に基づいて、上記目標を観測した観測値を含む複数の観測値のなかから、上記目標を観測した観測値を抽出し、
上記尤度算出部は、上記候補算出部が算出した上記他の状態量の推定値の複数の候補値それぞれについて、上記処理装置を用いて、上記相関部が仮定したモデル化誤差上限値と、上記相関部が抽出した観測値とに基づいて、上記他の状態量が上記候補値である尤度を算出し、
上記候補選択部は、上記処理装置を用いて、上記尤度算出部が算出した尤度に基づいて、上記候補算出部が算出した複数の候補値のなかから、上記他の状態量の推定値を選択し、
上記推定値更新部は、上記処理装置を用いて、上記相関部が仮定した第一の状態量の推定値と上記候補選択部が選択した他の状態量の推定値とを出発点として、上記目標の複数の状態量が推定値である尤度が高くなるよう、上記複数の状態量の推定値を更新し、
上記推定値出力部は、上記処理装置を用いて、上記推定値更新部が更新した上記複数の状態量の推定値を出力することを特徴とする推定装置。
【請求項2】
上記推定装置は、更に、収束判定部と、誤差上限値更新部と、詳細相関部と、詳細尤度算出部とを有し、
上記収束判定部は、上記処理装置を用いて、上記推定値更新部が更新した上記複数の状態量の推定値が収束したか否かを判定し、
上記誤差上限値更新部は、上記推定値更新部が上記複数の状態量の推定量を更新した場合に、上記処理装置を用いて、上記モデル化誤差上限値を小さくし、
上記詳細相関部は、上記処理装置を用いて、上記複数の状態量が、上記推定値更新部が更新した推定値であり、上記モデル化誤差が、上記誤差上限値更新部が小さくしたモデル化誤差上限値以下であるとの仮定に基づいて、上記複数の観測値のなかから上記目標を観測した観測値を抽出し、
上記詳細尤度算出部は、上記推定値更新部が更新した上記複数の状態量の推定値について、上記処理装置を用いて、上記誤差上限値更新部が小さくしたモデル化誤差上限値と、上記詳細相関部が抽出した観測値とに基づいて、上記複数の状態量が上記推定値である尤度を算出し、
上記推定値更新部は、上記推定値更新部が更新した上記複数の状態量の推定値が収束していないと上記収束判定部が判定した場合に、上記処理装置を用いて、更新した上記複数の状態量の推定値を出発点として、上記詳細尤度算出部が算出する尤度が更に高くなるよう、上記複数の状態量の推定値を更新し、
上記推定値出力部は、上記推定値更新部が更新した上記複数の状態量の推定値が収束したと上記収束判定部が判定した場合に、上記処理装置を用いて、上記推定値更新部が更新した上記複数の状態量の推定値を出力することを特徴とする請求項1に記載の推定装置。
【請求項3】
上記推定装置は、更に、収束判定部と、終了判定部と、誤差上限値更新部と、詳細相関部と、詳細尤度算出部とを有し、
上記収束判定部は、上記処理装置を用いて、上記推定値更新部が更新した上記複数の状態量の推定値が収束したか否かを判定し、
上記終了判定部は、上記推定値更新部が更新した上記複数の状態量の推定値が収束したと上記収束判定部が判定した場合に、上記処理装置を用いて、上記複数の状態量の推定を終了するか否かを判定し、
上記誤差上限値更新部は、上記複数の状態量の推定を終了しないと上記終了判定部が判定した場合に、上記処理装置を用いて、上記モデル化誤差上限値を小さくし、
上記詳細相関部は、上記処理装置を用いて、上記推定値更新部が更新した上記複数の状態量の推定値と、上記相関部が仮定し若しくは上記誤差上限値更新部が小さくしたモデル化誤差上限値とに基づいて、上記複数の観測値のなかから上記目標を観測した観測値を抽出し、
上記詳細尤度算出部は、上記推定値更新部が更新した上記複数の状態量の推定値について、上記処理装置を用いて、上記相関部が仮定し若しくは上記誤差上限値更新部が小さくしたモデル化誤差上限値と、上記詳細相関部が抽出した観測値とに基づいて、上記複数の状態量が上記推定値である尤度を算出し、
上記推定値更新部は、上記推定値更新部が更新した上記複数の状態量の推定値が収束していないと上記収束判定部が判定した場合、および、上記複数の状態量の推定を終了しないと上記終了判定部が判定した場合に、上記処理装置を用いて、更新した上記複数の状態量の推定値を出発点として、上記詳細尤度算出部が算出する尤度が更に高くなるよう、上記複数の状態量の推定値を更新し、
上記推定値出力部は、上記複数の状態量の推定を終了すると上記終了判定部が判定した場合に、上記処理装置を用いて、上記推定値更新部が更新した上記複数の状態量の推定値を出力することを特徴とする請求項1に記載の推定装置。
【請求項4】
上記詳細相関部は、上記処理装置を用いて、上記推定値更新部が更新した上記複数の状態量の推定値と、上記相関部が仮定し若しくは上記誤差上限値更新部が小さくしたモデル化誤差上限値とに基づいて、上記複数の観測値が上記目標を観測した観測値であるかと判定する範囲を算出するために使用するゲート行列を算出し、
上記終了判定部は、上記処理装置を用いて、上記詳細相関部が算出したゲート行列に基づいて、上記複数の状態量の推定を終了するか否かを判定することを特徴とする請求項3に記載の推定装置。
【請求項5】
上記終了判定部は、上記複数の状態量の推定を終了するか否かを判定した回数が所定の回数を超えた場合に、上記処理装置を用いて、上記複数の状態量の推定を終了すると判定することを特徴とする請求項3に記載の推定装置。
【請求項6】
上記誤差上限値更新部は、上記処理装置を用いて、次のいずれかの式の右辺を計算することにより、上記モデル化誤差上限値を更新することを特徴とする請求項2乃至請求項5のいずれかに記載の推定装置。
max←√[(amax−Δa
ただし、amaxは、上記モデル化誤差上限値を表わすベクトルのi番目の成分。Δaは、上記推定値更新部が更新した後の第一の状態量の推定値を表わすベクトルのi番目の成分と、上記推定値更新部が更新する前の第一の状態量の推定値を表わすベクトルのi番目の成分との差。
max←λ・amax
ただし、λは、モデル化誤差上限値の逓減率であり、0超1未満の実数。
【請求項7】
上記相関部及び上記詳細相関部は、上記処理装置を用いて、次のいずれかの式の右辺を計算することにより、ゲート行列を算出し、算出したゲート行列を使用して、上記複数の観測値が上記目標を観測した観測値であると判定する範囲を算出し、上記複数の観測値のうち、算出した範囲内に入る観測値を、上記目標を観測した観測値であると判定することを特徴とする請求項2乃至請求項6のいずれかに記載の推定装置。
S←R+Q
ただし、Sは、ゲート行列であり、上記観測値を表わすベクトルの成分の数を次数とする正方行列。Rは、観測雑音共分散行列であり、上記観測値を表わすベクトルの成分の数を次数とし、上記観測値を表わすベクトルの各成分の観測精度の二乗を対角成分とする対角行列。Qは、モデル化誤差共分散行列であり、上記観測値を表わすベクトルの成分の数を次数とし、上記モデル化誤差上限値を表わすベクトルの各成分から算出される値を対角成分とする対角行列。
S←R+Q+Q’
ただし、Q’は、量子化誤差共分散行列であり、上記観測値を表わすベクトルの成分の数を次数とする対称行列。
【請求項8】
上記推定装置は、複数の時刻における上記目標の位置を、上記目標を観測した観測値として入力し、所定の基準時刻における上記目標の位置及び速度及び加速度を、上記複数の状態量として推定し、
上記相関部は、上記処理装置を用いて、上記目標の加速度を上記第一の状態量として、上記加速度の推定値及びモデル化誤差上限値を仮定し、
上記相関部及び上記詳細相関部は、上記処理装置を用いて、次の式の右辺を計算することにより、上記モデル化誤差共分散行列の対角成分を算出することを特徴とする請求項7に記載の推定装置。
(t)←(amax・Δt/2)
ただし、q(t)は、時刻tにおけるモデル化誤差共分散行列Qのi行i列成分。amaxは、上記加速度のモデル化誤差上限値を表わすベクトルのi番目の成分。Δtは、上記基準時刻から上記時刻tまでの経過時間。
【請求項9】
上記推定装置は、複数の時刻における観測値を入力し、所定の基準時刻における複数の状態量を推定し、
上記相関部及び上記詳細相関部は、上記処理装置を用いて、上記基準時刻における上記目標の複数の状態量の推定値または推定値の候補に基づいて、上記複数の時刻それぞれにおける上記目標の複数の状態量を予測し、上記複数の観測値のうち次の不等式を満たす観測値を、上記目標を観測した観測値であると判定することを特徴とする請求項7または請求項8に記載の推定装置。
[z−h(x)]・S−1・[z−h(x)]≦d
ただし、zは、上記観測値を表わす縦ベクトル。hは、観測関数を表わす。xは、時刻tにおける上記目標の複数の状態量の予測値を表わす縦ベクトル。上付きのTは、行列の転置を表わす。Sは、上記時刻tにおけるゲート行列。上付きの−1は、逆行列を表わす。dは、ゲートの大きさを表わすパラメータであり、0超の実数。
【請求項10】
上記推定装置は、複数の時刻における観測値を入力し、所定の基準時刻における複数の状態量を推定し、
上記尤度算出部及び詳細尤度算出部は、上記処理装置を用いて、次の式の右辺を計算することにより、尤度を算出することを特徴とする請求項7乃至請求項9のいずれかに記載の推定装置。
【数1】

ただし、φ(Z|x)は、上記複数の状態量が上記推定値である尤度。lnは、自然対数を表わす。Pは、目標の探知確率を表わすパラメータであり、0超1未満の実数。βは、不要信号密度を表わすパラメータであり、0超の実数。detは、行列式を表わす。πは、円周率。Sは、時刻tにおけるゲート行列。mは、上記複数の観測値のうち、上記時刻tにおいて上記目標を観測した観測値であると、上記相関部または上記詳細相関部が判定した観測値の数。zk,jは、上記時刻tにおいて上記目標を観測した観測値であると、上記相関部または上記詳細相関部が判定した観測値を表わす縦ベクトル。hは、観測関数を表わす。xは、時刻tにおける上記目標の複数の状態量の予測値を表わす縦ベクトル。上付きのTは、行列の転置を表わす。上付きの−1は、逆行列を表わす。
【請求項11】
上記収束判定部は、上記処理装置を用いて、上記推定値更新部が上記複数の状態量の推定値を更新したことによる上記尤度の変化に基づいて、上記複数の状態量の推定値が収束したか否かを判定することを特徴とする請求項2乃至請求項10のいずれかに記載の推定装置。
【請求項12】
上記収束判定部は、上記複数の状態量の推定値が収束したか否かを判定した回数が所定の回数を超えた場合に、上記処理装置を用いて、上記推定値更新部が上記複数の状態量の推定値が収束したと判定することを特徴とする請求項2乃至請求項10のいずれかに記載の推定装置。
【請求項13】
上記推定値更新部は、上記処理装置を用いて、更新前よりも上記尤度が高くなる上記複数の状態量の推定値の差分を算出し、上記複数の状態量それぞれについて、算出した差分を上記状態量の推定値に加えることにより、上記複数の状態量の推定値を更新することを特徴とする請求項1乃至請求項12のいずれかに記載の推定装置。
【請求項14】
上記推定値更新部は、上記処理装置を用いて、最急降下法およびニュートン法および準ニュートン法および共役分配法のいずれかにより、更新前よりも上記尤度が高くなる上記複数の状態量の推定値の差分を算出することを特徴とする請求項13に記載の推定装置。
【請求項15】
データを処理する処理装置を有するコンピュータが実行することにより、上記コンピュータが請求項1乃至請求項14のいずれかに記載の推定装置として機能することを特徴とするコンピュータプログラム。
【請求項16】
データを処理する処理装置を有する推定装置が、目標の複数の状態量を推定する推定方法において、
上記複数の状態量のうち所定の第一の状態量を除く他の状態量について、上記処理装置が、上記他の状態量の推定値の複数の候補値を算出し、
算出した上記他の状態量の推定値の複数の候補値それぞれについて、上記処理装置が、上記第一の状態量が所定の推定値であり、上記第一の状態量が上記推定値であると仮定することによって生じるモデル化誤差が所定のモデル化誤差上限値以下であり、上記他の状態量が上記候補値であるとの仮定に基づいて、上記目標を観測した観測値を含む複数の観測値のなかから、上記目標を観測した観測値を抽出し、
算出した上記他の状態量の推定値の複数の候補値それぞれについて、上記処理装置が、仮定したモデル化誤差上限値と、抽出した観測値とに基づいて、上記他の状態量が上記候補値である尤度を算出し、
上記処理装置が、算出した尤度に基づいて、算出した複数の候補値のなかから、上記他の状態量の推定値を選択し、
上記処理装置が、仮定した第一の状態量の推定値と選択した他の状態量の推定値とを出発点として、上記目標の複数の状態量が推定値である尤度が高くなるよう、上記複数の状態量の推定値を更新し、
上記処理装置が、更新した上記複数の状態量の推定値を出力することを特徴とする推定方法。

【図1】
image rotate

【図2】
image rotate

【図3】
image rotate

【図4】
image rotate

【図5】
image rotate

【図6】
image rotate

【図7】
image rotate

【図8】
image rotate

【図9】
image rotate

【図10】
image rotate

【図11】
image rotate

【図12】
image rotate