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

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

中圖分類號：TH701 文獻標志碼：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ù)的重要組成部分，可以測量物場空間導電率/磁導率的分布情況^［1-2］。根據(jù)電磁感應原理，激勵源通入交變電流，物場空間產(chǎn)生交變磁場，即主磁場。由于物場空間介質(zhì)不同，由導電物體中感應的渦流與初級磁場疊加產(chǎn)生二次磁場。檢測線圈均勻分布在物場空間周圍，獲得不同投影方向的測量電壓。物場空間的電導率分布可以通過合適的圖像重建算法重構(gòu)。

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

1 EMT實驗室系統(tǒng)

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

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

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

其中A₁和A₂分別為檢測線圈2個端點處的矢量磁位值，表達式為

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

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

U=Sg（3）

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

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

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

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

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

搜索方向p_k定義為

其中r_k=U-g_k S為殘差向量，設置允許誤差e，當‖r_k ‖gt;e時，計算精確的搜索步長λ_k：

共軛梯度算法在迭代過程中對迭代的步長和約束條件要求較高，為了保證全局收斂性，一般會采用非精確的線性搜索方式確定搜索步長，常見的非精確線性搜索標準有Armijo線性搜索、Wolfe線性搜索和Golklstein線性搜索^［10］。β_k是共軛迭代參數(shù)，它的不同取法對應不同的共軛梯度算法，β_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算法的分子相同，而且各種算法具有不同的收斂特點和數(shù)值優(yōu)劣性質(zhì)。

3.1 譜共軛梯度算法

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

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

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

第2步：由Armijo搜索準則確定搜索步長λ_k使其滿足

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

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

第3步：計算

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

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

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

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

3.2 混合共軛梯度算法

PRP算法，HS算法和LS算法具有相同的r^T_k（r_k-r_k-1）項，數(shù)值表現(xiàn)較好，但是無法保證算法的充分下降性。而FR，CD和DY共軛梯度算法與合適的線搜索結(jié)合，收斂性較好，但是在實際計算中，其數(shù)值表現(xiàn)通常不如前3種方法。另外，F(xiàn)R算法存在的致命缺陷是一旦在某步生成小步長，后面將連續(xù)產(chǎn)生小步長。與此相反，HS算法和PRP算法則具有自動啟動趨勢，可以避免連續(xù)產(chǎn)生小步長^［11-14］。為了尋求具有全局收斂性且數(shù)值表現(xiàn)良好的算法，考慮將不同類型共軛梯度算法混合。在EMT技術(shù)中，結(jié)合FR算法和PRP算法的優(yōu)點，組合成一個收斂效果更好，數(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ù)，對式（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算法的基礎上，可以根據(jù)殘差自適應的調(diào)整大小。根據(jù)公式（7）可知，當出現(xiàn)小步長時，即r_k≈r_k-1，r^T_kr^T_k-1→1導致β^MCG_k→0，下一個搜索方向改成負梯度，保證了算法的良好收斂性和數(shù)值特征。

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

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

4 實驗室EMT系統(tǒng)圖像重建

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

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

相關系數(shù)（correlation-coefficient，CC）和圖像誤差（image error，IE）是評價EMT重建圖像質(zhì)量的重要指標^［15］，其中IE表達式為

CC表達式為

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

從表2可以看出，譜參數(shù)混合共軛梯度算法的圖像誤差值比其他2種對比算法小，說明其重建圖像的質(zhì)量較高，并且相關系數(shù)的絕對值較其他算法大，具有較好的相關度，逆問題的重建函數(shù)與真實的目標函數(shù)之間相關度較強，能比較準確地擬合真實目標函數(shù)，更接近真實的解?；旌瞎曹椞荻人惴ňC合了FR算法和PRP算法的優(yōu)點，成像質(zhì)量高于這2種算法，在運行速度和收斂性上都具有更大的優(yōu)勢。

5 結(jié) 語

在EMT技術(shù)中，為了提高共軛梯度算法重建圖像的質(zhì)量，從2個方面修正傳統(tǒng)的非線性共軛梯度算法：1）對搜索方向進行修正，通過給定搜索方向一個譜參數(shù)，控制殘差向量和上一步搜索方向的權(quán)重，來控制下一步的方向，通過對比收斂情況發(fā)現(xiàn)譜共軛梯度算法數(shù)值表現(xiàn)更好，大大改善了圖像的質(zhì)量；2）從共軛參數(shù)入手，結(jié)合FR算法快速收斂和PRP算法良好的數(shù)值表現(xiàn)的優(yōu)點，將2種算法按照一定比例混合得到混合共軛梯度算法，列出收斂圖和重建圖像，利用圖像誤差和相關度2個指標評價修正算法，發(fā)現(xiàn)譜參數(shù)混合共軛梯度算法的收斂性和成像質(zhì)量高于傳統(tǒng)算法，圖像重建質(zhì)量較高，能有效實現(xiàn)EMT實驗室系統(tǒng)被測物質(zhì)重現(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ù)理論的平面EMT圖像融合算法［J］.儀器儀表學報，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］.沈陽：東北大學，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］張靜.廣義隨機線性互補問題的數(shù)值求解方法［D］.天津：天津大學，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］.中國電機工程學報，2021，41（15）：5351-5361.

【責任編輯：封文江】