Геомагнетизм и аэрономия, 2023, T. 63, № 1, стр. 80-93

1. ВВЕДЕНИЕ

Все параметры ионосферы испытывают флуктуации относительно некоторого фона, за который обычно принимают спокойную медиану. Знание свойств ионосферы для спокойных условий важно при выделении и анализе возмущений любой природы, поскольку параметры ионосферы для спокойных условий служат фоном, относительно которого определяются эти возмущения. Флуктуации ионосферной плазмы обусловлены как постоянно существующими факторами (например, приливы, планетарные волны и внутренние гравитационные волны в атмосфере), так и эффектами, связанными с магнитосферными суббурями и бурями [Деминов, 2015]. Значительные по амплитуде флуктуации параметров ионосферы называют ионосферными возмущениями.

Основными механизмами ионосферных возмущений являются изменения состава нейтральной атмосферы, термосферные ветры и электродинамические процессы. Оценка и прогноз этих явлений в реальном времени все еще нуждаются в улучшении из-за ограниченности средств мониторинга и точных методов прогноза. Поэтому наиболее надежными средствами для их описания являются модели нейтральной атмосферы и ионосферы. Они находят широкое применение в космических навигационных и радиолокационных системах, которые требуют ионосферной коррекции для учета задержки сигнала или ошибок дальности. Широко известны 3-мерные эмпирические модели ионосферы, позволяющие выполнять долгосрочные (на месяцы и годы) прогнозы ключевых ионосферных параметров [Bilitza et al., 2017; Gulyaeva and Bilitza, 2012 и ссылки там]. Для более детальных краткосрочных ионосферных прогнозов (на часы и дни) применяются адаптивные модели, использующие ассимиляцию наблюдательных данных в реальном времени [Кринберг и др., 1986; Gulyaeva et al., 2013; Galkin et al., 2022]. В литературе представлены краткосрочные прогнозы карт полного электронного содержания GIM-TEC с применением нейронных сетей [Cesaroni et al., 2020; Liu et al., 2020]. Статистические методы успешно применяются для прогноза глобальных карт GIM-ТЕС с заблаговременностью от 1 до 10 дней [Lean, 2019; Iluore and Lu, 2022]. Однако для применения этих методов на практике требуется специальное вычислительное обеспечение, что не всегда возможно для заинтересованных пользователей.

Альтернативная возможность – адаптация медианной модели ионосферы по аэрономическим и динамическим параметрам, связанным с движением высокоширотных структур, выраженным аналитически, как это было принято в динамической модели критической частоты слоя F2 [Шубин и Деминов, 2019]. Преимущество такой коррекции медианной модели по сравнению с ассимиляцией текущих наблюдательных данных заключается в независимости от наличия этих данных в реальном времени.

Ошибки распространения радиосигнала в ионосфере пропорциональны общему содержанию электронов в ионосфере TEC на пути радиолуча. Поэтому удобная для пользователя модель TEC, охватывающая глобальный масштаб и все уровни солнечной и геомагнитной активности, может быть полезной в различных приложениях. В данной работе предложена аэрономическая и динамическая коррекция глобальной медианной модели полного электронного содержания GTEC (Shubin and Gulyaeva, 2022), подобно динамической модели критической частоты foF2 [Шубин и Деминов, 2019]. В разделе 2 представлены подробности аэрономической и динамической коррекции медианной модели GTEC. В разделе 3 приведены примеры применения адаптированной модели GDMTEC. В разделе 4 дано обсуждение результатов и в разделе 5 приведены выводы.

2. АДАПТАЦИЯ МЕДИАННОЙ МОДЕЛИ GTEC К ВОЗМУЩЕННЫМ УСЛОВИЯМ

Для построения адаптированной модели GDMTEC нами применен аэрономический подход к описанию ионосферной бури. Он основан на учете физических процессов, происходящих в термосфере. Впервые попытка использовать этот подход для моделирования foF2 была предпринята в работе [Шубин и Аннакулиев, 1995]. В дальнейшем с появлением более совершенных термосферных моделей типа NRLMSISE-00 [Picone et al., 2002] и более надежных лабораторных измерений констант скоростей основных ионно-молекулярных реакций [Hierl et al., 1997] данный подход получил свое развитие [Шубин и Аннакулиев, 1997; Аннакулиев и др., 2005]. Результаты его применения для прогноза foF2 изложены в работе [Шубин и Деминов, 2019]. Использование эмпирической модели термосферы NRLMSISE-00 [Picone et al., 2002] в качестве составной части полуэмпирической модели бури в ионосфере позволяет отразить важное свойство ионосферы, а именно, зависимость TEC(t) от предыстории изменения геомагнитной активности через соответствующую зависимость параметров термосферы.

В медианной модели GТЕС [Shubin and Gulyaeva, 2022] для того, чтобы отразить предысторию изменения TEC(t) мы ввели новый индекс солнечной активности в виде соотношения:

где интегральный индекс F10.7(τ = 0.96) представляет собой ряд суточных значений F10.7 с экспоненциальным сглаживанием, предшествующих заданному дню [Шубин и Деминов, 2019]. Индекс F10.7 – это поток солнечного радиоизлучения на длине волны 10.7 см в заданный день. Индексы F10.7 измеряются в единицах 10^–22 Вт/(Гц м²). Из τ = exp(–3/T) следует, что характерное время T индекса FF составляет около 27 дней, что соответствует одному обороту Солнца.

Модель GТЕС для геомагнитно-спокойных условий с Ap < 15 нТл, где Ap – суточный индекс геомагнитной активности, была построена на географической сетке с ячейками 5° по долготе и 2° по широте. Зависимость GТЕС от географических координат в фиксированный момент мирового времени UT для двух уровней солнечной активности – низкой L (F10.7 ≤ 80) и высокой H (F10.7 ≥ 120) представляется в виде:

где $\vartheta $ = 90° – φ; φ –широта; λ – долгота; $g_{n}^{m}$ и $h_{n}^{m}$ – коэффициенты разложения, а $P_{n}^{m}\left( {\cos \vartheta } \right)$ – присоединенные функции Лежандра. Для интерполяции по времени было использовано разложение Фурье с периодом T = 24 ч. Для пространственного разложения было использовано N = 15 по широте и M = 10 по долготе сферических гармоник. Количество коэффициентов в разложениях (2) определяется соотношением: J = M(2N – M + 1) + + N + 1. Таким образом, для каждого уровня солнечной активности, месяца и часа UT было получено J = 226 коэффициентов. Модель основана на анализе сферических гармоник глобальных карт полного электронного содержания ионосферы GIM-TEC, построенных в Лаборатории реактивного движения Калифорнийского технологического института, Пасадена, США (https://sideshow.jpl.nasa.gov/pub/iono_daily/).

Метод построения адаптированной модели GDMTEC аналогичен построению модели GDMF2 [Шубин и Деминов, 2019]. В ней модель отрицательной фазы ионосферной бури для критической частоты foF2 выражена в зависимости от спокойного фонового уровня с помощью коэффициента коррекции C_storm. Соответствующее соотношение для отрицательной фазы ионосферной бури полного электронного содержания имеет вид:

где C_storm = R(t)/R_q(t), с функциями из (4) представляет собой отношение параметров термосферы на высоте 300 км для геомагнитно-возмущенных R(t) и спокойных R_q(t) условий, а GТЕС – глобальная медианная модель (2) [Shubin and Gulyaeva, 2022]: Здесь n(O), n(N₂) и n(O₂) – концентрации нейтральных составляющих верхней атмосферы; m(O), m(N₂) и $m\left( {{{{\text{O}}}_{2}}} \right)$ – молекулярные веса нейтральных составляющих; β – коэффициент линейной рекомбинации ионов O⁺; k₁ и k₂ – коэффициенты ионно-молекулярных реакций O⁺ с N₂ и O₂ [Hierl et al., 1997]. Так как TEC(t) ~ ( foF2)²(t), показатель степени в выражении для R(t) уравнения (4) равен 1.30, в отличие от показателя 0.65, полученного для foF2 [Аннакулиев и др., 2005].

Для определения R(t) используется эмпирическая модель NRLMSISE-00 [Picone et al., 2002], которая является глобальной моделью термосферы, применимой при любом уровне солнечной и геомагнитной активности. В ней косвенно учтено, что основной причиной изменений параметров термосферы в периоды геомагнитных бурь является Джоулев нагрев термосферы высоких широт.

В субавроральной и авроральной ионосфере мы вводим динамические поправки к медианному распределению GТЕС. Тогда адаптированная модель GTEC будет иметь следующий вид:

Поправка к GТЕС в области главного ионосферного провала (ГИП) обозначена C_MIT, а в авроральном овале – C_AVR. Поправка C_AVR вычисляется по известным исправленным геомагнитным широтам аврорального пика – Ф_AVR и минимума ГИП – Ф_MIT. Модели положений минимума провала Ф_MIT, аврорального максимума Ф_AVR, а также аэрономические поправки C_MIT и C_AVR в уравнении (5) представлены в виде аналитических зависимостей от местного геомагнитного времени, геомагнитной активности и географической долготы. Выражения для ΔФ₁ имеют вид: (более подробно см. в статьях [Деминов и Шубин, 2018 и Шубин и Деминов, 2019]). Инвариантные геомагнитные широты в уравнениях аэрономических поправок находятся из модели IGRF-13 [Alken et al., 2021]. В качестве индикатора геомагнитной активности нами был введен средневзвешенный по времени индекс K *(t, τ) = 2.1ln(0.2a_p(t, τ) + 1.0), где коэффициент τ = 0.6 определен из условия минимума отклонения моделей от экспериментальных данных, a_p – 3-часовой индекс геомагнитной активности. Из τ = exp(–3/T) следует, что характерное время T для индекса K *(t, τ) составляет около 6 ч.

Величина C_storm также зависит от геомагнитной активности, но с характерным временем T ~ 11 ч, что следует из модели NRLMSISE-00 в уравнениях (3, 4). Таким образом, косвенно учитывается, что реакция аврорального пика ионизации на изменение геомагнитной активности происходит с двумя характерными временами. Вначале изменяется положение и величина этого пика, связанная с изменениями потока высыпающихся авроральных электронов и области локализации этих высыпаний, затем происходит дополнительное изменение величины аврорального пика ионизации, связанное с изменениями температуры и состава термосферы [Деминов и Шубин, 2019; Деминов и др., 2021]. Поправки C_AVR и C_MIT существенны только в области аврорального максимума ионизации и ГИП, соответственно они применяются в этих областях наряду с корректирующим выражением C_stormGТЕС в уравнении (5). Вне этих областей ${{С}_{{{\text{AVR}}}}}{\text{ехр}}\left( { - {\kern 1pt} {{{[{{({\text{Ф}} - {{{\text{Ф}}}_{{{\text{AVR}}}}})} \mathord{\left/ {\vphantom {{({\text{Ф}} - {{{\text{Ф}}}_{{{\text{AVR}}}}})} {\Delta {{{\text{Ф}}}_{{\text{1}}}}}}} \right. \kern-0em} {\Delta {{{\text{Ф}}}_{{\text{1}}}}}}]}}^{{\text{2}}}}} \right) \ll 1$ и C_MIT = 1, поэтому TEC(t) = C_storm GТЕС (уравнение 3) на широтах, расположенных к экватору от ГИП [Деминов и Шубин, 2018].

3. ПРИМЕРЫ ПРИМЕНЕНИЯ АДАПТИРОВАННОЙ МОДЕЛИ GDMTEC

Рассмотрим несколько примеров применения адаптированной модели GDMTEC в периоды геомагнитных возмущений. Для качественной и количественной оценки корректности GDMTEC нами были выбраны три геомагнитные бури в разные сезоны 2021 г: 19.02–24.02.2021 г. (Kp = 4.3 и Dst = –32 нТл), 27.08–01.09.2021 г. (Kp = 5.0 и Dst = –82 нТл) и 17.10–22.10.2021 г. (Kp = 4.0 и Dst = –55 нТл), где в скобках приведены максимальное значение Kp-индекса и минимальное значение Dst для каждой бури при низкой солнечной активности F10.7 от 71.3 до 86.0. Из этих данных следует, что анализируемые бури можно отнести к слабым и умеренным геомагнитным бурям.

Более детально параметры геомагнитной активности и ионосферы для геомагнитной бури 19.02–24.02.2021 г. представлены на рис. 1–3. На верхней панели рис. 1 представлены изменения геомагнитных индексов Kp и Dst с местным временем. На нижней панели точками показаны суточные вариации VTEC (вертикального полного содержания электронов) на трех станциях глобальной навигационной спутниковой системы (ГНСС). Нами были выбраны станции в разных долготных и широтных зонах. Станция Обнинск располагается на средних широтах. Станция Норильск относится к высокоширотной зоне и, как будет показано ниже во время выбранных для тестирования ионосферных бурь она находилась в авроральном овале. Станция Магадан лежит в субавроральной зоне вблизи минимума ГИП. Разработанная нами модель GDMTEC представлена сплошными кривыми, модель IRI-Plas [Gulyaeva and Bilitza, 2012], используемая без адаптации к текущим измерениям ТЕС, отмечена штриховыми линиями. Международная справочная модель ионосферы IRI, расширенная до плазмосферы – IRI-Plas, является одной из наиболее известных эмпирических моделей ионосферы и плазмосферы. Интегрирование электронной плотности в модели IRI-Plas производится от 60 км до 20200 км, которая примерно соответствует высоте орбиты ГНСС. Следует отметить, что при вычислениях VTEC в медианной модели IRI-Plas нами были использованы коэффициенты “URSI” с опцией “STORM”. На осях абсцисс отложено местное время LT, на осях ординат величины VTEC в единицах 1 TECU = 10¹⁶ эл/см². На рис. 1 вертикальными линиями выделены местные времена, для которых на рис. 2 для фиксированного времени 14 UT показаны широтные изменения VTEC на долготах, соответствующих трем выбранным наземным станциям ГНСС. На рис. 2 сплошные линии соответствуют модели GDMTEC, штриховые линии – модели IRI-Plas, точки – значения VTEC, полученные из GIM-TEC, крестики – значения VTEC на соответствующих станциях ГНСС. Значения VTEC на этом и последующих рисунках приведены с точностью ±1° по долготе и ±0.3h по UT. Отметим, что использованные для тестирования модели GDMTEC данные VTEC из GIM за 2021 г. не входили в построенную нами ранее медианную модель GTEC [Shubin and Gulyaeva, 2022]. Для наглядности на рис. 3 показаны положения анализируемых станций относительно границ аврорального овала ионизации и минимума ГИП в географических координатах в то же самое фиксированное мировое время, что на рис. 2. Эти границы являются элементами модели GDMTEC. Сплошной кривой обозначено положение минимума ГИП, штриховой кривой – положение экваториальной границы, точками – положения полярной границы аврорального овала, штрихпунктирной кривой – положение максимума ионизации в авроральном овале. Крестиками обозначены положения станций ГНСС, символом ⊕ обозначен Северный геомагнитный полюс. Видно, что согласно модели GDMTEC ст. Норильск в 14 UT оказывается в авроральной области, точнее, вблизи экваториальной границы этой области. Станция Магадан расположена в субавроральной области, ст. Обнинск – на средних широтах.

Аналогичные параметры геомагнитной активности, ионосферы, а также модели GDMTEC и IRI-Plas для бури 27.08–01.09.2021 г. представлены на рис. 4–6, а для бури 17.10–22.10.2021 г. – на рис. 7–9. Из данных на рис. 6 и 9 можно видеть, что ст. Норильск находилась почти в центре авроральной области, ст. Магадан – в области минимума ГИП для фиксированных моментов мирового времени в эти геомагнитные бури.

Для оценки эффективности модели GDMTEC были проведены ее сопоставления с величинами VTEC, полученными с помощью наземных приемников ГНСС и из GIM‑TEC, а также с моделью IRI-Plas. В качестве количественных оценок результатов сопоставления моделей с наблюдавшимися величинами VTEC используются среднеквадратическая ошибка RMSE (TECU) и среднее относительное отклонение MRD (%) модельных значений VTEC_mod от наблюдаемых данных VTEC_obs:

где N – число часовых значений VTEC наземных приемников ГНСС и GIM‑TEC, которые были использованы для оценок в каждой геомагнитной буре. Результаты представлены в табл. 1 и 2.

В таблице 1 приведены результаты расчета по модели IRI-Plas с двумя разными файлами IG_RZ.dat (с прогнозом до 11.2023 г.), созданными 09.2021 г. (IRI-Plas 1) и 02.2022 г. (IRI-Plas 2). Как видно из этой таблицы, ошибки с IG_RZ.dat от 09.2021 г. (IRI-Plas 1) существенно больше ошибок с IG_RZ.dat от 02.2022 г. (IRI-Plas 2). Файл IG_RZ.dat со скользящими средними за 12 месяцев индексами IG12 и R12 используется в моделях IRI-Plas [Gulyaeva and Bilitza, 2012] и IRI-2016 [Bilitza et al., 2017]. Наибольшие ошибки в моделях IRI-Plas 1 и IRI-Plas 2 достигаются для геомагнитной бури 17.10–22.10.2021 г., так как для этого периода индексы IG12 и R12 являются прогнозными. Результаты количественной оценки, как видно из табл. 1, свидетельствуют о том, что производительность модели GDMTEC, за исключением событий выделенных жирным шрифтом, лучше модели IRI-Plas для всех сезонов.

Количественная оценка точности широтного поведения моделей GDMTEC и IRI-Plas приведена в табл. 2.

В качестве примера на рис. 10 приведено глобальное распределения VTEC для 00 UT 20.02.2021 г., рассчитанное по GDMTEC. Пример показан для возмущенных условий (a_p = 15 нТл, Kp = 3.0). Белыми линиями показано асимметричное положение главного ионосферного провала в Северном и Южном полушариях.

4. ОБСУЖДЕНИЕ РЕЗУЛЬТАТОВ

Ранее в работе [Шубин и Деминов, 2019] было показано применение аэрономической и динамической коррекции foF2 во время геомагнитных бурь на средних и высоких широтах. В данном исследовании такая коррекция была осуществлена к VTEC. Применение адаптированной модели GDMTEC продемонстрировано на трех геомагнитных бурях в условиях низкой солнечной активности, в разные сезоны. Валидация модели GDMTEC была проведена на 3 наземных станциях ГНСС в разных долготных и широтных зонах: среднеширотной ст. Обнинск (55.1° N, 36.6° E), авроральной ст. Норильск (69.4° N, 88.4° E) и субавроральной ст. Магадан (59.6° N, 150.8° E).

Выше отмечалось, что для входных геофизических индексов IG12 и R12 в моделях IRI-Plas и IRI-2016 используется файл IG_RZ.dat. Этот файл обновляется примерно 2 раза в год. Поэтому эти модели применяются в основном для научных исследований или для долговременного прогноза параметров ионосферы в качестве климатических моделей. Преимущество введенного нами индекса FF (уравнение 1) состоит в том, что в нем учитывается предыстории изменения F10.7. Из этого следует, что построенная нами адаптированная модель GDMTEC может использоваться для краткосрочного прогноза VTEC.

В каждой геомагнитной буре, выбранной для проверки модели GDMTEC, для событий 14 UT 19.02.2022 г., 15 UT 27.08.2022 г. и 20 UT 17.10.2022 г. нами были проведены сравнения GDMTEC вдоль долгот соответствующих станциям ГНСС – Обнинск, Норильск, Магадан. Соответствующие широтные изменения модельных и экспериментальных VTEC показаны на рис. 2, 5 и 8. На рис. 3, 6 и 9 расположение этих станций относительно высокоширотных границ представлено в полярной системе координат. На этих рисунках нами была использована модель IRI-Plas с файлом IG_ RZ.dat от 02.2022 г.

Как видно из рис. 2 и 3 ст. Обнинск находилась на средних широтах, ст. Магадан – на экваториальной границе ГИП, а ст. Норильск – на экваториальной границе аврорального овала. При этом модель GDMTEC адекватно описывает широтные изменения как значений VTEC, полученных из GIM-TEC, так и VTEC на станциях ГНСС. Заметные расхождения наблюдаются в Южном полушарии и в области экваториальной аномалии. Модель IRI-Plas, за исключением 16 LT вблизи экватора, практически везде недооценивает VTEC.

Как изображено на рис. 5 и 6 в 15 UT 27.08.2022 г. ст. Магадан за счет смещения высокоширотных границ оказалась практически в минимуме ГИП, а ст. Норильск – в максимуме аврорального овала. Из рис. 5 следует, что модель GDMTEC более адекватно отражает поведение VTEC на долготах ст. Обнинск и ст. Магадан на средних широтах и переоценивает VTEC в области экватора. Модель IRI-Plas по-прежнему почти везде, за исключением области экватора на долготах Обнинска и Норильска, недооценивает VTEC.

Во время бури 17.10–22.10.2021 г. ст. Норильск остается вблизи максимума аврорального овала, ст. Магадан – на экваториальной границе ГИП и ст. Обнинск – в среднеширотной ионосфере. Модель GDMTEC достаточно точно оценивает величины VTEC на всех наземных станциях ГНСС, положение которых в пространстве и времени удовлетворяли условиям Δλ = ±1° и ΔUT = = 0.3 ч.

Как видно из табл. 2, ошибки GDMTEC значительно на 20–60% меньше, чем ошибки IRI-Plas, за исключением широтного хода VTEC в 21 LT 27.08.2022 г. на долготе Норильска (рис. 5). Соответствующие значения RMSE и MRD в табл. 2 выделены жирным шрифтом. Основные расхождения наблюдаются на субавроральных и авроральных широтах, где ионосфера наиболее изменчива.

Как следует из проведенной проверки модели GDMTEC, аэрономическая поправка позволяет адекватно моделировать отрицательные возмущения VTEC во время геомагнитных возмущений. Вместе с этим используемые в ней изменяющиеся границы аврорального овала и положения минимума ГИП во время геомагнитных возмущений достаточно точно описывают изменения VTEC на разных широтах и долготах, а также в разные сезоны. Следует отметить, что точность GDMTEC определяется также точностью используемой медианной модели GTEC в качестве фона для геомагнитно-спокойных условий и модели термосферы NRLMSISE-00 во время геомагнитных бурь.

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

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

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

Использование в качестве входных параметров модели GDMTEC интегральных индексов солнечной F10.7 и геомагнитной a_p активности, которые учитывают предысторию их изменения, позволят осуществлять краткосрочный прогноз TEC. При этом продолжительность и качество прогноза TEC будет определяться точностью прогнозов солнечных F10.7- и геомагнитных a_p-индексов.