本發(fā)明屬于電學(xué)層析成像技術(shù)領(lǐng)域，涉及利用稀疏正則化方法實(shí)現(xiàn)圖像重建及閾值迭代求解方法。

背景技術(shù)：

電學(xué)層析成像技術(shù)(Electrical Tomography，ET)出現(xiàn)于20世紀(jì)80年代后期，是一種基于電特性敏感機(jī)理的過程層析成像技術(shù)，它基于邊界測量值對測量區(qū)域介質(zhì)的電特性(電導(dǎo)率/介電系數(shù)/復(fù)導(dǎo)納/磁導(dǎo)率)分布信息進(jìn)行成像，進(jìn)而得到介質(zhì)的分布信息。電學(xué)層析成像技術(shù)主要包括電阻層析成像(Electrical Resistance Tomography，ERT)、電容層析成像(Electrical Capacitance Tomography，ECT)、電阻抗層析成像(Electrical Impedance Tomography，EIT)和電磁層析成像(Electrical Magnetic Tomography，EMT)。電學(xué)層析成像在多相流、地質(zhì)勘探以及醫(yī)學(xué)成像領(lǐng)域有廣泛的應(yīng)用前景，可以實(shí)現(xiàn)長期、持續(xù)監(jiān)測，實(shí)現(xiàn)功能性成像。

通過線性化處理，可以將電學(xué)層析成像中逆問題轉(zhuǎn)化為線性逆問題并基于靈敏度矩陣來求解。針對逆問題求解的病態(tài)性，可選取正則化方法處理。正則化算法通用目標(biāo)函數(shù)為：

其中，λ為正則化參數(shù)，用于調(diào)節(jié)誤差項(xiàng)與正則化項(xiàng)在目標(biāo)函數(shù)中的比重，G(x)為正則化項(xiàng)，根據(jù)具體應(yīng)用的需要可以設(shè)計(jì)為不同的形式。在電學(xué)層析成像領(lǐng)域中，經(jīng)典的Tikhonov算法，即以解的2范數(shù)作為正則化項(xiàng)的L₂正則化算法，目標(biāo)函數(shù)為是由Vauhkonen等人在1998年發(fā)表于《IEEE醫(yī)學(xué)成像》(IEEE Transactions on Medical Imaging)第17卷，第285-293頁，題為《基于電阻抗層析成像的Tikhonov正則化及先驗(yàn)信息選擇》(Tikhonov regularization and prior information in electrical impedance tomography)的文章中；近幾年，隨著壓縮傳感理論的發(fā)展，稀疏正則化算法重新受到多個領(lǐng)域的重視，成為了研究的熱點(diǎn)，其中以解的1范數(shù)作為正則化項(xiàng)的L₁正則化算法應(yīng)用最為廣泛，其目標(biāo)函數(shù)為：由Jin Bangti等人在2012年發(fā)表于《工程中的數(shù)值計(jì)算》(International Journal For Numerical Methods In Engineering)第89卷，第337-353頁，題為《基于稀疏正則化的電阻抗層析成像重建算法》(A reconstruction algorithm for electrical impedance tomography based on sparsity regularization)的文章中。采用L₂正則化方法求解逆問題所得的解會出現(xiàn)過光滑現(xiàn)象，重建圖像具有較大的偽影；而L₁正則化方法在圖像的邊界銳化和去偽影方面相比L₂正則化優(yōu)勢明顯，但由于其目標(biāo)函數(shù)的不可微性，求解起來相對困難，并且其對具有光滑物體分布的場域所求解時會出現(xiàn)過稀疏的問題，導(dǎo)致重建圖像小于場域物體的真實(shí)大小。但目前L₁正則化算法已成為稀疏化重建算法中最受歡迎的方法之一，不斷有新的求解算法出現(xiàn)，其中凸優(yōu)化類算法在求解L₁正則化問題中應(yīng)用最為廣泛。

凸優(yōu)化算法是將目標(biāo)函數(shù)的求解問題轉(zhuǎn)化為凸優(yōu)化問題來進(jìn)行求解的方法。常見的求解方法包括迭代閾值法、內(nèi)點(diǎn)法、映射梯度法等。其中迭代閾值算法由于其計(jì)算簡便、原理清晰而成為研究熱點(diǎn)，但同時也以迭代速度慢著稱。常用的迭代閾值算法包括迭代軟閾值算法、迭代硬閾值算法以及一些基于迭代閾值算法的改進(jìn)算法等。Amir Beck等人于2009年發(fā)表于《工業(yè)和應(yīng)用數(shù)學(xué)學(xué)會雜志：成像科學(xué)》(SIAM Journal on Imaging Sciences)第二卷，第183-202頁，題為《一種線性逆問題求解的快速迭代收縮閾值方法》(A fast iterative shrinkage-thresholding algorithm for linear inverse problems)提供了一種加速的軟閾值迭代方法來求解L₁正則化問題，提高了求解速度。

張玲玲等人于2011年在發(fā)表于《世界科學(xué)與工程學(xué)術(shù)學(xué)會：電路和系統(tǒng)處理》(WSEAS Transactions on Circuits and Systems)第10卷，第393-402頁，題為《一種用于電阻層析成像的快速迭代收縮閾值方法》(A Fast Iterative Shrinkage-Thresholding Algorithm for Electrical Resistance Tomography)的文章將該快速迭代閾值算法應(yīng)用于電阻層析成像中，并對成像結(jié)果進(jìn)行討論，驗(yàn)證了方法的有效性。

然而在現(xiàn)有研究中，利用基于L₁正則化的軟閾值迭代方法仍然如下缺點(diǎn)：(1)迭代速度慢；(2)重建圖像的效果雖然在去偽影和邊界銳化放面好于L₂正則化方法，但是依然存在所求解過于稀疏、重建圖像小于場域物體的真實(shí)分布大小的問題；(3)正則化參數(shù)的選取與場域分布密切相關(guān)，當(dāng)物體分布發(fā)生變化時，通常需要重新計(jì)算選擇更合適的正則化參數(shù)才能保證所求解逼近真實(shí)解，而且正則化參數(shù)的選取對圖像重建的質(zhì)量有很大的影響。但在現(xiàn)有研究中，基于迭代軟閾值的L₁正則化方法在進(jìn)行電學(xué)層析成像逆問題求解時，通常在整個迭代過程始終采用一個固定的正則化系數(shù)(λ值)，在實(shí)際應(yīng)用中，對于一些未知量經(jīng)常變化的應(yīng)用場景，重復(fù)計(jì)算正則化參數(shù)將使逆問題求解過程更加復(fù)雜，難以得到一個合適的正則化參數(shù)值。雖然目前很多學(xué)者已經(jīng)針對正則化系數(shù)選取問題做了很多研究，但是并沒有提出一種具有普遍適用性的最佳選取策略，多數(shù)正則化系數(shù)選取方法都只能在某些領(lǐng)域或某些情況下才能夠達(dá)到較好地應(yīng)用效果，導(dǎo)致算法的應(yīng)用范圍較小。所以研究的重點(diǎn)在于實(shí)現(xiàn)針對不同場域分布的正則化系數(shù)自適應(yīng)選取，在每一步迭代中自適應(yīng)地調(diào)節(jié)正則化參數(shù)，在一定范圍內(nèi)保證圖像重建的質(zhì)量不受初始正則化參數(shù)的影響，并且使重建圖像更接近于實(shí)物分布。

技術(shù)實(shí)現(xiàn)要素：

本發(fā)明的目的在于，基于一種新的閾值迭代方法，提出一個新的稀疏正則化目標(biāo)函數(shù)，并且自適應(yīng)地更新正則化參數(shù)，來解決電學(xué)層析成像的稀疏重建問題。本發(fā)明以Firm迭代閾值方法為基礎(chǔ)，求解新的稀疏正則化目標(biāo)函數(shù)，并根據(jù)場域分布，在每一步迭代計(jì)算中自適應(yīng)地更新正則化參數(shù)值，解決稀疏正則化參數(shù)難以選取、求解計(jì)算速度較慢的問題，在不改變正則化參數(shù)的前提下，仍然可以在不同場域分布的圖像重建中得到高質(zhì)量的圖像；同時也解決了L₁正則化圖像重建結(jié)果過于稀疏、重建圖像小于真實(shí)分布的問題，提高電學(xué)層析成像逆問題的求解精度和圖像質(zhì)量。本發(fā)明的技術(shù)方案如下：

一種基于Firm閾值迭代的電學(xué)層析成像稀疏重建方法，用于非連續(xù)介質(zhì)的電學(xué)層析成像，該方法將電學(xué)層析成像問題看作求解一個線性不適定問題Ax＝b，其中A為靈敏度矩陣，b為相對邊界測量值向量，x為需要求解的成像灰度值。

所建立的目標(biāo)函數(shù)為：其中τ為初始人為選擇的正則化系數(shù)，||·||為歐幾里得范數(shù)，ρ是每一步迭代中實(shí)際采用的正則化參數(shù)的大小；是基于Firm閾值迭代函數(shù)的新的稀疏正則化項(xiàng)。重建方法包含有以下步驟：

(1)根據(jù)被測場域，獲取重建所需的相對邊界測量值向量b和靈敏度矩陣A；邊界測量值的獲取，通常是指將被測對象置于電學(xué)層析成像測量系統(tǒng)中，被測場域外均勻分布n個電極，采用電流激勵電壓測量且激勵電極不測量的模式，采集循環(huán)激勵循環(huán)測量下各個電極上的邊界電壓，相對邊界測量值向量b為不含內(nèi)含物的空場邊界測量電壓向量b1和含有內(nèi)含物的有物場的邊界測量電壓向量b2之差。

(2)設(shè)置初始化參數(shù)：解初始值x₀＝A^Tb，稀疏度即非零元素個數(shù)K，初始正則化參數(shù)τ，常數(shù)t₁＝1，初始解x₀，迭代終止條件tol，解向量的元素總數(shù)N。

(3)對x_k進(jìn)行加速預(yù)處理，方式為其中預(yù)處理常數(shù)t₁＝1，t按照在每一步中進(jìn)行迭代更新，在接下來的運(yùn)算中用y_k代替x_k進(jìn)行迭代，以實(shí)現(xiàn)加速迭代求解的作用。

(4)根據(jù)上一循環(huán)中計(jì)算得到的正則化參數(shù)ρ，即本次閾值迭代計(jì)算中所選擇的閾值大小，對經(jīng)過加速處理后的上一步迭代結(jié)果進(jìn)行閾值迭代計(jì)算，代入迭代格式x_k+1＝S_ρ,λ(y_k-A^T(Ay_k-b))中，其中x_k+1＝S_ρ,λ(G(x_k))的閾值迭代函數(shù)算子為：將ρ代入進(jìn)行計(jì)算。

(5)計(jì)算本次迭代得到的解中非零元素的個數(shù)，更新稀疏度K，根據(jù)迭代解的稀疏度更新計(jì)算下一步迭代計(jì)算中使用的正則化參數(shù)ρ；正則化參數(shù)ρ的更新策略為：

(6)判斷迭代是否符合迭代終止條件，若是則迭代終止，進(jìn)行下一步操作；若否，則跳回第(3)步繼續(xù)迭代求解；

(7)根據(jù)最終求解所得灰度值，進(jìn)行成像。

本發(fā)明的有益效果是：基于一種新的稀疏正則化目標(biāo)函數(shù)，在計(jì)算過程中選擇與新目標(biāo)函數(shù)相對應(yīng)的Firm閾值函數(shù)進(jìn)行迭代求解，同時在每一步迭代過程中自適應(yīng)地根據(jù)解稀疏度來更新正則化參數(shù)，是一種新的圖像稀疏重建方法。方法有效地克服了傳統(tǒng)L₁正則化閾值迭代方法成像結(jié)果過于稀疏的缺點(diǎn)，提高了圖像重建的精度；也解決了正則化參數(shù)難以選取的問題，提高了算法的適用性和易用性，擴(kuò)展了稀疏重建方法的應(yīng)用范圍。

附圖說明

圖1為本發(fā)明的一種新的稀疏正則化目標(biāo)函數(shù)的電學(xué)層析成像重建方法的流程框圖；

圖2為本發(fā)明的電阻層析成像系統(tǒng)圓形單截面被測場域及電極分布；

圖3為本發(fā)明的實(shí)例，在選取一個小圓模型的真實(shí)分布時，傳統(tǒng)的基于L₁正則化的軟閾值迭代求解(FISTA)方法和基于新的稀疏正則化目標(biāo)函數(shù)的Firm閾值迭代(FIVTA)方法，分別在一組初始正則化參數(shù)條件下，圖像重建結(jié)果的示意圖，其中基于新的稀疏正則化目標(biāo)函數(shù)的Firm閾值迭代(FIVTA)方法即本專利發(fā)明的方法；

圖4為本發(fā)明的實(shí)例，在選取兩個小圓模型的真實(shí)分布時，傳統(tǒng)的基于L₁正則化的軟閾值迭代求解(FISTA)方法和基于新的稀疏正則化目標(biāo)函數(shù)的Firm閾值迭代(FIVTA)方法，分別在一組初始正則化參數(shù)條件下，圖像重建結(jié)果的示意圖，其中基于新的稀疏正則化目標(biāo)函數(shù)的Firm閾值迭代(FIVTA)方法即本專利發(fā)明的方法；

圖5為本發(fā)明的實(shí)例中，在兩個小圓模型中，傳統(tǒng)的基于L₁正則化的軟閾值迭代求解(FISTA)方法和基于新的稀疏正則化目標(biāo)函數(shù)的Firm閾值迭代(FIVTA)方法，分別在一組初始正則化參數(shù)(τ₁～τ₆)條件下，分別所對應(yīng)的迭代次數(shù)示意圖，其中初始正則化參數(shù)τ₁～τ₆是從大到小的排序；

圖中：

1、被測場域 2、電極

以上成像結(jié)果均在1％的隨機(jī)噪聲條件下測得。

具體實(shí)施方式

結(jié)合附圖和實(shí)施例對本發(fā)明的一種新的稀疏正則化目標(biāo)函數(shù)的電學(xué)層析成像重建方法加以說明。

本發(fā)明的一種新的稀疏正則化目標(biāo)函數(shù)的電學(xué)層析成像重建方法，以稀疏正則化為基礎(chǔ)，針對基于L₁正則化的軟閾值迭代方法所得解過稀疏以及參數(shù)難以選取的問題，結(jié)合一種基于基于新的稀疏正則化目標(biāo)函數(shù)的Firm閾值迭代方法，提出在迭代過程中利用前一步計(jì)算結(jié)果的稀疏度更新下一步迭代中的閾值，直至迭代終止的重建方法。

如圖1所示，為本發(fā)明的一種基于新的稀疏正則化目標(biāo)函數(shù)的電學(xué)層析成像重建方法流程圖。如圖2所示為電學(xué)層析成像之一的電阻層析成像系統(tǒng)圓形單截面被測場域及電極分布，采用16電極均勻分布在場域外壁。選取兩種典型的非連續(xù)介質(zhì)模型為實(shí)施例，場域內(nèi)物體真實(shí)分布如圖3左側(cè)一豎列所示。為了較好地體現(xiàn)本發(fā)明中基于新的稀疏正則化目標(biāo)函數(shù)(采用Firm閾值函數(shù)進(jìn)行迭代計(jì)算并在迭代中自適應(yīng)更新閾值大小)的方法與基于傳統(tǒng)L₁正則化(采用軟閾值函數(shù)進(jìn)行迭代計(jì)算并在迭代中不更新閾值大小)的不同，分別給出這兩種算法在一組不同初始正則化參數(shù)條件下的求解結(jié)果。實(shí)施例包括如下具體步驟：

一種基于Firm閾值迭代的電學(xué)層析成像稀疏重建方法，適用于非連續(xù)介質(zhì)的電學(xué)層析成像，該方法將電學(xué)層析成像問題看作求解一個線性不適定問題Ax＝b，其中A為靈敏度矩陣，b為相對邊界測量值向量，x為需要求解的成像灰度值。

閾值迭代方法是梯度下降方法的一種拓展，其內(nèi)在遵循的原理是一致的，均沿梯度的方向進(jìn)行迭代計(jì)算。和傳統(tǒng)Landweber迭代方法沿目標(biāo)函數(shù)的梯度方向求解相比，由于目標(biāo)函數(shù)中增加了一項(xiàng)正則化項(xiàng)，變成的形式，此時閾值迭代方法可以視為一步Landweber迭代之后再加上一個“閾值衰減”的步驟，沿著新的目標(biāo)函數(shù)的梯度方向進(jìn)行迭代求解。傳統(tǒng)的閾值迭代方法的迭代格式為x_k+1＝S_λ(x_k-tA^T(Ax_k-b))。

這種基于Firm閾值迭代方法的電學(xué)層析成像稀疏重建方法的目標(biāo)函數(shù)為：其中τ為初始人為選擇的正則化系數(shù)，||·||為歐幾里得范數(shù)，ρ是每一步迭代中實(shí)際采用的正則化參數(shù)的大小。是基于Firm閾值迭代函數(shù)得到的新的稀疏正則化項(xiàng)。

重建方法包含有以下步驟：

(1)針對非連續(xù)介質(zhì)的兩個典型模型，分別獲取各自重建所需的邊界測量值和靈敏度矩陣：

邊界測量值是將被測對象置于電學(xué)層析成像測量系統(tǒng)中，被測場域外均勻分布16個電極(如圖2所示)，采用電流激勵電壓測量且激勵電極不測量的模式，采集循環(huán)激勵循環(huán)測量下的邊界電壓，共獲得208個測量值；逆問題右端項(xiàng)b為不含內(nèi)含物的空場邊界電壓b1和含有內(nèi)含物的有物場的邊界測量電壓b2之差(即右端項(xiàng)相對邊界測量值b＝b1-b2)；

靈敏度矩陣是根據(jù)不含內(nèi)含物的空場的邊界測量電壓，結(jié)合靈敏度理論，計(jì)算靈敏度矩陣，計(jì)算公式為：

其中A_ij是第j個電極對對第i個電極對的靈敏度系數(shù)，φ_i,φ_j分別為第i個電極對及第j個電極對在激勵電流為I_i,I_j時場域電勢分布，x,y分布為場域的橫縱坐標(biāo)；

(2)設(shè)置初始化參數(shù)：

設(shè)置的初始化參數(shù)包括：解初始值x₀＝A^Tb，初始預(yù)估稀疏度(非零元素個數(shù))K＝N/5，正則化系數(shù)τ通常從開始從大到小選取以方便兩種算法進(jìn)行比較，t₁＝1，迭代終止條件tol＝2e^-2；

(3)對上一步得到的解x_k進(jìn)行預(yù)處理，目的是為了加速收斂，公式為：

其中，t_k的更新策略是：

(4)對預(yù)處理之后的解進(jìn)行閾值迭代計(jì)算，x_k+1＝S_ρ,λ(y_k-A^T(Ay_k-b))；其中

迭代公式的迭代格式為：

(5)計(jì)算非零元素的個數(shù)，更新稀疏度K，計(jì)算下一步迭代計(jì)算中使用的閾值ρ；閾值ρ的更新策略為：

(6)根據(jù)兩步迭代所求解之間殘差(||x_k+1-x_k||≤tol是否成立)來判斷迭代是否完成，若是則迭代終止，進(jìn)行下一步操作，若否，則跳回第(3)步繼續(xù)求解；

(7)根據(jù)最終求解所得灰度值，進(jìn)行成像。

圖3為本發(fā)明的實(shí)例，在選取一個小圓模型的真實(shí)分布時，傳統(tǒng)的基于L₁正則化的軟閾值迭代求解(FISTA)方法和基于新的稀疏正則化目標(biāo)函數(shù)的Firm閾值迭代(FIVTA)方法，分別在一組初始正則化參數(shù)條件下，圖像重建結(jié)果的示意圖；圖4為本發(fā)明的實(shí)例，在選取兩個小圓模型的真實(shí)分布時，傳統(tǒng)的基于L₁正則化的軟閾值迭代求解(FISTA)方法和基于新的稀疏正則化目標(biāo)函數(shù)的Firm閾值迭代(FIVTA)方法，分別在一組初始正則化參數(shù)條件下，圖像重建結(jié)果的示意圖；圖5為本發(fā)明的實(shí)例中，在兩個小圓模型中，傳統(tǒng)的基于L₁正則化的軟閾值迭代求解(FISTA)方法和基于新的稀疏正則化目標(biāo)函數(shù)的Firm閾值迭代(FIVTA)方法，分別在一組初始正則化參數(shù)(τ₁～τ₆)條件下，分別所對應(yīng)的迭代次數(shù)示意圖，其中初始正則化參數(shù)τ₁～τ₆是從大到小的排序；

可以看出，兩個模型中，在相同的正則化參數(shù)條件下，基于L₁正則化的軟閾值迭代求解(FISTA)方法的解過于稀疏，不能很好地描述物體的大小，相比實(shí)物偏?。换谛碌南∈枵齽t化目標(biāo)函數(shù)的Firm閾值迭代(FIVTA)方法更加準(zhǔn)確的描述了場域內(nèi)物體的分布，提高了逆問題求解精度和圖像重建質(zhì)量。同時，在兩個模型中還可以看出，隨著正則化參數(shù)選擇的不同，F(xiàn)IVTA(本算法)成像的結(jié)果穩(wěn)定性較高，基本不受初始正則化參數(shù)選擇的影響，在實(shí)際應(yīng)用中可以較好地解決正則化參數(shù)難以選取的問題。并且，本算法的迭代次數(shù)和計(jì)算時間相比傳統(tǒng)基于L₁正則化的軟閾值迭代求解(FISTA)方法都有較大的提高。

以上所述實(shí)施例為本發(fā)明的幾個較佳模型，本發(fā)明不局限于該實(shí)施例和附圖所公開的內(nèi)容。凡是不脫離本發(fā)明所公開的精神下完成的等效或修改，都在本發(fā)明保護(hù)的范圍。