Расчет случайных ошибок 3D реконструкции на основе аналитической модели Текст научной статьи по специальности «Физика»

Расчет случайных ошибок 3Б реконструкции на основе аналитической модели

Архипов Ю.Б., Иванова Ю.А.(и11а@а81гот1огт.ги) "Астроинформ СПЕ"

Для решения многих прикладных задач дистанционного зондирования немаловажную роль играет определение всех трех координат пространственного положения элементов наблюдаемой сцены по двум и более стереоснимкам. Примерами могут служить оценка рельефа верхней границы облачности тропических циклонов, шлейфов вулканических извержений, индустриальных и природных пожаров. Существуют также примеры из других прикладных областей. В подобных задачах требования к ошибкам знания геометрии съемки определяются конечными требованиями к допустимым ошибкам стереореконструкции. Рассматривается задача о геометрических ошибках 3D реконструкции сцены по двум и более изображениям. Считаем, что имеется 3 группы источников ошибок 3D реконструкции, по-разному влияющие на конечный результат:

1. ошибки положения камеры

2. ошибки ориентации камеры

3. ошибки изображения

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

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

Различные подходы к исследованию ошибок 3D-реконструкции были осуществлены в ряде работ [1, 2, 3, 4]. Первичные ошибки задачи постулировались гауссовыми с нулевым средним. Рассматривалась проективная геометрия построения изображения. В настоящей работе предлагается аналитическое решение задачи, обладающее большей наглядностью. Класс технических задач существенно расширен: геометрия построения изображения не обязана быть проективной (допускается сканирование с произвольным законом по пространству и времени). Для простоты наложено ограничение: распределение ошибок для каждого луча, соединяющего изображение точки и саму точку, должно обладать осевой симметрией, что обычно выполняется с достаточной для практических целей точностью. Не рассматривается корреляция между пространственными координатами близколежащих точек, что повысило бы относительную точность 3D-реконструкции рельефа. Последнее представляет собой отдельную серьезную задачу. Отметим, что существует класс прикладных задач, где такая корреляция отсутствует.

Постановка задачи

Рассмотрим камеру с произвольной внутренней геометрией. Общим случаем является приемная система с произвольным законом сканирования по пространству и времени. Частным случаем является, например, традиционная проективная камера. Точке трехмерного мира соответствует K сопряженных точек на K изображениях и K сопряженных лучей, исходящих из K центров проектирования (эквивалентная оптическая схема). Каждый луч можно характеризовать точкой начала и направляющим вектором.

Ошибки знания геометрии съемки разделяются на ошибки знания положения проективного центра и группу ошибок, приводящих к ошибке знания ориентации сопряженных лучей. Эта последняя группа ошибок состоит из:

- ошибок знания ориентации аппарата;

- ошибок союстировки осей аппарата и камеры;

- остаточных ошибок знания внутренней геометрии камеры;

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

В случае возможности подгонки параметров геометрии съемки по опорным точкам ошибки знания внешней геометрии заменяются на неустранимые ошибки подгонки. Фактически для целей исследования точности ЗБ-реконструкции группу ошибок, имеющих характер угловых ошибок, можно объединить в одну (двухкомпонентную) ошибку знания углового положения сопряженных лучей. В качестве исходных примем следующие допущения:

- ошибки знания точки начала луча по трем координатам имеют гауссово распределение с нулевым средним и одинаковыми дисперсиями ;

- ошибки знания направления луча будем характеризовать гауссовым законом распределения с нулевым средним и одинаковыми дисперсиями по координатам о2а точки

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

Тогда ошибки по координатам точки пересечения луча с перпендикулярной плоскостью, расположенной на расстоянии Ь от проективного центра, также будет распределена по нормальному закону с дисперсией

< = < + Ь2с2л , к = 1,К К, (1)

где К-количество изображений стереосъемки. При написании формулы (1) изначально

предполагалось, что диаграмма углового рассеяния, определяемая параметром аак,

достаточно узкая, что дает право пренебрегать изменением члена Ьк <ал на расстояниях порядка а кк.

Плотность распределения точки пересечения расчетного луча с плоскостью, перпендикулярной истинному лучу и находящейся на расстоянии Ьк от проективного центра записывается в виде

Рк( = 0, yt = 0) = -V ехр{- ^^ (2)

2па Ьк I 2а Ьк

где х^ yt - координаты точки пересечения плоскости истинным лучом.

На рисунке 1 для случая двух лучей штрихпунктирной линией изображены истинные лучи, точка пересечения которых является истинной точкой, прообразом сопряженных точек на снимках. Степень рассеяния расчетных лучей относительно истинных схематически изображена фигурами вращения, радиус сечения которых является среднеквадратическим отклонением точки пересечения лучом плоскости, расположенной на расстоянии Ьк от проективного центра.

Рис. 1

Основываясь на описанной выше модели явления выберем алгоритм вычисления (оценки) 3Б точки и найдем статистику ошибок применительно к этой оценке, исходя из того, что истинные сопряженные лучи неизвестны.

Нахождение апостериорной плотности вероятности 3D точки

Рассмотрим обратную задачу. Пусть теперь на рисунке 1 штрихпунктирными линиями изображены расчетные лучи, положение которых известно. Из-за присущих любой технической системе ошибок в измеренных, а также вычисленных по некоторым моделям величинах, положение расчетных лучей известно неточно. В частности, в трехмерном пространстве расчетные лучи, вообще говоря, не пересекаются.

Используя предыдущие определения, найдем статистику рассеяния истинных лучей

относительно расчетных. Проведем плоскость P, перпендикулярную (для определенности)

истинному лучу. Запишем для векторных координат t, m точек пересечения луча с

плоскостью P соотношение Байэса

p(m 11) p(t) p(t | m) = ——— — PV ' ' Jp(m|t) p(t) d(t)

или

p(m 11) p(t)

p(t Im) = ' \ (3)

p(m)

где p(t | m) - плотность вероятности истинного события t при условии наблюдаемого m. Ввиду малости поля ошибок можно считать априорную вероятность p(t) пренебрежимо мало изменяющейся в пределах этого поля. Аналогично, в пределах поля ошибок также мало изменяется априорная плотность события p(m), причем с достаточной точностью p(t) = p(m). Отсюда, учитывая (3), вытекает

p(t | m) = p(m | t) (4)

Перейдя в (4) к плотностям ошибок m - t p(t | m) = p(m - t | m), p(m | t) = p(m - t | t) придем к

p(m - t | m) = p(m - t | t) (5)

Последнее означает, что при узких полях ошибок условные распределения полей ошибок (5) практически совпадают.

Окончательно, учитывая центральную симметрию распределения (2) и соотношение (5), имеем для апостериорной вероятности распределения истинного луча

1ПЬк V lGLk J

где Xmk, ymk - координаты измеренного луча.

Далее будем искать оценку R для трехмерной истинной точки Ru и плотность распределения трехмерной ошибки р(е), где

е = R - Ru. (7)

Интересующая нас точка Ru реально может лежать только в узкой зоне сближения лучей. Максимальный линейный размер этой зоны при не слишком малых углах стереобазы имеет порядок максимальной из величин aLk. Учитывая малость входных ошибок, в зоне сближения лучей пренебрежем в формуле (1) изменением aLk от расстояния Lk. Перепишем формулу (1) в виде

= < + CWOk , к = 1,...K (8)

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

Переходя теперь к трехмерной апостериорной плотности распределения истинной точки, с учетом (6) и (8), имеем в зоне сближения лучей цилиндрическое распределение плотности

Рк(х,, | ^,^) = 4-^ехр(- (Х-^)2+ f"-^)21

или

Рк(R|mк) = Gk expj- 2^-J (9)

где

Ak, Gk - неизвестные постоянные коэффициенты (в пределах малой зоны сближения лучей); R - трехмерная координата точки пространства; mk - связанный вектор k-го измеренного луча; rk - расстояние от точки R до k-го луча; aLk - определяется формулой (8).

В силу независимости измерений плотность апостериорной вероятности равна произведению плотностей апостериорных вероятностей

p(R I m„m2,...m^) = Пpk(RI mk) (10)

к=1

Оценка (вычисление) 3D точки по принципу максимального правдоподобия

В соответствии с принципом максимального правдоподобия в качестве оценки истинной точки примем точку R, доставляющую максимум выражению (10). Для дальнейших оценок ограничимся случаем двух лучей. В этом случае формула (10) принимает вид

p(R I m1( m 2) = pj(R | m,) • pj(R | m 2) = Gfi.expl-^| (11)

На рис. 2 представлена геометрия съемки одной физической точки Rи в проекции на плоскость, параллельную расчетным лучам 1 и 2, соответствующим этой точке. Из-за присущих процессу ошибок лучи, вообще

говоря, не проходят через точку Rи и не

пересекаются между собой.

В зоне максимального сближения лучей в соответствии с предположением о малости ошибок изображены цилиндрические участки плотности вероятности нахождения истинной точки Rи для каждого луча.

Очевидно, максимум выражения (11) достигается на минимальном отрезке, соединяющем лучи. В самом деле, рассмотрим произвольную точку R, не лежащую на минимальном отрезке. Опустим из точки R перпендикуляр на минимальный отрезок. Rg -точка основания перпендикуляра. Очевидно

Р1(К | Ш1) < pl(Rg | Ш1>, Р2(Я | Ш1) < p2(Rg | Ш1)

Таким образом, точка максимума вероятности лежит на минимальном отрезке.

1

Рис.2

Итак, будем искать максимум апостериорной вероятности на минимальном отрезке. Для этого перепишем формулу (11) в виде:

>(г|ш 1 ,ш2)= Р1 (г|ш 1 )р2(гШ2)= О ■ exp

й

г + — 2

й

г--

2

2а:

2а2

(12)

где

О - постоянный коэффициент;

й

г - координата вдоль минимального отрезка, соединяющего лучи, причем г = - — на первом й

луче и г = — на втором луче; d -.расстояние между лучами.

Здесь и далее под о1 и о2 подразумеваются оы и оь2.

В соответствии с принципом максимума правдоподобия находим точку € максимума (12)

€=

1 +

а2

Ча1 У

-1

2

(13)

й

й

Как и следовало ожидать, при изменении о2/о1 от да до 0 € изменяется от--до —, т.е.

2 2

всегда находится между экспериментально измеренными лучами; при а2 =а1 точка максимума лежит посредине между лучами.

Плотность распределения вектора ошибок ЭБ реконструкции

Рассмотрим геометрию зоны пересечения лучей визирования в проекции на плоскость, содержащую истинные лучи визирования (рис. 3). Введем три естественные системы

2

2

2

координат: ХУЪ, Х1У121, Х2У2Ъ2, в которых описание явления имеет наиболее простой вид. Центры этих систем координат находятся в истинной точке Ми. Оси и Ё2 направлены на центры проектирования. Ось Z делит угол 0 между осями и Ё2 пополам. Оси X, Х1; Х2 лежат в плоскости истинных лучей.

ъ

Рис. 3

Рассмотрим расчетные лучи, характеризуемые векторами ш1 и ш2 (на том же рисунке). На минимальном отрезке, соединяющем эти лучи, находится оценка М истинной точки Ми.

В рамках предыдущих допущений расчетные лучи параллельны истинным и

характеризуются своим положением Хк, Ук и плотностями вероятностей

р.(|.. (14)

Задача состоит в нахождении плотности вероятности р(М), равной плотности вероятности вектора ошибки е = М - М и.

Обозначим Х1, У1, Х2, У2, реализации компонент первого и второго расчетных лучей в своих системах координат и выразим через них координаты (в системе ХУЪ) Х-, У, Ё расчетной

точки. Формулы для Х и Ё получаются из чисто геометрических соображений; для У - с учетом принятой оценки (5).

X = Х1 + Х 2

0 0 2cos— 2

а 2 а 2

у€=У1 +у (15>

а1 +а2 а1 +а2

€=Х1 - Х ^

0

2sin — 2

Очевидно, при принятых выше допущениях, математические ожидания компонент вектора Й = ((, У, *) равны нулю.

Найдем дисперсии и ковариации компонент вектора

Г \

ковариации: соу

используя свойство билинейности

^,^ь1п1 aibj соу(, п)• Имеем, с учетом независимости

У i, 1

смещений Х1, У1, Х2, У2

yar

_2

X = ст_ + а2

, 2 0 ' 4cos2 —

yar У=

соу1

((, Iе) = 0, соу((, *)--

22

22 -а2

2sin0

yar

22 е ст1 + ст2

*=7Ю

4sin2 —

(16)

соу1

= 0.

Распределение Й нормальное, как распределение линейных функций нормально распределенных независимых величин Х1, У1, Х2, У2, а именно, Й имеет 3-х мерное нормальное распределение с нулевым средним и ковариационной матрицей А = (ау)

А =

(„2.1 ст1 + ст2

4^2

22 ст1 -ст2

2sin0

2

22 ст,2 + ст:

2sin0

22 ст2 + ст22

. . 2 0 4sln2 — 2

(17)

Таким образом, получена ковариационная матрица ошибок (17) в системе координат X, У, Ъ (рис. 3), которую в дальнейшем будем называть специальной системой координат. Как видно из (17), в специальной системе координат отсутствует корреляция ошибок по У с ошибками по другим координатам. Назовем собственной системой координат такую систему, в которой корреляционная матрица ошибок имеет диагональный вид. Как и следовало ожидать, при ст1 = ст2 специальная система координат является собственной (см. выражение (17). Приведем теперь известные формулы, позволяющие пересчитывать корреляционную матрицу ошибок (17), полученную в специальной системе координат, к произвольной (декартовой) системе координат, которую далее будем для определенности называть местной.

Обозначим компоненты случайного вектора У в местной системе координат V и в специальной системе координат V. Пусть матрица поворота Я между этими системами координат определена формулой

V = ^ (18)

Выразим ковариационную матрицу вектора V

В = соу(У, V) (19)

через ковариационную матрицу А (17). Обозначая знаком Е математическое ожидание и используя линейность матричных операций, получим

В = Е((У - ЕУ)(У - EV)T) = Е((^ - ЕКу) - ЕКу)1) = - Е^ (V - Ev)TRT)

= - Е^ (V - Ev)T)RT = Я соу(^)ят = ЯАЯт

Окончательно

В = ЯАЯ (20)

Представление о величинах ошибок стереореконструкции в местной системе координат можно получить, вычислив СКЗ ошибок по осям местной системы координат

ст =

(21)

2

2

2

0

0

2

22 стст

1 ~ 2

0

0

0

Результаты и обсуждение

На рисунках 4-8 (см. приложение) представлены примеры расчетов ошибок стереореконструкции сцены, расположенной на характерной высоте 10 км (облачность) при съемке с двух точек, расположенных на высоте 500 км.

СКО ошибок положения точек съемки принято равным 12 м. В качестве СКО угловых ошибок принята достаточно произвольная величина 100 мкрад. При данном соотношении ошибки положения, угловой ошибки и расстояний съемки доминирующей (см. формулу (1)) является угловая ошибка.

Обсуждение зависимости угловой ошибки от многочисленных факторов, упомянутых в разделе "постановка задачи", выходит за рамки данной работы. Влияние кривизны Земли в данных результатах не отражено (хотя соответствующие расчеты также были проведены). Представление об ошибках в собственной системе координат дают значения ошибок по осям Х,У, Ъ в точке Х=У=0. На значения ошибок в остальных точках рассматриваемой плоскости оказывают влияние такие факторы, как изменение расстояний до точки визирования, изменение соотношения расстояний от точки визирования до спутников, изменение угла между двумя линиями визирования, а также фактор перераспределения ошибок между собственной и местной системами координат (см. формулу (20)). Последний фактор из перечисленных наглядно проявляет себя при рассмотрении зависимостей ошибок при малых углах стереобазы (рисунки 1а-3 а ), когда доминирующей ошибкой в специальной системе координат является Ъ-компонента ошибки. При этом основной вклад в Х и У ошибки в местной системе координат вносит ъ-ошибка специальной системы координат. Рисунок 86 иллюстрирует зависимость Ъ-ошибки от стереобазы для некоторых характерных точек. При малых стереобазах хорошо виден "эффект ножниц". Напомним, что в данной работе не исследовалось формирование входных ошибок. Учет увеличения угловой ошибки нахождения сопряженных точек при увеличении стереобазы привело к появлению на графиках рисунка 8б соответствующих минимумов.

Заключение

Разработан простой и физически прозрачный метод оценки ошибок стереореконструкции трехмерной сцены по двум ее изображениям. Метод нечувствителен к способу получения изображений (возможна любая комбинация изображений типа проективного, сканирующей линейки, развертки одного элемента изображения и т.д.). В качестве критерия выбора оценки 3Б-точки принят критерий максимума правдоподобия. Результаты расчетов иллюстрируют влияние геометрических параметров задачи на величину и распределение компонентов ошибки 3Б-реконструкции, в том числе на перераспределение компонентов ошибки между специальной и произвольной системами координат.

Дальнейшее развитие исследований предполагает как сравнение аналитических результатов данной работы с результатами статистического моделирования, так и количественное исследование влияния различных факторов на формирование входных ошибок. Особый интерес представляет синтез оптимального алгоритма 3Б-реконструкции на основе статистической модели рельефа, разработка которой представляет собой самостоятельную серьезную задачу.

Список литературы

1. Репин В.Г., Тартаковский Г.П. «Статистический синтез при априорной неопределенности и адаптация информационных систем», М., «Советское радио», 1977 г.

2. Бакут П.А., Шульц С.В. «Оценивание трехмерных координат точечного объекта по космическим стереонаблюдениям». Космические исследования , 1998, т.36, №4, с.407-416.

3. Jia Zong, Roger Davies, Jan-Peter Muller, David S. Diner. "Photogrammetric Retrieval of Cloud Advection and Top Height from the Multi-Angle Imaging Spectroradiometer (MISR)". Photogrammetric Ingineering & Remote Sensing, 2002, 68(8):821-829.

4. C. Vicent Tao, Yong Hu. "3D Reconstraction Methods Based on the Rational Function Model". Photogrammetric Ingineering & Remote Sensing, 2002, 68(7):705-714.

Приложение

Расчетные данные: высота спутников 500 км, высота объекта наблюдения 10 км, СКО положения спутников 12 м, СКО знания углового положения 100 мкрад.

а)

а)

500"

-500-

б)

2000 1000 500 1500 1000 2000

1 \ 1 Ист / I

Рис.4 Среднеквадратическая ошибка стереореконструкции по координате Х, м а) для стереобазы 6° (50 км), б) для стереобазы 60° (566 км)

500'

-50&

б)

Рис.5 Среднеквадратическая ошибка стереореконструкции по координате У, м а) для стереобазы 6° (50 км), б) для стереобазы 60° (566 км)

а) V,

500-

-500-

б)

Рис.6 Среднеквадратическая ошибка стереореконструкции по координате Ъ, м

а) для стереобазы 6° (50 км), б) для стереобазы 60° (566 км)

а) б)

Рис.7 Среднеквадратические ошибки стереореконструкции для стереобазы 60° (566 км) по координатам а) X, б) У, в) Ъ

а)

б)

сг, И

0 20 40 60 80 100

стереобата, градусы

Рис.8 Зависимость среднеквадратической ошибки стереореконструкции по координате Ъ от угловой стереобазы для трех характерных точек. а) выбранные точки, б)график зависимости