A Method for Predicting Geomagnetic Storms Based on Deep Learning Neural Networks Using Time Series of Matrix Observations of the URAGAN Muon Hodoscope

全文:

1. ВВЕДЕНИЕ

1.1. Общие положения, наблюдения мюонного годоскопа и геомагнитные бури

Прогнозирование геомагнитных бурь (ГМБ) является актуальной задачей геофизики, солнечно-земной физики и народного хозяйства. Проблема прогнозирования ГМБ к настоящему времени не решена исчерпывающим образом. Качество прогнозирования ГМБ зависит от вида исходных данных и применяемых математических технологий.

Используемая в рамках данной статьи экспериментальная информация, основывается на измерениях величин потоков мюонов (ПМ). Мюоны – это элементарные частицы, образующиеся в результате ядерных реакций между релятивисткими космическими протонами и отдельными атомами, которые входят в состав верхних слоев земной атмосферы [Мурзин, 2007; Астапов, 2014]. В меняющихся во времени величинах ПМ содержится информация о возможных экстремальных событиях в гелиосфере и потенциально возникающих ГМБ.

Измеренные ПМ формируются в матричные временные ряды наблюдений от мюонного годоскопа (МГ) УРАГАН, сконструированного в НИЯУ МИФИ [Yashin et al, 2015; Barbashina et al., 2008]. МГ-наблюдения берутся с сайта Experimental Complex NEVOD .

МГ является компьютеризованным прибором, измеряющим величины ПМ путем подсчета числа мюонов, попадающих на апертуру МГ для заданной системы телесных углов и установленного интервала дискретности по времени. В результате подсчета числа мюонов формируется временной ряд исходных минутных матриц МГ-наблюдений Y (i, j, T_0k, Tk), T – интервал дискретности по времени, k – временной индекс, обычно, T = 1 / 60 час (1 мин), T_0k – длительность интервалов регистрации, T_0k ≤ T. Интервалы T_0k не являются фиксированными и меняются в диапазоне ≈ секунд; исходные МГ-наблюдения состоят из случайных положительных целых чисел. Азимутальные и зенитные дискретные переменные (i, j), i = 1, ..., N₁,  j = 1, ..., N₂, N₁ = 90, N₂ = 76, определяют телесные углы (φ_i, ϑ_j), φ_i = Δφ (i – 1), ϑ_j = Δϑ (j – 1), Δφ = 1°, Δϑ = 4°, в которых производятся оценки величин ПМ. В ячейки с номерами (i, j) для МГ-матриц помещаются дискретно по времени количества собранных мюонов в виде целых чисел, соответствующих телесным углам (φ_i, ϑ_j).

Первичные минутные нормированные матричные МГ-наблюдения Y (i, j, Tk) формируются в виде временных рядов, определяющих количества зарегистрированных частиц, приведенных к одной секунде

$Y(i,j,Tk)=Y(i, j, T_{0k}, Tk)/T_{0k}, k=1,2,...$. (1)

Вторичные часовые нормированные матричные МГ-наблюдения формируются из минутных наблюдений (1) путем их усреднения на 60-ти минутных интервалах времени

$$\begin{gathered} Y(i,j,T_{0}n)=\sum _{k=k_1\left(n\right)}^{k_2\left(n\right)}Y(i,j,Tk)/60, \\ {k}_1\left(n\right)=1+60(n-1), \\ {k}_2\left(n\right)=60n, n=\mathrm{1,2,...}, T_0=60T. \end{gathered}$$ (2)

Экстремальные события в гелиосфере могут приводить к возникновению возмущений магнитного поля Земли. Геомагнитными бурями (ГМБ) считаются геомагнитные возмущения, имеющие относительную амплитуду колебаний больше некоторой заданной.

Геомагнитные индексы [Menvielle, 2011] служат своего рода индикаторами активности магнитного поля Земли.

Данная статья базируется на применении скалярных значений временных рядов Dst-индекса, которые берутся из сайта World Data Center of Geomagnetism. Dst-индекс измеряется в нанотеслах и формируется почасовым усреднениям значений меридианальных составляющих векторов напряженности геомагнитного поля для магнитных обсерваторий, расположенных по экватору Земли [Sugiura, 1991]. Cпокойным состояниям магнитосферы соответствует Dst-индекс со значениями в пределах (+20 ÷ –50) нТл. Для ГМБ принято считать, что Dst-индекс принимает значения в диапазоне (–50 ÷ –150) нТл. Значения Dst-индекса меньше, чем –150 нT относятся к исключительным случаям.

Достаточно распространенными для задач геофизики являются k-, kр- и aр-индексы с трехчасовым и суточным шагом дискретности. Однако, в нашем случае значения Dst-индекса в большей степени, чем перечисленные выше индексы, подходят для задачи прогнозирования ГМБ, вследствие их часовой дискретизации.

1.2. Обзор публикаций по прогнозированию геомагнитных бурь

Различаются долгосрочные, среднесрочные и краткосрочные прогнозирования ГМБ. В основном, прогнозирование производится в рамках различных технологий анализа космической погоды. Существует целый ряд служб различной ведомственной и государственной принадлежности, которые публикуют информацию по космической погоде, включая прогнозирование ГМБ.

В РФ в ИЗМИРАН активно работает Центр прогнозирования космической погоды [Гайдаш и др., 2016; Gaidash et al., 2017], поставляющий потребителям более 20 видов продуктов, связанных с краткосрочным-долгосрочным прогнозированием космической погоды и ГМБ. На сайте Институт Земного магнетизма, ионосферы и распространения радиоволн им. Н.В. Пушкова РАН предлагаются 3-х суточные прогнозы трехчасовых k- и kp-индексов, 8-ми суточные геомагнитные прогнозы на основе среднесуточных aр-индексов и вероятностей наибольших величин индексов в сутках, прогнозы геомагнитной активности на основе ар-индекса на 55 суток вперед, двухнедельные графики солнечных F10.7 и значенеий ар-среднесуточного индекса вместе с результатами их 24-часового прогнозирования, которые вычисляются на основе специально разработанных регрессионных моделей.

В институте космических исследований РАН в отделе физики космической плазмы реализуется проект по исследованию космической погоды [Плазменная гелиофизика, 2008], в котором разрабатываются, в том числе, методы и алгоритмы для осуществления краткосрочного прогнозирования ГМБ [Петрукович и др., 2016].

Целесообразно отметить работу Центра анализа космической погоды НИИЯФ им. акад. Д.В. Скобельцына МГУ и Службы анализа космической погоды ИПГ им. акад. Е.К. Федорова Росгидромета, занимающихся разработками в рассматриваемой предметной области

Заметной является деятельность Лаборатории рентгеновской астрономии ФИАН в решении проблемы прогнозирования ГМБ с привлечением технологий солнечной радиографии. На сайте Laboratory of X-ray приводятся вычисления kр-индекса на 3-х дневный и 27-дневный интервалы прогнозирования, которые осуществляются с учетом оценок параметров солнечного ветра.

Необходимо упомянуть подразделение Space Whether Prediction Center NOАA (National Оceanic and Аtmospheric Аdministration), CША, занимающееся вариантами прогнозирований ГМБ. На сайте NOAA помещаются данные по kр и aр-индексам для 1÷4-х недельных интервалов прогнозирований. Для решений задач прогнозирований ГМБ используются снимки от размещенных в космосе ультрафиолетового телескопа EIT (Extreme Ultraviolet Imaging Telescope) и солнечного коронографа LASCO (Large Angle and Spectrometric Coronagraph) National Оceanic, с помощью которых оцениваются характеристики нагрева солнечной короны и процессов в CME (Сoronal Мass Еjections), а также определяется влияние короны на солнечный ветер. Изображения LASCO используются для прогнозирующей модели WSA-Enlil National Оceanic, которая действует с 2011 года. Это крупномасштабная физическая гелиосферная модель предназначена для обеспечения заблаговременного предупреждения за 1–4 дня об изменениях в структурах солнечного ветра и направлений движений CME.

Для прямого прогнозирования ГМБ используется космический аппарат АСЕ (Advanced Composition Explorer), USA, National Оceanic, запущенный в 1997г, постоянно находящийся в точке Лагранжа, расположенной на расстоянии 1.44 млн. км от Земли на прямой, соединяющей Землю и Солнце. АСЕ непрерывно регистрирует параметры космического излучения и автоматически передает их Землю; параметры излучения размещаются на сайте NOAA для принятия последующих оперативных решений. На основе АСЕ осуществляется 1-часовое прогнозирование ГМБ.

Для получения необходимых результатов в рассматриваемой предметной области прогнозирования ГМБ применяется целый ряд математических подходов. Так, в публикациях [Белов, 2018; Григорьев, 2019], посвященных прогнозированию ГМБ и космической погоде, описывается метод глобальной съемки на основе данных сайта NMDB: The Neutron Monitor DataBase (НМ). В тезисах симпозиума ИЗМИРАН [Прогноз, 2023], содержится, фактически, обзор современных методов прогнозирования экстремальных гелиосферных событий.

НС широко используются в задачах распознавания-прогнозирования экстремальных событий в гелиосфере и магнитосфере [Бархатов, 2010; Бархатов и Ревунов, 2017]. Дополнения к указанным монографиям содержатся в материалах по НС на сайте нижегородской научно-исследовательской лаборатории физики Солнечно-Земных связей, где, в частности, описывается 12-часовое прогнозирование значений Dst-индекса на основе информации о параметрах околоземной плазмы и межпланетного магнитного поля.

Целый ряд публикаций, связанных с НС, Dst-индексом и ГМБ отличаются вариантами применяемых методов, программных продуктов и используемых источников информации. Указанные обстоятельства вносят значительное разнообразие в постановки задач прогнозирования.

В статьях [Pallochia, 2006; Gruet, 2018; Stepanova, 2000; Ефиторов, 2018; Dolenko, 2005; Широкий, 2015; Кристофоретти, 2022] помещены материалы, касающиеся исследований возможностей прогнозирований ГМБ, реализуемых с помощью различных вариантов НС.

Публикация [Белов и др., 2018] послужившая методической основой для данной статьи, посвящена задачам распознавания и прогнозирования ГМБ с применением МГ-наблюдений и наблюдений нейтронного монитора (НМ). Совместные МГ- и НМ-наблюдения формировались из скалярной переменной от МГ, вычисленной в виде усредненной суммы почасовых МГ-наблюдений, и скалярной переменной от НМ, вычисленной с помощью метода глобальной съемки и введеннных функций изотропных составляющих НM-наблюдений.

Приведенные здесь материалы по физике мюонов, конструкции МГ и сведений, касающихся вариантов задач прогнозирования ГМБ, легли в основу реализуемых в статье подходов. Для решения рассматриваемой задачи предложено использовать наблюдения мюонного годоскопа, ввиду их потенциальной перспективности, и применить НС глубокого обучения, которые являются эффективным инструментом анализа больших данных и в настоящее время широко применяются в научных и технических задачах прогнозирования [Агарвал, 2021; Goodfellow, 2016]. Для реализации прогнозирования здесь предложена математическая технология аппроксимацонной экстраполяции для матричных переменных.

Цель предлагаемой статьи состоит в разработке метода прогнозирования геомагнитных бурь на основе использоваания временных рядов матричных наблюдений МГ УРАГАН, временных рядов скалярных значений Dst-индекса, применения нейронных сетей глубокого обучения и реализации математической технологии аппроксимационного экстраполирования.

2. АНАЛИЗ НАБЛЮДЕНИЙ МЮОННОГО ГОДОСКОПА И Dst-ИНДЕКС

Все переменные, которые использовались для НС в данной статье, были подвергнуты синхронизации и дискретизации с часовым шагом в единой шкале времени UTC. Временной индекс n определял моменты дискретизации T₀n, T₀ = 1 час. Образовывались матричный временной ряд из часовых нормированных МГ-наблюдений Y_M (T₀n), состоящий из матриц Y (i, j, T₀n) (2) размерности (N₂, N₁), и скалярный временной ряд значений Dst-индекса Y_D (T₀n); начальный индекс n = 1, соответствовал моменту времени 01.01.2013, 00.00 часов, 1 ≤ n ≤ n_f, n_f = 78888.

Особенности МГ-наблюдений Y_M (T₀n) были проанализированы. На рис. 1 помещен пример 2D-изображения матрицы часовых нормированных МГ-наблюдений Y_M (T₀n₀), где индекс n₀ = 23424 относился к моменту времени 16.08.2015, 14 час. Из данного рисунка можно увидеть, что матрицы Y_M (T₀n₀) подвержены значительным модуляциям и шумам.

 

Рис. 1. 2D-изображение матрицы часовых нормированных МГ-наблюдений Y_M (T₀n₀).

 

Конструкционный фактор является доминирующим в модуляциях МГ-наблюдений. Элементы матриц для зенитных углов с индексами j = 25 ÷ 30 принимают максимальные значения, для углов с индексами j = 1 ÷ 2, 75 ÷ 76 – минимальные значения, которые различаются в несколько десятков раз. Так, математические ожидания (м.о.) m (j, n₀) и среднеквадратичные значения (с.к.з) σ (j, n₀)

$$\begin{gathered} m(j, n_0)=\frac{1}{N_1}\sum _{i=1}^{N_1}Y(i, j, T_0{n}_0), \\ {\sigma }^2(j,n_0)=\frac{1}{N_1-1}\sum _{i=1}^N{\left(Y(i,j,T_0{n}_0)-m(j,n_0)\right)}^2 \end{gathered}$$

для элементов Y (i, j, T₀n) с индексами j₁ = 3, j₂ = 27, j₃ = 74 принимают значения m (j₁, n₀) = 0.0643, m (j₂, n₀) = 0.3833, m (j₃, n₀) = 0.0082, σ (j₁, n₀) = 0.0242, σ (j₂, n₀) = 0.0565, σ (j₃, n₀) = 0.0073. Отношения вычисленных максимальных и минимальных м.о. и с.к.з. составляют величины ≈ 46.3 и ≈ 67.7 соответственно. Очевидно, эти величины указанных модуляций служат препятствием для прогнозирования ГМБ с учетом шумов и малых значений модуляций из-за гелиосферных экстремальных событий.

На рис. 2 помещены примеры графиков фрагментов исходных переменных, YM₀₁ (n) = YM₀₁ (i₀₁, j₀₁, T₀n), i₀₁ = 12, j₀₁ = 27, YM₀₂ (i₀₂, j₀₂, T₀n), i₀₂ = 22, j₀₂ = 15 и функции Dst-индекса Y_D (n) = Y_D (T₀n) для семимесячного временного участка 01.04.2013–31.10.2013. Отчетливо видно, что средние значения переменных YM₀₁, YM₀₂, для зенитных индексов j₀₁ = 27 и j₀₂ = 15 заметно различаются.

 

Рис. 2. Графики фрагментов исходных переменных YM₀₁, YM₀₂ и функции Dst-индекса Y_D для семимесячного временного участка.

 

Почти общепринято, что временные моменты времени T₀n, в которых реализуются ГМБ, определяются на основе сравнения текущих значений Dst-индекса Y_D (n) с задаваемым порогом Y_D0. Для нTл видно, на данном участке для 2013 г. имели место 13 ГМБ-событий.

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

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

Для пояснения особенностей возможных функциональных связей между МГ-наблюдениями и значениями Dst-индекса на участке с начальной – конечной датой 2013.06.20 – 2013.07.16 с укрупненным временным масштабом по сравнению с рис. 2 для YM₀₁ (n) и Y_D (n) с n₁ ≤ n ≤ n₂ и n₁ = 4101, n₂ = 4678, соответствовавших указанным датам, были произведены расчеты, результаты которых приведены на рис. 3: (а) – график функции Dst-индекса Y_D (n); (б) – график переменной YM₀₁ (n); (в) – график функции усреднения S_M (T₀n).

 

Рис. 3. Результаты расчетов переменных Y_D, YM₀₁ и функции усреднения S_M для анализа функциональных связей между МГ-наблюдениями и функцией Dst-индекса.

 

Из рис. 3а для Y_D (n) видно, что на этом участке произошли, по крайней мере 5 ГМБ, для которых имели место значительные понижения значения Dst-индекса. Однако, рассмотрение функции YM₀₁ (n) из рис. 3б не позволяет сделать заключение о видимой ее связи с функцией Y_D (n) – вполне просматриваемым понижениям значений Dst-индекса не в полной мере соответствовали какие-либо заметные изменения для функции YM₀₁ (n). Заключение о функциональной связи может быть получено на основе анализа функций Y_D (n) и S_M (T₀n), вычисленной на основе для МГ-наблюдений Y_M (i, j, T₀n) – рис. 3в.

$S_{Ì}\left(T_{0}n\right)=\frac{1}{N_1{N}_2}\sum _{i=1}^{N_1}\sum _{j=1}^{N_2}{Y}_M(i, j, T_{0}n), n_{01}\le n\le n_{02}$.

Видно, что понижениям функции Y_D (n) могут быть поставлены соответствующие понижения для функции S_M (T₀n). Так, ГМБ на рис. 3а происходили в моменты времени, которым были сопоставлены индексы n_GS1 = 4198, n_GS2 = 4316, n_GS3 = 4464, n_GS4 = 4550, n_GS5 = 4654. Понижения для S_M (T₀n) на рис. 3в происходили в моменты времени с индексами n_МГ1 = 4178, n_МГ2 = 4298, n_МГ3 = 4430, n_МГ4 = 4534, n_МГ5 = 4630. Нетрудно видеть из визуального рассмотрения рис. 3а и рис. 3в, что индексы n_МГ реализуются с опережениями по времени по отношению к индексам n_GS.

Графики рис. 3 позволили сделать вывод о существовании между значениями Dst-индекса Y_D (n) и МГ-наблюдениями – функциями Y_M (i, j, T₀n) (S_M (T₀n)) очевидной связи, однако, с достаточно сложным функциональным характером.

3. ПРОГНОЗИРОВАНИЕ ГЕОМАГНИТНЫХ БУРЬ НА ОСНОВЕ НЕЙРОННОЙ СЕТИ

3.1. Предварительная подготовка информации для нейронной сети

Предварительная подготовка временных рядов из МГ-наблюдений и Dst-индексов производилась для повышения эффективности работы НС [Кристофоретти, 2022]. Устранение низкочастотных аддитивных трендовых составляющих приближенно решалось путем нахождения м.о. m (Y_M (i, j)), m (Y_D) и вычислений центрированных переменных Y_MC (i, j, T₀n), Y_DC (n) для индексов 1 ≤ n ≤ n_f и смещенного для Dst-порога Y_{D0, C}

$$\begin{gathered} m\left(Y_M\right(i, j\left)\right)=\frac{1}{n_f}\sum _{n=1}^{n_f}{Y}_M(i, j, T_{0}n), i=\mathrm{1,} ...,N_1, \\ j=\mathrm{1,} ...,N_2, m\left(Y_D\right)=\frac{1}{n_f}\sum _{n=1}^{n_f}{Y}_D\left(n\right), \\ {Y}_{MC}(i,j,T_{0}n)=Y_M(i,j,T_{0}n)-m\left(Y_M\right(i,j\left)\right), \\ {Y}_{DC}\left(n\right)=Y_D\left(n\right)-m\left(Y_D\right), Y_{D\mathrm{0,}C}=Y_{D0}-m\left(Y_D\right). \end{gathered}$$

В центрированных переменных Y_MC (i, j, T₀n) устранение высокочастотных составляющих производилось с помощью цифрового низкочастотного фильтра [Filter Design, 2024] c нормированной частотой среза w_c; формировались отфильтрованные переменные, обозначаемые, как Y_{MC, Ф} (i, j, T₀n), 1 ≤ n ≤ n_f, i = 1, ..., N₁,  j = 1, ..., N₂.

На рисунке 4 помещены графики примеров фрагментов отфильтрованнных переменных Y_{M01C, Ф} = Y_{M01C, Ф} (i₁, j₁, T₀n), i₁ =12, j₁ = 5, Y_{M02C, Ф} = Y_{M02C, Ф} (i₂, j₂, T₀n), i₂ =50, j₂ = 25 с выбранной частотой среза w_c = 0.7 для четырехдневного временного интервала 01.04.2013−04.04.2013 г. в зависимости от времени T₀n.

Можно заключить, анализируя рис. 4, что переменные Y_{M01C, Ф}, Y_{M02C, Ф} хорошо приспособлены для экстраполирования: очевидно, что допустима их успешная аппроксимационная экстраполяция на основе параболических функций, по крайней мере, на 5−7 дискретных временных шагов с T₀ = 1 час.

 

Рис. 4. Графики фрагментов отфильтрованных переменных Y_{M01C, Ф}, Y_{M02C, Ф} для четырехдневного интервала.

 

3.2. Схема постановки задачи прогнозирования геомагнитных бурь

С учетом анализа наблюдений мюонного годоскопа, отмеченного в п. 2, вполне возможно осуществление прогнозирования ГМБ с помощью построения модели экстраполированных значений Dst-индекса Y_DE (n), которые сравниваются с задаваемым порогом Y_{D0, C}. На основе анализа неравенств Y_DE (n) ≥ Y_{D0, C} находятся индексы n, для которых выполняются рассматриваемые неравенства. Тем самым, определяются временные интервалы прогнозирования и моменты времени прогнозированных начал и окончаний ГМБ.

Матричные МГ-наблюдения Y_M (n) и скалярные значения Dst-индекса Y_D (n) рассматриваются на девятилетнем временном участке 01.01.2013−31.12.2021, на котором реализуются НС-этапы. На этапе обучения 01.01.2013−02.02.2020, 16 час., 1 ≤ n ≤ n_f0, n_f0 = 62128 используются МГ-наблюдения и значения Dst-индекса; результат этого этапа состоит в построении НС-модели на основе применения НС глубокого обучения. Используются МГ-наблюдения для участка 02.02.2020, 17 час. – 31.12.2020, n_f0 + 1 ≤ n ≤ n_f1, n_f1 = 70128; на этом этапе реализуется контроль этапа обучения – построенной НС-модели. На этапе валидации вычисляются экстраполированные МГ-наблюдения для участка 01.01.2021 – 31.12.2021, n_f1 + 1 ≤ n ≤ n_f, n_f = 78888 и находятся экстраполированные модельные оценки Dst-индекса Y_DE (n) с помощью формул экстраполяции и НС-модели на основе только Y_M (n); Dst-индекс применяется только для определения вероятностей правильных и ложных прогнозирований ГМБ.

Подалгоритм №1 этапа предварительной цифровой обработки для исходного матричного временного ряда МГ-наблюдений Y_M (n) и временного ряда Dst-индекса Y_D (n) осуществляет их фильтрацию с целью устранения низкочастотных составляющих и высокочастотных шумов. Результаты предварительной обработки обозначаются, как Y_{M1C, Ф}, Y_{M2C, Ф}, Y_D1C и Y_D2C; индекс 1 относится к этапу обучения, индекс 2 к этапу валидации построения модельных оценок Dst-индекса. Порог после предварительной обработки обозначается, как Y_{D0, C}.

Подалгоритм №2 реализует НС-обучение и НС-валидацию. На его вход подаются отьфильтрованный матричный временной ряд Y_{M1, C, Ф}, обучающий скалярный временной ряд Y_D1C, 1 ≤ n ≤ n_f0 и формируется НС-модель. Также в этом подалгоритме реализуется валидационный контроль качества НС-модели, n_f0 + 1 ≤ n ≤ n_f1.

Подалгоритмы №3 и №4 производили экстраполяцию МГ-наблюдений и вычисление экстраполированных модельных оценок Dst-индекса Y_{DCФ, E}; для этой цели использовались переменные Y_{M2C, Ф} на интервале n_f1 + 1 ≤ n ≤ n_f и НС-модель от этапа обучения.

Подалгоритм №5 реализует процедуру принятия решений для прогнозирования ГМБ путем сравнения вычисленных экстраполированных модельных оценок Dst-индекса Y_{DCФ, E} с порогом Y_{D0, C}. В этом подалгоритме используется Y_D2C для оценивания вероятностных характеристик прогнозирования.

3.3. Аппроксимационная экстраполяция матричных переменных мюонного годоскопа и прогнозирование геомагнитных бурь

Аппроксимационная экстраполяции на n_E шагов реализуется для переменных Y_{M2C, Ф} (i, j, n) относительно индексов n, удовлетворяющих неравенствам n_f1 + Δn − n_E + 1 ≤ n ≤ n_f − n_E и для всех сочетаний индексов i = 1, ..., N₁,  j = 1, ..., N₂. Производится построение системы аппроксимационных параболических моделей для текущего временного момента с индексом n с шагом экстраполяция на n_E шагов для индексов n + 1, n + 2, ..., n + n_E. Задается параметр n_d − количество индексов на интервале n, n − 1, ..., n − n_d + 1, на котором производятся построения аппроксимационных параболических моделей c₀ + c₁s + c₂s² с коэффициентами c^T = (c₀, c₁, c₂) и аппроксимационными функционалами для n_d точек, n − n_d + 1 ≤ s ≤ n, i = 1, ..., N₁,  j = 1, ..., N₂

$F\left(c,Y_{M2C,Ф}(i, j), n, n_d\right)=\sum _{s=n-n_d+1}^n{\left(Y_{M2C,Ф}(i, j, s)-c_0-c_{1}s-c_2{s}^2\right)}^2$.

Находятся векторы оптимальных коэффициентов аппроксимационных моделей c° (i, j, n, n_d) путем минимизации введенных квадратичных функционалов для всех сочетаний индексов i, j и заданных n, n_d

$c°(i, j, n, n_d)=\arg \left\{\underset{c}{\min }F\right(c,Y_{M2C,Ф}(i, j), n, n_d\left)\right\}$.

Образуются векторы с экстраполяцией Y_{M2C, ФE} (i, j, n, n_E, n_d) размерности Δn, i = 1, ..., N₁,  j = 1, ..., N₂, определенные в индексах n, n − 1, ..., n − n_E + 1, n − n_E, ..., n − Δn + 1.

Вектор Y_{M2C, ФE} (i, j, n, n_E, n_d) состоит из n_E координат, вычисляемых по формулам экстраполяции, и Δn − n_E координат, которые формируются основе координат от исходного вектора. Координаты промежуточного вектора ${\overline{Y}}_{M2C,ФE}(i, j, n, n_E, n_d)$ представляются следующими формулами

$$\begin{gathered} {\overline{Y}}_{M2C,ФE}(i, j, n+\mathrm{1,} n_E, n_d)= \\ =c_0^{\circ }(i, j, n, n_d)+c_1^{\circ }(i, j, n, n_d)\left(1\right)+c_2^{\circ }(i, j, n, n_d){\left(1\right)}^2, \\ {\overline{Y}}_{M2C,ФE}(i, j, n+2, n_E, n_d)= \\ =c_0^{\circ }(i, j, n, n_d)+c_1^{\circ }(i, j, n, n_d)\left(2\right)+c_2^{\circ }(i, j, n, n_d){\left(2\right)}^2, \\ {\overline{Y}}_{M2C,ФE}(i, j, n+2, n_E, n_d)= \\ =c_0^{\circ }(i, j, n, n_d)+c_1^{\circ }(i, j, n, n_d)\left(n_E\right)+c_2^{\circ }(i, j, n, n_d){\left(n_E\right)}^2, {\overline{Y}}_{M2C,ФE}(i, j, n, n_E, n_d)=Y_{M2C,Ф}(i, j, n, n_E), \\ {\overline{Y}}_{M2C,ФE}(i, j, n-\mathrm{1,} n_E, n_d)=Y_{M2C,Ф}(i, j, n-\mathrm{1,} n_E), \\ {\overline{Y}}_{M2C,ФE}(i, j, n-\Delta n+n_E-\mathrm{1,} n_E, n_d)=Y_{M2C,Ф}(i, j, n-\Delta n+n_E-\mathrm{1,} n_E). \end{gathered}$$

От промежуточных векторов ${\overline{Y}}_{M2C,ФE}(i, j, n+\mathrm{1,} n_E, n_d)$ осуществляется переход к экстраполированным векторам Y_{M2C, ФE} (i, j, n, n_E, n_d)

$$\begin{gathered} Y_{M2C,ФE}(i, j, n+\mathrm{1,} n_E, n_d)={\overline{Y}}_{M2C,ФE}(i, j, n+\mathrm{1,} n_E, n_d), \\ {Y}_{M2C,ФE}(i, j, n-n_E+\mathrm{2,} n_E, n_d)={\overline{Y}}_{M2C,ФE}(i, j, n+\mathrm{2,} n_E, n_d), \\ {Y}_{M2C,ФE}(i, j, n-n_E+1, n_E, n_d)={\overline{Y}}_{M2C,ФE}(i, j, n+1, n_E, n_d), \\ {Y}_{M2C,ФE}(i, j, n-n_E, n_E, n_d)={\overline{Y}}_{M2C,ФE}(i, j, n, n_E, n_d), \\ {Y}_{M2C,ФE}(i, j, n-n_E-1, n_E, n_d)={\overline{Y}}_{M2C,ФE}(i, j, n-1, n_E, n_d), \\ {Y}_{M2C,ФE}(i, j, n-∆n+1, n_E, n_d)={\overline{Y}}_{M2C,ФE}(i, j, ∆n+n_E-1, n_E, n_d), \end{gathered}$$

на основе которых собираются пакеты экстраполированных входных матриц Y_{M2C, ФE}. Далее, пакеты Y_{M2C, ФE} в скользящем режиме с единичным шагом поступают на вход системы НС-построения экстраполированных центрированнных модельных оценок Dst-индекса Y_DC (n, n_E, n_d) со смещением на n_E индексов. Оценки экстраполированных центрированных скалярных модельных оценок Dst-индекса без смещений Y_{DC, E} (n, n_E, n_d) представляются следующим соотношением

$Y_{DC,E}(n, n_d)=Y_{DC}(n-n_E, n_E, n_d)$. (3)

4. ОПЕРАЦИИ НЕЙРОННОЙ СЕТИ RESNET34 ДЛЯ ПРОГНОЗИРОВАНИЯ ГЕОМАГНИТНЫХ БУРЬ

На рис. 5 представляется укрупненная блок-схема вычислительных операций реализации НС-прогнозирования ГМБ.

 

Рис. 5. Блок-схема вычислительных операций нейросетевого прогнозирования ГМБ.

 

В блоке 1 производится предварительная обработка информации. На его вход подаются временные ряды пакетов Dst-векторов Y_D и матриц МГ-наблюдений Y_M. На выходе формируются: временные ряды пакетов центрированных Dst-векторов Y_D1C, центрированных отфильтрованных пакетов матриц МГ-наблюдений Y_{M1C, Ф} и временные ряды пакетов центрированных и отфильтрованных пакетов матриц МГ-наблюдений Y_{M2C, Ф}. В блоке 2 на основе Y_D1C, Y_{M1C, Ф} осуществляется процедуры обучения и валидации, в результате которых формируется НС-модель. В блоке 3 для переменных Y_{M2C, Ф} реализуется их экстраполяция на n_E шагов, выход этого блока обозначается, как Y_{M2C, ФE}. В блоке 4 на основе использования НС-модели и экстраполированных переменных Y_{M2C, ФE} вычисляются скалярные оценки экстраполированных Dst-индексов Y_{DCФ, E}. В блоке 5 с применением сравнений экстраполированных переменных Y_{DCФ, E} и порога Y_{D0, C} производится принятие решений о прогнозировании ГМБ; также в блок 5 поступает Y_D2C для оценивания вероятностей правильного и ложного прогнозирования.

Для решения сформулированной задачи прогнозирования ГМБ, связанной с цифровой обработкой большого объема данных, были использованы облачные технологии. Выполнение ресурсоемких вычислительных работ, которые требовались для данной НС глубокого обучения, осуществлялось с использованием вычислительных ресурсов вычислительного центра НИЯУ МИФИ.

Было произведено формирование варианта структуры НС глубокого обучения с учетом [Wintoft P.&Wik M, 2021]. Использовался язык программирования Python и библиотека модулей PyTorch. Проводились вычислительные эксперименты для рассмотрения работы модулей различных нейросетевых архитектур, таких как EfficientNet, VGG, DenseNet, Inception-v3 и ResNеt34 из библиотеки PyTorch. Для решения предложенной в статье задачи была взята модификация ResNet34, которая обеспечивала наилучшую точность и достаточно малое количество тренируемых параметров по сравнению с другими модулями. Дополнительный аргумент в пользу выбора модуля ResNet34 состоял в том, что при исследовании процедур обучения указанных НС, этот модуль занял второе место, работая без аугментации, т.е. без дополнения исходных реальных данных нереальными значениями.

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

Поскольку МГ-наблюдения занимали существенно большой объем памяти и модель Dst-индексов зависима от значительного числа параметров, то для решения поставленной задачи прогнозирования потребовались масштабные нейросетевые архитектуры, такие, как ResNet34. По мнению авторов, структуры НС из нескольких сверточных и полносвязных слоев обладали малыми возможностями и не в полной мере подходили для рассматриваемой задачи.

Анализ полученных результатов предварительных вычислительных экспериментов с НС ResNet34 и располагаемыми МГ-наблюдениями позволил назначить Δn размерность пакета входных матриц.

5. ПРИНЯТИЕ РЕШЕНИЯ ПО ПРОГНОЗИРОВАНИЮ ГЕОМАГНИТНЫХ БУРЬ, ВЫЧИСЛЕНИЕ ВЕРОЯТНОСТНЫХ ОЦЕНОК И ОЦЕНИВАНИЕ ИНТЕРВАЛОВ ПРОГНОЗИРОВАНИЯ

Правило принятия решений по прогнозированию ГМБ основывалось на процедуре бинарной классификации – сравнении вычисленных прогнозированных оценок Dst-индекса Y_DC (n, n_E, n_d) по формуле (3) из п. 3 с центрированным порогом Y_{D0, C}. Реализация сравнения производилась путем рассмотрения неравенств со значениями Dst-индекса и порогами для индексов n временного участка валидации

$$\begin{gathered} Y_{DC}(n,n_E,n_d)\le Y_{D\mathrm{0,}C}, \\ {n}_{f1}+\Delta n-n_E+1\le n\le n_f-n_E+1. \end{gathered}$$ (4)

С помощью неравенств (4) подсчитывались число правильных N_{R, ГМБ} (n_E, n_d, Y_{D0, C}) и ложных прогнозирований ГМБ N_{F, ГМБ} (n_E, n_d, Y_{D0, C}), поскольку на рассматриваемом участке валидации для каждого значения индекса  заранее известно случилось или не случилось для него событие ГМБ. Для N_{R, ГМБ} (n_E, n_d, Y_{D0, C}), N_{F, ГМБ} (n_E, n_d, Y_{D0, C}) записывались формулы, для которых использовалась функция Хевисайда: H (x) = 1, x ≥ 0, H (x) = 0, x < 0,

$$\begin{gathered} N_{R,ГМБ}(n_E, n_d, Y_{D\mathrm{0,}C})=\sum _{n=n_{f1}+\Delta n-n_E+1}^{n_f-n_E+1}H\left(Y_{D\mathrm{0,}C}-Y_{DC}\left(n\right)\right)H\left(Y_{D\mathrm{0,}C}-Y_{DC,E}(n,n_d)\right), \\ {N}_{R,ГМБ}(n_E, n_d, Y_{D\mathrm{0,}C})=\sum _{n=n_{f1}+\Delta n-n_E+1}^{n_f-n_E+1}H\left(Y_{DC}\left(n\right)-Y_{D\mathrm{0,}C}\right)H\left(Y_{D\mathrm{0,}C}-Y_{DC,E}(n,n_d)\right). \end{gathered}$$ (5)

Вычисление оценок вероятностей правильных и ложных прогнозирований ГМБ производилось путем рассмотрения неравенств для центрированных значений Dst-индекса на участке валидации

$Y_{DC}\left(n\right)\le Y_{D0,C}, n_{1f}+\Delta n-n_e+1\le n\le n_f-n_e+1$. (6)

С помощью неравенств (6) подсчитывалось число индексов с ГМБ N_{1, ГМБ} (Y_{D0, C}) и число индексов без ГМБ N_{0, ГМБ} (Y_{D0, C}),

$$\begin{gathered} N_{\mathrm{1,}ГМБ}\left(Y_{D\mathrm{0,}C}\right)=\sum _{n=n_{f1}+\Delta n-n_E+1}^{n_f-n_E+1}H(Y_{D\mathrm{0,}C}-Y_{D,C}(n)), \\ {N}_{\mathrm{0,}ГМБ}(Y_{D\mathrm{0,}C})=n_f-n{}_{f1}-\Delta n-N_{\mathrm{1,}ГМБ}. \end{gathered}$$ (7)

Вероятности β (n_E, n_d, Y_{D0, C}) правильных и α (n_E, n_d, Y_{D0, C}) ложных прогнозирований представлялись следующими формулами

$$\begin{gathered} \beta (n_E, n_d, Y_{D\mathrm{0,}C})=N_{R,ГМБ}(n_E, n_d, Y_{D\mathrm{0,}C})/N_{\mathrm{1,}ГМБ}\left(Y_{D\mathrm{0,}C}\right), \\ \alpha (n_E, n_d, Y_{D\mathrm{0,}C})=N_{F,ГМБ}(n_E, n_d, Y_{D\mathrm{0,}C})/N_{\mathrm{0,}ГМБ}\left(Y_{D\mathrm{0,}C}\right). \end{gathered}$$ (8)

Выражения (8) в статье использовались, как критерии качества прогнозирования ГМБ для бинарной классификации. Разумеется, вполне возможны другие критерии (метрики), такие, как матрица ошибок (Confusion matrix), меткость (Accuracy), точность (Precision), полнота (Recall), специфичность (Specificity), F1-мера (F1-score) и ее развитие метрика P4 [Loginom, 2023]. Критерии (формула 8) соответствовали известным англоязычным аналогам Fall-out и Recall и, в нашем случае, являлись более предпочтительными с физической точки зрения, чем выше перечисленные критерии.

Оценивание интервалов прогнозирования для ГМБ c номером r производилось на основе неравенств (6) и (4). Вычислялись временные границы ГМБ с применением формулы (6)

$n_1(r, Y_{D\mathrm{0,}C})\le n\le n_2(r, Y_{D\mathrm{0,}C}), r=\mathrm{1,} \mathrm{2,} ..., r_0$, (9)

где r₀ – число реальных ГМБ на участке валидации. Вычисление временных границ интервалов для прогнозированных ГМБ находилось с использованием выражения (4)

$n_1(s, n_E, n_d, Y_{D\mathrm{0,}C})\le n\le n_2(s, n_E, n_d, Y_{D\mathrm{0,}C}), s=\mathrm{1,} \mathrm{2,} ..., s_0$, (10)

где s₀ – число прогнозированных ГМБ, в общем случае имело место s₀ ≠ r₀.

Анализ взаимного расположения границ интервалов (n₁ (r, Y_{D0, C}), n₂ (r, Y_{D0, C})) и (n₁ (s, n_E, n_d, Y_{D0, C}), n₂ (s, n_E, n_d, Y_{D0, C})), для (9–10) позволил решить вопрос о прогнозировании ГМБ. Рассматривались для некоторого значения индекса r разности с индексами s = 1, 2, ..., s₀

$dn(r, s, n_E, n_d, Y_{D\mathrm{0,}C})=n_1(r, Y_{D\mathrm{0,}C})-n_1(s, n_E, n_d, Y_{D\mathrm{0,}C})$.

Полагалось, что к множеству S_p принадлежат индексы s, которые обеспечивали неравенство dn (r, s, n_E, n_d, Y_{D0, C}) ≥ 0. Вычислялось значение индекса s°, при котором достигался максимум положительной разности

$s^{\circ }=\arg \left\{\underset{s\in S_p}{\max }\right(dn(r,s,n_E,n_d,Y_{D\mathrm{0,}C})\}$.

Очевидно, что ГМБ с номером r для s° может считаться прогнозированной для индексов , которые удовлетворяли неравенству

$n_1(s^{\circ },n_E,n_d,Y_{D\mathrm{0,}C})\le n\le n_1(r,Y_{D\mathrm{0,}C})$. (11)

6. ЭКСПЕРИМЕНТАЛЬНОЕ ИССЛЕДОВАНИЕ ХАРАКТЕРИСТИК МЕТОДА ПРОГНОЗИРОВАНИЯ ГЕОМАГНИТНЫХ БУРЬ

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

Рассматривались МГ-наблюдения Y_{M2C, ФE} (n) и значения Dst-индекса Y_D2C (n) на шестимесячном участке 01.07.2021−31.01.2022 (75200 ≤ n ≤ 78888). На этом участке принимались параметры n_d = 5, n = 2, 4, 6 и формирровалась последовательность значений порогов Y_D0 (l) = Y_D0N + ΔY_D0 (l – 1), l = 1, ..., l_f, l_f = 9, Y_D0N = –70 нT, ΔY_D0 = 5 нTл. В данном случае, Y_D0 было заменено на Y_{D0, C} c учетом того, что эти параметры отличались на константу, Y_{D0, C} (l) = Y_D0 (l) = 5.78 нT. Вычислялись экстраполированные модельные оценки Dst-индексов Y_{DCФ, E}, (n, n_d, Y_{D0, C} (l)), которые сравнивались с порогом Y_{D0, C}.

Определялись по (4−8) оценки вероятностей β° (n_E, n_d, Y_D0) и α° (n_E, n_d, Y_D0). Для семейства функций с дискретными n_Em, m = 1, 2, 3, n_E1 = 2, n_E2 = 4, n_E3 = 6, чтобы обеспечить удобства визуализации, вводились смещения ΔY_D0 (m – 1) для оценок вероятностей

$$\begin{gathered} \alpha 0°(n_{Em}, n_d, Y_{D0})=\alpha °(n_{Em}, n_d, Y_{D0}-\Delta Y_{D0}(m-1\left)\right), -70nT\le Y_{D0}\le -30нT-\Delta Y_{D0}(m-1), \\ \beta 0°(n_{Em}, n_d, Y_{D0})=\alpha °(n_{Em}, n_d, Y_{D0}+\Delta Y_{D0}(m-1\left)\right), -70nT+\Delta Y_{D0}(m-1)\le Y_{D0}\le -30нT. \end{gathered}$$ (12)

На рис. 6 представлены семейства графиков, вычисленные на основе (12), вероятностных характеристик прогнозированния ГМБ: (а) – оценки вероятностей ложных прогнозирований α0° (n_E, n_d, Y_D0); (б) – оценки вероятностей правильных прогнозирования β0° (n_E, n_d, Y_D0). Вычисления оценок вероятностей произведены в зависимости от n_E, n_d = 5 и Y_D0, которые давали возможность оценить качество прогнозирования ГМБ. Так, при Y_D0 = –45нТл и n_E = 2, 4, 6 для диапазона значений вероятностей ложных прогнозирований α0° (n_E, n_d, Y_D0) = 0.001 ÷ 0.004, имел место диапазон вероятностей правильных прогнозирований β0° (n_E, n_d, Y_D0) = 0.85 ÷ 0.9; при Y_D0 = –35нТл и тех же выбранных n_E для диапазона значений вероятностей ложных прогнозирований α0° (n_E, n_d, Y_D0) = 0.027 ÷ 0.03 имел место диапазон вероятностей правильных прогнозирований β0° (n_E, n_d, Y_D0) = 0.98 ÷ 0.99.

 

Рис. 6. Графики оценок вероятностей ложного и правильного прогнозирования ГМБ.

 

Рассмотрение рис. 6а и рис. 6б с вычисленными оценками вероятностных характеристик, позволило заключить, что предложенный метод обеспечивал удовлетворительные для геофизической практики значения вероятностей правильных и ложных прогнозирований ГМБ.

Исходя из анализа вариантов задач прогнозирования, обычно, “удовлетворительными” для геофизической практики для обозначенных выше параметрах, связанных с рис. 6а, рис. 6б, могут быть приняты величины ограничений по вероятности ложного прогнозирования для α0° (n_E, n_d, Y_{D0, C}), в среднем, от ≈ 0.0025 до ≈ 0.028 и диапазону для вероятности правильного прогнозирования β0° (n_E, n_d, Y_{D0, C}), в среднем, от ≈ 0.85 до ≈ 0.99.

Производилось сравнение результатов методов прогнозирования ГМБ на основе матричных и скалярных временных рядов наблюдений от МГ с использованием НС. Очевидно, что матричные МГ-наблюдения явились более информативными по сравнению со скалярными МГ-наблюдениями. Так, из материалов настоящей статьи следовало, что для предлагаемого прогнозирования ГМБ с НС глубокого обучения при наилучших сочетаниях параметров для матричных МГ-наблюдений Y_D0 = –45 нT, n_d = 5, n_E = 6, вероятности ложных прогнозирований не превосходили значений 0.001 ÷ 0.004 и вероятности правильных прогнозирований составляли величины 0.85 ÷ 0.9. На основе материалов статьи [Гетманов, 2022], где рассмотрено прогнозирование ГМБ для скалярных МГ-наблюдений и сверточных НС, при наилучших сочетаниях параметров Y_D0 = –50 нT, n_d = 5, n_E = 5, вероятности ложных и правильных прогнозирований, соответственно, принимают значения 0.06 и 0.3.

Такое значительное отличие в вероятностных характеристиках прогнозирования для матричных и скалярных МГ-наблюдений связанно с тем, что в [Гетманов, 2022] рассматривается задача прогнозирования ГМБ при совместном использовании скалярных наблюдений от МГ и нейтронного монитора. Скалярные МГ-наблюдения, как показано в [Гетманов, 2022], явились малоэффективными для прогнозирования ГМБ. Матричные МГ-наблюдения, как было продемонстрировано в настоящей статье, позволяют обеспечить достаточно хорошие характеристики прогнозирования ГМБ.

6.2. Результаты оценивания интервалов прогнозирования геомагнитных бурь

Были рассмотрены Dst-индексы Y_D2C (n) и модельные оценки экстраполированных значений Dst-индекса Y_{DCФ, E} (n) на годовом временном участке 01.01.2021−31.12.2021 г.; в данном подпараграфе принято начало отсчета − 01.01.2021 г., 00:00 с начальным значением индекса n = 1. Вычислялись даты и значения индексов временных границ интервалов с точностью до часа для реальных ГМБ n₁ (r, Y_{D0, C}), n₂ (r, Y_{D0, C}), r = 1, 2, ..., r₀ (9) и временные границы для прогнозированных ГМБ n₁ (s, n_E, n_d, Y_{D0, C}), n₂ (s, n_E, n_d, Y_{D0, C}), s = 1, 2, ..., s₀ (10). В качестве регулирующих параметров принимались величины n_d, n_E. Назначалось исходное значение порога распознавания Y_D0 = –50 нTл; в результате центрирования было получено значение порога Y_{D0, C} = –44.22 нTл.

На рис. 7 представлен график функции Dst-переменной Y_D2C на рассматриваемом участке, пунктирной линией обозначено значение порога Y_{D0, C}.

 

Рис. 7. График Dst-переменной Y_D2C.

 

Видно по значениям Dst-индекса, что на указанном временном участке произошли 3 ГМБ. В табл. 1 приведены результаты вычислений дат начал-окончаний ГМБ и временных индексов границ n₁ (r, Y_{D0, C}), n₂ (r, Y_{D0, C}), r = 1, 2, ..., r₀, соответствующих графику (рис. 7), где r₀ = 3.

 

Таблица 1. Вычисления дат начал и окончаний реальных ГМБ и временных индексов границ

	Даты начал ГМБ	Даты окончаний ГМБ
1	27.08.2021 10:00	27.08.2021, 20:00	5722	5732
2	17.09.2021, 13:00	17.09.2021, 15:00	6229	6231
3	03.11.2021, 19:00	04.11.2021, 17:00	7363	7385

 

На рис. 8 для фиксированного знаачения n_d и параметров экстраполяции n_E помещены примеры графиков функций оценок экстраполированных переменных Y_{DCФ, E}: (а) – для n_d = 5, n_E = 2; (б) – для n_d = 5, n_E = 6; (в) – для n_d = 5, n_E = 10. На основе оценок переменнных Y_{DCФ, E} реализовывались нейросетевые прогнозирования; пунктирные линии определяли порог Y_{D0, C}.

 

Рис. 8. Графики экстраполированных переменных Y_{DCФ, E} для параметров экстраполяции n_E.

 

Приведенные здесь графики иллюстрировали важное обстоятельство, которое явилось отчетливо видимым: реализация правильных прогнозирований (Correct Prediction – CP) сопровождалась погрешностями – пропусками прогнозирований (Missed Prediction – MP) и ложными прогнозированями (False Prediction – FP). В табл. 2 представлены результаты расчетов прогнозирований в зависимости от значений n_E = 2, 4, 6, 8, 10 и n_d = 5 для дат начал и окончаний прогнозированных ГМБ и соответсвующих временных участков с границами n₁ (s, n_E, n_d, Y_{D0, C}), n₂ (s, n_E, n_d, Y_{D0, C}) в часах, s = 1, 2, ..., s₀.

 

Таблица 2. Вычисления дат начал и окончаний прогнозированных ГМБ и временных участков с границами n₁ (s, n_E, n_d, Y_{D0, C}), n₂ (s, n_E, n_d, Y_{D0, C}) в часах

	Параметр	Даты начал прогнозированных ГМБ	Даты окончаний прогнозированных ГМБ	n1 (s, nE, nd, YD0, C)	n2 (s, nE, nd, YD0, C)
1	2	27.08.2021, 02:00	27.08.2021, 11:00	5720	5729
2		03.11.2021, 11:00	04.11.2021, 05:00	7361	7369
1	4	27.08.2021, 05:00	27.08.2021, 16:00	5717	5728
2		03.11.2021, 15:00	04.11.2021, 08:00	7359	7376
1	6	27.08.2021, 03:00	27.08.2021, 15:00	5715	5727
2		03.11.2021, 13:00	04.11.2021, 05:00	7357	7373
1	8	26.08.2021, 22:00	27.08.2021, 14:00	5710	5726
2		03.11.2021, 11:00	04.11.2021, 03:00	7355	7371
3		05.11.2021, 22:00	05.11.2021, 08:00	7390	7392
1	10	14.08.2021, 23:00	15.08.2021, 01:00	5423	5425
2		26.08.2021, 22:00	26.08.2021, 14:00	5709	5726
3		17.09.2021, 05:00	17.09.2021, 07:00	6221	6223
4		17.09.2021, 09:00	17.09.2021, 13:00	6225	6229
5		03.11.2021, 10:00	04.11.2021, 02:00	7354	7370
6		04.11.2021, 21:00	05.11.2021, 01:00	7389	7393

 

Расчеты, помещенные в табл. 2, показали, что для n_E = 2, 6 реализовались 2 Correct Prediction – n_CP = 2 для ГМБ №1, 3 и 1 Missed Prediction – n_MP = 1 для ГМБ №2. Вероятно, ГМБ с №2 с малой амплитудой и длительностью 2 часа, воспринималась при нейросетевых расчетах, как шумовое возмущение, реализовавшееся, как пропуск прогнозирования. Расчеты для  показали реализацию n_CP = 2, n_MP = 1 и 4 False Preduction – n_FP = 4. Должно быть вполне очевидно, что качество прогнозирования определялось параметром экстраполяции n_E = 8, 10, в особенности для больших значений , когда, по вполне естественным причинам, нарушались условия успешных экстраполяций, что приводило к увеличению числа ложных прогнозирований.

Графическое сопоставление взаимного расположения временных интервалов с реальными и прогнозированными ГМБ, по данным, взятым из табл. 1, табл. 2 и определенным на основе соответствующих значений Dst-индекса может быть произведено с использованием рис. 9: (а) − соответствует ГМБ для даты 27.08.2021; (б) − соответствует ГМБ для даты 17.09.2021 г.

 

Рис. 9. Диаграммы сопоставления интервалов реальных и прогнозированных ГМБ для дат 27.08.2021, 17.09.2021.

 

На рис. 9а и 9б помещены диаграммы, позволившие графически определить взаимное расположение временных участков с реальными и прогнозированными ГМБ, по данным, взятым из табл. 1 и табл. 2, и определенных на основе соответствующих значений Dst-индекса.

Были рассмотрены ГМБ №1, 3 (табл. 1), для их анализа назначены 40-часовые временные интервалы 5705−5735, 7355−7385, на которых жирными линиями отмечены временные участки реальных ГМБ. Обычными линиями с индесами 1, 2, 3, 4, 5 отмечены временные участки прогнозированных ГМБ для n_E = 2, 4, 6, 8, 10.

Рассматривались участки прогнозирований ГМБ (11) на основе вычислений из табл. 2 разностей dn_p (r, s, n_E, n_d, Y_{D0, C}) = n₁ (r, Y_{D0, C}) − n₂ (s, n_E, n_d, Y_{D0, C}). Формировались оценки $d{n}_{p1}(n_E, n_d)=d{n}_{p1}(\mathrm{1,} s_1^{\circ }, n_E, n_d, Y_{D\mathrm{0,}C}), r=1$, $d{n}_{p3}(n_E, n_d)=d{n}_{p2}(\mathrm{1,} s_2^{\circ }, n_E, n_d, Y_{D\mathrm{0,}C}), r=3$. В табл. 3 для n_E = 1, 2, ..., 10 и n_d = 5 приведены результаты расчетов dn_p1 (n_E, n_d), dn_p3 (n_E, n_d) – оценок интервалов прогнозирований ГМБ и числа ложных прогнозирований n_FP (n_E, n_d) для n_E = 1 ÷ 10.

 

Таблица 3. Оценки dn_p1 (n_E, n_d), dn_p3 (n_E, n_d) – интервалов прогнозирований и n_FP (n_E, n_d) – числа ложных прогнозирований для ГМБ

nE	1	2	3	4	5	6	7	8	9	10
dnp1 (nE, nd)	1	2	3	5	6	7	8	12	12	13
dnp3 (nE, nd)	1	2	3	4	5	6	7	8	9	9
nFP (nE, nd)	0	0	0	0	0	0	0	1	2	4

 

Из табл. 3 следовало, что для шагов экстаполяции с n_E = 1 ÷ 3 прогнозирование ГМБ для временных дат 27.08.2021 г, 03.11.2021 г реализовывались с удовлетворительной вероятностью правильного прогнозирования ГМБ – оценки dn₁ (n_E, n_d), dn₃ (n_E, n_d) стабильно возрастали; при этом имел место один пропуск прогнозирования ГМБ для даты 17.08.2021 г. Для шагов экстрополяции n_E = 4 ÷ 7 вероятности правильного прогнозирования ГМБ уменьшались и параметры dn₁ (n_E, n_d), dn₃ (n_E, n_d) начинали незначительно флуктуировать. Для шагов экстраполяции n_E = 8 ÷ 10 вероятности правильного прогнозирования ГМБ далее уменьшались, реализовывалась значительная флуктуация оценок dn₁ (n_E, n_d), dn₃ (n_E, n_d) и одновременно появлялись множество ложных прогнозирований (см. рис. 8в).

Полученная картина прогнозирования ГМБ вполне правдоподобна с физической точки зрения − эффективность прогнозирования естественным образом ухудшалась при увеличении шагов экстраполяции.

Результаты произведенных расчетов, аналогичных помещенным в табл. 2, для n_d = 3, 5, 7 (индексы 1,2,3), связанные с построением оценок интервалов прогнозирований ГМБ dn_p1 (n_E, n_d) и числа ложных прогнозирований ГМБ n_FP (n_E, n_d) в зависимости от параметров n_E, n_d, представлены на рис. 10 и рис. 11. Была рассмотрена ГМБ для r = 1, произошедшая 27.08.2021. Для удобств анализа графиков были применены смещения по вертикальной оси на рис. 10, рис. 11, которые для индексов 2 и 3 составили величины 1 и 2.

 

Рис. 10. Оценки интервалов прогнозирований ГМБ dn_p1 (n_E, n_d).

 

Рис. 11. Оценки числа ложных прогнозирований ГМБ n_FP (n_E, n_d).

 

Можно заключить из представленных рисунков, что прогнозирование ГМБ по предлагаемому методу удовлетворительно осуществлялось для Y_D0 = −45 нT, n_E = 5 ÷ 6, n_d = 5 и, в среднем, принимало значение 4−5 часов до начала ГМБ на Земле, фиксируемой по значениям Dst-индекса.

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

1. Разработаннный статье метод прогнозирования ГМБ на основе матричных временных рядов наблюдений мюонного годоскопа УРАГАН с использованием НС глубокого обучения и алгоритма аппроксимационной экстраполяции оказался работоспособным и удовлетворительно эффективным, что свидельствует о том, что цель статьи достигнута.

2. В результате осуществления процедуры выбора НС-модуля, размерности пакета входных матриц МГ-наблюдений и частоты среза низкочастотного фильтра было установлено, что наиболее подходящим для рассматриваемой задачи прогнозирования ГМБ оказался программный модуль ResNet34 с оптимальной размерностью пакета входных матриц Δn = 60 и частотой среза w_c = 0.7.

2.1. Предложенный метод обеспечивал удовлетворительные для геофизической практики значения вероятностей правильных и ложных прогнозирований ГМБ. Так, при Y_D0 = −45 нTл, Y_D0 = −35 нTл, n_E = 2, 4, 6, n_d = 5 для значений вероятностей ложных прогнозирований α° (n_E, n_d, Y_{D0, C}) находившихся в диапазонах 0.001 ÷ 0.004, 0.027 ÷ 0.03, значения верооятностей правильных прогнозирований β° (n_E, n_d, Y_{D0, C}) находились в диапазонах 0.85 ÷ 0.9 и 0.98 ÷ 0.99.

2.2. Анализ расчетов, связанных с вычислением характеристик предложенного метода, позволил сделать вывод, что при параметрах Y_D0 = −45 нTл, n_E = 5 ÷ 6, n_d = 5 прогнозирование ГМБ осуществлялось с удовлетворительными вероятностными характеристиками и его время, в среднем, составляло 4–5 часов до начала ГМБ на Земле.

3. Произведенное экспериментальное исследование позволило установить, что предложенный в данной статье подход к прогнозированию ГМБ, существенно превысил по эффективности подход, описанный в публикации [Гетманов В., 2022, Геомагнетизм и Аэрономия], и который базировался на формировании из матричных временнных рядов МГ-наблюдений скалярных временных рядов, применении сверточных НС и апроксимационной экстраполяции.
4. Рассмотрение предложенного метода прогнозирования ГМБ на основе матричных временных рядов МГ-наблюдений с использованием НС глубокого обучения и аппроксимационной экстраполяции дало возможность сделать вывод о его перспективности для решения проблем солнечно-земной физики и задач народного хозяйства.

БЛАГОДАРНОСТИ

Авторы статьи выражают благодарность руководству Научно-образовательному Центру “НЕВОД” НИЯУ МИФИ за предоставленные для расчетов наблюдения мюонного годоскопа УРАГАН.

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

Работа выполнена в рамках государственного задания Геофизического центра РАН, утвержденного Министерством науки и высшего образования Российской Федерации.

作者简介

V. Getmanov

Geophysical Center of the Russian Academy of Sciences; Institute of Earth Physics of the Russian Academy of Sciences

编辑信件的主要联系方式.
Email: vgetm2015@yandex.ru
俄罗斯联邦, Moscow; Moscow

A. Gvishiani

Geophysical Center of the Russian Academy of Sciences; Institute of Earth Physics of the Russian Academy of Sciences

Email: a.gvishiani@gcras.ru
俄罗斯联邦, Moscow; Moscow

A. Solovyev

Geophysical Center of the Russian Academy of Sciences; Institute of Earth Physics of the Russian Academy of Sciences

Email: a.soloviev@gcras.ru
俄罗斯联邦, Moscow; Moscow

K. Zaytsev

MEPhI National Research Nuclear University

Email: kszajtsev@mephi.ru
俄罗斯联邦, Moscow

M. Dunaev

MEPhI National Research Nuclear University

Email: max.dunaev@mail.ru
俄罗斯联邦, Moscow

E. Yekhlakov

MEPhI National Research Nuclear University

Email: fruha1980@gmail.com
俄罗斯联邦, Moscow

参考

补充文件