eng
Выпуск #5/2017
Ю.А.Каламбет, Ю.П.Козьмин, А.C.Самохин
Фильтрация шумов. Сравнительный анализ методов
Фильтрация шумов. Сравнительный анализ методов
Просмотры: 6241
Предложена технология сравнения и отбора методов сглаживания в применении к определенным задачам. Проведена сравнительная оценка медианного сглаживания, скользящего среднего, сглаживания гауссианой, метода Савицкого-Голея и адаптивного сглаживания. Обсуждаются критерии выбора метода сглаживания для фильтрации шумов на хроматограммах.
УДК 543.08, 543.54; ВАК 02.00.02, 01.04.01
DOI: 10.22184/2227-572X.2017.36.5.88.101
УДК 543.08, 543.54; ВАК 02.00.02, 01.04.01
DOI: 10.22184/2227-572X.2017.36.5.88.101
ВВЕДЕНИЕ
Каждое измерение – это сумма полезного сигнала, случайных шумов и систематических погрешностей измерения. Шумы возникают в электронных системах регистрации, при изменении условий окружающей среды и т.д. Основная задача при сглаживании цифрового сигнала – максимально точная оценка полезного сигнала на основании массива измеренных данных. Существует мнение, что "сглаживание имеет лишь косметический эффект. Сглаженные данные не несут новой дополнительной информации" [1]. В свое время мы тоже внесли вклад в возможность расчета параметров пиков без фильтрации шумов [2]. Обсудим необходимость сглаживания позже, а прежде познакомимся с основными методами фильтрации, подробное обсуждение которых проведено в работах А.Фелингера [3] и Р.Г.Бреретона [4].
МЕТОДЫ СГЛАЖИВАНИЯ
Медианное сглаживание
В этом методе рассчитывают медианное среднее по 2k+1 точкам, и оно считается сглаженным значением в центральной k+1-й точке. Напомним, что медианное среднее находится в середине окна / интервала на k+1-й позиции (начиная с 1-й) в отсортированном массиве данных. Порядок – по возрастанию или убыванию – роли не играет, поскольку средний элемент в обоих случаях будет один и тот же. Этот метод очень эффективно убирает выбросы в данных (или, что то же самое, чрезвычайно устойчив), но потенциально меняет площадь и высоту пика. Тем не менее, медианная фильтрация по минимальному окну из трех точек в сочетании с некоторыми другими методами может давать неплохой результат.
Скользящее среднее
Самый простой метод сглаживания – скользящее среднее. В этом способе в окне из 2k+1 точек вычисляется среднее, которое считается сглаженным значением в центральной точке k+1. Таким образом, вес каждой точки равен 1/(2k+1) (рис.1). Частный случай скользящего среднего – сжатие данных (boxcar filtering), когда вместо нескольких последовательных точек берут их среднее значение. Такая процедура полезна, если исходно выбрана избыточная частота сбора данных, и ширина пика на половине высоты равна, к примеру, 200 точкам. Избыточная частота не улучшает точность анализа, а время обработки данных и требования к используемой памяти возрастают. При сжатии данных в N раз постоянная времени увеличивается в N, а величина шума уменьшается приблизительно в раз.
Вместо простого среднего можно использовать взвешенное среднее. Чаще всего, веса взвешенного среднего берутся распределенными по закону Гаусса (рис.1), что, как будет показано ниже, имеет большой смысл. Для того чтобы площадь пика не изменялась, сумма весов скользящего взвешенного среднего должна быть равна единице.
Метод Савицкого-Голея
Савицкий и Голей (Голай) [5] предложили способ фильтрации шумов, основанный на методе наименьших квадратов (МНК), который, несмотря на сложность модели, достаточно прост в расчетах. Суть заключается в построении по 2k+1 последовательным равноотстоящим точкам аппроксимирующего полинома n-й степени и использовании в качестве сглаженного значения величины полинома в k+1-й точке. Оказалось, что математически это значение вычисляется путем скользящего взвешенного среднего с весами точек, положительными в центре окна фильтрации и отрицательными на периферии (рис.1). Результаты фильтрации по Савицкому-Голею совпадают для полиномов степени 2n и 2n+1. К примеру, скользящее среднее и аппроксимация прямой (соответствующие полиному нулевой и первой степени) в качестве ответа дадут среднее значение сигнала – прямая всегда проходит через центр масс, равный среднему арифметическому, которое и выдается в качестве сглаженного значения аппроксимации прямой.
Линейные методы
Все методы, основанные на скользящем взвешенном среднем, являются линейными методами сглаживания. Это название происходит от математического термина "линейное преобразование" (в одномерном случае означающего примерно то же, что и "взвешенное среднее"), а не от аппроксимации функции линейной зависимостью. Более строго – линейным является метод, для которого результат сглаживания суммы сигналов совпадает с суммой результатов сглаживания сигналов: S(a+b)=S(a)+S(b), а также: S(k · a)=k · S(a). Не соответствующие этим условиям методы в литературе называются нелинейными, и важнейшим их представителем является медианное сглаживание.
Фильтр Калмана
Это онлайн-фильтр [6], то есть он ориентирован на фильтрацию значения в последней (по времени) из измеренных точек и использует при формировании нового значения предыдущий результат фильтрации (последнее свойство называется рекурсивностью). В простейшей реализации отфильтрованное значение складывается из взвешенной суммы последнего измеренного значения и предыдущего сглаженного. При не зависящих от времени постоянных коэффициентах этот метод носит название экспоненциального сглаживания. В общем случае фильтр строит некую самонастраивающуюся модель смены состояний системы и на ее основании производит сглаживание. Фильтр Калмана несимметричный, то есть использует только точки, измеренные ранее анализируемой, игнорируя более поздние. Этот факт можно рассматривать как достоинство (оперативное онлайн-сглаживание), а можно и как недостаток (используется только половина информации, доступной для оценки сглаженного значения). Фильтр может быть линейным или нелинейным.
Фильтрация с использованием
преобразования Фурье
Эти линейные фильтры имеют широкое распространение, но следует иметь в виду одно главное свойство большинства фурье-фильтров: они имеют аналог, который выглядит как скользящее взвешенное среднее. Причина в том, что произведение функций в фурье-пространстве превращается в свертку функций (эквивалентную скользящему взвешенному среднему) в исходном пространстве. Неприятное свойство фурье-преобразования состоит в резком уширении фурье-спектра при появлении таких особенностей сигнала, как скачок базовой линии или "срезанная" зашкаливанием детектора вершина пика. Это сильно затрудняет конструирование фильтров, базирующихся на фурье-спектре.
Фильтры, основанные на вейвлет-преобразовании
Вейвлет-преобразование по своим принципам очень похоже на преобразование Фурье, только в качестве объекта берется не вся хроматограмма, а только область пика, а в качестве базисных используются (в случае вейвлетов Эрмита) функция Гаусса и ее производные, а не синусы и косинусы. Такое разложение весьма подходит для отдельно стоящих пиков и не особенно хорошо сглаживает перекрывающиеся пики и базовую линию. Для применения вейвлет-преобразования чрезвычайно полезно заранее знать положение базовой линии и ожидаемую ширину пика.
Адаптивный метод фильтрации шумов
Несколько лет назад мы предложили адаптивный метод фильтрации шумов [7], работающий на принципах локальной минимизации доверительного интервала аппроксимации. Фильтр использует идею, впервые выдвинутую в работе [8]. Для каждой точки хроматограммы подбирается аппроксимационная модель, обеспечивающая минимальный доверительный интервал. Набор функций, по которым строится аппроксимация, может быть выбран исследователем в зависимости от специфики задачи. Мы демонстрировали работоспособность алгоритма на примере полиномиальной аппроксимации.
Для применения адаптивной фильтрации требуется оценка погрешности измерений. В тех случаях, когда она невозможна или неправильна (например, детектор выдает существенно сглаженный сигнал), метод работает плохо. В качестве результата сглаживания в каждой точке выдается два числа: сглаженное значение и его доверительный интервал.
Адаптивное сглаживание нелинейно, так как оптимальная аппроксимация каждой точки зависит от сигнала. Если рассматривать сигнал с шумом как сумму идеального сигнала и шума, то при их раздельной фильтрации сигнал не изменится, а шум может быть аппроксимирован постоянной так, что он в каждой точке уменьшится в раз. Если же сглаживать сумму сигнала и шума, то аппроксимация константой невозможна, и шумоподавление будет зависеть от локального профиля сигнала.
Адаптивный фильтр отличается от известного подхода Колмогорова-Винера [9, 10] тем, что он не пытается найти общую модель, оптимизирующую методом МНК всю экспериментальную кривую. Напротив, он ищет подходящую локальную аппроксимацию, минимизирующую погрешность в каждой точке. Общей модели при этом не существует, поскольку для каждой точки она может оказаться своей. Адаптивный фильтр не подразумевает требований непрерывности и гладкости изучаемой экспериментальной зависимости, которые являются обязательным условием для подхода Колмогорова-Винера.
Новый метод с временно закрытым алгоритмом работы следует в обозримом будущем ожидать в программе "МультиХром".
Особенности методов сглаживания
Для методов, основанных на скользящем взвешенном среднем, возможно сглаживание по четному числу точек, но его результат в большинстве случаев соответствует середине окна, то есть абсцисса при симметричном окне размером 2k попадает в середину между k-й и k+1-й точками разрядной сетки. Поскольку обычно задача фильтрации шумов не подразумевает сдвига разрядной сетки, такой вариант сглаживания крайне редко используется на практике.
Массив исходных (сырых) данных должен быть больше, чем окно фильтрации. Методы, основанные на скользящем взвешенном среднем с окном 2k+1, не могут вычислить аппроксимированное значение в первых и последних k точках массива данных. Для оценки значений в этих точках приходится привлекать дополнительные (выходящие за рамки основного метода сглаживания) соображения.
В случае многократного сглаживания скользящим взвешенным средним отметим интересный эффект, следующий из центральной предельной теоремы теории вероятностей. Если веса скользящего взвешенного среднего неотрицательны и их распределение имеет дисперсию σ2, то, вне зависимости от метода сглаживания, по мере роста числа N проходов, результат будет приближаться к однократному сглаживанию скользящим взвешенным средним с весами, распределенными по закону Гаусса, и дисперсией N ∙ σ2. К примеру, (N > 10)-кратное сглаживание скользящим средним с окном 2k+1 практически эквивалентно однократному сглаживанию гауссианой с дисперсией N(k+1)2/6.
В любом случае, сочетание нескольких проходов взвешенным средним с разным распределением весов будет эквивалентно одному проходу с весами, которые легко вычисляются по весам использованных фильтров. Достаточно сначала выполнить свертку по фильтрам, получив комбинированный, которым затем и выполнять сглаживание функции. Сочетание разных способов сглаживания может иметь смысл, если принципы их работы различны: например, сочетание медианного сглаживания и взвешенного среднего при наличии выбросов в исходных данных.
КАК СРАВНИВАТЬ МЕТОДЫ СГЛАЖИВАНИЯ
Наша предметная область – хроматография, поэтому проведем обработку фрагмента данных, содержащего отдельно стоящий пик характерной формы. Предположим, что частота оцифровки достаточна для вычисления не только площади пика [11], но и других его параметров [12]. Технология сравнения методов сглаживания линейных фильтров, использованная в работе [1] не годится, так как адаптивный фильтр не линеен.
При сравнении способов сглаживания в первую очередь нужно обращать внимание на качество сглаживания хроматограммы в области базовой линии. Чем лучше аппроксимирована базовая линия вне пика, тем более точен ее прогноз для пика. Величина площади практически не зависит от того, как сглажена область пика, поскольку все методы, основанные на скользящем взвешенном среднем, оставляют площадь неизменной. Однако, на нее существенное влияние оказывает то, как проведена базовая линия.
Программа "МультиХром" позволяет проводить обработку многоканальных данных. Мы сгенерировали несколько примеров (табл.1), каждый представляет собой 100-канальную хроматограмму, во всех каналах присутствует пик (рис.2), форма которого описывается экспоненциально модифицированной гауссианой (ЭМГ) [13, 14]:
,
где hG – высота, μG – положение вершины, σG – стандартное отклонение немодифицированной гауссианы; τ – постоянная времени модифицирующей экспоненты; erfcx() – шкалированная дополнительная функция ошибок [15].
Этот пик асимметричный: коэффициент τ/σG = 3, что соответствует коэффициенту асимметрии 2,8. На каждом канале к исходному сигналу был добавлен нормально распределенный шум (рис.2б). Шумы на разных каналах смоделированы независимо, хотя и имели одинаковые параметры генерации. Величина сигмы добавленного шума (E0 = 3 333 единиц) для двух примеров соответствует пределу определения (отношение сигнал / шум, S/N = 10), то есть без фильтрации при большем шуме пик уже будет формально непригоден для количественного расчета. Для третьего примера без фильтрации шумов пик соответствует пределу обнаружения (S/N = 3).
Первый и третий примеры демонстрируют работу алгоритмов фильтрации шумов в условиях недостаточной и избыточной частоты измерений: в первом сигма гауссианы была выбрана равной единице, в третьем сигма гауссианы была увеличена до двадцати, а отношение сигнал / шум уменьшено до трех.
Исследуемую процедуру сглаживания применяли к каждому из ста каналов хроматограммы. Из сглаженной хроматограммы поканально вычитали исходную – без шума (рис.2); в результате получали разностную 100-канальную хроматограмму. В дополнение к 100 профилям исходных данных рассчитывали еще четыре, так называемые "вычисляемые" каналы (рис.3).
1. Es – средний (по 100 каналам) канал разностной хроматограммы – средняя разностная хроматограмма. Отклонение Es от нуля характеризует систематическую погрешность сглаживания (рис.3а), индекс "s" происходит от слова systematic. Систематическая погрешность – это мера смещенности сглаженного значения. Ее нельзя исключить, так как форма сигнала заранее не известна. Она должна учитываться в суммарной погрешности оценки сглаженного значения.
,
где ytrue(i) – истинное значение интенсивности в i-й точке незашумленной хроматограммы; y(i, j) – значение интенсивности в i-й точке хроматограммы на j -м канале после применения процедуры сглаживания; N – общее число точек;
2. Er – стандартное отклонение интенсивности, рассчитанное из 100 сглаженных хроматограмм (СКО). Er является характеристикой случайной погрешности сглаживания (рис.3б), индекс "r" происходит от слова random:
;
3. Et – среднее стандартное отклонение сигнала от нуля, что приблизительно равно квадратному корню из суммы квадратов значений в каналах Es и Er (рис.3в), индекс "t" происходит от слова total:
;
4. R – отношение Er/Et, нормированное на 100 (процент случайной погрешности). По этому отношению можно судить о значимости систематического отклонения в окончательном результате сглаживания (рис.3г):
.
Чем меньше погрешности в первых трех каналах, тем лучше сглаживание в соответствующей точке хроматограммы. Для методов, основанных на скользящем среднем, канал Es, соответствующий бесконечно большой статистике (не 100, а бесконечное число каналов), можно получить, применив соответствующий фильтр к незашумленному пику с последующим вычитанием исходной хроматограммы.
Имея такие графики, можно проводить визуальный и количественный анализ результатов сглаживания. Для количественных оценок мы выбрали два интервала на оси абсцисс: область пика и область базовой линии (рис.2). Область пика (ОП) – тот диапазон номеров точек, для которого отклик исходного сигнала без шума превышает заданную величину (0,1% высоты). Эта область расширена таким образом, чтобы включать возмущения сигнала вне исходного пика, вызванные фильтрацией. В примере 2 мы считали областью пика диапазон от 125-й до 210-й точки. Область базовой линии (ОБ – от 300-й до 500-й точки) – достаточно протяженный интервал точек вне пика, находящийся в промежутке между концами пика и хроматограммы на расстоянии около 100 точек от каждого. В остальных примерах область пика соответствующим образом сужается (для узкого пика) или расширяется (для широкого пика).
Мы выбрали четыре ключевые характеристики метода сглаживания:
1. BNR (baseline noise reduction) – величина, равная отношению исходного стандартного отклонения шума E0 к среднеквадратичному шуму сглаженных данных в области базовой линии; NОБ – число точек в области базовой линии:
;
2. Etmax – максимальное значение среднеквадратичной погрешности в области пика:
;
3. Etav – среднеквадратичное шумоподавление в области пика; Nоп – число точек в области пика:
;
4. Rmin – минимальная величина процента случайной погрешности в области пика:
.
Особенности сглаживания по Савицкому-Голею
В методе Савицкого-Голея существует проблема – провал перед началом пика (рис.4а). Она характерна для полиномиального сглаживания начального участка крутого склона пика и заметна при визуальном исследовании сглаженной хроматограммы, если "глубина" провала превышает 3σ остаточного случайного шума. В противном случае провал маскируют "естественные" шумы базовой линии (рис.4б). Такой эффект может появиться, когда есть точки с отрицательным весом (рис.1). Для скользящего среднего и сглаживания гауссианой таких сложностей в принципе быть не может. Провал существенно затрудняет работу алгоритмов поиска границ пиков и искажает результаты количественной обработки данных. Появление провала перед пиком – только признак проблем. Если построить для этого случая график систематической погрешности (рис.5), то окажется, что максимальное искажение формы пика превышает абсолютную величину провала более чем в пять раз. Максимум амплитуды искажений всех исследованных методов почти всегда приходится на вершину пика.
Случайная погрешность метода Савицкого-Голея (как и всех других методов, основанных на скользящем среднем) для области пика и для базовой линии одинакова. Если мы представим сигнал с шумом как сумму исходного сигнала и шума, то вследствие линейности фильтра случайная погрешность сглаженного значения определяется целиком и полностью шумовой составляющей. Поскольку размеры окна сглаживания не меняются, то и шумоподавление во всех точках одинаково.
Особенности
адаптивного метода фильтрации шумов
Представленный подход к анализу ошибок аппроксимации был применен к адаптивному методу фильтрации шумов. Характерные графики вычисляемых каналов для примера 2 приведены на рис.6. В отличие от методов сглаживания, основанных на скользящем среднем, случайная погрешность перестает быть постоянной по хроматограмме; она возрастает вблизи вершины пика и на крутых склонах. Это явление вызвано автоматической "подстройкой" ширины окна аппроксимации под данные – на базовой линии окно аппроксимации широкое, поэтому случайная погрешность уменьшается, а вблизи вершины пика окно узкое, случайная погрешность увеличивается.
Высокая случайная погрешность в области пика является платой за несмещенность оценки. Адаптивный фильтр просто не может подобрать приемлемую полиномиальную аппроксимацию этой области с большим шумоподавлением. "Хорошее" отношение систематической и случайной погрешностей в области пика является следствием не только уменьшения величины систематической погрешности в адаптивном методе по сравнению с методом Савицкого-Голея, но и менее "агрессивной" фильтрации окрестностей вершины пика.
СРАВНЕНИЕ МЕТОДОВ НА ОДНОМ ГРАФИКЕ
Все методы, работающие как скользящее взвешенное среднее, содержат параметр, называемый "окно сглаживания". Понятие "окно сглаживания" для адаптивного метода не существует или сильно отличается от других; главное, от чего зависят результаты фильтрации – оценка погрешности измерений исходных данных. В качестве параметра адаптивного фильтра, определяющего подавление шумов на базовой линии, использовали максимальный разрешенный размер аппроксимационного окна.
Зависимость шумоподавления от окна соответствующего метода приведена на рис.7. Для облегчения графического сопоставления результатов фильтрации размер окна исключили из рассмотрения; вместо него по оси абсцисс рис.8 откладывали уровень шумоподавления на базовой линии (BNR). Такой вариант представления данных сокращает число параметров и упрощает процедуру сравнения разных методов.
Применимость метода сглаживания зависит от задачи, которую он выполняет. Если требуется несмещенная оценка параметров пика, то систематическая погрешность не должна изменять форму пика в большей степени, чем это делает остаточная случайная погрешность. В этом случае может быть использован простой и наглядный критерий выбора рабочей области фильтра (условия приемлемого искажения пика): процент случайной погрешности для любой точки пика не должен падать ниже 50. При таком условии систематическая погрешность не превышает стандартного отклонения случайной погрешности ни в одной точке пика; это область приемлемого (систематического) искажения пика. Соответственно, при доле случайной погрешности менее 50% искажения пика можно считать неприемлемыми. Пороговая величина может быть и другой, учитывающей особенности аналитической методики [16]; конкретное значение порога не меняет общего подхода к задаче фильтрации шумов.
Несложно увидеть, что метод Савицкого-Голея – лучший из основанных на скользящем среднем, поскольку с его помощью удается достичь наибольшего шумоподавления при условии приемлемого искажения пика (рис.8).
Возникает вопрос: можно ли работать при больших искажениях данных, когда систематическая погрешность существенно превышает случайную? Ответ зависит от того параметра пика, который мы пытаемся оценить. Если речь идет о площади пика, то при известной базовой линии она остается почти неизменной при любом методе сглаживания, основанном на скользящем среднем. Если провал перед пиком отсутствует, погрешность площади главным образом зависит от того, насколько хорошо подавлены шумы базовой линии. Соответственно, несмотря на искажение формы пика, целевой параметр – площадь – при увеличении окна сглаживания определяется с большей точностью. Если же нужно определить высоту или ширину на половине высоты, то результат будет неудовлетворительным. Правда, если истинная высота не нужна, а параметры сглаживания одинаковы в градуировочных графиках и в анализе, то за счет уменьшения случайной погрешности высоты сглаженной хроматограммы количественная оценка концентрации по высоте может быть достаточно точной, несмотря на искажение пика. В этом смысле сглаживание гауссианой предпочтительнее метода Савицкого-Голея: ширина окна не имеет "верхнего порога", связанного с провалом перед пиком. Следует помнить, что фильтрация увеличивает ширину пиков (при сглаживании гауссианой дисперсия пика равна сумме дисперсий пика и гауссианы), что уменьшает разрешение пиков и оценку эффективности хроматографической колонки. Высота пика, сглаженного скользящим взвешенным средним, при проценте случайной погрешности менее 50 становится более воспроизводимой, но смещенной оценкой истинной величины.
Поскольку наша реализация адаптивного метода основана на полиномиальной аппроксимации, естественно сравнить его с методом Савицкого-Голея, тем более что он оказался лучшим из методов скользящего среднего (см. рис.8). В методе Савицкого-Голея неприемлемое искажение пика в примере 2 появляется при подавлении шумов на базовой линии BNR ≈ 2,5 (место пересечения кривой Rmin с горизонталью на уровне 50%). В адаптивном методе искажение такого же порядка возникает при подавлении шумов на базовой линии BNR ≈ 4,7 и незначительно изменяется при еще большем шумоподавлении. Таким образом, отношение S/N адаптивного метода превышает этот параметр оптимально настроенного метода Савицкого-Голея:
• для первого примера – в 4,4/1,4 = 3,15 раза или на 20 lg(3,15) ≈ 9,9 дБ;
• для второго примера – в 4,7/2,5 = 1,88 раза или на 20 lg(1,88) ≈ 5,5 дБ.
Для третьего примера S/N не изменяется, если допустить увеличение систематической погрешности до удвоенной остаточной случайной (процент случайной – 33), то 11/6 = 1,83 · 20 lg(1,83) ≈ 5,3 дБ; такое же значение наблюдается для нового метода без ослабления критерия.
Для некоторых примеров есть небольшая область, где метод Савицкого-Голея дает лучшую аппроксимацию области пика, чем адаптивный. Она соответствует достаточно малым (менее 2,5) коэффициентам шумоподавления на базовой линии и весьма узка.
При малой ширине пика (пример 1) происходит ожидаемое – ни один из методов обработки данных, основанных на скользящем среднем, не попадает в область приемлемого искажения пика, за исключением фильтра Савицкого-Голея по пяти точкам. Этот фильтр обеспечивает падение уровня шумов базовой линии примерно в 1,4 раза. Адаптивный фильтр продолжает выдавать приемлемое шумоподавление вплоть до падения шумов базы в пять раз. При этом реального уменьшения максимальной погрешности внутри пика нет ни у одного метода – все значения коэффициентов шумоподавления E0/ Etmax находятся ниже единицы, то есть максимальная суммарная погрешность в области пика для всех методов больше исходной. Тем не менее, применение адаптивного метода дает выигрыш в погрешности вычисления площади и высоты за счет хорошего шумоподавления на базовой линии – она оказывается проведенной более точно.
Наиболее актуальна фильтрация шумов при работе вблизи пределов определения и обнаружения [17]; именно такие примеры мы анализировали в настоящей работе (табл.2). Особенно заметно преодоление пиком (пример 3) предела определения (S/N = 10), как в случае адаптивного, так и метода Савицкого-Голея, тогда как без фильтрации шумов этот пик находился на грани предела обнаружения (S/N = 3). Результаты адаптивного фильтра и метода Савицкого-Голея для этого пика совпадают вплоть до BNR = 6, далее шумоподавление в области пика адаптивного и нового фильтров заметно лучше, чем у метода Савицкого-Голея.
Предельно допустимое шумоподавление Савицкого-Голея снижается с уменьшением ширины пика и с ростом отношения сигнал / шум. При фильтрации шумов реальной хроматограммы вполне возможна ситуация, когда в ней присутствуют пики разной ширины и высоты. В этом случае приходится выбирать параметр сглаживания таким, чтобы для самого узкого и высокого пика искажения были в рамках приемлемости. Это касается всех алгоритмов, но для адаптивного и нового фильтров можно добиться того, чтобы шумоподавление на базовой линии было существенно (в 2–3 раза) выше, чем в случае применения метода Савицкого-Голея.
Главная трудность в применении адаптивного фильтра состоит в том, что для его нормальной работы требуется знание погрешности каждого измерения. В рамках нашей модели эта погрешность считается одинаковой для всех точек хроматограммы (гомоскедастичной). Реальные погрешности для разных точек могут отличаться, тогда модель адаптивного метода необходимо приспособить к неодинаковому (гетероскедастичному) шуму [18], что заметно усложнит вычисления. Обычно считается, что погрешность измерения зависит от уровня сигнала. Поскольку основной эффект сглаживания достигается за счет фильтрации шума базовой линии, а не области пика, использование адаптивным фильтром оценки шума базовой линии, как описано в [2, 7], дает приемлемые результаты, а погрешность базовой линии вполне может считаться гомоскедастичной. Ее можно оценить по всей хроматограмме, но при условии, что ранее к данным не применяли методы сглаживания. Многие производители оборудования встраивают в детекторы собственные программы подавления шумов, существенно искажающие полученные экспериментальные результаты. Тогда оценить погрешности возможно, установив определенные настройки детектора, например, минимально возможную постоянную времени.
Формальное фармакопейное определение отношения сигнал / шум [17, 19] в случае применения адаптивного метода способно вырасти в 2–5 раз в зависимости от формы сигнала и качества базовой линии. Максимальное доступное адаптивному методу шумоподавление может оказаться ниже тех чисел, которые мы получили на графиках, поскольку оно зависит от сглаживаемых данных. Если продолжительных участков базовой линии нет, то потенциальные возможности адаптивного метода могут оказаться нереализованными.
Новый метод фильтрации шумов является достойным соперником адаптивного фильтра, поскольку, как видно из представленных в этой работе результатов (см. рис.8), по всем критериям работает устойчивее адаптивного фильтра.
ВЫВОДЫ
Предложен подход к сравнению методов сглаживания, основанный на моделировании измерительной системы и ее шумов. Модельная система позволяет построить кривые, отражающие характеристики метода сглаживания: систематическую, (остаточную) случайную и суммарную погрешности, а также долю случайной погрешности в суммарной. На основании кривых рассчитываются интегральные характеристики метода сглаживания, для хроматографического сигнала – это шумоподавление на базовой линии, минимальное и среднее шумоподавление по пику, а также минимальный процент случайной погрешности в области пика. Критерии, основанные на анализе интегральных параметров, помогают выбрать метод, наиболее подходящий для конкретной системы.
ЛИТЕРАТУРА
1. Enke C.G., Nieman T.A. Signal-to-noise ratio enhancement by least-squares polynomial smoothing // Anal. Chem. 1976. 48. P. 705A–712A. doi:10.1021/ac50002a007.
2. Каламбет Ю.А., Михайлова К.В. Оценка величины шума и ее использование при обработке хроматографического сигнала // Лабораторный Журнал. 2002. 1. С. 32–35. http://multichrom.ru/Docs/ots-vel-shuma.pdf.
3. Felinger A. Data analysis and signal processing in chromatography // Elsevier, 1998. V. 21. 413 p.
4. Brereton R.G. Chemometrics: Data Analysis for the Laboratory and Chemical Plant. Wiley, 2003. 504 p. doi:10.1198/tech.2004.s738.
5. Savitzky A., Golay M.J.E. Smoothing and Differentiation of Data by Simplified Least Squares Procedures // Anal. Chem. 1964. 36. P. 1627–1639. doi:10.1021/ac60214a047.
6. Kalman R.E. A New Approach to Linear Filtering and Prediction Problems // J. Basic Eng. 1960. 82. P. 35–45. doi:10.1115/1.3662552.
7. Каламбет Ю., Мальцев С., Козьмин Ю. Фильтрация шумов: окончательное решение проблемы // АНАЛИТИКА. 2011. 1. P. 50–55.
8. Ariyur K.B., Jelinek J. Aeroengine Prognostics via Local Linear Smoothing, Filtering and Prediction // SAE Trans. J. Aerosp. 2004. 113. P. 1773–1780. doi:10.4271/2004-01-3160.
9. Колмогоров А.Н. Интерполирование и экстраполирование стационарных случайных последовательностей // Изв. АН СССР. Сер. Матем. 1941. 5. С. 3–14.
10. Wiener N. The Extrapolation, Interpolation and Smoothing of Stationary Time Series, New York: Wiley, 1949. 176 p.
11. Kalambet Y., Samokhin A., Kozmin Y. Finally, how many ponts per peak?, in: X Winter Symp. Chemom., Samara, 2016. doi:10.13140/RG.2.2.28145.28006.
12. Wahab M.F., Dasgupta P.K., Kadjo A.F., Armstrong D.W. Sampling frequency, response times and embedded signal filtration in fast, high efficiency liquid chromatography: A tutorial // Anal. Chim. Acta. 2016. 907. P. 31–44. doi:10.1016/j.aca.2015.11.043.
13. Каламбет Ю., Михайлова К., Козьмин Ю., Нагаев И. Восстановление хроматографических пиков при помощи экспоненциально-модифицированной гауссианы // АНАЛИТИКА. 2012. № 4. С. 62–67.
14. Kalambet Y., Kozmin Y., Mikhailova K., Nagaev I., Tikhonov P. Reconstruction of chromatographic peaks using the exponentially modified Gaussian function // J. Chemom. 2011. 25. P. 352–356. doi:10.1002/cem.1343.
15. Cody W.J. Rational chebyshev approximations for the error function // Math. Comput. 1969. 23. P. 631–637. doi:10.1090/S0025-5718-1969-0247736-4.
16. Дворкин В.И. Метрология и обеспечение качества химического анализа // М: Издательство МИТХТ, 2014. 416 c.
17. Государственная фармакопея Российской Федерации. XIII издание, Федеральная электронная медицинская библиотека, 2015.
18. Каламбет Ю.А., Мальцев С.А., Козьмин Ю.П. Доверительные интервалы метода взвешенных наименьших квадратов и стратегия градуировки // Заводская Лаборатория. Диагностика Материалов. 2015. 81. С. 69–76. http://www.multichrom.ru/Docs/zav-lab15-1-1-81.pdf.
19. Ettre L.S. Nomenclature for chromatography (IUPAC Recommendations 1993) // Pure Appl. Chem. 1993. 65. P. 819–872. doi:10.1351/pac199365040819.
Каждое измерение – это сумма полезного сигнала, случайных шумов и систематических погрешностей измерения. Шумы возникают в электронных системах регистрации, при изменении условий окружающей среды и т.д. Основная задача при сглаживании цифрового сигнала – максимально точная оценка полезного сигнала на основании массива измеренных данных. Существует мнение, что "сглаживание имеет лишь косметический эффект. Сглаженные данные не несут новой дополнительной информации" [1]. В свое время мы тоже внесли вклад в возможность расчета параметров пиков без фильтрации шумов [2]. Обсудим необходимость сглаживания позже, а прежде познакомимся с основными методами фильтрации, подробное обсуждение которых проведено в работах А.Фелингера [3] и Р.Г.Бреретона [4].
МЕТОДЫ СГЛАЖИВАНИЯ
Медианное сглаживание
В этом методе рассчитывают медианное среднее по 2k+1 точкам, и оно считается сглаженным значением в центральной k+1-й точке. Напомним, что медианное среднее находится в середине окна / интервала на k+1-й позиции (начиная с 1-й) в отсортированном массиве данных. Порядок – по возрастанию или убыванию – роли не играет, поскольку средний элемент в обоих случаях будет один и тот же. Этот метод очень эффективно убирает выбросы в данных (или, что то же самое, чрезвычайно устойчив), но потенциально меняет площадь и высоту пика. Тем не менее, медианная фильтрация по минимальному окну из трех точек в сочетании с некоторыми другими методами может давать неплохой результат.
Скользящее среднее
Самый простой метод сглаживания – скользящее среднее. В этом способе в окне из 2k+1 точек вычисляется среднее, которое считается сглаженным значением в центральной точке k+1. Таким образом, вес каждой точки равен 1/(2k+1) (рис.1). Частный случай скользящего среднего – сжатие данных (boxcar filtering), когда вместо нескольких последовательных точек берут их среднее значение. Такая процедура полезна, если исходно выбрана избыточная частота сбора данных, и ширина пика на половине высоты равна, к примеру, 200 точкам. Избыточная частота не улучшает точность анализа, а время обработки данных и требования к используемой памяти возрастают. При сжатии данных в N раз постоянная времени увеличивается в N, а величина шума уменьшается приблизительно в раз.
Вместо простого среднего можно использовать взвешенное среднее. Чаще всего, веса взвешенного среднего берутся распределенными по закону Гаусса (рис.1), что, как будет показано ниже, имеет большой смысл. Для того чтобы площадь пика не изменялась, сумма весов скользящего взвешенного среднего должна быть равна единице.
Метод Савицкого-Голея
Савицкий и Голей (Голай) [5] предложили способ фильтрации шумов, основанный на методе наименьших квадратов (МНК), который, несмотря на сложность модели, достаточно прост в расчетах. Суть заключается в построении по 2k+1 последовательным равноотстоящим точкам аппроксимирующего полинома n-й степени и использовании в качестве сглаженного значения величины полинома в k+1-й точке. Оказалось, что математически это значение вычисляется путем скользящего взвешенного среднего с весами точек, положительными в центре окна фильтрации и отрицательными на периферии (рис.1). Результаты фильтрации по Савицкому-Голею совпадают для полиномов степени 2n и 2n+1. К примеру, скользящее среднее и аппроксимация прямой (соответствующие полиному нулевой и первой степени) в качестве ответа дадут среднее значение сигнала – прямая всегда проходит через центр масс, равный среднему арифметическому, которое и выдается в качестве сглаженного значения аппроксимации прямой.
Линейные методы
Все методы, основанные на скользящем взвешенном среднем, являются линейными методами сглаживания. Это название происходит от математического термина "линейное преобразование" (в одномерном случае означающего примерно то же, что и "взвешенное среднее"), а не от аппроксимации функции линейной зависимостью. Более строго – линейным является метод, для которого результат сглаживания суммы сигналов совпадает с суммой результатов сглаживания сигналов: S(a+b)=S(a)+S(b), а также: S(k · a)=k · S(a). Не соответствующие этим условиям методы в литературе называются нелинейными, и важнейшим их представителем является медианное сглаживание.
Фильтр Калмана
Это онлайн-фильтр [6], то есть он ориентирован на фильтрацию значения в последней (по времени) из измеренных точек и использует при формировании нового значения предыдущий результат фильтрации (последнее свойство называется рекурсивностью). В простейшей реализации отфильтрованное значение складывается из взвешенной суммы последнего измеренного значения и предыдущего сглаженного. При не зависящих от времени постоянных коэффициентах этот метод носит название экспоненциального сглаживания. В общем случае фильтр строит некую самонастраивающуюся модель смены состояний системы и на ее основании производит сглаживание. Фильтр Калмана несимметричный, то есть использует только точки, измеренные ранее анализируемой, игнорируя более поздние. Этот факт можно рассматривать как достоинство (оперативное онлайн-сглаживание), а можно и как недостаток (используется только половина информации, доступной для оценки сглаженного значения). Фильтр может быть линейным или нелинейным.
Фильтрация с использованием
преобразования Фурье
Эти линейные фильтры имеют широкое распространение, но следует иметь в виду одно главное свойство большинства фурье-фильтров: они имеют аналог, который выглядит как скользящее взвешенное среднее. Причина в том, что произведение функций в фурье-пространстве превращается в свертку функций (эквивалентную скользящему взвешенному среднему) в исходном пространстве. Неприятное свойство фурье-преобразования состоит в резком уширении фурье-спектра при появлении таких особенностей сигнала, как скачок базовой линии или "срезанная" зашкаливанием детектора вершина пика. Это сильно затрудняет конструирование фильтров, базирующихся на фурье-спектре.
Фильтры, основанные на вейвлет-преобразовании
Вейвлет-преобразование по своим принципам очень похоже на преобразование Фурье, только в качестве объекта берется не вся хроматограмма, а только область пика, а в качестве базисных используются (в случае вейвлетов Эрмита) функция Гаусса и ее производные, а не синусы и косинусы. Такое разложение весьма подходит для отдельно стоящих пиков и не особенно хорошо сглаживает перекрывающиеся пики и базовую линию. Для применения вейвлет-преобразования чрезвычайно полезно заранее знать положение базовой линии и ожидаемую ширину пика.
Адаптивный метод фильтрации шумов
Несколько лет назад мы предложили адаптивный метод фильтрации шумов [7], работающий на принципах локальной минимизации доверительного интервала аппроксимации. Фильтр использует идею, впервые выдвинутую в работе [8]. Для каждой точки хроматограммы подбирается аппроксимационная модель, обеспечивающая минимальный доверительный интервал. Набор функций, по которым строится аппроксимация, может быть выбран исследователем в зависимости от специфики задачи. Мы демонстрировали работоспособность алгоритма на примере полиномиальной аппроксимации.
Для применения адаптивной фильтрации требуется оценка погрешности измерений. В тех случаях, когда она невозможна или неправильна (например, детектор выдает существенно сглаженный сигнал), метод работает плохо. В качестве результата сглаживания в каждой точке выдается два числа: сглаженное значение и его доверительный интервал.
Адаптивное сглаживание нелинейно, так как оптимальная аппроксимация каждой точки зависит от сигнала. Если рассматривать сигнал с шумом как сумму идеального сигнала и шума, то при их раздельной фильтрации сигнал не изменится, а шум может быть аппроксимирован постоянной так, что он в каждой точке уменьшится в раз. Если же сглаживать сумму сигнала и шума, то аппроксимация константой невозможна, и шумоподавление будет зависеть от локального профиля сигнала.
Адаптивный фильтр отличается от известного подхода Колмогорова-Винера [9, 10] тем, что он не пытается найти общую модель, оптимизирующую методом МНК всю экспериментальную кривую. Напротив, он ищет подходящую локальную аппроксимацию, минимизирующую погрешность в каждой точке. Общей модели при этом не существует, поскольку для каждой точки она может оказаться своей. Адаптивный фильтр не подразумевает требований непрерывности и гладкости изучаемой экспериментальной зависимости, которые являются обязательным условием для подхода Колмогорова-Винера.
Новый метод с временно закрытым алгоритмом работы следует в обозримом будущем ожидать в программе "МультиХром".
Особенности методов сглаживания
Для методов, основанных на скользящем взвешенном среднем, возможно сглаживание по четному числу точек, но его результат в большинстве случаев соответствует середине окна, то есть абсцисса при симметричном окне размером 2k попадает в середину между k-й и k+1-й точками разрядной сетки. Поскольку обычно задача фильтрации шумов не подразумевает сдвига разрядной сетки, такой вариант сглаживания крайне редко используется на практике.
Массив исходных (сырых) данных должен быть больше, чем окно фильтрации. Методы, основанные на скользящем взвешенном среднем с окном 2k+1, не могут вычислить аппроксимированное значение в первых и последних k точках массива данных. Для оценки значений в этих точках приходится привлекать дополнительные (выходящие за рамки основного метода сглаживания) соображения.
В случае многократного сглаживания скользящим взвешенным средним отметим интересный эффект, следующий из центральной предельной теоремы теории вероятностей. Если веса скользящего взвешенного среднего неотрицательны и их распределение имеет дисперсию σ2, то, вне зависимости от метода сглаживания, по мере роста числа N проходов, результат будет приближаться к однократному сглаживанию скользящим взвешенным средним с весами, распределенными по закону Гаусса, и дисперсией N ∙ σ2. К примеру, (N > 10)-кратное сглаживание скользящим средним с окном 2k+1 практически эквивалентно однократному сглаживанию гауссианой с дисперсией N(k+1)2/6.
В любом случае, сочетание нескольких проходов взвешенным средним с разным распределением весов будет эквивалентно одному проходу с весами, которые легко вычисляются по весам использованных фильтров. Достаточно сначала выполнить свертку по фильтрам, получив комбинированный, которым затем и выполнять сглаживание функции. Сочетание разных способов сглаживания может иметь смысл, если принципы их работы различны: например, сочетание медианного сглаживания и взвешенного среднего при наличии выбросов в исходных данных.
КАК СРАВНИВАТЬ МЕТОДЫ СГЛАЖИВАНИЯ
Наша предметная область – хроматография, поэтому проведем обработку фрагмента данных, содержащего отдельно стоящий пик характерной формы. Предположим, что частота оцифровки достаточна для вычисления не только площади пика [11], но и других его параметров [12]. Технология сравнения методов сглаживания линейных фильтров, использованная в работе [1] не годится, так как адаптивный фильтр не линеен.
При сравнении способов сглаживания в первую очередь нужно обращать внимание на качество сглаживания хроматограммы в области базовой линии. Чем лучше аппроксимирована базовая линия вне пика, тем более точен ее прогноз для пика. Величина площади практически не зависит от того, как сглажена область пика, поскольку все методы, основанные на скользящем взвешенном среднем, оставляют площадь неизменной. Однако, на нее существенное влияние оказывает то, как проведена базовая линия.
Программа "МультиХром" позволяет проводить обработку многоканальных данных. Мы сгенерировали несколько примеров (табл.1), каждый представляет собой 100-канальную хроматограмму, во всех каналах присутствует пик (рис.2), форма которого описывается экспоненциально модифицированной гауссианой (ЭМГ) [13, 14]:
,
где hG – высота, μG – положение вершины, σG – стандартное отклонение немодифицированной гауссианы; τ – постоянная времени модифицирующей экспоненты; erfcx() – шкалированная дополнительная функция ошибок [15].
Этот пик асимметричный: коэффициент τ/σG = 3, что соответствует коэффициенту асимметрии 2,8. На каждом канале к исходному сигналу был добавлен нормально распределенный шум (рис.2б). Шумы на разных каналах смоделированы независимо, хотя и имели одинаковые параметры генерации. Величина сигмы добавленного шума (E0 = 3 333 единиц) для двух примеров соответствует пределу определения (отношение сигнал / шум, S/N = 10), то есть без фильтрации при большем шуме пик уже будет формально непригоден для количественного расчета. Для третьего примера без фильтрации шумов пик соответствует пределу обнаружения (S/N = 3).
Первый и третий примеры демонстрируют работу алгоритмов фильтрации шумов в условиях недостаточной и избыточной частоты измерений: в первом сигма гауссианы была выбрана равной единице, в третьем сигма гауссианы была увеличена до двадцати, а отношение сигнал / шум уменьшено до трех.
Исследуемую процедуру сглаживания применяли к каждому из ста каналов хроматограммы. Из сглаженной хроматограммы поканально вычитали исходную – без шума (рис.2); в результате получали разностную 100-канальную хроматограмму. В дополнение к 100 профилям исходных данных рассчитывали еще четыре, так называемые "вычисляемые" каналы (рис.3).
1. Es – средний (по 100 каналам) канал разностной хроматограммы – средняя разностная хроматограмма. Отклонение Es от нуля характеризует систематическую погрешность сглаживания (рис.3а), индекс "s" происходит от слова systematic. Систематическая погрешность – это мера смещенности сглаженного значения. Ее нельзя исключить, так как форма сигнала заранее не известна. Она должна учитываться в суммарной погрешности оценки сглаженного значения.
,
где ytrue(i) – истинное значение интенсивности в i-й точке незашумленной хроматограммы; y(i, j) – значение интенсивности в i-й точке хроматограммы на j -м канале после применения процедуры сглаживания; N – общее число точек;
2. Er – стандартное отклонение интенсивности, рассчитанное из 100 сглаженных хроматограмм (СКО). Er является характеристикой случайной погрешности сглаживания (рис.3б), индекс "r" происходит от слова random:
;
3. Et – среднее стандартное отклонение сигнала от нуля, что приблизительно равно квадратному корню из суммы квадратов значений в каналах Es и Er (рис.3в), индекс "t" происходит от слова total:
;
4. R – отношение Er/Et, нормированное на 100 (процент случайной погрешности). По этому отношению можно судить о значимости систематического отклонения в окончательном результате сглаживания (рис.3г):
.
Чем меньше погрешности в первых трех каналах, тем лучше сглаживание в соответствующей точке хроматограммы. Для методов, основанных на скользящем среднем, канал Es, соответствующий бесконечно большой статистике (не 100, а бесконечное число каналов), можно получить, применив соответствующий фильтр к незашумленному пику с последующим вычитанием исходной хроматограммы.
Имея такие графики, можно проводить визуальный и количественный анализ результатов сглаживания. Для количественных оценок мы выбрали два интервала на оси абсцисс: область пика и область базовой линии (рис.2). Область пика (ОП) – тот диапазон номеров точек, для которого отклик исходного сигнала без шума превышает заданную величину (0,1% высоты). Эта область расширена таким образом, чтобы включать возмущения сигнала вне исходного пика, вызванные фильтрацией. В примере 2 мы считали областью пика диапазон от 125-й до 210-й точки. Область базовой линии (ОБ – от 300-й до 500-й точки) – достаточно протяженный интервал точек вне пика, находящийся в промежутке между концами пика и хроматограммы на расстоянии около 100 точек от каждого. В остальных примерах область пика соответствующим образом сужается (для узкого пика) или расширяется (для широкого пика).
Мы выбрали четыре ключевые характеристики метода сглаживания:
1. BNR (baseline noise reduction) – величина, равная отношению исходного стандартного отклонения шума E0 к среднеквадратичному шуму сглаженных данных в области базовой линии; NОБ – число точек в области базовой линии:
;
2. Etmax – максимальное значение среднеквадратичной погрешности в области пика:
;
3. Etav – среднеквадратичное шумоподавление в области пика; Nоп – число точек в области пика:
;
4. Rmin – минимальная величина процента случайной погрешности в области пика:
.
Особенности сглаживания по Савицкому-Голею
В методе Савицкого-Голея существует проблема – провал перед началом пика (рис.4а). Она характерна для полиномиального сглаживания начального участка крутого склона пика и заметна при визуальном исследовании сглаженной хроматограммы, если "глубина" провала превышает 3σ остаточного случайного шума. В противном случае провал маскируют "естественные" шумы базовой линии (рис.4б). Такой эффект может появиться, когда есть точки с отрицательным весом (рис.1). Для скользящего среднего и сглаживания гауссианой таких сложностей в принципе быть не может. Провал существенно затрудняет работу алгоритмов поиска границ пиков и искажает результаты количественной обработки данных. Появление провала перед пиком – только признак проблем. Если построить для этого случая график систематической погрешности (рис.5), то окажется, что максимальное искажение формы пика превышает абсолютную величину провала более чем в пять раз. Максимум амплитуды искажений всех исследованных методов почти всегда приходится на вершину пика.
Случайная погрешность метода Савицкого-Голея (как и всех других методов, основанных на скользящем среднем) для области пика и для базовой линии одинакова. Если мы представим сигнал с шумом как сумму исходного сигнала и шума, то вследствие линейности фильтра случайная погрешность сглаженного значения определяется целиком и полностью шумовой составляющей. Поскольку размеры окна сглаживания не меняются, то и шумоподавление во всех точках одинаково.
Особенности
адаптивного метода фильтрации шумов
Представленный подход к анализу ошибок аппроксимации был применен к адаптивному методу фильтрации шумов. Характерные графики вычисляемых каналов для примера 2 приведены на рис.6. В отличие от методов сглаживания, основанных на скользящем среднем, случайная погрешность перестает быть постоянной по хроматограмме; она возрастает вблизи вершины пика и на крутых склонах. Это явление вызвано автоматической "подстройкой" ширины окна аппроксимации под данные – на базовой линии окно аппроксимации широкое, поэтому случайная погрешность уменьшается, а вблизи вершины пика окно узкое, случайная погрешность увеличивается.
Высокая случайная погрешность в области пика является платой за несмещенность оценки. Адаптивный фильтр просто не может подобрать приемлемую полиномиальную аппроксимацию этой области с большим шумоподавлением. "Хорошее" отношение систематической и случайной погрешностей в области пика является следствием не только уменьшения величины систематической погрешности в адаптивном методе по сравнению с методом Савицкого-Голея, но и менее "агрессивной" фильтрации окрестностей вершины пика.
СРАВНЕНИЕ МЕТОДОВ НА ОДНОМ ГРАФИКЕ
Все методы, работающие как скользящее взвешенное среднее, содержат параметр, называемый "окно сглаживания". Понятие "окно сглаживания" для адаптивного метода не существует или сильно отличается от других; главное, от чего зависят результаты фильтрации – оценка погрешности измерений исходных данных. В качестве параметра адаптивного фильтра, определяющего подавление шумов на базовой линии, использовали максимальный разрешенный размер аппроксимационного окна.
Зависимость шумоподавления от окна соответствующего метода приведена на рис.7. Для облегчения графического сопоставления результатов фильтрации размер окна исключили из рассмотрения; вместо него по оси абсцисс рис.8 откладывали уровень шумоподавления на базовой линии (BNR). Такой вариант представления данных сокращает число параметров и упрощает процедуру сравнения разных методов.
Применимость метода сглаживания зависит от задачи, которую он выполняет. Если требуется несмещенная оценка параметров пика, то систематическая погрешность не должна изменять форму пика в большей степени, чем это делает остаточная случайная погрешность. В этом случае может быть использован простой и наглядный критерий выбора рабочей области фильтра (условия приемлемого искажения пика): процент случайной погрешности для любой точки пика не должен падать ниже 50. При таком условии систематическая погрешность не превышает стандартного отклонения случайной погрешности ни в одной точке пика; это область приемлемого (систематического) искажения пика. Соответственно, при доле случайной погрешности менее 50% искажения пика можно считать неприемлемыми. Пороговая величина может быть и другой, учитывающей особенности аналитической методики [16]; конкретное значение порога не меняет общего подхода к задаче фильтрации шумов.
Несложно увидеть, что метод Савицкого-Голея – лучший из основанных на скользящем среднем, поскольку с его помощью удается достичь наибольшего шумоподавления при условии приемлемого искажения пика (рис.8).
Возникает вопрос: можно ли работать при больших искажениях данных, когда систематическая погрешность существенно превышает случайную? Ответ зависит от того параметра пика, который мы пытаемся оценить. Если речь идет о площади пика, то при известной базовой линии она остается почти неизменной при любом методе сглаживания, основанном на скользящем среднем. Если провал перед пиком отсутствует, погрешность площади главным образом зависит от того, насколько хорошо подавлены шумы базовой линии. Соответственно, несмотря на искажение формы пика, целевой параметр – площадь – при увеличении окна сглаживания определяется с большей точностью. Если же нужно определить высоту или ширину на половине высоты, то результат будет неудовлетворительным. Правда, если истинная высота не нужна, а параметры сглаживания одинаковы в градуировочных графиках и в анализе, то за счет уменьшения случайной погрешности высоты сглаженной хроматограммы количественная оценка концентрации по высоте может быть достаточно точной, несмотря на искажение пика. В этом смысле сглаживание гауссианой предпочтительнее метода Савицкого-Голея: ширина окна не имеет "верхнего порога", связанного с провалом перед пиком. Следует помнить, что фильтрация увеличивает ширину пиков (при сглаживании гауссианой дисперсия пика равна сумме дисперсий пика и гауссианы), что уменьшает разрешение пиков и оценку эффективности хроматографической колонки. Высота пика, сглаженного скользящим взвешенным средним, при проценте случайной погрешности менее 50 становится более воспроизводимой, но смещенной оценкой истинной величины.
Поскольку наша реализация адаптивного метода основана на полиномиальной аппроксимации, естественно сравнить его с методом Савицкого-Голея, тем более что он оказался лучшим из методов скользящего среднего (см. рис.8). В методе Савицкого-Голея неприемлемое искажение пика в примере 2 появляется при подавлении шумов на базовой линии BNR ≈ 2,5 (место пересечения кривой Rmin с горизонталью на уровне 50%). В адаптивном методе искажение такого же порядка возникает при подавлении шумов на базовой линии BNR ≈ 4,7 и незначительно изменяется при еще большем шумоподавлении. Таким образом, отношение S/N адаптивного метода превышает этот параметр оптимально настроенного метода Савицкого-Голея:
• для первого примера – в 4,4/1,4 = 3,15 раза или на 20 lg(3,15) ≈ 9,9 дБ;
• для второго примера – в 4,7/2,5 = 1,88 раза или на 20 lg(1,88) ≈ 5,5 дБ.
Для третьего примера S/N не изменяется, если допустить увеличение систематической погрешности до удвоенной остаточной случайной (процент случайной – 33), то 11/6 = 1,83 · 20 lg(1,83) ≈ 5,3 дБ; такое же значение наблюдается для нового метода без ослабления критерия.
Для некоторых примеров есть небольшая область, где метод Савицкого-Голея дает лучшую аппроксимацию области пика, чем адаптивный. Она соответствует достаточно малым (менее 2,5) коэффициентам шумоподавления на базовой линии и весьма узка.
При малой ширине пика (пример 1) происходит ожидаемое – ни один из методов обработки данных, основанных на скользящем среднем, не попадает в область приемлемого искажения пика, за исключением фильтра Савицкого-Голея по пяти точкам. Этот фильтр обеспечивает падение уровня шумов базовой линии примерно в 1,4 раза. Адаптивный фильтр продолжает выдавать приемлемое шумоподавление вплоть до падения шумов базы в пять раз. При этом реального уменьшения максимальной погрешности внутри пика нет ни у одного метода – все значения коэффициентов шумоподавления E0/ Etmax находятся ниже единицы, то есть максимальная суммарная погрешность в области пика для всех методов больше исходной. Тем не менее, применение адаптивного метода дает выигрыш в погрешности вычисления площади и высоты за счет хорошего шумоподавления на базовой линии – она оказывается проведенной более точно.
Наиболее актуальна фильтрация шумов при работе вблизи пределов определения и обнаружения [17]; именно такие примеры мы анализировали в настоящей работе (табл.2). Особенно заметно преодоление пиком (пример 3) предела определения (S/N = 10), как в случае адаптивного, так и метода Савицкого-Голея, тогда как без фильтрации шумов этот пик находился на грани предела обнаружения (S/N = 3). Результаты адаптивного фильтра и метода Савицкого-Голея для этого пика совпадают вплоть до BNR = 6, далее шумоподавление в области пика адаптивного и нового фильтров заметно лучше, чем у метода Савицкого-Голея.
Предельно допустимое шумоподавление Савицкого-Голея снижается с уменьшением ширины пика и с ростом отношения сигнал / шум. При фильтрации шумов реальной хроматограммы вполне возможна ситуация, когда в ней присутствуют пики разной ширины и высоты. В этом случае приходится выбирать параметр сглаживания таким, чтобы для самого узкого и высокого пика искажения были в рамках приемлемости. Это касается всех алгоритмов, но для адаптивного и нового фильтров можно добиться того, чтобы шумоподавление на базовой линии было существенно (в 2–3 раза) выше, чем в случае применения метода Савицкого-Голея.
Главная трудность в применении адаптивного фильтра состоит в том, что для его нормальной работы требуется знание погрешности каждого измерения. В рамках нашей модели эта погрешность считается одинаковой для всех точек хроматограммы (гомоскедастичной). Реальные погрешности для разных точек могут отличаться, тогда модель адаптивного метода необходимо приспособить к неодинаковому (гетероскедастичному) шуму [18], что заметно усложнит вычисления. Обычно считается, что погрешность измерения зависит от уровня сигнала. Поскольку основной эффект сглаживания достигается за счет фильтрации шума базовой линии, а не области пика, использование адаптивным фильтром оценки шума базовой линии, как описано в [2, 7], дает приемлемые результаты, а погрешность базовой линии вполне может считаться гомоскедастичной. Ее можно оценить по всей хроматограмме, но при условии, что ранее к данным не применяли методы сглаживания. Многие производители оборудования встраивают в детекторы собственные программы подавления шумов, существенно искажающие полученные экспериментальные результаты. Тогда оценить погрешности возможно, установив определенные настройки детектора, например, минимально возможную постоянную времени.
Формальное фармакопейное определение отношения сигнал / шум [17, 19] в случае применения адаптивного метода способно вырасти в 2–5 раз в зависимости от формы сигнала и качества базовой линии. Максимальное доступное адаптивному методу шумоподавление может оказаться ниже тех чисел, которые мы получили на графиках, поскольку оно зависит от сглаживаемых данных. Если продолжительных участков базовой линии нет, то потенциальные возможности адаптивного метода могут оказаться нереализованными.
Новый метод фильтрации шумов является достойным соперником адаптивного фильтра, поскольку, как видно из представленных в этой работе результатов (см. рис.8), по всем критериям работает устойчивее адаптивного фильтра.
ВЫВОДЫ
Предложен подход к сравнению методов сглаживания, основанный на моделировании измерительной системы и ее шумов. Модельная система позволяет построить кривые, отражающие характеристики метода сглаживания: систематическую, (остаточную) случайную и суммарную погрешности, а также долю случайной погрешности в суммарной. На основании кривых рассчитываются интегральные характеристики метода сглаживания, для хроматографического сигнала – это шумоподавление на базовой линии, минимальное и среднее шумоподавление по пику, а также минимальный процент случайной погрешности в области пика. Критерии, основанные на анализе интегральных параметров, помогают выбрать метод, наиболее подходящий для конкретной системы.
ЛИТЕРАТУРА
1. Enke C.G., Nieman T.A. Signal-to-noise ratio enhancement by least-squares polynomial smoothing // Anal. Chem. 1976. 48. P. 705A–712A. doi:10.1021/ac50002a007.
2. Каламбет Ю.А., Михайлова К.В. Оценка величины шума и ее использование при обработке хроматографического сигнала // Лабораторный Журнал. 2002. 1. С. 32–35. http://multichrom.ru/Docs/ots-vel-shuma.pdf.
3. Felinger A. Data analysis and signal processing in chromatography // Elsevier, 1998. V. 21. 413 p.
4. Brereton R.G. Chemometrics: Data Analysis for the Laboratory and Chemical Plant. Wiley, 2003. 504 p. doi:10.1198/tech.2004.s738.
5. Savitzky A., Golay M.J.E. Smoothing and Differentiation of Data by Simplified Least Squares Procedures // Anal. Chem. 1964. 36. P. 1627–1639. doi:10.1021/ac60214a047.
6. Kalman R.E. A New Approach to Linear Filtering and Prediction Problems // J. Basic Eng. 1960. 82. P. 35–45. doi:10.1115/1.3662552.
7. Каламбет Ю., Мальцев С., Козьмин Ю. Фильтрация шумов: окончательное решение проблемы // АНАЛИТИКА. 2011. 1. P. 50–55.
8. Ariyur K.B., Jelinek J. Aeroengine Prognostics via Local Linear Smoothing, Filtering and Prediction // SAE Trans. J. Aerosp. 2004. 113. P. 1773–1780. doi:10.4271/2004-01-3160.
9. Колмогоров А.Н. Интерполирование и экстраполирование стационарных случайных последовательностей // Изв. АН СССР. Сер. Матем. 1941. 5. С. 3–14.
10. Wiener N. The Extrapolation, Interpolation and Smoothing of Stationary Time Series, New York: Wiley, 1949. 176 p.
11. Kalambet Y., Samokhin A., Kozmin Y. Finally, how many ponts per peak?, in: X Winter Symp. Chemom., Samara, 2016. doi:10.13140/RG.2.2.28145.28006.
12. Wahab M.F., Dasgupta P.K., Kadjo A.F., Armstrong D.W. Sampling frequency, response times and embedded signal filtration in fast, high efficiency liquid chromatography: A tutorial // Anal. Chim. Acta. 2016. 907. P. 31–44. doi:10.1016/j.aca.2015.11.043.
13. Каламбет Ю., Михайлова К., Козьмин Ю., Нагаев И. Восстановление хроматографических пиков при помощи экспоненциально-модифицированной гауссианы // АНАЛИТИКА. 2012. № 4. С. 62–67.
14. Kalambet Y., Kozmin Y., Mikhailova K., Nagaev I., Tikhonov P. Reconstruction of chromatographic peaks using the exponentially modified Gaussian function // J. Chemom. 2011. 25. P. 352–356. doi:10.1002/cem.1343.
15. Cody W.J. Rational chebyshev approximations for the error function // Math. Comput. 1969. 23. P. 631–637. doi:10.1090/S0025-5718-1969-0247736-4.
16. Дворкин В.И. Метрология и обеспечение качества химического анализа // М: Издательство МИТХТ, 2014. 416 c.
17. Государственная фармакопея Российской Федерации. XIII издание, Федеральная электронная медицинская библиотека, 2015.
18. Каламбет Ю.А., Мальцев С.А., Козьмин Ю.П. Доверительные интервалы метода взвешенных наименьших квадратов и стратегия градуировки // Заводская Лаборатория. Диагностика Материалов. 2015. 81. С. 69–76. http://www.multichrom.ru/Docs/zav-lab15-1-1-81.pdf.
19. Ettre L.S. Nomenclature for chromatography (IUPAC Recommendations 1993) // Pure Appl. Chem. 1993. 65. P. 819–872. doi:10.1351/pac199365040819.
Отзывы читателей
