CC BY
124
УДК 539.107
С. В. Молчанов, О. В. Рубан, И. Г. Мершиев,
Г. В. Мозжухин, Г. С. Куприянова
ПРИМЕНЕНИЕ АДАПТИВНОЙ ВЕЙВЛЕТ-ФИЛЬТРАЦИИ ДЛЯ ДЕТЕКТИРОВАНИЯ ЗАШУМЛЕННЫХ СИГНАЛОВ ЯДЕРНОГО КВАДРУПОЛЬНОГО РЕЗОНАНСА
Представлен новый подход к обработке зашумленного сигнала ядерного квадрупольного резонанса, содержащего радиочастотные помехи, с низким отношением сигнал/шум. Подход основан на предварительном статистическом анализе шумового сигнала по данным коэффициентов вейвлет-разложения и использовании адаптивно-порогой техники. Важным является автоматизированный выбор порога на основе параметров, характеризующих статистические свойства шума.
The article presents a new approach to signal processing with low signal-to-noise ratio. The approach is based on preliminary statistical analysis of the noise based on wavelet transform coefficients and the adaptive and threshold methods. Special attention is paid to the automatic detection of the threshold on the basis of statistical noise properties parameters.
Ключевые слова: ядерный квадрупольный резонанс, вейвлет, шум.
Keywords: nuclear quadruple resonance, wavelet, noise.
Введение
В последнее время одним из перспективных направлений прикладных исследований является разработка физических методов спектроскопии ядерного квадрупольного резонанса (ЯКР) и ядерного магнитного резонансов (ЯМР) для детектирования квадрупольных ядер для целей идентификации запрещенных веществ, исследования лекарственных препаратов, определения внутренних напряжений в материалах [1; 2; 3]. Результат детектирования зависит от множества факторов. Наиболее важными являются: частотный диапазон спектра ЯКР детектируемого вещества, влияющий на интенсивность сигнала ЯКР, примеси в исследуемом материале, вызывающие искажения сигнала, передвижения датчика, а также помехи от различных источников внешних сигналов. Ядерный квадру-польный резонанс ядер азота 14N широко применяется для детектирования запрещенных веществ и лекарственных препаратов. Частотный диапазон спектра ЯКР азотсодержащих соединений лежит в области 0,1—5 МГц, что приводит к относительно слабой интенсивности сигналов ЯКР.
Одной из проблем, возникающих при обработке сигналов ЯКР 14N, является повышение соотношения сигнал/шум и выделение слабого сигнала на фоне сильных внешних помех. Традиционно с этой целью исполь-
Вестник Российского государственного университета им. И. Канта. 2009. Вып. 4. С. 71-80.
зуют метод многократного накопления зашумленных сигналов, который теоретически для случая с белым шумом позволяет увеличить отношение сигнал/шум в -/м раз, где N — число накоплений [4]. Но в последнее время для регистрации сигналов ЯКР стали применяться высокодобротные датчики со значениями добротности Q > 10000 [5], что приводит к резкому сокращению полосы пропускания и, следовательно, к появлению коррелированных составляющих шума и к нестационарному поведению сигнала. Использование традиционных методов фильтрации не позволяет решить две наиболее важные проблемы, с которыми приходится сталкиваться при распознавании полезного сигнала, — это подавление паразитных сигналов и внешних сигналов, источники которых имеют основную частоту, близкую к частоте ЯКР, и увеличение отношения сигнал/шум в присутствии коррелированных шумов. Применением двухканальной системы в спектрометре ЯКР, один из каналов которого регистрирует только внешнюю помеху, удалось добиться подавления внешней помехи и повысить эффективность применения процесса накопления [6]. Это позволило идентифицировать сигнал от 150-граммового образца на расстоянии 15 см при 3000 накоплениях в отсутствие экранирования приемной системы. Значительного улучшения отношения сигнал/шум удалось достичь также за счет использования методики очищения сигнала от шума, основанной на вейвлет-фильтрации с применением алгоритма Кадзофа [5; 7]. Однако анализ большого числа экспериментальных данных показал, что эта методика не позволяет устранить паразитные сигналы, возникающие в режиме высокой добротности. Вейвлет-фильтр не отличает коррелированную помеху от полезного сигнала. А процесс накопления приводит к усилению и помехи и сигнала, что создает значительные трудности при идентификации полезного сигнала.
В данной работе предлагается новый подход, основанный на преимуществах вейвлет-анализа зашумленного сигнала с помехой и использовании адаптивно-пороговой техники, дающей возможность установить порог на основе параметров, характеризующих статистические свойства сигнала с шумом. Адаптивные устройства обработки данных отличаются наличием определенной связи параметров передаточной функции с параметрами входных, выходных и прочих дополнительных сигналов. Эти свойства адаптивных систем позволяют приемной системе «самой» настраиваться на оптимальную обработку сигналов. В основном адаптивная фильтрация применяется для очистки данных от нестабильных мешающих сигналов и шумов, перекрывающихся по спектру со спектром полезных сигналов. Также адаптивная фильтрация может применяться, когда полоса мешающих частот неизвестна, переменна и не может быть задана заранее для расчета параметрических фильтров. Эти свойства делают адаптивные системы наиболее подходящими для обработки сигналов с малым отношением сигнал/шум [8]. Однако результаты применения адаптивной фильтрации для анализа сигналов в значительной мере зависят как от отношения сигнал/ шум, так и от свойств самого сигнала. Кроме того, наиболее слабым местом является выбор порога при необходимости устранения помехи и диагностики полезного сигнала. Поэтому в
данной работе уделено основное внимание исследованию статистических характеристик реальных шумовых сигналов и сигналов ЯКР в приемном тракте спектрометра после цифровой регистрации.
Статистический анализ шумового сигнала
С целью выяснения свойств и особенностей ЯКР сигнала, регистрируемого в режиме с высокой добротностью, были изучены статистические свойства шума. Как правило, считается, что шум является белым гауссовым. Поэтому в данной работе проводится сравнение статистических свойств шумового сигнала со свойствами модельного белого шума. Экспериментальный шум исследовался при комнатной температуре на спектрометре ядерного квадрупольного резонанса на базе консоли ЯМР/ЯКР Аполло производства компании Tecmag (США). Несущая частота выбиралась равной частоте резонанса ЯКР 4641 кГц на ядрах азота 14М в нитрате натрия МаМО2. Датчик не содержал образца. Полоса пропускания цифрового приемника составляла 200 кГц. Использовалось квадратурное детектирование. При сборе данных исследуемый сигнал был оцифрован 1024 точками с интервалом между ними 10 мкс. Для изучения статистических свойств проводилось разложение шумового сигнала с помощью вейвлет-функций и исследовалось поведение коэффициентов разложения в разных полосах частот. Оценивались следующие статистические параметры: среднее значение Мх = (X XI)/п, дисперсия Бх = (2(х1 - Мх)2)/п, асимметрия Бх = (Х((х - Мх)/о)3)/п и эксцесс Кх = (2((х1 - Мх)/о)4)п - 3). Здесь х1 выборка значений реального шумового сигнала, о = (Бх)1/2 — стандартное отклонение. Выбиралась также — последовательность такой же длины модельного шумового сигнала, который предполагался белым гауссовым шумом. Результаты статистического анализа реального шума сравнивались с поведением модельного шума. Для обоих сигналов проведена нормировка и центрирование по формуле: хпогт = (х-Мх)/тах( |х-Мх |).
В качестве базиса разложения использовались функции Добеши второго порядка и пять уровней разложения ) [9; 10]. Для каждой )-й последовательности (то есть для каждой полосы разложения) коэффициентов разложения реального шума ут) и белого шума ( = 1...5) вычислялись статистические параметры. Псевдочастоты Б, соответствующие полосам, рассчитывались с помощью нахождения центральной частоты вейвлета. Результаты вычислений дисперсии, среднего значения, асимметрии и эксцесса вейвлет-коэффициентов в зависимости от полосы частот представлены на рисунке 1.
Рис. 1. Зависимости дисперсии, среднего значения, асимметрии и эксцесса вейвлет-коэффициентов от полосы частот для модельного белого шума и реального экспериментального шума, регистрируемого в ЯКР спектрометре
В большинстве методов выделения сигнала из шума предполагается, что шум является гауссовым. Поэтому важным вопросом, который решался в данной работе, — это выяснить, в какой мере шумовой сигнал ЯКР можно рассматривать как гауссовый шум. Изучив статистику реально наблюдаемого шума с помощью вейвлет-преобразования, мы убедились в том, что такое предположение в данном случае неверно.
Анализ показывает, что с уменьшением полосы наблюдения расхождения в значениях дисперсия и средние значения для белого и реального шума возрастают. Необходимо отметить, что известные алгоритмы шумоподавления основаны на сокращении количества шумовых коэффициентов в зависимости от некоторого порога, значение которого выбирается в зависимости от величины дисперсии коэффициентов на каждом уровне разложения. Например, в случае жесткого шумоподавления сохраняются неизменными все коэффициенты уровня ], большие или равные порогу, который вычисляется по формуле: г = гп = Су]21о%п, где с — дисперсия шума, а остальные обращаются в
ноль. Большой разброс значений дисперсии для реального шума для различных полос частоты не позволяет применять традиционные пороговые функции для оптимального шумоподавления. В этом случае необходимо учитывать значение дисперсии вейвлет-коэффициентов для каждого масштаба вейвлет-преобразования, повышая, таким образом, эффективность обработки зашумленных сигналов.
Кроме того, анализ показал значительное отклонение асимметрии и эксцесса от характеристик модельного белого шума. Поэтому в данной работе была проведена также проверка распределения на нормаль-
ность с помощью критерия Шапиро — Уилка [11]. Этот критерий наиболее подходит в данном случае: он используется для относительно небольших выборок данных (до 4000 отсчетов), хорошо подходит для асимметричных распределений. Критерий Шапиро — Уилка проверяет нулевую гипотезу о том, что выборка Хі...Хп распределена нормально. Тестовая статистика рассчитывалась по формуле
W = ■
2 а,х.
(і)
2 (Хі -< X >) і=1
Здесь Х(1) — статистика 1-го порядка, то есть 1-е наименьшее число в выборке; <х> — выборочное среднее; определяются выражением
(а1,~, ап) =
шт¥-1
4штУ 1У 1ш
где т = (т^.. тп)Т; Ш!... тп — ожидаемые значения порядковых статистик независимых и одинаково распределенных величин, относящихся к нормальному распределению, а V — ковариационная матрица этих порядковых статистик. Чем ближе значение Ш к единице, тем больше вероятность того, что нулевая гипотеза верна (рис. 2).
Рис. 2. Статистика распределения белого и реального шума и оценка нормальности с помощью критерия Шапиро — Уилка На рисунке 2 представлен график «квантиль-квантиль», показывающий отклонение наблюдаемого распределения от нормального. Здесь Ъ — квантиль нормального распределения — число стандартных отклонений от среднего в точке, которая должна принадлежать нормальному распределению с соответствующим средним и стандартным отклонением. Если выборка распределена нормально, точки образуют прямую линию. Отклонение от прямой указывает на ненормальное распределение.
В данном случае плотность распределения наблюдаемого шума существенно отличается от нормального закона и имеет длинный «хвост». Кроме того, полученные результаты коррелируют с гисто-
2
і=1
граммой распределения реального шумового сигнала, представленной на рисунке 3.
Хпогш
Рис. 3. Гистограмма распределения экспериментального шумового сигнала и распределение белого шума.
Статистический анализ свойств экспериментального шумового сигнала позволяет применить нам два подхода для повышения эффективности методики очищения сигнала от шума и повышения точности диагностики полезного сигнала.
Первый подход основан на применении фильтра, который строится на основе статистических свойств исследуемого сигнала. В данном случае целесообразно применять медианный фильтр, который эффективно ликвидирует «хвосты» в распределении (рис. 2) и преобразует форму распределения к гауссовой. С целью подбора размера окна для медианного фильтра было исследовано, как меняются асимметрия Бх и эксцесс Кх распределения в зависимости от окна медианного фильтра а (рис. 4). Здесь а — число точек, по которым производится усреднение.
Основываясь на том, что асимметрия и эксцесс для нормального распределения равны нулю, размер окна выбран равным а=5. Результаты применения медианного фильтра с окном а=5 представлены на рисунке 5. На рисунке показано, как выглядит гистограмма распределения и график «квантиль-квантиль» в результате обработки шумового сигнала оптимальным медианным фильтром. Видно, что распределение теперь близко к нормальному. Коэффициент Ш критерия Шапиро
— Уилкса близок к 1.
Рис. 4. Зависимость асимметрии и эксцесса от размера окна для экспериментального шумового сигнала
700
■0.45 -0.4 -0.35 -03 -0.25 -0.2 -0.15 -0.1 -0Ю5 0 0.05 0.1 0.15 0.2 0.25 0.3 035 0.4 0.45
хгюгш
а
Рис. 5. Гистограмма распределения экспериментального шумового сигнала и распределение белого шума (а) и статистика распределения белого и реального шума после применения медианного фильтра (б)
Применение медианного фильтра для предварительной обработки шумового сигнала позволяет более эффективно провести накопление
сигнала и повысить отношение сигнал/шум. Однако в ряде случаев выбор размера окна оказывается проблематичным в силу резких колебаний кривых эксцесса и асимметрии. Отсутствие ярко выраженного минимума этих кривых приводит к неоднозначности процесса выбора параметра а и не всегда удается оптимизировать этот процесс и добиться удовлетворительного результата. Более эффективным оказался второй подход, основанный на адаптивной фильтрации сигнала. Отличие данного метода от метода очищения сигнала от шума, который предлагается в пакетах MATLAB, заключается в том, что в зависимости от хода кривой дисперсии в каждой частотной полосе (рис. 1) выбирается свой порог. На рисунке 6 показана предлагаемая блок-схема адаптивной вейвлет-фильтрации. Второй подход фильтрации сигнала выполняется следующим образом:
— сбор данных образцового сигнала, разложение (декомпозиция) сигнала по базисным вейвлет-функциям;
— вычисление порогой функции для каждой полосы разложения в зависимости от статистических свойств сигнала по данным дисперсии в
соответствии с формулой Т^п, j = Ujyl2 log n , где aj = Dxj — дисперсия
вейвлет-коэффициентов в j-й полосе;
— пороговая обработка коэффициентов;
— обратное вейвлет-преобразование.
принятым
сигнал
1 1
пряиое адаптивная обратное
вейвлет пороговая вейвлет-
преобразование п обработка п преобразование
выходном
сигнал
Рис. б. Блок-схема адаптивной вейвлет-обработки сигнала Результаты и обсуждение
Изучение статистических свойств шумовых сигналов, регистрируемых в ЯКР спектрометре, датчик которого настроен на резонансную частоту, но не содержал образца, показало, что коррелированные шумы приводят к тому, что свойства шумового сигнала сильно отклоняются от свойств белого шума. Это приводит к тому, что накопление не является эффективным. Разработанный медианный фильтр, параметры окна которого подбираются в зависимости от поведения эксцесса и симметрии, позволяет во многом улучшить характеристики реального шума и увеличить отношение сигнал/шум при накоплении сигнала. Исследование большого числа сигналов показало, что более эффективна методика обработки сигналов с помехой — это адаптивная вейвлет-фильтрация. Данный подход был применен для обработки сигнала ЯКР МаМ02 на частоте 4641 кГц в присутствии внешней помехи при добротности контура 0=200. После предварительного анализа сигнала, полученного без образца, данные вводились в разработанную программу для вейвлет-разложения и вычисления дисперсии и среднего.
При обработке образцового сигнала эти данные позволяют автоматизировать процесс выбора порога в каждой полосе разложения. Пример обработки зашумленного сигнала с помощью математического пакета МЛТЬЛБ и предлагаемой методики представлен на рисунке 7.
Аотн. ^
■0.8 '-------1-------1-------1--------1-------1-------1-------1--------1-------1-------1
0 100 200 300 400 500 600 700 800 900 1000
X
а
б
Рис. 7. Результат обработки зашумленного сигнала ЯКР с помощью: а — пакетного вейвлет-преобразования; б — адаптивной вейвлет-фильтрации, при которой в зависимости от оценки шумовой дисперсии для каждого масштаба разложения устанавливается свой порог шумоподавления
Список литературы
1. Latosinska J. N. Nuclear Quadrupole Resonance spectroscopy in studies of biologically active molecular systems — a review / / Journal of Pharmaceutical and Biomedical Analysis. 2005. З8. Р. 577—587.
2. Miller J. B., Barrall G. A. Explosives Detection with Nuclear Quadrupole Resonance // American scientist. 2005. Vol. 93 (1). Р. 50.
3. King J. D., A. De Los Santos. Development and evaluation of magnetic resonance technologies, particular NMR, for detection of explosives. Appl. Magn. Reson. 2004. 25. Р. 535—565.
4. Cancino-De-Greiff H. F., Ramos-Garcia R., Lorenzo-Ginori J. V. Concepts in Magnetic Resonance. 2002. Vol. 14(6). Р. 388—401.
5. He D. F., Tachiki M. Itozaki H. 14N nuclear quadrupole resonance of p-nitrotoluene using a high-Tc rf SQUID, Supercond. Sci. Technol. 2007. 20. Р. 232 — 234.
6. Мозжухин Г. В., Куприянова Г. С. и др. Детектирование сигналов импульсного квадрупольного резонанса в условиях сильных помех // Вестник Российского государственного университета им. И. Канта. 2007. Вып. 3. С. 54.
7. Куприянова Г. С., Мозжухин Г. В. и др. Методы обработки сигналов магнитного резонанса в системах диагностики материалов в условиях технологических помех // Магнитный резонанса и его приложения. 5-я зимняя молодежная школа-конференция. 1 — 5 декабря 2008. Санкт-Петербург. 2008. А-183.
8. Chang S. G., Yu B. et al. Spatially adaptive wavelet thresholding with context modeling for image denoising // IEEE Trans. ImageProcessing. Sept. 2000. Vol. 9. P. 1522 — 1531.
9. Дьяконов В. П. MATLAB 6.5 SP1/7 +Simulink 5/6. Обработка сигналов и проектирование фильтров. М., 2005.
10. Chipman H. A., Kolaczyk E. D., McCulloch R. E. Adaptive Bayesian wavelet shrinkage // J. Amer. Statist. Assoc. 1997. Vol. 92, № 440. P. 1413 — 1421.
11. Кобзарь Ф. И. Прикладная математическая статистика. М., 2006.
Сведения об авторах
Г. В. Мозжухин — канд. физ.-мат. наук, доц., РГУ им. И. Канта, mgeorge@gazinter.net
С. В. Молчанов — ст. преп., РГУ им. И. Канта.
Г. С. Куприянова — д-р физ.-мат. наук, проф., РГУ им. И. Канта.
Authors
G. Mozzhukhin — Dr., IKSUR, mgeorge@gazinter.net.
S. Molchanov — IKSUR.
G. Kupriyanova — Prof., IKSUR.
