基于SOA優(yōu)化的電磁層析成像圖像重建研究

王亞?wèn)|，劉 澤，霍繼偉，趙鵬飛，王成飛，袁 偉

(1.北京交通大學(xué) 電子信息工程學(xué)院，北京 100044;2.中國(guó)鐵道科學(xué)研究院集團(tuán)有限公司 標(biāo)準(zhǔn)計(jì)量研究所，北京 100081)

0 引言

電磁層析成像(EMT,electromagnetic tomography)是一種近二、三十年發(fā)展起來(lái)的基于電磁感應(yīng)原理的過(guò)程層析成像技術(shù)[1]，其目的是研究具有電磁特性的物質(zhì)在空間的分布。典型的EMT系統(tǒng)主要由4個(gè)部分組成：傳感器陣列(即檢測(cè)激勵(lì)線圈)、激勵(lì)信號(hào)產(chǎn)生和分配單元、數(shù)據(jù)采集和調(diào)制單元以及圖像重建和信息獲取單元。通過(guò)對(duì)激勵(lì)線圈通入交變電流，在被測(cè)物場(chǎng)中產(chǎn)生激勵(lì)磁場(chǎng)；當(dāng)被測(cè)物場(chǎng)內(nèi)有具有導(dǎo)磁或?qū)щ姷奈矬w放入時(shí)，原有磁場(chǎng)分布將會(huì)改變，線圈之間的感性耦合也將隨之改變[2-3]。通過(guò)切換不同的激勵(lì)線圈，可獲得多個(gè)投影方向的檢測(cè)信息，再運(yùn)用合適的圖像重建算法就可以重建出放入物場(chǎng)中物體的電導(dǎo)率或者磁導(dǎo)率的分布圖，從而獲得放入物體的信息。由此可見(jiàn)，EMT具有非侵入、非接觸和快速性的特點(diǎn)，目前在生物醫(yī)學(xué)、無(wú)損檢測(cè)、鋼軌探傷、異物檢測(cè)和工業(yè)煉鋼等領(lǐng)域均已展開(kāi)相應(yīng)研究，具有很好的應(yīng)用以及發(fā)展前景[4-6]。

在EMT系統(tǒng)中，圖像重建算法是EMT技術(shù)研究中的重點(diǎn)和難點(diǎn)，圖像重建算法的好壞將直接影響最終成像的準(zhǔn)確性。目前的圖像重建算法主要有線性反投影(LBP)[7]、Tikhonov正則化[8]、Landweber迭代[9]等算法。在傳統(tǒng)的EMT系統(tǒng)中，邊界測(cè)量值的微小變化就會(huì)導(dǎo)致求解出的物場(chǎng)內(nèi)的介質(zhì)分布有很大變化，導(dǎo)致求解過(guò)程不穩(wěn)定，因此EMT技術(shù)本身具有病態(tài)性，嚴(yán)重影響了重建圖像的質(zhì)量。目前常采用截?cái)嗥娈愔邓惴ê蚑ikhonov正則化算法來(lái)改善靈敏度矩陣的病態(tài)性[10-11]，但都存在參數(shù)選取問(wèn)題。Tikhonov正則化算法主要依賴于正則化因子，但是正則化因子的選取較為復(fù)雜，在實(shí)際應(yīng)用中往往依據(jù)經(jīng)驗(yàn)進(jìn)行選取，具有較大的局限性。如果正則化參數(shù)選擇較小，原始的問(wèn)題就不能得到很好地近似，逆問(wèn)題的病態(tài)性的改進(jìn)效果不顯著；如果正則化參數(shù)選擇過(guò)大，雖然對(duì)誤差的敏感性會(huì)減少，但最終解會(huì)與正常值相差較大，達(dá)不到效果。

電磁層析成像的逆問(wèn)題的病態(tài)性的根源是靈敏度矩陣的條件數(shù)較大，本文通過(guò)利用人群搜索算法(SOA，seeker optimization algorithm)對(duì)物場(chǎng)為空時(shí)的靈敏度矩陣進(jìn)行優(yōu)化，能夠明顯降低其條件數(shù)。將優(yōu)化過(guò)后的靈敏度矩陣應(yīng)用于Landweber迭代算法進(jìn)行圖像重建，仿真實(shí)驗(yàn)結(jié)果表明，該方法能夠明顯提高圖像重建的效果以及精度。

1 理論基礎(chǔ)

1.1 電磁層析成像數(shù)學(xué)模型

電磁層析成像的研究包括正問(wèn)題以及逆問(wèn)題。正問(wèn)題是具有電磁特性的物質(zhì)在敏感場(chǎng)中的分布位置已知，求解其對(duì)激勵(lì)場(chǎng)的調(diào)制的結(jié)果；逆問(wèn)題則相反，是利用在正問(wèn)題中已經(jīng)得出的調(diào)制結(jié)果，來(lái)求解具有電磁特性的物質(zhì)的在物場(chǎng)的分布。在EMT中，其調(diào)制結(jié)果往往表現(xiàn)為邊界測(cè)量值，即線圈的感應(yīng)電壓。二者的關(guān)系可由圖1所示。

圖1 正問(wèn)題與逆問(wèn)題關(guān)系圖

EMT正問(wèn)題中的模型中線圈感應(yīng)電壓可用式(1)描述：

Udec=f(Be(x,y),Bo(x,y),μ(x,y),σ(x,y),ε(x,y))

(1)

式中,Udec為測(cè)量電壓，Be為激勵(lì)主磁場(chǎng)，Bo為由渦流產(chǎn)生的次級(jí)磁場(chǎng)，σ為電導(dǎo)率，μ為磁導(dǎo)率，ε為介電常數(shù)。

EMT的逆問(wèn)題就是利用已知的靈敏度矩陣S以及邊界測(cè)量值U來(lái)解出灰度矩陣G，再根據(jù)所求出的灰度矩陣G即可對(duì)靈敏場(chǎng)中被測(cè)物體進(jìn)行圖像重建。由于檢測(cè)線圈獲得的邊界測(cè)量值與物場(chǎng)空間之間的關(guān)系是非線性的，因此EMT的逆問(wèn)題可近似歸結(jié)為線性方程組式(2)的求解[12]。

U=S·G

(2)

式中,U∈RM為歸一化之后的邊界測(cè)量值矩陣，S∈RM×N為歸一化后的靈敏度矩陣，G∈RN為歸一化之后的灰度矩陣。M為獨(dú)立測(cè)量值的個(gè)數(shù)，N是逆問(wèn)題剖分單元數(shù)。

1.2 EMT靈敏度矩陣分析

在EMT逆問(wèn)題的成像計(jì)算中，靈敏度矩陣作為先驗(yàn)信息，反應(yīng)了特定的激勵(lì)檢測(cè)對(duì)于不同位置單元介質(zhì)變化的敏感程度，是圖像重建的基礎(chǔ)。在實(shí)際應(yīng)用中，獲取靈敏度矩陣主要有3種方法：模型擾動(dòng)法[13]、測(cè)量擾動(dòng)法[14]和場(chǎng)量提取法[15]。因模型擾動(dòng)法及測(cè)量擾動(dòng)法計(jì)算量大且較為繁瑣，因此在本文中，將采用簡(jiǎn)便易行，計(jì)算速度快的場(chǎng)量提取法。

采用場(chǎng)量提取法進(jìn)行靈敏度矩陣的求取并畫(huà)出典型激勵(lì)檢測(cè)組合的靈敏度分布如圖2所示。

圖2 典型靈敏度矩陣分布圖

1.3 人群搜索算法應(yīng)用分析

電磁層析成像中的圖像重建屬于第一類Fredholm積分方程問(wèn)題，其解具有不適定性，這主要表現(xiàn)在圖像重建的過(guò)程中靈敏度矩陣的病態(tài)性。條件數(shù)可以用來(lái)衡量靈敏度矩陣的病態(tài)程度，因此，可以采取一定的算法令靈敏度矩陣的條件數(shù)達(dá)到最小，從而改善其病態(tài)性，即轉(zhuǎn)化為一個(gè)優(yōu)化問(wèn)題。

人群搜索算法[16]能夠直接模擬人的智能搜索行為，立足傳統(tǒng)的直接搜索算法，以搜尋隊(duì)伍為種群，搜尋者位置為候選解，通過(guò)模擬人類搜尋“經(jīng)驗(yàn)梯度”和不確定推理進(jìn)行搜索方向和搜索步長(zhǎng)的更新，并不斷更新位置，完成對(duì)問(wèn)題的最優(yōu)求解。

人群搜索算法主要包含搜索步長(zhǎng)、搜索方向以及個(gè)體位置的更新[17]，為了簡(jiǎn)化模型，提高運(yùn)算效率，在本文中，對(duì)目標(biāo)函數(shù)值進(jìn)行排序時(shí)采用線性隸屬函數(shù)，定義如式(3)和式(4)所示：

(3)

Uij=Ui+(1-Ui)·rand,j=1,2,,D

(4)

式中,Ui是目標(biāo)函數(shù)i的隸屬度；Uij是j維搜索空間目標(biāo)函數(shù)值i的隸屬度；P為種群大?。籇為搜索空間維數(shù)；Ki為種群函數(shù)值按照降序排列之后的序列編號(hào)；rand()能夠生成0～1之間的隨機(jī)數(shù)，用來(lái)模仿人的隨機(jī)行為。

1.4 Landweber迭代算法

由于邊界測(cè)量值與敏感場(chǎng)中的電導(dǎo)率分布是非線性的關(guān)系，因此通過(guò)任何非迭代的算法都幾乎不可能找到精確解。目前Landweber迭代算法具有良好的成像質(zhì)量以及精度，是目前應(yīng)用最廣泛的圖像重建算法之一。

從原理上來(lái)講，Landweber算法是最速下降法的一種變形，其本質(zhì)是通過(guò)迭代來(lái)計(jì)算靈敏度矩陣的廣義逆矩陣[18]。一般地，Landweber迭代法更多的用來(lái)解決Fredholm積分方程的典型病態(tài)問(wèn)題。

本文將采用LBP重建的圖像作為迭代的初值，Landweber迭代算法如式(5)所示：

(5)

式中,β是增益因子，用來(lái)控制收斂速度，gk+1為第k步迭代的圖像灰度值，g0為利用線性反投影算法計(jì)算得到的灰度初始值，S為靈敏度矩陣，U為測(cè)量電壓值。

2 基于人群搜索算法優(yōu)化的EMT圖像重建

基于人群搜索算法優(yōu)化的電磁層析圖像重建的流程主要分為兩個(gè)部分：第一部分是利用人群搜索算法對(duì)利用場(chǎng)量提取法得到的靈敏度矩陣進(jìn)行預(yù)處理，降低其條件數(shù)；第二部分則是根據(jù)優(yōu)化之后靈敏度矩陣?yán)脠D像重建算法進(jìn)行成像。針對(duì)第二部分，本文所采用圖像重建算法為L(zhǎng)andweber迭代算法，為了便于表示，采用此方式進(jìn)行圖像重建結(jié)果的算法命名為O_Landweber迭代算法，流程如下：

1)參數(shù)的初始化，在人群搜索算法中，設(shè)定種群規(guī)模為sizepop，迭代次數(shù)為maxgen，空間維數(shù)為D，最大隸屬度值Umax，最小隸屬度值Umin，權(quán)重最大值Wmax，權(quán)重最小值Wmin，種群個(gè)體最大值為popmax，種群個(gè)體最小值為popmin，EMT逆問(wèn)題剖分網(wǎng)格數(shù)為N，Landweber迭代算法的迭代次數(shù)為V。

2)按照1.2節(jié)的方法計(jì)算靈敏度矩陣S，維度為D*N。

3)按照參數(shù)初始化中定義的種群中個(gè)體的最大值和最小值，完全隨機(jī)化生成初始種群T，維度為1*D。

4)根據(jù)人群搜索算法以適應(yīng)度值最小為最優(yōu)進(jìn)行搜索的特點(diǎn)[19]，對(duì)文獻(xiàn)[19]中提出的適應(yīng)度值的求取公式進(jìn)行取負(fù)操作，并且以矩陣T中的元素為對(duì)角線元素生成對(duì)角陣T1，然后計(jì)算種群中的個(gè)體的適應(yīng)度值。定義如式(6)和式(7)：

T1=diag(T)

(6)

(7)

其中:P是懲罰函數(shù)，定義如式(8)所示：

(8)

式中,cond(T1)表示T1的條件數(shù)。

5)按照1.3節(jié)的方法對(duì)搜索步長(zhǎng)、搜索方向進(jìn)行更新。

6)如果達(dá)到迭代次數(shù)，進(jìn)行下一步，否則返回第(4)步。

7)我們依據(jù)所改善的靈敏度矩陣的條件數(shù)占初始靈敏度條件數(shù)的比例來(lái)最終衡量靈敏度矩陣的優(yōu)化效果，定義如式(9)所示，delta越大表明優(yōu)化效果越好。在人群搜索算法中，本文采取雙停止策略，即若連續(xù)循環(huán)3次delta的變化量都不超過(guò)0.1%，或者已經(jīng)循環(huán)優(yōu)化了20次，則進(jìn)入(8)，否則進(jìn)入(4)繼續(xù)更新T。

(9)

式中,T1best表示該次循環(huán)后種群中的最佳個(gè)體。

8)在得到最終優(yōu)化的T1best之后，根據(jù)式(10)更新靈敏度矩陣以及測(cè)量電壓。本文將按照線性反投影算法得到的初值進(jìn)行Landweber迭代算法的計(jì)算，將更新后的靈敏度矩陣以及測(cè)量電壓應(yīng)用到式(5)所定義的Landweber迭代算法中進(jìn)行迭代，直到達(dá)到迭代次數(shù)V時(shí) 停止，進(jìn)行圖像重建。

(10)

式中,U1和S1分別表示優(yōu)化之后的測(cè)量電壓以及靈敏度矩陣。

基于SOA優(yōu)化的EMT圖像重建的流程圖如圖3所示。

圖3 基于SOA優(yōu)化的EMT圖像重建流程圖

3 仿真實(shí)驗(yàn)與分析

為了驗(yàn)證本文所提出的利用人群搜索算法改善靈敏度矩陣的病態(tài)程度，然后進(jìn)行圖像重建的有效性，主要選取了6種物場(chǎng)分布，采用傳統(tǒng)的8線圈EMT系統(tǒng)結(jié)構(gòu)進(jìn)行仿真[20-21]。為了比較不同算法的仿真效果，本文采用LBP，Tikhonov正則化，Landweber迭代算法以及第2節(jié)提出的O_Landweber迭代算法進(jìn)行成像。

3.1 仿真實(shí)驗(yàn)環(huán)境設(shè)定

與二維仿真相比，三維仿真更能夠模擬實(shí)際空間的復(fù)雜情況，準(zhǔn)確地表達(dá)出各物體之間相互關(guān)系，得到的仿真結(jié)果也就更為準(zhǔn)確，故本文采用電磁場(chǎng)有限元仿真軟件對(duì)仿真模型進(jìn)行三維建模。利用Matlab來(lái)控制在仿真過(guò)程中激勵(lì)線圈的切換，激勵(lì)電流的設(shè)置以及敏感場(chǎng)中物理量的獲取。

在建模時(shí)設(shè)定敏感場(chǎng)半徑為50 mm，線圈半徑為16 mm，激勵(lì)頻率為200 kHz，激勵(lì)電流為10 A，線圈材質(zhì)為銅。在進(jìn)行實(shí)驗(yàn)時(shí)，所放入的圓柱形物質(zhì)為半徑為8 mm的銅棒，其余為空氣，本文所建立的模型如圖4所示。

圖4 三維仿真模型

設(shè)定電磁層析成像的剖分單元格數(shù)N為322，如圖5所示；在人群搜索算法中，設(shè)定種群規(guī)模為sizepop為100，迭代次數(shù)為maxgen為50，空間維數(shù)D為28，最大隸屬度值Umax為0.950 0，最小隸屬度值Umin為0.011 1，權(quán)重最大值Wmax為0.9，權(quán)重最小值Wmin為0.1，種群個(gè)體最大值為popmax為1，種群個(gè)體最小值popmin為0；設(shè)定Landweber迭代算法的迭代次數(shù)V為1 000次[22]。

圖5 逆問(wèn)題網(wǎng)格剖分結(jié)果

3.2 仿真實(shí)驗(yàn)結(jié)果與分析

本文的仿真是在Matlab以及電磁場(chǎng)有限元仿真軟件中進(jìn)行的，利用Matlab來(lái)控制有限元仿真軟件的仿真進(jìn)程、仿真數(shù)據(jù)獲取以及圖像重建。仿真過(guò)程主要分為3個(gè)部分：第一部分為靈敏度矩陣的計(jì)算，第二部分為邊界測(cè)量值獲取，第三部分為圖像重建。首先獲取靈敏度矩陣，計(jì)算每個(gè)線圈在單位電流激勵(lì)下在被測(cè)物場(chǎng)中產(chǎn)生的磁場(chǎng)，然后確定激勵(lì)和檢測(cè)線圈的位置，從場(chǎng)量中提取靈敏度[15]，從而得到整個(gè)物場(chǎng)的靈敏度矩陣；獲取邊界測(cè)量值則采用單線圈激勵(lì)的策略，在電磁場(chǎng)有限元仿真軟件中，對(duì)3.1節(jié)中建立的模型中的8個(gè)銅線圈依次激勵(lì)，共可得到8×8個(gè)測(cè)量值；最后依據(jù)靈敏度矩陣和邊界測(cè)量值采用不同的圖像重建算法在MATLAB中進(jìn)行圖像重建。

本文通過(guò)采用循環(huán)利用人群搜索算法來(lái)改善靈敏度矩陣的病態(tài)程度，適應(yīng)度值和靈敏度矩陣條件數(shù)的收斂曲線如圖5所示，本次實(shí)驗(yàn)共循環(huán)了11次，每次循環(huán)次數(shù)maxgen為50，為了將所有的循環(huán)過(guò)程都呈現(xiàn)在一副圖片上，設(shè)定橫坐標(biāo)為總迭代次數(shù)，共迭代11*50=550代，縱坐標(biāo)為每次迭代所對(duì)應(yīng)的適應(yīng)度值以及靈敏度矩陣條件數(shù)。

圖6 適應(yīng)度值和靈敏度矩陣條件數(shù)收斂曲線

由圖6可知，根據(jù)適應(yīng)度值收斂曲線可以看出，迭代次數(shù)小于300時(shí)，隨著迭代次數(shù)的增加適應(yīng)度值下降明顯，當(dāng)?shù)螖?shù)大于300時(shí)，適應(yīng)度值下降趨于平緩，逐漸收斂，其表現(xiàn)為靈敏度矩陣條件數(shù)也趨于收斂。根據(jù)靈敏度矩陣條件數(shù)收斂曲線可以看出，用人群搜索算法對(duì)靈敏度矩陣進(jìn)行優(yōu)化，能夠有效降低其條件數(shù)。在迭代次數(shù)小于100時(shí)靈敏度矩陣條件數(shù)下降明顯，當(dāng)?shù)螖?shù)大于300時(shí)，適應(yīng)度值趨于穩(wěn)定，與適應(yīng)度值收斂曲線表現(xiàn)基本一致。

靈敏度矩陣條件數(shù)的改善量收斂曲線如圖7所示。

圖7 靈敏度矩陣條件數(shù)改善量收斂曲線

由圖7可知，在人群搜索算法循環(huán)迭代的過(guò)程中，靈敏度矩陣條件數(shù)改善量隨循環(huán)次數(shù)的增加不斷增加；在增加的過(guò)程中，曲線的斜率逐漸減小，表明靈敏度矩陣條件數(shù)的改善速率不斷減小，數(shù)值差異也不斷縮小，即可以停止迭代優(yōu)化。結(jié)合圖6來(lái)看靈敏度矩陣的最初的條件數(shù)為632.621 8，優(yōu)化之后降至77.009 9，共降低了87.83%，有效地改善了靈敏度矩陣的病態(tài)程度仿真得到的圖像重建結(jié)果如表1所示。

表1 仿真圖像重建結(jié)果

如表1所示，第一列為仿真的物體的真實(shí)分布圖像，第二列為采用LBP算法重建所得到的圖像，第三列為采用Tikhonov正則化算法重建所得到的圖像，第四列為采用Landweber迭代算法重建所得到的圖像，第五列為本文提出O_Landweber迭代算法重建所得到的圖像。相對(duì)于LBP算法以及Tikhonov正則化算法，通過(guò)Landweber迭代算法以及本文所提出O_Landweber算法重建所得到的圖像能夠較為準(zhǔn)確的重建出物體所在的位置，并且能夠減少偽影。由表中第四列Landweber算法重建的圖像可知，在物體周圍存在較多的偽影，物體的位置還不夠精確，本文提出的O_Landweber迭代算法能夠進(jìn)一步減少物體周圍的偽影，還原物體的實(shí)際位置。

為了以定量的方式來(lái)對(duì)以上4種算法的成像質(zhì)量以及成像精度進(jìn)行比較，本文采用圖像相對(duì)誤差(IE,image error)以及相關(guān)系數(shù)(CC,correlation coefficient)作為評(píng)價(jià)指標(biāo)進(jìn)行比較[23]。圖像相關(guān)系數(shù)是衡量變量之間線性相關(guān)程度的指標(biāo)；圖像誤差是重建圖像與原始圖像之間的差異程度。若將采用兩種不同成像算法同一種物場(chǎng)分布的仿真結(jié)果進(jìn)行比較，如果其中一種仿真圖像的圖像相對(duì)誤差更小，而相關(guān)系數(shù)更大，那么就說(shuō)明該仿真結(jié)果的圖像重建的質(zhì)量最好。定義如式(11)和(12)所示：

(11)

(12)

表2和表3給出了6中不同分布下4種不同的算法所得到的圖像重建結(jié)果的圖像誤差以及圖像相關(guān)系數(shù)。

表2 不同圖像重建算法對(duì)各類型分布的圖像相關(guān)系數(shù)

表3 不同圖像重建算法對(duì)各典型分布的圖像誤差模型

由表2、表3可知，當(dāng)被測(cè)物場(chǎng)中含有一個(gè)物體時(shí)，本文提出的O_Landweber迭代算法相比Landweber迭代算法的圖像相關(guān)系數(shù)能夠提升22%以上，圖像誤差能夠降低56%以上；當(dāng)被測(cè)物場(chǎng)中含有兩個(gè)物體時(shí)，不論是圖像相關(guān)系數(shù)還是圖像誤差，改善幅度均沒(méi)有一個(gè)物體時(shí)明顯，但仍有提升，圖像相關(guān)系數(shù)能夠提升8%以上，圖像誤差能夠降低28%以上。本文提出的O_Landweber迭代算法相對(duì)于LBP、Tikhonov正則化以及Landweber迭代算法的圖像相關(guān)系數(shù)更大，圖像誤差更小，這表明其所成圖像與原始分布的圖像誤差較小且線性相關(guān)程度較高，與仿真圖像重建結(jié)果表現(xiàn)的基本一致。

4 結(jié)束語(yǔ)

由于EMT自身特性，導(dǎo)致靈敏度矩陣往往都具有病態(tài)性。為了進(jìn)一步提高電磁層析成像的成像質(zhì)量，本文提出先通過(guò)人群算法降低靈敏度矩陣病態(tài)程度，然后再進(jìn)行圖像重建來(lái)提高成像質(zhì)量。通過(guò)循環(huán)利用人群搜索算法對(duì)靈敏度矩陣進(jìn)行不斷優(yōu)化的方法，使靈敏度矩陣的條件數(shù)能夠在迭代300代之后趨于收斂，相對(duì)于初始靈敏度矩陣條件數(shù)能夠降低87.83%，靈敏度矩陣病態(tài)程度得到明顯改善。仿真實(shí)驗(yàn)表明：與傳統(tǒng)的直接用于成像的電磁層析成像算法相比，本文所提出的O_Landweber迭代算法能夠?qū)D像相關(guān)系數(shù)最高提升45%，降低圖像誤差最高75%；與LBP算法、Tikhonov正則化算法，Landweber迭代算法相比，本文所提出的算法能夠明顯消除偽影，提高成像質(zhì)量，進(jìn)一步精確物體的位置。