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

РОССИЙСКАЯ АКАДЕМИЯ НАУК СИБИРСКОЕ ОТДЕЛЕНИЕ ВЫЧИСЛИТЕЛЬНЫЙ ЦЕНТР

УДК 518.517.948+533.9+539.143.43+615.47

МЕТОДЫ И АЛГОРИТМЫ ТОМОГРАФИИ ГАЗА И ПЛАЗМЫ

Диссертация на соискание ученой степени доктора физико-математических наук в форме научного доклада

05.13.16 - применение вычислительной техники, математического моделирования и математических методов в научных исследованиях

РГ6 о/}

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

2 О

■1

Пи колов Валерий Владимирович

Новосибирск - 1996

Работа выполнена в Институте теоретической и прикладной механики СО РАН (г.Новосибирск)

Официальные оппоненты: член-корреспондент РАН,

доктор физико-математических наук, профессор Э.П.Кругляков; доктор технических наук, профессор В.П.Пяткин^ доктор физико-математических наук, профессор В.Е.Куницын

Ведущая организация: Институт математики СО РАН

Защита состоится "28 " января 1997 года в 1<£ час. на заседании Специализированного совета Д 002.10.02 при Вычислительном центре Сибирского отделения РАН по адресу: 630090, Новосибирск, проспект академика Лаврентьева, С.

С докладом можно ознакомиться в читальном зале библиотеки ВЦ СО РАН (Новосибирск, пр. акад. Лаврентьева, 6) Автореферат разослан " /1996 г.

Ученый секретарь Специализированного совета к.т.н.

Г.И.Забиняко

1 Общая характеристика работы

1.1 Актуальность работы

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

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

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

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

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

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

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

1.2 Цель работы

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

1.3 Научная новизна

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

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

• развит модифицированный итерационный алгоритм класса Гершберга-Папулиса для задач двумерной и трехмерной томографии и найден параметр нелинейного управления качеством реконструкции в условиях малого числа ракурсов;

• создан новый итерационный алгоритм по гладкому восполнению двумерной и трехмерной функций, заданных на нерегулярной сетке, и на этой основе построен итерационный метод восполнения двумерных синограмм;

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

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

1.4 На защиту выносятся:

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

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

3. Новый алгоритм восполнения синограмм и результаты его применения для томографической диагностики в видимом диапазоне плазмы токамака RTP (Ньювэхайн, Нидерланды).

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

5. Разработанное программное обеспечение в виде пакета прикладных программ библиотечного типа по малоракурсной томографии газа и плазмы TOPAS.

1.5 Практическая ценность работы

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

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

- разработан алгоритм фильтрации и обратного проецирования для решения нелинейной задачи двумерной томографии в условиях слабой рефракции;

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

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

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

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

- разработан метод электронно-пучковой томографии гиперзвукового обтекания моделей разреженным потоком газа.

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

Результаты работы внедрены и используются в различных организациях, в том числе, в Физическом институте РАН (Москва), Московском институте электронной техники (Зеленоград), Институте физики плазмы (Ньювэхайн, Голландия), Эйндховенском технологическом университете (Эйндховен, Голландия), Институте физики плазмы (Юлих, Германия), Институте физики низкотемпературной плазмы (Грейфсвальд, Германия) и др.

Одномерный блок (ELLIPS) пакета томографии газа и плазмы TOPAS сдан в 1979 г. в Госфонд алгоритмов и программ СССР (per. номер П004148, см.[11]).

Имеются акты о внедрении пакета прикладных программ TOPAS и его отдельных частей в Институте общей физики (Москва), Физическом институте им.П.Н.Лебедева (Москва), Институте оптики атмосферы СО РАН (Томск), Московском авиационном институте (Москва), предприятии п/я Р-6324, Всесоюзном научно-исследовательском институте оптико - физических измерений (Москва), Институте физики АН Литвы (Вильнюс).

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

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

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

Основные результаты работы изложены в 4-х монографиях (с соавторами).

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

Всесоюзных семинарах по комплексам программ математической физики (Новосибирск, 1982; Ростов-на-Дону; 1990); Международной школе-семинаре по методам оптимизации и их приложениям (Иркутск, 1989); Всесоюзных и международных конференциях по некорректно поставленным задачам (Фрунзе, 1979; Москва, 1990); Всесоюзных и международных конференциях по методам аэрофизических исследований (Новосибирск, 1973, 1976, 1982, 1994, 1996); Всесоюзных и международном симпозиумах по вычислительной томографии (Новосибирск, 1983; Куйбышев, 1985; Киев, 1987; Ташкент, 1989; Звенигород, 1991; Новосибирск, 1993); Всесоюзных совещаниях по диагностике высокотемпературной плазмы (Харьков, 1986; Минск, 1990); Всесоюзных конференциях по генераторам низкотемпературной плазмы (Фрунзе, 1974, 1983; Алма-Ата, 1977; Новосибирск, 1980, 1989); Всесоюзных конференциях "Обработка изображений и дистанционные исследования" (Новосибирск, 1987, 1990); Международных конференциях по явлениям в ионизованных газах (Berlin, 1977; Minsk, 1981; Düsseldorf, 1983; Budapest, 1985); Международных конференциях по управляемому синтезу и физике плазмы (Innsbruck, 1992; Lisboa, 1993; Montpellier, 1994); Ежегодных европейских конференциях по исследованиям атмосферы (Iíiruna, Sweden, 1992; 1993); Международном совещании по реконструкции полностью трехмерных изображений в радиологии и ядерной медицине (Aix-les-Bains, 1995);

семинарах Института ядерной физики СО РАН, Института математики СО РАН, Международного томографического центра (Новосибирск), Института теоретической и прикладной механики СО РАН, Физического института РАН (Москва), Института атомной энергии (Москва), Института прикладной геофизики (Москва), Всесоюзного научно-исследовательского института оптико-физических измерений (Москва), Московского института электронной техники (Зеленоград), Института оптики атмосферы СО РАН (Томск), Института физики плазмы (Ньюв-эхайн, Голландия), Института физики плазмы (Юлих, Германия), Института физики низкотемпературной плазмы (Грейфсвальд, Германия), Центральной электронной лаборатории (Юлих, Германия), Энндховенского технологического университета (Эйндховен, Голландия), Объединенного европейского токамака (Абингдон, Великобритания).

1.7 Публикации

Основные материалы научного доклада опубликованы в 67 печатных работах, в том числе в 4 монографиях.

2 Основное содержание работы 2.1 Введение

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

Запишем уравнение переноса излучения вдоль прямой L:

MOW (1)

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

/с» /оо

Ле(1)ехр[-/ к(1')М'). (2)

■оэ JI

Интегральное излучение плазмы (проекция) регистрируется вдоль

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

В приближении оптически тонкой плазмы (к = 0) регистрируемое излучение описывается классическим интегральным преобразованием Радона (Natterer, 1990) (далее неизвестные локальные коэффициенты эмиссии будем обозначать через д(х,у), а регистрируемое интегральное излучение - через /(я,4))=

/со

g(.x,y)dl. (3)

-оо

Измерив в диагностических экспериментах интегральное излучение для набора углов {&},£ € [0,180°) и набора координат {р;}, можно получить оценку искомого решения :

5=я;1(/ + ч). (4)

Здесь R'1 - приближение к обратному преобразованию Радона, ij - шумовая составляющая проекционных данных.

Если пренебречь шумом, то для решения уравнения (3) можно получить формулу, известную как обратное преобразование Радона:

д(х,у) =--— Г Г НрЛНр (

2тг2 Jo J-oo (р — р0)2 К )

где pa — — isini-f- усos£ = rsin(y — i), а (г,у) - полярные координаты точки (х,у).

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

- малость числа направлений (ракурсов) наблюдений К;

- малость числа отсчетов ( по р ) TV;

- угловые ограничения в регистрации проекций;

- возможное рассогласование проекций;

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

- наличие шумов с различной статистикой в измерительных каналах;

- появление в некоторых случаях непрозрачных включений внутри плазменных объектов;

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

Методы получения проекционных данных в ПТ обычно являются дальнейшим развитием классических способов получения диагностической информации о плазме. В основном, это эмиссионные (различного спектрального диапазона - от инфракрасного до мягкого рентгеновского), трансмиссионные: интерферометрнческие (также в разных областях электромагнитных волн), дефлектомегрические (спекло- шлирен- и теневые методы), корпускулярные и т.п. методы (Кузнецов, Щеглов; Лохте-Хольтгревен, Подгорный, Пятницкий, Хаддлстоун). Все они, конечно же, получили свое развитие, прежде всего, для одноракурсных систем измерений. При определенных ограничениях на класс изучаемых плазменных образований даже такие системы сбора данных можно отнести к томографическим, если они позволяют по измеренной информации определить локальные характеристики плазмы в некоторых ее сечениях, а иногда - и во всем объеме.

2.2 Одномерные алгоритмы томографии

Случаи редукции двумерных задач томографии к одномерным

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

Такое сведение зачастую возможно при введении предположений о той или иной степени симметрии объекта: круговой, эллиптической, наличия плоскостей симметрии и т.д.

Рассмотрим, что даст при решении задач томографии априорная информация о характеристиках изолиний искомой функции.

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

Более сложные геометрии, допускающие все же редукцию двумерного интегрального уравнения к одномерному - это случай изолиний, заданных в виде смещенных эллипсов, треугольных изолиний, семейства вложенных парабол высших порядков, и в весьма общем случае - выпуклых изолиний без самопересечений, допускающих введение обобщенного радиуса [2, 8, 9, 12, 24].

2.2.1 Одноракурсная томография: Инверсия Абеля

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

Я{г) =ffo(| г |) =5о(г).

Тогда все проекции /(р) становятся эквивалентными и независящими от угловой переменной

/(р) = Л(1Р1)н/„(Р).

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

/<*> ГСО ГОп(г)

да(г)с1у = 2/ -^ЩЩйт. (6)

■оо Jc у/Г* — X

Полученое уравнение относительно до(т) является интегральным уравнением типа Вольтерра 1-го рода, в данном частном случае имеющим название уравнения Абеля.

Известно, что данное уравнение имеет аналитическое обращение:

Яо(г) = _1 г=-1-4- г лщ^ „х. (7)

7Г Л у/X2 — 7*2 1ХГ аг ¿г уж2 — г2

Отсюда легко заметить природу некоректности задачи обращения интегрального уравнения Абеля: необходимо дифференцировать зашумленные экспериментальные измерения, а также преодолевать сингулярность в интеграле в его нижнем пределе. В силу того, что уравнение Абеля встречается как первое приближение к описанию физических явлений самой разнообразной природы, литература по методам его инверсии очень обширна (см. обзоры в [2, 4]). В целом, эти методы можно разбить на два больших класса. Первый класс - это представление искомой функции в виде разложения в ряд по каким либо базисным функциям с последующим решением получившейся системы линейных алгебраических уравнений (СЛАУ) относительно коэффициентов такого разложения. В частном случае коэффициентами разложения могут быть и сами значения искомой функции {д,} в узлах сетки по радиусу:

= (8)

э

Ко второму классу алгоритмов относятся различные численные реализации использования формул инверсии (7).

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

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

д = (А+\УЛ +аП)-1 A+Wf, (9)

где Д+ - транспонированная матрица А, а диагональная матрица IV состоит из обратных дисперсий с~2, оцениваемых для вектора измерений в каждой точке. Матрица 12 вносит в задачу априорную информацию о степени гладкости искомого решения, и определяется из аппроксимации его интегральной нормы.

Наиболее сложная проблема при использовании решения по (9) - это определение параметра регуляризации а. В методе невязки этот параметр определяется из согласования нормы невязки с уровнем шумов (Морозов, 1967, 1970):

М-Ада\\ (10)

где <52 - это норма экспериментальных шумов (но только случайных, без систематических). Обычно этот метод дает несколько переглаженное решение д.

В методе статистической регуляризации имеется несколько подходов к выбору этого параметра (Воскобойников, 1984; Преображенский, Пикапов, 1982). Неплохо зарекомендовал себя итерационный метод:

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

a = N [5р{П(Л+и'Л + аПГ,} + (а»,П<0]~1. (12)

Метод статистической регуляризации позволяет получить и теоретическую оценку ошибки восстановленного радиального распределения в виде его дисперсии в каждой точке:

<г2 = (А+ XV А + «П):'. (13)

2.2.2 Эллиптические изолинии

Естественным обобщением задачи реконструкции объекта по изолиниям является переход от системы линий уровня, аппроксимируемых окружностями, к системе эллиптических изолиний. Изолинии в виде софокусных или смещенных один относительно другого эллипсов характерны для ряда режимов работы стеллараторов и токамаков (Silver, Roney, 1980), винтовых дуг и т.п. Кроме того, как отмечают авторы работы (Кузнецов, Щеглов, 1974), нередко целесообразным бывает наклонное микроволновое зондирование столба плазмы, позволяющее повысить точность восстановления 2D распределения электронной плотности плазмы.

Рис. 1: Система эллиптических изолиний со смещениями, s - угол веерного сканирования.

В этом случае, если столб обладает осевой симметрией, эллиптические изолинии появляются автоматически. Рассмотрим систему изолиний в виде смещенных эллипсов в условиях веерного сканирования излучения оптически тонкого объекта. Такая задача решалась нами в [9]. В более общем случае уравнение изолиний запишется в виде (см. рис.1):

[х-А(р)\2/а* + [у-В(р)Г/Ь2 (14)

где р имеет смысл обобщенной радиальной переменной, пробегающей значения от р = 1 (внешний эллипс) до р = О (точка сгущения (а,/?)); а и 6 - полуоси эллипса.

Уравнение (14) записано в систем*; координат А'ОК, повернутой на угол ^ относительно лабораторной системы отсчета; <1 - сдвиг детектора Г) в системе Х'ОУ.

Примем для простоты, что центры эллипсов смещаются закону:

Л(р) = <*(1 - p.), Bip) = /3(1 - р).

по линеиному

(15)

После ряда выкладок, схема которых тложена в [8, 9], приходим к следующей связи проекционных данных / (£) с восстанавливаемым распределением д(р):

/(«= f K(l,t)u{p)dp, (16)

где

л'(/>,£) =

ш2 =

Р =

Л{ =

Pi =

Р2 =

/>(ш2 - р?) +р(р — RçsinQ /Я)~{р -р1){р + р'г,)Р

2 аЪ т2

a2cos2 (ip -i- £) +Ьг sin2 (у +£), a cos (ip + Ç) — l) sin (ip + £), Ло +dci<7i,

\p- RfSin îl(»7i

|p-^.sb^Kmip)-1.

(17)

В выражениях для р1 и р2 верхние знаки берутся для луч-сумм, проходящих правее точки сгущения; нижние знаки — в противоположном случае. Результаты работы [9] соответствуют отсутствию сдвига внешнего эллипса, т.е. <I — 0; расчетам Сильвера и Рони отвечает еще более частный случ.-.м отсутствия внутрених смещений эллипсов, а = ¡3 = 0.

l'J

Переход к случаю параллельных лучей наблюдения производится заменой: В.(я{гг £ —► х, 4- £) —* где О — угол поворота внешнего эллипса относительно лабораторной системы координат.

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

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

Как показано автором в [2] хорошо поддается уточнению, например, положение точки сгущения.

2.2.3 Изолинии произвольной формы

В работах автора [8, 9, 12, 24] рассматривались различные возможности более общего задания вида изолиний и соответствующих постановок обратных задач локальной диагностики газовых и плазменных объектов, включая проблему редукции задачи для оптически плотной плазмы к оптически тонкому слою.

Если ограничиться классом замкнутых выпуклых изолиний без самопересечений Ф(х,у,р) = 0, то дифференциал пути <1з вдоль луча зрения, направленного под углом £ к оси ординат при веерном сканировании (см. рис.1), можно представить в виде

¿х _ Г ¿р, <1х/<1р> О,

dS~sini ~ I Q-{p,i)dp, dx/dp<0. (18)

Здесь Q+ и Q~ - функции, определенные зависимостью х = х(р), а параметр р однозначно характеризует заданную систему изолиний Ф = О и тем самым - восстанавливаемое распределение д(р). Допустим, что функция a;(t) двузначная и выполняются соотношения

Q+ = q + Q, Q-=q-Q. (19)

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

1(() = Ц е(х, у) ехр [- jf ' к(х, у) ds'] ds, (20)

где, как и прежде, с(х,у) - локальные коэффициенты эмиссии плазмы, а к(х, у) - ее локальные коэффициенты абсорбции. С учетом (18) выражение

(20) принимает вид

Щ) = ехр (- £ Я+н ар) р е [<?+ ехр ¿р')

— ехр ( / (?~кс(р')] (1р,

ехр

Р1

Р1

где Р1(£) есть решение уравнения

(22)

В частности, для д = 0 получается более простое соотношение

(23)

которое является обобщением уравнения Фримена - Каца. Оптически тонкому слою (к = 0) соответствует интегральное уравнение

В частном случае эллиптических изолиний из (24) сразу же следует (16). Аналогично можно рассмотреть иные конфигурации изолиний.

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

2.3 Двумерные задачи томографии

Ранее уже упоминалась формула обратного преобразования Радона, на основе которой легко получить так называемые сверточные алгоритмы КТ (Хермен, 1983). Эти алгоритмы дают неплохие результаты в медицинской томографии, где число направлений наблюдений обычно весьма велико. Однако в задачах плазменной томографии для обычно малого числа ракурсов особое значение преобретают итерационные методы, позволяющие получать приближенные оценки томограмм путем последовательного введения в алгоритм априорной физической информации.

Рассмотрим некоторые из переспективных алгоритмов сверточного и итерационного типов.

2.3.1 Сверточные алгоритмы

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

(24)

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

з(х,у) =

— Г dí Г dp 2тг2 Jo J-™(p~po)2

2ж Jo J-c

Р о

(Р ~ Ро) —х sin Í + у cos£.

(25)

(26)

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

2.3.2 Алгоритм ГИСЭвг

Устойчивость численным реализациям приведенных формул повышает их регуляризация в смысле регуляризации сингулярных интегралов:

9(х,У) = -г^г Г d£ Р^[/(р04-р,£)+/(ро-р,О' 2тг2 Jo Jo р2

1 [" roo dp

--/ di -f-

2тг Jo Jo p

df(Po + P,Z) S/(p0-p,0

dp

dp

■2/(p„,€)](27)

(28)

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

Алгоритм RICSS2 (Radon Inversion with Cubic Spline Smoothing) подробно описан в работах [3, 22]. Он основан на использовании регуляризо-ванного варианта аналитического преобразования Радона (28). Внутреннее интегрирование с обобщенной функцией 1 /р ведется здесь с обходом особенности в нуле. Радиус обхода t¡ подбирается в численном эксперименте. Оценим погрешность такого обхода для заданной величины r¡.

Запишем выражение (28) в виде:

Ir» /-»о dp д(х,у) = -—- / di i — 2тг Jo Jn p

где

n 2;rA Jo

0/(PO+P,O 3/(PO-P,Í)

dp

Op

+ Д, = g„ + Д„, (29)

" dp га/(ро + р,0 _ ¿>/(po -p,0 p [ dp dp

(30)

MM

Если функция на всем интервале исследования удовлетворяет

условию Гелдера с показателем d:

dñPut) e/(Pj,í)

Sp

9р

<CiÍpi -píI', de (o,i),

(31)

то получаем следующую оценку для погрешности Д, замены функции д(х,у) на д„(х,у):

|Д,1

0/(Ро+Р,Г) ö/(p0 -р,£")

Jo р [ ri dp . 2d J

< С, —(2pY = -¡C*V •

Ja p d

dp

(32)

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

Оценка дисперсии решения ад алгоритмом ШС882 получена в [22]:

81тг2(М - 1)

8 R2

(33)

[(N-1)*

Здесь т*„ - дисперсии первой и второй производных сигнала (по р) в точке ра, М - число ракурсов, N - число отсчетов для одного ракурса и величина ¿2:

! J 0V-3)/2 l (A'-J)/2 г -4- _ у — +16 у -

(N - 1У 4 ¿i i2 ¿i (2i - l)2

(34)

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

2.3.3 Итерационный алгоритм Гершберга-Папулиса

Одним из весьма эффективных итерационных алгоритмов для параллельной схемы регистрации проекций является алгоритм Гершберга - Папули-са (Gerchberg-Papoulis, далее - алгоритм GP) [45], основанный на использовании поочередного исправления очередной оценки томограммы, а затем ее фурье-спектра. В работах (Défrise, De Mol, 1983; Tarn, 1983) этот метод изучался для задачи с ограниченным угловым диапазоном, £ £ [£„„„,£,„„]. Нами [45] была разработана регуляризованная версия алгоритма MGP для задач малоракурсной томографии и приведены результаты исследования его свойств для трех ракурсов.

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

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

Упрощенно итерационные шаги регуляризовапного алгоритма MGP можно представить так [45].

Шаг 1. По известному набору проекций получается набор их одномерных фурье-образов (с помощью быстрого преобразования Фурье (БПФ)). По теореме о центральном сечении эти величины дают в полярной системе координат в частотной плоскости значения фурье-образов томограммы на (К — 1) лучах д(£ + к/2,и). На этом этапе значения спектральных амплитуд зануляются вне данных лучей.

Шаг 2. Выполняется обратное двумерное БПФ для получения оценки томограммы д(я, у).

Шаг 3. Вносится априорная информация о положительности решения и его пространственной ограниченности (кругом радиуса R).

Шаг 4. Выполняется прямое двумерное БПФ от оценки томограммы д(х,у) после шага 3. Значения спектра на лучах, определенных на шаге 1, заменяются на величины, определенные там же непосредственно по имеющимся проекциям. Вносится также априорная информация о спектре объекта, проводится его сглаживание, согласованное с уровнем шумов в измерениях.

Шаг 5. Проверяются критерии окончания итерационного процесса. Если они не выполняются, то осуществляется переход на шаг 2..

Применение здесь же учета положительности томограммы на основе нелинейного обратного проецирования [45] позволяет ускорить, в ряде случаев, процесс сходимости итераций и повысить точность конечного результата. Кроме того, на шаге 2 необходимо осуществлять интерполяцию в спектральной области, с переходом от полярной к декартовой системе координат. Для выбранного круга частот |i/| < vx применяется билинейная интерполяция, а для частот [ь»| ~> i>\ - интерполяция "ближайшего соседа". Последняя выполняется лишь для точек, попадающих в некоторую Ь - окрестность луча (£ -)- 7т/2,и). Постепенное уменьшение этой полосы "влияния" с ростом номера итерации также дает возможность ускорения итерационного процесса.

На рис.2, 3 проиллюстрирована работа этого алгоритма в вычислительном эксперименте [45].

На графиках рис.2 приведены восстановленные томограммы алгоритмом MGP (Modified Gerchberg - Papoulis algorithm) двух моделей. Первая модель - рис.2,а - две гауссианы, центры которых помещены в точки (—0.25,0) и (0.25,0); вторая модель (рис. 2,Ь) - те же гауссианы, но с центрами в точках (0,-0.25) и (0,0.25). Параметры N = 33, К = 4, число

1.00 0.30 0.00 -0.50 -1.00

9

ш

0.53 0.00 -0.50 -1.00

-1.00 -0.50 С.80 0.50 1.00

.871,>00

..000,-18

-1.00 -0.50 О.ОВ 0.50 1.00

Ри(. 2: Результаты воссз ановления алгоритмом МСР двух фантомов: "горизонтального" (а), и "вертикального" (Ь).

итераций £ = 60. Погрешность Д1 = 13%(а) и 29%(Ь). На рис.3,а показан ход ошибок решения и зависимости от номера итерации Ь. Кривая 1 соответствует "горизонтальной" модели рис.2,а, а кривая 2 - "вертикальной" модели рис. 2,Ь.

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

2.3.4 Критерий информативности проекций

Разница в качестве восстановления "горизонтальной" и "вертикальной" моделей (рис. 2), а также отличия между кривыми 1 и 2 рнс.З требуют отдельного объяснения. Прежде всего, отметим разницу между проекциями этих моделей в системе трехракурсной томографии. В обоих случаях направления наблюдения располагались под углами £ = 0°,60° и 120° относительно оси ОХ. "Горизонтальная" модель отличается в данном случае от "вертикальной" тем, что для нее в число зарегистрированных проекций

50

30

20

b

0

30 L 60

0 10 И 50 С

1'цс. Чааисимость погрешности реконструкции "горизонтальной" (1) и "пертикяльной"(2) моделей от номера итераций К — 4 (а) и 7 (Ь).

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

В работе (И.Г. Казанцев, 198В) введен критерий информативности проекции (для томограммы, заданной внутри единичного круга):

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

2.3.5 Алгоритм синтеза проекций

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

(35)

в наших работах [55, 57, 58, 52]. В работах (Prince, Willsky, 1990) также исследовались близкие методы, в частности, вводились ограничения в виде моментных законов сохранения в пространстве синограммы. Однако в них при этом рассматривались лишь параллельные проекции, и неполнота проекционных данных вводилась лишь по угловой переменной, т.е. по прицельному параметру эти проекции считались заданными полностью. Ниже будут расмотрены условия совместности, которым должны удовлетворять синограммы следуя основным свойствам преобразования Радона.

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

Изложим основные принципы метода синтеза проекций (МСП) на конкретном примере томографической системы регистрации излучения видимого диапазона (ViLT), работающей на голландском токамаке RTP [55].

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

аП '*"" ° " b "

Рис. 4:

a) "Идеальная" схема расположения детекторов в пространстве синограммы для веерной и параллельной схем измерений. Точка схождения для веерного пучка находится на расстоянии 3R от начала координат.

b) Расположение детекторов для реальной экспериментальной системы ViLT на токамаке RTP [57].

На рис.4,а приведено положение детекторов (в координатах (£,р)) для "идеального" случая 5 регистрирующих камер с веерным обзором, с 17 детекторами каждая (показаны ромбиками), и равномерным расположением этих камер в диапазоне углов от 0° до 180°. Вершина каждого веера начинается на расстоянии D = 3R от начала системы координат (R -радиус томограммы). Крестиками показаны положения детекторов для параллельной системы регистрации (11 камер, по 27 детекторов). Далее условно систему регистрации для рис.4,а будем называть идеальной, а

систему с рис.4,Ь - реальной (речь идет, разумеется, лишь об описании системы томографической регистрации ЛЧЬТ; для других установок "реальным" будет иное расположение детекторов на плоскости синограммы, например, возможно и расположение аналогичное рис. 4,а).

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

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

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

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

Процедуру интерполяции запишем как обратную задачу, а именно, интерполяции на систему реальных измерений (£"',р™) значений функции, заданной в декартовой системе координат р*). Этой процедуре соответствует алгебраическое уравнение:

г = Е (зб)

Здесь вектор /г - это оценка синограммы в тех ее точках (в общем случае расположенных нерегулярно), где на плоскости синограммы расположены детектора, а - это значения синограммы в соседних узлах регулярной декартовой сетки, окружающих измерительную точку (£(т,р]"). Для билинейной интерполяции достаточно взять всего 4 узла, для интерполяции с более высокой точностью число членов в этом уравнении может быть больше.

Если повторить подобную процедуру для всех измерительных точек в плоскости томограммы, то мы получим систему линейных алгебраических уравнений (СЛАУ):

А? = Г, (37)

где в матрице А собраны все элементы ат, влияющие на вычисление интерполированного значения f.

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

Предлагается решать систему (37) регуляризованным вариантом метода Качмажа, который широко известен в компьютерной томографии как метод ART. Обобщение данного метода интерполяции на трехмерный вариант развит в нашей работе [62].

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

Шаг 1. Для измеренных значений синограммы f методом ART находится первое приближение функции fd.

Шаг 2. В полученную оценку синограммы вносится вся известная в пространстве Радона априорная информация.

Шаг 3. По значениям проинтерполированных проекций методом FBP находится первое приближение для томограммы.

Шаг 4. Теперь уже в пространстве томограммы также вводится вся доступная физическая априорная информация.

Шаг 5. Вычисляется преобразование Радона от найденной томограммы, и далее осуществляется переход на Шаг 2.

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

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

- гладкость функций, вводимая сглаживанием;

- нулевые граничные условия;

- свойства симметрии и периодичности синограммы для оптически прозрачной плазмы;

- закон сохранения энергии и центра масс для отдельных параллельных проекций.

Сглаживание проводится либо скользящим усреднением в окне 3X3, либо регуляризованными сглаживающими кубическими сплайнами, которые легко адаптируются к заданному уровню измерительных шумов [4].

В цикле проведенных расчетов [1] рассматривались две схемы измерения проекций - веерная ("идеальная") и система УИ/Г, используемая для эмиссионной томографии плазмы голландского токамака 11ТР[55]. Для веерной системы предполагается, что система состоит из 5 камер, с 15 детекторами каждая; при интерполяции в систему параллельных проекций используется 11 ракурсов и 27 детекторов на каждом (рис.4,а). Для системы V] 1Л' 5 камер с 16 детекторами дают расположение точек в пространстве синограммы, приведенное на рис.4,Ь.

__Г7 Ч ' ' \ ■

Рис. 5: а)7\ггних томограмма модели -V 1. Ь) Точная синограмма модели Л' 1.

В качестве одной из моделей точной томограммы для исследования свойств алгоритма восполнения синограмм выбрано асимметричное кольцо. На рис.5,а изображена точная томограмма модели N 1, а на рис.5,Ь -соответствующая ей синограмма.

Для иллюстрации характера зависимости погрешности реконструкции синограммы и томограммы от номера итерации при фиксированном параметре релаксации Л на рис.6 приведены нормы погрешностей такой реконструкции для описываемой модели Л^ 1, А = 2.0. Данные кривые показывают, что достижение минимума (или неубывание) нормы невязки может служить критерием для автоматического окончания ■ 1е •» >» " ™ " итерационного процесса. В алгоритме

МСП используется именно этот критерий останова итерационного процес-

Рис. С: _ са.

Зависимость норм погрешностей ре- .,

конструкции для Л = 2.0 синограмм Восстановление методом фильтрации и (2), томограмм (1) и норм нсвячки для обратного проецирования (ГВР, веер-

синограммы (3) от номера итерации Ь ный варнант) дает близкие К МСП ре-(в процентах). „ '

чультаты. Погрешность реконструкции

томограммы при этом равна Д'^"1 = 37.82%, что даже хуже на 2%, чем лучший результат алгоритма МСП. Погрешность отклонения точной синограммы,-и полученной по восстановленной томограмме - Д"" = 11.79%.

Главное же достоинство алгоритма МСП перед И$Р (и многими другими, созданными для веерной томографии) в том, что у него шире область применений - он работает при нерегулярно расположенных детекторах в пространстве синограммы, когда алгоритмы ГВР просто неприменимы.

2.3,6 Томография объектов с непрозрачными включениями

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

Алгоритмические разработки, относящиеся к двумерным задачам с учетом возможной экранировки части луч-сумм, осуществлялись многими авторами (Cho, Vest, Oppenheim, Lewitt, Bates). В этих алгоритмах используется ряд характерных допущений: границей тела, как правило, считается окружность, начало координат совмещается с ее центром и т.п. При такого рода подходах к задаче довольно естественно также применение метода послойного расщепления.

Рассмотренный далее алгоритм [1] отличается от упомянутых в нескольких отношениях. Расположение непрозрачного включения в пределах финитного носителя искомого 2D распределения g(r) может быть произвольным. Возможна любая, но заранее известная выпуклая граница включения независимо от условий гладкости контура (например, треугольник, прямоугольник и т.д.). Наконец, допустимо различное число включений внутри рассматриваемой области.

Схематично геометрия задачи представлена на рис.7. На этом рисунке показаны два включения: круглое (1) и квадратное (2), а также одна из 'f (^/проекций /¿(р) с двумя "провалами", соответствующим затененным участкам. Подобная ситуация реализуется, например, при интерферометри-ческих исследованиях обтекания непрозрачных моделей газовым потоком. Каждой точке вне включений (скажем,

__В) отвечают определенные обл;

Рис. 7: Схема томографии объекта с непрозрачны- .

ми включениями сти потери информации; на

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

возрастает до 7Г.

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

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

Шаг 2. По совокупности дополненных проекций восстанавливается 2£> объект г) в первом приближении.

Шаг 3. Корректируется решение д^(т) на основе имеющейся априорной информации, в том числе выполняется процедура искусственного "просветления" включений.

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

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

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

Некоторые результаты применения описанного алгоритма показаны на рис.8 [2]. Здесь изображены два центральных сечения д(х,0) и д(0,у) двумерного объекта д(х,у), представляющего собой смещенную относительно начала координат гауссиану, которая нормирована на единицу в максимуме. Непрозрачным включением был круг радиусом 0.15 с центром в начале координат. В процессе "просветления" включение заменялось константой, равной 0.25. На томограмме выбирался растр 41 X 41; уровень "шума" в проекционных данных считался равным 5%; число ракурсов К = 7. Шаг 2 в итерационном процессе реализовался на основе алгоритма реконструкции Ев 1, описанного в [19]. Погрешность восстановления оказалась порядка 10%. Оценки по нескольким модельным экспериментам показали, что, если сравнить погрешность восстановления д'1'(г) в первой итерации с погрешностью после 10 - 15 итераций, то последняя, как правило, ниже первой в 2 - 3 раза.

0Л5

ais

wo

eso

0.5 10 x -lo -0.5 o

b

(.0 - e.5 o

0.5 vo

F'nc. S: Реконструкция неоднородного 2D объекта g(x,y), содержащего в центральной области осесимметричное непрозрачное включение (положение включения показано штрих - пунктирными линиями). Сплошные линии - исходный объект; штриховые - результат восстановления после 5 итераций. Сечения: а - jf = О, b - х = О

2.3.7 Методы реконструкции томограмм с учетом рефракции зондирующего излучения

Интерферометрия газовых или плазменных сред является широко распространенным диагностическим средством в современных исследованиях (Островская, Островский, Vest). Такой эксперимент позволяет визуализировать интегральный набег фазы световой волны, прошедшей сквозь изучаемую оптическую неоднородность. Переход от интегральных характеристик объекта к его локальным свойствам требует применения соответствующих методов томографической интерпретации. Последние наиболее полно развиты в линейном приближении [1, 2], когда пренебрегают явлениями рефракции. Учет рефракции приводит к нелинейным задачам томографии, методы решения которых еще недостаточно развиты. Интерес же к подобным постановкам возрастает все больше и больше (Aben, Andrienko, Dolovich, Gladwell, Lira, Vest, Шарафутдинов). В данном разделе описан метод определения возмущения показателя преломления по измеренным набегам фаз зондирующей волны, развитый в работе [39], и учитывающий реальные траектории лучей в объекте.

Рассмотрим задачу восстановления изменения показателя преломления Дп(х,у) в некоторой выделенной плоскости по измеренным набегам фаз Д <р световой волны:

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

(38)

2G

\

\

\

\

\

\

\

Рис. 9: Схема зондирования плоской волной оптической неоднородности.

Если семейство лучей, вдоль которых выполняется интегрирование в (38), суть прямые линии, то в классической задаче томографии уравнение (38) представляет собой линейное преобразование Радона (ЛПР) от функции 27гДп/А:

/(£>?)=# о{д(х,у)} = / д(х,у)6(р + хвгп£-усо5$)(1х(1у. (39)

оо 1 — оо

Здесь Но — интегральный оператор Радона, д(х,у) = Лп(х,у) — искомая функция, f = АД уз/2зг - проекционные данные.

В общем случае криволинейных лучей имеем запись нелинейной задачи интегральной геометрии (Лаврентьев, Романов, Шишатский, 1980)

где Д — оператор криволинейного интегрирования, который логично назвать обобщенным преобразованием Радона (ОПР). Нелинейность задачи (40) заключается в том, что траектории интегрирования заранее неизвестны и зависят от искомого показателя преломления п(г).

Решение уравнения (39) для случая прямолинейных лучей дается обратным линейным преобразованием Радона (ОЛПР):

которое осуществляется в два этапа. Первый этап - это свертка проекции

Щд} = /,

(40)

9 = Щ1{П = в0{ПР)},

(41)

/К,Р):

/<(Ро) = (/ * Л) = №,р) Цр - р0)с1р,

(42)

где ядро Н(р) в общем случае есть обобщенная функция —1/[2тг2(р — Ро)2]> Ро ~ —ж ягп£ + усо«£. Второй этап ОЛПР - это реализация оператора обратного проецирования В0'.

Если эффекты рефракции существенны, то использование алгоритма ОЛПР (42)-(43) не дает требуемой точности. Действительно, в этом случае вклад в проекцию /(£, р) в точке Р дается не интегрированием функции д(х,у) вдоль точек прямой р = —xstn$ + ycos£, а интегрированием вдоль траектории ст(г), определяемой уравнением эйконала

В (Cha, Vest, 1981) предложен итерационный алгоритм определения поля показателя преломления га(г) из (40) (далее этот алгоритм обозначается сокращением CV). В основу этого алгоритма положено построение оценок таких оптических путей Fl(£,p), которые получились бы от данного объекта п(г) в случае прямолинейных траекторий лучей. Последовательное нахождение этих оценок для линейной части проекций F1 позволяет строить очередные приближения для искомой функции применением обычного ОЛПР: д' — Я"1^'}. Для определения "псевдопроекций" R{g'} в (Cha, Vest, 1981) применялся алгоритм лучевой трассировки.

Для улучшения свойств сходимости изложенного алгоритма и ускорения времени счета на ЭВМ была разработана его модификация CV1 [39]. Прежде всего в схему итерационного процесса были введены параметры релаксации 0 < а, < 1, которые при a¡ ^ 1 делают схему более устойчивой, и возможная ее расходимость (наблюдаемая иногда как в (Cha, Vest, 1981), так и в наших численных экспериментах) проявляется при большем номере итерации. Кроме того, оператор Rр1 реализован в варианте обратного проецирования с фильтрацией (ОПФ) вида (41)-(42), тогда как в (Cha, Vest, 1981) используются полиноминальные разложения и суммирование рядов, что ставит непростую проблему выбора конечного числа членов ряда.

Далее описаны свойства алгоритма CV 1 в сравнении с новым методом, названным методом обращенной волны (ARW).

Метод обращенной волны

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

В первую очередь при вычислении поправок g¡ к д0 в уравнениях СVI заменим постоянный оператор R'1 на набор операторов ft"1, более точно аппроксимирующих линейную часть нелинейного оператора R'1. Поскольку в заданную точку Р проекции f(p) информация из объекта приходит по криволинейной траектории Lp, а не по прямым Ln (рис.9), то и при вычислении вклада проекции в решение д(х,у) на этапе обратного проецирования логично использовать величину /(р), а не /(pi), как это делается в методе Ча - Веста и его модификации СVI.

(43)

(VS)2 =п2(х,у).

(44)

Как уже отмечалось, составная часть ОЛПР - этап обратного проецирования В о — трактуется как "растяжение" значений отфильтрованных проекций /{(р) обратно тому направлению, вдоль которого они регистрировались. Образно это также можно представить как зеркальное отражение /¡(р) от плоского зеркала Ох назад к входной плоскости И (рис.9). Обобщим эту операцию зеркального отражения, и вместо плоского зеркала подставим зеркало обращения водного фронта (Зельдович, Пилипецкий, Шкунов, 1985), которое "растянет" проекцию /((р) (или отфильтрованную проекцию /((р)) вдоль реальных траекторий £г.

Тогда оператор обратного распространения волны определится как [39]:

ВП{Д$,Р)}= Г</Ш<г(г„)]. (45)

Jo

В этой формуле для заданного направления $ и фиксированной точки г0, в которой вычисляется результат действия оператора Вц, аргумент и определяется как координата той точки Р на экране Ог, в которую пришел луч из точки Го-

Вторая часть процедуры инверсии обобщенного преобразования Радона Л'1 - это этап фильтрации. Для обобщения этого оператора также воспользуемся наглядными эвристическими представлениями [30]. На рис.10,а условно изображена функция д(г), представляющая собой константу внутри круга малого радиуса и ее отфильтрованную проекцию /((р). Поскольку выбранная £ ^функция д(г) есть приближение к ^ ¿-функции, результат фильтрации ее проекции дает приближенное представление ядра свертки /1(р). На рис.10,а изображен частный случай известного ядра Шеппа-Логана. Результат обратного проецирования В0{/} показан здесь же, при-полосы" функ-

Рис. 10: Характеристическая функция круга и ее фильтрованные проекции в алгоритме Шеппа-Логана (а); нелинейная филь- чем положительные

трация импульснои проекции в направлении •* т / \ { = 0 и ее "след" на томограмме, получен- «ии П(Р) на ПОЛе томограммы Дании оператором обратного распространения ны прямой штриховкой, а отрица-волны Вц (Ь).

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

необходимо стремиться к такой же компенсации и для "криволинейных" полос. Этого частично можно добиться, если внести оператор нелинейног о проецирования Т, КОЮрЫЙ переводит ПЛОСКОС ТЬ ( = 1.1 в t = f3 по криволинейным траектория.".!, и ею обращение, выполняющее обратный переход (см. рис.10,Ь). В случае бе;рефракционного приближения эти операторы тождественно единичные. Т 1 = Т ГГ Е. В об/нем же случае эти операторы позволяют перевести друг в друга начальную и конечную точки траектории <г(г). Обозначим через \\{р). Тогда ''новую" проекцию ij>i(p) уже можно фильтровать обычным образом

<?<№) = Q{Vf(p)> = -(2л-3)-1 I >.-.Jr)h(p-f,„)dp, (-10)

Для последующего обращения волнового фронта необходимо в полученных 1 i'f(p) выполнить обратный переход к коордпна i ам в плоскости < = i2, т.е. получить ip((Tp). Toi да краткая запись оператора есть

Л,(г) - ÄiTi-J*),

Vdp) = Q/t гг. ;'.(<}. (.IT)

В этих выражениях индекс (/ — I) означает зависимость оператора от поля показателя преломления п, t (rj с прелыдуmn и шага итерационного процесса

»..¡(г) — fl'-'(r) -t- II».

Итак, окончательно основное уравнение итерации и алгоритме обра-щенпок полны (AUW) вьплндцт г;;" j.v юш и -а ofipa'-coM:

= ,z' -' + a, н;1 [f- Щ;,'-1 4- п„} 'г «„{»о}}- (4Ь)

Регуляризация задачи oeyrm ст влиется сглаживанием проекций кубическими сплайнами и прекращением процесса (48), согласованным с уровнем

!!!} Л'.JH В I'ptJC*K'miCliiiii>IX Д.'Ш'.'ЬГХ.

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

2.4 Трехмерные чадачи томографии

Введение

Наиболее сложными (как в экспериментальном, так и в вычислительном плане) являются задачи трехмерной эмиссионной томографии, в которых собственное излучение всею объема плазмы регистрируется двумерными дегск! орали [31].

0("мчно число доступных ракурсов наблюдения в таких экспериментах весьма невелико (порядка 3-10), и реконструкция томограмм объекта

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

В работах автора [31, 37, 42, 48, 49, 50], [51, 53, 54, СО. 61, 62, 64] разработаны новые алгоритмы трехмерной томографии, в которых развиты и обобщены методы, созданные нами рацее для задач двумерной томографии. К таким алгоритмам, в частности,относятся регулярнзовашплй вариант ОПФ, итерационный модифицированный а.т-Риг. 11: Схема регистрации двумерных проекции п горитм Гсршберга - IJany-задаче трехмерной эмиссионном томографии лиса (3£)Д/GP), ART

M ART и MENT.

Далее кратко упомянут алгоритм ОПФ, испытанный нами в реальном эксперименте, а также дано описание итерационного алгоритма 3DMGI'. Для последнего даются определенные в численном моделировании точностные характеристики.

Рассмотрим постановку задачи трехмерной эмиссионной томографии плазмы (рис.II).

Здесь g (г) - искомое трехмерное распределение локальных коэффициентов эмиссии плазмы, /е1>5(и,г>) - дну мерная проекция, т.е. интегральная светимость плазмы, зарегистрированная каким-либо двумерным детектором излучения (например, фотопленкой). Пренебрегая эффектами рефракции, рассеяния и самопоглощения, для параллельных систем регистрации излучения получаем связь между интегральными и локальными характеристиками плазмы в виде интегрального уравнения ::учевого Р - преобразования (Natterer, 1986):

f.,.(u,v)- I g(r)Jl,

(49)

где М - дифференциал вдоль прямой, характеризуемой направлением . и координатами (м,и) в плоскости Ь'0\'.

Задачей трехмерной томографии является решение уравнения (49) относительно д(г) при шнестных наборах функции (и, у),г =

Методам нахождения подобных решений в настоящее время посвящено уже значительное число работ, особенно связанных с медицинской томографией и промышленной дефектоскопией (см. обзоры и основополагающие работы в Ре1сПсатр, Спи^еа!, ОиИЬе^, Кп^ббоп, Е<Пю]т, Сгап1ши),

Petersson, Tuy, Minerbo, Herman, Орлов).

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

Здесь через q обозначен вектор точки в плоскости детектирования двумерной проекции /(ч).

В случаях со сложным расположением приемников излучения в пространстве, алгоритм обратного проецирования с фильтрацией, опирающийся на формулу (50), может дать приемлемую оценку плазменного распределения. Подобный алгоритм применялся нами в работе [31] для восстановления распределения рентгеновской эмиссии плазмы в исследованиях по лазерному термоядерному синтезу на установке "Кальмар".

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

2.4.1 Модифицированный алгоритм Гершберга-Папулиса

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

Пусть д(г) = у(х,у,г) - функция трех переменных, fs¡v(u,v) - ее параллельная проекция в направлении Обозначим также их Фурье -образы, соответственно, и /(</„, Тогда значение функции <7(1^, 1>у, Уг) на плоскости,проходящей через начало координат и имеющей вектор нормали совпадает с функцией f(^>u,^'„). Допустим, имеется N проекций f|l¡.í (u,v). Согласно проекционной теореме, их Фурье - образы совпадают с сечениями трехмерного Фурье - образа решения д(и) плоскостями с теми же направлениями векторов нормалей и проходящими через начало координат.

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

(50)

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

д, = Ф.Г-'д,.

&

Фl,Fg¡-l, г >

= {Л • ^ (эй

где Р и Р1'1, соответственно, операторы прямого и обратного трехмерного преобразования Фурье, фх и - операторы внесения априорной информации в пространственной и частотной областях, соответственно, д0 - Фурье-образ, полученный в результате начальной экстраполяции.

Данная трехмерная модификация алгоритма Гершберга-Папулиса является обобщением соответствующей двумерной модификации, развитой автором в работе [45].

Регуляризация итерационного процесса

В предыдущем разделе предполагалось, что восстановление происходит по идеальным данным, в том смысле, что Фурье - образы проекций в плоскости дают сечение Фурье - образа точного решения. Однако, проекционные данные, полученные в ходе какого-либо эксперимента, содержат некоторые погрешности - шумы измерений. В настоящее время разработано множество методов для подавления шумов как в пространственной, так и в частотной областях. Шумоподавление в пространственной области часто основано на том, что предполагается известным распределение шумов. Обычно его предполагают гауссовским с нулевым средним. В алгоритм ЗБМСР входит преобразование Фурье, которое при численной реализации заменяется дискретным преобразованием Фурье и, в конечном итоге, приводит к суммированию ряда. В работе (Тихонов, Арсенин, 1986) указан метод устойчивого суммирования рядов Фурье с неточно заданными коэффициентами ряда, который заключается в следующем. Пусть дан ряд Фурье некоторой функции

£(*)=£<:-.¥>«(*), (52)

в котором коэффициенты с„ известны приближенно. Регуляризованный метод суммирования заключается в том, что вместо ряда (52) суммируется ряд с другими коэффициентами с„ = с„/(1+а£„), где £„ - числа, порядок роста которых на бесконечности не ниже, чем 2"+л, Л > О, а - параметр регуляризации. При этом получается функция да, являющаяся приближением д. Величина а определяется, исходя из погрешности исходных данных, например, по критерию невязки:

Ус2 -= <т2, (53)

^ "(1+ «*„)'

где <т2 - суммарная дисперсия сигнала.

Вычислительный эксперимент

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

а) какова точность реконструкции изображения при малом числе ракурсов наблюдения;

б) какова чувствительность такого решения к шумам измерений;

в) как влияют на итоговое изображение отличия в информативности проекций, полученных в разных направлениях;

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

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

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

Для контроля сходимости итерационного процесса рассмотрим следующую норму ошибки решения:

Я 91

Ю02,%. (54)

Здесь и далее дЕ - точное решение, д''' - оценка решения, полученная на 1-оп итерации.

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

( 4 1п 2 (к — х0)2 41п 2 (у — уо)2 41п 2 (г - 20)2 ] д(х, у,г) = а ехр |------—-----1 ,

(55)

где а,/3,7 - полные полуширины гауссиан по координатным осям.

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

При этом величина Д1 стремится к постоянному пределу, который в значительной степени зависит от полярного угла под которым располагаются проекции (в численных экспериментах все проекции имели одинаковый полярный угол, а по азимутальному углу они были распределены равномерно) и, в меньшей степени, - от числа проекций. Сходимость алгоритма можно про-Рис. 12: График зависимости нормы отклоне- следить, например, на рис.12.

ния Дт от номера итерации X для модели N 1. « _ „ _________

* 11 Здесь показан результат восста-

новления гауссианы с параметрами а — /3 = 7 = 0.4, х0 = у0 = ^а = 0.0,

6050 40 20 го 10 о

л.

ю ¿

где ®0,Уо,20 20

Ш

I ' I

а

-X.

10

15

20 I

Рис.

СЛОЯ

координаты центра гауссианы.

Этот фантом в дальнейшем будем называть моделью N 1. На рис.12 приведено отклонение Дд, причем верхняя кривая соответствует углу в = 30°; нижняя -__в = 60°, К = 20 проекциям. Видно, что при уменьшении полярного угла отклонение полученного решения от точного заметно возрастает.

В численном моделировании изучено влияние толщины слоя "влияния" 6 на точность реконструкции. В местах изменения толщины слоя "влияния" на первых итерациях четко видны изломы зависимости величины погрешности решения от номера итерации I. Преимущество, получаемое при постоянно уменьшающейся толщине слоя "влияния", подтверждается рис.13(а,Ь). На рис.13(а) дана ошибка восстановления модели N 1, а на рис.13(Ь) - модели ./V 2: две гауссианы, разведенные по оси х: ж01 = 0.5, х02 = —0.5 (в обоих случаях в = 60°, 20 проекций). Здесь и далее

изменяющейся толщины на норму отклонения :

13: Влияние "влияния" 8 а)модель N 1; Ь)модель N 2.

стрелка отмечает номер итерации, после которой толщина слоя уменьшается, £1 - значение этой величины в шагах сетки (далее все величины толщины слоя "влияния" будут измеряться в этих единицах), 60 - вели-

чина слоя "влияния" при первичной экстраполяции.

На рис.13,(а) непрерывной линией нанесена кривая для итерационного процесса, при котором использовалась постоянная толщина слоя "влияния", равная 1, а штриховой - кривая для переменной величины слоя "влияния" ¿о = 1,¿i = 0.5,62 = 0.25. На рис.13,(Ь) кривая 1 - постоянный слой, 2 - переменный: = 0.125,= 0.0625.

2.4.2 Проблемы программной реализации алгоритмов томографии

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

Далее очень сжато описан пакет TOPAS — MICRO (см. [1, 2]), ориентированный на получение томографических изображений широкого круга объектов как в численном, так и в реальном физическом эксперименте. Разработка этого пакета начата в 1977 году, и на первых порах он включал в себя алгоритмы одномерной и двумерной томографии. Позже были разработаны блоки трехмерной томографии.

Пакет написан на языке Фортран-77, использовался ранее на ЭВМ БЭСМ-6, позже был адаптирован на персональном компьтере РС386/486. Комплекс построен на принципах модульного структурного программирования, с момента создания постоянно совершенствуется и пополняется новыми функциями и возможностями [14, 30].

Модульный анализ задач томографии.

При создании проблемно - ориентированных пакетов важную роль играет модульный анализ (Яненко, Карначук, Коновалов, 1977; Ковеня, Яненко, 1981).

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

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

Программы, ориентированные в значительной степени на обработку эксперимента, обязательно должны содержать блоки, по шоляюнше проводить вычислительный эксперимент (ВЭ). Это необходимо _члч н-еет.>-ринкего анализа характеристик новых алгоритмов, включаемых ц паком, а также для выполнения немаловажного этапа оптимизации томографических установок или схем проведения реального эксперимента {l'O).

Конкретная модульная струк rvpa ТОР AS — MICRO складывалась пит влиянием специфики математического аппарата КТ. а также прикладных интересов ангора в области диагностики газа и плазмы.

Основные стадии решения задачи КТ, очевидно, следующие: построение сеток по угловой и пространствен ной переменным для проекций; ввод (режим РЗ), ;шбо вычисление (режим ВЭ) соответствующих \»тим сеткам проекционных данных; расчет двумерной или трехмерной сеток, нл которых необходимо реконструировать томограмму ; наконец, непосредственно стадия самой реконструкции; далее следуют постобработка и интерпретация т омограмм, а также их визуализация.

Архитектура пакета.

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

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

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

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

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

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

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

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

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

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

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

Результаты использования пакета TOPAS — MICRO в реальных физических исследованиях, главным образом плазменных и газодинамических, описаны в наших работах [1, 4, 10, 13, 17, 18, 21, 23, 29, 50, 51], [54, 55, 58, 63].

3 Основные результаты и выводы

1. Разработан метод статистической регуляризации для решения зада-

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

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

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

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

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

6. Разработаны двумерные и трехмерные нелинейные модификации алгоритма Гершберга-Папулиса с максимальным учетом априорной информации об искомой томограмме.

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

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

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

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

Основные публикации по теме научного доклада

[1] Пикало« В.В., Мельникова Т.С. (1!И)5) Томография плазмы. Новосибирск: Наука, 229 с.

[2] Пикало« В.В., Преображенский Н.Г. (1987) Реконструктивная томография н газодинамике и физике плазмы. Новосибирск: Наука. -230

г.

[3] Жуков М.Ф. и др.(среди них - Никалов В.В.) (1987) Теория термической э.чсктрадуговой плазмы. Ч.И. Нестационарные процессы в термической плазме. Новосибирск: Наука. -202 с.

[1] Прсображснский II.Г., Пикало» В.В. (1982) Неустойчивые задачи диа-гнем гики плазмы. Новосибирск: Наука. -238 с.

[5] Пикалон В.В., Преображенский Н.Г. О преобразовании Абеля при топографической интерферометрии точечного взрыва /,' ФГВ. -1974. -Т.10, N 6. -С.923-930.

[в] Пикалон В В., Преображенский Н.Г. О восстановлении локальных характеристик плач мы а условиях ограниченной экспериментальной информации // Опт. спектр. -1070. -Т.40. N 6. -С.1094-1096,

[7j Пикалон В.В., Преображенскнй Н.Г., Гиппиус К.Ф., Колесников В.Н. К вопросу о надежности определения локальных значении интенсивности излучения плазмы // Физика плазмы. -1978. -Т.4, N 4. -С.926-930.

¡8; Пнкалоа В.В. Некорректные задачи локальной оптической диагностики газовых и плазменных объектов произвольной конфигурации // Инверсия Абеля и ее обобщения. -1978. Новосибирск: ИТПМ СО АН СССР, -С.25-67.

[!)] Пикллов В.В., Преображенский Н.Г. Локальная диагностика газовых и плазменных объектов с заданными изолиниями // Опт. спектр. -1979. -Т.46, N 1. -С.209-211.

[10] Климкин В.Ф., Пикалов В.В. О возможностях микроинтерферометрии при исследовании нестационарных процессов // 11МТФ. -1979. -N 3. -С.14-26.

[11] Пикалов B.B. Программа восстановления локальных коэффициентов эмиссии эллиптического плазменного объекта (ELLIPS) (ИТПМ, 1979, 28 е.), П004148 // Алгоритмы и программы. Информационный бюллетень ГФАП. -1980. -N 2(34). -С.48.

[12] Мельникова Т.С., Пикалов В.В., Преображенский Н.Г. О локальной диагностике оптически плотной асимметричной плазмы // Опт. спектр. -1980. -Т.48, N 3. -С.474-479.

[13] Melnikova T.S., Pickalov V.V. Temperature field measurements of electric arc plasma in longitudinal magnetic field // Beitr. Plasmaphysik. -1982. -Vol.22, N 2. -P.171-180.

[14] Пикалов B.B. Комплекс программ томографической реконструкции локальных характеристик плазмы // Комплексы программ математической физики (VII СКПМФ). Новосибирск: ИТПМ СО АН СССР, 1982. -С.317 - 319.

[15] Пикалов В.В. Модификации сверточного алгоритма оптической томографии газового потока // ШВсес. школа по методам аэрофизических исследований (сб. докладов). 4.2. -1982. Новосибирск: ИТПМ СО АН СССР, -С.232 - 235.

[16] Пикалов В.В. Инверсия Радона для треугольных плазменных конфигураций // Материалы III Всес. школы - семинара: Современные методы магнитного удержания, нагрева и диагностики плазмы. Харьков: ХФТИ, 1982. -С.160 - 162.

[17] Черевко A.C., Пикалов В.В., Тагильцев А.П., Юделевич И.Г., Энгель-шт B.C., Жеенбаев Ж.Ж. Изучение температурного поля плазменной струи двухструйного плазмотрона // ЖПС. -1983. -Т.38, N 3. -С.497-499.

[18] Гиппиус Е.Ф., Илюхин Б.И., Лунин H.H., Пикалов В.В. Особенности восстановления радиальных распределений интенсивностей излучения плазмы стеллараторов // Физика плазмы. -1983. -Т.9, N 5. -С.1111-1115.

[19] Мельникова Т.С., Пикалов В.В. Исследование плазмы с помощью плазменного томографа. -Новосибирск, 1983. -47 с. -(Препринт/ ИТ СО АН СССР; N 99-83).

[20] Пикалов В.В. Проблемы разработки и сопровождения программного обеспечения задач плазменной томографии // Всес. симпоз. вычисл. томографии. Тезисы докл. Новосибирск: ВЦ СО АН СССР, 1983. -С.150 - 151.

[21] Klimkin V.F., Pickalov V.V., Soloukhin R.I. Optical interferometry of spherical shock waves // Shock waves, explosions, and detonations. (Progr. in astronautics and aeronautics). -1983. -Vol.87. N.Y: AIAA, -P.196-204.

[22] Pickalov V.V., Melnikova T.S. Tomographic measurements of plasma temperature fields // Beitr. Plasmaphysik. -1984. -Vol.24, N 4. -P.417-430.

[23] Melnikova T.S., Pickalov V.V. Tomographic measurements of temperature fields in non-stationary arc plasma // Beitr. Plasmaphysik. -1984. -Vol.24, N 5. -P.431-445.

[24] Пикалов B.B., Преображенский H.Г. Атомно - спектральная томография // Изв. АН СССР, сер. физич. -1984. -Т.48, N 7. -С.1289-1296.

[25] Yudelevich I.G., Cherevko A.S., Engelsht V.S., Pikalov V.V., Tagiltsev A.P., Zheenbajev Zh.Zh. A two-jet plasmatron for the spectrochemical analysis of geological samples // Spectrochimica Acta (B). -1984. -Vol.39, N 6. -P.777-785.

[26] Жуков М.Ф., Мельникова T.C., Пикалов B.B. Томография низкотемпературной плазмы // Изв. СО АН СССР, сер. техн. наук. -1984. -N 10, вып.2. -С.39-46.

[27] Пикалов В.В., Преображенский Н.Г. Вычислительная томография в газовой динамике и физике плазмы // Некорректные задачи математической физики и анализа. -1984. Новосибирск: Наука, -С.106-111.

[28] Мельникова Т.С., Пикалов В.В. Эмиссионная томография нестационарной плазмы // ТВТ. -1984. -Т.22, N 4. -С.625-633.

[29] Мельникова Т.С., Пикалов В.В. Влияние расхода газа на температурные распределения дуговой плазмы // ТВТ. -1984. -Т.22, N 6. -С. 10671075.

[30] Пикалов В.В. Пакет прикладных программ, ориентированный на задачи вычислительной томографии // Вопросы реконструктивной томографии. -1985. Новосибирск: ВЦ СО АН СССР, -С.132-135.

[31] Денисов В.И., Захаренков Ю.А., Кологривов A.A., Пикалов В.В., Преображенский Н.Г., Рупасов A.A., Склизков Г.В., Шарапова Н.В., Ши-канов A.C. Методы эмиссионной и интерферометрической томографии лазерной плазмы // Краткие сообщения по физике. -1985. -N 12. -С.17-21.

[32] Baev V.K., Bazhaikin A.N., Buzukov A.A., Pickalov V.V. The spatial temperature distribution in nonstationary burning spray of pulverized liquid fuel // Archivum Combustion. -1986. -Vol.6, N 2. -P.115-124.

[33] Мельникова Т.С., Пикалов В.В. Томографическая диагностика низкотемпературной плазмы // Изв. СО АН СССР, сер. техн. наук. -1987. -N 11, вып 3. -С.60-68.

[34] Пикалов В.В. Математическое моделирование в задачах физической томографии // III Всес. симпоз. вычисл. томографии. Киев: Наукова Думка, 1987. -С.153-154.

[35] Вишняков Г.Н., Левин Г.Г., Пикалов В.В. Методы нормировки данных при количественной интерпретации томографических интерфе-рограмм // Опт. спектр. -1987. -Т.62, N 6. -C.1361-13G6.

[36] Климкин В.Ф., Пикалов В.В. О методе расчета распределения показателя преломления при интерферометрии сферических объектов // Опт. спектр. -1988. -Т.64, N 1. -С.159-164.

[37] Vishnyakov G.N., Drozhbin Yu.A., Levin G.G., Pickalov V.V., Ponomaryev A.M., Trofimenko V.V., Ushakov L.S. Chronotomography: a new approach to high speed frame recording // High speed Photography, Videography, and Photonics IV (Proceed. SPIE, vol.693). 1986. -P.170-173.

[38] Пикалов В.В. Методы оптической томографии плазмы // Оптическая томография (Тез. докл. Всес. семинара). Таллин: ИК АН ЭССР, 1988. -С.136-140.

[39] Пикалов В.В. Восстановление томограммы прозрачной неоднородности методом обращенной волны // Опт. спектр. -1988. -Т.65, N 4. -С.956-962.

[40] Аксенов В.П., Пикалов В.В. Томографическая реконструкция критериев качества лазерных пучков в атмосфере // Оптика атмосферы. -1988. -Т.1, N 10. -С.30-35.

[41] Бронников А.В., Воскобойников Ю.Е., Пикалов В.В., Шарапова Н.В. Сравнение алгоритмов восстановления по веерным проекциям // Электронное моделирование. -1988. -Т.10, N 6. -С.81-86.

[42] Пикалов В.В. Многомерная томография: алгоритмы и приложения // IV Всес. симпоз. вычисл. томографии. Тез. докл. 4.1. Новосибирск: ИМ СО АН СССР, 1989. -С.127-128.

[43] Пикалов В.В. Возможности пакета программ TOPAS для задач компьютерной томографии // Математическая обработка и интерпретация результатов физических экспериментов. -1989. М.: Энергоатом-издат, -С.57-63.

[44] Аксенов В.П., Пикалов В.В. Томографическое восстановление пространственно - энергетических параметров лазерных пучков // Квантовая электроника. -1990. -Т.17, N 2. -С.167-172.

[45] Melnikova T.S., Pickalov V. Computer-aided plasma tomography // High temperature dust-laden jets in plasma technology: Proc. Intern. Workshop, 6-8 Sept. 1988, Novos.,USSR. Eds: O.P.Solonenko, A.I.Fedorchenko. -1990. Utrecht: The Netherlands, VSP, -P.257-282.

[46] Zakharenkov Yu.A., Kosterin A.V., Shikanov A.S., Pikalov V.V., Preobrazhensky N.G. Laser produced plasma refractometry // Laser and Particle Beams. -1990. -Vol.8, N 1-2. -P.339 - 342.

[47] Пикалов В.В. Физические принципы и вычислительные методы компьютерной томографии // Математическое моделирование в машиностроении. Секция I. Общие вопросы. Тезисы докл. Первой Всес. школы - конф. Куйбышев: КуАИ, 1990. -С.35-36.

[48] Баландин A.JI., Пикалов В.В. Трехмерный алгоритм максимума энтропии в малоракурсной оптической томографии потоков. // Оптические методы исследования потоков. Тезисы докл. I Всес. конф. Новосибирск: ИТ СО АН СССР, 1991. -С.83-84.

[49] Пикалов В.В. Трехмерная томография: предобработка, реконструкция и визуализация на персональной ЭВМ // V Всес. симпоз. по вычислит, томографии. Тезисы докладов. М.: ВНИИФТРИ, 1991. -С.222-223.

[50] Веретенников А.В., Кошевой М.О., Панферов Н.В., Пикалов В.В., Ру-пасов А.А., Семенов О.Г., Шиканов А.С. Эмиссионная томография плазмы микропинчевого разряда // Физика плазмы. -1992. -Т.18, N 2. -С.249 - 252.

[51] Алпатов В.В., Левин Г.Г., Пикалов В.В., Романовский Ю.А. Оптическая томография искусственных образований в околоземной среде // Космические исследования. -1993. -Т.31, вып.1. -С.121 - 134.

[52] Pickalov V.V., Ingesson L.C., Donné A.J.H., Schram D.C. Tomography algorithms for visible light tomography on RTP // 1992 Intern. Confer, on Plasma Phys., Innsbruck, June 29 - July 3. Abstracts, vol.l6C, Pt.II, Innsbruck: EPS, 1992. -P.1143-1146.

[53] Glazkov V.N., Kazantsev I.G., Pikalov V.V. Fast algebraic reconstruction algorithm in three-dimensional tomography // Analytical methods for optical tomography, Proc. SPIE. -1992. -Vol.1843. -P.58-61.

[54] Alpatov V.V., Bulygin F.V., Levin G.G., Pikalov V.V., Romanovsky Yu. A. Optical tomography techniques applied to study geophysical artificial structures // Analytical methods for optical tomography, Proc. SPIE. -1992. -Vol.1843. -P.302-314.

[55] Ingesson L.C., Pickalov V.V., Donné A.J.H., Schram D.C. and RTP-team. First results with the visible light tomography system on RTP. // Proc. 20th EPS Conf. Controlled Fusion and Plasma Physics, Lisboa, 26-30 July 1993. Part III. 1993. -P.1147-1150.

[56] Balandin A.L., Likhachov A.V., Panferov N.V., Pickalov V.V., Rupasov A.V., Shikanov A.S. Emission microtomography of plasma // Analytical methods for optical tomography, Proc. SPIE. -1992. -Vol.1843. -P.68-82.

[57] Ingesson L.C., Pickalov V.V., Donné A.J.H., Schram D.C. and RTP-team. Interpretation of measurements from visible light tomography system on RTP // 21st EPS Conf. Contr. Fusion and Plasma Phys. Abstracts. Montpellier. EPS, 1994. -P.409.

[58] Ingesson L.C., Pickalov V.V., Donne A.J.H. and RTP - team. First tomographic reconstructions and a study of interference filters for visible light tomography on ЯТР. // Rev. Sci. lustrum. -1995. -V0I.6G, N 1. -P.622-624.

[59] Alpatov V.V., Pickalov V.V., Likhachov A.V. Informational analysis of auroral tomograph // Auroral Tomography Workshop. March 9-11, 1993. Kiruna. Sweden. IRF Technical Report 213. 1993. -P.35 - 46.

[60] Likhachov A.V., Pickalov V.V. Subpixel resolution in 3D cone - beam microtomography // Nucl. Instrum. Meth. Phys. Res. (A). -1995. -Vol.359, N 1-2. -P.370-375.

[61] Likhachov A.V., Pickalov V.V. A modification of ART method for cone-beam tomography of high space resolution // Computerized tomography. Proc. Fourth Intern. Sympos., Novosibirsk, Russia. Utrecht: VSP, 1995. -P.309-317.

[62] Likhachov A.V., Pickalov V.V. A new approach to the problem of 3D interpolation from arbitrary set of points // Computerized tomography. Proc. Fourth Intern. Sympos., Novosibirsk, Russia. Utrecht: VSP, 1995. -P.318-323.

[63] Alpatov V.V., Likhachov A.V., Pickalov V.V., Romanovsky Yu.A. Optical tomography of natural and artificial optical disturbances in the near terrestrial space // Computerized tomography. Proc. Fourth Intern. Sympos., Novosibirsk, Russia. Utrecht: VSP, 1995. -P.12-24.

[64] Лихачев А.В., Пикалов В.В. Частотная фильтрация в алгебраических алгоритмах трехмерной томографии // Автометрия. -1995. -N 4. -С.83-89.

[65] Likhachov A.V., Pickalov V.V. Frequency filtration in the algebraic algorithms of 3D-tomography // Proceed. 1995 Intern. Meet. Fully Three-Dimensional Image Reconstruction in Radiol. Nucl. Medicine, 4-6 July 1995, Aix-les-Bains. Grenoble: LETI, 1995. -P.103-107.

[66] Likhachov A.V., Maslov A.A., Mironov S.G., Pickalov V.V. A tomographic investigation of the density field in hypersonic flows // ICMAR'96, Proceed., Pt.II. Novosibirsk: ITAM, 1996. -P.169-174.

[67] Лихачев A.B., Пикалов B.B. Об одной постановке задачи хронотомо-графии // Опт. спектр. -1996. -N 4. -С.581-589.