CN109374567A - 基于有序子集快速收敛的太赫兹成像方法 - Google Patents

Abstract

本发明公开了一种基于有序子集快速收敛的太赫兹成像方法，主要解决现有技术成像速度较慢的问题。其实现过程为：固定太赫兹时域光谱THz‑TDS系统不动，通过三维扫描架系统控制目标物体移动，通过扫描架系统中的转台控制目标旋转的方式录取数据；对录取数据采用频域能量求和的方法进行预处理获得预处理后数据；利用比尔定律对预处理后数据进行处理得到衰减系数线积分数据；对衰减系数线积分数据根据子集均衡条件划分有序子集；对有序子集数据采用二次收敛步长算法进行成像，并通过循环迭代更新图像直到达到收敛要求，实现太赫兹成像。本发明加快了图像收敛速度，提高了成像效率，可用于生物组织、医学场景中的太赫兹成像。

Description

基于有序子集快速收敛的太赫兹成像方法

技术领域

本发明属于成像技术领域，更进一步涉及一种太赫兹成像方法，可用于生物组织、医学场景中的太赫兹成像。

背景技术

太赫兹波通常是指频率范围在0.1THz到10THz区间的电磁波，位于毫米波与红外光之间，具有穿透性、大带宽、频谱指纹识别性、安全无损的优点。作为一种新的技术，太赫兹技术在过去的几十年中取得了突破性进展，越来越受到研究人员和工业界的重视，在安全检查、医学、生物、国防等领域有着广泛的引用。

太赫兹成像能够获得物体的多层信息，是太赫兹技术发展的一个重要方向。基于最大似然的期望最大化算法MLEM是一种经典成像方法，具有较好的成像质量，从原理上来讲成像结果优于滤波反投影算法。有序子集OSEM算法的提出大幅度提高了收敛速度，因此OSEM算法及其各种变形算法得到了广泛地应用。但是太赫兹固有的特性增加了算法的复杂性，而且利用OSEM算法进行太赫兹成像仍然需要多次迭代，降低了成像效率。因此太赫兹成像中迭代过程耗时较高，收敛速度需要进一步提高。

天津大学申请的专利“一种太赫兹波计算机辅助层析成像装置及方法”中采用基于高斯波束的滤波反投影算法进行了三维太赫兹成像。高斯波束光斑大小在瑞利距离内可以接近其聚焦光斑的大小，因此限定了物品的大小。同时，采用滤波反投影算法容易出现伪影，降低图像质量。而为了提升成像质量需要采集大量的数据，会造成成像效率的降低。

发明内容

本发明的目的在于针对上述现有技术的不足，提出一种基于有序子集的快速收敛的太赫兹成像方法，以降低时间成本，提高成像效率。

本发明的技术思路是通过将有序子集和二次步长收敛成像方法的结合，在保证成像质量的同时进一步提高算法的收敛速度，从而提高成像效率。其实现方案包括如下：

(1)通过太赫兹时域光谱THz-TDS系统录取数据：

将TDS系统固定不动，通过三维扫描架系统实现目标物体相对于TDS系统的移动，通过扫描架系统中的转台实现目标的旋转；设初始角度为0度，角度采样间隔为Δθ，旋转次数为N，得到总的成像角度大小为θ＝NΔθ；

每旋转每一个角度，通过扫描架控制目标在垂直于波束的平面内进行水平移动，设移动间隔为Δl，水平移动次数为M，每移动一次则测量一次数据，每次测量数据的长度为K，得到录取的数据大小为N×M×K；

(2)对录取数据进行预处理：即先对数据进行快速傅里叶变换FFT，将数据从时域变换到频域，再用500-520GHz的太赫兹频段信号的积分和作为最终的信号，得到大小为N×M的预处理后数据S；

(3)根据比尔定律，对预处理后数据S进行归一化后取对数，再取负值，得到衰减系数线积分数据：B＝-log(S/max(S))；

(4)将衰减系数线积分数据B根据子集均衡条件按照角度划分为如下有序子集：

其中，k为第k个子集，k＝1,2...N_sub，N_sub为划分的子集个数，H为每个子集包含的角度的个数：H＝N/N_sub，θ_k,h为第k个子集中第h个角度，θ_k,h＝(k-1)Δθ+(h-1)N_subΔθ，h＝1,2...H，B_θk,h为B中对应角度为θ_k,h的数据；

(5)对子集数据利用二次步长收敛的期望最大化算法进行成像：

(5a)设定图像的初始值:其中J＝M²为图像像素个数；

(5b)利用第1个有序子集数据更新图像，得到更新后的图像

x⁽¹⁾＝x⁽⁰⁾·D_c

其中，x⁽¹⁾为更新后的图像，D_c为修正因子；

(5c)将重建出的x⁽¹⁾作为初始值代入下一个有序子集中，重复步骤(5b)，循环更新图像直到做完第N_sub个子集，完成本次迭代。

(6)将以上一次迭代的结果作为下一次迭代初值，重复(5b)到(5c)进行新一轮迭代，直到得到符合收敛要求的结果为止。

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

1.本发明由于采用了有序子集的数据划分方法，克服了现有技术采用最大似然期望值最大法收敛较慢的问题，提高了图像的收敛速度；

2.本发明由于采用了二次步长收敛成像方法，即使在数据量较少的情况下也可以保证成像质量，同时进一步提升收敛速度，提高成像效率。

附图说明

图1为本发明的实现流程图；

图2为本发明中用于数据采集的THz-TDS系统结构示意图；

图3为本发明仿真使用的原图；

图4为用本发明和常规方法的仿真成像结果图；

图5为分别用本发明和常规方法得到的仿真结果的相关系数变化曲线；

图6为本发明实测物品；

图7为用本发明和常规方法分别得到的实测数据成像结果图；

图8为分别用本发明和常规方法得到的实测结果的相关系数变化曲线。

具体实施方式

下面结合附图对本发明的实施例和效果做进一步的描述。

参照图1，本发明的实现步骤如下：

步骤1，数据采集。

数据采集可以采用太赫兹连续波系统或太赫兹时域光谱THz-TDS系统进行。但是由于太赫兹连续波系统只能采集到信号的强度信息，不能获得相位信息，无法获得某一频率或某一段频率的信号，而太赫兹时域光谱THz-TDS系统不仅能够获得物品的幅度信息还能获得相位信息，因此本实例采用太赫兹时域光谱THz-TDS系统录取数据。

如图2所示，太赫兹时域光谱THz-TDS系统的太赫兹信号是由激光激发，其中一路经反光镜传输到THz发射极，再经透镜聚焦到样品上，透过样品的太赫兹信号再经透镜传到THz探测器。另一路信号经光学扫描延迟线调整后传输到THz探测器作为参考信号，两路信号在THz探测器进行混频后再经过电流前置放大器输出到数字信号处理设备。

本步骤的数据采集实现如下：

首先，保持太赫兹时域光谱THz-TDS系统固定不动，通过扫描架系统中的转台实现目标的旋转，每旋转一个角度，通过扫描架控制目标在垂直于波束的平面内进行水平移动，设移动间隔为Δl，水平移动次数为M，每移动一次则录取一次数据，每次录取数据的长度为K，得到录取数据的大小为M×K，对应的成像区域大小为L×L,其中L＝MΔl；

一共旋转N次，得到成像角度的大小为θ＝NΔθ，总的录取数据大小为N×M×K，其中，Δθ为角度采样间隔。

步骤2，对录取数据进行预处理。

由于太赫兹时域光谱THz-TDS系统录取的数据包含多种信息，需要在成像之前进行预处理，现有的预处理方法包括：最大值法、最小值法、峰峰值法、时间延迟法、时域信号能量法、频率选取法、频域能量求和法；

图像分辨率与成像频率有关，因此为获得较高的分辨率需要利用较高的频率进行成像；信号信噪比对成像质量有重要影响，为了提高信号的信噪比，采用能量求和的方法进行成像，因此本发明采用频域能量求和法进行预处理，即先对数据进行快速傅里叶变换FFT将数据从时域变换到频域，再采用500GHz-520GH的太赫兹频段的积分和作为最终的信号，得到大小为N×M的预处理后数据S。

步骤3，计算衰减系数线积分数据。

本实例需要利用衰减系数的线积分进行成像，因此需要计算吸收率线积分数据:

根据比尔定律有关系式：S＝S₀ exp(-B)，其中S₀为S的最大值，B为衰减系数线积分，exp(·)表示以自然常数为底的指数函数；

利用上式对预处理后数据S进行归一化后取对数，再取负值，得到衰减系数线积分数据：

B＝-ln(S/max(S))。

步骤4，将衰减系数线积分数据B划分为有序子集。

将衰减系数线积分数据B根据子集均衡条件按照角度划分有序子集，其实现如下：

(4a)根据旋转角度的个数计算需要划分的子集个数：N_sub＝<sqrt(N/2)>，式中sqrt表示开平方，<·>表示取最接近N的因子的操作；

(4b)根据子集均衡条件，将旋转角度划分为N_sub个有序子集，其中第k个有序子集对应的角度为：{θ_k,h}＝{(k-1)Δθ+(h-1)N_subΔθ}，k＝1,2,...N_sub，h＝1,2...H，H＝N/N_sub为每个子集包含的角度的个数；

(4c)根据第k个有序子集对应的角度，得到第k个有序子集数据为：

其中为B中对应角度为θ_k,h的数据。

步骤5，对有序子集数据进行成像。

现有成像方法有最大似然期望值最大MLEM算法，二次步长收敛的期望最大化算法，基于罚函数的期望值最大法，基于贝叶斯公式的期望值最大法。

最大似然期望值最大MLEM算法采用一次步长，收敛速度较慢；基于罚函数的期望值最大法和基于贝叶斯公式的期望值最大法能够防止迭代次数过多导致的图像退化问题，但是增加了额外的计算量，导致时间开销较大；而二次步长收敛的期望最大化算法采用了二次步长收敛，比最大似然期望值最大MLEM算法采用一次步长收敛更快，而且在迭代次数较少的情况下可以忽略图像退化的问题，因此为了加快收敛速度，提高成像效率，本实例采用二次步长收敛的期望最大化算法进行成像，其实现如下：

(5a)设初始化图像为x⁰，其中x⁰中的元素为:J＝M²为图像的像素个数；

(5b)利用第1个有序子集数据更新初始化图像x⁰，得到更新后的图像x¹：

(5b1)计算图像与波束模型的卷积：u＝conv(x⁰,g)，其中conv表示卷积操作，g为高斯波束模型：式中，I₀为波束束腰位置的密度，w₀为波束束腰的半径，λ为波长，r为距离Z轴的距离，w(e)为距离Z轴零点位置为e处的波束半径：为瑞利距离；

(5b2)计算第1个有序子集数据所对应的估计值：

其中，I＝H·M为每个有序子集中包含的射线数，u_j为卷积后图像u中的第j个像素，a_ij为图像中第j个像素对第i射线的权重，如果射线穿过该像素则权重为1，否则为0；

(5b3)计算更新步长：式中，p_i为第1个有序子集中第i个衰减系数线积分数据；

(5b4)计算像素的修正值：

(5b5)计算图像的修正因子：D_c＝Deconv(C,g)，式中Deconv表示解卷积操作，C为由C_j组成的大小为M×M的修正矩阵；

(5b6)利用修正因子对初始图像进行修正，得到更新后的图像：x⁽¹⁾＝x⁽⁰⁾·D_c。

(5c)将修正后的图像x⁽¹⁾作为初始值代入下一个有序子集中，重复步骤(5b1)到(5b6)，循环更新图像直到做完第N_sub个有序子集，完成第一轮迭代。

步骤6，以上一次迭代的结果作为下一次迭代的初值，重复(5b)到(5c)进行新一轮迭代，直到符合收敛要求的结果为止，完成太赫兹成像。

下面结合仿真结果和实测结果对本发明的效果做进一步的说明：

1、仿真实验条件：

本发明的仿真实验是在主频2.1GHz、内存64GB的Dell计算机和MATLAB R2017a的环境中进行编程实现的。

仿真使用的原图像如图3，其大小为128*128像素。

2、仿真内容与结果：

仿真1，设成像角度范围为180°，角度采样间隔Δθ＝1°，有序子集个数N_sub＝9，水平移动次数M＝185，分别采用本发明和现有技术中的有序子集最大似然估计OSEM算法对图3进行100次的迭代成像仿真，结果如图4所示，其中：

图4(a)为用有序子集最大似然估计OSEM算法分别迭代1、5、10、20、50、100次的成像结果仿真结果；

图4(b)为用本发明分别迭代1、5、10、20、50、100次的成像结果仿真结果；

从图4可以看出本发明的收敛速度比有序子集最大似然估计OSEM算法更快，因此本发明方法的成像效率更高。

仿真2，对用本发明方法和有序子集最大似然估计OSEM算法得到的100幅图像分别计算相邻图像间的相关系数，结果如图5所示，其中实线代表本发明方法，虚线代表有序子集最大似然估计OSEM算法。

图像间的相关系数可以表示图像的变化快慢，相关系数曲线越陡峭说明收敛速度越快，从图5看出显然红色曲线更加陡峭，说明本发明方法的收敛速度更快，成像效率更高。

仿真3，设成像角度范围为180°，角度采样间隔Δθ＝2°，有序子集个数N_sub＝6，水平移动次数M＝30，水平移动间隔Δl＝2mm，分别采用本发明和现有技术中的有序子集最大似然估计OSEM算法对图6所示的两支笔进行迭代100次的实测数据成像，结果如图7所示，其中：

图7(a)为用有序子集最大似然估计OSEM算法分别迭代1、5、10、20、50、100次的实测数据成像处理结果；

图7(b)为用本发明分别迭代1、5、10、20、50、100次的实测数据成像处理结果；

图7中，采用OSEM算法的结果第一幅图比较模糊，而本发明方法的结果第一张图像就很清晰，说明本发明方法的收敛速度更快，因此本发明方法的成像效率更高。

仿真4，对用本发明方法和有序子集最大似然估计OSEM算法得到的100幅实测图像分别计算相邻图像间的相关系数，结果如图8所示，其中实线代表本发明方法，虚线代表有序子集最大似然估计OSEM算法。

从图8中可见实线更陡峭，说明本发明方法的收敛速度更快，成像效率更高。

以上仿真实验结果和实测数据处理结果均验证了本发明的正确性，有效性和可靠性。

Claims (3)

1.基于有序子集快速收敛的太赫兹成像方法，包括如下：

(1)通过太赫兹时域光谱THz-TDS系统录取数据：

将TDS系统固定不动，通过三维扫描架系统实现目标物体相对于TDS系统的移动，通过扫描架系统中的转台实现目标的旋转；设初始角度为0度，角度采样间隔为Δθ，旋转次数为N，得到总的成像角度大小为θ＝NΔθ；

每旋转一个角度，通过扫描架控制目标在垂直于波束的平面内进行水平移动，设移动间隔为Δl，水平移动次数为M，每移动一次则录取一次数据，每次录取数据的长度为K，得到录取的数据大小为N×M×K；

(2)对录取数据进行预处理：即先对数据进行快速傅里叶变换FFT，将数据从时域变换到频域，再用500-520GHz的太赫兹频段信号的积分和作为最终的信号，得到大小为N×M的预处理后数据S；

(3)根据比尔定律，对预处理后数据S进行归一化后取对数，再取负值，得到衰减系数线积分数据：B＝-ln(S/max(S))；

(4)将衰减系数线积分数据B根据子集均衡条件按照角度划分为如下有序子集：

其中，k为第k个子集，k＝1,2...N_sub，N_sub为划分的子集个数，H为每个子集包含的角度的个数：H＝N/N_sub，θ_k,h为第k个子集中第h个角度，θ_k,h＝(k-1)Δθ+(h-1)N_subΔθ，h＝1,2...H，为B中对应角度为θ_k,h的数据；

(5)对有序子集数据利用二次步长收敛的期望最大化算法进行成像：

(5a)设初始化图像为x⁰，其中x⁰中的元素为:其中J＝M²为图像像素个数；

(5b)利用第1个子集数据更新图像，得到更新后的图像

x⁽¹⁾＝x⁽⁰⁾·D_c

其中，x⁽¹⁾为更新后的图像，D_c为修正因子；

(5c)将重建出的x⁽¹⁾作为初始值代入下一个子集中，重复步骤(5b)，循环更新图像直到做完第N_sub个子集，完成本次迭代。

(6)将以上一次迭代的结果作为下一次迭代初值，重复(5b)到(5c)进行新一轮迭代，直到得到符合收敛要求的结果为止。

2.根据权利要求1所述的方法，其中(5b)中的为修正因子D_c，其计算如下：

(5b1)计算图像与波束模型的卷积：u＝conv(x⁰,g)，其中conv表示卷积操作，x0为初始化后的图像，g为高斯波束模型：式中，I₀为波束束腰位置的密度，w₀为波束束腰的半径，λ为波长，r为距离Z轴的距离，w(e)为距离Z轴零点位置为e处的波束半径： 为瑞利距离；

(5b2)计算第1个有序子集数据所对应的估计值：

其中，I＝H·M为有序子集中包含的射线数，u_j为卷积后图像u中的第j个像素，a_ij为图像中第j个像素对第i射线的权重，如果射线穿过该像素则权重为1，否则为0；

(5b3)计算更新步长：式中，p_i为有序子集中第i个衰减系数线积分数据；

(5b4)计算像素的修正值：

(5b5)计算图像的修正因子：D_c＝Deconv(C,g)，式中Deconv表示解卷积操作，C为由C_j组成的大小为M×M的修正矩阵。

3.根据权利要求1所述的方法，其中(4)中将衰减系数线积分数据B根据子集均衡条件按照角度进行有序子集划分，其实现如下：

(4a)根据旋转角度的个数计算需要划分的子集个数：N_sub＝＜sqrt(N/2)>，式中sqrt表示开平方，<·>表示取最接近N的因子的操作；

(4b)根据子集均衡条件，将旋转角度划分为N_sub个有序子集，其中第k有序个子集对应的角度为：{θ_k,h}＝{(k-1)Δθ+(h-1)N_subΔθ}，h＝1,2...H；

(4c)根据第k个子集对应的角度，得到第k个有序子集数据为：其中为B中对应角度为θ_k,h的数据。