Определение характеристик ИНЧ-волн, наиболее сильно реагирующих на незначительные изменения электронной плотности ионосферы в области высоких широт
Автор: Ахметов О.И., Мингалев И.В., Мингалев О.В., Суворова З.В., Белаховский В.Б., Черняков С.М.
Журнал: Солнечно-земная физика @solnechno-zemnaya-fizika
Статья в выпуске: 4 т.5, 2019 года.
Бесплатный доступ
Методами численного эксперимента исследуется задача определения характеристик электромагнитных волн УНЧ-диапазона (0.3-3 кГц), регистрируемых на уровне приземного слоя и несущих максимальное количество информации о состоянии волновода Земля-ионосфера. Проанализировано влияние горизонтальной пространственной структуры электронной плотности волновода Земля-ионосфера на особенности распространения электромагнитных волн. Выявлены характеристики, позволяющие регистрировать их инструментальными методами наблюдений в условиях слабо возмущенной ионосферы. Профили концентрации, используемые в численных экспериментах, получены по данным радара частичных отражений Полярного геофизического института, расположенного на радиофизическом полигоне «Туманный» Мурманской области (69.0° N, 35.7° E), с использованием модели IRI2016 во время солнечной вспышки 15 марта и последующей магнитной бури 17 марта 2013 г. Используемая в работе модель распространения электромагнитных сигналов разработана авторами...
Унч-волны, численное моделирование, ионосфера
Короткий адрес: https://sciup.org/142222503
IDR: 142222503 | DOI: 10.12737/szf-54201911
Текст научной статьи Определение характеристик ИНЧ-волн, наиболее сильно реагирующих на незначительные изменения электронной плотности ионосферы в области высоких широт
Одними из основных задач радиофизики в области распространения электромагнитных сигналов в атмосфере Земли являются поиск методов обеспечения стабильной связи или радиолокации на различных частотах; определение условий, когда эта связь невозможна; прогнозирование исключения таких условий и поиск мешающих связи или радиолокации факторов. Диагностика D-области ионосферы важна для многих физических приложений. В настоящее время основным отработанным средством такой диагностики являются радары частичных отражений. Ионозонды и радары некогерентного рассеяния формально захватывают D-область, но имеют там точность хуже, чем радары частичных отражений [Акимов, Калинин, 2017] . Состояние D-слоя ионосферы влияет на характеристики электромагнитных сигналов в диапазонах ИНЧ (0.3–3 кГц) и ОНЧ (3–30 кГц). Однако остается вопрос извлечения максимального количества информации о состоянии D-слоя из наземных измерений электромагнитных сигналов в указанном диапазоне.
В работе с помощью численного моделирования изучается влияние небольших (3–5 %) изменений концентрации электронов в D-слое ионосферы на характеристики электромагнитного сигнала на частоте 1500 Гц. Небольшие изменения концентрации электронов выбраны авторами с целью определения границ возможности диагностики D-слоя ионосферы сигналами ИНЧ-диапазона. Показано, что эти изменения не влияют на амплитуду отдельных компонент электромагнитного поля сигнала на уровне поверхности Земли, оказывают незначительное влияние на отношение амплитуд напряженностей электрического и магнитного полей и заметное влияние на разность фаз между напряженностями электрического и магнитного полей. Также исследовано влияние горизонтальных размеров области повышенной концентрации на перечисленные характеристики сигнала.
Для численного моделирования распространения сигналов ИНЧ-диапазона авторы использовали ранее разработанную ими модель распространения электромагнитных сигналов в системе литосфера— атмосфера—ионосфера. В модели численно интегрируется по времени система уравнений Максвелла на регулярной пространственной сетке. Для этого используется явная схема расщепления по пространственным направлениям и физическим процессам с противопотоковой аппроксимацией пространственных производных (метод Годунова с коррекцией потоков). Эта схема является консервативной, монотонной, имеет второй порядок точности по времени и третий по пространственным переменным.
В представленных численных экспериментах авторы используют прямоугольную область моделирования и декартову систему координат, пренебрегая кривизной сферических слоев литосферы, атмосферы и ионосферы по следующим причинам. Наибольший из горизонтальных размеров области моделирования примерно в шесть раз меньше радиуса Земли, а длина волны (~200 км) сравнима с размерами области моделирования. При таком соотношении масштабов сигнал распространяется горизонтально без многократных отражений от ионосферы и литосферы. Целью моделирования является исследование влияния локальных неоднородностей концентрации электронов в D-области ионосферы на характеристики электромагнитного сигнала. Плоская геометрия области моделирования для всех представленных численных экспериментов будет давать одинаковую малую ошибку, которая не окажет влияния на обсуждаемые закономерности в характеристиках сигнала.
ПОСТАНОВКА ЗАДАЧИ И ЗАДАНИЕ СРЕДЫ
В представленной работе авторы подробно рассматривают влияние незначительного повышения электронной плотности в нижней ионосфере в области высоких широт на прохождение электромагнитных волн ИНЧ-диапазона при различной конфигурации горизонтальных профилей электронной концентрации. Параметры среды задавались следующим образом: концентрация нейтральных частиц — с помощью эмпирической модели NRLMSISE-00; концентрация и температура электронов — с помощью эмпирической модели IRI 2016 с коррекцией по записям средневолнового радиолокатора вертикального излучения для исследования нижней ионосферы, расположенного на радиофизическом полигоне «Туманный» Полярного геофизического института в Мурманской области (69.0° N, 35.7° E) [Терещенко и др., 2003] , во время и после солнечной вспышки 15.03.2013, а именно 15.03.2013 в 09:00 UT (возмущенные условия, вспышка), 16.03.2013 в 09:00 UT (спокойные условия), 17.03.2013 в 06:00 UT (возмущенные условия, магнитная буря) и 16.03.2013 в 06:00 UT (спокойные условия). Вертикальные профили электронной плотности в указанных условиях показаны на рис. 1. Коррекция заключалась в перенесении тонкой структуры D- и E-областей, полученной по данным радиолокатора, на профили, полученные с помощью эмпирических моделей IRI 2016 и NRLMSISE.
Магнитная буря 17 марта 2013 г. была вызвана солнечной вспышкой М-класса 15 марта и началась с приходом к магнитосфере межпланетной ударной волны, которая проявилась на земной поверхности

Рис. 1. Аппроксимированные высотные профили концентрации электронов 15.03.2013 в 09:00 UT ( а ); 16.03.2013 в 09:00 UT ( б ); 17.03.2013 в 06:00 UT ( в ) и 16.03.2013 в 06:00 UT ( г )

Рис. 2. Горизонтальные профили концентрации электронов для эксперимента по распространению электромагнитной волны из спокойной области в возмущенную ( а ); из возмущенной области в спокойную ( б ); сквозь полосу повышенной концентрации электронов ( в )
в виде SSC (sudden storm commencement) в 06:00 UT. В этот момент скорость солнечного ветра резко увеличилась с ~400 до ~650–700 км/с, B z -компонента ММП стала отрицательной, что обеспечило постоянное поступление энергии в магнитосферу. Индекс SYM-H , характеризующий интенсивность магнитной бури, упал до –100 нТл и оставался на этом уровне. Авроральный AE -индекс резко возрос до ~1000 нТл и оставался повышенным. АЕ -индекс показывал в 16:00 UT появление другой авроральной активизации (увеличение AE до ~2500 нТл).
Несмотря на высокий класс вспышки, радар частичных отражений на полигоне «Туманный» и модель IRI показали незначительное повышение концентрации электронов на высотах D–E-областей ионосферы во время последующей магнитной бури. Авторы полагают, что такие условия хорошо подходят для разделения технически наблюдаемых и ненаблюдаемых характеристик электромагнитных волн и выделения тех из них, которые могут быть использованы для создания методов косвенных исследований ионосферы.
Профиль проводимости литосферы был взят из работы [Korja et al., 2002] , проводимость у поверхности составляла 2∙10–5 См.
Пространственная горизонтальная структура электронной плотности тоже может оказывать влияние на получаемые результаты. Для обнаружения этого влияния были проведены три серии экспериментов с горизонтальными профилями электронной плотности разных типов. Первая серия моделирует переход из спокойной области в область повышенной концентрации электронов. Горизонтальный профиль электронной концентрации показан на рис. 2, а. Вторая серия моделирует переход из области повышенной концентрации в спокойную область (профиль на рис. 2, б). Третья серия — прохождение полосы повышенной концентрации (рис. 2, в). Для всех серий переход осуществлялся по функции Гаусса. При этом в первых двух изменялся горизонтальный градиент для обнаружения его влияния в области перехода на распространение электромагнитных волн, в третьей изменялась ширина полосы с 0.25λ до λ, где λ — длина волны. Магнитное поле было взято вертикальным, что характерно для высоких широт, его индукция составляла 5.3^10-5 нТл.
ОБЛАСТЬ МОДЕЛИРОВАНИЯ, ПРОСТРАНСТВЕННАЯ СЕТКА И ИСТОЧНИК СИГНАЛА
Во всех поставленных в работе численных экспериментах область моделирования представляла собой параллелепипед с основанием 512x1280 км, высотой в атмосфере 150 км и глубиной в литосфере 50 км. Шаги сетки по горизонтали составляли 2 км, по вертикали в атмосфере 2 км и 0.5 км в литосфере. Шаг интегрирования по времени составлял 4^10-6 с.
На всех границах области моделирования применялись специальные условия свободного ухода волны для минимизации влияния нефизических отражений на результаты моделирования. Авторы применяли для моделирования постоянно действующий гармонический источник сигнала, а результаты анализировались при прекращении всех процессов, связанных с прохождением первичного фронта сигнала (35-й период колебаний). Поэтому, несмотря на почти 100 % отражение сигнала от высокопроводящих литосферы и ионосферы, находящихся внутри области моделирования, для исключения любой возможности накопления ошибок на вертикальных гранях расчетной области применялись такие же условия, как и на горизонтальных.
В качестве источника сигнала во всех представленных экспериментах используется поверхность одной из граней расчетной области. Такой источник сигнала позволяет задать не только амплитуду сигнала во времени, но и ее распределение в пространстве; задержками можно сформировать фронт волны необходимой формы подобно тому, как это делается в плоской эквидистантной фазированной антенной решетке [Воскресенский и др., 2012] . В экспериментах, обсуждаемых в данной работе, моделировалась плоская волна частотой 1500 Гц (длина волны ~200 км), излучаемая в область под прямым углом к плоскости источника.
ОПИСАНИЕ МОДЕЛИ
Используемая в представленной работе модель распространения электромагнитных сигналов в волноводе Земля—ионосфера с учетом инерции электронов в ионосфере построена на численном интегрировании по времени системы уравнений Максвелла на регулярной пространственной сетке с помощью явной схемы. Используется расщепление по пространственным направлениям и по физическим процессам. В схеме чередуются шаги распространения по каждому направлению и шаги затухания и вращения электрического поля сигнала. Также на шагах распространения по пространственным направлениям используется противопотоковая аппроксимация пространственных производных (метод Годунова с коррекцией потоков). Схема является консервативной, монотонной, имеет второй порядок точности по времени и третий по пространственным переменным.
Схема численного интегрирования уравнений Максвелла в литосфере
Пусть r =( x , y , z ) — декартовы координаты; t — время; E ( r , t ), D ( r , t ), H ( r , t ) и B ( r , t ) — напряженность и индукция электрического и магнитного полей; j ( r , t ) — плотность тока в момент времени t в точке с радиус-вектором г . Рассмотрим уравнение Фарадея и уравнение Максвелла в системе СИ
^ B = - rot E ( r , t ), ^ D = rot H ( r , t ) - j ( r , t ). (1)
о t d t
В литосфере эти уравнения замыкаются законом Ома в форме
j ( r , t ) = о ( г ) Е ( г , t ),
где o( r ) — скалярная проводимость среды, и материальными уравнениями
D ( r , t ) = S o S ( r ) E ( r , t ), B ( r , t ) = P 0 P ( r ) H ( r , t ), (3)
где e(r) и p(r) — безразмерные относительные ди- электрическая и магнитная проницаемости среды в низкочастотном пределе; е0 и р0 — электрическая и магнитная проницаемости вакуума.
В представленных расчетах мы полагали, что в литосфере p( r )=1, однако мы сохраним р в дальнейших формулах, поскольку используемая методика позволяет рассматривать общий случай с переменной p( r ).
Пусть c 0 = 1/ ^ £ 0 p 0 — скорость света в вакууме;
c(r) = c0 / V£(r)p(r) — скорость света в среде в точке с радиус-вектором r. Систему (1)-(3) можно записать в виде ав dE c2 +
= -rotE, = -rot д t д t £
B
Ц
- —E.
£ 0 £
Для численного интегрирования системы (4) применим метод расщепления по физическим процессам. Общий шаг интегрирования разделяется на два подшага. Один из них — это подшаг распространения, на котором интегрируется система уравнений dB д t
dE с22
=-rotE, = —rot II dt £(pj
второй — подшаг затухания сигнала, на котором аналитически интегрируется система уравнений dE _
--=--E dt£ по формулам E(r, t + t) = E(r, t) exp I - I. Пра-
( £0£(r)J вильное чередование подшагов расщепления обеспечивает второй порядок точности по времени.
Рассмотрим численное интегрирование системы
(5). Введем перенормированные поля
B(r, t) = c0B(r, t), E(r, t) = c0 / c(r)E(r, t), а также вектор M(r) = Vc в представленных расчетах, равный M(r) = Vc + cVц(r)/ ц(г) в общем случае. Умножая уравнение Фарадея на c0, получим уравнение
-
8 В
— = - rot(c E ). (6)
-
8 t
Подстановка введенных обозначений в уравнение Максвелла приводит это уравнение к виду dE = rot ( c — )-[m X — ].
Введем матрицы Rx, Ry, Rz и нулевую мат- рицу O:
" 0 0 |
0 ) |
f 0 0 1 ) |
||
R + = |
0 0 |
- 1 |
, R У = |
0 0 0 |
+ 0 1 |
0 у |
+- 1 0 0 , |
||
' 0 - 1 |
0' |
f 0 0 0 ) |
||
R z = |
1 0 |
0 |
, O = |
0 0 0 . |
+ 0 0 |
0 |
+ 0 0 0 v |
Оператор rot в декартовых координатах можно представить в виде
-
da da
rota = Rx — + Ry — + Rz —.
8x dy
Введем шестимерные вектор-столбец переменных u = ( B x , B y , B z , E x , E y , E z ) T и вектор F = ( F , F 2, F 3, F 4, F 5, F 6 ) T , компоненты которого заданы формулами
F j = F 2 = F 3 = 0, ( F 4 , F 5 , F 6 ) T = -[ m x — ] . (10)
Используя обозначения (8), а также формулы (9) и (10), систему уравнений (6), (7) можно представить в виде векторного уравнения du d(cu) d(cu) d( cu)
+ A x + A y + A z = F , (11)
8t 8x dy 8z где Ax, Ay, Az — симметричные матрицы размером 6X6, заданные формулами
^
A x
O R x
(- R x O )
^

О R y
(- R y O J
^
A z
'•> ^,
O R z
(- R z O J
Векторное уравнение (11) задает шестимерную линейную гиперболическую систему уравнений первого порядка, записанную в консервативной форме. Ее правая часть F линейно зависит от компонент вектора u. Для численного интегрирования таких систем разработано достаточно много разностных схем, в том числе схемы повышенного порядка точности, которые применяются для уравнений газовой динамики. Наиболее полное описание этих схем содержится в монографиях [Куликовский и др., 2012; Бисикало и др., 2013].
Используя метод расщепления по пространственным направлениям, можно построить явные монотонные схемы численного интегрирования системы (11), которые сводятся к последовательному интегрированию одномерных по пространству гиперболических систем уравнений. На каждом временном шаге нужно последовательно проинтегрировать три системы уравнений, например в таком порядке:
d u ' d ( c u )
+ A x— 8 t 8 x
= F , ^ A , ^ = f , x 8 t 8 y y
d um
---+
8 t
A z
d ( c "') = F d z z
При этом правые части систем (13) выбираются так, чтобы они не изменялись на своем шаге расщепления и удовлетворяли равенству
F = F x + F y + F z .
Эти условия будут выполнены, если задать их по формулам
T
0,0,0,0, - MB. , mb \ , zx yx
T
0,0,0, - M z B y ,0, - M x B y ) ,
T
0,0,0, - M y B z , M x B z ,0 ) .
На каждом шаге расщепления для двух компонент магнитного поля и двух компонент электрического поля, ортогональных направлению шага, рассчитывается распространение сигнала конечно-разностным способом. В качестве начальных условий для каждой системы уравнений в (13) берутся значения, рассчитанные на предыдущем шаге расщепления. Сохранить второй порядок аппроксимации по времени в схеме расщепления можно, если циклически изменять порядок выполнения шагов расщепления, например, выполняя сначала шаги по пространственным направлениям в следующем порядке: xyz , yxz , zxy , xzy , yzx , zyx . Обоснование этого утверждения содержится, например, в монографиях [Бисикало и др., 2013; Ковеня и др., 1981] .
Пусть заданы равномерная сетка по времени tn=10+nт, где т — шаг по времени, и равномерная пространственная сетка в декартовых координатах, целые и полуцелые узлы которых заданы соотношениями x = x0 + hxx, y, = y 0 + jhy , zk = z 0 + khz, r, j, k =( x,, y,, zk), xi+1/2 = x0 + (i + 1Z 2) hx, y,+1/2 = y0 + (j + 1Z 2) hy , zk+1/2 = z 0 +( k + 1/2) hz , в которых hx, hy, hz — шаги сетки по осям X, Y, Z соответственно. Будем использовать для значений функции f в узлах сетки обозначения f (rj,k, tn) = f ”,,k.
Рассмотрим послойный переход для первой системы из (13), которая интегрируется на шаге расщепления по направлению х . В этом случае представим заданную первой формулой в (10) матрицу
Ax в виде Ax = QЛQ , где Л — диагональная матрица, на диагонали которой стоят собственные числа матрицы Ax; Q — матрица, столбцы которой есть правые собственные векторы матрицы Ax, опре
- 1
деленные с точностью до множителя; Q — матрица, обратная Q. Строки матрицы Q 1 — левые соб ственные векторы матрицы Ax. Эти матрицы постоянны. Их можно взять в виде
( 0 0 0 0 0 0 )
0 1 0 0 0 0
Л =
0 - 1 0 0
0 0 0 0
0 0 0 0 1 0
( 0 0
0 0
0 - 1 )
(1 0 0 0 00
0 1/2 0 00 1/2
0 0 1/2 0 1/20
0 0 0 1 00
0 0 -1/2 0 1/20
( 0 - 1/2
0 0
0 1/2 )
0 0 0
0 0 ^
0 1 0 0
0 - 1
-•---1
Q
0 0
0 0
0 0
. 0 1
1 0 - 1 0
0 1
0 0
1 0 1 0
0 0 0
С помощью введенных матриц Л , Q и Q 1 первую систему уравнений в (13) представим в виде векторного уравнения
1 д ( c u ) д x
du '"
+ Q Л Q д t
= F x .
Введем вектор-столбец характеристических переменных (линеаризованные инварианты Римана)
1 по направлению x , заданный формулой w = Q u .
Компоненты вектора w заданы формулами
—-—~—•■ -"^—^ '•—~—^ .--——^ --—^—■ w1 Bx , w2 = By - Ez , w3 = Bz - Ey , w4 = Ex, w5 = Bz + Ey , w6 = By + Ez .
- 1
Умножим уравнение (15) слева на матрицу Q .
В результате получим векторное уравнение
^ w + Л ^w = Q - 1 F , (16)
д t дx x которое является системой из шести скалярных уравнений для характеристических переменных.
В итоге численное интегрирование первой системы из (13) сводится к численному интегрированию четырех независимых одномерных по пространству уравнений переноса (16) для характеристических переменных w2 = By - Ez , w3 = Bz - Ey , w5 = Bz + Ey, w6 = By + Ez .
Это интегрирование выполняется с помощью явной противопотоковой схемы, имеющей третий порядок точности по пространству и второй — по времени. Эта схема подробно описана в работах [Мингалев и др., 2018а, б] .
Рассмотрим кратко основные идеи, лежащие в основе численной схемы.
Схема численного интегрирования уравнений Максвелла в ионосфере
Для ионосферы мы полагали, что безразмерная относительная магнитная проницаемость среды p ( r ) = 1 и выполняется формула B ( r , t ) = ц 0 H ( r , t ). Также мы полагали, что имеет место поляризация плазмы D ( r , t ) = s 0 E ( r , t ) + P ( r , t ), где P ( r , t ) — вектор поляризации, причем плотность тока поляризации д Р / д t = j совпадает с полной плотностью тока в плазме. Уравнение Фарадея и уравнение Максвелла в этом случае принимают вид
— = - rot E , — = c 0'rot B - — j . (17)
д t , д t 0 s 0
Система (17) замыкается уравнением для плотности тока электронов, вызванного полем сигнала:
fj =- V ej -^e [ jX b^ro^E, д t где ve — частота столкновений электронов, Qe = e |Ввн| / me — гирочастота электронов, Ввн — внешнее геомагнитное поле, b = B / B , вн вн
® 2 = e 2 n e / ( m e e 0 ) — квадрат плазменной частоты электронов.
В ионосфере, как и в литосфере, применяется метод расщепления по физическим процессам. Правильное чередование подшагов расщепления обеспечивает второй порядок точности по времени. На подшаге распространения учитывается только распространение сигнала и численно интегрируется система уравнений
^ B = - rot E , — = c 22 rot B (18)
д t д t 0
с помощью той же схемы, что и система (5) в литосфере.
На подшаге затухания и вращения сигнала магнитное поле сигнала не изменяется и учитываются только затухание электрического поля и его вращение за счет внешнего геомагнитного поля. При этом в каждой точке расчетной сетки аналитически интегрируется система уравнений
| E = - — j, = - V e j a [ j X b ] + ro 2 S o E .
д t ^ o д t
Эта система разделяется на две независимых системы для продольных электрического поля E = ( b , j ) и тока j\ = ( b , j )
d E 2 . d j _ 2
, = j\\ , — = - V e j\\ + ® e E 0 E || (20)
dt Eo dt и поперечных электрического поля E± — E - bЕц и тока L = j — bj|| :
dE± = _2 • dt Eo
± ,
^ j T = - V e j ± - Q e [ L X b ] + ® 2 e 0 E ± .
d t
Системы (20) и (21) являются автономными линейными системами однородных дифференциальных уравнений (ОДУ) с постоянными коэффициентами. Решения задачи Коши для этих систем выражаются достаточно громоздкими аналитическими формулами.
В той части ионосферы, где выполнено неравенство ve t> 20, системы (20) и (21) заменяются на более простые уравнения и систему, дЕ\ д t
—
2 ® 2F
—Е\\
v e
5 E ± _ 2 ® 2 v e д t v 2 + Q 2
E ± — ^^ [ E ± x b ] v
которые также имеют аналитические решения задачи Коши.
Изложенная схема позволяет учесть частотную дисперсию и моделировать распространение широкополосных сигналов произвольной формы.
Условия на границе области моделирования
В модельных экспериментах на всех внешних границах было задано граничное условие свободного ухода волны, что достигалось обнулением входящего из-за пределов области моделирования потока. Такой подход позволяет получить низкие коэффициенты отражения плоской волны от границ области моделирования: для волн, падающих под углом от 80° до 90°, отношение амплитуды отраженной волны к амплитуде падающей плоской волны не превышает 0.01. При угле падения 60° это отношение составляет уже ~0.05, при угле 45° — ~0.16, при угле 27° — ~0.33, а при угле 18.4° — ~0.43 [Мингалев и др., 2018а]. Для сравнения: метод FDTD (finite-differences time-domain method) [Yee, 1966] при применении простых граничных условий, таких как условия Мура [Mur, 1981] и Лиао [Liao, 1984], дает отражения порядка 0.1…1 %, но только при падении волны на границу под прямым углом. При падении под острым углом коэффициент отражения растет до 100 % при падении по касательной. Однако при использовании непрерывно действующего источника даже столь малых отражений, которые порождает применяемая схема, достаточно для накопления ошибок в области моделирования, поэтому возникает необходимость использования методов подавления, подобных PML (perfectly matched layer), использующихся в FDTD-моделях [Berenger, 1994]. Именно такой тип источника применялся в представленных авторами экспериментах, что привело к необходимости адаптации и применения метода PML. Разделение схемы по пространственным переменным и физическим процессам позволяет применять профиль электрических и магнитных потерь, предложенный Беренгером, непосредственно к потокам противопотоковой схемы на границе области моделирования. Геометрический профиль потерь внутри отдельного слоя имеет вид
Р ( r) = — ^° c 0 ln( g )n ln( R o ) g* r / " ^ ) , (22) 2 A x ( g N — 1)
где g — коэффициент геометрической прогрессии; Δ x — шаг по пространству; c 0 — скорость света; N — номер PML-слоя, считая от интерфейса счетного региона и границы; r — расстояние от границы; R 0 — коэффициент отражения от первого слоя. В представленных численных экспериментах авторы используют профиль потерь, рассчитанный по формуле (22) со следующими параметрами: R 0 =0.01 (1 %), коэффициент прогрессии g =2.15, количество слоев 14. Несмотря на то, что коэффициент отражения от первого слоя не меньше, чем характерный для данной схемы при обнулении исходящих потоков на углах падения 80°–90°, а на практике даже больше вследствие отражений от последующих слоев, основным преимуществом метода PML является его крайне слабая зависимость от угла прихода электромагнитной волны. Данную особенность демонстрирует и адаптированный для противопотоковой схемы вариант.
РЕЗУЛЬТАТЫ И ОБСУЖДЕНИЕ
Авторами были проведены численные эксперименты семи типов: четыре эксперимента — с модельной областью, однородной по горизонтали и вертикальными профилями электронной плотности, как на рис. 1; три серии экспериментов — с переходом в направлении распространения волны, представленные на рис. 2. Модельный сигнал, исследуемый для всех представленных в работе численных экспериментов, задавался синусоидальной функцией с частотой 1500 Гц. Авторы не обсуждают переходные процессы, связанные с началом счета и прохождением фронта волны через области с не нулевым горизонтальным градиентом концентрации электронов. Все обсуждаемые далее результаты получены для устоявшегося режима (30-й –35-й периоды колебаний от начала моделирования). В связи с этим анализировать непосредственно временные формы компонент электромагнитного сигнала нецелесообразно, так как они идентичны исходному сигналу. Далее авторы сравнивают и обсуждают отдельные характеристики электромагнитного сигнала, такие как амплитуда и фаза, а также их производные.

Рис. 3. Разность амплитуд в процентах для электрического ( а ) и магнитного полей ( б ) для случаев распространения сигнала, когда вся область имела профиль электронной концентрации, характерный для спокойных условий 16.03.2013 в 06:00 UT и когда горизонтальный профиль концентрации электронов имел переход из спокойной области в область с повышенной концентрацией электронов
На рис. 3, а изображена разность амплитуд в процентах основной компоненты электрического поля волны для случаев распространения сигнала, когда вся область имела профиль электронной концентрации, характерный для спокойных условий 16.03.2013 в 06:00 UT и когда горизонтальный профиль концентрации электронов имел переход из спокойной области в область с повышенной концентрацией электронов, которая имела профиль, характерный для условий магнитной бури 17.03.2013 в 06:00 UT. На рис. 3, б изображена разность амплитуд основной компоненты магнитного поля волны в тех же экспериментах. На рис. 4, а , б изображена разность амплитуд основной компоненты электрического и магнитного полей волны для случаев распространения сигнала, когда вся область имела профиль электронной концентрации, характерный для условий магнитной бури 17.03.2013 в 06:00 UT и когда горизонтальный профиль концентрации электронов имел переход из области с повышенной концентрацией электронов в спокойную область, которая имела профиль, характерный для 16.03.2013, 06:00 UT. Видно, что распространение через области как положительного, так и отрицательного горизонтального градиентов электронной плотности оказывает влияние на амплитудные характеристики электромагнитных волн тем сильнее, чем выше градиент в области перехода. Однако в случае слабых возмущений электронной плотности изменения амплитуд основных компонент электромагнитного поля не достаточны для регистрации методами наземных наблюдений.
Анализ прохождения электромагнитной волной полосы повышенной концентрации электронов, расположенной перпендикулярно направлению распространения и параллельно плоскости волны, не выявил значительного влияния на амплитудные характеристики сигнала, независимо от соотношения ширины полосы и длины волны. Очевидно, что отдельная регистрация компонент электромагнитного поля не может быть использована для анализа слабых возмущений концентрации электронов в ионосфере, так как амплитудный отклик на них не превышает значений фонового шума.
Ситуация несколько улучшается в случае совместного анализа электрических и магнитных компонент электромагнитного поля. На рис. 5 показаны графики волнового сопротивления среды, рассчитываемого как отношение напряженностей магнитного и электрического полей для случаев горизонтально однородной среды с вертикальными профилями всех типов, представленных на рис. 1, и для случаев прохождения полосы повышенной концентрации электронов шириной λ и 0.25λ. Видно, что такой подход позволяет достаточно четко отделять случаи с наличием областей возмущений электронной плотности в ионосфере на радиотрассе от случаев спокойных условий, но слабо реагирует как на величину возмущений электронной плотности, так и на горизонтальную структуру данных возмущений. Для реализации данного подхода при анализе данных регистрирующая аппаратура в случае плоской волны, распространяющейся в плоскости волновода Земля—ионосфера, должна позволять фиксировать

Рис. 4. Разность амплитуд в процентах для случаев распространения сигнала, когда вся область имела профиль электронной концентрации условий магнитной бури 17.03.2013 в 06:00 UT и когда горизонтальный профиль концентрации электронов имел переход из области с повышенной концентрацией электронов в спокойную область: а — напряженность электрического поля; б — напряженность магнитного поля

Рис. 5. Волновое сопротивление среды
вертикальную компоненту электрического поля и две горизонтальные компоненты магнитного поля. В случае волн другого типа или произвольного направления распространения необходима фиксация всех шести компонент электромагнитного поля.
На рис. 6 показаны разности фаз напряженностей магнитного и электрического полей для случаев горизонтально однородной среды с вертикальными профилями всех типов, представленными на рис. 1, и для случаев прохождения полосы повышенной концентрации электронов шириной λ и 0.25λ. Видно, что фазовые характеристики наиболее сильно дифференцированы в зависимости от величины возмущения электронной плотности как в случае изменения вертикального профиля концентрации электронов, так и в случае горизонтального профиля в сравнении с ранее рассмотренными. Очевидно, что остается неопределенность, какого типа возмущение выз-

Рис. 6. Разность фаз между напряженностями электрического и магнитного полей
вало изменение фазы в каждом конкретном случае. Также остается вопрос зависимости сдвига фазы от формы профиля электронной концентрации, а также ее величины для волн различных диапазонов. Эти и многие другие вопросы требуют дополнительной проработки, выходящей за рамки представленной работы.
Для проведения анализа фазовых характеристик необходимо при регистрации компонент электромагнитного поля осуществлять точную привязку по времени регистрируемых электрических и магнитных компонент. В случае создания сети приемников, что может помочь с идентификацией горизонтальных ионосферных неоднородностей, необходимо обеспечить точную привязку по времени регистрируемых компонент поля к мировому времени, как, например, это делается в работе [Пильгаев и др., 2008] .
ЗАКЛЮЧЕНИЕ
Методами численного эксперимента авторами выполнена задача определения характеристик электромагнитных волн ИНЧ-диапазона, регистрируемых на уровне приземного слоя, несущих максимальное количество информации о состоянии волновода Земля—ионосфера.
Выявлено влияние горизонтальных градиентов концентрации электронов на амплитуду электромагнитных сигналов ИНЧ-диапазона при переходе между областями с различными характеристиками волновода Земля—ионосфера. Показано, что в случае незначительных изменений электронной плотности изменения амплитуды электромагнитного сигнала незначительны и не превышают 1 %.
Выявлено существенное влияние вида вертикального профиля электронной концентрации на волновое сопротивление среды и разность фаз между магнитной и электрической компонентами. Показано, что в случаях увеличения концентрации электронов в D-области волновое сопротивление среды увеличивается, а компонента электромагнитного поля E начинает запаздывать относительно H . По мнению авторов, отставание фазы электрического поля от магнитного обусловлено инерцией электронов и их столкновениями с нейтральными частицами ионосферы. Вызванная электрическим полем сигнала плотность электрического тока влияет на сдвиг фаз между моментами достижения максимальных значений электрическим и магнитным полями сигнала в фиксированной точке ионосферы.
На основании проведенных расчетов можно сделать вывод, что минимально пригодные для мониторинга состояния ионосферы наземные системы регистрации ИНЧ-волн должны содержать в своем составе три регистратора: один для вертикальной электрической компоненты и два для горизонтальных магнитных компонент электромагнитного поля. При этом важно обеспечить качественную синхронизацию между регистраторами и мировым временем.
Исследование выполнено при поддержке гранта Российского научного фонда (проект № 18-77-10018).
Список литературы Определение характеристик ИНЧ-волн, наиболее сильно реагирующих на незначительные изменения электронной плотности ионосферы в области высоких широт
- 1. Акимов В.Ф., Калинин Ю.К. Введение в проектирование ионосферных загоризонтных радиолокаторов. М.: Техносфера, 2017. 492 с.
- 2. Бисикало Д.В., Жилкин А.Г., Боярчук А.А. Газодинамика тесных двойных звезд. М.: Физматлит, 2013. 632 с.
- 3. Воскресенский Д.И., Степаненко В.И., Филиппов В.С. Устройства СВЧ и антенны. Проектирование фазированных антенных решеток. 4-е изд. М.: Радиотехника, 2012. 744 с.
- 4. Ковеня В.М., Яненко Н.Н. Метод расщепления в задачах газовой динамики. Новосибирск: Наука, 1981. 304 с.
- 5. Куликовский А.Г., Погорелов Н.В., Семенов А.Ю. Математические вопросы численного решения гиперболических систем уравнений. 2-е изд. М.: Физматлит, 2012. 656 с.
- 6. Мингалев И.В., Мингалев О.В., Ахметов О.И., Суворова З.В. Явная схема расщепления для уравнений Максвелла // Математическое моделирование. М.: РАН, 2018а. Т. 30, № 12. С. 17-38. DOI: 10.31857/S023408790001934-1.
- 7. Мингалев О.В., Мингалев И.В., Мельник М.Н., Ахметов О.И. и др. Новый метод численного интегрирования системы Власова-Максвелла // Математическое моделирование. М.: РАН. 2018б. Т. 30, № 10. С. 21-43. DOI: 10.31857/ S023408790001919-4.
- 8. Пильгаев С.В., Ахметов О.И., Филатов М.В., Федоренко Ю.В. Универсальное устройство синхронизации данных от GPS-приемника // Приборы и техника эксперимента. М.: РАН, 2008. № 3. С. 175-176.
- 9. Терещенко В.Д., Васильев Е.Б., Овчинников Н.А., Попов А.А. Средневолновый радиолокатор Полярного геофизического института для исследования нижней ионосферы // Техника и методика геофизического эксперимента, Апатиты: Изд. Кольского научного центра РАН, 2003. С. 37-46.
- 10. Berenger J.-P. A perfectly matched layer for the absorption of electromagnetic waves // J. Computational Phys. 1994. V. 114, N 2. P. 185-200. DOI: 10.1006/jcph.1994.1159.
- 11. Korja T., Engels M., Zhamaletdinov A.A., et al. Crustal conductivity in Fennoscandia - a compilation of a database on crustal conductance in the fennoscandian shield // Earth Planets Space. 2002. V. 54, N 5. P. 535-558. DOI: 10.1186/BF03353044.
- 12. Liao Z.P., Wong H.L., Yang B.P., Yuan Y.F. A transmitting boundary for transient wave analyses // Scienta Sinica (series A). 1984. V. 27, N 10. P. 1063-1076.
- 13. Mur G. Absorbing boundary conditions for the finite-difference approximation of the time domain electromagnetic field equations // IEEE Electromagnetic Compatibility. 1981. V. 23, N 4. P. 277-382. DOI: 10.1109/TEMC.1981.303970.
- 14. Yee Kane. Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media // IEEE Transactions on Antennas and Propagation. 1966. V. 14. P. 302-307. DOI: 10.1109/TAP.1966.1138693.