(57)【要約】【課題】源信号の状態が時間変化する場合であっても、信号特徴の状態の特定を可能とする推定システムを提供する。【解決手段】１以上の素子と、プロセッサと、メモリと、源信号の状態と状態間の時間的遷移を表す源信号状態モデルを物体種別ごとに格納したデータベースと、前記素子から出力された入力信号と、前記源信号状態モデルとを比較して源信号の状態を推定して状態推定結果を出力する状態推定部と、前記状態推定結果を提示する結果出力部と、を有する。【選択図】図２
2016065764.tif 000002
【特許請求の範囲】
【請求項１】
１以上の素子と、
プロセッサと、
メモリと、
源信号の状態と状態間の時間的遷移を表す源信号状態モデルを物体種別ごとに格納したデータベースと、
前記素子から出力された入力信号と、前記源信号状態モデルとを比較して源信号の状態を推定して状態推定結果を出力する状態推定部と、
前記状態推定結果を提示する結果出力部と、
を有することを特徴とする推定システム。
【請求項２】
請求項１に記載の推定システムであって、
前記素子から出力された入力信号に基づいて信号源の方向を推定し、方向推定結果として出力する方向推定部と、
複数の時刻における前記方向推定結果から、同時に存在しうる各物体の軌跡からなる集合を仮説として生成し、当該仮説の集合を算出する軌跡推定部と、をさらに有し、
前記状態推定部は、
前記仮説の集合の内の各軌跡に対して源信号の状態を推定し、前記推定された状態推定結果の尤もらしさに応じて前記仮説の集合の内の軌跡を選択し、当該選択された軌跡を出力することを特徴とする推定システム。
【請求項３】
請求項２に記載の推定システムであって、
前記状態推定部は、
前記選択された軌跡と、前記状態推定結果に応じて物体数を決定することを特徴とする推定システム。
【請求項４】
請求項２に記載の推定システムであって、
前記結果出力部が、前記選択された軌跡を出力することを特徴とする推定システム。
【請求項５】
請求項２に記載の推定システムであって、
前記状態推定部は、
前記選択された軌跡に基づいて各時刻に物体間が近接しているか否かを判定し、前記物体間が近接していない時刻の入力信号に応じて前記源信号状態モデルを更新することを特徴とする推定システム。
【請求項６】
１以上の素子と、プロセッサと、メモリとを備えた推定装置で、信号の状態を推定する推定方法であって、
前記推定装置が、前記素子から出力された入力信号を受け付ける第１のステップと、
前記推定装置が、源信号の状態と状態間の時間的遷移を表す源信号状態モデルを物体種別ごとに格納したデータベースから読み込んだ前記源信号状態モデルと、前記素子から受け付けた入力信号とを比較して、源信号の状態を推定して状態推定結果を出力する第２のステップと、
前記推定装置が、前記状態推定結果を提示する第３のステップと、
を含むことを特徴とする推定方法。
【請求項７】
請求項６に記載の推定方法であって、
前記推定装置が、前記素子から受け付けた前記入力信号に基づいて信号源の方向を推定し、方向推定結果として出力する第４のステップと、
前記推定装置が、複数の時刻における前記方向推定結果から、同時に存在しうる各物体の軌跡からなる集合を仮説として生成し、当該仮説の集合を算出する第５のステップと、をさらに含み、
前記第２のステップは、
前記仮説の集合の内の各軌跡に対して源信号の状態を推定し、前記推定された状態推定結果の尤もらしさに応じて前記仮説の集合の内の軌跡を選択し、当該選択された軌跡を出力することを特徴とする推定方法。
【請求項８】
請求項７に記載の推定方法であって、
前記第２のステップは、
前記選択された軌跡と、前記状態推定結果に応じて物体数を決定することを特徴とする推定方法。
【請求項９】
請求項７に記載の推定方法であって、
前記第３のステップが、
前記選択された軌跡を出力することを特徴とする推定方法。
【請求項１０】
請求項７に記載の推定方法であって、
前記第２のステップは、
前記選択された軌跡に基づいて各時刻に物体間が近接しているか否かを判定し、前記物体間が近接していない時刻の入力信号に応じて前記源信号状態モデルを更新することを特徴とする推定方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、信号の状態を推定する推定システムおよび方法に関する。
【背景技術】
【０００２】
信号の状態を推定する技術として、特許文献１が知られている。この特許文献１には、「予め用意した移動体の観測時刻毎の観測値を元に、該移動体の位置、速度、加速度及び角速度の少なくともひとつを含む移動経路の諸元を推定する方法であって、移動体が一種類の運動により移動した区間の複数からなる移動経路の諸元を推定する方法において、移動経路の諸元の要素と、区間毎の基本状態変数との対応関係を示すδ、及び、移動経路の諸元の次数Ｐを、その区間の運動の種類に応じて区間毎に求める段階と、前記観測値を元に予め求めた区間毎の基本状態遷移行列Φと、前記対応関係δとに応じた区間毎の状態遷移行列Φ￣を求める段階と、全区間における前記次数Ｐ及び前記状態遷移行列Φ￣に応じた計画係数データ行列Ｆを求める段階と、前記計画係数データ行列Ｆ、及び、前記観測値に基づく観測データ列Ｌに応じて、移動経路の諸元の推定値β＾を求める推定段階とを含むことを特徴とする」と記載されている。
【先行技術文献】
【特許文献】
【０００３】
【特許文献１】特開２００７−０３３２５２号公報
【発明の概要】
【発明が解決しようとする課題】
【０００４】
水中パッシブソーナーや、空中音響監視システムなど、物体から発せられる信号（源信号）を観測し、その物体の源信号の信号特徴の状態を推定するシステムは広く存在する。ここでの「信号特徴」とは、源信号が音響信号である場合は、たとえば源信号にフーリエ変換を施した後の周波数領域信号を表すベクトル、もしくは、源信号に離散コサイン変換やウォルシュ・アダマール変換などの基底変換を施した後の信号を表すベクトル、もしくは、源信号のケプストラムや源信号のウェーブレット変換後の信号を表すベクトルなどを意味する。
【０００５】
また、物体は同時に信号特徴の「状態」と呼ぶラベルを一つ持っており、その信号特徴の「状態」が対応する確率分布や確率過程が信号特徴を生成するものとする。確率分布や確率過程はたとえば多次元混合ガウス分布（ＧＭＭ）やガウス過程（ＧＰ）である。
【０００６】
上記特許文献１には、複数の運動状態がある場合に物体を追跡する方法が記載されている。しかし、源信号の信号特徴の状態が時間変化する場合に、その信号特徴の状態を特定する方法は開示されていない。このように上記従来のシステムでは、源信号の信号特徴の状態が時間変化する場合に信号特徴の状態を特定できない、という問題があった。
【０００７】
そこで、本発明は、上記問題点に鑑みて成されたもので、源信号の信号特徴の状態が時間変化する場合であっても、信号特徴の状態の特定を可能とする推定システムを提供する。
【課題を解決するための手段】
【０００８】
本発明は、１以上の素子と、プロセッサと、メモリと、源信号の状態と状態間の時間的遷移を表す源信号状態モデルを物体種別ごとに格納したデータベースと、前記素子から出力された入力信号と、前記源信号状態モデルとを比較して源信号の状態を推定して状態推定結果を出力する状態推定部と、前記状態推定結果を提示する結果出力部と、を有する。
【発明の効果】
【０００９】
本発明によれば、源信号の信号特徴の状態が時間変化する場合であっても、信号特徴の状態の特定が可能である。上記した以外の課題、構成及び効果は、以下の実施形態の説明により明らかにされる。
【図面の簡単な説明】
【００１０】
【図１】本発明の第１の実施例を示し、推定システムのハードウェア構成の一例を示すブロック図である。
【図２】本発明の第１の実施例を示し、推定システムの機能の一例を示すブロック図である。
【図３】本発明の第１の実施例を示し、推定システムで行われる処理の一例を示すフローチャートである。
【図４】本発明の第１の実施例を示し、結果出力部がユーザインタフェース部に推定結果を提示した例を示す画面イメージである。
【図５Ａ】本発明の第１の実施例を示し、源信号状態モデルデータベースのデータ構造の一例を示す図である。
【図５Ｂ】本発明の第１の実施例を示し、源信号状態モデルデータベースのデータ構造の一例を示す図である。
【図６】本発明の第１の実施例を示し、状態推定処理の一例を示すフローチャートである。
【図７】本発明の第２の実施例を示し、推定システムのハードウェア構成の一例を示すブロック図である。
【図８】本発明の第２の実施例を示し、推定システムの機能の一例を示すブロック図である。
【図９】本発明の第２の実施例を示し、推定システムで行われる処理の一例を示すフローチャートである。
【図１０】本発明の第２の実施例を示し、結果出力部がユーザインタフェース部に推定結果を提示した例を示す画面イメージである。
【図１１】本発明の第２の実施例を示し、状態推定処理の一例を示すフローチャートである。
【図１２】本発明の第３の実施例を示し、推定システムの機能の一例を示すブロック図である。
【図１３】本発明の第３の実施例を示し、推定システムで行われる処理の一例を示すフローチャートである。
【図１４】本発明の第３の実施例を示し、結果出力部がユーザインタフェース部に推定結果を提示した例を示す画面イメージである。
【発明を実施するための形態】
【００１１】
以下、本発明の実施形態を添付図面に基づいて説明する。
【実施例１】
【００１２】
本実施例１では、源信号（物体から発せられる信号）の信号特徴の状態が時間変化する場合であっても、信号特徴の状態の特定が可能な、推定システム１００の例について説明する。本実施例１の推定システムは例えば、水中パッシブソーナーや空中音響監視システムである。なお、前記課題で述べたように、源信号は、物体から発せられる信号であり、信号特徴は、源信号が音響信号である場合は、たとえば源信号にフーリエ変換を施した後の周波数領域信号を表すベクトル等を示す。
【００１３】
図１は、本実施例の推定システム１００のハードウェア構成の一例を示すブロック図である。
【００１４】
本実施例１の推定システム１００は、ハイドロホンやマイクロホンなどの素子１０１−１と、ＡＤ変換器１０２と、中央演算装置１０３と、ユーザインタフェース部１０４と、記憶媒体１０５と、揮発性メモリ１０６から構成されている。なお、揮発性メモリ１０６に代わって不揮発性メモリを用いてもよい。
【００１５】
素子１０１＿１は、水中あるいは空中の音波をアナログ電圧値に変換し、ＡＤ変換器１０２に送る。ＡＤ変換器１０２は、アナログ電圧値をデジタル信号に変換し、中央演算装置１０３に送る。中央演算装置１０３は、入力されたデジタル信号を入力信号として揮発性メモリ１０６（または記憶媒体１０５）に格納する。
【００１６】
揮発性メモリ１０６（または記憶媒体１０５）に格納した入力信号に基づいて、中央演算装置１０３が所定の推定処理を行う。この処理は、記憶媒体１０５に格納されたプログラムを、揮発性メモリ１０６に読み込んで中央演算装置１０３が実行する。
【００１７】
推定処理で出力される結果は、ユーザインタフェース部１０４が提示する。ユーザインタフェース部１０４は、モニタであってもよく、ネットワーク経由で接続された別のＰＣを経由したモニタ、マウス、キーボードなどであってもよく、モニタの代替としてプリンタなどの出力装置を用いてもよい。
【００１８】
図２は、本実施例１の推定システム１００の機能ブロック図である。
【００１９】
信号取得部２０４は、ＡＤ変換器１０２からチャンネルの入力信号を入力し、周波数領域信号を出力する。状態推定部２０２は、周波数領域信号と複数の源信号状態モデルを入力して、物体種別の推定結果を算出する。そして、結果出力部２０３は、物体種別の推定結果と警告表示を出力する。
【００２０】
信号取得部２０４は、予め設定されたサンプリングの時間間隔δ秒ごとに、ＡＤ変換器１０２から入力信号ｚ（ｔ）を受信し、入力信号ｚ（ｔ）をバッファリングする。ただし、入力信号ｚ（ｔ）は時刻ｔ秒での入力信号を意味する。
【００２１】
信号取得部２０４は、バッファリングされた過去Ｍ回にわたっての入力信号系列［ｚ（ｔ−（Ｍ−１）δ），ｚ（ｔ−（Ｍ−２）δ），…，ｚ（ｔ−δ），ｚ（ｔ）］に対して、窓関数の乗算と短時間フーリエ変換を行い、周波数領域信号ｘ（τ）を次式により計算する。
ｘ（τ）＝［ｘ（τ，１），…，ｘ（τ，Ｆ）］
ただし、δはサンプリングの時間間隔、Ｆは短時間フーリエ変換の周波数番号の個数、τは窓関数の乗算と短時間フーリエ変換の対象時間の範囲（ｔ−（Ｍ−１）δからｔまでの時刻帯）を一つの時間フレームとみなしたときのフレーム番号を表す。
【００２２】
状態推定部２０２は、源信号状態モデルデータベース２０１から読み出した複数の源信号状態モデルと、周波数領域信号ｘ（τ）＝［ｘ（τ，１），…，ｘ（τ，Ｆ）］に基づいて、状態推定処理を実行し、各種別の物体である確率からなる推定結果Ｑ（τ）を次式により出力する。
Ｑ（τ）＝［Ｑ（τ，１），…，Ｑ（τ，Ｋ）］
ただし、Ｑ（τ，ｋ）はフレーム番号τにおいてｋ番目の種別の物体が存在する確率を表し、Ｋは物体の種別の個数を表す。
【００２３】
結果出力部２０３は、推定結果Ｑ（τ）＝［Ｑ（τ，１），…，Ｑ（τ，Ｋ）］をユーザに提示する。提示方法は、たとえば、ユーザインタフェース部１０４を介した画像情報提示、点字ディスプレイによる提示、音声による提示、プリンタを介した画像情報の印刷などである。
【００２４】
信号取得部２０４と、状態推定部２０２と、結果出力部２０３の各機能部はプログラムとして揮発性メモリ１０６にロードされる。
【００２５】
中央演算装置１０３は、各機能部のプログラムに従って処理することによって、所定の機能を提供する機能部として稼働する。例えば、中央演算装置１０３は、信号取得プログラムに従って処理することで信号取得部２０４として機能する。他のプログラムについても同様である。さらに、中央演算装置１０３は、各プログラムが実行する複数の処理のそれぞれの機能を提供する機能部としても稼働する。計算機及び計算機システムは、これらの機能部を含む装置及びシステムである。
【００２６】
信号取得部２０４や状態推定部２０２等の各機能を実現するプログラム、テーブル等の情報は、ストレージサブシステムや不揮発性半導体メモリ、ハードディスクドライブ、ＳＳＤ（Ｓｏｌｉｄ  Ｓｔａｔｅ  Ｄｒｉｖｅ）等の記憶デバイス、または、ＩＣカード、ＳＤカード、ＤＶＤ等の計算機読み取り可能な非一時的データ記憶媒体に格納することができる。
【００２７】
図３は、本実施例１の推定システム１００で行われる処理の一例を示すフローチャートである。
【００２８】
まず、本実施例１の推定システム１００は、ステップＳ３０１で、停止操作が行われたか否かを判断し、測定停止操作が行われていなければステップＳ３０２以降の処理を実行する。測定停止操作が行われていれば推定システム１００は処理を終了する。
【００２９】
次に、ステップＳ３０２の信号取得処理では、信号取得部２０４が、入力信号ｚ（ｔ）をバッファリングし、バッファリングした入力信号系列［ｚ（ｔ−（Ｍ−１）δ），ｚ（ｔ−（Ｍ−２）δ），…，ｚ（ｔ−δ），ｚ（ｔ）］に対して窓関数の乗算と短時間フーリエ変換を行い、周波数領域信号ｘ（τ）＝［ｘ（τ，１），…，ｘ（τ，Ｆ）］を出力する。
【００３０】
次に、ステップＳ３０３の状態推定処理では、状態推定部２０２が、信号状態モデルデータベース２０１から読み出した複数の源信号状態モデルと、周波数領域信号ｘ（τ）＝［ｘ（τ，１），…，ｘ（τ，Ｆ）］に基づいて、状態推定処理を実行する。
【００３１】
状態推定処理は、物体種別推定結果Ｑ（τ，ｋ）からなる推定結果Ｑ（τ）＝［Ｑ（τ，１），…，Ｑ（τ，Ｋ）］を出力する。
【００３２】
次に、ステップＳ３０４の結果出力処理では、結果出力部２０３が、推定結果を出力し、Ｓ３０１に戻って上記処理を繰り返す。
【００３３】
図４は、本実施例の結果出力部２０３が、ユーザインタフェース部１０４を用いて、推定結果を提示した画面イメージである。ユーザインタフェース部１０４に出力される画面イメージ４００は、物体種別推定結果４１０と警告表示４２０の２つの領域を含む。
【００３４】
前述のとおり、本実施例１での推定結果は、物体種別推定結果からなる。物体種別推定結果は、各種別の物体である確率からなる。この提示により、ユーザは、特定の種別の物体が存在するか否かを容易に知ることができるので、この情報に基づいて、次の行動を決定することができるという効果がある。また、結果出力部２０３は、予め設定した物体種別の一致度（確率）が所定の閾値を超えていた場合、警告表示として“要注意”などのメッセージを出力し、警告表示４２０内に表示する。ユーザは警告表示４２０を参照することで、脅威の度合を知ることができる。
【００３５】
図５Ａおよび図５Ｂは、源信号状態モデルデータベース２０１のデータ構造の一例を示す図である。源信号状態モデルデータベース２０１は、物体種別kごとに異なった源信号状態モデルｕ（ｋ）を割り当て、データベースに格納している。それぞれの源信号状態モデルｕ（ｋ）が隠れマルコフモデル（ＨＭＭ：Hidden Markov Model）であり、Ｓ（ｋ）個の状態ｓ（ｋ，１），…，ｓ（ｋ，Ｓ（ｋ））を持つ。
【００３６】
源信号状態モデルｕ（ｋ）は、物体種別ｋがｉ番目の状態ｓ（ｋ，ｉ）において生成する源信号の周波数スペクトルの確率分布ｄｉｓｔ（ｋ，ｉ）、および、物体種別ｋが次の時刻フレームにｉ番目の状態ｓ（ｋ，ｉ）からｊ番目の状態ｓ（ｋ，ｊ）に遷移する遷移確率ｑ（ｋ，ｉ，ｊ）を持つ。
【００３７】
周波数スペクトルの確率分布ｄｉｓｔ（ｋ，ｉ）はたとえば混合ガウスモデル（ＧＭＭ）であり、ｄｉｓｔ（ｋ，ｉ）はＣ（ｋ，ｉ）個のクラスタ［ｃ（ｋ，ｉ，１），…，ｃ（ｋ，ｉ，Ｃ（ｋ，ｉ））］を持つ。ｃ＿ｌ番目のクラスタｃ（ｋ，ｉ，ｃ＿ｌ）は、下記（１）式のスペクトルの平均μ（ｋ，ｉ，ｃ＿ｌ）
【００３８】
【数１】
2016065764.tif 000003
と、下記（１）式のスペクトルの共分散行列Σ（ｋ，ｉ，ｃｌ）
【００３９】
【数２】
2016065764.tif 000004
により定義される。
【００４０】
例えば、図５Ａに示す、物体種別「ダイバー」（ｋ＝１）の信号状態モデルｕ（ｋ＝１）は、「吸気」に相当する「状態１」（ｓ（ｋ＝１，ｉ＝１））と、「呼気」に相当する「状態２」（ｓ（ｋ＝１，ｉ＝２））の２個の状態を有する。「状態１」の周波数スペクトル情報は１０ｋＨｚ以上の比較的高い周波数の成分を多く含むスペクトルであり、「状態２」の周波数スペクトル情報は１０ｋＨｚ以下の比較的低い周波数の成分に偏ったスペクトルである。
【００４１】
また、「ダイバー」が「呼気」と「吸気」を行う確率は各時刻で独立ではなく、ある時刻の状態が「呼気」である場合には、前後の時刻も「呼気」である可能性が高いこと、ある時刻の状態が「吸気」である場合には、前後の時刻も「吸気」である可能性が高いこと、その上で、ある遷移確率で「呼気」から「吸気」、「吸気」から「呼気」に状態を遷移すること、という時間的な性質が隠れマルコフモデルにより表現される。
【００４２】
例えば、図５Ｂに示す、物体種別「小型船舶」の信号状態モデルは、その物体の各航行速度に応じた「低速」に相当する「状態１」、「中速」に相当する「状態２」、「高速」に相当する「状態３」、「最高速」に相当する「状態４」の４個の状態を有する。
【００４３】
このように複数の信号状態を有するモデル化により、船舶のように同じ物体であっても航行速度が大きい場合と航行速度が小さい場合で異なる周波数特性を持つといった性質を表すことができる。また、船舶の速度は各時刻で独立ではなく、速度が時間的に或る程度連続すること、すなわち、前時刻の速度と近い速度に相当する状態へと遷移する確率が高いという時間的な性質が、隠れマルコフモデルにより表現される。
【００４４】
図６は、本実施例１の図３のステップＳ３０３に示した状態推定処理の一例を示すフローチャートである。
【００４５】
ステップＳ６０１のモデル比較処理では、状態推定部２０２が、前記の周波数領域信号ｘ（τ）＝［ｘ（τ，１），…，ｘ（τ，Ｆ）］をバッファリングし、周波数−時間の２次元信号（スペクトログラム）Ｘ（τ）＝［ｘ（τ−Ｔ＋１），…，ｘ（τ−１），ｘ（τ）］に変換する。
【００４６】
ただし、Ｔはバッファリングするフレーム数を意味する。状態推定部２０２は、このスペクトログラムＸ（τ）に対し、源信号状態モデルデータベース２０１から読みだした各物体種別ｋの源信号状態モデルｍ（ｋ）との時系列マッチングを行う。
【００４７】
時系列マッチングはたとえばＶｉｔｅｒｂｉアルゴリズムを用いて実行する。このように、スペクトログラムに対し、複数の状態を有し、かつ、上述した時間的な連続性を表現する源信号状態モデルとの時系列マッチングを行うことで、ダイバーのように「呼気」と「吸気」の複数の源信号状態を持つ物体種別や、船舶のように航行速度によって複数の源信号状態を持つ物体種別であっても、状態推定が可能である。
【００４８】
また、源信号状態モデルｍ（ｋ）が隠れマルコフモデルであるので、時系列マッチングの結果として、各物体種別ｋに対応する源信号状態モデルｍ（ｋ）から前記スペクトログラムＸ（τ）が生成される対数尤度Ｌｓ（τ，ｋ）が算出される。
【００４９】
【数３】
2016065764.tif 000005
【００５０】
上記（３）式により、Ｌｓ（τ，ｋ）の指数をとれば、フレームτ−Ｔ＋１〜τの間に存在する音源が物体種別ｋである確率ｐｓ（τ，ｋ）が分かる。そして、Ｑ（τ，ｋ）＝ｐｓ（τ，ｋ）を総合指標として、
推定結果Ｑ（τ）＝［Ｑ（τ，１），…，Ｑ（τ，Ｋ）］
として出力する。
【００５１】
本実施例１は、以上の構成により、源信号の状態が時間変化する場合であっても、信号特徴の状態を特定することが可能となる。ただし、上記実施例１では、物体の対象数が１で、かつ、軌跡を特定しない場合に物体の状態を特定し、物体の種別を推定することができる。
【実施例２】
【００５２】
本実施例２では、源信号の信号特徴の状態が時間変化する場合であっても、信号特徴の状態の特定、および、物体の軌跡の追跡が可能な、推定システム１００の例を説明する。
【００５３】
図７は、本実施例２の推定システム１００のハードウェア構成の一例を示すブロック図である。本実施例２の推定システム１００は、ハイドロホンやマイクロホンなどの素子１０１＿１〜１０１＿Ｎ、複数チャンネルＡＤ変換器１０２Ａ、中央演算装置１０３、ユーザインタフェース部１０４、記憶媒体１０５、揮発性メモリ１０６から構成されている。
素子１０１＿１〜１０１＿Ｎは、水中あるいは空中の音波をアナログ電圧値に変換し、複数チャンネルＡＤ変換器７０２に送る。複数チャンネルＡＤ変換器１０２Ａは、複数チャンネルのアナログ電圧値を複数チャンネルのデジタル信号に変換し、中央演算装置１０３に送る。なお、その他のハードウェアの構成は、前記実施例１と同様であり、同一のものに同一の符号を付して重複した説明は省略する。
【００５４】
図８は、本実施例２の推定システム１００の機能の一例を示すブロック図である。
【００５５】
前記実施例１の図２と異なる機能要素は、方向推定部８０１と軌跡推定部８０２を加えた点が主な相違点である。信号取得部８０５がＡＤ変換器１０２Ａから複数チャンネルの入力信号を入力し、複数チャンネルの周波数領域信号を出力する。複数チャンネルの周波数領域信号は方向推定部８０１と、状態推定部８０３へ出力される。
【００５６】
方向推定部８０１は、物体の方向推定結果を算出して軌跡推定部８０２へ出力する。軌跡推定部８０２は、方向推定結果に基づいて軌跡の仮説集合を算出し、状態推定部８０３へ出力する。状態推定部８０３は、複数チャンネルの周波数領域信号と軌跡仮説集合と複数の源信号状態モデルを入力して、推定結果を算出する。そして、結果出力部８０４は、前記実施例１の物体種別の推定結果と警告表示に加えて、物体の軌跡と、音源分離結果を出力する。
【００５７】
信号取得部８０５は、予め設定されたサンプリングの時間間隔δ秒ごとに、複数チャンネルＡＤ変換器１０２Ａから複数チャンネル入力信号ｚ（ｎ，ｔ）を受信し、ｚ（ｎ，ｔ）をバッファリングする。ただし、ｚ（ｎ，ｔ）は時刻ｔ秒での素子ｎの入力信号を意味する。
【００５８】
バッファリングされた過去Ｍ回にわたっての入力信号系列［ｚ（ｎ，ｔ−（Ｍ−１）δ），ｚ（ｎ，ｔ−（Ｍ−２）δ），…，ｚ（ｎ，ｔ−δ），ｚ（ｎ，ｔ）］に対して、窓関数の乗算と短時間フーリエ変換を行い、複数チャンネル周波数領域信号ｘ（ｎ，τ）＝［ｘ（ｎ，τ，１），…，ｘ（ｎ，τ，Ｆ）］を計算し、出力する。
【００５９】
ただし、δはサンプリングの時間間隔、Ｆは短時間フーリエ変換の周波数番号の個数、τは窓関数の乗算と短時間フーリエ変換の対象時間の範囲（ｔ−（Ｍ−１）δからｔまでの時刻帯）を一つの時間フレームとみなしたときのフレーム番号を表す。
【００６０】
方向推定部８０１は、複数チャンネル周波数領域信号ｘ（ｎ，τ）を入力として、方向推定処理を実行し、時刻τの方向推定結果ｙ（τ）＝［ｙ（τ，１），…，ｙ（τ，Ｄ）］を出力する。ただし、Ｄは離散化した方向番号の個数である。
【００６１】
軌跡推定部８０２は、複数時刻の方向推定結果［ｙ（τ−（Ｂ−１）），ｙ（τ−（Ｂ−２）），…，ｙ（τ）］に基づいて、軌跡推定処理を実行し、軌跡仮説集合ｈ（τ）＝［ｈ（τ，１），ｈ（τ，２），…，ｈ（τ，Ｈ）］を出力する。ただし、Ｈは仮説集合の要素数である。また、ｈ（τ，ｂ）はｂ番目の軌跡仮説を意味する。軌跡仮説ｈ（τ，ｂ）は音源方向の時系列であり、
ｈ（τ，ｂ）＝［θ（ｂ，τ−（Ｂ−１）），θ（ｂ，τ−（Ｂ−２）），…，θ（ｂ，τ）］
と表す。
【００６２】
状態推定部８０３は、源信号状態モデルデータベース２０１から読み出した複数の源信号状態モデルと、軌跡仮説集合ｈ（τ）、および、複数チャンネル周波数領域信号ｘ（ｎ，τ）に基づいて、状態推定処理を実行する。
【００６３】
状態推定処理は、各種別の物体である確率Ｑ（τ，ｋ）からなる物体種別推定結果Ｑ（τ）＝［Ｑ（τ，１），…，Ｑ（τ，Ｋ）］と軌跡Ｔｒａｊ（τ）の組からなる推定結果（Ｑ（τ），Ｔｒａｊ（τ））を出力する。
【００６４】
ただし、軌跡Ｔｒａｊ（τ）は、軌跡仮説集合に基づいて推定された、軌跡仮説と同様の方向の時系列であり、［θ￣（τ−（Ｂ−１）），θ￣（τ−（Ｂ−２）），…，θ￣（τ）］と表す。
【００６５】
結果出力部８０４は、推定結果（Ｑ（τ），Ｔｒａｊ（τ））をユーザに提示する。提示方法は、たとえば、ユーザインタフェース部１０４を介した画像情報提示、点字ディスプレイによる提示、音声による提示、プリンタを介した画像情報の印刷などである。
【００６６】
図９は、本実施例２の推定システム１００で行われる処理の一例を示すフローチャートである。
【００６７】
まず、本実施例の推定システム１００は、ステップＳ３０１で、停止操作が行われたか否かを判断し、測定停止操作が行われていなければステップＳ９０１以降の処理を実行する。測定停止操作が行われていれば推定システム１００は処理を終了する。
【００６８】
次に、ステップＳ９０１の信号取得処理では、信号取得部８０５が、複数チャンネル入力信号ｚ（ｎ，ｔ）をバッファリングし、バッファリングした入力信号系列［ｚ（ｎ，ｔ−（Ｍ−１）δ），ｚ（ｎ，ｔ−（Ｍ−２）δ），…，ｚ（ｎ，ｔ−δ），ｚ（ｎ，ｔ）］に対して窓関数の乗算と短時間フーリエ変換を行い、複数チャンネル周波数領域信号ｘ（ｎ，τ）＝［ｘ（ｎ，τ，１），…，ｘ（ｎ，τ，Ｆ）］を出力する。
【００６９】
次に、ステップＳ９０１の方向推定処理では、方向推定部８０１が、複数チャンネル周波数領域信号ｘ（ｎ，τ）に基づいて音源方向の推定を行う。音源方向の推定には、遅延和アレー法、最小分散ビームフォーマ法、ＧＣＣ−ＰＨＡＴ法、ＳＲＰ−ＰＨＡＴ法、ＭＵＳＩＣ法、ＥＳＰＲＩＴ法、ＳＰＩＲＥ法など周知または公知の技術が存在するが、いずれを用いても良い。
【００７０】
音源方向の推定結果は、方向ヒストグラムと呼ばれる、各時刻における各方向の音源の音量の大きさを表す２次元の信号ｙ（τ，θ）として得られる。たとえば、最小分散ビームフォーマ法の場合、下記の（４）式により、方向ヒストグラムを計算する。
【００７１】
【数４】
2016065764.tif 000006
【００７２】
ただし、ａ（ｆ，θ）＝［ａ（ｆ，θ，１），…，ａ（ｆ，θ，Ｎ）］は、周波数ｆにおける方向θのステアリングベクトルを表し、Ｒ（ｆ，τ）は周波数ｆにおける素子間共分散行列を表す。Ｒ（ｆ，τ）は、次の（５）式で計算される。
【００７３】
【数５】
2016065764.tif 000007
ただし、γは０以上１以下の定数である。
【００７４】
上記（４）式のａ（ｆ，θ）の各要素ａ（ｆ，θ，ｎ）は素子配置に基づくシミュレーションによって、次の（６）式で計算される。
【００７５】
【数６】
2016065764.tif 000008
【００７６】
ただし、η（ｎ）は素子１と素子ｎの間での源信号の到達時間の差であり、ＦｒｅｑＭａｘは周波数番号Ｆに対応する周波数である。
【００７７】
次に、ステップＳ９０３の軌跡推定処理では、軌跡推定部８０２が、複数時刻の方向推定結果［ｙ（τ−（Ｂ−１）），ｙ（τ−（Ｂ−２）），…，ｙ（τ）］を入力として、とりうると考えられる軌跡の集合である軌跡仮説集合ｈ（τ）＝［ｈ（τ，１），ｈ（τ，２），…，ｈ（τ，Ｈ）］を出力する。
【００７８】
具体的には、まず、各時刻τの方向推定結果の方向ヒストグラムｙ（τ）＝［ｙ（τ，１），…，ｙ（τ，Ｄ）］に対してｋ−ｍｅａｎｓなどのクラスタリングアルゴリズムを実行する。
【００７９】
クラスタリングの結果得られる代表点θ（τ，α＝１），…θ（τ，α＝Ａ）を音源方向の候補とみなす。αはクラスタの番号であり、Ａはクラスタ数を表す定数である。各時刻τのある音源方向候補θ（τ，α）と前の時刻τ−１のある音源方向候補θ（τ−１，β）との間をつないでいくと、全時間にわたっての音源方向の系列を得ることができる。
【００８０】
この系列一つを軌跡仮説と呼び、ｂ番目の軌跡仮説をｈ（τ，ｂ）と表す。そして、とりうる全ての場合の軌跡仮説ｈ（τ，ｂ）を列挙した集合を、軌跡仮説集合ｈ（τ）＝［ｈ（τ，１），ｈ（τ，２），…，ｈ（τ，Ｈ）］として出力する。
【００８１】
次に、ステップＳ９０４の状態推定処理では、状態推定部８０３が、信号状態モデルデータベース２０１から読み出した複数の源信号状態モデルと、軌跡仮説集合ｈ（τ）、および、複数チャンネル周波数領域信号ｘ（ｎ，τ）＝［ｘ（ｎ，τ，１），…，ｘ（ｎ，τ，Ｆ）］に基づいて、状態推定処理を実行し、推定した軌跡Ｔｒａｊ（τ）と、その軌跡の物体の源信号を推定した音源分離後の周波数領域信号ｓ（τ）＝［ｓ（τ，１），…，ｓ（τ，Ｆ）］と、各種別の物体である確率（物体種別推定結果と呼ぶ）Ｑ（τ）＝［Ｑ（τ，１），…，Ｑ（τ，Ｋ）］からなる、推定結果（Ｔｒａｊ（τ），ｓ（τ），Ｑ（τ））を出力する。
【００８２】
次に、ステップＳ９０５の結果出力処理では、結果出力部８０４が、前記の推定結果を出力し、ステップＳ３０１に戻って上記処理を繰り返す。
【００８３】
図１０は、本実施例２の結果出力部８０４が、ユーザインタフェース部１０４を用いて、推定結果（Ｔｒａｊ（τ），ｓ（τ），Ｑ（τ））を提示した画面イメージ４００Ａである。ユーザインタフェース部１０４に出力される画面イメージ４００Ａは、物体種別推定結果４１０と、警告表示４２０と、軌跡４３０及び音源分離結果４４０の４つの領域を含む。
【００８４】
前述のとおり、推定結果は、推定した軌跡Ｔｒａｊ（τ）を、例えば、２次元のマップで表示する軌跡４３０と、その軌跡の物体の源信号を推定した音源分離後信号ｓ（τ）をグラフとして表示する音源分離結果４４０と、物体種別推定結果Ｑ（τ）を物体の名称と一致度（確率）で表示する物体種別推定結果４１０からなる。物体種別推定結果Ｑ（τ）は、各種別ｋの物体である確率Ｑ（τ，ｋ）からなる。
【００８５】
推定した軌跡と、その軌跡の物体の源信号を推定した音源分離結果と、物体種別推定結果を同時に提示することで、ユーザは、特定の種別の物体が、特定の軌跡で移動しているという情報を容易に知ることができる。ユーザはこの情報に基づいて、次の行動を決定することができるという効果がある。また、前記実施例１と同様に、警告表示４２０の内容に基づいて、ユーザは脅威の度合を知ることができる。
【００８６】
図１１は、本実施例２の図９に示したステップＳ９０４で行われる状態推定処理の一例を示すフローチャートである。
【００８７】
まず、ステップＳ１１０１の運動尤度計算処理では、状態推定部８０３が、軌跡仮説集合ｈ（τ）＝［ｈ（τ，１），ｈ（τ，２），…，ｈ（τ，Ｈ）］の中の全ての軌跡仮説ｈ（τ，ｂ）に対し、ｈ（τ，ｂ）を構成する方向の時系列［θ（τ−（Ｂ−１）），θ（τ−（Ｂ−２）），…，θ（τ）］に基づいて、仮説ｂが正しい場合の真の音源位置の時系列ｔｒａｊ（τ，ｂ）＝［θ￣（τ−（Ｂ−１）），θ￣（τ−（Ｂ−２）），…，θ￣（τ）］を推定する。
【００８８】
具体的には、軌跡仮説を構成する音源方向θ（τ）を観測信号として、カルマンフィルタを実行することで、真の音源方向θ￣（τ）を計算する。たとえば、カルマンフィルタのシステムモデルと観測モデルは、それぞれ、次の（７）式及び（８）式で表される。
【００８９】
【数７】
2016065764.tif 000009
【００９０】
また、Δτは時間フレームの間隔を表し、Ｒξは或る共分散行列を表し、σψは或る分散値を表す。
【００９１】
【数８】
2016065764.tif 000010
【００９２】
ここで、カルマンフィルタの代わりにパーティクルフィルタなど他の追跡アルゴリズムを用いてもよい。ステップＳ１１０１の運動尤度計算処理ではさらに、次の（９）式により、ｈ（τ，ｂ）＝［θ（τ−（Ｂ−１）），θ（τ−（Ｂ−２）），…，θ（τ）］が観測された条件の下で、仮説ｂが正しい場合の真の音源位置の時系列ｔｒａｊ（τ，ｂ）が［θ￣（τ−（Ｂ−１）），θ￣（τ−（Ｂ−２）），…，θ￣（τ）］である尤度（運動尤度）Ｌｍ（τ，ｂ）を計算する。
【００９３】
【数９】
2016065764.tif 000011
【００９４】
次に、ステップＳ１１０２の信号分離処理では、状態推定部８０３が、軌跡仮説ｈ（τ，ｂ）ごとに、推定された音源位置の時系列ｔｒａｊ（τ，ｂ）＝［θ￣（τ−（Ｂ−１）），θ￣（τ−（Ｂ−２）），…，θ￣（τ）］と複数チャンネル周波数入力信号ｘ（ｎ，τ）に基づいて、軌跡仮説ｈ（τ，ｂ）ごとの各時刻τのθ￣（τ）方向から到来する音のみを抽出する。
【００９５】
音源分離には、たとえば、遅延和アレー法、最小分散ビームフォーマ法、最大ＳＮＲビームフォーマ法、局所ガウスモデリング法などの周知または公知の技術を用いればよい。ステップＳ１１０２の信号分離処理では、軌跡仮説ｂごとに対応する方向θ￣（τ）から到来した音の周波数領域信号ｓ（ｂ，τ）＝［ｓ（ｂ，τ，１），…，ｓ（ｂ，τ，Ｆ）］を出力するものとする。
【００９６】
状態推定部８０３が、予め各方向θの音を抽出する音源分離フィルタをシミュレーションにより学習しておけば、音源分離処理の実行時には、指定された方向に対応する音源分離フィルタを読み出すことで、低計算コストで移動音源に対する音源分離が可能である。
【００９７】
次に、ステップＳ１１０３のモデル比較処理では、状態推定部８０３が、前記の軌跡仮説ｂごとに対応する周波数領域信号ｓ（ｂ，τ）＝［ｓ（ｂ，τ，１），…，ｓ（ｂ，τ，Ｆ）］をバッファリングし、軌跡仮説ｂごとの周波数−時間の２次元信号（スペクトログラム）Ｓ（ｂ，τ）＝［ｓ（ｂ，τ−Ｔ＋１），…，ｓ（ｂ，τ−１），ｓ（ｂ，τ）］に変換する。ただし、Ｔはバッファリングするフレーム数を意味する。
【００９８】
このスペクトログラムＳ（ｂ，τ）に対し、源信号状態モデルデータベース２０１から読みだした各物体種別ｋの源信号状態モデルｍ（ｋ）との時系列マッチングを行う。時系列マッチングはたとえばＶｉｔｅｒｂｉアルゴリズムを用いて実行する。このように、スペクトログラムに対し、複数の状態を有し、かつ、上述した時間的な連続性を表現する源信号状態モデルとの時系列マッチングを行うことで、ダイバーのように「呼気」と「吸気」の複数の源信号状態を持つ物体種別や、船舶のように航行速度によって複数の源信号状態を持つ物体種別であっても、状態推定が可能である。
【００９９】
また、源信号状態モデルｍ（ｋ）が隠れマルコフモデルであるので、時系列マッチングの結果として、各物体種別ｋに対応する源信号状態モデルｍ（ｋ）から前記スペクトログラムＳ（ｂ，τ）が生成される対数尤度Ｌｓ（ｂ，τ，ｋ）が算出される。
【０１００】
次の、（１０）式により、Ｌｓ（ｂ，τ，ｋ）の指数をとれば、軌跡仮説ｂが正しい場合にフレームτ−Ｔ＋１〜τの間に存在する音源が物体種別ｋである確率ｐｓ（ｂ，τ，ｋ）が得られる。
【０１０１】
【数１０】
2016065764.tif 000012
【０１０２】
そして、次の、（１１）式により、フレームτ−Ｔ＋１〜τの間でどの軌跡仮説ｂが尤もらしいかの確率に対応するΩ（τ，ｂ）を計算する。
【０１０３】
【数１１】
2016065764.tif 000013
そして、次の、（１２）式により、まずΩ（τ，ｂ）に基づいてフレームτ−Ｔ＋１〜τの間の最尤の軌跡仮説ｂ＾を計算し、最尤の軌跡仮説ｂ＾が正しい場合に物体の種別がｋである確率Ｑ（τ，ｋ）＝ｐｓｓ（ｂ＾，τ，ｋ）を計算し、各種別の物体である確率からなる物体種別推定結果Ｑ（τ）＝［Ｑ（τ，１），…，Ｑ（τ，Ｋ）］を計算する。
【０１０４】
【数１２】
2016065764.tif 000014
【０１０５】
誤った軌跡仮説ｂにおいては，いずれの物体種別ｋに対するｐｓ（ｂ，τ，ｋ）も小さくなるため、（１１）式のΩ（τ，ｂ）が小さくなる。そのため、（１２）式においてΩ（τ，ｂ）が最大となるｂを選択することにより、誤った軌跡仮説が棄却できる。
【０１０６】
これにより、正しい軌跡仮説のみに絞りこむことができ、源信号の状態が時間変化する場合であっても高い追跡精度が得られる。また、音源の物体種別を推定する効果があるので、推定した物体種別に応じた行動を行うことができる。
【０１０７】
最後に、最尤の軌跡仮説ｂ＾を用いて、推定軌跡Ｔｒａｊ（τ）＝ｔｒａｊ（τ，ｂ＾）、および、音源分離後の周波数領域信号ｓ（τ）＝ｓ（ｂ＾，τ）を計算する。そして、それらの組として推定結果（Ｔｒａｊ（τ），ｓ（τ），Ｑ（τ））を出力する。
【０１０８】
本実施例２は、以上の構成により、源信号の状態が時間変化する場合であっても、信号特徴の状態の特定、および、物体の軌跡の追跡が可能となる。
【実施例３】
【０１０９】
本実施例３では、物体が複数存在し、源信号の信号特徴の状態が時間変化する場合であっても、信号特徴の状態の特定、および、軌跡の追跡が可能な、推定システム１００の例を説明する。
【０１１０】
以下は実施例２と異なる点について説明する。本実施例３の推定システム１００のハードウェア構成は実施例２に示した図７と同様である。
【０１１１】
図１２は、本実施例３の推定システム１００の機能ブロック図である。図１２において、実施例２の図８と異なる機能要素は、軌跡推定部１２０１と状態推定部１２０２と結果出力部１２０３である。
【０１１２】
軌跡推定部１２０１は、方向推定結果に加えて物体数を入力とし、複数物体の軌跡仮説集合を出力する。状態推定部１２０２は、複数物体の軌跡仮説集合と複数の源信号状態モデルを入力として、複数物体の推定結果を出力する。結果出力部１２０３は複数物体の推定結果を入力として、複数の推定結果を出力する。
【０１１３】
軌跡推定部１２０１は、複数時刻の方向推定結果［ｙ（τ−（Ｂ−１）），ｙ（τ−（Ｂ−２）），…，ｙ（τ）］に基づいて、軌跡推定処理を実行し、複数物体軌跡仮説集合ｈ（τ）＝［ｈ（τ，１），ｈ（τ，２），…，ｈ（τ，Ｈ）］を出力する。
【０１１４】
ただし、Ｈは複数物体仮説集合の要素数である。また、ｈ（τ，ｂ）はｂ番目の複数物体軌跡仮説を意味する。
【０１１５】
複数物体軌跡仮説ｈ（τ，ｂ）は各物体の音源方向の組の時系列であり、
ｈ（τ，ｂ）＝［［θ１（ｂ，τ−（Ｂ−１）），…，θＧ（ｂ，τ−（Ｂ−１））］，…，［θ１（ｂ，τ），…，θＧ（ｂ，τ）］］
と表す。ただし、Ｇは物体数である。
【０１１６】
状態推定部１２０２は、源信号状態モデルデータベース２０１から読み出した複数の源信号状態モデルと、複数物体軌跡仮説集合ｈ（τ）、および、複数チャンネル周波数領域信号ｘ（ｎ，τ）に基づいて、状態推定処理を実行する。
【０１１７】
状態推定処理は、各物体ｇが各種別の物体である確率Ｑ（τ，ｇ，ｋ）からなる物体種別推定結果Ｑ（τ）＝［［Ｑ（τ，１，１），…，Ｑ（τ，１，Ｋ）］，…，［Ｑ（τ，Ｇ，１），…，Ｑ（τ，Ｇ，Ｋ）］］と、軌跡Ｔｒａｊ（τ）の組からなる推定結果（Ｑ（τ），Ｔｒａｊ（τ））を出力する。
【０１１８】
ただし、軌跡Ｔｒａｊ（τ）は、複数物体軌跡仮説集合に基づいて推定された、複数物体軌跡仮説と同様の方向の時系列であり、［［θ１￣（τ−（Ｂ−１）），…，θＧ￣（τ−（Ｂ−１））］，…，［θ１￣（τ），…，θＧ￣（τ）］］と表す。
【０１１９】
結果出力部１２０３は、推定結果（Ｑ（τ），Ｔｒａｊ（τ））と警告表示をユーザに提示する。提示方法は、たとえば、ユーザインタフェース部１０４を介した画像情報提示、点字ディスプレイによる提示、音声による提示、プリンタを介した画像情報の印刷などである。
【０１２０】
図１３は、本実施例３の推定システム１００で行われる処理の一例を示すフローチャートである。前記実施例２で示した図９のフローチャートとの違いは、ステップＳ１３０１、Ｓ１３０２、Ｓ１３０３を新たに加えた点である。以下では、前記実施例２の図９との相違点について説明する。
【０１２１】
ステップＳ１３０１では、状態推定部１２０２が、物体数Ｇを１に初期化する。ステップＳ１３０２では、物体数Ｇが所定の閾値ＧＭａｘ以下であればステップＳ９０１〜Ｓ９０４、Ｓ１３０３以降の処理を実行し、物体数Ｇが所定の閾値ＧＭａｘを超えていればＳ９０５を実行する。
【０１２２】
ステップＳ１３０３では、物体数Ｇに１を加えてステップＳ１３０２を実行する。このループにより、Ｇ＝１，…，ＧＭａｘのそれぞれの物体数で事後確率を評価し、事後確率が最大となる物体数における最尤の推定結果を出力できる。
【０１２３】
本実施例３のステップＳ９０３で行われる軌跡推定処理を説明する。軌跡推定部１２０１が、複数時刻の方向推定結果［ｙ（τ−（Ｂ−１）），ｙ（τ−（Ｂ−２）），…，ｙ（τ）］を入力として、とりうると考えられる軌跡の集合である軌跡仮説集合ｈ（τ）＝［ｈ（τ，１），ｈ（τ，２），…，ｈ（τ，Ｈ）］を出力する。
【０１２４】
具体的には、まず、各時刻τの方向推定結果の方向ヒストグラムｙ（τ）＝［ｙ（τ，１），…，ｙ（τ，Ｄ）］に対してｋ−ｍｅａｎｓなどのクラスタリングアルゴリズムを実行する。
【０１２５】
クラスタリングの結果得られる代表点θ（τ，α＝１），…θ（τ，α＝Ａ）を音源方向の候補とみなす。ここで、αはクラスタの番号であり、Ａはクラスタ数を表す定数である。
【０１２６】
各時刻τにおいて、Ｇ個の物体のそれぞれに対応する音源方向候補θ（τ，α）を選んだ組み合わせγ（τ，ｆ）＝［θ（τ，αｆ１），…，θ（τ，αｆＧ）］を列挙し、複数音源候補γ（τ）＝［γ（τ，１），…，γ（τ，Ｆ）］とする。
【０１２７】
ただし、θ（τ，αｆｇ）は，ｆ番目（１≦ｆ≦Ｆ）の組み合わせのｇ番目（１≦ｆ≦Ｇ）の物体の音源方向を表し、Ｆは全ての組み合わせの個数を表す。ただし、同じτ，ｆにおいて互いに異なる物体ｇに対し、異なるαｆｇ、および、異なるθ（τ，αｆｇ）を割り当てるような組み合わせに限定する。
【０１２８】
各時刻τのある複数音源候補γ（τ，ｆ１）と前の時刻τ−１のある複数音源候補γ（τ？１，ｆ２）との間をつないでいくと、全時間にわたっての複数音源方向の系列を得ることができる。
【０１２９】
この系列の一つを複数物体軌跡仮説と呼び、ｂ番目の複数物体軌跡仮説をｈ（τ，ｂ）として表す。そして、とりうる全ての場合の複数物体軌跡仮説ｈ（τ，ｂ）を列挙した集合を、複数物体軌跡仮説集合ｈ（τ）＝［ｈ（τ，１），ｈ（τ，２），…，ｈ（τ，Ｈ）］として出力する。
【０１３０】
本実施例３のステップＳ９０４で行われる状態推定処理について図１１を参照しながら説明する。
【０１３１】
まず、ステップＳ１１０１の運動尤度計算処理で、実施例３の状態推定部１２０２は、複数物体軌跡仮説集合ｈ（τ）＝［ｈ（τ，１），ｈ（τ，２），…，ｈ（τ，Ｈ）］の中の全ての複数物体軌跡仮説ｈ（τ，ｂ）に対し、ｈ（τ，ｂ）の各物体ｇを構成する方向の時系列［θｇ（τ−（Ｂ−１）），θｇ（τ−（Ｂ−２）），…，θｇ（τ）］に基づいて、仮説ｂが正しい場合の真の音源位置の時系列ｔｒａｊｇ（τ，ｂ）＝［θｇ￣（τ−（Ｂ−１）），θｇ￣（τ−（Ｂ−２）），…，θｇ￣（τ）］を推定する。
【０１３２】
具体的には、軌跡仮説を構成する音源方向θｇ（τ）を観測信号として、実施例２と同じくカルマンフィルタを実行することで、真の音源方向θｇ￣（τ）を計算する。
【０１３３】
ステップＳ１１０１の運動尤度計算処理ではさらに、次の（１３）式により、
【０１３４】
【数１３】
2016065764.tif 000015
【０１３５】
ｈ（τ，ｂ）が観測された条件の下で、仮説ｂが正しい場合の真の音源位置の時系列ｔｒａｊｇ（τ，ｂ）が［θｇ￣（τ−（Ｂ−１）），θｇ￣（τ−（Ｂ−２）），…，θ￣（τ）］である尤度（運動尤度）Ｌｍ（τ，ｂ）を計算する。
【０１３６】
次に、ステップＳ１１０２の信号分離処理では、状態推定部１２０２が、複数物体軌跡仮説ｈ（τ，ｂ）ごとに、ｈ（τ，ｂ）の各物体ｇの推定音源位置の時系列ｔｒａｊｇ（τ，ｂ）＝［θｇ￣（τ−（Ｂ−１）），θｇ￣（τ−（Ｂ−２）），…，θｇ￣（τ）］と複数チャンネル周波数入力信号ｘ（ｎ，τ）に基づいて、複数物体軌跡仮説ｈ（τ，ｂ）の物体ｇごとの各時刻τのθｇ￣（τ）方向から到来する音のみを抽出する。
【０１３７】
音源分離には、たとえば、遅延和アレー法、最小分散ビームフォーマ法、最大ＳＮＲビームフォーマ法、局所ガウスモデリング法などを用いればよい。
【０１３８】
ステップＳ１１０２の信号分離処理は、軌跡仮説ｂごと物体ｇごとに対応する方向θｇ￣（τ）から到来した音の周波数領域信号ｓ（ｂ，ｇ，τ）＝［ｓ（ｂ，ｇ，τ，１），…，ｓ（ｂ，ｇ，τ，Ｆ）］を出力するものとする。
【０１３９】
次に、ステップＳ１１０３のモデル比較処理では、前記の軌跡仮説ｂ、物体ｇごとに対応する周波数領域信号ｓ（ｂ，ｇ，τ）＝［ｓ（ｂ，ｇ，τ，１），…，ｓ（ｂ，ｇ，τ，Ｆ）］をバッファリングし、仮説ｂ、物体ｇごとの周波数−時間の２次元信号（スペクトログラム）Ｓ（ｂ，ｇ，τ）＝［ｓ（ｂ，ｇ，τ−Ｔ＋１），…，ｓ（ｂ，ｇ，τ−１），ｓ（ｂ，ｇ，τ）］に変換する。ただし、Ｔはバッファリングするフレーム数を意味する。
【０１４０】
このスペクトログラムＳ（ｂ，ｇ，τ）に対し、源信号状態モデルデータベース２０１から読みだした各物体種別ｋの源信号状態モデルｍ（ｋ）との時系列マッチングを行う。
【０１４１】
時系列マッチングはたとえばＶｉｔｅｒｂｉアルゴリズムを用いて実行する。このように、スペクトログラムに対し、複数の状態を有し、かつ、上述した時間的な連続性を表現する源信号状態モデルとの時系列マッチングを行うことで、ダイバーのように「呼気」と「吸気」の複数の源信号状態を持つ物体種別や、船舶のように航行速度によって複数の源信号状態を持つ物体種別であっても、状態推定が可能である。
【０１４２】
また、源信号状態モデルｍ（ｋ）が隠れマルコフモデルであるので、時系列マッチングの結果として、各物体種別ｋに対応する源信号状態モデルｍ（ｋ）から前記スペクトログラムＳ（ｂ，ｇ，τ）が生成される対数尤度Ｌｓ（ｂ，ｇ，τ，ｋ）が算出される。
【０１４３】
【数１４】
2016065764.tif 000016
【０１４４】
上記（１４）式により、Ｌｓ（ｂ，ｇ，τ，ｋ）の指数をとれば、軌跡仮説ｂが正しい場合に、フレームτ−Ｔ＋１〜τの間に存在する音源のうちの物体ｇが物体種別ｋである確率ｐｓ（ｂ，ｇ，τ，ｋ）が得られる。そして、次の（１５）式により、
【０１４５】
【数１５】
2016065764.tif 000017
【０１４６】
フレームτ−Ｔ＋１〜τの間でどの軌跡仮説ｂが尤もらしいかの確率に対応するΩ（τ，ｂ）を計算する。そして、次の（１６）式により、
【０１４７】
【数１６】
2016065764.tif 000018
【０１４８】
まず、状態推定部１２０２は、Ω（τ，ｂ）に基づいてフレームτ−Ｔ＋１〜τの間の最尤の軌跡仮説ｂ＾を計算し、最尤の軌跡仮説ｂ＾が正しい場合に物体ｇの物体の種別がｋである確率Ｑ（τ，ｇ，ｋ）＝ｐｓ（ｂ＾，ｇ，τ，ｋ）を計算し、各物体ｇが各種別ｋである確率からなる物体種別推定結果Ｑ（τ）＝［［Ｑ（τ，１，１），…，Ｑ（τ，１，Ｋ）］，…，［Ｑ（τ，Ｇ，１），…，Ｑ（τ，Ｇ，Ｋ）］］を計算する。
【０１４９】
誤った軌跡仮説ｂにおいては，各物体ｇのいずれの物体種別ｋに対するｐｓ（ｂ，ｇ，τ，ｋ）も小さくなるため、上記（１５式）のΩ（τ，ｂ）が小さくなる。そのため、（１６）式においてΩ（τ，ｂ）が最大となるｂを選択することにより、誤った軌跡仮説が棄却できる。
【０１５０】
これにより、状態推定部１２０２は、正しい軌跡仮説のみに絞りこむことができ、源信号の状態が時間変化する場合であっても高い追跡精度が得られる。また、音源の物体種別を推定する効果があるので、推定した物体種別に応じた行動を行うことができる。また、本実施例は、複数の物体の軌跡仮説の検証の中に、源信号の状態が時間変化するモデルを導入しているので、源信号の状態が時間変化する複数の物体が、近接または交差するような場合でも、高精度な追跡が可能である。
【０１５１】
次に、状態推定部１２０２は、最尤の複数物体軌跡仮説ｂ＾を用いて、複数物体推定軌跡Ｔｒａｊ（τ）＝［ｔｒａｊ１（τ，ｂ＾），…，ｔｒａｊｇ（τ，ｂ＾），…，ｔｒａｊＧ（τ，ｂ＾）］、および、複数物体の音源分離後の周波数領域信号ｓ（τ）＝［ｓ（ｂ＾，１，τ），…，ｓ（ｂ＾，ｇ，τ），…，ｓ（ｂ＾，Ｇ，τ）］を計算する。
【０１５２】
そして、状態推定部１２０２は、それらの組として、物体数がＧの場合の推定結果Ｗ（Ｇ）＝（Ｔｒａｊ（τ），ｓ（τ），Ｑ（τ））をメモリに記憶する。最後に、図１３のＳ１３０１〜Ｓ１３０３の前回までのループでの物体数Ｇ＝１，…，ＧＭａｘでそれぞれ計算したΩ（τ，ｂ）が最大であるＧについて推定結果Ｗ（Ｇ）を出力する。これにより、存在する物体数を自動的に推定することができる。
【０１５３】
また、ステップＳ１１０３のモデル比較処理では、状態推定部１２０２が、最後に複数物体推定軌跡Ｔｒａｊ（τ）＝［ｔｒａｊ１（τ，ｂ＾），…，ｔｒａｊｇ（τ，ｂ＾），…，ｔｒａｊＧ（τ，ｂ＾）］に応じて源信号状態モデルデータベース２０１の学習を行う。
【０１５４】
すなわち、状態推定部１２０２は、複数物体推定軌跡Ｔｒａｊ（τ）が含む、ある物体ｇの推定軌跡ｔｒａｊｇ（τ，ｂ＾）＝［θｇ￣（τ−（Ｂ−１）），θｇ￣（τ−（Ｂ−２）），…，θｇ￣（τ）］が、物体ｇ以外の他のいずれの物体ｇ’の推定軌跡ｔｒａｊｇ’（τ，ｂ＾）＝［θｇ’￣（τ−（Ｂ−１）），θｇ’￣（τ−（Ｂ−２）），…，θｇ’￣（τ）］とも接近しない場合にのみ、源信号状態モデルデータベース２０１の学習を行う。すなわち、全てのτ’＝τ−（Ｂ−１），…，τ、全てのｇ’（ｇ’≠ｇ）に対して、｜θｇ￣（τ’）−θｇ’￣（τ’）｜＞εである場合のみ、物体ｇがいずれの物体とも接近していないと判断する。
【０１５５】
そして、物体ｇに対応する音源分離後周波数領域信号ｓ（ｂ＾，ｇ，τ）を入力としてＢａｕｍ−Ｗｅｌｃｈなどの学習アルゴリズムによりＨＭＭを学習し、源信号状態モデルデータベース２０１に格納する。このように物体が近接しない場合のみ源信号状態モデルを更新して学習を行うことで、源信号状態モデルの高精度な学習が可能である。
【０１５６】
図１４は、本実施例の結果出力部１２０３が、ユーザインタフェース部１０４を用いて、推定結果（Ｔｒａｊ（τ），ｓ（τ），Ｑ（τ））を提示した画面イメージ４００Ｂである。ユーザインタフェース部１０４に出力される画面イメージ４００Ｂは、物体種別推定結果４１０と、警告表示４２０と、軌跡４３０及び音源分離結果４４０の４つの領域を含む。
【０１５７】
前述のとおり、推定結果は、推定した複数物体軌跡Ｔｒａｊ（τ）と、それらの軌跡の各物体の源信号を推定した音源分離後信号ｓ（τ）と、各物体の物体種別推定結果Ｑ（τ）からなる。物体種別推定結果Ｑ（τ）は、各物体ｇの各種別ｋの物体である確率Ｑ（τ，ｇ，ｋ）からなる。
【０１５８】
複数の物体が同時に存在する場合であっても、各物体に対して推定した軌跡と、その軌跡の物体の源信号を推定した音源分離結果と、物体種別推定結果を同時に提示することで、ユーザは、特定の種別の物体が、特定の軌跡で移動しているという情報を容易に知ることができる。ユーザはこの情報に基づいて、次の行動を決定することができるという効果がある。
【０１５９】
また、本実施例３は、複数の物体の軌跡仮説の検証の中に、源信号の状態が時間変化するモデルを導入しているので、この例のように、源信号の状態が時間変化する複数の物体が交差するような場合でも、高精度な追跡が可能である。
【０１６０】
本実施例３は、以上の構成により、物体が複数存在し、源信号の状態が時間変化する場合であっても、信号特徴の状態の特定、および、軌跡の追跡が可能である。
【０１６１】
上記実施例１〜３の推定システム１００は、推定装置であってもよい。また、複数チャンネルＡＤ変換器１０２Ａと、中央演算装置１０３の間にネットワークなどの通信手段を含んでいても良い。
【０１６２】
なお、本発明は上記した実施例に限定されるものではなく、様々な変形例が含まれる。例えば、上記した実施例は本発明を分かりやすく説明するために詳細に説明したものであり、必ずしも説明した全ての構成を備えるものに限定されるものではない。また、ある実施例の構成の一部を他の実施例の構成に置き換えることが可能であり、また、ある実施例の構成に他の実施例の構成を加えることも可能である。また、各実施例の構成の一部について、他の構成の追加・削除・置換をすることが可能である。
【０１６３】
また、上記の各構成、機能、処理部、処理手段等は、それらの一部又は全部を、例えば集積回路で設計する等によりハードウェアで実現してもよい。また、上記の各構成、機能等は、プロセッサがそれぞれの機能を実現するプログラムを解釈し、実行することによりソフトウェアで実現してもよい。各機能を実現するプログラム、テーブル、ファイル等の情報は、メモリや、ハードディスク、ＳＳＤ（Ｓｏｌｉｄ  Ｓｔａｔｅ  Ｄｒｉｖｅ）等の記録装置、または、ＩＣカード、ＳＤカード、ＤＶＤ等の記録媒体に置くことができる。
【０１６４】
また、制御線や情報線は説明上必要と考えられるものを示しており、製品上必ずしも全ての制御線や情報線を示しているとは限らない。実際には殆ど全ての構成が相互に接続されていると考えてもよい。
【符号の説明】
【０１６５】
１０１＿１〜１０１＿ｎ  素子
１０２、１０２Ａ  ＡＤ変換器
２０１  源信号状態モデルデータベース
２０２、８０３、１２０２  状態推定部
２０３、８０４、１２０３  結果出力部
２０４  信号取得部
８０１  方向推定部
８０２、１２０１  軌跡推定部
【図１】
2016065764.tif 000019
【図２】
2016065764.tif 000020
【図３】
2016065764.tif 000021
【図４】
2016065764.tif 000022
【図７】
2016065764.tif 000023
【図８】
2016065764.tif 000024
【図９】
2016065764.tif 000025
【図１１】
2016065764.tif 000026
【図１２】
2016065764.tif 000027
【図１３】
2016065764.tif 000028
【図５Ａ】
2016065764.tif 000029
【図５Ｂ】
2016065764.tif 000030
【図６】
2016065764.tif 000031
【図１０】
2016065764.tif 000032
【図１４】
2016065764.tif 000033