Акустический журнал, 2022, T. 68, № 5, стр. 562-570

ВВЕДЕНИЕ

В ряде медицинских приложений, связанных с локальной гипертермией тканей организма человека, важно с надлежащей точностью осуществить безболезненные измерения глубинной температуры. Возможно, магнитно-резонансная термометрия в будущем позволит решить эту задачу [1]. Однако этот метод требует дорогого оборудования, обученного персонала и специально подготовленных помещений. Поэтому наличие альтернативных методов крайне актуально, если они обеспечивают надлежащую точность. Для медицинских приложений ошибка определения температуры не должна превышать 0.5–1 К, пространственное разрешение должно быть не хуже 5 мм.

Предлагается использовать пассивную акустическую термометрию, физической основой которой является регистрация теплового акустического излучения объекта [2–4]. Измерения шумового сигнала приводят к значительному времени интегрирования: в мегагерцовом диапазоне для получения требуемой точности необходимо усреднять сигнал в течение 30–50 с. Чтобы снизить это время без потери точности, предполагается при восстановлении температуры использовать уравнение теплопроводности с кровотоком (уравнение Пеннеса [5]). В акустотермографии этот подход уже исследовался теоретически в различных модификациях [6, 7]. В предлагаемом алгоритме восстанавливается не сама температура, а параметры уравнения теплопроводности, при этом искомые параметры не определяются каждый раз заново, а уточняются по мере измерений. Такой подход уже рассмотрен, в том числе и экспериментально. При нагреве модельного объекта – говяжьей печени – восстанавливались температуропроводность и источник нагрева [8], при охлаждении нагретых модельных объектов – пластилина и тефлона – восстанавливались температуропроводность и начальная температура [9, 10]. Новым является учет в этой схеме кровотока в теле человека.

В работе будет рассмотрен одномерный случай, с одной стороны, наиболее физически прозрачный, с другой стороны, наиболее некорректный. Если восстанавливать профиль температуры по результатам измерения теплового излучения на разных частотах, то различить расположенные на одной глубине узкий сильно нагретый слой и широкий, но слабо нагретый, – невозможно при существующей точности измерений.

Отметим работу [11] в области СВЧ-термометрии (математически СВЧ- и акустотермометрия описываются схожими уравнениями), где для восстановления температуры было использовано уравнение теплопроводности (без кровотока) и рассмотрен внешний нагрев объекта, что соответствует учету граничных условий. В данной работе предлагается рассмотреть локальный нагрев в глубине тела человека, т.е. ввести в уравнение теплопроводности источник.

Таким образом, в работе предложен алгоритм восстановления профиля температуры при локальном нагреве ткани с учетом уравнения теплопроводности с кровотоком. Предложенный алгоритм опробован экспериментально.

МАТЕРИАЛЫ И МЕТОДЫ

Расчетная модель. Рассмотрим одномерную модель локальной гипертермии мягких тканей человека (см. рис. 1а). Изменения температуры $T$ подчиняются уравнению теплопроводности с учетом кровотока [5]:

где ${{a}^{2}}$ – коэффициент температуропроводности, ${{\eta }}$ – кровоток, ${{T}_{0}}$ = 37°С – температура притекающей крови, $S\left( x \right)$ – источник нагрева, не меняющийся во времени, $x$ – координата, направленная вглубь тела, $t$ – время. Будем считать, что исследуемая среда однородна по температуропроводности и кровотоку ${{a}^{2}} = $ const, ${{\eta }} = $ const. Начальное распределение температуры: $T\left( {t = 0,x} \right) = {{T}_{0}}$, граничные условия: $T\left( {t,x = 0} \right)$ $ = T\left( {t,x = \infty } \right) = {{T}_{0}}$.

Акустояркостная температура, измеряемая датчиком, расположенным на границе тела при $x = 0$, рассчитывается по формуле [12]:

где ${{\gamma }}$ – коэффициент поглощения ультразвука по интенсивности. Будем считать, что исследуемая среда однородна по поглощению ${{\gamma }} = $ const.

Умножим обе части уравнения (1) на выражение ${{\gamma }}\exp \left( { - {{\gamma }}x} \right)$ и проинтегрируем по х. При этом первое слагаемое правой части (1) дважды проинтегрируем по частям. В результате, с учетом выражения (2), получим:

Предположим, что источник нагрева расположен достаточно глубоко и время нагрева невелико, так что влиянием источника на температурный градиент на поверхности тела можно пренебречь: ${{\partial T\left( 0 \right)} \mathord{\left/ {\vphantom {{\partial T\left( 0 \right)} {\partial x \approx 0}}} \right. \kern-0em} {\partial x \approx 0}}$.

Окончательно получаем обыкновенное дифференциальное уравнение относительно акустояркостной температуры

решение которого определяется выражением

В случае слабого кровотока ${{\eta }} = 0$ при небольшом времени нагрева $t \ll {1 \mathord{\left/ {\vphantom {1 {{{\gamma }^{2}}{{a}^{2}}}}} \right. \kern-0em} {{{\gamma }^{2}}{{a}^{2}}}}$ акустояркостная температура линейно растет со временем:

Алгоритм восстановления температуры. Исходя из выражений (5) и (6), можно предложить алгоритм восстановления меняющегося во времени профиля температуры: по временной зависимости измеряемой акустояркостной температуры определяются параметры уравнения теплопроводности и затем, с помощью этого уравнения, вычисляется профиль температуры в каждый момент времени. Предположим, что коэффициент температуропроводности, профиль источника и два его параметра: ширина d (взятая по половине от максимального значения) и положение x₀, а также коэффициент поглощения ультразвука известны. Необходимо найти неизвестные параметры уравнения теплопроводности – кровоток ${{\eta }}$ и амплитуду источника Q.

Для примера рассмотрим источник в форме гауссиана

Тогда интеграл, представленный в правой части выражения (4), равен (при ${{x}_{0}} \gg d$)

Для восстановления температуры предлагается следующий алгоритм (см. рис. 2). Входные параметры алгоритма: ${{a}^{2}}$, d, ${{x}_{0}}$, ${{\gamma }}$ и экспериментальная временная зависимость акустояркостной температуры ${{T}_{A}}\left( t \right)$, которая аппроксимируется параболой ${{T}_{A}}\left( t \right) \approx {{T}_{0}} + {{k}_{1}}t + {{k}_{2}}{{t}^{2}}$, проходящей через точку ${{T}_{A}}\left( 0 \right) = {{T}_{0}}$, ${{k}_{1}}$ и ${{k}_{2}}$ – рассчитываемые коэффициенты. Если полученная парабола вогнута или близка к прямой, т.е. ${{k}_{2}} \geqslant 0$, что говорит о слабом кровотоке ${{\eta }} \approx 0$, то акустояркостная температура аппроксимируется прямой ${{T}_{A}}\left( t \right) \approx {{T}_{0}} + {{k}_{3}}t$. При этом выходные параметры $Q = \frac{{{{k}_{3}}\exp \left( {\gamma {{x}_{0}} - {{{{d}^{2}}{{\gamma }^{2}}} \mathord{\left/ {\vphantom {{{{d}^{2}}{{\gamma }^{2}}} 8}} \right. \kern-0em} 8}\ln 2} \right)}}{{\gamma d}}$ и ${{\eta }}$ = 0. Если полученная парабола выпуклая, т.е. ${{k}_{2}} < 0$, то акустояркостная температура аппроксимируется выражением ${{T}_{A}}\left( t \right) \approx {{T}_{0}} + \frac{{{{k}_{4}}}}{{\eta *}}\left[ {1 - \exp \left( { - \eta *t} \right)} \right]$. При этом $Q = \frac{{{{k}_{4}}\exp \left( {\gamma {{x}_{0}} - {{{{d}^{2}}{{\gamma }^{2}}} \mathord{\left/ {\vphantom {{{{d}^{2}}{{\gamma }^{2}}} 8}} \right. \kern-0em} 8}\ln 2} \right)}}{{\gamma d}}$ и $\eta = \eta {\text{*}} + \,\,{{\gamma }^{2}}{{a}^{2}}$. По параметрам $Q$ и ${{\eta }}$ по уравнению (1) рассчитывается меняющееся во времени температурное распределение.

Для оценки точности алгоритма определяются погрешности восстановления максимальной температуры ${{T}_{{{\text{max}}}}} = T\left( {{{x}_{0}}} \right)$ и размера нагретой области h, взятого в точках, где увеличение температуры равно половине от максимального $T\left( {{{x}_{0}} - {h \mathord{\left/ {\vphantom {h 2}} \right. \kern-0em} 2}} \right) - {{T}_{0}} = $ $T\left( {{{x}_{0}} + {h \mathord{\left/ {\vphantom {h 2}} \right. \kern-0em} 2}} \right) - {{T}_{0}} = $ ${{\left( {{{T}_{{\max }}} - {{T}_{0}}} \right)} \mathord{\left/ {\vphantom {{\left( {{{T}_{{\max }}} - {{T}_{0}}} \right)} 2}} \right. \kern-0em} 2}$. Чтобы определить эти погрешности, на точное значение временной зависимости акустояркостной температуры накладывается нормально распределенная погрешность со стандартным отклонением ${{\delta }}{{T}_{A}}$. Эта погрешность определяется пороговой чувствительностью акустотермографа. По заданной с погрешностью акустояркостной температуре с помощью приведенного выше алгоритма определяются параметры $Q$ и ${{\eta }}$, по уравнению (1) рассчитывается температура и определяются максимальная температура и размер нагретой области. Указанная процедура повторялась 1000 раз. Это дало возможность статистически определить погрешности восстановления (стандартные отклонения) для максимальной температуры и размера нагретой области.

Экспериментальное опробование алгоритма. Алгоритм восстановления температуры был опробован при нагреве модельного объекта – пластизоля (пластизоль прозрачный, твердость 15–17, “Альпина Пласт”, Клин, Россия), акустические и теплофизические свойства которого близки свойствам мягких тканей тела человека [13, 14]. В качестве нагревателя использовалось помещенное в пластизоль, на расстоянии 50 ± 3 мм от поверхности, сопротивление 400 Ом, на которое подавали напряжение 18 В. Температура пластизоля контролировалась двумя цифровыми термометрами DS18S20P (Maxim Integrated, Сан-Хосе, США) с точностью 0.3 К. Датчики располагались на расстояниях 20 ± 3 и 41 ± 3 мм от границы пластизоля (см. рис. 1б).

Для измерений теплового акустического излучения был использован многоканальный акустотермограф [15–17], разработанный в ИПФ РАН А.Д. Мансфельдом и в настоящее время дорабатываемый Р.В. Беляевым (полоса пропускания 1.6–2.5 МГц, пороговая чувствительность при времени интегрирования 10 с – 0.2 К). Принимаемые акустические сигналы преобразовывались в электрические, которые усиливались, проходили через квадратичный детектор и усреднялись в течение 30 мс. С выходов многоканального акустотермографа сигналы подавались на 14-ти разрядный многоканальный АЦП Е14-140 (ЗАО “L-Card”, Москва, Россия) с частотой дискретизации 1 КГц на один канал. Разработанная программа проводила дальнейшее усреднение данных.

До измерений акустический датчик находился в держателе, который являлся акустическим черным телом, при комнатной температуре. Измерения теплового акустического излучения пластизоля проводились в течение 10–15 с. Для акустического контакта на поверхность объекта наносили гель для УЗИ “Медиагель”. После измерений датчик возвращали в держатель до следующего измерения. Интервал между измерениями составлял около 1 мин.

РЕЗУЛЬТАТЫ

Выбор параметров. В норме кровоток в разных тканях организма меняется в широких пределах [14] от 28–38 мл/(кг мин) или (4.7–6.3) × 10^–4 1/с в жировой ткани–скелетных мышцах до 1050–5000 мл/(кг мин) или (1.75–8.33) × 10^–2 1/с в печени–щитовидной железе. Коэффициент температуропроводности меняется в тканях незначительно [14], например, при 37°С температуропроводность в печени составляет 0.141, а в жировой ткани 0.131 мм²/с. Изменения этого параметра при изменении температуры также невелики: например, в печени – 0.28% на градус. Для оценки размера источника воспользуемся следующими соображениями. Предположим, что нагрев осуществляется с помощью введенного в организм световода, по которому передается ИК излучение. Это излучение поглощается в глубине организма, где и происходит нагрев. Согласно [14] характерная глубина проникновения (расстояние, на котором интенсивность света уменьшается в е раз) в печени на длинах волн 800 и 1000 нм составляет 1 и 0.5 мм, соответственно, а в почках на тех же длинах волн – 2.7 и 1.6 мм соответственно. Отметим, что источник при поглощении излучения не является гауссианом, однако расчеты показывают, что форма источника незначительно влияет на форму температурного распределения. Частотную зависимость коэффициента ослабления ультразвука в мягких тканях принято задавать эмпирической формулой: ${{\gamma }} = a{{f}^{b}}$ [14]. Например, для печени и мышц бедра коэффициент b = 1, коэффициент a = 0.162 и 0.128 1/(мм МГц), соответственно. В этом случае, на частоте 2 МГц коэффициент ослабления равен 0.0324 и 0.0256 1/мм соответственно.

После анализа литературы для рассчитываемой модели были выбраны следующие параметры уравнения теплопроводности: ${{x}_{0}}$ = 25 мм, d = 4 мм, a² = 0.14 мм²/с, ${{\gamma }}$ = 0.03 1/мм.

На рассчитанную акустояркостную температуру накладывалась нормально распределенная случайная ошибка, соответствующая пороговой чувствительности акустотермометра 0.2 К за 10 с [18].

Влияние кровотока на температуру показано на рис. 3, где представлены профили температуры в разных тканях с разным кровотоком при источнике в форме гауссиана. Пусть задача врача – нагреть некоторую область за 5 мин до оптимальной температуры, например, до 43°С. Отметим, что в случаях без кровотока и для скелетных мышц в представленном масштабе температурные профили практически совпадают. Ширина нагретой области при увеличении кровотока более чем на два порядка (скелетные мышцы–щитовидная железа) уменьшается незначительно – на 6 мм. При значительном кровотоке ширина температурного распределения стремится к ширине источника.

Расчет акустояркостной температуры. Временные зависимости акустояркостной температуры, рассчитанные по формуле (2), в случае отсутствия кровотока и при кровотоке как в печени показаны на рис. 4. Для примера показаны 6 реализаций акустояркостной температуры с погрешностью ${{\delta }}{{T}_{A}}~$ = 0.2 К. Из рис. 4а видно, что акустояркостная температура без кровотока растет линейно, при этом рост максимальной температуры замедляется, однако увеличивается размер нагретой области. Акустояркостная температура определяется интегралом от распределения глубинной температуры (т.е. пропорциональна произведению максимальной температуры на размер нагретой области) и поэтому на данном временном промежутке линейно увеличивается во времени.

Восстановление температурного распределения. Первый раз восстановление производилось через 50 с после начала нагрева и повторялось каждые 50 с в течение пяти минут. Если в первый раз были использованы данные только за первые 50 с нагрева, то в дальнейшем это время увеличивалось: во второй раз восстанавливали температуру за 100 с, в третий – за 150 с и т.д. В последнем восстановлении использовались все полученные за 5 мин данные.

Задача восстановления максимальной температуры без кровотока является линейной, и погрешность ее решения ${{\delta }}{{T}_{{{\text{max}}}}}$ пропорциональна выражению $\delta {{T}_{{\max }}}\sim {{\delta {{T}_{A}}} \mathord{\left/ {\vphantom {{\delta {{T}_{A}}} {\sqrt {\sum {{t}_{i}}^{2}} }}} \right. \kern-0em} {\sqrt {\sum {{t}_{i}}^{2}} }}$, где ${{t}_{i}}$ – моменты измерения акустояркостной температуры. Таким образом, для ее вычисления не требовалось рассчитывать реализации случайного процесса. Тем не менее, это было сделано для того, чтобы удостовериться в предлагаемом подходе перед его использованием для случая сильного кровотока. Погрешность восстановления температуры показана на рис. 4а. При интегрировании данных за 50 с погрешность составила 0.8 К, за 150 с – снизилась до 0.4 К, при интегрировании за все время измерений – составила 0.25 К.

Без кровотока ширина температурного распределения зависит от времени, но при заданном времени не зависит от амплитуды источника. Поэтому восстановленный размер температурного распределения совпадает с исходным.

Восстановление при кровотоке как в печени. Задача восстановления максимальной температуры и ширины нагретой области при наличии кровотока является нелинейной и была решена численно. Погрешность восстановления максимальной температуры показана на рис. 4б. При интегрировании данных за 50 с погрешность составила 1.2 К, за 150 с – снизилась до 0.45 К, при интегрировании за все время измерений составила 0.27 К.

Эксперимент с пластизолем. Пластизоль нагревали в течение 33 мин. Измеренные значения акустояркостной температуры и данные электронных термометров (до, во время и после нагрева) представлены на рис. 5. Исходная температура пластизоля была равна 20.1°С. Отметим, что акустояркостная температура и температура, измерявшаяся на расстоянии 20 мм от поверхности пластизоля (т.е. в 30 мм от центра нагрева), не сразу начали снижаться после окончания нагрева. По результатам измерений был определен коэффициент k₃ = 0.0016 К/с и рассчитана амплитуда источника Q = 0.05 К/с (при ${{a}^{2}} = $ 0.15 мм²/с [13], d = 2 мм, ${{x}_{0}} = $ 50 мм, ${{\gamma }} = $ 0.03 1/мм [13]). С помощью уравнения (1) были рассчитаны временные зависимости температуры в местах расположения электронных термометров: температурные диапазоны в пределах ± 3 мм показаны на рис. 5. Из графиков видно, что предложенный алгоритм работает.

ОБСУЖДЕНИЕ

Из трех параметров источника: глубины залегания, ширины и амплитуды восстанавливался только последний. Это связано с тем, что остальные параметры можно определить предварительно, до начала нагрева. Если нагрев осуществляется с помощью введенного в организм световода, по которому передается ИК излучение, то положение кончика световода контролируется с помощью стандартного медицинского УЗИ. Ширину источника можно измерить предварительно в модельном эксперименте. Теоретически в эксперименте на модели можно измерить и амплитуду источника (в этом случае в предлагаемом алгоритме остается только определить кровоток), но реально эта характеристика зависит от состояния кончика световода, от того, с какой тканью он контактирует. При нагреве может образоваться обгоревшая ткань, которая препятствует прохождению излучения из световода в окружающую ткань, что невозможно проконтролировать заранее [19].

Выбор максимальной температуры связан с использованием гипертермии в онкологии. В работе [20] для оценки тепловой дозы в условиях, которые клинически значимы для биологического эффекта, была предложена эталонная температура 43°C. В этом случае тепловая доза, выраженная в минутах, равна времени нагрева.

В качестве параметров восстановления взяты максимальная температура и размер нагретой области. Это связано с тем, что именно эти параметры удобно использовать при контроле гипертермии: врачам, в первую очередь, важна информация о размере и температуре нагретой области, но не о распределении температуры.

Отметим фундаментальность полученного результата: как известно, по данным об акустояркостной температуре в данный момент времени невозможно восстановить температурное распределение – по точке нельзя восстановить функцию. Это можно сделать, если измерена частотная зависимость акустояркостной температуры, которая связана с зависимостью коэффициента поглощения от частоты [21]. В данном исследовании проблема решается другим способом: рассматривается временная зависимость акустояркостной температуры, и по ней восстанавливаются параметры уравнения теплопроводности. После интегрирования уравнения восстанавливается температурное распределение, меняющееся во времени. Требование, чтобы температура удовлетворяла уравнению теплопроводности, приводит к повышению точности восстановления. Каждый раз (через 50, 100, …, 300 с) искомые параметры не определяются заново, а лишь уточняются (см. [5, 7]). Практически увеличивается время интегрирования, но при этом не страдает временное разрешение – получать новые результаты можно через любой промежуток времени. Это преимущество достигается за счет требования, чтобы параметры уравнения теплопроводности оставались постоянными. Отметим, что данное требование, вообще говоря, неприменимо к кровотоку – при нагреве ткани он увеличивается. В этом случае необходимо найти неизменный параметр меняющегося кровотока, например, считать, что кровоток растет линейно [22]. Это является темой отдельного исследования.

Отметим, что форма источника не является критичной для восстановления температуры. Если сравнить источники прямоугольной формы и гауссиан, расположенные на одной глубине, имеющие одинаковую ширину и площадь по кривой источника, то при пятиминутном нагреве, при одной и той же максимальной температуре 43°C, размеры нагретых областей различаются меньше, чем на 1 мм, а акустояркостные температуры различаются менее, чем на 0.03°С.

Рассмотрим пределы применимости модели, связанные с граничными условиями. Существенное ограничение – нулевой градиент температуры на границах области. Это означает, что все тепло остается в области, что можно считать справедливым только в начале нагрева. Таким образом, при слабом кровотоке модель можно считать адекватной, пока акустояркостная температура растет линейно со временем.

Укажем ограничения модели, связанные с акустической неоднородностью среды. Поглощение ультразвука различно в разных тканях организма. Если изменения коэффициента поглощения коррелируют с изменениями температуры, то акустояркостная температура меняется. Например, если в нагретой области поглощение больше среднего, то это увеличит акустояркостную температуру. Однако, согласно [14] температурная зависимость коэффициента ослабления ультразвука в мягких тканях при температуре ~37°C на частоте ~2 МГц практически отсутствует: температурный коэффициент равен 0 ± 0.4 (м К)^–1. Если же поглощение не зависит от температуры, то пространственная неоднородность поглощения не приводит к существенным изменениям акустояркостной температуры из-за интегрального характера этой характеристики.

В различных мягких тканях организма скорость звука различна. Согласно [14] скорость звука в печени и в мышечной ткани меняется в диапазоне 1568–1593 м/с или на ±1.6%. Эти вариации скорости приводят к изменению аппаратной функции (диаграммы направленности) приемного устройства, что необходимо учесть в 3D модели.

Отметим иллюстративный характер экспериментальной апробации алгоритма. Условия эксперимента не полностью соответствовали параметрам численного расчета: в шесть раз большее время нагрева; локальный источник – не нагретый слой, параллельный поверхности; наличие отражения от сопротивления, помещенного в пластизоль. Однако, задача экспериментальной проверки точности предложенного алгоритма не стояла. Актуальным является использование предложенного подхода к решению 3D обратной задачи, причем с учетом аппаратной функции акустического датчика [23]. Это – тема будущего исследования.

ЗАКЛЮЧЕНИЕ

Предложенный алгоритм позволяет восстановить температуру и размер нагретой области с точностью, приемлемой для медицинских приложений. Эта точность достигается через две минуты после начала нагрева, после чего температурное распределение можно контролировать практически непрерывно. Этот результат достигается благодаря использованию уравнения теплопроводности. Проведена экспериментальная проверка алгоритма.

Работа выполнена при поддержке РФФИ (грант № 20-02-00759), а также в рамках государственного задания Института радиотехники и электроники им. В.А. Котельникова РАН (№ государственной регистрации АААА-А19–119041590070-01).