Разработка модели распределения плотности токов при возбуждении ионосферы высокочастотным облучением Текст научной статьи по специальности «Электротехника, электронная техника, информационные технологии»

ЭЛЕКТРОННЫЕ И ЭЛЕКТРОМАГНИТНЫЕ

УСТРОЙСТВА

УДК 51.74: 533.9: 537.5: 537.876

А. Ю. Гришенцев, А. Г. Коробейников

РАЗРАБОТКА МОДЕЛИ РАСПРЕДЕЛЕНИЯ ПЛОТНОСТИ ТОКОВ ПРИ ВОЗБУЖДЕНИИ ИОНОСФЕРЫ ВЫСОКОЧАСТОТНЫМ ОБЛУЧЕНИЕМ

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

Ключевые слова: ионосфера, распространение электромагнитных колебаний, математическое моделирование.

Введение. Моделирование протекания в ионосфере динамических процессов, связанных с возбуждением электромагнитного поля, актуально для радиосвязи, метеорологии и ряда других научно-практических направлений [1—8]. Моделирование распространения радиоволн в ионосфере является сложной многопараметрической задачей, пока не имеющей универсального решения [1—11]. Задача усложняется особенностями плазменного состояния газа, пространственной и динамической неоднородностью ионосферы, неоднородностью геомагнитного поля, наличием собственных токов в ионосфере и др. Модель, рассмотренная в настоящей работе, предназначена для численного анализа и изучения переходных процессов, а также пространственного распределения плотности тока в ионосфере при искусственном облучении высокочастотным радиоизлучением.

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

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

Область ионосферы разделена по горизонтали на К коаксиальных слоев (от нуля до К-1) поверхностями сечения, имеющими форму усеченного сферой конуса с шаровым сегментом в основании и относительно малым углом при вершине (такие усеченные конусы, при некоторых допущениях, можно считать цилиндрами); по вертикали — на М слоев (от нуля до М-1) поверхностями сечения, имеющими форму шарового сегмента. Таким образом, ионосфера дробится на множество коаксиальных цилиндрических трубок общим числом МЯ=(М+1)К,

каждая из которых характеризуется радиусом гтк, толщиной Агтк и высотой Ит, где

т=0, 1, 2, ..., М-1 и к=0, 1, 2, ..., К-1. Область тропосферы и стратосферы возможно объединить в единый слой М и разделить на К-1 коаксиальных трубок высотой ЬМ, которые сверху примыкают к соответствующим трубкам из слоя ионосферы М-1, а снизу опираются на поверхность Земли.

Выбор значений Ит для последующего моделирования принимается с учетом изменения следующих величин как функций высоты: значение электронной плотности Ые, эффективного числа упругих столкновений электронов с нейтральными частицами и ионами, температуры (в том числе температурных инверсий), химического состава и плотности. Значения Агтк и Ятк выбираются так, чтобы можно было оценить неравномерность распределения плотности тока по ионосфере, вызванную проявлением поверхностного эффекта, т.е. с учетом глубины проникновения 5 электромагнитной волны:

где ю — круговая частота, ц — абсолютная магнитная проницаемость вещества, у — удельная объемная электропроводность вещества, X — длина волны. Наиболее значительные изменения с высотой (И) аргументов функции (1) будут наблюдаться для у в зависимости от значений Ые и Уей- [12]. При выборе Агтктребуется оценивать уточненное значение 5 с наибольшими значениями у для слоев ионосферы, в которых необходимо учитывать неравномерность распределения плотности тока. Окончательное уточнение геометрических параметров модели производится исходя из требуемой точности и доступных вычислительных ресурсов, а также с учетом особенностей расчета эквивалентных электрических параметров.

Эквивалентная схема замещения. Эквивалентная схема (рис. 2) получается в результате замены каждого дискретного элемента модели набором КЬС с определенными параметрами, значения которых могут динамически меняться во времени. В данной схеме элементы имеют коаксиальную структуру, каждый (см. рис. 1) обладает собственной индуктивностью Ь(тК+к) и активным сопротивлением Я(тК+к), где т=0, ..., М-1, к=0, ..., К-1. Кроме того, присутствуют межслойные электрические взаимодействия благодаря емкостным и индуктивным связям, индуктивно в модели взаимосвязаны все катушки, между соседними слоями присутствуют емкостные связи С(т(К-1)+к). Помимо емкостных связей, обусловленных токами смещения, присутствуют горизонтальные межслойные связи, обусловленные токами проводимости Я(т(К-1)+к). Такой подход позволяет рассматривать модель в виде схемы, состоя-

ние. 1

(1)

щей из обычных Л£С-элементов, и анализировать переходные процессы и устоявшиеся режимы стандартными методами электротехники.

R(K-2+NR) С(К-2) R(K-1)

L(K-1)

R(m(K-1)+K-2+NR) C(m(K-1)+K-2)

R(mK+K-1)

1 R(m(K-1)+k+NR) tC(m(K-1)+k)

R(mK+k)

R(0+NR)

C(Q)

R(0)

L(0)

I

•R(m(K-1)+NR) • C(m(K-1))

I R(mK)

L(mK+K-1){

L(mK+k)

L(mK)

R(M(K-1)+K-2+NR) C(M(K-1)+K-2) R(MK+K-1)

R(MK+K-1) C((M+2)(K-1))

C((M+1)(K-1)+k)

1 R(M(K-1)+NR) ,C(M(K-1)) R(MK)

;L(MK)

C((M+1)(K-1))

T—• • • I—I—I

Значение: NR=(M+1)K.

Рис. 2

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

Для составления системы дифференциальных уравнений первого порядка каждый элемент цепи выделим в отдельную ветвь. При отнесении ветвей графа схемы к дереву или связям будем придерживаться правила получения нормального дерева и нормального подграфа связей. К ветвям дерева последовательно отнесем ветви с источниками ЭДС, затем — ветви с конденсаторами, ветви с резисторами и в последнюю очередь — ветви с катушками индуктивности. В качестве связей (хорд) сначала выделим источники тока (в приведенном варианте схемы отсутствуют), затем — индуктивные элементы, резистивные ветви и в последнюю очередь — ветви с конденсаторами. Граф, полученный таким образом, называют нормальным [13].

Общее число узлов в графе схемы [(2M +1) K + 2], общее число ветвей [[ (4M +1) +1 - 2M ].

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

минимальный набор полей: порядковый номер строки (далее соответствует номеру ветви), номера узлов стока и истока, тип и номер элемента в соответствии с электрической схемой. Далее путем перебора строк DATA получим матрицу A, строки которой соответствуют номерам узлов, а столбцы — номерам ветвей. Номера ветвей соответствуют номеру элемента в массиве DATA.

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

Размерность матрицы A [(2M + 1)K + 2 - 1]х[K(4M +1) +1- 2M], где [(2M + 1)K + 2 -1] — число узлов, [[ (4M +1) +1 - 2M ] — число ветвей.

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

А = ||Ад Ас ||. (2)

С учетом того что таблица DATA была получена в соответствии с правилами формирования нормального графа, матрица Ад будет содержать набор данных, соответствующий нормальному графу. В полученных матрицах дерева Ад и связей Ас номера ветвей могут быть расположены не по порядку, что может вызывать некоторые неудобства при хранении в памяти ЭВМ и вычислениях, привести нумерацию к упорядоченному виду можно с помощью упорядоченного перебора и присвоения новых номеров ветвям. В таблице DATA для возможности взаимообратных преобразований необходимо добавить поле, содержащее новую нумерацию ветвей. В результате стандартных преобразований матриц получим систему дифференциальных уравнений, которая в совокупности с начальными условиями представляет задачу Коши, такую систему можно решать различными методами численного интегрирования. Точность и устойчивость решения будут во многом зависеть от выбранного шага интегрирования, при этом слишком малые значения могут существенно увеличить время расчетов. При формировании начальных условий необходимо учитывать следующие правила: iL (t - 0) = iL (t + 0) и uC (t - 0) = uC (t + 0), где iL — ток через индуктивности в зависимости от

времени t, uC — напряжение на конденсаторе в зависимости от времени t.

Расчет эквивалентных электрических параметров геометрической модели. Расчет эквивалентных значений активных сопротивлений при движении тока в вертикальном направлении (см. рис. 1) для дискретных элементов модели с номерами k и m можно осуществить по формуле:

RmK +k =-2-^Л-' (3)

Ykm 2nrmk mk

Межслойное сопротивление при движении тока в горизонтальном направлении можно рассчитать следующим образом. Рассчитаем элементарное приращение межслойного активного сопротивления для одного (внутреннего) полуслоя:

dR = 1 dr

aRm( K -1)+k+NR = JTi ,

Ymk 2n(rmk + dr)hm

далее запишем в интегральном виде:

гтк +"

Аг.

тк

ёг

Ят (К -1)+к+Ж

тк

Аг

тк

У тк

2лк„

Аг

тк

тк

гтк +

ёг

тк

У тк 2пкт (2гтк + Агтк )

для внешнего слоя на основании аналогичных преобразований получаем

(К-1)+к+Ж

1

Аг

т( к+1)

Ут(к+1) 2пкт (2гт(к +1) +Агт(к+1))

Полное межслойное сопротивление можно вычислить как

1 Агтк , 1

Аг

я

т( к+1)

т( К -1)+к+Ж

-. (4)

Утк 2пкт (2гтк + Агтк ) Ут(к +1) 2пкт (2гт(к +1) + Агт(к +1)) При решении задач на расчет индуктивности [14] принято рассматривать случай постоянного тока и низкой частоты, когда поперечные линейные размеры проводника меньше, или значительно меньше, глубины проникновения тока 5 (1); случай весьма высокой частоты, когда поперечные линейные размеры проводника больше или равны 5. При машинном расчете значения индуктивности (собственные и взаимные) целесообразно хранить в двумерном, квадратном, массиве ШБ размером [(М +1)К]х[(М + 1)К ]; элементы с одинаковыми номерами строк и столбцов будут образовывать множество собственных индуктивностей, остальные элементы будут соответствовать взаимным индуктивностям, причем ШБ[/, -] = ШБ[-,/]. Рассмотрим некоторые примеры, используемые в модели для расчета индуктивностей.

В общем случае собственную индуктивность коаксиального проводника приближенно можно определить как

ШБ[тК + к, тК + к ] = Ц

'тК +к

Д0 К

2п

(1П-

2И

Д

гтк +

А — 1 —1П с),

тк Д0

2

(5)

где ^0 — магнитная постоянная, с — табличная константа [14], зависящая от отношения внутреннего радиуса г

Аг Аг

-"2" к внешнему г + .

В пределах одного горизонтального слоя взаимная индуктивность рассчитывается для коаксиальных проводников в случае низкой частоты:

Г V

М = ^0 2п

1 Ч Д 1П —+ —

Р Д0 ШБ[/, - ] = IND[ -,/] = М

2 2 \ 5 - Ч )

1 5 1П—+-

п

(р2-пг

■1П Р -1

п 2

(

1 +

п

Л

2 2 5 ~ч

2 2 р -п

(6)

Аг,. Аг,.

где п = г,--, р = г, +--

г 2 г 2

Аг,.

Ч = г,

2

Аг,

5 = г, +--

- 2

причем г, > г,, а р = ч .

Расчет взаимной индуктивности между М элементами из различных горизонтальных слоев как между полыми цилиндрическими проводниками с общей осевой линией рассмотрен в [14], в общем случае значение М является функцией от линейных размеров г,, Аг,, к,, г-, Аг-, к- дискретных проводников - и расстоянием между их центрами х:

2

(7)

\М _ /(Г; , А/; , И, Ту, Агу, Иу, X), ШБ [/, у ]_ IND [у,; ]_ М. Рассчитаем межслойную электрическую емкость Стк как совокупность двух последо вательно включенных емкостей

Cmk = 2ns0Smk ■

hm

f rmk + Armk Л

и C'mk = 2ns0 S

hm

0b mk+1 '

ln

mk

ln

mk+1

+ Ar

mk +1

mk+1

С'тк — емкость внутреннего полуслоя, С\"пк — емкость внешнего полуслоя по отношению к площади, разделяющей слои:

С' С''

С _ тк^тк (о)

Стк ~ с' + С" ' ^ '

Стк Стк

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

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

СПИСОК ЛИТЕРАТУРЫ

1. Rycroft M. J. Solar-terrestrial energy program: handbook of ionospheric models // J. Atmospheric and Solar-Terrestrial Physics. Elsevier Science Publishing Company, Inc., 1998. P. 403—404.

2. Guide to Reference and Standard Ionosphere Models. American Institute of Aeronautics and Astronautics, 1999. 55 p.

3. Rasmussen C. E., Schunk R. W. A three-dimensional time-dependent model of the plasmasphere // J. Geophys. Res. 1990. Vol. 95. P. 6133—6144.

4. Смирнов В. М., Смирнова Е. В. Реконструкция пространственно-временной структуры ионосферы по данным спутниковых наблюдений // Современные проблемы дистанционного зондирования Земли из космоса. 2008. Т. 5. С. 561—566.

5. Думин Ю. В. Концентрация носителей заряда в метастабильной плазме с сильной кулоновской неидеальностью // Прикладная физика. 1999. № 5. С. 18—21.

6. Richmond A. D., Ridley E. C., Roble R. G. A thermosphere/ionosphere general circulation model with coupled electrodynamics // Geophys. Res. Lett. 1992. Vol. 19, N 6. P. 601—604.

7. Голыгин В. А., Сажин В. И., Унучков В. Е. Коррекция модели ионосферы по данным о максимально-применимых частотах реперных радиолиний // Исследовано в России. 2006 [Электронный ресурс]: <http://zhurnal.ape.relarn.ru/articles/2006/255.pdf>.

8. Намгаладзе А. А., Юрик Р. Ю. Математическое моделирование возмущений верхней атмосферы Земли [Электронный ресурс]: <http://www.rfbr.ru/pics/28326ref/file.pdf>.

9. Намгаладзе А. А., Мартыненко О. В., Волков М. А., Намгаладзе А. Н., Юрик Р. Ю. Математическое моделирование крупномасштабных возмущений верхней атмосферы Земли // Моделирование процессов в верхней полярной атмосфере. Апатиты: ПГИ КНЦ РАН, 1998. С. 167—249.

10. Гусятинский И. А., Немировский А. С., Соколов А. В., Троицкий В. Н. Дальняя тропосферная радиосвязь. М.: Связь, 1968. 248 с.

Дискретный интегродифференцирующий фильтр

47

11. Хлыбов Е. С., Гаврилов Б. Г., Егоров Д. Е. Модельные исследования ионосферных токов, электрических полей и ускорений заряженных частиц // Сб. тр. XLVII науч. конф. МФТИ. М.: Изд-во МФТИ, 2002.

12. Калинин А. И., Черенкова Е. Л. Распространение радиоволн и работа радиолиний. М.: Связь, 1971. 450 с.

13. Демерчан К. С., Нейман Л. Р., Коровкин Н. В., Чечурин В. Л. Теоретические основы электротехники. СПб: Питер, 2006. 576 с.

14. Калантаров П. Л., Цейтлин Л. А. Расчет индуктивностей. Л.: Энергоатомиздат, 1986. 488 с.

Сведения об авторах

канд. техн. наук; Санкт-Петербургский государственный университет информационных технологий, механики и оптики, кафедра проектирования компьютерных систем; E-mail: tigerpost@ya.ru д-р техн. наук, профессор; Санкт-Петербургский государственный университет информационных технологий, механики и оптики, кафедра проектирования компьютерных систем

Поступила в редакцию 16.03.10 г.

УДК 621.396:681.323

С. И. Зиатдинов

ДИСКРЕТНЫЙ ИНТЕГРОДИФФЕРЕНЦИРУЮЩИЙ ФИЛЬТР С ПРОМЕЖУТОЧНЫМ СГЛАЖИВАНИЕМ ОТСЧЕТОВ ВХОДНОГО СИГНАЛА

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

Ключевые слова: дискретизация сигнала, промежуточное суммирование, дискретная фильтрация, отношение сигнал/шум.

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

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

Суть предлагаемого метода фильтрации заключается в том, что из последовательности отсчетов входного сигнала х[1] формируются промежуточные суммы из т отсчетов

т-1

хЕ [п] _ £ х[пЬ -1 ],

;_0

Алексей Юрьевич Гришенцев Анатолий Григорьевич Коробейников

Рекомендована кафедрой проектирования компьютерных систем