УДК 519.6:551.510.413.5
МЕТОД ЗАПОЛНЕНИЯ ПРОПУСКОВ В ДАННЫХ КРИТИЧЕСКОЙ ЧАСТОТЫ НА ОСНОВЕ СОВМЕЩЕНИЯ МЕТОДА ВЕЙВЛЕТ-ПРЕОБРАЗОВАНИЯ И АВТОРЕГРЕССИОННЫХ МОДЕЛЕЙ
О.В. Мандрикова1, 2, Н.В. Глушкова1, 2
1Камчатский государственный технический университет, Петропавловск-Камчатский, 683003;
2Институт космофизических исследований и распространения радиоволн ДВО РАН, с. Паратунка, Камчатский край, 684034 e-mail: [email protected] e-mail: oksanam 1 @mail.kamchatka. ru
Описан метод многокомпонентного моделирования ионосферных данных, основанный на совмещении вейвлет-преобразования с моделями авторегрессии - проинтегрированного скользящего среднего, который позволяет выполнить анализ данных и построить прогноз. Для оценки метода и основанных на нем алгоритмов моделирования использовались данные критической частоты fF2, записанные на станции «Паратунка» (п-ов Камчатка). При моделировании данных выделены особенности, связанные с солнечной активностью, а также возникающие в периоды сильных землетрясений на Камчатке.
Ключевые слова: вейвлет-преобразование, модель авторегрессии - проинтегрированного скользящего среднего, критическая частота, аномалии.
Filling gaps method in critical frequency data on the based combination of the wavelet transforming method and the autoregressive models. O.V. Mandricova1, 2, N.V. Glushkova1,2 ^Kamchatka State Technical University, Petropavlovsk-Kamchatsky, 683003; institute of Cosmophysical Researches and Radio Wave Propagation, Paratunka, Kamchatka)
Method of multicomponent modeling of ionospheric data, based on combination of wavelet-transformation with autoregressive integrated moving average models is described. It allows to perform data analysis and build a prediction. To evaluate the method and modeling algorithms based on it we used data of critical frequency f0F2, recorded at the station “Paratunka” (Kamchatka Peninsula). While data modeling features associated with solar activity which appear in periods of strong earthquakes in Kamchatka are marked.
Key words: wavelet transform, autoregressive model, critical frequency, anomalies.
Введение
Одной из важных задач обработки и анализа ионосферных данных является задача контроля состояния ионосферы и автоматическое выделение и интерпретация аномалий, возникающих в ионосферной плазме в периоды возмущений [1]. Сложная структура регистрируемых ионосферных параметров и отсутствие формальной модели их описания делает поставленную задачу весьма сложной. Предметом данных исследований являются регистрируемые временные ряды критической частоты ионосферы f0F2. В данныхfF2 наблюдаются аномалии, длительностью от нескольких десятков минут до нескольких часов [2], возникающие в периоды повышенной солнечной активности. В сейсмоактивных областях эти аномалии могут возникать накануне сильных землетрясений [3]. Сложная структура аномалий не позволяет использовать для их выделения и анализа традиционные методы анализа временных рядов [4]. На основе совместного применения вейвлет-преобразования с моделями авторегрессии - проинтегрированного скользящего среднего (АРПСС) авторами предложен метод многокомпонентного моделирования данных f0F2, выполнения их прогноза и обнаружения аномальных особенностей.
В работе выполнено моделирование данных f0F2, регистрируемых на Камчатке, выделены характерные составляющие и изучена их структура. На основе прогнозирования данных выделены аномалии, возникающие в данных в периоды ионосферных возмущений. Статистический анализ показал, что эти аномалии могут возникать в периоды повышения солнечной активности и накануне сильных землетрясений.
Описание метода
Рассмотрим замкнутое пространство V = clos 2 (2j ф(2 t - к)): к е Z) с разрешением j = 0,
J L (R)
порожденное скэйлинг-функцией феХ2(Л) [5], где L2(R) - пространство Лебега. На основе кратномасштабных разложений до уровня m и основываясь на результатах работы [6], получим представление данных в виде:
— m
/о(0 = 2 (g [2 t] + е [2 t]) + f [2—mt ] , (нижний предел у), (1)
j=—1
где f [2—mt]е Vm,g\lJt] еЖ , W - пространство с разрешением j, порожденное вейвлет-
базисом ^ k(t) = 2J 2 Y(2Jt — к) ; детализирующие компоненты g[2jt] = 2dк*¥jk(t), где
к
коэффициенты разложения dJk = ^/, ^ ) ; аппроксимирующая компонента
f [ 2mt ] = 2 С— m, кф—m, к (t), где коэффициенты разложения С_ m, =( f, ф—m, к) ; e [ 2Jt ] - шумовые
к
составляющие (предполагается, что шум белый), разрешение j. В работе [7] показано, что подавление шума может быть выполнено на основе применения пороговой функции
| х,если|X > T
P (X) Чп ' (2)
10,если|X < T
для каждой детализирующей компоненты g [2t] в (1). Порог T = о2, о2 - дисперсия шума. Следуя работе [5], дисперсию шума с2 можно оценить на основе соотношения: о2 « Med (/,^ Л , где Med - медиана, j = —1, N - длина компоненты.
I' j, I Ю< к <N
Для идентификации модели применим следующие операции:
1. Используя традиционный подход [5], определим модели из класса моделей АРПСС для аппроксимации каждой из компонент f [2—mtJ и g [2J t], j = —1,m и оценим параметры моделей.
2. Выполним диагностические проверки полученных моделей компонент. Если погрешность модели компоненты удовлетворяет требованию, то будем считать, что данная компонента описывает характерные особенности структуры данных.
3. Объединим полученные представления в общую многокомпонентную конструкцию:
f(t) = 22 j (t j (t), (3)
J=1,M к=1,Nj
pj hJ
где j(t) = 2 YJл<к-I (t) — ’2^^jna%-„(t), - параметры авторегрессии компоненты с номером
l=1 n=1
= V^P%(0, t - коэффициенты разложения компоненты с номером М , рМ - порядок авторегрессионной модели компоненты с номером м , VVj - оператор взятия разности порядка Vj, 0^ - параметры скользящего среднего модели компоненты с номером м ,
- порядок модели скользящего среднего компоненты с номером М, а^ - ошибки модели компоненты с номером М , М - количество моделируемых компонент, NМ - длина компоненты с номером м , ЪМк - базис компоненты с номером м , j - разрешение.
Процедура выделения аномалий в компонентах разрешения j может быть построена на анализе остаточных ошибок моделей компонент при выполнении операции прогнозирования:
1. Прогнозирование значения s^ , q > 1 определяет прогноз s^k в момент t = к с
упреждением q. Значение s^ на основе модели (3) определяется как p^j Щ
j +q (t) =2 ^“"к +q— (t)——! 0>Г,.к + q—n (t )•
l=1
n=1
2. Остаточные ошибки компоненты модели с номером ц масштаба (—у) определяются как разность между прогнозными и фактическими значениями данных в момент времени ? = к + д :
aj ,к+д (?) = ^у',к+д,прогноз. (?) — ^ ,к+д, фактич.
(?).
3. Обнаружение аномалии в компоненте с номером ц масштаба (—у) можно выполнить на основе проверки условия:
1 иУ 2
°и, = 77£(а'‘+. (')) > Т-, ■ (4)
и г=1
где ТА - некоторое наперед заданное пороговое значение, определяющее наличие в данных аномалии масштаба (—у), UJ - длина окна наблюдения на масштабе (—у) .
Результаты экспериментов В экспериментах использовались часовые данные /Р2 за период 2001-2011 гг. С целью уменьшения погрешности получаемых результатов при моделировании были выбраны периоды с наименьшим количеством пропусков. В качестве базисных функций использовались ортогональные вейвлеты Добеши порядка 3, которые обеспечивают наименьшую погрешность аппроксимации данных /0р2 [8]. На основе кратномасштабных разложений до уровня т = 3 включительно было получено представление данных в виде (1). Уровень разложения определялся статистически и основывался на результатах работы [9]. Далее, на основе операции (2), были подавлены шумовые составляющие 1], результаты оценки дисперсии шума представлены
в табл. 1. Анализ результатов табл. 1 показывает, что уровень шума носит случайный характер.
Таблица 1
Результаты оценки дисперсии шума в данных /0П
Анализируемый период Дисперсия шума
2002 год 01.01.2002-31.03.2002 0,1929
01.04.2002-31.05.2002 0,1511
01.06.2002-31.08.2002 0,154
01.09.2002-30.11.2002 0,1851
01.12.2002-13.12.2002 0,1945
2006 год 01.01.2006-31.03.2006 0,1634
01.04.2006-31.05.2006 0,1269
01.06.2006-31.08.2006 0,1487
01.09.2006-30.11.2006 0,1676
01.12.2006-13.12.2006 0,186
2011 год 01.12.2010-31.03.2011 0,1674
01.04.2011-31.05.2011 0,2746
01.06.2011-31.08.2011 0,1511
01.09.2011-02.11.2011 0,1542
Полученные при идентификации близкие значения параметров моделей компонент позволили определить для данных/Р2 общую модель вида:
= Е Е <,* (0Ь% (і), <к (0 = (1 + 0,9В)2 (1 - Б)<$'к (0 + а$Л (і),
^=1>2 к=1, Щ
где Б1 Юзк(ґ) = Юз*_7(0, аЦк(ґ) - остаточные ошибки модели компоненты с номером з .
В табл. 2, на примере 2002 г., показаны параметры моделей компонент, полученные при идентификации для данных различных сезонов.
Таблица 2
_________________________________Результаты моделирования данных /ц/2________________________________
Анализируемый период Параметры моделей аппроксимирующих компонент Параметры моделей детализирующих компонент
Значение первого параметра Значение второго параметра Значение первого параметра Значение второго параметра
01.01.2002-31.03.2002 -0,9875 -0,9918 -0,9942 -0,9689
01.04.2002-31.05.2002 -0,8424 -0,8451 -0,998 -0,9068
01.06.2002-31.08.2002 -0,7019 -0,685 -0,9635 -0,9153
01.09.2002-30.11.2002 -0,9643 -0,9679 -1,021 -0,9591
Результаты статистики показали, что частота появления аномалий и их интенсивность
зависят от уровня солнечной активности (рис. 1, 2). Данный факт согласуется с результатами работы [1], где представлены исследования ионосферных возмущений в верхней атмосфере Земли на основе ОР8-мониторинга. Сопоставление полученных результатов с данными каталога землетрясений показывает, что в сейсмически спокойные периоды времени возрастание ошибок моделей наблюдается во время магнитных бурь, особенно для аппроксимирующих компонент. В периоды повышенной сейсмической активности (рис. 1, 3) характер процесса меняется, и существенное увеличение ошибок моделей наблюдается в периоды возникновения сильных землетрясений. Возрастание ошибок весной (март, апрель 2002 г. и 2011 г.), возможно, связано с переходными процессами в ионосфере, характерными для данного периода времени.
Рис. 1. Результат моделирования ионосферных данных за период 01.01.2002-4.12.2002. Стрелками отмечены моменты возникновения землетрясений
Рис. 2. Результат моделирования ионосферных данных за период 15.01.2011-25.06.2011.
Стрелками отмечены моменты возникновения землетрясений В табл. 3, на примере 2011 г., также показаны параметры моделей компонент, полученные при идентификации для данных различных сезонов после их вейвлет-восстановления. Результаты моделирования каждой компоненты и процесс их совмещения в общую многокомпонентную конструкцию показаны на рис. 3, 4. Анализ полученных результатов подтверждает эффективность предлагаемого метода.
Таблица 3
Результаты моделирования восстановленных компонент
Анализируемый период Восстановленная аппроксимирующая компонента Восстановленная детализирующая компонента 3-го уровня разложения Восстано детализи компо 2-го уровня вленная эующая нента разложения
Первый параметр Второй параметр Первый параметр Второй параметр Первый параметр Второй параметр
01.12.2010-31.05.2011 1,009 -0,2668 0,8176 -0,3443 0,3353 -0,6614
01.04.2011-31.05.2011 1,011 -0,266 0,8169 -0,3446 0,3193 -0,6778
15.01 29.01 12.02 26.02 12.03 26.03 09.04 23.04 07.05 21.05 04.06 18.06
Рис. 3. Результаты моделирования ионосферных данных за период 15.01.2011-25.06.2011. Сигнал критической частоты (черным), прогнозные значения (серым). Стрелками отмечены моменты возникновения землетрясений
Рис. 4. Результаты моделирования ионосферных данных за период 15.01.2011-29.01.2011. Сигнал критической частоты (черным), прогнозные значения (серым)
Выводы. В работе предложен метод многокомпонентного моделирования данных ионосферы, основанный на совмещении конструкции кратномасштабного анализа и авторегрессионных моделей. Выполнено моделирование и анализ данных критической частоты ионосферы f0F2, регистрируемых на Камчатке. В результате моделирования были изучены аппроксимирующая и детализирующая компоненты данных и выделены аномалии, возникающие в периоды повышенной солнечной или сейсмической активности на Камчатке.
Работа поддержана грантом Президента Российской Федерации МД-2199.2011.9, грантом РФФИ - ДВО РАН №11-07-98514-р_восток_а и грантом «У.М.Н.И.К.» - № 9633р/14207 от 30.08.2011.
Данные сейсмического каталога любезно предоставлены Камчатским филиалом геофизической службы РАН (г. Петропавловск-Камчатский).
Литература
1. Афрамович Э.Л., Перевалова Н.П. GPS-мониторинг и верхней атмосферы Земли. -Иркутск: ГУ НЦ РВХ ВСНЦ СО РАМН, 2006. - 480 с.
2. Дёмин М.Г. Ионосфера Земли. Плазменная гелиогеофизика. - М.: Физматлит, 2008. -Т. II. - С. 92-163.
3. Липеровская Е.В., Липеровский В.А., Похотелов О.А. О возмущениях в F-области ионосферы перед землетрясениями // Геофизические исследования. - 2006. - № 6. - С. 51-58.
4. Марпл.-мл. С.Л. Цифровой спектральный анализ и его приложения: Пер. с англ. -М.: Мир, 1990.
5. MallatS. A Wavelet tour of signal processing [пер. с анг.]. - М.: Мир, 2005. - 671 с.
9
6. Мандрикова О.В., Глушкова Н.В. Многокомпонентное моделирование и анализ аппроксимирующих компонент критической частоты ^2 на основе вейвлет-преобразования и моделей авторегрессии // Междунар. конф. по мягким вычислениям и измерениям ^СМ'2011). -СПб.: Изд-во СПбГЭТУ «ЛЭТИ», 2011. - Т. 2. - С. 139-143.
7. Мандрикова О.В., Горева Т.С. Метод идентификации структурных компонентов сложного природного сигнала на основе вейвлет-пакетов // Цифровая обработка сигналов. - М., 2010. - № 1.
- С. 45-50.
8. Мандрикова О.В., Полозов Ю.А. Критерии выбора вейвлет-функции в задачах аппроксимации природных временных рядов сложной структуры // Информационные технологии.
М., 2012. - № 1. - С. 31-36.
9. Мандрикова О.В., Глушкова Н.В. Метод моделирования данных критической частоты на основе совмещения вейвлет-преобразования и моделей авторегрессии - проинтегрированного скользящего среднего // Научные ведомости Белгородского государственного университета. -Белгород, 2011.- № 19. - С. 15-21.