CC BY
44
НАУЧНО-ТЕХНИЧЕСКИИ ВЕСТНИК ИНФОРМАЦИОННЫХ ТЕХНОЛОГИИ, МЕХАНИКИ И ОПТИКИ сентябрь-октябрь 2017 Том 17 № 5 ISSN 2226-1494 http://ntv.i1mo.ru/
SCIENTIFIC AND TECHNICAL JOURNAL OF INFORMATION TECHNOLOGIES, MECHANICS AND OPTICS September-October 2017 Vol. 17 No 5 ISSN 2226-1494 http://ntv.ifmo.ru/en
УДК 535.338.1, 519.6; 519.853.6
ИССЛЕДОВАНИЕ ПОГРЕШНОСТЕЙ НЕКОТОРЫХ МЕТОДОВ РАЗДЕЛЕНИЯ ПЕРЕКРЫВАЮЩИХСЯ СПЕКТРАЛЬНЫХ ЛИНИЙ В УСЛОВИЯХ ВОЗДЕЙСТВИЯ ПОМЕХ В.С. Сизикова, А.В. Лавров3 a Университет ИТМО, Санкт-Петербург, 197101, Российская Федерация Адрес для переписки: sizikov2000@mail.ru
Информация о статье
Поступила в редакцию 10.07.17, принята к печати 21.08.17 doi: 10.17586/2226-1494-2017-17-5-879-889 Язык статьи - русский
Ссылка для цитирования: Сизиков В.С., Лавров А.В. Исследование погрешностей некоторых методов разделения перекрывающихся спектральных линий в условиях воздействия помех // Научно-технический вестник информационных технологий, механики и оптики. 2017. Т. 17. № 5. С. 879-889. doi: 10.17586/2226-1494-2017-17-5-879-889
Аннотация
Предмет исследования. Рассмотрена одна из актуальных задач спектроскопии - разделение (сепарация) близких спектральных линий. Метод. Задача решена математическим и компьютерным путем. Выполнена минимизация функционала невязки между измеренным и рассчитанным спектрами. При этом линии-компоненты смоделированы гауссианами, и задача сводится к отысканию их параметров. Основной результат. Для минимизации функционала предложена модификация метода координатного спуска с использованием способа сужающихся ограничений. Для сглаживания и дифференцирования зашумленных экспериментальных спектральных данных предложено использовать сплайны. Разработано программное обеспечение на MatLab и выполнена обработка ряда спектров. Практическая значимость. Разработанная методика может быть использована для восстановления тонкой структуры спектров и, тем самым, для повышения разрешающей способности спектрометров. Ключевые слова
разделение (сепарация) линий спектра, минимизация функционала невязки, измеренный и рассчитанный спектры, гауссианы, метод координатного спуска с ограничениями, сплайны, программное обеспечение, MatLab Благодарности
Работа выполнена при поддержке РФФИ (грант № 13-08-00442).
STUDY OF ERRORS OF SOME METHODS FOR SEPARATING OVERLAPPED SPECTRAL LINES UNDER NOISE EFFECT V.S. Sizikov3, A.V. Lavrov3
3 ITMO University, Saint Petersburg, 197101, Russian Federation Corresponding author: sizikov2000@mail.ru
Article info
Received 10.07.17, accepted 21.08.17 doi: 10.17586/2226-1494-2017-17-5-879-889 Article in Russian
For citation: Sizikov V.S., Lavrov A.V. Study of errors of some methods for separating overlapped spectral lines under noise effect. Scientific and Technical Journal of Information Technologies, Mechanics and Optics, 2017, vol. 17, no. 5, pp. 879-889 (in Russian). doi: 10.17586/22261494-2017-17-5-879-889
Abstract
Subject of Research. We consider one of the actual problems of spectroscopy, that is, the separation of close spectral lines. Method. The problem is solved by a mathematical (computer) way, namely, by minimizing a functional of discrepancy between the measured and calculated spectra. In this case, the lines (components) are modeled by Gaussians and the problem is reduced to localization of their parameters. Main Results. To minimize the functional we propose the coordinate descent method modification with the use of decremental constraints technique. To smooth and differentiate noisy experimental spectral data, we suggest using splines. The software on MatLab is developed and a number of spectra are processed. Practical Relevance. The developed technique can be used to restore the fine structure of spectra and, thereby, to increase the resolving power of spectrometers. Keywords
spectrum lines separation, discrepancy functional minimization, measured and calculated spectra, Gaussians, coordinate
descent method with constraints, splines, software, MatLab
Acknowledgements
This work was supported by the Russian Foundation for Basic Research (RFBR), grant No. 13-08-00442.
Введение
В ряде работ ([1-8] и др.) рассматривается характерный для спектроскопии случай, когда отдельные линии перекрываются в той или иной степени вследствие их близости или из-за наличия магнитного или электрического поля, в результате чего каждая линия расщепляется на ряд перекрывающихся компонент (эффект Зеемана или Штарка) и создается сверхтонкая структура линии. Эффект перекрытия спектральных линий в спектре имеет место также, например, для колебательных спектров вещества в конденсированном состоянии [7]. Перекрытие создается за счет того, что взаимодействие молекул в жидкости ведет к уширению линий в спектре и, как следствие, к их частичному наложению.
Возникает актуальная задача разделения (сепарации) отдельных линий-компонент путем математической и компьютерной обработки спектра. Данную задачу сформулируем следующим образом.
Формулировка задачи
Пусть дан суммарный измеренный спектр 2(Х), где X - длина волны, равный сумме N отдельных линий-компонент с профилями интенсивности 2, (А), , = , т.е.
N
2 (А) = £ 2, (А), а <А< Ь , (1)
j=l
где , - номер линии, [а,Ь] - пределы суммарного спектра (ограниченная протяженность интервала регистрации). Требуется на основе суммарного (обычно зашумленного) измеренного спектра 2(Х) восстановить профили отдельных линий-компонент 2, (А), а также определить их число N. При этом возможно
моделировать профили линий некоторыми функциями (гауссианами, лоренцианами и т.д.) и использовать дополнительную информацию о линиях (начальные приближения параметров линий, диапазоны изменения параметров и т.д.).
Строго говоря, данная задача является некорректной (может не существовать решения, решение может быть неединственным или неустойчивым) [9-12]. Однако за счет использования дополнительной информации задача может стать промежуточной между корректными и некорректными задачами [12] и может иметь корректное решение, даже если не использовать, например, метод регуляризации Тихонова [9-11].
Оценка числа компонент
Важным и изначальным является вопрос о числе компонент-линий N, входящих в суммарный спектр 2(Х). Во многих публикациях ([1, 4, 6, 7] и др.) для оценки числа компонент N (а также координат максимумов линий А) используется метод производных суммарного спектра обычно 1-го,
2-го и 4-го порядков (т.е. I = 1, 2 или 4). Этот метод был предложен в работе [1], где была использована производная лишь 1-го порядка, а затем метод распространен на производные различных порядков [4, 6, 7].
Суть метода производных или дифференцирования состоит в следующем.
Моделируем каждую линию 2, (А) гауссианой1:
2 (X) = А ехр
Г (А-А,.)2 Л ' 2 ст2
(2)
где А, - амплитуда линии, А, - координата максимума (ее центральная длина волны), с, - среднеквадра-тическое отклонение (СКО), причем т, = %/21п2 ст, = 1,1773 ст, - полуширина гауссианы по относительному уровню 0,5 [6, 13, 14].
Производная 1-го порядка от 2 (А) по X равна (ср. [1])
2' (А) = - А, ехр
(а-а, )21 (а-а,
2ст2
(3)
а производная 2-го порядка равна (ср. [6])
2" (А) = - 4ехр
(А-А)2)( (А-А)
2ст?
- +1
(4)
Из формул (3) и (4) следует, что
1 Изложенная далее методика может быть распространена также на случай моделирования линии лоренцианой, экспонентой и т.д.
zj (xy) = 0, (5)
zj (X..) = - Ay/ст) , (6)
т.е. нулевое значение 1-й производной указывает на Xy (на максимумy-й компоненты-линии), а отрицательный экстремум 2-й производной также указывает на линию. При этом можно получить оценку амплитуды Ay (при наличии оценки стy):
Ay =a2y\zj (Xy )|. (7)
В работах [6, 7] и др. используется также производная 4-го порядка (ее положительные экстремумы указывают на линии). При этом в работе [6] предложена оригинальная методика, согласно которой вычисляются свертки 2-й и 4-й производных с отдельно взятым пиком (линией) типа (2). Это позволило оценить ст. и Ay.
Однако при обработке реального (измеренного) спектра необходимо численно дифференцировать его. При этом спектр обычно зашумлен, и его численное дифференцирование, тем более вычисление производной 4-го порядка, будет отягощено значительными погрешностями. Можно отфильтровать шумы, например, фильтром Савицкого-Голея (с помощью m-функции sgolayfilt.m в системе MatLab) [15, С. 193]. Можно также аппроксимировать спектр сглаживающим сплайном [16, С. 42; 17, С. 223; 18; 19] и дифференцировать сплайн или использовать метод производных без применения производной 4-го порядка.
Заметим также, что на практике отдельные линии zy (X) не могут быть дифференцированы в силу их недоступности, а может быть дифференцирован лишь суммарный спектр z(X) согласно (1). В результате этого параметры отдельных линий (в первую очередь Xy) будут определяться по производным (см. (3)-(7)) с погрешностями. Поэтому нужно использовать далее более точный метод сепарации линий (Нелдера-Мида и др.), при этом значения Xy, A., ст. и N, оцененные по производным, полагать в качестве начальных приближений.
Следующие методы могут быть использованы для более точной сепарации линий: симплекс-метод Нелдера-Мида [10, 20, 21], метод Фурье-самодеконволюции спектра [3, 7], метод «полуслепой» спектральной деконволюции [8], метод генетических алгоритмов [5, 7, 21, 22] и др. Метод Нелдера-Мида нами использован далее. Что касается метода Фурье-самодеконволюции [3], то в нем для разрешения перекрывающихся линий используется аподизация - усечение интерферограммы, по которой с помощью преобразования Фурье вычисляется спектр в Фурье-спектрометрах [23]. За счет аподизации ширины линий искусственно уменьшаются (до 5 раз), т.е. истинные профили линий-компонент искажаются ради их разрешения.
В данной работе для начальной оценки числа компонент N и координат их максимумов X (а также амплитуд A) используются производные суммарного спектра z лишь 1-го и 2-го порядков. При этом для повышения точности дифференцирования в случае зашумленности спектра предлагается использовать его сплайн-аппроксимацию. Затем параметры линий-компонент уточняются путем минимизации функционала невязки модифицированным методом координатного спуска (с использованием сужающихся ограничений на параметры) и для сравнения методом Нелдера-Мида (с заданием начальных приближений параметров). При этом аподизация и регуляризация не используются1.
Прямая задача
Начнем с прямой задачи - формирования или моделирования перекрывающихся линий-компонент, их суммарного спектра и добавления шумов.
Полагаем, что спектр z(X) есть сумма N линий-компонент согласно (1). Моделируем каждую линию Zy (X) гауссианой согласно (2). Полагаем, что суммарный спектр z(X) может быть зашумлен:
Z(x) = z (X) + noise, (8)
где noise - адаптивный шум2. На рис. 1 приведен случай трех линий в виде гауссиан и суммарный спектр незашумленный и с шумом.
1 Не используется также аппаратная функция спектрометра, дополнительно сглаживающая спектр. Ее влияние, связанное с решением интегрального уравнения, будет рассмотрено в отдельной публикации.
2 Шум может быть гауссов, рэлеевский, равномерный и т.д. Однако мы будем формулировать методику обработки зашумленного спектра, чувствительную не столько к типу шума, сколько к таким его характеристикам, как математическое ожидание (МО) и, особенно, СКО Ъг. Это будет достигнуто за счет использования сглаживающих сплайнов и ограничений на решение.
7
нТ 6 л
¡а о
О „• 4
2 В
2 л
Is.3
о
550
560
ч
/ / ч 4
/
t 1
.................... J X 3
/ \ / /
□X / у' \ V п
' р п U i 1 П п ~ 1
570
580
А, нм
590
600
Рис. 1. Исходные спектры. 1, 2, 3 - истинные перекрывающиеся линии (гауссианы) zi(A), z2(A), Z3(A);
4 - суммарный спектр z(A) = zi(A) + Z2(A) + Z3(A) - незашумленный (непрерывная линия) и зашумленный
(квадратики) 5-процентным шумом (5z = 0,1)
В примере на рис. 1 параметры трех линий A., A ., a., j = 1,N , где N = 3, объединены в единый вектор p длиной 3N = 9:
p = [A1, А1, a1, Аз, А2, a2, A3, А3, a3]. (9)
Точные (exact) параметры трех линий заданы равными pe = [3,7; 565; 3,6; 6,1; 574; 4,2; 2,8; 583; 3,9], причем значения А в (1)-(4), (8) заданы на дискретной сетке узлов: Ai = a + h(i - 1), где i = 1,n - номер узла дискретизации, a = Amin = 550 нм, b = Amax = 600 нм, h = ДА = 1 нм - шаг сетки, n = (b - a)/h + 1 = 51 -число узлов, приходящихся на интервал регистрации [a, b], а на контур каждой линии в пределах [А . - 3a ., А . + 3a . ] приходится соответственно 21, 37 и 23 отсчета.
Суммарный спектр z(A) рассчитан по формуле (1) при А = Ai, i = 1,n . Рис. 1 показывает, что отдельные линии-компоненты z. (А) в суммарном спектре z(A) проявляются слабо. К значениям z(A) были добавлены случайные нормальные погрешности 5z (см. (8)) с СКО, равным 8z = 0,1, что соответствует примерно 5% погрешности z (довольно значительная погрешность).
В работе [6] введено удобное понятие степени наложения пиков-линий (в наших обозначениях): D = т/Л , (10)
где т - полуширина пиков, а Д - расстояние между соседними пиками. Однако формула (10) применима лишь для случая, когда полуширины соседних пиков одинаковы. А для случая, вообще говоря, разных полуширин пиков предлагаем следующую формулу, как обобщение формулы (10):
(11)
D
j j+1
/л .
j j+1 / j j+1
J = 1, N - 1 ,
где т,,,+1 = (т, + т,+^/2 - средняя полуширина соседних (, и,+1) пиков, а Л,,+1 = А,.+1 - А,. - расстояние
между соседними пиками.
Однако формула (11) применима, строго говоря, для случая одинаковых амплитуд соседних пиков: А, = А,+1. Предлагаем следующую формулу для случая, вообще говоря, разных А, и А,+1 как обобщение формулы (11):
dj j+1 =■
j j+1
Л
j j+1
! ¡A. - A..+1 1 + —---—
. Aj + Aj+1 )
J = 1, N - 1.
(12)
Если А, = А,+1, то формула (12) переходит в формулу (11). А если, например, А,+1 << А , то
О,.+1 = 2 Т,.+1!Л,.+1 . Формула (12) показывает, что чем шире пики, чем они ближе друг к другу и чем
больше отличаются их амплитуды, тем больше степень наложения пиков О и тем сложнее будет разделить пики при решении обратной задачи (см. далее). Используем также величину
1 N - 1
D=N-Г £ dj j+■
(13)
среднее значение степени наложения.
В таблице в варианте 1 приведены данные прямой задачи по рис. 1.
0
Номер вари- анта Тип линий, величина шума 5г Точные параметры трех линий Степени наложения О23 и О (см. (12), (13)) Параметр Р сглаживания спектра г(Х) сплайном На что указывают производные г'(Х) и г''(Х) после сглаживания г(Х) сплайном Число итераций I и погрешность Ър
в методе координатного спуска (КС) (^ задает нач. шаг) в методе Нелдера-Мида (НМ)
1 Линии -гауссианы, 5г = 0,1 (5% шум) 3,7; 565; 3,6; 6,1; 574; 4,2; 2,8; 583; 3,9 (см. рис. 1) 0,635; 0,726; 0,681 (см. рис. 1) 0,4 (см. г(Х) на рис. 2, а) на X = 574; на X = 563, 573 и 583 (см. рис. 2) I = 18, Ър = 0,0103 при 5 = 100 (рис. 3); I = 27, Ър = 0,0137 при 5 = 10 I = 1098, Ър = 0,0697 (см. рис. 4)
2 Линии -гауссианы, 5г = 0,05 (2,5% шум) 3,7; 566; 3,6; 6,1; 574; 4,2; 2,8; 582; 3,9 0,714; 0,817; 0,766 0,5 на X = 573; на X = 565, 574 и 583 I = 18, Ър = 0,0433 при 5 = 100; I = 27, Ър = 0,0408 при 5 = 10 I = 1245, Ър = 0,0685
3 Линии -гауссианы, 5г = 0,1 (5% шум) 3,7; 565; 3,6; 3,5; 574; 4,2; 2,8; 583; 3,9 0,524; 0,589; 0,556 0,2 на X = 566, 573; на X = 564, 574 и 583 I = 18, Ър = 0,0075 при 5 = 100; I = 27, Ър = 0,0122 при 5 = 10 I = 1146, Ър = 0,0718
4 Линии -гауссианы, 5г = 0,05 (2,5% шум) 3,7; 567; 3,6; 3,5; 574; 4,2; 2,8; 581; 3,9 (см. рис. 5) 0,674; 0,757; 0,715 (см. рис. 5) 0,5 (см. г(Х) на рис. 6) на X = 569; на X = 565, 575 и 580 (см. рис. 6) I = 18, Ър = 0,0351 при 5 = 100 (рис. 7); I = 27, Ър = 0,0356 при 5 = 10 I = 1068, Ър = 0,0445 (см. рис. 8)
Таблица. Зависимость сепарации трех линий от различных факторов
Обратная задача
Под обратной задачей подразумевается задача сепарации (разделения) линий. Обратная задача является более сложной и важной, чем прямая задача.
Сначала оценим число компонент N в примере на рис. 1, для чего используем метод производных. Численно, например, с помощью т-функции Шй\т находим производные 1-го порядка г' (X) и 2-го порядка г "(X) по X суммарного спектра г^). На рис. 2, а, приведен спектр г^), зашумленный 5-процентным шумом. Первая производная г ' (X) от зашумленного спектра г(к) дает очень неустойчивый результат (пунктир на рис. 2, б). Еще более неустойчивый результат дает вторая производная г"(Х) (пунктир на рис. 2, в).
Чтобы сделать процедуру дифференцирования устойчивой, мы аппроксимировали зашумленный спектр Г(Х) сглаживающим кубическим сплайном [16, 17] с помощью т-функции СБарБ.т (ср. [18, 19]) и затем дважды дифференцировали сплайн. Заметим, что степень сглаживания сплайна регулируется параметром сглаживания Р 6 [0, 1], причем при Р = 0 сплайн получается максимально гладким, а при Р = 1 сглаживание отсутствует, и сплайн становится интерполирующим.
Первая производная г' (X)
0,
Исходная функция г(Х)
7 6 5 4 3 2 1 0
550 560
0,6 0,4 0,2 0 -0,2 -0,4 -0,6
: II ¿н
11 1 1\ '.......
т
У 1 п г Я', 11
V 1 т
1 1 " !
V |/|
570 580 590 600 X, нм а
550 560
570 580 590 600
X, нм
б
0,15 0,1 0,05 0
-0,05 -0,1 -0,15 -0,2 -0,25
Вторая производная г" (xx)
550 560
570 580 590 600
X, нм
Рис. 2. Метод производных. Исходная функция ¿(Л) незашумленная (непрерывная линия), зашумленная (квадратики) и сглаживающий сплайн (точки) (а); 1-я производная ¿(Л) от зашумленной функции ¿(Л) (пунктирная линия) и от сплайна (непрерывная линия) (б); 2-я производная ¿ (Л) от зашумленной ¿(Л) (пунктирная линия) и от сплайна (непрерывная линия) (в)
в
Данный параметр Р близок к 1 - а, где а - параметр регуляризации [17, С. 223]. Чем больше Р (а, значит, чем меньше а), тем меньше гладкость сплайна, и наоборот. Параметр Р можно выбирать способом подбора такого значения Р, при котором производные г '(А) и г"(А) имеют умеренную гладкость, причем 1-я производная дает четкий нулевой отсчет (или несколько отсчетов), а 2-я производная дает отрицательные экстремумы в количестве ~ N где N - предполагаемое количество линий.
На рис. 2, а, точками представлен сглаживающий сплайн, а на рис. 2, б, в, непрерывными линиями - производные г "(А) и г" (А) от сплайна. При этом способом подбора выбрано умеренное значение параметра сглаживания: Р = 0,4. Первая производная от сплайна имеет одно нулевое значение (см. (5)) при X = 574 нм, а вторая производная - три отрицательных экстремума (см. (6)) при X = 563, 573 и 583 нм. Сделан вывод, что число линий N = 3 (см. столбец 6 варианта 1 в таблице), и их центральные длины волн оценены как А, = 563; 573,5 и 583 нм.
Уточнение параметров линий. После получения оценки количества компонент N на основе 1-й и 2-й производных переходим к уточнению параметров компонент.
Компоненты г. (А), , = , суммарного спектра г(А) моделируем гауссианами согласно (2), у каждой - по три искомых параметра: амплитуда А., координата максимума А, и СКО с,..
Далее 3N параметров объединяем в единый вектор р длиной 3N по образцу (9). Параметры р3,
3 = 1,3N, находим путем минимизации функционала невязки:
Р = £ 6 - ) (14)
I=1
с ограничениями на параметры р3 в виде
Pmin/ ^ Pj ^ PmaaxJ, J = UN , (15)
где zt = z(Ai) - измеренные значения суммарного спектра, а zi = zt (p) - рассчитанные значения суммарного спектра согласно (1) при z(A) = z(Ai), а также согласно (2) при z(А) = z} (Ai) и согласно (9).
Для минимизации функционалов с ограничениями существует ряд методов: проекции градиента, условного градиента, сопряженных градиентов, наискорейшего спуска, оврагов, хорд, градиентного спуска и др. [10, 21, 24]. В данной работе для минимизации функционала (14) мы использовали разработанную нами модификацию метода координатного спуска [21, 25, 26]. В модификации метода координатного спуска (the method of coordinate descent) не используются градиенты, а вводятся ограничения на решение вида (15). Ограничения не позволяют выходить решению за пределы «коридора», даваемого неравенствами (15), тем самым, обеспечивая устойчивость и сходимость решения.
Приводим псевдокод модифицированного метода координатного спуска с использованием ограничений вида (15). Псевдокод оформлен в стиле m-функции системы MatLab, но не является m-функцией, а дает изложение метода, удобное для компьютерной реализации (в принципе, на любом языке программирования).
function [p, Fmin, Niter] = CoDesc(F, pmin, pmax, e, s, z, x) % F - минимизируемый функционал,
% pmin, pmax - массивы левых и правых ограничений для параметров, % e - задаваемая погрешность определения параметров (массив), % s - количество начальных шагов между pmin и pmax, % z, x - формальные параметры, передаваемые в F
% (в данном случае z - экспериментальный спектр z^), а x - сетка значений А) N ^ length(pmin);
h ^ (pmax - pmin) / s; d ^ h; % h и d - массивы начальных шагов p ^ (pmin + pmax) / 2; % массив начальных параметров % Итерации: Niter ^ 0;
while (any( | d | > e / 9 )) % пока хотя бы один из шагов больше предельного for J = 1 to N do % J - номер параметра
Niter ^ Niter + 1;
f ^ F(p, z, x); % вычисление функционала while ( | hj | > e. / 3 )
p.j ^ p. + h/;
if (p. < pminj) then % в случае выхода за левое ограничение
p. ^ pmin^
hj ^ - hj / 3; % уменьшение шага и смена его знака continue; % переход на while end if
if (pj > pmaxj) then % в случае выхода за правое ограничение
Pj ^ pmaxj;
hj ^ - hj / 3; % уменьшение шага и смена его знака continue; % переход на while end if
g ^ f % сохранение предыдущего значения функционала f ^ F(p, z, x); % вычисление нового значения функционала if (f > g ) then hj ^ - hj / 3; % уменьшение шага и смена его знака end if end while
hj ^ sign(hj) х d/; % восстановление шага со знаком end for
h ^ h / 9; d ^ h; % уменьшение всех шагов end while Fmin ^ f; return
Заметим, что не всегда можно задать достаточно узкий «коридор». В этом случае предлагается задавать сначала широкие ограничения (15), а затем их итеративно сужать (снова обращаясь к функции CoDesc), но так, чтобы решение не выходило из «коридора» (15). Данный способ называется способом сужающихся ограничений (the way of decremental constraints) [26]. В работе [26] он продемонстрировал свою эффективность в обратной задаче геофизики (гравиметрии).
Ниже мы использовали модифицированный метод координатного спуска в соединении со способом сужающихся ограничений для уточнения параметров линий.
Численные иллюстрации
В рамках системы программирования MatLab были разработаны программы для сепарации перекрывающихся линий. Выше в прямой задаче (а также на рис. 1, 2 и в начале таблицы) уже приведены исходные данные варианта 1 в дискретном виде и результаты дифференцирования суммарного спектра z(X).
Вариант 1 (а также другие варианты) описывает, например, физический процесс расщепления одиночной спектральной линии на зеемановский триплет в магнитном поле [27, С. 199], или расщепление линии в электрическом поле (эффект Штарка) [27, С. 857], или масс-спектр как суперпозицию пиков гауссовой формы [6], или эффект перекрытия линий колебательного спектра вещества в конденсированном состоянии [7], когда взаимодействие молекул ведет к уширению, например, трех линий и, как следствие, к их частичному наложению.
Далее выполнено уточнение параметров трех линий-гауссиан путем минимизации функционала (14) модифицированным методом координатного спуска (КС) с ограничениями типа (15): pmin = [3,3; 560; 3,3; 5,6; 572; 4; 2,6; 580; 3,7], pmax = [4; 568; 3,9; 6,5; 576; 4,3; 3,1; 586; 4,3]. Методом КС выполнено 18 итераций и получено следующее решение: p = [3,715; 565,0; 3,666; 6,068; 574,0; 4,163; 2,862; 582,9; 3,880]. При этом достигнуто следующее значение функционала: F = 0,457.
Среднеквадратическая погрешность решения p вычислена по формуле [26]
" 3N
' (16)
Sp =
J 3N "I1/2
— X wj (PJ " pj)
3N j=1
где PJ - вычисленные значения параметров линий, pJ - точные параметры (известные лишь при обработке модельных спеKтров), а WJ - вес^ равные ^ = 1 рт, J , прИЧем рти J = (ртт J + ртах J )/2 . Веса ^ введены потому, что искомые параметры (амплитуда, длина волны, полуширина линии) имеют, вообще говоря, различную физическую размерность и различный порядок величины, а введение весов делает слагаемые wJ р - )2 в (16) близкими по величине и безразмерными.
Погрешность Ър решения р методом КС по формуле (16) получилась равной 0,0103 при 5 = 100 и 0,0137 при 5 = 10.1 На рис. 3 и в таблице в варианте 1 представлен результат восстановления трех линий-гауссиан модифицированным методом КС.
1 В методе КС число s задает начальные шаги: hj = (pmaxj - pminj)/s.
Для сравнения было выполнено уточнение параметров трех линий-гауссиан путем минимизации функционала (14) также симплекс-методом Нелдера-Мида (НМ) [7, 20, 21], реализованным в т-функции 1т1П8еагсЬ.т [28, С. 399]. Специфика метода НМ состоит в том, что он не использует ограничения, как модифицированный метод КС, а использует начальные приближения для параметров. Кроме того, метод НМ (как и КС) не требует вычисления градиентов.
В методе НМ использовано начальное приближение для вектора р согласно (9) в виде: р = [3; 563; 3,3; 6; 573,5; 4,4; 3; 583; 3,8]. Методом НМ было выполнено 1098 итераций и получено уточненное решение: р = [4,053; 565,3; 3,728; 5,981; 574,0; 3,755; 3,215; 582,5; 3,890]. При этом Р = 0,372, а погрешность решения 5р = 0,0697 (примерно в 6-7 раз больше, чем методом КС). На рис. 4 и в таблице в варианте 1 представлен результат восстановления трех линий-гауссиан методом НМ.
Если же начальное приближение дляр в методе НМ положить равнымр = [3,5; 567; 3,5; 6; 573; 4,4; 3; 583; 3,8], а именно, приблизить р к точному решению, то получим более точное восстановление: 5р = 0,0314, I = 1237. Однако такое начальное приближение на практике трудно получить. Методу КС эта зависимость погрешности 5р от начального приближения менее свойственна, так как в нем задается сужающийся диапазон значений параметров (15).
7 6 5 4 3 2 1
л
&
и
<и
С
о
и
3 со
нн В
й Ё О
&
V
о
о
й
си
/V" // \\
Г ! г / / ч \ \\
/ \ / \ * \
X г
/ \ / » / V % \
1 5 « ' Г / V 4 Ьа а ■ ■ [
550
560
570
580
590
600
X, нм
Рис. 3. Нахождение параметров трех линий-гауссиан модифицированным методом координатного спуска.
Непрерывная линия - точный суммарный профиль г(Л), квадратики - использованный зашумленный профиль, пунктирные линии - найденные профили трех линий, кружки - суммарный найденный профиль,
штрих-пунктир - точные профили линий
7 6 5 4 3 2 1
А
£
И
<1>
С
о <и
А) и
м
нн В
Й Ё О
&
V
о
о
й
си
/ " ч / >• // VI А >1 \
/7 \\ л.
ч\ ¡1 \\ \а
ж
// \
/ У / \\ N V а
1 □ и 1 1 1 1 " | " п ~ 1
550
560
570
580
590
600
X, нм
Рис. 4. Нахождение параметров трех линий-гауссиан путем минимизации функционала Г методом Нелдера-Мида. Обозначения такие же, как на рис. 3
Были решены также другие примеры-варианты (см. таблицу и рис. 5-8).
Пояснения к таблице. Вариант 2 в таблице отличается от варианта 1, рассмотренного выше, меньшим расстоянием между линиями А и, как следствие, большей степенью наложения О, в результате чего погрешность 5р методом КС в варианте 2 увеличилась по сравнению с вариантом 1.
Вариант 3 отличается от варианта 1 тем, что амплитуда А2 средней линии близка амплитудам Л\ и А3, в результате степень наложения О уменьшилась и погрешность 5р методом КС понизилась.
Вариант 4 характерен еще большим сближением линий и еще большей степенью наложения О, в результате погрешность 5р стала больше, чем в вариантах 1 и 3.
Замечание. В варианте 1 (см. таблицу) погрешность ър методом КС равна ~ 0,01, а методом НМ ~ 0,07, т.е. примерно в 7 раз больше. В то же время значение функционала невязки Р методом КС равно ~ 0,45, а методом НМ ~ 0,37, т.е. меньше, чем методом КС. Этот неожиданный результат объясня-
0
0
ется тем, что при решении некорректных задач (а задача сепарации, строго говоря, является некорректной) малость невязки не является, вообще говоря, критерием малости погрешности решения [10, С. 224].
л &
ё и ° «
2 В ^ £
О с
X о
4,5 4
3,5 3
2,5 2
1,5 1
0,5 0
550
1
/ , V
1 2 \
/ \ /
V л 3
/ \
- / ''' -
у у V '■■.................. N Ч,
560
570 580
X, нм
590
600
Рис. 5. Исходные спектры (гауссианы) в варианте 4 (см. таблицу) Исходная функция z(X) Первая производная z'(X) Вт°рая пр°говодная г"(Х)
4,5 4
3,5 3
2,5 2
1,5 1
0,5
5°50
0,6 0,4 0,2 0 -0,2 -0,4
0,1 0,05 0
-0,05 -0,1 -0,15
560 570 580 590 600 X, нм
550 560
570 580 590 600
X, нм
б
550 560
570 580 590 600
X, нм
Рис. 6. Метод производных в варианте 4 (см. таблицу). Исходная функция г(А) незашумленная (непрерывная линия), зашумленная (квадратики) и сглаживающий сплайн (точки) (а); 1-я производная г'(А) от зашумленной функции г(А) (пунктирная линия) и от сплайна (непрерывная линия) (б); 2-я производная 2"(А) от зашумленной г(А) (пунктирная линия) и от сплайна (непрерывная линия) (в)
4,5 4
3,5
л &
и
С ^ ' о £ 3
9 "2,5
В я
В ° 2 н & 2
£ 1,5 1
0,5 0
й си
я*
I'' / -\ /;' ч //'
\,у \\ /'.
X «
/ » \ а \
/ / \ !
/ / V' \\
ч
550
560
570
580
590
600
X, нм
а
в
Рис. 7. Нахождение параметров трех линий-гауссиан модифицированным методом координатного спуска
в варианте 4 (см. таблицу)
Важными являются также вопросы о влиянии на сепарацию ошибочности оценки числа линий N о моделировании линий лоренцианами и другими функциями, о дополнительном сглаживании суммарного спектра под воздействием аппаратной функции спектрального прибора [13, 14, 29] и др. Эти вопросы могут быть предметом изучения в отдельных статьях.
550 560 570 580 590 600
А, нм
Рис. 8. Нахождение параметров трех линий-гауссиан симплекс-методом Нелдера-Мида (НМ) в варианте 4 (см. таблицу)
Заключение
В работе рассмотрен вопрос о разделении (сепарации) близких спектральных линий на основе измеренного зашумленного суммарного спектра z(А). Число линий-компонент N оценено по производным 1-го и 2-го порядков z'(А) и z''(А). При этом зашумленный спектр z(А) аппроксимирован сглаживающим
сплайном, что существенно снизило погрешность дифференцирования z(А).
Далее линии-компоненты смоделированы гауссианами. Параметры линий определены путем минимизации функционала невязки между измеренным и рассчитанным спектрами z(А). Для минимизации функционала использован модифицированный метод координатного спуска (КС) с сужающимися ограничениями. Разработаны программы для среды Ма1ЬаЪ, выполнена сепарация ряда спектров. Произведено сравнение результатов методом КС с результатами методом Нелдера-Мида (НМ), показавшие, что погрешность определения параметров линий методом КС в среднем меньше погрешности методом НМ.
Разработанная и представленная в работе методика может быть использована для разрешения близких линий в спектре, для восстановления тонкой структуры отдельных линий и т.д., тем самым, для повышения разрешающей способности спектрометров за счет математической и компьютерной обработки спектров.
д а н и
(D
О D (D
Й рц
4,5 4
3,5 3
м о Ç
эт 2,5 Я '
о 2 & 2
1,5 1
0,5 0
\
! \ 4 4 \
\ »/ J \\ ' ;;Д
К АД \\
/ / \ \ % /А ч
/ /' W
/ V
/х // .......... \\ Ni^p
Литература
1. Giese A.T., French C.S. The analysis of overlapping spectral absorption bands by derivative spectrophotometry // Applied Spectroscopy. 1955. V. 9. N 2. P. 78-96. doi: 10.1366/000370255774634089
2. Краулиня Э.К., Лиепа С.Я., Пикалов В.В., Скудра А.Я. К проблеме исследования атомной сенсибилизированной флуоресценции по контурам спектральных линий // Некорректные обратные задачи атомной физики / Сб. статей под ред. Н.Г. Преображенского. Новосибирск: Изд-во ИТПМ, 1976. С. 61-72.
3. Kauppinen J.K., Moffatt D.J., Mantsch H.H., Cameron D.G. Fourier self-deconvolution: a method for resolving intrinsically overlapped bands // Applied Spectroscopy. 1981. V. 35. N 3. P. 271-276. doi: 10.1366/0003702814732634
4. Михайленко В.И., Михальчук В.В. Методы разложения спектров с неразрешенной структурой // Журнал прикладной спектроскопии. 1987. Т. 46. № 4. С. 535-543.
5. Goldberg D.E. Genetic Algorithms in Search, Optimization and Machine Learning. Addison-Wesley, 1989. 412 p.
6. Манойлов В.В., Заруцкий И.В. Возможности алгоритма сверток с производными для оценки параметров масс-спектров, содержащих наложившиеся пики // Научное приборостроение. 2009. Т. 19. № 4. С. 103-108.
7. Курчатов И. С. Компьютерный анализ колебательных спектров водно-спиртовых растворов. [Электронный ресурс]. Режим доступа: genphys.phys.msu.ru/rus/diploma/diploma2010/Kurchatov.pdf (дата обращения: 31.07.2017)
8. Yan L., Liu H., Zhong S., Fang H. Semi-blind spectral deconvolution with adaptive Tikhonov regularization // Applied Spectroscopy. 2012. V. 66. N 11. P. 1334-1346. doi:
References
1. Giese A.T., French C.S. The analysis of overlapping spectral absorption bands by derivative spectrophotometry. Applied Spectroscopy, 1955, vol. 9, no. 2, pp. 78-96. doi: 10.1366/000370255774634089
2. Kraulinya E.K., Liepa S.I., Pickalov V.V., Scudra A.I. On the Investigation of Atomic Sensibilized Fluorescence by Spectrum Line Profiles. In Ill-Posed Inverse Problems in Atomic Physics. Ed. N.G. Preobrazhensky. Novosibirsk, IPAM Publ., 1976, pp. 61-72. (In Russian)
3. Kauppinen J.K., Moffatt D.J., Mantsch H.H., Cameron D.G. Fourier self-deconvolution: a method for resolving intrinsically overlapped bands. Applied Spectroscopy, 1981, vol. 35, no. 3, pp. 271-276. doi: 10.1366/0003702814732634
4. Mikhailenko V.I., Mikhal'chuk V.V. Method of expanding spectra with unresolved structure. Journal of Applied Spectroscopy, 1987, vol. 46, no. 4, pp. 327-335. doi: 10.1007/BF00660037
5. Goldberg D.E. Genetic Algorithms in Search, Optimization and Machine Learning. Addison-Wesley, 1989, 412 p.
6. Manoylov V.V., Zarutsky I.V. Capability of the algorithm on the base convolution processing signals form for the estimation of mass-spectra peak parameters in multiplets. Scientific Instrumentation, 2009, vol. 19, no. 4, pp. 103-108. (In Russian)
7. Kurchatov I.S. Computer analysis of vibrational spectra of water-alcohol solutions. genphys.phys.msu.ru/rus/diploma/diploma2010/Kurchatov.pdf (accessed 31.07.2017)
8. Yan L., Liu H., Zhong S., Fang H. Semi-blind spectral deconvolution with adaptive Tikhonov regularization. Applied Spectroscopy, 2012, vol. 66, no. 11, pp. 1334-1346. doi:
10.1366/11-06256
9. Тихонов А.Н., Арсенин В.Я. Методы решения некорректных задач. 3-е изд. М.: Наука, 1986. 288 с.
10. Верлань А.Ф., Сизиков В.С. Интегральные уравнения: методы, алгоритмы, программы. Киев: Наук. думка, 1986. 544 с.
11. Engl H., Hanke M., Neubauer A. Regularization of Inverse Problems. Dordrecht: Kluwer, 1996. 328 p.
12. Петров Ю.П., Сизиков В.С. Корректные, некорректные и промежуточные задачи с приложениями. СПб.: Политехника, 2003. 261 с.
13. Сизиков В.С., Кривых А.В. Восстановление непрерывных спектров методом регуляризации с использованием модельных спектров // Оптика и спектроскопия. 2014. Т. 117. № 6. С. 1040-1048. doi: 10.7868/S0030403414110166
14. Sizikov V., Sidorov D. Discrete spectrum reconstruction using integral approximation algorithm // Applied Spectroscopy. 2017. V. 71. N 7. P. 1640-1651. doi: 10.1177/0003702817694181
15. Дьяконов В., Абраменкова И. MATLAB. Обработка сигналов и изображений. СПб.: Питер, 2002. 608 с.
16. Воскобойников Ю.Е., Преображенский Н.Г., Седельников А.И. Математическая обработка эксперимента в молекулярной газодинамике. Новосибирск: Наука, 1984. 240 с.
17. Сизиков В.С. Математические методы обработки результатов измерений. СПб.: Политехника, 2001. 240 с.
18. Сизиков В.С. Инфракрасная томография горячего газа: математическая модель активно-пассивной диагностики // Научно-технический вестник информационных технологий, механики и оптики. 2013. № 6(88). С. 1-17.
19. Sizikov V.S., Evseev V., Fateev A., Clausen S. Direct and inverse problems of infrared tomography // Applied Optics. 2016. V. 55. N 1. P. 208-220. doi: 10.1364/A0.55.000208
20. Nelder J.A., Mead R. A simplex method for function minimization // The Computer Journal. 1965. V. 7. N 4. P. 308313. doi: 10.1093/comjnl/7.4.308
21. Kincaid D., Cheney W. Numerical Analysis: Mathematics of Scientific Computing. 3rd ed. Providence: AMS, 2009. 788 p.
22. Holland J.H. Adaptation in Natural and Artificial Systems. Cambridge: MIT Press, 1992. 232 p.
23. Белл Р.Дж. Введение в Фурье-спектроскопию. М.: Мир, 1975. 382 с.
24. Химмельблау Д. Прикладное нелинейное программирование. М.: Мир, 1975. 536 с.
25. Дьяконов В.П. Справочник по алгоритмам и программам на языке бейсик для персональных ЭВМ. М.: Наука, 1989. 240 с.
26. Sizikov V., Evseev V. Bars and spheroids in gravimetry problem. [Электронный ресурс]. Режим доступа: https://arxiv.org/abs/1604.06927, 2016.
27. Физический энциклопедический словарь / Под ред. А.М. Прохорова. М.: Советская энциклопедия, 1984. 944 с.
28. Дьяконов В.П. MATLAB 6: Учебный курс. СПб.: Питер, 2001. 592 с.
29. Раутиан С.Г. Реальные спектральные приборы // УФН. 1958. Т. 66. № 3. С. 475-517.
Авторы
Сизиков Валерий Сергеевич - доктор технических наук, профессор, профессор, Университет ИТМО, Санкт-Петербург, 197101, Российская Федерация, sizikov2000@mail.ru Лавров Алексей Валерьевич - старший преподаватель, Университет ИТМО, Санкт-Петербург, 197101, Российская Федерация, lavrov@corp.ifmo.ru
10.1366/11-06256
9. Tikhonov A.N., Arsenin V.Ya. Solutions of Ill-Posed Problems. New York, Wiley, 1977.
10. Verlan' A.F., Sizikov V.S. Integral Equations: Methods, Algorithms, Programs. Kiev, Naukova Dumka, 1986, 544 p. (In Russian)
11. Engl H., Hanke M., Neubauer A. Regularization of Inverse Problems. Dordrecht, Kluwer, 1996, 328 p.
12. Petrov Yu.P., Sizikov V.S. Well-Posed, Ill-Posed, and Intermediate Problems with Applications. Leiden-Boston, VSP, 2005, 234 p.
13. Sizikov V.S., Krivykh A.V. Reconstruction of continuous spectra by the regularization method using model spectra. Optics and Spectroscopy, 2014, vol. 117, no. 6, pp. 10101017. doi: 10.1134/S0030400X14110162
14. Sizikov V., Sidorov D. Discrete spectrum reconstruction using integral approximation algorithm. Applied Spectroscopy, 2017, vol. 71, no. 7, pp. 1640-1651. doi: 10.1177/0003702817694181
15. D'yakonov V., Abramenkova I. MATLAB. Processing of Signals and Images. St. Petersburg, Piter Publ., 2002, 608 p. (In Russian)
16. Voskoboinikov Yu.E., Preobrazhensky N.G., Sedel'nikov A.I. Mathematical Processing of Experiment in Molecular Gas Dynamics. Novosibirsk, Nauka Publ., 1984, 240 p. (In Russian)
17. Sizikov V.S. Mathematical Methods for Processing the Results of Measurements. St. Petersburg, Polytechnika Publ., 2001, 240 p. (In Russian)
18. Sizikov V.S. Infrared tomography of hot gas: mathematical model of active-passive diagnosis. Scientific and Technical Journal of Information Technologies, Mechanics and Optics, 2013, no. 6, pp. 1-17. (In Russian)
19. Sizikov V.S., Evseev V., Fateev A., Clausen S. Direct and inverse problems of infrared tomography. Applied Optics, 2016, vol. 55, no. 1, pp. 208-220. doi: 10.1364/A0.55.000208
20. Nelder J.A., Mead R. A simplex method for function minimization. The Computer Journal, 1965, vol. 7, no. 4, pp. 308-313. doi: 10.1093/comjnl/7.4.308
21. Kincaid D., Cheney W. Numerical Analysis: Mathematics of Scientific Computing. 3rd ed. Providence, AMS, 2009, 788 p.
22. Holland J.H. Adaptation in Natural and Artificial Systems. Cambridge, MIT Press, 1992, 232 p.
23. Bell R.J. Introductory Fourier Transform Spectroscopy. NYLondon, Academic Press, 1972, 382 p.
24. Himmelblau D.M. Applied Nonlinear Programming. NY, McGray-Hill, 1972, 416 p.
25. D'yakonov V.P. Handbook on Algorithms and Programs in Basic Language for PC. Moscow, Nauka Publ., 1989, 240 p. (In Russian)
26. Sizikov V., Evseev V. Bars and Spheroids in Gravimetry Problem. https://arxiv.org/abs/1604.06927, 2016
27. Encyclopedic Dictionary of Physics. Ed. A.M. Prokhorov. Moscow, Sovetskaya Entsiklopediya Publ., 1984. 944 p.
28. D'yakonov V.P. MATLAB 6: Training Course. St. Petersburg, Piter Publ., 2001, 592 p. (In Russian)
29. Rautian S.G. Real spectral apparatus. Soviet Physics Uspekhi, 1958, vol. 66(1), no. 2, pp. 245-273. doi: 10.1070/PU1958v001n02ABEH003099
Authors
Valery S. Sizikov - D.Sc., Full Professor, ITMO University, Saint Petersburg, 197101, Russian Federation, sizikov2000@mail.ru
Aleksei V. Lavrov - senior lecturer, ITMO University, Saint Petersburg, 197101, Russian Federation, lavrov@corp.ifmo.ru
