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

ISSN 0868-5886 НАУЧНОЕ ПРИБОРОСТРОЕНИЕ, 2012, том 22, № 3, c. 30-35

= РАБОТЫ ДЛЯ МАСС-СПЕКТРОМЕТРИИ =

УДК 621.38

© В. В. Манойлов, Л. В. Новиков

ОЦЕНКА ПАРАМЕТРОВ МАСС-СПЕКТРОМЕТРИЧЕСКОГО

ПИКА В ДУБЛЕТЕ

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

Кл. сл.: обработка сигналов, метод независимых компонент, собственные векторы матриц, оценка параметров масс-спектрометрических пиков

ВВЕДЕНИЕ

В ряде практических приложений элементного, изотопного и других методов анализа требуется разделить компоненты дублета — смеси двух на-ложившихся пиков, амплитуда одного из которых в десятки раз превосходит второй. Оценка параметров таких пиков: в мультиплетах, и в частности в дублетах, основана на математических методах восстановления сигналов, искаженных аппаратной функцией. Среди этих методов, применяемых в приборостроении, наиболее известными являются метод Тихонова [1], метод максимального правдоподобия [2, 3], метод деконволюции [4, 5], метод на основе фильтров Винера [6], метод сингулярных разложений (SVD) [7], метод наименьших квадратов [8], метод сверток [9, 10], метод независимых компонент (ICA) [11] и др. С целью подавления шума во всех этих операциях в последнее время стали использовать избыточное вейвлет-преобразование с последующей пороговой обработкой [12]. В настоящей работе на основе идеи метода ICA предлагается метод разделения пиков известной формы в присутствии шумов.

ТЕОРИЯ

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

xj — Cs s2

(i)

где С — произвольная положительная константа, определяющая амплитуду главного пика, С > 10, и может быть значительно больше 10; s1, s2 — векторы-столбцы отсчетов размерности N главного (s1) и искомого (s2) пиков. Эти пики положительны, а суммарный вектор образует вектор исходных данных, как правило искаженный аддитивным шумом. Пусть известна (с некоторой точностью) форма первого пика ^, что при известном его положении позволяет сформировать вектор-столбец модели s0 = s1. Требуется, используя эту модель и исходный сигнал ж1, выделить из дублета компоненту s2. Решение, т. е. оценку в2, будем искать в виде

s2 — axj + bs0.

(2)

С учетом выражения (2), полагая что это — точное решение, образуем скалярные произведения исходных данных:

(s2,xj) — a(xjT • xj) + b(sj • xj) — kna + k12b, ( s2 , s0 )— a (xj • s0 )+ Ь ( sj • s0 ) — k2ja + k22b,

(3)

которые удобно представить в матричном виде как:

( T Л

sj • x

T

Vs2 • s0 J

( k К jj

V k2j

k„ Y a^

22 J

v b J

— K

fa^

VbJ

(4)

Построим двухстолбцовую матрицу x—(xj s0 ).

Используя выражение (3) легко убедиться в том, что элементы матрицы К образованы произведе-

нием xT • x . Введем далее вектор e1 =

С e >

V e21 J

ложим, что

(a\ V b J

= d

e

V e21 J

станта. Имея ввиду (1) и (2), положим a = de11 = 1.

Следовательно, b = de21 = e21/e11 и

a

v b J

1

e

V e21 J

K

a

VbJ

= ± K e

e

"11 V e21 J "11 Следовательно, для оценки S2 из (2) будем иметь

=^ л

e

e

V e 21 J

(5)

s2 = x1 +-

e^

-sn.

(6)

и

k22 Л

C (sT • so ) + (sT • S0 )

(sT • so )- Л

Величина AJ изменяется от A0 = — с шагом

и по-

где d — некоторая кон-

Используя последнее выражение, покажем, что вектор е1 является собственным вектором матрицы К, для которого должно выполняться равенство К е1 = Л1е1, где Л — собственное число. Действительно, продолжая равенства (4), имеем

-ДА до A0m" , при которой норма отрицательной составляющей оценки s2- достигает минимума, а именно min^ { ||s2-||/||s2||] . Величину шага удобно выбирать из условия ДА0 = а\|sj ^Ц/ЦхЦ при

а = 0.1,...,0.5, чтобы обеспечить плавное достижение искомого минимума независимо от амплитуды главного пика.

ОПИСАНИЕ АЛГОРИТМА

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

а) формирование вектора-столбца исходных данных x1;

б) формирование вектора-столбца модели s0;

в) формирование матрицы x = (x1 s0);

в) вычисление матрицы

K = ( XT • X ):

Чтобы определить ошибку оценки §2 и, следовательно, справедливость предположения (2), получим выражение для отношения е21/е11 . Из выражения (5) с учетом (4) будем иметь

^21 е11 + (к22 — Л ) е21 = 0

' t sT ^

х 1 • х 1 so • x1

T T

Vx1 • so so • so J

г) вычисление собственных векторов e =

= (e^ e2 ) =

(e e >

e11 12

V e21 e22 J

(и собственных значений

(7)

Вычисление Л показывает, что Л <(8^ • 80), поэтому

— = —С0, причем С0 > С . е11

Подставляя выражение (1) в (6), с учетом последнего вывода можно убедиться в том, что компенсация главного пика при вычислении оценки §2 происходит с избытком, т. е. в §2 на месте главного пика будет всегда присутствовать отрицательная составляющая.

С целью устранения этой компоненты вводится итерационная процедура

§2' = Х1 + А—1 So, У = 1,2,..., Ушах .

Л = (А1, Л2 ) ) матрицы К;

д) начало итерационного цикла у = 1 от

А0 = е21 /еи до Д0ш,х;

е) вычисление оценки искомого вектора

§у = х1 + 180;

2 1 0 0

ж) подавление шума в оценке §у;

з) вычисление на каждой итерации функциона-

ла s ( j) = Iis-7

2- / 2

и) вычисление минимума е по нескольким итерациям. В качестве условия минимума принимается выполнение неравенства

(s(j)-s(j-2)V2 < р ,

где Р — порог, величина которого определяется уровнем шума (например, Р = а / N + 5 • 10—2), а — среднеквадратические отклонения (СКО) шума;

к) если условие минимума не выполняется — переход к п. (д), иначе — к п. (л);

л) оценка положения и амплитуды пика § на итерации у — 2.

e

11

e

e

11

к

e

e

32

В. В. МАНОЙЛОВ, Л. В. НОВИКОВ

МОДЕЛИРОВАНИЕ

Проверка качества разделения компонент дублета по предлагаемому алгоритму выполнена в среде МАТЛАБ на двух моделях исходного сигнала = Сб1 + s2, где векторы s1 и s2 построены

а) по формуле ассиметричной гауссианы (задний фронт в 2 раза круче переднего);

б) в виде пика с плоской вершиной.

Модель пиков асимметричной гауссианы строилась по формуле:

^ ) = А ехр | - (V - tl )2/2»1),

(8)

где А1 — амплитуда, ^ — положение, цк — среднеквадратическая ширина пика, причем ¡к = ц1 при V < t1, ¡к = ц при V > t1 и ц> ц .

Из отсчетов функции 51 (V) формируется вектор s1, причем компонента этого вектора 51i = 51 (iДt), где интервал Дt выбирается из условия 50-60 отсчетов на пик. По формуле (2) формируется также вектор s2 с параметрами А2, t2, ц2к и строится модель исходных данных = Cs1 + s2 с добавлением белого шума с дисперсией а2 при отношении сигнал/шум, определяемом как А2 /а .

На рис. 1 показан исходный сигнал х1 (V) на сетке V = -100Дt ^ 100Дt при С=10, ¡=10, А=1 t1 = 0, ¡ц =у[2 ц; А2 = 1, t2 = 15, ц2 = -у/2 ц и отношении сигнал/шум 10. Из рисунка видно, что главный пик полностью скрывает искомый пик s2 (V) . Это обстоятельство, а также наличие шума усложняют решение задачи разделения пиков.

Построим модель главного пика дублета s0 (V) , используя формулу (8) с параметрами А0 = 1 = 0 ¡0 ц, сформируем матрицу х и вычислим матрицу К. Вычислим далее собственные векторы (и собственные значения) этой матрицы, используя команду eig и выполним операции в соответствии с пп. (д),..., (л) предлагаемого алгоритма.

Результат выделения искомого пика из дублета для приведенных выше исходных данных показан на рис. 1. Число итераций при этом оказалось равным 16. Итерационный процесс, который характеризуется зависимостью е (j), проиллюстрирован на рис. 2.

По результатам моделирования вычислялись среднеквадратические отклонения (СКО) и смещения оценок параметров искомого пика по М реализациям шума (М = 20):

Рис. 2. Кривая сходимости итерационного процесса е (1) = 11$1 у||$1 || для бигауссовой модели при отношении сигнал/шум 10

Рис. 1. Бигауссова модель дублета.

х1 (V) — исходные данные; (V) — искомый пик (штриховая линия); £2 (V) — оценка пика (сплошная линия)

СКО положения

ХОи - ?2 )7(м -!)

смещения положения ?2 -12

1/2

где

- -_1V" ?2 - М ^ 4>т '

1 у± т-1

СКО амплитуды

Г \1 (Ат - А )2/( М - 1)

2,т

2 I т-1

смещения амплитуды

где

_ 1 М л

А - М ^ ^ ;

^ ^ т-1

СКО восстановления формы

1

-VI|§7 -!

12

А-2 А-2 •

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

Алгоритм был апробирован при отношении сигнал/шум от 10 и более, при С -10 и 100. Результаты сведены в табл. 1.

Модели пиков с плоской вершиной строятся из двух пиков гауссовой формы с / - 4 Д?, расстояние между которыми такое, что образуется плоская вершина. Ширина таких пиков выбиралась равной 30 Д?. На рис. 3 показан исходный сигнал

х1 (?) из пиков с плоской вершиной на сетке ?--100Д * 100Д при С -10; А1 -1, - 0, А2 -1, ?2 -10Д?, ширине пиков 30 Д?. Фронты пиков — гауссовой формы при / - 4Д? и отношении сигнал/шум 10. Из рисунка видно, что главный пик полностью скрывает искомый пик s2 (?) . Применение предлагаемой вычислительной процедуры позволило решить задачу оценки параметров скрытого в дублете пика s2 (?) , численные результаты которой приведены в табл. 2. Результат выделения искомого пика из дублета для приведенных выше исходных данных проиллюстрирован на рис. 3.

т

Табл. 1. Погрешность оценки параметров пика бигауссовой формы в дублете при и - 10Д, А1 -1, - 0, /1 - 72 и; А2 -1, ?2 -15 Д, /2 -л/2/

Превышение главного пика Отношение сигнал/шум Погрешность оценки положения Д? Погрешность оценки амплитуды в % СКО восстановления

СКО Смещение СКО Смещение

С=10 10 1.4 -1.65 7.2 -0.5 0.13

50 0.65 -1.8 3.2 -1 0.05

Без шума 0 -2 0 -1 0.004

С=100 10 1.7 -1.4 9 0 0.15

50 0.65 -1.8 3.2 -1 0.05

Без шума 0 -2 0 -1 0.004

34

В. В. МАНОИЛОВ, Л. В. НОВИКОВ

*1(0, ^2(0

1.2 1

0.8 0.6 0.4 0.2 0

МО

(Л -/-

----------------5 2 (/) Г 11 1

11 || 11 и

II 1 1 1\ \ • \ V.

Рис. 3. Модель дублета с плоской вершиной: х1 ) — исходные данные;

s2 (/) — искомый пик (штриховая линия); ) — оценка пика

(сплошная линия)

-50

50

100

Табл. 2. Погрешность оценки параметров пика плоской формы в дублете при А1 = 1, ? = 0, А2 = 1, ?2 = 10А?.

Ширина пиков — 30 А, фронты пиков — гауссовой формы при ¡л = 4Аt

Превышение главного пика Погрешность оценки положения А? Погрешность оценки амплитуды в % ■ « о

* ^ э! оа СКО Смещение СКО Смещение СКО восстан ления

£ £ о 5

С=10 10 2.1 0.55 6 6 0.118

50 1 0 3.2 1.2 0.035

Без шума 0 0 0 -1 0.015

С=100 10 2.05 0.6 7 7 0.13

50 1.1 0 3.5 1.2 0.04

Без шума 0 0 0 1.2 0.02

При моделировании изменялись положения главного и искомого пиков от 5 А? до 35 А? для бигаусса и от 10А? до 35А? для пика с плоской вершиной, т. е. от полного наложения до практически полного разделения. При максимальном разделении главного и исходного пиков для обеих моделей погрешности оценки положения и отклонения формы практически не изменились. Погрешность оценки амплитуды при этом уменьшилась на порядок.

ЗАКЛЮЧЕНИЕ

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

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

Авторы благодарят Тимченко Н.А. за активное участие в обсуждении результатов работы.

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

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

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

3. Gelfgat V.I., Kosarev E.L., Podolyak E.R. Program for signal recovery from noisy using the maximum likelihood principle // Comput. Phys. Commun. 1993. V. 74. P. 335.

4. Погуляй А.В., Васильев Ю.В., Туймедов Г.М., Мазу-нов В.А. Повышение разрешения спектрометра методами деконволюции // Химия и компьютерное моделирование. Бутлеровские сообщения. 1999. № 2. Приложение: Материалы VI Всероссийской конференции "Структура и динамика молекулярных систем".

5. Погуляй А.В., Аблазимов Р.Р., Туймедов Г.М., Мазу-нов В.А. Сравнительный анализ методов деконво-люции // Химия и компьютерное моделирование. Бутлеровские сообщения. 2001. № 6. Приложение: Материалы VIII Всероссийской конференции "Структура и динамика молекулярных систем".

6. Бардин Б.В., Белов В.Д., Новиков Л.В., Чижов Ю.В. Использование оптимального фильтра Винера для деконволюции электронных спектров // Научное приборостроение. 1999. Т. 9, № 1. С. 53-59.

7. Shrager R.I. Chemical transition measured by spectra

resolved using singular value decomposition // Che-mom. Intell. Lab. Syst. 1986. V. 1, N 1. P. 59-70.

8. Bardin B.V., Belov V.D., Mamro N.V. et al. A Computerized control and data handling system of a multifunctional electron spectrometer and software for processing complex electron spectra // Instruments and Experimental Techniques. 1999. V. 42, N 2. P. 205211.

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

10. Манойлов В.В., Заруцкий И.В. Возможности алгоритма сверток с производными для оценки параметров масс-спектров, содержащих наложившиеся пики // Научное приборостроение. 2009. Т. 19, № 4. С. 103-109.

11. Hyvarinen Aapo. Survey on independent component analysis // Neural Computing Surveys. 1999. N 2. P. 94-128. URL: (http ://www.icsi.berkeley.edu).

12. Князева Т.Н., Новиков Л.В., Орешко Н.И. Удаление нестационарного шума из экспериментальных данных // Научное приборостроение. 2008. Т. 18, № 2. С. 61-65.

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

Контакты: Манойлов Владимир Владимирович, [email protected]

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

PARAMETER ESTIMATION FOR THE MASS SPECTROMETRIC

PEAK IN THE DOUBLET

V. V. Manoilov, L. V. Novikov

Institute for Analytical Instrumentation of RAS, Saint-Petersburg

We propose a method of estimating the position and the intensity of the peak in the doublet, mostly hidden peak, whose intensity is greater than the required 10 or more times in the presence of noise. The position of the main peak at exactly this set, its shape is known with some accuracy, but the amplitude is unknown. The method allows to estimate the parameters of the hidden peak at a signal / noise ratio of ten and above, and the distance between the peaks of the half-width at half height, and more.

Keywords: signal processing, method of independent components, eigenvectors of matrices, estimation of parameters of mass spectrometric peaks