基于SOA優化的電磁層析成像圖像重建研究

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

(1.北京交通大學 電子信息工程學院，北京 100044;2.中國鐵道科學研究院集團有限公司 標準計量研究所，北京 100081)

0 引言

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

在EMT系統中，圖像重建算法是EMT技術研究中的重點和難點，圖像重建算法的好壞將直接影響最終成像的準確性。目前的圖像重建算法主要有線性反投影(LBP)[7]、Tikhonov正則化[8]、Landweber迭代[9]等算法。在傳統的EMT系統中，邊界測量值的微小變化就會導致求解出的物場內的介質分布有很大變化，導致求解過程不穩定，因此EMT技術本身具有病態性，嚴重影響了重建圖像的質量。目前常采用截斷奇異值算法和Tikhonov正則化算法來改善靈敏度矩陣的病態性[10-11]，但都存在參數選取問題。Tikhonov正則化算法主要依賴于正則化因子，但是正則化因子的選取較為復雜，在實際應用中往往依據經驗進行選取，具有較大的局限性。如果正則化參數選擇較小，原始的問題就不能得到很好地近似，逆問題的病態性的改進效果不顯著；如果正則化參數選擇過大，雖然對誤差的敏感性會減少，但最終解會與正常值相差較大，達不到效果。

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

1 理論基礎

1.1 電磁層析成像數學模型

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

圖1 正問題與逆問題關系圖

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

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

(1)

式中,Udec為測量電壓，Be為激勵主磁場，Bo為由渦流產生的次級磁場，σ為電導率，μ為磁導率，ε為介電常數。

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

U=S·G

(2)

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

1.2 EMT靈敏度矩陣分析

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

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

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

1.3 人群搜索算法應用分析

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

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

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

(3)

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

(4)

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

1.4 Landweber迭代算法

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

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

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

(5)

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

2 基于人群搜索算法優化的EMT圖像重建

基于人群搜索算法優化的電磁層析圖像重建的流程主要分為兩個部分：第一部分是利用人群搜索算法對利用場量提取法得到的靈敏度矩陣進行預處理，降低其條件數；第二部分則是根據優化之后靈敏度矩陣利用圖像重建算法進行成像。針對第二部分，本文所采用圖像重建算法為Landweber迭代算法，為了便于表示，采用此方式進行圖像重建結果的算法命名為O_Landweber迭代算法，流程如下：

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

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

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

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

T1=diag(T)

(6)

(7)

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

(8)

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

5)按照1.3節的方法對搜索步長、搜索方向進行更新。

6)如果達到迭代次數，進行下一步，否則返回第(4)步。

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

(9)

式中,T1best表示該次循環后種群中的最佳個體。

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

(10)

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

基于SOA優化的EMT圖像重建的流程圖如圖3所示。

圖3 基于SOA優化的EMT圖像重建流程圖

3 仿真實驗與分析

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

3.1 仿真實驗環境設定

與二維仿真相比，三維仿真更能夠模擬實際空間的復雜情況，準確地表達出各物體之間相互關系，得到的仿真結果也就更為準確，故本文采用電磁場有限元仿真軟件對仿真模型進行三維建模。利用Matlab來控制在仿真過程中激勵線圈的切換，激勵電流的設置以及敏感場中物理量的獲取。

在建模時設定敏感場半徑為50 mm，線圈半徑為16 mm，激勵頻率為200 kHz，激勵電流為10 A，線圈材質為銅。在進行實驗時，所放入的圓柱形物質為半徑為8 mm的銅棒，其余為空氣，本文所建立的模型如圖4所示。

圖4 三維仿真模型

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

圖5 逆問題網格剖分結果

3.2 仿真實驗結果與分析

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

本文通過采用循環利用人群搜索算法來改善靈敏度矩陣的病態程度，適應度值和靈敏度矩陣條件數的收斂曲線如圖5所示，本次實驗共循環了11次，每次循環次數maxgen為50，為了將所有的循環過程都呈現在一副圖片上，設定橫坐標為總迭代次數，共迭代11*50=550代，縱坐標為每次迭代所對應的適應度值以及靈敏度矩陣條件數。

圖6 適應度值和靈敏度矩陣條件數收斂曲線

由圖6可知，根據適應度值收斂曲線可以看出，迭代次數小于300時，隨著迭代次數的增加適應度值下降明顯，當迭代次數大于300時，適應度值下降趨于平緩，逐漸收斂，其表現為靈敏度矩陣條件數也趨于收斂。根據靈敏度矩陣條件數收斂曲線可以看出，用人群搜索算法對靈敏度矩陣進行優化，能夠有效降低其條件數。在迭代次數小于100時靈敏度矩陣條件數下降明顯，當迭代次數大于300時，適應度值趨于穩定，與適應度值收斂曲線表現基本一致。

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

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

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

表1 仿真圖像重建結果

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

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

(11)

(12)

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

表2 不同圖像重建算法對各類型分布的圖像相關系數

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

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

4 結束語

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