引言

磁共振眼球成像，已经广泛的应用于临床诊断中，例如用于葡萄膜黑色素瘤、视网膜母细胞瘤等眼部肿瘤的检测^[1-4].无需介入即可呈现物质内部结构，同时不同于光学成像方法，无需清晰的光线通路，可观测肿瘤的大小、位置、及邻近组织解剖结构情况等优点，使得磁共振成像（Magnetic Resonance Imaging, MRI）成为眼球研究的有力工具^[5-7].

然而，相比于光学成像方法，磁共振眼球成像也有它自身的一些缺点，例如分辨率较低、信噪比不足、容易产生运动伪影等.在MRI中，分辨率、信噪比和成像时间是3个相互制约的因素，它们不可能同时达到最优状态.对于眼球这个特殊的组织，它有较为精细的结构，需要较高的分辨率才可观察清楚.要想得到高分辨率和高信噪比的磁共振眼球图像，就需要较长的扫描时间.而人的眼球是一个易于活动的器官，很容易产生运动，较长的扫描时间会不可避免的使磁共振眼球图像产生运动伪影^[8]，降低图像质量.在闭上眼睛时，眼球会产生不由自主的全方位的运动，既有平面内的运动，又有跨越平面的运动.

为了解决磁共振眼球成像运动伪影的问题，一些学者提出了不同的策略，最常见的是“周期性眨眼的序列”^{[2, 9, 10]}.然而，这种方法需要受试者非常精确的配合.受试者在接受扫描时需要将眼睛目光汇聚在一个固定的位置，根据系统自动发出的提示音，每间隔一个短暂的时间眨眼一次，然后继续将目光盯住那个固定的位置.其中眨眼的瞬间机器是暂停扫描的.这种成像方法不经过训练很难一次成功.当受试者配合有偏差时，这种方法就会使图像质量明显下降.另外，近期还有学者在研究磁共振眼球成像方法时，让受试者睁开眼睛盯住一个位置，用探针检测眨眼，在不成像的眼睛涂上含有氧化铁的睫毛膏，当系统检测到受试者眨眼时，重新采集被破坏的K空间数据^[11].这种方法抑制了眨眼伪影，但同样需要受试者的配合，同时还需要额外的硬件设备并将其集成到系统中，增加了成像的复杂性.

现阶段大多数磁共振运动伪影校正方法只对整体的刚性运动有效，对于类似眼球这种存在局部的全方位运动的特殊部位，有学者在解决运动伪影的问题时提出了一种“快速磁共振成像三维体重建”（Snapshot magnetic resonance imaging with Volume Reconstruction, SVR）的方法，并将其应用到孕妇体内的胎儿大脑MRI中^[12].这种方法使用单层高速成像序列扫描以“冻结”胎儿大脑的运动.经过连续的移位采集以覆盖整个成像区域并确保足够的采样密度.使用6自由度刚性配准进行体积内的单层图像的回顾性对准.然后，应用基于B样条散射数据内插（B-spline Scattered Data Interpolation, SDI）的数据融合算法以生成最终的三维胎儿大脑图像.这种成像方法可以有效的抑制胎儿大脑的运动伪影，但对于眼球成像，其采集分辨率远高于胎儿大脑，以至于成像序列的单层采集时间无法足够短，进而难以“冻结”眼球的运动，部分图像会被运动破坏导致后续的配准不准确.尽管有磁共振并行成像方法可以加快扫描速度，但并行成像重建引起的信噪比损失会使一些图像信号淹没在噪声之中，同样使配准有时难以准确实现.这种方法并不适合眼球这一部位.

针对磁共振眼球成像中存在的这些问题，本文研究并实现了一种以超分辨率理论为基础的新的磁共振眼球成像方法.超分辨率，是一种利用一系列低分辨率图像生成一幅高分辨率图像的技术.在磁共振成像中，超分辨率更多的被应用于层厚方向分辨率的提升^{[13, 14]}.有的学者研究了平面内的磁共振超分辨率，对于静止的物体，通过平面内成像视野区域（Field Of View, FOV）的亚像素位移，实现平面内的图像分辨率的提升^[15].这一方法刚一提出便受到了争议^[16].由于平面内成像视野区域的亚像素位移对应于将采样后的K空间信号进行线性相位调制，不同的具有亚像素位移的成像视野区域所获取的磁共振图像的信号除了噪声之外并没有不同之处，不能提供用于平面内超分辨率的额外的信息.然而，另有学者指出，如果物体本身在多次扫描之间发生了运动，新的额外的信息会被引入^[17-19]，使平面内的超分辨率成为可能，但是其额外引入的信息量的大小未能被精确确定^[19].

本文提出了一种新的磁共振平面内超分辨率的策略，实现了平面内的超分辨率重建，并将这一技术运用于磁共振眼球成像中.由于在磁共振成像过程中，平面内的频率编码方向和相位编码方向可以有不同的采集分辨率.在采集磁共振信号时，可以在同一成像视野区域内采集两组数据，它们的相位编码方向正交，每组数据的频率编码方向和相位编码方向具有不同的采集分辨率，使得这两组数据在不同的方向上分辨率互补.另外通过旋转成像视野区域使得各方向上的采集分辨率都得到互补，运用合适的超分辨率算法，可使平面内的分辨率得到提升.

为了验证本文算法的有效性，我们先做了仿体的测试实验，产生了较好的效果，之后又在人体上采集了眼球数据，根据本文提出的新的磁共振眼球成像方法，受试者在接受扫描时，闭上眼睛，不需要做额外的配合工作.对每一个受试者，采集6组数据，这6组数据共享完全相同的参数，只是成像视野区域相对于中心分别旋转了一定的角度.这6组数据可分为3对，其中每对的两组数据在相位编码方向正交，频率编码方向和相位编码方向上的采集分辨率互补.每组数据都含有若干动态图像.由于眼球的不自主运动是刚性的，对于这些动态的数据，可以对眼球的运动做平面内的配准，然后丢弃一部分配准系数较低的帧.配准是超分辨率的一个必要步骤，由于在磁共振眼球成像时，层内的眼球运动可以通过配准来描述.而当不同的动态图像之间发生层间的运动时，会在眼球的不同层面成像，导致配准系数降低.显然层间的运动会降低配准系数，因此丢弃一部分配准系数较低的图像可以适当的减少层间的运动伪影.此外，由于采集每一帧图像的动态时间不可能足够短，导致部分动态图像会在一定程度上被运动破坏，这些运动伪影较大的帧也会因为配准系数较低而被丢弃.之后，本文通过最大后验概率（Maximum a Posteriori Probability, MAP）超分辨率算法^[20-22]提升层内的分辨率并且抑制了运动伪影.根据本文的采集策略，多组数据的分辨率可以互补，生成一个更高分辨率的图像.另一方面，眼球的不自主运动也可以提供超分辨率所需的亚像素位移.实验结果表明，本文方法在不需要受试者做额外配合、也不需要额外的磁场探针等硬件设备的情况下，可以得到更清晰的磁共振眼球图像，抑制了运动伪影的同时，提高了图像的分辨率，使图像的质量得到了提升.

1 理论

基于超分辨率重建算法，是将一系列同一场景但有差异的低分辨率图像合成为一幅高分辨率图像.1984年，Tsai和Huang^[23]开创了由低分辨率图像序列复原单幅高分辨率图像的工作，给出了频域的复原方法.该方法计算简单，但模型中没有考虑光学系统点扩散函数（Point Spread Function, PSF）、运动模型和观测噪声等因素的影响，且只能处理图像整体的运动，对于眼球这种局部运动无能为力.空域法是当前的主流方法之一，将复杂的全局和局部运动、光学模糊、空间可变PSF，以及非理想采样并入其中，并且具有很强的包含空域先验约束的能力.

对于空域法，可以设定一个由高分辨率到低分辨率的图像观测模型，即假设采集到的图像是由高分辨率图像变换得到^[24].这些变换包括几何平移、旋转、图像PSF模糊、降采样等操作.现假设高分辨率图像为X，低分辨率图像为y_k，k，1，2，3，….对于每幅低分辨率图像y_k，k，1，2，3，…，可视为由高分辨率图像X，经过刚性几何变换G_k，PSF模糊B_k，降采样D_k，再加上噪声V_k而得到，即

$ {{y}_{k}}={{D}_{k}}{{B}_{k}}{{G}_{k}}X+{{V}_{k}}, k=1, 2, 3, \cdots $

(1)

令${{W}_{k}}={{D}_{k}}{{B}_{k}}{{G}_{k}}$，可得

$ {{y}_{k}}={{D}_{k}}{{B}_{k}}{{G}_{k}}X+{{V}_{k}}={{W}_{k}}X+{{V}_{k}}, k=1, 2, 3, \cdots $

(2)

根据这一模型，超分辨率重建就转换为求解X的数学问题，这是一个不适定问题，有许多求解方法.本文采用了MAP法^[20-22]求解.MAP法是空域法中统计方法的典型代表，它把超分辨率问题看成一个统计估计问题，在已知一系列低分辨率图像序列的前提下，使出现高分辨率图像的后验概率达到最大，再利用贝叶斯公式，得

$ {{X}_{MAP}}=\underset{X}{\mathop{\arg \max }}\, P(X|{{y}_{1}}, {{y}_{2}}, {{y}_{3}}, \cdots )=\underset{X}{\mathop{\arg \max }}\, P(X|Y)=\underset{X}{\mathop{\arg \max }}\, \frac{P(Y|X)P(X)}{P(Y)} $

(3)

将(3)式中目标函数取对数，不改变其极值点，得：

$ {{X}_{MAP}}=\underset{X}{\mathop{\arg \max }}\, [\ln P(Y|X)+\ln P(X)-\ln P(Y)]=\underset{X}{\mathop{\arg \max }}\, [\ln P(Y|X)+\ln P(X)] $

(4)

假设图像噪声具有多变量高斯分布，所有低分辨率图像之间的噪声变量是相互独立的，于是

$ P({{V}_{k}})=P({{y}_{k}}-{{W}_{k}}X) $

(5)

根据Gauss随机场模型，在给定原始图像条件下，由于随机噪声的影响，采集观测到的图像数据y_k可视为一个随机场.可以认为该随机场的均值是W_kX.按照标准多变量高斯模型，给定X条件下y_k的分布应该符合下式：

$ P({{y}_{k}}|X)=P({{y}_{k}}-{{W}_{k}}X)=P({{V}_{k}}) $

(6)

因此

$ P({{y}_{k}}|X)=\frac{1}{\sqrt{{{(2\pi )}^{L}}\det |{{Q}_{\eta }}|}}\exp \left[-\frac{1}{2}{{({{y}_{k}}-{{W}_{k}}X)}^{T}}Q_{\eta }^{-1}({{y}_{k}}-{{W}_{k}}X) \right] $

(7)

其中L是图像y_k的大小（像素点个数）；${{Q}_{\eta }}=E[\eta {{\eta }^{T}}]={{\sigma }^{2}}I$，表示噪声的协方差矩阵，其中η表示低分辨率图像的噪声向量，而I则是单位矩阵.观测图像的这个条件分布模型是众多图像复原算法的基础.

(4) 式中P(X)可视为高分辨率图像X的先验约束项，不同类型的图像先验知识可以被用作估计P(X).Schultz和Stevenson基于统计估计框架，提出的一种具有边缘保持能力的Huber-Markov模型在MAP超分辨率算法中得到了广泛的应用^[21]，其数学表达式为：

$ P(X)=\frac{1}{Z}\exp \left\{ -\frac{1}{2\beta }\sum\limits_{c\in C}{\rho }[d_{c}^{t}(X), \alpha] \right\} $

(8)

在(8)式中，Z为归一化常数；β为Gibbs先验温度参数；c是一个局部像素组，称为团块，而图像上所有的团块构成集合C；$d_{c}^{t}(X)$是空域活动性测度，在高分辨率的平滑区域$d_{c}^{t}(X)$的值较小，而在边缘区域，$d_{c}^{t}(X)$的值较大.这里空域活动性测度即是计算高分辨率图像中每一个像素的二阶微分，即二阶邻域系统：

$ d_{m, n, 1}^{t}(X)={{X}_{m, n-1}}-2{{X}_{m, n}}+{{X}_{m, n+1}} $

(9)

$ d_{m, n, 2}^{t}(X)=0.5{{X}_{m+1, n-1}}-{{X}_{m, n}}+0.5{{X}_{m-1, n+1}} $

(10)

$ d_{m, n, 3}^{t}(X)={{X}_{m-1, n}}-2{{X}_{m, n}}+{{X}_{m+1, n}} $

(11)

$ d_{m, n, 4}^{t}(X)=0.5{{X}_{m-1, n-1}}-{{X}_{m, n}}+0.5{{X}_{m+1, n+1}} $

(12)

其中，${X_{m, n}}$表示图像中(m, n)位置的像素值，$d_{m, n, 1}^t(X)$、$d_{m, n, 2}^t(X)$、$d_{m, n, 3}^t(X)$、$d_{m, n, 4}^t(X)$分别对应于图像水平、对角、垂直和反对角方向的变化程度.

图像数据边缘像素点的概率被Huber边缘惩罚函数^{[21, 25]}所控制，其定义为：

$ \rho (x, \alpha ) = \left\{ \begin{array}{l} {x^2}, \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ |x|\alpha \\ 2\alpha |x| - {\alpha ^2}, \ \ |x| > \alpha \end{array} \right. $

(13)

(13)式中，α为阈值参数，控制先验模型的不连续性大小.将(7)~(13)式、代入(4)式中，得到高分辨率图像的最大后验概率估计为：

$ {X_{MAP}} = \mathop {\arg \min }\limits_X \left\{ {\sum\limits_{k = 1}^N {||{y_k} - {W_K}X||_2^2} + \gamma \sum\limits_{(m, n)} {\sum\limits_{r = 1}^4 \rho } [d_{m, n, r}^t(X), \alpha]} \right\} $

(14)

其中γ是P(X)先验知识项的权重系数.

2 材料与方法

实验包括仿体实验及在体实验两部分，仿体成像实验对象为圆柱水模，在体成像实验部位为志愿者眼部区域，感兴趣区域（Region Of Interest, ROI）为眼球.所有的扫描实验均在Philips 3.0 T Achieva TX MRI设备（Philips Healthcare, Best, The Netherlands）上完成.

2.1 仿体实验

对圆柱水模，我们在横断面先后采集了6组数据，其中前两组的成像视野区域绕着FH轴相对区域中心旋转了0°，中间两组旋转-30°，最后两组旋转30°.在每对数据中，相位编码方向相互正交，相位编码方向的分辨率设置为1 mm，而频率编码方向的分辨率为4 mm，层厚为4 mm.采用2D T1FFE序列，采集一层图像，成像视野区域为256×256 mm²，翻转角（Flip Angle, FA）为5°，脉冲间隔重复时间/回波时间（Repetition Time/Echo Time, TR/TE）= 2.5 ms/2.1 ms，半傅里叶采集K空间52.5%的数据.

对于PSF，本文选用高斯函数，使其半高宽等于图像像素点的单位大小.P(X)先验知识项的权重系数γ设置为0.002.在求解X_MAP时，使用梯度下降法迭代求解.

$ {X_{n + 1}} = {X_n} + \left\{ {\sum\limits_{k = 1}^N {W_k^T({y_k} - {W_K}{X_n}) + } \gamma \sum\limits_{(m, n)} {\sum\limits_{r = 1}^4 {{{(d_{m, n, r}^t)}^T}\rho '[d_{m, n, r}^t({X_n})} }, \alpha]} \right\} $

(15)

其中$\rho '(x, \alpha )$是$\rho (x, \alpha )$的导函数，定义为

$ \rho '(x, \alpha ) = \left\{ \begin{array}{l} 2x, \ \ \ \ \ \ \ \ \ \ \ \ \ \ |x|\geqslant\alpha \\ 2{\rm{sign}}(x)\alpha, \ \ |x| > \alpha \end{array} \right. $

(16)

需要注意的是，所有的图像在计算时都需要按照先列后行的顺序拉成一个向量，因而变换矩阵W_k是一个大型稀疏矩阵，它作用于图像向量时，相当于对图像做某种线性变换，最终求得的X_MAP也是一个向量，还需重新变形成对应大小的图像矩阵.

2.2 在体实验 2.2.1 准备工作

为了得到更高质量的磁共振眼球图像，我们定制了一个特殊的线圈，如图 1所示.该线圈为4通道眼球线圈，每个眼球各有两个通道.在扫描前，我们先将受试者成像的眼部盖上一层湿的纱布以减小与空气交界处的场不均匀性伪影^{[2, 26, 27]}.带上眼球线圈之后，我们用绑带将线圈绑在面部上方，增加稳定性，使其更加贴近眼部，增加图像的信噪比，如图 2所示.

2.2.2 数据采集

对于眼球，和水模数据的采集一样，我们在横断面先后采集了6组数据，同样对成像视野区域绕着FH轴方向相对区域中心进行了旋转操作.在每对数据中，相位编码方向相互正交.为了保持每一帧动态的图像具有一定的信噪比，相位编码方向的分辨率设置为0.2 mm，而频率编码方向的分辨率为0.6 mm，层厚为3 mm.所有的数据采集都使用2D T1FFE序列，仅采集一层，成像视野区域为192×192 mm²，翻转角为7°，TR/TE = 8.5 ms/4.1 ms，半傅里叶采集K空间52.5%的数据以加快动态采集时间.每组数据含有20帧，每一帧的动态时间为4.3 s，总共的采集时间为8 min 36 s.频率编码方向上较低的分辨率可以增加动态图像的信噪比，而动态扫描可以很大程度上抑制运动伪影.每一帧的动态时间需要尽可能的短，以“冻结”眼球的运动，部分未被“冻结”的眼球图像会在随后的处理中被丢弃.受试者在接受扫描时只需将眼睛闭上，不需要做任何睁眼眨眼等配合的工作.由于扫描是动态的，我们可以观察到眼球的不自主运动.

2.2.3 预处理工作

由于在数据采集时，平面内的分辨率较高，采集时间较短，导致图像有较低的信噪比.此外，由于眼球不自主运动的存在（包括平面内的运动和穿越平面的运动），部分图像帧被严重的破坏，这些帧的信号强度明显异常.在数据预处理部分，我们将剔除掉这部分信号强度异常的帧.剔除的比例为10%，对于每组20帧的动态扫描，保留其中的18帧数据.

在剔除信号异常的数据时，对于每组数据，我们先计算出每一帧的图像的信号强度，将所有帧的信号强度按顺序保存为一个向量$S = {[{s_1}, {s_2}, {s_3}, \cdots]^T}$.之后，根据每一帧的信号强度与总体信号强度平均值的偏离程度，去除掉偏离程度最大的那一帧数据.剔除一帧数据后，重新计算所有帧的平均信号强度，再剔除掉信号强度与总体信号强度平均值的偏离程度最大的那一帧数据.如此重复下去，直到剔除完预先设定好的图像帧的数量.

然后，我们对噪声的方差进行估计.假设图像噪声部分的实部和虚部符合独立的、有着相同的方差的高斯分布，那么其幅值符合瑞利分布，设方差为σ².根据数学理论可知，由两个独立的、均值相同的、方差为σ²的瑞利分布的差产生的随机变量，其均值为0、方差为2σ².对于每一组数据，若干动态图之间的噪声的幅值是相互独立的.在成像条件相同的条件下，其瑞利分布的均值和方差均大体一致.我们利用帧差法的思想，将每组数据若干帧动态的图像两两相减，试图通过相减去除掉图像本身的信号，而只留下背景噪声部分.两两相减后共得到若干幅帧差图，计算所有的帧差图的方差，找到其中最小的那一个，其对应的方差即近似为原始图像的噪声方差的2倍.因为方差最小的帧差图，已经通过相减几乎完全消除了图像信号本身的作用，而只保留了噪声部分.而方差较大的帧差图，则由于眼球运动的原因，遗留下了部分图像本身信号的作用，导致帧差图方差较大.

根据之前估计的图像噪声的方差σ²，我们使用BM3D算法^[28]对每一幅动态图进行不同噪声方差估计下的去噪（协同滤波时选择合适的阈值）.由于去噪会在一定程度上影响图像边缘部分的高频信息，我们将保留所有的方差估计下的去噪结果，在之后选择适当的去噪图像进行后续的计算.

2.2.4 图像配准

如果眼球的运动只发生在平面内，那么由于眼球的球体部分的运动是刚性的，所以可以用刚性的配准描述眼球在平面内的运动情况.然而，穿越平面的运动也时常存在，导致图像前后采集的不一致性.另一方面，由于图像分辨率较高，导致动态采集的周期不能足够短，在采集某一帧数据时，眼球的运动仍可以破坏在这一帧数据，使其产生运动伪影.鉴于上述情况，我们将根据配准的结果丢弃一部分数据，以减少层内和穿越层厚方向的运动伪影.

图像配准是超分辨率算法中关键的一步，其配准的精度对最终结果有至关重要的影响，必须要达到亚像素级别.在这一步中，对于每一组数据，我们都需要先选取一个参考帧.参考帧的眼球应是几乎没有运动伪影且恰好在所成像的平面内，几乎没有发生穿越平面的运动.由于眼球的不自主运动像是一种围绕某个中心的运动，对所有的动态图像做平均所得到的眼球图像可以作为匹配模板的近似估计.我们以图像矩阵的相关系数为标准，选取与平均图像相似度最大的帧作为配准的参考图像，将其他所有的图像与参考图像做配准.对于二维平面的刚体配准，有两个方向的平移量和旋转角度共3个自由度.

为了提高配准的精度和速度，我们选择了两种配准相似度度量标准，分别是相关系数和绝对差.这两种标准各有优劣，我们试图参考两者综合的结果作为配准的最终结果，以弥补单一标准下的不稳定性.此外，我们将图像进行双线性插值放大，以提高亚像素级别的配准精度.另一方面，我们将智能计算中的粒子群算法^{[29, 30]}应用到配准过程中，以加快其计算速度及准确性.由于配准问题的规模和复杂度并不非常大，考虑到计算时间的问题，种群粒子数和迭代次数都不必过大，本文把两者都设置为35.对于3个自由度的刚性配准（旋转角度和x、y方向的位移），粒子的长度设置为3.考虑到眼球的运动在一定的范围内，粒子的范围设定也需要恰当，本文设置为$\left[{-{{10}^ \circ }, -\frac{{{N_2}}}{{15}}, -\frac{{{N_1}}}{{15}}} \right] \sim \left[{{{10}^ \circ }, \frac{{{N_2}}}{{15}}, \frac{{{N_1}}}{{15}}} \right]$，N₁和N₂是图像两个维度的大小.粒子运动的最大速度为$\left[{{2^ \circ }, \frac{{{N_2}}}{{75}}, \frac{{{N_1}}}{{75}}} \right]$.速度压缩因子为1，学习因子为${c_1} = {c_2} = 1.494$.惯性权重则由初始值0.9在迭代过程中线性递减为0.000 2，使迭代最后阶段算法的精度得以提升.

考虑到对眼球做二维的刚体变换时，ROI会发生改变，用模板图像去配准待配准图像时，待配准图像的ROI会对相似性度量标准产生影响.因此，我们需要在正式配准前先固定待配准图像的ROI，然后再做亚像素配准.具体的配准步骤如下.

1.选取好参考图像，选用噪声方差估计为0.8σ²的去噪结果，将眼球球体部分非常仔细地手动画入ROI，并计算ROI的几何中心的坐标.其中画ROI时需要将图像放大32倍以上，以保证其准确性.

2.对每一帧待配准图像，选用噪声方差估计为0.8σ²的去噪结果，利用粒子群算法寻找使相关系数达到最大的配准参数以及对应的待配准图像的ROI.此时的配准参数结果并不十分准确，只是用于计算待配准图像的ROI.

3.对每一帧待配准图像，选用噪声方差估计为0.8σ²的去噪结果，固定之前计算得到的待配准图像的ROI，利用粒子群算法寻找使相关系数达到最大的配准参数.

4.对每一帧待配准图像，选用噪声方差估计为0的去噪结果（即没有去噪），固定之前计算得到的待配准图像的ROI，利用粒子群算法寻找使绝对差达到最小的配准参数.

5.将第3步和第4步计算得到的配准参数取平均值.

之前提到，由于运动伪影的存在，并不是所有的眼球图像都有较好的质量.在配准工作结束之后，我们需要继续丢弃一部分数据.对于每一组数据，我们将每一帧眼球图像的配准相关系数按降序排序，将配准绝对差按升序排序.根据两个排序结果，可以将每一帧图像确定一个平均排名，显然排在第一位的应当是参考模板图（参考图的配准相关系数为1，绝对差为0）.最后我们丢弃那些平均排名较低的图，可以在一定程度上抑制运动伪影.这里，我们将之前保留的每组的18帧数据剔除其中排名最靠后的5帧，最终每组数据留下13帧图像.

2.2.5 超分辨率

对于平面内的PSF，我们选用高斯函数，使PSF的半高宽等于高分辨率像素点的单位大小.根据之前计算得到的配准结果和对PSF的估计，我们已经获知了每一幅图像的几何变形参数G_k和PSF模糊B_k，D_k也是已知的.P(X)先验知识项的权重系数γ=0.01.对于X_MAP的求解，使用与仿体实验相同的方法.

另外需要说明的是，对每组数据分别做完配准和超分辨率之后，得到了6幅不同方向上分辨率互补的图像，这6幅图像还需要继续用同样的方法做配准和超分辨率，生成一幅最终的磁共振眼球图像.本文方法的流程图如图 3所示.

3 结果与讨论 3.1 仿体实验

在人体眼球实验之前，我们在仿体上作了超分辨率的测试.依照之前所叙述的采集策略，我们获得了若干幅低分辨率图像，并直接采集了高分辨率的图像作为参考图，如图 4所示.

仿体实验结果表明，我们提出的采集策略配合MAP超分辨率算法确实可以提高图像分辨率，成功的实现了磁共振平面内的超分辨率.

3.2 在体实验 3.2.1 噪声能量的估计及去噪的结果

我们用帧差法将每一组若干动态图像两两相减后，得到了若干帧差图像，然后对这些图像的方差由小到大升序排序.图 5的(a)、(b)、(c)、(d)分别显示了某组数据中排序为1、60、120、180的帧差图，其方差依次为225.12、249.14、282.37、336.96.可以看出，方差小的帧差图(a)，其图像信号被最大程度的消除掉，背景噪声得到了保留.由于眼球不自主运动的存在，部分图像相减后信号未被消除，导致方差较大.我们计算图 5(a)的方差就可以近似估计出图像背景噪声的能量.

在得到了噪声的能量（方差）大小之后，我们使用BM3D算法^[28]对图像进行去噪.某一帧图像去噪的结果如图 6所示.这里去噪只是为配准做准备，提高配准的相关系数及准确性，而不是用去噪后的图像做超分辨率.

3.2.2 配准及异常数据的剔除

我们任意地选取了一组数据，将两种不同的相似度度量下的配准参数绘制成曲线图（图 7）.可以看出，相关系数和绝对差用作相似度度量时，其配准结果比较吻合，但依然有误差.我们将两组配准参数取平均值，以降低单一标准下的误差.

根据本文的方法，对每一组数据，最先根据信号强度丢弃20幅动态图像中的2幅，然后根据配准结果，再丢弃配准系数排名较靠后的5幅，最后剩下13幅图像用于超分辨率.我们任意选取一组数据，将20帧图像根据配准排序结果展示如下（图 8），其中最后2幅图像是根据信号强度丢弃的.

3.2.3 超分辨率结果

由图 9的结果可以看出，本文的算法可以使采集到的磁共振眼球图像质量得到提升，抑制了运动伪影的同时，提高了图像分辨率.特别是眼球边缘处，本文的方法可以更清晰的观察到眼球最外层的膜结构.

3.3 讨论

本文从眼球成像的实际问题出发，为解决磁共振眼球成像中信噪比、分辨率和运动伪影三者之间的矛盾，提出了一种新的眼球图像获取方法.首先在成像时将受试者眼部盖上湿的纱布以抑制场不均匀性伪影.使用一种特制的眼球线圈，采集多组分辨率在不同方向上互补的动态眼球图像，再通过预处理，配准，异常数据丢弃，超分辨率等步骤，使磁共振眼球图像质量得到了提升.在整个图像数据采集过程中，受试者不需要做任何睁眼眨眼等配合工作，也不需要额外的磁场探针等硬件辅助，降低了成像的难度，减少了受试者的不舒适性.

需要说明的是，本文的方法只是为磁共振眼球成像提供了一种新的思路，其中的成像参数与数据处理中的参数并不是绝对的.在实际应用时，可以通过增加采集时间换取更高分辨率更高信噪比的眼球图像.异常数据丢弃的比例也是可以调整的，当丢弃数据过少时，会导致运动伪影的增加，丢弃数据过多时，会导致信噪比不足.此外，超分辨率算法中的先验知识权重系数也应选取为恰当适中的值.

最后要指出本文的算法所展现的效果来自3个方面共同作用的结果：一是运动伪影的消除；二是眼球不自主运动为磁共振平面内超分辨率提供的额外信息；三是不同方向互补的分辨率所带来的效果.这三个方面都为磁共振眼球图像质量的提升做出了贡献.

4 结论

本文提出了一种新的眼球图像获取方法，在受试者不需要做任何配合工作的情况下，使磁共振眼球图像的运动伪影减少，分辨率提高，图像质量得到了提升.这一新的思路为临床上眼球疾病的诊断提供了新的方法，具有一定的实用价值.