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.

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.

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.

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. 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. 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. 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. 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. Construction of a digital elevation model and a normalized DTM.

6. Creation of orthorectified image.

7. Classification of areas of vegetation, as well as buildings and structures

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:

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

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

where D _GC is the matrix of shifts of the corresponding points on stereo images, generated by the GC algorithm,

|| and || _L2 is the norm of the vector a in the space L ₂ ;

p _L - points of the first image of the stereopair,

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.

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

точек 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) - 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 - 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 is the set of points located on the ray S _i , from a neighborhood of a fixed radius centered at the point p _L ,

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:

in fig. 1 shows the proposed algorithm for constructing a DTM based on space stereo survey data;

in fig. 2 shows the proposed algorithm for stereo matching of image fragments of a stereo pair;

in fig. 3 shows an example of epipolarly aligned fragments of stereo images;

in fig. 4 is a diagram of the location of key points for calculating the DAISY descriptor;

;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

;

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.

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

- a priori known DTM with a resolution of at least 150 meters per pixel;

- 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. Image segmentation of a stereopair (pos. 2 of Fig. 1) by brightness by anisotropic median filtering to remove noise and small brightness spikes, including:

- 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. Splitting stereo images into smaller fragments of the same size (pos. 3 of Fig. 1);

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

- 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, the geometry of the placement of points is shown in Fig. four;

- 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;

- 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;

- 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;

- combine the obtained values into a vector of dimension 200 (25 points with 8 gradient directions in each);

- normalize the resulting vector, which is considered a feature vector (descriptor) DAISY (p) of the point p of the image.

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

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 ) - 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 ) - 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 ) - 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 ) is the set of points in the vicinity of the point p _L of a fixed radius,

V is the weighting coefficient of the penalty function,

T(⋅) - a function that returns 1 if the condition in brackets is true, and zero otherwise,

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

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;

сдвигов

соответственных точек второго изображения, находящихся в одноименных с 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)

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

:- 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) - 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 - 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 ) is the set of points in a neighborhood of a fixed radius centered at the point p _L of the ray S _i ,

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 ;

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

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

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

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.

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.

Attachment 1

GC search algorithm for shift matrix between corresponding points of epipolar aligned stereo images

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.

The energy function E _GC (D _GC ) has the form:

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 ) - 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 ) - 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 ) - set of points in the vicinity of the point p _L of fixed radius,

V is the weighting coefficient of the penalty function,

T(⋅) - a function that returns 1 if the condition in brackets is true, and zero otherwise,

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.

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.

Appendix 2

SGM search algorithm for shift matrix between points of epipolar aligned images

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

- 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 ),

сдвигов

соответственных точек второго изображения, находящихся в одноименных с 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)

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

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 - 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 ) 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 ) is the set of points of the second image in the vicinity of the point p _Rs of a fixed radius,

T(⋅) - a function that returns 1 if the condition in brackets is true, and zero otherwise,

I _L (q _L ) ^- brightness of the first image at point p _L ,

I _R (q _R ) - brightness of the second image at point q _R ,

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 ,

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

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

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

Claims (39)

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 ), 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 ) - 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; the value of E _{GC_smth} (D _GC ) is given by

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 ) - 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 ) - 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 ) - set of points in the vicinity of the point p _L of fixed radius, V is the weighting coefficient of the penalty function, T(⋅) - a function that returns 1 if the condition in brackets is true, and zero otherwise, 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: 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

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; 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

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

for the s-th version of the shift vector

determine the value of the energy function

according to the formula:

where

- 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

- 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) - 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 - 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 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

repeat the steps of defining the vector

for each ray S _i ; 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; 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