一种基于伪极坐标TV最小化直线轨迹CT图像重建方法

摘要

本发明公开了一种基于伪极坐标TV最小化直线轨迹CT图像重建方法，克服了现有技术中，直线轨迹计算机断层成像（linearcomputedtomography,LCT）技术的有限角度图像重建的问题。该发明包含以下步骤——步骤1：建立TV最小化重建模型；步骤2：利用ADM最小化TV模型；步骤3：利用PPFFT实现图像空-频域变换；步骤4：实现并运行算法，获得重建图像。该LCT重建技术基于交替方向法设计了TV最小化模型的求解算法，具有稳定的收敛性；并且，由于采用了伪极快速傅里叶变换，该算法具有优异的重建精度和计算效率。基于伪极坐标TV最小化LCT图像重建技术，在LCT技术投入实用化中具有重要意义。

说明书

技术领域

该发明涉及一种CT图像重建方法，特别是涉及一种基于伪极坐标TV最小化直线轨迹CT图像重建方 法。

背景技术

传统CT技术一般采用旋转轨迹(如圆轨迹，螺旋轨迹等)获取不同角度下的投影数据，需要对物体 或者扫描系统做旋转操作。在一些应用场合，如工业在线检测、海关集装箱检查、机场、车站的行李安全 检查等，这些场合对扫描速率要求较高，基于旋转轨迹的CT适用性并不佳。目前最新发展起来的一种采 用大张角射线源、大面积阵列探测器和直线扫描轨迹技术的新型X射线断层成像系统——直线轨迹CT (linear computed tomography,LCT)受到了广泛关注。在LCT扫描过程中，探测器和射线源保持相对 静止，待检测物体相对于探测器和射线源作直线运动，通过采集不同位置上的投影来重建物体的断层图像。 与传统CT成像方式相比，直线轨迹CT能够实现快速检测，在对扫描速度要求较高的应用场合具有重要的 实用价值。

然而，受射线源张角和探测器阵列尺寸限制，LCT投影数据采集通常限定在一个有限的范围内(小于 180°)，是一个典型的有限角度问题。有限角度扫描无法提供满足经典解析型重建算法所要求的投影数据 集，因此解析重建算法难以获得高质量的重建图像。图像重建质量是整个LCT系统投入实际应用的关键指 标，因此LCT图像重建技术成为了近年来CT成像领域研究的热点问题。基于压缩感知理论和正则化方法 发展起来的重建算法，能够对较低采样或有限角度下投影数据进行有效重建。基于压缩感知理论发展起来 的TV最小化重建算法基于图像内在的梯度稀疏特性，并利用该特性提升了有限角度重建的精度。另一方 面，LCT的投影频域采样具有伪极坐标分布的特性；该特性表明，基于伪极坐标的快速傅里叶变换(pseudo  polar fast Fourier transform,PPFFT)能够有效避免传统频域重建中的傅里叶域的插值问题，可以获 得更高的频域重建精度。因此，怎样结合TV模型和LCT投影采样频域分布特性设计高效高精度、具有稳 定收敛性的图像重建算法对于LCT技术的发展关键问题，并且具有重要意义。

针对直线轨迹计算机断层成像(linear computed tomography,LCT)技术的有限角度图像重建问题， 根据LCT投影频域采样分布特性，提出基于伪极坐标总变分(total variation,TV)最小化图像重建技 术。该技术分为投影数据预处理，建立重建模型，设计求解算法三部分：首先对采集到的投影数据进行频 域变换得到频域观测数据，然后基于图像梯度稀疏特性建立TV最小化重建优化模型，最后采用交替方向 法和伪极坐标快速傅里叶变换(PPFFT)设计求解算法。

对于采样投影数据不足的重建问题，通常采用迭代类算法进行处理。传统的迭代算法主要有基于数 据域的变换类算法、代数重建技术以及统计迭代类算法。虽然传统的迭代算法较解析类算法有较好重建质 量，但是对于有限角度的重建质量通常难以满足实用。近年来，基于压缩感知理论的重建算法分析了图像 内在的稀疏特性，并利用该特性提升了有限角度重建的精度。2007年，高河伟等人针对LCT的有限角度问 题结合GPEL Gerchberg-Papoulis-type extrapolation using Linogram)和TV正则化技术提出了GPEL-TV 算法，在抑制噪声的同时具有一定的边缘保持作用。在已有的凸集投影(projection onto convex sets, POCS)算法基础上，潘晓川等人以待重建图像TV最小化作为目标提出了ASD-POCS (adaptive-steepest-descent-projection onto convex sets)算法，针对不完全数据CT图像重建问 题取得了优于传统迭代类算法的效果。

采用交替方向法(alternating direction method,ADM)的优化框架的重建技术在TV最小化模型 的求解中表现出了优异性能。2011年，Vandeghinste等人针对TV最小化重建模型，将split-Bregman优 化算法应用于稀疏角度采样的CT图像重建中，取得了较高精度的重建。2013年，Zhang等人针对LCT有 限角度问题，采用各向异性TV优化模型重建提出了交替方向TV最小化(alternating direction TV  minimization,ADTVM)重建算法。该算法采用交替方向法(alternating direction method,ADM)，将 原优化问题的增广拉格朗日函数拆分为两个具有解析解的子问题，求解算法较为高效稳定，对有限角度问 题的解决起到了一定的推动作用。

LCT投影采样频域分布具有特殊性的性质：所有的采样点分布于以原点为中心的同心矩阵上，并且沿 中心成对称辐射状，每条过中心的片段所在直线呈等斜率间隔分布。LCT采样的频域分布类似于极坐标分 布但又相区别，故称其为伪极坐标分布。该种频域分布表明，采用基于伪极坐标的PPFFT技术能够实现图 像空间和采样频域空间的无插值精确快速变化。受频域重建启发，基于TV最小化重建模型，采用ADM并 结合PPFFT技术是设计高效高精度LCT重建算法的有效方法。

发明内容

本发明克服了现有技术中直线轨迹计算机断层成像(linearcomputedtomography,LCT)技术的有 限角度图像重建的问题，提供一种使用效果较好的基于伪极坐标TV最小化直线轨迹CT图像重建方法。 本发明的技术解决方案是，提供一种基于伪极坐标TV最小化直线轨迹CT图像重建方法，包括以下步骤：

步骤1：建立TV最小化重建模型；

步骤2：利用ADM最小化TV模型；

步骤3：利用PPFFT实现图像空-频域变换；

步骤4：实现并运行算法，获得重建图像。

所述建立TV最小化重建模型包括关于变量l的傅里叶变换：即获得函数f(x,y)的二维傅里叶变换； 每个探元之间的长度间隔为Δt，且都能接收到光源的辐射；设光源有2N个采样位置，其索引为n，采样 步进为Δl；从而ω_x和σ分别离散化表达为于是图像重建可以表达为寻找如下线性方程组 的解：

$F_{p}f=\hat{f},$

，其中，为观测到的傅里叶采样数据，为向量化的离散物体函数；组合系数矩阵 表示对f作离散伪极傅里叶变换。建立LCT的TV最小化重建模型，也就是求解如下带约 束条件的优化问题：

f^*＝arg min||f||_TV；

s.t.F_pf＝F(p),

其中，为沿图像i方向的差分算子，F(p)表示对投影数据做一维傅里 叶变换。

所述的利用ADM最小化TV模型包括：

(1)建立无约束重建模型：采用ADM框架给出其求解算法；令D_if＝z_i，则该TV重建 模型可以转化为：

$\min \underset{i}{\Sigma }{\left|\right|z_i\left|\right|}_1+\frac{\lambda }{2}{\left|\right|F_{p}f-\hat{f}\left|\right|}_2^2,$

s.t.D_if＝z_i,

其中为保真项因子，用来调整目标函数中观测数据的一致性程度；式为带约束优化模型，为了将其 转化为无约束优化问题，利用增广拉格朗日函数进行如下转化：

$\underset{f,z_i}{\min }{L}_A(f,z_i)=\frac{\lambda }{2}{\left|\right|F_{p}f-\hat{f}\left|\right|}_2^2+\underset{i}{\Sigma }(u_i^T(D_{i}f-z_i)+\frac{{\rho }_i}{2}{\left|\right|D_{i}f-z_i\left|\right|}_2^2+{\left|\right|z_i\left|\right|}_1),$

其中标量二次项惩罚系数，为乘子；

(2)利用ADM求解模型：基于ADM框架采用分离变量的方法将转化为两个单目标优化问题，也即

$\left(\begin{array}{c}\underset{f}{\min }{L}_{z_i}\left(f\right)=\frac{\lambda }{2}{\left|\right|F_{p}f-\hat{f}\left|\right|}_2^2+\underset{i}{\Sigma }(u_i^T(D_{i}f-z_i)+\frac{{\rho }_i}{2}{\left|\right|D_{i}f-z_i\left|\right|}_2^2)\\ \underset{z_i}{\min }{L}_f\left(z_i\right)=u_i^T(D_{i}f-z_i)+\frac{{\rho }_i}{2}{\left|\right|D_{i}f-z_i\left|\right|}_2^2+{\left|\right|z_i\left|\right|}_1\end{array}\right).$

采用ADM分别对子问题求取其最小值点，于是，基于ADM的求解迭代公式为：

$\left(\begin{array}{c}{f}^{(k+1)}={(\lambda F_p^H{F}_p+\underset{i}{\Sigma }{\rho }_i{D}_i^T{D}_i)}^{+}(\lambda F_p^H\hat{f}-\underset{i}{\Sigma }{D}_i^T(u_i^T-{\rho }_i{z}_i^{\left(k\right)}))\\ z_i^{(k+1)}=\max \left\{\right|D_i{f}^{\left(k\right)}+\frac{u_i^{\left(k\right)}}{{\rho }_i}|-\frac{1}{{\rho }_i},0\}·\mathrm{sgn}(D_i{f}^{\left(k\right)}+\frac{u_i^{\left(k\right)}}{{\rho }_i})\\ u_i^{(k+1)}=u_i^{\left(k\right)}+{\rho }_i(D_i{f}^{(k+1)}-z_i^{(k+1)})\end{array}\right)..$

所述的利用PPFFT实现图像空-频域变换包括：

(1)伪极傅里叶变换定义：LCT重建中所应用的F_pf及定义为：

$\left(\begin{array}{c}{F}_{p}f\stackrel{\Delta }{=}\frac{\Delta l^2}{2\pi }\underset{k_1=-N}{\overset{N-1}{\Sigma }}\underset{k_2=-N}{\overset{N-1}{\Sigma }}f(k_1\mathrm{\Delta l},k_2\mathrm{\Delta l})e^{-j(\frac{2\mathrm{\pi n}}{2N}{k}_1+\frac{2\mathrm{\pi m}}{2N}{k}_2·\frac{\mathrm{n\Delta t}}{D})}\\ =\hat{f}(\frac{\mathrm{n\pi }}{\mathrm{\Delta lN}},\frac{\mathrm{m\pi }}{\mathrm{\Delta lN}}·\frac{\mathrm{n\Delta t}}{D}),\end{array}\right)$

及

$\left(\begin{array}{c}{F}_p^H\hat{f}\stackrel{\Delta }{=}\frac{\mathrm{\Delta t\Delta l}}{4{\left(\mathrm{N\Delta l}\right)}^2}\underset{n=-N}{\overset{N-1}{\Sigma }}\underset{m=-M}{\overset{M-1}{\Sigma }}\hat{f}(\frac{\mathrm{n\pi }}{\mathrm{\Delta lN}},\frac{\mathrm{m\pi }}{\mathrm{\Delta lN}}·\frac{\mathrm{n\Delta t}}{D})e^{j(\frac{2\mathrm{\pi n}}{2N}{k}_1+\frac{2\mathrm{\pi m}}{2N}{k}_2·\frac{\mathrm{n\Delta t}}{D})}\\ =\tilde{f}(k_1\mathrm{\Delta l},k_2\mathrm{\Delta l}).\end{array}\right)$

现对F_pf分析并推导快速计算方法，F_pf实际上可简写为：

$\hat{f}[n,m]\stackrel{\Delta }{=}\frac{1}{2\pi }\underset{k_1=-N}{\overset{N-1}{\Sigma }}\underset{k_2=-N}{\overset{N-1}{\Sigma }}f(k_1,k_2)\exp (-j\frac{2\mathrm{\pi n}{k}_1}{2N}-j\frac{2\mathrm{\pi m}{k}_2}{2N}·\mathrm{n\alpha }),$

其中为方便计算且不失一般性，已令Δt＝1,Δl＝1；

(2)对阵列f(k₁,k₂)的每列做2N点一维快速傅里叶变换：展开二重求和式，可以得到：

$\left(\begin{array}{c}\hat{f}[n,m]=\frac{1}{2\pi }\underset{k_2=-N}{\overset{N-1}{\Sigma }}\exp (-j\frac{2\pi }{2N}·\mathrm{m\alpha }·n{k}_2)\underset{k_1=-N}{\overset{N-1}{\Sigma }}f(k_1,k_2)\exp (-j\frac{2\pi }{2N}n{k}_1)\\ =C·\underset{k_2=-N}{\overset{N-1}{\Sigma }}\exp (-j\frac{2\pi }{2N}·\mathrm{m\alpha }·n{k}_2)\underset{k_1=-N}{\overset{N-1}{\Sigma }}f(k_1,k_2)\exp (-j\frac{2\pi }{2N}n(k_1+N)).\end{array}\right)$

其中，为与求和过程无关的常量；分析内层关于变量k₁的求和式：

${\hat{f}}_1[n,k_2]=\underset{k_1=0}{\overset{2N-1}{\Sigma }}f(k_1,k_2)\exp (-j\frac{2\pi }{2N}n{k}_1),$

其计算过程实际为针对阵列f(k₁,k₂)的每列做2N点一维DFT，其计算可使用FFT技术实现快速计算；

(3)对阵列的每行做2N点线性调频Z变换：外层对变量k₂的求和：

$\hat{f}[n,m]=\frac{1}{2\pi }\exp \left(j\frac{2\pi }{2N}\mathrm{nN}\right)\underset{k_2=-N}{\overset{N-1}{\Sigma }}{\hat{f}}_1[n,k_2]\exp (-j\frac{2\mathrm{\pi m}{k}_2}{2N}·\mathrm{n\alpha }).$

此式中若nα＝1，则该式实际为对阵列的每行做2N点一维DFT，此时快速计算仍可以利用FFT 实现；当nα≠1时，式实际上是对序列做线性调频Z变换(chirp-Z transform,CZT)；令nα＝η， 将式简写为：

${\hat{f}}_n\left[m\right]=\underset{k_2=-N}{\overset{N-1}{\Sigma }}{\hat{f}}_{1,n}\left[k_2\right]\exp (-j\frac{2\mathrm{\pi m}{k}_2}{2N}·\eta ).$

注意到$2m{k}_2=m^2+k_2^2-{(m-k_2)}^2,$将该关系式代入式可得：

$\left(\begin{array}{c}{\hat{f}}_n\left[m\right]=\underset{k_2=-N}{\overset{N-1}{\Sigma }}{\hat{f}}_{1,n}\left[k_2\right]\exp (-j\frac{\mathrm{\pi \eta }}{2N}·(m^2+k_2^2-{(m-k_2)}^2))\\ =\exp (-j\frac{\mathrm{\pi \eta }}{2N}{m}^2)\underset{k_2=-N}{\overset{N-1}{\Sigma }}{\hat{f}}_{1,n}\left[k_2\right]\exp (-j\frac{\mathrm{\pi \eta }}{2N}·k_2^2)\exp (j\frac{\mathrm{\pi \eta }}{2N}·{(m-k_2)}^2),\end{array}\right)$

若定义则上式可以进一步改写为：

${\hat{f}}_n\left[m\right]=s\left[m\right]\underset{k_2=-N}{\overset{N-1}{\Sigma }}{\hat{f}}_{1,n}\left[k_2\right]·s\left[k_2\right]·s[m-k_2].$

式表明可由三步操作算得：首先用序列s[k₂]乘以得到然后利用 与s[k₂]计算卷积，最后再乘以s[m]即可得到在的计算过程中，卷积运算占 用了大部分计算时间，利用一维FFT可以替代此过程：分别先对与s[k₂]计算一维FFT，对 应相乘后做一维IFFT即可(为了避免FFT带来的周期循环效应，在做一维FFT之前应先对各序列补零至 二倍长度)；利用FFT可以使原卷积运算O(N²)复杂度降为N log N。

所述的实现并运行算法，获得重建图像包括：

输入数据p,初始化λ,ρ_i＞0.且设定最大迭代次数N，且令k＝0.

(1)数据预处理：

$\hat{f}(\omega ,-\mathrm{\omega t}/D)\leftarrow \frac{D/\sqrt{2\pi }}{\sqrt{D^2+t^2}}F\left(p(l-t,t)\right)$

执行如下迭代过程：

(2)更新f，公式如下：

$f^{(k+1)}\leftarrow {(\lambda F_p^H{F}_p+\underset{i}{\Sigma }{\rho }_i{D}_i^T{D}_i)}^{+}(\lambda F_p^H\hat{f}-\underset{i}{\Sigma }{D}_i^T(u_i^T-{\rho }_i{z}_i^{\left(k\right)}))$

(3)更新z_i，公式如下：

$z_i^{(k+1)}\leftarrow \max \left\{\right|D_i{f}^{\left(k\right)}+\frac{u_i^{\left(k\right)}}{{\rho }_i}|-\frac{1}{{\rho }_i},0\}·\mathrm{sgn}(D_i{f}^{\left(k\right)}+\frac{u_i^{\left(k\right)}}{{\rho }_i})$

(4)更新u_i，公式如下：

$u_i^{(k+1)}\leftarrow u_i^{\left(k\right)}+{\rho }_i(D_i{f}^{(k+1)}-z_i^{(k+1)})$

(5)k←k+1

若“未达到迭代设定的最大次数,”则返回第2步，否则停止；其中，F_p和由PPFFT实现图像空- 频域变换实现快速计算；迭代达到最大迭代轮数退出时，f^(N)即为输出重建图像。

与现有技术相比，本发明具有以下优点：利用伪极傅里叶变换和交替方向法实现了LCT无插值频域 重建高精度和高效重建。该LCT重建技术基于交替方向法设计了TV最小化模型的求解算法，具有稳定的 收敛性；并且，由于采用了伪极快速傅里叶变换，该算法具有优异的重建精度和计算效率。基于伪极坐标 TV最小化LCT图像重建技术，在LCT技术投入实用化中具有重要意义。

附图说明

图1是本发明一种基于伪极坐标TV最小化直线轨迹CT图像重建方法的流程示意图；

图2是本发明一种基于伪极坐标TV最小化直线轨迹CT图像重建方法中LCT扫描示示意图；

图3是本发明一种基于伪极坐标TV最小化直线轨迹CT图像重建方法中LCT投影几何模型示意图；

图4是本发明一种基于伪极坐标TV最小化直线轨迹CT图像重建方法中体模与重建结果的示意图；

图5是本发明一种基于伪极坐标TV最小化直线轨迹CT图像重建方法中重建均方根误差随迭代轮数 的下降曲线图。

具体实施方式

附图说明中标号1是探测器，2是被检测物体。

下面结合附图和具体实施方式对本发明基于伪极坐标TV最小化直线轨迹CT图像重建方法的技术方 案作进一步说明：

步骤1：建立TV最小化重建模型；步骤2：利用ADM最小化TV模型；步骤3：利用PPFFT实现图像空 -频域变换；步骤4：实现并运行算法，获得重建图像。

如图2所示，在LCT扫描过程中，探测器和射线源保持相对静止，待检测物体相对于探测器和射线源 作直线运动。如图3所示，点O为物体坐标系O-xy原点，其在光源运动直线上的投影为O'。为分析方 便引入一个等价虚拟探测器，该探测器与x轴重合，并且其中心为光源S在x轴上的垂足O"。在扫描过 程中，光源S的位置用其与O'的距离l标记；探测器探元B的位置用其与O"的距离t标记。

的投影数据p(l,t)可表示为：且并对q(l,t)做关于变量l的傅里叶变换：

为函数f(x,y)的二维傅里叶变换。

步骤1：建立TV最小化重建模型：

(1)设探测器共有2M个探元，其索引为m，每个探元之间的长度间隔为Δt，且都能接收到光源的 辐射；设光源有2N个采样位置，其索引为n，采样步进为Δl。从而ω_x和σ分别离散化表达为和 于是，

(或f(k₁,k₂))为二维连续函数f(x,y)的离散化表达。图像重建可以进一步表达为寻找如下线 性方程组的解：

$F_{p}f=\hat{f},---\backslash *\mathrm{MERGEFORMAT}\left(5\right)$

其中，为观测到的傅里叶采样数据，为向量化的离散物体函数。组合系数矩阵 表示对f作离散伪极傅里叶变换。建立LCT的TV最小化重建模型，也就是求解如下带约束 条件的优化问题：

f^*＝arg min||f||_TV；

                      \*MERGEFORMAT(6)

s.t.F_pf＝F(p),

其中，为沿图像i方向的差分算子，F(p)表示对投影数据做一维 傅里叶变换。

步骤2：利用ADM最小化TV模型：

(1)建立无约束重建模型：采用ADM框架给出其求解算法。令D_if＝z_i，则该TV重建 模型可以转化为：

$\min \underset{i}{\Sigma }{\left|\right|z_i\left|\right|}_1+\frac{\lambda }{2}{\left|\right|F_{p}f-\hat{f}\left|\right|}_2^2,---\backslash *\mathrm{MERGEFORMAT}\left(7\right)$

s.t.D_if＝z_i,

其中为保真项因子，用来调整目标函数中观测数据的一致性程度。式为带约束优化模型，为了 将其转化为无约束优化问题，利用增广拉格朗日函数进行如下转化：

$\underset{f,z_i}{\min }{L}_A(f,z_i)=\frac{\lambda }{2}{\left|\right|F_{p}f-\hat{f}\left|\right|}_2^2+\underset{i}{\Sigma }(u_i^T(D_{i}f-z_i)+\frac{{\rho }_i}{2}{\left|\right|D_{i}f-z_i\left|\right|}_2^2+{\left|\right|z_i\left|\right|}_1),---\backslash *\mathrm{MERGEFORMAT}\left(8\right)$

其中标量二次项惩罚系数，为乘子。

(2)利用ADM求解模型：基于ADM框架采用分离变量的方法将转化为两个单目标优化问题，也即

$\left(\begin{array}{c}\underset{f}{\min }{L}_{z_i}\left(f\right)=\frac{\lambda }{2}{\left|\right|F_{p}f-\hat{f}\left|\right|}_2^2+\underset{i}{\Sigma }(u_i^T(D_{i}f-z_i)+\frac{{\rho }_i}{2}{\left|\right|D_{i}f-z_i\left|\right|}_2^2)\\ \underset{z_i}{\min }{L}_f\left(z_i\right)=u_i^T(D_{i}f-z_i)+\frac{{\rho }_i}{2}{\left|\right|D_{i}f-z_i\left|\right|}_2^2+{\left|\right|z_i\left|\right|}_1\end{array}\right).---\backslash *\mathrm{MERGEFORMAT}\left(9\right)$

采用ADM分别对子问题求取其最小值点，于是，基于ADM的求解迭代公式为：

$\left(\begin{array}{c}{f}^{(k+1)}={(\lambda F_p^H{F}_p+\underset{i}{\Sigma }{\rho }_i{D}_i^T{D}_i)}^{+}(\lambda F_p^H\hat{f}-\underset{i}{\Sigma }{D}_i^T(u_i^T-{\rho }_i{z}_i^{\left(k\right)}))\\ z_i^{(k+1)}=\max \left\{\right|D_i{f}^{\left(k\right)}+\frac{u_i^{\left(k\right)}}{{\rho }_i}|-\frac{1}{{\rho }_i},0\}·\mathrm{sgn}(D_i{f}^{\left(k\right)}+\frac{u_i^{\left(k\right)}}{{\rho }_i})\\ u_i^{(k+1)}=u_i^{\left(k\right)}+{\rho }_i(D_i{f}^{(k+1)}-z_i^{(k+1)})\end{array}\right).---\backslash *\mathrm{MERGEFORMAT}\left(10\right)$

步骤3：利用PPFFT实现图像空-频域变换：

(1)伪极傅里叶变换定义：LCT重建中所应用的F_pf及定义为：

$\left(\begin{array}{c}{F}_{p}f\stackrel{\Delta }{=}\frac{\Delta l^2}{2\pi }\underset{k_1=-N}{\overset{N-1}{\Sigma }}\underset{k_2=-N}{\overset{N-1}{\Sigma }}f(k_1\mathrm{\Delta l},k_2\mathrm{\Delta l})e^{-j(\frac{2\mathrm{\pi n}}{2N}{k}_1+\frac{2\mathrm{\pi m}}{2N}{k}_2·\frac{\mathrm{n\Delta t}}{D})}\\ =\hat{f}(\frac{\mathrm{n\pi }}{\mathrm{\Delta lN}},\frac{\mathrm{m\pi }}{\mathrm{\Delta lN}}·\frac{\mathrm{n\Delta t}}{D}),\end{array}\right)---\backslash *\mathrm{MERGEFORMAT}\left(11\right)$

及

$\left(\begin{array}{c}{F}_p^H\hat{f}\stackrel{\Delta }{=}\frac{\mathrm{\Delta t\Delta l}}{4{\left(\mathrm{N\Delta l}\right)}^2}\underset{n=-N}{\overset{N-1}{\Sigma }}\underset{m=-M}{\overset{M-1}{\Sigma }}\hat{f}(\frac{\mathrm{n\pi }}{\mathrm{\Delta lN}},\frac{\mathrm{m\pi }}{\mathrm{\Delta lN}}·\frac{\mathrm{n\Delta t}}{D})e^{j(\frac{2\mathrm{\pi n}}{2N}{k}_1+\frac{2\mathrm{\pi m}}{2N}{k}_2·\frac{\mathrm{n\Delta t}}{D})}\\ =\tilde{f}(k_1\mathrm{\Delta l},k_2\mathrm{\Delta l}).\end{array}\right)---\backslash *\mathrm{MERGEFORMAT}\left(12\right)$

现对F_pf分析并推导快速计算方法，F_pf实际上可简写为：

$\hat{f}[n,m]\stackrel{\Delta }{=}\frac{1}{2\pi }\underset{k_1=-N}{\overset{N-1}{\Sigma }}\underset{k_2=-N}{\overset{N-1}{\Sigma }}f(k_1,k_2)\exp (-j\frac{2\mathrm{\pi n}{k}_1}{2N}-j\frac{2\mathrm{\pi m}{k}_2}{2N}·\mathrm{n\alpha }),---\backslash *\mathrm{MERGEFORMAT}\left(13\right)$

其中为方便计算且不失一般性，已令Δt＝1,Δl＝1。

(2)对阵列f(k₁,k₂)的每列做2N点一维快速傅里叶变换：展开二重求和式，可以得到：

$\left(\begin{array}{c}\hat{f}[n,m]=\frac{1}{2\pi }\underset{k_2=-N}{\overset{N-1}{\Sigma }}\exp (-j\frac{2\pi }{2N}·\mathrm{m\alpha }·n{k}_2)\underset{k_1=-N}{\overset{N-1}{\Sigma }}f(k_1,k_2)\exp (-j\frac{2\pi }{2N}n{k}_1)\\ =C·\underset{k_2=-N}{\overset{N-1}{\Sigma }}\exp (-j\frac{2\pi }{2N}·\mathrm{m\alpha }·n{k}_2)\underset{k_1=-N}{\overset{N-1}{\Sigma }}f(k_1,k_2)\exp (-j\frac{2\pi }{2N}n(k_1+N)).\end{array}\right)---\backslash *\mathrm{MERGEFORMAT}\left(14\right)$

其中，为与求和过程无关的常量。分析内层关于变量k₁的求和式：

${\hat{f}}_1[n,k_2]=\underset{k_1=0}{\overset{2N-1}{\Sigma }}f(k_1,k_2)\exp (-j\frac{2\pi }{2N}n{k}_1),---\backslash *\mathrm{MERGEFORMAT}\left(15\right)$

其计算过程实际为针对阵列f(k₁,k₂)的每列做2N点一维DFT，其计算可使用FFT技术实现快速计算。

(3)对阵列的每行做2N点线性调频Z变换：外层对变量k₂的求和：

$\hat{f}[n,m]=\frac{1}{2\pi }\exp \left(j\frac{2\pi }{2N}\mathrm{nN}\right)\underset{k_2=-N}{\overset{N-1}{\Sigma }}{\hat{f}}_1[n,k_2]\exp (-j\frac{2\mathrm{\pi m}{k}_2}{2N}·\mathrm{n\alpha }).---\backslash *\mathrm{MERGEFORMAT}\left(16\right)$

此式中若nα＝1，则该式实际为对阵列的每行做2N点一维DFT，此时快速计算仍可以利用 FFT实现。当nα≠1时，式实际上是对序列做线性调频Z变换(chirp-Z transform,CZT)。令 nα＝η，将式简写为：

${\hat{f}}_n\left[m\right]=\underset{k_2=-N}{\overset{N-1}{\Sigma }}{\hat{f}}_{1,n}\left[k_2\right]\exp (-j\frac{2\mathrm{\pi m}{k}_2}{2N}·\eta ).---\backslash *\mathrm{MERGEFORMAT}\left(17\right)$

注意到$2m{k}_2=m^2+k_2^2-{(m-k_2)}^2,$将该关系式代入式可得：

$\left(\begin{array}{c}{\hat{f}}_n\left[m\right]=\underset{k_2=-N}{\overset{N-1}{\Sigma }}{\hat{f}}_{1,n}\left[k_2\right]\exp (-j\frac{\mathrm{\pi \eta }}{2N}·(m^2+k_2^2-{(m-k_2)}^2))\\ =\exp (-j\frac{\mathrm{\pi \eta }}{2N}{m}^2)\underset{k_2=-N}{\overset{N-1}{\Sigma }}{\hat{f}}_{1,n}\left[k_2\right]\exp (-j\frac{\mathrm{\pi \eta }}{2N}·k_2^2)\exp (j\frac{\mathrm{\pi \eta }}{2N}·{(m-k_2)}^2),\end{array}\right)---\backslash *\mathrm{MERGEFORMAT}\left(18\right)$

若定义则上式可以进一步改写为：

${\hat{f}}_n\left[m\right]=s\left[m\right]\underset{k_2=-N}{\overset{N-1}{\Sigma }}{\hat{f}}_{1,n}\left[k_2\right]·s\left[k_2\right]·s[m-k_2].---\backslash *\mathrm{MERGEFORMAT}\left(19\right)$

式表明可由三步操作算得：首先用序列s[k₂]乘以得到然后利用 与s[k₂]计算卷积，最后再乘以s[m]即可得到在的计算过程中，卷积运算占 用了大部分计算时间，利用一维FFT可以替代此过程：分别先对与s[k₂]计算一维FFT，对 应相乘后做一维IFFT即可(为了避免FFT带来的周期循环效应，在做一维FFT之前应先对各序列补零至 二倍长度)。利用FFT可以使原卷积运算O(N²)复杂度降为N log N。

步骤4：实现并运行算法，获得重建图像：

现给出基于伪极傅里叶变换总变分最小化技术，获得LCT重建图像的具体方法如下：

输入数据p,初始化λ,ρ_i＞0.且设定最大迭代次数N，且令 k＝0.

(1)数据预处理：

$\hat{f}(\omega ,-\mathrm{\omega t}/D)\leftarrow \frac{D/\sqrt{2\pi }}{\sqrt{D^2+t^2}}F\left(p(l-t,t)\right)$

执行如下迭代过程：

(2)更新f，公式如下：

$f^{(k+1)}\leftarrow {(\lambda F_p^H{F}_p+\underset{i}{\Sigma }{\rho }_i{D}_i^T{D}_i)}^{+}(\lambda F_p^H\hat{f}-\underset{i}{\Sigma }{D}_i^T(u_i^T-{\rho }_i{z}_i^{\left(k\right)}))$

(3)更新z_i，公式如下：

$z_i^{(k+1)}\leftarrow \max \left\{\right|D_i{f}^{\left(k\right)}+\frac{u_i^{\left(k\right)}}{{\rho }_i}|-\frac{1}{{\rho }_i},0\}·\mathrm{sgn}(D_i{f}^{\left(k\right)}+\frac{u_i^{\left(k\right)}}{{\rho }_i})$

(4)更新u_i，公式如下：

$u_i^{(k+1)}\leftarrow u_i^{\left(k\right)}+{\rho }_i(D_i{f}^{(k+1)}-z_i^{(k+1)})$

(5)k←k+1

若“未达到迭代设定的最大次数,”则返回第2步，否则停止。其中，F_p和由步骤4实现快速 计算。迭代达到最大迭代轮数退出时，f^(N)即为输出重建图像，重建结果如图3所示，误差下降曲线如图 5所示。