技术领域：

本发明主要涉及重离子及质子治疗中布拉格峰(Bragg峰)展宽的模型建立及参数 计算的方法。

背景技术：

重离子及质子是目前放射治疗领域内公认的最有效的粒子，这是因为它具有显著优 于其它普通射线如X射线的生物物理学特性(Kraft G.Tumor therapy with heavy charged particles.Progress in particle and nuclear physics，2000，45：s473-544)。荷能重离子贯穿靶物 质时，主要是通过与靶原子核外电子的碰撞损失其能量，随离子能量的降低，这种碰撞 的几率增大。因此，离子在接近其射程末端时损失其大部分初始动能，形成一个高剂量 的能量损失峰，这就是Bragg峰，在其射程末端之后，即Bragg峰之后很少有剂量吸收； 同时离子在其入射通道上损失的能量较小，因而形成一个相对低剂量的坪区。这就是重 离子束特有的倒转深度剂量分布。利用这种倒转的深度剂量分布可以将其Bragg峰区瞄 准肿瘤，而使其前后及周围的健康组织受损很小。然而，肿瘤的大小通常要比Bragg峰 宽要宽，因此需要调节Bragg峰的位置，从而使峰区能够覆盖整个瘤区，使肿瘤均匀受 照。而离子束Bragg峰位的深度可以通过改变入射离子束的能量来调节，因而治疗时 Bragg峰位可精确地调整在肿瘤靶区上。

关于Bragg峰展宽的装置，国外曾有报道利用可变水柱，二进制过滤器，双楔形可 变吸收装置，环状楔形吸收装置(Wt.Chu，B.A.Ludewigt，and T.R.Renner，instrumentation for treatment of cancer using proton and light-ion beams)以及楔形过滤器等来改变入射离 子束的能量，从而实现Bragg峰位的深度调节即Bragg峰的展宽，前三种装置对Bragg 峰进行非连续的移动，第四种装置虽然能够连续移动Bragg峰，但未给出其设计参数的 计算方法，而第五种装置虽然也能展宽Bragg峰，但不易做到平滑展宽。

发明内容：

本发明的目的在于避免现有技术的不足之处，而提供一种重离子及质子治疗中 Bragg峰展宽的计算方法及其旋转降能装置。

本发明的目的可以通过采用以下技术方案来实现，一种重离子及质子治疗中Bragg峰 展宽的计算方法，其主要特点包括有如下步骤：

A.建立数学模型：

${\int }_0^T\mathrm{RBE}\left(\mathrm{LET}\right[x+{\int }_0^{T}v\left(t\right)\mathrm{dt}\left]\right)·D\left(t\right)·\mathrm{LET}[x+{\int }_0^{T}v\left(t\right)\mathrm{dt}]\mathrm{dt}=\theta \left(x\right)---\left(0\right)$

其中：T为布拉格峰(Bragg峰)移动一个周期所用的时间；v(t)为Bragg峰移动的 速度对时间的函数；LET(x)为粒子通路x上的传能线密度函数；RBE(LET(x))为相对 生物学效应，这里的数据引用的是12C6+相对于X射线的V79细胞的10％存活的相对生 物学效应；x为粒子注入深度；D(t)为单位时间内来自加速器的粒子数；函数9(x)为 t∈[0，T]时生物有效剂量分布，即为治癌所要求的Bragg峰展开后的形态。

B.用蒙特卡罗法求解模型。

C.将上述计算结果A*＝(a1 *，a2 *，a3 *，…，an *)作为设计重离子及质子治疗中Bragg峰 展宽降能装置的参数。

所述的重离子及质子治疗中Bragg峰展宽的计算方法，还包括：

$D\left(t\right)=\frac{F\left(t\right)}{I_{\mathrm{ionenergy}}}$

其中：D(t)为单位时间粒子通量；F(t)为单位时间内来自加速器的粒子能量通量； Ilonenergy为到达降能设备处的单粒子能量。

所述的重离子及质子治疗中Bragg峰展宽的计算方法，还包括： 当D(t)为常数时，

${\int }_0^T\mathrm{RBE}\left(\mathrm{LET}\right[x+{\int }_0^{T}v\left(t\right)\mathrm{dt}\left]\right)·\mathrm{LET}[x+{\int }_0^{T}v\left(t\right)\mathrm{dt}]\mathrm{dt}=\frac{1}{D\left(t\right)}·\theta \left(x\right)---\left(1\right)$

所述的重离子及质子治疗中Bragg峰展宽的计算方法，还包括：

$\begin{array}{cc}s\left(t\right)={\int }_0^{t}v\left(t\right)\mathrm{dt}& t\in [0,T]\end{array}----\left(2\right)$

s(t)为Bragg峰移动的距离时间函数。

由方程(1)和(2)得，

${\int }_0^T\mathrm{RBE}\left(\mathrm{LET}\right[x+s\left(t\right)\left]\right)·\mathrm{LET}[x+s\left(t\right)]\mathrm{dt}=\frac{1}{D}·\theta \left(x\right)---\left(3\right)$

所述的重离子及质子治疗中Bragg峰展宽的计算方法，还包括：

$s\left(t\right)=\mathrm{Broad}·\frac{a_1·t+t^{1.5}+a_2·t^2+a_3·t^3+...+a_n{·t}^n}{a_1·T+T^{1.5}+a_2·T^2+a_3·T^3+...+a_n·T^n}---\left(4\right)$

其中：Broad为Bragg峰的预期展宽宽度，α为待定参数(i＝1，2，…，n)。

所述的重离子及质子治疗中Bragg峰展宽的计算方法，还包括：Bragg峰展宽区的 相对标准误差：

$\mathrm{Relstdev}=\frac{1}{\frac{1}{n}\underset{1}{\overset{n}{\Sigma }}\theta \left(x_i\right)}·{\left\{\frac{1}{n-1}\underset{1}{\overset{n}{\Sigma }}{[\theta \left(x_i\right)-\frac{1}{n}\underset{1}{\overset{n}{\Sigma }}\theta \left(x_i\right)]}^2\right\}}^{\frac{1}{2}}---\left(5\right)$

其中：n为计算时Bragg展宽区所选计算用的点的数量。Relstdev：Bragg峰展宽区的 相对标准误差。

所述的重离子及质子治疗中Bragg峰展宽的计算方法，还包括步骤B用蒙特卡罗法 求解模型：

利用蒙特卡罗(Monte carlo)法求出使相对标准误差(Relstdev)为最小的待定系数 αi，其基本思想是：在估计的区域{(a1，a2，…，an)|ai∈(ci，bi)，i＝1，2，…，n}内随机取若干实验点， 然后从试验点中找出最小点。实验点的第i个分量ai服从[ci，bi]内的均匀分布。即：

ai＝ci+R(bi-ci)，  i＝1，2，…，n，其中ci，bi为ai取值的下限和上限；

R为一个[0，1]均匀分布的随机数。设：

P：实验点总数；

MaxP：最大实验点总数；

Broad：Bragg峰展宽宽度；

A*：迭代产生的最优点{ai *}，i＝1，2，…，n，

A：临时产生的点{ai}，i＝1，2，…，n，

Q：迭代产生的Relstdev()最小值，其初值为计算机所能表示的最大数，

Min：给定的Relstdev()的最大取值。

对迭代步骤做进一步的详述：

步骤1：初始化：给定Broad，T，MaxP；P＝0，Q大正数，

步骤2：P＝P+1，

步骤3：比较P与MaxP的大小，若P大于MaxP，退出迭代过程，输出结果，否 则，随机产生ai＝ci+R(bi-ci)，i＝1，2，…，n，

步骤4：计算相对标准误差Relstdev()，若Relstdev()大于Q，继续步骤(2)，否则， A*＝A，Q＝Relstdev()，

步骤5：比较Q与Mm，若Q大于Min，继续步骤(2)，否则退出迭代过程，输出 结果。

步骤6：结果输出：输出X*，Q，P，θ(x)。

一种重离子及质子治疗中Bragg峰展宽的旋转降能装置，包括有电机(1)通过支 架(3)固定，电机轴(2)通过在电机轴(2)上设有的键(7)与旋转降能盘(6)固连， 其主要特点是在旋转降能盘(6)的中部设有圆盘(61)，圆盘(61)的周围为降能盘(62)， 降能盘(62)的表面设有连续的曲面。

所述的重离子及质子治疗中Bragg峰展宽的旋转降能装置，所述的降能盘(62)的 表面连续的曲面为函数：

$s\left(t\right)={\int }_0^{t}v\left(t\right)\mathrm{dt\; t}\in [0,T].$

现用函数将各部分结构描述如下：

当r∈[r1，r2]时，

$\left\{\begin{array}{c}x=r·\cos \left(\alpha \right)\\ y=r·\sin \left(\alpha \right)\\ z=C+C_m·s\left(t\right)\end{array}\right.$

此实体为降能盘(62)的区域，当转子旋转时粒子沿轴方向通过该区域，达到降能 的目的。

当r∈[r0，r1)时，z＝z0，其实体为旋转降能盘(6)的中部设有的圆盘(61)，起支撑 转子的作用。

当r∈(0，r0)时，z＝0，此孔为电机轴(2)与旋转降能盘(6)连接孔，其作用是固 定转子于转动装置上。

这里，x，y为旋转降能盘(6)的横坐标和纵坐标，z为实体的高度。r，α为这些参 数方程的参数，其中r表示半径，α为转子所在极坐标系xy平面上的极角。这里参数C 为常数，有一定厚度C的介入，可以方便加工。Cm为粒子穿过该降能材料相当于穿越 单位长度水的厚度。

所述的重离子及质子治疗中Bragg峰展宽的旋转降能装置，所述的旋转降能盘(6) 的材料为有机玻璃或其它固体材料。若用有机玻璃则Cm＝0.862。

本发明的有益效果是，这一模型的建立可以实现Bragg峰的连续移动，使Bragg峰 展宽区水平平滑，实现在治癌过程中肿瘤区接收的能量达到完全均匀分布。

附图说明：

图1为本发明的参数计算程序框图。

图2为本发明Bragg峰展宽的装置实施例的结构示意图。

图3为本发明应用例Bragg峰展宽的效果示意图。

具体实施方式：

以下结合附图所示之最佳实施例作进一步详述：

入射的粒子：12C6+，

单粒子能量：100MeV/N，

靶：水(人体组织大部分为水，所以这里根据惯例以水作为靶代替人体做计算)。

实施例1：一种重离子及质子治疗中Bragg峰展宽的计算方法，包括有如下步骤：

A.建立数学模型：

所述的重离子及质子治疗中Bragg峰展宽的计算方法，主要包括：

${\int }_0^T\mathrm{RBE}\left(\mathrm{LET}\right[x+{\int }_0^{T}v\left(t\right)\mathrm{dt}\left]\right)·D\left(t\right)·\mathrm{LET}[x+{\int }_0^{T}v\left(t\right)\mathrm{dt}]\mathrm{dt}=\theta \left(x\right)---\left(0\right)$

其中：T为布拉格峰(Bragg峰)移动一个周期所用的时间(这里取T＝1s)；v(t)为 Bragg峰移动的速度对时间的函数；；LET(x)为粒子通路x上的传能线密度函数；， RBE(LET(x))为相对生物学效应，这里的数据引用的是12C6+相对于X射线的V79细胞 的10％存活的相对生物学效应。x为粒子注入深度；D(t)为单位时间内来自加速器的 粒子数；函数θ(x)为t∈[0，1]时生物有效剂量分布，即治癌所要求的Bragg峰展开后的 形态，先缓慢上升，然后迅速上升到某一点，从该点开始θ(x)为一常数，经过一段距离 后，θ(x)急剧下降；

所述的重离子及质子治疗中Bragg峰展宽的计算方法，还包括：

$D\left(t\right)=\frac{F\left(t\right)}{I_{\mathrm{ionenergy}}}$

其中：D(t)为单位时间粒子通量；F(t)为单位时间内来自加速器的粒子能量通量； Ilonenergy为到达降能设备处的单粒子能量。

所述的重离子及质子治疗中Bragg峰展宽的计算方法，还包括有：

${\int }_0^1\mathrm{RBE}\left(\mathrm{LET}\right[x+{\int }_0^{1}v\left(t\right)\mathrm{dt}\left]\right)·\mathrm{LET}[x+{\int }_0^{1}v\left(t\right)\mathrm{dt}]\mathrm{dt}=\frac{1}{D\left(t\right)}·\theta \left(x\right)---\left(1\right)$

其中：在加速器稳定的条件下D(t)通常为常数。当D(t)取常数时，因此由(0)可得(1)。

所述的重离子及质子治疗中Bragg峰展宽的计算方法，还包括：

$\begin{array}{cc}s\left(t\right)={\int }_0^{t}v\left(t\right)\mathrm{dt}& t\in \left[\mathrm{0,1}\right]\end{array}---\left(2\right)$

s(t)为Bragg峰移动的距离时间函数：由方程(1)和(2)得，

${\int }_0^1\mathrm{RBE}\left(\mathrm{LET}\right[x+s\left(t\right)\left]\right)·\mathrm{LET}[x+s\left(t\right)]\mathrm{dt}=\frac{1}{D\left(t\right)}·\theta \left(x\right)---\left(3\right)$

所述的重离子及质子治疗中Bragg峰展宽的计算方法，还包括：

$s\left(t\right)=\mathrm{Broad}·\frac{a_1·t+t^{1.5}+a_2·t^2+a_3·t^3}{a_1+1+a_2+a_3}---\left(4\right)$

其中：Broad为Bragg峰的预期展宽宽度，αi为待定参数(这里i＝1，2，3)；

所述的重离子及质子治疗中Bragg峰展宽的计算方法，还包括：Bragg峰展宽区的 相对标准误差：

$\mathrm{Relstdev}=\frac{1}{\frac{1}{n}\underset{1}{\overset{n}{\Sigma }}\theta \left(x_i\right)}·{\left\{\frac{1}{n-1}\underset{1}{\overset{n}{\Sigma }}{[\theta \left(x_i\right)-\frac{1}{n}\underset{1}{\overset{n}{\Sigma }}\theta \left(x_i\right)]}^2\right\}}^{\frac{1}{2}}---\left(5\right)$

其中：n为计算时Bragg峰展宽区所选计算用的点的数量，这里取n＝87。Relstdev： Bragg峰展宽区的相对标准误差。所述的重离子及质子治疗中Bragg峰展宽的计算方法， 其特征还包括所述公式(5)为优化函数求解所述公式(4)的待定系数αi(i＝1，2，3)。

B.模型求解

利用蒙特卡罗(Monte carlo)法求出使相对标准误差(Relstdev)为最小的待定系数 αi(i＝1，2，3)，

见图1，对迭代步骤做进一步的详述：

步骤1：初始化Broad，T，MaxP，Q为Broad＝1，T＝1，MaxP＝27000，Q＝1。当αi∈[0，1]， min＝0.01时，

步骤2：P＝P+1，

步骤3：比较P与MaxP的大小，若P大于MaxP，退出迭代过程，输出结果，否 则，随机产生ai＝ci+R(bi-ci)，i＝1，2，…，n，

步骤4：计算相对标准误差Relstdev()，若Relstdev()大于Q，继续步骤(2)，否则， A*＝A，Q＝Relstdev()；

步骤5：比较Q与Min，若Q大于Min，继续步骤(2)，否则退出迭代过程，输出 结果。

步骤6：结果输出：输出X*，Q，P，θ(x)。

求得：α1＝0.071627，α2＝0.29017，α3＝0.001832，相对标准误差为Relstedev＝Q＝0.012

通常当相对误差小于5％时就可以用于治疗.。峰坪比为3.5。

因此可得，

$s\left(t\right)=\frac{0.071627·t+t^{1.5}+0.29017·t^2+0.001832·t^3}{0.071627+1+0.29017+0.001832}=\frac{0.071627·t+t^{1.5}+0.29017·t^2+0.001832·t^3}{1.363629}$

C.设备设计

见图2，一种重离子及质子治疗中Bragg峰展宽的旋转降能装置，包括有电机1通 过支架3固定，电机轴2通过在电机轴2上设有的键7与旋转降能盘6固连，其主要特 点是在旋转降能盘6的中部设有圆盘61，圆盘61的周围为降能盘62，降能盘62的表面 设有连续的曲面。

所述的重离子及质子治疗中Bragg峰展宽的旋转降能装置，所述的降能盘62的表 面连续的曲面为函数：

$s\left(t\right)={\int }_0^{t}v\left(t\right)\mathrm{dt\; t}\in [0,T].$

这里，根据兰州重离子加速器治癌终端为例进行设计，该终端束流直径为45mm，这 里取r0＝6mm，r1＝115mm，r2＝175mm，因此，设计的展宽装置降能区的径向宽度为：r1- r2＝60mm，大于45mm，这样能保证束斑能落在整个降能区。所述的旋转降能盘6的材 料为有机玻璃，因此Cm＝0.862。参数C根据加工条件选取，这里取理想值C＝0。

当r∈[115，175]时，

$\left\{\begin{array}{c}\\ x=r·\cos \left(\alpha \right)\\ y=r·\sin \left(\alpha \right)\\ z=0.862·\frac{0.071627·t+t^{1.5}+0.29017·t^2+0.001832·t^3}{1.363629}\end{array}\right.$

此实体为降能所使用的区域，当转子旋转时粒子沿轴方向通过该区域，达到降能的 目的。

此实体为降能盘62的区域，当转子旋转时粒子沿轴方向通过该区域，达到降能的目 的。降能盘62表面为连续的曲面，0-π/2，π-3π/2为上升曲面；π/2-π，3π/2 -2π为下降曲面。

当r∈[6mm，115mm)时，z＝10mm，其实体为旋转降能盘6的中部设有的圆盘61， 起支撑转子的作用。

当r∈(0，6mm)时，z＝0mm，此孔为电机轴2与旋转降能盘6连接孔，其作用是固 定转子于转动装置上。

在使用时，病人躺在圆盘下，束流线通过降能盘照射到人体癌细胞处。粒子沿轴方 向穿过变化的曲面体，旋转降能盘6以大于150转/分的速度旋转。此值无需计算。

展宽的结果见图3，其中RBE为加入相对生物学效应的展宽结果，Physical为展宽 后物理剂量分布。Pristine为单能Bragg峰。展宽后的结果在图中已归一化，即，使入射能 量为相对值1。