Задачи на сфере: обратная геодезическая задача

Эта страница является черновиком статьи.

Обратная геодезическая задача — это нахождение начального направления и расстояния между двумя точками с известными координатами.

Введение

В качестве модели Земли принимается сфера с радиусом R, равным среднему радиусу земного эллипсоида. Введём следующие обозначения:

• φ — географическая широта,
• λ — географическая долгота,
• α — азимут дуги большого круга,
• σ — сферическое расстояние (длина дуги большого круга, выраженная в долях радиуса шара).

Линейное расстояние по дуге большого круга s связано со сферическим расстоянием σ формулой s = R σ.

Прямая и обратная геодезические задачи являются важными элементами более сложных геодезических задач.

Постановка задачи

x

Исходные данные

координаты пунктов Q₁ и Q₂ на сфере — φ₁, λ₁ и φ₂, λ₂.

Определяемые величины

расстояние между пунктами σ, начальный азимут направления α₁ в точке Q₁ на пункт Q₂.

На рисунке синим цветом выделены заданные элементы сферического треугольника, красным цветом неизвестные.

Алгоритм

Существует огромное множество подходов к решению поставленной задачи. Рассмотрим простой и надёжный векторный метод.

Последовательность решения:

1. преобразовать углы (λ₂ − λ₁) и φ₂ в декартовы координаты,
2. развернуть координатные оси на угол 90° − φ₁ вокруг оси Y,
3. преобразовать декартовы координаты в сферические,

Пример реализации алгоритма в виде функции языка Си:

/*
 * Решение обратной геодезической задачи
 *
 * Аргументы исходные:
 *     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;
}

В этом примере вычисление разности долгот заменено ещё одним вращением.

Следует заметить, что алгоритм решения обратной задачи является «перевёртышем» алгоритма прямой задачи.

Преобразование сферических координат в декартовы

${\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: ${\begin{array}{rcl}x'&=&-z\cos \phi _{1}+x\sin \phi _{1}y'&=&yz'&=&z\sin \phi _{1}+x\cos \phi _{1}\end{array}}$

Преобразование декартовых координат в сферические

Невозможно разобрать выражение (неизвестная функция «\begin{array}»): {\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°.

Альтернативные методы

Некоторые элементы альтернативных методов решения обратной задачи представлены в статье [[Вычисление_расстояния_и_начального_азимута_между_двумя_точками_на_сфере Вычисление расстояния и начального азимута между двумя точками на сфере]]

Ссылки

• Вычисление расстояния и начального азимута между двумя точками на сфере
• Степанов Н. Н. Сферическая тригонометрия
• Earth radius