2. 中国科学院大学 北京 100049;
3. 上海艾普强粒子设备有限公司 上海 201800
2. University of Chinese Academy of Sciences, Beijing 100049, China;
3. Shanghai APTRON Particle Equipment Corporation, Shanghai 201800, China
与传统的光子放疗相比，质子放疗由于其Bragg峰的存在，可以为治疗提供更好的剂量分布，从而提高肿瘤控制率，降低副作用[1]。目前质子放疗技术的照射方式主要有散射和笔型束扫描两种方式。相比于前者，后者的束流利用率更高，产生更低的中子剂量，并且无需患者定制准直器和补偿器，提供更好的适形照野。笔型束扫描被认为是质子放疗的发展方向，在建的质子放疗装置基本都使用笔型束扫描方式。笔型束治疗对加速器及照射系统提出了极高的要求。剂量分布能否达到临床要求，依赖于点扫描中每个点的位置精度和剂量精度。因此笔型束扫描治疗系统的调试和研究，需要测量质子束的位置和流强。
上海市战略性新兴产业重大项目：首台国产质子治疗示范装置由中国科学院上海应用物理研究所研发，装置采用同步加速器，其中固定束治疗室采用笔束扫描主动方式照射。目前对笔束照射位置的精确测量主要使用胶片、位置电离室、荧光靶等方法。Lynx-2D是一款由IBA公司开发的商用荧光靶探测器，它可以提供在线的测量，及时分析每次的测量结果。尤其在治疗头的调试中，极大提高工作效率[2]。
本文主要介绍了IBA公司商用Lynx-2D在国内首台质子放疗示范装置点扫描系统中的使用情况。由于Lynx-2D自身软件无法获取束斑照射点的位置，我们基于MATLAB开发了相应的分析工具，分别采用高斯分布拟合法和灰度重心法定义照射图像束斑的位置点，并测试了两种方法的有效性。
1 质子笔型束扫描照射质子笔束扫描的照射方式如图 1所示。首先将照射肿瘤靶区沿深度方向分为多层，每一层对应一个束流能量，通过控制加速器的引出能量，使束流到达不同层。在同一能量层内，由不同的照射点组成。通过扫描磁铁控制束流的横向偏转，实现对预定位置的照射，一般从所需最高能量开始照射，直到所有能量层上的照射点都照射完毕。中间还有剂量电离室监控照射剂量，条带电离室监控束流的照射位置。
|
图 1 质子笔束扫描示意图 Figure 1 Schematic of proton pencil beam scanning |
图 2为Lynx-2D (IBA)实物图。Lynx-2D主要由闪烁体屏、反射镜、CCD (Charge Coupled Device)相机组成。闪烁体屏的厚度为0.4 mm，屏的有效面积为30 cm×30 cm，有效分辨率为0.5 mm。Lynx-2D在出厂前进行了标定，并提供标定功能。所获的数据在内部已经做了背景去噪、滤波等处理，输出值为10 bit的像素值。如图 3所示，左边部分表示治疗头结构，右边部分为Lynx-2D结构。当质子束流经过扫描磁铁，偏转打到Lynx-2D闪烁体屏时，与其发生相互作用。闪烁屏受粒子辐射后产生光子（波长540 nm），发光强度与入射粒子的通量在一定范围内近似成线性关系。发出的荧光经过镜子反射后，在CCD相机中成像。通过对相机记录下的束斑图像的分析，就可以获得束斑的位置及其他性质。
|
图 2 Lynx-2D实物 Figure 2 Lynx-2D |
|
图 3 点扫描束斑位置测量 Figure 3 Position measurement of point scanning beam |
使用中主要控制Lynx-2D的光圈大小和曝光时间。光圈和曝光时间的设定值主要影响每次照射接受的粒子数，如果粒子数超过一定阈值，就会引起输出值过饱和。对于Lynx-2D而言，一定范围内，照射的粒子数与输出值近似成线性关系。
如图 4所示，意大利的CNAO (Centro Nazionale di Adroterapia Oncologica)质子中心在束流能量为118.19 MeV下，测量了不同光圈下粒子数和输出值的关系[3]。从图 4可以看出，在光圈值分别为20、50、70、100的情况下，粒子数与输出值几乎都为线性关系。
|
图 4 Lynx-2D不同光圈下粒子数与输出值对应关系 Figure 4 The relationship of particles and output value in different iris size |
在点扫描照射系统调试中，Lynx-2D最主要的功能就是确定束斑的照射位置。光斑的位置点识别定位在激光精密测量中运用十分普遍，常用的方法有高斯分布拟合法、灰度重心法等[4]。因此我们主要参考这两种算法模型来分别定义Lynx-2D的图像光斑位置点，该位置点就代表束斑的位置。
3.1 高斯分布拟合法拟合法是通过对束斑离散像素点的灰度值进行拟合，从而得到一个能代表所有像素点特征的曲线或者曲面。高斯分布拟合法认为光斑的光强在任意一个垂直于光束的截面(x, y)上都满足高斯分布。因此建立二维高斯模型：
| $f(x,y) = G{\rm{exp}}[ - \frac{{{{(x - {x_0})}^2}}}{{2\sigma _x^2}} - \frac{{{{(y - {y_0})}^2}}}{{2\sigma _x^2}}]$ | (1) |
式中：f(x, y)为在截面(x, y)处的光强；G为光强幅值；(x0, y0)为曲面的极值点，可以看成为图像光斑的位置点；σx、σy是两个方向上的标准差。对式(1)两边取对数并化简为多项式可得：
| $Z = A{x^2} + B{y^2} + Cx + Dy + F$ | (2) |
可以获得如下的参数对应关系：
| $\left\{ \begin{array}{l} A = - 1/2\sigma _x^2\\ B = - 1/2\sigma _y^2\\ C = {x_0}/\sigma _x^2\\ D = {y_0}/\sigma _y^2\\ F = {\rm{In}}G - x_0^2/2\sigma _x^2 - y_0^2/2\sigma _y^2\\ Z = {\rm{In}}f(x,y) \end{array} \right.$ | (3) |
取残差：
| ${{\varepsilon }_{i}}=\left( A{'}_{i}^{2}+By{'}_{i}^{2}+Cx{{'}_{i}}+Cy{{'}_{i}}+F \right)-Z{{'}_{i}}$ | (4) |
| $Q = {\rm{min}}\sum {{{\left( {\left( {A_i^{'}2 + By_i^{'}2 + Cx_i^{'} + Cy_i^{'} + F} \right) - Z_i^{}} \right)}^2}} $ | (5) |
根据最小值条件，分别令：
| $\frac{{\partial Q}}{{\partial A}} = \frac{{\partial Q}}{{\partial B}} = \frac{{\partial Q}}{{\partial C}} = \frac{{\partial Q}}{{\partial D}} = \frac{{\partial Q}}{{\partial F}} = 0$ | (6) |
对式(5)按照式(6)分别求导，罗列方程组，最终A、B、C、D、F系数可以通过式(7)中的矩阵方程求解。
| $\left[ {\begin{array}{*{20}{c}} {x{'}_1^2}&{y{'}_1^2}&{{{x{'}}_1}}&{{{y{'}}_1}}&1\\ {x{'}_2^2}&{y{'}_2^2}&{{{x{'}}_2}}&{{{y{'}}_2}}&1\\ {x{'}_3^2}&{y{'}_3^2}&{{{x{'}}_3}}&{{{y{'}}_3}}&1\\ \vdots&\vdots&\vdots&\vdots&\vdots \\ {x{'}_n^2}&{y{'}_n^2}&{{{x{'}}_n}}&{{{y{'}}_n}}&1 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} A\\ B\\ C\\ D\\ F \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{{Z'}_1}}\\ {{{Z'}_2}}\\ {{{Z'}_3}}\\ \vdots \\ {{{Z'}_n}} \end{array}} \right]$ | (7) |
由A、B、C、D、F的系数，我们可以得到x、y方向的标准差σx、σy，由此我们可以确定束斑的大小，同时也能得到极值点(x0, y0)，由此确定束斑的位置点。高斯分布拟合法计算量较大，当束斑的畸变较大时对位置点影响较大。
3.2 灰度重心法灰度重心法是以图像光斑灰度值为权重值，然后求光斑质心。假设Lynx-2D获得的图像由M×N个像素点组成，每个像素点的坐标为(i, j)。I(i, j)表示束斑图像中第i行、第j列像素点的灰度值，则束斑的位置点可由式(8)获得：
| ${x_0} = \frac{{\sum\limits_{i = 1}^M {\sum\limits_{j = 1}^M {i \times W(i,j)} } }}{{\sum\limits_{i = 1}^M {\sum\limits_{j = 1}^M {W(i,j)} } }},\;\;\;\;\;\;\;{y_0} = \frac{{\sum\limits_{i = 1}^M {\sum\limits_{j = 1}^M {j \times W(i,j)} } }}{{\sum\limits_{i = 1}^M {\sum\limits_{j = 1}^M {W(i,j)} } }}$ | (8) |
式中：W(i, j)为各个坐标点的权重值，取W(i, j)=I(i, j)，即将图像中各像素点的灰度值看成权重值，将(x0, y0)作为束斑的位置点。
实际的照射图像上可能存在多个束斑，而上面介绍的两种方法只针对识别单个束斑的位置。因此借助图像处理的方法，先对图像做二值化、联通域标定等操作进行区域分割，获得每个束斑的局部区域，在每个特定区域中，对束斑进行单独求位置[5-6]。每个点的局部位置加上该区域的位置，就可以得到该点在整幅图像中的位置。
在此基础上利用MATLAB开发了相应的分析工具。如图 5所示，主要包括Lynx-2D数据的可视化展示、单个束斑的x、y方向上的profile、不同能量半高宽(Full Width at Half Maximum, FWHM)变化等功能[7]。
|
图 5 MATLAB Graphy User Interface分析工具 Figure 5 MATLAB Graphy User Interface analyze tool |
由于Lynx-2D获得的图像存储是10 bit，MATLAB中的图像存储是8 bit，所以在图像的处理过程中还需要做一个数值转换。
4 测量结果 4.1 位置点稳定性验证由于光斑并不是呈理想分布，因此两种方法定义的位置点并不是对应同一个点。为了验证定义的位置点是否可靠，即所定义的位置点是否能如实反映束斑的移动。我们将Lynx-2D放置在可以二维移动的平台上，如图 6所示，x、y方向移动精度均小于0.01 mm。以x方向为例，在束流轨道矫正，确保束流垂直打到Lynx-2D平面后，得到图像1，束斑位置点P1(x1, y1)。然后控制移动平台使Lynx-2D在x方向上移动距离为L，y方向保持不变，再照射得到图像2，束斑位置点坐标P2(x2, y2)。
|
图 6 移动平台 Figure 6 Moving platform |
这个过程中，束流没有偏转，即束流从治疗头出来的位置是不变的，仅移动Lynx-2D，且两个束斑是分别在两幅图像上。P1、P2都是图像上通过算法获得的坐标位置点，距离L通过移动平台获取，可以认为是束斑的移动距离，在x方向上，所测距离与实际移动的距离的差值为：
| $\Delta D = \left| {{x_2} - {x_1}} \right| - L$ | (9) |
为此如图 7所示，在感兴趣区域内x、y方向分别照射了21个点，Lynx-2D每次等间距移动10 mm，移动范围为-100~+100 mm。以上过程束流保持不变，束流轨道没有加磁场偏转，仅保持Lynx-2D移动。Lynx-2D的光圈值为50，曝光时间为2000 ms，束流能量为70 MeV。
|
图 7 照射示意图 Figure 7 Projection schematic |
图 8表示高斯分布法和灰度重心法求束斑位置的结果。横坐标表示照射点的相对间隔位置，纵坐标表示图像上所测两点距离与实际移动距离的差值ΔD（式(9)），每次都取相邻两点。可见两种方法定位精度都在Lynx-2D的有效分辨率内。
|
图 8 高斯拟合法(a)和灰度重心法(b)求束斑位置 Figure 8 Position coordinate with Gaussian fitting method (a) and grayscale center of gravity method (b) |
在此基础上，总共测了5组数据，共105个点，如表 1所示，对两种方法求得的ΔD做了均值、标准差分析，以及比较了单个束斑的平均计算时间。从表 1中发现，两者的均值和标准差相差不大，单个束斑位置点的计算时间灰度重心法要好于高斯拟合法。计算环境采用Intel(R) Core(TM) i7-6600U CPU @2.60 GHz 2.80 GHz，MATLABR2010b 64 bit。
|
表 1 高斯拟合法和灰度重心法计算结果对比 Table 1 Results comparison between Gaussian fitting method and grayscale center of gravity method |
所有光学相机镜头都存在畸变的问题，畸变属于成像的几何失真，它是由于镜头平面上不同区域对影像的放大率不同而形成的画面扭曲变形现象。这种变形的程度从画面中心至画面边缘依次递增，在画面边缘的畸变更为明显。对于Lynx-2D而言，成像的范围并不是很大，只要确保在感兴趣的范围内不存在几何失真即可。
几何失真可以通过照射一个束斑矩阵来判断，每个束斑中心点实际间距100 mm，如果在图像上分析得到的两个点的距离也为100 mm，那么可以判断不存在几何失真的情况。图 9显示了9个束斑的照射位置分布，根据图片信息获得的束斑点间距与实际照射位置间距基本相同，因此并没有几何失真的影响。
|
图 9 束斑位置分布(70 MeV) Figure 9 Image of spots pattern (70 MeV) |
图 10表示将这9个点累加后在x和y方向上的profile分布情况。纵坐标表示像素值的相对值，横坐标表示位置。从图 10可以看出，得到的束斑在x方向上的分布比y方向上的要宽。profile两侧底部的凸起是由背景噪声引起的。
|
图 10 x方向(a)和y方向(b)的profile分布 Figure 10 The profile distribution of the x direction (a) and the y direction (b) |
Lynx-2D作为测量束斑位置分布的工具，具有在线测量的优势，在治疗头的调试中极大提高了工作效率。对于束斑位置定位除了上面介绍的两种算法，还可以结合其他算法，以提高定位精度。后续还可以将Lynx-2D的软件控制系统集成到加速器的控制系统中，这样可以和束流的引出关断更加同步。除此之外，Lynx-2D也可以作为剂量验证的工具[8]，相关的功能在后续使用中开发。
| [1] |
刘世耀. 质子和重离子治疗及其装置[M]. 北京: 科学出版社, 2012, 274-276. LIU Shiyao. The proton and heavy ions therapy and its facilities[M]. Beijing: Science Press, 2012, 274-276. |
| [2] |
Tamborini A, Raffaele L, Mirandola A, et al. Development and characterization of a 2D scintillation detector for quality assurance in scanned carbon ion beams[J]. Nuclear Instruments and Methods, 2016, 815(1): 23-30. |
| [3] |
Russo S, Mirandola A, Molinelli S, et al. Charaterization of a commercial scintillation detector for 2D dosimetry in scanned proton and carbon ion beams[J]. Medical Physics, 2017, 34(1): 48-54. |
| [4] |
孔兵, 王昭, 谭玉山. 基于圆拟合的激光光斑中心检测算法[J]. 红外与激光工程, 2002, 31(3): 275-279. KONG Bing, WANG Zhao, TAN Yushan. Algorithm of laser spot detection based on circle fitting[J]. Infrared and Laser Engineering, 2002, 31(3): 275-279. |
| [5] |
Gonzalez R C, Woods R E. 数字图像处理[M]. 2版. 阮秋琦, 阮宇智, 译. 北京: 电子工业出版社, 2007. Gonzalez R C, Woods R E. Digital image processing[M]. 2nd ed. RUAN Qiuqi, RUAN Yuzhi, transl. Beijing: Publishing House of Electronics Industry, 2007. |
| [6] |
章毓晋. 图像工程-图像分析[M]. 2版. 北京: 清华大学出版社, 2005. ZHANG Yujin. Image engineering (1)-image analysis[M]. 2nd ed. Beijing: Tsinghua University Press, 2005. |
| [7] |
罗军辉. MATLAB 7.0在图像处理中的应用[M]. 北京: 机械工业出版社, 2005. LUO Junhui. MATLAB 7.0 application in the image processing[M]. Beijing: Mechanical Industry Press, 2005. |
| [8] |
Lin L Y, Ainsley C G, Mertens T, et al. A novel technique for measuring the low-dose envelope of pencil-beam scanning spot profiles[J]. Physics in Medicine and Biology, 2013, 58(14): N171-N180. |