CC BY
99
УДК 539.107
С. В. Молчанов, Г. В. Мозжухин
ДЕТЕКТИРОВАНИЕ СИГНАЛОВ ЯДЕРНОГО КВАДРУПОЛЬНОГО РЕЗОНАНСА С ОГРАНИЧЕННОЙ ВЫБОРКОЙ ДАННЫХ
Рассмотрен метод выделения сигналов импульсного ядерного квадрупольного резонанса (ЯКР) из шумов в условиях малого отношения сигнал/шум в детекторах азотсодержащих соединений. Для анализа экспериментальных данных сигналов ЯКР использовалась авторегрессионная модель с применением метода максимальной энтропии. Полученные данные показывают эффективность этой методики в условиях ограничений времени детектирования и цифровой выборки сигнала в детекторах ЯКР.
The denoising method of the pulse nuclear quadrupole resonance (NQR) signals was investigated for the case of the small signal to noise ratio. The maximum entropy method of the autoregressive model was applied for the analysis of the experimental NQR data. The obtained results show the efficiency of proposed method in the conditions of the limited sampling rate for the limited time of the detection and digital signal sampling for the NQR detectors.
Ключевые слова: ядерный квадрупольный резонанс, авторегрессионная модель, метод максимальной энтропии.
Key words: nuclear quadrupolar resonance, autoregressive model, maximum entropy method.
Введение
В последние годы получило развитие применение метода ядерного квадрупольного резонанса для детектирования запрещенных азотсодержащих веществ, в том числе и взрывчатых, в сканерах багажа и детекторах мин [1]. Основной способ детектирования состоит в применении импульсного ЯКР с последующей цифровой обработкой полученных сигналов. По сравнению с обычными измерениями в импульсном методе ЯКР образец представлен незначительным количеством вещества (10 — 20 г) в большом объеме резонансной катушки (до 160 л) [2]. Кроме того, зачастую устройство находится в среде с большим количеством радиопомех, например в аэропортах. Для взрывчатых веществ характерные частоты ЯКР располагаются в диапазоне от 0,4 до 5,3 МГц. Причем для наиболее распространенного взрывчатого вещества, такого как тринитротолуол, частота обнаружения равна 842 кГц, а частота детектирования наиболее доступного для террористов нитрата аммония — 423 кГц. Как известно, в импульсном методе ЯКР интенсивность сигнала пропорциональна квадрату частоты детектирования сигнала ~ v2 [3]. Таким образом, обнаружение взрывчатых веществ по частоте ЯКР представляет сложную задачу в основном из-за слабого сигнала ЯКР и плохого отношения «сигнал/шум», обычно много меньшего І.
В качестве основного импульсного метода для регистрации сигналов ЯКР в сканерах багажа применяются многоимпульсные последовательности радиочастотных импульсов. В окнах наблюдения между радиочастотными импульсами регистрируются сигналы ЯКР в форме сигналов свободной индукции (ССИ) и спинового эха (ССЭ) в зависимости от выбранной модификации многоимпульсной последовательности. Одни из основных требований к регистрации сигналов в ЯКР сканерах — ограниченное время регистрации и необходимость отделения стабильной внешней помехи от сигнала ЯКР. Для увеличения количества зарегистрированных окон наблюдения используют минимальное время между импульсами в серии порядка 0,5 — 2 мс. При использовании цифровых методов регистрации время дискретизации составляет At = 1 - 20 мкс для используемых типов амтілигудно-цифровых преобразователей. Спектральная полоса цифровых приемников определяется как В -1/(2At) [4].
Для сокращения полосы приемника применяют большие значения времен дискретизации, что помимо сокращения полосы приемника может вести к значительным колебаниям амплитуды сигнала [5]. Наименьшая частота Av, которая может быть разрешена в спектре в импульсной фурье-спектроскоггии определяется соотношением Av<l/Ta, где Та — длительность зарегистрированного временного сигнала [6]. Стремление уменьшить временной интервал для
регистрации сигнала и сузить полосу пропускания для уменьшения входных шумов в цифровом приемнике приводит к тому, что регистрация сигнала происходит при условии, близком Ду = 1/Та.
При наблюдении сигналов свободной индукции, спинового эха с низким отношением «сигнал/шум» в ЯКР экспериментаторы используют традиционный метод выделения сигналов из шума — временное накопление. Эффективность этого метода зависит от когерентности наблюдаемых выборок данных, особенно при большом количестве накоплений N и пропорционально в случае некогерентного накопления получаем иную зависимость, так как случайные изменения фазы сигнала приводят к проигрышу по мощности в Q = exp(o2) раз, где о2 — дисперсия фазы [7]. В этих случаях в качестве оптимального обнаружителя для детерминированных и квазидетерминированных сигналов применяется квадратурный детектор, построенный по двухканальной схеме, что приводит к децимации выборки в два раза на выходе детектора. Это приводит в большинстве случаев к неудовлетворительным результатам применения фурье-анализа и непараметрических методов оценивания спектральной плотности мощности сигналов, так как разрешающая способность этих методов ДF зависит от длительности выборки Т = Foкнo/ ЛF, где F0кн0 — оконная функция [8]. Для преодоления указанных недостатков обычно используют параметрические методы, позволяющие выбрать функцию, аппроксимирующую реальный набор данных. Применение параметрических методов предполагает знание априорной информации, достаточной для построения параметрической модели исследуемого процесса. Если такая информация отсутствует, то применение параметрических методов может быть связано со значительными ошибками, возникающими из-за неадекватного выбора модели.
Авторегрессионные методы [9] позволяют улучшать или сохранять высокое разрешение без значительного ухудшения устойчивости спектральных оценок по сравнению с непараметрическими методами спектрального оценивания, недостатки которых обусловлены искажающим воздействием просачивания по боковым лепесткам из-за неизбежного взвешивания в них конечных последовательностей данных, а использование окна позволяет ослабить влияние боковых лепестков, но лишь за счет ухудшения спектрального разрешения. Таким образом, существующая методика регистрации сигналов ЯКР в детекторах-сканерах азотсодержащих соединений происходит в условиях ограниченной выборки точек окна наблюдения, что снижает спектральное разрешение метода и тем самым понижает достоверность обнаружения вещества. В данной работе исследуется эффективность использования параметрических методов для определения частоты сигнала ЯКР и возможность повышения частотного разрешения при ограниченной выборке.
Методы обработки данных. Авторегрессионная модель
В авторегрессионных методах при построении АР-модели данных, которая устанавливает функциональную связь между данными выборки сигнала, предпринимается попытка разложения данных на компоненты с помощью затухающих сигналов (метод Прони) [7]. Это улучшает временную локализацию и соответственно разрешение при оценке спектральной плотности мощности. В АР-модели выборка сигнала, как правило, представляется временным рядом и составляющей белого шума:
р
х(п) = - X а(к)х(п - к) + е(п), (1)
к=1
где p — порядок модели, e(n) — составляющая белого шума или ошибка, a(k) — коэффициенты временного ряда. Уравнение (1) можно представить в другом виде, используя оператор задержки z-k, тогда
х(п) = - X а(к)г~кх(п) + е(п), (2)
к=1
которое преобразуется к виду
;(п> . (з)
1+ Т. а(к)г к
ы\
Последнее уравнение эквивалентно z-преoбразoванию цифрового фильтра с бесконечной импульсной характеристикой. Тогда (3) можно представить как
где H(z) = 1
Р V
1+ Е a(k)z
к=1
Частотную характеристику фильтра H(o) получим, используя z = eimt, тогда
Н(со) =-----------------. (5)
1+ I а(к)е~]М к=1
Таким образом, ЛР-модель состоит из линейного рекурсивного фильтра H(a) порядка p с поданным на вход белым шумом e(n). Тогда спектральная плотность мощности полученного процесса будет равна произведению квадрата модуля частотной характеристики этого фильтра на спектральную плотность мощности белого шума
Р(а>) = \Н(со)\2 Pnolse(co) = \Щсо)\2 *2nolse(n) = ^ок>) . (6)
к=1
Здесь учтено, что (0}) = апыэе > причем это значение можно выразить через параметры
модели и, таким образом, определить спектральную плотность мощности. При расчете параметров АР-модели необходимо минимизировать ошибку е(и). Выражая среднеквадратическую ошибку с помощью соотношения
Е(п)=^ Уе2(п), (7)
N „
решаем задачу минимизации E(n), воспользовавшись уравнением
р
е(п) = х(п) + X а(к)х(п -к). (8)
к=1
Это приводит к решению уравнений Юла-Уокера, которые в матричной форме можно представить в следующем виде:
а(к) = -Етхх(к-тхх(к), (9)
где а(к) — вектор параметров модели при 1 <к<р; н!^(к - /) — автокорреляционная
симметричная матрица Теплица; Яхх (к) — вектор значений корреляционной функции Яхх (к) для к, 1 < к < р. Основная сложность при решении уравнении (9) состоит в оценке корреляционных функций, получаемых путем временного усреднения и учета краевых эффектов при построении корреляционной матрицы, которые выполняются следующими методами: автокорреляционным Левинсона — Дурбина, ковариационным, модифицированным ковариационным и методом Берга. Каждый из этих методов обладает достоинствами и недостатками, которые проявляются при спектральной оценке реальных сигналов. Нахождение параметров модели a(k) из (9) после подстановки в (6), учитывая, что Е(п)= сг^ж позволяет получить авторегрессионный спектр плотности мощности. Вычисление E(n) происходит по имеющейся выборке x(n), где п изменяется от 1 до N в случае использования непараметрических методов это эквивалентно оконному взвешиванию данных, приводящее к спектральному сглаживанию и, как следствие, к понижению разрешения. Для авторегрессионного метода это условие не является ограничением для спектрального разрешения, что позволяет эффективно использовать метод для обработки ограниченного набора данных. Для оптимальной аппроксимации данных в АР-модели необходимо выбрать порядок модели, позволяющий устранить появление ложных пиков и смещение частотных максимумов, с помощью критериев, предложенных Экейком. В случае ограниченного набора данных применяется информационный критерий Экейка:
А1С(р) = ЫЫЕ(р) + 2р. (10)
Развитием авторегрессионных методов стало применение принципа максимальной энтропии Берга (ММЭ) [9], при этом выбранный временной ряд должен отличаться максимальной хаотичностью, а корреляционная функция сигнала с ограниченной выборкой п должна совпадать
с корреляционной функцией выбранного временного ряда. При этом энтропия сигнала оценивается как
0,5
н[р«т= J inр»Ш(о, (и)
-0,5
где спектральная плотность мощности временного ряда:
со
Pxx(w) = TRxx(n)e-JG,n • (12)
-СО
Оценка спектра мощности с максимальной энтропией для ряда конечной длины —р <п< р приводит к следующему выражению:
РГ(а) = -— -----------• (13)
£ к(п)е-]0т
п=-р
Если знаменатель в (13) представить в виде симметричного полинома р 1
Y. k(n)z-n =—A(z)A(z-1), (14)
п=-р СТ
где о-2 — нормирующий множитель, а A(z) — полином порядка р, заданный как
A(z) = l + a1z~1 + ... + apz~p, (15)
то из (13) и (14) спектр мощности может быть выражен как
2
Сах(®) = —-—Г- (16)
A(z)A(z-1)
Соотношение (16) показывает, что оценка спектра мощности с максимальной энтропией — это спектр мощности АР-модели.
Моделирование и исследование параметрических методов
В нашей работе предварительно была исследована эффективность методов (ковариационного, модифицированного ковариационного и метода Берга) при расчетаx спектральной плотности мощности и определении частоты. Для этого были разработаны программы в среде MATLAB 7.0.1 (release 14). При одинаковьік вxодныx условияx проводилась оценка спектральной плотности мощности зашумленного сигнала (а на ее основе — частоты сигнала) 100 раз подряд. Если разность известной частоты чистого сигнала и оцененной частоты зашумленного сигнала составила меньше 3 % от значения частоты чистого сигнала, то такая оценка считалась достоверной. На основе этик данныч оценивалась эффективность. Анализ показал, что при анализе длинныx последовательностей отсчетов все методы дают практически одинаковые результаты, различия начинают проявляться в случае коротка сигналов. При использовании авторегрессионныx методов важно правильно выбрать порядок авторегрессионной модели — он должен быть в два раза большим числа синусоидальныx колебаний, которые предположительно содержатся в анализируемом сигнале. Метод Берга и модифицированный ковариационный метод позволили получить высокую разрешающую способность при анализе короткиx сигналов. Однако положения спектральныx пиков оказываются зависящими от начальныx фаз синусоид и при анализе суммы синусоид с шумом спектральные пики получаются слегка смещенными. Наилучшие результаты спектрального разрешения при ограниченной выборке сигнала с низким отношением сигнал/шум были достигнуты при использовании метода максимальной энтропии. На рисунке 1 показаны графики зависимости определения частоты сигнала ССИ от отношения «сигнал/шум» (SNR), где SNR варьировался от 0,01 до 1 с шагом 0,01. Видно, что при 0,1 < SNR < 1 пик сигнала xорошо различим, а при SNR < 0,1 он почти сливается с шумом и точное значение частоты определить практически невозможно.
Рисунок 2 отображает зависимость эффективности метода максимальной энтропии от SNR (от 0,01 до 0,3) и числа отсчетов сигнала N (от 1000 до 100 с шагом 100 плюс еще одно значение — 50 точек). Сильные спады заметны при значенияx SNR < 0,05 и N < 200.
Frequency, Hz
Рис. 1. Зависимость оцениваемой частоты в Гц от отношения «сигнал/шум» (SNR) в относительных единицах для 1000 отсчетов сигнала ССИ
Efficiency of the Maximum of Entropy method
Рис. 2. Зависимость эффективности метода в процентах от отношения «сигнал/шум» в относительных единицах и числа отсчетов N
Эксперимент
Экспериментальная установка собрана на основе консоли Аполло производства компании Tecmag (США). В качестве усилителя мощности использован линейный усилитель мощности 500 Вт производства компании TOMCO (Австралия). В состав датчика спектрометра вxодят высокодобротный резонансный контур, низкошумящий предусилитель (коэффициент усиления не менее 30 дб, э. д. с. шума не более 2,7 нВ) и система гашения «звона» после действия мощного радиочастотного импульса в резонансном контуре. Система гашения «звона» позволяет за счет изменения добротности резонансного контура сократить время нечувствительности приемника до 500 мкс. В последовательном резонансном контуре использовался переменный конденсатор 5 — 400 пф, 7,5 КВ, 20 А, производства фирмы Jennings, а также радиочастотная катушка с образцом. Добротность нагруженного последовательного контура, измеренного по ширине линии стандартным способом, составляла при измеренияx от 100 до 150. Образец представляет чешуированный тринитротолуол (TNT). Вес используемого образца составляет 100 г.
Регистрация осуществлялась на резонансной частоте 842 кГц с применением ограниченной выборки данных, 200 точек с интервалом 5 мкс, минимальное частотное разрешение при этом может составлять до 1 кГц, количество накоплений при этом равно 4000.
При комнатной температуре в спектре ЯКР образца имеется дублет 842 и 844 кГц, который разрешим при увеличении интервала окна наблюдения Та [10]. Фурье-спектр сигнала показан на рисунке 3. Как мы видим, разрешить тонкую структуру линий ТМТ не удалось. Это объясняется недостаточным количеством точек при данном частотном спектре. Так как полный частотный спектр лежит в диапазоне до 200 кГц, а частотное разрешение при преобразовании Фурье примерно около 1 кГц. Однако реализация в цифровом приемнике возможно только при числе отсчетов 2П, поэтому в данном случае количество отсчетов в частотной области составляет 256. Остальные точки заполняются нулями. При этом несколько увеличивается временной интервал окна. Хотя заполнение нулями и является стандартной процедурой в фурье-спектроскопии, оно не улучшает разрешающую способность преобразования по данной последовательности данных. Дополнение нулями позволяет получить интерполированное преобразование более сглаженной формы, а также точность оценивания частоты спектральных пиков [5]. В данном случае в области нахождения дублета мы имеем три точки, что не позволяет нам сделать эффективное разрешение. Использование авторегрессионного метода приводит к расширению эффективного окна наблюдения. Например, 256 обработанных точек с помощью авторегрессионного метода при интервале между точками 5 мкс дают эффективное окно наблюдения Тэ = 1280 мкс, что приводит к улучшению разрешения до 0,78 кГц.
Рис. 3. Необработанный спектр, полученный из временного окна наблюдения 1 мс, 200 точек с интервалом 5 мкс.
Для преобразования Фурье использовано 256 точек с заполнением нулям и дополнительных точек (показана часть спектра)
В результате применения метода максимальной энтропии при обработке сигналов ЯКР, полученных от образца ТМТ на резонансной частоте 842 кГц при N = 200 отсчетов (ограниченная выборка) на рисунке 4, представлен спектр обработанного сигнала.
Рис. 4. Спектр, полученный преобразованием Фурье из сигнала, обработанного авторегрессионным методом, использовано 256 точек (показана часть спектра)
Из рисунков 3 и 4 видно, что метод ММЭ позволяет получить хорошее спектральное разрешение и выделить спектральные составляющие сигнала TNT на частотах 844 и 848 кГц при ограниченной выборке, что подтверждает эффективность метода и результаты, полученные с помощью моделирования.
Заключение
В данной работе авторегрессионные методы были применены для обработки сигналов ЯКР в низкочастотной области. В качестве объекта выбран дублет в спектре TNT 842, 844 кГц с окном наблюдения 1 мс. Предварительно ковариационный метод, модифицированный ковариационный и метод Берга (ММЭ) были проанализированы при обработке модельных сигналов. Наилучшие результаты были получены с помощью метода максимальной энтропии Берга, который и был применен для обработки реальных сигналов. Хотя ММЭ уже применялся в ЯКР, однако его использование было строго ограничено случаем отношения сигнал/шум > 1 [11]. В нашей работе анализ проведен для случаев отношения сигнал/шум < 1. Также следует отметить, что помимо некоторого увеличения отношения «сигнал/шум», расширение эффективного временного окна наблюдения с помощью авторегрессионного метода дает возможность улучшить разрешение спектра, что не позволяет сделать простое заполнение нулями дополнительных точек в импульсной фурье-спектроскопии. При этом следует не забывать, что минимальная различаемая частота в спектре определяется соотношением Av<l/Ta.
Список литературы
1. Miller J. B., Barrall G. A. Explosives Detection with Nuclear Quadrupole Resonance. American Scientist, vol.93, 50-57, 2005.
2. G. V. Mozzhukhin, S. V. Molchanov, G. S. Kupriyanova, A. V. Bodnya, V. V. Fedotov, Hao Guoxin, Jin Yanbo, Ren Tianliang, Zhang Guojin. The Detection of Industrial Explosives by the Quadrupole Resonance Method: Some Aspects of the Detection of Ammonium Nitrate and Trinitrotoluene; In:Explosives Detection using Magnetic and Nuclear Resonance Techniques. Series: NATO Science for Peace and Security Series. Subseries: NATO Science for Peace and Security Series B: Physics and Biophysics. Fraissard, Jacques; Lapina, Olga (Eds.) Springer 2009, 295 p., 231 —244.
3. Сафин И. А., Осокин Д. Я. Ядерный квадрупольный резонанс в соединениях азота. М., 1977. С. 20.
4. NTNMR Reference Manual. Tecmag inc. URL: http://www.tecmag.com
5. Gregorovic A., Apih T. TNT detection with 14N NQR: Multipulse sequences and matched filter // Journal of Magnetic Resonance, 198 (2009), 215 — 221
6. P. A. Angelidis. Spectrum Estimation and the Fourier Transform in Imaging and Spectroscopy. Concepts in Magnetic Resonance, Vol. 8(5) 339 — 381 (1996).
7. Марпл.-мл. С. Л. Цифровой спектральный анализ и его приложения / пер. с англ. М., 1990.
8. Сергиенко А. Б. Цифровая обработка сигналов. СПб., 2002.
9. Burg I. P. Maximum Entropy Spectral Analysis. Proc. 37th Meeting of the Society of Exploration Geophysicists. Oklahoma City, Okla., October 1967.
10. Mozzhukhin G. V., Bodnya A. V., Fedotov V. V. et al. The Detection Of The Trinitrotoluene By Pure Nuclear Quadrupole Resonance. Appl. Magnetic Resonance, 29, 2005, 293 — 298.
11. СинявскийН.Я. Известия вузов. Физика. 1990. № 12. С. 91—92.
Об авторах
С. В. Молчанов — ст. преп., РГУ им. И. Канта.
Г. В. Мозжухин — канд. физ.-мат. наук, доц., РГУ им. И. Канта.
Authors
S. Molchanov — IKSUR.
G. Mozzhukhin — Dr., IKSUR.
