陳曉艷，房曉東

(天津科技大學(xué)電子信息與自動(dòng)化學(xué)院，天津 300222)

一種新的正則化圖像重建算法及參數(shù)優(yōu)化

陳曉艷，房曉東

(天津科技大學(xué)電子信息與自動(dòng)化學(xué)院，天津 300222)

針對(duì)電學(xué)層析成像技術(shù)的反問(wèn)題求解，提出一種新的正則化圖像重建算法及正則參數(shù)選擇方法．采用最小二乘擬合法構(gòu)建正則參數(shù)與評(píng)價(jià)參數(shù)的Matlab數(shù)學(xué)模型，取極值獲得正則參數(shù)值；利用COMSOL Multiphysics建立仿真場(chǎng)域，最后利用Matlab求解反問(wèn)題，重建圖像．通過(guò)對(duì)重建圖像及評(píng)價(jià)參數(shù)的對(duì)比分析，驗(yàn)證了方法的有效性．結(jié)果表明，所用方法及正則參數(shù)能夠有效改善反問(wèn)題的病態(tài)性，使反問(wèn)題的求解質(zhì)量最優(yōu)．

正則化；反問(wèn)題；圖像重建

電學(xué)層析成像技術(shù)(electrical tomography，ET)是過(guò)程層析成像技術(shù)(process tomography，PT)的一類，已被廣泛應(yīng)用于醫(yī)學(xué)診斷、工業(yè)測(cè)量等諸多領(lǐng)域．它在被測(cè)對(duì)象上按一定規(guī)律安放若干個(gè)電極，以非接觸或非侵入的方式測(cè)量被測(cè)對(duì)象的邊界電信息，然后利用圖像重建算法重建出被測(cè)對(duì)象內(nèi)部的電特性分布，實(shí)現(xiàn)對(duì)被測(cè)對(duì)象的可視化測(cè)量[1–2]．ET的反問(wèn)題求解具有高度的病態(tài)性，即使在測(cè)量數(shù)據(jù)中有很小的測(cè)量誤差，也會(huì)給反問(wèn)題的求解帶來(lái)災(zāi)難性的影響[3]．如何采取有效手段改善反問(wèn)題的病態(tài)性，從而獲得符合條件的滿意解，國(guó)內(nèi)外學(xué)者對(duì)此進(jìn)行了長(zhǎng)期的探索與研究．

正則化方法是解決反問(wèn)題病態(tài)性必不可少的方法，從某種角度，它可被視為保持解的存在性、唯一性及穩(wěn)定性的一個(gè)約束[4]．在20世紀(jì)60年代中期，蘇聯(lián)科學(xué)院院士Tikhonov最先提出求解不適定問(wèn)題的正則化方法[5]，為改善反問(wèn)題的病態(tài)性提供了解決的依據(jù)．2004年，Daubechies等[6]首次給出了稀疏正則化方法與迭代伸縮算法的理論分析，并指出其求解具有稀疏性的線性反問(wèn)題的有效性．2011年，彭黎輝[7]用Tikhonov正則化Landweber迭代方法對(duì)病態(tài)反問(wèn)題求解進(jìn)行了研究．正則化方法的求解質(zhì)量依賴于正則參數(shù)的選擇，目前主要還是依靠經(jīng)驗(yàn)來(lái)選取[8]．本文提出一種新的正則化算法，采用對(duì)角陣Q

代替單位陣I以施加阻尼作用，利用最小二乘擬合建立了正則參數(shù)和評(píng)價(jià)參數(shù)的數(shù)學(xué)模型，取極值獲得正則參數(shù)值，并優(yōu)化獲得最佳的正則參數(shù)，最后，通過(guò)重建圖像質(zhì)量和評(píng)價(jià)參數(shù)兩方面對(duì)方法的有效性進(jìn)行了驗(yàn)證．

1 算法原理

式中：U為被測(cè)對(duì)象的邊界電壓分布；J為靈敏度矩陣；σ為被測(cè)對(duì)象內(nèi)部的電特性分布矩陣(介電常數(shù)、電導(dǎo)率、磁導(dǎo)率)．ET反問(wèn)題求解，即已知被測(cè)對(duì)象的邊界電壓分布U及靈敏度矩陣J，求解對(duì)象內(nèi)部的電特性分布σ．若采用最小二乘法求解，可表示為

ET問(wèn)題的數(shù)學(xué)模型可表示為

令P=JTJ，由于矩陣P具有很高的條件數(shù)，導(dǎo)致求逆運(yùn)算中噪聲信號(hào)被放大，造成解的不穩(wěn)定．因此必須采取正則化手段，保證解的存在性、唯一性及穩(wěn)定性．

Tikhonov正則化算法是取與P同階次的單位陣I作為阻尼項(xiàng)，對(duì)P施加阻尼作用之后，經(jīng)過(guò)求逆運(yùn)算，從而進(jìn)行反問(wèn)題求解：

由于P=JTJ，J為任意矩陣時(shí)，P為對(duì)稱方陣，取P的每一行向量或列向量的最大元素均可構(gòu)成對(duì)角陣Q(Q的主對(duì)角元素非0，其他元素為0)．利用Q對(duì)P施加阻尼作用，從而獲得新的正則化算法：

在式(3)—式(4)中，α 為正則參數(shù)，通過(guò)調(diào)整α可改變阻尼作用的強(qiáng)弱．由于Q與P量級(jí)相同，在選取α 時(shí)，一般只需令其為0～1之間．相比于其他正則化算法，α 的范圍小，便于選擇，因而容易得到較為理想的求解結(jié)果．

2 評(píng)價(jià)參數(shù)

顯然，不同正則參數(shù)α 的選取，影響阻尼作用的強(qiáng)弱，會(huì)導(dǎo)致反問(wèn)題的求解結(jié)果σ 的不同．選取3個(gè)參數(shù)作為正則參數(shù)α 對(duì)σ 影響的評(píng)價(jià)參數(shù)．

2.1 誤差總和

定義誤差總和

式中：Ne為重建圖像的單元數(shù)；σi為重建電特性分布；σi*為真實(shí)電特性分布．E越小反問(wèn)題求解質(zhì)量越好．

2.2 相關(guān)系數(shù)

定義相關(guān)系數(shù)

2.3 結(jié)構(gòu)相似度

定義結(jié)構(gòu)相似度

3 正則參數(shù)優(yōu)化

由于正則參數(shù)α 會(huì)影響ET反問(wèn)題的最終結(jié)果，3個(gè)評(píng)價(jià)參數(shù)E、r、S也有差異．根據(jù)式(4)選取0～1范圍內(nèi)的α 進(jìn)行ET反問(wèn)題求解，相應(yīng)結(jié)果見(jiàn)表1．

由表1數(shù)據(jù)可知：當(dāng)α的值由0到1逐漸增大時(shí)，誤差總和E先減小后增大，相關(guān)系數(shù)r、結(jié)構(gòu)相似

度S先增大后減小，這表明反問(wèn)題求解質(zhì)量隨α值的增加先提高后降低．

3.1 模型建立

為得到理想的α值，采用最小二乘擬合的方法，以y＝f(x)的形式描述表1中各評(píng)價(jià)參數(shù)值隨α變化的規(guī)律．利用Matlab曲線擬合工具箱，取α＝0.1～1.0，步長(zhǎng)為0.01，采用四階擬合建立數(shù)學(xué)模型：

采用歸一化均方誤差(normalized mean square error)表征擬合精度，其表達(dá)式為如

3.2 α 求解

通常，α 取值越大，抗干擾性越強(qiáng)，問(wèn)題求解的穩(wěn)定性越高；但是α 過(guò)大會(huì)影響模型精度，覆蓋模型的有效信息．簡(jiǎn)言之，如果一個(gè)模型中α取值較小并且模型穩(wěn)定性好，則效果最好．分析表2中數(shù)據(jù)，當(dāng)α取0.403,1時(shí)，3個(gè)評(píng)價(jià)參數(shù)均較為理想．

4 圖像重建

利用有限元仿真軟件COMSOL Multiphysics 3.5,a對(duì)場(chǎng)域建模．圓形場(chǎng)域半徑為15,cm，場(chǎng)域表面安放16個(gè)電極，電極長(zhǎng)1,cm，電極寬度的占空比為20%，背景電導(dǎo)率為0.01,S/m，目標(biāo)電導(dǎo)率為1,S/m，激勵(lì)電極選取鈦電極(Titanium beta-21S)，激勵(lì)電流密度為1.392,mA．采用電流激勵(lì)–電壓測(cè)量、相鄰激勵(lì)–相鄰測(cè)量的模式．為了體現(xiàn)實(shí)際情況下的干擾，使結(jié)果更具一般性，在仿真數(shù)據(jù)中加入信噪比為42,dB的噪聲，根據(jù)不同的α值，重建圖像．場(chǎng)域原始模型和重建圖像如圖2所示．圖2(a)中白色區(qū)域?yàn)槟繕?biāo)區(qū)域，黑色區(qū)域?yàn)楸尘皡^(qū)域．重建圖像的灰度代表對(duì)應(yīng)的電導(dǎo)率數(shù)值．

根據(jù)圖2的重建圖像結(jié)果可知：當(dāng)α＝0，即不進(jìn)行正則化作用時(shí)，成像效果很不理想，場(chǎng)域中無(wú)法識(shí)別目標(biāo)的位置；當(dāng)α＝0.01時(shí)，成像效果有所改善，能夠分辨出目標(biāo)的大概位置，但偽影較多，目標(biāo)對(duì)比度較低；當(dāng)α＝0.403,1時(shí)，成像效果進(jìn)一步改善，偽影減少，目標(biāo)對(duì)比度增加，目標(biāo)大小與真實(shí)大小基本吻合，參照?qǐng)D像右側(cè)的色度條，此時(shí)目標(biāo)區(qū)域的灰度、背景區(qū)域的灰度所對(duì)應(yīng)的電導(dǎo)率值與模型的設(shè)定值0.01,S/m、1,S/m最為接近；當(dāng)α＝10時(shí)，偽影增加，目標(biāo)邊界模糊，成像質(zhì)量下降．表3給出了圖2對(duì)應(yīng)的評(píng)價(jià)參數(shù)值．

由表3可知：當(dāng)α＝0.01時(shí)，與α＝0(即不采用正則化方法)相比，各參數(shù)值趨于改善，這是由于進(jìn)行了正則化，在阻尼項(xiàng)的作用下，增強(qiáng)了反問(wèn)題求解時(shí)矩陣求逆的穩(wěn)定性；當(dāng)α＝0.403,1時(shí)，各參數(shù)值進(jìn)一步改善，求解質(zhì)量提高，這是由于此時(shí)的正則參數(shù)增加，正則化效果加強(qiáng)，模型穩(wěn)定性進(jìn)一步增強(qiáng)；當(dāng)α＝10時(shí)，各參數(shù)的質(zhì)量下降，求解質(zhì)量下降，這是由于當(dāng)正則參數(shù)過(guò)大時(shí)，過(guò)于強(qiáng)調(diào)模型的穩(wěn)定性和抗擾性，忽視模型自身的有效信息，導(dǎo)致了圖形偽影顯著，邊界模糊．表3的參數(shù)數(shù)據(jù)與圖2的重建圖像結(jié)果一致，且驗(yàn)證了算法具有一定的抗干擾能力．

5 結(jié) 語(yǔ)

本文針對(duì)二維ET反問(wèn)題求解，提出一種新的正則化圖像重建算法，采用最小二乘擬合法建立數(shù)學(xué)模型，取極值獲取了正則參數(shù)值．通過(guò)求解，從重建圖像效果及評(píng)價(jià)參數(shù)兩方面證實(shí)了方法的有效性．

不同于一般正則化算法采用單位陣作為阻尼項(xiàng)，本研究提出的正則化圖像重建算法，結(jié)合問(wèn)題自身的靈敏度矩陣，采用矩陣P的每一行向量或列向量的最大元素構(gòu)成對(duì)角陣Q進(jìn)行阻尼約束，從而使阻尼項(xiàng)Q較I更加貼近模型本身，與模型的相關(guān)性更強(qiáng)，與一般的正則化方法相比，大大縮小了α的選擇區(qū)間；采用最小二乘擬合建立數(shù)學(xué)模型，利用數(shù)學(xué)方法，通過(guò)取極值獲得正則參數(shù)α的值，使其在取值范圍內(nèi)最優(yōu)，與一般的正則化方法相比，避免了憑經(jīng)驗(yàn)選取α的盲目性及不確定性．這將為今后進(jìn)一步研究反問(wèn)題求解奠定基礎(chǔ)．

本研究的建模方法采用基于Matlab的最小二乘擬合法，未討論其他形式的擬合方法，今后將圍繞不同的方法建立數(shù)學(xué)模型，進(jìn)行對(duì)比研究．

［1］ Yue S H，Wu T，Cui L J，et al. Clustering mechanism for electric tomography imaging[J]. Science China：Information Sciences，2012，55(12)：2849–2864.

［2］ 何泳成. 電學(xué)層析成像圖像重建算法研究及軟件系統(tǒng)設(shè)計(jì)[D]. 天津：天津大學(xué)，2010.

［3］ Colton D，Kress R. Inverse acoustic and electromagnetic scattering theory[M]. 2nd ed. Berlin：Springer，1998.

［4］ Polydorides N. Image reconstruction algorithms for soft

field tomography[D]. Manchester，UK：University of Manchester Institute of Science and Technology，2002.［5］ 呂琪. 不適定問(wèn)題的迭代正則化方法研究[D]. 武漢：武漢理工大學(xué)，2012.

［6］ Daubechies I，Defrise M，De Mol C. An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[J]. Communications on Pure and Applied Mathematics，2004，57(11)：1413–1457.

［7］ 彭黎輝. 用于病態(tài)反問(wèn)題求解的正則化方法算例分析[J]. 云南民族大學(xué)學(xué)報(bào)：自然科學(xué)版，2011，20(2)：79–85.

［8］ 周旭勝，寇戈. 基于MATLAB的EIT圖像重構(gòu)算法研究[D]. 南京：南京理工大學(xué)，2010.

［9］ 吳克堅(jiān). 有限元近似誤差對(duì)EIT正逆問(wèn)題影響的定量研究及算子分解方法的應(yīng)用[D]. 西安：第四軍醫(yī)大學(xué)，2012.

責(zé)任編輯：常濤

A New Regularization Image Reconstruction Algorithm and its Parameter Optimization

CHEN Xiaoyan，F(xiàn)ANG Xiaodong
(College of Electronic Information and Automation，Tianjin University of Science & Technology，Tianjin 300222，China)

Focused on solving the inverse-problem of electrical tomography，a new regularization image reconstruction algorithm and a method to select the regularization parameters were proposed. A Matlab mathematic model of regularization parameter and evaluation parameter was constructed by the least square fit method．The values of regularization parameter can be obtained by taking the extreme．COMSOL Multiphysics was adopted to construct the simulation field，and Matlab was used to solve the inverse-problem and reconstruct the image. The effectiveness of the method was verified by contrasting the reconstructed image and the evaluation parameters. The result indicates that the method and regularization parameters can improve the ill-pose of the inverse-problem effectively and make the solution to the inverse problem optimal.

regularization；inverse-problem；image reconstruction

TP391.9

A

1672-6510(2014)06-0074-04

10.13364/j.issn.1672-6510.2014.06.014

2014–04–10；

2014–07–08

國(guó)家自然科學(xué)基金資助項(xiàng)目(61301246)；天津市自然科學(xué)基金資助項(xiàng)目(12JCYBJC19300)

陳曉艷（1973—），女，四川成都人，教授，cxywxr@tust.edu.cn.