CN100578546C - 基于复二维奇异谱分析的磁共振部分k数据图像重建方法 - Google Patents

Abstract

本发明涉及一种基于复二维奇异谱分析的磁共振部分K数据图像重建方法，包括从磁共振成像扫描仪中预设的图像空间范围中采集部分K数据、对该部分K数据进行补零成像处理、根据补零成像得到的近似图像及部分K数据进行模型参数估计、利用复系数加权二维奇异函数的磁共振图像的数学模型和复二维奇异谱分析模型进行磁共振复图像的重构。采用了该发明的基于复二维奇异谱分析的磁共振部分K数据图像重建方法，节省了扫描时间，图像具有高信噪比、高分辨率和高精确度，为医学核磁共振成像检测提供了高质量的可靠图像信息，适用范围较为广泛，给人们带来了很大便利，为医学成像检测技术的进一步发展和大范围普及应用奠定了坚实的理论和实践基础。

Description

基于复二维奇异谱分析的磁共振部分K数据图像重建方法

技术领域

本发明涉及医学成像检测技术领域，特别涉及快速磁共振成像技术领域，具体是指一种基于复二维奇异谱分析的磁共振部分K数据图像重建方法。

背景技术

随着现代医学技术的不断发展，核磁共振成像(MRI)技术已经成为医学成像检测领域中不可或缺的手段。磁共振信号空间(原始数据空间)称为K空间，即为傅里叶变换空间，K空间采样信号经过傅里叶反变换后再取模，即得到核磁共振(MR)图像。实现快速成像的方法主要有提高硬件设备档次(比如提高主磁场强)，采用快速策略(如FLASH，EPI等等)、部分K空间扫描(如半谱扫描，非中心对称扫描等等)以及非直角网格扫描(如螺旋扫描)等方法。其中一维部分K空间数据成像可以在硬件和扫描方式不变的前提下，成倍地提高扫描速度(请参阅文献：P.Margosian，F.Schmitt，D.Purdy，“Faster MR Imaging：Imaging with Halfthe Data，”Health Care Instrum.，vol.1，pp.195-197，1986.，和J.van Cuppen and A.van Est，“Reducing MR imaging time by one-sided reconstruction，”Mag.Reso.Imag.，vol.5，pp.526-527，1987.)。目前流行的策略是基于相位校正的部分K空间数据图像重构，典型的方法有半谱POCS、相位校正共轭对称法和HM法等等(请参阅文献：E.M.Haacke，E.D.Lindskog，W.Lin，″A fast，iterative partial Fourier technique capable of local phase recovery，″J.Magn.Reson.，vol.92，pp.126-145，1991.，V.A.Stenger，D.C.Noll，F.E.Boada，“Partial k-space reconstruction for 3Dgradient echo functional MRI：A comparison of phase correction methods，”Magn.Reson.Med.，vol.40，pp.481-490，1998.，D.C.Noll，D.G.Nishimura，A.Macovski，“Homodyne detection inmagnetic resonance imaging，”IEEE Trans.Med.Imag.，vol.10，PP.154-163，1991.，和G.McGibney，M.R.Smith，S.T.Nicholas and A.Crawley，”Quantitative evaluation of several partial-Fourierreconstruction algorithms used in MRI，”Magn.Reson.Med.，vol.30，pp.51-59，1993.)。

与此同时，二维部分K空间数据扫描可以比一维部分K空间数据扫描节省一半扫描时间，但是其成像方法只有补零法(FZI法)，即未采集的K空间数据补零后，再用离散傅里叶反变换成像。由于补零法成像总是有伪影，这是由于该方法是利用牺牲图像质量来换取成像速度的。这些缺点所造成的失真现象足以使临床诊断医生产生误诊，以致始终阻碍着它们进入医学临床应用的大门，这样就给人们的工作和生活带来了很大的不便，并在一定程度上限制了医学成像检测技术的进一步发展。

因此，提高部分K空间数据成像速度的关键问题是要寻找一种新的图像表示方法，即一种新的图像模型，在这种模型下，可以用较少的变量来表示图像。本发明的发明人曾经在文献中(具体请参阅文献：Luo JH，Zhu YM，MR image reconstruction from truncated k-space usinga layer singular point extraction technique，IEEE TRANSACTIONS ON NUCLEAR SCIENCE 51(1)：157-169 Part 1 FEB 2004)指出，任何一个实数字信号都可以用奇异函数的加权和表示。因此，这样的模型适合实信号的一维部分K空间数据以及在二维K空间中只有一个方向的K空间数据被截断条件下的图像重构问题。

发明内容

本发明的目的是克服了上述现有技术中的缺点，提供一种能够提高重建磁共振复图像的速度、有效降低图像误差、提高磁共振图像质量、高效实用、工作性能稳定可靠、适用范围较为广泛的基于复二维奇异谱分析的磁共振部分K数据图像重建方法。

为了实现上述的目的，本发明的基于复二维奇异谱分析的磁共振部分K数据图像重建方法如下：

该基于复二维奇异谱分析的磁共振部分K数据图像重建方法，其主要特点是，所述的方法包括以下步骤：

(1)从磁共振成像扫描仪中预设的图像空间范围中采集部分K数据G(k_x，k_y)，所述的K数据G(k_x，k_y)即为离散傅里叶变换空间数据G(k_x，k_y)；

(2)对该部分K数据进行补零成像处理；

(3)根据补零成像得到的近似图像及部分K数据进行模型参数估计；

(4)根据模型参数估计的结果，利用复系数加权二维奇异函数的磁共振图像的数学模型和复二维奇异谱分析模型进行磁共振复图像的重构。

该基于复二维奇异谱分析的磁共振部分K数据图像重建方法的复系数加权二维奇异函数的磁共振图像数学模型为：

$g(x,y)=\underset{i=1}{\overset{Q}{\mathrm{\Σ}}}{a}_i{u}_{{\mathrm{\δ}}_x}(x-x_i,y-y_i)$ 0≤x，y＜N；

其中，g(x，y)，0≤x，y＜N为像素为N×N的二维磁共振图像的复图像信号，(x_i，y_i)为g(x，y)的第i个奇异点，a_i为该第i个奇异点上的复奇异值，Q为g(x，y)的奇异点的个数，

为以(x_i，y_i)为奇异点的二维奇异函数。

该基于复二维奇异谱分析的磁共振部分K数据图像重建方法的复二维奇异谱分析模型为：

$G(k_x,k_y)=\underset{i=1}{\overset{Q}{\mathrm{\Σ}}}{a}_i{U}_{{\mathrm{\δ}}_{x,i}}(k_x,k_y),$ 0≤k_x，k_y＜N；

其中，G(k_x，k_y)＝DFT[g(x，y)]， $U_{{\mathrm{\δ}}_{x,i}}(k_x,k_y)=\mathrm{DFT}\left[u_{{\mathrm{\δ}}_x}(x-x_i,y-y_i)\right],$ DFT[·]为二维离散付里叶变换算子。

该基于复二维奇异谱分析的磁共振部分K数据图像重建方法的进行补零成像处理包括以下步骤：

(11)将K空间Ω划分为有数据子空间Ω_k和无数据子空间Ω_k；

(12)将Ω_k子空间中的数据用零代替，Ω_k空间保持不变，并根据以下公式分别得到G(k_x，k_y)和

补零成像后的频谱数据

和

$\tilde{G}(k_x,k_y)=\left\{\begin{array}{cc}G(k_x,k_y),& k_x,k_y\∈{\mathrm{\Ω}}_k\\ 0,& k_x,k_y\∈{\stackrel{\‾}{\mathrm{\Ω}}}_k\end{array}\right.,$ ${\tilde{U}}_{{\mathrm{\δ}}_{x,i}}(k_x,k_y)=\left\{\begin{array}{cc}{U}_{{\mathrm{\δ}}_{x,i}}(k_x,k_y),& k_x,k_y\∈{\mathrm{\Ω}}_k\\ 0,& k_x,k_y\∈{\stackrel{\‾}{\mathrm{\Ω}}}_k\end{array}\right.;$

(13)根据以下公式计算G_Δ(k_x，k_y)和

$U_{{\mathrm{\Δ\δ}}_{x,i}}(k_x,k_y)=\mathrm{DFT}\left[{\mathrm{\Δu}}_{{\mathrm{\δ}}_x}(x-x_i,y-y_i)\right];$

$G_{\mathrm{\Δ}}(k_x,k_y)=\underset{i=1}{\overset{Q}{\mathrm{\Σ}}}{a}_i{U}_{{\mathrm{\Δ\δ}}_{x,i}}(k_x,k_y),$ 0≤k_x，k_y＜N；

其中，

为g(x，y)的二维奇异函数

在y方向的差分；

(14)根据以下公式分别得到G_Δ(k_x，k_y)和补零成像后的频谱数据

和

${\tilde{G}}_{\mathrm{\Δ}}(k_x,k_y)=\left\{\begin{array}{cc}{G}_{\mathrm{\Δ}}(k_x,k_y),& k_x,k_y\∈{\mathrm{\Ω}}_k\\ 0,& k_x,k_y\∈{\stackrel{\‾}{\mathrm{\Ω}}}_k\end{array}\right.,$ ${\tilde{U}}_{{\mathrm{\Δ\δ}}_{x,i}}(k_x,k_y)=\left\{\begin{array}{cc}{U}_{{\mathrm{\Δ\δ}}_{x,i}}(k_x,k_y),& k_x,k_y\∈{\mathrm{\Ω}}_k\\ 0,& k_x,k_y\∈{\stackrel{\‾}{\mathrm{\Ω}}}_k\end{array}\right.;$

(15)根据以下公式分别得到

和

${\tilde{u}}_{{\mathrm{\δ}}_x}(x-x_i,{y-y}_i)=\mathrm{IDFT}\left[{\tilde{U}}_{{\mathrm{\δ}}_{x,i}}(k_x,k_y)\right];$

$\mathrm{\Δ}{\tilde{u}}_{{\mathrm{\δ}}_x}(x-x_i,y-y_i)=\mathrm{IDFT}\left[{\tilde{U}}_{{\mathrm{\Δ\δ}}_{x,i}}(k_x,k_y)\right];$

其中，IDFT[·]为二维离散付里叶反变换算子；

(16)根据以下公式分别得到g(x，y)和g(x，y)在y方向的差分Δg(x，y)的补零成像后的复图像信号

和

$\tilde{g}(x,y)=\underset{i=1}{\overset{Q}{\mathrm{\Σ}}}{a}_i{\tilde{u}}_{{\mathrm{\δ}}_x}(x-x_i,y-y_i),$ 0≤x，y＜N；

$\mathrm{\Δ}\tilde{g}(x,y)=\underset{i=1}{\overset{Q}{\mathrm{\Σ}}}{a}_i\mathrm{\Δ}{\tilde{u}}_{{\mathrm{\δ}}_x}(x-x_i,y-y_i),$ 0≤x，y＜N。

该基于复二维奇异谱分析的磁共振部分K数据图像重建方法的进行模型参数估计包括以下步骤：

(21)根据

和

使用二维层析法进行g(x，y)的奇异点集的估计处理；

(22)根据g(x，y)的奇异点集进行相应的奇异值的估计处理；

(23)将所得到的g(x，y)的奇异点集和相应的奇异值作为模型参数估计的结果返回。

该基于复二维奇异谱分析的磁共振部分K数据图像重建方法的使用二维层析法进行g(x，y)的奇异点集的估计处理包括以下步骤：

(1′)在

中，找到

绝对值最大，即

的位置坐标(x_i，y_i)，并将(x_i，y_i)加入奇异点队列Q中；

(2′)根据以下公式计算

$\mathrm{\Δ}\tilde{g}(x,y)=\mathrm{\Δ}\tilde{g}(x,y)-\mathrm{\α\Δ}{\tilde{u}}_{{\mathrm{\δ}}_x}(x-x_i,y-y_i),$ 0≤x，y＜N；

其中 $\mathrm{\α}=\frac{\left|\mathrm{\Ω}\right|\mathrm{\Δ}\tilde{g}(x_i,y_i)}{\left|{\mathrm{\Ω}}_k\right|},$ |Ω|为空间Ω中元素的个数，|Ω_k|为空间Ω_k中元素的个数；

(3′)判断 $\underset{0\&le;x,y<N}{\max }\left(\right|\mathrm{\&Delta;}\tilde{g}(x,y)\left|\right)>T$ 是否成立，如果是，则返回上述步骤(1′)，其中T为系统中预设的与噪声相关的阈值；

(4′)反之则将该奇异点队列Q作为奇异点集{(x₁，y₁)，(x₂，y₂)，...，(x_Q，y_Q)}输出。

该基于复二维奇异谱分析的磁共振部分K数据图像重建方法的进行g(x，y)的奇异点集的奇异值的估计处理可以包括以下步骤：

(31)根据以下公式计算

$\tilde{g}(x,y)=\mathrm{IDFT}\left[\tilde{G}(k_x,k_y)\right];$

(32)根据奇异点集{(x₁，y₁)，(x₂，y₂)，...，(x_Q，y_Q)}选取奇异函数

其中i＝1～Q，并按照以下公式计算：

$U_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y)=\mathrm{DFT}\left[u_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)\right];$

${\tilde{u}}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)=\mathrm{IDFT}\left[{\tilde{U}}_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y)\right];$

(33)根据以下公式联列奇异函数方程：

$\tilde{g}(x,y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i{\tilde{u}}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i),$ 0≤x，y＜N；

(34)用伪逆矩阵法解出所述的奇异函数方程，得到一个最小误差解，获得Q个复数奇异值{a₁，a₂，...，a_Q}并输出。

该基于复二维奇异谱分析的磁共振部分K数据图像重建方法的进行g(x，y)的奇异点集的奇异值的估计处理也可以包括以下步骤：

(41)根据奇异点集{(x₁，y₁)，(x₂，y₂)，...，(x_Q，y_Q)}选取奇异函数

其中i＝1～Q；

(42)根据以下公式计算

$U_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y)=\mathrm{DFT}\left[u_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)\right],$ 其中i＝1～Q；

(43)根据以下公式联列奇异谱方程：

$G(k_x,k_y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i{U}_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y),$ k_x，k_y∈Ω_k；

(44)用伪逆矩阵法解出所述的奇异谱方程，得到一个最小误差解，获得Q个复数奇异值{a₁，a₂，..，a_Q}并输出。

该基于复二维奇异谱分析的磁共振部分K数据图像重建方法的进行磁共振复图像的重构可以为：

基于模型参数估计的结果奇异点{(x₁，y₁)，(x₂，y₂)，..，(x_Q，y_Q)}和奇异值{a₁，a₂，...，a_Q}，根据以下公式重构所述的复图像信号g(x，y)：

$g(x,y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i{u}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)$ 0≤x，y＜N；

或者也可以包括以下步骤：

(51)基于模型参数估计的结果奇异点{(x₁，y₁)，(x₂，y₂)，...，(x_Q，y_Q)}和奇异值{a₁，a₂，...，a_Q}，根据以下公式重构所述的付里叶谱数据G(k_x，k_y)：

$G(k_x,k_y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i{U}_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y),$ 0≤k_x，k_y＜N；

(52)根据以下公式得到所述的复图像信号g(x，y)：

g(x，y)＝IDFT[G(k_x，k_y)]，0≤x，y＜N。

采用了该发明的基于复二维奇异谱分析的磁共振部分K数据图像重建方法，由于首先从实际磁共振设备中根据预设的图像空间范围采集部分K空间数据，然后对该部分K数据进行补零成像处理，接着根据该经过补零成像得到的近似图像及部分K数据信息进行模型参数估计，最后根据模型参数估计的结果，利用复系数加权二维奇异函数的磁共振图像的数学模型和复二维奇异谱分析模型进行磁共振复图像的重构，从而相对于一维磁共振部分K空间数据图像重建过程大大节省了扫描时间，实现了快速成像，同时确保了图像的高信噪比、高分辨率和高精确度；而且相比较现有技术中的二维部分K空间数据图像重构方法，能够克服补零法成像中所存在的伪影，有效降低图像误差，精确显示原磁共振图像，为医学核磁共振成像检测提供了高质量的可靠图像信息；同时，本发明的方法高效实用，工作性能稳定可靠、适用范围较为广泛，给人们的工作和生活带来了很大的便利，并且也为医学成像检测技术的进一步发展和大范围普及应用奠定了坚实的理论和实践基础。

附图说明

图1为本发明的奇异点为(58，36)的奇异函数图象。

图2a、2b、2c分别为物体模拟磁共振成像试验中的二维复图像信号g(x，y)、K数据图像G(k_x，k_y)和g(x，y)在y方向的差分Δg(x，y)的图像。

图2d、2e、2f分别为使用本发明的复二维奇异谱分析(2DSSA)方法对图2a、2b、2c进行补零成像后的

和

的图像。

图3a、3b分别为使用本发明的复二维奇异谱分析(2DSSA)方法对g(x，y)的二维奇异函数

进行截断频谱的补零成像后奇异函数

的实部函数和虚部函数图像。

图3c、3d分别为使用本发明的复二维奇异谱分析(2DSSA)方法对g(x，y)的二维奇异函数在y方向的差分

进行截断频谱的补零成像后的

的实部函数和虚部函数图像。

图4为本发明的基于复二维奇异谱分析的磁共振部分K数据图像重建方法的工作过程原理示意图。

图5为仿真实验中由仿真图像调制而成的相位缓慢变化的复数图像的相位图。

图6a、6b、6c分别为仿真实验中的原始图像、使用现有技术中的FZI法对图6a的图像进行重构后的图像和使用本发明的复二维奇异谱分析(2DSSA)方法对图6a的图像进行重建后的图像。

图7为实际人体磁共振成像试验中使用本发明的复二维奇异谱分析(2DSSA)方法和现有技术中的ZFI方法分别对部分K空间数据进行重构的图像与完全K空间谱数据重构的图像之间的标准差STD随部分K数据的中心点变化状况示意图。

图8为实际人体磁共振成像试验中使用本发明的复二维奇异谱分析(2DSSA)方法和现有技术中的ZFI方法分别对三维K空间切片中部分图像数据进行重构的图像与该对应切片中完全K空间谱数据重构的图像之间的标准差STD随不同切片的变化状况示意图。

图9a为实际人体磁共振成像试验中使用三维K空间第88幅横切片中完全K空间谱数据重构的图像。

图9b为该实际人体磁共振成像试验中三维K空间第88幅横切片中部分K空间数据示意图。

图9c为使用现有技术中的FZI方法根据图9b中的数据进行重构后的图像。

图9d为使用本发明的复二维奇异谱分析(2DSSA)方法根据图9b中的数据进行重构后的图像。

具体实施方式

为了能够更清楚地理解本发明的技术内容，特举以下实施例详细说明。

本发明是从实际磁共振设备中预设的图像空间范围中采集部分K数据，并采用二维直角网格形式的部分K空间数据成像方法，因而称运用复二维奇异谱分析方法直接重构磁共振图像方法为复二维奇异谱分析方法(2DSSA，2 Dimension Singular Spectrum Analysis)。

在阐述本发明的整体工作过程及工作原理之前，为了更加明确其技术含义，首先需要给出以下定义：

定义1：给定实的或复的一个数字信号，其差分不为零的点为奇异点，奇异点上的差分值为奇异值，奇异值可以是实数也可以是复数。

定义2：实数字信号w(x)，x＝0，1，...，N-1的有一个唯一奇异点，且奇异值为1，则称w(x)为奇异函数。

为了适应二维部分K空间数据的图像重构，我们定义二维奇异函数为

$u_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)=\mathrm{\&delta;}(x-x_i)u(y-y_i)...\left(1\right)$

其中δ(x-x_i)和u(y-y_i)分别是单位脉冲和单位阶跃函数：

$\mathrm{\&delta;}(x-x_i)=\left\{\begin{array}{cc}1,& x=x_i\\ 0,& x\&NotEqual;x_i\end{array}\right.,$ $u(y-y_i)=\left\{\begin{array}{cc}1,& y\&GreaterEqual;y_i\\ 0,& y<y_i\end{array}\right..$

如奇异点为(58，36)的奇异函数请参阅图1所示。

对于任意二维像素为N×N的磁共振图像g(x，y)，0≤x，y＜N，其复二维奇异函数磁共振图像模型为：

$g(x,y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i{u}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)$ 0≤x，y＜N    ......(2)

即，g(x，y)都可以用

的复数加权和表示。其中，(x_i，y_i)和a_i分别称为第i个奇异点(以适应磁共振图像多为复数的情况)和第i个奇异点的奇异值，Q为奇异点数。对公式(2)两边对y差分得，则图像g(x，y)的y方向差分Δg(x，y)表示为：

$\mathrm{\&Delta;g}(x,y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i{\mathrm{\&Delta;u}}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)$ ......(3)

$=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i\mathrm{\&delta;}(x-x_i)\mathrm{\&delta;}(y-y_i)$ 0≤x，y＜N

由上式可见，复二维奇异函数磁共振图像模型模型参数奇异点是y方向差分不为零的位置，而奇异值就是差分值。

对公式(2)两边同取离散付里叶变换，就得磁共振K空间G(k_x，k_y)的复二维奇异谱

分析图像模型：

$G(k_x,k_y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i{U}_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y),$ 0≤k_x，k_y＜N    ......(4)

其中G(k_x，k_y)＝DFT[g(x，y)]， $U_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y)=\mathrm{DFT}\left[u_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)\right].$ 即磁共振K空间数据可由二维奇异函数的谱函数的复系数加权和表示：其中DFT[·]表示二维离散付里叶变换算子。同理，对公式(3)两边同取付里叶变换有：

$G_{\mathrm{\&Delta;}}(k_x,k_y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i{U}_{{\mathrm{\&Delta;\&delta;}}_{x,i}}(k_x,k_y),$ 0≤k_x，k_y＜N    ......(5)

其中，G_Δ(k_x，k_y)＝DFT[Δg(x，y)]， $U_{{\mathrm{\&Delta;\&delta;}}_{x,i}}(k_x,k_y)=\mathrm{DFT}\left[{\mathrm{\&Delta;u}}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)\right].$

类似地也可以定义x方向的二维奇异函数

本发明中所讨论的一切理论方法都可以类同地应用于加权和表示图像的数学模型，以后不再说明。

如果仅已知一个部分K空间数据，要获得模型的奇异点和奇异值参数一般方法是困难的。本发明的策略是先对部分K空间数据进行补零法成像，然后用层析法从这图像找出奇异点，最后用奇异谱分析法获得奇异值。

为此，本发明需要先把K空间Ω划分为有数据子空间Ω_k和无数据子空间Ω_k。在补零法成像中，将Ω_k空间数据用零代替，Ω_k空间保持不变。这样G(k_x，k_y)、

G_Δ(k_x，k_y)和采用补零法成像后的频谱数据表示为：

$\tilde{G}(k_x,k_y)=\left\{\begin{array}{cc}G(k_x,k_y),& k_x,k_y\&Element;{\mathrm{\&Omega;}}_k\\ 0,& k_x,k_y\&Element;{\stackrel{\&OverBar;}{\mathrm{\&Omega;}}}_k\end{array}\right....\left(6\right)$

${\tilde{U}}_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y)=\left\{\begin{array}{cc}{U}_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y),& k_x,k_y\&Element;{\mathrm{\&Omega;}}_k\\ 0,& k_x,k_y\&Element;{\stackrel{\&OverBar;}{\mathrm{\&Omega;}}}_k\end{array}\right....\left(7\right)$

${\tilde{G}}_{\mathrm{\&Delta;}}(k_x,k_y)=\left\{\begin{array}{cc}{G}_{\mathrm{\&Delta;}}(k_x,k_y),& k_x,k_y\&Element;{\mathrm{\&Omega;}}_k\\ 0,& k_x,k_y\&Element;{\stackrel{\&OverBar;}{\mathrm{\&Omega;}}}_k\end{array}\right....\left(8\right)$

${\tilde{U}}_{{\mathrm{\&Delta;\&delta;}}_{x,i}}(k_x,k_y)=\left\{\begin{array}{cc}{U}_{{\mathrm{\&Delta;\&delta;}}_{x,i}}(k_x,k_y),& k_x,k_y\&Element;{\mathrm{\&Omega;}}_k\\ 0,& k_x,k_y\&Element;{\stackrel{\&OverBar;}{\mathrm{\&Omega;}}}_k\end{array}\right....\left(9\right)$

相应地，可以将部分K空间数据补零法成像后重构的信号表示为：

$\tilde{g}(x,y)=\mathrm{IDFT}\left[\tilde{G}(k_x,k_y)\right]...\left(10\right)$

${\tilde{u}}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)=\mathrm{IDFT}\left[{\tilde{U}}_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y)\right]...\left(11\right)$

$\mathrm{\&Delta;}\tilde{g}(x,y)=\mathrm{IDFT}\left[{\tilde{G}}_{\mathrm{\&Delta;}}(k_x,k_y)\right]...\left(12\right)$

$\mathrm{\&Delta;}{\tilde{u}}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)=\mathrm{IDFT}\left[{\tilde{U}}_{{\mathrm{\&Delta;\&delta;}}_{x,i}}(k_x,k_y)\right]...\left(13\right)$

其中IDFT[·]表示二维离散付里叶反变换算子。对公式(4)两边同时截断频谱，并用补零法成像，则补零法的图像函数

也可以用截断频谱的补零法奇异函数

的加权和表示：

$\tilde{g}(x,y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i{\tilde{u}}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i),$ 0≤x，y＜N    ......(14)

同样对于公式(5)，也可以变成为：

$\mathrm{\&Delta;}\tilde{g}(x,y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i\mathrm{\&Delta;}{\tilde{u}}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i),$ 0≤x，y＜N    ......(15)

请参阅图2a～2f所示，其中清晰地给出了G(k_x，k_y)与

g(x，y)与

以及Δg(x，y)与

的比较。

再请参阅图3a～3d所示，其中清晰地给出了Ω_k∶0≤k_x，k_y＜72的情形下

的实部函数和虚部函数以及

的实部函数和虚部函数的图像。

接下来需要考查

请参阅图3c所示，在原奇异点位上的有最大的函数值，而纵横两方向上有吉布斯效应。而

上的吉布斯效应是由全体

的吉布斯效应共同作用的结果。所以，

最大的地方最有可以是奇异点。据此，二维层析法求取奇异点的算法如下：

第一步，在

中找到

绝对值最大，即

的位置坐标(x_i，y_i)(即奇异点)，并加入奇异点队列Q中；

第二步，根据以下公式计算：

$\mathrm{\&Delta;}\tilde{g}(x,y)=\mathrm{\&Delta;}\tilde{g}(x,y)-\mathrm{\&alpha;\&Delta;}{\tilde{u}}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i),$ 0≤x，y＜N......(16)

其中 $\mathrm{\&alpha;}=\frac{\left|\mathrm{\&Omega;}\right|\mathrm{\&Delta;}\tilde{g}(x_i,y_i)}{\left|{\mathrm{\&Omega;}}_k\right|},$ |Ω|和|Ω_k|分别表示空间Ω和Ω_k中元素的个数；

第三步，判断 $\underset{0\&le;x,y<N}{\max }\left(\right|\mathrm{\&Delta;}\tilde{g}(x,y)\left|\right)>T$ 是否成立，其中T为系统中预设的与噪声相关的给定的阈值，如果成立，则返回第一步继续找奇异点；反之将该奇异点队列Q作为奇异点集{(x₁，y₁)，(x₂，y₂)，...，(x_Q，y_Q)}输出。

经过上述的求取奇异点的步骤后，求取该奇异点集所对应的奇异值的方法有以下两种：

(1)在图像空间中，根据上述的g(x，y)的奇异点集就确定了关于g(x，y)的补零成像后的奇异函数

集，奇异点集的奇异值可根据公式(14)或(15)确定。其算法如下：

第一步，求取 $\tilde{g}(x,y)=\mathrm{IDFT}\left[\tilde{G}(k_x,k_y)\right];$

第二步，根据奇异点集{(x₁，y₁)，(x₂，y₂)，...，(x_Q，y_Q)}选取奇异函数，i＝1～Q，并按照以下公式计算：

$U_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y)=\mathrm{DFT}\left[u_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)\right],$

${\tilde{u}}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)=\mathrm{IDFT}\left[{\tilde{U}}_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y)\right];$

第三步，联列奇异函数方程：

$\tilde{g}(x,y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i{\tilde{u}}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i),$ 0≤x，y＜N

用伪逆矩阵法解出所述的奇异函数方程，得到一个最小误差解，确定奇异值{a₁，a₂，...，a_Q}。

(2)在K空间中，根据上述的g(x，y)的奇异点集就确定了关于g(x，y)的奇异谱函数

集，利用公式(15)可计算出奇异点集对应的奇异值，其算法如下：

第一步，按奇异点集{(x₁，y₁)，(x₂，y₂)，...，(x_Q，y_Q)}选取奇异函数

第二步，i＝1～Q进行计算：

$U_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y)=\mathrm{DFT}\left[u_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)\right];$

第三步，联列奇异谱方程：

$G(k_x,k_y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i{U}_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y),$ k_x，k_y∈Ω_k    ......(17)

用伪逆矩阵法解出所述的奇异谱方程，得到一个最小误差解，确定奇异值{a₁，a₂，...，a_Q}。

最后是根据上述的奇异点和奇异值信息进行图像重构，有以下两种方法：

(1)根据奇异点{(x₁，y₁)，(x₂，y₂)，...，(x_Q，y_Q)}和奇异值{a₁，a₂，...，a_Q}，按公式(2)重构图像g(x，y)：

$g(x,y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i{u}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)$ 0≤x，y＜N

(2)按公式(4)重构未采集的K空间数据G(k_x，k_y)：

$G(k_x,k_y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i{U}_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y),$ 0≤k_x，k_y＜N；

然后再用离散傅里叶反变换成像：

g(x，y)＝IDFT[G(k_x，k_y)]，0≤x，y＜N。

请参阅图4所示，该基于复二维奇异谱分析的磁共振部分K数据图像重建方法，包括以下步骤：

(1)从磁共振成像扫描仪中预设的图像空间范围中采集部分K数据G(k_x，k_y)；

(2)对该部分K数据进行补零成像处理，包括以下步骤：

(a)将K空间Ω划分为有数据子空间Ω_k和无数据子空间Ω_k；

(b)将Ω_k子空间中的数据用零代替，Ω_k空间保持不变，并根据以下公式分别得到G(k_x，k_y)和

补零成像后的频谱数据

和

$\tilde{G}(k_x,k_y)=\left\{\begin{array}{cc}G(k_x,k_y),& k_x,k_y\&Element;{\mathrm{\&Omega;}}_k\\ 0,& k_x,k_y\&Element;{\stackrel{\&OverBar;}{\mathrm{\&Omega;}}}_k\end{array}\right.,$ ${\tilde{U}}_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y)=\left\{\begin{array}{cc}{U}_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y),& k_x,k_y\&Element;{\mathrm{\&Omega;}}_k\\ 0,& k_x,k_y\&Element;{\stackrel{\&OverBar;}{\mathrm{\&Omega;}}}_k\end{array}\right.;$

(c)根据以下公式计算G_Δ(k_x，k_y)和

$U_{{\mathrm{\&Delta;\&delta;}}_{x,i}}(k_x,k_y)=\mathrm{DFT}\left[{\mathrm{\&Delta;u}}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)\right];$

$G_{\mathrm{\&Delta;}}(k_x,k_y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i{U}_{{\mathrm{\&Delta;\&delta;}}_{x,i}}(k_x,k_y),$ 0≤k_x，k_y＜N；

其中，

为g(x，y)的二维奇异函数

在y方向的差分；

(d)根据以下公式分别得到G_Δ(k_x，k_y)和

补零成像后的频谱数据

和

${\tilde{G}}_{\mathrm{\&Delta;}}(k_x,k_y)=\left\{\begin{array}{cc}{G}_{\mathrm{\&Delta;}}(k_x,k_y),& k_x,k_y\&Element;{\mathrm{\&Omega;}}_k\\ 0,& k_x,k_y\&Element;{\stackrel{\&OverBar;}{\mathrm{\&Omega;}}}_k\end{array}\right.,$ ${\tilde{U}}_{{\mathrm{\&Delta;\&delta;}}_{x,i}}(k_x,k_y)=\left\{\begin{array}{cc}{U}_{{\mathrm{\&Delta;\&delta;}}_{x,i}}(k_x,k_y),& k_x,k_y\&Element;{\mathrm{\&Omega;}}_k\\ 0,& k_x,k_y\&Element;{\stackrel{\&OverBar;}{\mathrm{\&Omega;}}}_k\end{array}\right.;$

(e)根据以下公式分别得到

和

${\tilde{u}}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)=\mathrm{IDFT}\left[{\tilde{U}}_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y)\right];$

$\mathrm{\&Delta;}{\tilde{u}}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)=\mathrm{IDFT}\left[{\tilde{U}}_{{\mathrm{\&Delta;\&delta;}}_{x,i}}(k_x,k_y)\right];$

其中，IDFT[·]为二维离散付里叶反变换算子；

(f)根据以下公式分别得到g(x，y)和g(x，y)在y方向的差分Δg(x，y)的补零成像后的复图像信号

和

$\tilde{g}(x,y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i{\tilde{u}}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i),$ 0≤x，y＜N；

$\mathrm{\&Delta;}\tilde{g}(x,y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i{\mathrm{\&Delta;}\tilde{u}}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i),$ 0≤x，y＜N；

(3)根据该经过补零成像处理的部分K数据信息进行模型参数估计，包括以下步骤：

(a)根据

和

使用二维层析法进行g(x，y)的奇异点集的估计处理，包括以下步骤：

(i)在

中，找到绝对值最大，即

的位置坐标(x_i，y_i)，并将(x_i，y_i)加入奇异点队列Q中；

(ii)根据以下公式计算

$\mathrm{\&Delta;}\tilde{g}(x,y)=\mathrm{\&Delta;}\tilde{g}(x,y)-\mathrm{\&alpha;\&Delta;}{\tilde{u}}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i),$ 0≤x，y＜N；

其中 $\mathrm{\&alpha;}=\frac{\left|\mathrm{\&Omega;}\right|\mathrm{\&Delta;}\tilde{g}(x_i,y_i)}{\left|{\mathrm{\&Omega;}}_k\right|},$ |Ω|为空间Ω中元素的个数，|Ω_k|为空间Ω_k中元素的个数；

(iii)判断 $\underset{0\&le;x,y<N}{\max }\left(\right|\mathrm{\&Delta;}\tilde{g}(x,y)\left|\right)>T$ 是否成立，如果是，则返回上述步骤(i)，其中T为系统中预设的与噪声相关的阈值；

(iv)反之则将该奇异点队列Q作为奇异点集{(x₁，y₁)，(x₂，y₂)，...，(x_Q，y_Q)}输出。

(b)根据g(x，y)的奇异点集进行相应的奇异值的估计处理，可以包括以下步骤：

(i)根据以下公式计算

$\tilde{g}(x,y)=\mathrm{IDFT}\left[\tilde{G}(k_x,k_y)\right];$

(ii)根据奇异点集{(x₁，y₁)，(x₂，y₂)，..，(x_Q，y_Q)}选取奇异函数其中i＝1～Q，并按照以下公式计算：

$U_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y)=\mathrm{DFT}\left[u_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)\right];$

${\tilde{u}}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)=\mathrm{IDFT}\left[{\tilde{U}}_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y)\right];$

(iii)根据以下公式联列奇异函数方程：

$\tilde{g}(x,y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i{\tilde{u}}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i),$ 0≤x，y＜N；

(iv)用伪逆矩阵法解出所述的奇异函数方程，得到一个最小误差解，获得Q个复数奇异值{a₁，a₂，...，a_Q}并输出；也可以包括以下步骤：

(i)根据奇异点集{(x₁，y₁)，(x₂，y₂)，...，(x_Q，y_Q)}选取奇异函数

其中i＝1～Q；

(ii)根据以下公式计算

$U_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y)=\mathrm{DFT}\left[u_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)\right],$ 其中i＝1～Q；

(iii)根据以下公式联列奇异谱方程：

$G(k_x,k_y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i{U}_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y),$ k_x，k_y∈Ω_k；

(iv)用伪逆矩阵法解出所述的奇异谱方程，得到一个最小误差解，获得Q个复数奇异值{a₁，a₂，...，a_Q}并输出；

(c)将所得到的g(x，y)的奇异点集和相应的奇异值作为模型参数估计的结果返回；

(4)根据模型参数估计的结果，利用复系数加权二维奇异函数的磁共振图像的数学模型和复二维奇异谱分析模型进行磁共振复图像的重构；其中，该复系数加权二维奇异函数的磁共振图像数学模型为：

$g(x,y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i{u}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)$ 0≤x，y＜N；

其中，g(x，y)，0≤x，y＜N为像素为N×N的二维磁共振图像的复图像信号，(x_i，y_i)为g(x，y)的第i个奇异点，a_i为该第i个奇异点上的复奇异值，Q为g(x，y)的奇异点的个数，

为以(x_i，y_i)为奇异点的二维奇异函数；

该复二维奇异谱分析模型为：

$G(k_x,k_y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i{U}_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y),$ 0≤k_x，k_y＜N；

其中，G(k_x，k_y)＝DFT[g(x，y)]， $U_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y)=\mathrm{DFT}\left[u_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)\right],$ DFT[·]为二维离散付里叶变换算子；

该磁共振复图像的重构可以为：

基于模型参数估计的结果奇异点{(x₁，y₁)，(x₂，y₂)，...，(x_Q，y_Q)}和奇异值{a₁，a₂，...，a_Q}，根据以下公式重构所述的复图像信号g(x，y)：

$g(x,y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i{u}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)$ 0≤x，y＜N；

或者也可以包括以下步骤：

(a)基于模型参数估计的结果奇异点{(x₁，y₁)，(x₂，y₂)，...，(x_Q，y_Q)}和奇异值{a₁，a₂，...，a_Q}，根据以下公式重构所述的付里叶谱数据G(k_x，k_y)：

$G(k_x,k_y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i{U}_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y),$ 0≤k_x，k_y＜N；

(b)根据以下公式得到所述的复图像信号g(x，y)：

g(x，y)＝IDFT[G(k_x，k_y)]，0≤x，y＜N。

以下是使用本发明的方法进行的仿真实验：

首先应用计算机仿真数据进行算法测试，仿真实验用的图像是一幅灰度范围为0～255，图像尺寸为128×128。考虑到磁共振图像大多数情况下，在图像区域的相位是缓慢变化，变化范围一般在[0°，360°]之内(如果相位变化频率超过这个范围，可以通过K空间中心点平移方法进行校正)；所以，实验中，取相位变化范围为0°～360°。

仿真实验一：噪声的敏感性实验

将仿真图像调制上成一幅相位缓慢变化的复数图像，其相位图请参阅图5所示，再加入0均值的标准差分别是1、2、3、4、5、6、7、8、9的高斯白噪声，作为原始图像。然后对这原始图像进行傅里叶变换，再取K空间数据Ω_k＝{-24＜k_x，k_y＜40}，用本发明的方法进行图像重构，最后计算重构图像与原始图像的标准差。实验结果如表1所示。

STD	1	2	3	4	5	6	7	8	9
										2DSSA	1.3582	2.3576	3.6088	4.2717	4.3534	5.5481	6.0089	7.2448	7.6032
FZI	9.9308	9.9842	9.9928	10.051	10.142	10.148	10.187	10.315	10.332

表1.随噪声变化的重构精度比较

从表1中的数据结果可以看出，本发明的2DSSA方法重构图像的标准差(STD)逼近加入噪声的STD，噪声对本发明的2DSSA算法影响很小，并且对于高噪声图像有提高信噪比的功效。这主要是由于2DSSA通过层析法以避免因噪声引入虚假奇异点，通过求解伪逆方程进一步滤去层析中混入的虚假奇异点，从而大大抑制了噪声对2DSSA的影响。2DSSA因截去含噪声分量高的高频K空间，而后又重构还原了高频信号分量，从而提高了信噪比。而在现有技术的FZI方法中，在低噪声段，其STD主要来自于截断伪迹影响，在高噪声段受到截断伪迹和噪声的双重影响，所以都有较高的STD，重构质量相对较差。

仿真实验二：复数图像相位变化敏感性实验

在上述的仿真图像上加入标准差为5的均值为0的高斯噪声，对加噪声后的仿真图像进行相位调制构成测试用的原始图像，调制的相位变化范围分别是40°、80°、120°、160°、200°、240°、280°、320°、360°，然后对这原始图像进行傅里叶变换，再取K空间数据Ω_k＝{-24＜k_x，k_y＜40}，用本发明的方法进行图像重构，最后计算重构图像与原始图像的标准差。实验结果如表2所示。

PHASE	40	80	120	160	200	240	280	320	360
										2DSSA	4.6015	4.5875	4.5502	4.7380	5.1505	4.7065	4.115	4.2908	4.3534
FZI	10.157	10.025	10.167	10.604	10.125	10.213	10.213	10.186	10.142

表2.随图像相位变化的重构精度比较

由表2的数据可知，缓慢的相位变化对2DSSA重构精度没有任何相关性，也就是本发明的2DSSA方法和现有技术的ZFI都不受复数图像的缓慢相位变化的影响。由于缓慢变化的相位变化产生图像差分大小，相对于噪声产生的差分变化来说要小轻微得多。在层析法、伪逆求解奇异值方法的压制虚假奇异点的作用下，相位变化不易产生对2DSSA重构精度的影响。但是实际的相位变化并非如上述实验中那么简单，相位突变情况时有发生。

请参阅图6所示，其中给出一个仿真图像重构的例子。Ω_k＝{-24＜k_x，k_y＜40}，其加性高斯零均值噪声的STD为5。在FZI方法的图中，截断伪迹随处可见，2DSSA方法几乎没有截断伪迹，与原始图像比较仅仅有部分地方的细节变得略微模糊。

以下是使用本发明的方法进行实际MRI数据的测试实验。

实验三：分析算法的部分K数据采集的空间差异性

实验用的实际磁共振图像是一幅灰度范围为0～255，图像尺寸为256×176。用于实验的部分K空间大小为128×88。实验过程设计如下：

让部分K数据空间的中心点，从原点开始，以水平方向步距4，纵向步距6，向右下角移动，并每移动一次用本发明的2DSSA方法和现有技术的ZFI方法重构图像，把重构的图像和完全K空间数据重构的图像，进行比较，并给出了标准差STD随部分K数据的中心点变化情况，请参阅图7所示。从图7中可以发现：

第一，本发明的2DSSA方法比较现有技术的ZFI有好得多的重构精度，这说明2DSSA能较好重构未取到的部分频谱数据；

第二，在同样部分K空间数据条件下，采集空间的中心点位置稍偏离原点，对本发明的2DSSA方法的提取奇异和奇异值有益，但不宜过分不对称，否则会使频谱能量丢失过多，而引入较大误差。

实验四：2DSSA对图像解剖结构的敏感性

实际磁共振数据卷为256×256×276的三维K空间数据，实验设计如下：

根据上述实验三的结果，取部分K数据采集范围为Ω_k＝{-28＜k_x＜60，-42＜k_y＜86}，然对各片图像使用本发明的2DSSA方法和现有技术的ZFI方法进行重构，最后将重构图像与完全谱数据重构图像进行比较，并给出了随不同横切片的STD变化情况，结果如图8所示，其中表明本发明的2DSSA方法和现有技术的ZFI方法一样受图像解剖结构的影响，但2DSSA方法总是有比ZFI有更高重构图像精度。

为了从视觉直观地观察本方法效果，请参阅图9a～9d所示，其中给出了图像卷中第88幅横切面图像，其中图9a为完全K空间数据所作的图像，图9b为部分K空间示意图，图9c为根据图9b中的部分K空间的数据用现有技术中的ZFI方法所重构的图像，图9d为根据图9b中的部分K空间的数据用本发明的2DSSA方法所重构的图像。从图9c和图9d的直观比较可以看出，在现有技术的ZFI方法中，截断伪迹随处可见，而本发明的2DSSA方法的图像中几乎难寻这种截断伪迹，可以说本发明的2DSSA方法是行之有效的方法。

从以上可以看出，本发明的方法的基本思想是：首先给出二维奇异谱分析图像重构模型，为解决磁共振图像的相位问题引入复数加权系数，运用层析法、伪逆矩阵确定奇异值方法，从而能够较好的抑制噪声，并消除了相位变化对重构图像质量所造成的的不利影响。在仿真实验和实际磁共振图像数据的测试实验中，本发明的2DSSA方法都表现出比现有技术ZFI好得多的结果。本发明的2DSSA方法把目前研究的一维部分K空间数据重构问题推广到二维中去，这将对解决部分K空间磁共振数据的图像重建问题提供了一种新的思考方法。

采用了上述的基于复二维奇异谱分析的磁共振部分K数据图像重建方法，由于首先从实际磁共振设备中根据预设的图像空间范围采集部分K空间数据，然后对该部分K数据进行补零成像处理，接着根据补零成像得到的近似图像及部分K数据进行模型参数估计，最后根据模型参数估计的结果，利用复系数加权二维奇异函数的磁共振图像的数学模型和复二维奇异谱分析模型进行磁共振复图像的重构，从而相对于一维磁共振部分K空间数据图像重建过程大大节省了扫描时间，实现了快速成像，同时确保了图像的高信噪比、高分辨率和高精确度；而且相比较现有技术中的二维部分K空间数据图像重构方法，能够克服补零法成像中所存在的伪影，有效降低图像误差，精确显示原磁共振图像，为医学核磁共振成像检测提供了高质量的可靠图像信息；同时，本发明的方法高效实用，工作性能稳定可靠、适用范围较为广泛，给人们的工作和生活带来了很大的便利，并且也为医学成像检测技术的进一步发展和大范围普及应用奠定了坚实的理论和实践基础。

在此说明书中，本发明已参照其特定的实施例作了描述。但是，很显然仍可以作出各种修改和变换而不背离本发明的精神和范围。因此，说明书和附图应被认为是说明性的而非限制性的。

Claims (9)

1、一种基于复二维奇异谱分析的磁共振部分K数据图像重建方法，其特征在于，所述的方法包括以下步骤：

(1)从磁共振成像扫描仪中预设的图像空间范围中采集部分K数据G(k_x，k_y)，所述的K数据G(k_x，k_y)即为离散傅里叶变换空间数据G(k_x，k_y)；

(2)对该部分K数据进行补零成像处理；

(3)根据补零成像得到的近似图像及部分K数据进行模型参数估计；

(4)根据模型参数估计的结果，利用复系数加权二维奇异函数的磁共振图像的数学模型和复二维奇异谱分析模型进行磁共振复图像的重构。

2、根据权利要求1所述的基于复二维奇异谱分析的磁共振部分K数据图像重建方法，其特征在于，所述的复系数加权二维奇异函数的磁共振图像数学模型为：

$g(x,y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i{u}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)$ 0≤x，y＜N；

其中，g(x，y)，0≤x，y＜N为像素为N×N的二维磁共振图像的复图像信号，(x_i，y_i)为g(x，y)的第i个奇异点，a_i为该第i个奇异点上的复奇异值，Q为g(x，y)的奇异点的个数，为以(x_i，y_i)为奇异点的二维奇异函数。

3、根据权利要求2所述的基于复二维奇异谱分析的磁共振部分K数据图像重建方法，其特征在于，所述的复二维奇异谱分析模型为：

$G(k_x,k_y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i{U}_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y),$ 0≤k_x，k_y＜N；

其中，G(k_x，k_y)＝DFT[g(x，y)]， $U_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y)=\mathrm{DFT}\left[u_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)\right],$ DFT[·]为二维离散付里叶变换算子。

4、根据权利要求3所述的基于复二维奇异谱分析的磁共振部分K数据图像重建方法，其特征在于，所述的进行补零成像处理包括以下步骤：

(11)将K空间Ω划分为有数据子空间Ω_k和无数据子空间Ω_k；

(12)将Ω_k子空间中的数据用零代替，Ω_k空间保持不变，并根据以下公式分别得到G(k_x，k_y)和补零成像后的频谱数据和

$\tilde{G}(k_x,k_y)=\left\{\begin{array}{cc}G(k_x,k_y),& k_x,k_y\&Element;{\mathrm{\&Omega;}}_k\\ 0,& k_x,k_y\&Element;{\stackrel{\&OverBar;}{\mathrm{\&Omega;}}}_k\end{array}\right.,$ ${\tilde{U}}_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y)=\left\{\begin{array}{cc}{U}_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y),& k_x,k_y\&Element;{\mathrm{\&Omega;}}_k\\ 0,& k_x,k_y\&Element;{\stackrel{\&OverBar;}{\mathrm{\&Omega;}}}_k\end{array}\right.;$

(13)根据以下公式计算G_Δ(k_x，k_y)和

$U_{{\mathrm{\&Delta;\&delta;}}_{x,i}}(k_x,k_y)=\mathrm{DFT}\left[{\mathrm{\&Delta;u}}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)\right];$

$G_{\mathrm{\&Delta;}}(k_x,k_y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i{U}_{{\mathrm{\&Delta;\&delta;}}_{x,i}}(k_x,k_y),$ 0≤k_x，k_y＜N；

其中，

为g(x，y)的二维奇异函数

在y方向的差分；

(14)根据以下公式分别得到G_Δ(k_x，k_y)和

补零成像后的频谱数据

和

${\tilde{G}}_{\mathrm{\&Delta;}}(k_x,k_y)=\left\{\begin{array}{cc}{G}_{\mathrm{\&Delta;}}(k_x,k_y),& k_x,k_y\&Element;{\mathrm{\&Omega;}}_k\\ 0,& k_x,k_y\&Element;{\stackrel{\&OverBar;}{\mathrm{\&Omega;}}}_k\end{array}\right.,$ ${\tilde{U}}_{{\mathrm{\&Delta;\&delta;}}_{x,i}}(k_x,k_y)=\left\{\begin{array}{cc}{U}_{{\mathrm{\&Delta;\&delta;}}_{x,i}}(k_x,k_y),& k_x,k_y\&Element;{\mathrm{\&Omega;}}_k\\ 0,& k_x,k_y\&Element;{\stackrel{\&OverBar;}{\mathrm{\&Omega;}}}_k\end{array}\right.;$

(15)根据以下公式分别得到

和

${\tilde{u}}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)=\mathrm{IDFT}\left[{\tilde{U}}_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y)\right];$

$\mathrm{\&Delta;}{\tilde{u}}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)=\mathrm{IDFT}\left[{\tilde{U}}_{{\mathrm{\&Delta;\&delta;}}_{x,i}}(k_x,k_y)\right];$

其中，IDFT[·]为二维离散付里叶反变换算子；

(16)根据以下公式分别得到g(x，y)和g(x，y)在y方向的差分Δg(x，y)的补零成像后的复图像信号

和

$\tilde{g}(x,y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i{\tilde{u}}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i),$ 0≤x，y＜N；

$\mathrm{\&Delta;}\tilde{g}(x,y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i\mathrm{\&Delta;}{\tilde{u}}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i),$ 0≤x，y＜N。

5、根据权利要求4所述的基于复二维奇异谱分析的磁共振部分K数据图像重建方法，其特征在于，所述的进行模型参数估计包括以下步骤：

(21)根据

和

使用二维层析法进行g(x，y)的奇异点集的估计处理；

(22)根据g(x，y)的奇异点集，构造二维奇异谱方程，进行相应的奇异值的估计处理；

(23)将所得到的g(x，y)的奇异点集和相应的奇异值作为模型参数估计的结果返回。

6、根据权利要求5所述的基于复二维奇异谱分析的磁共振部分K数据图像重建方法，其特征在于，所述的使用二维层析法进行g(x，y)的奇异点集的估计处理包括以下步骤：

(1′)在

中，找到

绝对值最大，即 $\underset{0\&le;x,y<N}{\max }\left(\right|\mathrm{\&Delta;}\tilde{g}(x,y)\left|\right)$ 的位置坐标(x_i，y_i)，并将(x_i，y_i)加入奇异点队列Q中；

(2′)根据以下公式计算

$\mathrm{\&Delta;}\tilde{g}(x,y)=\mathrm{\&Delta;}\tilde{g}(x,y)-\mathrm{\&alpha;\&Delta;}{\tilde{u}}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i),$ 0≤x，y＜N；

其中 $\mathrm{\&alpha;}=\frac{\left|\mathrm{\&Omega;}\right|\mathrm{\&Delta;}\tilde{g}(x_i,y_i)}{\left|{\mathrm{\&Omega;}}_k\right|},$ |Ω|为空间Ω中元素的个数，|Ω_k|为空间Ω_k中元素的个数；

(3′)判断 $\underset{0\&le;x,y<N}{\max }\left(\right|\mathrm{\&Delta;}\tilde{g}(x,y)\left|\right)>T$ 是否成立，如果是，则返回上述步骤(1′)，其中T为系统中预设的与噪声相关的阈值；

(4′)反之则将该奇异点队列Q作为奇异点集{(x₁，y₁)，(x₂，y₂)，...，(x_Q，y_Q)}输出。

7、根据权利要求6所述的基于复二维奇异谱分析的磁共振部分K数据图像重建方法，其特征在于，所述的进行g(x，y)的奇异点集的奇异值的估计处理包括以下步骤：

(31)根据以下公式计算

$\tilde{g}(x,y)=\mathrm{IDFT}\left[\tilde{G}(k_x,k_y)\right];$

(32)根据奇异点集{(x₁，y₁)，(x₂，y₂)，...，(x_Q，y_Q)}选取奇异函数

其中i＝1～Q，并按照以下公式计算：

$U_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y)=\mathrm{DFT}\left[u_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)\right];$

${\tilde{u}}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)=\mathrm{IDFT}\left[{\tilde{U}}_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y)\right];$

(33)根据以下公式联列奇异函数方程：

$\tilde{g}(x,y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i{\tilde{u}}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i),$ 0≤x，y＜N；

(34)用伪逆矩阵法解出所述的奇异函数方程，得到一个最小误差解，获得Q个复数奇异值{a₁，a₂，...，a_Q}并输出。

8、根据权利要求6所述的基于复二维奇异谱分析的磁共振部分K数据图像重建方法，其特征在于，所述的进行g(x，y)的奇异点集的奇异值的估计处理包括以下步骤：

(41)根据奇异点集{(x₁，y₁)，(x₂，y₂)，...，(x_Q，y_Q)}选取奇异函数

其中i＝1～Q；

(42)根据以下公式计算

$U_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y)=\mathrm{DFT}\left[u_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)\right],$ 其中i＝1～Q；

(43)根据以下公式联列奇异谱方程：

$G(k_x,k_y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i{U}_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y),$ k_x，k_y∈Ω_k；

(44)用伪逆矩阵法解出所述的奇异谱方程，得到一个最小误差解，获得Q个复数奇异值{a₁，a₂，...，a_Q}并输出。

9、根据权利要求7或8所述的基于复二维奇异谱分析的磁共振部分K数据图像重建方法，其特征在于，所述的进行磁共振复图像的重构为：

基于模型参数估计的结果奇异点{(x₁，y₁)，(x₂，y₂)，...，(x_Q，y_Q)}和奇异值{a₁，a₂，...，a_Q}，根据以下公式重构所述的复图像信号g(x，y)：

$g(x,y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i{u}_{{\mathrm{\&delta;}}_x}(x-x_i,y-y_i)$ 0≤x，y＜N；

或者包括以下步骤：

(51)基于模型参数估计的结果奇异点{(x₁，y₁)，(x₂，y₂)，...，(x_Q，y_Q)}和奇异值{a₁，a₂，...，a_Q}，根据以下公式重构所述的付里叶谱数据G(k_x，k_y)：

$G(k_x,k_y)=\underset{i=1}{\overset{Q}{\mathrm{\&Sigma;}}}{a}_i{U}_{{\mathrm{\&delta;}}_{x,i}}(k_x,k_y),$ 0≤k_x，k_y＜N；

(52)根据以下公式得到所述的复图像信号g(x，y)：

g(x，y)＝IDFT[G(k_x，k_y)]，0≤x，y＜N。

Images