張 培,呂晶晶,李 媛
(中北大學(xué)儀器科學(xué)與動態(tài)測試教育部重點實驗室,山西太原 030051)
通過矩量法將上述兩個卷積方程轉(zhuǎn)化為代數(shù)方程組的形式:
方程(3)、(4)僅描述了物體在某一方向入射波作用下接受到的散射波,為獲得足夠多的信息,需要從K個方向發(fā)射超聲波,由此得到K×M個方程,為使上述方程組可解,通
80年代初,Johnson將微波技術(shù)中矩量法引入超聲成像領(lǐng)域[1],矩量法不再局限于Born近似的限制,可以方便地引入被測物體的某些先驗知識(這些先驗知識是克服任何逆散射問題的不適定性和對噪聲敏感性的關(guān)鍵),可以對大物體、高對比度、強散射的物體進行成像,因此適應(yīng)范圍更廣[2]。
不適定性是圖像重建逆問題求解中非常普遍的情況,本文采用Tikhonov和TSVD正則化方法解決其不適定性,兩種正則化算法的成功應(yīng)用依賴于正則化參數(shù)λ的合適選取。其中容差原理法確定正則化參數(shù)需要知道數(shù)據(jù)項中噪聲大小,而廣義交叉驗證法操作較為繁瑣且結(jié)果不穩(wěn)定。本文基于L-曲線法,提出以解的范數(shù)和殘差變化量的加權(quán)形式作為確定正則化參數(shù)的依據(jù),在迭代過程根據(jù)問題不適定性程度,自適應(yīng)地調(diào)整搜索范圍。通過兩種正則化技術(shù)和優(yōu)化策略來實現(xiàn)在空間域內(nèi)的圖像重構(gòu)。
當(dāng)用超聲波從某一方向作用于物體的被成像截面時,超聲波與物體內(nèi)部介質(zhì)相互作用產(chǎn)生的全場滿足非齊次亥姆霍茲方程[3],其解可用Lippmann-Schwinger積分方程表示:
式中:p(→)和pin分別為全場和入射場;G(→-) 為自由空間的格林函數(shù);O(=k20[n2(→)-1]為物體內(nèi)部介質(zhì)聲學(xué)特性的未知函數(shù)。
在實際應(yīng)用中,我們只能測到物體的散射場和入射場,將式(1)轉(zhuǎn)化為散射場形式:常KM>N,即要求方程組為超定的。
通過矩量法將上述兩個卷積方程轉(zhuǎn)化為代數(shù)方程組的形式:
方程(3)、(4)僅描述了物體在某一方向入射波作用下接受到的散射波,為獲得足夠多的信息,需要從K個方向發(fā)射超聲波,由此得到K×M個方程,為使上述方程組可解,通
這個超定方程組的求解具有較強的非線性。這里采用Born迭代方法[4]。具體描述如下:
①先假設(shè)全場等于入射場P(t)=P(in),由方程(4)求得未知函數(shù)的初始解O0,常稱其為Born逆解;
②然后將Ok代入方程(3)求更接近實際的物體內(nèi)部的全場P(kt),k表示迭代次數(shù);
③由第②步求得的全場P(kt),代入方程(4),計算散射場P(ks)=DOP(kt),并求計算散射場與測量散射場的差ΔP(ks)=P(ks)-P(s),若‖ΔP(ks)‖小于事先給定的δ,則停止迭代,否則,轉(zhuǎn)④;
④根據(jù)散射場的改變量 ΔP(ks),由方程 ΔP(ks)=DΔOkP(kt求未知函數(shù)的改變量ΔOk,然后賦Ok+1=Ok+ΔOk作為未知函數(shù)的新值;
⑤由新的未知函數(shù)Ok+1(rn),返回步驟②,進一步求更近似的全場。
在迭代過程中,將遇到不適定性逆散射方程的求解問題,而對不適定性問題的正則化技術(shù)的選擇是否得當(dāng),將直接影響迭代算法的穩(wěn)定性。所以,正則化技術(shù)是圖像重建問題的關(guān)鍵。
算法1:用Tikhonov方法計算Ax=b的算法如下[6]。
①矩陣A的奇異值分解:[U,S,V]=svd(A),其中S=diag(σ1,σ2,…σn)
整數(shù)λ稱為正則化參數(shù)。
③作圖(lg‖xλ‖2),lg‖Axλ-b‖2)),確定曲線拐角點,拐點對應(yīng)的xλ是正則解xtik。
算法2:用TSVD方法計算Ax=b的算法如下[7]。
①矩陣 A 的奇異值分解:[U,S,V]=svd(A)。
②對1≤k≤T,計算解的范數(shù):
③作圖(lg(‖xλ‖2),lg(‖Axλ-b‖2)),確定曲線拐角點,拐點對應(yīng)的xk是正則解xtsvd。
研究表明,對離散不適定問題,lg(‖xk‖2),lg(‖Axk-b‖2)的關(guān)系曲線總是呈現(xiàn)L-曲線形狀,該關(guān)系曲線有一個明顯的拐角,最優(yōu)正則化參數(shù)位于L-曲線拐點附近[8]。因此,可通過尋找L-曲線拐點作為最優(yōu)正則化參數(shù)。
實驗裝置我們采用圓環(huán)形結(jié)構(gòu),發(fā)射器(同時也是接收器)等間隔地放在圍繞物體的圓環(huán)上,半徑取20×λ,背景介質(zhì)為水。采用F=200 kHz的超聲波作為入射波,超聲波在水中的波速約為C0=1 500 m/s,波長約為7.5 mm,對應(yīng)的波數(shù)為k0=0.837 758 mm-1。按照矩量法,我們對包含物體的正方形區(qū)域進行均勻采樣,采樣間隔均為d=λ/10=0.75 mm,采樣點的個數(shù)為35×35,劃分后單元的總數(shù)為N=1 225,且在每個小單元上做內(nèi)接圓,其半徑為a=d/2=0.375 mm。換能器共有50個,可工作于發(fā)射、接收兩種模式,每次僅有一個換能器發(fā)射,其他換能器接收來自物體的散射波。
如圖1所示,實驗選用的樣品如下:對比度為20%的人工繪制的圖形,輪廓信息較為豐富。對于已知像函數(shù)O(→r),先由全場方程(3)求出ROI中的全場分布,然后通過散射場方程(4)求出物體外的散射場P(s)分布,再使用BI進行圖像重建。兩種方法的重建結(jié)果如圖2和圖3所示:
圖1 目標原始函數(shù)
圖2 采用Tikhonov迭代4次結(jié)果
圖3 TSVD正則化方法迭代18次的結(jié)果
圖4 Tikhonov方法相對誤差曲線
由圖4可以看出采用Tikhonov方法時,相對誤差在第4次相對誤差最小,當(dāng)超過第4次后,相對誤差變大。圖5可以看出采用TSVD方法迭代到第18次時最小誤差恒定不變,因此在第18次中止迭代。對比兩種方法,TSVD方法明顯好于的Tikhonov方法重建的結(jié)果。
圖5 TSVD方法相對誤差曲線
本文運用Born迭代算法解決了方程的非線性問題,對于離散不適定問題,采用Tikhonov和TSVD兩種正則化方法分別對方程進行了求解。經(jīng)實驗驗證:TSVD正則化方法明顯優(yōu)于Tikhonov正則化方法。超聲逆散射成像技術(shù)是一個比較復(fù)雜的系統(tǒng)工程,算法的穩(wěn)定性及速度只是問題的一個方面。目前尚沒有成熟的理論和技術(shù)可用于產(chǎn)品的自動檢測,對該方面的一些關(guān)鍵技術(shù)進行研究具有一定的學(xué)術(shù)價值和廣闊的應(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].信號處理,2009,25(2):169-173.
[3]劉超,汪元美.超聲層析成像的理論與實現(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]韓旭,劉杰,李偉杰.時域內(nèi)多源動態(tài)載荷的一種計算反求技術(shù)[J].力學(xué)學(xué)報,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曲線法的電容層析成像正則化參數(shù)優(yōu)化[J].天津大學(xué)學(xué)報,2006,39(3):306-309.