CN113230549A - 医用直线加速器光子源模型函数的构建方法 - Google Patents

医用直线加速器光子源模型函数的构建方法 Download PDF

Info

Publication number
CN113230549A
CN113230549A CN202110645477.5A CN202110645477A CN113230549A CN 113230549 A CN113230549 A CN 113230549A CN 202110645477 A CN202110645477 A CN 202110645477A CN 113230549 A CN113230549 A CN 113230549A
Authority
CN
China
Prior art keywords
photons
fluence
function
distribution
bremsstrahlung
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202110645477.5A
Other languages
English (en)
Other versions
CN113230549B (zh
Inventor
白彦灵
崔振国
王业伟
刘奇
冯丽娜
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Harbin Engineering University
Harbin Medical University
Original Assignee
Harbin Medical University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Harbin Medical University filed Critical Harbin Medical University
Priority to CN202110645477.5A priority Critical patent/CN113230549B/zh
Publication of CN113230549A publication Critical patent/CN113230549A/zh
Application granted granted Critical
Publication of CN113230549B publication Critical patent/CN113230549B/zh
Priority to US17/828,070 priority patent/US20220401754A1/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61N—ELECTROTHERAPY; MAGNETOTHERAPY; RADIATION THERAPY; ULTRASOUND THERAPY
    • A61N5/00—Radiation therapy
    • A61N5/10—X-ray therapy; Gamma-ray therapy; Particle-irradiation therapy
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61N—ELECTROTHERAPY; MAGNETOTHERAPY; RADIATION THERAPY; ULTRASOUND THERAPY
    • A61N5/00—Radiation therapy
    • A61N5/10—X-ray therapy; Gamma-ray therapy; Particle-irradiation therapy
    • A61N5/103—Treatment planning systems
    • A61N5/1031—Treatment planning systems using a specific method of dose optimization
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61N—ELECTROTHERAPY; MAGNETOTHERAPY; RADIATION THERAPY; ULTRASOUND THERAPY
    • A61N5/00—Radiation therapy
    • A61N5/10—X-ray therapy; Gamma-ray therapy; Particle-irradiation therapy
    • A61N5/1001—X-ray therapy; Gamma-ray therapy; Particle-irradiation therapy using radiation sources introduced into or applied onto the body; brachytherapy
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61N—ELECTROTHERAPY; MAGNETOTHERAPY; RADIATION THERAPY; ULTRASOUND THERAPY
    • A61N5/00—Radiation therapy
    • A61N5/10—X-ray therapy; Gamma-ray therapy; Particle-irradiation therapy
    • A61N5/1048—Monitoring, verifying, controlling systems and methods
    • A61N5/1064—Monitoring, verifying, controlling systems and methods for adjusting radiation treatment in response to monitoring
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/25—Design optimisation, verification or simulation using particle-based methods
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16H—HEALTHCARE INFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR THE HANDLING OR PROCESSING OF MEDICAL OR HEALTHCARE DATA
    • G16H20/00—ICT specially adapted for therapies or health-improving plans, e.g. for handling prescriptions, for steering therapy or for monitoring patient compliance
    • G16H20/40—ICT specially adapted for therapies or health-improving plans, e.g. for handling prescriptions, for steering therapy or for monitoring patient compliance relating to mechanical, radiation or invasive therapies, e.g. surgery, laser therapy, dialysis or acupuncture
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16H—HEALTHCARE INFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR THE HANDLING OR PROCESSING OF MEDICAL OR HEALTHCARE DATA
    • G16H40/00—ICT specially adapted for the management or administration of healthcare resources or facilities; ICT specially adapted for the management or operation of medical equipment or devices
    • G16H40/60—ICT specially adapted for the management or administration of healthcare resources or facilities; ICT specially adapted for the management or operation of medical equipment or devices for the operation of medical equipment or devices
    • G16H40/63—ICT specially adapted for the management or administration of healthcare resources or facilities; ICT specially adapted for the management or operation of medical equipment or devices for the operation of medical equipment or devices for local operation
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16H—HEALTHCARE INFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR THE HANDLING OR PROCESSING OF MEDICAL OR HEALTHCARE DATA
    • G16H50/00—ICT specially adapted for medical diagnosis, medical simulation or medical data mining; ICT specially adapted for detecting, monitoring or modelling epidemics or pandemics
    • G16H50/50—ICT specially adapted for medical diagnosis, medical simulation or medical data mining; ICT specially adapted for detecting, monitoring or modelling epidemics or pandemics for simulation or modelling of medical disorders
    • H—ELECTRICITY
    • H05—ELECTRIC TECHNIQUES NOT OTHERWISE PROVIDED FOR
    • H05H—PLASMA TECHNIQUE; PRODUCTION OF ACCELERATED ELECTRICALLY-CHARGED PARTICLES OR OF NEUTRONS; PRODUCTION OR ACCELERATION OF NEUTRAL MOLECULAR OR ATOMIC BEAMS
    • H05H9/00—Linear accelerators
    • H05H9/04—Standing-wave linear accelerators
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61N—ELECTROTHERAPY; MAGNETOTHERAPY; RADIATION THERAPY; ULTRASOUND THERAPY
    • A61N5/00—Radiation therapy
    • A61N5/10—X-ray therapy; Gamma-ray therapy; Particle-irradiation therapy
    • A61N5/103—Treatment planning systems
    • A61N5/1031—Treatment planning systems using a specific method of dose optimization
    • A61N2005/1034—Monte Carlo type methods; particle tracking
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61N—ELECTROTHERAPY; MAGNETOTHERAPY; RADIATION THERAPY; ULTRASOUND THERAPY
    • A61N5/00—Radiation therapy
    • A61N5/10—X-ray therapy; Gamma-ray therapy; Particle-irradiation therapy
    • A61N2005/1085—X-ray therapy; Gamma-ray therapy; Particle-irradiation therapy characterised by the type of particles applied to the patient
    • A61N2005/1091—Kilovoltage or orthovoltage range photons
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61N—ELECTROTHERAPY; MAGNETOTHERAPY; RADIATION THERAPY; ULTRASOUND THERAPY
    • A61N5/00—Radiation therapy
    • A61N5/10—X-ray therapy; Gamma-ray therapy; Particle-irradiation therapy
    • A61N2005/1092—Details
    • A61N2005/1095—Elements inserted into the radiation path within the system, e.g. filters or wedges
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/08—Probabilistic or stochastic CAD
    • G—PHYSICS
    • G21—NUCLEAR PHYSICS; NUCLEAR ENGINEERING
    • G21G—CONVERSION OF CHEMICAL ELEMENTS; RADIOACTIVE SOURCES
    • G21G1/00—Arrangements for converting chemical elements by electromagnetic radiation, corpuscular radiation or particle bombardment, e.g. producing radioactive isotopes
    • G21G1/04—Arrangements for converting chemical elements by electromagnetic radiation, corpuscular radiation or particle bombardment, e.g. producing radioactive isotopes outside nuclear reactors or particle accelerators
    • G21G1/12—Arrangements for converting chemical elements by electromagnetic radiation, corpuscular radiation or particle bombardment, e.g. producing radioactive isotopes outside nuclear reactors or particle accelerators by electromagnetic irradiation, e.g. with gamma or X-rays

Landscapes

  • Health & Medical Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Biomedical Technology (AREA)
  • Public Health (AREA)
  • General Health & Medical Sciences (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Pathology (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Animal Behavior & Ethology (AREA)
  • Radiology & Medical Imaging (AREA)
  • Veterinary Medicine (AREA)
  • Medical Informatics (AREA)
  • Physics & Mathematics (AREA)
  • Primary Health Care (AREA)
  • Epidemiology (AREA)
  • Theoretical Computer Science (AREA)
  • Plasma & Fusion (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Computer Hardware Design (AREA)
  • General Physics & Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Databases & Information Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Geometry (AREA)
  • Evolutionary Computation (AREA)
  • Business, Economics & Management (AREA)
  • General Business, Economics & Management (AREA)
  • Surgery (AREA)
  • Urology & Nephrology (AREA)
  • Radiation-Therapy Devices (AREA)
  • Analysing Materials By The Use Of Radiation (AREA)

Abstract

本发明涉及辐射剂量测量技术领域，提出了一种医用直线加速器光子源模型函数的构建方法，用于放射治疗方案中射线的剂量计算，该加速器的治疗光子束的源模型包含原射线光子和散射线光子的源模型两部分。这两部分的源模型函数中的物理参量为粒子的发射点位置坐标、单位动量矢量三维正交方向上的投影值、粒子的能量。利用该模型函数能够准确计算出在任意相空间平面上的光子注量信息、能谱信息、光子的单位动量方向信息，本发明涉及的源模型构建方法和思路，适用于放射治疗中所使用的加速器的各种标称能量的光子束源模型构建，具有建模参数少，函数中各个参数的物理意义明确、使用方便，能够用数学解析方法如实复现辐射场光子分布信息的优点。

Description

医用直线加速器光子源模型函数的构建方法
技术领域
本发明涉及辐射剂量测量技术领域，尤其涉及一种医用直线加速器光子源模型函数的构建方法。
背景技术
放射治疗是依赖于高能电离射线治疗肿瘤的方法，它和外科手术治疗方法、内科药物治疗方法一起构成了肿瘤治疗的三大主流手段。在放射治疗中，肿瘤细胞的杀灭概率和正常组织的损伤程度都与射线剂量的高低有关，因此，放射治疗在照射方案设计时，需要将射线剂量从不同方向上集中照射到肿瘤上，分散正常组织的剂量，最大可能提高肿瘤区域的剂量、降低正常组织器官的剂量、尽可能减少正常组织受照体积。因此，在放射治疗方案设计阶段时，剂量计算的准确性很重要。
放射治疗已经有了近百年历史，它的剂量计算方法或计算模型有很多，也在不断进步和发展。今天临床上应用最广泛的剂量计算都是数学解析方法，诸如笔形束算法、AAA算法、XB算法以及迭代卷积算法等等，但是这些算法在非均匀组织中的情况下，特别是密度差异比较大的区域内，计算的误差较大，误差达到4％-17％。
而蒙特卡洛剂量计算方法是目前公认的准确度、精度最高的计算方法，它基于统计物理的抽样数学方法，通过统计物理的理论、原子物理和核物理理论的数学方法可以精准模拟出粒子在介质中的输运过程，并可以记录粒子在输运过程中的各种信息。蒙特卡洛的剂量计算过程按照物理阶段可以分成两个阶段，包括：(1)模拟追踪打靶的电子，计算追踪电子与靶原子相互作用、相互作用过程中产生的轫致辐射光子以及轫致辐射光子被均整器修整的过程，记录修整后的治疗射线束的光子粒子信息；(2)追踪治疗射线束的光子在病人体内的输运过程，并记录输运过程中的信息。在实际治疗场景下，治疗的病人可以轮换更替变化，但医用直线加速器光子源装置是不变的，所以，预处理阶段在加速器机头内部空间的合适位置处选择一个垂直射线束的平面(一般选择加速器射线束均整齐下方)，记录到达这个平面上的治疗射线束光子的信息，这些光子信息包括了其在该平面上的位置坐标、运动动量方向、能量，记录的光子信息被存储在一个文件中，这种文件也称之为相空间文件，这个平面也称之为相空间平面。当需要做病人剂量计算的时候，这个相空间文件就作为粒子的相空间平面源输入文件直接用于剂量计算，可省去第一阶段加速器机头内的模拟过程，节省大量计算时间、提高工作效率。
但由于相空间文件必须记录足够多的粒子(数十亿)才能保证相空间平面的粒子注量的方差足够小，从而会导致文件的尺寸过大，占据很大的硬盘空间，而且在剂量计算中，从硬盘把粒子信息导入到计算程序中也十分的耗时；在相空间平面上的粒子注量的固有统计方差也会传递到后续的剂量计算过程中；另外，相空间文件一旦形成，记录的粒子数目、分布尺寸不能改变，不能跟随照射野而改变。因此，目前用数学函数描述方法来替代相空间是最好的选择，通过数学函数中的参量准确地表达粒子源物理特征，包括源发射的粒子能量特征、粒子发射点的空间位置、粒子的动量方向分布特征，这个数学函数直接被用作蒙特卡洛剂量计算粒子输入源的描述函数，就不再需要相空间文件，节省了时间、方便了计算，能显著改善了计算效率和方差。
现有技术中，供蒙特卡洛剂量计算使用的医用直线加速器光子剂量计算的源模型函数的建立方法：是利用靶位置的光子以及相空间的光子落点位置建立基于空间几何投影关系，然后去建立光子物理参数的数学形式，引入光子环的概念来建立源模型函数。
建立源模型函数过程中，来自于钨靶位置处的光子被称为原射线光子，这些原射线光子有不同的发射方向(即，动量方向)，按照动量方向对这些原射线光子分组后，它们在任意的相空间平面上会形成环带分布，这个特点就是所谓的光子环。因此，通过蒙特卡洛计算获取的相空间平面上的光子，按照动量方向分组后，会出现一环一环的高亮分布区，分布在这个高亮环形区域内的光子为原射线光子，而剩余的为散射线光子，没有聚焦区域，如图1所示。
如图2所示，高能电子束流引出后垂直击打薄钨靶，产生X线光子，任意有效位置点处的光子数量与击打该点的电子束流强度(数量)成正比。建立的数学模型如下：
Figure BDA0003109382900000021
Figure BDA0003109382900000022
Figure BDA0003109382900000023
(xs，ys)表示为光子在靶平面的位置坐标，r表示光子在相空间平面的落点距离中心的距离；(u，v，w)是光子的单位动量方向的三个分量投影；δs是击靶的电子束流分布方差。δ(u-u0)或δ(v-v0)为冲击函数，(u0，v0)为落到辐射场平面上的光子动量方向分量；zphsp为相空间平面(辐射场平面)到源的距离。
散射光子源模型函数：由于散射光子的动量方向是杂乱无章的，但是光子的动量分量满足u2+v2+w2＝1,是一个球面，所以根据散射光子按照动量分布分组后散射光子在半球面上分布如图3所示。按照其概率分布特点进行了方向(u和v)的拟合，不同的w所对应的散射光子的分布概率都是一个高斯函数，这个半球面上的散射光子可以拟合出很多高斯函数。基于拟合结果以及散射光子的角度，在相空间平面的半径区域，散射光子的动量分布为：
Figure BDA0003109382900000031
其中，α和β分别为动量方向投影u和v对应的方向角度。
但是，在上述模型函数中，至少存在以下缺点：
(1)有源平面的光子数量分布信息为
Figure BDA0003109382900000032
靶发射出来的光子经过均整器过滤后，光子注量的分布函数已经与靶发射的光子注量不同了，因此不能由单一高斯函数描述它的分布；
(2)光子在动量方向信息采用光子分布半径方式提供，均匀分布在这个环带上的光子具有相同的动量方向投影w。但由于加速器机头内部不可或缺的均整器的存在，使得在靶平面上的韧致辐射光子的注量分布与受均整器影响后的原射线光子注量分布形式不同；
(3)可确定光子的发射方向或在计算平面上的落点位置，但是不能提供原射线光子在这个方向的注量概率，或者沿光子动量方向上的光子注量概率，因此不能用于蒙特卡洛剂量剂量中的源模型；
(4)受均整器作用后的原射线光子分布函数不是一个简单的二维高斯函数，在本模型中，源平面上的光子注量分布函数表达如果是加速器靶产生光子分布，那么缺少了均整器的作用；如果是受均整器影响之后的原射线光子分布形式表达，这个函数不是简单高斯函数，因为靶产生的韧致辐射光子经过均整器之后的粒子注量分布已经不再是钟形，不能如实反映光子的注量分布信息，不适合计算复杂照射野、或细致野的剂量；
(5)该模型函数描述散射光子的函数很复杂，通过函数计算获取散射光子的位置、动量方向分量信息繁琐；
(6)上述模型函数缺失了不同环形带上粒子分布权重信息。
发明内容
因此，针对上述现有技术中存在的技术问题，本发明提出了一种新的构建医用直线加速器光子源模型函数的方法，用于放射治疗方案中射线的剂量计算，具有建模参数少，函数中各个参数的物理意义明确、使用方便，能够在任意位置高度的辐射场内如实复现光子注量分布信息、光子的单位动量方向信息、能量信息的优点。
具体的，主要通过以下技术方案来实现：
一种医用直线加速器光子源模型函数的构建方法，医用直线加速器光子源模型包括原射线光子源模型和散射线光子源模型两部分，该方法包括：
开启医用直线加速器，用产生的高能电子束轰击所述医用直线加速器的X线靶，将初始相空间平面上记录的靶面上产生的轫致辐射光子信息保存在第一相空间文件中，以及，将相空间平面上记录的经均整滤过器后的治疗光子束光子信息保存在第二相空间文件中；
通过光子逆向飞行计算，将第一相空间文件和第二相空间文件转换为靶平面上的相空间文件，分别记为B-PhSp和R-PhSp；将相空间文件R-PhSp中记录的光子分离为原射线光子和散射线光子，得到两个新的相空间文件分别记为p-PhSp和s-PhSp；
根据所述轫致辐射光子注量分布函数、靶面上注量点发射韧致辐射光子在动量分量上的注量概率分布函数和靶面上注量点发射韧致辐射光子在动量分量方向上被所述均整滤过器吸收的注量概率分布函数构建原射线光子源模型函数；
根据所述散射线光子注量分布函数和靶面上各注量点发射的散射光子在动量分量上的平均注量概率分布函数构建散射线光子源模型函数。
优选地，将相空间文件R-PhSp中记录的光子分离为原射线光子和散射线光子，具体包括：相空间文件R-PhSp记录的光子中，将位置坐标位于初级准直器的圆锥形孔径在靶底端面上的几何投影区域之内的光子记为原射线光子，所述原射线光子可沿着初级准直器孔径回到所述初始相空间平面上；将分布在所述几何投影区域之外的光子记为散射线光子。
优选地，将所述相空间文件B-PhSp和相空间文件p-PhSp中所记录的光子的位置坐标参量置于靶平面中心点(0，0，0)，其中靶平面设为0参考平面。其等价于击靶的电子数量提高数倍时靶面中心点产生的韧致辐射光子数量和此时原射线光子束中来自靶中心心点的光子数量。
优选地，加速器光子源模型函数表达的光子是光子束的能谱中的一个谱线能量的光子，所以，用能量宽度BinE将连续能谱的韧致辐射光子分成若干能量组，用Ei表示，i为自然数，表示第i组光子的能量分布中值为Ei，亦可称之为第i条谱线能量Ei，该组光子的能量分布宽度为
Figure BDA0003109382900000041
优选地，所述构建原射线光子源模型函数的步骤，具体包括：
所述轫致辐射光子注量分布函数为：
Figure BDA0003109382900000042
以及，
Figure BDA0003109382900000043
Figure BDA0003109382900000051
其中，
Figure BDA0003109382900000052
为
Figure BDA0003109382900000053
在NB(xs，ys)的分支比；
Figure BDA0003109382900000054
表示能量为Ei的韧致辐射光子总计数，NB表示所有能量的韧致辐射光子的总计数；NB0表示为靶中心点产生的韧致辐射光子计数；
根据韧致辐射光子在动量分量方向的注量分布，并拟合得到韧致辐射光子在动量分量方向的注量分布函数为：
Figure BDA0003109382900000055
其中，
Figure BDA0003109382900000056
表示在靶中心点产生的第i能量组的韧致辐射光子沿动量的分量方向的注量分布；
Figure BDA0003109382900000057
和
Figure BDA0003109382900000058
为韧致辐射光子中能量为Ei的光子沿动量分量方向上的注量分布标准差，
Figure BDA0003109382900000059
和
Figure BDA00031093829000000510
是函数表达式的系数，等于0度角方向上发射的韧致辐射光子计数，韧致辐射光子飞行的动量方向用单位动量在三维直角坐标系上的投影值表示，即(px，py，pz)，并且
Figure BDA00031093829000000511
B表示韧致辐射光子，Ei表示第i条谱线能量，m表示动量(momentum)，即，
Figure BDA00031093829000000512
和
Figure BDA00031093829000000513
为动量分布的两个标准差。
因为靶面上各个注量点产生的韧致辐射光子的角注量概率一样，则，靶面上注量点发射韧致辐射光子在动量分量上的注量概率分布函数为：
Figure BDA00031093829000000514
其中，
Figure BDA00031093829000000515
和
Figure BDA00031093829000000516
是韧致辐射光子沿动量分量方向的注量概率分布函数表达式的系数，根据
Figure BDA00031093829000000517
和
Figure BDA00031093829000000518
的归一计算得到；
根据轫致辐射光子沿动量的分量方向的注量分布与原射线光子沿动量的分量方向的注量分布之差，得到轫致辐射光子在动量分量方向上被所述均整滤过器吸收的注量分布，并拟合得到轫致辐射光子在动量分量方向上被所述均整滤过器吸收的注量分布函数：
Figure BDA00031093829000000519
其中，
Figure BDA00031093829000000520
表示在动量的分量方向(px，py)上第i条谱线能量的韧致辐射光子中被均整滤过器吸收的注量分布；
Figure BDA00031093829000000521
和
Figure BDA00031093829000000522
分别为上述两个二维高斯函数的标准差；
Figure BDA0003109382900000061
和
Figure BDA0003109382900000062
分别是上述两个高斯函数表达式的系数；
将任一注量点处第i条谱线能量的韧致辐射光子总计数做数学归一处理，得到靶面上注量点发射韧致辐射光子在动量分量方向上被所述均整滤过器吸收的注量概率分布函数为：
Figure BDA0003109382900000063
其中，
Figure BDA0003109382900000064
是韧致辐射光子被均整滤过器吸收的概率；
Figure BDA0003109382900000065
和
Figure BDA0003109382900000066
分别是该概率分布模型函数的两个系数，根据
Figure BDA0003109382900000067
和
Figure BDA0003109382900000068
的数学归一计算获得；
根据所述靶面上注量点发射韧致辐射光子在动量分量上的注量概率分布函数与靶面上注量点发射韧致辐射光子在动量分量方向上被所述均整滤过器吸收的注量概率分布函数之差，得到靶面上注量点出射的原射线光子在动量分量上的注量概率分布函数，因此，再结合所述轫致辐射光子注量分布函数，得到原射线光子源模型函数为：
Figure BDA0003109382900000069
其中：
Figure BDA00031093829000000610
以及
Figure BDA00031093829000000611
以及
Figure BDA00031093829000000612
其中，原射线光子源模型位于0参考平面上，且原射线光子飞行的动量方向用单位动量在三维直角坐标系上的投影值表示，即(px，py，pz)，并且
Figure BDA00031093829000000613
因此，略去zs和pz；上述函数式中各个参量均由拟合函数的过程中给出，
Figure BDA00031093829000000614
表示韧致辐射光子在动量分量上的注量概率分布；
Figure BDA00031093829000000615
是韧致辐射光子被均整滤过器吸收的概率；
Figure BDA0003109382900000071
为韧致辐射光子在动量分量上的注量概率分布函数的系数，
Figure BDA0003109382900000072
为韧致辐射光子在动量分量上的注量概率分布函数的分布标准差；
Figure BDA0003109382900000073
是韧致辐射光子在动量分量上的注量概率函数的系数，
Figure BDA0003109382900000074
是韧致辐射光子被均整滤过器吸收的概率函数的标准差，(f-absorbed)是被均整滤过器吸收的标识。
优选地，所述构建散射线光子源模型函数的步骤，具体包括：
将散射光子的相空间文件s-PhSp中所记录的散射光子的分布区域进行离散化处理，形成像素格，并以像素格为单位获取0参考平面上的散射光子的注量分布，用双高斯函数多项式对所述散射光子的注量分布进行拟合，得到散射光子的注量分布函数为：
Figure BDA0003109382900000075
其中，(xs，ys)为0参考平面上的任一点位置坐标，
Figure BDA0003109382900000076
是在注量点(xs，ys)的散射光子数，ai、bi是能量为Ei的散射光子的注量分布的拟合函数的系数；
Figure BDA0003109382900000077
是能量为Ei的散射光子的注量分布的拟合函数的分布标准差；参量下标中的s表示散射线标识；
靶所在的平面上各注量点发射的散射光子在动量分量上的注量分布的拟合函数为：
Figure BDA0003109382900000078
上式中，参量a'i、b'i、c'i、d'i为动量分量上的注量分布的拟合函数的系数；
Figure BDA0003109382900000079
为动量分量上的注量分布的拟合函数的标准差。各个参量的下标中，i表示第i条谱线能量组标识；
则，根据散射光子在动量分量上的注量分布函数得到散射光子在动量分量上的平均注量概率分布函数为：
Figure BDA00031093829000000710
其中，函数定义域r可用下面数学方法给出：
Figure BDA00031093829000000711
散射光子沿动量分量方向的平均注量概率分布函数中的参量ai、bi、ci、di、
Figure BDA0003109382900000081
由拟合过程给出，其中ai、bi、ci、di为函数的系数，
Figure BDA0003109382900000082
为沿动量分布的分布标准差，这些参量与能量密切相关，可以对能量拟合，分布形成各个参量与能量的关系函数；
根据散射光子的注量分布函数和靶面上各注量点发射的散射光子在动量分量上的平均注量概率分布函数构建散射光子源模型函数为：
Figure BDA0003109382900000083
本发明相较于现有技术具有以下有益效果：
本发明公开的模型函数的建模参数少，所表达的物理意义清晰，各个参数的物理意义明确、使用方便，直接明了地包含了加速器光子治疗束中光子完备的特征值信息，这些信息是光子的发射初始点位置坐标、粒子的单位动量方向在正交坐标系中的投影以及粒子的能量；利用数学函数直观地反映了在源平面上光子注量的空间分布函数、光子注量的角度分布函数，利用该模型函数能够直接地、准确地计算出加速器的治疗光子束在任意相空间平面上的光子注量信息、能谱信息、光子的单位动量方向信息，具有用数学解析方法如实复现辐射场光子分布信息的优点；可根据需要任意设定模拟计算的光子样本量；采用本模型函数计算任意相空间平面上的光子分布不存在统计涨落方差，这将大幅改善剂量计算的方差和收敛情况；本模型函数存储占用字节数极小，甚至可以忽略，因此有利于传播使用。
本发明涉及的源模型构建方法和思路，适用于放射治疗中所使用的各种标称能量的加速器的光子源模型构建，可用于加速器光子治疗束的蒙特卡洛剂量计算及相关产品开发。
附图说明
1、图1示出了现有源模型函数建立过程中的原射线光子模型，a示出了相空间平面上形成的光子分布环，b示出了相空间平面上光子位置以及动量方向；
2、图2为现有源模型函数建立过程中的原射线光子模型的几何示意图；
3、图3为现有源模型函数中根据散射光子按照动量分布分组后散射光子在半球面上的分布图；
4、图4为本发明实施例中韧致辐射光子的动量角分布剖面图；
5、图5为本发明实施例中被均整器吸收的光子在靶平面上的动量角分布剖面图；
6、图6为本发明实施例中散射光子在靶平面上的空间注量分布剖面图；
7、图7为本发明实施例中散射光子在靶平面上的动量角分布剖面图；
8、图8为本发明实施例中韧致辐射光子在动量分量上的注量概率分布剖面图与模型函数计算结果比较；
9、图9为本发明实施例中韧致辐射光子在动量分量上被均整器吸收的注量概率分布与模型函数计算结果的比较；
10、图10为本发明实施例中散射光子在0平面上的空间注量分布与模型函数计算结果比较；
11、图11为本发明实施例中散射光子在动量分量上的概率注量分布与模型函数计算结果比较；
12、图12为本发明实施例中初始韧致辐射光子在动量分量方向上分布的模型函数的系数及标准差与光子能量关系；
13、图13为本发明实施例中被均整器吸收的韧致辐射光子在动量分量方向上分布的模型函数的系数及标准差与光子能量关系；
14、图14为本发明实施例中散射光子源模型函数的系数与光子能量关系；
15、图15为本发明实施例中医用直线加速器产生治疗束光子的物理流程图；
16、图16为本发明实施例中医用直线加速器光子源模型函数的构建方法的流程图。
具体实施方式
为了使本领域技术人员更清楚的理解本发明的核心思想，下面将结合附图对其进行详细的说明。
本发明所提出的医用直线加速器光子源模型函数的构建方法主要用于放射治疗方案中射线的剂量计算，以及用于蒙特卡洛方法的计算粒子源。所涉及的医用直线加速器包括X线靶、射线源装置、加速管、均整滤过器以及准直器系统，其中准直器系统包括初级准直器系统和二级准直器系统。所述医用直线加速器产生治疗束光子的物理流程如图15所示：通过射线源装置的电子枪灯丝加热后发射的热电子被聚焦成直径1-2mm的电子束流，将电子束流引入到加速管中并被微波电场加速，加速到预定能量高度后引出，被引出出射窗的高能电子束轰击加速器的x线钨靶(厚度被视为薄靶，厚度不及1mm，可为0.89mm)，在靶平面上产生大量韧致辐射光子，韧致辐射光子束经过初级准直器的锥形孔、均整滤过器后，成为治疗所用的光子束。钨靶镶嵌在铜质靶衬中，铜质靶衬的冷却系统将带走电子束流轰击靶的过程中产生的大量热量。需要说明的是，上述加速器的治疗所用的光子束中既包含大量来自靶的光子，这部分光子即称之为原射线光子；也含有来自韧致辐射光子与靶衬、射线束均整滤过器相互作用过程中不断地与均整器发生相互作用而产生光电效应、康普顿效应、对效应以及瑞利散射效应，导致这部分初始轫致辐射光子在经过均整器的过程中，其初始运动方向和能量发生了改变，从而产生的大量散射光子。原射线光子可在空间上的靶平面上形成聚焦斑点，即加速器的靶面；散射光子没有聚焦斑点。
本发明的医用直线加速器的光子源模型函数的构建采用双源模型方式，所述光子源模型包含两部分：原射线光子源模型和散射光子源模型，但这两部分的源模型函数的参量为光子的状态信息，即光子的位置坐标(xs，ys，zs)、单位动量在正交坐标系下的投影(px，py，pz)、能量E，其中，光子的位置坐标采用直角坐标系；光子的动量方向采用光子单位动量在直角坐标系中的投影值刻画，即动量分量，表示为px、py、pz,且
Figure BDA0003109382900000101
原射线光子源模型与散射光子源模型之间存在明确的、可获知的比例关系。由于加速器的钨靶产生的巨量韧致辐射光子分布在一个连续能谱范围内，可离散化为一系列谱线，每个谱线能量用Ei表示，每一条谱线中有大量光子。
如图16所示为本发明公开的医用直线加速器光子源模型函数的构建方法流程示意图，并结合图4-14对构建医用直线加速器光子源模型函数的方法进行详细说明。其中，医用直线加速器以4MV光子束为例，关于坐标系：正交直角坐标系(x,y,z,)，z轴方向与靶中心点射线束一致，0参考平面为靶的底面(即光子出射面)；光子的飞行运动方向采用直角坐标系表达，坐标系方向平行靶的坐标系，光子采用单位动量方向在三个坐标方向上投影(px，py，pz)表示光子的方向，且有
Figure BDA0003109382900000102
pz＝1表示该光子的发射方向平行于射线束中心轴(垂直于治疗平面)，pz＝0表示该光子的发射方向在0平面上水平发射。
本发明实施例中的光子均按照能谱进行说明，每一条能谱线的能量宽度为BinE，每一条能谱线为一个能量组，将韧致辐射光子分成若干能量组，用Ei表示，i为自然数，表示第i组光子的能量分布中值为Ei，该组光子的能量分布宽度为
Figure BDA0003109382900000103
离散化能谱中，各个能量带的光子源模型函数形式统一，其函数中各个参量均与能量高度相关，是能量的函数。这种离散化能谱的源模型函数提供了光子的能量信息。
本发明的医用直线加速器的光子源模型函数包括所述原射线光子源模型函数和所述散射光子源模型函数，包括：
第一部分：构建原射线光子源模型函数
(1)、开启医用直线加速器，用产生的高能电子束轰击所述医用直线加速器的X线靶，将初始相空间平面上记录的靶面上产生的轫致辐射光子信息保存在第一相空间文件中，以及，将相空间平面上记录的经均整滤过器后的治疗光子束光子信息保存在第二相空间文件中。
为了获取正确的实验数据，需要在所述电子束流中心轴上适当位置选择一个与该中心轴垂直的平面，通过蒙特卡洛模拟计算获得经过该平面的所有光子信息，包括光子在该平面上的落点位置信息、单位动量分量、能量以及粒子类型，该平面为相空间平面，相空间平面的位置可根据实际应用需求选定，本发明不作具体的限定。
相空间平面安放在所述射线均整器的下方，紧靠均整器底端面的位置。在相空间平面上需要记录两类光子信息，一类是初始韧致辐射光子的光子信息，另一类是治疗射线束的光子信息。
用蒙特卡洛方法在初始相空间平面上记录初始韧致辐射光子的光子信息，只保留加速器的钨靶、初级准直器的部件，移除均整滤过器和二级准直器的部件；用蒙特卡洛方法在相空间平面上记录治疗射线束的光子信息，保留加速器的钨靶、靶衬、初级准直器、均整滤过器，移除加速器的二级准直器的部件描述。
需要说明的是，在初始相空间平面上记录的韧致辐射光子信息的相空间文件，本发明实施例记为第一相空间文件，在相空间平面上记录治疗射线束的光子信息的相空间文件，记为第二相空间文件。
(2)通过光子逆向飞行计算，将第一相空间文件和第二相空间文件转换为靶平面上的相空间文件，分别记为B-PhSp和R-PhSp；将相空间文件R-PhSp中记录的光子分离为原射线光子和散射线光子，得到两个新的相空间文件分别记为p-PhSp和s-PhSp。
为了将原射线光子源模型与散射光子源模型的位置统一到加速器的靶平面所在的平面上(即0参考平面)，必须将在上述两个相空间平面上记录的这两类光子(即韧致辐射光子和治疗光子束光子)经过光子逆向飞行计算处理，将它们分别转换到0参考平面上，相应的得到两个数据文件，即第一相空间文件转换为靶平面上的相空间文件，记为相空间文件B-PhSp，第二相空间文件转换为靶平面上的相空间文件，记为R-PhSp。根据在0参考平面上的韧致辐射光子数据文件可以得到韧致辐射光子在0参考平面上的注量分布。
在0参考平面上的治疗射线束的光子数据文件包含了原射线光子和散射线光子，需要做原、散射线分离处理，以便获得原射线光子束在0参考平面上的注量分布map矩阵。在相空间平面上所记录治疗射线束的光子，含有始发于加速器靶面的原射线光子以及均整滤过器中产生的散射线光子。原射线光子和散射线光子分离后，得到两个新的对应的相空间文件分别记为p-PhSp和s-PhSp。
在一个优选地实施例中，将相空间文件R-PhSp中记录的光子分离为原射线光子和散射线光子，具体可为：
相空间文件R-PhSp记录的光子中，将位置坐标位于初级准直器的圆锥形孔径在靶底端面上的几何投影区域之内的光子记为原射线光子，所述原射线光子可沿着初级准直器孔径回到所述初始相空间平面上；将分布在所述几何投影区域之外的光子记为散射线光子，而散射线光子的始发位置不在靶面上，而是分布在靶所在平面上的一个较大范围内，且散射线光子不能通过初级准直器孔径回到原来的相空间平面上。
(3)构建轫致辐射光子注量分布函数。
需要说明的是，在0参考面上任意一点的初始光子数量与击打该靶点的电子数量成正比，而真实的电子束流是具有确定半高宽的二维高斯分布。轰击靶平面的电子束流在靶入射面上的注量分布为二维高斯分布，该高斯函数中电子分布的全峰半高宽FWHM是已知的，决定了加速器的光子源尺寸，可来源于加速器厂商发布信息，或来自文献或本发明中推荐的典型值，电子束流的分布标准差δ＝FWHM/2.335；其能量也是二维高斯分布，决定了韧致辐射光子的能谱以及角分布特性。靶平面上任意一点产生的初始韧致辐射光子的产额正比于入射该点的电子数量(即：注量)，所以，在x线靶底端面发射的初始韧致辐射光子的产额为高斯分布，此高斯函数半高宽等于打靶电子束流的注量分布半高宽。
因此，加速器的靶产生的韧致辐射光子的注量分布的标准差δ和分布函数形式是已知的，为二维高斯函数，光子的总计数由模拟计算所需的光子总数决定，如下式所示，靶平面上任意一点产生的初始韧致辐射光子的注量分布函数可表示为：
Figure BDA0003109382900000121
且
Figure BDA0003109382900000122
其中，σ为分布标准差；(xs，ys)为靶平面上某一个点的坐标，下标s表示该坐标在靶平面上；NB0为靶面中心点的光子注量(产额)，NB(xs，ys)为靶平面上任一点(xs，ys)处的光子注量，下标B是轫致辐射光子，下标B0表示为韧致辐射光子在靶中心点的光子计数。
则离散化能谱后，加速器x线靶的靶面产生的韧致辐射光子注量分布函数为：
Figure BDA0003109382900000123
以及，
Figure BDA0003109382900000124
其中，
Figure BDA0003109382900000125
为
Figure BDA0003109382900000126
在
Figure BDA0003109382900000127
的分支比；
Figure BDA0003109382900000128
表示能量为Ei的韧致辐射光子总计数，NB表示所有能量的韧致辐射光子的总计数，NB0表示为靶中心点产生的韧致辐射光子计数；下标B是韧致辐射标识，Ei是第i条谱线的能量标识。
(3)根据韧致辐射光子在动量分量方向的注量分布，并拟合得到韧致辐射光子在动量分量方向的注量分布函数。
在一个优选地实施例中，将所述相空间文件B-PhSp中所记录的光子的位置坐标参量置于靶平面中心点(0，0，0)，其中靶平面设为0参考平面。
针对不同的能量，将X线靶上总的光子分成若干个能量组，例如分成20个能量组，各个能量组产生的韧致辐射光子在0参考面上按照(x,y,)的二维分布剖面图如图4所示，可知，在靶面上任意一点的电子注量的能谱相同，因此，在靶平面上任意一点产生的初始韧致辐射光子的能谱以及角分布特性相同。因此，将靶平面上产生的全部韧致辐射光子的位置坐标归于靶平面中心，相当于打靶电子束流增加n倍后在靶中心点产生的韧致辐射光子数量，而又不改变该靶面中心点处的初始韧致辐射光子的能谱和角分布特性，等价于加速器出射窗引出的电子束全部轰击靶中心产生的效果，可以有效地提高靶平面上光子数量的记录效率，从而提高工作效率。
因此，将靶面产生的所有韧致辐射光子的位置至于靶面中心，做靶面中心点的韧致辐射光子的沿动量的分量方向的注量分布并做函数拟合，进而得到韧致辐射光子沿动量的分量方向上的光子概率分布及函数拟合，该概率分布是指在该点产生的韧致辐射光子中沿某个方向发射的概率，因此，在靶面上任一点产生的韧致辐射光子的沿动量的分量方向上的光子概率分布函数也相同。
因此，获得韧致辐射光子沿动量的分量方向的注量分布，如图4所示，并拟合函数为：
Figure BDA0003109382900000131
其中，
Figure BDA0003109382900000132
表示在靶中心点产生的第i能量组的韧致辐射光子沿动量的分量方向的注量分布；
Figure BDA0003109382900000133
和
Figure BDA0003109382900000134
为韧致辐射光子中能量为Ei的光子沿动量分量方向上的注量分布标准差，
Figure BDA0003109382900000135
和
Figure BDA0003109382900000136
是函数表达式的系数，等于0度角方向上发射的韧致辐射光子计数，韧致辐射光子飞行的动量方向用单位动量在三维直角坐标系上的投影值表示，即(px，py，pz)，并且
Figure BDA0003109382900000137
其中，B表示韧致辐射光子，Ei表示第i条谱线能量，m表示动量(momentum)，即，
Figure BDA0003109382900000138
和
Figure BDA0003109382900000139
为动量分布的两个标准差。
(4)根据靶面上注量点发射韧致辐射光子在动量分量上的注量分布函数得到注量概率分布函数。
因为靶面上各个注量点产生的韧致辐射光子的角注量概率一样，则将所述靶面上注量点发射韧致辐射光子在动量分量上的注量分布函数进行拟合，得到靶面上注量点发射韧致辐射光子在动量分量上的注量概率分布函数为：
Figure BDA0003109382900000141
其中，
Figure BDA0003109382900000142
和
Figure BDA0003109382900000143
是韧致辐射光子沿动量分量方向的注量概率分布函数表达式的系数，根据
Figure BDA0003109382900000144
和
Figure BDA0003109382900000145
的归一计算得到。
(5)根据轫致辐射光子沿动量的分量方向的注量分布与原射线光子沿动量的分量方向的注量分布之差，得到轫致辐射光子在动量分量方向上被所述均整滤过器吸收的注量分布，并拟合得到轫致辐射光子在动量分量方向上被所述均整滤过器吸收的注量分布函数。
将所述轫致辐射光子的相空间文件B-PhSp及原射线光子的相空间文件中所记录的光子位置坐标置于靶中心(0，0，0)，并分别得到轫致辐射光沿动量的分量方向的注量分布和原射线光子沿动量的分量方向的注量分布B-map(如图4所示)和p-map，轫致辐射光子沿动量的分量方向的注量分布与原射线光子沿动量的分量方向的注量分布之差，即为在动量分量方向上被射线均整器吸收的韧致辐射光子的注量分布，如图5所示。
对所述轫致辐射光子在动量分量方向上被所述均整滤过器吸收的注量分布进行函数拟合，得到轫致辐射光子在动量分量方向上被所述均整滤过器吸收的注量分布函数：
Figure BDA0003109382900000146
其中，
Figure BDA0003109382900000147
表示在动量的分量方向(px，py)上第i条谱线能量的韧致辐射光子中被均整滤过器吸收的注量分布；
Figure BDA0003109382900000148
和
Figure BDA0003109382900000149
分别为上述两个二维高斯函数的标准差；
Figure BDA00031093829000001410
和
Figure BDA00031093829000001411
分别是上述两个高斯函数表达式的系数。
(6)将任一注量点处第i条谱线能量的韧致辐射光子总计数做数学归一处理，得到靶面上任一注量点发射韧致辐射光子在动量分量方向上被所述均整滤过器吸收的注量概率分布函数：
Figure BDA0003109382900000151
其中，
Figure BDA0003109382900000152
是韧致辐射光子被均整滤过器吸收的概率；
Figure BDA0003109382900000153
和
Figure BDA0003109382900000154
分别是该概率分布模型函数的两个系数，根据
Figure BDA0003109382900000155
和
Figure BDA0003109382900000156
的数学归一计算获得。
(7)构建原射线光子源模型函数。
需要说明的是，上述加速器靶发射的韧致辐射光子注量分布函数提供了光子发射点的坐标位置信息(xs，ys，zs)和产额信息NB(xs，ys)；靶面上注量点发射韧致辐射光子在动量分量上的注量概率分布函数与靶面上注量点发射韧致辐射光子在动量分量方向上被加速器均整器吸收的注量概率分布函数之差为该注量点出射的原射线束光子在动量分量上的注量概率，它提供了光子的动量方向信息(px，py，pz)和出现光子的概率。
根据所述靶面上注量点发射韧致辐射光子在动量分量上的注量概率分布函数与靶面上注量点发射韧致辐射光子在动量分量方向上被所述均整滤过器吸收的注量概率分布函数之差，得到靶面上注量点出射的原射线光子在动量分量上的注量概率分布函数，因此，再结合所述轫致辐射光子注量分布函数，得到原射线光子源模型函数为：
Figure BDA0003109382900000157
其中：
Figure BDA0003109382900000158
以及
Figure BDA0003109382900000159
以及
Figure BDA00031093829000001510
其中，原射线光子源模型位于0参考平面上，且原射线光子飞行的动量方向用单位动量在三维直角坐标系上的投影值表示，即(px，py，pz)，并且
Figure BDA00031093829000001511
因此，略去zs和pz。
上述函数式中各个参量均由拟合函数的过程中给出，
Figure BDA0003109382900000161
表示韧致辐射光子在动量分量上的注量概率分布；
Figure BDA0003109382900000162
是韧致辐射光子被均整滤过器吸收的概率；
Figure BDA0003109382900000163
为韧致辐射光子在动量分量上的注量概率分布函数的系数，
Figure BDA0003109382900000164
为韧致辐射光子在动量分量上的注量概率分布函数的分布标准差；
Figure BDA0003109382900000165
是韧致辐射光子在动量分量上的注量概率函数的系数，
Figure BDA0003109382900000166
是韧致辐射光子被均整滤过器吸收的概率函数的标准差，(f-absorbed)是被均整滤过器吸收的标识。
需要说明的是，本发明医用直线加速器治疗束中原射线光子源模型是关于光子信息参量xs、ys、zs、px、py、pz、Ei的函数；其中，在靶平面上产上的韧致辐射光子注量分布的标准差σ已知、模型函数中的参量
Figure BDA0003109382900000167
由拟合函数中给出。由于加速器靶平面的位置高度是已知固定值，因此，此处可略去zs的书写；另外，光子的飞行运动方向采用直角坐标系表达，坐标系方向平行靶的坐标系，光子采用单位动量方向在三个坐标方向上投影(px，py，pz)表示光子的方向，有
Figure BDA0003109382900000168
pz＝1表示该光子的发射方向平行于射线束中心轴(垂直于治疗平面)，pz＝0表示该光子的发射方向在0平面上水平发射，因此，pz的书写也可略去。
对于加速器治疗光子束中的原射线光子，xs、ys定义域满足：
Figure BDA0003109382900000169
即光子源的几何尺寸，在该模型中所描述原射线光子的发射点位置在源上。
第二部分：构建散射光子源模型函数
完整的散射线的光子源模型函数含两个组成部分，即在0参考平面上散射光子注量分布函数和某一注量点发射的散射光子在动量分量上的光子注量概率分布函数，共同构建散射光子源模型函数。
具体包括以下步骤：
(1)根据散射光子在0参考平面的注量分布拟合，得到散射光子的注量分布函数散射光子的注量分布函数。
需要说明的是，将散射光子的相空间文件s-PhSp中所记录的散射光子的分布区域进行离散化处理，形成像素格，并以像素格为单位获取0参考平面上的散射光子的注量分布，用双高斯函数多项式对所述散射光子的注量分布进行拟合，得到散射光子的注量分布函数为：
Figure BDA00031093829000001610
其中，(xs，ys)为0参考平面上的任一点位置坐标，
Figure BDA0003109382900000171
是在注量点(xs，ys)的散射光子数，ai、bi是能量为Ei的散射光子的注量分布的拟合函数的系数；
Figure BDA0003109382900000172
是能量为Ei的散射光子的注量分布的拟合函数的分布标准差；参量下标中的s表示散射线标识。
(2)根据靶面上各注量点发射的散射光子在动量分量上的注量分布函数得到散射光子在动量分量上的平均注量概率分布函数。
由于在0平面上的各个注量点的散射光子在动量分量方向上的注量概率分布存在轻度各向异性，因此，本发明实施例中以若干注量点的注量概率分布的平均作为各个注量点的散射光子在动量分量方向上的注量概率的代表，这种近似方式不影响模型的准确性和精度。
因此，靶面上各注量点发射的散射光子在动量分量上的注量分布的拟合函数为：
Figure BDA0003109382900000173
上式中，参量a'i、b'i、c'i、d'i为动量分量上的注量分布的拟合函数的系数；
Figure BDA0003109382900000174
为动量分量上的注量分布的拟合函数的标准差。各个参量的下标中，i表示第i条谱线能量组标识。
则，根据散射光子在动量分量上的注量分布函数得到散射光子在动量分量上的平均注量概率分布函数为：
Figure BDA0003109382900000175
其中，函数定义域r可用下面数学方法给出：
Figure BDA0003109382900000176
散射光子沿动量分量方向的平均注量概率分布函数中的参量ai、bi、ci、di、
Figure BDA0003109382900000177
Figure BDA0003109382900000178
由拟合过程给出，其中ai、bi、ci、di为函数的系数，
Figure BDA0003109382900000179
为沿动量分布的分布标准差，这些参量其与能量密切相关，可以对能量拟合，分别形成各个参量与能量的关系函数。参量小标的小写字母s是散射光子标识，其它下标意义同上。
(3)根据散射光子的注量分布函数和靶所在平面上各注量点发射的散射光子在动量分量上的平均注量概率分布函数构建散射光子源模型函数：
Figure BDA0003109382900000181
上述散射光子注量分布函数提供了散射光子源中光子的初始坐标(xs，ys，zs)和该注量点的散射光子计数Nscatter(xs，ys)；散射光子在动量分量上的注量概率分布函数提供了该注量点发射的散射光子的动量方向信息(px，py，pz)和概率。
以上对本发明实施例进行了详细介绍，本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处，综上所述，本说明书内容不应理解为对本发明的限制。

Claims (6)

1.一种医用直线加速器光子源模型函数的构建方法，其特征在于，医用直线加速器光子源模型包括原射线光子源模型和散射线光子源模型两部分，该方法包括：
开启医用直线加速器，用产生的高能电子束轰击所述医用直线加速器的X线靶，将初始相空间平面上记录的靶面上产生的轫致辐射光子信息保存在第一相空间文件中，以及，将相空间平面上记录的经均整滤过器后的治疗光子束光子信息保存在第二相空间文件中；
通过光子逆向飞行计算，将第一相空间文件和第二相空间文件转换为靶平面上的相空间文件，分别记为B-PhSp和R-PhSp；将相空间文件R-PhSp中记录的光子分离为原射线光子和散射线光子，得到两个新的相空间文件分别记为p-PhSp和s-PhSp；
根据所述轫致辐射光子注量分布函数、靶面上注量点发射韧致辐射光子在动量分量上的注量概率分布函数和靶面上注量点发射韧致辐射光子在动量分量方向上被所述均整滤过器吸收的注量概率分布函数构建原射线光子源模型函数；
根据所述散射线光子注量分布函数和靶所在的平面上各注量点发射的散射光子在动量分量上的平均注量概率分布函数构建散射线光子源模型函数。
2.如权利要求1所述的一种医用直线加速器光子源模型函数的构建方法，其特征在于，将相空间文件R-PhSp中记录的光子分离为原射线光子和散射线光子，具体包括：相空间文件R-PhSp记录的光子中，将位置坐标位于初级准直器的圆锥形孔径在靶底端面上的几何投影区域之内的光子记为原射线光子，所述原射线光子可沿着初级准直器孔径回到所述初始相空间平面上；将分布在所述几何投影区域之外的光子记为散射线光子。
3.如权利要求1或2所述的一种医用直线加速器光子源模型函数的构建方法，其特征在于，将所述相空间文件B-PhSp和相空间文件p-PhSp中所记录的光子的位置坐标参量置于靶平面中心点(0，0，0)，其中靶平面设为0参考平面。
4.如权利要求3所述的一种医用直线加速器光子源模型函数的构建方法，其特征在于，用能量宽度BinE将连续能谱的韧致辐射光子分成若干能量组，用Ei表示，i为自然数，表示第i组光子的能量分布中值，该组光子的能量分布宽度为
Figure FDA0003109382890000011
5.如权利要求4所述的一种医用直线加速器光子源模型函数的构建方法，其特征在于，所述构建原射线光子源模型函数的步骤，具体包括：
所述轫致辐射光子注量分布函数为：
Figure FDA0003109382890000021
以及，
Figure FDA0003109382890000022
其中，
Figure FDA0003109382890000023
为
Figure FDA0003109382890000024
在NB(xs，ys)的分支比；
Figure FDA0003109382890000025
表示能量为Ei的韧致辐射光子总计数，NB表示所有能量的韧致辐射光子的总计数；NB0表示为靶中心点产生的韧致辐射光子计数；σ为靶面产生初始韧致辐射光子分布的标准差，即，击靶的电子束流的分布标准差，以及，击靶电子束流的半高宽FWHM＝2.335*σ；
根据韧致辐射光子在动量分量方向的注量分布，并拟合得到韧致辐射光子在动量分量方向的注量分布函数为：
Figure FDA0003109382890000026
其中，
Figure FDA0003109382890000027
表示在靶中心点产生的第i能量组的韧致辐射光子沿动量的分量方向的注量分布；
Figure FDA0003109382890000028
和
Figure FDA0003109382890000029
为韧致辐射光子中能量为Ei的光子沿动量分量方向上的注量分布标准差，
Figure FDA00031093828900000210
和
Figure FDA00031093828900000211
是函数表达式的系数，等于0度角方向上发射的韧致辐射光子计数，韧致辐射光子飞行的动量方向用单位动量在三维直角坐标系上的投影值表示，即(px，py，pz)，并且
Figure FDA00031093828900000212
B表示韧致辐射光子，Ei表示第i条谱线能量，m表示动量(momentum)，即，
Figure FDA00031093828900000213
和
Figure FDA00031093828900000214
为动量分布的两个标准差；
因为靶面上个点的电子束流注量的能谱相同，所以靶面上个点产生的韧致辐射光子的动量分量上的注量概率分布与靶中心点的相同，则，靶面上注量点发射韧致辐射光子在动量分量上的注量概率分布函数为：
Figure FDA00031093828900000215
其中，
Figure FDA00031093828900000216
和
Figure FDA00031093828900000217
是韧致辐射光子沿动量分量方向的注量概率分布函数表达式的系数，根据
Figure FDA00031093828900000218
和
Figure FDA00031093828900000219
的归一计算得到；
Figure FDA00031093828900000220
表示韧致辐射光子在动量分量上的注量概率分布；
将所述相空间文件B-PhSp及p-PhSp所记录的光子沿着光子单位动量分量做光子通量分布，B-PhSp光子与p-PhSp光子的通量分布差即为韧致辐射光子被加速器均整滤过器吸收的光子在沿单位动量分量方向上的光子分布，记为
Figure FDA0003109382890000031
则根据轫致辐射光子沿动量的分量方向的注量分布与原射线光子沿动量的分量方向的注量分布之差，得到轫致辐射光子在动量分量方向上被所述均整滤过器吸收的注量分布，并拟合得到轫致辐射光子在动量分量方向上被所述均整滤过器吸收的注量分布函数：
Figure FDA0003109382890000032
其中，
Figure FDA0003109382890000033
表示在动量的分量方向(px，py)上第i条谱线能量的韧致辐射光子中被均整滤过器吸收的注量分布；
Figure FDA0003109382890000034
和
Figure FDA0003109382890000035
分别为上述两个二维高斯函数的标准差；
Figure FDA0003109382890000036
和
Figure FDA0003109382890000037
分别是上述两个高斯函数表达式的系数；
将任一注量点处第i条谱线能量的韧致辐射光子总计数做数学归一处理，得到靶面上注量点发射韧致辐射光子在动量分量方向上被所述均整滤过器吸收的注量概率分布函数为：
Figure FDA0003109382890000038
其中，
Figure FDA0003109382890000039
是韧致辐射光子被均整滤过器吸收的概率；
Figure FDA00031093828900000310
和
Figure FDA00031093828900000311
分别是该概率分布模型函数的两个系数，根据
Figure FDA00031093828900000312
和
Figure FDA00031093828900000313
的数学归一计算获得；
根据所述靶面上注量点发射的韧致辐射光子在动量分量上的注量概率分布函数与靶面上注量点发射的韧致辐射光子在动量分量方向上被所述均整滤过器吸收的注量概率分布函数之差，得到靶面上注量点出射的原射线光子在动量分量上的注量概率分布函数，因此，再结合所述轫致辐射光子注量分布函数，得到原射线光子源模型函数为：
Figure FDA00031093828900000314
其中：
Figure FDA00031093828900000315
以及
Figure FDA0003109382890000041
以及
Figure FDA0003109382890000042
其中，原射线光子源模型位于0参考平面上，且原射线光子飞行的动量方向用单位动量在三维直角坐标系上的投影值表示，即(px，py，pz)，并且
Figure FDA0003109382890000043
因此，略去zs和pz；上述函数式中各个参量均由拟合函数的过程中给出，
Figure FDA0003109382890000044
表示韧致辐射光子在动量分量上的注量概率分布；
Figure FDA0003109382890000045
是韧致辐射光子被均整滤过器吸收的概率；
Figure FDA0003109382890000046
为韧致辐射光子在动量分量上的注量概率分布函数的系数，
Figure FDA0003109382890000047
为韧致辐射光子在动量分量上的注量概率分布函数的分布标准差；
Figure FDA0003109382890000048
是韧致辐射光子在动量分量上的注量概率函数的系数，
Figure FDA0003109382890000049
是韧致辐射光子被均整滤过器吸收的概率函数的标准差，(f-absorbed)是被均整滤过器吸收的标识。
6.如权利要求4或5所述的一种医用直线加速器光子源模型函数的构建方法，其特征在于，所述构建散射线光子源模型函数的步骤，具体包括：
将所述散射光子的相空间文件s-PhSp中所记录的散射光子的分布区域进行离散化处理，形成像素格，并以像素格为单位获取0参考平面上的散射光子的注量分布，用双高斯函数多项式对所述散射光子的注量分布进行拟合，得到散射光子的注量分布函数为：
Figure FDA00031093828900000410
其中，(xs，ys)为0参考平面上的任一点位置坐标，
Figure FDA00031093828900000411
是在注量点(xs，ys)的散射光子数，ai、bi是能量为Ei的散射光子的注量分布的拟合函数的系数；
Figure FDA00031093828900000412
是能量为Ei的散射光子的注量分布的拟合函数的分布标准差；参量下标中的s表示散射线标识；
靶所在的平面上各注量点发射的散射光子在动量分量上的注量分布的拟合函数为：
Figure FDA00031093828900000413
Figure FDA0003109382890000051
上式中，参量a'i、b'i、c'i、d'i为动量分量上的注量分布的拟合函数的系数；
Figure FDA0003109382890000052
为动量分量上的注量分布的拟合函数的标准差；各个参量的下标中，i表示第i条谱线能量组标识；
则，根据散射光子在动量分量上的注量分布函数得到散射光子在动量分量上的平均注量概率分布函数为：
Figure FDA0003109382890000053
其中，函数定义域r可用下面数学方法给出：
Figure FDA0003109382890000054
散射光子沿动量分量方向的平均注量概率分布函数中的参量ai、bi、ci、di、
Figure FDA0003109382890000055
由拟合过程给出，其中ai、bi、ci、di为函数的系数，
Figure FDA0003109382890000056
为沿动量分布的分布标准差，这些参量与能量密切相关，可以对能量拟合，分别形成各个参量与能量的关系函数；
根据散射光子的注量分布函数和靶所在的平面上各注量点发射的散射光子在动量分量上的平均注量概率分布函数构建散射光子源模型函数为：
Figure FDA0003109382890000057
CN202110645477.5A 2021-06-10 2021-06-10 医用直线加速器光子源模型函数的构建方法 Active CN113230549B (zh)

Priority Applications (2)

Application Number Priority Date Filing Date Title
CN202110645477.5A CN113230549B (zh) 2021-06-10 2021-06-10 医用直线加速器光子源模型函数的构建方法
US17/828,070 US20220401754A1 (en) 2021-06-10 2022-05-31 Method for constructing photon source model function of medical linear accelerator

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110645477.5A CN113230549B (zh) 2021-06-10 2021-06-10 医用直线加速器光子源模型函数的构建方法

Publications (2)

Publication Number Publication Date
CN113230549A true CN113230549A (zh) 2021-08-10
CN113230549B CN113230549B (zh) 2021-12-31

Family

ID=77139443

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110645477.5A Active CN113230549B (zh) 2021-06-10 2021-06-10 医用直线加速器光子源模型函数的构建方法

Country Status (2)

Country Link
US (1) US20220401754A1 (zh)
CN (1) CN113230549B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117523138A (zh) * 2023-11-06 2024-02-06 南京航空航天大学 构造肿瘤拟合曲面的方法、电阻抗成像逆问题求解方法、检测系统

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103091693A (zh) * 2013-01-14 2013-05-08 中国科学院合肥物质科学研究院 一种基于二维平面剂量判断医用直线加速器稳定性的方法
CN104043203A (zh) * 2014-06-12 2014-09-17 合肥工业大学 基于医用直线加速器照射束特征建立照射源模型的方法
CN105204056A (zh) * 2015-09-29 2015-12-30 广州瑞多思医疗科技有限公司 放射治疗中在线剂量监测和验证的方法
CN106215331A (zh) * 2016-07-06 2016-12-14 中国科学院合肥物质科学研究院 一种基于x射线成像平板探测器的快速剂量验证方法
CN110237445A (zh) * 2019-07-05 2019-09-17 北京理工大学 基于epid的在体三维剂量监测及验证方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103091693A (zh) * 2013-01-14 2013-05-08 中国科学院合肥物质科学研究院 一种基于二维平面剂量判断医用直线加速器稳定性的方法
CN104043203A (zh) * 2014-06-12 2014-09-17 合肥工业大学 基于医用直线加速器照射束特征建立照射源模型的方法
CN105204056A (zh) * 2015-09-29 2015-12-30 广州瑞多思医疗科技有限公司 放射治疗中在线剂量监测和验证的方法
CN106215331A (zh) * 2016-07-06 2016-12-14 中国科学院合肥物质科学研究院 一种基于x射线成像平板探测器的快速剂量验证方法
CN110237445A (zh) * 2019-07-05 2019-09-17 北京理工大学 基于epid的在体三维剂量监测及验证方法

Also Published As

Publication number Publication date
CN113230549B (zh) 2021-12-31
US20220401754A1 (en) 2022-12-22

Similar Documents

Publication Publication Date Title
US10576303B2 (en) Methods and systems for beam intensity-modulation to facilitate rapid radiation therapies
Vassiliev et al. Validation of a new grid-based Boltzmann equation solver for dose calculation in radiotherapy with photon beams
Lovelock et al. A Monte Carlo model of photon beams used in radiation therapy
Schach von Wittenau et al. Correlated histogram representation of Monte Carlo derived medical accelerator photon‐output phase space
US9583302B2 (en) Convergent photon and electron beam generator device
Sikora et al. A virtual photon source model of an Elekta linear accelerator with integrated mini MLC for Monte Carlo based IMRT dose calculation
Sikora et al. A virtual source model of electron contamination of a therapeutic photon beam
Dipuglia et al. Validation of a Monte Carlo simulation for microbeam radiation therapy on the imaging and medical beamline at the Australian Synchrotron
EP3302699A1 (en) Method of selecting beam geometries
Naqvi et al. Convolution/superposition using the Monte Carlo method
CN113230549B (zh) 医用直线加速器光子源模型函数的构建方法
Bencheikh et al. Dosimetry investigation and evaluation for removing flattening filter configuration of linac: Monte Carlo study
US8280003B2 (en) Method for calculating head scatter phase space for radiation treatment using a multi-leaf collimator with dynamic jaws
Sikora Virtual source modelling of photon beams for Monte Carlo based radiation therapy treatment planning
Keshmiri et al. A high‐resolution dose calculation engine for X‐ray microbeams radiation therapy
CN116920285A (zh) 辐射治疗计划中的监测单元优化约束
Chiavassa et al. Validation of the analytical irradiator model and Monte Carlo dose engine in the small animal irradiation treatment planning system µ-RayStation 8B
US9227083B2 (en) Device and method for radiation dosimetry
Grzanka Modelling beam transport and biological effectiveness to develop treatment planning for ion beam radiotherapy
De Simoni et al. FRED: a fast Monte Carlo code on GPU for quality control in Particle Therapy
Chaves et al. Basic dosimetry of radiosurgery narrow beams using Monte Carlo simulations: A detailed study of depth of maximum dose
Harris Development and validation of an electron Monte Carlo model for the Elekta Infinity accelerator
Paganetti Monte Carlo simulations
US11235174B2 (en) Minibeam radiotherapy device
Fix et al. Treatment planning aspects and Monte Carlo methods in proton therapy

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant