アダプティブ空間フィルター
Part II
Adaptive beamformers for MEG
source-space analysis： Part II
2017-11-0２ 改訂版

関原謙介
株式会社 シグナルアナリシス

参考書
本スライドは以下の参考書 ”Adaptive Spatial Filters for …”の第４章から第９章までについての解説である．

www.electromagneticbrainimaging.infoに正誤表などのこの書籍に対する情報がアップロードされている．

MEGデータの適用に際し生じる問題
Minimum-variance filter
フィルター重み：R −1l(r )w(r ) = Tl (r )R −1l(r )1sˆ(r, t )w=(r )Rw(r ) T
出力パワー：=l (r )R −1l(r )2TSarvas公式を使ってリードフィールドを計算する場合，リードフィー
ルドは球中心でゼロとなる．したがって，フィルター重みや出力パワーは球中心で発散してしまう．

Center-of-the-head artifactセンサーアレイ感度の不均一さ，特に，球中心で感度ゼロとなることから生じるアーチファクトへの対応．

Sarvas公式を用いたコンピュータシミュレーション
BTI148チャンネルセンサーアレイ

Minimum-variance filterによる再構成結果

コンピュータシミュレーション：ソース配置

仮定した球中心の位置にアーチファクトを生じる．

解決策１：Array-gain constraint minimum variance filterw(r )TT=argmin w Rw subject to w l(r )l (r )w

空間フィルターのゲインをリードフィールドのノルムに等しくする． l(r ) は位置 r におけるセンサーアレイのゲインであるので，この制約条件により，センサーのゲインに等しい空間フィルターを求めることになる．

ユニットゲイン制約は，反対に，センサーアレイの感度にかかわらず，フィルターのゲインを１とする方法である．

R −1 l(r )フィルター重み： w(r ) =[ lT (r )R −1l(r )]l(r ) = l(r ) / l(r )T
1(r )l(r )l2
出力パワー： =sˆ(r, t )=l T (r )R −1l(r ) lT (r )R −1l(r )ベクトルソースに対するArray-gain minimum variance filter:
Scalar formulation（スカラー型拡張）
ηT LT (r )L(r )η
ηopt = arg max T T
であるので−1ηη L (r )R L(r )η
T−1

行列 [ L (r )L(r )][ LT (r )R −1 L(r )] の最小固有値を λmin

対応する固有ベクトルをフィルター重み： w(r ) =

出力パワー：

umin とすれば：
R −1 [L(r )umin / L(r )umin ]λmin
2sˆ(r, t )=1λminベクトルソースに対するArray-gain minimum variance filter:
Vector formulation（ベクトル型拡張）
W (r ) = [ wx (r ), wy (r ), wz (r )] として以下の最適化から重み行列を求める
W (r )

TT=argmin W (r )RW (r ) subject to W (r )L(r )

L(r )

W

重み行列： W (r ) = R −1 L (r )[ LT (r )R −1 L (r )]−1 LT (r ) を得る．
出力パワー：
sx (r, t )2 + sy (r, t )2 + sz (r, t )2 =tr [W T RW ]= tr [ LT (r )R −1 L (r )]−1  (r ) = L(r ) / L(r )ただし： L解決策2：Weight normalized minimum-variance filter
Minimum-variance filterの重みをnormalizeして用いる．
Tw (r )
R −1l(r )ˆ(srt=y(t ) として出力を求める,)w(r ) = Tを用いて−1w(r )l (r ) R l (r )w(r )=w(r )

w(r ) =w(r )=wT (r )w(r )R −1l(r )T−2l (r ) R l (r )出力パワー：

R −1l(r )T−2l (r ) R l (r )であるので：として計算したフィルター重みを用いることに等しいsˆ(r, t )2lT (r )R −1l(r )= Tl (r )R −2 l(r )と表される．

Weight normalized minimum-variance filter-続き
Weight normalized minimum-variance filterは以下の最適化から導くことができる．

w(r )

TT=w
Rw
subjecttow
w 1argminw

そのため，unit noise-gain minimum-variance filterとも呼ばれる．
（世の中では，noise-normalized beamformerと呼ばれる場合もある．）
この空間フィルターはminimum-varianceフィルターと比べて，高い空間分解能を与える．ただし，再構成結果の定量性に不安定さがあり，積極的な使用は進めない．
Signal Processing分野ではBorgiotti-Kaplan beamformerと呼ばれ，古くから知られている．
Borgiotti, G., and L. Kaplan. "Superresolution of uncorrelated interference sources by using
adaptive array techniques." IEEE Transactions on Antennas and Propagation 27.6 (1979):
842-845.

ベクトルソースに対するWeight normalized minimum-variance filter :Scalar formulation（スカラー型拡張）
ηT LT (r )R −1 L(r )η
ηopt = arg max T Tであるので−2ηη L (r )R L(r )η
T

行列 [ L (r )R−1L(r )]−1 [ LT (r )R −2 L(r )] の最小固有値を λmin
umin とすれば： ηopt = u min である．

対応する固有ベクトルをフィルター重み：

出力パワー：

R −1 L(r )uminw(r ) =uTmin LT (r )R −2 L(r )umin2sˆ(r, t )=1λminベクトルソースに対するWeight normalized minimum-variance filter :Vector formulation（ベクトル型拡張）以下の最適化を用いる

wxT
wTx ly (r ) 0,=
wTx lz (r ) 0
to wTx wx 1,=
argmin wx Rwx subject =wx=
wy argmin wTy Rwy subject to =wTy lx (r ) 0,=
wTy wy 1, =wTy lz (r ) 0wy=
wz argmin wTz Rwz subject to =wTz lx (r ) 0, =wTz ly (r ) 0,=
wTz wz 1wz

フィルター重み行列は以下で与えられる．
1 / ϒ1,1
[ wx (r ), wy (r ), wz (r )] = R −1 L(r )[ LT (r )R −1 L(r )]−1  0
 00
1 / ϒ 2,20

ここで，  = [ LT (r )R −1 L(r )]−1[ LT (r )R −2 L(r )][ LT (r )R −1 L(r )]−1

0
1 / ϒ 3,3 
0Sarvas公式を用いたコンピュータシミュレーション
BTI148チャンネルセンサーアレイ

コンピュータシミュレーション：ソース配置

ユニットゲインビームフォーマー

アレイゲインビームフォーマーソース配置スカラー
アレイゲイン

ソースタイムコース

スカラー
ノイズゲインベクトル
アレイゲイン
センサータイムコース

ベクトル
ノイズゲイン補足：再構成結果の次元について
シングルソースの仮定： y(t ) = s(r, t )l(r )TT=sˆ(r, t ) w==(r ) y(t ) s(r, t )w(r )l(r ) s(r, t )unit-gain filter：

再構成結果は電流の次元を持つ．
TTarray-gain filter：=sˆ(r, t ) w==(r ) y(t ) s(r, t )w(r )l(r ) s(r, t ) l(r )再構成結果は磁場の次元を持つ．
unit noise-gain filter：sˆ(r, t ) w==(r ) y(t ) s(r, t )TlT (r )R −1l(r )T−2l (r ) R l (r )再構成結果は磁場の次元を持つ．

≈ s (r , t ) l (r )array gain, unit-noise gain filterの再構成結果の解釈について
array-gain filterの再構成結果： sˆ(r, t ) = s(r, t ) l(r ) はセンサーアレイが検出する磁場強度の「ものさし」で測ったソース強度と解釈できる．

unit noise-gain filterの再構成結果： sˆ(r, t ) = s(r, t )1
は，センサーノイズの分散をσ2として，σlT (r )R −1l(r )lT (r )R −2 l(r )sˆ(r, t ) = s(r, t )lT (r )R −1l(r )σ lT (r )R −2 l(r )がフィルタ出力のSN比と解釈できる．したがって，再構成結果 sˆ(r, t ) / σ
画像として表示すれば，これは信号源推定のSN比を画像化したものになる．このようにある種の「統計量」を画像化することをStatisticalParametric Mappingと呼ぶ．

ビームフォーマーのベクトル・スカラー型の導出の詳細については，「Adaptive
spatial filter..」の第４章に記載がある．また，以下の論文がnoise-gain
constraint vector beamformerに関するオリジナル論文である．

また，スカラー型ビームフォーマーについては以下の論文が最初の提案である．

Biomag 96: Proceedings of the Tenth International Conference on Biomagnetism,
C. J. Aine eds. Feb. 1996, Page(s): 338-341, Springer-Verlag, New York.

環境ノイズ(low-rank interference)，およびブレインノイズ(high-rank interference)のビームフォーマー出
力への影響

Performance degradation due to interferences
Measurements in non-idealistic, real-life conditions
contains various types of interferences.
Beamformerに対するinterference（干渉信号）の影響を議論する．
2種類の干渉信号について議論する•Low-rank interferences•High-rank interferencesLow-rank interference
人工的な原因で発生する干渉信号

External disturbance with artificial origins•Power-line leakage•Additive base-line drift•Artificial noise from electrical appliances⇓
modeled as low rank signals
このような妨害信号の振る舞いは少数の“固有成分”によって記述できる

Low-rank interference—further assumptions
信号磁場↓
Additivity： =y(t )干渉・妨害信号↓Q∑ l(rq )sq (t ) + ε + d (t )q =1Covariance matrix relationshipTotal covariance↓

No correlation：=R

RS +ε + RD↑
Signal-plus-sensor-noise covarianceInterference covariance

RD = d (t )d T (t )PD

RD = ∑ λj uj uTj
j =1

When RD ∝ uuT （RDをランク１の行列と仮定）
Spatial-filter outputs:lT (r )R −1 y(t ) lT (r )( RS +ε + RD )−1 y(t )=lT (r )R −1l(r ) lT (r )( RS +ε + RD )−1 l(r )sˆ(r, t )QT−1S +εl (r ) Rl(rq )= ∑ sq (t ) Tq =1l (r )RS−+1 ε l(r )[1 −cos(l(r ), u | RS−+1 ε )cos(u, l(rq ) | RS−+1 ε )−1S +εcos(l(r ), l(rq ) | R)

[1 − cos2 (l(r ), u | RS−+1 ε )]This part is equal to 1≈lT (r )RS−+1 ε y(t )](when u is very different from l(r) )lT (r )RS−+1 ε l(r )Spatial-filter outputs when the interference does not exist

assumed source waveform
1148-channel sensor array

0-1
10-5-11}three sources-10100200
latency (ms)300400z (cm)0-10-15-50y (cm)5sphere for the forward modelingSimulated magnetic recordings

No interference

Periodicinterference

Periodicinterferencewith random
phasesLinear trendwith random
inclinations

Randomlinear trendwith slow
wave activity

Visualization of the first eigenvectors of the interferences

typical lead fieldof a brain source

PeriodicinterferenceLinear trendwith random
inclinations

Periodicinterference withrandom phases

Random linear
trend with slow
wave activity

Reconstruction results

No disturbance

PeriodicinterferenceLineartrend withrandominclinations

Periodicinterferencewith random
phasesRandomlineartrend withslow wave
activity

Reconstructed
time courses

10-1

No disturbance000

Periodic latency (ms)interferenceLineartrend withrandominclinations

100200300400Periodicinterferencewith random
phasesRandom linear
trend with slow
wave activity

Beamformerを用いてDBSを装着した患者の脳磁界からの信号源推定が行われている．

参考文献
Beamformer再構成に対するLow-rank interferenceの影響については「Adaptive spatial filter..」の第7章で議論されている．また，以下がオリジナル論文である．

High-rank interference
signal source

neurophysiological noise source

P∞

j =1

ky(t ) = ∑ l(rj )s(rj , t ) + ∑ l(rk )(rk , t ) + ε(t )(rk , t ) 皮質にランダムに分布したインコヒーレントな微弱ソースとしてモデル化できる．

The rank of y(t) is high.
•de Munck et al., IEEE Trans. Biomed. Eng., 39, 791-804, 1992.
•Valdes et al., Brain Topography, 4, 309-319, 1992.
•Lutkenhoner, J. Appl. Phys., 75, 7204-7210, 1994.

Spontaneous
brain activity

誘発脳磁界計測に重畳する自発脳磁界はBrain noiseと呼ばれる

Influence of background source activity(Influence of brain noise)
Computer-generated recordings→
Brain noise

Brain-noise added recordings→

Background activity causes a severe blur.

CorreltioncoefficientDual-condition experiments
Task: y(t ) = yS (t ) + yI (t ) + ε
Control: y=(t ) yI (t ) + ε
CCovariance matrix relations

R RS + Ri +ε
Task: =Control: RC = Ri +ε
Problem:• Brain noiseはソースの存在領域から発生する．除去にはコントロールデータ（Brain noiseのみを含み，関心信号を含まないデータ）が必要である．

Prewhitening estimation of signal covariance−1/2−1/2taskとcontrolデータから R =RC RRC を計算する．

R = RS + Ri +ε ⇒ R =R S + I

この計算をノイズの白色化（prewhitening）と呼ぶ

Q

R S の固有値展開を R S = ∑ γ j uj uTj とすれば，R の固有値展開は：
j =1

Q

M

QT
R = ∑ γ u u + ∑ u u =∑ (γ j + 1)uj uj .
Tj j jj=
1

Tj jj=
Q +1

j=1となる．したがって：

R の信号レベル固有値に対する固有ベクトルが R S の固有ベクトルに等しい
信号レベルとノイズレベル固有値に対する閾値の設定が必要

Prewhitening beamformingR =RC−1/2 RRC−1/2 を計算する
R の信号レベル固有値に対する固有ベクトルから R S を推定する：
QR S = ∑ γ j' uj uTj
j =1

ˆ を推定する：R S から，prewhiteningの演算を逆に行い RSQ

T 
1/2
1/2
'1/2ˆ
=RR=RRRuuRγ∑SS CC j j j  C
 j =1

1/2
Cˆ を用いてbeamformingを行う：
RS
ˆ −1 ˆy (t )l T (r ) R
TSS
ˆ
;
()
[][]sˆ(r, t ) =ytu
uu
uy(t )
S1Q
1QT−1
ˆ l (r )l (r ) RS
ブレインノイズの影響を効率よく除去できるものの，Qの決定が難しい．

Auditory evoked field

Conventional

73ms
Prewhitening90msConventional
Prewhitening参考文献
Prewhitening beamformerについては「Adaptive spatial filter..」の第8章で議論されている．また，以下の2報の論文がオリジナル論文である．

自発脳磁界解析への応用：
刺激に対してロックしていない脳活動（Inducedactivity）の検出と解析

Induced activity
Task-related modulation of oscillatory brain activity（タスクによって引き起こされる脳律動の変調成分）
Induced activity:•Stimulus-evoked but not phase-locked to the stimulus•frequency specific↓Event-related (spectral) power change
Power decrease → Event-related desynchronization (ERD)
Power increase → Event-related synchronization (ERS)
G. Pfurtscheller and F. H. Lopes da Silva: “Event-related EEG/MEG synchronization anddesynchronization: basic principles,” Clin. Neurophysiol., 1997

Induced activity may not be detected in averaged results
Stimulus onset

Motor-evoked MEG

trial 1100-trial-averaged results

trial 2

Subject A

trial 3trial K

Subject B

average
Induced activity has time jitter
No clear premotor response exists in average waveformCovariance matrices are computed from non-averaged trials, andaveraged across epochs
Stimulus onset
trial 1

⇒ R1

trial 2⇒ R2Ptrial 3⇒ R3

R = ∑ Rjj =1

trial K⇒ RPThis R contains unwanted backgroundinterferenceDual-condition experiments
Task: y(t ) = yS (t ) + yI (t ) + ε
Control: y=(t ) yI (t ) + ε
CCovariance matrix relations

R RS + Ri +ε
Task: =Control: RC = Ri +ε
Problem:
この場合， yS (t ) およびを持つ．

yI (t ) とも自発脳磁解であり，同程度の強度

Prewhitening beamforming
Task: =R RS + Ri +εControl: RC = Ri +ε⇒ Calculate R =RC−1/2 RRC−1/2

Signal covariance estimation
Q1/2T 
1/2
1/2 ˆ
=RS RC U S U S ( R − I ) =RCRC  ∑ γ j uj u  RC1/2
 j =1↑U S = [ u1, , uQ ] : R の信号レベル固有ベクトルyS (t ) および yI (t ) とも自発脳磁解であり，同程度の強度を持つため，Q を決めることが難しい．(の固有値スペクトルは明瞭な閾値を持たない)
R

Robinson’s F ratio method
Define s(r, t ) : source image from y(t )sC (r, t ) : source image from yC (t )Rtotal = R + RCsˆ(r, t )F (r ) =−1−1lT (r )Rtotaly(t )lT (r )RtotalyC (t )
=and sˆC (r, t )−1−1Tl (r )Rtotal l(r )lT (r )Rtotall (r )s(r, t )2 − sC (r, t )2sC (r, t )2S. E. Robinson and J. Vrba, Recent Advances in Biomagnetism, pp. 302-305, 1999

Frequency-domain beamformer(Narrow-band beamforming)
Because the induced signal is frequency specific, another strategy
that further reduces the brain-noise-influence is to use the weighttuned to the frequency of the induced signal.

Frequency-specific weight (Narrow-band weight)R( f )

Γ ( f )Γ H ( f ) where y(t ) ←→Γ(f)F .T .
T

w(r, f ) = R( f ) l(r ) / [ l (r )R( f ) l(r )]R( f ) : Cross-spectral matrix calculated from the frequency
band of the induced signal.

Time-frequency-domain beamformer
Time-frequency-specific weighty(t ) ←→Γ (τ , f )Windowed F.T.
Induced activityが含まれている時間周波数領域の中心時間をτ, 中心周波数をfとすれば：R(τ , f ) = Γ (τ , f )Γ H (τ , f )Tw(r, τ , f ) = R(τ , f ) l(r ) / [ l (r )R(τ , f ) l(r )]ここで，記号 ⋅ はエポックでのアベレージを表す．連続的に計測した
自発脳磁解をエポック化し，エポックでのアベレージを計算する

Results of frequency-specific-prewhitening filter for handmotor measurement
Voluntary right-finger movement (every 10 sec)

CT⇓
RLA

RL

TF maps of data from all sensors
A localized source is found near contra-lateral M1 area

PComparison between prewhitening and F-image methods
MRI overlay

Subject A

Subject B

Prewhitening results

F image results

Spike-locked dual-state adaptive spatial filter
発作間欠期のてんかん性スパイクの信号源推定をベータおよびガンマ帯のパワーチェンジから行う．
60分のcontinuous resting state MEGからスパイクを含む複数個のデータ部位を選び出す．

ControlとTask共分散行列を以下のように計算
Control (TC)

Task (T)RC = ∑∑ R( f ,TC )R= ∑∑ R( f ,T )spikes f ∈FW

spikes f ∈FW-700-500-100100

FW :12-55Hz

A. G. Guggisberg et al., NeuroImage Vol.39, 2008

Results from patients #1
Time-frequency maps

Source image

Surgical resection area

外科的切除部位はMEGのみでなく，他の検査も併用し総合的に決められる．
したがって切除部位がてんかん性スパイクの信号源位置推定の「正解」を与えると考え，切除部位と信号源推定結果を比較することで，信号源推定の正しさを評価する．

Results from patients #2
この患者の場合にはbeamformerの結果は外科的切除部位と一致しているが，ダイポールクラスターの位置は切除部位の位置からは大きく外れている

Time-frequency maps

Source image

Surgical resection area

Dual-state spatial filter をさらに拡張することで，時間周波数領域全体に対して信号源再構成が可能となるControl window

Task windowfrequency

f j +1fj( f j , tc )( f j , tk )f j −1
tktk +1

tk +2

tk + 3tcs(r, f j , tk ) : 時間周波数領域 ( f , t ) に対する再構成結果
jk

A set of s(r, f=1,2,...
=: k 1,2,... represents
j , tk ) : j
time-frequency reconstruction of source activity.

Five dimensional imaging from hand-motor data
Beta-band activity
LRLRLRLRLRLL

R

R

High-gamma-band activity
LR

L

R

LR

Comparison with ECoG results
Right finger (RD2) movement

Time-frequency beamformingについては以下の論文に詳細な記載がある．