Обработка дискретных спектров с помощью алгоритма интегральной аппроксимации Текст научной статьи по специальности «Математика»

УДК 535.33+517.615.5

ОБРАБОТКА ДИСКРЕТНЫХ СПЕКТРОВ С ПОМОЩЬЮ АЛГОРИТМА ИНТЕГРАЛЬНОЙ АППРОКСИМАЦИИ А.В. Кривых, В.С. Сизиков

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

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

Введение

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

Под спектром и(у) подразумевается зависимость интенсивности излучения и от частоты V. Спектры бывают непрерывные, дискретные, полосатые и комбинированные. Для разложения излучения в спектр и его регистрации используются оптические спектральные приборы [2, С. 703].

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

Целью данной работы является разработка методики и программного обеспечения для обработки дискретных спектров (решения обратной задачи спектроскопии) - восстановления истинного дискретного спектра по измеренному спектру и известной аппаратной функции спектрометра путем решения системы линейно-нелинейных уравнений (часть неизвестных входит в систему линейно, а часть - нелинейно) алгоритмом интегральной аппроксимации в рамках системы Ма1ЬаЪ.

Задача восстановления спектра называется обратной задачей спектроскопии [1, 3, 4, 6, 7] или задачей редукции к идеальному спектру [1, 4, 7] и является одним из вариантов редукционной проблемы Рэ-лея [4, 7].

Постановка задачи

Измеренный спектрометром (например, интерферометром Фабри-Перо [2, 5]) спектр и(у) обычно отличается от истинного спектра г(у), во-первых, большей сглаженностью (не разрешены близкие линии -результат воздействия аппаратной функции спектрального прибора К (V, V') [3-5, 7]), а, во-вторых, зашумленностью (слабые линии «тонут» в шуме - результат случайных погрешностей измерений [1, 3, 4, 7]).

Можно дать следующее определение аппаратной функции [1; 2, С. 704; 4-7]. Аппаратной функцией, АФ (спектральной чувствительностью, функцией пропускания, частотной характеристикой), спектрометра К (V, V') называется реакция спектрометра (в виде измеренной интенсивности) на дискретную линию единичной интенсивности и частоты V' при настройке спектрометра на частоту V .

Фиксируя V и изменяя V', получим некоторую зависимость К(V, V') в виде кривой (рис. 1). Аналогичные кривые получим для других значений V. В результате получим двухмерную функцию

К(V, V') .

1 уТ\

к^, v')

v' v v, v'

Рис. 1. Зависимость К (V, V') при некотором фиксированном V

Чем шире K(v, v'), тем более заглаженным будет измеренный спектр u(v) по сравнению с истинным спектром z(v). В случае дискретного спектра, когда искомый спектр z(v) состоит из отдельных почти монохроматических спектральных линий, характеризуемых их частотами и интенсивностями, задача восстановления истинного спектра описывается следующими соотношениями:

n _

^K(vI,vj) Zj + F = u(vi), i = 1,m , с <vi < d , (1)

j=1

где zj - амплитуда (интенсивность)j-й линии; vj - ее частота; n - число линий; vi - дискретный отсчет частоты настройки спектрометра v; m - число отсчетов; [с,d] - диапазон частот; u(vi) = u(vi) + 5w(vi), 5u - случайная компонента шума измерений; F - детерминированная компонента шума (фон).

В (1) известны u(vi), K(v, v' ), vi, с, d, m, а искомыми являются z j , vj , n, F. Соотношения (1)

образуют систему линейно-нелинейных уравнений (СЛНУ), поскольку часть неизвестных (z j и F) входит линейно, а часть (v j) - нелинейно.

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

Для решения системы (1) можно воспользоваться методами, предназначенными для решения СЛНУ, например, методом Прони [10], алгоритмом Пиблза-Берковича [11], алгоритмом Фальковича-Коновалова [12], но они либо ориентированы на специальный тип СЛНУ, либо оказываются весьма неточными, либо являются слишком громоздкими.

Можно использовать так называемый метод переменных проекций (the variable projection method) Голуба-Хегланда-Муллена [13], в котором также решается СЛНУ, однако для отыскания частот используется нелинейный метод (типа Гаусса-Ньютона).

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

Для эффективного решения СЛНУ (1), учитывающего ее специфику, воспользуемся алгоритмом интегральной аппроксимации [4, 7, 14-16], который учитывает особенности этой системы и который уже продемонстрировал свою эффективность в обработке сигналов [15]. Согласно данному алгоритму, реализуется следующая последовательность действий.

1. Решается интегральное уравнение (ИУ) Фредгольма I рода

b

J K (v, v') z(v') dv' = V(v), с <v< d (2)

a

методом регуляризации Тихонова [4, 17] с заниженным значением параметра регуляризации а (это необходимо для разрешения близких линий). В результате будет получено решение za (v' ), в котором могут разрешиться близкие линии, но из-за пониженности а также возникнут ложные флуктуации-линии.

2. В полученном решении za (v' ) на основе дополнительной информации выделяется ограниченное количество (L < N) наиболее мощных максимумов, причем N задается так, чтобы N > n , где n - предполагаемое число линий. Фиксируются частоты наиболее мощных максимумов v j , j = 1, L .

3. Решается уточняющая система линейных алгебраических уравнений (СЛАУ)

L ~ _

^K(vi, vj )Vj + F = V(vi), i = 1,m , с <vi < d , (3)

j=1

методом наименьших квадратов Гаусса относительно L, интенсивностей vz j и фона F .

4. Оставляются лишь те линии, значения интенсивностей vz j которых преодолели некоторый априори

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

Достоинством алгоритма является то, что наиболее сложная часть задачи - определение значений нелинейно входящих параметров (частот спектральных линий vv j ) - решается линейно, а именно, путем

решения линейного ИУ (2).

Метод регуляризации Тихонова

Задача решения уравнения (2) является некорректной [3, 4, 7, 17] (если решать уравнение (2), например, методом квадратур, то в качестве решения получим так называемую «пилу» [4, С. 182; 7, С. 205] - крайне неустойчивое решение). Исходя из этого, для его устойчивого решения необходимо применение устойчивых методов, например, метода регуляризации Тихонова [3, 4, 7, 9, 17]. Применительно к интегральному уравнению Фредгольма I рода (типа (2)) ь

Ау = |К(х,5)у(5)= /(х), с < х < ё ,

а

метод регуляризации Тихонова сводится к решению интегрального уравнения

ь

ауа (Г) +1Я(Г, 5) уа (5) = ¥(/), а < I < Ь ,

а

где а > 0 - параметр регуляризации, а новое ядро и новая правая часть равны

Я(/,5) = Я(5,0 = |К(х,0К(х,5)сЬс, ¥(0 = |К(х,0/(х)ёх .

с с

Численный пример

В рамках системы Ма1ЬаЪ7 разработано программное обеспечение для восстановления дискретных спектров, реализующее алгоритм интегральной аппроксимации, а также решен модельный пример [4, С. 89; 7, С. 220].

В нем истинный спектр задавался в виде семи дискретных спектральных линий с амплитудами (в условных единицах) 21 = 4,4, 12 = 4,6 , 23 = 1,1, 24 = 3,2, 25 = 3,2 , = 2,8, 27 = 3,6 и частотами (также в условных единицах) у1 = 2,28, у'2 = 2,36, у3 = 2,95 , у'4 = 3,02 , у'5 = 3,56, у6 = 3,64 , у7 = 3,69 .

АФ спектрометра задавалась частотно-неинвариантной (ширина К уменьшается с увеличением частоты настройки спектрометра V) функцией

к(у, V') = 0,9 -ехр

£ е

( (У-У')2 2ст 2(1 - 0,16у)

где ст = 0,05 . Значение детерминированной компоненты шума (фона) было взято ¥ = 0,2, а случайная компонента шума измерений имела среднеквадратическое отклонение (СКО) [18], равное 0,05 (2%).

На рис. 2 представлены истинный дискретный (линейчатый) спектр х(у), состоящий из 7 линий, измеренный (экспериментальный) спектр без шума и(у) и с шумом и (у) , а также АФ спектрометра К (у, у') на низкой и высокой частотах. Видим, что истинный спектр г(у) содержит близкие линии (две

слева, две посередине и три справа), которые в измеренном спектре и(у) не разрешаются.

Применив алгоритм интегральной аппроксимации для восстановления истинного спектра, сначала решаем ИУ (2) методом регуляризации Тихонова при а = 10-6. Полученное регуляризованное решение 2а (у) приведено на рис. 3. В нем разрешились все истинные спектральные линии, однако появилось много ложных линий (максимумов).

Взяв в регуляризованном решении 2а (у) первые Ь = 12 наиболее мощных максимумов, решаем уточняющую СЛАУ (3) относительно Ь +1 = 13 неизвестных (12 амплитуд ~ и фона ¥) методом наименьших квадратов Гаусса. На рис. 3 отмечены полученные значения ~ . Видим, что все ложные максимумы получили отрицательные значения или значения, близкие к нулю (в качестве барьера для отфильтро-вывания ложных линий использовался порог 2, равный 20% от полученного значения фона ¥, т.е. 2 = 0,2¥), а истинные максимумы получили значения ~ , весьма близкие к точным значениям амплитуд

В результате можно констатировать, что в модельном примере все 7 спектральных линий разрешились и с приемлемой точностью определились их частоты V- и интенсивности ~ , причем ни одна линия

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

I 2

<и

2

1 Г

| 3

1

1

5

4 \ / 1 \ /\

у

0

2

2,2

2,4

2,6

3,2

3,4

3,6

4

2,8 3 Частота, у. е.

Рис. 2. Численный пример. Прямая задача. По оси абсцисс - частота V , по оси ординат - г, и, и и К (в условных единицах): 1 - истинный дискретный спектр 2^); 2 - заглаженный спектр и (у); 3 - заглаженный и зашумленный спектр м(у); 4 - К(2,1, V ') - АФ на низкой частоте; 5 - К(3,9, V ') - АФ на высокой частоте

5

4

3

1

<и

о

0

1

со ^

0

1

ф

Л /' ! ! V ! \ !

н-

Л /. * \ \;

_!_к_'.__

-1

2 2,2 2,4 2,6 2,8 3 3,2 3,4 3,6 3,8 4

Частота, у.е.

Рис. 3. Численный пример. Обратная задача. По оси абсцисс - частота V , по оси ординат - г (в условных единицах): 1 - истинный спектр 2^) (вертикальные сплошные линии); 2 - регуляризованное решение 2а (V) (пунктир); 3 - восстановленный спектр (V) (вертикальные штрих-пунктирные линии)

Заключение

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

2

5

3

4

3

2

1

1

0

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

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

Работа выполнена при финансовой поддержке РФФИ (грант № 09-08-00034).

Литература

1. Раутиан С.Г. Реальные спектральные приборы // Успехи физических наук. - 1958. - Т. 66. - Вып. 3. -С. 475-517.

2. Физический энциклопедический словарь / Гл. ред. А.М. Прохоров. - М.: Сов. энциклопедия, 1984. - 944 с.

3. Кочиков И.В., Курамшина Г.М., Пентин Ю.А., Ягола А.Г. Обратные задачи колебательной спектроскопии. - М.: Изд-во МГУ, 1993. - 204 с.

4. Сизиков В.С. Математические методы обработки результатов измерений. - СПб: Политехника, 2001. - 240 с.

5. Ландсберг Г.С. Оптика: Учебное пособие для вузов. - 6-е изд. - М.: Физматлит, 2006. - 848 с.

6. Сизиков В.С., Кривых А.В. Использование способа моделирования при решении обратной задачи спектроскопии методом регуляризации // Изв. вузов. Приборостроение. - 2011. - Т. 54. - № 9. - С. 44-51.

7. Сизиков В.С. Обратные прикладные задачи и MatLab. - СПб: Лань, 2011. - 247 с.

8. Химельблау Д. Прикладное нелинейное программирование. - М.: Мир, 1975. - 536 с.

9. Гончарский А.В., Черепащук А.М., Ягола А.Г. Численные методы решения обратных задач астрофизики. - М.: Наука, 1978. - 336 с.

10. Кей С.М., Марпл С.Л. Современные методы спектрального анализа (обзор) // Труды Ин-та инж. по электротехнике и радиоэлектрон. - 1981. - Т. 69. - № 11. - С. 5-51.

11. Пиблз, Беркович. Многолучевой моноимпульсный радиолокатор // Зарубежн. радиоэлектроника. -1969. - № 10.

12. Фалькович С.Е., Коновалов Л.Н. Разрешение неизвестного числа сигналов // Радиотехника и электроника. - 1982. - Т. 27. - № 1. - С. 92-97.

13. Mullen K.M., van Stokkum I.H.M. The variable projection algorithm in time-resolved spectroscopy, microscopy and mass spectrometry applications // Numerical Algorithms. - 2009. - V. 51. - № 3. - P. 319-340.

14. Сизиков В.С. О моделировании некоторых некорректных задач с использованием принципов подобия // Электрон. моделирование. - 1981. - № 6. - С. 3-8.

15. Сизиков В.С. Обобщенный метод редукции измерений. I, III // Электрон. моделирование. - 1991. -Т. 13. - № 4. - С. 7-14; № 6. - С. 3-9.

16. Верлань А.Ф., Сизиков В.С., Мосенцова Л.В. Метод вычислительных экспериментов для решения интегральных уравнений в обратной задаче спектроскопии // Электрон. моделирование. - 2011. -Т. 33. - № 2. - С. 3-12.

17. Верлань А.Ф., Сизиков В.С. Интегральные уравнения: методы, алгоритмы, программы. - Киев: Нау-кова думка, 1986. - 544 с.

18. Дайнеко М.В., Сизиков В.С. Восстановление смазанных под углом и зашумленных изображений без учета граничных условий // Научно-технический вестник СПбГУ ИТМО. - 2010. - № 4 (68). - С. 28-32.

Кривых Александр Владимирович Сизиков Валерий Сергеевич

Санкт-Петербургский государственный университет информационных технологий, механики и оптики, студент, krivykh1987@mail.ru Санкт-Петербургский государственный университет информационных технологий, механики и оптики, доктор технических наук, профессор, sizikov2000@mail.ru