КОНФОРМНОЕ ОТОБРАЖЕНИЕ КАК МЕТОД УЧёТА НЕРОВНОСТЕЙ ПОДСТИЛАЮЩЕЙ ПОВЕРХНОСТИ С МАЛОЙ КРУТИЗНОЙ ПРИ ЧИСЛЕННОМ РЕШЕНИИ ПАРАБОЛИЧЕСКОГО УРАВНЕНИЯ Текст научной статьи по специальности «Математика»

УДК 537.874.3

А.В. Новиков, Ю.П. Акулиничев

Конформное отображение как метод учёта неровностей

подстилающей поверхности с малой крутизной

при численном решении параболического уравнения

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

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

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

В настоящее время имеется потребность в программных продуктах, позволяющих численно оценивать зоны радиопокрытия для заданных условий, ключевыми из которых являются рельеф, проводящие и диэлектрические свойства подстилающей (земной) поверхности. Сложность решения задачи распространения радиоволн (РРВ) связана с произвольностью рельефа. Одним из простейших способов представления рельефа поверхности является ступенчатая аппроксимация («stair case» [1]). Вторым шагом к уточнению описания рельефа является его кусочно-линейная аппроксимация («shift map» [1]). Две перечисленные методики представления рельефа содержат не вполне адекватное описание его профиля ввиду отсутствия кривизны у линейно'-ломаных функций (она равна либо нулю, либо бесконечности). Тем не менее эти методы дают качественно схожие результаты (см. разд. «Результаты расчетов») и им нет равных по вычислительным затратам и простоте учёта рельефа.

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

В итоге задача РРВ в однородной среде над неровной земной поверхностью преобразуется в задачу РРВ над плоской поверхностью в неоднородной среде (атмосфере). Это возможно благодаря свойству конформных отображений сохранять форму оператора Лапласа, который в новой системе координат будет с точностью до якобиана преобразования, записан так же, как и в старой. Эта идея не нова [3], однако результатов сравнительных расчётов применительно к распространению радиоволн в неоднородной среде в литературе нами пока не найдено. Исключение составляет упоминание Б.Д. Тартаковским в [4] о работе В.И. Аксёнова, в которой автор «предложил приближённый метод расчета амплитуды электромагнитной волны, зеркально отражённой от периодически неровной диэлектрической поверхности, который основан на замене такой поверхности «эквивалентной» слоисто-неоднородной средой». Эта работа ограничена случаем, когда период неровностей меньше длины волны.

В предлагаемом методе требуется представление профиля рельефа усеченным рядом Фурье.

Метод конформного отображения для решения двумерного параболического уравнения

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

Считаем также, что источник поля является линейным и расположен параллельно оси Oz. Тогда волны, распространяющиеся в такой среде, являются цилиндрическими, и для их описания достаточно использовать двумерное уравнение Гельмгольца для напряжённости поля E(x, у):

д2 Е д2 Е 2

—-+—•+к2е(х, у)Е = 0, (1)

дх2 дУ2

где к = 2п / X - волновое число.

Для решения задачи распространения радиоволн вместо уравнения Гельмгольца используют ПУ, которое следует из уравнения (1) после введения функции

и(х,у) = Е(х,у)ехр(-гкх), (2)

которая является медленно меняющейся функцией по х по сравнению с функцией электрического (или магнитного) поля Е(х, у).

Каждую точку с координатами (х, у) можно представить в виде комплексного числа г = х + гу. Отображение старой системы координат г в новую систему w = и + IV задаётся некоторой функцией w = /(г). Преобразование координат /(г) является конформным, если функция /(г) удовлетворяет уравнению Лапласа [5]

д2 //дх2 + д2//ду2 = 0. (3)

Примером конформного отображения является функция г'а-1п(г / га), которая отображает кольцо

а2 < х2 + у2 < (а + й)2

в прямоугольник

-па < и < па, 0 < V < й.

Такое логарифмическое отображение используется при аналитическом решении ПУ для случая

распространения радиоволн над полупроводящей Землёй с рельефом, имеющим постоянную кри-

визну [2]. В этом случае используется модифицированный индекс преломления М, в котором, помимо информации о рельефе, содержится информация о параметрах высотного профиля показателя преломления стандартной атмосферы.

Цель данной статьи - отыскать подходящее конформное отображение для произвольных неровностей земной поверхности, таких, профиль которых может задаваться набором отсчётов высоты.

Представим некоторое отображение / (г) в виде суммы периодических функций комплексной переменной N

/ (г)=г + Е /п(г)=и (х у)+Н х у), (4)

п=1

где /п (г) = -гап ехр(гКпг); Кп = 2пп / Б; Б - длина трассы (дальность); Б / п - пространственный период (вдоль оси Ох) функции/п (г),

N

и( х, у) = х + Е Ап ехр(-Кпу)т(Кпх+фп ), (5)

п=1

N

у)=у -Е Ап ехр (—Кп у) (Кпх+фп). (6)

п=1

Количество слагаемых N и их весовые коэффициенты ап = Ап ехр(гфп) подбираются так, чтобы кривая V = 0 с необходимой точностью совпадала с функцией §(х). Каждая из функций (5) и (6) удовлетворяет уравнению Лапласа (3), поэтому функция (4) является аналитической, а найденное преобразование координат - конформным.

Для упрощения дальнейшего анализа следует заметить, что ПУ является малоугловым приближением к уравнению Гельмгольца, поэтому максимальные наклоны земной поверхности у(х) при V = 0 должны быть малы

\йуЩтах = V(1 - 0 <<1, 1 = ЕАпКп . (7)

п

Найдем значения коэффициентов Ап и фп. Из выражения (6) следует, что уравнению V = 0 соответствует уравнение

Рис. 1. Неровная земная поверхность g(x) в декартовой системе координат

N

у =Е Ап ехр(-Кпу)е08(Кпх + фп). (8)

п=1

Первым приближением к решению уравнения (8) относительно у(х) будет функция

N

у(х) =Е Ап 008((пх + Фп), (9)

п=1

что справедливо при выполнении условия

? <<1, (10)

где следует напомнить, что числа Ап и Кп принадлежат множеству действительных неотрицательных чисел. Условие (10) не противоречит условию (7), т.к. оно слабее, поэтому для случая малых

углов наклона земной поверхности коэффициенты Ап и фп определяются из разложения функции ё(х) в ряд Фурье

N

ё(х) =Е Ап 008((пх + Фп). (11)

п=1

Функция ё(х) должна иметь нулевой средний уровень. Итак, при заданном массиве высот земной поверхности g = {ёг}, г = 0...2М-1, необходимо после центрирования вектора g методом прямого дискретного преобразования Фурье найти вектор с и вычислить значения коэффициентов

Ап = 2|сп|, Фп = аг§(сп), п = 1..М . (12)

Чем больше число учитываемых гармоник ряда Фурье, тем выше точность представления рельефа поверхности, но требуется больше времени.

И, наконец, запишем выражение для якобиана преобразования (4)

3 (х, у)=(|и) +(£) = (1+ис )2 + иД (13)

N N

где ис = Е АпКп ехр (- Кп у) 008 {Кп х + Фп ); из = Е АпКп ехр(-Кпу)81п (Кпх + Фп ).

п=1 п=1

При выполнении условия (7) якобиан (13) можно записать в более простом виде

3 (х, у) »1 + 2ис (14)

2

с ошибкой порядка 2? ^ 1.

После замены переменных в уравнении Гельмгольца (1) в соответствии с (5) и (6) имеем

Л-+а2Е+к21М_Е=0. (15)

ди2 ду2 3 (х у)

Переходя в (15) от уравнения Гельмгольца к ПУ, получим

2/кди+д^ + к2Г1ий-11и = 0. (16)

ди дv2 13 (х, у) )

Уравнение (16) записано в прямоугольной системе координат (и, V), поэтому его можно численно решать обычными для него способами (схема Кранка-Николсон, метод преобразования Фурье) [1].

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

ди / дv + аи = 0,

где коэффициент а равен

[г'к%/ е-1 - для горизонтальной, гк/ л/е-1 - для вертикальной поляризации.

а =

\гк

Результаты расчётов

Результаты сравнительного расчёта высотного профиля множителя ослабления для двух трасс (рис. 2) показаны на рис. 3.

Рис. 2. Треугольный и синусоидальный профили трасс РРВ

Здесь сокращение «Conf. map» означает метод конформного преобразования, «app. Jacob.» -использовалось приближённое вычисление якобиана J по (14), а «exact. Jacob.» - точное выражение (13) для J. «Shift map» - метод пошагового преобразования системы координат с кусочно-линейной аппроксимацией профиля рельефа. Длина волны - 10 см, высота антенны излучателя - 2 м, диаграмма направленности - гауссовская, шириной 16° (чтобы смоделировать всенаправленность для высот 0...150 м), длина трассы - 1 км, высота области расчёта - 150 м. Холм треугольной формы высотой 5 м и с основанием в 640 м расположен по центру трассы. Поляризация горизонтальная, коэффициент отражения - минус 1. Число гармоник для аппроксимации треугольного рельефа - 64 (СКО ошибки - около 2-3 см). Расчёт вёлся по схеме Кранка-Николсон с шагом сетки 8,8 см по высоте и 21,5 см по дальности (1710x4640 отсчётов). Влияние искусственно введенной диэлектрической проницаемости учитывалось применением метода расщепления [1].

- 30

-20

- 10

-20

- 10

Множитель ослабления, дБ

---Conf. map (app. Jacob.) ----------Conf. map (exact. Jacob.) ♦♦♦ Shift map

Рис. 4. Высотные профили множителя ослабления: слева - для профиля трассы треугольной формы,

справа - для синусоидальной

Результаты расчёта поля по методу «stair case» на рис. 3 не отображены из-за небольшой разницы (менее 1 дБ) по сравнению с результатами для метода «shift map».

Результаты оценки времени расчёта для всех кривых на рис. 3 приведены в таблице.

Время расчёта профилей множителя ослабления, с,

Профиль поверхности: Exact. Jacob. App. Jacob. Shift map

Синусоидальной формы 4 3 2

Треугольной формы 51 36 1

Заключение. Два профиля земной поверхности (см. рис. 2) имеют одинаковую максимальную крутизну (около 0,9°). Различия в профилях поля (см. рис. 3) при использовании точного (13) и приближённого (14) выражений для вычисления якобиана менее заметны для синусоидальной поверхности, для которой малый параметр t равен 0,016 против 0,05 для трассы с треугольным профилем. Также при аппроксимации синусоидального профиля рельефа использовались всего два члена ряда

Фурье против 64, что позволило сократить время расчёта на порядок (см. таблицу). Различия двух кривых («shift map» и «conf. map (exact. Jacob.)») по множителю ослабления на рис. 3 для двух трасс не превышают 2 дБ в области полутени (менее 20 м по высоте). Большие различия в минимумах (10 - 30 дБ), видимо, обусловлены тем, что метод конформного преобразования более корректно учитывает кривизну рельефа (правда, из-за математических трудностей решения (8) относительно y(x), это справедливо только для случая малых углов). Методы «stair case» и «shift map» не ограничены углами рельефа (исключая углы, близкие к прямому), правда, из-за ограничений ПУ это достоинство нивелируется. Хотя для широкоугольных вариантов ПУ есть шанс использовать, например, метод «shift map». При этом, однако, необходимо уменьшить шаг сетки по высоте в соответствии с теоремой отсчётов.

По затратам времени более предпочтителен метод «shift map». Достоинством метода конформного отображения является возможность сравнения влияния неровностей рельефа и неоднородностей тропосферы на уровень поля путём их пересчёта в неоднородности показателя преломления. Показано, что два метода («shift map» и «конформное отображение») дают достаточно близкие результаты.

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

Литература

1. Mireille Levy. Parabolic equation methods for electromagnetic wave propagation / The Institution of Electrical Engineers. - London: United Kingdom, 2000. - 336 p.

2. Фок В. А. Проблемы дифракции и распространения радиоволн. - М.: Сов. радио, 1970. - 520 с.

3. Долуханов М.П. Распространение радиоволн. - М.: Связь, 1972. - 336 с.

4. Тартаковский Б. Д. Второй Всесоюзный симпозиум по дифракции волн // Успехи физических наук. - 1962. - Т. LXXVIII, вып. 2. - С. 701-721.

5. Лаврентьев М.А. Методы теории функций комплексного переменного / М.А. Лаврентьев, Б.В. Шабат. - М.: Наука, 1965. - 716 с.

Новиков Анатолий Викторович

Ассистент каф. радиотехнических систем ТУСУРа

Тел.: 8-952-155-96-99

Эл. почта: anovikov@ms.tusur.ru

Акулиничев Юрий Павлович

Д-р техн. наук, профессор каф. радиотехнических систем ТУСУРа

Тел.: 8-913-806-80-01

Эл. почта: ayp@rts.rk.tusur.ru

Novikov A.V., Akulinichev Yu.P.

Conformal mapping as a method for the numerical solution of the parabolic wave equation in the presence of small-angle terrain irregularities

Conformal mapping which allows remapping of terrain irregularities into equivalent refractive index of the troposphere has been found. Comparative results of the numerical calculation of the propagation factor height profile are given for two surface models using two methods. Comparison was made with the method of piecewise-linear approximation of the terrain profile. It is shown that the conformal mapping method is suitable only for terrain profiles well approximated by a truncated Fourier series. In the absence of exact solutions for the considered paths (the triangular and sinusoidal hills) it should be noted only that both methods give the propagation factor height the profiles which are qualitatively similar (the difference is not more than 2 dB in penumbra and the shadow region).

Keywords: refractive index, conformal mapping, rough ground surface, radio propagation, parabolic wave equation.