SignalAnalysis
株式会社シグナルアナリシス 技術資料 本文書の著作権は株式会社シグナルアナリシスに帰属します．

信号部分空間の推定と  株式会社

  シグナルアナリシス½

信号部分空間

センサーアレイによる磁場計測

 アルゴリズム

関原謙介

生体磁気イメージングでは信号計測を多数のセンサー，すなわちセンサーアレイによる同時計測で行うのが普通である．
列ベクトル

：



番目のセンサーの時刻 での出力を  で表せば，センサーアレイ全体の出力は       










 で表す．ここで， はセンサー数を表し，上付きの  は行列の転置を意味する． 次元空間における位置ベクトルを       と表せば，位置  における信号源電流 は  次元の列ベクトル             

½ 信号源  電流とは神経や筋肉の活動で生じた電気生理学的な電流を意味する．

で表される．r: r = [x, y, z] と表し，位置 r における信号源電流強度をスカラー s(r, t) で表し，向きは3 次元の単位ベクトル η(r) = [ηx (r), ηy (r), ηz (r)]T で表す．ここで，上付きの T は行列の転置を意味する．したがって，s(r, t) = s(r, t)η(r) の関係がある．
位置 r に単位強度の信号源が存在すると仮定する．この信号源ベクトルが x 方向を向いていた場合のxym 番目のセンサー出力を lm(r)，y 方向を向いていた場合の出力を lm(r)，および z 方向を向いていた場
zxyz
合の出力を lm(r) と表記し，m 番目の行が，行ベクトル [lm(r), lm(r), lm(r)] に等しい M ×3 の行列 L(r)を定義する．この L(r) はリードフィールド行列と呼ばれ，センサーアレイの位置 r における感度を表す．
位置 r において η(r) 方向を向いた単位強度の信号源によるセンサーアレイの出力は l(r) = L(r)η(r)で表される．この列ベクトル l(r) をリードフィールドベクトルと呼ぶ．l(r) はセンサーアレイの位置 rにおける η(r) 方向の感度を現すものである．
リードフィールドベクトルを用いると信号源分布 s(r, t) とセンサー出力 y(t) の関係は

∫y(t) =L(r)η(r)s(r, t) dr + ε = y S (t) + ε.

(3)Ωと表すことができる．ここで，右辺の積分は信号源が存在する可能性のある 3 次元の領域 Ω について行う．この Ω を信号源空間 (source space) と呼ぶ．信号成分を表すベクトル y S (t) を信号ベクトルと呼び∫∫

y S (t) =L(r)η(r)s(r, t) dr =Ωl(r)s(r, t) dr(4)Ωである．また，M × 1 の列ベクトル ε は信号成分に加法的に加わるノイズを表す．ここで，ε は平均ゼロ，分散 %2 I の正規分布：

p(ε) = N (ε|0, %2 I)(5)を仮定する．ここで I は単位行列を表す．

1.2

低ランク信号モデリングと信号およびノイズ部分空間信号源分布 s(r, t) が Q 個の離散的な信号源で構成されていて，それぞれの位置を r 1 , . . . , r Q ，強度をs1 (t), . . . , sQ (t) で表すとする．すると信号源分布は，
s(r, t) =Q∑

sq (t)δ(r − r q )(6)

q=1と表されるので，式 (4) に代入すれば信号ベクトル y S (t) は∫y S (t) =L(r)ΩQ∑

sq (t)η q δ(r − r q ) dr =q=1Q∑
q=12sq (t)lq(7)と表される．η 1 , . . . , η Q は離散信号源の向きを表す．またここで，lq は q 番目の信号源に対するリードフィールドベクトル，lq = L(r q )η q である．ここで，信号源の数 Q はセンサー数 M よりも小さい，すなわち Q < M と仮定する．この仮定が成り立つ信号を低ランク信号．計測データを低ランク信号でモデル化することを低ランク信号モデリングと呼ぶ．
式 (7) を，時刻表記 (t) を省略して表すと，

yS =Q
∑

sq lq = s1 l1 + · · · sQ lQ(8)q=1であり，この式は信号ベクトル y S がリードフィールドベクトル l1 , · · · , lQ の線形和で表されること，すなわち，y S は l1 , · · · , lQ の張る空間の要素であることを示している．ここで l1 , · · · , lQ の張る空間を信号部分空間と呼び ES と定義する．すなわち，

ES = span{l1 , . . . , lQ }(9)であり，span{· · · } は括弧の中に記載されたベクトルの張る空間を意味する．式 (8) は明らかにy S ∈ ES(10)であることを示している．
リードフィールドベクトルは線形独立であるので, 式 (9) より信号部分空間 ES の次元は Q である．
今，l1 , . . . , lQ を列に持つ行列を H ：H = [l1 , . . . , lQ ]とすれば，低ランク信号 M > Q の仮定の下で，以下の関係：

xT H = 0(11)を満たす M 次元列ベクトル x が必ず存在する．このようなベクトルの集合：

EN = {x|xT H = 0}(12)をノイズ部分空間と呼ぶ．ノイズ部分空間 EN の次元は M − Q であり，ES と EN は直行補空間をなし，ES ∪ EN = <Mおよび

ES ∩ EN = ∅となることを示すことができる．ここで，<M は M 次元実ベクトル全体の集合である．

3

1.3

信号部分空間の推定

信号部分空間は式 (9) で定義されるが，この式を直接用いて信号部分空間を求めようとすれば信号源リードフィールドベクトルが既知，すなわち信号源の位置と向きが既知でなければ求めることができない．
しかし，信号源についてのそのような情報は実際の応用では入手できないのが普通である．そのような場合でも，信号部分空間はセンサーアレイの（時空間）データ行列から推定できる．以下，信号部分空間の推定について述べる．
センサーアレイの時系列出力 y(t1 ), . . . , y(tK ) が求まっているとして，データ行列 B を，

B = [y(t1 ), . . . , y(tK )] = [y 1 , . . . , y K ](13)と定義する．ここでは簡便さのため y(tj ) を（すなわち B の j 番目の列を）y j と表記した．また，Kは計測時間点の数であり，K > M とする．さらに信号データ行列 B S を

B S = [y S (t1 ), . . . , y S (tK )] = [y S1 , . . . , y SK ](14)と定義し，B S の j 番目の列を y Sj と表す．すると式 (3) のデータモデルはデータ行列を用いて，

B = BS + Bε(15)と表せる．ここで B ε は t1 , . . . , tK でのノイズベクトルを各列にもつノイズ行列である．先に，信号デー
タ行列 B S の列ベクトルは信号部分空間の要素であること (y Sj ∈ ES ) を示したが，さらに，B S の列空間は信号部分空間に等しいこと，つまり

ES = span(B S ) = span([y S1 , . . . , y SK ])(16)の関係もある．
S式 (16) によれば B S の列空間：span([y S1 , . . . , y K ]) を推定することにより信号部分空間を推定できる．この推定は以下のように行う．まず，信号データ行列 B S の特異値分解を（M < K を考慮して）
γ1 0 · · ·0 vT 1 
 0 γ2 · · ·0  ... B S = [u1 , . . . , uM ]  .
(17).. . .
..  .
.
.
.  .
 vTM
0 0 . . . γMと表す．ここで Λ = diag(γ1 , . . . , γM ) は特異値 γj (j = 1, . . . , M ) を対角成分にもつ対角行列で，特異値は大きさの順に番号付けされているとする．u1 , . . . , uM が空間方向特異値ベクトル，v 1 , . . . , v M が時4間方向特異値ベクトルである．B S はランク Q の行列であるので，Q 個のゼロでない特異値を持つ．つまり，B S の特異値分解は









B S = [u1 , . . . , uQ , uQ+1 , . . . , uM ] 


0·······..

····

γQ···..
.

··0····..

0·····0

γ10..
.

.

.

0v T1..
.

0 ..  
Tv

.
Q T
v·   Q+1 

..

.
0 


TvM0(18)となる．この特異値分解においてノンゼロ特異値 γ1 , . . . , γQ に対応した Q 個の特異値ベクトルのスパン，span([u1 , . . . , uQ ]) が B S の列空間，すなわち，信号部分空間 ES に等しい．つまり，特異値ベクト
ル u1 , . . . , uQ は信号部分空間の正規直交基底である．
ここで，信号データ行列 B S は未知量であるのでその特異値ベクトルを直接求めることはできない．
しかし，ノイズ ε が式 (5) に従いその共分散行列が %2 I であれば，この特異値ベクトル u1 , . . . , uQ はデータ行列 B の Q 個の最大特異値に対応した特異値ベクトルに漸近的に（データ量無限大の極限で）等しくなることを以下のように示すことができる．今，B ，B S ，B ε ，のサンプル共分散行列を

1BB T
K
1
RS = B S B TSK
1
Rε = B ε B TεK
R=と定義する．式 (15) より，信号とノイズタイムコースが無相関ならばR = RS + Rεの関係がある．データ量無限大の極限でRε = %2 I

5(19)(20)(21)であるので，この仮定の下で

R = RS + %2 I
= [u1 , . . . , uQ , uQ+1 , . . . , uM ] 


0·······..

····γQ + %2···..
.

··%2····..

0·····0γ1 + %20..
.

.

.

0v T1..
.





0 ..  
T.   vQ  (22)
 T
 v·   Q+1 
.


..
0 

T2v%Mを得る．上式は R の最大 Q 個の固有値に対応した固有ベクトルは u1 , . . . , uQ に等しいこと，すなわち，B S のノンゼロレベル固有値に対応した特異値ベクトル u1 , . . . , uQ はデータ行列 B の Q 個の最大特異値に対応した特異値ベクトルに等しいことを示している．したがって，（データ量無限大の場合には）信号部分空間の推定解 EbS はEbS = span([u1 , . . . , uQ ])(23)

として求めることができる．ここで，u1 , . . . , uQ はデータ行列 B の Q 個の最大特異値に対応した特異値ベクトルである．データ量が有限の場合には，式 (23) の EbS は信号部分空間の最尤推定解であることを証明できる．フォーマルな証明は第 1.6 節に記載されている．

1.4

Signal space projection (SSP) algorithm

1.4.1

センサーノイズ除去への応用

Signal space projection（SSP）アルゴリズムは signal subspace の性質を積極的に使ってセンサーノイズや妨害信号の除去（低減）を目指すアルゴリズムである．まず，センサーノイズの除去について述べる．データモデル

y(t) = y S (t) + ε(24)のもとで計測データ y(t) の SN 比の改善を目指す．ここでセンサーノイズはやはり式 (5) にある正規ノイズを仮定する．denoising の方法は直裁である．信号部分空間へのプロジェクター P が得られればこれをデータ行列 B に乗ずることによりノイズを低減できる．プロジェクター P としては信号部分空間基底ベクトルの最尤推定解である B の特異値ベクトルから求めるのが最も簡単である．実際，式 (53)6

よりプロジェクターを空間領域の信号部分空間の最尤推定解を用いてP = [u1 , . . . , uQ ][u1 , . . . , uQ ]T(25)とすれば，これをデータ行列 B に乗じて


QM
∑
∑
b S = P B = [u1 , . . . , uQ ][u1 , . . . , uQ ]T Bγj uj v Tj  =γj uj v Tj
j=1(26)j=1を得る．上式はデータ行列の小さな特異値をゼロとおい求まる結果に等しい．データ行列を特異値分解し，そのノイズレベル特異値をゼロとおいてノイズの低減を図ろうとする方法はよく知られており，SVDフィルターと呼ばれている．つまり，信号部分空間の最尤推定結果をもちいてプロジェクターを構成し，計測データを信号部分空間に投影することによりノイズ低減を計る方法は SVD フィルターとして知られている方法に等しい．

1.5

妨害信号 (interference) 除去への応用

次に，計測信号に重畳する妨害信号 (interference) を計測データから取り除こうとするアルゴリズムについて述べる．計測データ y(t) に対してデータモデル

y(t) = y S (t) + y I (t) + ε(27)を仮定する．ここで，y I (t) は計測データに含まれる妨害信号の成分である．妨害信号源が P 個の離散的な信号源で表されると仮定する．すると，第 p 番目の妨害信号源の時刻 t での強度を σp (t)，リードフィールドを ξ p として，妨害信号 y I (t) はy I (t) =P∑

σp (t)ξp(28)p=1と表すことができる．式 (9) と同様に妨害信号部分空間（interference subspace）EI は

EI = span{ξ1 , . . . , ξ P }(29)と定義される．
今，妨害信号源の位置や向きが既知であるとして，何らかの方法でリードフィールド ξ p が推定可能であるとすれば，データを妨害信号部分空間へ投影するプロジェクター P I は，行列 H I：H I = [ξ 1 , . . . , ξ P ]と定義して，()−1P I = H I H TI H I
H TI
7(30)

で与えられる．したがって，データを妨害信号部分空間の直交補空間へ投影すれば(I − P I )y I (t) = y I (t) − P I y I (t) = y I (t) − P IP∑
p=1σp (t)ξ p
σ(t)
 1

P()−1∑..
TT
=y(t)−σp (t)ξ p = 0= y I (t) − H I H I H IHI HI .
I
p=1σP (t)(31)b S (t) をであるので，信号ベクトルの推定結果 yb S (t) = (I − P I )y(t) = (I − P I )y S (t) + (I − P I )y I (t) + (I − P I )ε = y S (t) − P I y S (t) + ε0 (32)yとして求めることができる．ここで，ε0 = P I ε とした2 ．式 (32) から，この投影により妨害信号 y I (t)が除去されることがわかる．信号成分 y S (t) の受ける影響は式 (32) の右辺第 2 項で評価できる．この項は信号源のリードフィールド lq と妨害信号源のリードフィールド ξ p の直交性が高ければ，すなわち，内積 lTq ξ p が小さければ小さくなり，プロジェクターの信号成分への影響は小さく抑えられることはあ
きらかである．
ここで，妨害磁場源のリードフィールドに関して何の情報も無ければ（すなわち，妨害信号源の位置や向きに関し情報が無ければ）第 1.3 で述べた妨害信号部分空間の最尤推定解を用いる．このとき妨害信号のみを含み（関心）信号を含まないデータ

ȳ(t) = y I (t) + ε(33)の存在を仮定する．何らかの方法でこの ȳ(t) が得られれば，ȳ(t) から作られる時空間データ行列をB̄ = [ȳ(t1 ), . . . , ȳ(tK )](34)として，この行列を特異値展開（式 (17)）して得られる空間方向の特異値ベクトル u1 , . . . , uP が妨害信号部分空間基底ベクトルの最尤推定解である．したがって，妨害信号部分空間へのプロジェクター P I は

P I = [u1 , . . . , uP ][u1 , . . . , uP ]Tとして，データを妨害信号部分空間の直交補空間へ投影すれば，妨害信号 y I (t) は y I (t) =と展開できるため，(I − P I )y I (t) = 0
2 ノイズ "0 は p("0 ) = N ("0 |0, (I − PI)%2 ) で表される正規ノイズである．

8(35)
∑P

j=1 cj ujであるので，b S (t) = (I − P I )y(t) = y S (t) − P I y S (t) + ε0y(36)の成立は簡単に示すことができる．ここでも，信号成分 y S (t) の受ける影響は基底ベクトル ui (i =1, . . . , P ) と信号源のリードフィールドベクトル lq （q = 1, . . . , Q）との直交性に依存する．

1.6

信号部分空間の最尤推定

ここでは式 (23) を証明する．式 (23) を改めて書くとEbS = span{u1 , . . . , uQ }(37)

である．ここで，改めて，データモデルはy(t) = y S (t) + εであり，信号成分 y S (t) はy S (t) =Q∑(38)

sq (t)lq(39)p(ε) = N (ε|0, %2 I)(40)

q=1と表される．ここで，ノイズ ε はと仮定する．
第 1.1 節で述べた信号およびノイズ部分空間の議論より，信号ベクトル y S (t) に直交する (M − Q)個の線形独立なベクトルが存在する．それらを xj （j = 1, . . . , M − Q）とする．すなわち，

xTj y S (t) = 0,j = 1, . . . , M − Q(41)とする．式 (13) で定義された時空間データ行列 B に対するデータ尤度は，ノイズの仮定−式 (40) のもとで，log p(B) = −K∑[y(tk ) − y S (tk )]T [y(tk ) − y S (tk )](42)

k=1となる．ただし定数倍は無視した．したがって，もし xj （j = 1, . . . , M − Q）が既知の場合，信号ベク
トル y S (t) は式 (41) の制約条件の基で，式 (42) の対数尤度を最大とする y S (t) として求まる．この制約付き最適化は以下のラグランジアン：

L[y S (tk ), κ(tk )] = −K∑[y(tk ) − y S (tk )] [y(tk ) − y S (tk )] + 2
Tk=1K
∑
k=19

Aκ(tk )y S (tk )(43)

を無制約で最大化する問題に等しい．ただしここで，A = [xj , . . . , xM −Q ](44)

κ(tk ) = [κ1 (tk ), . . . , κM −Q (tk )]T

(45)

である．
式 (43) の y S (tk ) に関する微分係数は∂L[y S (tk ), κ(tk )]= −2[y(tk ) − y S (tk )] + 2Aκ(tk )∂y S (tk )であるので，右辺をゼロとおき

y S (tk ) = y(tk ) − Aκ(tk )(46)を得る．式 (41) は以下の式：

AT y S (tk ) = AT [y(tk ) − Aκ(tk )] = 0と等価であるので，

κ(tk ) = (AT A)−1 AT y(tk )(47)を得る．上式を式 (46) に代入しなおすと[]y S (tk ) = I − A(AT A)−1 AT y(tk ) = [I − Π A ] y(tk )(48)を得る．ただし，Π A = A(AT A)−1 AT である．ここで，Π A は行列 A の列ベクトルのスパン，すなわちノイズ部分空間へのプロジェクターである．
式 (47) は，行列 A が既知の場合の y S (tk ) の最尤推定解である．y S (tk ) が式 (47) に示す値をとる時の対数尤度関数 log p(B) の値は，log p(B) = −K∑y T (tk )Π A y(tk ) = −K tr [Π A R]

(49)

k=1

を得る．ただし，R は標本協共分散行列：

1 ∑R=y(tk )y T (tk )KKk=1である．式 (48) に示す y S (tk ) は A が既知のときの最尤推定解である．A が未知のときは式 (49) に示す対数尤度関数をさらに A で最大化する．式 (49) の対数尤度関数を最大とする A が A の最尤推定解10である．式 (49) の対数尤度関数は以下の上限を持つことが知られている．

log p(B) = −KK∑

tr [Π A R] ≤ −KM∑

λj(50)

j=Q+1

k=1

ここで，λj は標本協共分散行列 R の j 番目固有値である3 ．また，この上限はプロジェクター Π A が

Π A = [uQ+1 , . . . , uM ][uQ+1 , . . . , uM ]T(51)であるときに達成される．ここで，uQ+1 , . . . , uM は固有値 λQ+1 , . . . , λM に対応した R の固有ベクトルである．したがって，ノイズ部分空間の基底ベクトルの最尤推定解は uQ+1 , . . . , uM であり，

EbN = span{uQ+1 , . . . , uM }(52)である．ノイズ部分空間と信号部分空間は直交補空間をなすため，EbS = span{u1 , . . . , uQ }が成立する．

3 固有値は大きさの順に番号付けされているとする．

11(53)