基于AR模型的HIFU無(wú)損測(cè)溫方法研究

王亞，傅偉

（江蘇財(cái)經(jīng)職業(yè)技術(shù)學(xué)院 江蘇 淮安 223003）

超聲聚焦是一種先進(jìn)的無(wú)創(chuàng)傷手術(shù)，可在體外對(duì)生病部位做手術(shù)，不用開(kāi)刀，痛苦小；經(jīng)超聲照射后所形成的損傷區(qū)可準(zhǔn)確預(yù)測(cè)；組織凝固、變性的機(jī)制主要為超聲熱效應(yīng)，也伴有空化效應(yīng)，對(duì)人體無(wú)輻射；局部升溫的損傷區(qū)凝固后對(duì)熱吸收率加大，從而使周圍組織的溫度無(wú)變化，因而對(duì)焦斑以外的組織無(wú)損傷。

超聲聚焦需要進(jìn)行溫度監(jiān)控。病灶處溫度的高低以及溫度保持時(shí)間的長(zhǎng)短，都將直接影響照射的效果。長(zhǎng)期以來(lái)，缺乏一種無(wú)損二維溫度反饋被認(rèn)為是超聲聚焦應(yīng)用于深部熱療領(lǐng)域的主要障礙[1-2]。目前的治療反饋仍然靠熱電偶陣列獲得。雖然熱電偶是通用的，但它也是有損的，而且只能測(cè)量預(yù)先安置的傳感器所在位置（可能不符合實(shí)際感興趣組織的位置）的溫度，而且不適用于高強(qiáng)度的超聲外科手術(shù)。它們還有在手術(shù)中被損壞和充當(dāng)不想要的超聲散射體的潛在危險(xiǎn)。超聲無(wú)損測(cè)溫正好可以同時(shí)實(shí)現(xiàn)組織成像定位以及溫度監(jiān)控這兩個(gè)功能，充分滿足超聲聚焦的實(shí)際需要。因?yàn)樵\斷超聲與超聲聚焦同為超聲，所以可以非常方便的結(jié)合為一體。

高強(qiáng)度聚焦超聲治療技術(shù)是將高強(qiáng)超聲能量聚集于設(shè)定的焦點(diǎn)處，當(dāng)它作用于生物組織結(jié)構(gòu)的某處時(shí)，局部產(chǎn)生瞬間高溫，使癌細(xì)胞的蛋白質(zhì)達(dá)到凝固溫度（65℃以上）使其致死，以此達(dá)到治療的目的[3]。高強(qiáng)度聚焦超聲的焦點(diǎn)區(qū)的高強(qiáng)能量除破壞病變組織外，同樣也能破壞病變周圍的正常組織。為了殺死癌細(xì)胞而又不損傷正常組織，對(duì)組織溫度進(jìn)行測(cè)控非常重要，它直接影響到熱療的效果。

1 AR模型算法理論

1.1 生物組織離散隨機(jī)介質(zhì)模型

在超聲脈沖波掃描下，生物組織（如：肝臟、脾臟、腎臟等）具有如圖1所示的離散隨機(jī)介質(zhì)的散射模型。

在這個(gè)模型中有以下幾個(gè)假設(shè)：

圖1 超聲脈沖波作用下生物組織的離散隨機(jī)介質(zhì)模型Fig.1 Ultrasonic pulse wave discrete random medium model of biological tissue

1）生物組織和介質(zhì)，可以被認(rèn)為是一種半規(guī)則分布的散射體，這些散射體之間具有一定的間距d，稱之為平均散射元間距；

2）組織的平均散射間距d隨組織溫度的變化而變化，可以認(rèn)為d是組織內(nèi)溫度和組織的熱膨脹系數(shù)α的函數(shù)；

3）聲速c是溫度的函數(shù)，在水和大多數(shù)生物組織中，c隨溫度的升高而變大；在脂肪組織中，c隨溫度的升高而變小。

4）定義基頻f1=c/2d，它是由在介質(zhì)中以聲速c下兩個(gè)散射體間兩路行程（2d）來(lái)計(jì)算的。通過(guò)f1的變化就可以觀察到溫度T的變化（通過(guò)d和c的溫度相關(guān)性）。

1.2 接收信號(hào)模型

超聲信號(hào)在離散隨機(jī)介質(zhì)的生物組織中傳播，若為脈沖信號(hào)，假定組織介質(zhì)每個(gè)離散體的空間變化脈沖響應(yīng)為p（t），則超聲回波信號(hào)可以疊加表示為：

其中N為分辨單元內(nèi)散射元的個(gè)數(shù)，t為時(shí)間，Zi和Ai為傳感器成像軸第i個(gè)散射體的位置和幅度值，ci是介質(zhì)中隨空間變化的聲速，τi為傳感器的感應(yīng)面到第i個(gè)散射體回波信號(hào)通過(guò)并返回所用的時(shí)間（總行程為2Zi），如果介質(zhì)中隨空間變化的聲速能表示出來(lái)，則τi可以寫成：

隨空間變化的聲速為 c（ζ），方程（1）和（2）中包含了用 Pi（t）考慮介質(zhì)衰減和診斷超聲衍射的響應(yīng)。在頻域中方程（1）變?yōu)椋?/p>

其中 Pi（f）是 Pi（t）的傅立葉變換，它表示了第 i個(gè)散射中心的頻率相關(guān)的x（t）的頻譜密度函數(shù)（PSD）為：

其中E[*]是所期待的算子。

1.3 由PSD得到平均散射間距

平均散射體間隔對(duì)Sx的影響最簡(jiǎn)單情況是在實(shí)驗(yàn)區(qū)域內(nèi)一些散射中心以間隔d做規(guī)則分隔，此時(shí)方程（6）則變?yōu)椋?/p>

顯然當(dāng) fk=kc/2d，k=1，2，…∞ 時(shí)，方程（7）中的指數(shù)部分會(huì)變成1，一些頻譜線將在fk處疊加，也就是說(shuō)PSD會(huì)發(fā)生相關(guān)諧振現(xiàn)象。當(dāng)離散散射體之間不是均勻分布但有一定的規(guī)則性時(shí)，PSD中的局部峰值將由散射體中心的平均間距d決定。散射體間隔分布方差越大，則頻譜在f1，f2等處的諧振峰就越寬，越不確定。PSD中一些諧振峰值將會(huì)出現(xiàn)在f1和它的諧波 fk=k f1，k=1，2…∞ 處。

方程（6）和方程（7）顯示了 x（t）的 PSD 中包含了散射體分布的信息，因此 Sx（f）是d和fk的度量，可以間接的測(cè)量出溫度信息，典型的診斷傳感器用中心位于傳感器中心頻率處的高斯型P（f）來(lái)表示，只有出現(xiàn)在診斷傳感器的通帶內(nèi)f1的諧波在 Sx（f）中才能看到。

由上述的理論研究可以看到超聲無(wú)損測(cè)溫的本質(zhì)是利用回波信號(hào)的頻譜分析得出超聲體模中離散散射體的平均間隔d，然后根據(jù)d與溫度的相關(guān)性得到人體組織內(nèi)的溫度信息。因此可以模擬一個(gè)散射體的數(shù)學(xué)模型，假設(shè)散射體是平均分布的，如圖2所示。

圖2 離散介質(zhì)數(shù)學(xué)模型Fig.2 Mathematical model of discrete media

利用上面的數(shù)學(xué)模型可以得到超聲回波信號(hào)為:

其中 ak=fk（α，c），ak是組織的衰減系數(shù) α、聲速 c 的函數(shù)。

1.4 自回歸頻譜分析

在進(jìn)行從超聲回波信號(hào)中獲得體模溫度變化之前，必須首先計(jì)算它的頻譜，能夠產(chǎn)生高分辨率PSD-Sx（f）最簡(jiǎn)單也運(yùn)用最廣泛的是自回歸（AR）模型[4-6]，若要跟蹤和識(shí)別基頻f1和其諧振峰中的介質(zhì)溫度變化信息，需要一個(gè)高分辨率的PSD。AR-PSD是由AR參數(shù)所計(jì)算的，假定一個(gè)離散超聲回波信號(hào)為x（n），模擬經(jīng)過(guò)一個(gè)AR變換過(guò)程。一般來(lái)說(shuō)，一個(gè)AR過(guò)程是由下面的自相關(guān)序列所定義的：

其中 a（l），l=1，2…p 是 AR 模型系數(shù)參數(shù)，且定義 a（0）=1。p是模型的階數(shù)，σ2是AR模型的激發(fā)白噪聲方差。k=0時(shí)，δ（k）=1；k≠0 時(shí)，δ（k）=0。 為了計(jì)算 PSD，需要一種具有高分辨率和具有較強(qiáng)諧振細(xì)節(jié)的PSD算法，因?yàn)橛捎跍囟纫鸬膄1和其諧振的變化相當(dāng)微小。通常低于或相當(dāng)于標(biāo)準(zhǔn)傅立葉變換分析方法的細(xì)節(jié)限度。

方程（9）被稱為Yule-Walker方程，解方程中的AR參數(shù)有很多種方法，一種常用的方法就是Levinson-Durbin算法[7]，一旦AR參數(shù)算出之后，x（n）的AR-PSD可表示成：

其中 A（ f）=A（z）是在 z=exp（j2π f）下求得的，而 A（z）是AR參數(shù)的Z變換：

可見(jiàn)A（f）的零點(diǎn)將會(huì)在Sx（f）中形成波峰，功率譜密度中的某些峰值將相應(yīng)于用AR過(guò)程來(lái)模擬的由散射體間平均散射距離所形成的波峰。特別的是，在Sx（f）中這些諧振現(xiàn)象將會(huì)發(fā)生在f1和它所有的諧振峰中。

1.5 求解AR譜參數(shù)

要想得到較好的AR－PSD圖形，就必須準(zhǔn)確求得方程（9）中的 AR參數(shù) a[k]和 σ2。 在本文利用 Levinson_Durbin算法。設(shè)采集的超聲回波信號(hào)為離散時(shí)間信號(hào)x（n），則計(jì)算步驟如下:

1）根據(jù){X0，X1，…XN－1}計(jì)算 X（n）的自相關(guān)函數(shù)估計(jì) Rxx（k）。 （k=0，1，…Pmax），Pmax表示能達(dá)到的最大階數(shù)，可先估計(jì)。

3）利用下式遞推 k=2，3…Pmax時(shí)的各階 AR 系數(shù)

按這種遞推計(jì)算從P=2進(jìn)行下去，直到所計(jì)算得到的預(yù)測(cè)誤差的方差小于預(yù)先設(shè)定的門限值時(shí)或?yàn)V波器的階數(shù)達(dá)到預(yù)先確定的最大階數(shù)時(shí)。

整個(gè)AR譜估計(jì)法的流程如圖3所示。

圖3 AR譜估計(jì)流程Fig.3 AR spectral estimation process

2 仿真實(shí)驗(yàn)

為了驗(yàn)證超聲無(wú)損測(cè)溫系統(tǒng)仿真工作的正確性，進(jìn)行了肝組織樣品的實(shí)驗(yàn)。采用雙頻超聲脈沖反射模式，電路可以同步產(chǎn)生符合超聲測(cè)溫要求的雙超聲脈沖，并且能夠成功接收放大比較微弱的超聲回波信號(hào)。

2.1 實(shí)驗(yàn)過(guò)程及操作

建立一個(gè)實(shí)際的HIFU無(wú)損測(cè)溫實(shí)驗(yàn)系統(tǒng)，實(shí)驗(yàn)裝置的總體框圖如圖4所示。

圖4 實(shí)驗(yàn)裝置總體框圖Fig.4 Overall block diagram of the experimental apparatus

組織樣品采用的是新鮮的離體豬肝，在組織樣品被加熱（升溫）后，將超聲波回波信號(hào)作為時(shí)間的函數(shù)以20 MHz的采樣頻率采集1 000個(gè)數(shù)據(jù)點(diǎn)。整個(gè)實(shí)驗(yàn)過(guò)程中使用同一聚焦式超聲探頭，焦距設(shè)為67 mm，根據(jù)人體組織離散散射體模型理論，用肝組織模型中疊層厚度模擬肝組織中的平均散射體間隔d。豬肝的平均散射間距d近似0.11 μm，由平均均散射體間隔d與溫度T的相關(guān)性及d與基頻f的相關(guān)性，對(duì)于肝組織模型的實(shí)驗(yàn)，由回波信號(hào)計(jì)算回波信號(hào)的PSD諧振峰來(lái)代替溫度的計(jì)算，所用的加熱方法是直接加熱方法，在這種加熱情況下，肝組織樣品和超聲探頭二者間用充滿水的透明塑料薄膜隔開(kāi)，聚焦超聲探頭的焦點(diǎn)控制在豬肝樣品體表下的10 mm處，研究超聲波聲速在熱水及組織樣品中的溫度相關(guān)性以及超聲探頭本身的溫度相關(guān)性，發(fā)現(xiàn)在實(shí)驗(yàn)的小溫度變化（一般小于10℃）范圍內(nèi)，探頭的溫度相關(guān)性和超聲波聲速的溫度相關(guān)性可以忽略不計(jì)。并且在實(shí)驗(yàn)的溫度范圍內(nèi)設(shè)定超聲波在水中和組織中的聲速為1 500 m/s。

2.2 實(shí)驗(yàn)數(shù)據(jù)分析

將新鮮的離體豬肝經(jīng)過(guò)上述的實(shí)驗(yàn)，采樣時(shí)取得超聲回波信號(hào)時(shí)域圖如圖5所示。

圖5 超聲回波信號(hào)Fig.5 Ultrasonic echo signals

對(duì)采樣得到的20組超聲回波信號(hào)實(shí)行AR譜分析，取數(shù)據(jù)窗口長(zhǎng)度為2.0 mm，選取AR模型階次P=38，求得1 000點(diǎn)的回波信號(hào)PSD圖形，其中頻率分布為3～6 MHz，得到的一組數(shù)據(jù)的PSD圖形如圖6所示。

由圖6的AR－PSD圖形中可以觀察到，肝組織基頻f1的諧振峰fk在這個(gè)頻域內(nèi)形成的波峰還是非常明顯的，但有些頻譜圖像并不清楚，這時(shí)只能通過(guò)估計(jì)推算出在這段頻率內(nèi)出現(xiàn)的fk的位置，使用f1=fk+1－fk就可得到不同溫度下的f1。將頻譜圖中出現(xiàn)的波峰找出后求出f1的平均值，最后通過(guò)Δf和ΔT的溫度相關(guān)性求得了溫度變化。

圖6 由超聲回波信號(hào)得出的PSD圖形Fig.6 Obtained by the ultrasonic echo signal PSD graphics

圖7是對(duì)20組數(shù)據(jù)進(jìn)行AR-PSD分析后的頻譜圖。

圖7 由20組數(shù)據(jù)得到的PSD圖Fig.7 PSD FIG from 20 sets of data obtained

由圖7可以看出，每一組PSD圖像中的諧振波峰都有一點(diǎn)小的移動(dòng)，這正是所需要的頻譜移動(dòng)圖像，從中便可以很容易的提取出溫度信息。

3 結(jié) 論

本文描述了利用診斷脈沖超聲波無(wú)損檢測(cè)人體組織內(nèi)部溫度技術(shù)的算法理論研究及計(jì)算機(jī)仿真。理論研究已經(jīng)證明能從超聲回波信號(hào)的頻譜中提取出溫度信息。而且從PSD中提取的頻移和溫度變化的這種線性關(guān)系已由計(jì)算機(jī)仿真及實(shí)驗(yàn)得到定量的證明。

[1]任新穎，吳水才，曾毅.基于組織超聲回波時(shí)移的無(wú)創(chuàng)測(cè)溫技術(shù)及系統(tǒng)設(shè)計(jì)[J].北京生物醫(yī)學(xué)工程，2005，19（6）:42-46.REN Xin-ying，WU Shui-cai，ZENG Yi.Based on ultrasonic echo shifts of noninvasive temperature measurement technology and system design[J].Beijing Journal of Biomedical Engineering，2005，19（6）:42-46.

[2]夏雅琴，賈麗芹，彭見(jiàn)曙，等.一種用于超聲測(cè)溫的雙超聲脈沖發(fā)射接收電路的設(shè)計(jì) [J].中國(guó)醫(yī)療器械雜志，2002，26（4）:58-61.XIA Ya-qin，JIA Li-qin，PENG Jian-shu， et al.A method for ultrasonic measurement of double ultrasound pulse transmitting and receiving circuits design[J].Chinese Journal of Medical Instrumentation，2002，26（4）:58-61.

[3]ZU Wen-qia.Method for measuring the temperature in the body of human or animal with acoustic inversion[M].United States Patent Aplication Publication，2005.

[4]吳熙，錢盛友，孫福成，等.基于B超圖像處理的無(wú)損測(cè)溫方法研究[J].計(jì)算機(jī)工程與應(yīng)用，2007，43（32）:178-179.WU Xi，QIAN Sheng-you，SUN Fu-cheng，et al.Based on B ultrasound image processing study of noninvasive temperature measurement[J].Computer Engineering and Applications，2007，43（32）:178-179.

[5]宋平，錢盛友，馮艷玲，等.基于超聲散射回波分析的測(cè)溫方法研究[J].微計(jì)算機(jī)信息，2007，23（11-1）:267-269.SONG Pin，QIAN Sheng-you，F(xiàn)ENG Yan-ling，et al.Based on ultrasonic scattering echo analysis of temperature measurement method research[J].Micro Computer Information，2007，23（11-1）:267-269.

[6]劉正，候樹(shù)勛，任東風(fēng).高強(qiáng)度聚焦超聲治療過(guò)程中測(cè)溫技術(shù)的研究進(jìn)展[J].中國(guó)骨腫瘤骨病，2010，9（1）:84-87.LIU Zheng，HOU Shu-xun，REN Dong-fen.High intensity focused ultrasound in the treatment of process measurement technology research progress[J].Chinese Journal of Bone Tumor，2010，9（1）:84-87.

[7]Arthur RM，Straube WL，Trobaugh JW，Moros EG.Noninvasive estimation ofhyperthermia temperatures with ultrasound[J].International Journal of Hyperthermia，2005，21（6）:589-600.