PBX稀疏CT投影PICCS圖像重建算法

張建偉，張才鑫，陳 華，張 韜，4，劉 晨，劉豐林

（1.中國工程物理研究院化工材料研究所，四川 綿陽 621999；2.重慶大學(xué)光電技術(shù)及系統(tǒng)教育部重點(diǎn)實(shí)驗(yàn)室，重慶 400044；3.重慶大學(xué)工業(yè)CT 無損檢測教育部工程研究中心，重慶 400044；4.哈爾濱工業(yè)大學(xué)計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，黑龍江 哈爾濱150001）

1 引言

高聚物黏結(jié)炸藥（Polymer Bonded Explosive，PBX）是一種由炸藥晶體、黏結(jié)劑等構(gòu)成的多組份復(fù)合材料，因其性能優(yōu)良廣泛應(yīng)用于武器裝備。PBX 在加工成型過程中，不可避免地存在孔洞、氣泡等初始損傷［1］，這些損傷在應(yīng)力影響、存放環(huán)境等多因素相互作用下，易形成不同程度、不同形態(tài)、不同擴(kuò)展方式的裂紋［2-3］，影響 PBX 的力學(xué)性能、爆轟性能以及安全性能［4-5］。因此，對PBX 炸藥件初始損傷以及裂紋進(jìn)行檢測，將有利于保證PBX 性能以及運(yùn)輸儲(chǔ)存過程中的安全性。當(dāng)前，對于PBX 炸藥件初始損傷以及裂紋檢測主要采用的是計(jì)算機(jī)斷層成像（Computed tomography，CT）技術(shù)，在實(shí)際檢測中，CT 掃描過程需要采集一定數(shù)量投影數(shù)據(jù)，這個(gè)過程需要較長時(shí)間，檢測效率低。因此，亟需研究一種提高檢測效率的方法。稀疏投影CT 掃描的方式可減少投影采集時(shí)間提高檢測效率。但稀疏CT 投影是典型的不完備CT 投影，從不完備投影集重建圖像存在噪聲大、偽影嚴(yán)重、細(xì)節(jié)模糊等問題，影響對檢測物體內(nèi)部結(jié)構(gòu)和缺陷的觀察、判斷。因此，在稀疏投影提高檢測效率的前提下，選用合適的重建算法來保證重建圖像質(zhì)量至關(guān)重要。

CT 投影重建方法可分為兩類，一類是傳統(tǒng)解析算法，其最具代表性的是濾波反投影重建算法（Filtered back projection，F(xiàn)BP），但它對投影數(shù)據(jù)的完備性要較求高；另一類是迭代重建算法，其代表性的方法有代數(shù)重 建 算 法（Algebra Reconstruction Technique，ART）［6］和聯(lián)合代數(shù)重建算法（Simultaneous Algebra Reconstruction Technique，SART）［7］。 ART 和 SART算法放寬了對投影數(shù)據(jù)完備性的要求，但計(jì)算量大，重建速度慢，噪聲較大。壓縮感知（compressed sensing，CS）［8］理論的出現(xiàn)為 CT 圖像重建提供了新的思路。Sidky［9］利用圖像在梯度變換域內(nèi)具有稀疏性這一先驗(yàn)知識，提出TVM-POCS（Total Variation Minimization-Projection Onto Convex Sets，TVM-POCS）重建方法，并成功用于稀疏角圖像重建，但在重建圖像輪廓附近存在嚴(yán)重的滑坡偽影。之后對上述算法進(jìn)行改進(jìn)提出自適應(yīng)最速梯度下降-凸集投影算法（Adaptive Steepest Descent Projection Onto Convex Sets，ASD-POCS）［10］，以罰函數(shù)的形式增強(qiáng)解的穩(wěn)定性，有效抑制了條形偽影，但收斂速度慢。Chen 等［11］通過在重建模型中引入先驗(yàn)圖像，提出了基于先驗(yàn)圖像約束的壓縮感知圖像重建方法（Prior Image constrained compressed Sensing，PICCS），與上述算法相比，該方法收斂速度快，噪聲小，它通過引入先前掃描的CT 圖像作為先驗(yàn)信息，可以大幅提高重建圖像質(zhì)量，使從更稀疏的投影集中恢復(fù)高質(zhì)量圖像成為可能。

為了提高檢測效率，本研究利用稀疏采樣的CT 掃描方式獲得投影數(shù)據(jù)，采用PICCS 重建算法對稀疏投影數(shù)據(jù)重建抑制噪聲，減少偽影，采用準(zhǔn)靜態(tài)單軸壓縮TATB 基PBX 的損傷試件驗(yàn)證采樣方式和重建算法，研究內(nèi)容包括：（1）CT 重建模型和PICCS 算法重建模型；（2）準(zhǔn)靜態(tài)單軸壓縮TATB 基PBX 的損傷試件的稀疏投影CT 掃描重建實(shí)驗(yàn)，以及對完備投影集FBP 重建結(jié)果與稀疏投影集SART、PICCS 兩種算法重建結(jié)果進(jìn)行對比研究。

2 PICCS圖像重建算法

CT 重建的數(shù)學(xué)模型［12］可以歸結(jié)為：

式中，A=(amn)表示系統(tǒng)矩陣M×N，M代表投影個(gè)數(shù)，N代表重建圖像的像素個(gè)數(shù)表示投影向量表示重建圖像向量，等式（1）可以通過迭代進(jìn)行求解。但重建稀疏投影數(shù)據(jù)時(shí)，M通常是比N小，等式（1）是欠穩(wěn)定系統(tǒng)，有無窮多組解。

正則化方法的提出有利于模型（1）從無窮多組解中找出滿意的最優(yōu)解，利用正則化方法，模型（1）可轉(zhuǎn)化為如下最優(yōu)化問題［11］：

PICCS 算法重建模型［11］可表達(dá)為如下：

式中，fP表示高質(zhì)量先驗(yàn)圖像向量，TV是稀疏算子，α(0 ≤α≤1)是平衡參數(shù)，由先驗(yàn)圖像質(zhì)量確定，先驗(yàn)圖像質(zhì)量越好，α越接近1。在PICCS 算法中，先驗(yàn)圖像的好壞將決定重建圖像的質(zhì)量，為了獲得高質(zhì)量的先驗(yàn)圖像，一般采用濾波反投影（FBP）算法對較完備的投影進(jìn)行重建。

采用PICCS 算法重建模型，稀疏算子選用全變分（Total Variation，TV）變換，TV 求解的是圖像梯度 1范數(shù)，其表達(dá)式［9］如下：

式中，I和J分別是重建圖像的長和寬，單位為 pixel；fi，j是重建圖像柵格(i，j)坐標(biāo)上的像素值。

3 準(zhǔn)靜態(tài)單軸壓縮試驗(yàn)與CT掃描重建實(shí)驗(yàn)

3.1 試驗(yàn)樣品

檢測樣品為Φ10 mm×15 mm 的 TATB 基 PBX 圓柱形試件，由中國工程物理研究院化工材料研究所提供。

3.2 準(zhǔn)靜態(tài)單軸壓縮試驗(yàn)

準(zhǔn)靜態(tài)單軸壓縮試驗(yàn)在Instro 5582 萬能材料試驗(yàn)機(jī)上完成，將檢測樣品放置于材料試驗(yàn)機(jī)上進(jìn)行單軸壓縮試驗(yàn)，位移控制速率0.05 mm·min-1，準(zhǔn)靜態(tài)加載至破壞。力學(xué)試驗(yàn)結(jié)束后將樣品從材料試驗(yàn)機(jī)取下，進(jìn)行離線CT 掃描重建實(shí)驗(yàn)。

3.3 稀疏投影CT掃描重建實(shí)驗(yàn)

為了驗(yàn)證PICCS 算法對于高聚物粘結(jié)炸藥不完備投影重建的可行性，分別采集不同稀疏程度的不完備投影集，用SART、PICCS 二種重建算法進(jìn)行重建，將它們與完備投影集的FBP 重建結(jié)果對比，最終得出PICCS 重建算法在滿足對檢測物體缺陷和內(nèi)部結(jié)構(gòu)觀察、判斷的前提下，可以用稀疏幾倍的投影集重建出滿足要求的圖像，進(jìn)而提高幾倍檢測效率。

CT 掃描重建實(shí)驗(yàn)所用設(shè)備如下，電腦：Intel（R）Core（TM）i7-67 CPU @ 3.40GHz、8.00GB 內(nèi)存；編程 環(huán) 境 ：MATLAB R2017a；CT 設(shè)備 ：nanoVoxel 設(shè)備（中國天津三英精密工程研究中心）。

實(shí)驗(yàn)掃描參數(shù)如表1。對其進(jìn)行準(zhǔn)靜態(tài)單軸壓縮試驗(yàn)，然后對試件進(jìn)行稀疏投影CT 掃描重建實(shí)驗(yàn)觀察其內(nèi)部裂紋分布。試件壓縮后采集數(shù)字X 線攝影圖（Digital Radiography，DR）如圖1。SART 算法初始圖像為零，PICCS 的先驗(yàn)圖像是通過掃描破壞前的試件采用SART-TV 迭代500 次得到，如圖2。

表1 掃描參數(shù)Table 1 Scanning parameters

圖1 CT 重建 實(shí) 驗(yàn) PBX 試件 DR 圖FFiigg.1 Digital radiography（DR）image of PBX Specimen reconstructed by CT scanning

圖2 試件破壞前CT 掃描重建得到的先驗(yàn)圖像FFiigg.2 Specimen image prior to its failure reconstructed by CT scanning

在實(shí)驗(yàn)中，把720 個(gè)投影FBP 重建結(jié)果作為參照圖像，稀疏投影CT 掃描分別均勻采集3 組不同投影數(shù)n=360，240，180 的投影集進(jìn)行重建實(shí)驗(yàn)，角度增量分別是d=1°，1.5°，2°，分別對應(yīng) FBP 投影集稀疏倍數(shù)m=2，3，4，SART 迭代 100 次收斂，PICCS 迭代 20 次收斂。為了清晰地觀察重建圖像中細(xì)小裂紋，把重建結(jié)果適當(dāng)放大，如圖3 所示。從圖3 可以看出，SART 重建算法受噪聲影響較大，在稀疏3 倍時(shí)，細(xì)小裂紋被噪聲淹沒，但稀疏2 倍時(shí)，圖3（b）中箭頭所指細(xì)小裂紋可以被SART 重建出來，所以SART 可以用稀疏2 倍的投影集重建滿足要求的圖像；而PICCS 重建算法受噪聲影響較小，重建圖像比較平滑，在稀疏4 倍時(shí)，重建圖像出現(xiàn)少量偽影，箭頭所指裂紋被偽影覆蓋，無法滿足檢測要求，但稀疏3 倍時(shí)，PICCS 算法可以重建出箭頭所指的細(xì)小裂紋，圖像滿足科研觀察要求。主觀視覺判斷PICCS 算法用稀疏3 倍的投影重建結(jié)果與720 個(gè)投影FBP 重建結(jié)果相近。每個(gè)算法重建出滿足要求的圖像所需的投影數(shù)量以及所需的時(shí)間如表2 所示。雖然迭代重建時(shí)間要比FBP 重建時(shí)間長，但是隨著計(jì)算機(jī)硬件的提升以及GPU 并行計(jì)算的出現(xiàn)，迭代重建時(shí)間會(huì)大大減少。

表2 各算法重建結(jié)果滿足要求所需時(shí)間Table 2 Time required for the reconstruction of each algorithm to meet the requirements

為了能夠更清晰地看到PBX 細(xì)小裂紋的重建情況，更加準(zhǔn)確地判斷SART、PICCS 算法能夠利用稀疏幾倍的投影重建出滿足要求的圖像，分別把圖3 中m=3，4 時(shí)（d）（e）（f）（g）（h）（i）矩形區(qū)域放大，如圖 4所示。圖4 可以清晰地看到m=3 時(shí)，SART 重建圖像中箭頭所指裂紋被噪聲淹沒，不利于觀察和判斷，而PICCS 算法在m=3 時(shí)可以清晰重建出細(xì)小裂紋，但m=4 時(shí)，PICCS 算法重建圖像中箭頭所指裂紋模糊不清，所以，主觀視覺判斷PICCS 算法用稀疏3 倍的投影可以重建出滿足要求的圖像。

圖3 完備投影FBP 和稀疏投影不同算法重建結(jié)果對比Fig.3 Comparisons of reconstructed results between complete projection FBP and sparse projections algorithms

圖4 對應(yīng)圖3 中m=3，4 時(shí)矩形區(qū)域，不同方法重建結(jié)果局部放大圖Fig.4 Local enlargement of reconstructed results by different methods correspondingly to the rectangular region of m=3，4 in Figure 3

為了定量比較1/3 稀疏投影PICCS 和SART 的重建結(jié)果與720 個(gè)投影FBP 重建結(jié)果，分別計(jì)算它們的均方根誤差（Root Mean Squared Error，RMSE）［13］、峰值信噪比（Peak Signal to Noise Ratio，PSNR）［13］，表達(dá)式分別如式（5）（7）所示。

式 中 ，RMSE 是 觀 測 值Xobs，i與 真 值Xmodel，i偏 差 的 平 方和觀測次數(shù)n比值的平方根。RMSE 值越小，重建結(jié)果越精確。

式中，MSE 表示當(dāng)前圖像X和參考Y的均方誤差（Mean Square Error）［13］，I和J分別是圖像的高度和寬度，n為每像素的比特?cái)?shù)，一般取8，即像素灰階數(shù)為256，PSNR 單位是dB，數(shù)值越大表示失真越小，重建圖像質(zhì)量越好。

對于重建圖像，感興趣的是重建物體內(nèi)部結(jié)構(gòu)和缺陷，而重建圖像中重建物體之外的信息不關(guān)心，所以為了更精確的進(jìn)行定量計(jì)算，把1/3 稀疏投影PICCS和SART 的重建結(jié)果與720 個(gè)投影FBP 重建結(jié)果進(jìn)一步處理后僅保留重建物體信息，定量計(jì)算區(qū)域如圖3（g）虛線圓區(qū)域。在計(jì)算過程中，以720 個(gè)投影FBP 重建結(jié)果為參考圖像，計(jì)算結(jié)果如表3 所示。從計(jì)算結(jié)果得出，用1/3 稀疏投影重建，PICCS 重建結(jié)果比SART 重建結(jié)果更加接近參考圖像。

表3 720 個(gè)投影FBP 重建結(jié)果為參考圖像的定量計(jì)算結(jié)果Table 3 Quantitative calculation results of 720 projection FBP reconstruction results as reference image

上述圖像質(zhì)量評價(jià)指標(biāo)都需要參考圖像，為了進(jìn)一步分析重建結(jié)果，選取一個(gè)無參考圖像評價(jià)指標(biāo)信噪比（Signal to Noise Ratio，SNR）［13］，表達(dá)式如式（8）所示，GROI和SROI分別是圖像中感興趣區(qū)域的灰度均值和噪聲的方差（感興趣區(qū)域的方差），圖像信噪比越高圖像質(zhì)量越好，噪聲影響越小。為了比較1/3 稀疏投影 PICCS 和 SART 與 720 個(gè)投影 FBP 的重建結(jié)果質(zhì)量，分別取圖3（d）（e）（f）中 A 部分作為感興趣區(qū)域，計(jì)算結(jié)果如表4 所示，SART 算法用240 個(gè)投影重建的結(jié)果受噪聲影響較大，而PICCS 用240 個(gè)投影重建結(jié)果由于PICCS 算法引入高質(zhì)量先驗(yàn)圖像的約束抑制了噪聲，噪聲影響較小，重建圖像質(zhì)量較好。

表4 SNR 計(jì)算結(jié)果（各算法感興趣區(qū)域A）Table 4 SNR calculation results（region of interest A by different algorithms）

4 結(jié)論

（1）針對CT 檢測PBX 效率低的問題，采用稀疏投影CT 掃描的方法，即通過減少投影采集數(shù)量提高檢測效率。在稀疏投影CT 掃描重建實(shí)驗(yàn)中，通過采集1/3 稀疏的TATB 基PBX 的損傷試件投影進(jìn)行圖像重建，雖然迭代算法重建時(shí)間比FBP 長，但是隨著計(jì)算機(jī)硬件的提升以及GPU 并行計(jì)算的出現(xiàn)，迭代重建時(shí)間會(huì)大大減少，本文方法檢測效率還會(huì)進(jìn)一步提升，所以本文最終檢測效率提高2 倍左右。

（2）1/3 稀疏的PBX 投影是典型的不完備投影，本研究采用的PICCS 重建算法由于先驗(yàn)圖像的約束，重建過程中能抑制噪聲，減少偽影。實(shí)驗(yàn)結(jié)果表明，用1/3 稀疏的PBX 投影進(jìn)行重建，與經(jīng)典SART 重建算法相比，PICCS 重建算法可以重建出細(xì)小裂紋，重建結(jié)果與完備投影FBP 重建結(jié)果質(zhì)量更加接近，滿足對檢測物體內(nèi)部結(jié)構(gòu)和缺陷的觀察、判斷。

（3）CT 廣泛應(yīng)用于PBX 缺陷檢測以及細(xì)觀內(nèi)部呈現(xiàn)，所以對于稀疏投影CT 掃描重建實(shí)驗(yàn)中采集1/3稀疏投影檢測效率提高2 倍左右，因研究內(nèi)容以及精度要求不同而不同。