盛亮，張艷紅，高建鵬，李陽(yáng)，段寶軍，黑東煒

（1 西北核技術(shù)研究所強(qiáng)脈沖輻射環(huán)境模擬與效應(yīng)國(guó)家重點(diǎn)實(shí)驗(yàn)室，西安710024）

（2 清華大學(xué)工程物理系，北京100084）

0 引言

近二十年來(lái)，激光、Z 箍縮驅(qū)動(dòng)聚變點(diǎn)火技術(shù)取得了長(zhǎng)足進(jìn)步［1-2］，氘氚燃燒區(qū)的空間分布反映了所有點(diǎn)火物理過(guò)程的最終狀態(tài)，精確診斷燃燒區(qū)空間結(jié)構(gòu)對(duì)于改進(jìn)壓縮對(duì)稱(chēng)性技術(shù)措施、校驗(yàn)物理模型、評(píng)估點(diǎn)火成敗等具有重要意義，不同能量的中子圖像能夠直接給出燃燒區(qū)以及周?chē)镔|(zhì)在高壓縮比條件下的空間分布，因此中子成像技術(shù)是研究脈沖射線(xiàn)源性能的核心技術(shù)［3-5］。呂敏院士、王奎祿研究員等在20世紀(jì)80年代建立了基于厚針孔、射線(xiàn)圖像轉(zhuǎn)換屏、選通型像增強(qiáng)器、高速圖像記錄等關(guān)鍵技術(shù)的中子、伽馬成像系統(tǒng)，這些基本技術(shù)思路幾十年來(lái)沒(méi)有發(fā)生根本性的變化，但針對(duì)測(cè)量對(duì)象特征的不同進(jìn)行了一些重要技術(shù)改進(jìn)，如針對(duì)低強(qiáng)度脈沖中子源，開(kāi)展了半影孔、環(huán)形孔成像與圖像復(fù)原技術(shù)研究［6-8］，發(fā)展建立了多分幅相機(jī)、條紋相機(jī)、光纖-光電倍增管陣列等三種不同的時(shí)空分布信息記錄系統(tǒng)［9，10］。單光軸中子成像獲得的二維圖像是輻射區(qū)域中子的沿著特定方向線(xiàn)積分投影，無(wú)法直接反映氘氚燃燒區(qū)域內(nèi)部的三維空間分布，如氘氚燃燒區(qū)體積等重要物理參數(shù)需要根據(jù)旋轉(zhuǎn)對(duì)稱(chēng)假設(shè)給出。最近在高能量密度物理研究領(lǐng)域，數(shù)值模擬程序由二維向三維發(fā)展，美國(guó)國(guó)家點(diǎn)火裝置（National Ignition Fusion，NIF）團(tuán)隊(duì)的研究結(jié)果表明，三維模擬程序計(jì)算的中子產(chǎn)額與實(shí)驗(yàn)結(jié)果最為接近，三維程序給出的燃燒區(qū)分布需要直接的實(shí)驗(yàn)結(jié)果驗(yàn)證［11］。借鑒計(jì)算機(jī)層析成像思想，通過(guò)多角度成像的方式可以給出氘氚區(qū)三維分布信息，但由于脈沖射線(xiàn)源一般具有自輻射強(qiáng)、持續(xù)時(shí)間短、射線(xiàn)強(qiáng)度變化范圍大等特點(diǎn)，這些特點(diǎn)使得不同光軸之間像素級(jí)精確匹配、高精度時(shí)間關(guān)聯(lián)、相機(jī)靈敏度差異校正等存在較大困難，一般建立的光軸數(shù)量非常有限（≤5），如美國(guó)國(guó)家點(diǎn)火裝置上歷經(jīng)十余年，完成了三條中子成像光軸建設(shè)［12］。這種基于非常有限光軸中子圖像重建三維空間分布是一件極具挑戰(zhàn)性的工作，有解空間大、存在偽影、無(wú)法保持邊緣特征的問(wèn)題。中國(guó)工程物理研究院江少恩等在1997年進(jìn)行了ICF 聚變等離子體的X 射線(xiàn)三維成像技術(shù)探索，建立了基于代數(shù)重建（Algebraic Reconstruction Technique，ART）的三維圖像復(fù)原算法［13］；美國(guó)VOLEGOV P L 等建立了最大期望（Expectation Maximization，EM）算法圖像重建方法［14］，以及三維中子圖像球諧函數(shù)分解（Spherical Harmonic Decomposion，SHD）、柱諧函數(shù)分解重建方法（Cylindrical Harmonic Deomposition，CHD）［15-16］，這種函數(shù)分解解析方法具有重建速度快、重建結(jié)果在某些限定條件下解唯一確定的優(yōu)點(diǎn)，但其表達(dá)能力受投影角度數(shù)量的限制較大。

本文以輻射源強(qiáng)度空間分布為高斯函數(shù)的橢球源為研究對(duì)象，分別介紹并實(shí)現(xiàn)了球諧函數(shù)分解和EM迭代兩種三維圖像重建算法。與NIF 團(tuán)隊(duì)重建算法設(shè)計(jì)思路不同，根據(jù)慣性約束聚變靶燃燒區(qū)結(jié)構(gòu)具有一定球?qū)ΨQ(chēng)性這一物理特點(diǎn)，設(shè)計(jì)了一種將球諧函數(shù)分解與EM 迭代相結(jié)合的重建算法。該方法首先以球諧函數(shù)分解對(duì)輻射源進(jìn)行三維重建，然后將其重建結(jié)果作為EM 算法的初值進(jìn)行迭代重建，其最終重建結(jié)果可以看作是EM 迭代算法在球諧函數(shù)分解重建結(jié)果約束下尋找到的最優(yōu)解。

1 球諧函數(shù)分解三維圖像重建算法

橢球脈沖射線(xiàn)源及x=0（短軸中心位置）處切片圖如圖2（a）、（b）所示，本文給出的所有切片圖都為同一位置。在球諧函數(shù)分解解析求解的過(guò)程中，為了減少積分投影圖像中高頻噪聲對(duì)結(jié)果的影響，在積分投影圖像的頻域中增加了低通高斯濾波器，重建圖像及切片圖如圖2（c）、（d）所示。為了驗(yàn)證噪聲對(duì)球諧函數(shù)分解解析方法重建結(jié)果的影響，對(duì)橢球源的積分投影圖像添加兩種水平的高斯噪聲，噪聲的標(biāo)準(zhǔn)差分別為積分投影圖像最大值的2%和5%，對(duì)于噪聲投影圖像重建結(jié)果分別如圖2（e）、（f）和圖2（g）、（h）所示。通過(guò)與無(wú)噪聲重建結(jié)果的切片圖比較可以看到，在有噪聲的情況下，解析方法重建結(jié)果中心部分產(chǎn)生了較大畸變，而邊緣部分的重建偽影部分受噪聲影響不大。

圖2 球諧函數(shù)分解三維圖像重建結(jié)果Fig.2 The three-dimensional image reconstructed results by spherical harmonic decomposition

無(wú)噪聲重建結(jié)果的積分投影圖像與理想的橢球源積分投影圖像（149°，307°）的差值|Δ|及Δ的一維分布直方圖如圖3（a）、（b）所示，其中Δ為重建投影圖像與理想投影圖像之間的強(qiáng)度差，由于大部分差值數(shù)據(jù)為Δ=0，為了更好顯示差值分布，在進(jìn)行直方圖統(tǒng)計(jì)時(shí)沒(méi)有將Δ=0 的結(jié)果統(tǒng)計(jì)在內(nèi)。2%和5%噪聲水平的重建結(jié)果差值及直方圖分別如圖3（c）、（d）和圖3（e）、（f）所示?？梢钥吹诫S著噪聲水平的增加，非零差值的直方圖分布變的更寬，理想情況下差值圖像的直方圖分布的最可幾值在Δ=0 附近，當(dāng)噪聲為5%時(shí)，差值圖像的直方圖分布較理想情況下出現(xiàn)了較大偏移，高斯擬合結(jié)果給出的最可幾值出現(xiàn)在Δ=0.38 附近。

圖3 不同噪聲水平下重建三維圖像投影圖像與理想投影圖像之間的差值及其直方圖Fig.3 The difference image and its histogram between the projection image of reconstructed 3D image and ideal 3D image with different noise

2 EM 迭代算法

在中子成像中，物像距遠(yuǎn)大于輻射源自身尺寸，受厚針孔視角限制，射線(xiàn)經(jīng)過(guò)厚針孔的有效發(fā)射角度遠(yuǎn)小于1。成像系統(tǒng)獲取的中子圖像是具有三維空間結(jié)構(gòu)的脈沖輻射源沿著某一方向的積分投影圖像，其矩陣形式為

式中，投影數(shù)據(jù)p為M×1 的向量，源強(qiáng)度f(wàn)為N×1 的向量，投影矩陣A為M×N的矩陣。EM 算法是一種統(tǒng)計(jì)類(lèi)代數(shù)迭代重建方法，可用于處理投影角度少、數(shù)據(jù)不完備的問(wèn)題，較ART、MART 等代數(shù)類(lèi)迭代算法具有更好的重建效果。受中子成像探測(cè)數(shù)量限制和投影噪聲的影響，EM 算法隨著迭代次數(shù)增加，迭代噪聲增強(qiáng)，易引起圖像質(zhì)量退化［20］。

EM 迭代算法的基本思想為：假設(shè)射線(xiàn)源放射強(qiáng)度f(wàn)服從泊松分布，利用觀(guān)測(cè)投影數(shù)據(jù)p，通過(guò)EM 迭代尋找使似然函數(shù)的期望最大化的參量f，似然函數(shù)定義為

式中，Ai是A矩陣第i行。在EM 算法中，依賴(lài)于一組完全觀(guān)測(cè)數(shù)據(jù){zij}，數(shù)據(jù)的每個(gè)元素表示探測(cè)像素i探測(cè)到第j個(gè)輻射源體素內(nèi)發(fā)射的光子數(shù)，服從Poisson 分布

式中，aij是A矩陣第i行、j列元素，fj表示向量f第j列的元素值，投影數(shù)據(jù)p與觀(guān)測(cè)數(shù)據(jù)z的關(guān)系為

式中，pi表示p第i行的元素值。結(jié)合式（8）和（10），目標(biāo)函數(shù)為

采用EM 算法分為兩個(gè)步驟。

E-step：已知當(dāng)前源圖像fk，投影值p，觀(guān)測(cè)數(shù)據(jù)的條件期望為

將式（11）的目標(biāo)函數(shù)的隨機(jī)變量zij用期望值代替，則有

M-step：使條件期望最大化，對(duì)目標(biāo)函數(shù)求偏導(dǎo)。

令式（14）為零，得到經(jīng)典的EM 算法迭代格式［21-22］

將每個(gè)投影方向數(shù)據(jù)作為一個(gè)子集，采用有序子集最大期望重建算法（Ordered Subset Expectation Maximization，OSEM）可以加速算法收斂速度［21-22］，本文中的算例迭代5 次即可收斂至兩次迭代差別可忽略的程度。在利用EM 算法進(jìn)行圖像重建時(shí)易產(chǎn)生高頻噪聲，可以利用中值濾波（Median Filter，MF）、全變差最?。═otal Variation Minimum，TVM）等約束方法對(duì)迭代數(shù)據(jù)進(jìn)行平滑約束。文獻(xiàn)［14］引入圖像空間分布的先驗(yàn)信息，使用Gibbs 能量分布函數(shù)對(duì)EM 算法的迭代過(guò)程進(jìn)行約束，提高了算法的噪聲抑制能力。EM算法獲得的是局部最優(yōu)解，通常在利用EM 迭代算法時(shí)，初值設(shè)為0 或者1。如果能夠?qū)M 算法的初值限定在一定范圍內(nèi)，就有可能獲得與受初值約束的局部最優(yōu)解。大多數(shù)對(duì)EM 算法解空間的約束依賴(lài)于投影圖像數(shù)據(jù)的先驗(yàn)假設(shè)，如高斯先驗(yàn)、脈沖噪聲先驗(yàn)、結(jié)構(gòu)化先驗(yàn)等［23］。激光、Z 箍縮驅(qū)動(dòng)慣性約束聚變研究中，為了獲得更高的能量轉(zhuǎn)換效率，靶丸的內(nèi)爆壓縮過(guò)程盡可能為對(duì)稱(chēng)壓縮，其結(jié)構(gòu)上具有一定的球?qū)ΨQ(chēng)性，基于此類(lèi)物理先驗(yàn)信息可以幫助我們限定EM 類(lèi)算法的初值，將重建結(jié)果限定在物理結(jié)構(gòu)先驗(yàn)約束范圍內(nèi)，將球諧函數(shù)分解重建結(jié)果作為EM 迭代算法的初值是ICF 中子圖像三維重建的一種自然選擇。當(dāng)然，根據(jù)重建對(duì)象的不同，EM 迭代算法的初值約束可以有不同的選擇，比如輻射源具有一定柱對(duì)稱(chēng)性，可以利用阿貝爾變換（Abel Transform）或者文獻(xiàn)［16］提出的柱諧函數(shù)分解重建結(jié)果作為EM 迭代算法的初值。

圖4（a）～（j）為無(wú)噪聲情況下圖1所示5 個(gè)投影角度橢球源OSEM、MF-OSEM、球諧函數(shù)分解初值約束的MF-OSEM 算法三維圖像重建結(jié)果其切片圖。圖4（b）中OSEM 算法重建圖像給出的切片圖存在較多迭代偽影，包含輻射源內(nèi)部的高頻噪聲和邊緣噪聲，圖4（d）中MF-OSEM 算法的結(jié)果對(duì)高頻噪聲進(jìn)行了較好的抑制；MF-OSEM、球諧函數(shù)初值約束的MF-OSEM 重建三維圖在與理想的橢球源積分投影圖像（149°，307°）的差值及其直方圖分別如圖5（a）、（b）和圖5（c）、（d）所示。從直方圖中看球諧函數(shù)約束的MF-OSEM算法差值圖像較MF-OSEM 算法絕對(duì)值更小，分布也更為均勻，受投影角度的影響也更小一些，直方圖的強(qiáng)度分布都集中在Δ=0 附近。圖5（e）、（f）和圖5（g）、（h）分別為添加2%和5%噪聲之后，球諧函數(shù)初值約束的MF-OSEM 算法重建投影圖像與理想投影圖像之間的差值及其直方圖，與球諧函數(shù)分解結(jié)果類(lèi)似，隨著噪聲水平的提高，差值圖像的絕對(duì)值也變大，直方圖分布也變的更寬，與圖3（e）差值圖像邊緣與理想圖像差別較大不同，差值圖像圖5（g）分布較為均勻，圖5（h）的直方圖高斯擬合給出的最可幾值出現(xiàn)在Δ=0.29附近。

圖1 橢球輻射源5 個(gè)角度投影Fig.1 The 5 projection views of the ellipsoidal radiation source

圖4 EM 類(lèi)算法三維圖像重建結(jié)果Fig.4 The EM algorithm reconstructed results

圖5 不同噪聲水平下球諧函數(shù)初值約束的MF-OSEM 算法重建三維圖像投影圖像與理想投影圖像之間的差值及其直方圖Fig.5 The difference image and its histogram between the projection image of the MF-OSEM reconstructed 3D image with initial values constrained by SHD and ideal 3D image with different noise

3 重建結(jié)果量化比較

在重建圖像結(jié)果評(píng)價(jià)中，通常采用一些量化指標(biāo)進(jìn)行評(píng)價(jià)。本文采用峰值信噪比（Peak-Signal to Noise Ratio，PSNR）、均方根誤差（Root Mean Squared Error，RMSE）、KL 偏差（Kullback-Leibler divergence）［24］等3 種重建圖像質(zhì)量評(píng)價(jià)方法對(duì)幾種圖像重建結(jié)果進(jìn)行定量比較。

設(shè)理想圖像為X，重建圖像為Y，P×Q為圖像大小，PSNR 定義為

Xmax為圖像X的最大值（本文中為1.0）。RMES 定義為

KL 偏差DKL定義為

從式（16）和（17）中可以看出，重建結(jié)果評(píng)價(jià)指標(biāo)中PSNR 越大，RMSE 越小。KL 偏差越小，表示重建圖像與理想圖像之間越接近。文中涉及到的幾種重建算法量化比較結(jié)果如表1所示。

從表1 的結(jié)果可以看出，隨著噪聲水平的提高，所有算法重建結(jié)果都會(huì)變差，其中OSEM 算法由于在重建過(guò)程中沒(méi)有對(duì)噪聲進(jìn)行抑制，其重建結(jié)果對(duì)噪聲最為敏感。球諧函數(shù)分解方法雖然依賴(lài)于投影角度數(shù)量，在本文5 個(gè)投影角度情況下，對(duì)橢球源依然保持了較好的重建精度。MF-OSEM 算法和基于球諧函數(shù)分解的MF-OSEM 算法整體效果較好，采用球諧函數(shù)重建結(jié)果作為MF-OSEM 的迭代初值可以把EM 類(lèi)算法的迭代結(jié)果限制在初值附近。在本文的量化評(píng)價(jià)指標(biāo)下，該方法的重建結(jié)果與理想圖像之間的差別最小，對(duì)噪聲適應(yīng)性也較好。球諧函數(shù)分解算法為解析方法，占用計(jì)算資源較小，在處理器intel（R）core（TM）i7-4790 CPU@3.6GHz，內(nèi)存為24GB 的計(jì)算資源上完成一次計(jì)算的時(shí)間約為6 分鐘。EM 類(lèi)迭代算法占用了大部分計(jì)算時(shí)間，幾種EM 類(lèi)迭代算法所需時(shí)間差別不大。本文中EM 算法迭代5 次左右即可收斂，所需時(shí)間約為1 455 s（～24 min）。需要說(shuō)明的是，本文計(jì)算的三維結(jié)構(gòu)總像素?cái)?shù)為100×100×100，如果要重建更多像素?cái)?shù)圖像，則需要在較高性能服務(wù)器上進(jìn)行計(jì)算以縮短計(jì)算時(shí)間。

表1 幾種三維重建算法效果比較Table 1 The comparison of these three-dimensional image reconstruction algorithms

4 結(jié)論

三維圖像數(shù)據(jù)是校驗(yàn)慣性約束聚變靶物理設(shè)計(jì)的重要結(jié)果，極少角度投影三維圖像重建存在解空間過(guò)大的問(wèn)題。本文建立了強(qiáng)度空間分布為高斯函數(shù)理想橢球源5 個(gè)投影角度球諧函數(shù)分解、EM 迭代等三維重建算法。球諧函數(shù)分解解析重建方法表達(dá)能力受投影角度數(shù)量的限制，EM 迭代算法只具有局部最優(yōu)解，基于慣性約束聚變靶為具有一定球?qū)ΨQ(chēng)性的射線(xiàn)源的認(rèn)識(shí)先驗(yàn)，將球諧函數(shù)分解解析重建結(jié)果作為EM 迭代算法的初始值，從而把EM 迭代算法重建結(jié)果約束在球諧函數(shù)分解解析結(jié)果附近。重建結(jié)果及量化評(píng)價(jià)指標(biāo)表明，該算法具有良好的效果，對(duì)噪聲的適應(yīng)性也較好，可作為類(lèi)似輻射源極少角度投影三維重建的基準(zhǔn)比較算法。下一步，將在實(shí)驗(yàn)室構(gòu)建相應(yīng)的多光軸成像系統(tǒng)，通過(guò)實(shí)際實(shí)驗(yàn)結(jié)果來(lái)進(jìn)一步校驗(yàn)算法的有效性。