(57)【特許請求の範囲】
【請求項１】
複数の移動対象が存在する動画像に関する観測データが入力として与えられる場合に、確率的な生成モデルを用いて、前記複数の移動対象に関する真値の推定量を示す隠れ状態量の時間発展関数と観測データを計算するための観測関数とを示すダイナミクスの、パターン数および各特徴を推定することで、前記複数の移動対象の隠れ状態量を推定する動画像処理装置による動画像処理方法であって、
前記動画像処理装置は、前記観測データ、前記移動対象それぞれの隠れ状態量、前記移動対象それぞれに関する前記観測データと前記隠れ状態量それぞれとの間の関係を示す複数の隠れ変数、前記隠れ変数それぞれの重み、および、前記ダイナミクスのハイパーパラメータ、を記憶する記憶手段と、演算手段と、を備えており、
前記演算手段は、
隠れ変数集合推定部が、前記観測データに基づいて、前記複数の隠れ変数、および、前記隠れ変数それぞれの重みを推定するステップと、
隠れ状態推定部が、前記観測データ、現在計算している時刻の１つ前の時刻の前記隠れ状態量、および、前記隠れ変数集合推定部が推定した前記複数の隠れ変数に基づいて、前記隠れ状態量を推定するステップと、
リサンプリング部が、前記隠れ変数集合推定部が推定した前記複数の隠れ変数、および、前記隠れ状態推定部が推定した前記隠れ状態量に基づいて、前記隠れ変数の一推定結果とそれに付随した隠れ状態量の推定結果を示すパーティクルを、前記隠れ変数の重みが大きいほど多く選択されるようにサンプリングするステップと、
ハイパーパラメータ更新部が、前記観測データ、前記リサンプリング部によるサンプリング結果に含まれる前記隠れ変数および前記隠れ状態量、ならびに、現在計算している時刻の１つ前の時刻の前記ハイパーパラメータに基づいて、前記ハイパーパラメータをベイズ推定により更新するステップと、
を繰り返し実行することで、前記隠れ状態推定部による前記複数の移動対象の隠れ状態量の推定を行う
ことを特徴とする動画像処理方法。
【請求項２】
前記隠れ変数集合推定部は、前記観測データに基づいて、前記複数の隠れ変数、および、前記隠れ変数それぞれの重みを推定するステップにおいて、
隠れ変数決定部が、前記１つ前の時刻の前記隠れ状態量、前記１つ前の時刻の前記複数の隠れ変数、および、前記１つ前の時刻の前記隠れ変数それぞれの重みに基づいて、前記パーティクルごとに、前記複数の隠れ変数を推定するステップと、
隠れ変数重み計算部が、前記隠れ変数決定部が推定した前記複数の隠れ変数に基づいて、前記パーティクルごとに、前記隠れ変数それぞれの重みを推定するステップと、
を実行することを特徴とする請求項１に記載の動画像処理方法。
【請求項３】
前記隠れ変数決定部は、前記１つ前の時刻の前記隠れ状態量、前記１つ前の時刻の前記複数の隠れ変数、および、前記１つ前の時刻の前記隠れ変数それぞれの重みに基づいて、前記パーティクルごとに、前記複数の隠れ変数を推定するステップにおいて、
前記１つ前の時刻の前記隠れ状態量、前記１つ前の時刻の前記複数の隠れ変数、および、現在計算している時刻の前記観測データの入力を受け付け、
移動対象数決定部が、前記複数の隠れ変数のうち、ある時刻において前記移動対象それぞれが前記動画像中に存在するか否かを表す前記隠れ変数を決定するステップと、
動作パターン番号決定部が、前記複数の隠れ変数のうち、前記移動対象の動作パターンの識別子を表す前記隠れ変数を決定するステップと、
動作パターンパラメータ決定部が、前記動作パターンのパラメータを決定するステップと、
移動対象対応関係決定部が、前記複数の隠れ変数のうち、前記観測データと前記移動対象の対応関係を表す前記隠れ変数を決定するステップと、
を実行することを特徴とする請求項２に記載の動画像処理方法。
【請求項４】
複数の移動対象が存在する動画像に関する観測データが入力として与えられる場合に、確率的な生成モデルを用いて、前記複数の移動対象に関する真値の推定量を示す隠れ状態量の時間発展関数と観測データを計算するための観測関数とを示すダイナミクスの、パターン数および各特徴を推定することで、前記複数の移動対象の隠れ状態量を推定する動画像処理装置であって、
前記観測データ、前記移動対象それぞれの隠れ状態量、前記移動対象それぞれに関する前記観測データと前記隠れ状態量それぞれとの間の関係を示す複数の隠れ変数、前記隠れ変数それぞれの重み、および、前記ダイナミクスのハイパーパラメータ、を記憶する記憶手段と、
前記観測データに基づいて、前記複数の隠れ変数、および、前記隠れ変数それぞれの重みを推定する隠れ変数集合推定部、
前記観測データ、現在計算している時刻の１つ前の時刻の前記隠れ状態量、および、前記隠れ変数集合推定部が推定した前記複数の隠れ変数に基づいて、前記隠れ状態量を推定する隠れ状態推定部、
前記隠れ変数集合推定部が推定した前記複数の隠れ変数、および、前記隠れ状態推定部が推定した前記隠れ状態量に基づいて、前記隠れ変数の一推定結果とそれに付随した隠れ状態量の推定結果を示すパーティクルを、前記隠れ変数の重みが大きいほど多く選択されるようにサンプリングするリサンプリング部、
前記観測データ、前記リサンプリング部によるサンプリング結果に含まれる前記隠れ変数および前記隠れ状態量、ならびに、現在計算している時刻の１つ前の時刻の前記ハイパーパラメータに基づいて、前記ハイパーパラメータをベイズ推定により更新するハイパーパラメータ更新部、
を有し、前記隠れ変数集合推定部、前記隠れ状態推定部、前記リサンプリング部、および、前記ハイパーパラメータ更新部による処理を繰り返すことで、前記隠れ状態推定部による前記複数の移動対象の隠れ状態量の推定を行う演算手段と、
を備えていることを特徴とする動画像処理装置。
【請求項５】
前記隠れ変数集合推定部は、
前記１つ前の時刻の前記隠れ状態量、前記１つ前の時刻の前記複数の隠れ変数、および、前記１つ前の時刻の前記隠れ変数それぞれの重みに基づいて、前記パーティクルごとに、前記複数の隠れ変数を推定する隠れ変数決定部と、
前記隠れ変数決定部が推定した前記複数の隠れ変数に基づいて、前記パーティクルごとに、前記隠れ変数それぞれの重みを推定する隠れ変数重み計算部と、
を備えることを特徴とする請求項４に記載の動画像処理装置。
【請求項６】
前記隠れ変数決定部は、
前記１つ前の時刻の前記隠れ状態量、前記１つ前の時刻の前記複数の隠れ変数、および、現在計算している時刻の前記観測データの入力を受け付け、
前記複数の隠れ変数のうち、ある時刻において前記移動対象それぞれが前記動画像中に存在するか否かを表す前記隠れ変数を決定する移動対象数決定部と、
前記複数の隠れ変数のうち、前記移動対象の動作パターンの識別子を表す前記隠れ変数を決定する動作パターン番号決定部と、
前記動作パターンのパラメータを決定する動作パターンパラメータ決定部と、
前記複数の隠れ変数のうち、前記観測データと前記移動対象の対応関係を表す前記隠れ変数を決定する移動対象対応関係決定部と、
を備えることを特徴とする請求項５に記載の動画像処理装置。
【請求項７】
コンピュータを請求項４ないし請求項６のいずれか一項に記載の動画像処理装置として機能させるための動画像処理プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、複数の移動対象が存在する動画像などのように、複数対象の隠れ状態量（真値の推定値）から出力される情報によって構成される時系列データが入力として与えられる場合に、複数対象の時間発展（物理量の値が時々刻々変化すること）を制御するダイナミクスパターンの数と各パターンの特徴を推定し、未知数の追跡対象の隠れ状態量を精度良く推定する技術に関する。
【背景技術】
【０００２】
複数の移動対象の追跡は動画像処理の研究において重要な問題の一つである。追跡とは、動画像中から検出された追跡対象の正しい位置や状態量を、その時間変化をモデリングすることによって連続的に推定することである。結果として、動画像から各対象の位置（状態）の時空間軌跡を得ることができる。動画像中の複数対象を同時に追跡することができれば、行動パターンの認識や状況の自動監視、動画像のコンテンツ理解などさまざまな用途に役立てることが可能となる。
【０００３】
従来の複数対象追跡アルゴリズム（例えば非特許文献１，２参照）は、追跡すべき追跡対象の数が未知、あるいは時間的に変動する場合においても各対象の状態軌跡を推定することが可能である。特に、非特許文献２の技術によれば、追跡対象の数の変化も確率的にモデル化しており、動画像に限らず多くの時系列データに関して追跡処理を行うことが可能である。
【０００４】
ほとんどの既存の複数対象追跡手法（非特許文献１，２など）では、追跡対象の時間発展を単一のダイナミクスモデルでモデル化している。すなわち、動画像中のあらゆる対象が同一のダイナミクスパターン（動作パターン）で移動すると仮定している。しかし、多くの移動物体がシーン中で観察される場合、各対象の特徴に応じた、複数の異なるダイナミクスが存在するという仮定の方がより自然で正しいモデル化であると予想される。なお、本明細書において、ダイナミクスとは、「隠れ状態量の時間発展関数」および「観測量を計算するための観測関数」のことを指し、詳細については後記する。
【０００５】
このような追跡モデルの設計には、ｉ）ダイナミクスパターンの数、そして、ｉｉ）各パターンの特徴、を事前に知る必要がある。これらの情報を人手によって推定することは困難なため、機械学習による解決法が必要となる。ｉ）の問題は、時系列データのクラスタリング手法（例えば非特許文献３，４など）を用いて解決できる。ｉｉ）の問題は、時系列データにおけるモデルパラメータの推定問題として解決される（例えば非特許文献５，６）。
【非特許文献１】B. Leibe, N.Cornelis, K. Cornelis  and L. Van Gool, “Dynamic 3D Scene Analysis from a Moving Vehicle”, in Proceedings of the IEEE International Conference on Computer Vision and Pattern Recognition, 2007.
【非特許文献２】S. Saerkkae, A. Vehtari and J. Lampinen, “Rao-Blackwellized Particle Filter for Multiple Target Tracking”, Information Fusion, Vol. 8, N0. 1， pp・2-15, 2007.
【非特許文献３】M. Ramoni and P. Sebastiani and P. Cohen, “Bayesian Clustering by Dynamics”, Machine Learning, Vol. 47, No. 1, pp. 91-121 ， 2002.
【非特許文献４】P. Smyth, “Clustering Sequences with Hidden Markov Models”，in Advances in Neural Information Processing Systems 9, pp. 648-654, 1997.
【非特許文献５】F. Caron， M. Davy， A. Doucet， E. Duflos and P. Vanheeghe，“Bayesian Inference for Linear Dynamic Models with DirichletProcess Mixtures”, IEEE Transactions on Signal Processing, To appear, 2007.
【非特許文献６】M. J. Cassidy and W. D. Penny, “Bayesian Nonstationary Autoregressive Models for Biomedical Signal Analysis”, IEEE Transactions on Biomedical Engineering, Vol. 49, No. 10, PP. 1142-1152, 2002.
【発明の開示】
【発明が解決しようとする課題】
【０００６】
しかしながら、既存のクラスタリング、あるいはパラメータ推定手法は、対象ごとに分離された観測系列（時系列）データにしか用いることができない。すなわち、複数の追跡対象が混在する観測系列データに対して、これらの従来手法を用いてパターンを学習することは不可能である。
【０００７】
そこで、本発明は、前記した問題を解決するためになされたものであり、複数の移動対象に対する追跡をより高精度で行うことを目的とする。
【課題を解決するための手段】
【０００８】
前記課題を解決するために、本発明は、複数の移動対象が存在する動画像に関する観測データが入力として与えられる場合に、確率的な生成モデルを用いて、前記複数の移動対象に関する真値の推定量を示す隠れ状態量の時間発展関数と観測データを計算するための観測関数とを示すダイナミクスの、パターン数および各特徴を推定することで、前記複数の移動対象の隠れ状態量を推定する動画像処理装置による動画像処理方法であって、前記動画像処理装置は、前記観測データ、前記移動対象それぞれの隠れ状態量、前記移動対象それぞれに関する前記観測データと前記隠れ状態量それぞれとの間の関係を示す複数の隠れ変数、前記隠れ変数それぞれの重み、および、前記ダイナミクスのハイパーパラメータ、を記憶する記憶手段と、演算手段と、を備えており、前記演算手段は、隠れ変数集合推定部が、前記観測データに基づいて、前記複数の隠れ変数、および、前記隠れ変数それぞれの重みを推定するステップと、隠れ状態推定部が、前記観測データ、現在計算している時刻の１つ前の時刻の前記隠れ状態量、および、前記隠れ変数集合推定部が推定した前記複数の隠れ変数に基づいて、前記隠れ状態量を推定するステップと、リサンプリング部が、前記隠れ変数集合推定部が推定した前記複数の隠れ変数、および、前記隠れ状態推定部が推定した前記隠れ状態量に基づいて、前記隠れ変数の一推定結果とそれに付随した隠れ状態量の推定結果を示すパーティクルを、前記隠れ変数の重みが大きいほど多く選択されるようにサンプリングするステップと、ハイパーパラメータ更新部が、前記観測データ、前記リサンプリング部によるサンプリング結果に含まれる前記隠れ変数および前記隠れ状態量、ならびに、現在計算している時刻の１つ前の時刻の前記ハイパーパラメータに基づいて、前記ハイパーパラメータをベイズ推定により更新するステップと、を繰り返し実行することで、前記隠れ状態推定部による前記複数の移動対象の隠れ状態量の推定を行うことを特徴とする。
【０００９】
かかる発明によれば、複数の隠れ変数および隠れ変数それぞれの重みの推定、隠れ状態量の推定、隠れ変数の重み等に基づいたパーティクルのサンプリング、ハイパーパラメータの更新、の各処理を繰り返し実行することで、複数の移動対象に対する追跡をより高精度で行うことができる。
【００１０】
また、本発明は、隠れ変数集合推定部が、前記観測データに基づいて、前記複数の隠れ変数、および、前記隠れ変数それぞれの重みを推定するステップにおいて、隠れ変数決定部による、前記１つ前の時刻の前記隠れ状態量、前記１つ前の時刻の前記複数の隠れ変数、および、前記１つ前の時刻の前記隠れ変数それぞれの重みに基づいて、前記パーティクルごとに、前記複数の隠れ変数を推定するステップと、隠れ変数重み計算部による、前記隠れ変数決定部が推定した前記複数の隠れ変数に基づいて、前記パーティクルごとに、前記隠れ変数それぞれの重みを推定するステップと、を実行することが望ましい。
【００１１】
かかる発明によれば、隠れ変数決定部がパーティクルごとに複数の隠れ変数を、また、隠れ変数重み計算部がパーティクルごとに隠れ変数それぞれの重みを、それぞれ推定することで、パーティクルの数だけ異なる推定仮説（および各重み）をもつことができる。したがって、移動対象が急激に（瞬間的に）運動方向を変化させた場合などでも、対応する推定仮説を使うことで、複数の移動対象をより高精度で追跡することができる。
【００１２】
また、本発明は、前記隠れ変数決定部が、前記１つ前の時刻の前記隠れ状態量、前記１つ前の時刻の前記複数の隠れ変数、および、前記１つ前の時刻の前記隠れ変数それぞれの重みに基づいて、前記パーティクルごとに、前記複数の隠れ変数を推定するステップにおいて、前記１つ前の時刻の前記隠れ状態量、前記１つ前の時刻の前記複数の隠れ変数、および、現在計算している時刻の前記観測データの入力を受け付け、移動対象数決定部による、前記複数の隠れ変数のうち、ある時刻において前記移動対象それぞれが前記動画像中に存在するか否かを表す前記隠れ変数を決定するステップと、動作パターン番号決定部による、前記複数の隠れ変数のうち、前記移動対象の動作パターンの識別子を表す前記隠れ変数を決定するステップと、動作パターンパラメータ決定部による、前記動作パターンのパラメータを決定するステップと、移動対象対応関係決定部による、前記複数の隠れ変数のうち、前記観測データと前記移動対象の対応関係を表す前記隠れ変数を決定するステップと、を実行することが望ましい。
【００１３】
かかる発明によれば、移動対象数決定部が移動対象の存在を表す隠れ変数を決定することで、移動対象の増減に対応することができる。また、動作パターン番号決定部が移動対象の動作パターンの識別子を表す隠れ変数を決定することで、移動対象が動作パターンを切り替えながら行動する状況に対応することができる。さらに、動作パターンパラメータ決定部が動作パターンのパラメータを決定することで、動作パターンの初期パラメータが不適切であることによって追跡精度が低下する事態を回避できる。また、移動対象対応関係決定部が観測データと移動対象の対応関係を表す隠れ変数を決定することで、観測データと移動対象の対応関係を把握することができる。
【００１４】
また、本発明に係る動画像処理プログラムは、コンピュータを動画像処理装置として機能させることを特徴とする。このような構成により、このプログラムをインストールされたコンピュータは、このプログラムに基づいた各機能を実現することができる。
【発明の効果】
【００１５】
本発明によれば、複数の移動対象に対する追跡をより高精度で行うことができる。
【発明を実施するための最良の形態】
【００１６】
以下、図面を参照（言及図以外の図も適宜参照）して、本発明を実施するための最良の形態（以下、「実施形態」という。）について詳細に説明する。本実施形態では、動画像中に存在する複数の対象（移動対象または追跡対象ともいう。）を同時に追跡する手法について説明する。なお、前記したように、ダイナミクスとは、「隠れ状態量（後記するｘ（ｔ））の時間発展関数（後記する式（２０））」および「観測量（後記するｙ（ｔ））を計算するための観測関数（後記する式（２１））」のことを指す。また、ダイナミクスの学習とは、式（２０）、式（２１）に必要となるパラメータ（後記するξ，ψ）を推定することである。
【００１７】
従来手法と本実施形態による手法（以下、本手法という。）について改めて説明すると、従来手法では、複数対象の追跡に先んじて対象の状態変化を表現するダイナミクスモデルを決定しておく必要があるが、一般に対象のダイナミクスの数や各ダイナミクスのパラメータの値は未知であり、これらを自動的に推定して設定することは困難であった。一方、対象のダイナミクスを推定する手法も存在するが、これらの手法は複数対象が存在するデータに適用不能であった。本実施形態では、複数対象が存在するデータに対して、確率的な生成モデルを用いてダイナミクスパターンの数とパラメータを推定すると同時に未知数の対象に対する追跡を行うことができる手法について説明する。
【００１８】
（動画像処理装置の構成）
まず、本実施形態の動画像処理装置の構成について説明する。図１は、本実施形態に係る動画像処理装置の構成を模式的に示す機能ブロック図である。動画像処理装置１は、例えば、ＣＰＵ（Central Processing Unit）、ＲＡＭ（Random Access Memory）、ＲＯＭ（Read Only Memory）、ＨＤＤ（Hard Disk Drive）、入出力インタフェースなどから構成されるコンピュータ装置である。図１に示すように、動画像処理装置１は、演算手段２（例えばＣＰＵとＲＡＭ）、記憶手段３（例えばＲＯＭとＨＤＤ）、キーボードやマウスなどの入力手段４、および、液晶表示機などの出力手段５を備え、それらがバスライン６で接続されている。なお、ここでは動画像処理装置１の構成の概要について説明し、詳細については図２Ａ以下で説明する。
【００１９】
記憶手段３は、動画像に関する観測量（観測データまたは観測値ともいう。）３１、動画像におけるそれぞれの移動対象に関する真値の推定量を示す複数の隠れ状態量３２、それぞれの移動対象に関する観測データと隠れ状態量（隠れ状態ともいう。）との間の関係を示す複数の隠れ変数３３、隠れ変数それぞれの重み３４、および、ダイナミクスのハイパーパラメータを含む各種パラメータ３５、を記憶する（詳細は後記）。なお、記憶手段３は、特に図示しないが、演算手段２の動作プログラムなども記憶する。また、以下において、演算手段２が記憶手段３に対して各データの読み書きを行う場合、「記憶手段３に対して」などの記載は省略する。
【００２０】
演算手段２は、隠れ変数集合推定部２１、隠れ状態推定部２２、リサンプリング部２３、および、ハイパーパラメータ更新部２４を備えている。
【００２１】
隠れ変数集合推定部２１は、観測データに基づいて、複数の隠れ変数、および、それぞれの隠れ変数の重みを推定するものであり、隠れ変数決定部２１１と隠れ変数重み計算部２１２とを備えている。
【００２２】
隠れ変数決定部２１１は、１つ前の時刻の隠れ変数、複数の隠れ状態量、および、それぞれの隠れ変数の重みに基づいて、パーティクルごとに、隠れ変数を推定するものであり、移動対象数決定部２１１１、動作パターン番号決定部２１１２、動作パターンパラメータ決定部２１１３、および、移動対象対応関係決定部２１１４を備えている。
【００２３】
移動対象数決定部２１１１は、ある時刻において各移動対象が動画中に存在するか否かを表す隠れ変数を決定する。動作パターン番号決定部２１１２は、移動対象の動作パターンの識別子（識別番号など）を表す隠れ変数を決定する。
【００２４】
動作パターンパラメータ決定部２１１３は、動作パターンのパラメータを決定する。移動対象対応関係決定部２１１４は、観測データと移動対象の対応関係を表す隠れ変数を決定する。
【００２５】
隠れ変数重み計算部２１２は、隠れ変数決定部２１１が推定した隠れ変数に基づいて、パーティクル（詳細は後記）ごとに、それぞれの隠れ変数の重みを推定する。
【００２６】
隠れ状態推定部２２は、観測データ、１つ前の時刻の隠れ状態量、および、隠れ変数集合推定部２１が推定した複数の隠れ変数に基づいて、隠れ状態量を推定する。
【００２７】
リサンプリング部２３は、隠れ変数集合推定部２１が推定した複数の隠れ変数、および、隠れ状態推定部２２が推定した隠れ状態量に基づいて、隠れ変数の一推定結果とそれに付随した隠れ状態量の推定結果を示すパーティクルの重みの分布から事後分布の高い部分を推測してサンプリングを複数回行うリサンプリングを行う。
【００２８】
ハイパーパラメータ更新部２４は、リサンプリング部２３によるリサンプリング結果に基づいて、ダイナミクスのパラメータを更新する。
【００２９】
（計算処理の概要と目的）
次に、本手法の計算処理の概要とその目的について説明する。図２Ａは、複数対象追跡のイメージ図である。観測データの一部は追跡すべき対象と関係ないノイズであることもある。本手法では、複数の追跡対象の状態ベクトルｘ（ｔ）（例えば２次元の位置ベクトル）を、観測された時系列データ（観測データ）Ｙ（ｔ）＝｛ｙ（ｔ）｝，ｔ＝１，２，・・・，から推定することが目標である。
【００３０】
ｘ（ｔ）は隠れ状態量（すなわち直接の観測は不可能）であり、また、計算を簡単にするため、ｘ（ｔ）にはマルコフ性（将来における事象の起こる確率は、現在の状態だけから決まり、過去の状態に依存しないという性質）を仮定する。なお、ｘ（ｔ）は直接の観測が不可能である、というのは、あるデータが観測されたとしてもそのデータはノイズによって真値からずれている可能性があり正確なデータとは限らない、という意味である。つまり、隠れ状態量とはこの真値（またはその推定値）のことを指す。そして、真値は、観測不能であるが、本手法によれば高精度で推定することができる。
【００３１】
さらに、観測データと隠れ状態量の間の関係を表現するために隠れ変数φ（ｔ）を定義する。実際の計算では、各時刻ｔにおける隠れ状態量ｘ（ｔ）を推定するために、観測データＹ（ｔ）と時刻ｔ−１での隠れ変数φ（ｔ−１）および隠れ状態量ｘ（ｔ−１）を用いて、隠れ変数φ（ｔ）の値の候補を確率的に生成し、各候補値のもとでの隠れ状態量ｘ（ｔ）を推定する。複数の隠れ変数と隠れ状態量の組み合わせ候補を、その確からしさで重みづけすることによって各時刻での隠れ状態量の確率分布を推定することを繰り返す。
【００３２】
一例をあげて説明する。今、目的がビデオ画像中の人物追跡であるとすると、隠れ状態量は各時刻における人物の位置、あるいは速度ベクトルである。観測量は動画像データあるいはそこから抽出された特徴量ベクトルの時系列データである。隠れ変数は、各時刻におけるシーン（動画像）中の人物数、各人物の動き方パターン、観測量と人物の対応関係（あるピクセルあるいは特徴ベクトルがどの人物に対応するか）といった情報である。各時刻ｔにおいて得られた観測情報と時刻ｔ−１での人物数や各人物の位置といった情報から、時刻ｔにおける隠れ変数と各人物の状態の候補値を複数推定し、その確からしさの重みづけによって隠れ状態量の分布を推定することを繰り返す。
【００３３】
（生成モデル）
続いて、本手法の計算に用いられる生成モデルについて説明する。図２Ｂは、本手法の計算に用いられる生成モデルのグラフィカルモデルである。図２Ｂでは、確率変数間の依存関係と、それらの変数から時系列データが生成される過程を表現している。生成されるデータは追跡対象の状態（位置）を表すｘ（ｔ）と、状態量にノイズの重畳した観測量であるｙ（ｔ）である。この生成モデルはダイナミクスパラメータの混合モデルに従う複数対象の行動データを表現できる。
【００３４】
複数の追跡対象が存在するため、隠れ状態量ｘ（ｔ）はｉ、観測量ｙ（ｔ）はｍの添字（インデックス）を持つ。ｉ番目の対象（隠れ状態量）をｘｉ（ｔ）、ｍ番目の観測量をｙｍ（ｔ）とする。Ｎｔは時刻ｔにおける追跡対象の数、Ｍｔは観測量の数を表す。複数の追跡対象（隠れ状態量）と観測量が存在する場合、各対象ごとの隠れ状態量の推定には、どの対象がどの観測量を出力したのか（対応するのか）、という対応関係を解決しなければならない。この問題はData Associationと呼ばれており、隠れ状態量の正確な推定のためには重要な問題である。この生成モデルでは、Data Association問題を解決するため、ｃおよびｊという隠れ変数を導入している。
【００３５】
各隠れ変数について説明する。まず、ｃｉ（ｔ）は追跡対象ｉのbirth(addition)とdeath(deletion)を表す変数である。すなわち、どれだけの対象が現在現れ、あるいは消えたかを説明する。ｃｉ（ｔ）＝１のとき、ｉ番目の対象は生存、すなわちシーン中に存在することを表し、逆にｃｉ（ｔ）＝０の場合は対象がシーン外に存在することを示す。時刻ｔに新しい追跡対象が生成された場合には、新しいインデックスｉ＾（本明細書において、記号「＾」は直前の文字の上に位置する記号であるものとする。）を導入して、ｃｉ＾（ｔ）＝１とする。ｊｍ（ｔ）は、観測量ｙｍ（ｔ）をその出力元である対象ｘｉ（ｔ）と関連づけるdata association変数である。ｊｍ（ｔ）＝ｉはｉ番目の追跡対象ｘｉ（ｔ）がｍ番目の観測量ｙｍ（ｔ）を出力したことを示す。
【００３６】
パラメータは添字ｋでインデキシングされた集合｛ξｋ｝、｛ψｋ｝となっている。ｘ（ｔ）（ｙ（ｔ））はパラメータξｋ（ψｋ）で特徴づけられた分布に従って生成される。本モデルでは状態空間モデルとしてKalman filter（カルマンフィルタ。例えば「北川源四郎、“時系列解析入門”、岩波書店、pp. 125-141, 209-222、2005.」（以下、「非特許文献７」という。）参照）を用いるので、ξｋとψｋはその正規分布パラメータである。なお、ここではKalman filterを用いたが、これに限定されず、一般に正規分布で特徴づけられるモデルが利用可能である。
【００３７】
Ｇ０およびＨ０はNormal Inverse Wishart distributions（ＮＩＷ）（例えば「渡部洋、“ベイズ統計学入門”、福村出版、1999.」参照）とし、そのパラメータ（ハイパーパラメータ）をそれぞれθξ、θψとする。ここでＮＩＷ分布を用いたのは正規分布に対する共役分布であり計算が簡略化できるためであり、その他の分布を用いることも可能である。
【００３８】
各時刻ｔにおいて、まず非特許文献５の手法に従ってｉ番目の追跡対象を支配するダイナミクスパターンのインデックスｚｉ（ｔ）を生成する。この過程はパラメータγを持つChinese Restaurant Process（ＣＲＰ）（例えば「D. Blackwell and J. MacQueen, “Ferguson Distributions via Polyaurnschemes”, The Annals of Statistics, Vol. 1, pp. 353-355, 1973.」（以下、「非特許文献８」という。）および「C. Kemp, J. B. Tenenbaum, T. L. Griffiths, T. Yamada and N. Ueda, “Learning Systems of Concepts With An Infinite Relational Model”, in Proceedings of 21st National Conference on Artificial Intelligence, 2006.」参照）とπでパラメタライズされるベルヌイ試行でモデル化する。
【００３９】
ＣＲＰはDirichlet Process Mixture（ＤＰＭ）の実現例の一つであり、クラスタクリングの事前分布を与える。ＣＲＰを用いることでダイナミクスのパターン数（クラスタ数）を事前に固定することなく、ダイナミクスのパターン化（クラスタリング）をモデル化可能である。ｚｉ（ｔ）＝ｋが生成されたら、それに対応したパラメータξｋ、ψｋを用いて隠れ状態量と観測量が生成される。このモデルの特徴は、パラメタライズされた分布を表現するためにＤＰＭを用いた点である。ＤＰＭは、ノンパラメトリックベイズモデルの一つで、混合数が無限の混合モデルと理解できる。ＤＰＭを用いることのメリットは、従来の混合モデルの推定のように各混合コンポーネントの性質（パラメータ）推定だけでなく、混合数についてもベイズ的枠組で同時に推定できる点にある。
【００４０】
以上で、ｘｉ（ｔ）、ｙｍ（ｔ）の生成に必要なすべての隠れ変数が準備された。隠れ状態量ｘｉ（ｔ）は、対象ｉがシーン中に存在する（ｃｉ（ｔ）＝１）場合のみ生成される。生成過程はインデックスｚｉ（ｔ）＝ｋで決定されたξｋとψｋで特徴づけられる。続いて、観測量ｙｍ（ｔ）の生成のために、data association変数ｊｍ（ｔ）を生成する。ｍ番目の観測量ｙｍ（ｔ）はｊｍ（ｔ）＝ｉのときｘｉ（ｔ）と観測分布のパラメータψｚｉ（ｔ）によって生成される。このプロセスをすべてのｔについて繰り返すことで、複数の異なるダイナミクスで支配された、複数対象の行動軌跡の重ね合わせからなる、一般の複数対象の行動時系列データを生成することが可能である。
【００４１】
（計算アルゴリズム）
次に、本手法の計算アルゴリズムについて説明する。最終的な目標は、各時刻において事後分布ｐ（ｘ（ｔ）｜Ｙ（ｔ））を計算することであり、例えば次の式（１）によって計算できる。
【数１】
0004972016.tif 000002
【００４２】
図３は、本手法の全体の処理を示すフローチャートである。つまり、図３は、隠れ状態量、および隠れ変数を推定するアルゴリズムの概略を示している。図１の動画像処理装置１の演算手段２は、時刻ｔ＝１（ステップＳ１）から時刻ｔ＝Ｔ（ステップＳ８でＮｏ）まで、時刻を「１」ずつ更新しながら（ステップＳ９）、ステップＳ２〜Ｓ７の処理を行う。
【００４３】
各時刻ｔにおいて、演算手段２がその時刻における観測値を取得（入力）した（ステップＳ２）のち、隠れ変数集合推定部２１は各隠れ変数（集合）を推定する（ステップＳ３：詳細は図４で後記）。この際、複数の推定結果をそれらの重みとともに保持することで、隠れ変数の推定分布を表現する。次に、隠れ状態推定部２２はその隠れ変数のもとで想定される最良の状態量を計算する（時刻ｔの隠れ状態（量）を推定）（ステップＳ４：詳細は図１１で後記）。これらの推定された状態量は隠れ変数の結果に依存する。前記したように、隠れ変数の一推定結果と、それに付随した隠れ状態量の推定結果を「パーティクル」と呼ぶ。
【００４４】
続いて、リサンプリング部２３は、パーティクルの重みに従って、リサンプリングと呼ばれる処理を行う（ステップＳ５：詳細は図１２で後記）。リサンプリングとは、パーティクルおよびその重みの分布から、事後分布の高い部分を推測し、そこからサンプリングを何度も行うことである。ハイパーパラメータ更新部２４は、このリサンプリングの結果を受け、この時刻の観測データと、ステップＳ４で推定された隠れ状態量に基づいて、ダイナミクスのパラメータであるハイパーパラメータを更新する（ステップＳ６：詳細は図１３で後記）。演算手段２は、時刻ｔの隠れ状態量、隠れ変数およびハイパーパラメータを出力し（ステップＳ７）、例えば出力手段５に追跡結果表示を行う。
【００４５】
図４は、図３のステップＳ３の詳細を示すフローチャートであり、隠れ変数（集合）の推定方法を示すフローチャートである。ここでは、ｓ＝１，・・・,Ｓ個の隠れ変数の予測値と重みを計算する。まず、時刻ｔにおける各種隠れ変数をまとめて、φ（ｔ）＝｛｛ｃｉ（ｔ）｝,｛ｊｍ（ｔ）｝,｛ｚｉ（ｔ）｝,｛ξｋ（ｔ）｝,｛ψｋ（ｔ）｝｝とする。
【００４６】
ここで、各隠れ変数について簡単に説明する。ｃｉ（ｔ）は、時刻ｔにおいてｉ番目の追跡対象がシーン中に存在するかどうかを表す変数である。動画像における人物追跡の例では、時刻ｔにおいて動画像のシーン内に人物が存在するか否かを表現している。ｃ（ｔ）＝｛ｃｉ（ｔ）｝,ｉ＝１,・・・,Ｎｔの値によって現在の追跡対象数が決定できる。
【００４７】
ｊｍ（ｔ）は、時刻ｔにおいてｍ番目の観測量を生成した追跡対象の番号を表す。動画像における人物追跡の例では、例えばｍ番目のピクセルがどの人物に対応するか、というような情報である。ｊ（ｔ）＝｛ｊｍ（ｔ）｝,ｍ＝１,・・・Ｍｔの値は各追跡対象の隠れ状態量の推定値に大きな影響を与える。
【００４８】
ｚｉ（ｔ）は、時刻ｔにおいてｉ番目の追跡対象の状態量変化を決定するダイナミクスのパターン番号である。動画像における人物追跡の例では、各人物の時刻ｔにおける動きのパターン、例えば走っている、右向きに歩く、というような情報である。ｚ（ｔ）＝｛ｚｉ（ｔ）｝，ｉ＝１，・・・，Ｎｔの値も、各追跡対象の隠れ状態量の推定値、およびｊ（ｔ）の決定に影響を与える。
【００４９】
ξｋ（ｔ）およびψｋ（ｔ）は、時刻ｔにおいてｋ番目のダイナミクスパターンを表す。より正確には、各パターンを特徴づけるパラメータである。動画像における人物追跡の例では、人物の位置ベクトルの時間変化量、あるいは観測ノイズの平均や分散に対応する。
【００５０】
隠れ変数集合推定部２１は、時刻ｔにおける観測量ｙ（ｔ）を入力し（ステップＳ３１）、ｓ＝１（ステップＳ３２）からｓ＝Ｓまで（ステップＳ３６でＮｏ）、ｓの値を「１」ずつ更新しながら（ステップＳ３７）、ステップＳ３３〜Ｓ３５の処理を行う。
【００５１】
ステップＳ３３において、隠れ変数集合推定部２１は、時刻ｔ−１における隠れ状態量ｘ（ｓ）（ｔ−１）、隠れ変数φ（ｓ）（ｔ−１）、重みｗ（ｓ）（ｔ−１）を入力する。
【００５２】
ステップＳ３４において、隠れ変数決定部２１１は、時刻ｔにおけるｓ番目の隠れ変数φ（ｓ）（ｔ）を決定する（詳細は図５で後記）。具体的には、ｓ番目のパーティクルに対してφ（ｓ）（ｔ）を次の式（２）に示すようにサンプリングすることで決定する。
【数２】
0004972016.tif 000003
【００５３】
ステップＳ３５において、隠れ変数重み計算部２１２は、ｓ番目の隠れ変数の重みｗ（ｓ）（ｔ）を計算する（詳細は図１０で後記）。具体的には、各々のパーティクルは、式（３）に示すように、ｗ（ｓ）（ｔ）によって重みづけられる。
【数３】
0004972016.tif 000004
次の式（４）〜（８）に示すように、条件付き分布ｐ（φ（ｔ）｜Φ（ｔ−１））を定義できる。
【数４】
0004972016.tif 000005
【００５４】
ステップＳ３６でＮｏの場合、隠れ変数重み計算部２１２は、時刻ｔにおけるＳ個の隠れ変数φ（ｔ）と重みｗ（ｔ）を出力する（ステップＳ３８）。
【００５５】
図５は、図４のステップＳ３４の詳細を示すフローチャートであり、隠れ変数の予測値の計算方法を示すフローチャートである。図５は、前記した提案分布からのφ（ｓ）（ｔ）のサンプリング過程であり、次の式（９）〜（１３）のように表記できる。つまり、図５は、ｑ（・）からのサンプリングを表現している。
【数５】
0004972016.tif 000006
【００５６】
隠れ変数決定部２１１は、時刻ｔ−１における隠れ状態量ｘ（ｔ−１）、隠れ変数φ（ｔ−１）を入力し（ステップＳ３４１）、時刻ｔにおける観測量ｙ（ｔ）を入力する（ステップＳ３４２）。
【００５７】
ステップＳ３４３において、移動対象数決定部２１１１は、追跡対象の生成と消滅を表す隠れ変数ｃ（ｔ）を決定する（詳細は図６で後記）。
【００５８】
ステップＳ３４４において、動作パターン番号決定部２１１２は、追跡対象の動作パターンの番号を表す隠れ変数ｚ（ｔ）を決定する（詳細は図７で後記）。
【００５９】
ステップＳ３４５において、動作パターンパラメータ決定部２１１３は、動作パターンのパラメータを表すξ（ｔ），ψ（ｔ）を決定する（詳細は図３，４，５で後記）。
【００６０】
ステップＳ３４６において、移動対象対応関係決定部２１１４は、観測値と追跡対象の対応関係を表す隠れ変数ｊ（ｔ）を決定する（詳細は図９で後記）。
【００６１】
ステップＳ３４７において、隠れ変数決定部２１１は、時刻ｔにおける隠れ変数φ（ｔ）を出力する。
【００６２】
図６は、図５のステップＳ３４３の詳細を示すフローチャートであり、式（９）を用いた追跡対象のbirth(addition)およびdeath(deletion)を表現する隠れ変数ｃ（ｔ）のサンプリングである。
【００６３】
まず、図６の処理の概要について説明する。非特許文献２において、ｃ（ｔ）の時間発展は２段階のベルヌイ試行としてモデル化されており、本モデルでもそれを踏襲する。すなわち、ｉ）前時刻においてシーン中に存在する追跡対象（ｃｉ（ｔ−１））＝１）は確率Ｐｄでシーンから消える（ｃｉ（ｔ）＝０）。そうでなければｃｉ（ｔ）＝１とする。ｉｉ）新たな追跡対象が確率Ｐｂで生成される。各時刻で新たに生成される追跡対象の数を「１」に制限（非特許文献２参照）すると、全追跡対象についてｃ（ｔ）の計算を行ったときに、その値を得る確率は式（１４）のようになる。
【００６４】
【数６】
0004972016.tif 000007
ここで、ｎｓはシーン中に現存している追跡対象の数、ｎｄは時刻ｔにシーンから消えた対象数、そして（ｎｂ＝｛０，１｝）が新たに生成した追跡対象の数である。
【００６５】
具体的には、移動対象数決定部２１１１は、まず、ｃ（ｔ−１）を入力し（ステップＳ６０１）、ｉ＝１（ステップＳ６０２）からｉ＝Ｎｔ−１まで（ステップＳ６０７でＮｏ）、ｉの値を「１」ずつ更新しながら（ステップＳ６０８）、ステップＳ６０３〜Ｓ６０６の処理を行う。
【００６６】
移動対象数決定部２１１１は、ｃｉ（ｔ−１）＝１か否か判断し（ステップＳ６０３）、Ｙｅｓの場合は時刻ｔにおいて追跡対象ｉが存在するか否か判断する（ステップＳ６０４）。ステップＳ６０４でＹｅｓの場合はｃｉ（ｔ）に１を代入する（ステップＳ６０６）。ステップＳ６０３とステップＳ６０４のいずれかでＮｏの場合は、ｃｉ（ｔ）に０を代入する（ステップＳ６０５）。
【００６７】
ステップＳ６０９において、移動対象数決定部２１１１は、新しい追跡対象が発生したか否か判断し、Ｙｅｓの場合はＮｔにＮｔ−１＋１を代入し（ステップＳ６１０）、ｃＮｔ（ｔ）に１を代入し（ステップＳ６１１）、Ｎｏの場合はＮｔにＮｔ−１を代入する（ステップＳ６１２）。最後に、移動対象数決定部２１１１は、ｃ（ｔ）を出力する（ステップＳ６１３）。
【００６８】
図７は、図５のステップＳ３４４の詳細を示すフローチャートであり、式（１１）を用いたｚ（ｔ）のサンプリングである。まず、図７の処理の概要について説明する。ｚ（ｔ）の生成にはChinese Restaurant Process（ＣＲＰ）（非特許文献８参照）を利用する。ＣＲＰとは、ノンパラメトリックベイズモデルの一種であるDirichlet Process Mixture（ＤＰＭ）の実現例で、サンプルの分割（クラスタリング）の事前分布である。今、ｋでインデックスされた混合分布の要素（クラスタ）集合があり、ｉ番目のサンプルのクラスタインデックスをｚｉと書くことにする。ｚ１からｚｉ−１が与えられたとき、ｚｉの分布は次の式（１５）のようになる。
【数７】
0004972016.tif 000008
【００６９】
ここで、ｍｋはｋ番目のクラスタの大きさである。このサンプリングプロセスをｚｉ〜ＣＲＰ（γ）と書く。式（１５）からすぐ分かるようにＣＲＰから生成されるクラスタリングではサイズが大きい少数のクラスタが生成されやすい。これはクラスタリングの事前分布として適切な性質である。ＣＲＰによって生成されたクラスタリングに適切な尤度関数を組み合わせることで、ｚ（ｔ）の事後分布も導くことができる。
【００７０】
本モデルでは、各時刻で最良のクラスタリング推定を得るため、式（１５）をオンラインで更新する。Ｚ（ｔ−１）を事前の情報として、時刻ｔにおけるｚｉ（ｔ）のサンプリングを次の式（１６）のように行う。つまり、ここで、「オンライン」とは、「インクリメンタル」の意味であり、式（１５）の計算は、実際には式（１６）のように時刻ｔ−１までのＺの推定結果を用いて計算する、ということである。
【数８】
0004972016.tif 000009
【００７１】
ｍｋ（ｔ−１）は時刻ｔ−１までに得られたｋ番目のクラスタの大きさである。また、｜Ｚ（ｔ−１）｜は時刻ｔ−１におけるクラスタのサイズ総数、すなわちｍｋ（ｔ−１）の総和である。
【００７２】
本モデルではＣＲＰによって各対象のダイナミクスを変更することができるが、そういったダイナミクスの変更は一定の頻度でしか発生しないという場合が考えられる。ダイナミクス変更の頻度もモデル化するため新しいパラメータπを導入し、確率πでダイナミクスの変更を行うこととする。最終的に式（１１）は次のようにモデル化される。まず、式（１７）に示すように、式（１１）を対象ｉで分解する。
【数９】
0004972016.tif 000010
【００７３】
その上で、各ｚｉ（ｔ）を次の式（１８）、式（１９）のようにサンプリングする。
【数１０】
0004972016.tif 000011
【００７４】
具体的には、動作パターン番号決定部２１１２は、まず、ｚ（ｔ−１）、ｃ（ｔ）を入力し（ステップＳ７０１）、ｉ＝１（ステップＳ７０２）からｉ＝Ｎｔまで（ステップＳ７０９でＮｏ）、ｉの値を「１」ずつ更新しながら（ステップＳ７０４）、ステップＳ７０３〜Ｓ７０８の処理を行う。
【００７５】
動作パターン番号決定部２１１２は、ｃｉ（ｔ）＝１か否か判断し（ステップＳ７０３）、Ｎｏの場合はステップＳ７０４に進み、Ｙｅｓの場合は追跡対象ｉの動作パターンが変化したか否か判断する（ステップＳ７０５）。ステップＳ７０５でＹｅｓの場合は、ＣＲＰにしたがって動作パターンのインデックスｋを決定し（ステップＳ７０６）、ｚｉ（ｔ）にｋを代入する（ステップＳ７０７）。ステップＳ７０５でＮｏの場合は、ｚｉ（ｔ）にｚｉ（ｔ−１）を代入する（ステップＳ７０８）。最後に、動作パターン番号決定部２１１２は、ｚ（ｔ）を出力する（ステップＳ７１０）。
【００７６】
図８は、図５のステップＳ３４５の詳細を示すフローチャートであり、Kalman filterの正規分布パラメータであるξ（ｔ）、ψ（ｔ）をそれぞれ式（１２）、式（１３）からサンプリングするプロセスである。
【００７７】
まず、図８の処理の概要について説明する。先に述べたように、Kalman filterでなく、他のモデルでも正規分布で特徴づけられるモデルならば同様のモデル化が可能である。時刻ｔにおいてパーティクルｓが推定するクラスタ数をｋｔとする。追跡対象の隠れ状態量ｘｉ（ｔ）および観測量ｙｍ（ｔ）がｋ番目のダイナミクスクラスタから生成される場合（すなわちｊｍ（ｔ）＝ｉ  and  ｚｉ（ｔ）＝ｋの場合）、それらの生成過程はKalman filterを用いて次の式（２０）、式（２１）のようになる。
【数１１】
0004972016.tif 000012
【００７８】
ここで、ｆとｈは正規分布ノイズを持つ線形モデルである。システムノイズの平均と共分散行列をξｋ（ｔ）＝｛ｑ，Ｑ｝、観測ノイズの平均と共分散行列をψｋ（ｔ）＝｛ｒ，Ｒ｝とする。これらのパラメータ推定のため、パラメータの事前分布にNormal Inverse Wishart distribution（ＮＩＷ）（非特許文献５参照）を導入する。各クラスタのシステムノイズ、観測ノイズそれぞれにＮＩＷを設定して、パラメータをθｋ＝｛θξｋ，θψｋ｝と表す。
【００７９】
時刻ｔ−１までのデータに基づいた事後分布から、時刻ｔにおける真のパラメータξｋ，ψｋの推定値ξｋ（ｔ）とψｋ（ｔ）をサンプリングする。時刻ｔにおけるサンプリングには、時刻ｔ−１におけるθの事後推定値θｋ（ｔ−１）＝｛θξｋ（ｔ−１），θψｋ（ｔ−１）｝を用いる。式（１２）、式（１３）をそれぞれｐ（ξ（ｔ）｜θξ（ｔ−１））＝Πｋ（ξｋ（ｔ）｜θξｋ（ｔ−１）），ｐ（ψ（ｔ）｜θψ（ｔ−１））＝Πｋ（ψｋ（ｔ）｜θψｋ（ｔ−１））のように分解し、各要素を次の式（２２）、式（２３）のように表す。
【数１２】
0004972016.tif 000013
【００８０】
具体的には、動作パターンパラメータ決定部２１１３は、まず、ｚ（ｔ）、θ（ｔ−１）、θ（０）を入力し（ステップＳ８０１）、ｋ＝１（ステップＳ８０２）からｋ＝Ｋｔまで（ステップＳ８０８でＮｏ）、ｋの値を「１」ずつ更新しながら（ステップＳ８０９）、ステップＳ８０３〜Ｓ８０７の処理を行う。
【００８１】
動作パターンパラメータ決定部２１１３は、ｋは時刻ｔに生成された新しいインデックスか否か判断し（ステップＳ８０３）、Ｙｅｓの場合はξｋ（ｔ）をＮＩＷ（θξｋ（０））からサンプリングし（ステップＳ８０４）、ψｋ（ｔ）をＮＩＷ（θψｋ（０））からサンプリングする（ステップＳ８０５）。ステップＳ８０３でＮｏの場合はξｋ（ｔ）をＮＩＷ（θξｋ（ｔ−１））からサンプリングし（ステップＳ８０６）、ψｋ（ｔ）をＮＩＷ（θψｋ（ｔ−１））からサンプリングする（ステップＳ８０７）。最後に、動作パターンパラメータ決定部２１１３は、ξ（ｔ）、ψ（ｔ）を出力する（ステップＳ８１０）。
【００８２】
図９は、図５のステップＳ３４６の詳細を示すフローチャートであり、式（２５）によるdata association変数Ｊのサンプリングであるが、ここでは式（２４）についても同時に説明する。
【００８３】
まず、図９の処理の概要について説明する。各観測量が独立であると仮定し、次の式（２４）のように分解する。
【数１３】
0004972016.tif 000014
【００８４】
本モデルではｐ（ｊｍ（ｔ））を一様分布、すなわち一般にdata associationに関して事前知識（情報）がない状況を想定する。ただし、後記する実験では事前知識を導入した例を示す。
【００８５】
Ｊの探索空間が広い場合には、疑似的に尤度を導入することで外れ値を除外することも可能である（例えば「M. K. Pitt and N. Shephard, “Filtering via simulation: Auxiliary particle filters”, Journal of the American Statistical Association， Vol. 94, N0. 446， pp. 590-599， 1999.」参照）。後記する実験では、ｊのproposal distributionとして次の式（２５）のような分布を用いた（非特許文献２参照）。
【数１４】
0004972016.tif 000015
【００８６】
この式は、“代表的な”状態ベクトルｘ＾ｉ（ｔ）における尤度を使って事後分布に近い分布を計算しようとしている。本モデルではｘ＾ｉ（ｔ）として、ｘｉ（ｔ）予測分布（フィルタリング前）の平均値を利用した。
【００８７】
具体的には、移動対象対応関係決定部２１１４は、まず、時刻ｔ−１における隠れ状態量ｘ（ｔ−１）を入力し（ステップＳ９０１）、時刻ｔの観測値ｙ（ｔ）を入力し（ステップＳ９０２）、Ｊ（ｔ−１）、ｃ（ｔ）、ｚ（ｔ）、ξ（ｔ）、ψ（ｔ）を入力する（ステップＳ９０３）。
【００８８】
続いて、移動対象対応関係決定部２１１４は、ｍ＝１（ステップＳ９０４）からｍ＝Ｍｔまで（ステップＳ９１２でＮｏ）、ｍの値を「１」ずつ更新しながら（ステップＳ９１３）、ステップＳ９０５〜Ｓ９１１の処理を行う。
【００８９】
また、移動対象対応関係決定部２１１４は、ｉ＝１（ステップＳ９０５）からｉ＝Ｎｔまで（ステップＳ９１０でＮｏ）、ｉの値を「１」ずつ更新しながら（ステップＳ９０７）、ステップＳ９０６〜Ｓ９０９の処理を行う。
【００９０】
移動対象対応関係決定部２１１４は、ｃｉ（ｔ）＝１か否か判断し（ステップＳ９０６）、Ｎｏの場合はステップＳ９０７に進み、Ｙｅｓの場合はｚｉ（ｔ）、ξ（ｔ）、ψ（ｔ）からｘ＾ｉを計算し（ステップＳ９０８）、ｑ（ｊｍ（ｔ）＝ｉ｜Ｊ（ｔ−１），ｙｍ（ｔ））を計算する（ステップＳ９０９）。
【００９１】
ステップＳ９１１において、移動対象対応関係決定部２１１４は、ｑ（ｊｍ（ｔ）｜Ｊ（ｔ−１），ｙｍ（ｔ））に従ってｊｍ（ｔ）を決定する。最後に、移動対象対応関係決定部２１１４は、ｊ（ｔ）を出力する（ステップＳ９１４）。
【００９２】
次に、図４のステップＳ３５の処理、つまり、隠れ変数の重み（以下、単に「重み」ともいう。）を計算する処理について説明する。図１０は、図４のステップＳ３５の詳細を示すフローチャートである。重みの計算式は既に式（３）に示した通りである。
【００９３】
具体的には、隠れ変数重み計算部２１２は、まず、時刻ｔ−１における隠れ状態量ｘ（ｔ−１）と重みｗ（ｔ−１）を入力し（ステップＳ１００１）、時刻ｔにおける隠れ変数φ（ｔ）と観測値ｙ（ｔ）を入力する（ステップＳ１００２）。
【００９４】
続いて、隠れ変数重み計算部２１２は、ｑ（φ（ｔ）｜Φ（ｔ−１），Ｙ（ｔ））を計算し（ステップＳ１００３）、ｐ（φ（ｔ）｜Φ（ｔ−１））を計算する（ステップＳ１００４）。
【００９５】
その後、隠れ変数重み計算部２１２は、ｐ（ｙ（ｔ）｜φ（ｔ））を計算し（ステップＳ１００５）、ｗ（ｔ）を計算する（ステップＳ１００６）。最後に、隠れ変数重み計算部２１２は、ｗ（ｔ）を出力する（ステップＳ１００７）。
【００９６】
続いて、図３のステップＳ４の処理、つまり、ステップＳ３で得られた隠れ変数をもとに隠れ状態量を推定する処理について説明する。図１１は、図３のステップＳ４の詳細を示すフローチャートである。
【００９７】
まず、図１１の処理の概要について説明する。パーティクルフィルタの利用によって、求める事後分布は次の式（２６）のように近似される。
【数１５】
0004972016.tif 000016
【００９８】
今、隠れ変数Φ（ｔ）（ｓ）＝｛φ（ｔ）（ｓ）｝を生成したので、これはｐ（ｘ（ｔ），Φ（ｔ）（ｓ）｜Ｙ（ｔ））を計算することに相当する。計算方法自体は、一般的なKalman filterによる状態推定手法（非特許文献７）と同じである。
【００９９】
具体的には、隠れ状態推定部２２は、まず、時刻ｔ−１における隠れ状態量ｘ（ｔ−１）を入力し（ステップＳ１１０１）、時刻ｔにおける隠れ変数φ（ｔ）と観測量ｙ（ｔ）を入力する（ステップＳ１１０２）。
【０１００】
続いて、隠れ状態推定部２２は、ｓ＝１（ステップＳ１１０３）からｓ＝Ｓまで（ステップＳ１１１５でＮｏ）、ｓの値を「１」ずつ更新しながら（ステップＳ１１１６）、ステップＳ１１０４〜Ｓ１１１４の処理を行う。
【０１０１】
また、隠れ状態推定部２２は、ｉ＝１（ステップＳ１１０４）からｉ＝Ｎｔまで（ステップＳ１１１４でＮｏ）、ｉの値を「１」ずつ更新しながら（ステップＳ１１０６）、ステップＳ１１０５〜Ｓ１１１３の処理を行う。
【０１０２】
隠れ状態推定部２２は、ｃ（ｓ）ｉ（ｔ）＝１か否か判断し（ステップＳ１１０５）、Ｎｏの場合はステップＳ１１０６に進み、Ｙｅｓの場合はｚ（ｓ）ｉ（ｔ）＝ｋで選択されるξ（ｓ）ｋ（ｔ）、ψ（ｓ）ｋ（ｔ）をカルマンフィルタにセットする（ステップＳ１１０７）。
【０１０３】
その後、隠れ状態推定部２２は、ｍ＝１（ステップＳ１１０８）からｍ＝Ｍｔまで（ステップＳ１１１２でＮｏ）、ｍの値を「１」ずつ更新しながら（ステップＳ１１１０）、ステップＳ１１０９、Ｓ１１１１の処理を行う。
【０１０４】
隠れ状態推定部２２は、ｊ（ｓ）ｍ（ｔ）＝ｉか否か判断し（ステップＳ１１０９）、Ｎｏの場合はステップＳ１１１０に進み、Ｙｅｓの場合はｙ（ｓ）ｍ（ｔ）を用いてｘ（ｓ）ｉ（ｔ）をフィルタリングする（ステップＳ１１１１）。
【０１０５】
ステップＳ１１１２でＮｏの場合、隠れ状態推定部２２は、各ｍに対するフィイルタリング結果の平均をｘ（ｓ）ｉ（ｔ）として設定する。
【０１０６】
ステップＳ１１１５でＮｏの場合、隠れ状態推定部２２は、時刻ｔにおける隠れ状態量ｘ（ｔ）を出力する（ステップＳ１１１７）。
【０１０７】
次に、図３のステップＳ５の処理について説明する。この処理は、重みの大きい、すなわち事後分布の高いパーティクル（隠れ変数および隠れ状態量）を選択的にコピーし、重みの小さいパーティクルと入れ換えることで、より確からしいパーティクル集合を保持するために行う。図１２は、図３のステップＳ５の詳細を示すフローチャートである。なお、図１２の処理（計算）は、従来のパーティクルフィルタにおけるリサンプリングと同じである（非特許文献７参照）。
【０１０８】
リサンプリング部２３は、重みｗ（ｔ）、隠れ変数Φ（ｔ）、隠れ状態量Ｘ（ｔ）を入力し（ステップＳ１２０１）、パーティクルの重みｗ（ｔ）を正規化する（ステップＳ１２０２）。
【０１０９】
その後、リサンプリング部２３は、ｓ＝１（ステップＳ１２０３）からｓ＝Ｓまで（ステップＳ１２０８でＮｏ）、ｓの値を「１」ずつ更新しながら（ステップＳ１２０９）、ステップＳ１２０４〜Ｓ１２０７の処理を行う。
【０１１０】
リサンプリング部２３は、ｕ∈[０，１]を生成し（ステップＳ１２０４）、Σｌ＝１ｓ＾−１ｗ（ｔ）（ｌ）（「ｌ＝１」から「ｓ＾−１」までのｗ（ｔ）（ｌ）の総和。以下同様）＜ｕ≦Σｌ＝１ｓ＾ｗ（ｔ）（ｌ）となるｓ＾を見つける（ステップＳ１２０５）。
【０１１１】
リサンプリング部２３は、Φ＾（ｓ）（ｔ）にΦ（ｓ＾）（ｔ）を、Ｘ＾（ｓ）（ｔ）にＸ（ｓ＾）（ｔ）をそれぞれ代入し（ステップＳ１２０６）、ｗ＾（ｓ）（ｔ）に１／ｓを代入する（ステップＳ１２０７）。
【０１１２】
ステップＳ１２０８でＮｏの場合、リサンプリング部２３は、重みｗ＾（ｔ）、隠れ変数Φ＾（ｔ）、隠れ状態量Ｘ＾（ｔ）を出力する（ステップＳ１２１０）。
【０１１３】
続いて、図３のステップＳ６の処理の詳細について説明する。図１３は、図３のステップＳ６の詳細を示すフローチャートである。
【０１１４】
まず、図１３の処理の概要について説明する。図１３の処理は、サンプリングされたξｋ（ｔ），ψｋ（ｔ）を用いて対象のトラッキングを実行し、ハイパーパラメータθのオンライン推定を行うものである。θξｋ（ｔ）の推定は次の式（２７）〜式（２９）の漸化式に従う。
【数１６】
0004972016.tif 000017
【０１１５】
ここで、θξｋ（０）はθξｋ（ｔ）の初期値である。共役性からθξｋ（ｔ）をＰ（ξｋ｜Ｘ（ｔ），θξｋ（０））＝ＮＩＷ（ξｋ；θξｋ（ｔ））となるように定義できて、式（３０）に示すように、ハイパーパラメータθξｋ（ｔ）の推定値はθξｋ（ｔ−１）と推定された隠れ状態量ｘ（ｔ）の尤度を用いて計算できることが示された。
【数１７】
0004972016.tif 000018
【０１１６】
同様にθψｋ（ｔ）に関しても次の式（３１）で推定できる。
【数１８】
0004972016.tif 000019
【０１１７】
これらの更新済みハイパーパラメータθ（ｔ）は次の時刻において式（２２）、式（２３）に用いられる。実際には、ｘ（ｔ）およびｙ（ｔ）は複数の追跡対象ｉおよび観測データｍからなっているので、次の式（３２）、式（３３）に示すように、ｋに関連するデータの集合を用いる。
【数１９】
0004972016.tif 000020
【０１１８】
式（３２）において、「s.t.」は「such that」の略記であり、後ろに続く条件に適合するものだけに適用する、という意味である。つまり、式（３２）は、式（３０）中のｘ（ｔ）はｘ[ｋ]（ｔ）に読み替えて計算する（ただし、ｘ[ｋ]（ｔ）とは、全追跡対象ｉの状態量ｘｉ（ｔ）の中でｚｉ（ｔ）＝ｋであるｉを持つｘｉ（ｔ）の集合である）という意味である。
【０１１９】
具体的には、ハイパーパラメータ更新部２４は、まず、時刻ｔにおける隠れ変数Φ（ｔ）、隠れ状態量ｘ（ｔ）、観測量ｙ（ｔ）を入力する（ステップＳ１３０１）。
【０１２０】
続いて、ハイパーパラメータ更新部２４は、ｓ＝１（ステップＳ１３０２）からｓ＝Ｓまで（ステップＳ１３２０でＮｏ）、ｓの値を「１」ずつ更新しながら（ステップＳ１３２１）、ステップＳ１３０３〜Ｓ１３１９の処理を行う。
【０１２１】
また、ハイパーパラメータ更新部２４は、ｋ＝１（ステップＳ１３０３）からｋ＝Ｋｔまで（ステップＳ１３１８でＮｏ）、ｋの値を「１」ずつ更新しながら（ステップＳ１３１９）、ステップＳ１３０４〜Ｓ１３１７の処理を行う。
【０１２２】
さらに、ハイパーパラメータ更新部２４は、ｘ[ｋ]とｙ[ｋ]に｛ｎｕｌｌ｝を代入した（ステップＳ１３０４）後、ｉ＝１（ステップＳ１３０５）からｉ＝Ｎｔまで（ステップＳ１３１５でＮｏ）、ｉの値を「１」ずつ更新しながら（ステップＳ１３０７）、ステップＳ１３０６〜Ｓ１３１４の処理を行う。
【０１２３】
ハイパーパラメータ更新部２４は、ｃ（ｓ）ｉ（ｔ）＝１か否か判断し（ステップＳ１３０６）、Ｎｏの場合はステップＳ１３０７に進み、Ｙｅｓの場合はｚ（ｓ）ｉ（ｔ）＝ｋか否か判断する（ステップＳ１３０８）。
【０１２４】
ステップＳ１３０８において、ハイパーパラメータ更新部２４は、Ｎｏの場合、ステップＳ１３０７に進み、Ｙｅｓの場合、ｘ[ｋ]に｛ｘ[ｋ]，ｘ（ｓ）ｉ（ｔ）｝を代入し（ステップＳ１３０９）、ｍ＝１（ステップＳ１３１０）からｍ＝Ｍｔまで（ステップＳ１３１４でＮｏ）、ｍの値を「１」ずつ更新しながら（ステップＳ１３１２）、ステップＳ１３１１、Ｓ１３１３の処理を行う。
【０１２５】
ハイパーパラメータ更新部２４は、ｊ（ｓ）ｍ（ｔ）＝ｉか否か判断し（ステップＳ１３１１）、Ｎｏの場合はステップＳ１３１２に進み、Ｙｅｓの場合はｙ[ｋ]に｛ｙ[ｋ]，ｙ（ｓ）ｍ（ｔ）｝を代入する（ステップＳ１３１３）。
【０１２６】
ステップＳ１３１５でＮｏの場合、ハイパーパラメータ更新部２４は、ｘ[ｋ]を用いてθξｋ（ｓ）（ｔ）に代入したθξｋ（ｓ）（ｔ−１）を出力し（ステップＳ１３１６）、ｙ[ｋ]を用いてθψｋ（ｓ）（ｔ）に代入したθψｋ（ｓ）（ｔ−１）を出力する（ステップＳ１３１７）。
【０１２７】
このように、本実施形態の動画像処理装置１によれば、確率的な生成モデルを用いてダイナミクスパターンの数とパラメータを推定すると同時に未知数の移動対象に対する追跡を行うことができる。
【０１２８】
なお、動画像処理装置１の最大の特徴は、複数の対象の追跡と、それぞれの対象のダイナミクスを同時に推定できる点である。具体的には、図２Ｂ、図５、式（４）〜式（１３）などに示すように、ｃ，ｚ，ξ，ψ，ｊという５つの隠れ変数を推定することにより、複数対象追跡とダイナミクス推定を同時に計算できる。
【０１２９】
一方、従来技術では、｛ｃ，ｊ｝だけ（複数の追跡対象があるという情報）、もしくは｛ｚ，ｘ，ｙ｝だけ（複数のダイナミクスがあるという情報）しか推定しないため、本実施形態と比較した場合に推定されない変数による情報が欠落する。
【０１３０】
また、本実施形態の動画像処理装置１によれば、経時とともに追跡対象の数が変わっても、その変化を推定して正しく追跡を行うことができる。これは、式（１４）において時刻ｔでの対象の生成と消滅を表すｃ（ｔ）を推定する処理によって実現できる。
【０１３１】
さらに、本実施形態の動画像処理装置１によれば、推定されるダイナミクスの数を自動的に決定することもできる。これは、式（１５）の計算方法をオンライン計算用に変更した式（１６）でクラスタの数を増減させることによって実現できる。
【０１３２】
また、パーティクルを用いない（例えばカルマンフィルタだけを用いる）従来手法では、各時刻において最善と思われる隠れ変数の推定値を一つだけ求める。一方、本手法ではカルマンフィルタとパーティクルフィルタを組み合わせており、パーティクルごとに個別に隠れ変数の推定を行うため、パーティクルの数だけ異なる推定仮説をもつことになる。
【０１３３】
この違いは、例えば追跡対象が急激に（瞬間的に）運動方向を変化させた場合などに結果として現れる。パーティクルを用いない従来手法の場合には、それまでの動きに基づいて最善と思われる運動方向だけで推定を行うので、急激な変化には対応しにくく、追跡に失敗する可能性が高くなる。一方、パーティクルを用いる本手法の場合には、いくつも仮説を生成しておくので、そのうちの一つでも運動方向の急変を予測していれば引き続き追跡できる可能性が高くなる。このようにいくつも仮説があるとどれを信頼してよいかの判断が必要になるので、本手法では、隠れ変数の重み（≒パーティクルの重み）ｗ（ｔ）で、それぞれの仮説の信頼度を表現する。そして、急激な方向転換のように起こりにくい事象の推定値には小さい重みをつけておけば、通常の（予測通りの）状態での追跡では重要視されない一方で、実際にそのような仮説が起こったときに対応する余地を残すことができる。
【０１３４】
さらに、各隠れ変数による作用と効果について改めて説明する。まず、移動対象の存在を表す隠れ変数（動画像中に存在するか否かを表す隠れ変数）（ｃ（ｔ））を推定することで、移動対象の増減（画面からの消失および画面への登場）を表現できる。この変数が存在しない場合には、追跡対象数は全データ（時間）を通じて一定であると仮定したことになるので、そうでないデータの推定には適用できないモデルとなる。
【０１３５】
また、動作パターンの識別子を表す隠れ変数（ｚ（ｔ））を推定することで、各追跡対象が動作パターンを切り替えながら行動する状況を表現できる。この変数が存在しない場合には、追跡対象は常に一定の動作パターンで行動すると仮定したことになるので、その仮定が正しくない場合には追跡自体が困難となり、例えば表３に示した対数尤度のような追跡精度の指標が劣化する可能性が高くなる。
【０１３６】
さらに、動作パターンのパラメータ（ξ、ψ）を推定することで、このパラメータの初期設定（値）にかかわらず、追跡精度を向上することができる。つまり、このパラメータを推定しない場合には、前記したような複数の動作パターンを許容するモデルを用いたとしても、そのモデルの正確さは動作パターンの初期パラメータで完全に決定されてしまい、その初期設定が良くない場合には追跡精度の向上が図りにくくなる。
【０１３７】
また、観測データと移動対象の対応関係を表す隠れ変数（ｊ（ｔ））を推定することで、複数対象の追跡を行うことができる。この隠れ変数の推定は、複数対象の追跡を行う場合には必須である。この隠れ変数を推定しない場合には、複数の追跡対象が存在する場合に正しい推定を行うことができない、つまり、単一の追跡対象しか追跡できない。
【０１３８】
なお、本発明は前記実施形態に限定されるものではなく、その趣旨を変えない範囲で実施することができる。つまり、ハードウェア、ソフトウェアの具体的な構成について、本発明の主旨を逸脱しない範囲で適宜変更が可能である。
【０１３９】
（実験例）
次に、人工データと実動画データを用いて本手法の性能を確認した実験について説明する。
＜比較手法＞
本手法の比較対象として隠れ変数の数を減じた２モデルを用意し、隠れ変数導入の効果を評価する。これら３モデルの特徴を表１に示した。
【表１】
0004972016.tif 000021
【０１４０】
第１の比較手法は、Saerkkaeらのモデル（非特許文献２参照）に似たモデルである。このモデルはダイナミクスのハイパーパラメータθ（ｔ）の更新（式（３０）、式（３１））を行わず（no）、常にデフォルトの初期値θ（０）を用いる。さらに、ＣＲＰを用いたダイナミクスのクラスタリング（式（１６）、式（１８）、式（１９））も行わない（no）。結果としてすべての追跡対象に対してξ（ｔ）とψ（ｔ）は常に唯一の初期分布からサンプリングされる。
【０１４１】
第２の比較手法は、θ（ｔ）のオンライン推定（式（３０）、式（３１））を行うが（yes）、ｚｉ（ｔ）に対するクラスタリング（式（１６）、式（１８）、式（１９））を行わない（no）。従って、推定されるハイパーパラメータは追跡対象ごとに独立に学習される。このモデルでは、最初のモデルにくらべ複数の追跡対象が異なるダイナミクスを持つことを表現できる点で複雑な（移動対象が多い）場合に対応できるようになっているが、パラメータのパターン化（クラスタリング）が行われない。
【０１４２】
最後に、本手法では、前記した通り、ダイナミクスのクラスタリング（式（１６）、式（１８）、式（１９））を行い（yes）、ハイパーパラメータをオンラインで推定（式（３０）、式（３１））する（yes）。第２の比較手法との相違点は、各追跡対象が共有した運動パターンを時間的に切替えながら軌跡を生成するようにモデル化している点にある。
【０１４３】
〔人工データを用いた実験〕
この実験では、［０：２００］×［０：２００］の仮想２次元空間中の移動質点のトラッキングとクラスタリングをタスクとする。データの例を図１４に示す。図１４は、人工データに対応する図であり、（ａ）が観測データ（「○」で示す）、（ｂ）が正解データである。図１４（ｂ）において、「●」が追跡対象を示し、「○」がノイズを示す。追跡対象（質点）の隠れ状態量は、各対象の位置を表す２次元の実数ベクトルとする。また、観測量も２次元実数ベクトルである。これらはノイズによって劣化した対象位置の情報を表す。なお、これらのベクトルは整数ベクトルでも機能する。各対象は次の式（３４）、式（３５）に示すランダムウォークモデルに従う。
【数２０】
0004972016.tif 000022
【０１４４】
ここで、Ｎ（・）は正規分布を表していて、ｖ（ｔ）およびｗ（ｔ）はそれぞれ正規分布からサンプリングされる。本実験では、各対象は表２に示す４つのダイナミクスパターンをランダムに切り替えながら軌跡と出力を生成するものとした。なお、｛｝Ｔは転置した行列またはベクトルを意味し、diag｛｝は対角行列を意味する。
【表２】
0004972016.tif 000023
【０１４５】
全パターンにおいてｒは零ベクトルである。ＮＩＷの初期値θ（０）は、ξｋ（ｔ）およびψｋ（ｔ）の平均値が次の式（３６）の値となるように設定する。
【数２１】
0004972016.tif 000024
【０１４６】
その他、追跡対象の軌跡推定に必要なパラメータを列挙する。時系列データは３００ステップである。パーティクル数はＳ＝３００、ＣＲＰのconcentration parameterはγ＝２とした。各時刻においてＰｂ＝０．１の確率で新しい追跡対象が発生するものとする。シーン中に存在する追跡対象（ｃｉ（ｔ−１）＝１）は式（３７）で計算される確率（消滅確率）Ｐｄでシーンから消滅するものとする。なお、ｔｎはｊｍ（・）＝ｉとなる観測量ｙｍが少なくとも一つ存在した最後の時刻を表す。また、λ＝０．１とした。
【数２２】
0004972016.tif 000025
【０１４７】
続いて、実験結果を検討する。各時刻におけるデータ対数尤度の平均値を表３に示す。
【表３】
0004972016.tif 000026
【０１４８】
人工データの列の値（実動画像の列の値も同様）は、推定された隠れ変数（すなわちモデルそのもの）から観測データが生成される可能性を表すもので、値が大きい程モデルの適合度が高いことを示す。表３から、Kalman filterのパラメータ推定およびそのクラスタリングが対数尤度を向上させることが確認できる。これは、本手法がオンラインのトラッキングとクラスタリングによってより良いモデルとパラメータを獲得したことを示している。
【０１４９】
図１５は、本手法（π＝０．２）における、正規分布ノイズのパラメータｑの最終時刻における分布図（人工データ実験における速度ノイズの平均値分布図）である。ｑは追跡対象の平均速度バイアスを表している。生成された人工データは表２に示した通り４パターンのｑをもっている。図１５における各プロット点は一つのダイナミクスクラスタのｑの平均値を表している。サイズの数字は各クラスタのデータサイズ（どれだけ多くの対象がそのクラスタに属したか）を表す。大きく４つの主要なクラスタ（ｓｉｚｅ＞１００）が獲得され、それらはデザイン（設定）された「正解（ground truth）」に近い値を得た。すなわち、本手法はダイナミクスパターンのクラスタリングとそれらのパラメータ推定に成功したといえる。
【０１５０】
〔実動画像データによる実験〕
次に、実動画像データを用いた実験について説明する。デジタルカメラで撮影された３２０×２４０ピクセルの動画像を用いて、シーン中の歩行者から抽出される特徴点の追跡とクラスタリングを行う。
【０１５１】
目標となる特徴点は次のようにして抽出した。最初に背景差分（特定の対象を特定の条件で撮影した背景画像と、別途同一条件で撮影した観測画像を比較すること）を行い、二値化によって前景ピクセルを抽出する。続いて、前景部分中から黒色のピクセルだけを選び、それらをmean shiftクラスタリングで少数のクラスタ中心へ量子化する。これらのクラスタ中心座標を特徴点とする。これら特徴点の位置ベクトル（２次元のベクトル）が観測量となる。この位置ベクトルは、一般に実数ベクトルであるが、整数ベクトルとしても機能する。おおよそ一人の歩行者から１〜３個の特徴点が抽出された。
【０１５２】
本実験ではdata association分布（式（２４））に色モデルを設定した（非特許文献１参照）。モデルに用いるのは８ビットのＲＧＢ（Red Green Blue）ヒストグラムである。各特徴点の周辺ピクセルにおけるヒストグラムと、前時刻に得られた追跡対象のヒストグラムの間でBhattacharya係数を計算する。対象ｉのヒストグラムは、ｊｍ（ｔ−１）＝ｉとなった観測点のヒストグラムの平均とする。ｊｍ（ｔ）＝ｉとなる確率は次の式（３８）で計算する。
【数２３】
0004972016.tif 000027
【０１５３】
ここでσは事前に設定した定数パラメータ（σ＝０．１５）で、ｄ（ｍ，ｉ）はｍ番目の観測データとｉ番目の追跡対象のヒストグラム間のBhattacharya coefficientである。
【０１５４】
状態空間モデルは先の人工データ実験と同じものを用いる。ＮＩＷの初期値θ０に関しては、ξｋ（ｔ）およびψｋ（ｔ）の平均値が次の式（３９）の値になるように設定した。
【数２４】
0004972016.tif 000028
【０１５５】
他のパラメータ設定を列挙する。動画像のフレーム数は２００フレームである。この２００フレームは３０ＦＰＳ（Frame Per Second）でキャプチャされた１０００フレームから間引いて抽出した。パーティクルの総数はＳ＝５００個、ＣＲＰのconcentration parameterはγ＝０．１とした。新規追跡対象の生成確率はＰｂ＝０．１とし、シーン中の対象が消滅する確率は式（３７）で計算された確率Ｐｄである（λ＝０．１とする）。
【０１５６】
続いて実験結果を検討する。前記したように、表３の実動画像の列は平均対数尤度を示す。人工データの実験の場合と同様、本手法によって対数尤度が向上することが確認できる。次に、平均速度のバイアスであるｑの分布を図１６と図１７に示す。それぞれπ＝０．１としたときの１０番目、１９９番目（最終）フレームでの分布である。
【０１５７】
実データであるため「正解（ground truth）」は存在しないが、動画像データを観察すると、ほとんどの歩行者は画面内の上下方向にのみ動いていることがわかった。また、ごくたまに店に入るため左右方向に方向転換する。以上のことから、基本的に動きのパターンは上下方向の２種であることが予想される。入力フレーム数が少ない状態（図１６）では、まだそういったトレンドを推測することができないため、各対象ごとに個別のダイナミクスパターンを設定して尤度を上げようとする。しかし、情報が蓄積されると（図１７）、最終的には上下方向２種類のダイナミクスクラスタのみが存在するようになった。
【０１５８】
最後に、トラッキング結果の連続スナップショット（８枚）を図１８に示す。「＃０」などは時系列に沿ったフレーム番号を示す。図１８は、ショッピングモールの吹き抜け部分を上から撮影したものである。各画面の左右両端部分はショッピングモールの２階廊下部分、画面中央部分は１階通路である。１階の各店は２階廊下部分の下に存在し、それらの入り口では光が漏れている。実際には、１階通路において、画面右側を歩く人のほとんどは画面上方向に歩き、画面左側を歩く人のほとんどは画面下方向に歩いている。
【０１５９】
対象の推定位置は矩形で示している。なお、１人の人物に対して、矩形の追跡マークが複数ついている部分もあるが、これは人物を直接追跡しているのではなく、人物とおぼしき部分から黒色領域（髪の部分やズボンの部分など）を個別に抽出し、その中心点を追跡しているからであり、追跡結果は正しいものとなる。
【０１６０】
それぞれの矩形の上部に付した数字はダイナミクスクラスタのインデックス（ｚ）に対応する。上へ動くバイアスを獲得したクラスタは「１」、そして下向きのダイナミクスは「０」に対応している。フレーム１００までは上向きに歩いていた画面右側の二人が、フレーム１２５では店に入るために反転している。本手法では、それを検知して対応する特徴点のダイナミクスパターンを上向き「１」から下向き「０」に切り替えた。つまり、本手法によって複数の移動対象に対する追跡をより高精度で行っていることが確認できる。
【図面の簡単な説明】
【０１６１】
【図１】本実施形態に係る動画像処理装置の構成を模式的に示す機能ブロック図である。
【図２Ａ】複数対象追跡のイメージ図である。
【図２Ｂ】本手法の計算に用いられる生成モデルのグラフィカルモデルである。
【図３】本手法の全体の処理を示すフローチャートである。
【図４】図３のステップＳ３の詳細を示すフローチャートである。
【図５】図４のステップＳ３４の詳細を示すフローチャートである。
【図６】図５のステップＳ３４３の詳細を示すフローチャートである。
【図７】図５のステップＳ３４４の詳細を示すフローチャートである。
【図８】図５のステップＳ３４５の詳細を示すフローチャートであり
【図９】図５のステップＳ３４６の詳細を示すフローチャートであ
【図１０】図１０は、図４のステップＳ３５の詳細を示すフローチャートである。
【図１１】図３のステップＳ４の詳細を示すフローチャートである。
【図１２】図３のステップＳ５の詳細を示すフローチャートである。
【図１３】図３のステップＳ６の詳細を示すフローチャートである。
【図１４】人工データに対応する図であり、（ａ）が観測データ、（ｂ）が正解データである。
【図１５】本手法（π＝０．２）における、正規分布ノイズのパラメータｑの最終時刻における分布図である。
【図１６】実動画像を用いた実験における平均速度のバイアスであるｑの分布図である（１０番目のフレーム）。
【図１７】実動画像を用いた実験における平均速度のバイアスであるｑの分布図である（１９９番目のフレーム）。
【図１８】実動画像を用いた実験におけるトラッキング結果の連続スナップショット（８枚）である。
【符号の説明】
【０１６２】
１      動画像処理装置
２      演算手段
３      記憶手段
４      入力手段
５      出力手段
６      バスライン
２１    隠れ変数集合推定部
２２    隠れ状態推定部
２３    リサンプリング部
２４    ハイパーパラメータ部
２１１  隠れ変数決定部
２１２  隠れ変数重み計算部
２１１１  移動対象数決定部
２１１２  動作パターン番号決定部
２１１３  動作パターンパラメータ決定部
２１１４  移動対象対応関係決定部
【図１】
0004972016.tif 000029
【図２Ａ】
0004972016.tif 000030
【図２Ｂ】
0004972016.tif 000031
【図３】
0004972016.tif 000032
【図４】
0004972016.tif 000033
【図５】
0004972016.tif 000034
【図６】
0004972016.tif 000035
【図７】
0004972016.tif 000036
【図８】
0004972016.tif 000037
【図９】
0004972016.tif 000038
【図１０】
0004972016.tif 000039
【図１１】
0004972016.tif 000040
【図１２】
0004972016.tif 000041
【図１３】
0004972016.tif 000042
【図１４】
0004972016.tif 000043
【図１５】
0004972016.tif 000044
【図１６】
0004972016.tif 000045
【図１７】
0004972016.tif 000046
【図１８】
0004972016.tif 000047