Радиотехника и электроника, 2020, T. 65, № 5, стр. 434-441

ВВЕДЕНИЕ

Численному моделированию удаленной регистрации в магнитосфере и на земной поверхности низкочастотных волн от искусственных возмущений в ионосфере уделяется большое внимание (см., например, [1–3]). Это связано с вводом в строй ряда мощных наземных коротковолновых (КВ) передатчиков, излучение которых модулировано на низкой частоте (НЧ). Излучение таких передатчиков (“нагревных стендов”) интенсивно поглощается в нижней ионосфере, приводя к возникновению эффекта нелинейной демодуляции и излучению НЧ-волн в ионосферу и волновод “Земля-ионосфера” [4]. Во всем мире функционирует более десятка нагревных стендов, наиболее мощными и оснащенными научной аппаратурой из которых являются стенды EISCAT (Норвегия) и HAARP (Аляска) [5, 6].

К настоящему времени можно считать, что теория распространения волн очень и сверхнизких частот (ОНЧ и СНЧ) в околоземном пространстве, в частности в волноводе (резонаторе) “Земля–ионосфера”, в основном построена [7–10]. Соответствующие алгоритмы часто физически недостаточно наглядны и требуют значительных усилий для интерпретации результатов. По этой причине задачу можно разбить на связанные части, для каждой из которых строятся относительно простые численные алгоритмы. Это относится, прежде всего, к алгоритмам построения лучевых траекторий в ионосфере, магнитосфере [11–14]. Даже в волноводе “Земля–ионосфера”, где на низких частотах приближение геометрической оптики формально не выполнимо вблизи верхней “стенки”, возможно построение эффективных численных алгоритмов в рамках “метода скачков” [15–17]. При этом коэффициенты отражения от нижней неоднородной анизотропной ионосферы необходимо вычислять независимо.

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

Цель данной работы – представить описание одного из возможных алгоритмов решения такой задачи. В качестве иллюстрации проведены расчеты в случае СНЧ (50…1000 Гц), которым в последнее время уделяется повышенное внимание для диагностики состояния земной коры дистанционными методами и в задачах связи.

1. ОСНОВНЫЕ СООТНОШЕНИЯ

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

Для расчетов напряженности поля используем уравнения Максвелла совместно с моделью среды в виде неоднородной холодной многокомпонентной магнитоактивной плазмы. Такая среда описывается тензором диэлектрической проницаемости ${\hat {\varepsilon }}$ (магнитная проницаемость ионосферной плазмы равна магнитной проницаемости вакуума μ₀). Тензор ${\hat {\varepsilon }}$ в системе координат, в которой ось z направлена вдоль внешнего геомагнитного поля, имеет вид [20]

где

ν_k – эффективная частота соударений частиц сорта k; f – частота волны, ω = 2πf, i – мнимая единица. Суммирование ведется по сорту k заряженных частиц, составляющих плазму, с учетом знака заряда в величинах Y_k:

где f_Hk и f_pk – гиро- и плазменная частота частицы сорта k: k = 1 (электроны), 2, … (ионы). Для электронов имеем

где N_e измеряется в см^–3; e, m_e – заряд и масса электрона; ε₀ – диэлектрическая проницаемость вакуума. В максимуме слоя F2 ионосферы имеем днем f_pe ≈ 10 МГц, ночью – f_pe ≈ 2…5 МГц в зависимости от условий (времени суток, сезона, географического положения и активности Солнца):

В пределах ионосферы пространственная структура вектора геомагнитного поля $\overrightarrow {{{B}_{0}}} $ хорошо описывается моделью точечного магнитного диполя, расположенного вблизи центра Земли с осью, наклоненной под некоторым углом к оси вращения Земли [21]. По аналогии с географическими координатами вводятся понятия геомагнитной широты Φ и долготы Λ. На высоте h и широте Φ значение гирочастоты можно оценить так:

где R₀ – средний радиус Земли (~6370 км). Для ионов сорта k в формулах (3), (4) следует использовать соответствующие значения концентрации N_k и массы m_k.

Пусть начало правой декартовой прямоугольной системы координат расположено в точке O на поверхности земли, а ось z направлена вертикально вверх. Плоская ЭМ-волна с волновым вектором $\vec {k}$ падает снизу на горизонтальный неоднородный по оси z слой плазмы под углом θ с вертикалью. Верхняя и нижняя границы слоя соответствуют высотам z₁ и z₂. Направим ось х так, чтобы вектор $\vec {k}$ находился в плоскости XZ. Направление локального геомагнитного поля определяется углом ψ вектора ($ - \overrightarrow {{{B}_{0}}} $) с вертикалью и задаваемым углом азимута A, который отсчитывается от оси х до направления проекции $\overrightarrow {{{B}_{0}}} $ на земную поверхность. Отсчет идет против часовой стрелки, если смотреть на поверхность сверху вниз. Величина угла ψ определяется из соотношения ctgψ = 2 tg Φ [21]. Мы будем использовать для описания направления $\overrightarrow {{{B}_{0}}} $ его направляющие косинусы {l, m, n}, где

Такой выбор системы отсчета позволяет записать любую переменную волновую величину $\vec {V}\left( {x,y,z,t} \right)$ в виде

где k₀ = ω/c – волновое число в вакууме; Γ = sin θ. Опуская везде данный экспоненциальный множитель, мы можем превратить систему уравнений в частных производных в систему обыкновенных дифференциальных уравнений первого порядка [18]:

где $\vec {e}$ – вектор-столбец вида {E_x, –Ey, Z₀H_x, Z₀H_y}, искусственно образованный из компонент E_x, E_y электрического, и компонент H_x, H_y магнитного поля волны, причем компоненты магнитного поля для уравнивания размерностей умножены на импеданс вакуума Z₀ = (µ₀/ε₀)^1/2 = 120π. Компоненты 4 × 4 матрицы $\hat {T}$ выражаются через компоненты 3 × 3 матрицы электрической восприимчивости $\hat {M} = \hat {\varepsilon } - \hat {1}$ ($\hat {1}~$ – единичная матрица):

где C = cos θ.

В выбранной выше системе координат легко показать, что элементы матрицы $\hat {M}$ имеют вид

где параметры S, P, D определены в уравнениях (1) с учетом ионного состава ионосферного слоя и эффективной частоты соударений частиц в слое.

Недостающие в уравнении (6) компоненты E_z и H_z находятся явно из уравнений Максвелла:

Для выделения падающих и отраженных волн в волноводе при z < z₂ запишем соотношения между компонентами полей в этой области:

где верхний знак соответствует восходящей волне, а нижний – нисходящей.

Система (6) состоит из четырех обыкновенных линейных дифференциальных уравнений первого порядка. Система однородна, т.е. имеет четыре линейно независимых решения ${{\vec {e}}_{j}}$ (j = 1…4), причем любое другое решение (6) есть линейная комбинация этих решений. Легко показать, что решения ${{\vec {e}}_{j}}$ в однородной среде соответствуют классическим магнитоионным модам, условно называемым обыкновенной и необыкновенной волнами, которые распространяются вверх и вниз [18].

Предположим, что при z > z₁ параметры ионосферы мало меняются на длине волны. В этом случае при падении волны на слой снизу на верхней границе слоя возникают только обыкновенная и необыкновенная волны, направленные вверх. Пусть это будут волны j = 1, 2. При z = z₁ решение (6) может быть представлено в виде $\vec {e}$ = a₁${{\vec {e}}_{1}}$ + a₂${{\vec {e}}_{2}}$, где a_1,2 – комплексные коэффициенты. В силу линейности и однородности системы (6) подобное представление решения справедливо для любой высоты z < z₁, т.е. коэффициенты a_1, 2 не зависят от z. Следовательно, необходимо найти только два линейно независимых решения (6), соответствующих восходящим обыкновенной и необыкновенной волнам. Для этого необходимо вычислить начальные значения ${{\vec {e}}_{{1,2}}}$(z₁) в однородной среде выше неоднородного слоя. Для существования таких решений требуется выполнение известного условия det($\hat {T}$ – q$\hat {1}$) = 0, где коэффициент q описывает изменение плоской волны вдоль оси z как exp(–ik₀qz). Зная значения q_1,2 на высоте z₁, можем проинтегрировать (6) вниз для соответствующих собственных векторов ${{\vec {e}}_{{1,2}}}$(z) матрицы $\hat {T}$.

Раскрывая определитель с учетом значений элементов матрицы $\hat {T}$ (7), получим известное комплексное алгебраическое уравнение четвертой степени относительно q (“квадрик Букера”) [18]. На рис. 1 схематически представлено расположение корней квадрика q_j (j = 1…4) на комплексной плоскости в случае низких частот, когда различие корней для четырех магнитоионных мод наиболее наглядно. Для восходящих волн Re(q) > 0. Из рисунка видно, что для “обыкновенной” волны |Im(q₁)| $ \gg $ Re(q₁), и эта волна относительно быстро затухает в неоднородной среде. Для “необыкновенной” волны, наоборот, |Im(q₂)| $ \ll $ $ \ll $ Re(q₂), затухание волны относительно невелико. По этой причине свистящие атмосферики распространяются на большие расстояния из полушария в полушарие в ионосфере и магнитосфере. Поэтому эту НЧ-моду часто называют “свистовой”.

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

1. На поверхности q в комплексном пространстве строим около начала координат правильный треугольник с характерным размером много меньше показателя преломления волны ${{\vec {e}}_{2}}$ на высоте z₁. В качестве характерного размера выбираем, например, величину радиуса описанной окружности.

2. Проверяем значения q в вершинах треугольника и находим наименьшее из них.

3. Переносим центр треугольника в эту вершину и уменьшаем его характерный размер.

При циклическом повторении процедуры 2, 3 треугольник неизбежно “падает” в один из корней, останавливаясь при достижении требуемой точности.

4. Степень полинома понижается на единицу, и процедуры 1–4 повторяются до вычисления всех четырех корней.

Этот простой алгоритм использовался нами в данной задаче продолжительное время и ни разу не выдал ошибочный результат.

Для численного интегрирования (6) используем стабильный метод численного интегрирования с переменным шагом и контролем ошибки интегрирования Рунге–Кутта–Мерсона [22], хотя возможно эффективное применение и других методов [23].

Начальные значения ${{\vec {e}}_{{1,2}}}$(z₁) для интегрирования системы уравнений (6) находим как собственные векторы матрицы $\hat {T}$ на высоте z₁ при известных собственных значениях q_j:

Линейная система уравнений (11) легко решается с точностью до постоянного множителя, например, E_y:

где индекс j не приводится,

На низких частотах обыкновенная волна сильно затухает вдоль оси z, следовательно, при интегрировании сверху вниз обыкновенная волна ${{\vec {e}}_{1}}$(z) быстро нарастает в сравнении с необыкновенной волной ${{\vec {e}}_{2}}$(z) (рис. 1). Из-за конечного числа разрядов в мантиссе вещественных чисел на ЭВМ этот эффект сильного нарастания одного из решений может привести к потере их линейной независимости, поэтому ее необходимо периодически восстанавливать в процессе интегрирования уравнений (6), пользуясь тем фактом, что ортогональные векторы всегда линейно независимы [24]. При достижении амплитуды решения ${{\vec {e}}_{1}}$(z) некоторого заданного порога, образуем вектор ${{\vec {e}}_{{20}}}$ = ${{\vec {e}}_{2}}$ + $\xi {{\vec {e}}_{1}}$ и потребуем, чтобы скалярное произведение ($\vec {e}_{1}^{*}$, ${{\vec {e}}_{{20}}}$) было равно 0. Отсюда находим множитель

где звездочка означает комплексное сопряжение. Поскольку вектор ${{\vec {e}}_{{20}}}$ есть линейная комбинация решений (6), то он тоже является решением, линейно независимым с ${{\vec {e}}_{1}}$. Отметим, что в сильно неоднородной среде понятие магнитоионной моды теряет свой изначальный физический смысл. В результате процедуры численного интегрирования, описанной выше, получаем при z = z₂ два формально правильных линейно независимых решения (6). Однако в результате применения процедуры ортогонализации эти решения в совокупности не имеют какого-либо физического смысла. С их помощью можно сконструировать все характеристики волнового поля в пустом пространстве при z ≤ z₂ [23].

Во-первых, можно разделить волны ниже слоя на восходящие и нисходящие части. Из выражений (7) и (8) легко получить значения q = ±C, причем, верхний знак соответствует волнам вверх ($\vec {U}$), а нижний – вниз ($\vec {D}$), т.е.

Возьмем электрическое поле $\vec {E}$, тогда $\vec {E}$ = $\vec {U}$ + $\vec {D}$ и ${\text{d}}\vec {E}$/dz = –ik₀C$\vec {U}$ + ik₀C$\vec {D}$. Выражение для производной можно также получить из уравнения (6). Получаем с учетом соотношений (10):

где C = cos θ.

Во-вторых, можно показать, что линейная комбинация ${{\vec {U}}_{{\text{п}}}}$ = ${{\vec {U}}_{2}}$ + $\zeta {{\vec {U}}_{1}}$, где ζ = –($\vec {U}_{1}^{*}$${{\vec {U}}_{2}}$)/($\vec {U}_{1}^{*}$${{\vec {U}}_{1}}$) представляет собой “проникающую” волну, которая порождает наименее затухающую необыкновенную волну выше уровня z₁. У линейно независимых решений ${{\vec {e}}_{1}}$ и ${{\vec {e}}_{{\text{п}}}}$ уравнения (6) ниже высоты z₂ есть простой физический смысл: они имеют соответственно наименьший и наибольший коэффициенты прохождения по мощности через неоднородный плазменный слой z₂ – z₁, причем ${{\vec {e}}_{{\text{п}}}}$ = ${{\vec {e}}_{2}}$ + $\zeta {{\vec {e}}_{1}}$.

В-третьих, векторы ${{\vec {e}}_{1}}$ и ${{\vec {e}}_{{\text{п}}}}$ можно использовать как базис для получения характеристик известной плоской волны, падающей на слой снизу (коэффициенты отражения, прохождения, в том числе по мощности). Традиционно рассматриваются падающие плоские волны двух независимых поляризаций. Вертикальная поляризация означает, что вектор $\vec {E}$ напряженности электрического поля волны лежит в плоскости (x, z), в которой лежит и вектор $\vec {k}$. В случае горизонтальной поляризации вектор $\vec {E}$ перпендикулярен этой плоскости, т.е. параллелен оси y. Обозначим эти векторы как ${{\vec {E}}_{\parallel }}$ и ${{\vec {E}}_{ \bot }}$ соответственно. Поскольку волны этих условно выбранных поляризаций должны быть линейными комбинациями базисных решений (${{\vec {e}}_{1}}$, ${{\vec {e}}_{{\text{п}}}}$), то из простейших геометрических соображений легко получить коэффициенты разложений ${{\vec {E}}_{\parallel }}$ и ${{\vec {E}}_{ \bot }}$ для волн вверх (↑) и вниз (↓), которые будут, очевидно, одинаковы для этих направлений: $\vec {E}_{\parallel }^{ \uparrow }$ = A_∥${{\vec {U}}_{1}}$ + B_∥${{\vec {U}}_{{\text{п}}}}$,

Если задана амплитуда E_∥ падающей волны вертикальной поляризации, то структуру поля ниже слоя можно определить из соотношений (13). Аналогично, если задана амплитуда поля E_⊥ падающей волны горизонтальной поляризации, то

Если на анизотропный слой падает волна, например, горизонтальной поляризации, то отразятся волны и горизонтальной, и вертикальной поляризаций. Поэтому принято представлять коэффициент отражения в виде 2 × 2 матрицы $\hat {R}$ [18]. Приведем определения элементов матрицы:

Из соотношений (10), (13), (14) и (15) получим

Определим поляризации ρ базовых решений:

поскольку ${{\vec {e}}_{1}}$, ${{\vec {e}}_{{\text{п}}}}$ ортогональны, то |ρ_п||ρ₁| = 1 и arg ρ₁ = arg ρ_п + π. Из разложений (13) и (14) получим

Коэффициенты прохождения по напряженности электрического поля волны легко получить, учитывая тот факт, что на уровень z₁ проходит только одна проникающая волна. Получаем два коэффициента:

Эти выражения можно переписать в более простом виде:

В заключение описательной части приведем выражения для коэффициентов прохождения по мощности. Мощность падающей на слой проникающей волновой моды вдоль оси z составляет

Выше слоя при z = z₁ вертикальная часть мощности формируется начальным полем свистовой моды ${{\vec {e}}_{{\text{п}}}}$:

Вертикальный коэффициент прохождения по мощности D_пz определяем так:

Определим полный коэффициент прохождения по мощности D_п:

где относительные “боковые” потоки мощности s_{x, y} на высоте z₁ вычисляются как отношения соответствующих компонент вектора Умова–Пойнтинга к его z-компоненте:

причем, компоненты напряженности поля волны E_z и H_z определяются из соотношений (9а) и (9б).

3. РЕЗУЛЬТАТЫ РАСЧЕТОВ

Далее везде будем использовать численные значения параметров плазмы, которые соответствуют приполярной ионосфере на широте 60° в равноденственный период при умеренной солнечной активности [25]. Процент ионных составляющих в слое [z₁, z₂] приведен в табл. 1.

Профили N_e(z) и ν_e(z) аппроксимировались кубическими сплайнами. Такой метод гарантирует отсутствие “паразитных” частичных отражений в сильно неоднородном слое нижней ионосферы. Для дневной ионосферы выбирали слой 50…100 км; для ночной ионосферы – 75…150 км. На верхней границе слоя плазменная частота для электронов днем равна 1268.2 кГц, ночью – 1750.6 кГц, гирочастоты электронов составляют днем 1507.1 кГц, ночью – 1472.7 кГц. Для ионов соответствующие частоты приведены в табл. 1. Частоты нижнего гибридного резонанса f_LHR равны соответственно 4.95 и 5.15 кГц.

Расчеты проводились для значений θ = 10°, A = = 180° (направление на магнитный юг) и частот в интервале 50…1000 Гц. Выбор этих величин определялся условиями работы станции HAARP. Однако следует учесть, что длина волны на частоте 50 Гц составляет в атмосфере 6000 км. При длине окружности экватора 40 000 км на ней укладывается шесть целых длин волн. Такой волновой процесс можно лишь условно считать “локальным”, а дальнейшее понижение частоты требует рассмотрения сферической полости “Земля–ионосфера” в виде резонатора, а не волновода. Если учесть, однако, что в нижней ионосфере на высоте z₁ показатель преломления волн относительно велик, то получим величину длины волны в среде ~30…40 км. Такую величину вполне можно считать “локальной”.

На рис. 2 приведены расчетные зависимости показателя преломления n свистовой моды от частоты при z = z₁ для дневной и ночной моделей ионосферы. Величина n для ночи выше дневных значений, поскольку ночная высота z₁ значительно больше дневной.

На рис. 3 представлена частотная зависимость полного коэффициента прохождения по мощности D (фактически соответствует величине D_п в тексте) для дневной и ночной моделей ионосферы. Для диапазона сверхдлинных дневные значения D, как правило, заметно меньше ночных из-за столкновительного поглощения волн на высотах 50…70 км в дневное время. Из рис. 3 видно, что подобная закономерность начинает развиваться только на частотах f > 500 Гц. На меньших частотах ионные резонансы в дневное время формируют характерный максимум прохождения на частотах f ~ 30 Гц, что соответствует примерно половине гиро-частоты наилегчайшего иона N⁺. Отметим, что величина D на частотах 0.2…1 кГц не мала ни днем, ни ночью, и почти достигает значения 0.4. Без учета влияния геомагнитного поля величина D равнялась бы нулю.

На рис. 4а приведены зависимости от частоты модулей компонентов матрицы $\hat {R}$ коэффициента отражения (см. (15), (16)) для дня (сплошные кривые) и ночи (пунктир). Приведены зависимости от f величин R₁₁ (кривые 1 и 4), R₁₂ (кривые 2 и 5) и R₂₂ для дня (кривая 3). Для выбранного значения азимута из рисунка видно (кривые 1, 3), что значения R₁₁ и R₂₂ очень близки, или практически совпадают (как и значения не диагональных элементов R₁₂ и R₂₁, из которых кривая приведена только для величины R₁₂). По этой причине на рис. 4б изображены только зависимости аргументов R₁₁(f) – кривые 1, 4; и аргументов R₁₂(f) – кривые 2 и 5.

Из рис. 4 видно, что в выбранном интервале частот преобладают диагональные элементы $\hat {R}$, причем в ночное время аргумент элементов матрицы относительно стабилен. В дневное время с ростом частоты при переходе из области ионных колебаний в область электронных волн аргумент диагональных элементов заметно меняется (кривая 1 на рис. 4б). Кривая 3 для аргумента R₂₂ на рис. 4б отсутствует, поскольку практически совпадает с кривой 1. Мы не приводим аналогичные зависимости поляризации отраженной волны в целом, поскольку в используемом интервале частот ионосфера играет роль почти идеального поляризационного фильтра, поэтому поляризация отраженной волны является круговой.

ЗАКЛЮЧЕНИЕ

Приведено описание алгоритма вычисления характеристик прохождения низкочастотных электромагнитных волн через магнитоактивную плоскослоистую нижнюю ионосферу с параметрами, аналогичными параметрам ионосферы над станцией HAARP, Аляска. Алгоритм основан на процедуре численного интегрирования волновых уравнений для случая падения плоской волны на ионосферу снизу с учетом ионного состава плазмы. В качестве иллюстрации работы алгоритма вычислены основные характеристики прохождения волны через нижнюю ионосферу в интервале частот 0.05…1 кГц.

Коэффициент прохождения по мощности указанных НЧ-волн достаточно велик и находится в интервале 0.2…0.4 для большей части используемых частот. В дневное время вблизи частоты ~300 Гц коэффициент прохождения в дневное время имеет максимум и превышает ночные значения примерно в два раза. Во всем исследованном интервале частот поляризация отраженной волны является круговой с хорошей степенью точности.