CC BY
8
Вычислительные технологии, 2020, том 25, № 1, с. 5-38. © ИВТ СО РАН, 2020 Computational Technologies, 2020, vol. 25, no. 1, pp. 5-38. © ICT SB RAS, 2020
ISSN 1560-7534 elSSN 2313-691X
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ
D01:10.25743/ICT.2020.25.1.002
Малоракурсная реконструкция светимости плазмы для сферического токамака
А. Н. Баженов1'2, П. А. Злтылкин1,2
1 Физико-технический институт им. А.Ф. Иоффе РАН, Санкт-Петербург, Россия
2СПбПУ им. Петра Великого, Санкт-Петербург, Россия
Контактный автор: Баженов Александр Н., e-mail: a_bazhenov@inbox.ru
Поступила 28 мая 2019 г., доработана 11 ноября 2019 г., принята в печать 22 ноября 2019 г.
Публикация посвящена применению методов вычислительной геометрии, интервального анализа и линейного программирования к задачам физики управляемого термоядерного синтеза. Рассмотрены геометрические аспекты проблемы, получены проекции светимостей различных объемов сферического токамака на плоскость матричного детектора, изучены изображения предполагаемых макроскопических структур и микроскопических включений. Для набора модельных распределений светимости объема токамака поставлена задача восстановления сигнала. Решение получено с использованием задач линейного программирования.
Ключевые слова: интервальная система линейных алгебраических уравнений, допусковое множество, мера вариабельности, распознающий функционал, линейное программирование.
Цитирование: Баженов А.Н., Затылкин П.А. Малоракурсная реконструкция светимости плазмы для сферического токамака. Вычислительные технологии. 2020; 25(1):5-38.
Введение
Работа посвящена применению методов вычислительной геометрии и интервального анализа к задачам физики плазмы. Производится реконструкция светимости объема плазмы сферического токамака "Глобус-М" ФТИ им. А.Ф. Иоффе РАН [1] в области мягкого рентгеновского излучения.
Для физики управляемого термоядерного синтеза реконструкция светимости объема плазмы токамака включает как восстановление пространственной светимости плазмы в отдельные моменты времени разряда, так и ее эволюцию во времени. Такая реконструкция дает информацию о транспорте примесей в плазме и может использоваться совместно с другими диагностиками плазмы для ее всестороннего исследования. Обычно такие задачи изучаются с использованием нескольких детекторов в виде камер-обскур, расположенных в различных ракурсах по отношению к экваториальной плоскости камеры токамак. В основу алгоритма обработки пространственно-временных зависимостей в работе [2] положен метод Кормака (преобразование Радона), использующий разложение в ряд Фурье по угловой переменной и ортогональным полиномам Цернике — по радиальной. Получены реконструкции, воспроизводящие физически значимые топологические структуры плазмы.
Диагностика мягкого рентгеновского излучения сферического токамака "Глобус-М" ФТИ им. А.Ф. Иоффе РАН имеет ряд особенностей. Для регистрации излучения в течение ряда лет использовалась одна камера-обскура с ориентацией нормали в экваториальной плоскости токамака так, чтобы спроецировать на детектор как можно больший объем токамака. При таком расположении детектора задача реконструкции плохо обусловлена из-за значимой мультиколлинеарности хорд наблюдения. Геометрия сферического токамака существенно отличается от классической тороидальной геометрии. Форма сепаратрисы (последняя замкнутая поверхность магнитного потока, условная граница плазмы) напоминает букву "Б". Структура поверхностей равного магнитного потока очень далека от набора концентрических окружностей, и аналитические решения на основе элементарных (гармонических, эллиптических, полиномов) функций в духе [2] вряд ли реализуемы. Размеры вакуумной камеры токамака "Глобус-М" имеют масштаб 1 м, размер плазмы в меридиональном сечении составляет около 50 см. Размер одного пространственного элемента детектора — немного менее 1 мм, а самого детектора 3 см, поэтому различные неопределенности в геометрических данных и вычислениях имеют существенно большее значение, чем на крупных установках.
В диссертации [3] впервые на сферических токамаках (Глобус-М) получено восстановление двумерного распределения радиационных потерь плазмы в меридиональном сечении. В этой работе создан и применен алгоритм для реконструкции двумерного профиля радиационных потерь на основе хордовых измерений. Разбиение области выполнено с помощью прямоугольной сетки, что делает затруднительным разрешение областей светимости малого объема.
В связи с тем, что на токамаке "Глобус-М" накоплено большое количество данных, несколько десятков тысяч разрядов, есть мотивация к обработке этого материала как для задач осмысления уже полученной информации, так и для дальнейших исследований на модернизированном токамаке "Глобус-М2" [4].
0.1. Постановка задачи определения светимости элемента плазмы
в конкретный момент времени и временной эволюции светимости по хордовым наблюдениям
С математической точки зрения нахождение светимости различных областей объема плазмы по данным матричного детектора представляет собой обратную задачу. Точнее, такие задачи относятся к интегральной геометрии. В книге [5] дается определение таких задач как восстановление функции на некотором линейном пространстве по множеству значений этой функции на заданном семействе вложенных в это пространство многообразий.
За время разряда детектор производит последовательно большое количество измерений, при этом светимость плазмы непрерывно меняется. Фиксируем момент времени р = 1, 2,... М, N ~ 105. Пусть П — объем плазмы внутри токамака. Значения светимости I(х, у, г) плазмы положительны в объеме П и равны нулю вне этого объема:
Детектор представляет собой набор пикселей, которые расположены в виде двумерной структуры, так что они образуют "столбцы" и ортогональные им "строки". В геометрии камеры-обскуры излучение плазмы попадает в пиксели через диафрагму С пре-
> 0, если (х,у,г) € П, 0, если (х, у, г) </ П.
(1)
небрежимо малого диаметра. Через "столбцы" детектора и диафрагму можно провести плоскости Рг, г = 1, 2,... ,п. В плоскости Рг светимость I(х,у,г) принимает значения I(хг,уг), где (хг,уг) — координаты точек в рг. В каждой плоскости pi через конкретные пиксели в столбце детектора и диафрагму проходят лучи Ь^, ] = 1, 2,... ,т. Плазма токамака очень разрежена и не поглощает собственного излучения. В связи с этим светимость, регистрируемая детектором, от малого участка й1 луча вносит вклад й1 ■ 1г(хг, Уг). Уравнение прямой на плоскости можно записать в следующем виде:
х соъ(ф) + у 8т(<^) = р. (2)
Здесь |р| — длина перпендикуляра, опущенного из начала координат на прямую; — угол между осью х и перпендикуляром. Семейство прямых на плоскости представляет собой двухпараметрическое семейство Ь(р,р), р Е К, р Е (0, 2ж). Светимость, регистрируемая детектором, — интеграл по прямой Ь:
У I(х,у)М = /(р,р). (3)
Нр,1Р)
Формулу (3) называют преобразованием Радона функции I(х,у) (формула (5.1.2) в [5]).
Перейдем к дискретному варианту расположения пикселей детектора. Светимость, регистрируемая пикселем детектора, — интеграл по прямой Ь^:
(хг,уг)М = ¡гз. (4)
Здесь /у — данные пикселей детектора; г = 1, 2,... ,п — номер плоскости; ] = 1, 2,... , т — номер луча.
Перейдем к дискретизации объема П. Рассмотрим разбиение объема П на совокупность непересекающихся объемов П:
я
П= У Пк, П П П = 0, г = ]. (5)
к=1
Припишем каждому объему П значение светимости . В этом случае система (4) примет вид
я
1к = . (6)
к=1
Здесь — длина отрезка луча ] в плоскости г, пересекающего объем П. Общее
число уравнений в системе (6) равно пт, число неизвестных — д. В случае д < пт система матричных уравнений переопределена.
Чтобы привести матричное уравнение к виду, который удобно решать методами линейной алгебры, произведем перенумерование пикселей детектора в одномерном виде. Пусть в — номер уравнения в системе (6) и в = 1, 2,... пт. Тогда система уравнений (6) примет вид системы линейных алгебраических уравнений (СЛАУ)
я
А1к31к = Л
к=1
(7)
или в матричной форме
LI = I (8)
с матрицей L = Alsk. Уравнение (8) описывает распределение светимости в I-объеме плазмы в момент времени tp. Совокупность величин I(t) = {Ikописывает эволюцию светимости плазмы в течение конкретного разряда.
0.2. Интервальный анализ и его применение к решению СЛАУ
Для реализации этой задачи разработан математический аппарат для двумерной реконструкции светимости с учетом неопределенностей входных данных методами интервального анализа, система обозначений приводится в публикации [6].
Согласно [7], классическая интервальная арифметика IR — алгебраическая система, образованная интервалами х = [х, х] С R так, что для любой арифметической операции " к" из множества { + , — , • , / } результат операции между интервалами определяется как
X к у = { Х*у | X Е Х,у Е у }. (9)
Развернутые формулы для арифметических операций выглядят следующим образом:
х + у = [х + у, х + у], (10a)
х — у = [х — у, х — у~\, (10b)
х • у = [min{xy, ху, ху, ху}, max{xy, ху, ху, ху}], (10c)
х/у = х • [ 1/у, 1/у ] для у^ 0. (10d)
Действующий IEEE-стандарт для интервальных вычислений принят в 2015 г. [8].
Интервальные системы линейных алгебраических уравнений (ИСЛАУ) записываются как
ап^1 + ai2^2 + ... аыХ1 = 6Ь
О21Ж1 + а22Х2 + ... а2пХ2 = Ь2, (11)
+ Q"m2%2 + . . .
Кратко систему (11) можно записать следующим образом:
Ах = Ь, (12)
где А = (aij) — m х n-матрица, b — m-вектор правой части. Арифметические операции с интервалами — это операции классической интервальной арифметики [7].
Далее будут использоваться обозначения: mid а = —(а + а) — середина интервала,
rad а = —(а — а) — радиус интервала.
Решение системы уравнений (12) образуют множества решений, наиболее популярными из которых являются объединенное множество решений, образованное решениями всех точечных систем Ах = Ь:
Euni(A, b) = {х Е Rra | Ах = b для некоторых А Е А и b Е Ь} , (13)
и допусковое множество решении
ЕЫ(А, Ь) = [х е Кга | Ах е Ь для любых А е А}. (14)
Принадлежность к допусковому множеству можно тестировать путем вычисления знака распознающего функционала. Его выражение имеет вид
Tol (х) = Tol (х, A, b) = min < rad bi —
1<i<m '
Xj
.
:15)
mid bi — ^^ a,i
3 = 1
Принадлежность x Е Etoi(A, b) равносильна Tol (x; A,b) > 0, т. е. допусковое множество решений интервальной линейной системы Ах = b есть множество уровня функционала Tol
Stol(A, b) = {х Е R | Tol (х; А, Ь) > 0}. (16)
В тексте под "интервализацией" точечного математического объекта а понимается введение его нижней и верхней границ а и а: так что интервальный объект а понимается как а = [а, а] := {х Е R | а < х < а}. Иначе а = mid а + [— 1,1] • rad а. В качестве mid а в таком случае используется точечное значение а.
0.3. Дополнение к постановке задачи
Основанием для расчета являются данные плоского мозаичного детектора рентгеновского излучения в геометрии камеры-обскуры с шагом по времени порядка 1 мкс и данные по восстановлению границы плазмы с шагом по времени, равным 1 мс. Эти границы плазмы используются для постановки геометрической части задачи. Сначала производится разбиение расчетной области (5) на криволинейные элементы в соответствии с представлениями о структуре плазмы. Предполагается, что в рамках одного элемента Qk разбиения светимость постоянна. Для построения матрицы СЛАУ (8) вычисляются длины хорд наблюдения элементов детектора (пикселей). Элементом матрицы является длина Alkij отрезка хорды наблюдения L^ от пикселя в конкретном элементе Q.
Зарегистрированные в каждом пикселе детектора величины излучения формируют правую часть СЛАУ. Данные детектора снимаются гораздо чаще, чем реконструируется геометрия плазмы. В связи с этим можно представить правую часть интервальным образом, задав нижнюю и верхнюю границы сигналов в каждом пикселе на небольшом отрезке времени. Таким образом, для текущей конфигурации плазмы и конкретного момента времени формировалась интервальная система линейных алгебраических уравнений (12).
Алгебраическая часть задачи состоит в получении ИСЛАУ и в случае ее неразрешимости — регуляризации задачи. Решение отыскивается методами линейного программирования как точка, требующая наименьшего уширения интервалов неопределенности правой части ИСЛАУ по норме 11 при ограничениях неотрицательности решения (см. подразд. 3.4). Этот вид регуляризации дает гарантированную разрешимость ИСЛАУ.
Информация о пространственном положении плазмы, сепаратрисы и магнитной оси, данные детектора зависят от момента времени разряда. Приведем последовательность действий при получении и исследовании решения для единичного момента разряда. Для изучения эволюции светимости объема плазмы в течение разряда выбирается часть длительности разряда, из базы данных извлекаются входные данные и производятся необходимые вычисления для единичного момента разряда.
1. Входные данные
Входные данные для расчета:
1. Геометрическая информация о конструктивных элементах: чертежи вакуумной камеры токамака, устройство и расположение матричного детектора — взаимное расположение входной диафрагмы и плоскости, в которой расположены чувствительные элементы детектора (далее — пиксели), размеры и расположение пикселей.
2. Расчетные координаты геометрии плазмы: области положения границы плазмы и координаты точки, задающей положение "магнитной оси" — условной оси вращения плазмы в полоидальном направлении.
Входные данные делятся на общие геометрические данные, которые определяются устройством токамака "Глобус-М", и данные, относящиеся к конкретному разряду и моментам времени внутри него.
1.1. Общие геометрические данные
Общая геометрическая информация: чертежи вакуумной камеры токамака, устройство и расположение матричного детектора, размеры и расположение чувствительных элементов детектора — была неизменна с момента запуска токамака "Глобус-М" до начала его реконструкции в 2018 г.
Классические токамаки имеют вид открытого тора с большим отношением большого радиуса к малому. Иначе говоря, они имеют большое "горло". В конце XX в. получила популярность идея более компактной установки, по форме близкой к сфере, — сферо-мак. "Глобус-М" [1] ФТИ им. А.Ф. Иоффе РАН — первый российский токамак такого типа. Всего в мире около десяти установок вида сферомак. На рис. 1 представлена пространственная модель вакуумной камеры токамака "Глобус-М". Форма токамака в основном формируется пересечением сферы и двух цилиндров. Внутренний цилиндр —
Рис. 1. Пространственная модель вакуумной камеры токамака "Глобус-М" и криволинейный элемент разбиения объема плазмы
Fig. 1. The spatial model of the vacuum chamber of the "Globus-M" tokamak and a curved element for dividing the plasma volume
токопровод, внешний цилиндр служит для размещения оборудования. На рис. 1 оранжевым цветом представлен элемент разбиения объема токамака на элементы с симметрией по азимутальному углу.
1.2. Матричный детектор мягкого рентгеновского излучения
Для характеризации радиационных потерь плазмы используется мягкое рентгеновское излучение частично ионизованных примесей, само же рабочее вещество плазмы (водород, дейтерий) полностью ионизовано и в этой области спектра излучать не может. Матричный детектор мягкого рентгеновского излучения подробно описан в [3]. Расположение детектора (вид сверху) представлено на рис. 2, по осям показаны размеры, иллюстрирующие геометрию токамака и детектора.
На рис. 2, а показаны внешний диаметр вакуумной камеры, центральный соленоид (серый цвет), хорды наблюдения — пунктир. Детектор расположен за пределами вакуумной камеры и в масштабе левого рисунка очень мал. На рис. 2, б представлены элементы матричного детектора. Границы конструктива показаны окружностями, чувствительные элементы (пиксели) — квадратами. Пиксели лежат в одной плоскости, нормаль которой направлена в центральную область вакуумной камеры. Хорды наблюдения для пикселей проходят через апертуру камеры-обскуры, которая показана окружностью.
Рис. 2. Токамак "Глобус-М": a — вид сверху; б — детектор и хорды наблюдения: С — апертура детектора, плоскости P1, P8, P16 — сечения, используемые для построения матрицы длин хорд Fig. 2. Tokamak "Globus-M": a — view from above; б — detector and observation chords: С — detector aperture, planes P1, P8, P16 — sections used to construct a matrix of chord lengths
a
Если центр координат (0, 0) расположить в геометрическом центре токамака, то координаты проекций крайних точек детектора на экваториальную плоскость составят соответственно (0, -0.708) и (0.0288, -0.719), размеры в метрах с точностью до 1 мм. Координаты апертуры камеры-обскуры (0.0290, -0.677, 0).
Хорды наблюдений различных пикселей неравнозначны ввиду того, что соленоид непрозрачен для излучения. Относительно малая часть хорд (48 из 256) заканчивается на центральном столбе, остальные — на стенке вакуумного объема. Этот факт учитывается в подразд. 2.7 при вычислении элементов геометрической матрицы.
Детектор находится в экваториальной плоскости, область наблюдения симметрична относительно экватора. На рис. 3, а показано расположение детектора относительно камеры токамака в меридиональном сечении, на рис. 3, б — область детектора. Координаты апертуры камеры-обскуры в системе координат (R, Z) составляют (-0.677, 0). Также здесь штриховыми линями показаны хорды наблюдения в меридиональном сечении. Хорды наблюдения образуют углы до 0.33 рад (18.9 град.) относительно экваториальной плоскости. На рис. 4 представлено расположение пикселей в плоскости детектора. Они образуют структуру, близкую к регулярной по горизонтали, и практически регулярную по вертикали.
Данные детектора хранятся в базе данных токамака "Глобус-М" в виде трехмерного массива:
b = b(i,j,t). (17)
Рис. 3. Детектор и хорды наблюдения (а), меридиональная плоскость (б) Fig. 3. Detector and observation chords (а), meridional plane (б)
□ □ [ □ □ [ 3 □ □ D □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ E □ HE
Ы Ы 1 □ □ [ □ □ [ □ □ [ □ □ [ □ □ [ d Ш Ш D □ □ D □ □ D □ □ D □ □ D □ □ LLI U U □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ Ы Ы □ □ □ □ □ □ □ □ □ □ Ы Ы Ы □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ EJ EJ E □ HE □ HE □ HE □ HE □ HE
□ □ [ □ □ [ □ □ [ □ □ [ □ □ [ П П [ D □ □ D □ □ D □ □ D □ □ 3 □ □ Л П П □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ П П П □ □ □ □ □ □ □ □ □ □ П n □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ П П n □ HE □ HE □ HE □ □ E El El E п п г
□ □ [ □ □ [ D □ □ D □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ □ HE □ HE
О 0.005 0.01 0.015 0.02 0.025 0.03
у, m
Рис. 4. Детектор мягкого рентгеновского излучения. 16 х 16 пикселей размером 0.88 х 1.23 мм Fig. 4. Soft X-ray detector. 16 х 16 pixels in size 0.88 х 1.23 mm
Здесь i,j = 1, 2,..., 16 — индексы пикселей, t — индекс момента времени. Размер массива t ~ 105 для каждого разряда.
1.3. Расчетная форма границы плазмы в токамаке
В разделе физики плазмы, относящемся к магнитным ловушкам, считается, что внешняя граница плазмы соответствует последней замкнутой поверхности потока магнитного поля, удерживающего ее. Эту границу называют сепаратрисой, на ней давление газа плазмы уравновешено давлением магнитного поля. За пределами сепаратрисы плазмы нет, но есть остаточное давление газа, в котором вещество находится не в состоянии плазмы.
На рис. 5 показаны различные сечения токамака "Глобус-М" [1]: стенки вакуумной камеры — сплошная линия, сепаратрисы (границы плазмы) и магнитной оси (условный центр тока плазмы) — прямоугольники.
На установках типа токамак принято хранить данные о геометрическом месте сепаратрисы в так называемых g-файлах, являющихся результатом работы программы EFIT (Equilibrium Fitting) [9] реконструкции поверхностей магнитного потока. В ^-файлах содержится различная информация: RBDRY, ZBDRY — (R, Z)-координаты точек сепаратрисы, NDRY — число точек сепаратрисы, rdim, zdim — размер сетки по радиусу и вертикали, R, Z — координаты точек сетки по радиусу и вертикали, flux — двумерный массив магнитного потока.
0.4
0.2
I
□
I
/
-0.4
-0.2
0
0.2
0.4
0.6
0.8
R, m
Рис. 5. Форма плазмы в токамаке. Сплошная линия — вакуумная камера, прямоугольники — граница плазмы, рассчитанная на прямоугольной сетке EFIT
Fig. 5. Plasma form in tokamak. Solid line — vacuum chamber, rectangles — plasma boundary calculated on a rectangular grid
В окрестности минимума распределения flux находится "магнитная ось". Область "магнитной оси" можно считать "центром" плазмы, вокруг которой происходит вращение плазмы в полоидальном направлении.
2. Детали геометрических вычислений
Форма поверхностей плазмы сферического токамака "Глобус-М" существенно отлична от "классической" тороидальной. Строго говоря, даже равновесная поверхность не имеет тороидальной (азимутальной) симметрии. Вместе с тем степень асимметрии вдоль азимутального угла не очень велика: неоднородность магнитного поля в этом направлении можно из геометрических соображений оценить в пределах 10 %. Гораздо большую роль может играть асимметрия между внешней и внутренней поверхностями. Для полноценной трехмерной постановки задачи необходимо соотнести степень дискретизации и количество уравнений, равное 256. При разбиении токамака на 8-16 частей по азимуту остается 16-32 степени свободы в меридиональном сечении. В предположении равноценности направлений z и г имеем, в лучшем случае, возможность разбиения на восемь частей с размером порядка 10 см, что выглядит весьма грубым представлением для токамака "Глобус-М".
В связи с этим на данном этапе задача ставится в двумерной постановке. Объем плазмы мыслится как набор вложенных поверхностей постоянного давления, на которых свойства плазмы постоянны.
2.1. Проекция окружности на плоскость детектора
Центральные проекции круга на плоскость в последние годы тщательно изучаются применительно к задачам компьютерного зрения. Например, в работе R. Matsuoka и S. Maruyama [10] исследован вопрос об искажении эллиптической формы проекции круга в общей геометрии расположения проецируемого круга S0, центра проекции точки С и плоскости Р. Причина искажения — в разном удалении различных точек круга S0 от центра проекции и вследствие этого в разном увеличении / для его разных точек.
На рис. 6 показано, что проекция центра проецируемого круга и центр эллипса не совпадают. В публикации [10] представлены формулы, позволяющие вычислить величину смещения. Теоретически можно было бы вычислить проекции для всего набора окружностей, заполняющих объем токамака, и на основании этих данных построить геометрическую модель задачи. В современных компьютерных пакетах, например OpenGL, подобные расчеты выполняются с использованием проекционных матриц [11].
В настоящей статье мы не используем такой подход по нескольким причинам. Вычисления по формулам в духе [10, 11] весьма громоздки, что плохо с учетом неизбежного распространения ошибок исходных данных. Нам представляется более простым вычисление средствами векторной алгебры.
На рис. 7 приведена центральная проекция окружности, лежащей в плоскости, ортогональной оси токамака, на плоскость детектора. Такими окружностями равной светимости можно заполнить весь объем токамака (см. рис. 1). Здесь О — геометрический центр токамака, S — окружность нулевой толщины, С — входная диафрагма камеры-обскуры, Р — плоскость детектора, Sp — проекция окружности на плоскость детектора.
е: eccentricity of projection
Рис. 6. Проекция круга на плоскость при центральной проекции (рис. 1 из [10]) Fig. 6. The projection of a circle onto a plane with a central projection (Fig. 1 from [10])
S(A с - p
r Ö 1 d\ \..... Si ; p
к С Sp
Рис. 7. Расположение детектора относительно токамака, пространственное представление и вид сверху. Проекция круглой области токамака на плоскость детектора Fig. 7. The location of the detector relative to the tokamak, spatial representation and top view. The projection of the circular tokamak area onto the detector plane
Опишем движение точки М по окружности и получим множество точек S для центральной проекции. Обозначения соответствуют рис. 7.
OM = OA + г(р), (18)
г(р) = R0z х M. (19)
Здесь OM — вектор, описывающий точки окружности S; OA — центр окружности, Roz — матрица вращения точки М вокруг оси Oz.
Для каждой точки М найдем вектор n« единичной длины, проходящий через нее и апертуру детектора С. Зададим уравнение точек прямой как
MC = M + t • nM. (20)
Здесь t — параметр, описывающий удаление точки прямой от точки М.
Вычислим координаты пересечения прямой MC, проходящей через точки М и С, и плоскости детектора Р(N, P). Обозначения плоскости Р(N, P) соответствуют плоскости с нормалью N = (Nx,Ny,NZ), проходящей через точку P = (Рх,Ру,PZ). Уравнение плоскости
N • (R - P) = 0. (21)
Здесь R = (х, у, z).
Приравнивание R координате точек прямой MC в уравнении плоскости (21) позволяет вычислить значение параметра t:
t = n^P-M). (22)
n„•N v '
Окончательно подстановка значения t из формулы (22) в выражение (20) дает координаты проекции окружности в плоскости детектора.
На рис. 8 F — проекция наиболее удаленной от С точки окружности в токамаке, N — наиболее близкой; L,R — точки окружности с координатами х = 0 и у = ±R.
Из набора окружностей можно получить разбиение объема токамака на непересекающиеся поверхности, задав в качестве вершин множества точки сетки меридионального сечения.
Геометрический смысл величины параметра abs(i) и его поведение в зависимости от угла ip могут быть поняты из построения (см. рис. 7). Размах изменения abs(t) соответствует соотношению большой и малой полуосей эллипса, являющегося проекцией круга на плоскость детектора. Эксцентриситет m определяется в первую очередь соотношением расстояний \\NC|| и \\FC|| между диафрагмой детектора и ближайшей и наиболее удаленной к ней точками окружности S. С учетом малости размеров детектора по отношению к размеру токамака имеем
^ у/(г - Хс)2 + Н2 m & . (23)
у/(г + Хс)2 + Н2 v 7
Здесь Н = \\ОА\\, г = \\АМ\\, Хс = \\ОС\\ (см. рис. 7). Формула (23) носит иллюстративный характер и пригодна в случае направления нормали детектора к центру токамака. Характерный вид графика параметра t приведен на рис. 9. Для "границы" плазмы, показанной на рис. 5, график параметра m представлен на рис. 10.
Рис. 8. Вид токамака сверху (a): F,L,R,N — характерные точки окружности, С— апертура детектора; проекция круглой области токамака на плоскость детектора: FN, LR — проекции диаметров окружности (б)
Fig. 8. Top view of a tokamak (a): F,L,R,N — characteristic points of a circle, С — detector aperture; projection of the circular tokamak area onto the detector plane: FN, LR — projection of circle diameters (б)
Рис. 9. График параметра Формула (22). Рис. 10. График параметра т для "границы"
Обозначения точек как на рис. 8
плазмы, показанной на рис. 5. Формула (23)
Fig. 9. The graph of the parameter t as for Fig. 10. Parameter graph m for the "boundary" formula (22). Notation for points are the same of the plasma shown in fig. 5. Formula (23) as in fig. 8
a
Максимальное значение abs t достигается на большой полуоси эллипса. При повороте направления нормали по отношению к оси токамака эксцентриситет эллипса увеличивается и в конечном счете график переходит через параболу в гиперболу. При этом abs t неограниченно возрастает.
В подразд. 3.2 обсудим характер фигур, представленных на рис. 8.
2.2. Влияние конечных размеров светящегося кольца на вид его центральной проекции на плоскость
Рассмотрим влияние конечных размеров светящегося кольца на вид его центральной проекции на плоскость. В таком случае кольцо будет открытым круглым тором. Используем интервальную арифметику согласно формулам (10). Пусть светящееся кольцо имеет радиус R =0.3 ми расположено на расстоянии Z = 0.1 м выше экваториальной плоскости. Радиус кольца положим равным г = 0.005 м.
В формуле (19) положим радиус Roz интервальной величиной R в классической интервальной арифметике:
R = [R, R]. (24)
Здесь R, R — соответственно нижняя и верхняя границы R:
R = mid R — rad R, R = mid R + rad R, mid R = R.
Конкретно
R = 0.295, R = 0.305, mid R =0.3.
Применим формулы (20) и (22) для построения проекции "объемной" окружности (круглого тора) на плоскость. Технически использование интервальной арифметики (10) производилось в системе Octave с применением библиотеки interval [12].
На рис. 11 представлена проекция на плоскость кольца круглого открытого тора с большим и малым радиусами 0.3 и 0.005 м соответственно. Характерной особенностью проекции является различная передача близкой и удаленной от центра проекций точек круглого тора. Удаленные области имеют малый угловой размер видимости, а ближние — существенно больше. Соответственно также различаются и их проекции. Если светимость тора постоянна на его теле, то плотности освещенности проецируемых областей будут находиться в обратном отношении: проекции удаленных областей будут ярче ближних.
Приведем пример. Зададимся конкретными параметрами кольца круглого открытого тора с большим и малым радиусами 0.3 и 0.005 м соответственно. Вычислим t по формулам (20) и (22) как интервальную величину.
На рис. 12 представлены графики параметров t и rad t. Для интервальной величины t изображены sup t и inf t. Что касается радиуса параметра rad t, то его величина изменяется почти на порядок, отношение максимального радиуса к минимальному составляет ^ 8.6. Если же неопределенность имеет не радиус, а высоту исходного кольца, то интервализации при расчетах подвергается координата Z. При rad Z = 0.005 м в формуле (19) получаем фигуру на рис. 13.
Характерной особенностью проекции такого типа является разность в "размывании" вдоль оси Z для близкой и удаленной от центра проекций точек цилиндра. Элементы разбиения на рис. 16 (см. подразд. 2.6) имеют протяженность как по радиальной координате, так и по оси Z, так что их проекции имеют черты фигур как на рис. 11, так и рис. 13. В подразд. 3.2 на рис. 25 представлена подобная структура.
Рис. 11. Проекция круглой области токамака на плоскость детектора с учетом конечности радиуса. Формулы (20) и (22). Система Octave и библиотека interval
Fig. 11. The projection of the circular tokamak region onto the plane of the detector, taking into account the finiteness of the radius. Formulas (20) and (22). Octave system and Interval library
Рис. 12. Графики параметров t (сверху) и rad t (снизу). Система Octave и библиотека interval Fig. 12. Graphs of parameters t (above) and rad t (below). Octave system and Interval library
Пусть интересующая нас структура токамака образует в нем окружность в вертикальной плоскости. Анализ таких структур необходим для понимания временной эволюции светимости, например движения светящегося объема в полоидальном направлении. В плоскости детектора проекция такой области представлена на рис. 14.
Рис. 13. Проекция цилиндрической круглой области токамака на плоскость детектора Fig. 13. The projection of the cylindrical circular area of the tokamak on the plane of the detector
Рис. 14. Проекция круглой области токамака в вертикальной плоскости на плоскость детектора
Fig. 14. The projection of the circular tokamak area in the vertical plane onto the plane of the detector
Основная черта проекции, как и в предыдущих примерах, — большая разница в радиусах точек проекций для более близких к детектору областей. Это приводит к проекции ближней на большее число пикселей, чем дальней, и соответственно, к меньшей яркости при прочих равных условиях (см. ниже рис. 23 и 24).
Итак, наблюдаемое проявление конечной величины светящихся объектов при центральной проекции в случае проекции объема на плоскость — значительная разница в яркости проекций близких и удаленных объектов. В сочетании с существенно более высоким перспективным искажением более близких объектов это приводит к очень значительной разнице между отображением идентичных структур в детекторе в зависимости от их расположения по отношению к его апертуре, как видно из рис. 22.
Рассмотренные примеры иллюстрируют различное размывание изображений границ проекций цилиндрической и круглой областей токамака (рис. 13 и 14). Причиной размывания изображений является увеличение угла наблюдения структуры по мере приближения к апертуре камеры-обскуры, при этом выполнение расчетов в интервальной арифметике представляет естественный и удобный способ описания таких явлений.
2.3. Разбиение меридионального сечения на области равной светимости (построение сетки)
Зададимся вопросом о типе разбиения расчетной области (сетки). Это первый шаг в процессе нахождения решений. Полезны следующие свойства сетки:
• адаптивность к форме сепаратрисы и области магнитной оси;
• легкость построения и изменения для массовых расчетов, например для временной эволюции;
• легкость объединения групп элементов сетки (ячеек) и их измельчения по мере необходимости подобно тому, как это делается при расчетах по методу конечных элементов и при обработке видеоизображений;
• возможность простой параметрической модификации для привлечения данных
других диагностик и модельных представлений. Должны ли элементы расчетной области описываться аналитически или же их вычислительно выгодней задавать алгоритмически? В данной статье мы рассматриваем все поверхности как расчетные множества точек без попыток привлечения аналитических описаний.
2.4. Построение разбиения области пространства между замкнутыми кривыми и радиальное разбиение
Для каждого сечения плазмы плоскостью необходимо провести его разбиение на элементы, в которых светимость полагается постоянной. Формализуем задачу. Пусть заданы множество, ограничивающее область пространства Si, и другая область внутри нее So. Задача: построить разбиение пространства множеством кривых S2, являющихся "средними" поверхностями обоих множеств, с выделением опорных направлений.
В g-файлах, являющихся результатом работы программы EFIT (Equilibrium Fitting) [9] для реконструкции поверхностей магнитного потока, содержится двумерный массив магнитного потока Ф. Задаваясь набором констант С
Ф = С, (25)
получаем набор вложенных поверхностей с возрастающим расстоянием от магнитной оси до сепаратрисы. Далее необходимо выбрать радиальные направления. Характерная форма сепаратрисы, представленная на рис. 5, свидетельствует о наличии особых
0.6
-0.6
-0.6 -0.4 -0.2 0 0.2 0.4 0.6
Г, m
Рис. 15. Оечение объема плазмы в меридиональном сечении и система областей постоянной светимости
Fig. 15. The plasma volume cross section in the meridional section and the system of regions of constant luminosity
направлений в меридиональном сечении токамака, причиной которых является текущая конфигурация магнитных полей. Объективную информацию об этих направлениях дают дифференциальные свойства сепаратрисы.
Теорема о четырех вершинах [13] утверждает, что функция кривизны простой замкнутой гладкой плоской кривой имеет четыре локальных экстремума (в частности, по меньшей мере два локальных максимума и два локальных минимума). Название теоремы отражает соглашение называть экстремальные точки функции кривизны вершинами.
Опираясь на информацию о кривизне поверхности, можно построить вложенные поверхности следующим образом. Пусть найдены четыре особые точки. Соединив их с центром, получим четыре сектора. Внутри каждого сектора проводим дополнительное разбиение. Такой тип сетки похож на паутину. Подобное построение используется в программе Spider [14], разработанной для моделирования равновесия и эволюции плазмы токамака.
На рис. 15 представлен пример разбиения плазмы в меридиональной плоскости. Первые шесть внутренних замкнутых областей получены из условия (25). Помимо этого добавлена еще одна поверхность за сепаратрисой для учета светимости примесей вне объема плазмы вблизи сепаратрисы. По форме эта поверхность подобна сепаратрисе. Радиальных разбиений каждой вложенной поверхности 26, так что всего элементов сетки 26 х 8 = 208.
В немеридиональных сечениях объема токамака плоскостями расчет сетки производится c сечений исходного разбиения плоскостями, как показано в подразд. 2.5.
2.5. Сечения аксиально симметричной фигуры плоскостями
Постановка задачи сечения аксиально симметричной фигуры плоскостями: необходимо спроецировать двумерную фигуру S, заданную в координатах X и Z, на плоскость х = Н, где Н — расстояние от центра токамака до плоскости сечения.
Для реализации этого проецирования воспользуемся формулами поворота вокруг оси Oz в предположении, что все точки фигуры (х, 0, z) совершают поворот на один и тот же угол а. Координаты точек (x',y',z') фигуры Sout после поворота задаются следующим образом:
Здесь а — угол поворота точки поверхности с координатой х в плоскость, где х' = Н. Таким образом, получаем плоские фигуры S^. Предположим, что множество точек, построенных по формуле (26), таково, что
х • cos(a) — у • sin(a) х • sin(a) + у • sin(a)
z
\
(26)
minх > 0,
$out
(27)
иначе говоря, фигура вращения, образованная вращением вокруг оси Z, содержит "горло". При этом имеются несколько вариантов связности сечений такой фигуры вращения плоскостью х = Н. Если расстояние Н от оси Z меньше либо равно наименьшей
координате по х фигуры, то сечение является двусвязной фигурой. Если Н больше минимальной и меньше либо равно максимальной координатам, то сечение — односвязная фигура. Третий вариант имеет место, когда Н больше максимальной координаты по х, в этом случае секущая плоскость не пересекает фигуру.
Сечения круглого тора хорошо изучены и описываются аналитически [15]. Эти кривые (овалы Кассини, леминиската Бернулли) использовались как тестовые примеры для отладки и тестирования алгоритмов.
2.6. Хорды наблюдения для пикселей детектора
Финальным этапом геометрических расчетов является построение хорд наблюдения Ьг для каждого пикселя детектора мягкого рентгеновского излучения и вычисление длины хорды Д/у в каждой из областей разбиения .
Это наиболее трудоемкая часть вычислений, которая влечет и основные возможные неточности в дальнейших вычислениях. Причиной этих неточностей являются неопределенности в исходных данных как формы сепаратрисы, вычисленной на прямоугольной сетке, так и геометрии детектора: положения его входной диафрагмы, наклона к экваториальной плоскости и оси камеры токамака.
По постановке задачи искомые переменные пропорциональны длинам хорд в объемах, на которые разбивается плазма. Ниже представлен алгоритм нахождения длин хорд с учетом нахождения групп по 16 пикселей детектора в одной плоскости.
Примеры хорд наблюдения в плоскости 8, которой соответствует расстояние от центра токамака Н = 0.295, представлены на рис. 16. Из всего разбиения сечения то-камака плоскостью 8 показаны сечение вакуумной камеры и элемент разбиения с номером 161.
На рис. 17 приведены два фрагмента рис. 16, показывающие пересечения элемента разбиения с номером 161 хордами наблюдения. Область ячейки с номером 161 в данной плоскости пересекается четырьмя хордами в части токамака, ближней к детектору, и двумя хордами в части токамака, дальней от детектора.
Алгоритм нахождения матрицы Ь^ длин хорд
Вход
Форма сепаратрисы — массив Б = (КВВЕУ, ЕВВВУ) 2Б-распределение магнитного потока /1их(К, Z) Выход
Матрица длин хорд Ь^ в областях равной светимости х^ Алгоритм
Разбиение меридионального сечения 5 на области ] = 1,..., 208 равной светимости (построение сетки).
Нахождения плоскостей Рг, г = 1,..., 16. Построение сечений Si, г = 1,..., 16.
Построение набора хорд ] = 1,..., 16 для пикселей детектора в каждом сечении Sí. Нахождение длин хорд Ь^ в каждой из областей ].
m х п=256 х 208 section=8 Н= 0.295 cell=161
b
......¿(J -1 - - - -" г----- -- ----
г— -- г - - -
\ 1 - - k.y-
i i
-0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6
Г, m
Рис. 16. Хорды наблюдения в плоскости Н = 0.295 и 161-й элемент разбиения
Fig. 16. Chords of observation in the plane H = 0.295 and the 161st element of the partition
Рис. 17. Увеличенные области расположения 161-го элемента разбиения и его пересечения хордами наблюдения
Fig. 17. Enlarged areas of the 161st element of the partition and its intersection by observation chords
2.7. Матрица СЛАУ
Вычисленные хорды наблюдений для отдельных пикселей позволяют построить матрицу СЛАУ. В общем случае пространственного расположения детектора относительно объема плазмы матрица длин хорд в элементах не имеет симметрий. Однако в случае установки детектора параллельно оси токамака элементы детектора расположены вертикально и в каждом из 16 столбцов 16 хорд наблюдения лежат в одной плоскости, как на рис. 4. В таком случае структуру матрицы длин хорд можно представить в блочном виде:
Ь =(¿1 ,¿2 ,...,¿16 )Т. (28)
Всего число блоков Ь^ равно 16 — по числу столбцов детектора. Число хорд наблюдения равно тоже 16 по числу строк детектора. В развернутом виде блоки Ь^ в формуле (28) имеют вид
Сечение
Li —
L16 —
16
Хорда 1 2
15
16
Слой 1
«1,1 «2,1
«15,1 \«16,1
«2,/
«15,/ «16,/
. ^ Слой 8
«1,п-/+1 - -«2,п-/+1 - -
«15,п-/+1 - -«16,п-/+1 - -
• «1,п
• «2,п
- «15,п
- «16,п)
241
242
255
256
«241,1 - - - «241,/ - - - «241,п-/+1 - - - «241 ,п
«242,1 - - - «242,/ - - - «242,п-/+1 - - - «42,2п
«255,1 - - - «255,/ - - - «255,п-/+1 - - - «255,п
\«256,1 - - - «256,/ - - - «256,п-/+1 - - - «256,п^
(29)
Х1
Xi
Хп-1+1
Рис. 18. Структура матрицы длин хорд: точки — ненулевые элементы матрицы, линии — плоскости сечения и слои разбиения
Fig. 18. The structure of the matrix of lengths of chords: points are nonzero elements of the matrix, lines are plane sections and layers of partition
1
Рис. 19. Структура матрицы длин хорд, часть рис. 18. Показаны одна плоскость сечения 8 и один слой разбиения 7; квадраты — ненулевые элементы матрицы
Fig. 19. The structure of the matrix for lengths of chords, part of fig. 18. One plane of section 8 and one layer of partition 7 are shown; squares are nonzero matrix elements
На рис. 18 по горизонтали штриховыми линиями выделены блоки 16 плоскостей, по вертикали - 8 слоев по 26 ячеек. Точками показаны ненулевые элементы матрицы независимо от их величины. Нумерация ячеек идет по сворачивающейся спирали, содержащей 8 витков по 26 ячеек. Нумерация плоскостей, в каждой из которых находится по 16 хорд наблюдения, идет от периферии токамака с большими расстояниями Н между плоскостью сечения и его центром до малых расстояний. Например, для сечения 8, показанного на рис. 16, сечение будет содержать шесть ненулевых элементов для столбца, соответствующего элементу разбиения 161, с искомой светимостью х161. Данный элемент находится в слое 7 от внешней границы сетки и на пятом радиусе.
На рис. 19 видна особенность структуры матрицы для центральной проекции с центром в экваториальном сечении токамака: половина хорд наблюдения идет выше экваториального сечения токамака, половина — ниже, как это показано на рис. 3.
3. Исследование геометрической природы задачи и модельные вычислительные эксперименты
Данные на детекторе представляют собой центральную проекцию тороида на плоскость детектора. Рассмотрим, как проецируется светимость объема токамака с формой, близкой к сферической, на плоскость детектора.
3.1. Макроскопические структуры
Исследование модельных прямых задач имеет целью понять возможные типы решений и научиться строить алгоритмы их распознавания по данным детектора. Наиболее естественная пара структур: ядро и оболочка. По аналогии с астрономической тер-
минологией свечение оболочки называют также "коронарным". Сначала рассмотрим проекцию свечения небольшой центральной области плазмы токамака на плоскость детектора — "изображение" этого объема. На рис. 20 показаны распределение светимости в меридиональном сечении токамака и ее проекция на плоскость детектора.
Рис. 20. Проекция на детектор — свечение "ядра"
Fig. 20. The projection onto the detector — the luminosity of the "core"
Рис. 21. Проекция на детектор — свечение второго слоя внутреннего объема
Fig. 21. The projection onto the detector — the luminosity of the second layer of the internal volume
Рис. 22. Проекция на детектор — свечение последнего слоя внутреннего объема
Fig. 22. The projection onto the detector — the luminosity of the final layer of the internal volume
По аналогии с применяемым в оптике и обработке изображений термином "функция распределения точки" (point spread function, PSF) изображение объема на детекторе можно назвать функцией распределения объема (volume spread function, VSF).
Каковы характерные особенности "изображения" ядра плазмы? Распределение яркости симметрично относительно проекции экваториального сечения, что является следствием расположения детектора и нормали к его входной диафрагме в этой плоскости. Максимум яркости смещен относительно геометрического центра детектора (столбец 12, строки 8 и 9), в целом распределение яркости проекции имеет тип "комета". Оно соответствует изображению круглого тора относительно небольшого диаметра. На рис. 20 также можно видеть ослабление яркости пикселей для столбцов 1-3. Причина — геометрическая тень центрального соленоида, расположенного в геометрическом центре токамака, как показано на рис. 1.
Если вместо центрального объема ("ядра") рассмотреть последовательность концентрических объемов, то их изображение будет представлять "расплывание" изображения по плоскости детектора с потерей светимости в месте первичного положения максимума.
На рис. 21 и 22 показаны распределение светимости в меридиональном сечении то-камака и ее проекция на плоскость детектора для второго и последнего (см. рис. 5) слоев внутреннего (ограниченного сепаратрисой) объема соответственно. Поверхность
Рис. 23. Проекция на детектор — свечение внешней части сепаратрисы
Fig. 23. The projection onto the detector — the luminosity of the outer part of the separatrix
Рис. 24. Проекция на детектор — свечение внутренней части сепаратрисы
Fig. 23. The projection onto the detector — the luminosity of the inner part of the separatrix
максимальной яркости теперь окружает геометрическое место максимума от самого внутреннего слоя и приобретает тип "пустотелой кометы". Можно увидеть на рис. 21 увеличенное изображение части рис. 20.
При дальнейшем удалении светящегося объема от "центра" тороида (магнитной оси) происходит продолжающееся расплывание изображения с его одновременным увеличением с центром в изображении магнитной оси. На всех изображениях проявляется затенение области проекции центрального соленоида (столбцы 1-3).
Рассмотрим теперь "коронарный" тип светимости, выделив область небольшой "толщины" непосредственно за сепаратрисой. Выделим в этом слое внешнюю и внутреннюю части. Изображение светимости внешней и внутренней частей сепаратрисы представлено на рис. 23 и 24. Поскольку внешнюю часть сепаратрисы пересекают все хорды наблюдения, проекция ее свечения не имеет структуры, первые три столбца затенены. В противоположность внешней части изображение внутренней части имеет четкую геометрическую границу за четвертым столбцом матрицы детектора.
Из-за практически порядкового соотношения областей освещенности проекций внешней и внутренней частей сепаратрисы в изображении проявляется прежде всего внутренняя ее часть, а внешняя выглядит как фон. На рис. 14 приведен иллюстративный пример, поясняющий этот факт.
3.2. Микроскопические структуры
Рассмотрим теперь вопрос о том, как проявляется свечение областей малого размера. Будем пользоваться результатами, представленными в подразд. 2.1. Пусть вся область плазмы не излучает, за исключением одной ячейки. Для определенности возьмем ячейку, расположенную в седьмом слое (от стенок камеры) и на девятом радиусе в разбиении по углу. Для краткости такое расположение в меридиональной плоскости будем обозначать (layer, radial direction), в данном примере (7, 9).
Изображение свечения этой ячейки на детекторе показано на рис. 25. Изображение четырехугольного кольцевого объекта развернуто в фигуру, образованную суммой проекций эллипсов, которые являются результом проекций окружностей согласно формулам (20) и (22) — рис. 8. Яркость каждого пикселя определяется попадающей в него долей светимости объема. Зависит от длины хорд наблюдений в выбранном объеме, т. е.
Рис. 25. Проекция на детектор — свечение ячейки (7, 9)
Fig. 25. Projection onto the detector — luminosity of the cell (7, 9)
не только от факта нахождения светящегося объема на пути хорды наблюдения, но и от длины хорды в этом объеме. Помимо этого различные участки светящейся поверхности дают различные проекции на плоскость детектора из-за перспективного искажения, что иллюстрируют рис. 11 и 13. Эти факторы учтены при построении матрицы (29).
3.3. Геометрические структуры для тестирования решений
Для тестирования техники решения полезно сформировать набор задач, на которых можно установить корректность получаемых результатов. Первые два примера: "пирамида" и "спираль" — относятся к регулярным гладким решениям. "Пирамида" — монотонное изменение светимости от слоя к слою (рис. 26). Данное распределение моделирует ситуацию с преобладающим свечением ядра, убывающим по слоям, с резкими границами светимости между слоями. Распределение светимости в проекции на детектор практически центрально симметрично, с затенением от соленоида токамака в левой части.
"Спираль" — монотонное изменение светимости от ячейки к ячейке по номеру (рис. 27). Данное модельное распределение является вариантом распределения "пирамида" с большим количеством уровней светимости, равным числу ячеек в разбиении расчетной области. Распределение светимости в проекциии на детектор практически
Рис. 26. Проекция на детектор — свечение типа "пирамида"
Fig. 26. Projection onto the detector — luminosity of the "pyramid" type
Рис. 27. Проекция на детектор — свечение типа "спираль"
Fig. 27. Projection onto the detector — luminosity of the "spiral" type
Рис. 28. Проекция на детектор — свечение типа "лес"
Fig. 28. Projection onto the detector — luminosity of the "wood" type
смещено вверх, связано с порядком обхода ячеек и затенением от соленоида токамака в левой части.
"Лес" — случайно расположенные в меридиональном сечении токамака светящиеся ячейки (рис. 28). Этот тип моделирует худший для разрешимости случай распределения светимости.
Наиболее проблемной из представленных моделей с вычислительной точки зрения является задача со случайными данными. Ее важность состоит в выяснении того, насколько техника решения способна разрешить непредсказуемое распределение светящихся примесей в плазме.
3.4. Решение модельных задач
Данный подраздел представляет краткое изложение алгебраической части решения задачи и введен в настоящую публикацию для иллюстрации процесса тестирования полученных геометрических результатов.
Рассмотрим типовые задачи, предложенные в подразд. 3.3. С математической точки зрения эмиссионная томография представляет собой обратную задачу — нахождение светимости различных областей объема плазмы по данным матричного детектора. Решение задачи реконструкции светимости с учетом неопределенностей входных данных ведется методами интервального анализа с использованием метода распознающего функционала, упомянутого в подразд. 0.2. Пусть имеем ИСЛАУ
Ах = Ь, (30)
где А — матрица длин хорд в областях различной светимости, Ь — данные детектора, х — значения светимости. Правая часть формируется согласно предложенной в работе [16] процедуре "естественной интервализации". Применительно к задаче реконструкции светимости процедура состоит в нахождении интервала Ь(г):
b(k,t) =
min b, max b
teT(k) teT(k)
(31)
Здесь к = к (г,]) — индекс пикселя матрицы детектора — см. формулу (17) в одномерном представлении в соответствии со структурой матрицы (29); Т(к) — некоторая
окрестность времени, определение которой требует отдельного исследования, связанного с темпом изменчивости светимости на коротких промежутках времени. Эта изменчивость различна для разных разрядов и разных временных участков каждого разряда.
Для получения физически состоятельного решения необходимо как нахождение непустого допускового множества решений, так и получение неотрицательных компонент решения. Реализация метода распознающего функционала tolsovty [17] на основе г-алгоритма [18] решает задачу безусловной максимизации, что не гарантирует неотрицательности компонент решения. Исследования конкретных задач показали, что условие неотрицательности решения в подавляющем числе случаев не выполняется.
Решение ИСЛАУ получаем с помощью метода "ведущего центра" — строим как ¿1-центр информационного множества с помощью метода и программного обеспечения, разработанных С.И. Жилиным [19]. В этом подходе ищется решение задачи линейного программирования специального вида, ограничения в которой описывают подвижные границы сжимаемого в точку информационного множества и могут быть дополнены условием неотрицательности на компоненты вектора решений. Этот подход, по существу, производит 11 -регуляризацию задачи с ограничением на неотрицательность компонент вектора решения.
Постановка задачи линейного программирования выглядит следующим образом [19]. Для условия
А • х С Ь, (32)
где А — точечная матрица правой части ИСЛАУ, Ь — интервальный вектор правой части ИСЛАУ, ставится задача линейного программирования (ЗЛП):
п
inN Wi, (33)
7Л -^
min
x,w
i=1
mid bi — Wi • rad bi < (A • x)i < mid bi + Wi • rad bi, (34)
Xj > 0, j = 1,... ,m, (35)
Wi > 0, г = 1,... ,n. (36)
Здесь Wi — масштабирующие множители для радиусов правой части "по уравнениям". Идея введения масштабирующих множителей состоит в выяснении того, насколько для конкретного уравнения необходимо увеличить радиус интервала правой части для удовлетворения условия (34).
Технически процедура вычисления строится следующим образом. Пусть А — матрица размером п х т, что соответствует СЛАУ с п уравнениями и т неизвестными. Строятся матрица С ограничений ЗЛП и расширенный вектор г правой части для уравнения (32):
с f А —diag(rad &)\ ^ f mid & \ (37)
\—А —diag(rad b)J ' у—mid bj '
Здесь diag(f) — диагональная матрица со значениями на диагонали, равными компонентам вектора v, размер матрицы ограничений С — (т + п) х (т + п), размер вектора искомого решения т + п.
Таким образом, ставится задача минимизации суммы весов Wi при неопределенностях компонент правой части rad b. Веса Wi при начале вычислений устанавливаются равными 1, т.е. ищется решение, помещающееся в интервальной правой части. При необходимости Wi увеличиваются или уменьшаются, показывая своей величиной
количественную степень совместности (при < 1) или несовместности (при > 1) конкретного уравнения. По условию (34) имеем
Tol
ЗЛП
= 0,
(38)
и (согласно [7], подраздел 6.5 "Полное исследование разрешимости") допусковое множество Stoi (А, Ь) (16) решения ИСЛАУ (32) содержит единственную точку, на которой достигается условие (38).
На рис. 29 представлены графики для типа модельного распределения "лес".
Относительное количество восстановленных переменных, существенно отличающихся от модельных, составляет не более 5 %. Получение приемлемого решения для случайного распределения светимости дает основание полагать возможность качественного решения и для других задач.
Анализ величин компонент вектора w в методе "ведущего центра" дает инструмент анализа качества построения решаемой ИСЛАУ. Наличие значимого количества компонент Wi > 1 дает основание полагать наличие неточности в определении элементов матрицы и необходимости уточнения процедуры ее построения, начиная с исходной посылки: определения положения границы плазмы — сепаратрисы.
В качестве меры вариабельности решения ИСЛАУ в [20] предложена следующая конструкция ive (interval variability of the estimate):
ive(A; b) = \fn ( min cond2 A ) • max Tol •
W A J R"
II arg тахЖп Tol II
(39)
В выражении (39) минимум числа обусловленности min^ ^ cond2 А находится среди всех реализаций матрицы ИСЛАУ, п — размерность вектора х, ||arg maxRn Tol ||2 —
Рис. 29. Проекция на детектор — свечение типа "лес". Модельное распределение (синяя линия) и решение техникой из [19] (красная линия)
Fig. 29. Projection onto the detector — luminosity of the "wood" type Model distribution (blue line) and solution by the technique from [19] (red line)
максимальное значение распознающего функционала, ||arg max^n Tol ||2, ||Ь||2 — нормы векторов аргумента экстремума распознающего функционала и правой части.
Оценка величин в (39) дает значение у/п (min^ ^ cond2 А} — 103 для разбиений расчетной области с п — 102. Величина min^ ^ cond2 А — 102 найдена методами диагональной оптимизации [21]. Интервализация матрицы А была проведена с использованием радиусов ее элементов, равных 1 мм: rad а^ = 10-3 Vi,j. Процедура интервализации определена в подразд. 0.2. Полученные оценки дают возможность судить об изменении величины множества решений при ненулевых значениях maxRn Tol:
/ . ,ч 3 ^ , II arg maxRn Tol IL , ,
ive(A; b) - 103 ■ max Tol ■ ^——R-^. (40)
||ащ тахКп То11|2
Ввиду того что величина -гт—п-~ I, выражение (40) можно использовать
П П2
в качестве оценки влияния величины maxgn Tol на размеры допускового множества.
2
Заключение
Реконструкция объемной светимости плазмы токамака — многостадийное исследование. Настоящая работа имеет дело с геометрическими аспектами проблемы и решением ряда модельных задач. В постановке предполагается, что плазма имеет азимутальную пространственную симметрию с известной границей в меридиональном сечении.
В работе представлен алгоритм нахождения проекций светимости объемов плазмы на плоскость матричного детектора. Исследован набор модельных прямых задач для построения алгоритмов их распознавания по данным детектора. Получены изображения предполагаемых макроскопических структур и микроскопических включений. Получены решения для набора тестовых задач, что показывает возможность успешной реконструкции для реальных данных.
Методической основой работы является применение методов интервального анализа для решения геометрических и алгебраических задач. Такой подход позволяет элегантно с минимальным количеством вычислительных затрат получать качественные и количественные результаты с учетом неопределенности входных данных. Исследование геометрических структур с интервальными входными данными служит основой для понимания характера проекций макроскопических структур светимости объема на плоскость детектора. Алгебраическая разрешимость исследуется в интервальной постановке. Результатом исследования является информация о наличии неопределенностей в геометрических данных и связанных с ними вычислениях за счет получения результатов о светимости плазмы путем решения задач линейного программирования.
Развитие исследований предполагает следующие направления. По мере получения решений, соответствующих принятым физическим представлениям, можно модифицировать способ разбиения, сгущая сетку в областях со значимыми величинами неодно-родностей и разрежая ее в областях с более равномерным распределением светимости. Интересна также трехмерная постановка задачи с решением недоопределенной системы уравнений. При этом решение двумерной задачи дает разумное первое приближение. Исследования в области интервальной регуляризации показывают [20], что таким образом достижимо получение содержательных результатов и для плохообусловленных задач. Также интересно расширить число методов оптимизации распознающего функционала.
Благодарности. Авторы выражают благодарность С.П. Шарому за общие консультации по методам интервального анализа, С.И. Жилину — за консультации и помощь
в расчетах, а также участникам российского вебинара по интервальному анализу за
конструктивные обсуждения математических аспектов исследования.
Список литературы
[1] Гусев В.К., Голант В.Е., Гусаков Е.З., Дьяченко В.В., Ирзак М.А., Минаев В.Б., Мухин Е.Е. и др. Сферический токамак Глобус-М. Журнал технической физики. 1999; 69(9):58-62.
[2] Днестровский Ю.Н., Лядина Е.С., Саврухин П.В. Пространственно-временная томографическая задача диагностики плазмы. № 5201/15. М.: ИАЭ; 1990: 46.
[3] Сладкомедова А.Д. Исследование радиационных потерь плазмы сферического токамака "Глобус-М": Дис. ... канд. физ.-мат. наук. СПб.: ФТИ им. А.Ф. Иоффе РАН; 2017:18.
[4] Уникальная научная установка (УНУ) "Сферический токамак Глобус-М". http://globus. rinno.ru/pages/status_unu-5.html
[5] Кабанихин С.И. Обратные и некорректные задачи. Учебник для студентов высших учеб. заведений. Новосибирск: Сиб. науч. изд-во; 2009: 457.
[6] Kearfott B., Nakao M., Neumaier A., Rump S., Shary S.P., Van Hentenryck P.
Standardized notation in interval analysis. Computational Technologies. 2010; 15(1):7-13.
[7] Шарый С.П. Конечномерный интервальный анализ. Новосибирск: Изд-во "XYZ"; 2019: 617. Адрес доступа: http://www.nsc.ru/interval/Library/InteBooks/SharyBook. pdf
[8] IEEE Std 1788-2015 — IEEE standard for interval arithmetic. Available at: https:// standards.ieee.org/standard/1788-2015.html
[9] Lao L.L., John H. St., Stambaugh R.D., Kellman A.G., Pfeiffer W. Reconstruction of current profile parameters and plasma shapes in tokamaks. Nuclear Fusion. 1985; 25(11): 1611-1622.
[10] Matsuoka R., Maruyama S. Eccentricity on an image caused by projection of a circle and a sphere. XXIII ISPRS Congress, 12-19 July 2016, Prague, Czech Republic. ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences. 2016; (III-5):19-26.
[11] Kessenich J., Sellers G., Shreiner D. The OpenGL® Programming Guide. 9th edition. Available at: http://www.opengl-redbook.com/
[12] Interval package. Available at: https://wiki.octave.org/Interval_package
[13] DeTurck D., Gluck H., Pomerleano D., Shea Vick D. The four vertex theorem and its converse. Notices of the American Mathematical Society. 2007; 54(2):192-207.
[14] Иванов А.А., Мартынов А.А., Медведев С.Ю. и др. Вычислительный код SPIDER.
Математическое моделирование равновесия и эволюции плазмы токамака. Вопросы атомной науки и техники. Сер.: Термоядерный синтез. 2014; 37(1):80-96.
[15] Савелов A.A. Плоские кривые. Систематика, свойства, применение: Справочное руководство. М.: Гос. изд-во физ.-матем. лит-ры; 1960: 296.
[16] Баженов А.Н., Затылкин П.А. Применение интервального подхода для спектрального анализа. Исследование разрешимости задачи с использованием распознающего функционала. Измерительная техника. 2019; (2):6-12.
[17] ШШарый С.П., ШШарая И.А. Распознавание разрешимости интервальных уравнений и его приложения к анализу данных. Вычислительные технологии. 2013; 18(3):80-109.
[18] Стецюк П.И. Субградиентные методы ralg5 и ralg4 для минимизации овражных выпуклых функций. Вычислительные технологии. 2017; 22(2):127—149.
[19] Zhilin S.I. Simple method for outlier detection in fitting experimental data under interval error. Chemometrics and Intelligent Laboratory Systems. 2007; 88(1):60-68.
[20] Ш!арый С.П. О мере вариабельности оценки параметров в статистике интервальных данных. Вычислительные технологии. 2019; 24(5):90-108.
[21] Сергеев Я.Д., Квасов Д.Е. Диагональные методы глобальной оптимизации. М.: Физ-матлит; 2008: 352.
Вычислительные технологии, 2020, том 25, № 1, с. 5-38. © ИВТ СО РАН, 2020 ISSN 1560-7534
Computational Technologies, 2020, vol. 25, no. 1, pp. 5-38. © ICT SB RAS, 2020 elSSN 2313-691X
MATHEMATICAL MODELLING
DOI:10.25743/ICT.2020.25.1.002 Small angle reconstruction of plasma luminosity for spherical tokamak
Bazhenov, Alexander N.1,2, Zatylkin, Pavel A.1,2
:A.F. Ioffe Physical-Technical Institute of the RAS, St. Petersburg, 194021, Russia 2SPbPU Peter the Great, St. Petersburg, 194021, Russia
Corresponding author: Bazhenov, Alexander N., e-mail: a_bazhenov@inbox.ru
Received May 28, 2019, revised November 11, 2019, accepted November 22, 2019
Abstract
The problems of reconstruction of plasma luminosity are important for physics and technology of power plants-tokamaks. The Globus-M research tokamak obtained a large amount of data using a matrix detector in pinhole camera geometry.
From the mathematical point of view, finding the luminosity for different regions of the plasma volume according to the matrix detector is an inverse problem related to the field of integral geometry. An essential feature of the particular task is the use of a single fixed camera with a small viewing angle. In this regard, application of methods of harmonic analysis of data is not enough.
The paper investigates the geometric aspects of the problem. In the general view, a three-dimensional object is projected onto a two-dimensional plane through a diaphragm. Under the assumption of azimuthal symmetry, there is a central projection of the luminosity of the body of rotation onto a flat matrix detector. The initial information for the calculation is the plasma boundary obtained from magnetic sensors. There is no reliable information about the internal structure of the plasma, so its division into regions of the equal luminosity is not unambiguous.
The paper presents an algorithm for finding the projections of the luminosity of plasma volumes on the plane of the matrix detector. A set of model direct problems for the construction of algorithms for their recognition according to the detector data was investigated. Images of supposed macroscopic structures and microscopic inclusions were obtained.
The methodological basis of the work is the use of interval analysis methods for solving geometric and algebraic problems. This approach allows obtaining qualitative and quantitative results that takes into account the uncertainty of the input data with the minimum amount of computational costs. Algebraic solvability is investigated in the interval formulation using response functionality.
Solutions for a set of test problems are obtained, which demonstrate the availability of successful reconstruction for real data. An important result of the study is an information about the presence of uncertainties in geometric data and related calculations by obtaining results about the luminosity of the plasma by solving linear programming problems.
Keywords: interval system of linear algebraic equations, tolerance set, measure of variability, recognizing functional, linear programming.
Citation: Bazhenov, A.N., Zatylkin, P.A. Small angle reconstruction of plasma luminosity for spherical tokamak. Computational Technologies. 2020; 25(1):5-38. (In Russ.)
Acknowledgements. The authors express their gratitude to S.P. Shary for general advice on interval analysis methods, to S.I. Zhilin for advice and assistance in calculations, and the participants of the Russian webinar on interval analysis for the constructive discussions of mathematical aspects of the study.
References
1. Gusev V.K., Golant V.E., Gusakov E.Z., D'yachenko V.V., Irzak M.A., Minaev V.B., Mukhin E.E. et al. Spherical tokamak Globus-M. Technical Physics. 1999; 69(9):58-62. (In Russ.)
2. Dnestrovskiy Yu.N., Lyadina E.S., Savrukhin P.V. Space-time tomographic problem of plasma diagnostics. No. 5201/15. Moscow: IAE; 1990: 46. (In Russ.)
3. Sladkomedova A.D. Issledovanie radiatsionnykh poter' plazmy sfericheskogo tokamaka "Globus-M" [Study of radiation losses for plasma of a spherical tokamak Globus-M]. Sankt-Peterburg: Fiz.-tekhn. in-t im. A.F. Ioffe RAN; 2017: 18. (In Russ.)
4. Unique scientific instrument (UNU) "Spherical tokamak Globus-M". Available at: http://globus. rinno.ru/pages/status_unu-5.html (In Russ.)
5. Kabanikhin S.I. Obratnye i nekorrektnye zadachi. Uchebnik dlya studentov vysshikh uchebnykh zavedeniy [Inverse and ill-posed Problems]. Novosibirsk: Sibirskoe nauchnoe izdatel'stvo; 2009: 457. (In Russ.)
6. Kearfott B., Nakao M., Neumaier A., Rump S., Shary S.P., Van Hentenryck P. Standardized notation in interval analysis. Computational Technologies. 2010; 15(1):7-13.
7. Shary S.P. Konechnomernyy interval'nyy analiz [Finite-dimensional interval analysis]. Novosibirsk: Izdatel'stvo "XYZ"; 2019: 617. Available at: http://www.nsc.ru/interval/Library/InteBooks/ SharyBook.pdf (In Russ.)
8. IEEE Std 1788-2015 — IEEE standard for interval arithmetic. Available at: https://standards. ieee.org/standard/1788-2015.html
9. Lao L.L., John H. St., Stambaugh R.D., Kellman A.G., Pfeiffer W. Reconstruction of current profile parameters and plasma shapes in tokamaks. Nuclear Fusion. 1985; 25(11):1611-1622.
10. Matsuoka R., Maruyama S. Eccentricity on an image caused by projection of a circle and a sphere. XXIII ISPRS Congress, 12-19 July 2016, Prague, Czech Republic. ISPRS Annals of the Photogram-metry, Remote Sensing and Spatial Information Sciences. 2016; (III-5):19-26.
11. Kessenich J., Sellers G., Shreiner D. The OpenGL® Programming Guide. 9th edition. Available at: http://www.opengl-redbook.com/
12. Interval package. Available at: https://wiki.octave.org/Interval_package
13. DeTurck D., Gluck H., Pomerleano D., Shea Vick D. The four vertex theorem and its converse. Notices of the American Mathematical Society. 2007; 54(2):192-207.
14. Ivanov A.A., Martynov A.A., Medvedev S.Yu. et al. The spider code. Mathematical modeling of tokamak plasma equilibrium and evolution. Problems of Atomic Science and Technology. Ser.: Thermonuclear Fusion. 2014; 37(1):80-96. (In Russ.)
15. Savelov A.A. Ploskie krivye. Sistematika, svoystva, primenenie: Spravochnoe rukovodstvo [Plane curves. Systematics, properties, application (reference guide)]. Moscow: Gosudarstvennoe Izdatel'stvo Fiziko-matematicheskoy Literatury; 1960: 296. (In Russ.)
16. Bazhenov A.N., Zatylkin P.A. Interval approach for spectral analysis. Use of recognizing functional to study of the solvability of a problem. Mesaurement Techniques. 2019; (2):6-12. (In Russ.)
17. Shary S.P., Sharaya I.A. Recognizing solvability of interval equations and its application to data analysis. Comput. Technologies. 2013; 18(3):80-109. (In Russ.)
38
A. H. BaxeeHOB, n. A. 3aTbij[KHH
18. Stetsyuk P.I. Subgradient methods ralg5 and ralg4 for minimization of ravine-like convex functions. Comput. Technologies. 2017; 22(2):127-149. (In Russ.)
19. Zhilin S.I. Simple method for outlier detection in fitting experimental data under interval error. Chemometrics and Intelligent Laboratory Systems. 2007; 88(1):60-68.
20. Shary S.P. On a variability measure for estimates of parameters in the statistics of interval data. Comput. Technologies. 2019; 24(5):90-108. (In Russ.)
21. Sergeev Ya.D., Kvasov D.E., Diagonal methods of global optimization Moscow: Fizmatlit; 2008: 352. (In Russ.)
