ВЫЧИСЛЕНИЕ МНОЖИТЕЛЯ ОСЛАБЛЕНИЯ НАД ЗЕМНОЙ ПОВЕРХНОСТЬЮ МЕТОДОМ ПАРАБОЛИЧЕСКОГО УРАВНЕНИЯ

В. В. Ахияров

МГТУ им. Н.Э. Баумана

Получена 26 января 2011 г.

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

Ключевые слова: распространение радиоволн, параболическое уравнение.

Abstract. The solution of a parabolic equation for prediction of the field strength above the earth, allowing for the vertical atmosphere stratification and the geometry of terrain, is considered. The numerical algorithm based on split-step Fourier transform is shown. Solutions for well-known diffraction problems such as spherical earth and a wedge are obtained. The results for real terrain profiles are produced.

Keywords: radio wave propagation, parabolic equation.

Метод параболического уравнения (ПУ) широко используется в акустике, электродинамике и оптике для расчета волновых полей в неоднородных средах. Идея метода принадлежит М.А. Леонтовичу, который преобразовал исходную краевую задачу для уравнения Гельмгольца (эллиптического типа) к параболическому уравнению. Далее В.А. Фок получил аналитические решения ПУ над сферической Землей для линейной модели атмосферы и приземного волноводного канала [1]. В настоящее время существует возможность решения ПУ численными методами для произвольного высотного профиля показателя преломления атмосферы и любой геометрии рельефа [2]. Рассматриваемая задача считается скалярной (см. рис.1), т.е. при горизонтальной поляризации излучения отличной от нуля является поперечная компонента электрического поля:

а при вертикальной поляризации решение ищется для поперечной компоненты магнитного поля:

Функция должна удовлетворять уравнению Гельмгольца:

, (1)

где – волновое число, - высотный профиль показателя преломления атмосферы.

Рис.1. Геометрия задачи.

Полагая зависимость поля от времени в виде , выделим в быстро осциллирующий множитель:

и получим уравнение Гельмгольца для амплитуды поля:

. (2)

Далее представим (2) в виде произведения двух дифференциальных операторов [2]:

, (3)

где:

. (4)

В (3) первый множитель дает уравнение для волны, распространяющейся в направлении оси Х, а второй – для волны, имеющей противоположное направление.

Дифференциальный оператор (4) можно аппроксимировать различными способами, однако проще всего ограничиться первыми двумя членами ряда Тейлора:

. (5)

Подстановка (5) в уравнение

(6)

приводит к ПУ для распространения в направлении оси Х:

. (7)

Выражение (5) справедливо в малоугловом (параксиальном) приближении, которое выбирается из условия (см. рис.1), поэтому ПУ (7) является малоугловым приближением уравнения (6). Другим способом аппроксимации оператора (4) является выражение [3]:

, (8)

где соответствует оператору (4) для свободного пространства. После подстановки правой части (8) в уравнение (6) получим:

. (9)

Для численного решения ПУ (7) и (9) может использоваться либо конечно-разностная методика, либо подход, основанный на вычислении прямого и обратного преобразования Фурье:

, (10.а)

, (10.б)

где , а максимальная высота Z_max связана со значением и размером преобразования Фурье L критерием Найквиста [5]:

, (11)

который позволяет определить шаг по высоте и дальности .

Метод решения (7) и (9), основанный на паре преобразований Фурье (10.а, б), заключается в следующем: на дальности x поле разлагается в угловой спектр плоских волн и умножается на передаточную функцию слоя пространства. Далее вычисляется обратное преобразование Фурье, соответствующее распределению поля по высоте на дальности и результат умножается на дополнительный фазовый множитель, учитывающий рефракцию радиоволн в атмосфере Земли. Передаточная функция определяется выражением:

для уравнения (7) и

для уравнения (9).

Таким образом, на каждом шаге численного решения необходимо вычислять выражение:

, (12.а)

при решении (7) или

, (12.б)

при решении (9).

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

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

Для учета кривизны Земли в уравнениях (7) и (9) необходимо выполнить замену модифицированным показателем преломления , где a – радиус Земли. Тогда с учетом выражения

получим решение ПУ (9)

(13)

для расчета напряженности поля над сферической Землей.

Поскольку мало отличается от единицы, на практике используется модуль приведенного показателя преломления или индекс рефракции:

. (14)

Для того, чтобы учесть влияние рельефа на напряженность поля, используются методы ступенчатого моделирования, кусочно-линейной аппроксимации и конформного преобразования [2]. В данной работе используется метод ступенчатого моделирования, при этом трасса распространения аппроксимируется вертикальными «ступеньками», как это показано на рис.2. При распространении над горизонтальными участками S₁ поле вычисляется с использованием выражений (12.а, б), а на ступеньках S₂ результаты численного решения приравниваются к нулю (см. рис.3).

Рис.2. Моделирование геометрии рельефа.

Рис.3. Алгоритм вычисления амплитуды поля.

При реализации вычислительных алгоритмов, основанных на преобразованиях Фурье (10.а, б), используется метод зеркальных изображений, т.е. область (см. рис.1) отображается в нижнюю полуплоскость относительно оси X. Ограничение области расчетов приводит к отражению волн от искусственных границ на высоте , поэтому необходимо подавление поля, выходящего за пределы области расчетов. Для алгоритмов, основанных на прямом и обратном преобразовании Фурье по поперечной координате, подавление искусственных отражений осуществляется с помощью поглощающего или идеально согласованного слоя, при этом решение ищется в полосе , где .

Поглощающий слой представляет собой пространственный фильтр, реализованный на основе весового окна Хэмминга аналогично тому, как это делается при спектральной обработке сигналов. Идеально согласованный слой обычно рассматривается как среда с потерями, в которой амплитуда поля экспоненциально убывает. Таким образом, с точки зрения реализации вычислительных алгоритмов, ограничение области расчетов поглощающим или идеально согласованным слоем эквивалентно применению весового окна [6]. Следует отметить, что наличие рельефа и сферичность Земли не влияют на алгоритм подавления поля в области .

С использованием рассмотренных алгоритмов были получены решения для сферической Земли, линейной модели атмосферы и приземного волноводного канала. Источник излучения моделировался линейной апертурой размером , что соответствует шагу по вертикальной координате. На рис.4 представлены результаты расчетов множителя ослабления для линейной модели атмосферы с градиентом показателя преломления , что соответствует эквивалентному радиусу Земли . Исходные данные выбраны следующими: длина волны , высота источника , горизонтальная координата x является расстоянием по дуге сферической земли. Сравнение расчетов методом ПУ для граничного условия Дирихле на земной поверхности с аналитическим решением по дифракционной формуле В.А. Фока [1] для Земли с реальными электрическими параметрами показало их полное соответствие [6]. Поэтому при численном решении данной задачи представление о граничных условиях Дирихле вполне оправдано.

Рис.4. Множитель ослабления над сферической земной поверхностью.

Далее были выполнены расчеты для модели приземного волноводного канала, который определяется безразмерной функцией , связанной с индексом рефракции [1]:

где .

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

, (15)

где y_i – приведенная высота точки инверсии, y_l – параметр, , (z – высота над Землей).

На рис.5 представлен высотный профиль индекса рефракции для следующих исходных данных [1]: , , , на рис.6 –значения множителя ослабления в диапазоне высот от 0 до 75 м на дальности до 500 км. Высоты источника выбраны равными (рис.6.а) и (рис.6.б), значение соответствует высоте точки инверсии на профиле индекса рефракции. Из приведенных рисунков видно, что при высоте источника 9,31 м внутри слоя инверсии наблюдается большая напряженность поля, которая обусловлена более высокой кривизной профиля на малых высотах. Необходимо отметить, что полученные результаты полностью соответствуют аналитическому решению данной задачи [7].

Рис.5. Высотный профиль индекса рефракции.

Рассмотрим результаты решения задачи дифракции для клина, синусоидальной поверхности и реальных профилей рельефа с использованием алгоритма ступенчатого моделирования (см. рис.2). Как и в предыдущем случае, размер апертуры источника излучения равен .

а)

б)

Рис.6. Множитель ослабления в приземном тропосферном волноводе:

а – ; б – .

На рис.7 представлены результаты расчетов множителя ослабления над клином длиной 10 км и высотой 200 м при и высоте источника . На рис.8 показано сравнение численного и аналитического результатов решения данной задачи при тех же исходных данных и высоте точки наблюдения, равной высоте источника. Сплошная линия на данном рисунке получена в результате численного решения ПУ, точки – аналитический расчет с использованием интерференционной формулы в освещенной области и геометрической теории дифракции в зоне тени. Рис.9 иллюстрирует результаты решения ПУ для синусоидального профиля трассы протяженностью 10 км и высоте препятствий, равной высоте источника при .

Рис.7. Дифракция на клине.

Рис.8. Сравнение результатов решения ПУ (сплошная линия)

с аналитическими расчетами (точки).

Рис.9. Дифракция на синусоидальном профиле.

Рассмотрим результаты расчетов для реальных профилей земной поверхности. На рис.10 показана геометрия рельефа, на рис.11 представлены результаты расчетов множителя ослабления при следующих исходных данных: длина волны , высота источника равна высоте точки наблюдения – . Сплошная линия на данном рисунке получена путем решения интегрального уравнения [8], точки соответствуют параболическому уравнению. Поскольку в методе интегральных уравнений изменение показателя преломления с высотой не учитывается, исходный профиль был изогнут по дуге радиуса , что соответствует нормальным условиям рефракции.

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

Рис.10. Профиль рельефа.

Рис.11. Дистанционная зависимость множителя ослабления.

Сплошная кривая – решение интегрального уравнения,

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

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

Рис.12. Множитель ослабления над земной поверхностью.

Следует отметить, что в вычислительном отношении рассмотренные алгоритмы являются очень эффективными: использование ноутбука с процессором Core 2 Quad Q9000 (тактовая частота 2 ГГц) и оперативной памятью объемом 4 ГБ позволяет получить любой из представленных в работе результатов за время, не превышающее нескольких минут. Ключевым моментом столь высокой скорости расчетов является использование критерия Найквиста (5), который дает возможность выбрать интервал дискретизации существенно больше длины волны. Например, для малоуглового приближения и размера преобразования Фурье при длине волны получим интервал дискретизации , что соответствует разбиению трассы длиной 10 км на интервалов. С использованием данных параметров представленные на рис.7 и рис.9 результаты можно получить менее чем за одну минуту.

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

Литература

1. Фок В.А. Проблемы дифракции и распространения электромагнитных волн. М.: Сов. радио. 1970.

2. Levy М.F. Parabolic equation method for electromagnetic wave propagation. London. IEE. 2000. 336 p.

3. Levy М.F. Perfectly matched layer truncation for parabolic wave equation models // Proc. R. Soc. Lond. A. 2001. pp. 2609-2624.

4. Иванов В.К., Шаляпин В.Н., Левадный Ю.В. Рассеяние ультракоротких волн на тропосферных флуктуациях в приводном волноводе // Известия вузов. Радиофизика. 2009. т.LII. №4. С.307-317.

5. Sevgi L., Uluisik C., Akleman F. A MATLAB-based two-dimensional parabolic equation radiowave propagation package. – IEEE Antennas and Propagation magazine, 2005, vol. 47, no.4, pp.164-175.

6. Ахияров В.В. Метод параболического уравнения в теории дифракции. – Успехи современной радиоэлектроники. 2010. №9. c.72-80.

7. Ахияров В.В., Чернавский С.В. Использование численных методов для изучения условий распространения радиоволн // Радиотехника. 2011. №10. с.101-110.

8. Ахияров В.В. Методы численного решения задачи дифракции радиоволн над земной поверхностью // Электромагнитные волны и электронные системы. 2010. т.15. №3. С.39-46.