Продолжение
вот этого.
Наконец-то доказал свою формулу для оценки шума. Разумеется, всё это уже было выведено в мат. статистике, в регрессионном анализе, но у нас его не преподавали, была теория вероятностей, линейная алгебра, и даже на старших курсах - теория случайных процессов, и как-то так вышло, что эту тему я аккуратно в своё время обошёл по кругу со всех сторон, да так и не изучил. Потом несколько раз порывался, но тяжело это, очень специфическая штука, и очень всё абстрактно.
Впервые НА ПРАКТИКЕ мне понадобилась идемпотентная матрица, за это надо будет выпить...
Под катом - моё доказательство, применительно конкретно к видеоизмерителю параметров сближения, с примерами матриц из прошлого поста.
Так уж исторически у меня повелось, что я каждую точку (каждый уголковый отражатель) рассматривал по отдельности. Говорил: вот координаты точки на фотоприёмной матрице, (fny, fnz), а вот такие мы ожидаем: (gny, gnz), а вот матрица M размером 6×2, по сути, набор из 6 двумерных векторов, показывающих, куда поползёт точка на фотоприёмной матрице, если начать менять один из 6 параметров сближения.
Но на деле вот эта размерность 2 (поскольку точки регистрируются на фотоприёмной матрице) оказывается не особенно полезной. Можно от этого абстрагироваться, не пытаться группировать все входные данные по точкам. Просто мы заявляем: у нас имеется N "каналов измерения", которые каким-то известным образом реагируют на изменение D параметров сближения. И не столь важно, что каналы 1, 3, 5, и так далее - это y-координата точек, а каналы 2, 4, 6, и так далее - это z-координаты. (у нас было определено, что ось X идёт вдоль оптической оси объектива, поэтому фотоприёмная матрица занимает плоскость YZ). В случае "стереорежима", например, у нас на каждую точку может получиться по 4 "измерительных канала" - это координаты точки в "левом глазу" и в "правом глазу", и далее решать будем в точности так же!
Сейчас мы, по крайней мере, для дальнейших выкладок, "освобождаем" одну размерность. Теперь у нас будет вектор невязок H размером N×1, где N - количество "измерительных каналов". В "игрушечном видеоизмерителе" из прошлого поста было N=3, поскольку у нас было 3 точки, но мы работали с фотоприёмной ЛИНЕЙКОЙ, т.е для каждой точки измеряли лишь одну координату. В настоящем видеоизмерителе, работающем по 8 уголковым отражателям, будет N=16, поскольку для каждой точки измеряется две координаты. Невязка - это разность между измеренными координатами и "модельными", которые мы ожидали получить при текущих параметрах сближения.
В "игрушечном видеоизмерителе", вектор невязок в начале работы равнялся:
Далее, у нас будет одна-единственная матрица частных производных, M, размером D×N, где D - количество параметров сближения, которые мы измеряем. В "игрушечном видеоизмерителе" мы рассматривали случаи D=0 (здесь его рассматривать не будем!), D=1, D=2 и затем, мельком, D=3, при котором задача перестаёт быть оптимизационной и даёт нулевую неустранимую невязку. Наконец, в реальном видеоизмерителе, D=6.
Ранее у нас было матриц Mn по количеству точек, но эти "старые" матрицы мы "утрамбовали" в столбцы, и из этих столбцов собрали одну большую матрицу. Каждый её столбец выражает, как соответствующий канал измерения отреагирует на изменения каждого из параметров сближения.
При измерении 1 параметра сближения в "игрушечном видеоизмерителе", мы получали матрицу:
При измерении 2 параметров сближения мы получали матрицу:
А для 3 параметров сближения мы матрицу не строили, но она имеет размеры 3×3.
Далее, у нас появляется ненормированная оценка сдвига:
Поскольку матрица M имеет размер D×N, а невязка H - размеры N×1, то ненормированная оценка сдвига K имеет размеры D×1, т.е по количеству параметров сближения. К каждому из них мы и посчитали "поправку", которая уже должна вести в правильную сторону (если перекрёстные связи не чересчур велики), но масштаб может гулять в огромных пределах...
Потом мы вводим нормирующую матрицу J:
её размеры: D×D, она симметричная и, как минимум, неотрицательно определённая, а если задача имеет решение - положительно определённая. Контрпример, когда матрица так и останется неотрицательно определённой, с нулевым определителем: попытка измерить "параметры сближения" до звёзд. Реакция на "приближение к звёздам" на 1 миллиметр будет нулевой, насколько будет хватать точности, т.е параметр дальности у нас "повиснет в воздухе". Более реалистичный пример: "плоская" мишень, например, мы нарисовали на стене 5 точек, и поставили видеоизмеритель строго перпендикулярно на приличчном расстоянии. Тогда реакция на "поворот стены" будет в первом приближении нулевой. А вот для "объёмных" мишеней, где элементы не лежат в одной плоскости, такой проблемы нет.
И наконец, мы получаем оценку вектора параметров:
Как видно, исходный вектор невязок мы домножаем на хитрючую матрицу, построенную из матрицы M. Подобное выражение мы уже встречали при
выводе аффинного алгоритма, называется это псевдообратной матрицей Мура-Пенроуза, выползает в регрессионном анализе буквально повсюду.
Но на этом мы не останавливаемся! Теперь нам хочется посчитать новые невязки, после того, как мы подправили исходные параметры сближения на вектор s. В реальной задаче мы это делаем по исходным формулам, ДО ЛИНЕАРИЗАЦИИ. Например, у нас проекция точки на плоскость записывалась как fy=y/x. Вот уточнив значения x,y конкретной точки в системе координат видеоизмерителя, мы честно выполняем это деление.
Для оценки шумов, однако, мы можем остаться в рамках линеаризированной задачи. В таком случае мы можем понять, насколько уменьшатся невязки, с помощью матрицы M. Первая её строка - это куда уползут ожидаемые значения (координаты точек), если первый параметр сближения увеличить на единицу. Вторая строка - если второй параметр сближения увеличить на единицу, и так далее. И тогда нетрудно заметить, что суммарное изменение координат точек после поправки параметров сближения на вектор s - это:
Эту поправку надо ВЫЧЕСТЬ из исходного вектора невязок, поскольку невязка - это разница между измеренным и МОДЕЛЬНЫМ. Это у нас как раз поправка к МОДЕЛЬНОМУ значению, поэтому знак минус. Новый вектор невязок:
E - это единичная матрица. Получилось выражение ещё интереснее, зависящее исключительно от матрицы M. Размер матрицы, которая в скобках: N×N.
Назовём эту матрицу как-нибудь, например, Q:
Попробуем возвести её в квадрат:
Как видим, матрица осталась той же самой! Такие матрицы называют идемпотентными.
В общем-то, мы должны были это ожидать! Это означает, что если мы попытаемся провести вторую итерацию данного алгоритма, она уже ничего не изменит! В линейной задаче итерации не нужны, всё, что можно подправить, подправляется в один заход. Это когда исходная задача нелинейная, у нас по мере приближения к ответу будет корректироваться матрица M, но как только мы достаточно близко подошли к ответу, нелинейностями уже можно пренебречь, и тогда ответ перестаёт меняться, сколько итераций мы не возьмём.
Чтобы не быть голословными, приведём примеры этих матриц для "игрушечного видеоизмерителя". По сути, мы "ручками" посчитали их, только не изобразили в явном виде. Первый раз мы записывали неустранимые невязки при измерении ОДНОГО параметра сближения:
То же самое можно записать так:
а вектор из 3 случайных величин - это наши исходные невязки, H.
Ещё, кстати, матрица Q является симметричной, это легко проверяется при попытке оттранспонировать выражение для неё.
А вот такая матрица получается после измерения двух параметров сближения:
В прошлый раз я умудрился не заметить, что все 3 строчки отличаются только множителем! Первая с третьей вообще идентичны, а вторая получается из них умножением на (-2). Ранг матрицы: 1. А у предыдущей матрицы ранг 2, там сложение первых двух строк даст третью со знаком минус.
Идемпотентность можно проверить "прямой подстановкой", это забавно :) На удивление приятный получился пример, не зря за "игрушечный видеоизмеритель" взялся. Всё-таки надо эти матрицы ручками "потрогать", чтобы хорошо всё понять, а ковырять матрицы 6×6, или, того хуже, 16×16 (характерные размеры для "взрослого" видеоизмерителя) - удовольствие ниже плинтуса...
Мы близки к цели. Осталось к этому делу привязать нашу сумму квадратов неустранимых невязок...
По сути, нам надо возвести в квадрат все компоненты матрицы - и сложить их! Так выходит, поскольку когда мы возводим в квадрат линейную комбинацию от независимых случайных величин с последующим взятием мат. ожидания, зануляется всё, кроме квадратов, т.е когда случайная величина умножается на саму себя.
Когда мы
выводили аффинный алгоритм, мы сталкивались с похожей задачей, и поняли, как "в матричных понятиях" записать эту сумму квадратов компонентов матрицы. Для этого умножаем матрицу на транспонированную - и берём след от получившейся матрицы!
Покажем, как это работает для произвольной матрицы 3×3:
Поскольку наша матрица Q симметричная, её умножение на транспонированную - то же самое, что возведение в квадрат. А поскольку она идемпотентная, то и возводить в квадрат её не требуется, она останется той же самой!
Остаётся лишь взять след! Поехали:
Погодите, что-то не то... Вот же у нас выше матрицы перед глазами, где след отнюдь не нулевой... Да и нулевой след означал бы нулевую сумму квадратов невязок, а это очевидно не так... Всё дело в том, что единичные матрицы E бывают разные! Изначальная матрица, след от которой мы берём, имеет размер N×N, поэтому след единичной матрицы в левой части выражения: N.
Но когда мы решили переставить местами матрицы в правой части (это разрешённая операция при взятии следа), мы от матрицы N×N перешли к матрице D×D! Поэтому и единичная матрица там имеет след D, а вовсе не N. Перепишем:
Иными словами, если каждый компонент исходной невязки имел дисперсию σ2, то сумма квадратов невязок после коррекции параметров сближения будет равна (N-D)σ2.
Так что моя формула
была верна! (правда, сейчас я под N понимаю количество каналов измерения, а тогда - количество точек, из-за чего в той формуле было 2N. И вместо произвольного количества параметров сближения D я тогда брал конкретное число 6)
В том, что след матрицы Q равен N-D, можно убедиться на матрицах из нашего "игрушечного видеоизмерителя". В первой из них D=1 (измеряем лишь один параметр), и след равен двум. Во второй, D=2, и след равен единице. Ещё можно убедиться, что след такой идемпотентной матрицы равен её рангу.
И наконец, заметим ещё одну вещь: если N=D, то матрица Q обратится в ноль! Доказывается это в одну строчку. Вот это выражение:
в общем случае упростить нельзя. Мы не можем "разбить" центральную скобочку, отдельно инвертируя M и MT, по той причине, что матрица M обычно ПРЯМОУГОЛЬНАЯ, размера D×N, обратную можно брать только к квадратной матрице! Но если N=D, значит, нам повезло:
Выражение свелось к единичной матрице. Соответственно, Q обращается в ноль.
Это и понятно: у нас хватает степеней свободы, чтобы полностью устранить любые невязки. Или, скорее, нам совсем "впритык" хватает данных, чтобы посчитать параметры сближения, но уже совсем не остаётся избыточности, чтобы убедиться в согласованности этих данных, оценить шумы измерения. В принципе, уже исходя из следа матрицы, равного N-D, можно было заключить, что матрица Q становится нулевой. Ведь N-D - это сумма квадратов всех её элементов, поэтому в ноль это выражение может обратиться лишь для нулевой матрицы. Или вспомнить, что ранг идемпотентной матрицы равен её следу, так что нулевой след - это и нулевой ранг, а это и есть нулевая матрица.
Сейчас мне вся эта математика понадобилась, чтобы оценить точность "контрольных" измерений с помощью теодолита. Причём не просто "для себя", а чтобы и метрологи не придрались, и приглашённые ими эксперты откуда-нибудь из ВНИИФТРИ, которые будут аттестовывать эти все методики... Потому как мало сделать прибор - ты ещё попробуй докажи, что он работает как положено!
Безусловно, велосипед изобретаю. По крайней мере, когда мы обсуждали эту проблему с "Промышленной геодезией", они говорили, что это для них стандартная задача, со стандартными решениями. У лазерного трекера чуть ли не соответствующий "режим измерений" есть, говоришь ему - вот есть такой объект, есть эдакий, измерь их и скажи, как они взаимно расположены. А у нас тут, ВНЕЗАПНО, завалялся тахеометр, надо будет ещё в нём покопаться как следует, вдруг он тоже шибко умный, всё за нас сделает? Боюсь, впрочем, что нет: координаты всех уголковых отражателей в пространстве он, конечно, выдаст, а вот превращать это в параметры сближения всё равно придётся "ручками", в смысле, какими-то своими вычислительными методами, а их, хочешь не хочешь, надо будет аттестовать. А если всё равно такая волокита - тахеометр мне задачу ничуть не облегчит. Пожалуй, точность повысит, но мне и точности теодолита за глаза. Впрочем, тахеометр "потискать" всяко неплохо бы...