Способ построения цифровой модели поверхности по данным космической стереосъемки

Изобретение относится к области вычислительной техники для анализа стереоизображений при построении цифровой модели поверхности (ЦМП) по данным космической стереосъемки. Технический результат заключается в повышении точности и детальности построения цифровой модели поверхности по данным космической стереосъемки. Технический результат достигается за счет сопоставления стереоизображений с нахождением матрицы сдвигов точек второго изображения относительно точек первого изображения в два этапа. На первом этапе построение матрицы сдвигов выполняется через нахождение минимума модифицированной энергетической функции алгоритмом GraphCut (GC), а на втором нахождение методом Semi-global matching (SGM) точной матрицы сдвигов в каждой точке р первого изображения по минимуму модифицированной энергетической функции для каждого i-го луча, исходящего из точки р. Для каждой точки р изображения стереопары формируется вектор-признак, называемый дескриптором DAISY(р). Координаты дескриптора DAISY(p) являются двумерными свертками модулей градиентов яркости по 8 направлениям с гауссовыми ядрами различных радиусов, вычисленных в точке р и еще в 24-х точках, равномерно распределенных по трем концентрическим окружностям с центром в точке р. 6 ил.

 

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

Построение цифровой модели поверхности (ЦМП) по данным космической стереосъемки является одной из актуальнейших задач анализа данных дистанционного зондирования Земли, поскольку позволяет извлекать из снимков высотную информацию, представляющую собой ценный источник информации для различных прикладных задач. Однако стереобработка космических изображений сопряжена с рядом фундаментальных трудностей, таких как: высокая изменчивость поля яркости, вызванная атмосферными и иными явлениями; малый размер объектов интереса, сопоставимый с геометрическими искажениями, вызванными разницей углов съемки; сложная траектория движения космического аппарата и специфический способ получения снимков; существенный объем данных для обработки.

Точность ЦМП определяется точностью определения высот объектов местности (высот рельефа и относительных высот зданий и сооружений) и точностью положения объектов в плане. Детальность ЦМП определяется наличием относительно небольших по площади (порядка десятков кв. метров) и по высоте (порядка единиц метров) объектов, выступающих на фоне рельефа или на фоне плоскостей более крупных объектов. Точность и детальность ЦМП может быть измерена путем сравнения с эталонной моделью поверхности, построенной по результатам воздушного лазерного сканирования и имеющей более высокую (сантиметровую) точность в плане и по высоте.

Важнейшим вопросом при создании способа построения ЦМП является выбор алгоритма стереосопоставления. В подавляющем большинстве современных способов и пакетов программ авторы используют алгоритм [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 состоит в определении для каждой точки первого и второго изображения вектор-признака DAISY(p) (дескриптора DAISY [5]) и замене в выражении для энергетической функций E_GC(D_CC), приведенной в приложении 1, меры близости E_GC__data(D_GC) на

где D_GC - матрица сдвигов соответственных точек на стереоснимках, формируемая алгоритмом GC,

||а||_L2 - норма вектора а в пространстве L₂;

p_L - точки первого снимка стереопары,

p_R(p_L,D_GC) - точки второго снимка стереопары, соответствующие точкам p_L и находящиеся с ними в одной строке.

Модификация алгоритма SGM состоит в замене обоих слагаемых в приведенной в приложении 2 энергетической функции : первого - меры близости точек p_R второго изображения к точкам p_L первого изображения на , и второго слагаемого , представляющего собой штраф за несовпадение сдвига в точке p_L и сдвигов в соседних точках q, находящихся в окрестности точки p_L фиксированного радиуса, на модификацию . Модифицированная энергетическая функция имеет вид

где - вектор сдвигов точек второго изображения p_R относительно точек p_L первого, находящихся на луче S_i, направленном в точку р первого изображения, для которой ищется сдвиг,

- модифицированная мера близости точек p_R второго изображения относительно точек p_L первого изображения, находящихся на луче S_i, направленном в точку р первого изображения,

- суммарный штраф, модифицированный за счет введения штрафных коэффициентов C_S(p_L,q) за несовпадение сдвигов в точках p_L, расположенных на луче S_i, и сдвигов в соседних с p_L точках q из окрестности фиксированного радиуса,

C_S(p_L, q) - коэффициент сегментации, который равен 1, если яркости точек p_L и q в сегментированном изображении не равны, и 3 - в противном случае,

- сдвиги точек второго изображения, соответствующих точкам p_L и q из первого изображения, находящихся на луче S_i,

p_L ∈ S_i - точки первого изображения, расположенные на луче S_i,

С₁ и С₂ - штрафы за несовпадение сдвигов в соседних точках,

q ∈ N(p_L)∈S_i - множество точек, расположенных на луче S_i, из окрестности фиксированного радиуса с центром в точке p_L,

Т(⋅) - функция, которая возвращает 1, если условие в скобках верно, и ноль в противном случае.

Комбинирование модифицированных алгоритмов заключается в том, что по матрице которая доставляет минимум функции , находят грубые оценки сдвигов соответственных точек второго изображения относительно точек первого. Далее локально в точке р и входящих в нее лучей S_i в интервалах сдвигов, которые определяются сдвигами (где p_L ∈ S_i, записанными в матрице , по минимумам энергетических функций для всех лучей S_i находят в точке р первого изображения точную оценку сдвига соответственной ей точки второго изображения.

Более подробно сущность предлагаемого способа поясняется дальнейшим описанием и чертежами:

на фиг. 1 представлен предлагаемый алгоритм построения ЦМП по данным космической стереосъемки;

на фиг. 2 представлен предлагаемый алгоритм стереосопоставления фрагментов изображений стереопары;

на фиг. 3 представлен пример эпиполярно выравненных фрагментов стереоснимков;

на фиг. 4 представлена схема расположения ключевых точек для вычисления дескриптора DAISY;

на фиг. 5 представлена геометрическая связь точки p_L на луче S_i с точками p_R в зоне поиска сдвига точки соответствия ;

на фиг. 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) в последовательности:

- определяют диапазон поиска сдвигов точек соответствия по цифровой матрице рельефа (ЦМР) как интервал

где k_p - коэффициент перевода 1 м высоты местности в пиксели смещения на стереопаре, вычисляемый индивидуально для каждого фрагмента по RPC-коэффициентам снимков,

- минимальная и максимальная высота рассматриваемой местности по ЦМР в пределах фрагмента;

- находят методом оптимальную матрицу сдвигов такую, что значение каждого элемента матрицы принадлежит интервалу , то есть , и для которой достигается минимум модифицированной энергетической функции :

где - мера близости точек p_L первого изображения точкам p_R второго изображения,

p_L=(x_L,y_L) - точки первого снимка стереопары,

p_R(p_L,D_GC)=(x_R,y_R)=(x_L+D_GC(p_L),y_L) - точка второго снимка стереопары, сдвиг которой относительно точки p_L первого снимка равен D_GC(p_L) и находится в пределах интервала ,

D_GC(p_L) - сдвиг точки второго изображения относительно точки p_L первого изображения, находящихся в одноименных строках изображений и в пределах интервала ,

N(p_L)- множество точек в окрестности точки p_L фиксированного радиуса,

V - весовой коэффициент штрафной функции,

Т(⋅) - функция, которая возвращает 1, если условие в скобках верно, и ноль в противном случае,

С - весовой коэффициент штрафной функции для условия D_GC(p_L)=∅ (∅ - пустое множество) означающего, что сдвиг точки p_R относительно точки p_L не определен.

5.3 Выполнение второго этапа стереосопоставления (точного сопоставления фрагментов разбиения) с помощью модифицированного алгоритма SGM (поз. 5.3 фиг. 2) в последовательности:

определяют диапазон поиска сдвигов для каждой точки снимка p_L как интервал

в точку р первого изображения от края изображения проводят 16 равномерно разнесенных по углу лучей (i=1…16). В дальнейшем тексте точки первого изображения, находящиеся на S_i луче, обозначены p_L. Расположение точек р и p_L на луче S_i показано на фиг. 5;

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

- рассматривают множество вариантов вектора сдвигов точек p_Rs второго изображения относительно точек p_L ∈ S_i первого изображения, координаты которого - последовательные значения (из заранее заданного диапазона)

- для s-того варианта вектора сдвигов определяют значение энергетической функции :

где - мера близости всех точек p_L первого изображения, расположенных на луче S_i, к точкам p_Rs второго изображения, имеющих сдвиги

C_S(p_L,q) - коэффициент сегментации, который равен 1, если яркости точек p_L и q в сегментированном изображении не равны, и 3 - в противном случае,

- модифицированный за счет введения коэффициентов C_S(p_L,q) суммарный штраф за несовпадение сдвигов в точках p_L, расположенных на луче S_i, и сдвигов в соседних с p_L точках q ∈ S_i из окрестности фиксированного радиуса,

p_L ∈ S_i - точки первого изображения, расположенные на луче S_i,

С₁ и С₂ - штрафы за несовпадение сдвигов в соседних точках,

q ∈ N(p_L) - множество точек в окрестности фиксированного радиуса с центром в точке p_L луча S_i,

Т(⋅) - функция, которая возвращает 1, если условие в скобках верно, и ноль в противном случае;

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

повторяют шаги определения вектора для каждого луча S_i;

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

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

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

6. Построение ЦМП фрагментов изображения в географических координатах (поз. 6 фиг. 1) в последовательности:

вычисляют высоту местности в пикселе р по найденным сдвигам , базовой высоте фрагмента и коэффициенту k_p перевода метров высоты в пиксели изображения:

по матрице h(p) вычисляют географические координаты точек поверхности Н(Х, Y) с помощью RPC-коэффициентов для того, чтобы получить ЦМП, привязанную к местности.

7. Объединение ЦМП всех фрагментов разбиения Н(Х, Y) в итоговую ЦМП (поз. 7 фиг. 1) и выполнение медианной фильтрации полученной матрицы высот.

Для проверки работоспособности способа был проведен вычислительный эксперимент по построению цифровой модели местности по данным космической стереосъемки. Было выполнено сравнение качества ЦМП, построенной предлагаемым способом, с ЦМП прототипа, использующего сопоставление стереоснимков по алгоритму GC. Сравнение проводилось по значениям среднеквадратического отклонения оценок высоты ЦМП от эталонной ЦМП, построенной по данным воздушной лазерной съемки. Сравнение производилось по областям, не отмеченным в ходе предварительной классификации как водная поверхность и облачный или растительный покров. Полученное в прототипе среднеквадратичное значение ошибки (СКО) 3,74 м было улучшено до 2,93 м., что составляет уменьшение СКО ошибки на 21,6%. При этом для полученной ЦМП также были оценены точность по высоте (среднее абсолютное отклонение высот объектов ЦМП от высот объектов на эталоне), которая составила 1,81 м и точность в плане (отклонение границ объектов ЦМП - строений от границ этих же объектов на эталоне), которая составила 1,87 м.

На фиг. 6 показаны фрагменты матрицы высот рельефа, вычисленных предложенным методом, и соответствующие им результаты лидарной съемки. Визуальное сравнение показывает, что вычисленные матрицы высот рельефа качественно близки к результатам лидарной съемки.

Представленный способ построения ЦМП может найти свое применение для решения широкого круга задач, которые удобно разбить на три направления.

Первым важным направлением является формирование эталонных ЦМП по космическим стереоснимкам, необходимых для навигации беспилотных летательных аппаратов. Вне зависимости от типа задачи, данный способ обладает рядом важных преимуществ в сравнении с построением ЦМП по аэросъемке: более широкое покрытие; возможность получения данных на территории, где использование беспилотных летательных аппаратов невозможно; меньшая стоимость; высокая оперативность получения требуемой информации.

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

Третье направление - природный и экологический мониторинг, инвентаризация сельскохозяйственных угодий, оценка динамики изменения площади лесных массивов, создание тематических природных карт.

Приложение 1

Алгоритм GC поиска матрицы сдвигов между соответственными точками эпиполярно выравненных стереоизображений

Задачу поиска матрицы сдвигов D между точками первого p_L и второго p_R эпиполярно выровненных стереоизображений решают через нахождение минимума «энергетической» функции E_GC(⋅). Значения элементов матриц D для алгоритмов GC и SGM определяются по-разному, поэтому в дальнейшем тексте обозначаются D_GC для алгоритма GC и D_SGM для алгоритма SGM.

Энергетическая функция E_GC(D_GC) имеет вид:

E_GC(D_GC)=E_GC__data(D_GC)+E_GC__smth(D_GC),

где - мера близости точек p_L первого изображения точкам p_R второго изображения,

- суммарный штраф за несовпадение сдвигов в точках p_L и сдвигов в соседних с p_L точках q из окрестности фиксированного радиуса,

p_L=(.x_L,y_L) - каждая точка первого снимка стереопары,

p_R(p_L,D_GC)=(x_R,y_R)=(x_L+D_GC(p_L),y_L) - точка второго снимка стереопары, сдвиг которой относительно точки p_L первого снимка равен D_GC(p_L),

- средние яркости первого и второго изображений стереопары в окрестностях точек p_L и p_R соответственно,

N(p_L) - множество точек в окрестности точки p_L фиксированного радиуса,

V - весовой коэффициент штрафной функции,

Т(⋅) - функция, которая возвращает 1, если условие в скобках верно, и ноль в противном случае,

С - весовой коэффициент штрафной функции для условия D_GC(p_L)=∅ (∅ - пустое множество) означающего, что сдвиг точки p_R относительно точки p_L не определен.

Минимум функции E_GC(D_GC) в алгоритме GC ищут для всех точек изображения одновременно специальным численным методом, который относится к классу глобальных алгоритмов.

Приложение 2

Алгоритм SGM поиска матрицы сдвигов между точками эпиполярно выравненных изображений

В алгоритме SGM поиск минимума функции E_SGM(⋅), в отличие от алгоритма GC, осуществляется не глобально для всех точек сразу, а для каждой точки р первого изображения отдельно и независимо от других, поэтому алгоритм относится к классу полу-глобальных. Для этого

- в точку р первого изображения от края изображения проводят 16 равномерно разнесенных по углу лучей S_i, в дальнейшем тексте точки первого изображения, находящиеся на луче S_i, обозначаем p_L=(x_L,y_L),

- для i-того луча находят вектор сдвигов соответственных точек второго изображения, находящихся в одноименных с p_L строках, следующим образом:

рассматривают множество вариантов вектора сдвигов точек p_Rs второго изображения относительно точек p_L ∈ S_i первого изображения, координаты которого - последовательные значения (из заранее заданного диапазона)

для s-того варианта вектора сдвигов определяют значение энергетической функции

где - мера близости всех точек p_L первого изображения, расположенных на луче S_i, к точкам p_Rs второго изображения, имеющих сдвиги

- суммарный штраф за несовпадение сдвигов в точках p_L, расположенных на луче S_i, и сдвигов в соседних с p_L точках q ∈ S_i из окрестности фиксированного радиуса,

p_L ∈ S_i - точки первого изображения, расположенные на луче S_i,

С₁ и С₂ - штрафы за несовпадение сдвигов в соседних точках,

q_L ∈ N(p_L) - множество точек первого изображения в окрестности точки p_L фиксированного радиуса,

q_R ∈ N(p_Rs) - множество точек второго изображения в окрестности точки p_Rs фиксированного радиуса,

Т(⋅) - функция, которая возвращает 1, если условие в скобках верно, и ноль в противном случае,

I_L(q_L) ^- яркость первого изображения в точке p_L,

I_R(q_R) - яркость второго изображения в точке q_R,

N - количество точек в заданной окрестности;

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

- повторяют шаги определения вектора для каждого луча S_i,

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

- выбирают в векторе первую координату, отвечающую точке р первого изображения. Значение сдвига, записанное в этой точке, является искомым точным значением сдвига соответственной точки второго изображения;

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

Литература

[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.

Способ построения цифровой модели поверхности (ЦМП) по данным космической стереосъемки включает получение пары стереоснимков местности в панхроматическом диапазоне и RPC-коэффициентов математического преобразования координат точек (пикселей) космического изображения на точки местности и обратно; разбиение стереоснимков на одинаковые по размеру меньшие фрагменты; эпиполярное выравнивание фрагментов изображений стереопары в панхроматическом диапазоне через повороты изображений относительно центров фрагментов снимков так, чтобы смещения, вызванные сменой ракурса при получении стереопары, лежали вдоль одной строки; выполнение процедуры стереосопоставления, в результате которой получают матрицу сдвигов соответственных точек; процедура стереосопоставления включает в себя: нахождение оптимальной матрицы сдвигов для которой достигается минимум энергетической функции E_GC(D_GC)

E_GC(D_GC)=E_{GC_data}(D_GC)+E_GC__smth(D_GC),

где E_{GC_data}(D_GC) - мера близости точек p_L первого изображения к точкам p_R второго изображения,

E_GC__smth(D_GC) - суммарный штраф за несовпадение сдвигов в точках p_L и сдвигов в соседних с p_L точках q из окрестности фиксированного радиуса;

значение E_{GC_smth} (D_GC) определяется выражением

где p_L=(x_L, y_L) - каждая точка первого снимка стереопары,

p_R(p_L,D_DC)=(x_R, y_R)=(x_L+D_GC(p_L),y_L) - соответственная точка второго снимка стереопары, сдвиг которой относительно точки p_L первого снимка равен D_GC(p_L),

D_GC(p_L) - сдвиг соответственной точки второго изображения относительно точки p_L первого изображения, находящихся в одноименных строках изображений,

N(p_L) - множество точек в окрестности точки p_L фиксированного радиуса,

V - весовой коэффициент штрафной функции,

Т(⋅) - функция, которая возвращает 1, если условие в скобках верно, и ноль - в противном случае,

С - весовой коэффициент штрафной функции для условия D_GC(p_L)=∅ (∅ - пустое множество), означающего, что сдвиг точки p_R относительно точки p_L не определен;

строят ЦМП фрагментов изображения в географических координатах в последовательности:

вычисляют высоту местности h(p) в пикселе р по найденным сдвигам и параметрам съемки; вычисляют по матрице h(p) географические координаты точек поверхности H(X,Y) с помощью RPC-коэффициентов для того, чтобы получить ЦМП, привязанную к местности; объединяют ЦМП всех фрагментов разбиения H(X,Y) в итоговую ЦМП и выполняют медианную фильтрацию полученной матрицы высот; отличающийся тем, что из панхроматических стереоснимков, принимаемых для построения ЦМП, в результате предварительной обработки исключены пиксели, соответствующие местоположению облаков и водной поверхности; в состав исходной информации для построения ЦМП включена априорно известная цифровая модель рельефа (ЦМР) с разрешением не ниже 150 метров на пиксель; перед разбиением стереоснимков выполняют сегментацию изображений стереопары методом анизотропной медианной фильтрации для удаления шума и мелких выбросов яркости, сегментация изображений включает: выпуск из каждой точки р изображения 16 лучей длиною 7 пикселей, вычисление дисперсии яркости изображения вдоль каждого луча, вычисление медианы яркости точек вдоль луча с наименьшей дисперсией, присвоение точке р панхроматического изображения найденного значения медианы; в операции стереосопоставления перед нахождением оптимальной матрицы сдвигов для каждой точки р изображений стереопары вычисляют вектор-признак точки, называемый дескриптором DAISY(p), и диапазон поиска сдвигов соответственных точек по ЦМР; дескриптор DAISY (р) вычисляют по алгоритму: вокруг точки р панхроматического изображения строят набор из трех концентрических окружностей радиусами 3, 6 и 9 пикселей; на каждой из этих окружностей выбирают по 8 точек, равномерно распределенных вдоль окружности; в исходной точке р, а также в 8 точках, выбранных на окружности с радиусом 3 пикселя, выполняют двумерную свертку с гауссовым ядром модулей градиентов яркости по 8 направлениям, радиус гауссова ядра 2 пикселя; в каждой из 8 точек, выбранных на окружности с радиусом 6 пикселей, выполняют свертку с гауссовым ядром модулей градиентов яркости по 8 направлениям, радиус гауссова ядра 3 пикселя; в каждой из 8 точек, выбранных на окружности с радиусом 9 пикселей, выполняют свертку с гауссовым ядром модулей градиентов яркости по 8 направлениям, радиус гауссова ядра 4 пикселя; объединяют полученные значения в вектор размерности 200; нормализуют полученный вектор, который считают дескриптором DAISY(р) точки р изображения; диапазон поиска сдвигов соответственных точек по ЦМР вычисляется как интервал

где k_p - коэффициент перевода 1 м высоты местности в пиксели смещения на стереопаре, вычисляемый для каждого фрагмента по RPC-коэффициентам снимков,

- минимальная и максимальная высота рассматриваемой местности по ЦМР в пределах фрагмента;

вычисления энергетической функции E_GC(D_GC) выполняют для пар точек первого p_L и второго р_R изображения, сдвиг между которыми D_GC(p_L) находится в пределах интервала после нахождения оптимальной матрицы сдвигов ищется матрица точных значений в локальных областях каждой точки р первого снимка модифицированным алгоритмом SGM в последовательности: определяют диапазон поиска сдвигов для каждой точки снимка p_L как интервал

в точку р первого изображения от края изображения проводят 16 равномерно разнесенных по углу лучей S_i (i=1…16); в дальнейшем тексте по алгоритму SGM точки первого изображения, находящиеся на луче S_i, обозначены p_L; для i-того луча находят вектор сдвигов соответственных точек второго изображения, находящихся в одноименных с p_L строках, для этого рассматривают множество вариантов вектора сдвигов точек p_Rs второго изображения относительно точек p_L ∈ S_i первого изображения, координаты которого - последовательные значения

для s-того варианта вектора сдвигов определяют значение энергетической функции по формуле:

где

- модифицированная мера близости всех точек p_L первого изображения, расположенных на луче S_i, к точкам p_Rs второго изображения, имеющим сдвиги

- модифицированный за счет введения коэффициентов C_S(p_L, q) суммарный штраф за несовпадение сдвигов в точках p_L, расположенных на луче S_i, и сдвигов в соседних с p_L точках q ∈ S_i из окрестности фиксированного радиуса,

C_S(p_L,q) - коэффициент сегментации, который равен 1, если яркости точек p_L и q в сегментированном изображении не равны, и 3 - в противном случае,

p_L ∈ S_i - точки первого изображения, расположенные на луче S_i,

С₁ и С₂ - штрафы за несовпадение сдвигов в соседних точках,

q ∈ N(p_L) ∩ S_i - множество точек луча S_i в окрестности фиксированного радиуса с центром в точке p_L,

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

повторяют шаги определения вектора для каждого луча S_i;

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

выбирают в векторе первую координату, отвечающую точке р первого изображения, значение сдвига, записанное в этой точке, является искомым точным значением сдвига соответственной точки второго изображения; повторяют вычисление сдвига для всех точек р первого изображения и формируют для этого изображения результирующую матрицу сдвигов , в которой на месте яркости I_L(р) исходного изображения стоит значение ; в операциях построения ЦМП фрагментов изображения в географических координатах значение h(р) - высоты местности в пикселе р определяют по найденным сдвигам базовой высоте фрагмента и коэффициенту k_p перевода метров высоты в пиксели изображения: