基于相位衬度成像的X射线生物特征识别仪

许明伟 薛艳玲 陈荣昌 李 可 杜国浩邓 彪 谢红兰 肖体乔

1（中国科学院上海应用物理研究所 上海201800）

2（上海光源科学中心中国科学院上海高等研究院 上海201204）

3（中国科学院大学 北京100049）

国家生物安全对于濒危物种的保护、防止外来生物入侵和维持生态系统的稳定至关重要［1-2］，而海关是维护生物安全的第一道防线。X射线因具有强穿透性的特点被海关部门广泛用于安全检测工作中［3］，2008年诞生的直线加速器Linatron Mi6能产生3.5 MeV和6 MeV两种能量的X射线对货物进行透视检测［4］；为了穿透厚度超过300 mm的集装箱钢板，2012年诞生了能产生6 MeV和9 MeV高能X射线的IDK-6/9 MeV线性电子加速器［5］；2017年又出现了可产生光子能量5 MeV、6 MeV和9 MeV，焦点尺寸2 mm的X射线直线加速器，并且通过人工智能算法自动标注X射线图像中的目标物［6-7］。以上这些装置基于吸收衬度和双能成像透视货物，大幅提高了工作效率［8-9］。但是生物样品主要由C、H、O、N等轻元素组成，不同组织之间密度差异很小［10-14］，通过上述高能量、大焦点X射线光源和基于吸收衬度的成像方法无法对弱吸收生物组织清晰成像，很难用于生物样品的品种和产地鉴别，现有的基于吸收衬度成像的安检设备也无法满足海关生物安全检测的需求。

近年来基于高亮度同步辐射光源发展起来的X射线相位衬度成像方法，利用空间相干X射线在样品内部传播时的相位变化成像，对于硬X射线来说，它的复折射率中代表相位的虚部比代表吸收的实部高两个量级以上，因而相位衬度成像的灵敏度远高于吸收衬度成像，特别适合于弱吸收样品成像。研发基于相位衬度成像的X射线生物特征识别仪，可以解决现有X射线成像安检仪无法实现弱吸收样品成像的难点问题。X射线相位衬度成像主要是基于高亮度、空间相干性好的同步辐射光源发展起来的［15］，其实现成像的方式主要包括干涉法［16］、衍射增强法［17］、光栅微分法［18］和同轴法［19］等实现方式，其中干涉法和衍射增强法对X射线光子的利用效率较低［20-22］。光栅微分法可降低对X射线光源空间相干性的要求，但需要复杂、精密的光栅及调节机构［18］，且光子利用效率低、数据采集时间长。同轴轮廓法基于准相干X射线在自由空间的传播实现相位衬度成像，具有光子利用效率高和光路简单的特点［19，23］，特别适用于高效率和高性价比X射线成像检测仪的研发。

本文报道了一种基于相位衬度成像的X射线生物特征识别仪。针对海关检测要求成像速度快的特点，采用光子利用效率高的X射线同轴轮廓相衬成像方法实现高效数据采集。海关检测样品种类多、尺寸和形状不统一，仪器结构设计中通过调节射线源、样品和探测器三者之间的距离，可选择合适的成像放大倍率，适应不同样品检测的需求。为保证检测环境的辐射安全、样品的切换及调整方便，采用防护棚屋屏蔽、门禁安全联锁和样品台自动控制的系统设计。此外，为解决二维投影像存在的信号堆叠和很难用于生物特征识别的问题，发展基于相衬成像的显微电子计算机断层扫描（Computed Tomography，CT）及其三维显微结构定量分析方法，以实现不同产地生物样品的识别。

1 成像原理及系统设计

1.1 X射线相衬成像原理

当X射线穿过物质时，除了会发生X射线的吸收外，还会伴随着X射线的折射现象，通常使用复折射率表示X射线与物质的相互作用［22，24-25］：

式中：δ、β≪1，当X射线在物质中传播时，其位相因子中呈现出一个衰减因子：

式中：n为X射线的复折射率；k=2πλ为波矢；λ为入射X射线的波长；z为沿光轴的传输距离。复折射率表达式中的δ项为相位项，β项为吸收项，相位项通常比吸收项高两个量级以上。这意味着物体对X射线相位的调制作用比吸收更加明显，利用X射线透过样品后携带的相位信息可以得到弱吸收样品精细结构的高灵敏成像［16，26-27］。

由于相位的记录是通过局域光干涉实现的，因此相位衬度成像对X射线光源的空间相干性要求较高［16］。理想的平行光或点光源是完全空间相干的，但实际上任何光源均有一定的相干空间尺度［28］。在非点源照明的情况下，物平面上的相干长度可表示为［19，29］：

式中：λ为光源波长；R1为光源到样品的距离；σ为光源点尺寸。光源尺寸越小，产生的X射线空间相干性越好。生物样品检测成像所需的光子能量通常为20 keV，对应的波长为λ=0.062 nm。当光源点尺寸σ=2μm、R1=0.1m时，相干长度d=3.1μm，相干长度d还会随着R1的增大而增大。

样品厚度和成分的变化导致X射线波前在通过样品后的形状发生变化。X射线透过样品后的相位变化φ可表示为［29］：

同轴相衬成像原理可由菲涅尔衍射来解释，根据菲涅耳衍射基尔霍夫公式的一阶近似，对于点光源照射纯相位样品的情况下，像平面的光强分布为［19，29］：

式中：z为近似光轴；k=2πλ为波矢；re为经典电子半径；ρ为样品的电子密度；R1为样品到光源的距离；R2为样品到探测器的距离；M=(R1+R2)R1为成像系统的几何放大倍率。像平面上的光强分布正比于电子密度的二阶微分，与光子能量无关［19］，这意味着可以实现多色光相衬成像。成像的空间分辨率很大程度上取决于光源的尺寸［29-30］，因此为同时保证相位衬度和高空间分辨，本仪器研发采用微焦点X光管作为成像系统的射线源。

归因于空间相干X射线的局域干涉，同轴相衬成像只能直接获得样品密度突变界面的边缘增强轮廓像。要获得样品的密度分布，还需发展相应的相位恢复算法。目前相位恢复算法主要有：基于Born近似和Rytov近似的算法［31］、强度传播方程（Transport of Intensity Equation，TIE）算法［32-33］和衬度传递函数（Contrast Transfer Function，CTF）算法［34］等。由于这些算法引入了一系列近似和假设，使得相位恢复结果不稳定。此外，其它一些相位恢复 算 法：如GSF（Gerchberg-Saxton-Fienup）算法［35-36］、Yang-Gu算 法［37］、Input-Output算 法［38］、Steepest-Descent算法［35］、Conjugate-Gradient算法［35］以及Wu等提出的迭代算法［39-40］，虽然稳定性更好，但运算量大，效率低。基于相位-吸收二元性（Phase-Attenuation Duality，PAD）的算法，因其单距、非迭代和稳定性好等优点，被广泛用于提取同轴相衬像的相位信息［30，41-42］。基于以上分析，采用PAD算法实现同轴相衬像的重构是X射线生物特征识别仪研发的最佳选择。

像平面的光强分布是经样品调制后的X射线波前在像平面干涉的结果，为了从单张投影图的光强分布中提取相位信息，假设样品为均质的，从平行光照明的强度输运方程（Transport-of-Intensity Equation，TIE）出发，结合菲涅尔衍射积分公式可以得到适用于提取X光管成像系统相位信息的PAD算法公式［30，43］：

式中：位置矢量r⊥垂直于光轴方向z；T(r⊥)为均匀物体在像平面上的投影厚度；Iin为入射光的平均强度；F{}和F-1{}分别代表傅里叶变换和傅里叶逆变换。

由式（5）可知：成像衬度随着R2的增大而增大［29］，考虑到该系统为基于锥形束的投影成像方式，成像的几何放大率M受R1和R2的共同影响。因此，在实际成像过程中应在保证所需放大倍率M的情况下尽可能增大R2以获得最佳衬度。

1.2 X射线成像系统

为提供同轴相衬成像所要求的准相干X射线，该系统采用微焦点透射式钨靶X光管（X-ray WorX GmbH，XWT-225-THE Plus），管 电 压 范 围20~225 kV，管电流范围20~1 000μA，焦点尺寸2μm，最大靶功率50 W，X射线张角160°，射线源可连续4 h提供稳定的X射线束。

为满足海关生物安全识别工作中样品尺寸差异大的问题，该系统采用基于GdOS闪烁体的互补金属氧化物半导体（Complementary Metal Oxide Semiconductor，CMOS）X射 线 平 板 探 测 器（Teledyne DALSA，Rad-icon 2329）实现大视场高分辨成像［44］，探测器单像素尺寸49.5μm，像素阵列4 608×5 890，成像面积228 mm×291mm。

图1为研制的X射线生物特征识别仪，其中图1（a）为成像光路示意图。该系统中X射线源、样品台和探测器同轴放置，其中X射线源固定，样品台可实现X、Y、Z三维平动和W一维旋转的四维运动。同时，样品台和探测器可沿X方向移动，以调节光源到样品的距离（Source-Object Distance，SOD）和光源到探测器的距离（Source-Detector Distance，SDD）来任意调节成像系统的几何放大倍率。SOD最小为10 mm，SDD最大为1 570 mm，此时为最大几何放大倍率157×，对应的有效像素尺寸0.32μm。根据成像条件的不同，曝光时间在数百毫秒到十几秒之间。图1（b）为仪器外观照片，辐射防护棚屋将整个X射线成像系统封闭起来，并配有相应的安全联锁系统，以保证测试过程中的辐射安全。配置了一台高性能计算机，可实现光源、样品台与及探测器的棚屋外控制及数据采集和处理。图1（c）为识别仪器的内部结构照片，仪器门可完全打开，以方便维护。前门配铅玻璃观察窗，以方便测试过程中实时观察样品及部件状态。该系统通过以太网（Ethernet）进行硬件控制与数据传输，整套软件运行于Windows操作系统，操作简单，便于海关部门对采集到的数据进行管理和对硬件的升级改造。此外，制定发展了完整的数据处理、设备维护和操作的使用规范，方便相关人员的使用。

图1 仪器的成像光路图(a)、仪器外观(b)及仪器内部照片(c)Fig.1 Diagram of the imaging optical path in the instrument(a),an appearance photo of the instrument(b),a photo inside the instrument(c)

由于二维投影成像存在信号堆叠问题，仅从二维投影图难以获取生物样品的特征结构信息。针对生物样品三维空间结构信息获取和分析的需求，将CT成像技术引入到海关生物特征识别中，通过对生物样品相位衬度三维信息的定量分析，实现同种、不同产地生物样品的识别。为此，在DALSA GigE软件开发包（Software Development Kit，SDK）的基础上，结合KOHZU电机软件控制接口，开发出具有自定义曝光时间、单张拍摄、单张保存、实时预览、电机控制以及CT扫描功能的软件（图2），使该仪器具备拍摄CT的能力。通过该软件对样品进行CT扫描，可快速无损获取被检测样品的微米级三维空间结构，有效提高生物安全检测的准确率。

图2 仪器CT扫描程序流程图Fig.2 Flow chart of the CT scanning program of the instrument

2 测试结果与分析

2.1 投影成像测试

首先对该系统进行了二维投影成像测试，采用分辨率靶（JIMART RC-02）测试该系统的空间分辨率，采用小鱼和香樟叶测试生物样品成像效果，实验条件见表1。成像时有效像素尺寸的选择取决于成像所需的分辨率和视场大小：为了测试该系统的空间分辨率，采用0.50μm的有效像素对分辨率靶成像；考虑到样品较大，采用8.15μm的有效像素对小鱼成像，4.95μm的有效像素对香樟叶成像。曝光时间的选择取决于成像信噪比，成像信噪比与到达像素点的光子数成正比，而到达像素点的光子数与光源靶功率、有效像素尺寸的平方、探测器闪烁体对不同能量X射线的响应系数成正比，与SDD的平方成反比，且与样品材质和厚度相关，因此实际的曝光时间要从以上几个方面来综合考虑。分辨率靶成像测试采用的靶功率小、SDD大、有效像素尺寸小，需要5 s的曝光时间才能得到高信噪比图像；小鱼成像测试采用的靶功率大、SDD小、有效像素尺寸大，小鱼对X射线的调制能力强，0.7 s的曝光时间即可得到高信噪比图像；香樟叶成像测试采用的管电压小、靶功率小，探测器闪烁体对低能X射线的响应系数小、样品薄，需要10 s的曝光时间才能得到高信噪比图像。测试结果如图3所示。综合图3（a）的投影像及图3（b）的强度分布轮廓图，可以认为该成像系统的空间分辨率为3μm。图3（c）小鱼二维投影像说明该系统可清晰分辨小鱼体内不同层次的骨骼、鱼肠内食物和鱼鳍，考虑到样本已脱水，这应该是相位衬度和吸收衬度共同作用的混合衬度成像。根据得到的样本内部结构信息，可以用来推断小鱼样本的种类及相关信息，从而用于样本识别。图3（d）所示香樟叶二维投影像表明该系统可清晰分辨叶片上的各级叶脉，这属于相位衬度占主导地位的成像，从而验证了所研制的X射线生物特征识别仪具备较好的相位衬度成像能力。上述投影成像测试结果表明，该系统可用于简单结构生物样品或薄样本的生物特征识别。

表1 二维成像测试条件Table 1 Test conditions of the two dimensional imaging

图3 分辨率靶投影像(a)，分辨率靶投影像标记处的强度分布轮廓图(b)，小鱼投影像(c)，香樟叶投影像(d)Fig.3 Projection image of the resolution target(a),intensity profile at the mark in Fig.3a(b),projection image of the tiddler(c),projection image of the cinnamomum camphora leaves(d)

2.2 CT成像测试

相衬CT功能主要用于弱吸收样品的三维显微结构成像，适用于厚样品原位无损测试及定量分析识别。为测试其相位衬度成像效果，主要采用不含水样品测试，以降低吸收衬度的影响。采用小鱼样品和药物微丸进行三维成像测试，脱水小鱼样品属相衬与吸收的混合衬度成像，药物微丸由于均为低Z元素组成且尺度小，可看成纯相位样品。测试过程中，通过CT扫描软件采集不同角度下的投影图像，使用PAD算法进行相位恢复，采用基于极坐标变换的方法去除切片的环状伪影［45］，使用Avizo 9.7.0（Thermo Fisher）软件对切片进行三维渲染。小鱼测试结果见图4，其测试条件为管电压50 kV、管电流700μA、靶功率21.3 W、曝光时间700 ms、有效像素尺寸8.15μm。图4（a、b）是小鱼吸收衬度成像结果，图4（a）为小鱼骨骼及肠道内食物残渣的三维结构，图4（b）给出了小鱼腹部方向的视图，可清晰观察到尾椎骨的整体显微结构，且可明显看出侧部鱼骨断裂和错位的情况，应该是样品干燥后挤压所致。图4（c、d）是小鱼样品的相位衬度成像结果，图4（c）为小鱼体表的相位衬度成像三维结构，可观察到小鱼体表鱼鳞的排列结构，图4（d）给出小鱼腹部方向的剖视图，可清晰观察到鱼肠的形状、空间位置以及鱼肠内食物残渣的具体分布情况。图4表明该系统具有较高的密度分辨能力，适用于对生物样品进行高衬度成像。图5是纯相位样品药物微丸的三维显微结构，测试条件为管电压50 kV、管电流500μA、靶功率16 W、曝光时间3 s、有效像素尺寸3 μm。样品准备过程中将含微丸的药片压碎，挑选出药物微丸，为CT成像方便，将微丸粘在塑料吸管上。图5（a）是吸管上药物微丸的三维空间分布以及表面形貌结构，图5（b）为微丸及吸管的横截面三维视图，可清晰观察到药物微丸及吸管内部存在多层复杂空腔结构。图5（c）为单个微丸的三维显微结构，可见微丸从内到外依次由空腔、药物层、隔离层、肠衣层组成，各层之间存在一定的间隙。由于微丸及吸管样品在管电压50 kV的X射线照射下可看做纯相位物体，从而验证了所研制的仪器具有较高的相位衬度成像分辨能力。

图4 小鱼骨骼及肠道内食物残渣的三维结构(a)，腹部方向视图(b)，体表结构(c)，纵向剖视图(d)Fig.4 Three-dimensional structure of the skeleton and food residues of the tiddler(a),orientation view of the abdomen(b),surface structure(c),longitudinal cutaway view(d)

图5 药物微丸粘附于塑料吸管的三维空间分布及表面形貌结构(a)，微丸及吸管截面三维视图(b)，单个微丸三维显微结构(c)Fig.5 Three-dimensional spatial distribution and surface morphology of pellets adhered to a pipette(a),three-dimensional crosssection view of pellets and pipette(b),three-dimensional microstructure of a single pellet(c)

2.3 生物样本识别测试

对于较难根据外观鉴别的生物样本，采用显微CT实现其内部结构的无损成像，并进一步定量分析，可以实现同种、不同产地生物样本的识别。选取不同产地的核桃作为测试样本，对X射线而言这是一个混合衬度样品，核桃衣可以作为纯相位物体。由于核桃被坚硬的外壳包裹，常规的检测方法很难区分不同产地的核桃，通过研制的X射线生物特征识别仪对核桃进行CT扫描，再对切片进行定量分析可达到无损鉴别核桃产地的目的。测试样本分别为产自云南省楚雄彝族自治州（Chuxiong）、浙江省杭州市（Hangzhou）和河南省洛阳市（Luoyang）三地的核桃，每个产地各取3个样本，测试条件均为管电压80 kV、管电流100μA、靶功率6.6 W、曝光时间500 ms、有效像素尺寸19.09μm。三维成像结果见图6，图6（a、b）为核桃的三维体积渲染图，从图6（b）核桃的横截面视图可以看出，识别仪可以实现核桃壳、核桃衣和核桃肉的清晰成像。图6（c）为核桃衣相衬CT成像结果，其中浅色部分为核桃衣，深色实体部分为核桃壳，说明识别仪可以对核桃衣这类纯相位物体进行成像。图6（d）为核桃仁的三维渲染图，表明该识别仪可无损获取样本内部复杂的三维显微结构，通过对切片进行阈值分割实现核桃壳和核桃肉的结构信息提取，并对核桃样品进行进一步的图像定量分析。

图6 核桃的三维成像体积渲染图，核桃整体(a)，核桃横截面(b)，核桃衣(c)，核桃仁(d)Fig.6 3D volume renderings of the walnut,the whole body(a),a cross-section view(b),walnut coat(c),walnut kernel(d)

通常仅从形态差异较难实现不同产地样本的识别，依据三维图像的定量分析可以获取更多的样品特征信息。对核桃样品特征作了定量分析（图7），结果表明：云南省楚雄彝族自治州所产核桃的平均外壳厚度为914.5μm，果仁平均体积占比为36.0%；浙江省杭州市所产核桃的平均外壳厚度为1 114.9μm，果仁平均体积占比为29.7%；河南省洛阳市所产核桃的平均外壳厚度为1 543.2μm，果仁平均体积占比为33.9%。结合各产地的经纬度：楚雄（101°30′E，25°02′N），杭州（120o10′E，30o15′N），洛阳（112°26′E，34°43′N），从图7（b）可知，核桃壳厚度与产地的纬度存在明显的正相关性。通过对CT扫描产生的三维数据进行定量分析可以看出，不同产地的核桃外壳平均厚度和果仁体积占比有明显差异，因此能以此来无损鉴定核桃的产地。尽管实际物种鉴别中还需要更多的检测样品，但初步的测试结果表明，研制的生物特征识别仪可为生物样本鉴别提供更多有用的信息。基于该仪器建立的检测方法，可推广到其它生物物种的识别或显微鉴定中，从而提高海关生物安全检测的准确率。

图7 核桃样品特征的定量分析结果 (a)果仁体积占比，(b)外壳平均厚度Fig.7 Quantitative analysis results of characteristics of walnut samples(a)Proportion of kernel volume,(b)Average shell thickness

3 结语

根据海关生物安全检测对于数据采集效率高、成像速度快的需求，成功研发了一台基于同轴轮廓相衬成像的X射线生物特征识别仪。为了解决海关检测样品种类复杂、尺寸和形状差异大的问题，将微焦点X光管与CMOS大面阵X射线平板探测器相结合，利用锥形X射线束投影成像的几何放大作用，通过调整样品、探测器与光源的相对距离，实现成像系统放大倍率的调节以适应实际应用场景。此外，该方法也降低了高分辨成像对探测器像素尺寸的要求。结合相衬显微CT成像技术，该识别仪可用于定量分析生物样品的三维结构特征，实现不同产地生物样品的识别，从而提升海关生物安全检测工作的准确率。分辨率靶测试结果表明：采用单像素尺寸49.5μm的大面积平板探测器，实现了3μm的空间分辨率。生物样品测试结果表明：相比目前海关通用的基于吸收衬度成像的X射线检测手段，新研发的基于相位衬度成像的X射线生物特征识别仪具有更高的空间分辨和密度分辨能力，更适合用于弱吸收样品的成像。此外，研发的识别仪通过以太网进行硬件控制、数据采集与处理，便于海关部门对采集到的数据进行管理和对硬件的升级改造。该系统还可用于中药材显微鉴定、结构药剂学、聚合物材料等领域的研究，为相关学科的发展提供一个有效的相衬成像显微CT检测平台。

致谢 感谢中国科学院上海药物研究所张继稳研究员课题组提供的药物微丸样品以及中国科学院上海应用物理研究所骆钧尧在CT扫描软件编写过程中提供的帮助。