Large-scale city building fine satellite three-dimensional reconstruction method based on prior guidance iteration

Abstract

The invention discloses a large-scale city building fine satellite three-dimensional reconstruction method based on prior guidance iteration, which comprises the following steps of: acquiring three kinds of data of a to-be-constructed area, namely an image sequence of a panchromatic satellite optical image, a DEM (Digital Elevation Model) and building vector data with building floor height information, selecting one image as a reference image, and selecting other images as source images; using an angle screening method to obtain a candidate image sequence for subsequent calculation; re-sampling gray information of the reference image to an object space to obtain a gray base map, registering building vector data with floor height information with the gray base map, and fusing with a DEM (Digital Elevation Model) re-sampled to the object space to generate elevation priori; step-by-step downsampling is carried out on a candidate image sequence, an image pyramid is constructed, then an improved object space semi-global dense matching method is used for carrying out layer-by-layer iterative matching on the image pyramid from top to bottom, and a reconstructed DSM is obtained.

Description

Technical Field

[0001] This invention belongs to the technical field of 3D reconstruction, specifically relating to a method for fine satellite 3D reconstruction of large urban buildings based on prior-guided iteration. Background Technology

[0002] With the acceleration of global urbanization, the rapidly growing urban population poses challenges to urban design and planning. 3D models generated from the 3D reconstruction of urban areas provide accurate spatial data support for refined urban management, holding significant value in modern urban governance and planning. Among these, automatic 3D reconstruction methods based on satellite remote sensing optical imagery are widely used to generate urban 3D models to aid urban management due to the inherent advantages of low acquisition cost and wide coverage of optical imagery, attracting widespread attention from the photogrammetry and remote sensing communities. Although the development of traditional reconstruction methods in the field of satellite 3D reconstruction has made it possible to automatically produce geographic 3D products such as digital surface models (DSMs) based on stereo satellite optical imagery, complex urban scenes, due to the presence of numerous high-rise buildings and their complex and dense building topology, generally exhibit more severe shadows, occlusions, and complex topological structures compared to natural terrain scenes. This leads to a significant decrease in reconstruction quality for traditional reconstruction methods in urban scenes.

[0003] Although existing technologies have improved upon traditional semi-global stereo matching algorithms in computer vision to address the characteristics of remote sensing imagery and constructed 3D reconstruction frameworks for satellite imagery, severe shadows and occlusions caused by densely packed, structurally complex, and topologically detailed building clusters in urban scenes, combined with the pixel mixing phenomenon commonly found in satellite imagery, lead to typical problems in 3D reconstruction results, such as building structural adhesion, elevation deviation, blurred edges, and reconstruction voids. For example, in 2016, Ghuffar et al. proposed a semi-global dense matching satellite stereo digital surface model (DSM) generation method based on object space, such as... Figure 1 As shown, the main problems are as follows: 1. Ghuffar et al.'s method employs a constant penalty term in the cost aggregation stage of object-side semi-global dense matching. This means applying a uniform smoothing constraint to any disparity variation across all pixel neighborhoods. While this fixed penalty strategy effectively suppresses noise in flat or texture-continuous regions, it lacks adaptive adjustment capabilities for image structural features. Because the fixed penalty cannot adjust the constraint strength based on image texture or gradient changes, the method's performance in urban areas is limited. It is particularly prone to over-smoothing at abrupt elevation changes such as building rooftops and facades, and buildings and ground surfaces, leading to blurred building edges, facade collapse, or structural incompleteness.

[0004] 2. Ghuffar et al.'s method used cost aggregation based on the smoothness assumption in the cost aggregation stage of object-side semi-global dense matching. However, due to the prevalent occlusion, shadows, and pixel mixing phenomena in urban areas, the initial cost value may contain significant errors. These errors accumulate gradually because they cannot be effectively identified in subsequent aggregation processes, ultimately causing the reconstruction results to fail to accurately reflect the true structural characteristics of ground features, resulting in reduced reconstruction quality of building edges and inaccurate estimation of building heights.

[0005] 3. Ghuffar et al.'s method uses the matching results of low-resolution pyramids to limit the disparity search range of high-resolution pyramids, thereby accelerating the efficiency of dense matching algorithms when performing dense matching on large-scale images. However, in practical applications in urban scenes, dense urban scenes often have characteristics such as small spacing between buildings, complex building topology, and rich details. In such cases, multi-layer downsampling of the image will result in the loss of a large amount of detail and exacerbate the uncertainty of grayscale information from mixed pixels, causing problems such as building adhesion, missing building bodies, and loss of building details, which affect the reconstruction quality. Summary of the Invention

[0006] This invention provides a method for fine satellite 3D reconstruction of large urban buildings based on prior-guided iteration, in order to solve the technical problems of fixed penalty terms, poor adaptability, and obvious errors in existing schemes when performing dense matching.

[0007] To achieve the above objectives, the present invention provides the following technical solution: A method for high-resolution satellite 3D reconstruction of large-scale urban buildings based on prior-guided iteration includes the following steps: Step 1: Obtain three types of data for the area to be constructed: panchromatic satellite optical imagery, DEM, and image sequence of building vector data with building height information. Select one image as the reference image and the others as source images. Then, use the angle filtering method to obtain candidate image sequence for subsequent calculations. Step 2: Resample the grayscale information of the reference image to the object space to obtain a grayscale base map. Then, register the building vector data with floor height information with the grayscale base map, and then fuse it with the DEM resampled to the object space to generate an elevation prior. Step 3: By downsampling the candidate image sequence step by step, an image pyramid is constructed. Then, an improved object-side semi-global dense matching method is used to iteratively match the image pyramid layer by layer from top to bottom to obtain the reconstructed DSM. The improved object-side semi-global dense matching method uses CENSUS to calculate the cost value, and then sequentially uses the building prior-guided cross-arm cost aggregation method and the grayscale base map-guided scan line cost aggregation method to perform two-stage aggregation optimization of the cost value. At the same time, when performing layer-by-layer iterative matching on the image pyramid, the elevation prior is used to initialize the elevation search range of the top layer, and the elevation search range of the remaining layers is constrained by the matching results of the upper layer to constrain the elevation search range of the lower layer.

[0008] Furthermore, when constructing the image pyramid, the source image, reference image and corresponding object space to be matched are downsampled step by step. Except for the full-size layer at the bottom, the sampling interval of each layer is performed according to the following relationship. image: , ; Object space: , , ; in, and V represents the horizontal and vertical displacement step size of the image, respectively. , Y, Z represents the displacement step size in the latitude, longitude, and height directions in the object space, respectively. Indicates the layer number; For the full-size layer (layer 0), the image is not downsampled, and the displacement step size relationship in the object space is set to... , , .

[0009] Furthermore, during the layer-by-layer iterative matching, an improved object-side semi-global dense matching method is used for the current layer of the image pyramid. Perform matching to obtain the elevation value corresponding to the horizontal coordinate v of each object point. Together they form an elevation map Then check the elevation map Morphological erosion and morphological dilation were performed to obtain elevation values. and Then the next layer of the image pyramid The search range corresponding to the horizontal coordinate v is expressed as ),in, This represents the elevation displacement step size. The step size method coefficients are... For full-size layers, the search range is represented as follows: ),in, This is the elevation value after matching the previous layer adjacent to the full-size layer.

[0010] Furthermore, the search range at the top of the image pyramid is based on the elevation prior definition, which includes ground elevation. Building outline buffer zone elevation and building elevation within the building outline Three types, Ground elevation minimum elevation and highest elevation The definition is as follows: in, Elevation values ​​from reference DEM imagery; Building outline buffer zone elevation minimum elevation and highest elevation The definition is as follows: in, Based on reference building height; Building elevation within the building outline buffer zone minimum elevation and highest elevation The definition is as follows: Furthermore, for the improved object-space semi-global dense matching method, the following formula is used to calculate the cost. , in, This represents the transformation value corresponding to the grayscale value within the projection window on the reference image. This represents the transformation value corresponding to the grayscale value within the projection window on the i-th source image, where N is the total number of source images. This represents the cost calculated for the i-th source image and the reference image. The weighted average of all cost values ​​is the final cost value, where u and v represent the vertical and horizontal coordinates in the image, and X, Y, and Z represent the coordinates of any point in the object space in the longitude, latitude, and elevation directions, respectively. Using cross-arm cost aggregation, given a point V in the object space. When it finds the endpoint When the arm stops in a certain direction, the stopping rule at that point is as follows: in, This represents the color difference between two object points; the corresponding grayscale information comes from the grayscale base image. and Indicates the grayscale threshold. This represents stricter grayscale constraints. Represents the spatial distance between two object points. and Indicates the distance threshold. This represents the previous object point, taking the left arm as an example. Indicates the roof area. Indicates the ground area; The calculation formula for the scan line cost aggregation is as follows: in, Let V represent the coordinates of a point V in object space. Equivalent to the horizontal coordinates (X, Y) of the object point, where Z represents the elevation coordinates of the object point. This represents the aggregated cost of the scan line at the assumed elevation Z for point V. Indicates the object point after the cross arm cross cost aggregation. Given the cost at elevation Z, r defines the convergence direction. and This represents the penalty term, which satisfies... , express The initial values ​​correspond to small and large parallax changes, respectively. and Representing the object point in object space respectively and The horizontal coordinates of the grayscale base image The corresponding grayscale value in the image.

[0011] Furthermore, when obtaining the candidate image sequence using the angle screening method in step one, the image with the smallest solar zenith angle and the imaging angle closest to the vertical viewing angle is selected from the image sequence as the reference image, and the remaining images are the source images. Then, the cross angle and solar angle difference between the reference image and each source image are calculated one by one. Finally, the source images with cross angle and solar angle differences greater than a set threshold and the reference image are selected to form the candidate image sequence.

[0012] Furthermore, the intersection angle between the reference image and any source image is calculated using the following formula. Difference with the sun angle , in, and These represent the angles of the principal point vectors of the reference image and the source image, respectively. and These represent the solar zenith angle and solar azimuth angle of the reference image, respectively. and These represent the solar zenith angle and solar azimuth angle of the source image, respectively.

[0013] Furthermore, in step two, the building vector data with floor height information is first divided into multiple small blocks, and then local affine or homography transformation is used to register each small block with the corresponding area of ​​the grayscale base map. Then, the registration result is rasterized and resampled together with the DEM into the object space for superposition and fusion to generate elevation prior.

[0014] Compared with the prior art, the beneficial effects of the present invention are: 1. A two-stage cost aggregation strategy is introduced. Before the traditional scanline cost aggregation, a cross-arm cross aggregation is added, which enables the construction of cross-shaped support regions in the local support region, achieving local cost smoothing with edge preservation and effectively avoiding false propagation caused by crossing strong edges. Then, the second-stage scanline cost aggregation is carried out on this basis, using global smoothness constraints to supplement the energy consistency across regions.

[0015] 2. Using building outline priors to guide the first stage of cross-arm aggregation, the building roof outline data registered with the grayscale base map is sampled into the object space as a terrain prior and threshold, which together with the cross-arm range constraint to reduce the gradual accumulation of errors during cost aggregation and improve the reconstruction quality of building edges and height.

[0016] 3. The resampled orthophoto is used as a grayscale base map to dynamically guide the second-stage scan line aggregation penalty term. By resampling the orthophoto to the object space, a dynamic correlation between pixel intensity gradient and penalty term is established, which solves the problem of insufficient adaptability of fixed penalty parameters at building edges and improves the reconstruction quality of building edges.

[0017] 4. Using building height and DEM-constructed coarse terrain elevation constraints as the highest-level pyramid elevation range constraints to guide semi-global dense matching of objects, the problem of building loss and adhesion in the three-dimensional matching reconstruction results of urban areas is effectively improved, and the outlier noise of building edges in the reconstruction results is effectively suppressed by the range limitation. Attached Figure Description

[0018] Figure 1 This is a schematic diagram illustrating the method of Ghuffar et al. in the prior art. Figure 2 This is a schematic diagram of the overall framework of the present invention; Figure 3 This is a schematic diagram of the terrain prior initialization elevation constraint of the present invention; Figure 4 This is a schematic diagram illustrating the main stages of the improved object-side semi-global dense matching method of the present invention. Figure 5 This is a distribution diagram of the experimental area in a specific embodiment of the present invention; Figure 6 This is a comparative schematic diagram of DSM images obtained by reconstructing region A using various methods in a specific embodiment of the present invention. In this diagram, (A) is an orthophoto satellite image, (A-1) is the ground truth of GT, (A-2) is the ASP (MGM) result, (A-3) is the S2P result, (A-4) is the S2P (M) multi-view fusion result, and (A-5) is the result of the method of the present invention. Figure 7 This is a comparative schematic diagram of DSM images obtained by reconstructing region B using various methods in a specific embodiment of the present invention. In this diagram, (B) is an orthophoto satellite image, (B-1) is the ground truth of GT, (B-2) is the ASP (MGM) result, (B-3) is the S2P result, (B-4) is the S2P (M) multi-view fusion result, and (B-5) is the result of the method of the present invention. Figure 8 This is a comparative schematic diagram of DSM images obtained by reconstructing region C using various methods in a specific embodiment of the present invention. In this diagram, (C) is an orthophoto satellite image, (C-1) is the ground truth GT value, (C-2) is the ASP (MGM) result, (C-3) is the S2P result, (C-4) is the S2P (M) multi-view fusion result, and (C-5) is the result of the method of the present invention. Figure 9This is a comparative schematic diagram of DSM images obtained by reconstructing region D using various methods in a specific embodiment of the present invention. In this diagram, (D) is an orthophoto satellite image, (D-1) is the ground truth of GT, (D-2) is the ASP (MGM) result, (D-3) is the S2P result, (D-4) is the S2P (M) multi-view fusion result, and (D-5) is the result of the method of the present invention. Figure 10 This is a comparative schematic diagram of DSM images obtained by reconstructing region E using various methods in a specific embodiment of the present invention. In this diagram, (E) is an orthophoto satellite image, (E-1) is the ground truth of GT, (E-2) is the ASP (MGM) result, (E-3) is the S2P result, (E-4) is the S2P (M) multi-view fusion result, and (E-5) is the result of the method of the present invention. Figure 11 is a comparative schematic diagram of DSM images obtained by reconstructing the region AE using various methods in a specific embodiment of the present invention. Detailed Implementation

[0019] To make the technical means, creative features, objectives and effects of this invention easier to understand, the following embodiments, in conjunction with the accompanying drawings, specifically illustrate the method for fine satellite 3D reconstruction of large-scale urban buildings based on prior guided iteration. It should be noted that the description of these embodiments is for the purpose of helping to understand this invention, but does not constitute a limitation of this invention.

[0020] Based on the semi-global dense matching satellite stereo digital surface model (DSM) generation method framework disclosed by Ghuffar et al., this invention proposes a high-resolution satellite 3D reconstruction method for large-scale urban buildings based on prior-guided iteration, such as... Figure 2As shown, this method consists of three main parts. The first part preprocesses the input multi-view images, mainly including the screening of multi-view images and the registration of prior information. Through preprocessing, we obtain a high-quality candidate image sequence, which is used for subsequent dense matching of images with grayscale base maps, building outline priors, and terrain elevation priors. The second part proposes an object-side semi-global dense matching method based on urban features. This method introduces a two-stage cost aggregation strategy. In the first stage, cross-arm aggregation is performed under the guidance of building outline priors to maintain local structural consistency. In the second stage, scanline cost aggregation is performed under the dynamic guidance of grayscale base maps, and an adaptive penalty term is set, thereby fundamentally improving the accuracy and stability of dense matching in urban areas. The third part proposes a pyramid iterative matching method based on terrain prior constraints to accelerate the object-side semi-global dense matching method proposed in the second part. By introducing building terrain information, improving the transmission method of prior results between pyramid layers, and introducing a full-size matching layer, the limitations of traditional pyramid iterative matching in complex urban areas are effectively alleviated, and the integrity and accuracy of 3D reconstruction are further improved.

[0021] Specifically as follows: 1. Image preprocessing 1.1. Screening of multi-view images This invention utilizes three types of data: panchromatic satellite optical imagery, DEM imagery, and building vector data with building height information. The panchromatic satellite optical imagery consists of 32 images taken by the three-line array cameras of the ZY3-01 and ZY3-03 satellites. Each image is accompanied by a corresponding rational polynomial coefficient (RPC) file for establishing image geometric positioning. The three-line array cameras provide a forward-looking angle of +22°, a frontal angle of 0°, and a back-looking angle of -22°. The ground resolution of the frontal camera on the ZY3-01 satellite is 2.1 meters, and the ground resolution of the frontal and back-looking cameras is 3.5 meters. The ground resolution of the frontal and back-looking cameras on the ZY3-03 satellite is further upgraded to 2.5 meters compared to the ZY3-01 satellite. The panchromatic satellite optical imagery was captured between February 2020 and July 2025, with an average cloud cover of less than 20%, covering the entire Hong Kong Special Administrative Region and parts of Guangdong Province, China.

[0022] The DEM imagery, with a resolution of 5 meters, was provided by the Hong Kong Lands Department and comprehensively covers the entire Hong Kong Special Administrative Region. The building vector data, including floor height information, comes from the Hong Kong Lands Department's Hong Kong regional building vector data provided by the Hong Kong Spatial Data Platform. All permanent buildings or structures in this data are displayed as polygons, and building attributes include building height, building type, and building name. Detailed data is shown in Table 1.

[0023] Table 1 Image Data Information and Sources Satellite imagery acquisition is constrained by orbital mechanics, exhibiting two typical observation modes: Imagery from the same orbit has short temporal phase (intervals within minutes) and a fixed cross-angle; imagery from different orbits may have significant temporal phase differences (from days to years) and highly variable cross-angles. These temporal phase differences and cross-angle angles affect the quality of dense matching. In urban areas, considering the mutual occlusion between buildings, relatively small cross-angles typically produce higher-quality DSMs. Simultaneously, solar angle differences play a crucial role in urban reconstruction, as they can project different directional illuminations onto anisotropic / mirror surfaces, thus affecting the degree and extent of shadow occlusion in the image. Therefore, to avoid these factors in multi-view imagery interfering with our reconstruction results, we considered cross-angle and solar angle differences to screen multi-view imagery. Specifically, we adopted the method of Qin et al., which comprehensively considered cross-angle and solar angle differences to score the images, selecting images with scores exceeding a set threshold to form a candidate image sequence for matching.

[0024] Since the solar zenith angle and solar azimuth angle both affect the solar angle difference, before screening, the image with the smallest solar zenith angle and the imaging angle closest to the vertical viewing angle is selected as the reference image from the image sequence. All other images are used as source images. Considering that the resource series images used in our experimental data include downward-looking viewing angles that provide good vertical observation conditions, the downward-looking image with the best observation conditions is used as the reference image for the entire matching stage. The intersection angle and solar angle difference between the reference image and each source image are calculated one by one using the following formula: in, and These are the principal point vector angles of the reference image and the source image, respectively. and These are the solar zenith angle and solar azimuth angle from the reference image, respectively. and These are the solar zenith angle and solar azimuth angle of the source image, respectively.

[0025] 1.2. Registration of Prior Information First, we establish an accurate representation of ground features in the object space. We utilize the grayscale information of the reference image, which reflects the reflectivity or radiance of the ground surface to the sensor's sensitive wavelength bands. Compared to oblique imagery, using a downward-looking image as the reference image provides a more stable geometric reference in terms of occlusion and projection distortion. However, due to terrain undulations and lens imaging distortion, some geometric offsets still appear in the image. Therefore, we use orthorectification to eliminate the tilt displacement of most tall buildings, thereby ensuring the accurate positioning of planar features. The corrected reference image contains geographic coordinate information and can be used to resample to the object space to generate a grayscale base map. Next, we register the building vector data with floor height information with the grayscale base map, and simultaneously resample the DEM to the object space to fuse and generate an elevation prior.

[0026] The registration process is based on affine transformation, estimating the building's vector data using feature points between the building vector data and the grayscale base image. However, off-axis imaging and the elevation differences of buildings in urban areas cause building tilt (non-planar distortion) in the image. Off-axis images reveal roofs and walls, with buildings tilting outwards and occluding elements, resulting in displacements of different directions and magnitudes at different heights. A single affine transformation cannot simultaneously satisfy this across the entire image. To mitigate this issue, we divide the building vector data into small blocks and register each block with the corresponding image block in the grayscale base image. The registered building vectors are then rasterized based on their contained building outline and floor height information and resampled into the object space to obtain the building roof footprint map and floor height map. Local affine or homography transformations are used during registration to partially improve the registration error; for local offsets that cannot be eliminated, fine-tuning is performed through manual registration. Finally, the ground elevation information from the DEM and the building elevation information from the floor height map (estimated by multiplying the floor height by the average floor height by 3) are superimposed to obtain the elevation prior. Although we employed a series of methods during the registration process to maintain spatial consistency between ground feature information and the real world, some local biases inevitably remained and could not be completely corrected. It is worth noting that subsequent matching methods incorporated this prior information as a non-strict soft constraint in stereo matching; therefore, the impact of a small amount of residual error on the final reconstruction result is limited.

[0027] 2. Improved object-space semi-global dense matching The improved object-space semi-global dense matching is consistent with the object-space semi-global matching method published by Ghuffar et al. in terms of overall process. First, an object space O(X,Y,Z) needs to be constructed. The space enclosed by X and Y in the object space represents the geographical range of the area to be reconstructed, and Z represents the elevation range of the area to be reconstructed. An object point in the object space is projected onto the reference image and the source image after being corrected by the regional network stereo adjustment. The gray value of the projected position is used for comparison to calculate the cost. Second, due to the widespread noise and interference in the image, the initial cost calculated in this process cannot be directly used to calculate the optimal elevation. It is necessary to optimize the cost through subsequent cost aggregation. Our improvement mainly focuses on the cost aggregation stage. Cost aggregation is divided into two stages: cross-arm cost aggregation and scanline cost aggregation. The building roof outline information registered to the object space and the orthophoto are used as the building roof footprint map and grayscale base map, respectively, to guide the two-stage cost aggregation. Finally, the optimal elevation value is selected using the WTA algorithm to obtain the elevation map Z(v), and a DSM image is generated by adding coordinates, datum and other information. Figure 4 The main process of the improved object-side semi-global dense matching method is demonstrated.

[0028] 2.1 Cost Calculation In this invention, we employ window-based SCENSUS cost as the cost estimate. The matching window is defined within the object space, centered on the object block (X, Y, Z) and parallel to the reference plane of the object space. During projection, the window is reprojected onto both the reference and source images. The grayscale information within the window is used to calculate the cost, as follows: Figure 4 As shown in (a), CENSUS exhibits good radiometric robustness and outperforms traditional intensity-based matching methods when handling illumination variations, making it widely used for cost calculation. However, CENSUS cost calculation primarily targets paired windows, which may lead to multiple windows (greater than 2) in multi-view matching. We extend this to the multi-window case and propose the SCENSUS cost. The SCENSUS cost calculates the Hamming distance between the transformed values ​​of the reference image and the transformed values ​​of the source image, and then takes a weighted average of all Hamming distance results. The cost calculation formula is defined as follows: in, This represents the transformation value corresponding to the grayscale value within the projection window on the reference image. This represents the transformation value corresponding to the grayscale value within the projection window on the i-th source image, where N is the total number of source images. This represents the cost calculated for the i-th source image and the reference image. The weighted average of all cost values ​​represents the final cost. This projection method gives the algorithm stronger rotation and scale invariance compared to traditional semi-global matching methods. Furthermore, the reprojection of the matching window during the reprojection process can break through the limitations of traditional integer pixels and generate sub-pixel level coordinates in the image. This characteristic provides an important foundation for obtaining high-precision 3D results in subsequent steps.

[0029] 2.2 Cost Aggregation of Crossarms Based on Building Awareness Mei et al. proposed a cross-arm aggregation method to support matching costs in the support region, thereby reducing matching ambiguity and noise in the initial cost. The cross-arm aggregation is based on a simple but effective assumption: adjacent pixels with similar colors should have similar differences. Although this algorithm was tested on the Middlebury Benchmark dataset and achieved better matching results than direct cost aggregation, when applied to 3D reconstruction of satellite imagery, the large GSD of satellite imagery makes it more susceptible to noise and mixed pixels, making it difficult to achieve accurate range constraints using only a simple threshold. Therefore, considering the characteristics of urban areas, this invention samples building roof contour data registered with the grayscale base map into the object space as a terrain prior, which, along with the threshold, acts on the cross-arm range constraint. The process is as follows: Figure 4 As shown in (b). The study uses cross-arm aggregation for both the roof and ground areas. The outline range, as the constraint of the highest floor, effectively avoids the cost of different floor heights being included, which could lead to incorrect floor height estimation. We use cross-arm aggregation in the object space. Given an object point V in the object space, when it finds its endpoint... When the arm stops in a certain direction, the rule for stopping at that point is modified as follows: in, This represents the color difference between two points, and the corresponding grayscale information comes from the grayscale base image. and This is the grayscale threshold. This represents stricter grayscale constraints; It represents the spatial distance between two points. and Distance threshold; For object point The previous object point (taking the left arm as an example); For the roof area, The area is the ground.

[0030] 2.3 Grayscale background-guided scanline cost aggregation Building upon cross-arm aggregation, we employ scanline aggregation. To further alleviate matching ambiguity and ensure energy consistency across regions, an optimizer with both smooth constraints and moderate parallelism is needed. The process is as follows: Figure 4 As shown in (c). This invention employs the same method as OSGM and the semi-global constraint MVLL method, referencing the cost aggregation method of semi-global matching, to decompose the global smoothness constraint into a one-dimensional optimization problem involving multiple paths. However, both OSGM and the semi-global constraint MVLL method use fixed penalty values ​​for edge penalties. This can lead to either oversensitivity or oversmoothing in application scenarios such as urban building reconstruction where elevation changes abruptly. To address the needs of the scenario, we sample orthophotos with geographic information into the object space as a base map. The grayscale values ​​of this base map provide effective edge change information for dynamically adjusting the penalty value, ensuring the robustness of the algorithm in discontinuous scenarios. The cost aggregation formula is as follows: in, Let V represent the coordinates of a point V in object space, which is equivalent to (X,Y,Z) mentioned earlier. This represents the aggregated cost of the scan line at the assumed elevation Z for point V. The object point after cross-arm intersection cost aggregation Given the cost at elevation Z, r defines the convergence direction. and This represents the penalty term, which satisfies... , for The initial values ​​correspond to small and large parallax changes, respectively. and Let these represent the object points with horizontal coordinates in object space. )and grayscale base image The corresponding grayscale value in the image.

[0031] In this study, we consider aggregating the cost of implementation from eight directions. Selecting eight directions can improve the ability to capture complex textures and edge continuity in images while maintaining operational efficiency. The cost aggregation function is defined as follows: 3. Pyramid Iterative Matching Method Based on Terrain Prior Constraints Dense matching, as the most computationally intensive step in traditional satellite 3D reconstruction, faces severe computational and stability challenges in large-format imagery and densely built-up urban scenes due to the massive number of pixels, wide search range, and mismatches caused by weak textures and occlusion. Therefore, reducing computational complexity while maintaining accuracy is crucial for improving the efficiency and reliability of large-scale urban 3D reconstruction. To address these challenges, pyramid iterative matching significantly reduces computational load while maintaining accuracy through multi-scale matching from coarse to fine. The coarse-resolution layer provides stable initial disparity, while the high-resolution layer completes fine matching within a constrained search range, thus becoming a commonly used acceleration strategy in object-side semi-global dense matching frameworks. However, directly applying the traditional pyramid method to urban scenes still faces significant limitations. Urban scenes are characterized by high building density, large differences in building scale, and rich structural details. Oversampling can lead to the complete erasure of small-volume buildings, reduced building spacing, and even adhesion, resulting in the destruction of building outlines. Furthermore, compared to mountainous terrain, urban elevation exhibits discontinuous abrupt changes, and its elevation column changes often lack approximation at the distance scale. Therefore, once the coarse layer suffers from the loss of architectural details due to downsampling, the error will accumulate in subsequent layers and constrain the elevation range within the error interval, seriously affecting the accuracy of the final reconstruction result. This invention makes three improvements to the traditional pyramid framework for urban building scenarios: 1) It uses DEM images and building outlines to construct an elevation prior to initialize the coarse layer parallax, reducing the number of pyramid layers and enhancing matching stability, suppressing anomalies caused by occlusion, shadows, and mixed pixels; 2) It introduces noise culling and filling of voids in the building structure after the coarse layer parallax to ensure the structural consistency transmitted downwards and avoid error amplification; 3) It adds a full-resolution matching layer to refine the elevation interval and improve matching accuracy.

[0032] 3.1. Constructing an image pyramid When constructing the image pyramid, it is necessary to downsample the source image, reference image, and corresponding object space to be matched. Except for the full-size layer at the bottom, the sampling intervals of the other layers maintain the following relationship: image , ; Object Space , , , in, and V represents the horizontal and vertical displacement steps of the image, respectively. , Y, Z represents the displacement step size in the latitude and longitude directions and the height direction in the object space, respectively. Indicates the layer number; For full-size layers, the image is not downsampled, and the displacement step size in the object space is set to be consistent with the horizontal direction of the upper layer, but twice that in the elevation direction. , , .

[0033] 3.2 Interlayer transfer First, the search range for the top layer of the image pyramid is initialized using elevation priors constructed from DEM and building floor height maps. Elevation priors based on finite-resolution DEMs and estimated floor height elevations are insufficiently accurate. To avoid excessive interference from elevation priors in the matching process, a dual optimization strategy is adopted. On the one hand, relatively broad offset values ​​are set for the elevation priors of the ground and roof, appropriately relaxing their constraints and reducing their coercive influence on the matching results. On the other hand, a buffer is constructed in the building edge region, and a wider elevation search range is used to ensure the autonomy of matching in the building edge region. The buffer is obtained using the XOR operation result after building contour dilation and erosion, as follows: Figure 3 As shown.

[0034] We set the ground elevation as a priori elevation. Building outline buffer zone elevation and building elevation within the building outline Three types, including ground elevation. minimum elevation and highest elevation The definition is as follows: in, The initial elevation range of the ground is constructed by adding an elevation offset to the DEM elevation values ​​for non-building areas.

[0035] Building buffer zone elevation minimum elevation and highest elevation The definition is as follows: in, Based on the reference building height, the buffer has the largest elevation search range to avoid excessive interference from contour data matching, ensuring that the edges of the reconstructed building come from the aggregated results.

[0036] Building elevation within the building buffer zone minimum elevation and highest elevation The definition is as follows: Similar to the definition of ground elevation, the initial elevation range within the building buffer is constructed by adding an elevation offset to the building height.

[0037] Stereo matching is performed at each layer, with each run based on image and object space at different scales within the layer, executed sequentially from coarse to fine. Except for the highest and lowest layers (full-size layers), inter-layer results are passed as follows: An improved semi-global dense object-space matching method is used to match the current layer of the image pyramid. Perform matching to obtain the elevation value corresponding to the horizontal coordinate v of each object point. Together, they form an elevation map Z(v). Then, morphological erosion and morphological dilation are performed on the elevation map Z(v) to obtain the elevation values. and Then the next adjacent layer The search range corresponding to the horizontal coordinate v can be expressed as: ),in, This represents the elevation displacement step size. The step size method coefficients.

[0038] During the transfer process, the matching results from the upper layer are used to construct the initial search values ​​and elevation constraints for the lower layer. If significant noise exists in the upper layer, the search range may deviate from the actual terrain, preventing the lower layer from obtaining the correct solution. Such errors accumulate between layers, affecting the overall reconstruction accuracy. Although semi-global matching has some robustness in urban scenes, errors caused by occlusion, shadows, and mixed pixels are more prominent at coarse resolutions. Downsampling may introduce new mixed pixels or texture blurring, further exacerbating noise. Furthermore, some errors that could be recovered at high-resolution layers with rich textures and clear edges cannot be recovered if the upper-layer matching range systematically shifts, as the lower layer will lose the possibility of searching for the true solution. Therefore, effective removal of upper-layer errors and the use of structured strategies to restore architectural details are crucial for overall reconstruction accuracy.

[0039] In urban area reconstruction, elevation anomalies typically manifest as local outliers, including fine-grained impulse noise and large-scale structural deviations caused by occlusion or pixel mixing. To address this characteristic, this invention proposes an error removal and structured filling strategy based on the coarse-layer dense matching results of the pyramid iterative framework: First, median filtering is used to suppress isolated impulse noise; second, eight-neighbor connected domain analysis is used to detect and remove large-scale outlier regions, while cross-validation of matching results using left-right consistency checks improves the robustness of the removal. After noise removal, void areas inevitably appear in the coarse-layer elevation. Given that urban building roofs are flat and have significant abrupt changes from the ground surface, direct mean filling easily leads to cross-structure mixing and edge blurring. This invention introduces prior guidance for filling based on building contours, using only pixels within the same building unit to complete elevation restoration, effectively maintaining the integrity of the roof plane and edge clarity, reducing error propagation between pyramid layers, and playing a crucial positive role in high-resolution matching.

[0040] 3.3. Full-size layer matching To further improve the elevation accuracy of pyramid iterative matching in urban areas, this invention adds a full-size layer below the bottom layer (i.e., the highest resolution object space) of the traditional pyramid. The elevation sampling interval of the initial object space at the bottom of the original pyramid is... x、 y、2 This resolution is sufficient to obtain structurally complete, geometrically correct, but still relatively coarse elevation results. Using this coarse-resolution matching result, the present invention constructs its corresponding initial elevation search range within the full-size layer; specifically, the horizontal intervals are consistent with the upper object space. x and y, elevation interval adopted z.

[0041] The elevation search range for the full-size layer is derived from the offset of the upper-layer elevation map, eliminating the need for morphological dilation and erosion operations. Simultaneously, the WTA algorithm is used to select the optimal elevation value as the final DSM output. For the elevation value Z, sub-pixel interpolation is performed using the neighborhood height and cost to obtain sub-pixel precision elevation values, increasing the continuity of the matching results. Finally, Equation 9 is used to determine the optimal elevation value. The selection strategy involves choosing a coordinate system, where c0 represents the cost corresponding to the current elevation value Z, and c1 and c2 represent the costs corresponding to elevation values ​​Z-1 and Z+1, respectively. Finally, a DSM image is generated by adding coordinates, datum information, and other details. This selection strategy avoids the enormous computational cost of full-space search, enabling full-size layers to stably converge to accurate reconstruction results at finer elevation intervals.

[0042] By adding this full-size layer to the image pyramid structure, the entire matching framework further improves the detail recovery capability of elevation estimation without significantly increasing the computational burden, and can effectively cope with the challenges of scenarios with frequent elevation changes and rich structural details in urban areas. To demonstrate the effectiveness of the proposed technology, this invention selected Lantau Island (green), Hong Kong Island (blue), and the Hong Kong Peninsula (red) within the Hong Kong Special Administrative Region as the experimental scope. From these, 25 representative regions of interest were carefully selected as experimental areas for DSM generation. The quality difference between the proposed method and the currently mainstream traditional open-source semi-global methods for generating DSMs was compared. The distribution of the experimental areas is shown below. Figure 5 As shown, these areas cover more than ten typical large-scale urban infrastructure scenarios, each with unique architectural structural features and reconstruction challenges, enabling a comprehensive evaluation of the algorithm's performance under different conditions. Furthermore, considering the presentation needs during qualitative analysis, five of the most representative areas were selected from the chosen 25 regions to demonstrate the DSM and mesh reconstruction effects, such as... Figure 5 As shown. The specific characteristics of the five most representative regions are as follows: 1) Area A: This area centers on the Hong Kong Velodrome (an oval-shaped stadium) and Tseung Kwan O Sports Ground (a crescent-shaped stadium), surrounded by a dense cluster of high-rise buildings. The challenge in this scenario lies in managing the occlusion relationship between the complex curved structure and the tall buildings.

[0043] 2) Area B: This area includes the Hung Hom Coliseum (a square stadium), Hung Hom Station (an integrated transportation hub composed of multiple rectangular units), and the International Metropolitan Shopping Centre (a cluster of high-rise buildings), surrounded by the Hong Kong Polytechnic University campus buildings and dense low- to mid-rise commercial buildings. The challenge of this scene lies in handling the mixture of various geometric shapes and the mutual occlusion of the high-density building cluster.

[0044] 3) Area C: This area centers around the Hong Kong Cultural Centre (inverted cone-shaped twin towers + wave-shaped roof), and includes semi-enclosed arc-shaped commercial buildings, the Canton Road boutique street commercial complex (diamond-shaped + cylindrical high-rise buildings), and surrounding dense low-rise buildings. The challenge of this scene lies in handling the reconstruction of the irregular curved surface structure.

[0045] 4) Area D: This area mainly consists of the rectangular plant buildings and cylindrical storage tanks of the Qingshan Power Plant. The difficulty in its reconstruction lies in dealing with the large-scale regular geometry and repetitive structure of industrial buildings.

[0046] 5) Area E: This area centers around the Hong Kong Convention and Exhibition Centre (a wing-shaped, arched roof exhibition hall), and includes the Grand Hyatt and Renaissance hotels (enclosed high-rise buildings) as well as a dense cluster of mixed-use commercial and residential buildings. The challenge in this scenario lies in managing the large-span spatial structure and the spatial hierarchy of the building complex.

[0047] To comprehensively evaluate the reconstruction performance of the proposed method in diverse urban built environments, this invention conducted experimental analyses at both qualitative and quantitative levels. Qualitatively, five representative regions were selected, and their reconstructed digital surface models (DSMs) and 3D mesh models were displayed to visually assess elevation representation capabilities and structural reconstruction quality. Quantitatively, error and completeness statistics were performed on the DSMs of the five regions (Table I), and further, error and completeness statistics were conducted on all 25 regions of interest (Table II) to comprehensively reflect the reconstruction accuracy and robustness of the algorithm in large-scale urban scenarios.

[0048] Table 1 compares the average accuracy and completeness of DSMs generated using different methods in five of the most representative sub-regions, presented by region partitioning. Area A: This area mainly consists of crescent-shaped and oval-shaped stadiums and high-rise buildings, with significant shadows near the high-rises and oval stadiums. All methods used in the comparison within this area... Figure 6 Both methods were able to largely reconstruct the terrain changes and the overall outline of the stadium building, but significant differences were observed when dealing with surrounding high-rise buildings. ASP's reconstruction in this area was unsatisfactory; its reconstructed stadium outline exhibited significant boundary expansion, a rough roof surface, and outlier protrusions. Furthermore, densely distributed high-rise buildings showed large areas of indistinguishable, fused forms, while buildings in sparse areas showed structural gaps and impaired spatial continuity. S2P's reconstruction was slightly better than ASP's, preserving the structure of some buildings, but the fused forms prevented it from accurately representing the true form of the buildings. Although S2P(M) employed a multi-view fusion strategy to improve reconstruction quality, the strategy lacked robustness, resulting in the adoption of large parallax values ​​in some areas. Figure 5In the previous method, the elevation of the hilltop area was significantly overestimated, deviating severely from the actual terrain. Furthermore, several high-rise buildings were not successfully reconstructed, demonstrating the instability of the fusion strategy in highly complex areas. In contrast, the method proposed in this invention shows significant advantages in the reconstruction of the stadium and surrounding high-rise buildings. The DSM results exhibited consistency between height continuity and the actual elevation, with clear boundaries and natural transitions in the roof structure. Particularly in the area of ​​the upper left corner of the stadium affected by self-shadowing, other methods showed uncertainty leading to edge defects, while the method of this invention maintained structural integrity, successfully reconstructing a complete elliptical outline and demonstrating good shadow adaptability. In terms of quantitative evaluation, as shown in Table I, the method of this invention achieved an average MAE of 4.174m and an RMSE of 5.920m in this area, significantly outperforming the results of ASP and S2P. Although S2P(M) achieved the highest value in terms of integrity, its average MAE and RMSE were 7.016m and 7.530m respectively, reflecting that its accuracy was severely affected by fusion errors.

[0049] Area B: This area is dominated by a square stadium and its surrounding high-density cluster of low- to mid-rise buildings. The buildings are densely packed, with a complex topological structure, exhibiting significant mutual shading, shadow shading, self-shading, and parallel shadowing, constituting a typical urban dense area reconstruction scenario. In the reconstruction of this area, the ASP (Average Spanning Floor) performed poorly overall, with significant elevation errors, allowing only a vague identification of some major building structures. For example... Figure 7 As shown, densely distributed low- and mid-rise building areas generally exhibit excessively high elevation values ​​and overlapping building forms, while high-rise building areas suffer from significantly underestimated reconstructed elevations, making it impossible to discern their actual heights and resulting in severe structural information loss. S2P shows some improvement in geometric representation compared to ASP, accurately reconstructing the structure of the central stadium and some surrounding low-rise buildings. The low- and mid-rise buildings in the lower left corner, with their outlines identifiable, still exhibit local overlap and elevation underestimation issues. The stadium's main body elevation is accurately reconstructed, with clear edge outlines, and no complete loss of high-rise buildings. However, in the high-density building area in the upper right corner, S2P also demonstrates poor geometric differentiation capabilities, exhibiting significant elevation anomalies and structural blurring similar to ASP, resulting in a substantial decrease in reconstruction quality. S2P(M) improves some edge representation by introducing multi-view information, particularly in the boundary outline representation of the area surrounding the stadium. However, from... Figure 7It can be observed that, limited by the insufficient robustness of the fusion strategy, this method exhibits elevation overestimation in several areas. In the lower left corner, where there are low-rise buildings, multi-view fusion does indeed provide more accurate elevation restoration compared to binocular S2P. However, in flat areas, such as road areas, the elevation of buildings and the ground is not effectively separated, resulting in blurred and blended surface details. The problem persists in the densely populated upper right corner, showing no significant improvement. In contrast, the method proposed in this invention demonstrates superior overall performance in this area. The elevation representation in the stadium area is similar to the S2P result, but the boundaries are... Figure 7 The resulting structures exhibit sharper and more regular geometric shapes, demonstrating a stronger ability to preserve details. The adhesion problem of low-rise buildings in the lower left corner is significantly alleviated, and the elevation estimation is closer to the actual ground conditions, although there is still a certain degree of underestimation. In the dense area of ​​the upper right corner, some buildings show adhesion and error accumulation, but the elevation deviation is smaller than that of S2P(M), and the building shapes are clearer. In terms of quantitative assessment, as shown in Table I, the method of this invention obtained an average MAE of 3.931m and an RMSE of 4.993m in this area, with accuracy superior to ASP and S2P. Although S2P(M) achieved a reconstruction integrity rate of 89.04%, due to the lack of a robust fusion strategy for high-density urban building scenarios, its average MAE and RMSE increased by 0.495m and 1.056m respectively compared to S2P.

[0050] Area C: This area contains a large number of high-rise commercial buildings and low- to mid-rise buildings. These buildings exhibit complex, irregularly shaped curved structures and are closely arranged, yet the relationships between them are relatively distinct. The problems are similar to those in Area B, but the average building height in this area is higher than in Area B, resulting in more severe shading and obstruction caused by the building height. In the reconstruction results for this area, the ASP (Average Scaling Average) remained unsatisfactory. Figure 8 As shown, high-rise buildings generally exhibit significantly overestimated elevation values, while mid- and low-rise buildings are severely integrated with roads, and in some areas, structural information is completely missing, making it difficult to identify the true outlines of features from the results. S2P shows a significant improvement in reconstruction quality compared to ASP. It can not only accurately restore the overall form of most high-rise commercial buildings but also capture the elevation distribution of mid- and low-rise buildings to some extent. However, due to the unresolved issue of building boundary integration in densely populated areas, [further issues remain]. Figure 8The grid structure presented in the image shows incomplete building structures, making it difficult to accurately reflect the independence between buildings. While S2P(M) enhances the representation of some edge contours by introducing a multi-view fusion strategy, it still suffers from problems similar to those in regions A and B: insufficient robustness of the fusion process, leading to a general overestimation of building elevations in high-density areas. Incorrect parallax information is incorporated into the fusion process, causing elevation values ​​in some areas to deviate from actual ground heights, affecting the accuracy of the overall structural representation. In contrast, the method proposed in this invention performs better in this region. In the area on the left, dominated by high-rise commercial buildings, our reconstruction results effectively distinguish individual buildings, with clear structural edges and accurate building mass representation, reproducing building outlines and facade changes more meticulously than S2P. On the right side of the region, although some low- and mid-rise buildings experience localized information loss due to occlusion and reflection issues, and some building structures are excessively blended into the ground, the overall elevation estimation is still superior to other methods, especially demonstrating stronger suppression of adhesion and structural misjudgment. As shown in Table I, the quantitative evaluation results indicate that the method of this invention achieved the best error performance among all methods in this region, with an average MAE of 3.397m and an RMSE of 4.984m. S2P achieved similar reconstruction accuracy, with an average MAE of 3.498m, but it was slightly inferior in terms of boundary representation and structural connectivity in dense regions. Although S2P(M) achieved a 5.28% higher integrity rate than the method of this invention, its average MAE and RMSE reached 3.930m and 5.466m respectively, which were 0.533m and 0.482m higher than the method of this invention. This further indicates that its accuracy was significantly affected during the fusion process, making it unable to balance integrity and geometric accuracy.

[0051] Region D: This region primarily consists of industrial buildings of the factory type, characterized by large-scale, regularly distributed rectangular roofs. These roofs often lack sufficient texture information or exhibit large areas of repetitive textures, posing a significant challenge to stereo matching-based 3D reconstruction methods, particularly in accurately estimating elevations and distinguishing building boundaries. In the reconstruction results for this region, ASP showed improved performance compared to its performance in areas with dense high-rise buildings, and was able to largely restore the main building structures. However, auxiliary buildings, due to their strong texture continuity with the main factory buildings, were often incorrectly merged and reconstructed, resulting in blurred, blocky structures. Figure 9As can be observed, the ground area reconstructed by ASP suffers from weak texture and a lack of effective parallax information, leading to the smoothing of details such as pipes and roads, resulting in a systematic overestimation of ground elevation and difficulty in distinguishing road structures. S2P performs similarly to ASP in this area, reconstructing the main form of the factory building while also exhibiting problems of attached auxiliary buildings and missing ground information. This indicates that the method's performance in handling large-scale, low-texture areas is still limited by the insufficient robustness of traditional matching strategies. After fusing multi-view information, S2P(M) corrects the elevation anomalies present in S2P in some building roof areas, making the elevation estimation closer to the actual terrain. However, the adhesion of auxiliary buildings and the loss of ground details are still not fundamentally resolved, indicating that its multi-view fusion strategy still suffers from insufficient reliability in scenarios with repetitive textures and boundary continuity. In contrast, the method proposed in this invention demonstrates stronger building separation and elevation recovery capabilities. Figure 9 As can be observed, our results accurately distinguish between the main factory building and auxiliary buildings, with clear boundaries and distinct structures, significantly improving the problem of easy adhesion in low-texture scenes using traditional methods. Simultaneously, in texture-scarce areas such as rooftops, our method maintains the continuity and stability of elevation estimation, reconstructing the roof planar features relatively completely and avoiding the generation of abnormal elevation points. From a quantitative evaluation perspective, as shown in Table I, our method achieves the best error performance in terms of average MAE and RMSE in this region. Particularly in terms of RMSE, it is 0.77m lower than the second-ranked method, demonstrating a significant accuracy advantage. Regarding completeness, the difference between our method and S2P(M) is small, only 0.56% lower.

[0052] Region E: This region contains an exhibition hall with a complex, wing-shaped, arched roof structure and its surrounding densely packed high-rise buildings. The scene exhibits occlusions similar to Region A, and due to its complex roof structure, it also suffers from penumbra shadows. It serves as an ideal sample for evaluating the performance of 3D reconstruction algorithms in handling scenes with complex curved buildings and vertical elevation differences. In the reconstruction of this region, ASP performed poorly, especially in handling the complex roof structure of the exhibition hall. The wing-shaped roof, which should have had a distinct arc and undulation, was... Figure 10The structure was simplified to an arc-shaped ridge on the ground, lacking three-dimensionality and structural information, making it virtually impossible to identify its architectural attributes. Furthermore, the reconstruction of surrounding high-rise buildings failed to effectively restore the original structure, exhibiting severe deficiencies in elevation information and structural damage, consistent with the performance of ASP in the previous areas. S2P achieved better results in the reconstruction of the exhibition hall roof, largely restoring its arc-shaped structural features. The elevation estimation was natural and continuous, closely resembling the curvature changes of a real roof. However, in the high-rise building area, due to the dense building distribution and weak facade texture features, the S2P results showed significant adhesion. Although some elevation information was recovered, the lack of clear boundary segmentation resulted in a generally poor overall structural representation, failing to accurately reflect the actual building outline. S2P(M), utilizing a multi-view fusion strategy, further improved the reconstruction performance in the exhibition hall roof area. The elevation changes of its arc-shaped roof were smoother and more consistent with the real curved surface structure, with a more continuous and realistic geometric expression. However, the elevation offset problem during multi-view fusion was not completely resolved, resulting in a systematic overestimation of the roof elevation. The method proposed in this invention demonstrates comprehensive advantages in this area. Compared to S2P, our method can also accurately reconstruct the elevation undulations of the pavilion's arched roof, but it offers clearer and sharper edge contours and better structural continuity. In the reconstruction of surrounding high-rise buildings, our method continues its advantages in the aforementioned area, effectively separating dense buildings, reducing adhesion, achieving three-dimensional separation between buildings, and providing more accurate elevation estimation, complete morphology, and good spatial consistency. Quantitative results, as shown in Table I, demonstrate that our method achieves the best reconstruction accuracy in this area, with an average MAE of 2.218m and an RMSE of 5.794m, representing improvements of 3.725m and 0.260m respectively compared to the worst-performing ASP, demonstrating significant advantages. Regarding reconstruction completeness, our method achieves 76.215%, only about 2.13% lower than the multi-view fusion method S2P(M), maintaining high accuracy while ensuring good reconstruction coverage.

[0053] To further verify the universality and robustness of the method of this invention in large-scale complex urban environments, we conducted a systematic quantitative evaluation on all 25 regions of interest, and statistically analyzed the reconstruction accuracy and completeness of the digital surface model (DSM) in different regions. The results are shown in Table II. Furthermore, we compared the DSM reconstruction results for different aberrations using different methods, such as... Figure 11As shown, all areas are divided into three main sections based on geographical location: Area I (green) is the outlying islands, Area II (blue) is Hong Kong Island, and Area III (red) is the Kowloon Peninsula. Area I mainly consists of large-scale, low-rise buildings such as airport terminals and amusement parks, with minimal topographic relief and strong roof structure continuity. In this area, the average reconstruction error (MAE) and RMSE of our method are 3.437 meters and 4.301 meters, respectively, which are only 0.317 meters and 0.019 meters different from the suboptimal method S2P, demonstrating considerable geometric accuracy. In terms of completeness, although slightly lower than S2P's 86.36%, it still reaches 84.09%, showing that our method has good adaptability and efficient elevation recovery capabilities in large-scale building reconstruction scenarios. Area II contains a large number of densely built-up mid-to-high-rise buildings and large-scale infrastructure along the coast, and is a typical high-density urban core area with complex structure and severe obstruction. In this region, our method demonstrates a significant advantage in accuracy. Compared to the second-best algorithm, the average MAE and RMSE are improved by 0.652 meters and 1.093 meters, respectively, fully demonstrating that our proposed method has superior matching and elevation recovery capabilities when dealing with high-rise, densely built-up buildings. ASP has the lowest reconstruction accuracy in this region, consistent with its analysis in Table I. Although the multi-view fusion strategy S2P(M) achieved the highest integrity rate in this region (5.45% higher than our method), its average MAE and RMSE are 1.574 meters and 1.394 meters higher than our method, respectively, indicating that its fusion mechanism suffers from geometric distortion in complex building scenes. Region III, located in the Kowloon Peninsula, features the most diverse building types, including mid-to-high-rise residential areas near the coastline, dense low-rise building clusters, large-scale transportation hubs, and industrial facilities. Some areas are built on hillsides, resulting in significant interference from natural topography. In such highly heterogeneous mixed scenarios, our method still exhibits good stability and accuracy, with average MAE and RMSE of 4.311 meters and 5.505 meters, respectively. In terms of reconstruction completeness, it differs from S2P(M) by only 0.95%, further validating the robustness and generalization ability of our method in environments with multi-source occlusion, complex textures, and drastic geometric changes. Overall, the evaluation results shown in Table II fully demonstrate that the method proposed in this invention has excellent geometric accuracy and high completeness in handling different types of urban structures, especially in high-rise dense areas and large-scale facility scenes, where it has a greater reconstruction advantage than traditional binocular and multi-view fusion methods.

[0054] Table 2 compares the average accuracy and completeness of DSM generation using different methods across twenty-five regions (based on three major regions), presented by region. It is important to note that the schemes and arrangements of this application shown in the exemplary embodiments are merely exemplary. Although only a few embodiments are described in detail in this disclosure, those who consult this disclosure will readily understand that many modifications are possible (e.g., variations in various parameter values ​​(temperature, power, humidity, etc.), installation arrangements, names, colors, logical orders, etc.) without substantially departing from the novel teachings and advantages of the subject matter described in this application. Therefore, all such modifications are also included within the scope of the invention, and the order or sequence of any process or method steps may be changed or rearranged according to alternative embodiments. In the claims, any "apparatus plus function" clause is intended to cover the structure described in this invention for performing the function, and not only structural equivalents but also equivalent structures. Other substitutions, modifications, changes, and omissions may be made in the design, operation, and arrangement of the exemplary embodiments without departing from the scope of the invention. Therefore, the invention is not limited to the particular embodiments but extends to various modifications that still fall within the scope of the appended claims.

[0055] Furthermore, in order to provide a concise description of exemplary embodiments, not all features of actual embodiments (i.e., those features that are not relevant to the best mode of carrying out the invention as currently considered, or those features that are not relevant to implementing the invention) may be omitted.

[0056] It should be understood that numerous specific implementation decisions can be made during the development of any practical implementation, such as in any engineering or design project. Such development efforts may be complex and time-consuming, but for those skilled in the art who benefit from this disclosure, the development effort will be a routine work of design, manufacturing, and production without requiring much experimentation.

[0057] It should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit it. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can be made to the technical solutions of the present invention without departing from the spirit and scope of the technical solutions of the present invention, and all such modifications or substitutions should be covered within the scope of the claims of the present invention.

Claims

1. A method for fine satellite 3D reconstruction of large-scale urban buildings based on prior-guided iteration, characterized in that... Includes the following steps: Step 1: Obtain three types of data for the area to be constructed: panchromatic satellite optical imagery, DEM, and image sequence of building vector data with building height information. Select one image as the reference image and the others as source images. Then, use the angle filtering method to obtain candidate image sequence for subsequent calculations. Step 2: Resample the grayscale information of the reference image to the object space to obtain a grayscale base map. Then, register the building vector data with floor height information with the grayscale base map, and then fuse it with the DEM resampled to the object space to generate an elevation prior. Step 3: By downsampling the candidate image sequence step by step, an image pyramid is constructed. Then, an improved object-side semi-global dense matching method is used to iteratively match the image pyramid layer by layer from top to bottom to obtain the reconstructed DSM. The improved object-side semi-global dense matching method uses CENSUS to calculate the cost value, and then sequentially uses the building prior-guided cross-arm cost aggregation method and the grayscale base map-guided scan line cost aggregation method to perform two-stage aggregation optimization of the cost value. At the same time, when performing layer-by-layer iterative matching on the image pyramid, the elevation prior is used to initialize the elevation search range of the top layer, and the elevation search range of the remaining layers is constrained by the matching results of the upper layer to constrain the elevation search range of the lower layer.

2. The method for fine satellite 3D reconstruction of large urban buildings based on prior-guided iteration as described in claim 1, characterized in that: When constructing the image pyramid, the source image, reference image and corresponding object space to be matched are downsampled step by step. Except for the full-size layer at the bottom, the sampling interval of each layer is performed according to the following relationship. image: , ; Object space: , , ; in, and V represents the horizontal and vertical displacement step size of the image, respectively. , Y, Z represents the displacement step size in the latitude, longitude, and height directions in the object space, respectively. Indicates the layer number; For the full-size layer (layer 0), the image is not downsampled, and the displacement step size relationship in the object space is set to... , , .

3. The method for fine satellite 3D reconstruction of large-scale urban buildings based on prior-guided iteration as described in claim 2, characterized in that: During the layer-by-layer iterative matching, an improved object-side semi-global dense matching method is used for the current layer of the image pyramid. Perform matching to obtain the elevation value corresponding to the horizontal coordinate v of each object point. Together they form an elevation map Then check the elevation map Morphological erosion and morphological dilation were performed to obtain elevation values. and Then the next layer of the image pyramid The search range corresponding to the horizontal coordinate v is expressed as ),in, This represents the elevation displacement step size. The step size method coefficients are... For full-size layers, the search range is represented as follows: ),in, The elevation value is the one that matches the previous layer adjacent to the full-size layer.

4. The method for fine satellite 3D reconstruction of large-scale urban buildings based on prior-guided iteration as described in claim 3, characterized in that: The search range at the top of the image pyramid is based on an elevation prior definition, which includes ground elevation. Building outline buffer zone elevation and building elevation within the building outline Three types, Ground elevation minimum elevation and highest elevation The definition is as follows: in, Elevation values ​​from reference DEM imagery; Building outline buffer zone elevation minimum elevation and highest elevation The definition is as follows: in, Based on reference building height; Building elevation within the building outline buffer zone minimum elevation and highest elevation The definition is as follows:

5. The method for fine satellite 3D reconstruction of large-scale urban buildings based on prior-guided iteration as described in claim 1, characterized in that: For the improved semi-global dense matching method for object space, the following formula is used to calculate the cost. , in, This represents the transformation value corresponding to the grayscale value within the projection window on the reference image. This represents the transformation value corresponding to the grayscale value within the projection window on the i-th source image, where N is the total number of source images. This represents the cost calculated for the i-th source image and the reference image. The weighted average of all cost values ​​is the final cost value, where u and v represent the vertical and horizontal coordinates in the image, and X, Y, and Z represent the coordinates of any point in the object space in the longitude, latitude, and elevation directions, respectively. Using cross-arm cost aggregation, given a point V in the object space. When it finds the endpoint When the arm stops in a certain direction, the stopping rule at that point is as follows: in, This represents the color difference between two object points; the corresponding grayscale information comes from the grayscale base image. and Indicates the grayscale threshold. This represents stricter grayscale constraints. Represents the spatial distance between two object points. and Indicates the distance threshold. This represents the previous object point, taking the left arm as an example. Indicates the roof area. Indicates the ground area; The calculation formula for the scan line cost aggregation is as follows: in, Let V represent the coordinates of a point V in object space. Equivalent to the horizontal coordinates (X, Y) of the object point, where Z represents the elevation coordinates of the object point. This represents the aggregated cost of the scan line at the assumed elevation Z for point V. Indicates the object point after the cross arm cross cost aggregation. Given the cost at elevation Z, r defines the convergence direction. and This represents the penalty term, which satisfies... , express The initial values ​​correspond to small and large parallax changes, respectively. and Representing the object point in object space respectively and The horizontal coordinates of the grayscale base image The corresponding grayscale value in the image.

6. The method for fine satellite 3D reconstruction of large-scale urban buildings based on prior-guided iteration as described in claim 1, characterized in that: In step one, when the candidate image sequence is obtained by using the angle screening method, the image with the smallest solar zenith angle and the imaging angle closest to the vertical viewing angle is selected from the image sequence as the reference image, and the remaining images are the source images. Then, the cross angle and solar angle difference between the reference image and each source image are calculated one by one. Finally, the source images with cross angle and solar angle differences greater than a set threshold and the reference images are selected to form the candidate image sequence.

7. The method for fine satellite 3D reconstruction of large-scale urban buildings based on prior-guided iteration as described in claim 6, characterized in that: The intersection angle between the reference image and any source image is calculated using the following formula. Difference with solar angle , in, and These represent the angles of the principal point vectors of the reference image and the source image, respectively. and These represent the solar zenith angle and solar azimuth angle of the reference image, respectively. and These represent the solar zenith angle and solar azimuth angle of the source image, respectively.

8. The method for fine satellite 3D reconstruction of large urban buildings based on prior-guided iteration as described in claim 1, characterized in that: In step two, the building vector data with floor height information is first divided into multiple small blocks. Then, local affine or homography transformation is used to register each small block with the corresponding area of ​​the grayscale base map. The registration result is then rasterized and resampled together with the DEM into the object space for overlay and fusion to generate the elevation prior.

Images