Вычисление пересечения двух кругов?

29

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

Например, учитывая:

  • Широта / долгота (37,673442, -90,234036) радиус 107,5 нм
  • Шир. / Долг. (36.109997, -90.953669) Радиус 145 НМ

Я должен найти две точки пересечения с одной из них (36.948, -088.158).

Было бы тривиально легко решить это на плоской плоскости, но у меня нет опыта решения уравнений на несовершенной сфере, такой как поверхность Земли.

Будет
источник
1
Если все ваши радиусы будут такими маленькими (менее нескольких километров), тогда Земля в этом масштабе практически плоская, и вы можете также выбрать точную, простую проекцию и выполнить обычные евклидовы вычисления. Удостоверьтесь, что вы вычисляете пересечение более чем на три знака после запятой - неточность в третьем знаке после запятой настолько же велика, что и любой из ваших радиусов!
whuber
1
Я должен был добавить единицы измерения, эти радиусы указаны в м. Миль, так что это все еще небольшое расстояние относительно поверхности Земли, но больше, чем несколько километров. Как этот масштаб влияет на искажение? Я пытаюсь найти решение с точностью менее 1 нм, поэтому оно не должно быть сверхточным. Благодарность!
Уилл
Это все хорошо знать, потому что это показывает, что вы можете использовать сферическую модель Земли - более сложные эллипсоидальные модели не нужны.
whuber
@whuber Означает ли это, что проблему можно переформулировать так: найти пересечение трех сфер, где одна из сфер является землей, а две другие центрированы в точках с их соответствующими радиусами?
Кирк Куйкендалл
@ Кирк Да, это способ сделать это, предполагая сферическую модель земной поверхности. После некоторых предварительных расчетов это сводится к частному случаю задачи трилатерации в 3D. (Расчеты необходимы для преобразования расстояния вдоль сферических дуг в расстояния вдоль сферических хорд, которые становятся радиусами двух меньших сфер.)
whuber

Ответы:

21

Это не намного сложнее в сфере, чем в самолете, как только вы поймете, что

  1. Рассматриваемые точки - это взаимные пересечения трех сфер: сфера с центром под местоположением x1 (на поверхности земли) данного радиуса, сфера с центром под местоположением x2 (на поверхности земли) данного радиуса и сама земля , которая представляет собой сферу с центром в O = (0,0,0) данного радиуса.

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

Следовательно, задача сводится к пересечению линии со сферой, что легко.


Вот подробности. Входными данными являются точки P1 = (lat1, lon1) и P2 = (lat2, lon2) на поверхности земли, рассматриваемые как сфера, и два соответствующих радиуса r1 и r2.

  1. Преобразовать (широта, долгота) в (x, y, z) геоцентрические координаты. Как обычно, потому что мы можем выбрать единицы измерения, в которых Земля имеет единичный радиус,

    x = cos(lon) cos(lat)
    y = sin(lon) cos(lat)
    z = sin(lat).
    

    В этом примере P1 = (-90,234036 градусов, 37,673442 градусов) имеет геоцентрические координаты x1 = (-0,00323306, -0,7915, 0,61116) и P2 = (-90,953669 градусов, 36,109997 градусов) имеет геоцентрические координаты x2 = (-0,0134475, -0. 0,589337).

  2. Преобразуйте радиусы r1 и r2 (которые измеряются вдоль сферы) в углы вдоль сферы. По определению, одна морская миля (NM) составляет 1/60 градуса дуги (что равно pi / 180 * 1/60 = 0,0002908888 радиан). Поэтому, как углы,

    r1 = 107.5 / 60 Degree = 0.0312705 radian
    r2 = 145 / 60 Degree = 0.0421788 radian
    
  3. Геодезический круг радиуса r1 вокруг x1 является пересечением поверхности Земли с евклидовой сферы радиуса греха (r1) с центром в Cos (r1) * x1.

  4. Плоскость, определяемая пересечением сферы радиуса sin (r1) вокруг cos (r1) * x1 и земной поверхности, перпендикулярна x1 и проходит через точку cos (r1) x1, откуда ее уравнение равно x.x1 = cos (r1) («.» представляет обычное скалярное произведение ); аналогично для другого самолета. На пересечении этих двух плоскостей будет уникальная точка x0, представляющая собой линейную комбинацию x1 и x2. Запись x0 = a x1 + b * x2, два плоских уравнения

    cos(r1) = x.x1 = (a*x1 + b*x2).x1 = a + b*(x2.x1)
    cos(r2) = x.x2 = (a*x1 + b*x2).x2 = a*(x1.x2) + b
    

    Используя тот факт, что x2.x1 = x1.x2, который я напишу как q, решение (если оно существует) дается

    a = (cos(r1) - cos(r2)*q) / (1 - q^2),
    b = (cos(r2) - cos(r1)*q) / (1 - q^2).
    

    В текущем примере я вычисляю a = 0,973503 и b = 0,0260194.

    Очевидно, нам нужно q ^ 2! = 1. Это означает, что x1 и x2 не могут быть ни одной, ни антиподальной точкой.

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

    n = x1~Cross~x2
    

    выполняет задание, если n отлично от нуля: еще раз, это означает, что x1 и x2 не совпадают и не диаметрально противоположны. (Мы должны позаботиться о том, чтобы вычислить перекрестное произведение с высокой точностью, поскольку оно включает в себя вычитания с большим количеством отмен, когда x1 и x2 находятся близко друг к другу.) В этом примере n = (0.0272194, -0.00631254, -0.00803124) ,

  6. Поэтому мы ищем до двух точек вида x0 + t * n, которые лежат на поверхности земли: их длина равна 1. Эквивалентно, их квадратная длина равна 1:

    1 = squared length = (x0 + t*n).(x0 + t*n) = x0.x0 + 2t*x0.n + t^2*n.n = x0.x0 + t^2*n.n
    

    Член с x0.n исчезает, потому что x0 (будучи линейной комбинацией x1 и x2) перпендикулярен n. Два решения легко

    t = sqrt((1 - x0.x0)/n.n)
    

    и его негатив. Еще раз требуется высокая точность, потому что когда x1 и x2 близки, x0.x0 очень близко к 1, что приводит к некоторой потере точности с плавающей запятой. В этом примере t = 1,07509 или t = -1,07509. Поэтому две точки пересечения равны

    x0 + t*n = (0.0257661, -0.798332, 0.601666)
    x0 - t*n = (-0.0327606, -0.784759, 0.618935)
    
  7. Наконец, мы можем преобразовать эти решения обратно в (широта, долгота), преобразовав геоцентрические (x, y, z) в географические координаты:

    lon = ArcTan(x,y)
    lat = ArcTan(Sqrt[x^2+y^2], z)
    

    Для долготы используйте обобщенный арктангенс, возвращающий значения в диапазоне от -180 до 180 градусов (в вычислительных приложениях эта функция принимает в качестве аргументов как x, так и y, а не только отношение y / x; иногда ее называют «ATan2»).

    Я получаю два решения (-88.151426, 36.989311) и (-92.390485, 38.238380), показанные на рисунке желтыми точками.

3D фигура

На осях отображаются геоцентрические (x, y, z) координаты. Серый участок - это часть земной поверхности от -95 до -87 градусов по долготе, от 33 до 40 градусов по широте (отмечена одной градусной сеткой). Земная поверхность была сделана частично прозрачной, чтобы показать все три сферы. Правильность вычисленных решений очевидна тем, как желтые точки находятся на пересечениях сфер.

Whuber
источник
Билл, это круто. Одно пояснение, которое вы можете добавить, основываясь на том, кто пытался его реализовать. На шаге 2 вы не даете явное преобразование из градусов в радианы.
Джерси Энди
@ Джерси Спасибо за предложенное редактирование. Я немного изменил его, чтобы избежать избыточности и сохранить формулы как можно более четкими. Прочитав ветку, на которую вы ссылаетесь, я также вставил ссылку, чтобы объяснить точечный продукт.
whuber
8

Эллипсоидальной случай:

Эта проблема является обобщением проблемы нахождения морских границ, определенных как «срединные линии», и по этой теме имеется обширная литература. Мое решение этой проблемы заключается в использовании эквидистантной азимутальной проекции:

  1. Угадай в точке пересечения
  2. Спроецируйте две базовые точки, используя эту предполагаемую точку пересечения в качестве центра равноудаленной азимутальной проекции,
  3. Решите проблему пересечения в 2-мерном проецируемом пространстве.
  4. Если новая точка пересечения находится слишком далеко от старой, вернитесь к шагу 2.

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

Формулы приведены в разделе 14 геодезических на эллипсоиде вращения . Эллипсоидальная эквидистантная азимутальная проекция предоставлена GeographicLib . Версия MATLAB доступна в геодезических проекциях для эллипсоида .

cffk
источник
+1 Это удивительная статья: ваше скромное описание здесь не оправдывает себя.
whuber
См. Также мою более короткую статью о геодезических "Алгоритмы для геодезических" dx.doi.org/10.1007/s00190-012-0578-z (скачать бесплатно!), А также исправления и дополнения для этих статей geographiclib.sf.net/geod-addenda.html
cffk
1

Вот некоторый код R, чтобы сделать это:

p1 <- cbind(-90.234036, 37.673442) 
p2 <- cbind(-90.953669, 36.109997 )

library(geosphere)
steps <- seq(0, 360, 0.1)
c1 <- destPoint(p1, steps, 107.5 * 1852)
c2 <- destPoint(p2, steps, 145 * 1852)

library(raster)
s1 <- spLines(c1)
s2 <- spLines(c2)

i <- intersect(s1, s2)
coordinates(i)

#        x        y
# -92.38241 38.24267
# -88.15830 36.98740

s <- bind(s1, s2)
crs(s) <- "+proj=longlat +datum=WGS84"
plot(s)
points(i, col='red', pch=20, cex=2)
Роберт Хейманс
источник
1

Исходя из ответа @ whuber , вот код Java, который полезен по двум причинам:

  • в нем рассказывается о ArcTan (для Java и, возможно, для других языков?)
  • он обрабатывает возможные крайние случаи, включая случай, не упомянутый в ответе @ whuber.

Он не оптимизирован и не завершен (я пропустил такие очевидные классы, как Point), но должен сделать свое дело.

public static List<Point> intersection(EarthSurfaceCircle c1, EarthSurfaceCircle c2) {

    List<Point> intersections = new ArrayList<Point>();

    // project to (x,y,z) with unit radius
    UnitVector x1 = UnitVector.toPlanar(c1.lat, c1.lon);
    UnitVector x2 = UnitVector.toPlanar(c2.lat, c2.lon);

    // convert radii to radians:
    double r1 = c1.radius / RadiusEarth;
    double r2 = c2.radius / RadiusEarth;

    // compute the unique point x0
    double q = UnitVector.dot(x1, x2);
    double q2 = q * q;
    if (q2 == 1) {
        // no solution: circle centers are either the same or antipodal
        return intersections;
    }
    double a = (Math.cos(r1) - q * Math.cos(r2)) / (1 - q2);
    double b = (Math.cos(r2) - q * Math.cos(r1)) / (1 - q2);
    UnitVector x0 = UnitVector.add(UnitVector.scale(x1, a), UnitVector.scale(x2, b));

    // we only have a solution if x0 is within the sphere - if not,
    // the circles are not touching.
    double x02 = UnitVector.dot(x0, x0);
    if (x02 > 1) {
        // no solution: circles not touching
        return intersections;
    }

    // get the normal vector:
    UnitVector n = UnitVector.cross(x1, x2);
    double n2 = UnitVector.dot(n, n);
    if (n2 == 0) {
        // no solution: circle centers are either the same or antipodal
        return intersections;
    }

    // find intersections:
    double t = Math.sqrt((1 - UnitVector.dot(x0, x0)) / n2);
    intersections.add(UnitVector.toPolar(UnitVector.add(x0, UnitVector.scale(n, t))));
    if (t > 0) {
        // there's only multiple solutions if t > 0
        intersections.add(UnitVector.toPolar(UnitVector.add(x0, UnitVector.scale(n, -t))));
    }
    return intersections;
}

Также, что важно, обратите внимание на использование atan2- это обратное тому, что вы ожидаете от ответа @ whuber (я не знаю почему, но это работает):

    public static Point toPolar(UnitVector a) {
        return new Point(
                Math.toDegrees(Math.atan2(a.z, Math.sqrt(a.x * a.x + a.y * a.y))),
                Math.toDegrees(Math.atan2(a.y, a.x)));          
    }
ianhoolihan
источник
0

Рабочий код 'R' для ответа @wuhber.

P1 <- c(37.673442, -90.234036)
P2 <- c(36.109997, -90.953669) 

#1 NM nautical-mile is 1852 meters
R1 <- 107.5
R2 <- 145

x1 <- c(
  cos(deg2rad(P1[2])) * cos(deg2rad(P1[1])),  
  sin(deg2rad(P1[2])) * cos(deg2rad(P1[1])),
  sin(deg2rad(P1[1]))
);

x2 <- c(
  cos(deg2rad(P2[2])) * cos(deg2rad(P2[1])),
  sin(deg2rad(P2[2])) * cos(deg2rad(P2[1])),
  sin(deg2rad(P2[1]))
);

r1 = R1 * (pi/180) * (1/60)
r2 = R2 * (pi/180) * (1/60)

q = dot(x1,x2)
a = (cos(r1) - cos(r2) * q) / (1 - q^2)
b = (cos(r2) - cos(r1) * q)/ (1 - q^2)

n <- cross(x1,x2)

x0 = a*x1 + b*x2


t = sqrt((1 - dot(x0, x0))/dot(n,n))

point1 = x0 + (t * n)
point2 = x0 - (t * n)

lat1 = rad2deg(atan2(point1[2] ,point1[1]))
lon1= rad2deg(asin(point1[3]))
paste(lat1, lon1, sep=",")

lat2 = rad2deg(atan2(point2[2] ,point2[1]))
lon2 = rad2deg(asin(point2[3]))
paste(lat2, lon2, sep=",")
Шри
источник
-1

Если один из кругов - это Нортстар, то с юнит-сферой есть самый простой способ.

Вы можете измерить свою широту с Nortstar. Тогда у вас есть относительная позиция в этой сфере. v1 (0, sin (la), cos (la)) Вы знаете положение (угол) другой звезды (star2) из ​​альманаха. v2 (sin (lo2) * cos (la2), sin (la2), cos (lo2) * cos (la2)) Его векторы. Из уравнения сферы.

lo2 - относительная долгота. Это неизвестно .

Угол между вами и звездой2 вы тоже можете измерить, (м) И вы знаете, что внутренним продуктом двух единичного вектора является cos (угол) между ними. cos (m) = точка (v1, v2) теперь вы можете рассчитать относительную долготу (lo2). lo2 = экос ((соз (м) -sin (ли) * sin (ла2)) / (соз (ли) * соз (ла2)))

В конце концов вы добавляете настоящую долготу звезды2 к lo2. (или sub, зависит от его на западной стороне от вас, или на востоке.) lo2 теперь ваша долгота.

Извините за мой английский, я никогда не изучаю этот язык.


2 вещи: Северная звезда означает Полярную звезду.

Другая. Поскольку угол измеряется относительно горизонтали, всегда требуется коррекция на 90 градусов. Это справедливо и для угла м.

PS: реальный угол означает: положение звезды - время коррекции.

docwho
источник
Не ясно, как это отвечает на вопрос.
whuber