用于活體溫度評估的非線性超聲熱應變模型?

孫彪 郭霞生 屠娟,2 章東,2

(1 南京大學物理學院 南京 210093)

(2 中國科學院聲學研究所 聲場聲信息國家重點實驗室 北京 100190)

0 引言

在腫瘤和高血壓等疾病的治療研究中，熱療被認為是具有重要價值的微創或非侵入性技術[1?2]。考慮到熱療的安全性和有效性，在治療過程中對生物組織的溫度變化進行監控很有必要。其中，超聲輻照腎周脂肪治療高血壓過程中的溫度監控是本文作者團隊近年來關注的重點。

與其他溫度成像方法相比，超聲測溫技術具有無電離輻射、靈活、易操作、低成本、便攜和兼容性好等優勢，被認為是一種有前景的溫度評估方法。其原理主要依賴于溫度變化引起超聲傳播過程的變化，具體包括散射點信號強度變化[3?4]、回波信號的譜峰偏移[5]等。其中，基于超聲回波偏移的熱應變測溫技術受到廣泛關注[6?8]。

超聲熱應變測溫基于組織熱膨脹和聲速隨溫度線性變化的假設。當組織溫度發生變化時，由于熱膨脹和聲速變化引發回波時移，該時移的軸向梯度(即熱應變)與溫度變化量存在線性關系，可用于溫度估計[9]。Simon等[10]對加熱的凝膠進行實驗，在4.22℃溫度變化范圍內實現的評估誤差小于0.5℃。Liu等[11]對豬心臟組織加熱，借助GPU運算實現了溫度的實時監控。但當組織溫度高于50℃(高于體溫13℃)時，組織的性質通常發生較大變化，該線性溫度評估模型不再適用。即對于生物組織，線性的聲速變化和熱膨脹僅在有限的溫度范圍適用[9]。因此，對于更大范圍的溫度評估，線性熱應變測溫模型需要改進。

對活體的溫度變化進行評估時，熱應變遠小于組織運動導致的機械應變，容易被后者掩蓋[9]。因此，有效抑制活體生理運動同樣是臨床應用之前要解決的關鍵問題。Daniels等[12]考慮活體呼吸的周期性，采取動態幀選取的方法來抑制呼吸運動，進而降低豬腎臟中溫度評估的誤差。Bayat等[13]用自適應運動補償的方法來消除呼吸脈搏對活體超聲評估的干擾。對于上述溫度評估方法，報道的溫度評估誤差為4℃左右，并不能滿足實際的要求。在之前的工作中，動態幀選取和自適應濾波結合的常系數溫度估計(Constant coefficient model with dynamic frame selection and adaptive filtering，簡稱CDA)算法可以很好地抑制活體中生理運動對溫度評估的干擾[14]。

本文進一步發展超聲熱應變測溫理論，提出一種非線性熱應變測溫模型，結合之前發展的運動抑制算法對活體溫度變化進行評估。以豬的腎周脂肪作為目標區域，進行溫度評估實驗。在3種不同溫度變化范圍的實驗中，以最大相關系數動態幀選取和自適應濾波相結合的方法對活體豬的呼吸心跳的干擾進行抑制，并取得了很好的效果。本算法將推動超聲測溫在臨床熱療中的應用，可顯著提高熱療的安全性和有效性。

1 溫度評估模型和算法

1.1 非線性熱應變模型

超聲成像中，軸向深度z處的散射子引起的時間延遲可以描述為聲速和溫度的函數：

其中，T0(ξ)是軸向位置ξ處的溫度，c[ξ,T0(ξ)]為軸向位置ξ處溫度為T0(ξ)時的聲速。

加熱時，目標區域組織發生熱膨脹。若溫度變化范圍較大，熱膨脹引起的體積變化與溫度之間的關系是非線性的[11]，聲路徑上兩個散射點之間的軸向距離l隨溫度變化的關系近似描述為

進一步，認為聲速也是溫度的非線性函數[11]，表示為

當組織溫度從T0變化為T時，z處回波的時移(回波時延的變化量)為

其中，tf(z)和ti(z)分別為T和T0時刻的回波時延。將方程(1)和方程(4)代入方程(5)，得到

將式(6)對z微分，并忽略高階小量，得到

進一步考慮到∣β2(z)δ2T(z)∣?1+β1(z)δT(z)，得到

假設加熱區的組織為單一組分，即區域內各處具有相同的基線溫度，同時考慮到[β1(z)δT(z)]2?1，由式(8)得到

其中，k1=α2(z)?β2(z)+β1(z)[α1(z)?β1(z)]和k2=α1(z)?β1(z)為非線性溫度評估模型的兩個測溫系數。

通過求解方程(9)，得到時移與溫度變化之間的關系：

此即為文中采用的非線性溫度評估模型。如忽略其中的二階小量，則方程(10)化簡為

與文獻[10,15]中描述的線性熱應變測溫模型一致。由此可見，對于小范圍的溫度變化的情況，式(10)中的二階變化量很小。因此，非線性模型也適用于小范圍的溫度變化情況。

1.2 CDA算法

與線性模型的溫度評估算法類似[10,15]，非線性溫度評估模型也是基于組織回波偏移估計獲得具體的溫度估計。而對于活體組織的回波偏移估計，應考慮活體生理運動的干擾。對本文考慮的靶區組織腎周脂肪而言，其生理運動主要由呼吸和心跳運動共同引起[9]。本文采用CDA算法來對二維溫度圖像中的噪聲信號進行抑制，該算法包含動態幀選取、熱應變計算和自適應濾波3部分。動態幀選取算法通過計算超聲射頻(Radio frequency,RF)圖像間的歸一化互相關系數，并據此選擇運動狀態相似的圖像，以抑制短周期噪聲(心跳)引起的干擾。自適應濾波算法則基于呼吸運動的全局性，用自適應濾波器(Adaptive filter,AF)對長周期的呼吸干擾進行運動補償。CDA算法流程圖如圖1所示，具體方法可以參考文獻[14]。該算法大致按如下3個步驟進行：

圖1 算法流程Fig.1 The f low diagram of the CDA algorithm

(1)對升溫過程的超聲RF圖像，以心跳周期為選幀范圍，使用動態幀選取算法挑選出最大相關性的RF圖像序列I1。在升溫前的RF圖像中，對I1逐幀挑選相關性最好的圖像作為非升溫的參照圖像序列I2。

(2)基于選取的RF圖像I1、I2，通過散斑追蹤計算RF圖像中的散射點位移分布，再依次應用軸向數字差分濾波器和軸向、橫向的低通濾波得到熱應變分布S1及S2。

(3)設計歸一化最小均方(Normalized least mean square,NLMS)自適應濾波器，通過S2所包含的機械應變信息，以步長0.1迭代以確定濾波器的系數。然后將濾波器應用于含噪聲的熱應變圖像S1，得到降噪后的熱應變分布，若已知溫度系數k1和k2，則可得到溫度估計結果。

1.3 實驗設置

實驗前將豬禁食12 h，完成清洗后，用異丙酚和七氟醚對豬進行麻醉。使用呼吸機對處于深度睡眠狀態下豬進行輔助呼吸。實驗裝置示意圖如圖2所示。對豬備皮，使用B超成像探頭對腎周脂肪定位，在平行于脊柱的位置使用高頻電刀開出長約5 cm的創口。使用128陣元的L12-5線陣探頭和便攜式超聲設備(定制，珠海醫凱超聲設備有限公司)對目標區域成像并采集RF數據。在B超成像畫面的引導下，將通過模具固定的微波消融針(ECO-200G，南京億高微波系統工程有限公司)和熱電偶(HYP2,OMEGA)沿成像平面的法向插入腎周脂肪。微波消融針的發射單元和熱電偶的感溫區域均在尖端，并位于成像平面上。微波消融針和熱電偶的間隔為1.3 mm。動物實驗準備、麻醉和手術得到了江蘇省人民醫院的支持，且遵循美國國家衛生研究院《實驗動物的護理和使用指南》(2011)，并經南京總醫院倫理委員會批準(批準號20140510)。

圖2 實驗裝置Fig.2 The experimental setup

實驗中，在加熱前采集至少兩個呼吸周期長度的非升溫RF信號。將微波消融針的發射功率設定在不同水平，在每個加熱過程，同步獲得熱電偶測量的溫度和B超圖像的RF數據。熱電偶采樣率20 Hz，RF數據采樣率40 MHz，幀率32 Hz。每幀RF圖像尺寸為2608(軸向)×128(橫向)像素。

2 結果

2.1 二維溫度分布

根據多次微波加熱活體豬腎周脂肪的實驗，采集37℃～50℃升溫過程的RF數據用以計算熱應變，根據熱電偶在RF圖像中的位置提取熱電偶處的熱應變估計值，將其和熱電偶測得的溫度曲線進行最小二乘擬合確定非線性模型溫度系數k1、k2。由此確定的非線性系數(二次項系數)為k1=0.0045，線性溫度系數(一次項系數)k2=0.0930。一組超聲CDA算法評估的結果如圖3所示，其中給出了靶區在4個不同時刻的二維溫度分布圖像。由圖3可知，隨著時間的增加，組織溫度逐漸升高，溫升區域范圍不斷擴大。對于真實生物組織傳熱過程，無論是熱傳導方程計算的模擬結果，還是紅外成像記錄的真實溫度擴散過程，加熱點的溫度變化都是由低到高，溫升范圍都是由中間向四周擴散。因此本文的超聲溫度評估算法的溫度和真實的物理過程具有很好的一致性。

圖3 二維溫度變化序列Fig.3 Two-dimensional temperature change sequence diagram

2.2 兩種不同模型的對比

圖4分別給出了對應3種不同的加熱功率(對應不同的溫升范圍)下，傳統線性模型和本文非線性模型的溫度評估結果，其中線性模型采用了式(11)描述的時移-溫度關系。對于小范圍的溫度變化(12℃，滿足線性模型溫升的范圍)，線性模型和非線性模型的最大溫度偏差都在2℃左右，在這種情況下，兩者的溫度評估精度是相當的；當溫度變化范圍提高到17℃時，線性模型溫度評估的最大評估誤差為3.04℃，而非線性模型的最大誤差為2.03℃。相較而言，非線性模型具有更高的溫度評估精度。隨著微波功率繼續增加，組織最大溫升達到25℃，線性模型的溫度評估最大誤差達到5.08℃，而非線性模型的溫度評估誤差最大為2.34℃，兩種模型的性能差距進一步擴大。因此，對于較大溫升范圍過程中的溫度評估，非線性模型更具有優勢。

圖4 線性模型與非線性模型升溫曲線對比Fig.4 Temperature curves of linear model and nonlinear model

3 討論

對于超聲評估活體溫度的算法，溫度評估模型是影響溫度評估準確性的重要因素。研究表明，無論是活體還是離體，只有在溫度變化較小的情況，線性熱應變模型的評估結果才會比較準確[8,16]。當生物體溫度高于50℃，基于回波偏移的線性熱應變模型無法對生物組織的升溫過程進行有效的評估[10]。對于大范圍的溫度變化，生物組織的體積和聲速隨溫度的變化存在非線性關系[17]。因此，本文提出將熱膨脹和聲速描述為溫度的非線性函數，為較大范圍溫度變化下的溫度評估問題提供了合適的解決方案。

需要指出的是，本文的模型僅考慮了軸向的回波時移，而實際組織受熱時溫度場是全向擴散的。事實上，當組織發生溫度變化時，由于熱膨脹和聲速變化的共同作用產生軸向的回波偏移。隨著加熱過程的進行，RF圖像在橫向上僅存在由于熱膨脹導致的圖像變形，在軸向上存在由于熱膨脹和聲速變化導致的圖像變化。但重要的是，熱膨脹形成的回波偏移遠小于聲速變化導致的偏移[9]，因此文中所得熱應變主要源于聲速變化的貢獻。盡管如此，如能同時從橫向考慮組織受熱導致的圖像變形，是完全可能進一步提高精度的。但這通常很困難，因為B超圖像的橫向分辨率一般遠低于軸向。

非線性熱應變模型可在較大的溫度變化范圍內保證溫度評估的準確性。本文的實驗結果表明，對于25℃的溫升，非線性模型的誤差精度基本可以滿足臨床的測溫需求。但對于活體脂肪組織，當熱療使其溫度達到70℃時，組織可能發生變性，非線性模型也就不再適用。如圖5所示，當溫度變化達到40℃以上時(組織溫度大于70℃)，雖然非線性模型的溫度評估誤差仍小于線性模型的評估誤差，但是其最大估計偏差已大于4℃，難以滿足臨床的要求。這可能與脂肪組織的嚴重變性(如液化)有關[18]。因此，對于變性的組織，無論線性還是非線性模型溫度評估模型都不再適用。鑒于此，Liu等[19]提出了自適應系數和紅外成像相結合的算法對消融治療中的溫度(治療溫度超過90℃)進行評估，取得了一定的效果。但該方法須借助紅外成像，難以應用在深部活體組織的測溫中，因而其應用受到限制。

圖5 溫度變化40℃的線性和非線性模型溫度估計Fig.5 Linear and non-linear model temperature estimation for temperature change of 40℃

對超聲熱應變測溫在活體中的應用，生理運動的影響必須考慮在內。此前有很多關于運動補償提高超聲測溫精度的研究。例如，基于在加熱之前采集一系列參考幀的運動補償算法[8,20]，利用靶區外的組織變形進行運動校正[13]等。這些算法在離體實驗中顯示出令人滿意的結果，但在活體中的結果并不理想。本文基于心跳運動的準周期性，通過動態幀的選取來降低心跳周期性的運動干擾。此外，對于活體呼吸(血流等)帶來的組織旋轉和形變[8]，簡單的平滑濾波和動態幀選取法并不足以解決運動步長問題，本文中采用自適應濾波得到了較好的結果。綜合非線性模型和活體運動抑制算法，本文中對25℃溫升的脂肪組織將測溫誤差控制在2.5℃以內，取得了令人滿意的效果。

4 結論

本文利用生物組織升溫過程中物理量變化的非線性，提出了用于活體溫度評估的非線性熱應變模型，對大范圍的溫度變化進行了有效的評估。采用動態選幀算法和自適應濾波相結合的方法來抑制活體的呼吸及心跳對溫度評估算法的干擾，并取得了很好的效果。實驗中以活體豬腎周脂肪為對象，對于小于25℃的溫升過程，測溫誤差不高于2.5℃。但對于更大范圍的溫度變化(大于40℃)，非線性溫度評估效果雖然比線性模型有明顯的提升，但還無法滿足臨床的需求。本文提出的非線性模型提升了超聲熱應變測溫的適用溫度范圍，拓展了超聲溫度評估技術的應用前景，有助于提高熱療相關的治療過程的安全性。