CC BY
73
УДК 551.4:681.3.06
T.B. Байкалова, U.U. Дубина, A.B. Евтюшкин, В.Н. Шишков
Алгоритм восстановления рельефа местности по стереопаре
Введение
В настоящее время процесс построения рельефа местности основан на хорошо разработанном и применяемом на практике стерео-фототопографическом методе создания топографических карт. Весь процесс работ осуществляется на одном из серийных универсальных фотограмметрических приборов типа СПР, СЦ и т.д. Измерения проводятся по стсреомодели сфотографированной местности, полученной по законам проективной геометрии с двух перекрывающихся снимков, что требует большого количества затрат рабочего времени, так как эта процессы не автоматизированы. В современных пакетах по обработке изображений для ЭВМ данный процесс всегда предполагает наличие базы данных высот характерных точек местности, а построение изолиний осуществляется методами аппроксимации и сглаживающих фильтров.
Алгоритм восстановления рельефа по цифровым моделям местности, представленный в данной статье, основан на построении марковского случайного поля и не имеет аналогов в широко распространенных пакетах по обработке изображений. Нахождение области перекрытая и базиса фотографирования производится автоматически поиском максимального коэффициента корреляции между снимками. Восстанавливаемое сечение зависит от качества аэрофотоснимков, масштаба съемки и точек оцифровки на дюйм. Разработанное программное обеспечение позволяет полностью автоматизировать весь процесс работ.
Алгоритмы построения цифровых
моделей меетноети
Построение карты высот местности по стереопаре требует идентификации на парных аэрофотоснимках одинаково расположенных точек местности. Для этого производится корреляционный анализ левого 5](.) и правого ¿>2 (•) снимков на окне фиксированного размера.
, . е - а ■ Ь / п
г (x1, у,, р) =-;
[(с -а2/п){й - Ь2/п)] ■ 1/2
a = X Si(i + j +
b =X S 2(i + X2' j + У2);
c=X Si(i+xi' j+yi) ■2;
d = X S +x2' J+У2) ■2;
e = X Si(i+xi'J+yi) ■S 2(i+x2'J+y 2);
x2=x 1 • .v„ • /»; 1/2=1/1 +1/0,
где суммирование ведется по фиксированному окну (переменные i и j); 5[(.) и S->(.) — значения яркостей снимков; п — число точек в окне; (,г-|,г/|), (Х2,У2) — координаты коррелируемых точек; .>•(,. у о — базовое относительное смещение снимков.
Для каждой точки первого снимка вычисляется коэффициент корреляции между ее окрестностью и окрестностью точки второго снимка для всего набора возможных параллаксов. Возможные параллаксы определяются по предварительной оценке минимальных и максимальных высот и разрешению сканера. Число возможных значений параллаксов можно оценить по следующей формуле:
N = (Hmax -H;1П)B ■ DPI H ■ M ■ 0.245 '
где H — высота съемки в метрах; М — масштаб съемки; В — базис в метрах; 1)14 — разрешение сканера точек/дюйм.
Тогда за искомое значение параллакса можно принять оценку р, обеспечивающую максимум корреляционной функции при фиксированном (Л'|.//| ).
Однако на реальных снимках существуют слабо коррелируемые участки (коэффициент корреляции 0.5 и меньше). Ситуация может осложняться наличием периодических структур, имеющих несколько корреляционных максимумов. Для таких ситуаций одного корреляционного критерия недостаточно. Дополнительным критерием может являться гладкость построенной модели рельефа.
Учет гладкости рельефа требует введения условной вероятности параллакса в зависимости от значений параллаксов в соседних точках, т.е. описания марковскими случайными полями (МСП) [2,3,4]. Один из вариантов построения марковской модели восстановления высот рассмотрен в [4].
Для построения МСП задастся вероятность параллакса р в точке (х,у), исходя из двух моделей, заданных в виде гиббсовских потенциалов Гч. и V;(:
exp(aVp(x y p) + bVs(x y p))
P( x y, p):
X exP(aVp (x У, pt) + bVs(x y, pi)
p
Гладкость рельефа обеспечивается гиббсов-ским потенциалом следующего вида:
К(х у, р) = X Л(Р' р' (х +'' у + ■/))'
Л(Р\, Р 2 ) = (Р\ - Р2 ) 2,
где р'(х,у) — значения параллаксов в точке (х.у): /. У принимают значения от -1 до +1 за исключением случая ¿=0 и ,)=0; \7р(х.у.р) — гиббсовский потенциал, зависящий только от коэффициента корреляции для значения параллакса р в точке (х.у).
Последовательная релаксация значений р в соответствии с заданными вероятностями и генерацией последовательности карт параллаксов р'(х.у) является марковским процессом. При постоянных а и 6, носящих название обратной температуры, марковский процесс является стационарным и имеет предельное распределение вероятностей, не зависящее от начального состояния.
Начальное состояние может быть любым, например, являться независимым шумом, но для ускорения сходимости марковского процесса К Предельному СОСТОЯНИЮ ПрСДПОЧТИТСЛЬНСЙ указывать в качестве начального некоторое приближенное значение, например, полученное в результате решения задачи без учета условия гладкости рельефа. В результате нескольких порождений р'(х.у) можно вычислить статистические оценки параллаксов для каждой точки снимка. Выбирая в каждом элементе растра в качестве решения значение, обеспечивающее максимум вероятностей, получаем решение исходной задачи.
Переход от значений параллаксов в каждой точке снимка к нормальным высотам и расчет прямоугольных координат проекции Гаусса Крюгера точек местности осуществляется по следующему алгоритму:
1. Угловые элементы внешнего ориентирования а |, Ю |, С I левого снимка принимают равными нулю. После этого определяют элементы взаимного ориентирования снимков а Ю [', «2', Ю2' и С 2 п0 известным в фотограмметрии формулам. Задача решается методом приближений. В первом приближении косинусы неизвестных принимают равными единице, а произведения синусов — равными нулю.
2. Дирекционный угол Т и угол наклона V базиса фотографирования получают из матрицы поворота:
А(Т V ) = А(а,Ю|,£,)* А(а,',£|'),
где элементами матриц А(«|,ЮьС|) и А(а[',С| ) являются направляющие косинусы, определяемые по формулам при (Х = а,\,
Ю = С = С| и а = ах\"ю =о, С = С
соответственно:
а | =со8а*со8%-8та*8тю*8т%, а 2=-сова* вт% - вта* вт ю* сов%, аз=-8та*со8Ю,
Ь]=С08Ю*8Н1%,
Ь2=С08Ю*С08%, ( 1 )
Ь^—этю,
С|=8та*со8%+со8а*8тю*8т%, С2=-8та*8т%+со8а*8тю*со8%,
СЗ=С08а*С08Ю.
3. Углы наклона второго снимка определяют из матрицы:
А(а2,ю2,%2) = А(т,п) * А(а2',ю2',х2'), где элементы матрицы А(а2',Ю2',%2') находят из формул (1) при
а=а2', ю=ю2', %=Х2-
4. Координаты правого центра проектирования вычисляют по формулам:
Л*2 Л.V| Их. Г*2 Г«, Ну. Х*2 /л-. И:.
где
Вх=В*со$п*со&т, Ву=В*со&п*$\т, В2=В*&тП~,
Лл-|. 2а>[ — координаты левого центра
проектирования; В — базис фотографирования.
5. Вычисляют координаты точек стереопары, соответствующие их значениям на горизонтальных снимках, используя формулы:
.Г2„ ;/*Л2\.Х2'.
где — пространственные координаты
точки левого снимка, определяемые по формулам:
Л" х"*п 1~у"*П'2:Г*ч:>. =х * Ъ |+'у * 6 2 -/* Ъ з, (2)
'/.' Д'*Г|~|/*Г2
в которых направляющие косинусы <1\.<1->.....'';(
являются элементами матрицы А(а|,Ю|>%|)> Г2', Ъ-2 — пространственные координаты точки правого снимка, определяемые по формулам (2), в которых направляющие косинусы являются элементами матрицы А(а2,Ю2,%2); / — фокусное расстояние.
6. Вычисляют фотограмметрические координаты точек цифровой модели, используя формулы:
Л Л"л-|—<!Л". Г......)л-|—<!)". 7.......7.нх-р-&7..
где (I V Л*.С|(|. (!) Л*,/и„ <17......
Л'= (Вх+х20*Вг/Л (.с,0-.г-20); р - величина параллакса, полученная из МСП.
Аналогично строят вторую цифровую модель по второй стереопаре и т.д. Последовательное соединение заключается в приведении масштаба последующей модели к масштабу предыдущей по связующим точкам, расположенным в зонах тройного перекрытия.
Эта математическая модель является свободной, и она не ориентирована относительно системы координат, используемой в геодезии.
7. Для внешнего ориентирования модели используют опорные точки. Обычно плановое положение опорных точек задают координатами
в проекции Гаусса-Крюгсра, а высотное — нормальными высотами. Поэтому предварительно переходят от прямоугольных координат в проекции Гаусса-Крюгера и нормальных высот к пространственным прямоугольным координатам. За начало пространственной прямоугольной системы X'YZ' принимают опорную точку, расположенную примерно в центре обрабатываемого участка. Определяют криволинейные геодезические координаты по прямоугольным координатам Гаусса-Крюгера, используя известные в сфсроидичсской геодезии формулы, а затем пространственные прямоугольные координаты Х\ F, Z'. В случаях когда ось абсцисс пространственной системы принимается параллельной плоскости осевого меридиана, систему XYZ' необходимо повернуть вокруг оси OZ' на угол Со, равный сближению осевого и проходящего через начало координат меридианов. После поворота на угол Со имеем:
А"=А" *cosG()-F *sin С0, Y= F:*cosC0+A":*sm С0. Z /.'.
По разностям геодезических и фотограмметрических координат определяют элементы внешнего ориентирования цифровой модели, которыми являются: Лп. F0, Z0 — пространственные прямоугольные координаты начала фотограмметрической системы, V — продольный угол наклона модели, X — по~ перечный угол наклона модели, V — угол поворота модели, I — знаменатель масштаба модели.
Связь между координатами точки X, Г, Z в
тгтг 171 Т 71 гтг I 1
системе A1Z и А , 1 , Z в фотограмметрической системе XYZ' выражается формулами:
A=Ao+<L¥'=Ao+ (а, *X'+a2*Y+aH*Z') %
Y= у а+d F=Yо +(b\ *X'+b2*Y+bH*Z') */. (3) Z=Z0+dZ'=Z0+ (с | *X'+c2*Y+cH*Z' ) */. где (IV. d)". dZ — приращения координат определяемой точки относительно начала координат; а [,...,— направляющие косинусы, вычисляемые по формулам (1) при a = V ,
W=V , X=X ■ Одна опорная точка дает возможность составить три уравнения с семыо неизвестными. Поэтому для ориентирования модели необходимо не менее трех опорных точек, из которых две должны иметь все три координаты, т.е. X, F, Z, а третья — только Z.
8. Координаты модели, исправленные за деформацию, можно определить, используя степенные полиномы.
9. Для перехода от фотограмметрической пространственной системы координат к координатам в проекции Гаусса-Крюгера и нормальным высотам используют известные в сфсроидичсской геодезии формулы.
Выбор подобной схемы расчетов позволяет легко модифицировать алгоритм путем ввода дополнительной модели, например, коррекции
по опорным точкам с известными координатами и высотами.
Программная реализация
Для практического использования описанного выше алгоритма построения изолиний рельефа по цифровой модели местности разработан комплекс программ для ЭВМ типа РС/АТ-486 и RS-6000 па языке С++. Этот комплекс состоит из программы расчета и построения карты параллаксов, полученных методом МСП, программы расчета элементов внешнего и взаимного ориентирования пары снимков, программы расчета координат и нормальных высот элементов изображения и его пересчета в картографическую проекцию Гаусса-Крюгера, программы оценки точности и коррекции полученных результатов.
Апробация
Разработанное программное обеспечение апробировалось на аэрофотоснимках территории, расположенной в Кемеровской области и примыкающей к Беловскому водохранилищу, где происходят интенсивные эрозионные процессы. Аэрофотоснимки сканировались на черно-белом полутоновом сканере с разрешением 300 dpi.
В результате была получена картосхема исследуемой территории. Основные контуры местности построены по фотопланам ВИСХАГИ, а изолинии рельефа — по продольным и поперечным параллаксам, полученным в области перекрытия пары аэроснимков (рис. 1). На данной картосхеме также изображены визуально отдешифрированные линии поверхностного смыва почвы, которые совпадают с характерными изгибами горизонталей. Расхождения в изображениях изолиний рельефа при совмещении полученной картосхемы с фотопланом достигают 1 мм в центре и 2 мм по краям рабочей площади.
Координатная и высотная привязка, а также оценка точности построения контуров и изолиний рельефа осуществлялась по 5 опорным точкам, равномерно расположенным в зоне перекрытия. Разность в значениях прямоугольных координат и нормальных высот точек, полученных в результате построения математической модели местности и измерения по топографической карте, не превышает 10 м в плане и 50 м по высоте, что соответствует точности карты масштаба 1:100000 на данную территорию.
Заключение
Описанная математическая модель и разработанное программное обеспечение позволяет восстанавливать изолинии и значения абсолютных высот местности по парс аэрофотоснимков, вычислять прямоугольные координаты исследуемого участка в проекции Гаусса-Крюгера с необходимой точностью для решения ряда тематических задач. Полученные цифровые карты в виде матрицы высот или векторных
изолиний позволяют существенно расширить мационных системах и технологиях опсратив-иримснснис аэро- и космоснимков в гсоинфор- ного цифрового картографирования.
Литература
1. Куштин И.Ф., Бруевич П.1Г., Лысков Г.Л. Снра-вочник техника фотограмметриста. М., 1988.
2. Гимельфарб Г.Л., Залесный A.B. Цифровая обработка изображений, представленных моделями марковских случайных нолей. Киев, 1991. 43 с. (Пренр./ЛН УССР 11и т кибернетики им. В.М.Глушкова; 91-42).
3. Залесный Л.В., Юшаков В.1Г. Эксперименты но генерации изображений // Теоретические и прикладные вопросы распознавания изображений. Киев, 1992. 15-26 с.
4. Gimelfarb G.L. Cibbs Random Fields and Compound Bayesian Decisions at the Lower Lever of Digital Image Proeessing//Pallern Recognition. 1991. N.l P.39-49.
