ВОССТАНОВЛЕНИЕ НЕПРЕРЫВНЫХ СПЕКТРОВ АДАПТИВНЫМ СПОСОБОМ ВЫЧИСЛИТЕЛЬНЫХ ЭКСПЕРИМЕНТОВ С РЕГУЛЯРИЗАЦИЕЙ
ВОССТАНОВЛЕНИЕ НЕПРЕРЫВНЫХ СПЕКТРОВ АДАПТИВНЫМ СПОСОБОМ...
УДК 535.338.1+519.642.3+519.6
ВОССТАНОВЛЕНИЕ НЕПРЕРЫВНЫХ СПЕКТРОВ АДАПТИВНЫМ СПОСОБОМ ВЫЧИСЛИТЕЛЬНЫХ ЭКСПЕРИМЕНТОВ С РЕГУЛЯРИЗАЦИЕЙ
В.С. Сизиков, А.В. Кривых
Рассмотрена обратная задача спектроскопии – восстановление непрерывных спектров путем математической обработки измеренных спектров, искаженных аппаратной функцией спектрометра и помехами. Задача сводится к решению интегрального уравнения Фредгольма I рода. Задача его решения некорректна, поэтому для получения устойчивого решения используется метод регуляризации Тихонова. При этом применен адаптивный способ вычислительных экспериментов, согласно которому, наряду с исходным спектром P, обрабатывается модельный спектр Q с задаваемым истинным спектром z и моделируемым измеренным спектром u с учетом дополнительной (априорной) информации об истинном спек-
тре P. Это позволяет выбрать параметр регуляризации . Предложенная методика может быть использована для по-
вышения разрешающей способности спектрометра. Приведены численные иллюстрации. Ключевые слова: непрерывный спектр, обратная задача спектроскопии, интегральное уравнение, метод регуляризации Тихонова, способ вычислительных экспериментов, повышение разрешающей способности спектрометра.
Введение Измеренный спектрометром (например, интерферометром Фабри–Перо) спектр u() (где – длина волны) обычно отличается от истинного спектра z() [1–8]. Это проявляется, во-первых, в большей сглаженности спектра u() по сравнению с z() , а именно, в спектре u() не разрешены близкие линии, сглажена тонкая структура спектральной линии, что является результатом воздействия аппаратной функции спектрального прибора [1–9]. Во-вторых, это проявляется в зашумленности спектра u() , а именно, слабые линии «тонут» в шуме, что является результатом погрешностей измерений [1–3], а также воздействия среды, через которую проходит излучение [10]. Дадим следующее определение аппаратной функции (АФ) [3, 6–8] (ср. [9, С. 32, 704]): аппаратной функцией K(, ) спектрометра называется его реакция (в виде измеренной интенсивности) на дискретную линию единичной интенсивности и длины волны при настройке спектрометра на длину волны . Форма аппаратной функции (ширина и т.д.) может заметно меняться с изменением длины волны настройки [min , max ] , где [min , max ] – диапазон длин волн изучаемой части спектра. Обычно с увеличением АФ становится шире, что характерно для широкополосной спектрометрии, например, изучения спектра звезды во всем видимом диапазоне. Если же АФ практически не изменяется при изменении , то АФ является разностной (инвариантной): K(,) K( ) , что имеет место, например, при изучении тонкой структуры отдельной линии [3, 6, 8], когда диапазон [min , max ] мал. На рис. 1 в качестве примера приведен смоделированный непрерывный измеренный спектр u() , сглаженный аппаратной функцией спектрометра K(, ) , а также зашумленный (и дискретизированный) измеренный спектр u() u() u (где u – шум) и АФ спектрометра, причем, поскольку в дан-
22 Научно-технический вестник информационных технологий, механики и оптики,
2013, № 3 (85)
В.С. Сизиков, А.В. Кривых
ном примере K(, ) – функция неразностная, то приведено два ее «сечения» (подробности примера см.
дальше). В принципе похожий вид может иметь непрерывный узкополосный спектр [6, С. 200], например, сверхтонкая структура отдельной линии, обусловленная магнитными или электрическими полями (эффект Зеемана или Штарка), а также тепловым уширением (эффект Доплера) [10], однако в этом случае диапазон [min , max ] мал, а АФ – разностная: K(, ) K( ) .
6
Ииннттееннссииввнноосстть,ь,у.у.е.е.
5 u()
4 u~()
3
2
1 K(620, )
K(485, ) 0
450 460 470 480 490 500 510 520 530 540 550 560 570 580 590 600 610 620 630 640 650
min
, , нм
max
Рис. 1. u() – измеренный без шума спектр; u() – измеренный зашумленный и дискретизированный
спектр; K ( , ) – АФ при некоторой длине волны настройки ; [min , max ] – широкий диапазон длин волн
Как будет видно далее, в примере на рис. 1 в измеренном спектре u() (тем более, в зашумленном
спектре u() ) не разрешены близкие линии и не выявлены слабые, причем этот эффект тем сильнее, чем
шире АФ K(, ) (а также чем выше уровень шумов), другими словами, чем меньше разрешающая спо-
собность спектрометра [1, 9]. В данной работе ставится известная обратная задача спектроскопии – задача восстановления ис-
тинного спектра z() по измеренному спектру u() и аппаратной функции K(, ) [1–8, 11–15]. Дан-
ная задача описывается интегральным уравнением (см. дальше), задача решения которого некорректна, поэтому его обычно решают методом регуляризации Тихонова. При этом важным является вопрос о вы-
боре параметра регуляризации . В данной работе предлагается новый адаптивный способ (вычисли-
тельных экспериментов) для выбора параметра .
Математическая формулировка обратной задачи спектроскопии
Рассмотрим случай непрерывного спектра, обычно характерного для веществ с повышенной плотностью (расплавленный жидкий металл, плазма и т.д.). Измеренная интенсивность u() при настройке
спектрометра на длину волны равна сумме (интегралу) по всем истинным интенсивностям z() с ве-
совой функцией K:
b
u() z() K (, ) d , a
где a min , b max , откуда, варьируя значение (т.е. выполняя сканирование по спектру) и учитывая
зашумленность спектра u() , получим:
b
K (, ) z() d u() , c d ,
a
где [c, d] – пределы изменения (обычно более широкие, чем [a,b] ).
(1)
В соотношении (1) известны (измерены или заданы) u() , K(, ) , a, b, c, d, а z() является ис-
комым истинным спектром. Соотношение (1) есть интегральное уравнение Фредгольма I рода, причем
Научно-технический вестник информационных технологий, механики и оптики, 2013, № 3 (85)
23
ВОССТАНОВЛЕНИЕ НЕПРЕРЫВНЫХ СПЕКТРОВ АДАПТИВНЫМ СПОСОБОМ...
K(, ) является ядром уравнения, u() – правой частью, а z() – искомой функцией. Если K(,) K( ) , то
K ( ) z() d u() , 0 .
(2)
0
Соотношение (2) есть интегральное уравнение Фредгольма I рода типа свертки на полуоси. Реше-
ние уравнения (1) или (2) дает возможность, в принципе, восстановить истинный спектр z() . Однако
задача решения уравнений (1) и (2) является некорректной (существенно неустойчивой) [2–4, 6, 8, 16]: если решать уравнение (1), например, методом квадратур, а уравнение (2) – методом преобразования Фурье (инверсной фильтрации), то в качестве решения получим так называемую «пилу» [3, 6, 8] – крайне неустойчивое решение. По этой причине для устойчивого решения этих уравнений необходимо применение таких методов, как регуляризация Тихонова [2–4, 6–8, 11–16], параметрическая фильтрация Винера [3, 6, 8, 16], итеративная регуляризация Фридмана [6, 8, 16] и др.
При обработке спектра в широком диапазоне длин волн следует учитывать изменение формы АФ
K(, ) с изменением длины волны настройки . При обработке же спектра в узкой полосе следует ис-
пользовать уравнение Фредгольма I рода с разностным ядром (ср. (2)):
b
K ( ) z() d u() , c d .
(3)
a
Задача решения уравнений (1)–(3) связана с задачей редукции к идеальному спектральному при-
бору [1–4, 9, 17] – с одним из вариантов редукционной проблемы Рэлея [3, 6, 8, 13]. Успешное решение
задачи редукции позволит путем математической обработки результатов измерений повысить разре-
шающую способность спектрального прибора. В настоящей статье воспользуемся методом регуляриза-
ции Тихонова. Что касается других устойчивых методов (фильтрации Винера, итеративных методов и
др.), то они изложены в различных публикациях ([3, 6, 8, 16] и др.) и также могут быть применены для
устойчивого восстановления спектров.
Краткая формулировка метода регуляризации Тихонова
Запишем уравнение (1) в виде
b
A z K (, ) z() d u() , c d ,
(4)
a
где A – оператор, соответствующий ядру K. Метод регуляризации Тихонова сводится к решению инте-
грального уравнения Фредгольма II рода
b
z (t) B(t, ) z () d U (t) , a t b , a
где 0 – параметр регуляризации, а новое ядро и новая правая часть равны
(5)
dd
B(t, ) B(,t) K (,t) K (, ) d , U (t) K (,t)u() d . cc
В таком варианте уравнение (5) обычно решается методом квадратур [3, 4, 6, 8, 16]. Если же рассматривать уравнение (2) или (3), то его решение методом регуляризации Тихонова будет включать преобразование Фурье и -регуляризацию (подробности см. в [3, 4, 6, 8, 11–16, 18–21]).
При этом важным является вопрос о выборе параметра регуляризации и об учете дополнитель-
ной (априорной) информации относительно искомого спектра z() . Существует ряд способов выбора
параметра регуляризации : способ невязки, обобщенный принцип невязки, метод перекрестной значимости, локальный регуляризующий алгоритм, способ подбора и др. [3, 4, 6, 8, 13, 16, 18–21].
Способ вычислительных экспериментов
В данной работе получает дальнейшее развитие способ вычислительных экспериментов для выбора параметра регуляризации (другие его названия – способ псевдообратного оператора, способ эталонных, или модельных примеров, способ моделирования) [3, 6, 8, 16, 22, 23]. Данный способ учитывает дополнительную (априорную) информацию об искомом спектре (оценку количества спектральных линий, их параметров и т.д.) и поэтому является интерактивным и адаптивным способом.
Кратко изложим способ вычислительных экспериментов. Рассмотрим операторное уравнение I рода: A z u (ср. (4)). Полагаем, что вместо точных u и K из-
вестны u и K такие, что || u u || , || A A || , где и – верхние оценки погрешностей по норме
правой части u и ядра K. При использовании метода регуляризации Тихонова решается операторное
24 Научно-технический вестник информационных технологий, механики и оптики,
2013, № 3 (85)
В.С. Сизиков, А.В. Кривых
уравнение z AT A z AT u (ср. (5)), где AT – транспонированный оператор. Обозначим z z z –
погрешность регуляризованного решения z , а z – точное решение (нормальное псевдорешение [16, 20, 21]). В работах [16, 22, 23] получена следующая оценка относительной погрешности регуляризованного решения по норме:
||
z || z ||
||
()
,
(6)
где
()
|| A
2
||
p p 1
.
(7)
Здесь отн отн , причем отн || u || и отн || A || – относительные погрешности исходных
данных; p || A ||2 , A – псевдообратный оператор: Au z [20, С. 184]. Функция () является верх-
ней огибающей для истинной относительной погрешности
отн
()
||
||
z z ||
||
.
(8)
В работах [16, 22] показано, что функция () имеет (единственный) минимум при условии
p (|| A || )2 27 16 1, 69 или || A || || A || 3 3 4 1, 30 . Согласно соотношениям (6)–(8), оценка
относительной погрешности || z || || z || регуляризованного решения z зависит от A и (точнее, от
произведения || A || ). По этой причине, если решается несколько задач (другими словами, обрабатыва-
ется несколько спектров) с одинаковыми A и , то для них оценки погрешности
отн ()
||
z || z ||
||
|| A~ ||
2
p p 1
будут одинаковыми.
Отсюда следует, что при решении некоторого исходного примера P (т.е. при обработке исходного
спектра uP ) с неизвестным решением (спектром) zP можно использовать результаты решения другого,
модельного, примера Q с известным (заданным) точным решением (спектром) zQ , причем с такими же
A и , что и в примере P. При этом при решении примера Q можно рассчитать функцию
отн () Q || z Q || || zQ || и по ней найти опт Q – оптимальное значение , при котором
отн () Q
min
.
Это
значение
опт Q
может быть использовано при решении исходного примера (спек-
тра) P. При этом необходимо также определить p || A ||2 . Оценка p может быть получена путем подбо-
ра такого значения p, при котором огибающая кривая () касается набора кривых отн ()Q (см. рис. 3).
Добавим, что для повышения эффективности изложенного способа модельный пример Q (или несколько примеров) должен содержать дополнительную информацию об исходном примере (спектре) P, а именно, оценку количества спектральных линий (максимумов) в искомом спектре zP , соотношений их
интенсивностей и значений их длин волн. Данную оценку должен делать опытный спектроскопист. Ис-
пользование такой информации в модельном примере Q позволит более удачно выбрать параметр .
Данный способ следует считать адаптивным и интерактивным способом.
Численная иллюстрация
В рамках системы программирования MATLAB 7 был разработан пакет программ для восстанов-
ления истинных непрерывных спектров z() путем численного решения интегрального уравнения
Фредгольма I рода методом регуляризации Тихонова с использованием способа вычислительных экспериментов.
Сначала был рассмотрен первый пример (рис. 1) – оригинал P, у которого известен зашумленный
измеренный спектр u() на равномерной сетке узлов min , min h,, max , где min 450 нм ; max 650 нм ; h const 1 нм – шаг дискретизации; n (max min ) h 200 – число шагов дискретизации по . Известна также аппаратная функция – дифракционная АФ Рэлея (ср. [1, 4, 5]) вида
K (, )
1 ()
sinc2
()
1 ()
sin ( ) ()2
( ) ()
,
(9)
Научно-технический вестник информационных технологий, механики и оптики, 2013, № 3 (85)
25
ВОССТАНОВЛЕНИЕ НЕПРЕРЫВНЫХ СПЕКТРОВ АДАПТИВНЫМ СПОСОБОМ...
где () – полуширина АФ по уровню 0, равная приблизительно ширине АФ по уровню 0,5, которую мы
положили равной () 8 min нм . Спектр полагается широкополосным (от фиолетового до красного), поэтому ширина АФ непостоянна, а именно, (min ) 8 нм , а (max ) 11, 55 нм , т.е. (max ) (min ) 1, 44 . При этом истинный спектр z() в примере P неизвестен.
Из рис. 1 видно, что измеренный спектр u() имеет довольно сложную структуру, а именно, со-
держит шесть явных флуктуаций, две из которых (при 525 нм и 620 нм ), скорее всего, состоят каждая из двух линий, но они не разрешились из-за того, что АФ имеет немалую ширину и, тем самым, ограничивает разрешающую способность спектрометра. Кроме того, есть намек на то, что при 507 нм и 543 нм имеются еще две слабые линии. Таким образом, все указывает на то, что на самом деле в спектре имеются не менее восьми спектральных линий. В связи с этим в качестве второго (модельного или эталонного) примера Q был составлен близкий к оригиналу P пример, истинный спектр zQ () которого состоит из 9 спектральных линий в виде гауссиан:
zQ () 2, 0 exp{[( 486) /10]2} 0, 4exp{[( 512) / 5]2}
8,5exp{[( 522) / 2]2} 9, 2exp{[( 530) / 2]2}
0,5exp{[( 542) / 5]2} 8, 2exp{[( 566) / 6]2}
2,5exp{[( 592) / 4]2} 4,5exp{[( 614) / 7]2}
3, 0 exp{[( 626) / 5]2}.
Измеренный спектр uQ () в примере Q был рассчитан согласно выражению
b
uQ () K (, ) zQ () d , c d , a
численно. При этом a = 460, b = 640, c = 450, d = 650 нм. Погрешности измеренной uP () были оценены примерно в 1%, что соответствует среднеквадра-
тическому отклонению СКО ≈ 0,02. По этой причине к значениям uQ () были добавлены случайные
нормальные погрешности с СКО от 0,01 до 0,025, что соответствует отн 0,5 1,25% , поскольку значе-
ние отн в исходном примере P известно неточно. АФ спектрометра в примере Q была взята в виде (9),
причем (поскольку АФ известна также неточно) () было взято равным () (8 ) min , где 0 0,3 , что соответствует отн 0 3% .
Далее модельный пример Q был решен методом квадратур с регуляризацией Тихонова с помощью
разработанной m-функции Tikh.m [6, С. 207] для ряда значений параметра регуляризации , и была
построена зависимость относительной погрешности регуляризованного решения z () по отношению к точному решению z() (см. (8)):
отн ()
z () z() z()
.
На рис. 2 представлены зависимости отн () для ряда погрешностей отн и отн . На рис. 2 пред-
ставлена также огибающая () (см. (7)), при построении которой было положено 102 и
|| A |||| A || || u ||L2 || z ||L2 0,82 . Для ряда значений p от 100 до 270 были рассчитаны кривые ()
(рис. 2). Было выбрано то значение p, при котором одна из кривых касается набора кривых отн () , а
именно, p = 100. Этому соответствует значение параметра регуляризации 103 . Из рис. 2 видно, что,
несмотря на разброс кривых отн () и () , значения p и, как следствие, определяются уверенно.
При значении 103 , выбранном с помощью решения модельного примера Q как вспомогательного восстановлен спектр в исходном примере P (рис. 3). Как видно из рис. 3, в примере P разрешились близкие линии и восстановились слабые линии, правда, на краях спектра проявился эффект Гиббса, однако в слабой форме (на уровне погрешностей метода). Аналогичные результаты получены для других, весьма различных, непрерывных спектров [3, 6–8, 14, 22, 23], т.е. изложенная в работе методика вычислительных экспериментов может быть использована для широкого класса спектров (с близкими линиями, со слабыми линиями, узкими и широкими линиями и т.д.).
26 Научно-технический вестник информационных технологий, механики и оптики,
2013, № 3 (85)
В.С. Сизиков, А.В. Кривых
1
0,9
0,8
0,7
0,6
отн ()
0,5
0,4
p=270 p=200
() p=100
0,3
0,2
0,1
0–9 –8 –7 –6 –5 –4 –3 –2 –1 0 1 lg
Рис. 2. Зависимости отн () для ряда погрешностей отн и отн и огибающие ()
2
9 11
8
7
интеИннстиевннсоисвтнь,осу.теь., у.е.
6 5 22
4
3
2 33
1
0
460 470 480 490 500 510 520 530 540 550 560 570 580 590 600 610 620 630 640 , , нм
Рис. 3. 1 – истинный спектр zP () ; 2 – измеренный спектр u() ; 3 – восстановленный спектр z P ()
Заключение
Практическое использование изложенной методики позволит повысить разрешающую способность спектрометра. Спектральный прибор может быть соединен с компьютером или со спецпроцессором с заложенным в него математическим и программным обеспечением, реализующим методы и численные алгоритмы решения обратной задачи спектроскопии. В результате такого комплексирования (соединения прибора с компьютером) можно разрешить близкие и выделить слабые линии спектров излучения (или поглощения), а именно, в физике – спектров газов, жидкостей, металлов, плазмы; в астрофизике – спектров звезд, планет, галактик, туманностей, комет; в металлургии – спектров расплавленных металлов в домнах; в геофизике – спектров залежей руд, минералов, нефти, газа и т.д.
Работа выполнена при поддержке РФФИ (грант № 13-08-00442).
Научно-технический вестник информационных технологий, механики и оптики, 2013, № 3 (85)
27
ВОССТАНОВЛЕНИЕ НЕПРЕРЫВНЫХ СПЕКТРОВ АДАПТИВНЫМ СПОСОБОМ...
Литература
1. Раутиан С.Г. Реальные спектральные приборы // Успехи физических наук. – 1958. – Т. 66. – Вып. 3. – С. 475–517.
2. Кочиков И.В., Курамшина Г.М., Пентин Ю.А., Ягола А.Г. Обратные задачи колебательной спектроскопии. – М.: Изд-во МГУ, 1993. – 204 с.
3. Сизиков В.С. Математические методы обработки результатов измерений. – СПб: Политехника, 2001. – 240 с.
4. Старков В.Н. Конструктивные методы вычислительной физики в задачах интерпретации. – Киев: Наукова думка, 2002. – 264 с.
5. Fleckl T., Jäger H., Obernberger I. Experimental verification of gas spectra calculated for high temperatures using the HITRAN/HITEMP database // J. Phys. D: Appl. Phys. – 2002. – V. 35. – P. 3138–3144.
6. Сизиков В.С. Обратные прикладные задачи и MatLab. – СПб: Лань, 2011. – 256 с. 7. Сизиков В.С., Кривых А.В. Применение способа эталонных примеров при решении обратной задачи
спектроскопии методом регуляризации // Изв. вузов. Приборостроение. – 2011. – Т. 54. – № 9. – С. 44– 51. 8. Сизиков В.С. Интегральные уравнения и MatLab в задачах томографии, иконики и спектроскопии. – Saarbrücken: LAP, 2011. – 252 c. 9. Физический энциклопедический словарь / Гл. ред. А.М. Прохоров. – М.: Сов. энциклопедия, 1984. – 944 с. 10. Ландсберг Г.С. Оптика: Учебное пособие для вузов. – 6-е изд. – М.: ФИЗМАТЛИТ, 2003. – 848 с. 11. Преображенский Н.Г., Седельников А.И. Оптимизация спектроскопических измерений на основе методов регуляризации // Журнал прикладной спектроскопии. – 1981. – Т. 35. – Вып. 4. – С. 592–599. 12. Брагинская Т.Г., Клюбин В.В. Решение обратной задачи спектроскопии оптического смещения методом регуляризации Тихонова. Препринт № 855. – Л.: ЛИЯФ, 1983. – 60 с. 13. Глазов М.В., Болохова Т.А. Решение редукционной проблемы Рэлея с использованием различных модификаций метода регуляризации // Оптика и спектроскопия. – 1989. – Т. 67. – Вып. 3. – С. 533– 537. 14. Кривых А.В., Сизиков В.С. Комплексированное восстановление непрерывных спектров с использованием псевдообратной матрицы // XLI Неделя науки СПбГПУ: материалы научно-практической конференции с международным участием. Ч. XIII. – СПб: Изд-во Политехн. ун-та, 2012. – С. 240–242. 15. Кривых А.В., Сизиков В.С. Обработка дискретных спектров с помощью алгоритма интегральной аппроксимации // Научно-технический вестник СПбГУ ИТМО. – 2011. – № 5 (75). – С. 14–18. 16. Верлань А.Ф., Сизиков В.С. Интегральные уравнения: методы, алгоритмы, программы. – Киев: Наукова думка, 1986. – 544 с. 17. Краулиня Э.К., Лиепа С.Я., Пикалов В.В., Скудра А.Я. К проблеме исследования атомной сенсибилизированной флуоресценции по контурам спектральных линий // Некорректные обратные задачи атомной физики / Сб. статей под ред. Н.Г. Преображенского. – Новосибирск: ИТПМ, 1976. – 133 с. 18. Воскобойников Ю.Е., Преображенский Н.Г., Седельников А.И. Математическая обработка эксперимента в молекулярной газодинамике. – Новосибирск: Наука, 1984. – 240 с. 19. Воскобойников Ю.Е., Мухина И.Н. Локальный регуляризирующий алгоритм восстановления контрастных сигналов и изображений // Автометрия. – 2000. – № 3. – С 45–53. 20. Морозов В.А. Регулярные методы решения некорректно поставленных задач. – М.: Наука, 1987. – 240 с. 21. Тихонов А.Н., Гончарский А.В., Степанов В.В., Ягола А.Г. Численные методы решения некорректных задач. – М.: Наука, 1990. – 232 с. 22. Сизиков В.С. Обобщенный метод редукции измерений // Электронное моделирование. – 1991. – Т. 13. – № 4. – С. 7–14. 23. Верлань А.Ф., Сизиков В.С., Мосенцова Л.В. Метод вычислительных экспериментов для решения интегральных уравнений в обратной задачи спектроскопии // Электронное моделирование. – 2011. – Т. 33. – № 2. – С. 3–12.
Сизиков Валерий Сергеевич Кривых Александр Владимирович
– Санкт-Петербургский национальный исследовательский университет информационных технологий, механики и оптики, доктор технических наук, профессор, sizikov2000@mail.ru
– Санкт-Петербургский национальный исследовательский университет информационных технологий, механики и оптики, аспирант, krivykh1987@mail.ru
28 Научно-технический вестник информационных технологий, механики и оптики,
2013, № 3 (85)
УДК 535.338.1+519.642.3+519.6
ВОССТАНОВЛЕНИЕ НЕПРЕРЫВНЫХ СПЕКТРОВ АДАПТИВНЫМ СПОСОБОМ ВЫЧИСЛИТЕЛЬНЫХ ЭКСПЕРИМЕНТОВ С РЕГУЛЯРИЗАЦИЕЙ
В.С. Сизиков, А.В. Кривых
Рассмотрена обратная задача спектроскопии – восстановление непрерывных спектров путем математической обработки измеренных спектров, искаженных аппаратной функцией спектрометра и помехами. Задача сводится к решению интегрального уравнения Фредгольма I рода. Задача его решения некорректна, поэтому для получения устойчивого решения используется метод регуляризации Тихонова. При этом применен адаптивный способ вычислительных экспериментов, согласно которому, наряду с исходным спектром P, обрабатывается модельный спектр Q с задаваемым истинным спектром z и моделируемым измеренным спектром u с учетом дополнительной (априорной) информации об истинном спек-
тре P. Это позволяет выбрать параметр регуляризации . Предложенная методика может быть использована для по-
вышения разрешающей способности спектрометра. Приведены численные иллюстрации. Ключевые слова: непрерывный спектр, обратная задача спектроскопии, интегральное уравнение, метод регуляризации Тихонова, способ вычислительных экспериментов, повышение разрешающей способности спектрометра.
Введение Измеренный спектрометром (например, интерферометром Фабри–Перо) спектр u() (где – длина волны) обычно отличается от истинного спектра z() [1–8]. Это проявляется, во-первых, в большей сглаженности спектра u() по сравнению с z() , а именно, в спектре u() не разрешены близкие линии, сглажена тонкая структура спектральной линии, что является результатом воздействия аппаратной функции спектрального прибора [1–9]. Во-вторых, это проявляется в зашумленности спектра u() , а именно, слабые линии «тонут» в шуме, что является результатом погрешностей измерений [1–3], а также воздействия среды, через которую проходит излучение [10]. Дадим следующее определение аппаратной функции (АФ) [3, 6–8] (ср. [9, С. 32, 704]): аппаратной функцией K(, ) спектрометра называется его реакция (в виде измеренной интенсивности) на дискретную линию единичной интенсивности и длины волны при настройке спектрометра на длину волны . Форма аппаратной функции (ширина и т.д.) может заметно меняться с изменением длины волны настройки [min , max ] , где [min , max ] – диапазон длин волн изучаемой части спектра. Обычно с увеличением АФ становится шире, что характерно для широкополосной спектрометрии, например, изучения спектра звезды во всем видимом диапазоне. Если же АФ практически не изменяется при изменении , то АФ является разностной (инвариантной): K(,) K( ) , что имеет место, например, при изучении тонкой структуры отдельной линии [3, 6, 8], когда диапазон [min , max ] мал. На рис. 1 в качестве примера приведен смоделированный непрерывный измеренный спектр u() , сглаженный аппаратной функцией спектрометра K(, ) , а также зашумленный (и дискретизированный) измеренный спектр u() u() u (где u – шум) и АФ спектрометра, причем, поскольку в дан-
22 Научно-технический вестник информационных технологий, механики и оптики,
2013, № 3 (85)
В.С. Сизиков, А.В. Кривых
ном примере K(, ) – функция неразностная, то приведено два ее «сечения» (подробности примера см.
дальше). В принципе похожий вид может иметь непрерывный узкополосный спектр [6, С. 200], например, сверхтонкая структура отдельной линии, обусловленная магнитными или электрическими полями (эффект Зеемана или Штарка), а также тепловым уширением (эффект Доплера) [10], однако в этом случае диапазон [min , max ] мал, а АФ – разностная: K(, ) K( ) .
6
Ииннттееннссииввнноосстть,ь,у.у.е.е.
5 u()
4 u~()
3
2
1 K(620, )
K(485, ) 0
450 460 470 480 490 500 510 520 530 540 550 560 570 580 590 600 610 620 630 640 650
min
, , нм
max
Рис. 1. u() – измеренный без шума спектр; u() – измеренный зашумленный и дискретизированный
спектр; K ( , ) – АФ при некоторой длине волны настройки ; [min , max ] – широкий диапазон длин волн
Как будет видно далее, в примере на рис. 1 в измеренном спектре u() (тем более, в зашумленном
спектре u() ) не разрешены близкие линии и не выявлены слабые, причем этот эффект тем сильнее, чем
шире АФ K(, ) (а также чем выше уровень шумов), другими словами, чем меньше разрешающая спо-
собность спектрометра [1, 9]. В данной работе ставится известная обратная задача спектроскопии – задача восстановления ис-
тинного спектра z() по измеренному спектру u() и аппаратной функции K(, ) [1–8, 11–15]. Дан-
ная задача описывается интегральным уравнением (см. дальше), задача решения которого некорректна, поэтому его обычно решают методом регуляризации Тихонова. При этом важным является вопрос о вы-
боре параметра регуляризации . В данной работе предлагается новый адаптивный способ (вычисли-
тельных экспериментов) для выбора параметра .
Математическая формулировка обратной задачи спектроскопии
Рассмотрим случай непрерывного спектра, обычно характерного для веществ с повышенной плотностью (расплавленный жидкий металл, плазма и т.д.). Измеренная интенсивность u() при настройке
спектрометра на длину волны равна сумме (интегралу) по всем истинным интенсивностям z() с ве-
совой функцией K:
b
u() z() K (, ) d , a
где a min , b max , откуда, варьируя значение (т.е. выполняя сканирование по спектру) и учитывая
зашумленность спектра u() , получим:
b
K (, ) z() d u() , c d ,
a
где [c, d] – пределы изменения (обычно более широкие, чем [a,b] ).
(1)
В соотношении (1) известны (измерены или заданы) u() , K(, ) , a, b, c, d, а z() является ис-
комым истинным спектром. Соотношение (1) есть интегральное уравнение Фредгольма I рода, причем
Научно-технический вестник информационных технологий, механики и оптики, 2013, № 3 (85)
23
ВОССТАНОВЛЕНИЕ НЕПРЕРЫВНЫХ СПЕКТРОВ АДАПТИВНЫМ СПОСОБОМ...
K(, ) является ядром уравнения, u() – правой частью, а z() – искомой функцией. Если K(,) K( ) , то
K ( ) z() d u() , 0 .
(2)
0
Соотношение (2) есть интегральное уравнение Фредгольма I рода типа свертки на полуоси. Реше-
ние уравнения (1) или (2) дает возможность, в принципе, восстановить истинный спектр z() . Однако
задача решения уравнений (1) и (2) является некорректной (существенно неустойчивой) [2–4, 6, 8, 16]: если решать уравнение (1), например, методом квадратур, а уравнение (2) – методом преобразования Фурье (инверсной фильтрации), то в качестве решения получим так называемую «пилу» [3, 6, 8] – крайне неустойчивое решение. По этой причине для устойчивого решения этих уравнений необходимо применение таких методов, как регуляризация Тихонова [2–4, 6–8, 11–16], параметрическая фильтрация Винера [3, 6, 8, 16], итеративная регуляризация Фридмана [6, 8, 16] и др.
При обработке спектра в широком диапазоне длин волн следует учитывать изменение формы АФ
K(, ) с изменением длины волны настройки . При обработке же спектра в узкой полосе следует ис-
пользовать уравнение Фредгольма I рода с разностным ядром (ср. (2)):
b
K ( ) z() d u() , c d .
(3)
a
Задача решения уравнений (1)–(3) связана с задачей редукции к идеальному спектральному при-
бору [1–4, 9, 17] – с одним из вариантов редукционной проблемы Рэлея [3, 6, 8, 13]. Успешное решение
задачи редукции позволит путем математической обработки результатов измерений повысить разре-
шающую способность спектрального прибора. В настоящей статье воспользуемся методом регуляриза-
ции Тихонова. Что касается других устойчивых методов (фильтрации Винера, итеративных методов и
др.), то они изложены в различных публикациях ([3, 6, 8, 16] и др.) и также могут быть применены для
устойчивого восстановления спектров.
Краткая формулировка метода регуляризации Тихонова
Запишем уравнение (1) в виде
b
A z K (, ) z() d u() , c d ,
(4)
a
где A – оператор, соответствующий ядру K. Метод регуляризации Тихонова сводится к решению инте-
грального уравнения Фредгольма II рода
b
z (t) B(t, ) z () d U (t) , a t b , a
где 0 – параметр регуляризации, а новое ядро и новая правая часть равны
(5)
dd
B(t, ) B(,t) K (,t) K (, ) d , U (t) K (,t)u() d . cc
В таком варианте уравнение (5) обычно решается методом квадратур [3, 4, 6, 8, 16]. Если же рассматривать уравнение (2) или (3), то его решение методом регуляризации Тихонова будет включать преобразование Фурье и -регуляризацию (подробности см. в [3, 4, 6, 8, 11–16, 18–21]).
При этом важным является вопрос о выборе параметра регуляризации и об учете дополнитель-
ной (априорной) информации относительно искомого спектра z() . Существует ряд способов выбора
параметра регуляризации : способ невязки, обобщенный принцип невязки, метод перекрестной значимости, локальный регуляризующий алгоритм, способ подбора и др. [3, 4, 6, 8, 13, 16, 18–21].
Способ вычислительных экспериментов
В данной работе получает дальнейшее развитие способ вычислительных экспериментов для выбора параметра регуляризации (другие его названия – способ псевдообратного оператора, способ эталонных, или модельных примеров, способ моделирования) [3, 6, 8, 16, 22, 23]. Данный способ учитывает дополнительную (априорную) информацию об искомом спектре (оценку количества спектральных линий, их параметров и т.д.) и поэтому является интерактивным и адаптивным способом.
Кратко изложим способ вычислительных экспериментов. Рассмотрим операторное уравнение I рода: A z u (ср. (4)). Полагаем, что вместо точных u и K из-
вестны u и K такие, что || u u || , || A A || , где и – верхние оценки погрешностей по норме
правой части u и ядра K. При использовании метода регуляризации Тихонова решается операторное
24 Научно-технический вестник информационных технологий, механики и оптики,
2013, № 3 (85)
В.С. Сизиков, А.В. Кривых
уравнение z AT A z AT u (ср. (5)), где AT – транспонированный оператор. Обозначим z z z –
погрешность регуляризованного решения z , а z – точное решение (нормальное псевдорешение [16, 20, 21]). В работах [16, 22, 23] получена следующая оценка относительной погрешности регуляризованного решения по норме:
||
z || z ||
||
()
,
(6)
где
()
|| A
2
||
p p 1
.
(7)
Здесь отн отн , причем отн || u || и отн || A || – относительные погрешности исходных
данных; p || A ||2 , A – псевдообратный оператор: Au z [20, С. 184]. Функция () является верх-
ней огибающей для истинной относительной погрешности
отн
()
||
||
z z ||
||
.
(8)
В работах [16, 22] показано, что функция () имеет (единственный) минимум при условии
p (|| A || )2 27 16 1, 69 или || A || || A || 3 3 4 1, 30 . Согласно соотношениям (6)–(8), оценка
относительной погрешности || z || || z || регуляризованного решения z зависит от A и (точнее, от
произведения || A || ). По этой причине, если решается несколько задач (другими словами, обрабатыва-
ется несколько спектров) с одинаковыми A и , то для них оценки погрешности
отн ()
||
z || z ||
||
|| A~ ||
2
p p 1
будут одинаковыми.
Отсюда следует, что при решении некоторого исходного примера P (т.е. при обработке исходного
спектра uP ) с неизвестным решением (спектром) zP можно использовать результаты решения другого,
модельного, примера Q с известным (заданным) точным решением (спектром) zQ , причем с такими же
A и , что и в примере P. При этом при решении примера Q можно рассчитать функцию
отн () Q || z Q || || zQ || и по ней найти опт Q – оптимальное значение , при котором
отн () Q
min
.
Это
значение
опт Q
может быть использовано при решении исходного примера (спек-
тра) P. При этом необходимо также определить p || A ||2 . Оценка p может быть получена путем подбо-
ра такого значения p, при котором огибающая кривая () касается набора кривых отн ()Q (см. рис. 3).
Добавим, что для повышения эффективности изложенного способа модельный пример Q (или несколько примеров) должен содержать дополнительную информацию об исходном примере (спектре) P, а именно, оценку количества спектральных линий (максимумов) в искомом спектре zP , соотношений их
интенсивностей и значений их длин волн. Данную оценку должен делать опытный спектроскопист. Ис-
пользование такой информации в модельном примере Q позволит более удачно выбрать параметр .
Данный способ следует считать адаптивным и интерактивным способом.
Численная иллюстрация
В рамках системы программирования MATLAB 7 был разработан пакет программ для восстанов-
ления истинных непрерывных спектров z() путем численного решения интегрального уравнения
Фредгольма I рода методом регуляризации Тихонова с использованием способа вычислительных экспериментов.
Сначала был рассмотрен первый пример (рис. 1) – оригинал P, у которого известен зашумленный
измеренный спектр u() на равномерной сетке узлов min , min h,, max , где min 450 нм ; max 650 нм ; h const 1 нм – шаг дискретизации; n (max min ) h 200 – число шагов дискретизации по . Известна также аппаратная функция – дифракционная АФ Рэлея (ср. [1, 4, 5]) вида
K (, )
1 ()
sinc2
()
1 ()
sin ( ) ()2
( ) ()
,
(9)
Научно-технический вестник информационных технологий, механики и оптики, 2013, № 3 (85)
25
ВОССТАНОВЛЕНИЕ НЕПРЕРЫВНЫХ СПЕКТРОВ АДАПТИВНЫМ СПОСОБОМ...
где () – полуширина АФ по уровню 0, равная приблизительно ширине АФ по уровню 0,5, которую мы
положили равной () 8 min нм . Спектр полагается широкополосным (от фиолетового до красного), поэтому ширина АФ непостоянна, а именно, (min ) 8 нм , а (max ) 11, 55 нм , т.е. (max ) (min ) 1, 44 . При этом истинный спектр z() в примере P неизвестен.
Из рис. 1 видно, что измеренный спектр u() имеет довольно сложную структуру, а именно, со-
держит шесть явных флуктуаций, две из которых (при 525 нм и 620 нм ), скорее всего, состоят каждая из двух линий, но они не разрешились из-за того, что АФ имеет немалую ширину и, тем самым, ограничивает разрешающую способность спектрометра. Кроме того, есть намек на то, что при 507 нм и 543 нм имеются еще две слабые линии. Таким образом, все указывает на то, что на самом деле в спектре имеются не менее восьми спектральных линий. В связи с этим в качестве второго (модельного или эталонного) примера Q был составлен близкий к оригиналу P пример, истинный спектр zQ () которого состоит из 9 спектральных линий в виде гауссиан:
zQ () 2, 0 exp{[( 486) /10]2} 0, 4exp{[( 512) / 5]2}
8,5exp{[( 522) / 2]2} 9, 2exp{[( 530) / 2]2}
0,5exp{[( 542) / 5]2} 8, 2exp{[( 566) / 6]2}
2,5exp{[( 592) / 4]2} 4,5exp{[( 614) / 7]2}
3, 0 exp{[( 626) / 5]2}.
Измеренный спектр uQ () в примере Q был рассчитан согласно выражению
b
uQ () K (, ) zQ () d , c d , a
численно. При этом a = 460, b = 640, c = 450, d = 650 нм. Погрешности измеренной uP () были оценены примерно в 1%, что соответствует среднеквадра-
тическому отклонению СКО ≈ 0,02. По этой причине к значениям uQ () были добавлены случайные
нормальные погрешности с СКО от 0,01 до 0,025, что соответствует отн 0,5 1,25% , поскольку значе-
ние отн в исходном примере P известно неточно. АФ спектрометра в примере Q была взята в виде (9),
причем (поскольку АФ известна также неточно) () было взято равным () (8 ) min , где 0 0,3 , что соответствует отн 0 3% .
Далее модельный пример Q был решен методом квадратур с регуляризацией Тихонова с помощью
разработанной m-функции Tikh.m [6, С. 207] для ряда значений параметра регуляризации , и была
построена зависимость относительной погрешности регуляризованного решения z () по отношению к точному решению z() (см. (8)):
отн ()
z () z() z()
.
На рис. 2 представлены зависимости отн () для ряда погрешностей отн и отн . На рис. 2 пред-
ставлена также огибающая () (см. (7)), при построении которой было положено 102 и
|| A |||| A || || u ||L2 || z ||L2 0,82 . Для ряда значений p от 100 до 270 были рассчитаны кривые ()
(рис. 2). Было выбрано то значение p, при котором одна из кривых касается набора кривых отн () , а
именно, p = 100. Этому соответствует значение параметра регуляризации 103 . Из рис. 2 видно, что,
несмотря на разброс кривых отн () и () , значения p и, как следствие, определяются уверенно.
При значении 103 , выбранном с помощью решения модельного примера Q как вспомогательного восстановлен спектр в исходном примере P (рис. 3). Как видно из рис. 3, в примере P разрешились близкие линии и восстановились слабые линии, правда, на краях спектра проявился эффект Гиббса, однако в слабой форме (на уровне погрешностей метода). Аналогичные результаты получены для других, весьма различных, непрерывных спектров [3, 6–8, 14, 22, 23], т.е. изложенная в работе методика вычислительных экспериментов может быть использована для широкого класса спектров (с близкими линиями, со слабыми линиями, узкими и широкими линиями и т.д.).
26 Научно-технический вестник информационных технологий, механики и оптики,
2013, № 3 (85)
В.С. Сизиков, А.В. Кривых
1
0,9
0,8
0,7
0,6
отн ()
0,5
0,4
p=270 p=200
() p=100
0,3
0,2
0,1
0–9 –8 –7 –6 –5 –4 –3 –2 –1 0 1 lg
Рис. 2. Зависимости отн () для ряда погрешностей отн и отн и огибающие ()
2
9 11
8
7
интеИннстиевннсоисвтнь,осу.теь., у.е.
6 5 22
4
3
2 33
1
0
460 470 480 490 500 510 520 530 540 550 560 570 580 590 600 610 620 630 640 , , нм
Рис. 3. 1 – истинный спектр zP () ; 2 – измеренный спектр u() ; 3 – восстановленный спектр z P ()
Заключение
Практическое использование изложенной методики позволит повысить разрешающую способность спектрометра. Спектральный прибор может быть соединен с компьютером или со спецпроцессором с заложенным в него математическим и программным обеспечением, реализующим методы и численные алгоритмы решения обратной задачи спектроскопии. В результате такого комплексирования (соединения прибора с компьютером) можно разрешить близкие и выделить слабые линии спектров излучения (или поглощения), а именно, в физике – спектров газов, жидкостей, металлов, плазмы; в астрофизике – спектров звезд, планет, галактик, туманностей, комет; в металлургии – спектров расплавленных металлов в домнах; в геофизике – спектров залежей руд, минералов, нефти, газа и т.д.
Работа выполнена при поддержке РФФИ (грант № 13-08-00442).
Научно-технический вестник информационных технологий, механики и оптики, 2013, № 3 (85)
27
ВОССТАНОВЛЕНИЕ НЕПРЕРЫВНЫХ СПЕКТРОВ АДАПТИВНЫМ СПОСОБОМ...
Литература
1. Раутиан С.Г. Реальные спектральные приборы // Успехи физических наук. – 1958. – Т. 66. – Вып. 3. – С. 475–517.
2. Кочиков И.В., Курамшина Г.М., Пентин Ю.А., Ягола А.Г. Обратные задачи колебательной спектроскопии. – М.: Изд-во МГУ, 1993. – 204 с.
3. Сизиков В.С. Математические методы обработки результатов измерений. – СПб: Политехника, 2001. – 240 с.
4. Старков В.Н. Конструктивные методы вычислительной физики в задачах интерпретации. – Киев: Наукова думка, 2002. – 264 с.
5. Fleckl T., Jäger H., Obernberger I. Experimental verification of gas spectra calculated for high temperatures using the HITRAN/HITEMP database // J. Phys. D: Appl. Phys. – 2002. – V. 35. – P. 3138–3144.
6. Сизиков В.С. Обратные прикладные задачи и MatLab. – СПб: Лань, 2011. – 256 с. 7. Сизиков В.С., Кривых А.В. Применение способа эталонных примеров при решении обратной задачи
спектроскопии методом регуляризации // Изв. вузов. Приборостроение. – 2011. – Т. 54. – № 9. – С. 44– 51. 8. Сизиков В.С. Интегральные уравнения и MatLab в задачах томографии, иконики и спектроскопии. – Saarbrücken: LAP, 2011. – 252 c. 9. Физический энциклопедический словарь / Гл. ред. А.М. Прохоров. – М.: Сов. энциклопедия, 1984. – 944 с. 10. Ландсберг Г.С. Оптика: Учебное пособие для вузов. – 6-е изд. – М.: ФИЗМАТЛИТ, 2003. – 848 с. 11. Преображенский Н.Г., Седельников А.И. Оптимизация спектроскопических измерений на основе методов регуляризации // Журнал прикладной спектроскопии. – 1981. – Т. 35. – Вып. 4. – С. 592–599. 12. Брагинская Т.Г., Клюбин В.В. Решение обратной задачи спектроскопии оптического смещения методом регуляризации Тихонова. Препринт № 855. – Л.: ЛИЯФ, 1983. – 60 с. 13. Глазов М.В., Болохова Т.А. Решение редукционной проблемы Рэлея с использованием различных модификаций метода регуляризации // Оптика и спектроскопия. – 1989. – Т. 67. – Вып. 3. – С. 533– 537. 14. Кривых А.В., Сизиков В.С. Комплексированное восстановление непрерывных спектров с использованием псевдообратной матрицы // XLI Неделя науки СПбГПУ: материалы научно-практической конференции с международным участием. Ч. XIII. – СПб: Изд-во Политехн. ун-та, 2012. – С. 240–242. 15. Кривых А.В., Сизиков В.С. Обработка дискретных спектров с помощью алгоритма интегральной аппроксимации // Научно-технический вестник СПбГУ ИТМО. – 2011. – № 5 (75). – С. 14–18. 16. Верлань А.Ф., Сизиков В.С. Интегральные уравнения: методы, алгоритмы, программы. – Киев: Наукова думка, 1986. – 544 с. 17. Краулиня Э.К., Лиепа С.Я., Пикалов В.В., Скудра А.Я. К проблеме исследования атомной сенсибилизированной флуоресценции по контурам спектральных линий // Некорректные обратные задачи атомной физики / Сб. статей под ред. Н.Г. Преображенского. – Новосибирск: ИТПМ, 1976. – 133 с. 18. Воскобойников Ю.Е., Преображенский Н.Г., Седельников А.И. Математическая обработка эксперимента в молекулярной газодинамике. – Новосибирск: Наука, 1984. – 240 с. 19. Воскобойников Ю.Е., Мухина И.Н. Локальный регуляризирующий алгоритм восстановления контрастных сигналов и изображений // Автометрия. – 2000. – № 3. – С 45–53. 20. Морозов В.А. Регулярные методы решения некорректно поставленных задач. – М.: Наука, 1987. – 240 с. 21. Тихонов А.Н., Гончарский А.В., Степанов В.В., Ягола А.Г. Численные методы решения некорректных задач. – М.: Наука, 1990. – 232 с. 22. Сизиков В.С. Обобщенный метод редукции измерений // Электронное моделирование. – 1991. – Т. 13. – № 4. – С. 7–14. 23. Верлань А.Ф., Сизиков В.С., Мосенцова Л.В. Метод вычислительных экспериментов для решения интегральных уравнений в обратной задачи спектроскопии // Электронное моделирование. – 2011. – Т. 33. – № 2. – С. 3–12.
Сизиков Валерий Сергеевич Кривых Александр Владимирович
– Санкт-Петербургский национальный исследовательский университет информационных технологий, механики и оптики, доктор технических наук, профессор, sizikov2000@mail.ru
– Санкт-Петербургский национальный исследовательский университет информационных технологий, механики и оптики, аспирант, krivykh1987@mail.ru
28 Научно-технический вестник информационных технологий, механики и оптики,
2013, № 3 (85)