Алгоритм восстановления координат и высотной привязки исторической топографической карты (с учётом того, что доступны ориентиры и несколько контрольных точек с современными координатами):
1) Подготовка изображения
- Отсканировать карту в высоком разрешении, сохранить геометрически неизменённой (без обрезки/растяжения).
- Определить координаты пикселей опорных точек на растровом изображении $(x_i,y_i)$.
2) Подбор и проверка опорных точек (GCP)
- Выбрать устойчивые ориентиры (перекрёстки дорог, углы зданий, церкви, мосты).
- Записать соответствующие современные координаты $(X_i,Y_i)$ (в одной системе координат).
- Проверить соответствие и метки (минимум >2–6, см. ниже). Исключить очевидные ошибочные GCP.
3) Выбор модели преобразования (по числу и характеру искажений)
- Ригидное (перенос+поворот+масштаб, 2D Helmert / similarity) — параметров 4. Минимум 2 неколлинеарных GCP. Формула:
$$X=s(\cos \theta x-\sin \theta y)+t_x,Y=s(\sin \theta x+\cos \theta y)+t_y.$$ - Аффинное (линейное искажение, сдвиг/скейл/сдвиг по осям/сдвиг по углу) — 6 параметров. Минимум 3 GCP:
$$X=a_{1}x+a_{2}y+a_3,Y=b_{1}x+b_{2}y+b_3.$$ Решается как линейная СЛАУ: собрать матрицу $A$ из $(x_i,y_i,1)$ и вектор наблюдений $L$.
- Гомография (проективное преобразование, перспективные искажения) — 8 параметров. Минимум 4 GCP.
- Нелинейные (полиномиум 2-го порядка, thin‑plate spline, радиальные базисы) — для учёта локальных деформаций; требует больше GCP (полином 2-го порядка по каждой координате — по 6 коэффициентов → минимум 6 GCP).
- Практика: начать с similarity → affine → homography/TPS, смотреть остатки и распределение векторов сдвига.
4) Построение и решение системы (метод наименьших квадратов)
- Для аффинной: для всех точек собрать уравнения
$$\{\begin{array}{l}{X}_i=a_1{x}_i+a_2{y}_i+a_3,\\ Y_i=b_1{x}_i+b_2{y}_i+b_3\end{array}$$ и решить $\mathbf{Ap}=\mathbf{L}$ методом наименьших квадратов (SVD).
- Для homography: решить нелинейную систему через нормализацию и линейное решение/итерации.
- Для similarity: можно аналитически оценить $s,\theta ,t_x,t_y$ или решать как линейную задачу после линейной параметризации.
5) Оценка качества и фильтрация выбросов
- Вычислить остатки $\Delta X_i=X_i-{\hat{X}}_i,\Delta Y_i=Y_i-{\hat{Y}}_i$.
- RMSE: $RMSE=\sqrt{\frac{1}{n}\sum _{i=1}^n((\Delta X_i{)}^2+(\Delta Y_i{)}^2)}$.
- Визуализировать векторы остатков; при больших локальных выбросах применить RANSAC для устойчивой оценки параметров или удалить проблемные GCP.
6) Применение трансформации к растровому изображению
- Использовать полученную модель для прямого преобразования пикселей в координаты и/или создать геореференцированный растр (например GDAL: задать GCP и выполнить warp).
- Для нелинейных деформаций использовать интерполяцию (cubic, TPS) при резэмплинге.
7) Высотная привязка
- Варианты в зависимости от наличия высот на карте:
a) Если карта содержит подписи высот отдельных точек/вершин контуров: сопоставить эти значения с современными точками $(Z_i)$. Подгонка линейного сдвига/масштаба:
$$Z=\alpha z_{map}+\beta$$ оценить $\alpha ,\beta$ по МНК при наличии подписанных высот.
b) Если имеются контурные линии (полилинии) с известными уровнями — создать цифровую модель рельефа (TIN/интерполяция): размечать линии по пикселям → векторизовать → построить DEM (методы: TIN, IDW, Kriging, B-spline).
c) Если есть лишь отдельные контрольные высоты $Z_i$ — интерполировать высоту в произвольных точках: например IDW
$$Z(p)=\frac{{\sum }_i{w}_i{Z}_i}{{\sum }_i{w}_i},w_i=\frac{1}{d(p,p_i{)}^p},$$ где обычно $p=2$.
- Учесть вертикальные системы (геоид/эллипсоид) — выполнить преобразование датумов при необходимости.
8) Проверка итоговой привязки
- Контроль: независимые контрольные точки (не использованные при построении) для валидации.
- Отчёт: RMSE по плану и высоте, карты векторов сдвига, графики остатков.
9) Рекомендации практические
- Минимум GCP: similarity — 2, affine — 3, homography — 4, полином 2-го порядка — 6; для TPS — желательно >10–15 ровно распределённых GCP.
- Если контрольных точек мало и карта сильно деформирована — применять локальные привязки: разбивать карту на пласты и гео-регионально выравнивать.
- Использовать инструменты: QGIS (Georeferencer), GDAL (gdal_translate + gdalwarp), специализированные пакеты (OpenCV для оценки матриц трансформации, SAGA/GRASS для DEM/контуров).
- Документировать все шаги, GCP и метрики качества.
Кратко: подобрать и проверить GCP → выбрать модель трансформации по характеру искажений → решить параметры МНК/устойчивым методом → оценить и удалить выбросы → применить трансформацию к растру → получить высоты через явные подписи/векторные контуры или интерполяцию, учитывая вертикальный датум.