基于L曲線(xiàn)法的超聲CT正則化參數(shù)優(yōu)化

張 培，呂晶晶，李 媛

（中北大學(xué)儀器科學(xué)與動(dòng)態(tài)測(cè)試教育部重點(diǎn)實(shí)驗(yàn)室，山西太原 030051）

通過(guò)矩量法將上述兩個(gè)卷積方程轉(zhuǎn)化為代數(shù)方程組的形式:

方程（3）、（4）僅描述了物體在某一方向入射波作用下接受到的散射波，為獲得足夠多的信息，需要從K個(gè)方向發(fā)射超聲波，由此得到K×M個(gè)方程，為使上述方程組可解，通

0 引言

80年代初，Johnson將微波技術(shù)中矩量法引入超聲成像領(lǐng)域［1］，矩量法不再局限于Born近似的限制，可以方便地引入被測(cè)物體的某些先驗(yàn)知識(shí)（這些先驗(yàn)知識(shí)是克服任何逆散射問(wèn)題的不適定性和對(duì)噪聲敏感性的關(guān)鍵），可以對(duì)大物體、高對(duì)比度、強(qiáng)散射的物體進(jìn)行成像，因此適應(yīng)范圍更廣［2］。

不適定性是圖像重建逆問(wèn)題求解中非常普遍的情況，本文采用Tikhonov和TSVD正則化方法解決其不適定性，兩種正則化算法的成功應(yīng)用依賴(lài)于正則化參數(shù)λ的合適選取。其中容差原理法確定正則化參數(shù)需要知道數(shù)據(jù)項(xiàng)中噪聲大小，而廣義交叉驗(yàn)證法操作較為繁瑣且結(jié)果不穩(wěn)定。本文基于L-曲線(xiàn)法，提出以解的范數(shù)和殘差變化量的加權(quán)形式作為確定正則化參數(shù)的依據(jù)，在迭代過(guò)程根據(jù)問(wèn)題不適定性程度，自適應(yīng)地調(diào)整搜索范圍。通過(guò)兩種正則化技術(shù)和優(yōu)化策略來(lái)實(shí)現(xiàn)在空間域內(nèi)的圖像重構(gòu)。

1 UCT圖像重建原理

當(dāng)用超聲波從某一方向作用于物體的被成像截面時(shí)，超聲波與物體內(nèi)部介質(zhì)相互作用產(chǎn)生的全場(chǎng)滿(mǎn)足非齊次亥姆霍茲方程［3］，其解可用Lippmann-Schwinger積分方程表示:

式中:p（→）和pin分別為全場(chǎng)和入射場(chǎng);G（→-） 為自由空間的格林函數(shù);O（=k20［n2（→）-1］為物體內(nèi)部介質(zhì)聲學(xué)特性的未知函數(shù)。

在實(shí)際應(yīng)用中，我們只能測(cè)到物體的散射場(chǎng)和入射場(chǎng)，將式（1）轉(zhuǎn)化為散射場(chǎng)形式:常KM＞N，即要求方程組為超定的。

通過(guò)矩量法將上述兩個(gè)卷積方程轉(zhuǎn)化為代數(shù)方程組的形式:

方程（3）、（4）僅描述了物體在某一方向入射波作用下接受到的散射波，為獲得足夠多的信息，需要從K個(gè)方向發(fā)射超聲波，由此得到K×M個(gè)方程，為使上述方程組可解，通

這個(gè)超定方程組的求解具有較強(qiáng)的非線(xiàn)性。這里采用Born迭代方法［4］。具體描述如下:

①先假設(shè)全場(chǎng)等于入射場(chǎng)P（t）=P（in），由方程（4）求得未知函數(shù)的初始解O0，常稱(chēng)其為Born逆解;

②然后將Ok代入方程（3）求更接近實(shí)際的物體內(nèi)部的全場(chǎng)P（kt），k表示迭代次數(shù);

③由第②步求得的全場(chǎng)P（kt），代入方程（4），計(jì)算散射場(chǎng)P（ks）=DOP（kt），并求計(jì)算散射場(chǎng)與測(cè)量散射場(chǎng)的差ΔP（ks）=P（ks）-P（s），若‖ΔP（ks）‖小于事先給定的δ，則停止迭代，否則，轉(zhuǎn)④;

④根據(jù)散射場(chǎng)的改變量 ΔP（ks），由方程 ΔP（ks）=DΔOkP（kt求未知函數(shù)的改變量ΔOk，然后賦Ok+1=Ok+ΔOk作為未知函數(shù)的新值;

⑤由新的未知函數(shù)Ok+1（rn），返回步驟②，進(jìn)一步求更近似的全場(chǎng)。

在迭代過(guò)程中，將遇到不適定性逆散射方程的求解問(wèn)題，而對(duì)不適定性問(wèn)題的正則化技術(shù)的選擇是否得當(dāng)，將直接影響迭代算法的穩(wěn)定性。所以，正則化技術(shù)是圖像重建問(wèn)題的關(guān)鍵。

2 正則化方法及參數(shù)的選取

算法1:用Tikhonov方法計(jì)算Ax=b的算法如下［6］。

①矩陣A的奇異值分解:［U，S，V］=svd（A），其中S=diag（σ1，σ2，…σn）

整數(shù)λ稱(chēng)為正則化參數(shù)。

③作圖（lg‖xλ‖2），lg‖Axλ-b‖2）），確定曲線(xiàn)拐角點(diǎn)，拐點(diǎn)對(duì)應(yīng)的xλ是正則解xtik。

算法2:用TSVD方法計(jì)算Ax=b的算法如下［7］。

①矩陣 A 的奇異值分解:［U，S，V］=svd（A）。

②對(duì)1≤k≤T，計(jì)算解的范數(shù):

③作圖（lg（‖xλ‖2），lg（‖Axλ-b‖2）），確定曲線(xiàn)拐角點(diǎn)，拐點(diǎn)對(duì)應(yīng)的xk是正則解xtsvd。

研究表明，對(duì)離散不適定問(wèn)題，lg（‖xk‖2），lg（‖Axk-b‖2）的關(guān)系曲線(xiàn)總是呈現(xiàn)L-曲線(xiàn)形狀，該關(guān)系曲線(xiàn)有一個(gè)明顯的拐角，最優(yōu)正則化參數(shù)位于L-曲線(xiàn)拐點(diǎn)附近［8］。因此，可通過(guò)尋找L-曲線(xiàn)拐點(diǎn)作為最優(yōu)正則化參數(shù)。

3 重建結(jié)果比較

實(shí)驗(yàn)裝置我們采用圓環(huán)形結(jié)構(gòu)，發(fā)射器（同時(shí)也是接收器）等間隔地放在圍繞物體的圓環(huán)上，半徑取20×λ，背景介質(zhì)為水。采用F=200 kHz的超聲波作為入射波，超聲波在水中的波速約為C0=1 500 m/s，波長(zhǎng)約為7.5 mm，對(duì)應(yīng)的波數(shù)為k0=0.837 758 mm-1。按照矩量法，我們對(duì)包含物體的正方形區(qū)域進(jìn)行均勻采樣，采樣間隔均為d=λ/10=0.75 mm，采樣點(diǎn)的個(gè)數(shù)為35×35，劃分后單元的總數(shù)為N=1 225，且在每個(gè)小單元上做內(nèi)接圓，其半徑為a=d/2=0.375 mm。換能器共有50個(gè)，可工作于發(fā)射、接收兩種模式，每次僅有一個(gè)換能器發(fā)射，其他換能器接收來(lái)自物體的散射波。

如圖1所示，實(shí)驗(yàn)選用的樣品如下:對(duì)比度為20%的人工繪制的圖形，輪廓信息較為豐富。對(duì)于已知像函數(shù)O（→r），先由全場(chǎng)方程（3）求出ROI中的全場(chǎng)分布，然后通過(guò)散射場(chǎng)方程（4）求出物體外的散射場(chǎng)P（s）分布，再使用BI進(jìn)行圖像重建。兩種方法的重建結(jié)果如圖2和圖3所示:

圖1 目標(biāo)原始函數(shù)

圖2 采用Tikhonov迭代4次結(jié)果

圖3 TSVD正則化方法迭代18次的結(jié)果

圖4 Tikhonov方法相對(duì)誤差曲線(xiàn)

由圖4可以看出采用Tikhonov方法時(shí)，相對(duì)誤差在第4次相對(duì)誤差最小，當(dāng)超過(guò)第4次后，相對(duì)誤差變大。圖5可以看出采用TSVD方法迭代到第18次時(shí)最小誤差恒定不變，因此在第18次中止迭代。對(duì)比兩種方法，TSVD方法明顯好于的Tikhonov方法重建的結(jié)果。

圖5 TSVD方法相對(duì)誤差曲線(xiàn)

4 結(jié)論

本文運(yùn)用Born迭代算法解決了方程的非線(xiàn)性問(wèn)題，對(duì)于離散不適定問(wèn)題，采用Tikhonov和TSVD兩種正則化方法分別對(duì)方程進(jìn)行了求解。經(jīng)實(shí)驗(yàn)驗(yàn)證:TSVD正則化方法明顯優(yōu)于Tikhonov正則化方法。超聲逆散射成像技術(shù)是一個(gè)比較復(fù)雜的系統(tǒng)工程，算法的穩(wěn)定性及速度只是問(wèn)題的一個(gè)方面。目前尚沒(méi)有成熟的理論和技術(shù)可用于產(chǎn)品的自動(dòng)檢測(cè)，對(duì)該方面的一些關(guān)鍵技術(shù)進(jìn)行研究具有一定的學(xué)術(shù)價(jià)值和廣闊的應(yīng)用前景。

［1］Tracy M L，Johnson S A.Inverse Scattering Solutions by a Sinc Basis，Multiple Source，Moment Method-Part II:Numerical Evaluations［J］.Ultrasonic Imaging，1983，5（4）:376-392.

［2］陶進(jìn)緒，張東文，葉寒生.頻域法超聲逆散射成像［J］.信號(hào)處理，2009，25（2）:169-173.

［3］劉超，汪元美.超聲層析成像的理論與實(shí)現(xiàn)［D］.浙江大學(xué)，生物醫(yī)學(xué)工程，博士論文，2003.

［4］Chew W C，Wang Y M.Reconstruction of Two-Dimensional Permittivity Distribution Using the Distorted Born Iterative Method［J］.IEEE Transactions on Medical Imaging，1990，9（2）:218-225.

［5］韓旭，劉杰，李偉杰.時(shí)域內(nèi)多源動(dòng)態(tài)載荷的一種計(jì)算反求技術(shù)［J］.力學(xué)學(xué)報(bào)，2009，41（4）:595-601.

［6］Hansen P C，O＇Leary D P.The Use of the L-Curve in the Regularization of Discrete Ill-Posed Problems［J］.SIAM Journal on Scientific Computing，1993，14（6）:1487-1503.

［7］Hansen P C，Jensen T K.An Adaptive Pruning Algorithm for the Discrete L-curve Criterion［J］.Journal of Computational and Applied Mathematics，2007，198（2）:483-492.

［8］王化祥，何永勃，朱學(xué)明.基于L曲線(xiàn)法的電容層析成像正則化參數(shù)優(yōu)化［J］.天津大學(xué)學(xué)報(bào)，2006，39（3）:306-309.