CC BY
52
РАДИОФИЗИКА
УДК 520
ПАРАМЕТРЫ АСТРОФИЗИЧЕСКИХ ОБЪЕКТОВ ПО ДАННЫМ О МОДУЛЯЦИИ ИНТЕНСИВНОСТИ ИХ ЭЛЕКТРОМАГНИТНОГО ИЗЛУЧЕНИЯ.
I. АЛГОРИТМЫ ОБРАБОТКИ ДАННЫХ НАБЛЮДЕНИЙ
© 2011 г. А.Г. Кисляков, Е.И. Шкелёв, С.Ю. Лупов, К.Г. Кислякова
Нижегородский госуниверситет им. Н.И. Лобачевского [email protected]
Поступила в редакцию 29.09.2010
Кратко описаны методы определения параметров астрофизических объектов на основе анализа спектрально-временных (СВА) особенностей модуляции интенсивности их излучения. Применён билинейный метод Вигнера - Виля (впервые в астрономии) в сочетании с методом Фурье-анализа со «скользящим окном», что позволило реализовать широкий диапазон время-частотного цифрового анализа (от ~10-4 до ~10 Гц). Приводятся основные алгоритмы методов СВА, а также процессов подготовки цифровых сигналов к анализу (например, частотной фильтрации, поиска подобных сегментов в сигнале и др.).
Ключевые слова: анализ Фурье, преобразование Вигнера - Виля, динамические спектры.
Введение
Электромагнитное излучение астрофизических объектов является основным источником информации о физических условиях и эволюции этих обычно недоступных для контактных измерений небесных тел. С течением времени растёт арсенал средств для наблюдений космического электромагнитного излучения (КЭИ), при этом увеличивается и количество исследуемых объектов. Изменяются интенсивность, спектр и временные характеристики наблюдаемого излучения, что требует адаптации методов регистрации и обработки сигналов к их особенностям. Следует учитывать и вариации сопутствующих естественному излучению различных мешающих сигналов, в основном инструментального происхождения. В последнем случае становится необходимой предварительная подготовка полученных сигналов к спектральновременному анализу (СВА), а иногда требуется и постановка контрольных измерений с целью отделения артефактов от истинных характеристик излучения.
В настоящей работе представлены методы определения физических характеристик источ-
ников КЭИ на основе СВА модуляции интенсивности этого излучения. Авторы применили, в частности, билинейный метод СВА Вигнера -Виля (ВВ; впервые в астрономии) в сочетании с методом Фурье-анализа со «скользящим окном» (СФА), что позволило реализовать широкий диапазон время-частотного анализа (исследованы квазипериодические модуляции с частотами от ~10-4 до ~10 Гц). СФА-ВВ широко применен в исследованиях солнечного излучения по данным наблюдений на таких радиотелескопах, как 14-м и 2.6-м радиотелескопы в Метсахови (Финляндия, ХТУ), 45-м радиотелескоп и радиогелиограф в Нобеяма (Япония) и др. (диапазон длин волн от X ~ 3 мм до X ~ 10 см). Обнаружена звёздная активность, аналогичная солнечной, при наблюдениях на 100-м радиотелескопе в Эффельсберге (Германия, X ~ 5 см). Исследования звёздной активности в дальнейшем развивались на основе данных наблюдений в УФ-диапазоне с борта различных ИСЗ (MOST, COROT и др.). Внеатмосферные наблюдения
КЭИ на волнах X ^ 102 - 104 м также стали основным источником данных о столь длинноволновом излучении вследствие его экранирования ионосферой Земли. При этом в число
объектов исследований вошли планеты-гиганты Солнечной системы, их крупнейшие спутники, а также плазма межпланетной среды (ИСЗ CASSINI, ULYSSES, STEREO и др.). Ниже приводятся краткое описание алгоритма СФА-ВВ и некоторые полученные результаты.
1. Этап подготовки данных наблюдений к обработке
Поскольку почти во всех практически значимых случаях астрофизические данные представляют собой дискретные ряды отсчётов интенсивности сигналов (за исключением фотографий в оптике), естественно применение цифровых методов их обработки. Одна из неизбежных особенностей цифровых методов заключается в ограниченности динамического диапазона получаемых данных. В результате из поля зрения могут исчезнуть относительно слабые, но представляющие интерес составляющие спектрально-временной картины события. При подготовке данных их подвергают следующим операциям: 1) вычитанию постоянной составляющей (или среднего); 2) вычитанию основного «тренда» изменения данных, если он является достаточно медленным по сравнению с представляющими интерес осцилляциями излучения; 3) аппроксимации ряда данных достаточно медленным полиномом (методом наименьших средних квадратов) с последующим вычитанием полученной кривой; 4) «нормировке» (см. ниже); 5) цифровой фильтрации доминирующих частотных составляющих, не представляющих интереса. Заметим, что цифровая фильтрация с помощью подходящих фильтров может избавить от необходимости перечисленных вычитаний. В качестве иллюстрации рассмотрим пример обработки профиля двойного всплеска солнечного радиоизлучения, наблюдавшегося с помощью 14-м радиотелескопа обсерватории Мет-сахови на частоте 37.5 ГГц [1].
Отсчёты интенсивности солнечного радиоизлучения производились периодически с интервалом At = 0.1 с. Всплески излучения инициированы солнечными вспышками в активной области [1], на которую ориентирована антенна радиотелескопа. Заметно сходство профилей наблюдавшихся всплесков: вслед за мощным кратковременным «предвестником» (рис. 1а) происходят постепенный рост и спад интенсивности радиоизлучения (всплески типа ПРС). В результате вычитания аппроксимирующей кривой (рис. 1б) выделяются мощные квазиперио-дические осцилляции в области «предвестни-
ков», в то же время можно отметить затухание осцилляций с уменьшением интенсивности всплеска и наличие интервалов быстрых флуктуаций излучения. Всё это говорит о сложном механизме возникновения и развития всплеска, рассмотрение которого выходит за рамки настоящей работы. Аналогичный результат (см. рис. 2) получается при фильтрации данных с помощью цифрового фильтра верхних частот (см. ниже). Заметно сглаживание флуктуаций вследствие инерционности фильтра. Существенно, что в результате частотной фильтрации теряется также информация об интенсивности флуктуаций радиоизлучения.
Перейдём теперь к рассмотрению метода «нормировки». Этот приём основан на использовании аналитического сигнала z(n) [2]: z(n) = s(n) + jsT (n) = M(n)exp(jy(n)), (1.1)
где s(n) - исследуемая цифровая реализация сигнала, sI.(n) - гильбертова сопряжённая функции s(n), а n - номер отсчёта. Огибающая (или модуль) M(n) и фаза y(n) аналитического сигнала определяются выражениями
M (n) = \s2(n) + s,?(n)J
/2
и y(n) = arctg
■sy (n) s(n)
«Нормировка» заключается в делении г(п) на М(п), при этом выравниваются амплитуды всех спектральных составляющих процесса ^(п). Заметим, что функция М(п) не обращается в ноль вследствие ортогональности функций s(n) и sг(n). При компьютерном анализе цифровых данных затруднено использование интегральной формы sг(n), т.к. эта функция является свёрткой функций s(n) и 1/п, поэтому знаменатель подынтегральной функции может обращаться в ноль. В связи с этим функция sг(t) обычно находится путём двойного преобразования Фурье. Сначала ищется спектр исследуемой реализации s(a)^s(n), а затем осуществляется переход sг(a)^sг(n), т.к. между спектрами дискретной последовательности и её гильбертовой сопряжённой существует простая связь: 5г(ю) = js(a) [2]. Заметим, что приём нормировки даёт хорошие результаты при анализе простых нестационарных процессов, при этом удаётся проследить вариации «мгновенной» частоты колебания.
Для иллюстрации сказанного рассмотрим пример спектров тестового сигнала:
w(t) = ехр[-0.01(? - 15)2] 8т(6.28? + 0.1?2 ). (1.2) График функции (1.2) представлен на рис. 3а.
Дискретная реализация этой функции при Дt = 0.01 с была подвергнута спектральному анализу с использованием алгоритмов Фурье
Рис. 1. а) Профиль солнечных всплесков 28.08.1990 на частоте 37 ГГц; б) флуктуации солнечного излучения 28.08.1990 (всплески вычтены)
Рис. 2. Флуктуации солнечного излучения 28.08.1990 после частотной фильтрации
(ДПФ) и Вигнера - Виля (ДВВ) (см. графики рис. 3б и 3в, соответственно; алгоритмы описаны в разделе 2). Очевидно, что преобразование ДВВ даёт более высокое разрешение на плоскости (га,0, чем ДПФ. На рис. 3г, д показано также влияние нормировки сигнала на его спектр. В обоих случаях спектры принимают вид простого ЛЧМ-сигнала, при этом преобразование ДВВ обеспечивает более высокое разрешение на плоскости (га, ?).
В работе [3] рассматриваются примеры более сложных тестовых сигналов, являющихся суперпозицией импульсных и гармонических сигналов, при этом неизбежно их взаимодействие (в случае преобразования ВВ) и появление ложных сигналов. Подобные тесты позволяют на практике отделять искомые сигналы от
артефактов, при этом полезно для контроля использовать метод ДПФ.
Необходимым тестом данных измерений является проверка их эквидистантности, что возможно при условии, что ряд отсчётов сопровождается сведениями о времени измерения. Такая информация может быть использована также для фильтрации мешающих импульсов. Данные наземных инструментов обычно эквидистантны и период следования отсчётов сохраняется достаточно хорошо. Имеется также возможность постановки контрольных измерений для выявления помех.
При наблюдениях с борта ИСЗ ситуация усложняется и необходим более тщательный контроль данных наблюдений. Рассмотрим пример обработки данных наблюдений УФ-из-
Рис. З. Иллюстрация действия «нормировки» на спектр модельного сигнала после его Фурье и ВВ преобразований
лучения звезды G7V и планеты Exo2b с борта ИCЗ COROT. Из рис. 4а видно, что в данных присутствуют импульсы с периодом следования —512 с, более мощные, чем отсчёты световой кривой (CK). Разработана программа «чистки» данных от этих импульсов инструментального происхождения. Рисунок 4б показывает поток данных после «чистки» с постоянным периодом отсчётов З2 с. Исправленная CK использовалась для определения модуляций её интенсивности вследствие вращения звезды, а также при транзитах планеты Exo2b перед диском звезды (см. [4-6], а также раздел 3).
Для полноты этого раздела приведём алгоритмы использовавшихся нами цифровых фильтров. Aлгоритм работы фильтров основан
на наложении гауссова окна на спектр сигнала, полученный с помощью дискретного преобразования Фурье (ДПФ), с последующим вычислением обратного ДПФ. Г ауссово окно определяется следующими выражениями:
■d • N
W-нч (i )= expj i 2 •^
i = 0,..., N/2, - для ФНЧ,
WB4 (i) = 1 - exp|- i2 2-1 • f • d • N i = 0,..., N/2, - для ФВЧ,
Wim (i) =
i = 0,..., N/2, - для ППФ, (i ) =1 - exp|-i = 0,..., N/2, - для ПЗФ,
где Л - частота среза для фильтров нижних и верхних частот, Л - частота дискретизации, /0 -центральная частота полосы пропускания (для полосно-пропускающего фильтра - ППФ) и запирания (для полосно-задерживающего фильтра - ПЗФ), ДЛ- полоса пропускания (задержания) для полосовых фильтров, ё - коэффициент децимации (см. также [2, 7]).
Приведённые выше фильтры не исчерпывают инструментарий, используемый при подготовке сигнала к анализу. Разработан пакет простых и эффективных программ (работают на Фортране или Си в зависимости от формата данных) для осреднения данных на скользящем интервале, для вычитания среднего или тренда, а также для других простых преобразований. Программы позволяют быстро обрабатывать файлы ёмкостью до 150 МБ и более (порядка 107 отсчётов).
2. Алгоритмы спектрально-временного анализа
В качестве примера на рис. 5а показан экран компьютера, иллюстрирующий процесс подготовки сигнала к СВА. Вслед за командой «Открыть» на экране высвечивается реализация сигнала или её часть, соответствующая положению курсора. При этом активируются кнопка «Набор» (выбор цифрового фильтра, децимации
2
Рис. 4. Чистка сигнала (световой кривой) от инструментальных помех
їда* к*. ,т, чюч • і
ШАГ ) 2ІІ Г. хн М Мка[Го :• СОХ.КЛК
Рис. 5. а) Окно СВА при подготовке сигнала к анализу; б) панель управления процессом СВА
и нормировки сигнала) и кнопки выбора анали- зации может быть сохранена в виде Ш-файла. зируемого интервала сигнала. После фиксации То же самое может быть сделано после цифро-этого интервала соответствующая часть реали- вой фильтрации или иной трансформации сиг-
нала. Трансформированный сигнал также показывается на дисплее компьютера. Нажатием кнопки «Анализ» переходим к процессу СВА (рис. 5б). Активированные кнопки указывают, что (после нажатия кнопки «Старт») произведён анализ спектра методом Фурье при окне Хеннинга размером 2048 точек. Прочие кнопки определяют выбор уровней отсчёта спектральной плотности (по модулю) и интервала анализируемых частот.
Для дискретных во времени сигналов si, состоящих из N отсчетов, для вычисления спектра Sk, также состоящего из N отсчетов, используется ДПФ [2, 8]:
N-1
Sk = Z wnd(i)Si
. 2nik
1 N
k = 0,1,...,N -1.
i=0
Для получения спектрально-временного распределения ДПФ вычисляется в скользящем окне теп^/). Используемые функции окон приводятся ниже (см. также [9]):
( ч 11, 0 < і < N
(і) = •< - прямоугольное окно;
wnd2 (i ) =
cos21 — I, 0 < i < N/2
2N) - окно Хеннинга;
! > N 2
1 -6•'тт] + 6•( 1^] , 0 <i <N/4
wnd3 (i) =
2 •! 1 -
2-і
N
з
N/4 < i < N/ 2 i > N/ 2
с 11-летним солнечным циклом. Результирующий динамический спектр представлен в левой части дисплея на рис. 5б. Сечение этого динамического спектра по вертикали (отмечено курсором на спектре) показано в правой части дисплея. Динамический спектр может быть сохранен в рабочей папке анализа в виде Ьшр-файла. Нажатие кнопок СТКЬ-ІКБ активирует сохранение в той же папке сечения динамического спектра в виде Ш-файла. Подробное рассмотрение вариаций ИКМА Солнца выходит за рамки настоящей работы.
Разработан специальный вариант СВА, позволяющий сохранять необходимое число сечений динамического спектра. Это число определяется децимацией динамического спектра, а также выбором интервалов анализируемого сигнала. Количество сохраняемых спектров достаточно велико (до 105-106 и более) и ограничивается только возможностями используемого графического редактора. В результате можно получить усреднённые по заданным интервалам времени динамические спектры, что даёт возможность выделения (в эти периоды времени) представляющих интерес слабых частотных составляющих.
При инициации процесса СВА методом Вигнера - Виля (кнопка Виг-Виль на рис. 5б) также используются подобные прежним установки: выбор интервала анализируемого сигнала, децимация, нормировка, уровни отсчёта спектральной плотности и т.п. Команда «Старт» даётся после всех необходимых установок. Алгоритм дискретного преобразования Вигнера -Виля (ДВВ) определяется выражением
2N-1
E (kAt, mAf) = 2 At Z
i=0
Sk+iSk-i єхРІ - І
. n im
N
- окно Парзена.
К достоинствам данного метода следует отнести его быстродействие, особенно если использовать для вычислений алгоритмы быстрого преобразования Фурье (БПФ) [2, 8]. К недостаткам можно отнести тот факт, что разрешение по частоте обратно пропорционально количеству отсчетов в окне, т.е. размер окна необходимо выбирать большего размера, что сказывается на разрешении по времени.
Сигнал, часть реализации которого можно видеть на рис. 5а, описывает суточные вариации индекса корональной магнитной активности (ИКМА) Солнца за период 1939-2005 гг. Можно видеть вариации этой активности, связанные
где £(кД?,тД/) - дискретное распределение Вигнера - Виля энергии сигнала на плоскости (?, /); S - последовательность комплексных отсчётов исследуемого процесса (аналитический сигнал), Д - период дискретизации, Д/ = = 1/(4МД?), к = 0, 1, 2,..., 2 N т = 0, 1, 2,..., 2К Исследование нестационарных процессов проводится обычно с использованием скользящего окна и БПФ. Выше уже приводился пример обработки нестационарного тестового сигнала, где показано преимущество преобразования ВВ в реализации более высокого разрешения на вре-мя-частотной плоскости (см. рис. 3а-д).
Заметим, что распределение ВВ является вещественным лишь при интегрировании непрерывных функций в бесконечных пределах. Дискретное преобразование в конечных пределах даёт в результате комплексную функцию
e
0
Е(Х/), которая в литературе получила название «псевдо-ВВ преобразования» [10]. В связи с этим приходится оперировать модулем функции Е(/ (или её реальной частью).
3. Программа «Подобные сегменты»
Программа SimilarSegments (Подобные сегменты, ППС) предназначена для обработки сигналов, содержащих повторяющуюся информацию. В отличие от гармонического анализа, такой подход позволяет выделять и исследовать периодические процессы с импульсами сложной или даже изменяющейся формы.
Для поиска в исследуемом сигнале х участков, наиболее похожих на заданный, необходимо найти все значения параметра п2, при которых величина
/=хк+,-(а+к -))2
В выражении (3.1): п1 - начало выделенного сегмента сигнала, п2 - начало предполагаемого похожего сегмента, M - количество отсчетов в выделенном сегменте,
М-1 ' X
М-1
X хщ+,- В х
Х„
А =
М
Если масштабирование искомого сегмента не требуется, то 5=1, в противном случае
М-1 М-1 М-1/ ч
ХП\ +, х^2 +,
В = ,=0___________,=о_______
М-1
МХ(
Хи1 +, ' Хп2 +, I
(М-1
х„
П2 +,
V ,=0
М-1
- м X х2
(3.1)
имеет локальные минимумы (чем меньше /, тем более найденный участок похож на заданный).
Масштабирование может потребоваться при анализе процессов с изменяющейся интенсивностью.
На рис. 6а показан дисплей компьютера при работе программы ПС. После открытия файла с сигналом (кнопка «Открыть») реализация сигнала может быть преобразована путём нажатия кнопок в верхнем левом углу экрана (действие
а)
б)
2
2
,=0
Рис. 6. а) Панель управления поиском подобных сегментов; б) средний период обращения планеты и средний транзитный импульс
кнопки обозначено стрелками). Масштабирование и выделение сегмента сигнала осуществляется с помощью мышки. На рис. 6а представлен результат выделения сегмента сигнала. Задано число искомых подобных сегментов (затем нажата кнопка «Поиск»).
В результате поиска найденные сегменты сохраняются (после команды «Сохранить») в рабочей папке программы под заданным именем в виде txt-файла.
Результаты анализа СК звезды G7V при вращении планеты Exo2b (данные ИСЗ COROT) приводятся на рис. 6б. На основании подобных данных можно судить как о сейсмической активности звезды, так и о характере её взаимодействия с планетой; всё это будет предметом другой работы.
Заключение
В заключение заметим, что обработка значительных по объёму данных происходит чаще всего в условиях, когда ни время, ни форма проявления ожидаемого феномена не определены априори. Надёжность результатов анализа, несомненно, увеличится при использовании итеративных методов анализа, допускающих модификацию параметров обработки и адаптацию методов путём многократного повторения ключевых процедур. Это могут быть процедуры как подготовительного, так и итогового этапа. В связи с этим приходится повторять одни и те же операции при различных начальных условиях и с разной степенью детализации исследуемого сигнала. Это неизбежно, например, при анализе нестационарных сигналов с изменяющимися во времени параметрами.
Детальное рассмотрение функционирования комплексной программы СВА, комбинирующей «скользящее» преобразование Фурье с преобразованием Вигнера - Виля, содержится в работе [11], где представлена модель системы с использованием терминологии и обозначений, принятых в унифицированном языке моделирования UML [12].
Исследование выполнено совместно с сотрудниками Института космических исследований АН Австрии М.Л. Ходаченко и Х. Ламме-
ром, обеспечившим доступ к данным измерений с ИСЗ.
Работа поддержана конкурсным контрактом КД НК-21П с Федеральным агентством образования и грантом № 228319 Европейского союза в рамках проекта EuroPlanet RI FP7.
Список литературы
1. Urpo S., Puhakka P., Oinaskallio M. et al. Selected Radio Maps and Major Solar Radio Flares Measured at Metsahovi in 1996-2001. 2003, HUT-MET-46.
2. Марпл-мл. С.Л. Цифровой спектральный анализ и его приложения / Пер. с англ. О.И. Хабарова и Г.А. Сидоровой. Под ред. И.С. Рыжака. М.: Мир, 1990. 584 с.
3. Шкелёв Е.И., Кисляков А.Г., Лупов С.Ю. // Изв. высш. уч. зав. Радиофизика. 2002. Т. 45(5). С. 433-442.
4. Khodachenko M.L., Kislyakov A.G., Weingrill J. et al. An Efficient Algorithm for Analysis of Stellar Light Curves and Detection of Transits // 1-st CoRoT International Symposium, Feb. 2-5, 2009, Paris, France.
5. Кисляков А.Г., Лупов С.Ю., Ходаченко М.Л. и др. Программы анализа вариаций света звёзд с целью поиска планет // Труды научной конференции по радиофизике. ННГУ, 2010.
6. Кисляков А.Г., Ходаченко М.Л., Ламмер Х. и др. Активность звезды G7V, связанная с планетой Exo2b // Труды научной конференции по радиофизике. ННГУ, 2010.
7. Блейхут Р. Быстрые алгоритмы цифровой обработки сигналов. М.: Мир, 1989.
8. Allen R.L., Mills D.W. Signal Analysis: Time, Frequency, Scale and Structure. Wiley-IEEE Press, 2004. 966 p.
9. Pollock D.S.G. A Handbook of Time-Series Analysis, Signal Processing and Dynamics. London: Academic Press, 1999. 782 p.
10. Коэн Л. // ТИИЭР. 1989. Т. 77. № 10. C. 72-120
11. Шкелёв Е.И., Лупов С.Ю. // Вестник Нижегородского университета им. Н.И. Лобачевского. Серия Радиофизика. 2004. № 1. С. 55-61.
12. Крэг Л. Применение UML шаблонов проектирования: Пер. с англ. 2-е изд. М.: Изд-во «Вильямс», 2002. 624 с.
PARAMETERS OF ASTROPHYSICAL OBJECTS ACCORDING TO THEIR ELECTROMAGNETIC EMISSION INTENSITY MODULATION DATA. I.
OBSERVATIONAL DATA PROCESSING ALGORITHMS
A.G. Kislyakov, E.I. Shkelev, S.Yu. Lupov, K.G. Kislyakova
A brief description has been given of parameter determination of astrophysical objects on the basis of spectraltime analysis (STA) of modulation peculiarities of their emission intensities. The application of the bilinear Wigner-Ville (WV) method (for the first time in astronomy) together with the sliding window Fourier-analysis (SWF) has made it possible to realize the broad-band time-frequency digital (from ~10-4 to ~10 Hz have been studied). Main algorithms of STA and digital signal processing before the analysis are presented (for instance, frequency filtering, search for similar signal segments, etc.).
Keywords: Fourier analysis, Ville -Wigner transform, dynamic spectra.
