CC BY
24
УДК 519.6: 551.510.413.5
ПРОГНОЗИРОВАНИЕ ДАННЫХ КРИТИЧЕСКОЙ ЧАСТОТЫ ИОНОСФЕРЫ НА ОСНОВЕ МНОГОКОМПОНЕНТНОЙ МОДЕЛИ
О.В. Мандрикова1,2, Н.В. Глушкова1,2
1Камчатский государственный технический университет, Петропавловск-Камчатский, 683003;
2Институт космофизических исследований и распространения радиоволн ДВО РАН, с. Паратунка, Камчатский край,684034
e-mail: nv. slushkova@yandex. ru e-mail:oksanam1 @mail.kamchatka. ru
Исследование направлено на разработку технологий и автоматизированных систем по изучению спокойного (фонового) хода ионосферы и выявлению аномалий в периоды ионосферных возмущений. Представлен метод моделирования и прогнозирования ионосферных данных, основанный на совмещении вейвлет-преобразования с моделями авторегрессии - проинтегрированного скользящего среднего, позволяющий выполнить анализ данных и построить их прогноз. На основе операции прогнозирования данных может быть решена задача заполнения в них пропусков с учетом суточного и сезонного хода. Анализ остаточных ошибок прогнозирования позволил выделить особенности, связанные с солнечной активностью, а также возникающие в периоды сильных землетрясений на Камчатке. Для оценки метода использовались данные критической частоты ионосферы, записанные на станции «Паратунка» (п-ов Камчатка, ИКИР ДВО РАН) за период 1969-2011 гг.
Ключевые слова: анализ аномалий, частоты ионосферы, многокомпонентное моделирование, вейлет-преобразование, авторегрессия.
Data forecasting of ionosphere critical frequency based on the on the basis of multicomponent model.
O.V. Mandricova1,2, N.V. Glushkova1,2 ^Kamchatka State Technical University, Petropavlovsk-Kamchatsky, 683003; 2Institute of Cosmophysical Researches and Radio Wave Propagation, Paratunka, Kamchatka, 684034)
The research is aimed at the development of technologies, automated systems for studing calm (background) ionosphere motion and detecting anomalies during ionospheric disturbances. The paper describes simulation technique and forecast method of ionospheric data which is based on combination of wavelet transform with autoregressive models - integrated sliding average; this method allows us to carry out data analysis and to make prognostication. Based on the data prognostication the problem of filling in gaps in them can be solved taking into account daily and seasonal trends. Analysis of the residual prediction error allowed to identify some peculiarities associated with solar activity and arrising during strong earthquakes in Kamchatka. To evaluate the method, ionosphere critical frequency data registered at “Paratunka” (Kamchatka peninsular, IKIR FEB RAS) were used in the period from 1969 to 2011.
Key words: analysis of anomalies, ionosphere frequency, multicomponent simulation, wavelet transformation, autoregressive.
Введение
Предметом исследований являются данные критической частоты ^2-слоя ионосферы, представленные в виде временного ряда. Одной из важных задач обработки и анализа ионосферных данных является задача контроля состояния ионосферы и автоматическое выделение и интерпретация аномалий, возникающих в периоды ионосферных возмущений [1-3]. Основные сложности решения поставленных задач связаны с исходной априорной неопределенностью знаний о структуре ионосферных данных, условий их формирования, а также отсутствием формальной модели их описания. Динамический режим ионосферы определяется распределением ее параметров по высоте, плотностью атмосферы, химическим составом и спектральными характеристиками солнечного излучения [1-3]. На фоне регулярных суточных и сезонных изменений в данных fF2 наблюдаются аномалии, длительностью от нескольких десятков минут до нескольких часов, обусловленные повышением солнечной и сейсмической активностями [3, 4]. Традиционные методы моделирования временных рядов не позволяют адекватно описывать характерные изменения параметров ионосферы и не обеспечивают эффективное выделение и интерпретацию аномалий [4-
6]. Предложенный в данной статье метод основан на применении вейвлет-преобразования [7, 8]. Ввиду наличия быстрых алгоритмов преобразования данных и большого разнообразия базисных функций с компактными носителями этот аппарат позволяет детально изучить внутреннюю структуру данных и выделить локальные особенности различной формы и длительности [7]. На основе совместного применения вейвлет-преобразования с моделями авторегрессии -проинтегрированного скользящего среднего (АРПСС), авторами статьи предложен метод моделирования и прогнозирования данных /ор2. Идентификация многокомпонентных моделей основана на применении операции кратномасштабного анализа и представлении данных в виде аппроксимирующей и детализирующих компонент и аппроксимации полученных компонент моделями АРПСС.
Построены многокомпонентные модели критической частоты ионосферы ^2-слоя по данным Камчатки, выполнено прогнозирование и анализ данных. Выделены аномалии, возникающие в периоды ионосферных возмущений. Анализ полученных результатов показал, что данные аномалии возникают в периоды повышенной солнечной активности и могут наблюдаться накануне сильных землетрясений на Камчатке.
Многокомпонентное моделирование, прогнозирование и анализ данных на основе совмещения кратномасштабного анализа и моделей АРПСС
Структура разложения пространства Лебега Ь2(Я), порожденная ортогональным вейвлетом ¥ е Ь2 (К), имеет вид [7, 8]:
®
Ь2(К) = £:=...Ф W-1 Ф Ж0 Ф Ф ..., (1)
1е2
где WJ := с1о^^;п е 2), '¥] п = 21 /2 ¥(2Ч - п), Ъ - множество целых чисел.
Функция / при этом представляется в виде суммы компонент:
V/ еЬда ! /(0 =......+ g-^) + g0(0 + gl(t) +., gj еWj, 1 е 2 . (2)
Каждая компонента g ■ в (2) имеет единственное представление в виде вейвлет-ряда,
последовательность коэффициентов которого дает локализованную спектральную информацию от / в 1 -й октаве (частотном диапазоне): gj (^ = £djn¥1п ^),где ¥'=|¥уи^2 - базис
п
пространства ^. Коэффициенты определяются из соотношения: = ^/,.
Коэффициенты ^ п будем рассматривать как результат отображения / в пространство с разрешением 1 .
Используя разложение Ь2 (К) в (1), мы имеем последовательность вложенных друг в друга замкнутых подпространств V, 1 е 2 , определенных формулой [8]:
VI =...Ф Wj-2 Ф Wj- , (3)
при этом с1о^2 (и 1е2 Vj ) = Ь (К). Порождает пространства (3) скэйлинг-функция фе Ь(К) [7, 8]. На основе соотношений (1), (3) получаем следующее разложение пространства Ь (К) :
Ь (К) = V Ф WJ Ф WJ+1 Ф ... (4)
Пусть /■ - некоторая проекция / на V- для фиксированного 1 е 2 . Пространство V- в этом случае рассматривается в качестве пространства выборки и /■ - измерения / в V.. Как доказано в [8],
Vj = Wj- + Vj- = .... = Wj- + .... + Wj-т + Vj-т Vm, то /■ имеет единственное кратномасштабное разложение (КМА) до уровня т :
fj (t) = gj-1(t) + gj-2(t) + ... + gj - m (t) + fj-m (t) .
(5)
где gj (t) eWj , fj -m (t) -m .
Тогда введя предположение, что регистрируемые данные содержат шумовую составляющую e(t), и используя соотношение (5) для f . (j = 0), исходную непрерывную функцию ft). полученную по совокупности регистрируемых дискретных данных, представим в виде
-m
f0 (t) = I (gj (t) + Є (t” + f-m (t) . (б)
j=-1
где /_т (1) е ^ ^ - пространство с разрешением порожденное вейвлет-базисом
¥],к(О = 2 /2*¥(2Ч - к) , аппроксимирующая компонента /-т (1) = ^ с_т>кф-т>к (1), где
коэффициенты разложения с_тк =( f, ф_тк ), детализирующие компоненты g (t) = I j (‘),
k
где коэффициенты разложения djk = ^ f, ; e (t) - шумовые составляющие (предполагается,
что шум белый), разрешение j (соответствует масштабу (—j)).
В работе [9] показано, что, имея представление данных в виде (6), подавление шума может быть выполнено на основе применения пороговой функции
| х, еслиX ^ T
PT (X) = 1п 1т (?)
10, еслщ X < T
для каждой детализирующей компоненты g. (t) , где порог T = О2 , О2 - дисперсия шума. Следуя работам [7, 9], дисперсию шума О можно оценить на основе соотношения:
а 2 * Med( f, Л ,
IV \0<k< N
где Med - медиана, j = — 1, N - длина компоненты.
Для выделения компонент конструкции (6), описывающих характерные особенности данных, и идентификации их параметров применим следующие операции.
1. Отобразим данные критической частоты в вейвлет-пространство до уровня m на основе операции (5).
2. Для каждой детализирующей компоненты gj (t) = I j *j. (t) на основе операции (7)
k
выполним подавление шума.
3. Восстановим каждую из полученных компонент f [2—mt] и g[2Jt], j = —1,mдо масштаба
j = 0, получим восстановленные компоненты вида: f0 (t) = I co,k Ф„ (t), g0j 40 = I dJXk (t),
k k
где j - разрешение компоненты до восстановления.
4. Используя традиционный подход [10], определим модели из класса моделей АРПСС для
аппроксимации каждой из полученных восстановленных компонент f0 (t) = I с0,kФо,k (t) и
к
й' '(1) = Е <Хк (1).
к
5. Выполним диагностические проверки полученных моделей компонент. Если диагностические проверки модели компоненты подтверждают ее адекватность, то будем считать, что модель компоненты готова к использованию и данная компонента является характерной.
6. Объединим модели выделенных характерных компонент в общую многокомпонентную конструкцию. Получим многокомпонентную модель вида:
к
/м) = X Хіо) ь%(і),
(8)
;=і,м »=і, щ
параметры авторегрессии компоненты с номером ц, юцк ^ ) = VV 1 Рцк ), р**. = су к , вц к = djk,ц = 2,М, - порядок авторегрессионной модели компоненты с номером ц, Л;ц,9^ -
порядок модели и параметры скользящего среднего модели компоненты с номером ц, аЦк -остаточные ошибки модели компоненты с номером ц, М - количество моделируемых
Щк = к, ц = 2,М - вейвлет-базис компоненты с номером ц, і - разрешение.
Если данные содержат аномалию, то произойдет изменение их структуры. Поэтому процедура выделения аномалий в компонентах разрешения і может быть построена на обработке остаточных
ошибок аЦк моделей компонент при выполнении операции прогнозирования. Данная процедура предполагает следующие операции:
1. Прогнозирование значения *, д > 1 определяет прогноз * цк в момент I = к с
упреждением д. Значение * Ц к+ч на основе модели (8) определяется как
2. Остаточные ошибки компоненты модели с номером р разрешения определяются как разность между прогнозными и фактическими значениями данных в момент времени t = к + q :
3. Обнаружение аномалии в компоненте с номером р разрешения 1 можно выполнить на основе проверки условия:
где Г - некоторое наперед заданное пороговое значение, определяющее наличие в данных аномалии разрешения 1, Ц - длина окна наблюдения для разрешения 1.
В экспериментах использовались часовые данные критической частоты f0F2 за период 1969-2011 гг. Учитывая сезонный характер ионосферного процесса, данные предварительно были разделены на сезоны и анализировались отдельно. В данной работе представлены результаты прогнозирования данных зимнего и летнего периодов времени. В силу физических и технических причин в регистрируемых данных содержатся пропуски, что существенно затрудняет процесс их моделирования и анализа. С целью уменьшения погрешности получаемых результатов были выбраны временные периоды, когда измерения у^2 велись без существенных пропусков.
В качестве базисных функций использовались ортогональные вейвлеты Добеши порядка 3, которые, как показала статистика, обеспечивают наименьшую погрешность аппроксимации данных у^2 [4, 11]. На основе кратномасштабных разложений до уровня т = 3 включительно было получено представление данных в виде
характерных компонент, Nц - длина компоненты с номером ц, Ь1 к = ф . к- скэйлинг-функция,
і=і
П=1
(9)
Результаты экспериментов
-3
/„(') =Х (gj (I) + е (I)) + /-з(1) ,
і=-1
где gJ ) = к(1) - детализирующие компоненты разложения, ^к =^ /.
к
ы<)=1 С3. £ф_3к(1) - аппроксимирующая компонента, с_тк = ^/.§_ткв;.) - шумовые
к
составляющие.
Уровень разложения определялся статистически и основывался на результатах работы [4], в которой показано, что исходные данные f0F2, аппроксимирующие компоненты 1-го и 2-го уровней разложения не могут быть аппроксимированы моделью АРПСС (наблюдалась существенная автокорреляция остатков модели). Наилучшие результаты при диагностике моделей выделенных компонент были получены для уровня разложения т = 3 .
Далее на основе операции (7) были подавлены шумовые составляющие в. (?), анализ
результатов оценки дисперсии шума показал, что уровень шума носит случайный характер [5].
В табл. 1 показаны параметры моделей восстановленных компонент, полученные на основе описанных выше операций (1)-(5). Результаты моделирования данных зимнего и летнего периодов времени показали, что все модели АРПСС восстановленных компонент 3-го уровня разложения имеют второй порядок. Близкие значения параметров зимнего и летнего периодов позволяют предположить, что параметры восстановленных компонент не зависят от сезона.
Таблица 1
Параметры моделей восстановленных компонент
Анализируемый период Аппроксимирующая компонента Детализирующая компонента 3-го уровня разложения Детализирующая компонента 2-го уровня разложения
первый параметр второй параметр первый параметр второй параметр первый параметр второй параметр
зимний период 04.01.70-05.02.70 1,01 -0,27 0,82 -0,34 0,38 -0,61
04.02.75-25.02.75 1,01 -0,27 0,83 -0,35 0,39 -0,63
07.12.79-22.12.79 1,01 -0,27 0,83 -0,35 0,4 -0,62
23.01.81-06.02.81 1,01 -0,27 0,76 -0,36 0,22 -0,78
07.02.83-23.02.83 1,01 -0,27 0,83 -0,34 0,33 -0,68
01.01.91-26.01.91 1,01 -0,27 0,81 -0,35 0,31 -0,69
12.01.92-05.02.92 1,01 -0,27 0,81 -0,35 0,21 -0,79
01.12.00-22.12.00 1,01 -0,27 0,83 -0,34 0,38 -0,61
24.01.02-25.02.02 1,01 -0,27 0,8 -0,35 0,34 -0,65
21.12.03-03.02.04 1,01 -0,27 0,82 -0,34 0,39 -0,6
10.12.10-31.12.10 1,01 -0,27 0,83 -0,34 0,3 -0,68
08.02.11-27.02.11 1,01 -0,27 0,81 -0,35 0,44 -0,47
летний период 31.05-02.07.1979 1,03 -0,25 0,79 -0,35 0,30 -0,71
22.05.-18.06.1984 1,01 -0,26 0,82 -0,34 0,40 -0,61
05.08-31.08.1987 1,01 -0,27 0,83 -0,34 0,42 -0,59
02.06-29.06.1989 1,02 -0,26 0,80 -0,35 0,39 -0,61
15.08-31.08.1992 1,01 -0,26 0,83 -0,34 0,35 -0,65
30.07-21.08.1999 1,02 -0,26 0,82 -0,34 0,43 -0,58
13.06-28.06.2000 1,01 -0,26 0,76 -0,36 0,43 -0,58
30.07.-17.08.2002 1,02 -0,26 0,83 -0,34 0,36 -0,64
27.06-10.07.2005 1,01 -0,29 0,83 -0,34 0,44 -0,57
Процедура выделения аномалий в полученных компонентах fF2 была основана на операциях (1)-(3) и определении дисперсии ошибки прогноза (величина Пи в соотношении (9)). В качестве
примера на рис. 1 представлены результаты прогнозирования каждой восстановленной компоненты с шагом упреждения q = 1,2,3,4,5 и процесс их совмещения в общую многокомпонентную конструкцию для временного периода 09.02.2011-27.02.2011, а также результаты расчета дисперсии ошибки прогноза в скользящем временном окне, равном 24 ч. Анализ графиков показал, что накануне сейсмического события 20.02 наблюдается увеличение дисперсии ошибок Ои .
Рис. 1. Результаты прогнозирования данных f0F2: черная линия - исходные данные f0F2, серая линия - прогноз восстановленных компонент, (а) - прогнозирование аппроксимирующей компоненты, (б) - прогнозирование аппроксимирующей и детализирующей компонент 3-го уровня разложения, (в) - прогнозирование аппроксимирующей компоненты и детализирующих компонент 3-го и 2-го уровня разложения
Выводы. В работе представлен метод моделирования и прогнозирования данных критической частоты ионосферы, основанный на совмещении конструкции кратномасштабного анализа и авторегрессионных моделей. Идентифицированы параметры многокомпонентной модели зимнего и летнего периодов времени данных foF2■ На основе полученной модели выполнено прогнозирование и анализ данных и изучена их внутренняя структура. Анализ результатов операции прогнозирования позволил обнаружить аномальные эффекты, которые обусловлены солнечной активностью, а также процессами в литосфере, формирующимися в периоды повышенной сейсмической активности (анализировались события энергетического класса с к>12 в
радиусе R~200 км от Петропавловска-Камчатского).
Работа поддержана грантом грантом РФФИ - ДВО РАН №11-07-98514-р_восток_а и грантом «У.М.Н.И.К.» - № 9633р/14207 от 30.08.2011.
Данные сейсмического каталога любезно предоставлены Камчатским филиалом геофизической службы РАН (г. Петропавловск-Камчатский).
Литература
1. Афраймович Э.Л., Перевалова Н.П. GPS-мониторинг верхней атмосферы Земли. -Иркутск: ГУ НУ РВХ ВСНЦ СО РАМН, 2006. - 480 с.
2. Дёмин М.Г. Ионосфера Земли. Плазменная гелиогеофизика. - М.: Физматлит, 2008. -T.II.
- С. 92-163.
3. Липеровская Е.В., Липеровский В.А., Похотелов О.А. О возмущениях в F-области ионосферы перед землетрясениями // Геофизические исследования. - 2006. - № 6. - С. 51-58.
4. Мандрикова О.В., Глушкова Н.В. Метод моделирования данных критической частоты на основе совмещения вейвлет-преобразования и моделей авторегрессии - проинтегрированного скользящего среднего // Научные ведомости Белгородского государственного университета. -2011. - № 19.- С. 59-63.
5. Глушкова Н.В., Мандрикова О.В. Обнаружение и анализ аномалий в данных критической частоты ионосферы на основе совмещения вейвлет-преобразования и авторегрессионных моделей // Международная молодежная конференция «Прикладная математика, управление и информатика» (ПМУИ-2012). - Белгород: Ид «Белгород», 2012. - Т. 1. - С. 355-360.
6. Март.-мл. С.Л. Цифровой спектральный анализ и его приложения: Пер. с англ. - М.: Мир, 1990. - 547 с.
7. MallatS. A Wavelet tour of signal processing [пер. санг.]. - М.: Мир, 2005. - 671 с.
8. Чуи К. Введение в вейвлеты: Пер. с англ. - М.: Мир. - 2001. - 412 с.
9. Мандрикова О.В., Горева Т.С. Метод идентификации структурных компонентов сложного природного сигнала на основе вейвлет-пакетов // Цифровая обработка сигналов. - 2010. -№ 1. - С. 45-50.
10. Бокс Дж., Дженкинс Г. Анализ временных рядов: Прогноз и управление. - М.: Мир, 1974. - 604 с.
11. Мандрикова О. В., Полозов Ю.А. Критерии выбора вейвлет-функции в задачах аппроксимации природных временных рядов сложной структуры // Информационные технологии.
- 2012. - № 1. - С. 31-36.
