Известия РАН. Серия физическая, 2019, T. 83, № 1, стр. 70-75

ВВЕДЕНИЕ

Акустическая термотомография является перспективным методом неинвазивного исследования мягких биотканей [1–3]. Она основана на регистрации приемниками, находящимися вне исследуемого объекта, всей совокупности акустических полей, приходящих от данного объекта. Это собственное тепловое акустическое излучение объекта (зависящее от произведения собственной температуры объекта $T(\vec {r})$ и коэффициента вязкости $b(\vec {r}),$ пропорционального коэффициенту поглощения) и фоновое тепловое акустическое излучение, а также аналогичные им по интенсивности и спектру искусственно создаваемые поля (так называемая “подсветка”) [4]. Кроме того, рассеяние перечисленных видов акустических полей на неоднородностях фазовой скорости звука $\upsilon (\vec {r}) \equiv \frac{{{{\omega }^{2}}}}{{c_{0}^{2}}} - \frac{{{{\omega }^{2}}}}{{{{c}^{2}}(\vec {r})}}$ (где $\omega $ – круговая частота) и неоднородностях коэффициента поглощения приводит к образованию так называемых рассеянных полей. Здесь ${{c}_{0}} \equiv {\text{const}}$ и $c(\vec {r})$ – скорость звука вне и внутри объекта, соответственно. Тем самым, регистрируемые приемниками сигналы несут в себе информацию сразу о $T(\vec {r}),$ $b(\vec {r}),$ $c(\vec {r})$ [3–6]. Для задач медицинской диагностики целесообразно обеспечить раздельное восстановление этих параметров, что и обсуждается в представляемой работе.

Методу акустической термотомографии присуще чрезвычайно низкое входное отношение сигнал/помеха. Это связано с тем, что полезным является сигнал от элемента разрешения с линейным размером не более чем несколько длин волн, в то время как в роли помехи выступают сигналы, излучаемые всей областью томографирования. Данное обстоятельство, в сочетании с желаемым разрешением по температуре в несколько десятых долей градуса, приводит к необходимости увеличения времени накопления сигнала и применения антенных решеток с фокусировкой для съема данных [3, 5, 7, 8].

Для раздельного восстановления $T(\vec {r}),$ $b(\vec {r})$ и $c(\vec {r})$ в корреляционных термотомографических схемах необходимо иметь несколько линейно независимых наборов экспериментальных данных [3] для разных ракурсов приема сигналов. Принятые сигналы с целью их дальнейшей корреляционной обработки предварительно представляются в комплексном аналитическом виде, попарно перемножаются, и результат перемножения усредняется по времени накопления. Тем самым формируются элементы выборочной функции-матрицы когерентности ${{W}_{{ik}}}{\text{(}}\tau {\text{),}}$ т.е. выборочной корреляционной функции при комплексном представлении сигналов, для каждой рассматриваемой пары приемников с номерами $(i,k);$ τ – дискретизованная временнáя задержка прихода сигнала на второй приемник данной пары, по сравнению с приходом на первый приемник. При бесконечно большом времени накопления значения ${{W}_{{ik}}}{\text{(}}\tau {\text{)}}$ совпадают с элементами ${{\Gamma }_{{ik}}}{\text{(}}\tau {\text{)}}$ пространственно-временнόй функции-матрицы взаимной когерентности поля.

НЕОДНОРОДНОСТЬ МАЛОГО ВОЛНОВОГО РАЗМЕРА

Пусть сначала объект представляет собой одиночную рефракционно-поглощающую неоднородность с малым волновым размером, находящуюся в однородной непоглощающей фоновой среде. В этом случае действительная и мнимая части матрицы ${{\Gamma }_{{ik}}}{\text{(}}\tau {\text{),}}$ “сфазированной” в центр неоднородности ${{\vec {r}}^{p}},$ линейно пропорциональны (с точностью до главных членов) пространственной плотности мощности источников соответствующего типа [3, 6]:

Коэффициенты $A\left( {b({{{\vec {r}}}^{p}})} \right)$ и $C\left( {\upsilon ({{{\vec {r}}}^{p}})} \right)$ линейно пропорциональны, соответственно, коэффициенту вязкости $b({{\vec {r}}^{p}})$ и значению функции $\upsilon ({{\vec {r}}^{p}}),$ характеризующей неоднородность скорости звука в данной точке. Значения $T_{{bg}}^{{(i)}}({{\vec {r}}^{p}})$ и $T_{{bg}}^{{(k)}}({{\vec {r}}^{p}})$ – это температуры тех частей присутствующего в точке ${{\vec {r}}^{p}}$ фонового акустического излучения, которые после выхода из этого элемента разрешения распространяются в направлении i-го или, соответственно, k-го приемников и регистрируются ими. Под температурой фонового акустического излучения ${{T}_{{bg}}}({{\vec {r}}^{p}})$ подразумевается температура внешнего (по отношению к рассматриваемому элементу разрешения с центром в точке ${{\vec {r}}^{p}}$) излучения, присутствующего в данной точке ${{\vec {r}}^{p}}$ и создаваемого всеми полями, кроме поля, излучаемого рассматриваемым элементом разрешения. Таким образом, температура фонового излучения в элементе разрешения с центром ${{\vec {r}}^{p}}$ различна, в общем случае, “с точки зрения” различных направлений приема, т.е. анизотропна. Здесь и далее рассматривается режим разностных временных задержек ${{\tau }_{{ik;p}}} \equiv {{t}_{k}}({{\vec {r}}^{p}}) - {{t}_{i}}({{\vec {r}}^{p}}),$ при котором фазирование осуществляется путем компенсации относительной разности времен распространения сигнала ${{t}_{i}}({{\vec {r}}^{p}})$ и ${{t}_{k}}({{\vec {r}}^{p}})$ из точки ${{\vec {r}}^{p}}$ до соответствующего приемника текущей пары $(i,k).$

Для рассматриваемого объекта в виде неоднородности с малым волновым размером (или в виде тонкого слоя [5, 6]) значения $T_{{bg}}^{{(i)}}({{\vec {r}}^{p}})$ и $T_{{bg}}^{{(k)}}({{\vec {r}}^{p}})$ известны и совпадают с температурой внешнего поля, падающего на такой объект снаружи с соответствующего направления. Тогда в двух уравнениях (1) при фиксированном ${{\vec {r}}^{p}}$ присутствуют три неизвестные величины: $T({{\vec {r}}^{p}}),$ $b({{\vec {r}}^{p}}),$ $c({{\vec {r}}^{p}});$ для их определения предлагалось создание дополнительной (надфоновой) чисто анизотропной подсветки [4, 6, 7]. В случае же произвольного объекта пространственное распределение фоновой температуры становится неизвестным, поскольку собственное тепловое излучение объекта (от всех элементов разрешения, кроме рассматриваемого) также начинает влиять на фоновое поле и его температуру. В любом случае для нахождения $T({{\vec {r}}^{p}}),$ $b({{\vec {r}}^{p}})$ и $c({{\vec {r}}^{p}})$ термотомограф должен последовательно обеспечивать три режима работы.

В “исходном” состоянии томографической системы (первый режим) падающее на объект тепловое излучение совместно с собственным излучением объекта создает некоторое исходное пространственное распределение фоновой температуры ${{T}_{{bg}}}({{\vec {r}}^{p}}) \equiv {{T}_{{bg{\text{0}}}}}({{\vec {r}}^{p}}),$ зависящее, в общем случае, от точки ${{\vec {r}}^{p}}.$ Это распределение ${{T}_{{bg{\text{0}}}}}({{\vec {r}}^{p}})$ может быть как анизотропным, т.е. $T_{{bg{\text{0}}}}^{{(i)}}({{\vec {r}}^{p}}) \ne T_{{bg{\text{0}}}}^{{(k)}}({{\vec {r}}^{p}}),$ так и изотропным. В соотношения (1) при этом входят значения $T_{{bg{\text{0}}}}^{{(i)}}({{\vec {r}}^{p}})$ и $T_{{bg{\text{0}}}}^{{(k)}}({{\vec {r}}^{p}}).$ В других двух режимах создаются, соответственно, два типа дополнительной внешней анизотропной подсветки для каждой пары приемников $(i,k)$ и фиксированного элемента разрешения с центром ${{{\mathbf{r}}}^{p}}$ [3, 4, 6, 7]. Подсветка типа I проходит через данный элемент разрешения и дополнительно подсвечивает только приемник i; при этом приемник k дополнительно не подсвечен. Тогда, с точки зрения i-го и k-го приемников, температура фонового излучения в рассматриваемом элементе разрешения составляет $T_{{bg}}^{{{\text{I}}(i)}}({{\vec {r}}^{p}})$ и $T_{{bg}}^{{{\text{I}}(k)}}({{\vec {r}}^{p}})$ соответственно:

Здесь $\delta T_{{bg\_an}}^{{{\text{I}}(i)}}({{\vec {r}}^{p}})$ – чисто анизотропная добавка к температуре исходного фонового излучения за счет дополнительной анизотропной подсветки типа I. Подсветка типа II имеет противоположную асимметрию: дополнительно подсвечен только приемник k, но не приемник i, т.е. Подстановка выражений (2) в (1), а также подстановка выражений (3) в (1) дает: Индексы I и II у элементов пространственно-временнόй матрицы когерентности $\Gamma _{{ik}}^{I}(\tau )$ и $\Gamma _{{ik}}^{{ II}}(\tau )$ означают тип подсветки. Для обеспечения приемлемой точности результата обработки удобно использовать суммарные и разностные комбинации $\Gamma _{{ik}}^{I} \pm \Gamma _{{ik}}^{{II}}$ [3, 5]:

В рассматриваемом случае одиночной неоднородности с малым волновым размером легко обеспечить изотропность исходного фонового излучения, т.е. $T_{{bg{\text{0}}}}^{{(i)}}({{\vec {r}}^{p}})$ = $T_{{bg{\text{0}}}}^{{(k)}}({{\vec {r}}^{p}}) \equiv {{T}_{{bg{\text{0}}}}}({{\vec {r}}^{p}}),$ а также одинаковую анизотропную добавку к фоновой температуре $\delta T_{{bg\_an}}^{{{\text{I}}(i)}}({{\vec {r}}^{p}})$ = $\delta T_{{bg\_an}}^{{{\text{I}}I(k)}}({{\vec {r}}^{p}}) \equiv \delta {{T}_{{bg\_an}}}({{\vec {r}}^{p}}),$ поскольку рассеяние и поглощение полей отсутствует вне неоднородности. Тогда измерение в первом режиме, которое проводится в условиях чисто изотропного фонового излучения с известной температурой ${{T}_{{bg{\text{0}}}}}({{\vec {r}}^{p}}),$ дает для элементов матрицы ${{\Gamma }_{{ik}}}(\tau ) \equiv \Gamma _{{ik}}^{{isotr}}(\tau ),$ согласно (1):

Для второго и третьего режимов выражения (6) и (7) упрощаются, выделяя, соответственно, только диссипативный или только рефракционный вклад:

Соотношения (8)–(10) приводят к системе из трех уравнений, однозначно разрешаемой относительно неизвестных $b({{\vec {r}}^{p}}),$ $\upsilon ({{\vec {r}}^{p}}),$ $T({{\vec {r}}^{p}})$ при известных (для одиночной малой неоднородности) ${{T}_{{bg{\text{0}}}}}({{\vec {r}}^{p}})$ и $\delta {{T}_{{bg\_an}}}({{\vec {r}}^{p}})\,:$

НЕОДНОРОДНОСТЬ ПРОИЗВОЛЬНОГО РАЗМЕРА

Процедура раздельного восстановления искомых характеристик объекта с произвольными размерами усложняется. Во-первых, вклад в матрицу когерентности ${{\Gamma }_{{ik}}}({{\tau }_{{ik;p}}}),$ сфазированную на точку $\vec {r} = {{\vec {r}}^{p}}$ при фиксированной паре приемников $(i,k),$ вносят термоакустические источники во всех тех элементах разрешения с индексом $p1,$ для которых временнáя задержка ${{\tau }_{{ik;p1}}}$ на распространение сигнала до пары приемников такая же, как и для элемента с индексом p, т.е. ${{\tau }_{{ik;p1}}} = {{\tau }_{{ik;p}}}.$ Местоположение и количество элементов $p1$ зависит от конкретной используемой схемы акустического термотомографирования и существенно уменьшается с применением фокусировки. Если их число отлично от единицы, то для повышения точности оценки изображений целесообразно применять корреляционную обработку, строящуюся на методе максимального правдоподобия [9]. Во-вторых, каждый элемент разрешения теперь располагается в среде, неоднородной по скорости звука и поглощению (которые, к тому же, неизвестны), а восстановление информативных для диагностики значений $b({{\vec {r}}^{p}}),$ $\upsilon ({{\vec {r}}^{p}})$ и $T({{\vec {r}}^{p}})$ должно сопровождаться сопутствующей оценкой вспомогательных величин – фоновой температуры и анизотропной добавки к ней. В случае объекта, волновой размер которого не мал, эти вспомогательные величины неизвестны и зависят как от рассматриваемой точки внутри области томографирования, так и от направления, с которого излучение падает на приемник, согласно (1)–(5).

Наиболее прост в алгоритмическом плане вариант фокусирующих схем типа рассмотренных в [3, 7, 8]. Поскольку в них каждый элемент разрешения объекта сканируется отдельно от других, то играет роль только выделенный элемент с центром ${{\vec {r}}^{p}},$ а вклад от остальных элементов разрешения декоррелирован для данной задержки ${{\tau }_{{ik;p}}}.$ Благодаря этому изображения получаются по аналогии с соотношениями (1) фазировкой отдельного элемента выборочной матрицы ${{W}_{{ik}}}{\text{(}}\tau = {{\tau }_{{ik;p}}}{\text{),}}$ где в роли $(i,k)$ выступает пара сопряженных приемников, принимающих сигнал от конкретного элемента разрешения с центром в точке ${{\vec {r}}^{p}}$ [3, 7, 8]. Термоакустическое излучение от каждого элемента разрешения проходит через объект с неоднородными пространственными распределениями амплитудного коэффициента поглощения $\alpha (\vec {r})$ и фазовой скорости $c(\vec {r}).$ Чтобы учесть поглощение сигнала вдоль оси лучевой трубки ${{\vec {r}}^{p}} \to {{\vec {R}}_{i}},$ идущей от точки ${{\vec {r}}^{p}}$ до соответствующего эффективного центра (с радиус-вектором ${{\vec {R}}_{i}}$) текущего i-го приемника, вводится амплитудный коэффициент проникновения излучения $\kappa ({{\vec {r}}^{p}},i) \cong \exp \left[ { - \int_{{{{\vec {r}}}^{p}} \to {{{\vec {R}}}_{i}}} {\alpha (\vec {r}\,{\text{'}})d{{l}_{{\vec {r}\,{\text{'}}}}}} } \right].$ Неоднородности скорости звука учитываются на этапе расчета временнόй разностной задержки ${{\tau }_{{ik;p}}} \cong \int_{{{{\vec {r}}}^{p}} \to {{{\vec {R}}}_{k}}} {\frac{1}{{c(\vec {r}{\kern 1pt} {\text{'}})}}d{{l}_{{r'}}}} $ – $\int\limits_{{{{\vec {r}}}^{p}} \to {{{\vec {R}}}_{i}}} {\frac{1}{{c(\vec {r}{\kern 1pt} ')}}d{{l}_{{\vec {r}'}}}} $ при формировании ${{\Gamma }_{{ik}}}(\tau = {{\tau }_{{ik;p}}})$ в (1) или при обработке реальных данных ${{W}_{{ik}}}(\tau = {{\tau }_{{ik;p}}}).$ Тогда обработка отдельно $\operatorname{Re} {{W}_{{ik}}}(\tau )$ и отдельно $\operatorname{Im} {{W}_{{ik}}}(\tau )$ дает, по аналогии с (1):

где $T_{{bg}}^{{(ik)}}({{\vec {r}}^{p}}) \equiv {{\left\{ {T_{{bg}}^{{(i)}}({{{\vec {r}}}^{p}}) + T_{{bg}}^{{(k)}}({{{\vec {r}}}^{p}})} \right\}} \mathord{\left/ {\vphantom {{\left\{ {T_{{bg}}^{{(i)}}({{{\vec {r}}}^{p}}) + T_{{bg}}^{{(k)}}({{{\vec {r}}}^{p}})} \right\}} 2}} \right. \kern-0em} 2}$ и $\delta T_{{bg\_an}}^{{(ik)}}({{\vec {r}}^{p}})$ ≡ $ \equiv T_{{bg}}^{{(i)}}({{\vec {r}}^{p}}) - T_{{bg}}^{{(k)}}({{\vec {r}}^{p}}).$

Для раздельного восстановления характеристик среды требуется обеспечить те же три набора данных, как и в случае неоднородности с малым волновым размером. В условиях изотропного фонового акустического излучения с температурой, близкой к средней температуре исследуемого объекта, т.е. к $ \approx 37^\circ {\text{C}} = 310\,\,{\text{K}}$ (в медицинских приложениях), собственное тепловое излучение биоткани будет вносить искажения порядка $ \cong {{10}^{{ - 2}}};$ следовательно, температура итогового фонового поля ${{T}_{{bg{\text{0}}}}}$ ≡ ${{T}_{{bg(31{\text{0)}}}}} = 310\,\,{\text{K}}$ будет практически однородной и изотропной для всех элементов разрешения томографируемой области. При этом формируется выборочная матрица когерентности ${{W}_{{ik}}}(\tau ) \equiv W_{{ik}}^{{(310)}}(\tau )$ – первый набор данных; результатом его обработки, в соответствии с (12), является изображение $\Theta _{{Re}}^{{\text{0}}}({{\vec {r}}^{p}})\,:$

Второй и третий наборы данных в виде выборочных матриц когерентности $W_{{ik}}^{{\text{I}}}(\tau )$ и $W_{{ik}}^{{{\text{II}}}}(\tau )$ формируются, соответственно, добавлением к исходному однородному и изотропному фоновому излучению с температурой 310 K дополнительной чисто анизотропной подсветки типа I и типа II. Такая дополнительная подсветка создает изначально вне объекта в соответствующем направлении заданную анизотропную добавку $\delta {{T}_{0}}$ к температуре ${{T}_{{bg{\text{(310)}}}}}$ ($\delta {{T}_{0}}$ одинаково по величине для обоих наборов данных), а в элементе разрешения с центром в точке ${{\vec {r}}^{p}}$ создает анизотропную добавку $\delta T_{{bg\_an}}^{{{\text{I}}(i)}}({{\vec {r}}^{p}})$ или $\delta T_{{bg\_an}}^{{{\text{II}}(k)}}({{\vec {r}}^{p}})$ к той же температуре ${{T}_{{bg{\text{ (310)}}}}}.$ Тогда фоновая температура в рассматриваемом элементе разрешения с точки зрения i-го и k-го приемников, по аналогии с (2) и (3), такова: $T_{{bg}}^{{{\text{I}}(i)}}({{\vec {r}}^{p}})$ = ${{T}_{{bg{\text{(310)}}}}} + \delta T_{{bg\_an}}^{{{\text{I}}(i)}}({{\vec {r}}^{p}}),$ $T_{{bg}}^{{{\text{I}}(k)}}({{\vec {r}}^{p}}) = {{T}_{{bg(310)}}};$ $T_{{bg}}^{{{\text{II}}(i)}}({{\vec {r}}^{p}}) = {{T}_{{bg{\text{(310)}}}}},$ $T_{{bg}}^{{{\text{II}}(k)}}({{\vec {r}}^{p}})$ = ${{T}_{{bg{\text{(310)}}}}} + \delta T_{{bg\_an}}^{{{\text{II}}(k)}}({{\vec {r}}^{p}}).$ С учетом этого, а также соотношений (12) и (13), по аналогии с (4) и (5) получается:

Обработка матриц $W_{{ik}}^{{\text{I}}} - W_{{ik}}^{{(310)}}$ и $W_{{ik}}^{{{\text{II}}}} - W_{{ik}}^{{(310)}}$ исключает в (14) и (15) влияние как собственной температуры объекта $T(\vec {r}),$ так и температуры исходного изотропного фонового излучения ${{T}_{{bg{\text{(310)}}}}},$ оставляя вклад только от дополнительной анизотропной подсветки.

Чтобы связать неизвестные анизотропные добавки $\delta T_{{bg\_an}}^{{{\text{I}}(i)}}({{\vec {r}}^{p}})$ и $\delta T_{{bg\_an}}^{{{\text{II}}(k)}}({{\vec {r}}^{p}})$ с $\delta {{T}_{0}},$ вводится, по аналогии с $\kappa ({{\vec {r}}^{p}},i),$ эффективный энергетический коэффициент проникновения дополнительного подсвечивающего внешнего излучения в точку ${{\vec {r}}^{p}}\,:$ ${{\kappa }_{{en}}}({{\vec {r}}^{p}},S)$ ≅ $\exp \left[ { - 2\int_{S \to {{{\vec {r}}}^{p}}} {\alpha (\vec {r}{\kern 1pt} {\text{'}})d{{l}_{{\vec {r}'}}}} } \right].$ Интегрирование ведется вдоль оси лучевой трубки, связывающей подсвечивающий излучатель ($S = {{S}^{I}}$ и $S = {{S}^{{II}}}$ для подсветки типа I и типа II, соответственно) и центр элемента разрешения ${{\vec {r}}^{p}}.$ Тогда $\delta T_{{bg\_an}}^{{{\text{I}}(i)}}({{\vec {r}}^{p}})$ = $ = \kappa _{{en}}^{I}({{\vec {r}}^{p}},{{S}^{I}})\delta {{T}_{0}},$ $\delta T_{{bg\_an}}^{{{\text{II}}(k)}}({{\vec {r}}^{p}})$ = $\kappa _{{en}}^{{II}}({{\vec {r}}^{p}},{{S}^{{II}}})\delta {{T}_{0}},$ и соотношения (14) и (15) приобретают вид:

где $\Theta _{{Re}}^{I}({{\vec {r}}^{p}}),$ $\Theta _{{Im}}^{I}({{\vec {r}}^{p}}),$ $\Theta _{{Re}}^{{II}}({{\vec {r}}^{p}}),$ $\Theta _{{Im}}^{{II}}({{\vec {r}}^{p}})$ – получаемые изображения.

Из (13), (16) и (17) следует, что для раздельного восстановления искомых характеристик среды достаточно построить изображение $\Theta _{{Re}}^{{\text{0}}}({{\vec {r}}^{p}})$ и построить изображения ${{\Theta }_{{Re}}}({{\vec {r}}^{p}})$ и ${{\Theta }_{{Im}}}({{\vec {r}}^{p}})$ только для одного типа дополнительной анизотропной подсветки. Тем не менее, при этом необходима очень высокая точность фазирования выборочных матриц, что требует априорного знания пространственных распределений $c(\vec {r})$ и $\alpha (\vec {r}).$ Однако на практике эти распределения заранее неизвестны, и именно их требуется определить наряду с температурными характеристиками. Фазировать выборочные матрицы приходится, в общем случае, на основе лишь предварительно оцененных значений этих распределений. Поэтому для обеспечения приемлемого качества получаемых изображений нужны корреляционные данные от обоих типов дополнительной анизотропной подсветки. Это позволяет, наподобие (11), использовать суммарные и разностные комбинации выборочных матриц когерентности двух типов, что усиливает полезный эффект и повышает точность восстановления. Так, из (16) и (17) вытекает, что двум комбинациям

соответствуют два изображения ${{\Theta }^{{sum}}}({{\vec {r}}^{p}})$ и ${{\Theta }^{{dif}}}({{\vec {r}}^{p}})\,:$ Изображения ${{\Theta }^{{sum}}}({{\vec {r}}^{p}})$ и ${{\Theta }^{{dif}}}({{\vec {r}}^{p}})$ отражают, соответственно, диссипативные и рефракционные свойства исследуемого объекта.

ЗАКЛЮЧЕНИЕ

Соотношения (13), (18) и (19) являются основными для определения собственной температуры и акустических характеристик на основе изображений $\Theta _{{Re}}^{{\text{0}}},$ ${{\Theta }^{{sum}}}$ и ${{\Theta }^{{dif}}},$ полученных в режиме разностных задержек. Тем не менее, даже наличие всех корреляционных данных, обеспечиваемых тремя описанными режимами томографирования, не избавляет от необходимости применения итерационной процедуры. На каждом итерационном шаге уточняются оценки искомых температурных и акустических характеристик, а также оцениваются коэффициенты κ и ${{\kappa }_{{en}}},$ необходимые для следующего шага [3].

Работа поддержана грантом РФФИ № 16-29-02097 офи_м.