Limiting accuracy of two-dimensional parabolic equation solution using split-step Fourier transform method.pdf Введение В настоящее время эффективность любой радиотехнической системы (РТС) в существенной степени зависит от условий распространения радиоволн (РРВ) в неоднородной среде. Поэтому развитие методов прогноза ожидаемых характеристик поля и соответствующей эффективности РТС является одной из актуальных задач в областях радиолокации, радионавигации и беспроводной связи. В данной статье не рассматриваются довольно простые, но недостаточно точные [1] эмпирические модели прогнозирования, такие, как модели Longley-Rice [2] и Hata-Davidson [3], рекомендация МСЭ-R P.1546 и др. Более перспективны детерминированные модели, позволяющие имитировать процессы РРВ практически в любом диапазоне радиоволн на конкретной трассе и учитывать рефракцию и дифракцию, отражения от препятствий на трассе и изменения коэффициента преломления среды в пространстве. Естественно, такие модели сложнее эмпирических и на их применение требуется значительно больше вычислительных затрат. Однако с развитием техники эта проблема постепенно теряет актуальность [4]. Наиболее популярной математической детерминированной моделью РРВ является параболическое уравнение (ПУ). Впервые оно было введено в работах Леонтовича и Фока по решению задач РРВ над сферической поверхностью Земли [5, 6]. Широкое распространение использования ПУ началось после того, как Hardin и Tappert импортировали преобразование Фурье с расщеплением из алгоритмов, применяемых при подводных акустических исследованиях [7], а Claerbout разработал конечно-разностную схему для геофизических приложений [8]. С тех пор способы применения ПУ в качестве модели атмосферного РРВ регулярно совершенствуются [4, 9] и используются на практике, например в приложениях прогнозирования радиолокационного покрытия [10]. Скалярное ПУ также широко применяется в сферах гидроакустики, оптики, сейсмологии [11] и др. Если волна является монохроматической, т.е. зависимость от времени имеет вид exp(-iωt), то распространяющееся в декартовом пространстве xOz электромагнитное поле заданной поляризации ψ(x, z) можно описать двумерным скалярным уравнением Гельмгольца. Если при этом поле распространяется преимущественно в направлении оси Ox+, то в поперечном направлении Oz оно изменяется значительно медленнее, чем вдоль трассы, поэтому удобно использовать представление ψ(x, z) = u(x, z) exp(ikx), где u(x, z) - медленно изменяющаяся в пространстве комплексная огибающая, k = 2π/λ - волновое число, λ - длина волны. Тогда уравнение Гельмгольца можно записать в виде [11] , (1) где , ε(x, z) - относительная диэлектрическая проницаемость среды. В практических задачах вторым множителем в (1) пренебрегают и используют первый, так как именно он описывает распространение поля в положительном направлении оси Ox. Существует несколько способов аппроксимации оператора Q. Вид выбранной аппроксимации определяет максимально допустимое значение βmax модуля угла β между вектором Пойтинга и осью Ox+, при котором ПУ (1) приемлемо. Для моделей Claerbout [12], Feit-Fleck [13] и Padé [14] значение этого угла может достигать 30-45. Однако наиболее распространена самая простая модель малоуглового (стандартного) ПУ (βmax < 10-15) . (2) Детальный вывод уравнения (2) приведен в [11] путем аппроксимации оператора Q первыми двумя членами ряда Тейлора по β и в [15] путем аппроксимации самого уравнения Гельмгольца с учетом того, что u(x, z) изменяется по переменной x значительно медленнее чем по z. В практических задачах используют численные методы решения ПУ, которые предполагают наложение на область расчета дискретной сетки с шагами дискретизации по дальности ∆x и поперек трассы (например, по высоте) ∆z. Существующие сегодня сеточные методы можно разделить на три группы: методы дискретного преобразования Фурье с расщеплением (ДПФ), методы конечных разностей (МКР) и конечных элементов (МКЭ) [16-19]. Несмотря на то, что для последних двух разработаны практически идеальные граничные условия (ГУ) [4], при прочих равных условиях их точность чаще всего меньше точности ДПФ [4]. К тому же при увеличении частоты излучения МКР и МКЭ требуют уменьшать шаги дискретизации расчетной сетки [20], что значительно повышает объем проводимых вычислений [4, 16]. Согласно методике ДПФ с расщеплением значения поля определяются лишь в узлах прямоугольной сетки с ячейками Δx•и Δz последовательно шаг за шагом по мере удаления от источника, т.е. находят решение дискретного аналога уравнения (2) , (3) где u(x) и u(x+∆x) - векторы-столбцы значений комплексной огибающей напряженности поля в узлах сетки при x и x+∆x соответственно; F и F-1 = FH - квадратные матрицы прямого и обратного ДПФ; K(x) - диагональная матрица коэффициентов передачи гармоник ряда Фурье на одном шаге по правилам РРВ в свободном пространстве (при ε(x, z) = 1) в слое (x, x+∆x); L(x) - диагональная матрица, вводимая для учета нижнего и верхнего граничных условий и называемая поглощающим слоем (ПС); E(x) - диагональная матрица с элементами Enn = 0.5 k ∆x[ε(x, ∆z∙n) - 1], введенная в соответствии с методом расщепления для учета неоднородностей среды на данном шаге. Уже несколько десятилетий используется упрощенный вариант (3), а именно считается, что матрицы L и K не зависят от расстояния x [21]. Более того, элементы матрицы K формально определяются из непрерывного ПУ (2), описывающего РРВ в среде, безграничной по оси Oz, , (4) хотя требование теоремы отсчетов при этом учитывается. Более корректно было бы задавать их путем сравнения дискретных спектров S(; x) = F u(x) реализаций комплексной огибающей напряженности поля протяженностью N•∆z в начале и в конце шага. При РРВ в тропосфере вдоль поверхности земли нижнее ГУ обычно вводится в виде смешанного ДПФ, которое позволяет успешно учитывать электрические свойства подстилающей поверхности трассы РРВ [22, 23]. Попытки же введения ГУ на границах со свободным пространством наталкиваются на существенные трудности, поскольку метод ДПФ уже обладает собственным граничным условием - это периодичность функции, описываемой рядом Фурье. В итоге и амплитудная, и фазовая характеристики такого пространственного фильтра при больших значениях угла β слабо удовлетворяют исходному уравнению (2) [24]. Поэтому пока единственным простым способом введения верхнего ГУ является применение поглощающих слоев (ПС). В их качестве, основываясь только на интуитивных предположениях, обычно используют степенные оконные функции, окна Хэннинга или Хэмминга [25]. Тем не менее, метод ДПФ до настоящего времени считается эталонным для поверки других приближенных методов [4, 9]. Причина состоит в том, что указанные недостатки уменьшаются до любой желаемой степени простым увеличением поперечного размера N•∆z области расчета. Что, разумеется, приводит к увеличению вычислительных затрат и снижению оперативности прогнозирования. Вероятно, этим и объясняется тот факт, что в доступной литературе не удалось найти результатов исследования потенциальной точности метода ДПФ при ограничении объема вычислений, тем более, что весьма сомнительно существование аналитических методов решения этой задачи. Цель данной работы - определение предела для среднеквадратической ошибки (СКО) расчета напряженности поля в однородной безграничной среде при использовании метода (3) без увеличения вычислительных затрат, т.е. фактически определение минимальной СКО расчета элементов дискретной функции Грина интегрального уравнения, соответствующего уравнению (2) [11]. По необходимости при этом также решается задача поиска оптимальных значений матриц K(x) и L(x) в (3) и исследуется характер их зависимости от геометрических параметров вычислительной схемы, в первую очередь, от расстояния x до источника (первичного или вторичного, для интегрального ПУ это не имеет значения). Результаты исследования и их обсуждение Решение задачи проводилось при следующих предположениях. 1. Рассматривалось РРВ в неограниченной среде, поэтому влияние подстилающей поверхности не учитывалось. Расчетная область, содержащая N ячеек по оси Oz, делилась на три части: нижнюю и верхнюю половины ПС, состоящая каждая из NL ячеек, и среднюю (рабочую) область из NW = N - 2NL ячеек, в которой поглощение отсутствует. 2. Поскольку область, занятая ПС, играет вспомогательную роль, то источник излучения и все неоднородности среды, которые формально являются вторичными источниками, располагались лишь в рабочей области. Только в этой области вычислялась и СКО расчета поля. 3. С методической точки зрения удобно было считать, что отсчеты первичного поля при x = 0 есть совокупность NW независимых гауссовских комплексных случайных величин с нулевым математическим ожиданием и единичной дисперсией. 4. Поскольку речь идет о потенциальной точности, сделано предположение, что на предыдущих шагах ошибки не возникали, поэтому эталонные поля в начале и в конце шага задавались по результатам тестового расчета. Для теста применялся тот же метод ДПФ с теми же источниками, но поперечная ширина области расчета выбрана достаточно большой N•∆z >> β•x так, чтобы значения поля на границах были практически равны нулю, что позволяет всегда использовать граничное условие Дирихле. 5. Оптимальные значения K и L находили методом поочередных последовательных приближений, при этом в качестве нулевого приближения использованы те значения, которые были получены на предыдущем шаге по х. Поскольку известно, что при x → 0 матрица K определяется формулой (4), а при x → ∞ имеем L ≈ I, где I - единичная матрица, то процедура расчета для контроля проводилась дважды: при увеличении и при уменьшении расстояния х. 6. Для поиска наилучших значений элементов матриц K и L на каждой итерации применялась стандартная процедура линейного оценивания по критерию минимума СКО [26]. Чтобы не повторять однотипных расчетов, для расстояния от источника и шага по дальности используем нормировку, основанную на теореме отсчетов, , , (5) причем во всех приводимых результатах принято ∆z βmax/ λ = 0.5. Итоговое значение относительной СКО расчета поля σ на дальности m ∆x определялось как , (6) где v(m ∆z) и u(m ∆z) - результаты тестового и проведенного расчетов соответственно. Зависимости величины СКО σ от нормированного расстояния до источника излучения xн при различных размерах рабочей области N ∆z и ширины поглощающих слоев NL ∆z представлены на рис. 1. Видно, что величина СКО σ с ростом xн увеличивается до максимального значения , (7) которое достигается на расстоянии x0 ≈ 0.5 (N - 2 NL) ∆xн, а последующее уменьшение значений можно аппроксимировать гиперболой . (8) Такой характер зависимости можно объяснить тем, что до расстояний x0 сферичность фазовых фронтов распространяющихся волн выражена еще довольно ярко, а на больших расстояниях любая сферическая волна преобразуется в отрезок плоской, поэтому представление поля в виде фиксированного числа плоских волн становится проще и точнее (рис. 1). Рис. 1. Зависимость величины относительной СКО σ от нормированного расстояния при ∆xн = 1, числе отсчетов поля по высоте N = 128 (а) и 256 (б), числе отсчетов ПС: кр. 1 - NL = N/16; кр. 2 - NL = N/8; кр. 3 - NL = N/4; кр. 4 - соответствующие аппроксимации (8) Характер зависимостей коэффициентов передачи гармоник (КПГ) ряда Фурье от параметров расчетной схемы следующий. До расстояний x0 аргументы КПГ определяются традиционной формулой (4), однако их модули с увеличением уменьшаются линейно, что полностью согласуется с результатами работ [24, 27]. На больших расстояниях xн >> x0 матрица K(x) близка к единичной, так как сферические волны от всех источников в дальней зоне вырождаются в одну плоскую. Рис. 2. Формы модулей (а) и аргументов (б) поглощающих слоев при N = 128, NL = 16, ∆xн = 1 на нормированных расстояниях от источника излучения: кр. 1 - xн = 64; кр. 2 - xн = 128; кр. 3 - xн = 256; кр. 4 - окно Хэмминга Получаемые в результате работы вышеописанного алгоритма поглощающие слои (рис. 2), во-первых, являются комплексными. Это существенное отличие полученных ПС от традиционно применяемых чисто действительных ПС в виде окон Хэннинга или Хэмминга (рис. 2, кривые 4). Во-вторых, формы оптимальных ПС изменяются с расстоянием от источника излучения x. При этом вносимое ими ослабление поля с ростом x уменьшается, что также объяснимо уменьшением сферичности фронтов распространяемых радиоволн с отдалением от источника излучения. Величина СКО σ расчета поля, естественно, зависит от величины шага дискретизации по дальности ∆xн. При малых ∆xн эта зависимость практически линейна. Так, на рис. 3 кривые 2, 4 с достоверностью 0.999 аппроксимируются линейными функциями и соответственно, т.е. ошибки на соседних шагах почти полностью коррелированы. Поэтому не следует делать поспешный вывод о том, что с уменьшением шага Δxн можно добиться сколь угодно малой величины СКО, поскольку на практике при этом возрастает количество шагов и, следовательно, количество накапливаемых коррелированных ошибок. Рис. 3. Зависимость величины относительной СКО σ от шага по дальности при N = 64 и NL = 8 на расстояниях xн = 32 (кр. 1, 2) и xн = 128 (кр. 3, 4) при использовании традиционных КПГ и окна Хэмминга (кр. 1, 3) и предложенного алгоритма (кр. 2, 4) Если же вместо оптимальных значений K и L использовать традиционные значения (4) и окно Хэмминга соответственно, то с увеличением ∆x величина СКО σ растет значительно быстрее, что является прямым подтверждением эффективности применения оптимальных КПГ и ПС. К сожалению, точность расчетов (7) на практике недостижима, так как на m-м шаге по дальности входным воздействием являются отсчеты поля, рассчитанные на (m-1)-м шаге, а не результаты тестового расчета. Если учесть этот фактор в вышеописанном алгоритме подбора КПГ и ПС, то, как показывают предварительные выборочные расчеты, абсолютные значения СКО естественно увеличатся (рис. 4). Однако характер изменения - рост до некоторого значения μ с дальнейшим уменьшением - сохраняется. Максимальное значение СКО можно оценить, согласно , (9) Рис. 4. Зависимость величины относительной СКО σ от нормированного расстояния при N = 64, ∆xн = 1, числе отсчетов в ПС NL = 4 (кр. 1, 4) и NL = 8 (кр. 2, 5) при использовании традиционной (кр. 1, 2) и предложенной (кр. 4, 5) схем расчета; кр. 3 - то же, что и кр. 2, при N = 1024 где θ - величина, определяемая в (7). При этом точность расчетов значительно выше той, что достигается с применением традиционных КПГ (4) и окна Хэмминга. Более того, даже довольно значительное увеличение размеров расчетной области в традиционной схеме [14] не приводит к выигрышу в точности по отношению к предложенному алгоритму (см. кривые 2 и 3 на рис. 4). Результаты даже этих небольших по объему исследований обнадеживают и указывают на перспективность проведения более масштабных исследований в предположениях, приближенных к условиям практического использования метода. Заключение Приведенные результаты показывают, что имеется возможность существенного повышения точности расчета напряженности поля в рамках традиционного метода ДПФ и даже в отдельных случаях без всяких дополнительных затрат, при этом лучшие результаты получаются, если учитывать дистанционную зависимость матриц КПГ и ПС в (3). Выполнить в неоднородной среде это можно лишь в случае, если от решения уравнения в частных производных (2) перейти к более приемлемой идеологии решения эквивалентного интегрального уравнения [11,15]. Естественно, это потребует некоторого дополнительного увеличения объема вычислений. Описанная методика моделирования и полученные результаты показывают, что эти данные действительно определяют потенциальную, хотя и далеко не всегда достижимую, точность расчета напряженности поля в среде, которая в среднем является однородной, поскольку любые случайные неоднородности в среднем приводят к расширению спектра волн и, как следствие, к увеличению СКО расчета. Если это условие нарушается, СКО может как увеличиться, так и уменьшиться. Например, если среднее значение диэлектрической проницаемости при удалении от оси Ox уменьшается по параболическому закону (такая среда называется гармоническим осциллятором [28]), даже без проведения расчетов можно утверждать, что это уменьшит СКО. Отметим, что мы представили данные лишь для «минимаксного» варианта: при наличии взаимной корреляции источников, либо различий их интенсивностей, потенциальная точность будет увеличиваться (пока проверено на отдельных примерах). В работе приведены результаты только для двумерного пространства, получить подобные для трехмерного не представляет принципиальных трудностей, но при этом существенно возрастет и без того большой объем вычислений при проведении моделирования.

Скачать электронную версию публикации