2018-5-10改訂版

空間フィルター（信号源の線形推定）
Non-Adaptive Spatial Filters

株式会社シグナルアナリシス

関原謙介

生体磁気再構成問題は劣決定な推定問題として定式化される
データベクトル y1(t ) y2 (t ) y(t ) =  yM (t )

y1 (t )y2 (t )yM (t )

脳活動によって生じた電流

M：センサー数いかにして，データベクトル y(t) と，脳活動（厳密には脳活動によって生じた電流と）の関連を表現するか？
ソースと呼ぶソース電流と計測データの関係を導く
センサーリードフィールドの導入

位置 r に 単位強度の電流源があり，x 方向を向いている場合のセンサー出力：

xxxlMxl2 l3ylMyl y l3y 2l1

l1xy 方向を向いている場合のセンサー出力：

zyyyxz 2l1zz

l1x ,l2x , ,lMxzlMzl z l3x

リードフィールド行列

l1y ,l2y , ,lMyz 方向を向いている場合のセンサー出力：

l1z ,l2z , ,lMz

以下の（M×３）の行列はリードフィールド行列と呼ばれ，センサーアレイの位置 r における感度を表す．

 l x (r ) l y (r ) l z (r ) 
11
1





L(r ) =  l jx (r ) l jy (r ) l jz (r )   
l x (r ) l y (r ) l z (r ) MM
M磁場の信号源は電流ベクトルであるが，まず，種々の空間フィルター法の導出・説明は信号源がスカラーであると仮定して行う．信号源がベクトル量である場合への拡張は後に議論する．
これは，脳内部の位置 r でソースの向き η(r ) が決まっていて，既知である場合に相当する．
この場合位置 r でのセンサーリードフィールドはリードフィールドベクトル l (r) ： l x (r ) l y (r ) l z (r )  η (r )
1
1
1
 x l(r ) = L(r )η(r ) =    ηy (r )
l x (r ) l y (r ) l z (r ) η (r )
MM
M z を用いて表す．

空間フィルターとは：基本的な定式化
weight vector
 y (t ) ↓1T==sˆ(r, t ) w=(r ) y(t ) [w1(r ), , wM (r )]  y (t )↑ M 
Estimate of source activityM∑ wm (r )ym (t )m =1

センサーデータに再構成する場所 r (フィルターポインティング部位) に依存した重みを与えてソース強度を推定する線形推定法．
フィルター重みは位置 r におけるソースの活動 s (r, t ) を選択的に計測するーvirtual (仮想)センサー

空間フィルターとは，任意の位置に局在した感度分布を持つ仮想センサーを信号処理により構成する技術と言う解釈も可能フィルター出力
空間フィルター出力： sˆ(r, t ) = w (r ) y(t )T

空間フィルター出力パワー：sˆ(r, t )2

データサンプル共分散行列
T↓=wT (r ) y(t )  wT (r ) y(t )
wT (r )Rw(r )空間フィルターノイズ出力パワー：

σ2out

Tw (=r )ε  w (r )ε 
w=(r ) εε w(r ) σTTTT

2w(r )

2↑

重みベクトルのノルム（大きさ）がフィルターのノイズゲインである

Non-adaptive spatial filter:w(r ) is data independent
重みは位置 r のリードフィールドにのみ依存する．
Adaptive spatial filter:w(r ) is data dependent
重みは位置 r のリードフィールドと計測データに依存する．
すなわち，計測データに応じたアダプティブな重みを構成できる．adaptive beamformerあるいは単にbeamformerと呼ばれる． beamformerについては別に説明する．

線形推定法と空間フィルター
ソースが存在する可能性のある領域（ソース空間）をボクセルで分割し，ソースの分布を離散化し，ボクセルにおけるソース活動の推定に線形推定（重み行列：W = [ w(r1 ), , w(rN )] )を用いれば

 s(r1 , t ) T==Wy(t )s(rN , t )

 wT (r1 ) = y(t )
 wT (r )
N 

 wT (r1 ) y(t ) 

 wT (r ) y(t )
Nとなる．
第 j 行を取り出せば，sˆ(rj , t ) = w (rj ) y(t ) であるので，rjの位置をポイントした空間フィルター表現になっている．
T

重み行列は全てのボクセルにおける空間フィルター： sˆ(r, t ) = w (r ) y(t )Tを並べたものである．
つまり，空間フィルターとは，「特別なある方法」を指すのではなく，全ての線形信号源推定法の総称である．

ミニマムノルム解
ミニマムノルム法は解 x の観測データ y に対する一致度：

F=y − Hx22を満たしながら，解のノルム： x を最小にすると言う基準を用いて

ˆ(t ) = H
ミニマムノルム解 x

T( HH T )−1 y(t ) を導いた．

ミニマムノルム解は，重み行列： W T = H T ( HH T )−1 として，xˆ(t ) = W T y(t ) と書ける．
ミニマムノルム解も（線形推定法であるので）空間フィルターとして表現できる．

Minimum-norm filter
全ボクセル位置でのlead field matrixrj

H =  l(r1 ), , l(rN )
データモデル： s (r )  1 y =  l(r1 ), , l(rN )    = Hx s(rN )
H

x

Gram matrix
G = ∫ l(r )lT (r )d r ≈ HH T

Minimum-norm solution
センサー感度の2次モーメント行列xˆ

TH=( HH T )−1 y H T G −1 yMinimum-norm filterrjMinimum-norm solution sˆ(r )   lT (r1 ) 
 1        sˆ(r )  =  lT (r )  G −1 yj j       sˆ(r )   lT (r ) N 
 N  Minimum-norm filtersˆ(r ) Solution at rj :T−1
ˆ(s rj ) = l (rj )G y


任意の位置 r において：

Tw=(r ) y; w(r ) G −1l(r )ミニマムノルム法の問題点-5

シングルソースの再構成-10

ミニマムノルム法は3次元のソース分布を再構成できない．特に，センサーから遠ざかる方向（ソースの深さ方向）に対して，精度が著しく劣る．
80年代から良く知られた事実

-15-5

再構成結果

0y (cm)5

ミニマムノルム解より，もっと性能のいい線形推定法（空間フィルター）は存在するのか??
あるとすると，どのようなものか？
空間フィルターの性能を評価する基準で一番大事なのは？

位置バイアス
再構成結果
ここに点状のソースを仮定したのにも係らず，再構成結果はかけ離れたものとなっている．

位置バイアス真の信号源位置と点像再構成結果の隔たりを位置バイアスとよぶ．
点状信号源の再構成結果（point spread function：PSF)が必ず信号源位置で最大となるピークを持つ．
この場合，再構成は位置バイアスを持たないという．

Point-spread function (reconstruction profile of a point source)(r , r0 )=wT (r )l (r0 )r0 に点状ソースが存在する場合

Location bias
点像分布関数

位置バイアスを持たない



点像分布関数は r0 でピークを持つ
(r0 , r0 ) > (r , r0 )r0

補足説明ー位置バイアス
位置バイアスは不偏性とは異なる．
不偏推定解：
真の信号源分布が再構成できる解推定解が位置バイアスを持たない．
点状の信号源に対して，真の信号源位置をピークにもつ再構成結果が得られる．

ミニマムノルム解から導かれる空間フィルター（１）Weight（noise） normalized minimum norm filterｄSPM（dynamic statistical parametric mapping）の名前で知られている．

w(r )−1−1−1G=l (r ) G l (r )G −1l(r )lT (r )G −2 l(r )重みの大きさが1になるように規格化．

ミニマムノルムフィルターの出力SN比に相当するフィルター出力を与える．
この方法は理論的には位置バイアスを持つ＊．実際の応用で位置バイアスは小さいためか， わりと人気のある手法．
A. M. Dale et al., Neuron,Vol.26, pp.55-67, 2000
＊理論的な証明は以下の文献にあり：
R. E. Greenblatt, et al., IEEE Transactions on Signal Processing, vol. 53, no. 9, pp. 3403-3412, 2005.

ミニマムノルム解から導かれる空間フィルター（２）sLORETA Standardized low-resolution tomographyw(r ) =G −1l(r )lT (r )G −1l(r )この規格化を「標準化」（standardization）と呼んでいる．この標準化に関して明確な解釈は知られていない．
sLORETAフィルターはノイズゼロの極限で位置バイアスを持たないことを理論的に示すことができる．

R. D. Pascual-Marqui, Methods and Findings in Experimental and ClinicalPharmacology,Vol.24, pp.5-12, 2002

Minimum-norm法を基にした空間フィルター（まとめ）Minimum-normdSPM

sLORETAw(r ) = G −1l(r )w(r ) =w(r ) =G −1l(r )lT (r )G −2 l(r )G −1l(r )lT (r )G −1l(r )これらは， w(r ) = ζ (r )G −1l(r ) の形で表すことができ，minimum-norm法をベースにした空間フィルターと考えることができる．

Minimum-norm法を基にした空間フィルター： w(r ) = ζ (r )G−1l (r )において，バイアスを持たないように係数 ζ (r ) を決める．

r と r’ の2点のフィルター出力について：
r に単位強度のソースがある場合

wT (r )l(r ) ≥ wT (r ')l(r )ζ (r )lT (r )G −1l(r ) ≥ ζ (r ')lT (r ')G −1l(r )r’ に単位強度のソースがある場合

r’rr’

rwT (r ')l(r ') ≥ wT (r )l(r ')ζ (r ')lT (r ')G −1l(r ') ≥ ζ (r )lT (r )G −1l(r ')ζ (r )lT (r )G −1l(r ) ≥ ζ (r ')lT (r ')G −1l(r )ζ (r ')lT (r ')G −1l(r ') ≥ ζ (r )lT (r )G −1l(r ')上の2式よりlT (r ')G −1l(r ') ζ (r ) lT (r ')G −1l(r )≥≥ TT−1l (r )G l(r ') ζ (r ') l (r )G −1l(r )Geometric meanを計算：

ζ (r )
ζ (r ')

lT (r ')G −1l(r ') lT (r ')G −1l(r )=lT (r )G −1l(r ') lT (r )G −1l(r )したがって， ζ (r ) = 1 /lT (r ')G −1l(r ')lT (r )G −1l(r )lT (r )G −1l(r ) とすれば，以下のsLORETAの重みを得る：
w(r ) =

G −1l(r )lT (r )G −1l(r )Spatial matched filter（最も簡単な空間フィルター）w(r ) = l(r ) / l(r )重みを単に再構成位置での規格化されたリードフィールドに等しいとする．

単位強度の信号源が1個，位置 r0 に存在する場合，
フィルター出力：=sˆ(t )wT y(t ) ∝lT (r )l(r0 )l(r ) l(r0 )y(t ) = l(r0 ) と表される．
l(r0 ) と l(r ) の内積

内積は2個のベクトルの類似度（spatial correlation）を表すので，再構成結果は信号源位置でピークを形成する．したがって位置バイアスを持たない．
Spatial matched filterはその「簡単さ」のため，いろいろな分野で独立に「発明されている」
レーダー・ソーナー信号処理： Delay-and-sum beamformer
CTスキャナー再構成： 逆投影

Bti-148センサーア
レイを仮定

点像分布関数
Spatial matched filter

Weight-normalized
minimum norm

-5-10-15-5

Minimum normsLORERA0y (cm)5グラム行列について考察
G = ∫ l(r )lT (r )d r :センサー感度分布の2次モーメント行列(i,j) 要素： [G ]i , j = ∫ li (r )l j (r )d r は i と j 番目のセンサーの感度分布の重なりを表す．
脳磁計センサーアレイ

CTスキャン（ランクMの行列）

1

0G≈


001

0




0

0


1


1

1G≈


どちらに近いかは1

Gのcondition number
により予測

11

1




1

1

1


（ランク1の行列）

Bti-148センサーアレイによるグラム行列の固有値スペクトル20個目で10-340個目で10-4程度まで小さくなる

コンディションナンバー
非常に大きい1
1G≈
グラム行列はどちらかと言うとに近い．
1G −1 の計算にregularizationが必要，すなわち，(G + γ I )−11  1
1  1  
1  1とする必要あり．

dSPM，sLORETAとも実際の使用においては正則化を用いる．

dSPM

w(r ) =(G + γ I )−1 l(r )lT (r )(G + γ I )−2 l(r )→γ →∞ l (r )l (r )sLORETAw(r ) =(G + γ I )−1 l(r )lT (r )(G + γ I )−1 l(r )→γ →∞ l (r )l (r )dSPM，sLORETAとも正則化定数を大きくしていくとSpatial matched filterに近づく

位置バイアス（まとめ）• Spatial matched filterとsLORETAの結果は位置バイアスを持たない．
• dSPM（Weight (noise)-normalized minimum-normフィルター）は，点状ソースに対して局在したピークを持つが，ピークは真のソース位置とは異なる位置で最大値となる，つまり位置バイアスを持つ．
• ミニマムノルム法は局在した解を持たず，ボクセル領域の上方に張り付いたような解となる．

ベクトル型空間フィルター
ソースが3次元電流ベクトルであることを考慮した拡張s (r, t )x
=sy (r, t )s (r, t ) z
 wTx (r )
T T Tw(r)b()w(r),w(r),w(r)y()
W(r ) y(t )=t=t y yz
 x
 wT (r )
 z 

ベクトル型空間フィルターではソース電流の x, y, z 成分をそれぞれ別のウエイトで再構成する．すなわち，重みは M ×3の行列となる．

ベクトル型Minimum-norm filter
全ボクセル位置でのlead field matrix
H =  L(r1 ), , L(rN ) 
Measurement equation s (r )  1 y =  L(r1 ),  , L(rN )     = Hx s (r ) N Hx

Minimum-norm solutionx

T
H=( HH T )−1 y H T G −1 yrj

ベクトル型 Minimum-norm filter ˆs (r )   LT (r1 ) 
 1       
ˆs (r )  =  LT (r )  G −1 yj j       
ˆs (r )   LT (r ) N 
 N  rj Solution at rj :ˆs (rj ) = LT (rj )G −1 y
任意の位置 r において：ベクトル型 Minimum-norm filter

Tˆs(r ) W=(r ) y; W (r ) G −1 L(r )=（参考）スカラー型： sˆ(r )T
w=(r ) y; w(r ) G −1l(r )ベクトル型weight-normalized filter (dSPM)−1

G −1 L(r )F

tr LT (r )G −2 L(r )−1

G −1l(r )−1
G=L(r ) G −1 L(r )W (r )（参考）スカラー型−1
G=l(r ) G −1l(r )w(r ){lT (r )G −2 l(r )重みが行列の場合のノルム{} tr {G L(r) G L(r)}= tr { L(r ) G L(r )}=W (r )F}

tr W (r )W (r )=TT

T−1−2A
補足 =F

2=∑ Ai , ji, j{}

tr AT A−1

ベクトル型 sLORETA

W (r ) = G L(r )  L (r )G L(r )
−1（参考）スカラー型

w(r ) =T

G −1l(r )lT (r )G −1l(r )ベクトル型 Spatial matched filter

W (r ) = L(r ) / L(r )（参考）スカラー型

w(r ) = l(r ) / l(r )

−1

−1/2

参考書
本書の第2章および第3章の内容を解説した．

参考文献