说明
转到末尾下载完整示例代码。
投影器与投影背景知识#
本教程介绍投影器与信号空间投影（SSP）的背景知识，涵盖投影器的加载与保存、向 Raw 对象添加和移除投影器、“已应用”与“未应用”投影器之间的区别，以及 MNE-Python 在哪些阶段会自动应用投影器。
首先导入所需的 Python 模块；同时定义一个简短的函数，以便更方便地绘制多张外观相似的图：
# Authors: The MNE-Python contributors.
# License: BSD-3-Clause
# Copyright the MNE-Python contributors.
import os
import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.mplot3d import Axes3D # noqa
from scipy.linalg import svd
import mne
def setup_3d_axes():
ax = plt.axes(projection="3d")
ax.view_init(azim=-105, elev=20)
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("z")
ax.set_xlim(-1, 5)
ax.set_ylim(-1, 5)
ax.set_zlim(0, 5)
return ax
什么是投影？#
用最通俗的话说，投影是一种把一组点转换为另一组点的运算，并且对结果点再次执行同样的投影不会再产生任何影响。举个简单的几何例子，设想三维空间中的点 \((3, 2, 5)\)。把该点投影到 \(x, y\) 平面上，效果很像太阳位于正上方时这个点投下的影子：
ax = setup_3d_axes()
# plot the vector (3, 2, 5)
origin = np.zeros((3, 1))
point = np.array([[3, 2, 5]]).T
vector = np.hstack([origin, point])
ax.plot(*vector, color="k")
ax.plot(*point, color="k", marker="o")
# project the vector onto the x,y plane and plot it
xy_projection_matrix = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 0]])
projected_point = xy_projection_matrix @ point
projected_vector = xy_projection_matrix @ vector
ax.plot(*projected_vector, color="C0")
ax.plot(*projected_point, color="C0", marker="o")
# add dashed arrow showing projection
arrow_coords = np.concatenate([point, projected_point - point]).flatten()
ax.quiver3D(
*arrow_coords,
length=0.96,
arrow_length_ratio=0.1,
color="C1",
linewidth=1,
linestyle="dashed",
)
说明
@ 符号表示 NumPy 数组上的矩阵乘法，自 Python 3.5 / NumPy 1.10 起引入。plot(*point) 这种写法利用 Python 的参数展开机制，把 point 的元素“解包”为函数的各个独立位置参数。换言之，plot(*point) 会展开为 plot(3, 2, 5)。
注意，我们使用矩阵乘法计算了点 \((3, 2, 5)\) 在 \(x, y\) 平面上的投影：
……并且对结果再次施加投影，得到的仍然只是这个结果本身：
从信息的角度看，这次投影取用点 \(x, y, z\)，并去掉了该点在 \(z\) 方向上位置远近的信息；现在我们所知道的只是它在 \(x, y\) 平面上的位置。此外，把投影矩阵作用于 \(x, y, z\) 空间中的任意点，都会将其归约为 \(x, y\) 平面上对应的点。描述这一过程的术语是子空间（subspace）：投影矩阵把原空间中的点投影到一个维数低于原空间的子空间中。我们的子空间之所以是 \(x,y\) 平面（而不是比如 \(y,z\) 平面），直接取决于投影矩阵中的具体取值。
示例：将投影用作降噪#
描述这种“信息损失”或“投影到子空间”的另一种说法是：投影降低了测量的秩（或“自由度”）——在这里是从 3 维降到 2 维。反过来说，如果你知道 \(z\) 方向上的测量成分只是测量方法带来的噪声，而你只关心 \(x\) 和 \(y\) 成分，那么把三维测量投影到 \(x, y\) 平面就可以视为一种降噪方式。
当然，如果所有测量噪声恰好都集中在 \(z\) 方向上，那确实非常幸运；这样你直接丢弃 \(z\) 成分即可，不必费心构造投影矩阵或做矩阵乘法。换一个情形：假设为了完成测量，你必须扣动测量设备上的触发器，而扣动触发器的动作会使设备轻微移动。如果你能测出扣动触发器如何影响设备位置，就可以“校正”实际测量，从而把扣动触发器的影响“投影掉”。这里假定触发器的平均效应是使测量设备移动 \((3, -1, 1)\)：
trigger_effect = np.array([[3, -1, 1]]).T
知道了这一点，我们就能计算一个与触发器效应正交的平面（这里用到的事实是：给定法向量 \((A, B, C)\)，过原点的平面方程为 \(Ax + By + Cz = 0\)），并把实际测量投影到该平面上。
由 trigger_effect 向量计算投影矩阵是通过奇异值分解（SVD）完成的；感兴趣的读者可以查阅网络资料或线性代数教材了解该方法的细节。有了投影矩阵，就可以对原始向量 \((3, 2, 5)\) 做投影以去除触发器的影响，然后再绘图：
# compute the projection matrix
U, S, V = svd(trigger_effect, full_matrices=False)
trigger_projection_matrix = np.eye(3) - U @ U.T
# project the vector onto the orthogonal plane
projected_point = trigger_projection_matrix @ point
projected_vector = trigger_projection_matrix @ vector
# plot the trigger effect and its orthogonal plane
ax = setup_3d_axes()
ax.plot_trisurf(x, y, z, color="C2", shade=False, alpha=0.25)
ax.quiver3D(
*np.concatenate([origin, trigger_effect]).flatten(),
arrow_length_ratio=0.1,
color="C2",
alpha=0.5,
)
# plot the original vector
ax.plot(*vector, color="k")
ax.plot(*point, color="k", marker="o")
offset = np.full((3, 1), 0.1)
ax.text(*(point + offset).flat, "({}, {}, {})".format(*point.flat), color="k")
# plot the projected vector
ax.plot(*projected_vector, color="C0")
ax.plot(*projected_point, color="C0", marker="o")
offset = np.full((3, 1), -0.2)
ax.text(
*(projected_point + offset).flat,
"({}, {}, {})".format(*np.round(projected_point.flat, 2)),
color="C0",
horizontalalignment="right",
)
# add dashed arrow showing projection
arrow_coords = np.concatenate([point, projected_point - point]).flatten()
ax.quiver3D(
*arrow_coords,
length=0.96,
arrow_length_ratio=0.1,
color="C1",
linewidth=1,
linestyle="dashed",
)
和之前一样，该投影矩阵会把 \(x, y, z\) 空间中的任意点映射到那个平面上，而一个点一旦被投影到该平面，再次施加投影就不会产生任何效果。因此可以清楚地看到，尽管投影后的点在 \(x\)、\(y\)、\(z\) 三个方向上都有变化，这组投影点实际上只有两个有效维度（即它们被约束在一个平面上）。
EEG 或 MEG 信号的投影运作方式非常相似：点 \(x, y, z\) 对应单个时间点上每个传感器的取值，而投影矩阵会随你想投影掉信号的哪些方面（即哪种噪声）而变化。唯一真正的区别在于，你处理的不再是单个三维点 \((x, y, z)\)，而是由 \(N\) 维“点”组成的时间序列（每个采样时刻一个点），其中 \(N\) 通常为几十到几百（取决于你的 EEG/MEG 系统有多少传感器）。幸运的是，由于投影是矩阵运算，即使信号有数百维、数万个时间点，计算也能非常快速地完成。
信号空间投影（SSP）#
上文提到，投影矩阵会随你想投影掉的噪声类型而变化。信号空间投影（SSP）[1] 通过比较包含与不包含目标信号的测量，来估计该投影矩阵应当是什么。例如，你可以额外采集“空房间”（empty room）记录，即在没有被试在场时记录传感器上的活动。通过观察空房间测量中各 MEG 传感器的空间活动模式，可以构造一个或多个 \(N\) 维向量，给出环境噪声在传感器空间中的“方向”（类似于上例中“触发器效应”的向量）。SSP 也常用于去除心跳和眼动伪迹——在这些情况下，噪声方向不是由空房间记录估计，而是通过检测伪迹、提取其周围的 epoch 并取平均来估计。示例参见 使用 SSP 修复伪迹。
一旦确定了噪声向量，就可以构造一个与它们正交的超平面，并建立投影矩阵，把实验记录投影到该超平面上。这样，测量中与环境噪声相关的成分就能被去除。同样应当清楚，投影会降低数据的维数——传感器信号的数量不变，但它们不再全部线性无关——不过传感器通常有几十到几百个，而要消除的噪声子空间只有 3–5 维，因此自由度的损失通常并不构成问题。
MNE-Python 中的投影器#
在我们的示例数据中，SSP 已经使用空房间记录完成，但投影器与原始数据一同存储，且尚未被应用（换句话说，投影器尚未激活）。下面加载示例数据并截取为 60 秒；从 read_raw_fif() 的输出中可以看到这些投影器：
sample_data_folder = mne.datasets.sample.data_path()
sample_data_raw_file = os.path.join(
sample_data_folder, "MEG", "sample", "sample_audvis_raw.fif"
)
raw = mne.io.read_raw_fif(sample_data_raw_file)
raw.crop(tmax=60).load_data()
Opening raw data file /home/circleci/mne_data/MNE-sample-data/MEG/sample/sample_audvis_raw.fif...
Read a total of 3 projection items:
PCA-v1 (1 x 102) idle
PCA-v2 (1 x 102) idle
PCA-v3 (1 x 102) idle
Range : 25800 ... 192599 = 42.956 ... 320.670 secs
Ready.
Reading 0 ... 36037 = 0.000 ... 60.000 secs...
在 MNE-Python 中，环境噪声向量使用主成分分析（principal component analysis，通常缩写为“PCA”）计算，因此 SSP 投影器通常命名为“PCA-v1”之类。（顺便一提，由于执行 PCA 的过程底层会用到奇异值分解，已发表论文中也常见“投影器使用 SVD 计算”这样的说法。）投影器存储在 raw.info 的 projs 字段中：
print(raw.info["projs"])
[<Projection | PCA-v1, active : False, n_channels : 102>, <Projection | PCA-v2, active : False, n_channels : 102>, <Projection | PCA-v3, active : False, n_channels : 102>]
raw.info['projs'] 是一个普通的 Python 列表，元素为 Projection 对象，因此可以通过索引访问单个投影器。Projection 对象本身类似于 Python 的 dict，所以可以用它的 .keys() 方法查看其中包含哪些字段（通常不需要直接访问其属性，但如有必要也可以这样做）：
first_projector = raw.info["projs"][0]
print(first_projector)
print(first_projector.keys())
<Projection | PCA-v1, active : False, n_channels : 102>
dict_keys(['desc', 'kind', 'active', 'data', 'explained_var'])
Raw、Epochs 和 Evoked 对象都有一个布尔型 proj 属性，用于指示对象中是否存有未应用/未激活的投影器。换言之，如果至少存在一个投影器且它们全部处于激活状态，则 proj 属性为 True。此外，每个单独的投影器也有一个布尔型 active 字段：
print(raw.proj)
print(first_projector["active"])
False
False
计算投影器#
在 MNE-Python 中，可以使用通用函数 mne.compute_proj_raw()、mne.compute_proj_epochs() 和 mne.compute_proj_evoked() 计算 SSP 向量。这些函数的一般假设是：传入的数据包含原始数据、epoch，或你想通过投影修复的伪迹的平均值。实践中，这通常涉及空房间记录的连续原始数据，或平均后的 ECG、EOG 伪迹。此外还提供了一组面向典型用例的高级便捷函数，用于计算投影向量，其中包括 mne.preprocessing.compute_proj_ecg() 和 mne.preprocessing.compute_proj_eog()，分别用于计算 ECG 和 EOG 相关的伪迹成分；这些用法的示例参见 使用 SSP 修复伪迹。若要把 EEG 参考信号作为投影器计算，可以使用 mne.set_eeg_reference() 函数；更多信息参见 设置 EEG 参考。
警告
最好仅在实际会使用的通道上计算投影器（例如排除坏通道）。这样可以确保投影向量保持正交归一化，并正确捕捉感兴趣的活动。
可视化投影器的效果#
通过比较应用与不应用投影器的图，可以看到投影器对测量信号的影响。默认情况下，raw.plot() 会在绘图前于后台应用投影器（但不会修改 Raw 对象）；可以用下面展示的布尔参数 proj 控制这一行为，也可以通过绘图窗口右下角的 Proj 按钮打开投影器界面，交互式地开启或关闭投影器。这里我们只查看磁强计，并取文件开头 2 秒的样本。
mags = raw.copy().crop(tmax=2).pick(picks="mag")
for proj in (False, True):
with mne.viz.use_browser_backend("matplotlib"):
fig = mags.plot(butterfly=True, proj=proj)
fig.subplots_adjust(top=0.9)
fig.suptitle(f"proj={proj}", size="xx-large", weight="bold")
Using matplotlib as 2D backend.
Using qt as 2D backend.
Using matplotlib as 2D backend.
Using qt as 2D backend.
投影器的其他可视化方法见教程 使用 SSP 修复伪迹。
加载与保存投影器#
除了降低环境噪声，SSP 还可用于其他类型的信号清理。你可能已经注意到，上一张图的磁强计信号中有两个大幅偏转没有被空房间投影器去除——那是被试心跳产生的伪迹。SSP 同样可以用来去除这些伪迹。示例数据中包含用于降低心跳噪声的投影器，它们被保存在与原始数据分开的文件中，可以用 mne.read_proj() 函数加载：
ecg_proj_file = os.path.join(
sample_data_folder, "MEG", "sample", "sample_audvis_ecg-proj.fif"
)
ecg_projs = mne.read_proj(ecg_proj_file)
print(ecg_projs)
Read a total of 6 projection items:
ECG-planar-999--0.200-0.400-PCA-01 (1 x 203) idle
ECG-planar-999--0.200-0.400-PCA-02 (1 x 203) idle
ECG-axial-999--0.200-0.400-PCA-01 (1 x 102) idle
ECG-axial-999--0.200-0.400-PCA-02 (1 x 102) idle
ECG-eeg-999--0.200-0.400-PCA-01 (1 x 59) idle
ECG-eeg-999--0.200-0.400-PCA-02 (1 x 59) idle
[<Projection | ECG-planar-999--0.200-0.400-PCA-01, active : False, n_channels : 203>, <Projection | ECG-planar-999--0.200-0.400-PCA-02, active : False, n_channels : 203>, <Projection | ECG-axial-999--0.200-0.400-PCA-01, active : False, n_channels : 102>, <Projection | ECG-axial-999--0.200-0.400-PCA-02, active : False, n_channels : 102>, <Projection | ECG-eeg-999--0.200-0.400-PCA-01, active : False, n_channels : 59>, <Projection | ECG-eeg-999--0.200-0.400-PCA-02, active : False, n_channels : 59>]
相应地，mne.write_proj() 函数可以把投影器以 .fif 格式保存到磁盘：
mne.write_proj('heartbeat-proj.fif', ecg_projs)
说明
按照约定，MNE-Python 期望投影器保存为以 -proj.fif（或 -proj.fif.gz）结尾的文件名；如果你不遵循这一建议，程序会发出警告。
添加与移除投影器#
上面打印从文件加载的 ecg_projs 列表时，它显示了两个针对梯度计的投影器（前两个，标记为“planar”）、两个针对磁强计的投影器（中间两个，标记为“axial”），以及两个针对 EEG 传感器的投影器（最后两个，标记为“eeg”）。我们可以用 add_proj() 方法把它们添加到 Raw 对象中：
6 projection items deactivated
要移除投影器，有相应的 del_proj() 方法，它会根据投影器在 raw.info['projs'] 列表中的索引将其移除。若要用新的投影器替换现有投影器，可使用 raw.add_proj(ecg_projs, remove_existing=True)。
要查看 ECG 投影器如何影响测量信号，我们可以再次分别绘制应用与不应用投影器时的数据（不过请记住，plot() 方法只是临时应用投影器以便可视化，并不会永久改变底层数据）。下面把前面创建的、仅含空房间 SSP 投影器的 mags 变量，与同时含有空房间和 ECG 投影器的数据进行比较：
mags_ecg = raw.copy().crop(tmax=2).pick(picks="mag")
for data, title in zip([mags, mags_ecg], ["Without", "With"]):
with mne.viz.use_browser_backend("matplotlib"):
fig = data.plot(butterfly=True, proj=True)
fig.subplots_adjust(top=0.9)
fig.suptitle(f"{title} ECG projector", size="xx-large", weight="bold")
Using matplotlib as 2D backend.
Using qt as 2D backend.
Using matplotlib as 2D backend.
Using qt as 2D backend.
投影器在何时被“应用”？#
默认情况下，从 Raw 数据创建分段数据（epoch）时会应用投影器，不过可以通过向 Epochs 构造函数传入 proj=False 来推迟应用。然而，即使投影器尚未被应用，mne.Epochs.get_data() 方法返回的数据仍会像是已应用投影器一样（但 Epochs 对象本身不会改变）。此外，如果数据未预加载（not preloaded），就无法应用投影器。如果数据是内存映射的（即未预加载），可以检查 _projector 属性，看数据载入内存后是否会应用投影器。
最后，在进行逆成像（即使用 mne.minimum_norm.apply_inverse()）时，投影器会被自动应用。在处理 Raw、Epochs 或 Evoked 对象时，也可以通过对象的 apply_proj() 方法手动应用投影器。默认情况下，apply_proj() 会应用附加到该实例上的所有投影器。若只想应用某个特定投影器或子集，可通过 projs= 传入所需的投影器。所选的投影器必须已附加到该实例上。例如，ecg_projs[2:4] 是两个磁强计（“axial”）ECG 投影器，因此我们可以只应用这两个：
raw_ecg_mag = raw.copy().apply_proj(projs=ecg_projs[2:4])
Created an SSP operator (subspace dimension = 2)
2 projection items activated
SSP projectors applied...
这样就无需临时移除其他已附加的投影器。
警告
请记住，投影器一旦被应用就无法撤销，因此在交互式/探索性分析中，最好在应用投影器之前先使用对象的 copy() 方法。
最佳实践#
一般而言，建议在从 Raw 数据创建 Epochs 时应用投影器。这样建议有两个原因：
数据被缩减为仅包含感兴趣的片段（即 epoch）之后再应用投影器，计算开销更小。
如果你对 epoch 使用基于幅值的剔除标准，最好在应用投影器之后根据信号进行剔除，因为投影器可能把某些 epoch 中的噪声降低到可接受的水平（从而增加可用 epoch 的数量，进而提高后续分析的统计功效）。
参考文献#
脚本总运行时间：（0 分钟 3.963 秒）