岳士弘,榮西拉,馬海濤,盧 劍
(天津大學(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)證.
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).
本文模糊算子是指基于模糊關(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)算能夠起到不同的作用.
本文提出算法的核心思想是根據(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)極大,而模糊模型能夠提供更好的選擇.
假設(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)率分布.
利用有限元仿真軟件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
本文將模糊理論中的模糊算子運(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)用范圍.