时间反转超声成像检测乳腺微钙化的仿真

吴水才，欧阳亚丽，吴薇薇，周著黄

(1.北京工业大学生命科学与生物工程学院，智能化生理测量与临床转化北京市国际科研合作基地，北京 100124；2.首都医科大学生物医学工程学院，北京 100069)

乳腺癌是影响中国女性身体健康的主要恶性肿瘤之一，在中国女性肿瘤患者中的发病率居首位，并呈现出逐年上升的趋势[1].乳腺微钙化是早期乳腺癌的主要影像学特征之一，有研究发现超过40%的乳腺癌病例以微钙化为第一特征[2]，并且95%的乳腺导管内癌是通过分析乳腺X线摄影所检测到的微钙化而诊断出[3].微钙化的检出对于早期诊断乳腺癌具有非常重要的意义.

乳腺微钙化通常指超声或乳腺X线摄影图像上直径为0.1～1.0 mm的微小钙化灶[4-5].根据化学成分的不同乳腺微钙化主要分为2种：一种以草酸钙为主，主要常见于乳腺良性病变；另一种以羟基磷灰石为主，在良性和恶性肿瘤中均可见[6].乳腺微钙化可能存在于乳腺肿瘤中，也可能孤立存在，必须根据其大小、数量、分布以及形态来确定乳腺微钙化是良性还是恶性[7].目前，乳腺X线摄影被认为是检测乳腺微钙化的金标准[8].然而，乳腺X线摄影使病人暴露于电离辐射中，对致密型的女性乳腺，尤其是年轻女性，诊断效果相对较差[9]，误诊率为15%～25%，造成不必要的活检[10-11].超声成像是一种非电离技术，是检测乳腺微钙化较为安全的方法，并能对致密乳腺组织进行有效成像.然而，传统超声成像(B超)受空间分辨率和散斑噪声等因素影响，很难检出乳腺微钙化[3].

时间反转(time reversal, TR)方法基于波动方程时反不变性原理，在对散射子位置先验知识未知的情况下，对散射子实现时间和空间上的能量聚焦，从而有效地抑制多径的干扰以及信号散射[12].在医学超声成像[13-14]、碎石术[15]、地下探测[16]及工业检测[17]等领域引起了广泛关注.Robert等[13]提出一种改进的时间反转算子分解(decomposition of the time reversal operator, DORT)算法，用于检测乳腺微钙化，提高了传统DORT的抗噪能力，但该方法只适用于聚焦发射，并且只有在介质中的散斑噪声很强时才有效.Labyed等[14]通过数值模拟研究发现，基于时间反转多信号分类(time reversal with multiple signal classification, TR-MUSIC)算法能有效估算散射子的密度和压缩比，可用于区分乳腺钙化与其他组织散射子.TR-MUSIC算法最早由Devaney[18]提出，通过将时间反转聚焦[19]与多信号分类法[20]相结合，实现超分辨率成像.

TR-MUSIC与DORT算法均基于时间反转技术，通过获取时间反转算子进行奇异值分解，划分信号子空间与噪声子空间，并且信噪子空间之间相互正交.DORT方法利用信号子空间产生后向传播的回传信号，实现选择性聚焦，但当多个散射子之间的距离较小时，很难进行成像.而TR-MUSIC方法利用噪声子空间与目标背景格林函数向量空间的正交特性，来定义成像伪谱，成像质量得到了显著提高.

然而，TR-MUSIC仅适用于目标数量少于换能器阵元数量，因此Labyed等[21]提出窗口化TR-MUSIC成像对其进行改进，通过将成像平面分割成多个子区域，并将TR-MUSIC算法应用于每个子区域，将所有子区域的图像组合在一起形成整个图像，使得目标数量大于换能器阵元数量时，图像质量也可以得到显著提高.Asgedom等[22]提出利用相位相干时间反转多信号分类(phase-coherent with multiple signal classification, PC-MUSIC)法对散射子进行检测与定位，提高了传统TR-MUSIC的纵向分辨率.Foroozan等[23]提出频率聚焦的TR-MUSIC成像算法，降低了TR-MUSIC的计算复杂度.

本文基于TR超声成像算法通过自适应聚焦实现超分辨率成像的特性，通过Field Ⅱ超声仿真软件[24]，设置成像区域，采用全矩阵采集回波信号的方式进行成像仿真，定量分析不同阵元数量、阵元间距以及阵元输出脉冲信号周期个数对成像分辨率的影响.设置点散射子模拟乳腺微钙化，分析不同阈值、不同频带宽度对TR超声成像的影响，从而选择出TR超声成像的最佳成像参数，提高成像质量.并通过对比传统超声成像、TR-MUSIC成像以及PC-MUSIC成像，验证TR超声成像对乳腺微钙化检测的有效性.

1 TR超声成像算法与相关原理

1.1 全矩阵数据采集理论

全矩阵采集(full matrix capture, FMC)技术是一种新的超声数据采集方法，源于合成孔径技术，采用单阵元发射、全部阵元接收[25].以具有N个阵元的超声探头为例，首先从第一个阵元开始激发，所有阵元接收回波信号(A-scan信号)，并进行存储，获得N列回波数据S1j，其中，j=1,2,…,N.然后激发第2个阵元，所有阵元再依次接收并保存接收数据，获得N列回波数据S2j，重复该发射接收过程，直到最后一个阵元激发并被所有阵元接收完成后结束，因而采集的超声阵列数据可包含N×N个时域信号，称为全矩阵数据[26]，如图1所示.

采用全矩阵数据采集方式获得的全矩阵数据，可保存为一个三维阵列

(1)

式中：i,j=1, 2, …,N；S(t,xi,xj)为第i阵元激发、第j阵元接收到的回波数据；(xi, 0)为发射阵元坐标；(xj, 0)为接收阵元坐标；t为时间.

相比于传统B超扫查方式获得的数据，FMC能够最大程度上获得探头中所有发射-接收对的A-scan数据结合，所携带的与点散射子相关的信息量最大，从而为后处理提供了更加完整的数据信息.全矩阵采集是进行数据后处理成像算法研究的第1步，也是超声阵列后处理检测方法的基础.

1.2 TR原理

TR法是通过激发阵列换能器某一阵元，发出声波信号，沿不同路径传播，当遇到点散射子时发生散射，散射信号再由所有阵元接收，将接收信号进行时间上的逆序操作，按照先到后发、后到先发的方式，分别由相应的阵元重新发射出去.根据波动方程的时反不变性原则，由不同阵元发出、沿不同路径传播的声波信号将同时到达声源位置，即在传播介质和阵列换能器先验知识未知的情况下，实现声束的自聚焦[19].

1.3 阵列响应矩阵与TR矩阵

假设使用N个阵元组成超声阵列换能器，对含有M个点散射子的非均匀介质进行检测，采用全矩阵数据采集方式依次激发每一个阵元，由其引起的背散射信号被阵列的N个阵元同时接收，并对接收信号进行一维傅里叶变换，即可获得角频率ω时的阵列响应矩阵Ki,j(ω)，其中，i,j=1, 2, …,N.玻恩近似下，阵列响应矩阵K[27]可表示为

(2)

式中：F(ω)为发射脉冲和阵列响应之间的传递函数；Φ(ω)为F(ω)的相位响应；γk(rm)为波动函数，是对散射子反射率的一种度量；i为虚数单位.

假设密度波动可以忽略不计，对于均匀激励的阵列换能器，方向向量GT为格林函数在换能器表面上的积分[28]

(3)

由于时域信号的TR对应频域信号的相位共轭，TR矩阵T可由阵列响应矩阵[29]

T=KHK=K*K

(4)

推得.式中H与“*”分别为矩阵的伴随算子与矩阵的共轭算子.

根据TR矩阵与阵列响应矩阵的相关性，多重信号分类算法可通过对阵列响应矩阵K进行奇异值分解进而成像.奇异分解表达式[30]表示为

K=U∑VH

(5)

U和V分别由左奇异向量uj和右奇异向量vj组成，根据奇异值大小，可将uj与vj划分为两部分：

U=[US|UN]=[u1,u2,…,uM|uM+1,…,uN]

(6a)

V=[VS|VN]=[v1,v2,…,vM|vM+1,…,vN]

(6b)

式中：US和VS为信号子空间(奇异值非零)；UN和VN为噪声子空间(奇异值为零).对于均匀介质，检测到的目标个数应等于非零奇异值个数[31].但在实际应用中，受背景噪声的影响，产生附加奇异值，使得信号子空间与噪声子空间难以区分，在这种情况下，可设置相应阈值进行划分，并将小于阈值的奇异值处理为零.

1.4 TR超声成像算法

1.4.1 TR-MUSIC成像

TR-MUSIC成像算法可分为2种，单频TR多信号分类(CF-TR-MUSIC)算法和多频TR多信号分类超声(MF-TR-MUSIC)成像算法.

CF-TR-MUSIC成像工作频率通常设定为换能器的中心频率ωc[32].根据信号子空间与噪声子空间的正交性，其在成像区域内任意一点的强度可表示为[33]

(7)

式中：PS(r,ωc)为成像区域任意一点r处的幅值；操作算子A(r,ωc)可由信号子空间US或VS计算得到

(8a)

(8b)

式中“‖‖”表示欧式范数.

MF-TR-MUSIC成像将给定的频率带宽Δω进行组分，选择其中的单点频率ω0为工作频率(ω0∈Δω)，构建阵列响应矩阵K，对K进行奇异值分解划分信噪子空间，根据式(7)可实现成像.将不同工作频率下得到的超声图像组合，即可得到MF-TR-MUSIC成像方法在成像区域内任意点r处的幅值.成像函数[23]可表示为

(9)

式中：Δω为频率带宽；Nω为单点频率个数.

1.4.2 PC-MUSIC成像

PC-MUSIC是对MF-TR-MUSIC的进一步修改，通过引入归一化混合阵列TR操作算子A，保留相位信息，操作算子A[23]表示为

(10)

2 仿真成像

2.1 基本仿真参数分析

参数的选择影响阵列的声学特性，若参数选择不当，阵列声场将会产生伪像，影响成像分辨率，因此，从声学角度讲，为了得到高分辨率的超声图像，首先需要选择合适的阵列参数.

本文通过Field Ⅱ超声仿真软件，以全矩阵采集方式进行超声仿真成像，分别定量分析阵元数量、阵元间距以及阵元输出超声脉冲周期个数对超声成像结果的影响，提高成像质量.

2.1.1 阵元数量分析

阵列包含阵元数量的多少决定了阵列的孔径宽度，阵元数目越多，即阵列孔径越宽.分别使用ele=16、32、64、128阵元的超声阵列，对一个理想点散射子(x=0,z=22)进行全矩阵采集成像，成像结果如图2所示.其中，图2(a)与(b)给出了阵元分别为16和128时的超声成像，以40 dB动态范围显示.图2(c)给出了z=22 mm处不同阵元数量下的横向强度曲线，图2(d)给出了x=0 mm处不同阵元数量下的纵向强度曲线.对比图2(a)(b)，可以看出图2(b)横向上伪影程度要小于图2(a).根据图2(c)(d)可发现，阵元数目主要影响图像横向分辨率，随着阵元数量的增多，横向分辨率越来越高，纵向分辨率变化不大，128阵元时获得了最好的成像结果.但若继续增加阵元数目，不仅系统硬件成本增高，由于采用全矩阵采集，数据量大，将会大大增加软件的计算复杂度，因此本文选择使用128阵元的阵列换能器.

2.1.2 阵元间距分析

线性阵列换能器的阵元间距同样是关系阵列声学性能的一个基本参数，等于阵元宽度和相邻阵元间隙之和.使用128阵元的超声阵列，分别取阵元间距kerf为2λ,λ,λ/2,λ/4,λ/8，其中λ为波长，对理想点散射子(x=0,z=22)进行仿真成像.成像结果如图3所示.图3(a)(b)分别给出了kerf=λ和kerf=λ/8时的超声成像，以40 dB动态范围显示.图3(c)为z=22 mm处不同阵元间距下的横向强度曲线，图3(d)为x=0 mm处不同阵元间距下的纵向强度曲线.比较图3(a)(b)，可看出阵元间距为λ时产生的伪像要大于阵元间距为λ/8时的超声图像.从图3(c)(d)可以看出，不同阵元间距下主瓣宽度基本无差异，但旁瓣随阵元间距的增大而增大，因此阵元间距越小，成像质量就会越好，当kerf=λ/8时，成像效果最好.

2.1.3 阵元输出脉冲周期个数分析

阵元输出超声信号是一个包含n个周期的脉冲信号.在与上述仿真条件相同的情况下，使用128阵元的超声阵列，设置阵元间距kerf=λ/8，分别取脉冲信号所包含的周期个数n=1, 2, 3, 4, 5，进行全矩阵采集超声成像，结果如图4所示.图4(a)(b)分别是脉冲信号包含周期个数n=2和n=5时的超声成像，以40 dB动态范围显示.图4(c)为z=22 mm处不同周期脉冲信号下的横向强度曲线，图4(d)为x=0 mm处不同周期脉冲信号下的纵向强度曲线.对比图4(a)(b)，可明显看出n=5时散射子纵向方向拉长，即纵向分辨率降低.从图4(c)可看出，脉冲信号周期个数n对主瓣影响不大，旁瓣随着周期个数n的增大而增大.从图4(d)可看出，随着周期个数n的增大，纵向强度曲线中的主瓣宽度整体随之增大，即纵向分辨力降低.但当n=2时，主瓣宽度最窄，即纵向分辨力最好.

根据上述分析结果，Field Ⅱ仿真基本参数设置如表1所示.

表1 Field Ⅱ仿真基本参数设置

2.2 传统超声成像

利用Field Ⅱ超声仿真软件进行成像仿真，利用Field Ⅱ超声仿真软件进行成像仿真，利用Field Ⅱ超声仿真软件进行成像仿真，在不考虑噪声的前提下，首先设置2个理想点散射子模拟乳腺微钙化，如图5所示.两点相距L，并将成像区域划分为343×343个网格，利用全矩阵采集可获得128×128条回波数据.设置2点散射子距离L分别为2λ,λ,λ/2，分别进行传统B超成像，并以40 dB动态范围显示.成像结果及在z=22 mm处与B超图像相对应的超声横向强度曲线如图6所示.

传统B超成像遵循瑞利准则[34]，即当介质中的相邻2点目标之间的距离小于分辨率极限时，成像系统将无法区分这2个点目标.根据瑞利准则，超声成像系统可分辨两点目标的最小距离l[35]

(11)

式中：λ为波长；θ在成像中心位置处，换能器边缘与其垂直方向的夹角.

由式(11)，可计算出l在成像深度为22 mm处的大小约为0.34 mm.根据图6，当两点散射子分别相距2λ和λ时，传统B超成像可分辨出两点目标.但当两点相距λ/2时，由于此时两点距离小于l，因此无法区分两点.

2.3 TR-MUSIC与PC-MUSIC成像

TR-MUSIC算法突破了瑞利准则的衍射极限，具有超分辨率成像特性[33].保持Field Ⅱ基本仿真参数不变，利用TR-MUSIC算法对介质中的2个点散射子成像，分别设置两散射子距离为2λ,λ,λ/2.对超声阵列数据中的时域信号进行傅里叶变换，选定阵列换能器的中心频率f0作为工作频率，提取工作频率处的频谱值，构建阵列响应矩阵K，对阵列响应矩阵进行奇异值分解，将奇异值进行归一化，得到归一化的奇异值，如图7所示.

TR-MUSIC通过设定阈值划分信噪子空间进行成像，非零奇异值的数量代表了图像中目标点的个数.从图7可以看出2个明显大于0的奇异值，但随着两散射子之间距离的变近，其中的一个奇异值迅速减小.因此，若阈值设置不当，可能成像时不能分辨出两点目标.通常选择使用最大奇异值的10%作为阈值.

若利用多频TR-MUSIC对目标成像，首先需要选择频带宽度.以MF-TR-MUSIC成像为例，模拟设置一个理想点散射子(x=0,z=22)，选择不同频率带宽BW，分别设置为BW=1，2，3，4，5，分析不同频带宽度下，点散射子的横向与纵向强度曲线，结果如图8所示.由于是理想环境，从图8可以看出，横向强度曲线与纵向强度曲线在不同频带宽度下均差异较小，但当BW=3时，横纵分辨率可达到最好.

基于以上分析，Field Ⅱ基本参数设置不变，设置2点散射子分别相距λ和λ/2，以中心频率3.5 MHz作为工作频率，最大奇异值的10%为阈值，多频成像频带宽度为3.0 MHz，并设置衰减系数为0.3 dB/(cm·MHz)，分别进行CF-TR-MUSIC成像、MF-TR-MUSIC成像以及MF-PC-MUSIC成像，成像结果如图9所示.

由图9可以看出，相较于图6中的传统B超成像，TR-MUSIC超声成像能够有效区分开距离相近的2个点散射子，并且伪影程度较小，显著提高了传统B超成像分辨力.CF-TR-MUSIC与MF-TR-MUSIC能够准确定位，但纵向分辨率不高，并且MF-TR-MUSIC的显示结果明显大于CF-TR-MUSIC显示的结果.而PC-MUSIC则有效提高了纵向分辨率，但散射子位置发生偏移，定位不准确.并且，结果显示，2点散射子相距λ/2时的PC-MUSIC成像比2点散射子相距λ时的散射子成像强度稍微减弱.

考虑到实际应用时随机噪声对成像的影响，因此在上述仿真中引入高斯白噪声，直接添加在超声回波信号上，噪声的信噪比SNR分别设置为0、10、20 dB，如图10所示.

设置2点散射子相距λ/2，分析不同SNR条件下高斯白噪声对传统B超成像、CF-TR-MUSIC成像、MF-TR-MUSIC成像以及MF-PC-MUSIC成像分辨率的影响，结果分别如图11～14所示.可见CF-TR-MUSIC成像、MF-TR-MUSIC成像以及MF-PC-MUSIC成像与未添加噪声时的成像效果基本无差异，而传统B超成像对噪声较为敏感.

有研究发现[33]，在散射子设定的实际位置处，操作算子A在带宽内的相位随频率的增大而呈上升趋势，而成像显示位置处，操作算子A在带宽内的相位保持不变，因此由于相位差异，点散射子位置发生偏移.因此，在已知散射子实际位置时，可对其进行相位补偿.基本仿真参数不变，在深度22 mm处分别设置一个理想点散射子和2个相距λ/2的点散射子，分别进行PC-MUSIC成像与PC-MUSIC相位补偿成像，设置高斯白噪声SNR=10 dB，结果如图15所示.其中，红色圆圈表示散射子设置的实际位置.由图15(a)和(c)可以发现，散射子成像显示位置与实际位置不吻合，图15(b)和(d)为经过相位补偿后的成像，散射子成像位置与实际位置重合，实现了散射子的准确定位.

3 讨论

本文通过Field Ⅱ超声仿真软件，设置点散射子模拟乳腺微钙化，分析影响成像分辨率的相关参数，对传统B超成像、TR-MUSIC成像、PC-MUSIC成像进行仿真研究.研究结果表明，基于TR的超声成像能够区分相邻较近的点散射体，有效抑制噪声，显著提高了传统B超成像分辨率.其中，CF-TR-MUSIC成像与MF-TR-MUSIC成像均能正确定位点散射子，但纵向方向延长，导致纵向分辨率降低.PC-MUSIC成像方法在其成像函数中保留了相位信息，在横向分辨率基本不变的情况下，有效提高了TR-MUSIC的纵向分辨率.但由于忽略了阵元间的相位响应，PC-MUSIC不能准确定位点散射子.经过对PC-MUSIC进行相位补偿，提高了PC-MUSIC的定位精度，实现点散射子的准确定位.但相位补偿需预先估计散射子的位置，可先通过TR-MUSIC对散射子进行定位，然后再通过PC-MUSIC进行相位补偿.但由于采用全矩阵数据采集方式，本身数据量就比较大，而多频TR超声成像需组合多个频率下获得的超声图像，在由TR-MUSIC成像确定散射子位置后，相位补偿的PC-MUSIC需对每个散射子进行重新定位，计算非常繁琐，实时性差，不利于实际应用.

在未来的工作中，仍需对TR超声成像算法做进一步的改进，提高定位精度，降低计算复杂性，并利用真实的超声数据进行实验验证.

4 结论

1)通过Field Ⅱ仿真研究发现，TR超声成像可提高传统B超成像检测乳腺微钙化点的成像分辨率和抗噪能力；

2)TR-MUSIC成像能准确定位点散射子，但纵向分辨率低；PC-MUSIC提高了TR-MUSIC的纵向分辨率，但不能准确定位，需进行相位补偿.