Розділ 8

Геодезія всерйоз

Справжня геодезична мережа порушує два наші зручні припущення: виміри мають різну точність, а зв'язок між вимірами й невідомими буває нелінійним. МНК справляється з обома, додавши ваги і ітерації Гаусса — Ньютона.

8.1. Два нових ускладнення

Проблема

Кілометровий нівелірний хід виміряно гірше, ніж 300-метровий. Чи справедливо ставитися до них однаково?
Віддалемір дає відстань від станції до точки: $d = \sqrt{(x - x_s)^2 + (y - y_s)^2}$. Невідомі $x$, $y$ сидять під коренем. Модель нелінійна, і жоден трюк з розділу 7 її не лінеаризує.

8.2. Ваги: кому довіряти більше

Якщо похибка $i$-го виміру має середню квадратичну величину $\sigma_i$, природно, щоб точніші виміри важили більше. Кожному виміру приписують вагу

$$p_i = \frac{\sigma_0^2}{\sigma_i^2}$$

і мінімізують зважену суму квадратів $\sum p_i v_i^2$. Константа $\sigma_0$ — похибка «еталонного» виміру вагою 1 (середня квадратична похибка одиниці ваги). Нормальні рівняння набувають вигляду

$$A^\mathsf{T}P A\,\hat{\boldsymbol\theta} = A^\mathsf{T}P\,\mathbf l, \qquad P = \operatorname{diag}(p_1, \dots, p_n).$$

Жодного нового методу тут немає: досить помножити кожен рядок $A$ і відповідний елемент $\mathbf l$ на $\sqrt{p_i}$ і розв'язати звичайний МНК. Саме так рахує код цього розділу.

У нівелюванні похибки з кожної станції накопичуються, тому дисперсія перевищення росте пропорційно довжині ходу: $\sigma_i^2 \propto L_i$. Отже, $p_i = 1/L_i$ (довжина в кілометрах). Повернімося до ходу з розділів 1–2.

Розв'язок зваженого МНК для замкненого ходу — поправки пропорційні довжинам: $v_i = -w\,L_i / \sum L_j$. Це і є варіант «пропорційно довжині» з розділу 1, який на практиці вважають правильним. У розділі 2 він програвав за звичайною сумою квадратів, бо тоді ми не враховували, що ходи різні. Кожен розв'язок найкращий за своїм критерієм, і таблиця це показує: без ваг мінімальна $\sum v^2$, з вагами — $\sum p v^2$.

8.3. Зрівнювання нівелірної мережі

Тепер повноцінна мережа: два репери з відомими висотами (Rp1 = 100,000 м, Rp2 = 103,500 м), три невідомі пункти A, B, C і сім виміряних ліній. Невідомих $k = 3$, вимірів $n = 7$, отже надлишкових $r = 4$. Кожна лінія дає рівняння $H_\text{куди} - H_\text{звідки} = h$ з вагою $1/L$.

Додати промах +25 мм у лінію 5 (C → B)

№	Лінія	h, м	L, км	v, мм

Результат зрівнювання — не лише висоти, а й три речі, яких немає в «ручних» способах.

Поправки до кожного виміру $v_i$. Геодезист завжди на них дивиться: поправки мусять бути малими й без системи.
Оцінка точності вимірів $\sigma_0 = \sqrt{\sum p_i v_i^2 / r}$. Ділимо на кількість надлишкових вимірів $r = n - k$, а не на $n$ (розділ 5). Якщо $\sigma_0$ значно більша за очікувану для приладу, у мережі промах.
Оцінка точності результатів: похибка кожної висоти $\sigma_{H_j} = \sigma_0\sqrt{Q_{jj}}$, де $Q = (A^\mathsf{T}PA)^{-1}$. Пункт C, до якого сходиться найбільше ліній, визначено найточніше. Звідки ця формула і що таке еліпс похибок, розповімо в розділі 9.

Увімкніть промах. $\sigma_0$ різко зросте, а найбільша поправка вкаже на лінію 5. Але зверніть увагу: поправки з'явилися і на сусідніх лініях. МНК «розмазує» один промах по всій мережі, намагаючись узгодити все з усім. Саме тому пошук промахів — окрема наука (розділ 10).

8.4. Нелінійні задачі: метод Гаусса — Ньютона

Трилатерація: положення точки визначають за виміряними відстанями до кількох станцій з відомими координатами. Так працює і віддалемірна засічка, і (з поправкою на годинник приймача) GNSS. Рівняння $d_i = \sqrt{(x - x_i)^2 + (y - y_i)^2}$ нелінійні, але ідея розв'язку проста.

Беремо наближення $(x_0, y_0)$. Геодезисти майже завжди знають його з карти чи попередніх вимірів.
Замінюємо кожну відстань її лінійною апроксимацією в околі наближення (дотичною площиною — перший член ряду Тейлора): $$d_i(x_0 + \Delta x,\, y_0 + \Delta y) \approx d_i^0 + \frac{x_0 - x_i}{d_i^0}\,\Delta x + \frac{y_0 - y_i}{d_i^0}\,\Delta y.$$ Коефіцієнти — компоненти одиничного вектора від станції до точки. Разом вони утворюють матрицю $J$ (якобіан).
Розв'язуємо звичайний лінійний МНК $J\,\Delta \approx \mathbf d - \mathbf d^0$ відносно поправок $\Delta$.
Оновлюємо наближення і повторюємо, доки поправки не стануть мізерними.

показати справжню точку

#	x	y	S	\|Δ\|

Що варто помітити:

Збіжність дуже швидка. Кількох ітерацій вистачає навіть з наближення за сотні метрів. Біля розв'язку кількість правильних цифр на кожному кроці приблизно подвоюється.
Кожен крок — звичайний МНК з розділів 4–6, з ваговою матрицею, якщо виміри різноточні. Уся «нелінійність» живе лише в тому, що матрицю $J$ перераховують на кожній ітерації.
Геометрія важлива. Поставте станції в ряд. Кола перетинаються у двох симетричних точках, і метод збігається до тієї, біля якої почав. Чаша з розділу 3 гарантувала одне дно лише для лінійних моделей. У нелінійних задач мінімумів може бути кілька, і наближення вирішує, в який із них ми потрапимо.

8.5. Зважений і нелінійний МНК у машинному навчанні

Зв'язок з ML

Ваги прикладів: параметр sample_weight у fit() більшості моделей sklearn — це та сама матриця $P$. Його використовують, коли дані мають різну надійність або коли класи незбалансовані.
Гетероскедастичність — статистична назва ситуації, коли шум різних спостережень різний. Ігнорувати її — те саме, що зрівнювати ходи різної довжини без ваг.
Нелінійний МНК: scipy.optimize.least_squares і curve_fit реалізують Гаусса — Ньютона з модифікацією Левенберга — Марквардта. Навчання нейромережі з MSE — теж нелінійний МНК, лише з мільйонами параметрів.
Звідки градієнти. Градієнт суми квадратів дорівнює $-2J^\mathsf{T}\mathbf r$. Зворотне поширення помилки (backpropagation) — ефективний спосіб обчислити саме цей добуток, не будуючи $J$ явно. А матрицю $J^\mathsf{T}J$ з Гаусса — Ньютона у глибокому навчанні називають Gauss–Newton matrix. На її наближеннях будують методи другого порядку.
Початкове наближення в нелінійному МНК — те саме, що ініціалізація ваг нейромережі: від неї залежить, у яку «яму» ландшафту втрат потрапить навчання.

🔍 Хочете більше подробиць? Корельовані виміри, Левенберг — Марквардт, дефект датуму, збіжність

Корельовані виміри

Якщо похибки вимірів корельовані (скажімо, кілька кутів з однієї станції), діагональної матриці ваг не досить. Беруть $P = \sigma_0^2\,C^{-1}$, де $C$ — повна коваріаційна матриця вимірів. Це узагальнений МНК (Ейткен, 1935). Геометрично ми просто змінюємо спосіб вимірювати відстань у просторі спостережень: «одиничне коло» стає еліпсом.

Левенберг — Марквардт

Якщо наближення погане, крок Гаусса — Ньютона буває завеликим і «перестрибує» розв'язок. Метод Левенберга — Марквардта розв'язує замість $J^\mathsf{T}J\,\Delta = J^\mathsf{T}\mathbf r$ систему

$$\bigl(J^\mathsf{T}J + \lambda I\bigr)\,\Delta = J^\mathsf{T}\mathbf r.$$

За малого $\lambda$ це Гаусс — Ньютон, за великого — обережний крок уздовж градієнта. $\lambda$ підлаштовують на кожній ітерації. Такий самий доданок $\lambda I$ ми побачимо в ridge-регресії (розділ 10), і це не збіг: в обох випадках він стримує надто великі значення.

Дефект датуму: навіщо відомі пункти

Приберіть з нівелірної мережі обидва репери. Перевищення визначають лише різниці висот, тож усі висоти можна зсунути на ту саму сталу, і жоден вимір цього не помітить. Стовпці $A$ стають лінійно залежними (їхня сума — нуль), і $A^\mathsf{T}PA$ вироджена (розділи 5–6). Геодезисти називають це дефектом датуму. Його усувають, фіксуючи висоту хоча б одного пункту, або шукають розв'язок з найменшою нормою (вільна мережа, розділ 6). Для планових мереж дефект більший: зсув за $x$, за $y$ і поворот, а без виміряних відстаней ще й масштаб.

Коли Гаусс — Ньютон збігається

Метод замінює функцію її лінійною апроксимацією, тож працює добре, коли модель в околі розв'язку «майже лінійна» і нев'язки малі. Для геодезичних задач із наближеними координатами з точністю до кількох відсотків відстані це майже завжди так. Практичне правило зупинки: поправки менші за десяту частку очікуваної точності (скажімо, 0,1 мм). Тоді й оцінку точності рахують з останнього $J$ так само, як для лінійної задачі.

8.6. Код

Python: повне зрівнювання нівелірної мережі з вагами й оцінкою точності та трилатерація Гауссом — Ньютоном з нуля (у коментарі — однорядковий відповідник зі scipy). JavaScript: зважений розподіл нев'язки й ітерації Гаусса — Ньютона з нормальними рівняннями 2×2.

8.7. Перевірте себе

У скільки разів вага однокілометрового нівелірного ходу більша за вагу чотирикілометрового?

Замкнений хід з ходів 1, 1, 2 і 1 км має нев'язку +10 мм. Яка поправка (у мм) дістанеться двокілометровому ходу?

Станція в $(0;\,0)$, наближення точки $(3;\,4)$. Чому дорівнює $\partial d/\partial x$ у цій точці?

Скільки невідомих у задачі визначення положення приймачем GNSS за псевдовідстанями?

Далі: чому саме квадрати

Ми вже рахуємо не лише розв'язок, а й його точність: $\sigma_0$, $\sigma_H$. Звідки ці формули? Чому ваги саме $1/\sigma^2$? І чому взагалі квадрати, а не модулі? Відповідь дає теорія ймовірностей: гауссів розподіл похибок, принцип максимальної правдоподібності й теорема Гаусса — Маркова. Це розділ 9, там же з'являться еліпси похибок.