A method for dense matching and 3D reconstruction of high-resolution remote sensing image mapping conversion

Abstract

The invention discloses a high-resolution remote sensing image mapping conversion dense matching and three-dimensional reconstruction method, comprising the steps of: receiving a first remote sensing image; orthorectifying 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; adaptively extending and blocking the third remote sensing image, and densely matching the blocked images to generate a corresponding disparity map; based on the disparity map and blocking information, restoring each pair of same-name image points to their coordinates on the orthorectified image, then inversely orthorectifying the same-name point pairs, restoring the point pair coordinates to the first high-resolution remote sensing image, and obtaining a three-dimensional point cloud formed by corresponding three-dimensional point coordinates; and rasterizing and sampling the three-dimensional point cloud to generate a high-resolution DSM.

Description

Technical Field

[0001] The invention belongs to the technical field of image processing, and in particular relates to a high-resolution remote sensing image mapping conversion dense matching and three-dimensional reconstruction method. Background Art

[0002] Laser detection and ranging (LiDAR) data and aerial remote sensing images have been widely used to generate large-scale digital surface models (DSM). Large-scale digital surface models (DSM) generated by high-resolution satellite imagery (HRSI) have the advantages of being comparable, cheaper, and more readily available.

[0003] In recent years, with the rapid development of earth observation technology, the earth observation system has been increasingly improved. Foreign high-resolution earth observation remote sensing satellite technology represented by IKONOS, WorldView, and Pleiades has become very mature, and high-resolution civil 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 (such as the three-linear array sensor of Ziyuan-3 can provide forward, downward, and rearward three-dimensional views), thereby providing an effective and reliable data source for three-dimensional information extraction. At present, the application of three-dimensional information extraction from high-resolution linear array satellite stereo observation data includes but is not limited to: digital surface model generation, urban three-dimensional modeling, natural disaster assessment, and target positioning and tracking. Summary of the invention

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

[0005] receiving a first remote sensing image;

[0006] orthorectifying the first remote sensing image to obtain a second remote sensing image;

[0007] resampling the second remote sensing image to obtain a third remote sensing image;

[0008] Adaptively expanding and dividing the third remote sensing image into blocks, and densely matching the divided images to generate corresponding disparity maps;

[0009] Based on the disparity map and the block information, each pair of image points with the same name is restored to its coordinates on the orthorectified image, and then the same-name point pairs are reversely orthorectified, and the point pair coordinates are restored to the first high-resolution remote sensing image to obtain a three-dimensional point cloud formed by the corresponding three-dimensional point coordinates;

[0010] The 3D point cloud is rasterized and sampled to generate a high-resolution DSM. BRIEF DESCRIPTION OF THE DRAWINGS

[0011] The above and other objects, features and advantages of the exemplary embodiments of the present invention will become readily understood by reading the detailed description below with reference to the accompanying drawings. In the accompanying drawings, several embodiments of the present invention are shown in an exemplary and non-limiting manner, in which:

[0012] Figure 1 A schematic diagram of an existing satellite image 3D reconstruction processing flow.

[0013] Figure 2 A schematic flow chart of a method for three-dimensional reconstruction of remote sensing images according to one embodiment of the present invention.

[0014] Figure 3 An example diagram of row and column differences of the same-name image points after orthorectification according to one embodiment of the present invention.

[0015] Figure 4 An example diagram of the cost accumulation process in the disparity space according to one embodiment of the present invention.

[0016] Figure 5 An example diagram of a remote sensing image according to one embodiment of the present invention.

[0017] Figure 6 An example diagram of a remote sensing image according to one embodiment of the present invention.

[0018] Figure 7 An example diagram of a remote sensing image according to one embodiment of the present invention.

[0019] Figure 8 An example diagram of a remote sensing image according to one embodiment of the present invention.

[0020] Fig. 9 An example of a comparison of remote sensing images according to one embodiment of the present invention. The left image is a DSM generated by an embodiment of the present invention, and the right image is a STRM 30m resolution DEM.

[0021] Fig.10 An example diagram of three-dimensional reconstruction of a remote sensing image according to one embodiment of the present invention. DETAILED DESCRIPTION

[0022] Most of the existing satellite image 3D reconstruction software uses an improved version of the semi-global matching (SGM) algorithm to perform dense image matching based on the epipolar constraint. A core technology for extracting 3D information from satellite images is the stereo dense matching technology. The epipolar image is used as the input data for dense matching. Its correctness directly affects the reliability of the matching result, and thus affects the accuracy of the 3D information. Therefore, correct epipolar image input is crucial for dense matching.

[0023] Unlike frame images, the epipolar lines of linear array images are hyperbolic, which brings great challenges to epipolar image generation. Traditional epipolar image generation methods are applicable to a single stereo observation mode (such as co-orbit or different orbits) and lack versatility. In addition, with the increasing number of earth observation satellites, the possibility of forming stereo observations between different platforms is increasing. Therefore, it is necessary to obtain an alternative epipolar image generation method that can be applied to co-orbit, different orbits and different-source stereo image pairs.

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

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

[0026] Currently, most of the software for 3D reconstruction of various satellite images follows the following processing flow:

[0027] Firstly, feature matching is performed on satellite stereo image pairs to obtain image points with the same name, so as to achieve geometric correction of the satellite camera model and obtain accurate satellite image orientation parameters.

[0028] Then, in order to solve the problem that the satellite image frame is too large and the core line is an approximate curve, the satellite image is processed in blocks:

[0029] First, perform epipolar correction on each pair of image blocks, then perform dense matching to obtain a disparity map, and then achieve forward intersection based on the disparity map to obtain the point cloud of each block;

[0030] Finally, each point cloud is fused to obtain the 3D point cloud of the entire area and the rasterized DSM product. Figure 1 shown.

[0031] It can be seen that in the existing scheme, the original satellite stereo image pair must be corrected for the core line before dense matching can be performed. However, due to the large frame size of satellite images, wide coverage, large terrain differences within the images, and differences in various sensor types, there is no universal core line constraint method that can be applied to most scenarios, and it cannot meet the needs when processing large-scale satellite data.

[0032] When densely matching based on epipolar-rectified images, the traditional tSGM algorithm only searches for points of the same name for pixels on the same epipolar line, and uses a multi-level pyramid approach to match from top to bottom, using the matching measurement information of the upper layer to search the search range of the lower pyramid. However, this method is not applicable to situations where the epipolar lines are not on a straight line, and the matching accuracy of epipolar-rectified images with unaligned epipolar lines is very poor.

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

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

[0035] In the first step, the original image is orthorectified based on the reference low-resolution DSM. This step converts the image points to the object coordinate system in longitude and latitude. Subsequently, resampling is performed on the basis of the orthorectified image to obtain a resampled image with the coordinates of the same-name image points constrained within the RPC accuracy range.

[0036] In the second step, the wide-width satellite images are adaptively expanded and divided into blocks based on the reference DEM elevation values, and dense matching is performed on the divided image pairs to generate corresponding disparity maps. The specific dense matching process adopts a two-dimensional window search strategy based on tSGM, and a Gaussian pyramid is used to use the results of the same-name point search in the upper layer to narrow the search range of the lower layer.

[0037] In the third step, each pair of same-name image points is restored to their coordinates on the orthorectified image based on the disparity map and block information, and then the same-name point pairs are reversely orthorectified based on the reference low resolution to restore the point pair coordinates to the original image. Then, RPC is used to perform forward intersection on the original image pair to obtain the corresponding three-dimensional point coordinates.

[0038] Finally, the 3D point cloud is subjected to post-processing operations such as denoising, fusion, and water removal based on the reference DEM, and 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 multiple DSM fusion operations to improve the DSM accuracy. The overall scheme of the disclosed embodiment is as follows Figure 2 As shown. Figure 2 In the case of satellite sensors with only two stereo images, the left and right images are equivalent to the rear view image and the front view image respectively. In actual reconstruction, only two images need to be reconstructed in three dimensions. For satellite sensors with three stereo images, two methods are generally used for reconstruction:

[0039] 1. Randomly select two of the three-view images for reconstruction, where the left and right images are equivalent to the two selected images.

[0040] 2. Use downward and forward views, downward and backward views to reconstruct respectively, obtain two DSMs and then perform DSM fusion.

[0041] Furthermore, for orthorectified image generation and blocking, 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 corrected for the nuclear line before stereo image matching, reducing the two-dimensional matching search to a one-dimensional search. Usually, the nuclear line of a high-resolution remote sensing image is similar to a hyperbola, but it can be approximately regarded as a straight line in a local range. However, on the one hand, in areas with large terrain fluctuations, even small enough blocks cannot constrain all image points with the same name to the same nuclear line. On the other hand, the current best method of extracting feature point pairs based on the image square for basic matrix estimation has a biggest problem, which is that it is difficult to ensure that the feature matching obtains the same-name points accurately enough, and even when using RANSAC for same-name point filtering, it is difficult to remove all erroneous point pairs. This makes it difficult to generate nuclear line images that meet the requirements for all terrains.

[0042] In view of this, the present invention uses rpc-based and reference low-resolution DEM to orthorectify the left and right images respectively. The image after orthorectification can ensure that the image points with the same name are within the neighborhood, which can also play a role in limiting the search range and avoid the problem that the subsequent dense matching algorithm fails due to the failure of the kernel line constraint. When the rpc accuracy is high enough, the search range can be limited to within dozens of pixels, and the same kernel line constraint must also be searched in dozens or even higher one-dimensional spaces in areas with large terrain undulations. The distance between the image points with the same name after orthorectification is as follows: Figure 3 .

[0043] like Figure 3 As shown in the figure, in the Loess Plateau area where the terrain undulation exceeds 1000m, the search range of the same-name image points after orthorectification is still within 5 pixels in rows and columns, which can fully ensure that the entire search range of the same-name image points is within 100 pixels.

[0044] After obtaining the orthorectified images corresponding to the left and right images, resampling is required to unify the resolution and frame size of the orthorectified images. The orthorectified images after resampling have the same geographic coordinate affine parameters, which greatly facilitates the subsequent search for overlapping areas of the left and right images. There is no need to find the corresponding overlapping areas based on rpc. Due to the error of rpc, the overlapping areas obtained in this way often need to be extended at the edge. In the process of resampling, the required DSM resolution is added as an adjustable parameter. According to the adjustable parameter, orthophotos of any resolution size can be generated, and based on this, DSMs of the same resolution can be generated in the subsequent process.

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

[0046] 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 orthophoto blocks and a certain overlap between blocks. By dividing the image into blocks, on the one hand, the disparity search can be limited to a smaller interval, which can reduce the mismatch caused by the large disparity search range in the repeated texture area. On the other hand, it is easy to implement multi-threaded processing, so that each pair of image blocks can be processed as an independent and non-interfering thread, which will greatly reduce the time required to generate the entire DSM.

[0047] Before performing dense matching, you can do the following three data preparations:

[0048] A. RPC adjustment (this process is not necessary). Regional block adjustment based on feature matching can improve the elevation accuracy obtained by front interaction.

[0049] B. Reference DEM of the corresponding area (required). Elevation constraints based on DEM can improve the accuracy of 3D reconstruction.

[0050] C. Even light and color (this processing is not necessary). For image pairs with inconsistent colors, color evenness can improve the accuracy of dense matching.

[0051] The above is the data preprocessing that can be performed, but image pairs without unnecessary preprocessing can also be densely matched, but the accuracy may be lower in some cases (such as large color differences in the image pairs).

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

[0053] 1) The basic SGM algorithm is to establish a global energy function and obtain the potential optimal disparity when the energy function reaches the minimum. The energy function is

[0054]

[0055] The global energy function E(D) consists of a data term and a smoothing term. The first term on the right side of equation (4.1) is the data term, which represents the sum of the matching costs of all pixels, C(p,D p ) indicates that the parallax of the left image pixel p is D p The second term is the smoothing term, which aims to increase the penalty number P1 for the pixel q whose disparity difference of the pixel p is equal to 1. This makes the disparity of adjacent pixels 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 returns 0. The third item is the smoothing item, which gives a larger penalty number P2 to the pixel q whose adjacent disparity difference of pixel p is greater than 1. Usually P2 changes with the change of image gradient, and P2 cannot be too large, otherwise the disparity edge will be over-smoothed.

[0056] 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, c is the number of image columns, and d max ,d min Represents the upper and lower bounds of the disparity search. After calculating the cost matrix, consider accumulating each point on the cost matrix from multiple directions (8 or 16 directions), such as Figure 4 as shown.

[0057] For pixel point p, its accumulation function in the r direction is

[0058]

[0059] 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 accumulation of the previous pixel cost of pixel p in the r direction (including the penalty coefficient); the last term minus the minimum cost value in the accumulation direction ensures Lr(p,d) <C max (p, d) + P2, this item has no effect on the selection of the optimal path. The final cost accumulation value can be obtained by summing the cost accumulation values ​​in different directions. The cost accumulation summation formula is

[0060]

[0061] After the accumulation of the cost matrix is completed, within the disparity search range, the winner-takes-all (WTA) algorithm is used to obtain the position with the minimum accumulation 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.

[0062]

[0063] 2) The improved algorithm based on tSGM. The motivation is that for large-format high-resolution remote sensing images, it is difficult to apply the semi-global matching algorithm when using the entire window as the disparity search range. On the one hand, the memory overhead required for the cost matrix calculation and cost accumulation process is relatively large, and the calculation time is also relatively long. On the other hand, since the matching cost is calculated for each pixel over the entire window, when there is a lot of low-texture or repetitive texture information in the image, the probability of false matching will increase.

[0064] For the tSGM algorithm, by performing semi-global matching layer by layer on the pyramid image, the disparity search boundary for the next layer of matching is restricted using the matching results of the current layer, thereby dynamically limiting the disparity search range for each pixel. When the tSGM algorithm determines the disparity search range of a pixel, it first searches for the maximum disparity d in the window around the disparity pixel D(x). max and the minimum disparity d min , then it determines whether the difference between the maximum and minimum disparities exceeds the maximum preset search range R. If d max - d min < R, at this time, 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 preset range, at this time, the disparity search range is defined as Equation (4.3.1-5).

[0065]

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

[0067] However, when applying the expansion of the one-dimensional search area in tSGM to the search window, the adopted strategy is different. The cost calculation steps are basically the same, which is to use the census algorithm to calculate the cost on each layer of the pyramid image. The difference is that for one-dimensional search to calculate the cost, it only needs to calculate the pixels within the search range on the same row, while for two-dimensional search, it needs to calculate the cost values for the entire search area. The same is that the finally obtained cost value results are all stored in a one-dimensional array.

[0068] The biggest difference is the cost aggregation stage. When obtaining the neighborhood pixel cost aggregation value of a given disparity, the neighborhood is no longer just the left and right neighbors along the one-dimensional direction, but includes eight neighbors in all directions around it. Therefore, the formula for cost aggregation should be changed to:

[0069]

[0070] The direction is the index that stores eight directions, namely upper left, upper, upper right, left, right, lower left, lower, and lower right. When performing cost aggregation calculations later and needing to obtain the cost aggregation value of adjacent pixels, the value of the corresponding position can be obtained directly according to the index.

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

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

[0073] At the same time, in order to avoid the search area being still large when it has been reduced, this paper adopts the following strategies for optimization:

[0074] a. At the top level of the pyramid, a two-dimensional pixel-by-pixel search is performed on the search image for each pixel in the reference image.

[0075] b. For pixels with high matching confidence in the previous layer, a fixed-size square window is directly given in the next layer for search. If the cost aggregation value of the same-name image points found is lower than the given threshold or the same-name image points are located on the boundary, the search range is expanded and the search is continued until the cost value of the same-name image points is lower than the threshold and does not fall on the boundary.

[0076] c. Use the search range of adjacent pixels to narrow the search domain of pixels in the maximum disparity range. This is because the unreliable boundary disparity in the previous step will cause the corresponding pixels in the next layer to use the maximum disparity search range.

[0077] d. No search is performed for areas where the mask is 0.

[0078] In order to further verify the technical effect of the embodiments of the present disclosure, a test example is given below for illustration.

[0079] First, the reconstruction results of the disclosed technical solution in different terrains. The 3D reconstruction pipeline of the disclosed invention generates a large number of DSM products in key terrain areas in the satellite stereo image pair test, as follows Figure 5 , 6, 7, and 8, from left to right are the left and right ortho-resampled images and the generated DSM.

[0080] Secondly, for the detailed results of the difficult terrain, the generated DSM details are also very significant compared with the STRM 30m resolution DSM product, as shown below Fig. 9 Shown is the DSM details of K2, the world's second highest peak, and its comparison with the STRM 30m resolution DEM.

[0081] from Fig. 9 It can be seen that the method proposed in the present disclosure can still generate a DSM with high precision and high resolution in a stereo image pair with a minimum elevation of 3500m, a maximum elevation of about 8600m, and a single pair of elevation differences of 5000m. It can be seen from the detail map that the DSM generated by the method proposed in the present disclosure is basically consistent with the SRTM 30m resolution DEM in terms of general texture and terrain trend, which fully verifies the correctness and reliability of the method. In terms of the fineness and clarity of the terrain, the method proposed in the present disclosure far exceeds the SRTM. On the one hand, this is because the method can generate a DSM product with a resolution basically consistent with the original image, and on the other hand, it also proves that the method is complete and accurate in reconstructing terrain at the pixel level.

[0082] In the reconstruction results of key areas, for areas composed of hundreds of images, this method adopts the method of pair-by-pair reconstruction and final fusion. Fig.10 The DSM of the Taihu Lake area was reconstructed from 120 images, including cities and mountainous areas. It can be seen that the texture of the entire area is clear and the reconstructed DSM is completely consistent with the orthophoto.

[0083] The quantitative comparison of the embodiments of the present disclosure is 50 stereo image pairs were selected from key areas of the national resource three-image system Comparative reconstruction and DSM elevation accuracy comparison results after evaluation using ICESAT-1 screened global laser points fruit.

[0084]

[0085] Appendix: The test time is obtained under the configuration of i511400CPU host + 16G running memory.

[0086] The terms and abbreviations involved in the present disclosure are explained as follows.

[0087] Digital Surface Model (DSM)

[0088] Semi-global Matching (SGM)

[0089] Digital Terrain Model (DTM)

[0090] Digital Elevation Model (DEM)

[0091] Digital Surface Model (DSM)

[0092] Digital Orthophoto Map (DOM)

[0093] Rational Polynomial Coefficients (RPC: Rational Polynomial Coefficients).

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

[0095] It should be understood that in the embodiments of the present invention, the term "and / or" is only a description of the association relationship of the associated objects, indicating that three relationships may exist. For example, A and / or B can represent: A exists alone, A and B exist at the same time, and B exists alone. In addition, the character " / " in this article generally indicates that the associated objects before and after are in an "or" relationship.

[0096] Those of ordinary skill in the art will appreciate that the units and algorithm steps of each example described in conjunction with the embodiments disclosed herein can be implemented in electronic hardware, computer software, or a combination of the two. In order to clearly illustrate the interchangeability of hardware and software, the composition and steps of each example have been generally described in terms of function in the above description. Whether these functions are performed in hardware or software depends on the specific application and design constraints of the technical solution. Professional and technical personnel can use different methods to implement the described functions for each specific application, but such implementation should not be considered to be beyond the scope of the present invention.

[0097] In the several embodiments provided in the present 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 schematic. For example, the division of the units is only a logical function division. There may be other division methods in actual implementation, such as multiple units or components can be combined or integrated into another system, or some features can be ignored or not executed. In addition, the mutual coupling or direct coupling or communication connection shown or discussed can be an indirect coupling or communication connection through some interfaces, devices or units, or it can be an electrical, mechanical or other form of connection.

[0098] In addition, each functional unit in each embodiment of the present invention may be integrated into one processing unit, or each unit may exist physically separately, or two or more units may be integrated into one unit. The above-mentioned integrated unit may be implemented in the form of hardware or in the form of software functional units.

[0099] If the integrated unit is implemented in the form of a software functional unit and sold or used as an independent product, it can be stored in a computer-readable storage medium. Based on this understanding, the technical solution of the present invention is essentially or the part that contributes to the prior art, 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 a number of instructions for a computer device (which can be a personal computer, a server, or a network device, etc.) to perform all or part of the steps of the method described in each embodiment of the present invention. The aforementioned storage medium includes: U disk, mobile hard disk, read-only memory (ROM, Read-Only Memory), random access memory (RAM, Random Access Memory), disk or optical disk and other media that can store program codes.

[0100] The above is only a specific embodiment of the present invention, but the protection scope of the present invention is not limited thereto. Any technician familiar with the technical field can easily think of various equivalent modifications or replacements within the technical scope disclosed by the present invention, and these modifications or replacements should be included in the protection scope of the present invention. Therefore, the protection scope of the present invention shall be based on the protection scope of the claims.

Claims

1. A method for three-dimensional reconstruction of remote sensing images, characterized in that: The following steps are involved: receiving a high-resolution first remote sensing image, the first remote sensing image comprising a forward-looking image, a downward-looking image, and a rearward-looking image, Based on the reference low-resolution DEM, orthorectifying the forward-view image, the downward-view image, and the rear-view image of the first remote sensing image to obtain a second remote sensing image; Correspondingly, the second remote sensing image includes a forward orthophoto, a downward orthophoto and a rearward orthophoto; resampling the second remote sensing image to obtain a third remote sensing image; Adaptively extending and partitioning the third remote sensing image based on the reference DEM elevation value, and densely matching the partitioned images to generate corresponding disparity maps, wherein the dense matching adopts a two-dimensional window search strategy based on tSGM; 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 the same-name point pairs are reversely orthorectified, and the point pair coordinates are restored to the first high-resolution remote sensing image. The first remote sensing image pair is forward intersected using the RPC model to obtain a three-dimensional point cloud formed by the corresponding three-dimensional point coordinates; The three-dimensional point cloud is rasterized and sampled to generate a high-resolution DSM, and the DSM is subjected to hole filling, bilateral filtering and / or multiple DSM fusion operations to improve the DSM accuracy.

2. The method according to claim 1, characterized in that: Before performing the dense matching, the following processing is performed on the divided image map, including: Perform elevation constraints based on DEM, and / or perform block adjustment based on feature matching, and / or perform uniform lighting and color matching for image pairs with inconsistent colors.

3. The method according to claim 2, characterized in that The 2D window search strategy based on tSGM includes: By establishing a global energy function, the potential optimal disparity is obtained when the energy function reaches the minimum. The global energy function is: The global energy function E(D) consists of a data term and a smoothing term. The first term on the right side of formula (1) is the data term, which represents the sum of the matching costs of all pixels. p ) indicates that the parallax of the left image pixel p is D p The second item is the smoothing item, which is to increase the penalty number P1 for the pixel q whose disparity difference value near the image point p is equal to 1. Where T[|D p -D q |=1] is a judgment function, which returns 1 when the judgment condition is true, otherwise it returns 0. The third item is a smoothing item. A larger penalty number P2 is given to the pixel q whose disparity difference of the pixel p is greater than 1. The strategy first calculates the matching cost of each pixel on the left image within the disparity search range to form a cost matrix with a dimension of r×c×(d max -d min ), where r is the number of image rows, c is the number of image columns, and d max d min represents the upper and lower bounds of the disparity search, After calculating the cost matrix, consider accumulating each point on the cost matrix from multiple directions. For pixel point p, its accumulation function in the r direction is In formula (2), 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 accumulation of the previous pixel cost of pixel p in the direction of r; the last term minus the minimum cost value in the accumulation direction ensures that Lr(p, d) <C max (p,d)+P2, The final cost accumulation value can be obtained by summing the cost accumulation values ​​in different directions. The cost accumulation summation formula is: After the accumulation of the cost matrix is ​​completed, the winner-takes-all (WTA) algorithm is used within the disparity search range to obtain the position with the minimum accumulation value, which is the optimal disparity of the current pixel. The formula for selecting the optimal disparity is shown in formula (4). D is the final disparity map.

4. The method according to claim 3, characterized in that Based on 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, thereby dynamically limiting the disparity search range of each pixel.

5. The method according to claim 4, characterized in that When determining the disparity search range of a pixel, First, search for the maximum disparity value d in the window around the disparity pixel D(x) max , minimum value d min ,Then determine whether the maximum and minimum parallax differences exceed the maximum preset search range R, If d max -d min <R, at this time the upper boundary range is set to t max = d max -D(x) + 2, the lower boundary range is set to t min = D(x) - d max + 2, if it exceeds the maximum preset range, at this time the parallax search range is defined as formula (5), After determining the disparity search boundary of the current layer of image, the upper and lower disparity boundaries are doubled as the disparity search range for matching the next layer of epipolar images, and then semi-global matching is performed layer by layer to obtain the final matching result.

6. The method according to claim 5, characterized in that In the two-dimensional window search strategy based on tSGM, in the cost aggregation stage, when obtaining the neighborhood pixel cost aggregation value of a given disparity, the neighborhood at this time includes eight neighborhoods in all directions of the surrounding area, and the formula for cost aggregation is: The direction is the index that stores eight directions, namely, upper left, upper, upper right, left, right, lower left, lower, and lower right. When performing cost aggregation calculations later and needing to obtain the cost aggregation value of adjacent pixels, the value of the corresponding position is obtained directly according to the index.

7. The method according to claim 6, characterized in that In the processing stage after selecting the optimal disparity, the following strategies are adopted for optimization when using biquadratic curve fitting for sub-pixel interpolation from the shaped disparity: (a) At the top level of the pyramid, a two-dimensional pixel-by-pixel search is performed on the search image for each pixel in the reference image; (b) For pixels with high matching confidence in the previous layer, a fixed-size square window is directly given in the next layer for search. If the cost aggregation value of the same-name pixel found is lower than the given threshold or the same-name pixel is located on the boundary, the search range is expanded and the search is continued until the cost value of the same-name pixel is lower than the threshold and does not fall on the boundary; (c) narrowing the search domain of pixels with the maximum disparity range by using the search range of adjacent pixels; (d) No search is performed for areas where the mask is 0.

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

Images