Известия РАН. Серия физическая, 2021, T. 85, № 3, стр. 366-371

Получение информации об ионосферно-магнитосферной плазме по наблюдениям полярных сияний

Б. В. Козелов 1*, А. В. Ролдугин 1

1 Федеральное государственное бюджетное научное учреждение Полярный геофизический институт
Апатиты, Россия

* E-mail: boris.kozelov@gmail.com

Поступила в редакцию 25.09.2020
После доработки 20.10.2020
Принята к публикации 27.11.2020

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

Аннотация

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

ВВЕДЕНИЕ

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

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

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

ОПРЕДЕЛЕНИЕ ВЫСОТЫ АВРОРАЛЬНОГО СВЕЧЕНИЯ ТРИАНГУЛЯЦИЕЙ

Первое большое исследование полярных сияний методом триангуляции было проведено еще Штормером в 1913 г. [1]. Было зарегистрировано более 300 авроральных структур, для которых получены высотные профили свечения. Как ни странно, в дальнейшем такой подход использовался редко, даже при переходе на более совершенные технологии регистрации изображений.

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

Необходимо отметить, что в системе ALIS [2] и некоторых других работах, например [3], реализован более сложный вариант триангуляционного подхода – пространственная томография области аврорального свечения при наблюдении из нескольких точек. Однако использование томографических алгоритмов восстановления накладывает дополнительные ограничения на качество изображений.

При использовании камер с небольшим полем зрения с небольшим пространственным разнесением и наблюдениях вблизи зенита можно избежать проблем идентификации структур: на камерах будет наблюдаться практически одна и та же структура, но смещенная на фоне удаленных звезд. По величине смещения можно определить высоту наблюдаемой структуры. Такой подход используется в системе авроральных камер MAIN (Multiscale Auroral Imaging Network), данные которой обсуждаются далее и описание которой приведено в работе [4]. Две камеры Guppy-1 и Guppy-2 в этой системе идентичны, имеют поле зрения по диагонали 18 градусов при угловом разрешении 0.038 градусов/пиксель и расположены на расстоянии 4.12 км друг от друга в направлении восток-запад. Камеры ориентированы на область вблизи магнитного зенита, и их работа синхронизована по сигналам GPS. Камеры снабжены одинаковыми стеклянными фильтрами, подавляющими красную часть спектра. Полоса пропускания и калибровка регистрируемых скоростей счета камер к абсолютным значениям приведены в работе [5]. Для анализа общей динамики полярных сияний будем использовать данные панхроматической (в широком спектральном диапазоне, охватывающем видимый свет) камеры всего неба, также входящей в систему MAIN.

На рис. 1 приведен пример регистрации тремя камерами системы MAIN свечения неба 30 марта 2017 с 23:00 UT до 23:30 UT. На верхней панели (рис. 1а) приведена кеограмма, построенная по данным камеры всего неба, т.е. зависимость от времени интенсивности свечения в проходящем через зенит сечении поля зрения камеры с севера на юг (здесь и далее изображение полярных сияний инвертировано, т.е. черное – интенсивность больше, белое – меньше). Из кеограммы видны основные морфологические особенности полярных сияний в рассматриваемый интервал времени: в течение всего интервала южнее камеры наблюдались пульсирующие сияния, в 23:02 UT на полюсной (северной) границе пульсирующих сияний появилась яркая дуга, дальнейшее развитие и расширение которой к северу привело после 23:10 UT к появлению активных ярких сияний в зените камеры. После 23:15 UT активизация сияний угасает, а зона пульсирующих сияний к 23:27 UT расширяется на все поле зрения.

Рис. 1.

Свечение неба, зарегистрированное камерами системы MAIN 30 марта 2017 с 23:00 UT до 23:30 UT: кеограмма для разреза в направлении север–юг, построенная по данным камеры всего неба (а); кеограммы для разреза вдоль эпиполярной линии восток–запад, построенные по данным камер Guppy-1 и Guppy-2 (б, в); зависимость эффективной высоты авроральных структур в поле зрения камер Guppy-1 и Guppy-2 (г), ось справа – энергия электронов, высотный профиль энерговыделения для которых имеет максимум на соответствующей высоте. Зенитные углы на панелях а–в в градусах.

Ниже на рис. 1б, 1в приведены кеограммы, построенные по данным разнесенных на 4.12 км камер Guppy-1 и Guppy-2, причем кеограммы строились по сечениям вдоль эпиполярной линии [6], лежащей в плоскости, проходящей через камеры и магнитный зенит. Постоянный фон, засветка от города, удален. Шкала абсолютных значений приведена для эмиссии 558 нм согласно данным калибровки [5] и соображению, что эта эмиссия составляет 35% интегральной интенсивности в сине-зеленой части спектра сияний. Видно, что структуры в сияниях повторяются на обеих кеограммах с некоторым пространственным смещением. По простым геометрическим соображениям, из треугольника с вершинами в наблюдаемой точке и в точках установки камер, это относительное угловое смещение может быть пересчитано в высоту структуры по формуле:

(1)
$h = \frac{{d~\cos \left( {w + f~dy} \right)}}{{\sin \left( {f~dx} \right)~}}.$
Здесь: w = 15° – зенитный угол центра изображения, dx – относительное смещение (параллакс) вдоль эпиполярной линии, в пикселах, f = 0.038 градусов/пиксел – угловое разрешение камер, dy – положение эпиполярной линии относительно центра изображения, пиксел, d = 4.12 км – расстояние между точками наблюдения.

На рис. 1г приведена зависимость “эффективной” высоты структур в поле зрения камер Guppy-1 и Guppy-2 в рассматриваемый интервал времени. “Эффективная” высота определялась по угловому смещению, которому соответствует максимум коэффициента корреляции по общему полю зрения для пары одновременных изображений. Так как контрастность структуры сияния на исходных изображениях не всегда достаточная, то предварительно изображения обрабатывались фильтром Собела.

ПОВЫШЕНИЕ КОНТРАСТА С ПОМОЩЬЮ ФИЛЬТРА СОБЕЛА

Фильтрация Собела является одним из типовых методов цифровой обработки изображений, и их реализация включена в основные математические пакеты [6, 7]. Также часто используемые фильтры Робертса и Превитта, в данной работе не рассматривались, так как для целей данной работы имеют недостатки по сравнению с фильтром Собела: фильтр Робертса более чувствителен к шуму, а фильтр Превитта дает большее сглаживание. Фактически, фильтр Собела – это оператор дискретного дифференцирования, приближенно вычисляющий градиент интенсивности изображения. Компоненты вектора градиента в каждой точке вычисляются путем свертки (сonvolution) с масками

(2)
${{X}_{{mask}}} = \left[ {\begin{array}{*{20}{c}} { - 1}&0&1 \\ { - 2}&0&2 \\ { - 1}&0&1 \end{array}} \right]~,\,\,\,\,{{Y}_{{mask}}} = \left[ {\begin{array}{*{20}{c}} 1&2&1 \\ 0&0&0 \\ { - 1}&{ - 2}&{ - 1} \end{array}} \right].$

Практический интерес представляет градиент на некотором характерном масштабе m пикселов, а не на минимальном масштабе при максимальном разрешении, при котором, на масштабе 1–2 пиксела, сказывается влияние звезд и случайных флуктуаций. Поэтому матрицы (2) масштабировались в m раз по обеим размерностям, при этом каждый элемент матрицы заменялся матрицей m × m, заполненной значением этого элемента. В данной работе использовалось значение = 8, т.е. масштаб по угловому разрешению 0.3°, что для типичных высот сияний (90–200 км) соответствует пространственным масштабам 0.5–1.0 км.

По полученным компонентам вектора градиента Gx и Gy можно получить абсолютную величину и направление градиента:

(3)
$A = \sqrt {G_{x}^{2} + G_{y}^{2}} ,\,\,\,\,~\varphi = {\text{arctan}}\left( {\frac{{{{G}_{y}}}}{{{{G}_{x}}}}} \right).$

На рис. 2 приведен пример обработки пары изображений с камер Guppy-1 и Guppy-2 с помощью описанной выше фильтрации. Компоненты вектора градиента Gx и Gy приведены красным и голубым цветом. На исходных изображениях присутствует диффузное свечение без очевидной структуры, однако после обработки видны волокнистые элементы, повторяющиеся на обоих снимках, что дает возможность проводить корреляцию и определить относительное смещение. Максимум корреляции для этой пары достигается при смещении 65 пикселей, что из формулы (1) соответствует высоте 92.3 км.

Рис. 2.

Пример фильтрации Собела для выделения структуры. Верхний ряд – исходные изображения, полученные камерами 31 марта 2017 в 00:33:30 UT. Нижний ряд – компоненты вектора градиента Gx и Gy приведены красным и голубым цветом.

ДИНАМИКА ВЫСОТЫ СИЯНИЙ И ЭНЕРГИИ ЭЛЕКТРОНОВ

Вернемся к обсуждению полученной высоты сияния, динамика которой в рассматриваемый интервал времени приведена на рис. 1г. Для сравнения на правой вертикальной оси рис. 1г отложена энергия моноэнергетического потока высыпающихся электронов с изотропным питч-угловым распределением, высотный профиль энерговыделения которого имеет максимум на соответствующей высоте [8]. Из рис. 1г видно, что первые 6 минут алгоритм триангуляции не дает устойчивого значения высоты, что вызвано помехами от облачности в поле зрения камеры Guppy-2. После 23:06 UT высота определяется уверенно, со сравнительно небольшим разбросом в соседние моменты времени. Исключениями являются небольшие интервалы вблизи 23:10 UT и 23:17 UT, когда полярные сияний “убегают” из поля зрения камер. Из сравнения с кеограммой, приведенной на рис. 1а и полученной по камере всего неба, можно выделить следующие особенности:

23:07–23:09 UT – появление параллельной слабой дуги немного к полюсу от яркой основной дуги. Высота этой дуги не ниже 95 км, что соответствует энергии электронов не более 30 кэВ.

23:10:30–23:16 UT – активизация дуги, прохождение небольшой омега-структуры в сияниях. Явно видно постепенное повышение эффективной высоты от 90 до 100 км, что может свидетельствовать о дисперсии по энергии инжектированных электронов – уменьшение энергии от 50 до 20 кэВ.

23:18-23:19 UT – пересечение поля зрения несколькими слабыми дугами. Высота свечения 102–104 км, энергия электронов 10–12 кэВ.

23:24–23:26 UT – волокна “черных сияний” (black aurora), высота 90 км, соответствующая энергия 50 кэВ.

В 23:26 UT и далее поле зрения камер заполнено пульсирующими сияниями, высота 92–100 км, энергия 20–40 кэВ, с крупномасштабной (~1 мин) модуляцией. Интересно, что увеличениям интенсивности свечения соответствует увеличение эффективной высоты, т.е. уменьшение средней энергии электронов.

ОЦЕНКА РАЗМЕРНОСТИ ПРОСТРАНСТВЕННО-ВРЕМЕННОЙ ДИНАМИКИ СИЯНИЙ

Визуально пространственно-временная динамика полярных сияний весьма сложна, может охватывать большой диапазон пространственных и временных масштабов [9], большой динамический диапазон. Тем не менее, во многих случаях эта динамика имеет внутреннюю организацию, не столь очевидную, но которая дает надежду на возможность ее простого описания, как систему сравнительно небольшой размерности [10, 11].

Покажем это для нашего события методом вложения [12], используя алгоритм Грассбергера–Прокаччи [13]. Ранее этот подход нами применялся при исследовании динамики полярных сияний по данным телевизионных наблюдений [10, 11]. Численный метод заключается в следующем. Из экспериментальных данных, измеренных во время исследуемого процесса, составляются эмпирические траектории системы в модельных фазовых пространствах разных размерностей d. Для траектории в каждом пространстве вычисляется корреляционный интеграл – вероятность Cd(ε) того, что расстояние между точками на траектории не превышает некоторой величины ε, и анализируется зависимость Cd(ε) от размерности вложения d и масштаба ε. Если динамика системы может быть описана небольшим числом дифференциальных уравнений, то с ростом d, существует такое de, что при d >de зависимость Cd(ε) практически не меняется, и на этой зависимости есть степенной участок Cd(ε) ~ εD.

Если исходный набор данных характеризует систему в (квази-) стационарном состоянии и d ≥ 2D + 1, то D является оценкой размерности аттрактора системы в ее фазовом пространстве. Если полученное значение D не является целым числом, то система имеет аттрактор, обладающий фрактальными свойствами. В нашем случае мы будем использовать такую оценку, как оценку размерности траектории системы, характеризующей степень сложности переходного процесса.

Для применения описанного численного алгоритма к набору данных в виде растровых изображений {I: IRN × RM} определим “расстояние” между изображениями, т.е определим метрическое пространство {I, dist}:

(4)
$\begin{gathered} dist~\left( {{\mathbf{I}}\left( {{{t}_{1}}} \right),{\mathbf{I}}\left( {{{t}_{2}}} \right)} \right) = \left| {{\mathbf{I}}\left( {{{t}_{1}}} \right) - {\mathbf{I}}\left( {{{t}_{2}}} \right)~} \right| = \\ = \sum\limits_{i = 1}^N {\sum\limits_{j = 1}^M {\left| {{{I}_{{i,j}}}\left( {{{t}_{1}}} \right) - {{I}_{{i,j}}}\left( {{{t}_{2}}} \right)} \right| \cdot {{2}^{{ - i - j}}}} } . \\ \end{gathered} $

В нашем случае растровые изображения имеют вид матриц 150 × 150 из центральной части поля зрения камеры всего неба, индексированных временем t с разрешением 1 секунда. Каждая такая матрица является “точкой” в метрическом пространстве с метрикой (4). Из последовательных d точек с постоянным шагом dt формируем “вектора” в пространстве вложения Id. Рассматривались d = 3, …, 10. Обоснование выбора шага dt опускаем для упрощения, для наших данных dt = 1 c.

Метрика пространствах вложения определяется аналогично:

(5)
$dist\left( {{{{\mathbf{I}}}^{d}}\left( {{{t}_{1}}} \right),{{{\mathbf{I}}}^{d}}\left( {{{t}_{2}}} \right)} \right) = \sum\limits_{k = 1}^d {dist({{I}_{k}}\left( {{{t}_{1}}} \right),{{I}_{k}}\left( {{{t}_{2}}~} \right))\cdot~{{2}^{{ - k}}}} .$

Тогда корреляционный интеграл определяется так:

(6)
${{C}_{d}}\left( {{\varepsilon ) = }\frac{1}{{{{N}^{2}}}}\sum\limits_{i = 1}^N {\sum\limits_{j = 1,j \ne i}^N {\theta (\varepsilon - dist({{{\mathbf{I}}}^{d}}\left( {{{t}_{i}}} \right),{{{\mathbf{I}}}^{d}}\left( {{{t}_{j}}} \right)} } } \right).$
Здесь θ – функция Хевисайда, N – число векторов в наборе данных.

Для анализа мы выбрали 15-минутные наборы последовательных изображений (900 шт. в наборе). Из них формировались “вектора” размерности d = 3 ,…, 10 и вычислялись корреляционные интегралы по (6). В нашем случае уже при d > 5 заметного изменения формы зависимости C(ε) не было, поэтому далее приводим результаты только для d = 6. На рис. 3а, 3б приведены результаты расчета корреляционного интеграла для интервалов 23:05–23:20 UT и 23:20–23:35 UT. На зависимостях C(ε) пунктиром выделены степенные участки ~εDc.

Рис. 3.

Результаты расчета корреляционного интеграла (6) при d = 6 для интервалов: 23:05–23:20 UT (а); 23:20–23:35 UT (б). Пунктиром выделены степенные участки ~εDc и приведены аппроксимации значений Dc.

Первый интервал, см. кеограмму на рис. 1, включает активизацию авроральной дуги на приполюсной границе пульсирующих сияний, прохождение небольшой омега-формы и расширению области активных сияний к полюсу. Несмотря на морфологически сложное строение структуры сияний, корреляционная размерность траектории оказывается небольшой Dc = 1.7, т.е. пространственно-временные изменения подчинены единой низко-размерной динамике. Близкие фрактальные значение корреляционной размерности характерны для траекторий систем вблизи фазовых переходов [14].

Во время второго анализируемого интервала, 23:20–23:35 UT, происходит постепенное заполнение всего поля зрения камеры пульсирующими сияниями: в начале интервала пульсирующие формы видны только в южной части, к концу интервала они занимают все поле зрения. Значение корреляционной размерности повышается до Dc = 2.8. Т.е. динамика становится более сложной, в ее формировании играет роль большее число степеней свободы. Следует отметить, что далее в интервале 23:35–23:50 UT, не приведенном здесь, когда пульсации практически неизменно сохраняются по всему небу, зависимость C(ε) не имеет явного степенного участка, а резко спадает к малым масштабам, что говорит о дальнейшем усложнении динамики сияний.

Такое изменение может быть объяснено следующим образом. Активизация сияний была связана с инжекцией энергичных частиц в магнитосфере [15], которая сформировала доминирующее крупномасштабное движение в плазме, что выразилось в малой размерности динамики сияний во время первого интервала. Далее, по мере ослабления возмущения, стал постепенно увеличиваться относительный вклад других локальных слабо коррелированных диссипационных процессов, таких как взаимодействие волн и частиц в отдельных дактах, о чем свидетельствуют пульсирующие формы в это время. Ранее в работе [10] было показано, что хотя динамика одного пульсирующего пятна в группе может быть близка периодической осцилляции (Dc = 2.0), для группы пятен пульсации могут быть слабо связаны между собой, что приводит большому значению Dc.

ЗАКЛЮЧЕНИЕ

Как видно из проведенных выше примеров, из данных оптических наблюдений можно получать полезные численные оценки, характеризующие динамические диссипативные процессы в магнитосферно-ионосферной системе. Конечно, этими примерами возможности оптических наблюдений не исчерпываются, и для полноценного использования их необходимо включать в комплексный анализ вместе с другими видами наблюдений. В частности, рассмотренный случай относится к одному из первых интервалов магнитного сопряжения японского спутника ERG(ARASE) с наземной аппаратурой ПГИ на Кольском полуострове. В работе [16] проведен одновременный анализ данных оптических наземных наблюдений камерой всего неба и данных регистрации на спутнике ERG электромагнитных волн ОНЧ диапазона, комплексный анализ которых будет продолжен в последующих работах.

Работа А.В. Ролдугина поддержана РФФИ (проект № 19-52-50025-ЯФ-а).

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

  1. Störmer C. // Kristiania, Geofysisk Publikationer 1. 1921. No. 5. P. 269.

  2. Steen Å., Brändström U. // STEP International Newsletter. 1993. V. 3. No. 5. P. 11.

  3. Frey S., Frey H.U., Carr D.J. et al. // J. Geophys. Res. 1996. V. 101. No. A10. Art. No. 21731.

  4. Kozelov B.V., Pilgaev S.V., Borovkov L.P., Yurov V.E. // Geosci. Instrum. Method. Data Syst. 2012. V. 1. P. 1.

  5. Kozelov B.V., Brändström B.U.E., Sigernes F. et al. // Proc. 36-th Ann. Seminar “Physics of Auroral Phenomena” (Apatity, 2013). P. 151.

  6. Визильтер Ю.В., Желтов С.Ю., Бондаренко А.В. и др. Обработка и анализ изображений в задачах машинного зрения: курс лекций и практических занятий. М.: Физматкнига, 2010. 672 с.

  7. Duda R., Hart P. Pattern classification and scene analysis. N.Y.: John Wiley and Sons, 1973. P. 271.

  8. Иванов И.Е., Козелов Б.В. Перенос электронных и протонно-водородных потоков в атмосфере Земли. Апатиты: Изд. КНЦ РАН, 2001. 260 с.

  9. Kozelov B.V., Uritsky V.M., Klimas A.J. // Geophys. Res. Lett. 2004. V. 31. Art. No. L20804.

  10. Kozelov B.V., Vjalkova N.Y. // Int. J. Geomagn. Aeron. 2005. V. 5. Art. No. GI3005.

  11. Kozelov B.V., Kozelova T.V., Kornilova T.A. // Proc. ICS-6 (Seattle, 2002). P. 432.

  12. Takens F. // Int. J. Bifurc. Chaos. 1993. V. 3. P. 241.

  13. Grassberger P., Procaccia I. // Phys. Rev. Lett. 1983. V. 50. P. 346.

  14. Sitnov M.I., Sharma A.S., Papadopoulos K. et al. // J. Geophys. Res. 2000. V. 105. Art. No. 12955.

  15. Лазутин Л.Л., Козелова Т.В., Мередит Н. и др.// Косм. исслед. 2007. Т. 45. № 2. С. 99; Lazutin L.L., Kozelova T.V., Meredith N.P. et al. // Cosm. Res. 2007. V. 45. No. 2. P. 89.

  16. Kawamura S., Hosokawa K., Kurita S. et al. // J. Geophys. Res. 2019. V. 124. No. 4. P. 2769.

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