Super large scene modeling method and system based on stereoscopic satellite image pre-positioning refinement

Abstract

This invention relates to the field of data modeling technology and discloses a method and system for modeling ultra-large scenes based on pre-processing of stereo satellite imagery. The method includes: performing joint regional network adjustment and epipolar resampling on stereo satellite image pairs based on orthophotos and aerospace radar elevation reference data to obtain standardized stereo epipolar images; generating global initial digital surface model raster data; using a pixel-level linked editing platform to complete defect identification and batch geometric correction of ground features, outputting a global corrected elevation raster; updating the original dense matching point cloud elevation based on the corrected elevation to obtain an optimized color point cloud; combining vector boundary constraints to perform constrained Deloni triangulation to construct a regular triangular mesh model; completing pixel-level texture fitting; and outputting a realistic 3D model after mesh simplification and texture compression. This method solves the problems of traditional elevation optimization lacking multi-view constraints from stereo imagery and relying solely on numerical elevation correction, which results in insufficient accuracy and difficulty in achieving fine correction of all ground features.

Description

Technical Field

[0001] This invention relates to the field of data modeling technology, specifically to a method and system for modeling ultra-large scenes based on pre-processing of stereo satellite imagery. Background Technology

[0002] As the scale of real-scene 3D construction continues to expand, the demand for large-scale, high-precision real-scene 3D models has significantly increased in fields such as land spatial planning and smart city governance. High-resolution stereo satellite imagery, with its advantages of wide coverage and low acquisition cost, has become the core data source for large-scene 3D modeling. The mainstream modeling process is based on digital photogrammetry, which sequentially completes stereo satellite image pair adjustment, epipolar resampling, and dense matching to generate a digital surface model. Then, triangular meshes are constructed and textures are mapped to obtain a 3D model, which can streamline the production of city-level and terrain-level scenes. However, satellite occlusion, imaging angle deviation, and matching algorithm errors can cause various distortions of ground features in the initial digital surface model, such as jagged building edges, undulating top surface elevations, abrupt changes in water body elevations, distorted forest terrain, and elevation holes in shadow areas. These errors are propagated throughout the entire modeling process, resulting in geometric fragmentation and texture misalignment in the finished model.

[0003] Currently, defect correction methods rely on manual revision of each region after the model is completed, which is costly, has a long production cycle, and makes it difficult to unify the quality of the entire model. Alternatively, elevation can be corrected using a general surface algorithm, which is only suitable for bare terrain and cannot distinguish the features of buildings, water bodies, and forest land. Moreover, this is a post-remedial measure and cannot prevent the transmission of errors from the source.

[0004] In summary, most automated elevation optimization methods lack multi-view geometric constraints from the original stereo imagery. They rely solely on the elevation data itself for numerical correction, resulting in limited accuracy improvement and failing to provide reliable geometric support for the refined correction of all ground features. Summary of the Invention

[0005] In view of this, the present invention provides a method and system for modeling ultra-large scenes based on pre-processing of stereo satellite imagery, in order to solve the problem that related automated elevation optimization methods generally lack multi-view geometric constraints of original stereo imagery, rely solely on the elevation data itself for numerical correction, and have limited accuracy improvement, making it difficult to provide reliable geometric support for the fine correction of all elements of the ground.

[0006] In a first aspect, this invention provides a method for modeling ultra-large scenes based on pre-processing of stereo satellite imagery. This method includes: acquiring true orthophoto data of the target area, original densely matched point clouds, stereo satellite image pair data, basic digital line vector data, orthophotos, and space radar elevation reference data; based on the orthophotos and space radar elevation reference data, performing joint regional network adjustment and epipolar resampling processing on the stereo satellite image pair data to obtain standardized stereo epipolar images; using a dense matching algorithm, generating global initial digital surface model raster data based on the standardized stereo epipolar images; and inputting the standardized stereo epipolar images, global initial digital surface model raster data, and basic digital line vector data into... The pre-built digital surface model and stereoscopic image pixel-level linkage editing platform are used to perform defect classification and identification processing and batch correction processing of geometric constraints of ground features, outputting global corrected digital surface model raster data; based on the global corrected digital surface model raster data, the elevation dimension data of the original dense matching point cloud is updated to generate an optimized color point cloud; using the basic digital line vector as boundary constraint, the constrained Deloni triangulation algorithm is used to construct a regular triangular mesh model based on the optimized color point cloud; pixel-level texture bonding processing is performed on the regular triangular mesh model based on the true orthophoto image to obtain a texture-bonded triangular mesh model; mesh simplification and texture compression are performed on the texture-bonded triangular mesh model in sequence to output a realistic 3D model.

[0007] The ultra-large scene modeling method based on stereo satellite imagery pre-processing provided in this embodiment first acquires true orthophoto data, original dense matching point clouds, stereo satellite image pairs, basic digital line vector data, orthophotos, and space radar elevation reference data for the target area. Based on the orthophotos and space radar elevation reference data, joint regional network adjustment and epipolar resampling processing are performed on the stereo satellite image pairs to generate standardized stereo epipolar images. In the image preprocessing stage, the unified calibration of the global spatial coordinate benchmark is completed, eliminating systematic errors and geometric offsets between multiple stereo images, and standardizing the epipolar geometric relationship between images, so that corresponding image points are on the same scan line. Second, by using a dense matching algorithm, initial digital surface model raster data for the entire region is generated based on the standardized stereo epipolar images. The surface elevation information of the entire region is calculated by matching corresponding points pixel by pixel. After rasterization interpolation, a regular grid elevation dataset is formed, which completely covers the top elevation morphology of various land features such as artificial buildings, natural water bodies, and forest vegetation in the target area, and truly reflects the surface elevation distribution in the initial matching state. Subsequently, by inputting standardized stereo epipolar imagery, global initial digital surface model raster data, and basic digital line vector data into a pre-built digital surface model and stereo image pixel-level linkage editing platform, defect classification and identification processing and batch correction processing of geometric constraints for different land features are carried out. The global corrected digital surface model raster data is output, and the treatment of elevation defects is completed in a concentrated manner before the triangulation mesh construction process. Based on the multi-view pixel information of stereo imagery, an intuitive geometric reference is provided for elevation correction. Differentiated constraint rules are set to match the morphological patterns of different land features, so that the corrected elevation data conforms to the inherent geometric characteristics of various land features, realizing the source treatment of elevation defects and avoiding the transmission and amplification of defects along the production chain. Furthermore, by updating the elevation dimension data of the original densely matched point cloud based on the raster data of the globally corrected digital surface model, an optimized color point cloud is generated. This retains the planar coordinate distribution and color attribute information of the original densely matched point cloud. The elevation values ​​are replaced and optimized point-by-point using planar coordinate mapping, eliminating the need to re-execute the globally corrected 3D dense matching operation. This improves the elevation accuracy of the point cloud while controlling computational resource consumption, ensuring complete adaptation of the point cloud data to the corrected elevation system, and maintaining the naturalness of the original point cloud's planar distribution and the authenticity of its color information. In addition, using basic digital line vectors as boundary constraints, a constrained Deloni triangulation algorithm is employed to construct a regular triangular mesh model based on the optimized color point cloud. The outlines of ground features are set as hard constraint boundaries for triangulation, ensuring that the triangulation process strictly follows the boundary direction of the ground features. This avoids edge fragmentation and misalignment caused by free triangulation, guaranteeing the regularity and fit of the triangular facets at the feature boundaries, and constructing a three-dimensional mesh skeleton with clear geometric boundaries and a regular shape.Then, by performing pixel-level texture bonding processing on the regular triangular mesh model based on true orthophotos, a texture-bonded triangular mesh model is obtained. A precise mapping relationship is established between mesh vertices and image pixels, accurately projecting the color information of the real image onto the 3D mesh surface. Combined with global texture color balancing and seam smoothing, color differences and splicing marks between different image data sources are eliminated, ensuring the visual continuity and realism of the 3D model surface texture, and ensuring precise correspondence between texture information and geometric structure. Finally, by sequentially performing mesh simplification and texture compression on the texture-bonded triangular mesh model, a realistic 3D model is output. The density of triangular faces and the volume of texture data are adjusted according to hierarchical standards. While preserving key geometric features and texture details, the overall data volume is optimized, forming standardized result data adaptable to different loading environments. This improves the scene adaptability and data universality of the final model result, meeting the data loading and rendering needs of different applications. By implementing this invention, the problem that related automated elevation optimization methods generally lack multi-view geometric constraints from original stereo images, relying solely on elevation data for numerical correction, resulting in limited accuracy improvement and difficulty in providing reliable geometric support for the refined correction of all elements of the terrain is addressed.

[0008] Secondly, this invention provides a large-scale scene modeling system based on pre-processing of stereo satellite imagery. The system includes: an acquisition module for acquiring true orthophoto data of the target area, original densely matched point clouds, stereo satellite image pairs, basic digital line vector data, orthophotos, and space radar elevation reference data; a sampling module for performing joint regional network adjustment and epipolar resampling processing on the stereo satellite image pairs based on the orthophotos and space radar elevation reference data to obtain standardized stereo epipolar images; a generation module for generating global initial digital surface model raster data based on the standardized stereo epipolar images using a dense matching algorithm; and a correction module for inputting the standardized stereo epipolar images, global initial digital surface model raster data, and basic digital line vector data into... The system integrates a pre-built digital surface model and a pixel-level linkage editing platform for stereoscopic images. It performs defect classification and identification processing, and batch correction processing based on geometric constraints of different features, outputting raster data of the globally corrected digital surface model. An optimization module updates the elevation dimension data of the original densely matched point cloud based on the raster data of the globally corrected digital surface model, generating an optimized color point cloud. A construction module uses basic digital line vectors as boundary constraints and employs a constrained Deloni triangulation algorithm to construct a regular triangular mesh model based on the optimized color point cloud. A bonding module performs pixel-level texture bonding processing on the regular triangular mesh model based on true orthophotos, obtaining a texture-bonded triangular mesh model. An output module sequentially performs mesh simplification and texture compression on the texture-bonded triangular mesh model, outputting a realistic 3D model. Attached Figure Description

[0009] To more clearly illustrate the specific embodiments of the present invention or the technical solutions in the prior art, the drawings used in the description of the specific embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are some embodiments of the present invention. For those skilled in the art, other drawings can be obtained from these drawings without creative effort.

[0010] Figure 1 This is a flowchart illustrating a specific example of a method for modeling ultra-large scenes based on pre-processing of stereo satellite imagery according to an embodiment of the present invention. Figure 2 This is a principle block diagram of a specific example of a large-scale scene modeling system based on pre-processing of stereo satellite imagery according to an embodiment of the present invention. Detailed Implementation

[0011] To address the technical problems mentioned in the background section, this application provides a method for modeling ultra-large scenes based on pre-processing of stereo satellite imagery, as detailed below. Figure 1 As shown, the method includes: Step S101: Obtain true orthophoto data, raw dense matching point cloud, stereo satellite image pair data, basic digital line vector data, orthophoto and space radar elevation reference data of the target area.

[0012] Specifically, true orthophoto data refers to orthophoto products that have undergone terrain correction and projection aberration correction, and eliminated the tilt and occlusion of ground features. These products possess accurate geographic coordinates and true top morphology of ground features. True orthophoto products covering the entire target area are retrieved from the already produced satellite imagery database. The coordinate references and coverage areas of each image are checked, and the temporal consistency and color uniformity are verified. Images with cloud or snow occlusion or imaging defects are removed. Qualified images are then stitched together to form a complete true orthophoto dataset covering the target area. Typical products include meter-level and sub-meter-level true orthophotos produced by domestic high-resolution optical satellites. Raw densely matched point clouds refer to three-dimensional point cloud data generated based on preliminary matching of stereo satellite imagery, containing the three-dimensional coordinates of the points and... Color information is obtained by retrieving point cloud results from the previous stereo image matching process. The coverage integrity and attribute completeness of each point cloud are checked, and gross errors significantly deviating from the surface morphology are removed. After unifying the coordinate reference, an original densely matched point cloud dataset is formed, which is directly integrated into subsequent elevation optimization and grid construction processes. Stereo satellite image pairs refer to stereo observation combinations composed of two satellite images of the same area from different imaging angles. This is the original data source for elevation calculation. Stereo image pairs covering the target area, whether from the same or different orbits, are selected from the satellite image archive. The imaging time, spatial resolution, and image overlap of each pair are checked. Image pairs meeting the overlap requirements and imaging quality are selected. After unifying the data format and coordinate reference, a stereo satellite image pair dataset is formed. Common data sources include domestically produced high-resolution satellite imagery products with stereo imaging capabilities; basic digital line vector data refers to vector datasets describing the spatial location and attribute information of various geographic features, including core geographic features such as buildings, water systems, roads, and vegetation. This involves retrieving basic geographic information vector results for the corresponding target area, extracting feature layers related to feature refinement, such as building outlines, water system boundaries, road boundaries, and vegetation extent, checking the geometric integrity and attribute accuracy of the features, and forming a basic digital line vector dataset after coordinate benchmark unification. This dataset is used for feature classification, labeling, and boundary constraint reference; orthophotos refer to geometrically calibrated standard orthophoto products with stable planar coordinate accuracy, serving as the benchmark data source for extracting planar control points. Officially released orthophotos covering the target area are retrieved, and the planar accuracy and feature clarity of each image are checked. The images are then stitched together to form a continuous orthophoto dataset covering the entire area. Commonly used products include orthophotos accompanying the national geographic census and annually updated orthophoto products. Aerospace radar elevation reference data refers to elevation data products generated based on synthetic aperture radar interferometry technology. It has uniform elevation accuracy across the entire area and serves as the benchmark data source for extracting elevation control points. Aerospace radar elevation data covering the target area are retrieved, and the coverage integrity and elevation consistency of the data are checked block by block. After the coordinate benchmark and elevation benchmark are converted and aligned, an aerospace radar elevation reference dataset is formed. Typical products include full-coverage elevation data products produced by domestic and foreign radar altimetry satellites.

[0013] Step S102: Based on orthophotos and aerospace radar elevation reference data, perform joint regional network adjustment and epipolar resampling on the stereo satellite image pairs to obtain standardized stereo epipolar images.

[0014] Specifically, step S102 includes: Step a1: Extract planar control points from orthophotos and extract elevation control points from aerospace radar elevation reference data.

[0015] Furthermore, planar control points refer to landmark points with precise planar coordinate information used to constrain the planar geometric accuracy of the image, while elevation control points refer to topographic points with precise elevation values ​​used to constrain the elevation calculation accuracy of the image. First, a preliminary screening of ground feature characteristics across the entire orthophoto is performed. The target area is then divided into grid cells according to the principle of uniform coverage. Within each cell, corner points of clearly defined and fixed-location features are selected as candidate planar control points. Typical selection targets include the turning points of rigid road intersections, the outer corner vertices of regular permanent buildings, the edge corners of paved plazas, and the endpoints of piers of large bridges. Seasonally changing farmland, temporary structures, and concealed points obscured by vegetation are avoided to ensure that the points can be stably identified in the stereo satellite imagery. In the aerospace radar elevation reference data, priority is given to selecting areas with stable terrain and gentle undulations to set up candidate elevation control points, avoiding areas with abrupt changes in terrain such as steep slopes, gullies, and deep pits. The elevation consistency of the surrounding area of ​​each point is verified, candidate points with abnormal elevation changes are eliminated, and alternative points in the same area are added to make up the number, ultimately forming a set of plane control points and a set of elevation control points that are evenly distributed and of reliable quality.

[0016] Step a2: Construct adjustment control datum based on horizontal control points and vertical control points.

[0017] Furthermore, the adjustment control datum refers to a spatial reference framework that provides unified coordinate constraints for the adjustment of stereo satellite image pairs within a regional network, used to unify the coordinate system and accuracy benchmark of multiple images. First, all planar and elevation control points are incorporated into a unified spatial coordinate system, unifying the coordinate projection method and elevation datum of the points, completing coordinate transformation and datum alignment for points from different sources. Control blocks are divided according to point distribution and image orbit coverage, establishing common connecting points between adjacent blocks to ensure coordinate continuity at block boundaries. Point accuracy verification is performed on all control points, calculating the relative coordinate and elevation deviations between adjacent control points. Abnormal control points with deviations exceeding reasonable ranges are identified and replaced with alternative points in the same area, then re-verified until the overall control network accuracy meets requirements. Combining the image overlap distribution and terrain complexity of the target area, the weight allocation of control points at different locations is adjusted to increase the constraint strength in image overlap areas and terrain-complex areas, ensuring uniform and stable constraint strength of the adjustment datum across the entire region, ultimately forming an adjustment control datum with complete coverage and consistent accuracy.

[0018] Step a3: Perform joint regional network adjustment on the stereo satellite image pairs based on the adjustment control benchmark, and identify and retain the stereo satellite image pairs that have passed the adjustment.

[0019] Furthermore, joint regional network adjustment refers to the adjustment operation that uses the adjustment control datum as a constraint to perform a comprehensive solution of the exterior orientation elements and densification point coordinates of multiple stereo satellite image pairs, eliminating systematic errors and geometric misalignments between images. First, the observation equations for regional network adjustment are constructed based on the rational function model parameters of the stereo satellite image pairs. The orbital attitude parameters and camera interior orientation elements of each image are used as unknowns to be solved. The plane control points and elevation control points in the adjustment control datum are used as mandatory constraint observations. Simultaneously, corresponding tie points are automatically extracted in the overlapping areas of adjacent images as additional constraints, constructing a set of regional network adjustment error equations covering the entire region. A bundle adjustment strategy was adopted, using the ray bundle of each image as the basic unit. A functional relationship between image point coordinates and object coordinates was established based on the collinearity condition equation. Global least squares iterative solution was performed on all unknowns. During the iteration, the weights of observations were dynamically adjusted according to the distribution of control point residuals and tie point residuals. Observations with large residuals were gradually reduced in weight until they were removed. After convergence, the optimized exterior orientation elements and the 3D coordinates of the global densification points for each image were obtained. After the adjustment was completed, accuracy assessment was performed. The planar and elevation residuals of each image at all control points and tie points were calculated for each image. The mean and standard error of the residuals for each image were statistically analyzed. The thresholds for acceptable adjustment were set as follows: the standard error of the planar residuals was no greater than 1.0 pixel, and the standard error of the elevation residuals was no greater than 0.5 meters. Stereo satellite image pairs that met the threshold requirements were deemed acceptable and retained. Image pairs that exceeded the threshold requirements were investigated for causes and either re-entered into the adjustment process or were removed. Finally, a set of stereo satellite image pairs with unified coordinate references and acceptable geometric accuracy was formed.

[0020] Step a4: Perform epipolar resampling on the adjusted stereo satellite image pairs to generate standardized stereo epipolar images.

[0021] Furthermore, epipolar resampling refers to the image preprocessing operation of rearranging pixels in the original stereo image based on epipolar geometry principles, so that corresponding points in the left and right images are in the same scan row. Based on the exterior orientation elements of the image obtained from adjustment, the epipolar geometric equations corresponding to the stereo image pair are derived, the epipolar direction and pixel mapping relationship of the left and right images are calculated, and a one-to-one mapping table of pixel coordinates of the left and right images is established. Row-by-row resampling is performed on the pixel array of the original stereo satellite image pair according to the epipolar direction. An adapted interpolation algorithm is used to calculate the gray value of the new pixel to ensure a smooth transition of grayscale information in the resampled image and avoid jagged grayscale distortion. After resampling, multiple sets of corresponding ground features are extracted in the entire image area for row alignment verification. The row number deviation of the corresponding points in the left and right images is checked. For areas where the deviation exceeds the allowable range, the epipolar parameters are corrected and resampled. The geometric consistency and grayscale quality of the epipolar image are verified scene by scene. After confirming that all requirements are met, a standardized stereo epipolar image with normal geometric relationships and strict row alignment of corresponding points is output to meet the input requirements of the subsequent dense matching process.

[0022] Step S103: Using a dense matching algorithm, generate global initial digital surface model raster data based on standardized stereo epipolar images.

[0023] Specifically, step S103 includes: Step b1: Perform pixel-by-pixel matching on the standardized stereo epipolar image using a multi-view dense matching algorithm to generate a dense point cloud with the same name.

[0024] Furthermore, the multi-view dense matching algorithm refers to an image matching algorithm that combines observation information from multiple overlapping stereo images to solve for corresponding image points pixel by pixel and calculate their 3D coordinates. Dense corresponding point clouds refer to a high-density set of 3D coordinate points obtained through pixel-by-pixel matching. First, based on the epipolar geometric constraints of standardized stereo epipolar images, the matching search range and matching cost calculation rules are set. A semi-global matching strategy is used to complete pixel-level matching cost aggregation, traversing all pixels row by row along the epipolar direction to solve for the disparity value corresponding to each pixel. During the matching process, the observation results of multiple overlapping images are integrated. The optimal value of multiple sets of matching results is taken as the final disparity for the multi-view overlapping area. Matching confidence thresholds are set for image edges and occluded areas, and pixel matching results with insufficient confidence are eliminated. Based on the image's interior and exterior orientation elements, the disparity value is converted into 3D spatial coordinates, and each point is simultaneously assigned corresponding image grayscale and color information, ultimately forming a dense corresponding point cloud covering the entire area. Typical matching results can completely restore the elevation morphology of features such as building rooftops, ground slopes, and vegetation canopies.

[0025] Step b2 involves performing statistical filtering and morphological filtering on the dense point cloud to obtain a denoised point cloud.

[0026] Furthermore, statistical filtering refers to a point cloud denoising method that removes outliers based on the statistical distribution characteristics of the point cloud neighborhood, while morphological filtering refers to a point cloud normalization method that removes discrete noise points based on morphological operation rules. First, statistical filtering is performed. A specified neighborhood space is defined for each point, and the mean and standard deviation of the spatial distances of all points within the neighborhood are calculated. Outliers deviating from the normal distribution range are identified based on statistical thresholds, and floating discrete points and edge mismatches caused by matching errors are identified as outliers and removed. After statistical filtering, morphological filtering is performed. A surface morphology reference surface is constructed based on morphological opening and closing operations. The elevation difference between the point cloud points and the reference surface is compared, and discrete noise points exceeding the reasonable elevation difference range are removed. Special attention is paid to matching anomalies caused by water surface reflections and building shadows. After both rounds of filtering, the spatial continuity of the point cloud is verified block by block, and effective points with excessively small holes are filled in. Finally, a denoised point cloud with complete noise removal and morphological continuity is obtained.

[0027] Step b3: Perform rasterization interpolation on the denoised point cloud of the same name to generate global initial digital surface model raster data. The grid spacing of the global initial digital surface model raster data is in the range of [0.5m, 2m].

[0028] Further, rasterization interpolation refers to the spatial interpolation process of converting discrete point clouds into regular grid elevation data. The initial global digital surface model raster data refers to the initial elevation dataset stored in regular grid form, covering the entire target area and containing the top elevation information of ground features. First, the regular grid spacing is set based on production accuracy requirements. The target area is then divided into global units according to the grid spacing, forming a fully covered grid array. For each grid unit, the elevation values ​​of all point cloud points within the unit are statistically analyzed, and the average elevation of the effective points is taken as the initial elevation value of the grid. For blank grid units lacking effective point clouds, an inverse distance weighted interpolation algorithm is used to interpolate and complete the data based on the effective grid elevation values ​​of the surrounding neighborhood. During interpolation, a neighborhood search range and weight decay rules are set, and the neighborhood range is appropriately reduced for areas such as building edges and water system boundaries to preserve the boundary features of ground features. After all grid calculations are completed, the continuity and coverage integrity of the global elevation data are verified, and the initial global digital surface model raster data with a unified format and consistent coordinate reference is output. The grid spacing of the initial digital surface model raster data for the entire domain ranges from [0.5m to 2m]. This range is determined by the spatial resolution of the stereo satellite imagery, the upper limit of the accuracy of dense matching, the representation requirements of all ground features, and the computational efficiency of the production process. It represents a reasonable range that balances accuracy and production efficiency. The ground resolution of mainstream high-resolution stereo satellite imagery is at the sub-meter level. Dense matching can stably output a planar accuracy that corresponds to the elevation accuracy. When the grid spacing is less than the lower limit of the range, it will not substantially improve the true geometric accuracy of the elevation data. Instead, it will introduce random noise from the matching process into the raster data, causing fragmented distortions at the edges of ground features. At the same time, it will significantly increase the data volume and computational cost of subsequent refinement and modeling stages, failing to bring commensurate accuracy gains. When the grid spacing exceeds the upper limit of the interval, the morphological features of small buildings, narrow waterways, and road details will be smoothly lost, and the local undulations of the terrain cannot be fully preserved. This fails to meet the requirements for feature reproduction in a full-element realistic 3D model and will also be incompatible with the resolution of the true orthophoto image used for subsequent texture mapping, causing geometric misalignment during texture mapping and affecting the visual effect and geometric accuracy of the final model. This interval is compatible with the matching capabilities of most mainstream high-resolution stereo satellite imagery, takes into account the morphological expression accuracy of various features such as man-made buildings, natural waterways, and woodland terrain, while keeping the data volume within a reasonable range to meet the computing resource allocation and production cycle requirements of large-scale production.

[0029] Step S104 involves inputting standardized stereo epipolar imagery, global initial digital surface model raster data, and basic digital line vector data into a pre-built digital surface model and stereo image pixel-level linkage editing platform. Defect classification and identification processing and batch correction processing of geometric constraints of ground features are performed, and global corrected digital surface model raster data is output.

[0030] Specifically, step S104 includes: The aforementioned pre-built digital surface model and stereoscopic image pixel-level linkage editing platform adopts a three-layer technical architecture: data layer, service layer, and application layer. The data layer is based on a hybrid storage of geospatial database and file system, uniformly managing multi-source data such as stereo epipolar imagery, digital surface model raster, and digital line vector, and supports reading and writing standard geographic data formats such as GeoTIFF, IMG, LAS, and SHP. The service layer provides three core computing services: spatial registration engine, defect identification engine, and batch correction engine. Each service exposes call interfaces through RESTful APIs, supporting parameterized configuration and task queue scheduling. The application layer provides four functional modules: data loading, visual browsing, interactive editing, and result export, supporting simultaneous display of multi-source data and pixel-level position synchronization. The interaction between the platform and this method is as follows: This method inputs standardized stereo epipolar imagery, global initial digital surface model raster data, and basic digital line vector data through the platform's data loading interface. After the platform automatically completes spatial geometric registration, it calls the defect identification engine to perform anomaly detection and classification labeling, and then calls the batch correction engine to perform geometric constraint correction according to the land cover category. Finally, it outputs global corrected digital surface model raster data through the results export interface. The entire process supports parameter configuration and progress callback.

[0031] Step c1 involves performing spatial geometric registration on the standardized stereo epipolar image, the global initial digital surface model raster data, and the basic digital line vector data to obtain a linked dataset.

[0032] Furthermore, spatial geometric registration refers to the preprocessing operation that unifies multi-source geographic data into the same spatial coordinate system to achieve complete spatial alignment of features with the same name. Linked datasets refer to multi-source data combinations that have one-to-one spatial correspondence and can be called synchronously. First, the standardized stereo epipolar image, the initial digital surface model raster data, and the basic digital line vector data are uniformly converted to the same legal spatial coordinate system, completing the alignment conversion between the coordinate projection method and the elevation datum. Then, corresponding ground features are uniformly selected as registration control points throughout the entire area. Typical control points include the center point of a cross-shaped road intersection, the outer corner vertex of a rectangular permanent building, the bend point of a main water system, and the corner points of a large hardened site, etc., which are ground features with clear characteristics and fixed locations. Based on the selected registration control points, spatial transformation parameters are calculated, and spatial transformation adjustments are made to the standardized stereo epipolar image, the initial digital surface model raster data, and the basic digital line vector data, respectively. The planar position deviation of the corresponding ground features is checked point by point, and the spatial transformation parameters are iteratively optimized until the deviation is within a reasonable range. Finally, the pixel positions of the stereo image, the grid positions of the digital surface model, and the feature positions of the vector data achieve pixel-level precise correspondence. At any spatial location, the grayscale color information, elevation values, and ground feature vector attribute information of the image can be retrieved simultaneously, forming a linked dataset with spatial alignment and attribute linkage.

[0033] Step c2: Perform defect classification and identification processing on the linked dataset to generate elevation anomaly area data with land feature category labels.

[0034] Specifically, step c2 includes: Step d1 involves using an edge detection algorithm and a deep learning semantic segmentation model to perform anomaly detection on the linked dataset and identify elevation anomaly areas. These areas include jagged edges and undulating top surfaces of buildings, abrupt changes in water elevation, uneven areas of forest canopies, and elevation voids caused by shadows.

[0035] Furthermore, edge detection algorithms refer to image detection algorithms that identify edge distortions by calculating changes in elevation gradients, while deep learning semantic segmentation models refer to intelligent recognition models that achieve automatic classification of ground features and anomaly identification based on labeled sample training. The deep learning semantic segmentation model adopts an encoder-decoder architecture. The encoder uses ResNet-50 as the backbone network to extract multi-scale image features and elevation features. The decoder restores spatial resolution through stepwise upsampling and skip connections, outputting a pixel-level classification result map of the same size as the input image. The model training uses a weighted combination of cross-entropy loss and Dice loss as the loss function, with a cross-entropy loss weight of 0.7 and a Dice loss weight of 0.3 to alleviate the imbalance problem of ground feature category samples. The training sample annotation standard is pixel-level manual annotation, with annotation categories including buildings, water bodies, woodlands, and shadow occlusion. The annotation boundaries are drawn along the actual contours of the ground features, and the size of a single annotation sample is not less than 1024×1024 pixels. The training set and validation set are divided in a 7:3 ratio. Data augmentation strategies such as random flipping, random cropping, and color dithering are used during the training process. The model convergence criterion is that the cross-intersection over union (CUI) of the validation set is not less than 0.85. First, the global elevation gradient value is calculated grid-by-grid based on the initial digital surface model raster data. Edge detection algorithms are used to capture edge stripes with abrupt elevation changes. Then, abnormal edges are screened by combining the normal morphological characteristics of ground features. Jagged abrupt stripes that do not conform to the regular shape of artificial buildings are identified as jagged areas of building edges, and areas with discrete elevation gradient distribution within the top surface are identified as undulating areas of building top surfaces. At the same time, the grayscale information of the standardized stereo epipolar image and the elevation information of the initial digital surface model raster data are synchronously input into the deep learning semantic segmentation model. The model automatically segments the spatial range of four types of ground features: buildings, water bodies, woodlands, and shadow occlusion. Then, combined with elevation statistical features, abnormal states within the corresponding ranges are identified. For example, elevation variance exceeding the normal water body fluctuation threshold within the water body range is identified as a water body elevation jump area. Local elevation variation coefficient exceeding the normal canopy fluctuation range in the woodland area is identified as a woodland canopy concave-convex area. The area of ​​building shade and mountain shadow without continuous effective elevation values ​​is identified as a shadow occlusion elevation void area. The integration rules for edge detection and semantic segmentation recognition results are as follows: First, the union of the two types of results is taken to determine the candidate anomaly range. Then, the intersection of the overlapping areas of the two types of results is taken as the high-confidence anomaly area for direct confirmation. For areas identified only by edge detection but not covered by semantic segmentation, a secondary verification is performed based on the prior morphological rules of the land cover in the area. Those that meet the morphological characteristics of the land cover anomaly are retained, and those that do not are removed. For areas identified only by semantic segmentation but not responded to by edge detection, a verification is performed based on the statistical characteristics of the elevation gradient in the area. Those with abnormal elevation gradients are retained, and those with normal elevation gradients are determined to be misidentified by semantic segmentation and removed. For areas where the two types of results conflict in category classification, the land cover category of semantic segmentation is used as the standard, and the anomaly type is confirmed by combining the elevation gradient characteristics.The final output shows the complete spatial extent of four types of elevation anomalies: jagged edges and undulating top surfaces of buildings, abrupt changes in water elevation, uneven areas of forest canopies, and areas with elevation voids due to shadows.

[0036] Step d2 involves reviewing and supplementing the elevation anomaly areas to obtain the reviewed elevation anomaly areas.

[0037] Furthermore, review and defect supplementation refer to automated defect quality control operations based on multi-view image cross-validation and multi-source data constraint verification. These operations are used to eliminate misjudgments and supplement omissions in the preliminary identification results of the algorithm, thereby improving the accuracy and completeness of abnormal area identification. A cross-validation benchmark is constructed based on the independent matching elevation results of multi-view standardized stereo epipolar imagery. For each preliminarily identified elevation anomaly area, the consistency of the independently calculated elevation data from multiple views is compared. Differentiated elevation deviation tolerance thresholds are set according to land cover categories: a multi-view elevation deviation of no more than 0.3 meters is considered consistent for building areas, no more than 0.2 meters for water areas, no more than 0.8 meters for forest areas, and no more than 0.5 meters for shaded areas. If the multi-view data all show that the elevation status conforms to the normal morphological characteristics of the land cover and the deviation is within the corresponding tolerance threshold range, it is determined to be a normal area misidentified by the algorithm and is automatically removed from the anomaly range. For example, areas where the inherent morphology of land cover, such as the edge of a steep slope formed by natural terrain or a local protrusion of a natural mountain, is misidentified as anomaly are automatically excluded from the anomaly range. A multi-scale detection strategy is employed to perform a secondary traversal detection on the linked dataset, adjusting the detection threshold and window size to identify subtle defects missed in the initial detection. Simultaneously, based on the feature boundary constraints of the basic digital line vector data, locations prone to matching distortions, such as the junctions of adjacent buildings, the edges of building eaves, and the boundaries of small ponds and dams, are located. Combining elevation gradient features and prior feature morphology rules, subtle defects such as adhesion distortion, elevation misalignment, and water surface disturbances missed by the initial annotation algorithm are automatically supplemented. The boundary contours of anomaly areas are automatically corrected simultaneously, ensuring that the boundaries of anomaly areas perfectly match the spatial distribution of actual defects. After completing full-domain cross-validation and supplementary detection, all valid anomaly areas are automatically merged, overlapping areas with duplicate annotations are eliminated, broken anomaly strips are completed, and the boundary morphology of anomaly areas is automatically standardized, resulting in a verified elevation anomaly area with accurate boundaries, no missed detections, and no misjudgments.

[0038] Step d3: Based on the land feature attributes, perform classification and labeling on the verified elevation anomaly areas to generate elevation anomaly area data with attached land feature category labels. The land feature category labels include building category labels, water body category labels, forest land category labels, and shadow occlusion category labels.

[0039] Furthermore, the land feature attributes originate from the feature classification attribute information attached to the basic digital line vector data. Classification labeling refers to the process of matching the corresponding land feature category and adding attribute tags to each elevation anomaly area. The verified elevation anomaly areas are spatially overlaid and matched with the basic digital line vector data. Based on the spatial correspondence, the classification attribute information of the corresponding land feature elements is associated with each elevation anomaly area. For anomaly areas spanning multiple land feature categories, they are classified according to the dominant land feature category with the largest proportion within the area, ensuring that each anomaly area corresponds to a unique land feature category. After completing the attribute association, corresponding land feature category tags are added to each elevation anomaly area, corresponding to building category tags, water body category tags, forest land category tags, and shadow occlusion category tags, respectively. A one-to-one mapping relationship is established between the anomaly areas and the subsequent corresponding correction rules, allowing subsequent correction processes to directly call the matched correction parameters and processing logic based on the category tags. Ultimately, this results in elevation anomaly area data with complete attributes and clear categories, attached with land feature category tags.

[0040] Step c3 involves performing batch correction processing on the elevation anomaly area data based on geometric constraints of ground features to generate raster data of the global corrected digital surface model.

[0041] Specifically, step c3 includes: Step e1: For the elevation anomaly areas with attached building category markers, perform edge straightening constraints and top surface elevation smoothing.

[0042] Furthermore, edge straightening constraint refers to the correction operation of straightening jagged distortion at the edges according to the regular geometric shape of the building, while top surface elevation smoothing refers to the correction operation of eliminating random elevation fluctuations on the building's top surface. Based on the building outline elements in the basic digital line vector data, the coordinates of each corner point and the direction parameters of each edge line of each building are extracted. The correction range of the building edge is defined based on the vector outline. For the elevation grid with jagged distortion within the range, the elevation values ​​are adjusted grid by grid according to the straight direction of the corresponding edge line, so that the building edge forms a continuous and straight geometric shape, and the geometric turning features of the building corner points are strengthened at the same time, eliminating the elevation tail phenomenon at the building base. For the undulating area of ​​the building's top surface, the elevation values ​​of the top surface grid are adjusted using a neighborhood weighted smoothing algorithm to eliminate random elevation fluctuations caused by dense matching, while retaining the reasonable elevation differences of the building's roof facilities, such as the elevation features of small structures such as roof stairwells, elevator machine rooms, and water tanks, avoiding excessive smoothing and loss of ground feature details. Finally, the elevation correction of all abnormal areas of buildings is completed.

[0043] Step e2: For areas with abnormal elevations marked by water body categories, perform elevation unification for static water areas and longitudinal slope constraint for flowing water areas.

[0044] Furthermore, the unification of elevation in static water bodies refers to the operation of correcting the elevation of static water bodies such as lakes, reservoirs, and ponds to a uniform horizontal elevation, while the longitudinal slope constraint in flowing water bodies refers to the correction operation of setting a gentle longitudinal slope according to the natural flow direction of the river. For abnormal areas with attached water body category markers, static and flowing water bodies are first distinguished based on the water system element attributes in the basic digital line vector data. For static water bodies such as lakes, reservoirs, and ponds, multiple connecting elevation points are uniformly sampled along the shoreline around the water body, and a uniform horizontal elevation value that matches the shoreline connection height is calculated. All elevation grids within the water body are uniformly assigned this elevation value to ensure that the water surface is flat and smoothly connected to the surrounding shoreline. For flowing water bodies such as rivers and ditches, multiple control sections are selected along the natural direction of the river, and the corresponding elevation values ​​at each section are extracted. A gentle longitudinal slope change rate that conforms to the natural hydrological law is set, and the elevation value of each elevation grid is calculated segment by segment along the flow direction to make the water surface elevation transition continuously and gently along the flow direction, completely eliminating elevation jumps and turbulent distortions.

[0045] Step e3: Perform canopy morphology filtering on the elevation anomalies of the attached forest land category markers.

[0046] Furthermore, canopy morphology filtering refers to a filtering and correction operation that suppresses local elevation irregularities caused by individual tree canopies while preserving the overall topographic trend of the forest land. For elevation anomaly areas with attached forest land category markers, an appropriate filter window size is set based on the tree species and average canopy size of the forest land. The filter window is appropriately increased in broad-leaved forest areas and appropriately decreased in coniferous forest areas. Morphological opening operations are used to first remove local elevation peaks caused by individual tree canopy protrusions, and then morphological closing operations are used to fill local depressions formed by gaps between tree canopies. During the filtering process, the filtering intensity is adaptively adjusted according to the topographic slope of the forest land area. The filtering intensity is appropriately reduced in mountainous areas with steep slopes to avoid destroying the macroscopic topographic undulation characteristics of the mountain slope. After filtering is completed, the spatial continuity of forest land elevation is verified block by block, fully preserving the macroscopic topographic morphology such as mountain slopes and valley orientations, so that the forest land surface presents a continuous and gentle overall canopy morphology, eliminating the fragmented irregularities and elevation disturbances caused by individual tree canopies.

[0047] Step e4: For elevation anomaly areas with attached shadow occlusion category markers, fill elevation holes using multi-view elevation information from stereo epipolar imagery.

[0048] Furthermore, reusing multi-view elevation information from stereo epipolar images refers to the operation of filling occlusion holes using elevation data extracted from stereo images at different imaging angles. For elevation anomaly areas marked with shadow occlusion categories, the spatial range corresponding to the hole is first located. Standardized stereo epipolar images covering all imaging angles within this spatial range are retrieved. The occlusion status of the corresponding position in each image is determined from each viewpoint. Unobstructed and clearly imaged viewpoints are selected, and the effective matching elevation values ​​at the corresponding positions are extracted. The extracted effective elevation data is then used to fill the elevation grid of the hole area. For areas that are completely occluded and for which no effective elevation can be extracted from any imaging angle, the elevation values ​​within the hole area are calculated using a trend interpolation algorithm based on continuous elevation data around the hole and along the natural direction of the surrounding terrain. This ensures a smooth transition between the filled elevation and the surrounding terrain, completely eliminating elevation holes and data gaps in the shadowed areas.

[0049] Step e5: Generate raster data for the global corrected digital surface model.

[0050] Furthermore, the fully corrected digital surface model raster data refers to the final elevation raster result with consistent accuracy and continuous morphology after the correction of all land features in the entire region. After completing the zonal correction of the four types of land features, a smooth transition process is applied to the boundary areas connecting different land features, such as the junctions between buildings and the surrounding ground, water bodies and shorelines, and forest land and open land. The elevation values ​​of the boundary areas are adjusted grid by grid to eliminate abrupt elevation changes caused by differences in different correction rules. Subsequently, a full-area integrity verification is carried out, checking for any uncovered blank areas and residual defects block by block according to the map sheet, and supplementing and correcting local problem points. Finally, the grid size, coordinate projection, elevation datum, and data storage format of the full-area data are unified to ensure that the full-area data standards are completely consistent, and the final output is fully corrected digital surface model raster data that completely covers the target area, has regular and natural land feature morphology, and consistent elevation accuracy.

[0051] Step S105: Update the elevation dimension data of the original dense matching point cloud based on the global correction digital surface model raster data to generate an optimized color point cloud.

[0052] Specifically, step S105 includes: Step f1: Extract the planar coordinate information and color attribute information of the original densely matched point cloud to generate basic attribute data of the point cloud.

[0053] Furthermore, point cloud basic attribute data refers to the set of point cloud attributes that retains the original point cloud planar spatial distribution and color features while removing the original elevation values. It is the core carrier for subsequent elevation updates. The process involves reading all attribute fields from the original densely matched point cloud, separating the planar coordinate components and elevation components from the 3D coordinate information, and extracting the x-coordinate and y-coordinate values ​​of each point as planar coordinate information. Simultaneously, the image color information corresponding to each point is extracted, including grayscale values ​​in the visible light band as color attribute information. During the extraction process, the number and spatial arrangement of points in the original densely matched point cloud are fully preserved. The planar position and color assignment of any point are not changed; only the original elevation values ​​are removed, without any interpolation or resampling of the planar distribution. For example, various points covering building rooftops, natural surfaces, vegetation canopies, and water surfaces retain the original matched planar positions and corresponding image colors, only removing original elevation values ​​with errors. This ultimately forms point cloud basic attribute data with complete attributes and unchanged distribution.

[0054] Step f2: Establish the planar coordinate mapping relationship between the global corrected digital surface model raster and the basic attribute data of the point cloud.

[0055] Furthermore, the planar coordinate mapping relationship refers to the positional correspondence and indexing relationship between the grid cells of the globally corrected digital surface model raster and the points in the basic attribute data of the point cloud, established based on the same planar coordinate system. This is used to quickly retrieve the corrected elevation value corresponding to each point in the point cloud. First, the coordinate reference and projection method of the globally corrected digital surface model raster and the basic attribute data of the point cloud are verified to confirm that they are in a completely consistent spatial reference system. Each point in the basic attribute data is traversed, and based on the planar coordinate value of the point, the row and column position corresponding to the point in the globally corrected digital surface model raster is calculated to determine the specific grid cell into which the point falls. For points falling on grid boundaries, the correspondence is determined using a neighboring grid weight allocation method. A corresponding grid index entry is established for each point, constructing a globally unified planar coordinate mapping index table. For example, point cloud points at the edge of a building correspond to the corrected elevation grid of the building edge, and water surface points correspond to the corrected elevation grid of the water body, ensuring that all points can be accurately matched to the corresponding elevation grid cells through the mapping relationship.

[0056] Step f3 involves replacing the elevation values ​​of the original densely matched point cloud point by point based on the planar coordinate mapping relationship to generate an optimized color point cloud.

[0057] Furthermore, the process iterates through all points in the basic attribute data of the point cloud, and based on the established planar coordinate mapping relationship, it retrieves the elevation values ​​of the corresponding raster grid of the global corrected digital surface model, directly replacing the original elevation values ​​of the corresponding points in the original densely matched point cloud. During the replacement process, the planar coordinate and color attribute information of the points are fully preserved; no changes are made to the planar position or color assignments, only the elevation dimension is updated. For points at the boundary of features, the replacement is completed after weighted interpolation calculations based on the elevation values ​​of multiple adjacent grids, eliminating elevation jumps caused by differences in correction rules for different land types. The entire process does not require re-performing the dense matching calculation of the stereo image; the point cloud elevation is optimized and updated solely through numerical replacement. For example, the elevation of a building edge point cloud that originally exhibited jagged distortion becomes straight and regular after replacement, and the elevation of a water surface point cloud that originally had disordered elevation becomes flat and continuous after replacement, ultimately forming an optimized color point cloud with natural planar distribution, accurate color information, and improved elevation accuracy.

[0058] Step S106: Using the basic digital line vector as boundary constraint, the constrained Deloni triangulation algorithm is adopted to construct a regular triangular mesh model based on the optimized color point cloud.

[0059] Specifically, step S106 includes: Step g1: Extract the feature outlines from the basic digital line vector data and set the feature outlines as hard constraint edges of the triangulation.

[0060] Furthermore, hard constraint edges refer to mandatory boundary lines that must be strictly followed during triangulation operations, prohibiting triangular facets from crossing or penetrating them. This ensures precise alignment between the triangular mesh and the boundaries of ground features. During execution, feature layers related to ground feature boundaries are selected from the basic digital line vector data. Closed contour lines of various core features are extracted, typically including the outer contour lines of permanent buildings, the shorelines of lakes and rivers, the edges of main roads, the boundaries of large woodlands, and the boundaries of large hardened areas. Topological preprocessing is performed on all extracted contour lines, systematically checking and eliminating topological issues such as overlapping lines, redundant nodes, self-intersections, and breakpoints. Node connection processing is applied to broken contours, deduplication and fusion processing is performed on overlapping lines, and geometric correction processing is performed on self-intersecting contours to ensure that each contour line has a continuous and closed geometric shape. After completing the topology verification, a mandatory constraint attribute mark is added to all contour lines, incorporating them into the constraint rule system of triangulation. This clarifies that the edges of triangular faces in the triangulation operation must not cross the constraint boundaries. For example, after setting the outer contour line of a building as a hard constraint edge, the triangular faces generated by the triangulation will not simultaneously cover the top surface of the building and the surrounding ground. After setting the shoreline of a water system as a hard constraint edge, the triangular faces will not simultaneously cross the water surface and the shoreline, thus ensuring the geometric regularity of the boundaries of features from a topological perspective.

[0061] Step g2: Using the optimized colored point cloud as nodes and hard-constrained edges as boundaries, perform constrained Deloni triangulation to obtain a two-dimensional triangular mesh.

[0062] Furthermore, constrained Deloni triangulation refers to a triangular mesh construction algorithm that, based on the empty circle characteristic of traditional Deloni triangulation, introduces mandatory boundary constraints while simultaneously considering the optimality of mesh shape and the fit of terrain boundaries. During execution, the optimized color point cloud is first projected onto a two-dimensional plane, forming a set of triangulation nodes composed of planar coordinate points. Simultaneously, all hard-constrained edges are incorporated into the constraint framework of the triangulation operation, and the topological skeleton of the constraint boundaries is pre-constructed, clarifying the node order and connection relationships of the constraint edges. A point-by-point insertion triangulation algorithm is used to add all nodes sequentially. Each time a node is inserted, the triangular mesh structure of the local region is updated, ensuring that the generated triangular patches always satisfy the Deloni empty circle criterion. At the same time, the spatial positional relationship between triangular edges and hard-constrained edges is strictly verified. If a triangular edge crosses the constraint boundary, the mesh structure is immediately adjusted through a local topological transformation to ensure that no triangular edge crosses the constraint boundary. For special locations such as corners of features and boundary junctions, the side lengths and angles of local triangular patches are optimized to avoid narrow and elongated triangles and deformed patches. For example, the triangular sides at building corners strictly follow the direction of the building corner, and the triangular sides at the shoreline of waterways are continuously and evenly distributed along the shoreline. The result is a two-dimensional triangular mesh with boundaries that accurately fit the outline of features and regular and uniform patch shapes.

[0063] Step g3 maps the two-dimensional triangular mesh to three-dimensional space to generate a regular triangular mesh model.

[0064] Furthermore, mapping a 2D triangular mesh to 3D space refers to assigning a corresponding elevation value and color attribute to each vertex of the planar mesh, restoring the planar topology to a 3D geometric form. This is the core step in constructing the geometric skeleton of a 3D model. During execution, all vertices of the 2D triangular mesh are traversed. Based on the planar coordinates of each vertex, the elevation values ​​and color attribute information of the corresponding points in the color point cloud are matched and optimized, expanding the 2D planar coordinates of the vertices to 3D spatial coordinates. Simultaneously, the color attributes of the corresponding vertices are associated with each triangular facet, establishing a correspondence between facets and color information. The mapping process completely preserves all topological connections formed by the 2D partitioning, without changing the node connection structure and boundary morphology of the triangular facets. Dimensionality is increased only through elevation assignment, avoiding boundary misalignment problems caused by changes in the topological structure. For the vertices at the boundaries of ground features, the elevation values ​​corresponding to the constraint boundaries are strictly used to ensure that the triangular patches on both sides of the boundary present elevation differences that conform to the characteristics of the ground features. For example, the inner vertices of the building boundary use the elevation of the building's top surface, while the outer vertices use the ground elevation, resulting in a clear difference in the building facade elevation after mapping. The vertices of the water surface area uniformly use the flat elevation, resulting in a continuous and gentle water surface shape after mapping. Finally, a regular triangular mesh model with clear geometric boundaries, regular ground feature shapes, and stable topological structure is generated.

[0065] Step S107: Perform pixel-level texture bonding processing on the regular triangular mesh model based on the true orthophoto image to obtain a texture-bonded triangular mesh model.

[0066] Specifically, step S107 includes: Step h1 establishes a texture coordinate mapping relationship for the regular triangular mesh model and generates mesh texture coordinate data.

[0067] Furthermore, the texture coordinate mapping relationship refers to the positional correspondence between 3D mesh vertices and 2D image pixels, used to determine the texture sampling area corresponding to each triangular facet surface in the true orthophoto. During execution, the spatial coordinate reference of the regular triangular mesh model and the true orthophoto is first verified to confirm that their coordinate projection methods are completely consistent with the planar coordinate system. All vertices of the regular triangular mesh model are traversed vertex by vertex. Based on the planar spatial coordinates of each vertex, the corresponding pixel position of the vertex in the true orthophoto pixel coordinate system is calculated, and this pixel position is recorded as the vertex's texture coordinate. The calculation process strictly follows the principle of one-to-one spatial correspondence: the vertex of the building top surface corresponds to the pixel position of the building top surface in the true orthophoto; the vertex of the water system shoreline corresponds to the pixel position of the water system boundary in the true orthophoto; and the vertex of the road edge line corresponds to the pixel position of the road edge in the true orthophoto, ensuring that the texture coordinates at the feature boundaries are accurately aligned with the geometric boundaries. After completing the texture coordinate calculation for all vertices, the texture coordinates of three vertices are associated with each triangular facet, clarifying the linear interpolation rules for the texture within the facet, ultimately forming mesh texture coordinate data where each vertex carries the corresponding image pixel position information.

[0068] Step h2: Match the pixel coordinates and mesh texture coordinates of the true radio image, perform pixel-level texture mapping, and obtain the initial texture mesh model.

[0069] Furthermore, pixel-level texture mapping refers to the process of accurately projecting the pixel color information of a true orthophoto onto a 3D mesh surface based on texture coordinate correspondence, generating a 3D model with the texture of the real image. During execution, the regular triangular mesh model is traversed triangular facet by triangular facet. Based on the texture coordinates of the three vertices of each facet, the pixel range covered by the facet on the true orthophoto is determined. A high-precision texture sampling algorithm is used to collect the color values ​​of the image pixel by pixel, and the collected color information is mapped to the surface of the corresponding triangular facet. The sampling process maintains pixel-level matching accuracy without blurring or downsampling. Triangular facets of building rooftops are accurately mapped to the pixel colors of the corresponding building rooftops in the true orthophoto; triangular facets of water areas are accurately mapped to the image colors of the corresponding water areas; and triangular facets of forest areas are accurately mapped to the image colors of the corresponding vegetation, ensuring a complete correspondence between texture information and geometric position. For triangular facets spanning multiple true orthophoto scenes, color information is sampled from the corresponding images and allocated to the corresponding areas of the facet according to spatial proportion. After completing the texture sampling and mapping of all facets, an initial texture mesh model with complete image texture is formed.

[0070] Step h3 involves performing texture color balancing and seam smoothing on the initial texture mesh model to obtain a texture-fitting triangular mesh model.

[0071] Furthermore, texture color uniformity refers to the process of unifying the global texture color benchmark and eliminating color differences between different image data sources, while seam smoothing refers to the process of eliminating texture gaps at the joints of adjacent images and improving visual continuity. During execution, global color statistics are first performed on the initial texture mesh model to calculate the brightness, contrast, and hue distribution parameters of different image coverage areas. Based on the global color benchmark, the color parameters of each area are adjusted to make the texture areas corresponding to images from different time phases and batches more consistent. For example, for adjacent areas covered by satellite images from different orbits, brightness and hue deviations are eliminated through color parameter adjustments. After global color uniformity is completed, the joint positions of all adjacent images are located, and transition processing ranges are defined on both sides of the joint line. Color gradient interpolation is used to adjust the texture color of the triangular faces within the range, ensuring a smooth transition of texture on both sides of the joint and eliminating harsh stitching marks. For areas where linear features such as roads and waterways cross joints, the texture transition effect along the feature lines is further optimized to ensure the visual continuity of the feature textures. After completing all processing, the overall texture quality is verified, and supplementary adjustments are made to areas with uneven local color and obvious seams. Finally, a texture-fitting triangular mesh model with uniform color and continuous and natural texture is obtained.

[0072] Step S108: Perform mesh simplification and texture compression sequentially on the texture-fitting triangular mesh model to output a realistic 3D model.

[0073] Specifically, mesh simplification refers to the geometric processing operation that reduces the total number of triangular faces while preserving key geometric features, thereby optimizing the overall data volume of the model. During execution, an edge-folding mesh simplification algorithm is employed. First, a global geometric feature analysis is performed on the texture-fitting triangular mesh model to identify key feature areas such as feature boundaries, terrain transition lines, and abrupt landform changes. Different simplification thresholds are set for regions with different feature levels. For regions with significant features such as building outlines, mountain ridgelines, valley bottoms, and road edges, the simplification level is reduced, retaining sufficient triangular faces to maintain geometric accuracy. For regions with minimal terrain undulation, such as open plains, large bodies of water, and gently sloping woodlands, the simplification level is increased, merging coplanar triangular faces and eliminating redundant faces. The simplification process strictly adheres to feature boundary constraints; all simplification operations do not cross feature outline boundaries and do not alter the overall topological connectivity of the model. In densely built-up urban areas, more facets are retained to maintain building angles and facades, while in large suburban farmlands and gently sloping mountainous areas, the number of facets is appropriately reduced. This effectively reduces the overall data volume while ensuring geometric accuracy, thus simplifying the mesh. Texture compression refers to image processing operations that encode and compress texture image data while preserving texture visual quality, thereby reducing the data storage volume. During execution, the texture data of the simplified texture-fitting triangular mesh model is first divided into blocks, and differentiated compression parameters are set based on the type of terrain feature and the richness of texture details. For areas with rich texture details and high recognition requirements, such as building rooftops, road intersections, and waterfronts, a low-loss compression mode is used to preserve clear texture details and color levels. For areas with high texture homogeneity, such as large woodlands, open water surfaces, and homogeneous farmland, a higher compression ratio encoding mode is used to reduce the data volume without affecting the overall visual effect. Throughout the compression process, the correspondence between texture coordinates and mesh vertices remains unchanged to avoid texture misalignment, stretching, distortion, and other problems. The building textures in the urban core area retain high clarity and detail, while the compression ratio in the contiguous forest areas in the suburbs is appropriately increased. After completing the full-area texture compression process, the storage format and coordinate reference of the model are unified, and the final output is a real-scene 3D model with accurate geometric shape, natural texture reproduction, and optimized data volume.

[0074] This embodiment provides a method for modeling ultra-large scenes based on pre-processing of stereo satellite imagery. By advancing the elevation defect management process to before triangulation mesh construction, the pre-processing step prevents matching errors from being amplified along the modeling chain at the source. Close technical synergy exists between each step: joint regional network adjustment and epipolar resampling unify the spatial coordinate benchmarks of multiple imagery scenes, providing geometric accuracy assurance for subsequent dense matching; the initial digital surface model generated by dense matching forms a pixel-level linkage with the stereo epipolar imagery, enabling defect identification to utilize both elevation gradient features and image texture features simultaneously; batch correction based on geometric constraints for different land features sets differentiated correction rules for the morphological characteristics of four types of land features: buildings, water bodies, woodlands, and shadows, achieving refined elevation management for all land features; the corrected elevation data updates the original point cloud in reverse, ensuring that the point cloud elevation accuracy is consistent with the corrected elevation system, providing high-quality node data for triangulation mesh construction; constrained Deloni triangulation uses land feature outlines as hard constraint boundaries, ensuring the regularity of the mesh geometric boundaries; pixel-level texture bonding and subsequent simplification and compression ensure the visual realism and data universality of the model. The above-mentioned steps are interconnected and together enable high-precision, automated, and mass production of ultra-large-scale real-scene 3D models, effectively solving the technical problems of existing automated elevation optimization methods lacking multi-view geometric constraints and having limited accuracy improvement.

[0075] The above are embodiments of the ultra-large scene modeling method based on pre-processing of stereo satellite imagery provided in this application. Other embodiments of ultra-large scene modeling based on pre-processing of stereo satellite imagery provided in this application are described below.

[0076] This invention also discloses a large-scale scene modeling system based on pre-processing of stereo satellite imagery, such as... Figure 2 As shown, the system includes: The acquisition module is used to acquire true orthophoto data, raw dense matching point cloud, stereo satellite image pair data, basic digital line vector data, orthophotos, and space radar elevation reference data of the target area.

[0077] The sampling module is used to perform joint regional network adjustment and epipolar resampling processing on the stereo satellite image pairs based on the orthophoto and the aerospace radar elevation reference data to obtain a standardized stereo epipolar image.

[0078] The generation module is used to generate global initial digital surface model raster data based on the standardized stereo epipolar image using a dense matching algorithm.

[0079] The correction module is used to input the standardized stereo epipolar image, the global initial digital surface model raster data, and the basic digital line vector data into a pre-built digital surface model and stereo image pixel-level linkage editing platform, perform defect classification and identification processing and batch correction processing of geometric constraints of ground features, and output global corrected digital surface model raster data.

[0080] The optimization module is used to update the elevation dimension data of the original dense matching point cloud based on the global corrected digital surface model raster data, and generate an optimized color point cloud.

[0081] The construction module is used to construct a regular triangular mesh model based on the optimized color point cloud, using the basic digital line vector as boundary constraints and employing the constrained Deloni triangulation algorithm.

[0082] The bonding module is used to perform pixel-level texture bonding processing on the regular triangular mesh model based on the true orthophoto image to obtain a texture-bonded triangular mesh model.

[0083] The output module is used to sequentially perform mesh simplification and texture compression on the texture-fitting triangular mesh model to output a realistic 3D model.

[0084] The technical solution of this invention has the following advantages: This embodiment provides a large-scale scene modeling system based on pre-processing of stereo satellite imagery. The functional modules form an ordered processing chain according to the data flow: the acquisition module uniformly accesses multi-source geographic data, providing standardized input for the entire process; the sampling module completes the adjustment of stereo image pairs and epipolar resampling to establish a unified spatial geometric benchmark; the generation module performs dense matching based on epipolar imagery to output an initial elevation raster; the correction module relies on a pixel-level linked editing platform to complete defect identification and batch correction of features, achieving pre-processing of elevation defects; the optimization module updates the corrected elevation data back to the original point cloud, generating a high-precision optimized color point cloud; the construction module performs triangulation with vector boundaries as constraints to form a regular 3D mesh skeleton; the bonding module completes pixel-level texture mapping and color balancing, outputting a mesh model with realistic textures; and the output module performs simplification and compression to generate the final standardized result. The modules are connected through standardized data interfaces, supporting parallel processing and pipeline operations, together forming a complete large-scale scene real-world 3D modeling technology system, realizing automation of the entire modeling process and pre-processing of elevation accuracy.

[0085] In this embodiment, the ultra-large scene modeling system based on stereo satellite image pre-processing is presented in the form of functional units. Here, a unit refers to an ASIC (Application Specific Integrated Circuit), a processor and memory that execute one or more software or fixed programs, and / or other devices that can provide the above functions.

Claims

1. A method for modeling ultra-large scenes based on pre-processing of stereo satellite imagery, characterized in that, The method includes: Acquire true orthophoto data, raw dense matching point cloud, stereo satellite image pair data, basic digital line vector data, orthophotos, and space radar elevation reference data for the target area; Based on the orthophoto and the aerospace radar elevation reference data, the stereo satellite image pair data is subjected to joint regional network adjustment and epipolar resampling to obtain a standardized stereo epipolar image. Using a dense matching algorithm, global initial digital surface model raster data is generated based on the standardized stereo epipolar image. The standardized stereo epipolar image, the global initial digital surface model raster data, and the basic digital line vector data are input into a pre-built digital surface model and stereo image pixel-level linkage editing platform to perform defect classification and identification processing and batch correction processing of geometric constraints of ground features, and output global corrected digital surface model raster data. Based on the global corrected digital surface model raster data, the elevation dimension data of the original dense matching point cloud is updated to generate an optimized color point cloud; Using the basic digital line vector as boundary constraints, a regular triangular mesh model is constructed based on the optimized color point cloud using the constrained Deloni triangulation algorithm. Based on the true orthogonal image, pixel-level texture bonding processing is performed on the regular triangular mesh model to obtain a texture-bonded triangular mesh model. The texture-fitting triangular mesh model is sequentially subjected to mesh simplification and texture compression to output a realistic 3D model.

2. The method according to claim 1, characterized in that, The process of performing joint regional network adjustment and epipolar resampling on the stereo satellite image pairs based on the orthophoto and the spaceborne radar elevation reference data to obtain a standardized stereo epipolar image includes: Extract planar control points from the orthophoto image, and extract elevation control points from the aerospace radar elevation reference data; An adjustment control datum is constructed based on the plane control points and the elevation control points; Based on the adjustment control benchmark, joint regional network adjustment is performed on the stereo satellite image pairs to identify and retain qualified stereo satellite image pairs; Perform epipolar resampling on the adjusted stereo satellite image pairs to generate standardized stereo epipolar images.

3. The method according to claim 2, characterized in that, The process of generating global initial digital surface model raster data based on the standardized stereo epipolar image using a dense matching algorithm includes: A multi-view dense matching algorithm is used to perform pixel-by-pixel matching on the standardized stereo epipolar image to generate a dense point cloud with the same name. Statistical filtering and morphological filtering are performed on the dense point cloud to obtain a denoised point cloud. Rasterization interpolation is performed on the denoised point cloud of the same name to generate global initial digital surface model raster data. The grid spacing of the global initial digital surface model raster data is in the range of [0.5m, 2m].

4. The method according to claim 3, characterized in that, The process involves inputting the standardized stereo epipolar image, the initial global digital surface model raster data, and the basic digital line vector data into a pre-built digital surface model and stereo image pixel-level linkage editing platform. Defect classification and identification processing, as well as batch correction processing based on geometric constraints of different features, are performed to output globally corrected digital surface model raster data, including: Spatial geometric registration is performed on the standardized stereo epipolar image, the global initial digital surface model raster data, and the basic digital line vector data to obtain a linked dataset; Defect classification and identification processing is performed on the linked dataset to generate elevation anomaly area data with land cover category labels; The data of the elevation anomaly area are subjected to batch correction processing based on geometric constraints of ground features to generate raster data of the global corrected digital surface model.

5. The method according to claim 4, characterized in that, The defect classification and identification process based on the linked dataset generates elevation anomaly area data with land cover category labels, including: An edge detection algorithm and a deep learning semantic segmentation model are used to perform anomaly detection on the linked dataset to identify elevation anomaly regions, which include jagged edges and undulating top surfaces of buildings, abrupt elevation changes in water bodies, uneven areas of forest canopies, and elevation voids caused by shadows. The elevation anomaly area is reviewed and defect supplemented to obtain the reviewed elevation anomaly area; Based on the land feature attributes, the verified elevation anomaly areas are classified and labeled to generate elevation anomaly area data with attached land feature category labels. The land feature category labels include building category labels, water body category labels, forest land category labels, and shadow occlusion category labels.

6. The method according to claim 4, characterized in that, The step of performing batch correction processing based on the elevation anomaly area data to generate global corrected digital surface model raster data includes: For areas with abnormal elevations and attached building category markers, edge straightening constraints and top surface elevation smoothing are applied. For areas with abnormal elevations marked by water body categories, implement elevation unification for static water areas and longitudinal slope constraint for flowing water areas. For areas with abnormal elevations and attached forest land category markers, canopy morphology filtering is applied. For elevation anomaly areas marked with attached shadow occlusion categories, multi-view elevation information from stereo epipolar imagery is reused to fill elevation voids; Generate raster data for a globally corrected digital surface model.

7. The method according to claim 6, characterized in that, The step of updating the elevation dimension data of the original densely matched point cloud based on the global corrected digital surface model raster data to generate an optimized color point cloud includes: Extract the planar coordinate information and color attribute information of the original densely matched point cloud to generate basic attribute data of the point cloud; Establish the planar coordinate mapping relationship between the global corrected digital surface model raster and the point cloud basic attribute data; Based on the aforementioned planar coordinate mapping relationship, the elevation values ​​of the original densely matched point cloud are replaced point by point to generate an optimized colored point cloud.

8. The method according to claim 7, characterized in that, The process of constructing a regular triangular mesh model based on the optimized color point cloud, using the basic digital line vector as boundary constraints and employing a constrained Deloni triangulation algorithm, includes: Extract the feature outlines from the basic digital line vector data, and set the feature outlines as hard constraint edges of triangulation; Using the optimized colored point cloud as nodes and the hard-constrained edges as boundaries, constrained Deloni triangulation is performed to obtain a two-dimensional triangular mesh; The two-dimensional triangular mesh is mapped to three-dimensional space to generate a regular triangular mesh model.

9. The method according to claim 8, characterized in that, The step of performing pixel-level texture fitting processing on the regular triangular mesh model based on the true orthophoto image to obtain a texture-fitted triangular mesh model includes: Establish texture coordinate mapping relationships for the regular triangular mesh model and generate mesh texture coordinate data; Match the pixel coordinates of the true radiometric image with the mesh texture coordinate data, perform pixel-level texture mapping, and obtain the initial texture mesh model; The initial texture mesh model is subjected to texture color uniformation and seam smoothing processing to obtain a texture-fitting triangular mesh model.

10. A large-scale scene modeling system based on pre-processing of stereo satellite imagery, characterized in that, The system includes: The acquisition module is used to acquire true orthophoto data of the target area, raw densely matched point cloud, stereo satellite image pair data, basic digital line vector data, orthophoto and aerospace radar elevation reference data; The sampling module is used to perform joint regional network adjustment and epipolar resampling on the stereo satellite image pair data based on the orthophoto and the aerospace radar elevation reference data to obtain a standardized stereo epipolar image. The generation module is used to generate global initial digital surface model raster data based on the standardized stereo epipolar image using a dense matching algorithm. The correction module is used to input the standardized stereo epipolar image, the global initial digital surface model raster data, and the basic digital line vector data into a pre-built digital surface model and stereo image pixel-level linkage editing platform, perform defect classification and identification processing and batch correction processing of geometric constraints of ground features, and output global corrected digital surface model raster data. The optimization module is used to update the elevation dimension data of the original dense matching point cloud based on the global corrected digital surface model raster data, and generate an optimized color point cloud; The construction module is used to construct a regular triangular mesh model based on the optimized color point cloud, using the basic digital line vector as boundary constraints and employing the constrained Deloni triangulation algorithm. The bonding module is used to perform pixel-level texture bonding processing on the regular triangular mesh model based on the true orthophoto image to obtain a texture-bonded triangular mesh model. The output module is used to sequentially perform mesh simplification and texture compression on the texture-fitting triangular mesh model to output a realistic 3D model.

Images