基于改进射线追踪法的超声ct成像系统
阅读说明：本技术 基于改进射线追踪法的超声ct成像系统 (Ultrasonic CT imaging system based on improved ray tracing method ) 是由 赵悦 袁于 沈毅 王艳 冯乃章 于 2019-10-23 设计创作，主要内容包括：本发明涉及一种基于改进射线追踪法的超声CT成像系统,包括全包围超声换能器、信号采集模块、数据处理及成像模块。数据处理及成像模块中的处理器执行储存在存储器中的计算机程序时实施步骤：基于全包围超声换能器扫查人体器官,测量声波的渡越时间；重建目标初始的声速分布；基于目标的声速分布进行B超成像；从B超图像中提取器官组织边界信息；迭代更新声波的传播路径；迭代更新目标的声速分布。本发明与传统成像方案相比,在保证超声CT成像精度的同时,减小计算量,提高重建速度,并改善B超图像的精度和超声CT成像结果中的边缘清晰度。(The invention relates to an ultrasonic CT imaging system based on an improved ray tracing method, which comprises a full-surrounding ultrasonic transducer, a signal acquisition module and a data processing and imaging module. The processor in the data processing and imaging module, when executing the computer program stored in the memory, performs the steps of: scanning human organs based on a fully-enclosed ultrasonic transducer, and measuring the transit time of sound waves; reconstructing initial sound velocity distribution of a target; b ultrasonic imaging is carried out based on the sound velocity distribution of the target; extracting organ tissue boundary information from the B ultrasonic image; iteratively updating the propagation path of the sound wave; the sound speed distribution of the target is iteratively updated. Compared with the traditional imaging scheme, the invention reduces the calculated amount, improves the reconstruction speed and improves the precision of the B-mode ultrasonic image and the edge definition in the ultrasonic CT imaging result while ensuring the ultrasonic CT imaging precision.)
技术领域
本发明涉及一种超声CT成像系统，具体涉及一种改进的射线追踪法，以及基于该算法的利用B超中提取的乳腺组织边界信息进行超声CT成像的系统。
背景技术
超声CT作为一种新兴的超声断层图像重建方法，其基本原理是利用超声波对***进行扫查，并从超声信号中提取渡越时间等有效信息，以对乳腺组织的声学性质如声速分布等进行重建。由于其具有分辨率高、无辐射等特点，在无损探伤、乳腺癌早期筛查等领域有着广泛的应用前景。
传统的超声CT成像过程可以被分为三个步骤：利用全包围探头对***进行扫查并从声信号中提取渡越时间；重建声在乳腺组织中的传播路径；基于渡越时间和声传播路径，求解成像区域内的声速分布。
上述步骤中，重建准确的声传播路径是一大难点。根据重建声路径的方法不同，超声CT的重建方法又可以分为两类。第一类方法中，超声波沿直线传播，其传播路径是从发射阵元到接收阵元的直线。该方法计算简单，但是没有考虑到声在乳腺组织中的折射、反射等现象，因此基于该方法的重建结果中有着边界模糊、形状畸变等问题。第二类方法对声的传播过程进行了更精确的建模，因此超声波沿曲线传播。该方法中，首先基于惠更斯原理计算声波到达成像区域内各点的渡越时间，然后利用梯度下降法找到渡越时间最短的路线作为声路径，这种方法被称为射线追踪法。与第一种方法相比，射线追踪法考虑了声在传播过程中的折射和衍射现象，其数学模型与真实模型更加接近，因此重建结果更为准确；其缺点是声传播路径需要在迭代过程中逐步优化，所需的计算量极为庞大，这降低了该方法的实用性。
与传统的超声成像(即B超成像)方式相比，超声CT的重建结果能直观地反应乳腺组织的声速分布。但是受重建原理所限，该成像方法对边界不敏感，重建结果中常存在着边界模糊和计算量大的问题。反之，B超成像利用回波对乳腺组织的边界进行检测，重建耗时短，边界明显，但是重建结果中斑点噪声较为严重。除此之外，在B超成像中存在着若干假设，如组织声速为1540m/s、超声波沿直线传播等，而实际的重建过程中这些假设并不严格成立，因此在重建结果中往往有着严重的形状畸变。
本发明提出了一种技术方案，其基于改进的射线追踪法，利用从B超中提取的乳腺组织边界信息计算超声波在***中的传播路径。该方案求得的超声波传播路径与传统的射线追踪法的结果接近，但是计算量减少了90％以上。基于该算法进行超声CT的重建，在保证重建结果高精度的同时，可以将重建时间减小一个量级以上；除此之外，本发明还能同时改善B超图像的精度和超声CT成像结果中的边缘清晰度。
发明内容
本发明的目的在于改进传统的射线追踪法，并基于该算法将B超图像中提取的乳腺组织边界信息用于超声CT成像。与基于传统射线追踪法的超声CT成像相比，在保证超声CT成像精度的同时，减小计算量，提高重建速度，并改善B超图像的精度和超声CT成像结果中的边缘清晰度。
本发明的技术方案第一方面为一种超声CT成像系统，包括：配置成环绕人体器官的已确定位置的全包围超声换能器；
信号采集模块，用于对全包围超声换能器接收到的声信号进行采样；
与所述信号采集模块连接的数据处理及成像模块，该数据处理及成像模块包括存储器和处理器，其特征在于，所述处理器执行储存在所述存储器中的计算机程序时实施的方法包括：
步骤一、基于全包围超声换能器扫查人体器官，测量声波的渡越时间；
步骤二、重建目标初始的声速分布；
步骤三、基于目标的声速分布进行B超成像；
步骤四、从B超图像中提取器官组织边界信息；
步骤五、迭代更新声波的传播路径；
步骤六、迭代更新目标的声速分布。
在一些实施例中，所述的步骤一包括：使全包围超声换能器发射和接收超声波，并基于“阈值-相位”法从记录的超声信号中提取渡越时间。
在一些实施例中，所述的步骤二包括：记声在器官组织中传播的路径为从发射阵元到接收阵元的直线，以L(0)表示；记声慢度为s；并且基于SART算法求解公式L(0)s(0)＝t，得到s的初始解s(0)。
在一些实施例中，所述的步骤三包括：基于全包围超声换能器，通过电子聚焦和延时求和的方式进行数字波束合成，对人体器官进行B超成像；在延时求和的过程中，基于步骤二或步骤五中重建得到的器官组织的声慢度分布计算各个发射阵元的发射延迟时间。
在一些实施例中，所述的步骤四包括：利用Sobel算子对B超成像结果中的每个像素点做卷积和运算，然后选取阈值以提取边缘；进行滤波操作后，生成器官组织的边界。
在一些实施例中，所述的步骤五包括：对声传播过程中在器官组织边界上发生折射和反射现象进行建模，基于声慢度s(r)计算超声波到达器官组织边界每个点的声达时间，并迭代更新超声波的在器官组织中的传播路径，得到L(r+1)。
在一些实施例中，所述的步骤六包括：迭代更新目标的声速分布，基于SART算法求解L(r+1)s(r+1)＝t，求得更新后的声慢度s(r+1)。
在一些实施例中，所述的全包围超声换能器包括由多个阵元组成的环状超声探头。所述的系统还包括升降装置，其中升降装置包括：支撑架，用于支撑整个全包围超声换能器及附属的扫查设备；内空腔，其中填充超声耦合剂以接纳人体器官；以及电机，用于带动全包围超声换能器沿人体器官的参考平面的法线方向进行线性运动。所述的系统，还包括控制模块，该控制模块一方面通过电子开关切换来控制全包围超声换能器的声波发射和接收，另一方面通过控制升降装置来调节换能器的高度，以实现对人体器官的三维扫查
本发明的技术方案第二方面为一种计算机可读存储介质，其上储存有计算机程序。所述计算机程序被处理器执行时实施的方法包括：
步骤一、基于全包围超声换能器扫查人体器官，测量声波的渡越时间；
步骤二、重建目标初始的声速分布；
步骤三、基于目标的声速分布进行B超成像；
步骤四、从B超图像中提取器官组织边界信息；
步骤五、迭代更新声波的传播路径；
步骤六、迭代更新目标的声速分布。
本发明与现有技术相比具有如下优点：
本发明改进了传统的射线追踪法，基于从B超图像中提取的乳腺组织边界重建声波的传播路径与声慢度分布，与传统的射线追踪法相比，在保证超声CT成像精度的同时，计算量减少了90％以上；同时，由于将B超图像中提取的乳腺组织边界应用于超声CT中，该方案的成像结果边界更为清晰；除此之外，基于超声CT重建的声慢度分布进行电子聚焦和波束合成，改善了B超的成像结果。
本发明提出的成像方案适用性强，成像效果好。该方案虽是基于全包围超声换能器提出，但其并非是必需的，只要使用位置精确可知的超声换能器对人体***采集足够多的透射式投影数据，例如使用若干线型超声换能器环绕人体***旋转并采集超声信号，即可根据本发明提出的重建方法进行超声CT成像；该方案同时改善了超声CT和B超的成像结果，其中，超声CT的成像结果与传统方案相比，在保证重建精度不下降的基础上，计算量减少了90％以上。
附图说明
图1为本发明的系统架构框图。
图2为本发明系统的实施方法的流程图。
图3为本发明系统中的全包围超声换能器进行超声CT成像的示意图。
图4为本发明中从声信号中提取渡越时间的示意图。
图5为本发明中声波沿直线或折线传播的示意图。
图6为本发明中重建乳腺组织声速分布使用的仿真数字标准体模。
图7为本发明中得到的乳腺组织声速分布初始值的示意图。
图8为本发明中最终重建得到的乳腺组织声速分布的示意图。
具体实施方式
下面结合实施例说明本发明的具体实施方式。本发明的超声CT成像系统可以应用于人体器官或者组织的成像，本文主要以***为例描述成像方案，但不作为成像对象的限制。
根据本发明的超声CT成像系统，其硬件架构可以包括：升降装置、全包围超声换能器、信号采集模块、控制模块、数据处理及成像模块，如图1所示。其中，升降装置包括支撑架、内空腔以及电机。全包围超声换能器位置可以包括由多个阵元组成的环状超声探头，且其空间位置预先确定。多个阵元可以以圆形的全包围方式布置，如图3和图5所示。支撑架用于支撑整个扫查设备(包括全包围超声换能器)，内空腔用于接纳人体***并盛放超声耦合剂(例如水)，电机与全包围超声换能器相连，用于控制换能器在竖直平面内的升降运动。全包围超声换能器环绕升降装置的内空腔，用于发射并接收超声波信号。数据采集模块对全包围超声换能器接收到的声信号进行采样，并将结果输出到数据处理及成像模块。控制模块一方面通过电子开关切换来控制全包围超声换能器的声波发射和接收，另一方面通过控制升降装置来调节换能器的高度，以实现对人体***的三维扫查。数据处理及成像模块负责对采集到的超声信号进行分析，并完成超声CT图像及B超图像的重建。
参照图2，数据处理及成像模块实施超声CT成像的方法主要包括以下步骤：首先基于全包围超声换能器扫查人体***，测量声波的渡越时间；其次重建目标初始的声速分布；第三，基于目标的声速分布进行B超成像；第四，从B超图像中提取乳腺组织边界信息；第五，迭代更新声波的传播路径；最后，迭代更新目标的声速分布。下面对各个步骤及运算过程进行详细描述。
步骤一：基于全包围超声换能器扫查人体***，测量声波的渡越时间
本步骤中，用全包围超声换能器环绕人体***，发射并接收超声波，并基于“阈值-相位”法从记录的超声信号中提取渡越时间。
具体步骤如下：
(1.1)在人体***置于阵元数为N的全包围超声换能器中央，且换能器与***浸没于水中后，使换能器各阵元依次发射声波；每个阵元发射声波时，所有阵元接收并记录声信号，由此记录总计N×N个声信号；
(1.2)设接收到的声信号以p＝f(τ)表示，其中p表示声压，τ表示声的传播时间；则单个超声信号对应的渡越时间t由下式得到：
其中，Δp是阈值，其值略高于记录的超声信号中的噪声最大值，τ-和τ+分别表示τ的左邻域和右邻域；该式的含义为：在接收到的超声信号第一次高于阈值后，取第一个极大值点作为特征点P，P所对应的时间即为该超声信号的声达时间t；从每个超声信号中提取渡越时间，构成一个N2×1的向量t：
其中，th表示从第h个发射-接收阵元对对应的声信号中提取的渡越时间。
步骤二：重建目标初始的声速分布
记声在乳腺组织中传播的路径为从发射阵元到接收阵元的直线，以L(0)表示；记声慢度(即声速度的倒数)为s；基于SART算法求解公式L(0)s(0)＝t，得到s的初始解s(0)。
具体步骤如下：
(2.1)将成像区域划分为M×M个网格，以直线连接发射阵元与接收阵元，并求出该直线在每个网格中的长度，由于共有N2个发射-接收阵元对，因此结果为一个N2×M2维的矩阵，记作
(2.2)重建成像区域内的声速分布：首先将声速sinit初始化为水的声速，然后基于SART算法求解L(0)s(0)＝t：
式中，λ是用于调控收敛快慢的参数；上式的解的形式为：
其中，表示第k个网格中的声慢度的初始值。
步骤三：基于目标的声速分布进行B超成像
该步骤中，基于全包围超声换能器，通过电子聚焦和延时求和的方式进行数字波束合成，对人体***进行B超成像；延时求和的过程中，基于步骤二或步骤五中重建得到的乳腺组织的声慢度分布计算各阵元的发射延迟时间：
式中，Tr表示发射超声波的阵元，Fo表示电子聚焦的焦点，
步骤四：从B超图像中提取乳腺组织边界信息
该步骤中，利用Sobel算子对B超成像结果中的每个像素点做卷积和运算，然后选取合适的阈值以提取边缘，进行简单的滤波操作后，即可得到乳腺组织的边界。
步骤五：迭代更新声波的传播路径
对声传播过程中在乳腺组织边界上发生折射和反射现象进行建模，基于声慢度s(r)计算超声波到达乳腺组织边界每个点的声达时间，并迭代更新超声波的在乳腺组织中的传播路径，得到L(r+1)。
对每个发射超声波的阵元Tr，接收超声波的阵元为Tr′，实施如下步骤：
(5.1)对乳腺组织边界上的每个点PA，计算超声波从Tr到达PA的声达时间的初始值，具体计算方式为：若Tr与PA之间的连线经过了乳腺组织边界，则
(5.2)以下式迭代更新从Tr到达乳腺组织边界上每个点PA的声达时间
式中，表示更新后的声达时间，PB表示乳腺组织边界上除PA外的任一点，
(5.3)计算从Tr到Tr′的最小的声传播时间μTr,Tr′：
该声传播时间对应的声路径即为更新后的声传播路径。
对每个发射-接收阵元对之间的声传播路径进行迭代更新，更新后的声波传播路径矩阵记为L(r+1)。
步骤六：迭代更新目标的声速分布
基于SART算法求解L(r+1)s(r+1)＝t，求得更新后的声慢度s(r+1)：
其中Δs为预先设定的小参量，若s(r+1)满足迭代停止条件|s(r+1)-s(r)|<Δs，则s(r +1)为最终输出的声慢度分布，否则将其作为先验信息输入步骤三，再次进行B超成像，直到求得的s(r+1)满足条件。
下面，结合具体的计算实例，通过基于MATLAB仿真平台，对本发明中超声CT成像系统的技术方案进行仿真验证及评价。实例中设计并使用的数字体模如图6所示。
执行步骤一：基于全包围超声换能器扫查人体***，测量声波的渡越时间
用全包围环形超声换能器对体模进行扫查并记录超声信号，再基于“阈值-相位”法提取渡越时间。
具体步骤如下：
(Z1.1)用阵元数为256的全包围的环形超声换能器对仿真数字体模进行扫查(参考图3或5)：换能器各阵元依次发射声波，每个阵元发射声波时，所有阵元接收并记录声信号，由此记录总计65536个声信号；
(Z1.2)接收到的声信号p＝f(τ)的波形如图4所示，其中p表示声压，在本实施例中p经过了归一化处理，即信号幅值绝对值的最大值为1，τ表示声的传播时间；由下式计算该超声信号对应的渡越时间t：
其中，Δp是阈值，其值略高于记录的超声信号中的噪声最大值，本实施例中取0.07，τ-和τ+分别表示τ的左邻域和右邻域。由上式计算得该超声信号对应的渡越时间为1.825微秒；从每个超声信号中提取渡越时间，构成一个65536×1的向量t。
执行步骤二：重建目标初始的声速分布
记声在乳腺组织中传播的路径为从发射阵元到接收阵元的直线，以L(0)表示；记声慢度(即声速度的倒数)为s；基于SART算法求解公式L(0)s(0)＝t，得到s的初始解s(0)。
具体步骤如下：
(Z2.1)将成像区域划分为256×256个网格，以直线连接发射阵元与接收阵元，并求出该直线在每个网格中的长度，由于共有65536个发射-接收阵元对，因此结果为一个65536×65536维的矩阵，记作
(Z2.2)重建成像区域内的声速分布：首先将声速sinit初始化为水的声速1515m/s，然后基于SART算法求解L(0)s(0)＝t：
式中，λ是用于调控收敛快慢的参数，上式的解是一个65536×1的向量，每个元素分别表示65536个网格内的声慢度值，再根据其重建出成像区域内的速度值分布，结果如图7所示。
执行步骤三：基于目标的声速分布进行B超成像
该步骤中，基于全包围超声换能器，通过电子聚焦和延时求和的方式进行数字波束合成，对人体***进行B超成像；延时求和的过程中，基于步骤二或步骤五中重建得到的乳腺组织的声慢度分布计算各阵元的发射延迟时间：
式中，Tr表示发射超声波的阵元，Fo表示电子聚焦的焦点，
执行步骤四：从B超图像中提取乳腺组织边界信息
该步骤中，利用Sobel算子对B超成像结果中的每个像素点做卷积和运算，然后选取合适的阈值以提取边缘，进行简单的滤波操作后，即可得到乳腺组织的边界。
执行步骤五：迭代更新声波的传播路径
对声传播过程中在乳腺组织边界上发生折射和反射现象进行建模，基于声慢度s(r)计算超声波到达乳腺组织边界每个点的声达时间，并迭代更新超声波的在乳腺组织中的传播路径，得到L(r+1)。
对每个发射超声波的阵元Tr，接收超声波的阵元为Tr′，实施如下步骤：
(Z5.1)对乳腺组织边界上的每个点PA，计算超声波从Tr到达PA的声达时间的初始值，具体计算方式为：若Tr与PA之间的连线经过了乳腺组织边界，则
(Z5.2)以下式迭代更新从Tr到达乳腺组织边界上每个点PA的声达时间
式中，
(Z5.3)计算从Tr到Tr′的最小的声传播时间μTr,Tr′：
该声传播时间对应的声路径即为更新后的声传播路径。
对每个发射-接收阵元对之间的声传播路径进行迭代更新，更新后的声波传播路径矩阵记为L(r+1)。
执行步骤六：迭代更新目标的声速分布
基于SART算法求解L(r+1)s(r+1)＝t，求得更新后的声慢度s(r+1)：
现有Δs＝2×10-6(m/s)-1为预先设定的小参量，若s(r+1)满足迭代停止条件|s(r+1)-s(r)|<Δs，则s(r+1)为最终输出的声慢度分布，否则将其作为先验信息输入执行步骤三，再次进行B超成像，直到求得的s(r+1)满足条件，本实施例中，迭代一次后s(r+1)收敛。
本实施例中，基于改进的射线追踪法重建的乳腺组织声速分布如图8所示，其重建结果与基于传统的射线追踪法的结果相比，重建图像的质量基本相当：以重建结果的声速分布与仿真体模的声速分布的均方根误差作为评价指标，两种方法的结果相差不超过2.5％；特别地，基于本发明提出的算法的重建结果比传统的射线追踪法的边界清晰度高出50％以上；除此之外，基于改进的射线追踪法重建乳腺组织的声速分布，与传统的射线追踪法相比，其计算量减少了90％以上，可以节省重建时间，在本实施例中，基于改进的射线追踪法重建乳腺组织的声速分布耗时952秒，而传统的射线追踪法耗时为12018秒。
应当认识到，本发明的实施例还可以由计算机硬件、硬件和软件的组合、或者通过存储在非暂时性计算机可读存储器中的计算机指令来实现或实施。所述方法可以使用标准编程技术-包括配置有计算机程序的非暂时性计算机可读存储介质在计算机程序中实现，其中如此配置的存储介质使得计算机以特定和预定义的方式操作-根据在具体实施例中描述的方法和附图。每个程序可以以高级过程或面向对象的编程语言来实现以与计算机系统通信。然而，若需要，该程序可以以汇编或机器语言实现。在任何情况下，该语言可以是编译或解释的语言。此外，为此目的该程序能够在编程的专用集成电路上运行。
进一步，上述方法可以在可操作地连接至合适的任何类型的计算平台中实现，包括但不限于个人电脑、迷你计算机、主框架、工作站、网络或分布式计算环境、单独的或集成的计算机平台、或者与带电粒子工具或其它成像装置通信等等。本发明的各方面可以以存储在非暂时性存储介质或设备上的机器可读代码来实现，无论是可移动的还是集成至计算平台，如硬盘、光学读取和/或写入存储介质、RAM、ROM等，使得其可由可编程计算机读取，当存储介质或设备由计算机读取时可用于配置和操作计算机以执行在此所描述的过程。此外，机器可读代码，或其部分可以通过有线或无线网络传输。当此类媒体包括结合微处理器或其他数据处理器实现上文所述步骤的指令或程序时，本文所述的发明包括这些和其他不同类型的非暂时性计算机可读存储介质。当根据本发明所述的方法和技术编程时，本发明还包括计算机本身。
以上所述，只是本发明的较佳实施例而已，本发明并不局限于上述实施方式，只要其以相同的手段达到本发明的技术效果，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明保护的范围之内。在本发明的保护范围内其技术方案和/或实施方式可以有各种不同的修改和变化。
- 上一篇：一种医用注射器针头装配设备
- 下一篇：超声探头