Laser-based building engineering flatness measuring method

Abstract

The invention provides a laser-based building engineering flatness measurement method, which comprises the following steps: acquiring initial position data and elevation information of each floor measurement point of a building, and generating a preliminary data set containing three-dimensional coordinates and elevation values; according to the preliminary data set, generating a station distribution matrix by adopting multi-station cooperative measurement; establishing a unified coordinate reference system by adopting an adjustment algorithm based on the station distribution matrix; the local coordinates of all floors are converted into an overall coordinate system through coordinate transformation, and a multi-layer measurement data set is generated; according to the multi-layer measurement data set, calculating an inter-floor elevation difference, and generating an absolute elevation reference value; generating floor flatness distribution data by adopting gridding sampling; and through fitting the plane equation, generating a building three-dimensional flatness distribution data set, calculating a flatness qualified rate and a comprehensive evaluation index, and outputting a flatness quality report.

Description

Technical Field

[0001] The present invention relates to the field of information technology, and in particular to a laser-based construction engineering flatness measurement method. Background Art

[0002] Flatness measurement in construction projects is a key technical step in ensuring the quality and safety of building structures, directly impacting the overall stability and service life of a building. As modern buildings grow taller and more complex, the requirements for flatness measurement accuracy and efficiency are becoming increasingly stringent, making technological innovation in this field particularly important.

[0003] Traditional methods for measuring building flatness primarily rely on conventional measuring equipment such as levels and total stations. These methods have significant shortcomings when faced with complex building environments. Measurement accuracy is significantly affected by environmental factors, the operation process is cumbersome and time-consuming, and it is difficult to achieve continuous measurement over a large area. This is particularly evident in vertical measurements of high-rise buildings, which exhibit significant limitations. The core challenge currently facing construction project flatness measurement stems from the technical complexity of multi-site collaborative measurement. In high-rise buildings or large construction projects, a single measurement site cannot cover the entire building, and multiple measurement sites are required to work together to obtain complete flatness data.

[0004] However, data coordination and synchronization between multiple sites present technical challenges, and the data obtained at each measurement point often suffers from inconsistent benchmarks. This inconsistent benchmark further complicates vertical through-the-building measurement. Specifically, in the entire measurement process, from the building foundation to the roof, the flatness data between floors lacks coherence and consistency. When the flatness data between building floors cannot be effectively connected, it is difficult to establish a complete building spatial coordinate system, which in turn affects the accurate assessment of the overall building quality. Therefore, how to build a laser measurement network that can achieve multi-site collaboration, establish a unified elevation benchmark, and ensure through-the-building flatness measurement from the foundation to the roof has become a key issue that needs to be addressed in laser-based construction project flatness measurement methods. Summary of the Invention

[0005] The present invention provides a laser-based construction engineering flatness measurement method, which mainly includes: Obtain the initial position data and elevation information of the measurement points on each floor of the building to generate a preliminary data set containing three-dimensional coordinates and elevation values; based on the preliminary data set, use multi-site collaborative measurement to generate a site distribution matrix; based on the site distribution matrix, use the adjustment algorithm to establish a unified coordinate reference system; use coordinate transformation to convert the local coordinates of each floor into the global coordinate system to generate a multi-layer measurement data set; based on the multi-layer measurement data set, calculate the elevation difference between floors and generate an absolute elevation reference value; use grid sampling to generate floor flatness distribution data; generate a three-dimensional building flatness distribution data set by fitting the plane equation, calculate the flatness pass rate and comprehensive evaluation indicators, and output a flatness quality report.

[0006] Furthermore, the multi-site collaborative measurement is used to generate a site distribution matrix, including: counting the number of measurement points in the preliminary data set, and if the number exceeds a preset threshold, using the Delaunay triangulation algorithm to calculate the rangefinder deployment position; synchronously collecting the three-dimensional coordinates and elevation values of each site to generate a sub-dataset; using the highest-precision measurement station as the base station, using the ICP algorithm to convert the sub-dataset to the base station coordinate system to generate a site distribution matrix; calculating the standard deviation of the coordinates of each point in the matrix, marking and reviewing outliers, and outputting the final measurement data set.

[0007] Furthermore, the unified coordinate reference system is established by using an adjustment algorithm, including: extracting three-dimensional coordinates from a site distribution matrix, and calculating relative position deviations between sites using a least squares adjustment algorithm; constructing a geometric constraint equation based on a triangle closure condition, and calculating a rotation matrix and a translation vector using SVD decomposition; converting local coordinates into a unified coordinate system, calculating the distance error between adjacent sites, and if the error exceeds a preset threshold, adjusting the site position using the LM algorithm, and outputting optimized coordinate reference data.

[0008] Furthermore, the local coordinates of each floor are converted into a global coordinate system through coordinate transformation, including: extracting local coordinates from laser measurement data, calculating offsets to generate an initial coordinate array; obtaining the principal axis direction of the point set through principal component analysis, constructing a rotation matrix and a translation vector, and generating a transformation matrix; applying the transformation matrix to transform the coordinates and calculating the standardized elevation value; if the elevation difference exceeds a preset threshold, using least squares optimization to adjust the coordinates and generate structured data including the site number, global coordinates and elevation value.

[0009] Furthermore, the calculation of the elevation difference between floors and the generation of an absolute elevation reference value includes: extracting adjacent floor point clouds from multi-layer measurement data, and using the ICP algorithm to align and calculate the elevation difference; calculating weights based on point cloud density and measurement accuracy, and accumulating the elevation difference using a weighted average algorithm; if the cumulative error exceeds a preset threshold, optimizing the elevation difference using the least squares method; detecting and replacing outliers using the DBSCAN algorithm, and generating a structured table containing floor numbers and absolute elevation values.

[0010] Furthermore, the grid sampling is used to generate floor flatness distribution data, including: using a laser plane scanner to perform grid sampling, obtain node coordinates and elevation values, and convert them into an absolute coordinate system; calculating the elevation deviation matrix, and using Kriging interpolation to optimize the deviation; calculating the local flatness index with the node as the center, and using the DBSCAN algorithm to detect abnormal coordinates; replacing abnormal values with the average value of the elevation deviation of adjacent nodes, and generating a flatness distribution matrix containing node coordinates and elevation deviations.

[0011] Furthermore, the generated site distribution matrix is specifically an M×4 site distribution matrix, where M is the total number of measurement points and the 4 columns correspond to the X, Y, and Z coordinates and elevation values, respectively. For each measurement point in the site distribution matrix, the standard deviation σ of its three-dimensional coordinates in multiple measurements is calculated. If the standard deviation of any coordinate component of a point exceeds 3σ, it is marked as a coordinate outlier. After manual review of all outliers, the final building measurement data set is output.

[0012] Furthermore, the relative position deviation between the stations is calculated using a least squares adjustment algorithm based on the three-dimensional coordinates and the elevation value of each measurement point in the station distribution matrix; Establish a set of geometric constraint equations based on triangle closure conditions and solve the unified coordinate transformation parameters; The rotation matrix R and translation vector T satisfy the equation ∑(R·Pi+T-Qi)² to be minimum, Pi is the local coordinate of site i, and Qi is the target coordinate.

[0013] Furthermore, the standardized elevation value difference between adjacent floors is extracted from the multi-layer measurement data set, and an elevation transfer path between floors is established by vertical laser projection; The weighted average elevation transfer algorithm is used to accumulate and calculate the total elevation difference from the base layer to the target layer layer by layer, where the elevation difference ΔH = ∑wi·ΔHi, wi is the weight of the i-th layer, and ΔHi is the elevation difference of the i-th layer.

[0014] Furthermore, according to the floor flatness distribution data matrix, a least squares plane fitting algorithm is used to calculate the best fitting plane equation; The plane equation is z=ax+by+c, a, b, and c are fitting parameters, and the distance deviation from each node to the fitting plane is calculated to obtain the root mean square error of the flatness.

[0015] The technical solution provided by the embodiment of the present invention may have the following beneficial effects: The invention discloses a laser-based method for measuring the flatness of construction projects. A laser rangefinder is used to collect position data of key points on multiple floors, generate a unified spatial coordinate system, and realize high-precision floor elevation transmission.

[0016] Gridded laser scanning is used to obtain floor plane elevation deviations, the least squares method is used to fit the flatness distribution, and an overall flatness model is constructed based on 3D point cloud fusion.

[0017] This invention effectively addresses the accuracy and efficiency challenges of large-scale building flatness measurements through multi-site collaborative measurement, coordinate transformation, and error adjustment. It ultimately generates a visual flatness quality report, including floor pass rates and overall evaluation indicators. This provides comprehensive and accurate data support for construction project quality control, significantly improving the automation and reliability of building flatness measurements. BRIEF DESCRIPTION OF THE DRAWINGS

[0018] Figure 1 The present invention is a flow chart of a laser-based construction engineering flatness measurement method.

[0019] Figure 2 Schematic diagram of a laser-based construction engineering flatness measurement method of the present invention.

[0020] Figure 3 This is another schematic diagram of a laser-based construction engineering flatness measurement method of the present invention.

[0021] Figure 4 This is another schematic diagram of a laser-based construction engineering flatness measurement method of the present invention.

[0022] Figure 5 This is another schematic diagram of a laser-based construction engineering flatness measurement method of the present invention. DETAILED DESCRIPTION

[0023] To help those skilled in the art better understand the technical solutions in this specification, the following will provide a clear and complete description of the technical solutions in the embodiments of this specification, in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of this specification, not all of them. All other embodiments derived by those skilled in the art based on the embodiments in this specification without creative effort shall fall within the scope of protection of this specification.

[0024] like Figure 1-5 In this embodiment, a laser-based construction engineering flatness measurement method may specifically include: In step S101, a laser rangefinder is used to collect initial location data and elevation information for key measurement points on each floor of the building, generating a preliminary dataset containing the site number, three-dimensional X, Y, and Z coordinates, and elevation values. If the number of measurement points exceeds a preset threshold of 100, a multi-site collaborative measurement mode is initiated. Data is collected synchronously by multiple rangefinders and integrated to create a site distribution matrix containing the three-dimensional coordinates and elevation values of all sites.

[0025] Key measurement points on each floor of the building are scanned using a laser rangefinder to obtain an initial dataset containing station numbers, 3D coordinates (X, Y, Z), and elevation information. The number of valid measurement points in the initial dataset is counted. If this number exceeds a pre-set threshold N, the multi-station collaborative measurement process begins. Based on the building floor plan, the Delaunay triangulation algorithm is used to calculate the optimal deployment locations for the rangefinders, ensuring that all key measurement points are within the coverage of at least two rangefinders. K rangefinders are placed according to the deployment plan, and 3D coordinates and elevation information are simultaneously collected at each station to generate K sub-datasets. The station with the highest accuracy is selected as the base station, and the ICP algorithm is used to transform the sub-datasets of the remaining stations into the base station's coordinate system. Spatial interpolation is performed on all transformed data to generate an M×4 station distribution matrix, where M is the total number of measurement points and the four columns correspond to the X, Y, Z coordinates, and elevation values. For each measurement point in the station distribution matrix, the standard deviation σ of its 3D coordinates across multiple measurements is calculated. If the standard deviation of any coordinate component of a point exceeds 3σ, it is marked as a coordinate outlier. After manual review of all outliers, the final building measurement dataset is output.

[0026] For example, when using a laser rangefinder to measure a five-story office building, acquiring an initial dataset is crucial. The laser rangefinder emits a laser beam and receives the reflected signal, calculating the three-dimensional coordinates and elevation information of key points. Assuming 10 key points are selected on each floor, such as wall corners and column centers, a data table containing the site number, X, Y, and Z coordinates, and elevation is generated after scanning. When counting the number of valid measurement points, assuming a total of 60 points, the preset threshold N is 50. Since 60 exceeds 50, the multi-site collaborative measurement process begins. This ensures sufficient data volume and avoids accuracy issues caused by insufficient coverage at a single site.

[0027] In one possible implementation, a Delaunay triangulation algorithm is used to optimize rangefinder deployment based on a building floor plan. The floor plan shows each floor as a rectangular area. The algorithm constructs a triangular mesh and calculates the rangefinder locations so that each key point is covered by at least two devices.

[0028] For example, for 10 key points on the first floor, the algorithm determined that rangefinders should be placed at two diagonal corners of a rectangular area to ensure unobstructed line of sight. This deployment maximizes coverage, reduces measurement blind spots, and improves data reliability.

[0029] Specifically, three rangefinders (K=3) were deployed, named Station A, Station B, and Station C, to synchronously collect data to generate three sub-datasets. Each sub-dataset contains the X, Y, and Z coordinates and elevation values of 60 points. Station A was selected as the base station because its location at the center of the building minimizes signal interference and maximizes accuracy. The ICP algorithm was used to convert the data from Stations B and C to the coordinate system of Station A. ICP iteratively calculates the minimum distance between point clouds and adjusts the coordinates of the remaining stations to achieve a unified coordinate system.

[0030] For example, the initial coordinates of a point at station B are (10.5, 20.3, 3.2), which are adjusted to (10.4, 20.2, 3.1) after conversion, and aligned with station A, eliminating the deviation between stations.

[0031] In one embodiment, spatial interpolation generates a site distribution matrix. Assuming that the total number of measurement points M is 60, the matrix dimension is 60×4, and each row stores the X, Y, Z, and elevation of a point.

[0032] For example, if the data for a point is (15.2, 25.1, 4.0, 3.9), kriging is used for interpolation to ensure that the data is smooth and reflects spatial continuity, improving the accuracy of subsequent analysis.

[0033] For example, the standard deviation σ of the site distribution matrix is calculated to assess measurement stability. For a certain point, the X coordinates of multiple measurements are 15.2, 15.3, and 15.1, with a calculated standard deviation σ of 0.1. If the standard deviation of the X coordinate of a point is 0.4, exceeding 3σ (0.3), it is marked as an outlier. This outlier may be caused by equipment obstruction or irregularities in the reflective surface, requiring manual review. During review, the surveyor checks the environment at the point and discovers that the glass curtain wall is causing the reflection anomaly. Remeasurement is then performed and the data is updated.

[0034] As you can see, marking and reviewing outliers improves the reliability of the dataset, ensuring accurate building measurements for subsequent structural analysis or construction verification. The final output dataset is a calibrated 60×4 matrix containing the precise coordinates and elevations of all key points, providing a reliable basis for building monitoring.

[0035] Step S102 uses the least squares adjustment algorithm to calculate the relative position deviations between sites based on the three-dimensional coordinates and elevation values of each measurement point in the site distribution matrix. A set of geometric constraint equations based on triangle closure conditions is established to solve the unified coordinate transformation parameters, where the rotation matrix R and the translation vector T satisfy the equation ∑(R Pi + T - Qi)² to a minimum, where Pi is the local coordinate of site i and Qi is the target coordinate. If the distance error between adjacent sites is greater than a preset threshold of 5 mm, the site positions are recalibrated, the site distribution matrix is updated, and a unified spatial coordinate reference system is generated.

[0036] Three-dimensional coordinates and elevation data are extracted from the station distribution matrix, and the coordinate data is normalized to eliminate dimensional differences. The normalized coordinates are input into the least-squares adjustment algorithm to construct the observation equation and the design matrix. The normal equations are solved to obtain a set of residual vectors for the relative positions between stations. A triangle closure error equation is established based on the residual vector set, and the modulus of the three-station vector cross product is defined as the closure condition. SVD decomposition is used to calculate the rotation matrix parameters, and the translation vector is solved by combining the barycentric coordinate difference to obtain a coordinate transformation parameter set containing the quaternion rotation and translation components. The coordinate transformation parameter set is applied to transform the local coordinates into a set of transformed coordinate points, and the Euclidean distance error between adjacent stations is calculated. If the error exceeds a preset threshold, the LM algorithm is used to iteratively adjust the station position, updating the station distribution matrix until all distance errors fall below the threshold. The optimized coordinate reference data is then output. Three-dimensional coordinates are extracted from the optimized coordinate reference data, and the standard deviation of each point relative to its neighbors is calculated as a spatial consistency indicator. After removing outliers exceeding three times the mean square error, kriging interpolation is used to complete the data for missing locations to generate the final spatial reference coordinate set. Verify the closure difference equations of all triangles in the final coordinate set, confirm that the vector cross product modulus length is less than the 1e-6 accuracy requirement, and output the verified reference coordinate system data.

[0037] For example, when using a laser rangefinder to process the three-dimensional coordinates and elevation data of a high-rise building, processing the station distribution matrix is a critical step. The station distribution matrix contains the three-dimensional coordinates and elevation values of all measured points. For example, the data for a certain point may be 10.2, 15.3, 5.1, and 4.9. Normalization eliminates dimensional differences. Assuming the X coordinate range is 0 to 50 meters, a point X = 10.2 is normalized to 0.204. After normalization, the data is input into the least squares adjustment algorithm to construct the observation value equation. The design matrix is based on the geometric relationships between stations, such as the distance constraint between stations A and B, and generates a matrix containing relative positions. After solving the normal equations, a set of residual vectors is obtained, which reflects the deviation between stations.

[0038] For example, the residual vectors from station A to station B are 0.1, 0.2, and 0.05, indicating that the position deviation is small.

[0039] In one possible implementation, the triangle closure error equation is used to verify the consistency of station locations. For example, using stations A, B, and C as an example, the vector cross product modulus must be close to zero.

[0040] For example, if the coordinates of station A are 10, 10, 5, station B is 12, 11, 5, and station C is 11, 12, 5, the modulus of the cross product is close to 0, satisfying the closure condition. SVD decomposition is used to calculate the rotation matrix, which is combined with the difference in the barycentric coordinates to determine the translation vector, generating the transformation parameters in quaternion form.

[0041] For example, the quaternion represents the rotation of station B relative to station A, and the translation vector is 0.3, 0.2, 0.1. After applying the transformation parameters, the local coordinates are transformed into the unified coordinate system, and a point at station B is transformed from 12.5, 11.3, 5.2 to 12.4, 11.2, 5.1.

[0042] Specifically, Euclidean distance error is used to assess conversion accuracy. Assume the distance error between stations A and B is 0.15 meters. A preset threshold of 0.1 meters triggers an iterative adjustment of the LM algorithm. After adjustment, the error is reduced to 0.08 meters, and the station distribution matrix is updated. Three-dimensional coordinates are extracted from the optimized coordinate reference data, and spatial consistency indicators are calculated.

[0043] For example, if the standard deviation between a point and its neighbors is 0.05 meters, which is less than three times the mean error of 0.15 meters, the point is retained. If the standard deviation of a point is 0.2 meters, it is marked as an anomaly, removed, and then supplemented using kriging interpolation.

[0044] For example, the coordinates of a void point are interpolated from the neighboring points 10.3, 15.4, and 5.0 to produce 10.4, 15.5, and 5.0. The final coordinate set verifies the triangle closure error, with the cross product modulus less than 1e-6, confirming the data accuracy.

[0045] For example, kriging interpolation considers spatial correlation, generating smoothed data based on a weighted average of neighboring points to ensure spatial continuity. The LM algorithm uses iterative optimization to reduce error accumulation and improve coordinate system consistency. These methods collectively improve the accuracy of the site distribution matrix, providing a reliable basis for building monitoring.

[0046] Step S103: According to the unified spatial coordinate reference system, the local coordinates of the laser measurement data of each floor are converted into the overall coordinate system of the building through a 4x4 coordinate transformation matrix, where the transformation matrix M is composed of a rotation matrix R and a translation vector T. The standardized elevation value of each measurement point is calculated, and the elevation of the base layer is used as the zero point to generate a multi-layer measurement data set, which includes the site number, overall coordinates and elevation value.

[0047] Extract local coordinate points from the laser survey data and construct an array containing the station number and local coordinates. Based on the preset origin coordinates, calculate the 3D coordinate offset of each point to generate an initial coordinate array. Extract the center of gravity coordinates from the initial coordinate array and calculate the translation vector T. Use principal component analysis to determine the principal axis orientations of the point set, decompose them into Euler angles about the Z, Y, and X axes, and construct a rotation matrix R. Combine R and T to generate a 4x4 transformation matrix M. Apply transformation matrix M to transform the local coordinate points into a global coordinate system, outputting a global coordinate array. Extract the X, Y, and Z components from the global coordinate array and, combined with the base layer elevation data, calculate the standardized elevation of each point. If the elevation difference between the elevation and the base layer exceeds a preset threshold, use the elevation difference as the error term and the coordinate offset as the optimization variable. Scipy's least_squares function performs a least-squares optimization to generate an optimized coordinate array. Combine the optimized coordinate array with the station number to construct a dictionary structure, with the station number as the key and the X, Y, and Z coordinates and elevation as the value. Extract all coordinate points from the dictionary and calculate the Euclidean distance between adjacent points. If the distance exceeds a preset threshold, gradient descent iterations are performed using Scipy's minimize function, using the squared distance deviation as the loss function and the coordinate point as the optimization variable. The result is an updated coordinate array. The elevation difference between each point in the updated coordinate array is compared with the base layer. If the difference exceeds the threshold and elevation data is missing, linear interpolation is performed between adjacent points using Scipy's interp1d function to generate the final coordinate dictionary. The output is structured data containing the station number, global coordinates, and elevation values, using the WGS84 datum.

[0048] For example, when processing the coordinates of high-rise building sites using laser measurement data, you can extract local coordinate point sets to construct an array containing the site numbers and coordinates. Suppose a building has 10 measurement points, numbered S1 through S10, with local coordinates in meters. For example, the coordinates of S1 are 10.5, 15.2, 5.0. Based on a preset origin coordinate (e.g., 0, 0, 0), the offsets for each point are calculated to generate an initial coordinate array.

[0049] For example, the offset of S1 is its own coordinates: 10.5, 15.2, 5.0. Next, the barycentric coordinates are calculated from the initial coordinate array. Assuming the barycentrics of the 10 points are 12.0, 14.0, and 5.5, the translation vector T is the negative of the barycentric coordinates: -12.0, -14.0, and -5.5. This method ensures that the point set is based on the barycentric coordinates, reducing the deviation of the coordinate transformation.

[0050] In a possible implementation, the principal axis directions of the point set are obtained through analysis and decomposed into Euler angles.

[0051] For example, the analysis of the principal axis direction of the point set shows a 30-degree rotation around the Z axis, a 10-degree rotation around the Y axis, and a 5-degree rotation around the X axis. This is how the rotation matrix R is constructed. Combining R and T generates a 4x4 transformation matrix M, which is used to transform the local coordinates into the global coordinate system.

[0052] For example, the S1 coordinates 10.5, 15.2, 5.0 may become -1.5, 1.2, -0.5 after transformation. This transformation unifies the coordinate base and facilitates subsequent analysis.

[0053] Specifically, after extracting the X, Y, and Z components from the global coordinate array, the normalized elevation is calculated based on the base layer elevation data. For example, if the base layer elevation is 0 meters, assume the S2 global coordinates are -0.8, 2.1, and 0.2, with an elevation difference of 0.2 meters. If the preset threshold is 0.1 meters, the S2 elevation difference exceeds the threshold and requires optimization. Using Scipy's least_squares function, with the elevation difference as the error term, the coordinate offset is adjusted to generate an optimized coordinate array.

[0054] For example, after S2 optimization, the coordinates may become -0.8, 2.1, and 0.1, and the elevation difference falls within the threshold.

[0055] For example, the optimized coordinate array is combined with the station number to generate a dictionary, where the corresponding values for the key S1 are -1.5, 1.2, -0.5, and 0.0. The coordinates are extracted from the dictionary and the Euclidean distance between adjacent points is calculated. For example, if the distance between S1 and S2 is 2.5 meters, if the preset threshold is 2.0 meters, the squared distance deviation is used as the loss function for any points exceeding the threshold, and Scipy's minimize function is used to iteratively adjust the loss.

[0056] For example, after adjustment, the S2 coordinates become -0.7, 2.0, 0.1, and the distance is reduced to 1.9 meters, which meets the requirements.

[0057] In one possible implementation, the elevation difference is checked for the updated coordinate array. If the S3 elevation difference is 0.3 meters and data is missing, the Scipy interp1d function is used to interpolate the S2 and S4 coordinates to 0.1 meters. The final coordinate dictionary is output in WGS84 datum, containing the site number, global coordinates, and elevation values. This multi-step optimization approach ensures coordinate data accuracy and consistency, providing reliable support for building monitoring.

[0058] Step S104 extracts the standardized elevation differences between adjacent floors from the multi-layer measurement data set. Vertical laser projections are used to establish inter-floor elevation transfer paths. A weighted average elevation transfer algorithm is then used to cumulatively calculate the total elevation difference from the base layer to the target layer, layer by layer. Here, elevation difference ΔH = ∑wi·ΔHi, where wi is the weight of the i-th floor and ΔHi is the elevation difference of the i-th floor. If the cumulative error exceeds a preset threshold of 10 mm, the error is redistributed using the least squares method to generate absolute elevation benchmarks for each floor.

[0059] From the multi-layered point cloud data, point cloud matching surfaces of adjacent floors were extracted. The point clouds of adjacent floors were registered using the ICP algorithm. The coordinate offset after registration was calculated as the raw elevation difference ΔHi, generating an initial elevation difference array. The PDAL tool was used to calculate the point cloud density grid for each floor. Combined with the total station measurement accuracy report, the weight coefficient wi was normalized by multiplying the density and accuracy by the inverse product to generate a weight array. The numpy.average function was used to perform a weighted calculation on the elevation difference array and the weight array, outputting the cumulative elevation difference ΔH. The preset tolerance threshold for elevation difference was 3σ, where σ is the nominal accuracy of the total station. When ΔH exceeded the threshold, the scipy.optimize.lsq_linear solver was used, using the raw measurement error of each layer as a constraint, to optimize and adjust the elevation difference ΔHi, generating an optimized elevation array. The elevation differences of each layer were extracted from the optimized array and added to the known elevation values of the base layer to obtain the absolute elevation reference value in the WGS84 coordinate system. The DBSCAN algorithm was used to detect the elevation differences between adjacent floors. The ε parameter was set to 2 times the mean square error, and discrete points were marked as outliers. For each outlier floor, the arithmetic mean of the elevations of the floors above and below it was linearly replaced to generate the final elevation dataset. The output was a structured table containing the floor number, WGS84 elevation value, and elevation benchmark value.

[0060] In a possible implementation, to extract point cloud matching surfaces of adjacent floors from point cloud data of a multi-story building, the preprocessing quality of the point cloud data must be ensured.

[0061] For example, a high-rise building has 20 floors. A laser scanner captures point cloud data for each floor. The point cloud contains feature points for the floor, ceiling, and walls. When extracting matching surfaces, flat areas, such as floor slabs, are preferred because their geometric features are stable and suitable for registration. Using the plane segmentation algorithm in the PCL library, with a distance threshold of 0.05 meters, a subset of the point cloud for each floor slab is extracted.

[0062] For example, the layer 5 point cloud subset contains 100,000 points, representing the floor surface, and is saved in PCD format.

[0063] Specifically, when using the ICP algorithm to register adjacent floor point clouds, it is necessary to set the initial registration parameters.

[0064] For example, the initial alignment of the 5th and 6th layer point clouds assumes a translation of 0 and a rotation of 0 degrees. The ICP algorithm calculates the optimal transformation matrix between the two layers of point clouds through iterative optimization and outputs a translation vector and a rotation matrix.

[0065] For example, after registration, the translation vectors of the 6th floor relative to the 5th floor are 0.02, 0.01, and 0.15 meters, indicating an elevation offset of 0.15 meters. This is stored in the initial array as the original elevation difference ΔHi. After all floors are registered, an array containing 20 elevation differences is generated.

[0066] For example, when calculating the point cloud density raster, use the PDAL tool to process each layer of the point cloud.

[0067] For example, the fifth layer of point cloud is rasterized with a grid size of 0.1 meters, resulting in a density of 1000 points per square meter. Based on the total station accuracy report, assuming a nominal accuracy of σ of 0.02 meters, the weight coefficient wi is calculated.

[0068] For example, the fifth layer has high density and high precision, so wi is 0.8; the sixth layer has low density and wi is 0.6. After generating the weight array, call the numpy.average function to calculate the weighted elevation difference ΔH.

[0069] For example, the weighted average elevation difference of 20 floors is 0.18 meters.

[0070] In one possible implementation, we check whether ΔH exceeds the threshold. Assume that 3σ is 0.06 meters and ΔH is 0.18 meters, which exceeds the threshold. We use the scipy.optimize.lsq_linear solver to optimize the elevation differences between layers, using the total station measurement error as a constraint.

[0071] For example, the elevation difference of layer 6 is adjusted from 0.15 meters to 0.10 meters to generate an optimized elevation array. The base layer elevation, such as 0 meters, is then added to obtain the absolute elevation in the WGS84 coordinate system.

[0072] For example, the elevation of the 6th floor is 5.10 meters.

[0073] Specifically, the DBSCAN algorithm is used to detect outliers, and ε is set to 2 times the mean error, that is, 0.04 meters.

[0074] For example, the 7th floor, with an elevation difference of 0.20 meters, is marked as a discrete point. The 7th floor outlier is replaced by the average of 5.10 meters on the 6th floor and 5.30 meters on the 8th floor. This generates a structured table containing floor numbers (e.g., F5), WGS84 elevations (e.g., 5.10 meters), and elevation benchmark values. This approach ensures elevation data consistency and supports the accuracy requirements of building monitoring.

[0075] In step S105, based on the absolute elevation reference value of each floor, a laser plane scanner is used to establish a uniformly distributed measurement grid within the floor plane using a 50 cm x 50 cm grid sampling method. Laser ranging is performed on each grid node to obtain the elevation deviation value and generate a two-dimensional distribution data matrix of floor flatness. The matrix records the coordinates and elevation deviation value of each node.

[0076] A laser planar scanner was used to perform fixed grid sampling of the floor plane, obtaining the plane coordinates x, y, and elevation z of each grid node. The coordinates of at least three reference points were measured using a total station. The scanner data was converted to an absolute coordinate system to generate a two-dimensional data matrix with elevations. The elevation z value was extracted from the two-dimensional data matrix and subtracted from the design elevation value to obtain an elevation deviation matrix. The root mean square error (RMS) of the elevation deviation matrix was calculated. If it exceeded a preset threshold of 0.5 mm, Kriging interpolation was performed using the pykrige library. The grid coordinates x, y, and elevation deviation were input and the output was an optimized elevation deviation matrix. Within the optimized elevation deviation matrix, a 3×3 window was created centered at each node. The standard deviation of the elevation deviation within the window was calculated as a local flatness indicator to generate a flatness distribution dataset. The sklearn.cluster.DBSCAN algorithm was used with a neighborhood radius of 0.2 m and a minimum sample size of 3 to detect elevation anomalies. For each elevation anomaly, the elevation deviations of the four adjacent grid nodes were extracted, and the arithmetic mean was calculated to replace the outlier value to generate a flatness correction dataset. Match the corrected data to the original grid coordinates and reconstruct the final elevation deviation distribution matrix using numpy.meshgrid.

[0077] For example, when using a laser plane scanner to perform fixed-grid sampling of floor plans, a high-precision laser scanner can be selected to ensure millimeter-level accuracy for the elevation z values of grid nodes. For example, for a floor slab on a 20-story building, a design grid size of 0.2 m x 0.2 m is used. A 100 x 100 grid matrix is generated, with each node recording the plane coordinates x, y, and elevation z. The scanner scans the floor surface at a fixed frequency, acquiring evenly distributed point cloud data, which is converted into coordinates and elevation values for grid nodes. This approach ensures comprehensive data coverage and is suitable for subsequent elevation deviation analysis.

[0078] In one possible implementation, the scanner data is converted to an absolute coordinate system by measuring the coordinates of at least three reference points using a total station. For example, three fixed columns at the corners of a floor are selected as reference points. Their WGS84 coordinates are (100.00, 200.00, 0.00), (100.00, 210.00, 0.00), and (110.00, 200.00, 0.00) meters, respectively. The total station has a measurement accuracy of 0.01 meters. Using a coordinate conversion algorithm, the scanner's relative coordinate system is aligned to the WGS84 coordinate system, generating a two-dimensional data matrix with absolute coordinates. This method ensures consistent reference of elevation data and facilitates cross-floor comparisons.

[0079] Specifically, the elevation value z is extracted from the two-dimensional data matrix and the design elevation value is subtracted to obtain the elevation deviation matrix. Assuming the design elevation is 0.00 meters and the measured elevation of a node is 0.03 meters, the deviation is 0.03 meters. If the elevation deviation matrix of the entire floor shows a root mean square error of 0.06 mm, which exceeds the preset threshold of 0.05 mm, further optimization is required. Kriging interpolation is performed using the pykrige library, and the grid coordinates xy and elevation deviation are input to generate a smooth elevation deviation matrix. This interpolation method uses spatial autocorrelation to reduce local elevation mutations and improve data smoothness.

[0080] For example, in the optimization of the elevation deviation matrix, a 3×3 window is established with each node as the center, and the standard deviation of the elevation deviation in the window is calculated as the local flatness index. Assuming that the elevation deviation in a node window is 0.02, 0.03, 0.04 mm, etc., and the standard deviation is 0.01 mm, it means that the local flatness is relatively high. If the standard deviation of a certain area reaches 0.05 mm, there may be construction defects. After generating the flatness distribution dataset, sklearn.cluster.DBSCAN is used to detect outliers, setting the neighborhood radius to 0.2 meters and the minimum number of samples to 3. If the elevation deviation of a node is 0.10 mm, it is far greater than the surrounding nodes and is marked as an anomaly. The deviations of its four adjacent nodes, such as 0.02, 0.03, 0.02, and 0.04 mm, are extracted, and the average value of 0.0275 mm is calculated to replace the outlier. This method effectively eliminates outliers and improves data reliability.

[0081] In one possible implementation, the corrected flatness data is matched to the original grid coordinates, and the final elevation deviation distribution matrix is reconstructed using numpy.meshgrid. For example, if a 100×100 deviation matrix is generated for a particular floor, and the deviation in a certain area is concentrated around 0.03 mm, indicating good flatness, this matrix can be used for building quality assessment and support subsequent construction adjustments.

[0082] In step S106, a least-squares plane fitting algorithm is used to calculate the best-fitting plane equation based on the two-dimensional floor flatness distribution data matrix. The plane equation is z = ax + by + c, with a, b, and c as fitting parameters. The distance deviation from each node to the fitting plane is calculated to obtain the root mean square error (RMS) of the flatness. If the error exceeds the preset building code threshold of 3 mm, the node area is marked as unqualified, and a floor flatness assessment dataset is generated, containing the distribution of qualified and unqualified areas.

[0083] The three-dimensional coordinates and distance deviation values of unqualified nodes were extracted from the flatness assessment dataset to generate a node deviation dataset. Based on the node deviation dataset, the k-nearest neighbor algorithm was used to obtain the neighborhood node set for each unqualified node, using Euclidean distance as the metric and a preset neighborhood radius. For each neighborhood node set, the inverse distance between each node and its neighbors was calculated as a weight. This weight was multiplied by the distance deviation value of the corresponding node, and the sum was divided by the weighted sum to obtain the corrected distance deviation value. If the corrected deviation value exceeded the preset building code threshold, the node coordinates were added to the persistent unqualified node set. Within the persistent unqualified node set, the DBSCAN algorithm was used, setting a neighborhood radius relevant to building codes and a minimum sample size of 5, to cluster unqualified areas. The arithmetic mean of the coordinates of all nodes in each cluster was calculated to obtain the cluster center coordinate set. For the cluster center coordinate set, a square grid with a side length twice the neighborhood radius was generated, with each center as the origin. The distance deviation values of all nodes within the grid were extracted and the arithmetic mean was calculated to generate the flatness correction distribution dataset.

[0084] For example, when extracting the 3D coordinates and distance deviation values of unqualified nodes from a flatness assessment dataset, data generated by a high-precision laser scanner can be used. Assuming a floor is a 20m x 20m rectangular area with a 0.5m x 0.5m grid size, a total of 1,600 nodes are generated. Each node records the 3D coordinates (x, y, z) and distance deviation value. Assuming the design elevation is 0m and the measured elevation of a node is 0.05m, the deviation value is 0.05m. If the building code requires a deviation threshold of 0.03m, nodes with a deviation exceeding 0.03m are marked as unqualified.

[0085] Preferably, the coordinates and deviation values of all unqualified nodes can be extracted through data screening to form a node deviation data set.

[0086] For example, the coordinates of a node are (5.0, 5.0, 0.05) with a deviation of 0.05 meters, and the coordinates of another node are (5.5, 5.0, 0.04) with a deviation of 0.04 meters. These nodes are recorded as unqualified.

[0087] In one possible implementation, the k-nearest neighbor algorithm is used to obtain the neighboring nodes of the unqualified node, with a neighborhood radius of 1 meter and calculation based on Euclidean distance. Assume that the coordinates of an unqualified node are (5.0, 5.0, 0.05), and its neighborhood contains four nodes with coordinates of (4.5, 5.0, 0.02), (5.5, 5.0, 0.04), (5.0, 4.5, 0.03), and (5.0, 5.5, 0.02). The weighted average deviation is calculated using the inverse of the distance as the weight.

[0088] For example, if the distances are 0.5 meters, 0.5 meters, 0.5 meters, and 0.5 meters, respectively, and the weight is 2, the weighted sum of the deviations is divided by the weighted sum, resulting in a corrected deviation of 0.0275 meters. If the deviation still exceeds 0.03 meters after correction, the node is added to the set of persistently unqualified nodes. This method smoothes the deviations through weighted averaging, reducing the impact of local anomalies.

[0089] For example, for a collection of consistently unqualified nodes, the DBSCAN algorithm is applied to cluster unqualified areas, with a neighborhood radius of 1 meter and a minimum sample size of 5. Suppose a region contains six unqualified nodes with coordinates such as (5.0, 5.0, 0.05) and (5.5, 5.0, 0.04), meeting the clustering criteria and forming a cluster of unqualified areas. The average coordinates of the nodes in the cluster are calculated to obtain the cluster center, for example, (5.2, 5.1, 0.045). This clustering method effectively identifies concentrated unqualified areas and facilitates the localization of construction problem areas.

[0090] In one possible implementation, a 2-meter square grid is generated, with the cluster center as the origin, covering the area surrounding the cluster center. For example, if the cluster center is (5.2, 5.1, 0.045), the grid range is (4.2, 4.1) to (6.2, 6.1), containing 16 nodes. The deviation values of these nodes are extracted, and the average deviation, for example, 0.035 meters, is calculated to form a flatness correction distribution dataset. This method, through grid analysis, refines the deviation distribution of unqualified areas, providing data support for subsequent construction optimization.

[0091] It should be noted that the generated flatness correction distribution dataset can be used to evaluate the construction quality of specific areas of a floor.

[0092] For example, an average deviation of 0.035 meters in a certain area indicates a minor flatness issue, which can be corrected through local grinding or filling. This approach ensures that data analysis is closely aligned with actual construction requirements, improving the practicality of floor flatness assessment.

[0093] In step S107, based on the floor flatness assessment dataset, a vertical penetration analysis algorithm based on 3D point cloud fusion is used to construct a 3D flatness model from the foundation to the roof. Kriging spatial interpolation is used to fill in data gaps in uncovered areas, where the interpolation function f(x, y, z) is based on the elevation values of adjacent points. If the flatness variation between adjacent floors exceeds a preset continuity constraint threshold of 5 mm, the outlier data point is corrected to generate a 3D building flatness distribution dataset.

[0094] 3D point cloud data was obtained from a floor flatness assessment dataset. ICP registration was performed using CloudCompare, and the point clouds of each floor were integrated to generate a unified 3D point cloud dataset for the building. Based on this dataset, the standard deviation of elevation within a 10 cm × 10 cm grid was calculated using a moving window method. A flatness model from the foundation to the roof was constructed, and a preliminary flatness distribution was generated. For areas with missing data in the preliminary distribution, elevation values were filled in using the ordinary kriging interpolation module in the GDAL library, based on an exponential variogram model, to generate a continuous flatness distribution. Elevation differences between adjacent floors were extracted from the continuous distribution. Points exceeding a preset 5 mm threshold were identified as outliers using the DBSCAN clustering algorithm. For each outlier point, the elevation difference between each point and its neighbor within a 20 cm radius was calculated. A weighted average correction was performed using the inverse of the Euclidean distance, and a set of corrected elevation values was generated. Based on the correction results, kriging interpolation was re-performed to update the continuous flatness distribution. Data for each floor was extracted from the updated distribution, and a polynomial regression fit was used to fit the inter-floor elevation difference curve. This generated a dataset for analyzing the overall building flatness.

[0095] For example, when acquiring 3D point cloud data based on a floor flatness assessment dataset, a high-precision laser scanner can be used to scan a building and generate point cloud data. Consider a five-story building, each measuring 30 meters by 30 meters, with a point cloud density of 100 points per square meter. CloudCompare is used to perform ICP registration to integrate the point cloud data from each floor. During the registration process, feature points in each layer's point cloud, such as wall corners and pillar edges, are selected to ensure millimeter-level point cloud alignment accuracy.

[0096] For example, the initial deviation between two point clouds is 2 cm. After ICP registration, the deviation is reduced to 1 mm, forming a unified building point cloud dataset. This method ensures the spatial consistency of the point cloud data and facilitates subsequent analysis.

[0097] In one possible implementation, a moving window method is used to calculate the standard deviation of elevation within a 10 cm x 10 cm grid. Assume that a grid on a particular floor contains 100 points, with elevation values ranging from -0.02 m to 0.03 m, and the calculated standard deviation is 0.015 m. By traversing all grids, a preliminary flatness distribution map is generated, reflecting areas of elevation fluctuation.

[0098] For example, if the standard deviation of a certain area is 0.02 meters, it indicates poor flatness and possible construction defects. This method can intuitively reflect the local flatness differences.

[0099] It is understandable that for areas with missing data, ordinary kriging interpolation using the GDAL library is used to fill in the elevation values. Assuming that 10% of the point cloud data is missing in a certain area on a certain floor, the elevation values of the neighboring points are interpolated based on the exponential variation function.

[0100] For example, if the coordinates of a missing point are 10.0, 10.0, and the elevations of the neighboring points are 0.01m and 0.02m, the interpolated elevation is 0.015m. This interpolation generates a continuous flatness distribution, ensuring data integrity.

[0101] For example, when extracting the elevation difference between adjacent floors, suppose the elevations of the corresponding grid points on two floors are 0.05 meters and 0.01 meters, respectively. The difference is 0.04 meters, exceeding the 5 mm threshold. The DBSCAN clustering algorithm is used with a neighborhood radius of 20 cm and a minimum sample size of 5 to mark the outlier points.

[0102] For example, if there are six outliers in a certain area with coordinates such as 10.0, 10.0, and 0.05, they are clustered to form an outlier cluster. This method effectively identifies concentrated outlier areas.

[0103] In one possible implementation, for a set of outliers, the elevation difference within a 20-cm neighborhood is calculated, with the inverse of the Euclidean distance as the weight.

[0104] For example, if the elevation of an outlier point is 0.05 meters, and the elevations of four neighboring points are 0.02 meters and 0.03 meters, respectively, and the distance between them is 15 centimeters, the weighting is 1 / 0.15. The weighted average of the corrected elevation is 0.035 meters. Re-run kriging interpolation to update the continuous flatness distribution and ensure the accuracy of the correction.

[0105] For example, based on the updated distribution, a polynomial regression is used to fit a curve of inter-floor elevation differences. For example, if a quadratic polynomial is fitted to inter-floor elevation difference data, the curve shows that the elevation difference in a certain area gradually increases from 0.01 to 0.04 meters, indicating possible structural settlement. The output data set for overall building flatness analysis includes the distribution of elevation differences for each floor, providing a basis for construction quality assessment. This method clearly reflects the overall building flatness trend.

[0106] Step S108: Based on the three-dimensional flatness distribution dataset of the building, the flatness qualification rate of each floor is calculated, which is defined as the ratio of the qualified area area to the total area, and the overall building flatness comprehensive evaluation index, which is defined as the weighted average of the root mean square error of all floors. A report on the overall flatness quality of the building is generated, using color coding. Green represents qualified areas, red represents unqualified areas, and yellow represents critical areas where the error is close to the 3 mm threshold. A visual building project flatness measurement result dataset is generated.

[0107] The elevation data for each floor is obtained from the three-dimensional building flatness distribution dataset. Using a 1m x 1m grid, each floor is divided into fixed-size grid cells. The mean elevation of each grid cell is calculated to obtain the floor elevation distribution. Based on the floor elevation distribution, the deviation of each grid cell's elevation from the mean elevation is calculated. If the deviation is less than 5 mm, the area is marked as qualified. If the deviation is equal to 5 mm, the area is marked as critical. If the deviation is greater than 5 mm, the area is marked as unqualified, resulting in the regional division results. Based on the regional division results, the number of grid cells with qualified areas on each floor is counted, and the ratio of the area of the qualified area to the total area is calculated to obtain the flatness pass rate for each floor. Based on the flatness pass rate of each floor, a weighted average method is used, with the weight being the current floor area divided by the total building area, to generate a comprehensive evaluation index for the overall building flatness. Based on the regional division results, a color-coded method is used, with qualified areas colored green, critical areas colored yellow, and unqualified areas colored red, to generate a visual flatness distribution dataset.

[0108] For example, when obtaining floor elevation data from a building's 3D flatness distribution dataset, a high-precision laser scanner can be used to generate point cloud data, and point cloud processing software can then be used to extract each floor's elevation information. For example, suppose a building has 10 floors, each measuring 50 meters by 50 meters. After preprocessing the point cloud data, an elevation dataset is generated. Using a 1 meter by 1 meter grid, each floor can be divided into 2500 grid cells, each containing sufficient point cloud data to calculate the mean elevation.

[0109] For example, the elevation values of the point cloud within a grid on a certain floor range from 0.01 to 0.03 meters, with a calculated mean of 0.02 meters. This gridding method facilitates subsequent statistical analysis and is suitable for flatness assessment of large floors.

[0110] In one possible implementation, the deviation of each grid cell's elevation from the average floor elevation is calculated. For example, if the average floor elevation is 0.025 meters and the mean elevation of a grid cell is 0.03 meters, the deviation is 0.005 meters, equal to 5 millimeters, and the area is marked as critical. If the mean elevation of another grid cell is 0.032 meters and the deviation is 0.007 meters, greater than 5 millimeters, the area is marked as unacceptable. This deviation analysis method intuitively reflects local flatness differences and facilitates the identification of problem areas.

[0111] For example, when counting the number of grid cells in the qualified area on each floor, we can traverse all grids and record the number of grid cells with a deviation of less than 5 mm. Assuming a floor has 2,000 grid cells, of which 1,800 have a deviation of less than 5 mm, the qualified area is 1,800 square meters, accounting for 72% of the total area. This statistical method clearly quantifies the flatness pass rate and provides a basis for construction quality assessment.

[0112] In one possible implementation, a weighted average method is used to calculate the overall building flatness evaluation index. For example, assume a 10-story building, each with an area of 2,500 square meters, for a total area of 25,000 square meters. If the pass rate for a particular floor is 80%, the weight is 2,500 / 25,000 = 0.1. The weighted pass rates for each floor are calculated sequentially and summed up to produce the overall index. This method comprehensively considers the differences in floor area, ensuring objective evaluation results.

[0113] For example, when visualizing a flatness distribution dataset, color coding is used. Acceptable areas are colored green, critical areas are colored yellow, and unacceptable areas are colored red. For example, suppose there are 10 unacceptable meshes in a certain area of a floor, concentrated in a corner. Marking them in red clearly indicates the construction defect areas. This visualization method makes it easy to quickly locate problem areas and improves analysis efficiency.

[0114] It can be understood that for critical areas, the distribution trend of their elevation deviations can be further analyzed.

[0115] For example, critical areas on a particular floor may be concentrated at the edges, potentially due to uneven formwork support during construction. By combining color coding with deviation analysis, the construction team can optimize processes accordingly. This multi-dimensional approach provides a comprehensive view of the building's flatness, supporting quality control.

[0116] For example, in the extended analysis, the elevation deviation trend between floors can be combined to determine whether there is overall settlement.

[0117] For example, a gradual increase in deviation at the edge of a floor, from 5 mm to 7 mm, could indicate localized settlement risk. By combining meshing and color coding, the construction team can quickly locate and implement reinforcement measures. This analytical approach, progressing from the local to the overall level, ensures a comprehensive assessment.

[0118] The above description is merely a preferred embodiment of one or more embodiments of this specification and is not intended to limit one or more embodiments of this specification. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of one or more embodiments of this specification shall be included in the scope of protection of one or more embodiments of this specification.

Claims

1. A laser-based method for measuring the flatness of a construction project, characterized in that: include: Obtain the initial position data and elevation information of the measurement points on each floor of the building and generate a preliminary data set containing three-dimensional coordinates and elevation values; Based on the preliminary dataset, a site distribution matrix was generated using multi-site collaborative measurements; Based on the station distribution matrix, an adjustment algorithm is used to establish a unified coordinate reference system; The local coordinates of each floor are converted into the global coordinate system through coordinate transformation to generate a multi-layer measurement data set; Based on the multi-layer measurement data set, the elevation difference between floors is calculated to generate the absolute elevation reference value; Grid sampling is used to generate floor flatness distribution data; By fitting the plane equation, a three-dimensional flatness distribution data set of the building is generated, the flatness qualification rate and comprehensive evaluation index are calculated, and a flatness quality report is output.

2. The method according to claim 1, wherein The method of generating a site distribution matrix by adopting multi-site collaborative measurement includes: Count the number of measurement points in the preliminary data set. If the number exceeds a preset threshold, the Delaunay triangulation algorithm is used to calculate the rangefinder deployment position. Synchronously collect the three-dimensional coordinates and elevation values of each site to generate sub-datasets; Taking the station with the highest accuracy as the base station, the ICP algorithm is used to transform the sub-datasets into the base station coordinate system to generate the station distribution matrix; Calculate the standard deviation of the coordinates of each point in the matrix, mark and review the outliers, and output the final measurement data set.

3. The method according to claim 1, wherein The method of establishing a unified coordinate reference system by using an adjustment algorithm includes: Extract three-dimensional coordinates from the station distribution matrix and use the least squares adjustment algorithm to calculate the relative position deviation between stations; Construct geometric constraint equations based on triangle closure conditions and use SVD decomposition to calculate rotation matrix and translation vector; The local coordinates are converted into a unified coordinate system, and the distance error between adjacent sites is calculated. If the error exceeds the preset threshold, the LM algorithm is used to adjust the site position and output the optimized coordinate benchmark data.

4. The method according to claim 1, wherein The method of converting the local coordinates of each floor into the global coordinate system through coordinate transformation includes: Extract local coordinates from laser measurement data and calculate offsets to generate an initial coordinate array; Obtain the principal axis direction of the point set through principal component analysis, construct the rotation matrix and translation vector, and generate the transformation matrix; Apply the transformation matrix to transform the coordinates and calculate the normalized elevation value; If the elevation difference exceeds the preset threshold, the coordinates are adjusted using least squares optimization to generate structured data containing the site number, overall coordinates, and elevation values.

5. The method according to claim 1, wherein The step of calculating the elevation difference between floors and generating an absolute elevation reference value includes: Extract adjacent floor point clouds from multi-layer measurement data and use ICP algorithm to calculate elevation difference; The weights are calculated based on the point cloud density and measurement accuracy, and the elevation differences are accumulated using a weighted average algorithm; If the cumulative error exceeds the preset threshold, the elevation difference is optimized using the least squares method; Outliers are detected and replaced using the DBSCAN algorithm to generate a structured table containing floor numbers and absolute elevation values.

6. The method according to claim 1, wherein The generation of floor flatness distribution data using grid sampling includes: A laser plane scanner is used to perform grid sampling, obtain node coordinates and elevation values, and convert them into an absolute coordinate system; Calculate the elevation deviation matrix and use Kriging interpolation to optimize the deviation; The local flatness index is calculated with the node as the center, and the DBSCAN algorithm is used to detect abnormal coordinates; The outliers are replaced by the average of the elevation deviations of adjacent nodes to generate a flatness distribution matrix containing node coordinates and elevation deviations.

7. The method according to claim 2, wherein The generated site distribution matrix is specifically an M×4 site distribution matrix, where M is the total number of measurement points and the 4 columns correspond to the X, Y, Z coordinates and elevation values respectively; For each measurement point in the site distribution matrix, calculate the standard deviation σ of its three-dimensional coordinates in multiple measurements; If the standard deviation of any coordinate component of a point exceeds 3σ, it is marked as a coordinate outlier; After manual review of all outliers, the final building measurement dataset is output.

8. The method according to claim 7, wherein Calculating relative position deviations between sites using a least squares adjustment algorithm based on the three-dimensional coordinates and elevation values of each measuring point in the site distribution matrix; Establish a set of geometric constraint equations based on triangle closure conditions and solve the unified coordinate transformation parameters; The rotation matrix R and translation vector T satisfy the equation ∑(R·Pi+T-Qi)² to be minimum, Pi is the local coordinate of site i, and Qi is the target coordinate.

9. The method according to claim 5, wherein The method further includes extracting standardized elevation value differences between adjacent floors from the multi-layer measurement data set, and establishing an elevation transfer path between the floors by vertical laser projection; The weighted average elevation transfer algorithm is used to accumulate and calculate the total elevation difference from the base layer to the target layer layer by layer, where the elevation difference ΔH = ∑wi·ΔHi, wi is the weight of the i-th layer, and ΔHi is the elevation difference of the i-th layer.

10. The method according to claim 6, wherein The method further includes calculating a best-fit plane equation using a least squares plane fitting algorithm based on the floor flatness distribution data matrix; The plane equation is z=ax+by+c, a, b, and c are fitting parameters, and the distance deviation from each node to the fitting plane is calculated to obtain the root mean square error of the flatness.

Images