Содержание к диссертации
Введение
1. Постановка задачи реконструкции распределений плотности в сечениях объектов 9
1.1. Принцип исследования 9
1.2. Преобразование Радона как основа математического описания задачи 13
1.3. Виды задач реконструкции распределений плотности 17
1.4. Использование преобразования Фурье 20
1.5. Дискретное преобразование Радона 22
1.6. Итерационный алгоритм определения длин пересечений прямой с элементами области сканирования 27
Выводы по главе 1 30
2. Усовершенствование известных методов реконструкции распределений плотности 31
2.1. Оператор обратного проецирования 31
2.2. Формулы обращения преобразования Радона 34
2.3. Метод свертки и обратной проекции и его модификации 38
2.4. Метод итераций Качмажа 43
2.5. Алгебраический алгоритм восстановления в форме, удобной для компьютерной реализации 44
Выводы по главе 2 46
3. Повышение устойчивости методов реконструкции распределений плотности и решение задач с неполными данными 48
3.1. Некорректные задачи 48
3.2. Использование метода регуляризации Тихонова для повышения устойчивости методов реконструкции 50
3.3. Способы решения задач с неполными данными 53
3.4. Способы включения дополнительной информации в методы реконструкции распределений плотности 54
3.5. Предлагаемый метод нахождения границ объекта 56
Выводы по главе 3 59
4. Численные эксперименты 60
4.1. Количественные оценки результатов реконструкции 60
4.2. Погрешность реконструкции при использовании метода регуляризации Тихонова 63
4.3. Сходимость алгебраического алгоритма восстановления 69
4.4. Реконструкция распределений плотности по ограниченному числу проекций 76
Выводы по главе 4 85
5. Реконструкция распределений плотности по данным компьютерного томографа РКТ-01 87
5.1. Преобразование веерных проекционных данных в параллельные проекционные данные 87
5.2. Определение параметров рентгеновского компьютерного томографа РКТ-01 90
5.3. Предварительные преобразования значений интенсивности 92
5.4. Анализ результатов численных экспериментов 93
Выводы по главе 5 96
Заключение 97
Литература 100
Приложение 107
- Преобразование Радона как основа математического описания задачи
- Формулы обращения преобразования Радона
- Использование метода регуляризации Тихонова для повышения устойчивости методов реконструкции
- Погрешность реконструкции при использовании метода регуляризации Тихонова
Введение к работе
Актуальность темы. Компьютерная томография (КТ) применяется для исследования внутренней структуры материальных объектов. В диссертационной работе задача определения внутренней структуры объектов рассматривается на примере рентгеновской трансмиссионной КТ.
Рентгеновская трансмиссионная КТ применяется в медицине для исследования внутренних органов [12, 18], в промышленности для контроля качества изделий [17], в геофизике для исследования мантии Земли [22], в физике для диагностики плазмы [10, 22] и т.д. Известно, что промышленные рентгеновские компьютерные томографы позволяют обнаружить дефекты размером порядка 1 мкм, благодаря чему удается повысить качество и надежность выпускаемой продукции. Срок окупаемости затрат на промышленные рентгеновские томографы во многих случаях в 5-10 раз меньше срока окупаемости технологического оборудования. К достоинствам рентгеновской трансмиссионной КТ относятся высокая скорость и невысокая стоимость исследования. В медицинских приложениях конкурентом рентгеновской трансмиссионной КТ является магниторезонансная томография, достоинством которой является отсутствие вредных воздействий на организм пациента. Однако магниторезонансная томография обладает рядом недостатков по сравнению с рентгеновской трансмиссионной КТ. К таким недостаткам относятся: позднее выявление гематом, невозможность исследования при наличии металлических имплантантов, сложность обеспечения искусственной вентиляции легких, большее влияние движений больного на качество реконструкции, невозможность исследования при наличии клаустрофобии [12, 60, 66]. Эти недостатки приводят к тому, что в некоторых случаях, например, в нейрохирургической и неврологической практике, рентге новская трансмиссионная КТ является предпочтительным методом исследования [12].
В рентгеновской трансмиссионной КТ внутренняя структура объекта характеризуется распределением плотности. В данном случае под плотностью понимается плотность электронов, т.е. количество электронов в единице объема. Хотя плотность объекта является функцией трех пространственных координат, обычно трехмерный объект представляют в виде набора сечений. Внутри каждого сечения плотность рассматривается как функция только двух переменных. Методика исследования такова: тонкий пучок рентгеновских лучей сканирует сечение объекта, испытывает частичное поглощение веществом, и изменение интенсивности рентгеновского излучения, являющееся интегральной характеристикой плотности исследуемого сечения, фиксируется детектором излучения. Затем с помощью компьютера по полученным данным производится реконструкция искомого распределения плотности в сечении объекта. Результат реконструкции выводится на дисплей в виде изображения.
Качество и скорость реконструкции распределений плотности зависят от конструктивных особенностей рентгеновского компьютерного томографа и от используемого математического аппарата. В связи с этим, одной из важных задач рентгеновской трансмиссионной КТ является совершенствование методов реконструкции распределений плотности. При этом возникают две фундаментальные проблемы. Во-первых, задача реконструкции распределений плотности является некорректной, а именно: небольшие погрешности в измеряемых данных могут привести к большим погрешностям в реконструируемой функции. Поэтому приходится прилагать усилия, направленные на повышение устойчивости методов реконструкции. И, во-вторых, для всех методов решения некорректных задач характерно противоречие между разрешающей способностью и устойчивостью. Данное противоречие состоит в том, что увеличение разрешающей способности (уменьшение систематической погрешности метода реконструкции) вызывает понижение устойчивости (увеличение случайной погрешности метода реконструкции), и наоборот. Поэтому необходимо иметь возможность устанавливать требуемое соотношение между устойчивостью и разрешающей способностью методов реконструкции распределений плотности.
В разработку методов реконструкции внутренней структуры объектов внесли свой вклад такие ученые как И. Радон, Р.Н. Брейсуэлл, Л.А. Шепп, Б.Ф. Логен, Ф. Наттерер, А.Н. Тихонов, В.Я. Арсенин, И.Н. Троицкий и многие другие. Данная диссертационная работа посвящена разработке усовершенствованных устойчивых методов и алгоритмов реконструкции распределений плотности на основе метода свертки и обратной проекции, алгебраического алгоритма восстановления, метода регуляризации Тихонова и др.
Цель диссертационной работы. Целью данной работы является совершенствование существующих методов реконструкции распределений плотности, направленное на повышение их устойчивости, разрешающей способности и скорости, а также разработка новых алгоритмов и программ на их основе. Повышение устойчивости и разрешающей способности методов реконструкции является важной задачей КТ, поскольку эти характеристики методов определяют погрешность реконструкции распределений плотности. При этом необходимо иметь возможность устанавливать требуемое соотношение между устойчивостью и разрешающей способностью. Увеличение скорости реконструкции позволяет увеличить размер реконструируемой матрицы плотности, что повышает пространственную разрешающую способность рентгеновского компьютерного томографа.
Задачи исследования. Для достижения поставленной цели, необходимо решить следующие задачи:
• Проведение анализа и сравнения существующих методов реконструкции распределений плотности в сечениях объектов.
• Применение методов регуляризации к существующим методам реконструкции распределений плотности для повышения их устойчивости.
• Исследование зависимости соотношения устойчивости и разрешающей способности от параметров методов реконструкции распределений плотности.
• Разработка методик нахождения дополнительной (априорной) информации о реконструируемой функции и включение найденной информации в методы реконструкции распределений плотности.
• Разработка дискретных алгоритмов и программ для реконструкции распределений плотности и их оптимизация, направленная на повышение скорости реконструкции и снижение объема требуемой машинной памяти.
• Нахождение количественных оценок результатов реконструкции и сравнение полученных усовершенствованных методов реконструкции распределений плотности с существующими методами.
Научная новизна:
• Разработаны усовершенствованные устойчивые методы реконструкции распределений плотности в сечениях объектов на основе метода свертки и обратной проекции и алгебраического алгоритма восстановления.
• Разработан итерационный алгоритм определения длин пересечений прямой с элементами области сканирования.
• Разработан метод нахождения границ реконструируемого объекта по исходным проекционным данным.
• Алгебраический алгоритм восстановления приведен к форме, удобной для компьютерной реализации.
• Предложено упрощение стабилизирующего множителя для метода регуляризации Тихонова, не снижающее эффективность метода.
• Показано, что решение задачи реконструкции распределений плотности методом р-фильтрации обратной проекции совпадает с решением уравнения Радона, приведенного к стандартному виду - интегральному уравнению Фредгольма I рода.
• На основе усовершенствованных методов реконструкции распределений плотности разработано программное обеспечение, позволяющее проводить численные эксперименты в области КТ.
• Исследована зависимость соотношения устойчивости и разрешающей способности метода свертки и обратной проекции с использованием метода регуляризации Тихонова от значений параметра и порядка регуляризации.
• Исследована зависимость соотношения устойчивости и разрешающей способности алгебраического алгоритма восстановления от количества итераций и значения параметра релаксации.
Практическая ценность. Полученные в результате проделанной работы усовершенствованные устойчивые методы реконструкции позволяют реконструировать распределения плотности с меньшими погрешностями, чем существующие методы. На основе данных методов создано программное обеспечение, с помощью которого были проведены эксперименты по реконструкции распределений плотности по данным, полученным с рентгеновского компьютерного томографа РКТ-01, расположенного во ВНИИ технической физики (Российский Федеральный ядерный центр, г. Снежинск). Результаты экспериментов позволили сделать вывод о том, что результаты данной работы могут быть использованы при создании программного обеспечения рентгеновских компьютерных томографов.
Преобразование Радона как основа математического описания задачи
Рассмотрим задачу реконструкции распределений плотности в сечениях объектов на примере рентгеновского компьютерного томографа с параллельной схемой сканирования (см. рис 1.4, а также рис. 1.2) [9, 24-26].
На рис. 1.4 показано сечение D исследуемого объекта, характеризуемое плотностью вещества (коэффициентом поглощения рентгеновских лучей) f(x,y), где х,у - неподвижная относительно объекта система декартовых координат (с ней совмещена система координат х ,у , необходимая далее). Расположенные на раме рентгеновские трубки излучают остронаправленные пучки рентгеновских лучей интенсивности I0, которые, пройдя через вещество и испытав частичное поглощение, регистрируются соответствующими детекторами (приемниками). Такой эксперимент проводится для ряда значений угла поворота рамы ф є [0,7і). Введем вращающуюся (неподвижную относительно рамы) систему декартовых координат s,l. Тогда можно обозначить через I(s, ф) значения интенсивности излучения, измеренные детекторами.
Поскольку система декартовых координат s,l повернута относительно системы декартовых координат х,у на угол Ф, то, используя формулы преобразования прямоугольных декартовых координат при повороте осей, получаем связь между системами координат s,l и х,у [6]: х = scosq -l sirup, у = s simp+ 1 cosy. Поэтому мы можем записать уравнение (1.5) как [4, 24] од Rf(x,y) = \f(scos(p-lsin(p,ssmq) + lcos(p)dl = q(s,q ) . (1-6) -со Уравнение (1.5) справедливо в случае отсутствия погрешностей в результатах измерений. Однако на практике невозможно точно измерить значения интенсивности излучения I(s, ф), поскольку, во-первых, число квантов, испускаемых источником излучения за единицу времени является случайной величиной, и, во-вторых, поглощение квантов веществом тоже происходит случайным образом [1, 18, 22]. Таким образом, значения интенсивности I(s, ф) неизбежно будут измеряться с погрешностями, т.е. измеряться будут не точные значения I(s, ф), а некоторые за-шумленные значения 7( ,Ф) = /( ,Ф) + 5/, (1.7) где 6/ - случайная погрешность измерения. Учитывая (1.7), приходим к преобразованию Радона в виде [11] q(s, p) = Rf(x,y) + i}(s,q ), (1.8) где г(.у,ф) - шумовая составляющая проекционных данных. Заметим, что интегральное уравнение Радона (1.5) имеет нестандартный вид, а именно: в нем отсутствует ядро и интеграл является одномерным, хотя задача является двумерной. Для решения такого уравнения требуется применение специальных методов, изложенных в главе 2.
Вычисление преобразования Радона функции f(x,y), т.е. определение функции q(s,q ), является прямой задачей КТ, а обращение преобразования Радона, т.е. реконструкция функции f(x,y) по известной функции q(s,(p), - обратной задачей. При решении гораздо более сложной обратной задачи возникают следующие трудности. Во-первых, обратная задача КТ является некорректной [9, 16, 22], т.е. небольшие погрешности в проекционных данных q(s,q ) могут привести к большим погрешностям в реконструируемой функции f(x,y). Поскольку на практике функция q(s,ф) всегда содержит погрешности, то необходимо прилагать усилия, направленные на повышение устойчивости методов реконструкции. Во-вторых, для всех методов решения некорректных задач характерно противоречие между разрешающей способностью и устойчивостью [3, 16], а именно: увеличение разрешающей способности (уменьшение систематической погрешности метода реконструкции) вызывает понижение устойчивости (увеличение случайной погрешности метода реконструкции), и наоборот.
Существует два вида задач реконструкции распределений плотности: задача с полными данными и задача с неполными данными [9, 10, 45, 54-56, 62]. В задаче с полными данными предполагается, что проекционные данные q(s, p) известны для всех прямых L(s,y), проходящих через исследуемое сечение. В задаче с неполными данными предполагается, что известен только неполный набор проекционных данных: значения функции I(s, ф) вдоль некоторых направлений измерены весьма неточно или не измерены вовсе. Если в таких случаях попытаться применить стандартные алгоритмы, то реконструированные распределения плотности оказываются искаженными сильными погрешностями. Задачи с неполными данными могут быть следующих типов:
Задача с ограниченным диапазоном углов возникает в КТ, если значения интенсивности излучения I(s, ф) могут быть измерены только в некотором ограниченном угловом интервале [10, 11, 39, 45, 52, 62]. Эта задача может быть проиллюстрирована схемой багажного компьютерного томографа (рис. 1.5) [54], который состоит из неподвижного источника рентгеновского излучения и неподвижного набора детекторов. Исследуемый объект лежит на движущемся конвейере, и плоскость досмотра проходит через плоский веерный пучок рентгеновских лучей. На рисунке 1.5 справа показаны рентгеновские лучи, проходящие через фиксированную точку плоскости досмотра при прохождении ее сквозь веерный пучок. Углом расхождения веерного пучка определяется ограниченный угловой интервал, в котором известны проекционные данные. Проекционные данные для других направлений, например, параллельных ленте конвейера, не могут быть получены. детееторы излучения
Внешняя задача возникает в том случае, когда исследуемая плоскость (сечение) содержит область, непрозрачную для рентгеновских лучей (рис. 1.6) [40, 54]. На практике такая задача может возникнуть, если исследуемая плоскость содержит металлическую структуру, например, хирургический зажим [35] или искусственный бедренный сустав [42]; такая же задача возникает в том случае, когда исследуемое сечение проходит через костное образование, скажем, позвоночник [46].
В каждой проекции имеется набор известных линейных интегралов, соответствующих прямым, проходящим через объект мимо непрозрачной области, и отсутствует набор линейных интегралов, соответствующих прямым, которые проходят через непрозрачную область. Кроме того, значения интенсивности излучения I(s, ф), измеренные вдоль непрозрачных для рентгеновских лучей структур, содержат большие погрешности из-за увеличения жесткости излучения [18, 46] и небольших движений объекта во время исследования. При решении внешней задачи необходимо доопределить проекционные данные q(s,(p), соответствующие прямым L(s,(p), проходящим через непрозрачную область.
Внутренняя задача возникает в КТ, если сканирование всего сечения необязательно или нежелательно [9]. Если интерес вызывает лишь небольшой участок тела (исследуемая область), то не стоит подвергать воздействию рентгеновского излучения остальную часть тела. В этом случае значения интенсивности излучения I(s, р) измеряются только вдоль прямых, которые пересекают исследуемую область или проходят вблизи нее.
Формулы обращения преобразования Радона
Теорема Брейсуэлла. Сформулируем теорему Брейсуэлла или проекционную теорему о центральном сечении, которая утверждает, что одномерное преобразование Фурье проекции некоторой функции двух переменных на некоторое направление совпадает с сечением двумерного преобразования Фурье этой функции плоскостью, направленной под тем же углом [24, 31].
Подставив данное выражение в (2.9), получим формулу (2.8). Оказалось, что решение уравнения Фредгольма I рода (1.9) методом двумерного преобразования Фурье совпадает с решением уравнения (1.5) методом р-фильтрации обратной проекции (2.8). Таким образом, решение (2.9) уравнения Фредгольма I рода (1.9) реализует формулу обращения (2.7). Однако при дискретизации формулы (2.8) возникают следующие трудности [65]. Во-первых, для вычисления S 2Bq(s,y) необходимо знать Bq(s,q ) не только в области реконструкции, но и во всей области определения. И, во-вторых, функция Ф2дО,ф) имеет особенность в точке (0,0). Из-за этих трудностей метод р-фильтрации обратной проекции не получил широкого распространения. Метод преобразования Фурье. Метод преобразования Фурье основан на использовании проекционной теоремы о центральном сечении. Запишем (2.4) в виде: F(a cos ф, со sin ф) = -= Q9 (ю), л/2я где F(cocoscp,(osin(p) = Ф2/(х,у), а ф(со) = Oq9(s). Применяя к обеим частям уравнения обратное преобразование Фурье (1.17), получаем [30, 33, 59, 72, 73] f(x,y) = -—-ті Ш j{ \,artcg 2- L Vco, . Однако при дискретной реализации алгоритм Фурье дает значительные пофешности в реконструируемой функции, поскольку один из шагов алгоритма Фурье требует интерполяции при переходе от полярных координат к декартовым координатам [9, 30]. На этом этапе возникают существенные сложности, связанные с тем, что количество отсчетов данных на единицу площади с удалением от начала координат уменьшается, вследствие чего пофешность интерполяции при больших значениях частоты неизбежно растет.
Метод свертки и обратной проекции является дальнейшей эволюцией метода преобразования Фурье и позволяет полностью избежать интерполяции за счет перехода от декартовых координат к полярным координатам в пространстве частот. Данный метод был впервые применен в радиоастрономии Брейсуэллом и Риддли [32]. Его применение к медицинским задачам восходит к работам Шеппа и Логена [69], а также Ра-мачандрана и Лакшминараянана [63]. Метод свертки и обратной проекции получил крайне широкое распространение, поскольку он не требует больших вычислительных затрат и его достаточно просто реализовать [43].
Таким образом, дискретный алгоритм принимает следующий вид: 1. Производится фильтрация исходных проекционных данных qrp в частотной области. Для каждой проекции с номером pe[o,P-l] вычисляется дискретное преобразование Фурье (1.23)-(1.24), производится фильтрация и вычисляется обратное дискретное преобразование Фурье (1.25)-(1.26). 2. Производится обратное распространение фильтрованных проек ционных данных pgrp по плоскости сканирования в соответствии с одним из дискретных вариантов формулы (2.1). Оценка количества требуемых операций составляет 0(N2 Р) или 0(N Р R), не считая операций, необходимых для выполнения прямого и обратного преобразований Фурье.
Выражение (2.20) представляет собой одну итерацию ААВ, в результате которой мы переходим от приближения Fk к приближению Fk+l. Для выполнения одной итерации требуется O(R-P-N) операций. Длины / можно вычислить с помощью алгоритма определения длин пересечений прямой с элементами области сканирования.
Формула (2.20) не является удобной с точки зрения компьютерной реализации, поскольку и распределение плотности, и исходные проекционные данные удобнее хранить в виде двумерных, а не одномерных массивов. Для того чтобы получить соответствующий вид формулы (2.20), вернемся к старым индексам г,р и п,т. Следует отметить, что процесс итераций сходится только в случае отсутствия погрешностей в проекционных данных qrp. При наличии же погрешностей ААВ обеспечивает лишь полусходимость вследствие некорректности задачи: на некотором шаге алгоритм дает удовлетворительное решение, но затем расходится [9]. Поэтому определенную сложность представляет определение количества итераций. Обычно момент остановки итераций определяют по невязке или по поправке [2]. Основным преимуществом ААВ является его универсальность. Этим алгоритмом можно воспользоваться при любой схеме сканирования, а также в задачах с неполными данными. Оценка количества требуемых операций для рассмотренного алгоритма составляет O(R-P-N-K). 1. Решение уравнения Радона, приведенного к интегральному уравнению Фредгольма I рода (1.9), методом двумерного преобразования Фурье (2.9) совпадает с решение уравнения (1.5) методом р-фильтрации обратной проекции (2.8). Таким образом, решение (2.9) реализует формулу обращения (2.7). 2. Для обеспечения сходимости метода свертки и обратной проекции (2.11)-(2.12) необходимо ввести предельную частоту сошах =л/Л. В результате интеграл (2.12) можно взять аналитически в виде (2.14). 3. Предложенная модификация дискретного алгоритма свертки и обратной проекции, полученного из уравнения (2.11), позволяет увеличить скорость реконструкции распределений плотности. 4. Алгебраический алгоритм восстановления в форме (2.20) приведен к форме (2.21), более удобной с точки зрения компьютерной реализации. Применение алгоритма определения длин пересечений прямой с элементами области сканирования при численной реализации формулы (2.21) позволяет увеличить скорость реконструкции и снизить объем требуемой машинной памяти. 5. Поскольку алгебраический алгоритм восстановления обладает лишь полусходимостью, то необходимо выбрать критерий, позволяющий определять момент остановки итераций.
Использование метода регуляризации Тихонова для повышения устойчивости методов реконструкции
В предыдущем параграфе были рассмотрены общие принципы регуляризации. В данном параграфе мы получим модифицированные формулы обращения преобразования Радона, использующие метод регуляризации Тихонова, откуда будут следовать конкретные методы реконструкции распределений плотности, использующие данный метод регуляризации. Тем не менее, формула (3.7) достаточно сложна, поэтому более эффективный путь - вычислять интеграл (3.6) численно по какой-либо квадратурной формуле [6], например, по формуле трапеций или по формуле левых прямоугольников [16].
Методы реконструкции распределений плотности, использующие метод регуляризации Тихонова, подавляют высокие частоты Фурье более аккуратно, чем методы реконструкции, использующие усечение по частоте. С одной стороны, высокие частоты наиболее сильно реагируют на погрешности в исходных данных, и поэтому их нужно подавлять, а, с другой стороны, от высоких частот зависит разрешающая способность методов реконструкции, поэтому это подавление должно быть умеренным [16, 67, 68]. Таким образом, с помощью задания значений параметра и порядка регуляризации можно устанавливать соотношение между разрешающей способностью (т.е. систематической погрешностью) и устойчивостью (т.е. случайной погрешностью) методов реконструкции распределений плотности [67, 68].
В главе 1 подробно описаны различные типы задач с неполными данными. Задача реконструкции распределений плотности по неполным данным представляет собой существенную трудность, поскольку реконструируемая функция не определяется однозначно неполным набором проекционных данных [9, 54]. В этом случае применение стандартных алгоритмов реконструкции приводит к тому, что реконструированное распределение плотности содержит существенные погрешности.
Существуют приемы, позволяющие снизить влияние неполноты проекционных данных на погрешность реконструкции. Одним из таких приемов является дополнение неполных проекционных данных [9, 50]. Данный подход заключается в том, что по известным проекционным данным происходит экстраполяция проекционных данных в исключенную область. Если существует более чем одно решение задачи экстраполяции, то выбирают решение с минимальной нормой. Использование описанного приема приводит к исчезновению наиболее сильных погрешностей и появлению новых, сравнительно небольших, но все же существенных [51].
Другим способом решения задач с неполными данными является применение метода реконструкции (например, алгебраического алгоритма восстановления), не требующего дополнения проекционных данных [10, 11, 48, 54, 56]. Следует отметить, что применение таких методов не гарантирует удовлетворительное решение задачи реконструкции распределений плотности.
Наконец, можно снизить погрешность реконструкции с помощью включения в метод реконструкции распределений плотности какой-либо дополнительной (априорной) информации о решении [10, 11, 39, 54, 55]. Для задач с неполными данными использование дополнительной информации оказывается особенно полезно, поскольку реконструкция по неполным данным - это не просто задача обращения, а задача аппроксимации и оптимизации. Следует отметить, что использование дополнительной информации снижает погрешность реконструкции и при решении задачи с полными данными. При этом включать найденную дополнительную информацию можно как в стандартные методы реконструкции (например, в метод свертки и обратной проекции), так и в методы, не требующие дополнения исходных данных.
Отметим, что все описанные приемы можно использовать как отдельно друг от друга, так и совместно в любой комбинации. Мы не будем использовать дополнение исходных проекционных данных. Вместо этого, будем заменять отсутствующие проекционные данные некоторыми значениями на основе найденной дополнительной информации.
Будем представлять дополнительную (априорную) информацию в виде ограничений, которым должно удовлетворять реконструированное распределение плотности. Приведем некоторые обычные ограничивающие условия и возможные способы их реализации [54]: 1. Неотрицательность: отрицательные значения распределения плотности заменяются нулевыми. 2. Известные границы объекта: значения распределения плотности, соответствующие элементам области сканирования вне границ объекта, заменяются нулевыми. 3. Известная область объекта: в известной области задаются известные значения распределения плотности. 4. Известный интервал значений: значения распределения плотности, выходящие за верхнюю границу интервала, заменяются максимальным значением, значения, выходящие за нижнюю границу, - минимальным значением. 5. Известное среднее значение в области: все значения в этой области изменяются на одну и ту же константу так, чтобы полученное среднее значение удовлетворяло ограничению. Отметим, что приведенные способы реализации ограничений не являются единственно возможными. Так, например, требование неотрицательности может быть реализовано путем замены отрицательных значений абсолютными значениями или средней величиной окружающих положительных значений.
Несмотря на универсальность данного результата, мы не будем использовать формулу (3.9) из-за ее излишней сложности. Вместо этого будем включать дополнительную информацию непосредственно в численные алгоритмы любым удобным способом.
Остается открытым вопрос о том, как найти какую-либо дополнительную информацию. К сожалению, в литературе этот вопрос практически не рассматривается. Зачастую предлагается предварительно реконструировать функцию каким-либо методом реконструкции, по полученным данным определить дополнительную информацию и включить ее в вычисления на следующем шаге [54]. Однако, как мы покажем в следующем параграфе, дополнительную информацию можно извлечь непосредственно из исходных проекционных данных.
1. Применение методов регуляризации позволяет повысить устойчивость методов реконструкции распределений плотности. При этом для итерационных методов реконструкции выбор параметра регуляризации сводится к выбору количества итераций.
2. Предложенное упрощение стабилизирующего множителя для метода регуляризации Тихонова не снижает эффективность метода.
3. С помощью задания значений параметра и порядка регуляризации для метода регуляризации Тихонова можно устанавливать требуемое соотношение между разрешающей способностью и устойчивостью методов реконструкции.
4. Включение дополнительной информации о решении в методы реконструкции при решении задач с неполными данными позволяет снизить погрешность реконструкции распределений плотности.
5. Хотя существует универсальное решение (3.9) уравнения (3.1) с учетом уравнения ограничений (3.8), данная формула излишне сложна, поэтому предпочтительнее включать найденную дополнительную информацию непосредственно в численные алгоритмы любым удобным способом.
6. Предложенный метод определения границ исследуемого объекта по исходным проекционным данным позволяет снизить погрешность и увеличить скорость реконструкции распределений плотности.
Погрешность реконструкции при использовании метода регуляризации Тихонова
Изучим вопрос о зависимости погрешности реконструкции при использовании метода регуляризации Тихонова от параметра и порядка регуляризации. Для реконструкции распределений плотности будем использовать метод свертки и обратной проекции в виде (3.5)-(3.6) и включать в него информацию о границе объекта (3.11), а также требование неотрицательности реконструируемого распределения плотности.
Решение в таком виде имеет следующие преимущества. Во-первых, вычислять значения плотности требуется не для всей области сканирования, а только для подмножества точек, не выходящих за границу объекта, что увеличивает скорость реконструкции. Во-вторых, для всех остальных точек значения плотности известны точно, что снижает погрешность реконструкции.
В таблице 4.1 приведены значения относительного СКО и относительной гладкости реконструированных распределений плотности в случае измерений без погрешностей, т.е. без шума в проекционных данных. Как видно, использование регуляризации в этом случае не сильно снижает погрешность реконструкции [28]. Так, максимальное снижение погрешности реконструкции составляет менее 1%. В данном случае единственное преимущество регуляризации состоит в том, что она по 64 зволяет получить распределение плотности с большей относительной гладкостью. В таблице 4.2 приведены значения относительного СКО и относительной гладкости, полученные при реконструкции распределений плотности при наличии в проекционных данных 5% шума с равномерным распределением. В этом случае применение метода регуляризации Тихонова оказывается полезным.
Еще одно преимущество регуляризации заключается в том, что, посредством выбора соответствующих значений параметра и порядка регуляризации, можно добиться требуемого соотношения между устойчивостью и разрешающей способностью метода реконструкции. Например, можно повысить значение относительной гладкости без существенного увеличения погрешности реконструкции (см. рис. 4.5). В таблице 4.3 приведены значения относительного СКО и относительной гладкости, полученные при реконструкции распределений плотности при наличии в проекционных данных 10% шума с равномерным распределением. Здесь мы наблюдаем похожую картину: регуляризация помогает уменьшить относительное СКО и увеличить относительную гладкость функции в одно и то же время, при этом разница между решением без регуляризации и с регуляризацией отличаются еще больше, чем в случае с 5% шумом. Решение без регуляризации дает погрешность около 25%, а применение регуляризации уменьшает эту погрешность почти вдвое.
Таким образом, применение метода регуляризации Тихонова позволяет повысить устойчивость методов реконструкции распределений плотности, а посредством выбора значений параметра и порядка регуляризации можно варьировать соотношение между устойчивостью и разрешающей способностью методов реконструкции распределений плотности. Таблица 4. Реконструкция с помощью метода свертки и обратной проекции по проекционным данным с 10% шумом
Графики зависимости значений относительной гладкости от параметра регуляризации для реконструкции с 5% шумом. Сравнение данных в таблицах 4.1-4.3 показывает, что относительная гладкость растет с увеличением значения параметра регуляризации (см. рис. 4.7), а относительное СКО убывает, пока не достигнет минимума, а затем начинает расти (см. рис. 4.6). Поэтому параметр регуляризации следует выбирать аккуратно. Анализируя таблицы 4.1-4.3 можно заключить, что:
1. С увеличением уровня шума увеличивается значение параметра регуляризации, соответствующее минимальной погрешности реконструкции.
2. Чем больше значение порядка регуляризации, тем меньше значение параметра регуляризации, соответствующее минимальному значению относительного СКО (см. рис. 4.6).
3. С увеличением уровня шума увеличивается минимально достижимое значение относительного СКО. Компьютерные эксперименты показывают, что минимально возможное значение относительного СКО достигается при р = 2 (см. рис. 4.6), а оптимальное соотношение между относительным СКО и относительной гладкостью - при р = 1 (см. рис. 4.6 и 4.7). В таблице 4.5 приведены значения относительного СКО и относительной гладкости, полученные при реконструкции распределений плотности при помощи ААВ по проекционным данным без погрешностей с использованием дополнительной информации. Как видно, использование дополнительной информации существенно увеличивает скорость сходимости ААВ. Например, при реконструкции без релаксации (Х = 1) и без использования дополнительной информации на 100-ой итерации погрешность реконструкции составляет 9%, а с использованием дополнительной информации — менее 1%, причем в этом случае значение относительной гладкости реконструированного распределения плотности равно значению относительной гладкости исходного распределения плотности. Кроме того, при использовании дополнительной информации убывание значения относительного СКО происходит монотонно при любых значениях параметра релаксации. Также следует отметить, что при больших значениях параметра релаксации наблюдается монотонный рост значения относительной гладкости реконструированного распределения плотности. В отличие от реконструкции по проекционным данным без погрешности, при реконструкции по зашумленным данным более предпочтительным оказывается выбор малого значения параметра релаксации. Это объясняется тем, что при малом значении параметра релаксации быстро реконструируется среднее значение функции, и, таким образом, к моменту начала расхождения итераций мы имеем минимальную погрешность реконструкции. Однако выбор момента остановки итераций, исходя из минимума относительного СКО, не всегда будет оптимальным. Например, рассмотрим рисунки 4.8 и 4.9. Рисунок 4.8 соответствует реконструкции с параметром релаксации X = 0.05 на 5-ой итерации. В этом случае мы получаем минимальное значение относительного СКО, но и большое значение относительной гладкости. Для более точной реконструкции границ областей предпочтительнее произвести реконструкцию с тем же параметром релаксации до 10-ой итерации (рис. 4.9). Несмотря на то, что в этом случае значение относительного СКО больше, чем на 5-ой итерации, значение относительной гладкости меньше, и, таким образом, выше разрешение. Отсюда можно сделать вывод, что момент остановки итераций предпочтительнее определять методом подбора, поскольку в этом случае можно добиться требуемого соотношения между разрешающей способностью и устойчивостью ААВ. Технически это можно реализовать следующим образом: выводить на дисплей изображение реконструированного распределения плотности после завершения каждой итерации, сохранять результат реконструкции в памяти и дать оператору возможность определять момент остановки итераций и просматривать изображения реконструированных распределений плотности, полученные на предыдущих шагах с тем, чтобы выбрать один или несколько наиболее подходящих результатов.
