Satellite photogrammetry point cloud optimization method for large-scale building area

Abstract

The invention discloses a satellite photogrammetry point cloud optimization method for a large-scale building area. Belongs to the technical field of remote sensing surveying and mapping, and particularly relates to the technical field of satellite photogrammetry point cloud optimization methods for large-scale building areas. And satellite photogrammetry point cloud optimization of a large-scale building area is realized. The method comprises the following steps: S1, based on a stereo image pair, generating a photogrammetric point cloud set # imgabs0 #, segmenting # imgabs1 # into a ground point set # imgabs2 # and a non-ground point set # imgabs3 # S2, and carrying out European clustering on # imgabs4 # to realize building monotonization: S3, carrying out facade removal based on an independent building object; s4, roof redundant points are constrained based on the building contour vector; s5, repairing ground holes based on grid projection to obtain final point cloud set # imgabs5 #

Description

Technical Field

[0001] The present invention relates to the technical field of remote sensing mapping, and in particular to the technical field of satellite photogrammetry point cloud optimization methods for large-scale building areas. Background Art

[0002] A realistic 3D model is a high-precision, highly realistic 3D digital model built based on real-world scene data. By integrating multi-source remote sensing data (aerial imagery, satellite imagery, and LiDAR point clouds) with ground-based measurement data, and leveraging photogrammetry, computer vision, and 3D reconstruction technologies, it creates a comprehensive, multi-level, and detailed representation of real-world geographic entities (topography, buildings, vegetation, and roads).

[0003] Compared with airborne point clouds and vehicle-mounted point clouds collected by LiDAR (Light Detection And Ranging), photogrammetry point clouds (PhoDAR) generated based on satellite images have become the main data source for city-level real-scene 3D modeling due to their wide coverage and relatively low acquisition cost. However, PhoDAR cannot penetrate obstacles like LiDAR and can only obtain information on the ground surface, vegetation canopy, and the top and partial side of buildings. Moreover, it is difficult to obtain complete side information of buildings due to the viewing angle of satellite or aerial images. In the forward intersection process of photogrammetry, due to the obstruction of the roof, multiple beams of light from one image intersect with the same beam of light from another image, generating multiple erroneous intersection points in one direction, and ultimately forming an erroneous "slope" effect on the photogrammetry point cloud, such as Figure 2 shown.

[0004] Therefore, in the photogrammetric point cloud generated based on satellite imagery, the sides of buildings in densely populated urban areas often appear tilted rather than vertical walls. This will cause distortion in the DSM (digital surface model) and the real-life 3D model, making it impossible to truly reflect the physical spatial relationship and undermining the usability of the real-life 3D model in urban planning and engineering design. Secondly, when texture mapping the real-life 3D model, incorrect intersection points in the photogrammetric point cloud will lead to confusion in the mapping coordinates, resulting in distorted texture stretching, dislocation, and fractures on the building surface, destroying the visual effect. Due to the large shooting area of satellite images, the problem of side occlusion exists in almost all buildings of a certain height, resulting in a huge workload for manually removing erroneous points, requiring the use of automated algorithms. Summary of the Invention

[0005] In response to the problems of data acquisition limitations, three-dimensional model geometric distortion, texture mapping quality impairment and low manual processing efficiency in the existing remote sensing mapping field, the present invention provides a satellite photogrammetry point cloud optimization method for large-scale building areas. By performing layered filtering and individualization on the point cloud, removing facades and constraining roofs of building targets, and repairing ground holes through grid discretization, the present invention achieves satellite photogrammetry point cloud optimization for large-scale building areas.

[0006] The method comprises the following steps:

[0007] S1. Generate photogrammetric point cloud based on stereo image pairs and will Segmentation into ground point sets and non-ground point sets

[0008] S2, yes Euclidean clustering is used to achieve building individualization;

[0009] S3. Remove the facade based on the independent building object to obtain the roof point cloud Will After Euclidean clustering, Z roof point cloud sets are obtained

[0010] S4, constraining rooftop redundant points based on building outline vector;

[0011] S5. Repair ground holes based on grid projection to obtain the final point cloud set

[0012] Further, the Segmentation into ground points and non-ground points This is done through the cloth filtering algorithm.

[0013] Further, the The specific method of using Euclidean clustering to realize building monomerization is as follows:

[0014] S31, will The initial state of is marked as "unvisited";

[0015] S32, will All non-ground points in are clustered into different clusters C k , and mark the status of the clustered non-ground point set as “visited”;

[0016] S33, iterative update C k , obtain M independent building objects Ω C ={C0,C1,…,C M}, k = [0, 1,…, M].

[0017] Furthermore, the iterative update C k Specifically: Among them, t represents the number of iterative updates, p j represents the jth non-ground point cloud, Indicates p j The point set within the neighborhood radius ∈; when When Ω is obtained, the iteration ends. C .

[0018] Further, Ω C Constraints need to be met: Ω C ={C k |N min ≤|C k |≤N max}, where |C k | represents cluster C k The number of points, N min Indicates C k The minimum number of points, N max Indicates C k The maximum number of points.

[0019] Furthermore, the facade removal based on independent building objects is performed by the formula: Calculate, where [m s ,m e ] represents the interval index of the roof slice, Represents the final roof interval, m represents the mth C k Slice, p i represents the i-th non-ground point cloud, z i Indicates p i Elevation slice in the z-axis direction, z min Indicates p i The minimum distribution range in the z-axis direction, Δz represents p i The difference between the maximum and minimum distribution ranges in the z-axis direction.

[0020] Furthermore, the roof redundancy points are specifically constrained based on the building outline vector: For constraints, optimize Get N optimized roof point cloud sets Among them, ρ represents the overlap rate between the building outline vector and the center of the roof point cloud, α represents the minimum overlap rate threshold, Represents the set of building outline vectors closest to the center of the roof point cloud, v l *Represents the building outline vector that has the highest overlap with the center of the roof point cloud.

[0021] Furthermore, the grid projection-based ground hole repair is specifically as follows:

[0022] S81, will and Projected onto the xOy plane and divided into a regular grid;

[0023] S82, traverse the regular grid to obtain the coordinates of the hole point (x h ,y h ,z h ), (x h ,y h ,z h ) Get a complete ground point cloud set

[0024] S83, will and Take the union and get The beneficial effects of the present invention are:

[0025] The method described in this paper first uses a cloth filtering algorithm to filter large-scale urban point clouds into ground point clouds and non-ground point clouds, centrally processing erroneous intersection points generated by non-ground point clouds. Secondly, Euclidean clustering is used to separate ground features, and density and height difference detection are used to distinguish between buildings and low-lying features (including bushes and trees). Building facades are then removed by slicing along elevations based on differentials, and roof shapes are constrained based on building vectors. Finally, based on the restored non-ground portions, areas missing ground points are filled in to obtain a complete optimized point cloud for the target scene.

[0026] The method described in the present invention can correct the urban building point cloud, thereby improving the precision and aesthetics of the produced DSM and real-scene three-dimensional model. It is easy to implement and deploy and can be applied in actual engineering. BRIEF DESCRIPTION OF THE DRAWINGS

[0027] Figure 1 This is a flow chart of the method according to an embodiment of the present invention;

[0028] Figure 2 The building point cloud according to the embodiment of the present invention generates an erroneous front intersection effect diagram;

[0029] Figure 3 Schematic diagram of the CSF ground filtering algorithm according to an embodiment of the present invention;

[0030] Figure 4 This is a schematic diagram of facade separation based on density slicing according to an embodiment of the present invention;

[0031] Figure 5 This is a rendering of the separation of the roof and facade according to an embodiment of the present invention;

[0032] Figure 6 1. A comparison diagram of a rooftop redundant point (DSM) and a remote sensing image according to an embodiment of the present invention, wherein (a) is a schematic diagram of a rooftop redundant point (DSM), and (b) is a schematic diagram of a remote sensing image.

[0033] Figure 7 Comparison diagram of roof constraint effects according to an embodiment of the present invention, (a) shows a vector diagram of a building, (b) shows a DSM diagram of the point cloud before constraint, and (c) shows a DSM diagram of the point cloud after constraint.

[0034] Figure 8 This is a rendering of the ground repair effect according to an embodiment of the present invention;

[0035] Figure 9 This is a schematic diagram of the white model effect of the three-dimensional model before the local point cloud optimization of the experimental area 1 according to the embodiment of the present invention;

[0036] Figure 10 This is a schematic diagram of the white model effect of the three-dimensional model after the optimization of the local point cloud of the experimental area 1 according to the embodiment of the present invention;

[0037] Figure 11 This is a schematic diagram of the white model effect of the three-dimensional model before the local point cloud optimization of the experimental area 2 according to the embodiment of the present invention;

[0038] Figure 12 This is a schematic diagram of the white model effect of the three-dimensional model after local point cloud optimization in the experimental area 2 according to an embodiment of the present invention. DETAILED DESCRIPTION

[0039] The technical solution of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the embodiments described are only some embodiments of the present invention, not all embodiments. All other embodiments obtained by ordinary technicians in this field based on the embodiments of the present invention without making any creative efforts shall fall within the scope of protection of the present invention.

[0040] This embodiment provides a satellite photogrammetry point cloud optimization method for large-scale building areas. This embodiment uses Jilin-1 Gaofen 04A satellite image data to complete regional block adjustment, epipolar image generation, dense matching and three-dimensional point cloud solution to obtain photogrammetry point clouds, and uses L1C fused images to obtain building vectors. Microsoft Studio 2022 is used as a simulation platform, and the PCL library is used to verify the method described in this embodiment. L1C fused images refer to satellite image products generated by subjecting satellite raw data (L0) to radiometric correction, geometric correction and multispectral / panchromatic band fusion processing.

[0041] The flow chart of the method is as follows Figure 1 As shown, the method includes the following steps:

[0042] S1. Generate photogrammetric point cloud based on stereo image pairs and will Segmentation into ground point sets and non-ground point sets

[0043] The relevant operations of step S1 are introduced with a specific example:

[0044] Generate a photogrammetric point cloud of the area using satellite image pairs taken at the same location that meet the base-to-height ratio requirements and have good image conditions. like Figure 2 As shown in Figure 2, the incorrect intersection points appear mainly on the building facade. The cloth simulation filter (CSF) algorithm is first used to filter the The basic idea of the cloth filtering algorithm is to invert the point cloud and create a uniform grid above the highest point of the point cloud. Assuming a piece of cloth is falling from above due to gravity, the cloth node stops moving when it collides with the point cloud, and the contact point is recorded.

[0045] like Figure 3 As shown, the relative positions between the nodes are constrained by the rigidity of the cloth to prevent excessive deformation. The gravity, collision, and internal force calculations are repeated until the cloth is stable or the maximum number of iterations is reached. After the iterative simulation of "cloth", the point cloud can be classified: if the distance from the point cloud to the corresponding cloth node is less than the threshold, it is marked as a ground point. Otherwise it is a non-ground point

[0046] S2, yes Euclidean clustering is used to achieve building individualization.

[0047] The relevant operations of step S2 are introduced with a specific example:

[0048] To detect false facades on every building, the first step is to singulate each building. This singulation relies on searching for connected components in a graph based on Euclidean distance. Its core concepts include neighborhood definition, iterative region growing, and cluster validity determination.

[0049] First, all non-ground points The initial state is marked as "unvisited", that is, set the array V:

[0050]

[0051] Among them, p irepresents the i-th non-ground point cloud, and Fasle indicates that the initial state is marked as “unvisited”.

[0052] In order to realize neighborhood search, for any point The point set within its neighborhood radius ∈:

[0053]

[0054] where ||·|| 2 is the Euclidean norm. In code deployment, the neighborhood search can be accelerated by KD-Tree. j represents the j-th non-ground point cloud.

[0055] Through the region growing process, All points in the cluster are clustered into different clusters C k And mark the clustered points as "True", that is, the access status:

[0056]

[0057] Next, cluster generation is achieved through iterative updates:

[0058]

[0059] Among them, ∪ means taking the union of the points of the t-1 cluster in the previous round and the new cluster in this round, and t means the number of iterations.

[0060] when When the cluster stops growing, we get Ω C .

[0061] Final cluster set Ω C The following point constraints must be met:

[0062] Ω C ={C k |N min ≤|C k |≤N max}

[0063] where |C k | represents cluster C k The number of points, N min Indicates C k The minimum number of points, N max Indicates C k The maximum number of points.

[0064] After clustering, we get M independent building objects Ω C ={C0,C1,…,C M}.

[0065] S3. Remove the facade based on the independent building object to obtain the roof point cloud Will After Euclidean clustering, Z roof point cloud sets are obtained

[0066] The relevant operations of step S3 are introduced with a specific example:

[0067] In this embodiment, the separation of the roof and the facade is achieved based on the following assumption: the roof points are concentrated at higher positions with a density significantly higher than that of the ground.

[0068] Input Ω in sequence C Each C k , that is, a single building point cloud, by analyzing the density distribution of a single building point cloud in the vertical direction, the building point cloud C is separated k Get roof point cloud and facade point cloud

[0069] Since the data format of the point cloud is an unordered three-dimensional point, this embodiment is based on the basic idea of differentiation and transforms the point cloud C into k Slicing is done to achieve regularization and discretization. First, calculate the distribution range of the point cloud on the z-axis:

[0070]

[0071] Among them, z i Indicates p i Elevation slices in the z-axis direction, represents z i Point cloud of a single building within the range, z min Indicates p i The minimum distribution range in the z-axis direction, z max Indicates p i The maximum distribution range in the z-axis direction.

[0072] Set the size of the point cloud slice along the z axis to Δz, and the total number of slices N can be calculated bin :

[0073]

[0074] Calculate the histogram H and convert z i In [z min +mΔz,z min +(m+1)Δz] i Divide into H(m):

[0075]

[0076] Where I(·) represents the characteristic function, which determines the elevation slice z i Whether it belongs to the mth slice, m means the mth C k slice.

[0077] Then, we distinguish roof and facade slices based on the following assumption: the density is higher than the minimum density threshold and it is considered as roof. The minimum density threshold is set as:

[0078]

[0079] Among them, λ represents the stretching coefficient, which controls the size of the threshold τ;

[0080] Traverse the histogram H and mark the set of slices that continuously satisfy H(m)≥τ:

[0081]

[0082] Among them, [m s ,m e ] represents the interval index of the detected roof slice, m s is the starting index, m e is the end index, each interval index corresponds to a height range along the z axis [z min +m s Δz,z nin +m e Δz].

[0083] In addition, due to the presence of subtle special structures on the top of the building, the density of these structures is not enough to constitute the roof, but they are part of the roof. Therefore, a peak expansion strategy is needed to find the global highest density interval:

[0084]

[0085] And expand to the top slice:

[0086]

[0087] Combine the general roof slice and the top roof slice to get the final roof interval:

[0088]

[0089] According to the interval index in point cloud C k Filter out the roof points and get the roof point cloud

[0090]

[0091] The facade point cloud is obtained by complement operation:

[0092]

[0093] However, some complex buildings may have complex structures, including multiple roofs and facades at varying elevations, making it difficult for subsequent design algorithms to separate these roofs and facades. Therefore, these buildings need to be re-single-unitized to ensure that a series of individual roof objects are ultimately obtained, facilitating the implementation of subsequent algorithms.

[0094] The specific implementation method is to convert the building point cloud C k The roof point cloud stripped out After another Euclidean clustering, the number of clusters obtained after clustering is determined. If the number of clusters exceeding the point threshold σ is greater than 1, it means that the building point cloud C is a composite structure containing multiple roofs and facades. After the second clustering, the roofs can be separated into independent objects, and finally Z roofs are obtained. The schematic diagram of vertical slicing along the z axis is as follows Figure 4 As shown, the separation effect of the roof and facade is as follows Figure 5 shown.

[0095] S4. Constrain the redundant roof points based on the building outline vector.

[0096] The relevant operations of step S4 are introduced with a specific example:

[0097] After removing the facade, only the roof remains in the building point cloud. However, the incorrect intersection points are not only presented as facade slopes, but also exist at the junction of the roof and the facade. The size of the building in such a point cloud will be much larger than the actual situation. Figure 6 As shown in (a) and (b), the DSM generated by the point cloud is significantly wider than the original remote sensing image in the north-south direction of the building.

[0098] Therefore, the shape of the roof must be strictly constrained, keeping its dimensions as close to the original image as possible. Many existing feature extraction algorithms can effectively extract building outlines from remote sensing images. This method leverages these algorithms to extract building outline vectors from the L1C fused image as ground truth, eliminating erroneous intersections of roof edges. Since vectors are two-dimensional data and do not contain elevation, this step first projects the three-dimensional points onto two dimensions using the RPC (Rational Polynomial Coefficients) model. This means that pixel coordinates are obtained given the longitude, latitude, and elevation.

[0099] Convert the three-dimensional coordinates of the point cloud under the UTM system into the ground point space coordinates (longitude Latitude λ and altitude h) and then normalized:

[0100]

[0101] in, represents the standard translation parameter, Represents the standard scale parameters. (U, V, W) represents the normalized ground point spatial coordinates. (c, r) represents the normalized image point coordinates, which are obtained by the following formula:

[0102]

[0103] Among them, each term is obtained as follows:

[0104] Numc(u,v,w)=a1+a2V+a3U+a4W

[0105] +a5VU+a6VW+a7UW+a8V 2

[0106] +a9U 2 +a 10 W 2 +a 11 UVW+a 12 V 3

[0107] +a 13 VU 2 +a 14 VW 2 +a 15 V 2 U+a 16 U 3

[0108] +a 17 UW 2 +a 18 V 2 W+a 19 U 2 W+a 20 W 3

[0109] Denc(u,v,w)=b1+b2V+b3U+b4W

[0110] +b5VU+b6VW+b7UW+b8V 2

[0111] +b9U 2 +b 10 W 2 +b 11 UVW+b 12 V 3

[0112] +b 13 VU 2 +b 14 VW 2 +B15 V 2 U+B 16 U 3

[0113] +b 17 UW 2 +b 18 V 2 W+b 19 U 2 W+b 20 W 3

[0114] Number(u,v,w)=c1+c2V+c3U+c4W

[0115] +c5VU+c6VW+c7UW+c8V 2

[0116] +c9U 2 +c 10 W 2 +c 11 UVW+c 12 V 3

[0117] +c 13 VU 2 +c 14 VW 2 +c 15 V 2 U+c 16 U 3

[0118] +c 17 UW 2 +c 18 V 2 W+c 19 U 2 W+c 20 W 3

[0119] Denr(u,v,w)=d1+d2V+d3U+d4W

[0120] +d5VU+d6VW+d7UW+d8V 2

[0121] +d9U 2 +d 10 W 2 +d 11 UVW+d 12 V 3

[0122] +d 13 VU2 +d 14 VW 2 +d 15 V 2 U+d 16 U 3

[0123] +d 17 UW 2 +d 18 V 2 W+d 19 U 2 W+d 20 W 3

[0124] Among them, a1~a 20 , b1~b 20 , c1~c 20 , d1~d 20 Collectively referred to as rational polynomial coefficients (RPCs).

[0125] Combining the above equations, we can obtain the original image point coordinates (C, R):

[0126]

[0127] Assemble roof point clouds After the above RPC model is converted into a two-dimensional point set of the roof After that, you need to work with the building outline vector collection Complete one-to-one matching. Find the center of the roof

[0128]

[0129] Among them, a represents the roof The a-th point in (c a ,r a ) represents the coordinates of the a-th standardized image point.

[0130] Find the center of a vector

[0131]

[0132] Where b represents a vector The b-th image point in (x b ,y b ) represents the vector coordinates of the b-th image point.

[0133] Based on the vector center point set, search for the nearest l candidate vectors for each roof center in turn Represents the building outline vector set closest to the center of the roof point cloud and performs polygon encirclement detection:

[0134]

[0135] Among them, I(·) represents the characteristic function, which determines the planar roof point cloud Is the g-th image point in the vector V l In (c g ,r g ) represents the coordinates of the g-th normalized image point, and calculates the overlap rate between the building outline vector and the center of the roof point cloud:

[0136]

[0137] in, Point cloud for the roof The number of points, S represents the set of points within the building outline vector, |S| represents the roof point cloud The number of points within the building outline vector. This method selects the vector with the best overlap as the best match:

[0138]

[0139] Among them, α is the minimum overlap rate threshold, that is, under the premise of meeting the minimum overlap, select The vector constraint solution with the highest overlap rate, Represents the building outline vector that has the highest overlap with the center of the roof point cloud.

[0140] According to the algorithm process, each roof point cloud that meets the overlap threshold is processed in turn, and finally N optimized roof clusters are obtained. After the roof is constrained, the generated DSM effect is as follows Figure 7 shown.

[0141] S5. Repair ground holes based on grid projection to obtain the final point cloud set

[0142] The relevant operations of step S5 are introduced with a specific example:

[0143] After the point cloud optimization of the building area is completed, due to the removal of the facade and part of the roof, these places will lose points on the xOy plane perpendicular to the z axis, resulting in holes in the generated DSM and other 3D products. With roof point cloud Project it onto the xOy plane and discretize it using a regular grid to detect the parts with missing points.

[0144] The first step is to divide the xOy plane into a regular grid. The grid index (c, d) is calculated as follows:

[0145]

[0146] (x, y, z) is the coordinate of the point cloud in the UTM coordinate system, Δ represents the length of the grid unit, x min Indicates the minimum value of the x-axis component of the point cloud in the UTM coordinate system, min Indicates the minimum value of the y-axis component of the point cloud in the UTM coordinate system. When calculating the grid index number for each point in turn, a mask is generated synchronously, marking the grid where the point exists as true. After the traversal is completed, all grid points marked as false can be retrieved in reverse order and determined to be holes. These empty grid points are calculated back to the UTM coordinate system based on the index number:

[0147] x h =c*Δ+x min ,y h =d*Δ+y min

[0148] Prioritize searching with radius r1 (x h ,y h )'s neighborhood point set

[0149]

[0150] in, Indicates p i 's coordinates.

[0151] like If the number of points is less than the neighborhood point set threshold β, the search radius is expanded to r2 to further search the neighborhood point set. For each hole (x h ,y h ) to perform such a neighborhood search and perform interpolation operations based on the neighborhood:

[0152]

[0153] Traversing the entire grid can obtain the coordinates of all hole points (x h ,y h ,z h ), merged into the ground point cloud In the process, a complete ground point cloud is obtained Finally, the ground point cloud and roof point cloud Take the union ∪ to get the final point cloud Effects such as Figure 8 shown.

[0154] For a clearer and more beautiful visualization, the effect comparison before and after the optimization of the point cloud is presented in the form of a 3D model white model. This embodiment conducts experiments in two regions, and the experimental results are as follows: Figure 9 、 10 , 11 and 12. The simulation results show that after the method of the present invention is used, a large number of error points in the point cloud are cleaned and optimized, which has practical value.

Claims

1. A satellite photogrammetry point cloud optimization method for large-scale building areas, characterized by: The method comprises the following steps: S1. Generate photogrammetric point cloud based on stereo image pairs and will Segmentation into ground point sets and non-ground point sets S2, yes Euclidean clustering is used to achieve building individualization; S3. Remove the facade based on the independent building object to obtain the roof point cloud Will After Euclidean clustering, Z roof point cloud sets are obtained S4, constraining rooftop redundant points based on building outline vector; S5. Repair ground holes based on grid projection to obtain the final point cloud set 2. The satellite photogrammetry point cloud optimization method for large-scale building areas according to claim 1 is characterized in that: The said Segmentation into ground points and non-ground points This is done through the cloth filtering algorithm.

3. The satellite photogrammetry point cloud optimization method for large-scale building areas according to claim 2 is characterized in that: The pair The specific method of using Euclidean clustering to realize building monomerization is as follows: S31, will The initial state of is marked as "unvisited"; S32, will All non-ground points in are clustered into different clusters C k , and mark the status of the clustered non-ground point set as "visited"; S33, iterative update C k , obtain M independent building objects Ω C ={C0,C1,…,C M }, k = [0, 1,…, M].

4. The satellite photogrammetry point cloud optimization method for large-scale building areas according to claim 3 is characterized in that: The iterative update C k Specifically: Among them, t represents the number of iterative updates, p j represents the jth non-ground point cloud, Indicates p j The point set within the neighborhood radius ∈; when When Ω is obtained, the iteration ends. C .

5. The satellite photogrammetry point cloud optimization method for large-scale building areas according to claim 4 is characterized in that: Ω C The constraints need to be met: Ω C ={C k |N min ≤|C k |≤N max }, where |C k | represents cluster C k The number of points, N min Indicates C k The minimum number of points, N max Indicates C k The maximum number of points.

6. The satellite photogrammetry point cloud optimization method for large-scale building areas according to claim 5, characterized in that: The facade removal based on independent building objects is performed by the formula: Calculate, where [m s ,m e ] represents the interval index of the roof slice, Represents the final roof interval, m represents the mth C k Slice, p i represents the i-th non-ground point cloud, z i Indicates p i Elevation slice in the z-axis direction, z min Indicates p i The minimum distribution range in the z-axis direction, Δz represents p i The difference between the maximum and minimum distribution ranges in the z-axis direction.

7. The satellite photogrammetry point cloud optimization method for large-scale building areas according to claim 6, characterized in that: The roof redundancy points are specifically constrained based on the building outline vector: For constraints, optimize Get N optimized roof point cloud sets Among them, ρ represents the overlap rate between the building outline vector and the center of the roof point cloud, α represents the minimum overlap rate threshold, Represents the set of building outline vectors closest to the center of the roof point cloud, Represents the building outline vector that has the highest overlap with the center of the roof point cloud.

8. The satellite photogrammetry point cloud optimization method for large-scale building areas according to claim 7, characterized in that: The method of repairing ground holes based on grid projection is specifically as follows: S81, will and Projected onto the xOy plane and divided into a regular grid; S82, traverse the regular grid to obtain the coordinates of the hole point (x h ,y h ,z j ), (x h ,y h ,z h ) Get a complete ground point cloud set S83, will and Take the union and get

Images