(57)【要約】【課題】実環境に近い条件下であっても、音源数を適切に推定する。【解決手段】音源数推定装置Ｂは、Ｎ個の音源からの信号が混合された混合信号をM個のマイクロホンで観測した観測信号から、各観測信号の時間周波数成分からなる観測信号ベクトルyτωに対応する特徴ベクトルxτωを抽出する。そして、音源数推定装置Ｂは、特徴ベクトルxτωを所定の確率モデルにあてはめ、各音源の尤度の時系列が周波数ビン間で同期しているほど高い評価値を与える評価関数を用いて、確率モデルのモデルパラメータを推定し、該モデルパラメータを用いて、観測信号ベクトルyτωが、Ｎ個の音源よりも多く設定された各クラスタに属する条件付き確率である事後確率を計算する。その後、音源数推定装置Ｂは、事後確率を用いて、各クラスタの事後確率の総和を算出し、各クラスタの事後確率の総和に基づいて、音源数を推定する。【選択図】図３
2016045225.tif 000002
【特許請求の範囲】
【請求項１】
kをクラスタのインデックスとし、τを時間フレームのインデックスとし、ωを角周波数とし、
Ｎ個の音源からの信号が混合された混合信号をM個のマイクロホンで観測した観測信号から、各観測信号の時間周波数成分からなる観測信号ベクトルyτωに対応する特徴ベクトルxτωを抽出する特徴抽出部と、
前記特徴ベクトルxτωを所定の確率モデルにあてはめ、各音源の尤度の時系列が周波数ビン間で同期しているほど高い評価値を与える評価関数を用いて、前記確率モデルのモデルパラメータを推定し、該モデルパラメータを用いて、前記観測信号ベクトルyτωが、前記Ｎ個の音源よりも多く設定された各クラスタに属する条件付き確率である事後確率を計算するモデル推定部と、
前記事後確率が有意な値をとるクラスタの個数を、音源数として推定する音源数推定部と、
を含み、
前記確率モデルは、各音源に関する特徴ベクトルxτωの分布の重み付き和で表される混合モデルであり、
前記確率モデルの混合重みは、時間フレームτに依存し、角周波数ωに依存しない重みであり、
前記確率モデルのモデルパラメータは、前記混合重みと、各音源に関する前記特徴ベクトルxτωの分布のパラメータである
ことを特徴とする音源数推定装置。
【請求項２】
請求項１に記載の音源数推定装置であって、
前記音源数推定部は、各クラスタの事後確率の総和を算出し、該各クラスタの事後確率の総和を２つにクラスタリングし、総和の大きい方のクラスタに属するクラスタの数を音源数として推定することを特徴とする音源数推定装置。
【請求項３】
請求項１に記載の音源数推定装置であって、
前記音源数推定部は、各クラスタの事後確率の総和を算出し、該各クラスタの事後確率の総和が所定の閾値以上であるかを判定し、前記所定の閾値以上である事後確率の総和に対応するクラスタの数を音源数として推定することを特徴とする音源数推定装置。
【請求項４】
請求項１から３のいずれか１つに記載の音源数推定装置であって、
クラスタkに関する前記特徴ベクトルxτωの分布は、平均方向をakωとし、密度パラメータをκkωとするワトソン分布であり、
クラスタkに関する前記特徴ベクトルxτωの分布のパラメータは、前記平均方向akωと前記密度パラメータκkωであることを特徴とする音源数推定装置。
【請求項５】
請求項１から４のいずれか１つに記載の音源数推定装置であって、
前記混合重みの事前分布は、クラスタkに依存しないハイパーパラメータφを各混合重みの指数とする前記混合重みについてのディリクレ分布であることを特徴とする音源数推定装置。
【請求項６】
請求項１から５のいずれか１つに記載の音源数推定装置であって、
前記モデル推定部は、クラスタkに関する前記特徴ベクトルxτωの分布とクラスタkの時間フレームτにおける混合重みとの積に基づいて、特徴ベクトルxτωが与えられたもとで、xτωに対応する前記観測信号ベクトルyτωがクラスタkに属する条件付き確率を計算する事後確率計算部と、
前記条件付き確率とクラスタkに依存しないハイパーパラメータφとに基づいて、前記混合重みを更新する混合重み更新手段と、
前記条件付き確率と前記特徴ベクトルxτωとに基づいて、クラスタkに対する相関行列Rkωを計算する相関行列更新手段と、
前記相関行列Rkωの正規化された主成分ベクトルを新たな値として前記平均方向akωを更新する平均方向更新手段と、
前記相関行列Rkωの最大固有値に基づいて、前記密度パラメータκkωを更新する密度パラメータ更新手段と、
周波数ビンごとに、前記評価関数が最大になるように、前記平均方向akωと前記密度パラメータκkωとをそれぞれ音源間で並べ替えるパーミュテーション解決手段と、
を含むことを特徴とする音源数推定装置。
【請求項７】
請求項６に記載の音源数推定装置であって、
γkτωを前記条件付き確率とし、αkτを前記混合重みとし、d(τ,ω)を前記観測信号ベクトルyτωに寄与するクラスタ番号とし、Fを周波数ビンの数とし、・Hを・のエルミート転置とし、λkωを前記相関行列Rkωの最大固有値とし、
前記事後確率計算部は、次式により前記条件付き確率を計算し、
【数３３】
2016045225.tif 000035
前記混合重み更新手段は、次式により求めたα’kτを新たな値として前記混合重みを更新し、
【数３４】
2016045225.tif 000036
前記相関行列更新手段は、次式により求めたR’kωを新たな値として前記相関行列Rkωを更新し、
【数３５】
2016045225.tif 000037
前記密度パラメータ更新手段は、次式により求めたκ’kωを新たな値として前記密度パラメータκkωを更新する
【数３６】
2016045225.tif 000038
ことを特徴とする音源数推定装置。
【請求項８】
請求項１から７のいずれか１つに記載の音源数推定装置であって、
ｎを音源の番号とし、＾Ｎを前記音源数推定部により推定された音源数とし、
前記音源数推定部は、前記事後確率が有意な値をとるクラスタの番号k(1)、k(2)、…、k(^N)を出力し、
前記事後確率が有意な値をとるクラスタの番号k(1)、k(2)、…、k(^N)と、前記各クラスタに属する条件付き確率である事後確率とを用いて、前記音源数推定部により推定された音源数に対応するマスクmｎτωを求めるマスク作成部と、
前記観測信号ベクトルyτωから前記マスクmｎτωを用いて時間周波数領域の分離音を計算する分離音作成部と、
を含むことを特徴とする音源数推定装置。
【請求項９】
kをクラスタのインデックスとし、τを時間フレームのインデックスとし、ωを角周波数とし、
Ｎ個の音源からの信号が混合された混合信号をM個のマイクロホンで観測した観測信号から、各観測信号の時間周波数成分からなる観測信号ベクトルyτωに対応する特徴ベクトルxτωを抽出するステップと、
前記特徴ベクトルxτωを所定の確率モデルにあてはめ、各音源の尤度の時系列が周波数ビン間で同期しているほど高い評価値を与える評価関数を用いて、前記確率モデルのモデルパラメータを推定し、前記モデルパラメータを用いて、前記観測信号ベクトルyτωが、前記Ｎ個の音源よりも多く設定された各クラスタに属する条件付き確率である事後確率を計算するステップと、
前記事後確率が有意な値をとるクラスタの個数を、音源数として推定するステップと、
を含み、
前記確率モデルは、各音源に関する特徴ベクトルxτωの分布の重み付き和で表される混合モデルであり、
前記確率モデルの混合重みは、時間フレームτに依存し、角周波数ωに依存しない重みであり、
前記確率モデルのモデルパラメータは、前記混合重みと、各音源に関する前記特徴ベクトルxτωの分布のパラメータである
ことを特徴とする音源数推定方法。
【請求項１０】
請求項１から８のいずれか１つに記載の音源数推定装置としてコンピュータを機能させるための音源数推定プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、音源数推定装置、音源数推定方法および音源数推定プログラムに関する。
【背景技術】
【０００２】
従来、音源分離技術として、クラスタリングに基づく音源分離を行う技術が知られている。例えば、音源数が既知であることを前提として、音源数と同数のクラスタを用いてクラスタリングを行う（例えば、非特許文献１参照）。１音源の特徴量は、各周波数内で同じ特徴量空間に集中しやすいため、このようなクラスタリングにより音源分離を行うことができる。
【０００３】
各音源の特徴量空間の形状は、周波数に依存するが、１音源信号は、全周波数で同時に立ちあがりやすく、且つ、立ち下がりやすい。すなわち、音源アクティビティは周波数間で同期する。非特許文献１に記載のクラスタリングでは、この音源アクティビティの同期性を適切にモデル化し、全体最適化に組み込むことにより、音源位置特徴量と音源アクティビティの同時クラスタリングを実現している。
【０００４】
ところが、非特許文献１に記載のクラスタリングでは、音源数が既知であることを前提としており、音源数が未知である場合には、音源数と同数のクラスタを用いてクラスタリングを行うことができない。一方、音源数を推定する技術が知られている。例えば、音源数を推定する技術として、「音源数≦マイクロホン数、かつ残響なし」であることを前提とし、音源数を推定する方法がある（例えば、非特許文献２を参照）。
【先行技術文献】
【非特許文献】
【０００５】
【非特許文献１】伊藤信貴、荒木章子、中谷智広“時変混合重みに基づくパーミュテーション問題のないクラスタリングベース音源分離”、社団法人電子情報通信学会、信学技報,  VOL.113,  NO.27,  2013年5月
【非特許文献２】MATI  WAX,THOMAS  KAILATH,“Detection  of  Signals  by  Information  Theoretic  Criteria”,  IEEE  TRANSACTIONS  ON  ACOUSTICS,  SPEECH,  AND  SIGNAL  PROCESSING,  VOL.  ASSP-33,  NO.2,  APRIL  1985
【発明の概要】
【発明が解決しようとする課題】
【０００６】
しかしながら、上記した従来の音源数を推定する技術では、「音源数≦マイクロホン数、かつ残響なし」であることを前提としており、実環境に近い条件下では、音源数を適切に推定することができないという課題があった。つまり、従来の音源数を推定する技術は、「音源数≦マイクロホン数、かつ残響なし」という理想的な条件でしか成立せず、例えば、「音源数＞マイクロホン数、または残響あり」というような状況である場合には、音源数を適切に推定することができない。
【課題を解決するための手段】
【０００７】
上述した課題を解決し、目的を達成するために、本発明の音源数推定装置は、kをクラスタのインデックスとし、τを時間フレームのインデックスとし、ωを角周波数とし、Ｎ個の音源からの信号が混合された混合信号をM個のマイクロホンで観測した観測信号から、各観測信号の時間周波数成分からなる観測信号ベクトルyτωに対応する特徴ベクトルxτωを抽出する特徴抽出部と、前記特徴ベクトルxτωを所定の確率モデルにあてはめ、各音源の尤度の時系列が周波数ビン間で同期しているほど高い評価値を与える評価関数を用いて、前記確率モデルのモデルパラメータを推定し、該モデルパラメータを用いて、前記観測信号ベクトルyτωが、前記Ｎ個の音源よりも多く設定された各クラスタに属する条件付き確率である事後確率を計算するモデル推定部と、前記事後確率が有意な値をとるクラスタの個数を、音源数として推定する音源数推定部と、を含み、前記確率モデルは、各音源に関する特徴ベクトルxτωの分布の重み付き和で表される混合モデルであり、前記確率モデルの混合重みは、時間フレームτに依存し、角周波数ωに依存しない重みであり、前記確率モデルのモデルパラメータは、前記混合重みと、各音源に関する前記特徴ベクトルxτωの分布のパラメータであることを特徴とする。
【発明の効果】
【０００８】
本発明によれば、実環境に近い条件下であっても、音源数を適切に推定することが可能であるという効果を奏する。
【図面の簡単な説明】
【０００９】
【図１】図１は、第一実施形態に係るモデル推定装置の機能構成を例示する図である。
【図２】図２は、第一実施形態に係るモデル推定装置の処理フローを例示する図である。
【図３】図３は、第二実施形態に係る音源数推定装置の機能構成を例示する図である。
【図４】図４は、第二実施形態に係る音源数推定装置の処理フローを例示する図である。
【図５】図５は、実際の音源数を超えるクラスタ数を設定した場合の事後確率の例を示す図である。
【図６】図６は、実験結果を示す図である。
【図７】図７は、音源とクラスタ数とを同じに設定した場合における特徴量空間にプロットされる音源位置特徴量を示す図である。
【図８】図８は、音源よりもクラスタ数を多く設定した場合における特徴量空間にプロットされる音源位置特徴量の例を示す図である。
【図９】図９は、音源よりもクラスタ数を多く設定した場合における特徴量空間にプロットされる音源位置特徴量の例を示す図である。
【図１０】図１０は、第三実施形態に係る音源数推定装置の機能構成を例示する図である。
【図１１】図１１は、第三実施形態に係る音源数推定装置の処理フローを例示する図である。
【図１２】図１２は、音源数推定方法を実行するコンピュータを示す図である。
【発明を実施するための形態】
【００１０】
以下、この発明の実施の形態について詳細に説明する。なお、図面中において同じ機能を有する構成部には同じ番号を付し、重複説明を省略する。
【００１１】
［発明のポイント］
この発明の音源数推定装置における技術のポイントは、音源数が未知である条件下において、音源数を推定することが可能である点である。詳細は後述するが、音源数よりも多いクラスタ数を設定して非特許文献１に記載のクラスタリングを行い、各クラスタの事後確率を計算し、事後確率が有意な値をとるクラスタの個数を、音源数として推定するものである。
【００１２】
［第一実施形態］
この発明の第一実施形態は、複数の音源からの信号を複数個のマイクロホンで観測し、各音源に関するモデルパラメータを推定するモデル推定装置である。
【００１３】
図１を参照して、第一実施形態のモデル推定装置Ａの機能構成例を説明する。モデル推定装置Ａは、周波数領域変換部１、特徴抽出部２及びモデル推定部３を含む。モデル推定部３は、事後確率計算部３１、パラメータ更新部３２及びパラメータ保持部３３を含む。パラメータ更新部３３は、混合重み更新手段３２１、相関行列更新手段３２２、平均方向更新手段３２３、密度パラメータ更新手段３２４及びパーミュテーション解決手段３２５を含む。
【００１４】
図２を参照して、モデル推定装置Ａの動作例を手続きの順に従って説明する。周波数領域変換部１へM個のマイクロホンにより観測された時間領域の混合信号~ytが入力される。時間領域の混合信号~ytは式(１)で定義される。
【数１】
2016045225.tif 000003
【００１５】
ここで、tは時間インデックスを表し、・T（上付き文字のT）はベクトル・の転置を表し、~ymtはm（1≦m≦M）番目のマイクロホンで観測された時間領域の混合信号を表す。
【００１６】
周波数領域変換部１は入力された時間領域の混合信号~ytから短時間フーリエ変換などにより時間周波数領域の観測信号ベクトルyτωを生成し出力する（ステップＳ１）。時間周波数領域の観測信号ベクトルyτωは式(２)で定義される。
【数２】
2016045225.tif 000004
【００１７】
ここで、τは時間フレームのインデックスを表し、ωは角周波数を表し、ymτωは混合信号~ymtの時間周波数領域での表現である。
【００１８】
特徴抽出部２は、周波数領域変換部１の出力する時間周波数領域の観測信号ベクトルyτωを入力とし、特徴ベクトルxτωを計算し出力する（ステップＳ２）。特徴ベクトルxτωの計算は、時間周波数領域の観測信号ベクトルyτωを正規化することにより行ってもよいし、時間周波数領域の観測信号ベクトルyτωを白色化した後に正規化することによって行ってもよいし、時間周波数領域の観測信号ベクトルyτωを正規化した後に白色化して再度正規化することにより行ってもよい。例えば、時間周波数領域の観測信号ベクトルyτωを正規化することにより特徴ベクトルxτωを計算する場合は式(３)により計算すればよい。
【数３】
2016045225.tif 000005
【００１９】
また、例えば、時間周波数領域の観測信号ベクトルyτωを白色化した後に正規化する場合には、以下のように特徴ベクトルxτωを計算すればよい。まず、時間周波数領域の観測信号ベクトルyτωを用いて、式(４)により時間周波数領域の観測信号ベクトルyτωの標本相関行列Rωｙを計算する。
【数４】
2016045225.tif 000006
【００２０】
ここで、Tはフレームの個数であり、・H（上付き文字のH）はエルミート（Hermite）転置である。
【００２１】
次に、標本相関行列Rωｙの固有値と固有ベクトルを計算する。計算した標本相関行列Rωｙの固有値を、大きい順に並べたものをσω1,σω2,…,σωMと表す。したがって、式(５)の関係が成り立つ。
【数５】
2016045225.tif 000007
【００２２】
ここで、標本相関行列Rωｙはエルミート行列であるから、固有値σω1,σω2,…,σωMはすべて実数であることに注意する。また、固有値σω1,σω2,…,σωMに対応し、正規直交系をなす標本相関行列Rωｙの固有ベクトルをuω1,uω2,…,uωMで表す。ここで、標本相関行列Rωｙはエルミート行列であるから、このような固有ベクトルが存在することに注意する。
【００２３】
次に、行列Σωを式(６)により求め、行列Uωを式(７)により求める。
【数６】
2016045225.tif 000008
【００２４】
次に、行列Uω及び行列Σωを用いて、時間周波数領域の観測信号ベクトルyτωを白色化したベクトルy’τωを式(８)により計算する。
【数７】
2016045225.tif 000009
【００２５】
最後に、次式のようにベクトルy’τωをそのノルムで正規化することにより、特徴ベクトルxτωを計算する。
【数８】
2016045225.tif 000010
【００２６】
モデル推定部３は、特徴ベクトルxτωを特徴ベクトルの分布を表す確率モデルに当てはめ、確率モデルを評価する所定の評価関数を用いて、クラスタリングに適した確率モデルのモデルパラメータを計算する。そして、モデル推定部３は、該モデルパラメータを用いて、観測信号ベクトルyτωが、Ｎ個の音源よりも多く設定された各クラスタに属する条件付き確率である事後確率を計算する。
【００２７】
音源の位置が固定の場合、各周波数ビンにおいて、特徴ベクトルxτωは理想的にはそれぞれの音源ごとに固有の値をとる。ただし、実際には、雑音・残響の影響、モデル化誤差などに起因する変動が存在するため、特徴ベクトルxτωは音源ごとにある値を中心にクラスタを形成する。そのため、この発明では、クラスタkに関する特徴ベクトルxτωの分布を、例えば、以下のようにワトソン（Watson）分布でモデル化する。
【００２８】
【数９】
2016045225.tif 000011
【００２９】
ここで、akωはクラスタkに関する特徴ベクトルの分布の中心を表し、平均方向（mean  orientation）と呼ばれ、κkωはクラスタkに関する特徴ベクトルの分布の広がりの小ささを表し、密度パラメータ（concentration  parameter）と呼ばれる。M(a,b,x)はクンマー（Kummer）関数である。クンマー関数についての詳細は「S.  Sra  and  D.  Karp,  “The  multivariate  Watson  distribution:  maximum-likelihood  estimation  and  other  aspects”,  arXiv:  1104.4422v2,  2012.（参考文献１）」を参照されたい。ここで、特徴ベクトルの分布が周波数ビンごとに定義されることに注意する。
【００３０】
音声をはじめとする多くの音源信号は、「音源信号の時間周波数変換の振幅値の時系列{|sｎτω|}τが、周波数ビン間で類似する」という共通振幅変調の性質をもつ（例えば、「G.  J.  Brown,  “Computational  Auditory  Scene  Analysis:  A  Representational  Approach”,  Ph.D.  thesis,  University  of  Sheffield,  1992.」を参照）。この発明では、この共通振幅変調の性質を、パーミュテーション問題を回避するための手掛かりとして利用できることに着目した。混合信号の各時間周波数成分に寄与する音源信号は高々一つであるという仮定（WDO（W-Disjoint  Orthogonality）性の仮定）にもとづき、この共通振幅変調の性質を、クラスタリングの枠組みにおいて利用しやすい表現で言い換えると、「観測信号に寄与するクラスタ番号の時系列{d(τ,ω)}τは、周波数ビン間で類似する」と言える。この発明では、この周波数ビン間での{d(τ,ω)}τの類似性を、「d(τ,ω)の事前分布P(d(τ,ω)=k)が、フレームτに依存（時変）し、周波数ビン（角周波数ω）にはよらない（周波数非依存）」とモデル化する。このような各音源信号に対する、周波数ビン間での振幅変調の共通性を利用することにより、パーミュテーションを引き起こさずに、クラスタリングできる。この事前確率をαkτにより表す。なお、αkτはΣk=1Kαkτ=1を満たす。ここでは、音源数は未知であることを前提とし、モデル推定装置Ａでは、音源数を超えるクラスタ数が設定されているものとする。
【００３１】
この事前確率は、１個の時間フレームごとに変化すると仮定してもよいし、複数の時間フレームからなるブロックごとに変化すると仮定してもよい。事前確率が１個の時間フレームごとに変化すると仮定する場合、任意のクラスタkと任意の時間フレームτに対して、αkτは独立変数であり、推定すべきパラメータである。
【００３２】
一方、事前確率が数個の時間フレームからなるブロックごとに変化すると仮定する場合、Bをブロックの総数とし、ブロック番号をb=1,2,…,Bとし、Jを各ブロック内における時間フレームの総数とし、各ブロック内における時間フレームの番号をj=1,2,…,Jとすると、τ=(b−1)×J+jと表せ、αk,(b−1)×J+j（j=1,2,…,J）は等しくなるから、推定すべきパラメータである混合重みは~αkb=αk,(b−1)×J+1により定義される~αkbである。以下では、特に断りのない限り、事前確率が１個の時間フレームごとに変化すると仮定する場合について説明する。
【００３３】
以上より、特徴ベクトルxτωの尤度関数は、式(１１)で表す混合モデルで与えられる。
【数１０】
2016045225.tif 000012
【００３４】
ここで、Θは、式(１２)に示すパラメータ集合である。
【数１１】
2016045225.tif 000013
【００３５】
ここで、{αkτ}kτは式(１３)により定義される。
【数１２】
2016045225.tif 000014
【００３６】
他の同様の記法もこれにならって定義される。以降では、αkτを混合重みと呼ぶ。混合重みαkτの事前分布として式(１４)に示すディリクレ（Dirichlet）分布を用いる。
【数１３】
2016045225.tif 000015
【００３７】
ここで、Γはガンマ関数であり、φはハイパーパラメータと呼ばれる。φの値を十分大きく定めることにより、混合重みαkτの変動を抑えることができる。φの値を微調整する必要はないが、例えば、φ=1,10,100,1000などの値を用いることができる。
【００３８】
混合重みαkτ以外のパラメータについては一様な事前分布を仮定する。したがって、p(Θ)=Πτp({αkτ}k)である。
【００３９】
モデル推定部３では、特徴ベクトルxτωを以上のようにモデル化された確率モデルに当てはめ、確率モデルを評価する所定の評価関数を用いて、事後確率及びクラスタリングに適したパラメータ集合Θを求める。
【００４０】
以下、モデル推定部３の各部の処理を詳細に説明する。モデル推定部３は、図１に示すとおり、事後確率計算部３１、パラメータ更新部３２及びパラメータ保持部３３を含む。モデル推定部３での処理に先立ち、パラメータ集合Θの初期値をパラメータ保持部３３に用意しておく（ステップＳ０）。この初期値は、例えば、αkτ=1/K、κkω=20とし、akωは{xτω}τωから無作為に選ぶことにより設定することができる。
【００４１】
事後確率計算部３１は、パラメータ保持部３３に記憶されたパラメータ集合Θから事後確率γkτω、すなわち特徴ベクトルxτωが与えられたもとでd(τ,ω)=kとなる条件付き確率を式(１５)により計算する（ステップＳ３１）。
【数１４】
2016045225.tif 000016
【００４２】
パラメータ更新部３２は、図１に示すとおり、混合重み更新手段３２１、相関行列更新手段３２２、平均方向更新手段３２３、密度パラメータ更新手段３２４及びパーミュテーション解決手段３２５を含み、現在のパラメータ集合Θを更新して新たなパラメータ集合Θ’を生成する（ステップＳ３２）。
【００４３】
混合重み更新手段３２１は、事後確率γkτωを用いて、式(１６)を計算することにより、混合重みαkτを新しい値α’kτに更新する。
【００４４】
【数１５】
2016045225.tif 000017
【００４５】
ここで、Fは周波数ビンの個数を表す。φ=1のとき、α’kτは全周波数ビンにわたる事後確率γkτωの平均値となることがわかる。φの増加とともに、α’kτは定数1/Kに近づく。
【００４６】
相関行列更新手段３２２は、特徴ベクトルxτωと事後確率γkτωを用いて、式(１７)を計算することにより、各クラスタkに対する相関行列Rkωを新しい値R’kωに更新する。
【数１６】
2016045225.tif 000018
【００４７】
平均方向更新手段３２３は、相関行列Rkωの正規化された主成分ベクトルとして、平均方向akωを新しい値a’kωに更新する。
【００４８】
密度パラメータ更新手段３２４は、相関行列Rkωの最大固有値λkωを用いて、密度パラメータкkωを式(１８)により新しい値к’kωに更新する。
【数１７】
2016045225.tif 000019
【００４９】
パーミュテーション解決手段３２５は、式(１９)〜(２１)に示すように、各周波数ビンにおいて、平均方向a’kωと密度パラメータк’kωを、事後確率p(Θ’|{xτω}τω)が最大になるように音源間で置換し、パーミュテーションを解決する（ステップＳ３２５）。
【００５０】
【数１８】
2016045225.tif 000020
【００５１】
なお、以上では、混合重みが１個の時間フレームごとに変化する場合の処理について説明したが、混合重みが複数の時間フレームからなるブロックごとに変化する場合は、混合重み更新手段３２１における混合重みαkτの更新式(１６)において、分子の事後確率γkτωの時間フレーム内の和を事後確率γkτωのブロック内の和に置き換え、分母のFをF×Jで置き換えればよい。一方、相関行列更新手段３２２、平均方向更新手段３２３、密度パラメータ更新手段３２４及びパーミュテーション解決手段３２５においては、混合重みが１個の時間フレームごとに変化する場合の処理と同一の処理を行えばよい。
【００５２】
以下、パラメータ更新部３２における各更新式の導出根拠を説明する。パラメータ更新はEM（Expectation-Maximization）アルゴリズムを導入して、それに基づき行う。なお、{d(τ,ω)}τωは、EMアルゴリズムにおける隠れ変数として扱う。
【００５３】
まず、MAP（Maximum  a  posteriori）推定のためのコスト関数L(Θ)は、式(２２)〜(２４)により与えられる。
【数１９】
2016045225.tif 000021
【００５４】
ここで、{xτω}τωは互いに独立であると仮定し、Θに依存しない定数項を無視した。この目的関数を式(２５)に示す制約条件のもとで最大化する。
【数２０】
2016045225.tif 000022
【００５５】
目的関数L(Θ)は、パーミュテーション問題がない場合に大きい値を取るため、L(Θ)の最大化によりパーミュテーション問題が回避できる。実際、式(２４)の第一項から分かるように、目的関数L(Θ)が大きくなるのは、混合重みαkτが大きい値をとるk、τに対し、クラスタkに対する尤度（もっともらしさ）p(xτω|d(τ,ω)=k,akω,κkω)が大きい場合である。したがって、L(Θ)の最大化により、クラスタkに対する尤度の時系列{p(xτω|d(τ,ω)=k,akω,κkω)}τが周波数ビン間で同期する。このことと、上述の「観測信号に寄与する音源インデックスの時系列{d(τ,ω)}τは、周波数ビン間で類似する」という性質を考え合わせると、L(Θ)はパーミュテーション問題がない場合に大きい値を取ることがわかる。EMアルゴリズムで用いる評価関数（Q関数）は式(２６)(２７)により与えられる。
【００５６】
【数２１】
2016045225.tif 000023
【００５７】
更新後のパラメータ集合Θ'は次式により定義される。
【数２２】
2016045225.tif 000024
【００５８】
Q関数を式(２５)の制約のもとで最大にするものとして導かれる。すなわち、混合重みαkτの新たな値α'kτを求める式(１６)は、ラグランジュ（Lagrange）の未定乗数法によって、式(２８)(２９)により導出される。
【数２３】
2016045225.tif 000025
ここで、μはラグランジュの未定乗数である。
【００５９】
平均方向の更新式は、
【数２４】
2016045225.tif 000026
を解くことで得られる。したがって、クーラン・フィッシャー（Courant-Fischer）の定理より、Rkωの最大固有値に対応する固有ベクトル（主成分ベクトル）によって、平均方向を更新すればよい。
【００６０】
また、密度パラメータの更新式(１８)については、まず∂Q/∂κkω=0より式(３１)を得る。
【数２５】
2016045225.tif 000027
【００６１】
ここで、
【数２６】
2016045225.tif 000028
であり、λkωは相関行列Rkωの最大固有値である。上式は、近似的に次の式(３２)のように解くことができる。
【００６２】
【数２７】
2016045225.tif 000029
【００６３】
パラメータ保持部３３は、パラメータ更新部３２での更新処理により得られたパラメータ集合Θ’を記憶する（ステップＳ３３）。また、事後確率計算部３１での次回の処理の際には、記憶したパラメータ集合Θ’をパラメータ集合として提供する。なお、第一の実施形態のモデル推定装置を用いて推定した、パラメータ集合Θを用いることで、音源定位、音源数が未知の条件下での音源分離、音源数と分離音の同時推定等を行うことが可能である。
【００６４】
ステップＳ３１からステップＳ３３までの処理は、事前に設定した最大反復回数max_iterに達するまで、またはパラメータ更新部３２における各パラメータの更新による変動幅が収束判定の閾値Δよりも小さくなるまで、反復して行う（ステップＳ９１）。最大反復回数max_iter及び閾値Δの具体的な値は、例えば、max_iter=100、Δ=10−10とすることができる。
【００６５】
ステップＳ９１において、モデル推定部３における処理が最大反復回数max_iterに達した場合、または各パラメータの更新による変動幅が閾値Δよりも小さくなった場合、モデル推定部３は反復終了後の事後確率γokτωを出力する。
【００６６】
［第二実施形態］
この発明の第二実施形態は、第一実施形態のモデル推定装置Ａを用いて音源数推定装置として構成した実施形態である。
【００６７】
図３を参照して、第二実施形態の音源数推定装置Ｂの機能構成例を説明する。音源数推定装置Ｂは、第一実施形態のモデル推定装置Ａの各部に加えて、音源数推定部４を含む。
【００６８】
音源数推定部４は、事後確率が有意な値をとるクラスタの個数を、音源数として推定する。具体的には、音源数推定部４は、事後確率計算部３１によって計算された各クラスタの事後確率の入力を受け付け、各クラスタの事後確率を用いて、各クラスタの事後確率の総和を算出する。例えば、音源数推定部４は、下記式（３３）を用いて、各クラスタの事後確率の総和を算出する。なお、各クラスタの事後確率の総和を算出することに限定させるものではなく、例えば、特定の周波数範囲に限った事後確率の部分和を算出してもよい。
【００６９】
【数２８】
2016045225.tif 000030
【００７０】
そして、音源数推定部４は、各クラスタの事後確率の総和を２つにクラスタリングし、総和の大きい方のクラスタに属するクラスタの数を音源数として推定する。
【００７１】
図４を参照して、音源数推定装置Ｂの動作例を手続きの順に従って説明する。ステップＳ０からステップＳ９１までの処理は第一実施形態のモデル推定装置Ａの動作例と同様であるので詳細な説明は省略する。
【００７２】
音源数推定部４は、各クラスタの事後確率を用いて、各クラスタの事後確率の総和を算出する（ステップＳ４１）。具体的には、音源数推定部４は、事後確率計算部３１によって計算された各クラスタの事後確率の入力を受け付け、各クラスタの事後確率を用いて、各クラスタの事後確率の総和を算出する。
【００７３】
そして、音源数推定部４は、各クラスタの事後確率の総和を２つにクラスタリングする（ステップＳ４２）。例えば、音源数推定部４は、各クラスタの事後確率の総和に対して、クラスタ数２のk-meansクラスタリングを適用してクラスタリングする。
【００７４】
続いて、音源数推定部４は、事後確率の総和の大きい方のクラスタに属するクラスタ数を音源数と推定する（ステップＳ４３）。例えば、音源数推定部４は、セントロイドがより大きいクラスタに属するクラスタの個数を、音源数として推定する。
【００７５】
ここで、音源数推定部４は、上記で例示したように、クラスタ数が「２」のk-meansクラスタリングを適用してクラスタリングする場合に限定されるものではなく、より簡易な処理として、所定の閾値を用いて、各クラスタの事後確率の総和に対して閾値処理を行ってもよい。具体的には、音源数推定部４は、各クラスタの事後確率の総和が所定の閾値以上であるかを判定し、所定の閾値以上である事後確率の総和に対応するクラスタの数を音源数として推定してもよい。クラスタ数が「２」のk-meansクラスタリングを適用してクラスタリングする方法は、所定の閾値以上である事後確率の総和に対応するクラスタの数を音源数として推定する方法と比べ、残響時間などの条件の変化に対して、より頑健であると期待される。
【００７６】
第一実施形態の説明では、音源数が未知であることを前提として、音源数を超えるクラスタ数を設定しているものとして説明した。このように、音源数を超えるクラスタを設定して、段落００４０〜００６５において説明した、パラメータの更新処理および事後確率の計算処理を行うことで、音源数と同数のクラスタのみ有意な事後確率をもち、残りのクラスタは小さい事後確率をもつこととなる。ここで、図５において、音源数が「３」であり、クラスタ数「４」の場合に、第二実施形態に記載のクラスタリング方法で得られた事後確率の例を示す。図５について、横軸が時間であり、縦軸が周波数であり、輝度の大きい点ほど、事後確率が大きいことを示している。
【００７７】
つまり、実際の音源数と同数のクラスタのみ、音源アクティビティが有意な値をもっていることがわかる。図５の例を用いて説明すると、「cluster1」、「cluster2」、「cluster4」の３つのクラスタは有意な事後確率をもち、「cluster3」は小さい事後確率であり、有意な事後確率をもっていないといえる。
【００７８】
このため、上述したように、音源数推定部４は、クラスタ数２のk-meansクラスタリングを適用してクラスタリングした場合には、「cluster1」、「cluster2」および「cluster4」の３つのクラスタの事後確率の総和と、「cluster3」のみの１つのクラスタの事後確率の総和とでそれぞれクラスタリングされる。そして、事後確率の総和の大きい方のクラスタに属するクラスタ数（要素数）が、「cluster1」、「cluster2」および「cluster4」の３つであるから、音源数が「３」と推定される。
【００７９】
［実験結果］
ここで、音源数推定装置Ｂによる音源数推定処理を実施した場合の実験結果について説明する。
【００８０】
図６に、音源数推定実験における正解率を示す。図６に示すように、φ＝６００とした場合、全ての条件において、１００％の正解率が得られた。また、φを１→６００と増加させることにより、正解率が向上する傾向があった。一方、φをさらに６００→１０００と増加させると、音源数Ｎ＝３、残響時間３７０ｍｓの条件での正解率が、１００％→８８％と減少した。
【００８１】
以上のように、従来は音源数の推定が困難であるような条件、例えば、「音源数>マイクロホン数、または残響あり」のような条件下でも、音源数を適切に推定することが可能である。
【００８２】
続いて、本発明の原理について説明する。まず、非特許文献１に記載のクラスタリングに基づく音源分離法について説明する。非特許文献１に記載のクラスタリングに基づく音源分離法は、音源数が既知であることを前提としており、音源数と同数のクラスタを用いてクラスタリングを行うものである。図７を参照して、上記の非特許文献１の技術について説明する。ここで、図７中の×印は、各時間周波数点における音源位置特徴量（例：マイクロホン間の位相差、時間差、振幅比）を表す。図７に示すように、同一音源の×印がそれぞれ円で囲まれており、１音源の特徴量は、各周波数内で、同じ特徴量空間に集中しやすい。各音源の特徴量空間の形状は、周波数に依存するが、１音源信号は、全周波数で同時に立ち上がりやすく、且つ、立ち下がりやすい。すなわち、音源アクティビティは周波数間で同期する。このようなクラスタリングでは、この音源アクティビティの同期性を適切にモデル化し、全体最適化に組み込むことにより、音源位置特徴量と音源アクティビティの同時クラスタリングを実現している。
【００８３】
上述のように、非特許文献１の技術では、音源数と同数のクラスタを用いてクラスタリングを行っていた。これを、音源数よりも多いクラスタを用いる構成に変更することは、これまでの技術常識に鑑みれば、容易に想到できるものではなかった。なぜならば、クラスタリングに基づく音源分離において、音源数よりも多いクラスタを用いると、通常、図８に示すように、１つの音源が複数のクラスタに分裂してしまうため、各クラスタが音源に対応しないからである。なお、図８は、音源数２、クラスタ数４の例である。
【００８４】
これに対して、本発明に至る研究の過程で、次の全く新しい知見が得られた。すなわち、非特許文献１のクラスタリングにおいて、音源数よりも多いクラスタを用いる構成に変更すると、上記の技術常識に反して、１音源が複数のクラスタに分かれるようなことは起こらず、１音源は１つのクラスタにまとまる。つまり、先行技術のクラスタリングでは、上述のように、音源アクティビティの同期性を適切にモデル化しているため、図８のようなことは起こらず、図９のように１音源は１つのクラスタにまとまる。なお、図９は、音源数２、クラスタ数３の例である。図９に示すように、音源１および音源２がそれぞれ一つのクラスタにまとまっている。また、残り一つのクラスタは小さい事後確率をもつクラスタとなる。
【００８５】
本発明では、以上の知見に基づき、非特許文献１のクラスタリングを、音源数よりも多いクラスタを用いる構成に変更する。また、本発明では、クラスタリング結果に基づき音源数を決定する際に、各クラスタの事後確率の総和を大小２クラスタにクラスタリングし、大きい方のクラスタの要素数として音源数を決定する。これにより、固定の閾値を用いる場合と比べて、より頑健に音源数を推定することができる。本発明の方法により、音源数が未知でも音源分離を行ったり、音源数を推定したりすることが可能になる。また、上記の実験結果からも分かるように、本発明の方法を用いて、極めて良好に音源数を推定できることが確認された。
【００８６】
また、本発明の音源数推定法では、「音源数＞マイクロホン数、または残響あり」の場合にも、音源数を適切に推定することができる。実際、実験により、「音源数＞マイクロホン数」の場合や、残響時間が比較的長い場合にも、本発明の音源数推定法により、きわめて良好に音源数を推定できることが示された。
【００８７】
上述してきたように、音源数推定装置Ｂは、kをクラスタのインデックスとし、τを時間フレームのインデックスとし、ωを角周波数とし、Ｎ個の音源からの信号が混合された混合信号をM個のマイクロホンで観測した観測信号から、各観測信号の時間周波数成分からなる観測信号ベクトルyτωに対応する特徴ベクトルxτωを抽出する。そして、音源数推定装置Ｂは、特徴ベクトルxτωを所定の確率モデルにあてはめ、各音源の尤度の時系列が周波数ビン間で同期しているほど高い評価値を与える評価関数を用いて、確率モデルのモデルパラメータを推定し、該モデルパラメータを用いて、観測信号ベクトルyτωが、Ｎ個の音源よりも多く設定された各クラスタに属する条件付き確率である事後確率を計算する。その後、音源数推定装置Ｂは、事後確率を用いて、各クラスタの事後確率の総和を算出し、各クラスタの事後確率の総和に基づいて、音源数を推定する。これにより、音源数推定装置Ｂでは、実環境に近い条件下であっても、音源数を適切に推定することが可能である。
【００８８】
また、音源数推定装置Ｂは、各クラスタの事後確率の総和を２つにクラスタリングし、総和の大きい方のクラスタに属するクラスタの数を音源数として推定する。音源数推定装置Ｂは、閾値を用いる場合と比べて、より頑健に音源数を推定することができる。
【００８９】
また、音源数推定装置Ｂは、各クラスタの事後確率の総和が所定の閾値以上であるかを判定し、所定の閾値以上である事後確率の総和に対応するクラスタの数を音源数として推定してもよい。この場合には、上記の２つにクラスタリングを行う場合と比べて、より簡易に音源数を推定することができる。
【００９０】
［第三実施形態］
この発明の第三実施形態は、第二実施形態に係る音源数推定装置Ｂの構成に音源分離部５および時間領域変換部６を追加した音源数推定装置Ｃとして構成した実施形態である。
【００９１】
図１０を参照して、第三実施形態の音源数推定装置Ｃの機能構成例を説明する。音源数推定装置Ｃは、第二実施形態の音源数推定装置Ｂの各部に加えて、音源分離部５及び時間領域変換部６を含む。音源分離部５は、マスク作成部５１及び分離音作成部５２を含む。
【００９２】
図１１を参照して、音源数推定装置Ｃの動作例を手続きの順に従って説明する。ステップＳ０からステップＳ４３までの処理は第二実施形態の音源数推定装置Ｂの動作例と同様であるので詳細な説明は省略する。音源数推定部４は、前記事後確率が有意な値をとるクラスタの番号k(1)、k(2)、…、k(^N)を出力し、音源分離部５での処理に供する。
【００９３】
音源分離部５は、音源数推定部４の出力する前記事後確率が有意な値をとるクラスタの番号k(1)、k(2)、…、k(^N)と、周波数領域変換部１の出力する混合音の時間周波数変換yτωと、事後確率計算部３１の出力する反復終了後の事後確率γokτωとを用いて、分離音の時間周波数変換^sｎτωを推定する（nは音源の番号）。
【００９４】
マスク作成部５１は、音源数推定部４の出力する前記事後確率が有意な値をとるクラスタの番号k(1)、k(2)、…、k(^N)と、事後確率計算部３１の出力する前記反復終了後の事後確率γokτωとを用いて、音源数推定部４により推定された音源数に対応するマスクmｎτωを求める（ステップＳ５１）。まず、マスク作成部５１は、反復終了後の事後確率γokτωのうち、有意な値である事後確率γｓｎτωをγｓｎτω＝γoｋ（ｎ）τωにより計算する。ここで、n=1、2、…、^Nであり、γの上付き文字のＳは、Significant（有意）を表す。次に、マスク作成部５１は、反復終了後の有意な値である事後確率γsｎτωを用いて、式(３４)により、アクティブな音源の番号dｓ(τ,ω)の推定値^dｓ(τ,ω)を計算する。
【数２９】
2016045225.tif 000031
【００９５】
次に、マスク作成部５１は、マスクmｎτωを式(３５)により計算する。
【数３０】
2016045225.tif 000032
【００９６】
なお、マスク作成部５１は、マスクmｎτωを式(３６)により求めてもよい。
【数３１】
2016045225.tif 000033
【００９７】
分離音作成部５２は、式(３７)により、マスクmｎτωを混合音の時間周波数変換y1τωに乗算し、分離音の時間周波数変換^sｎτωを計算する。これにより、周波数領域の観測信号を音源ごとに分離する（ステップＳ５２）。
【数３２】
2016045225.tif 000034
【００９８】
時間領域変換部６は、音源ｎごとに、時間周波数領域の分離信号^sｎτωを時間領域の分離信号~^sｎｔに変換して出力する（ステップＳ６）。
【００９９】
上述してきたように、音源数推定装置Ｃは、音源数が未知である場合であっても、音源数を推定した後に、音源分離技術を実現することができる。また、音源数推定装置Ｃは、パーミュテーション問題を生じず、二段階の処理を必要としない音源分離技術を実現することができる。これにより、例えば、音源位置などが時間的に変化する時変の環境での音声強調のためのオンライン音源分離を容易に実現することが可能となる。
【０１００】
［システム構成等］
また、図示した各装置の各構成要素は機能概念的なものであり、必ずしも物理的に図示の如く構成されていることを要しない。すなわち、各装置の分散・統合の具体的形態は図示のものに限られず、その全部または一部を、各種の負荷や使用状況などに応じて、任意の単位で機能的または物理的に分散・統合して構成することができる。例えば、モデル推定部３と音源数推定部４を統合してもよい。さらに、各装置にて行なわれる各処理機能は、その全部または任意の一部が、ＣＰＵおよび当該ＣＰＵにて解析実行されるプログラムにて実現され、あるいは、ワイヤードロジックによるハードウェアとして実現され得る。
【０１０１】
また、本実施例において説明した各処理のうち、自動的におこなわれるものとして説明した処理の全部または一部を手動的におこなうこともでき、あるいは、手動的におこなわれるものとして説明した処理の全部または一部を公知の方法で自動的におこなうこともできる。この他、上記文書中や図面中で示した処理手順、制御手順、具体的名称、各種のデータやパラメータを含む情報については、特記する場合を除いて任意に変更することができる。
【０１０２】
［プログラム］
また、上記実施形態に係る音源数推定装置Ｂ、Ｃが実行する処理をコンピュータが実行可能な言語で記述したプログラムを作成することもできる。この場合、コンピュータがプログラムを実行することにより、上記実施形態と同様の効果を得ることができる。さらに、かかるプログラムをコンピュータ読み取り可能な記録媒体に記録して、この記録媒体に記録されたプログラムをコンピュータに読み込ませて実行することにより上記実施形態と同様の処理を実現してもよい。以下に、音源数推定装置Ｂ、Ｃと同様の機能を実現する音源数推定プログラムを実行するコンピュータの一例を説明する。
【０１０３】
図１２は、音源数推定プログラムを実行するコンピュータを示す図である。図１２に示すように、コンピュータ１０００は、例えば、メモリ１０１０と、ＣＰＵ１０２０と、ハードディスクドライブインタフェース１０３０と、ディスクドライブインタフェース１０４０と、シリアルポートインタフェース１０５０と、ビデオアダプタ１０６０と、ネットワークインタフェース１０７０とを有する。これらの各部は、バス１０８０によって接続される。
【０１０４】
メモリ１０１０は、ＲＯＭ（Read  Only  Memory）１０１１およびＲＡＭ（Random  Access  Memory）１０１２を含む。ＲＯＭ１０１１は、例えば、ＢＩＯＳ（Basic  Input  Output  System）等のブートプログラムを記憶する。ハードディスクドライブインタフェース１０３０は、ハードディスクドライブ１０９０に接続される。ディスクドライブインタフェース１０４０は、ディスクドライブ１０４１に接続される。ディスクドライブ１０４１には、例えば、磁気ディスクや光ディスク等の着脱可能な記憶媒体が挿入される。シリアルポートインタフェース１０５０には、例えば、マウス１１１０およびキーボード１１２０が接続される。ビデオアダプタ１０６０には、例えば、ディスプレイ１１３０が接続される。
【０１０５】
ここで、図１２に示すように、ハードディスクドライブ１０９０は、例えば、ＯＳ１０９１、アプリケーションプログラム１０９２、プログラムモジュール１０９３およびプログラムデータ１０９４を記憶する。上記実施形態で説明した各データは、例えばハードディスクドライブ１０９０やメモリ１０１０に記憶される。
【０１０６】
また、音源数推定プログラムは、例えば、コンピュータ１０００によって実行される指令が記述されたプログラムモジュールとして、ハードディスクドライブ１０９０に記憶される。具体的には、上記実施形態で説明した音源数推定装置Ｂ、Ｃが実行する各処理が記述されたプログラムモジュールが、ハードディスクドライブ１０９０に記憶される。
【０１０７】
また、音源数推定プログラムによる情報処理に用いられるデータは、プログラムデータとして、例えば、ハードディスクドライブ１０９０に記憶される。そして、ＣＰＵ１０２０が、ハードディスクドライブ１０９０に記憶されたプログラムモジュール１０９３やプログラムデータ１０９４を必要に応じてＲＡＭ１０１２に読み出して、上述した各手順を実行する。
【０１０８】
なお、音源数推定プログラムに係るプログラムモジュール１０９３やプログラムデータ１０９４は、ハードディスクドライブ１０９０に記憶される場合に限られず、例えば、着脱可能な記憶媒体に記憶されて、ディスクドライブ１０４１等を介してＣＰＵ１０２０によって読み出されてもよい。あるいは、分散データ処理プログラムに係るプログラムモジュール１０９３やプログラムデータ１０９４は、ＬＡＮ（Local  Area  Network）やＷＡＮ（Wide  Area  Network）等のネットワークを介して接続された他のコンピュータに記憶され、ネットワークインタフェース１０７０を介してＣＰＵ１０２０によって読み出されてもよい。
【符号の説明】
【０１０９】
Ａ  モデル推定装置
Ｂ、Ｃ  音源数推定装置
１  周波数領域変換部
２  特徴抽出部
３  モデル推定部
３１  事後確率計算部
３２  パラメータ更新部
３２１  混合重み更新手段
３２２  相関行列更新手段
３２３  平均方向更新手段
３２４  密度パラメータ更新手段
３２５  パーミュテーション解決手段
３３  パラメータ保持部
４  音源数推定部
５  音源分離部
５１  マスク作成部
５２  分離音作成部
６  時間領域変換部
【図１】
2016045225.tif 000039
【図２】
2016045225.tif 000040
【図３】
2016045225.tif 000041
【図４】
2016045225.tif 000042
【図５】
2016045225.tif 000043
【図６】
2016045225.tif 000044
【図７】
2016045225.tif 000045
【図８】
2016045225.tif 000046
【図９】
2016045225.tif 000047
【図１０】
2016045225.tif 000048
【図１１】
2016045225.tif 000049
【図１２】
2016045225.tif 000050