馬敏，范廣永，孫穎

（中國民航大學(xué) 電子信息與自動(dòng)化學(xué)院，天津300300）

航空發(fā)動(dòng)機(jī)作為航空器的核心設(shè)備，其系統(tǒng)的結(jié)構(gòu)復(fù)雜且長期工作于高溫、高壓的惡劣環(huán)境下，限制了一些傳統(tǒng)設(shè)備對(duì)航空發(fā)動(dòng)機(jī)尾氣進(jìn)行監(jiān)測(cè)，基于電容層析成像技術(shù)（Electrical Capacitance Tomography，ECT）的航空發(fā)動(dòng)機(jī)尾氣檢測(cè)系統(tǒng)為解決此問題提供了可能性。

ECT技術(shù)的優(yōu)點(diǎn)在于非侵入性、安全性、響應(yīng)速度快、安裝便捷及成本低等，使得其在工業(yè)過程可視化檢測(cè)中成為一項(xiàng)極具發(fā)展?jié)摿Φ募夹g(shù)［1］。

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

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

1 ECT系統(tǒng)

1.1 系統(tǒng)組成

傳統(tǒng)的ECT系統(tǒng)主要由三大子系統(tǒng)組成：傳感器系統(tǒng)、數(shù)據(jù)采集系統(tǒng)和上位機(jī)成像系統(tǒng)。3個(gè)系統(tǒng)具有各自的功能。傳感器系統(tǒng)主要測(cè)量極板間的電容值，數(shù)據(jù)采集系統(tǒng)將測(cè)得的電容值傳輸?shù)缴衔粰C(jī)，上位機(jī)通過軟件及程序進(jìn)行圖像重建，這就是一個(gè)簡(jiǎn)單的ECT系統(tǒng)的成像過程。圖1為傳統(tǒng)ECT系統(tǒng)模型。

圖1 傳統(tǒng)ECT系統(tǒng)模型Fig.1 Traditional ECT system model

1.2 算法介紹

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

圖2 BICG算法計(jì)算流程Fig.2 BICG algorithm calculation flowchart

2 仿真設(shè)計(jì)

仿真使用COMSOL5.3軟件對(duì)12電極的ECT系統(tǒng)進(jìn)行建模，對(duì)系統(tǒng)進(jìn)行有限元分割，分割網(wǎng)絡(luò)設(shè)置為64×64，共有3 228個(gè)有效單元，通過有效網(wǎng)絡(luò)來求解場(chǎng)內(nèi)的電容值及靈敏度分布情況，即ECT技術(shù)的正問題。空?qǐng)霾牧系慕殡姵?shù)設(shè)為1，屏蔽罩及電極的介電常數(shù)為2.2，管道內(nèi)被測(cè)物體的介電常數(shù)設(shè)為4.2。圖3（a）為電極數(shù)目為12的傳統(tǒng)ECT系統(tǒng)模型，圖3（b）為通過有限元分割后得到的模型。

通過有限元方法分割完網(wǎng)格后，對(duì)場(chǎng)內(nèi)進(jìn)行計(jì)算，計(jì)算結(jié)果可通過靜電場(chǎng)電勢(shì)（見圖4（a））及場(chǎng)內(nèi)電場(chǎng)線分布情況（見圖4（b））呈現(xiàn)出來。

圖3 ECT系統(tǒng)仿真Fig.3 ECT system simulation

圖4 傳感器電勢(shì)分布Fig.4 Sensor potential distribution

3 圖像重建算法

實(shí)驗(yàn)通過5種圖像重建算法，即Tikhonov算法、Landweber算法、CG算法、BICG算法及本文改進(jìn)算法對(duì)各模型仿真進(jìn)行圖像重建［12-13］。

3.1 LBP算法

LBP算法基于以下2個(gè)理想條件：①靈敏度系數(shù)均勻分布在場(chǎng)內(nèi)；②電容值與靈敏度矩陣和介質(zhì)分布有關(guān)。

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

3.2 Landweber算法

Landweber算法計(jì)算過程如下：

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

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

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

3.3 CG算法

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

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

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

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

3）當(dāng)k＝0，1，2，…，n時(shí)，進(jìn)行循環(huán)計(jì)算。

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

轉(zhuǎn)到式（11）。

3.4 BICG算法

BICG算法計(jì)算過程如下：

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

2）給定參數(shù)r0＝C′－S′G0，p0＝r0，pt0＝rt0，rt0為r0的共軛矩陣；進(jìn)行循環(huán)迭代計(jì)算。

3）當(dāng)k＝0，1，2，…，n時(shí)，進(jìn)行循環(huán)計(jì)算。

轉(zhuǎn)到式（16）。

3.5 Tikhonov算法

正則化算法是以最小二乘準(zhǔn)則及平滑準(zhǔn)則作為理論依據(jù)計(jì)算的，目標(biāo)泛函為

式中：α為正則化參數(shù)；L為正則化算子。

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

其原理就是引入一個(gè)正則化參數(shù)，使原來不能求逆的矩陣可以求逆矩陣。設(shè)L為單位矩陣，可得標(biāo)準(zhǔn)Tikhonov算法的解為

正則化參數(shù)α的取值是至關(guān)重要的，直接影響到圖像重建的質(zhì)量，而且正則化參數(shù)α的選擇是很難確定的。本文實(shí)驗(yàn)選取的正則化參數(shù)為0.003。

3.6 改進(jìn)算法

本文改進(jìn)算法是在BICG算法基礎(chǔ)上加入了正則化的思想。傳統(tǒng)的BICG算法應(yīng)用時(shí)針對(duì)不同的流型所要的迭代次數(shù)會(huì)有很大的變化，為了改進(jìn)這一不足，加入了正則化的思想。將BICG算法經(jīng)過Tikhonov正則化后變?yōu)?/p>

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

4 實(shí)驗(yàn)與分析

4.1 圖像重建

運(yùn)行MATLAB 2014a軟件，首先，對(duì)建模完成的程序分別進(jìn)行空?qǐng)觥⑽飯?chǎng)、滿場(chǎng)電容值和靈敏度的計(jì)算，并進(jìn)行歸一化處理；其次，分別通過Tikhonov算法、Landweber算法、CG算法、BICG算法及改進(jìn)算法等5種算法對(duì)各模型進(jìn)行圖像重建。仿真模型選擇核心流、雙泡流、三泡流、層流和環(huán)流。圖像重建結(jié)果如圖5所示。

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

4.2 圖像重建分析

仿真實(shí)驗(yàn)中，通過計(jì)算重建圖像的圖像錯(cuò)誤率（IME）和圖像相關(guān)系數(shù)（CORR）來對(duì)重建圖像進(jìn)行評(píng)估，對(duì)重建圖像的效果進(jìn)行評(píng)價(jià)［14-15］。將Tikhonov算法、Landweber算法、CG算法、BICG算法及改進(jìn)算法的重建圖像的圖像錯(cuò)誤率和圖像相關(guān)系數(shù)進(jìn)行對(duì)比，分析改進(jìn)算法的成像優(yōu)缺點(diǎn)。圖像錯(cuò)誤率IME和圖像相關(guān)系數(shù)CORR的計(jì)算公式分別為

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

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

表1 IME計(jì)算值Tab le 1 Calcu lated values of IM E

實(shí)驗(yàn)中，除了計(jì)算圖像錯(cuò)誤率和相關(guān)系數(shù)，還對(duì)各算法成像時(shí)間和迭代算法的迭代次數(shù)做了統(tǒng)計(jì)。表3為各成像算法數(shù)據(jù)對(duì)比。通過多次調(diào)節(jié)對(duì)比，各迭代算法選取的迭代次數(shù)不同，選擇各算法最優(yōu)成像效果。通過選擇的迭代次數(shù)進(jìn)行成像時(shí)長對(duì)比。由數(shù)據(jù)可以直接看出，改進(jìn)算法迭代次數(shù)相對(duì)其他迭代類算法大大減少，而且成像時(shí)間也大大縮短。

表2 CORR計(jì)算值Table 2 Calculated values of CORR

表3 成像算法數(shù)據(jù)對(duì)比Table 3 Com parison of im aging algorithm data

5 結(jié) 論

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

2）針對(duì)不同流型所對(duì)應(yīng)的參數(shù)和迭代次數(shù)有所不同，而且變化幅度較大，在實(shí)際應(yīng)用中可能會(huì)受到一些約束，為此將BICG算法與正則化算法相結(jié)合，通過調(diào)節(jié)一個(gè)正則化參數(shù)便可得到較好的成像效果。

3）為使改進(jìn)算法在ECT技術(shù)和其他領(lǐng)域可以得到更廣泛、更有效的應(yīng)用，還需要對(duì)其進(jìn)行進(jìn)一步的研究和發(fā)展。