Способ деконволюции спектрометрической информации и обнаружения спектральных пиков Текст научной статьи по специальности «Электротехника, электронная техника, информационные технологии»

ISSN 0868-5886

НАУЧНОЕ ПРИБОРОСТРОЕНИЕ, 2017, том 27, № 2, c. 75-82

МАТЕМАТИЧЕСКИЕ МЕТОДЫ г

И МОДЕЛИРОВАНИЕ В ПРИБОРОСТРОЕНИИ

УДК 621.391.837: 681.3 © Б. В. Бардин

СПОСОБ ДЕКОНВОЛЮЦИИ СПЕКТРОМЕТРИЧЕСКОЙ ИНФОРМАЦИИ И ОБНАРУЖЕНИЯ СПЕКТРАЛЬНЫХ ПИКОВ

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

Кл. сл.: деконволюция спектра, обнаружение спектральных пиков, фильтр Винера

ВВЕДЕНИЕ

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

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

b

J h( x, y)s( y)dy + n( x) = u( x), (i)

a ^ '

с < x < d,

в котором h(x, y) — аппаратная функция, s(x) — искомая функция (сигнал), u(x) — зарегистрированный прибором спектр, n(x) — случайный шум измерения.

Если параметры прибора не зависят от времени, то аппаратная функция зависит только от разности аргументов

h(x, y) = h(x - y), (2)

и интеграл в уравнении (1) становится интегралом свертки или конволюции (convolution). Поэтому задачу восстановления сигнала принято называть деконволюцией.

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

Задачей расшифровки спектра, как правило, является определение двух параметров для каждого спектрального пика — положение на шкале и величина (амплитуда или площадь) пика. Такую информацию можно представить как набор дельта-функций с соответствующими параметрами, который можно трактовать как исходный спектр s(x) в выбранной математической модели. При этом аппаратная функция h(x, y) в уравнении (1) будет характеризовать отклонение зарегистрированного спектра от такого "истинного" спектра, определяемое как параметрами прибора, так и физическими свойствами процесса.

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

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

СОСТОЯНИЕ ПРОБЛЕМЫ ДЕКОНВОЛЮЦИИ

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

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

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

Для задачи восстановления сигналов в следующей постановке:

1) сигнал представляет стационарный процесс,

2) в качестве критерия оптимизации используется минимум среднего квадрата ошибки,

3) используется линейная обработка информации, —

оптимальным решением является фильтр Винера [1]

Ж (ю) = -

Н *(ю)

Н (ю)| +

Ры (ю) Ps (ю)

(3)

где Н(ю) Н*(ю) -

фурье-образ аппаратной функции сопряженная величина, а Р5! (ю)

и Ры (ю) — спектральные плотности энергии соответственно сигнала и шума. Аппаратная функция в рассматриваемой задаче соответствует по форме и ширине спектральному пику.

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

- фильтр Тихонова: не требует информации о спектральной плотности энергии сигнала [2];

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

плотности энергии шума [3];

- свертка сигнала со второй производной аппаратной функции: не требует информации о спектральной плотности энергии как сигнала, так и шума [4].

Таким образом, фильтр Винера реализует теоретический предел разрешения в поставленной задаче. Снятие последнего (3-го) ограничения в постановке задачи открывает широкое многообразие нелинейных методов обработки сигналов. При этом естественно возникает вопрос о теоретическом пределе разрешения в такой постановке задачи. Такой предел был найден Е.Л. Косаревым [5]. В этой же работе им было определено, что этот предел может быть реализован методом максимума правдоподобия. Естественно, что разрешение в данном случае получается существенно выше, чем при использовании фильтра Винера. Однако решение задачи в этом случае требует больших временных затрат, т. к. в процессе решения могут потребоваться тысячи и десятки тысяч итераций алгоритма, на каждой из которых выполняется большой объем вычислений.

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

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

ДЕКОНВОЛЮЦИЯ И ПОВЫШЕНИЕ РАЗРЕШЕНИЯ

Основой рассматриваемого способа является фильтр Винера (3). Основной трудностью применения фильтра Винера является необходимость знания спектральной плотности энергии сигнала и шума. Однако, как отмечено выше, в задачах спектрометрии искомым входным сигналом является последовательность дельта-функций, или импульсный процесс. Спектральная плотность такого процесса [6]

Р3 (ю) = 2va2,

(4)

т. е. является константой.

Здесь V — средняя частота событий (импульсов), а2 — значение среднего квадрата амплитуд импульсов, и эти характеристики обычно можно определить из входных экспериментальных данных.

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

¿(ю) = РМ/ Ps (ю), (5)

становится константой.

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

В противном случае, для определения точного выражения для фильтра Винера необходимо производить предварительные сложные исследования характеристик шумов для разных типов задач. Однако на практике чаще всего на основе предыду-

щего опыта в качестве Х(ш) подбирают некоторую константу. Типичные значения такой константы обычно лежат в пределах 10-10-8.

Реакция фильтра Винера (3)

^ «) = ^4 [ Ж (ю) • ^ [ (t) + и(t)]] (6)

на модельный спектр отображена на рис. 1, а, б.

Здесь: £г0^) — зарегистрированный спектр при отсутствии шумов, п(() — шум, и — соответственно прямое и обратное пре-

образования Фурье.

Модельный спектр на рис. 1 состоит из двух пачек — слившихся гауссовых пиков. Первая (левая) пачка состоит из двух пиков шириной 100 ед., амплитудой 1000 ед., и с расстоянием между пиками 80 ед. Вторая пачка состоит из трех пиков. Первые два пика имеют ширину 100 ед., амплитуду 1000 ед. и расстояние между этими пиками 65 ед. Третий пик имеет ширину 25 ед., амплитуду 250 ед., а расстояние между вторым и третьим пиками составляет 75 ед.

Рис. 1. Деконволюция спектрометрической информации.

Sг — входной регистрируемый сигнал; Sv — реакция фильтра Винера

Аппаратная функция принята соответствующей гауссовому пику шириной 100 ед. Параметр регуляризации (5) фильтра Винера принят равным Л(о) = 2.5 • 10~4. На рис. 1, а, шум отсутствует — п(() = 0. На остальных рисунках добавлен белый гауссов шум. Отношение сигнал/шум = 10 на уровне 2 с.к.о.

Реакция фильтра Винера на рис. 1, а, за пределами полезного сигнала — пика или пачки пиков, содержит колебания Гиббса, которые в сложных спектрах могут создавать проблемы обнаружения пиков — приводить к нахождению ложных пиков.

На рис. 1, б, к реакции фильтра добавляются шумы из той части спектра входных шумов, которые проходят через фильтр, т. е. совпадают по полосе частот с полезным сигналом и поэтому линейной фильтрацией их не устранить. Это дополнительно усугубляет проблему ложных пиков.

Для подавления нежелательных колебаний и шумов там, где их не должно быть — в отсутствие сигнала, в настоящей работе использована дополнительная операция умножения выхода фильтра Винера на входной сигнал

5v(t) = ^1 [Ж(о) • ^[£г0^) + п^)]] • ^с^) / Ап . (7)

Здесь Ап — нормирующий делитель сохранения размерности, в качестве его может быть использовано максимальное или среднее значение сигнала.

Поскольку на практике сигнал 5г0^) неизвестен, то для его оценки можно пропустить непосредственно зарегистрированный спектр ) = 5г0^) + п^) через сглаживающий фильтр с целью устранения высокочастотных шумов.

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

Реакция теперь уже нелинейного фильтра (7) отображена на рис. 1, в, из которого видно, что в данной модели колебания Гиббса и низкочастотные шумы подавлены практически полностью.

На рис. 1, г, показан выход операции повышения разрешения (7), когда в модель добавлена погрешность выполнения операции устранения фона. В этом случае колебания и шумы не подавляются полностью, что и свидетельствует о важности последней операции. Причем погрешность может быть как положительной (как на рисунке) — недокомпенсация фона, так и отрицательной — перекомпенсация.

ОБНАРУЖЕНИЕ СПЕКТРАЛЬНЫХ ПИКОВ

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

Рис. 2. Обнаружение спектральных пиков.

— абсолютная величина реакции фильтра Винера 5,; Т^ — фоновая составляющая пороговой кривой (8); Thl — "гауссовская" локальная составляющая (9) пороговой кривой; — результат пороговой обработки сигнала 15,1; — результат применения к операции двойного дифференцировая

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

На рис. 2, а, кривая это абсолютная величина £ на рис. 1, г.

Если средняя фоновая длина, занимаемая пиками на оси абсцисс графика меньше половины общей длины спектра, то фоновый уровень шумов может быть вычислен медианным фильтром с шириной окна, в два раза большей длины максимальной пачки пиков. Тогда фоновая составляющая пороговой кривой Т1 будет

т ¿Г) = med_Flltr (| 5У|, Ж )• К г

(8)

где — окно фильтра, Кг — коэффициент

величиной порядка 3-5.

Распространенной в спектрометрии помехой или типом шума являются ложные, сателлитные пики, которые могут возникать около "истинных" пиков. Так, например, в капиллярном электрофорезе, в секвенаторе ДНК-анализа такие ложные пики появляются при недостаточном подавлении перекрестного влияния частотных каналов флуоресценции. При обработке сателлитные пики могут появляться при недостаточном подавлении колебаний Гиббса. Существуют также другие специфические причины появления таких помех. В рассматриваемой модели (рис. 1) сателлитный пик имитируется последним — третьим пиком второй пачки.

Для подавления сателлитных пиков в пороговую кривую можно ввести локальную составляющую, которая вычисляется как реакция на сигнал

1^1 гауссового сглаживающего фильтра с шириной ядра, равной, например, удвоенной ширине искомых "истинных" пиков, соответствующих аппаратной функции, т. е. локальная составляющая

Щ(г) = gauss_Flltr(|5У|, Ж ^)•

К ь

(9)

где Ж ^ — ширина ядра фильтра, Кь — коэффициент порядка 0.5-2.0.

В целом пороговая кривая (рис. 2, а) определяется как

) =

\ТИ{^) > Th1(t): Т^Ц), [Thf(t) <ЩЦ): ).

(10)

Результат пороговой обработки (рис. 2, б) сигнала к „(¿)| определяется как

£ vh(t) =■

к v(t)| - т) > 0:

кv(t) - ) < 0:

£v(t)| - ТЪ^), 0.

(11)

Как видно из рис. 2, б, за счет введения локального порога (Кь= 1.2) сателлитный пик подавлен, и, кроме того, полностью разрешены пики первой пачки.

Сигнал £практически полностью очищен от шумов — высокочастотные шумы подавлены фильтром Винера, низкочастотные шумы отсечены пороговой обработкой. Это позволяет дополнительно повысить разрешение задачи обнаружения пиков при помощи простой, достаточно эффективной, но чувствительной к шумам операции — двойного дифференцирования. На рис. 1, в, представлена инвертированная 2-я производная £ла(0 от сигнала £А(0. При этом, поскольку шумы отсутствуют, 2-я производная вычислялась как 2-я разность. Как видно из рисунка, все полезные пики рассматриваемого модельного спектра разрешены.

Положение пика на шкале определяется как положение максимального значения сигнала £ла(0 в соответствующей области, где £Аа(0 превышает нуль на произвольно малый порог.

РЕЗУЛЬТАТЫ И ВЫВОДЫ

На рис. 3 в качестве примера показано обнаружение пиков для участка спектра капиллярного электрофореза для задачи ДНК-анализа методом секвенирования ДНК. На рисунке Рк — обнаруженные пики, причем амплитуда и ширина пиков показаны условно, т. к. в настоящей работе рассматривается только обнаружение пиков и определение их положения на шкале.

Рассмотрим вопрос предельного разрешения рассматриваемо метода и его сравнения с другими методами.

Рис. 3. Обнаружение пиков на участке спектра капиллярного электрофореза для задачи ДНК-анализа. £г — входной регистрируемый сигнал; Рк — спектральные пики сигнала

Разрешение удобно определить как

Rs = 5/W, (12)

где 5 — расстояние между соседними пиками, W— ширина пиков. Т. е. это минимально различимое относительное расстояние между пиками, и чем меньше Rs, тем выше разрешение.

Классическое определение разрешения, восходящее к Рэлею, не предполагает решения интегрального уравнения (1), т. е. не использует процедуру повышения разрешения, и за предельное разрешение принимается 5 = W, или Rs = 1. Но это фактически декларация, предназначенная для "ручного" восприятия информации. Например, в случае оптики зто соответствует визуальному различению функций s(j).

Для линейных методов деконволюции предельное разрешение можно приблизительно оценить по рис. 1, а, б, для первой пачки пиков величиной порядка 0.8. И хотя наложившиеся пики этой пачки достаточно хорошо разделились, появление колебаний Гиббса и усиление низкочастотных шумов в реальных сложных спектрах со значительным диапазоном амплитуд в пределах пачки приводят к большим сложностям по определению пороговой линии обнаружения пиков. Не самый сложный пример такого спектра приведен на рис. 3.

Степень повышения разрешения любым методом ограничивается уровнем шумов и зависит от формы аппаратной функции. В работе [5] при помощи модельного эксперимента методом максимального правдоподобия определен теоретический предел разрешения для гауссовской аппаратной функции и отношения сигнал/шум 45 дБ. В этой работе для характеристики повышения разрешения используется величина сверхразрешения sR, которая является обратной величиной Яз с некоторым коэффициентом. Перерасчет предельной величины sR дает Яз = 0.15.

Для сравнения рассматриваемого в настоящей работе метода с методом максимума правдоподобия был проведен модельный эксперимент по разделению и идентификации пиков в пачке из двух пиков шириной 100 ед. с таким же уровнем шума (с/ш = 45 дБ). Кроме того, поскольку заданный уровень шума низкий, с целью повышения разрешения параметр регуляризации снижен до величины Л = 2.5•Ю-8. Расстояние между пиками подобрано близкое к предельному — 35 ед.

На рис. 4 отображен процесс работы метода. Обозначения кривых на рисунке соответствуют обозначениям рис. 1.

Рис. 4. Предельный случай задачи обнаружения пиков 5г — входной регистрируемый сигнал; 5, — реакция фильтра Винера; |5,| — ее абсолютное значение; П — пороговая кривая; — выходной сигнал метода

Как видно из рис. 4, а (реакция фильтра Винера), заниженное значение параметра регуляризации X приводит к усилению (раскачиванию) низкочастотных шумов, однако они практически полностью подавляются операцией перемножения выхода фильтра на входной сигнал (рис. 4, б), и в результате пики разрешаются полностью (рис. 4, г).

Таким образом, разрешение рассматриваемого способа Rs = 0.35 в 2.3 раза хуже теоретического предела, однако время решения задачи этим способом на несколько десятичных порядков меньше времени решения методом максимума правдоподобия. Так, время решения задачи, представленной на рис. 1, 2, на современном, средней "мощности" персональном компьютере (процессор i3) составило 24 мс. Решение же задачи методом максимума правдоподобия, вероятно, займет время порядка минут или даже десятков минут. Во всяком случае эксперимент, описанный в работе [5], на среднем компьютере 80-х годов прошлого века длился сутки.

СПИСОК ЛИТЕРАТУРЫ

1. Wiener N. Extrapolation, interpolation and smoothing of stationary time series with engineering applications. N.Y., J. Wiley, 1950.

2. Тихонов А.Н., Арсенин В.Я. Методы решения некорректных задач. М.: Наука, 1986. 288 с.

3. Василенко Г.И. Теория восстановления сигналов. М.: Советское радио, 1979. 272 с.

4. Сирвидас С.И., Заруцкий И.В., Ларионов А.М., Маной-лов В.В. Обнаружение, разделение и оценка параметров масс-спектрометрических пиков методом свертки экспериментальных данных с производными гауссовых функций // Научное приборостроение. 1999. Т. 9, № 2. С. 71-75.

5. Косарев Е.Л. О пределе сверхразрешения при восстановлении сигналов // Радиотехника и электроника. 1990. Вып. 1. С. 68-87.

6. Букингем М. Шумы в электронных приборах и системах. М.: Мир, 1986. 398 с.

Институт аналитического приборостроения РАН, г. Санкт-Петербург

Контакты: Бардин Борис Васильевич, bardin.bv@mai1.ru

Материал поступил в редакцию: 20.03.2017

ISSN 0868-5886

NAUCHNOE PRIBOROSTROENIE, 2017, Vol. 27, No.2, pp. 75-82

WAY DECONVOLUTION SPECTROMETER INFORMATION AND DETECTION OF SPECTRAL PEAKS

B. V. Bardin

Institute for Analytical Instrumentation of RAS, Saint-Petersburg, Russia

Basis of operation of a deconvolution of a spectrum is filtering according to Winer. Multiplication of an exit of the filter of Winer by an input spectrum is made for suppression of fluctuations of Gibbs and low-frequency noise. The threshold curve determining operation of detection of peaks contains the background component dividing a useful signal and noise, and the local component separating noise satellite peaks. 2-fold differentiation of result of threshold processing is made for additional increase in resolution of detection of peaks.

Keywords: spectrum dekonvolyution, detection of spectral peaks, filter of Wiener

REFERENСES

1. Wiener N. Extrapolation, interpolation and smoothing of stationary time series with engineering applications. N.Y., J. Wiley, 1950.

2. Tihonov A.N., Arsenin V.Ya. Metody resheniya nekor-rektnyh zadach [Methods of the solution of incorrect tasks]. Moscow, Nauka Publ., 1986. 288 p. (In Russ.).

3. Vasilenko G.I. Teoriya vosstanovleniya signalov [Theory of restoration of signals]. Moscow, Sovetskoe radio Publ., 1979. 272 p. (In Russ.).

4. Sirvidas S.I., Zaruckij I.V., Larionov A.M., Manoylov V.V. [Detection, division and assessment of parameters of mass

Contacts: Bardin Boris Vasil'evich, bardin.bv@mail.ru

and spectrometer peaks by method of convolution of experimental data with derivatives of Gaussian functions]. Nauchnoe Priborostroenie [Scientific Instrumentation], 1999, vol. 9, no. 2, pp. 71-75. (In Russ.).

5. Kosarev E.L. [About a superpermission limit at restoration of signals]. Radiotekhnika i elektronika [Soviet journal of communications technology and electronics], 1990, no. 1, pp. 68-87. (In Russ.).

6. Bukingem M. Shumy v elektronnych priborach i siste-mach [Noise in electronic devices and systems]. Moscow, Mir Publ., 1986. 398 p.

Article received in edition: 20.03.2017