Using classification methods in forecasting the level of geomagnetic field disturbance based on the Kp-Index

Толық мәтін

1. ВВЕДЕНИЕ

Геомагнитные возмущения – это один из важнейших факторов космической погоды, роль которого будет возрастать с развитием космической отрасли и глобальной цифровой индустрии как на Земле, так и в околоземном космическом пространстве [Daglis, 2001; Cole, 2003; Vassiliadis, 2007; Kudela, 2013; Schrijver et al., 2015; McGranaghan et al., 2021]. На Земле магнитные бури могут вызывать сбои в работе радиосвязи, трубопроводов, линий электропередач и электрических сетей, а также вызывать проблемы со здоровьем человека [Лазутин, 2012; Qiu et al., 2015; Белаховский и др., 2023].

Помимо этого, геомагнитные возмущения опосредованно оказывают заметное влияние на состояние околоземного космического пространства за счет того, что после примерно половины магнитных бурь поток релятивистских электронов внешнего радиационного пояса Земли (РПЗ) увеличивается на порядок и более (см., например, [Kataoka et al., 2008; Мягкова и др., 2013]). Резкий рост потоков электронов внешнего РПЗ несет существенную опасность для аппаратуры, находящейся на борту космических аппаратов, поскольку экстремальные потоки электронов могут привести к сбоям в электронных микросхемах [Белов и др., 2004; Романова и др., 2005; Iucci et al., 2005].

Учитывая сказанное выше, задача прогнозирования геомагнитных возмущений имеет очень большое практическое значение.

Геомагнитную активность принято характеризовать специальными индексами. Одним из наиболее распространенных геомагнитных индексов является Кр-индекс, впервые введенный Юлиусом Бартельсом в 1939 году [Bartels et al, 1939; Bartels, 1949].

Как известно, Кр – планетарный индекс, характеризующий глобальную возмущенность магнитного поля Земли в трехчасовом интервале времени. Индекс Kp определяется как среднее значение уровней возмущения двух горизонтальных компонент геомагнитного поля, наблюдаемых на 13 магнитных обсерваториях, расположенных в субавроральной зоне между 48 и 63 градусами северной и южной геомагнитных широт. Для определения Kp-индекса используются стандартизованные значения локальных K-индексов этих 13 обсерваторий.

Также известно, что основной причиной, которая вызывает возмущения в магнитосфере Земли, являются процессы, происходящие на Солнце и в гелиосфере [Akasofu et al., 1974]. Источниками возмущений магнитосферы Земли могут служить как высокоскоростные потоки солнечного ветра (ВСП СВ) из корональных дыр, так и корональные выбросы массы (КВМ), сопровождающие солнечные вспышки и достигающие орбиты Земли (например, [Ермолаев и др., 2009]). ВСП СВ играют основную роль в периоды минимумов солнечной активности, когда нет ни вспышек, ни связанных с ними мощных КВМ. Существование связи интенсивности геомагнитных возмущений (а, следовательно, и амплитуды Кр-индекса) с параметрами солнечного ветра (СВ), в частности с его скоростью, а также с параметрами межпланетного магнитного поля (ММП), известна достаточно давно [Nishida, 1972]. Детально данная связь исследовалась в работах [Elliott et al., 2013; Boroyev et al., 2020] и др.

Кр-индекс имеет самую долговременную историю наблюдения после aa-индекса (первый доступен с 1932, а второй с 1868 г.). Это дает возможность исследовать связь геомагнитной активности с процессами в межпланетном пространстве, солнечном ветре и магнитосфере Земли и прогнозировать Кр на основе этой связи (например, [Wang et al., 2015]). Также наличие столь долговременной истории измерений Кр-индекса позволяет предположить, что другим достаточно перспективным подходом к его прогнозированию является применение методов машинного обучения, в том числе искусственных нейронных сетей [Boberg et al., 2000; Bala et al., 2012; Ji, E. Y. et al, 2013; Tan et al., 2018].

Авторами данной работы – сотрудниками лаборатории адаптивных методов обработки данных НИИЯФ МГУ – также достаточно давно ведутся работы по прогнозированию геомагнитных возмущений, а именно, амплитуды Dst-индекса, методами машинного обучения. Было показано, что лучшее качество прогноза Dst-индекса достигается при построении нейросетевой модели, использующей в качестве входных данных как историю Dst-индекса, так и параметры СВ (скорость) и ММП (компонента Bz) [Dolenko et al., 2005]. В относительно недавних работах лаборатории [Myagkova et al., 2017; Ефиторов и др., 2018; Мягкова и др., 2021] каждый пример также содержал среднечасовые значения нескольких основных параметров СВ и ММП и самого прогнозируемого параметра (в данном случае индекса Dst) за последние 24 ч, что позволило улучшить качество прогноза. Использование подобного подхода стало возможным только в последнее время, когда накопились достаточно длинные однородные временные ряды спутниковых измерений параметров СВ и ММП, полученные в экспериментах на космическом аппарате (КА) АСЕ (подробнее о входных данных см. ниже). Но поскольку в большинстве прикладных моделей, касающихся явлений космической погоды, требуется прогноз не Dst, а Кр- индекса, было принято решение распространить имеющиеся наработки на задачу прогнозирования Кр-индекса.

В настоящей работе исследуется задача классификации – прогнозирования категорий геомагнитных возмущений по уровню индекса Kp на разных горизонтах. Выбранные категории условно соответствуют состояниям “невозмущенная магнитосфера” (Kp от 0 до 2-), “слабые геомагнитные возмущения” (Kp от 2 до 3+) и “умеренные и сильные геомагнитные возмущения” (Kp от 4-).

В работе [Gadzhiev I. et al, 2023] исследовалась возможность классификации категорий индекса Kp с другим разбиением по классам, что привело к более выраженной несбалансированности классов – в результате качество работы моделей на последних горизонтах для наиболее важных для прогнозирования небольших классов получалось неудовлетворительным из-за недостатка данных для обучения. Поэтому было принято решение использовать указанные выше границы категорий.

Для классификации используются предыдущие значения Kp и двенадцати других физических величин, о которых будет подробно рассказано в следующем разделе. Классификация индекса Kp выполняется на восьми горизонтах с интервалом три часа (то есть первый прогноз категории выполняется на три часа вперед, а последний на сутки вперед). Очевидно, что сложность задачи классификации увеличивается с увеличением горизонта.

Таким образом, цель настоящей работы состоит в выработке метода построения модели машинного обучения, наилучшим образом выполняющей поставленную задачу классификации состояний магнитосферы Земли по категориям геомагнитной возмущенности, определяемым по уровню индекса Kp. Для достижения этой цели решаются следующие задачи: сравнение результатов различных методов машинного обучения; исследование возможности использования специальных методик для преодоления несбалансированности категорий и выбора оптимальных гиперпараметров моделей; сравнение результатов, получаемых при включении в обучающую выборку данных 23 и 24 циклов солнечной активности либо только 24 цикла; анализ входных признаков, оказавшихся наиболее существенными с точки зрения метода-победителя.

2. ИСТОЧНИКИ ДАННЫХ ДЛЯ ПРОГНОЗИРОВАНИЯ Кр-ИНДЕКСА

Магнитосфера Земли представляет собой сложную динамическую систему, состояние которой зависит как от внешнего воздействия (прежде всего со стороны СВ), так и от собственной предыстории. По этой причине в число входных признаков при использовании методов машинного обучения в качестве индикаторов внешнего воздействия включают обычно параметры СВ и ММП. Требование алгоритмов машинного обучения к стационарности или хотя бы квазистационарности используемых временных рядов приводит к необходимости использования данных космических аппаратов, находящихся на гало-орбитах вокруг точки Лагранжа L1 между Солнцем и Землей, влиянием изменения положения которых на результаты измерений в первом приближении можно пренебречь. Помимо этого, время подлета солнечного ветра от точки L1 до Земли составляет порядка 40 минут, что позволяет повысить качество прогнозирования при небольших горизонтах прогноза.

Ключевым условием эффективного поступления энергии СВ в магнитосферу Земли является отрицательное значение Bz-компоненты ММП (в системе GSM) [Калегаев В.В. и др.]. Существенное влияние на возможное возникновение геомагнитных возмущений оказывает также значение общей амплитуды ММП. By- и в особенности Bx-компонента ММП практически не оказывают влияния на возмущение магнитосферы Земли. Как правило, авторы настоящей статьи включают значения этих признаков в число входных признаков создаваемых моделей для того, чтобы проверить способность используемых методов отбора и/или оценки существенности входных признаков обнаруживать это.

Ключевыми параметрами СВ, оказывающими влияние на состояние магнитосферы Земли, являются скорость СВ и плотность протонов в СВ. Третьим по значимости среди параметров СВ принято считать его температуру, которая, впрочем, имеет высокую корреляцию со скоростью СВ. По этой причине в данной работе авторами было принято решение среди параметров СВ ограничиться скоростью и плотностью.

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

В то же время, для физических и геофизических исследований используется экваториальный индекс Dst, вычисляемый на основе измерений в приэкваториальной зоне и характеризующий динамику кольцевого тока. Специфика и динамика этого индекса таковы, что использование его значения в качестве входного признака при прогнозировании значения индекса Kp позволяет повысить качество прогноза. По этой причине он используется и в настоящей работе.

Иные признаки за пределами этого базового состава, хотя и используются в некоторых работах, на данный момент не стали общепринятыми. Из их числа авторы обычно используют гармонические функции времени с суточным и годовым периодами, позволяющие методу машинного обучения учесть корреляции изменений геомагнитных индексов с вращением Земли вокруг Солнца и вокруг своей оси (разумеется, подобные корреляции не обязаны иметь причинно-следственный характер).

В работе используются данные о параметрах СВ и ММП с КА ACE, [ACE Browse Hourly Averages 1997-2023], полученные при помощи приборов SWEPAM (Solar Wind Electron Proton Alpha Monitor) и MAG, а также значения геомагнитных индексов с сайта Мирового Центра данных в Киото [World Data Center for Geomagnetism, Kyoto, 1997-2023].

Приведем полный список используемых в данной работе переменных:

1. данные о величине ММП в системе единиц GSM – покомпонентно (Bx, By, Bz) и |B| (модуль ММП), нТл;
2. данные о параметрах плазмы СВ: скорость СВ (v), км/с и плотность протонов (np), см^-3;
3. данные о значении прогнозируемого индекса Кр;
4. данные о значении геомагнитного индекса Dst.
   Все перечисленные переменные используются в виде временных рядов с шагом в 3 часа.
5. Для учета изменений индекса Кр, связанных с вращением Земли вокруг своей оси и вокруг Солнца, на вход алгоритмов машинного обучения подается также информация о часе суток и часе года, представленная в виде значений синуса и косинуса с суточным (hourSin, hourCos) и годовым (daySin, dayCos) периодами.

Для обеспечения возможности учета истории изменения параметров в работе исследуется возможность использования т.н. погружения временных рядов, т.е. помимо текущих значений параметров, алгоритм машинного обучения принимает на вход данные о значениях указанных параметров за последние сутки с шагом в 3 часа.

3. ОБОСНОВАНИЕ ПОСТАНОВКИ ЗАДАЧИ КЛАССИФИКАЦИИ. ОБОСНОВАНИЕ ВЫБОРА ГРАНИЦ КАТЕГОРИЙ

Как правило, при рассмотрении задачи прогнозирования значений того или иного геомагнитного индекса используется регрессионная постановка задачи, в рамках которой прогнозируется конкретное значение индекса в тот или иной момент времени, а оценка качества решения задачи основана на той или иной мере погрешности прогноза, т.е. отличия спрогнозированного значения от реального. Поэтому постановка задачи прогнозирования уровня возмущенности геомагнитного поля, т.е. задачи классификации по определению попадания значения Kp-индекса в тот или иной диапазон, нуждается в дополнительном разъяснении.

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

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

Ясно, что успешность подхода, связанного с рассмотрением ХСС, напрямую связана с качеством решения задачи классификации, которое, в свою очередь, напрямую связано с адекватностью алгоритмов (правил) выделения ХСС. В данной работе мы рассматриваем самый простой вариант, который может не обязательно оказаться оптимальным – выделение ХСС только по уровню индекса Kp. При этом правильность установки границ между категориями является критически важной.

В своих предыдущих работах авторы использовали в качестве тренировочного набора данные, начиная с ноября 1997 года, то есть с начала работы КА АСЕ. Но проведенный анализ показывает, что уровень геомагнитной возмущенности в последнем из завершенных циклов солнечной активности (СА) – 24-м, существенно снизился. Особенно хорошо это заметно на экстремальных магнитных бурях с Кр = 9 – последняя такая буря наблюдалась в 2003 году, то есть почти двадцать лет назад. На рис. 1 представлено распределение Кр-индекса по его величине с 21.10.1997 г. по 07.01.2023 г., а также отдельно для 23, 24 и части 25 цикла СА. Из приведенных гистограмм видно, что число интервалов с Кр = 5 составляет 0.83% для полного массива данных, но если смотреть отдельно 23 и 24 цикл СА, то получается 1.2% для Кр = 5 в 23 цикле и 0.57% в 24 цикле. Для значений Кр, соответствующих невозмущенным периодам, ситуация обратная – например, для Кр = 1 – 10.25% для всего набора данных, 8.97% – для 23 цикла СА, 11.23% – для 24 цикла СА.

Из четвертого столбца рис. 1, где приведены распределения для текущего, 25 цикла СА, можно увидеть, что для упомянутых выше значений Кр доли составляют 0.41% и 1.29% соответственно, и все распределение в целом повторяет распределение 24-го цикла. Поэтому нами было принято решение создать два набора классификационных моделей: один построен на данных 23 и 24 циклов СА, другой – на данных только 24 цикла.

Диаграммы, приведенные на рис. 1, объясняют выбор границ для “спокойного”, “слабо возмущенного” и “возмущенного” периода. Если, как это было принято в 23 цикле, считать буревыми значения Кр, начиная с 5, то в текущем цикле СА мы получим всего 15 бурь. Поэтому были установлены следующие границы “категорий” индекса Kp – “спокойная магнитосфера” – от 0 до 2-, “слабо возмущенная магнитосфера” – от 2 до 3+, и “возмущенная магнитосфера” – 4- и выше.

 

Рис. 1. Гистограммы индекса Kp для разных циклов солнечной активности.

 

4. ОЦЕНКА КАЧЕСТВА КЛАССИФИКАЦИИ

Для построения моделей первого типа в качестве тренировочного набора был использован интервал значений всех переменных с 1997 по 2019 годы.

В качестве тестового набора для оценки качества работы всех моделей был использован интервал значений всех переменных (в том числе и индекса Kp) за 2020‒2022 гг. На рис. 2 изображена динамика индекса Kp для моделей первого типа с разбиением на тренировочный и тестовый наборы.

 

Рис. 2. Недельные максимумы индекса Kp с 1997 года с разбиением на тренировочный и тестовый наборы.

 

Для построения моделей второго типа в качестве тренировочного набора был использован интервал значений всех переменных с 2009 по 2019 годы.

Для оценки качества построенных моделей для каждого класса была выбрана метрика F₁, которая вычисляется по формуле:

F₁ = 2 (Precision × Recall) / (Precision + Recall),

где Precision (точность) – доля примеров, для которых модель правильно выбрала данную категорию, от всех примеров, для которых модель выбрала данную категорию; Recall (полнота) – доля примеров, распознанных моделью как пример данной категории, от всех примеров данной категории.

Для оценки результатов мы используем точность и полноту отдельно для каждого класса (категории). Однако в качестве финальной метрики качества классификации использовалась усредненная по всем классам метрика F₁.

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

5. ОПИСАНИЕ МЕТОДОВ ПРОГНОЗИРОВАНИЯ

В настоящей работе тестировались различные модели классификации для предсказания категории геомагнитного индекса Kp на двух уровнях погружения – 0 (на вход подается вектор признаков в предыдущий момент времени) и 8 (соответствует 24 часам, на вход подается вектор признаков с погружением за все предыдущие 8 отсчетов). Шаг временного ряда во всех случаях составляет 3 часа.

Перечислим методы, использованные для решения задачи классификации категории геомагнитного индекса Kp:

1. Две различных реализации алгоритма градиентного бустинга (GB) [Friedman, J. H., 2001] (CatBoost [Prokhorenkova, L., et al., 2018], LightGBM [Ke G., et al., 2017]), реализации которых могут существенно различаться в силу используемых эвристик, различных приближений, а также большого количества гиперпараметров, для которых в этих реализациях устанавливаются различные значения по умолчанию.
2. Случайный лес (RF) [Breiman, L., 2001].
3. Логистическая регрессия (LR) [Cox, D.R., 1958].
4. Линейная регрессия с регуляризацией: L1 (Lasso), L2 (Ridge) [Hoerl, A. E., Kennard, R. W., 1970] и L1+L2 (Elastic Net).
5. Многослойный персептрон (Perceptron) [Rumelhart, D.E., et al., 1986].
6. Дополнительно для уровня погружения 8 тестируются две рекуррентные нейронные сети – GRU [Cho, K., et al., 2014] и LSTM [Hochreiter S., Schmidhuber J., 1997]. Для GRU и LSTM на вход подавалась последовательность из векторов длиной 8, сформированная скользящим окном.
7. Тривиальная модель (TM) – инерционный прогноз, для которого прогнозируемое значение категории индекса Kp равно последнему известному значению.

В скобках даны кодовые обозначения моделей, которые будут использованы далее на рисунках. Для обучения всех нейронных сетей использовался метод оптимизации Adam [Kingma D.P., Ba J., 2015] и ранняя остановка по валидационному набору данных.

Для преодоления основной проблемы задачи – несбалансированности данных использована методика выравнивания классов-меньшинств за счет генерации синтетических примеров – SMOTE [Chawla N.V., et al., 2002]. Все перечисленные модели (кроме персептрона, GRU и LSTM) протестированы с использованием SMOTE и без для оценки применимости методики в задаче.

Для всех моделей (в конфигурациях со SMOTE и без) также был протестирована специальная методика кросс-валидации для временных рядов в целях выбора оптимальных гиперпараметров каждой модели (TSS). Осуществлялось разбиение тренировочного набора на блоки для кросс-валидации, после чего модель с лучшими параметрами обучалась заново на полном тренировочном наборе. Размер всех блоков при кросс-валидации был выбран одинаковым, для того чтобы значения метрики качества классификации можно было корректно усреднять по всем наборам (так как при изменении количества данных качество классификации в большинстве случаев должно изменяться). Общая схема кросс-валидации изображена на рис. 3.

 

Рис. 3. Общая схема кросс-валидации.

 

Количество разбиений при кросс-валидации K = 5.

Перечисленные модели протестированы со значениями гиперпараметров по умолчанию и после кросс-валидации.

Проверка применимости выбранной методики кросс-валидации и методики SMOTE представляет отдельный интерес, так как ранее не удалось подтвердить эффективность этих методов при ином разбиении индекса Kp на категории с более выраженной несбалансированностью классов [Gadzhiev I. et al, 2023].

Вывод об эффективности работы каждой из моделей строится на основании сравнения качества классификации по усредненной метрике F₁ с другими моделями и с тривиальной моделью.

6. РЕЗУЛЬТАТЫ

Таким образом, для каждого из рассматриваемых методов машинного обучения было построено по 8 моделей: с использованием и без использования SMOTE, с использованием и без использования подбора оптимальных параметров на основе кросс-валидации для временных рядов, с использованием и без использования погружения временного ряда на 8 отсчетов (24 часа). Следует также отметить, что для градиентного бустинга и для рекуррентных нейронных сетей использовались по две различных реализации. Для каждого метода мы приведем только наилучший результат. Отметим также, что все реализации линейной регрессии показали устойчиво более низкие результаты, чем логистическая регрессия, в связи с чем мы их не приводим.

На рис. 4 изображена метрика качества классификации F₁ для всех моделей первого типа (обученных на данных 23 и 24 циклов) на всех горизонтах. Последний столбец (TM) – метрика качества для тривиальной модели. Видно, что наилучшей модели превзойти тривиальную модель удается на всех горизонтах.

 

Рис. 4. Метрика F₁ на разных горизонтах для моделей первого типа.

 

Использование методики SMOTE (выравнивание представительности категорий) и подбор гиперпараметров при помощи кросс-валидации временных рядов дали положительный эффект на всех горизонтах. При этом на первом горизонте (3 часа) лучшие результаты показал градиентный бустинг в реализации LightGBM (F₁ = 0.677), на остальных горизонтах – случайный лес (метрика F₁ монотонно убывает от 0.589 для горизонта в 6 часов, 0.540 для горизонта в 9 часов и 0.500 для горизонта в 12 часов до 0.438 для горизонта в 24 часа).

Отметим, что использование погружения не дало прироста по качеству работы моделей. Среди возможных причин, почему качество могло не улучшиться, можно отметить следующие:

1. Слишком большая размерность входных данных, что особенно существенно для персептронов.
2. Характерное время изменения некоторых переменных (например, Bz) составляет меньше трех часов, что не может быть учтено при шаге временного ряда в 3 часа.

На следующих этапах исследования предполагается проверить последнее предположение путем использования при погружении для всех параметров, кроме самого индекса Kp, временных рядов с шагом в один час.

На рис. 5 представлены результаты, аналогичные результатам рис. 4, однако полученные для моделей второго типа, обученных только на данных 24 цикла СА. Для большинства горизонтов результаты моделей второго типа оказались несколько хуже, чем результаты моделей первого типа: метрика F₁ от 0.672 при горизонте в 3 часа и 0.576 при горизонте в 6 часов до 0.432 для горизонта в 24 часа.

 

Рис. 5. Метрика F₁ на разных горизонтах для моделей второго типа.

 

В то же время, для моделей второго типа справедливы два основных вывода, сделанных для моделей первого типа. Совместное использование методики SMOTE и выбор гиперпараметров моделей с помощью кросс-валидации временных рядов позволяет превзойти результаты тривиальной модели на всех горизонтах. При этом, использование погружения не дает прироста по качеству работы моделей.

Более интересными являются выводы, которые можно сделать при сравнении лучших результатов, полученных для моделей первого и второго типа (рис. 6). Видно, что для всех горизонтов прогнозирования модели второго типа уступили моделям первого типа. Это означает, что при переходе от моделей первого типа к моделям второго типа негативный эффект от двукратного уменьшения размера обучающей выборки оказался более существенным, чем позитивный эффект от приближения распределения данных в обучающей выборке к распределению данных в тестовом наборе и от устранения части данных, которые могут отличаться от данных тестового набора по динамике изменения основных признаков (если предположить, что магнитосфера Земли, представляющая собой сложную динамическую систему, во время разных циклов СА находится в разных характерных состояниях).

 

Рис. 6. Сравнение лучших результатов моделей первого и второго типа.

 

Стоит отметить, что этот результат согласуется с результатами, полученными авторами в их предыдущих работах [Dolenko et al., 2014; Доленко и др., 2016; Myagkova et al., 2016]. Предположение о том, что разбиение полного массива данных на части с построением отдельной модели для каждой из частей может привести к улучшению результатов прогнозирования по сравнению с единой моделью, построенной на полном массиве данных, проверялось несколькими способами: при разбиении массива данных по наличию или отсутствию геомагнитного возмущения, при его разбиении путем кластеризации, а также при его разбиении по уровню СА. Во всех случаях оказывалось, что единая модель, обученная на полном массиве данных с максимально возможным количеством примеров, показывала наилучшие результаты. В этом смысле результат, полученный в данной работе, согласуется с нашими предыдущими исследованиями. Однако исследования в этом направлении следует продолжить, прежде всего в направлении оптимизации методов выделения характерных скрытых состояний.

На рис. 7 представлена точность, а на рис. 8 – полнота модели-победителя первого и второго типа для каждого горизонта в разрезе по категориям. Видно, что при монотонном убывании обеих метрик для каждой из категорий с увеличением горизонта прогнозирования наблюдается также вполне ожидаемый эффект, связанный с разной представительностью категорий. Несмотря на использование методики SMOTE для частичного выравнивания представительности, точность и полнота для категорий 2 и 3, соответствующих возмущенной магнитосфере, оказываются кратно ниже соответствующих показателей для категории 1, представительность которой намного выше.

 

Рис. 7. Точность по категориям для моделей первого типа (слева) и второго типа (справа).

 

Рис. 8. Полнота по категориям для моделей первого типа (слева) и второго типа (справа).

 

При этом низкие значения точности (рис. 7) означают повышенное количество ошибок 1 рода (“ложных тревог”, т.е. ситуаций, когда модель ошибочно отнесла пример к категории, соответствующей возбужденному состоянию магнитосферы). Более критичными являются низкие значения полноты (рис. 8) – повышенное количество ошибок 2 рода (“пропусков цели”, т.е. ситуаций, когда модель не смогла правильно идентифицировать принадлежность примера к категории, соответствующей возбужденному состоянию магнитосферы). Наблюдаемая картина означает, что следует изменить относительные веса точности и полноты в используемой F-мере, т.е. перейти от метрики F₁, например, к метрике F₂ – как для оценки качества модели, так и для отбора и сравнения моделей. Авторы планируют сделать это в своей следующей работе.

Обращает на себя внимание существенное повышение значений полноты для категории 2, соответствующей возмущенной магнитосфере, при переходе от моделей первого типа к моделям второго типа, обученным только на данных 24-го солнечного цикла (за счет некоторого ухудшения этого показателя для категории 1 – невозмущенная магнитосфера, а также за счет некоторого ухудшения точности). Результаты моделей второго типа по показателю полноты на всех горизонтах превышают результаты тривиальной модели, что уже неплохо по той причине, что переходы между категориями случаются сравнительно редко, и поэтому результаты тривиальной модели оказываются достаточно высокими. Все это может свидетельствовать о различиях во взаимосвязях между процессами на Солнце, в солнечном ветре и магнитосфере Земли в “сильном” 23-м цикле СА, с одной стороны, и в “слабых” 24-м и 25-м циклах, с другой. Из этого следует, что необходимо более тщательно подходить к выбору временного диапазона данных, на котором строятся модели.

С точки зрения анализа взаимосвязей между параметрами СВ и ММП, значением Dst-индекса и вращением Земли вокруг своей оси и вокруг Солнца, с одной стороны, и возмущением магнитосферы в терминах Kp- индекса, с другой стороны, представляет интерес сравнение относительных значений значимости (существенности) входных признаков задачи, получаемых из модели-победителя.

На рис. 9 показана относительная значимость входных признаков для лучшей модели с горизонтом прогнозирования 3 часа (градиентный бустинг в реализации LightGBM с использованием SMOTE и подбора параметров на основе кросс-валидации). Распределения значимости, показанные на рис.9 в виде “ящиков с усами”, получены из анализа полного набора деревьев в методе градиентного бустинга. (Средняя линия “ящика” соответствует медианному значению значимости, нижняя и верхняя границы “ящика” – 1-му и 3-му квартилям, краевые засечки нижнего и верхнего “усов” – 5-му и 95-му процентилям.)

 

Рис. 9. Относительная значимость входных признаков для лучшей модели с горизонтом прогнозирования 3 часа (градиентный бустинг в реализации LightGBM с использованием SMOTE и подбора параметров на основе кросс-валидации).

 

Как видно из рис. 9, максимально значимым является сам индекс Кр, что полностью ожидаемо. Далее по значимости идут две фазы (cos и sin) гармонической функции, соответствующей часу суток, что является не настолько очевидным, но отражает влияние на планетарный индекс распределения наземных станций, которые используются для его определения, по земному шару. Далее по значимости следуют скорость СВ, а также амплитуда и Bz-компонента ММП, что также является объяснимым с точки зрения физики – к Земле приходит корональный выброс массы или высокоскоростной поток СВ, возрастают скорость и плотность, а в случае южного направления ММП развивается геомагнитное возмущение. Относительно низкое (шестое по списку) место Bz может быть объяснено тем, что в данном исследовании используется трехчасовое усреднение временных рядов, а южное (отрицательное) Вz, как правило, сохраняется в течение существенно меньшего времени. Низкие величины значимости для компонент ММП Bx и By полностью согласуются с физическими представлениями и подтверждают правильность методики проведения настоящего исследования.

7. ЗАКЛЮЧЕНИЕ

В данной работе были протестированы различные методы классификации при прогнозировании категории геомагнитного индекса Kp по предыдущим значениям самого индекса Kp, параметров СВ и ММП и амплитуды Dst-индекса, а также по значениям времени и даты прогнозирования. Было показано, что методика кросс-валидации временных рядов и методика SMOTE для преодоления несбалансированности классов позволяют превзойти качество предсказания тривиальной инерционной модели (значение прогноза равно предыдущему известному значению категории индекса Kp) на горизонте до 24 часов (с шагом 3 часа), что показывает эффективность применения данных методик в задаче классификации геомагнитного индекса Kp при использовании выбранного разбиения индекса на категории. Модели градиентного бустинга (реализация LightGBM) и случайного леса с использованием этих методик показывают лучшие результаты по сравнению со всеми остальными моделями.

Было показано, что использование погружения (подачи на вход модели предыдущих данных) на 24 часа не дает прироста качества работы моделей. Причиной этому может быть слишком большая размерность входных данных; следовательно, одним из направлений дальнейших исследований должно являться использование отбора признаков. Кроме того, следует протестировать использование часового шага при погружении некоторых переменных, для того чтобы точнее учесть характерное время их изменения.

Ограничение набора входных данных только данными 24 цикла СА, более близкого по статистическим показателям (доли категорий по количеству примеров, условные среднее и медианное значения индекса Kp) к тестовому набору данных (начало 25 цикла СА, 2020-2022 гг.) по сравнению с полным набором и тем более 23 циклом СА, не привело к улучшению результатов прогнозирования в терминах метрики F₁, использованной в данной работе как основная метрика для оценки качества модели. Однако при более детальном рассмотрении результатов в терминах точности и полноты оказалось, что при ограничении набора входных данных 24-м циклом солнечной активности полнота модели существенно увеличивается на всех горизонтах, превышая полноту тривиальной модели. Это может свидетельствовать о том, что различие внешних проявлений СА в 23 и 24 циклах означает изменение некоторых закономерностей и взаимосвязей между процессами на Солнце, в солнечном ветре и магнитосфере Земли.

В целом низкие значения полноты моделей-победителей на всех горизонтах прогнозирования свидетельствуют о необходимости перехода к F-метрике с иными весами точности и полноты по сравнению с метрикой F₁, использованной в настоящей работе.

Результаты анализа относительной значимости входных признаков для лучшей модели с горизонтом прогнозирования 3 часа (градиентный бустинг в реализации LightGBM с использованием SMOTE и подбора параметров на основе кросс-валидации) согласуются с физическими представлениями о существенности параметров СВ и ММП и других переменных для прогнозирования геомагнитных индексов. Была также обнаружена существенная взаимосвязь измеряемого уровня геомагнитного возмущения (измеряемого индекса Kp) с вращением Земли вокруг своей оси. Учитывая то обстоятельство, что при получении оценок значимости входных признаков задачи не использовалась какая-либо априорная информация об их физическом смысле, их совпадение с физическими представлениями свидетельствует об адекватности построенных моделей.

В контексте исследования подхода, связанного с выделением характерных скрытых состояний (ХСС) магнитосферы Земли как динамической системы и раздельным прогнозированием значения геомагнитного индекса в рамках каждого состояния, в ходе дальнейших исследований следует сравнить между собой качество прогнозирования значений индекса Kp методами машинного обучения следующими способами: а) при построении единой прогнозирующей модели без выделения ХСС; б) при построении отдельной прогнозирующей модели в рамках каждой целевой категории, определяемой текущим значением индекса Kp так, как это было сделано в настоящей работе; в) при построении отдельной прогнозирующей модели в рамках каждого состояния, определяемого алгоритмом классификации так, как это было сделано в настоящей работе. Полученные в настоящей работе результаты свидетельствуют о необходимости разработки и использования более изощренных методов выделения ХСС магнитосферы Земли.

ФИНАНСИРОВАНИЕ РАБОТЫ

Исследование выполнено за счет Российского научного фонда, грант № 23-21-00237, (https://rscf.ru/project/23-21-00237/).

Авторлар туралы

I. Gadzhiev

Skobeltsyn Institute of Nuclear Physics, Moscow State University

Хат алмасуға жауапты Автор.
Email: ismailgadzhievff@gmail.com
Ресей, Moscow

О Barinov

Skobeltsyn Institute of Nuclear Physics, Moscow State University

Email: obar@sinp.msu.ru
Ресей, Москва

I. Myagkova

Skobeltsyn Institute of Nuclear Physics, Moscow State University

Email: irina@srd.sinp.msu.ru
Ресей, Москва

S. Dolenko

Skobeltsyn Institute of Nuclear Physics, Moscow State University

Email: dolenko@srd.sinp.msu.ru
Ресей, Москва

Әдебиет тізімі

Қосымша файлдар