Отбраковка "выбросов" и оценка параметров масс-спектрометрических сигналов для прецизионного изотопного анализа Текст научной статьи по специальности «Компьютерные и информационные науки»

ISSNG868-5886

НАУЧНОЕ ПРИБОРОСТРОЕНИЕ, 2002, том 12, № 3, c. 38-46

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

УДК 621.384.8

© В. В. Манойлов, И. В. Заруцкий

ОТБРАКОВКА "ВЫБРОСОВ" И ОЦЕНКА ПАРАМЕТРОВ МАСС-СПЕКТРОМЕТРИЧЕСКИХ СИГНАЛОВ ДЛЯ ПРЕЦИЗИОННОГО ИЗОТОПНОГО АНАЛИЗА

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

обнаружения и оценки параметров одиночных пиков.

ВВЕДЕНИЕ

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

АЛГОРИТМЫ ОТБРАКОВКИ "ВЫБРОСОВ" В МАССИВАХ ДАННЫХ ИЗОТОПНЫХ МАСС-СПЕКТРОВ

Под "выбросами", или "ложными" элементами, понимают данные, сильно отличающиеся по величине от математического ожидания анализируемой выборки случайной величины. Плотность распределения данных случайных величин, в которых имеются "выбросы", / обычно представляется в виде распределения Тьюки, которое записывается в следующем виде:

fv = /1(1 -а) + f2a,

где /1 — плотность распределения данных, в которых отсутствуют "выбросы", /2 — плотность распределения данных, которые являются "выбросами", а — относительное количество "выбросов" в анализируемой выборке данных. Например, если в анализируемой выборке данных содержится 10 % "выбросов", то а = 0.1, если 5 %, то а = 0.05. В случае нормального распределения анализируемых данных

/V = Ш(т1,а1)(1 -а) + Ы(ш2,02)а, (1)

где /1 и /2 представляют собой функции Гаусса со средними значениями соответственно т1, т2 и средними квадратичными отклонениями ст1, а2. Очень часто "выбросами" являются данные, плотность распределения которых имеет среднее значение т2 = т1 , а средние квадратичные отклонения отличаются в к раз, то есть а2 = ко1 . При

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

Алгоритм с поочередным исключением измерений

Мы полагаем, что измерения являются линейной регрессией

Ж) = А + Вг. + Ж(г.), (/ = 1,...,п), (2)

где Ж (г.) — напряжение шума.

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

ходного данного вычисляется значение его среднего квадратического отклонения от линии регрессии. В случае превышения отклонением значения 35/^ - 1), где N — объем выборки, данное измерение заменяется значением, "лежащим"

на линии регрессии. После этого по полученным новым данным вычисляются коэффициенты А и В. Данный алгоритм внедрен в программное обеспечение масс-спектрометров МИ1201ИГ [1] в Институте геологии рудных месторождений и в масс-спектрометре МИ1201В в Объединенном Институте геологии, геофизики и минералогии СО РАН.

Алгоритм на основе метода квадрата

минимального медианного остатка

Для того, чтобы использовать метод квадрата минимального медианного остатка [1], предлагается рассматривать преобразованные измерения г(гг) = у(гг) - Вг{ и затем искать

/(В) = шт {шеф(г) - Дг^)..^))2]},

где Т (2(г1)... г(гп)) — оценка параметра А, минимизирующая квадрат медианного остатка. Для нахождения В в качестве первого его приближения возьмем медианную оценку производной исходных измерений у(гг), т. е.

(

В = med

Ж+і) - Ж-1)

у гг+1 гг—1 J

В дальнейшем выполняются операции по схеме оценки методом квадрата минимального медианного остатка для переменной

2(г. ) = У(гт ) — Щ — гт ^ где У(гт ) = ШеЙ(У ({г )) .

1. Строим вариационный ряд:

2(0 < 2(г2) < ... < 2(гк) < ... < 2(^2к+1); п = 2к+1.

2. Образуем следующие к разностей:

А1 = г(гк+2) — 2(г1);

А 2 = г(гк+3 ) — г(г2 );

Ак = г(г2к+1 ) — г(гк ).

3. Находим минимальную разность

А^0 <А Л < ••• <Аг .

4. Находим оценку параметра А

2(Ам+к+1 ) + 2(г Л0 )

А = ■

2

нее квадратичное отклонение 8 = 1.4835л/а2 .

6. Вычисляем абсолютные величины остатков: г(г) = г(гг) — А и проводим отбраковку "пог (г)

дозрительных" измерений. Если -------> 3 , то г —

8

отсчет "подозрительный". Этот отсчет заменяется отсчетом "лежащим" на линии регрессии для данного времени измерения. Полученный ряд измерений у(г)) обрабатываем по схеме наименьших

квадратов ^ (г(гг) — А — Вгг) ^ шт.

г

7. Находим новые оценки параметров А и В и окончательную оценку дисперсии.

Проверка данного алгоритма на моделях с различным отношением сигнал/шум показала результаты, представленные в табл. 1. На рис. 1 представлен пример восстановления линейного сигнала, в котором 50 % "выбросов".

Алгоритмы отбраковки "выбросов" в масс-спектрометрическом изотопном анализе веществ в твердой фазе

Отбраковка "выбросов" производится на следующих этапах проведения автоматизированного изотопного анализа веществ в твердой фазе:

1. при определении среднего значения по результатам измерений на каждом из пиков спектра;

2. при определении средних значений изотопных отношений по серии спектров;

3. при определении средних значений изотопных отношений по нескольким сериям спектров.

Табл. 1. Результаты отбраковки "выбросов" алгоритмом на основе метода квадрата минимального медианного остатка

СКО шума, % Заданное количество выбросов, % Отбракованное количество выбросов, %

0 11 11

20 40 40

0...5 50 50

5...10 44 44

5.10 50 48

0 2 О 50 46

33 25 25

5. Находим дисперсию шума измерений через медианное отклонение а2 = (А — г (г Л0))2 и сред-

Рис. 1. Влияние односторонних "выбросов" на параметры линейного сигнала. Тонкой линией показан искаженный "выбросами" сигнал, а жирной — сигнал, восстановленный после удаления "выбросов"

Для каждого из этих этапов предлагается простая и достаточно надежная процедура отбраковки "выбросов".

1. Строится вариационный ряд из обрабатываемых данных.

2. Из всех данных отбрасывается а % данных, находящихся в верхней и нижней частях вариационного ряда, где а — априорный процент "выбросов" в соответствии с плотностью распределения, описанной формулой (1).

3. По оставшимся данным определяется среднее значение т и среднее квадратичное отклонение 7.

4. Дальнейшая процедура отбраковки "выбросов" является стандартной по критерию к7 , где к определяется по таблицам распределения Стью-дента в зависимости от объема выборки.

На рис. 2 показывается, какие из обрабатываемых данных участвуют в определении среднего значения т и среднего квадратичного отклонения

7. Проверка работы данного достаточно простого алгоритма на масс-спектрометрах МИ 1201Т

и МИ1320 [2] в Всероссийском геологическом институте показала возможность отбраковки до 15 % "выбросов".

Рис. 2. Упорядоченная выборка обрабатываемых данных, построенная в порядке возрастания чисел. Темным фоном показаны числа, которые не участвуют в определении среднего и среднего квадратичного значений выборки

настройка на опорный пик

Время настройки на центр пика, выполняемой во всех масс-спектрометрах с дискретной ступенчатой разверткой, можно сократить при использовании описываемых ниже алгоритмов. Без применения этих алгоритмов время получения оценки положения середины пика при постоянной времени измерительной системы, равной 0.2 с, достигает 80 с (резистор в цепи обратной связи электро-ме1т2рического усилителя имеет сопротивление 10 Ом). В применяемых до сих пор алгоритмах необходимо проводить измерение 50-60 значений ионного тока в области предполагаемого положения центра пика при разных величинах напряженности магнитного поля [1]. Время, в течение которого удается проделать необходимое количество измерений, столь велико еще и по той причине, что для снижения влияния динамических ошибок перед измерением каждого нового значения ионного тока необходимо дождаться завершения переходных процессов в системах прибора. Необходимое для оценки положения центра пика количество исходных данных и, следовательно, общее время настройки можно существенно сократить, если найти модель пика, удачно (с точки зрения минимума средней квадратичной ошибки) описывающей пик. В этом случае можно оценить положение центра и другие параметры пика лишь по нескольким значениям ионного тока, примерно равномерно распределенным по области, где находится пик. А если при вычислениях учитывать динамические свойства измерительной системы, то становится возможным снимать значения (ионного тока), не дожидаясь окончания переходных процессов, что еще несколько снизит общее время оценки положения центра пика. Уменьшение времени работы программы позволяет снизить влияние нестабильных процессов в источнике ионов, ионно-оптической системе и в системе управления магнитным полем.

Модели пика для изотопных масс-спектрометров

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

) = ехр

и прямоугольной функции

2а2

гее* (V) =

1,

0,

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

(

Ф

V V

а

+ф

а

л

\

(3)

где Ф — функция Лапласа, А — масштабный множитель.

Графики функции i(t) (при 11 = 60.9, Ь = 108.5, А = 256.3, с = 7.5) и пика изотопа стронция 88

/с 88ч ->

(Ьг ) приведены на рис. 3.

Используя указанную модель, произведем оценку положения центра пика. Так как масс-спектрометр относится к классу так называемых линейных приборов [3], то выходной сигнал может быть представлен в виде свертки (при отсут-

ствии шума и при пренебрежении дискретным характером измерений на конечном интервале)

I(V) = | — т)і(т)<іт .

(4)

Здесь t и т — параметры времени при управлении напряженностью магнитного поля. Такой переход оказывается допустимым, т. к. напряженность магнитного поля в пределах одного пика можно считать линейной функцией от времени. Естественно, встает вопрос о восстановлении истинной зависимости i(t) по измеренным значениям I(^ решением уравнения (4).

указаны в тексте) и пика изотопа стронция 8г 88

Трудность состоит в следующем. Данные I ) реально регистрируются с некоторой погрешностью, а задача (4) относится к классу обратных, некорректных в классическом смысле задач [3], поэтому возможно получение неустойчивого, осциллирующего решения. Действительно, если формально задачу (4) решать методом так называемой деконволюции с использованием преобразования Фурье, то в результате при практических вычислениях получается осциллирующая кривая, не являющаяся искомым решением [4]. Для нахождения устойчивых решений обратных задач в вычислительной математике наиболее мощными методами являются регуляризация по Тихонову и оптимальная фильтрация Винера [5]. Поскольку в регистрируемых данных масс-спектрометра I (^) неизбежно присутствует шум, то требование выполнения точного равенства (4) не имеет смысла и целесообразно искать приближенное решение, удовлетворяющее уравнению (4) с допустимой точностью. Кроме того, иногда интеграл (4) удается взять в общем виде. Тогда наиболее простым и логичным методом решения данной задачи является подбор параметров для аналитического решения интеграла (4) методом наименьших квадратов.

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

спектрометра, апериодическую функцию

к(0 = е““П(0 ,

(5)

1,

t > 0, t < 0,

т. е. единичная функция

Хевисайда; а — параметр, характеризующий инерционность, имеет размерность с-1.

Для аппаратной функции вида (5) и модели пика (3) решением интеграла (4) является уравнение

) = А

-,-а Л

а

t - t■^

V V

а

42

-Ф

а

42

(6)

/у

Рис. 4. Графики функции /(/) (6) (при условиях, указанных в тексте) и пика изотопа стронция 8г 88

Невязки, вычисленные для модели (6) носят периодический характер, что указывает на наличие колебательной составляющей в аппаратной функции. Следующая модель аппаратной функции, более точно описывающая динамические процессы внутри измерительной системы, — уравнение свободных затухающих колебаний

(7)

График функции (6) (при А = 37.7, а = 0.09, t\ = 41, ^ = 88, а = 7.5) и пик 8г88 изображены на рис. 4.

Степень близости модели к экспериментальным данным определялась по значениям невязок и по критерию X Пирсона [6].

Оценка среднего значения невязок для модели (3) и модели (6) не превышает 1.5 % от площади пика. Проверка соответствия моделей (3) и (6) экспериментальных данных по критерию X Пирсона показала ее правдоподобность с вероятностью 0.95. Для проверки моделей (3) и (6) экспериментальные данные регистрировались соответственно с интервалом времени 1.0 и 0.2 секунды.

где ю — круговая частота затухающих колебаний.

Аналитическое выражение свертки функций (3) и (7) получается довольно громоздким. Поэтому модельный пик проще вычислять непосредственно по формуле (4). Поскольку функция свертки (3) и (7) получается достаточно гладкой и без особенностей, нет необходимости в использовании сложных адаптивных алгоритмов интегрирования. Для вычислений можно воспользоваться одной из многоточечных прямых формул интегрирования, например Уэддля [7]. График функции (4) (при А = 780, ^ = 30, Ь = 78, л = 7, а = 0.22, ю = 0.050) и пик 8г88, снятый с интервалом между измерениями в 0.2 с, представлены на рис. 5. На этом же рисунке показан пик, восстановленный рассматриваемым методом. На рис. 6 изображены модели пика при (А = 90, 8 = 40, t2 = 86, /л = 7, а = 0.22, ю = = 0.048) и пик 8г84, снятый с интервалом между отсчетами в 0.275 с. Можно видеть, какое сильное влияние на снимаемые данные оказывает инерционность системы регистрации. Сравнение этой модели с экспериментальными данными показало,

Рис. 5. График функции (4), снятый при указанных в тексте условиях, и пик изотопа стронция 8г 88

Рис. 6. Модель пика, вычисленная при указанных в тексте условиях, и пик изотопа стронция 8г84

что среднее значение невязок не превышает 4.8 % от площади пика. Так как модели пика представляют собой функции, нелинейные относительно своих параметров, то для их оценивания использовался итерационный метод Ньютона [8, 9].

Результаты проверки работы алгоритмов

Проверка работы алгоритмов показала:

1. Практическую независимость от шумов точности оценки параметров пиков с использованием модели формы в диапазоне сигнал/шум от 50 и выше.

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

3. Погрешность восстановления центра пика не превышает минимального шага системы управления разверткой.

ОБНАРУЖЕНИЕ ПИКОВ И ОЦЕНКА ИХ ПАРАМЕТРОВ В РЕЖИМЕ С НЕПРЕРЫВНОЙ РАЗВЕРТКОЙ

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

Теория и возможности метода обнаружения и оценки параметров масс-спектрометрических пиков методом свертки экспериментальных данных с производными гауссовых функций описаны в статье [10]. В настоящей работе приведем основные формулы и покажем на двух примерах возможности метода.

Частные производные второго и четвертого порядка для гассовых функций представляются в виде следующих функций:

<$2(0=(1 - ^2/^о2)ехр(-^2/2^о2), (8)

54(0 = (3 - 6Х2/цо2 + ^4/^о2)ехр(-^2/2^о2).

Эти функции "подобны" пику в том смысле, что имеют экстремумы при X = о.

Теперь рассмотрим свертки «^(Х) и ^(Х) с отдельно взятым пиком:

У 2 (т) = - Л1 ехр У 4 (т) = - Л1 ехр

2^2

( X2 ^

2М2 ,

52(т- X )ах, (9)

54(т- X ^Х. (1о)

Взяв интегралы, получим:

У 2 (т) =

= Лцп12 /[1 + (ц/цо)2]32 х х (1 -т2/(ц2 +Цо2)) х

х ехр(-т2/(ц2 +Цо2)); (11)

у 4 СО =

= А^л12 /[1 + (и /Яо)2]52 х х (з - бт2/(и2 +и2))+о4/(и2 + Ао)2 х х ехр(-т2/(и2 +ио2))- (12)

Из формул (11) и (12) следует:

— При т = 0 у2(т) и у4(т) имеют абсолютный максимум.

— У2(0)/у4(0) = 1 + (и + А))2

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

Таким образом, производится обнаружение пика в спектре. Факт обнаружения пика в данном случае аналогичен действию фильтра с выделением заданной частоты [11]. Традиционный фильтр, выделяющий полосу, основан на вычислении свертки исходного сигнала с синусоидальной функцией. В связи с ортогональностью и единственностью функций на основе производных четных порядков {5'„(^)}, в базисе этой системы может быть представлен сигнал исследуемого масс-спектра, аналогично представлению сигнала в традиционном гармоническом базисе. Выполнение операций свертки при п = 2 и п = 4 дает возможность выделить пик в заданной точке оси масс, если он там существует, точно так же, как свертка с синусоидальной функцией дает возможность выделить сигнал заданной частоты.

Затем по значению сверток в максимуме, используя формулы (11) и (12), можно вычислить полуширину и и амплитуду А обнаруженного пика. Обозначим значение свертки у2(т) в максимуме как С2, а значение свертки у4(т) в максимуме как С4, тогда:

А=Ао ■

^2

А=

1

С _ ^

Vе* ,

3

2

1+А

Ао

(13)

2 (14)

2пи

Предлагаемый метод проверен с помощью компьютерной программы.

Примеры работы метода

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

Табл. 2. Возможности алгоритма обнаруживать и оценивать параметры одиночного пика при различных отношениях сигнал / шум (с/ш)

Отношение с/ш Значение амплитуды пика в мВ Погрешность оценки в %

100 1000 0

10 1004 0.4

5 1021 2

1 1070 7

0.5 1160 16

-0.2 ---------------------1------------------1------------------1-------------------1------------------1------------------1------------------1-------------------1------------------1------------------

-10 -8 -6 -4 -2 0 2 4 6 8 10

Рис. 7. Исходный сигнал с соотношением сигнал / шум 100

Свертка гаусоианы с четвертой производной

-ДО--------------------------1-----------------------1------------------------1-----------------------1-----------------------1-----------------------1-----------------------1-----------------------

-20 -15 -10 -5 0 5 10 15 20

Рис. 8. Свертка сигнала (рис. 7) с четвертой производной гауссианы

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

На рис. 9, 10 представлены аналогичные графики для отношения сигнал / шум 0.5.

Гауссовы функции, с производными которых производятся свертки, должны иметь полуширину примерно в 1.5-2 раза меньше, чем полуширина пика в сигнале масс-спектрометра.

ЗАКЛЮЧЕНИЕ

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

исходный сигнал

-Е 1 1 1 1 1 1 1 1 1 1

-10 -В -Б -4 -2 Q 2 4 6 8 10

Рис. 9. Исходный сигнал с соотношением сигнал/шум 0.5

Свертка гауссианы с четвертой производное

Рис. 9. Свертка сигнала (рис. 8) с четвертой производной

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

1. Манойлов В.В., Аракелянц М. А., Чернышев И.В. Программное обеспечение для определения изотопного состава аргона в автоматизированном комплексе на базе масс-спектрометра МИ1201ИГ // Научное приборостроение. 1999. Т. 9, № 4. С. 84-87.

2. Костоянов А.И., Манойлов В.В., Ефис Ю.М., Родионов М. В. Масс-спектрометрический комплекс для определения изотопного состава трудноионизируемых металлов // ПТЭ. 2000. № 1. С. 101-103.

3. Разников В.В., Разникова М. О. Информационно-аналитическая масс-спектрометрия. М.: Наука, 1992. 248 с.

4. Шубин В.М., Манойлов В.В. и др. Измерительно-вычислительная система для оценки атомных долей изотопов металлов на твердофазном масс-спектрометре МИ1201 с непрерывной разверткой // Научное приборостроение. 2000. Т. 10, № 2. С. 63-67.

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

6. Справочник по вероятностным расчетам. М.: Воениздат, 1970. 53б с.

7. Дьяконов В.П. Справочник по MathCAD PLUS 6.0 PRO. M.: СК Пресс, 1997. 336 с.

8. Вержбицкий В.М. Численные методы (линейная алгебра и нелинейные уравнения). М.: Высш. шк., 2000. 266 с.

9. Чебраков Ю.В. Теория оценивания параметров

в измерительных экспериментах. СПб.:

СПбГУ, 1997. 300 с.

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

11. Макс Ж. Теория и техника обработки сигналов при физических измерениях. М.: Мир, 1983. 280 с.

Институт аналитического приборостроения РАН, л/г ^ п/: опт

„ Л г Материал поступил в редакцию 16.06.2002.

Санкт-Петербург ^ ^

OUTLIER REJECTION AND PARAMETER ESTIMATION OF MASS SPECTROMETRIC SIGNALS FOR HIGH PRECISION ISOTOPIC ANALYSIS

V. V. Manoilov, I. V. Zarutsky

Institute for Analytical Instrumentation RAS, Saint-Petersburg

The paper presents algorithms of signal processing for gas and solid mass spectrometers applied to isotopic analysis of substances in geochronology and metal analysis in industry. The algorithms perform outlier rejection, peak center adjustment, single peak detection and parameter estimation.