RU2778076C1 - Method for building a digital surface model based on data from space stereo imaging - Google Patents

Abstract

FIELD: computing technology.

SUBSTANCE: invention relates to the field of computing technology for analysing stereo images when building a digital surface model (DSM) based on data from space stereo imaging. The technical result is achieved by comparing stereo images and finding the shift matrix of the points of the second image relative to the points of the first image in two stages. At the first stage, the shift matrix

is built by finding the minimum of the modified energy function

by the GraphCut (GC) algorithm, and at the second stage, the exact shift matrix

is found by Semi-global matching (SGM) at each point p of the first image by the minimum of the modified energy function

for each i^th ray originating at point p. For each point p of the image of a stereo pair, an attribute vector is formed, called the DAISY(p) descriptor. The coordinates of the DAISY(p) descriptor are two-dimensional convolutions of modules of brightness gradients in 8 directions with Gaussian kernels of varying radii calculated at point p and at 24 more points evenly distributed over three concentric circumferences with a centre at point p.

EFFECT: increase in the accuracy and level of detail when building a digital surface model based on data from space stereo imaging.

1 cl, 6 dwg

Description

Изобретение относится к области картографии и может быть использовано для повышения точности и детальности построения цифровой модели поверхности по данным космической стереосъемки.The invention relates to the field of cartography and can be used to improve the accuracy and detail of building a digital surface model based on space stereo survey data.

Построение цифровой модели поверхности (ЦМП) по данным космической стереосъемки является одной из актуальнейших задач анализа данных дистанционного зондирования Земли, поскольку позволяет извлекать из снимков высотную информацию, представляющую собой ценный источник информации для различных прикладных задач. Однако стереобработка космических изображений сопряжена с рядом фундаментальных трудностей, таких как: высокая изменчивость поля яркости, вызванная атмосферными и иными явлениями; малый размер объектов интереса, сопоставимый с геометрическими искажениями, вызванными разницей углов съемки; сложная траектория движения космического аппарата и специфический способ получения снимков; существенный объем данных для обработки.The construction of a digital surface model (DSM) based on space stereo survey data is one of the most urgent tasks in the analysis of Earth remote sensing data, since it allows extracting altitude information from images, which is a valuable source of information for various applied problems. However, stereo processing of space images is associated with a number of fundamental difficulties, such as: high variability of the brightness field caused by atmospheric and other phenomena; small size of objects of interest, comparable to geometric distortions caused by the difference in shooting angles; complex trajectory of the spacecraft and a specific method of obtaining images; significant amount of data to process.

Точность ЦМП определяется точностью определения высот объектов местности (высот рельефа и относительных высот зданий и сооружений) и точностью положения объектов в плане. Детальность ЦМП определяется наличием относительно небольших по площади (порядка десятков кв. метров) и по высоте (порядка единиц метров) объектов, выступающих на фоне рельефа или на фоне плоскостей более крупных объектов. Точность и детальность ЦМП может быть измерена путем сравнения с эталонной моделью поверхности, построенной по результатам воздушного лазерного сканирования и имеющей более высокую (сантиметровую) точность в плане и по высоте.The accuracy of the DSM is determined by the accuracy of determining the heights of terrain objects (heights of the relief and the relative heights of buildings and structures) and the accuracy of the position of objects in the plan. The detail of the DMF is determined by the presence of relatively small in area (of the order of tens of square meters) and in height (of the order of meters) objects protruding against the background of the relief or against the background of the planes of larger objects. The accuracy and detail of the DTM can be measured by comparison with a reference surface model built from the results of airborne laser scanning, which has a higher (centimeter) accuracy in plan and height.

Важнейшим вопросом при создании способа построения ЦМП является выбор алгоритма стереосопоставления. В подавляющем большинстве современных способов и пакетов программ авторы используют алгоритм [1] Semi-Global Matching (SGM). Данный алгоритм совмещает в себе высокую вычислительную эффективность с не менее высоким качеством построенной карты глубины. Однако, эта эффективность достигается на стереопарах сантиметровой точности, получаемых при аэросъемке, которые обладают более детализированной и предсказуемой структурой по сравнению с космическими стереопарами. В то же время точность и детальность работы данного алгоритма на космических стереопарах, получаемых, например, системой Pleiades [2], имеющих пространственное разрешение порядка 0.5 м, недостаточны. Альтернативой данному методу является метод разрезов графа GraphCut (GC) [3], который благодаря своей глобальной структуре лучше подходит для космических снимков. Однако, он не обладает той же чувствительностью к малым изменениям глубины, то есть к микрорельефу, что и алгоритм SGM.The most important issue in creating a method for constructing a DTM is the choice of a stereo matching algorithm. In the vast majority of modern methods and software packages, the authors use the algorithm [1] Semi-Global Matching (SGM). This algorithm combines high computational efficiency with equally high quality of the constructed depth map. However, this efficiency is achieved on centimeter-precision stereopairs obtained from aerial photography, which have a more detailed and predictable structure compared to space stereopairs. At the same time, the accuracy and detail of the operation of this algorithm on space stereopairs obtained, for example, by the Pleiades system [2], which have a spatial resolution of about 0.5 m, are insufficient. An alternative to this method is the GraphCut (GC) graph cut method [3], which, due to its global structure, is better suited for satellite images. However, it does not have the same sensitivity to small depth changes, i.e. microrelief, as the SGM algorithm.

Описание алгоритмов GC и SGM нахождения сдвигов между соответственными точками эпиполярно выравненных изображений приведено в Приложениях 1 и 2 соответственно.The description of the GC and SGM algorithms for finding shifts between the corresponding points of epipolarly aligned images is given in Appendixes 1 and 2, respectively.

Прототипом предлагаемого способа построения ЦМП является работа [4], в которой предлагается способ построения ЦМП, а также приводятся результаты сравнения качества полученной ЦМП для различных алгоритмов стереосопоставления.The prototype of the proposed method for constructing a DTM is the work [4], which proposes a method for constructing a DTM, and also presents the results of comparing the quality of the obtained DTM for various stereo matching algorithms.

Представленная в прототипе последовательность операций состоит из следующих шагов.Presented in the prototype sequence of operations consists of the following steps.

1. Загрузка входных данных, состоящих из панхроматического снимка с разрешением не менее 1 м и 4-х канального мультиспектрального изображения, включающего красный, зеленый, синий и инфракрасный каналы яркости, с разрешением не менее 2 м, а также RPC коэффициентов математического преобразования координат точек (пикселей) космического изображения на точки местности и обратно.1. Loading input data consisting of a panchromatic image with a resolution of at least 1 m and a 4-channel multispectral image, including red, green, blue and infrared channels of brightness, with a resolution of at least 2 m, as well as RPC coefficients of mathematical transformation of point coordinates (pixels) of the space image to the terrain points and back.

2. Предварительная обработка исходных данных, включающая выбор соответственных точек для коррекции RPC-коэффициентов и повышение пространственного разрешения мультиспектральных изображений до разрешения панхроматического.2. Preliminary processing of the initial data, including the selection of the corresponding points for the correction of RPC coefficients and the increase in the spatial resolution of multispectral images to panchromatic resolution.

3. Эпиполярное выравнивание изображений стереопары, обеспечивающее размещение соответственных точек в одноименных строках выравненных изображений.3. Epipolar alignment of images of a stereopair, which ensures the placement of the corresponding points in the lines of the same name in the aligned images.

4. Нахождение матрицы D сдвигов пар соответственных точек на стереоснимках через сопоставление стереоизображений поверхности (решение задачи стереосопоставления). Матрица сдвигов D соответственных точек состоит из элементов D(p), каждый из которых представляет собой расстояние в пикселях между точкой р первого изображения, и точкой р+D(p) второго изображения стереопары, причем, поскольку предварительно выполняется эпиполярное выравнивание, то смещение возможно только вдоль строки.4. Finding the matrix D of shifts of pairs of corresponding points on stereo images by matching stereo images of the surface (solution of the stereo matching problem). The shift matrix D of the corresponding points consists of elements D(p), each of which represents the distance in pixels between the point p of the first image, and the point p + D(p) of the second image of the stereopair, and since the epipolar alignment is preliminarily performed, the shift is possible only along the line.

5. Построение цифровой модели рельефа и нормализованной ЦМП.5. Construction of a digital elevation model and a normalized DTM.

6. Создание ортотрансформированного изображения.6. Creation of orthorectified image.

7. Классификация областей растительности, а также строений и сооружений7. Classification of areas of vegetation, as well as buildings and structures

8. Извлечение и моделирование отдельных объектов (зданий и сооружений).8. Extraction and modeling of individual objects (buildings and structures).

Недостатком данного способа при формировании ЦМП по космическим снимкам является недостаточно высокая точность и детальность, достигаемая используемыми алгоритмами.The disadvantage of this method in the formation of the DTM on satellite images is not enough high accuracy and detail achieved by the algorithms used.

Целью изобретения является создание способа построения более точной и высоко детальной ЦМП по данным космической стереосъемки.The aim of the invention is to create a method for constructing a more accurate and highly detailed DTM based on space stereo survey data.

Для достижения поставленной цели повышения точности и детальности ЦМП, в отличие от прототипа, предлагается:To achieve the goal of improving the accuracy and detail of the DSM, in contrast to the prototype, it is proposed:

- учитывать априорно известную цифровую модель рельефа (ЦМР) с разрешением не ниже 150 метров на пиксель;- take into account a priori known digital elevation model (DEM) with a resolution of at least 150 meters per pixel;

- выполнять построение ЦМП по стереоснимкам, в которых априорно исключены из рассмотрения области, маскируемые облаками, а также водные поверхности;- to build a DTM based on stereo images, in which areas masked by clouds, as well as water surfaces are a priori excluded from consideration;

- выполнять сегментацию панхроматических изображений и учитывать ее при стереосопоставлении снимков таким образом, чтобы избежать значительных скачков высоты внутри сегментов;- perform segmentation of panchromatic images and take it into account when comparing stereo images in such a way as to avoid significant height jumps within segments;

- находить матрицу

сдвигов пар соответственных точек на стереоснимках (задача стереосопоставления) за счет модификации алгоритмов GC и SGM и их комбинирования.- find matrix

shifts of pairs of corresponding points on stereo images (stereo matching problem) by modifying the GC and SGM algorithms and combining them.

Модификация алгоритма GC состоит в определении для каждой точки первого и второго изображения вектор-признака DAISY(p) (дескриптора DAISY [5]) и замене в выражении для энергетической функций E_GC(D_CC), приведенной в приложении 1, меры близости E_GC__data(D_GC) наThe modification of the GC algorithm consists in determining for each point of the first and second images the feature vector DAISY(p) (DAISY descriptor [5]) and replacing the expression for the energy functions E _GC (D _CC ), given in Appendix 1, with the proximity measure E _{GC _data} ( _DGC ) on

где D_GC - матрица сдвигов соответственных точек на стереоснимках, формируемая алгоритмом GC,where D _GC is the matrix of shifts of the corresponding points on stereo images, generated by the GC algorithm,

||а||_L2 - норма вектора а в пространстве L₂;|| and || _L2 is the norm of the vector a in the space L ₂ ;

p_L - точки первого снимка стереопары,p _L - points of the first image of the stereopair,

p_R(p_L,D_GC) - точки второго снимка стереопары, соответствующие точкам p_L и находящиеся с ними в одной строке.p _R (p _L ,D _GC ) - points of the second image of the stereopair corresponding to the points p _L and being in the same line with them.

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

: первого - меры близости

точек p_R второго изображения к точкам p_L первого изображения на

, и второго слагаемого

, представляющего собой штраф за несовпадение сдвига

в точке p_L и сдвигов в соседних точках q, находящихся в окрестности точки p_L фиксированного радиуса, на модификацию

. Модифицированная энергетическая функция имеет видA modification of the SGM algorithm consists in replacing both terms in the energy function given in Appendix 2

: first - proximity measures

points p _R of the second image to points p _L of the first image on

, and the second term

, which is the penalty for mismatching the shift

at the point p _L and shifts at neighboring points q located in the vicinity of the point p _L of a fixed radius, for modification

. The modified energy function has the form

где

- вектор сдвигов

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

- shift vector

points of the second image p _R relative to the points p _L of the first, located on the ray S _i directed to the point p of the first image, for which the shift is being sought,

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

- a modified measure of the proximity of points p _R of the second image relative to the points p _L of the first image, located on the beam S _i directed to the point p of the first image,

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

is the total penalty modified by introducing penalty coefficients C _S (p _L ,q) for the mismatch between shifts at points p _L located on the ray S _i and shifts at points q adjacent to p _L from a neighborhood of a fixed radius,

C_S(p_L, q) - коэффициент сегментации, который равен 1, если яркости точек p_L и q в сегментированном изображении не равны, и 3 - в противном случае,C _S (p _L , q) - segmentation coefficient, which is equal to 1 if the brightness of the points p _L and q in the segmented image are not equal, and 3 otherwise,

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

- shifts of the points of the second image corresponding to the points p _L and q from the first image located on the ray S _i ,

p_L ∈ S_i - точки первого изображения, расположенные на луче S_i,p _L ∈ S _i - points of the first image located on the ray S _i ,

С₁ и С₂ - штрафы за несовпадение сдвигов в соседних точках,C ₁ and C ₂ - penalties for non-coincidence of shifts at neighboring points,

q ∈ N(p_L)∈S_i - множество точек, расположенных на луче S_i, из окрестности фиксированного радиуса с центром в точке p_L,q ∈ N(p _L )∈S _i is the set of points located on the ray S _i , from a neighborhood of a fixed radius centered at the point p _L ,

Т(⋅) - функция, которая возвращает 1, если условие в скобках верно, и ноль в противном случае.T(⋅) is a function that returns 1 if the condition in parentheses is true, and zero otherwise.

Комбинирование модифицированных алгоритмов заключается в том, что по матрице

которая доставляет минимум функции

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

(где p_L ∈ S_i, записанными в матрице

, по минимумам энергетических функций

для всех лучей S_i находят в точке р первого изображения точную оценку сдвига

соответственной ей точки второго изображения.The combination of the modified algorithms is that, according to the matrix

which delivers the minimum of the function

, find rough estimates of the shifts of the corresponding points of the second image relative to the points of the first. Further, locally at the point p and the rays S _i included in it in the intervals of shifts, which are determined by the shifts

(where p _L ∈ S _i , written in the matrix

, by minima of energy functions

for all rays S _i find at the point p of the first image an exact estimate of the shift

the corresponding point of the second image.

Более подробно сущность предлагаемого способа поясняется дальнейшим описанием и чертежами:In more detail, the essence of the proposed method is explained by the following description and drawings:

на фиг. 1 представлен предлагаемый алгоритм построения ЦМП по данным космической стереосъемки;in fig. 1 shows the proposed algorithm for constructing a DTM based on space stereo survey data;

на фиг. 2 представлен предлагаемый алгоритм стереосопоставления фрагментов изображений стереопары;in fig. 2 shows the proposed algorithm for stereo matching of image fragments of a stereo pair;

на фиг. 3 представлен пример эпиполярно выравненных фрагментов стереоснимков;in fig. 3 shows an example of epipolarly aligned fragments of stereo images;

на фиг. 4 представлена схема расположения ключевых точек для вычисления дескриптора DAISY;in fig. 4 is a diagram of the location of key points for calculating the DAISY descriptor;

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

;in fig. 5 shows the geometric connection of the point p _L on the ray S _i with the points p _R in the search area for the shift of the correspondence point

;

на фиг. 6 представлены для сравнения два ЦМП: первое изображение, полученное предлагаемым способом по данным космической стереосъемки, второе изображение (эталонное ЦМП) по данным лазерной аэросъемки.in fig. Figure 6 shows for comparison two DMFs: the first image obtained by the proposed method from space stereo survey data, the second image (reference DMF) from laser aerial survey data.

Предлагаемый способ построения ЦМП по данным космической стереосъемки в соответствии со схемой фиг. 1 и фиг. 2 производится в последовательности:The proposed method for constructing a DTM based on space stereo data in accordance with the scheme of Fig. 1 and FIG. 2 is produced in the sequence:

1. Получение входных данных (поз. 1 фиг. 1), включающих:1. Obtaining input data (pos. 1 of Fig. 1), including:

- пары стереоснимков местности в панхроматическом диапазоне с исключенными пикселями, соответствующими местоположению облаков и водной поверхности;- pairs of stereo images of the area in the panchromatic range with excluded pixels corresponding to the location of the clouds and the water surface;

- априорно известную ЦМР с разрешением не ниже 150 метров на пиксель;- a priori known DTM with a resolution of at least 150 meters per pixel;

- RPC-коэффициенты математического преобразования координат точек (пикселей) космического изображения на точки местности и обратно;- RPC-coefficients of mathematical transformation of the coordinates of points (pixels) of the space image to the points of the terrain and vice versa;

2. Сегментацию изображений стереопары (поз. 2 фиг. 1) по яркости методом анизотропной медианной фильтрации для удаления шума и мелких выбросов яркости, включающую:2. Image segmentation of a stereopair (pos. 2 of Fig. 1) by brightness by anisotropic median filtering to remove noise and small brightness spikes, including:

- выпуск из каждой точки изображения 16 лучей длиною 7 пикселей,- release from each point of the image 16 rays 7 pixels long,

- вычисление дисперсии яркости вдоль каждого луча,- calculation of the brightness dispersion along each beam,

- вычисление медианы яркости точек вдоль луча с наименьшей дисперсией,- calculation of the median of the brightness of points along the ray with the smallest dispersion,

- присвоение точке панхроматического изображения найденное значение медианы,- assigning a point of the panchromatic image to the found value of the median,

3. Разбиение стереоснимков на одинаковые по размеру меньшие фрагменты (поз. 3 фиг. 1);3. Splitting stereo images into smaller fragments of the same size (pos. 3 of Fig. 1);

4. Эпиполярное выравнивание фрагментов изображений стереопары в панхроматическом диапазоне через повороты изображений относительно центров снимков так, чтобы смещения, вызванные сменой ракурса при получении стереопары, лежали вдоль одной строки (поз. 4 фиг 1). Пример эпиполярного выравнивания приведен на фиг. 3, горизонтальные линии вдоль снимков соответствуют строкам выровненных изображений;Fig. 4. Epipolar alignment of image fragments of a stereopair in the panchromatic range through rotation of images relative to the centers of images so that the displacements caused by a change in angle when obtaining a stereopair lie along one line (pos. 4 of Fig. 1). An example of an epipolar alignment is shown in FIG. 3, the horizontal lines along the pictures correspond to the lines of the aligned pictures;

5. Выполнение процедуры стереосопоставления (поз. 5 фиг. 1), в результате которой получают матрицу сдвигов соответственных точек. Процедура стереосопоставления включает в себя следующие этапы.5. Execution of the stereo matching procedure (pos. 5 of Fig. 1), as a result of which a matrix of shifts of the corresponding points is obtained. The stereo matching procedure includes the following steps.

5.1 Вычисление дескриптора DAISY (р) для каждой точки р изображений стереопары (поз. 5.1 фиг. 2) по алгоритму:5.1 Calculation of the DAISY descriptor (p) for each point p of images of a stereopair (pos. 5.1 of Fig. 2) according to the algorithm:

- вокруг точки р панхроматического изображения строят набор из трех концентрических окружностей радиусами 3, 6 и 9 пикселей;- a set of three concentric circles with radii of 3, 6 and 9 pixels is built around the point p of the panchromatic image;

- на каждой из этих окружностей выбирают по 8 точек, равномерно распределенных вдоль окружности, геометрия размещения точек приведена на фиг. 4;- on each of these circles, 8 points are selected, evenly distributed along the circle, the geometry of the placement of points is shown in Fig. four;

- в исходной точке р, а также в 8 точках окружности радиусом 3 пикселя выполняют свертку с гауссовым ядром модулей градиентов яркости по 8 направлениям, радиус гауссового ядра 2 пикселя;- at the starting point p, as well as at 8 points of a circle with a radius of 3 pixels, convolution with the Gaussian kernel of the modules of brightness gradients in 8 directions is performed, the radius of the Gaussian kernel is 2 pixels;

- в каждой из 8 точек окружности радиусом 6 пикселей выполняют свертку с гауссовым ядром модулей градиентов яркости по 8 направлениям, радиус гауссового ядра 3 пикселя;- in each of the 8 points of the circle with a radius of 6 pixels, convolution is performed with the Gaussian kernel of the brightness gradient modules in 8 directions, the radius of the Gaussian kernel is 3 pixels;

- в каждой из 8 точек окружности радиусом 9 пикселей выполняют свертку с гауссовым ядром модулей градиентов яркости по 8 направлениям, радиус гауссового ядра 4 пикселя;- in each of the 8 points of the circle with a radius of 9 pixels perform convolution with the Gaussian kernel modules of brightness gradients in 8 directions, the radius of the Gaussian kernel is 4 pixels;

- объединяют полученные значения в вектор размерности 200 (25 точек по 8 направлений градиента в каждой);- combine the obtained values into a vector of dimension 200 (25 points with 8 gradient directions in each);

- нормализуют полученный вектор, который считают вектором-признаком (дескриптором) DAISY (р) точки р изображения.- normalize the resulting vector, which is considered a feature vector (descriptor) DAISY (p) of the point p of the image.

5.2 Выполнение первого этапа стереосопоставления с помощью модифицированного алгоритма GC (поз. 5.2 фиг. 2) в последовательности:5.2 Performing the first stage of stereo matching using the modified GC algorithm (pos. 5.2 of Fig. 2) in the sequence:

- определяют диапазон поиска

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

shifts of correspondence points in the digital relief matrix (DEM) as an interval

где k_p - коэффициент перевода 1 м высоты местности в пиксели смещения на стереопаре, вычисляемый индивидуально для каждого фрагмента по RPC-коэффициентам снимков,where k _p is the coefficient of conversion of 1 m of terrain height into offset pixels on a stereopair, calculated individually for each fragment according to the RPC coefficients of images,

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

- the minimum and maximum height of the area under consideration according to the DEM within the fragment;

- находят методом оптимальную матрицу сдвигов

такую, что значение

каждого элемента матрицы принадлежит интервалу

, то есть

, и для которой достигается минимум модифицированной энергетической функции

:- find the optimal shift matrix by the method

such that the value

each element of the matrix belongs to the interval

, that is

, and for which the minimum of the modified energy function is reached

:

где

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

- measure of proximity of points p _L of the first image to points p _R of the second image,

p_L=(x_L,y_L) - точки первого снимка стереопары,p _L =(x _L ,y _L ) - points of the first image of the stereopair,

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 _R (p _L, D _GC )=(x _R ,y _R )=(x _L +D _GC (p _L ),y _L ) - point of the second image of the stereopair, the shift of which relative to the point p _L of the first image is equal to D _GC ( p _L ) and is within the interval

,

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

,D _GC (p _L ) - shift of the point of the second image relative to the point p _L of the first image, located in the same image lines and within the interval

,

N(p_L)- множество точек в окрестности точки p_L фиксированного радиуса,N(p _L ) is the set of points in the vicinity of the point p _L of a fixed radius,

V - весовой коэффициент штрафной функции,V is the weighting coefficient of the penalty function,

Т(⋅) - функция, которая возвращает 1, если условие в скобках верно, и ноль в противном случае,T(⋅) - a function that returns 1 if the condition in brackets is true, and zero otherwise,

С - весовой коэффициент штрафной функции для условия D_GC(p_L)=∅ (∅ - пустое множество) означающего, что сдвиг точки p_R относительно точки p_L не определен.C is the weight coefficient of the penalty function for the condition D _GC (p _L )=∅ (∅ is an empty set) meaning that the shift of the point p _R relative to the point p _L is not defined.

5.3 Выполнение второго этапа стереосопоставления (точного сопоставления фрагментов разбиения) с помощью модифицированного алгоритма SGM (поз. 5.3 фиг. 2) в последовательности:5.3 Performing the second stage of stereo matching (exact matching of partition fragments) using the modified SGM algorithm (pos. 5.3 of Fig. 2) in the sequence:

определяют диапазон поиска

сдвигов для каждой точки снимка p_L как интервалdefine the search range

shifts for each point of the image p _L as an interval

в точку р первого изображения от края изображения проводят 16 равномерно разнесенных по углу лучей (i=1…16). В дальнейшем тексте точки первого изображения, находящиеся на S_i луче, обозначены p_L. Расположение точек р и p_L на луче S_i показано на фиг. 5;to the point p of the first image from the edge of the image, 16 beams evenly spaced along the angle (i=1...16) are carried out. In the following text, the points of the first image located on the S _i beam are denoted p _L . The location of the points p and p _L on the ray S _i is shown in FIG. 5;

для i-того луча находят вектор

сдвигов

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

соответственной точки p_R иллюстрирует фиг. 5) следующим образом:for the i-th ray find the vector

shifts

of the corresponding points of the second image, located in the same-named lines with p _L , (connection of points p _L with the search area

corresponding point p _R is illustrated in FIG. 5) as follows:

- рассматривают множество вариантов

вектора сдвигов

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

(из заранее заданного диапазона)- consider multiple options

shift vectors

points p _Rs of the second image with respect to points p _L ∈ S _i of the first image, whose coordinates are successive values

(out of predefined range)

- для s-того варианта вектора сдвигов

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

:- for the s-th version of the shift vector

determine the value of the energy function

:

где

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

where

- a measure of the proximity of all points p _L of the first image, located on the beam S _i , to the points p _Rs of the second image, having shifts

C_S(p_L,q) - коэффициент сегментации, который равен 1, если яркости точек p_L и q в сегментированном изображении не равны, и 3 - в противном случае,C _S (p _L ,q) - segmentation coefficient, which is equal to 1 if the brightness of the points p _L and q in the segmented image are not equal, and 3 otherwise,

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

- modified due to the introduction of coefficients C _S (p _L ,q) the total penalty for the mismatch between shifts at points p _L located on the ray S _i and shifts at points q ∈ S _i adjacent to p _L from a neighborhood of a fixed radius,

p_L ∈ S_i - точки первого изображения, расположенные на луче S_i,p _L ∈ S _i - points of the first image located on the ray S _i ,

С₁ и С₂ - штрафы за несовпадение сдвигов в соседних точках,C ₁ and C ₂ - penalties for non-coincidence of shifts at neighboring points,

q ∈ N(p_L) - множество точек в окрестности фиксированного радиуса с центром в точке p_L луча S_i,q ∈ N(p _L ) is the set of points in a neighborhood of a fixed radius centered at the point p _L of the ray S _i ,

Т(⋅) - функция, которая возвращает 1, если условие в скобках верно, и ноль в противном случае;T(⋅) - a function that returns 1 if the condition in brackets is true, and zero otherwise;

- находят методом динамического программирования среди множества векторов

) оптимальный вектор

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

- find by dynamic programming among a set of vectors

) optimal vector

shifts of the corresponding points of the second image, for which the minimum of the energy function is achieved

повторяют шаги определения вектора

для каждого луча S_i;repeat the steps of defining the vector

for each ray S _i ;

выбирают среди полученных 16-ти векторов-сдвигов

тот, для которого достигается наименьшее значение

choose among the obtained 16 shift vectors

the one for which the smallest value is reached

выбирают в векторе

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

соответственной точки второго изображения;choose in vector

the first coordinate corresponding to the point p of the first image, the shift value recorded at this point is the desired exact shift value

corresponding point of the second image;

повторяют вычисление сдвига

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

, в которой на месте яркости h(р) исходного изображения стоит значение

repeat the shift calculation

for all points p of the first image and form the resulting shift matrix for this image

, in which the brightness h(p) of the original image is replaced by the value

6. Построение ЦМП фрагментов изображения в географических координатах (поз. 6 фиг. 1) в последовательности:6. Construction of the DTM of image fragments in geographic coordinates (pos. 6 of Fig. 1) in the sequence:

вычисляют высоту местности в пикселе р по найденным сдвигам

, базовой высоте

фрагмента и коэффициенту k_p перевода метров высоты в пиксели изображения:calculate the height of the terrain in pixel p from the found shifts

, base height

fragment and the coefficient k _p of converting meters of height into pixels of the image:

по матрице h(p) вычисляют географические координаты точек поверхности Н(Х, Y) с помощью RPC-коэффициентов для того, чтобы получить ЦМП, привязанную к местности.the matrix h(p) calculates the geographical coordinates of the points of the surface H(X, Y) using RPC coefficients in order to obtain a DTM tied to the area.

7. Объединение ЦМП всех фрагментов разбиения Н(Х, Y) в итоговую ЦМП (поз. 7 фиг. 1) и выполнение медианной фильтрации полученной матрицы высот.7. Combining the DTM of all fragments of the partition H(X, Y) into the final DTM (pos. 7 of Fig. 1) and performing median filtering of the resulting DEM.

Для проверки работоспособности способа был проведен вычислительный эксперимент по построению цифровой модели местности по данным космической стереосъемки. Было выполнено сравнение качества ЦМП, построенной предлагаемым способом, с ЦМП прототипа, использующего сопоставление стереоснимков по алгоритму GC. Сравнение проводилось по значениям среднеквадратического отклонения оценок высоты ЦМП от эталонной ЦМП, построенной по данным воздушной лазерной съемки. Сравнение производилось по областям, не отмеченным в ходе предварительной классификации как водная поверхность и облачный или растительный покров. Полученное в прототипе среднеквадратичное значение ошибки (СКО) 3,74 м было улучшено до 2,93 м., что составляет уменьшение СКО ошибки на 21,6%. При этом для полученной ЦМП также были оценены точность по высоте (среднее абсолютное отклонение высот объектов ЦМП от высот объектов на эталоне), которая составила 1,81 м и точность в плане (отклонение границ объектов ЦМП - строений от границ этих же объектов на эталоне), которая составила 1,87 м.To test the operability of the method, a computational experiment was carried out to build a digital terrain model based on space stereo imaging data. The quality of the DTM built by the proposed method was compared with the DTM of the prototype, using the comparison of stereo images using the GC algorithm. The comparison was carried out using the values of the root-mean-square deviation of the DMF height estimates from the reference DMF constructed from airborne laser survey data. Comparisons were made for areas not identified during the preliminary classification as water surface and cloud or vegetation cover. Received in the prototype root-mean-square error (RMS) of 3.74 m was improved to 2.93 m, which is a decrease in the error RMS by 21.6%. At the same time, for the obtained DSM, the accuracy in height (the average absolute deviation of the heights of the DSM objects from the heights of objects on the standard) was also estimated, which amounted to 1.81 m and the accuracy in plan (deviation of the boundaries of the objects of the DSM - buildings from the boundaries of the same objects on the standard) , which amounted to 1.87 m.

На фиг. 6 показаны фрагменты матрицы высот рельефа, вычисленных предложенным методом, и соответствующие им результаты лидарной съемки. Визуальное сравнение показывает, что вычисленные матрицы высот рельефа качественно близки к результатам лидарной съемки.In FIG. Figure 6 shows fragments of the elevation matrix calculated by the proposed method and the corresponding results of the lidar survey. A visual comparison shows that the calculated elevation matrixes are qualitatively close to the results of the lidar survey.

Представленный способ построения ЦМП может найти свое применение для решения широкого круга задач, которые удобно разбить на три направления.The presented method of constructing a DMP can be used to solve a wide range of problems, which can be conveniently divided into three areas.

Первым важным направлением является формирование эталонных ЦМП по космическим стереоснимкам, необходимых для навигации беспилотных летательных аппаратов. Вне зависимости от типа задачи, данный способ обладает рядом важных преимуществ в сравнении с построением ЦМП по аэросъемке: более широкое покрытие; возможность получения данных на территории, где использование беспилотных летательных аппаратов невозможно; меньшая стоимость; высокая оперативность получения требуемой информации.The first important direction is the formation of reference DMF based on space stereo images necessary for the navigation of unmanned aerial vehicles. Regardless of the type of task, this method has a number of important advantages in comparison with the construction of a DTM based on aerial photography: wider coverage; the possibility of obtaining data on the territory where the use of unmanned aerial vehicles is impossible; lower cost; high efficiency in obtaining the required information.

Второе направление связано с управлением территориями, картографией, контролем темпов строительства и обнаружения объектов незаконного строительства, корректировкой кадастровой и градостроительной документации, оценкой геоморфологических изменений.The second direction is related to the management of territories, cartography, monitoring the pace of construction and detection of illegal construction sites, updating cadastral and urban planning documentation, and assessing geomorphological changes.

Третье направление - природный и экологический мониторинг, инвентаризация сельскохозяйственных угодий, оценка динамики изменения площади лесных массивов, создание тематических природных карт.The third direction is natural and environmental monitoring, inventory of agricultural lands, assessment of the dynamics of changes in the area of forests, creation of thematic natural maps.

Приложение 1Attachment 1

Алгоритм GC поиска матрицы сдвигов между соответственными точками эпиполярно выравненных стереоизображенийGC search algorithm for shift matrix between corresponding points of epipolar aligned stereo images

Задачу поиска матрицы сдвигов D между точками первого p_L и второго p_R эпиполярно выровненных стереоизображений решают через нахождение минимума «энергетической» функции E_GC(⋅). Значения элементов матриц D для алгоритмов GC и SGM определяются по-разному, поэтому в дальнейшем тексте обозначаются D_GC для алгоритма GC и D_SGM для алгоритма SGM.The task of finding the shift matrix D between the points of the first p _L and second p _R epipolarly aligned stereo images is solved by finding the minimum of the "energy" function E _GC (⋅). The values of the elements of the matrices D for the GC and SGM algorithms are defined differently, therefore, in the following text, they are denoted by D _GC for the GC algorithm and D _SGM for the SGM algorithm.

Энергетическая функция E_GC(D_GC) имеет вид:The energy function E _GC (D _GC ) has the form:

E_GC(D_GC)=E_GC__data(D_GC)+E_GC__smth(D_GC),E _GC (D _GC )=E _GC _ _data (D _GC )+E _GC _ _smth (D _GC ),

где

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

- measure of proximity of points p _L of the first image to points p _R of the second image,

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

- total penalty for non-coincidence of shifts at points p _L and shifts at points q adjacent to p _L from a neighborhood of a fixed radius,

p_L=(.x_L,y_L) - каждая точка первого снимка стереопары,p _L =(.x _L ,y _L ) - each point of the first image of the stereopair,

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 _R (p _L ,D _GC )=(x _R ,y _R )=(x _L +D _GC (p _L ),y _L ) - point of the second image of the stereopair, the shift of which relative to the point p _L of the first image is equal to D _GC ( p _L ),

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

- average brightness of the first and second images of the stereopair in the vicinity of the points p _L and p _R , respectively,

N(p_L) - множество точек в окрестности точки p_L фиксированного радиуса,N(p _L ) - set of points in the vicinity of the point p _L of fixed radius,

V - весовой коэффициент штрафной функции,V is the weighting coefficient of the penalty function,

Т(⋅) - функция, которая возвращает 1, если условие в скобках верно, и ноль в противном случае,T(⋅) - a function that returns 1 if the condition in brackets is true, and zero otherwise,

С - весовой коэффициент штрафной функции для условия D_GC(p_L)=∅ (∅ - пустое множество) означающего, что сдвиг точки p_R относительно точки p_L не определен.C is the weight coefficient of the penalty function for the condition D _GC (p _L )=∅ (∅ is an empty set) meaning that the shift of the point p _R relative to the point p _L is not defined.

Минимум функции E_GC(D_GC) в алгоритме GC ищут для всех точек изображения одновременно специальным численным методом, который относится к классу глобальных алгоритмов.The minimum of the function E _GC (D _GC ) in the GC algorithm is searched for for all points of the image simultaneously by a special numerical method, which belongs to the class of global algorithms.

Приложение 2Appendix 2

Алгоритм SGM поиска матрицы сдвигов между точками эпиполярно выравненных изображенийSGM search algorithm for shift matrix between points of epipolar aligned images

В алгоритме SGM поиск минимума функции E_SGM(⋅), в отличие от алгоритма GC, осуществляется не глобально для всех точек сразу, а для каждой точки р первого изображения отдельно и независимо от других, поэтому алгоритм относится к классу полу-глобальных. Для этогоIn the SGM algorithm, the search for the minimum of the function E _SGM (⋅), in contrast to the GC algorithm, is carried out not globally for all points at once, but for each point p of the first image separately and independently of the others, so the algorithm belongs to the class of semi-global ones. For this

- в точку р первого изображения от края изображения проводят 16 равномерно разнесенных по углу лучей S_i, в дальнейшем тексте точки первого изображения, находящиеся на луче S_i, обозначаем p_L=(x_L,y_L),- to the point p of the first image from the edge of the image, 16 evenly spaced rays S _i are carried out, in the further text the points of the first image located on the beam S _i are denoted by p _L =(x _L ,y _L ),

- для i-того луча находят вектор

сдвигов

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

shifts

corresponding points of the second image, located in the same-named lines with p _L , as follows:

рассматривают множество вариантов

вектора сдвигов

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

(из заранее заданного диапазона)considering many options

shift vectors

points p _Rs of the second image with respect to points p _L ∈ S _i of the first image, whose coordinates are successive values

(out of predefined range)

для s-того варианта вектора сдвигов

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

for the s-th version of the shift vector

determine the value of the energy function

где

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

where

- a measure of the proximity of all points p _L of the first image, located on the beam S _i , to the points p _Rs of the second image, having shifts

- суммарный штраф за несовпадение сдвигов

в точках p_L, расположенных на луче S_i, и сдвигов

в соседних с p_L точках q ∈ S_i из окрестности фиксированного радиуса,

- total penalty for mismatched shifts

at points p _L located on the ray S _i , and shifts

at points q ∈ S _i adjacent to p _L from a neighborhood of a fixed radius,

p_L ∈ S_i - точки первого изображения, расположенные на луче S_i,p _L ∈ S _i - points of the first image located on the ray S _i ,

С₁ и С₂ - штрафы за несовпадение сдвигов в соседних точках,C ₁ and C ₂ - penalties for non-coincidence of shifts at neighboring points,

q_L ∈ N(p_L) - множество точек первого изображения в окрестности точки p_L фиксированного радиуса,q _L ∈ N(p _L ) is the set of points of the first image in the vicinity of the point p _L of a fixed radius,

q_R ∈ N(p_Rs) - множество точек второго изображения в окрестности точки p_Rs фиксированного радиуса,q _R ∈ N(p _Rs ) is the set of points of the second image in the vicinity of the point p _Rs of a fixed radius,

Т(⋅) - функция, которая возвращает 1, если условие в скобках верно, и ноль в противном случае,T(⋅) - a function that returns 1 if the condition in brackets is true, and zero otherwise,

I_L(q_L) ^- яркость первого изображения в точке p_L,I _L (q _L ) ^- brightness of the first image at point p _L ,

I_R(q_R) - яркость второго изображения в точке q_R,I _R (q _R ) - brightness of the second image at point q _R ,

N - количество точек в заданной окрестности;N is the number of points in a given neighborhood;

находят методом динамического программирования среди множества векторов

оптимальный вектор

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

are found by dynamic programming among a set of vectors

optimal vector

shifts of the corresponding points of the second image, for which the minimum of the energy function is achieved

- повторяют шаги определения вектора

для каждого луча S_i,- repeat the steps of defining the vector

for each ray S _i ,

- выбирают среди полученных 16-ти векторов-сдвигов

тот, для которого достигается наименьшее значение

- choose among the received 16 shift vectors

the one for which the smallest value is reached

- выбирают в векторе

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

соответственной точки второго изображения;- choose in vector

the first coordinate corresponding to the p point of the first image. The shift value written at this point is the exact shift value you are looking for

corresponding point of the second image;

- повторяют вычисление сдвига

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

в которой на месте яркости I_L(р) исходного изображения стоит значение

- repeat the calculation of the shift

for all points p of the first image and form the resulting shift matrix for this image

in which the brightness I _L (p) of the original image is replaced by the value

ЛитератураLiterature

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

Claims (39)

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

для которой достигается минимум энергетической функции E_GC(D_GC)A method for constructing a digital surface model (DSM) based on space stereo survey data includes obtaining a pair of stereo images of the area in the panchromatic range and RPC coefficients for mathematical transformation of the coordinates of points (pixels) of the space image to the points of the area and vice versa; splitting stereo images into smaller fragments of the same size; epipolar alignment of image fragments of a stereopair in the panchromatic range through rotation of images relative to the centers of image fragments so that the displacements caused by a change in angle when obtaining a stereopair lie along one line; performing a stereo matching procedure, which results in a matrix of shifts of the corresponding points; stereo matching procedure includes: finding the optimal shift matrix

for which the minimum of the energy function E _GC (D _GC )

E_GC(D_GC)=E_{GC_data}(D_GC)+E_GC__smth(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 второго изображения,where E _{GC_data} (D _GC ) is a measure of proximity of points p _L of the first image to points p _R of the second image, E_GC__smth(D_GC) - суммарный штраф за несовпадение сдвигов в точках p_L и сдвигов в соседних с p_L точках q из окрестности фиксированного радиуса;E _GC _ _smth (D _GC ) - total penalty for non-coincidence of shifts at p _L points and shifts at q points adjacent to p _L from a neighborhood of a fixed radius; значение E_{GC_smth} (D_GC) определяется выражениемthe value of E _{GC_smth} (D _GC ) is given by

где p_L=(x_L, y_L) - каждая точка первого снимка стереопары,where p _L =(x _L, y _L ) - each point of the first image of the stereopair, p_R(p_L,D_DC)=(x_R, y_R)=(x_L+D_GC(p_L),y_L) - соответственная точка второго снимка стереопары, сдвиг которой относительно точки p_L первого снимка равен D_GC(p_L),p _R (p _L, D _DC )=(x _R , y _R )=(x _L +D _GC (p _L ),y _L ) - the corresponding point of the second image of the stereopair, the shift of which relative to the point p _L of the first image is equal to D _GC (p _L ), D_GC(p_L) - сдвиг соответственной точки второго изображения относительно точки p_L первого изображения, находящихся в одноименных строках изображений,D _GC (p _L ) - shift of the corresponding point of the second image relative to the point p _L of the first image, located in the same image lines, N(p_L) - множество точек в окрестности точки p_L фиксированного радиуса,N(p _L ) - set of points in the vicinity of the point p _L of fixed radius, V - весовой коэффициент штрафной функции,V is the weighting coefficient of the penalty function, Т(⋅) - функция, которая возвращает 1, если условие в скобках верно, и ноль - в противном случае,T(⋅) - a function that returns 1 if the condition in brackets is true, and zero otherwise, С - весовой коэффициент штрафной функции для условия D_GC(p_L)=∅ (∅ - пустое множество), означающего, что сдвиг точки p_R относительно точки p_L не определен;C is the weight coefficient of the penalty function for the condition D _GC (p _L )=∅ (∅ is an empty set), which means that the shift of the point p _R relative to the point p _L is not defined; строят ЦМП фрагментов изображения в географических координатах в последовательности:DTM of image fragments is built in geographic coordinates in the sequence: вычисляют высоту местности 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(р) точки р изображения; диапазон поиска

сдвигов соответственных точек по ЦМР вычисляется как интервалcalculating the terrain height h(p) in pixel p based on the found shifts and shooting parameters; calculate the geographic coordinates of the points of the surface H(X,Y) from the matrix h(p) using RPC coefficients in order to obtain a DTM tied to the area; combine the DTM of all fragments of the partition H(X,Y) into the final DTM and perform median filtering of the resulting DEM; characterized in that, as a result of preliminary processing, pixels corresponding to the location of clouds and the water surface are excluded from panchromatic stereo images taken to build the DSM; a priori known digital elevation model (DEM) with a resolution of at least 150 meters per pixel is included in the initial information for building a DEM; before splitting the stereo images, segmentation of images of a stereo pair is performed using anisotropic median filtering to remove noise and small brightness spikes, image segmentation includes: issuing 16 rays of 7 pixels long from each point p of the image, calculating the dispersion of the image brightness along each beam, calculating the median of the brightness of points along the beam with the smallest variance, assigning the point p of the panchromatic image of the found value of the median; in the stereo matching operation before finding the optimal shift matrix

for each point p of images of a stereopair, a feature vector of the point, called the DAISY(p) descriptor, and the search range are calculated

shifts of the corresponding points along the DEM; the descriptor DAISY (p) is calculated according to the algorithm: a set of three concentric circles with radii of 3, 6 and 9 pixels is built around the point p of the panchromatic image; on each of these circles, 8 points are selected, evenly distributed along the circle; at the starting point p, as well as at 8 points selected on a circle with a radius of 3 pixels, perform a two-dimensional convolution with a Gaussian kernel modules of brightness gradients in 8 directions, the radius of the Gaussian kernel is 2 pixels; at each of the 8 points selected on a circle with a radius of 6 pixels, convolution is performed with the Gaussian kernel of the brightness gradient modules in 8 directions, the radius of the Gaussian kernel is 3 pixels; in each of the 8 points selected on a circle with a radius of 9 pixels, perform convolution with the Gaussian kernel modules of brightness gradients in 8 directions, the radius of the Gaussian kernel is 4 pixels; combine the obtained values into a vector of dimension 200; normalizing the resulting vector, which is considered the descriptor DAISY(p) of the point p of the image; search range

shifts of the corresponding points along the DTM is calculated as an interval

где k_p - коэффициент перевода 1 м высоты местности в пиксели смещения на стереопаре, вычисляемый для каждого фрагмента по RPC-коэффициентам снимков,where k _p is the coefficient of conversion of 1 m of terrain height into offset pixels on a stereopair, calculated for each fragment by the RPC coefficients of images,

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

- the minimum and maximum height of the area under consideration according to the DEM within the fragment; вычисления энергетической функции E_GC(D_GC) выполняют для пар точек первого p_L и второго р_R изображения, сдвиг между которыми D_GC(p_L) находится в пределах интервала

после нахождения оптимальной матрицы сдвигов

ищется матрица точных значений

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

сдвигов для каждой точки снимка p_L как интервалcalculations of the energy function E _GC (D _GC ) are performed for pairs of points of the first p _L and second p _R images, the shift between which D _GC (p _L ) is within the interval

after finding the optimal shift matrix

looking for a matrix of exact values

in the local areas of each point p of the first image by the modified SGM algorithm in the sequence: determine the search range

shifts for each point of the image p _L as an interval

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

сдвигов

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

вектора сдвигов

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

to the point p of the first image from the edge of the image spend 16 evenly spaced along the angle of the rays S _i (i=1...16); in the following text, according to the SGM algorithm, the points of the first image located on the ray S _i are denoted p _L ; for the i-th ray find the vector

shifts

of the corresponding points of the second image, located in the lines of the same name with p _L , for this, a set of options is considered

shift vectors

points p _Rs of the second image with respect to points p _L ∈ S _i of the first image, whose coordinates are successive values

для s-того варианта вектора сдвигов

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

по формуле:for the s-th version of the shift vector

determine the value of the energy function

according to the formula:

гдеwhere

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

- modified measure of proximity of all points p _L of the first image, located on the beam S _i , to points p _Rs of the second image, having shifts

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

- modified due to the introduction of coefficients C _S (p _L , q) the total penalty for the mismatch of shifts at points p _L located on the ray S _i and shifts at points q ∈ S _i adjacent to p _L from a neighborhood of a fixed radius, C_S(p_L,q) - коэффициент сегментации, который равен 1, если яркости точек p_L и q в сегментированном изображении не равны, и 3 - в противном случае,C _S (p _L ,q) - segmentation coefficient, which is equal to 1 if the brightness of the points p _L and q in the segmented image are not equal, and 3 otherwise, p_L ∈ S_i - точки первого изображения, расположенные на луче S_i,p _L ∈ S _i - points of the first image located on the ray S _i , С₁ и С₂ - штрафы за несовпадение сдвигов в соседних точках,C ₁ and C ₂ - penalties for non-coincidence of shifts at neighboring points, q ∈ N(p_L) ∩ S_i - множество точек луча S_i в окрестности фиксированного радиуса с центром в точке p_L,q ∈ N(p _L ) ∩ S _i is the set of points of the ray S _i in a neighborhood of a fixed radius centered at the point p _L , находят методом динамического программирования среди множества векторов

оптимальный вектор

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

are found by dynamic programming among a set of vectors

optimal vector

shifts of the corresponding points of the second image, for which the minimum of the energy function is achieved

повторяют шаги определения вектора

для каждого луча S_i;repeat the steps of defining the vector

for each ray S _i ; выбирают среди полученных 16-ти векторов-сдвигов

тот, для которого достигается наименьшее значение

choose among the obtained 16 shift vectors

the one for which the smallest value is reached

выбирают в векторе

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

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

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

, в которой на месте яркости I_L(р) исходного изображения стоит значение

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

базовой высоте

фрагмента и коэффициенту k_p перевода метров высоты в пиксели изображения:choose in vector

the first coordinate corresponding to the point p of the first image, the shift value recorded at this point is the desired exact shift value

corresponding point of the second image; repeat the shift calculation

for all points p of the first image and form the resulting shift matrix for this image

, in which the brightness I _L (p) of the original image is replaced by the value

; in the operations of constructing a DTM of image fragments in geographic coordinates, the value h(p) - the height of the terrain in pixel p is determined by the found shifts

base height

fragment and the coefficient k _p of converting meters of height into pixels of the image:

Images