Формула для решения обратной геодезической задачи

Основные понятия геодезии

Геодезическая сеть – система взаимосвязанных геодезических пунктов, которые расположены на поверхности Земли и служат для определения географических координат точек, а также для проведения геодезических измерений и построения картографических материалов.

Прямая геодезическая задача – это задача определения геодезического положения одной или нескольких точек на поверхности Земли по известным параметрам (например, координаты, углы) других точек.

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

Координаты – числовые значения, позволяющие определить положение точки на поверхности Земли. Обычно задаются в виде географической широты и долготы, а также высоты над уровнем моря.

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

Картографические материалы – графические изображения, отображающие географические объекты на поверхности Земли. Могут быть представлены в виде карт, планов, схем и других географических материалов.

Геоинформационные системы – комплекс программных и аппаратных средств, предназначенных для сбора, хранения, обработки, анализа и визуализации геодезической информации.

Топография – раздел геодезии, занимающийся определением и описанием природных и искусственных объектов на поверхности Земли, их размерами, формой и положением.

Рельеф – характерный облик поверхности Земли, представляющий собой перепады высот и наклонов между точками.

Примеры

Решение прямой геодезической задачи

Для решения прямой геодезической задачи, неоходимо создать объект класса .

var directEllipsoid = new DirectProblemService(new Ellipsoid());
var directSpheroid = new DirectProblemService(new Spheroid());

В качестве параметра в конструктор следует передать объект реализующий интерфейс , в котором задаются полярный и экваториальный радиус, а так же коэффициент полярного сжатия.

Для решения прямой задачи вызвать метод , в который передать в качестве параметров, начальную точку, азимут = направление и расстояние:

var point1 = new Point(15, 25, 53, CardinalLongitude.W, 28, 7, 38, CardinalLatitude.N);
var azimuth = 21;
var distance = 2000;
var directAnswer = directEllipsoid.DirectProblem(point1, azimuth, distance);

Ответ содержит вторую точку ортодромии и обратный азимут .

Решение обратной геодезической задачи

Для решения обратной геодезической задачи, неоходимо создать объект класса .

var inverseEllipsoid = new InverseProblemService(new Ellipsoid());
var inverseSpheroid = new InverseProblemService(new Spheroid());

Для решения обратной задачи вызвать метод , в который передать в качестве параметров, две точки:

var point1 = new Point(15, 25, 53, CardinalLongitude.W, 28, 7, 38, CardinalLatitude.N);
var point2 = new Point(59, 36, 30, CardinalLongitude.W, 13, 5, 46, CardinalLatitude.N);
var inverseAnswer = inverseEllipsoid.OrthodromicDistance(point1, point2);

Ответ содержит прямой и обратный азимуты , , а также расстояние между точками .

Вычисление точки пересечения ортодромий

Для вычисления точки пересечения ортодромий, неоходимо создать объект класса .

var intersectEllipsoid = new IntersectService(new Ellipsoid());
var intersectSpheroid = new IntersectService(new Spheroid());

Для рассчёта вызвать метод , в который передать в качестве параметров, по две точки для каждой из двух ортодромий:

var point1 = new Point(22, 36, 30, CardinalLongitude.E, 13, 5, 46, CardinalLatitude.N);
var point2 = new Point(27, 25, 53, CardinalLongitude.E, 15, 7, 38, CardinalLatitude.N);
var point3 = new Point(20, 36, 30, CardinalLongitude.E, 17, 5, 46, CardinalLatitude.N);
var point4 = new Point(26, 25, 53, CardinalLongitude.E, 13, 7, 38, CardinalLatitude.N);
var intersectCoord = intersectEllipsoid.IntersectOrthodromic(point1, point2, point3, point4);

Ответом будет точка — объект класса , в котором определены долгота и широта, в десятичных градуса (,) или в радианах (,).

Вычисление широты по долготе или долготы по широте

Для рассчётов, неоходимо создать объект класса .

var intermediateEllipsoid = new IntermediatePointService(new Ellipsoid());
var ntermediateSpheroid = new IntermediatePointService(new Spheroid());

Для вычисления широты вызвать метод , в который передать значение долготы, для которого мы вычисляем широту, и две координаты характеризующие ортодромию.

var coord1 = new Point(10, 10);
var coord2 = new Point(30, 50);
var lat = intermediateEllipsoid.GetLatitude(20, coord1, coord2);

Для вычисления долготы вызвать метод , в который передать значение широты, для которого мы вычисляем долготу, и две координаты характеризующие ортодромию.

var coord1 = new Point(10, 10);
var coord2 = new Point(30, 50);
var lat = intermediateEllipsoid.GetLongitude(20, coord1, coord2);

В обоих случая ответом будет значение типа .

[править] Пример программной реализации

Исходники вышеприведённых функций можно найти в архиве Sph.zip в файле 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 buf1024;
  double pt12, pt22;
  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 = Radians(lat1);
    pt11 = Radians(lon1);
    pt2 = Radians(lat2);
    pt21 = 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 ;
}

В архиве Sph.zip этот код находится в файле inv.c. Создадим исполняемый модуль inv компилятором gcc:

$ gcc -o inv inv.c sph.c -lm

Впрочем, в архиве есть Makefile. Для MS Windows готовую программу inv.exe можно найти в архиве Sph-win32.zip.

Программа читает данные из стандартного ввода консоли и отправляет результаты на стандартный вывод. Для чтения и записи файлов используются символы перенаправления потока «>» и «<» соответственно. Из каждой строки ввода программа считывает координаты двух точек φ₁, λ₁, φ₂, λ₂, которые должны быть в градусах, решает обратную задачу и записывает в строку вывода α₁, α₂, s (азимуты прямого и обратного направлений в градусах; расстояние между пунктами в километрах, а точнее, в единицах, определённых константой A_E).

Создадим файл inv.dat, содержащий одну строку данных:

30 0 52 54

После запуска программы

$ inv < inv.dat

получим α₁, α₂, s:

44.804060 262.415109 5001.1309

В архиве Sph-py.zip находятся скрипты на языке Питон. Выполнение скрипта в командной консоли:

$ python inv.py inv.dat

[править] Алгоритм

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

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

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

Можно устранить второй пункт, если в первом заменить долготу λ₂ на разность долгот (λ₂ − λ₁).

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

/*
 * Решение обратной геодезической задачи
 *
 * Аргументы исходные:
 *     pt1  - {широта, долгота} точки Q1
 *     pt2  - {широта, долгота} точки Q2
 *
 * Аргументы определяемые:
 *     azi  - азимут начального направления
 *     dist - расстояние (сферическое)
 */
void SphereInverse(double pt1, double pt2, double *azi, double *dist)
{
  double x3, pt2;
 
  SpherToCart(pt2, x);			// сферические -> декартовы
  Rotate(x, pt11, 2);			// первое вращение
  Rotate(x, M_PI_2 - pt1, 1);	// второе вращение
  CartToSpher(x, pt);	     		// декартовы -> сферические
  *azi = M_PI - pt1;
  *dist = M_PI_2 - pt;
 
  return;
}

Следует заметить, что прямая и обратная задача математически идентичны, и алгоритмы их решения зеркально отражают друг друга.

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

В данном случае в качестве сферических координат φ, λ подставим φ₂, λ₂.

Реализация на Си:

/*
 * Преобразование сферических координат в вектор
 *
 * Аргументы исходные:
 *     y - {широта, долгота}
 *
 * Аргументы определяемые:
 *     x - вектор {x, y, z}
 */
void SpherToCart(double y, double x)
{
  double p;
 
  p = cos(y);
  x2 = sin(y);
  x1 = p * sin(y1);
  x = p * cos(y1);
 
  return;
}

Вращение вокруг оси

Представим оператор вращения вокруг оси X на угол θ в следующем виде:

Операторы вращения вокруг осей Y и Z получаются перестановкой символов.

Реализация вращения вокруг i-ой координатной оси на Си:

/*
 * Вращение вокруг координатной оси
 *
 * Аргументы:
 *     x - входной/выходной 3-вектор
 *     a - угол вращения
 *     i - номер координатной оси (0..2)
 */
void Rotate(double x, double a, int i)
{
  double c, s, xj;
  int j, k;
 
  j = (i + 1) % 3;
  k = (i - 1) % 3;
  c = cos(a);
  s = sin(a);
  xj = xj * c + xk * s;
  xk = -xj * s + xk * c;
  xj = xj;
 
  return;
}

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

В данном случае в роли сферических координат φ, λ окажутся углы (90° − σ), (180° − α₁).

Реализация на Си:

/*
 * Преобразование вектора в сферические координаты
 *
 * Аргументы исходные:
 *     x - {x, y, z}
 *
 * Аргументы определяемые:
 *     y - {широта, долгота}
 *
 * Возвращает:
 *     длину вектора
 */
double CartToSpher(double x, double y)
{
  double p;
 
  p = hypot(x, x1);
  y1 = atan2(x1, x);
  y = atan2(x2, p);
 
  return hypot(p, x2);
}

[править] Алгоритм

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

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

Если третий пункт пропустить, на выходе вместо долготы λ₂ получится разность долгот (λ₂ − λ₁), что упростит алгоритм. Останется только прибавить долготу первого пункта. С другой строны, благодаря третьему пункту долгота λ₂ всегда будет в диапазоне .

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

/*
 * Решение прямой геодезической задачи
 *
 * Аргументы исходные:
 *     pt1  - {широта, долгота} точки Q1
 *     azi  - азимут начального направления
 *     dist - расстояние (сферическое)
 *
 * Аргументы определяемые:
 *     pt2  - {широта, долгота} точки Q2
 */
void SphereDirect(double pt1, double azi, double dist, double pt2)
{
  double pt2, x3;
 
  pt = M_PI_2 - dist;
  pt1 = M_PI - azi;
  SpherToCart(pt, x);			// сферические -> декартовы
  Rotate(x, pt1 - M_PI_2, 1);	// первое вращение
  Rotate(x, -pt11, 2);		// второе вращение
  CartToSpher(x, pt2);	     		// декартовы -> сферические
 
  return;
}

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

В данном случае в качестве сферических координат φ, λ подставим углы (90° − σ), (180° − α₁).

Реализация на Си:

/*
 * Преобразование сферических координат в вектор
 *
 * Аргументы исходные:
 *     y - {широта, долгота}
 *
 * Аргументы определяемые:
 *     x - вектор {x, y, z}
 */
void SpherToCart(double y, double x)
{
  double p;
 
  p = cos(y);
  x2 = sin(y);
  x1 = p * sin(y1);
  x = p * cos(y1);
 
  return;
}

Вращение вокруг оси

Представим оператор вращения вокруг оси X на угол θ в следующем виде:

Операторы вращения вокруг осей Y и Z получаются перестановкой символов.

Реализация вращения вокруг i-ой координатной оси на Си:

/*
 * Вращение вокруг координатной оси
 *
 * Аргументы:
 *     x - входной/выходной 3-вектор
 *     a - угол вращения
 *     i - номер координатной оси (0..2)
 */
void Rotate(double x, double a, int i)
{
  double c, s, xj;
  int j, k;
 
  j = (i + 1) % 3;
  k = (i - 1) % 3;
  c = cos(a);
  s = sin(a);
  xj = xj * c + xk * s;
  xk = -xj * s + xk * c;
  xj = xj;
 
  return;
}

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

В данном случае в роли сферических координат φ, λ окажутся φ₂, λ₂.

Реализация на Си:

/*
 * Преобразование вектора в сферические координаты
 *
 * Аргументы исходные:
 *     x - {x, y, z}
 *
 * Аргументы определяемые:
 *     y - {широта, долгота}
 *
 * Возвращает:
 *     длину вектора
 */
double CartToSpher(double x, double y)
{
  double p;
 
  p = hypot(x, x1);
  y1 = atan2(x1, x);
  y = atan2(x2, p);
 
  return hypot(p, x2);
}

Практическое применение вычисления дирекционного угла

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

Прокладка трасс – вычисление дирекционного угла позволяет определить направление по которому должна пролегать трасса дороги, трубопровода, электропередачи и других инженерных сооружений

Это особенно важно при проектировании новых объектов или модернизации существующих.

Определение координат земных объектов – при обновлении карт или создании новых картографических материалов необходимо знать координаты каждого объекта. Вычисление дирекционного угла помогает не только определить координаты объекта, но также его относительное расположение по отношению к другим объектам.

Навигация – вычисление дирекционного угла позволяет определить путь и направление движения, а также точное положение в пространстве

Это особенно важно для мореплавателей, пилотов, автомобилистов и других людей, работающих в сфере транспорта.

Геофизические исследования – в геофизических исследованиях, например, при поиске месторождений полезных ископаемых, вычисление дирекционного угла позволяет определить точное направление взаимодействия и получить максимально точные результаты.

Архитектура и строительство – при проектировании зданий, сооружений или жилых комплексов очень важно определить направление, размеры и расположение объектов. Вычисление дирекционного угла помогает сделать эти расчеты и построить объекты в точности согласно заранее заданным параметрам.

Это лишь некоторые примеры практического применения вычисления дирекционного угла. Благодаря точным расчетам и использованию современных геодезических инструментов, мы можем получить максимально точные результаты и эффективно применять их в различных областях деятельности.