Несколько месяцев назад я наткнулся на книгу [1] австрийского историка, спелиолога и преподавателя университета на пенсии Х. Куша и др., в которой он и его команда утверждают, что вместе обнаружили доисторические и частично выгравированные каменные плиты на австрийских раскопках. с другими предметами, такими как статуэтки, керамика и органический материал. На одной из этих пластин (Steintafel IV, стр. 179, рис. 252) изображен ряд отверстий, соединенных линиями — классическая фигурка из палочек звездного созвездия Тельца, кажется. Рисунок 1 представляет собой цифровую реконструкцию исходного изображения и показывает 13 отмеченных точек и белый кружок, некоторые из которых соединены белыми линиями.

Хотя эта и другие пластины содержат весьма любопытные дополнительные изображения, исследовательская группа настаивает на их подлинности. Я спрашивал себя, могу ли я сделать вывод по рисунку созвездия, подлинное оно или нет. Ведь линии классических и современных изображений созвездий Тельца, как видно на Рис. 2, лишь немного отклоняются от Рис. 1, что само по себе довольно любопытно и подгонка кажется «слишком хорошей, чтобы быть правдой».

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

Как лучше всего совместить фигурку с ее галактическим аналогом и найти дату или период, для которого это сопоставление подходит лучше всего?

Поиск подходящего инструмента для работы

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

  • translation: перемещение набора точек в другое положение x/y,
  • вращение: вращение заданной точки на угол до 360°,
  • масштабирование: уменьшение или расширение точки, заданной постоянным коэффициентом,
  • отражение: отражает каждую точку набора точек на заданной оси.

Особенно простой и эффективный алгоритм был разработан Кабшем (1976 и 1978 гг.) и позднее усовершенствован Умеямой (1991 г.). Как отмечают другие, такие как T. Siipola на своем веб-сайте [3], алгоритм использует разложение по сингулярным значениям (SVD) ковариационной матрицы двух наборов точек и получает оптимальный результат с помощью чисто алгебраических операций. В дополнение к исходному пространственному преобразованию Кабша, поддерживающему только перенос и вращение, Умеяма добавил поддержку коэффициента масштабирования и исправил случайное нежелательное отражение . Между прочим, Сиипола также применяет базовую механику алгоритма Умеямы к двум заданным звездным созвездиям.

Пейринг — это забота

Следующий вопрос касался проблемы определения пар совпадающих точек между каменной плитой и настоящим звездным созвездием. У некоторых звезд были вполне очевидные аналоги, но как насчет более неоднозначных пар, скажем, если мне придется выбирать из множества вариантов рядом со звездами? Что ж, отмеченные пятна на каменной пластине имели почти однозначные аналоги среди 14 ярчайших звезд по своей видимой величине. Только главная звезда Альдебаран, представленная большим кругом, должна была быть принята в центре этого круга. Более строгий статистический анализ показал, что предложенное отображение действительно является наиболее вероятным.

Как я могу выразить свою неуверенность в отношении конкретных пар?

Хотя алгоритм Умеямы довольно популярен, в нем отсутствует поддержка отдельных весов для выражения уверенности в определенных сопоставлениях. Однако алгоритм может быть адаптирован для учета индивидуальных весов, см., например, реализацию в библиотеке PyTorch3D [4].

Роль весов

На рисунках 3 и 4 показаны результаты картирования положения 14 самых ярких звезд созвездия Тельца, как они были видны с Земли в 6000 году до н. э. Гравировка на каменной плите показана серым цветом, а фактическое звездное созвездие - синим. В то время как левое изображение использует одинаковые веса для всех точек, правое изображение использует в 100 раз больший вес для двух выбранных точек, чтобы указать на субъективно более высокую достоверность их сопоставления.

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

Правильная мера ошибки

Результирующее преобразование, примененное к первому набору точек, минимизирует средневзвешенную квадратичную ошибку (wMSE) между наборами точек. Каждая точечная ошибка вносит свой вклад в общую среднеквадратичную ошибку своим весом. Если набор исходных точек может быть перемещен, повернут и масштабирован таким образом, что он полностью перекрывает набор целевых точек, то результирующая ошибка равна нулю.

Полученные результаты

Чтобы найти даты, для которых (равновзвешенное) наложение приводит к минимальным ошибкам, я рассчитал положение звезд для интервала в 200 000 лет (100 000 лет до н. программное обеспечение [2]. На рис. 5 показано, как звезды меняют свое положение в течение длительного периода времени. Кроме того, графики показывают наилучшее визуальное совпадение между 25000 г. до н.э. и 25000 г. н.э.

Обувь подходит?

Чтобы оценить качество наложения, я нормализовал возвращаемые среднеквадратические ошибки до интервала [0;1]. На рисунках 6aи 6b мы видим, что гравюра на камне показывает наименьшие ошибки для временного периода 7000 г. до н.э. — 3000 г. н.э., с минимальным получено для 2000 г. до н.э., см. рис. 7.

Самый низкий диапазон ошибок в 5 % считается оптимальным в данном контексте. Эксперименты с разными весами точек существенно не изменили основной результат, полученный при равных весах.

Это подлинное?

Анализ показал, что можно предположить, что выгравированная фигурка действительно представляет созвездие Тельца, учитывая качество совпадения, которое может быть достигнуто с помощью операции наложения. На рисунке скорее всего изображено звездное созвездие Тельца, которое можно было наблюдать с Земли в течение последних 10 000 лет, охватывающих период от неолита до современной эпохи.

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

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

Взвешенный алгоритм Умеямы

Даны два набора m-мерных точек размера n в матричной форме, расположенных по столбцам:

или, в векторной форме

и вектор весов

Примените веса к каждой точке в обоих наборах и определите взвешенные центроиды наборов точек:

и соответствующая взвешенная ковариационная матрица

Затем вычислите его сингулярное разложение

и исправить отражения/зеркалирование, так как они нежелательны для наших целей

Для масштабирования требуется взвешенная дисперсия A.

Теперь компоненты преобразования можно вычислить следующим образом:

Матрица поворота

Масштаб фактор

Перевод вектор

Теперь преобразованные точки можно вычислить следующим образом:

где

Минимизированную средневзвешенную ошибку теперь можно вычислить следующим образом:

Аналогичный вывод, но без поддержки веса, вместе с реализацией на Python можно найти в [3].

Пример реализации в R

# Note: A, B contain row vectors. w is a row vector.
umeyama_weighted <- function(A, B, w) {
  m <- dim(A)[1]
  n <- dim(A)[2]
  
  if (missing(w)) w <- rep(1, each=n)
  W <- diag(c(w))
  sw <- sum(W)
  
  # determine weighted centroids
  Amu <- rowSums(A %*% W) / sw
  Bmu <- rowSums(B %*% W) / sw
  Ac <- (A - Amu) %*% W
  Bc <- (B - Bmu) %*% W
  
  C <- Ac %*% t(Bc) / sw # weighted covariance matrix
  
  Z <- svd(C) # singular value decomp. of weighted covariance matrix
  U <- Z$u
  D <- Z$d
  V <- Z$v
  
  S <- diag(m)
  S[-1] <- S[-1] * det(U %*% t(V)) # correct for reflections
  
  R <- V %*% S %*% t(U) # rotation matrix
  
  varA <- sum(Ac^2) / sw
  scale <- 1/varA * sum(D %*% S) # determine scale factor
  translate <- Bmu - (scale * R %*% Amu) # determine transl. vector
  
  Aprime <- scale * R %*% A + c(translate) # new point
  # determine weighted mean squared error
  err <- 1/sw * sum((B - Aprime)^2 %*% W) 
  
  return(list(c=scale, R=R, t=translate, err=err, P=Aprime))
}

[1]: Куш, Х., и Куш, И. (2021). Geheime Unterwelt: Auf den Spuren von Jahrtausende alten unterirdischen Völkern. Стокер. ISBN 978–3853653234.

[2]: Зотти Г., Хоффманн С. М., Вольф А., Шеро Ф. и Шеро Г. (2021). Смоделированное небо: Stellarium для исследований в области культурной астрономии. Журнал Skyscape Archaeology, 6 (2), 221–258. https://doi.org/10.1558/jsa.17822

[3]: Siipola, T. (7 июня 2022 г.). Выравнивание шаблонов точек с помощью алгоритма Кабша–Умеямыhttps://zpl.fi/aligning-point-patterns-with-kabsch-umeyama-algorithm/

[4] PyTorch3D (7 июня 2022 г.). Исходный код для pytorch3d.ops.points_alignment https://pytorch3d.readthedocs.io/en/latest/_modules/pytorch3d/ops/points_alignment.html