張凱洋,賴建偉,黃志勇
(海軍工程大學, 武漢 430030)
algorithm
固體火箭發(fā)動機因其結構簡單、使用維護方便,廣泛應用于軍事領域。與液體火箭發(fā)動機以及吸氣式發(fā)動機不同,固體發(fā)動機在生產(chǎn)過程中推進劑已直接澆筑于燃燒室內(nèi),無法使用拆解的方法對其進行全面檢查或檢測。因此,如何在不破壞發(fā)動機的前提下,對其內(nèi)部可能存在的缺陷進行快速、準確檢測,對于保持導彈武器作戰(zhàn)狀態(tài)、提高導使用安全性、降低日常檢測和維護成本具有重要作用。
目前,用于固體火箭發(fā)動機無損檢測的方法包括超聲、射線、激光、紅外、微波和ICT[1]。與其他無損檢測技術相比,ICT檢測技術具有直接得到物體內(nèi)部結構的數(shù)字圖像,檢測結果的存儲、傳輸、分析和處理方便等突出特點。目前,ICT技術已經(jīng)廣泛應用于各類工業(yè)領域,在對小型固體火箭發(fā)動機的檢測過程中也展現(xiàn)出了良好的檢測效果。但對于直徑在2 m以上的大型固體火箭發(fā)動機,采用ICT檢測方法,即使進行一個發(fā)動機斷層面的重建就需要耗費大量時間,而完成發(fā)動機整體結構檢測需要完成上千個斷層面的重建。ICT檢測技術應用于此類大型固體發(fā)動機檢測,這存在時間及人力成本過大的問題。而在此類大型工件的ICT檢測過程中,有時并不需要關心其內(nèi)部結構只需對工件外部環(huán)形區(qū)域進行檢測即可。史源源等[2]對固體火箭發(fā)動的切向CT重建技術進行了研究,采用基于待測物體先驗圖像的IRR迭代算法(iterative reconstruction reprojection,IRR)實現(xiàn)了切向投影數(shù)據(jù)的圖像重建。熊亮等研究了基于最大似然期望最大化的ML-EM算法(maximum likelihood expectation maximization,ML-EM)。曾理等[4]針對石油管道的成像問題通過引入圖像先驗信息,提出了基于凸集投影全變分最小化的TVM-POCS算法(total variation minimization- projection onto convex sets,TVM-POCS)得到了較好的成像效果。
傳統(tǒng)ICT檢測采用全斷面掃描成像,完成檢測需要掃描并重建固體火箭發(fā)動機整個斷面。與之相比,切向ICT檢測只采集穿過物體外部的射線投影數(shù)據(jù),在圖像重建過程中也僅重建待測物體外部環(huán)形區(qū)域。在對大型固體火箭發(fā)動機進行檢測時,采用切向ICT檢測,可以在滿足對脫粘等缺陷進行檢測的前提下降低對設備硬件條件的要求,提高檢測效率,因此具有較高的應用價值。
以窄扇束TR掃描ICT系統(tǒng)為例,其對固體火箭發(fā)動機進行檢測過程中采用平移/旋轉掃描方式。如圖1(a)所示,在進行TR掃描過程中,射線源與探測器以某一固定步長從待檢發(fā)動機一端掃描至另一端,完成一次掃描后射線源與探測器復位,檢測系統(tǒng)根據(jù)設定參數(shù)將待檢發(fā)動機旋轉一固定角度,完成旋轉后射線源與探測器再進行一次平移掃描,依此往復,直至得到待檢發(fā)動機360°范圍內(nèi)的全部掃描數(shù)據(jù)。而采用切向ICT檢測方法時,在射線源與探測器同步平移的過程中,不需要掃過整個發(fā)動機截面,只需要掃描發(fā)動機殼體部分。射線源與探測器同步平移一定步長后即復位,待檢發(fā)動機旋轉固定角度后射線源與探測器進行下一輪平移掃描,直至得到待檢發(fā)動機360°范圍內(nèi)的切向掃描數(shù)據(jù),如圖1(b)所示。
如圖1(a)陰影部分所示,對于傳統(tǒng)的TR全局掃描方法在圖像重建過程中得到的圖像為待測發(fā)動機完整的斷面圖像。而切向掃描方法只對待測發(fā)動機外部環(huán)形區(qū)域斷面圖像進行重建,重建區(qū)域如圖1(b)陰影部分所示。
圖1 窄扇束TR全局掃描/切向掃描示意圖
切向ICT檢測方法存在的主要問題在于如何在采集數(shù)據(jù)基礎上完成待測物體外部環(huán)形區(qū)域的圖像重建。由于相較于全局掃描,切向ICT檢測只掃描并采集待測物體的部分投影數(shù)據(jù),使得在此基礎上的圖像重建成為一類投影數(shù)據(jù)不完備條件下的圖像重建問題。其解決難點在于待測物體的投影數(shù)據(jù)是截斷的,雖然理論上存在有唯一解,但是由于計算機數(shù)值計算誤差以及噪聲的影響,重建圖像沿徑向的邊緣會出現(xiàn)嚴重的偽影。對于切向ICT圖像的重建,現(xiàn)有重建算法中較典型的包括:ART代數(shù)迭代算法(algebraic reconstruction techniques,ART),IRR迭代算法,ML-EM算法以及基于凸集投影全變分最小化的TVM-POCS算法。
圖像重建過程中可以將原始圖像視為一個n×n大小的灰度圖像矩陣F,將圖像矩陣F按行展開得到圖像向量f,圖像重建問題即求解如下方程組,即:
(1)
式(1)中:M為投影射線總數(shù);pj為第j條射線穿過待測物體后的投影值;wi, j為投影系數(shù)反映第i個像素對第j條射線的影響,在數(shù)值上可以取為第j條射線與第i個像素塊的相交長度。式(1)寫為矩陣形式,有:
Rf=P
(2)
式(2)中:R為M×N系統(tǒng)投影系數(shù)矩陣;f為由灰度圖像矩陣F展開后的圖像向量;P為由對應射線投影值組成的投影數(shù)據(jù)向量。由于單條射線相交的像素有限,因此系統(tǒng)投影系數(shù)矩陣R為一大型稀疏矩陣。因此,方程(2)在求解利用常規(guī)的矩陣理論難以進行求解。
ART迭代算法其核心思想就是將系統(tǒng)矩陣R按行拆開,將每一行看作一個n維超平面。然后從初始值開始,依次計算迭代值在每一個超平面上的垂直投影,直到迭代收斂至所有超平面的共同交點。利用ART迭代算法求解方程(1),每次迭代過程可以表示為:
(3)
當完成M次投影得到f(M)為一次完整迭代,之后將f(M)作為初始值代回,重復上述過程直到結果滿足終止條件。在ART算法的基礎上,Andersen等[5]提出的聯(lián)合代數(shù)重建算法SART(simultaneous algebraic reconstruction techniques,SART)。與ART算法不同,SART算法將初始值并行投影至所有超平面,然后計算所有投影落點的平均值,以此作為下一次投影的迭代值。相較于ART算法,SART算法增強了抗干擾能力,并擁有了進行并行計算的可能。
IRR 是一種基于濾波反投影的迭代算法,由于外部問題的投影數(shù)據(jù)不完備,因此需要引入虛擬探測器并對采集得到的投影數(shù)據(jù)進行補充[6]。IRR 算法的本質在于通過迭代,不斷改善對虛擬探測器投影數(shù)據(jù)的估計,使其更加接近真實情況。虛擬探測器布置示意圖如圖2。
圖2 虛擬探測器布置示意圖
設真實探測器得到的數(shù)據(jù)為P1,虛擬探測器的數(shù)據(jù)估計為P2,令P0={p|p∈P1∪P2}。將補完得到的投影數(shù)據(jù)代入式(2)中,即有:
Rf=P0
(4)
P0=Rf1
(5)
(6)
完成式(4)、式(5)、式(6),即完成一次迭代過程,由于每次迭代的投影數(shù)據(jù)均完全包含了原始數(shù)據(jù),因此每次迭代后結果不會更加偏離真實情況。重復以上迭代過程直到滿足終止條件得到重建圖像。
(7)
式(7)中,ni為測量過程中產(chǎn)生的噪聲信號。一般研究認為測量產(chǎn)生的誤差噪聲近似適用于泊松噪聲模型,即:
(8)
式(8)取對數(shù),得到極大似然代價函數(shù),即fi由pi得到的概率為:
(9)
在式(9)基礎上求得fi,使得ψ(f) 取得極值,求得的fi即為圖像重建的最優(yōu)估計值。上述方程使用解析方法難以求解,一般采用EM方法迭代求解。其迭代過程數(shù)學表示為:
(10)
由于ML-EM 算法在重建過程中既考慮了系統(tǒng)的物理模型,又考慮了觀測數(shù)據(jù)的統(tǒng)計特性,重建出的圖像要優(yōu)于 FBP重建出的圖像。雖然ML-EM算法在一定迭代中可以取得不錯的效果,但是在實際使用中卻存在一個問題,即當觀測數(shù)據(jù)被噪聲污染比較嚴重時,重建圖像的效果并不是完全隨著迭代的進行而越來越好,而是當?shù)螖?shù)達到一定值之后,重建圖像的質量會下降[6]。
TVM-POCS算法最早由Sidky等人應用到稀疏角CT圖像重建問題中。Zeng L等將此方法應用于石油管道切向CT投影數(shù)據(jù)的處理,提出了基于子區(qū)域平滑的凸集投影總變差最小化迭代算法。對于灰度圖像,其圖像的全變分定義為:
(11)
TVM-POCS算法的目的在于求得解,使得其在滿足式(2)基礎上使得圖像的全變分最小,即:
(12)
在實際求解過程中,由于測量誤差等因素的影響不一定存在完全滿足式(2)的求解結果。因此在實際求解過程中允許求解結果存在一定誤差,即:
(13)
TVM-POCS算法的實現(xiàn)包含2個步驟:POCS步主要保證重建結果滿足數(shù)據(jù)保真項;TVM步主要是使重建結果的TV減小。兩步交替迭代進行,最終得到重建圖像,求解過程如下。
步驟1數(shù)據(jù)初始化,令f(0)=0;
步驟2將初始值代入式(3)中,利用ART算法進行一次完整的迭代運算,得到fART;
步驟3對fART進行非負約束,即
(14)
步驟4計算增量因子d,即
(15)
步驟5計算圖像全變分梯度大小,圖像fART的全變分梯度大小為:
(16)
步驟6如式(17)所示,利用vi, j對圖像進行全變分最小化修正,其中α為調(diào)節(jié)因子,即
(17)
利用仿真模體的計算機仿真掃描數(shù)據(jù),可以驗證以上圖像重建算法在固體火箭發(fā)動機切向ICT檢測圖像重建過程中的有效性并比較以上算法各自的優(yōu)缺點。
本文仿真實驗采用Matlab進行編程,程序運行環(huán)境采用I7-9750H處理器,RAM為16 GB,利用SART算法進行迭代計算過程中調(diào)用GPU加速運算過程。仿真實驗主要參數(shù)如表1所示。
表1 切向ICT檢測仿真實驗參數(shù)設置
仿真模體及其缺陷設置如圖3所示,全局ICT檢測仿真掃描結果如圖4(a)所示,為驗證切向ICT重建算法,在斷面重建過程中只采用邊緣部分探測器采集的數(shù)據(jù),所采集數(shù)據(jù)如圖4(b)所示。
圖3 仿真模體及其缺陷設置示意圖
圖4 仿真掃描結果示意圖
利用圖4(b)所示的切向掃描數(shù)據(jù),分別采用ART、IRR、ML-EM及TVM-POCS等4種算法,對仿真模體外部環(huán)形區(qū)域進行重建,重建結果如圖5所示。
圖5 仿真數(shù)據(jù)重建效果圖
如圖6所示,將以上重建圖像中對應仿真模體的缺陷部分進行局部放大,可以較直觀地比較不同成像算法的成像效果。
從圖6中可以發(fā)現(xiàn),相較于IRR算法及ML-EM算法,ART算法和TVM-POCS算法在對切向ICT掃描數(shù)據(jù)進行圖像重建時效果更好。由于IRR算法引入了人為設定的虛擬探測器數(shù)據(jù),在環(huán)形重建區(qū)域的切向上產(chǎn)生了嚴重的偽影。切向偽影的存在使得采用IRR算法重建圖像難以對徑向拓展的缺陷進行有效檢測。ML-EM算法較好地還原了仿真模體內(nèi)的缺陷特征,但重建圖像存在嚴重的棋盤狀噪點,這容易導致被測物體內(nèi)可能的微小缺陷被噪聲掩蓋。
圖6 仿真模體缺陷成像效果圖
(18)
重建結果與原始圖像的均方誤差越大,說明重建結果的誤差越大。記圖像內(nèi)所有像素中的最大灰度值為fmax,重建結果與原始圖像的峰值信噪比定義為:
(19)
SSIM是一種衡量兩幅圖像相似度的指標,其計算公式為:
(20)
式(20)中:μ、μ′分別為F、F′的均值;σ2、σ′2分別為F、F′的方差;σ″為F和F′的協(xié)方差。設L為圖像的灰度范圍,則c1=(0.01L)2,c2=(0.03L)2。
當其他條件相同時,重建結果與原始圖像的均方誤差越小,峰值信噪比越大,結構相似性越大則說明重建效果越好。4種重建算法與原始圖像的均方誤差、峰值信噪比和結構相似性如表2所示。
表2 4種算法重建圖像的RMSE、PSNR和SSIM
通過比較4種算法所得到的重建圖像與原始圖像的RMSE、PSNR和SSIM數(shù)值可以發(fā)現(xiàn),TVM-POCS算法在以上3種評價指標上均優(yōu)于其余3種算法,這與圖5、圖6所示重建圖像的直觀比較效果是一致的。
對于切向ICT掃描得到的數(shù)據(jù),ART算法、IRR算法、ML-EM算法及TVM-POCS算法均可以較好地重建出被測物體斷面圖像。從圖6仿真模體內(nèi)的模擬脫粘圖像上看,切向ICT檢測可以滿足固體火箭發(fā)動機界面脫粘的檢測要求。
從上述4種切向ICT重建算法的仿真重建結果上看,TVM-POCS算法的圖像重建結果優(yōu)于其他3種算法。ART算法及ML-EM算法的主要缺點在于圖像噪點較多,影響對缺陷的識別判斷,IRR算法的主要缺點在于重建圖像在切向上偽影嚴重,由于脫粘在圖像上是切向延伸的,因此切向偽影將導致無法根據(jù)重建圖像準確計算脫粘缺陷的實際大小。
將ICT檢測技術應用于固體火箭發(fā)動機無損檢測,對于掌握發(fā)動機質量情況和變化規(guī)律,保持和提高部隊戰(zhàn)斗力水平具有重大意義。對于大型固體火箭發(fā)動機,切向ICT檢測對于發(fā)動機殼體各類復合材料之間脫粘問題的檢測是可行的。本文主要針對切向ICT圖像重建算法進行了梳理和研究,通過仿真實驗分析了各個算法的優(yōu)缺點,可為下一步切向ICT重建算法的研究提供參考。