CC BY
52
1
ОБЗОРНАЯ СТАТЬЯ
УДК 535.3+517.926.4+519.642.7
ИНФРАКРАСНАЯ ТОМОГРАФИЯ ГОРЯЧЕГО ГАЗА: МАТЕМАТИЧЕСКАЯ МОДЕЛЬ АКТИВНО-ПАССИВНОЙ ДИАГНОСТИКИ
В.С. Сизиков
Сизиков Валерий Сергеевич - доктор технических наук, кандидат физико-математических наук, профессор кафедры измерительных технологий и компьютерной томографии СПбНИУ ИТМО и Института международных образовательных программ при СПбГПУ. Области научных интересов - томография, обработка изображений, спектроскопия, гравиметрия, гидроакустика, астрофизика, прикладная математика, интегральные уравнения, некорректные задачи и регулярные методы их решения. Автор более 130 научных работ, из них 9 книг (монографии и учебники). Член Американского математического общества, референт журнала «Mathematical Reviews», член диссертационного совета при НИУ ИТМО. Научный руководитель шести кандидатских и консультант одной докторской диссертации.
Сформулированы основные задачи инфракрасной томографии горячего газа на примере пламени горелки (burner). Описан вариант, когда используются два режима диагностики пламени: активный (ON), с включенным источником просвечивания, и пассивный (OFF) - без источника. Выведены два интегральных уравнения относительно коэффициента абсорбции к и функции Планка B среды (по которой можно рассчитать температурный профиль среды Tg). В случае осевой симметрии и параллельного сканирования пламени уравнения преобразованы в одномерные сингулярные интегральные уравнения типа Абеля относительно к и B. Для их численного решения использованы обобщенный метод квадратур, сглаживание данных с помощью сплайнов и метод регуляризации Тихонова. Разработан пакет программ на MATLAB7. С его помощью выполнена обработка результатов экспериментальной диагностики пламени горелки при некотором волновом числе v в некотором сечении пламени. Особенностью методики является то, что она не требует специального определения к путем прямого измерения или с помощью базы данных, например HITRAN/HITEMP. Обзор предназначен для студентов, аспирантов, преподавателей и научных сотрудников, специализирующихся по томографии, интегральным уравнениям, некорректным задачам, прикладной математике и программированию.
Ключевые слова: ИК томография, активный и пассивный режимы диагностики, коэффициент абсорбции, температурный профиль, интегральные уравнения, осевая симметрия, параллельное сканирование.
Введение
Рассматривается задача инфракрасной (ИК) томографии (IR tomography) [1-4] горячего газа, согласно которой определяются характеристики газовой среды по результатам ее внешней диагностики посредством ИК излучения, проходящего через среду.
Области применения ИК томографии и решаемые задачи:
- определение распределения абсорбции (поглощения), эмиссии (излучения) и температуры горячих газов - пламени лабораторной горелки [5], плазмы [6-10], пламени котлов, топок и парогенераторов [11], горячих газовых потоков [12], например, газа, вытекающего из сопла ракеты;
- ИК спектроскопическая томография температуры и концентрации частиц применительно к диагностике горения и обжига пылевидного угля или биомассы [12];
- точечная (не по сечениям) диагностика, или лазерная термометрия пламени горелки на основе когерентного антистоксова Раман-рассеяния (coherent anti-Stokes Raman scattering, CARS) [5; 13, С. 303, 615];
- определение параметров слоистой атмосферы - температуры, давления, затухания, абсорбции, эмиссии и др. [14-16] в рамках спутниковой метеорологии [17];
- тепловая (ИК) томография неразрушающего контроля композиционных материалов - анализ изменения поверхностной температуры материалов во времени с помощью нейронной сети из персептронов, тепловой контроль углепластика в авиакосмической промышленности и др. [18];
- в биомедицинской оптике - зондирование биотканей (кожи и др.) излучением ближнего ИК диапазона [19];
- ИК томография времени жизни и диффузионной длины носителей заряда в слитках полупроводникового кремния [20].
Видим, что область применения ИК томографии довольно широкая и разнообразная. Однако мы остановимся, главным образом, на ИК томографии горячих газов, в которых имеют место абсорбция и
эмиссия (примеры: пламя лабораторной горелки, котлов, топок и парогенераторов, плазма, газ, вытекающий из сопла ракеты). Это - одна из задач физики горения (combustion) [4, 5, 21, 22]. При этом будем также иметь в виду примыкающую задачу определения параметров атмосферы на основе ИК радиации, эмитированной атмосферой и перехватываемой сенсорами спутников, как задачу удаленного восприятия (remote sensing) профиля атмосферной температуры [17]. Задачи ИК томографии горячих газов и определения параметров атмосферы имеют ту общность, что они все используют дифференциальное уравнение переноса излучения (в разных вариантах).
Объяснение различий РКТ и ИК томографии. Задача ИК томографии имеет как сходства, так и различия с рентгеновской компьютерной томографией (РКТ) [8, 23-28]. Как известно [13, С. 227, 638], область длин волн X рентгеновских лучей составляет примерно 10-5-102 нм, а ИК излучение занимает широкую область от 740 до 2 • 106 нм, но мы будем рассматривать более узкую область ИК излучения -от 2000 до 15000 нм. При этом будем рассматривать область температур среды Tg от 1000 до 2000°С. На
рис. 1 приведена функция Планка B для длин волн X от 0 до 15000 нм и для трех значений Tg . Из рис. 1 видно, что в области длин волн рентгеновских лучей функция Планка B имеет очень малые значения, а значит, излучательная способность е практически отсутствует. В области же длин волн ИК томографии функция Планка имеет повышенные значения, поэтому излучательная способность среды существенна. Это порождает следующие различия РКТ и ИК томографии.
0,7
0,6
о
Й «
cq tf
СП g
g
к й н С
к О
0,5
0,4
0,3
0,2
0,1
1
1
2
3
J У ---
1000 2000 3000
4000 5000 6000 7000 8000 9000 10000 11000 12000 13000 14000 15000 Длина волны X, нм
Рис. 1. Функция Планка газа B(k,Tg). 1 - Tg = 2000°С; 2 - Tg = 1500°C; 3 - Tg = 1000°C
В РКТ имеет место лишь поглощение (абсорбция) сканирующих лучей, проходящих через вещество; искомой функцией является одна функция - коэффициент абсорбции к; используется лишь активный режим (с включенным источником просвечивания); необходимо решать лишь одно уравнение - интегральное уравнение (ИУ) Радона относительно к.
В ИК томографии имеет место не только абсорбция, но и эмиссия газа, поэтому искомыми функциями являются две - коэффициент абсорбции к и коэффициент излучения е, связанный с функцией Планка среды B (которая связана с температурным профилем среды Tg); в идеале необходимы (как, например, в [4]) два режима измерений: активный (в [4] он назван режимом ON - с источником или с зеркалом [1]) и пассивный (в [4] он назван режимом OFF - без источника и без зеркала [1], когда источником излучения является сама среда), и нужно решать два интегральных уравнения (ИУ) относительно к и B. Однако часто [2, 3] рассматривают лишь режим ON и получают одно (дифференциальное или интегральное) уравнение относительно двух искомых функций, к и B. Чтобы преодолеть этот дефицит данных (одно уравнение и две искомых функции), нередко полагают [2], что коэффициент к измерен (некоторым образом) или используют базу данных HITRAN/HITEMP [29, 30] и подобные для определения к.
В данной работе для ИК томографии горячего газа используются два режима измерений - ON и OFF. Будут получены два ИУ относительно двух функций к и B. Заметим, что в работах [1, 4] и др. также сформулирована задача эмиссионно-абсорбционной томографии на основе измерений с внешним источником или зеркалом и без него (with and without an external source or mirror). В работе [1] для выполнения
0
диагностики осесимметричной плазмы использовано зеркало (которое предложил ранее Реагсе), расположенное с противоположной стороны от детектора (приемника) (рис. 2). Зеркало позволяет обойтись без источника, заменяя его. В дальнейшем вместо схемы с зеркалом стали использовать систему источник-приемник (как и в данной работе), при этом уравнения, полученные в [1] для схемы с зеркалом, оказались в общем и целом справедливыми для системы источник-приемник.
Ых)
(детектор)
т
У\(х)
| В(Т0) (источник) или ///,'// (зеркало)
Рис. 2. Прохождение луча через z-сечение пламени
В [4] отмечено, что разность ON и OFF измерений позволяет определить коэффициент абсорбции к, а затем определить функцию Планка B и температуру Tg . Приведены выражения для интенсивностей в
виде контурных интегралов и реконструктивные алгоритмы Рамачандрана-Лакшминарайяна и Шеппа-Логана. Однако в [4] методика изложена слишком кратко; при этом не использованы такие эффективные приемы, как сглаживание экспериментальных данных, регуляризация и т.д. (эти вопросы рассмотрены в других публикациях - см. дальше).
Заметим, что в обратной задаче диагностики плазмы [2-4, 6-10] используется как активный режим (просвечивание плазмы), так и пассивный режим, когда источником излучения является сама плазма, однако задача не доведена до двух интегральных уравнений. Также не доведена до двух уравнений методика диагностики слоистой атмосферы [14-17].
В данной работе эти вопросы рассмотрены объединенно и по-новому.
Основные соотношения
Уравнение переноса ИК излучения. Рассмотрим горелку (burner) [5], создающую ламинарное ровное пламя (рис. 3). Такая горелка часто используется для тестирования различных методов и алгоритмов [2-5].
Пусть через некоторое z-сечение пламени параллельно оси y (рис. 3) передается монохроматическое излучение абсолютно черного тела (а.ч.т.) 2hc2 v3
B(T>) =-
(1)
е л^/ад _1:
где к - постоянная Планка, с - скорость света в вакууме, V = 1/X - волновое число, X - длина волны, кв - постоянная Больцмана, Т0 - температура а.ч.т.
Излучение, посылаемое источником (а.ч.т), испытывает абсорбцию с коэффициентом абсорбции к = к (х, у). Кроме того, среда излучает с коэффициентом излучения е = е( х, у) = е(Т? (х, у)), где
Tg = Тв (х, у) - распределение температуры по г-сечению среды. Имеет место дифференциальное уравнение переноса излучения с учетом абсорбции и эмиссии, но без рассеяния [1-4, 8, 14, 31]:
Ш (х, у)
dy
■ = -к(x, y) I(x, y) + е(x, у),
(2)
где I(х, у) - искомая интенсивность излучения, причем к(х, у(х)) = 0 при у(х) < у1 (х) и у(х) > у2 (х). где у1(х) и у2(х) - границы среды (рис. 2).
Абсолютно
черное тело, Т0*900°С
В(То)
X ^
Горелка
Рис. 3. Горелка
Полагаем, что в каждой точке (х, у) имеет место локальное термодинамическое равновесие [31]. Это позволяет использовать закон Кирхгофа [31]: отношение излучательной и поглощательной способ-
з
ностей равно — = В(Т ), где В(Т ) - функция Планка для среды (ср. (1)):
к К К
ВТ) = В(Т (х, у)) =-
2Ис2 V3
,ИсЧкВТ£ (Х>У)
- 1
(3)
или е(х, у) = к(х, у) В(ТК (х,у)). В результате уравнение переноса излучения (2) примет вид (ср. [2, 9, 14])
(4)
^^ = к(х, у) [В(^ (х, у)) -1(х, у)]: ёу
где х играет роль параметра, т.е. уравнение (4) справедливо при каждом значении х (см. рис. 4, где представлена схема параллельного сканирования некоторого г-сечения пламени при некотором ракурсе ф).
Рис. 4. Параллельное прохождение лучей через г-сечение пламени
Аналитическое решение уравнения переноса излучения. Решение уравнения (4) для каждой точки (х, у) каждого луча имеет следующий вид (ср. [1, 2, 4, 9, 31; 32, С. 308]):
I (x, y) = B(T))exp
( у Л у ( у Л
_ | к(х, у') йу' + | к(х, у') ВТ (х, у'))ехр _{ к(х, у") йу" V у1 (х) / у1 (х) V у'
Для дальнейшей обработки более важна интенсивность на детекторе:
( у2( х) Л у2( х) ( у2( х) А
dy'.
Ir (x) = B(T)) exp
- J к(x, y) dy + J к(x, y) B(Tg (x, y))exp - J к(x, y') dy' dy. (5)
ч y1 (x) / Я(x) V у /
Соотношение (5) следует рассматривать как ИУ относительно двух функций к и B. Чтобы уравнений было также два, помимо активного режима ON с источником (1), дающего функцию IR, используем режим пассивной диагностики пламени OFF, когда нет источника и его роль играет само пламя. Обозначим измеренную детекторами интенсивность в режиме OFF через Ig (x). Выражение для Ig (x) имеет
вид (ср. [1, 4] и (5)):
y2(x) ( y2(x) ^
Ig (x) = J к(x, y) B(Tg (x, y)) exp - J к(x, y') dy' dy. (6)
y (x) V У /
Теперь задачу можно сформулировать следующим образом: по двум измеренным функциям IR и Ig нужно определить две функции к и B путем решения системы двух ИУ (5) и (6). Однако искомые функции к и B - это функции двух переменных x и y, а исходные (измеренные) функции IR и Ig - функции лишь одной переменной x. Исходя из этого, необходимо, как в РКТ, выполнить измерения под разными ракурсными углами ф и получить IR (x, ф) и Ig (x, ф) или рассмотреть, например, случай, когда изолинии к и B являются окружностями (случай осевой симметрии, см. дальше). В этом случае искомые функции будут функциями одной переменной: к = к (r), B = B(Tg (r)), где r - расстояние от оси пламени.
Данную задачу можно упростить. Обозначим через IT разность двух измеренных функций IR и Ig , т.е. IT (x) = IR (x) -Ig (x). Получим уравнение, описывающее процесс абсорбции без эмиссии в активном режиме:
( y2 (x)
it (x) = B(T)) exp
J к (x, y) dy
Я(x)
(7)
В результате получаем не систему двух ИУ, а два самостоятельных ИУ, решаемых последовательно: сначала решается уравнение (7) относительно к(х, у), а затем уравнение (6) относительно В(Т^, (х, у))
при уже найденной к(х, у). После вычисления В(Т^, (х, у)) = В(х, у) можно рассчитать температурный
профиль (см. (3)):
кв
Te (x, y) =
in 12hc!v_+1
(8)
В х, у)
При этом функции 1Т и должны быть получены, как сказано выше, под различными ракурсами Ф, т.е. получены 1Т (х, ф) и Iв (х, ф) или же рассмотрена схема осевой симметрии (в принципе, при одном ракурсе), когда 1т = 1т (х), ^ = 18 (х), к = к (г), В = В(Т^, (г)).
Случай осевой симметрии и параллельного сканирования
Рассмотрим случай, когда в каждом г-сечении пламени изолинии постоянных к и Т^, (а также е)
являются окружностями, т.е. имеет место круговая (осевая) симметрия. Кроме того, полагаем, что выполняется параллельное сканирование (рис. 5). Вариант осевой симметрии и параллельного сканирования рассмотрен в ряде работ ([1, 4, 8, 9, 21, 22, 33] и др.). В данной работе дается дальнейшее развитие этого варианта.
Определение коэффициента абсорбции к. Рассмотрим ИУ (7) относительно к(х, у). Учтем осевую симметрию пламени. Введем вместо у новую переменную г (расстояние от г-оси симметрии), и уравнение (7) примет вид сингулярного интегрального уравнения (СИУ) Абеля:
я
г
2J
л/Т2-
7 к(r)dr = q(x), 0 < x < R,
(9)
где
q( x) = - ln
it (x) B(T0)'
(10)
Рассмотрим вопрос о численном решении СИУ (9). Отметим следующие основные алгоритмы численного решения уравнения (9) (ср. [34]).
Ш
.(детектор)
ЩТ0) (источник)
Рис. 5. Осевая симметрия и параллельное сканирование пламени Алгоритм 1 (сдвиг сеток узлов) [35-38]. Вводятся дискретные сетки узлов по переменным г и х:
rf = jh, xi = r + Д, j, i = 0,1,
= R , где h = R/n - шаг дискретизации, Д - сдвиг между сетками,
равный Д = к/ 2 [35, С. 8; 37, С. 273], или Д е [0, к/ 2] [36, С. 27], или Де (0, к/ 2) [38]. Введение сдвига А позволяет устранить особенности при использовании численных методов, например, метода квадратур. Алгоритм 2 (варианты метода квадратур). В работе [35, С. 8] при построении метода дискретных
1 г 1
вихрей (одного из вариантов метода квадратур) интеграл в уравнении — I
2л •'-1
1 г i у(x)dx
= f (xo), xo e (-1,1)
(связанном с уравнением типа (9)) расписывается по формуле левых прямоугольников и используются сетки по x и x0 со сдвигом Д = h/2 . В результате получается система линейных алгебраических уравнений (СЛАУ) относительно значений уj = у(xj).
В работах [21, 33] изложен метод под названием «onion-peeling» («луковая кожура»), в котором область r e [0, R] аппроксимируется кольцами постоянных значений к между rj -Дг/2 и rj + Дг/2 для
каждого rj. Сетки узлов по r и x полагаются одинаковыми (без сдвига). Особенность метода состоит в
' rj +Дг /2 rdr
при r, - Дг /2 > xi берется аналитически и имеет конечное значение.
том,
что интеграл j
ТГ2^
Далее используется формула средних прямоугольников, и получается СЛАУ с верхней прямоугольной матрицей. Похожий метод изложен в работе [34].
Алгоритм 3 (аппроксимация решения) [6, С. 154]. В работах [36, С. 15, 27; 39, С. 149, 159] искомое решение к (г) и (или) д( х) аппроксимируется алгебраическим или тригонометрическим многочленом или полиномиальным сплайном, коэффициенты которого определяются из условия минимума невязки, что приводит к проекционному методу (Галеркина, коллокации и др.) и к решению СЛАУ.
В перечисленных алгоритмах имеет место саморегуляризация, причем в случае использования сдвига сеток величина сдвига А играет роль параметра регуляризации: чем ближе А к к/2 , тем устойчивее решение к (г), но меньше разрешающая способность численного метода, и чем ближе А к 0, тем менее устойчиво решение, но выше разрешение. Кроме того, во всех алгоритмах получаются СЛАУ с преобладающей (но не бесконечной) диагональю матрицы.
Отметим еще ряд алгоритмов. Уравнение (9), как известно, имеет аналитическое решение [6, 8, 21, 22, 33, 34]
1 R
к (r) = - -j
тг J
q'(x)
- lyjx2 - r2
dx, 0 < r < R .
(11)
Однако решение (11) содержит производную q'( x) от экспериментальной, а значит, зашумленной функции q(x), а задача дифференцирования является некорректной [40, С. 18]. Кроме того, интеграл в (11) является несобственным (сингулярным). Тем не менее, для вычисления решения согласно (11) предложен ряд следующих алгоритмов.
Алгоритм 4 (интерполяция и метод квадратур). В работе [33] для вычисления производной q'( x)
г r dx
использована интерполяция по трем, а также по двум соседним точкам. Интеграл в (11) I . вы-
^ Vx2 - r2
числяется аналитически. Похожий алгоритм изложен в работе [34] (использована обобщенная формула левых прямоугольников).
Алгоритм 5 (аппроксимация q (x)) [7, 41, 42]. Функция q( x) аппроксимируется линейной комбинацией сглаживающих полиномов (или сплайнов) на всем интервале x е [0,R]. Вычисляется q'(x) путем дифференцирования полинома (или сплайна). Решение к (r) согласно (11) вычисляется путем суммирования интеграла в (11) по отрезкам, что выполняется аналитически [7, С. 188-189].
Алгоритм 6 (без использования q'(x)). В работе [43] выражение (11) с помощью интегрирования по частям преобразовано в выражение, не содержащее производной q' (x) [6-8, 21]:
к =-1^q(R^ + Jx[q(x^dx^, 0 < r < R.
* I Vr2 - r2 r V(x2 - r2)3 J
Данный алгоритм реализован в работе [7, С. 217-220] с использованием кубического сплайна [44, P. 273] для q(x).
Алгоритм 7 (с использованием регуляризации). ИУ Абеля (9) вследствие сингулярности имеет, с одной стороны, сложности в реализации алгоритмов, а с другой стороны, обладает саморегуляризацией, в результате чего задача его решения является умеренно некорректной [21]. В ряде работ ([21, 22, 45] и др.) для повышения устойчивости алгоритмов использован метод регуляризации Тихонова [40, 46-48].
В данной работе рассматривается решение уравнения (9) относительно искомой функции к(r), а также численное вычисление к (r) согласно (11). Новым является комплексный подход - предварительное сглаживание экспериментальных функций I R и Ig , использование обобщенного метода квадратур (в модификации) для решения уравнения (9) и вычисления (11) и метода регуляризации Тихонова.
Обобщенный метод квадратур решения сингулярного интегрального уравнения
Рассмотрим СИУ (9). В работе [33] предложен, а в работе [21] реализован численный метод «onion-peeling», в котором использованы равномерные совпадающие сетки узлов по r и x и квадратурная формула средних прямоугольников для вычисления интеграла в (9). В работе [34] применены также равномерные совпадающие сетки по r и x и формула левых прямоугольников.
В данной работе используются неравномерные сетки и квадратурная формула левых прямоугольников, что порождает более общий и удобный алгоритм.
Итак, рассмотрим решение уравнения (9) на сетках узлов, вообще говоря, неравномерных и несовпадающих:
0 < x1 < r1 < x2 < r2 < x3 < r3 <... < xt < rt <... < xn < rn < R (12)
или rt = xt + A,, A( е [0,(xt+1 -xt )/2]. Здесь R = rmax - граничное значение, при котором к(R) = 0 . На каждом промежутке [r.., r.+1), j е [1, n -1] полагаем приближенно к(r) = к(rj) = = const. Получим
0+1
J~nr=T к (r) dr = (V rj+1 - x2 - j*2) к
r. \r x
7' (13)
; е [1, п _ 1], х < г) < г!+! < я, - обобщенную квадратурную формулу левых прямоугольников для сингулярности г/-%/г2 _ х2 , причем множители ^г}+1 _ х2 _ ^г2 _ х21 являются квадратурными коэффициентами этой формулы. Специфика
формулы состоит в том, что сингулярный интеграл | У+' Г = йг в^гчисляется аналитически точно и
г ыг2_х2
без особенности. Если же вычислять его численно по обычной квадратурной формуле левых прямоугольников, то при х = г возникает деление на ноль.
i £ [1, П - 1].
Интеграл в (9) есть сумма интегралов (13):
1~п=к (г) *=Е (V у - х2 ^л/г;гтх;г )■
В результате СИУ (9) преобразуется в
п-1 -- .-
Е А кУ = 2, '' е [1, П - 1] , где А = VгУ+1 - х, -\1 гУ- х,2
(14)
т.е. СЛАУ с верхней треугольной матрицей (п -1) х (п -1) относительно ку, у е [1, п -1]. Рекуррентное решение СЛАУ (14) имеет следующий вид:
к - gn-,/2 к - к + ( гп - П ^
кп-1~ ' кп~ кп-2 +
рп-1, п-1
П-1
2 - S Pjkj
' к-1 - kn-2),
к. _-
j-i+1
i - п - 2,П - 3,...,1,
(15)
где kj = к (г ), qt = q( xj). В решение (15) добавлено вычисление кп с помощью линейной экстраполяции [49, С. 213, ф. (8.80)]. В работе [33] также приведены формулы типа (13), но для случая равномерных и совпадающих сеток по г и x и с использованием формулы средних прямоугольников, при этом формулы типа (15) не приведены. Формулы (12)-(14) являются более общими и удобными, чем в работах [21, 33], а формула (15) - явной. Метод согласно (12)-(15) назван в работе [34] обобщенным методом квадратур решения СИУ (9).
Приведенный метод согласно (12)-(15) может быть применен также к численному вычислению решения к (г ) согласно (11) (ср. [33, 34]), однако в данной работе мы ограничимся численным вычислением к (г ) путем решения СИУ (9) согласно (12), (15).
Определение функции Планка В. Теперь рассмотрим соотношение (6) как ИУ относительно B(x, y) - B(T (x, y)) при уже найденной к(x, y). Соотношение (6) можно записать в виде
Ух( x )
i
У]( x)
J к(x,y)B(x,y)dy - Ig(x), -R < x < R.
где
( y2 ( x )
k(x, y) - к(x, y) exp
- J к(x, y') dy'
V y
(16)
(17)
- ядро интегрального уравнения, а В( х, у) = В(Т^ (х, у)) - искомая функция.
Учтем осевую симметрию и параллельное сканирование. При этом, чтобы учесть особенности записи интегралов в (16) и (17) при использовании полярной координаты г, запишем интеграл (16) в виде суммы двух интегралов (рис. 6):
У2 ( х ) 0 У2 ( х )
| к(х,у)В(х,у) dy = | к(х,у) В(х,у) dy + | к(х,у) В(х,у) dy . (18)
у, (х) у, (х) 0
y,(,x)Zy<i) 0 <, у < у,(.\)
а б
Рис. 6. Осевая симметрия и параллельное сканирование пламени:
y £ [Ух(x),0] (а) и y £ [0, y2 (x)] (б)
Для первого интеграла в (18) (т.е. при у < 0) имеем у = _4г2 _ х2 , ёу = _
тёг, а для вто-
рого интеграла (т.е. при у > 0 ) у = -у/г2 _ х2 , ёу = . г =ёг . Сумма этих интегралов равна
^- И!+! «г, .г
г' к (г')
г к (г)
,ехр ■!
! _! .г : в(г) а>
х х N гх2 ]
= ехР Нда Л- '[■!
г к (г)
тг7^
ехр|!
г г' к (г')
'2 _ х2
Лг'
г г' к (г')
г
ехр1 _!-/::
Лг'
В(г) Лг.
где
В результате уравнение (16) с учетом (19) запишется в виде СИУ типа ИУ Абеля:
л
!К(х,г) В(г)Лг = 0(х), 0 < х < Я , гк(г)
К(х,г) = (е"(хг5 + е_"(хг5
л/г2 _ х2 1
б( х) = / (х) е"( хя\
г г' к (г')
л Г г '
1( х, г ) = 1-;= Лг'
2 _ х2
ёг', 0 < х < Я, х < г < Я .
(19)
(20)
(21) (22) (23)
Видим, что, во-первых, само уравнение (20) является сингулярным, так как его ядро К (х, г) (21) содержит сингулярный множитель, и, во-вторых, интеграл (23) является сингулярным. Используем, как и при решении СИУ (9), обобщенный метод квадратур. В результате сингулярный интеграл (23) будет приближенно равен (ср. (13))
[0, у = ,,
■- ■-\ (24)
"(х,, г]) =
,е[1,п_1], у>!
Ё (у/г^ _jr7)kf, у > ,,
на сетках узлов (12), а СИУ (20) с учетом (22)-(24) преобразуется в СЛАУ
Ё рВу = й , , е [1, п _ 1],
где
Р = (Я^х к] .(е"г> + е_"г>):
(25)
(26)
е, = ^ (х,) - е"(я> . (27)
СЛАУ (25) имеет верхнюю треугольную матрицу (п _ 1) х (п _ 1). Ее рекуррентное решение (ср.
(15)):
В = , в = В 2 +
п_1 £> > п п_2
Рп _1,п _1
•(Ви_1 _ Ви_2),
В,. =-
у=,+1
р.
, = п _2,п _3,...,1,
(28)
где В1 = В(г), = е(х1), а Р у и евыражаются формулами (26) и (27). После получения В(г) можно найти температурный профиль (см. (8)):
Нс\/ кв
Те (г) =
Ы^СУ. + ! В(г)
(29)
Использование регуляризации
Задача решения интегральных уравнений I рода (7) и (6), а также их частных случаев (9) и (20) является, строго говоря, некорректной, несмотря на то, что интегральные уравнения (9) и (20) являются сингулярными, а сингулярные ИУ обладают, как известно [35, С. 194; 37, С. 495], свойством естествен-
г
ной регуляризации, или саморегуляризации, обусловленным тем, что при их численном решении квадратурами матрица получающейся СЛАУ имеет преобладающую диагональ. Кроме того, интегральные уравнения Вольтерра I рода, к которым относятся ИУ (9) и (20), занимают промежуточное положение между ИУ Фредгольма I рода и Вольтерра II рода [46, С. 110-114]. Другими словами, задача решения ИУ Вольтерра I рода может быть корректной (well-posed) или некорректной (ill-posed) в зависимости от того, какими методами и в каких пространствах она решается.
В данной работе выясняется, в какой степени метод обобщенных квадратур (без регуляризации, но со сглаживанием результатов измерений, см. дальше) делает задачу решения СИУ (9) и (20) устойчивой и насколько повышает ее устойчивость метод регуляризации Тихонова (ср. [10, 21, 22, 40, 45-48]). В задаче определения коэффициента абсорбции к запишем СЛАУ (14) в виде Ак = f, (30)
где
Aij = j Pj J > ^ f = qj2, i, j e [1,n -1]. (31)
I 0, j < i,
Согласно методу регуляризации Тихонова, вместо (30) с учетом (31) решается СЛАУ (а/ + ATA) ка = ATf, (32)
где а > 0 - параметр регуляризации, / - единичная матрица. Часто [10, 16, 21, 22, 45-47] регуляризован-ная СЛАУ записывается в виде
{allT + ATA) = ATf (33)
или
{clallt + ATA) = ATf , (34)
где матрица L (n -1) x n - дискретный лапласиан. Записи (33) и (34) имеют более высокий порядок регуляризации, чем (32).
Важным является вопрос о выборе параметра регуляризации а. Для выбора а обычно используются принцип невязки, критерий L-кривой, метод перекрестной значимости и др. [7, 10, 22, 46-48]. В данной работе мы остановимся на принципе невязки. Согласно нему, а выбирается согласно равенству
Aa- f II = 8 , (35)
где 5 - верхняя оценка шума по норме: 8 = ||8f||2. Здесь f = q/2, поэтому (см. (10)) 8f = ~~ = ~. Подробнее:
I JUI 8/„ 1 I
, (36)
8=2 IS
Ч ^
Its
где |8/T| = /Tm -IT s|, /Tm - измеренные (measured) значения /T , а /Ts - значения /T , полученные аппроксимирующим (сглаживающим) сплайном (spline) (см. рис. 9), которые предлагается использовать в качестве «точных» значений /т , поскольку точные значения /т неизвестны в случае реальных (не модельных) измерений.
Невязка - f || может быть записана как
,1/2
IK-f|L [(^ка )i - f, ]2 ^ =|S
Л V2
SlAjj(ка)j -f J I . (37)
Принцип невязки (35)-(37) эффективен, когда f зашумлена, а погрешность 5 оценена надежно. Если решается модельный пример с заданным точным коэффициентом абсорбции к, то а нужно выбирать по минимуму относительной погрешности решения (ср. [45])
||ка- ц
8а ' (38)
где к - точный (заданный) коэффициент абсорбции.
Обозначим через а(1 параметр регуляризации, выбранный по принципу невязки (discrepancy), а через а^ - по минимуму погрешности 8а (38).
В задаче определения функции Планка B запишем СЛАУ (25) в виде
AB = Q , (39)
где
* = 't
j ^ i, j <i,
i, j е [1, n -1].
Согласно методу регуляризации Тихонова, вместо (39) с учетом (40) решаем СЛАУ
(а/ + АТА) Ва = Ате .
(40)
(41)
Численные иллюстрации
Разработан пакет программ на MATLAB7 для решения задачи определения коэффициента абсорбции к и температурного профиля Tg по результатам активно-пассивной диагностики пламени горелки в
предположении осевой симметрии пламени и параллельного сканирования его. Приведем численные иллюстрации решения этой задачи.
На рис. 7 представлены экспериментальные (measured) функции IR m (x) и I (x), а также
IT,m (x) = IR,m (x) -1g,m (x) на следующей сетке узлов:
x = 0; 0,4; 0,8; 1,2; 1,6; 2; 2,4; 2,8; 3,2; 3,6; 3,8 см (42)
при v = 2271,496 см-1, T0 = 894,4°С, n = 12, HAB = 12 мм (height above burner).
Измерения выполнены в Laboratory Rise DTU, Denmark в рамках совместного проекта [12].
0,25 Г-
1Яш(х) (ON), Вт-м-1
в—
0,2
0,15
0,1
0,05
0,5
1,5
2,5
3,5
х, см
Рис. 7. Измеренные функции /Ят(х), /^т(х) и /Тт(х)
На рис. 8 показан коэффициент абсорбции кт(г) (пунктирная кривая), найденный в результате численного решения СИУ Абеля (9) обобщенным методом квадратур согласно (15) на сетке узлов (42) при г = х1 на основе измеренной (шеа8шеф функции /Т т (х).
Видим, что решение кт(г) имеет значительные (нереальные) флуктуации. Это обусловлено грубостью шага дискретизации по х (другими словами, малостью п), а также наличием погрешностей измерения функций /Ят (х) и /^,т (х) (и, как следствие, /Т,т (х)).
Рис. 8 показывает также резко пониженное значение кт(0) (на оси симметрии пламени). Понижение значения к (0) наблюдалось и в других работах, например, [4, 21], но в меньшей степени.
Мы добавили еще ограничение на решение кт (г): определение (корректировка) такого значения кт (0), при котором к'т (0) = 0 (равенство нулю производной при г = 0 следует из физических соображений). Численно такое значение кт (0) можно получить на основе полинома Лагранжа 2-й степени, проходящего через точки кт (0), кт (г2) и кт (г3), и формулы (2.65) работы [46, С. 121] при /'(а) = а = 0 :
к =
r3 к2 r2 к3 (r2 + r3)(r3 - r2)
(43)
0
0
1
2
3
4
где к1 = кш(гг), к2 =кш(г2), к3 =кш(г3). В частном случае, когда г2 -г1 = г3 -г2, имеем к1 = (4к2 -к3)/3. На рис. 8 точками отмечена часть кривой кш (г) вблизи г1 = 0 с найденным (скорректированным) согласно (43) значением к1 = кш(0). Видим, что скорректированное значение к1 значительно отличается от нескорректированного и является более правдоподобным.
0,16
0,14 0,12 0,1 0,08 0,06 0,04 0,02
0 0,5 1 1,5 2 2,5 3 3,5 4
г, см
Рис. 8. Коэффициент абсорбции кш(г) (без регуляризации) и кша (г) (с регуляризацией)
На рис. 8 приведены также значения коэффициента абсорбции кша (г) (непрерывная кривая), найденные в результате решения СЛАУ (30) методом регуляризации Тихонова согласно (33). При этом параметр регуляризации а выбран по принципу невязки согласно (35), причем оценка погрешности 5 рассчитана по формуле (36). Получилось: 5 = 0,037 , а = 10-0 09 (см. далее рис. 11). Отметим, что 5 оценивается надежно, если в качестве «точных» значений 1Т использовать 1Т 8 - значения 1Т , полученные сглаживающим сплайном (см. рис. 9). 0,25
0,2
0,15
0,1
0,05
0
Рис. 9. Измеренные функции 1Кш(х), I (х) и 1Тш(х) ( о); узлы х = 0(0,2)3,8 кубических сглаживающих сплайнов (• ); непрерывные кривые - сплайны
Рис. 8 показывает, что метод регуляризации Тихонова сам заметно сгладил флуктуации в к (г), хотя при решении СЛАУ (33) использовались несглаженные 1Т , а именно, 1Т ш . Точками отмечена часть кривой кша (г) вблизи оси симметрии с корректированным согласно (43) значением кша (0). Заключения о
полученном решении kma (r) (а также km (r)) будут сделаны дальше (после использования сглаживания экспериментальных данных).
Чтобы уменьшить шаг дискретизации по x, а также сгладить погрешности измерений функций IR m (x) и Ig m (x) (и, как следствие, IT m (x)), был использован кубический сглаживающий сплайн. Вопросу о сглаживании данных посвящено много разработок. Использованы аппроксимирующие (сглаживающие) сплайны [7, 8, 10, 49], приближение линейной комбинацией полиномов [7, 41, 42], метод регрессии [44, 50, 51] и др. Например, в работе [51] использована процедура сглаживания экспериментальных зашумленных данных путем так называемой локально взвешенной регрессии (LOESS), в которой каждая точка сглаженных данных формируется посредством весового среднего некоторого числа соседних точек. Как отмечено в [22], эта процедура не устраняет прямо некорректность томографической задачи, но понижает неустойчивость решения.
В данной работе нами был использован кубический сглаживающий (аппроксимационный) сплайн [7, С. 42-86, 188-193; 10, С. 164; 44, Р. 273; 49, С. 223]. На рис. 9 представлена аппроксимация функций ir m (x), Ig m (x) и ITm (x) кубическими сглаживающими сплайнами с помощью m-функции csaps.m (со сглаживающим параметром p = 0.97, что соответствует умеренному сглаживанию). При этом было введено дополнительное условие для всех трех функций: I'(0) = 0 .
0,16
0,14 0,12 0,1 0,08 0,06 0,04 0,02
ka(r), 1/с] 1 1
k(r), 1/см
| 1 1 р
0 0,5 1 1,5 2 2,5 3 3,5 4
г, см
Рис. 10. Коэффициент абсорбции к (г), полученный путем решения СИУ (9) обобщенным методом квадратур согласно (15), и ка(г), полученный путем решения СЛАУ (33) методом регуляризации Тихонова
0,06 0,05
0,04
0,03 0,02 0,01
0
-3,5 -3 -2,5 -2 -1,5 -1 -0,5 0 0,5 1
1ёа
Рис. 11. Невязка ||Ака - /|| без сглаживания (пунктирная линия) и со сглаживанием (непрерывная линия) измеренной функции 1Т т( х)
0,25
0,2
0,15
0,1
0,5
-<---_ Q(x)
0,5
1
2,5
1,5 2
х, см
Рис. 12. Правая часть Q(x), Вт/м
3,5
0,5
0,4
0,3
0,2
0,1
0
B(r) .........................П
B<x(r)
0,5
1
1,5 2 r, см
2,5
3
3,5
4
Рис. 13. Функция Планка B(r) (без регуляризации) и Bx (r) (с регуляризацией), Вт/м
2500
2000
1500
1000
500
ТН(Г)
Ч \
Tgx(r) \ Ч ч \ \ Nw \
\ \ v х \ \
\ X \ \ ч ч \ \
0
0,5 1 1,5 2 2,5 r, см
3
3,5 4
Рис. 14. Температурные профили Tg (r) (без регуляризации), Tga (r) (с регуляризацией)
и TH(r) (по методике CARS), К
0
3
4
Со сплайна были сняты сглаженные значения IT (x) в 20 узлах x = 0(0,2)3,8 и по ним снова решено ИУ (9) обобщенным методом квадратур согласно (15) (см. решение к (r ) на рис. 10), а также решена СЛАУ (30) методом регуляризации Тихонова согласно (33). Параметр регуляризации а выбран способом невязки (5 = 0,003, а = 10"2'1) (решение ка(r) - на рис. 10, а выполнение способа невязки - на рис. 11). При этом использованы ограничения: к'(0) = к'а (0) = 0 . Из рис. 10 видно следующее:
1. сглаживание сплайном измеренных функций IR m(x), I (x) и IT m(x) привело к (умеренному)
сглаживанию решений к (r ) и ка (r ) ;
2. решения к (r ) (без регуляризации) и ка (r ) (с регуляризацией) получились практически одинаковыми.
Это говорит о том, что решение СИУ (9) можно получить, используя сглаживание измеренных функций и обобщенный метод квадратур без регуляризации, и при этом решение будет близким к регу-ляризованному решению.
Следующий шаг - определение функции Планка В и температурного профиля Tg. На рис. 12
представлена правая часть Q(x) уравнения (20), рассчитанная согласно (22), (23).
Путем решения СИУ (20) обобщенным методом квадратур (согласно (28)) была найдена функция Планка B(r) для среды (пламени горелки), представленная на рис. 13. Была найдена также функция
Планка Ва (r) методом регуляризации Тихонова путем решения СЛАУ (41) (см. рис. 13). В заключение был рассчитан температурный профиль среды (осесимметричного пламени горелки) Tg (r ) по B(r ) (без регуляризации) и Tga (r ) по Ва (r ) (с регуляризацией) согласно (29) - см. рис. 14. Видим, что результаты
получились весьма близкими без регуляризации (но со сглаживанием и обобщенными квадратурами) и с регуляризацией, что подтверждает то, что задача решения сингулярных интегральных уравнений является умеренно некорректной и обладает свойством саморегуляризации.
На рис. 14 для сравнения приведен один из температурных профилей TH(r ) (Hartung's profile), полученных методом CARS [5] также для газовой горелки, но с несколько иными условиями. Видим заметные различия, но качественно картины похожие.
Заключение
По-новому изложена известная ИК методика определения коэффициента абсорбции к и функции Планка среды В (по которой можно определить температурный профиль Tg ). Активно-пассивная диагностика позволяет получить две экспериментальные функции и два новых интегральных уравнения относительно к и В. Рассмотрен случай осевой симметрии (и параллельного сканирования), для которого справедливы сингулярные интегральные уравнения (СИУ) Абеля. Сформулирован (в новой редакции) обобщенный метод квадратур их решения относительно к и B. Выполнена обработка экспериментальных данных в двух вариантах - путем решения СИУ обобщенным методом квадратур и методом регуляризации Тихонова, причем без сглаживания и со сглаживанием экспериментальных данных. Эта обработка подтвердила, что СИУ обладает саморегуляризацией, задача его решения является умеренно некорректной, а использование обобщенного метода квадратур с предварительным сплайн-сглаживанием экспериментальных данных позволяет получить практически такие же результаты, как и методом регуляризации Тихонова.
В дальнейшем будут рассмотрены варианты веерного сканирования и изолиний в виде эллипсов. Автор благодарен А. Фатееву и В. Евсееву за предоставленные экспериментальные данные (рис. 7) и за полезные обсуждения.
Работа выполнена при поддержке РФФИ (грант № 13-08-00442) и DTU, Denmark (Project No. 010246).
Литература
1. Porter R.W. Numerical solution for local émission coefficients in axisymmetric self-absorbed sources // SIAM Review. - 1964. - V. 6. - № 3. - P. 228-242.
2. Tourin R.H., Krakow B. Applicability of infrared emission and absorption spectra to determination of hot gas temperature profiles // Applied Optics. - 1965. - V. 4 - № 2. - P. 237-242.
3. Krakow B. Spectroscopic temperature profile measurements in inhomogeneous hot gases // Applied Optics. -1966. - V. 5. - № 2. - P. 201-209.
4. Hall R.J. and Bonczyk P.A. Sooting flame thermometry using emission/absorption tomography // Applied Optics. - 1990. - V. 29. - № 31. - P. 4590-4598.
5. Härtung G., Hult J., Kaminski C.F. A flat flame burner for the calibration of laser thermometry techniques // Measur. Sci. Technol. - 2006. - V. 17. - P. 2485-2493.
6. Преображенский Н.Г., Пикалов В.В. Неустойчивые задачи диагностики плазмы. - Новосибирск: Наука, 1982. - 238 с.
7. Воскобойников Ю.Е., Преображенский Н.Г., Седельников А.И. Математическая обработка эксперимента в молекулярной газодинамике. - Новосибирск: Наука, 1984. - 240 с.
8. Пикалов В.В., Преображенский Н.Г. Реконструктивная томография в газодинамике и физике плазмы. - Новосибирск: Наука, 1987. - 239 с.
9. Пикалов В.В., Мельникова Т.С. Томография плазмы (Низкотемпературная плазма, Т. 13). - Новосибирск: Наука. Сиб. изд. фирма РАН, 1995. - 229 с.
10. Старков В.Н. Конструктивные методы вычислительной физики в задачах интерпретации. - Киев: Наук. думка, 2002. - 264 с.
11. Kästner W. et al. Application of evolutionary algorithms to the optimization of the flame position in coal-fired utility steam generators / E. Hüllermeier et al. (Eds.) // Proc. 13th Intern. Conf. IPMU 2010. - Part I. -Berlin: Springer, 2010. - V. 80. - P. 722-730.
12. Evseev V., Fateev A., Sizikov V., Clausen S., Nielsen K.L. On the development of methods and equipment for 2D-tomography in combustion // Report on Annual meeting of Danish Physical Society, 21-22 June 2011. - 32 p.
13. Физический энциклопедический словарь / Гл. ред. А.М. Прохоров. - М.: Сов. энциклопедия, 1984. -944 с.
14. Goody R.M., Yung Y.L. Atmospheric Radiation. Theoretical Basis. - 2nd ed. - New York-Oxford: Oxford University Press, 1989. - 536 p.
15. Huang B., Smith W.L., Huang H.-L., Menzel W.P. A hybrid iterative method for ATOVS temperature profile retrieval // Techn. Proc. 9th Intern. TOVS Study Conf., Igls, Austria, 20-26 Feb. 1997. - P. 177-187.
16. Doicu A., Trautmann T., Schreier F. Numerical Regularization for Atmospheric Inverse Problems. - Berlin: Springer, 2010. - 431 p.
17. Menzel W.P. Applications with Meteorological Satellites. Techn. Document WMO/TD № 1078. - Univ. Wisconsin, 2001. - 242 p.
18. Вавилов В.П., Нестерук Д.А., Ширяев В.В., Иванов А.И., Swiderski W. Тепловая (инфракрасная) томография: терминология, основные процедуры и применение для неразрушающего контроля композиционных материалов // Дефектоскопия. - 2010. - № 3. - С. 3-15.
19. Зимняков Д.А., Тучин В.В. Оптическая томография тканей // Квантовая электроника. - 2002. - Т. 32. -№ 10. - С. 849-867.
20. Ахметов В. Д., Фадеев Н.В. Инфракрасная томография времени жизни и диффузионной длины носителей заряда в слитках полупроводникового кремния // Физика и техника полупроводников. - 2001. -Т. 35. - Вып. 1. - С. 40-47.
21. Daun K.J., Thomson K.A., Liu F., Smallwood G.J. Deconvolution of axisymmetric flame properties using Tikhonov regularization // Applied Optics. - 2006. - V. 45. - № 19. - P. 4638-4646.
22. Äkesson E.O., Daun K.J. Parameter selection methods for axisymmetric flame tomography through Tikhonov regularization // Applied Optics. - 2008. - V. 47. - № 3. - P. 407-416.
23. Наттерер Ф. Математические аспекты компьютерной томографии. - М.: Мир, 1990. - 288 с.
24. Суинделл Б., Уэбб С. Рентгеновская трансмиссионная томография // Физика визуализации изображений в медицине. - М.: Мир, 1991. - Т. 1. - С. 138-173.
25. Гуров И.П., Сизиков В.С., Щекотин Д.С. Методы восстановления изображений в рентгеновской томографии // Научно-технический вестник СПбГУ ИТМО. - 2003. - № 5 (11). - С. 97-104.
26. Марусина М.Я., Казначеева А.О. Современное состояние и перспективы развития томографии // Научно-технический вестник СПбГУ ИТМО. - 2007. - № 8 (42). - С. 3-13.
27. Сизиков В.С. Обратные прикладные задачи и MatLab. - СПб: Лань, 2011. - 256 с.
28. Симонов Е.Н. Физика визуализации изображений в рентгеновской компьютерной томографии. - Челябинск: Изд-во НИУ ЮУрГУ, 2013. - 550 с.
29. Fleck T., Jäger H., Obernberger I. Experimental verification of gas spectra calculated for high temperatures using the HITRAN/HITEMP database // J. Phys. D: Applied Physics. - 2002. - V. 35. - № 23. - P. 31383144.
30. Rothman L.S. et al. The HITRAN molecular spectroscopic database and HAWKS (HITRAN Atmospheric Workstation): 1996 edition // J. Quant. Spectrosc. Radiat. Transfer. - 1998. - V. 60. - № 5. - P. 665-710.
31. Radiative transfer [Электронный ресурс]. - Режим доступа: http://en.wikipedia.org/wiki/Radiative_transfer, свободный. Яз. англ. (дата обращения 10.08.2013).
32. Бронштейн И.Н., Семендяев К.А. Справочник по математике для инженеров и учащихся втузов. - 13-е изд. - М.: Наука, 1986. - 544 с.
33. Dasch C.J. One-dimensional tomography: a comparison of Abel, onion-peeling, and filtered backprojection methods // Applied Optics. - 1992. - V. 31. - № 8. - P. 1146-1152.
34. Сизиков В.С., Смирнов А.В., Федоров Б.А. Численное решение сингулярного интегрального уравнения Абеля обобщенным методом квадратур // Изв. вузов. Математика. - 2004. - № 8 (507). - С. 62-70.
35. Белоцерковский С.М., Лифанов И.К. Численные методы в сингулярных интегральных уравнениях. -М.: Наука, 1985. - 256 с.
36. Габдулхаев Б.Г. Прямые методы решения сингулярных интегральных уравнений первого рода. - Казань: Изд-во Казанс. ун-та, 1994. - 288 с.
37. Лифанов И.К. Метод сингулярных интегральных уравнений и численный эксперимент. - М.: ТОО «Янус», 1995. - 520 с.
38. Бойков И.В. Приближенные методы решения сингулярных интегральных уравнений. - Пенза: Изд-во Пенз. гос. ун-та, 2004. - 316 с.
39. Иванов В.В. Теория приближенных методов и ее применение к численному решению сингулярных интегральных уравнений. - Киев: Наук. думка, 1968. - 287 с.
40. Тихонов А.Н., Арсенин В.Я. Методы решения некорректных задач. - 3-е изд. - М.: Наука, 1986. -288 с.
41. Minerbo G.N., Levy M.E. Inversion on Abel's integral equation by means of orthogonal polynomials // SIAM J. Num. Anal. - 1969. - V. 9. - № 4. - P. 598-616.
42. Косарев Е.Л. О численном решении интегрального уравнения Абеля // Журн. вычисл. матем. и матем. физики. - 1973. - Т. 13. - № 6. - С. 1591-1596.
43. Deutsch M., Beniaminy I. Derivative-free inversion of Abel's integral equation // Appl. Phys. Lett. - 1982. -V. 41. - № 1. - P. 27-28.
44. Martinez W.L., Martinez A.R., Solka J.L. Exploratory Data Analysis with MATLAB. - 2nd ed. - Boca Raton: CRC Press, 2010. - 495 p.
45. Daun K.J. Infrared species limited data tomography through Tikhonov reconstruction // J. Quant. Spectrosc. Radiat. Transfer. - 2010. - V. 111. - № 1. - P. 105-115.
46. Верлань А.Ф., Сизиков В.С. Интегральные уравнения: методы, алгоритмы, программы. - Киев: Наук. думка, 1986. - 544 с.
47. Engl H.W., Hanke M., Neubauer A. Regularization of Inverse Problems. - Dordrecht: Kluwer, 1996. -328 p.
48. Hansen P.C. Discrete Inverse Problems: Insight and Algorithms. - Philadelphia: SIAM, 2010. - 213 p.
49. Сизиков В.С. Математические методы обработки результатов измерений. - СПб: Политехника, 2001.
- 240 с.
50. Корн Г., Корн Т. Справочник по математике для научных работников и инженеров. - М.: Наука, 1968.
- 720 с.
51. Cleveland W.S., Devlin S.J. Locally weighted regression: an approach to regression analysis by local fitting // J. Amer. Stat. Assoc. - 1988. - V. 83. - № 403. - P. 596-610.
Сизиков Валерий Сергеевич - Россия, Санкт-Петербург, Санкт-Петербургский национальный исследо-
вательский университет информационных технологий, механики и оптики, профессор; Институт международных образовательных программ при СПб ГПУ, профессор; доктор технических наук, профессор, [email protected]
