第5章: ICAによるアーティファクト除去¶
フィルタリング(第4章)は不要な周波数を除去します。しかし、まばたきや筋アーティファクトは、周波数の点で本物の脳信号と重なっているため、それらを壊さずに単純にフィルタで取り除くことはできません。ICA(独立成分分析)は代わりに、信号を空間的に、統計的に独立な発生源へと分離します。これにより、アーティファクトの発生源だけを特定して除去し、それ以外はすべてそのまま残すことができます。
import mne
import matplotlib.pyplot as plt
from pathlib import Path
# 日本語の文字がグラフ内で正しく表示されるよう、フォントを設定する
plt.rcParams["font.sans-serif"] = ["Yu Gothic", "Meiryo", "MS Gothic", "DejaVu Sans"]
plt.rcParams["axes.unicode_minus"] = False
mne.set_log_level("WARNING")
sample_folder = Path(mne.datasets.sample.data_path()) / "MEG" / "sample"
raw = mne.io.read_raw_fif(sample_folder / "sample_audvis_raw.fif", preload=True)
raw.pick(["eeg", "eog"])
raw.filter(l_freq=1.0, h_freq=40.0) # ICAは、あらかじめ1 Hz以上でハイパスフィルタされたデータで最もよく働く
raw
| General | ||
|---|---|---|
| Filename(s) | sample_audvis_raw.fif | |
| MNE object type | Raw | |
| Measurement date | 2002-12-03 at 19:01:10 UTC | |
| Participant | Unknown | |
| Experimenter | MEG | |
| Acquisition | ||
| Duration | 00:04:38 (HH:MM:SS) | |
| Sampling frequency | 600.61 Hz | |
| Time points | 166,800 | |
| Channels | ||
| EEG | and | |
| EOG | ||
| Head & sensor digitization | 146 points | |
| Filters | ||
| Highpass | 1.00 Hz | |
| Lowpass | 40.00 Hz | |
ICAの仕組み(概念的に)¶
頭皮の電極を、いくつかの独立した「スピーカー」が同時に鳴っている部屋の中の、複数のマイクだと考えてみてください。スピーカーの中には脳の発生源もあれば、あるものは「まばたき」、別のものは「筋緊張」かもしれません。それぞれの電極は、距離や幾何学的配置に応じて、すべてのスピーカーの異なる混合を拾います。ICAはどのスピーカーがどれかを知りませんが、できるだけ統計的に互いに独立した成分を探すことで、記録を(想定される)独立な元の発生源へと数学的に「アンミックス(逆混合)」できます。
成分に分解したら、それぞれを確認し、(その頭皮上のトポグラフィと時間経過に基づいて)どれが眼球/筋のアーティファクトらしいかを判断し、それだけを除去してから、それ以外のすべてを再び混ぜ合わせます。これはフィルタリングより外科的な手法です。周波数帯域ではなくアーティファクトの発生源そのものを除去するため、たまたまアーティファクトと同じ周波数を持つ本物の脳信号も生き残ります。
ica = mne.preprocessing.ICA(n_components=15, random_state=97, max_iter="auto")
ica.fit(raw)
ica.plot_components()
plt.show()
それぞれのサブプロットは、1つの成分の頭皮上のトポグラフィ(その成分がどのチャンネルに最も強く寄与しているか)を表します。典型的なまばたき成分は、最前部(Fp1/Fp2 付近)を中心にした、強く左右対称な塊のように見えます。まばたきはほとんどの記録において、前頭部で最も支配的な信号源だからです。
EOGチャンネルを使った自動検出¶
目視での確認も有効ですが、主観的で時間もかかります。もし記録に専用のEOGチャンネルが含まれている場合(このサンプルデータセットには含まれており、眼球運動を直接計測しています)、MNEはどのICA成分がそれと強く相関しているかを自動的に見つけられます。トポグラフィを目で確認するより、はるかに信頼性が高い方法です。
eog_indices, eog_scores = ica.find_bads_eog(raw)
print("まばたきに関連するとフラグ付けされた成分:", eog_indices)
ica.plot_scores(eog_scores)
plt.show()
まばたきに関連するとフラグ付けされた成分: [np.int64(0)]
フラグ付けされた成分をもう少し詳しく見てみましょう。そのトポグラフィ、時間経過、実際のEOG信号との相関を、plot_properties が1つの診断図にまとめてくれます。
if eog_indices:
ica.plot_properties(raw, picks=eog_indices)
plt.show()
成分を除外してクリーンアップを適用する¶
除去したい成分のインデックスを ica.exclude に設定し(自動検出の結果でも、自分で目視確認した結果でも構いません)、その後 ica.apply() を呼ぶと、それらの成分を差し引いた形で信号が再構成されます。
ica.exclude = eog_indices # または上で見た内容に基づいて、例えば [0, 3] のように手動で設定する
raw_clean = raw.copy()
ica.apply(raw_clean)
raw.plot(picks="eeg", n_channels=10, duration=10, scalings="auto", title="ICAクリーンアップ前")
raw_clean.plot(picks="eeg", n_channels=10, duration=10, scalings="auto", title="ICAクリーンアップ後")
plt.show()
EOGチャンネルがない場合¶
あなた自身の記録(32チャンネル、CPz を基準)には専用のEOGチャンネルがありません。これはよくあることです。その場合は次のようにします。
- 目視確認に頼ります。
ica.plot_components()とica.plot_sources(raw)(各成分の時間経過を表示するので、まばたきのようなスパイクを文字通り目で追うことができます)を使い、上で説明した前頭部トポグラフィの経験則を使って、まばたき/筋活動らしい成分を目で見つけます。 - あるいは、もっと単純な(そして大まかな)代替策として、エポック化した後に振幅ベースでエポックを除外する方法があります。いずれかのチャンネルのピーク・ツー・ピーク振幅が異常に大きいエポックをまるごと破棄するもので、これは第6章で2番目の防衛線として使います。
まとめ¶
- ICAは信号を独立な成分に分離し、そのうち眼球/筋アーティファクトらしいものだけを除去します。
- EOGチャンネルがあれば
ica.find_bads_eog()でこれを自動化でき、なければica.plot_components()/ica.plot_sources()を目視確認します。 - ICAは必ず、ハイパスフィルタ済み(1 Hz以上)のデータに対してフィットさせます。
次へ: 第6章 — イベントとエポック化