CC BY
46
НАУЧНЫЙ ВЕСТНИК МГТУ ГА серия Радиофизика и радиотехника
УДК 629 7 778
МАТЕМАТИЧЕСКАЯ МОДЕЛЬ КОСВЕННОЙ ОПТИМАЛЬНОЙ ПРИВЯЗКИ ИЗОБРАЖЕНИЙ БОРТОВОГО МНОГОКАДРОВОГО ДАТЧИКА К КАРТЕ МЕСТНОСТИ В ЗАДАЧАХ ДИСТАНЦИОННОГО ЗОНДИРОВАНИЯ ЗЕМЛИ
В.Н. СТЕПАШКИН, А.Е. ГРИЦЕНКО
Статья представлена доктором технических наук, профессором Логвиным А.И.
Представлены математические соотношения для расчета поправок к значениям параметров фотограмметрической модели многокадрового бортового датчика изображения, обеспечивающих оптимальную координатную привязку кадров без использования опорных ориентиров местности.
Введение
Предварительная обработка видеоинформации, получаемой с помощью бортовых датчиков изображения (ДИ), включает геометрическую коррекцию - трансформирование (ортотрансформирование) снимков, а также привязку изображений к карте местности. От успешного решения данной задачи зависит эффективность последующих процедур:
- определения координат изображений объектов в системе координат (СК) местности;
- сшивки кадров видеоряда вдоль маршрута полета с построением фотоплана (ортофотоплана) местности;
- получения стереоскопического изображения местности и стереовидео;
- комплексного дешифрирования изображений, получаемых различными датчиками одновременно в нескольких диапазонах спектра электромагнитного излучения и др.
Применительно к случаю малоразмерного беспилотного летательного аппарата (МБЛА) информация, получаемая бортовыми ДИ и датчиками пилотажно-навигационной системы (ПНС), имеет ряд особенностей:
- малая площадь захвата местности отельными кадрами;
- относительно большое число кадров на единицу длины маршрута при съемке как видеокамерой или тепловизионной системой, так и цифровым фотоаппаратом, а также значительная степень перекрытия кадров;
- невысокая точность измерения углового положения МБЛА бортовой ПНС;
- высокая динамика пилотажно-навигационных параметров (ПНП) МБЛА;
- сложность синхронизации и установления соответствия между моментами измерения пи-лотажно - навигационных параметров (ПНП) и моментами получения кадров изображения.
Для решения данных проблем и обеспечения возможности автоматического выполнения перечисленных выше задач предлагается на этапе предварительной обработки введение процедуры сквозной взаимной привязки изображений кадров видеоряда на основе комплексирования пилотажно-навигационной информации (ПНИ) с данными о рассогласовании положения, о параллаксах, отдельных точек (фрагментов) на смежных кадрах видеоряда, получаемых независимо от измерений ПНП.
1. Постановка задачи
Обозначения: Хла,Ула - прямоугольные координаты летательного аппарата в фотограмметрической СК (ФСК); 2ла - абсолютная высота полета; Нла - текущая относительная высота по-
лета ЛА над надирной точкой местности (Хн, Ун); Ж и W - вектор и модуль вектора путевой скорости ЛА; Wx, Wy, Wz - проекции вектора путевой скорости на оси ФСК; Хп,Уп,2п,Ип^х,п^у;пЖ;п - параметры Xла,Yла,Zла,Hла,Wx,Wy,Wz в момент формирования п-того кадра изображения; X, Y, Z=Z(X,Y) - прямоугольные координаты некоторой точки местности в ФСК, имеющей координаты (1,]) в СК изображения; Z=Z(X,Y) - цифровая модель рельефа (ЦМР); H(X,Y)=ZЛa-Z(X,Y) - высота полета по отношению к точке ^^^) местности; ила, ула, ула, Сла - продольный (тангаж), поперечный (крен) и курсовой углы, измеряемые датчиком углового положения и матрица направляющих косинусов данных углов; иу, уу, уу, Су — продольный, поперечный и курсовой углы установки датчика изображения относительно осей чувствительности датчика углового положения и матрица направляющих косинусов данных углов; и, у, у, С=||ер,ч||, р^=1,2,3 - углы ориентации ДИ относительно осей ФСК и матрица направляющих косинусов данных углов; I, I - размеры кадров; положение оптического центра (¡о, ]0) и эквивалентное фокусное расстояние; £,, - параметры внутреннего ориентирования снимков.
Постановка задачи на описательном уровне
Рассматривается участок маршрута полета МБЛА с установленным на его борту кадровым датчиком изображения, работающим в режиме с перекрытием кадров не менее 60%. Производится регистрация изображений, а также ПНП, причем последние регистрируются с погрешностями, обусловленными ограниченной точностью измерителя ПНП, рассогласованием моментов времени регистрации ПНП и кадров изображения.
В общем случае в наличие имеется ЦМР.
Требуется привязать каждый из снимков видеоряда к карте местности путем определения параметров функции пересчета координат точек изображения в координаты их прообразов на местности. При этом параметры определяются так, чтобы обеспечивалось в среднем:
- минимальное рассогласование пересчитываемых к местности координат соответственных точек по всем парам смежных кадров видеоряда;
- минимальное отклонение прообразов (Xn,н,Yn,н) надирных точек (¡н(п),]н(п)) снимков от линии пути ЛА.
Математическая постановка задачи
Исходные данные:
- последовательные кадры изображения, полученные с перекрытием в моменты времени 1п: {х(п)1,]}, 1=0, ..,1-1; ]=0,..,1-1; п=0,..,К;
- результаты измерений вектора ПНП, относимые к моментам ^ и выполненных с некоторой неизвестной погрешностью: (Xn=XЛa(tn), Yn=YЛa(tn), Нп=НЛа(1п), Wx,n=Wx(tn), Wy,n =Wy(tn),
Wz,n =Wz(tn), и(п)=ила(0, у (п)=УлаО-п), У (п)=Ула(0)т;
- наборы из Кп пар координат соответственных точек смежных кадров п и п+1: |(1(п)к,1, ](п)к,1),(1(п+1)к,0, ](п1)к,0), к=1,..,Кп}, п=0,..,К-1. Индекс «1» означает, что точка кадра имеет соответственную точку на следующем смежном кадре; индекс «0» означает, что точка кадра имеет соответственную точку на предыдущем смежном кадре;
- значения I, I, ¡0, ]0, 1^;
- установочные углы датчика изображения (иу, уу, уу);
- ЦМР Z=Z(X,Y).
Требуется получить:
А. Параметрически заданные математические процедуры ф^ фу для вычисления координат ^У) точки местности в заданной ФСК по координатам (¡(п),](п)) ее образа на п-ом изображении:
X=фx(i(n),j(п)/0(п)); Y=фy(i(n),j(п)/0(п)), (1.1)
где 0 - вектор параметров функций фx, фу; 0(п) - вектор значений параметров 0 применительно к п- му изображению;
Б. Математическую процедуру для определения параметров 0(п), обеспечивающих минимум функции ^(0(О),., 0(М):
IV — 1 Л. I
5(0(О|,..,0(~>)=
(ф (£.’, Л'1/0м )-Ф ( £Т, ІЇГ/ 0'"+,)))
+
(1.2)
=0 к=1
=1 ф ('кї, 1? / 0м )-Фу (С11, ІІТ/ 0(“*,)))
I[{фу ('?”)./0(П))-)2 +(фу (І'0,Л<“) /0(“))-)2).
М=П \ /
Геометрическая иллюстрация
На рис. 1. 1 показан набор исходных кадров видеоряда, на которых обозначены пары соответственных точек (светлые и темные кружки). Второй индекс «1» означает, что точка кадра имеет соответственную точку на следующем смежном кадре; индекс «0» означает, что точка кадра имеет соответственную точку на предыдущем смежном кадре. Также показано прогнозное положение надирных точек, координаты которых имеют индекс «н». При наличии параметрически заданных процедур фх, фу пересчета координат точек изображения в координаты их прообразов на местности можно отобразить области захвата, а также обозначенные на кадрах точки в контексте равномасштабного фотоплана.
2
Варьирование параметров 0 для каждого кадра позволяет трансформировать и смещать области захвата вместе с указанными на рис. 1. 1 точками. Оптимальными для каждого кадра значениями параметров являются такие значения, при которых сумма квадратов длин, показанных на рис. 1.1 отрезков по всему видеоряду, минимальна.
Порядок _решения задачи
На первом этапе, используя фотограмметрическую модель кадрового датчика изображения, требуется определить вид процедуры (1.1) и уточнить список параметров, подлежащих коррекции.
На втором этапе требуется решить задачу на определение экстремума функции многих переменных (1.2) и получить выражения для вычисления оптимальных параметров процедуры
(1.1) применительно к каждому кадру видеоряда. Одновременно требуется определить минимально необходимое число пар соответственных точек смежных кадров
На третьем этапе требуется определить поправку на смещение координат X,Y единую для всех кадров видеоряда и связанную с систематической погрешностью измерения координат в течение времени формирования видеоряда.
На рис. 2.1 изображена схема формирования кадра дискретного изображения. В случае дискретного изображения связь координат точек на местности и на изображении определяется модифицированными на случай цифрового изображения прямыми классическими фотограмметрическими уравнениями (ФГУ) следующего вида:
2. Предварительный переход к эквивалентным плановым снимкам
) С11(І0 - І) + С12('о - 0 - С1з/э .
С31(/0 - І) + С32('о - 0 - С3з/э ’
(2.1)
У - Ула = (7 (X, У) - 7ЛА )
С21 (/0 І) + С22('о '-) С2з1э
С31(/0 - І) + С32('о - 0 - С3з/э
Кадр
изображения
Связь координат изображения точки местности и координат ее прообраза на местности определяется обратными ФГУ следующего вида:
-• _ , = f C11(X - X ЛА ) + C21(Y - YIA ) - C31(Z - Z ЛА ) . (22)
J J° Ja Ci3(X - Xда ) + C23(Y - Yia ) - C33(Z - Zia )’
' - ' = f C12 (X - XЛА ) + C22 (Y - YIA ) - C32(Z - ZIA )
0 f э Ci3(X - Xia ) + C23(Y - Yia ) - C33(Z - Zia )’
где: ^^cosy-cosu; ^^siny-cosy-cosy-sinusiny; св—siny siny-cosy-sinucosy;
с^-sinycosu; ^^cosy-cos^^ny-sinu-siny; ^^siny-sinu-cosy-cosy siny; (2.3)
^^sinu; ^^cosu-siny; ^^cosu-cosy;
При этом входящие в (2.3) эквивалентные углы, определяющие ориентацию снимка в пространстве, вычисляются по формулам:
* *
C2 • * C32
y = arCtg —*^; u = arcsin C31; y = arCtg-,
C11 C33
где с*у - элементы матрицы С*=СлаСу.
Решение поставленной выше задачи затрудняется по следующим причинам:
- относительная сложность нелинейных ФГУ. Это затрудняет формирование целевой функции вида (1.1) и поиск ее экстремума;
- необходимость в случае выраженного рельефа применения итерационной процедуры определения X,Y. Это ставит под сомнение принципиальную возможность реализации самого подхода;
- снижение эффективности компьютерной поддержки поиска пар соответственных точек вследствие возможных существенных искажений изображения, связанных с угловой динамикой носителя.
Для преодоления указанных трудностей предлагается введение процедуры предварительного трансформирования (ортотрансформирования) кадров с целью их приведения к пространственно ориентированным эквивалентным плановым снимкам (ЭПС) на основе имеющейся пилотажной информации.
Эквивалентный горизонтальный снимок (ЭГС) - изображение {xir,jr}, полученное при строго нулевых углах ориентации снимка (u, y, y) в выбранной ФСК. Ортонормированный ЭГС -строго равномасштабный ЭГС (2.5). Для горизонтальной местности эти снимки идентичны, причем ФГУ имеют вид:
X - Xia =-( j, - А„)-**• Y - Yia =-( 'г - ',. К- («)
где k* - постоянный масштабный коэффициент, размерность м/пиксел.
В общем виде зависимость наклонного и горизонтального снимков можно записать следующим образом:
Ч ° ^ - 1.н = ^( i,j/p ); Ai ° .¡г - .¡г.н = .¡г( i,j/p ); (2 6)
Di ° i - io = i-1( i,j/p ); Aj ° j - jo = j-1( i,j/p ),
где, например, функции i^ij/p) и jr(i,j/p) образуются приравниванием правых частей (2.5) и
(2.1) в отсутствие выраженного рельефа (Z^Y^Z^H^onst) и путем итерационного опреде-
ления X, Y из уравнений (2.1) с последующей подстановкой их в равенства (2.5), в случае выраженного рельефа; p - вектор параметров преобразований: p =(u, y, y, Н)т.
ЭПС будем называть изображение, полученное по ФГМ вида (2.6), в которых вместо точ-
(2.7)
ных параметров р подставлены приближенные измерения ро=(ио, уо, Уо, Но)т параметров р:
Ч ° 1п - 1п.н = I(А 1, А]/р0 );
А1 ° 1 - 1о = 1 Г1 (А1 п ,А]п /р0 );
А)п ° ]п - ]п.н = ]Г (А ^ А]/Р );
А Л /р0
А] ° ] = Г(Ч,А]п/ро ) ,
где (1п,н, ]п,н) - положение надирной точки на плановом снимке.
Формулы (2.7) представляют собой математическую модель трансформирования исходного кадра в эквивалентный плановый снимок.
3. Определение вида целевой функции
Ниже представлены два возможных варианта аппроксимации (1.1).
А. Коррекция масштаба т, разворота в и сдвига (ЗХ=(х&, <У=(у&) ЭПС в масштабе местности
Коррекция определяется формулой:
✓ ^ \ ^ в _Бтв^(А/ ^ л
= 1 + т
V У ) V у )
Бт в соб в
+
1 )
(¿н у
или
(3.1)
V"
Х=Хн+ +аА1п-ЬА]п;
У=Ун+ Wy•8t +ЬА1п+аА]п,
где: а°шсо80; Ь°швт0; 8t - рассогласование моментов времени регистрации кадра изображения и используемого для этого кадра измерения ПНП.
Приравнивая правые части уравнений (3.1) для X и для У, записанных для пар соответст-
N-1
венных точек кадров п и п+1, получим набор из 21 Кп равенств для определения 3(№+1) пара-
„=о
метров 0(п)=(ап, Ьп, 5^)т (п=о,.. ,,К) (индекс “п” для краткости опущен):
Хп+ Wx,n•5tn +апА1(п)к,1 -ЬпА](п)к,1=Хп+1+ Wx,д+l•8t0+l +ап+1А1(п+1)к,о -Ьп+^+^о; (3.2)
Уп+ Wy,n•8tn +ЬпА1(п)к,1+апА](п)к,1=Уп+1+ Wy,n+l•8tn+l +Ьп+1А1п(п+1)к,о+ап+1А](п+1)к,о.
Так для N>3 для получения решения достаточно иметь не менее, чем по две пары соответственных точек на смежные кадры.
Для предотвращения накопления систематической погрешности в определении 8t в случае статистического решения задачи целесообразно для каждого п ввести ограничение вида: Хп,н=Хп, Уп,н=Уп, где (Хп,н,Уп,н) - расчетное положение прообраза надирной точки п-го кадра:
Хп,н=Хп+8^х,п; Уп,н=Уп+5^у,п.
Целевая функция (1.1) задачи оптимизации (3.2), (3.3) в этом случае имеет вид:
5(0(о),..,0от)=
(Х„ + 1¥х,Л + а„А4? - Ь, Д/<? - X,* - (.„А* - я„,Л4Г + Ь^АкГ)’ + (у,+(, А+ь„а,{,+а„ а/;„> - г„„ - (..„а* - ь^А/кГ - ч^А/к"4 )2
N
I (Я)2 ((2„+к).
(3.3)
N-1 К„
'II
„=о к =1
+
(3.4)
„=о
Получаемая в процессе минимизации относительно 0(о),..,0(^ система линейных уравнений легко решается методом Г аусса.
Б. Коррекция угловой ориентации ЭПС
По аналогии с (3.1) и с учетом (2.5) коррекция определяется формулами:
Х=Хн+Л|пХп-Л.]г.1, ■ к*,У=У„+ ^"Хп-А^ к*.
где А1Г.л, А]Г.л, вычисляются путем линеаризации модели 2.6, записанной для ЭПС:
Ді
(ДІп , ДІп Кп ) @ *г (Діп> ДІп 'К ) + Т! (Д/” ’ Дп іЯ" )
. Кк)
( Э/ (Діп , Д/п ІКп ) ^
К);
ж
(к) п
ж
(к)
У к=п
(3.5)
А]г.л = (Ап . А/п /рп > @ ]г (А/п . А/п /р ) + I
к
Рп =(8и, 8у, 8у, Ге, к*,Н+8Н, 1п,н+81. ]п,н+8])т; роп =(о, о, о, Ге, к*, н, 1п,н, ]п,н)т;
8рп =(8и, 8у, 8у, о, о, 8Н, 81(8и, 8у, 8у), 8](8и, 8у, 8у))т.
При формировании вектора параметров 0 для п-го кадра следует принять во внимание, что погрешности 81, 8] определяются углами (8и, 8у, 8у), поэтому целесообразно их не учитывать: 0(п)= (8и(п), 8у(п), 8у(п), 8Н~(п), 8t(n))т. Таким образом:
Н * ( Ді
к-
п| 1+
V ^ I Н *
ёН Л ДіпДІП о, ^ ДІП
. -ёь- — ёу-
Ді
Л
1 + ^П Г2
Г У У
ёу
(3.6)
1г.л
Н
(
к
ДІП (Л , ёН Л Діп ДІП
1 + -
V Ге V Н *
К
ёу + -Д-П- ёу-/е
1+
д/2 Л л
І ёь
Г2
Ге У у
Приравнивая правые части уравнений (3.5) для Х и для У, записанных для пар соответст-
N-1
венных точек кадров п и п+1, пол учим набор из 21 Кп равенств для определения 5(№+1) пара-
„=о
метров 0(п)= (8и(п), 8у(п), 8у(п), 8Н~(п), 8t(n))т ("=0,...^> (индекс “п” для краткости опущен и введены нормированные координаты А1~°А1/Ге; А]-°А]/:Ге):
Xn+8t(n)W(n)x - Н*[А]~(п)к 1(1+8Н~(п)) -Л1~(п)к 1А]~(п)к 18у(п)+А1 ~(п)к 18у(п) - (1+А]~(п)2к 1>8и(п)] = =Хп+1+8^п+1%(п+1)х - Н*[А]~(п+1)к,о(1+8Н~(п+1)) - А1-(п+1)к,оА]-(п+1)к,о8у(п+1)+А1-(п+1)к,о8у(п+1) -
(1+ЛГ(п+1\о)8и(п+1)];
Уn+8t(n)W(n)y - Н*[А1~(п)к 1(1+8Н~(п)) - А1~(п)к 1А]~(п)к 18и(п) - ЛГ(п)к 18у(п) - (1+АГ(п)\ 1>8у(п)] = =Уп+1+8t(n+1)W(n+1)y - Н*[АГ(п+1)к,о(1+8Н~(п+1)) - А1-(п+1)к,оА]-(п+1)к,о8и(п+1) - А]-(п+1)к,о8у(п+1) -
- (1+АГ(п+1 )2к,о)8у(п+1)], к=1,..,Кп, "=0,..,N-1.
Целевая функция (1.1) в этом случае имеет вид (символы “~” для краткости опущены):
N-1 К„
^(0(о),..,0(^)= I I (Хп+81;(п^х,п-Н*[А](п)к, 1 (1+8Н(п))-А1(п)к, 1 А](п)к, 18у(п)+А1(п)к, 18у(п)- (1+А](п)2кд)8и(п)]-„=о к=1
-Xn+1-8t(n+1>W(n+1)x+H*[Аj(n+1)k,0(1+8H(n+1))-Аi(n+1)k,0Аj(n+1)k,08У:n+1)+Аi(n+1)k,08y(n+1) - (1+А](п+1)2к,о)8и(п+1)])2+
N-1 К„
+ II (У"+81(п>№у,п-Н*[А1'(п)кД(1+8Н'(п)) -ЛГ(п)к,1Л]~(п)к,18и(п) -А]~(п)кД8у(п) - (1+Л1(п)2к,1)8У:п)]-Уп+1-„=о к=1
-8t(n+1>W(n+1)y+H*[Аi(n+1)k,0(1+8H(n+1))-Аi(n+1)k,0Аj(n+1)k,08u(n+1)-Аj(n+1)k,08y(n+1)- (1 +Л1(п+1)2к,0)8у(п+1)] )2+
N 2
+т (*«) (г,2,,+ж, і).
,=0
к
о
Для достижения определенности системы достаточно находить по три пары соответственных точек для каждой пары смежных кадров видеоряда.
Заключение
Особенности представленного подхода к привязке:
- в отличие от традиционного подхода, когда опорные точки определяются на снимке и на карте, здесь опорные (соответственные) точки определяются на смежных снимках. В связи с этим:
- обеспечивается высокая надежность автоматического определения соответственных точек корреляционно - экстремальным методом;
- не требуется наличия подробной карты с опорными ориентирами для каждого снимка видеоряда;
- информация о параметрах привязки снимков к местности берется из данных о сквозной взаимной привязке кадров в контексте фотоплана. Это позволяет повысить точность привязки, учитывая то обстоятельство, что относительные измерения пространственного положения ЛА могут быть значительно более точными, чем абсолютные;
- использование надирных точек кадра в процедуре не допускает накопление ошибок, связанных с дискретностью изображения.
ЛИТЕРАТУРА
Лобанов А.Н. Фотограмметрия. - М.: Недра, 1984.
THE MATHEMATICAL MODEL FOR INDIRECT OPTIMAL FITTING OF AIRBORNE MULTIFRAMING OPTICAL SENSOR IMAGES TO THE MAP IN REMOTE SENSING OF THE EARTH
Stepashkin V.N., Gritsenko A.E.
The mathematic relations is presented for corrections to photogrammetric model parameters values for airborne multiframing optical sensor image provided optimal coordinate binding for each frame of image videoset without using of reference ground points.
Сведения об авторах
Степашкин Владимир Николаевич, 1959 г.р., окончил МФТИ (1983), ВВИА им. проф. Н Е. Жуковского (1988), кандидат технических наук, доцент кафедры приборного, высотного и оптикоэлектронного авиационного оборудования ВВИА им. проф. Н.Е.Жуковского, автор более 100 научных работ, область научных интересов - дистанционное зондирование Земли, обработка многомерных сигналов (изображений), распознавание образов, фотограмметрия.
Гриценко Александр Евгеньевич, 1974 г.р., окончил Иркутское ВВАИУ (1996), адъюнкт кафедры приборного, высотного и оптико-электронного авиационного оборудования ВВИА им. проф. Н.Е. Жуковского, область научных интересов - цифровая обработка изображений, мониторинг земной поверхности.
