基于壓縮感知理論的高能閃光照相密度反演方法

蘆存博,盛云霄

(1.北京計(jì)算機(jī)技術(shù)及應(yīng)用研究所，北京 100854; 2.陸軍裝甲兵學(xué)院 演訓(xùn)中心，北京 100072)

高能閃光照相是在核禁試的情況下，進(jìn)行核武器研究的主要手段之一。高能閃光照相能夠驗(yàn)證以前在核試驗(yàn)情況下形成的核武器編碼，并能對(duì)不完善的編碼進(jìn)行修正，以之來推動(dòng)在核禁試情況下的核武器研究。

高能閃光照相成像系統(tǒng)的結(jié)構(gòu)由X射線源、準(zhǔn)直器、客體、圖像接收系統(tǒng)組成。準(zhǔn)直器主有兩個(gè)作用：(1)降低被照客體的動(dòng)態(tài)范圍，便于圖像接收系統(tǒng)接收?qǐng)D像；(2)降低成像平面上散射X射線強(qiáng)度，以便提高客體中心區(qū)域的對(duì)比度。閃光照相的圖像接收系統(tǒng)除了氯化銀的膠片外，還有CCD相機(jī)、磷探測(cè)器等。在閃光照相中，由光源發(fā)出的X射線穿透照相客體后到達(dá)圖像接收系統(tǒng)，可記錄X射線的空間強(qiáng)度分布，形成閃光照相圖像。

高能閃光照相利用X射線的穿透能力以及它與物質(zhì)相互作用的性質(zhì)，根據(jù)成像平面上單位面積的X射線能量空間分布，通過圖像分析來提取客體的結(jié)構(gòu)信息和密度信息，進(jìn)而確定客體的幾何性質(zhì)和物理性質(zhì)，從而對(duì)客體及其內(nèi)部結(jié)構(gòu)進(jìn)行定量和物理診斷。

利用高能閃光照相獲得客體的幾何界面位置以及空間密度分布對(duì)研究客體在爆轟過程中的行為具有重要意義。在核物理領(lǐng)域，為了研究爆炸物體在某一時(shí)刻的物理狀態(tài)，例如爆轟過程中物質(zhì)密度的變化，常使用高能X射線照射物體，并根據(jù)投影信息重建密度分布等內(nèi)部結(jié)構(gòu)信息。高能X射線照相診斷的客體信息 “隱藏”于高能X射線照相得到的投影圖像中。高能閃光照相屬于投影成像技術(shù)，密度圖像重建技術(shù)是高能閃光照相圖像處理技術(shù)中的重要研究?jī)?nèi)容[1]。

高能閃光照相圖像重建問題是一個(gè)典型的高維線性反演問題?，F(xiàn)有的圖像重建算法大致有3種[2]：(1)解析算法，例如Abel變換、Radon變換、濾波反投影(Filter Back Projection，F(xiàn)BP)等；(2)迭代算法，例如代數(shù)重建技術(shù)(Algebraic Reconstruction Technique，ART)、約束共軛梯度(Constrained Conjugate Gradient，CCG)法、偏微分方程(Partial Differential Equation，PDE)約束法等；(3)統(tǒng)計(jì)算法，例如極大似然模型和期望最大化(Maximum Likelihood-Expectation Maximum，ML-EM)方法、馬爾可夫鏈蒙特卡羅(Markov Chain Monte Carlo，MCMC)方法。

在高能閃光照相中，每次實(shí)驗(yàn)只有一個(gè)或兩個(gè)方向的投影數(shù)據(jù)，因此要求密度反演的客體必須具有軸對(duì)稱性。這是因?yàn)椋绻腕w是軸對(duì)稱的，那么在與對(duì)稱軸垂直的平面上從不同方向得到的射線投影圖像是等效的，從而只需1幅投影圖像即可進(jìn)行密度重建工作。如果客體不具有軸對(duì)稱性，就需要多次實(shí)驗(yàn)。因此，需要研究少數(shù)投影數(shù)據(jù)條件下的密度反演問題。但是國(guó)內(nèi)在非軸對(duì)稱客體的密度反演技術(shù)方面的研究較少。

壓縮感知(Compressive Sensing，CS)[3-9]由于考慮了信號(hào)的稀疏性或可壓縮性而成為一種有效的信源處理技術(shù)，其核心是利用測(cè)量矩陣實(shí)現(xiàn)了對(duì)原始信號(hào)從高維空間到低維空間的線性映射，然后在重建算法中利用信號(hào)的稀疏性/可壓縮性實(shí)現(xiàn)了原始信號(hào)的高概率精確重建。在高能閃光照相下，待重建客體的物質(zhì)密度圖像可以看做是一種簡(jiǎn)單的三維圖像。通常情況下，組成客體的物質(zhì)數(shù)量不會(huì)太多，并且物質(zhì)的分布有一定規(guī)律，因此在表征客體的三維密度圖像中，不同區(qū)域之間相同密度元素值的相似性可以在圖像重建算法設(shè)計(jì)中作為先驗(yàn)信息加以利用，其可以用CS的思想表征。

針對(duì)高能閃光照相密度重建算法，國(guó)內(nèi)在軸對(duì)稱客體方面開展了大量數(shù)值模擬工作[1-2， 10-12]。但是目前還沒有關(guān)于非軸對(duì)稱客體方面的工作，對(duì)于基于CS理論的重建算法的研究也較少。利用先驗(yàn)信息約束重建圖像并結(jié)合CS理論而形成的重建算法將開辟一條閃光照相密度反演的新路徑。利用這種算法獲得高質(zhì)量的重建圖像，對(duì)提高核武器初級(jí)診斷閃光實(shí)驗(yàn)的客體界面和密度測(cè)量的精度具有重要的現(xiàn)實(shí)意義。

現(xiàn)有的利用CS思想的全變差(Total Variation，TV)類算法[13-14]雖然考慮了圖像的局部相似性，但沒有考慮圖像的非局部相似性，圖像的先驗(yàn)稀疏信息利用還不充分，圖像的重建精度還有進(jìn)一步提升的空間。為了解決這一問題，本文將組稀疏模型集成于TV框架之下，提出了一種基于組稀疏正則化(Group-Sparsity Regularization，GSR)的全變分重建技術(shù)(TV-GSR)，充分利用圖像的先驗(yàn)稀疏信息，利用客體的上、下、左、右4點(diǎn)對(duì)稱性降低圖像重建的規(guī)模，增加了重構(gòu)精度，加快了重建速度。仿真實(shí)驗(yàn)表明，文中提出的TV-GSR算法同時(shí)考慮了客體圖像的局部相似性和非局部自相似性，充分利用了圖像的先驗(yàn)稀疏信息，提升了圖像在無(wú)噪聲和有噪聲情況下的重建精度，對(duì)于高能閃光圖像和紋理細(xì)節(jié)豐富的CT(Computed Tomography)圖像都有良好的效果，具有普適性。

1 閃光照相圖像重建的數(shù)學(xué)模型

在高能閃光照相中，雖然客體處于爆轟壓縮的高速運(yùn)動(dòng)狀態(tài)，但是由于X光持續(xù)時(shí)間特別短，在此時(shí)間內(nèi)客體的壓縮運(yùn)動(dòng)位移及密度變化非常小，因此可以認(rèn)為在X光持續(xù)時(shí)間內(nèi)的客體處于靜止?fàn)顟B(tài)[15]?？紤]單能、無(wú)模糊、無(wú)散射的理想簡(jiǎn)化情況，閃光照相的投影圖像可以認(rèn)為是沿一組平行光束進(jìn)行采集的。閃光照相成像過程可以簡(jiǎn)述為：由光源發(fā)出的X光穿透照相客體后到達(dá)圖像探測(cè)平面，由圖像接收器件記錄X光的空間強(qiáng)度分布，形成照相圖像。閃光照相成像過程可以表示為

(1)

式中，I(x，y)為X光的透射強(qiáng)度；I0為X光的入射強(qiáng)度；μ(x，y，z)為照相客體質(zhì)量吸收系數(shù)的空間分布；ρ(x，y，z)為照相客體的空間密度分布；z為X光的傳播方向。令p(x，y)=ln[I0/I(x，y)]表示X光穿透客體的光程，則式(1)可以變換為

(2)

式(2)即為閃光照相投影成像模型，其矩陣形式為

p=Aμ1

(3)

2 組稀疏模型

(4)

式中，操作符/表示兩個(gè)矩陣元素之間的除法；1Ps×m為與fGk相同大小的全1矩陣。

令DGk為每一個(gè)結(jié)構(gòu)組fGk的自適應(yīng)字典。在此模型中，定義e為f的估計(jì)，在每次迭代中都可以得到e。令eGk為結(jié)構(gòu)組fGk的相應(yīng)估計(jì)。一旦得到eGk，對(duì)其應(yīng)用奇異值分解(Singular Value Decomposition，SVD)，可得

(5)

(6)

式中，DG表示所有DGk的級(jí)聯(lián)；αG表示所有αGk的級(jí)聯(lián)。

3 TV-GSR算法的具體描述

TV-GSR算法包括兩個(gè)模塊，分別為TV算法模塊和GSR算法模塊，這兩個(gè)模塊是相互獨(dú)立的，級(jí)聯(lián)到一起的。將GSR模型集成于TV框架之下，對(duì)每次迭代TV的重建結(jié)果使用提出的GSR模型進(jìn)行正則化處理，最后把GSR正則化的結(jié)果作為下一次TV迭代的初始圖像，這樣循環(huán)往復(fù)直到滿足停止準(zhǔn)則為止。本文中TV算法使用文獻(xiàn)[18]中的算法，下面重點(diǎn)介紹GSR算法模塊。

GSR的優(yōu)化問題可以表示為

(7)

(8)

但是，式(8)很難求解，因?yàn)閘0范數(shù)優(yōu)化問題是非凸的。本文用分裂Bregman迭代(Split Bregman Iteration， SBI)算法[17]求解這個(gè)問題?？紤]如下約束優(yōu)化問題

(9)

根據(jù)SBI算法，式(9)中的最小化問題可以被分裂成下列3個(gè)子問題。

(10)

(11)

bt+1=bt-(xt+1-Gyt+1)

(12)

對(duì)于SBI算法，將相關(guān)的初始值設(shè)置為t=0，選擇μ>0，b0=0，x0=0，y0=0。式(10)～式(12)循環(huán)往復(fù)直到滿足停止準(zhǔn)則為止。

(13)

(14)

(15)

然后，式(8)所示的最小化問題被轉(zhuǎn)換成關(guān)于f和αG的兩個(gè)子問題。對(duì)于給定的αG，式(13)所示的關(guān)于f的子問題是一個(gè)嚴(yán)格的二次凸優(yōu)化問題，其可以被定義為

(16)

通過設(shè)置P1(ft+1)的梯度為0，可以獲得式(16)的閉環(huán)求解，其為

(17)

式中，E是單位矩陣。對(duì)于給定的f，式(14)所示的關(guān)于αG的子問題可以被定義為

(18)

式中，et+1=ft+1-bt。根據(jù)文獻(xiàn)[19]，式(18)可以被轉(zhuǎn)化為

為減少道岔故障對(duì)有軌電車運(yùn)營(yíng)帶來的影響，一般宜在 3 個(gè)方向設(shè)置單渡線供道岔發(fā)生故障后臨時(shí)折返。雙Y道岔在各種故障情況下臨時(shí)運(yùn)行情況如圖 4 所示。

(19)

式中，ξ=(λ×Ps×m×n)/(μ×N)。然后，式(18)轉(zhuǎn)化為關(guān)于所有結(jié)構(gòu)組fGk的n個(gè)子問題。由于fGk=DGkαGk，eGk=DGkβeGk，每一個(gè)結(jié)構(gòu)組fGk的最小化問題可以被定義如下

(20)

因此，式(20)的一個(gè)閉環(huán)求解為

(21)

式中，hard(·)為硬閾值函數(shù)。一旦對(duì)于所有結(jié)構(gòu)組的αGk被計(jì)算出來，關(guān)于αG的子問題的最終求解也能被確定。

鑒于目前武器診斷的應(yīng)用實(shí)際，客體被假設(shè)滿足上、下、左、右4點(diǎn)對(duì)稱(參照方向?yàn)槌跏嫉纳渚€方向)。重建三維客體只需要重建出上半部分的客體即可，下半部分的客體可以利用上下對(duì)稱性得到。

TV-GSR算法的具體實(shí)現(xiàn)可描述為：利用TV對(duì)每次迭代得到的數(shù)據(jù)進(jìn)行正向投影到圖像域，再對(duì)TV重建的相應(yīng)二維圖像切片的左半部分圖像和右半部分圖像進(jìn)行加權(quán)平均，使用提出的GSR模型對(duì)加權(quán)平均后的圖像進(jìn)行正則化處理，最后利用左右對(duì)稱性把GSR正則化的結(jié)果恢復(fù)出原始大小的圖像，并將此圖像作為下一次TV迭代的初始圖像，這樣循環(huán)往復(fù)直到滿足停止準(zhǔn)則為止。

4 實(shí)驗(yàn)仿真及分析

本章節(jié)通過簡(jiǎn)單三維客體和復(fù)雜CT圖像上的數(shù)值仿真，對(duì)比了TV-GSR算法、TV[18]算法、SART算法和SART-GSR[17]算法的密度重建精度。投影方式采用多角度的投影，在0°～180°之間均勻取值，例如當(dāng)投影圖像為3幅時(shí)，投影角度為0°、60°、120°。在無(wú)噪聲重建中，對(duì)比算法為TV算法、SART算法和SART-GSR算法；在有噪聲重建中，投影數(shù)據(jù)被高斯噪聲污染，對(duì)比算法為TV算法。

4.1 簡(jiǎn)單三維客體實(shí)驗(yàn)

圖1 原始三維客體密度圖像的二維展開圖Figure 1. Two-dimensional expanded structure of original three-dimensional object

實(shí)驗(yàn)1在無(wú)噪聲的情況下，當(dāng)投影角度的數(shù)量變化時(shí)，不同重建算法的重建均方根誤差對(duì)比曲線如圖2所示。

圖2 不同投影角度數(shù)量下的重建均方根誤差Figure 2. Reconstruction root-mean-square error under different number of projection angle

如果客體整體的重建密度誤差在5%以下就滿足重建精度要求，7個(gè)方向和6個(gè)方向的投影圖像對(duì)于TV算法和TV-GSR算法分別來說基本是最佳的投影圖像數(shù)目。GSR的加入在保證成像質(zhì)量的情況下，有效改善了TV算法的性能，可以減少一個(gè)方向的投影圖像，從而降低輻射劑量。

實(shí)驗(yàn)2在有噪聲的情況下，當(dāng)投影角度的數(shù)量為5和9時(shí)，不同噪聲強(qiáng)度比例下的TV和TV-GSR重建算法的重建均方根誤差對(duì)比曲線如圖3(a)所示；當(dāng)投影角度的數(shù)量變化時(shí)，不同噪聲強(qiáng)度比例下的TV-GSR重建算法的重建均方根誤差對(duì)比曲線如圖3(b)所示。仿真圖中算法標(biāo)注括號(hào)里面的數(shù)字表示投影角度的數(shù)量。

(a)

(b)圖3 不同噪聲強(qiáng)度比例下的重建均方根誤差(a)投影角度的數(shù)量為5和9的情況 (b)不同投影角度數(shù)量的情況Figure 3. Reconstruction root-mean-square error under different percentage of noise strength(a)The case that the number of projection angle equals 5 and 9 (b)The case that the number of projection angle differs

從圖3可以看出，與TV算法相比，TV-GSR算法改善了客體的抗噪聲性能，甚至在高斯噪聲強(qiáng)度比例高達(dá)50%的情況下，客體整體的重建密度誤差仍然在10%以下。而對(duì)于TV算法來說，在投影角度的數(shù)量為5時(shí)，當(dāng)高斯噪聲強(qiáng)度比例為15%時(shí)，客體整體的重建密度誤差超過10%；而對(duì)于投影角度的數(shù)量為9時(shí)，當(dāng)高斯噪聲強(qiáng)度比例為20%時(shí)，客體整體的重建密度誤差超過10%。

4.2 復(fù)雜CT圖像實(shí)驗(yàn)

待重建客體為大小為128×128的CT二維圖像切片“pelvic_image_up”，其滿足左右對(duì)稱性，如圖4所示。

圖4 CT圖像切片F(xiàn)igure 4. CT image slice

在計(jì)算重建均方根誤差式時(shí)，相應(yīng)的Di為像素的真值，Ci為像素的重建值。在無(wú)噪聲的情況下，當(dāng)投影角度的數(shù)量變化時(shí)，對(duì)于算法TV-GSR、TV、SART和SART-GSR，不同投影角度數(shù)量下的重建均方根誤差對(duì)比曲線如圖5所示。在有噪聲的情況下，當(dāng)投影角度的數(shù)量為8時(shí)，不同噪聲強(qiáng)度比例下的TV和TV-GSR重建算法的重建均方根誤差對(duì)比曲線如圖6所示。

圖5 無(wú)噪聲情況下的重建均方根誤差Figure 5. Reconstruction root-mean-square error in noiseless scenario

圖6 有噪聲情況下的重建均方根誤差Figure 6. Reconstruction root-mean-square error in noisy scenario

從圖5和圖6可以看出，對(duì)于復(fù)雜CT圖像，在無(wú)噪聲和有噪聲的情況下，GSR的加入不僅有助于客體重建精度的改善，還改善了客體的抗噪聲性能。

從上述簡(jiǎn)單三維客體和復(fù)雜CT圖像上的數(shù)值仿真可以看出，提出的TV-GSR算法同時(shí)考慮了客體圖像的局部相似性和非局部自相似性，充分利用了圖像的先驗(yàn)稀疏信息，提升了圖像在無(wú)噪聲和有噪聲情況下的重建精度，對(duì)于高能閃光圖像和紋理細(xì)節(jié)豐富的CT圖像都有良好的效果，具有普適性。

5 結(jié)束語(yǔ)

持續(xù)提升閃光照相密度反演的精度是圖像分析研究人員的努力方向。本文提出了一種基于GSR的TV重建技術(shù)(TV-GSR)，對(duì)TV算法模塊和GSR 算法模塊進(jìn)行循環(huán)交替迭代，直到滿足停止準(zhǔn)則為止。本文采用上下對(duì)稱性來降低分層重建的規(guī)模；采用左右對(duì)稱性來降低GSR正則化的復(fù)雜度。通過仿真實(shí)驗(yàn)，定量給出了無(wú)噪聲情況下投影圖像數(shù)目和有噪聲情況下噪聲強(qiáng)度比例對(duì)密度重建精度的影響規(guī)律。基于紋理細(xì)節(jié)相對(duì)豐富的復(fù)雜CT圖像的實(shí)驗(yàn)結(jié)果也證明了引入GSR有助于改善客體重建精度。本文提出的TV-GSR算法同時(shí)考慮了客體圖像的局部相似性和非局部自相似性，提升了圖像在無(wú)噪聲和有噪聲情況下的重建精度，對(duì)于高能閃光圖像和復(fù)雜CT圖像都有良好的效果，具有普適性。