阅读说明：本技术 一种基于费马原理的超声ct图像重建方法及系统 (Ultrasonic CT image reconstruction method and system based on Fermat principle ) 是由 尉迟明 丁明跃 方小悦 武云 宋俊杰 周亮 张求德 周权 于 2019-08-27 设计创作，主要内容包括：本发明属于功能成像技术领域,公开了一种基于费马原理的超声CT图像重建方法及系统,其中基于费马原理的超声CT声速重建方法包括步骤：(1)渡越时间的提取；(2)准备进入迭代；(3)根据有限差分方法计算从每个发射阵元出发到成像区域每个像素点的时间τ,并根据费马原理计算胖射线路径,其中胖射线路径的宽度随迭代次数的增加而收窄；(4)反问题的求解,并更新声速值；(5)判断迭代是否终止。本发明通过对方法整体流程进行改进,通过路径的优化,基于费马原理,尤其是变参数的、路径由宽至窄的变路径宽度的胖射线路径,使得本发明中超声CT声速及衰减系数重建方法及系统具有快速、稳定、成像效果更好的特点。(The invention belongs to the technical field of functional imaging and discloses an ultrasonic CT image reconstruction method and system based on the Fermat principle, wherein the ultrasonic CT sound velocity reconstruction method based on the Fermat principle comprises the following steps: (1) extracting the transit time; (2) preparing to enter iteration; (3) calculating the time tau from each emission array element to each pixel point of an imaging region according to a finite difference method, and calculating a fat ray path according to the Fermat principle, wherein the width of the fat ray path is narrowed along with the increase of the iteration times; (4) solving an inverse problem and updating a sound velocity value; (5) it is determined whether the iteration is terminated. The method and the system for reconstructing the ultrasonic CT sound velocity and the attenuation coefficient have the characteristics of rapidness, stability and better imaging effect by improving the overall process of the method, optimizing the path and based on the Fermat principle, particularly the fat ray path with the variable parameter and the variable path width from wide to narrow.)

一种基于费马原理的超声CT图像重建方法及系统

技术领域

本发明属于功能成像技术领域，更具体地，涉及一种基于费马原理的超声CT图像重建方法及系统，属于超声断层成像中透射式超声成像模式，可实现超声CT声速及衰减系数的重建。

背景技术

超声CT是一种断层成像模式，采用超声探头发射和接收信号，产生反射数据和透射数据，利用这些数据重建出不同模态的超声断层图像，以便更好地观测物体内部的参数信息。超声CT检测具有价格低廉、对人体无辐射等优点，随着探头加工工艺和计算机高性能运算的快速发展，近些年来超声断层成像技术又再次成为了工业领域和医学领域的研究热点。

利用超声CT系统采集的数据可以重建出声速、衰减、密度等多种参数信息，属于功能像的领域。在这里，我们采用透射数据。有研究表明，在病变初期，病变组织功能参数的变化要早于结构变化。因此，超声CT功能成像对病变的早期诊断具有重要辅助意义。

以声速这一功能参数为例，超声CT声速重建方法包括基于射线理论的重建算法和基于波动理论的重建算法。基于波动理论的方法成像分辨能力更好，但是容易受微小误差扰动，因而鲁棒性不高，并且运算量很大，对系统精度和数据的性噪比要求更高，目前还不适用于实际应用。基于射线理论的重建算法模型更简单，稳定性更高，并且运算量较小，目前来看是一类高效的、稳定的、适用于临床的声速重建方法。

近些年来国内外对超声CT声速重建方法的研究逐渐活跃起来。基于射线理论的研究领域，射线路径的追踪方法众多，针对不同应用场景适用性不同，成像效果也不同。有采用直线路径的，也有采用折线路径的。采用直线路径是模仿X线CT重建方法，理论上来说，超声的传播更近似于折线。但是，由于临床超声探头的带宽有限，实际上超声的传播路径也并非是一条折线。因而路径的追踪是该类方法的难点。

发明内容

针对现有技术的以上缺陷或改进需求，本发明的目的在于提供一种基于费马原理的超声CT图像重建方法及系统，其中通过对方法的整体流程进行改进，通过路径的优化，基于费马原理，尤其是变参数的、路径由宽至窄的变路径宽度的胖射线路径，一方面，不仅考虑射线路径上的点对接收信息的影响，同时还考虑射线领域上的其它点对接收信息的影响，利用特定的路径追踪过程，使得本发明中超声CT声速及衰减系数重建方法及系统具有快速、稳定、成像效果更好的特点，另一方面，尤为重要的是，实现了低分辨率到高分辨率加深重建的方法，使重建过程更加稳定，重建分辨率更高，减少重建伪影。

为实现上述目的，按照本发明的一个方面，提供了一种基于费马原理的超声CT声速重建方法，其特征在于，包括以下步骤：

(1)渡越时间的提取：

针对待测对象，采用环形阵列探头，按发射阵元预先设定的顺序，基于其中一个发射阵元，采集满足预先设定接收阵元角度要求范围内全部接收阵元接收到的透射波数据；接着，按接收阵元预先设定的顺序，对于其中每一条透射波数据，在透射波数据上预估的渡越时间点附近选取匹配窗，窗长度预先设定为N，接着，以当前检索点为分界点，将窗分为两段，利用公式(1)计算当前检索点的AIC值，然后，对窗内的点依次检索，记AIC值最小的点为渡越时间点；

AIC(k)--k*log(var(d(1,k)))+(N-k-1)*log(var(d(k+1，N))) (1)

该公式(1)中，var代表方差，d(1，k)表示第一个点到第k个点的数据；

如此得到针对同一个发射阵元、且按接收阵元顺序对应的一系列渡越时间，并最终由此得到针对全部发射阵元的、且按发射阵元及接收阵元顺序对应的渡越时间总集合，将该渡越时间总集合列向量化即可得到渡越时间列向量[T_tof]；

(2)给定预先设定的声速值作为初始声速值，记迭代次数为1，进入第一次迭代；

(3)根据有限差分方法计算从每个发射阵元出发到成像区域每个像素点的时间τ，并根据费马原理计算胖射线路径，具体是：

将成像区域划分成m₁×m₂个网格，其中m₁、m₂均为预先设定的正整数；

然后，按所述发射阵元预先设定的顺序，对于其中一个发射阵元，将与该发射阵元相对应的、且满足预先设定接收阵元角度要求范围内的全部接收阵元，按所述接收阵元预先设定的顺序，选取其中一个接收阵元，得到一对发射阵元-接收阵元；接着，就这一对发射阵元-接收阵元计算声波从该发射阵元到该接收阵元经过的胖射线路径；具体的，采用有线差分法计算每一个发射阵元发射声波，传播到成像区域每个像素点的时间τ，以成像区域的每一个像素点P点为对象，若满足

τ_SP+τ_RP-τ_SR≤Δt (2)

则，该网格的胖射线路径值为有效值，记为预先设定的值；否则为无效值，记为0；从而得到这一对发射阵元-接收阵元下、对应于成像区域各个网格的一系列胖射线路径值，接着再将这一系列胖射线路径值行向量化，即可得到胖射线路径行向量；

所述公式(2)中，S代表发射阵元，R代表接收阵元，S和R互为发射接收，τ_SP代表声波由S出发、到达P的时间，τ_RP代表声波由R出发、到达P的时间，τ_SR代表声波由S出发、到达R的时间；Δt为按迭代次数从预先设定的时间阈值数列{Δt}中选取的第i个时间阈值，i与所述迭代次数在数值上相等，该时间阈值数列{Δt}中的时间阈值是按由大到小的顺序排列，且该时间阈值数列{Δt}中的任意一个时间阈值均与探头中心频率f相关；

如此即可得到针对同一个发射阵元、且按接收阵元顺序对应的一系列胖射线路径行向量，并最终得到针对全部发射阵元的、且按发射阵元及接收阵元顺序对应的胖射线路径行向量总集合；将每个行向量视为一个元素整体，将该胖射线路径行向量总集合中的所有元素排列成一列，即可得到行位置与发射阵元及接收阵元顺序相对应的胖射线路径总矩阵[L]；

(4)反问题的求解，并更新声速值：

基于所述步骤(1)得到的渡越时间列向量[T_tof]和所述步骤(3)得到的胖射线路径总矩阵[L]构建如式(3)所示的路径-慢度-时间方程组：

[L][S]＝[T_tof] (3)

该式(3)中，[S]是总行数为m₁×m₂、总列数为1的待求解的慢度列向量；

然后，采用随机梯度下降法解方程组求解得到[S]，接着对[S]中的每个元素取倒数，即可得到针对待测对象的速度重建值，然后以该速度重建值更新声速值；同时更新迭代次数，使迭代次数加1；

(5)判断迭代次数是否小于预先设定的迭代次数阈值，若小于，则重复步骤(3)和步骤(4)；若不小于，则终止迭代，并输出最终的声速值；

此外，所述预先设定的时间阈值数列{Δt}中时间阈值的总个数大于等于所述预先设定的迭代次数阈值。

作为本发明的进一步优选，所述步骤(1)中，所述预先设定接收阵元角度要求范围是以这一发射阵元正对面为基准、以所述环形阵列，探头的圆心为圆心、圆心角角度从-α到+α共计2α角度范围内的接收阵元，其中，α为预先设定的正角度，并且2α满足90°～270°。

作为本发明的进一步优选，所述步骤(2)中，是以水中的声速值作为所述初始声速值；

所述步骤(3)中，网格的胖射线路径值中，有效值是记为1。

按照本发明的另一方面，本发明提供了一种基于费马原理的超声CT衰减系数重建方法，其特征在于，包括以下步骤：

(1)能量参数的提取：

针对待测对象，采用环形阵列探头，按发射阵元预先设定的顺序，基于其中一个发射阵元，采集满足预先设定接收阵元角度要求范围内全部接收阵元接收到的透射波的能量参数，并将这些能量参数按接收阵元预先设定的顺序排列，如此得到针对同一个发射阵元、且按接收阵元顺序对应的一系列能量参数；并最终由此得到针对全部发射阵元的、且按发射阵元及接收阵元顺序对应的能量参数总集合，将该参量参数总集合列向量化即可得到能量参数列向量[P]；

(2)给定预先设定的衰减系数值作为初始衰减系数值，记迭代次数为1，进入第一次迭代；

(3)根据有限差分方法计算从每个发射阵元出发到成像区域每个像素点的时间τ，并根据费马原理计算胖射线路径，具体是：

将成像区域划分成m₁×m₂个网格，其中m₁、m₂均为预先设定的正整数；

然后，按所述发射阵元预先设定的顺序，对于其中一个发射阵元，将与该发射阵元相对应的、且满足预先设定接收阵元角度要求范围内的全部接收阵元，按所述接收阵元预先设定的顺序，选取其中一个接收阵元，得到一对发射阵元-接收阵元；接着，就这一对发射阵元-接收阵元计算声波从该发射阵元到该接收阵元经过的胖射线路径；具体的，是以成像区域的每一个网格中心点P点为对象，若满足

τ_SP+τ_RP-τ_SR≤Δt (2)

则，该网格的胖射线路径值为有效值，记为预先设定的值；否则为无效值，记为0；从而得到这一对发射阵元-接收阵元下、对应于成像区域各个网格的一系列胖射线路径值，接着再将这一系列胖射线路径值行向量化，即可得到胖射线路径行向量；

所述公式(2)中，S代表发射阵元，R代表接收阵元，S和R互为发射接收，τ_SP代表声波由S出发、到达P的时间，τ_RP代表声波由R出发、到达P的时间，τ_SR代表声波由S出发、到达R的时间；Δt为按迭代次数从预先设定的时间阈值数列{Δt}中选取的第i个时间阈值，i与所述迭代次数在数值上相等，该时间阈值数列{Δt}中的时间阈值是按由大到小的顺序排列，且该时间阈值数列{Δt}中的任意一个时间阈值均与探头中心频率f相关；

如此即可得到针对同一个发射阵元、且按接收阵元顺序对应的一系列胖射线路径行向量，并最终得到针对全部发射阵元的、且按发射阵元及接收阵元顺序对应的胖射线路径行向量总集合；将每个行向量视为一个元素整体，将该胖射线路径行向量总集合中的所有元素排列成一列，即可得到行位置与发射阵元及接收阵元顺序相对应的胖射线路径总矩阵[L]；

(4)反问题的求解，并更新衰减系数值：

基于所述步骤(1)得到的能量参数列向量[P]和所述步骤(3)得到的胖射线路径总矩阵[L]构建如式(4)所示的路径-衰减-能量参数方程组：

[L][A]＝[P] (4)

该式(4)中，[A]是总行数为m₁×m₂、总列数为1的待求解的衰减系数列向量；

然后，采用随机梯度下降法解方程组求解得到[A]，这样即可得到针对待测对象的衰减系数重建值，然后以该衰减系数重建值更新衰减系数值；同时更新迭代次数，使迭代次数加1；

(5)判断迭代次数是否小于预先设定的迭代次数阈值，若小于，则重复步骤(3)和步骤(4)；若不小于，则终止迭代，并输出最终的衰减系数值；

此外，所述预先设定的时间阈值数列{Δt}中时间阈值的总个数大于等于所述预先设定的迭代次数阈值。

作为本发明的进一步优选，所述步骤(1)中，所述预先设定接收阵元角度要求范围是以这一发射阵元正对面为基准、以所述环形阵列探头的圆心为圆心、圆心角角度从-α到+α共计2α角度范围内的接收阵元，其中，α为预先设定的正角度，并且2α满足90°～270°；优选的，2α＝180°。

作为本发明的进一步优选，所述步骤(2)中，是以水中的衰减系数值作为所述初始衰减系数值；

所述步骤(3)中，网格的胖射线路径值中，有效值是记为1。

按照本发明的又一方面，本发明提供了利用上述基于费马原理的超声CT声速重建方法的超声CT图像重建方法，其特征在于，该方法利用上述基于费马原理的超声CT声速重建方法，还包括以下步骤：

(6)成像：将所述步骤(5)得到的最终的声速值排列成二维矩阵，形成m₁×m₂的声速值二维矩阵，然后基于该矩阵得到二维像素图，该二维像素图中的每个像素与声速值相对应；

所述二维像素图是对声速值进行对数压缩、灰度映射、并最终显示得到。

按照本发明的又一方面，本发明提供了利用上述基于费马原理的超声CT衰减系数重建方法的超声CT图像重建方法，其特征在于，该方法利用上述基于费马原理的超声CT衰减系数重建方法，还包括以下步骤：

(6)成像：将所述步骤(5)得到的最终的衰减系数值排列成二维矩阵，形成m₁×m₂的衰减系数值二维矩阵，然后基于该矩阵得到二维像素图，该二维像素图中的每个像素与衰减系数值相对应；

所述二维像素图是对衰减系数值进行对数压缩、灰度映射、并最终显示得到。

按照本发明的再一方面，本发明提供了一种基于费马原理的超声CT声速重建系统，其特征在于，该系统包括：

渡越时间提取模块用于：针对待测对象，采用环形阵列探头，按发射阵元预先设定的顺序，基于其中一个发射阵元，采集满足预先设定接收阵元角度要求范围内全部接收阵元接收到的透射波数据；按接收阵元预先设定的顺序，对于其中每一条透射波数据，在透射波数据上预估的渡越时间点附近选取匹配窗，窗长度预先设定为N，以当前检索点为分界点，将窗分为两段，利用公式(1)计算当前检索点的AIC值，对窗内的点依次检索，记AIC值最小的点为渡越时间点；

AIC(k)＝k*log(var(d(1,k)))+(N-k-1)*log(var(d(k+1，N))) (1)

该公式(1)中，var代表方差，d(1，k)表示第一个点到第k个点的数据；

从而得到针对同一个发射阵元、且按接收阵元顺序对应的一系列渡越时间；并最终得到针对全部发射阵元的、且按发射阵元及接收阵元顺序对应的渡越时间总集合，将该渡越时间总集合列向量化即可得到渡越时间列向量[T_tof]；

迭代计算功能模块包括胖射线路径计算模块和反问题求解模块，其中，

所述胖射线路径计算模块用于：根据有限差分方法计算从每个发射阵元出发到成像区域每个像素点的时间τ，并根据费马原理计算胖射线路径，具体是：

将成像区域划分成m₁×m₂个网格，其中m₁、m₂均为预先设定的正整数；

按所述发射阵元预先设定的顺序，对于其中一个发射阵元，将与该发射阵元相对应的、且满足预先设定接收阵元角度要求范围内的全部接收阵元，按所述接收阵元预先设定的顺序，选取其中一个接收阵元，得到一对发射阵元-接收阵元；就这一对发射阵元-接收阵元计算声波从该发射阵元到该接收阵元经过的胖射线路径；具体的，是以成像区域的每一个网格中心点P点为对象，若满足

τ_SP+τ_RP-τ_SR≤Δt (2)

则，该网格的胖射线路径值为有效值，记为预先设定的值；否则为无效值，记为0；从而得到这一对发射阵元-接收阵元下、对应于成像区域各个网格的一系列胖射线路径值，将这一系列胖射线路径值行向量化，从而得到胖射线路径行向量；

所述公式(2)中，S代表发射阵元，R代表接收阵元，S和R互为发射接收，τ_SP代表声波由S出发、到达P的时间，τ_RP代表声波由R出发、到达P的时间，τ_SR代表声波由S出发、到达R的时间；Δt为按迭代次数从预先设定的时间阈值数列{Δt}中选取的第i个时间阈值，i与所述迭代次数在数值上相等，该时间阈值数列{Δt}中的时间阈值是按由大到小的顺序排列，且该时间阈值数列{Δt}中的任意一个时间阈值均与探头中心频率f相关；

从而得到针对同一个发射阵元、且按接收阵元顺序对应的一系列胖射线路径行向量，并最终得到针对全部发射阵元的、且按发射阵元及接收阵元顺序对应的胖射线路径行向量总集合；将每个行向量视为一个元素整体，将该胖射线路径行向量总集合中的所有元素排列成一列，得到行位置与发射阵元及接收阵元顺序相对应的胖射线路径总矩阵[L]；

所述反问题求解模块用于：基于所述渡越时间列向量[T_tof]和所述胖射线路径总矩阵[L]构建如式(3)所示的路径-慢度-时间方程组：

[L][S]＝[T_tof] (3)

该式(3)中，[S]是总行数为m₁×m₂、总列数为1的待求解的慢度列向量；

采用随机梯度下降法解方程组求解得到[S]，对[S]中的每个元素取倒数，从而得到针对待测对象的速度重建值，该速度重建值即为更新后的声速值；

判断是否终止的功能模块用于：判断迭代次数是否小于预先设定的迭代次数阈值，若小于则重复迭代；若不小于，则终止迭代，并输出最终的声速值；

此外，所述预先设定的时间阈值数列{Δt}中时间阈值的总个数大于等于所述预先设定的迭代次数阈值。

按照本发明的最后一方面，本发明提供了一种基于费马原理的超声CT衰减系数重建系统，其特征在于，该系统包括：

能量参数提取模块用于：针对待测对象，采用环形阵列探头，按发射阵元预先设定的顺序，基于其中一个发射阵元，采集满足预先设定接收阵元角度要求范围内全部接收阵元接收到的透射波的能量参数，并将这些能量参数按接收阵元预先设定的顺序排列，得到针对同一个发射阵元、且按接收阵元顺序对应的一系列能量参数；并最终得到针对全部发射阵元的、且按发射阵元及接收阵元顺序对应的能量参数总集合，将该参量参数总集合列向量化即可得到能量参数列向量[P]；

迭代计算功能模块包括胖射线路径计算模块和反问题求解模块，其中，所述胖射线路径计算模块用于：根据有限差分方法计算从每个发射阵元出发到成像区域每个像素点的时间τ，并根据费马原理计算胖射线路径，具体是：

将成像区域划分成m₁×m₂个网格，其中m₁、m₂均为预先设定的正整数；

按所述发射阵元预先设定的顺序，对于其中一个发射阵元，将与该发射阵元相对应的、且满足预先设定接收阵元角度要求范围内的全部接收阵元，按所述接收阵元预先设定的顺序，选取其中一个接收阵元，得到一对发射阵元-接收阵元；就这一对发射阵元-接收阵元计算声波从该发射阵元到该接收阵元经过的胖射线路径；具体的，是以成像区域的每一个网格中心点P点为对象，若满足

τ_SP+τ_RP-τ_SR≤Δt (2)

则，该网格的胖射线路径值为有效值，记为预先设定的值；否则为无效值，记为0；从而得到这一对发射阵元-接收阵元下、对应于成像区域各个网格的一系列胖射线路径值，将这一系列胖射线路径值行向量化，得到胖射线路径行向量；

所述公式(2)中，S代表发射阵元，R代表接收阵元，S和R互为发射接收，τ_SP代表声波由S出发、到达P的时间，τ_RP代表声波由R出发、到达P的时间，τ_SR代表声波由S出发、到达R的时间；Δt为按迭代次数从预先设定的时间阈值数列{Δt}中选取的第i个时间阈值，i与所述迭代次数在数值上相等，该时间阈值数列{Δt}中的时间阈值是按由大到小的顺序排列，且该时间阈值数列{Δt}中的任意一个时间阈值均与探头中心频率f相关；

从而得到针对同一个发射阵元、且按接收阵元顺序对应的一系列胖射线路径行向量，并最终得到针对全部发射阵元的、且按发射阵元及接收阵元顺序对应的胖射线路径行向量总集合；将每个行向量视为一个元素整体，将该胖射线路径行向量总集合中的所有元素排列成一列，即可得到行位置与发射阵元及接收阵元顺序相对应的胖射线路径总矩阵[L]；

反问题求解模块用于：基于所述能量参数列向量[P]和所述胖射线路径总矩阵[L]构建如式(4)所示的路径-衰减-能量参数方程组：

[L][A]＝[P] (4)

该式(4)中，[A]是总行数为m₁×m₂、总列数为1的待求解的衰减系数列向量；

采用随机梯度下降法解方程组求解得到[A]，从而得到针对待测对象的衰减系数重建值，该衰减系数重建值即为更新后的衰减系数值；

判断是否终止的功能模块用于：判断迭代次数是否小于预先设定的迭代次数阈值，若小于则重复迭代；若不小于，则终止迭代，并输出最终的衰减系数值；

此外，所述预先设定的时间阈值数列{Δt}中时间阈值的总个数大于等于所述预先设定的迭代次数阈值。

通过本发明所构思的以上技术方案，与现有技术相比，由于基于费马原理，不仅考虑射线路径上的点对接收信息的影响，还考虑射线领域上的其它点对接收信息的影响，得到胖射线路径总矩阵，路径的追踪过程简单快捷，且成像分辨能力更高。另一方面，尤为重要的是，本发明基于变参数的、路径由宽至窄的变路径宽度的胖射线路径，实现了低分辨率到高分辨率加深重建的方法，使重建过程更加稳定，重建分辨率更高，减少重建伪影。

针对某一对发射阵元-接收阵元，判断声波从该发射阵元到该接收阵元经过的胖射线路径时，具体是基于公式(2)，当成像区域的成像点满足该公式(2)时，则判断该成像点位于胖射线路径上，该成像点的胖射线路径值被标记为有效值，否则记为无效值0。并且，该有效值可预先统一设定为1，这样针对全部成像点将得到非0即1的胖射线路径矩阵，便于后续计算处理。

τ_SP+τ_RP-τ_SR≤Δt (2)

每次迭代时，式(2)中Δt的取值将发生变化，迭代次数越大，Δt越小，也就是说，Δt在迭代的初期是比较大的值，随着迭代加深，Δt逐渐减小；例如，Δt可以从(1/f)、…、(1/3f)、…、(1/5f)、…、(1/10f)、…依次降低。Δt越大，胖射线路径越宽，相反地，Δt越小，胖射线路径越窄。由于对于不同成像对象、不同成像系统，路径的宽窄直接影响成像的分辨率，本发明中这种Δt随着迭代的加深逐渐减小的变化过程，实质对应了一种变参数(变路径宽度)的胖射线路径的图像重建方法，是一种从低分辨率到高分辨率加深重建的方法，使重建过程更加稳定，重建分辨率更高，减少重建伪影。并且，本发明所采用的从宽路径到窄路径的迭代，有利于避免陷入局部最优解；本发明预先设定的迭代次数阈值可以是大于等于10的正整数(如50次、100次等)；若采用一次迭代并直接取Δt最小，那么得到的将是局部最优解，而本发明通过缓慢变化能够避免陷入局部最优解。

进一步的，本发明通过选取满足预先设定接收阵元角度要求范围内全部接收阵元接收到的透射波数据，能进一步对透射波数据的准确性进行把控。考虑到发射阵元附近的接收阵元主要接收的是反射信号，针对任意一个发射阵元，本发明优选采集以发射阵元正对面为基准、以所述环形阵列探头的圆心为圆心、圆心角角度从-α到+α共计2α角度范围内的接收阵元接收得到的透射波数据，2α优选满足90°～270°(即角度范围最少不少于90°，最多不超过270°)，从而再利用发射阵元、接收阵元预先设定的顺序(如逆时针顺序或顺时针顺序)，得到元素行位置与发射阵元及接收阵元顺序相对应的渡越时间列向量[T_tof](或能量参数列向量[P])及胖射线路径总矩阵[L]。

此外，本发明采用随机梯度下降法来求解该方程组，随机梯度下降法作为一种迭代求解方程组的方法，介于梯度下降法和牛顿法之间，是一种良好的迭代优化方法。本发明中，由于超声CT重建数据量极大，计算量极大，采用随机梯度下降法求解方程组，在迭代优化求解过程中，每次迭代随机选取一个样本来对待重建参数进行更新，使得优化过程加快，从而能大大减小运算负担。这里所说的一个样本，即为反问题求解步骤中公式(3)或公式(4)方程组中的某一个方程。

总体来说，本发明提出的基于费马原理的声速及衰减系数重建方法及系统的优点有以下几点：1、基于费马原理，路径的追踪过程简单快捷；2、基于费马原理的路径追踪，使成像分辨能力更高；3、基于变路径的胖射线路径的图像重建，能够实现从低分辨率到高分辨率加深重建，使重建过程更加稳定，重建分辨率更高，减少重建伪影；4、判断路径是否有效时，优选采用非0即1，给后续计算带来方便；5、采用随机梯度下降法求解反问题，避免了海森矩阵的计算，运算量小，求解精度高。

附图说明

图1是AIC法提取渡越时间示意图。

图2是本发明中基于费马原理的胖射线路径示意图。

图3是本发明所采用的仿真数据真实模型示意图。

图4是本发明方法重建出的声速图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。

本发明中基于费马原理的超声CT图像重建方法及系统，不仅考虑到射线路径上的点对该信息具有影响，并且射线领域上的其它点对接收信息也具有影响。同时，胖射线路径采用变参数(变路径宽度)思想，路径宽度随着迭代的加深(即，随着迭代次数的增多)，不断收窄，能够实现从低分辨率到高分辨率加深重建。

以声速重建为例，本发明中基于费马原理的声速重建方法，概括来说，包括以下步骤：

(1)提取渡越时间。

(2)给定初始声速值，如水的声速值，进入第一次迭代。

(3)根据有限差分方法计算从每个发射阵元出发到成像区域每个像素点的时间τ，根据费马原理计算胖射线路径。

(4)反问题的求解，更新声速值。

(5)判断迭代次数是否小于给定的迭代次数阈值，若小于，回到步骤(3)，否则终止迭代，输出声速值。

具体的：

超声CT系统采用环形阵列探头，换能器阵元呈环形均匀分布。采集数据利用单发全接的模式，即一个阵元发射，所有阵元接收，依次使所有阵元发射，假设有T个阵元，则可采集到T²条数据记录。

首先对数据进行渡越时间(tof)的提取。渡越时间即波形到达接收位置的时间，该时间是接收波形中首达波形的起跳点时间。本发明采用Akaike信息准则(AIC)提取。

接着计算声波从发射阵元到接收阵元经过的射线路径。本发明基于费马原理，因此路径不是射线，而是一个区域。首先根据费马原理，判断成像点是否在路径区域，并存储在路径矩阵当中。最后是建立并求解路径-慢度-时间方程组，即反问题的建立与求解。慢度是声速的倒数，本发明中的慢度指待重建对象的慢度，路径矩阵可以是非零即一的矩阵，若该点在路径上，即为1，否则为0(1表示成像点在路径区域，0表示成像点不在路径区域)。时间是待重建对象的渡越时间。

并且，在路径处理过程中，每次迭代时，随着迭代次数的不同，胖射线路径宽窄不同；迭代次数越大，路径越窄，由此实现变路径(变路径宽度)的胖射线路径，从而对应一种从低分辨率到高分辨率加深重建的方法，使重建过程更加稳定，重建分辨率更高，减少重建伪影。

以下为具体实施例：

实施例1

以体模为待测对象，本发明中基于费马原理的声速重建方法是将超声CT系统采集到的透射数据进行适当的处理，最后重建得到超声声速图像；具体处理步骤包括渡越时间的提取、迭代开始的准备、迭代过程中胖射线路径的计算、反问题的求解、迭代是否终止的判断、以及进一步的声速图像的显示。

本发明采集数据采用环形阵列探头，利用单发全接的模式，即，一个阵元发射，所有阵元接收，依次使所有阵元发射，假设有T个阵元，则采集到T²条数据记录。透射波主要由发射阵元位置对面的阵元接收，因此本发明采用发射阵元对面的180度范围内阵元(即，发射阵元正对面±90度、共计180度范围内、共计T/2个阵元；当然也可以采用其他预先设定的接收阵元角度要求范围，如其他满足90°～270°的具体角度范围要求，以90°的角度要求范围为例，发射阵元正对面±45度、共计90度范围内的阵元)接收到的数据重建声速，即，从T²条数据记录中选取T×T/2条数据。

具体处理时，首先是采用AIC法提取待重建对象数据的渡越时间。

AIC法的算法流程如下：首先在数据上预估的渡越时间点附近选取合适的窗，窗长度为N。以当前检索点为分界点，将窗分为两段，利用公式(1)计算当前检索点的AIC值。var代表方差，d(1，k)表示第一个点到第k个点的数据。对窗内的点依次检索，记AIC值最小的点为渡越时间点。

AIC(k)=k*log(var(d(1,k)))+(N-k-1)*log(var(d(k+1，N))) (1)

这样针对T×T/2条数据，就可以得到T×T/2个渡越时间，将这些渡越时间排列成一个列向量，即可得到总行数为(T×T/2)、总列数为1的渡越时间列向量。

接着，准备进入迭代：给定预先设定的声速值作为初始声速值，如水中的声速值，并记迭代次数为1，进入第一次迭代。

接着，根据有限差分方法计算从每个发射阵元出发到成像区域每个像素点的时间τ，并根据费马原理计算胖射线路径，即，计算声波从发射阵元到接收阵元经过的胖射线路径。本发明基于费马原理，即认为不仅射线路径上的点对该信息具有影响，射线领域上的其它点对接收信息也具有影响。路径不是一条射线，而是一个区域。该区域的判断基于费马原理。如图2所示，是胖射线路径的示意图，S为发射阵元，R为接收阵元，S和R互为发射接收。P为声场中的任意一点，首先计算S发射声波，到达P点的时间，接着计算R发射，到达P点的时间。若

τ_SP+τ_RP-τ_SR≤Δt (2)

则P点在路径区域内(式中，Δt为一个依据经验预先选取的时间阈值)。同时，每次迭代时，Δt的取值不同，实现变路径(变路径宽度)的胖射线路径。Δt在迭代的初期是比较大的值，随着迭代加深，Δt逐渐减小。例如可以从(1/f)、…、(1/3f)、…、(1/5f)、...(1/10f)、…依次降低，f指探头中心频率。Δt越大，胖射线路径越宽，相反地，Δt越小，胖射线路径越窄。具体实施时，可以如下：Δt为按迭代次数从预先设定的时间阈值数列{Δt}中选取的第i个时间阈值，i与所述迭代次数在数值上相等，该时间阈值数列{Δt}中的时间阈值是按由大到小的顺序排列。本发明采取的是一种变参数(变路径宽度)的胖射线路径的图像重建方法，是一种从低分辨率到高分辨率加深重建的方法，使重建过程更加稳定，重建分辨率更高，减少重建伪影。从宽路径到窄路径的迭代，有利于避免陷入局部最优解)。由于对于不同成像对象、不同成像系统，路径的宽窄直接影响成像的分辨率，因此{Δt}可以靠经验预先选取。

将成像区域划分成m×m个网格(m为正整数，数值可以预先设定；当然，成像区域长度方向上划分的网格数与其宽度方向上划分的网格数两者可以不同，例如，也可以将成像区域划分成m₁×m₂个网格，m₁、m₂不相等；与现有技术常规设置相似，m₁×m₂优选与渡越时间总个数相等，允许一定偏差)，每个网格是否属于满足费马原理、属于路径区域，则可以将网格的中心点基于公式(2)来判断。判断属于路径区域的话，则可以将该网格的值记为1；判断不属于路径区域的话，则可以将该网格的值记为0；由此得到一个元素值非0即1的总列数为(m×m)的行向量。这样针对所有发射阵元-接收阵元，可以得到总行数为T×T/2、总列数为m×m的路径矩阵L，由此完成路径计算。

在路径计算完成之后，构建路径-慢度-时间方程组，S是慢度(S为总行数满足(m×m)的列向量)，T_tof是渡越时间列向量(该列向量的总行数满足(T×T/2))，L是路径矩阵。

LS＝T_tof (3)

采用随机梯度下降法解方程组，输出S，取倒数，则为此次迭代得到的体模的速度重建值，从而完成第1次迭代。

此次迭代完毕后，将得到的速度重建值作为更新后的声速值，同时更新迭代次数，使迭代次数加1，重复上述迭代过程的具体操作进行再次迭代，直到迭代次数大于等于预先设定的迭代次数阈值，此时再终止迭代。

其中，最初预估的渡越时间点，可参考现有技术中的常规操作，例如，可以根据水的理论声速来计算大概的渡越时间点，或者利用肉眼在波形上观察渡越时间点所在的大概位置。窗长度N的选择也是根据波形的实际形状，经验选取，如果不合适可以适当增大或减小。

实施例2

基于与实施例1中重建声速相似的处理原理及处理步骤，可以重建衰减系数。

以同样的环形阵列探头为例，首先，利用单发全接的模式，采集得到T²条能量参数数据记录。透射波主要由发射阵元位置对面的阵元接收，因此本发明采用发射阵元对面的180度范围内阵元(即，发射阵元正对面±90度、共计180度范围内、共计T/2个阵元)接收到的能量参数数据重建衰减系数，即，从T²条数据记录中选取T×T/2条能量参数数据，得到能量参数的列向量P，P的总行数满足(T×T/2)。

然后，与实施例1中声速重建相似，给定预先设定的衰减系数值作为初始衰减系数值，记迭代次数为1，进入第一次迭代；迭代过程中，按迭代次数选择预先设定的Δt时间阈值参与构建公式(2)，并得到胖射线路径矩阵L。与实施例1相似，实施例2中也是每次迭代时，Δt的取值不同，实现变路径(变路径宽度)的胖射线路径。Δt在迭代的初期是比较大的值，随着迭代加深，Δt逐渐减小。

接着，构建如公式(4)所示的路径-衰减-能量参数方程组：

LA＝P (4)

其中，A为待求解的衰减系数；然后，采用随机梯度下降法解方程组，求解得到一维含有(m×m)个元素的向量A；这样即可直接得到待测对象的衰减系数重建值，从而完成第1次迭代。

此次迭代完毕后，将得到的衰减系数重建值作为更新后的衰减系数值，同时更新迭代次数，使迭代次数加1，重复上述迭代过程操作进行再次迭代，直到迭代次数大于等于预先设定的迭代次数阈值，此时再终止迭代。

本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。