引言

在磁共振成像（MRI）系统中，谱仪是主要的核心部件，对成像功能与质量有着至关重要的影响.它的主要功能是控制成像系统、完成脉冲序列的运行、产生梯度信号、发射射频信号，以及接收和处理MRI回波信号.目前，国内外很多实验室都在研制谱仪^[1-8].梯度是成像的关键环节，丰富多样的脉冲序列的区别主要表现在谱仪所发出的三路梯度波形上，其在时间与幅度上的细小偏差都会对磁共振图像造成显著影响.通常，在谱仪研制与序列调试中，通过观察示波器波形或重建图像来判断梯度问题^{[9, 10]}.然而，受显示效果以及存储深度的影响，示波器观察梯度波形难以捕捉细节.另一方面，基于重建图像的判断需要丰富的经验，某些情况下难以判断.因此，本文提出通过采集磁共振梯度波形信号，进行波形脉冲特征点提取的方法，以分析谱仪序列运行中梯度波形是否符合设计要求，从而快速判断谱仪梯度硬件电路与脉冲序列是否存在问题.

当前，随着电子技术与器件的发展，高速采集与数据处理技术在硬件及软件等方面已逐步成熟. National Instructions公司（NI）的LabVIEW软件工具在数据采集方面具有强大的优势，在通信、遥感、生物医疗等领域都有很广泛的应用^[11-14].

采用虚拟仪器LabVIEW软件控制高速采集卡DAQ-2005设计实现多路采集系统，对谱仪的梯度输出进行采集.通过对波形数据进行直方图统计、滤波、差分计算等分析，提取出波形的特征点，这些特征点包含时间与幅度信息.利用实验室自主研制的谱仪对该方法进行了多次实验，证明了该方法的有效性，从而为谱仪研制和脉冲序列开发提供了一种辅助测试手段.

1 实验原理与设计

基于LabVIEW软件开发平台和高速数据采集卡DAQ-2005，本文设计了LabVIEW采集系统的前面板和程序框图，完成谱仪3个方向梯度，包括选层梯度G_z、频率编码梯度G_x和相位编码梯度G_y，以及1个数字输出波形信号的采集.设计控制模块完成对采集卡的启动、采样率、停止等参数功能的配置；设计采样模块实现波形的触发采样、连续采样、实时显示、同步存储等功能.波形数据文件进行格式转换后，进行数据分析，包括直方图统计、滤波、差分计算等分析过程.整体结构框图如图 1所示.

1.1 硬件部分

采集系统运用的是凌华科技（ADLINK）DAQe-2000系列的数据采集卡DAQ-2005.它是一款多通道、多功能的高性能数据采集卡，采用32位外部设备互连总线接口（PCI），通过高速度和高精度的模数转换（A/D）芯片，可实现高速、高精度同步数据采集.其可配置4个差分的模拟量输入（AI）通道达到最大的降噪，4个通道最高采样率可同时达到500 KHz；配置16位8 K先入先出（FIFO）存储器；配备模拟及数字触发功能、24通道输入输出端口（I/O）及2通道16位定时器/计数器.

1.2 软件部分

采集系统采用了NI公司的LabVIEW V2014作为软件开发平台.它具有直观多样的图形控件、友好的人机界面、强大的数据显示与处理和仪器控制能力，以及丰富的数据库函数等优点.

如图 2所示，利用LabVIEW建立虚拟仪器面板，加载数据采集卡ADLINK DAQPilot ExpressVI的驱动程序，对硬件进行自检和调试，对外接信号进行A/D转换，存入至DATA FIFO中，量化后的16位信号通过PCI总线存入计算机内存中.采集参数，如采样通道选择和采样率调整，是通过DAQPilot Ⅵ配置实现.实验中采用连续采样和触发采样两种模式.连续采样模式由While循环控制，采集启动与停止动作通过点击按钮或自定义固定时间参数完成.触发采样模式的触发信号由外接谱仪的输出提供.波形在采样过程中都外接谱仪的时钟信号.

采集程序将波形数据及其文件属性一起进行同步存储. 1次完整的成像过程持续数分钟，完成1次采样将产生10亿字节（GB）量级的数据文件.因此同步存储采用数据管理系统格式（Technical Data Management System，TDMS）以满足波形存储空间及快速写入与读取速度的需求. TDMS文件是以二进制方式存储数据；采用3层逻辑结构，分别为文件、通道组和通道；存储速度可达600 MB/s，能满足采集系统的需要.

图 2是梯度波形连续采样程序框图. 图 2中While循环控制整个采样过程的运行与停止；通过two.tsk的DAQPilot Ⅵ程序实现波形采样与输出；由波形图表进行实时显示，同时将4个通道的波形数据同步存储至名为datafile的TDMS文件中.

完成波形采样后，对TDMS的波形数据文件进行格式转换，保留通道的数据即波形数据，去除文件及通道组的数据即相关属性信息.

1.3 特征点提取

设定各个波形脉冲发生变化的点为特征点.通过以下的数据处理过程提取出特征点，并与预先设定的特征点进行比对，确定梯度波形发生问题的位置与时刻，最终判断硬件电路与脉冲序列是否存在问题.波形特征点提取过程如图 3所示.

第一步骤：对波形幅值做直方图统计分析.

直方图是概括梯度波形幅值属性的图形化方法，它以脉冲幅值为底边，以该幅值频数为高度.在直方图中，选定频数最多的直方图块（低电平值）进行均值计算，计算得到的幅值作为波形调零的参数.调零操作使得波形在零附近震荡，并且不改变其波形变化状态.选定最小脉冲计算3 db值，最小脉冲是直方图中最靠近低电平值的脉冲. 3 db值是参照滤波器通带截止与阻带截止概念所设定的值. 3 db值=0.707×E1，E1为最小脉冲的幅值. 3 db值作为阈值来判断脉冲的状态，在本文中，脉冲上升至3 db值大小时属于进入高电平状态；脉冲下降至3 db值大小时属于高电平结束状态.

第二步骤：滤波计算.

波形在采集过程中或多或少混叠了来自系统内部或外界的噪声干扰.因此需要进行噪声去除以保证后续波形特征点提取的正确性及波形平滑.中值滤波对脉冲噪声有良好的滤除作用，在滤除噪声的同时，能够保护信号的边缘使之不被模糊.选用中值滤波让脉冲序列的每一个点用该点邻域中各点的中位数代替，从而去除脉冲孤立的噪声点.中值滤波的表达式如(1)式所示：

$ {{s}_{1}}(i)=\text{median}[{{s}_{0}}(i-M), \ \cdot \cdot \cdot, \ {{s}_{0}}(i), \ \cdot \cdot \cdot, \ {{s}_{0}}(i+M)] $

(1)

其中$ {{s}_{0}}(i) $表示采集的原始梯度波形信号，$ {{s}_{1}}(i) $表示中值滤波信号，（2M+1）表示滤波窗的窗口长度.中值滤波后，进一步进行线性滤波以减少高斯噪声的影响，如(2)式所示：

$ {{s}_{2}}(i)=\sum\limits_{j=-N}^{N}{{{w}_{j}}{{s}_{1}}(i+j)} $

(2)

其中$ {{s}_{2}}(i) $表示线性滤波信号.为了避免边缘模糊，选择矩形窗. N表示矩形窗的长度，设置为3，$ {{w}_{j}} $表示矩形窗的系数.

对滤波数据取绝对值，即$ {{s}_{2}}=|{{s}_{2}}(i)| $.

第三步骤：一阶差分计算.

在物理上一阶差分表达物体运动的平均速度，在信号上表示幅值在时间上的变化率.对滤波数据进行一阶差分计算，获取相邻时间点脉冲的变化情况.选择两个方向的差分：前向差分和逆向差分，分别如(3)式和(4)式所示.

$ {{s}_{3}}(i)={{s}_{2}}(i+\delta )-{{s}_{2}}(i) $

(3)

$ {{s}_{4}}(i)={{s}_{2}}(i)-{{s}_{2}}(i-\delta ) $

(4)

其中δ表示差分间隔.在谱仪输出中，由于数字输出信号的突变性，δ一般取值较小，而3个梯度波形信号的变化较为平稳，δ取值相对较大.

第四步骤：特征点提取.

二阶差分表现了物体的变化趋势.当某个点的二阶差分达到最值时说明它是一个极值点，表明一阶差分值的符号正相反.若一阶差分值的符号左正右负表明此点是极大值点，反之为极小值点，体现在波形上表明抖动变化最为明显.二阶差分表达式如(5)式所示.

$ \nabla (i)=\left[\sum\limits_{i=-N}^{N}{{{s}_{3}}}(i) \right]/2N-\left[\sum\limits_{i=-N}^{N}{{{s}_{4}}}(i) \right]/2N $

(5)

N表示邻域大小.选择该点一阶差分值N邻域内的均值作为该点的差分值进行计算，使得二阶差分结果变化更为剧烈. N在数字输出信号时取值较小，可为1；而在梯度波形信号时取值相对较大，可为5.

二阶差分值的最值查找是根据阈值划分搜索范围.由阈值确定查找的起始时间，再分别向前和延后1 ms的时间范围内查找最大值和最小值.通过此方法找到的最值的位置就是脉冲抖动最为剧烈的点，即定义的脉冲特征点.

图 4是特征点提取过程的示意图.

其中图 4(a)、4(b)和4(c)分别是原始数据、调零后的数据和滤波数据，它们实际是不一样的，由于数据量与动态范围的问题，在matlab显示结果较为相似. 图 4(d)为差分计算过程，包括前向差分、逆向差分和二阶差分. 图 4(e)先由阈值判断出脉冲是上升还是下降状态，并确定查找的起始时间；然后向前和向后方向查找，如箭头方向所示.若脉冲处于上升状态，二阶差分值的最小值对应脉冲变化的起始位置，最大值对应脉冲变化的结束位置；若脉冲处于下降状态，二阶差分值的最大值对应脉冲变化的起始位置，最小值对应脉冲变化的结束位置；提取出的脉冲特征点用黑色“×”符号表示.

2 结果与讨论

选层梯度的波形特征点提取如图 5所示. 图 5(a)显示了整个波形的直方图. 图 5(b)选择8个连续脉冲为例，提取其特征点.这里，相位编码步为256，平均扫描次数（NS）为2，脉冲间隔（TR）为500 ms，采样率为500 kHz.其中直方图为了显示效果，对Y轴（数量）求取了对数. 表 1给出了图 5所提取出的选层梯度特征点的时间与幅值信息.

图 6给出了相位编码梯度的特征点提取结果示例.同样地，图 6(a)显示整个波形的直方图. 图 6(b)以8个连续脉冲为例，提取其特征点.其中直方图为了显示效果对Y轴（数量）求取了对数. 表 2给出了图 6所提取出的相位编码梯度特征点的时间与幅值信息.

当梯度受到了较为严重的外界噪声干扰时，运用同样的方法对读梯度进行处理. 图 7给出了读梯度（频率编码梯度）的特征点提取结果示例，图 7(a)显示了整个波形的直方图，图 7(b)以8个连续脉冲为例，提取其特征点.其中直方图为了显示效果对Y轴（数量）求取了对数. 表 3给出了图 7所提取出的读梯度特征点的时间与幅值信息.

图 7的直方图分布中除3个主要脉冲分量外，还包括其他分量，而且其他分量的数量较多，不可忽视的.但是，根据序列的发生特点，读梯度应该只具有3个脉冲分量，说明读梯度在采集的过程中受到较为严重的干扰.

我们设定脉冲幅值的预期误差为幅值的2%，脉冲时间差的预期误差为10 μs.从表 1、表 2、表 3可以看出，选层、读以及相位编码梯度的幅值误差为±0.02 V，脉冲的时间误差不超过2 μs，表明了所提取出的特征点的时间与幅值均在设计预期的误差范围内，这证明了该方法的有效性，可用该方法评估谱仪输出梯度和序列运行的准确度.

3 结论

基于LabVIEW软件和DAQ-2005数据采集卡，采集了MRI谱仪输出的梯度波形，对采集数据进行直方图统计、滤波、差分计算等分析，提取出波形的特征点，能够快速地判断谱仪梯度发生是否存在问题，以及确定发生问题的位置.这种方法准确度高、速度快，为实验室自研谱仪设计与脉冲序列开发提供了一种非常有效的辅助方法.