СПОСОБ ПОСТРОЕНИЯ ЦИФРОВОЙ МОДЕЛИ ПОВЕРХНОСТИ ПО ДАННЫМ КОСМИЧЕСКОЙ СТЕРЕОСЪЕМКИ
METHOD FOR BUILDING A DIGITAL SURFACE MODEL BASED ON DATA FROM SPACE STEREO IMAGING

Изобретение относится к области картографии и может быть использовано для повышения точности и детальности построения цифровой модели поверхности по данным космической стереосъемки.
Построение цифровой модели поверхности (ЦМП) по данным космической стереосъемки является одной из актуальнейших задач анализа данных дистанционного зондирования Земли, поскольку позволяет извлекать из снимков высотную информацию, представляющую собой ценный источник информации для различных прикладных задач. Однако стереобработка космических изображений сопряжена с рядом фундаментальных трудностей, таких как: высокая изменчивость поля яркости, вызванная атмосферными и иными явлениями; малый размер объектов интереса, сопоставимый с геометрическими искажениями, вызванными разницей углов съемки; сложная траектория движения космического аппарата и специфический способ получения снимков; существенный объем данных для обработки.
Точность ЦМП определяется точностью определения высот объектов местности (высот рельефа и относительных высот зданий и сооружений) и точностью положения объектов в плане. Детальность ЦМП определяется наличием относительно небольших по площади (порядка десятков кв. метров) и по высоте (порядка единиц метров) объектов, выступающих на фоне рельефа или на фоне плоскостей более крупных объектов. Точность и детальность ЦМП может быть измерена путем сравнения с эталонной моделью поверхности, построенной по результатам воздушного лазерного сканирования и имеющей более высокую (сантиметровую) точность в плане и по высоте.
Важнейшим вопросом при создании способа построения ЦМП является выбор алгоритма стереосопоставления. В подавляющем большинстве современных способов и пакетов программ авторы используют алгоритм [1] Semi-Global Matching (SGM). Данный алгоритм совмещает в себе высокую вычислительную эффективность с не менее высоким качеством построенной карты глубины. Однако, эта эффективность достигается на стереопарах сантиметровой точности, получаемых при аэросъемке, которые обладают более детализированной и предсказуемой структурой по сравнению с космическими стереопарами. В то же время точность и детальность работы данного алгоритма на космических стереопарах, получаемых, например, системой Pleiades [2], имеющих пространственное разрешение порядка 0.5 м, недостаточны. Альтернативой данному методу является метод разрезов графа GraphCut (GC) [3], который благодаря своей глобальной структуре лучше подходит для космических снимков. Однако, он не обладает той же чувствительностью к малым изменениям глубины, то есть к микрорельефу, что и алгоритм SGM.
Описание алгоритмов GC и SGM нахождения сдвигов между соответственными точками эпиполярно выравненных изображений приведено в Приложениях 1 и 2 соответственно.
Прототипом предлагаемого способа построения ЦМП является работа [4], в которой предлагается способ построения ЦМП, а также приводятся результаты сравнения качества полученной ЦМП для различных алгоритмов стереосопоставления.
Представленная в прототипе последовательность операций состоит из следующих шагов.
1. Загрузка входных данных, состоящих из панхроматического снимка с разрешением не менее 1 м и 4-х канального мультиспектрального изображения, включающего красный, зеленый, синий и инфракрасный каналы яркости, с разрешением не менее 2 м, а также RPC коэффициентов математического преобразования координат точек (пикселей) космического изображения на точки местности и обратно.
2. Предварительная обработка исходных данных, включающая выбор соответственных точек для коррекции RPC-коэффициентов и повышение пространственного разрешения мультиспектральных изображений до разрешения панхроматического.
3. Эпиполярное выравнивание изображений стереопары, обеспечивающее размещение соответственных точек в одноименных строках выравненных изображений.
4. Нахождение матрицы D сдвигов пар соответственных точек на стереоснимках через сопоставление стереоизображений поверхности (решение задачи стереосопоставления). Матрица сдвигов D соответственных точек состоит из элементов D(p), каждый из которых представляет собой расстояние в пикселях между точкой р первого изображения, и точкой р+D(p) второго изображения стереопары, причем, поскольку предварительно выполняется эпиполярное выравнивание, то смещение возможно только вдоль строки.
5. Построение цифровой модели рельефа и нормализованной ЦМП.
6. Создание ортотрансформированного изображения.
7. Классификация областей растительности, а также строений и сооружений
8. Извлечение и моделирование отдельных объектов (зданий и сооружений).
Недостатком данного способа при формировании ЦМП по космическим снимкам является недостаточно высокая точность и детальность, достигаемая используемыми алгоритмами.
Целью изобретения является создание способа построения более точной и высоко детальной ЦМП по данным космической стереосъемки.
Для достижения поставленной цели повышения точности и детальности ЦМП, в отличие от прототипа, предлагается:
- учитывать априорно известную цифровую модель рельефа (ЦМР) с разрешением не ниже 150 метров на пиксель;
- выполнять построение ЦМП по стереоснимкам, в которых априорно исключены из рассмотрения области, маскируемые облаками, а также водные поверхности;
- выполнять сегментацию панхроматических изображений и учитывать ее при стереосопоставлении снимков таким образом, чтобы избежать значительных скачков высоты внутри сегментов;
- находить матрицу  сдвигов пар соответственных точек на стереоснимках (задача стереосопоставления) за счет модификации алгоритмов GC и SGM и их комбинирования.
GC
CC
GC
data
GC
Модификация алгоритма GC состоит в определении для каждой точки первого и второго изображения вектор-признака DAISY(p) (дескриптора DAISY [5]) и замене в выражении для энергетической функций E(D), приведенной в приложении 1, меры близости E_(D) на
GC
где D - матрица сдвигов соответственных точек на стереоснимках, формируемая алгоритмом GC,
а
L2
2
|||| - норма вектора а в пространстве L;
L
p - точки первого снимка стереопары,
R
L
GC
L
p(p,D) - точки второго снимка стереопары, соответствующие точкам p и находящиеся с ними в одной строке.
R
L
L
L
Модификация алгоритма SGM состоит в замене обоих слагаемых в приведенной в приложении 2 энергетической функции : первого - меры близости  точек p второго изображения к точкам p первого изображения на , и второго слагаемого , представляющего собой штраф за несовпадение сдвига  в точке p и сдвигов в соседних точках q, находящихся в окрестности точки p фиксированного радиуса, на модификацию . Модифицированная энергетическая функция имеет вид
R
L
i
где  - вектор сдвигов  точек второго изображения p относительно точек p первого, находящихся на луче S, направленном в точку р первого изображения, для которой ищется сдвиг,
R
L
i
- модифицированная мера близости точек p второго изображения относительно точек p первого изображения, находящихся на луче S, направленном в точку р первого изображения,
S
L
L
i
L
- суммарный штраф, модифицированный за счет введения штрафных коэффициентов C(p,q) за несовпадение сдвигов в точках p, расположенных на луче S, и сдвигов в соседних с p точках q из окрестности фиксированного радиуса,
S
L
L
C(p, q) - коэффициент сегментации, который равен 1, если яркости точек p и q в сегментированном изображении не равны, и 3 -  в противном случае,
L
i
- сдвиги точек второго изображения, соответствующих точкам p и q из первого изображения, находящихся на луче S,
L
i
i
p ∈ S - точки первого изображения, расположенные на луче S,
1
2
С и С - штрафы за несовпадение сдвигов в соседних точках,
L
i
i
L
q ∈ N(p)∈S - множество точек, расположенных на луче S, из окрестности фиксированного радиуса с центром в точке p,
Т(⋅) - функция, которая возвращает 1, если условие в скобках верно, и ноль в противном случае.
i
L
i
i
Комбинирование модифицированных алгоритмов заключается в том, что по матрице  которая доставляет минимум функции , находят грубые оценки сдвигов соответственных точек второго изображения относительно точек первого. Далее локально в точке р и входящих в нее лучей S в интервалах сдвигов, которые определяются сдвигами  (где p ∈ S, записанными в матрице , по минимумам энергетических функций  для всех лучей S находят в точке р первого изображения точную оценку сдвига  соответственной ей точки второго изображения.
Более подробно сущность предлагаемого способа поясняется дальнейшим описанием и чертежами:
на фиг. 1 представлен предлагаемый алгоритм построения ЦМП по данным космической стереосъемки;
на фиг. 2 представлен предлагаемый алгоритм стереосопоставления фрагментов изображений стереопары;
на фиг. 3 представлен пример эпиполярно выравненных фрагментов стереоснимков;
на фиг. 4 представлена схема расположения ключевых точек для вычисления дескриптора DAISY;
L
i
R
на фиг. 5 представлена геометрическая связь точки p на луче S с точками p в зоне поиска сдвига точки соответствия ;
на фиг. 6 представлены для сравнения два ЦМП: первое изображение, полученное предлагаемым способом по данным космической стереосъемки, второе изображение (эталонное ЦМП) по данным лазерной аэросъемки.
Предлагаемый способ построения ЦМП по данным космической стереосъемки в соответствии со схемой фиг. 1 и фиг. 2 производится в последовательности:
1. Получение входных данных (поз. 1 фиг. 1), включающих:
- пары стереоснимков местности в панхроматическом диапазоне с исключенными пикселями, соответствующими местоположению облаков и водной поверхности;
- априорно известную ЦМР с разрешением не ниже 150 метров на пиксель;
- RPC-коэффициенты математического преобразования координат точек (пикселей) космического изображения на точки местности и обратно;
2. Сегментацию изображений стереопары (поз. 2 фиг. 1) по яркости методом анизотропной медианной фильтрации для удаления шума и мелких выбросов яркости, включающую:
- выпуск из каждой точки изображения 16 лучей длиною 7 пикселей,
- вычисление дисперсии яркости вдоль каждого луча,
- вычисление медианы яркости точек вдоль луча с наименьшей дисперсией,
- присвоение точке панхроматического изображения найденное значение медианы,
3. Разбиение стереоснимков на одинаковые по размеру меньшие фрагменты (поз. 3 фиг. 1);
4. Эпиполярное выравнивание фрагментов изображений стереопары в панхроматическом диапазоне через повороты изображений относительно центров снимков так, чтобы смещения, вызванные сменой ракурса при получении стереопары, лежали вдоль одной строки (поз. 4 фиг 1). Пример эпиполярного выравнивания приведен на фиг. 3, горизонтальные линии вдоль снимков соответствуют строкам выровненных изображений;
5. Выполнение процедуры стереосопоставления (поз. 5 фиг. 1), в результате которой получают матрицу сдвигов соответственных точек. Процедура стереосопоставления включает в себя следующие этапы.
5.1 Вычисление дескриптора DAISY (р) для каждой точки р изображений стереопары (поз. 5.1 фиг. 2) по алгоритму:
- вокруг точки р панхроматического изображения строят набор из трех концентрических окружностей радиусами 3, 6 и 9 пикселей;
- на каждой из этих окружностей выбирают по 8 точек, равномерно распределенных вдоль окружности, геометрия размещения точек приведена на фиг. 4;
- в исходной точке р, а также в 8 точках окружности радиусом 3 пикселя выполняют свертку с гауссовым ядром модулей градиентов яркости по 8 направлениям, радиус гауссового ядра 2 пикселя;
- в каждой из 8 точек окружности радиусом 6 пикселей выполняют свертку с гауссовым ядром модулей градиентов яркости по 8 направлениям, радиус гауссового ядра 3 пикселя;
- в каждой из 8 точек окружности радиусом 9 пикселей выполняют свертку с гауссовым ядром модулей градиентов яркости по 8 направлениям, радиус гауссового ядра 4 пикселя;
- объединяют полученные значения в вектор размерности 200 (25 точек по 8 направлений градиента в каждой);
- нормализуют полученный вектор, который считают вектором-признаком (дескриптором) DAISY (р) точки р изображения.
5.2 Выполнение первого этапа стереосопоставления с помощью модифицированного алгоритма GC (поз. 5.2 фиг. 2) в последовательности:
- определяют диапазон поиска  сдвигов точек соответствия по цифровой матрице рельефа (ЦМР) как интервал
p
где k - коэффициент перевода 1 м высоты местности в пиксели смещения на стереопаре, вычисляемый индивидуально для каждого фрагмента по RPC-коэффициентам снимков,
- минимальная и максимальная высота рассматриваемой местности по ЦМР в пределах фрагмента;
- находят методом оптимальную матрицу сдвигов  такую, что значение  каждого элемента матрицы принадлежит интервалу , то есть , и для которой достигается минимум модифицированной энергетической функции :
L
R
где  - мера близости точек p первого изображения точкам p второго изображения,
L
L
L
p=(x,y) - точки первого снимка стереопары,
R
L,
GC
R
R
L
GC
L
L
L
GC
L
p(pD)=(x,y)=(x+D(p),y) - точка второго снимка стереопары, сдвиг которой относительно точки p первого снимка равен D(p) и находится в пределах интервала ,
GC
L
L
D(p) - сдвиг точки второго изображения относительно точки p первого изображения, находящихся в одноименных строках изображений и в пределах интервала ,
L
L
N(p)- множество точек в окрестности точки p фиксированного радиуса,
V - весовой коэффициент штрафной функции,
Т(⋅) - функция, которая возвращает 1, если условие в скобках верно, и ноль в противном случае,
GC
L
R
L
С - весовой коэффициент штрафной функции для условия D(p)=∅ (∅ - пустое множество) означающего, что сдвиг точки p относительно точки p не определен.
5.3 Выполнение второго этапа стереосопоставления (точного сопоставления фрагментов разбиения) с помощью модифицированного алгоритма SGM (поз. 5.3 фиг. 2) в последовательности:
L
определяют диапазон поиска  сдвигов для каждой точки снимка p как интервал
i
L
L
i
в точку р первого изображения от края изображения проводят 16 равномерно разнесенных по углу лучей (i=1…16). В дальнейшем тексте точки первого изображения, находящиеся на S луче, обозначены p. Расположение точек р и p на луче S показано на фиг. 5;
L
L
R
для i-того луча находят вектор  сдвигов  соответственных точек второго изображения, находящихся в одноименных с p строках, (связь точек p с зоной поиска  соответственной точки p иллюстрирует фиг. 5) следующим образом:
Rs
L
i
- рассматривают множество вариантов  вектора сдвигов  точек pвторого изображения относительно точек p ∈ S первого изображения, координаты которого - последовательные значения  (из заранее заданного диапазона)
- для s-того варианта вектора сдвигов  определяют значение энергетической функции :
L
i
Rs
где  - мера близости всех точек p первого изображения, расположенных на луче S, к точкам p второго изображения, имеющих сдвиги
S
L
L
C(p,q) - коэффициент сегментации, который равен 1, если яркости точек p и q в сегментированном изображении не равны, и 3 - в противном случае,
S
L
L
i
L
i
- модифицированный за счет введения коэффициентов C(p,q) суммарный штраф за несовпадение сдвигов в точках p, расположенных на луче S, и сдвигов в соседних с p точках q ∈ S из окрестности фиксированного радиуса,
L
i
i
p ∈ S - точки первого изображения, расположенные на луче S,
1
2
С и С - штрафы за несовпадение сдвигов в соседних точках,
L
L
i
q ∈ N(p) - множество точек в окрестности фиксированного радиуса с центром в точке p луча S,
Т(⋅) - функция, которая возвращает 1, если условие в скобках верно, и ноль в противном случае;
- находят методом динамического программирования среди множества векторов ) оптимальный вектор  сдвигов соответственных точек второго изображения, для которого достигается минимум энергетической функции
i
повторяют шаги определения вектора  для каждого луча S;
выбирают среди полученных 16-ти векторов-сдвигов  тот, для которого достигается наименьшее значение
выбирают в векторе  первую координату, отвечающую точке р первого изображения, значение сдвига, записанное в этой точке, является искомым точным значением сдвига  соответственной точки второго изображения;
повторяют вычисление сдвига  для всех точек р первого изображения и формируют для этого изображения результирующую матрицу сдвигов , в которой на месте яркости h(р) исходного изображения стоит значение
6. Построение ЦМП фрагментов изображения в географических координатах (поз. 6 фиг. 1) в последовательности:
p
вычисляют высоту местности в пикселе р по найденным сдвигам , базовой высоте  фрагмента и коэффициенту k перевода метров высоты в пиксели изображения:
по матрице h(p) вычисляют географические координаты точек поверхности Н(Х, Y) с помощью RPC-коэффициентов для того, чтобы получить ЦМП, привязанную к местности.
7. Объединение ЦМП всех фрагментов разбиения Н(Х, Y) в итоговую ЦМП (поз. 7 фиг. 1) и выполнение медианной фильтрации полученной матрицы высот.
Для проверки работоспособности способа был проведен вычислительный эксперимент по построению цифровой модели местности по данным космической стереосъемки. Было выполнено сравнение качества ЦМП, построенной предлагаемым способом, с ЦМП прототипа, использующего сопоставление стереоснимков по алгоритму GC. Сравнение проводилось по значениям среднеквадратического отклонения оценок высоты ЦМП от эталонной ЦМП, построенной по данным воздушной лазерной съемки. Сравнение производилось по областям, не отмеченным в ходе предварительной классификации как водная поверхность и облачный или растительный покров. Полученное в прототипе среднеквадратичное значение ошибки (СКО) 3,74 м было улучшено до 2,93 м., что составляет уменьшение СКО ошибки на 21,6%. При этом для полученной ЦМП также были оценены точность по высоте (среднее абсолютное отклонение высот объектов ЦМП от высот объектов на эталоне), которая составила 1,81 м и точность в плане (отклонение границ объектов ЦМП - строений от границ этих же объектов на эталоне), которая составила 1,87 м.
На фиг. 6 показаны фрагменты матрицы высот рельефа, вычисленных предложенным методом, и соответствующие им результаты лидарной съемки. Визуальное сравнение показывает, что вычисленные матрицы высот рельефа качественно близки к результатам лидарной съемки.
Представленный способ построения ЦМП может найти свое применение для решения широкого круга задач, которые удобно разбить на три направления.
Первым важным направлением является формирование эталонных ЦМП по космическим стереоснимкам, необходимых для навигации беспилотных летательных аппаратов. Вне зависимости от типа задачи, данный способ обладает рядом важных преимуществ в сравнении с построением ЦМП по аэросъемке: более широкое покрытие; возможность получения данных на территории, где использование беспилотных летательных аппаратов невозможно; меньшая стоимость; высокая оперативность получения требуемой информации.
Второе направление связано с управлением территориями, картографией, контролем темпов строительства и обнаружения объектов незаконного строительства, корректировкой кадастровой и градостроительной документации, оценкой геоморфологических изменений.
Третье направление - природный и экологический мониторинг, инвентаризация сельскохозяйственных угодий, оценка динамики изменения площади лесных массивов, создание тематических природных карт.
Приложение 1
Алгоритм GC поиска матрицы сдвигов между соответственными точками эпиполярно выравненных стереоизображений
L
R
GC
GC
SGM
Задачу поиска матрицы сдвигов D между точками первого p и второго p эпиполярно выровненных стереоизображений решают через нахождение минимума «энергетической» функции E(⋅). Значения элементов матриц D для алгоритмов GC и SGM определяются по-разному, поэтому в дальнейшем тексте обозначаются D для алгоритма GC и D для алгоритма SGM.
GC
GC
Энергетическая функция E(D) имеет вид:
GC
GC
GC
data
GC
GC
smth
GC
E(D)=E_(D)+E_(D),
L
R
где  - мера близости точек p первого изображения точкам p второго изображения,
L
L
- суммарный штраф за несовпадение сдвигов в точках p и сдвигов в соседних с pточках q из окрестности фиксированного радиуса,
L
L
L
p=(.x,y) - каждая точка первого снимка стереопары,
R
L
GC
R
R
L
GC
L
L
L
GC
L
p(p,D)=(x,y)=(x+D(p),y) - точка второго снимка стереопары, сдвиг которой относительно точки p первого снимка равен D(p),
L
R
- средние яркости первого и второго изображений стереопары в окрестностях точек p и p соответственно,
L
L
N(p) - множество точек в окрестности точки p фиксированного радиуса,
V - весовой коэффициент штрафной функции,
Т(⋅) - функция, которая возвращает 1, если условие в скобках верно, и ноль в противном случае,
GC
L
R
L
С - весовой коэффициент штрафной функции для условия D(p)=∅ (∅ - пустое множество) означающего, что сдвиг точки p относительно точки p не определен.
GC
GC
Минимум функции E(D) в алгоритме GC ищут для всех точек изображения одновременно специальным численным методом, который относится к классу глобальных алгоритмов.
Приложение 2
Алгоритм SGM поиска матрицы сдвигов между точками эпиполярно выравненных изображений
SGM
В алгоритме SGM поиск минимума функции E(⋅), в отличие от алгоритма GC, осуществляется не глобально для всех точек сразу, а для каждой точки р первого изображения отдельно и независимо от других, поэтому алгоритм относится к классу полу-глобальных. Для этого
i
i
L
L
L
- в точку р первого изображения от края изображения проводят 16 равномерно разнесенных по углу лучей S, в дальнейшем тексте точки первого изображения, находящиеся на луче S, обозначаем p=(x,y),
L
- для i-того луча находят вектор  сдвигов  соответственных точек второго изображения, находящихся в одноименных с p строках, следующим образом:
Rs
L
i
рассматривают множество вариантов  вектора сдвигов  точек pвторого изображения относительно точек p ∈ S первого изображения, координаты которого - последовательные значения  (из заранее заданного диапазона)
для s-того варианта вектора сдвигов  определяют значение энергетической функции
L
i
Rs
где  - мера близости всех точек p первого изображения, расположенных на луче S, к точкам p второго изображения, имеющих сдвиги
L
i
L
i
- суммарный штраф за несовпадение сдвигов  в точках p, расположенных на луче S, и сдвигов  в соседних с pточках q ∈ S из окрестности фиксированного радиуса,
L
i
i
p ∈ S- точки первого изображения, расположенные на луче S,
1
2
С и С - штрафы за несовпадение сдвигов в соседних точках,
L
L
L
q ∈ N(p) - множество точек первого изображения в окрестности точки pфиксированного радиуса,
R
Rs
Rs
q ∈ N(p) - множество точек второго изображения в окрестности точки pфиксированного радиуса,
Т(⋅) - функция, которая возвращает 1, если условие в скобках верно, и ноль в противном случае,
L
L
L
-
I(q)  яркость первого изображения в точке p,
R
R
R
I(q) - яркость второго изображения в точке q,
N - количество точек в заданной окрестности;
находят методом динамического программирования среди множества векторов  оптимальный вектор  сдвигов соответственных точек второго изображения, для которого достигается минимум энергетической функции
i
- повторяют шаги определения вектора  для каждого луча S,
- выбирают среди полученных 16-ти векторов-сдвигов  тот, для которого достигается наименьшее значение
- выбирают в векторе  первую координату, отвечающую точке р первого изображения. Значение сдвига, записанное в этой точке, является искомым точным значением сдвига  соответственной точки второго изображения;
L
- повторяют вычисление сдвига  для всех точек р первого изображения и формируют для этого изображения результирующую матрицу сдвигов  в которой на месте яркости I(р) исходного изображения стоит значение
Литература
[1] H. Hirschmuller. Stereo Processing by Semi-Global Matching and Mutual Information // IEEE Transactions on Pattern Analysis and Machine Intelligence, Volume: 30, Issue: 2, Feb. 2008. P. 328-341.
[2] https://racurs.ru/ers-data/optical-images/pleiades/
[3] V. Kolmogorov, R. Zabih, Computing visual correspondence with occlusions using graph cuts. In: Proc. Int'l. Conf. Computer Vision, pp. 5 08-515 (2001) DOI: 10.1109/ICCV.2001.937668
[4] T. Kraub, M. Lehner, P. Reinartz. Generation of coarse 3D models of urban areas from high resolution stereo satellite images // The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences. Vol. XXXVII. Part Bl. Beijing 2008. P. 1091-1098.
[5] Tola E. A DAISY: An Efficient Dense Descriptor Applied to Wide-Baseline Stereo / Tola E., Lepetit V., Fua P. // IEEE Transactions on Pattern Analysis and Machine Intelligence. - 2010 - Vol. 32. Is. 5. P. 815-830.