Численное моделирование характеристик радиосигналов в локально возмущенной неоднородной анизотропной ионосфере земли

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

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

Для компьютерного моделирования характеристик неоднородной среды, в которой распространяется электромагнитная волна, была, выбрана, полуэмпирическая модель распределения электронной концентрации в ионосфере IRI — 2001 [1].

Для построения зависимости плазменной частоты от пространственных координат использовались данные, получаемые в автоматическом режиме с интернет-сайта.

Данные, полученные с помощью интернет-ресурса, были сглажены с помощью кубических сплайнов [2] с окном от 0.1 до 0.25 МГц по плазменной частоте для обеспечения непрерывности первых производных от плазменной частоты по пространственным координатам.

3.    Математическая модель ионосферного возмущения

Модель возмущения, вызываемого бегущей вертикально вверх акустической волной, задавалась следующей формулой:

N = N q (1 + ^ К (H,т,t_m)).

V s

Здесь N и No — электронная концентрация возмущенной и невозмущенной ионосферы в точке с координатами Н (высота над поверхностью Земли) и дальностью от источника радиоизлучения т + ту. При этом т отсчитывается по горизонтали в плоскости трассы от центра возмущения, расположенного на дальности т у вдоль поверхности Земли в той же плоскости. Скорость звука V_s в (1) вычисляем по известной формуле для идеального газа:

V s =

■ -R g Т А2

Ц

Здесь R g — универсальная газовая постоянная, у — показатель адиабаты, который можно принять равным 7/5 (двухатомный газ). Зависимость молярной массы ц и температуры Т в (2) от высоты взяты по данным международной стандартной модели нейтральной атмосферы. Величина V_e в (1) характеризует амплитуду возмущения и имеет размерность скорости.

Безразмерный коэффициент К(H,т,t_m) характеризует положение центра возмущения в зависимости от времени Т_т и характер изменения интенсивности возмущения по мере удаления от центра возмущения. Он задается формулой

• ■ А1     ⁽ 2L V       ⁽ Нс - Н 1 А

  К (H,T,t m )= F

  
  

• [Ho ^exp   Uo/ J ^exp   ( H2 ) J .

  
  

Здесь Н _с = V_s • t_m — высота центра возмущения, поднимающегося вертикально вверх со скоростью V_s;

• h = Н — Y_c~ смещение относительно центра возмущения по высоте;

Но — полутолщина возмущения в вертикальной плоскости;

Ro = то + Нс1ду — «текущая» полуширина возмущения в горизонтальной плоскости,

• то — полуширина возмущения при t_m = 0.

• у — угол, характеризующий рост горизонтальных размеров возмущения при его подъеме вверх.

В формуле (3) множитель F (h/Но) задает форму волнового возмущения. При выполнении расчётов в качестве F использовались функции F1 и F2:

• ^F1 = [7 ⁺ 71] ^exP [" (770) 1 ,  ^F2 = ^k1 (^exP [" (770 ⁺ -2) 1 ^{— k}2 ^exP ["^k3 (7) 1} .

4. Моделирование характеристик радиосигналов

Здесь ki = ±1, k2 и кз — положительные коэффициенты, 71 и 72 — параметры в формулах.

В рамках метода бихарактеристик [3] задача определения амплитудно-фазовой структуры волнового поля для квазимонохроматического сигнала сводится к решению системы из шести дифференциальных уравнений:

с2т,   ш де ш ^кх 2 дк _ж

• с² 1      ш де

• ш ^к У   2 дк _у

Ах_

Ат =

АУ = Г

Ат

1   Гс- । ^ш де 1   ¹           Ак _х _ ш де Щ , ш де 1   ¹                /гр

• ^{£ +} 2 дш ,         Ат = 2 дх ^{Е +} 2 дш ,

].Г ш8е I   '        ^   ^шд Г ш8е 1   '

• ^{ь +} 2 дш ,         Ат = 2 ду ^{Ь +} 2 дш ,

dz

Дт

— Г С2В     ^{ш д}    [ шд^ 1

= ш ^кг 2 8k_z\ I^{е +} 2 8ш J

Дк ₂

Дт

= ш 1 [е +21 1^-1

с начальными условиями Коши:

ж|_Г0=0 = $0, к_х\_Г0=о = cosa cos^, < y|ro=о = Уо, к_у |_Го=о = sin a cos Д

, гк=о = z q , к_г|го=о = sin Д.

Здесь г = ($, у, z) — это декартовые координаты текущей точки на траектории луча;

• к = (к_х, к у , k_z ) — волновой вектор в той же точке;

• ш — круговая частота сигнала в радианах; с — скорость света в вакууме;

• т — параметр размерности времени вдоль траектории;

• ($о,Уо,2о) — координаты источника излучения;

• а, Д — угловые координаты траектории;

е = е ( т, к,ш) — эффективная диэлектрическая проницаемость среды, в которой распространяется волна.

Для контроля точности интегрирования системы дифференциальных уравнений (5) -

• (7) используется дисперсионное соотношение:

к_х + ку + к2 — ( ^Ш ) е(г, к,ш) = 0.                            (8)

Для магнитоактивной среды без учета влияния поглощения на форму лучевых траекто рий эффективная диэлектрическая проницаемость может быть представлена в следующем виде [4-6]:

е ( г, к) = 1 —

2X (1 — X )

2(1 — X) — УТ ± ^УТ + 4(1 — X)2Уь2

Здесь X = (шм/ш}²', ш2 = Х/е о ш — квадрат угловой плазменной частоты среды; епт — заряд и масса электрона; Ео — электрическая постоянная;

• X — концентрация электронов;

У ь = У cos 9, У т = У sin 9, где У = (ш_нш);

• ш н = (|е|/т)Во — угловая гирочастота в радианах;

Во — модуль индукции внешнего магнитного поля;

• 9 — угол между направлением вектора к и вектора Во.

Отметим, что в формуле (9) знак «+» соответствует обыкновенной волне, а знак «-» необыкновенной.

Для построения континуума лучевых траекторий в фазовом пространстве был разработан на языке СД комплекс программ ВЕАМН. При создании программы, выполняющей численное интегрирование системы обыкновенных дифференциальных уравнений (5) - (7), использовался модифицированный метод Кутта-Мерсона.

Рассмотрим результаты численного моделирования лучевого распространения для обыкновенной волны в ионосфере Земли, когда диэлектрическая проницаемость среды описывается формулой (9).

При выполнении расчётов использовалась декартовая система координат, для которой в точке расположения передатчика ось У направлена по местной вертикали вверх. Начало системы координат находится на поверхности Земли. Ось X выбрана в плоскости большого круга (плоскость, проходящая через центр Земли, передатчик и приёмник). Ось Z дополняет оси X и У до правой тройки. Частота передатчика, расположенного на поверхности Земли, была выбрана равной 14 МГц. Диаграмма направленности передатчика — это узкий лепесток в азимутальной плоскости, сориентированный в плоскости большого круга. Расчёты выполнялись для 14 часов московского декретного времени для сезона с высокой солнечной активностью.

На рисунках 1 и 2 представлены проекции лучевой структуры обыкновенной волны, первоначально излучаемой точечным источником в плоскости большого круга, соответственно на плоскость (X ,Y) (рис. 1) и плоскость (X , Z) (рис. 2). Комплекс программ ВЕАМН, выполняя интегрирование лучевых траекторий с автоматическим выбором шага интегрирования по переменной т, обеспечивал выполнение дисперсионного соотношения (8) с точность нс хуже чем 10—12.

На рис. 1 видно, что все лучи отражаются либо от ионосферного слоя Е, либо от слоя F1, либо от слоя F2. Учитывая малость отклонений лучей от плоскости большого круга, на основе анализа рисунка 1 можно заключить, что лучи, отражающиеся от слоёв F1 л F2. образуют сложную каспонднуто структуру «бабочка» — волновая катастрофа А5. Лучи, отражающиеся от слоя Е, образуют каспоидную структуру «остриё» — волновая катастрофа А3.

Рис. 1

Рис. 2

Согласно методу бихарактеристик [3, 7] волновое поле определяется по формуле

Д                     1 / 2

Г ^ом(х,у,г) = С-^гсом £ А(« к ,Р к ) IJ I х ^к=1                  ^,                                  (10)

X exp {—G(aj,Рк) + г [ЛФ(«к,Рк) - Mj^ } , где N — число лучей, приходящих в точку с координатами (х, у, г);

а к , Рк ~ угловые координаты луча с номером к',

Mj — число каустических точек вдоль лучевой траектории с номером к:

Сге<,м = 103 • (^) ; Л = 10ш/3 — числовые константы:

J — якобиан, который определяется по следующей формуле:

J=

х_а х Д К_х У а Ур ^К У ^ a Р р ^K

(И)

Отметим, что при перемещении вдоль траектории луча якобиан меняет знак после прохождения каждой каустической точки.

Ф(а, Р ) — это набег фазы, определяемый как интеграл вдоль лучевой траектории.

г ( а,р

Ф(а,Р )= j

kdx

= [ ^-^Щ) _dT

(^ + ^    • 2)    ’

где T (а, Р) — значение независимой переменной t в точке с ко ординатами (х,у,гУ

G(a, Р ) — интегральное поглощение вдоль этой же траектории [5];

А(а,Р ) — коэффициент, определяемый по диаграмме направленности передающей антенны.

Расчет поля по формуле (10) возможен всюду за исключением областей, где якобиан (11) близок к нулю. Поверхность, на которой этот якобиан равен нулю, является огибающей семейства лучевых траекторий в пространстве (х, у, г) и называется каустикой. В окрестности каустик и каустической тени метод геометрической оптики неприменим. Для расчета волнового поля в этих областях можно использовать метод канонического оператора Маслова [7-9].

5.    Метод Доплера

Доплеровское смещение частоты мы вычисляем как первую производную по времени (например, по tm) от набега фазы вдоль лучевой траектории, связывающей приёмник и передатчик [5]:

, = - ± .dl

^D    2л dt_m

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

Зафиксировав время t_m и определив для него параметры ионосферы, мы находим пару значений (t_m, Ф) для каждой лучевой траектории, соединяющей приёмник и передатчик. Проведя аналогичные расчёты для моментов времени t_m — dT nt+ + dT, мы получаем возможность вычислить производную по времени t_m от фазы по формуле

^D = - 4^dT [Ф(tm + dT) — Ф(tm — dT)]> то есть вычислить доплеровское смещение частоты для одной моды принимаемого сигнала.

Пример расчёта вариаций доплеровского смещения частоты, вызванных акустическим возмущением, перемещающимся вертикально вверх в атмосфере Земли, представлен на рисунке 3. Расчёты выполнены для F = F2 при ki = 1, k2 = 1.9, кз = 4 и 72 = -0.1. При таких параметрах характерные горизонтальные и вертикальные размеры неоднородности, вызванной акустическим возмущением, составили 10-40 км, а отклонение электронной концентрации от её значения в невозмущенной ионосфере были около 1-2%. Частота излучаемого сигнала — 12 мегагерц. По горизонтальной оси отложено время t_m в минутах, а по вертикальной оси доплеровское смещение частоты в герцах fo- Время t_m = 0 соответствует 14 часам московского декретного времени. Ненулевые значения J d до и после прохождения акустического возмущения обусловлены регулярными суточными вариациями ионосферы, не связанными с прохождением возмущения. Эффект акустической волны «накладывается» на доплеровское смещение частоты, обусловленное естественными вариациями интегральной электронной плотности ионосферы вдоль траектории распространения волны с течением времени.

Следует отметить, что на исследуемой трассе отчетливо воспроизводится характерная «седлообразная форма» зависимости доплеровского смещения частоты от времени, что достаточно хорошо согласуется с результатами экспериментальных измерений, описанных в работах [11-13].

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

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

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

Следует также отметить, что по данным работы [11] некоторые характерные изменения в ионосфере, связанные с сейсмической активностью, могут наблюдаться задолго до регистрации самого землетрясения (период времени порядка 2-3 дней). Данный эффект позволяет использовать метод доплеровского зондирования для упреждающего прогноза возможных негативных событий.

Работа выполнена при поддержке ФЦП «Научные и научно-педагогические кадры инновационной России» на 2009-2013 годы.