王亞東，劉 澤，霍繼偉，趙鵬飛，王成飛，袁 偉

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

0 引言

電磁層析成像(EMT,electromagnetic tomography)是一種近二、三十年發(fā)展起來的基于電磁感應(yīng)原理的過程層析成像技術(shù)[1]，其目的是研究具有電磁特性的物質(zhì)在空間的分布。典型的EMT系統(tǒng)主要由4個(gè)部分組成：傳感器陣列(即檢測激勵(lì)線圈)、激勵(lì)信號產(chǎn)生和分配單元、數(shù)據(jù)采集和調(diào)制單元以及圖像重建和信息獲取單元。通過對激勵(lì)線圈通入交變電流，在被測物場中產(chǎn)生激勵(lì)磁場；當(dāng)被測物場內(nèi)有具有導(dǎo)磁或?qū)щ姷奈矬w放入時(shí)，原有磁場分布將會改變，線圈之間的感性耦合也將隨之改變[2-3]。通過切換不同的激勵(lì)線圈，可獲得多個(gè)投影方向的檢測信息，再運(yùn)用合適的圖像重建算法就可以重建出放入物場中物體的電導(dǎo)率或者磁導(dǎo)率的分布圖，從而獲得放入物體的信息。由此可見，EMT具有非侵入、非接觸和快速性的特點(diǎn)，目前在生物醫(yī)學(xué)、無損檢測、鋼軌探傷、異物檢測和工業(yè)煉鋼等領(lǐng)域均已展開相應(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)中，邊界測量值的微小變化就會導(dǎo)致求解出的物場內(nèi)的介質(zhì)分布有很大變化，導(dǎo)致求解過程不穩(wěn)定，因此EMT技術(shù)本身具有病態(tài)性，嚴(yán)重影響了重建圖像的質(zhì)量。目前常采用截?cái)嗥娈愔邓惴ê蚑ikhonov正則化算法來改善靈敏度矩陣的病態(tài)性[10-11]，但都存在參數(shù)選取問題。Tikhonov正則化算法主要依賴于正則化因子，但是正則化因子的選取較為復(fù)雜，在實(shí)際應(yīng)用中往往依據(jù)經(jīng)驗(yàn)進(jìn)行選取，具有較大的局限性。如果正則化參數(shù)選擇較小，原始的問題就不能得到很好地近似，逆問題的病態(tài)性的改進(jìn)效果不顯著；如果正則化參數(shù)選擇過大，雖然對誤差的敏感性會減少，但最終解會與正常值相差較大，達(dá)不到效果。

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

1 理論基礎(chǔ)

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

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

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

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

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

(1)

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

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

U=S·G

(2)

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

1.2 EMT靈敏度矩陣分析

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

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

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

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

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

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

人群搜索算法主要包含搜索步長、搜索方向以及個(gè)體位置的更新[17]，為了簡化模型，提高運(yùn)算效率，在本文中，對目標(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為種群大??；D為搜索空間維數(shù)；Ki為種群函數(shù)值按照降序排列之后的序列編號；rand()能夠生成0～1之間的隨機(jī)數(shù)，用來模仿人的隨機(jī)行為。

1.4 Landweber迭代算法

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

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

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

(5)

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

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

基于人群搜索算法優(yōu)化的電磁層析圖像重建的流程主要分為兩個(gè)部分：第一部分是利用人群搜索算法對利用場量提取法得到的靈敏度矩陣進(jìn)行預(yù)處理，降低其條件數(shù)；第二部分則是根據(jù)優(yōu)化之后靈敏度矩陣?yán)脠D像重建算法進(jìn)行成像。針對第二部分，本文所采用圖像重建算法為Landweber迭代算法，為了便于表示，采用此方式進(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ǎ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]，對文獻(xiàn)[19]中提出的適應(yīng)度值的求取公式進(jìn)行取負(fù)操作，并且以矩陣T中的元素為對角線元素生成對角陣T1，然后計(jì)算種群中的個(gè)體的適應(yīng)度值。定義如式(6)和式(7)：

T1=diag(T)

(6)

(7)

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

(8)

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

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

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

7)我們依據(jù)所改善的靈敏度矩陣的條件數(shù)占初始靈敏度條件數(shù)的比例來最終衡量靈敏度矩陣的優(yōu)化效果，定義如式(9)所示，delta越大表明優(yōu)化效果越好。在人群搜索算法中，本文采取雙停止策略，即若連續(xù)循環(huán)3次delta的變化量都不超過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)更新靈敏度矩陣以及測量電壓。本文將按照線性反投影算法得到的初值進(jìn)行Landweber迭代算法的計(jì)算，將更新后的靈敏度矩陣以及測量電壓應(yīng)用到式(5)所定義的Landweber迭代算法中進(jìn)行迭代，直到達(dá)到迭代次數(shù)V時(shí) 停止，進(jìn)行圖像重建。

(10)

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

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

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

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

為了驗(yàn)證本文所提出的利用人群搜索算法改善靈敏度矩陣的病態(tài)程度，然后進(jìn)行圖像重建的有效性，主要選取了6種物場分布，采用傳統(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)確，故本文采用電磁場有限元仿真軟件對仿真模型進(jìn)行三維建模。利用Matlab來控制在仿真過程中激勵(lì)線圈的切換，激勵(lì)電流的設(shè)置以及敏感場中物理量的獲取。

在建模時(shí)設(shè)定敏感場半徑為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ǎng)格剖分結(jié)果

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

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

本文通過采用循環(huán)利用人群搜索算法來改善靈敏度矩陣的病態(tài)程度，適應(yīng)度值和靈敏度矩陣條件數(shù)的收斂曲線如圖5所示，本次實(shí)驗(yàn)共循環(huán)了11次，每次循環(huán)次數(shù)maxgen為50，為了將所有的循環(huán)過程都呈現(xiàn)在一副圖片上，設(shè)定橫坐標(biāo)為總迭代次數(shù)，共迭代11*50=550代，縱坐標(biāo)為每次迭代所對應(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ù)收斂曲線可以看出，用人群搜索算法對靈敏度矩陣進(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)迭代的過程中，靈敏度矩陣條件數(shù)改善量隨循環(huán)次數(shù)的增加不斷增加；在增加的過程中，曲線的斜率逐漸減小，表明靈敏度矩陣條件數(shù)的改善速率不斷減小，數(shù)值差異也不斷縮小，即可以停止迭代優(yōu)化。結(jié)合圖6來看靈敏度矩陣的最初的條件數(shù)為632.621 8，優(yōu)化之后降至77.009 9，共降低了87.83%，有效地改善了靈敏度矩陣的病態(tài)程度仿真得到的圖像重建結(jié)果如表1所示。

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

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

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

(11)

(12)

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

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

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

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

4 結(jié)束語

由于EMT自身特性，導(dǎo)致靈敏度矩陣往往都具有病態(tài)性。為了進(jìn)一步提高電磁層析成像的成像質(zhì)量，本文提出先通過人群算法降低靈敏度矩陣病態(tài)程度，然后再進(jìn)行圖像重建來提高成像質(zhì)量。通過循環(huán)利用人群搜索算法對靈敏度矩陣進(jìn)行不斷優(yōu)化的方法，使靈敏度矩陣的條件數(shù)能夠在迭代300代之后趨于收斂，相對于初始靈敏度矩陣條件數(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)一步精確物體的位置。