CN106725347A - 基于随机变量交替方向乘子法荧光分子断层成像重建方法 - Google Patents

Abstract

本发明公开了一种基于随机变量交替方向乘子法的荧光分子断层成像重建方法，利用随机对偶坐标下降方法的随机性和可分解性将交替方向乘子方法转换成随机参数形式，然后通过随机参数的交替方向乘子方法的参数交替更新来求解，其实现步骤如下：(1)大规模荧光数据采集；(2)建立测量数据与目标分布之间的线性关系；(3)转换线性关系为凸优化问题；(4)使用随机参数的对偶坐标下降方法与交替方向乘子法来交替求解凸优化问题，获得目标分布图和重建时间。在保证重建图像质量的同时，有效提高了荧光分子断层成像的重建效率，在医学分子影像、重建方法等领域有重要的应用价值。

Description

基于随机变量交替方向乘子法荧光分子断层成像重建方法

技术领域

本发明属于分子影像领域，涉及一种基于随机变量交替方向乘子法的荧光分子断层成像重建方法。

背景技术

荧光分子断层成像(以下简称FMT)是近年发展起来的一种新型成像模态。它利用外部光源激发荧光探针(荧光团、荧光染料等)使其发射光子，利用荧光采集装置收集荧光信号，结合数学模型，利用反演算法可在体重建荧光团的位置和浓度分布。通常采用电荷耦合器件(charge coupled device，CCD)全角度非接触式成像系统来采集老鼠荧光投影数据，FMT重建是个高度病态逆问题，大规模的多投影荧光数据能够降低FMT的病态性。但是在使用大规模荧光投影数据进行重建时会占用大量的内存，并且花费大量时间。

发明内容

本发明的目的是克服上述FMT重建时的病态性以及提高重建效率，提出了一种基于随机变量交替方向乘子法的重建方法。本发明采用多点激发和多角度测量，获得了大规模的测量数据，在重建过程中，在光传输模型和Robin边界条件下，结合有限元方法将表面测量所得的光子能量信息以及所需重建的荧光目标，用一个系统矩阵建立线性关系，由于在荧光分子断层成像中，荧光分子探针在生物组织中的分布相对稀疏，根据压缩感知理论，可以将该线性关系转换为带惩罚项的凸优化问题，采用随机变量交替方向乘子方法来求解，从而获得重建的荧光目标三维分布和浓度。

为实现上述目的，本发明的具体步骤如下：

(1)获得测量数据

主要包括：

1.利用激发光源对固定在旋转台上的重建目标进行360度的透射式断层扫描；

2.使用光学检测仪器获得测量数据，同时获取表面测量得到的荧光信息Φ。

(2)获得重建目标的结构信息以及光学参数信息

(3)在扩散近似模型和Robin边界条件下，结合有限元方法，以重建目标的结构信息和光学参数为先验信息，用系统矩阵A建立表面测量得到的荧光信息Φ和所要重建的荧光目标X的线性关系。

(4)根据压缩感知理论，将上述线性关系转化为“损失函数+正则化项”的凸优化问题：

其中：z₁,z₂,...,z_n∈R^p，B∈R^p×d，w为系统矩阵，线性约束条件：Zx+By＝0。

(5)对于步骤(4)中的凸优化问题，采用随机变量交替方向乘子方法来进行迭代求解，该方法结合了随机参数的对偶坐标下降方法的可分解性与交替方向乘子法。其中，随机对偶坐标上升方法将系统矩阵分成K个子矩阵，每次迭代只选择其中一个子矩阵，使每次迭代只在随机选择的子矩阵中进行，利用随机对偶坐标下降法的随机性和可分解性结合到交替方向乘子法中，同时引入mini-batch方法加速收敛，使子矩阵转换成mini-batch I_k(k∈{1,...,K})，进行迭代计算，本方法迭代形式为：

y^(t)←q^(t)-prox(q^(t)|nψ(ρηB·)/(ρηB))

w^(t)←w^(t-1)-γρ{n(Zx^(t)+By^(t))-(n-n/K)(Zx^(t-1)+By^(t-1))}

(6)通过以上方法，x和y的值交替更新，以此来对w进行迭代求解，设置迭代次数T，令t＝1,…T，迭代结束后输出结果w^(T)，完成重建过程，最后输出重建结果。

(7)设置对比方法，对于步骤(4)中的凸优化问题，采用线性交替方向乘子方法进行求解，线性交替方向乘子方法迭代形式为：

w^(t)←w^(t-1)-γρ{Zx^(t)+By^(t)}

设置迭代次数T，令t＝1,…T，迭代结束后输出结果w^(T)，完成对比方法的重建过程，最后输出重建结果。

8)显示结果，将重建结果和成像目标的解剖结构进行图像融合，在Tecplot软件中进行展示；并对比本发明的方法和对比方法的重建时间和重建精度。

采用上述方法后，本发明具有以下优点：

第一，本发明采用多点激发，获取了大规模的多投影荧光数据，能有效降低FMT重建问题病态性，提高重建图像的精度。

第二，本发明将重建问题转化为“损失函数+正则化项”的凸优化问题，针对大规模的荧光数据所产生复杂系统矩阵，采用随机参数交替方向乘子方法分解系统矩阵为多个子矩阵，使得求解过程不需要占用大量的内存，实现加速重建。

第三，本发明在采用随机参数交替方向乘子方法的同时，又采用了mini-batch方法，进一步加速收敛，减少重建所需要花费的时间

附图说明

图1为本发明的荧光分子断层成像重建方法的流程图；

图2为本发明中用于仿真实验的数字鼠模型及激发点分布；

图3为本发明的重建算法获得的重建结果及荧光分布图；

其中图(a)和图(c)分别显示为本发明的重建结果和在Z＝16.4mm截面的荧光目标荧光产额值分布图；图(b)和图(d)为线性交替方向乘子方法的重建结果在Z＝16.4mm截面的荧光目标荧光产额值分布图；图(e)为两种方法重建结果的定量分析。

具体实施方式

下面结合附图具体实施例来详细描述本发明的技术方案。

结合附图1对本发明的步骤做进一步的描述：

(1)大规模荧光数据采集

1.利用激发光源对固定在旋转台上的重建目标进行360度的透射式断层扫描；

透射式断层成像，将激光器与光学检测仪器放置在成像目标的两侧，激光照射重建目标激发荧光团发出荧光，荧光穿透成像目标被激光器对面的光学检测仪器检测到。

多角度透射式断层扫描，用电脑控制旋转台等间隔旋转一定角度，本发明中为10度，激光器发射点状激光照射成像目标，转一个角度激发一次，这样就进行了多次激发，从而实现了多角度的透射式成像。

2.使用光学检测仪器获得测量数据，获得荧光信息Φ；

在步骤1中，激光器照射一次成像目标，光学检测仪器就采集一组荧光信号，得的一组测量数据，多角度激发对应产生多组测量数据，将数据应用非接触式光学断层成像方法中描述的生物体表面三维能量重建技术获取成像目标体表面的三维荧光数据分布。

(2)基于扩散近似模型和Robin边界条件，结合有限元方法，将重建目标的结构信息和光学参数信息作为先验信息，用系统矩阵A建立表面测量得到的荧光信息Φ和所要重建的荧光目标X的线性关系。

AX＝Φ

其中，X表示需要重建的荧光目标；Φ表示表面测量得到的荧光信息；A为一个大小为m*n的系统矩阵。系统矩阵A包含了在前向问题求解过程中所得到的每个激发点探测到的结点的荧光强度。

(3)FMT系统中，荧光分子探针在生物组织中的分布相对较为稀疏，因此根据压缩感知理论，可建立“损失函数+正则化项”的凸优化模型对上述线性关系进行转化求解。

线性约束条件：Zx+By＝0。

(4)对于步骤(3)中的凸优化问题，本发明采用随机变量交替方向乘子方法来进行迭代求解，该方法结合了随机参数的对偶坐标上升方法的可分解性与交替方向乘子法。包括：

1.随机对偶坐标下降方法，利用此方法来解决凸优化问题的对偶问题：

针对复杂系统矩阵A，首先将系统矩阵分解成K个子空间，因此目标函数也被分解成了K个子问题，以均匀概率1/K选取一个子空间进行迭代。

2.交替方向乘子方法，利用交替方向乘子法对步骤(4)中的凸优化问题进行迭代求解，其迭代形式为

w^(t)←w^(t-1)-γρ{Zx^(t)+By^(t)}

将随机对偶坐标下降方法的可分解性代入交替方向乘子法，形成本发明的随机参数的交替方向乘子方法，并引入mini-batch方法加速收敛，将子矩阵转换成mini-batchI_k(k∈{1,...,K})，每个mini-batch被选中概率统一为1/K，使用mini-batch方法进行迭代时，参数更新规则如下：

w^(t)←w^(t-1)-γρ{n(Zx^(t)+By^(t))-(n-n/K)(Zx^(t-1)+By^(t-1))}

其中，代入到交替方向乘子方法中，迭代形式转换成：

y^(t)←q^(t)-prox(q^(t)|nψ(ρηB·)/(ρηB))

w^(t)←w^(t-1)-γρ{n(Zx^(t)+By^(t))-(n-n/K)(Zx^(t-1)+By^(t-1))}

(5)通过以上方法，x和y的值交替更新，以此来对w进行迭代求解，设置迭代次数T，令t＝1,…T，迭代结束后输出结果w^(T)，完成重建过程，最后输出重建结果，得到重建的荧光目标X和重建时间。

(6)设置对比方法，对于步骤(4)中的凸优化问题，采用步骤(4)中第2点的线性交替方向乘子方法进行迭代求解。设置迭代次数T，令t＝1,…T，迭代结束后输出结果w^(T)，完成重建过程，最后输出重建结果，得到重建的荧光目标X和重建时间。

(7)结果展示及评价，将重建结果和成像目标的解剖结构进行图像融合，在Tecplot软件中进行展示；同时引入位置误差L作为重建结果的评价标准以及重建时间作为本发明采用的方法的评价标准。位置误差即为重建目标中心和真实目标之间的欧氏距离。位置误差越小，说明重建结果与真实位置越接近。

附图2(a)为用于仿真实验的数字鼠躯干三维模型图，包括的器官有：心脏，肺，肝脏，胃，肾脏，荧光目标，图2(b)为激发点分布图。

附图3为基于本发明方法的重建结果。图(a)展示为本发明方法的重建结果，图(b)为线性交替方向乘子方法的重建结果，图(c)为本发明重建在z＝16.4mm处的截面，图(d)为线性交替方向乘子法重建在z＝16.4mm处的截面，图(e)为两种方法的重建的结果分析。

在一种实施例中，荧光目标真实位置中心点坐标为(11.9mm，6.40mm，16.40mm)，重建结果中心点坐标为(11.83mm，6.30mm，15.96mm)。位置误差L＝0.46mm，本发明的重建方法的重建时间为28.4s。对比实验的重建结果中心坐标为(11.64mm，6.36mm，16.12mm)，位置误差L＝0.38，重建时间为147s，证明，基于本发明的FMT重建，位置误差差别不大，但重建时间是常用的线性交替方向乘子方法的1/5左右，是一种有效的荧光分子断层成像重建方法。

以上所述，仅为本发明中的一个具体实例，但本发明的保护范围并不局限于此，任何熟悉该技术的人在本发明所揭露的技术范围内，可理解想到的变换或替换，都应涵盖在本发明的包含范围之内，因此，本发明的保护范围应该以权利要求书的保护范围为准。

Claims (1)

1.一种基于随机变量交替方向乘子法的重建方法，其特征在于：包括以下步骤：

(1)多点激发和多角度测量获得大规模荧光数据：

①利用激发光源对固定在旋转台上的重建目标进行360度的透射式断层扫描；

透射式断层成像，将激光器与光学检测仪器放置在成像目标的两侧，激光照射重建目标激发荧光团发出荧光，荧光穿透成像目标被激光器对面的光学检测仪器检测到。

多角度透射式断层扫描，用电脑控制旋转台等间隔旋转一定角度，本发明中为10度，激光器发射点状激光照射成像目标，转一个角度激发一次，这样就进行了多次激发，从而实现了多角度的透射式成像。

②使用光学检测仪器获得测量数据，获得荧光信息Φ：在步骤①中，激光器照射一次成像目标，光学检测仪器就采集一组荧光信号，得的一组测量数据，多角度激发对应产生多组测量数据，将数据应用非接触式光学断层成像方法中描述的生物体表面三维能量重建技术获取成像目标体表面的三维荧光数据分布。

(2)获得重建目标的结构信息以及光学参数信息：基于扩散近似模型和Robin边界条件，结合有限元方法，将重建目标的结构信息和光学参数信息作为先验信息，用系统矩阵A建立表面测量得到的荧光信息Φ和所要重建的荧光目标X的线性关系：

AX＝Φ

其中，X表示需要重建的荧光目标；Φ表示表面测量得到的荧光信息；A为一个大小为m*n的系统矩阵。系统矩阵A包含了在前向问题求解过程中所得到的每个激发点探测到的结点的荧光强度。

(3)根据压缩感知理论，将上述线性关系转化为“损失函数+正则化项”的凸优化问题：FMT系统中，荧光分子探针在生物组织中的分布相对较为稀疏，因此根据压缩感知理论，可建立“损失函数+正则化项”的凸优化模型对上述线性关系进行转化求解：

$\underset{w\∈R^d}{min}\frac{1}{n}\underset{i=1}{\overset{n}{\Σ}}{f}_i\left({z_i}^{T}w\right)+\mathrm{\ψ}\left(B^{T}w\right)$

线性约束条件：Zx+By＝0。

(4)对于步骤(3)中的凸优化问题，采用随机变量交替方向乘子方法来进行迭代求解，该方法结合了随机参数的对偶坐标上升方法的可分解性与交替方向乘子法。包括：

I、随机对偶坐标下降方法，利用此方法来解决凸优化问题的对偶问题：

$\underset{y\∈R^d}{\underset{x\∈R^n}{min}}\{-\frac{1}{n}\underset{i=1}{\overset{n}{\Σ}}{f}_i^{*}\left(x_i\right)+{\mathrm{\ψ}}^{*}\left(\frac{y}{n}\right)|Zx+By=0\}$

针对复杂系统矩阵A，首先将系统矩阵分解成K个子空间，因此目标函数也被分解成了K个子问题，以均匀概率1/K选取一个子空间进行迭代。

II、交替方向乘子方法，利用交替方向乘子法对步骤(4)中的凸优化问题进行迭代求解，其迭代形式为

$y^{\left(t\right)}\&LeftArrow;\underset{y}{\arg min}\{{\mathrm{n\&psi;}}^{*}\left(\frac{y}{n}\right)-<w^{(t-1)},{\mathrm{Zx}}^{(t-1)}+By>+\frac{\mathrm{\&rho;}}{2}||{\mathrm{Zx}}^{(t-1)}+By|{|}^2\}$

$x^{\left(t\right)}\&LeftArrow;\underset{x}{\arg min}\{{\&Sigma;}_{i=1}^n{f}_i^{*}\left(x_i\right)-<w^{(t-1)},Zx+{\mathrm{Bt}}^{\left(t\right)}>+\frac{\mathrm{\&rho;}}{2}||Zx+{\mathrm{By}}^{\left(t\right)}|{|}^2\}$

w^(t)←w^(t-1)-γρ{Zx^(t)+By^(t)}

将随机对偶坐标下降方法的可分解性代入交替方向乘子法，形成本发明的随机参数的交替方向乘子方法，并引入mini-batch方法加速收敛，将子矩阵转换成mini-batch I_k(k∈{1,...,K})，每个mini-batch被选中概率统一为1/K，使用mini-batch方法进行迭代时，参数更新规则如下：

$q^{\left(t\right)}=y^{(t-1)}+\frac{B^T}{{\mathrm{\&rho;\&eta;}}_B}\{w^{(t-1)}-\mathrm{\&rho;}({\mathrm{Zx}}^{(t-1)}+{\mathrm{By}}^{(t-1)})\}$

$y^{\left(t\right)}\&LeftArrow;\underset{y}{\arg min}\{{\mathrm{n\&psi;}}^{*}\left(\frac{y}{n}\right)-<w^{(t-1)},{\mathrm{Zx}}^{(t-1)}+By>+\frac{\mathrm{\&rho;}}{2}||{\mathrm{Zx}}^{(t-1)}+By|{|}^2\}$

w^(t)←w^(t-1)-γρ{n(Zx^(t)+By^(t))-(n-n/K)(Zx^(t-1)+By^(t-1))}

其中，代入到交替方向乘子方法中，迭代形式转换成：

$q^{\left(t\right)}=y^{(t-1)}+\frac{B^T}{{\mathrm{\&rho;\&eta;}}_B}\{w^{(t-1)}-\mathrm{\&rho;}({\mathrm{Zx}}^{(t-1)}+{\mathrm{By}}^{(t-1)})\}$

y^(t)←q^(t)-prox(q^(t)|nψ(ρηB·)/(ρηB))

$x_I^{\left(t\right)}\&LeftArrow;prox(x_I^{(t-1)}+\frac{Z_I^T}{{\mathrm{\&rho;\&eta;}}_{Z,I}}\{w^{(t-1)}-\mathrm{\&rho;}({\mathrm{Zx}}^{(t-1)}+{\mathrm{By}}^{\left(t\right)})\left\}\right|\frac{{\mathrm{\&Sigma;}}_{i\&Element;I}{f}_i^{*}}{{\mathrm{\&rho;\&eta;}}_{Z,I}})$

w^(t)←w^(t-1)-γρ{n(Zx^(t)+By^(t))-(n-n/K)(Zx^(t-1)+By^(t-1))}

通过以上方法，x和y的值交替更新，以此来对w进行迭代求解，设置迭代次数T，令t＝1,…,T，迭代结束后输出结果w^(T)，完成重建过程，最后输出重建结果，得到重建的荧光目标X和重建时间。

(5)设置对比方法，对于步骤(4)中的凸优化问题，采用步骤(4)中II的线性交替方向乘子方法进行迭代求解。设置迭代次数T，令t＝1,…,T，迭代结束后输出结果w^(T)，完成重建过程，最后输出重建结果，得到重建的荧光目标X和重建时间。

(6)结果展示及评价：将重建结果和成像目标的解剖结构进行图像融合，在Tecplot软件中进行展示；同时引入位置误差L作为重建结果的评价标准以及重建时间作为评价标准，其中位置误差即为重建目标中心和真实目标之间的欧氏距离。