A method and system for generating a digital elevation model (DEM) from stereo satellite imagery

Abstract

The application provides a method and system for generating a digital elevation model DEM from stereo satellite images, acquiring same-orbit forward-looking and rear-view panchromatic and multispectral satellite images, correcting RPC system errors of the panchromatic stereo image pair through regional network adjustment, implementing orthographic resampling on the panchromatic and multispectral images after adjustment to obtain L2-level stereo images; then fusing the orthographic multispectral and panchromatic images to generate high-resolution color DOM, simultaneously generating DSM through dense matching and spatial triangulation, and superimposing the DOM and the DSM to generate three-dimensional color point cloud; respectively adopting a deep network to perform semantic segmentation on the three-dimensional point cloud and the DOM, and correcting the optimized classified point cloud through two-dimensional and three-dimensional fusion rules; performing elevation setting, height reduction processing and bottom elevation interpolation for water bodies, woodland and building areas respectively, and finally extracting ground point interpolation to generate high-precision DEM. The application fully considers multi-class feature semantic constraints, and the accuracy in complex terrain areas is significantly improved compared with traditional algorithms.

Description

Technical Field

[0001] This invention relates to the fields of photogrammetry and remote sensing, and more specifically, to a method and system for generating a DEM from stereo satellite imagery. Background Technology

[0002] With the rapid development of industries such as digital villages, smart water conservancy, and low-altitude economy, the demand for DEM data production with large coverage, high precision, and high timeliness is increasing rapidly. DEM data covering the entire province, city, county, and township is an important basic data for building a high-precision and highly realistic digital geospatial foundation.

[0003] Currently, DEM data collection and production methods include ground-based altimetry, aerial photogrammetry, and lidar mapping. However, due to limitations in data collection, production results, and efficiency, these methods cannot meet the requirements for large-scale, full-area DEM production.

[0004] With the increasing availability of stereo satellite imagery, technologies and methods for DEM production based on remote sensing stereo image pairs have emerged, making rapid and low-cost large-scale DEM production increasingly possible. However, due to limitations in the data acquisition accuracy of stereo satellite imagery and bottlenecks in stereo modeling technology, it is currently impossible to directly and automatically generate large-scale, high-precision DEM data based on stereo imagery. Existing methods propose DEM generation based on GF7 stereo image pairs, but technically, they only generate 3D dense point clouds through multi-view stereo image matching and directly generate DEM data, without considering the influence of ground features such as buildings, vegetation, and water bodies on terrain features. This results in the inability to obtain high-precision, realistic terrain in areas with buildings and vegetation. Furthermore, multi-view stereo matching methods may fail to match corresponding points in areas with weak texture, such as water surfaces, failing to guarantee the accuracy and completeness of the terrain in water areas.

[0005] In summary, although various satellites have relatively mature stereo imaging capabilities, there is currently a lack of technologies and methods for generating high-precision DEMs based on stereo satellite imagery. Existing automated DEM generation technologies generally lack the calculation and inference of the actual ground elevation of features such as buildings, vegetation, and water bodies in 3D scenes. After automatically generating DEM data, a large amount of manual editing and correction is required, which cannot directly meet the requirements of various industries for digital geospatial terrain bases in terms of accuracy, timeliness, and cost. Summary of the Invention

[0006] This invention addresses the technical problems existing in the prior art by providing a method and system for generating a DEM from stereo satellite imagery.

[0007] According to a first aspect of the present invention, a method for generating a DEM from stereo satellite imagery is provided, comprising:

[0008] Acquire stereo satellite imagery data from the same orbit, the stereo satellite imagery data including forward-looking and backward-looking panchromatic images, forward-looking and backward-looking multispectral images, and corresponding rational polynomial coefficient RPC files;

[0009] The stereo region network adjustment is performed on the front-view and rear-view panchromatic images to correct the systematic errors of the original RPC parameters and obtain the corrected RPC model.

[0010] Based on the corrected RPC model and the publicly available DEM data, orthorectified resampling is performed on the forward-looking and backward-looking panchromatic images and the forward-looking and backward-looking multispectral images, respectively, to obtain orthorectified forward-looking and backward-looking panchromatic images and orthorectified forward-looking and backward-looking multispectral images.

[0011] The orthophoto and backview multispectral images are fused with the orthophoto and backview panchromatic images to obtain a high-resolution color orthophoto image DOM.

[0012] Based on the corrected RPC model, dense image matching and spatial triangulation are performed on the orthophoto and back-view panchromatic images to generate a digital surface model (DSM).

[0013] The digital surface model (DSM) and the orthophoto DOM are spatially overlaid. The elevation values ​​of each grid point of the digital surface model (DSM) are sampled and the corresponding color values ​​are sampled from the orthophoto DOM to construct a three-dimensional color point cloud.

[0014] A first deep learning network is used to perform semantic segmentation of the three-dimensional color point cloud to obtain three-dimensional semantic labels; a second deep learning network is used to perform feature segmentation of the orthophoto DOM to obtain two-dimensional semantic labels.

[0015] Based on preset two-dimensional and three-dimensional fusion rules, the three-dimensional semantic labels and the two-dimensional semantic labels are fused and corrected to obtain an optimized point cloud containing land feature classification labels;

[0016] Based on the optimized point cloud's feature classification labels and the orthophoto DOM's feature segmentation results, the point cloud elevations in the water area, forest area, and building area are respectively subjected to targeted correction processing.

[0017] Ground points are extracted from the corrected point cloud, and interpolation calculations are performed based on the ground points to generate a digital elevation model (DEM) within the coverage area of ​​the stereo satellite image data.

[0018] Based on the above technical solution, the present invention can also be improved as follows.

[0019] Optionally, stereo region network adjustment is performed on the forward-looking and backward-looking panchromatic images to correct the systematic errors of the original RPC parameters, resulting in a corrected RPC model, including:

[0020] Feature point extraction and homonym matching are performed on the front-view and rear-view panchromatic images to obtain the set of connection points;

[0021] Based on the normalized image coordinate output of the original RPC model, a two-dimensional affine transformation parameter is introduced to construct an image square error compensation model.

[0022] Using the image plane coordinates of the connection point on the foresight and backsight panchromatic images as the observed values, and the two-dimensional affine transformation correction parameters of the foresight and backsight panchromatic images as the unknowns to be solved, the least squares error equation is constructed and iteratively adjusted to solve it.

[0023] Substitute the obtained two-dimensional affine transformation correction parameters back into the original RPC model, and refit the corrected RPC model.

[0024] Optionally, the first deep learning network is a Point Transformer network, and the three-dimensional semantic labels include at least ground points, building points, water surface points, and vegetation points; the second deep learning network is a Vitae network, and the two-dimensional semantic labels include at least water bodies, buildings, woodlands, and grassy areas.

[0025] Optionally, the step of fusing and correcting the three-dimensional semantic labels and the two-dimensional semantic labels based on preset two-dimensional and three-dimensional fusion rules to obtain an optimized point cloud containing land cover classification labels includes:

[0026] The 2D semantic segmentation results are projected onto the 3D point cloud space, so that each 3D point simultaneously obtains a 3D semantic label. and two-dimensional semantic tags ;

[0027] When the three-dimensional semantic label With the two-dimensional semantic tags When consistent, retain the three-dimensional semantic tags. As the final semantic label of a 3D point;

[0028] When the three-dimensional semantic label With the two-dimensional semantic tags When there is inconsistency, the three-dimensional semantic labels within the preset neighborhood of the three-dimensional points are statistically analyzed respectively. With the two-dimensional semantic tags Support rate: Select the semantic label with the higher support rate as the final semantic label of the 3D point.

[0029] When the two-dimensional semantic label The three-dimensional semantic tags are labeled as water bodies. When tagged as other categories, the two-dimensional semantic tag is used. As the final semantic label of a 3D point;

[0030] When the three-dimensional semantic label The two-dimensional semantic tags are labeled as buildings. When marked as a non-building, if the elevation value of a 3D point is higher than the elevation values ​​of surrounding ground points and the difference between the two exceeds a preset threshold, then the 3D semantic tag will be... As the final semantic label for a 3D point.

[0031] Optionally, targeted correction processing can be performed on the point cloud elevation within the water body area, including:

[0032] Based on the optimized point cloud's land cover classification labels and the water body vector mask of the orthophoto DOM, a three-dimensional point set within the water body area is extracted;

[0033] The statistical mean of the elevation of the three-dimensional point set within the water body area is calculated as the average elevation of the water body.

[0034] The elevation values ​​of all points within the water body area are uniformly updated to the average elevation of the water body.

[0035] An edge transition zone is constructed using the water body vector boundary, and weighted elevation smoothing is applied to the edge points within the edge transition zone:

[0036] ;

[0037] in, These are the original elevation values ​​of the edge points. The average elevation of the water body. This represents the smoothed elevation value of the edge points.

[0038] Optionally, targeted correction processing can be performed on the point cloud elevation within the forest area, including:

[0039] Based on the land cover classification labels of the optimized point cloud and the forest vector mask of the orthophoto DOM, a three-dimensional point set is extracted within the forest area;

[0040] An elevation histogram is constructed for the three-dimensional point set, and a threshold segmentation method is used to divide the three-dimensional point set into a low-level suspected ground point set and a high-level tree canopy point set;

[0041] Extract a preset proportion of low points from the set of suspected low-level ground points to calculate the benchmark ground elevation, and calculate the average tree height for the set of high-level tree canopy points;

[0042] For each point in the high-level tree canopy point set, the average tree height is subtracted from its original elevation to restore the true surface elevation of the forest land.

[0043] Optionally, targeted correction processing can be performed on the point cloud elevation within the building area, including:

[0044] Based on the building element vector mask of the orthophoto DOM, extract the polygon boundaries of each building area;

[0045] For each building area, a buffer zone is constructed by extending the polygon boundary outward by a preset distance. Points classified as ground level within the buffer zone are extracted as a set of reference points for the elevation of the building's bottom.

[0046] The base elevation of the building bottom is determined based on the elevation statistics of the set of building bottom elevation reference points.

[0047] At the DSM grid location of the building area, using the building's bottom elevation as a reference, a three-dimensional ground point at the building's bottom is generated using an interpolation method, replacing the original elevation value of the building area.

[0048] Optionally, the step of extracting ground points from the corrected point cloud and performing interpolation calculations based on the ground points to generate a digital elevation model (DEM) within the coverage area of ​​the stereo satellite image data includes:

[0049] From the point cloud after land cover classification and targeted correction, the original bare land points with semantic labels of ground are extracted, and the water body ground points after leveling, the forest land ground points after elevation reduction, and the building bottom ground points generated by interpolation are merged to form a complete ground point set.

[0050] Using the planar coordinates of each point in the ground point set as vertices, an irregular triangular mesh is constructed using a triangulation algorithm;

[0051] At each node of the target DEM grid, linear interpolation is performed based on the irregular triangular network to calculate the elevation value of each node, generating regular grid DEM raster data covering the range of the stereo satellite image.

[0052] According to a second aspect of the present invention, a system for generating a DEM from stereo satellite imagery is provided, comprising:

[0053] The acquisition module is used to acquire stereo satellite image data in the same orbit. The stereo satellite image data includes forward-looking and backward-looking panchromatic images, forward-looking and backward-looking multispectral images, and corresponding rational polynomial coefficient RPC files.

[0054] The correction module is used to perform stereo region network adjustment on the front-view and rear-view panchromatic images, correct the systematic errors of the original RPC parameters, and obtain the corrected RPC model.

[0055] The resampling module is used to perform orthorectified resampling on the forward-looking and backward-looking panchromatic images and the forward-looking and backward-looking multispectral images based on the corrected RPC model and the publicly available DEM data, respectively, to obtain orthorectified forward-looking and backward-looking panchromatic images and orthorectified forward-looking and backward-looking multispectral images.

[0056] The first fusion module is used to fuse the orthophoto forward-looking and backward-looking multispectral images with the orthophoto forward-looking and backward-looking panchromatic images to obtain a high-resolution color orthophoto image DOM.

[0057] The first generation module is used to perform dense image matching and spatial triangulation on the orthophoto and back-view panchromatic images based on the corrected RPC model to generate a digital surface model (DSM).

[0058] The construction module is used to spatially overlay the digital surface model (DSM) with the orthophoto DOM, sample the elevation values ​​of each grid point of the digital surface model (DSM) and sample the corresponding color values ​​from the orthophoto DOM to construct a three-dimensional color point cloud;

[0059] The semantic segmentation module is used to perform semantic segmentation of the three-dimensional color point cloud using a first deep learning network to obtain three-dimensional semantic labels; and to perform feature segmentation of the orthophoto DOM using a second deep learning network to obtain two-dimensional semantic labels.

[0060] The second fusion module is used to perform fusion correction on the three-dimensional semantic labels and the two-dimensional semantic labels based on preset two-dimensional and three-dimensional fusion rules to obtain an optimized point cloud containing land feature classification labels.

[0061] The correction processing module is used to perform targeted correction processing on the point cloud elevation in the water area, forest area and building area based on the land cover classification labels of the optimized point cloud and the feature segmentation results of the orthophoto DOM.

[0062] The second generation module is used to extract ground points from the corrected point cloud, perform interpolation calculations based on the ground points, and generate a digital elevation model (DEM) within the coverage area of ​​the stereo satellite image data.

[0063] According to a third aspect of the present invention, a computer-readable storage medium is provided having a computer program stored thereon, which, when executed by a processor, implements a method for generating a DEM from stereo satellite imagery.

[0064] This invention provides a method and system for generating a DEM from stereo satellite imagery. It acquires panchromatic and multispectral satellite imagery from the same orbit, performs regional network adjustment to correct RPC system errors on the panchromatic stereo image pairs, and then performs orthorectification resampling on the adjusted panchromatic and multispectral images to obtain L2-level stereo imagery. The orthorectified multispectral and panchromatic images are then fused to generate a high-resolution color DOM, while a DSM is generated through dense matching and spatial triangulation. The DOM and DSM are then overlaid to generate a three-dimensional color point cloud. A deep network is used to perform semantic segmentation on the three-dimensional point cloud and DOM, and an optimized classification point cloud is obtained through two-dimensional and three-dimensional fusion rule correction. For water bodies, forest areas, and building areas, elevation leveling, elevation reduction, and bottom elevation interpolation are performed respectively. Finally, ground point interpolation is extracted to generate a high-precision DEM. This DEM fully considers the semantic constraints of various land features such as water bodies, vegetation, and buildings, and significantly improves accuracy in complex terrain areas compared to traditional single algorithms. Attached Figure Description

[0065] Figure 1 A flowchart of a method for generating a digital elevation model (DEM) from stereo satellite imagery, provided in one embodiment of the present invention;

[0066] Figure 2 This is a schematic diagram of the corresponding connection points of the front-view and rear-view panchromatic images according to an embodiment of the present invention;

[0067] Figure 3 This is a schematic diagram of a system structure for generating a digital elevation model (DEM) from stereo satellite imagery, provided in one embodiment of the present invention.

[0068] Figure 4 This is a schematic diagram of the hardware structure of a possible computer-readable storage medium provided by the present invention. Detailed Implementation

[0069] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention. In addition, the technical features of the various embodiments or individual embodiments provided by the present invention can be arbitrarily combined with each other to form feasible technical solutions. Such combinations are not constrained by the order of steps and / or structural composition patterns, but must be based on the ability of those skilled in the art to implement them. When the combination of technical solutions is contradictory or cannot be implemented, it should be considered that such a combination of technical solutions does not exist and is not within the scope of protection claimed by the present invention.

[0070] Currently, the main technical means for producing and acquiring high-precision DEM terrain data are relatively traditional methods such as ground altimetry, aerial photogrammetry, and lidar mapping. These traditional DEM production and acquisition methods have high costs for both field and office work, long acquisition and production cycles, and limited capacity for large-scale terrain production, making it difficult to meet the requirements for large-scale basic terrain acquisition and updating.

[0071] In addition, with the development of satellite constellations, stereo satellite imagery data has gradually become more widespread and has gradually become an important source of surveying and mapping geographic information data that is timely and easy to obtain. However, there is still a lack of mature technical solutions for the automated production of high-precision DEMs based directly on stereo satellite imagery.

[0072] This invention primarily addresses the challenge of efficiently and automatically generating full-area DEM terrain data based on stereo satellite imagery. Furthermore, considering the complexity and diversity of ground features, it employs AI-driven 3D scene feature perception results as guiding constraints. For areas such as buildings, vegetation, and water bodies, it fully considers the two- and three-dimensional spatial relationships between features and infers and constrains the building base height, actual ground points in forest areas, and the planarity of water body terrain. This achieves high-precision DEM generation covering all elements including ground, buildings, vegetation, and water bodies. The accuracy is significantly higher than traditional methods using single stereo matching to generate DEMs, greatly reducing the cost of subsequent manual DEM terrain correction.

[0073] Figure 1 This diagram illustrates a method for generating a digital elevation model (DEM) from stereo satellite imagery according to an embodiment of the present invention. Figure 1 As shown, the method includes the following steps:

[0074] Step 1: Acquire stereo satellite imagery data from the same orbit. The stereo satellite imagery data includes forward-looking and backward-looking panchromatic images, forward-looking and backward-looking multispectral images, and corresponding rational polynomial coefficient RPC files.

[0075] First, L1B-level co-orbit forward-looking and backward-looking stereo satellite imagery data is acquired, including four images: forward-looking and backward-looking panchromatic images and forward-looking and backward-looking multispectral images. In one embodiment, the ground sampling distance (GSD) of the panchromatic image is approximately 0.25m, and the GSD of the multispectral image is approximately 1.0m. The L1B-level imagery is also accompanied by rational polynomial coefficient (RPC) files.

[0076] Step 2: Perform stereo region network adjustment on the front-view and rear-view panchromatic images to correct the systematic errors of the original RPC parameters and obtain the corrected RPC model.

[0077] Understandably, stereo region network adjustment is performed on the fore-see and back-see panchromatic images acquired in step S1 to complete the relative orientation of the stereo image pairs and eliminate systematic errors in the RPC parameters of the original L1B fore-see and back-see images. The specific implementation steps of stereo region network adjustment are as follows:

[0078] (1) Feature point extraction and matching: Feature points were extracted from the front-view and back-view panchromatic images using the SIFT algorithm. Corresponding points were matched using the FLANN fast approximate nearest neighbor search algorithm. Mismatched points were removed using the RANSAC algorithm to obtain a uniformly distributed set of connection points. Let i = 1, 2, ..., N, where N is a positive integer. A schematic diagram of the corresponding connection points of the forward and backward panchromatic images is shown below. Figure 2 As shown.

[0079] (2) Construction of RPC image-side affine error compensation model: The original RPC model of the satellite is subject to systematic geometric deviations due to orbit and attitude disturbances. This invention introduces two-dimensional affine transformation parameters for error compensation based on the normalized image coordinate output of the original RPC. The two-dimensional affine transformation parameters include row direction parameter set and column direction parameter set. The formula for simultaneous calculation of row and column direction parameter correction is as follows:

[0080] ;

[0081] In the formula, a0, a1, a2, b0, b1, b2 are the two-dimensional image-side affine transformation correction parameter set; where a0, a1, a2 construct the row coordinate correction equation, and b0, b1, b2 construct the column coordinate correction equation; r and c are the normalized image-side row and column coordinates obtained from the ground coordinates of the original RPC model; r' and c' are the normalized image point coordinates after affine system error compensation.

[0082] (3) Adjustment solution: The image plane coordinates of the stereo tie points on the foresight and backsight panchromatic images are used as adjustment observations. The two-dimensional affine transformation correction parameters of the two images, namely a0, a1, a2, b0, b1, b2, are used as unknowns to be solved. The standard least squares error equation of photogrammetry is constructed as follows:

[0083] ;

[0084] In the formula, V is the observation residual vector, and B is the coefficient matrix of the error equation. The vector represents the affine parameter corrections, and l is a constant term vector, obtained by subtracting the observed coordinates of the image points from the coordinates calculated by the original RPC model. The least squares criterion is used for iterative solution. The convergence criterion is that the sum of the squared differences between (r', c') after all affine corrections and (r, c) in the original RPC model is less than a preset convergence threshold. After iteration, the obtained affine correction parameters are substituted back into the original RPC polynomial to refit and output the complete set of RPC coefficients after system error correction. To ensure the stability and orientation accuracy of the stereo adjustment solution, a uniform grid constraint is used to arrange stereo connection points in the overlapping area of ​​the two images, with a total of no less than 25 effective connection points across the entire domain. After adjustment, the error per unit weight image point is controlled to be less than 0.5 pixels.

[0085] Step 3: Based on the corrected RPC model and the publicly available DEM data, orthorectify the forward-looking and backward-looking panchromatic images and the forward-looking and backward-looking multispectral images respectively to obtain orthorectified forward-looking and backward-looking panchromatic images and orthorectified forward-looking and backward-looking multispectral images.

[0086] Understandably, this step acquires publicly available DEM data (SRTM), performs orthorectification resampling on the multispectral imagery and the adjusted panchromatic imagery obtained in step S2, eliminates topographic projection differences, and obtains L2-level foresight and backwardsight stereo satellite imagery data. The specific implementation steps of orthorectification resampling are as follows:

[0087] (1) Obtain publicly available SRTM DEM data as a reference elevation surface. The spatial resolution of publicly available SRTM DEM data is about 30m, the elevation accuracy is better than 16m, and the coverage range meets the global land area.

[0088] (2) For each pixel in the target orthophoto (forward and backward multispectral images, and adjusted forward and backward panchromatic images), let its ground coordinates be (X,Y). Based on the corrected RPC model and the elevation value H provided by the SRTM DEM, calculate the corresponding image coordinates (r,c) in the original image. First, take the elevation H of the SRTM DEM at (X,Y) as the initial elevation value, substitute it into the RPC model to calculate the image coordinates (r,c), and then obtain the corresponding gray value in the original image through bilinear interpolation and assign it to the corresponding pixel in the orthophoto.

[0089] (3) Perform the above orthophoto resampling process on the multispectral image and the adjusted front and back panchromatic image respectively to eliminate the projection difference caused by terrain undulation and obtain L2 level front and back panchromatic and front and back multispectral orthophoto images.

[0090] Step 4: Fuse the orthophoto and backview multispectral images with the orthophoto and backview panchromatic images to obtain a high-resolution color orthophoto image DOM.

[0091] Understandably, this step performs image fusion on the orthophoto panchromatic and multispectral images obtained in step S3. The multispectral color image is upsampled to the resolution of the panchromatic image, outputting high-resolution color orthophoto data (DOM). The image fusion uses the Gram-Schmidt (GS) orthogonalization fusion algorithm, and the specific implementation steps are as follows:

[0092] (1) Upsample the low-resolution multispectral image to the resolution of the panchromatic image and resample it using bicubic interpolation to ensure that the multispectral image and the panchromatic image are spatially registered after upsampling.

[0093] (2) The multispectral bands upsampled to panchromatic resolution are assigned weighting coefficients according to the spectral response interval of the satellite multispectral payload, and the equivalent low-resolution panchromatic luminance component PAN is generated by band weighted averaging. low It is used to match the spectral brightness characteristics of real high-resolution panchromatic images.

[0094] (3) Perform Gram-Schmidt (GS) orthogonal transformation fusion processing: Perform forward GS orthogonal decomposition only on the upsampled multispectral dataset, and obtain (GS1, GS2, ..., GS) after decomposition. n After decomposition, the first orthogonal component GS1 and the simulated equivalent panchromatic PAN low The brightness characteristics are consistent, mainly containing brightness information, with the remaining components being pure spectral residual components after removing the brightness information; a high-resolution true panchromatic image (PAN) is used. high Replace the first orthogonal component GS1 to form the decomposed component (PAN). high ,GS2,...,GS n Then, the GS inverse orthogonal transform is performed to reconstruct the multispectral bands, resulting in a high-resolution fused color image. The simplified calculation formula for reconstructing the multispectral bands is:

[0095] ;

[0096] In the formula, The original multispectral i-th band GS orthogonal higher-order residual components; For the fused high-resolution i-th multispectral band; coef i The fused regression coefficients for the i-th band are obtained by combining the low-resolution original multispectral bands with the equivalent simulated panchromatic PAN. low The solution is obtained by least squares linear regression.

[0097] (4) The fused orthophoto DOM has both the high spatial resolution of the panchromatic image and the rich color information of the multispectral image. The resolution of the orthophoto DOM is consistent with that of the panchromatic image.

[0098] Step 5: Based on the corrected RPC model, perform dense image matching and spatial triangulation on the orthophoto and backview panchromatic images to generate a digital surface model (DSM).

[0099] Understandably, based on the RPC model corrected using stereo image pairs, a standard satellite stereo mapping process is employed. Dense parallax fields are obtained through dense image matching, and surface elevation is calculated using spatial triangulation methods to generate a digital surface model (DSM). The specific implementation steps for DSM generation are as follows:

[0100] (1) Epipolar Image Generation: Based on the corrected RPC model, the fore- and back-view panchromatic images are resampled into epipolar image pairs so that corresponding points are located on the same epipolar line, reducing the two-dimensional search problem to a one-dimensional search.

[0101] (2) Dense image matching: The SGM semi-global matching algorithm is used to match the orthophoto and foresight panchromatic stereo image pairs pixel by pixel and generate a disparity map D(i,j).

[0102] (3) Aerial triangulation: Using the corrected RPC model and dense parallax field, the three-dimensional ground coordinates of each matching point pair are calculated by forward intersection. Let the image coordinates of corresponding points on the fore-and-after images be (r1,c1) and (r2,c2) respectively. The RPC models of the two images are combined to form a system of equations. The ground point coordinates (Lat,Lon,H) are solved by the least squares method. After point-by-point calculation, the digital surface model (DSM) is generated by gridding according to the rules. The grid spacing is set to 1m to 2m.

[0103] Step 6: Spatially overlay the digital surface model (DSM) with the orthophoto DOM, sample the elevation values ​​of each grid point of the digital surface model (DSM) and sample the corresponding color values ​​from the orthophoto DOM to construct a three-dimensional color point cloud.

[0104] Understandably, the digital surface model (DSM) and orthophoto (DOM) are spatially overlaid, and 3D elevation point sampling and RGB color assignment are performed to output 3D color point cloud data. The specific implementation steps for generating 3D color point clouds are as follows:

[0105] (1) Elevation point sampling: Each grid node of the DSM regular grid is taken as a three-dimensional point, and its plane coordinates (X,Y) and elevation value Z are obtained.

[0106] (2) Color assignment: For each 3D point (X,Y,Z), bilinear interpolation sampling is performed on the DOM image based on its planar coordinates (X,Y) to obtain the corresponding RGB color value (Color). red Color green Color blue), and assign it to the three-dimensional point.

[0107] (3) Point cloud construction: Combine all 3D points with RGB color information into a 3D color point cloud dataset P = {(X i , Y i Z i , R i G i B i )}, i=1,2,...,M, where M is the total number of DSM grid points. The planar coordinates of the point cloud adopt the same projected coordinate system as DSM and DOM to ensure the spatial correspondence between the 3D point cloud and the 2D image. The generated 3D color point cloud contains both geometric elevation information and spectral color information, providing rich multimodal feature input for subsequent AI semantic segmentation.

[0108] Step 7: Use a first deep learning network to perform semantic segmentation of the three-dimensional color point cloud to obtain three-dimensional semantic labels; use a second deep learning network to perform feature segmentation of the orthophoto DOM to obtain two-dimensional semantic labels.

[0109] Understandably, a 3D point cloud semantic segmentation model based on the Point Transformer network model framework is constructed to support the semantic segmentation of ground points, building points, water surface points, and vegetation points. This model is then used to perform semantic segmentation and classification of the 3D color point cloud obtained in step S6. The specific implementation steps of the Point Transformer-based 3D point cloud semantic segmentation model are as follows:

[0110] (1) Point Transformer self-attention feature extraction: Point cloud coordinates and RGB colors are normalized to the [-1,1] interval. The multi-layer encoder extracts geometric and color features through neighborhood self-attention. The core attention calculation formula is:

[0111] ;

[0112] Q is the query feature matrix, representing the local features extracted from each point in the current point cloud; K is the key feature matrix, representing the features of all reference points in the neighborhood; V is the value feature matrix, corresponding one-to-one with K, storing the effective geometric, color, and elevation features of the neighborhood points. T d is the transpose of matrix K. k The channel dimension of the K feature vectors is used. Based on the spatial distance weighted neighborhood attention weights, the extraction of elevation gradient and terrain change features is enhanced. Softmax() is a normalized activation function that converts the similarity into attention weights between 0 and 1. All weights are summed to 1. Attn(Q,K,V) finally outputs the weighted fused point cloud features.

[0113] (2) Overall loss function of the model: Multi-class cross-entropy loss is used in combination with elevation smoothing regularization term to constrain the classification consistency of adjacent points and suppress minor noise.

[0114] Loss total =Loss crossentropy +λLoss reg

[0115] Among them, Loss total The total model loss is the core optimization objective of backpropagation during network training. crossentropy The cross-entropy loss for multi-class classification is the main loss for semantic segmentation. reg The elevation smoothing regularization loss is used, and λ is the weight hyperparameter of the regularization term, which is used to balance the magnitude of the two losses. The hyperparameter λ is fixed at 0.1. The decoder outputs four types of land feature classification labels: 0=ground point, 1=building point, 2=water surface point, and 3=vegetation point.

[0116] (3) Batch point cloud inference output: The three-dimensional color point cloud is divided into blocks and fed into the trained Point Transformer model. Semantic labels are output point by point and written into the point cloud classification field to obtain the initial three-dimensional classification point cloud.

[0117] Simultaneously, a DOM image feature segmentation model based on the Vitae network model framework is constructed, supporting the segmentation and extraction of water bodies, roads, buildings, woodlands, and grasslands. This constructed DOM image feature segmentation model is then used to segment and extract features from the DOM data obtained in step S4. The specific implementation steps of the Vitae-based DOM image feature segmentation model are as follows:

[0118] (1) Vitae window attention image segmentation: Taking a three-channel color DOM as input, layered window attention extracts two-dimensional multi-scale visual features and outputs six semantic masks: water, road, building, woodland, grassland and bare land; attention is calculated only within the window, reducing the computational power consumption of large-scale image reasoning.

[0119] (2) Joint loss function for image segmentation: To address the imbalance problem of small target samples, a joint loss function of Dice loss and cross-entropy is adopted:

[0120] ;

[0121] ;

[0122] Where P is the prediction mask and G is the ground truth mask. img The overall joint loss of the DOM image segmentation model is the optimization objective of the network backpropagation; Loss cross∑P is the pixel-level cross-entropy loss, the basic classification loss; Dice is the Dice similarity coefficient, which measures the degree of overlap between the predicted mask and the ground truth mask; P is the predicted semantic mask (binary map of roads, water bodies, buildings, woodlands, etc.) output by the model; G is the manually labeled ground truth standard mask; ∑P∩G is the total number of pixels intersecting the predicted region and the real feature region; ∑P is the sum of all pixels in the predicted mask; ∑G is the sum of all pixels in the ground truth mask.

[0123] (3) Vector feature generation: The model inference outputs the raster mask of each feature, and the contour extraction algorithm vectorizes it to generate a set of closed vector surfaces. The vectors are completely matched with the spatial coordinates of DOM and DSM, and serve as the constraint boundary for subsequent elevation correction.

[0124] Step 8: Based on the preset two-dimensional and three-dimensional fusion rules, the three-dimensional semantic labels and the two-dimensional semantic labels are fused and corrected to obtain an optimized point cloud containing land feature classification labels.

[0125] Understandably, this step constructs semantic segmentation and fusion processing rules for 2D and 3D spatial features, corrects the semantic segmentation results of the 3D RGB color point cloud and the feature segmentation results of the DOM image, and outputs the corrected feature point cloud classification results. The specific implementation steps of the semantic segmentation and fusion processing rules for 2D and 3D spatial features are as follows:

[0126] (1) Spatial mapping: The DOM 2D segmentation results are projected onto the 3D point cloud space through spatial mapping relationships, so that each 3D point simultaneously has the 3D semantic label L predicted by the Point Transformer network. 3D The two-dimensional semantic label L predicted by the Vitae network 2D .

[0127] (2) Fusion Correction Rules: The following decision rules are used to fuse and correct the two sets of segmentation results: Rule 1: When L 3D With L 2D If they are consistent, keep L directly. 3D As the final classification result of 3D points. Rule 2: When L 3D With L 2D When there is inconsistency, the support rates of the two sets of labels within the K nearest neighbors (K=20) of the 3D point are calculated separately, and the label with the higher support rate is taken as the final label. Rule 3: For weakly textured regions such as water bodies, the semantic label noise is relatively large because 3D point cloud stereo matching is prone to failure. 2D Marked as water body and L 3D When categorized as other, L should be used first. 2D The judgment result. Rule 4: For building areas, when L 3D Marked as building and L 2DWhen marked as a non-building point, check if the elevation value of that point is significantly higher than that of surrounding ground points (elevation difference greater than 3m). If so, retain L. 3D The building's tag.

[0128] (3) Post-processing filtering: The fused classification results are smoothed by majority voting filtering. The point set within a 2m radius of each 3D point is taken and its category is uniformly updated to the category with the highest frequency in the neighborhood. Isolated noise points are eliminated and the corrected feature point cloud classification results are output.

[0129] Step 9: Based on the land feature classification labels of the optimized point cloud and the feature segmentation results of the orthophoto DOM, the point cloud elevations in the water area, forest area and building area are respectively corrected.

[0130] Understandably, this step involves targeted correction of the point cloud elevations within the identified water, woodland, and building areas.

[0131] In this process, the water body feature areas extracted from the DOM imagery are subjected to average elevation calculation. The Z-values ​​of the point cloud elevations within the water body area are then updated and leveled based on the average elevation to ensure the consistency and planarity of the 3D point cloud heights within the water body area. The specific implementation steps for water body elevation leveling are as follows:

[0132] (1) Water body point cloud extraction: Based on the corrected and classified point cloud output by step S8, combined with the DOM water body vector mask, all three-dimensional points within the mask boundary are extracted to form the water body point set P. water .

[0133] (2) Calculation of the average elevation of the water body: For the water body point set P water Elevation values ​​Z of all points i Statistical analysis was performed, and the average elevation H was calculated after removing outlier points using the truncated mean method. water .

[0134] (3) Elevation leveling operation: The elevation Z values ​​of all point clouds within the water body area are uniformly updated to the average elevation H. water This operation ensures the consistency and flatness of the 3D point cloud height within the water body area, eliminating abnormal fluctuations in water body elevation caused by 3D matching failure.

[0135] (4) Water body edge transition processing: An edge transition zone is constructed by buffering the water body vector boundary outward by a distance equal to the resolution of a panchromatic image grid. Weighted elevation smoothing is applied to the boundary points within the edge transition zone to avoid abrupt elevation changes at the water body edge after leveling. Let the original elevation of the edge point be Z. orig The average elevation of the water body is H water The formula for calculating the elevation after boundary point smoothing is as follows:

[0136] .

[0137] in, This represents the smoothed elevation value of each boundary point within the edge transition zone.

[0138] For forest feature areas extracted from DOM imagery, the average tree height is calculated by combining sparse ground points and the lowest points of trees in the forest. The Z-value of the point cloud elevation in the forest area is then reduced to obtain a 3D point cloud that approximates the true ground elevation of the forest. The specific steps for forest elevation reduction are as follows:

[0139] (1) Extraction of point cloud of the whole forest area: Based on the corrected classification point cloud output by S8, combined with the forest vector mask output by DOM in step S7, all three-dimensional points within the boundary of the forest mask are extracted to form the forest point set P. forest At the same time, two-dimensional and three-dimensional labels are superimposed to constrain and filter points whose three-dimensional labels are vegetation.

[0140] (2) Stratification of ground low point and tree canopy high point: For P forest Elevation Z of all points within i An elevation histogram is constructed, and Otsu's maximum inter-class variance thresholding method is used for adaptive segmentation to divide the point set into a low-level suspected ground point set P. low High-level tree canopy point set P high .

[0141] (3) Solving for the average tree height of the region: from the low-level point set P low Extract the lowest 10% elevation quantiles and calculate the mean as the forest land baseline ground elevation H. ground For the tree canopy point set P high The average elevation of all tree crown points within the area is taken as the mean. Finally, the average tree height H is obtained by subtracting the forest land reference ground elevation from the elevation of all tree crown points within the area and then calculating the mean. tree .

[0142] (4) Canopy point elevation reduction and boundary truncation treatment: only for P high The inner canopy point is subjected to a height reduction operation, which subtracts the area's average tree height H from the original elevation. tree To restore the true ground elevation under the forest canopy. Furthermore, for elevation values ​​lower than P after the elevation reduction. low Points with elevation values, denoted by H ground Truncation corrections should be made to avoid excessive reduction in elevation.

[0143] For the building feature area extracted from the DOM image, ground points within a 0.5m buffer zone around the building boundary are used as references for the building's bottom elevation. A 3D ground point cloud of the building's bottom is then interpolated to ensure the accuracy of the building's bottom elevation. The specific steps for interpolating the building's bottom elevation are as follows:

[0144] (1) Building region extraction: Extract the polygon boundaries of the building feature region {B} based on the DOM segmentation results of step S7. k Let k = 1, 2, ..., K, where K is the number of buildings. Perform a morphological closing operation on each building polygon to fill the internal holes.

[0145] (2) Buffer construction: For each building boundary polygon, extend outward by 0.5m to construct a buffer. Extract the points classified as ground class in step S8 within the buffer as the building bottom elevation reference point set R. k .

[0146] (3) Estimation of the building base elevation: For the reference point set R k The elevation values ​​in the data are statistically analyzed, and the median is taken as the reference elevation H of the building's base. base Using the median instead of the mean can effectively resist the influence of noise points and outliers that may exist in the buffer zone. When there are fewer than 3 reference points, the elevation of the nearest ground point around the building is used as a substitute.

[0147] (4) Point cloud interpolation generation at the bottom of the building: At the DSM grid location within the building area, using H base Using the baseline elevation, the inverse distance weighted interpolation (IDW) method is used to generate three-dimensional ground points at the bottom of the building.

[0148] (5) Replace the original building area elevation value with the ground point generated by interpolation at the bottom of the building to ensure the rationality and smoothness of the building bottom elevation.

[0149] Step 10: Extract ground points from the corrected point cloud, perform interpolation calculations based on the ground points, and generate a digital elevation model (DEM) within the coverage area of ​​the stereo satellite image data.

[0150] Understandably, ground points are extracted from the corrected point cloud classification results, and DEM interpolation is performed using these ground points to generate the final DEM data within the coverage area of ​​the stereo image. This DEM generation process fully utilizes the perceptual constraint correction of multiple ground features in the 3D scene, and has a significant advantage over single DEM generation algorithms in terms of terrain accuracy in areas with buildings, water bodies, and dense vegetation cover. The specific implementation steps of DEM interpolation generation are as follows:

[0151] (1) Ground point extraction and merging: From the point cloud classification results after correction, the original bare ground points with the semantic label of ground are extracted, and the ground points of water bodies after leveling in step S9, the ground points of forest land after elevation reduction, and the ground points of building bottoms generated by interpolation are merged to form a complete ground point set P. ground .

[0152] (2) Construction of irregular triangular network (TIN): with Pground Planar coordinates (X) i , Y i Using vertices 1, 2, and 3, an irregular triangular network (TIN) is constructed using the Delaunay triangulation algorithm. The Delaunay triangulation satisfies the empty circle property, meaning that the circumcircle of each triangle does not contain any other points, thus ensuring the local optimality and terrain fidelity of the triangular network.

[0153] (3) DEM grid interpolation: at each node (X) of the target DEM grid g , Y g At point ), the elevation value is calculated using the TIN linear interpolation method.

[0154] (4) For edge regions outside the TIN coverage area, the nearest neighbor assignment method is used to supplement them.

[0155] The final output is a regular gridded DEM raster data covering the area of ​​the stereo satellite imagery. The grid spacing is set to twice the GSD of the panchromatic imagery to balance the accuracy and data size of the generated DEM. Taking the 4D Gaofen-1 03 stereo satellite imagery as an example, the panchromatic imagery resolution is 0.25m, and the corresponding DEM grid spacing is set to 0.5m.

[0156] See Figure 3 This invention illustrates a system for generating a DEM from stereo satellite imagery according to an embodiment of the present invention. The system includes:

[0157] The acquisition module 301 is used to acquire stereo satellite image data in the same orbit. The stereo satellite image data includes forward-looking and backward-looking panchromatic images, forward-looking and backward-looking multispectral images, and corresponding rational polynomial coefficient RPC files.

[0158] Correction module 302 is used to perform stereo region network adjustment on the front-view and rear-view panchromatic images, correct the systematic errors of the original RPC parameters, and obtain the corrected RPC model;

[0159] The resampling module 303 is used to perform orthorectified resampling on the forward-looking and backward-looking panchromatic images and the forward-looking and backward-looking multispectral images based on the corrected RPC model and the publicly available DEM data, respectively, to obtain orthorectified forward-looking and backward-looking panchromatic images and orthorectified forward-looking and backward-looking multispectral images.

[0160] The first fusion module 304 is used to fuse the orthophoto forward-looking and backward-looking multispectral images with the orthophoto forward-looking and backward-looking panchromatic images to obtain a high-resolution color orthophoto image DOM.

[0161] The first generation module 305 is used to perform dense image matching and spatial triangulation on the orthophoto and back-view panchromatic images based on the corrected RPC model to generate a digital surface model (DSM).

[0162] The construction module 306 is used to spatially overlay the digital surface model (DSM) and the orthophoto DOM, sample the elevation values ​​of each grid point of the digital surface model (DSM) and sample the corresponding color values ​​from the orthophoto DOM to construct a three-dimensional color point cloud.

[0163] The semantic segmentation module 307 is used to perform semantic segmentation of the three-dimensional color point cloud using a first deep learning network to obtain three-dimensional semantic labels; and to perform feature segmentation of the orthophoto DOM using a second deep learning network to obtain two-dimensional semantic labels.

[0164] The second fusion module 308 is used to perform fusion correction on the three-dimensional semantic label and the two-dimensional semantic label based on the preset two-dimensional fusion rules to obtain an optimized point cloud containing land feature classification labels.

[0165] The correction processing module 309 is used to perform targeted correction processing on the point cloud elevation in the water area, forest area and building area based on the land feature classification label of the optimized point cloud and the feature segmentation result of the orthophoto DOM.

[0166] The second generation module 3010 is used to extract ground points from the corrected point cloud, perform interpolation calculations based on the ground points, and generate a digital elevation model (DEM) within the coverage area of ​​the stereo satellite image data.

[0167] It is understood that the system for generating a digital elevation model (DEM) from stereo satellite imagery provided by this invention corresponds to the method for generating a digital elevation model (DEM) from stereo satellite imagery provided in the foregoing embodiments. The relevant technical features of the system for generating a digital elevation model (DEM) from stereo satellite imagery can be referred to the relevant technical features of the method for generating a digital elevation model (DEM) from stereo satellite imagery, and will not be repeated here.

[0168] Please see Figure 4 , Figure 4 This is a schematic diagram illustrating an embodiment of a computer-readable storage medium provided by the present invention. (See diagram below.) Figure 4 As shown, this embodiment provides a computer-readable storage medium 400 storing a computer program 411, which, when executed by a processor, implements the steps of a method for generating a digital elevation model (DEM) from stereo satellite imagery.

[0169] The method and system for generating a digital elevation model (DEM) from stereo satellite imagery provided in this invention have the following advantages:

[0170] 1. Based on the stereo image pairs, the digital surface model (DSM) generated is spatially fused and color space mapped with the color DOM to obtain RGB point cloud data, which serves as the basis for subsequent high-precision 3D ground feature semantic segmentation.

[0171] 2. The latest Point Transformer network framework is used for 3D point cloud ground feature classification, which has higher accuracy and faster efficiency than the traditional CSF ground point filtering algorithm and can accurately separate ground points;

[0172] 3. For weak texture areas of water bodies, a combination of DOM full-element segmentation and 3D point cloud semantic segmentation is adopted to effectively solve the problem of abnormal water body point elevation caused by failure of 3D stereo disparity matching.

[0173] 4. A combined approach of DOM full-element extraction and 3D point cloud semantic segmentation is adopted, which fully leverages the advantages of point clouds and images in 3D scene element classification, further improving the accuracy of point cloud classification and reducing erroneous extraction.

[0174] 5. Under the unified constraint of spatial consistency of all elements in the 3D scene, an adaptive DEM generation strategy was adopted for the ground, water surface, forest area and building bottom area. Compared with the traditional method of interpolating and generating DEM by only relying on ground point classification, the global accuracy is significantly improved.

[0175] 6. The entire DEM generation process is driven by an AI model, and the process is fully automated without human intervention, which verifies that high-precision full-domain DEM data can be produced based on a large range of stereo image pairs.

[0176] 7. The entire process fully utilizes AI intelligent perception and constraint correction of multiple features in the 3D scene, which has a significant advantage over single DEM generation algorithms in terms of terrain accuracy in areas with buildings, water bodies, and dense vegetation coverage.

[0177] It should be noted that the descriptions of each embodiment in the above embodiments have different focuses. For parts that are not described in detail in a certain embodiment, please refer to the relevant descriptions in other embodiments.

[0178] Those skilled in the art will understand that embodiments of the present invention can be provided as methods, systems, or computer program products. Therefore, the present invention can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, the present invention can take the form of a computer program product embodied on one or more computer-usable storage media (including, but not limited to, disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0179] This invention is described with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of the invention. It will be understood that each block of the flowchart illustrations and / or block diagrams, and combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, special-purpose computer, embedded computer, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, generate instructions for implementing the flowchart illustrations. Figure 1 One or more processes and / or boxes Figure 1 A device that provides the functions specified in one or more boxes.

[0180] These computer program instructions may also be stored in a computer-readable storage medium that can direct a computer or other programmable data processing device to function in a particular manner, such that the instructions stored in the computer-readable storage medium produce an article of manufacture including instruction means, which are implemented in a process Figure 1 One or more processes and / or boxes Figure 1 The function specified in one or more boxes.

[0181] These computer program instructions may also be loaded onto a computer or other programmable data processing equipment to cause a series of operational steps to be performed on the computer or other programmable equipment to produce a computer-implemented process, thereby providing instructions that execute on the computer or other programmable equipment for implementing the process. Figure 1 One or more processes and / or boxes Figure 1 The steps of the function specified in one or more boxes.

[0182] Although preferred embodiments of the invention have been described, those skilled in the art, upon learning the basic inventive concept, can make other changes and modifications to these embodiments. Therefore, the appended claims are intended to be interpreted as including both the preferred embodiments and all changes and modifications falling within the scope of the invention.

[0183] Obviously, those skilled in the art can make various modifications and variations to this invention without departing from its spirit and scope. Therefore, if these modifications and variations fall within the scope of the claims of this invention and their equivalents, this invention also intends to include these modifications and variations.

Claims

1. A method for generating a digital elevation model (DEM) from stereo satellite imagery, characterized in that, include: Acquire stereo satellite imagery data from the same orbit, the stereo satellite imagery data including forward-looking and backward-looking panchromatic images, forward-looking and backward-looking multispectral images, and corresponding rational polynomial coefficient RPC files; The stereo region network adjustment is performed on the front-view and rear-view panchromatic images to correct the systematic errors of the original RPC parameters and obtain the corrected RPC model. Based on the corrected RPC model and the publicly available DEM data, orthorectified resampling is performed on the forward-looking and backward-looking panchromatic images and the forward-looking and backward-looking multispectral images, respectively, to obtain orthorectified forward-looking and backward-looking panchromatic images and orthorectified forward-looking and backward-looking multispectral images. The orthophoto and backview multispectral images are fused with the orthophoto and backview panchromatic images to obtain a high-resolution color orthophoto image DOM. Based on the corrected RPC model, dense image matching and spatial triangulation are performed on the orthophoto and back-view panchromatic images to generate a digital surface model (DSM). The digital surface model (DSM) and the orthophoto DOM are spatially overlaid. The elevation values ​​of each grid point of the digital surface model (DSM) are sampled and the corresponding color values ​​are sampled from the orthophoto DOM to construct a three-dimensional color point cloud. A first deep learning network is used to perform semantic segmentation of the three-dimensional color point cloud to obtain three-dimensional semantic labels; a second deep learning network is used to perform feature segmentation of the orthophoto DOM to obtain two-dimensional semantic labels. Based on preset two-dimensional and three-dimensional fusion rules, the three-dimensional semantic labels and the two-dimensional semantic labels are fused and corrected to obtain an optimized point cloud containing land feature classification labels; Based on the optimized point cloud's feature classification labels and the orthophoto DOM's feature segmentation results, the point cloud elevations in the water area, forest area, and building area are respectively subjected to targeted correction processing. Ground points are extracted from the corrected point cloud, and interpolation calculations are performed based on the ground points to generate a digital elevation model (DEM) within the coverage area of ​​the stereo satellite image data. The step of fusing and correcting the three-dimensional semantic labels and the two-dimensional semantic labels based on preset two-dimensional and three-dimensional fusion rules to obtain an optimized point cloud containing land cover classification labels includes: The 2D semantic segmentation results are projected onto the 3D point cloud space, so that each 3D point simultaneously obtains a 3D semantic label. and two-dimensional semantic tags ; When the three-dimensional semantic label With the two-dimensional semantic tags When consistent, retain the three-dimensional semantic tags. As the final semantic label of a 3D point; When the three-dimensional semantic label With the two-dimensional semantic tags When there is inconsistency, the three-dimensional semantic labels within the preset neighborhood of the three-dimensional points are statistically analyzed respectively. With the two-dimensional semantic tags Support rate: Select the semantic label with the higher support rate as the final semantic label of the 3D point. When the two-dimensional semantic label The three-dimensional semantic tags are labeled as water bodies. When tagged as other categories, the two-dimensional semantic tag is used. As the final semantic label of a 3D point; When the three-dimensional semantic label The two-dimensional semantic tags are labeled as buildings. When marked as a non-building, if the elevation value of a 3D point is higher than the elevation values ​​of surrounding ground points and the difference between the two exceeds a preset threshold, then the 3D semantic tag will be... As the final semantic label for a 3D point.

2. The method according to claim 1, characterized in that, The front-view and rear-view panchromatic images are subjected to stereo region network adjustment to correct the systematic errors of the original RPC parameters, resulting in a corrected RPC model, including: Feature point extraction and homonym matching are performed on the front-view and rear-view panchromatic images to obtain the set of connection points; Based on the normalized image coordinate output of the original RPC model, a two-dimensional affine transformation parameter is introduced to construct an image square error compensation model. Using the image plane coordinates of the connection point on the foresight and backsight panchromatic images as the observed values, and the two-dimensional affine transformation correction parameters of the foresight and backsight panchromatic images as the unknowns to be solved in the image square error compensation model, the least square error equation is constructed and iteratively adjusted to solve it. Substitute the obtained two-dimensional affine transformation correction parameters back into the original RPC model, and refit the corrected RPC model.

3. The method according to claim 1, characterized in that, The first deep learning network is a PointTransformer network, and the three-dimensional semantic labels include at least ground points, building points, water surface points, and vegetation points; the second deep learning network is a Vitae network, and the two-dimensional semantic labels include at least water bodies, buildings, woodlands, and grassy areas.

4. The method according to claim 1, characterized in that, Targeted correction processing is performed on the point cloud elevation within the water body area, including: Based on the optimized point cloud's land cover classification labels and the water body vector mask of the orthophoto DOM, a three-dimensional point set within the water body area is extracted; The statistical mean of the elevation of the three-dimensional point set within the water body area is calculated as the average elevation of the water body. The elevation values ​​of all points within the water body area are uniformly updated to the average elevation of the water body. An edge transition zone is constructed using the water body vector boundary, and weighted elevation smoothing is applied to the edge points within the edge transition zone: ； in, These are the original elevation values ​​of the edge points. The average elevation of the water body. This represents the smoothed elevation value of the edge points.

5. The method according to claim 1, characterized in that, Targeted correction processing of point cloud elevations within forest areas, including: Based on the land cover classification labels of the optimized point cloud and the forest vector mask of the orthophoto DOM, a three-dimensional point set is extracted within the forest area; An elevation histogram is constructed for the three-dimensional point set, and a threshold segmentation method is used to divide the three-dimensional point set into a low-level suspected ground point set and a high-level tree canopy point set; Extract a preset proportion of low points from the set of suspected low-level ground points to calculate the benchmark ground elevation, and calculate the average tree height for the set of high-level tree canopy points; For each point in the high-level tree canopy point set, the average tree height is subtracted from its original elevation to restore the true surface elevation of the forest land.

6. The method according to claim 1, characterized in that, Targeted correction processing is performed on the point cloud elevation within the building area, including: Based on the building element vector mask of the orthophoto DOM, extract the polygon boundaries of each building area; For each building area, a buffer zone is constructed by extending the polygon boundary outward by a preset distance. Points classified as ground level within the buffer zone are extracted as a set of reference points for the elevation of the building's bottom. The base elevation of the building bottom is determined based on the elevation statistics of the set of building bottom elevation reference points. At the DSM grid location of the building area, using the building's bottom elevation as a reference, a three-dimensional ground point at the building's bottom is generated using an interpolation method, replacing the original elevation value of the building area.

7. The method according to claim 1, characterized in that, The step of extracting ground points from the corrected point cloud, performing interpolation calculations based on the ground points, and generating a digital elevation model (DEM) within the coverage area of ​​the stereo satellite image data includes: From the point cloud after land cover classification and targeted correction, the original bare land points with semantic labels of ground are extracted, and the water body ground points after leveling, the forest land ground points after elevation reduction, and the building bottom ground points generated by interpolation are merged to form a complete ground point set. Using the planar coordinates of each point in the ground point set as vertices, an irregular triangular mesh is constructed using a triangulation algorithm; At each node of the target DEM grid, linear interpolation is performed based on the irregular triangular network to calculate the elevation value of each node, generating regular grid DEM raster data covering the range of the stereo satellite image.

8. A system for generating a DEM from stereo satellite imagery, characterized in that, include: The acquisition module is used to acquire stereo satellite image data in the same orbit. The stereo satellite image data includes forward-looking and backward-looking panchromatic images, forward-looking and backward-looking multispectral images, and corresponding rational polynomial coefficient RPC files. The correction module is used to perform stereo region network adjustment on the front-view and rear-view panchromatic images, correct the systematic errors of the original RPC parameters, and obtain the corrected RPC model. The resampling module is used to perform orthorectified resampling on the forward-looking and backward-looking panchromatic images and the forward-looking and backward-looking multispectral images based on the corrected RPC model and the publicly available DEM data, respectively, to obtain orthorectified forward-looking and backward-looking panchromatic images and orthorectified forward-looking and backward-looking multispectral images. The first fusion module is used to fuse the orthophoto forward-looking and backward-looking multispectral images with the orthophoto forward-looking and backward-looking panchromatic images to obtain a high-resolution color orthophoto image DOM. The first generation module is used to perform dense image matching and spatial triangulation on the orthophoto and back-view panchromatic images based on the corrected RPC model to generate a digital surface model (DSM). The construction module is used to spatially overlay the digital surface model (DSM) with the orthophoto DOM, sample the elevation values ​​of each grid point of the digital surface model (DSM) and sample the corresponding color values ​​from the orthophoto DOM to construct a three-dimensional color point cloud; The semantic segmentation module is used to perform semantic segmentation of the three-dimensional color point cloud using a first deep learning network to obtain three-dimensional semantic labels; and to perform feature segmentation of the orthophoto DOM using a second deep learning network to obtain two-dimensional semantic labels. The second fusion module is used to perform fusion correction on the three-dimensional semantic labels and the two-dimensional semantic labels based on preset two-dimensional and three-dimensional fusion rules to obtain an optimized point cloud containing land feature classification labels. The correction processing module is used to perform targeted correction processing on the point cloud elevation in the water area, forest area and building area based on the land cover classification labels of the optimized point cloud and the feature segmentation results of the orthophoto DOM. The second generation module is used to extract ground points from the corrected point cloud, perform interpolation calculations based on the ground points, and generate a digital elevation model (DEM) within the coverage area of ​​the stereo satellite image data. The step of fusing and correcting the three-dimensional semantic labels and the two-dimensional semantic labels based on preset two-dimensional and three-dimensional fusion rules to obtain an optimized point cloud containing land cover classification labels includes: The 2D semantic segmentation results are projected onto the 3D point cloud space, so that each 3D point simultaneously obtains a 3D semantic label. and two-dimensional semantic tags ; When the three-dimensional semantic label With the two-dimensional semantic tags When consistent, retain the three-dimensional semantic tags. As the final semantic label of a 3D point; When the three-dimensional semantic label With the two-dimensional semantic tags When there is inconsistency, the three-dimensional semantic labels within the preset neighborhood of the three-dimensional points are statistically analyzed respectively. With the two-dimensional semantic tags Support rate: Select the semantic label with the higher support rate as the final semantic label of the 3D point. When the two-dimensional semantic label The three-dimensional semantic tags are labeled as water bodies. When tagged as other categories, the two-dimensional semantic tag is used. As the final semantic label of a 3D point; When the three-dimensional semantic label The two-dimensional semantic tags are labeled as buildings. When marked as a non-building, if the elevation value of a 3D point is higher than the elevation values ​​of surrounding ground points and the difference between the two exceeds a preset threshold, then the 3D semantic tag will be... As the final semantic label for a 3D point.

9. A computer-readable storage medium having a computer program stored thereon, characterized in that, When executed by a processor, the program implements the steps of the method according to any one of claims 1 to 7.

Images