固體火箭發動機射線檢測仿真系統①

任曉雙，王 冬，陳 平

(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，與實際膠片測量的黑度更加接近，在可檢測范圍內更準確。該仿真系統為固體火箭發動機射線檢測的工藝制定和優化搭建了一個穩定的仿真實驗平臺，具有很大的實用價值。