基于QR分解的對稱共軛梯度法成像算法*

馬世文，王化祥

(天津大學電氣與自動化工程學院，天津 300072)

電阻層析成像技術(Electrical Resistance Tomography，ERT)［1］是近年來發展起來的新型測量技術，由于該技術具有非侵入、無輻射、響應快、結構簡單以及成本低廉等優點，在醫學臨床監護和工業測量等領域具有廣闊的應用前景［2－4］。

如圖1所示，電阻層析成像系統由三個部分構成:(1)獲取被測物場信息的空間敏感陣列電極［5］，在交流電壓/電流激勵下，形成一個可從不同角度掃描被測對象的空間敏感場，物場內部介質分布或結構的運動變化對敏感場產生調制作用，使敏感陣列輸出相應的信號;(2)數據采集與處理單元，實時快速地采集空間敏感陣列輸出的反映被測物場介質分布狀態的大量測量數據，完成相應的解調、濾波處理;(3)圖像重建與物場特征參數提取單元，運用圖像重建算法，依據處理后的數據，獲得被測物場的圖像及其變化的時間歷程，實現被測物場內部介質分布可視化，并提取出與介質分布相對應的被測物場的特征信息。

圖1 ERT系統示意圖

為了便于ERT逆問題的求解，通過靈敏場的分析可將ERT的非線性物理模型線性化并歸一化［6］，得出近似的線性基本方程如下

式中，z為n×1階的測量數據矢量矩陣;S為n×m階的靈敏度系數矩陣(歸一化后的雅可比矩陣);g為m×1階的灰度矢量矩陣;n代表獨立測量電極對數;m代表場域的剖分單元數。

根據Geselowitz Sensitivity理論，ERT系統的近似敏感場分布如式［7－8］

其中，φi和φj分別為第i極板和第j極板施加電壓V，而其余極板接地的電勢分布。

共軛梯度法是一種基于最優搜索方向的迭代算法，同其他算法相比，它的收斂速度快，成像質量高。本文在共軛梯度法的基礎上研究了對稱共軛梯度法，進而通過QR分解實現了對稱共軛梯度法的單步成像法。

1 共軛梯度法

共軛梯度法(Conjugate Gradient，CG)是共軛方向類算法中極為重要的一種方法。CG方法適合于系數矩陣為對稱正定的情形。如果不考慮舍入誤差，理論上它能保證最多迭代n步(n為方程組的階數)，便可求得方程組Ax=b的精確解［9－10］。

定義:設A為一個N×N階對稱正定矩陣，p，q為兩個N維向量，如果

則稱向量p和q為相互A共軛。

定理:設A為一N×N階對稱正定矩陣，p1，…，pN為A共軛的N維向量，則此向量系必為線性獨立［11］。

上述定理給出了A共軛向量系的線性獨立性，而共軛梯度法的思想是只要依次沿著n個非零共軛方向進行搜索，便能達到其極小點，因此關鍵是如何構造A的共軛搜索方向。

在ERT圖像重建中，由于靈敏度矩陣S既非正定又不對稱，為保證重建算法具有較好的收斂性，必須對靈敏度矩陣進行正規化。

在式(1)兩端同乘以靈敏度矩陣的轉置ST，則

令 S'=STS，z'=STz。

取初始灰度值 g0為零向量，r0=z'－S'g0，p0=r0，依次按如下步驟循環迭代:

式中，pk為第k次修正方向，而各修正方向之間滿足共軛關系

2 對稱共軛梯度法SCG

ERT逆問題的求解具有不適定性，因此為得到逆問題的有效解，需補充解的附加約束條件以及先驗知識。Tikhonov方法［12］是一種有效解決病態逆問題的方法，它已廣泛應用于ERT的圖像重建。

式(4)經過Tikhonov正則化后可轉變為

式中:α為正則化因子，本文取0.01;I為同維單位矩陣。

本文在Tikhonov正則化的基礎上提出了對稱共軛梯度法［13］。對稱共軛梯度法的每一步的迭代矩陣為對稱矩陣。其具體實現過程如下:

令 S'=STS+αI，z'=STz，任取 g0，令 r0=z'－S'g0，p0=r0，迭代步驟如下:

將式(19)整理，有

3 基于QR分解的對稱共軛梯度法

為了快速得到高質量的重建圖像，在此提出一種新型的基于QR分解的單步對稱共軛梯度成像法RQSCG。系數矩陣STS+αI為非奇異的矩陣，則存在唯一的正交矩陣Q和上三角矩陣R［14－15］，使得

經過QR分解后，式(10)可等價變化為

經過系數矩陣的QR分解后，則可利用對稱共軛梯度法對式(23)進行求解。獲得滿意解g*RQ后，利用g*=QTg*RQ獲得需要的灰度值向量。其中系數矩陣的QR分解可離線計算。

使用對稱共軛梯度法對方程組求解時，設方程迭代解的初始值g(0)RQ，經過k步迭代后的解為g(k)RQ。設系數矩陣SRQ的最大特征值為β，最小特征值為α。根據共軛梯度法的求解原理和數學歸納法可以得出方程經過k步迭代后解的收斂性公式如下［16］

根據QR分解的原理可以得出SRQ≈I，所以β≈α，因此經過單步計算方程解得收斂性公式為

由上式可知，方程在經過單步計算后必收斂，從而理論上證明了基于QR分解的對稱共軛梯度法的單步成像的合理性。

4 實驗仿真與結果分析

為判定算法的有效性，設計停止迭代誤差為

設定停止迭代誤差限值為1×10－6，則共軛梯度法、對稱共軛梯度法和基于QR分解的對稱共軛梯度法的成像誤差和迭代次數以及達到停止迭代誤差的時間如表1、表2、表3所示。

表1 三種算法的迭代步數

表2 三種算法的迭代誤差

表3 三種算法的迭代時間 單位:ms

由表1、2、3可見，對稱共軛梯度算法與共軛梯度算法相比，明顯減少了迭代次數，具有更快的收斂速度。而基于QR分解的對稱共軛梯度法，僅需要一步迭代，其迭代誤差的數量級便可達到10－15，遠優于本文所設定的停止迭代誤差限值。

仿真實驗中，在Intel Core 2 CPU，2GB內存的PC機上，利用COMSOL軟件和Matlab軟件，采用16電極仿真模型，管道內徑為200.5 mm，管壁厚度6.5 mm。三種流型的管內剖分單元數分別為2 148、3 460、1 584。油井內的動態情況通常是油、氣、水混合流體在鋼質套管內流動，井內流體來自于地層孔間縫隙，地層油和天然氣的電導率一般為10－9S·m－1～10－16S·m－1，而地層水由于絕大部分含有一定濃度的鹽離子，其電導率一般為 50 S·m－1～0.5 S·m－1。仿真實驗中介質電導率設為35.869 S·m－1。仿真原圖如表4所示，仿真結果如表5所示，其中共軛梯度法和對稱共軛梯度法各迭代500次，基于QR分解的對稱共軛梯度法迭代1次完成。

表4 仿真原型

表5 三種算法的仿真結果

5 結論

本文提出的SCG算法明顯提高了成像實時性，改善了重建圖像質量。在此基礎上提出RQSCG，實現了一步迭代成像。并對多相流ERT系統進行圖像重建，其成像質量和重建速度與共軛梯度法算法進行比較。仿真實驗結果表明，SCG算法、RQSCG算法對重建圖像質量有明顯的改進，而RQSCG算法的實時性明顯優于CG和SCG算法。

［1］馬藝馨，徐苓安，姜常珍.電阻層析成像技術的研究［J］.儀器儀表學報，2001，22(2):195－198，213.

［2］BrownB H.MedicalImpedance Tomography and Process Impedance Tomography:a Brief Review［J］.Meas.Sci.Technol.，2001，12:991－996.

［3］Brown B H，Barber D C.Electrical Impedance Tomography:the Construction and Application of Electrical Impedance Images［J］.Medical Progress through Technology，1987，13:69－75.

［4］張雪輝，王化祥.電容層析成像數字化系統設計［J］.傳感技術學報，2007，20(8):1826－1830.

［5］孫青，王化祥.一種新型網絲傳感器優化設計［J］.傳感技術學報，2010，23(4):465－470.

［6］Yang Wuqiang，Peng Lihui.Image Reconstruction Algorithms for Electrical Capacitance Tomography［J］.Measurement Science and Technology，2003，14:1－13.

［7］唐磊.電學成像系統圖像重建算法研究及可視化軟件設計［D］.天津:天津大學，2008

［8］Vauhkonen M，Vadasz D，Kaipio J P，et al.Electrical Impedance Tomography with Basis Constraints［M］.Inverse Problems，1997，13:523－530.

［9］Hanke M.Conjugate Gradient Type Methods for Ill Posed Problems［M］.Scientific ＆ Technical，Longman，Harlow，1995.

［10］Rajen Manicon Murugan.An Improved Electrical Impedance Tomography(EIT)Algorithm for the Detection and Diagnosis of Early Stages ofBreastCancer［D］.PHD thesis，University of Manitoba，1999.

［11］鄧乃揚.無約束最優化計算方法［M］.北京:科學出版社，1982.

［12］WH Winston.Tikhonov A N and Arsenin V Y，Solutions of illposed problems［M］.V.H.Winston ＆ Sons，a Division of Scripta Technica，Inc，1977.

［13］劉東毅，邵琛.下降對稱的 Polak-Ribiere-Polyak共軛梯度法［J］.天津大學學報，2010，43(4):367－372.

［14］戴華.矩陣論［M］.北京:科學出版社，2001

［15］和斌濤.初等變換與矩陣的QR分解的關系［J］.科學技術與工程，2009，9(21):6484－6485.

［16］李愛芹.線性方程組的迭代解法［J］.科學技術與工程，2007，7(14):3357－3364.