автореферат диссертации по информатике, вычислительной технике и управлению, 05.13.18, диссертация на тему:Численное моделирование процессов в лазерной оптической медицинской томографии

Московский физико-технический институт (государственный университет)

РГб од

На правах рукописи

Верблюденко Павел Александрович

УДК 53.072; 53:681.3 517.958:536.71

Численное моделирование процессов в лазерной оптической медицинской томографии

05.13.18 — "Теоретические основы математического моделирования, численные методы и комплексы программ"

АВТОРЕФЕРАТ диссертации на соискание ученой степени кандидата физико-математических наук

Москва 1997 год

Работа выполнена в Московском физико-техническом институте

Научный руководитель — доктор физико-математических наук,

профессор Сон Эдуард Евгеньевич

Консультант — кандидат физико-математических наук

Калюжный Владимир Викторович

Официальные оппоненты —

академик, доктор физико-математических наук Гуляев Юрий Васильевич кандидат физико-математических наук Сергиевская Алла Леонидовна

Ведущая организация — Международный лазерный центр МГУ

Защита состоится "_"_1997 г. в_часов на заседании диссертационного совета К 063.91.03 при Московском физико-техническом институте (государственном университете) по адресу: 141700 Московская обл., г. Долгопрудный, Институтский пер. 9, МФТИ.

С диссертацией можно ознакомиться в библиотеке МФТИ

Автореферат разослан "_"_1997 г.

Ученый секретарь диссертационного совета

_ д. ф-м. н., проф. Самыловский Александр Иванович

ОБЩАЯ ХАРАКТЕРИСТИКА РАБОТЫ

Эффективность лечебно-профилактической работы в современной медицине определяется прежде всего совершенством диагностической техники. Главная цель научных исследований в этой области — это разработка методов и аппаратных средств, позволяющих обеспечить максимальный объем точной информации о состоянии и поведении отдельных органов, тканей и всего организма человека в целом при минимизации вредных для организма последствий проведения обследования.

Одно из перспективных направлений развития медицинской диагностики заключается в использовании различных видов излучений (электромагнитное, СВЧ, УЗ и другие). Основное внимание при развитии методов лучевой диагностики должно уделяться корректной обработке и верификации получаемых результатов при одновременном соблюдении всех необходимых критериев безопасности организма и отдельных органов.

Методы лучевой диагностики различаются по степени и сложности математической (в частности компьютерной) обработки получаемых с их помощью данных. Для одних методов (например рентгенография) характерным является то, что полученные в процессе обследования данные о пациенте фактически являются конечными и не подвергаются в дальнейшем какой-либо существенной математической обработке. Однако развитие вычислительной математики и компьютерной техники привели к тому, что появилась возможность подвергать полученные при обследовании данные самой разнообразной математической обработке и, таким образом, существенно увеличить количество получаемой при обследовании информации. Эта возможность была успешно использована в большинстве методов лучевой диагностики, среди которых по распространенности и перспективности выделяются томографические методы (от греч. tomos — часть, слой и grapho — писать), основанные на получении изображения внутреннего строения объекта в тонком слое.

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

Совершенствование лучевой диагностики связано с поиском источника такого детектируемого излучения, которое как можно лучше удовлетво-

ряло бы требованиям по информативности и безопасности обследования организма. Здесь широкое поле для исследований представляют оптические методы диагностики, в которых в качестве детектируемого излучения используется излучение в дальнем, ближнем ИК и видимом диапазонах. Излучение этого диапазона, являясь по своей сути тепловым излучением, не ионизирует ткани, хорошо переносится организмом человека. При этом оптические свойства биотканей очень чувствительны не только к их структуре, но и к функциональному состоянию (например кровенаполнен-ности), что говорит о возможности широкого применения и новых приложениях этого метода.

Систематическое теоретическое обоснование переноса оптического излучение через рассеивающие среды дано в работах A.Ishimaru. Сегодня теоретическим и экспериментальным исследованиям этого процесса посвящены многочисленные работы таких зарубежных авторов, как B.C. Wilson, S. L. Jacques, B.Chance, R.R. Alfano, D.A. Benaron, D.K. Stevenson, W. Cheong, S. A. Prahl, A. J. Weich, L.O.Svaasand, J.Beuthan. В нашей стране и странах СНГ проблемам лазерной диагностики и оптической томографии посвящены работы В.В.Тучина, А.В.Приезжева, B.C. Летохова,

A.Я.Хайруллиной и многих других исследователей. Систематическое изложение математических методов обработки томографических данных применительно к рентгеновской томографии дано в работах G.T.Herman. Методы решения прямых и обратных задач переноса излучения глубоко исследованы в работах отечественных авторов: А.Н.Тихонова

B.Я.Арсенина Г.И.Марчука, В.И.Лебедева, Н.Н.Яненко, A.A.Самарского, Ю.В.Ракитского, A.B. Гончарского.

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

Целью диссертации является разработка и обоснование эффективного численного метода решения краевой задачи для уравнения лучистого переноса в биотканях, позволяющего достаточно полно учесть эффекты мно-

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

Для достижения поставленной цели решались следующие задачи:

1. Проведен анализ существующих томографических методов медицинской диагностики, их достоинств и недостатков, сделан вывод о путях дальнейшего развития методов оптической томографии.

2. Разработан эффективный численный метод решения краевой задачи для уравнения лучистого переноса, учитывающий эффекты многократного рассеяния, применимый для решения прямых и обратных задач лучистого переноса в различных средах, в том числе, и в биотканях.

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

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

5. Проведено расчетно-параметрическое исследование с целью определения возможности использования оптического излучения для медицинской диагностики нейрофизиологических процессов.

6. Сформулирована обратная задача реконструкции изображения по регистрации оптического излучения с учетом рассеяния. Решена модельная задача реконструкции изображения.

Научная новизна работы заключается в следующем.

1. Разработан оригинальный численный метод решения краевой задачи для уравнения лучистого переноса в различных средах, в том числе в биотканях.

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

3. В результате численного моделирования и расчетно-параметричес-ких исследований показана возможность нейрофизиологических исследований с помощью оптических томографических методов.

4. Сформулирован численный метод решения обратной задачи оптической томографии с учетом рассеяния излучения.

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

Основные положения, выносимые на защиту:

1. численный метод для решения прямой и обратной задачи;

2. результаты анализа теоретических и экспериментальных данных;

3. результаты расчетно-параметрических исследований.

Апробация работы.

Результаты диссертации представлялись на международной конференции "Применение лазеров в биологии и медицине" (Киев, октябрь 1995г.), международной конференции "International workshop on advanced electronics technology '95" (Москва, ноябрь-декабрь 1995г), XXIV Всероссийской школе-симпозиуме по когерентной оптике и голографии (Москва, январь 1996г), научной конференции ФАРШ МФТИ (Москва, май 1996г), XXXIX юбилейной научной конференции МФТИ "Современные проблемы фундаментальной и прикладной физики и математики" (Москва, ноябрь 1996г). По теме диссертации опубликовано 8 работ, список которых дан в конце автореферата.

Объем работы. Диссертационная работа изложена на 138 страницах машинописного текста, состоит из: Введения, 3-х глав, Заключения и Приложения, включает 56 рисунков и список основной использованной литературы из 78 наименований отечественных и зарубежных авторов на 7 страницах.

СОДЕРЖАНИЕ РАБОТЫ

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

В Главе 1 настоящей работы сделан обзор современного состояния томографических способов лучевой диагностики, приведены сведения об их принципиальных схемах, физических основах функционирования, областях применимости, проанализированы достоинства и недостатки. Пока-

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

В настоящее время можно выделить несколько основных способов медицинской томографии: рентгеновскую, позитронно-эмиссионную, одно-фотонную эмиссионную, ЯМР-томографию, ультразвуковую томографию. Каждый из этих способов, как правило, диагностирует свой круг заболеваний и можно выделить целый ряд органов и заболеваний которые недостаточно хорошо диагностируются этими методами. Идеальным способом диагностики мог бы стать безопасный недорогой метод с высоким пространственным и временным разрешением, позволяющий выполнять диагностические исследования в реальном масштабе времени. Широкое поле для исследований представляет использование в качестве детектируемого излучения — лазерного излучения в дальнем, ближнем ИК и видимом диапазонах.

Лазерное излучение уже сегодня используется для решения многочисленных задач медицинской диагностики и для терапевтических целей. Для определения внутренней структуры используется излучение видимого и инфракрасного диапазонов, которое не ионизирует ткани и достаточно хорошо переносится в больших дозах. По имеющимся сведениям, излучение этих диапазонов безопасно и по другим параметрам. Оптические источники недороги, легко транспортируются, могут использоваться длительное время. Спектр измеряемых параметров в случае оптического источника довольно широк.

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

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

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

Хотя универсальный оптический томограф для обследования всего организма еще не создан, есть примеры разработок приборов для диагностики отдельных органов и тканей, основанные на частотном методе. Так в апреле 1996 года в Далласе (США) на 27-ой конференции по раку груди фирмой Imaging Diagnostic Systems Inc. был продемонстрирован томографический лазерный маммограф, начались его клинические испытания. Годом раньше, в августе 1995 года фирма Humphrey Instruments (подразделение известной фирмы Carl Zeiss Inc.) начала внедрение и продажу оптического когерентного офтальмологического томографа.

В нашей стране исследования в области лазерной медицинской и биологической диагностики ведутся практически с момента создания оптических квантовых генераторов. Тем не менее, задача создания прибора для получения изображений внутренних органов человека с использованием оптических источников излучения в явном виде до последнего времени не ставилась. Широкое развитие получили другие методы диагностики, основанные на использовании рассеяния света (нефелометрия), дифракции, интерференции, спектроскопии, флуоресценции и т.д. С начала 1990-х годов интерес к оптической томографии растет. Исследуются возможности оптической диагностики многослойных структур, моделирующих органы человека, акустооптической томографии неоднородных сред, других методов оптической интроскопии.

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

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

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

В Главе 2 строится модель распространения излучения в сильнорас-сеивающих средах для плоскопараллельной геометрии. Для численной реализации этой модели разрабатывается специальный метод матриц переноса, с его помощью решается прямая задача распространения излучения, анализируются экспериментальные данные, численные результаты других авторов. Определяются основные свойства биотканей и их особенности. Определен диапазон длин волн, оптимальный для использования в оптическом томографе, коэффициенты отражения, поглощения, рассеяния, фазовые функции. Проводится решение задачи распространения излучения в многослойной ткани и расчетно-параметрическое исследование для моделирования оптической томографии.

Рассмотрим уравнение переноса:

При решении задачи о переносе излучения возникают две проблемы. Во-первых, это разработка эффективной численной схемы решения инте-гродифференциального уравнения переноса при наличии острых распределений индикатриссы и сильном рассеянии. Во-вторых, удовлетворение граничных условий, особенность которых заключается в том, что часть граничных условий задана на одной границе (для входящего излучения), а часть — на другой (для выходящего излучения). Оптимальным может считаться метод, решающий обе поставленные проблемы.

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

Используем особенности нашей задачи, состоящие в следующем. Во-первых, это — независимость коэффициентов задачи от интенсивности

<Л(г, т)

т-

йг

-1

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

Фоновым излучением пренебрегаем. Рассмотрим интенсивность падающего излучения 1(г,т) как совокупность компонент так, что направлению распространения каждой компоненты s¡ соответствует своя интенсивность Ц и направляющий косинус ть так что: 1(г,т)={11(г)=1(г,т1), 1=1Х...М

Уравнение переноса будет справедливо для каждой ¿-ой компоненты. Если теперь интеграл в правой части аппроксимировать суммой по компонентам, то получим:

Щ = -И/;«

Это уравнение представимо через матричные операторы:

¿г

Здесь матрица А — матрица переноса, она имеет следующие компоненты:

А = Ш = ~ 5,у)}

Рассмотрим уравнение переноса в матричном виде. Если компоненты матрицы переноса не зависят от пространственных координат в среде (среда однородная), то тогда матричное уравнение имеет решение, выражаемое через матричную экспоненту:

1(У)=ехр[А,г]1о

Реальная среда, рассматриваемая в данной работе существенно неоднородна. Введем пространственную дискретизацию. Положим, что направление т/ \г\ совпадает с направлением оси ОЪ. Разобьем путь луча на 5 отрезков, считая, что компоненты матрицы переноса А постоянны на каждом из этих отрезков. Для я-ого отрезка решение будет иметь вид:

13=ехр[А3^]13.1,

где к3 — длина 5-ого отрезка.

Из этой реккурентной формулы вытекает решение транспортного уравнения в приближении матриц переноса:

1 = Пехр[АЛ]1„

х

Полученная формула позволяет находить решение во всей интересующей области. Удовлетворение граничных условий сводится к одноразовому решению линейной алгебраической системы уравнений. Необходимость в любых итерациях отсутствует. Таким образом, можно считать, что построенная схема удовлетворяет сформулированным выше требованиям к оптимальной схеме для решения лучистого уравнения переноса.

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

*02

=£><3=

Граничные условия будут выглядеть следующим образом:

1о = Уо,}, 1=0,

I = п,п + к I

Здесь к компонент с п-ной по я+й-ую — это те компоненты, которые лежат в правой полуплоскости от обеих границ. Сначала вычислим матрицу переноса Р, которая в случае отсутствия фонового излучения равна:

Ф = Пехр[АД]

Затем матрица переноса и вектора интенсивностей представляются в блочном виде (с учетом найденных граничных условий):

Ф

Теперь можно записать следующее:

I = По, Р,10 = 0, Р210 = 12

Видно, что для того, чтобы решить матричное уравнение переноса нужно найти неизвестные компоненты 101 вектора 10, а затем найти 12.

ф2 , I = . I» = I»:

ф, 0 1о,

Сравнивая данные о — эффективном коэффициенте ослабления, для разных тканей и длин волн in vivo и in vitro, можно заметить, что: значения эффективного показателя ослабления ¿1ец не сильно отличаются для измерений in vivo и in vitro, хотя в случае in vitro наблюдается слабая тенденция к его увеличению (не более, чем в 1.5 — 2 раза); абсолютные значения не велики и лежат в пределах от 1.0 до 20.0 см*1, что говорит о небольшом диффузионном ослаблении, глубина проникновения излучения составляет 0.05 — 1.00 см; в диапазоне 600 — 1000 нм четко прослеживается наличие "окна прозрачности". Для биотканей характерной особенностью является то, что коэффициент рассеяния у них существенно больше коэффициента поглощения, причем отличие может достигать нескольких порядков. Так, например, у мышцы кролика в "окне прозрачности" поглощение составляет порядка 1см-1, в то время как рассеяние в том же диапазоне свыше 300 см-1. Для показателей поглощения и рассеяния печени крысы в "окне прозрачности" характерны значения 3.0 — 7.0 и 100.0 — 150.0 см-1 соответственно. В общем же если коэффициент поглощения биотканей лежит в пределах 0.5 — 10.0 см-1, то для коэффициента рассеяния характерны числа 10 — 500 см-1.

Вид индикатриссы заранее, как правило, не известен, поэтому при решении как прямой, так и обратной задачи приходится прибегать к различной ее аппроксимации, используя априорную информацию о свойствах исследуемых тканей. В данном случае такой информацией является значение параметра анизотропии g (g = 0.68 — 0.98), характерное значение которого равно 0.9. Близкое к единице значение g говорит о сильной вы-тянутости индикатриссы вперед.

Важная задача, возникающая при обработке полученных с помощью оптических источников изображений — это интерпретация полученной картины распределения интенсивности, идентификация наблюдаемых на ней объектов. В качестве примера решения такой задачи рассмотрим фотографии пальца человека, сделанные на длине волны 788 нм (рис. 1). На фотографии видны параллельные вертикальные размытые полосы. Для того, чтобы понять их происхождение, решим уравнение переноса излучения методом матриц переноса.

Сначала профотометрируем снимок пальца. Для того, чтобы при этом избежать ошибки в определении месторасположения темных полос и величины интенсивности в этих местах, усредним фотометрируемые интенсивности по вертикали в прямоугольнике, выделенном на рис. 1 светлыми линиями. Результаты фотометрирования представлены на рис. 2.

Далее, воспользуемся моделью сустава пальца, которая дает нам представление о расположении слоев ткани в его поперечном сечении (см. рис. 3). При этом положим для ткани кости = 1.5 см-1, - 180 см-1, g = 0.97; для кожи //а = 7 см-1, /а. = 130 см*1, g = 0.93; для синовиальной оболочки /1а = 2.4 см-1, = 200 см1, g = 0.95 и проведем расчеты интенсивности прошедшего через поперечное сечение пальца излучения методом матриц переноса. Получим кривую распределения интенсивности поперек пальца (см. рис. 2). Из сопоставления расчета с данными фотометрии видно, что темные полосы на фотографии образованы синовиальной оболочкой, суммарный коэффициент ослабления (¿ц = ц3 + ps) которой больше, чем у кости и кожи. Согласие расчета с экспериментом вполне удовлетворительное, если не принимать во внимание краевые эффекты.,

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

Нейрофизиология — деликатная область, малоизученная, труднодоступная с помощью ионизирующих излучений или инвазивного вмешательства. На основании обзора томографических методов можно утверждать, что оптическое излучение в ИК диапазоне — почти идеальное средство исследования головного мозга, при условии решения всех математических проблем обработки сигналов. В настоящей работе для иллюстрации возможностей применения разработанного метода матриц переноса рассматривается модельная задача прохождение излучения через многослойную ткань мозга. Целью численного расчета было определение чувствительности характеристик прошедшего излучения (его интенсивности) к изменениям состояния исследуемой ткани.

Целый ряд мозговых процессов (кровоизлияние, кровенаполненность, прохождение нервных импульсов и т.д.) сопровождается изменением оптических свойств определенных участков мозга. Для анализа возможности регистрации подобных изменений рассмотрим простейшую многослойную модель тканей мозга человека. В этой модели слоями располагались ткани кожи, кости, белого вещества, серого вещества, и снова белого вещества, кости и кожи. Были приняты следующие значения оптических свойств:

кожа — /¿j = 0.7 см1, - 130 см1, g = 0.93; кость — //а = 1.8 см-1, /4 = 160 см"1, g = 0.97; белое вещество — - 1.6 см1, = 51 см1, g = 0.96; серое вещество — д, = 2.6 см-1, /а. = 60 см-1, g = 0.88.

Для моделирования изменений состояния исследуемой ткани в ней размещали элемент, моделирующий сосуд с кровью, и изменяли его местоположение, кровенаполненность (концентрация воздуха изменялась при этом от.0 до 15% объема), размер (от 1 до 10% всей толщины многослойной ткани). При этом оптические свойства крови принимались следующими: //а = 2.84 см1, ^ = 505 см"1, g = 0.992.

Такая постановка задачи соответствует актуальной медицинской задаче ■— исследованию кровенаполненности кровеносных сосудов, расположенных внутри головного мозга человека. Адекватная диагностика кровенаполненности способствует повышению эффективности соответствующих хирургических операций и терапии.

Результаты расчетов представлены на рис. 5-7. Ввдно, что модель использованная для расчетов, позволяет учесть изменения в прошедшей интенсивности, связанные с положением сосуда, хотя согласно этой модели изменения и незначительны — доли процента. Гораздо большая чувствительность наблюдается к изменению кровенаполненности, при ее уменьшении на десять процентов интенсивность прошедшего излучения увеличивается на пять процентов. Зависимость интенсивности прошедшего излучения от размера сосуда носит экспоненциальный характер.

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

Глава 3 посвящена решению обратной томографической задачи. Обосновываются приближения, сделанные при решении. Метод матриц переноса обобщается для решения обратной задачи определения оптических свойств биоткани по данным сканирования. Осуществляется выбор алгоритма реконструкции, построение фантома, предлагается схема сканирования. Для разработанного фантома и предложенной схемы сканирования решается прямая задача прохождения излучения через ткань, а затем, на основании полученных проекционных данных, решается обратная задача восстановления изображения фантома. Обсуждается область применимости предложенного метода решения обратной задачи и дальнейшие этапы создания оптического томографа.

Аппроксимируя индикатриссу рассеяния можно получить два варианта упрощенного уравнения переноса, соответственно для транспортного приближения:

ёа1{Г'то) = -Ц,1(г,тй) + - е)1{г,т0)

аг

и приближения дельта-функцией:

= -м,1(г, т0) + мМ', Щ)

аг

Здесь - /а — альбедо единичного рассеяния.

Видно, что во втором и третьем случаях возможно представление уравнения переноса в простом виде:

Л{г,т0)

—-— = -Ы(г, т0 ) а!

Здесь Ь е — эффективный коэффициент ослабления, который равен соответственно для транспортного приближения и приближения дельта-функцией: ИеЦ1=Иа/МеЦ2=Ма/ё- Таким образом, возможно аналитическое решение уравнения переноса следующего вида:

1(г,т0) =/(г0,т0) ехр[- ¡леии2г\

Разобьем восстанавливаемое изображение на ячейки. Для определенности представим наше изображение в виде квадрата, разбитого на ячейки. Оптическая плотность изображения представится в виде вектора Ь, где Ь/ — оптическая плотность /-ой ячейки. Тогда основное уравнение компьютерной томографии можно представить в виде произведения матричных операторов:

Щ ■ Ь; = и,- или К • В = и

Здесь Е — это матрица сканирования — элемент равен расстоянию, которое г'-ая прямая сканирования проходит в /-ой ячейке изображения; К; — это обработанная информация о прошедшем излучении полученная от ¿-ой прямой сканирования.

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

Если матрица сканирования известна нам с точностью Л, а правая часть и с точностью 8, то данная задача может быть некорректной (корректность понимается в смысле Адамара), т.е. не выполняется условие непрерывности обратного к интегральному оператора (даже если точное решение этой задачи с приближенной правой частью и оператором и существует, то оно может быть неустойчиво к малым возмущениям исходных данных). Алгоритмы решения таких уравнений хорошо разработаны в работах А.Н.Тихонова и основаны на методе регуляризации.

Решение В, которое минимизирует регуляризующий функционал удовлетворяет уравнению Эйлера:

ВДйВ - аВ = О

Таким образом приходим к системе уравнений: АаВ = АВ + аСВ = Р

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

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

Выше уравнение переноса было упрощено с использованием плоскопараллельной геометрии и некоторых простых приближений индикатрис-сы. Для того, чтобы полнее учесть эффекты рассеяния, используем для решения обратной томографической задачи метод матриц переноса. Рассмотрим матрицу переноса А.

А = -^А, -Ц,А2 =Ц,Аз

здесь

= алМ'^) ^ =^А -А

л, л2 — л3 - л, л2

т. т, 2

Решение транспортного уравнения методом матриц переноса:

1 = Пехр[АД]1»

s

Рассмотрим матричную экспоненту:

exptAAJsS^-

/-о I ■

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

exp[As/zs] = Е + A shs

и, соответственно, для произведения матричных экспонент: П«р[ал] = Е+5>Л

s s

Решение транспортного уравнения будет выглядеть так:

i = i, +5>.М.

j

Теперь предположим, что путем введения некоторых "усредненных" значений альбедо и среднего косинуса рассеяния или подбором определенной аппроксимации индикатриссы нам удалось сформировать матрицу А3 не зависящую от пространственных координат. Тогда

As = щ, А3 и: I -10 = (Н-ВМо

Здесь компоненты вектора В — значения во всех ячейках восстанавливаемого изображения, а компоненты вектора H — расстояния, пройденные лучем сканирования в каждой такой ячейке. Заметим, что произведение векторов (НВ) есть некоторое число, не зависящее от номера компоненты вектора интенсивности I. Поэтому для каждой из этих компонент можно записать:

h - I0i = (H-B)-/w или (/,- - I0i)/ I0i = Н-В

Обозначим (/,■ - I0i)/ loi - Uh тогда C/(- = H-B

Пусть для решения томографической задачи реализуется схема сканирования по К направлениям, тогда для каждой прямой сканирования:

Ui k = %В или U; = R-B

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

Теперь предположим, что нам известно "усредненное" значение не среднего косинуса рассеяния и альбедо, как в прошлый раз, а только "усредненное" значение среднего косинуса. В этом случае сформировать независимую от пространственных координат матрицу А3 мы не можем, зато можем сформировать матрицы А] и А2.

Рассуждая похожим образом, в этом случае также можно получить уравнение переноса для каждой компоненты вектора интенсивности:

S S

или для компоненты i вектора интенсивности Uoi = Uu ВГН - U2i В2-Н

и в обезразмеренном виде:

(Uoi/Un) = BrH - (U2l/U,i)BrH = Н (Bi - Шц/U,,)В2)

Здесь U0i - Ii - ¡oi . Aiij ^v = A2ij h\. Обозначая для fe-ой

прямой сканирования

Uik = W0i/UH)k и Bf = Bi - (U2i/UH)B2 ,

окончательно получим

Ui k = Hft- Bi или U, = R- B(,

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

Bs = juts (1/2 w0s - U2i/Un)

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

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

Решим прямую задачу прохождения луча света через построенный фантом. Для фантома выберем поперечное сечение головного мозга человека (см. рис. 8). Возьмем небольшой фрагмент этого сечения, профото-метрируем, получим фантом для восстановления (см. рис. 9). Огтические свойства головного мозга для диапазона длин волн 0.8 — 0.9 мкм возьмем: для белого вещества д = 52.6 см-1, д = 0.96, для серого вещества /4 = 62.8 см"1, g = 0.88

Применим транспортное приближение индикатриссы. На рис. 4 показано, к чему приводит замена коэффициентов поглощения и рассеяния на коэффициент Ь в уравнении переноса. Видно, что качественная картина распределения прошедшего излучения для такого грубого приближения, которым по существу является транспортное, сохраняется достаточно хорошо и, следовательно, оно может быть использовано для оценочных расчетов.

Уровень погрешностей в операторе примем равным 10"10, а в правой части 10"8. Результаты решения прямой задачи будут служить исходными данными для обратной задачи, которую решим описанным выше матричным алгоритмом.

Результаты решения обратной задачи приведены на рис. 10. Сравнивая рис. 9 и рис. 10 можно сказать, что все основные объекты, которые мы видим на изображении фантома — "выпуклости" и "вогнутости" — хорошо просматриваются и на реконструированном изображении. Это говорит о том, что после реконструкции информация о расположении в объекте белого и серого вещества сохранилась, несмотря на небольшие различия (порядка 4-5%) в их оптических свойствах. Видно, что матричный алгоритм обеспечивает приемлемое качество восстановления исходного фантома, относительная ошибка при этом составляет около 2%.

ОСНОВНЫЕ ВЫВОДЫ И РЕЗУЛЬТАТЫ ДИССЕРТАЦИОННОЙ РАБОТЫ

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

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

Новизна работы заключается в создании метода матриц переноса для решения прямой и обратной задач переноса излучения, позволяющего: разработать на его основе эффективную численную схему решения инте-гродифференциального уравнения переноса при наличии острых распределений индикатриссы и сильном рассеянии; удовлетворить граничным условиям, особенность которых заключается в том, что часть граничных условий задана на одной границе (для входящего излучения), а часть — на другой (для выходящего излучения); учесть рассеяние при решении обратной томографической задачи. Проведено сравнение расчетов методом матриц переноса с имеющимися теоретическими результатами. Показана эффективность работы метода, позволяющего некоторым образом учитывать все особенности задачи. Проведен предварительный анализ возможности использования данного численного метода в нейрофизиологии.

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

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

ОСНОВНЫЕ ПУБЛИКАЦИИ ПО ТЕМЕ ДИССЕРТАЦИИ:

1. Верблюденко П.А. Оптические свойства биотканей в видимом и ближнем ИК поддиапазонах длин волн./Применение лазеров в биологии и медицине./ / Сборник научных докладов,-тезисов и методик по лазерной медицине. Часть 1. Материалы международной конференции. Киев, 11-14 октября 1995 года, с.23

2. Верблюденко П.А. Сон Э.Е. Один алгоритм решения обратной задачи прохождения лазерного излучения через биоткани./Фундаментальные исследования в области прикладной физики и математики в технических вузах России.// Сборник аннотаций и научных статей. Том 1. -М., Государственный Комитет РФ по,высшему образованию, 1995.

3. PA.Verbludenko, E.E.Son. About the challenge of creating the laser medical tomograph. Optical properties of biological tissues in IR range / International workshop on advanced electronics technology '95. - Presidium of Russian academy of Science, Moscow, Russia, p.43

4. PA.Verbludenko, E.E.Son. About the challenge of creating the laser medical tomograph. Method of solving the invers problem of laser light transmission through biological tissue. / International workshop on advanced electronics technology '95. - Presidium of Russian academy of Science, Moscow, Russia, p.44

5. ПА.Верблюденко. Моделирование оптического компьютерного медицинского томографа / / Современные проблемы фундаментальной и прикладной физики и математики. Тезисы докладов XXXIX юбилейной научной конференции МФТИ. Выпуск 3. Механика, прикладная физика, химическая физика, физика живых систем. 1996, с. 13.

6. ПА.Верблюденко. Численное моделирование процессов в лазерной оптической медицинской томографии: Препринт МФТИ №97/2 — М., 1997, 26 с.

7. ПА.Верблюденко, Э.Е.Сон. Восстановление изображений в лазерной медицинской томографии. Труды XXIV Школы-симпозиума по когерентной оптике и голографии. Долгопрудный, 1996.

8. ПА.Верблюденко. Томография. В уч. пособии А.А.Веденова. Физика XX век. — М.: Химия и Жизнь, 1997г., с. 310-320. (в печати)

РИСУНКИ

Рисунок 1. Картина прохождения света длиной волны 788 нм сквозь палец

человека.

Рисунок 2. Распределение интенсивности света в прямоугольнике на рис. 1. Квадраты — результат расчета.

СИНОВИАЛЬНАЯ ОБОЛОЧКА

Рисунок 3. Модель сустава пальца (поперечное сечение).

0.5 1 1.5 Расстояние, отн. ед.

---Расчет с рис. 2

а Без поглощения

▲ Без рассеяния

-х Транспортное приближение

Рисунок 4. Сравнение расчетов при различных условиях.

2

0.3024

й 0.3023

Ф 0.3022 я

о 0.3021 ь~ 0.302

0

§ 0.3019

ш

у 0.3018

1 0.3017 к 0.3016

0.3015

0.4 0.5 0.6

Положение сосуда, отн. ед.

Рисунок 5. Зависимость относительной интенсивности прошедшего излучения от положения сосуда.

Я 0.37

0%

15%

5% 10%

Кровенаполненность

Рисунок 6. Зависимость относительной интенсивности прошедшего излучения от кровенаполненности.

Размер сосуда

Рисунок 7. Зависимость относительной интенсивности прошедшего излучения от размера сосуда.

Рисунок 8. Поперечное сечение головного мозга человека и его небольшой

фрагмент

Рисунок 10. Восстановленный фантом.