中华放射医学与防护杂志  2018, Vol. 38 Issue (7): 541-546   PDF    
基于蒙特卡罗和CT数据体素模型的双能X射线吸收法腰椎骨密度成像参数优化
李师 , 皮一飞 , 霍万里 , 汪志 , 陈志 , 徐榭     
230027 合肥, 中国科学技术大学核科学技术学院
[摘要] 目的 研究双能X射线吸收法(DXA)腰椎成像过程并优化成像参数,为提高双能X射线骨密度仪的成像质量提供理论依据。方法 采用蒙特卡罗方法,构建基于计算机断层扫描数据(CT)的体素模型,利用模拟产生的高低能X射线能谱,分别照射体素模型腰椎部位得到骨密度图像,计算每幅图像的品质因数(FOM)。根据FOM大小优化高低能管电压组合、过滤铜片厚度以及入射粒子数比值。结果 FOM在管电压组合为75与200 kVp达到最低值1.59×10-2。当过滤铜片厚度由0 mm增加至3 mm时,FOM数值从6.30×10-2下降到1.87×10-2,且下降趋势逐渐变缓。随着低能与高能的入射粒子数比从1:3增加到19:1,FOM先降低后升高,在3:1时达到最低值1.40×10-2结论 低能管电压为70~85 kVp,高能管电压为160~200 kVp,产生高能能谱时额外添加过滤铜片厚度0.3 mm,入射粒子数比值(低能/高能)为1~5时,能在获得较好的DXA腰椎图像质量的同时使患者受到较小的辐射剂量。
[关键词] 蒙特卡罗模拟     计算机体素模型     双能X射线     骨密度成像     成像参数优化    
Optimization of DXA lumber spine bone imaging parameters based on Monte Carlo method and CT voxel phantom
Li Shi, Pi Yifei, Huo Wanli, Wang Zhi, Chen Zhi, Xu Xie     
School of Nuclear Science and Technology, University of Science and Technology of China, Hefei 230027, China
Fund programs: Key Research Development Program of China (2017YFC0107500); National Natural Science Foundation of China (113751182, 113751181)
Corresponding author: Chen Zhi, Email:zchen@ustc.edu.cn
[Abstract] Objective To study the lumber spine imaging process of dual-energy X-ray absorptiometry (DXA) and parameters used to optimize the image quality. Methods A computational voxel phantom was constructed from patient computed tomography (CT) data. Using the Monte Carlo radiation transport method, a dual energy x-ray beam was simulated to scan the phantom of lumbar spine to generate a bone density image. The Figure of Merit (FOM) of each image was claculated. Parameters including the combination of the high and low energy tube voltage, the thickness of Cu filter, and the ratio of two beam energy incident photon number were optimized, which based on FOM. Results FOM reaches a minimum of 1.59×10-2 with the tube voltage combination of 75 and 200 kVp. With the thickness of the Cu filter from 0 mm to 3 mm, FOM decreases from 6.30×10-2 to 1.87×10-2, showing a gradually slow-down trend. With the incident photon number ratio (low energy/high energy) increasing from 1:3 to 19:1, FOM decreases firstly and then increases, reaching a minimum of 1.40×10-2 at 3:1. Conclusions According to the simulation results, the combinations of low tube voltage from 70 kVp to 85 kVp and high tube voltage from 160 kVp to 200 kVp, 0.3 mm Cu filter and beam incident photon number ratio from 1 to 5 can yield the best lumbar spine image quality with the lowest patient dose.
[Key words] Monte Carlo     Computational voxel phantom     Dual-energy X-ray     Bone density imaging     Imaging parameters optimization    

双能X射线骨密度仪是目前最常用的骨质疏松诊断仪器[1-2],其原理是采用双能X射线吸收法(dual-energy X-ray absorptiometry,DXA)即高能和低能两种X射线穿过人体组织[3],根据不同能量X射线在不同组织中的衰减系数不同的特点,计算出骨骼的面密度(即单位投影面积的物体质量)并得到骨密度图像。近年来,随着基于DXA骨密度图像评估骨质疏松和骨折风险的新方法不断发展,对DXA成像质量的要求越来越高[4-6]。然而制约DXA图像质量的一个重要因素是非最优化的成像参数。本研究构建了基于CT数据的体素模型,利用蒙特卡罗通用软件MCNPX[7]模拟DXA腰椎成像过程,通过改变不同双能X射线管电压组合,添加不同厚度的过滤铜片以及不同入射粒子数比值,在尽可能低的剂量范围内,以期获得DXA腰椎成像的最优化参数。

材料与方法

1.蒙特卡罗方法模拟成像过程

(1) 体素模型的构建:选取1名77岁男性的CT数据,编写MATLAB脚本,读取CT的DICOM文件,获取HU值。采用Schneider方法对HU值进行划分[8];并根据HU值与物质密度的转换关系和HU值与材料组成的对应关系,计算出每个像素的密度与材料组成,将其分别归类到24种组织材料,每个材料都有其独有的元素组成和密度;按照像素在CT上的空间位置拼接每个体素,所有体素填充构成可用于模拟计算的体素模型[9]

(2) 蒙特卡罗模拟扫描过程:模拟骨密度仪电压切换法[10],由SpekCalc软件计算获得钨阳极X射线管能谱[11-12]。光子从平面面源出射,沿着正位后前方向照射体素模型的腰椎区域[13],射野大小为10 cm×18 cm。照射区域的人体厚度约为22 cm,处于正常腰椎厚度范围内,可以忽略能谱硬化现象的影响[14]。探测器分为200×360个理想探测单元,每个探测单元尺寸为0.5 mm×0.5 mm×15 cm。使用计数卡F4记录每个探测单元光子通量(光子强度,photons/cm2)和相对误差,使用计数卡F6记录整个体模吸收剂量,单位MeV/(g·photon)换算为μGy/photon。

2. DXA骨密度图像的计算:DXA骨密度计算公式[2, 15]为:

$ {m_{\rm{b}}} = \frac{{\ln \frac{{I_0^L}}{{{I^L}}} - {R_{\rm{s}}} \cdot \ln \frac{{I_0^H}}{{{I^H}}}}}{{\tau _{\rm{b}}^l - {R_{\rm{s}}} \cdot \tau _{\rm{b}}^h}} $ (1)

式中,mb为骨骼的面密度,g/cm2ILIH分别为低能和高能的透射线强度,I0LI0H分别为低能和高能的入射线强度;τblτbh分别为低能和高能X射线在骨骼中的质量衰减系数,τslτsh分别为低能和高能X射线在软组织中的质量衰减系数,cm2/g;${R_{\rm{s}}} = \frac{{\tau _{\rm{s}}^l}}{{\tau _{\rm{s}}^h}}$,为脂肪与肌肉含量的测量值,须提前被定义以计算骨密度图像中每个像素点数值,通过扫描不含骨骼的软组织获得Rs[15]。若mb=0,Rs为:

$ {R_{\rm{s}}} = \frac{{ - \ln \frac{{I_{\rm{s}}^L}}{{I_0^L}}}}{{ - \ln \frac{{I_{\rm{s}}^H}}{{I_0^H}}}} $ (2)

式中,IsLIsH分别为低能和高能X射线照射去除骨骼只存在软组织体模的透射线强度。双能X射线能谱分别照射完整体模、空气、去除骨骼的软组织体模,即双能能谱组合得到6幅投影图像,对应IiI0iIsii=L, H。将这6幅图像的每个像素数据对应代入公式(2)和(1)中,其中不同能量X射线在骨骼中的质量衰减系数从NIST数据库中查表获得[16],在临床成像情况下,通过能谱的平均能量插值得到能谱质量衰减系数等效值,从而计算得到被照射区域的DXA骨密度图像,图中每个像素的数据值代表了该点的骨密度值。

3.骨密度图像噪声、吸收剂量与品质因数的计算

(1) 骨密度图像噪声:噪声是衡量图像质量的重要指标,在此只考虑由探测器记录的X射线强度统计特性引起的量子噪声,该噪声符合高斯分布,随着曝光量的增加而减小。对公式(1)进行误差分析[17-19],并根据MCNPX中相对误差R定义—计数结果平均值标准差与计数结果平均值之比,得到由ILIHI0LI0H统计误差导致的骨密度图像噪声,用骨密度值标准差σmb表示:

$ \sigma _{{m_{\rm{b}}}}^2 = {\left( {\frac{1}{{\tau _{\rm{b}}^l - \tau _{\rm{b}}^h{R_{\rm{s}}}}}} \right)^2}\left( {R_{I_0^L}^2 + R_{{I^L}}^2} \right) + {\left( {\frac{{ - {R_{\rm{s}}}}}{{\tau _{\rm{b}}^l - \tau _{\rm{b}}^h{R_{\rm{s}}}}}} \right)^2}\left( {R_{I_0^H}^2 + R_{{I^H}}^2} \right) $ (3)

为了方便比较分析,只需计算每幅DXA骨密度图像的噪声σ2mb平均值。获得Rs均值后,由模拟结果得到相对误差R2I0LRIL2R2I0HRIH2平均值。因为能谱的有效能量相当于具有相同衰减效果的某一单能光子射束的辐射能量[20],在计算噪声的平均值时,通过能谱的有效能量插值得到质量衰减系数等效值。

(2) 吸收剂量:用整个体模吸收剂量D表示患者所受的辐射剂量,因为MCNPX所得结果为相对剂量,为了得到总吸收剂量D(μGy),计算如下:

$ D = {D_1}{N_1} + {D_{\rm{h}}}{N_{\rm{h}}} $ (4)

式中,DlDh分别为低能和高能照射时的吸收剂量,μGy/photon;NlNh分别为低能和高能照射的入射粒子数。

(3) 品质因数FOM:最优化的DXA成像参数能够在保证目标噪声的前提下,使得患者所受的辐射剂量尽可能小,品质因数(FOM)是同时评价图像质量及辐射剂量的综合指标,但在不同文献中定义有所不同,本研究根据Bolin和Preuss提出的方法[21-22],将品质因数FOM表示为:

$ F = \sigma _{{m_{\rm{b}}}}^2D $ (5)

吸收剂量D与入射粒子数成正比例关系,当入射粒子数足够多时,图像噪声σ2mb与入射粒子数成反比例关系,可知FOM与入射粒子数无关。

本研究针对常用的管电压取值区间40~200 kVp进行优化,并模拟大范围的入射粒子数比以及添加常见过滤铜片厚度观察其对品质因数FOM的影响规律。

结果

1.骨密度平均值计算:为了验证模拟是否正确,在获得骨密度图像后,计算了骨密度平均值,并与临床值比较。利用原始CT的HU值>300的像素投影,将其灰度叠加投影在一个平面上,划分骨骼和软组织区域,定位相关区域坐标,将骨密度的图像叠加在相应位置,如图 1所示。不同成像参数下骨密度计算平均值大约在0.6~1.5 g/cm2范围内,与临床值范围基本符合[23],并且由于图像噪声基于骨密度值标准差,表明图像噪声越低的同时骨密度值也越准确,管电压组合为80和200 kVp时,图像噪声较低,骨密度计算平均值为0.991 g/cm2,相比文献中70~79岁中国正常男性0.964 g/cm2[23],相对误差仅为0.028,验证了所建模型和模拟方法的正确性。

图 1 骨密度图像叠加在CT的灰度投影相应位置示意图 Figure 1 A bone density image superimposed on the corresponding position of the gray scale projection of CT

2.双能X射线管电压组合优化:低能管电压40~100 kVp,间隔5 kVp,共13组低能管电压数据,过滤片2.7 mm Al;高能管电压110~200 kVp,间隔10 kVp,共10组高能管电压数据,过滤片2.7 mm Al+ 0.3 mm Cu;入射粒子数比1 :1。不同高低能管电压组合得到130幅骨密度图像,计算每幅图的FOM。将FOM与不同管电压组合作等高线图,如图 2所示,低能取值在70~80 kVp,高能取值在160~200 kVp时,FOM值较小。由图 3图 4可知,低能管电压为75 kVp,高能管电压为200 kVp时,FOM达到最低值1.59×10-2

图 2 管电压组合与品质因数关系图 Figure 2 Correlation of tube voltage combination with FOM

图 3 高能管电压固定为200 kVp与不同能量的低能管电压组合时,图像噪声、吸收剂量、品质因数变化规律图  A.图像噪声;B.吸收剂量;C.品质因数FOM Figure 3 Correlation of tube voltage combination with image noise, absorbed dose and FOM at combinations of 200 kVp high-energy tube voltage with different low-energy tube voltages  A. Image noise; B. Absorbed dose; C. FOM

图 4 低能管电压固定为75 kVp与不同能量的高能管电压组合时,图像噪声、吸收剂量、品质因数变化规律图   A.图像噪声;B.吸收剂量;C.品质因数FOM Figure 4 Correlation of tube voltage combination with image noise, absorbed dose and FOM at combinations of 75 kVp low-energy tube voltage with different high-energy tube voltages  A. Image noise; B. Absorbed dose; C. FOM

3.过滤铜片厚度优化:为了减少高低能谱之间的重叠,在产生高能能谱时额外添加一定厚度的过滤铜片,将高能能谱预硬化,过滤软X射线部分[15]。产生高能能谱时,本研究在固定2.7 mm过滤铝片基础上添加的不同厚度的过滤铜片,管电压组合为55和120 kVp,入射粒子数比1 :1。由表 1可知,随着过滤铜片厚度增加,吸收剂量升高,图像噪声和FOM下降。

表 1 不同厚度的过滤铜片对应的吸收剂量、图像噪声和品质因数 Table 1 Absorbed dose, image noise, and FOM for different thicknesses of the Cu filter

4.入射粒子数比优化:添加过滤铜片0.3 mm,管电压组合为70与200 kVp,保持总入射粒子数始终为1×109,定义入射粒子数比r=Nl/Nhr值分别为1 :3、1 :1、3 :2、3 :1、4 :1、9 :1、19 :1。如表 2所示,随着比值增加,吸收剂量降低,图像噪声和FOM先降低后升高,当r为3 :2时,图像噪声最小,当r为3 :1,FOM值最小。

表 2 不同入射粒子数比对应的吸收剂量、图像噪声和品质因数 Table 2 Absorbed dose, image noise, and FOM in relation to different incident photon number ratios

讨论

当管电压取值在常用诊断能量区间40~200 kVp,正位腰椎成像情况下,模拟结果显示:低能管电压为70~85 kVp,高能管电压为160~200 kVp;产生高能能谱时额外添加0.3 mm过滤铜片;入射粒子数比r值为1~5,能在获得较好的DXA图像质量的同时使患者受到较小的辐射剂量。不同成像参数下图像质量差别较大,与模拟的最差参数相比,理想参数的选择使FOM值降低约14倍。

管电压组合优化结果显示,吸收剂量随着管电压取值增大而升高,当高能管电压固定时,低能管电压为80 kVp时图像噪声最小,75 kVp时FOM值最小;当低能管电压固定时,范围在200 kVp以下时,高能管电压取值越高,图像噪声和FOM值越低。双能X射线骨密度仪生成厂家美国Hologic公司,通过电压切换法产生70 kVp与140 kVp的能谱,美国Lunar公司、Norland公司以及Sopha公司通过K边缘过滤法产生的能谱为40 keV与70 keV[13],这说明目前大部分商用的骨密度仪能量选择不能达到最理想的成像效果。

为了分离双能X射线,在原有过滤铝片基础上添加过滤铜片,本研究结果显示,仅添加0.1 mm铜片,图像噪声和FOM明显下降;当过滤铜片选择0.3 mm时图像质量已明显提高,继续增加过滤铜片厚度,FOM值下降缓慢但是患者受到的辐射剂量却迅速增加。结果表明过滤铜片对图像质量有显著影响,但过度增加铜片厚度图像质量并不会有太大的提升反而会使患者受到更多不必要的辐射剂量。

相比高能X射线,低能X射线衰减更加严重,到达探测器的透射线更少,所以统计误差更大,入射粒子数比值r < 1,意味着更少的低能X射线,这种统计误差将大幅度增加,因此r值越小FOM值越大;r>1,意味着更多的低能X射线,这种统计误差将减小,但当r值过大时,高能的统计误差增大,FOM值升高。所以,r取值不应 < 1且不应过大,结果显示r值为1~5时FOM较小。本研究使用蒙特卡罗方法模拟骨密度成像,并计算骨密度值,将其与临床值相比较,验证了模拟的正确性。由模拟结果直接得到患者全身吸收剂量,根据MCNPX模拟结果误差定义计算得到图像噪声,从而获得FOM,优化了成像参数。相比以往研究中使用X射线平均能量来代替X射线能谱以及基于简单的几何模型,本文模拟了双能X射线能谱照射并基于真实人体CT数据构建的体素模型,所得结果更符合实际情况。

由于FOM还与组织的厚度和成分等因素有关[23],以上结果仅针对标准人体厚度的正位腰椎扫描,对于其他体型患者(肥胖患者、小孩等)以及其他人体部位(侧位腰椎、前臂骨等)在骨密度扫描时的参数选择会有所变化。另外本研究为理论模拟研究,只考虑了统计特性引起的噪声,未考虑实际噪声情况,但所得结果可为临床实际应用提供理论依据。

志谢 林卉、裴曦提供蒙特卡罗模拟的指导和论文修改意见,在此表示感谢
利益冲突 全体作者未因进行该研究而接收任何不正当的利益,在此对研究的独立性和科学性予以保证
作者贡献声明 李师设计研究方案、处理分析数据并撰写论文;皮一飞构建体素模型、编写骨密度值计算程序;霍万里、陈志、徐榭负责指导研究思路;汪志提供CT数据
参考文献
[1]
Kanis JA, McCloskey EV, Johansson H, et al. A reference standard for the description of osteoporosis[J]. Bone, 2008, 42(3): 467-475. DOI:10.1016/j.bone.2007.11.001
[2]
Blake GM, Fogelman I. Technical principles of dual energy x-ray absorptiometry[J]. Semin Nucl Med, 1997, 27(3): 210-228. DOI:10.1016/S0001-2998(97)80025-6
[3]
梁昊, 马成兴, 朱俊杰. DEXA骨密度仪中的图像处理与定量计算[J]. 核电子学与探测技术, 2002, 22(3): 263-264.
Liang H, Ma CX, Zhu JJ. Image processing and quantitative computation of DEXA bone density[J]. Nucl Electron Detection Technol, 2002, 22(3): 263-264. DOI:10.3969/j.issn.0258-0934.2002.03.021.issn.0258-0934.2002.03.021
[4]
Silva BC, Leslie WD, Resch H, et al. Trabecular bone score:a noninvasive analytical method based upon the DXA image[J]. J Bone Miner Res, 2014, 29(3): 518-530. DOI:10.1002/jbmr.2176
[5]
Jager PL, Jonkman S, Koolhaas W, et al. Combined vertebral fracture assessment and bone mineral density measurement:a new standard in the diagnosis of osteoporosis in academic populations[J]. Osteoporos Int, 2011, 22(4): 1059-1068. DOI:10.1007/s00198-010-1293-3
[6]
Winzenrieth R, Michelet F, Hans D. Three-dimensional (3D) microarchitecture correlations with 2D projection image gray-level variations assessed by trabecular bone score using high-resolution computed tomographic acquisitions:effects of resolution and noise[J]. J Clin Densitom, 2013, 16(3): 287-296. DOI:10.1016/j.jocd.2012.05.001
[7]
Pelowitz DB. MCNPX user's manual version 2. 7. 0[M]. Los Alamos: Los Alamos National Laboratory, 2011.
[8]
Schneider W, Bortfeld T, Schlegel W. Correlation between CT numbers and tissue parameters needed for Monte Carlo simulations of clinical dose distributions[J]. Phys Med Biol, 2000, 45(2): 459-478. DOI:10.1088/0031-9155/45/2/314
[9]
冯铓, 张练, 陈志, 等. DICOM2 Phantom: 从CT的DICOM图像到蒙卡程序MCNP的转换软件的初步开发[C]. 合肥: 第一届合肥国际放射医学物理论坛, 2015.
Feng M, Zhang L, Chen Z, et al. DICOM2 Phantom: initial development of conversion software from CT's DICOM image to Monte Carlo program MCNP[C]. Hefei: The First Hefei International Forum for Radiological Medical Physics (HIFRMP), 2015.
[10]
李婵, 李亮, 陈志强. 双能X射线骨密度仪技术进展综述[J]. CT理论与应用研究, 2014, 23(5): 717-730.
Li C, Li L, Chen ZQ. Review of developments of the dual energy X-ray absorptiometry techniques[J]. CT Theory Appl, 2014, 23(5): 717-730.
[11]
Poludniowski GG. Calculation of x-ray spectra emerging from an x-ray tube. Part Ⅱ. X-ray production and filtration in x-ray targets[J]. Med Phys, 2007, 34(6): 2175-2186. DOI:10.1118/1.2734726
[12]
Poludniowski G, Landry G, DeBlois F, et al. SpekCalc:a program to calculate photon spectra from tungsten anode x-ray tubes[J]. Phys Med Biol, 2009, 54(19). DOI:10.1088/0031-9155/54/19/N01
[13]
Adams J. Single-and dual-energy:X-ray absorptiometry[J]. Eur Radiol, 1997, 7(0938-7994): S20-S31. DOI:10.1007/pl00006861
[14]
Blake GM, McKeeney DB, Chhaya SC, et al. Dual energy x-ray absorptiometry:the effects of beam hardening on bone density measurements[J]. Med Phys, 1992, 19(2): 459-465. DOI:10.1118/1.596834
[15]
International Atomic Energy Agency. Dual energy X ray absorptiometry for bone mineral density and body composition assessment[M]. Vienna: IAEA, 2011.
[16]
Petoukhova AL, van Egmond J, Eenink MG, et al. The ArcCHECK diode array for dosimetric verification of HybridArc[J]. Phys Med Biol, 2011, 56(16): 5411-5428. DOI:10.1088/0031-9155/56/16/021
[17]
Shaw CC, Gur D. Comparison of three different techniques for dual-energy subtraction imaging in digital radiography:a signal-to-noise analysis[J]. J Digit Imaging, 1992, 5(4): 262-270. DOI:10.1007/bf03167808
[18]
Michael GJ, Henderson CJ. Monte Carlo modelling of an extended DXA technique[J]. Phys Med Biol, 1998, 43(9): 2583-2596. DOI:10.1088/0031-9155/43/9/011
[19]
Lemacks MR, Kappadath SC, Shaw CC, et al. A dual-energy subtraction technique for microcalcification imaging in digital mammography--a signal-to-noise analysis[J]. Med Phys, 2002, 29(8): 1739-1751. DOI:10.1118/1.1494832
[20]
涂振邦, 陳為立, 朱健豪, 等. 低能量X射線射質之研究[J]. 中華放射線醫學雜誌, 2004, 29(4): 185-190.
Tu ZB, Chen WL, Zhu JH, et al. The study of the beam quality of the low energy X-rays[J]. Chin J Radiol, 2004, 29(4): 185-190.
[21]
Bolin FP, Preuss LE. Optimum photon energies for the measurement of bone mineral and fat fractions[J]. Br J Radiol, 1978, 51(609): 750. DOI:10.1259/0007-1285-51-609-750
[22]
Michael GJ, Henderson CJ. Simulation studies of optimum energies for DXA:dependence on tissue type, patient size and dose model[J]. Australas Phys Eng Sci Med, 1999, 22(4): 126-35.
[23]
梅铁成, 高玉红, 梅菊红, 等. 骨折病人与正常人骨密度的对比分析[J]. 中国骨质疏松杂志, 2005, 11(2): 210-211.
Mei TC, Gao YH, Mei JH, et al. Analysis of bone mineral density in fractured patients and healthy subject[J]. Chin J Osteopor, 2005, 11(2): 210-211. DOI:10.3969/j.issn.1006-7108.2005.02.020.issn.1006-7108.2005.02.020