Известия РАН. Серия физическая, 2021, T. 85, № 3, стр. 372-377

ВВЕДЕНИЕ

В настоящее время численное моделирование распространения коротких радиоволн (КВ) является одним из главных инструментов интерпретации и анализа ионограмм вертикального и наклонного зондирования ионосферы. На основе численного моделирования распространения КВ проводится валидация ионосферных моделей, исследуется динамика ионосферной плазмы, а также осуществляется прогноз оптимальных параметров радиосвязи. Для задач имитационного моделирования широко используется приближение геометрической оптики [1], основанное на решении уравнения эйконала. Наиболее распространенным и широко используемым подходом к решению уравнения эйконала является метод бихарактеристик [2]. Метод бихарактеристик позволяет быстро и с высокой точностью рассчитывать лучевые траектории в трехмерной неоднородной анизотропной плазме, поэтому используется в многочисленных исследованиях [3–5]. Для решения бихарактеристической системы лучевых уравнений необходимо задать начальные условия, которые в двухточечной (краевой) задаче заранее неизвестны. Наиболее распространенным подходом к решению краевой задачи с заданными положениями передатчика и приемника является численное решение системы лучевых уравнений в сочетании с методом последовательных приближений или пристрелки, где направление запуска систематически уточняется до тех пор, пока луч не попадет в указанную точку приема с заданной точностью. До недавнего времени пристрелка являлась единственным широко используемым подходом для решения краевой задачи о расчете лучевых траекторий КВ. Существуют вариации оптимизации такого алгоритма. Например, в работе [6] авторы рассматривают применение симплексной минимизации методом Нелдера–Мида [7] относительно азимута и угла возвышения. На основе метода пристрелки выполняется численное моделирование ионограмм наклонного зондирования [6, 8]. Однако в условиях нелинейной краевой задачи метод пристрелки является слабоустойчивым, особенно для лучей с большой расходимостью [9, 10], в результате чего решение граничной задачи требует значительных временных и вычислительных затрат [11, 12]. Поэтому поиск альтернативных методов решения граничной задачи является актуальной задачей.

Другим перспективным методом решения краевой задачи является вариационный подход, в основе которого лежит принцип Ферма о стационарности фазового пути луча. Главным достоинством такого подхода является точное выполнение граничных условий. Ранее были предложены реализации вариационного метода расчета лучевых траекторий радиоволн [9, 10, 13, 14]. Наиболее значимый вклад в развитие и применение вариационного метода для расчета радиотрасс оказала работа Коулмана [15]. Предложенный Коулманом комплексный вариационный подход позволяет использовать простую минимизацию для поиска “верхних” (лучи Педерсена) и трансионосферных лучей и численное решение вариационного уравнения для “нижних” лучей. Кроме того, реализация вариационного метода Коулмана учитывает анизотропию ионосферы. Однако представленные выше реализации вариационного подхода эффективны только в определении верхних и трансионосферных лучей, фазовый путь которых соответствует минимуму функционала Ферма. В то же время определение односкачковых нижних, многоскачковых, волноводных радиолучей оптимизационными методами на основе вариационного подхода на протяжении долгого времени представляло проблему. Другие трудности реализации вариационного метода были связаны с многолучевостью ионосферного распространения. Для определения множества лучевых траекторий в неоднородной среде требуется организовать выборку начальных приближений луча, что также вызывает трудности. В работе [16] представлены решения важнейших проблем вариационного метода: найден способ определения нижних лучей, решены проблемы многолучевости и подбора начальных приближений на основе “глобальной оптимизации”.

Ранее нами был рассмотрен вопрос, как изменение количества точек вдоль трассы, построенной вариационным методом, влияет на угол возвышения α. Значения полученных α использовались в качестве входных параметров для модели, основанной на методе бихарактеристик [17]. Было исследовано, как меняется решение для каждого α при изменении параметра нормализованного значения начального шага для длины вдоль трассы Δτ по модели, основанной на методе бихарактеристик. Было определено, что оптимальным является Δτ = 0.1, когда получаемые траектории методом бихарактеристик совпадают или достаточно близко проходят с траекторией, полученной вариационным методом.

Целью данной работы является сравнительный анализ времени и точности решений краевой задачи о расчете радиолучей с зафиксированными точками передачи и приема традиционным методом пристрелки и перспективным вариационным методом. Данная работа позволит оценить преимущества и применимость различных подходов для расчета КВ радиотрасс.

МЕТОД ПРИСТРЕЛКИ

В качестве модели, основанной на решении уравнения эйконала методом бихарактеристик, использовалась [4, 18]. Для нахождения лучей, приходящих в точку приема, использовалась реализация алгоритма пристрелки, подробно описанная в статье [19], которая состоит в итеративном изменении азимута и угла места таким образом, чтобы максимально приблизить к приемнику точку прихода сигнала на Землю.

ВАРИАЦИОННЫЙ МЕТОД

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

Здесь интегрирование производится вдоль некоторой кривой $\gamma $, задающей траекторию луча, которая соединяет точки A и B, $n\left( {\vec {r}} \right)$ – показатель преломления в точке $\vec {r} = \left( {x,y,z} \right)~$, лежащей на кривой $\gamma $, и $dl$ – элемент длины вдоль $\gamma $. Поскольку положение первой и последней точки зафиксировано, задача о нахождении траектории радиолуча сводится к поиску экстремумов (стационарных точек) функции $S\left( {\vec {r}} \right) = S\left( {{{{\vec {r}}}_{1}},{{{\vec {r}}}_{2}}, \ldots ,{{{\vec {r}}}_{P}}} \right)$, где P – число “подвижных” точек, задающих кусочно-линейное представление луча. Ранее было показано, что радиолучи могут соответствовать двум типам стационарных точек [20]: минимумам и седловым точкам первого порядка. Для поиска стационарных точек различного типа был предложен метод обобщенной силы [16].

Поиск минимума в рамках метода обобщенной силы, основан на модификации отрицательного градиента функционала $S\left( {\vec {r}} \right)$. Полученный фиктивный вектор силы $\vec {F} = \left( {{{{\vec {F}}}_{1}},~{{{\vec {F}}}_{2}}, \ldots ,~{{{\vec {F}}}_{P}}~} \right)$, действующий на каждую точку кусочно-линейной кривой $\vec {r}$ задает направление оптимизации лучевой траектории. Реализация процесса оптимизации представляет собой последовательные смещения траектории от некоторого начального положения ${{\vec {r}}^{{{\text{initial}}}}}$ к оптимальной конфигурации ${{\vec {r}}^{{{\text{final}}}}}$, соответствующей нулю вектора силы $\vec {F}$. В нашей работе сходимость к оптимуму осуществляется методом проецирования скорости [21]. В результате оптимальная конфигурация луча соответствует минимуму функционала $S\left( {\vec {r}} \right)$ и удовлетворяет принципу Ферма.

В рамках метода обобщенной силы для идентификации седловых точек первого порядка дополнительно используется численный расчет гессиана $H$ (матрицы вторых производных) целевой функции $S\left( {\vec {r}} \right)$. Оценка знакопостоянства матрицы Гессе позволяет определить направления и скорость возрастания/убывания функции $S\left( {\vec {r}} \right)$. В результате вектор силы $\vec {F}$ может быть скорректирован в направлении поиска седловых точек. Важно отметить, что реализация поиска седловых точек отличается от прямой минимизации только способом расчета силы $\vec {F}$. Именно это обстоятельство позволило реализовать единый подход к поиску, как минимумов, так и седловых точек первого порядка, основанный на различных модификациях многомерного вектора $\vec {F}$. Кроме того, переходные свойства седловых точек позволили на основе метода обобщенной силы создать алгоритм систематического поиска лучей без необходимости подбора начальных приближений. Данный алгоритм получил название глобальная оптимизация и подробно описан в [16].

ПОСТАНОВКА ЗАДАЧИ

Для сравнения двух численных методов была выбрана изотропная модель ионосферы, где профиль электронной концентрации N_e зависит только от высоты, $h$, над поверхностью Земли [3]:

где h₀ = 300 км – высота максимума N_e на высотах ионосферы, h_m = 150 км – полутолщина слоя, N_max = 1 ⋅ 10⁶ см^–3 – значение N_e в максимуме параболического слоя. Критическая плазменная частота такой среды приблизительно равна 9 МГц. Используемые модель среды и численные модели распространения радиоволн учитывают сферичность Земли. Для выполнения численного эксперимента была выбрана радиотрасса Ловозеро (68.0° с.ш., 35.02° в.д.) – Горьковская (60.27° с.ш., 29.38° в.д.). Максимально применимая частота (МПЧ) этой радиотрассы при распределении электронной концентрации, описанному выше, оценена на основе трассировки лучей методом бихарактеристик и составляет 13 МГц. Численный эксперимент проводился для радиоволны с частотой 12 МГц, что соответствует диапазону между критической и МПЧ рассматриваемой радиотрассы. Для выполнения численного расчета и нахождения всех решений мы применили метод бихарактеристик совместно с пристрелкой и вариационный метод.

РЕЗУЛЬТАТЫ

Процесс поиска решений методом пристрелки состоял из двух этапов. На первом этапе проводилась трассировка в заданном по дуге большого круга азимутальном направлении на точку приема. В данном численном эксперименте указанный способ трассировки обоснован изотропным приближением среды и отсутствием горизонтальных градиентов электронной концентрации. Трассировка выполнялась для разных углов места α, изменяемых с шагом Δα = 1°, и позволила определить два решения (необходимых для запуска второго этапа – алгоритма последовательных приближений), когда дальность между точкой излучения и точкой прихода больше и меньше расстояния между двумя заданными точками передачи и приема. На рис. 1 показана иллюстрация работы метода пристрелки как для нижнего (а), так и для верхнего (б) лучей. После работы второго этапа искомый луч попадает в некоторую окрестность, радиус которой задается как требуемая точность расчетов (0.5 км для наших расчетов).

Расчет траекторных характеристик вариационным методом осуществлялся в соответствии с алгоритмом глобальной оптимизации (см. рис. 2). Начальное приближение задавалось в виде прямой линии, соединяющей положения передатчика и приемника. На первом этапе глобальной оптимизации методом обобщенной силы было определено первое решение – нижний луч, соответствующий седловой точке первого порядка. Переход от начального приближения к первому решению осуществлялся кусочно-линейной кривой с дискретизацией D = 5. Для увеличения точности расчетов найденного нижнего луча, дискретизация траектории последовательно увеличивалась с 5 до 61 точек с шагом 10 точек. Максимальное значение дискретизации равно числу точек траектории, полученной пристрелкой и методом бихарактеристик (см. табл. 1). Результаты численных расчетов нижнего луча вариационным подходом приведены в табл. 2 . На втором этапе глобальной оптимизации осуществлялся переход от найденной седловой точки (нижний луч) к ближайшему минимуму. Оптимизация методом обобщенной силы определила локальный минимум, соответствующий верхнему лучу. Аналогично предыдущему этапу дискретизация последовательно увеличивалась с 5 до 113 точки с шагом 10 точек.

В результате численного эксперимента были рассчитаны основные параметры траекторных характеристик радиотрассы Ловозеро–Горьковская методом пристрелки и вариационным методом. Анализ результатов показывает полное согласие дистанционных и временных параметров решений, полученных двумя методами. Разница между значениями группового и геометрического путей не превышает 0.5 и 3 км для нижнего и верхнего лучей соответственно. Полное согласие наблюдается по параметрам времени задержки и максимальной высоты отражения лучей. Разница между найденными значениями параметра угла места для нижнего и верхнего лучей не превышает 0.9° и 0.2° соответственно. Указанные различия в значениях параметров могут быть объяснены особенностями исследуемых моделей, осуществляющих численные расчеты лучей, а также наличием у метода пристрелки области “попадания”, с заданным конечным радиусом. В результате этого, полученные пристрелкой лучи не всегда могут точно попасть в точку приема (см. табл. 1).

Одним из ключевых параметров численного эксперимента является время расчета лучей. Моделирование проводилось на ЭВМ с процессором Intel Core i7-7700 с тактовой частотой 3.6 ГГц и объемом оперативной памяти 8 ГБ. Полученные результаты показывают, что первый этап метода пристрелки наиболее времязатратный. Согласно табл. 1, время трассировки до нахождения “вилки” для нижнего луча при шаге Δα = 1° занимает порядка 75–90% общего времени поиска решения, а для верхнего – порядка 50–90%. Следовательно, время, требуемое для нахождения верхнего луча на втором этапе, может быть сопоставимо со временем трассировки. Увеличение шага Δα приведет к снижению времени трассировки. Однако шаг не должен быть слишком большим, иначе некоторые решения могут быть потеряны. Отметим, что время расчета зависит от безразмерного параметра Δτ, влияющего на дискретизацию луча через пространственные градиенты показателя преломления. Как показано в табл. 1 изменение параметра Δτ в диапазоне 0.5 до 0.01 приводит к значительному росту временных затрат: от 16.8 с до 15.2 мин для нижнего луча и от 35.6 с до 18.9 мин для верхнего луча. Одновременно наблюдается значительное увеличение дискретизации и точности расчета луча. В данном случае оптимальное значение параметра Δτ составляет 0.1.

Рассмотрим временные затраты расчетов, выполненных вариационным методом (см. табл. 2 и 3). Диапазон времени расчета верхнего луча с дискретизацией от 5 до 113 точек составляет 0.1–3.1 с. Суммарное время определения верхнего луча в результате последовательного увеличения дискретизации составило 21.5 с. Однако, согласно табл. 2 время расчета нижнего луча значительно увеличивается. Для дискретизации луча от 5 до 61 точки время расчета составляет 3.1–29.8 с. Суммарное время определения нижнего луча в результате последовательного увеличения дискретизации составило 2.9 мин, что сопоставимо с результатами по методу пристрелки (включая время предварительной трассировки). В данном случае увеличение времени расчета связано с необходимостью численного расчета матрицы вторых производных, что неизбежно ведет к росту вычислительных затрат.

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

ЗАКЛЮЧЕНИЕ

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

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

Авторы выражают благодарность Бессарабу П.Ф. за неоценимый вклад в развитие вариационного метода и многолетнее сотрудничество. Авторы признательны Бессарабу Ф.С. и Карпову И.В. за внимание и полезные дискуссии. Численное моделирование выполнено при финансовой поддержке гранта Президента Российской Федерации для молодых ученых № МК-2584.2019.5 (Котова Д.С., Носиков И.А.). Разработка и развитие метода обобщенной силы и глобальной оптимизации выполнено при финансовой поддержке Российского научного фонда (проект № 17-77-20009 – Клименко М.В.).