Содержание к диссертации
Введение
1 Обзор систем прогнозирования характеристик СВЧ радиоволн 18
1.1 Современные системы прогнозирования 18
1.2 Обзор существующих численных методов 20
1.3 Стандартное и широкоугольное параболические уравнения .. 21
1.4 Особенности методов решения параболического уравнения в частотной области 25
1.5 Особенности методов решения параболического уравнения в пространственной области 27
1.6 Погрешности при использовании метода расщепления 28
1.7 Задачи исследований 30
2 Граничные условия при решении двумерного параболического уравнения над неровной земной поверхностью 31
2.1 Метод пошагового преобразования системы координат 31
2.2 Метод конформного отображения 35
2.3 Результаты сравнительного численного расчёта поля методами пошагового преобразования координат и конформного отображения 42
2.4 Сравнение влияния неоднородностей подстилающей поверхности и тропосферы на характеристики поля 45
2.5 Частичная оптимизация формы и параметров дискретного поглощающего слоя 51
2.6 Дискретное прозрачное граничное условие в пространственной области (схема Кранка-Николсон) 59
2.7 Дискретное прозрачное граничное условие в спектральной области (метод преобразования Фурье) 73
3 Разработка макета системы оперативного прогнозирования характеристик радиосигналов 80
3.1 Схема Кранка-Николсон и метод преобразования Фурье 81
3.2 Средства сбора данных о состоянии нижнего слоя тропосферы вблизи трассы распространения и о земной поверхности 97
3.3 Программный комплекс для расчёта основных характеристик электромагнитного поля 99
4 Экспериментальная проверка результатов моделирования
4.1 Одновременные измерения метеопараметров тропосферы (2 120 м) и множителя ослабления ЮЗ
4.2 Измерения матричных импульсных характеристик канала распространения радиоволн 112
Заключение 129
- Стандартное и широкоугольное параболические уравнения
- Результаты сравнительного численного расчёта поля методами пошагового преобразования координат и конформного отображения
- Дискретное прозрачное граничное условие в пространственной области (схема Кранка-Николсон)
- Программный комплекс для расчёта основных характеристик электромагнитного поля
Стандартное и широкоугольное параболические уравнения
Как уже было сказано во введении, основным современным способом расчёта поля СВЧ радиоволн, учитывающим неоднородности среды распространения, является метод ПУ. Это уравнение является линейным неоднородным дифференциальным уравнением в частных производных второго порядка, относящееся к классу параболических уравнений [16, и. 10.3-3]. Данный класс предполагает знание только одного начального условия, которое в рассматриваемом случае определяет источник излучения (антенна). Существующие аналитические методы решения ПУ [16, и. 10.3 - 10.5] не годятся, т.к. практическую ценность система прогнозирования приобретает, в частности, тогда, когда она не ограничена известными классами подстилающих поверхностей и неоднородностей среды распространения. Поэтому главная роль переходит к численным методам, пик развития которых приходится на 60-80-е г.г. XX столетия, т.е. после появления первых ЭВМ.
Сутью численных методов является переход от непрерывных функций к сеточным [17]. Это же относится и к их производным требуемого порядка. При таком раскладе неизбежно возникает методическая погрешность расчётов, какой бы разрядностью не обладал вычислительный процессор. Эта погрешность тем меньше, чем меньше наибольший объём ячейки расчётной сетки при условии, что разрядность процессора достаточно велика. При этом возрастает объём вычислений. В каждом случае необходимо индивидуально оценивать число операций и ошибку рассматриваемого численного метода при определённой разрядности процессора. После введения сетки (сеточных функций, определённых на ней) заданное дифференциальное уравнение можно решать двумя споеобами.
Первым способом является непосредственная замена непрерывной функции и её производных выбранными конечно-разностными аналогами, составление системы линейных алгебраических уравнений (СЛАУ) и нахождение алгоритма её решения. Отличия в алгоритмах решения СЛАУ приведут к различным численным методам: методу исключения Гаусса, методу прогонки как частному случаю метода Гаусса, итерационным методам.
Второй способ отличается от первого тем, что вместо составления СЛАУ используется некоторое функциональное преобразование, переводящее соответствующее разностное уравнение в так называемую спектральную область (преобразование Фурье, -преобразование). После этого преобразования разностное уравнение переходит в алгебраическое, которое решается соответствующими методами. Затем найденный спектр решения преобразуется в исходную пространственную область. Или, что чаще всего, заранее находят выражение для коэффициента передачи, которое используют в качестве множителя при пересчёте спектра исследуемого поля на следующий щаг по определённой координате. К спектральному преобразованию К предъявляется требование существования обратного преобразования КГ1, т.е. Гг[К [м]] = и.
Таким образом, делаем вывод, что выбор численных методов решения ПУ невелик. Это или метод прогонки как имеющий наименьшее количество операций (используется преимущество ленточной структуры матрицы СЛАУ), или метод преобразования Фурье, достоинством которого является возможность использования оптимизированных для распараллеливания БПФ графических процессоров (GPU). Более подробно эти методы рассмотрены в разд. 1.4-1.5 и 3.2 данной работы.
В данном разделе речь пойдёт о выводе ПУ из уравнения Гельмгольца и о соответствующих ограничениях применимости ПУ.
Исторически первые работы, посвящённые методу ПУ, вышли под руководством М.А. Леонтовича и В.А. Фока в 40-е годы прошлого столетия. В этих работах показано, что задачу о РРВ над сферической полупроводящей поверхностью можно решить сравнительно простым методом ПУ. До решения поставленной задачи методом ПУ, Фок нашёл решение, следующее из медленно сходящегося ряда для вектора Герца, который получается из исходных уравнений Максвелла (рещение Ми задачи дифракции на шаре). Фок использовал представление данного ряда контурным интегралом, в котором заменил функции Ганкеля и Бесселя их асимптотиками через им введённые функции Эйри, используя для этого «большой параметр» задачи (отношение радиуса кривизны земной поверхности к длине волны поля). В совместных с Леонтовичем работах Фок показал, что найденное им решение можно получить методом ПУ. Это решение можно использовать для любых удалений от передатчика: для малых удалений выражение переходит в интерференционные формулы, а для больших удалений - в одночленную формулу Б.А. Введенского.
Таким образом, метод ПУ является универсальным методом решения задачи распространения радиоволн для случаев, когда радиус кривизны подстилающей поверхности и размеры неоднородностей тропосферы много больше длины волны. При выполнении этих условий можно пренебречь продольной диффузией поля по сравнению с поперечной. При этом из физических соображений следует, что угол рассеяния будет мал (например, дифракция света на больших по сравнению с длиной волны препятствиях [9]). Приведём ниже известный [9, 14] вывод ПУ из уравнения Гельмгольца. При этом сначала найдём вид стандартного (малоуглового) ПУ, а затем покажем путь перехода к вариантам широкоугольного.
Результаты сравнительного численного расчёта поля методами пошагового преобразования координат и конформного отображения
Рельеф земной поверхности, показанный на рис. 2.9, является шероховатой, в целом плоской, поверхностью, характер отражений от которой зависит от угла скольжения. При малых углах преобладает зеркальное отражение, что мы и наблюдаем (кривая С на рис. 2.9) до высот 20 м в виде двулучевой интерференционной картины, а при больших углах индикатриса рассеяния расширяется, что приводит к искажению интерференционной картины поля на остальных высотах.
Таким образом, даже в том случае, когда земная поверхность является почти плоской и слабошероховатой, её влияние на характеристики поля сопоставимо с влиянием неоднородностей индекса преломления тропосферы, СКО которых на два порядка превышает реально наблюдаемые значения. То есть на коротких трассах прямой видимости и при малых углах рассеяния можно не учитывать влияние случайных неоднородностей индекса преломления тропосферы (речь идёт о диапазоне частот 100 МГц ... 30 ГГц).
Что касается больших углов РРВ (более 10 градусов), то теоретических результатов оценки влияния неоднородностей тропосферы нет по причине нерешённости уравнения (2.16) относительно двух зависимых переменных Ап и ф„ (квадратурных составляющих Acosq „ и Ansinq n). В этом случае необходимо решить эквивалентное матричное уравнение
Численное решение ПУ ееточным методом предполагает ограниченность области расчёта, что требует задания ГУ. Нижняя граница (подстилающая поверхность) определяется ГУ Леонтовича [35, 39], которое является локальным, т.е. граничное поле на текущем шаге по дальности не зависит от граничных значений поля на предыдущих щагах. Верхнюю границу области расчёта как правило можно считать открытой, т.к. при решении задачи РРВ в естественных средах верхнюю границу обычно выбирают так, чтобы область расчёта включала все наиболее значимые изменения показателя преломления. Открытость границы означает её прозрачность, т.е. распространяющаяся волна не должна от неё отражаться. Это требование приводит к нелокальности верхнего ГУ [28-30], что вызывает трудности нахождения и численной реализации оператора нелокального ГУ. Если данный оператор удаётся найти, то возникает проблема сокращения вычислительных затрат, т.к. его прямая реализация для двумерной задачи требует 0(N2) количества операций (дискретная свёртка).
Проблема прозрачных ГУ для схемы К-Н и двумерного ПУ решена в разд. 2.6. Для метода преобразования Фурье (т.е. в спектральной области) проблема нелокального ГУ до конца (с оптимизацией вычиелительных затрат) не решена. Решение проблемы нелокального ГУ, найденное в разд. 2.6, соверніенно не применимо к трёхмерному ПУ. Поэтому альтернативным вариантом для метода преобразовании Фурье является техника поглощающих слоёв, которая по сути своей является эмпирической.
Прямоугольная область расчёта с поглощающим слоем на верхней границе Параметры слоёв были подобраны так, чтобы ослабление на краях слоя не было равным нулю, т.е. чтобы амплитуда поглощения постепенно уменьшалась от единицы до заданной величины (слабый слой). Это требование определилось исходя из предварительных численных расчётов поля и сравнения с эталоном. Под эталоном в данном случае понимается аналитическое решение двумерного ПУ в свободном пространстве (приложение В).
Для численных расчётов был выбран дискретный точечный источник поля (с ДН прямоугольной формы; см. введение или Приложение В), расположенный в центре рабочей области (рис. 2.10). Ширина ДН по соответствующей координате определяется как 2аіС8Іп(Я./2Дг),рад. (2.35) где Az - шаг расчётной сетки по высоте (ширине). Ограничение шага Az снизу расстоянием полуволны связано с требованием для широкоугольного ПУ соблюдения углов РРВ от -90 до +90 градусов. Эти ограничения снимаются только при переходе к уравнению Гельмгольца (на другой качественный уровень), что в данной работе не имеет смысла.
Дальность (длина трассы РРВ) х нормировалась как xX/Az2, (2.36) где X - длина волны источника. Шаг Az и высота рабочей области фиксировались, поэтому число отсчётов по высоте также было неизменным. Для вычисления эталонного поля использовались результаты из приложения В. За меру ошибки было принято среднеквадратическое отклонение от эталона. Шаг по высоте Az = 1 м, шаг по дальности Дх = 5 м. Длина волны X = 0,1 м. Число отсчётов по высоте N = 256. Толщина слоя выбиралась из ряда: 16, 32, 64.
В табл. 2.2 для примера приведены результаты оптимизации слоя линейной формы (2.34), толщиной 16. Оптимизация заключается в подборе амплитуды слоя А при фиксировании числа отсчётов по дальности и толщины слоя, так, что обеспечивается минимум среднеквадратической ошибки. Параметры у, 5 и р пояснены ниже.
Дискретное прозрачное граничное условие в пространственной области (схема Кранка-Николсон)
Таким образом, метод прогонки позволяет обратить трёхдиагональную матрицу, используя 0(М) операций (примерно 5М) и найти вектор поля на текущем шаге по дальности i + 1 по вектору из предыдущего шага / .
Если прогонку начинать сверху (левая прогонка, [17]), то в качестве начального значения для прогоночного коэффициента ам можно взять предельное значение коэффициента а,- при итерациях (3.6), т.к. можно предположить, что прогонка вниз идёт из бесконечности и значение этого коэффициента уже установилось. Решая квадратное уравнение (3.6) относительно установившегося значения а, найдём единственный корень, по модулю меньший единицы ос=1 + (l-ylu )/2g. (3.8) Тем самым мы сэкономим часть времени, взяв сц.= const. Требование ос 1 выполняется для любого g (и даже комплексного) - это необходимо и достаточно для устойчивости метода прогонки [17]. Коэффициент передачи по угловому спектру для схемы К-Н l + 2g [l-cos(2npAz)] sinp К(Р) l+2g[l-cos(2np Z)], Р X (3-9) где (З - угол распространения плоской волны exp[/(хcos(3 + 7sinp)] по отношению к оси дальности х (рис. 3.1); р - пространственная частота. Коэффициент (3.9) есть отношение угловых спектров вектора поля на текущем шаге по дальности г + 1 и вектора на предыдушем шаге / . Коэффициент передачи (3.9) получается после 2-преобразования (3.4) по переменной у и нахождения отношения lfM (z) / ІҐ г). Далее переменная z-преобразoвaния заменяется экспонентой Фурье ехр(2тфДг).
Модуль коэффициента передачи (3.9) равен единице, а фаза зависит от пространственной частоты следуюшим образом Ф(р) = -2arctg 4Im(g)sm (yAz) 11+4Re( )sin (npAz) Формулы (3.9) и (3.10) пригодятся при сравнении точности схемы К-Н и метода преобразования Фурье (разд. 3.3), коэффициент передачи которого можно считать эталоном для свободного пространства. Показано [26], что точность схемы К-Н можно повысить за счёт введения реальной части в коэффициент g из (3.4). Это обсуждается в разд. 3.3.
Неоднородности тропосферы (коэффициент преломления п) учитываются с помощью методики расщепления [разд. 1.6], которая сводится к простому умножению вектора (матрицы) поля на комплексную экспоненту. Метод преобразования Фурье Методику численного решения ПУ с помощью преобразования Фурье опишем для двумерного случая, после чего обобщим её на трёхмерный.
Коэффициент передачи (3.9) для схемы К-Н был получен как следствие формул для конечно-разностной схемы (по сути это коэффициент передачи некоторого цифрового фильтра). Отсюда и ограничение по точности схемы К-Н, т.к. желательно получить коэффициент передачи для непрерывного ПУ. Метод преобразования Фурье позволяет преодолеть это ограничение, хотя если потом при решении использовать ДПФ, то - преодолеть только частично. Однако вклад в ошибку дискретности преобразования Фурье незначителен по сравнению с ошибкой дискретности схемы К-Н.
Коэффициент передачи для метода преобразования Фурье является точным и находится с помощью решения непрерывного однородного ПУ в свободном пространстве посредством интеграла Фурье по поперечной координате (высоте)
Формула (3.21) определяет прямое смешанное преобразование Фурье. Для расчёта поля и(х+Ах, z) угловой спектр V(p) умножаетея на коэффициент передачи (3.18), после чего от полученного произведения находится обратное смешанное преобразование Фурье. Выражение (3.21) должно быть дополнено вычислением дополнительного коэффипиента, без которого обратное преобразование в общем случае найти невозможно (см. ниже).
Для получения формулы обратного смешанного преобразования Фурье найдем из уравнения (3.20) поле u(z)
Таким образом, из (3.25) следует, что для нахождения обратного смешанного преобразования Фурье необходимо знать дополнительную константу С, поэтому в формулу прямого преобразования (3.21) необходимо добавить её вычисление. Слагаемое Cexp(-az) в (3.25) при Re(a) 0 (случай вертикальной поляризации при РРВ) определяет поверхностную быстро затухающую волну (аналогия с диэлектрическим волноводом).
Программный комплекс для расчёта основных характеристик электромагнитного поля
Для расчёта поля по методу ПУ была выбрана прямоугольная трёхмерная область, нижняя грань которой совпадает с минимальной высотой рельефа (64 м над уровнем моря), а две боковых удалены на 200 м от передатчика вправо и влево. Высота расчётной области - 320 м, дальность - 7000 м, ширина - 400 м. По числу отсчётов это 1344, 2000 и 1152 соответственно. Высота антенны передатчика над подстилающей поверхностью -2 м, приёмных антенн - 4,4 м. В качестве модели ДН передающей антенны была взята гауссовская ДН, шириной 3. Коэффициент отражения от подстилающей поверхности - минус единица. Рефракция -нормальная (высотный градиент - минус 40 У-ед./км).
Расчёт был сделан для диапазона частот 9640 ... 9960 МГц с шагом 10 МГц. В результате чего был получен комплексный спектр (рис. 4.19), состоящий из 33 отсчётов.
Рассчитанная по методу ПУ амплитудно-частотная характеристика семикилометрового канала РРВ в диапазоне 9640 ... 9960 МГц Путём вычисления обратного преобразования Фурье была оценена ИХ канала с разрешением 3,1нс (рис. 4.20). Частотная характеристика канала (рис. 4.19) перед вычислением ИХ центрировалась по частоте (восстановление аналитического сигнала по его спектру). Наблюдаемый (рис. 4.19) теоретический завал (3 - 4 дБ) на частоте 9960 МГц по сравнению с частотой 9640 МГц объясняется интерференцией, а точнее изменением основного периода интерференции, который рассчитывается для плоской горизонтальной трассы по известным формулам 123 7Д = Xid/ilh,) = 0,0311203-7000 / (2-2) 54,46 м, Т2 = X2d/(2hs) = 0,0301205-7000 / (2-2) 52,71 м, на величину Тх-Тг 1,75 м.
При расчёте ИХ учитывался лишь достаточно гладкий рельеф почвы. Более мелкие неоднородности растительного происхождения и даже линия электропередачи, пересекающая трассу на расстоянии 15 м от передатчика, учтены не были в связи с отсутствием корректной методики.
Диаграмма множителя ослабления (для частоты 9640 МГц) в плоскости «дальность-высота», проходящей через центр расчётной области, показана на рис. 4.21. передатчик - слева на высоте 16 м, приёмник - справа на высоте 68,4 м; рельеф - светло-коричневый; дальность - от О до 7000 м; высота -от О до 80 м
Из рис. 4.21 следует, что поле в области приёмника формируется, в основном, за счёт прямых лучей и лучей, отражённых от области, удалённой от передатчика на 500 ... 1000 м. На раестояниях от передатчика более 1000 м рельеф практически не влияет. Чтобы определить область влияния по поперечной координате, следует проанализировать диаграмму множителя оелабления в плоекости «дальность--координата» для выеоты, например, 15,4 м (рис. 4.22).
Анализ рис. 4.22 показывает, что эта облаеть ограничена величиной ±20 м по 7-координате. Также можно отметить поперечную асимметрию поля, связанную с соответствующей асимметрией рельефа. То, что метод ПУ учитывает боковые отражения, можно показать на примере расчёта поля с рельефом в виде продольного клина, раеположенного по центру трассы, при условии, что передатчик расположен, например, на 20 м правее этого центра (рис. 4.23 - 4.24). - методом привязного аэрологического зондирования проведены измерения высотных профилей индекса преломления приземного слоя тропосферы (высоты от 2 до 120 м). Измеренные профили близки к линейным со значениями градиентов в диапазоне от -0,07 до -0,03 N-ед./м и максимальными отклонениями от линейной аппроксимации ±4 N-ед. Проведённые измерения показали, что привязное аэрологическое зондирование возможно только при скорости ветра не более 1 м/с. - из анализа экспериментальных импульсных характеристик канала РРВ и характеристик, полученных путём моделирования (метод ПУ), следует, что и та, и другая характеристики указывают на наличие сигналов, рассеянных объектами вблизи трассы, но расчётная характеристика слабо учитывает отражения от более далёких объектов. Вероятно, это объясняется тем, что при расчёте учитывался лишъ достаточно гладкий рельеф почвы без учёта более мелких неоднородностей растительного происхождения (лес, кустарники).
Временным препятствием при внедрении подобных систем прогнозирования являются высокие вычислительные затраты, что снижает оперативность прогнозирования, но в ближайшем будущем эта трудность будет быстро преодолеваться. Например, обзоры развития профессиональных вычис-лительных видеокарт по технологии CUDA свидетельствуют о повышении скорости вычислений в 10 - 20 раз.
Например, текущий расчёт (в частности, рис. 4.19) проводился с шагом по высоте и Y-координате около 25 - 35 см при длине волны 3 см. Это возможно благодаря тому, что по методу НУ вычисляется комплексная огибающая высокочастотного поля и расчёт ведётся в небольшом угловом секторе (не более 15 градусов). При этом размер матрицы поля составил 1152x1344, а число шагов по дальности - 2000. Время, затраченное программой, составило 10 минут на двух ядрах Core2Duo 2,4GHz с общей загрузкой CPU на уровне 80%. Объём используемой памяти - 24 МБ (одинарная точность с плавающей точкой). Для сравнения был проведён подобный расчёт на восьми ядрах двухпроцессорного сервера на CPU Xeon 2,2GHz с общей загрузкой CPU 55% - время расчёта составило 4 минуты, т.е. всего в 2,5 раза меньше (неполная загрузка ядер связана с реализацией многопоточности в процедурах БПФ FFTW 3.2.2 [43]).
Если же теперь считать поле с шагом полуволны, то размер матрицы придётся увеличить до 20000x25000, т.е. в 320 раз. Время расчёта увеличится, по крайней мере, во столько же раз (53 часа для Core2Duo).
Для учёта значимых особенностей рельефа необходимо комбинирование методов расчёта. Метод НУ пригоден для учёта плавных неровностей рельефа и плавных изменений показателя преломления тропосферы (как правило до 1 - 2 км по высоте). При этом дальность должна быть такой, чтобы угол рассеяния был малым (не более 15 градусов для малоуглового ПУ и 45 градусов для некоторых широкоугольных вариантов). По сути, методика комбинирования лучевых методов и метода ПУ уже разработана и реализована [12], но дополнения расчётов методами ГТД или ГО для явно выделенных дискретных отражателей (строения, опоры и Т.П.), по всей видимости, никем сделано не было.
