Геомагнетизм и аэрономия, 2020, T. 60, № 5, стр. 588-599

Выделение вихревых токовых структур в ионосфере и оценка их параметров по наземным магнитным данным

В. Е. Чинкин 1*, А. А. Соловьев 12, В. А. Пилипенко 12

1 Геофизический центр Российской академии наук (ГЦ РАН)
г. Москва, Россия

2 Институт физики Земли им. О.Ю. Шмидта Российской академии наук (ИФЗ РАН)
г. Москва, Россия

* E-mail: v.chinkin@gcras.ru

Поступила в редакцию 03.12.2019
После доработки 31.01.2020
Принята к публикации 21.05.2020

Полный текст (PDF)

Аннотация

Предложена система обработки данных 2D-сети магнитных станций с целью выделения вихревых токовых структур в ионосфере и оценки их характерных параметров. В качестве примера методика применена к анализу структуры дневных конвективных вихрей TCV (traveling convection vortices) по данным арктических станций. Задача решается с применением методов оптимизации различных функций, полученных на основе пространственной интерполяции и последующей регуляризации данных. Разработанный подход дает возможность не только автоматически находить несколько вихревых структур, но и определять текущие значения их характерных параметров: пространственную структуру продольных токов и групповую скорость горизонтального распространения вихря вдоль ионосферы.

1. ВВЕДЕНИЕ

Перенос энергии электромагнитных возмущений из магнитосферы Земли к ионосфере происходит преимущественно за счет продольных (вдоль силовых линий магнитосферного магнитного поля) токов. В высоких широтах, где геомагнитное поле В почти вертикально, наземное магнитное возмущение создается системой вихревых холловских токов, возбуждаемых продольными токами. Вихревые ионосферные структуры могут иметь самые различные масштабы, в зависимости от типа магнитосферного возмущения: от планетарных масштабов при воздействии межпланетной ударной волны [Fujita et al., 2005] до малых кинетических масштабов в авроральных структурах [Chmyrev et al., 1988]. В локализованных структурах, связанных с вихревыми токами, сосредоточена основная энергия нестационарных магнитосферно-ионосферных возмущений. Именно они ответственны за всплески интенсивности геоиндуцированных токов в протяженных линиях электропередач [Ngwira et al., 2015; Belakhovsky et al., 2019]. Возможность автоматического выделения локализованных вихревых структур и определение их параметров по данным наземных магнитометров является крайне важной задачей.

Как правило, наличие вихревых структур в возмущениях геомагнитного поля определяется визуально по картине эквивалентных ионосферных токов, построенной по магнитным данным (например, [Engebretson et al., 2013]). В данной работе предложен новый метод, который не только дает возможность автоматически находить даже несколько одновременных вихревых структур в ионосфере по данным 2D-сети магнитометров, но и определять их характерные параметры.

Для апробации метода выделения вихревых структур рассмотрены конвективные холловские вихри (в английской литературе “TCV – traveling convection vortices”) [Glassmeier, 1992]. Конвективные вихри в ионосфере являются хорошо известным и наглядным проявлением импульсного воздействия солнечного ветра на магнитосферу [Kivelson and Southwood, 1991; Sibeck and Korotova, 2000; Kataoka et al., 2002]. Наземным отражением подобных событий служат изолированные импульсные возмущения магнитного поля с характерными временными масштабами ~10 мин и амплитудами до первых сотен nT (в англоязычной литературе MIE – “magnetic impulse events”) [Luhr and Brawert, 1994; Lanzerotti et al., 1987, 1990]. Спорадические MIE/TCV-импульсы являются одним из нестационарных явлений, характерных для дневной высокоширотной области. В качестве примера будет рассмотрен типичный дневной конвективный вихрь по данным арктических станций.

2. МОДЕЛЬ ИОНОСФЕРНОГО КОНВЕКТИВНОГО ВИХРЯ

Продольные токи j| |, протекающие вдоль силовых линий магнитосферного магнитного поля между магнитосферой и ионосферой, замыкаются в проводящем Е-слое ионосферы на систему поперечных педерсеновских JP и холловских JH токов

$J = {{J}_{P}} + {{J}_{H}} = {{\Sigma }_{P}}E + {{\Sigma }_{H}}[e \times E].$
Здесь J – горизонтальный ионосферный ток, проинтегрированный по толщине проводящего слоя ионосферы; Е – горизонтальное электрическое поле в ионосфере; е – нормаль к ионосфере. Соотношение между педерсеновским и холловским токами определяется отношением интегральных (по толщине Е-слоя ионосферы) педерсеновской ΣP и холловской ΣH проводимостей. С математической точки зрения, педерсеновский JP и холловский JH токи являются безвихревой (дивергентной) и бездивергентной (вихревой) составляющими вектора горизонтального тока J. Используя двумерные операторы Rot A = rotzA и Div A = ∂Ax/∂x +Ay/∂y, эти свойства на плоскости ионосферы можно представить в виде Rot JP = 0 и Div JH = 0. В высоких широтах, где геомагнитное поле практически вертикально и j||jz, выполняется условие непрерывности токов
${{j}_{z}} = - {\text{Div}}{{J}_{P}}.$
Это соотношение показывает, что продольный ток jz в горизонтально однородной ионосфере замыкается только педерсеновским током. При вертикальном геомагнитном поле суммарный магнитный эффект под ионосферой от продольного и педерсеновских токов равен 0, поэтому наземное магнитное возмущение создается только системой ионосферных холловских токов. Таким образом, эквивалентные токи в горизонтально-однородной ионосфере, определяемые по наземным магнитным данным, являются холловскими токами. В области ионосферы с горизонтально однородными интегральными проводимостями
(1)
${{j}_{z}} = ({{{{\Sigma }_{P}}} \mathord{\left/ {\vphantom {{{{\Sigma }_{P}}} {{{\Sigma }_{H}}}}} \right. \kern-0em} {{{\Sigma }_{H}}}}){\text{Rot}}{{J}_{H}}.$
Отношение ΣНР ~ 1.5–2.0 мало меняется даже при сильных вариациях ионосферной проводимости. Источник вихревых систем, создающий импульс TCV на земной поверхности, обычно представляют в виде либо (а) азимутально симметричной токовой трубки, либо (б) пары противоположно направленных токов [Glassmeier, 1992]. Нами используется модель (а), изображенная на рис. 1. На рисунке 2 качественно продемонстрированы графики пространственной структуры горизонтальной Br(R) и вертикальной Bz(R) составляющих магнитного поля, создаваемых однородной токовой трубкой, которое наблюдалось бы на наземных станциях под горизонтально-однородной ионосферой в высоких широтах. Педерсеновские токи растекаются симметрично от центра втекающего продольного тока, но магнитного отклика на Земле они не возбуждают, так как он компенсируется полем, создаваемым продольным током (Bφ= 0).

Рис. 1.

Схематическое изображение продольных токов (J||), втекающих и вытекающих из ионосферы, и их замыкания на поперечные педерсеновские (JP) и холловские (JH) токи в проводящем E-слое ионосферы для магнитосферного возмущения в виде азимутально симметричной токовой трубки.

Рис. 2.

Безразмерные радиальные функции Br(R) и Bz(R) составляющих магнитного поля и производной ∂zBz(R) для модели вихря, показанной на рис. 1. Отношение высоты ионосферы к радиусу контура h0/Rc = 3/16. Ось y начала координат совпадает с центром вихря, штриховые линии обозначают его границы.

Качественно представим, что на расстоянии порядка радиуса Rc величина педерсеновских токов спадает до пренебрежимо малой величины. В свою очередь холловские токи заполняют область от r = 0 до r = Rc. Создаваемое ими магнитное поле Br имеет биполярную форму (обращается в 0 под центром вихревой токовой системы). По компоненте Bz модельный TCV-импульс представляет собой униполярное возмущение магнитного поля. Качественно картина магнитного возмущения, представленная на рис. 1, соответствует численной модели TCV [McHenry and Clauer, 1987]. Ниже опишем методику приближенного определения положения центров подобных локализованных токовых структур. В качестве оценки характерного размера токового контура, берется линия ∂Bz/∂z = 0, хотя в дальнейшем ∂Bz/∂z не вычисляется, а приводится здесь лишь для демонстрации физического смысла данной модели TCV.

3. ПОСТРОЕНИЕ ВЕКТОРНОГО ПОЛЯ ЭКВИВАЛЕНТНОГО ИОНОСФЕРНОГО ТОКА

В качестве исходных используются данные 2D магнитной сети, представляющие собой компоненты горизонтального вектора магнитной индукции Bi(k) = (Xi(k), Yi(k))T, i = 1, , Nst, где k – номер обрабатываемой минуты, зарегистрированные на Nst станций. Для удобства географические координаты станций (lati, loni)T записываются в системе коширота-долгота (ψi, φi)T = (π/2 – lati, loni)T. В данных могут содержаться пропуски во времени; в случае, если пропуск меньше 10 с, этот пропуск во времени линейно интерполируется, если данные пропущены на концах временного интервала, то интерполяция производится методом ближайшего значения. В исходных данных также устраняются постоянная составляющая, суточная вариация и продолжительные длительные вариации. Фильтрация данных осуществляется с помощью фильтра высоких частот с конечной импульсной характеристикой. Вид магнитной записи с удаленным низкочастотным трендом показан на рис. 3. Рассматриваемый в качестве примера импульс TCV, зарегистрированный 31 января 1997 г. на меридиональном профиле магнитных станций, показан на рис. 4.

Рис. 3.

Примеры фильтрации и удаления тренда магнитограмм по трем компонентам станции PBQ (Пост-де-ля-Белен, Канада) за 31 января 1997 г.: Ориг. – исходный сигнал, Фильт. – отфильтрованный сигнал.

Рис. 4.

Импульс TCV, зарегистрированный на меридиональном профиле магнитных станций (геомагнитные компоненты X, Y, Z). У правой оси ординат даны коды станций и их геомагнитные координаты.

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

${{({{x}_{i}}_{{,k}},{{y}_{i}}_{{,k}})}^{T}} = M({{\varphi }_{0}} + {{\varphi }_{k}}){{({{\psi }_{i}}{\text{cos}}({{\varphi }_{i}}),{{\psi }_{i}}{\text{sin}}({{\varphi }_{i}}))}^{T}},$
где $M(\varphi )$ = $\left( {\begin{array}{*{20}{c}} {\cos (\varphi )}&{ - \sin (\varphi )} \\ {\sin (\varphi )}&{\cos (\varphi )} \end{array}} \right)$ – матрица поворота; φ0 = 2π(t0/1440) – угол, соответствующий повороту Земли на момент времени t0, с которого начинаются ряды анализируемых данных, φk = 2π(k – 1)/ 1440, верхний индекс Т означает транспонированную матрицу. Для рассматриваемых региональных возмущений с масштабами менее 103 km, сферичность Земли можно не учитывать.

Составляющие вектора индукции магнитного поля i-й станции в k-й момент времени Bi(k) записаны в полярной системе: ось Xφ, ψ направлена на полюс, базисный вектор оси Yφ, ψ определяется по стандартному правилу ey = ez × ex. Тогда для каждого момента времени k координаты преобразуются с помощью матрицы перехода Ui,k)

$U({{\varphi }_{{i,k}}}) = \left( {\begin{array}{*{20}{c}} { - \cos ({{\varphi }_{{i,k}}})}&{ - \sin ({{\varphi }_{{i,k}}})} \\ { - \sin ({{\varphi }_{{i,k}}})}&{\cos ({{\varphi }_{{i,k}}})} \end{array}} \right),$
$\begin{gathered} {\text{cos(}}{{\varphi }_{i}}{{_{,}}_{k}}{\text{)}} = {{{{x}_{i}}(k)} \mathord{\left/ {\vphantom {{{{x}_{i}}(k)} {{{{[{{x}_{i}}{{{_{,}}}_{k}} + {{y}_{i}}{{{_{,}}}_{k}}]}}^{{{1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0em} 2}}}}}}} \right. \kern-0em} {{{{[{{x}_{i}}{{{_{,}}}_{k}} + {{y}_{i}}{{{_{,}}}_{k}}]}}^{{{1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0em} 2}}}}}}, \\ {\text{sin(}}{{\varphi }_{i}}{{_{,}}_{k}}{\text{)}} = {{{{y}_{i}}{{{_{,}}}_{k}}} \mathord{\left/ {\vphantom {{{{y}_{i}}{{{_{,}}}_{k}}} {{{{[{{x}_{i}}{{{_{,}}}_{k}} + {{y}_{i}}{{{_{,}}}_{k}}]}}^{{{1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0em} 2}}}}.}}} \right. \kern-0em} {{{{[{{x}_{i}}{{{_{,}}}_{k}} + {{y}_{i}}{{{_{,}}}_{k}}]}}^{{{1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0em} 2}}}}.}} \\ \end{gathered} $

По множеству точек данных {xi,k, yi,k, Bi(k)} стандартной Matlab-функцией griddata решается задача интерполяции с неравномерно расположенными точками данных с применением кубических полиномов. Регулярная сетка задается следующим образом: xi = xi – 1 + Δx, yj = yj – 1+ Δy, x0= xmin, y0= ymin, шаг по x и y, соответственно, Δx = (xminxmax)/(Nx – 1), Δy = (yminymax)/(Ny – 1), где xmin, xmax, ymin, ymax – границы рассматриваемого участка; Nx, Ny – число узлов сетки. Таким образом, для каждого момента времени k получаются матрицы: Bk(i, j) = (Bx,k(i, j), By,k(i, j))T, i = 1, , Nx, j = 1, , Ny. В данной работе из-за удобства преобразования координат при вращении использовались географические координаты станций. Однако в любой момент координаты могут быть преобразованы в геомагнитные.

Матрицы Bk(i, j) содержат в себе дискретное 2D векторное поле вариаций магнитного поля B(t, x, y) на поверхности Земли. Поле эквивалентного ионосферного тока J можно представить путем поворота векторов наземного поля B на угол π/2 вокруг базисного вектора ez. Таким образом, мы вводим дискретное плоское векторное поле эквивалентного ионосферного тока Jk(i, j) = = (–By,k(i, j), Bx,k(i, j))T для каждого шага k. Векторное поле эквивалентного ионосферного тока J для момента времени 31 января 1997 г., 14:53 UT, показано на рис. 5.

Рис. 5.

Плоское векторное поле эквивалентного ионосферного тока J, построенное для события 31 января 1997 г. 14:53 UT до коррекции. Направление на Солнце – слева, длина стрелок зависит от интенсивности магнитного возмущения. Жирные стрелки отображают реальные вектора эквивалентного тока в точках расположения станций.

Однако полученное таким образом векторное поле J может не быть бездивергентным. Поэтому, чтобы построить ассоциированное с холловскими токами бездивергентное поле J°, восстановим плоскую безвихревую составляющую магнитного поля. Для этого нам необходимо представить анализируемое дискретное магнитное поле B(x, y) как градиент от потенциальной функции и определить эту функцию. Для этого введем функцию ρ(xi, yj), которую можно интерпретировать как “заряд” поля горизонтальных компонент:

(2)
$\rho (x,y) = {\text{Div}}{{B}_{ \bot }}.$
Применим к ρ(xi, yj) сглаживающий фильтр Гаусса для устранения шумов, появившихся при дифференцировании. Выразим потенциал магнитного поля по следующей упрощенной схеме:
${{\varphi }_{k}}(x,y) = \frac{1}{{4\pi }}\sum\limits_{i,j}^{{{N}_{x}},{{N}_{y}}} {\left\{ {\begin{array}{*{20}{c}} {\frac{{\rho ({{x}_{i}},{{y}_{j}})}}{{\sqrt {{{{({{x}_{i}} - x)}}^{2}} + {{{({{y}_{j}} - y)}}^{2}}} }}\Delta x\Delta y} \\ {0,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,(x,y) = ({{x}_{i}},{{y}_{j}}).} \end{array}} \right.} $
Тогда скорректированное векторное поле можно вычислить следующим образом:
$B_{ \bot }^{^\circ } = {\text{Grad}}(C\varphi ).$
При одновременном анализе множества событий описанный подход может приводить к занижению оценки модуля поля B по сравнению с его реальными значениями. Для частичного устранения этого недостатка введена дополнительная масштабирующая константа С, относительно которой решается дополнительная оптимизационная задача:
$\sum\limits_{i = 1}^{{{N}_{{st}}}} {{{{\left| {{{B}_{i}}(k) - {{{\left. {{\text{Grad}}(C\varphi )} \right|}}_{{{{x}_{i}}(k),{{y}_{j}}(k)}}}} \right|}}^{2}}} \to \min .$
Пример скорректированного бездивергентного поля эквивалентного ионосферного тока J°, рассчитанного через $B_{ \bot }^{^\circ },$ показан на рис. 6 (на рис. 5 оно нескорректированное). Это полученное бездивергентное поле мы ассоциируем с JH:
(3)
${\mathbf{J}}^\circ \approx {{\mu }_{0}}{{{\mathbf{J}}}_{H}}.$
Именно это поле будет использовано в дальнейшем для выделения вихрей и символ “°” будет опускаться.

Рис. 6.

Пример выделения вихревых структур для события 31 января 1997 г. 14:53 UT: отображение плоского скорректированного векторного поля эквивалентного ионосферного тока J и величины плотности модуля продольного тока |Jz|. Направлению вихря против часовой стрелки соответствует маркер кружочка его центра и вытекающий ток jz < 0, направлению по часовой стрелке соответствуют маркер звездочи центра и втекающий ток jz > 0.

5. ОПРЕДЕЛЕНИЕ СТРУКТУРЫ И ХАРАКТЕРИСТИК ВИХРЕЙ

5.1. Поиск центра вихря

При автоматическом определении центров вихрей центр магнитного возмущения предлагается оценивать по экстремумам следующей функции:

(4)
$G(x,y) = {\text{Rot}}({{\mathbf{J}} \mathord{\left/ {\vphantom {{\mathbf{J}} {\left| J \right|}}} \right. \kern-0em} {\left| J \right|}}).$
Экстремумы функции (4) соответствуют максимумам продольного тока, втекающего или вытекающего в ионосферу, согласно (1). Причем, из связи ротора и циркуляции следует, что она же является максимальной циркуляцией минимальных контуров. Нормировка введена для того, чтобы сравнительно небольшие вихри не терялись на фоне больших структур. Поэтому каждый экстремум функции (4) может быть признан характерной внутренней точкой некоторого токового контура в ионосфере.

Численное вычисление функции G связано с вычислением дивергенции поля B/|B|, которое выполняется стандартной Matlab-функцией divergence. Так как после дифференцирования значительно усиливаются низкочастотные шумы, G сглаживается пространственным фильтром Гаусса. Поиск экстремумов сглаженной функции (4) выполняется методом перебора с условием

$\begin{array}{*{20}{c}} {G({{x}_{i}},{{y}_{{j + 1}}})\left\langle {G({{x}_{i}},{{y}_{j}})} \right\rangle G({{x}_{i}},{{y}_{{j - 1}}})} \\ {G({{x}_{{i - 1}}},{{y}_{j}})\left\langle {G({{x}_{i}},{{y}_{j}})} \right\rangle G({{x}_{{i + 1}}},{{y}_{j}})} \end{array}.$
Чтобы исключить большинство случайных максимумов, применен критерий на среднеквадратическое отклонение
(5)
$\left| {{\text{Rot}}({\mathbf{J}})} \right| > \chi \sqrt {M[{\text{Rot}}{{{({\mathbf{J}})}}^{2}}]} ,$
где M[Rot(J)2] – среднее квадратов роторов по всему полю; χ – настраиваемый параметр. Эмпирически установлено, что наилучший результат достигается при χ = 2. Вместо формулы (5) можно было бы использовать жесткое задание уровня отсечки |Rot(J)| > K, где K – значение отсечки, но такая структура не отражала бы особенности двумерной поверхности, которая задается |Rot(J)|, и параметр K пришлось бы подбирать для каждого события свой. В предложенном нами подходе параметр χ намного меньше подвержен перенастройке. Таким образом, можно получить первое оценочное множество центров вихрей ri(k) = = (xi(k), yi(k))T. Пример данной техники можно видеть на рис. 6, где были распознаны 4 вихревые структуры. Кружочком и звездочками отмечены центры с втекающим и вытекающим продольными токами, соответственно.

5.2. Оценка масштаба и интенсивности вихрей

Пусть точка r0 = (x0, y0)T – найденный центр вихря, контур Γ(R) – окружность радиуса R с центром в точке r0, а S – поверхность, образованная этим контуром. Если рассмотреть циркуляцию вектора J по контуру Γ(R)

$\begin{gathered} C(\Gamma (R)) = \oint\limits_\Gamma {Jdl} = \\ = \int\limits_S {{\text{Rot}}J} dS \sim \int\limits_S {{\text{Div}}{{B}_{ \bot }}} = - \int\limits_S {\frac{{\partial {{B}_{z}}}}{{\partial z}}\,} dS, \\ \end{gathered} $
то для плоского приближения ей соответствует частная производная потока Φ по перпендикулярному поверхности Земли направлению через поверхность S:
$C(\Gamma (R)) = - \frac{{\partial \Phi }}{{\partial z}}.$
Считая эти формулы физической интерпретацией метода, подчеркнем, что Φ не вычисляется, а выкладки приведены для указания связи с нашей численной моделью, показанной на рис. 2. Масштаб вихревой структуры можно оценить путем нахождения экстремума циркуляции:
$\frac{{\partial C(\Gamma (R))}}{{\partial R}} = 0.$
Значение максимальной циркуляции эквивалентного тока будет соответствовать масштабу максимального изменения потока магнитной индукции в вертикальном направлении. Необходимо подчеркнуть, что такую систему легко обобщить для случая любого многопараметрического контура. В качестве первого приближения естественно предложить те контуры, которые формируют сечения Земли, так как для их масштабирования требуется лишь один параметр. На рисунке 7 продемонстрированы графики C(Г(R)) для соответствующих вихрей, которые были определены на предыдущем шаге (рис. 6). На них видны экстремумы, по которым оцениваются масштабы этих вихрей.

Рис. 7.

Зависимость нормированной на радиус земли циркуляции от радиуса для каждого распознанного вихря на рис. 6. Координаты вихрей указанны в системе коширота–долгота (ψ, φ).

На практике циркуляция вычисляется численным интегрированием контурного интеграла

$\begin{gathered} \oint\limits_\Gamma {Jdl} = \oint\limits_\Gamma {[{{B}_{x}}dx + {{B}_{y}}dy]} \approx \\ \approx R\sum\limits_{i = 1}^N {[{{B}_{y}}({{x}_{i}},{{y}_{i}})\sin ({{\varphi }_{i}}) + {{B}_{x}}({{x}_{i}},{{y}_{i}})\cos ({{\varphi }_{i}})]\Delta \varphi } , \\ \end{gathered} $
где xi = R cos(φi), yi = R sin(φi); значения Bx(xi, yj), By(xi, yj) находятся билинейным интерполированием матрицы Bk(xi, yj); φi = Δφ(i – 1); Δφ = 2π/N – длина шага; N – число разбиений, определяющее точность вычисления. Плотность продольного тока jz, характеризующую интенсивность рассматриваемого вихря в заданный момент времени, будем численно искать в узлах сетки, исходя из выражений (1) и (3):
${{j}_{z}}(x,y) = ({1 \mathord{\left/ {\vphantom {1 {{{\mu }_{0}}}}} \right. \kern-0em} {{{\mu }_{0}}}}){\text{Rot}}{\mathbf{J}}.$
Так как отношение проводимостей изменяется в течение суток слабо, мы полагаем ΣPH = 1. Пример данной техники оценки масштаба можно видеть на рис. 6 и 7.

5.3. Определение траектории вихря

С помощью методики оценки масштабов вихрей каждому центру вихря ri поставлен в соответствие его масштаб ri = (xi(k), yi(k))TRi(k), где i – порядковый номер вихря в момент времени k. Для того чтобы связать разные кадры наблюдений, введем метрику между распознанными структурами с индексом i момента времени k и точкой j момента времени k – 1:

$\rho (i,k{\text{|}}\,j,k - 1) = \exp \left( { - \frac{{{{\Delta }_{{i,j}}}(k)}}{{{{\sigma }_{1}}}} - \frac{{\Delta {{R}_{{i,j}}}(k)}}{{{{\sigma }_{2}}}}} \right),$
где
$\begin{gathered} {{\Delta }_{{i,j}}}(k) = \left| {{{r}_{i}}(k) - {{r}_{j}}(k - 1)} \right|, \\ \Delta {{R}_{{i,j}}}(k) = \left| {{{R}_{i}}(k) - {{R}_{j}}(k - 1)} \right|, \\ \end{gathered} $
σ1, σ2 – эмпирически подбираемые параметры. Тогда индекс ближайшей к ( j, k – 1) точке будет определяться как $i_{j}^{^\circ }(k)$ = $\arg \left\{ {\mathop {\max }\limits_i (\rho (i,k{\text{|}}{\kern 1pt} j,k - 1))} \right\}.$ Кроме того, на этом этапе возможно дополнительно исключить из рассмотрения точки, которые находятся слишком далеко от остальных согласно введенной метрике. Связные точки, отражающие траекторию движения одного вихря, должны удовлетворять условию:
$\rho \left( {i_{j}^{^\circ }(k),k{\text{|}}\,j,k - 1} \right) < {{d}_{0}},$
где d0 – эмпирический коэффициент, зависящий от σ1 и σ2. Реализацию данной процедуры можно наблюдать на рис. 8. Теперь можно говорить о том, что мы можем проводить наблюдения от момента зарождения вихря до момента его исчезновения или до момента, когда структура выйдет из области видимости сети магнитометров. Пороговый уровень, когда можно говорить о зарождении/исчезновении вихря, определяется исходя из соотношения (5): он уникален для каждого момента времени и определяется исходя из общей геомагнитной обстановки и эмпирически подобранного параметра χ.

Рис. 8.

Пример определения траекторий движения вихревых структур во время события 31 января 1997 г. Линии указывают перемещение центров двух вихрей (с номерами 2 и 4 на рис. 6) за наблюдаемый период, начало движения помечено кружочком. На рисунке указаны время зарождения Tb, время существования D и длина пути L каждого вихря. Для наглядности приведено положение вихрей на 14:55 UT, обозначения аналогичны рис. 6. На рисунке отмечены IAGA-коды станций.

6. ДИНАМИЧЕСКИЕ ХАРАКТЕРИСТИКИ ВИХРЕЙ ВО ВРЕМЯ СОБЫТИЯ 31.01.1997 г.

Продемонстрируем работу предложенной методики для оценки динамических параметров дневных TCV-вихрей за 31 января 1997 г. Соответствующий импульс (см. рис. 4) представляет собой типичное TCV-событие в предполуденные часы. Исходными данными служат геомагнитные наблюдения станций сетей CARISMA и MACCS, полученные через портал SuperMAG. Предварительная обработка данных включает фильтрацию, интерполяцию и коррекцию.

На рисунке 8 продемонстрированы траектории движения центров двух наиболее продолжительных возмущений, имевших место с 14:50 до 14:58 UT, с указанием времени зарождения каждого вихря, длительности и пройденного расстояния за время его существования. Характерное время существования анализируемых вихрей составляет ~10 min. Из рисунка видно, что вихри за время своего существования смещаются к ночной стороне и к полюсу Земли. Анти-солнечное направление распространения вызвано их “сдуванием” с дневной стороны потоком солнечного ветра. Характерные скорости смещения центров вихрей составляют ~3.9 km/s и ~1.3 km/s. Для сравнения, линейная скорость вращения Земли на рассматриваемых широтах (~70° N) составляет порядка 0.16 km/s, что свидетельствует о том, что дрейфующие в западном направлении вихри на утренней стороне движутся относительно поверхности Земли.

Графики на рис. 9 показывают изменения во времени сглаженных характеристик вихревых возмущений: значение продольного тока Jz в центре каждого вихря, оценочный радиус вихря и циркуляцию тока по этому радиусу. Характерный размер вихря составляет ~500 km, порядок величины продольного тока до ~0.2 A/km2 при величинах возмущения геомагнитного поля ~150 nT.

Рис. 9.

Эволюция характерных параметров двух наиболее продолжительных вихрей (с номерами 2 и 5 на рис. 6) за 31 января 1997 г., включая интенсивность продольного тока Jz в центре вихря, радиус вихря и циркуляцию тока по радиусу. По оси абсцисс указано время в минутах с момента зарождения каждого вихря. Типы линий графиков соответствуют типам траекторий вихрей на рис. 8.

Представленные графики отчетливо демонстрируют прямую зависимость между величиной продольного тока Jz и циркуляцией по контуру (рис. 9, первый и третий графики). Необходимо отметить, что определять параметры вихря на последних минутах жизни зачастую сложно из-за низких значений продольного тока. Кроме того, иногда вихрь пропадает из зоны видимости, входя в слепую зону сети станций. Так, в рассматриваемом случае отрицательный вихрь (штрихпунктирные линии на рис. 8 и 9) пропадает из видимости почти на пике величины продольного тока.

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

7. ОБСУЖДЕНИЕ

Реализованный в работе метод выделения вихревых структур является простой альтернативой сложному методу сферически-гармонического анализа [Amm and Viljanen, 1999]. При сравнении полученных результатов можно заметить лишь незначительные отличия результатов, даже с учетом того, что изложенный метод не учитывает сферичность Земли. Еще одно заметное преимущество заключается в том, что алгоритм не требует большого числа станций по всей земной сфере, и для его работы достаточно только данных в области решения.

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

Наблюдения на региональных сетях магнитометров показывают наличие эффектов горизонтального распространения TCV вдоль ионосферы. Обычно для оценки эффектов распространения геомагнитных возмущений сопоставляют магнитограммы станций, расположенных либо вдоль геомагнитного меридиана, либо параллели. В долготном направлении наблюдаемые сдвиги по времени между импульсами на разнесенных станциях, как правило, соответствовали фазовым скоростям ~10 km/s в анти-солнечном направлении. В проекции на экваториальную плоскость магнитосферы это значение примерно соответствует характерной скорости обтекания флангов магнитопаузы потоком плазмы магнитослоя. В широтном направлении наблюдаемые сдвиги во времени указывали на кажущееся распространение к полюсу с фазовой скоростью ~1 km/s. Широтное распространение связано с тем, что сигнал распространяется вдоль магнитных силовых линий от экваториальной области магнитосферы к высокоширотной ионосфере. Поскольку время распространения растет с увеличением широты Φ, на земной поверхности будет наблюдаться запаздывание, соответствующее кажущемуся распространению к полюсу. Предложенная же методика определяет движение центра вихря, т.е. определяет групповую скорость возмущения. Для анализируемого события скорость движения вдоль географической широты составляет 4.1 km/s, вдоль географической долготы – 0.3 km/s, т.е. оказывается меньше типичных фазовых скоростей.

Сейчас можно сказать, что даже с учетом сделанных приближений существующий алгоритм справляется со своей задачей вполне удовлетворительно. В частности, метод может быть использован для анализа большого числа событий за длительный интервал времени с последующим статистическим анализом [Clauer and Petrov, 2002].

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

Представленный набор алгоритмов выполняет целостную автоматическую обработку данных, начиная с анализа исходных измерений магнитного поля Земли до конечного результата вычислений ряда динамических характеристик вихревых токовых структур. Предложенный метод позволяет численно определять текущие параметры вихревых структур: положение центра и размер вихревого образования, траекторию и скорость движения, а также направление и интенсивность продольного тока. Работа всей системы продемонстрирована на примере изучения TCV, зарегистрированного на 2D-сети станций CARISMA и MACCS 31.01.1997 г. Предложенные алгоритмы позволяют получать вполне удовлетворительные результаты и могут быть широко использованы исследователями в области солнечно-земной физики и геофизики.

Список литературы

  1. Amm O., Viljanen A. Ionospheric disturbance magnetic field continuation from the ground to the ionosphere using spherical elementary current systems // Earth Planets Space. V. 51. P. 431–440. 1999.

  2. Belakhovsky V., Pilipenko V., Engebretson M., Sakharov Ya., Selivanov V. Impulsive disturbances of the geomagnetic field as a cause of induced currents of electric power lines // J. Space Weather Spac. V. 9. № A18. 2019. https://doi.org/10.1051/swsc/2019015

  3. Chmyrev V.M., Bilichenko S.V., Pokhotelov O.A., Marchenko V.A., Stenflo L. Alfven vortices and related phenomena in the ionosphere and magnetosphere // Physica Scripta. V. 38. № 6. P. 841–854. 1988.

  4. Clauer C.R., Petrov V.G. A statistical investigation of traveling convection vortices observed by the west coast Greenland magnetometer chain // J. Geophys. Res. V. 107. A7. 2002. https://doi.org/10.1029/2001JA000228

  5. Engebretson M.J., Yeoman T.K., Oksavik K. et al. Multi-instrument observations from Svalbard of a traveling convection vortex, electromagnetic ion cyclotron wave burst, and proton precipitation associated with a bow shock instability // J. Geophys. Res. V. 118. № 6. P. 2975–2997. 2013. https://doi.org/10.1002/jgra.50291

  6. Fujita S., Tanaka T., Motoba T. A numerical simulation of the geomagnetic sudden commencement: 3. SC in the magnetosphere-ionosphere compound system // J. Geophys. Res. V. 110. № A11203. 2005. https://doi.org/10.1029/2005JA011055

  7. Glassmeier K.-H. Traveling magnetospheric convection twin-vortices: observations and theory // Ann. Geophysicae. V. 10. № 8. P. 547–565. 1992.

  8. Kataoka R., Fukunishi H., Lanzerotti L.J. et al. Traveling convection vortices induced by solar wind tangential discontinuities // J. Geophys. Res. V. 107. № A12. P. 1455. 2002. https://doi.org/10.1029/2002JA009459

  9. Kivelson M.G., Southwood D.J. Ionospheric travelling vortex generation by solar wind buffeting of the magnetosphere // J. Geophys. Res. V. 96. № A2. P. 1661–1667. 1991.

  10. Lanzerotti L.J., Wolfe A., Trivedi N., Maclennan C.G., Medford L.V. Magnetic impulse events at high latitudes: Magnetopause and boundary layer plasma processes // J. Geophys. Res. V. 95. № A1. P. 97–107. 1990. https://doi.org/10.1029/JA095iA01p00097

  11. Lanzerotti L.J., Hunsucker R.D., Rice D. et al. Ionosphere and ground-based response to field-aligned currents near the magnetospheric cusp regions // J. Geophys. Res. V. 92. № A7. P. 7739–7743. 1987. https://doi.org/10.1029/JA092iA07p07739

  12. Luhr H., Blawert A. Ground signatures of traveling convection vortices // Geophys. Monogr. Ser. V. 81. P. 231. 1994. https://doi.org/10.1029/GM081p0231

  13. McHenry M.A., Clauer C.R. Modeled ground magnetic signatures of flux transfer events // J. Geophys. Res. V. 92. № A10. P. 11231–11240. 1987. https://doi.org/10.1029/JA092iA10p11231

  14. Ngwira C.M., Pulkkinen A.A., Bernabeu E., Eichner J., Viljanen A., Crowley G. Characteristics of extreme geoelectric fields and their possible causes: Localized peak enhancements // Geophys. Res. Lett. V. 42. № 17. P. 6916–6921. 2015. https://doi.org/10.1002/2015GL065061

  15. Sibeck D.G., Korotova G.I. Testing models for traveling convection vortices: Two case studies // Geophys. Res. Lett. V. 27. № 3. P. 325–328. 2000. https://doi.org/10.1029/1999GL003700

Дополнительные материалы отсутствуют.