電容成像雙共軛梯度圖像重建改進算法

馬敏，范廣永，孫穎

（中國民航大學 電子信息與自動化學院，天津300300）

航空發動機作為航空器的核心設備，其系統的結構復雜且長期工作于高溫、高壓的惡劣環境下，限制了一些傳統設備對航空發動機尾氣進行監測，基于電容層析成像技術（Electrical Capacitance Tomography，ECT）的航空發動機尾氣檢測系統為解決此問題提供了可能性。

ECT技術的優點在于非侵入性、安全性、響應速度快、安裝便捷及成本低等，使得其在工業過程可視化檢測中成為一項極具發展潛力的技術［1］。

1984年，Brown對電學成像技術進行了綜述，并且提出了反投影算法［2］；1985年，Murai和Kagawa［3］根據Geselowitz靈敏度原理將靈敏度系數法引入到電學成像中；1998年，Vauhkonen等［4］提出了卡爾曼濾波法，并將其應用到電學層析成像中；1999年，楊五強等［5］提出了改進的Landweber迭代算法，并成功應用到電學成像領域；2002年，Kim等［6］研究出一種擴展的卡爾曼濾波算法；2003年，楊五強等［7］對ECT技術中常用的算法進行了綜述；2005年，王化祥等［8］提出了基于靈敏度矩陣的奇異值分解理論及正則化共軛梯度（Conjugate Gradient，CG）法；2007年，王化祥等［9］提出了基于自適應網格的改進總變差算法。

本文將雙共軛梯度（Biconjugate Gradient，BICG）算法引入到ECT技術中。為了得到更好的效果，對該算法做了一些改進，將BICG算法與正則化算法相結合，并與傳統的迭代類算法進行比較［10-11］。

1 ECT系統

1.1 系統組成

傳統的ECT系統主要由三大子系統組成：傳感器系統、數據采集系統和上位機成像系統。3個系統具有各自的功能。傳感器系統主要測量極板間的電容值，數據采集系統將測得的電容值傳輸到上位機，上位機通過軟件及程序進行圖像重建，這就是一個簡單的ECT系統的成像過程。圖1為傳統ECT系統模型。

圖1 傳統ECT系統模型Fig.1 Traditional ECT system model

1.2 算法介紹

CG算法一般適用于求解大型線性方程，屬于最速下降法的一種改進，誤差不計的條件下，CG算法在N步內就可以達到收斂效果。雖然CG算法在有限的步驟內就可到達收斂，但對于求解出的矩陣G來說，N的取值是非常大的。因此，本文在經典CG算法的基礎上引進了BICG算法來提高計算速度。BICG算法使用2組共扼向量作為搜索方向，而CG算法僅使用1組共軛向量作為搜索方向，因此BICG算法的速度得到大幅度提高。如果系數矩陣A為實對稱陣，且求解節點數較多時，BICG算法收斂速度將會更快。通過實驗仿真和數據對比，驗證了BICG算法效果明顯高于其他幾種迭代類算法。圖2為BICG算法的計算流程。

圖2 BICG算法計算流程Fig.2 BICG algorithm calculation flowchart

2 仿真設計

仿真使用COMSOL5.3軟件對12電極的ECT系統進行建模，對系統進行有限元分割，分割網絡設置為64×64，共有3 228個有效單元，通過有效網絡來求解場內的電容值及靈敏度分布情況，即ECT技術的正問題。空場材料的介電常數設為1，屏蔽罩及電極的介電常數為2.2，管道內被測物體的介電常數設為4.2。圖3（a）為電極數目為12的傳統ECT系統模型，圖3（b）為通過有限元分割后得到的模型。

通過有限元方法分割完網格后，對場內進行計算，計算結果可通過靜電場電勢（見圖4（a））及場內電場線分布情況（見圖4（b））呈現出來。

圖3 ECT系統仿真Fig.3 ECT system simulation

圖4 傳感器電勢分布Fig.4 Sensor potential distribution

3 圖像重建算法

實驗通過5種圖像重建算法，即Tikhonov算法、Landweber算法、CG算法、BICG算法及本文改進算法對各模型仿真進行圖像重建［12-13］。

3.1 LBP算法

LBP算法基于以下2個理想條件：①靈敏度系數均勻分布在場內；②電容值與靈敏度矩陣和介質分布有關。

式中：c為電容值矩陣；s為靈敏度矩陣；g為介質分布向量；sT為s的轉置矩陣，用來代替s的逆矩陣。

3.2 Landweber算法

Landweber算法計算過程如下：

把式（3）展開為函數f（g）的最速下降方向，即梯度為

以LBP算法作為圖像初值，Landweber算法為

設λmax為sTs的最大特征值，α≤2λmax保證Landweber算法收斂。

3.3 CG算法

CG算法計算過程為：在式（7）的兩邊乘以S的轉置矩陣ST，得到式（8）。

式中：S為靈敏度矩陣；G為介電常數矩陣；C為電容值。

1）記初始灰度值為向量G0，設為零向量。

2）給定參數令STS＝S′，STC＝C′，得到新的未知方程S′G＝C′。r0＝C′－S′G0，p0＝r0；進行循環迭代計算。

3）當k＝0，1，2，…，n時，進行循環計算。

式中：pk為第k次的搜索方向；rk為殘差向量。

轉到式（11）。

3.4 BICG算法

BICG算法計算過程如下：

1）記初始灰度值為向量G0，設為零向量。

2）給定參數r0＝C′－S′G0，p0＝r0，pt0＝rt0，rt0為r0的共軛矩陣；進行循環迭代計算。

3）當k＝0，1，2，…，n時，進行循環計算。

轉到式（16）。

3.5 Tikhonov算法

正則化算法是以最小二乘準則及平滑準則作為理論依據計算的，目標泛函為

式中：α為正則化參數；L為正則化算子。

對于給定的正則化算子L，選取合適的正則化參數α＞0，最小化方程，可以得到Tikhonov算法的解為

其原理就是引入一個正則化參數，使原來不能求逆的矩陣可以求逆矩陣。設L為單位矩陣，可得標準Tikhonov算法的解為

正則化參數α的取值是至關重要的，直接影響到圖像重建的質量，而且正則化參數α的選擇是很難確定的。本文實驗選取的正則化參數為0.003。

3.6 改進算法

本文改進算法是在BICG算法基礎上加入了正則化的思想。傳統的BICG算法應用時針對不同的流型所要的迭代次數會有很大的變化，為了改進這一不足，加入了正則化的思想。將BICG算法經過Tikhonov正則化后變為

式中：y0＝STx0，x0為正則化向量。選取合適的x0，令y0＝STx0，S′＝STS＋y0yT0，取初始灰度值G0為零向量，r0＝C′－S′G0，p0＝r0，其循環迭代過程與BICG算法相同。

4 實驗與分析

4.1 圖像重建

運行MATLAB 2014a軟件，首先，對建模完成的程序分別進行空場、物場、滿場電容值和靈敏度的計算，并進行歸一化處理；其次，分別通過Tikhonov算法、Landweber算法、CG算法、BICG算法及改進算法等5種算法對各模型進行圖像重建。仿真模型選擇核心流、雙泡流、三泡流、層流和環流。圖像重建結果如圖5所示。

圖5 各模型重建圖像Fig.5 Reconstructed images of each model

4.2 圖像重建分析

仿真實驗中，通過計算重建圖像的圖像錯誤率（IME）和圖像相關系數（CORR）來對重建圖像進行評估，對重建圖像的效果進行評價［14-15］。將Tikhonov算法、Landweber算法、CG算法、BICG算法及改進算法的重建圖像的圖像錯誤率和圖像相關系數進行對比，分析改進算法的成像優缺點。圖像錯誤率IME和圖像相關系數CORR的計算公式分別為

由式（23）和式（24）計算得IME和CORR，分別如表1和表2所示。

通過數據對比可以看出，使用CG算法得到核心流模型的成像效果最佳，改進算法效果也相差不大，但BICG算法迭代次數少，成像時間快，效率會更高一些。一般使用核心流模型驗證傳感器系統的基本性能。反觀其他模型的圖像重建數據對比，BICG算法和改進算法成像效果優于其他成像算法，尤其在環流模型中，其成像效果遠遠超出其他算法。

表1 IME計算值Tab le 1 Calcu lated values of IM E

實驗中，除了計算圖像錯誤率和相關系數，還對各算法成像時間和迭代算法的迭代次數做了統計。表3為各成像算法數據對比。通過多次調節對比，各迭代算法選取的迭代次數不同，選擇各算法最優成像效果。通過選擇的迭代次數進行成像時長對比。由數據可以直接看出，改進算法迭代次數相對其他迭代類算法大大減少，而且成像時間也大大縮短。

表2 CORR計算值Table 2 Calculated values of CORR

表3 成像算法數據對比Table 3 Com parison of im aging algorithm data

5 結 論

1）本文提出了將BICG算法應用到ECT技術，而且圖像重建的分辨率更高，誤差減小，通過仿真實驗可以發現BICG算法的可行性，不同的流型都具有較好的成像效果，得到更高質量的重建圖像，大幅度縮短重建圖像的成像時間，可以更好地應用于在線監測。

2）針對不同流型所對應的參數和迭代次數有所不同，而且變化幅度較大，在實際應用中可能會受到一些約束，為此將BICG算法與正則化算法相結合，通過調節一個正則化參數便可得到較好的成像效果。

3）為使改進算法在ECT技術和其他領域可以得到更廣泛、更有效的應用，還需要對其進行進一步的研究和發展。