CN116894923A - High-resolution remote sensing image mapping conversion dense matching and three-dimensional reconstruction method - Google Patents

Abstract

The invention discloses a high-resolution remote sensing image mapping conversion dense matching and three-dimensional reconstruction method, which includes the steps of: receiving a first remote sensing image; performing orthorectification on the first remote sensing image to obtain a second remote sensing image; and performing orthorectification on the second remote sensing image. The remote sensing image is resampled to obtain the third remote sensing image; the third remote sensing image is adaptively expanded into blocks, and the block images are densely matched to generate the corresponding disparity map; based on the disparity map and block information, each Restore the image point with the same name to its coordinates on the orthorectified image, then perform inverse orthorectification on the point pair with the same name, restore the coordinates of the point pair to the first high-resolution remote sensing image, and obtain the corresponding three-dimensional point coordinates formed by Three-dimensional point cloud; raster sampling of three-dimensional point cloud to generate high-resolution DSM.

Description

A high-resolution remote sensing image mapping conversion dense matching and three-dimensional reconstruction method

Technical field

The invention belongs to the field of image processing technology, and particularly relates to a high-resolution remote sensing image mapping conversion dense matching and three-dimensional reconstruction method.

Background technique

Laser detection, radar ranging (LiDAR) data and aerial remote sensing images have been widely used to generate large-scale digital surface models (DSM, Digital Surface Model). Large-scale digital surface models (DSM) generated from high-resolution satellite imagery (HRSI) have the advantage of being comparable, cheaper, and more readily available.

In recent years, with the rapid development of earth observation technology, the earth observation system has become increasingly perfect. Foreign high-resolution earth observation remote sensing satellite technologies represented by IKONOS, WorldView, and Pleiades have become very mature, and high-resolution civilian mapping satellites represented by Ziyuan-3 have also reached the world's advanced level. The high-resolution linear array satellites listed above can form stereoscopic observations (for example, the Ziyu-3 linear array sensor can provide front-view, downward-view and rear-view three-view stereo), thus providing an effective and reliable data source for three-dimensional information extraction. . Currently, the applications of 3D information extraction from high-resolution linear array satellite stereoscopic observation data include but are not limited to: digital surface model generation, urban 3D modeling, natural disaster assessment, and target positioning and tracking.

Contents of the invention

One embodiment of the present invention is a three-dimensional reconstruction method of high-resolution remote sensing images, which uses high-resolution remote sensing image mapping conversion dense matching and three-dimensional reconstruction. Includes the following steps:

Receive the first remote sensing image;

Perform orthorectification on the first remote sensing image to obtain a second remote sensing image;

Resample the second remote sensing image to obtain the third remote sensing image;

Adaptively expand the third remote sensing image into blocks, and perform dense matching on the effects of the blocks to generate corresponding disparity maps;

Based on the disparity map and block information, each pair of image points with the same name is restored to its coordinates on the orthorectified image, and then inverse orthorectification is performed on the same-named point pair to restore the coordinates of the point pair to the first high-resolution remote sensing image. above, obtain the three-dimensional point cloud formed by the corresponding three-dimensional point coordinates;

Raster sampling of 3D point clouds generates high-resolution DSM.

Description of the drawings

The above and other objects, features and advantages of exemplary embodiments of the present invention will become apparent upon reading the following detailed description with reference to the accompanying drawings. In the drawings, several embodiments of the invention are shown by way of example and not by way of limitation, in which:

Figure 1 is a schematic diagram of an existing three-dimensional reconstruction processing flow chart of satellite images.

Figure 2 is a schematic flowchart of a three-dimensional reconstruction method for remote sensing images according to one embodiment of the present invention.

Figure 3 is an example diagram of the row and column difference of the same image point after orthorectification according to one embodiment of the present invention.

Figure 4 is an example diagram of a cost accumulation process in disparity space according to one embodiment of the present invention.

Figure 5 is an example of a remote sensing image according to one embodiment of the present invention.

Figure 6 is an example of a remote sensing image according to one embodiment of the present invention.

Figure 7 is an example of a remote sensing image according to one embodiment of the present invention.

Figure 8 is an example of a remote sensing image according to one embodiment of the present invention.

Figure 9 is an example diagram of a remote sensing image comparison according to one embodiment of the present invention. The picture on the left is the DSM generated by the embodiment of the present invention, and the picture on the right is the STRM 30m resolution DEM.

Figure 10 is an example diagram of three-dimensional reconstruction of a remote sensing image according to one embodiment of the present invention.

Detailed ways

Most existing satellite image 3D reconstruction software uses an improved version of the semi-global matching (SGM) algorithm to perform dense image matching based on epipolar line constraints. A core technology for extracting three-dimensional information from satellite images is stereoscopic dense matching technology. As the input data of dense matching, the correctness of the epipolar image directly affects the reliability of the matching results, and thus the accuracy of the three-dimensional information. Therefore, correct epipolar line image input is crucial for dense matching.

Different from the frame image, the epicenter of the linear array image is hyperbolic, which brings great challenges to the generation of epipolar images. The traditional epipolar image generation method is suitable for a single stereoscopic observation mode (such as same orbit or different orbit) and lacks versatility. In addition, with the increasing number of Earth observation satellites, the possibility of stereoscopic observation between different platforms is increasing. Therefore, an alternative core line that can be applied to co-orbital, out-of-orbit and heterosource stereo pairs is obtained. Image generation methods became a need.

A stereo matching algorithm has been proposed. Its four steps mainly include: matching cost calculation, cost aggregation, disparity calculation and disparity refinement. From the perspective of matching strategy, stereo matching can be divided into: local stereo matching algorithm, semi-global stereo matching algorithm and global stereo matching algorithm.

The local stereo matching algorithm has strong real-time performance, but for pixels with discontinuous parallax at the edge, it is easy to cause mismatching, resulting in low matching accuracy, and the algorithm is easily affected by noise. The global algorithm has high matching accuracy, but the space and time complexity of this type of algorithm is high, and it has been very difficult to put it into practical engineering applications until now. The semi-global algorithm has good matching speed and accuracy and can meet the real-time requirements, so it is widely used in stereo matching. This includes the proposed semi-global matching algorithm (Semi-global Matching, SGM).

At present, various software for 3D reconstruction of satellite images mostly follow the following processing flow:

First, feature matching is performed on satellite stereo image pairs to obtain image points with the same name, in order to achieve geometric correction of the satellite camera model and obtain accurate satellite image orientation parameters;

Then, in order to solve the problem that the satellite image frame is too large and the epipolar line is an approximate curve, the satellite image is processed into blocks:

First, epipolar line correction is performed on each pair of image blocks, and then dense matching is performed to obtain a disparity map. Then, based on the disparity map, forward intersection is performed to obtain the point cloud of each block;

Finally, each point cloud is fused to obtain the three-dimensional point cloud of the entire area and the rasterized DSM product. As shown in Figure 1.

It can be seen that in the existing solution, the original satellite stereo image pairs need to be corrected for epipolar lines before dense matching can be carried out. However, due to the large frame and wide coverage of satellite images, excessive terrain differences within the images, and differences in various sensor types, etc. The problem is that there is no universal kernel line constraint method that can be applied to most scenarios and cannot meet the needs when processing large-scale satellite data.

When performing dense matching based on epipolar line corrected images, the traditional tSGM algorithm only searches for identical points on pixels on the same epipolar line. It uses a method of constructing a multi-level pyramid to match from top to bottom, and uses the matching measurement information of the upper layer to Corso searches the lower pyramid's search range, but this method is not applicable to situations where the epipolar lines are not on a straight line. For epipolar lines that are not aligned, the corrected image matching accuracy is very poor.

Therefore, the purpose of this disclosure is to solve the problem that the epipolar line constraint method mentioned above is difficult to apply to the three-dimensional reconstruction of satellite images in complex scenes, and the problem of poor matching accuracy caused by the inability to perform window search in subsequent dense matching.

According to one or more embodiments, a three-dimensional reconstruction method based on orthorectification and window disparity search is provided. The following is a brief introduction to the process.

The first step is to perform orthorectification on the original image based on the reference low-resolution DSM. This step converts the image points to the object coordinate system in units of longitude and latitude. Subsequently, resampling is performed on the basis of the orthorectified image. A resampled image whose coordinates of the same image point are constrained within the rpc accuracy range.

In the second step, the large-width satellite image is adaptively expanded into blocks based on the reference DEM elevation value, and the corresponding disparity map is generated by dense matching for the block image pairs. The specific dense matching process uses a two-dimensional window based on tSGM. The search strategy adopts Gaussian pyramid and uses the same-name point search results in the upper layer to narrow the search range in the lower layer.

The third step is to restore each pair of same-named image points to their coordinates on the orthorectified image based on the disparity map and block information, and then perform an inverse orthorectification process on the same-named point pair based on the reference low resolution, and convert the point pair coordinates Restore to the original image, and then use RPC to perform forward intersection on the original image pair to obtain the corresponding three-dimensional point coordinates.

Finally, the 3D point cloud is subjected to post-processing operations such as denoising, fusion, and water body removal based on the reference DEM. Then the 3D point cloud is rasterized and sampled to generate a high-resolution DSM. Finally, the DSM is subjected to hole filling, bilateral filtering, and multiplexing. DSM fusion operation improves DSM accuracy. The overall scheme of the embodiment of the present disclosure is shown in Figure 2. In Figure 2, for a satellite sensor with only two stereoscopic images, the left and right images are equivalent to the rear-view image and the front-view image respectively. In actual reconstruction, only three-dimensional reconstruction of the two images is required. For satellite sensors with three stereoscopic images, two methods are generally used for reconstruction:

1. Select any two images containing three-view images for reconstruction. The left and right images here are equal to the two selected images.

2. Use the downward view, the forward view, and the downward view and the backward view respectively to reconstruct, and then perform DSM fusion after obtaining two DSMs.

Furthermore, for orthorectified image generation and segmentation, in the current processing flow, in order to improve the accuracy of image matching and reduce the time required for matching, the image needs to be epipolarly corrected before stereoscopic image matching, and the two-dimensional matching is The search is reduced to a one-dimensional search. Normally, the epipolar line of high-resolution remote sensing images is similar to a hyperbola, but it can be approximately regarded as a straight line in a local range. But on the one hand, in areas with large terrain undulations, even if it is small enough Blocking also cannot constrain image points with the same name to the same core line. On the other hand, the currently best method of extracting feature point pairs based on image cubes for basic matrix estimation has one of the biggest problems, which is that it is difficult to ensure feature matching. Obtaining the same-name points is accurate enough that even when using RANSAC to filter the same-name points, it is difficult to remove all wrong point pairs. This makes it difficult to generate epipolar line images that meet the requirements for all terrains.

In view of this, this disclosure uses RPC-based and reference low-resolution DEM to orthorectify the left and right images respectively. The orthorectified image can ensure that the image point with the same name is within the neighborhood, which can also limit the search range. function, and avoids the problem that the failure of the kernel line constraint makes the subsequent dense matching algorithm invalid. When the RPC accuracy is high enough, the search range can be limited to dozens of pixels, and the same epipolar line constraint also needs to be searched in dozens or even higher one-dimensional spaces in areas with large terrain undulations. The distance between the same image points after orthorectification is shown in Figure 3.

As shown in Figure 3, in the Loess Plateau area with terrain undulations exceeding 1000m, the search range of the image point with the same name after orthorectification is still within 5 pixels in the row and row, which can fully ensure that the entire search range of the image point with the same name is within 100 pixels. Within.

After obtaining the orthorectified images corresponding to the left and right images, resampling is needed to unify the resolution and frame size of the orthorectified images. The orthorectified images after resampling have the same geographical coordinate affine parameters. This is It provides great convenience for finding overlapping areas in subsequent left and right images. There is no need to find the corresponding overlapping area based on rpc. In view of the error of rpc, the overlapping area obtained in this way often needs to be edge expanded. In the resampling process, the required DSM resolution is added as an adjustable parameter. According to the adjustable parameters, an orthoimage of any resolution size can be generated, and based on this, a DSM of the same resolution is generated in the subsequent process.

Then, masks will be generated for the orthorectified resampled image pairs respectively. The purpose of generating the masks is to reduce the search for invalid areas in the subsequent dense matching process, which will greatly reduce the matching time.

Finally, the generated orthorectified resampled image and its mask are divided into blocks with extended areas to ensure that there are no gaps between the disparity maps generated by adjacent orthorectified image blocks and that a certain overlap is guaranteed between blocks. .By dividing the image into blocks, on the one hand, the parallax search can be limited to a smaller interval, which can reduce the mismatching caused by the excessive parallax search range in repeated texture areas. On the other hand, it facilitates multi-threading processing, making the processing Each pair of image blocks can be used as an independent thread without interfering with each other, which will greatly reduce the time required to generate the entire DSM.

Before performing dense matching, you can make the following three data preparations:

A. RPC adjustment (this process is not necessary). Regional network adjustment based on feature matching can improve the elevation accuracy obtained by front-end interaction.

B. Reference dem of the corresponding area (required). Dem-based elevation constraints can improve the accuracy of three-dimensional reconstruction.

C. Light and color uniformity (this process is not necessary). For image pairs with inconsistent colors, color uniformity can improve the accuracy of dense matching.

The above is possible data preprocessing, but image pairs without unnecessary preprocessing can also be densely matched, but the accuracy may be worse in some cases (such as the color difference between image pairs is large).

The following describes the tSGM algorithm based on window search in the embodiment of the present disclosure.

1) The basic algorithm of SGM is to establish a global energy function to obtain the potential optimal parallax when the energy function reaches a minimum. The energy function is

The global energy function E(D) consists of data items and smoothing items. The first item on the right side of equation (4.1) is the data item, which represents the sum of the matching costs of all pixels, and C(p,D _p ) represents the left image pixel p The matching cost with the right image point when the disparity is D _p [24]. The second term is the smoothing term, whose purpose is to increase the penalty number P ₁ for the pixel q with a disparity difference equal to 1 near the image point p. This makes the adjacent pixels The disparity of the points is as consistent as possible, where T[|D _p -D _q |=1] is a judgment function, which returns 1 when the judgment condition is true, otherwise it is 0. The third item is the smoothing item, for the pixel point p Pixels q with adjacent parallax differences greater than 1 are given a larger penalty number P ₂ . Usually P ₂ changes with the change of the image gradient, and P ₂ cannot be too large, otherwise the parallax edge will be over-smoothed.

The SGM algorithm first calculates the matching cost of each pixel on the left image within the disparity search range to form a cost matrix. The dimension of the cost matrix is r×c×(d _max -d _min ), where r is the number of image rows and c is the image column. numbers, d _max and d _min represent the upper and lower bounds of disparity search. After calculating the cost matrix, consider accumulating each point on the cost matrix from multiple directions (8 or 16 directions), as shown in Figure 4.

For a pixel point p, its accumulation function in the r direction is

In the formula, the first term on the right side represents the matching cost when the disparity of pixel p is d; the second term represents the minimum value of the previous pixel cost accumulation (including the penalty coefficient) of pixel p in the r direction; the last term minus The minimum cost value in the accumulation direction ensures that Lr(p,d)<C _max (p,d)+P _2. This term has no impact on the selection of the optimal path. It is calculated by accumulating the cost values in different directions. After summing, the final cost accumulation value can be obtained. The cost accumulation summation formula is

After completing the accumulation of the cost matrix, use the Winner Takes All (WTA) algorithm within the disparity search range to obtain the minimum position of the accumulated value, which is the optimal disparity of the current pixel. The formula for optimal disparity selection is shown in Equation (4.4), and D is the final disparity map

2) An improved algorithm based on tSGM. The motivation is that for large-format, high-resolution remote sensing images, it is difficult to apply a semi-global matching algorithm when using the entire window as the parallax search range. On the one hand, the memory overhead required for cost matrix calculation and cost accumulation process is relatively large, and the calculation time is relatively long; on the other hand, since each pixel calculates the matching cost on the entire window, when there are many low textures or Repeating texture information will increase the probability of mismatching.

For the tSGM algorithm, semi-global matching is performed layer by layer on the pyramid image, and the matching results of this layer are used to limit the disparity search boundary of the next layer of matching, thereby dynamically limiting the disparity search range of each pixel. The tSGM algorithm determines the disparity search of the pixel range, first search for the maximum value d _max and minimum value d _min of disparity in the window around the disparity pixel D(x), and then determine whether the maximum and minimum disparity difference exceeds the maximum preset search range R. If d _max -d _min <R, the upper boundary range is set to t _max =d _max -D(x)+2, and the lower boundary range is set to t _min =D(x)-d _max +2. If it exceeds the maximum Default range, at this time the parallax search range is defined as formula (4.3.1-5)

After determining the disparity search boundary of this layer of images, the upper and lower disparity boundaries are doubled as the disparity search range for epipolar line image matching of the next layer, and then semi-global matching is performed layer by layer to obtain the final matching result.

To apply the expansion of the one-dimensional search area in tSGM to the search window, the strategies adopted are different. The cost calculation steps are basically the same. They use the census algorithm to calculate the cost on each layer of the pyramid image. The difference is that , the one-dimensional search calculation cost only needs to calculate the pixels within the search range on the same line, while the two-dimensional search needs to calculate the cost value for the entire search area. The same is that the final cost value results are stored in a one-dimensional array.

The biggest difference is the cost aggregation stage. When obtaining the neighborhood pixel cost aggregation value for a given disparity, the neighborhood at this time is no longer only the left and right neighborhoods along the one-dimensional direction, but includes all surrounding directions. of eight neighborhoods, so the formula for cost aggregation should be changed to:

The direction stores indexes in eight directions, pointing to the eight directions of upper left, upper, upper right, left, right, lower left, lower, and lower right. When performing subsequent cost aggregation calculations, it is necessary to obtain the cost aggregation value of adjacent pixels. , the value of the corresponding position can be obtained directly based on the index.

In the process of selecting the optimal disparity, formula 4.1.1-4 is still used for calculation, and the value with the smallest cost aggregation value is selected from the search window corresponding to each pixel.

The disparity post-processing stage is basically the same, except that when performing sub-pixel interpolation from the plastic disparity, a biquadratic curve can be used for fitting, making the obtained disparity value more accurate.

At the same time, in order to avoid that the search area is still large even though it has been reduced, this article adopts the following strategies for optimization:

a. At the top level of the pyramid, perform a two-dimensional pixel-by-pixel search on the search image for each pixel on the reference image.

b. For pixels with high matching confidence in the previous layer, a fixed-sized square window is directly given in the next layer to search. If the cost aggregation value of the found image point with the same name is lower than the given threshold or the image point with the same name is located on the boundary , expand the search range and search again until the cost value of the same image point is lower than the threshold and does not fall on the boundary.

c. Use the search range of adjacent pixels to narrow the search domain of pixels with the maximum disparity range. This is because the boundary disparity considered untrustworthy in the previous step will cause the corresponding pixels in the next layer to use the maximum disparity search range.

d. Do not search for areas with a mask of 0.

In order to further verify the technical effects of the embodiments of the present disclosure, test examples are given below for illustration.

The first is the reconstruction results of the disclosed technical solution in different terrains. The three-dimensional reconstruction pipeline of this disclosure generated a large number of DSM products of key terrain areas in the satellite stereo image pair test, as shown in Figures 5, 6, 7, and 8 below. From left to right, the left and right ortho-resampled images and generated DSM.

Secondly, for the detailed results of difficult and difficult terrain. In terms of the details of the generated DSM, the effect of comparing STRM 30m resolution DSM products is also very significant. Figure 9 below shows the DSM details of K2, the world’s second highest peak, and the comparison with STRM 30m resolution DEM.

As can be seen from Figure 9, the method proposed in this disclosure can still achieve high-precision and high-resolution DSM generation in a stereo pair with a minimum elevation of 3500m, a maximum elevation of about 8600m, and a single pair of height differences of 5000m. From the detailed figure It can be seen that the DSM generated by the method proposed in this disclosure is basically consistent with the SRTM 30m resolution DEM in terms of rough texture and terrain direction, which fully verifies the correctness and reliability of this method, while in the fineness and clarity of the terrain In terms of accuracy, the method proposed in this disclosure far exceeds SRTM. On the one hand, this method can generate a DSM product with a resolution that is basically consistent with the original image. On the other hand, it also proves that this method is complete in reconstructing terrain at the pixel level. and precise.

In the reconstruction results of key areas, for an area composed of hundreds of images, this method uses pairwise reconstruction and final fusion. Figure 10 below shows the DSM of the Taihu Lake area reconstructed from 120 images, including cities and In the mountainous area, it can be seen that the texture of the entire area is clear, and the reconstructed DSM is completely consistent with the orthophoto.

The quantitative comparison of the embodiments of this disclosure is the DSM elevation accuracy comparison result after comparative reconstruction of 50 scene stereo pairs in selected key areas of the National Resource Three Images and evaluation using ICESAT-1's filtered global laser points .

Attachment: The test time is obtained under the i511400 CPU host + 16G running memory configuration.

The terms and abbreviations used in this disclosure are explained below.

Digital Surface Model (DSM, Digital Surface Model)

Semi-global Matching (SGM)

Digital terrain (or ground) model (DTM, Digital Terrain Model, abbreviation DTM)

Digital Elevation Model (DEM)

Digital Surface Model (DSM)

Digital Orthophoto Map (DOM)

Rational Polynomial Coefficients (RPC: Rational Polynomial Coefficients).

tSGM algorithm (SGM with a tube shaped disparity range, abbreviated as tSGM)

It should be understood that in the embodiment of the present invention, the term "and/or" is only an association relationship describing associated objects, indicating that three relationships can exist. For example, A and/or B can mean: A exists alone, A and B exist simultaneously, and B exists alone. In addition, the character "/" in this article generally indicates that the related objects are an "or" relationship.

Those of ordinary skill in the art can appreciate that the units and algorithm steps of each example described in conjunction with the embodiments disclosed herein can be implemented with electronic hardware, computer software, or a combination of both. In order to clearly illustrate the relationship between hardware and software Interchangeability, in the above description, the composition and steps of each example have been generally described according to functions. Whether these functions are performed in hardware or software depends on the specific application and design constraints of the technical solution. Skilled artisans may implement the described functionality using different methods for each specific application, but such implementations should not be considered to be beyond the scope of the present invention.

In the several embodiments provided in this application, it should be understood that the disclosed systems, devices and methods can be implemented in other ways. For example, the device embodiments described above are only illustrative. For example, the division of the units is only a logical function division. In actual implementation, there may be other division methods. For example, multiple units or components may be combined or can be integrated into another system, or some features can be ignored, or not implemented. In addition, the coupling or direct coupling or communication connection between each other shown or discussed may be an indirect coupling or communication connection through some interfaces, devices or units, or may be electrical, mechanical or other forms of connection.

In addition, each functional unit in various embodiments of the present invention can be integrated into one processing unit, or each unit can exist physically alone, or two or more units can be integrated into one unit. The above integrated units can be implemented in the form of hardware or software functional units.

If the integrated unit is implemented in the form of a software functional unit and sold or used as an independent product, it may be stored in a computer-readable storage medium. Based on this understanding, the technical solution of the present invention is essentially or contributes to the existing technology, or all or part of the technical solution can be embodied in the form of a software product, and the computer software product is stored in a storage medium , including several instructions to cause a computer device (which can be a personal computer, a server, or a network device, etc.) to execute all or part of the steps of the method described in various embodiments of the present invention. The aforementioned storage media include: U disk, mobile hard disk, read-only memory (ROM, Read-Only Memory), random access memory (RAM, Random Access Memory), magnetic disk or optical disk and other media that can store program code.

The above are only specific embodiments of the present invention, but the protection scope of the present invention is not limited thereto. Any person familiar with the technical field can easily think of various equivalent methods within the technical scope disclosed in the present invention. Modifications or substitutions shall be included in the protection scope of the present invention. Therefore, the protection scope of the present invention should be subject to the protection scope of the claims.

Claims (9)

1. The three-dimensional reconstruction method of the remote sensing image is characterized by comprising the following steps of:

receiving a first remote sensing image;

carrying out orthorectification on the first remote sensing image to obtain a second remote sensing image;

resampling the second remote sensing image to obtain a third remote sensing image;

performing adaptive expansion partitioning on the third remote sensing image, and performing dense matching on the partitioned image to generate a corresponding parallax image;

based on the parallax map and the blocking information, restoring each pair of homonymous image points to the coordinates of the homonymous image points on the orthographic correction image, performing inverse orthographic correction on the homonymous image point pairs, restoring the point pair coordinates to the first high-resolution remote sensing image, and obtaining a three-dimensional point cloud formed by the corresponding three-dimensional point coordinates;

rasterizing the three-dimensional point cloud generates a high resolution DSM.

2. The method of claim 1, wherein the first remote sensing image is a high resolution remote sensing image.

3. The method of claim 2, wherein the first remote sensing image comprises a front view image, a bottom view image, and a back view image.

4. A remote sensing image three-dimensional reconstruction method according to claim 3, wherein the forward-looking image, the downward-looking image and the backward-looking image of the first remote sensing image are orthorectified based on a reference low resolution DEM.

5. The method of claim 4, wherein the second remote sensing image comprises a front-looking orthographic image, a lower-looking orthographic image, and a rear-looking orthographic image.

6. The method of claim 1, wherein the third remote sensing image is adaptively expanded and segmented based on a reference DEM elevation value.

7. The method for three-dimensional reconstruction of remote sensing images according to claim 6, wherein the dense matching adopts a two-dimensional window search strategy based on tSGM.

8. The method of claim 1, wherein after recovering the point-to-coordinate system onto the first high resolution remote sensing image, the first remote sensing image pair is subjected to forward intersection by using an RPC model to obtain the corresponding three-dimensional point coordinate system.

9. A computer program product comprising a computer program, characterized in that the computer program is executed by a processor to implement the method of any one of claims 1-8.