Геодезические системы пространственных координат

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

Рассматриваются преобразования между пространственными координатными системами. Приводится пример программной реализации на языке Питон.

Земной эллипсоид

Земным эллипсоидом называется эллипсоид вращения, поверхность которого по форме и размерам довольно близка к поверхности геоида.

Поверхность эллипсоида образуется вращением эллипса вокруг его малой оси.

Эллипс обычно определяется размером его большой полуоси a и сжатием f. Реже вместо сжатия задаётся размер малой полуоси b:

$f={\frac {a-b}{a}}\qquad b=a(1-f)$

В теории и практике вычислений широко используются такие параметры, как полярный радиус кривизны поверхности c, первый эксцентриситет e и второй эксцентриситет e′:

$c={\frac {a}{1-f}}\qquad e^{2}=f(2-f)\qquad e'^{2}={\frac {e^{2}}{1-e^{2}}}$

Пример функции Питона, вычисляющей по a и f параметры b, c, e и e′:

def initSpher(a, f):
    b = a * (1. - f)
    c = a / (1. - f)
    e2 = f * (2. - f)
    e12 = e2 / (1. - e2)
    return (b, c, e2, e12)

Системы координат

Рассмотрим следующие системы координат.

1. Геоцентрические декартовы прямоугольные координаты:
   • начало координат находится в центре эллипсоида,
   • ось z расположена вдоль оси вращения эллипсоида и направлена в северный полюс,
   • ось x лежит в пересечении экватора и начального меридиана,
   • ось y лежит в пересечении экватора и меридиана с долготой L = 90°.
2. Система геодезических координат:

   геодезическая широта B

   угол между нормалью к поверхности эллипсоида и плоскостью экватора,

   геодезическая долгота L

   двугранный угол между плоскостями данного и начального меридианов,

   геодезическая высота H

   кратчайшее расстояние до поверхности эллипсоида.

3. Топоцентрические декартовы прямоугольные координаты:
   • начало координат находится в некоторой точке Q₀ (B₀, L₀, H₀) над эллипсоидом,
   • ось z расположена вдоль нормали к поверхности эллипсоида и направлена вверх,
   • ось x расположена в плоскости меридиана и направлена на север,
   • ось y перпендикулярна к осям x и z и направлена на восток.

Преобразования координат

Переход от геодезических координат к геоцентрическим

Это преобразование выполняется по следующим формулам:

${\begin{array}{rcl}x&=&(N+H)\cos B\cos L\\y&=&(N+H)\cos B\sin L\\z&=&(N+H-e^{2}N)\sin B\end{array}}$

Здесь N — так называемый радиус кривизны первого вертикала:

$N={\frac {a}{1-e^{2}\sin ^{2}B}}={\frac {c}{1+e'^{2}\cos ^{2}B}}$

Реализация на Питоне:

def fromLatLong(lat, lon, h, a, f):
    b, c, e2, e12 = initSpher(a, f)
    cos_lat = math.cos(lat)
    n = c / math.sqrt(1. + e12 * cos_lat ** 2)
    p = (n + h) * cos_lat
    x = p * math.cos(lon)
    y = p * math.sin(lon)
    z = (n + h - e2 * n) * math.sin(lat)
    return (x, y, z)

Переход от геоцентрических координат к геодезическим

Проще всего вычисляется долгота:

$\operatorname {tg\,} L={\frac {y}{x}}$

Сложнее с определением широты и высоты. Существует множество способов решения этой непростой задачи. Воспользуемся итеративным методом Боуринга.

В начале находится предварительная оценка параметра θ:

$\operatorname {tg\,} \theta ={\frac {z}{r}}\left(1+e'^{2}{\frac {b}{\rho }}\right)(1-f)$

Здесь ρ — геоцентрический радиус-вектор, r — расстояние от оси вращения эллипсоида:

$\rho ={\sqrt {x^{2}+y^{2}+z^{2}}}\qquad r={\sqrt {x^{2}+y^{2}}}$

Затем вычисляется тангенс широты:

$\operatorname {tg\,} B={\frac {z+e'^{2}b\sin ^{3}\theta }{p-e^{2}a\cos ^{3}\theta }}$

и параметр θ уточняется:

$\operatorname {tg\,} \theta =(1-f)\operatorname {tg\,} B$

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

В конце определяется высота:

$H={\frac {r}{\cos B}}-N={\frac {z}{\sin B}}-N(1-e^{2})$

def toLatLong(x, y, z, a, f):
    b, c, e2, e12 = initSpher(a, f)
    p = math.hypot(x, y)
    if p == 0.:
        lat = math.copysign(math.pi / 2., z)
        lon = 0.
        h = math.fabs(z) - b
    else:
        t = z / p * (1. + e12 * b / math.hypot(p, z))
        for i in range(2):
            t = t * (1. - f)
            lat = math.atan(t)
            cos_lat = math.cos(lat)
            sin_lat = math.sin(lat)
            t = (z + e12 * b * sin_lat ** 3) / (p - e2 * a * cos_lat ** 3)
        lon = math.atan2(y, x)
        lat = math.atan(t)
        cos_lat = math.cos(lat)
        n = c / math.sqrt(1. + e12 * cos_lat ** 2)
        if math.fabs(t) <= 1.:
            h = p / cos_lat - n
        else:
            h = z / math.sin(lat) - n * (1. - e2)
    return (lat, lon, h)

Переход от геодезических координат к топоцентрическим. Обратная пространственная задача

Постановка задачи: начало топоцентрической системы координат задано точкой Q₀ (B₀, L₀, H₀); по геодезическим координатам точки Q (B, L, H) вычислить её топоцентрические координаты x, y, z.

Конформное преобразование между двумя декартовыми прямоугольными системами координат всегда может быть представлено последовательностью сдвигов и вращений координатной системы. Данное преобразование можно реализовать по следующему алгоритму:

• по геодезическим координатам точки B, L, H вычислить её геоцентрические координаты x, y, z,
• сместить начало координат вдоль оси z на величину e² N₀ sin B₀ (т.е. до вершины конуса, образованного нормалями, лежащими на параллели с широтой B₀),
• повернуть систему координат вокруг оси z на угол L₀, чтобы ось x оказалась в плоскости меридиана точки Q₀,
• повернуть систему координат вокруг оси y на угол 90° − B₀, чтобы ось z совпала с нормалью к поверхности эллипсоида в точке Q₀,
• сместить начало координат вдоль оси z на величину N₀ + H₀,
• изменить знак x на противоположный.

Реализация алгоритма:

def inverse3d(lat0, lon0, h0, lat, lon, h, a, f):
    b, c, e2, e12 = initSpher(a, f)
    sin_lat = math.sin(lat0)
    n = a / math.sqrt(1. - e2 * sin_lat ** 2)
    x, y, z = fromLatLong(lat, lon, h, a, f)
    z = z + e2 * n * sin_lat
    x, y = rotate(x, y, lon0)
    z, x = rotate(z, x, math.pi / 2. - lat0)
    z = z - (n + h0)
    x = -x
    return (x, y, z)

Для реализации алгоритма понадобится функция вращения:

def rotate(x, y, a):
    c, s = math.cos(a), math.sin(a)
    return (x * c + y * s, -x * s + y * c)

Рассмотренная задача является разновидностью обратной геодезической задачей в пространстве. Вместо декартовых прямоугольных топоцентрических координат может требоваться вычисление каких-то других топоцентрических величин, например, полярных координат «дальность-азимут-зенитное расстояние», варианты могут быть разные. Однако в большинстве случаев сначала находятся x, y, z, по которым и находятся искомые значения.

Программы

import math

def rotate(x, y, a):
    c, s = math.cos(a), math.sin(a)
    return (x * c + y * s, -x * s + y * c)

def forward3d(lat1, lon1, h1, x, y, z, a, f):
    b, c, e2, e12 = initSpher(a, f)
    sin_lat = math.sin(lat1)
    n = a / math.sqrt(1. - e2 * sin_lat ** 2)
    x = -x
    z = z + (n + h1)
    z, x = rotate(z, x, lat1 - math.pi / 2.)
    x, y = rotate(x, y, -lon1)
    z = z - e2 * n * sin_lat
    return toLatLong(x, y, z, a, f)

def inverse3d(lat1, lon1, h1, lat2, lon2, h2, a, f):
    b, c, e2, e12 = initSpher(a, f)
    sin_lat = math.sin(lat1)
    n = a / math.sqrt(1. - e2 * sin_lat ** 2)
    x, y, z = fromLatLong(lat2, lon2, h2, a, f)
    z = z + e2 * n * sin_lat
    x, y = rotate(x, y, lon1)
    z, x = rotate(z, x, math.pi / 2. - lat1)
    z = z - (n + h1)
    x = -x
    return (x, y, z)

# -*- coding: utf-8 -*-
from sys import argv
import math
import latlong

script, fn = argv

a, f = 6378136., 1./298.25784 # ПЗ-90

lat0, lon0, hgt0 = math.radians(65.), math.radians(45.), 500.

fp = open(fn, 'r')
for line in fp:
    x, y, z = map(float, line.split(" "))
    lat, lon, hgt = latlong.forward3d(lat0, lon0, hgt0, x, y, z, a, f)
    print "%.8f %.8f %.3f" % (math.degrees(lat), math.degrees(lon), hgt)
fp.close()

# -*- coding: utf-8 -*-
from sys import argv
import math
import latlong

script, fn = argv

a, f = 6378136., 1./298.25784 # ПЗ-90

lat0, lon0, hgt0 = math.radians(65.), math.radians(45.), 500.

fp = open(fn, 'r')
for line in fp:
    lat, lon, hgt = map(float, line.split(" "))
    lat = math.radians(lat)
    lon = math.radians(lon)
    x, y, z = latlong.inverse3d(lat0, lon0, hgt0, lat, lon, hgt, a, f)
    print "%.3f %.3f %.3f" % (x, y, z)
fp.close()

Расчёты

-40000 30000 0

64.63992461 45.62743323 695.578

-40000.000 30000.000 0.000

Ссылки

• Robert Burtch, A Comparison of Methods Used in Rectangular to Geodetic Coordinate Transformations, 2006