OPM 脑磁数据实操：从原始.mat 文件到 MNE Raw 与 BIDS 标准化全流程

如何从一个 OPM-MEG 原始数据文件夹出发，逐步拆解 MATLAB .mat 字段，构造 MNE Raw 对象，并把数据整理成可验证的 BIDS 结构
Getting started with QuSpin OPM data

Gitee项目源码 multimodal-bids-eeg-meg-squid

1. 先看

在这个项目里，EEG、SQUID-MEG、OPM-MEG 三类数据的原始格式完全不同。

EEG 是 BrainVision 三件套：

.vhdr / .vmrk / .eeg

SQUID 是 Ricoh/KIT 系统文件：

.con

而 OPM 数据是 MATLAB .mat：

run01.meg.mat
run01.eeg.mat

这个数据结构或许有规范，但是还不够规范
数据矩阵、传感器位置、传感器方向、trigger、采样率、通道名、MRI 信息都散落在不同字段里。

所以 OPM 转换的核心问题：

我怎么确认哪些字段是数据？
哪些字段是通道信息？
哪些字段是传感器几何？
哪些字段应该进 MNE Raw？
哪些字段应该写入 BIDS sidecar？

以OPM中听觉任务为例补充一下实操过程

2. 再细细看

原始数据是这样的结构：

OPM-SQUID-EEG-dataset/
  002/
    OPM/
      Auditory/
        run01.meg.mat
        run01.eeg.mat
        run02.meg.mat
        run02.eeg.mat
      Motor/
      Rest/
      Somatosensory/
  005/
  006/
  093/
  095/
  EMPTY/
    OPM/
      run01.meg.mat

这里第一眼要注意几件事

第一，run01.meg.mat 和 run01.eeg.mat 是成对出现的。meg.mat 是 OPM-MEG 主数据，eeg.mat 在数据集的介绍文档中有提到，是眼电通道

• 记得细细查找数据集是否有说明书

第二，同一个任务有多个 run。例如 auditory 有 run01、run02。这些不能直接合并成一个文件，而应该在 BIDS 里保留为不同的 run-01、run-02

第三，EMPTY/OPM/run01.meg.mat 是空房间数据，同样在文档中有介绍

我最终希望整理成这样的 BIDS 结构：

The Brain Imaging Data Structure

bids_dataset/
  sub-002/
    ses-01/
      meg/
        sub-002_ses-01_task-auditory_acq-opm_run-01_meg.fif
        sub-002_ses-01_task-auditory_acq-opm_run-01_meg.json
        sub-002_ses-01_task-auditory_acq-opm_run-01_channels.tsv
        sub-002_ses-01_task-auditory_acq-opm_run-01_events.tsv
        sub-002_ses-01_task-auditory_acq-opm_run-01_coordsystem.json

  sub-emptyroom/
    ses-01/
      meg/
        sub-emptyroom_ses-01_task-noise_acq-opm_run-01_meg.fif

这里有一个关键设计：OPM 和 SQUID 都属于 MEG，所以都放在 meg/ 目录下，用 acq-opm 和 acq-squid 区分设备。

3. 细看 .mat 字段

拿到 .mat 文件后，需要先写一个检查脚本（AI🏃‍跑即可）：

inspect_mat.py

它的作用是递归输出 MATLAB 文件里的变量、嵌套 struct、数组形状、dtype 和数值预览。

使用方式：

python inspect_mat.py "...\093\OPM\Auditory\run01.meg.mat"

一个 OPM .meg.mat 文件里可以看到这些顶层变量：

meg_parm
bexp
refmg
bexp_ext
CoordType
MEGinfo
Measurement
pick
Qpick
ref_pick
ref_Qpick
PositionFile

或者使用MATLAB打开查看，哦这才是聪明的人类应该做的

这些字段大致可以分为四类。

3.1 主数据字段

bexp: shape=(30, 854000)

显然这是 OPM-MEG 主数据矩阵，
这说明有 30 个 OPM 通道，854000 个采样点。

如果采样率是 2000 Hz，则记录时长约为：

854000 / 2000 = 427 s

3.2 trigger 字段

bexp_ext: shape=(1,854000)

这个字段一般是额外通道，比如 trigger pulse。
也就是说，它和主数据有相同的采样点数，但只有一个额外通道。

在构造 MNE Raw 时， 要设置为 stim 类型，作为触发标志

3.3 通道和采样信息

核心字段是：

MEGinfo

里面包含：

SampleFreq
Nsample
Nchannel
MEGch_name
ChannelInfo
ExtraChannelInfo
Trial
MEG_ID

这些字段决定了 MNE Info 对象怎么建。

例如：

SampleFreq = 2000
Nchannel = 30
MEGch_name = M001_y, M001_z, ...
ExtraChannelInfo.Channel_name = trigger_pulse

3.4 传感器位置和方向

OPM 转换最重要的几何字段是：

pick
Qpick
CoordType
PositionFile

其中：

pick  -> 每个 OPM 通道的位置
Qpick -> 每个 OPM 通道的敏感轴方向

这两个字段决定了每个通道在 MNE 里的 loc 应该怎么写。

如果只把数据矩阵放进 Raw，而不处理 pick/Qpick，后续就没法做空间分析、传感器图、源定位了

4. 字段拆完以后，怎么构造 MNE Raw

OPM .mat 转换的核心逻辑可以概括成：

MAT 字段
  ↓
数据矩阵 bexp
  ↓
通道名 / 采样率 MEGinfo
  ↓
传感器位置方向 pick / Qpick
  ↓
mne.create_info
  ↓
mne.io.RawArray
  ↓
写出 FIF

关键代码逻辑大致如下。

4.1 读取采样率、通道名、通道数

def meg_info_from_mat(mat):
    meg_info = mat["MEGinfo"]
    sfreq = float(field(meg_info, "SampleFreq"))
    n_channels = int(field(meg_info, "Nchannel"))
    ch_names = [str(name) for name in as_list(field(meg_info, "MEGch_name"))]
    device = str(field(meg_info, "device"))
    return meg_info, sfreq, n_channels, ch_names, device

这里要注意：不要依赖 MATLAB struct 的字段顺序，要按字段名读取。否则不同文件或不同 MATLAB 保存版本可能会让脚本不可用

4.2 检查数据矩阵形状

meg_data = np.asarray(mat["bexp"], dtype=float)

if meg_data.shape[0] != n_channels:
    raise ValueError("bexp channel count does not match MEGinfo.Nchannel")

if meg_data.shape[1] != int(field(meg_info, "Nsample")):
    raise ValueError("bexp sample count does not match MEGinfo.Nsample")

让AI审查一下，顺手加几个检查点

数据矩阵通道数 == MEGinfo.Nchannel
数据矩阵采样点数 == MEGinfo.Nsample

4.3 把 trigger channel 追加到 Raw

bexp_ext作为额外通道拼到数据矩阵后面：

all_data = meg_data
all_names = ch_names
all_types = ["mag"] * n_channels

bexp_ext = mat.get("bexp_ext")
if bexp_ext is not None and np.size(bexp_ext) > 0:
    bexp_ext = np.asarray(bexp_ext, dtype=float).reshape(1, -1)
    all_data = np.vstack([all_data, bexp_ext])
    all_names = all_names + ["trigger_pulse"]
    all_types = all_types + ["stim"]

这样最后 MNE 里不仅有 OPM-MEG 通道，也有 trigger 通道，后续要用来提取事件：

mne.find_events(raw, stim_channel="trigger_pulse")

4.4 创建 MNE Info 和 RawArray

info = mne.create_info(
    ch_names=all_names,
    sfreq=sfreq,
    ch_types=all_types,
)

raw = mne.io.RawArray(all_data, info)

到这一步，只是得到了一个“有数据、有通道名、有采样率”的 Raw，还不够。OPM 的传感器位置和方向还需要写进 raw.info["chs"][i]["loc"]。

4.5 写入 OPM 传感器位置和方向

MNE 每个通道的 loc 是一个 12 维数组。对 OPM 来说，至少要写入：

loc[0:3]  -> 传感器位置
loc[3:6]  -> X轴方向 (ex)
loc[6:9]  -> y轴方向 (ey)
loc[9:12] -> z轴方向 (ez) 敏感轴向

loc 的核心就是用一个位置和三个互相垂直的方向向量，完整地定义了一个传感器在空间中的六自由度姿态
Three degrees of freedom for position, and three for orientation）

def make_opm_loc(position, sensitive_axis):
    """将 OPM 的位置和敏感轴方向转换为 MNE 的 12 维 loc 数组"""
    loc = np.zeros(12)

    # 位置：直接写入
    loc[0:3] = np.asarray(position, dtype=float)

    # 敏感轴方向 (Z轴)：归一化
    ez = np.asarray(sensitive_axis, dtype=float)
    ez /= np.linalg.norm(ez)

    # 构造X轴 ：用Z轴与世界Z轴叉乘，得到水平切向
    ex = np.cross(ez, [0, 0, 1])
    # 如果敏感轴恰好指向正上方或正下方，叉乘结果为0，改用Y轴做回退
    if np.linalg.norm(ex) < 1e-6:
        ex = np.cross(ez, [0, 1, 0])
    ex /= np.linalg.norm(ex)

    # 构造Y轴：用Z轴与X轴叉乘，补全右手系
    ey = np.cross(ez, ex)

    loc[3:6] = ex
    loc[6:9] = ey
    loc[9:12] = ez
    return loc

然后对每个 OPM 通道写入：

for idx in range(n_channels):
    raw.info["chs"][idx]["loc"] = make_opm_loc(pick[idx], Qpick[idx])

Raw 不只是一个时间序列矩阵，而是带有传感器空间信息的 MEG 数据对象。

踩坑记录：coil_type 未填充导致 HFC 无效

我一开始对 coil_type 和 loc 的理解也不够清楚。后来发现，这两个字段会直接影响 OPM 后续空间去噪方法，比如 HFC

在完成上述 loc 填充后，我跑了一下 HFC（Homogeneous Field Correction）

结果：完全没有效果。数据没有任何改善，HFC 根本没被执行

问题定位
经过排查，问题出在 coil_type 字段上。

raw.info['chs'][idx] 中有一个关键字段叫 coil_type。它告诉 MNE 这个通道是什么类型的传感器——是磁力计、梯度计，还是别的什么

mne.preprocessing.compute_proj_hfc 在执行时，会根据 coil_type 来选择正确的物理模型和线圈定义。

而我当时用 mne.create_info 创建 OPM 通道时，coil_type 默认是 0（未定义）。
HFC 拿到一个 coil_type=0 的通道，不知道该按什么模型处理，于是静默跳过了——不报错，也不生效。

但 coil_type 不是 .mat 里直接有一个字段叫“线圈类型”然后拿来填。它是这么来的：

MEGinfo.device = QZFM
数据集/文档说明是 QuSpin OPM
MNE 里有对应的 FIFFV_COIL_QUSPIN_ZFOPM_MAG

info["chs"][i]["loc"] = make_opm_loc(pick_pos[i], pick_ori[i])
info["chs"][i]["coord_frame"] = FIFF.FIFFV_COORD_HEAD
info["chs"][i]["coil_type"] = FIFF.FIFFV_COIL_QUSPIN_ZFOPM_MAG

5. 哪些信息放进 Raw，哪些信息放进 sidecar

Converting an example EEG dataset for sharing in BIDS

这是做 BIDS 时很容易混乱的地方

MNE 需要直接使用的信息 -> 放进 Raw / Info
说明性、追踪性、原始来源信息 -> 放进 JSON sidecar
逐通道表格信息 -> 放进 channels.tsv
事件信息 -> 放进 events.tsv
坐标系说明 -> 放进 coordsystem.json

5.1 放进 Raw / Info 的信息

数据矩阵
通道名
通道类型
采样率
line frequency
传感器 loc
bad channel 标记

5.2 放进 _meg.json 的信息

TaskName
Manufacturer
SamplingFrequency
PowerLineFrequency
RecordingDuration
MEGChannelCount
TriggerChannelCount
SourceFile
PositionFile
MEG_ID
MRI_ID
Vcenter
Vradius
AssociatedEmptyRoom

例如：

{
  "TaskName": "auditory",
  "Manufacturer": "QuSpin",
  "SamplingFrequency": 2000,
  "RecordingType": "continuous",
  "SourceFile": ".../run01.meg.mat",
  "AssociatedEmptyRoom": "sub-emptyroom/ses-01/meg/..."
}

5.3 放进 _channels.tsv 的信息

每个通道一行：

name
type
units
sampling_frequency
status
status_description

OPM 主通道写成：

MEGMAG

trigger 通道写成：

TRIG

5.4 放进 _events.tsv 的信息

事件表来自 trigger 通道：

onset
duration
trial_type
value
sample

其中：

onset = sample / sfreq

这样后续 epoch 时就可以直接根据 events 表来做。

5.5 放进 _coordsystem.json 的信息

OPM 的坐标系来自原始 MAT：

CoordType = SPM_Right_m

所以我写入：

{
  "MEGCoordinateSystem": "Other",
  "MEGCoordinateUnits": "m",
  "MEGCoordinateSystemDescription": "SPM_Right_m"
}

6. Empty-room 怎么处理

MEG 数据里 empty-room 很重要，尤其是 OPM。它可以用于：

环境噪声评估
noise covariance
去噪前后比较
后续 HFC/回归方法参考

所以不要将 empty-room 混进某个被试，而是单独放在：

sub-emptyroom/ses-01/meg/

例如：

sub-emptyroom_ses-01_task-noise_acq-opm_run-01_meg.fif

然后在每个实验 run 的 _meg.json 里记录：

{
  "AssociatedEmptyRoom": "sub-emptyroom/ses-01/meg/sub-emptyroom_ses-01_task-noise_acq-opm_run-01_meg.fif"
}

这样后续要计算 empty-room covariance，就能追踪到对应噪声文件。

7. 为什么不能把多个 run 合并

同一个被试可能做了两次 auditory，比如：

run01.meg.mat
run02.meg.mat

我的处理方式是保留为两个 run：

run-01
run-02

而不是合并。

原因有三

第一，每个 run 可能有不同的噪声状态、头动状态、坏道状态。

第二，后续 QC 时可以按 run 比较质量。

第三，BIDS 本身就鼓励用 run 区分同一任务的重复采集。

所以合理结构是：

sub-093_ses-01_task-auditory_acq-opm_run-01_meg.fif
sub-093_ses-01_task-auditory_acq-opm_run-02_meg.fif

而不是：

sub-093_ses-01_task-auditory_acq-opm_meg.fif

8. 最终转换流程

最后 OPM 转换整理成批处理脚本：

batch_opm_to_bids.py

整体流程是：

扫描 source dataset
  ↓
找到每个 subject/task/run 的 runXX.meg.mat
  ↓
读取 MAT 字段
  ↓
验证 bexp / MEGinfo shape 是否一致
  ↓
构造 MNE RawArray
  ↓
写入 OPM 传感器 loc
  ↓
追加 trigger 通道
  ↓
写出 BIDS FIF
  ↓
写 meg.json / channels.tsv / events.tsv / coordsystem.json
  ↓
记录 conversion log

转换完成后，OPM 部分的结果是：

OPM experimental runs: 32
OPM empty-room runs: 1
Total OPM FIF files: 33

9. 验证

写出 BIDS 文件不等于转换成功。我还需要验证：

FIF 能不能被 MNE 读取
channels.tsv 行数是否等于 raw channel 数
events.tsv 是否存在
coordsystem.json 是否存在
AssociatedEmptyRoom 是否指向真实文件

所以拷打AI：

validate_opm_bids.py

验证结果：

[1/5] bids-validator: PASS
[2/5] mne.read_raw_bids: 33/33 OK
[3/5] channels.tsv 行数 == raw.nchan: 33/33 OK
[4/5] events.tsv 存在且 onset 在 [0, duration] 内: 32/32 OK
[5/5] AssociatedEmptyRoom 指向真实文件: 32/32 OK

再用 mne_bids.read_raw_bids() 做综合验证，确认这些文件不只是能被 read_raw_fif() 打开，而是真的能作为 BIDS 数据进入后续 pipeline。

总结

一份 OPM 数据从是有格式到 开源标准 的链路
不同设备基本逻辑都是这样，往上套一套，AI🏃‍♀️跑一跑就行咧