李化欣

(中北大學(xué) 信息與通信工程學(xué)院， 山西 太原 030051)

0 引 言

光學(xué)CT (光干涉CT、 光偏折CT等)是從多個方向上的投影數(shù)據(jù)進(jìn)行二維或三維被測場分布的反演. 與干涉CT相比， 光偏折 CT具有寬動態(tài)測試范圍， 適用于惡劣環(huán)境， 因此被廣泛用于測量復(fù)雜流場[1-2]. 莫爾偏折層析技術(shù)測試動態(tài)范圍大， 設(shè)備簡單、 對機械穩(wěn)定性要求低， 該技術(shù)適用于高溫高速流場的測量[3-5]. 然而， 偏轉(zhuǎn)層析成像的數(shù)學(xué)方面是復(fù)雜的， 正是這種復(fù)雜性阻礙了實際重建算法的發(fā)展. 計算偏折層析成像技術(shù)可分為變換方法和級數(shù)展開法. 當(dāng)重建使用全場視角范圍的大量等間隔投影數(shù)據(jù)時， 諸如濾波反投影方法的變換類算法可以獲得更好的結(jié)果， 但這種方法難以處理不完全數(shù)據(jù). 相比之下， 級數(shù)展開法[6]， 尤其是代數(shù)重建技術(shù)[7](ART)更加實用， 因為它可以靈活地處理不完全的數(shù)據(jù). 通過積分變換轉(zhuǎn)換投影數(shù)據(jù)后， 以ART算法為代表的迭代算法， 對于不完全投影數(shù)據(jù)， 極少投影數(shù)據(jù)， 抑制噪聲等方面比系列展開法更適合解決實際問題的重建. 近些年來出現(xiàn)的子集迭代重建算法[8]， 定向代數(shù)重建技術(shù)[9]、 自適應(yīng)聯(lián)合代數(shù)重建技術(shù)[10]等算法不斷改進(jìn)著重建算法的性能. 但是這些都是CT的一般重建算法， 不能反映光偏折CT的基本特征. 將光偏折投影數(shù)據(jù)沿垂直于投影方向的積分可以獲到光程差投影， 然后通過使用重建算法重建被測場. 但是， 在極少采樣數(shù)據(jù)的情況下， 這類重建算法也不能重建出滿足高精度指標(biāo)的溫度場. 在極少采樣情況下研究CT重建過程時， 壓縮感知(CS)對于信號處理是一種顛覆性的信息處理理論， 它開闊了研究人員的研究思路. 從Donoho和Candes提出的CS理論可以看出： 若信號或者圖像在某一特定的變換域存在稀疏表示， 那么可以適當(dāng)?shù)貎?yōu)化從極少的投影數(shù)據(jù)重建出原來的信號或者圖像[11]. 通過原始信號或者圖像的隨機投影而獲得投影數(shù)據(jù)， 該投影數(shù)據(jù)遠(yuǎn)小于奈奎斯特采樣定理所需的投影數(shù)據(jù). 這樣可以有效降低采樣系統(tǒng)的復(fù)雜性， 節(jié)約系統(tǒng)成本， 加快信號或者圖像的重建速度. 壓縮感知的理論已經(jīng)在數(shù)字圖像重建[12-13]、 圖像去除噪聲等領(lǐng)域得到了應(yīng)用， 而且該特性非常適合用于測量不穩(wěn)定[14]、 復(fù)雜流場的極度欠采樣重建.

本文針對偏折層析成像中極度欠采樣條件下重建精度較低的問題， 結(jié)合 CS理論， 將被測場的全變分作為稀疏性先驗?zāi)Ｐ停?利用最速下降法來調(diào)整全變分， 使得溫度場的全變分最小化. Sidky[15]等人提出的自適應(yīng)梯度下降法， 設(shè)置了很多參數(shù)來確定步長. 而本文在研究傳統(tǒng)ART-TV算法的基礎(chǔ)上， 稍微改進(jìn)了算法的步長， 并采用一致性約束步長作為梯度下降法的步長， 使梯度自適應(yīng)下降. 并把該方法用于重建模擬溫度場并分析重建誤差， 驗證算法的有效性. 對于稀疏投影數(shù)據(jù)的被測場重建， 本文改進(jìn)的算法重建質(zhì)量和抗噪聲能力比ART-TV算法更好.

1 原 理

1.1 偏折信息的轉(zhuǎn)化

如圖 1 所示， 設(shè)被測場某截面的折射率分布為

(1)

式中：n0為環(huán)境折射率,f(x,y)為相對折射率分布. 如圖 1 所示將坐標(biāo)系x-y逆時針旋轉(zhuǎn)θ角后得到坐標(biāo)系x′-y′. 當(dāng)平行于x′軸的入射光線經(jīng)過被測場時， 光沿著y′方向的偏折角為

(2)

該光線通過被測場的光程差為

(3)

式中：l是光線經(jīng)過被測場的光程. 當(dāng)光以小角度偏轉(zhuǎn)時， 光可以近似為被測場區(qū)域中的直線.

(4)

因此， 積分運算可以將偏折投影數(shù)據(jù)轉(zhuǎn)化為光程差投影數(shù)據(jù),即

(5)

圖 1 光線偏折原理與坐標(biāo)系示意圖Fig.1 Schematic diagram of the beam deflection and coordinate system

將重建區(qū)域劃分為k=N×N個網(wǎng)格， 將每個網(wǎng)格中的折射率變化設(shè)置為常數(shù)， 并設(shè)置總投影方向為q， 每一個方向的采樣數(shù)為r， 這樣總采樣數(shù)為m=q×r. 通過(5)可以將偏折投影方程轉(zhuǎn)換為光程差投影方程.

(6)

式中：Aij是第i條光線通過第j個網(wǎng)格的長度；pi是通過偏折角轉(zhuǎn)化得到的第i條光線的投影數(shù)據(jù).

1.2 ART-TV算法

壓縮感知從通過隨機投影獲得的少量投影數(shù)據(jù)中重建信號. 重建是在滿足觀測值的條件下找到最稀疏解的過程， 圖像或分布矩陣用f表示， 那么投影觀測值表示為P=Af，Ψ代表某一稀疏變換， 投影矩陣用A表示， 此處的Ψ和A是不相干的. 則通過求解式(7)就可重建原始分布.

(7)

式中零范數(shù)‖·‖0表示非零元素的數(shù)量. 式(7)是典型的非凸優(yōu)化問題， 因為此問題的求解過程相當(dāng)困難， 所以Candés和Donoho 提出了采用范數(shù)l1代替l0范數(shù)[16-17], 式(7)變?yōu)?/p>

(8)

根據(jù)壓縮感知的原理以及偏折角的變換， 用有限差分作為稀疏變換， 結(jié)合物理量(溫度、 密度)分布總變差的稀疏性， 通常通過總變差TV(Total Variation)來衡量， 并且式(8)中的‖Ψf‖1就成為‖f‖TV， 如式(9)所示

(9)

其中

式中：f表示折射率的分布矩陣；s,t是矩陣f的坐標(biāo)序號； 矩陣A是光線通過網(wǎng)格的長度. ART-TV算法的流程如下：

1) 取初始值f(0)=0， 從已知角度投影數(shù)據(jù)進(jìn)行ART迭代， 得到初始被測場的分布

j=1,2,…,k.

(10)

(11)

3) TV梯度下降

初始化TV梯度下降法

f(TV-GRAD)=f(ART-POS),

(12)

d=‖f(0)-f(ART-POS‖.

(13)

TV求導(dǎo)

(14)

(15)

式中：α為梯度下降法的松弛因子；d表示在正約束之后得到的折射率與不受約束的折射率的差值.

4) 將TV下降處理后的值賦值給ART作為迭代的初值， 進(jìn)行下一次迭代

f(0)=f(TV-GRAD).

當(dāng)d的差異足夠小或一定次數(shù)的迭代時， 迭代結(jié)束.

1.3 一致性控制步長算法

本文的算法是對TV最小化過程中梯度下降法中的步長進(jìn)行一致性控制， 可以根據(jù)被測場的重建質(zhì)量自動更新迭代步長， 并且可以確保梯度下降方向為被測場收斂的方向. 在初始迭代時， 投影數(shù)據(jù)存在很大的不一致性， 并且‖Af-P‖的值比較大， 所以可以增大迭代步長來使得f在可行域中； 當(dāng)‖Af-P‖的值比較小時， 保證在約束的可行域中找到最優(yōu)解. 所以我們希望通過數(shù)據(jù)的不一致性來確定梯度下降法中的步長. 一致性控制步長是基于被測場的原始投影數(shù)據(jù)與重建的投影數(shù)據(jù)之間的差， 以給出新的控制迭代步長.

在本方法中， 式(15)中的常數(shù)α可以由一致性控制步長函數(shù)η(n)來代替. 在本方法中， 修改后的梯度下降法為

(16)

其中

(17)

一致性控制步長函數(shù)η(n)在梯度下降的過程中不變， 只在ART開始迭代更新時修改. 式中f(0)是初始值， 在ART-TV開始迭代時n=1， 分子和分母的值一樣，η(n)=1， 隨著迭代次數(shù)的增多Af和P值逐漸接近， 因此， 步長從1開始， 并隨著迭代次數(shù)的增加而逐漸減小.

2 實驗與結(jié)果分析

2.1 溫度場的評價標(biāo)準(zhǔn)

溫度場重建的質(zhì)量好壞可以通過均方根誤差(Root Mean Square Error, RMSE)來衡量， RMSE越小說明重建效果越好. 用峰值誤差(Peak Value Error, PVE)來評價溫度場重建峰值的誤差， PVE越小說明重建的峰值越好.

(18)

(19)

2.2 仿真實驗

在本文的模擬仿真實驗中， 采用滿足高斯分布的三峰溫度場函數(shù)

0

(20a)

t(x,y)=0 其他.

(20b)

如圖 2 所示， 設(shè)網(wǎng)格數(shù)為30×30， 重建被測場區(qū)域的實際長度為30 cm×30 cm. 模擬環(huán)境溫度為0 ℃， 實際上三個峰值的溫度值分別為150, 200和250 ℃.

圖 2 三峰高斯模擬溫度場Fig.2 Three peak Gauss simulation of temperature field

在本文中， 采樣以0°～180°的范圍內(nèi)的22.5°的間隔進(jìn)行， 并且獲得8個方向上的投影數(shù)據(jù)， 并且每個方向上的樣本數(shù)量為40. 初始值為0， ART迭代算法中松弛因子λ=0.25， 并且在20次迭代后得到折射率的分布. 溫度t和折射率n之間的關(guān)系可以通過G-D(Gladston-Dale)公式來確定.

(21)

為了研究有噪聲條件下本文算法的重建質(zhì)量， 對帶有白噪聲的偏折角投影數(shù)據(jù)進(jìn)行了仿真. 圖 3 為兩種算法的重建結(jié)果， 圖 4 為兩種算法在投影數(shù)據(jù)疊加了白噪聲后的重建結(jié)果. 表 1 為兩種算法重建后的誤差.

圖 3 無噪聲的三峰溫度場重建結(jié)果Fig.3 Three peak temperature field reconstructed without noise

圖 4 有噪聲的三峰溫度場重建結(jié)果Fig.4 Three peak temperature field reconstructed with noise

表 1 兩種算法的重建誤差

從圖 2 和圖 3 的重建結(jié)果可以看出， 在8個投影方向的極度不完全數(shù)據(jù)的條件下， 不論在有噪聲的情況下還是在沒有噪聲的情況下， 都可以看出本文算法重建出的表面質(zhì)量明顯都優(yōu)于ART-TV算法. 而且本文算法可以更有效地重構(gòu)峰值. 從表1的重建誤差中可以看出， 迭代相同的次數(shù)， 不論是否存在噪聲， 本文算法的RMSE更小， 說明本算法的收斂速度更快， 誤差更小. 本文算法的誤差PVE比ART-TV算法要小， 說明本文的算法可以更好、 更有效地重構(gòu)出溫度場的峰值. 從有噪聲的誤差中可以看出本文算法的誤差RMSE和PVE都比ART-TV算法要小， 說明本文算法的抗噪聲性能更好一些.

2.3 液化石油氣火焰溫度場的重建

為了驗證本算法的有效性， 將該算法用于液化石油氣火焰溫度場的重建. 實驗數(shù)據(jù)由青島科技大學(xué)流場光測實驗室提供， 用莫爾偏折儀來獲取多方向的莫爾條紋圖， 實驗裝置見文獻(xiàn)[4].

圖 5 6個不同方向的莫爾條紋圖Fig.5 Morie deflectograms at six view angles

其中激光波長632.8 nm， 使用兩個完全相同的Ronchi光柵， 莫爾條紋圖成像在屏上， 然后再用CCD采集讀入計算機.

在重建過程中將重建區(qū)域劃分成30×30的網(wǎng)格， 使用本文提出的算法， 參數(shù)松弛因子取值為λ=0.55， 環(huán)境溫度為20°， 使用從到等間隔的6個投影方向的投影數(shù)據(jù)重建， 每個投影方向下的采樣數(shù)為150， 用本文改進(jìn)的算法迭代20次后， 得到溫度場的分布圖， 如圖 6 所示.

圖 6 重建溫度場Fig.6 Reconstructed temperature field

3 結(jié) 論

本文將偏折角轉(zhuǎn)化的迭代算法與TV范數(shù)相結(jié)合， 在傳統(tǒng)ART-TV算法的基礎(chǔ)上， 對梯度下降法的步長稍作調(diào)整， 采用一致性控制步長， 使梯度下降. 通過對三峰溫度場的仿真實驗， 通過重建結(jié)果可以看出， 本文的算法重建的溫度場平滑度優(yōu)于ART-TV算法， 在峰值的構(gòu)建方面表現(xiàn)出優(yōu)越性， 并且可以更好地抑制噪聲. 實驗結(jié)果表明， 本文的算法對于極少投影方向的溫度場的重建是有效的.