论文学习：Good practice for conducting and reporting MEG research

部分术语：

Magnetoencephalographic (MEG) 脑磁图

electrooculogram (EOG) 眼电图

electrocardiogram (ECG) 心电图

electromyogram (EMG) 肌动电流图

Magnetic Resonance Imaging (MRI) 磁共振成像

Event-Related Potential (ERP) 事件相关电位

Event-Related Field (ERF) 事件相关磁场

Inter-stimulus intervals (ISI) 刺激间隔

finite impulse response (FIR) 有限脉冲响应

fast Fourier transform (FFT) 快速傅里叶变化

infinite impulse response (IIR) 无限脉冲响应

General Linear Model (GLM) 广义线性模型

论文学习：Good practice for conducting and reporting MEG research

两本作者推荐的刊物：

Hansen et al., 2010

Henson, R.N., Flandin, G., Friston, K.J., Mattout, J., 2010. A parametric empirical Bayesian framework for fMRI-constrained MEG/EEG source reconstruction. Hum. Brain Mapp. 31, 1512–1531.

Picton et al., 2000

Picton, T.W., Bentin, S., Berg, P., Donchin, E., Hillyard, S.A., Johnson, R., Miller, G.A., Ritter, W., Ruchkin, D.S., Rugg, M.D., Taylor, M.J., 2000. Guidelines for using human event-related potentials to study cognition: recording standards and publication criteria. Psychophysiology 37, 127–152.

1 数据获取 (Data acquisition)

1.1 考虑事项

• 对系统设备的考虑
  • 运行前检查其设置，如是否关闭了传感器等
  • 关注系统的运行状态
• 对参与者的考虑
  • 移除任何可能干扰测量的事物
  • 向参与者介绍要做的任务
  • 在进行过程中给予参与者指示
• 对实验设计的考虑
  • 要足以应对各种突发情况，收集足量的数据
• 对常规要求的考虑
  • 足够高的采样率 (sampling rate) 和滤波器通带 (filter pass-band)
  • 记录补充数据 (supplementary data)
  • 等等

1.2 建议和说明

对系统设备的建议

• 系统在实验前后要有两分钟的 “empty room” 的记录时间
  • 也就是说要在实验前后记录2分钟环境噪音
• 为了获得最佳性能，有时需要对单个传感器进行调整
• 准备一本日志

对参与者的建议

• 鉴别不适合的参与者
  • 参与者的任务：深呼吸5-10s，当MEG设备提示时眨眼，每秒张闭嘴一次
• 关于头部定位线圈的放置和头部的数字化
  • 建议在覆盖传感器的同时在参与者的头上也附上定位线圈
  • 避免线圈位置完全对称
  • 建议对头部进行数字化
• 将杜瓦瓶 (dewar) 置于参与者上方
  • 建议参与者尽可能地靠近MEG传感器
  • 建议参与者尽可能地避免头部活动

对实验设计的建议

• 进行简单的光谱分析及存储原始数据和光谱，以供以后比较
• 建议同时采集EOG (水平和垂直方向) 和ECG
• 在实验开始前后记录至少2分钟 “resting state”
  • 记录的时候必须保证MEG系统和所有需要进行实验的设备已经全部准备和安置完毕
• 根据经验，抽样频率应该比感兴趣的最高频率高3-4倍

对常规要求的建议

• 如果对源的定位是计划分析的一部分，则需要获取MEG研究参与者的个体结构/解剖学的磁共振成像 (anatomical MRI)
• 为了实验的准确性，建议通过目视检查进行对齐

1.3 报告撰写

需要在报告中呈现的内容

• 被试者的准备
  • 穿无磁性的衣服
  • 移除金属物品
• 要点和受试者可能的任务
• MEG
  • 系统 (供应商、型号、采集软件的版本)
  • 传感器类型
  • 采样频率
  • 通频带 (采集滤波器的截止频率)
• MRI
  • 系统 (供应商、型号)
  • 序列 (sequence)
  • 立体像素的大小 (vexel size)
• EOG/ECG/EMG
  • 电极的类型、数量、位置和阻抗，放大器 (amplifier) 类型
• 头部形状
  • 用于数字化头部形状的系统
• 头部运动
  • 对头部活动和最大活动的记录
  • 如果不能连续记录，则应说明运行前后头位置的最大差异
• 参与者的位置
  • 说明参与者是坐位还是仰卧位
  • 说明参与者的头部在MEG头盔中的位置
• 外部刺激和记录设备
  • 与MEG触发器相关的状态延迟
  • 描述数据修正过程等

2 数据预处理 (Data preprocessing)

2.1 背景

失真 (artifact) 的种类

• 系统相关的失真 (system related artifacts)
  • 如SQUID jumps，传感器问题
• 外界引起的失真 (external artifacts)
  • 如人体自身的磁场和外界的噪音
• 生理导致的失真 (physiological artifacts)
  • 如眼动、眨眼和肌肉活动等

两种处理失真的常用方法

• **目视检查 (visual inspection) **以识别被污染的数据段
• 信号处理技术 (signal-processing techniques) 以减少不是由大脑活动引起的失真部分

2.2 建议和说明

• 为了识别失真和检查失真修正方法的准确性，强烈建议对数据进行目视检查

• 处理系统相关的失真

  • 必须要在进行下一步分析之前排除平坦且嘈杂的通道 (channel)

  • 可以通过对线路频率及其谐波处拟合和减去sine和cosine函数 (Mitra and Pesaran, 1999) ，或使用ICA (Barbati et al., 2004) 来减少线路干扰 (line interference)

    Line interference can be reduced by fitting and subtracting sine and cosine functions (Mitra and Pesaran, 1999) at the line frequency and its harmonics, or with ICA (Barbati et al., 2004).

  关于第二小点我不是很理解，故附上了原文，并在这里给出作者参考的文献：

  Mitra, P.P., Pesaran, B., 1999. Analysis of dynamic brain imaging data. Biophys. J. 76, 691–708.

  Barbati, G., Porcaro, C., Zappasodi, F., Rossini, P.M., Tecchio, F., 2004. Optimisation of an independent component analysis approach for artifact identification and removal in magnetoencephalographic signals. Clin. Neurophysiol. 115, 1220–1232.

• 排除环境噪音

  • 可以应用 ICA , PCA , 基于回归的方法和 SSS 等来排除环境噪音

• 修正生理导致的失真

  • 可以使用信号处理技术来抑制 EOG/ECG 失真
  • 修正头部运动

• 使用自动或半自动的方法 (automatic or semi-automatic method) 或者目视检查来排除失真

• 注意每一步骤的顺序

  1. 排除系统相关的失真
  2. 排除线路噪音 (line noise)
  3. 排除环境噪音
  4. 排除生理噪音

• 其他

  • 当想要进行后续的源重建时，正演模型 (forward model) 中还需要考虑使用子空间投影或回归技术从数据中去除空间分量或减去加权参考信号
  • 多通道脑磁图信号中的低频调制可能是由呼吸相关的失真引起的
    • 可以使用应用于原始、连续记录的高通滤波器(截止频率约为1hz)来缓解
  • 刺激装置，尤其是神经电刺激器，会造成类似脉冲的失真
    • 使用PCA事件相关的修正 (PCA event-related correction)
    • 如果失真发生在大脑活动之前，可以切除被污染的片段并插入省略的样本
  • 失真可能有属于生理相关频带范围内的频率成分。肌肉失真 (muscle artifacts) (例如，来自颈部或颌部肌肉) 以高频振荡爆发的形式出现在MEG信号中

2.3 报告撰写

需要在报告中呈现的内容

• 每个预处理步骤的所有参数和这些步骤的顺序
• 说明坏的MEG传感器数量和是否有用插值代替的数据
• 滤波器的类型、顺序、截止频率和应用方向
• ICA算法、输入数据、估计的组件数量、被排除的组件数量和选择这些组件的标准
• 用于后续分析的纯净试验 (clean trials) 或数据片段的数量
• 排除数据的标准 (尤其是基于目视检测时)
• 指出是否使用头部运动补偿算法以及使用哪种算法
• 对于自适应干扰抑制方法如SSS的时间扩展(Taulu and Simola, 2006)，指出相关参数
  • 对于tSSS，指出相关窗口的持续时间和相关极限

3 ERFs分析 (Analysis of event-related fields)

3.1 背景

对事件的分析，通常采用的方法是通过一个刺激或者行为诱发大脑反应，并对每个事件的数据取平均值进行分析。

• 对EEG的分析方法：Event-Related Potential (ERP) 事件相关电位
• 对MEG的分析方法：Event-Related Field (ERF) 事件相关磁场

然而，这种方法是有局限的，将环境活动简单归为噪音的方法如今已不再适用，大量实验表明大脑对事件的反应是多因素相互作用的结果。比如，不仅仅是感官输入，还有感官与大脑动态的变化等等因素也会影响大脑的反应。

3.2 建议和说明

• 将实验次数减到最少

  • MEG响应存在实验和对象间的可变性
  • 总的来说，实验的数量的确定需要有一个权衡。
    • 其最小值取决于生理相关的大脑区域，如细胞密度、种类和位置。
    • 其最大值取决于受试者能够按需要的标准进行任务的时长 (需要保持稳定的头部位置和避免眨眼等)

• 在事件相关的范例中，在刺激间隔 (ISI) 中引入随机元素

  • 减少期望对诱发反应的影响
  • 避免混叠或累积的周期性干扰(如线路噪声)对平均值的影响

• 通过EOG和EMG测量来排除高频振荡和肌肉失真 (muscle artifact) 污染

• 应用滤波器来提取感兴趣的预定义频带内的信号

  • 建议对较长的数据集应用数字滤波器

• 基线矫正 (baseline correction)

  • 对连续原始MEG图使用低截止频率的高通滤波器可以有效地消除偏移，但不建议对短的诱发数据使用
  • 从每个通道的信号中减去特定时间间隔内数据的平均值

• 对传感器级数据进行分组分析

  • 在没有使用运动补偿技术等的情况下，不建议在传感器级别进行大量 (跨科目) 或间歇的平均
  • 即使没有头部位置校正，原始的传感器级的分析也是对数据的无假设分析，通常会在保守(假阴性)方面出错

  这里我也不是很理解

3.3 报告撰写

在报告中应该呈现的内容

• 进行的刺激的明确时间表

  • 包括ISIs，刺激类型和持续时长

• 每个实验条件和所有受试者使用的epoch数量

• 所用的基线矫正技术和它的参数

• 滤波器

  • 包括滤波器类型，频率截止，顺序/核心长度

• 控制肌肉失真对高频(>40 Hz)成分的可能污染

• 用于估计峰值延迟、振幅和此分析的时间窗口的标准 (如果使用的话)

• 如果合并了不同类型的传感器(如磁力计和梯度仪)，需要指定使用的方法

• 如果进行了降维 (PCA，ICA，传感器选择) ，说明选择组件/传感器的步骤和标准

4 光谱分析 (Spectral analysis)

4.1 背景

光谱分析的目的是将信号转变为时间频域 (time frequency domain)

• 参数方法 (parametric method)
  • AR-models
• 非参数方法 (non-parametric method)
  • 贝叶斯变化
  • wavelet 变化
  • Hilbert 变化

4.2 建议和说明

• 可以使用补零来实现光谱插值，以便能够对实验中不同长度的光谱表示进行平均

• 为了避免零填充可能带来的光谱泄漏效应，建议在光谱分析之前进行基于试验的平均减法

• 建议对所有频率使用相同长度的数据段，或分别分析频率范围

• 为了达到最佳的灵敏度，建议在更高的频率上增加频率平滑

  • 电生理数据中不同频率分量的带宽随频率的增加而增加，且峰值频率也不同
  • 30Hz以下 几乎不平滑
  • 30-50Hz 普通的平滑
  • 50Hz以上 更大程度的平滑

• MEG 数据很难确保一定准备，受试者可能头动等，为了方便解释光谱的数量，通常需要对一个特定频率进行矫正

  • 方法一：使用一个基线时间窗或基线任务
  • 方法二：基于功率水平的差异，或功率对数的差异，这相当于计算相对变化。此外，还可以考虑考虑基线中的方差信息，并计算z-score。最后，可以计算相对功率，其中每个频带的功率由总功率归一化得到

• 对用于基线归一化的时间窗的选择

  • 刺激前选择基线：避免将基线演唱到刺激开始时
  • 基线起始时间应避免包括前一次刺激的结束时间

• 接近刺激开始的基线偏移可能会降低时频图中与刺激相关的低频分量的幅度

• 对限带效应峰值的频率的考虑需要将固有频率分辨率也加入考虑

  The interpretation of the frequency of the peak of a band-limited effect needs to take the intrinsic frequency resolution into account.

• 小波 (wavelet) 的频率特性 (frequency specificity) 需要为算法的精确性让步

• 时域中的短暂失真或早或晚会对光谱估计产生影响

• 建议使用至少为最低频率小波窗口长度一半的试验数据进行填充来避免边缘数据的丢失

• 在解释光谱分析的结果时需要考虑MEG传感器的类型

4.3 报告撰写

在报告中需要有所呈现的内容

普通要求 (general)

• 使用的时间窗口 (time-window) 的长度
• 重叠窗口的数量
• 如果应用了基线估计的标准化，说明时间区间和标准化的流程
• 展现结果的数据必须清晰地描述所用的变化 (如log-transform)
  • 表格的每一个值、单元和刻度都需要提供名字
  • 图片尽可能展示方差信息

特定方法的要求

• 非参数方法 (non-parametric methods)
  • 具体说明插零的数量
  • 具体说明用于控制频谱泄漏的锥形策略 (tapering strategy)
• 参数方法 (parametric methods)
  • 用于估计AR模型的算法
  • AR模型的顺序和用来决定这个顺序的方法
  • 详细说明在估计AR模型之前的预处理步骤
    • 总平均值的去除 (??? 不是很理解)
    • 数据标准化

5 源重建 (Source reconstruction)

5.1 背景

源定位的目的是在传感器级上重构信号的神经源，所有的源重建方法都需要对正演模型 (forward model) 进行描述，即详细说明在已知位置和方向上具有统一的强度的神经源产生的磁场分布的模型。

在定义了体积导体模型之后，通常需要指定一个源模型，该模型用有限数量的参数逼近潜在的连续电流密度分布。

5.2 建议和说明

• 在扫描过程中准确的配准和受试者的头部保持不动非常有利于源重建的质量

• 在准备受试者时，特别注意头部定位线圈的放置和基准位置的说明；

  • 同时要让受试者保持一个舒适的位置，这样有利于其头部保持不动

• 源本身的本地化 / 提取源时间序列可以当作源定位过程的一部分或者单独的步骤

  这里不确定主语是哪一个所以两个都写了

• 在源级重建时间序列的过程中，要确保用于估计源时间序列的线性加权不仅适用于其位置的活动，而且它们是无偏的

  • 在比较两种情况时，为了得到一个无偏的线性映射，需要将这两个情况都包含到反演过程中去
  • 我们应该只比较在反演参数估计的同一时频窗口内的神经活动的估计

• 当将源重构约束为一组 ROIs 或预先定义的源模型时，并非来自建模区域的大脑活动可能会被错误地投射到这些源上

• 在源估计中需要考虑降维

  • 消除数据的自由度
  • 减小传感器协方差矩阵的秩

• 空间范围估计通常需要使用不同源模型的特定算法

5.3 报告撰写

在报告中需要有所呈现的内容

普通要求

• 头部的体积导体模型的细节和用于导程电场 (lead field) 计算的算法
• 选择一个特定的源分析方法的基本原理
• 源重建算法的细节和其所有可调参数
  • 时间间隔
  • 正则化的选择和分析
  • 基本原理
• 所使用源定位技术的缺点和对结果的影响
• 如果坐标呗连接到大脑结构上，需要说明lookup表的来源

特定方法的要求

波束形成 (Beamforming)

• 用于波束形成器输出的协方差矩阵和度量的条件数量
• 如果对重建的源图像进行了噪声归一化，应提供如何计算噪声水平的详细信息
• 关于已扫描的偶极网格 (dipole grid) 的细节
  • 规则体积网格的体素大小
  • 皮层表的节点数
• 指定所有协方差矩阵计算的细节
  • 应用于数据的时间窗口、过滤器和其他预处理

分布式源模型 (Distributed source models)

• 关于已扫描的偶极网格的细节
  • 规则体积网格的体素大小
  • 皮层薄片的节点数量
  • 偶极位向的约束

偶极子定位 (Dipole fitting)

• 研究者的选择对用偶极拟合方法得到的结果很重要，需要详细说明
  • 选择原因
  • 偶极子的数量
  • 所使用的时间窗口 (单个延迟，多个延迟)
  • 所使用的偶极子模型 (移动、旋转或固定偶极子)
  • 是否某些参数是固定的而其他参数是优化的
• 所使用的成本函数与优化算法

6 连通性分析 (Connectivity analysis)

6.1 背景

MEG连通性分析的目的是研究大脑区域之间的信息流和相互作用。对于连通性分析的方法孰好孰坏，目前尚未达成共识。在估计MEG数据的连通性时电磁场的扩散 (与电流有关的磁场可以无限延伸的现象) 为核心问题，这意味着我们实际测出来的传感器级生理电信号是会相互干扰的。

避免因磁场扩散而对生理相互作用作出错误解释的一种方法是量化只能由非瞬时相互作用解释的连通性，减少磁场扩散影响的推荐方法是在底层源级而不是在传感器级分析连接性。

• 两步方法 (two-step approaches)
  • 首先执行源时间过程的估计，然后将这些位置的波束形成滤波器应用到传感器信号中，以获得源时间程。
• 一步方法 (one-step approaches)
  • 基于一个完整的用于MEG信号产生的模型。

6.2 建议和说明

• 连通性分析应该在源级 (source level) 进行

  • 相比传感器级更简单且更容易与其他研究结合

  • 减轻了电磁场扩散引起的混合问题所带来的解释困难

    补充优点：源数据的数量会比传感器信号要少得多

• 如果连通性分析仍然在传感器级别进行，那么应该注意，磁场扩散的影响取决于MEG传感器的类型

• 建议通过对比两个实验的条件 (或刺激前和刺激后时间窗口，或一个条件和替代数据) 来进行源空间的连通性分析

  • 进一步减少磁场扩散的影响
  • 如果进行了对比，建议匹配两种条件下的数据样本总数

• 当使用使用限带振荡相位 ( a band-limited oscillation) 估计的连通性度量时，建议仔细分析,选择足够窄的频带

• 估计的质量取决于数据的信噪比 (signal-to-noise ratio) ，噪声高的信号的相位估计精度较低

  • 报告的相位同步差应该伴随着信号幅值的报告

• 与大脑活动的任何其他指标一样，连通性测量也需要进行统计测试

  • 大多数连通性测试是片面的估计量，其数据质量取决于数据数量和信噪比
  • 测试用的人造数据可能会破坏连通性
  • 因为要研究大量的潜在连接，对于连通性分析而言，大量比较的问题通常更为严重

• 如果使用了不适当的替代数据来在单个情况下测试连接性的重要性，可能会产生虚假的结果

• 去除空间分量或减去加权参考信号会影响MEG主传感器之间的相关性，这可能会混淆对后续连通性结果的解释

• 当连通性的估计涉及到一个自回归模型或一组嵌入参数用于状态空间重构时，模型或参数的错误说明可能导致假阳性结果和反向连接

• 用相干虚部评估耦合有一定的局限性

  • 其对真正的生理零相位耦合不敏感
  • 比较时可能会出现问题，解释连贯的虚部变化并不是无关紧要的
  • 两个或多个具有不同方向的独立源可能过于接近以至于无法将它们与具有多个方向成分的单一源区分开来
  • 来自心脏的电信号导致了非零的虚相干

• 未观测到的源可能会影响对观测到的连通性模式的解释

• 当使用两步方法 (tow-step approaches) 来估计源级的连通性时，只分析所使用的重建方法的源时间过程的特征

• 在定向测量中，对观测到的依赖关系与零滞后依赖关系 (zero-lag dependencies) 的检验可以表明场扩散的存在

6.3 报告撰写

需要在报告中呈现的内容：

普通要求

• 说明预处理步骤的细节
  • 滤波器 (因为过滤器的应用可以改变连通性分析的结果)
• 描述用于处理诸如心跳信号或眨眼等失真的程序
• 报告用于估计连通性的数据量
• 考虑电磁场扩散的潜在影响下讨论对结果的解释
  • 这对于传感器级的分析尤为重要
• 讨论对估计的连通值在不同环境下的解释 (不同的信噪比和源功率) 并报告不同情况间的功率变化。

特定方法的要求

• 基于自回归模型的连通性测度 (Granger因果关系和派生测度，以及传递熵)

  Autoregressive model based connectivity measures (Granger causality and derived measures, and transfer entropy)

  • 详细的自回归模型拟合程序 (特别是选择的模型顺序)

• 动态因果模型 (Dynamic causal modeling)

  • 建议选择考虑评估的特定模型

    • 为什么在模型空间中包含特定的源区域
    • 为什么考虑特定的连通方案

    参考文献：

    Stephan, K.E., Penny, W.D., Moran, R.J., den Ouden, H.E., Daunizeau, J., Friston, K.J., 2010. Ten simple rules for dynamic causal modeling. Neuroimage 49, 3099–3109.

7 统计 (Statistics)

7.1 背景

统计分析 (statistical analysis) 的目的是量化从数据中估计的特征的可靠性 (不无视噪音和不同实验对象的差异) 。基本上，先进的MEG数据分析管道中的许多步骤都带有可能需要量化的不确定性。在这里，我们将主要考虑与假设检验相关的统计问题，这些问题解决了关于激发数据 (motivated data) 获取和实验设计的神经科学问题。

几种提出关键问题的测试

• 参数方法或无参方法
  • 参数方法依赖于关于数据特定参数的概率分布的假设
  • 无参方法使用重采样策略从数据中估计概率分布或置信区间
• 单变量或多变量
  • 多变量测试 (multivariate test) ：在单个测试中对多个特征 (变量) 进行一次测试
  • 大量单变量测试 (mass univariate tests) ：使用单独的单个测试对每个变量进行测试
• 基于经典推理或贝叶斯推理
  • 经典推理于检验数据的一个正式假设 (formal hypothesis) 相关
  • 贝叶斯统计可以提供具有一定规模的效应的概率，但需要先验概率分布

需要判断统计推断需要在传感器还是源空间中进行

• 传感器级分析是相对自由和稳健的假设，但假设(反演)是正确的

• 源级分析将更加敏感。

  在某些特定情况下，在源空间中直接执行统计信息可能会更好

  • 当传感器数据没有在多个重复中对齐时
  • 当源的方向在不同的对象之间有很大差异时
  • 对于多模式的数据集 (mutimodal datasets) 而言
  • 当源空间中存在先验时
  • 当需要对大脑的特定区域进行推断时

7.2 建议和说明

• 提前决定感兴趣的关键实验比较

  • 避免在决定测试什么比较之前查看数据的问题

• 使用参数或非参数统计的决定应该基于对数据分布的检查或知识

  • 除非有足够的重复，否则对分布假设的正式检验是不敏感的，所以重复很少时更建议使用非参数方法

• 当基本假设得到满足(或数据被转换以满足假设)时，参数方法总是比非参数方法更敏感

  • 当与这些假设存在偏差时，非参数方法将更加稳健

• RFT和MCP的非参数方法都提供了对测试统计数据的特定特性进行推断的能力

• 提前决定所要探索的科学问题是否与时空模式的识别或效应在空间/时间上的定位和映射有关

• 在单变量框架内，多个变量仍然可以作为建模的因素

• 从理论/经验上证明关于数据分布的假设，同时在使用任何数据测试之前检车数据值的分布

• 考虑多重比较——跨空间、时间、频率以及实验设计中的多重对比，确保任何校正多重比较的方法的假设都是满足的

• 由于传感器和源空间分析的灵敏度不同，有时会出现一种特殊效果在其中一种分析中很显著而在另一种分析中则不然

• 预先限制统计检验的范围是提高其灵敏度的一个有效且强力的方法

  • 只对没有预先假设的数据维度使用多个测试，并将这些测试的数量限制到最低限度

• 明确推理类型，不管基于FWE，FDR，或随机/固定的对参与者的效应分析

• 源空间的质量单变量方法 (mass univariate method) 引起了与MCP处理方法有关的问题

  • 诸如RFT等方法的假设可能会被违反，使用RFT对经典p值的修正将过于保守，建议使用无参的MCP方法或者贝叶斯推理

7.3 报告撰写

在报告中需要呈现的内容

• 描述使用的完整统计模型

  • 对于 (参数化的) 广义线性模型 (GLM)
    • 使用的设计矩阵
    • 自由度
    • 纠正误差中的任何相关性
  • 对于非参数测试
    • 所使用的技术类型(随机化、排列、自举)
    • 重新取样的数量

• 说明先验比较和后验比较，以及多次比较的纠正方法

• 影响大小及其统计

  • 诸如平均子等的参数估计