CN105997153A - 一种超声ct的成像方法 - Google Patents

Abstract

本发明公开了一种超声CT的成像方法。包括将成像区域网格化为M个网格点；环形探头以成像区域为中心360°旋转，获得Q*N*K组数据；然后获得该数据对应的渡越时间向量T₀以及直线路径向量L₀；最后迭代计算，获得所述M个网格点对应的速度分布V，并将所述速度分布V归一化，完成成像。本发明通过优化环形探头中发射阵元以及接收阵元的采集次数和采集角度，并对获得的数据进行处理与图像重建，从而使重建出来的速度分布更精确，以提高成像的信噪比，提高系统的成像质量。

Description

一种超声CT的成像方法

技术领域

本发明属于超声断层成像领域，更具体地，涉及一种超声CT的成像方法。

背景技术

超声CT成像，即超声断层成像，具有超声非侵入、无辐射和成本低的优势，以断层成像的方式进行扫描，在乳腺癌等疾病早期诊断方面具有广阔的应用前景。

超声CT成像分为反射式成像和透射式成像，其中，采用声速的透射式成像是通过重建乳腺内组织的声速分布和衰减系数进行成像。超声CT声速成像利用首达波到达时间，重建乳腺内组织的声速分布，从而正确地区分良性肿瘤与癌。因为和正常组织与良性肿瘤相比，癌具有较高的声速。

目前由在乳腺癌早期检测领域走在前列的美国Karmanos癌症中心研制的环形阵列的超声CT成像扫描器CURE可以得到亚毫米级高分辨率的图像。这种扫描器每次使用单阵元向环形阵列的中心发射超声波，所有阵元接收，直到每个阵元循环发射一次，所以可以接收到多个角度的超声信号，得到斑点噪声少的高分辨率图像。虽然这种扫描器实现了对成像物体的360度扫描，但由于环形探头之间存在间隙，从而无法获取到更多的成像目标的信息，进而声速重建时由于方程组的解并非唯一，所以求解出来的声速分布跟实际的声速分布差别很大，重建图像的质量不高。

发明内容

针对现有技术的以上缺陷或改进需求，本发明提供了一种超声CT的成像方法，其目的在于优化环形探头中发射阵元以及接收阵元的采集次数和采集角度，并对获得的数据进行处理与图像重建，从而使重建出来的速度分布更精确，以提高成像的信噪比，提高了系统的成像质量。

为实现上述目的，按照本发明的一个方面，提供了一种超声CT的成像方法，包括以下步骤：

(1)对成像区域进行网格化处理后获得M个网格点；

(2)环形探头以成像区域为中心360°旋转，获得所述M个网格点对应的Q*N*K组数据；其中，表示采集次数，N为环形探头的发射阵元的数量，K为每个发射阵元对应的接收阵元的数量，[.]表示向下取整；

(3)获得所述Q*N*K组数据对应的渡越时间向量T₀以及直线路径向量L₀；

(4)选取初始慢速S₀，以T₀为初始渡越时间、L₀为初始传播路径，迭代求解方程组L_f-1*dS_f＝dT_f，dT_f＝T_f－T_f-1，S_f＝S_f-1+dS_f，直至dT_f≤渡越时间扰动变化阈值ξ，获得实际传播路径L_f；

(5)获得所述M个网格点对应的速度分布V，并将所述速度分布V归一化，完成成像；其中，速度分布V＝1/s，实际慢速s＝T₀/L_f。

优选地，在所述步骤(1)中，所述网格点的边长为ε₁～ε；其中，ε＝max(ε₁,ε₂),ε₁为超声波的半波长，ε₂为成像区域中病变组织尺寸的1/4，以保证成像速度的同时，保持成像区域中的病变组织的清晰度。

优选地，在所述步骤(2)中，所述环形探头旋转的角度为α为与N互质的整数。

优选地，在所述步骤(2)中，K为3～(N/2+1)。

优选地，在所述步骤(4)中，所述初始慢速为水中的声速。

优选地，在所述步骤(4)中，所述渡越时间扰动变化阈值ξ为e^－2～e^－4。

优选地，所述步骤(4)具体包括以下子步骤：

(4-1)令迭代次数f＝1，初始慢速S_f-1＝S₀，初始渡越时间T_f-1＝T₀，初始传播路径L_f-1＝L₀；

(4-2)根据方程组L_f-1*dS_f＝dT_f，S_f＝S_f-1+dS_f，dT_f＝T_f－T_f-1，获得迭代渡越时间，迭代慢速，迭代传播路径L_f以及传播路径渡越时间扰动变化dT_f；

(4-3)判断所述渡越时间扰动变化dT_f是否小于渡越时间扰动变化阈值ξ，是则获得实际传播路径为迭代传播路径L_f，否则f＝f+1，返回步骤(4-2)。

优选地，在所述步骤(5)中归一化的方法为其中，G(i)表示归一化后的速度分布，V(i)表示网格点i的速度分布，i为1～M，min(V(i))为V(i)的最小值，max(V(i))为V(i)的最大值。

优选地，在所述步骤(5)后还包括：将所述速度分布V归一化，并将归一化后的速度分布映射为灰度值。

总体而言，通过本发明所构思的以上技术方案与现有技术相比，具有下列有益效果：

1、由于本发明的成像方法，使得环形探头从多角度采集超声数据，从面能比传统超声获取更加全面的成像数据，弥补了阵元间存在空隙而无法全面扫描到成像区域的问题，有效缓解了声速重建中逆问题的不适定性，从而提高了信噪比和速度重建的精度。

2、由于环形探头多次旋转扫描的次数过多，会造成信息的冗余，重建时求解逆问题的复杂度大大提高，重建效率较低等问题；而扫描次数过少，则会造成成像不清晰；本发明根据求解慢速时的不适定性，获得了最优的采集次数，在提高重建精度的前提下，最大限度地保持了数据的采集效率。

附图说明

图1是实施例1的环形阵列示意图；

图2是实施例1的基于环形阵列多次旋转扫描模式的超声CT声速重建流程图；

图3a是实施例1的环形阵列多次旋转扫描前重建的超声CT声速图像；

图3b是实施例1的环形阵列旋转一次扫描后重建的超声CT声速图像。

具体实施方式

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

本发明提供了一种超声CT的成像方法，包括以下步骤：

(1)将成像区域网格化为M个网格点；网格点的边长越小，成像精度越高，然而会影响其计算效率，因此网格点的边长最小为超声波的半波长，最大不超过超声波的半波长和成像区域中病变组织尺寸的1/4的最大值；

(2)环形探头以成像区域为中心360°旋转，获得所述M个网格点对应的Q*N*K组数据；其中，表示采集次数，N为环形探头的发射阵元的数量，K为3～(N/2+1)，为每个发射阵元对应的接收阵元的数量，通常接收阵元优选为以发射阵元的对面阵元为中心呈π/4～π多个阵元，以在保证成像精度的同时，尽可能的减少计算量，[.]表示向下取整；在实践中，通常令环形探头以固定角度旋转Q－1次，α为与N互质的整数，为操作简易考虑，可令α＝1；

(3)利用加权AIC(Akaike information criterion)渡越时间提取法、采用初至波方法(the first break method)或AIC和归一化相结合法获得所述Q*N*K组数据对应的渡越时间向量T₀以及直线路径向量L₀；

(4-1)令迭代次数f＝1，初始慢速S_f-1＝S₀，初始渡越时间T_f-1＝T₀，初始传播路径L_f-1＝L₀；初始慢速选取任意正数即可，为减少迭代次数，通常以水中的声速为初始慢速；

(4-2)根据方程组L_f-1*dS_f＝dT_f，S_f＝S_f-1+dS_f，dT_f＝T_f－T_f-1，获得迭代渡越时间，迭代慢速，迭代传播路径L_f以及传播路径渡越时间扰动变化dT_f；

(4-3)判断所述渡越时间扰动变化dT_f是否小于渡越时间扰动变化阈值ξ，是则获得实际传播路径为迭代传播路径L_f，否则f＝f+1，返回步骤(4-2)；渡越时间扰动变化阈值ξ为e^－2～e^－4，要求的成像精度越高可将该阈值设置得越小；

(5)获得所述M个网格点对应的速度分布V，完成成像；其中，速度分布V＝1/s，实际慢速s＝T₀/L_f；

(6)将所述速度分布V归一化，并将归一化后的速度分布映射为灰度值；归一化可利用公式其中，G(i)表示归一化后的速度分布，V(i)表示网格点i的速度分布，i为1～M，min(V(i))为V(i)的最小值，max(V(i))为V(i)的最大值；映射的方法通常采用线性映射或动态范围压缩法；以线性映射为例，对于一个灰度级为256的图像，则令max(G(i))对应映射后的灰度的最大值255，min(G(i))对应映射后的灰度的最小值，其它的速度分布则按该比例映射；

当采用线性映射时，也可以直接将归一化前的速度分布映射为灰度值，即灰度值[.]表示向下取整。。

实施例1

本发明提出的探头多次旋转扫描模式的乳腺超声CT的声速成像方法，通过环形探头作一定次数的旋转运动，以此来获得美国Karmanos癌症中心的超声CT成像扫描器CURE无法实现的阵元间间隙所对应的成像物体的信息。同时通过规划最优旋转次数，降低声速重建时逆问题的不适定性，重建出高质量的乳腺超声CT声速图像(本实施例中重建出的成像目标的声速与实际声速的相对误差减小了一倍)。

一种基于环形阵列多次旋转扫描模式的超声CT声速成像方法，其步骤包括最优旋转扫描次数确定、环形阵列多次旋转扫描下的数据采集、子孔径数据选取、渡越时间向量和初始路径矩阵整合、折线路径更新、声速图像评价标准、求逆问题，解速度分布以及图像显示步骤。

(1)子孔径数据的采集

(1-1)规划最优旋转扫描次数

由于多次旋转扫描方法使超声环形探头覆盖到更多的成像区域，因此，这种扫描方法，不仅可以采集到更丰富的信息，而且可以抑制阵元间的稀疏性造成的伪信号，从而提高图像的重建质量。理论上来说，只要每次旋转的角度足够小，旋转扫描次数越多，采集到的信息更全面、更丰富，但是如果扫描次数过多，会造成信息的冗余；采集时间过长，采集效率比较低；重建时求解逆问题的复杂度大大提高，重建效率较低等问题。所以如何选取最优的旋转扫描次数，是本发明优先考虑的一项重要内容。本发明利用声速重建中逆问题的不适定性，基于最大程度地降低逆问题的不适定性来规划最优的旋转扫描次数。

因此为了准确地求解成像物体的声速分布，精确地重建图像，我们尽可能地让路径矩阵的秩等于慢度向量s的维数M，而通过多次旋转环形探头，可以产生更多地从发射阵元到接收阵元的路径，从而路径矩阵的秩可以增大至或接近M。因此，获得最优的旋转扫描次数为：则采集的数据的次数为Q次；

将成像区域网格化，分成M＝(m-1)*(n-1)个网格。其中，M为成像区域剖分的网格数。理论上剖分的网格尺寸可以无限接近0。虽然网格尺寸越小，求解出来的声速分布越精确，但是网格越密，计算量越大，所以最佳的M由成像区域和计算效率共同决定。一般地，病变组织的尺寸较大时，剖分不要求那么密，网格数M可适当降低，只要满足网格尺寸小于病变组织尺寸的1/4即可；当未知成像区域中病变组织尺寸时，剖分的网格大小不超过(λ/2)*(λ/2)，λ为超声波波长。N为环形探头的阵元数量，K为子孔径的阵元数量；本实施例中，N＝72,K＝9,M＝64*64。；因为声速重建是利用透射信号成像，子孔径即为发射阵元对面阵元的邻近阵元，子孔径接收到的信号主要包含透射信息，理论上K最小为3，最大为(N/2+1)。如果K过大，不仅图像的旁瓣水平会很高，重建的计算量也会提高；但K过小，能够利用到的信息较少，降低重建的精度。[.]表示向下取整。

(1-2)数据采集

利用N阵元环形探头，循环采集成像物体与环形探头的阵元一一对应的N组数据；沿逆时针方向(或顺时针方向)旋转角度且a与N互质，再采集成像物体的N组数据；直至沿逆时针方向采集第Q次数据，每个网格共获得Q*N*K组子孔径阵元接收到的信号。

(2)渡越时间向量和初始路径矩阵整合：

将采集到的Q+1组子孔径数据分别利用目前美国Karmanos癌症中心使用的，Cuiping Li等提出的加权AIC(Akaike information criterion)渡越时间提取法(除此之外，也可以采用初至波方法(the first break method)，Xiaolei Qu提出的AIC和归一化互相关结合的方法(AICNCC)等)计算每批数据的渡越时间向量，维数N*K，并按对应的旋转扫描的先后顺序放在总的渡越时间向量矩阵T₀中，维数N*K*(Q+1)，完成渡越时间向量的整合。同时，每组子孔径数据对应的初始路径矩阵计算是借助成像区域网格化完成的。将成像区域网格化，分成m行n列，这样整个成像区域分为M个网格，M＝(m-1)*(n-1)。由于每个网格都对应某一个速度，而速度的倒数是慢度，所以整个成像区域的速度分布可以由一个M*1的慢度向量s来描述。计算环阵上每个阵元在网格中的位置(第几行第几列)；假设超声信号从发射阵元沿直线路径传播至接收阵元，根据路径经过的网格位置及在网格中的路径长度表示每个发射阵元到接收阵元的直线路径，得到每组子孔径数据对应的直线路径矩阵，维数为N*K行M列；最后将这Q+1个路径矩阵也按对应的旋转的先后顺序放在总的直线路径矩阵L₀中，维数为N*K*(Q+1)行，M列，这里的直线路径矩阵是步骤五中折线路径更新的初始路径矩阵。

(3)计算速度分布V

(3-1)折线路径更新：

超声CT的声速重建问题实质为求解如下一类线性方程组：Ls＝T，其中，s为M*1的慢度向量，可反映整个成像目标的速度分布，L是超声波从发射阵元到发射阵元对面的子孔径的接收阵元的折线路径矩阵，T为子孔径中所有从发射阵元到接收阵元超声波的渡越时间向量。

本发明由于规划了最优扫描次数，从而降低了线性方程组的不适定性，从而得到方程组更准确的解。因为声速成像中，我们只利用到了子孔径的信息，所以网格剖分充分密的情况下线性方程组(1)是个欠定问题，即N*K＜M。

而由数值分析可知，设A为p阶方阵，非齐次线性方程组Ax＝b……(2)有唯一解等价于系数矩阵A的秩rank(A)＝p，即系数矩阵的秩等于方程组中未知数的个数。这里A对应重建中的路径矩阵L，x对应重建中的慢度向量s，b对应重建时的实际的渡越时间向量T₀，系数矩阵A的秩对应线性无关的从发射阵元到接受阵元路径的个数。

原理是从发射与接收点之间的假想初始路径开始，根据最小走时准则对路径进行扰动，从而求出接收点处的走时及射线路径。步骤四中计算的总的直线路径矩阵L₀为这里折线路径更新的初始路径矩阵。对由阵元i发射并且被阵元j接收的超声波，它的传播路径L_ij，慢度向量S_M×1和渡越时间T_ij满足：L_ij(1×M)*S_M×1＝T_ij……(4)；这样所有的发射接收阵元组合的路径矩阵L和渡越时间T满足以下等式：L*s＝T……(5)；假设初始的慢度为S₀(任意正数即可，但为了减少迭代次数，加快迭代速度，一般由声波在水中的速度表示)，将初始路径矩阵L₀和S₀向量相乘得到计算的渡越时间T₁，然后将T₁和实际的渡越时间向量作差，得出渡越时间扰动，即dT₁＝T₁-T₀。假设在慢度分布微小扰动时，路径矩阵没有发生变化。解方程组L₀*dS₁＝dT₁，得慢度扰动dS₁，又由S₁＝S₀+dS₁，得更新后的慢度S₁。根据慢度S₁，由最短路径算法重新规划每个发射阵元到接收阵元的最短路径，然后选取每个发射阵元到子孔径中阵元的路径，组成N*K行，M列的路径矩阵，记为L₁。再将L₁和S₁相乘得到计算的渡越时间T₂，计算新的渡越时间扰动，即dT₂＝T₂-T₁，同样利用在慢度分布微小扰动时，路径矩阵没有发生变化的假设，解方程组L₀*dS₂＝dT₂，得慢度扰动dS₂，又由S₂＝S₁+dS₂，得更新后的慢度S₂，重复之前的计算，不断更新折线路径，不断更新慢度分布。

(3-2)：声速图像评价标准：该评价标准指的是步骤五中每次更新折线路径时得到的渡越时间扰动变化dT_f的范数上限ξ，计算的渡越时间和实际的渡越时间之差称为渡越时间的扰动。渡越时间扰动变化的上限是一个很小的数，ξ一般小于e^-2(在本实施例中ξ＝e^-2)，当渡越时间扰动变化不超过上限时，这时就认为达到我们所要求的精度，就可以进行下一步图像显示；如果渡越时间扰动变化超过规定的上限，则继续更新折线路径，更新f-1次，直到扰动变化dT_f不超过上限，这时的折线路径矩阵记为L_f，对应的慢度为S_f。

(3-3)求逆问题，解速度分布：通过步骤五和步骤六中得到的最终的折线路径矩阵L_f和步骤四中得到的渡越时间向量T₀，通过一种成熟的方法-TVAL3方法(其他还有非线性共轭梯度方法，联合代数重建算法等)求解逆问题：L_f*s＝T₀……(6)得出慢度分布s，利用V＝1/s，将s转化成速度向量V，维数为(m-1)*(n-1)。

(4)图像显示：该步骤包括将步骤七中得到的速度向量归一化及灰度映射，最终得到超声CT图像。为了提高图像的对比度，我们采用将速度向量V归一化。将V的最小值min(V)，最大值max(V)取出来，然后按照(7)中函数将其映射到[0,1]，并最终生成(m-1)*(n-1)的超声CT图像。灰度映射采用简单的线性映射(其他还有动态范围压缩法等)，即成比例地将最小的速度映射到0，将最大的速度映射到255或511(分别对应图像显示中常用的256灰度级或512灰度级的图像)。

$f\left(V\right)=\frac{V-min\left(V\right)}{max\left(V\right)-\min \left(V\right)}---\left(7\right)$

从图3中我们可以看出旋转扫描一次对声速重建图像的影响，可以看出，经一次扫描，声速重建后的图像质量即获得了极大的提升。

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

Claims (9)

1.一种超声CT的成像方法，其特征在于，包括以下步骤：

(1)对成像区域进行网格化处理后获得M个网格点；

(2)环形探头以成像区域为中心360°旋转，获得所述M个网格点对应的Q*N*K组数据；其中，表示采集次数，N为环形探头的发射阵元的数量，K为每个发射阵元对应的接收阵元的数量，[.]表示向下取整；

(3)获得所述Q*N*K组数据对应的渡越时间向量T₀以及直线路径向量L₀；

(4)选取初始慢速S₀，以T₀为初始渡越时间、L₀为初始传播路径，迭代求解方程组L_f-1*dS_f＝dT_f，dT_f＝T_f－T_f-1，S_f＝S_f-1+dS_f，直至dT_f≤渡越时间扰动变化阈值ξ，获得实际传播路径L_f；

(5)获得所述M个网格点对应的速度分布V，完成成像；其中，速度分布V＝1/s，实际慢速s＝T₀/L_f。

2.如权利要求1所述的成像方法，其特征在于，在所述步骤(1)中，所述网格点的边长为ε₁～ε；其中，ε＝max(ε₁,ε₂),ε₁为超声波的半波长，ε₂为成像区域中病变组织尺寸的1/4。

3.如权利要求1所述的成像方法，其特征在于，在所述步骤(2)中，所述环形探头旋转的角度为α为与N互质的整数。

4.如权利要求1所述的成像方法，其特征在于，在所述步骤(2)中，K为3～(N/2+1)。

5.如权利要求1所述的成像方法，其特征在于，在所述步骤(4)中，所述初始慢速为水中的声速。

6.如权利要求1所述的成像方法，其特征在于，在所述步骤(4)中，所述渡越时间扰动变化阈值ξ为e^－2～e^－4。

7.如权利要求1所述的成像方法，其特征在于，所述步骤(4)具体包括以下子步骤：

(4-1)令迭代次数f＝1，初始慢速S_f-1＝S₀，初始渡越时间T_f-1＝T₀，初始传播路径L_f-1＝L₀；

(4-2)根据方程组L_f-1*dS_f＝dT_f，S_f＝S_f-1+dS_f，dT_f＝T_f－T_f-1，获得迭代渡越时间，迭代慢速，迭代传播路径L_f以及传播路径渡越时间扰动变化dT_f；

(4-3)判断所述渡越时间扰动变化dT_f是否小于渡越时间扰动变化阈值ξ，是则获得实际传播路径为迭代传播路径L_f，否则f＝f+1，返回步骤(4-2)。

8.如权利要求1所述的成像方法，其特征在于，在所述步骤(5)中归一化的方法为其中，G(i)表示归一化后的速度分布，V(i)表示网格点i的速度分布，i为1～M，min(V(i))为V(i)的最小值，max(V(i))为V(i)的最大值。

9.如权利要求1-8中任意一项所述的成像方法，其特征在于，在所述步骤(5)后还包括：将所述速度分布V归一化，并将归一化后的速度分布映射为灰度值。