3-1-9
空間勾配スペクトルのベイズモデリングによる複数音源の瞬時定位 ∗◎植野夏樹（東大院情報理工，NTT CS 研）, 亀岡弘和（NTT CS 研）

1

はじめに波源定位は，レーダやソナーといった幅広い応用を有している。特に，小さいアレイで，移動する波源を瞬時に定位し追跡できるようにすることは重要課題である。波源定位問題に対する従来法としては, Multiple Signal Classification (MUSIC) 法 [1]，Generalized Cross-Correlation methods with Phase Transform (GCC-PHAT) 法 [2], 波源拘束偏微分方程式に基づく手法 [3, 4, 5] などがある。
MUSIC 法や GCC-PHAT 法は，音源に対し平面波を仮定し各音源のセンサ間での到来時間差を定位の手がかりとするため，一般にアレイサイズは大きい方が有利となる。また，いずれもセンサアレイの受信信号間の自己相関関数や相互相関関数といった，統計量に基づく手法であるため，音源を高い精度で定位するためには観測時間幅を十分長く取る必要がある。
このため，これらの手法は小さいアレイサイズと瞬時的な観測のみによる波源定位には必ずしも向いていない。一方，波源拘束偏微分方程式に基づく手法
は，各時刻ごとに成立する音響信号の時空間偏微分方
程式を元に音源定位を行うもので, 理論的には瞬時の小領域観測のみで波源定位を行うことが可能である。
ただし，この手法は単一波源に対して成立する方程式をベースとしているため，雑音や複数の点音源が存在する場合などのように観測音響信号が偏微分方程式から逸脱する場合に脆弱であるという欠点を有している。
そこで本研究では，本枠組においてランダム雑音や複数の点音源が存在する場合への拡張を可能にするため，波源拘束偏微分方程式を基にしたアレイ観測信号の確率モデルを導き，ベイズ推論に基づく波源定位アルゴリズムを導出する。また，数値シミュレーションによって提案手法の有効性を検証する。

2

音源定位問題の定式化

2.1 波源拘束偏微分方程式
観 測 点 の 基 準 と な る 位 置 ベ ク ト ル をT
r = (x, y, z) ∈ R3 とし，単一波源の位置ベクT
トルを r0 = (x0 , y0 , z0 ) ∈ R3 とする。波源の信号を g(t)，音速を c とし，単一点波源からの球面波伝播を仮定すると観測点における観測値は()
1∥r − r0 ∥f (r, t) =g t−(1)∥r − r0 ∥cと表される。ただし，√∥r − r0 ∥ = (x − x0 )2 +(y − y0 )2 +(z − z0 )2 (2)である。観測点から波源方向へ向かう単位ベクトルを n とすると，f (r, t) の空間微分は，()
1∥r − r0 ∥∇f (r, t) =gt−n(3)∥r − r0 ∥2c()1∥r − r0 ∥+ġ t −n (4)c∥r − r0 ∥c∗

Fig. 1点音源 r0 から観測点 r へ到来する球面波。となる。また，f (r, t) の時間微分は()
1∥r − r0 ∥f˙(r, t) =ġ t −∥r − r0 ∥cとなるので，g を消去して，}{
1 ˙
1f (r, t) + f (r, t) n
∇f (r, t) =Rc(5)(6)のように，観測信号とその時間・空間微分のみを含む方程式を立てることができる。ただし，R = ∥r − r0 ∥は観測点から波源までの距離である。この式を音源拘束式と呼ぶ [3, 4, 5]。以上のように音源拘束式は，任意の音源信号波形で成り立つ，音源の位置と空間の場の一意な関係を記述する偏微分方程式である。

2.2 波源拘束式に基づく音響信号の確率モデル化
観測信号の空間微分を空間差分で近似するため，以下では図 2 に示すマイクロホンアレーを仮定する。ただし，マイクロホンの配置は観測信号の空間微分を空間差分で近似できるものであれば良く，以下の理論は図 2 の配置に限らない。図 2 のマイクロホンアレーの場合，7 本のマイクロホンを用いて各時刻 tj で，基準点における信号 f0,j およびその各方向の空間差分
f −ff −ff −ffx,j = 1,j2D 2,j , fy,j = 3,j2D 4,j , fz,j = 5,j2D 6,j を得ることができる。ただし，j は標本時刻のインデックスを表す。このとき式 (6) は()
11fi,j =f0,j + f˙0,j ni(7)Rcと表せる。ただし，i = x, y, z で，nx , ny , nz はそれ
ぞれ n の x, y, z 方向の成分である。式 (7) の左辺を右辺に移項すると()
11 ˙−fi,j +f0,j + f0,j ni = 0(8)Rcが得られる。ここで，f0,j , fx,j , fy,j , fz,j を窓関数で窓掛けして取得された信号とする。切り出し区間の両端点の影響を無視できるものとすると，式 (8) は周波数領域で()1ωm−Fi,m ++jni F0,m = 0(9)Rc

Multiple source localization with Bayesian modeling of spatial gradient spectra. by UENO, Natsuki(University of Tokyo, NTT), KAMEOKA, Hirokazu (NTT)

日本音響学会講演論文集- 479 -2017年3月

に従う。det Cm (θ) = −1 ̸= 0（Cm (θ) は正則）であるので，fm はfm = Cm (θ)−1 ϵm(20)と表される。よって式 (19) より，fm の確率密度関数p(fm |θ) = NC (fm ; 0, Cm (θ)−1 Σm Cm (θ)−H ) (21)Fig. 2 観測信号 f の空間微分を取得するためのアレ
f −fイ幾何の例。例えば x 方向の f の空間微分は 1,j2D 2,jで近似できる。と表される。ただし，F0,m , Fx,m , Fy,m , Fz,m は f0,j ,fx,j , fy,j , fz,j の離散 Fourier 変換であり，m は離散周波数インデックスである。
式 (9) の右辺は差分近似に伴う誤差により実際には必ずしも厳密に 0 にはならない。そこで，式 (8) の右辺を()1ωm−Fi,m ++jni F0,m = ϵi,m(10)Rcのように誤差変数 ϵx,m , ϵy,m , ϵz,m に置き換え，これらを平均が 0 で互いに独立な正規確率変数（複素正規分布に従う確率変数）
2ϵi,m ∼ NC (ϵi,m ; 0, σm)(11)と仮定する。また，観測点における観測信号の各周波2数成分を平均が 0，分散が σ0,mの正規確率変数とする。これは，F0,m = ϵ0,m(12)2ϵ0,m ∼ NC (ϵ0,m ; 0, σ0,m)(13)と仮定することに相当する。
ここで，Fx,m , Fy,m , Fz,m , F0,m を並べたベクトルと ϵx,m , ϵy,m , ϵz,m , ϵ0,m を並べたベクトルを

fm = (Fx,m , Fy,m , Fz,m , F0,m ) ∈ C
Tϵm = (ϵx,m , ϵy,m , ϵz,m , ϵ0,m )T ∈ C44

θ̂ = argmax log p(f |θ)θ∏p(f |θ) =p(fm |θ)(15)(16)の形で書ける。ただし θ = {R, n} であり，Cm (θ) は
−1( R1 + j ωcm )nx−1( R1 + j ωcm )ny  (17)Cm (θ) = −1 ( R1 + j ωcm )nz 
1(18)の複素正規分布ϵm ∼ NC (ϵm ; 0, Σm )日本音響学会講演論文集(23)

mを解く最尤推定問題に帰着する。

2.3 混合音信号の確率モデル化
音声信号や楽音など実世界の音響信号の多くは時間周波数成分がスパースである。従って，複数の音源が同時に混在する場合でも，各時間周波数点では高々一つの音源のみが支配的であると仮定できる場合が多い。この音源の時間周波数成分のスパース性の仮定と以上の f の確率モデル化に基づき，音源が複数個存在する場合，および雑音が存在する場合の観測信号の確率分布を導くことができる。
信号の切り出しフレームの時刻のインデックスをl，音源インデックスを k = 0, . . . , K とし，k = 0 は雑音，k ̸= 0 は点音源に対応するものとする。また，
点音源 k の位置を θ(k) = {R(k) , n(k) } とする。ここで，雑音を含む全音源の時間周波数成分のスパース性を仮定し，周波数 m，時刻 l において zm,l 番目の音源のみが非零のパワーをもち，それ以外の音源のパワーを 0 とする。このとき，所与の zm,l の下での観測信号の時間周波数成分とその空間差分 ym,l ∈ C4（以後，観測信号）の条件付き確率密度関数は(k)p(ym,l |zm,l = k, ϑ) = NC (ym,l ; 0, Λm,l )(24)で与えられる。ここで，{(k)Cm (θ(k) )−1 Σm,l Cm (θ(k) )−H (k ̸= 0)(k)
Λm,l =Γm,l(k = 0)
(25)
(k)

(19)

(k)

222Σm,l = diag(σm, σm, σm, σ0,m,l 2 )(26)である。また，Γm,l は雑音の時間周波数成分の分散共分散行列で，周波数にのみ依存する正規化分散共分散行列モデル Wm と時刻にも依存する雑音のパワー
2νm,lの積で与えられる。式 (11), (13) より，ϵm は平均が 0，分散共分散行列 Σm が2222Σm = diag(σm, σm, σm, σ0,m)(22)(14)とする。式 (10) はCm (θ)fm = ϵmを得る。以上より，観測スペクトルおよびその空間差分 f = {fm }m が与えられた下で，単一音源を定位する問題は，

2Γm,l = νm,lWm(27)で表されるものとする。拡散雑音を仮定する場合のWm の設定方法については [6] を参照されたい。ϑ はす(k)2べての未知パラメータ {θ(k) }k , {σm}m , {σ0,m,l }k,m,l ,2{νm,l}m,l を表す。zm,l の事前確率を p(zm,l = k) =(k)αm,l とすると，観測信号 y = {ym,l }m,l の確率密度- 480 -

2017年3月

関数（ϑ の尤度関数）はp(y|ϑ) =K∏∑

2.5 変分推論アルゴリズム
以上の生成モデル化により，p(ym,l |zm,l = k, ϑ)p(zm,l = k)p(Y |N, ρ, Z, λ, ζ, γ)∞∏ ∏(k) 1[zm,l =k]=NC (ym,l ; 0, Λm,l )m,l k=0=K
∏∑

(k)

(k)

αm,l NC (ym,l ; 0, Λm,l )(28)

m,l k=0∞∏ ∏m,l k=0

p(Z|V ) =と書ける。以上より，複数の音源と雑音が存在する場
合の各音源の位置 θ̂ = {θ̂(1) , . . . , θ̂(K) } を推定する問
題は，観測信号 y が与えられた下でm,l k=01[zπk m,lp(V ) = Beta(v0 ; αb , βb )

=k]
∞
∏(43)
(44)

Beta(vk ; 1, βs )(45)k=1ϑ̂ = argmax log p(y|ϑ)(29)

ϑを解く最尤推定問題に帰着する。この最適化問題の大域解は解析的に解くことはできないが，Z = {zm,l }m,lを潜在変数とした Expectation-Maximization (EM)
アルゴリズムにより停留点を探索することができる。

2.4 ノンパラメトリックベイズモデリング
実環境では音源数が未知の場合が多い。上記の定式化では音源数 K を既知と仮定したが，音源数を仮定せずともモデルの複雑度を適応させながら実際の音源数に合わせて音源定位を行えるようにすることが望ましい。そこで，本章では上記の生成モデルを [7]と同様に Dirichlet 過程混合モデルに拡張する。k = 0を雑音，k = 1, . . . , ∞ を点音源のインデックスとし，zm,l を加算無限次元の離散分布 π0 , π1 , . . . , π∞ に従って生成される確率変数とする。
zm,l ∼ Categorical(π0 , π1 , . . . , π∞ )(30)

こ こ で ，点 (m, l) に お い て 雑 音 が 支 配 的 で あ る（zm,l = 0 となる）確率 π0 を，超パラメータ αb , βbのベータ分布に従って生成される変数

π 0 = v0v0 ∼ Beta(αb , βb )(31)(32)とし，点音源 k ̸= 0 が支配的である（zm,l = k となる）確率 πk を，棒折り過程 [8] に従って決まる変数

π k = vkk−1∏

(1 − vk )(33)

j=0vk ∼ Beta(1, βs )(34)とする。以上のプロセスで生成される π1 , π2 , . . . の期
待値は，k が大きいほど指数的に小さくなる傾向を持
つため，大きい k に対応した音源ほどアクティブになる確率が低くなる。よって，観測信号からパラメータを推論する際，必要最小限の音源インデックス数の混合モデルで観測信号を説明しようとする効果がもたらされる。以上の生成モデルにおいて，全未知変数
Θ は以下となる。

Θ ={N, ρ, Z, V, λ, ζ, γ}
N = {n(k)(35)}k(k)= 1/R(k))

(37)
(38)
(39)

2
λ = {λm } (λm = 1/σm)

(40)

(k)
ζ = {ζm,l }(41)

γ = {γm,l }日本音響学会講演論文集

p(N ) =p(ρ) =(k)(k)2(ζm,l = 1/σ0,m,l )2(γm,l = 1/νm,l)∞∏vMF(n(k) ; ξ (k) , κ(k) )(46)

N (ρ(k) ; µ(k) , η (k) )(47)

Gamma(λm ; αλ,m , βλ,m )

(48)

k=1∞
∏
k=1p(λ) =∏m,lp(ζ) =∞∏ ∏(k)(49)

Gamma(γm,l ; αγ,m,l , βγ,m,l )(50)

m,l k=1p(γ) =∏(k)(k)Gamma(ζm,l ; αζ,m,l , βζ,m,l )

m,lと仮定すると，観測信号 Y = {ym,l }m,l と未知変数
Θ の同時分布p(Θ, Y ) = p(Y |N, ρ, Z, λ, ζ, γ)p(Z|V )p(V )p(N )p(ρ)p(λ)p(ζ)p(γ)(51)

を具体的に記述することができる。ただし，vMF(·)
は von Mises-Fisher 分 布 ，N (·) は 実 正 規 分 布 ，Gamma(·) はガンマ分布を表す。
Θ の事後分布 p(Θ|Y ) = p(Θ, Y )/p(Y ) を解析的に得ることは難しいが，変分推論法に基づき近似分布を反復計算により得ることができる。変分推論は，∫q(Θ)dΘ = 1 を満たす非負の変関数 q(Θ) を，事後分
布 p(Θ|Y ) との間の Kullback-Leibler (KL) ダイバー
ジェンス
∫q(Θ)
dΘ (52)KL(q(Θ)∥p(Θ|Y )) = q(Θ) log
p(Θ|Y )が小さくなるように求める方法である。q(Θ) をq(Θ) = q(N )q(ρ)q(Z)q(V )q(λ)q(ζ)q(γ)(53)のように近似できると仮定すると，q(N ), q(ρ), q(Z),q(V ), q(λ), q(ζ), q(γ) について反復的に式 (52) を最小化することで p(Θ|Y ) の近似分布を得ることができる。また，q(Z) に関して，以下の打ち切り近似(36)

ρ = {ρ }k (ρZ = {zm,l }m,lV = {vk }k(k)と書ける。また，その他の事前分布 p(V )，p(N )，p(ρ)，p(λ)，p(ζ) 及び p(γ) をq(zm,l = K ∗ + 1) = · · · = q(zm,l = ∞) = 0(54)を行う。この近似はモデルの複雑度（音源数）を固定したことを意味するのではなく q の関数空間をある範囲に限定したことを意味する。従って，p(Θ|Y ) をできるだけ良く近似したければ K ∗ を大きくして q(Θ)
が取りうる範囲を広くすれば良い。

(42)- 481 -

2017年3月

(a) 残響時間 0.5 s，観測長2779 ms(b) 残響時間 0.5 s，観測長
1665 ms(c) 残響時間 0.2 s，観測長2779 ms(d) 残響時間 0.2 s，観測長
1665 ms

Fig. 3 音源位置とマイクロホン位置。赤バツ印が音源位置。青マル印がマイクロホンアレーの中心位置
導出は省略するが，式 (52) を各 q について最小化する更新式は以下の形となる。

q(N ) =∗K
∏vMF(n(k) ; ξ̄ (k) , κ̄(k) )(55)

N (ρ(k) ; µ̄(k) , η̄ (k) )(56)(z(57)

k=1q(ρ) =∗K∏

k=1q(Z) =∏m,lq(V ) =m,lφ̄m,l∗K∏)

Fig. 4(k)(k)Beta(vk ; β̄0 , β̄1 )(58)Gamma(λm ; ᾱλ,m , β̄λ,m )(59)

k=0q(λ) =∏m,l∗q(ζ) =

∏ K
∏(k)

(60)

Gamma(γm,l ; ᾱγ,m,l , β̄γ,m,l )(61)

m,l k=1q(γ) =

∏(k)

(k)

Gamma(ζm,l ; ᾱζ,m,l , β̄ζ,m,l )

m,l3

数値実験

提案手法の性能を検証するため，残響環境下の複数音源定位の数値シミュレーションを行った。今回，
x 方向，y 方向の 2 次元モデルを使用した。部屋サイズは 6 m×10 m×4 m とし，中心に 0.03 m 間隔で 7つのマイクロホンを図 2 のように配置した。壁面の反射係数は 0.7308 及び 0.4566（Sabine の残響公式による残響時間がそれぞれ 0.5 s 及び 0.2 s に相当）とした。マイクロホンのサンプリング周波数は 32 kHz，短時間 Fourier 変換のフレーム長は 64 ms（オーバーラップは 32 ms)，観測信号全体の長さは 2779 ms 及び 1665 ms とした。音源数は 3 個，音源の位置は部屋の中心からそれぞれ (1, 0, 0) m，(−0.5, 0.87, 0) m，(−0.5, −0.87, 0) m とした。音源信号は話速バリエーション型音声データベース（SRV-DB）を利用した。
提案法は変分推論に基づく手法（VBEM），正しい音源数（3 個）を仮定した EM アルゴリズムによる手法（EM1），誤った音源数（6 個）を仮定した EM アルゴリズムによる手法（EM2）を評価し，従来法である MUSIC 法との比較を行った。EM アルゴリズムに
22
おいては定常雑音を仮定し νm,l= ν̃mとした。どの手法においても，検出閾値をもとに音源の有無及び方向を推定し，そのうち真の音源方向との角度の誤差が±τ 以内のものがあれば正解，なければ脱落誤り，検出された音源方向のうちどの真の音源方向にも属さ

日本音響学会講演論文集

実験結果ないものを誤挿入として，F 尺度を算出した。各 τ について，検出閾値を変化させたときの最も高い F 尺度のプロットを図 4 に示す。多くの場合提案法，特に変分推論に基づく手法が他の手法より高精度な定位を達成していることが確認できる。

4結論

本研究は，波源拘束変微分方程式を基に音響信号の空間勾配スペクトルをベイズモデリングすることで小領域・瞬時観測による複数音源の定位を可能にする手法を提案し，数値実験で有用性を示した。
謝辞 本研究は JSPS 科研費 26730100 の助成を受けて行われた。

参考文献[1] R. O. Schmidt, IEEE Trans. AP, 34 (3), 276–280,
1986.
[2] C. H. Knapp et al., IEEE Trans. ASSP, 24(4),320–327, 1976.
[3] 藤田他, 音講論 (秋), 3-1-3, 483–484, 2006.
[4] S. Ando et al., Proc. ICSV14, 2007.
[5] 小山他, 音講論 (秋), 2-8-20, 679–682, 2008.
[6] 鈴木&亀岡, 音講論 (春), 3-3-4, 615–618, 2016.
[7] H. Kameoka et al., IEICE Trans. FE, E96-A(10),1928–1937, 2013.
[8] J. Sethuraman, Statistica Sinica, 4, 639–650,
1994.

- 482 -2017年3月