Developing space-time model of the main geomagnetic field based on satellite data for the epoch 2015–2020

Abstract

Based on measurements of SWARM satellites, a model of the main geomagnetic field for the 2020 epoch and a model of annual field changes until 2025 were built. These models were included in the international model of the main geomagnetic field IGRF-13. The method of data selection and model calculation is described.

Full Text

1. Введение

Интерес к изучению магнитного поля Земли (МПЗ) и построению его моделей начался с внедрением магнитного компаса в практику мореплавания. Основными методами в изучении МПЗ являются непосредственные наблюдения над пространственным распределением магнитного поля и его вариациями на поверхности Земли и в околоземном пространстве. Уже в самом начале организации магнитных наблюдений особое внимание было обращено на проведение наблюдений в акваториях морей и океанов. В течение первых двух столетий (1500—1700 гг.) после открытия Колумбом склонения — отклонения магнитной стрелки от направления на север — измерения на морях и океанах носили случайный характер. Первая систематическая съемка склонения и наклонения была организована в конце XVII в. английским адмиралтейством. В конце XVIII и в течение всего XIX в. кроме склонения на судах стали проводиться и наблюдения за наклонением и горизонтальной составляющей магнитного поля. Важным моментом в истории развития науки о земном магнетизме после обнаружения пространственного изменения МПЗ явились систематические измерения магнитного склонения. Анализ таких измерений позволил Гильберту [Гильберт, 1956] предложить для описания МПЗ модель однородно намагниченного шара.

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

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

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

2. Аналитическое представление главного геомагнитного поля

Основными источниками данных для моделирования главного магнитного поля в последние десятилетия стали специальные низкоорбитальные спутники для магнитной съемки на протяженном временном интервале: MAGSAT, CHAMP. Проект SWARM [Friis-Christensen et al., 2006], запущенный в конце 2013 г. и включающий три специальных спутника для точных измерений магнитного поля, предоставляет ученым единый набор данных для всего модельного периода 2015—2020 гг., что позволяет создавать модели, основанные только на данных SWARM. Впервые данные SWARM были использованы для построения моделей МПЗ через 2 и 4 года после их запуска [Olsen et al., 2016; Terence et al., 2018].

Общий подход к созданию глобальных магнитных моделей разработан уже давно [Chapman and Bartels, 1940; Alken et al., 2021a] и заключается в минимизации разности экспериментально измеренных значений компонент магнитного поля и теоретических значений, рассчитанных по коэффициентам разложения магнитного поля по сферическим гармоникам. Согласно теории Гаусса потенциал магнитного поля, создаваемый источниками, расположенными внутри сферы радиуса r, можно описать общей формулой:

$\begin{array}{c}U\left(r,\theta ,\lambda \right)=\\ =a\sum _{n=1}^{\infty }\sum _{m=0}^n\left\{\left[g_n^m\cos m\lambda +h_n^m\sin m\lambda \right]\cdot {\left(\frac{a}{r}\right)}^{n+1}\right\}\times \\ \times P_n^m\left(\cos \theta \right)+\\ +\left\{\left[q_n^m\cos m\lambda +s_n^m\sin m\lambda \right]\cdot {\left(\frac{r}{a}\right)}^n\right\}\times \\ \times P_n^m\left(\cos \theta \right),\end{array}$ (1)

где a — средний радиус Земли;

n — степень сферических гармоник;

r — геоцентрическое расстояние;

m — порядок сферических гармоник;

$g_n^m$, $h_n^m$, $q_n^m$, $s_n^m$ — сферические гармонические коэффициенты;

$P_n^m\left(\cos \theta \right)$ — присоединенные полиномы Лежандра в нормировке Шмидта.

Если мы строим модель магнитного поля только внутренних источников, то коэффициенты $q_n^m$ и $s_n^m$, описывающие поля, создаваемые источниками вне сферы, нас не интересуют, и члены с ними отбрасываются.

Составляющие напряженности поля по осям координат X, Y, Z находятся дифференцированием выражения (1) по соответствующей координате:

$Z=B_r=\frac{-\partial U}{\partial r}$; $X=B_{\theta }=\frac{-\partial U}{r\partial \theta }$; $Y=B_{\lambda }=\frac{-1}{r\sin \theta }\frac{\partial U}{\partial \lambda }$.

Если имеются измерения МПЗ в ряде точек Bi, то, рассматривая модельное поле как функцию коэффициентов $g_n^m$ и $h_n^m$, можно определить среднеквадратичное отклонение измеренных значений от модельных Bm:

$err=\sum {(Bi-Bm)}^2.$ (2)

Рассматривая err как функцию гармонических коэффициентов $g_n^m$ и $h_n^m$, находим значения коэффициентов, дающих минимальное значение отклонения.

Для учета изменчивости коэффициентов во времени используется два метода. Возможно выбрать для построения модели небольшой интервал времени, в течение которого коэффициенты можно считать постоянными, а затем, исследуя временные изменения самих коэффициентов, определить их вековой ход. Другой метод — записать коэффициенты в виде:

$gi=g{i}_0+k*\left(t-t_0\right),$

где gi₀ — значение коэффициента в момент времени t₀; t — текущее время. Достоинством первого метода является меньшее количество неизвестных, а недостатком — меньшее количество точек измерения (более короткий временной интервал для анализа). Во втором случае интервал сбора данных может быть бо́льшим, но и количество неизвестных существенно больше. Кроме того, первый метод не зависит от характера изменения коэффициентов, его можно описать любой удобной функцией, а во втором случае предполагается только линейный характер изменения.

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

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

Рис. 1. Покрытие измерениями SWARM поверхности земного шара за одни сутки.

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

Существует множество различных подходов к расчету коэффициентов гармоник, и международная модель геомагнитного поля International Geomagnetic Reference Field (IGRF) представляет собой комбинацию 9 моделей, разработанных независимо друг от друга. Такой подход позволяет проводить взаимное сравнение полученных моделей и с помощью их комбинаций получать наиболее достоверные значения гармонических коэффициентов. Детально процедура создания модели IGRF из отдельных субмоделей описана в [Alken et al., 2021a].

3. Выбор данных

Для расчета модели мы использовали векторные 1-sec данные со всех трех спутников SWARM за период январь 2014 г. — август 2019 г. из (https://earth.esa.int/web/guest/swarm/data-access), краткое описание методики расчетов было приведено в [Petrov and Bondar, 2021]. Предварительно данные были проверены, и все ошибочные данные были удалены. Данные обсерваторий и пунктов векового хода не использовались. Основная причина отказа от использования наземных данных — неоднородное пространственное распределение обсерваторий по земной поверхности.

Одним из общих методов уменьшения влияния ионосферных и магнитосферных токов является выбор таких временных интервалов, когда эти токи минимальны. Как правило, выбор основывается на различных индексах геомагнитной и солнечной активности, местном времени и состоянии освещенности ионосферы Солнцем. Каждая исследовательская группа использовала свой собственный метод уменьшения влияния ионосферных и магнитосферных токов и набор индексов [Olsen et al., 2000; Finlay et al., 2015], которые могли частично совпадать. В некоторых случаях вводились свои индексы [Chambodut et al., 2010] или даже вообще предварительный отбор данных не производился [Бондарь и др., 2006; Головков и Зверева, 1998]. Для выделения спокойных интервалов мы использовали следующий набор индексов: планетарный индекс Kp (http://wdc.kugi.kyotou.ac.jp/kp/index.html), индекс ASY/SYM и индекс авроральной геомагнитной активности AE (https://wdc.kugi.kyoto-u.ac.jp/aeasy/index.html). Для 2018—2019 гг., когда окончательные значения AE-индекса отсутствовали, рисунки предварительных значений AE из (https://wdc.kugi.kyoto-u.ac.jp/ae_realtime/index.html) были оцифрованы.

На первом этапе отбирались только 3-часовые интервалы с Кp менее 1.5 для текущего интервала и Кp менее 2.5 для предыдущего 3-часового интервала. Это позволило выбрать только спокойные интервалы с отсутствием заметной геомагнитной активности.

Затем из выбранных по Kp интервалов интервалы с вариацией SYMH более 20 нТл также были удалены. Эта процедура исключала влияние кольцевого тока и токов магнитопаузы. Окончательно при расчете из отобранных интервалов использовались только точки со значением AE менее 100 нТл. Это ограничивало влияние ионосферных и продольных токов.

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

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

4. Создание модели

Чтобы получить равномерное распределение точек в пространстве, поверхность Земли была разделена на шестиградусные (~660 км) широтные зоны, и каждый пояс был разделен на ячейки по долготе около 660 км шириной, так что в каждой зоне было несколько ячеек. Экваториальный пояс имел 60 ячеек, а самый высокоширотный — только 3, но площадь каждой ячейки была примерно одинаковой. Всего получилось 1146 ячеек. Затем все отобранные данные помещались в соответствующие ячейки в соответствии с координатами спутника. Обычно одни сутки наблюдений (без учета отбора данных по уровню геомагнитной активности) всех трех спутников позволяют заполнить примерно 70% ячеек, а двухсуточные данные — 95%. Принимая во внимание выбор данных только для спокойных интервалов, мы использовали 11-дневный интервал для расчета отдельной модели. Чтобы получить лучшее временное разрешение при моделировании, мы рассчитывали модели для каждого пятого дня, используя пять дней до и пять дней после даты. После заполнения ячеек использовался медианный фильтр для получения только одного значения для каждой компоненты, для каждой ячейки и для каждого спутника, поэтому в большинстве случаев для моделирования использовались 9 точек в каждой ячейке. Для нахождения коэффициентов вычислялись производные уравнения (2) по неизвестным $g_n^m$ и $h_n^m$, и они приравнивались нулю, далее полученная система линейных уравнений решалась стандартным методом Гаусса.

Коэффициенты любой сферической гармоники можно нанести на график как функцию даты. На рис. 2 показан пример такого графика. Точки, для которых было заполнено менее 90% всех ячеек, отбрасывались и в дальнейшем анализе не участвовали. Выбросы имеют несколько причин. Первая — ошибки в исходных данных SWARM. Предварительный отбор данных делался, но удалить все ошибки в таких больших массивах данных сложно. Вторая — вариации магнитного поля и, соответственно, получаемых моделей, связанные с изменением внешних и ионосферных источников. При построении модели использовались данные только спокойных интервалов, но какие-то возмущения могли оставаться и в такие моменты. Третья и, мы считаем, основная причина — неравномерное распределение данных по поверхности земного шара. Практически всегда в данных имеются более или менее длительные пропуски или данные исключались согласно описанным выше критериям, что приводило к отсутствию данных для отдельных регионов и изменению модели.

Рис. 2. Вариации коэффициентов сферических гармоник $g_1^0$ , $g_1^1$ , $h_1^1$  за период 2014—2020 гг.

5. Модель ИЗМИРАН для 2020 года

Описанный выше метод позволил рассчитать и построить графики изменения гармонических коэффициентов с начала 2014 г. до конца августа 2019 г. (модель требовалось представить до конца сентября) и экстраполировать их на 01.01.2020 г. Если коэффициент изменяется линейно, то эта экстраполяция не вызывает проблем. Однако при наличии заметной нелинейности для таких коэффициентов могут быть получены различные результаты, зависящие от используемой модели экстраполяции и интервала, выбранного для экстраполяции.

На рис. 2 показано изменение $g_1^0$, $g_1^1$ и $h_1^1$. Видно, что $g_1^0$ и $h_1^1$ изменяются практически линейно, аппроксимация данных параболой также показывает, что для $g_1^0$ и $h_1^1$ подгонка линейна (квадратичный член мал), но для $g_1^1$ нелинейность явно видна и линейная аппроксимация плохо подходит. Многие другие коэффициенты гармоник также демонстрируют видимую нелинейность, иногда очень сильную (см. $g_2^2$, $h_3^3$ на рис. 3). Некоторые гармонические коэффициенты даже не могут быть описаны параболической аппроксимацией и требуют кубической (рис. 3).

Рис. 3. Изменение коэффициентов сферических гармоник $g_2^2$  и $h_3^3$  по период 2014—2020 гг. Сплошные линии показывают параболическую аппроксимацию данных, пунктирные — кубическую.

Соответственно, использование разных методов аппроксимации даже при интерполяции на 4 месяца вперед может дать различные значения. Например, для показанной на рис. 3 гармоники $h_3^3$ различные методы аппроксимации при экстраполяции даже на 3 месяца дают результаты от –546 до –542 нТл. Для получения модели на 01.01.2020 г. мы решили использовать данные только за последние 2 года. Значения коэффициентов рассчитывались как параболически-экстраполированные данные за интервал 01.09.2017 г. — 31.08.2019 г. Для каждой сферической гармоники использовалась трехпроходная параболическая аппроксимация данных. После первой аппроксимации точки с отклонением более двух стандартных отклонений от полученной параболы были исключены, и следующая аппроксимация производилась на основе оставшихся данных. Эту процедуру повторяли 3 раза. На рис. 4 показаны все имеющиеся точки, но точки выше и ниже пунктирных линий в процессе аппроксимации были исключены, и окончательно парабола подбиралась только по точкам между пунктирными линиями.

Рис. 4. Пример отбора данных для вычисления гармонических коэффициентов модели.

Окончательная модель IGRF-2015 (до 2020 г. модель на 2015 г. считалась предварительной, как и модель 2020 г. будет считаться предварительной до 2025 г.) рассчитывалась таким же образом из аппроксимации для интервала 01.01.2014 г. — 31.12.2015 г., и в качестве модели использовались расчетные значения для 01.01.2015 г. Погрешности коэффициента гармоник вычислялись как среднеквадратичное отклонение от аппроксимирующей кривой, деленное на квадратный корень из числа использованных точек.

Основную сложность при построении моделей внутренних источников геомагнитного поля представляет исключение из данных в первую очередь магнитных полей ионосферных и затем магнитосферных токов. Поскольку величина ионосферных токов сильно зависит от проводимости ионосферы, определяемой в первую очередь освещенностью, то их влияние должно быть существенно больше в летнем полушарии, и их неполное ограничение должно привести к годовой вариации в гармонических коэффициентах. Рисунки 2 и 3 показывают, что в основных гармониках такой вариации нет. Просмотр графиков всех остальных гармоник также показал отсутствие таких вариаций, за исключением гармоник $g_2^0$, $g_3^0$, $g_4^0$. Наибольшую амплитуду годовая вариация имеет в гармонике , показанной на рис. 5.

Рис. 5. Изменения коэффициента $g_4^0$  разложения магнитного поля по гармоническим функциям за период 2014— 2020 гг.

Амплитуда этой вариации не превышает 1 нТл, а использование для определения параметров аппроксимации интервала точно в два года исключило ее влияние на значение $g_4^0$. Таким образом, можно утверждать, что влияние ионосферных токов очень мало, а методика практически исключает их влияние на получаемую модель.

6. Построение прогнозной модели годовой вариации

Модель IGRF требует в качестве значений годовой вариации магнитного поля (secular variation — SV) линейную аппроксимацию на период 2020—2025 гг., что легко получить в случае линейного изменения гармонического коэффициента, например $g_1^0$ или $h_1^1$. Гораздо сложнее оценить SV в случае нелинейных данных. Для некоторых гармоник (например, $h_3^3$, $h_2^2$, $g_3^2$) аппроксимации за 2015—2019 гг. и 2018—2019 гг. были существенно различны (SV значительно изменяются за пятилетний интервал), и экстраполяция двухлетнего набора данных на 3 года вперед очень ненадежна. В итоге SV 2022.5 были рассчитаны на основе 5-летнего (08.2014 г. — 08.2019 г.) параболического приближения. Этот выбор приводит к большей ошибке в начале периода 2020—2025 гг. по сравнению с оценкой за последние 2 года, но, поскольку нет никаких знаний о возможном характере изменений годовой вариации в будущем, мы полагаем, что такая консервативная оценка приведет к меньшей величине возможной ошибки.

7. Оценка точности модели

Любая модель, построенная по экспериментальным данным, неизбежно содержит ошибки, связанные с методикой построения модели, ошибкой исходных данных и методикой отбора данных. Если теоретические основы моделирования (методика представления потенциального поля в виде разложения по сферическим гармоникам) сомнений не вызывают, то остальные пункты не имеют однозначных критериев, и, естественно, модели, построенные различными авторами, несколько отличаются. В [Alken, 2021b] рассмотрены все 12 использованных для построения IGRF-моделей и показано, что разброс в значениях сферических гармонических коэффициентов, использованных при построении модели IGRF, имеет наибольшую величину в 3.5 нТл для $g_1^0$ и уменьшается до 0.1 нТл с ростом номера гармоник. Там же рассмотрено отклонение значения вертикальной компоненты Z отдельных моделей от усредненной модели IGRF, максимальное отклонение у большинства моделей имеет величину ~20 нТл и только у одной достигает 50 нТл. Построенная нами разность для Z-компоненты, рассчитанной по модели ИЗМИРАН и IGRF, имеет экстремальные значения ±17 и среднеквадратичное отклонение 5 нТл. Для X: соответственно –10; 8; 3, для Y: –10; 11; 3.5. В (https://www.ngdc.noaa.gov/IAGA/vmod/igrfhw.html ) оценена среднеквадратичная погрешность компонент магнитного поля, связанная с ошибками коэффициентов современных моделей IGRF величиной в 10 нТл, наши оценки дают даже меньшую величину. Сравнение модельных данных с реальными измерениями в обсерваториях большого смысла не имеет, так как реальные данные содержат значительный вклад коровых источников (аномальное поле), но в ряде работ такое сравнение делалось, например [Ciarán and Beggan, 1971; Golovkov et al., 1997], и для IGRF-13 получено стандартное отклонение 87 нТл для X, 73 нТл для Y, 114 нТл для Z, для магнитного склонения — 0.4°. Поскольку модели типа IGRF определяют низшие гармоники моделей высокого разрешения (EMM, HGDM и др.), в которых поля коровых источников в значительной степени учтены, именно погрешность компонент магнитного поля, связанная с ошибками коэффициентов современных моделей IGRF, в значительной степени определяет и точность моделей высокого разрешения. Практически значение магнитного склонения D наиболее важно, на рис. 6 показано распределение разности склонения, рассчитанной по модели ИЗМИРАН и IGRF. Экстремальные значения разности D равны –13 и 12, стандартное отклонение — 1.6 мин, причем экстремальные значения наблюдаются на высоких (>70°) широтах, а на территории России эта разность не превышает двух минут.

Рис. 6. Изолинии распределения разности D-компонент, рассчитанной по модели ИЗМИРАН и IGRF.

8. Заключение

В конце 2019 г. модели главного магнитного поля Земли были созданы несколькими организациями независимо и различными методами. В декабре 2019 г. рабочей группой V Международной ассоциации по геомагнетизму и аэрономии эти модели были объединены в единую модель, ставшую стандартной моделью IGRF-2020 (или IGRF-13, по порядковому номеру модели) и действующую на период до 2025 г. ИЗМИРАН представил свои варианты для моделей DGRF-2015, IGRF-2020 и SV-2022.5, они соответствовали требованиям рабочей группы и были включены в итоговую модель IGRF-13. Сравнение модели ИЗМИРАН с усредненной моделью IGRF показывает, что различие между ними не значимо и модель ИЗМИРАН, рассчитанная на промежуточные даты между 2020—2025 гг., может применяться вместо модели IGRF.