«МЕТОДЫ И АЛГОРИТМЫ РЕШЕНИЯ ЗАДАЧИ СТРУКТУРНОГО СИНТЕЗА СИСТЕМЫ ИСТОЧНИКОВ И ДЕТЕКТОРОВ ЗОНДИРУЮЩЕГО ИЗЛУЧЕНИЯ ...»

Нижегородский государственный университет им. Н.И. Лобачевского

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

Горшков Антон Валерьевич

МЕТОДЫ И АЛГОРИТМЫ РЕШЕНИЯ ЗАДАЧИ СТРУКТУРНОГО СИНТЕЗА СИСТЕМЫ ИСТОЧНИКОВ И ДЕТЕКТОРОВ ЗОНДИРУЮЩЕГО ИЗЛУЧЕНИЯ

05.13.01 – «Системный анализ, управление и обработка информации

(в наук

е и промышленности)» (технические науки)

Диссертация на соискание ученой степени кандидата технических наук

Научный руководитель д.т.н. Гергель Виктор Павлович

Научный консультант к.ф.-м.н. Кириллин Михаил Юрьевич Нижний Новгород -

СОДЕРЖАНИЕ

ВВЕДЕНИЕ

ПРОБЛЕМА МОДЕЛИРОВАНИЯ В ЗАДАЧАХ ОПТИЧЕСКОЙ

1.

ДИФФУЗИОННОЙ СПЕКТРОСКОПИИ

1.1. МЕТОД ОПТИЧЕСКОЙ ДИФФУЗИОННОЙ СПЕКТРОСКОПИИ (ОДС)

1.2. ОБЩАЯ ПОСТАНОВКА ЗАДАЧ СТРУКТУРНОГО АНАЛИЗА И СИНТЕЗА СИСТЕМЫ

ИСТОЧНИКОВ И ДЕТЕКТОРОВ ЗОНДИРУЮЩЕГО ИЗЛУЧЕНИЯ

1.3. ПРИМЕНЕНИЕ МЕТОДОВ МОНТЕ-КАРЛО ДЛЯ РЕШЕНИЯ УРАВНЕНИЯ ТЕОРИИ

ПЕРЕНОСА ИЗЛУЧЕНИЯ

1.3.1. КРАТКИЙ АНАЛИЗ ПОДХОДОВ К РЕШЕНИЮ УРАВНЕНИЯ ТЕОРИИ ПЕРЕНОСА

ИЗЛУЧЕНИЯ.

1.3.2. МЕТОД МОНТЕ-КАРЛО В ЗАДАЧЕ РАСПРОСТРАНЕНИЯ СВЕТА.

1.3.3. ПРОБЛЕМА ВЫБОРА ОПТИЧЕСКИХ ПАРАМЕТРОВ МОДЕЛИРОВАНИЯ.

1.3.4. ПРОБЛЕМА ОПИСАНИЯ СЛОЖНОЙ ГЕОМЕТРИИ.

1.3.5. СПОСОБЫ ПОВЫШЕНИЯ ЭФФЕКТИВНОСТИ МЕТОДА МОНТЕ-КАРЛО.

1.3.6. РЕАЛИЗАЦИЯ МЕТОДА МОНТЕ-КАРЛО НА ВЫЧИСЛИТЕЛЯХ С ПАРАЛЛЕЛЬНОЙ

АРХИТЕКТУРОЙ.

1.3.7. МЕТОД МОНТЕ-КАРЛО В ЗАДАЧЕ ОДС.

1.4. ПОСТАНОВКА ЗАДАЧИ ДИССЕРТАЦИОННОГО ИССЛЕДОВАНИЯ

МОДИФИЦИРОВАННЫЙ МЕТОД МОНТЕ-КАРЛО ДЛЯ РЕШЕНИЯ ЗАДАЧ

2.

СТРУКТУРНОГО АНАЛИЗА И СИНТЕЗА СИСТЕМЫ ИСТОЧНИКОВ И

ДЕТЕКТОРОВ ЗОНДИРУЮЩЕГО ИЗЛУЧЕНИЯ

2.1. ПОСТАНОВКА ЗАДАЧИ

2.2. МОДИФИЦИРОВАННЫЙ МЕТОД МОНТЕ-КАРЛО

2.2.1. ОБЩЕЕ ОПИСАНИЕ

2.2.2. ИНИЦИАЛИЗАЦИЯ ФОТОНА.

2.2.3. ПОЛУЧЕНИЕ ДЛИНЫ СВОБОДНОГО ПРОБЕГА ФОТОНА.

2.2.4. ПЕРЕМЕЩЕНИЕ ФОТОНА.

2.2.5. ОБНОВЛЕНИЕ ТРАЕКТОРИИ ДВИЖЕНИЯ ФОТОНА.

2.2.6. ОБНОВЛЕНИЕ ВЕСА ФОТОНА.

2.2.7. ВЫЧИСЛЕНИЕ НАПРАВЛЕНИЯ ДВИЖЕНИЯ ФОТОНА.

2.2.8. ПРИМЕНЕНИЕ МЕТОДА СУЩЕСТВЕННОЙ ВЫБОРКИ.

2.2.9. ВЫБОР ВЕКТОРА ПРИТЯЖЕНИЯ.

2.2.10. ПЕРЕСЕЧЕНИЕ ФОТОНА С ГРАНИЦЕЙ СЛОЯ.

2.2.11. ЗАВЕРШЕНИЕ МОДЕЛИРОВАНИЯ ФОТОНА.

АВТОМАТИЗИРОВАННАЯ СИСТЕМА ПОДДЕРЖКИ СТРУКТУРНОГО

3.

АНАЛИЗА И СИНТЕЗА СИСТЕМЫ ИСТОЧНИКОВ И ДЕТЕКТОРОВ

ЗОНДИРУЮЩЕГО ИЗЛУЧЕНИЯ

3.1. ТРЕБОВАНИЯ К СОСТАВУ ОБЕСПЕЧЕНИЙ АВТОМАТИЗИРОВАННОЙ СИСТЕМЫ ................

3.2. ВЫСОКОУРОВНЕВАЯ АРХИТЕКТУРА ПРОГРАММНОГО ОБЕСПЕЧЕНИЯ

АВТОМАТИЗИРОВАННОЙ СИСТЕМЫ

3.3. ОБЩАЯ СХЕМА ФУНКЦИОНИРОВАНИЯ ПРОГРАММНОГО ОБЕСПЕЧЕНИЯ

3.4. МЕТОДИКА УПРАВЛЕНИЯ СИСТЕМОЙ ИСТОЧНИКОВ И ДЕТЕКТОРОВ С

ИСПОЛЬЗОВАНИЕМ АВТОМАТИЗИРОВАННОЙ СИСТЕМЫ

3.5. ОПТИМАЛЬНЫЙ АЛГОРИТМ ПОИСКА ПЕРЕСЕЧЕНИЙ

3.5.1. ПОСТАНОВКА ЗАДАЧИ.

3.5.2. СУЩЕСТВУЮЩИЕ АЛГОРИТМЫ ПОИСКА ПЕРЕСЕЧЕНИЙ.

3.5.3. ВЫБОР ОПТИМАЛЬНОГО АЛГОРИТМА.

3.6. ПАРАЛЛЕЛЬНЫЕ АЛГОРИТМЫ ДЛЯ РАЗЛИЧНЫХ АРХИТЕКТУР

3.6.1. ПАРАЛЛЕЛЬНЫЙ АЛГОРИТМ ДЛЯ МНОГОЯДЕРНЫХ ЦЕНТРАЛЬНЫХ ПРОЦЕССОРОВ...... 3.6.2. ПАРАЛЛЕЛЬНЫЙ АЛГОРИТМ ДЛЯ СОПРОЦЕССОРОВ INTEL XEON PHI

3.6.3. ПАРАЛЛЕЛЬНЫЙ АЛГОРИТМ ДЛЯ ГРАФИЧЕСКИХ ПРОЦЕССОРОВ.

3.6.4. ПАРАЛЛЕЛЬНЫЙ АЛГОРИТМ ДЛЯ СИСТЕМ С РАСПРЕДЕЛЕННОЙ ПАМЯТЬЮ................. АНАЛИЗ И ТЕСТИРОВАНИЕ МОДИФИЦИРОВАННОГО МЕТОДА МОНТЕКАРЛО

4.1. ПОДБОР ПАРАМЕТРОВ ПРЕДЛОЖЕННОГО МЕТОДА

4.2. ПРОВЕРКА КОРРЕКТНОСТИ ПРЕДЛОЖЕННОГО МЕТОДА

4.2.1. СРАВНЕНИЕ С РЕЗУЛЬТАТАМИ ДИФФУЗИОННОЙ ТЕОРИИ.

4.2.2. СРАВНЕНИЕ С РЕЗУЛЬТАТАМИ СТАНДАРТНОГО МЕТОДА МОНТЕ-КАРЛО.

4.3. АНАЛИЗ ВЛИЯНИЯ ПАРАМЕТРОВ МЕТОДА НА ТОЧНОСТЬ РАСЧЕТОВ

4.3.1. АНАЛИЗ ВЛИЯНИЯ КОЛИЧЕСТВА МОДЕЛИРУЕМЫХ ФОТОНОВ НА ТОЧНОСТЬ

РАСЧЕТОВ.

4.3.2. АНАЛИЗ ВЛИЯНИЯ МИНИМАЛЬНОГО ВЕСА ФОТОНА НА ТОЧНОСТЬ РАСЧЕТОВ............

4.4. АНАЛИЗ УСТОЙЧИВОСТИ МЕТОДА К ИЗМЕНЕНИЮ ВХОДНЫХ ПАРАМЕТРОВ

АЛГОРИТМА

4.4.1. АНАЛИЗ УСТОЙЧИВОСТИ МЕТОДА ПРИ ИСПОЛЬЗОВАНИИ РАЗНЫХ

ПОСЛЕДОВАТЕЛЬНОСТЕЙ ПСЕВДОСЛУЧАЙНЫХ ЧИСЕЛ.

4.4.2. АНАЛИЗ УСТОЙЧИВОСТИ МЕТОДА ПРИ ИЗМЕНЕНИИ ОПТИЧЕСКИХ ХАРАКТЕРИСТИК

БИОТКАНИ

4.5. МОДЕЛИРОВАНИЕ ОПТИЧЕСКОЙ ФУНКЦИОНАЛЬНОЙ ДИАГНОСТИКИ ГОЛОВНОГО

МОЗГА

4.5.1. АНАЛИЗ ЗАВИСИМОСТИ СИГНАЛА ОТ РАССТОЯНИЯ «ИСТОЧНИК ДЕТЕКТОР» В ТРЕХ

ГЕОМЕТРИЯХ.

4.5.2. ИССЛЕДОВАНИЕ ГЛУБИНЫ ЗОНДИРОВАНИЯ С ПОМОЩЬЮ КАРТ РАССЕЯНИЯ.............. 4.5.3. АНАЛИЗ РАСПРЕДЕЛЕНИЯ ФОТОНОВ ПО ДЛИНАМ ПРОБЕГА.

4.5.4. СРАВНЕНИЕ РЕЗУЛЬТАТОВ МОДЕЛИРОВАНИЯ С МОДЕЛЬНЫМ ЭКСПЕРИМЕНТОМ........ ЗАКЛЮЧЕНИЕ

СПИСОК ЛИТЕРАТУРЫ

ВВЕДЕНИЕ

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

Ряд методов оптической биомедицинской диагностики уже прошел стадию предклинических испытаний и успешно применяется в клинической практике. Одним из таких методов является оптическая когерентная томография (ОКТ), обеспечивающая имиджинг внутренней структуры биотканей с высоким пространственным разрешением (до единиц микрон) [1, 2]. Однако, глубина видения для данного метода существенно ограничена многократным рассеянием зондирующего излучения в биотканях (так, для кожи человека она составляет не более 2 мм), поэтому основные области применения ОКТ связаны с офтальмологией [3] и диагностикой поверхностных тканей [4].

Для диагностики биотканей на больших глубинах необходимо применять методы, для которых информативным является многократно рассеянное излучение. Одним из таких методов является оптическая диффузионная спектроскопия (ОДС), предоставляющая широкие возможности для неинвазивной диагностики [5, 6]. Метод основан на регистрации многократно рассеянного объектом зондирующего излучения на нескольких длинах волн, определяемых спектрами поглощения исследуемых компонент организма.

Применение метода оптической диффузионной спектроскопии позволяет решать такие задачи, как диагностика раковых опухолей [7, 8], в частности, рака груди; мониторинг активности зон коры головного мозга [9]; планирование фотодинамической терапии [10, 11]; мониторинг состояния пациента при хирургическом вмешательстве [12]; определение состояния кожных покровов [13]; и др.

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

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

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

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

Объектом исследования в данной работе является система источников и детекторов, взаимодействующих посредством переноса зондирующего излучения в сложных гетерогенных средах, а предметом – методы и алгоритмы Монте-Карло моделирования для решения задач структурного анализа и синтеза такого рода систем.

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

В соответствии с целью диссертационной работы поставлены следующие задачи:

1. Разработка более эффективного с точки зрения быстродействия модифицированного метода Монте-Карло для решения задач структурного анализа и синтеза системы источников и детекторов зондирующего излучения;

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

3. Исследование эффективности разработанного метода и параллельных алгоритмов на примере решения задачи моделирования ОДС;

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

5. Применение полученного программного комплекса для моделирования функциональной диагностики головного мозга человека методом ОДС.

Методы исследования.

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

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

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

2. Предложены новые параллельные алгоритмы обработки информации о траекториях движения фотонов в сложных гетерогенных средах, предназначенные для исполнения на многоядерных процессорах, графических процессорах, ускорителях Intel Xeon Phi и гетерогенных кластерных системах.

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

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

Структура и объем работы. Диссертация состоит из введения, четырех глав, заключения и списка литературы. Общий объем работы составляет 137 страниц, включая 70 рисунков и 8 таблиц. Список литературы включает 93 наименования.

Содержание работы.

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

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

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

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

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

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

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

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

В §1.3 дано описание методов Монте-Карло для решения уравнения теории переноса излучения, а так же краткий обзор других подходов к решению данной задачи.

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

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

Еще одним широко распространенным классом методов, которые с успехом применяются для решения уравнения переноса, являются сеточные методы, и в частности метод конечных разностей. Наиболее известный в данной области программный пакет – NIRFAST (Near Infrared Fluorescence and Spectral Tomography), подходит для решения целого ряда задач реконструкции в области оптики биотканей, в частности, задач обнаружения и локализации раковых опухолей.

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

Методами, которые лишены описанных выше недостатков, являются методы МонтеКарло для моделирования распространения зондирующего излучения. Однако при использовании этих методов возникает другая сложность – для проведения численного моделирования требуется обработка значительного объема информации (необходимо моделировать траектории движения 108 - 1010 пакетов частиц), а значит и существенные вычислительные ресурсы.

В §1.4 приведена постановка цели диссертационного исследования: создание инструментариев поддержки структурного анализа и синтеза системы источников и детекторов, взаимодействующих посредством переноса зондирующего излучения в сложных гетерогенных средах.

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

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

В работе представлена обобщенная модель системы источников и детекторов зондирующего излучения. В основе лежит модель переноса фотонов в плоскопараллельных гетерогенных биологических тканях, применяемая в программном пакете MCML (Monte Carlo Modeling of Light transport).

Исследуемый объект рассматривается в трехмерном пространстве и состоит из набора слоев. Каждый слой описывает определенный тип биологической ткани, обладающей набором оптических характеристик. Оптические характеристики постоянны в рамках слоя. Например, если моделировать перенос излучения в голове человека, то в ней можно выделить такие слои, как кожа головы, жировая ткань, череп, цереброспинальная жидкость, серое и белое вещество головного мозга.

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

Источник излучения представляет собой бесконечно тонкий луч фотонов и описывается направлением и положением в трехмерном пространстве.

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

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

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

Идея метода Монте-Карло в данной задаче состоит в случайной трассировке набора фотонов в среде. Фотоны объединяются в пакеты, каждый пакет обладает весом. Далее понятия «фотон» и «пакет фотонов» будут отождествляться. Начинает движение пакет фотонов от источника излучения. Далее на каждом шаге трассировки случайным образом определяется его направление, величина смещения и поглощенный вес. Моделирование пакета завершается либо при его поглощении средой (когда вес пакета становится меньше минимального), либо если он вылетает за границы исследуемого объекта.

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

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

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

Данная идея была известна ранее, однако разработанные на ее основе методы поддерживали работу только с однородным плоскопараллельным и полубесконечным слоем.

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

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

В §3.1 приведены требования к функциональности АС и ее составу как системы обеспечений (технического, информационного, математического, алгоритмического, программного и методического).

В §3.2 описана высокоуровневая архитектура программного обеспечения АС. Разработанный программный комплекс содержит модуль моделирования распространения зондирующего излучения применительно к задачам анализа системы источников и детекторов, модуль для генерации трехмерных триангулированных поверхностей, модуль для трехмерной визуализации триангулированных поверхностей и модуль для визуализации результатов моделирования. Для создания программного кода использовались языки программирования C++ и C#, общий объем кода – 12 750 строк.

В §3.3 описана общая схема функционирования программного обеспечения АС.

В §3.4 приведена методика управления системой источников и детекторов зондирующего излучения с применением АС.

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

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

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

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

В §4.1 приведен анализ параметров модифицированного метода Монте-Карло, а также приведены рекомендации по выбору их конкретных значений.

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

Основное достоинство модифицированного метода Монте-Карло состоит в том, что для получения достаточного объема статистики для целевого детектора требуется выполнять трассировку значительно меньшего числа фотонов, чем при использовании стандартного метода Монте-Карло. Время моделирования при этом уменьшается на порядок.

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

В §4.4 приведен анализ устойчивости метода к изменению последовательности псевдослучайных чисел и оптических параметров среды.

В §4.5 приведен анализ результатов моделирования оптической функциональной диагностики головного мозга человека с использованием разработанной АС.

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

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

В заключении изложены основные научные и практические результаты диссертационной работы.

Апробация результатов и публикации.

Материалы диссертации опубликованы в 12 работах. Из них 5 статей в рецензируемых журналах («Вестник ННГУ» [14, 15], «Journal of Computational Science» [16], «Journal of Selected Topics in Quantum Electronics» [17], « Procedia Computer Science» [18]), из которых 4 входит в Перечень ВАК российских рецензируемых журналов, 6 работ, представляющие собой публикации в трудах конференций, 1 работа в сборнике «Суперкомпьютерное моделирование в науке, образовании и промышленности». Получено 1 свидетельство о государственной регистрации программы для ЭВМ №2013611386 «Моделирование распространения излучения в биотканях для задач оптической диффузионной спектроскопии головного мозга человека» от 9 января 2013 г.

Результаты диссертации обсуждались на семинарах кафедры Математического обеспечения ЭВМ факультета Вычислительной математики и кибернетики ННГУ и докладывались на следующих конференциях:

Международная конференция «Высокопроизводительные параллельные вычисления на кластерных системах» (Пермь, 2010);

Всероссийские конференции «Высокопроизводительные параллельные вычисления на кластерных системах» (Нижний Новгород, 2011 – 2013);

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

Международная конференция «High Performance Computing and Simulation»

(Амстердам, Нидерланды, 2012);

Международная конференция «Topical problems of Biophotonics» (Нижний

1. ПРОБЛЕМА МОДЕЛИРОВАНИЯ В ЗАДАЧАХ

ОПТИЧЕСКОЙ ДИФФУЗИОННОЙ СПЕКТРОСКОПИИ

1.1. Метод оптической диффузионной спектроскопии (ОДС) В медицинской практике часто возникают задачи определения биомедицинского состояния различных типов биологической ткани [19]. Это, например, задача диагностики и лечения раковых опухолей, в частности, рака груди [7, 8]; задача мониторинга степени оксигенации зоны коры головного мозга для определения ее активности [9]; планирование фотодинамической терапии [10, 11]; мониторинг состояния пациента при хирургическом вмешательстве [12]; определение состояния кожных покровов [13]; и др. Одним из подходов, который уже доказал свою эффективность при решении данного класса задач, является оптическая диффузионная спектроскопия.

Оптическая диффузионная спектроскопия (ОДС) – это метод, предназначенный для определения биомедицинского состояния биологической ткани посредством измерения ее параметров поглощения и рассеяния при использовании различных длин волн зондирования [5, 6].

Спектроскопия поглощения – подход, обычно применяемый в лабораторных условиях для измерения концентраций и молекулярного состояния абсорбирующих молекул. А метод ОДС позволяет применять спектроскопию поглощения для исследования биологических тканей человека. Основное различие этих методов состоит в том, что в случае применения обычной спектроскопии поглощения исследуемый объект должен быть нерассеивающим, а так же находиться в стеклянной или кварцевой кювете для анализа. В случае же ОДС, кюветой может являться тело человека, что позволяет выполнять биохимический анализ “in vivo”.

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

Одним из возможных применений метода ОДС на практике является функциональная диагностика головного мозга человека. Основная задача функциональной диагностики мозга – выявление активности той или иной зоны коры головного мозга человека “in vivo”.

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

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

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

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

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

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

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

1.2. Общая постановка задач структурного анализа и синтеза системы источников и детекторов зондирующего излучения Рассмотрим задачу анализа системы источников и детекторов зондирующего излучения подробнее.

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

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

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

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

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

Целью моделирования является возможность оценки следующих величин, результатов анализа, которыми являются:

Объемное распределение поглощенной интенсивности излучения в среде;

Интенсивность рассеянного назад излучения на детекторах;

Фотонные карты траекторий для каждого детектора.

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

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

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

1.3. Применение методов Монте-Карло для решения уравнения теории переноса излучения 1.3.1. Краткий анализ подходов к решению уравнения теории переноса излучения.

1. Аналитическое решение уравнения переноса излучения.

Законы распространения излучения в различных средах описываются в теории переноса излучения [5, 6]. В классической теории переноса излучения основным понятием является лучевая интенсивность. Она определяет средний поток энергии через площадку, сосредоточенный в телесном угле вблизи направления в интервале частот :

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

Основное дифференциальное уравнение этой теории называется уравнением переноса излучения. Стационарное уравнение теории переноса излучения для монохроматического света имеет вид:

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

Здесь предполагается, что внутри среды источников излучения нет.

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

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

Так, в случае слабого рассеяния, когда рассеивающая среда является разреженной, а рассеивающий объем не слишком велик, применяется приближение первого порядка теории переноса излучения. В этом случае решение уравнения переноса излучения можно получить итерационным путем [21]. При этом предполагается, что полная интенсивность, падающая на частицы, приближенно равна ослабленной падающей интенсивности, которая известна. Так называемое решение в приближении первого порядка справедливо для оптически тонких и слабо рассеивающих сред, когда интенсивность прошедшей волны (когерентной составляющей), описывается законом Бугера. А в случае очень узкого падающего пучка (например, луча лазера), приближение первого порядка справедливо и для более плотных тканей.

С ростом оптического пути диаграмма направленности рассеянного излучения становится почти изотропной. Для описания такого излучения можно использовать диффузионное приближение, которое позволяет значительно упростить исходное уравнение переноса излучения и свести его к решению дифференциальных уравнений [5, 6]. Это позволяет эффективно использовать диффузионное приближение как способ учета многократного рассеяния.

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

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

Во многих задачах лазерной диагностики биологических систем падающую волну можно представить в виде узкого коллимированного пучка. Это справедливо, например, для определения концентрации бактерий в жидкости при помощи лазерного луча [24], или для определения насыщения крови кислородом при помощи волоконно-оптического катетера [25]. В случае крупномасштабных рассеивающих неоднородностей, размеры которых велики по сравнению с длиной волны излучения, индикатриса рассеяния отлична от нуля в узком секторе углов в направлении распространения вперед. Это позволяет существенно упростить описание остронаправленных волновых пучков путем перехода к малоугловому приближению [26]. Уравнение переноса в малоугловом приближении допускает точное решение, однако оно выглядит довольно громоздко и оказывается не слишком удобным для численных расчетов. Более простое приближение можно получить, перейдя от малоуглового приближения к приближению, соответствующему диффузии рассеянного излучения по угловым переменным. Малоугловое приближение пригодно для не слишком больших трасс распространения, на которых пучок остается достаточно сильно сколлимированным.

Часто интерес представляет не сама лучевая интенсивность, а лишь интегралы от, дающие энергетические характеристики поля излучения. Если освещение диффузно, и среда достаточно замутнена, то экспериментальные результаты хорошо описывает двухпотоковая теория Кубелки и Мунка. Эта теория основана на модели двух световых потоков, распространяющихся в прямом и обратном направлениях. Однако данная теория неприменима для описания падающего на среду коллимированного пучка. В этом случае следует использовать четырехпотоковую теорию, которая учитывает в дополнение к диффузионным потокам два коллимированных лазерных пучка, внешний падающий и отраженный от задней поверхности образца [21, 27]. Семипотоковая теория является простейшим трехмерным представлением падающего лазерного пучка и рассеяния излучения в полубесконечной среде [28].

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

2. Сеточные методы.

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

Существует два основных подхода к решению уравнения переноса в плоскопараллельных средах методом конечных разностей [29].

Первый основан на решении дифференциальной формы первого порядка уравнения переноса. Этот метод известен как теория дискретного пространства переноса излучения (Discrete Space theory of Radiative Transfer, DS) [30, 31]. Он определяет алгоритм вычисления внутренних и исходящих полей излучения. Разностная схема метода имеет первый порядок сходимости. Данная схема может применяться для решения уравнения переноса в случае сферической симметрии [32]. Дополнительно она может использоваться для решения линейных задач в расширяющейся среде в неподвижной и в совместно движущейся системе координат [33, 34], а также для решения задачи поляризованного поля излучения [35, 36].

Второй подход опирается на решение дифференциальной формы второго порядка для уравнения переноса излучения. В частности, в работе [37] предложена центральная разностная схема, дающая второй порядок сходимости. Далее она была улучшена до четвертого порядка сходимости путем использования Эрмитовой схемы [38]. Данный метод широко применяется в различных приложениях, например, для моделирования атмосферы [39]. Схема второго порядка сходимости также используется для решения задач, например, задачи переноса поляризованного излучения.

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

Наиболее известным пакетом решения уравнения теории переноса излучения, который основан на применении сеточных методов, является пакет NIRFAST (Near Infrared Fluorescence and Spectral Tomography) [23, 40, 41]. Он предназначен для моделирования распространения света ближнего инфракрасного диапазона в неоднородных биологических тканях и включает в себя решение как прямой задачи моделирования распространения излучения, так и задачи восстановления оптических параметров ткани.

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

Программный пакет NIRFAST хорошо подходит для решения целого ряда задач реконструкции в области оптики биотканей, в частности, задач обнаружения и локализации раковых опухолей [42].

Следует отметить, что помимо решения прямой и обратной задачи моделирования распространения излучения, NIRFAST может решать проблему получения реальной геометрии биологических тканей на основании результатов компьютерной (КТ) и магнитнорезонансной (МРТ) томографии [41].

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

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

1.3.2. Метод Монте-Карло в задаче распространения света.

Применение подхода на основе метода Монте-Карло для решения задачи моделирования распространения света – еще один способ численного решения уравнения переноса излучения. Использование метода Монте-Карло в данной задаче основано на вероятностной природе распространения излучения в средах.

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

Идея метода состоит в трассировке набора фотонов в среде. Для улучшения качества метода рассматриваются не отдельные фотоны, а так называемые пакеты фотонов. Каждому пакету ставится в соответствие определенный начальный вес. Величина веса в начальный момент времени должна определяться, исходя из мощности источника излучения. Обычно, для упрощения расчетов, начальный вес принимают за единицу. Далее понятия «фотона» и «пакета фотонов» будут отождествляться.

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

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

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

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

Дополнительными результатами моделирования в работе [44] являются:

интенсивность рассеянного назад излучения на верхней границе слоя;

интенсивность рассеянного вперед излучения на нижней границе слоя.

Таким образом, все фотоны можно разделить на три типа (рисунок 1.4):

1. Фотоны, поглощенные внутри слоя;

2. Диффузионно-отраженные фотоны, вылетевшие за пределы слоя сверху;

3. Прошедшие фотоны, вылетевшие за пределы слоя снизу.

Следует отметить и процедуру завершения процесса трассировки. В работе [44] моделирование пакета фотонов завершается тогда, когда вес пакета становится меньше определенной достаточно малой величины (например, 0.0001). При этом применяется подход, который носит название «рулетка». Он состоит в том, что при достижении достаточно малого веса, пакету дается шанс 1 из продолжить движение в слое с весом.

Дальнейшее развитие метода Монте-Карло в задаче моделирования распространения излучения велось в направлении обобщения этого метода на неоднородные среды. Так в работе [45] описан алгоритм трассировки пакета фотонов в многослойной среде, где каждый слой обладает своим набором оптических характеристик. Каждый слой бесконечен по двух направлениям и имеет только ширину. Соответственно граница слоя представляет собой две параллельные плоскости.

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

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

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

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

Отметим, что описываемые здесь подходы реализованы в рамках программного пакета для моделирования распространения излучения в гетерогенных средах MCML (Monte Carlo Modeling of Light transport) [45], который является наиболее известной реализацией метода Монте-Карло в рассматриваемой задаче. Алгоритм, реализованный в пакете, является базой для построения более эффективных либо более общих методов моделирования распространения света другими исследователями.

1.3.3. Проблема выбора оптических параметров моделирования.

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

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

Но даже если принять описанную выше оптическую модель биоткани как адекватную для рассматриваемого класса задач, остается актуальной проблема выбора конкретных оптических характеристик для каждого из слоев ткани. Иными словами, определение оптических характеристик для каждого конкретного типа биоткани – это отдельная довольно сложная задача. И на текущий момент она не решена полностью даже для набора хорошо известных типов тканей. Есть лишь ряд приближенных значений, полученных путем проведения реальных физических экспериментов над исследуемым типом ткани и последующего моделирования распространения света с целью решения обратной задачи для определения оптических характеристик среды [46, 47]. Причем в работах разных исследователей обычно даются разные (хотя и обычно близкие) значения этих характеристик для одного и того же типа биологической ткани.

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

1.3.4. Проблема описания сложной геометрии.

Одним из направлений исследований в области Монте-Карло моделирования распространения излучения в различных средах являются способы описания геометрии сложных объектов. В первых работах, посвященных описанию метода Монте-Карло для решения уравнения переноса излучения, границы объектов и их слоев описывались с помощью плоскостей [43, 44, 45]. Основное достоинство плоскопараллельной геометрии состоит в легкости определения факта пересечения фотона с границей. Основной недостаток – плохая аппроксимация реальных объектов. Например, в задаче моделирования распространения света в голове человека геометрия различных тканей головы далека от плоскопараллельной. А значит, требуется развитие существующих методов с тем, чтобы они имели возможность моделирования объектов с достаточно сложной геометрией.

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

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

1. Аналитически заданные поверхности.

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

Такой подход используется, например, в работе [48], где производится сравнение метода Монте-Карло с сеточным методом конечных элементов. В качестве тестового объекта там рассматривается объем, состоящий из двух вложенных друг в друга цилиндров разных диаметров.

Другой пример – аппроксимации геометрии головы человека набором плоскостей [49]. Внешний вид используемой в работе геометрии приведен на рисунке 1.5.

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

Рисунок 1.5. Аппроксимации геометрии головного мозга человека Достоинствами такого подхода являются простота и скорость поиска пересечений фотона с границами слоя и малый объем памяти, требуемый для хранения информации о геометрии объекта. Как результат, такой подход представления границ обеспечивает высокую скорость работы метода. Отметим, что перечисленные достоинства имеют место только в случае, когда граница слоя состоит из небольшого числа поверхностей. Это условие определяет и основной недостаток такого подхода: недостаточную степень аппроксимации реальной геометрии объекта при использовании небольшого числа простых поверхностей.

2. Воксельный подход.

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

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

С использованием данного подхода был разработан программный пакет для моделирования распространения излучения в биологических тканях с произвольной геометрией слоев [52]. Данный пакет был протестирован на корректность путем сравнения результатов его работы с диффузионной теорией, в результате чего была установлена применимость воксельного подхода на практике.

Рисунок 1.6. Воксельная геометрия головы человека, полученная из данных МРТ Впоследствии программный код на основе данного подхода был применен к задаче моделирования оптической диффузионной спектроскопии для функциональной диагностики мозга [53]. При этом геометрия головы человека, используемая для моделирования, была получена из данных МРТ.

Отметим, что процесс формирования воксельной геометрии исследуемого объекта сам по себе является сложной вычислительной задачей. На текущий момент предложены алгоритмы построения такой геометрии по данным МРТ и КТ. Например, в работе [54] приводится как алгоритм построения воксельной модели головы человека по данным МРТ и КТ, так и метод моделирования распространения излучения в полученной модели. Тем самым реализуется возможность проведения моделирования функциональной диагностики мозга с учетом физиологических особенностей каждого конкретного пациента.

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

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

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

3. Триангулированные поверхности.

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

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

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

Рисунок 1.8. Алгоритм поиска пересечений с помощью проекции [56] Так, в работе [56] поиск пересечения состоит из трех этапов (рисунок 1.8). На первом этапе поверхность и вектор направления движения фотона проецируются на двумерную плоскость. На втором этапе плоскость покрывается равномерной сеткой, после чего ищется пересечение двумерного вектора направления фотона с ячейками этой сетки. Это позволяет выделить те области, в рамках которых пересечение может иметь место. Третий этап включает в себя тест на пересечение прямой с треугольниками, попавшими в найденные ранее области. Данный алгоритм преследует целью сокращение числа тестов пересечения прямой с треугольником за счет отсечения заведомо лишних тестов.

Другой способ сокращения числа тестов на пересечение состоит в рекурсивном разбиении пространства на подобласти [57]. В результате такого разбиения получается октодерево, каждый узел которого описывает некоторый параллелепипед в исследуемой области и содержит либо набор треугольников, если это листовой узел, либо ссылки на узлыпотомки (рисунок 1.9). Это позволяет осуществлять поиск пересечения фотона с поверхностью следующим образом. На первом этапе ищется пересечение фотона с узлами дерева (соответствующими им параллелепипедами). Если найдено пересечение с некоторым узлом, поиск продолжается для узлов-потомков. На втором этапе, когда найден листовой узел дерева, с которым пересекается фотон, ищется пересечение со всеми треугольниками поверхности, которые содержит узел.

Рисунок 1.9. Алгоритм поиска пересечения путем разбиение двумерного Отметим, что в обоих описанных выше алгоритмах делается попытка сокращения числа тестов на пересечение фотона с треугольником. При этом в качестве самого теста используется алгоритм, описанный в работе [58].

1.3.5. Способы повышения эффективности метода Монте-Карло.

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

В данном пункте обсуждаются различные способы повышения эффективности метода Монте-Карло в задаче моделирования распространения света.

Первая группа методов основана на применении метода понижения дисперсии, рассматриваемого в рамках теории методов Монте-Карло. Идея данного метода состоит в увеличении точности алгоритма при сохранении постоянного объема статистических данных.

1. Схема расчетов с использованием весов.

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

Описанный подход является одним из способов понижения дисперсии и носит название схемы расчетов с использованием весов (в англоязычной литературе «survival biasing» или «implicit capture»). Применение данного подхода описано, например, в работах [43, 44, 45].

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

Другой подход – случайный обрыв траектории (в англоязычной литературе «Russian roulette technique»). При достижении определенного малого веса с вероятностью пакет фотонов продолжает движение с весом либо с вероятностью его трассировка прекращается. Подход используется в тех случаях, когда вес пакета фотонов становится достаточно малым и продолжать рассчитывать его траекторию невыгодно, так как его вклад в результат становится невелик, а вычислений много; и в то же время пренебречь таким весом нельзя [59]. Описанный метод завершения процесса трассировки используется, например, в работах [44, 45].

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

2. Метод существенной выборки.

Метод существенной выборки (в англоязычной литературе «importance sampling») – один из подходов к уменьшению дисперсии случайной величины при работе с методами Монте-Карло. Идея подхода основана на наблюдении, что некоторые значения случайной величины в процессе моделирования имеют большую значимость для оцениваемой функции, чем другие. Если эти «более значимые» значения будут появляться в процессе выбора случайной величины чаще, дисперсия оцениваемой функции уменьшится. Следовательно, базовая методология данного метода заключается в выборе распределения, которое способствует появлению «более значимых» значений случайной величины. Такое «смещенное» распределение изменяет оцениваемую функцию, если применяется прямо в процессе расчета. Однако, если результат расчета перевзвешивается согласно этому смещенному распределению, то гарантируется, что новая оцениваемая функция не будет смещена [59].

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

Так, в работе [60] рассматривается задача моделирования сигнала оптической когерентной томографии в однородных биологических тканях. При этом важным моментом там является разделение сигнала оптической когерентной томографии (ОКТ) на два типа.

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

Для более точного определения сигнала ОКТ первого типа в работе применяется метод смещенного угла (в англоязычной литературе «biased angle») [61]. Идея применения данного подхода состоит в следующем. При достижении фотоном нужной глубины угол его рассеяния меняется на противоположный, вес пакета фотонов при этом корректируется для получения несмещенной оценки искомой случайной величины. Такое вмешательство приводит к развороту фотона в сторону детектора, тем самым повышая вероятность учета его излучения (рисунок 1.10).

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

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

Рисунок 1.10. Метод смещенного угла в задаче определения сигнала ОКТ [60] Описанный в работе [62] алгоритм моделирования предназначен для работы с однородными плоскопараллельными средами и позволяет получить те же по точности результаты, что и стандартный метод Монте-Карло, но на порядок быстрее. Более того, данный подход представляется применимым и для задач оптической диффузионной спектроскопии.

Основной недостаток предлагаемого метода состоит в том, что процесс трассировки каждого пакета фотонов опирается на результаты моделирования предыдущего пакета, тем самым теряется важное свойство, обычно присущее алгоритмам на базе метода Монте-Карло – независимость фотонов друг от друга. А значит, выполнить распараллеливание данного метода представляется затруднительным.

В работе [63] решается задача вычисления коэффициента пропускания излучения через толстый однородный слой. Геометрия слоя является плоскопараллельной. Источник излучения расположен на верхней границе слоя и направлен перпендикулярно его границе. Детектор имеет ограниченные размеры и располагается на противоположной границе слоя соосно источнику излучения.

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

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

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

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

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

Таким образом, алгоритм, описанный в работе [64], наиболее подходит для решения задач в области оптической диффузионной спектроскопии в случае однородных плоскопараллельных тканей. При этом алгоритм позволяет существенно сократить время моделирования для получения результатов, сравнимых с результатами обычного метода Монте-Карло.

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

3. Другие методы повышения эффективности алгоритма моделирования.

Наряду с применением методов понижения дисперсии, исследователи пытаются снизить время моделирования распространения излучения методом Монте-Карло другими способами.

В процессе трассировки пакета фотонов выполняется определенный набор вычислений, в основе которого лежит работа с логарифмическими и тригонометрическими функциями. В работе [65] все такие функции были аппроксимированы полиноминальными и рациональными версиями, что позволило сократить время вычислений примерно в 4 раза на процессоре семейства Intel Pentium 4. Аппроксимации были подобраны таким образом, чтобы погрешность вычисления каждой функции не превышала 1%, при этом общая погрешность результатов моделирования не превышает 2%.

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

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

Третий подход, описанный в работе [67], предлагает использование модификации стандартного метода Монте-Карло для решения уравнения переноса излучения, которая направлена на возможность моделирования движения фотонов во времени. Принципиальным ограничением описываемого в работе подхода также является ограниченная область применения. Алгоритм работает только для бесконечных и полубесконечных сред.

1.3.6. Реализация метода Монте-Карло на вычислителях с параллельной архитектурой.

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

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

Очевидным подходом к ускорению алгоритма моделирования является использование возможностей современных многоядерных процессоров [55, 68, 69]. В силу высокой степени параллелизма задачи, применение многоядерных CPU позволяется существенно сократить время моделирования. Эффективность распараллеливания в этом случае составляет более 90% [68, 69].

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

В рассматриваемой задаче обмен данными между процессорами нужно осуществлять всего два раза – в начале и в конце процесса моделирования. Это означает, что алгоритм решения уравнения переноса излучения методом Монте-Карло должен показывать хорошую эффективность распараллеливания и близкую к линейной масштабируемость на системах с распределенной памятью, что подтверждается соответствующими результатами исследователей [68, 69].

Заметим, что наиболее эффективные параллельные алгоритмы моделирования распространения излучения для CPU сочетают в себе как работу с многопроцессорными системами, так и работу с многоядерными CPU в рамках одного узла кластера [69].

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

Например, в работе [70] алгоритм моделирования распространения излучения в многослойных биологических тканях с плоскопараллельной геометрией слоев был реализован на ПЛИС (программируемая логическая интегральная схема). Область применения полученного решения – фотодинамическая терапия раковых опухолей. Использование ПЛИС позволило сократить время вычисления в 28 раз по сравнению с однопоточной CPU версией. При этом ПЛИС потребляет в 25 раз меньше мощности.

Существенным недостатком предлагаемого решения является сложность процесса реализации алгоритма на ПЛИС. Так разработка описанной выше интегральной схемы заняла 12 человеко-месяцев. А в случае реализации алгоритма моделирования для гетерогенных сред с произвольной геометрией слоя с силу существенного усложнения задачи длительность процесса разработки значительно увеличится.

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

Технические особенности переноса алгоритма для моделирования распространения света на GPU для случая однородной плоскопараллельной полубесконечной среды описаны в работе [71]. Ускорение GPU версии с одинарной точностью по сравнению с однопоточной CPU программой двойной точности здесь составляет 1080 раз. Следует отметить, что приведенные в работе оценки ускорения не следует рассматривать в качестве ориентира для задач моделирования распространения излучения методом Монте-Карло. Вопервых, сравнение производилось для версий с разной точностью вычислений, а вовторых, решалась задача моделирования в достаточно простом частном случае – для полубесконечной однородной среды.

Реализация алгоритма моделирования для более общего случая – трехмерной гетерогенной среды с произвольной геометрией приведена в работе [72]. Геометрия слоев исследуемого объекта описывается набором вокселей. Ускорение GPU версии относительно однопоточной CPU реализации составляет от 75 до 300 в зависимости от постановки задачи. Однако в работе опять же происходит сравнение версии для графического процессора с одинарной точностью и версии для центрального процессора двойной точности, что не совсем корректно.

Результаты переноса алгоритма моделирования для случая гетерогенных сред с произвольной геометрией слоя с триангуляцией границ описаны в работе [56]. Здесь производилось сравнение GPU реализации и однопоточной CPU версии с одинаковой точностью вычислений. В работе продемонстрировано ускорение GPU реализации в 4 – 11 раз в зависимости от параметров задачи, что позволяет говорить о графических процессорах как полезном инструменте для ускорения вычислений рассматриваемого алгоритма.

В работе [73] описывается реализация алгоритма моделирования для распределенных систем (имеются в виду не кластерные системы, а набор отдельных компьютеров, объединенных глобальной сетью, например Internet). Предполагается, что каждый узел распределенной системы имеет в своем распоряжении NVidia GPU с поддержкой вычислений общего назначения. В качестве технологии построения распределенной системы рассматривается одноранговая (P2P) сеть. В работе показано, что использование P2P сетей в качестве дополнительного уровня для распараллеливания процесса моделирования переноса фотонов может существенно сократить время вычислений, особенно для действительно больших задач.

1.3.7. Метод Монте-Карло в задаче ОДС.

Метод оптической диффузионной спектроскопии, описанный в §1.1, позволяет определять биомедицинское состояние биологической ткани посредством измерения ее параметров поглощения и рассеяния. Обычно необходимо проводить исследование определенной области биоткани, для чего используется нужная конфигурация источника и детекторов. Задача моделирования – определить требуемую конфигурацию, исходя из оптических характеристик исследуемого объекта и расположения нужной области.

Таким образом, метод моделирования распространения излучения в различных средах для решения задачи ОДС должен обеспечивать следующие возможности:

Работу с биологическими тканями, которые относятся к классу сильно рассеивающих сред;

Поддержку достаточно сложной геометрии слоев исследуемого объекта;

Возможность моделирования сигнала на детекторах;

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

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

Сеточные методы моделирования переноса фотонов позволяют решать значительно более общие задачи и удовлетворяют почти всем требованиям из приведенных выше.

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

Единственным подходом, который удовлетворяет всем выше приведенным требованиям, является моделирование методом Монте-Карло. Применимость данного подхода для решения задач оптической диффузионной спектроскопии подтверждается рядом исследователей [54, 74, 75].

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

Описанный подход представлен, например, в работе [74], где в качестве объекта исследования рассматривается голова человека. Показано, что степень проникновения света достаточна для проведения мониторинга зон коры головного мозга.

Следующий шаг решения задачи ОДС – проверка и, при необходимости, корректировка найденных на предыдущем этапе параметров конфигурации «источник-детектор»

путем проведения дополнительного моделирования. Область мониторинга задается в качестве одного или нескольких дополнительных слоев. Проводится моделирование с найденными ранее параметрами при изменении оптических свойств области мониторинга и анализируется изменение сигнала на детекторе. Если сигнал меняется существенно – в рамках данного этапа задача считается решенной.

Заключительный этап – проверка полученных результатов путем построения тестовой установки и проведения реальных физических экспериментов.

Важным аспектом при проведении функциональной диагностики методом ОДС для конкретного пациента является учет его физиологических особенностей, в частности, геометрии исследуемых тканей. При проведении функциональной диагностики головного мозга существует возможность получать геометрию тканей головы с помощью МРТ и КТ, и уже с учетом полученных данных проводить моделирование. Такой подход описан в работах [54, 75].

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

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

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

Задачами диссертационного исследования являются:

1. Разработка более эффективного с точки зрения быстродействия модифицированного метода Монте-Карло для решения задач структурного анализа и синтеза системы источников и детекторов зондирующего излучения;

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

3. Исследование эффективности разработанного метода и параллельных алгоритмов на примере решения задачи моделирования ОДС;

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

5. Применение полученного программного комплекса для моделирования функциональной диагностики головного мозга человека методом ОДС.

2. МОДИФИЦИРОВАННЫЙ МЕТОД МОНТЕ-КАРЛО ДЛЯ

РЕШЕНИЯ ЗАДАЧ СТРУКТУРНОГО АНАЛИЗА И

СИНТЕЗА СИСТЕМЫ ИСТОЧНИКОВ И ДЕТЕКТОРОВ

ЗОНДИРУЮЩЕГО ИЗЛУЧЕНИЯ

2.1. Постановка задачи Рассмотрим постановку задач структурного анализа и синтеза системы источников и детекторов зондирующего излучения в контексте моделирования оптической диффузионной спектроскопии (ОДС).

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

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

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

Используемые в диссертации модель описания биологических тканей и модель распространения зондирующего излучения являются обобщением соответствующих моделей, описанных в работах [43, 45, 54, 56, 64], с целью наиболее эффективного решения задач оптической диффузионной спектроскопии. Опишем их подробнее.

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

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