固體火箭發(fā)動機射線檢測仿真系統(tǒng)①

任曉雙，王 冬，陳 平

(1.中北大學(xué) 信息探測與處理山西省重點實驗室，太原 030051；2.中國人民解放軍96630部隊，北京 102206)

0 引言

射線檢測中的仿真技術(shù)是計算機圖形學(xué)、數(shù)學(xué)、無損檢測和軟件設(shè)計等方面的綜合應(yīng)用。該技術(shù)是對實際射線檢測系統(tǒng)及檢測過程中的原理進行建模，并利用計算機軟件模擬計算來獲取與實際檢測相同的結(jié)果[1]。利用計算機CAD輔助軟件對固體火箭發(fā)動機進行建模，設(shè)計射線檢測工藝方案，確定相應(yīng)照射角度及相關(guān)參數(shù)。這樣可縮短檢測工藝的編制周期，節(jié)省大量人力、物力資源，且能對當(dāng)前射線檢測工藝進行優(yōu)化和分析；依據(jù)仿真技術(shù)獲取的理論透照圖像，實現(xiàn)復(fù)雜結(jié)構(gòu)的檢測，進而排除結(jié)構(gòu)因素對檢測人員評片的影響。

目前，國內(nèi)外研究人員對射線檢測仿真軟件進行了研究。美國Iowa州立大學(xué)無損評估中心的XRSIM仿真系統(tǒng)支持使用者自定義射線源、檢測對象和膠片等的相關(guān)參數(shù)[2]，該系統(tǒng)獲得了美國國家工業(yè)標(biāo)準(zhǔn)技術(shù)局基金項目的支持。BAM是德國的一個仿真系統(tǒng)[3]，該系統(tǒng)在起初的版本中只支持導(dǎo)入體素描述的檢測對象文件，目前新版本支持CAD軟件建立的各種格式的實體模型。歐洲各國聯(lián)合開發(fā)了RADICAD系統(tǒng)，該系統(tǒng)可導(dǎo)入BRL-CAD軟件包設(shè)計的CAD模型[4]。CIVA是法國原子能委員會開發(fā)的一套軟件[5]，最初用于多種檢測(超聲、渦流和射線檢測)數(shù)據(jù)的識別與處理，現(xiàn)在可實現(xiàn)這些檢測系統(tǒng)的仿真計算。國內(nèi)各大高校和研究機構(gòu)也研究了射線仿真技術(shù)，但這些單位重點研究透射成像仿真算法，沒有對仿真系統(tǒng)的物理特性深入探討，目前沒有成熟的射線檢測計算機仿真軟件。

為解決在固體火箭發(fā)動機射線檢測中，缺少先驗知識、工藝參數(shù)難確定、耗時長、耗費大等問題，本文從軟件工程角度出發(fā)，采用面向?qū)ο蟮哪K化設(shè)計模式，通過可視化的軟件編程框架研究開發(fā)了固體火箭發(fā)動機射線檢測仿真系統(tǒng)，同時針對固體火箭發(fā)動機模型尺寸大和結(jié)構(gòu)復(fù)雜特點，進行仿真計算耗時長的問題，提出基于三維空間投影約束的求交加速算法，并利用GPU并行技術(shù)實現(xiàn)硬件加速[6]。

1 仿真系統(tǒng)總體設(shè)計

1.1 理論基礎(chǔ)

X射線衰減定律和光線追蹤技術(shù)[7]是本固體火箭發(fā)動機射線檢測仿真系統(tǒng)的理論基礎(chǔ)。假設(shè)理想狀態(tài)下的點射線源發(fā)射出一系列錐形狀的射線，擊穿發(fā)動機打在背后的感光膠片上。在該過程中射線與發(fā)動機物質(zhì)發(fā)生相互作用而產(chǎn)生衰減，經(jīng)過衰減后的X光子數(shù)目N(E)遵循比爾定律[8]：

(1)

式中N0(E)為在每單位空間立體角上射線源發(fā)出含有總能量為E的光子總數(shù)目；ΔΩ為射線源和探測器探元相應(yīng)的立體角度；u(x，y，z，E)為E能量下的發(fā)動機空間坐標(biāo)(x，y，z)處的X射線線性衰減系數(shù)；l為X射線穿過發(fā)動機時的衰減路徑。

從理論上分析，X射線檢測的物理系統(tǒng)仿真主要包括X射線源產(chǎn)生的多譜仿真、待檢測對象模體仿真和膠片成像仿真。投影成像算法利用光線追蹤法對待檢測模型求交計算得到衰減距離，然后根據(jù)式(1)計算出經(jīng)過衰減后膠片上的能量值。

1.2 基本框架

本文的仿真系統(tǒng)利用了模塊化的設(shè)計思想，依據(jù)實際射線檢測系統(tǒng)各部分的物理原理來設(shè)計相應(yīng)模塊。仿真系統(tǒng)的基本框架如圖1所示。整個仿真系統(tǒng)由射線源仿真模塊、檢測對象仿真模塊、膠片特性仿真模塊、投影成像算法模塊以及信息交互模塊構(gòu)成。前4個模塊主要是完成高能X射線的透照仿真，其中包括高能射線特性、透照參數(shù)、透照工藝、仿真模體、膠片黑度仿真等，這4個模塊之間通過信息交互模塊來傳遞信息。

圖1 仿真系統(tǒng)的基本框架Fig.1 The basic framework of the simulation system

2 仿真系統(tǒng)物理建模

2.1 射線源仿真

X射線源的仿真表征主要由幾何信息和能譜信息組成，本文采用的是9 MeV直線加速器。直線加速器在理想情況下是一個點光源，但實際上具有一定的尺寸和形狀，焦點的形狀取決于直線加速器的靶材設(shè)計，依據(jù)實際檢測系統(tǒng)中2 mm焦點尺寸的加速器，采用圓形離散化方法表示加速器焦點，如圖2所示，在2 mm直徑的圓內(nèi)分布5個點源，其中1個點源位于圓心位置，其余4個點源呈90°等角度布于圓周上。

圖2 直線加速器示意圖Fig.2 Linear accelerator schematic

X射線源的能量信息包括劑量值、能譜分布、方向劑量分布等參數(shù)。在很大程度上，射線源光譜的能量分布對仿真計算投影圖像的質(zhì)量有影響。其分布在數(shù)學(xué)形式上是連續(xù)的，但是現(xiàn)階段幾乎不使用模擬的數(shù)學(xué)表達式對其進行刻畫。本文將這種連續(xù)分布的特性進行分階段離散化處理，采用該方法在軟件層面上有利于對其建立數(shù)學(xué)模型，另外可在編程方面使實現(xiàn)簡單。對分段能量的X光子進行相應(yīng)計算獲取投影圖像，最終融合各能量段的投影圖像生成最終投影結(jié)果。

使用蒙特卡羅方法計算的X射線能譜如圖3所示。在本文中將能譜(0～9 MeV)離散為16個能區(qū)，每個能譜區(qū)根據(jù)積分面積計算出該能量區(qū)間下對投影結(jié)果的貢獻百分比，作為多能投影融合的能譜權(quán)值，材料在每個能譜區(qū)的衰減系數(shù)來自美國國家標(biāo)準(zhǔn)技術(shù)局(NIST)數(shù)據(jù)庫。

圖3 9 MeV直線加速器X射線能譜圖Fig.3 9 MeV linear accelerator X-ray energy spectrum

2.2 固體火箭發(fā)動機模型仿真

固體火箭發(fā)動機模型應(yīng)該由幾何信息和材質(zhì)信息兩部分組成，其中幾何信息包括發(fā)動機模型的物理位置及相應(yīng)的尺寸等；材質(zhì)信息包括發(fā)動機模型的材料密度、X射線線性衰減系數(shù)等。

固體火箭發(fā)動機是一個多零件、多材質(zhì)、多類型的檢測樣本，本文采用SolidWorks軟件[9]建立多材質(zhì)模型STL樣本來表示幾何信息。圖4為固體火箭發(fā)動機建立的簡化三層結(jié)構(gòu)模型，從左到右分別是外殼、絕熱層、藥柱。

圖4 Solid Works軟件設(shè)計的固體火箭發(fā)動機簡化模型Fig.4 Solid rocket motor simplified model by Solid Works

STL樣本采用三角面片來描述模型的表面結(jié)構(gòu)，包括三角面片的頂點坐標(biāo)和面法向量信息。同樣的STL樣本模型，采用三角面片數(shù)量的不同將導(dǎo)致模型的精確度不同，三角面片個數(shù)越多，模型就越精細(xì)，但是相應(yīng)的計算時間也將增加。本文所有實驗，簡化三級固體火箭發(fā)動機模型采用103 682個三角面片描述其幾何信息。

SolidWorks軟件建立的多材質(zhì)模型STL文件無法存入材料信息，本文采用如下方法來完成模型幾何信息和材質(zhì)信息的相互匹配：每導(dǎo)入一個固體火箭發(fā)動機STL配件由用戶在材料信息表中選擇材料信息和相關(guān)的衰減指數(shù)，與模型幾何信息數(shù)據(jù)一同寫入仿真系統(tǒng)，從而將表示多材質(zhì)的固體火箭發(fā)動機STL模型裝配到仿真系統(tǒng)中進行計算。

2.3 膠片特性仿真

本文采用膠片黑度值[10]是來表征膠片的特性仿真，因此確定仿真膠片的特性曲線是至關(guān)重要的。在實際檢測中測定膠片特性曲線[11]，普遍使用的方法是，使用X射線源，固定管電壓、使用濾波板，以保證輻射質(zhì)量，然后再固定管電流、焦距、增感鉛屏厚度等條件，僅改變曝光時間t這一參數(shù)，測定膠片的底片黑度D與曝光量E(表示管電流乘以曝光時間t)的關(guān)系，為顯示斜率用D-lgE表示膠片特性曲線，該曲線近似為一條簡單的上直斜線。在實際的固體火箭發(fā)動機射線檢測中，射線源參數(shù)采用的是劑量值，而當(dāng)管電壓和濾波等其他條件固定時，吸收劑量值和曝光量成正比，因此認(rèn)為膠片的黑度和吸收劑量也近似成正比，于是本文采用D-lgK特性曲線[12]來表征膠片特性仿真，其中K為達到黑度D時的吸收劑量，D-lgK曲線也近似為一條上斜直線，故計算仿真黑度值：

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

(2)

式中a為近似膠片特性直線斜率；b為矯正參數(shù)；c為截距。這3個參數(shù)由膠片特性所決定，可根據(jù)膠片特性曲線調(diào)整。

3 投影成像算法及加速

國內(nèi)外在射線檢測仿真系統(tǒng)上成像方法主要有兩種：一種是使用體素模型進行光線追蹤求交算法，但是當(dāng)精度較高、模型體積大、數(shù)據(jù)量大時，仿真速度很慢，不符合固體火箭發(fā)動機模型的應(yīng)用；另一種是使用表面模型進行光線求交算法，計算效率要比體素模型高，但是針對固體火箭發(fā)動機模型體積大數(shù)據(jù)量大的問題，求交速度依然是個瓶頸。本文重點研究基于STL表面模型的高效求交算法，提高射線與大數(shù)據(jù)量樣本的求交速度。

3.1 成像算法描述

由X射線源發(fā)射出一系列錐形的X射線，穿透固體火箭發(fā)動機樣本模型，撞擊到固體火箭發(fā)動機后面的面陣探測器的各個探元上，該過程中X射線穿透固體火箭發(fā)動機并與其發(fā)生相互作用而衰減，這是投影成像的基本思想。

投影成像算法是仿真系統(tǒng)的核心，由以下三層循環(huán)組成：(1)射線源點陣循環(huán)，對每個離散的焦點計算投影；(2)固體火箭發(fā)動機模型的零件循環(huán)，通過X射線與每個零件相交計算，求出衰減距離；(3)像素矩陣循環(huán)，對探測器每個像素計算經(jīng)過模型衰減后的X射線能量值。

在每次循環(huán)中首先要求出射線與固體火箭發(fā)動機模型三角形的交點，再對交點排序求出模型的衰減距離，其中關(guān)鍵的求交方法如下：

射線直線方程參數(shù)形式：

(3)

式中 (x0，y0，z0)為射線上一點；(m，n，p)為射線的方向向量。

三角形所在平面方程參數(shù)形式：

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

(4)

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

先判斷射線直線與三角平面是否平行，如果不平行將直線方程式(3)代入平面方程式(4)，求出交點，然后判斷交點是否在三角形內(nèi)。方法如下：

對于平面內(nèi)任意一點如圖5所示，都可由如下方程來表示：

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

(5)

圖5 點與三角形的關(guān)系示意圖Fig.5 Schematic of the relationship between points and triangles

如果讓P位于三角形ABC內(nèi)部(如圖5所示)，u和v必須滿足如下三個條件：

3.2 三維空間投影約束的求交加速算法

射線與STL模型求交要求遍歷模型中的所有三角面片來完成直線與三角形求交運算，導(dǎo)致仿真計算耗時長的問題。針對這一問題，黃魁東等[13]采用將模型八叉樹剖分的思想來加速射線與模型的求交運算，但是針對尺寸大、數(shù)據(jù)量大、結(jié)構(gòu)復(fù)雜的固體火箭發(fā)動機模型，八叉樹剖分模型需要耗很長時間，這對于仿真系統(tǒng)的實際應(yīng)用需求無法滿足。

本文對模型采用三維空間平面投影約束的方法快速濾除絕大多數(shù)與射線不相交的三角面片，從而加速遍歷運算。將射線和三角面片頂點從三維空間中投影到任意兩個二維平面上，這樣形成了二維平面上3個點和線的關(guān)系，如圖6所示。判斷3個投影頂點是否在射線投影直線的同一側(cè)，如果在同一側(cè)說明射線與該三角面片肯定不相交，增加這個約束條件可避免射線與絕大多數(shù)不相交三角面片的求交計算，可節(jié)約大量的計算時間。

圖6 二維投影平面上3個點和線的關(guān)系示意圖Fig.6 The relationship between three points and line on

三維空間投影約束的求交加速算法如下：

針對模型所有三角面片，首先分別確定射線和三角面片3個頂點在XOY平面和XOZ平面的垂直投影直線方程和點坐標(biāo)；然后判斷這3個點是否在投影直線的同一側(cè)，若都不在同一側(cè)則計算出空間射線方程和三角面片的平面方程；其次根據(jù)計算出的射線方程和平面方程判斷射線是否與平面相交，若相交則求出交點坐標(biāo)；最后判斷交點是否在三角面片內(nèi)，若交點在三角面片內(nèi)則表示射線與模型相交，記錄相交點坐標(biāo)信息用以后續(xù)計算最終的射線衰減距離。

對一條射線與模型所有三角面片遍歷求交點并排序計算出衰減距離消耗的時間做實驗，三維空間投影約束前后加速實驗結(jié)果見表1，實驗結(jié)果表明加約束條件后可以提高10倍左右的加速比。

表1 約束算法加速效果對比

3.3 CUDA并行加速

在射線投影計算過程中涉及到射線源不同離散焦點位置、探測器不同探元位置的射線與所有三角面片的求交運算，同時為了保證仿真投影結(jié)果的精度，要求仿真探測器分辨率大，計算量龐大。經(jīng)過實驗測試在CPU Intel Core i7-6700K(八核，主頻4.00 GHz)的計算機上對模型所有三角面片遍歷加約束條件后計算出衰減距離一次需要10 ms左右，生成4096×4096大小的投影，理論估計需要233.02 h，因此加約束條件后計算速度的提升依然無法滿足實際應(yīng)用。考慮到在投影計算過程中每條射線計算均是獨立的，本文利用CUDA技術(shù)對核心算法進行加速，CUDA技術(shù)[14]是英偉達公司使用顯卡GPU來輔助電腦CPU分擔(dān)計算任務(wù)的技術(shù)，充分利用顯卡的高性能并行計算能力，可提高算法的運算速度。CUDA線程并行是GPU細(xì)粒度的并行，線程并行運行的核函數(shù)都是相同的，CUDA流并行可實現(xiàn)對同一個核函數(shù)傳遞不同的參數(shù)，實現(xiàn)任務(wù)級別的并行。利用流并行機制可實現(xiàn)不同固體火箭發(fā)動機模型零件級別的并行，線程并行實現(xiàn)探測器探元級別上的并行計算。

對比CPU在不同GPU設(shè)備上做了計算速度實驗，實驗參數(shù)：射線源離散個數(shù)為5個；探測器大小為4096×4096，像素規(guī)格0.1 mm；殼體，27 840個三角面片；絕熱層，38 400個三角面片；推進劑，37 442個三角面片；總計103 682個三角面片。實驗結(jié)果見表2。

表2 CUDA加速效果對比

表2結(jié)果表明，在不同GPU硬件加速情況下可達到500～1000倍的加速比，GPU的性能越好，加速比將越大。基于CUDA 的算法加速，使得該方案能滿足于實際固體火箭發(fā)動機射線檢測工程的速度需求。

4 軟件編程與實驗結(jié)果

4.1 系統(tǒng)開發(fā)

根據(jù)固體火箭發(fā)動機射線檢測仿真系統(tǒng)的功能分析和總體設(shè)計，系統(tǒng)主控平臺采用微軟公司MFC以文檔為中心的程序架構(gòu)思想，在Windows 7操作系統(tǒng)上，使用Visual studio 2013結(jié)合OpenGL三維可視化庫[15]編寫了固體火箭發(fā)動機射線檢測仿真系統(tǒng)主界面，見圖7。

圖7 仿真系統(tǒng)主界面Fig.7 Simulation system main interface

4.2 實驗結(jié)果

固體火箭發(fā)動機簡化三層模型殼體、絕熱層、藥柱前肩位置的衰減距離可視化結(jié)果見圖8。從仿真結(jié)果上看，對模型的物理結(jié)構(gòu)仿真非常準(zhǔn)確。

圖8 殼體、絕熱層、藥柱的衰減距離可視化圖Fig.8 Attenuation distance visualization of case,insulation and grain

根據(jù)實際試驗測出0.02 Gy吸收劑量對應(yīng)膠片2黑度，本文采取式(2)中a=1，c=0，計算得出b=5000。不同劑量下的黑度仿真可視化結(jié)果及黑度值分布曲線見圖9。

在相同條件下對仿真膠片與實際膠片測量的黑度值作了對比實驗。建立的簡化固體火箭發(fā)動機仿真模型與真實發(fā)動機的尺寸等參數(shù)有所差別，定量誤差分析要求黑度值的測量位置參數(shù)一致，因此本文采用定性分析黑度值仿真結(jié)果。參見圖10。

(a)劑量為3 Gy

(b)劑量為10 Gy

(c)劑量為25 Gy

(a)固體火箭發(fā)動機實際膠片與投影仿真結(jié)果

(b)劑量為10 Gy在相同位置處

如圖10所示，仿真黑度分布曲線趨勢與實際膠片測得黑度曲線非常相似，在10 Gy劑量條件下仿真黑度數(shù)值大小與實際系統(tǒng)同劑量下膠片測得的黑度數(shù)值相近。在實際檢測固體火箭發(fā)動機缺陷時更關(guān)注絕熱層的黑度值，實驗結(jié)果表明在絕熱層的黑度仿真與實際測量曲線基本吻合，最大差值為0.23，與文獻[10]的黑度值最大誤差0.3相比，本文的黑度值仿真在檢測范圍內(nèi)更準(zhǔn)確，證明了本文黑度仿真模型的可靠性，絕熱層黑度對比結(jié)果如圖11所示。

圖11 絕熱層處仿真與實際黑度值對比Fig.11 Comparison of simulation and actual blackness values at the thermal insulation

5 結(jié)論

本文圍繞固體火箭發(fā)動機的射線檢測流程和原理，設(shè)計了非點光源與多能譜的固體火箭發(fā)動機射線檢測仿真系統(tǒng)。在保證投影精度的前提下，通過加三維空間投影約束的STL模型求交算法和GPU并行技術(shù)實現(xiàn)了固體火箭發(fā)動機投影仿真與黑度仿真的工程應(yīng)用。實驗結(jié)果表明，該方案較之于CPU下STL模型求交算法有500～1000倍左右的加速比，仿真膠片與實際膠片黑度最大差值為0.23，與實際膠片測量的黑度更加接近，在可檢測范圍內(nèi)更準(zhǔn)確。該仿真系統(tǒng)為固體火箭發(fā)動機射線檢測的工藝制定和優(yōu)化搭建了一個穩(wěn)定的仿真實驗平臺，具有很大的實用價值。