基于模糊算子的電學(xué)層析成像算法

岳士弘，榮西拉，馬海濤，盧 劍

(天津大學(xué)電氣自動(dòng)化與信息工程學(xué)院，天津 300072)

電學(xué)層析(electrical tomography，ET)是一種新型、先進(jìn)、非侵入式的可視化測量技術(shù)，它具有成本低、響應(yīng)快、安全性高等優(yōu)點(diǎn)，在工業(yè)和醫(yī)學(xué)領(lǐng)域具有廣闊的應(yīng)用前景[1-2]．一個(gè)典型的ET 系統(tǒng)主要通過被測場域邊界電極陣列的輪流激勵(lì)和測量；測量值通過信號(hào)數(shù)據(jù)采集單元；最后上傳到上位機(jī)進(jìn)行圖像重建與顯示可視化．一個(gè)ET 系統(tǒng)結(jié)構(gòu)簡單易于實(shí)現(xiàn)．然而，由于ET 過程固有的“欠定性”和“軟場效應(yīng)”[3-4]，其成像的空間分辨率很低，難以滿足很多應(yīng)用中的基本需求．在大多數(shù)情況下，ET 成像算法依賴于一個(gè)基于目標(biāo)函數(shù)的優(yōu)化問題．從最先提出的核心優(yōu)化算法——高斯-牛頓(Gaussian-Newton，GN)法[5-7]、吉洪諾夫正則化(Tikhonov regularization，TR)[8]法到最近的基于Lp范數(shù)各種變形形式[9-11]以及Bregman 散度[12]，它們都有各自可應(yīng)用的范圍以及明顯的局限性．如何面向一個(gè)具體的應(yīng)用對象，采用有效的方法提高ET 的空間分辨率和穩(wěn)定性，幾十年來一直是該領(lǐng)域一個(gè)重點(diǎn)和難點(diǎn)問題．

針對上述問題，本文基于模糊算子[13]提出一個(gè)一般的ET 成像算法．它具有兩個(gè)優(yōu)勢：①具有近百個(gè)功能各異的模糊算子可供使用者選擇，特別是目前ET 中目標(biāo)函數(shù)常用的加權(quán)和形式只不過是模糊算子的一個(gè)特例；②一個(gè)易于實(shí)現(xiàn)的非傳統(tǒng)優(yōu)化的求解方法；比較而言，很多已有的ET 優(yōu)化算法往往要求可導(dǎo)或至少是連續(xù)的，面臨著如何優(yōu)化和實(shí)現(xiàn)的問題，而本文方法只要求目標(biāo)函數(shù)在可行域上每點(diǎn)有唯一的函數(shù)值就可以實(shí)現(xiàn)優(yōu)化；這種要求更加符合工程實(shí)際．最后，通過重建圖像質(zhì)量和評價(jià)參數(shù)對方法的有效性進(jìn)行了驗(yàn)證．

1 研究背景與相關(guān)工作

1.1 電學(xué)層析原理

ET 技術(shù)主要分為電阻層析成像[14]、電容層析成像[15]和電磁層析成像[16]．這3 類技術(shù)雖存在一定差異，其重構(gòu)原理是相似的．本文僅以二維電阻層析成像為例進(jìn)行說明．

將調(diào)查的場域采用有限元刨分和離散化，并將邊界測量與內(nèi)部場域關(guān)系進(jìn)行線性化后得

式中：m 為被測場域內(nèi)的獨(dú)立測量數(shù)；n 為像素點(diǎn)個(gè)數(shù)；Jm×n為雅可比矩陣，表示測量電極對每個(gè)像素點(diǎn)單位電導(dǎo)率增量的響應(yīng)程度或靈敏度；σ 為被測場域的電導(dǎo)率分布．記U為邊界測量構(gòu)成的向量，則歸一化式(1)后得到

該方法求解速度快并且易于實(shí)現(xiàn)．但是方程(2)對應(yīng)的逆問題求解本質(zhì)上是非線性、欠定和病態(tài)的，因此LBP 算法很難獲得準(zhǔn)確的解．大多數(shù)ET 算法把問題轉(zhuǎn)化為以下優(yōu)化問題，即

式中||·||為任意范數(shù)．過去大多數(shù)研究集中在使用不同形式的范數(shù)對式(4)進(jìn)行求解；但另一種改變是式(4)左邊增加限制項(xiàng)，其中最典型的是TR 算法，它具有下列形式，即

式中λ 為正則化系數(shù)．

這些方法對于具體應(yīng)用往往有特有的優(yōu)勢．但是，這些方法的使用往往帶有一定的假設(shè)且泛化能力不強(qiáng)難以移植[18]．此外，對于任意一個(gè)目標(biāo)函數(shù)如何求解也是關(guān)鍵問題[19]，克服這些問題是當(dāng)前ET 研究領(lǐng)域的重點(diǎn)和難點(diǎn)．

1.2 模糊算子

本文模糊算子是指基于模糊關(guān)系的合成運(yùn)算[20].模糊算子種類豐富而多樣，在模糊推理過程中，算子選取的不同，會(huì)導(dǎo)致推理方法的可適用范圍和推理結(jié)果的不同[21]．對于任意a，b∈[0，1]，模糊算子T(a，b)(常稱為T模)必須滿足3 個(gè)條件：

常用的模糊算子有以下幾種：

上述“max”算子就是在兩個(gè)變量之間取較大的一個(gè)，這個(gè)算子應(yīng)用最為簡單和方便，但容易丟失信息．相比于max 算子，式(7)中其他算子能更多地保留運(yùn)算前的信息．如T2和T4實(shí)際上是求兩個(gè)向量的夾角正弦和正切值等；它們對于模糊推理和模糊運(yùn)算能夠起到不同的作用．

2 基于模糊三角算子的圖像重建算法

2.1 基本動(dòng)機(jī)和模型

本文提出算法的核心思想是根據(jù)測量值找到一個(gè)合適的模糊算子，通過這個(gè)模糊算子代替式(2)中的求和運(yùn)算，實(shí)現(xiàn)下述一般化轉(zhuǎn)換，即

式中：角標(biāo)k和j對應(yīng)第k個(gè)方程(測量)以及第j個(gè)變量，[skj]∈S；“·”表示對于任何兩項(xiàng)使用模糊三角算子做運(yùn)算．使用式(8)的主要?jiǎng)訖C(jī)包括以下幾個(gè)方面．

(1) 由于軟場效應(yīng)和非線性，S 中每個(gè)元素skj所反映的任意第k個(gè)測量值uk對于第j個(gè)像素σj具有高度不確定性，從而看成模糊量是必要和可行的；將S 中的值離散和模糊化之后，一定程度上能反映靈敏度系數(shù)本身不準(zhǔn)確性以及測量值中噪聲的影響，具有較好的可解釋性．

(2) 大量具有不同可適用范圍和功能的模糊算子為式(8)的表達(dá)和拓展提供了豐富的選擇．實(shí)際上，式(2)相當(dāng)于模糊算子中取T3形式的一個(gè)特殊情況；相對于一般的工程需求和復(fù)雜工況，有潛力實(shí)現(xiàn)更好和更有效的計(jì)算．

(3) 式(8)中取模糊算子構(gòu)造的任何形式目標(biāo)函數(shù)都必須通過優(yōu)化才能實(shí)現(xiàn)電學(xué)層析過程．在第2.2節(jié)將使用一個(gè)一般的優(yōu)化方法使式(8)的求解和實(shí)現(xiàn)成為可能．

(4) 工程中很多先驗(yàn)信息是用語言規(guī)則表示的.模糊系統(tǒng)正是基于經(jīng)驗(yàn)和規(guī)則進(jìn)行推理和分析的，能夠較好地將先驗(yàn)信息和工程實(shí)踐經(jīng)驗(yàn)導(dǎo)入模糊系統(tǒng)，從而抑制噪聲并提高系統(tǒng)的魯棒性．特別是對于ET系統(tǒng)面臨的具有不一致、不確定和不完整特點(diǎn)的工程問題，精確模型及其優(yōu)化往往耗費(fèi)極大，而模糊模型能夠提供更好的選擇．

2.2 實(shí)現(xiàn)算法

假設(shè)ET 的邊界測量值已知且靈敏度矩陣S 已經(jīng)預(yù)先確定．基于一般的模糊算子對應(yīng)的ET 任務(wù)是如何求解最優(yōu)電導(dǎo)率的分布，使得式(8)對應(yīng)的目標(biāo)函數(shù)誤差最小化．

假設(shè)過程的m個(gè)測量值uk的被歸一化為一組離散值a1，a2，…，am，滿足0＜a1＜a2＜…＜am-1＜am＜1，那么求解式(8)右邊的目標(biāo)函數(shù)等價(jià)于近似求解一組σ1，σ2，…，σn的值，使得

式中b1，b2，…，bm是在0，a1，a2，…，am之間的任意插值，滿足0＜b1＜a1＜b2＜…＜am＜bm． 對于滿足此條件的ai和bi和任意的模糊算子 “·”，必然存在一個(gè)嚴(yán)格增函數(shù)f[22]，滿足f(0)＝0，且

式(10)的意義在于：對于式(9)中任意模糊算子的運(yùn)算結(jié)果，都可以通過一個(gè)一元增函數(shù)f及其簡單的加法運(yùn)算代替來計(jì)算．

如果f(br)的值均為已知，那么為了求出σj的值，必須根據(jù)式(10)求解未知量f(skjσj)，j＝1，2，…，n．這類不等式有很多求解方法，本文使用一個(gè)易于實(shí)現(xiàn)和操作的糾錯(cuò)算法[23]來實(shí)現(xiàn)．

式中[s/2]、[(s＋1)/2]中的符號(hào)[·]表示向下取整，即不超過s/2 以及(s＋1)/2 的最大整數(shù)．

如果設(shè)Qk=(q1, q2, …,q2m)，Rk=(r1, r2, …, r2m)，W =(u1, u2, u3, …,u2m)，Es表示第s列為1 其余列均為0 的1×2m單位矩陣，式(12)、式(16)和式(17)的求解可以被寫為

求式(18)的一個(gè)一般解W*的迭代步驟為

式中：μ 是一個(gè)權(quán)重系數(shù)，若無先驗(yàn)信息可取為1；Ai是將Qk、Rk、Es按照迭代循環(huán)排列后的向量；||Ai||是取定的矩陣標(biāo)準(zhǔn)范數(shù)．通常序列Ai的迭代可以是無限的，但該算法必然會(huì)在有限的次數(shù)內(nèi)收斂．事實(shí)上，式(19)的迭代次數(shù)不超過上限[23]

式中f-1為單調(diào)函數(shù)f的逆算子．為了最終求出各個(gè)像素的σ分布值，必須將單調(diào)函數(shù)f具體表達(dá)出來．

考慮到ET 過程的近似性，可以與已知的測量值uk構(gòu)造數(shù)據(jù)點(diǎn)列(uk，f(uk))進(jìn)行樣條插值[24]，i＝1，2，…，m，設(shè)獲得的插值函數(shù)為F(x)，而作為一個(gè)單調(diào)單元函數(shù)，它的逆函數(shù)F-1(x)值即為F(x)．它也是

f--1(x)的近似表達(dá)式，這里取為F(x)的倒數(shù)．因此，F(xiàn)-1(x)的表達(dá)式為

最后，各個(gè)像素的σ的分布值為

由于靈敏度系數(shù)skj的值是完全確定的，式(25)作為選擇ET 中關(guān)鍵的并算子一般通式．

當(dāng)一組模糊算子可以適用和選擇時(shí)，由式(23)可以定義不同的結(jié)果，為了檢測哪個(gè)算子是最適合的結(jié)果，可以通過下面的式子求解，即

使得式(26)獲得最小值的模糊三角算子為最優(yōu)解．

最后，當(dāng)任意一個(gè)模糊三角算子給定后，基于上述各個(gè)公式，可以得出以下關(guān)于電學(xué)層析的新優(yōu)化算法基本實(shí)施步驟．

步驟1獲得邊界電壓測量數(shù)據(jù)并歸一化后升序排列為u1，u2，…，um；

步驟2確定空場的靈敏度系數(shù)skj，k＝1，2，…，m；j＝1，2，…，n；

步驟3根據(jù)式(16)和式(17)計(jì)算值；

步驟4根據(jù)式(19)求解不等式(18)得到一個(gè)最優(yōu)解W*；

步驟 5根據(jù)式(22)求解f(skjσ)和f(bj)；

步驟 6根據(jù)式(25)計(jì)算電導(dǎo)率分布．

3 仿真實(shí)驗(yàn)和結(jié)果

利用有限元仿真軟件COMSOL 對場域建模，圓形場域半徑為10 cm，場域表面分別安放16 個(gè)電極，電極長為1 cm，背景電導(dǎo)率為1 S·m，目標(biāo)電導(dǎo)率為0.01 S·m．原始模型和重建圖像如表1 所示．采用常用的相鄰激勵(lì)和相鄰測量方法，可獲得的獨(dú)立樣本數(shù)為208 個(gè)，為此場域刨分為812 個(gè)單元．

表1 中第1 列對應(yīng)原型，其中紅色像素構(gòu)成目標(biāo)，藍(lán)色區(qū)域像素對應(yīng)背景．這些原型又兩兩分為3組：1 組為兩個(gè)目標(biāo)離散的模型；2 組為人體胸部截面模型；3 組為兩個(gè)目標(biāo)連續(xù)的模型；這些原型具有不同的代表性．第2 到第6 列為式(10)中采用模糊算子T2、T4、T5以及使用TR 得到ET 圖像；其中T3相當(dāng)于使用了LBP 算法；而TR 算法的結(jié)果受正則化參數(shù)λ 的影響．為此，實(shí)驗(yàn)中遍歷了參數(shù)λ在一定范圍內(nèi)的所有可能值，并取了對應(yīng)值下最佳空間分辨率的重構(gòu)結(jié)果進(jìn)行比較．

本實(shí)驗(yàn)實(shí)施過程中做以下設(shè)定．

(1) 歸一化處理．對靈敏度系數(shù)矩陣S 與測量值U 相乘作為整體進(jìn)行歸一化，即先把靈敏度系數(shù)矩陣的每一行都乘以該行對應(yīng)的最大值，然后再進(jìn)行歸一化．由于作為求解依據(jù)的測量值不是完全準(zhǔn)確的而是存在誤差，對于模擬氣泡分布比較特殊的樣本，求解得到的W 是有可能出現(xiàn)個(gè)別元素有異常值，因此需要足夠數(shù)量和不同分布的樣本來避免這些偶然因素．

(2) 離散化的尺度．為使求解不等式不過于復(fù)雜同時(shí)照顧實(shí)際問題中的精度要求，aj和bj各使用30 個(gè)不同的但等間隔的取值，即aj＝1/30，1/15，1/10，2/15，…，1；bj＝1/60，1/20，1/12，7/60，…，59/60．考慮算法的邊值的范圍，要求aj和bj的每個(gè)值不能等于 0 或等于 1，最后一個(gè) a 值可改為0.999 9．這相當(dāng)于輸出的所有灰度值可能為30 個(gè)不同的離散值，即把連續(xù)的灰度值轉(zhuǎn)化為30 個(gè)不同的灰度等級，這對于僅有812 個(gè)像素的ET 應(yīng)用于實(shí)際問題是足夠的．

(3)W的初值為全1 向量．雖然理論上求解W的糾錯(cuò)算法必然會(huì)收斂，但實(shí)際收斂的速度因不同樣本而異可能會(huì)十分緩慢．為保證效率，規(guī)定序列Aj經(jīng)過1 000 個(gè)完整周期之后，必須終止．同時(shí)，給每一次W的修正量取不同的修正權(quán)重μ值，隨著序列Aj的循環(huán)次數(shù)不斷增加，每次檢測收斂條件時(shí)，W的修正量不斷減少．第i 次循環(huán)的權(quán)重μi為100/(i＋99)．例如序列Aj已經(jīng)在最后一個(gè)周期時(shí)，W的修正能力只有W1的修正能力的1/10．

(4) 閾值設(shè)定．為了避免異常值，先計(jì)算前50個(gè)樣本求解得到的W向量的每一個(gè)元素的平均值和標(biāo)準(zhǔn)差E，此后對于每一個(gè)新樣本得到的新的W向量，計(jì)算每一個(gè)元素和該元素平均值的偏差，如果超過5E 則視為無效值．

表1 顯示了仿真實(shí)驗(yàn)所測試的6 個(gè)原型在不同的模糊算子下的成像結(jié)果，按照3 組不同的類型具體說明如下．

(1) 對于一組中目標(biāo)離散分布的模型，可以直觀地觀察到T2和T5目標(biāo)成像更加接近于原型，T4次之而T3(LBP)成像結(jié)果總體是最差的，這是由于LBP不包含任何迭代過程，成像結(jié)果最差，因此，T2和T5更具有分辨離散分布目標(biāo)的能力．

表1 使用模糊算子T2、T3、T4、T5 的ET成像對比Tab.1 Comparison of ET images by using fuzzy operators T2，T3，T4，and T5

對于TR 算法而言，其重構(gòu)圖像中的目標(biāo)偽跡很小，不僅目標(biāo)位置正確而且邊緣清楚．

(2) 對于2 組中肺部界面模型而言，由于設(shè)置了3 種電導(dǎo)率分布的目標(biāo)，它們分別對應(yīng)肺部組織(包含2 個(gè)和4 個(gè)結(jié)節(jié))、心臟和外部骨骼，導(dǎo)致ET 成像結(jié)果都不太理想；各個(gè)組織的邊界都不夠清楚，幾乎無法發(fā)現(xiàn)心臟；但是，模擬醫(yī)學(xué)診斷過程的結(jié)節(jié)能夠一定程度上呈現(xiàn)．

(3) 對于最后兩個(gè)目標(biāo)聯(lián)系的界面模型而言，分辨率是有邊界的清晰度決定；其中，第2 個(gè)模型對應(yīng)有兩層邊界；從ET 結(jié)果圖中可以看出，基于T3的算法實(shí)際輸出與原型比較更加接近，所有邊界基本上可以辨識(shí)出來．這說明T3(LBP)在重構(gòu)真實(shí)的邊界時(shí)是有優(yōu)勢的．比較而言．基于模糊算子T4的結(jié)果空間分辨率比T3有所提高，但仍然存在著大量偽跡；TR 重構(gòu)的界面是最不清楚的．

上述4 個(gè)模糊算子對應(yīng)的空間分辨率可以用相關(guān)系數(shù)λR(correlation coefficient)和RE(relative error，相對誤差er)指標(biāo)來進(jìn)一步定量說明和評價(jià)，分別定義為

式中：G和G′分別表示重建圖像和實(shí)際圖像的灰度值向量；σi和σi*分別為σ 和σ*的第i 個(gè)像素；σ0和σ0*分別為σ 和σ*的平均值；n 為表示像素總數(shù)．er越小表明ET 算法重構(gòu)的圖像越接近真實(shí)圖像，從而成像分辨率越高；er值在極端情況下不小于0．λ指標(biāo)越大表示重建圖像的灰度與實(shí)際灰度相關(guān)程度越大，重建圖像的質(zhì)量越高．

表2 顯示了3 個(gè)算法的λ值與er值．顯然，模糊算子的λR值與er值顯示結(jié)果與它們的可視化圖像是一致的．依據(jù)不同算子和模型，er值最大可以減小0.146 2，且平均減小0.093 1．總體而言，當(dāng)目標(biāo)離散時(shí)，T3對應(yīng)結(jié)果是最劣的．而T5目標(biāo)函數(shù)比模糊算子的目標(biāo)函數(shù)值大近一個(gè)數(shù)量級．當(dāng)目標(biāo)連續(xù)時(shí)，T3對應(yīng)結(jié)果是有優(yōu)勢的．對于像模擬人肺部斷面這樣的復(fù)雜分布的情況，幾種算子各有特點(diǎn)；這說明不同的模糊三角3 種有各自可適用的范圍．TR 算法由于帶有約束項(xiàng)，對于離散的小目標(biāo)辨識(shí)能力很高；但是對于連續(xù)的目標(biāo)邊緣(如兩種介質(zhì)邊界面)則幾乎失效，而T3(LBP)算法則正好相反．使用Matlab 中的Tic 和Toc 函數(shù)對6 個(gè)模型統(tǒng)計(jì)得：T2～T5的平均運(yùn)行時(shí)間基本是一致的，分別為0.201 s、0.223 s、0.230 s和0.219 s，TR 由于要求逆矩陣時(shí)間為0.332 s，略長于其他方法．

表2 使用模糊算子T2、T3、T4、T5 和TR的精度比較Tab.2 Comparison of accuracy by using fuzzy operators T2，T3，T4，T5，and TR

4 結(jié) 語

本文將模糊理論中的模糊算子運(yùn)用在ET 核心方程組的求解上．首先從ET 系統(tǒng)的欠定性和病態(tài)性出發(fā)，分析了ET 圖像重建中采用模糊算子代替矩陣乘法中的求和運(yùn)算的理論依據(jù)與合理性，然后針對ET 系統(tǒng)的特點(diǎn)，提出了一種可代替?zhèn)鹘y(tǒng)求和運(yùn)算的模糊算子的算法，并給出了采用模糊算子的ET 成像方法，最后通過仿真實(shí)驗(yàn)，證明了不同的模糊三角算子對于不同的可視化目標(biāo)有各自適用范圍，從而拓寬了ET 技術(shù)的可應(yīng)用范圍并提高了針對性．特別情況下，已有的ET 算法所常用的加和形式只不過是本文算法取特殊模糊算子的一個(gè)特例，證明本文算法具有一定泛化能力．

考慮到模糊算子的一般性和普遍性，本文僅為提高ET 成像的空間分辨率提供一條新的途徑或者說打開一個(gè)新的窗口；所使用的實(shí)驗(yàn)驗(yàn)證也是有限的.本文將來的研究聚焦在真實(shí)和復(fù)雜環(huán)境中，探討本文提出方法的可應(yīng)用范圍．