CN110044937B - 一种基于压缩感知的ct成像方法和装置 - Google Patents

Abstract

本发明涉及一种基于压缩感知的CT成像方法和装置，方法具体包括以下步骤：CT机采用角度压缩比进行物体扫描，获取投影数据；采集投影数据输入至CT逆成像模型；基于ADMM算法求解CT逆成像模型；根据求解结果构建二维图像。与现有技术相比，本发明基于lasso模型构建了全新的CT逆成像模型，因为该模型具有稀疏性，所以CT机能够采用一定的角度压缩比进行投影采集，显著提高了CT机的扫描速度，减少数据的采集量，使其能够满足运动性大的器官的扫描。同时，本发明采用ADMM算法对CT逆成像模型进行求解，减少了数据计算量，减轻计算机系统的计算负荷，极大提高了成像速度。

Description

一种基于压缩感知的CT成像方法和装置

技术领域

本发明涉及一种CT成像技术领域，尤其是涉及一种一种基于压缩感知的CT成像方法和装置

背景技术

CT(Computed Tomography)可以在不破坏样品的情况下，利用样品对X射线能量的吸收特性，对生物组织和工程材料的样品进行成像来获取样品内部的结构信息。当X射线扫描人体组织时，部分射线被组织吸收，部分射线穿过人体被检测器接收，形成接受信号。因为人体各种组织的疏密程度不同，X射线的穿透能力不同，穿透人体组织的距离不同，所以探测器接收到的射线就有了差异。将所接收到的信号，转换为数字信息后计算机进行处理，显示出断层图像。

CT机采取旋转/平移方式，在病人的两侧分别装有一个X线管和一个晶体检测器，在对X线管所产生的光束进行准直后，使光束与检测器成为一条连线，当作扫描运动时X线管和检测器作平移运动扫描，然后，X线管和检测器绕病人旋转1°，再作一次平移扫描，如此重复，总共旋转180°，这样，就能从不同的方向上采集某一扫描部位的投影信号。但是现有的CT机存在以下问题：

(1)随着CT机在医学上的广泛应用，由于高剂量X射线对人体的照射，可能引起基因突变、诱发癌症等。但降低剂量就会增加CT成像的统计噪声，对诊断造成困难。

(2)CT成像的扫描过程需要花费时间，要3.5～6分钟才完成一次扫描，同时成像时间长，无法对运动性大的器官进行扫描。

(3)扫描数据量大，对计算机内存要求较高。

(4)减少扫描角度是解决上述问题的合理方法，但传统的重构算法(例如滤波反投影算法)重构的图像退化严重，未能考虑到数据的稀疏性等。

发明内容

本发明的目的就是为了克服上述现有技术存在的缺陷而提供一种基于压缩感知的CT成像方法和装置。

本发明的目的可以通过以下技术方案来实现：

一种基于压缩感知的CT成像方法，具体包括以下步骤：

S1、CT机采用角度压缩比进行物体扫描，获取投影数据；

S2、采集投影数据输入至CT逆成像模型，该CT逆成像模型表达式如下：

式中，E为第二矩阵，x为求解量，b为投影数据中的射线吸收量，δ为预设参数；

S3、基于ADMM算法求解CT逆成像模型；

S4、根据求解结果构建二维图像。

进一步地，所述的第二矩阵E获取步骤如下：

A1、根据CT机参数获取第一矩阵D(α,k)，其中，α为CT机的扫描角度，k为CT机发出射线的数量，第一矩阵D(α,k)表示α角度,第k条射线穿过扫描区域距离矩阵；

A2、将矩阵D(α,k)按列连接转换成行向量d(α,k)；

A3、按照α和k的顺序，将行向量d(α,k)重新组合成第二矩阵E，其表达式为：

进一步地，步骤S3中，基于ADMM算法求解CT逆成像模型步骤为：

首先，将CT逆成像模型表达式转化为：

s·t·x＝y；

x，y∈Rⁿ；

然后采用如下迭代公式进行求解：

式中，ρ为惩罚参数，τ为正则化变量。

进一步地，步骤S1中，角度压缩比大于等于4。

一种CT成像装置，所述的装置包括处理器以及存储器，所述处理器调用存储器中的程序，用于实现以下步骤：

S1、CT机采用角度压缩比进行物体扫描，获取投影数据；

S2、采集投影数据输入至CT逆成像模型，该CT逆成像模型表达式如下：

式中，E为第二矩阵，x为求解量，b为投影数据中的射线吸收量，δ为预设参数；

S3、基于ADMM算法求解CT逆成像模型；

S4、根据求解结果构建二维图像。

进一步地，所述的第二矩阵E获取步骤如下：

A1、根据CT机参数获取第一矩阵D(α,k)，其中，α为CT机的扫描角度，k为CT机发出射线的数量，第一矩阵D(α,k)表示α角度,第k条射线穿过扫描区域距离矩阵；

A2、将矩阵D(α,k)按列连接转换成行向量d(α,k)；

A3、按照α和k的顺序，将行向量d(α,k)重新组合成第二矩阵E，其表达式为：

进一步地，步骤S3中，基于ADMM算法求解CT逆成像模型步骤为：

首先，将CT逆成像模型表达式转化为：

s.t.x＝y；

x,y∈Rⁿ；

然后采用如下迭代公式进行求解：

式中，ρ为惩罚参数，τ为正则化变量。

进一步地，步骤S1中，角度压缩比大于等于4。

与现有技术相比，本发明具有以下优点：

(1)本发明利用CT投影信号的稀疏性，构建了基于lasso模型全新的CT逆成像模型。

(2)该模型使得能够采用一定的角度压缩比进行投影采集，显著提高了CT机的扫描速度，减少数据的采集量，减小了射线可能对人体的伤害，使其能够满足运动性大的器官的扫描。

(3)对采集数据预先进行筛选后，仅对有效数据进行处理，减少了数据计算量，减轻计算机系统的计算负荷。

(4)本发明采用ADMM算法对CT逆成像模型进行求解。由于该算法在求解CT成像模型具有解析解，所以极大提高了计算速度。

附图说明

图1为本发明的流程示意图。

图2为CT成像原理示意图。

图3为扫描区域划分示意图。

具体实施方式

下面结合附图和具体实施例对本发明进行详细说明。本实施例以本发明技术方案为前提进行实施，给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述的实施例。

本实施例提供了一种基于压缩感知的CT成像方法，具体包括以下步骤：

S1、CT机采用大于等于4的角度压缩比进行物体扫描，获取投影数据；

S2、采集投影数据输入至CT逆成像模型，该CT逆成像模型表达式如下：

式中，E为第二矩阵，x为求解量，b为投影数据中的射线吸收量，δ为预设参数；

S3、基于ADMM算法求解CT逆成像模型：

首先，将CT逆成像模型表达式转化为：

s.t.x＝y；

x,y∈Rⁿ；

然后采用如下迭代公式进行求解：

式中，ρ为惩罚参数，τ为正则化变量；

S4、根据求解结果采用现有的算法构建二维图像。

该方法的理论说明如下：

一、CT成像的数学描述

如图1所示，CT机的工作原理为：平行入射的X射线垂直于探测器平面，等距排列的每一个探测器单元都是接收器。X射线的发射器和探测器相对位置固定不变，整个CT机绕某个固定旋转中心旋转180(α)次。每旋转一次都会得到521(k)组等距单元探测器接收到的衰减后X射线，通过增益处理得到180组信息，再通过计算机对采集到的180组信息进行处理、运算和反演，从而得到物体的横断图。

CT机从本质意义上看就是建立CT系统的数学模型，根据接收信息来求得未知介质的对X射线的吸收状态，判断未知介质是否发生病变等。CT机实现上述功能的关键在于求得函数关系，建立成像模型。

如图2所示，将含有该介质的正方形扫描区域划分为若干个小的像素方格，可看作矩阵X，即：

其中，x_ij表示第i行第j列的像素格的吸收率。

CT机的任意一束X射线在穿过介质时，关注该射线穿过的像素格的距离。设矩阵D如下：

元素d_ij表示射线穿过第i行第j列的像素格的距离。该像素方格吸收率矩阵x(i,j)，则该像素格的吸量为d_ijx_ij。利用上述的CT成像模型，则某射线对扫描区域的所有像素格的吸收量为：

该射线在扫描区域的全部吸收量b(采集量)为所有像素格上吸收量的和，即：

显然，b也就是矩阵D和X的点乘所得到的矩阵的F₁范数。

二、CT成像的数学模型

CT机共有k＝1,2,…,512共512条射线l，且从角度α＝1,2,…,180度，共180个角度实现对扫描区域的扫描。记b(α,k)表示以角度α，第k条射线l(α,k)穿过扫描区域的吸收量。该吸收量在考虑了CT机的吸收衰减后，与公式(4)计算的结果一致。

根据CT机的扫描结果得到b(α,k)，以及CT机的默认参数——角度为α和k条射线，获取射线穿过扫描区域像素格的距离矩阵为D(α,k)，则得到数学模型：

D(a，k)⊙X＝Ab(α，k)；

α＝1,2,…，180，k＝1，2，…512； (5)

其中，矩阵D(α,k)表示α角度，第k条射线l(α,k)穿过扫描区域距离矩阵，X表示扫描区域的吸收率矩阵(即求解量)，b(α,k)表示α角度，第k条射线l(α,k)穿过扫描区域的吸收量(即采集量)，λ表示吸收率的衰减系数，⊙表示矩阵点乘运算。

公式(5)本质上是一个线性方程组，共包含180×512＝92160个方程，n²个未知数。求解该线性方程组即可求得每个像素格的吸收率x_ij，从而能够进行图像的构建。

三、压缩感知的数学模型

压缩感知，又称压缩采样。作为一个新的采样理论，它通过开发信号的稀疏特性，在远小于奈奎斯特(Nyquist)采样率的条件下，用随机采样获取信号的离散样本，然后通过非线性重建算法完美的重建信号。压缩感知理论的核心思想主要包括两点。第一个是信号的稀疏结构。另外一点是不相关特性。稀疏信号的有用信息的获取可以通过一个非自适应的采样方法将信号压缩成较小的样本数据来完成。

压缩感知方法抛弃了当前信号采样中的冗余信息。它直接从连续时间信号变换得到压缩样本，然后在数字信号处理中采用优化方法处理压缩样本。这里恢复信号所需的优化算法常常是一个已知信号稀疏的欠定线性逆问题。

设信号x∈Rⁿ是r-稀疏的，即x的n个分量中最多有r个是非零的，这里r＜＜n。测量矩阵A_m×n是按标准正态分布N(0,1随机产生的矩阵。由测量矩阵A对数据x采样，可得到m(m＜n)维的测量数据b，即：

b_m×l＝A_m×nx_n×1；

通过欠采样数据b来恢复信号x是一个求解不确定线性方程组的病态问题，其结果通常有无穷多解。也就是说，根据观测信号b，是无法确定原始信号x的。但根据Cand`es，Donoho和陶哲轩创立的压缩感知理论，如果信号x∈Cⁿ本身是稀疏信号，则可以通过求解如下的约束优化问题来重构信号[1,2]。

s.t.Ax＝b，x∈Cⁿ；

其中，||x||₀表示向量x的l₀范数，即向量x的非零分量的数量。

但是公式(6)存在难以优化的问题，随着信号长度n的增加，算法的复杂度呈几何级数增加。解决该问题的关键结果是压缩感知理论的重要结论：即如果信号x是稀疏的，则优化问题(6)的稀疏解x可以以概率1通过求解一个凸优化问题得到[1,2]，即：

s·t.Ax＝b，x∈Cⁿ。

在建立模型之初，为了尽量减小因缺少自变量而出现的模型偏差，通常会选择尽可能多的自变量。然而，建模过程需要寻找对因变量最具有解释力的自变量集合，就需要通过自变量的选择来提高模型的解释性和预测精度。指标选择在统计建模过程中是极其重要的问题，Lasso(The Least Absolute Shrinkage and Selectionator operator)模型则是一种能够有效实现指标集合精简的估计方法。

Lasso模型的基本思想是在回归系数的绝对值之和小于一个常数的约束条件下，使残差平方和最小化，从而能够产生某些严格等于0的回归系数，得到解释力较强的指标集合。在应用中，常常遇到特征数量大于训练样本数数量的工程问题，目标是找到一个简单的模型来拟合数据。A＝D向量，标准Lasso模型为

其中，C＞0是一个正则化的回归参数。

四、基于压缩感知的CT逆成像模型

模型(5)如下所示：

D(α，k)⊙X＝λb(α，k)；

α＝1，2，…，180，k＝1，2，…，512；

将矩阵D(α，k)按列连接转换成行向量d(α，k)，即：

d(α，k)＝[d₁₁，d₂₁，…，d_n1，d₁₂，d₂₂，…，dn₂，…，d_1n，d_2n，…，d_nn]；

按照α，k的顺序，将按照上面公式得到行向量组成矩阵E，即：

将矩阵X按列转换成向量x，将b(α，k)按照α，k的顺序，排列成向量b，则模型(5)等价于线性方程组Dx＝b。

通常情况下，扫描区域里大多数像素格不含有吸收组织，所以模型(5)可转换成如下的Lasso模型求解：

其中，δ＞0为正则化参数，E∈R^p×n，x∈Rⁿ，b∈R^p。

五、求解压缩感知的数学模型

交替方向乘子法(Alternating Direction Method of Multipliers，ADMM)是一种适用于求解可分凸优化问题的计算方法，可进行分布式计算，特别适用于解决统计学习，压缩感知，机器学习等理论中的优化问题。交替方向乘子法通过分解协调过程，可以将大的全局优化问题分解为多个较小、较容易求解的子问题，并通过协调子问题的解而得到全局问题的最优解。

交替方向乘子法主要用于如下形式的求解线性等式约束的可分凸优化问题，其形式如下：

其中A∈R^p×n和B∈R^p×m是已知的矩阵，c∈R^p是已知的向量，函数f：Rⁿ→R∪{+∞}和g：R^m→R∪{+∞}是正常的(proper)、闭的凸函数，

和

是闭凸集。

问题(10)中目标函数的拉格朗日函数和增广拉格朗日函数分别为：

L₀(x，y，z)＝f(x)+g(y)+z^T(Ax+By-c) (11)

和

其中，x，y是原始变量，z是对偶变量，ρ＞0是惩罚参数。

求解问题(10)的交替方向乘子法的迭代公式如下：

CT逆成像模型(9)可转化为

s.t.x＝y；

x，y∈Rⁿ.

优化模型(13)的增广拉格朗日函数为

求解优化模型问题(13)的交替方向乘子法的迭代公为

上述交替方向乘子法具有解析解迭代公式：

其中，ρ＞0是惩罚参数，

五、求解过程数值仿真

本节给出基于随机数据和实际脑电信号的两个仿真实验。所有数值仿真实验的结果都是基于如下的软硬件配置：Intel(R)Core(TM)i7-8650U CPU频率1.90G Hz，32GB物理内存，MATLAB版本为R2017b，操作系统为Window 10，64位操作系统。

仿真1：

对64×64，角度压缩比4，运用现有的MATLAB自带的CVX优化，运行时间60.61秒。

仿真2：

对64×64，角度压缩比4，运用本发明的ADMM优化，运行时间32.37秒。

从以上仿真中可明确采用ADMM算法具有优越性。

以上详细描述了本发明的较佳具体实施例。应当理解，本领域的普通技术人员无需创造性劳动就可以根据本发明的构思作出诸多修改和变化。因此，凡本技术领域中技术人员依本发明的构思在现有技术的基础上通过逻辑分析、推理或者有限的实验可以得到的技术方案，皆应在由权利要求书所确定的保护范围内。

Claims (4)

1.一种基于压缩感知的CT成像方法，其特征在于，具体包括以下步骤：

S1、CT机采用角度压缩比进行物体扫描，获取投影数据；

S2、采集投影数据输入至CT逆成像模型，该CT逆成像模型表达式如下：

式中，E为第二矩阵，x为求解量，b为投影数据中的射线吸收量，δ为预设参数；

S3、基于ADMM算法求解CT逆成像模型；

S4、根据求解结果构建二维图像；

所述的第二矩阵E获取步骤如下：

A1、根据CT机参数获取第一矩阵D(α,k)，其中，α为CT机的扫描角度，k为CT机发出射线的数量，第一矩阵D(α,k)表示α角度,第k条射线穿过扫描区域距离矩阵；

A2、将矩阵D(α,k)按列连接转换成行向量d(α,k)；

A3、按照α和k的顺序，将行向量d(α,k)重新组合成第二矩阵E，其表达式为：

步骤S3中，基于ADMM算法求解CT逆成像模型步骤为：

首先，将CT逆成像模型表达式转化为：

s.t.x＝y；

x,y∈Rⁿ；

然后采用如下迭代公式进行求解：

式中，ρ为惩罚参数，τ为正则化变量。

2.根据权利要求1所述的基于压缩感知的CT成像方法，其特征在于，步骤S1中，角度压缩比大于等于4。

3.一种CT成像装置，其特征在于，所述的装置包括处理器以及存储器，所述处理器调用存储器中的程序，用于实现以下步骤：

S1、CT机采用角度压缩比进行物体扫描，获取投影数据；

S2、采集投影数据输入至CT逆成像模型，该CT逆成像模型表达式如下：

式中，E为第二矩阵，x为求解量，b为投影数据中的射线吸收量，δ为预设参数；

S3、基于ADMM算法求解CT逆成像模型；

S4、根据求解结果构建二维图像；

所述的第二矩阵E获取步骤如下：

A1、根据CT机参数获取第一矩阵D(α,k)，其中，α为CT机的扫描角度，k为CT机发出射线的数量，第一矩阵D(α,k)表示α角度,第k条射线穿过扫描区域距离矩阵；

A2、将矩阵D(α,k)按列连接转换成行向量d(α,k)；

A3、按照α和k的顺序，将行向量d(α,k)重新组合成第二矩阵E，其表达式为：

步骤S3中，基于ADMM算法求解CT逆成像模型步骤为：

首先，将CT逆成像模型表达式转化为：

s.t.x＝y；

x,y∈Rⁿ；

然后采用如下迭代公式进行求解：

式中，ρ为惩罚参数，τ为正则化变量。

4.根据权利要求3所述的CT成像装置，其特征在于，步骤S1中，角度压缩比大于等于4。

Images