CC BY
22
УДК 621.371
Е.Б. Ипатов1, Е. А. Палки»}’2, В. И. Чивилёв1, Д. Е. Ипатов1
1 Московский физико-технический институт (государственный университет)
2 Российский новый университет
Численное моделирование характеристик радиосигналов в локально возмущенной неоднородной анизотропной ионосфере Земли
Исследуется задача о моделировании влияния воздействия акустической волны на характеристики радиосигнала, распространяющегося на коротковолновой трассе в неоднородной анизотропной ионосфере Земли. Результаты расчётов показывают наличие доплеровского смещения частоты при прохождении радиосигнала через движущееся акустическое возмущение, что подтверждается экспериментально.
Ключевые слова: ионосфера, регулярно неоднородная среда, распространение радиоволн, численное моделирование, метод бихарактеристик, каустические особенности.
1. Введение
Данная работа посвящена исследованию влияния акустических атмосферных возмущений на характеристики квазимонохроматического импульса, распространяющегося в ионосфере Земли. В основе исследования положен метод численного моделирования, причём используемая модель ионосферы позволяет выбирать актуализированные параметры с учётом солнечной активности, географических координат и магнитного поля Земли. Такая постановка определяется необходимостью решения прикладных задач радионавигации, радиолокации, а также потребностями мониторинга различных геофизических процессов, протекающих в земной коре и нейтральной атмосфере, соприкасающейся с ионосферой.
Чувствительность ионосферной плазмы к различным геофизическим процессам может быть использована для выявления новых физических закономерностей, для решения ряда прикладных задач, связанных с хозяйственной деятельностью человека и безопасностью его существования.
2. Модель ионосферы, используемая в расчётах
Для компьютерного моделирования характеристик неоднородной среды, в которой распространяется электромагнитная волна, была выбрана полуэмпирическая модель распределения электронной концентрации в ионосфере 1Ш — 2001 [1].
Для построения зависимости плазменной частоты от пространственных координат использовались данные, получаемые в автоматическом режиме с интернет-сайта modelweb.gsfc.nasa.gov/models/iri.html для трассы, задаваемой двумя точками на поверх-
ности Земли. Для исследования была выбрана среднеширотная трасса, проходящая над территорией Российской Федерации. Географические координаты передатчика: 55.8° с.ш. и 37.6° в.д.; приёмника: 38° с.ш. и 5° в.д.
Данные, полученные с помощью интернет-ресурса, были сглажены с помощью кубических сплайнов [2] с окном от 0.1 до 0.25 МГц по плазменной частоте для обеспечения непрерывности первых производных от плазменной частоты по пространственным координатам.
3. Математическая модель ионосферного возмущения
Модель возмущения, вызываемого бегущей вертикально вверх акустической волной, задавалась следующей формулой:
N _ N0(1 + ^К(Н, г,іт)).
V*
(1)
Здесь N и N0 — электронная концентрация возмущенной и невозмущенной ионосферы в точке с координатами Н (высота над поверхностью Земли) и дальностью от источника радиоизлучения г + г ь- При этом г отсчитывается по горизонтали в плоскости трассы от центра возмущения, расположенного на дальности гь вдоль поверхности Земли в той же плоскости. Скорость звука У3 в (1) вычисляем по известной формуле для идеального газа:
1/2
(2)
Здесь Кд — универсальная газовая постоянная, 7 — показатель адиабаты, который можно принять равным 7/5 (двухатомный газ). Зависимость молярной массы у и температуры Т в (2) от высоты взяты по данным международной стандартной модели нейтральной атмосферы. Величина Уе в (1) характеризует амплитуду возмущения и имеет размерность скорости.
Безразмерный коэффициент К(Н,г,Ьт) характеризует положение центра возмущения в зависимости от времени Тт и характер изменения интенсивности возмущения по мере удаления от центра возмущения. Он задается формулой
к
н0
ехр
— (І)
ехр
(3)
Здесь Нс = ■ Ьт — высота центра возмущения, поднимающегося вертикально вверх
со скоростью
к = Н — Ус — смещение относительно центра возмущения по высоте;
Но — полутолщина возмущения в вертикальной плоскости;
= го + Hc^gр — «текущая» полуширина возмущения в горизонтальной плоскости,
Го — полуширина возмущения при = 0,
р — угол, характеризующий рост горизонтальных размеров возмущения при его подъеме вверх.
В формуле (3) множитель Р(Н/Но) задает форму волнового возмущения. При выполнении расчётов в качестве Р использовались функции ^ и
*1 =
I \^и/ I у^и / у/ I
(4)
Здесь &1 = ±1 к2 и кз — положительные коэффициенты, 71 и 72 — параметры в формулах.
ї +7/ ехр - ( Щ) , ^2 _ &1 |ехр - (її+< — к2ехр —к:і (ї)
4. Моделирование характеристик радиосигналов
В рамках метода бихарактеристик [3] задача определения амплитудно-фазовой структуры волнового поля для квазимонохроматического сигнала сводится к решению системы из шести дифференциальных уравнений:
йх
(ІТ
йу
йт
с2 и ш дг
ш 2 дкх
' + ш де_ ] 1
Ь + 2 дш\
' [■
є + шдє_ 1 1
Ь + 2 дш\
йкх + шдє_ ] 1
йт 2 дх I + 2 дш\
ш дє 1 1
йку ш д£ Г _ + ш дє_
йт 2 ду V + 2 дш
(5)
(6)
У
dz
dr
с2 к ш де
шк* 2 dkz
Г£ I 1 1 dkz _ шд§_ \£ + ^de_"I 1
\£ + 2 дш\ , dr _ 2 dz L£ + 2 дш\ V Ч
с начальными условиями Коши:
х\Г0=о _ хо, кх1го=о _cosacos3, у\г 0=о _ Уо, ку |г0=о _ sin a cos 3, z\r0=о _ -го, kzlr0=о _ sin3.
Здесь г _ (х, у, z) — это декартовые координаты текущей точки на траектории луча; к _ (кх, ку, kz) — волновой вектор в той же точке;
ш — круговая частота сигнала в радианах; с — скорость света в вакууме;
(хо, Уо, %о) — координаты источника излучения;
a, 3 — угловые координаты траектории;
е _ е(г, к, ш) — эффективная диэлектрическая проницаемость среды, в которой распространяется волна.
Для контроля точности интегрирования системы дифференциальных уравнений (5) -(7) используется дисперсионное соотношение:
кХХ + ку2 + kl2 — е(г,к,ш) = 0. (8)
Для магнитоактивной среды без учета влияния поглощения на форму лучевых траекторий эффективная диэлектрическая проницаемость может быть представлена в следующем виде [4-6]:
е( г,к) = 1----------------2Х (1 — X) . ^
2(1 — X) — ±^У^ + 4(1 — X)2Fl2
Здесь X = (ш^/ш)2; ш2м = N/е0т — квадрат угловой плазменной частоты среды;
е и т — заряд и масса электрона; £о — электрическая постоянная;
N — концентрация электронов;
Уь = У cos0, Уг = У 81п 0, где У = (шдш);
шн = (|е|/т)В0 — угловая гирочастота в радианах;
Во — модуль индукции внешнего магнитного поля;
в — угол между направлением вектора к и вектора Во-
Отметим, что в формуле (9) знак «+» соответствует обыкновенной волне, а знак «-» необыкновенной.
Для построения континуума лучевых траекторий в фазовом пространстве был разработан на языке С# комплекс программ ВЕАМН. При создании программы, выполняющей численное интегрирование системы обыкновенных дифференциальных уравнений (5) - (7), использовался модифицированный метод Кутта-Мерсона.
Рассмотрим результаты численного моделирования лучевого распространения для обыкновенной волны в ионосфере Земли, когда диэлектрическая проницаемость среды описывается формулой (9).
При выполнении расчётов использовалась декартовая система координат, для которой в точке расположения передатчика ось У направлена по местной вертикали вверх. Начало системы координат находится на поверхности Земли. Ось X выбрана в плоскости большого круга (плоскость, проходящая через центр Земли, передатчик и приёмник). Ось Z дополняет оси X и У до правой тройки. Частота передатчика, расположенного на поверхности Земли, была выбрана равной 14 МГц. Диаграмма направленности передатчика — это узкий лепесток в азимутальной плоскости, сориентированный в плоскости большого круга. Расчёты выполнялись для 14 часов московского декретного времени для сезона с высокой солнечной активностью.
На рисунках 1 и 2 представлены проекции лучевой структуры обыкновенной волны, первоначально излучаемой точечным источником в плоскости большого круга, соответственно
на плоскость (Х,У) (рис. 1) и плоскость (Х^) (рис. 2). Комплекс программ ВЕАМН, выполняя интегрирование лучевых траекторий с автоматическим выбором шага интегрирова-
не хуже чем 10-12.
На рис. 1 видно, что все лучи отражаются либо от ионосферного слоя Е, либо от слоя ^1, либо от слоя Е2. Учитывая малость отклонений лучей от плоскости большого круга, на основе анализа рисунка 1 можно заключить, что лучи, отражающиеся от слоёв и ^2, образуют сложную касиоидную структуру «бабочка» — волновая катастрофа А$.
Е
катастрофа А3.
ОН--------------1--------------1-------------1--------------1--------------1-------------1------------Ч- ’
400 500 600 700 800 900 1000 1100
Рис. 1
Рис. 2
Согласно методу бихарактеристик [3, 7] волновое поле определяется по формуле
N
игео™(х,у, г) = с™* Е А(ак,$к) к=1
cos(«fc)
1/2
(10)
х ехр{-С(ак,Рк) + * [ЛФ(ак, Зк) - Мк%]} ,
где N — число лучей, приходящих в точку с координатами (х,у, г); а к, Зк — угловые координаты луча с номером к;
Мк — число каустических точек вдоль лучевой траектории с номером к; Сгеом = 103 ■ (^) ; Л = 10ш/3 — числовые константы;
J — якобиан, который определяется по следующей формуле:
у» у'а К
г'а Кг
(11)
Отметим, что при перемещении вдоль траектории луча якобиан меняет знак после прохождения каждой каустической точки.
Ф(а, З) — это набег фазы, определяемый как интеграл вдоль лучевой траектории.
Т(а,')
Ф(а, З) = / к(1х =
ш(е -км ■1)
йт,
(12)
где Т(а, З) — значение независимой переменной т в точке с координатами (х,у, г);
С(а,З) — интегральное поглощение вдоль этой же траектории [5];
А(а, З) — коэффициент, определяемый по диаграмме направленности передающей антенны.
Расчет поля по формуле (10) возможен всюду за исключением областей, где якобиан (11) близок к нулю. Поверхность, на которой этот якобиан равен нулю, является огибающей семейства лучевых траекторий в пространстве (х,у, г) и называется каустикой. В окрестности каустик и каустической тени метод геометрической оптики неприменим. Для расчета волнового поля в этих областях можно использовать метод канонического оператора Маслова [7-9].
х
J
о
5. Метод Доплера
Доплеровское смещение частоты мы вычисляем как первую производную по времени (например, по Ьт) от набега фазы вдоль лучевой траектории, связывающей приёмник и передатчик [5]:
= 1_ йф
° 2ж Мт'
Здесь Ф — набег фазы, который определяется по формуле (12) как интеграл вдоль траектории, связывающей приёмник и передатчик.
Зафиксировав время и определив для него параметры ионосферы, мы находим пару значений (£т, Ф) для каждой лучевой траектории, соединяющей приёмник и передатчик. Проведя аналогичные расчёты для моментов времени — йТ и + йТ, мы получаем возможность вычислить производную по времени от фазы по формуле
= — 47ГйТ[Ф(^т + йТ) — Ф(— йТ^
то есть вычислить доплеровское смещение частоты для одной моды принимаемого сигнала.
Пример расчета вариаций доплеровекого смещения частоты, вызванных акустическим возмущением, перемещающимся вертикально вверх в атмосфере Земли, представлен на рисунке 3. Расчёты выполнены для Р = Р2 при к1 = 1, к2 = 1.9, кз = 4 и 72 = —0.1. При таких параметрах характерные горизонтальные и вертикальные размеры неоднородности, вызванной акустическим возмущением, составили 10 40 км, а отклонение электронной концентрации от её значения в невозмущенной ионосфере были около 1 2%. Частота излучаемого сигнала — 12 мегагерц. По горизонтальной оси отложено время в минутах, а по вертикальной оси доплеровское смещение частоты в герцах /д. Время Ьт = 0 соответствует 14 часам московского декретного времени. Ненулевые значения /д до и после прохождения акустического возмущения обусловлены регулярными суточными вариациями ионосферы, не связанными с прохождением возмущения. Эффект акустической волны «накладывается» на доплеровское смещение частоты, обусловленное естественными вариациями интегральной электронной плотности ионосферы вдоль траектории распространения волны с течением времени.
1,2
0,8'
0,4
-0,4'
-0,8
-1,2
-1,6
{ ГЦ
Л
/ '
/
5 / 7 < Э 1 1 13 1
1т, МИН
Рис. 3
Следует отметить, что на исследуемой трассе отчетливо воспроизводится характерная «седлообразная форма» зависимости доплеровекого смещения частоты от времени, что достаточно хорошо согласуется с результатами экспериментальных измерений, описанных в работах [11 13].
6. Заключение
Разработанный комплекс алгоритмов расчета пространственно-временной структуры радиосигнала позволяет проводить полномасштабное численное моделирование экспериментов по распространению сложных радиосигналов в ионосферной плазме. Производительность предлагаемых алгоритмов позволяет применять их также для решения таких задач, как распространение широкополосных и еверхширокополоеных сигналов, формирование полей ближней зоны в неоднородных средах, диагностика ионосферной плазмы импульсами специальной формы.
Для оценки принципиальной возможности обнаружения ионосферных возмущений, вызванных акустической волной от наземных или подземных источников (землетрясения, взрывы), было проведено компьютерное моделирование, для которого численными методами на основе метода бихарактеристик было рассчитано доплеровское смещение частоты при распространении электромагнитной волны через возмущенную ионосферу. Расчёты показали, что величину доплеровекого смещения определяет не столько величина изменения электронной концентрации, сколько вариации градиентов электронной концентрации, обусловленные возмущением.
Следует также отметить, что по данным работы [П] некоторые характерные изменения в ионосфере, связанные с сейсмической активностью, могут наблюдаться задолго до регистрации самого землетрясения (период времени порядка 2-3 дней). Данный эффект позволяет использовать метод доплеровекого зондирования для упреждающего прогноза возможных негативных событий.
Работа выполнена при поддержке ФЦП «Научные и научно-педагогические кадры инновационной России» на 2009-2013 годы.
Литература
1. Bilitza D. International Reference Ionosphere // Radio Sci. — 2001. — V. 36, N 2. — P. 261— 275.
2. Стечкин С. Б., Субботин Ю. Н. Сплайны в вычислительной математике. — М.: Наука, 1976.
3. Лукин Д. С., Спиридонов Ю. Г. Применение метода характеристик для численного решения задач распространения радиоволн в неоднородной и нелинейной среде // Радиотехника и электроника. — 1969. — Т. 14, № 9. — С. 1673-1677.
4. Гинзбург В. Л. Распространение электромагнитных волн в плазме. — М.: Наука, 1967.
5. Дэвис К. Радиоволны в ионосфере. — М.: Мир, 1973.
6. Александров А. Ф., Богданкевич Л. С., Рухадзе А. А. Основы электродинамики плазмы: учебник для физических специальностей университетов. — М.: Высш. ж., 1988.
7. Маслов В. П., Федорюк М. В. Квазиклассическое приближение для уравнений квантовой механики. — М.: Наука, 1976.
8. Ипатов Е. Б., Лукин Д. С., Палкин Е. А. Численная реализация метода канонического оператора Маслова в задачах распространения коротких радиоволн в ионосфере Земли // Известия вузов. Радиофизика. — 1990. — Т. 33, Л*8 5. — С. 562-573.
9. Ipatov Е.В., Lukin D.S. and Palkin E. A. Maslov canonical operator in problems of numerical simulation of diffraction and propagation of waves in inhomogeneous media // Soviet journal of numerical analysis and mathematical modeling. VNU Science press BV. — 1990. - V. 5, N 6. - P. 465-488.
10. Ипатов Е.Б., Глушнёв С. Ф., Чивилёв В. И., Палкин Е.А. Моделирование дисперсионных эффектов радиосигналов в неоднородной ионосфере Земли // Радиотехника и электроника. - 2003. — Т. 48, № 12. — С. 1436-1442.
11. Таращук Ю. Е., Нагорский П. Н., Борисов Б. Б. и др. Нестационарные процессы в ионосфере Земли и их влияние на распространение коротких радиоволн // Томск: изд. ТРУ. - 1986. - С. 163.
12. Georges Т. Н. HF Doppler studies of traveling ionospheric disturbances // J. Atm. and Terr. Phvs. - 1968. - V. 30. - P. 1441-1481.
13. Bartning Group delay variations of oblique incidence resulting from medium scale traveling ionospheric disturbances // Radio Sci. — 1973. — V. 13. — P. 147-157.
Поступила в редакцию 15.06.2011.
