CC BY
23
УДК 551.510
ИССЛЕДОВАНИЕ АЛГОРИТМА ЧАСТОТНО-ВРЕМЕННОГО АНАЛИЗА РЯДОВ ИЗМЕРЕНИЙ СПЕКТРОВ ИЗЛУЧЕНИЯ НОЧНОГО МЕЗОСФЕРНОГО ОЗОНА В МИЛЛИМЕТРОВОМ
ДИАПАЗОНЕ ВОЛН С. Б. Розанов, А. С. Завгородний, А. Н. Игнатьев
Для изучения вариаций излучения ночного мезосфер-ного озона, наблюдаемых на миллиметровых волнах с поверхности Земли, разработан и испытан алгоритм частотно-временного анализа нерегулярных по времени рядов данных. Приведены некоторые результ,ат,ы численных и натурных экспериментов.
Ключевые слова: озон, миллиметровые волны, мезосфера, спектрометр, выборка, частотный анализ, окно наблюдений, численное моделирование.
Введение. Атмосферный озон защищает жизнь на Земле от ультрафиолетового излучения Солнца, влияет на тепловой баланс и динамику атмосферы, на климат [1]. В ночное время содержание озона выше 50 км, в мезосфере, значительно увеличивается. Наземные измерения на миллиметровых (ММ) волнах позволяют получать данные о вертикальном распределении озона (ВРО) на высотах до 90-100 км [2, 3]. В ходе работ, проведенных в 1999-2006 гг. в ФИАН и ИФА РАН, были обнаружены быстрые вариации излучения ночного озона с периодами от нескольких минут до нескольких часов [2, 3]. Создание в ФИАН нового автоматизированного озонометра ММ диапазона с высокой чувствительностью [4] дало возможность более детально исследовать эти вариации. Целями данной работы были выбор и исследование характеристик алгоритма частотно-временного анализа неравномерных по времени рядов измерений спектров ночного озона на ММ волнах, позволяющих выявлять колебания интенсивности (яр-костной температуры) излучения озона, определять их периоды и длительности.
Аппаратура и методика наблюдений. Единичные спектры излучения озона (линия с центральной частотой f0 = 142.175 ГГц, длина волны 2.1 мм) накапливаются в течение
ФИАН, 119991 Россия, Москва, Ленинский пр-т, 53; e-mail: [email protected].
100 с и записываются каждые 110 с с помощью 50-канального фильтрового анализатора спектра АС-50 с разрешением 0.1 МГц в центре линии [4]. Накопление спектров периодически, обычно раз в 25 минут, прерывается для калибровки озонометра и уточнения его параметров. В результате данные наблюдений представляют собой неравномерные по времени ряды значений яркостных температур в каналах АС-50.
Алгоритм частотно-временного анализа. В [3] показано, что ВРО в мезосфере может быть определено по разностям яркостных температур линии озона, измеренных в каналах АС-50 с отстройками частоты 0, 50, 150 и 250 кГц от центра линии. Ряды разностей яркостных температур единичных спектров и являются объектом частотно-временного анализа. Наиболее известные методы такого анализа: оконное преобразование Фурье (ПФ) [5] и вейвлет-преобразование [6], применимы только к равномерным выборкам. Интерполяция неравномерных данных на равномерную сетку может, как отмечалось в [7], приводить к потере в спектре мощности сигнала высоких частот, приближающихся к половине частоты дискретизации. Задача адаптации аппарата вейвлет-преобразования для анализа неравномерных выборок, а главное, последующего определения статистических характеристик получаемых результатов весьма трудоемка, её решение в рамках данной работы было признано нецелесообразным.
Метод частотного анализа неравномерных рядов данных, свободный от недостатков ПФ, был предложен Ломбом [8] и Скарглом [7] с использованием результатов Барнинга [9]. Пусть выборка состоит из N измерений X? величины X (£), полученных в моменты времени £ = , ] = 1,..., N, ¿3-+1 > . Периодограмма Ломба (спектр мощности) величины X(£) определяется как [8]:
Яьх (ш)
1 - <
2
N
^ (X? — XX) сое ш(£? — т) 3=1
N
+
N
(X? — X) вт ш(£? — т)
3 = 1
N
^ сое2 — т)
3=1
^ вт2 ш(£? — т)
3=1
(1)
где ш - круговая частота, X - среднее значение X(£), временной сдвиг т определяется как
N
^ вт (2ш£?)
^(2шт) = 33=-
^ сов(2ш£3-)
3=1
(2)
В [7, 8] показано, что периодограмма Ломба (1) инвариантна по отношению к сдвигу шкалы времени, а её статистические свойства для неравномерных выборок оказывают-
2
2
ся такими же, как свойства преобразования Фурье для равномерных выборок. Так, если X(¿) - случайная гауссова величина с дисперсией аХ, то плотность вероятности того, что (ьх(ш) = гаХ, равна р(г) = е-* (порог г > 0), а среднее по частоте значение QLX (ш) равно аХ [7]. Т.к. на результаты измерений накладывается шум аппаратуры, требуется оценить вероятность появления в результатах расчёта (ьх (ш) шумового выброса, превышающего уровень гаХ на одной из М контролируемых частот. В [7, 10] показано, что эта вероятность "ложной тревоги" равна
Ро = 1 - (1 - е-*)м.
Вероятность р0 на практике выбирается достаточно малой, р0 < 0.1. В этом случае
р0 и Ме-г. (3)
Для выборок с небольшим нарушением регулярности, как в нашем случае, принимают М и N [10]. Тогда из (3) получаем
г и 1п(^ро). (4)
Частотный анализ методом Ломба позволяет решить задачу обнаружения колебательных процессов в рядах данных, однако для изучения эволюции наблюдаемых явлений нужен частотно-временной анализ, при котором часть выборки выделяется скользящим окном наблюдений. При использовании окна наблюдений вместо исследуемой величины X(¿) берётся её произведение на некоторую весовую функцию (в нашей работе - гауссовой формы) с полушириной Тш, а в формуле (4) общее число точек выборки N заменяется на их эффективное количество в пределах окна.
Численные эксперименты. Алгоритм, созданный для анализа результатов измерений спектров излучения ночного озона, базируется на соотношениях (1) и (2) с учётом рекомендаций работ [8, 10]. В качестве величины X(¿) используются измеряемые в градусах Кельвина (К) разности яркостных температур линии озона в центральных каналах фильтрового анализатора спектра [4]. Для проверки разрешения алгоритма по частоте и времени использовался тестовый сигнал в виде суммы трёх гармоник с амплитудами 10, 5 и 7 К на частотах 1, 2 и 3 мГц соответственно. Каждая из гармоник "включалась" на своем интервале времени, интервалы частично перекрывались. Отсчеты времени выборки с количеством элементов N = 1000 формировались как сумма регулярной (с периодом 110 с) и случайной (с СКО а^ = 5 с) составляющих. В результате обработки такого сигнала с использованием окон с полушириной Т,ш = 5 — 30 минут
были определены моменты "включения" и "выключения" гармоник, оказавшиеся в хорошем согласии с характеристиками исходного сигнала. Средние уровни спектральных составляющих гармоник соотносятся примерно как 4:1:2 ^ 102 : 52 : 72, т.е. соответствуют квадратам амплитуд (мощностям) гармоник.
Чтобы смоделировать выделение гармонического сигнала из шумов аппаратуры, использовались тестовые сигналы
Х1(Ь) = в 8т(2п^) + СдЦ)
(5)
в виде суммы гармоники с частотой = 1 мГц и амплитудой и центрированного гауссова шума д(£) с единичной дисперсией, масштабируемого множителем О =1.6 К, соответствующим шуму аппаратуры при времени накопления спектра 100 с. Число элементов выборки N = 300 соответствует типичному числу единичных спектров озона в течение зимней ночи. На рис. 1(а) показаны периодограммы Ломба (1) чисто шумового сигнала, когда = 0, для одной реализации шума (тонкая линия) и среднее по 10 реализациям (толстая линия). Горизонтальные линии на рис. 1(а) отмечают вероятности "ложной тревоги" р0 = 0.5, 0.1, 0.05, 0.01. Видно, что для единичной реализации шума максимальные пики лежат между значениями 0.5 и 0.1 вероятности р0 . Поэтому
Рис. 1: (а) периодограммы Ломба (1) для центрированного гауссова шума: для одной реализации шума (тонкая линия) и средний спектр по 10 реализациям (толстая линия); (б) периодограмма Ломба для гармонического сигнала с частотой = 1 мГц и амплитудой Бг = 0.5 К в присутствии гауссова шума с дисперсией (О2 = 1.6К2).
регистрация гармонической компоненты сигнала X1 (t) в нашем случае может считаться достаточно надёжной при превышении ею уровня, соответствующего вероятности p0 = 0.05. Среднее по частоте значение мощности для единичной реализации шума составило 2.572, для среднего спектра по 10 реализациям - 2.492, обе эти величины близки к теоретическому значению G2 = 2.562. На рис. 1(б) приведена периодограмма Ломба (1) для сигнала Xi(t) при амплитуде гармоники S1 = 0.5 К. Горизонтальные линии соответствуют значениям p0 = 0.05, 0.01. Видно, что спектральная составляющая на частоте 1 мГц превышает шумы и достигает порога обнаружения, соответствующего вероятности p0 = 0.05 (обведено кружком). При этом отношение амплитуды сигнала к СКО шума составляет S1/G ~ 0.3.
В настоящее время алгоритм частотно-временного анализа используется для обработки экспериментальных массивов спектров излучения ночного мезосферного озона. Предварительный анализ выявил колебания разностей яркостных температур в центре линии с периодами от 5 до 40 мин и "временами жизни" процессов от 1 до 8 периодов.
зЗаключение. Обоснован выбор алгоритма частотно-временного анализа неравномерных по времени рядов измерений спектров излучения ночного мезосферного озона в линии 142.175 ГГц, основанного на методе периодограмм Ломба. Для выявления колебаний с "временами жизни" меньше длительности сеансов наблюдений использованы скользящие временные окна гауссовой формы. Численные эксперименты показали корректность работы алгоритма. При обработке экспериментальных данных выявлены вариации излучения озона с периодами от 5 до 40 мин, которые могут быть связаны с распространением внутренних гравитационных волн (волн плавучести).
Работа поддержана программами ОФН РАН "Радиоэлектронные методы в исследованиях природной среды и человека" и "Новые источники миллиметрового и терагер-цового излучения и их перспективные приложения". Авторы выражают благодарность д.ф.-м.н. С. В. Соломонову и к.ф.-м.н. Е. П. Кропоткиной за полезные обсуждения работы.
ЛИТЕРАТУРА
[1] Scientific Assessment of Ozone Depletion: 2014. Report No. 55 (WMO, Geneva, 2014), http://www.esrl.noaa.gov/csd/assessments/ozone/2014/report.html.
[2] С. В. Соломонов, Е. П. Кропоткина, А. И. Семенов, Краткие сообщения по физике ФИАН, № 10, 31 (2001).
[3] А. Н. Игнатьев, Радиометрия атмосферного озона и окиси хлора на миллиметровых волнах. Дисс. канд. физ.-мат. наук (МГУ, Москва, 2006).
[4] С. Б. Розанов, А. С. Завгородний, С. В. Логвиненко и др., Известия вузов. Радиофизика 54, (8-9), 708 (2011).
[5] А. Б. Сергиенко, Цифровая обработка сигналов. 2-е изд. (Питер, Санкт-Петербург, 2006).
[6] В. П. Дьяконов, Вейвлеты. От теории к практике (Солон-Р, Москва, 2002).
[7] J. D. Scargle, Astrophys. Space Sci. J. 263, 835 (1982).
[8] N. R. Lomb, Astrophys. Space Sci. 39, 447 (1976).
[9] F. J. M. Barning, Bulletin of the Astron. Institutes of the Netherlands 17, 22 (1963). [10] Spectral Analysis of Unevenly Sampled Data, in Numerical Recipies, 3rd Edition, Ed.
by W.H. Press et al., (Cambridge, 1986), Chap. 13, p. 685.
Поступила в редакцию 7 июня 2016 г.
