王朝華 陳德峰 劉榮海

（1．國網(wǎng)河南省電力公司電力科學(xué)研究院，河南 鄭州 450052；2．檢測成像北京高校工程研究中心，北京 100048；3．云南電網(wǎng)有限責(zé)任公司電力科學(xué)研究院，云南 昆明 650000）

0 引言

CT系統(tǒng)工作時，X射線被前準(zhǔn)直器準(zhǔn)直成一個狹窄的射束，通過被檢測的工件后被探測器組接收，得到多個角度下射線的強(qiáng)度信號，這些電信號經(jīng)過光電倍增管或光二極管放大后，再經(jīng)過模數(shù)轉(zhuǎn)換器轉(zhuǎn)換成數(shù)字信號，通過計算機(jī)運(yùn)用一定的數(shù)學(xué)方法處理后，在顯示器上得到一個重建的、完整的二維層析圖像。由于在同一斷面上各點(diǎn)材料的密度不同，因此各點(diǎn)的線性衰減系數(shù)也不相同，這樣，吸收X射線的量也不同，這些差異可在圖像上顯示出來，為研究、判斷分析提供了可靠的依據(jù)。

國內(nèi)外現(xiàn)有工業(yè)X射線CT還無法實現(xiàn)對高速旋轉(zhuǎn)物體（例如高速運(yùn)轉(zhuǎn)中的航空發(fā)動機(jī)、高速離心機(jī)以及燃?xì)廨啓C(jī)等）的CT成像[1-2]。有關(guān)高速旋轉(zhuǎn)物體CT成像的研究尚處在探索階段。2003年提出了基于隨機(jī)正交投影識別和重排的高速旋轉(zhuǎn)物體CT成像方法，并對該方法進(jìn)行了仿真數(shù)據(jù)驗證，但是該方法要求所采集的數(shù)據(jù)仍是準(zhǔn)靜止且不發(fā)生數(shù)據(jù)混疊的。2014年相關(guān)學(xué)者分析了研制航空發(fā)動機(jī)在線成像CT系統(tǒng)面臨的難點(diǎn)問題，給出了系統(tǒng)設(shè)計、數(shù)據(jù)采集方案以及圖像重建方法，提出了降低由數(shù)據(jù)混疊引起的圖像模糊方法，并構(gòu)建高能CT系統(tǒng)進(jìn)行模擬成像測試。高速旋轉(zhuǎn)物體CT成像的主要難點(diǎn)是受X射線源流強(qiáng)和探測器數(shù)據(jù)采集速度的限制，無法在物體高速旋轉(zhuǎn)的條件下采集準(zhǔn)靜止的CT數(shù)據(jù)，因此經(jīng)典的CT成像模型不再適用，如果仍采用經(jīng)典CT成像模型將導(dǎo)致重建圖像在沿旋轉(zhuǎn)方向出現(xiàn)模糊。

1 數(shù)學(xué)模型

CT常用算法大體上可分為3類，第一類為系列重建法或稱系列（級數(shù)）展開法，例如聯(lián)合迭代重建法；第二類為變換重建法，例如傅里葉變換重建法等；第三類為其他重建法。目前，變換重建法中的卷積反投影法應(yīng)用較多，常用的有平行束卷積反投影法和扇形束卷積反投影法。

該文在若干假設(shè)下，針對電網(wǎng)設(shè)備高速旋轉(zhuǎn)物體CT成像給出了數(shù)據(jù)采集時間的選取準(zhǔn)則，建立了射線混疊的CT成像數(shù)學(xué)模型，并提出了相應(yīng)的圖像迭代重建算法。當(dāng)數(shù)據(jù)存在噪聲時，射線混疊的CT成像數(shù)學(xué)模型具有高病態(tài)性。經(jīng)模擬CT數(shù)據(jù)和實際CT數(shù)據(jù)驗證，該文所提出的算法在數(shù)值上具有較好的魯棒性，可以正確重建高速旋轉(zhuǎn)物體的CT圖像，且不會在旋轉(zhuǎn)方向產(chǎn)生圖像模糊。

每個CT掃描數(shù)據(jù)Ij可以看作是入射流強(qiáng)I由被測物體x經(jīng)過一定調(diào)制過程Fj得到的數(shù)據(jù)，即Ij=Fj（I，x）；而CT圖像重建過程則可以看作是由一系列關(guān)系Ij=Fj（I，x）解調(diào)出x（j∈J，J是掃描數(shù)據(jù)指標(biāo)集）。在此意義上，該文的思想與利用編碼孔徑和利用準(zhǔn)直器提高分辨率的思想有相似之處，而不同之處在于相應(yīng)的調(diào)制過程和解調(diào)方法，類似的思想還有基于虛擬焦點(diǎn)的超分辨重建方法。該文對高速旋轉(zhuǎn)物體CT成像做了一些較為理想的假設(shè)，而實際系統(tǒng)可能會有旋轉(zhuǎn)不勻速、旋轉(zhuǎn)軸晃動的情況，但相應(yīng)的CT成像問題本質(zhì)上仍是調(diào)制解調(diào)問題，只是在調(diào)制關(guān)系Ij=Fj（I，x，Θ）中增加了參數(shù)族Θ，因此需要同時研究Θ的測量或估計方法和圖像重建（解調(diào)）算法。

假設(shè)：被檢測物體是勻速轉(zhuǎn)動的；旋轉(zhuǎn)軸不發(fā)生晃動；探測器的采樣時間是連續(xù)可選的；X射線是單能的；忽略多次散射光子以及電子對淹沒產(chǎn)生的光子。在該假設(shè)下，高速旋轉(zhuǎn)物體斷層CT成像如公式（1）所示。

2 高速轉(zhuǎn)動數(shù)據(jù)快速處理方法

高速轉(zhuǎn)動圖像重建具有數(shù)據(jù)量大的特點(diǎn)，因此在工程應(yīng)用中考慮重建算法的計算效率也是CT技術(shù)的一個重要方面。目前對反投影部分的CUDA加速策略主要如下：1）在線程分配時考慮反投參數(shù)的Z軸相關(guān)性。每個線程各自完成對應(yīng)Z軸上一系列點(diǎn)的重建任務(wù)，這種線程分配方式可以減少GPU的重復(fù)計算量。2）紋理存儲器的使用。紋理存儲器能夠被緩存，因此使用紋理存儲器可以提高訪問速度，并且紋理存儲器帶有優(yōu)化過的尋址系統(tǒng)，支持方便、高效的二維插值操作，可用來完成反投影中的雙線性插值運(yùn)算。3）對全局存儲器釆用合并訪問優(yōu)化。在滿足合并訪問條件的情況下，只需要一次傳輸便可以處理來自一個half-warp的16個線程對全局內(nèi)存的訪問請求。4）對三角函數(shù)計算的優(yōu)化。在GPU中完成一次三角函數(shù)運(yùn)算需要32個時鐘周期，而一次乘加只需要4個時鐘周期，因此在計算三角函數(shù)時釆用遞歸方法，將不同角度下的三角函數(shù)運(yùn)算轉(zhuǎn)換為一次乘加運(yùn)算，便可在一定程度上縮短GPU的運(yùn)算時間。由于GPU架構(gòu)本身在高性能計算方面的優(yōu)勢，在應(yīng)用上述優(yōu)化策略后， 便可獲得CPU上百倍的加速比。在上述并行策略的基礎(chǔ)上，做了更細(xì)致的優(yōu)化，使優(yōu)化后的加速效果獲得大幅提升。

與上述優(yōu)化策略相比，該文所提出的方法與現(xiàn)有方法的不同主要體現(xiàn)在以下3個方面：1）線程分配方式。2）通過使用常數(shù)存儲器存儲計算反投影參數(shù)時的中間量，減少GPU的重復(fù)計算量。3）通過優(yōu)化一個kernel中反投影的張數(shù)來減少訪問全局存儲器的次數(shù)。

3 數(shù)值試驗

首先，采用仿真數(shù)據(jù)驗證該文提出的數(shù)據(jù)采集和重建算法的有效性。其次，利用實際工業(yè)CT系統(tǒng)采集的數(shù)據(jù)驗證該文算法對多能CT數(shù)據(jù)的有效性。

3.1 仿真數(shù)據(jù)驗證

根據(jù)以下條件模擬高速旋轉(zhuǎn)物體的扇形束CT數(shù)據(jù)：射線源為能量8 MeV的單能射線源；探測器為線陣列探測器，由512個探測器單元組成，各單元的尺寸為0.60 mm。被測物體模型為如圖1（a）所示的鐵材質(zhì)工件，環(huán)形的外徑為208.30 mm，內(nèi)徑為158.70 mm，其最長弦為134.91 mm，鐵對8 MeV的光子的線性衰減系數(shù)為0.023 63 mm-1。射線源到旋轉(zhuǎn)中心的距離為1 000.00 mm，射線源到探測器中心的距離為1 200.00 mm。射線混疊數(shù)據(jù)由旋轉(zhuǎn)17圈（6 120 °）采集的720 個投影數(shù)據(jù)得到，即每旋轉(zhuǎn)8.5 °采集1幅投影數(shù)據(jù)。掃描數(shù)據(jù)中加入泊松噪聲的方法如下：由已知光子流強(qiáng)I0和被測物體模型，根據(jù)公式（6）數(shù)值計算；再以為泊松分布的均值產(chǎn)生一個隨機(jī)值，并用該隨機(jī)值替換Ik，j。 由此可知，數(shù)據(jù)的噪聲水平與I0和物體吸收率是相關(guān)的。對給定物體來說，I0越小，數(shù)據(jù)噪聲水平越高。

圖1 理想仿真數(shù)據(jù)以及重建圖像

分別考察重建算法對無噪聲數(shù)據(jù)和有噪聲數(shù)據(jù)的重建效果。圖1分別為用3種方法對無噪聲的混疊掃描數(shù)據(jù)重建的CT圖像，重建圖像矩陣為512×512。使用FBP算法直接重建的圖像如圖1（a）所示，第 4.1 節(jié)所述的直接解調(diào)算法重建的圖像如圖1（b）所示，第3.2節(jié)提出的迭代算法重建的圖像如圖1（c）所示。圖1（a）在旋轉(zhuǎn)方向存在明顯的圖像模糊，遠(yuǎn)離旋轉(zhuǎn)中心的像素模糊程度相應(yīng)加重。在無噪聲的情況下，由直接解調(diào)方法和該文提出的迭代算法均能重建滿意的圖像。

3.2 實際數(shù)據(jù)驗證

該文在射線混疊CT數(shù)學(xué)模型中假設(shè)射線是單能的，然而工業(yè)CT的射線源所發(fā)出的射線是多能的。為驗證該文方法對多能射線的有效性，用工業(yè)CT設(shè)備掃描的數(shù)據(jù)進(jìn)行試驗驗證。采用延長探測器采樣時間的方法來模擬掃描數(shù)據(jù)混疊。試驗所用的工業(yè)CT系統(tǒng)探測器單元個數(shù)為2 604個，探測器單元寬度為0.25 mm，探測器的采樣時間為60 ms。射線源焦點(diǎn)到轉(zhuǎn)臺中心的距離為600.00 mm，射線源焦點(diǎn)到探測器的距離為1 139 mm。將葉輪放置于轉(zhuǎn)臺旋轉(zhuǎn)中心處進(jìn)行CT掃描。轉(zhuǎn)臺共旋轉(zhuǎn)7圈（即2 520°），探測器的采樣角度數(shù)為1 800 個，即每旋轉(zhuǎn)1.4 °采集1幅投影數(shù)據(jù)。葉輪實物照片和掃描數(shù)據(jù)如圖2所示。在葉輪中加入4根鉛筆芯作為標(biāo)記物，以考察重建圖像沿旋轉(zhuǎn)方向的模糊程度。

圖2 葉輪實物照片和混疊掃描數(shù)據(jù)

圖3（a）、 圖3（c）和 圖3（e）分別為直接FBP算法、直接解調(diào)算法和該文迭代重建算法重建的圖像，圖3（b）、圖3（d）和圖3（f）分別為它們的局部放大圖像。由直接FBP算法重建的圖像發(fā)生沿旋轉(zhuǎn)方向的模糊，直接解調(diào)算法已難以重建出葉輪的信息，而該文迭代重建算法仍可以重建出質(zhì)量較好的圖像。

圖3 葉輪實際掃描數(shù)據(jù)重建圖像

4 結(jié)語

該文在若干假設(shè)下，針對高速旋轉(zhuǎn)物體的CT成像問題建立了射線混疊的CT數(shù)學(xué)模型，提出了CT采樣時間的選取準(zhǔn)則及其相應(yīng)的迭代重建算法，并通過仿真數(shù)據(jù)和實際數(shù)據(jù)驗證了該方法的可行性。盡管該方法做了射線單能的假設(shè)，但是實際工業(yè)CT數(shù)據(jù)驗證表明，將該方法用于多能數(shù)據(jù)也具有很好的近似性。該文是對高速旋轉(zhuǎn)物體CT的初步研究，在物體轉(zhuǎn)速有擾動以及旋轉(zhuǎn)軸有晃動等情況下，還需要進(jìn)一步深入研究如何重建高速旋轉(zhuǎn)物體的CT圖像。