陳國(guó)平，李怡萱

(重慶郵電大學(xué)電工理論與新技術(shù)實(shí)驗(yàn)室，重慶400065)

0 引言

微波致熱超聲掃描成像技術(shù)(microwave-induced thermo-acoustic tomography，MITAT)是生物醫(yī)學(xué)成像中頗有前途的一種技術(shù)，也是醫(yī)學(xué)診斷的一個(gè)重要工具，它兼有電磁參數(shù)高對(duì)比度和超聲波長(zhǎng)高分辨率的優(yōu)點(diǎn)。MITAT采用調(diào)制的微波脈沖照射生物組織［1］，微波脈沖提供高對(duì)比度的照射信號(hào)，由于腫瘤和背景介質(zhì)的吸收系數(shù)存在差異，因此，由熱效應(yīng)產(chǎn)生的熱膨脹使腫瘤組織產(chǎn)生高幅度的超聲波［2-4］，從內(nèi)向外傳播直至被組織周圍的傳感器陣列接收。由此可知腫瘤是作為MITAT的聲源，本文對(duì)腫瘤組織的定位其實(shí)是解決源反演的問(wèn)題。

目前大多數(shù)研究方法是基于模型的反演，通常使用均勻或隨機(jī)非均勻模型為背景。時(shí)間反轉(zhuǎn)鏡(time reversal mirror，TRM)方法被證明比后向背影(back projection，BP)方法具有更高的對(duì)比度［5］，與BP相比，TRM技術(shù)可以有效地克服多徑效應(yīng)，消除復(fù)雜環(huán)境對(duì)信號(hào)造成的時(shí)延，它對(duì)非均勻媒質(zhì)中的目標(biāo)具有良好的成像定位能力。信號(hào)在非均勻媒質(zhì)中傳播時(shí)，時(shí)間反轉(zhuǎn)鏡可以實(shí)現(xiàn)有效能量的相關(guān)疊加，提高信噪比，最終達(dá)到最佳的空時(shí)匹配濾波。若信號(hào)傳播的環(huán)境不變，接收器陣列發(fā)射的時(shí)反信號(hào)在源發(fā)射原始信號(hào)的位置重新聚焦，并且介質(zhì)環(huán)境越復(fù)雜，重聚焦的信號(hào)質(zhì)量越高，因此大大提高了對(duì)目標(biāo)物體的定位精度［6-8］。然而TRM的實(shí)質(zhì)是通過(guò)接收信號(hào)特性和物理機(jī)制以及環(huán)境模型來(lái)實(shí)現(xiàn)對(duì)目標(biāo)的定位要求，但是通常生物組織中密度和聲速等參數(shù)的分布情況是未知的，MITAT中熱聲源(腫瘤組織)信息也是未知的，因此，生物組織環(huán)境的建模成為腫瘤定位一個(gè)重要的前置過(guò)程。

較為先進(jìn)的快速行進(jìn)方法(fast marching method，F(xiàn)MM)首先假設(shè)一個(gè)環(huán)境分布，根據(jù)測(cè)量數(shù)據(jù)得到的渡越時(shí)間(time of fly，TOF)和路徑計(jì)算得到的TOF運(yùn)用同時(shí)代數(shù)重建技術(shù)(simultaneous algebraic reconstruction technique，SART)對(duì)假設(shè)的環(huán)境分布進(jìn)行迭代更新，從而重建真實(shí)的環(huán)境分布［9］。在一定條件下TR-FMM能較好地估測(cè)出信道參數(shù)的分布，也能解決TRM技術(shù)的難點(diǎn)，但是在非均勻的生物組織中使用TR-FMM算法只考慮了幾何光學(xué)和費(fèi)馬原理，沒(méi)有考慮衍射方面，并且TR-FMM必須知道天線的起始位置，有時(shí)也不能滿足復(fù)雜生物環(huán)境中腫瘤這一類主動(dòng)源的高精度定位需求。

為了進(jìn)一步提高M(jìn)ITAT對(duì)生物組織中腫瘤定位的精確度，本文提出TR-adjoint方法，將時(shí)間反轉(zhuǎn)技術(shù)和adjoint理論結(jié)合應(yīng)用在生物組織腫瘤定位的研究中。該方法使用腫瘤作為主動(dòng)源發(fā)射信號(hào)，天線獲得傳播結(jié)果，其中包括生物組織的非均勻介質(zhì)分布信息，然后用adjoint方法提取相關(guān)信息，建立并迭代更新合成的聲速模型，使合成腫瘤信號(hào)與觀測(cè)腫瘤信號(hào)之間的差異最小化，更新后的模型可用于TRM定位腫瘤源。這是一種僅通過(guò)觀測(cè)接收到的信號(hào)，就能自適應(yīng)求解復(fù)雜傳播環(huán)境信道信息的新算法，具有較強(qiáng)的研究?jī)r(jià)值。

1 TR-adjoint原理及應(yīng)用

MITAT中影響目標(biāo)成像的因素有聲密度、聲速度以及發(fā)射源的位置，通常聲密度的影響非常小，可以被忽略，本文僅討論聲速度的影響。TR-adjoint方法的主要思想是假設(shè)一個(gè)介質(zhì)參數(shù)模型，根據(jù)此模型在接收器處獲得合成信號(hào)，因此，我們可以根據(jù)真實(shí)的和假設(shè)的介質(zhì)參數(shù)模型得到觀測(cè)信號(hào)和合成信號(hào)之間的差異，這種差異可以通過(guò)走時(shí)的失配函數(shù)來(lái)量化，然后根據(jù)這個(gè)差異用共軛梯度法迭代更新假設(shè)的模型參數(shù)，直至使觀測(cè)信號(hào)和合成信號(hào)之間的差異在一定準(zhǔn)則下達(dá)到最小，從而獲得理想的參數(shù)模型，即足夠接近真實(shí)情況，然后對(duì)腫瘤源進(jìn)行準(zhǔn)確定位。

1.1 走時(shí)失配函數(shù)分析

我們的目的是最小化二乘走時(shí)失配函數(shù)［10］，即最小化合成數(shù)據(jù)和真實(shí)數(shù)據(jù)之間的差異，關(guān)鍵是計(jì)算出衡量差異大小的失配函數(shù)的靈敏度內(nèi)核。定義走時(shí)失配函數(shù)為

(1)式中:m表示模型參數(shù)，非均勻生物組織的介質(zhì)參數(shù)可以用m(ρ，c)表示，ρ是聲密度，c是聲速分布;N為接收天線總數(shù);Tr(m)為在假設(shè)的模型m下第r根接收天線預(yù)測(cè)合成信號(hào)的走時(shí);Torbs為觀測(cè)信號(hào)的走時(shí)。失配函數(shù)的梯度為

(2)式中，δTr是走時(shí)Tr的Fréchet導(dǎo)數(shù)。這里Fréchet導(dǎo)數(shù)定義為位函數(shù)對(duì)空間任意一點(diǎn)介質(zhì)某個(gè)物理量變化的靈敏度函數(shù)，根據(jù)合成的和觀測(cè)的波形互相關(guān)有如下表達(dá)［11-12］，

(3)式中:wr是時(shí)間窗函數(shù)，取決于模型參數(shù)的變化范圍，以確保穩(wěn)定的互相關(guān)，比如最高和最低的聲波速度可以確定最短和最長(zhǎng)的時(shí)窗;T表示一個(gè)周期，si(xr，t)表示第i個(gè)接收天線獲得的合成信號(hào);δsi是由模型擾動(dòng)δm引起的合成波場(chǎng)的變化。Nr是歸一化因子，定義為

(4)式并中，δsi是由模型擾動(dòng)δm引起的合成波場(chǎng)的變化。將δsi的玻恩近似形式［13］代入到(3)式和

(2)式中得:

(5)式中，sj(x'，t')表示si(xr，t)的時(shí)間反轉(zhuǎn)傳播;si(xr，t)是可以通過(guò)接收器觀測(cè)到的信號(hào)。接收器用時(shí)反信號(hào)作為虛擬發(fā)射源，計(jì)算觀測(cè)數(shù)據(jù)與合成數(shù)據(jù)時(shí)反的差異，即產(chǎn)生adjoint波場(chǎng)。利用格林張量的相互作用Gij(x，x';t－t')=Gji(x'，x;t－t')，定義走時(shí)adjoint波場(chǎng)為

因此，走時(shí)的adjoint源定義為

根據(jù)下面各向同性的介質(zhì)模型參數(shù)之間的關(guān)系［14］:

(8)式中:κ是體積模量;μ是剪切模量，因?yàn)镸ITAT源反演方法通常是在聲學(xué)近似下進(jìn)行的，往往假設(shè)生物組織的傳播介質(zhì)是各向同性且靜止的，產(chǎn)生的壓力波是無(wú)旋流的，所以，可以忽略剪切模量μ(x)。(8)式中參數(shù) δjk，δlm，δjl，δkm，δjm，δkl具體說(shuō)明見(jiàn)文獻(xiàn)［14］。另外，密度參數(shù)被認(rèn)為是均勻分布的。因此，將(8)式代入到(5)式中，則失配函數(shù)的擾動(dòng)可以簡(jiǎn)化用體積模量來(lái)表示為

根據(jù)TR-adjoint的目的是估計(jì)生物組織環(huán)境的聲速分布，則聲速內(nèi)核Kc用體積模量?jī)?nèi)核Kκ表示為

在對(duì)電力系統(tǒng)變電站運(yùn)行管理的過(guò)程中，應(yīng)該結(jié)合國(guó)家的相關(guān)要求，與電氣企業(yè)的發(fā)展現(xiàn)狀，做好科學(xué)合理的規(guī)劃工作，從而保證電力變電系統(tǒng)的順利運(yùn)行。在這個(gè)基礎(chǔ)上，應(yīng)該落實(shí)各工作人員的工作職責(zé)，要求每一個(gè)工作人員都需要遵守相關(guān)的制度，按時(shí)完成工作任務(wù)，做好變電運(yùn)行的保證工作。

(11)式中，s(x，T－t)·s(x，t)實(shí)際上是合成波場(chǎng)和伴隨波場(chǎng)中每個(gè)時(shí)間幀互相關(guān)的總和，以上理論分析出了聲速c(x)與走時(shí)失配函數(shù)的Fréchet導(dǎo)數(shù)的關(guān)系，可知由當(dāng)前模型的波場(chǎng)和adjoint波場(chǎng)之間的相互作用得到參數(shù)靈敏度內(nèi)核，隨后，再用共軛梯度法不斷迭代更新假設(shè)的模型，直到找到最佳的參數(shù)模型為止。

1.2 共軛梯度法

本文應(yīng)用非線性共軛梯度法來(lái)迭代參數(shù)模型，該方法開(kāi)始于一個(gè)初始的均勻模型m0，上文中提到非均勻生物組織的介質(zhì)參數(shù)可以用m(ρ，c)表示，m0即表示介質(zhì)參數(shù)密度和聲速是均勻分布的初始模型，在該初始模型下進(jìn)行前向模擬產(chǎn)生合成信號(hào)。由(1)式計(jì)算出初始的失配函數(shù)χ0，由(7)式計(jì)算出走時(shí)的adjoint源，adjoint源在與前向模擬相同的模型中進(jìn)行時(shí)間反轉(zhuǎn)傳播，由(11)式計(jì)算出失配函數(shù)梯度的靈敏度內(nèi)核Kκ，對(duì)不同位置的接收天線得到的初始Kκ進(jìn)行求和，將此和定義為g0，g0=表示失配函數(shù)的初始梯度。并設(shè)最初的共軛梯度搜索方向等于負(fù)的失配函數(shù)初始梯度，即p0=－g0，如果‖p0‖＜ε，其中ε是一個(gè)極小的值，那么模型m0則是要求的模型。反之，

①進(jìn)行線性搜索獲得標(biāo)量vk，使的值最小化，其中:

②更新參數(shù)模型:mk+1=mk+vkpk，然后計(jì)算

③更新共軛梯度的搜索方向p:

2 TR-adjoint的仿真實(shí)例

為了驗(yàn)證我們所述的TR-adjoint方法對(duì)生物組織腫瘤定位的有效性，我們進(jìn)行了如下的仿真實(shí)驗(yàn)。在仿真過(guò)程中，我們主要考慮生物組織中聲速模型估計(jì)的準(zhǔn)確性及腫瘤定位的精確性。

圖1為假設(shè)的均勻介質(zhì)參數(shù)模型、主動(dòng)源和接收器陣列的模擬設(shè)置圖。仿真區(qū)域采用20×20像素的完全匹配層(perfectly matched layer，PML)吸收邊界，總像素?cái)?shù)為250×250，模擬的生物組織大小為150 mm×150 mm。如圖1所示，傳播介質(zhì)密度為1 000 kg/m3，傳播介質(zhì)聲速為1 500 m/s。200個(gè)接收器均勻地分布在半徑為60 mm的陣環(huán)上，我們假設(shè)主動(dòng)源(腫瘤組織)位于離中心位置在X軸正方向20個(gè)像素，發(fā)出的信號(hào)強(qiáng)度歸一化并且波形是高斯脈沖。

圖1 假設(shè)參數(shù)模型、主動(dòng)源和接收器陣列設(shè)置Fig.1 Supposed model parameters，active source and receiver array

根據(jù)初始假設(shè)的均勻介質(zhì)參數(shù)模型開(kāi)始應(yīng)用本文所提出的TR-adjoint方法，每迭代一次，模型就被更新一次，根據(jù)更新的模型經(jīng)過(guò)前向模擬后即更新了合成信號(hào)，然后由合成信號(hào)重新計(jì)算(9)式的聲速內(nèi)核和(7)式的伴隨源，直到觀測(cè)信號(hào)和更新后的合成信號(hào)之間的走時(shí)失配函數(shù)達(dá)到收斂。真實(shí)的介質(zhì)參數(shù)模型、主動(dòng)源和接收器陣列設(shè)置如圖2所示。假設(shè)的參數(shù)模型經(jīng)過(guò)12次迭代后的仿真結(jié)果如圖3所示。

圖2為真實(shí)的介質(zhì)參數(shù)模型、主動(dòng)源和接收器陣列的模擬設(shè)置圖。仿真區(qū)域采用20×20像素的完全匹配層吸收邊界，總像素?cái)?shù)為250×250，模擬的生物組織的大小為150 mm×150 mm。如圖2所示，傳播介質(zhì)密度為1 000 kg/m3，傳播介質(zhì)聲速為1 500 m/s，為驗(yàn)證重建的生物組織的聲速度分布的有效性，設(shè)置2個(gè)半徑為15 mm，聲速為1 650 m/s的高速區(qū)分布在均勻介質(zhì)區(qū)域(167，167)，(83，83)處，另外2個(gè)半徑為15 mm，聲速為1 350 m/s的低速區(qū)位于仿真區(qū)域(83，167)，(167，83)處，200個(gè)接收器均勻的分布在半徑為60 mm的陣環(huán)上。主動(dòng)源腫瘤組織位于仿真區(qū)域的中心位置(125，125)，發(fā)出的信號(hào)強(qiáng)度歸一化并且波形是高斯脈沖。

圖2 真實(shí)參數(shù)模型，主動(dòng)源和接收器陣列設(shè)置Fig.2 Real model parameters，active source and receiver array

圖3 假設(shè)的參數(shù)模型經(jīng)過(guò)12次迭代后的仿真結(jié)果Fig.3 Results of supposed model with 12 iterations

圖3是經(jīng)TR-adjoint方法得到的生物組織聲速模型，從圖3中可以明顯看出經(jīng)過(guò)12次迭代，假設(shè)的參數(shù)模型已經(jīng)極其接近圖2中的真實(shí)模型。

圖4顯示了每次迭代更新后失配函數(shù)的大小，失配函數(shù)用來(lái)衡量真實(shí)與合成信號(hào)之間差異的大小，因此失配函數(shù)的值越小，表明差異越小，當(dāng)失配函數(shù)達(dá)到收斂時(shí)即差異達(dá)到最小，此時(shí)假設(shè)的模型最接近真實(shí)情況，從圖4中可以看到，在第4次或第6次迭代后失配函數(shù)的大小已經(jīng)下降了92%，12次迭代后失配函數(shù)達(dá)到收斂，大小幾乎接近于零，可見(jiàn)此時(shí)的模型即為所要求的最佳聲速模型。最后在此模型下運(yùn)用TRM技術(shù)對(duì)腫瘤組織進(jìn)行定位，結(jié)果如圖5所示。圖2中真實(shí)發(fā)射源處于仿真區(qū)域(125，125)，仿真結(jié)果中腫瘤組織位于(124.745 3，124.745 4)，由此可知定位誤差為0.255 mm。比在生物組織常規(guī)環(huán)境的TRM定位技術(shù)提高將近60多倍。

圖4 迭代次數(shù)與失配函數(shù)大小關(guān)系Fig.4 Relationship between the number of iterations and the size of the misfit function

圖5 腫瘤組織定位結(jié)果Fig.5 Tumor cell localization results

在實(shí)現(xiàn)TR-Adjoint算法時(shí)，為了避免在仿真運(yùn)行中存儲(chǔ)所有s場(chǎng)與s+場(chǎng)的時(shí)間步長(zhǎng)，可以通過(guò)在特定的一個(gè)步長(zhǎng)中使用TR方法得到s(x，t)場(chǎng)與s(x，T－t)場(chǎng)的相互作用。在一個(gè)時(shí)不變的媒質(zhì)中檢測(cè)收集信號(hào)時(shí)，TR方法可以在T－t時(shí)間步長(zhǎng)時(shí)重現(xiàn)波場(chǎng)，這些過(guò)程可通過(guò)把時(shí)反信號(hào)作為初始條件來(lái)實(shí)施。另外，算法需要使用正反向模擬波場(chǎng)，計(jì)算量較大，為提高仿真效率，算法使用基于快速傅立葉變換(fast Fourier transform，F(xiàn)FT)的偽譜時(shí)域算法(pseudo spectral time domain，PSTD)計(jì)算內(nèi)核以實(shí)現(xiàn)波場(chǎng)的精確計(jì)算。共軛梯度迭代算法實(shí)現(xiàn)簡(jiǎn)單，內(nèi)存需求小，在全波形反演中得到廣泛應(yīng)用。

對(duì)于定位一個(gè)主動(dòng)發(fā)射源的目標(biāo)物體來(lái)講，在64位的Windows 7系統(tǒng)平臺(tái)上，利用Matlab 7.1軟件，其中計(jì)算機(jī)的配置為4 GB的內(nèi)存，雙CPU，運(yùn)用adjoint算法對(duì)乳腺組織中的腫瘤定位技術(shù)僅僅需要10 min的時(shí)間。如果在實(shí)驗(yàn)以及仿真中有更加精良的硬件設(shè)備，那么定位精度將會(huì)更加準(zhǔn)確。所有的這些都表明adjoint對(duì)生物組織腫瘤這一類主動(dòng)源定位具有很大的價(jià)值。

3 結(jié)束語(yǔ)

在介質(zhì)參數(shù)模型未知的MITAT環(huán)境中應(yīng)用TR-adjoint方法，可提高生物組織中腫瘤定位的精度。仿真結(jié)果表明，TR-adjoint方法不僅可以較逼真地重建介質(zhì)參數(shù)模型，而且可以對(duì)腫瘤進(jìn)行高精度定位，定位誤差僅為0.255 mm，比傳統(tǒng)TRM方法有較明顯的定位優(yōu)勢(shì)。此外，該方法計(jì)算波場(chǎng)時(shí)采用全波模擬技術(shù)，可以滿足MITAT對(duì)成像的有限頻率特性要求。TR-adjoint方法顯著的優(yōu)點(diǎn)是僅通過(guò)接收器接收到的信號(hào)就可以估計(jì)和更新所有的模型參數(shù)，進(jìn)而對(duì)腫瘤源進(jìn)行高精度定位，因此在實(shí)際應(yīng)用中具有相當(dāng)大的價(jià)值及前景。

［1］KELLNBERGER S，HAJIABOLI AR，RAZANSKY D，et al.Near-field thermoacoustic tomography of small animals［J］.Phys.Med Biol，2011，56(11):3433-3444.

［2］LIU L，XU Y，WANG L V.Transcranial Thermo-acoustic Tomography:A Comparison of Two Imaging Algorithms［J］.Medical Imaging，IEEE Transactions on，2013，32(2):289-294.

［3］CHEN GuoPing，ZHAO Zhiqin，GONG Wei，et al.The Development of the Microwave-Induced Thermo-acoustic Tomography Prototype System［J］.Chinese Science Bulletin，2009，52(12):1786-1789.

［4］XIONG W，BAUER D R，WITTE R，et al.Microwave-Induced Thermo-acoustic Imaging Model for Potential Breast Cancer Detection［J］.IEEE Trans Biomed Eng，2012，59(10):2782-2791.

［5］XU M，WANG L V.Pulsed-microwave-induced thermoacoustic tomography:Filtered backprojection in a circular measurement configuration［J］.Medical Physics，2002，29(8):1661.

［6］BAL G，RYZHIK L.Time reversal and refocusing in random environment［J］.SIAM Appl Math，2003(63):1375-1498.

［7］MOSZKOWSKI Ben.Compositional Reasoning using Intervals and Time Reversal［J］.Annals of Mathematics and Artificial Intelligence，2014，7(1):175-250.

［8］XIA Yunlong，F(xiàn)U Yongqing.Refocusing with Time Reversal Mirror in Random Medium［J］.Procedia Engineering，2012，(29):2600-2604.

［9］WANG Jingguo，ZHAO Zhiqin，SONG J，et al.An Iterative Method for Correcting the Acoustic Refraction in Heterogeneous Tissue for Microwave-Induced Thermo-Acoustic Tomography［J］.Progress in electromagnetics research，2012(130):225-240.

［10］周小娟，李春曉.基于偏最小二乘分析和稀疏表示的目標(biāo)跟蹤算法［J］.重慶郵電大學(xué)學(xué)報(bào):自然科學(xué)版，2014，(1):104-110.

ZHOU Xiaojuan，LI Chunxiao.A partial least squares analysis and sparse representation of the target tracking algorithm based on［J］.Journal of Chongqing University of Posts and Telecommunications:Natural Science Edition，2014，(1):104-110.

［11］MARQUERING H，DAHLEN F，NOLET G.Three-dimensional sensitivity kernels for finite-frequency traveltimes:the banana-doughnut paradox［J］.Geophysical Journal International，1999，137(3):805-815.

［12］DAHLEN F，HUNG S H，NOLET G.Fréchet kernels for finite frequency traveltimes—I.Theory［J］.Geophysical Journal International，2000，141(1):157-174.

［13］HUDSON J A，HERITAGE J R.The use of the Born approximation in seismic scattering problems［J］.Geophysical Journal of the Royal Astronomical Society，1981，66(1):221-240.

［14］FICHTNER A，BUNGE H P，IGEL H.The adjoint method in seismology:I.Theory［J］.Physics of the Earth and Planetary Interiors，2006，157(1):86-104.