Задачи на сфере: обратная геодезическая задача
Обратная геодезическая задача — это нахождение начального направления и расстояния между двумя точками с известными координатами.
Введение
В качестве модели Земли принимается сфера с радиусом R, равным среднему радиусу земного эллипсоида. Введём следующие обозначения:
- φ — географическая широта,
- λ — географическая долгота,
- α — азимут дуги большого круга,
- σ — сферическое расстояние (длина дуги большого круга, выраженная в долях радиуса шара).
Линейное расстояние по дуге большого круга связано со сферическим расстоянием формулой Невозможно разобрать выражение (SVG (MathML можно включить с помощью плагина для браузера): Недопустимый ответ («Math extension cannot connect to Restbase.») от сервера «https://wikimedia.org/api/rest_v1/»:): {\displaystyle s = R σ} .
Прямая и обратная геодезические задачи являются важными этапами в решении более сложных геодезических задач.
Постановка задачи
- Исходные данные
- координаты пунктов Q₁ и Q₂ на сфере — φ₁, λ₁ и φ₂, λ₂.
- Определяемые величины
- расстояние между пунктами σ, начальный азимут направления α₁ в точке Q₁ на пункт Q₂.
На рисунке синим цветом выделены заданные элементы сферического треугольника, красным цветом неизвестные.
Алгоритм
Существует огромное множество подходов к решению поставленной задачи. Рассмотрим простой и надёжный векторный метод.
Последовательность решения:
- преобразовать углы (λ₂ − λ₁) и φ₂ в декартовы координаты,
- развернуть координатные оси на угол 90° − φ₁ вокруг оси Y,
- преобразовать декартовы координаты в сферические,
Пример реализации алгоритма в виде функции языка Си:
/*
* Решение обратной геодезической задачи
*
* Аргументы исходные:
* pt1 - {широта, долгота} точки Q₁
* pt2 - {широта, долгота} точки Q₂
*
* Аргументы определяемые:
* azi - азимут начального направления
* dist - расстояние (сферическое)
*/
void SphereInverse(double pt1[], double pt2[], double *azi, double *dist)
{
double x[3], pt[2];
SpherToCart(pt2, x); // сферические -> декартовы
Rotate(x, pt1[1], 2); // (первое вращение)
Rotate(x, M_PI_2 - pt1[0], 1); // (второе) вращение
CartToSpher(x, pt); // декартовы -> сферические
*azi = M_PI - pt[1];
*dist = M_PI_2 - pt[0];
return;
}
В этом примере вычисление разности долгот заменено ещё одним вращением.
Следует заметить, что алгоритм решения обратной задачи является «перевёртышем» алгоритма прямой задачи.
Преобразование сферических координат в декартовы
Невозможно разобрать выражение (Ошибка преобразования. Сервер («https://wikimedia.org/api/rest_») сообщил: «Cannot get mml. TeX parse error: Double subscripts: use braces to clarify»): {\displaystyle {\begin{array}{rcl}x&=&\cos \phi _{2}\cos _{(}\lambda _{2}-\lambda _{1})\\y&=&\cos \phi _{2}\sin _{\lambda }_{2}\\z&=&\sin \phi _{2}\end{array}}}
Реализация на Си:
/*
* Преобразование сферических координат в вектор
*
* Аргументы исходные:
* y - {широта, долгота}
*
* Аргумены определяемые:
* x - вектор {x, y, z}
*/
void SpherToCart(double y[], double x[])
{
double p;
p = cos(y[0]);
x[2] = sin(y[0]);
x[1] = p * sin(y[1]);
x[0] = p * cos(y[1]);
return;
}
Вращение вокруг оси
Собственно работа векторного метода выполняется вращением вокруг координатной оси Y:
Преобразование декартовых координат в сферические
Невозможно разобрать выражение (SVG (MathML можно включить с помощью плагина для браузера): Недопустимый ответ («Math extension cannot connect to Restbase.») от сервера «https://wikimedia.org/api/rest_v1/»:): {\displaystyle \begin{array}{rcl} α₁ & = & 180^{\circ} − \operatorname{arctg} \frac{y}{x} \\ σ & = & 90^{\circ} − \operatorname{arctg} \frac{z}{\sqrt{x^2 + y^2}} \end{array}}
/*
* Преобразование вектора в сферические координаты
*
* Аргументы исходные:
* x - {x, y, z}
*
* Аргументы определяемые:
* y - {широта, долгота}
*
* Возвращает:
* длину вектора
*/
double CartToSpher(double x[], double y[])
{
double p, z;
z = x[2];
p = sqrt(x[0] * x[0] + x[1] * x[1]);
y[1] = atan2(x[1], x[0]);
y[0] = atan2(z, p);
return sqrt(p * p + z * z);
}
Пример программной реализации
Исходники вышеприведённых функций можно найти в архиве в файле sph.c. Кроме того, в файл sph.h включены следующие определения:
#define A_E 6371.0 // радиус Земли в километрах
#define Degrees(x) (x * 57.29577951308232) // радианы -> градусы
#define Radians(x) (x / 57.29577951308232) // градусы -> радианы
Теперь напишем программу, которая обращается к функции SphereInverse для решения обратной задачи:
#include <stdio.h>
#include <stdlib.h>
#include "sph.h"
int main(int argc, char *argv[])
{
char buf[1024];
double pt1[2], pt2[2];
double lat1, lon1, lat2, lon2, azi1, azi2, dist;
while (fgets(buf, 1024, stdin) != NULL) {
sscanf(buf, "%lf %lf %lf %lf", &lat1, &lon1, &lat2, &lon2);
pt1[0] = Radians(lat1);
pt1[1] = Radians(lon1);
pt2[0] = Radians(lat2);
pt2[1] = Radians(lon2);
SphereInverse(pt2, pt1, &azi2, &dist);
SphereInverse(pt1, pt2, &azi1, &dist);
printf("%f\t%f\t%.4f\n", Degrees(azi1), Degrees(azi2), dist * A_E);
}
return 0;
}
В архиве этот код находится в файле inv.c. Создадим исполняемый модуль inv компилятором gcc:
$ gcc -o inv inv.c sph.c -lm
Для MS Windows готовую программу inv.exe можно найти в архиве.
Программа читает данные со стандартного ввода консоли и отправляет результаты на стандартный вывод. Для чтения и записи файлов используются символы перенаправления потока «>» и «<» соответственно. Из каждой строки ввода программа считывает координаты двух точек φ₁, λ₁, φ₂, λ₂, которые должны быть в градусах, решает обратную задачу и записывает в строку вывода α₁, α₂, s (азимуты прямого и обратного направлений в градусах и расстояние между пунктами в единицах, определённых константой A_E).
Создадим файл inv.dat, содержащий одну строку данных:
30 0 52 54
После запуска программы
$ inv < inv.dat
получим α₁, α₂, s:
44.804060 262.415109 5001.1309
Решение обратной задачи средствами PROJ.4
В пакет PROJ.4 входит программа geod, предназначенная для решения прямых и обратных геодезических задач на сфере. Так выглядит команда обработки файла inv.dat:
$ geod +a=6371000 -I -f "%f" -F "%.4f" +units=km inv.dat
Параметр +a определяет радиус сферы, -I — решение обратных задач, -f — формат вывода угловых величин, -F — формат вывода длин линий, +units — единица измерения расстояний. В результате получим идентичный результат:
44.804060 -97.584891 5001.1309
Отличие значений α₂ на 360° объясняется тем, что inv выводит азимуты в диапазоне от 0° до 360°, а geod от −180° до +180°.
Альтернативные методы
Некоторые элементы альтернативных методов решения обратной задачи представлены в статье [[Вычисление_расстояния_и_начального_азимута_между_двумя_точками_на_сфере Вычисление расстояния и начального азимута между двумя точками на сфере]]