一种超声ct的成像方法
【专利摘要】本发明公开了一种超声CT的成像方法。包括将成像区域网格化为M个网格点；环形探头以成像区域为中心360°旋转，获得Q*N*K组数据；然后获得该数据对应的渡越时间向量T0以及直线路径向量L0；最后迭代计算，获得所述M个网格点对应的速度分布V，并将所述速度分布V归一化，完成成像。本发明通过优化环形探头中发射阵元以及接收阵元的采集次数和采集角度，并对获得的数据进行处理与图像重建，从而使重建出来的速度分布更精确，以提高成像的信噪比，提高系统的成像质量。
【专利说明】
一种超声CT的成像方法
技术领域
[0001 ]本发明属于超声断层成像领域，更具体地，涉及一种超声CT的成像方法。
【背景技术】
[0002] 超声CT成像，即超声断层成像，具有超声非侵入、无辐射和成本低的优势，以断层 成像的方式进行扫描，在乳腺癌等疾病早期诊断方面具有广阔的应用前景。
[0003] 超声CT成像分为反射式成像和透射式成像，其中，采用声速的透射式成像是通过 重建乳腺内组织的声速分布和衰减系数进行成像。超声CT声速成像利用首达波到达时间， 重建乳腺内组织的声速分布，从而正确地区分良性肿瘤与癌。因为和正常组织与良性肿瘤 相比，癌具有较高的声速。
[0004] 目前由在乳腺癌早期检测领域走在前列的美国Karmano s癌症中心研制的环形阵 列的超声CT成像扫描器CURE可以得到亚毫米级高分辨率的图像。这种扫描器每次使用单阵 元向环形阵列的中心发射超声波，所有阵元接收，直到每个阵元循环发射一次，所以可以接 收到多个角度的超声信号，得到斑点噪声少的高分辨率图像。虽然这种扫描器实现了对成 像物体的360度扫描，但由于环形探头之间存在间隙，从而无法获取到更多的成像目标的信 息，进而声速重建时由于方程组的解并非唯一，所以求解出来的声速分布跟实际的声速分 布差别很大，重建图像的质量不高。
【发明内容】
[0005] 针对现有技术的以上缺陷或改进需求，本发明提供了一种超声CT的成像方法，其 目的在于优化环形探头中发射阵元以及接收阵元的采集次数和采集角度，并对获得的数据 进行处理与图像重建，从而使重建出来的速度分布更精确，以提高成像的信噪比，提高了系 统的成像质量。
[0006] 为实现上述目的，按照本发明的一个方面，提供了一种超声CT的成像方法，包括以 下步骤：
[0007] (1)对成像区域进行网格化处理后获得Μ个网格点；
[0008] (2)环形探头以成像区域为中心360°旋转，获得所述Μ个网格点对应的Q*N*K组数 据;其中，0 ，表示采集次数，N为环形探头的发射阵元的数量，每个发射阵元 对应的接收阵元的数量，[.]表示向下取整；
[0009] (3)获得所述Q*N*K组数据对应的渡越时间向量To以及直线路径向量L〇;
[0010] (4)选取初始慢速SQ，以To为初始渡越时间、L〇为初始传播路径，迭代求解方程组 Lf-wdSf = dTf，dTf = Tf - Tf-i，Sf = Sf-ddSf，直至dTf <渡越时间扰动变化阈值ξ，获得实际传 播路径Lf;
[0011] (5)获得所述Μ个网格点对应的速度分布V，并将所述速度分布V归一化，完成成像； 其中，速度分布V = 1 /s，实际慢速s = TQ/Lf。
[0012]优选地，在所述步骤⑴中，所述网格点的边长为εχ~ε ;其中，ε ε2)，￡丄为 超声波的半波长，ε2为成像区域中病变组织尺寸的1/4,以保证成像速度的同时，保持成像 区域中的病变组织的清晰度。
[0013] 优选地，在所述步骤（2)中，所述环形探头旋转的角度为#?，α为与Ν互质的整 数。
[0014] 优选地，在所述步骤(2)中，Κ为3~(Ν/2+1)。
[0015] 优选地，在所述步骤(4)中，所述初始慢速为水中的声速。
[0016] 优选地，在所述步骤(4)中，所述渡越时间扰动变化阈值ξ为eT2~eT4。
[0017] 优选地，所述步骤(4)具体包括以下子步骤：
[0018] (4-1)令迭代次数f = 1，初始慢速Sf-i = So，初始渡越时间Tf-i = To，初始传播路径 Lf-i - Lo ；
[0019] (4-2)根据方程组 Lf-AdSFdThSFSf-i+dSndTFTf -Tf-!，获得迭代渡越时间， 迭代慢速，迭代传播路径Lf以及传播路径渡越时间扰动变化dTf;
[0020] (4-3)判断所述渡越时间扰动变化dTf是否小于渡越时间扰动变化阈值ξ，是则获 得实际传播路径为迭代传播路径Lf，否则f = f+l，返回步骤(4-2)。
[0021] 优选地，在所述步骤(5)中归一化的方法戈
，其中，G(i) 表示归一化后的速度分布，V(i)表示网格点i的速度分布，i为1~M，min(V(i))为V(i)的最 小值，max(V(i))为V(i)的最大值。
[0022] 优选地，在所述步骤(5)后还包括:将所述速度分布V归一化，并将归一化后的速度 分布映射为灰度值。
[0023] 总体而言，通过本发明所构思的以上技术方案与现有技术相比，具有下列有益效 果：
[0024] 1、由于本发明的成像方法，使得环形探头从多角度采集超声数据，从面能比传统 超声获取更加全面的成像数据，弥补了阵元间存在空隙而无法全面扫描到成像区域的问 题，有效缓解了声速重建中逆问题的不适定性，从而提高了信噪比和速度重建的精度。
[0025] 2、由于环形探头多次旋转扫描的次数过多，会造成信息的冗余，重建时求解逆问 题的复杂度大大提高，重建效率较低等问题;而扫描次数过少，则会造成成像不清晰;本发 明根据求解慢速时的不适定性，获得了最优的采集次数，在提高重建精度的前提下，最大限 度地保持了数据的采集效率。
【附图说明】
[0026] 图1是实施例1的环形阵列示意图；
[0027] 图2是实施例1的基于环形阵列多次旋转扫描模式的超声CT声速重建流程图；
[0028] 图3a是实施例1的环形阵列多次旋转扫描前重建的超声CT声速图像；
[0029] 图3b是实施例1的环形阵列旋转一次扫描后重建的超声CT声速图像。
【具体实施方式】
[0030] 为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对 本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并 不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要 彼此之间未构成冲突就可以相互组合。
[0031] 本发明提供了一种超声CT的成像方法，包括以下步骤：
[0032] (1)将成像区域网格化为Μ个网格点；网格点的边长越小，成像精度越高，然而会影 响其计算效率，因此网格点的边长最小为超声波的半波长，最大不超过超声波的半波长和 成像区域中病变组织尺寸的1 /4的最大值；
[0033] (2)环形探头以成像区域为中心360°旋转，获得所述Μ个网格点对应的Q*N*K组数 据;其中
，表示采集次数，N为环形探头的发射阵元的数量，K为3~(N/2+1)，为 每个发射阵元对应的接收阵元的数量，通常接收阵元优选为以发射阵元的对面阵元为中心 呈V4~31多个阵元，以在保证成像精度的同时，尽可能的减少计算量，[.]表示向下取整;在 实践中，通常令环形探头以固定角度#?旋转Q - 1次，α为与N互质的整数，为操作简易考 虑，可令α = 1;
[0034] (3)利用加权AIC(Akaike information criterion)渡越时间提取法、采用初至波 方法(the first break method)或AIC和归一化相结合法获得所述Q*N*K组数据对应的渡 越时间向量Το以及直线路径向量Lo;
[0035] (4-1)令迭代次数f = 1，初始慢速Sf-1 = So，初始渡越时间Tf-ι = Το，初始传播路径 LhiLo;初始慢速选取任意正数即可，为减少迭代次数，通常以水中的声速为初始慢速；
[0036] (4-2)根据方程组 Lf-AdSf = dTf，Sf = Sf-i+dSf，dTf = Tf - Tf-!，获得迭代渡越时间， 迭代慢速，迭代传播路径Lf以及传播路径渡越时间扰动变化dTf;
[0037] (4-3)判断所述渡越时间扰动变化dTf是否小于渡越时间扰动变化阈值ξ，是则获 得实际传播路径为迭代传播路径L f，否则f = f+l，返回步骤(4-2);渡越时间扰动变化阈值ξ 为eT2~eT4,要求的成像精度越高可将该阈值设置得越小；
[0038] (5)获得所述Μ个网格点对应的速度分布V，完成成像;其中，速度分布V=l/s，实际 慢速 s = To/Lf;
[0039] (6)将所述速度分布V归一化，并将归一化后的速度分布映射为灰度值；归一化可 利用公￥
其中，G(i)表不归一化后的速度分布，V(i)表不网格点 i的速度分布，i为1~111^11(￥(1))为￥(1)的最小值，111&以￥(1))为￥(1)的最大值;映射的方 法通常采用线性映射或动态范围压缩法；以线性映射为例，对于一个灰度级为256的图像， 则令max(G(i))对应映射后的灰度的最大值255，11^11(6(1))对应映射后的灰度的最小值，其 它的速度分布则按该比例映射；
[0040] 当采用线性映射时，也可以直接将归一化前的速度分布映射为灰度值，即灰度值
[·]表示向下取整。。
[0041 ] 实施例1
[0042] 本发明提出的探头多次旋转扫描模式的乳腺超声CT的声速成像方法，通过环形探 头作一定次数的旋转运动，以此来获得美国Karmanos癌症中心的超声CT成像扫描器CURE无 法实现的阵元间间隙所对应的成像物体的信息。同时通过规划最优旋转次数，降低声速重 建时逆问题的不适定性，重建出高质量的乳腺超声CT声速图像(本实施例中重建出的成像 目标的声速与实际声速的相对误差减小了一倍）。
[0043] -种基于环形阵列多次旋转扫描模式的超声CT声速成像方法，其步骤包括最优旋 转扫描次数确定、环形阵列多次旋转扫描下的数据采集、子孔径数据选取、渡越时间向量和 初始路径矩阵整合、折线路径更新、声速图像评价标准、求逆问题，解速度分布以及图像显 示步骤。
[0044] (1)子孔径数据的采集
[0045] (1-1)规划最优旋转扫描次数
[0046] 由于多次旋转扫描方法使超声环形探头覆盖到更多的成像区域，因此，这种扫描 方法，不仅可以采集到更丰富的信息，而且可以抑制阵元间的稀疏性造成的伪信号，从而提 高图像的重建质量。理论上来说，只要每次旋转的角度足够小，旋转扫描次数越多，采集到 的信息更全面、更丰富，但是如果扫描次数过多，会造成信息的冗余;采集时间过长，采集效 率比较低;重建时求解逆问题的复杂度大大提高，重建效率较低等问题。所以如何选取最优 的旋转扫描次数，是本发明优先考虑的一项重要内容。本发明利用声速重建中逆问题的不 适定性，基于最大程度地降低逆问题的不适定性来规划最优的旋转扫描次数。
[0047] 因此为了准确地求解成像物体的声速分布，精确地重建图像，我们尽可能地让路 径矩阵的秩等于慢度向量s的维数M，而通过多次旋转环形探头，可以产生更多地从发射阵 元到接收阵元的路径，从而路径矩阵的秩可以增大至或接近M。因此，获得最优的旋转扫描 次数为:……，则采集的数据的次数为Q次；
[0048] 将成像区域网格化，分成M = (m-1) * (η-1)个网格。其中，Μ为成像区域剖分的网格 数。理论上剖分的网格尺寸可以无限接近0。虽然网格尺寸越小，求解出来的声速分布越精 确，但是网格越密，计算量越大，所以最佳的Μ由成像区域和计算效率共同决定。一般地，病 变组织的尺寸较大时，剖分不要求那么密，网格数Μ可适当降低，只要满足网格尺寸小于病 变组织尺寸的1/4即可；当未知成像区域中病变组织尺寸时，剖分的网格大小不超过(λ/2)* (λ/2)，λ为超声波波长。Ν为环形探头的阵元数量，Κ为子孔径的阵元数量;本实施例中，Ν = 72，Κ = 9，Μ = 64*64。；因为声速重建是利用透射信号成像，子孔径即为发射阵元对面阵元的 邻近阵元，子孔径接收到的信号主要包含透射信息，理论上Κ最小为3,最大为(Ν/2+1)。如果 Κ过大，不仅图像的旁瓣水平会很高，重建的计算量也会提高;但Κ过小，能够利用到的信息 较少，降低重建的精度。[.]表示向下取整。
[0049] (1-2)数据采集
[0050]利用Ν阵元环形探头，循环采集成像物体与环形探头的阵元一一对应的Ν组数据； 沿逆时针方向（或顺时针方向）旋转角度·^α，且a与Ν互质，再采集成像物体的Ν组数据;直 至沿逆时针方向采集第Q次数据，每个网格共获得Q*N*K组子孔径阵元接收到的信号。
[0051 ] (2)渡越时间向量和初始路径矩阵整合：
[0052]将采集到的Q+1组子孔径数据分别利用目前美国Karmanos癌症中心使用的， Cuiping Li等提出的加权AIC(Akaike information criterion)渡越时间提取法(除此之 外，也可以采用初至波方法(the first break method)，Xiaolei Qu提出的AIC和归一化互 相关结合的方法(AICNCC)等)计算每批数据的渡越时间向量，维数N*K，并按对应的旋转扫 描的先后顺序放在总的渡越时间向量矩阵To中，维数N*K*(Q+1 )，完成渡越时间向量的整 合。同时，每组子孔径数据对应的初始路径矩阵计算是借助成像区域网格化完成的。将成像 区域网格化，分成m行η列，这样整个成像区域分为Μ个网格，M = (m-1) * (η-1)。由于每个网格 都对应某一个速度，而速度的倒数是慢度，所以整个成像区域的速度分布可以由一个M*1的 慢度向量s来描述。计算环阵上每个阵元在网格中的位置(第几行第几列）；假设超声信号从 发射阵元沿直线路径传播至接收阵元，根据路径经过的网格位置及在网格中的路径长度表 示每个发射阵元到接收阵元的直线路径，得到每组子孔径数据对应的直线路径矩阵，维数 为N*K行Μ列；最后将这Q+1个路径矩阵也按对应的旋转的先后顺序放在总的直线路径矩阵 Lo中，维数为N*K*(Q+1)行，Μ列，这里的直线路径矩阵是步骤五中折线路径更新的初始路径 矩阵。
[0053] (3)计算速度分布V
[0054] (3-1)折线路径更新：
[0055]超声CT的声速重建问题实质为求解如下一类线性方程组:Ls = T，其中，s为M*1的 慢度向量，可反映整个成像目标的速度分布，L是超声波从发射阵元到发射阵元对面的子孔 径的接收阵元的折线路径矩阵，T为子孔径中所有从发射阵元到接收阵元超声波的渡越时 间向量。
[0056] 本发明由于规划了最优扫描次数，从而降低了线性方程组的不适定性，从而得到 方程组更准确的解。因为声速成像中，我们只利用到了子孔径的信息，所以网格剖分充分密 的情况下线性方程组(1)是个欠定问题，即N*K<M。
[0057] 而由数值分析可知，设A为p阶方阵，非齐次线性方程组Ax = b……（2)有唯一解等 价于系数矩阵A的秩rank(A)=p，即系数矩阵的秩等于方程组中未知数的个数。这里A对应 重建中的路径矩阵L，X对应重建中的慢度向量s，b对应重建时的实际的渡越时间向量To，系 数矩阵A的秩对应线性无关的从发射阵元到接受阵元路径的个数。
[0058]原理是从发射与接收点之间的假想初始路径开始，根据最小走时准则对路径进行 扰动，从而求出接收点处的走时及射线路径。步骤四中计算的总的直线路径矩阵Lo为这里 折线路径更新的初始路径矩阵。对由阵元i发射并且被阵元j接收的超声波，它的传播路径 Lij，慢度向量Smxi和渡越时间Tij满足:Lij(ixM)*SMxi = Tij......⑷；这样所有的发射接收阵兀 组合的路径矩阵L和渡越时间T满足以下等式:L*s = T……（5);假设初始的慢度为So(任意 正数即可，但为了减少迭代次数，加快迭代速度，一般由声波在水中的速度表示），将初始路 径矩阵Lo和So向量相乘得到计算的渡越时间!^，然后将1^和实际的渡越时间向量作差，得出 渡越时间扰动，即cmiMo。假设在慢度分布微小扰动时，路径矩阵没有发生变化。解方程 组=cm，得慢度扰动dSi，又由Si=So+dSi，得更新后的慢度Si。根据慢度Si，由最短路 径算法重新规划每个发射阵元到接收阵元的最短路径，然后选取每个发射阵元到子孔径中 阵元的路径，组成N*K行，Μ列的路径矩阵，记为U。再将LdPSi相乘得到计算的渡越时间T 2， 计算新的渡越时间扰动，即dTsiTs-Ti，同样利用在慢度分布微小扰动时，路径矩阵没有发 生变化的假设，解方程组Lo*dS 2 = dT2，得慢度扰动dS2，又由S2 = Si+dS〗，得更新后的慢度S2， 重复之前的计算，不断更新折线路径，不断更新慢度分布。
[0059] (3-2):声速图像评价标准:该评价标准指的是步骤五中每次更新折线路径时得到 的渡越时间扰动变化dT f的范数上限ξ，计算的渡越时间和实际的渡越时间之差称为渡越时 间的扰动。渡越时间扰动变化的上限是一个很小的数，ξ-般小于eT 2(在本实施例中| = e _2)，当渡越时间扰动变化不超过上限时，这时就认为达到我们所要求的精度，就可以进行下 一步图像显示;如果渡越时间扰动变化超过规定的上限，则继续更新折线路径，更新f-Ι次， 直到扰动变化dT f不超过上限，这时的折线路径矩阵记为Lf，对应的慢度为&。
[0060] (3-3)求逆问题，解速度分布:通过步骤五和步骤六中得到的最终的折线路径矩阵 Lf和步骤四中得到的渡越时间向量Το，通过一种成熟的方法-TVAL3方法(其他还有非线性共 辄梯度方法，联合代数重建算法等)求解逆问题:L f*S = T〇……（6)得出慢度分布s，利用V = 1/8，将8转化成速度向量￥，维数为(111-1)*(11-1)。
[0061] (4)图像显示:该步骤包括将步骤七中得到的速度向量归一化及灰度映射，最终得 到超声CT图像。为了提高图像的对比度，我们采用将速度向量V归一化。将V的最小值min (V)，最大值max(V)取出来，然后按照(7)中函数将其映射到[0，1]，并最终生成(m-l)*(n-l) 的超声CT图像。灰度映射采用简单的线性映射(其他还有动态范围压缩法等），即成比例地 将最小的速度映射到〇,将最大的速度映射到255或511 (分别对应图像显示中常用的256灰 度级或512灰度级的图像）。
[0062]
[0063]从图3中我们可以看出旋转扫描一次对声速重建图像的影响，可以看出，经一次扫 描，声速重建后的图像质量即获得了极大的提升。
[0064]本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以 限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含 在本发明的保护范围之内。
【主权项】
1. 一种超声CT的成像方法，其特征在于，包括以下步骤： (1) 对成像区域进行网格化处理后获得M个网格点； (2) 环形探头以成像区域为中心360°旋转，获得所述M个网格点对应的Q*N*K组数据;其 中，表示采集次数，N为环形探头的发射阵元的数量，K为每个发射阵元对应的 接收阵元的数量，[.]表示向下取整； (3) 获得所述Q*N*K组数据对应的渡越时间向量To以及直线路径向量L0; (4) 选取初始慢速So，以To为初始渡越时间、Lo为初始传播路径，迭代求解方程组Lf-AdSf = dTf，dTf = Tf - Tf-! ,Sf = Sf-WdSf，直至dTf<渡越时间扰动变化阈值ξ，获得实际传播路径 Lf; (5) 获得所述M个网格点对应的速度分布V，完成成像;其中，速度分布V= Ι/s，实际慢速 s = To/Lf 〇2. 如权利要求1所述的成像方法，其特征在于，在所述步骤（1)中，所述网格点的边长为 ει~ε ;其中，E=Iiiax(EU2) J1为超声波的半波长，ε2为成像区域中病变组织尺寸的1/4。3. 如权利要求1所述的成像方法，其特征在于，在所述步骤(2)中，所述环形探头旋转的 角度戈，α为与N互质的整数。4. 如权利要求1所述的成像方法，其特征在于，在所述步骤(2)中，K为3~(Ν/2+1)。5. 如权利要求1所述的成像方法，其特征在于，在所述步骤(4)中，所述初始慢速为水中 的声速。6. 如权利要求1所述的成像方法，其特征在于，在所述步骤(4)中，所述渡越时间扰动变 化阈值ξ为e_2~e_4。7. 如权利要求1所述的成像方法，其特征在于，所述步骤(4)具体包括以下子步骤： (4-1)令迭代次数f = 1，初始慢速Sf-i = So，初始渡越时间Tf-i = To，初始传播路径Lf-i = Lo； (4-2)根据方程组Lf-AiSf = dTf ,Sf = Sf-WdSf，dTf = Tf- Τη，获得迭代渡越时间，迭代 慢速，迭代传播路径Lf以及传播路径渡越时间扰动变化dTf; (4-3)判断所述渡越时间扰动变化dTf是否小于渡越时间扰动变化阈值ξ，是则获得实际 传播路径为迭代传播路径Lf，否则f = f+l，返回步骤(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归一化，并将归一化后的速度分布映射为灰度值。
【文档编号】A61B8/15GK105997153SQ201610559079
【公开日】2016年10月12日
【申请日】2016年7月15日
【发明人】尉迟明, 丁明跃, 娄翠娟, 王珊珊, 宋俊杰, 李春雨, 方小悦, 彭杨
【申请人】华中科技大学