固体火箭发动机射线检测仿真系统①

任晓双，王 冬，陈 平

(1.中北大学 信息探测与处理山西省重点实验室，太原 030051；2.中国人民解放军96630部队，北京 102206)

0 引言

射线检测中的仿真技术是计算机图形学、数学、无损检测和软件设计等方面的综合应用。该技术是对实际射线检测系统及检测过程中的原理进行建模，并利用计算机软件模拟计算来获取与实际检测相同的结果[1]。利用计算机CAD辅助软件对固体火箭发动机进行建模，设计射线检测工艺方案，确定相应照射角度及相关参数。这样可缩短检测工艺的编制周期，节省大量人力、物力资源，且能对当前射线检测工艺进行优化和分析；依据仿真技术获取的理论透照图像，实现复杂结构的检测，进而排除结构因素对检测人员评片的影响。

目前，国内外研究人员对射线检测仿真软件进行了研究。美国Iowa州立大学无损评估中心的XRSIM仿真系统支持使用者自定义射线源、检测对象和胶片等的相关参数[2]，该系统获得了美国国家工业标准技术局基金项目的支持。BAM是德国的一个仿真系统[3]，该系统在起初的版本中只支持导入体素描述的检测对象文件，目前新版本支持CAD软件建立的各种格式的实体模型。欧洲各国联合开发了RADICAD系统，该系统可导入BRL-CAD软件包设计的CAD模型[4]。CIVA是法国原子能委员会开发的一套软件[5]，最初用于多种检测(超声、涡流和射线检测)数据的识别与处理，现在可实现这些检测系统的仿真计算。国内各大高校和研究机构也研究了射线仿真技术，但这些单位重点研究透射成像仿真算法，没有对仿真系统的物理特性深入探讨，目前没有成熟的射线检测计算机仿真软件。

为解决在固体火箭发动机射线检测中，缺少先验知识、工艺参数难确定、耗时长、耗费大等问题，本文从软件工程角度出发，采用面向对象的模块化设计模式，通过可视化的软件编程框架研究开发了固体火箭发动机射线检测仿真系统，同时针对固体火箭发动机模型尺寸大和结构复杂特点，进行仿真计算耗时长的问题，提出基于三维空间投影约束的求交加速算法，并利用GPU并行技术实现硬件加速[6]。

1 仿真系统总体设计

1.1 理论基础

X射线衰减定律和光线追踪技术[7]是本固体火箭发动机射线检测仿真系统的理论基础。假设理想状态下的点射线源发射出一系列锥形状的射线，击穿发动机打在背后的感光胶片上。在该过程中射线与发动机物质发生相互作用而产生衰减，经过衰减后的X光子数目N(E)遵循比尔定律[8]：

(1)

式中N0(E)为在每单位空间立体角上射线源发出含有总能量为E的光子总数目；ΔΩ为射线源和探测器探元相应的立体角度；u(x，y，z，E)为E能量下的发动机空间坐标(x，y，z)处的X射线线性衰减系数；l为X射线穿过发动机时的衰减路径。

从理论上分析，X射线检测的物理系统仿真主要包括X射线源产生的多谱仿真、待检测对象模体仿真和胶片成像仿真。投影成像算法利用光线追踪法对待检测模型求交计算得到衰减距离，然后根据式(1)计算出经过衰减后胶片上的能量值。

1.2 基本框架

本文的仿真系统利用了模块化的设计思想，依据实际射线检测系统各部分的物理原理来设计相应模块。仿真系统的基本框架如图1所示。整个仿真系统由射线源仿真模块、检测对象仿真模块、胶片特性仿真模块、投影成像算法模块以及信息交互模块构成。前4个模块主要是完成高能X射线的透照仿真，其中包括高能射线特性、透照参数、透照工艺、仿真模体、胶片黑度仿真等，这4个模块之间通过信息交互模块来传递信息。

图1 仿真系统的基本框架Fig.1 The basic framework of the simulation system

2 仿真系统物理建模

2.1 射线源仿真

X射线源的仿真表征主要由几何信息和能谱信息组成，本文采用的是9 MeV直线加速器。直线加速器在理想情况下是一个点光源，但实际上具有一定的尺寸和形状，焦点的形状取决于直线加速器的靶材设计，依据实际检测系统中2 mm焦点尺寸的加速器，采用圆形离散化方法表示加速器焦点，如图2所示，在2 mm直径的圆内分布5个点源，其中1个点源位于圆心位置，其余4个点源呈90°等角度布于圆周上。

图2 直线加速器示意图Fig.2 Linear accelerator schematic

X射线源的能量信息包括剂量值、能谱分布、方向剂量分布等参数。在很大程度上，射线源光谱的能量分布对仿真计算投影图像的质量有影响。其分布在数学形式上是连续的，但是现阶段几乎不使用模拟的数学表达式对其进行刻画。本文将这种连续分布的特性进行分阶段离散化处理，采用该方法在软件层面上有利于对其建立数学模型，另外可在编程方面使实现简单。对分段能量的X光子进行相应计算获取投影图像，最终融合各能量段的投影图像生成最终投影结果。

使用蒙特卡罗方法计算的X射线能谱如图3所示。在本文中将能谱(0～9 MeV)离散为16个能区，每个能谱区根据积分面积计算出该能量区间下对投影结果的贡献百分比，作为多能投影融合的能谱权值，材料在每个能谱区的衰减系数来自美国国家标准技术局(NIST)数据库。

图3 9 MeV直线加速器X射线能谱图Fig.3 9 MeV linear accelerator X-ray energy spectrum

2.2 固体火箭发动机模型仿真

固体火箭发动机模型应该由几何信息和材质信息两部分组成，其中几何信息包括发动机模型的物理位置及相应的尺寸等；材质信息包括发动机模型的材料密度、X射线线性衰减系数等。

固体火箭发动机是一个多零件、多材质、多类型的检测样本，本文采用SolidWorks软件[9]建立多材质模型STL样本来表示几何信息。图4为固体火箭发动机建立的简化三层结构模型，从左到右分别是外壳、绝热层、药柱。

图4 Solid Works软件设计的固体火箭发动机简化模型Fig.4 Solid rocket motor simplified model by Solid Works

STL样本采用三角面片来描述模型的表面结构，包括三角面片的顶点坐标和面法向量信息。同样的STL样本模型，采用三角面片数量的不同将导致模型的精确度不同，三角面片个数越多，模型就越精细，但是相应的计算时间也将增加。本文所有实验，简化三级固体火箭发动机模型采用103 682个三角面片描述其几何信息。

SolidWorks软件建立的多材质模型STL文件无法存入材料信息，本文采用如下方法来完成模型几何信息和材质信息的相互匹配：每导入一个固体火箭发动机STL配件由用户在材料信息表中选择材料信息和相关的衰减指数，与模型几何信息数据一同写入仿真系统，从而将表示多材质的固体火箭发动机STL模型装配到仿真系统中进行计算。

2.3 胶片特性仿真

本文采用胶片黑度值[10]是来表征胶片的特性仿真，因此确定仿真胶片的特性曲线是至关重要的。在实际检测中测定胶片特性曲线[11]，普遍使用的方法是，使用X射线源，固定管电压、使用滤波板，以保证辐射质量，然后再固定管电流、焦距、增感铅屏厚度等条件，仅改变曝光时间t这一参数，测定胶片的底片黑度D与曝光量E(表示管电流乘以曝光时间t)的关系，为显示斜率用D-lgE表示胶片特性曲线，该曲线近似为一条简单的上直斜线。在实际的固体火箭发动机射线检测中，射线源参数采用的是剂量值，而当管电压和滤波等其他条件固定时，吸收剂量值和曝光量成正比，因此认为胶片的黑度和吸收剂量也近似成正比，于是本文采用D-lgK特性曲线[12]来表征胶片特性仿真，其中K为达到黑度D时的吸收剂量，D-lgK曲线也近似为一条上斜直线，故计算仿真黑度值：

D=a×lg(b·K)+c

(2)

式中a为近似胶片特性直线斜率；b为矫正参数；c为截距。这3个参数由胶片特性所决定，可根据胶片特性曲线调整。

3 投影成像算法及加速

国内外在射线检测仿真系统上成像方法主要有两种：一种是使用体素模型进行光线追踪求交算法，但是当精度较高、模型体积大、数据量大时，仿真速度很慢，不符合固体火箭发动机模型的应用；另一种是使用表面模型进行光线求交算法，计算效率要比体素模型高，但是针对固体火箭发动机模型体积大数据量大的问题，求交速度依然是个瓶颈。本文重点研究基于STL表面模型的高效求交算法，提高射线与大数据量样本的求交速度。

3.1 成像算法描述

由X射线源发射出一系列锥形的X射线，穿透固体火箭发动机样本模型，撞击到固体火箭发动机后面的面阵探测器的各个探元上，该过程中X射线穿透固体火箭发动机并与其发生相互作用而衰减，这是投影成像的基本思想。

投影成像算法是仿真系统的核心，由以下三层循环组成：(1)射线源点阵循环，对每个离散的焦点计算投影；(2)固体火箭发动机模型的零件循环，通过X射线与每个零件相交计算，求出衰减距离；(3)像素矩阵循环，对探测器每个像素计算经过模型衰减后的X射线能量值。

在每次循环中首先要求出射线与固体火箭发动机模型三角形的交点，再对交点排序求出模型的衰减距离，其中关键的求交方法如下：

射线直线方程参数形式：

(3)

式中 (x0，y0，z0)为射线上一点；(m，n，p)为射线的方向向量。

三角形所在平面方程参数形式：

A(x-x1)+B(y-y1)+C(z-z1)=0

(4)

式中 (x1，y1，z1)为平面上一点(三角形任意一顶点)；(A，B，C)为三角平面的法向量。

先判断射线直线与三角平面是否平行，如果不平行将直线方程式(3)代入平面方程式(4)，求出交点，然后判断交点是否在三角形内。方法如下：

对于平面内任意一点如图5所示，都可由如下方程来表示：

P=A+u·(C-A)+v·(B-A)

(5)

图5 点与三角形的关系示意图Fig.5 Schematic of the relationship between points and triangles

如果让P位于三角形ABC内部(如图5所示)，u和v必须满足如下三个条件：

3.2 三维空间投影约束的求交加速算法

射线与STL模型求交要求遍历模型中的所有三角面片来完成直线与三角形求交运算，导致仿真计算耗时长的问题。针对这一问题，黄魁东等[13]采用将模型八叉树剖分的思想来加速射线与模型的求交运算，但是针对尺寸大、数据量大、结构复杂的固体火箭发动机模型，八叉树剖分模型需要耗很长时间，这对于仿真系统的实际应用需求无法满足。

本文对模型采用三维空间平面投影约束的方法快速滤除绝大多数与射线不相交的三角面片，从而加速遍历运算。将射线和三角面片顶点从三维空间中投影到任意两个二维平面上，这样形成了二维平面上3个点和线的关系，如图6所示。判断3个投影顶点是否在射线投影直线的同一侧，如果在同一侧说明射线与该三角面片肯定不相交，增加这个约束条件可避免射线与绝大多数不相交三角面片的求交计算，可节约大量的计算时间。

图6 二维投影平面上3个点和线的关系示意图Fig.6 The relationship between three points and line on

三维空间投影约束的求交加速算法如下：

针对模型所有三角面片，首先分别确定射线和三角面片3个顶点在XOY平面和XOZ平面的垂直投影直线方程和点坐标；然后判断这3个点是否在投影直线的同一侧，若都不在同一侧则计算出空间射线方程和三角面片的平面方程；其次根据计算出的射线方程和平面方程判断射线是否与平面相交，若相交则求出交点坐标；最后判断交点是否在三角面片内，若交点在三角面片内则表示射线与模型相交，记录相交点坐标信息用以后续计算最终的射线衰减距离。

对一条射线与模型所有三角面片遍历求交点并排序计算出衰减距离消耗的时间做实验，三维空间投影约束前后加速实验结果见表1，实验结果表明加约束条件后可以提高10倍左右的加速比。

表1 约束算法加速效果对比

3.3 CUDA并行加速

在射线投影计算过程中涉及到射线源不同离散焦点位置、探测器不同探元位置的射线与所有三角面片的求交运算，同时为了保证仿真投影结果的精度，要求仿真探测器分辨率大，计算量庞大。经过实验测试在CPU Intel Core i7-6700K(八核，主频4.00 GHz)的计算机上对模型所有三角面片遍历加约束条件后计算出衰减距离一次需要10 ms左右，生成4096×4096大小的投影，理论估计需要233.02 h，因此加约束条件后计算速度的提升依然无法满足实际应用。考虑到在投影计算过程中每条射线计算均是独立的，本文利用CUDA技术对核心算法进行加速，CUDA技术[14]是英伟达公司使用显卡GPU来辅助电脑CPU分担计算任务的技术，充分利用显卡的高性能并行计算能力，可提高算法的运算速度。CUDA线程并行是GPU细粒度的并行，线程并行运行的核函数都是相同的，CUDA流并行可实现对同一个核函数传递不同的参数，实现任务级别的并行。利用流并行机制可实现不同固体火箭发动机模型零件级别的并行，线程并行实现探测器探元级别上的并行计算。

对比CPU在不同GPU设备上做了计算速度实验，实验参数：射线源离散个数为5个；探测器大小为4096×4096，像素规格0.1 mm；壳体，27 840个三角面片；绝热层，38 400个三角面片；推进剂，37 442个三角面片；总计103 682个三角面片。实验结果见表2。

表2 CUDA加速效果对比

表2结果表明，在不同GPU硬件加速情况下可达到500～1000倍的加速比，GPU的性能越好，加速比将越大。基于CUDA 的算法加速，使得该方案能满足于实际固体火箭发动机射线检测工程的速度需求。

4 软件编程与实验结果

4.1 系统开发

根据固体火箭发动机射线检测仿真系统的功能分析和总体设计，系统主控平台采用微软公司MFC以文档为中心的程序架构思想，在Windows 7操作系统上，使用Visual studio 2013结合OpenGL三维可视化库[15]编写了固体火箭发动机射线检测仿真系统主界面，见图7。

图7 仿真系统主界面Fig.7 Simulation system main interface

4.2 实验结果

固体火箭发动机简化三层模型壳体、绝热层、药柱前肩位置的衰减距离可视化结果见图8。从仿真结果上看，对模型的物理结构仿真非常准确。

图8 壳体、绝热层、药柱的衰减距离可视化图Fig.8 Attenuation distance visualization of case,insulation and grain

根据实际试验测出0.02 Gy吸收剂量对应胶片2黑度，本文采取式(2)中a=1，c=0，计算得出b=5000。不同剂量下的黑度仿真可视化结果及黑度值分布曲线见图9。

在相同条件下对仿真胶片与实际胶片测量的黑度值作了对比实验。建立的简化固体火箭发动机仿真模型与真实发动机的尺寸等参数有所差别，定量误差分析要求黑度值的测量位置参数一致，因此本文采用定性分析黑度值仿真结果。参见图10。

(a)剂量为3 Gy

(b)剂量为10 Gy

(c)剂量为25 Gy

(a)固体火箭发动机实际胶片与投影仿真结果

(b)剂量为10 Gy在相同位置处

如图10所示，仿真黑度分布曲线趋势与实际胶片测得黑度曲线非常相似，在10 Gy剂量条件下仿真黑度数值大小与实际系统同剂量下胶片测得的黑度数值相近。在实际检测固体火箭发动机缺陷时更关注绝热层的黑度值，实验结果表明在绝热层的黑度仿真与实际测量曲线基本吻合，最大差值为0.23，与文献[10]的黑度值最大误差0.3相比，本文的黑度值仿真在检测范围内更准确，证明了本文黑度仿真模型的可靠性，绝热层黑度对比结果如图11所示。

图11 绝热层处仿真与实际黑度值对比Fig.11 Comparison of simulation and actual blackness values at the thermal insulation

5 结论

本文围绕固体火箭发动机的射线检测流程和原理，设计了非点光源与多能谱的固体火箭发动机射线检测仿真系统。在保证投影精度的前提下，通过加三维空间投影约束的STL模型求交算法和GPU并行技术实现了固体火箭发动机投影仿真与黑度仿真的工程应用。实验结果表明，该方案较之于CPU下STL模型求交算法有500～1000倍左右的加速比，仿真胶片与实际胶片黑度最大差值为0.23，与实际胶片测量的黑度更加接近，在可检测范围内更准确。该仿真系统为固体火箭发动机射线检测的工艺制定和优化搭建了一个稳定的仿真实验平台，具有很大的实用价值。