Разработка автоматизированной системы для анализа геомагнитных вариаций на основе вейвлет-технологий Текст научной статьи по специальности «Компьютерные и информационные науки»

структуры и моделирование 2014. №4(32). С. 185-197

УДК 004.415.28:550.385.1

РАЗРАБОТКА АВТОМАТИЗИРОВАННОЙ СИСТЕМЫ ДЛЯ АНАЛИЗА ГЕОМАГНИТНЫХ ВАРИАЦИЙ НА ОСНОВЕ ВЕЙВЛЕТ-ТЕХНОЛОГИЙ

С.Н. Верзунов

ст. преподаватель каф. ИВТ, аспирант, e-mail: [email protected]

Н.М. Лыченко д.т.н., проф. кафедры ИВТ, e-mail: [email protected]

Кыргызско-Российский славянский университет, г. Бишкек

Аннотация. Разработаны методика и программные средства для анализа геомагнитных вариаций на базе вейвлет-преобразования. На примере анализа магнитных вариаций показана эффективность разработанной методики. В частности было выяснено, что для периода магнитных бурь характерно появление квазипериодических вариаций магнитного поля с периодами в 120 и 30 мин. Программные средства разработаны на языке Python 3, являются кроссплатформенными и имеют широкий функционал.

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

1. Предмет и задачи разработки

Изучению вариаций геомагнитного поля Земли посвящено множество работ, при этом в основном для анализа измерительной информации традиционно использовался аппарат Фурье-преобразования. Однако в силу явной нестационарности магнитовариационного сигнала Фурье-анализ не способен дать достаточно полной информации о сигналах. Учесть особенности нестационарных временных рядов возможно, используя вейвлет-технологии, применению которых в анализе геомагнитных данных в последнее время посвящено значительное количество работ [1, 2, 3].

Недостаток предложенных в настоящее время методик анализа временных рядов заключается в том, что они не учитывают особенностей магнитовариа-ционных данных — достаточно большой размерности (порядка 216 в типичном случае) и наличия пропусков, возникающих вследствие сбоев применяемого оборудования, а также помех в каналах передачи данных. В отличие от методов, предложенных в работах [2, 4], в настоящей работе предлагается использовать интерполяцию для заполнения небольших пропусков, всегда имеющихся в реальных данных такого рода, и Фурье-анализ (БПФ), поскольку он эффективен при исследовании данных большой размерности. В методе, предложенном в работе [1], для анализа магнитных вариаций используется дискретное

вейвлет-преобразование, которое позволяет анализировать сигнал только по масштабам и сдвигам, равным степени двойки. Однако для магнитовариацион-ного сигнала непрерывное вейвлет-преобразование более удобно; его некоторая избыточность, связанная с непрерывным изменением масштабного коэффициента и параметра сдвига, является при этом положительным качеством, так как позволяет более полно и чётко представить и проанализировать содержащуюся в данных информацию, поэтому в настоящей работе предложено использовать непрерывное вейвлет-преобразование.

Цель данной работы — разработать алгоритм и программные средства для анализа геомагнитных вариаций на основе вейвлет-преобразования, позволяющие сделать цифровую обработку временных рядов геофизических данных доступной широкому кругу учёных, занимающихся анализом вариаций магнитного поля Земли. Разработанная автоматизированная система предназначена, в первую очередь, для просмотра и вейвлет-анализа вариаций одной из составляющих геомагнитного поля, однако, она может быть использована для анализа любых одномерных данных, представленных в формате CSV (Comma Separated Value — значения, разделённые запятой). В качестве исходных данных для отладки и иллюстрации работы автоматизированной системы использованы данные интерактивного ресурса по солнечно-земной физике (Space Phisics Interactive Data Resource, SPIDR). Разработанный алгоритм анализа геомагнитных вариаций включает два этапа: предобработка измерительной информации (интерполирование, исключение тренда) и собственно анализ на основе Фурье и (или) вейвлет-преобразований.

2. Разработка алгоритма анализа геомагнитных данных

2.1. Данные для анализа

Как было сказано выше, данные для анализа предоставляет интерактивный ресурс по солнечно-земной физике (Space Phisics Interactive Data Resource, SPIDR, http://spidr.ngdc.noaa.gov/spidr), который проектировался как распределённая сеть синхронных баз данных, позволяющая получать, визуализировать и моделировать данные по космической погоде при помощи сети Интернет. В рассматриваемой автоматизированной системе пользователь может сделать запрос на сервер SPIDR, выбрав интервал времени, необходимые данные и список станций или спутников (параметры данных). Заполненные формы обрабатываются на сервере при помощи Java-сервлетов, которые распределяют запросы на серверы базы данных, сервер посылает результат пользователю [6]. На рис. 1 показана запись модуля вектора магнитных вариаций с минутным осреднением из обсерватории Новосибирск в период мощной магнитной бури 20 ноября 2003 года, полученная путём запроса к базе данных SPIDR.

2.2. Интерполяция исходных данных

Так как исходные данные (временной ряд) вследствие несовершенства используемой аппаратуры и каналов передачи данных содержат пропуски, необхо-

хЮ4

95 94

—

93

92

91

90

89

88

87 86

2003-09-29 2003-10-12 2003-10-25 2003-11-07 2003-11-20 2003-12-03 2003-12-16 2003-12-29

итс

Рис. 1. Исходные данные — вариации модуля вектора напряженности геомагнитого поля в период с октября по декабрь 2003 года, записанные на станции Новосибирск. По оси абсцисс отложены данные; по оси ординат — напряженность, нТ (http://spidr.ngdc.noaa.gov/ spidr/servlet/ImageServlet?width=800&height=4 00&representation= UNDEF&elem=F, 1,Geom,NVS_min&dateFrom=20031001&dateTo=20031231)

дима их интерполяция. По графику исходного ряда видно, что характеристики сигнала не остаются постоянными, изменяется его форма, частота и амплитуда. В связи с этим для аппроксимации ряда /п = /(¿п), заданного в моменты ¿0,г1,...,гп значениями /0, /1,..., /п, используется функция ^(¿п) (кубический сплайн), определённая на интервале г е [Ьг-1,1^] как

<Рг(г) = аг + ъг(г - г—) + ф - и-1)2 + ¿г(г - г-)3, (1)

где аг,Ъг ,сг, ¿г — коэффициенты сплайна, г = 1,2, 3 ,...,п — номер сплайна.

Коэффициенты сплайна находятся по методу подгонки, подробно описанному в работе [6]. Согласно этому методу коэффициент аг = /г-1, коэффициент Ъг = (/г - /г-1)Мг(сг+1 + 2сг)кг/3, где кг = ¿г - ¿г-1, 1 < г < п, ¿г = (сг+1 - сг)/(3кг). То есть остальные коэффициенты выражаются через коэффициент сг, а сам он находится по формуле:

сг = кг - /гсг+1, (2)

где 1г,кг — специальные подгоночные коэффициенты, которые определяются соотношениями:

г г - кг~ 1 кг-1 т кг

кг =-т-т—, 1г =-т-, (3)

вг - кг-1кг-1 вг - кг—11г—1

где вг = 2(кг-1 - кг), при этом к1 = 0,11 = 0.

В процессе выполнения прямого хода метода подгонки вычисляются все подгоночные коэффициенты 1г,кг по рекуррентному соотношению (3). Далее последовательно применяется формула (2), которая позволяет найти значения

всех коэффициентов вг, зная которые легко вычислить все остальные коэффициенты сплайнов.

После этого интерполирующую функцию можно рассчитать с помощью соотношения (1) в любой момент £ на интервале \Ь0,Ьп]. После интерполяции временной ряд задан значениями функции, следующими друг за другом с постоянным шагом д£:

2.3. Исключение линейного тренда

В настоящем алгоритме обработки предлагается использовать как вейвлет-, так и Фурье- преобразование, поскольку для некоторых интервалов временных рядов Фурье-анализ также вполне информативен. Поэтому, т.к. преобразование Фурье очень чувствительно к линейному тренду, а вейвлеты первого порядка не менее чувствительны, его необходимо исключить. Это можно сделать путём замены исходных значений ряда (4) первыми последовательными разностями.

Пусть yt = yt + et; yt = a + bt, где et; — отклонение от тренда.

Тогда at = yt - yt-1 = a + bt + et - (a + b(t - 1) + et-i) = b + (et - et-i).

Коэффициент b — константа, которая не зависит от времени. При наличии линейной тенденции остатки et достаточно малы и носят случайный характер. Поэтому первые разности уровней ряда at не зависят от переменной времени, их можно использовать для дальнейшего анализа.

Коэффициент b — среднее значение ряда определяется по формуле:

тогда центрированный ряд определяется: хг = Дг — Ь.

На рис. 2б показан центрированный и освобождённый от линейного тренда ряд. Однако так как линейный тренд очень мал, результат будет заметен только при внимательном сравнении периодограмм в области низких частот.

2.4. Построение периодограммы

Периодограмма вычисляется с использованием быстрого преобразования Фурье сигнала х длиной N по формуле:

при помощи алгоритма Кули-Тьюки [7]. После нахождения преобразования Фурье (5) периодограмма вычисляется следующим образом:

fk = f (tk),tk = kAt, k = 0,1,..., N - 1.

(4)

Xm = / y ~2n^N/2 i ^ / y x2ra+1 aNm/2, m = 0, 1,...,N - 1 (5)

n=0 n=0

Dj = —[(reXj)2 + (imXj)2], j = 0,1..., N/2

1

Таблица 1. Часто употребляемые вейвлет-базисы

Название вейвлет-базиса Определение Фурье-образ

DOG (derivative of a Gaussian — производная функции Гаусса) Г t2 1 ^m(t) = (—1)mdm e—-2 , где = dm [...] /dtm,m > 1 Фт(к) = m(ik)me—"2"

Морли -2 ф(Ь) = ejk0te—T- J(t) = 0(k)e— (k 2k0) , где 0(k) — функция Хевисайда

Пауля im V« = (m +1)(1 — a)m+i ФсО = 0(k)kme—k

При нечётном N вычисления продолжают до N/2 + 0,5. Отсчёты периодограммы соответствуют частотам ^ = Ау],] = 0,1,..., N/2, где Аи = ^^

2.5. Построение скалограммы

Для вычисления вейвлет-преобразования можно воспользоваться формулой:

1 м-1 / + -ъ

<«.">=псм) £ л -

где

N— 1 /t ,

2 V a

п{а,Ъ) =

k=0

Вейвлет-функция (6) вычисляется для значений аргументов а^ и Ъ, где i = = 0,1,...,Na — 1; Ъ = 0,1,...,Nb — 1 [8]. Используя (6), теперь можно ввести оценку локального спектра энергии — скалограмму [9]:

S(aj, bj) = |ЖА(аг,Ъг)|2.

Параметр a — масштаб, он определяет размер вейвлета. Параметр Ъ — сдвиг, он задаёт временную локализацию вейвлета, а символом * обозначена процедура комплексного сопряжения материнской (базисной) вейвлет-функции. В табл. 1 приведены наиболее часто употребляемые вейвлет-базисы [10].

Выбор функции ф(Ь), базисного вейвлета, как правило, определяется тем, какую информацию необходимо извлечь из сигнала. Каждый вейвлет имеет характерные особенности во временном и в частотном пространстве, поэтому иногда с помощью разных вейвлетов можно полнее выявить и подчеркнуть те или иные свойства анализируемого сигнала. Например, MHAT-вейвлет (Mexican

2

hat — «Мексиканская шляпа»), получается из DOG при m =2. Он имеет узкий энергетический спектр и два равных нулю момента, хорошо приспособлен для анализа сложных сигналов, так как коэффициенты W(a,b) зависят от малого интервала области частот вейвлета [11].

Другим часто применяемым базисом является вейвлет Пауля. Чем больше m, тем больше нулевых моментов имеет этот вейвлет. Вейвлет Морли представляет собой плоскую волну, модулированную гауссианом единичной ширины. С увеличением k0 возрастает частотная избирательность базиса, но ухудшается временное разрешение вейвлет-преобразования. На рис. 2г, е показаны скалограммы, выполненные вейвлетом Морле в разных диапазонах частот при k0 = 5. Необходимо отметить также, что у вейвлета Морле равен нулю только нулевой момент, поэтому он чувствителен к линейному тренду. Также своими особенностями обладает вейвлет Хаара.

В разработанной автоматизированной системе исследователю предоставляется возможность выбора любого из описанных базисных вейвлетов.

2.6. Построение скейлограммы

Оценка глобального спектра энергии иллюстрируется построением скейло-граммы

где N — число точек, по которому осуществляется осреднение [10]. Скей-лограмма является прямым аналогом сглаженной периодограммы в Фурье-анализе.

2.7. Построение вейвлет-скелета

Широкие контуры линий близких по частоте на рис. 2г мешают проследить за эволюцией их частот во времени. Чтобы нивелировать влияние контуров, необходимо построить вейвлет-скелет (скелетон), выделив те точки скалограм-мы, в которых она имеет максимумы по масштабу а и сдвигу Ь [1]:

где Бц = Б(аг,Ъг). На скелетонах рис. 3 хорошо заметно различие между периодическими компонентами сигнала, идущими параллельно оси времени (3а) и случайными, идущими параллельно оси периодов (3б).

(6)

Рис. 2. Скейлограмма а) и скалограмма б), выполненная вейвлетом Морле во всем возможном диапазоне частот; скейлограмма в) и скалограмма г) в диапазоне частот до 700 мин.; периодограмма д) центрированного и освобожденного от линейного тренда е)

a)

б)

Рис. 3. Скелетон, полученный при помощи вейвлета Морле во всем возможном диапазоне периодов а) и в диапазоне до 700 часов б)

3. Блок-схема алгоритма анализа геомагнитных данных

Как следует из вышеизложенного, алгоритм анализа геомагнитных вариаций на основе вейвлет-технологий можно представить в виде следующей диаграммы деятельности (Activity Diagramm), выполненной в среде UML-моделирования Umbrello UML Modeller (https://umbrello.kde.org/) (рис. 4).

Для того чтобы начать анализ, пользователь должен задать параметры данных — временной интервал, компоненту вектора напряжённости магнитного поля, название станции, интервал дискретизации и (или) имя файла. В поле эти данные автоматически загружаются из хранилища SPIDR или сразу импортируются из файла по желанию пользователя. Для хранения данных используется файл в формате CSV, что позволяет анализировать не только магнитовариационные данные, но и любые другие одномерные временные ряды, единственным условием является соблюдение требований формата CSV. Следующий шаг — интерполяция данных и построение графика интерполированного временного ряда. На основании полученного графика пользователь должен принять решение, содержит ли исследуемый сигнал линейный тренд, и в этом случае дать команду системе исключить его. Необходимым этапом анализа является построение скалограммы. Она даёт возможность оценить степень нестационарности сигнала. В случае если на исследуемом участке характеристики сигнала (амплитуда, период) остаются приблизительно постоянными, то для его анализа достаточно построить периодограмму, она должна отражать все значимые особенности сигнала. В противном случае требуется построение скейлограммы, которая даёт возможность определить, содержит ли сигнал близкие по частоте компоненты, что говорит о необходимости построения вейвлет-скелета. Некоторые исследователи полагают, что скелетон не только чётко и без лишних деталей визуализирует структуру анализируемого сигнала, но и содержит всю информацию о нем [11]. Однако если вейвлет-скелет не по-

Рис. 4. Алгоритм анализа геомагнитных вариаций на основе вейвлет-технологий

казывает нужной информации, это означает, что требуется изменить параметры вейвлет-преобразования — тип или порядок вейвлета, а, возможно, и диапазон исследуемых частот. Эта возможность также поддерживается алгоритмом.

Таким образом, процесс анализа полностью контролируется пользователем, он управляет всеми значимыми параметрами, от которых зависит результат работы системы, которая, в свою очередь, берет на себя наиболее трудоёмкие и рутинные операции — загрузку данных, построение периодограммы, скелетона и другие.

4. Общая структура автоматизированной системы

Программная часть автоматизированной системы написана на языке Python 3 и использует следующие внешние компоненты без изменения исходного кода (рис. 5). Numpy — реализует численные алгоритмы и работу с матрицами. Scipy — библиотека научных алгоритмов для Python, из этого модуля импортируются функция исключения линейного тренда и гамма-функция для вейвлет-базиса Пауля. Matplotlib 1.3 используется для построения графиков. Qt 4.8 применяется для построения графического интерфейса пользователя

Рис. 5. Структура автоматизированной системы

(GUI — Graphic User Interface) версии 4.8, так на момент написания программы версия Matplotlib 1.3 не поддерживает Qt 5 версии. PyQt — обеспечивает интерфейс для использования Qt 4 совместно с языком Python. Все эти модули поставляются с исходным кодом по той или иной свободной лицензии и доступны для использования в научных и коммерческих проектах.

Для упрощения разработки и сопровождения код программы разбит на несколько логических модулей. GUI — содержит формы и классы, отвечающие за построение графического интерфейса пользователя. Wawelet — компонент, содержащий вейвлеты, и класс, необходимый для выполнения непрерывного вейвлет-преобразования. Модуль Painter — отвечает за построение скалограм-мы, скейлограммы и периодограммы. Loader — содержит классы, необходимые для загрузки данных из сети Интернет и импорта их в программу. Этот компонент по интерфейсу CGI (Common Gateway Interface — общий интерфейс шлюза) связывается с ресурсом SPIDR, о котором подробно рассказано выше.

Отличительными особенностями программных средств является достаточно широкий набор базовых вейвлетов и регулируемых параметров преобразования, графический интерфейс пользователя, кроссплатформенность (благодаря использованию интерпретатора языка Python программные средства могут работать под управлением широкого круга операционные систем, не требуя при этом установки таких тяжеловесных и требующих приобретения лицензии сред, как, например, Matlab). Все функции, выполняющие длительные математические расчёты, выполняются в фоновом режиме и отделены от интерфейса пользователя, что позволило достичь хорошей отзывчивости пользовательского интерфейса.

5. Обсуждение результатов анализа

При помощи описанной выше автоматизированной системы анализа геомагнитных вариаций на основе вейвлет-технологий был исследован временной ряд полного вектора напряженности магнитного поля в период значительных

магнитных возмущений октября-ноября 2003 года по данным обсерватории Новосибирск (рис.1). В периодограмме (рис. 2а) ряда, освобождённого от линейного тренда (рис. 2б), можно видеть 4 значимые концентрации мощности с периодами 1440, 720, 120 и 30 минут. Периодограмма сильно изрезана, что свидетельствует о значительной стохастичности, присущей процессу геомагнитной активности. Спектральная полоса с периодом 1440 мин. (24 ч) — так называемая суточная вариация, является следствием вращения Земли вокруг своей оси. Вариация с периодом в 720 мин. называется лунной и возникает из-за приливного действия Луны. Причина остальных двух вариаций лежит в солнечно-земных взаимодействиях, однако, их природа и механизм возникновения до конца еще не известны. Кроме того, используя только периодограмму, невозможно выяснить, какие из этих вариаций характерны для периода магнитной бури, а какие проявляются и в магнитоспокойные дни. Следовательно, в алгоритме анализа геомагнитных вариаций оправдано использование вейвлет преобразования.

Из рис. 2г, е видно, что суточные и лунные вариации остаются постоянными. В небольших пределах изменяется их амплитуда. Напротив, 30 и 120-минутные вариации наблюдаются преимущественно в магнитовозмущен-ные дни, причем 30-минутные пульсации наблюдаются только во время значительных магнитных бурь, например, таких как 23 октября и 20 ноября 2003 года.

Заключение

Таким образом, адекватность и эффективность разработанных методики и программных средств показана на примере анализа магнитных вариаций в октябре-ноябре 2003 года. В частности, было выяснено, что для периода магнитных бурь характерны появления квазипериодических вариаций магнитного поля с периодами в 120 и 30 мин.

Разработанная автоматизированная система предназначена, в первую очередь, для просмотра и анализа вариаций одной из составляющих геомагнитного поля. Однако она может быть использована для исследования любых одномерных данных, представленных в формате CSV. Процесс анализа полностью контролируется пользователем, он управляет всеми значимыми параметрами, от которых зависит результат работы системы, которая, в свою очередь, берет на себя наиболее трудоёмкие и рутинные операции — загрузку данных, построение периодограммы, скелетона и другие. Отличительными особенностями разработанных программных средств, доступных по адресу https://github.com/abalckin/tesla, является достаточно широкий набор базовых вейвлетов и регулируемых параметров преобразования, графический интерфейс пользователя, кроссплатформенность.

Литература

1. Мандрикова О.В., Смирнов С.Э., Соловьев И.С. Метод определения Индекса геомагнитной активности на основе вейвлет-пакетов // Геомагнетизм и аэрономия. 2012. Том 52, № 1.

2. Иванов В.В. Вейвлет-анализ как инструмент анализа геомагнитных вариаций: теоретические основы и экспериментальные результаты // Доклады научной конференции Базы данных, инструменты и информационные основы полярных геофизических исследований. URL: www.izmiran.ru/polar2012/abstracts/tez_ ivanovvv.doc

3. Мандрикова О.В., Соловьев И.С. Вейвлет-технология анализа вариаций геомагнитного поля // Доклады 13-й Международной конференции Обработка и передача измерительной информации. URL: www.autex.spb.su/dspa/dspa2011/ dspa2011-6-2.doc

4. Витязев В.В. Вейвлет-анализ временных рядов: Учеб.пособие. СПб.: Изд-во С.-Петерб. ун-та, 2001.

5. Web Services Guide. REST API v1, v2 Space Physics Interactive Data Resource. URL: http://spidr.ngdc.noaa.gov/spidr/docs/SPIDR.REST.WSGuide. en.pdf (Version 3.0, May 2010).

6. Мудров А.Е. Численные методы для ПЭВМ на языках Бейсик, Фортран и Паскаль. Томск: МП «РАСКО», 1991.

7. Cooley, James W., John W.T. "An algorithm for the machine calculation of complex Fourier series". Math. Comput. 1965. 19: 297-301.

8. Foster G. Wavelets for period analysis of unevenly sampled time series // Astron. J. Vol. 112, N. 4 P. P. 1709-1729.

9. Scargle J.D. Wavelet and OtherMulti-resolution methods for Time Series Analysis. Statistical Challenges in Modern Astronomy II / Ed GJ. Babu and E.D. Ferigelson. N.Y. Springer-Verlag. P. 333-347.

10. A Practical Guide to Wavelet Analysis Christopher Torrence, Gilbert P. Compo. Amer. Meteor. Soc. / In Bulletin of the American Meteorological Society. Vol. 79, No. 1. (1 January 1998).

11. Астафьева Н.М. Вейвлет-анализ: основы теории и примеры применения // Успехи физических наук. Ноябрь 1996. Том 166, № 11,

DEVELOPMENT OF AUTOMATED SYSTEM FOR ANALYSIS GEOMAGNETIC VARIATIONS BASED ON WAVELET-TECHNOLOGIES

S.N. Verzunov

Senior Teacher, postgraduate, e-mail: [email protected] N.M. Lychenko

Doctor of Technical Sciences, e-mail: [email protected]

Kyrgyz-Russian Slavic University, Bishkek

Abstract. The methods and software tools for the analysis of geomagnetic variations on the basis of wavelet transform are developed. An analysis of magnetic variations demonstrated the efficacy of the developed technique. In particular, it was found that for the period of magnetic storms the appearance of quasi-periodic magnetic field variations with periods of 120 and 30 min is typical. Software tools are developed in Python 3, they are cross-platform and have a wide functional.

Keywords: wavelet transform, digital processing of non-stationary non-uniform time-series, analysis of geomagnetic variations, software.