CC BY
18ISSN 0868-5886
НАУЧНОЕ ПРИБОРОСТРОЕНИЕ, 2009, том 19, № 3, c. 35-40
ТЕОРЕТИЧЕСКИЕ ИССЛЕДОВАНИЯ
УДК 621.391.26
© В. В. Манойлов, И. В. Заруцкий
ИССЛЕДОВАНИЕ АЛГОРИТМОВ ОТБРАКОВКИ ВЫБРОСОВ В МАСС-СПЕКТРОМЕТРИЧЕСКИХ СИГНАЛАХ
Описываются результаты исследования нескольких типов алгоритмов отбраковки выбросов и сглаживания в масс-спектрометрических сигналах. Исследование проведено с помощью вычислительного эксперимента с использованием компьютерных моделей масс-спектрометрических сигналов, искаженных шумами и выбросами. Показаны преимущества модифицированного для масс-спектрометрии алгоритма Rousseeuw, основанного на вычислении квадратов медиан разности текущих отсчетов сигналов.
Кл. сл.: методы обработки сигналов, масс-спектрометрия, статистический анализ данных, сглаживание и фильтрация сигналов
ПОСТАНОВКА ЗАДАЧИ
Одной из важных задач обработки сигналов масс-спектрометрического анализа является удаление выбросов, которое имеет существенное значение для повышения точности измерений. Мы наблюдаем мультимодальную функцию (спектр) в аддитивной смеси со стационарным шумом и выбросами (ложными наблюдениями). Требуется освободиться от выбросов и оценить дисперсию (среднеквадратическое отклонение) шума при минимально возможном искажении мультимодаль-ной функции.
Под "выбросами", или "ложными элементами", понимают данные, сильно отличающиеся по величине от математического ожидания анализируемой выборки случайной величины. Плотность распределения данных случайных величин, в которых имеются выбросы ^ обычно представляется в виде распределения Тьюки, которое записывается в следующем виде:
х = /(1 - а) + £а,
где f — плотность распределения данных, в которых отсутствуют выбросы; У2 — плотность распределения данных, которые являются выбросами, а — относительное количество выбросов в анализируемой выборке данных. Например, если в анализируемой выборке данных содержится 10 % выбросов, то а = 0.1, если 5 %, то а = 0.05. В случае нормального распределения анализируемых данных
X = N(т,СТ1Х1 - а) + N(т2,а2)а, где f1 и У2 представляют собой функции Гаусса
со средними значениями соответственно т1, т2 и средними квадратичными отклонениями ст1, сг2. Очень часто выбросами являются данные, плотность распределения которых имеет среднее значение т2 = т1, а средние квадратичные отклонения отличаются в k раз, т. е. ст2 = . При этом число k может быть существенно больше 3. Для отбраковки выбросов наиболее простым и достаточно надежным методом является метод "3 а" (трех сигм). При наличии в выборке данных выбросов оценка значения величины а может быть искажена. Описываемые ниже алгоритмы показывают, каким образом в задачах масс-спектрометри-ческого анализа производится оценка значения величины а .
В работе Rousseeuw Р.Л. [1] дается изящный и простой метод отыскания а , основанный на определении mmmed(ri2), где ri 2— квадрат разности экспериментальных наблюдений в выборке заданного размера N [у1...yi...ум.], называемой скользящим окном; М = minmed(ri2) является статистикой, для которой наблюдение у{ = а + М является верным (не выбросом); а — оценка функции (в работе [1] константы), на которую наложены шум и выбросы. Для нахождения minmed(ri2) необходимо:
- построить вариационный ряд, т. е. упорядочить наблюдения по величине
У1 < У2 < . . . < ум ;
- составить разности
Дй=У[м/2]+й- Уй, й = 1, 2,..., [N/2];
- найти среди них минимальную Дй0
36
В. В. МАНОЙЛОВ, И. В. ЗАРУЦКИЙ
и тогда соответствующая этой разности полусумма (у[лг/2]+м + ук0)/2 будет являться оценкой функции а , а разность (уЛ/2]+м - ук0) является оценкой статистики Ы, которая реализует minmed(гI2).
Все же основная идея метода минимума квадрата медиального отклонения для нашей задачи должна быть в определенной степени доработана. Основная идея доработки метода минимума квадрата медианного отклонения для масс-спектро-метрии принадлежит Сирвидасу С.И. [2]. Суть доработки состоит в том, что масс-спектр является быстроизменяющейся функцией, следовательно, нужна оценка текущего значения первой и второй производной и, кроме того, требуется дополнение в виде процедуры типа скользящего окна, т. е. внедрение локального сглаживания [2] и определение оптимальных размеров скользящего окна по критерию минимальной суммарной дисперсии. Ниже описывается доработанный метод, представленный в виде алгоритма, и результаты исследования его возможностей в сравнении с другими алгоритмами, решающими задачу отбраковки выбросов.
АЛГОРИТМ МЕТОДА МИНИМУМА КВАДРАТА МЕДИАННОГО ОТКЛОНЕНИЯ, МОДИФИЦИРОВАННЫЙ ДЛЯ ОБРАБОТКИ МАСС-СПЕКТРОМЕТРИЧЕСКИХ СИГНАЛОВ
1. Ввод данных.
2. Формируем массив х1, 7 = [1, п].
3. Цикл 7 от 1 до п.
4. Цикл м (ширина окна) от 3 до wmsx (м — нечетное).
5. y. = x[i + j], j = 1,w .
6. Вычисляем dy. =
yj+2 - yj 2
7. Строим вариационный ряд dyj .
8. Вычисление у1 = meddУj..
9. Вычисляем d2 у. = -——
10. Строим вариационный ряд d2у. .
11. Вычисляем y2 = med d2y. .
12. Вычисляем zi = y. - y:( j - h) - y
2 (j - h)2
где h =
13. Строим вариационный ряд .
14. Вычисляем такое, что 5[+ к +1]
: min(Z[ j + h +1] - Z[ j]); j = [1, h].
Z[ jo + h +1] + Z[ j0]
15. Вычисляем =-
16. Вычисляем
2
=
(Z[ jo + h +1] - Z[ jo])2 1483
Г] ~ам\
17. Проверим условие --- < 3.
18. Если п. 17 выполнен, заносим индекс . в массив I.
19. Вычисляем — число элементов в массиве I.
20. Вычисляем b = — V z .
w S J
21. Вычисляем г = V (5 - Ь )2.
ч / : V . м /
22. Для всех . £ I и 1 < 7 + . < Л увеличиваем элемент массива J[7 +.] на 1.
23. Повторим шаги 5-22 для следующих м .
24. Определяем ч0 такое, что гщ = т1п гм .
25. Значение Ьм, соответствующее гмт1п, заносим в новый массив "чистых" данных X. = Ь .
7 м0
26. Повторяем шаги 4-26 для следующих I.
Поясним некоторые шаги алгоритма.
Шаг 12. Для работы данного метода в масс-спектрометрии необходимо оценить значение функции в данной точке. Имея значения первой и второй производных в каждой точке, можно воспользоваться разложением Тейлора
X X2
у( х) = у( хо) + у'( хо)1! + у"С xo)^!...,
откуда значение у( х0) есть разность
X X2
у(xo) = у (х) - у'( - у "(:с0)2\.
Таким образом, получаем оценку функции в окне размеромм (м = 2k +1 для k = 0,±1,±2,±3...±k ):
5 2
а = 5(5) = ул+5 - у""5)5 - у"5) — -const.
Шаг 16. Оценка среднего квадратичного значения а основывается на том факте, что интервал [а - Ы, а + Ы] содержит половину наблюдений у{ ("верные" наблюдения) и соответственно интервал [-Ы, Ы] содержит половину остатков г . Из теории вероятности известно, что радиус интервала, в который попадает половина значений случайной величины — вероятное отклонение
2
2
(В.О.) от математического ожидания нормальной величины с дисперсией а2 равен 0.6745 а, т. к.
значение функции ошибок erf(0.6745 / V2) = 0.5, что означает, что вероятность
Р(|х| < 0.6745а) = 0.5 при а = 1.
M
Отсюда оценка а : а =-= 1.483M .
0.6745
Адаптивность алгоритма заключается в том, что из всех полученных значений bw при различной ширине окна W в качестве оценки берется значение с индексом W0, которое доставляет
min J(y - bw )2 .
w i
После окончания работы алгоритма в массиве J для каждого элемента массива исходных данных хранится количество его упоминаний в качестве выброса i.
Ниже описываются результаты исследования доработанного для масс-спектрометрических сигналов метод отбраковки выбросов.
МОДЕЛИРОВАНИЕ "ЗАГРЯЗНЕННОГО" СИГНАЛА И ОБРАБОТКА ДАННОГО СИГНАЛА С ПРИМЕНЕНИЕМ РАЗЛИЧНЫХ МЕТОДОВ ОТБРАКОВКИ ВЫБРОСОВ
В основе модели лежит гауссова функция вида
120 14Ü 360 1В0 200 220 240 260 200
Рис. 1. Исходный сигнал с ложными выбросами
Параметры функции вводятся в программу с последующей записью результатов. Далее моделируется случайный шум с нормальным распределением с заданным СКО. В основе модели создания шума лежит функция формирования случайного числа. Суммируются 12 случайных чисел, из суммы вычитается 6, полученное число умножается на СКО — результатом является шум, который добавляется к значениям функции уг-. Далее увеличиваем i на 1 и повторяем процедуру.
В основе модели создания выбросов лежит функция формирования случайного числа. Из датчика случайных чисел берется число и, если оно меньше, чем процент выбросов, то к "зашумлен-ной" уже функции при текущем значении i прибавляет выброс с заданной амплитудой. После этого производится "очистка" сигнала при помощи метода параметрического сглаживания. Параметры модели гауссовой функции:
А=30; ¡л = 3 ; t = 0.
Для проведения математического эксперимента на языке программирования С++ написана программа. В качестве функций для тестирования выбраны следующие три варианта:
1) исходная гауссова функция (параметры указаны выше);
2) исходная функция с наложенным шумом и
3) с выбросами.
СКО шума — 0.0005, % выбросов (% выбр.) и их амплитуда (Ампл.) вводились в программу в диалоговом режиме.
Рис. 2. Сигнал после отбраковки ложных выбросов алгоритмом минимизации квадрата медианного остатка (МКМО)
о
38
В. В. МАНОЙЛОВ, И. В. ЗАРУЦКИЙ
........ ........ ........ ........ .......п
-------- -------- -------- -------- -------- -------- --------
120 Щ 130 130 200 2 2D 240 260 280
Рис. 3. Сигнал после отбраковки ложных выбросов алгоритмом медианы в скользящем окне (МСО) — искажена вершина
120 140 160 ISO 2Ш 220 240 2Ю. 2Ш
Рис. 4. Сигнал после отбраковки ложных выбросов цифровым фильтром Чебышева — искажена форма
Средние значения невязок от "идеального" сигнала для сигналов, обработанных разными алгоритмами, при различных количествах и амплитуде ложных выбросов
Параметры сигнала Методы обработки
% выбросов Амплитуда выбросов Мкмо Мсо Сглаживание Фильтр Чебышева
25 20 2.34 2.97 9.22 4.88
25 30 2.52 2.98 9.47 12.31
25 40 2.30 2.93 10.19 12.33
35 20 2.32 2.95 6.51 10.24
35 30 2.42 2.50 9.00 12.96
35 40 2.31 3.20 12.70 12.27
40 20 2.34 2.68 6.21 10.31
40 30 2.40 2.70 8.09 11.43
40 40 2.74 2.84 13.32 13.25
АНАЛИЗ РЕЗУЛЬТАТОВ МОДЕЛИРОВАНИЯ
Результаты эксперимента (графики) — моделирование функции Гаусса, наложение шумов и выбросов, фильтрация зашумленной функции — представлены на рис. 1-4 и в таблице. Обозначены:
Мкмо — алгоритм минимизации квадрата медианного остатка,
Мсо — алгоритм медианы в скользящем окне, Сглаживание — сглаживание в скользящем окне суммированием с "весами",
Фильтр Чебышева — цифровой фильтр Чебы-шева.
Из данных в таблице видно, что по сравнению с другими методами Мкмо имеет наименьшие значения средних значений невязок от "идеального" сигнала. Поэтому можно сделать вывод о том, что этот метод меньшего всего искажает "идеальный" сигнал.
Метод Мкмо оставляет на графике четко различимый пик, позволяющий дальнейший анализ причин его возникновения. Этот пик в дальнейшем может быть сглажен при повторном применении метода Мкмо.
ЗАКЛЮЧЕНИЕ
Рассмотренные выше алгоритмы, основанные на поиске минимума квадрата медианного остатка, являются универсальными для решения большинства задач обработки информации масс-спектро-метрического изотопного анализа. Для ряда частных задач могут быть использованы более простые алгоритмы, описанные в работах [3, 4, 5, 6,
7].
СПИСОК ЛИТЕРАТУРЫ
1. Rousseeuw P.J. Least Median of Squares Regression // Journal of the American Statistical Association. 1984. V. 79. P. 871-880.
2. Sirvidas S.I., Manoylov V.V. The Software Outliers' Eliminator and Noise Smoother for Spectral Data // 7th International school-seminar on automation and computing in science, engineering and industry. РАН, МГУ. Ялта, 1996. P. 32.
3. Манойлов В.В., Заруцкий И.В. Отбраковка "выбросов" и оценка параметров масс-спектро-метрических сигналов для прецизионного изотопного анализа // Научное приборостроение. 2002. Т. 12, № 3. С. 67-70.
4. Манойлов В.В., Заруцкий И.В. Алгоритмы первичной обработки масс-спектрометри-
ческих сигналов для прецизионного изотопного анализа // Вопросы атомной науки и техники. Серия "Техническая физика и автоматизация". Научно-технический сборник. Вып. 56. М.: Мин-во РФ по атомной энергии, Центральный научно-исследовательский институт информации и технико-экономических исследований в атомной науке и технике, 2002. С.52-74.
5. Манойлов В.В., Заруцкий И.В. Алгоритмы обработки масс-спектрометрических сигналов для изотопного и химического анализа // Труды LVII Научной сессии, посвященной Дню радио. М.: Российское научно-техническое общество радиотехники, 2002. Т. 1. С. 274277.
6. Манойлов В.В., Аракелянц М.М., Чернышев И.В., Сердюк Г.И. Программное обеспечение измерительно-вычислительной системы на основе IBM/PC-AT для определения возраста геологических образцов калий-аргоновым методом на масс-спектрометре МИ-1201ИГ // 12-й Международный симпозиум по проблемам модульных информационно-вычислительных систем и сетей. Тезисы докладов 3.7. РАН, МГУ. М., СПб., 1997. С. 43.
7. Манойлов В.В., Аракелянц М.М., Чернышев И.В. Программное обеспечение для определения изотопного состава аргона в автоматизированном комплексе на базе масс-спектрометра МИ1201ИГ // Научное приборостроение. 1999. Т. 9, № 4. С. 84-95.
Институт аналитического приборостроения РАН, Санкт-Петербург
Материал поступил в редакцию 19.05.2009.
THE INVESTIGATION OF ALGORITHMS OF OUTLIERS' ELIMINATOR IN MASS-SPECTRAL SIGNALS
V. V. Manoylov, I. V. Zarutsky
Institute for Analytical Instrumentation RAS, Saint-Petersburg
The results of investigation of some types of algorithms of outliers eliminator and smoothing in mass-spectral signal are discussed. The investigation was performed with the help of computational experiment on the
40
B. B. MAHOHHOB, H. B. 3APYU,KHH
base models of mass-spectral signals, distorted by noise and outliers. The advantages of the modified Rous-seeuw's method for mass-spectral analysis are shown. This method is based on the least median of squares regression.
Keywords: method for signal treatment, mass-spectometry, statistical data analysis, signal smoothing and filtration
HAyHHOE nPHEOPOCTPOEHHE, 2009, tom 19, № 3
