Содержание к диссертации
Введение
Глава 1. Обзор математических моделей, методов решения ОЗЭ и особенностей формирования биопотенциалов мозга 9
1.1. Математические модели для решения прямой задачи ЭЭГ 9
1.2. Методы решения ОЗЭ 34
Глава 2. Решение ОЗЭ при однодипольном моделировании источника ЭЭГ и одномоментном наборе исходных данных 39
2.1. Однородная модель 39
2.2. Неоднородная модель 44
2.2.1. Разработка оптимального алгоритма поиска минимума целевого функционала 51
2.2.2. Методы поиска минимума функционала с использованием базы готовых решений прямой задачи 65
2.2.3. Результаты исследований на основе однодипольной локализации 72
Глава 3. Решение ОЗЭ при двухдипольном моделировании источника ЭЭГ и одномоментном наборе исходных данных 85
3.1. Постановка задачи для однородной и неоднородной моделей головы 86
3.2. Выбор метода поиска целевого функционала 87
3.3. Результаты экспериментов с двухдипольной моделью источника 93
Глава 4. Способы повышения степени адекватности моделей для решения ОЗЭ 110
4.1. Решение прямой задачи для однородного сфероида 110
4.2. Уточнение модели источника ЭЭГ с использованием квадрупольнои аппроксимации генератора электрической активности 116
4.3. Использование временной информации ЭЭГ активности при решении ОЗЭ 132
4.3.1. Математические модели пространственно-временной активности 134
4.3.2. Определение количества доминантных источников при пространственно-временном моделировании активности мозга 137
4.3.3. Постановка и решение ОЗЭ для пространственно-временных моделей 140
4.4. Статистическая обработка результатов однодипольной локализации. Способ повышения точности локализации 146
Глава 5. Результаты клинических испытаний разработанных методик локализации источников электрической активности мозга 153
Общие выводы 161
Список литературы 164
Приложение 171
- Разработка оптимального алгоритма поиска минимума целевого функционала
- Результаты исследований на основе однодипольной локализации
- Уточнение модели источника ЭЭГ с использованием квадрупольнои аппроксимации генератора электрической активности
- Статистическая обработка результатов однодипольной локализации. Способ повышения точности локализации
Введение к работе
Правильная и своевременная диагностика состояния и заболеваний головного мозга всегда была и остается до настоящего времени одной из ключевых задач современной медицины, поскольку этот орган выполняет огромное количество жизненно важных функций. Однако, несмотря на большое количество существующих методов диагностики мозга, он остается наименее изученной частью человеческого организма. Это связано, прежде всего, со сложностью его строения, которая определяется большим количеством нейронов в мозге, сложностью взаимодействий между нервными центрами в различных функциональных состояниях мозга, отсутствием достоверных знаний о механизмах генеза ЭЭГ и пространственно-временных взаимосвязях в нейронных сетях головного мозга, с отсутствием единой точки зрения о природе, структуре и характере источников потенциалов во многих случаях патологий конкретных пациентов. По этой причине изучению мозга сегодня уделяется большое внимание. С этим неразрывно связана потребность в достоверных и точных методах изучения и измерения его функциональных характеристик, поскольку эти методы являются инструментарием для изучения мозга.
Метод электроэнцефалографии занимает важное место среди таких неинвазивных методов функциональной диагностики состояния мозга, как магнитоэнцефалография (МЭГ), эхоэнцефалография (ЭхоЭГ), реоэнцефалография (РЭГ), доплерография. Несмотря на появление таких диагностических методов, как компьютерная томография (КТ), магнитно-резонансная томография (МРТ), позитронно-эмиссионная томография (ПЭТ), интерес врачей и исследователей к электроэнцефалографии и методам ее анализа в последнее время возрастает. Это связано, прежде всего, с тем, что методы нейровидения основаны на применении вредных для мозга излучений. Различие диагностических результатов методов нейровидения и электроэнцефалографии заключается в том, что в первом случае это получение информации о морфологии очагов патологии, а во втором -информации о функциональных очагах. Электроэнцефалография незаменима на ранних стадиях формирования очагов патологии, когда органические изменения незначительны настолько, что не выявляются методами нейровидения, а биоэлектрическая активность содержит информацию о наличии патологии. Кроме того, преимуществами электроэнцефалографии перед методами нейровидения являются: существенно более низкая стоимость аппаратуры и, следовательно, обследования, безвредность ввиду отсутствия вредных для мозга излучений (рентгеновское в случае КТ и магнитное в случае МРТ), высокое временное разрешение.
Среди методов анализа данных электроэнцефалографии (корреляционный, спектральный, когерентый, частотный, фазово-частотный, картирование) в последние годы наибольший интерес ученых привлекает методика трехмерной локализации источников электрических потенциалов, генерируемых мозгом [53,54,68,83,94], которая позволяет получить информацию о положении функциональных очагов. Эта методика основана на решении обратной задачи электроэнцефалографии (ОЗЭ), для которой исходными данными являются потенциалы, регистрируемые на поверхности скальпа, а результатом -пространственные характеристики источника-генератора этих потенциалов. Это математическая задача, для решения которой необходимы значительные вычислительные ресурсы ЭВМ.
В 70-80-е годы ввиду недостаточных вычислительных возможностей ЭВМ, доступных в клинической практике, программы для локализации источников электрической активности мозга были основаны на грубых упрощенных моделях [10,28,52,92] и приближенных методах решения ОЗЭ [61,69], которые позволяли с учетом допустимых затрат времени лишь приблизительно (с точностью до 1ч-2 см) оценивать характер распределения очагов патологической активности. Это приводит к получению недостоверных результатов и дискредитирует метод в глазах врачей. В настоящее время в связи с широким внедрением в клиническую практику высокоскоростных ЭВМ должны ставиться качественно другие цели - разработка методов локализации высокой точности, работающих в режиме реального времени, на основе исследования математических моделей и методов решения ОЗЭ с их помощью. Получение достоверных и проверенных знаний о механизмах генеза ЭЭГ, разработка адекватных математических моделей электрической активности мозга являются актуальной задачей. Сегодня приоритет нужно отдавать использованию в практике моделей, обеспечивающих наличный уровень знаний о механизмах генеза ЭЭГ, и повышению точности методов решения, что обусловлено достаточным уровнем развития вычислительной техники, использующейся в широкой клинической практике.
Поэтому целью данной работы является разработка эффективных методов решения проблемы локализации источников электрической активности головного мозга в норме и при патологии, что требует разработки адекватных биообъекту математических моделей, в том числе моделей электрической активности источников биопотенциалов, обеспечивающих высокую точность, а также реализация этих методик в виде эффективных компьютерных программ. Можно выделить два главных направления для достижения поставленной цели: исследование и разработка математических моделей головы, моделей электрической активности и моделей источников биопотенциала, являющихся наиболее точными с точки зрения имеющихся знаний о строении и функционировании мозга; разработка эффективных методов решения ОЗЭ и реализация их в виде программного обеспечения.
В рамках обозначенных направлений задачами данной работы являются:
Исследование существующих моделей, использующихся при локализации источников электрической активности мозга, и возможностей их совершенствования с учетом покровов мозга, а так же с учетом реальной геометрии черепа.
Разработка эффективного метода решения ОЗЭ для модели одного и двух подвижных дипольных источников.
Разработка метода решения ОЗЭ на основе мультипольного моделирования источников электрического потенциала.
Разработка и исследование метода решения ОЗЭ для пространственно-временной модели электрической активности мозга.
Исследование влияния количества регистрирующих электродов на точность локализации.
Реализация оптимального метода решения ОЗЭ в виде программного обеспечения в автоматизированном программно-аппаратном комплексе для его использования в клинической практике.
В работе получены следующие новые результаты:
Разработана дипольно-квадрупольная модель источника биопотенциала мозга, являющаяся более адекватной, чем дипольная модель и позволяющая повысить точность локализации корковых источников потенциала на 0.5-^1 см в сравнении с дипольной локализацией.
Разработана сфероидальная однородная модель головы, более адекватная, чем сферическая однородная.
Впервые установлены характерные особенности структуры функционала, определяющего решение ОЗЭ, и наличие его дополнительных ложных минимумов. Разработан метод поиска глобального минимума, соответствующего точному решению задачи локализации. Исследована зависимость структуры функционала от количества и расположения электродов.
Разработан алгоритм решения прямой задачи для неоднородной сферической модели, основанный на базе компонент скальповых потенциалов.
5. Разработан алгоритм определения количества доминантных источников при пространственно-временном моделировании электрической активности мозга, представляющий собой развитие метода главных компонент.
Практическая значимость работы заключается в следующем:
Разработаны эффективные алгоритмы решения ОЗЭ для однодипольной, двухдипольной, дипольно-квадрупольной моделей источника и пространственно-временной модели электрической активности мозга, отличающиеся от известных аналогов адекватностью и точностью.
Исследована взаимосвязь количества электродов с точностью локализации. Предложено использование плотных электродных сеток при анализе электроэнцефалограмм методами пространственной локализации.
Разработан метод статистической обработки результатов однодипольной локализации, позволяющий повысить точность и степень достоверности заключения по результатам проведенного анализа.
Создана программа локализации источников электрической активности мозга. Она является составной частью аппаратно-программного комплекса "Нейрон-Спектр" ("НейроСофт", г. Иваново), использующегося в клинической практике для регистрации и анализа биопотенциалов головного мозга.
Результаты диссертации внедрены в клиническую практику Ивановской и Владимирской областных клинических больниц, городской больницы № 7 г. Иваново, Центрального института экспертизы трудоспособности инвалидов (г. Москва) и других лечебных учреждений (около 10 внедрений).
Основное содержание диссертации заключается в пяти главах.
В первой главе приведен подробный аналитический обзор литературы по теме реконструкции пространственных источников биопотенциала в мозге. Отдельно обсуждаются существующие математические модели и методы решения ОЗЭ. Рассмотрены базовые принципы возникновения биопотенциалов мозга, механизм статистической суммации потенциалов отдельных нейронов и некоторые статистические закономерности их распространения. Вторая глава представляет собой подробное исследование однодипольной модели источника потенциала в однородной и неоднородной сферических моделях головы. В ней описаны сами математические модели, проанализированы методы минимизации целевого функционала и представлены результаты численных исследований созданной методики однодипольной локализации для однородной и неоднородной моделей головы. В третьей главе проведено исследование двухдипольной модели источника в однородной и неоднородной сферических моделях головы. В ней представлены результаты исследования структуры функционала в случаях стандартной ЭЭГ и при увеличенном количестве электродов. Четвертая глава содержит описание и результаты численного исследования пространственно-временных моделей электрической активности мозга. В ней описана дипольно-квадрупольная модель источника потенциала и методика локализации таких источников, а также приведен сравнительный анализ структуры функционала в случаях дипольного и дипольно-квадрупольного источника. Кроме того, в этой главе описаны сфероидальная однородная модель головы и метод статистического уточнения результатов локализации. В пятой главе представлены и проанализированы результаты клинических испытаний разработанных методик локализации. В общих выводах по работе перечислены основные результаты и выводы.
Разработка оптимального алгоритма поиска минимума целевого функционала
В связи с этим становится очевидной необходимость определения количества источников при использовании пространственно-временной модели. Методы, предназначенные для определения количества источников известного сигнала, разработаны и используются в предложениях, не связанных с электрофизиологией [91,57]. В [91] такой подход представлен в рамках алгоритма классификации множественного сигнала (Multiple Signal Classification - MUSIC). Автор использует этот метод для определения параметров множественных волновых фронтов, регистрируемых антенной, на основе измеряемых потенциалов. В [79] алгоритм MUSIC, включающий в себя методику определения количества источников сигнала, как составную часть, распространен на анализ данных МЭГ. При этом авторы используют многодипольную модель генерации ЭЭГ, то есть одновременно несколько активных диполей всех трех вышеописанных типов (подвижные, с фиксированным положением, с фиксированным положением и ориентацией). Методика накладывает ограничения на количество диполей в модели, которое не может превышать количество датчиков, регистрирующих магнитное поле на поверхности скальпа. Для анализа отбираются несколько одномоментных сечений МЭГ. Формируется матрица исходных данных, в которой столбцы - потенциалы магнитного поля, записанные в момент времени ten точках регистрации, строки - значения потенциала магнитного поля в і точке регистрации в m последовательные моменты времени. Количество дипольных источников этих исходных данных предлагается оценивать, как равное числу собственных значений матрицы исходных данных, превышающих некоторое пороговое значение. Ключевым предложением в этом методе является то, что разные диполи генерируют линейнонезависимые компоненты сигнала во всех анализируемых срезах. При помощи разделения собственных чисел производится выделение из исходных данных полезного сигнала и шума. Очевидно, что такое разделение возможно лишь при качественной записи МЭГ. В [78,96] для определения количества источников МЭГ используется аналогичная методика, отличающаяся тем, что авторы используют сингулярное разложение матрицы исходных данных вместо разложения по собственным значениям. Способ численной оценки значения, разделяющего полезный сигнал и шум, рассматривается в [57] в приложении к проектированию антенн. В [63] сделана попытка преодоления ограничения, накладываемого на количество локализуемых источников, принятого в [79]. В экспериментах автор использует однородную сферическую модель и 37 датчиков магнитного поля. Эксперименты показывают неработоспособность такого способа определения количества доминантных источников в случае отсутствия четкой границы между "большими" и "малыми" собственными числами. Этот недостаток свойственен и методу MUSIC. Метод MUSIC не работоспособен когда: исходные данные содержат сильный шум, затрудняющий их разделение на сигнал и шум; взятые для анализа срезы МЭГ сильно коррелируют друг с другом; локализуемые источники расположены вблизи друг друга.
На аналогичной процедуре определения количества источников основывается метод анализа главных компонентов (Principal Components Analysis - РСА). Впервые этот метод применен в работе [78], в которой он используется для локализации источников зрительных вызванных потенциалов (ВП). Известно, что у этих ВП есть несколько различных источников. Соответственно, однодипольный анализ в данном случае не отразит точной картины генерации ВП. Авторы применяют метод РСА на трехслойной сферической неоднородной модели с 24-электродной сеткой. Основой алгоритма РСА является выделение главных компонентов из матрицы исходных данных U размерности NxM, где N - количество электродов в схеме записи ВП, М — количество временных срезов, взятых для анализа. Количество главных компонентов, вносящих основной вклад в общий сигнал, определяет количество активных источников. В экспериментах по локализации источников зрительных ВП авторами установлено, что не более двух диполей генерируют 95% общего сигнала. Остальные 5% относятся к шуму. Такое отношение сигнал/шум может быть получено только при большом количестве усреднений (более 200) для записи ВП. Главные компоненты матрицы U определяются сингулярным разложением где V - ортонормальная матрица, описывающая временные функции, соответствующие главным компонентам, Р - ортонормальная матрица, описывающая пространственное распределение главных компонентов, S - диагональная матрица сингулярных чисел, W -матрица "коэффициентов загрузки". Каждый элемент WL. представляет собой вклад к главной компоненты в потенциал, записанный j электродом. Величины сингулярных чисел интерпретируются как мощности соответствующих главных компонентов. Количество главных компонентов/? определяется уровнем полезного сигнала а где Skk - сингулярные числа матрицы U. Количество главных компонентов и есть количество активных источников. Определение параметров диполей производится методом наименьших квадратов для минимизации функционала ошибки. Недостатком данного метода является невозможность во многих случаях четко отделить существенные и несущественные главные компоненты, что приводит к ошибке определения количества источников и лишает смысла дальнейшее решение ОЗЭ.
В работах [96,75] предприняты попытки усовершенствовать метод РСА, а в частности его этап определения количества доминантных источников ЭЭГ. Авторы [96] отмечают схожесть алгоритмов MUSIC и РСА на этом этапе. Однако они работают не всегда правильно, особенно в случаях, когда невозможно отделить несущественные собственные числа и соответствующие им главные компоненты.
В ранних работах [28,74,95] предлагается определять количество источников на основе визуального анализа потенциальных карт на поверхности скальпа. Однако, в этих исследованиях рассматривается не более двух источников и очевидно, что такой подход не работает в случае наличия мультидипольной активности. В то же время метод картирования потенциала и других характеристик электрической активности мозга - это признанный и широко использующийся в клинической практике инструмент. В дальнейшем речь пойдет о потенциальных картах, хотя кроме них при анализе ЭЭГ используются частотные и амплитудные. Карта позволяет визуально оценить характер распределения потенциалов на поверхности скальпа. Точность карты потенциалов зависит от плотности электродной сетки, с помощью которой она вычисляется. Достоверно известны значения потенциала лишь в точках регистрации. Значения в промежуточных точках вычисляются с помощью различных методов интерполяции: линейной [86], поверхностными сплайнами [86], сферическими сплайнами [85]. В подавляющем большинстве случаев ЭЭГ записывается по стандартной международной схеме "10 - 20" [57]. В этом случае количество электродов не превышает 24, а среднее значение между любыми двумя ближайшими 6 см. Общий и достаточно подробный характер карты представим с помощью ограниченного количества электродов. Однако, результаты экспериментов [81] показывают, что значение картируемой величины (в данном случае - значений лапласиана) в промежуточных точках зависят от плотности электродной сетки. Следовательно, никакой метод построения карты не дает возможности сколько-нибудь точно аппроксимировать значения картируемой величины и особенно в случае разреженных электродных сеток. Несмотря на это, точности традиционной сплайн-интерполяции достаточно для визуально анализируемых карт. Однако, для методов трехмерной локализации источников электрической активности точность карты скальповых потенциалов имеет важнейшее значение. Некоторые исследователи [79,64,51] считают, что в отсутствии плотных сеток электродов можно заменить их картой, полученной интерполяцией.
Результаты исследований на основе однодипольной локализации
Канал регистрирует разность потенциалов записываемых активным и пассивным (референтным) электродами. Набор описаний пар электродов представляет собой схему отведения (записи) ЭЭГ. В клинической практике записываются ЭЭГ с количеством каналов 8-21. В данной работе эксперименты проводились с использованием 19-канальной схемы отведения с объединенным ушным электродом, принятым в качестве референтного. Это означает, что для 19 активных электродов пассивным является некий гипотетический электрод, потенциал которого вычисляется как среднее арифметическое потенциалов, регистрируемых ушными электродами. Потенциал, регистрируемый им, получается путем усреднения потенциалов 2-х реальных ушных электродов. Такая схема выбрана ввиду наибольшей распространенности в медицинской практике.
Процесс решения задачи минимизации описывается обобщенным алгоритмом (Рис.2.5), каждый этап которого представляет собой отдельную исследовательскую подзадачу.
Первым шагом алгоритма оптимизации является выбор начальных условий или точки старта. Проведенные численные эксперименты позволили установить, что скорость сходимости итерационного процесса минимизации в большой степени зависит от этого выбора, причем близость точки старта к решению не гарантирует высокой скорости сходимости и точности решения. Определить и обосновать правила выбора начальных условий позволяет знание характера распределения значения функционала (31) внутри области, моделирующей мозг. На Рис.2.6 (а-г) представлены трехмерные графики, представляющие собой поверхности, описывающие распределение функционала по сферам нескольких фиксированных радиусов (R=4.5H-7.5 СМ С шагом 1 см). На этих графиках сферы разных радиусов представлены в виде прямоугольных разверток. Координаты у,х развертки соответствуют угловым сферическим координатам 6, СС. Графики показывают распределение функционала внутри модели головы для случая, когда дипольныи источник с координатами (-4.38, -3.71, 1.79) и тангенциальной ориентацией расположен в однородной модели головы.
На Рис. 2.4 представлено распределение целевого функционала на поверхности постоянного радиуса вдоль линий проходящих через дипольныи источник при 19 и 64-канальной регистрации ЭЭГ. Важно отметить, что в случае 64-канальной ЭЭГ щелевидность минимума существенно меньше, а также меньше ложных минимумов и они ближе к решению. Щелевидность определяется по разности между минимальным и максимальным значением второй производной по различным направлениям в точке минимума. В результате анализа характера распределения функционала для различных положений дипольного источника сделаны следующие выводы: 1. На поверхности с R 4 см отчетливо виден один глобальный минимум с пологим характером. 2. Характер рельефа поверхностей усложняется по мере приближения R к радиусу модели. Это выражается в появлении нескольких локальных максимумов. 3. Локальные максимумы поверхностей больших радиусов находятся под электродами и имеют максимальные значения на поверхностях с R = 7.5 -9 см. 4. В случаях расположения дипольного источника вблизи поверхности модели могут обнаруживаться локальные минимумы. 5. Глобальный минимум имеет форму впадины вытянутой вдоль радиуса модели независимо от ориентации диполя. На основании проведенного анализа определяется область наилучшего исходного положения диполя. Глобальный минимум отчетливо проявляется на слоях с R 1.5 см. и там доминирует. Таким образом, в качестве оптимальных начальных условий предпочтительно задавать источники с R Зсм, то есть в той области, где функционал изменяется плавно, возможно направление только на один минимум и нет сложных поверхностей перегиба. Это положение получило практическое подтверждение в результате численных экспериментов, в которых точки старта находились в различных областях модели для одинаково расположенных дипольных источников. Следующим шагом алгоритма (Рис.2.5) является определение проекций момента диполя на декартовые оси координат или ориентации момента и его модуля. Эту подзадачу можно не решать, если включить эти параметры источника в число переменных модели. Однако такой подход приводит к увеличению количества степеней свободы в 2 раза. В случае однодипольного моделирования источника поля такое увеличение (от 3 до 6) не оказывает существенного влияния на качество решения ОЗЭ при использовании 19 электродной схемы записи ЭЭГ. С увеличением количества дипольных источников модели до двух количество переменных - 12 - становится сравнимым с количеством исходных данных - 19, что существенно затрудняет поиск минимума F (31). Поэтому более предпочтительным является вычисление проекций дипольного момента исходя из условия минимальности функционала при фиксированных координатах источника. Для этого функционал F (31) в виде Значения проекций момента диполя, получаемые решением СЛАУ (32) минимизируют функционал (31) при заданных координатах диполя. Таким образом, количество переменных модели уменьшается в два раза и становится равно 3 или 6 для случаев 1 или 2 дипольных источников соответственно. Наиболее важным этапом алгоритма (Рис.2.5) является определение направления, в котором выбирается следующее положение диполя. Способ вычисления вектора направления - это основной фактор, определяющий скорость и точность решения задачи минимизации. Простейшее решение этой подзадачи реализуется в виде нового созданного в результате численных экспериментов алгоритма без вычисления производных, названного "клетка". Он особенно эффективен, когда нет возможности вычислять производные целевой функции в явном виде. Выбор оптимального направления происходит следующим образом: 1. вокруг центральной точки строится клетка из 6 точек, каждая из которых получается из центральной путем приращения одной из ее координат со знаком "+" или "-"; 2. во всех 7 точках вычисляются значения функционала невязки; 3. за направление минимизации принимается вектор, соединяющий центральную клетки и точку в которой функционал минимален. Критерием остановки этого метода минимизации является событие, когда функционал минимален в центральной точке клетки. В рамках алгоритма "клетка" необходимо выбрать оптимальный шаг клетки, то есть приращение координат, определяющее размеры области охватываемой клеткой, а также определить степень необходимости применения методов одномерной оптимизации. При выборе шага клетки необходимо учитывать следующее: 1. Большой шаг клетки уменьшает точность решения и увеличивает вероятность пропуска минимума в случае сложного характера поверхности. 2. Маленький шаг клетки снижает скорость процесса оптимизации и таким образом делает предпочтительным использование метода одномерной минимизации. Эксперименты показали, что оптимальной является двухэтапная методика, при которой на первом этапе шаг "клетки" составляет 5 мм, а на втором - 1-2 мм. Для решения данной задачи были разработаны и протестированы модификации алгоритма "клетка": 1. Вокруг центральной точки строится куб с 27 узлами. В узлах куба вычисляются значения функционала, по которым строится интерполирующий полином. Точка, минимизирующая этот полином в пределах куба принимается центром следующего куба. 2. Вокруг центральной точки достраиваются не 6, а 14 точек, образующих решетку куба. В остальном алгоритм не отличается от исходного. С помощью этих модификаций алгоритма удается лишь незначительно повысить точность решения задачи оптимизации, но достигается это за счет увеличения времени решения.
Уточнение модели источника ЭЭГ с использованием квадрупольнои аппроксимации генератора электрической активности
Из всех возможных условий достижения минимума при решении данной задачи проверялись только эти три как наиболее приемлемые по ее условиям. В ходе численных экспериментов установлено, что первый критерий неудовлетворителен из-за разнообразия типов конфигурации поверхности функционала в окрестности его минимума. Это затрудняет определение предельного значения модуля вектора направления минимизации.
В некоторых ситуациях (особенно при расположении источников близко к центру) область минимума обширная и плоская. В такой ситуации модуль вектора направления минимизации становится близким к 0 задолго до минимума (5-7 мм).
Второй критерий, также как и первый, не работает корректно в случае плоского обширного минимума, но он универсален для всех методов определения направления минимизации. Это условие достижения минимума следует очень осторожно применять в ситуации, когда минимум имеет щелевидный характер. В таком случае градиентный метод, например, дает направление поперек щели, и шаг может оказаться очень малым, хотя реально минимум далеко. Подобные ситуации имеют место в случае расположения дипольных источников на внешних радиусах модели. Учитывая вышеизложенные соображения и результаты экспериментов, решено комбинировать критерии минимальности длины шага и убывания функционала. Условие остановки по длине шага (нескольких шагов) полезно тем, что позволяет избежать большого количества итераций, существенно не уточняющих решение.
Отдельно следует упомянуть ситуацию, когда текущее положение дипольного источника близко к границе мозга, а вектор направления минимизации направлен за пределы модели. В таком случае минимум считается достигнутым в последней точке. В многодипольной модели прекращается движение только одного этого диполя.
В процессе поиска минимума функционала (31) многократно решается прямая задача, то есть вычисляются значения скальповых потенциалов, генерируемых ТДИ в различных точках внутри области моделирующей мозг. Этот процесс достаточно прост и относительно быстр в случае расчетов с использованием однородной шаровой модели. Однако, решение этой же задачи для более адекватной реальности трехслойной неоднородной шаровой модели требует существенно большего времени ввиду отсутствия соответствующих конечных формул (см. раздел 2.1.2). В таких условиях решение ОЗЭ по алгоритму (Рис.2.5) требует больших вычислительных затрат. Они еще больше возрастают при использовании методов выбора направления минимизации, требующих вычисления производных, так как в этом случае необходимо использовать конечно-разностные аналоги производных. Поэтому при расчетах по неоднородной модели может применяться метод, основанный на использовании предварительно созданной базы решений прямой задачи. В базе хранятся значения (pxi, Pyi, Pzi, где - порядковый номер положения, в котором находится дипольный источник, генерирующий эти потенциалы. Положения ДИ определяются узлами пространственной решетки, расположенной внутри области, моделирующей мозг. Эта решетка может соответствовать декартовой либо сферической системам координат. Недостаток сферической решетки в ее неравномерной плотности в зависимости от радиуса. В тоже время декартова решетка неравномерна вблизи границы области, моделирующей мозг. Поиск минимума функционала с использованием базы производится простым перебором. Недостатком такого метода является необходимость создания очень "подробной" базы.
В противном случае, если узлы решетки базы расположены далеко друг от друга, может происходить потеря решения (Рис. 2.9), что было подтверждено экспериментальным путем. Для неоднородной шаровой модели была расчитана база со сферической решеткой (RMo3ra= 8 см, Ячерепа = 8.7 см, Рчжальпа = 9.2 см). В решетке по радиусу 8 узлов, по углу в - 24, по углу а - 48 (всего 9217 узлов). При некотором положении тестового диполя (6.7; 1.3; 1.35) с ориентацией вдоль оси ОХ функционал достигал своего минимального значения в точке X , отстоящей от тестового диполя на 2.8 см. При этом ближайший к тестовому диполю узел находится на расстоянии 0.5 см от него. Для проверки гипотезы наличия ложного минимума из всех трех точек (тестовая, узел ближайший к тестовой точке, узел X ) проводился градиентный поиск с фиксированным шагом 0.1 см. Из точки, в которой расположен ТД, смещения не произошло. При старте из узла ближайшего к ТД метод остановился в 0.17 см от ТД. При старте из узла X метод в качестве минимума определил точку удаленную от ТД на 2.95 см. Результаты этого эксперимента подтверждают наличие локального минимума целевого функционала и доказывают необходимость использования более подробной базы потенциалов. Однако увеличение размеров базы влечет за собой увеличение времени затрачиваемого на вычисление функционала во всех узлах. При количестве узлов базы превышающем 12000-15000 более эффективными по скорости оказываются методы без использования базы. В такой ситуации возможно применение алгоритма перебора, но с двухступенчатой базой. База первой ступени строится по редкой решетке, второй - по частым решеткам расположенным вокруг узлов решетки базы первой ступени. После определения узла с Fmjn по базе первой ступени расчеты продолжаются с использованием более подробной базы второй ступени, соответствующей узлу с Fmjn. При таком подходе удается сократить время расчетов, но вероятность пропуска точки минимума остается и даже возрастает, так как решетка базы первой ступени должна быть достаточно редкой, а также существенно возрастает объем базы (за счет базы второй ступени). В проведенных экспериментах для построения базы первой ступени использовалась решетка с 2250 узлами, вокруг каждого из которых создавалась база второй ступени из 1000 узлов. В результате объем базы второго уровня составил 108 Мбайт. Расчеты показали, что при использовании такого алгоритма нередки пропуски минимума еще на этапе перебора первой базы. Возможно использование следующей модификации алгоритма с двухступенчатой базой. После окончания перебора по первой базе вокруг узла Fmjn выстраивается кубическая решетка из 27 узлов с Fmjn в центре, для чего могут быть использованы ближайшие узлы решетки базы (Рис. 2.10).
Статистическая обработка результатов однодипольной локализации. Способ повышения точности локализации
Следующим шагом по проверке точности алгоритма решения прямой задачи стало сравнение значений функционала, вычисляемых по двум различным алгоритмам: 1. Алгоритм, в котором диполь моделируется двумя точечными зарядами, расположенными на расстоянии hd друг от друга. 2. Алгоритм, в котором используется математический диполь, то есть компоненты (рх,(ру,(pz вычисляются путем дифференцирования выражения для потенциала (30). В процессе экспериментов с алгоритмом 1 база диполя hd менялась hd = 10", 10" , 10 , 10"10 см. При этом установлено, что точность решения ОЗЭ меняется несущественно при уменьшении hd ниже 10"5. Далее было проведено сравнение функционалов, вычисленных по алгоритму 1 с hd = 10"5 и по алгоритму 2. Характер изменения этих функционалов полностью совпадает, а конкретные значения различаются на доли процента. В итоге решено в дальнейшем использовать алгоритм 2, так как при этом не нужно дважды вычислять потенциал в точке наблюдения, как в алгоритме 1. Это позволяет повысить точность решения ОЗЭ, а также сократить время ее решения в 1.9 раза.
В результате вышеописанных экспериментов сделан вывод о доброкачественности и естественности появления нескольких возвышенностей на поверхностях функционала (Рис.2.6). Вполне очевидно, что такие неравномерности рельефа поверхности и наличие локальных минимумов приведут к усложнению решения задачи поиска минимума и, следовательно, к увеличению времени. В подходе, описанном в [10,28], эта проблема обходится за счет замены неоднородной модели "эквивалентной" однородной с большим радиусом. В этом случае неравномерности рельефа поверхности растягиваются и отодвигаются вслед за границей модели на большие радиусы. Естественно, в подобной ситуации процедура поиска минимума поверхности упрощается, но следствием этого достижения являются снижение степени адекватности модели и точности решения. В дальнейшем при исследовании поведения функционала проводились эксперименты не с 19, а с 64 каналами ЭЭГ. На Рис.2.7 и Рис.2.8 можно видеть поверхности функционала в случаях 32-х и 64-канальных ЭЭГ при тех же исходных данных, как при 19-канальной ЭЭГ на Рис.2.6. электродов уменьшаются размеры максимумов, располагающихся под электродами, однако соответственно увеличивается их количество. Одновременно с этим сужается область глобального минимума, что позволяет более точно выделять его. Наконец, при использовании 230-электродной регулярной (Рис.2.13) сетки замечено, что максимумы отсутствуют, либо незаметны даже на поверхностях, близких к поверхности с Rj. Характер поверхностей существенно проще, чем в случае 19-канальной ЭЭГ. Появление на поверхности Рис.2.13 (б) двух минимумов вместо одного объясняется тем, что источник с тангенциальной ориентацией расположен под несколькими электродами, распложенными цепочкой. Поэтому эти два минимума есть один разделенный пополам несколькими максимумами, слившимися в одну гряду.
В результате сделан вывод о том, что сложный характер поверхности является следствием недостатка исходных данных (количества электродов). Также установлено, что возвышенности поверхности и электроды имеют приблизительно одинаковые угловые координаты. Неполное совпадение этих координат объясняется ориентацией дипольного источника. Таким образом, сложность рельефа поверхности вблизи Rj является основной трудностью при решении ОЗЭ в случае однодипольной модели источника.
Поиск оптимального метода минимизации функционала (31) был начат с тестирования методов, использующих готовые базы прямых решений (п.О). Наиболее показательными являются результаты решения ОЗЭ с использованием трех различных баз: 1. База, построенная на основе решетки в сферических координатах. Параметры решетки: AR = у5 = 1.56см, А0 = {5, Аа = 2 - Таким образом, размеры самого большого элемента решетки 1.63x1.63x1.56 см. Количество узлов в базе 2250. 2. База, построенная на основе решетки в сферических координатах. Параметры решетки: AR = \ =0.97см, Ав = п/2А, Аа = 2%&. Таким образом, размеры самого большого элемента решетки 1x1x0.97 см. Количество узлов в базе 9216. 3. База, построенная на основе решетки в декартовых координатах. Параметры решетки: Ах = Ay = Az, z — 4 см . Размеры элемента решетки 0.5x0.5x0.5 см. Количество узлов в решетке 13400. Способ решения с использованием двухступенчатой базы (п.0) после первых экспериментов отмечен, как неподходящий ввиду того, что точность была низкой даже при расположении тестового дипольного источника на і? 4.5 см, и при этом время расчетов было больше, чем при одноступенчатой базе. Схема экспериментов описывается следующим образом: 1. Инициализация значений координат (xd,yd,zd) и проекций момента тестового дипольного источника. 2. Вычисление значений потенциала в точках регистрации (решение прямой задачи для заданных параметров источника). 3. Поиск минимума функционала (31) для посчитанных в пункте 2 значений потенциалов. Во всех вычислениях фигурируют не значения потенциалов, а разности потенциалов: потенциала активного электрода и потенциала референтного электрода. Потенциал референтного электрода в описываемых в этой работе экспериментах вычисляется, как среднее арифметическое потенциалов ушных электродов А1 и А2. Таким образом, количество исходных данных представляется разностями потенциалов 19 каналов ЭЭГ. В экспериментах с базами 1,2,3 ТД устанавливались в трех различных сферических слоях 2cM /?rf 4cM, 4см 7? 6см, 6см /? / 7.5см. для выяснения зависимости точности локализации от глубины положения ТД и от используемой базы потенциалов. Координата z ТД ограничивалась снизу значением —4 см, так как в нижних отделах головы мозга нет и, следовательно, решение ОЗЭ, полученное в этой области, не имеет важного физиологического значения. Во всех экспериментах ориентация ТД вдоль оси ОХ. В таблицах 3-5 приведены результаты этих экспериментов.
Числа в этих таблицах обозначают количество (в % от общего количества экспериментов) попаданий решений ОЗЭ, полученных с помощью алгоритмов с тремя базами разной плотности, в различные окрестности тестового источника при его расположении в шаровых слоях разной глубины. Считается, что остальные решения попали точно в тестовый источник. Все эксперименты проводились на ПК с процессором Pentium 166, RAM 64. Длительность перебора баз с целью определения узла, в котором функционал (31) минимален, составляет 2.5 с, 9 с, 14 с для 1,2,3 баз соответственно. Очевидно, что перебор по базе 1 приводит к большему количеству ошибок, чем с случае двух других баз. При этом нет зависимости точности решения от радиуса, как в случае более подробной базы 2. Это объясняется большим шагом решетки базы 1. В случае базы 2 ошибок меньше в глубине модели, то есть там, где решетка плотнее. Ошибки при использовании базы 3 меньше по абсолютной величине, но несколько больше по количеству, чем ошибки получаемые при использовании базы 2 (в таблице 5 не отражены ошибки 1..5 мм, так как декартова сетка решетки имеет шаг 5 мм и такая погрешность может восприниматься не как ошибка решения, а как погрешность метода). Наилучшим вариантом по дисперсии ошибок локализации является база 3. Однако, учитывая отмеченное выше время перебора базы, оптимальным выбором является база 2. Кроме того, абсолютную величину ошибок в этом случае легко снизить, применяя интерполяцию функционала по его значениям вокруг узла с минимумом с помощью метода описанного в п.0. С целью сокращения времени решения задачи минимизации протестированы следующие методы: 1. Метод направленного перебора - "клетка".
