電磁層析成像技術(shù)中的譜參數(shù)混合共軛梯度法

摘 要：電磁層析成像技術(shù)是基于電磁感應(yīng)原理的過程檢測(cè)技術(shù)，圖像重建是解決其反問題的關(guān)鍵。基于共軛梯度算法，修正搜索方向和迭代參數(shù)，推導(dǎo)出基于譜參數(shù)的混合共軛梯度算法，以提高圖像重建質(zhì)量和收斂性。介紹了電磁層析成像系統(tǒng)的實(shí)驗(yàn)室模型和靈敏度分布矩陣，作為反問題基礎(chǔ)；對(duì)搜索方向進(jìn)行修正，提出一個(gè)新的譜參數(shù)，用來控制新舊搜索方向之間的角度，根據(jù)每一步迭代的殘差結(jié)果，尋找最優(yōu)解；改善共軛參數(shù)，分析非線性共軛梯度算法中FR算法（fletcher-reeves method）和PRP算法（polak-ribiere-polyak method）的各自優(yōu)勢(shì)，將它們按照一定比例進(jìn)行混合，得到一種新的混合共軛梯度算法；將譜參數(shù)混合共軛梯算法應(yīng)用到電磁層析成像實(shí)驗(yàn)室系統(tǒng)中，構(gòu)建了3種典型實(shí)驗(yàn)室模型，對(duì)比譜參數(shù)混合共軛梯度算法與傳統(tǒng)算法并作出評(píng)價(jià)。實(shí)驗(yàn)結(jié)果表明：譜參數(shù)混合共軛梯度算法的重建圖像的質(zhì)量更高，具有較好的數(shù)值表現(xiàn)；混合共軛梯度算法結(jié)合了FR算法和PRP算法的優(yōu)點(diǎn)，收斂速度比PRP算法快，成像質(zhì)量比其他算法高。

關(guān) 鍵 詞：電磁層析成像技術(shù); 共軛梯度; 圖像重建關(guān) 鍵 詞：分子篩; 鈦硅沸石; 無溶劑法; 硅鈦摩爾比

中圖分類號(hào)：TH701 文獻(xiàn)標(biāo)志碼：A

doi：10.3969/j.issn.1673-5862.2024.03.005

Spectral parameters hybrid conjugate gradient method in electromagnetic tomography technology

CUI Song^1，2， LYU Yan^1，2， CHEN Lanfeng^1，2LI Liu¹， LUO Yue¹， WANG Chengyu¹， AN Zhigang¹， WANG Yanli²

（1. College of Physical Science and Technology， Shenyang Normal University， Shenyang 110034， China）（1. College of Physical Science and Technology， Shenyang Normal University， Shenyang 110034， China; 2. College of Grain Science and Technology， Shenyang Normal University， Shenyang 110034， China）

Abstract：Electromagnetic tomography is a kind of process detection technology based on the principle of electromagnetic induction. Image reconstruction is the key to solve its inverse problem. Based on the conjugate gradient algorithm， the search direction and iteration parameters are modified. A hybrid conjugate gradient algorithm based on spectral parameters is derived to improve the quality and convergence of image reconstruction. The laboratory model and sensitivity distribution matrix of the electromagnetic tomography technology system are introduced as the basis of the inverse problem. A new spectral parameter conjugate gradient algorithm is proposed to modify the search direction， which is used to control the angle between the old and new search directions. The search direction is determined according to the residual of each step in order to find the optimal solution. Combining the advantages of FR algorithm and PRP algorithm in nonlinear conjugate gradient algorithm， they are mixed by certain proportion to obtain a new hybrid conjugate gradient algorithm. The spectral parameter hybrid conjugate gradient method is applied to the electromagnetic tomography technology laboratory system， three typical laboratory models are constructed， and the spectral parameter hybrid conjugate gradient algorithm and the traditional algorithm are compared and evaluated. The experimental results show that the spectral parameter hybrid conjugate gradient algorithm has higher quality of reconstructed image and better numerical performance， and it combines the advantages of FR algorithm and PRP algorithm， with faster convergence speed and higher imaging quality than other algorithms.

Key words：electromagnetic tomography technology; conjugate gradient algorithm; image reconstruction

電磁層析成像技術(shù)（electromagnetic tomography technology，EMT）是電層析成像技術(shù)的重要組成部分，可以測(cè)量物場(chǎng)空間導(dǎo)電率/磁導(dǎo)率的分布情況^［1-2］。根據(jù)電磁感應(yīng)原理，激勵(lì)源通入交變電流，物場(chǎng)空間產(chǎn)生交變磁場(chǎng)，即主磁場(chǎng)。由于物場(chǎng)空間介質(zhì)不同，由導(dǎo)電物體中感應(yīng)的渦流與初級(jí)磁場(chǎng)疊加產(chǎn)生二次磁場(chǎng)。檢測(cè)線圈均勻分布在物場(chǎng)空間周圍，獲得不同投影方向的測(cè)量電壓。物場(chǎng)空間的電導(dǎo)率分布可以通過合適的圖像重建算法重構(gòu)。

EMT技術(shù)包括正問題和反問題研究。正問題是在已知物場(chǎng)空間分布時(shí)，獲得并分析檢測(cè)值的過程，比如，建立物理模型、提取檢測(cè)值、構(gòu)建靈敏度矩陣等。而圖像重建是反問題的主要內(nèi)容，目前主要的圖像重建算法有線性反投影算法（linear back projection，LBP）算法、Tikhonov正則化算法、迭代算法、共軛梯度算法及神經(jīng)網(wǎng)絡(luò)方法等^［3-5］。

1 EMT實(shí)驗(yàn)室系統(tǒng)

8線圈EMT實(shí)驗(yàn)室系傳感器陣列如圖1所示，8個(gè)線圈均勻、同高度排列在圓柱物場(chǎng)壁周圍，交替作為激勵(lì)和檢測(cè)，采用同一套激勵(lì)源輪流通過SMA（sub-miniature-A，反極性公頭）接口與激勵(lì)線圈連接形成激勵(lì)磁場(chǎng)。物場(chǎng)空間為圓柱形，管壁（屏蔽層）采用厚5mm的鋁材質(zhì)，激勵(lì)頻率為10MHz，實(shí)驗(yàn)裝置實(shí)物圖如圖1（a）所示。圖1（b）是線圈傳感器陣列的截面圖，其中R₁=200mm為物場(chǎng)空間半徑，R₂=225mm為屏蔽層內(nèi)徑，R₃=300mm為屏蔽層外徑。單個(gè)線圈傳感器的結(jié)構(gòu)如圖1（c）所示，激勵(lì)/檢測(cè)線圈的組成為12匝Φ0.5mm實(shí)心銅線，Φ50mm PVC絕緣骨架^［6］。

2 EMT技術(shù)反問題研究

根據(jù)檢測(cè)線圈獲得的檢測(cè)值來重構(gòu)物場(chǎng)空間電導(dǎo)率/磁導(dǎo)率分布的過程稱為電磁層析成像技術(shù)的反問題，它是正問題的逆過程。檢測(cè)值與物場(chǎng)空間分布之間非線性關(guān)系可以表達(dá)為

其中A₁和A₂分別為檢測(cè)線圈2個(gè)端點(diǎn)處的矢量磁位值，表達(dá)式為

式中：φ為通過檢測(cè)線圈的磁通量；μ，σ分別表示磁導(dǎo)率、電導(dǎo)率；n表示檢測(cè)線圈匝數(shù)；l表示沿管道軸向的長度；ω為角頻率。

由式（1）和（2）獲得的檢測(cè)值數(shù)量遠(yuǎn)遠(yuǎn)小于被測(cè)物場(chǎng)空間所描述的像素點(diǎn)，因此，該類方程是欠定的。為了解決這類方程的軟場(chǎng)問題，一般使用數(shù)值方法將非線性問題線性化，將檢測(cè)值近似為

U=Sg（3）

其中U∈R^M為檢測(cè)電壓矩陣;S∈R^M×N為空間靈敏度矩陣；g∈R^N為灰度矩陣；M是獨(dú)立測(cè)量值個(gè)數(shù)（M=n×n，n為線圈個(gè)數(shù)）；N是采用有限元方法時(shí)的剖分單元數(shù)。

靈敏度矩陣S_ij （k）是利用有限元方法（finite element method，F(xiàn)EM）剖分出的網(wǎng)格中第k個(gè)單元充滿被測(cè)介質(zhì)時(shí)，第i個(gè)線圈作激勵(lì)，第j個(gè)線圈的檢測(cè)值。獲取靈敏度矩陣的方法有實(shí)驗(yàn)擾動(dòng)法、模擬擾動(dòng)法和場(chǎng)量公式提取法^［7-8］。模擬擾動(dòng)法將物場(chǎng)空間剖分成若干有限元單元，在每個(gè)單元加入擾動(dòng)介質(zhì)，每個(gè)線圈輪流作為激勵(lì)源，遍歷所有單元獲得每個(gè)有限元的靈敏度，構(gòu)成靈敏度矩陣。

3 基于譜參數(shù)的混合共軛梯度法

正則化方法是解決過程層析技術(shù)反演問題的主要手段，共軛梯度算法是一種特殊的正則化方法 ^［9］。根據(jù)變分原理和最速下降思想，在每個(gè)迭代點(diǎn)g_k處，搜索方向關(guān)于對(duì)稱正定矩陣A共軛，滿足

p^T_iAp_j=0，（i≠j）（4）

搜索方向p_k定義為

其中r_k=U-g_k S為殘差向量，設(shè)置允許誤差e，當(dāng)‖r_k ‖gt;e時(shí)，計(jì)算精確的搜索步長λ_k：

共軛梯度算法在迭代過程中對(duì)迭代的步長和約束條件要求較高，為了保證全局收斂性，一般會(huì)采用非精確的線性搜索方式確定搜索步長，常見的非精確線性搜索標(biāo)準(zhǔn)有Armijo線性搜索、Wolfe線性搜索和Golklstein線性搜索^［10］。β_k是共軛迭代參數(shù)，它的不同取法對(duì)應(yīng)不同的共軛梯度算法，β_k可表示為

根據(jù)β_k參數(shù)取值不同定義共軛梯度算法為FR算法、PRP算法、HS算法（hestenes-stiefel method）、CD算法（conjugate descent method）、DY算法（dai-yuan method）和LS算法（liu-storey method）^［11］。從公式上看，PRP算法，HS算法，LS算法的分母相同，而FR算法，CD算法和DY算法的分子相同，而且各種算法具有不同的收斂特點(diǎn)和數(shù)值優(yōu)劣性質(zhì)。

3.1 譜共軛梯度算法

在解決EMT技術(shù)的逆問題時(shí)，基于Armijo搜索準(zhǔn)則，設(shè)計(jì)更適合EMT系統(tǒng)θ_k參數(shù)的FR譜共軛梯度算法，表達(dá)式為

θ_k譜參數(shù)由搜索方向p_（k-1）和殘差向量‖r_k‖組成，具體的算法步驟如下：

第1步：給定初始值g₀=0，δ=0.6，σ=0.4，ε=0.0001，k=1，計(jì)算r₁=Ag₀-b，p₁=-r₁，若‖r_k‖≤ε，停止迭代，否則進(jìn)入下一步。

第2步：由Armijo搜索準(zhǔn)則確定搜索步長λ_k使其滿足

ψ（g_k+λ_kp_k）≤ψ（g_k）+δλ_kr^T_kp_k

ψ（g_k+λ_kp_k）≥ψ（g_k）+σλ_kr^T_kp_k（9）

第3步：計(jì)算

g_k+1=g_k+λ_kp_k，k≥1（10）

第4步：計(jì)算r_k+1，若‖r_k+1‖≤ε，則停止迭代；否則循環(huán)k=k+1。

第5步：計(jì)算β_k，θ_k和搜索方向p_k=-θ_kr_k+β_kp_k-1。

基于FR算法設(shè)計(jì)的譜共軛梯度算法引入了2個(gè)方向參數(shù)變量，通過相鄰2次的殘差向量改變前后2次搜索方向的比值，從而得到新的搜索方向。此種算法不僅具有良好的數(shù)值表現(xiàn)，還有較快的收斂速度。

3.2 混合共軛梯度算法

PRP算法，HS算法和LS算法具有相同的r^T_k（r_k-r_k-1）項(xiàng)，數(shù)值表現(xiàn)較好，但是無法保證算法的充分下降性。而FR，CD和DY共軛梯度算法與合適的線搜索結(jié)合，收斂性較好，但是在實(shí)際計(jì)算中，其數(shù)值表現(xiàn)通常不如前3種方法。另外，F(xiàn)R算法存在的致命缺陷是一旦在某步生成小步長，后面將連續(xù)產(chǎn)生小步長。與此相反，HS算法和PRP算法則具有自動(dòng)啟動(dòng)趨勢(shì)，可以避免連續(xù)產(chǎn)生小步長^［11-14］。為了尋求具有全局收斂性且數(shù)值表現(xiàn)良好的算法，考慮將不同類型共軛梯度算法混合。在EMT技術(shù)中，結(jié)合FR算法和PRP算法的優(yōu)點(diǎn)，組合成一個(gè)收斂效果更好，數(shù)值表現(xiàn)優(yōu)秀的混合共軛梯度算法。將β^PRP_k和β^FR_k按照一定的權(quán)重相加，得到混合共軛梯度算法的β^MCG_k參數(shù)，其具體形式如下：

β^MCG_k=μ_kβ^PRP_k+（1-μ_k）β^FR_k（11）

定義μ_k=‖r_k‖‖r_k-1‖，將前后2次殘差向量的夾角作為參數(shù)，對(duì)式（11）化簡為

β^MCG_k=μ_kβ^PRP_k+（1-μ_k）β^FR_k=μ_k‖r_k‖²-μ_kr^T_kr_k-1+（1-μ_k）‖r_k‖²‖r_k-1‖²=β^FR_k1-r^T_kr^T_k-1 （12）

混合共軛梯度算法的β^MCG_k是在FR算法的基礎(chǔ)上，可以根據(jù)殘差自適應(yīng)的調(diào)整大小。根據(jù)公式（7）可知，當(dāng)出現(xiàn)小步長時(shí)，即r_k≈r_k-1，r^T_kr^T_k-1→1導(dǎo)致β^MCG_k→0，下一個(gè)搜索方向改成負(fù)梯度，保證了算法的良好收斂性和數(shù)值特征。

設(shè)計(jì)實(shí)驗(yàn)室EMT系統(tǒng)的3種典型模型，如圖2所示，保證3種模型硬件參數(shù)相同，采用不同共軛梯度算法進(jìn)行圖像重建，得出誤差隨迭代次數(shù)變化的趨勢(shì)，分析算法的收斂性，迭代誤差收斂圖如圖3所示。其中藍(lán)色曲線表示PRP算法的收斂情況，綠色曲線表示FR算法的收斂情況，紅色曲線表示混合共軛梯度算法的收斂情況。

由圖3可知： PRP算法的變化較慢，收斂速度較慢，但是趨于真實(shí)值之后會(huì)比較穩(wěn)定；FR算法收斂較快，且存在小步長的情況，整體的數(shù)值表現(xiàn)略弱；而混合共軛梯度算法的收斂速度較FR算法略慢些，但是收斂較平滑，誤差值變化范圍較小，沒有出現(xiàn)振蕩的情況，其收斂情況較FR算法和PRP算法優(yōu)化。

4 實(shí)驗(yàn)室EMT系統(tǒng)圖像重建

為了測(cè)試不同共軛梯度算法在EMT實(shí)驗(yàn)系統(tǒng)中的成像質(zhì)量，將被測(cè)物體分別放置于靠近線圈端，即高靈敏度位置，如圖2（a）中模型；多物體分布成像效果模型，如圖2（b）和（c）所示。被測(cè)物體選用趨膚效應(yīng)較強(qiáng)的銅介質(zhì)，直徑為15mm的銅棒（磁導(dǎo)率μ≈1，電導(dǎo)率γ=1.7×10^-8Ω·m）。基于實(shí)驗(yàn)室EMT系統(tǒng)3種典型模型，分別采用3種共軛梯度算法FR算法、PRP算法和基于譜函數(shù)的混合共軛梯度MCG算法進(jìn)行圖像重建，成像結(jié)果如圖4所示。表1列出了3種算法針對(duì)3種典型模型的迭代次數(shù)和運(yùn)行時(shí)間。

由圖4可知，譜參數(shù)混合共軛梯度算法的重建圖像清晰度更高且邊界清晰，能夠更準(zhǔn)確地反映出擾動(dòng)目標(biāo)的位置，結(jié)合表1的收斂時(shí)間和迭代次數(shù)數(shù)據(jù)可以證明，混合算法不僅具有相對(duì)較快的收斂速度，且成像精度大大提高，更適合用于實(shí)時(shí)檢測(cè)系統(tǒng)中，如醫(yī)學(xué)成像或者高鐵軌道探測(cè)等對(duì)成像精度要求較高的領(lǐng)域。

相關(guān)系數(shù)（correlation-coefficient，CC）和圖像誤差（image error，IE）是評(píng)價(jià)EMT重建圖像質(zhì)量的重要指標(biāo)^［15］，其中IE表達(dá)式為

CC表達(dá)式為

其中：_i為重建圖像中每個(gè)剖分單元的灰度值；g為原始物場(chǎng)的二維圖像的灰度值；g＾_和_i分別為_i和g的平均值。IE越小，說明重建圖像的灰度值越接近原始圖像的值，即圖像的精度越高。CC越大，說明圖像與原始圖像的相關(guān)度越高，圖像質(zhì)量越高。3種典型模型重建圖像評(píng)價(jià)結(jié)果見表2。

從表2可以看出，譜參數(shù)混合共軛梯度算法的圖像誤差值比其他2種對(duì)比算法小，說明其重建圖像的質(zhì)量較高，并且相關(guān)系數(shù)的絕對(duì)值較其他算法大，具有較好的相關(guān)度，逆問題的重建函數(shù)與真實(shí)的目標(biāo)函數(shù)之間相關(guān)度較強(qiáng)，能比較準(zhǔn)確地?cái)M合真實(shí)目標(biāo)函數(shù)，更接近真實(shí)的解。混合共軛梯度算法綜合了FR算法和PRP算法的優(yōu)點(diǎn)，成像質(zhì)量高于這2種算法，在運(yùn)行速度和收斂性上都具有更大的優(yōu)勢(shì)。

5 結(jié) 語

在EMT技術(shù)中，為了提高共軛梯度算法重建圖像的質(zhì)量，從2個(gè)方面修正傳統(tǒng)的非線性共軛梯度算法：1）對(duì)搜索方向進(jìn)行修正，通過給定搜索方向一個(gè)譜參數(shù)，控制殘差向量和上一步搜索方向的權(quán)重，來控制下一步的方向，通過對(duì)比收斂情況發(fā)現(xiàn)譜共軛梯度算法數(shù)值表現(xiàn)更好，大大改善了圖像的質(zhì)量；2）從共軛參數(shù)入手，結(jié)合FR算法快速收斂和PRP算法良好的數(shù)值表現(xiàn)的優(yōu)點(diǎn)，將2種算法按照一定比例混合得到混合共軛梯度算法，列出收斂圖和重建圖像，利用圖像誤差和相關(guān)度2個(gè)指標(biāo)評(píng)價(jià)修正算法，發(fā)現(xiàn)譜參數(shù)混合共軛梯度算法的收斂性和成像質(zhì)量高于傳統(tǒng)算法，圖像重建質(zhì)量較高，能有效實(shí)現(xiàn)EMT實(shí)驗(yàn)室系統(tǒng)被測(cè)物質(zhì)重現(xiàn)。

參考文獻(xiàn)：

［1］KRAJCINOVIC D，F(xiàn)ONSEKA G U.The continuous damage theory of brittle materials［J］.J Appl Mech，1981，48（4）：809-824.

PEYTO N A J，YU Z Z，LYON G，et al.An overview of electromagnetic inductance tomography：Description of three different systems［J］.Meas Sci Technol，1996，7（3）：261-271.

［2］WATSON S，WILLIAMS R J，GRIFFITHS H，et al.A transceiver for direct phase measurement magnetic induction tomography［C］//Proceedings of the 23rd IEEE Annual EMBS International Conference，Istanbul：IEEE Press，2001：3182-3184.

［3］徐笑，黃云志，韓亮.基于陣列旋轉(zhuǎn)和改進(jìn)證據(jù)理論的平面EMT圖像融合算法［J］.儀器儀表學(xué)報(bào)，2022，43（10）：136-144.

［4］WU X J，ZHAO Q，GAO M Y，et al.Image reconstruction algorithm of electromagnetic tomography based on fractional Kalman filter［J］.Flow Meas Instru，2022，86（6）：121-128.

［5］HAN M，CHENG X L，XUE Y Y.Comparison with reconstruction algorithms in magnetic induction tomography［J］.Physiol Meas，2016，37（5）：683-697.

［6］李柳.電磁層析成像技術(shù)的研究［D］.沈陽：東北大學(xué)，2013.

［7］XIONG H L，XU L G.Electromagnetic tomography：An analytical approach to detectability limits and sensitivity maps［J］.Pro Nat Sci，2000，10（7）：549-550.

［8］LIU Z，XUE F Q，YANG G，et al.Experimental measurement method of sensitivity matrix for electromagnetic tomography［J］.Sci Instrument Chinese J，2013，34（9）：1982-1988.

［9］HAO N T.Tikhonov regularization algorithm for pseudomonotone variational inequalities［J］.Acta Math Viet，2006，31（3）：283-289.

［10］張靜.廣義隨機(jī)線性互補(bǔ)問題的數(shù)值求解方法［D］.天津：天津大學(xué)，2020.

［11］LI Z，ZHOU W，LI D H.A descent modified Polak–Ribière–Polyak conjugate gradient method and its global convergence［J］.Num Analy Ima Jour，2006，26（4）：629-640.

［12］YUAN G，LU J，WANG Z.The modified PRP conjugate gradient algorithm under a non-descent line search and its application in the Muskingum model and image restoration problems［J］.Soft Comput，2021，25（8）：5867-5879.

［13］NA H，MA C.The modified conjugate gradient methods for solving a class of generalized coupled Sylvester-transpose matrix equations［J］.Math Appl Comput，2014，67（8）：1545-1558.

［14］FANG M，WANG M，SUN M，et al.A modified hybrid conjugate gradient method for unconstrained optimization［J］.J Math-UK，2021，2021（1）：5597863.

［15］霍繼偉，劉澤，王亞東，等.平面電磁層析成像鋼軌探傷［J］.中國電機(jī)工程學(xué)報(bào)，2021，41（15）：5351-5361.

【責(zé)任編輯：封文江】