超声系统前端理论与模拟仿真-续
作者:蒋志强
??本人同意他人对我的文章引用,但请在引用时注明出处,谢谢.作者:蒋志强
前言
????????近期整理了一下彩超前端及波束合成相关的内容,很早以前已经有过一次,这次把其它的内容总结一下,又做了不少仿真和实验,以便大家能直觀的理解,希望对朋友们有所帮助。
内容覆盖
? ? ? ? 发射,声场,声场传播,声场交互,Grating lobe, Side lobe,PSF,回波信号接收,接收聚焦,孔径与变迹,MLA,MLT,RTB回溯波束合成,RTB问题与修正,Diversed?Wave(发散波)Plane Wave(平面波)?Focused Wave(聚焦波),HFPS(高帧率成像),合成孔径(Synthetic Aperture),大孔径发射,小孔径发射,系统设计与工程考量,其它...
导读说明
本文以叙述说明的方式科普相关知识,涉及到相关知识点的段落时,留意前面用蓝色加粗字体提示。为了形象直观的描述问题,我做了大量的模拟仿真实验,目标读者是相关实验室同学或新进入彩超领域的同事理解,大家先记住结论,很多关键细节的内容我无法在文中展开。大家在一线工程实践中,自然就明白了。
如果你是彩超系统工程师(不管是偏FPGA或偏SW的系统工程师),欢迎你留言提出各种意见与建议。
正文内容
发射,声场
超声探头本身包含数以百计的element,我们在发射声波时可以逐个控制阵元的发射延时和电压,从而控制声场。概念上,波束合成包含两部分,一个是发射,另一个是接收;?
我们可以抽象为以下传递函数? Beam = Func_rx(Interaction(Func_tx(signal),tissue))
Func_tx是发射部分,它直接决定了声场的情况。我仿真了一组声场,我们先直观的感受一下,声压在组织中的传播所形成的声场,见下面的视频
我的声场仿真动画
如果视频点不开,就看图吧,从左看到右,靠脑补动起来
我的声场仿真不动的画
?这组仿真,是一个64阵元的小探头(按常规的相控阵配置的),声场区域长宽都是4厘米,发射聚焦在2厘米的深度,2个Cycle的激励,2MHz的发射频率。实际上声场本身是机械波的传播过程,在传播过程中会与组织发生交互,在各个位置上不断发生反射,折射,散射,衍射。上面大家看到的声场中,没有任何反射等交互,是因为我们仿真的是在水中传播的情况,没有任何的物质。但我们可以看到,由于控制了发射延迟,让所有的阵元的声波同时到达2厘米处,在2厘米处波束很窄很集中,在焦点前后逐渐远离发射焦点的地方,波束逐渐变宽。
Grating Lobe
在上面的例子中,声场传播需要一段时间,为了简单便捷的表示声场,我们一般通过在空间上计算beam profile来描述波束。现在我们来仿真一个,同样的相控阵探头,3MHz的发射,但是不再是正向下,而是向左偏转45度。
?我们向左偏转45度的beam profile仿真
为了便于直观查看,我调整到了一个合适的动态范围。细心的朋友能留意到,波束右侧有一个小的能量分布,那部分杂波能量我们叫为Grating Lobe(栅瓣),它会带来grating lobe artifact。如我们仿真的情况,你扫查的是左侧45度,而如果grating lobe处有组织的强反射,这部分信号将混入到我们后面的接收波束合成中,导致右侧栅瓣的组织出现在左侧45度的地方。grating lobe与发射频率,探头阵元的间距相关。一般情况下,只要它偏得比较远,在成像区域外,我们就可以不去理会它,否则你要考虑一下,这个特定的探头换个发射频率或者缩小成像区域。
一副二维图通常需要 重复我们仿真的那种发射若干次(下图左绿线所示),才能完成一副成像,如下图绿线所示。
?别人的相控阵图像
我们模拟一下上图相控阵探头的二维扫查,直观感受声场能量分布变化;
我们在2D扫查中的beam profile模拟仿真动画
如果视频点不开,还是看图吧,靠脑补动起来
?仿真完整一副2D相控扫查beam profile不动的画
线阵探头的情况也类似,只是线阵扫查时,发射一般不偏转,或者偏转较小(如空间复合,穿刺增强时偏转较大是比较特殊的情况,穿刺增强可能会引入较大grating lobe伪影)。
下面,我仿真一个128晶元的线阵,3.8厘米宽度的线阵探头,进行完整一帧图像扫查的声场情况;
?顶部白线区域表示探头孔径,绿线位置表示发射时,实际激励晶元的位置。
?注意对比一下相控阵探头偏转时的grating lobe,线阵就会小不少。
聚焦波/平面波/发散波?
目前为止我们的发射都是聚焦在FOV(Field of View)内,实际上我们的发射聚焦至少可以有三种以上:
1. 焦点在FOV内:聚焦波(Focused Wave)
2. 焦点在无穷远(即无焦点):平面波(Plane Wave)
3. 焦点在探头后面:发射波(Diversed?Wave)
4. 其它
我在相同的探头,相同的孔径,用相同的激励波形,以不同聚焦方式,仿真对比一下相同的FOV下的声场,如下图
?几种发射,用在成像上各有不同的特性。我先暂时用聚焦波为大家说明,因为目前仍有很多产品使用聚焦波,稍后在波束合成时,再用平面波与发散波进行一下对比分析。
By now,我们控制了基本的发射,有了声场,我们来看看接收后的情况。我们能够接收到信号,来源于扫查目标中有声阻抗的变化,如果一个无穷大的水槽,我们将什么信号都收不到,所以我们模拟几个靶点。
我们这次模拟的是一把128阵元,宽度为3.8厘米的线阵,7MHz的发射频率,中央正下放配置几个理想靶点(体积无限小,足够的强幅度,类似于数学中的delta function),对于接收到的信号,逐点计算延迟曲线进行接收聚焦,如下图所示
扫查示意图(左)模拟动态接收聚焦成像(右)
声场传播,声场交互
当声场中有了反射点时,声场就会与这些反射点发生交互,常见的反射,衍射,折射等。下面我用一个192晶元的线阵外加一个反射点模拟一下这个交互过程,用红色标注处是反射点位置。大家留意声场的交互过程
回波信号接收
?上面的交互过程的反射信号,被探头接收转换为电信号后,就是我们说的通道数据。下图是上组实验中接收到的信号;
怕大家看不清楚,红色框内,我拉大了显示; 每个通道的信号,基本长得和激励发射时的波形一样。实际上这是不对的,因为彩超系统不是线性系统,且声音在传播过程中也是非线性的,波形也会逐步发生畸变,所以我们的仿真不是太真,不过基本够用。你一认真,就输了。
Side Lobe
我们再回到靶点体模图,注意到有点靶点两侧有小的震荡,这是我故意配置动态范围比较夸张,为了说明Side?Lobe的影响。其实各个靶点都有,只是明显层度而已,在这个动态范围下比较明显。如果不想出现得太明显的Side?Lobe改改配置,就不那么明显了。
接收聚焦,孔径
在发射端我只给了一个延迟曲线和焦点,但接收我们可以根据深度变化,每个点的延迟曲线都可以单独计算,所以叫动态接收聚焦(Anyway,在工程上常常为了省资源,分成很多小段,一小段使用一组接收延迟曲线),由于越远的地方信号越弱,越近的地方信号越强(但杂波多),我们一般保持成像的深度与使用的接收阵元的个数(探头孔径)的比值为一个常数,即F#,焦径比,焦点与孔径大小的比值。
距离越远,成像需要开的孔径越大。大家还记得2019年与2021年EHT观察距离我们5500万光年的M87星系中心黑洞时,雷达的理论抽象孔径虽然达到了1.2万公里(即地球直径),但F#仍高达4.3X10^16,看不清是正常的,要知道我们上面成像仿真中F#=1,如下图
成像系统F#比较
?回到我们的超声成像,右侧图像基本反映了我们的实际扫查结构,但大家请留意:除了发射焦点深度(2.5cm)上聚焦比较理想,其它深度都比较差,我们预期远离焦点位置会差一些,但结果也太差了。主要原因是,我们还没有做遍迹Apodization;
变迹(Apodization)
所谓Apodization是指探头在发射/接收信号时,给与不同阵元不同的权重,从而使得波束合成的效果更好。根本原因是,阵元具有指向性,垂直阵元方向上信号不管发射还是接收端SNR最高,于是我们赋予更高的权重。我们可以在发射做Apo,也可以在接收做Apo,或者两者都做,我在Rx时对比一下有无Apo的区别,见下图
接收动态聚集无Apo(左)接收动态聚集有Apo(右)
怕看不清楚,放大些,怼近看
?放大看 接收动态聚集无Apo(左)接收动态聚集有Apo(右)
不仅发射,接收具有指向性,声音作为一种波,在与组织交互过程中也存在指向性,而并非均匀的球面波,而是与发射波阵面的法线夹角成负相关。?
Beam = Func_rx(Interaction(Func_tx(signal),tissue))
我们再回看波束的抽象函数,三个函数都具有方向性,给予方向上的权重,从本质上符合超声系统的成像原理;
点扩散函数(PSF)
Apodization实际是以Main Lobe稍微变宽为代价,换来一级旁瓣大幅压制,but anyway还是一笔划算的compromise。上面以理想靶点为输入,得到的图像,我们常称为PSF(Point Spread Function)。PSF是泛指任意成像系统,当输入为一个delta function时,它的输出结果成像所对应的抽象函数,它描述了成像系统本身,如下图所示
PSF示意
在超声成像系统的发展过程中,我们总希望成像更快更好。在更快的方向上,如果整个二维成像平面的线数一定的情况下,我们需要等待声音信号回来,这是物理限制,于是发展出了MLA,MLT等方法,下面我们通过模拟仿真来介绍一下;
多线接收(MLA)
MLA即Multiple Line Acquisition,之前我们是发射一次,接收一条线,MLA是发射一次,接收多条线。But how, MLA的做法非常简单粗暴,只发射一次,但按照新的位置进行接收布线,并按新的布线位置计算接收延迟曲线。以MLA=4为例,如下图所示
? ? ? 左1:常规接收布线位置与beam profile关系;? ? ? ? ? 左2:常规接收布线位置与声场关系;
右1:MLA4接收布线位置与beam profile关系;????????右2:MLA4接收布线位置与声场关系;
对于一次发射后,接收到的信号放在FIFO里,现在使用两种不同的策略的接收Delay Curve进行波束合成,左侧是常规合成一条线的方式,右侧是合成四条线的方式,如下图所示;
?左图:一条线的接收延迟曲线? ? ? ? ? ? ? ? ? ? ? ? 右图:MLA4时的接收延迟曲线
MLA由于忽略了声场能量的不同,以及实际发射声场波阵面相位的差异,成像的效果可以预料一定会比单线要差。所以我们再强调一下,MLA的图像质量退化,来源于:
1.声场横向位置的能量差异;
2.声场波阵面的行程相位差异;
回头我们再来从这两方面修正它。我们先看看MLA图像质量下降会有多大呢,我们模拟仿真一下,让接收布线相同情况下,一组波束合成使用单发射单接收一条线,另一组单发射多接收(MLA=4),如下图所示
左图:MLA=1,Tx=100,Rx=100? ? ? ? ? ? ? ? ? ? ? ? 右图:MLA=4,Tx=25,Rx=100?
我们再怼近看,在MLA=4时,靶点会更宽一些?Is that true ??
?我先给结论:
1. 图像品质会下降;
2.靶点宽窄会变化,取决于所处声场的位置,有可能窄,有可能宽;
为了直观的描述第二个结论,我把探头平移一下,模拟手持探头向左平移,看一下它的变化
Ultrasound MLA4 movement
平移探头动的画
如果视频点不开,还是看下面不动画,自己脑补动起来
上排:NO MLA(MLA=1)? ? ? ? ? ? ? ? ? ? 下排:MLA=4
?我们可以看到MLA=4,在相对运动时,靶点明显忽宽忽窄,扭来扭去,这就是我上面提到的第二点,也经常叫做MLA block?artifact;而在MLA=1,在相对运动时,靶点时稳定的(除了在探头边上,由于接收孔径不对称造成少量轻微扭曲);
对于这类MLA artifact,一般都会做一些Interp或STB进行补偿,但补偿的结果是以牺牲横向分辨率换取artifact不可见。这个我应该也模拟一下,因为我很懒,就不做了。
只有靶点,不方便我们系统对比MLA的实际影响,我做一个完整的体模比较,有靶点,组织,及对比分辨率靶球,如下图所示
?????左图:常规体模 Tx=160,MLA=1,? Rx=160? ? ? ? 右图:常规体模 Tx=40,MLA=4,? Rx=160?
?从这个对比,大家就很好理解为什么多MLA的图像下降叫做block artifact,非常形象。
发射回溯波束合成(RTB)
下面我对比一下效果好一些的修正MLA artifact的方式:Recursive Transmit Beamformation(RTB),这种方式能帮助在做多MLA时,不牺牲横向分辨率。RTB是95年左右Freeman与O'Donell老师和Pai-Chi Li?老师最早在论文里提出的思路,后面两位老师做彩超的朋友应该都比较熟悉。大概二十年前(2003还是2004有点记不清了),Jensen老师(学超声的朋友应该都很了解)首先在实验室里,搭建了以此思路并完成合成孔径的实时平台(当时GPU还没有发展起来,FPGA也很难做到,印象中是50块Xeon CPU并连的工程方案,现在想起来还是要给他们实验室点个大大的赞)。
RTB的原理很简单,由于多MLA问题本质是忽略了不同接收线在声场中波阵面的不同相位,所以解决方式就是不要忽略它。我们还是以MLA=4为例(其它情况类似),我画个示意图如下
?我画的个示意图
黑色弧线表示某时刻声场波阵面的同相位处,我们这一次发射,准备接收信号以后合成四条扫描线,即MLA=4,图中标注L1-L4。这四条接收线,使用同样的接收延迟曲线Delay Curve 1-4。但严格来讲,该时刻延迟曲线,只对图中的红色点是正确的,因为波阵面正处于红点的中心。只有当波阵面处于黄色点中心时,该延迟曲线才应该去合成L2,L3上的黄色数据点,同理只有当波阵面处于蓝色点中心时,该延迟曲线才应该用作合成L1,L4上的蓝色数据点。我们应该在不同时刻的波阵面下去使用该延迟曲线,我再画个示意图如下
我再画个示意图
如图所示,我们需要让波阵面回退到黄色位置时,去合成那L2,L3的那两个点,需要波阵面回退到蓝色位置时,去合成L1,L4的那两个点。我们只要在接收的通道数据里提早或延后一些,就等效于波阵面不同的时刻,实现了波阵面回退或前移。这就匹配RTB名字的意思了,即回溯发射波束合成。
为了对比更明显一些,我以MLA=8的情况进行比较,如下图所示
(PS:下面这个图有问题,但我换电脑时这组数据找不到,这组实验我不打算去改正了。大家正好可以从结果上来找找错。文章后面部分,我会换新数据再做)
MLA=8 comparison between without?RTB and with RTB?
?看上去似乎差别不大,但留意焦点位置,在2.5 cm深度位置的变化。在越远离焦点的位置,波阵面会越来越接近平面波,所以回溯的差别不大,在越接近焦点的地方,差别会越大。这个里面,有一些特殊的处理,我就不展开了。我做成视频,切换着对比,我们就看得更明显一些,如下视频
MLA8 with and without RTB
MLA=8?without?RTB and with RTB 切换
由于聚焦波在焦点附近无声区的存在,如果保持接收布线密度不变的情况下进一步加大MLA数量,则会让X形的黑色区域更大更明显。当然我们可以在后面做合成孔径(Synthetic Aperture)时,系统工程师费点心(其实也不是太恼火)做一些特殊处理,把这个问题修正一下。
当然这是聚焦波特有的需要特别处理的地方,平面波和发散波由于声场内没有聚焦,不存在无声区,做高倍数的MLA时,工程上更简单点。
HFPS(高帧率成像), 平面波,发散波
(由于换新电脑,之前一些资料找不到了,换个更炫酷的体模来仿真)
下面我做个极端例子,使用128晶元线阵探头,平面波发射,只做一次发射打开全部晶元,波束合成时,接收MLA=512,覆盖整个成像区域,情况如下?
?相同的情况,换作发散波,单次发射,打开全部晶元,MLA=512,覆盖全区域,结果如下
左:体模? ? ? ? 中:单次全孔径平面波? ? ? ? 右:单次全孔径发散波
?扩散波与平面单次发射,全区域成像结果比较类似。我们可以发现,这两种发射由于成像区域内没有焦点,也就不存在无声区,即使我做到512的MLA覆盖4厘米的宽度,也不会出现暗黑的无声区域。另外,大家应该留意到MLA的BLOCK格子也不会存在。因为声压全场统一,不会由于不同的MLA线处声压强度,导致MLA线实际位置横向偏移,同时这种极端扫查情况下,帧速率(FPS)可以做到每秒几千甚至上万。可以把这样的优势用在需要高帧速的场景,更精细的运动追踪,时间解析度更高的造影等等。
但缺点也很明显,grating lobe伪影明显,信噪比低,分辨率差等。在实际系统设计上,我们肯定使用更小的孔径,更少的MLA,多次发射,降低FPS,来平衡图像品质。?
合成孔径(Synthetic Aperture)
与平面波相比,扩散波非常类似,唯一的不同是非全孔径发射时,在后续的合成孔径处理中,平面波需要扭转方向,来提供不同波阵面的方面,而扩散波不必如此,它天生确保了波阵面的方向不同,有利于合成孔径。
在具体说合成孔径之前,大家先看效果,对比一下,有个直观印象先。下边右图是进行23组合成孔径的效果,偏转角度从-15度至15度,发射是全孔径激励的平面波的成像结果,中间是单次平面波无偏转,全孔径激励的平面波的成像结果,左图是原始体模;
合成孔径对超声图像的提升是显而易见的。为了更直观的表达,下面的动画展示了进行31次合成孔径的处理过程中带来的变化,由于旋转了31个角度,波阵面各有不同,所以每个角度实际上也进行了类似聚焦波的RTB处理,效果如下:
Synthetic Aperture即合成孔径,是成像系统中一个比较广泛应用的成像技巧;虽然我们这里将用到超声成像里,但实际在雷达,射电望远镜,声呐等成像系统中广泛应用;
合成孔径示意
在我们成像系统中将发射或接收孔径不同的N组信号(或者发射与接收都不同的N组信号),进行合成,用以最终成像,就可以被称为合成孔径;合成孔径最大的优势在提升信噪比,付出的代价是牺牲时间分辨率,受组织运动影响明显;另外由于靶点反射信号形态与波阵面的法向量相关,几个不同角度的信号合成,本质上能起到spatial compounding的效果,或者说spatial compounding也可以归为合成孔径的范畴,这将对收缩靶点的旁瓣,提升信噪比有良好的帮助;
合成孔径在彩超产品上的实际情况
在常规成像过程中,一般不会像上面实验那样打开全孔径,一般都是小孔径多次重叠发射,较小的MLA(大部分情况不超过32),合适的合成次数(部分情况不超过16),降低FPS(大部分场景用不上几千帧的FPS,除了几个特别应用,例如shear wave追踪,运动可视化或HFPS的造影),以减小side lobe以及运动对合成孔径带来的影响。
做在产品中,一般有两条技术路线:
1. FPGA工程路线
? ? ? ? 优势:减小了软件端的压力,避免了PCI-E传输的压力和成本。在通过WIFI传输的手持类系统中,有着得天独厚的低成本低数据传输压力的独特优势;
? ? ? ? 缺点:可扩展性,可调试性,后期高端应用的扩展性弱;
2. SW工程路线
? ? ? ? 优势:可扩展性,可调试性,高端功能的扩展(如通道内运动追踪补偿的波束合成,自适应计算的波束合成),打开了后期中高端产品的性能上限;
? ? ? ? 缺点:成本一般高于FPGA方案,对数据传输接口数据率有较高要求,对高频超高频探头的应用不友好(例如CMUT/PMUT等高频/超高频探头);
从长远来看,两种技术路线各有长短,我们需要用发展的眼光来看待。
FPGA的扩展性,调试性会随着技术的发展而进步,有一个接近软件端的可扩展可调试性也不是没可能;而对软件端而言,随着电脑主机的技术发展,导致性能逐年提升的同时成本每年都在下降,若有一天综合成本接近FPGA方案的低成本,我也不会惊讶;
对自主研发的厂商来说,团队同时拥有两套路线,将大大拓展其未来的商业利益;
结束语
行业展望
这篇概述文章,在一年多以前我就写得差不多了,由于时间原因完全整理好,就一直没有发出来。
几个月前,离开了上一家工作的超声团队,正好有时间来整理一下资料给发出来了。希望能够帮到刚毕业或还在实验室学习的彩超领域的同学们。
我个人对国内彩超自主研发团队是充满信心的。从2000年初到现在,已经有一批自主研发的工程师成长了近二十年,逐步形成了研发梯队。在彩超领域早已不再是仰望国外的年代了,甚至在彩超领域的个别项目上已经局部领先欧美日的研发。
彩超这个复合学科,高技术门槛的领域,会不会像手机,家电,太阳能面板,高压输电等领域一样,最终由中国人全部自主且领先?这个答卷,需要年轻一代自主研发人员与我们这批老腊肉来共同作答。
至少,我对年轻一代是充满信心的!
忙点自己的开发
工作这么些年,离职后也闲不下来,我发现自己一个人在家开发点东西,比在大厂内工作效率高了至少一个数量级,每天都思如泉涌。
我打算设计开发一套彩超系统链路工程化的软件框架,在设计架构上希望至少不弱于行业内大厂的水准(例如不弱于GE挪威团队Steen N. Eric几人开发设计的cSound框架)。Flag似乎立得有点高,但不push自己一把,自己也不知道自己在系统与工程两方面真正的limit在哪里,荒废了难得的“空闲”时间。
要保持在心流状态进行研发,也必须处在自己能力边缘,高标准是必要条件。运气使然,我应该是国内唯一一个全面深入研究过cSound框架与工程实现所有细节,并在工程上进行大规模拆解重组的工程师。它的设计优势不只是在波束合成这一小块,而是在整体架构链路的设计和工程架构上,而它的短板如果没有在研发的第一线做过几次自主研发框架的搭建,踩过几次坑,也不可能想到,更无法谈超越了。我恰好同时具备以上条件的工程师,只能说是运气好?
构思了一阵,目前开发已经开始一段时间了,研发团队暂时就我一个人。GPT和Bard也参与讨论分析很多,也算一个吧(GPT比Bard更靠谱一些)。AI在综合性系统性的问题,设计层面架构层面的问题? 还比较白痴,但特定技术细节方面的专业深度让我刮目相看,给我很多启发,在调试中找错的能力也不错(表扬一下这免费的优质劳动力)
但千头万绪太多事要做,感受一下我日常开发中的桌面混乱状态
好的方面是,整个可运行的架构缓慢但稳健的逐步成型了;不好的方面是,我得一会儿做理论分析,一会儿去模拟仿真测试一下,一会儿思考软件框架写设计文档,一会儿做具体代码调试BUG,一会儿设计定义资源配置文件,反复切换任务经常自己把自己搞蒙了。
但我已经慢慢适应了在这种混乱中乱而有序的推进开发,我把这叫做有序的混沌开发(WTF)。当这一套框架开发到一定阶段时,如果不出意外,我会把它可运行的软件包放出来与大家共享。
开发到完成基本链路功能算结束(常规聚焦/发散/平面波的合成孔径肯定要完成,基本B/C/D串联也需要完成,CPU/GPU/AVX配置都要能实时运行,效率要满足正常产品需求),余下的肯定要和团队一起推进。毕竟要完成整个体系是一个团队的工作量。独行快,但众行远。
那时正好竞业期也结束了,重新出发,岂不完美?
在2023年底,预祝大家来年好运,也祝自己的开发在新年一切顺利!
本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。 如若内容造成侵权/违法违规/事实不符,请联系我的编程经验分享网邮箱:veading@qq.com进行投诉反馈,一经查实,立即删除!