JP4568223B2 - 推定装置 - Google Patents

Description

本発明は、観測データから隠れ状態を推定する推定装置に関する。

従来より、運転中のドライバの頭部姿勢など、推定対象（ここではドライバ）を拘束することなく推定対象の状態（ここでは頭部姿勢）を検出したい場合に、カメラで撮影された画像を観測データとして用いることが行われている。

しかし、この場合、観測データである画像から直接取得されるのは、画素毎の輝度や色情報に過ぎず、従って、画像から頭部姿勢を直接検出することはできない。
このような観測データ（画像）から直接検出することができない推定対象の状態（例えば頭部姿勢）のことを隠れ状態と呼び、この推定対象の隠れ状態の事後確率分布を観測データから算出する手法として、時系列ベイズ推定が知られている。

時系列ベイズ推定では、種々のモデル変数（隠れ状態変数を含む）の分布を求める際に、複雑な積分計算を必要とするため、一般的には、解析的に解くことができないが、モデル変数の確率分布がガウス分布に従い、モデルを規定するシステム方程式が線形性を有することを仮定することによって、推定対象の隠れ状態の事後確率分布の算出に伴う積分計算を解析的に行うことを可能としたカルマンフィルタを用いて、高速に時系列ベイズ推定を実現する手法が知られている（例えば、非特許文献１参照。）。

しかし、ドライバの頭部姿勢を推定することを考えた場合、実際の環境下では、頭部の複雑な動きや、直射日光や西日、街灯などによる照明条件の急激な変化、ドライバの手や帽子などの遮蔽物の存在などに基づく複雑な外乱（非ガウス性雑音）が観測データに加わる。つまり、モデル変数の確率分布がガウス分布に従わなかったり、システム方程式が線形から外れたりするため、カルマンフィルタにて、安定した推定精度を得ることは一般に困難である。

また、観測データとして画像を用いた場合、一般に、扱うモデル変数の次元が高次元となる場合が多く、上述した積分計算のために、非常に膨大な処理負荷が発生し、実時間での実現が非常に困難であるという問題がある。

これに対して、非ガウス性雑音が加わった観測データから隠れ状態の時系列推定や予測を近似的な時系列ベイズ推定によって行うパーティクルフィルタが知られている。
このパーティクルフィルタは、有限のパーティクルによって各時刻での隠れ状態の事前確率分布及び事後確率分布を離散的に表現し、それを用いて時系列推定や予測を行うものである。

但し、パーティクルフィルタによる事後確率分布の近似精度を高めるためには、パーティクル数を増加させる必要があり、一方、実時間処理のためにはパーティクル数を少なく抑える必要がある。つまり、パーティクルフィルタを用いた場合、精度と処理時間との間には、トレードオフの関係がある。

そこで、状態変数を解析的に積分可能なものとそうでないものに分け，解析的に積分不可能なものに対してのみパーティクルフィルタを用いることにより、パーティクルフィルタで推定する次元を小さくでき、比較的少ないパーティクル数によって推定精度の向上と計算量の抑制を両立させるRao-Blackwellisationという手法が知られている（例えば、非特許文献２参照。）。

画像処理の分野でもこのRao-Blackwellizationの適用例があり、一般に推定対象の隠れ状態と観測データの中間に推定対象の隠れ状態と観測データの双方と因果関係を持つ中間の隠れ状態（例えば，複数の顔の特徴点の画像面上の座標）を設けることで、線形ガウス過程を仮定することができる推定対象の隠れ状態と中間の隠れ状態とから構成される確率システム（上位層）と線形ガウス過程を仮定することができない中間の隠れ状態と観測可能な状態とから構成される確率システム（下位層）とに分け、上位層における時系列推定では例えばカルマンフィルタを用い、下位層における時系列推定ではパーティクルフィルタを用いることにより、パーティクルフィルタで推定する次元を小さくして、推定精度の向上と処理時間の抑制を両立させる手法が知られている（例えば、非特許文献３参照。）。

なお、Rao-Blackwellisationの適用には、状態変数の一部が解析的に積分可能である必要があるが、これは線形ガウス過程モデルに対するカルマンフィルタに限定されるものではない。（例えば混合正規分布モデル，隠れマルコフモデル，ディリクレ過程モデルなどでもよい。）
A. Blake, R. Curwen and A. Ziseerman ： "A frame work for spatio-temporal control in the tracking of visual contours", International Journal of Computer Vision, 11, 2, pp.127-145 (1993)． G. Casella and C. Robert ： "Rao-blackwellisation of sampling schmes", Biometrika, 83, 1, pp.81-94 (1996)． A. Doucet, S. Godsil and C. Andrieu ： "On sequential monte carlo sampling methods for bayesian filtering", Statistics and Computiong, 10, 3, pp.197-208 (2000)．

しかし、非特許文献３に記載の手法では、下位層においてパーティクルフィルタで推定した中間の隠れ状態の推定結果（中間の隠れ状態の事後確率分布）を推定対象の隠れ状態に対する観測データとして用いて上位層で推定対象の隠れ状態の推定を行っているため、下位層で中間の隠れ状態の推定精度が低下すると、上位層での推定対象の隠れ状態の推定精度も低下しまい、しかも、下位層で中間の隠れ状態の推定精度が一旦低下すると、これを正常な状態に復元することが一般に困難であるという問題があった。

即ち、例えば、下位層の推定に用いるパーティクルフィルタが、ドライバの顔のある特徴点が画像（観測データ）上に投影された座標（中間の隠れ状態）を時系列推定するように構成されていると仮定した場合、この特徴点と類似した部分（偽特徴点という）が観測データ中に存在すると、この偽特徴点は特徴点と共に大きな尤度を持つことになる。このため、特徴点が不規則で急激な動きをする等した場合に、中間の隠れ状態の事前確率分布及び事後確率分布を離散近似するパーティクル群が特徴点から外れて偽特徴点に捕らわれてしまうことがある。このようにパーティクル群が偽特徴点に一旦捕らわれてしまうと、以後、パーティクルフィルタでは、その誤った特徴点の座標の推定結果（中間の隠れ状態の事後確率分布）から生成した予測（中間の隠れ状態の事前確率分布）に従うパーティクル群を用いて、次回の観測及び推定を行うことになるため、パーティクル群が偽特徴点を追跡する異常な状態から抜け出すことが困難になってしまうのである。

そこで、本発明は、非ガウス性雑音が存在する環境でも、精度良くしかもリアルタイムに時系列ベイズ推定の良い近似計算を行うことが可能な推定装置を提供することを目的とする。

ここで、図１７，図１８は、以下で使用する用語の定義の理解を容易にするための説明図である。
上記目的を達成するためになされた本発明の推定装置では、推定対象の隠れ状態予測手段が、当該手段での時刻ｔ−１の推定情報と、推定対象の隠れ状態の時間変化に伴う動きのモデルを表す推定対象の隠れ状態の状態遷移構造とに基づいて、時刻ｔの推定対象の隠れ状態の状態予測情報を生成する。

すると、推定対象の隠れ状態の状態遷移構造による中間の隠れ状態予測手段が、推定対象の隠れ状態予測手段によって生成された時刻ｔの推定対象の隠れ状態予測情報と，推定対象の隠れ状態と前記中間の隠れ状態との間の因果関係とに基づいて，時刻ｔの推定対象の隠れ状態の状態遷移構造による中間の隠れ状態予測情報を生成する。

そして、中間の隠れ状態の尤度観測手段が、後述する混合手段によって生成された時刻ｔの中間の隠れ状態予測情報と、時刻ｔの観測可能な状態の観測データとに基づいて，時刻ｔの中間の隠れ状態の尤度を生成する。

すると、中間の隠れ状態推定手段が、後述する混合手段によって生成された時刻ｔの中間の隠れ状態予測情報と，中間の隠れ状態の尤度観測手段によって生成された時刻ｔの中間の隠れ状態の尤度とに基づいて時刻ｔの中間の隠れ状態推定情報を生成する。

そして、推定対象の隠れ状態の尤度観測手段が、中間の隠れ状態推定手段によって生成された時刻ｔの中間の隠れ状態推定情報と，推定対象の隠れ状態予測手段によって生成された時刻ｔの推定対象の隠れ状態予測情報とに基づいて，時刻ｔの推定対象の隠れ状態の尤度を生成する。

すると、推定対象の隠れ状態推定手段が、推定対象の隠れ状態推定手段によって生成された時刻ｔの推定対象の隠れ状態予測情報と，推定対象の隠れ状態の尤度観測手段によって生成された時刻ｔの推定対象の隠れ状態の尤度とに基づいて時刻ｔの推定対象の隠れ状態推定情報を生成する。

なお、推定対象の隠れ状態の状態遷移構造による中間の隠れ状態予測手段が、時刻ｔの推定対象の隠れ状態予測情報と、推定対象の隠れ状態と中間の隠れ状態との因果関係とに基づいて中間の隠れ状態を予測した、時刻ｔの推定対象の隠れ状態の状態遷移構造による中間の隠れ状態予測情報を生成し、中間の隠れ状態の状態遷移構造による中間の隠れ状態予測手段が、時刻ｔ−１の中間の隠れ状態推定情報と、中間の隠れ状態の状態遷移構造とに基づいて中間の隠れ状態を予測した、時刻ｔの中間の隠れ状態の状態遷移構造による中間の隠れ状態予測情報を生成し、混合手段が、時刻ｔの推定対象の隠れ状態の状態遷移構造による中間の隠れ状態予測情報と、時刻ｔの中間の隠れ状態の状態遷移構造による中間の隠れ状態予測情報とを、予め設定された混合比に従って混合することで時刻ｔの中間の隠れ状態予測情報を生成する。

つまり、中間の隠れ状態推定手段が使用する時刻ｔの中間の隠れ状態予測情報には、その中間の隠れ状態推定手段での時刻ｔ−１での推定結果（時刻ｔ−１での中間の隠れ状態推定情報）だけでなく、推定対象の隠れ状態推定手段での時刻ｔ−１での推定結果（時刻ｔ−１での推定対象の隠れ状態推定情報）も反映されている。

そして、時刻ｔ−１の中間の隠れ状態推定情報から中間の隠れ状態の状態遷移構造によって導かれる時刻ｔの中間の隠れ状態の状態遷移構造による中間の隠れ状態予測情報と、時刻ｔ−１の推定対象の隠れ状態推定情報から推定対象の隠れ状態の状態遷移構造と、推定対象の隠れ状態と中間の隠れ状態との因果関係とからによって導かれる時刻ｔの推定対象の隠れ状態の状態遷移構造による中間の隠れ状態予測情報とは異なっているため、時刻ｔの中間の隠れ状態の状態遷移構造による中間の隠れ状態予測情報と時刻ｔの推定対象の隠れ状態の状態遷移構造による中間の隠れ状態予測情報とが同時に大きく精度を低下させてしまう従来手法の問題を解決することができる。

従って、本発明の推定装置によれば、何等かの原因（観測データが画像、中間の隠れ状態が画面上に投影された顔の特徴点の座標である場合は、顔の特徴点の複雑な動きや照明の変動、遮蔽物による顔の特徴点の隠れなど）で、時刻ｔの推定対象の隠れ状態の状態遷移構造による中間の隠れ状態予測情報及び時刻ｔの中間の隠れ状態の状態遷移構造による中間の隠れ状態予測情報のうち一方の精度が大きく低下したとしても、両者を混合することで生成される時刻ｔの中間の隠れ状態予測情報は、その影響が緩和されたものとなり、実環境で発生する状況の変化や大きな外乱の影響を除去することができ、頑強で正確な推定を実現することができる。

ところで、中間の隠れ状態推定手段は、時刻ｔの中間の隠れ状態推定情報と、その時刻ｔの中間の隠れ状態推定情報の生成に使用した時刻ｔの定対象の隠れ状態の状態遷移構造による中間の隠れ状態予測情報及び時刻ｔの中間の隠れ状態の状態遷移構造による中間の隠れ状態予測情報とに基づき、時刻ｔの中間の隠れ状態推定情報を最も良く近似する時刻ｔの推定対象の隠れ状態の状態遷移構造による中間の隠れ状態予測情報と時刻ｔの中間の隠れ状態の状態遷移構造による中間の隠れ状態予測情報との混合比を求める混合比算出手段を備え、この混合比算出手段にて算出された混合比を、時刻ｔ＋１の中間の隠れ状態予測情報を生成する際に使用するように構成されていることが望ましい。

なお、時刻ｔの中間の隠れ状態推定情報を最も良く近似する混合比を求めた場合、時刻ｔの推定対象の隠れ状態の状態遷移構造による中間の隠れ状態予測情報及び時刻ｔの中間の隠れ状態の状態遷移構造による中間の隠れ状態予測情報のうち予測精度の高い方がより高い比率で混合されるような混合比となる。

つまり、本発明の推定装置によれば、時刻ｔの推定対象の隠れ状態の状態遷移構造による中間の隠れ状態予測情報と時刻ｔの中間の隠れ状態の状態遷移構造による中間の隠れ状態予測情報との混合比を、予測精度の良好な側がより大きな割合で混合されるように、状況に応じて動的に変化させることになるため、より頑強で正確な推定を実現することができる。

なお、混合比算出手段では、混合比の算出に、オンライン化したＥＭアルゴリズムを用いることで，統計的に望ましい推定値をオンラインで得ることができ，変動するデータに追従することができるため，本推定装置に更なる頑健性を与える。

また、例えば、推定対象の隠れ状態推定手段は、推定対象の隠れ状態の確率分布が線形性，ガウス性を有し、中間の隠れ状態推定手段は、中間の隠れ状態の確率分布が非線形性，非ガウス性を有するものとして構成されていてもよい。

そして、この場合、中間の隠れ状態推定手段は、中間の隠れ状態推定情報の生成にパーティクルフィルタを用いること、また、推定対象の隠れ状態推定手段は、推定対象の隠れ状態の推定にカルマンフィルタを用いることが望ましい。

なお、観測データが顔面を撮影した画像であり、推定対象の隠れ状態が頭部姿勢を表すアフィンパラメータであり、中間の隠れ状態が顔面の特徴を表す顔の特徴点の画像上での座標である場合には、本発明の推定装置を、ドライバの頭部姿勢を推定して、居眠りやわき見運転を防止する運転支援システム等を構成する際に好適に用いることができる。

以下に本発明の実施形態を図面と共に説明する。
図１は、本発明を適用した画像処理装置１の全体構成を示すブロック図である。なお、画像処理装置１は、車両に搭載され、ドライバの頭部を撮影した画像から、ドライバのわき見や居眠り等を検出するための情報である頭部姿勢を時系列的に求める処理を実行する装置である。

図１に示すように、画像処理装置１は、ドライバシートの前面上部又は下部（例えば、メータ内）に設置され、ドライバの頭部を含む画像を取得するカメラ２と、カメラ２が取得した画像を取り込むビデオキャプチャーボード３と、初期化指令が入力された場合に、ビデオキャプチャーボード３を介して取り込んだ入力画像ｚ_t （ｔは時刻を表す）から、ドライバの顔の特徴を表す複数（Ｎ個）の特徴点の座標（以下、顔特徴点という）ｘ＿｛ｎ｝（ｎ＝１，２，…Ｎ）を抽出する初期特徴点抽出部４と、初期特徴点抽出部４が抽出した顔特徴点ｘ＿｛ｎ｝のそれぞれについて、その顔特徴点ｘ＿｛ｎ｝を中心とする予め設定された大きさの画像をテンプレートｔｐⁿ として切り出して記憶するテンプレート生成・記憶部５と、初期特徴点抽出部４が抽出した顔特徴点ｘ＿｛ｎ｝間の幾何学的拘束条件を記述した初期顔モデル（観測行列Ｃａ）を生成して記憶する観測行列生成・記憶部６と、入力画像ｚ_t ，テンプレートｔｐⁿ ，観測行列Ｃａに基づいて、時系列ベイズ推定を行うことにより、頭部姿勢を表すアフィンパラメータａ_t を推定する推定部８と、推定部８が推定したアフィンパラメータａ_t に従って、頭部姿勢θ_x，θ_y，θ_zを算出する頭部姿勢算出部９から構成される。

このうち、初期特徴点抽出部４では、図２（ａ）に示すように、顔特徴点として、目尻，目頭，口元、鼻（鼻の穴，鼻の中心，鼻の左右端）などを抽出する。但し、本実施形態では、左右目尻、左右目頭，左右口元，鼻の中心の７（＝Ｎ）個の点を抽出するものとする。また、初期化指令は、ドライバが正面を向き、顔特徴点の隠れのない画像をカメラ２で撮影した時に、ドライバの操作によって入力されるものとする。

また、観測行列生成・記憶部６では、図２（ｂ）に示すように、顔特徴点ｘ＿｛１｝〜ｘ＿｛Ｎ｝（線図の各頂点）が一つの平面上にあるとした平面顔モデルを仮定して、顔特徴点ｘ＿｛ｎ｝間の位置関係（相対位置，距離）を表す観測行列Ｃ_a を生成する。なお、平面顔モデルを仮定するのは、計算を簡単にするためである。具体的には、図３に示すように、画面上の座標を表す水平方向の座標軸をｕ軸、垂直方向の座標軸をｖ軸、３次元空間の水平方向の座標軸をｘ軸、垂直方向の座標軸をｙ軸、光軸方向（奥行き方向）の座標軸をｚ軸とし、初期特徴点抽出部４で抽出されたｎ番目の顔特徴点の座標を（ｃｕ_n，ｃｖ_n）で表すものとして、観測行列Ｃ_a は、（１）式にて表される。

テンプレート生成・記憶部５が生成するテンプレートｔｐⁿ は、具体的には、図２（ａ）中の点線白枠にて示すようなものとなる。但し、図では、左口元、鼻の中心、左目尻についてのみ示す。

推定部８で推定されるアフィンパラメータａ_t は、実世界（三次元のｘｙｚ空間）での回転による座標変換を、射影面である画像面（二次元のｕｖ平面）でのアフィン変換による座標変換で表した場合に得られる変換行列の要素ａ₁〜ａ₆を縦に並べて得られるベクトルである。

即ち、ある特徴点の三次元座標を（ｘ，ｙ，ｚ）^T 、この特徴点をｘ軸を中心に角度θｘ，ｙ軸を中心に角度θｙ，ｚ軸を中心に角度θｚだけ回転させた移動後の座標を（ｘ’，ｙ’，ｚ’）^T とすると、移動（回転）前の三次元座標と移動後の三次元座標の関係は（２）式で表される。

これら三次元座標（ｘ，ｙ，ｚ）^T ，（ｘ’，ｙ’，ｚ’）^T を画像面上に射影した二次元座標、すなわち画像上での特徴点の座標を（ｕ，ｖ）^T ，（ｕ’，ｖ’）^T とした場合、移動前の二次元座標と移動後の二次元座標の関係は、アフィンパラメータａ_t を用いて（３）式で表される。

そして、頭部姿勢算出部９では、頭部姿勢情報として、頭部の３軸角度θｘ，θｙ，θｚを（４）〜（６）式を用いて算出する。

つまり、（ｘ’，ｙ’，ｚ’）^T の画像面への正射影と、（ｕ’，ｖ’）^T とが一致すると仮定すると、（２）（３）式から（７）〜（１０）式が導かれ、更に、これら（７）〜（１０）式から（４）〜（６）式が導かれる。

このように、時刻ｔにおける頭部姿勢を表すアフィンパラメータａ_t が得られれば、平面顔モデルの初期位置（即ち、観測行列Ｃ_a ）から、時刻ｔの全ての顔特徴点の位置Ｐ_t （＝Ｃ_aａ_t ）を予測することができる。

また、時刻ｔにおけるアフィンパラメータａ_t は、顔特徴点ｘ＿｛１｝〜ｘ＿｛Ｎ｝の初期位置ｘ＿｛１｝(０)〜ｘ＿｛Ｎ｝(０)と時刻ｔでの位置ｘ＿｛１｝(ｔ)〜ｘ＿｛Ｎ｝(ｔ)との対応関係を３組以上特定できれば求めることができる。従って、アフィンパラメータａ_t が求められると、その求めたアフィンパラメータａ_t から、逆に、時刻ｔにおける画像に映っていない隠れた顔特徴点の位置を予測することも可能となる。

以下では、本発明の主要部である推定部８について詳述する。
ここで、図４は、推定部８が扱う状態空間モデルを表すグラフである。
図４に示すように、推定部８は、時刻ｔにおける顔モデルの姿勢を表すアフィンパラメータａ_t を隠れ状態として求める上位層と、時刻ｔにおける顔特徴点群の位置ｘ(ｔ)(＝（ｘ＿｛１｝(ｔ), ｘ＿｛２｝(ｔ), …，ｘ＿｛Ｎ｝(ｔ)）^T)を隠れ状態として求める顔特徴点ｘ＿｛ｎ｝毎に設けられたＮ個の下位層（図では一つだけ示す）とを備えており、時刻ｔまでに入力された入力画像系列ｚ_1:t から、これら隠れ状態ａ_t ，ｘ(ｔ) を、時系列ベイズ推定により推定するものである。

なお、時系列ベイズ推定は、図５に示すように、状態変数をすべて確率分布として扱い、隠れ状態に関する時刻ｔ−１における推定結果（推定分布）から時刻ｔにおける予測分布を求めると共に、時刻ｔにおける観測データから、検出すべき隠れ状態であるもっともらしさ（尤度）を求め、これら予測分布と尤度を考慮して時刻ｔにおける推定結果（推定分布）を得るという手順を繰り返すことで、逐次的に隠れ状態を推定するものである。

つまり、推定部８では、入力画像系列（観測データ）ｚ_1:t に基づいて、アフィンパラメータ（頭部姿勢の隠れ状態）ａ_t を推定した事後確率分布（推定分布）ｐ（ａ_t ｜ｚ_1:t ）を求めることになり、これを数式で記述すると（１１）（１２）式で表される。

ここで、ｐ（ａ_t ｜ｚ_1:t-1 ）はアフィンパラメータａ_t の事前確率分布（予測分布）、ｐ（ｘ(ｔ)｜ａ_t ，ｚ_1:t-1 ）は顔特徴点群の位置ｘ(ｔ)の事前確率分布（予測分布）、ｐ（ｚ_t ｜ａ_t ，ｚ_1:t-1 ），ｐ（ｚ_t ｜ｘ(ｔ)）は尤度を表す。

そして、（１２）式の部分を下位層がパーティクルフィルタを用いて推定し、（１１）式の部分を上位層がカルマンフィルタを用いて推定する。
なお、あるパラメータｆの確率分布がガウス分布（正規分布）に従う場合、その確率分布は、平均をμ，分散をΣ（標準偏差σ）とすると、（１３）式で表すことができる。つまりこの場合、パラメータの確率分布の計算は、実際には、その平均μと分散Σを求めれば十分である。

次に、図６は推定部８の具体的な構成を表すブロック図である。
図６に示すように、推定部８は、顔特徴点ｘ＿｛ｎ｝毎に設けられたＮ個の追跡器ＢＫ_n からなり、それぞれが一つの顔特徴点ｘ＿｛ｎ｝を追跡して、時刻ｔにおけるその顔特徴点ｘ＿｛ｎ｝の位置ｘ＿｛ｎ｝(ｔ)を推定した確率分布をガウス近似した確率分布を規定する特徴点推定分布パラメータ（平均、標準偏差）を生成する顔特徴点追跡部１０と、各追跡器ＢＫ_n にて算出された特徴点推定分布パラメータに基づき、カルマンフィルタを用いてアフィンパラメータａ_t 及びアフィンパラメータの予測値ａ_t+1，その分散Ｖ_t+1 （以下、予測値等という）を算出するアフィンパラメータ算出部３０と、アフィンパラメータ算出部３０にて算出されるアフィンパラメータの予測値等を保持する予測値保持部４０と、予測値保持部４０にて保持された時刻ｔ−１において算出されたアフィンパラメータの予測値等に基づいて、時刻ｔの顔特徴点群の位置ｘ(ｔ) (＝（ｘ＿｛１｝(ｔ), ｘ＿｛２｝(ｔ), …，ｘ＿｛Ｎ｝(ｔ)）^T)を予測した確率分布ｐ（ｘ(ｔ)｜ａ_t ）を規定する上位予測分布パラメータ（平均値，分散）を算出して、追跡器ＢＫ_n のそれぞれに供給する予測分布パラメータ算出部４１とを備えている。

つまり、アフィンパラメータ算出部３０が上述の上位層、顔特徴点追跡部１０を構成する各追跡器ＢＫ_n が上述の下位層に相当する。
なお、アフィンパラメータ算出部３０で使用するカルマンフィルタでは、アフィンパラメータａ_t の状態遷移過程として、アフィンパラメータａ_t の動きモデルが等速直線運動であること、システムノイズが標準正規分布（ガウス分布でμ＝０，Σ＝１）の分散を線形変換したガウス分布に従うことを表す（１４）式に示す関係式を仮定し、また、下位層にて推定される顔特徴点群の位置ｘ(ｔ) の期待値ベクトルｆ_t の観測過程として、（１５）式に示す線形の関係式を仮定している。ここで、このような観測モデルを仮定するのは、カルマンフィルタの適用が可能となるためには、システムが線形で表現されることと、状態及び尤度の確率分布がガウス分布で表されることとが必要なためである。そして、このガウス分布で近似された顔特徴点群の位置の期待値ベクトルｆ_t はアフィンパラメータａ_t を推定する際の観測情報に相当する。

但し、Ａ_a は動きモデル（等速直線運動）、Ｂ_a はシステムノイズの分散の大きさ（予測からの外れ具合）、Ｃ_a は観測行列（顔特徴点座標とアフィンパラメータとの間の幾何関係）、Ｄ_a は観測ノイズの分散の大きさ（観測値に加わっているノイズの大きさ）を表す行列である。また、ｗ_a,t ，ｖ_a,t は、標準正規分布に従うノイズを表すベクトルである。

ここで図７は、アフィンパラメータ算出部３０の具体的な構成を表すブロック図である。
図７に示すように、アフィンパラメータ算出部３０は、顔特徴点追跡部１０から、顔特徴点ｘ＿｛１｝〜ｘ＿｛Ｎ｝毎に供給される特徴点推定分布パラメータ（平均、標準偏差）に基づいて、時刻ｔにおける各顔特徴点の位置ｘ＿｛１｝(ｔ)〜ｘ＿｛Ｎ｝(ｔ)の確率分布の，期待値（平均）を配列してなる顔特徴点群の位置の期待値ベクトルｆ_t 、標準偏差を配列してなる雑音行列Ｄ_a,t を生成する算出器３１と、予測値保持部４０を介して供給される時刻ｔ−１におけるアフィンパラメータの推定値に基づく予測値等（以下単にアフィンパラメータ予測値等という）、及び算出器３１にて生成された雑音行列Ｄ_a,t 、観測行列生成・記憶部６に記憶された観測行列Ｃ_a に基づいてカルマンゲインＫ_t を算出するカルマンゲイン算出部３２と、カルマンゲインＫ_t ，アフィンパラメータ予測値等，顔特徴点群の位置の期待値ベクトルｆ_t ，観測行列Ｃ_a に基づいて、事後確率分布を表現するアフィンパラメータ及びその分散の推定値（以下単にアフィンパラメータ推定値等という）ａ_t ，Ｖ_t を推定するアフィンパラメータ推定部３３と、アフィンパラメータ推定部３３にて推定されたアフィンパラメータ推定値等ａ_t ，Ｖ_t に基づいて、事前確率分布を表現する次時刻ｔ＋１におけるアフィンパラメータ予測値等を算出するアフィンパラメータ予測部３４とからなる。

なお、算出器３１が生成する顔特徴点群の位置の期待値ベクトルｆ_t は（１６）式、雑音行列Ｄ_a,t は（１７）式で表される。また、ｄｉａｇは、各成分を対角成分として持つ対角行列を表す。

また、カルマンゲイン算出部３２ではカルマンゲインＫ_t を（１８）式により、アフィンパラメータ推定部３３ではアフィンパラメータの推定値ａ_t を（１９）式、アフィンパラメータの分散の推定値Ｖ_t を（２０）式により算出する。

更に、アフィンパラメータ予測部３４では、アフィンパラメータの予測値を（２１）式、アフィンパラメータの分散の予測値を（２２）式により算出する。但し、各式において、ハットが付されたパラメータａ_t ，Ｖ_t は推定値を表し、バーが付されたパラメータａ_t ，Ｖ_t は予測値を表す（以下同様）。

また、（１３）式を用いることにより、アフィンパラメータａ_t の事後確率分布ｐ（ａ_t ｜ｚ_1:t ）は（２３）式、事前確率分布ｐ（ａ_t+1 ｜ｚ_1:t ）は（２４）式で表される。

ところで、アフィンパラメータの推定値ａ_t を算出する（１９）式は、顔特徴点群の位置ｘ(ｔ)の観測値（ｆ_t ）と予測値（Ｃ_a ａ_t ）との差をカルマンゲインＫ_t で乗じた補正値を、予測値に加算したものを推定値ａ_t とすることに相当する。

このとき、カルマンゲインＫ_t は、顔特徴点群の位置ｘ(ｔ)の観測値と予測値との混合比を制御するものと見なすことができ、予測値の分散が大きいほど（又は観測値の分散が小さいほど）、カルマンゲインＫ_t は大きな値となり、逆に、予測値の分散が小さいほど（又は観測値の分散が大きいほど）カルマンゲインＫ_t は小さな値となる。

なお、予測値や観測値の分散が大きいということは、予測値や観測値の信頼度が低い（想定したモデルから外れ易い，観測値に大きなノイズが加わり易い）ことを意味し、予測値や観測値の分散が小さいということは、予測値や観測値の信頼度が高い（想定したモデルに近い動きをし易い，観測値にノイズが加わり難い）ことを意味する。

つまり、アフィンパラメータ算出部３０では、予測値より観測値の信頼度が高い時には、カルマンゲインＫ_t が大きくなり、観測値を重視した推定値を生成し、逆に、観測値より予測値の信頼度が高い時には、カルマンゲインＫ_t が小さくなり、予測値を重視した推定値を生成することになる。

次に、上位予測分布パラメータ算出部４１では、顔特徴点の位置ｘ(ｔ)の上位予測分布ｐ（ｘ(ｔ)｜ａ_t ）は（２５）式で与えられるため、この式に従って、上位予測分布パラメータである平均μ及び分散Σを算出する。

次に、顔特徴点追跡部１０を構成する追跡器ＢＫ_n について説明する。
なお、追跡器ＢＫ_n は、パーティクルフィルタによって一つの顔特徴点ｘ_n を追跡するものである。ここで、パーティクルフィルタの動作概要を、図８に示す説明図を参照して説明する。

図８に示すように、パーティクルフィルタでは、推定したい対象の状態（顔特徴点の位置）の実現値（画像上の座標）をパーティクルで表し、予測，観測，リサンプルという三つの過程を繰り返し実行することで、時系列的にパーティクルの分布を求める。なお、パーティクルフィルタでは、カルマンフィルタとは異なり、ガウス分布に限らず任意の確率分布を対象とすることができる。

まず、予測過程では、推定したい対象の動きを考慮して、状態空間（ここでは画像上）での各パーティクルの状態（ここでは画像上での位置）を遷移させ、推定したい対象の動きに加わるノイズを考慮してパーティクルを確率的に散らすことで，推定したい対象が存在しそうな位置にパーティクル群を配置する。これにより、推定したい対象の状態を予測した確率分布が、パーティクル群の位置と数により離散的かつ近似的に表現されることになる。

次に、観測過程では、各パーティクルの状態が推定したい対象の状態であることの尤度（ここでは顔特徴点を表すテンプレートとの正規化相関値）を、尤度に応じてパーティクルの重みを算出する。これにより、推定したい対象の状態の確率分布が、重みのあるパーティクルによって表現されることになる。

また、リサンプル過程では、重みの小さいパーティクルを消去すると共に重み大きいパーティクルを重みのない複数のパーティクルに増殖させることにより、推定したい対象の状態の確率分布を、重みのあるパーティクルによる表現から、重みのないパーティクルにより表現に確率的に変換する。

そして、このリサンプル（リサンプル過程で生成）されたパーティクルを用いて上述の予測過程を行うことにより、パーティクルで表された推定したい対象の状態の確率分布が時系列的に求められることになる。

ここで図９は、顔特徴点ｘ＿｛ｎ｝を追跡する追跡器ＢＫ_n の具体的な構成を示すブロック図である。
図９に示すように、追跡器ＢＫ_n は、パーティクルのリサンプル値（後述する）に基づいて顔特徴点を予測した確率分布（以下、下位予測分布という）に従ったパーティクルを生成する下位予測分布生成部１１と、上位予測分布パラメータ算出部４１にて生成された顔特徴点についての上位予測分布パラメータ（顔特徴点座標の平均，分散）、及び下位予測分布生成部１１にて生成されたパーティクルに基づいて、上位予測分布と下位予測分布とを予め設定された混合比α_a,t ，α_z,t （０≦α_a,t ≦１，０≦α_z,t ≦１，α_a,t ＋α_z,t ＝１）にて混合してなる混合分布π（ｘ(ｔ)｜α_t ）に従ったパーティクルを生成する混合分布生成部１４とを備えている。

また、追跡器ＢＫ_n は、混合分布生成部１４が生成するパーティクルの重みを、観測データｚ_t ，テンプレートｔｐⁿ を用いて算出する観測部２２と、混合分布生成部１４が生成するパーティクルの位置、及び観測部２２が算出するパーティクルの重みに基づいて、顔特徴点の位置を推定した確率分布（以下、推定分布という）を規定する特徴点推定分布パラメータ（平均，標準偏差）を算出する推定分布パラメータ算出部２５と、観測部２２が算出するパーティクルの重みに従って、混合分布生成部１４が生成したパーティクルをリサンプルするリサンプル部２６と、リサンプル部２６がリサンプルしたパーティクル、上位予測分布パラメータ算出部４１にて算出された上位予測分布パラメータ、下位予測分布生成部１１にて生成されたパーティクルに基づいて、混合分布生成部１４が次回使用する混合比α_a,t ，α_z,t を生成する混合比設定部１８とを備えている。

このうち、予測分布生成部１１は、リサンプル部２６にてリサンプルされたＮｐ個のパーティクルを保持するリサンプルパーティクル保持部１２と、リサンプルパーティクル保持部１２に保持された前回の推定結果（推定分布）に従ったパーティクルを遷移／拡散させることで、下位予測分布に従ったパーティクルを生成するパーティクル遷移／拡散部１３とからなる。

なお、パーティクル遷移／拡散部１３では、各パーティクルに（２６）式を作用させる。

つまり、行列Ａ_f は、顔特徴点のダイナミクスとして等速直線運動を仮定したことを意味するものであり、図１０に示すように、パーティクル群の位置を前々回（時刻ｔ−２）の推定位置から前回（時刻ｔ−１）の推定位置までの間の移動量だけ遷移（シフト）させることに相当し、また、行列Ｂ_f は、移動量にランダム（白色ガウス）ノイズを付加することを意味するものであり、パーティクル群の位置を拡散させることに相当する。

混合分布生成部１４は、上位予測分布パラメータμ_n,t ，σ_n,t に従って、上位予測分布ｐ（ｘ(ｔ)｜ａ_t ）に従ったα_a,t ×Ｎｐ個のパーティクルを発生させるパーティクル生成部１５と、予測分布生成部１１にて生成された下位予測分布ｐ（ｘ(ｔ)｜ａ_t-1 ，ｚ_1:t-1 ）に従ったＮｐ個のパーティクルの中から、α_z,t ×Ｎｐ個のパーティクルを確率的に選択するパーティクル選択部１６と、パーティクル生成部１５にて生成されたパーティクル、及びパーティクル選択部１６にて選択されたパーティクルを混合して、上位予測分布と下位予測分布との混合分布π（ｘ(ｔ)｜α_t ）≒ｐ（ｘ(ｔ)｜ａ_t ，ｚ_1:t-1 ）に従ったＮｐ個のパーティクルを出力する混合部１７とからなる。

なお、パーティクル生成部１５は、具体的には、平均及び分散（又は標準偏差）を入力することで、その平均及び分散（又は標準偏差）で特定される正規分布に従った乱数を発生する周知の乱数発生器を用いればよい。

観測部２２は、混合分布生成部１４にて生成されたパーティクルのそれぞれについて、入力画像ｚ_t から、パーティクルが表す位置を中心に、テンプレートｔｐⁿ と同じ大きさの画像（以下、切出画像という）Ｓ_t を切り出す画像切出部２３と、画像切出部２３が切り出した切出画像Ｓ_t のそれぞれについて、切出画像Ｓ_t とテンプレートｔｐⁿ との類似度を表す正規化相関値を求め、その求めた正規化相関値に従って各パーティクルの重みｗ_t を設定する重み算出部２４とからなる。

なお、重み算出部２４では、（２７）式により重みｗ_t を決定し、このパーティクルの重みと位置から求められる確率分布が今回（時刻ｔ）の推定分布（事後確率分布）ｐ（ｘ(ｔ)｜ａ_t ，ｚ_1:t ）となる。

なお、図８及び図１１中では、パーティクルの重みｗ_t を濃淡で示す。
つまり、真の観測過程（（１２）式におけるｐ（ｚ_t ｜ｘ(ｔ)））は未知であるため，パーティクルの重みｗ_t の計算には何らかのモデルが必要となり、ここでは、テンプレートｔｐⁿ と切出画像Ｓ_t との正規化相関を用いてパーティクルの重みｗ_t を決定している。このため、重み算出部２４におけるパーティクルの重みｗ_t の計算が、陰に観測過程を規定していることになる。

推定分布パラメータ算出部２５は、パーティクルの重みｗ_t を考慮して、混合分布生成部１４が生成したパーティクルの位置の平均（期待値）、及び標準偏差（又は分散）を算出する。このように、重みを持つパーティクルにて表現される確率分布を、その平均と標準偏差で表現されるガウス分布で表すということは、この非ガウス性を有する確率分布を、ガウス性を有する確率分布で近似することを意味する。そして、アフィンパラメータ算出部３０では、上述したように、この平均と標準偏差が示すガウス分布を用いて顔特徴点の尤度を求め、アフィンパラメータの推定に用いる。

リサンプル部２６は、位置と重みにより確率分布を表現するパーティクルを、位置と数により確率分布を表現するパーティクルに確率的に変換する。具体的には、パーティクルの持つ重みの大きさに応じて、重みの小さいパーティクルを消去し、重みの大きいパーティクルを増殖させる。

これは、図８及び図１２に示すように、重みの大きいパーティクルが存在する位置、即ち、推定分布が高い値を持つ位置に、パーティクルを確率的に寄せることに相当する。このようにしてリサンプルされたパーティクルは、顔特徴点の事後確率分布を表すものとして、予測分布生成部１１と混合比設定部１８とに供給される。

次に、混合比設定部１８は、予測分布生成部１１にて生成された下位予測分布に従うパーティクルに基づき、下位予測分布を規定する下位予測分布パラメータ（平均，分散）を算出する下位予測分布パラメータ算出部１９と、上位予測分布パラメータ算出部４１にて算出された上位予測分布パラメータ、下位予測分布パラメータ算出部１９にて算出された下位予測分布パラメータ、リサンプル部２６にて生成された推定分布に従うパーティクルに基づき、次回の推定で使用する混合比α_a,t+1 ，α_z,t+1 を算出する混合比算出部２０と、混合比算出部２０にて算出された混合比α_a,t+1 ，α_z,t+1 を保持して、前回の推定結果から算出された混合比を混合分布生成部１４に供給する混合比保持部２１とからなる。

このうち、混合比算出部２０では、推定分布（リサンプルされたパーティクル）と、その生成に用いられた二つの予測分布とに基づき、二つの予測分布から生成される混合分布が推定分布と最も近くなるように求めた混合比を、次回の推定で使用する混合比として算出する。

具体的には、二つの分布の距離を測る指標として、（２８）式に示す周知のカルバックライブラーダイバージェンスを使用し、これを最小化する混合比αを求める。この問題は、結局、（２９）式に示す対数尤度Ｌ（α）を最大にすること、即ち、（３０）式に示すように、対数尤度Ｌ（α）の微分値がゼロとなるαを求めるという問題に帰着する。

実際には、過去の不正確な推定による影響を徐々に忘却するための忘却係数をλｓ（０≦λｓ≦１）として、この忘却係数λｓを組み込んだ（３１）式に示す対数尤度を用いる。

しかし、（３１）式の最大化は、非線形な形となり計算が困難であるため、ここでは（３１）式のαに関する最大化と等価となる（３６）式のαに関する最大化をオンラインＥＭアルゴリズムを用いて行う。この（３６）式は（３１）式で示した対数尤度の、指標変数ｍ（注目するパーティクルが上位と下位のどちらの予測分布からサンプルされたかを示す隠れ変数である）の予測事後分布についての期待値（期待対数尤度）である。

このＥＭアルゴリズムでは、まず、現在設定されている（今回の推定で使用した）混合比α_t ＝｛α_a,t ，α_z,t ｝を用いて、上位予測分布及び下位予測分布について（３７）（３８）式で表される値（図１３参照）を、リサンプルされたパーティクルのそれぞれについて算出する（Ｅ−ＳＴＥＰ）。

次に、Ｅ−ＳＴＥＰでの算出値と、（４０）〜（４３）式を用いることにより混合比α_a ，α_z を算出し、この混合比α_a ，α_z を次回の混合分布を生成する際に用いる混合比α_a,t+1 ，α_z,t+1 とする（Ｍ−ＳＴＥＰ）。

なお、（４０）〜（４３）式は、（３６）式の右辺の確率分布を表す二つの項のうち、前側（式中左側）の項に、Ｅ−ＳＴＥＰでの算出結果を代入することで期待対数尤度の一部を定数化し、その一部が定数化された期待対数尤度を用いて（３９）式を、α_a,t+1 ＋α_z,t+1 ＝１という制約下で解くことで求められる。

このようにして求められた混合比α_a,t+1 ，α_z,t+1 は、下位予測分布（顔特徴点の推定分布）に基づく予測精度が高い（即ち、顔特徴点の追跡が上手くいっている）ときは、下位予測分布（特徴点による予測）を重視したものとなり、顔特徴点の推定分布に基づく予測精度が低い（即ち、顔特徴点の追跡が上手くいっていない）ときは、上位予測分布（顔モデルによる予測）を重視したものとなるようにされている。

以上説明したように、本実施形態の画像処理装置１では、画面上における顔特徴点ｘ＿｛ｎ｝の追跡（位置の推定）をパーティクルフィルタを用いて行う下位層と、下位層で得られた顔特徴点群の位置ｘ(ｔ)の推定結果（尤度）から顔モデルへのアフィン変換過程をカルマンフィルタを用いて推定する上位層とを有する階層的なモデルを用いて、入力画像からアフィンパラメータａ_tを推定するようにされている。

また、本実施形態の画像処理装置１では、下位層での顔特徴点ｘ＿｛ｎ｝の推定に使用する予測分布として、上位層でのアフィンパラメータの推定結果ａ_t に基づく上位予測分布と下位層での顔特徴点ｘ＿｛ｎ｝の推定結果に基づく下位予測分布とを混合してなる混合分布を用いており、しかも、下位予測分布の予測精度が高い時には下位予測分布を重視し、下位予測分布の予測精度が低い時には上位予測分布を重視した混合分布が生成されるように混合比を動的に変化させ、精度の低下した階層を互いに補い合うようにされている。

このように、本実施形態の画像処理装置１によれば、非ガウス性ノイズが重畳しやすい顔特徴点ｘ＿｛ｎ｝の追跡（位置の推定）にのみパーティクルフィルタを用いているため、パーティクルフィルタで使用するパーティクル数を抑制することができると共に、実環境において発生する状況の変化や大きなノイズの影響を最小限に抑えることができる。

つまり、下位層では、個々の顔特徴点ｘ＿｛ｎ｝の位置を高い精度で推定できるため、精度の高い下位予測分布を得ることができるが、その反面、ノイズ（手やめがねによる顔特徴点の隠れ／顔向きの大きな変化／照明条件の変化など）の影響を受けた場合には、推定結果（ひいては下位予測分布）が大きく外れやすいという特徴がある。一方、上位層では、平面顔モデルを仮定したことに基づく誤差が存在し、また、一つ顔特徴点での精度低下が、全ての顔特徴点の上位予測分布に影響するため、アフィンパラメータの推定結果から個々の顔特徴点を高精度に予測することは困難であるが、その反面、顔特徴点群の位置ｘ(ｔ) を全体として推定するため、顔特徴点の一部がノイズの影響を大きく受けていても、その結果として得られる上位予測分布が大きく外れたものとなってしまうことがないという特徴がある。このように相補的な性質を持つ二つの予測分布を混合した混合分布は、二つの予測分布の長所を併せ持つことになる。

従って、手や飲み物等により顔特徴点の一部が隠されたり、運転行動における高速な頭部姿勢の変化（例えば、正面方向からルームミラー方向に顔を向け、再び正面に向き直るなど）により、顔特徴点が予め仮定した動きモデル（等速直線運動）とは大きく異なる動きをする等して、下位予測分布の推定精度が大きく低下したとしても、混合分布を使用する下位層での顔特徴点の推定精度が大きく低下することがない。その結果、本実施形態の
画像処理装置１によれば、上位層でのアフィンパラメータの推定、ひいてはアフィンパラメータに基づくドライバの頭部姿勢の算出を高い信頼性，頑健性をもって行うことができる。

また、本実施形態の画像処理装置１では、オンライン化したＥＭアルゴリズムを用いて混合比の算出を行うため、動的に変化する実際の環境への対処が可能となっている。
ここで、画像処理装置１の動作をシミュレーションによって確認した結果を示す。

なお、テンプレートｔｐⁿ には、実際の車内でドライバーが正面を向いている時に撮影した画像（正面顔画像）から、人手によって切り出した両眼・鼻孔画像を用いた。
また、シミュレーション画像は、図３に示すように、各顔部品（両眼，鼻）は、単一平面上にあるものとし、この平面がＹ軸を中心とした回転角（ヨー角）がθ_yだけ回転するものとして、次のように作成した。

即ち、頭部を左右に振っている人の頭部姿勢のヨー角を、実際の車内で連続的に計測し、その計測結果に基づいて、図１４に示すように、各時刻毎に回転後の座標を透視変換により画像面上に射影することで得られる画像面上の両眼中心および鼻の位置に、正面顔画像から切り出したテンプレートを配置したものをシミュレーション画像とした。つまり、シミュレーション画像では、各時刻のヨー角に応じて画像上での両眼中心・鼻の位置は変化するが、各顔部品は常に正面向きになっている。

更に、車内における遮蔽（手，飲み物などによる顔特徴点の隠れ）を模擬するために、シミュレーション画像には、観測ノイズの大きな大ノイズ領域Ｒを設けた。但し、この大ノイズ領域Ｒは、頭部姿勢が右向きのときに、右眼尻がちょうど隠れるような大きさ，位置に設定した。

また各顔特徴点を追跡する追跡器ＢＫ_n が用いるパーティクルの個数Ｎｐは、実時間処理が十分に実現可能な５０個とし、混合比α_a,t ，α_z,t を算出する際に用いるオンラインＥＭアルゴリズムの忘却係数λｓは、０．５とした。更に、（２６）式におけるｂ_u，ｂ_vは、ｂ_u ＝ｂ_v ＝７、（２７）式におけるＤ_f は、Ｄ_f ＝０．１５とした。

上述のようにして作成したシミュレーション画像列に対して、本実施形態の画像処理装置１を適用し、推定されたアフィンパラメータａ_t から復元した頭部姿勢のヨー角を図１５に示す。この図１５では、推定結果のヨー角が、正面方向を表す０°を中心にプラス方向，マイナス方向の両方に変化しており、頭部を左右に振る様子がよく推定されていることがわかる。

また、図１６（ａ）には、シミュレーションを５０回試行し、その結果に基づいて頭部姿勢の平均推定誤差（平均二乗誤差）と平均誤差標準偏差を求めた結果を示す。なお、ここでは比較のために、上位層の混合比α_a,t を１（即ち、上位予測分布のみ使用）及び０（即ち、下位予測分布のみを使用）に固定して、同様のシミュレーションを行った場合の結果も示す。

また、図１６（ｂ）には、各顔特徴点について算出された混合比α_a,t を示す。
但し、ここでは、混合比α_a,t の最大値を０．８、最小値を０．２とした。図１６（ｂ）中のグレー領域は、顔特徴点が図１４の大ノイズ領域に隠れていることを示す。

図示されているように、右眼，右眼尻が大ノイズ領域に隠れる期間では、上位層の混合比α_a,t が大きくなって、上位予測分布（顔モデルに基づく予測）を重視した混合分布が生成されるように変化している。つまり、画像処理装置１では、実環境における複雑な外乱に対応して、その時の状況に即した混合比を算出（ひいては混合分布を生成）し、外乱に対して頑強で正確な顔特徴点の自動追跡を実現していることがわかる。

なお、本実施形態において、アフィンパラメータ予測部３４が「推定対象の隠れ状態予測手段」、上位予測分布パラメータ算出部４１が「推定対象の隠れ状態の状態遷移構造による中間の隠れ状態予測手段」、観測部２２が「中間の隠れ状態の尤度観測手段」、推定分布パラメータ算出部２５が「中間の隠れ状態推定手段」、算出器３１が「推定対象の隠れ状態の尤度観測手段」、カルマンゲイン算出部３２及びアフィンパラメータ推定部３２が「推定対象の隠れ状態推定手段」、下位予測分布生成部１１が「中間の隠れ状態の状態遷移構造による中間の隠れ状態予測手段」、混合分布生成部１４が「混合手段」、混合比設定部１８が「混合比算出手段」に相当する。
［他の実施形態］
以上、本発明の一実施形態について説明したが、本発明は上記実施形態に限定されるものではなく、本発明の要旨を逸脱しない範囲において、様々な態様にて実施することが可能である。

例えば、上記実施形態では、画像処理装置１を入力画像に基づいてドライバの頭部状態を推定するために使用しているが、これに限定されるものではなく、入力画像以外の観測データを用いて、ドライバの頭部状態以外の隠れ状態を推定するために使用してもよい。

なお、応用例としては顔の追跡に限らず、階層表現を持ったいかなる剛体の追跡でも良い。さらに一般的には、推定分布や予測分布が解析的に積分可能となる変数の階層関係が保証されていれば、追跡対象を剛体に限る必要もなく、またデータを画像データに限る必要もない。

実施形態の推定装置の全体構成を示すブロック図。 顔特徴点、平面顔モデルを示す説明図。 推定装置に適用する座標系を示す説明図。 推定部が扱う状態空間モデルを表すグラフ。 時系列ベイズ推定の動作概略を示す説明図。 推定部の構成を示すブロック図。 アフィンパラメータ算出部の構成を示すブロック図。 パーティクルフィルタの動作概要を示す説明図。 追跡器の構成を示すブロック図。 パーティクルフィルタの動作（遷移／拡散）の具体例を示す説明図。 パーティクルフィルタの動作（観測）の具体例を示す説明図。 パーティクルフィルタの動作（リサンプル）の具体例を示す説明図。 混合比の算出に必要なパラメータの意味を示す説明図。 シミュレーションに使用した画像の構成を示す説明図。 シミュレーションによるヨー角の算出結果を示すグラフ。 シミュレーションの結果に基づいて算出した頭部姿勢の誤差を示す表、及びＥＭアルゴリズムによる顔特徴点毎の混合比の算出結果を示すグラフ。 状態遷移構造に関わる用語を理解するための説明図。 状態間でやりとりされる情報、及び各手段との関係を示す説明図。

符号の説明

１…画像処理装置、２…カメラ、３…ビデオキャプチャーボード、４…初期特徴点抽出部、５…テンプレート生成・記憶部、６…観測行列生成・記憶部、８…推定部、９…頭部姿勢算出部、１０…顔特徴点追跡部、１１…下位予測分布生成部、１２…リサンプルパーティクル保持部、１３…パーティクル遷移／拡散部、１４…混合分布生成部、１５…パーティクル生成部、１６…パーティクル選択部、１７…混合部、１８…混合比設定部、１９…下位予測分布パラメータ算出部、２０…混合比算出部、２１…混合比保持部、２２…観測部、２３…画像切出部、２４…重み算出部、２５…推定分布パラメータ算出部、２６…リサンプル部、３０…アフィンパラメータ算出部、３１…算出器、３２…カルマンゲイン算出部、３３…アフィンパラメータ推定部、３４…アフィンパラメータ予測部、４０…予測値保持部、４１…上位予測分布パラメータ算出部、ＢＫ_n …追跡器。

Claims (6)

1. 観測可能な状態から推定対象の隠れ状態を時系列に推定する推定装置であって、
   前記推定対象の隠れ状態から観測可能な状態が観測される過程が階層構造を有し、
   前記階層構造は、前記推定対象の隠れ状態と前記観測可能な状態の間に、前記推定対象の隠れ状態と前記観測可能な状態の両者と因果関係を持つ中間の隠れ状態を持ち，前記推定対象の隠れ状態が時間変化に伴って状態遷移する構造を持つ階層構造であって、
   前時刻（以下、時刻ｔ−１とする）の推定対象の隠れ状態推定情報と推定対象の隠れ状態の状態遷移構造とに基づいて、現時刻ｔ（以下、時刻ｔとする）の推定対象の隠れ状態予測情報を生成する推定対象の隠れ状態予測手段と、
   前記時刻ｔの推定対象の隠れ状態予測情報と、前記推定対象の隠れ状態と前記中間の隠れ状態との間の因果関係とに基づいて、時刻ｔの推定対象の隠れ状態の状態遷移構造による中間の隠れ状態予測情報を生成する推定対象の隠れ状態の状態遷移構造による中間の隠れ状態予測手段と、
   時刻ｔの観測可能な状態の観測データと、後述する混合手段によって生成された時刻ｔの中間の隠れ状態予測情報とに基づいて、時刻ｔの中間の隠れ状態の尤度を生成する中間の隠れ状態の尤度観測手段と、
   前記時刻ｔの中間の隠れ状態予測情報と、前記時刻ｔの中間の隠れ状態の尤度と、に基づいて時刻ｔの中間の隠れ状態推定情報を生成する中間の隠れ状態推定手段と、
   前記時刻ｔの中間の隠れ状態推定情報と、前記時刻ｔの推定対象の隠れ状態予測情報とに基づいて、時刻ｔの推定対象の隠れ状態の尤度を生成する推定対象の隠れ状態の尤度観測手段と、
   前記時刻ｔの推定対象の隠れ状態予測情報と、前記時刻ｔの推定対象の隠れ状態の尤度と、に基づいて時刻ｔの推定対象の隠れ状態推定情報を生成する推定対象の隠れ状態推定手段と、
   を備え、
   前記階層構造に、前記中間の隠れ状態にも時間変化に伴って状態遷移する構造と、
   時刻ｔ−１の中間の隠れ状態推定情報と中間の隠れ状態の状態遷移構造とに基づいて、時刻ｔの中間の隠れ状態の状態遷移構造による中間の隠れ状態予測情報を生成するところの、中間の隠れ状態の状態遷移構造による中間の隠れ状態予測手段と、
   前記時刻ｔの推定対象の隠れ状態の状態遷移構造による中間の隠れ状態予測情報と、前記時刻ｔの中間の隠れ状態の状態遷移構造による中間の隠れ状態予測情報を予め設定された混合比に従って混合することで前述の時刻ｔの中間の隠れ状態予測情報を生成する混合手段と、
   を設けたことを特徴とする推定装置。
2. 前記中間の隠れ状態推定手段は、
   前記時刻ｔの中間の隠れ状態推定情報と、該時刻ｔの中間の隠れ状態推定情報の生成に使用した前記時刻ｔの推定対象の隠れ状態の状態遷移構造による中間の隠れ状態予測情報及び前記時刻ｔの中間の隠れ状態の状態遷移構造による中間の隠れ状態予測情報とに基づき、前記時刻ｔの中間の隠れ状態推定情報を最もよく近似する前記時刻ｔの推定対象の隠れ状態の状態遷移構造による中間の隠れ状態予測情報と前記時刻ｔの中間の隠れ状態の状態遷移構造による中間の隠れ状態予測情報との混合比を求める混合比算出手段を備え、該混合比算出手段にて算出された混合比を、次時刻（以下、時刻ｔ＋１とする）の中間の隠れ状態予測情報を生成する際に使用することを特徴とする請求項１に記載の推定装置。
3. 前記混合比算出手段は、前記混合比の算出に、オンライン化したＥＭアルゴリズムを用いることを特徴とする請求項２に記載の推定装置。
4. 前記中間の隠れ状態推定手段は、前記時刻ｔの中間の隠れ状態推定情報の生成に、パーティクルフィルタを用いることを特徴とする請求項１〜３のいずれかに記載の推定装置。
5. 前記推定対象の隠れ状態推定手段は、前記時刻ｔの推定対象の隠れ状態推定情報の生成に、カルマンフィルタを用いることを特徴とする請求項１〜４のいずれかに記載の推定装置。
6. 前記観測可能な状態は、顔面を撮影した画像であり、前記推定対象の隠れ状態は、頭部姿勢を表すアフィンパラメータであり、前記中間の隠れ状態は、顔面の特徴的な部位の画像面上での位置を表す複数の特徴点の画像面上での座標であることを特徴とする請求項１〜５のいずれかに記載の推定装置。