3-3-4
波源拘束差分方程式に基づく音響信号の確率モデル化と複数音源定位アルゴリズム ∗☆鈴木惇 (東大院・情報理工), 亀岡弘和 (東大院・情報理工 NTT CS 研)

1

はじめに

波源定位は，レーダやソナーといった幅広い応用を有する。特に，小さいアレイで，移動する波源を瞬時に定位し追跡可能にすることは重要課題である。波源定位問題に対する従来法としては, Multiple SignalClassification (MUSIC) 法 [1]，Generalized CrossCorrelation methods with Phase Transform (GCCPHAT) 法 [2], 波源拘束偏微分方程式に基づく手法[3, 4, 5] などがある。
MUSIC 法や GCC-PHAT 法は，センサアレイの受信信号間の自己相関関数や相互相関関数といったセンサ間の到来時間差に由来する統計量に基づく手法であるため，音源を高い精度で定位するためには大きいアレイサイズと十分長い観測時間幅を必要とする。
一方，波源拘束偏微分方程式に基づく手法は，各時刻ごとに成立する音響信号の時空間偏微分方程式を元に音源定位を行うもので, 理論的には瞬時の小領域観測のみで波源定位を行うことが可能だが，単一波源に対して成立する方程式を基にしているため，複数の音源を同時に定位することはできない。また，雑音が存在する場合など，観測音響信号が偏微分方程式から逸脱する場合に脆弱であるという欠点を有する。
そこで本研究は，小領域・瞬時観測による複数音源の波源定位を可能にする手法の提案を目的として，波源拘束偏微分方程式の周波数領域表現を基にした音響信号およびその空間差分の確率分布を構築することで複数の音源と加法雑音が存在する場合の混合信号およびその空間差分の確率分布を記述し，それに基づき複数の音源の波源定位を行うための ExpectationMaximization (EM) アルゴリズムを導出する。また、シミュレーションで提案手法の有用性を検証する。

2手法

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

Fig. 2 観測信号 f の空間微分を取得するためのアレイ幾何の例。例えば x 方向の f の空間微分は f1,t −f2,tで近似できる。となる。また，f (r, t) の時間微分は()
1|r − r 0 |f˙(r, t) =ġ t −|r − r 0 |cとなるので，g を消去して，{}
11∇f (r, t) =f (r, t) + f˙(r, t) nRc(3)(4)のように，観測信号とその時間・空間微分のみを含む方程式を立てられる。ただし，距離 R = |r − r 0 | である。この式を波源拘束式と呼ぶ [3, 4, 5]。以上のように波源拘束式は，任意の音源信号波形で成り立つ，音源の位置と空間の場の一意な関係を記述する偏微分方程式である。

2.2 波源拘束式に基づく音響信号の確率モデル化
観測信号の空間微分を空間差分で近似する場合を考える。簡単のため，以下では図 2 のようなマイクロホンアレイの配置を仮定する。このとき 7 本のマイクロホンを用いて各時刻 tl で，基準点における信f −f
号 f0,l およびその各方向の空間差分 fx,l = 1,l2D 2,l ,f −f
f −ffy,l = 3,l2D 4,l , fz,l = 5,l2D 6,l を得ることができる。
ただし，l は離散時刻のインデックスを表す。基準点における観測信号の時間微分を時間差分で近似する
ことにすると，式 (4) は)(1 f0,l − f0,l−11f0,l +nξ(5)fξ,l =RcTGenerative modeling of acoustic signals for multiple sound source localization based on spatio-temporal
partial diﬀerence equation. by SUZUKI, Atsushi (Graduate School of Information Science and Technology, The University of Tokyo), KAMEOKA, Hirokazu (Graduate School of Information Science andTechnology, The University of Tokyo NTT Communication Science Laboratories, NTT Corporation)

日本音響学会講演論文集- 615 -

2016年3月

と表せる。ただし，ξ = x, y, z で，nx , ny , nz はそれ
ぞれ n の x, y, z 方向の成分，T はサンプリング周期である。式 (5) の左辺を右辺に移項すると()
11 f0,l − f0,l−1−fξ,l +f0,l +nξ = 0(6)RcTで与えられる。ただし blockdiag{} は括弧内の行列を対角に並べたブロック対角行列。式 (9), (10) より，ϵは平均が 0，分散共分散行列が Σ}{Σ = blockdiag Σ0 , . . . ΣL/2{ 2}(17)222Σm = diag σx,m, σy,m, σz,m, σ0,mが得られる。ここで，f0,l , fx,l , fy,l , fz,l を窓関数で窓掛けして取得された信号とする。切り出し区間の両端点の影響を無視すると，式 (6) は周波数領域で()2πm1
1 − e−j L−Fξ,m ++nξ F0,m = 0(7)RcTの複素正規分布と表される。ただし，F0,m , Fx,m , Fy,m , Fz,m は f0,m ,fx,m , fy,m , fz,m の離散 Fourier 変換であり，m は離散周波数インデックスである。
式 (7) の右辺は雑音の存在や差分近似に伴う誤差により実際には必ずしも厳密に 0 にはならない。そこ
で，式 (6) の右辺を)(2πm1 − e−j L
1+nξ F0,m = ϵξ,m (8)−Fξ,m +RcTのように誤差変数 ϵx,m , ϵy,m , ϵz,m に置き換え，これらを平均が 0 で互いに独立な正規確率変数（複素正規分布に従う確率変数）
2ϵξ,m ∼ NC (ϵξ,m ; 0, σξ,m)(9)と仮定する。また，観測点における観測信号の各周波2の正規確率変数とす
数成分を平均が 0，分散が σ0,mる。これは，F0,m = ϵ0,m ,2ϵ0,m ∼ NC (ϵ0,m ; 0, σ0,m)(10)と仮定することに相当する。ここで，[]Tf m = Fx,m Fy,m Fz,m F0,m ∈ C4 ,[]Tϵm = ϵx,m ϵy,m ϵz,m ϵ0,m ∈ C4(11)(12)とし，これらを縦に連結して[]TTf = fT∈ C2L+4.
.
.
f0
L/2[]TT.
.
.
ϵϵ = ϵT∈ C2L+40
L/2(13)(15)(18)に従う。det A(θ) = (−1)L/2+1 ̸= 0（A(θ) は正則）であるので，f は

f = A(θ)−1 ϵ(19)と表され，式 (18) より，

f ∼ NC (f ; 0, A(θ)−1 ΣA(θ)−H )(20)が言える。従って，観測信号およびその空間差分が与えられた下での最尤音源位置 θ̂ はθ̂ = argmax log p(f |θ)(21)

θp(f |θ) = NC (f ; 0, A(θ)−1 ΣA(θ)−H )(22)により得られる。

2.3 複数音源の定位アルゴリズム
以上の f の確率モデル化により，音源が複数個存在する場合，および雑音が存在する場合の観測信号の確率分布が導ける。音源インデックスを k とし，音源 k に由来する観測信号の成分，音源位置パラメー
タをそれぞれ f (k) , θ(k) とする。また，f (k) の周波数(k)m の成分エネルギーを σ0,m 2 とする。式 (23) より，
f (k) ∼ NC (f (k) ; 0, A(θ(k) )−1 Σ(k) A(θ(k) )−H ) (23)となる。また，加法雑音をυ とし，観測信号を y =∑ (k)(1)f+υとする。
f,.
. . , f (K) , υ が互いに独立kであれば，観測信号 y は y ∼ p(y|θ)

p(y|θ))(∑=NC y; 0,A(θ(k) )−1 Σ(k) A(θ(k) )−H + Γ=∏m(14)と表記すると，式 (8) はA(θ)f = ϵϵ ∼ NC (ϵ; 0, Σ)(kNC y m ; 0,∑(k) −H
C m (θ(k) )−1 Σ(k)) + Γmm C m (θ)k(24)に従う。ただし，Γ は υ の分散共分散行列である。
以上より，複数の音源と雑音が存在する場合の各音源の最尤音源位置 θ̂ = {θ̂(1) , . . . , θ̂(K) } は，観測信号TT
y = (y T0 , . . . , y L/2 ) が与えられた下での形で書ける。ただし，θ = {R, n} であり，A(θ) は{}θ̂ = argmax p(y|θ)(25)A(θ) = blockdiag C 0 (θ), . . . , C L/2 (θ)θ
−j 2πmを解くことにより得られる。
−1( R1 + 1−e cT L )nx
2πmy を不完全データ，x = (f (1) T , . . . , f (K) T , υ T )T を

1−e−j L1
4×4−1(R +cT 2πm )ny  ∈ C
C m (θ) = 
完全データと見なすと，以上の最尤推定問題に対し
−j−1 ( R1 + 1−e cT L )nz 
Expectation-Maximization (EM) アルゴリズムを適1用できる。完全データ対数尤度 log p(x|θ) は(16)log p(x|θ)= − log | det πΛ| − xH Λ−1 x(26)日本音響学会講演論文集- 616 -

2016年3月


O と書かれる。ただし，A(θ(1) )−1 Σ(1) A(θ(1) )−H{}..
(k) 2(k)(k) 2(k) 2(k) 2.
Σ=diagσ,σ,σ,σ(35)Λ=mx,my,mz,m0,mA(θ(K) )−1 Σ(K) A(θ(K) )−H(k)(k)OΓ であり，Φm は Φ の 4 × 4 のブロック対角成分。
M ステップでは Q(θ, θ′ ) ができるだけ大きくなる(27)(k)ように n(k) , R(k) , σ0,m 2 , Γ を更新する。Q(θ, θ′ ) をで与えられるので，y が与えられた下での log p(x|θ)(k)最大にする n(k) , R(k) , σ0,m 2 , Γ の同時最適解を解析の x に関する条件付き期待値（Q 関数）は的に求めることは難しいが，座標勾配法によりそれQ(θ, θ′ )ぞれの変数に関して Q(θ, θ′ ) が最大となるように反復′=E[log p(x|θ)|y, θ ](28)更新することで p(y|θ) を局所最大化することができ′−1Hる（EM アルゴリズムでは，M ステップで補助関数= − log | det πΛ| − tr(Λ E[xx |y, θ ])(29)が単調に増大することが保証されていれば収束性は∑θ=−log |Σ(k) | − log |Γ| − tr(Γ−1 E[vv H |y, θ ′ ])保証される）。以下に M ステップの更新方法を示す。
k1) n(k) の更新式∑(k) (k)H′(k) H (k)−1(k)n(k) は単位ベクトルなので，∥n(k) ∥22 = 1 の下で−tr(A(θ ) ΣA(θ )E[f f|y, θ ])tr(A(θ(k) )H Σ(k) −1 A(θ(k) )Φ(k) ) ができるだけ小さくkなるように n(k) を更新する。この制約つき最適化問
xで与えられる。ただし，= は x に関係する項のみに
題は，目的関数ついての等号を意味する。この関数が増大するようL(n(k) ) = tr(A(θ(k) )H Σ(k)−1 A(θ(k) )Φ(k) )(36)に θ を更新するステップ（M ステップ）と，更新した(k)(k)′′∑θ を θ′ に代入し，E[f f H |y, θ ] と E[vv H |y, θ ]−1=tr(C m (θ(k) )H Σ(k)C m (θ(k) )Φ(k)mm )を計算するステップ（E ステップ）を繰り返すことで
mp(y|θ) を局所的に最大にする θ を求められる。
(37)
完全データ x と不完全データ y の関係は][を n(k)T n(k) − 1 = 0 の条件のもと最小化する問題と(30)
y = Hx,
H= I· · I} I| ·{z等価である。この問題は二次制約付き二次最適化問K題であり、Lagrange 未定乗数法などで解ける。と書けるので，E[x|y, θ ′ ] と E[xxH |y, θ ′ ] はそれぞれ
2) 音源距離 R(k) の更新式
ρ(k) = R1(k) とする。上記と同様，C m (θ(k) ) は
E[x|y, θ ′ ] = ΛH T (HΛH T )−1 y
E[xxH |y, θ ′ ] = Λ − ΛH T (HΛH T )−1 HΛ+ E[x|y, θ ′ ]E[x|y, θ ′ ]H
 (1)
Φ∗

.
..
=

Φ(K)
∗

 (31)
Φ(K+1)

で与えられる。以上より，以下のアルゴリズムを得る。

1. θ を初期設定する。
2. （E ステップ）θ を θ ′ に代入し，式 (31) によりΦ(1) , . . . , Φ(K) , Φ(K+1) を計算する。
3. （M ステップ）下式により θ を更新する。
θ ← argmax Q(θ, θ′ )
θ∑Q(θ, θ′ ) = −log |Σ(k) | − log |Γ|(32)(38)
C m (θ(k) ) = ρ(k) V (k) + P (k)
m


(k)
0 0 0 nx
(k) (39)V (k) = 0 0 0 ny(k) 0 0 0 nz 0 0 00
2πm(k)
1−e−j L−1nx
cT 2πm(k) 
1−e−j L(k)−1ny
 (40)
cT 2πmPm = (k) 
1−e−j L
−1nz cT1と，ρ(k) に依存する項としない項に分解できるので，
Q(θ, θ′ ) の ρ(k) に関する偏微分を 0 と置くことで，∑−1 (k) (k)tr(ℜ(V (k)T Σ(k)P m Φm ))m

ρ(k) = − m∑−tr(A(θ(k) )H Σ(k)−1 A(θ(k) )Φ(k) )k− tr(Γ−1 Φ(K+1) )(33)

2.4 M ステップ更新式
A(θ(k) ) は C 1 (θ(k) ), . . . , C L (θ(k) ) を対角に並べたブロック対角行列なので，(k) H(k)−1(k)(k)

tr(A(θ ) ΣA(θ )Φ )∑(k) H (k) −1
tr(C m (θ ) Σm C m (θ(k) )Φ(k)
=m )

を得る。
(k)
3) σ0,m 2 の更新式
−1
上記と同様，Σ(k)を
m−1
Σ(k)= G(k)mm +1(k)σ0,m 2E 4,4(42)(k)と，σ0,m 2 に依存する項としない項に分解できるので，Q(θ, θ′ ) の σ0,m 2 に関する偏微分を 0 と置くことで，(k)

(34)

m

日本音響学会講演論文集

(41)

mk

∑

−1 (k) (k)tr(V (k)T Σ(k)V Φm )
m- 617 -(k)σ0,m 2 = [Φ(k)m ]4,4(43)2016年3月

Fig. 3Fig. 5 音源位置とマイク位置。赤バツ印が音源位置(別々に実験)。青マル印がマイク位置。

フレーム幅が L = 64 の場合の定位精度。

EM アルゴリズムは，望ましくない局所解に陥ることがあり、それを避けるために 360 点のランダムな方向から最尤なものを初期解する手法 (図中の “Proposed360”) についても実験を行った。初期解を吟味して与えた場合，提案手法は比較手法と比べ，高精度な定位が可能な場合があり，MUSIC 法や GCC-PHATが使えない短時間測定からの定位問題で有用であることが分かる。
Fig. 44

フレーム幅が L = 256 の場合の定位精度。

を得る。
4) 雑音分散共分散行列 Γ の更新式
雑音の分散共分散行列を}{2W L/2
Γ = blockdiag ν02 W 0 , . . . , νL/2

本研究は，波源拘束偏微分方程式の確率モデル化に基づき小領域・瞬時観測による複数音源の定位を可能にする手法を提案し，数値実験で有用性を示した。

(44)のように，正規化分散共分散行列モデル W m と周波22数 k の成分エネルギー νmの積で表し，νmを推定すべき変数とする。ここで W m は空間無相関モデルや，マイクの 2 点間の空間相関係数がηd (ω) =sin ωdc
ωdc(45)で与えられる拡散音場モデルなどから導かれる定数行列。Q(θ, θ′ ) の ν 2 に関する偏微分を 0 と置くことで，2νm=1(K+1)
tr(W −1)m Φm
4(46)を得る。

3数値実験

提案手法の性能を検証するため，残響環境下の 1 音源定位の数値シミュレーションを行った。今回，x 方向，y 方向の 2 次元モデルを使用した。部屋サイズは6m×10m×8m とし，中心に 0.1m 間隔でマイクを 5つ配置した。壁面の反射係数は 0.01, 0.5, 0.8（残響の影響の大きさに相当）の 3 条件とした。サンプリン
グ周波数は 16000, フレーム幅は L = 64, 256 とした。
音源の配置は図 5 の通り。それぞれの条件下で 40 種類の音源を発させ、到来角度を推定した際の誤差の二乗和平均の平方根 (deg) を評価した。通常の EM アルゴリズムによる手法が図中の “Proposed” である。
単一音源のみの存在を仮定した尤度関数（式 (22)）を最大にする方法（従来法に相当・図中の “OneSrc”）と比較している。図 3, 4 に実験結果を示す。

日本音響学会講演論文集

結論

謝辞 本研究は JSPS 科研費 26730100 の助成を受けて行われた。

参考文献[1] R. O. Schmidt, “Multiple emitter location andsignal parameter estimation,” IEEE Transactions on Antennas and Propagation, vol. 34,
no. 3, pp. 276–280, 1986.
[2] C. H. Knapp et al., “The generalized correlation method for estimation of time delay,” IEEE
Trans. ASSP, vol. 24, no. 4, pp. 320–327, 1976.
[3] 藤田悠哉, 小野順貴, 安藤繁, “有限時間窓と離散フーリエ変換の利用を可能にする音源定位の高速厳密解法とその実験” 日本音響学会 2006 年秋季
研究発表会講演論文集, 3-1-3, pp. 483–484, Sep.
2006.
[4] S. Ando, N. Ono, T. Nara, “Direct algebraic
method for sound source localization with finest
resolution both in time and frequency,” in Proc.
ICSV14, Jul. 2007.
[5] 小山翔一, 栗原徹, 安藤繁, “偏微分方程式の空間荷重積分による瞬時音源定位,” 日本音響学会 2008 年
秋季研究発表会講演論文集, 2-8-20, pp. 679–682,
Sep. 2008.

- 618 -2016年3月