基于纹影法的聚焦超声声场重建算法研究

陆彦邑，刘俏俏，赵纯亮，沈 勇，王 华

引言

纹影法作为一种非侵入式测量方式，具有直接、快速、实时等优势［1］，已经在低频声波［2］、超声波［3－5］以及聚 焦 超 声［6］的 检 测 中 得 到 了 很 多 的 应用。根据声波与光波的作用规律，声波在透明介质中发射时，用一束激光垂直穿过声波，介质密度的改变会引起折射率发生改变，形成稳定的相位光栅，导致光波相位变化［7］。据此，利用纹影成像技术便可得到相应的声场纹影图像，从而重建出声场声压分布。然而现有的纹影法并不能直接利用纹影图像重建得到声压分布图像，原因是纹影系统中所用的空间滤波器多为刀口［8］、针孔［9］、或档板［10］等，这些空间滤波技术得到的图像光强信息与声压梯度有关［11］，重建后的图像为声压梯度图像，并不是声压分布图像。而利用Zernike相衬技术［12］则可以将光强与声压的梯度关系转变为光强与声压的关系，从而对纹影图像进行重建得到声压分布图像，具有简单、直接、可靠的特点。基于此，本文的纹影系统利用相衬滤波器得到声场纹影图像，经反投影重建算法从光强图像重建出声压分布图像。首先研究凹球壳聚焦超声换能器的声场分布特性；而后利用激光纹影系统采集声场的实时影像，并借鉴CT图像的反投影重建算法对聚焦超声声场进行重建，得到声场声压分布。

1 纹影系统

1.1 实验装置

纹影法是研究透明介质中声场的一种常见光学方法，实验装置如图1所示［13］。Laser为He－Ne激光仪；EL为扩束透镜，L1为准直透镜，L2、L3分别为变换透镜和成像透镜；SF为空间相位滤波器；FT为聚焦超声换能器，固定于装有脱气水的水槽中；CCD为摄像机，用于获取声场图像。

图1 纹影系统实验装置图Fig.1 Experimental setup of schlieren system

1.2 实验原理及方法

根据声光效应，透明介质中，超声波的传播引起介质密度及折射率变化，形成稳定的相位光栅，引起光线偏转。当光线在均匀介质中传播的时候，介质中折射率任一点都一样，光线沿直线传播。而当超声波作用于均匀介质时，声压导致介质密度变化，从而引起介质折射率改变。此时，光线在传播路径上相位发生改变，不再沿直线传播。

用激光器投射出激光，垂直穿过超声束。当无超声作用时，平行光聚焦在SF上；当发射超声时，引起传播路径上介质折射率改变，形成相位光栅，平行光经相位光栅衍射，衍射光通过SF滤波，由CCD获取声场实时影像。

激光相位改变φL（t）与声压改变p（t）的关系为［14］：

式中：μOP为介质压光系数，水中取1.51×10－10m2／N；λL为 激 光 波 长，He－Ne 激 光 波 长 为632.8nm；L为激光穿过声场时的有效作用宽度。本文所用聚焦超声换能器为重庆海扶医疗科技股份有限公司生产的CZF型超声波妇科治疗仪，经计算，该换能器L数量级为10－4m，声压数量级可达106Pa。计算可得φL（t）≪1弧度，因此该过程可以用空间滤波原理中的Zernike相衬技术来表示，即不可见的相位分布可以转化为可见的光强分布［12］：

φ（x，y）为经过声波调制后的光波相位分布［15］：

式中：P（x，y，z，t）为声压分布；ke0为真空波数；pop为压电光学常数。根据该式，相位分布为声压在激光传播路径上的线积分。可见，纹影系统得到的纹影图像光强也可以表示为声压在传播路径上的线积分：

该过程与CT图像的投影过程类似。

将超声换能器固定于装满脱气水的水槽中，打开激光仪，驱动换能器发射超声波。此时，便能在屏幕上看到超声场的实时显示。依次改变换能器驱动功率大小，超声场亮度也随之发生变化。

2 超声场图像重建算法

2.1 CT图像重建理论

根据著名的Lambert－Beer定律，可得到非均匀介质中，X射线投射后的投影这p为

式中：I0为入射X射线强度；I为出射X射线强度；衰减系数μ（l）是随路径l连续变化的函数。

对比（4）式与（5）式，可以看出，纹影系统与CT扫描系统具有相似的数学表达。因此，可以将光强分布看作是该方向上的声场投影值，从而借助CT图像的重建方法来重建超声场声压分布。

2.2 直接反投影法重建声场

直接反投影重建算法［16］的基本内容是，在对某体层一个方向的扫描完成后，以得到的投影为灰度值沿着扫描路径经过的体素回抹到体素对应的像素上，改变方向后的多次扫描形成多次回抹，同一像素上多次回抹的灰度值累加即完成图像重建。

如图2所示，可以把激光发射装置看作由n×n个小发射器组成，由此在屏幕上可得到n×n个像素的光强分布图像。而实际上，根据CT重建算法，我们只需要其中一列光强分布即可重建出这列分布的激光穿过的声压截面。将该列数据沿着激光路径均匀地回抹到要重建的声压截面上，旋转激光器，可得到多个角度的光强分布，将所有角度的同一位置光强分布都回抹到重建图像上，再除以角度数，即可得到相应的声压截面分布。

图2 重建算法原理图Fig.2 Schematic of reconstruction algorithm

本文中，由于无法旋转纹影系统以获得多个角度的投影值，只能得到如图2所示的一个角度下的投影值。因此，需要寻找一种模型来获得其他各个角度与该角度下的投影值的关系，从而得到多个角度的投影值重建声场。基于声场原本具有的声压分布，以此为模型即可得到一个角度与其他各个角度投影值的关系。方法如下：1）通过瑞利积分得到与纹影系统对应的声压截面的归一化声场分布图［17］；2）将该图在图2所示方向上定为第1个方向（水平方向），把该方向图像的行向量相加得到第1列数据（初始数据）；3）顺时针旋转图像，再在水平方向重复步骤2），依次旋转359°，得到359列数据；4）用步骤3）得到的359列数据分别除以第1列数据即可得到水平方向投影值与旋转后投影值的关系。

用纹影系统采集一幅声场光强图，取图像中最大值那一列数据，该列数据即对应图2中水平位置投影所得投影值。将该列数据分别乘以上述第4步里所得投影关系数据，便得到了360个角度下的所有投影数据。将所有投影数据回抹为声压对应的像素值，再除以角度数，即可完成声场重建。算法流程图如图3所示。

图3 重建算法流程图Fig.3 Flow chart of reconstruction algorithm

3 结果与讨论

为了确定仿真图像精度，先对CCD采集的纹影图像的像素大小进行标定，方法为获取已知直径为1mm的短棒的纹影图像，计算其直径所占像素值大小。经标定，本文所用CCD相机采集的图像20个像素宽度为1mm。采集一幅声场纹影图像，由于图像中感兴趣区域为声焦域附近，因此截取包含焦域在内的声场轴向4mm、横向4mm范围做重建，该图像大小记为M×N像素值。

3.1 仿真结果

对声场进行仿真时，同样只计算包含焦域在内的声场轴向4mm、横向4mm范围，也使仿真图像大小为M×N像素值以便于之后进行重建。使用Matlab进行仿真实验，所得理论归一化声压分布的x－o－y 面如图4（a）所示；根据（4）式在z轴上对理论声场积分可得纹影图像理论光强分布，如图4（b）所示；用上文所述算法得到的声压分布重建图像如图4（c）所示。对比可知，重建算法能成功地重建出声场声压分布。

图4 仿真图像Fig.4 Simulated images

3.2 不同电功率下纹影图像

本文采用重庆海扶医疗科技股份有限公司生产的CZF型超声波妇科治疗仪作为聚焦超声发射装置，其换能器为中间开孔型凹球壳聚焦超声换能器，频率为9.1MHz，焦距为12mm。调节换能器驱动电压与电流可得不同电功率下的纹影图像，如图5所示，下方为换能器。可以看出，随着功率的增大声场图像亮度变大，声场强度增强，焦域大小也随之变大，焦域形态更加清晰。用其余驱动功率下纹影图像减去功率为0时的背景图像即可重建。分析图像后可知，焦域呈椭球形，宽度为10－4m数量级，长度为10－4m～10－3m数量级。

3.3 重建结果

图6为利用直接反投影法重建所得用灰度级表示的声压分布。由于功率为7W时的声场较弱，纹影图像与减去背景后的图像中均无明显声焦域形态，因此不对其进行重建。

图5 不同驱动电压下声场纹影图像Fig.5 Schiliren images of different driving voltages

图6 不同电压下的重建结果图Fig.6 Results of reconstruction with different driving voltages

根据其余结果可以看出，基于理论值得到的重建声场与纹影图像吻合度较高。进一步分析其声轴与横向声压分布，如图7所示，可见随着电功率增大声场变强，与纹影图像变化趋势吻合。

图7 声压分布Fig.7 Acoustic pressure distributions with different driving voltages

以声压下降3dB为界限，可得图7中各图的声焦域轴向与横向像素值大小，继而根据标定结果得到其具体长度与宽度。表1将结果与理论图结果进行对比发现，12W时横向最接近理论声场，30W时轴向最接近理论声场。表1表明重建声场轴向与横向比值要比理论声场轴向与横向比值要小，也就是说，重建后声场横向变宽了。这可能有以下几个原因：1）纹影图像本身轴向与横向比值就要比理论声场轴向与横向比值要小，所以基于纹影图像的重建结果也如此；2）换能器工艺设计和环境影响使得声场本身就与理论声场存在差别，这也是进行声场测量的主要原因；3）本文使用直接反投影算法，使用此算法重建CT图像时会出现星状伪迹，因此从图7可看出重建声场也不可避免地出现了伪迹。

表1 不同电压下声焦域大小与理论值对比结果Table 1 Focal region size of different voltages compared with the theoretical value

4 结论

利用光学纹影法检测聚焦超声声场，基于CT成像系统与纹影系统的相似性，本文提出了一种借鉴CT图像重建算法来重建超声场声压分布图像的方法。通过实验得到了不同电功率下的聚焦超声声压分布图像，结果图像与理论图像的对比论证了该方法的可行性。本研究为纹影系统用于聚焦声场检测提供了实验依据，也为声焦域分析提供了理论基础。

［1］ Jiang Xueping，Chen Xi，Qian Menglu.Theoretical and experimental investigation of imaging the acoustic fields by schlieren techniques［J］.Technical A－coustics，2011，30（05）：1－4.

姜学平，程茜，钱梦騄.纹影法对声场成像的理论和实验研究［J］.声学技术，2011，30（05）：1－4.

［2］ Torras－Rosell A，Barrera－Figuera S，Jacobsen F.Sound field reconstruction using acousto－optic tomography［J］.Acoustical Society of America，2012，131（5）：3786－3793.

［3］ Kudo N，Sanbonmatsu Y，Shimizu K.Microscopic visualization of high－frequency ultrasound fields using a new method of schlieren photography［C］.US：IEEE，2010：829－832.

［4］ Chinnery P A，Humphrey V F，Beckett C.The schlieren image of two－dimensional ultrasonic fields and cavity resonances［J］.Acoustical Society of A－merica，1997，101（1）：250－256.

［5］ Zhu Weimin.A study on measurement of ultrasound field based on optical technology［J］.Measurement Technique，2013，5：16－18.

朱卫民.一种基于光学技术的超声场测量方法研究［J］.计量技术，2013，5：16－18.

［6］ Remenieras J P，Matar O B，Calle S，et al.Acoustic pressure measurement by acousto－optic tomography［C］.US：IEEE，2001：505－508.

［7］ Settles G S.Schlieren and shadowgraph techniques：visualizing phenonmena in transparent media［M］.New York：Springer，2001：25－28.

［8］ Moller D，Degen N，Dual J.Schlieren visualization of ultrasonic standing waves in mm－sized chambers for ultrason.ic particle manipulation［J］.Journal of Nanobiotechnology，2013，11（1）：1－5.

［9］ Unverzagt C，Olfert S，Henning B.A new method of spatial filtering for schlieren visualization of ultrasound wave fields［J］.Physics Procedia，2010，3（1）：935－942.

［10］Zhu Guozhen，Lu Kean，Fu Deyong，et al.Experi－ments on two kinds of threshold for the acoustic pressure gradient of a schlieren system［J］.Measurement Science and Technology，2002，13（4）：483－487.

［11］Brownlee C，Pegoraro V，Shankar S，et al.Physically－based interactive flow visualization based on schlieren and interferometry experimental techniques［J］.IEEE Transactions on Visualization and Computer Graphics，2011，17（11）：1574－1586.

［12］Su Xianyu，Li Jitao，Cao Yiping，et al.Information optics［M］.2nd ed.Beijing：Science Press，2011：207－208.

苏显渝，李继陶，曹益平，等.信息光学［M］.2版.北京：科学出版社，2011，207－208.

［13］Shan Zijuan，Wang Dingxing，Li Zhengzhi.Properties of a laser schlieren system［J］.Acta Optica Sinica，1984，4（10）：880－886.

单子娟，王定兴，李正直.一种激光纹影仪的光学特性［J］.光学学报，1984，4（10）：880－886.

［14］Harvey G，Gachagan A，Mcnab A.Ultrasonic field measurement in test cells combining the acousto－optic effect，laser interferometry ＆tomography［J］.Proc.IEEE Ultrason.Symp，2004，2：1038－1041.

［15］Holm A，Persson H W，Lindstrom K.Optical diffraction tomography of ultrasonic fields with Algebraic Reconstruction Techniques［J］.Proc.IEEE Ultrason.Symp，1990，2：685－688.

［16］Huang Liyu.Basic principle of medical imaging［M］.Beijing：Publishing House of Electronics Industry，2009：99－102.

黄力宇.医学成像的基本原理［M］.北京：电子工业出版社，2009：99－102.

［17］Zhang Dejun.High intensity focused ultrasound transducer［J］.Chinese Journal of Ultrasound Diagnosis，2000，1（2）：1－4.

张德俊.高强度聚焦超声换能器［J］.中国超声诊断杂志，2000，1（2）：1－4.