Конформное преобразование

2020-03-22T13:57:07Z

ArmRus: опечатка названия столбца расчетов

{{Статья|Опубликована|helmert2d}}

{{Аннотация|Конформное преобразование на плоскости широко используется в геодезии при создании местных координатных систем на небольшие территории, ограниченные размерами населённого пункта.}}

== Введение ==

Следующие формулы связывают координаты точек ''x'', ''y'', заданные в местной системе координат (МСК), и координаты ''X'', ''Y'', заданные в государственной системе координат (ГСК):

<math>\begin{align}
X & = X_0 + m \big[ (x - x_0) \cos \theta - (y - y_0) \sin \theta \big] \\
Y & = Y_0 + m \big[ (x - x_0) \sin \theta + (y - y_0) \cos \theta \big]
\end{align}</math>

где ''m'' — масштабный множитель, ''θ'' — угол разворота, ''X''0, ''Y''0, ''x''0, ''y''0 — координаты одного из геодезических пунктов в ГСК и МСК, как правило. Этот набор параметров называется «ключ».

Исходный материал для определения ключа — пары координат пунктов геодезической сети, полученные из независимого уравнивания одних и тех же измерений в МСК и в ГСК<ref>ГКИНП (ОНТА)-01-271-03 Руководство по созданию и реконструкции городских геодезических сетей с использованием спутниковых систем ГЛОНАСС/GPS.</ref>. В зависимости от класса пунктам (вернее, соответствующим парам уравнений) назначаются веса ''p''.

== Алгоритм нахождения параметров ==

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

<math>\begin{align}
X' & = a_0 + a_1 x - b_1 y \\
Y' & = b_0 + b_1 x + a_1 y
\end{align}</math>

Определению подлежат четыре параметра: ''a''0, ''b''0, ''a''1, ''b''1,.

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

=== Шаг 1: вычисление взвешенных средних ===

<math>
x_c = \frac{\sum\limits_{i=1}^n p_i x_i}{\sum\limits_{i=1}^n p_i} \quad
y_c = \frac{\sum\limits_{i=1}^n p_i y_i}{\sum\limits_{i=1}^n p_i} \quad
X_c = \frac{\sum\limits_{i=1}^n p_i X_i}{\sum\limits_{i=1}^n p_i} \quad
Y_c = \frac{\sum\limits_{i=1}^n p_i Y_i}{\sum\limits_{i=1}^n p_i}
</math>

=== Шаг 2: перенос осей в центр масс ===

<math>
\Delta x_i = x_i - x_c \quad
\Delta y_i = y_i - y_c \quad
\Delta X_i = X_i - X_c \quad
\Delta Y_i = Y_i - Y_c
</math>

=== Шаг 3: вычисление ''a''1 и ''b''1 ===

{| border="0" cellspacing="16"
|-
| <math>S_1 = \sum\limits_{i=1}^n p_i \Delta X_i \Delta x_i</math>
| <math>S_2 = \sum\limits_{i=1}^n p_i \Delta Y_i \Delta y_i</math>
|-
| <math>S_3 = \sum\limits_{i=1}^n p_i \Delta Y_i \Delta x_i</math>
| <math>S_4 = \sum\limits_{i=1}^n p_i \Delta X_i \Delta y_i</math>
|-
| colspan="2" align="center" | <math>S_5 = \sum\limits_{i=1}^n p_i \left( \Delta x_i^2 + \Delta y_i^2 \right)</math>
|-
| <math>a_1 = \frac{S_1 + S_2}{S_5}</math>
| <math>b_1 = \frac{S_3 - S_4}{S_5}</math>
|}

=== Шаг 4: вычисление ''a''0 и ''b''0 ===

<math>\begin{align}
a_0 & = X_c - a_1 x_c + b_1 y_c \\
b_0 & = Y_c - b_1 x_c - a_1 y_c
\end{align}</math>

=== Шаг 5: вычисление невязок ===

<math>
v_{xi} = X_i - X_i' \quad
v_{yi} = Y_i - Y_i'
</math>

Невязки позволяют выявить точки, плохо укладывающиеся в полученную модель. Классическая процедура удаляет такие «отлетающие» точки, после чего вычисление параметров повторяется без них. Робастные процедуры переназначают веса уравнениям в соответствии с невязками, и повторное вычисление повторяется с полным набором точек при том, что «отлетающие» влияют на результат незначительно.

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

=== Шаг 6: вычисление ключа ===

Вычислим масштабный множитель и угол разворота:

<math>
m = \sqrt{a_1^2 + b_1^2} \qquad
\theta = \operatorname{arc\,tg} \frac{b_1}{a_1}
</math>

Выберем ''j''-й пункт с малыми невязками по возможности в середине массива точек. Его координаты в обеих системах ''Xj'', ''Yj'', ''xj'', ''yj'' становятся параметрами ''X''0, ''Y''0, ''x''0, ''y''0.

== Пример вычисления параметров ==

Даны координаты четырёх пунктов ''x'', ''y'' в исходной системе, ''X'', ''Y'' в конечной системе с весами ''p'':

{| class="wikitable"
|-
! ''i'' !! ''xi'' !! ''yi'' !! ''Xi'' !! ''Yi'' !! ''pi''
|- align="right"
| 1 || 1334.71 || 285.94 || 83477.64 || 87377.60 || 1.0
|- align="right"
| 2 || 563.67 || −5197.34 || 82557.14 || 81916.51 || 1.0
|- align="right"
| 3 || 4444.27 || 1153.79 || 86610.19 || 88160.39 || 1.0
|- align="right"
| 4 || −252.07 || 2881.90 || 81962.05 || 90016.34 || 1.0
|- align="right"
| colspan="5" | ∑''p'' = || 4.0
|}

{| class="wikitable"
|-
! ''i''
! ''pixi''
! ''piyi''
! ''piXi''
! ''piYi''
|- align="right"
| 1 || 1334.71 || 285.94 || 83477.64 || 87377.60
|- align="right"
| 2 || 563.67 || −5197.34 || 82557.14 || 81916.51
|- align="right"
| 3 || 4444.27 || 1153.79 || 86610.19 || 88160.39
|- align="right"
| 4 || −252.07 || 2881.90 || 81962.05 || 90016.34
|- align="right"
| ∑ = || 6090.58 || −875.71 || 334607.02 || 347470.84
|}

Взвешенные средние:

{| class="wikitable"
|-
! ''xc'' !! ''yc'' !! ''Xc'' !! ''Yc''
|- align="right"
| 1522.645 || −218.9275 || 83651.755 || 86867.71
|}

Перенос осей в центр масс:

{| class="wikitable"
|-
! ''i'' !! ∆''xi'' !! ∆''yi'' !! ∆''Xi'' !! ∆''Yi''
|- align="right"
| 1 || −187.935 || 504.8675 || −174.115 || 509.89
|- align="right"
| 2 || −958.975 || −4978.4125 || −1094.615 || −4951.20
|- align="right"
| 3 || 2921.625 || 1372.7175 || 2958.435 || 1292.68
|- align="right"
| 4 || −1774.715 || 3100.8275 || −1689.705 || 3148.63
|}

{| class="wikitable"
|-
! ''i''
! ''pi'' ∆''xi'' ∆''Xi''
! ''pi'' ∆''yi'' ∆''Yi''
! ''pi'' ∆''xi'' ∆''Yi''
! ''pi'' ∆''yi'' ∆''Xi''
! ''pi'' (∆''xi''2 + ∆''yi''2)
|- align="right"
| 1 || 32722.3 || 257426.9 || −95826.2 || −87905.0 || 290210.8
|- align="right"
| 2 || 1049708.4 || 24649116.0 || 4748077.0 || 5449445.0 || 25704224.1
|- align="right"
| 3 || 8643437.7 || 1774484.5 || 3776726.2 || 4061095.5 || 10420246.0
|- align="right"
| 4 || 2998744.8 || 9763358.5 || −5587920.9 || −5239483.7 || 12764744.5
|- align="right"
| ∑ = || 12724613.2 || 36444385.8 || 2841056.2 || 4183151.8 || 49179425.3
|-
!   !! ''S''1 !! ''S''2 !! ''S''3 !! ''S''4 !! ''S''5
|}

Параметры конформного преобразования:

{| border="0"
|-
|
{| class="wikitable" style="float:left"
|-
! ''a''1 !! ''b''1
|- align="right"
| 0.99978799 || −0.02728978
|}
{| class="wikitable" style="float:right"
|-
! ''a''0 !! ''b''0
|- align="right"
| 82135.407 || 87128.144
|}
|}

Невязки:

{| class="wikitable"
|-
! ''i'' !! ''vxi'' !! ''vyi''
|- align="right"
| 1 || 0.002 || 0.001
|- align="right"
| 2 || 0.016 || −0.013
|- align="right"
| 3 || −0.032 || −0.016
|- align="right"
| 4 || 0.013 || 0.028
|}

Масштаб и разворот:

{| class="wikitable"
|-
! ''m'' !! ''θ''
|- align="right"
| 1.00016037 || −1°33′48.72″
|}

Сконструируем ключ на основе первого геодезического пункта:

{| class="wikitable"
|-
! ''X''0 !! ''Y''0 !! ''x''0 !! ''y''0 !! ''m'' !! ''θ''
|- align="right"
| 83477.64 || 87377.60 || 1334.71 || 285.94 || 1.00016037 || −1°33′48.72″
|}

== Пример программной реализации ==

Вот листинг простой утилиты на языке C:

<syntaxhighlight lang="c">
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

int main(int argc, char *argv[])
{
char buf[1024], name[32];
double x[2], y[2], wgt;
double xc[2], yc[2];
double dx[2], dy[2], dz[2];
double a[2][2], scale, rotation;
double s[5];
int i;
FILE *fp0, *fp1;

if (argc < 4) {
printf("usage: findkey input-file key-file var-file\n");
exit(EXIT_FAILURE);
}

if ((fp0 = fopen(argv[1], "r")) == NULL) {
printf("can't open %s\n", argv[1]);
exit(EXIT_FAILURE);
}

/* подсчитать сумму координат */
for (i = 0; i < 5; i++)
s[i] = 0.;
while (fgets(buf, 1024, fp0) != NULL) {
sscanf(buf, "%s %lf %lf %lf %lf %lf",
name, &x[0], &x[1], &y[0], &y[1], &wgt);
s[0] += x[0] * wgt;
s[1] += x[1] * wgt;
s[2] += y[0] * wgt;
s[3] += y[1] * wgt;
s[4] += wgt;
}
rewind(fp0);

/* найти центр масс */
for (i = 0; i < 2; i++) {
xc[i] = s[i] / s[4];
yc[i] = s[2 + i] / s[4];
}

/* подсчитать сумму произведений */
for (i = 0; i < 5; i++)
s[i] = 0.;
while (fgets(buf, 1024, fp0) != NULL) {
sscanf(buf, "%s %lf %lf %lf %lf %lf",
name, &x[0], &x[1], &y[0], &y[1], &wgt);
/* вычислить разности */
dx[0] = x[0] - xc[0];
dx[1] = x[1] - xc[1];
dy[0] = y[0] - yc[0];
dy[1] = y[1] - yc[1];
/* суммировать */
s[0] += dx[0] * dy[0] * wgt;
s[1] += dx[1] * dy[1] * wgt;
s[2] += dx[0] * dy[1] * wgt;
s[3] += dx[1] * dy[0] * wgt;
s[4] += (dx[0] * dx[0] + dx[1] * dx[1]) * wgt;
}
rewind(fp0);

/* найти первичные параметры */
a[1][0] = (s[0] + s[1]) / s[4];
a[1][1] = (s[2] - s[3]) / s[4];
a[0][0] = yc[0] - a[1][0] * xc[0] + a[1][1] * xc[1];
a[0][1] = yc[1] - a[1][1] * xc[0] - a[1][0] * xc[1];

/* найти вторичные параметры */
scale = hypot(a[1][0], a[1][1]);
rotation = atan2(a[1][1], a[1][0]);

/* вывести параметры в файл ключа */
if ((fp1 = fopen(argv[2], "w")) == NULL) {
printf("can't create %s\n", argv[2]);
exit(EXIT_FAILURE);
}
fprintf(fp1, "%.3f\n", a[0][0]);
fprintf(fp1, "%.3f\n", a[0][1]);
fprintf(fp1, "%.12f\n", a[1][0]);
fprintf(fp1, "%.12f\n", a[1][1]);
fprintf(fp1, "%.12f\n", scale);
fprintf(fp1, "%+.10f\n", rotation * 180. / M_PI);
fclose(fp1);

/* вывести данные вместе с невязками */
if ((fp1 = fopen(argv[3], "w")) == NULL) {
printf("can't create %s\n", argv[3]);
exit(EXIT_FAILURE);
}
while (fgets(buf, 1024, fp0) != NULL) {
sscanf(buf, "%s %lf %lf %lf %lf %lf",
name, &x[0], &x[1], &y[0], &y[1], &wgt);
/* "наблюдённые" dx[], dy[] */
dx[0] = x[0] - xc[0];
dx[1] = x[1] - xc[1];
dy[0] = y[0] - yc[0];
dy[1] = y[1] - yc[1];
/* "вычисленные" dy[] */
dz[0] = a[1][0] * dx[0] - a[1][1] * dx[1];
dz[1] = a[1][1] * dx[0] + a[1][0] * dx[1];
fprintf(fp1, "%s %.3f %.3f %.3f %.3f %g %.3f %.3f\n",
name, x[0], x[1], y[0], y[1], wgt, dy[0] - dz[0], dy[1] - dz[1]);
}
fclose(fp1);

fclose(fp0);

return 0;
}
</syntaxhighlight>

Сохраним код в файл '''findkey.c'''. Исполняемый модуль можно создать, например, компилятором '''gcc''':

<syntaxhighlight lang="bash">
$ gcc -o findkey findkey.c -lm
</syntaxhighlight>

Для MS Windows можно загрузить уже скомпилированную [http://gis-lab.info/other/findkey.zip программу].

Утилита '''findkey''' запускается в командной строке с именами трёх файлов в качестве аргументов: входной файл данных, выходные файл параметров и файл невязок.

Подготовим файл данных '''data.dat''', структура которого соответствует таблице исходных данных, т. е. содержит в колонках названия пунктов, входные координаты ''x'', ''y'', выходные координаты ''X'', ''Y'', веса:

<pre>
1 1334.71 285.94 83477.64 87377.60 1.0
2 563.67 -5197.34 82557.14 81916.51 1.0
3 4444.27 1153.79 86610.19 88160.39 1.0
4 -252.07 2881.90 81962.05 90016.34 1.0
</pre>

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

<syntaxhighlight lang="bash">
$ findkey data.dat key.dat var.dat
</syntaxhighlight>

получим параметры '''key.dat'''

<pre>
82135.407
87128.144
0.9997879942
-0.0272897781
1.0001603698
-1.56353244
</pre>

и невязки '''var.dat'''

<pre>
1 1334.710 285.940 83477.640 87377.600 1 0.002 0.001
2 563.670 -5197.340 82557.140 81916.510 1 0.016 -0.013
3 4444.270 1153.790 86610.190 88160.390 1 -0.032 -0.016
4 -252.070 2881.900 81962.050 90016.340 1 0.013 0.028
</pre>

== Заключение ==

Положенные в основу статьи формулы называются в учебниках геодезии «неполными» в противоположность «полным». Дело в том, что при большом удалении объекта от осевого меридиана исходной проекции Гаусса-Крюгера возникает значительный градиент масштаба изображения в направлении запад-восток. Чтобы компенсировать возникающие при этом специфические искажения конформного отображения, в «полные» выражения добавляют необходимые члены разложений формул проекций. Разумеется, такой подход несовершенен, как любые ограниченные разложения. В статье [http://gis-lab.info/qa/local-cs.html Добавление местной координатной системы в GIS] в качестве альтернативы предлагается переход от ГСК к проекции с нулевым градиентом масштаба в центре объекта или вблизи него, что делает «полные» формулы ненужными.

== Примечания ==

<references />

== Ссылки ==

* [http://www.google.ru/url?q=http://communities.bentley.com/products/geospatial/desktop/m/geospatial_desktop_gallery/206218/download.aspx&sa=U&ei=YG47UdKNPOeG4gT92oD4Ag&ved=0CDgQFjAJ&usg=AFQjCNFbjYGDGtp3cTXaKerd9LuE8l9_zQ Convert Local Coordinate Systems to Standard Coordinate Systems, McCoy J., 2012]
* [http://user.gs.rmit.edu.au/rod/files/publications/COTRAN_1.pdf Coordinate transformations in surveying and mapping, Deakin R.E., 2004]
* [http://kartoweb.itc.nl/geometrics/coordinate%20transformations/coordtrans.html Coordinate transformations, Knippers R., 2009]
* [http://gis-lab.info/qa/affine-math.html Аффинные преобразования - математика]
* [http://gis-lab.info/qa/local-cs.html Добавление местной координатной системы в GIS]
* [http://gis-lab.info/qa/polynom.html Полиномиальные преобразования]
* [http://gis-lab.info/qa/polynom-calc-examples.html Полиномиальные преобразования - примеры реализации]
* [http://gis-lab.info/qa/rmse.html Среднеквадратичная ошибка (RMSE)]

GIS-Lab - Вклад [ru]

Конформное преобразование