An engineering earthwork volume measurement method based on unmanned aerial vehicle machine vision

Abstract

The application discloses an engineering earthwork measurement method based on unmanned aerial vehicle machine vision, and particularly relates to the technical field of engineering measurement and digital modeling, and comprises the following steps: based on double-phase images obtained by an unmanned aerial vehicle, a time sequence relationship is established through time information and spatial pose information, pixel-level time sequence consistency analysis is first carried out to form a static reliable image dataset; then, multi-view three-dimensional reconstruction is carried out, and a three-dimensional point cloud model is generated by constraining the ground height mutation area, candidate areas of earthwork edges are extracted based on the surface morphology characteristics of the soil body, and finally, the real boundary of the closed earthwork is determined through normal section stability analysis, and the earthwork volume calculation is completed; through the time sequence analysis of the double-phase images, the application eliminates dynamic interference, reconstructs stable point clouds based on the static reliable images, extracts the real earthwork boundary, and accurately closes the cutting area, so that the stability and accuracy of the engineering earthwork calculation are improved, the overall improvement of the complex terrain applicability and the engineering measurement reliability level is enhanced.

Description

Technical Field

[0001] This invention relates to the field of engineering surveying and digital modeling technology, and more specifically, to a method for measuring engineering earthwork volume based on UAV machine vision. Background Technology

[0002] In the construction of civil engineering, municipal construction, mining, and site leveling projects, the measurement of earthwork volume is a crucial technical foundation for construction planning, progress control, and project settlement. With the development of drone technology and machine vision technology, using drones to acquire images of engineering sites and then performing earthwork cutting and filling analysis or volume difference calculation through 3D reconstruction has gradually become a common technique in the field of engineering surveying.

[0003] Existing methods for measuring earthwork volume in engineering projects based on UAV machine vision typically generate a 3D terrain model through single or multiple image acquisitions, then determine the earthwork edges based on surface height differences or geometric abrupt changes, and perform earthwork cut-fill analysis or volume difference calculations on this basis. However, in actual engineering applications, construction site environments are often complex and variable, and earthwork edges are not regular, flat geometric interfaces, but are often located in complex terrain areas with significant local undulations and irregular shapes, making edge determination methods that rely solely on height differences or geometric abrupt changes unreliable and difficult to apply consistently.

[0004] During construction, non-target objects such as construction equipment, temporarily stockpiled gravel and construction waste, weed cover, and various temporary structures are common. These non-target objects tend to adhere to the soil in space, creating visual adhesion in images. This leads to non-soil objects being mistakenly included in the soil model in the 3D reconstruction results. Existing technologies often lack effective constraints on changes in the time dimension, making it difficult to distinguish between the true surface morphology and dynamic disturbances that appear or disappear within a short period. Consequently, the edges of earthworks are blurred, distorted, or incorrectly extended in the 3D model.

[0005] Traditional methods for determining earthwork boundaries are prone to failure, leading to boundary offsets, incomplete boundary closures, or boundaries being mixed with non-target bodies, resulting in systematic deviations during volume calculations. These deviations are particularly pronounced in large-scale projects or those with complex terrain, affecting the accuracy of earthwork volume calculations to misleading construction decisions, causing imbalances between cut and fill, and even posing safety hazards. Therefore, this invention proposes a method for measuring earthwork volume based on UAV machine vision to address these problems. Summary of the Invention

[0006] To achieve the above objectives, the present invention provides the following technical solution:

[0007] A method for measuring engineering earthwork volume based on UAV machine vision includes the following steps:

[0008] The UAV takes at least two consecutive images of the target engineering area along the same preset flight path, and records the corresponding time information and spatial pose information for each frame of the image, forming a dual-temporal image dataset with a continuous temporal relationship.

[0009] Based on the dual-temporal image dataset, a temporal consistency analysis of pixel changes between images at different times is performed. Non-fixed target areas that have undergone positional or morphological changes in the time dimension are removed from the dual-temporal image dataset, resulting in a static and reliable image dataset that only reflects the static state of the surface.

[0010] Multi-view 3D reconstruction is performed using a static reliable image dataset. During the reconstruction process, edge continuity constraints are applied to areas with abrupt changes in surface height, so that the reconstructed 3D point cloud maintains continuous point cloud distribution and stable normal changes in the earthwork edge area, generating a 3D point cloud model containing point-level reliability information.

[0011] In the three-dimensional point cloud model, the point cloud is separated based on the surface morphology features of the soil. By jointly analyzing the surface roughness, particle size characteristics, normal distribution and height variation trend, non-soil targets that are spatially attached to the soil are separated from the soil surface, and continuous earthwork edge candidate regions are extracted based on the separation results.

[0012] Multiple profiles are constructed along the normal direction of the candidate earthwork edge area. Stability analysis is performed on the profile changes at different times. Edge positions with consistent position changes over time are retained, edge offsets caused by local undulating terrain are suppressed, the closed true boundary of earthwork is determined, and the 3D point cloud model is clipped based on the true boundary of earthwork. Earthwork cutting and filling analysis or volume difference calculation is performed in the clipped area, and the corresponding engineering earthwork volume results are output.

[0013] In a preferred embodiment, the camera pose parameters and imaging scale parameters are kept consistent during image acquisition.

[0014] In a preferred embodiment, timing consistency analysis refers to:

[0015] Geometric registration is performed on the dual-temporal image dataset based on temporal and spatial pose information to establish a one-to-one pixel correspondence across temporal phases;

[0016] For the corresponding pixels, calculate the brightness difference, gradient difference, and texture difference, and normalize and fuse each difference to form the pixel variation distribution;

[0017] Based on the pixel variability distribution, and constrained by the distance relationship between spatially adjacent pixels, adjacent pixels whose variability satisfies the preset continuous distribution conditions are merged to form multiple spatially continuous variability regions.

[0018] In a preferred embodiment, excluding non-surface fixed target areas means:

[0019] The centroid position and region contour of the candidate change region are matched in two time phases, and the displacement distance and morphological change amount across time phases are calculated.

[0020] Candidate change areas whose displacement distance or morphological change exceeds the corresponding preset screening threshold are identified as non-fixed target areas on the ground surface.

[0021] Non-surface fixed target areas are mapped to generate image culling masks. Based on the culling masks, the corresponding image areas are removed, and neighboring static pixels are used for interpolation filling and boundary smoothing to form a static reliable image dataset.

[0022] In a preferred embodiment, interpolation filling and boundary smoothing include the following steps:

[0023] Determine the spatial extent of the area to be removed, and select pixels that maintain a static surface state in images at different times around the boundary of the area as neighboring static pixels;

[0024] Based on the spatial distance relationship between the neighboring static pixels and the pixels inside the eliminated region, the pixel values ​​of the neighboring static pixels are weighted linearly calculated, and the calculation result is used as the filling value of the pixels inside the eliminated region.

[0025] Within the boundary zone between the filled area and the original image, the variation range of adjacent pixel values ​​is constrained to ensure a continuous transition of pixel values ​​along the spatial direction, forming a smooth image boundary.

[0026] In a preferred embodiment, the logic for generating a 3D point cloud model containing point-level reliability information is as follows:

[0027] Based on the spatial pose information of images from different perspectives in the static reliable image dataset, the pixel coordinate difference of the same feature point under each perspective is calculated, and the average value of the pixel coordinate difference of each perspective is used as the disparity consistency measure. Feature points that meet the consistency threshold are selected to generate the initial 3D point cloud.

[0028] Calculate the height difference between adjacent points in the initial 3D point cloud, and identify areas where the height difference is greater than a preset height change threshold as areas of abrupt change in surface height.

[0029] In areas with abrupt changes in surface height, the spatial distance and normal angle between adjacent points are calculated. Points whose spatial distance exceeds a preset spacing threshold or whose normal angle exceeds a preset angle threshold are identified as outliers and removed. At the same time, new points are inserted between the retained points at fixed intervals to complete local densification and reconstruction.

[0030] For each 3D point, calculate the average reprojection error at each viewpoint, and count the number of viewpoints at which the 3D point was successfully observed. The weighted result of the average reprojection error and the number of viewpoints is used as point-level reliability information and written into the 3D point cloud model.

[0031] In a preferred embodiment, point cloud separation and extraction of candidate regions for earthwork edges includes:

[0032] Four types of features are calculated for each point using a fixed radius neighborhood: surface roughness, particle size features, normal distribution, and height variation trend. The four types of features are normalized and weighted to obtain the soil similarity score for each point. Points with soil similarity scores below the separation threshold and spatially connected are identified as non-soil targets and removed from the 3D point cloud model.

[0033] Calculate the local point density gradient and normal mutation intensity for the removed soil point set, and mark the points whose point density gradient exceeds the candidate threshold or whose normal mutation intensity exceeds the candidate threshold as edge points;

[0034] Spatial connectivity closure processing is performed on the edge points to output continuous candidate regions for earthwork edges.

[0035] Surface roughness is taken as the average distance from neighboring points to the fitted plane; particle scale features are taken as the equivalent diameter of the neighboring point cloud after spatial connectivity and aggregation; normal distribution is taken as the variance of the angle between the neighboring normals; and height variation trend is taken as the average height difference along the main direction of the neighborhood.

[0036] In a preferred embodiment, determining the true boundary of the earthwork refers to:

[0037] A profile sampling line is established based on the normal direction of each edge point in the candidate area of ​​the earthwork edge, and the intersection sequence of the two-phase three-dimensional point cloud is extracted on the profile sampling line at a fixed step size to form a profile surface.

[0038] For each set of profiles, calculate the difference in normal projection distance and its mean square value at the intersection points of the edges. Intersection points with mean square values ​​lower than the stability threshold are taken as stable edge locations. Median constraints are applied to the stable edge locations along the boundary tangent to suppress the offset caused by local fluctuations.

[0039] Connect the stable edge positions according to spatial adjacency and perform closure verification to obtain the closed earthwork true boundary. Then, cut the 3D point cloud model with the earthwork true boundary to obtain the cut point cloud. Calculate the cut-fill volume or volume difference within the cut point cloud with the design reference plane or historical reference plane as a reference, and output the engineering earthwork volume result.

[0040] In a preferred embodiment, performing a closure check means:

[0041] Connect the stable edge positions sequentially according to spatial order to form the boundary path;

[0042] Calculate the spatial distance between the beginning and end positions along the boundary path and the overall circumferential direction consistency of the path;

[0043] When the distance between the beginning and end of the path is less than a preset distance threshold and the circumferential direction of the path remains consistent, the boundary path is determined to be a closed boundary.

[0044] If the path does not meet the closure condition, the corresponding boundary path is rejected for use in 3D point cloud clipping and volume calculation.

[0045] The technical effects and advantages of this invention are as follows:

[0046] This invention utilizes unmanned aerial vehicles (UAVs) to acquire images of a target engineering area at least twice along a pre-set, identical flight path, constructing a dual-temporal image dataset with a continuous temporal relationship. It then performs temporal consistency analysis on pixel changes between images from different times, effectively identifying and eliminating non-fixed surface target areas that have undergone positional or morphological changes over time. This ensures that the data entering subsequent processing reflects only the static state of the surface. In real-world engineering scenarios, construction equipment, personnel activities, temporary stockpiles, and rapidly changing surface cover are easily misidentified as terrain components, interfering with 3D reconstruction and earthwork volume calculation. This invention, through temporal consistency constraints, eliminates these dynamic factors at the data source, enabling the static, reliable image dataset to accurately reflect surface morphological changes themselves. This avoids misjudging short-term changes as earthwork changes, providing more stable and reliable basic data support for engineering earthwork volume measurement.

[0047] This invention performs multi-view 3D reconstruction based on a static, reliable image dataset. During the reconstruction process, it applies edge continuity constraints to areas with abrupt changes in surface height, ensuring that the generated 3D point cloud maintains continuous point cloud distribution and stable normal changes in the earthwork edge region. Since earthwork edges are often accompanied by significant height abrupt changes and geometric discontinuities, traditional reconstruction methods easily produce problems such as sparse, broken, or directionally disordered point clouds in these areas, leading to blurred or offset boundaries. This invention applies targeted constraints to areas with abrupt changes in surface height during the 3D reconstruction stage, ensuring that the geometric structure of the earthwork edge region is stably represented at the point cloud level. This reduces the amplification effect of reconstruction errors at the edges, providing a more consistent geometric basis for subsequent point cloud-based soil surface morphology analysis and boundary determination.

[0048] This invention separates point clouds based on soil surface morphology features within a 3D point cloud model and combines this with stability analysis of the normal profile of candidate earthwork edge regions. This effectively suppresses edge offsets caused by local undulating terrain, accurately determines the closed true boundary of earthwork, and then performs earthwork cut-and-fill analysis or volume difference calculation. Under complex engineering terrain conditions, irregular surface undulations and local disturbances can easily lead to unstable edge determination, resulting in an overly large or small trimmed area, affecting the reliability of volume calculation results. This invention first extracts continuous candidate earthwork edge regions, then constructs a profile in the normal direction and performs stability analysis over time, retaining only edge positions with consistent location changes, thereby obtaining a closed and reliable true boundary of earthwork. Based on this boundary, the 3D point cloud model is trimmed before performing earthwork cut-and-fill analysis or volume difference calculation, ensuring that the final output of the engineering earthwork volume is highly consistent with the actual soil area, thus enhancing the practical value of engineering earthwork measurement results in construction management, engineering accounting, and decision support. Attached Figure Description

[0049] To facilitate understanding by those skilled in the art, the present invention will be further described below with reference to the accompanying drawings;

[0050] Figure 1 This is a schematic diagram of a method for measuring engineering earthwork volume based on UAV machine vision, as described in this invention. Detailed Implementation

[0051] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative effort are within the scope of protection of the present invention.

[0052] Reference Figure 1 The following examples were obtained:

[0053] Example 1: A method for measuring engineering earthwork volume based on UAV machine vision, comprising the following steps: The UAV performs at least two consecutive image acquisitions of the target engineering area along a preset identical flight path, and records the corresponding time information and spatial pose information for each frame of the acquired image, forming a dual-temporal image dataset with a temporal continuity relationship; its significance lies in ensuring the consistency of observation conditions through identical flight paths and continuous image acquisition, and using time information and spatial pose information to organize the images into a comparable temporal continuous data basis, providing a unified coordinate reference and traceable data source for subsequent cross-temporal pixel correspondence, change judgment and dynamic interference identification.

[0054] By performing temporal consistency analysis on pixel changes between images from different time periods using a dual-temporal image dataset, non-fixed surface target areas that have undergone positional or morphological changes over time are removed from the dual-temporal image dataset, resulting in a static and reliable image dataset that only reflects the static state of the surface. The significance of this is that temporal consistency analysis isolates non-fixed surface target areas that change over time from the images, preventing them from being mistakenly treated as surface features in 3D reconstruction. This reduces false terrain, edge adhesion, and volumetric errors caused by dynamic factors from the source, making the static and reliable image dataset more realistically reflect the static state of the surface and laying a reliable data foundation for the stable generation of subsequent 3D point cloud models.

[0055] Multi-view 3D reconstruction is performed using a static, reliable image dataset. During the reconstruction process, edge continuity constraints are applied to areas with abrupt changes in surface height. This ensures that the reconstructed 3D point cloud maintains continuous point cloud distribution and stable normal variation in the earthwork edge region, generating a 3D point cloud model containing point-level reliability information. The significance lies in improving the stability of the reconstruction input using a static, reliable image dataset. By applying edge continuity constraints to areas with abrupt changes in surface height, it suppresses breaks, misconnections, and noise diffusion at the edges, ensuring continuous point cloud distribution and stable normal variation in the earthwork edge region. This provides a 3D point cloud model with higher geometric quality for subsequent soil surface morphology analysis and edge extraction. Simultaneously, the point-level reliability information provides quantifiable evidence for subsequent separation processing, edge selection, and result reliability assessment.

[0056] In a 3D point cloud model, point clouds are separated based on soil surface morphology features. Through joint analysis of surface roughness, particle size characteristics, normal distribution, and height variation trends, non-soil targets that are spatially attached to the soil are separated from the soil surface. Based on the separation results, continuous candidate regions for soil edges are extracted. The significance of this method lies in addressing the edge adhesion problem caused by the spatial attachment of soil and non-soil targets in complex construction site environments. By using joint analysis of surface roughness, particle size characteristics, normal distribution, and height variation trends, a basis for judging soil morphology is established. Non-soil targets that affect edge determination and volume calculation are separated from the soil surface, reducing the risk of edge distortion or erroneous expansion. Furthermore, continuous candidate regions for soil edges are extracted from the purified soil point cloud, providing a convergent search range and clear edge candidates for subsequent stability analysis.

[0057] Multiple profiles are constructed along the normal direction of the candidate earthwork edge area. Stability analysis is performed on the profile changes at different times, retaining edge positions with consistent changes over time, suppressing edge shifts caused by local undulating terrain, determining the closed true earthwork boundary, and cropping the 3D point cloud model based on the true earthwork boundary. Earthwork cut-fill analysis or volume difference calculation is performed within the cropped area to output the corresponding engineering earthwork volume results. The significance lies in transforming edge determination into a comparable profile change problem by constructing multiple profiles along the normal direction of the candidate earthwork edge area. Stability analysis is then performed on the profile changes at different times to filter out edge drifts caused by local undulating terrain, point cloud noise, or residual interference, retaining edge positions with consistent changes over time to form a closed true earthwork boundary. Cropping the 3D point cloud model with the closed boundary ensures the integrity and clear boundaries of the volume calculation area. Finally, earthwork cut-fill analysis or volume difference calculation is performed within the cropped area to obtain engineering earthwork volume results that are closer to the actual soil range, reducing systematic bias and improving the reliability of engineering decision-making and settlement basis.

[0058] In one embodiment, taking the target engineering area as an example, a rectangular area with a length of 480 meters and a width of 320 meters is taken as the operation boundary. The flight altitude is set to 90 meters, the forward overlap rate is 80%, and the lateral overlap rate is 70%. The flight path is fixed and saved in the form of a waypoint sequence, with the waypoint spacing set to 25 meters and the turning radius set to 12 meters. At the same time, the camera tilt angle is set to 90 degrees vertical top view, the shutter speed is 1 / 1000 second, the ISO is 200, and the focal length is 24 mm to ensure that the imaging geometry at the same waypoint is repeatable. After the flight path is fixed, the flight path is used as the sole source of the flight path for the next two consecutive image acquisitions, thereby providing a repeatable trajectory basis for the UAV to perform at least two consecutive image acquisitions of the target engineering area along the preset same flight path.

[0059] The first round of continuous image acquisition is performed, with time and spatial pose information written frame by frame. After takeoff, the UAV enters a fixed flight path and continuously captures images at a fixed trigger interval of 0.8 seconds, resulting in, for example, 620 frames per round. Each frame is simultaneously written with its capture time information (time resolution no less than 0.01 seconds) and spatial pose information (including 3D position coordinates and attitude angles, with 3D position coordinate accuracy no less than 0.2 meters and attitude angle accuracy no less than 0.5 degrees). Simultaneously, an index table is created for the first round of acquired images in chronological order. The index table fields include image number, time information, spatial pose information, and waypoint number, ensuring that each frame has a traceable spatiotemporal identifier, providing the first temporal data source for the subsequent formation of a dual-temporal image dataset with a continuous temporal relationship.

[0060] After the first round of acquisition, without changing the flight path and imaging geometry, a second round of continuous image acquisition was performed, and the camera attitude parameters and imaging scale parameters were locked for consistency. The second round of acquisition started 90 seconds after the first round, still using the same waypoint sequence, the same flight altitude of 90 meters, the same waypoint spacing of 25 meters, and the same turning radius of 12 meters, with a shooting trigger interval of 0.8 seconds, resulting in 615 frames of images. To maintain consistency of camera attitude parameters and imaging scale parameters during image acquisition, before the second round took off, the camera tilt angle was locked at 90 degrees, the focal length was locked at 24 mm and automatic zoom was disabled, the image resolution was locked at 5472×3648 pixels and resolution auto-switching was disabled, and the in-flight attitude deviation alarm threshold was set to 1.0 degree. When the attitude angle deviated from the threshold, the marker position of that frame was immediately recorded and a reshoot at the same waypoint was triggered. Each frame of the second round of images was also written with time information and spatial pose information, and a second round index table was established in chronological order, so that the two rounds of images formed a one-to-one correspondence in terms of flight path, attitude, and imaging scale, providing a strict prerequisite for subsequent bi-temporal comparison.

[0061] The acquisition results from the two rounds were paired and aggregated according to time and spatial pose to form a temporally continuous bi-temporal image dataset, and consistency verification was performed. The first and second round index tables were paired according to waypoint number first and time information second. The spatial pose information difference of each pair of paired images was verified. The position difference was controlled within 0.5 meters and the attitude angle difference was controlled within 1.0 degree. Paired images exceeding the threshold were removed or replaced by re-enactment. The retained image pairs were written into a unified dataset directory structure. The directory structure was grouped by waypoint number and sorted by time information within the group. After completion, a finished dataset containing 600 image pairs was obtained. This finished dataset directly corresponds to the bi-temporal image dataset with a temporally continuous relationship in terms of organization. By verifying the consistency of attitude angle difference and resolution, the implementation effect of maintaining the consistency of camera attitude parameters and imaging scale parameters was further guaranteed, thus providing a highly reliable input for subsequent temporal consistency analysis and 3D reconstruction.

[0062] In one embodiment, a temporal consistency analysis is performed on pixel changes between images from different times based on a dual-temporal image dataset. First, geometric registration is performed on the dual-temporal image dataset based on temporal and spatial pose information to establish a one-to-one pixel correspondence across temporal phases. Using image pairs under the same waypoint number as processing units, the temporal and spatial pose information of the two images is read. The three-dimensional position coordinates and attitude angles in the spatial pose information are used to construct the geometric mapping relationship between the two images. Reprojection is then performed on a unified reference plane to align the two images at the same scale and orientation. For example, the reference plane resolution is set to 0.03 meters per pixel on the ground, and the two images are resampled to a unified 4096×4096 pixel grid, ensuring that the same ground location falls within the same pixel index range in both images. After geometric registration, a one-to-one pixel correspondence across temporal phases is established for each reference plane pixel, and the validity flag of the correspondence is recorded for subsequent difference calculations to remove pixels exceeding the overlapping area.

[0063] Based on the one-to-one pixel correspondence across time phases, brightness difference, gradient difference, and texture difference are calculated for corresponding pixels, and each difference is normalized and fused to form a pixel variation distribution. For any pixel index, the brightness values ​​of two frames are extracted and the brightness difference is calculated, where the brightness difference is the absolute difference between the brightness values ​​of the two frames. Simultaneously, the gradient magnitude is calculated within a 3×3 neighborhood centered on the pixel, and the gradient difference is the absolute difference between the gradient magnitudes of the two frames. Furthermore, the texture descriptor is calculated within a 7×7 neighborhood centered on the pixel. The texture descriptor can be taken from the contrast index in the local gray-level co-occurrence relationship statistics, and the texture difference is the absolute difference between the contrast indices of the two frames. For example, if the brightness difference of a pixel is 18... If the gradient difference is 6 and the texture difference is 0.12, then the three types of differences are normalized according to preset upper limits: the normalization upper limit for brightness difference is 60, the normalization upper limit for gradient difference is 20, and the normalization upper limit for texture difference is 0.30, resulting in normalization results of 0.30, 0.30, and 0.40, respectively. These are then fused according to weights of 0.4, 0.3, and 0.3, resulting in a pixel variation of 0.33. This process is repeated for all valid pixels to form a pixel variation distribution covering the reference plane.

[0064] Based on the pixel variability distribution, and constrained by the distance relationship between spatially adjacent pixels, adjacent pixels whose variability satisfies the preset continuous distribution conditions are merged to form multiple spatially continuous variability regions. Using the pixel variability distribution as input, the distance relationship between spatially adjacent pixels is set to eight-neighbor adjacency, with a distance threshold of one pixel interval. The preset continuous distribution conditions are set as follows: the difference in variability between adjacent pixels does not exceed 0.08 and the pixel variability is not less than 0.25. When traversing the pixel variability distribution, pixels that satisfy both spatial adjacency and continuous distribution conditions are merged into the same set to form variability region candidates. For example, if 240 pixels in a local area satisfy the conditions of pixel variability not less than 0.25 and adjacent variability difference not exceeding 0.08, these pixels are merged to form a spatially continuous variability region. Simultaneously, variability region candidates with fewer than 30 pixels are eliminated to avoid interference from scattered changes caused by noise in subsequent judgments.

[0065] The structured output of spatially continuous change regions provides input boundaries for subsequent removal of non-surface fixed target areas from the dual-temporal image dataset. For each spatially continuous change region, the region's bounding rectangle, area, and centroid are calculated, and the region's pixel set is written into a change region list. The change region list fields include region number, area, centroid location, bounding rectangle, and pixel index set. For example, the output of change region number 12 shows an area of ​​2.16 square meters, a centroid location at reference plane coordinates (1860, 2145), and a bounding rectangle with a width of 1.8 meters and a height of 1.2 meters, and is written into the corresponding pixel index set. In this way, temporal consistency analysis not only constructs the pixel variability distribution but also elevates the change information from the pixel level to a traceable spatially continuous change region, providing stable and reusable data objects for subsequent processing.

[0066] In one embodiment, to remove non-surface fixed target areas that undergo positional or morphological changes over time from the dual-temporal image dataset and obtain a static, reliable image dataset that only reflects the static state of the surface, the centroid positions and regional contours of candidate change areas in the dual-temporal images are matched, and the displacement distance and morphological changes across time phases are calculated. Using the candidate change areas output from temporal consistency analysis as input, for each candidate change area, a set of regional contour boundary points is extracted from the first temporal image, and the centroid position is calculated. In the second temporal image, a search window is constructed centered on the centroid position, for example, with a side length of 120 pixels, and candidate contours that satisfy the regional contour similarity constraint are extracted within the search window. The regional contour similarity constraint compares the radial distance sequence from the contour points to the centroid, and the average absolute difference of the radial distance sequence is taken as the contour difference value. When the contour difference value is less than 6 pixels, contour matching is completed, and the centroid position of the second temporal image is obtained. The displacement distance across time phases is taken as the Euclidean distance between the two centroid positions. For example, if the centroid of the first time phase is (1860, 2145) and the centroid of the second time phase is (1878, 2136), the displacement distance is 20.1 pixels. The morphological change is taken as the relative rate of change of the area of ​​the two regions. For example, if the area of ​​the first time phase is 240 pixels and the area of ​​the second time phase is 330 pixels, the morphological change is |330−240| / 240=0.375.

[0067] Candidate change regions whose displacement distance or morphological change exceeds the corresponding preset screening threshold are identified as non-fixed surface target regions. To ensure consistency between the judgment criteria and the imaging scale, pixel displacement is first converted into ground displacement based on the imaging scale parameters of the dual-temporal image dataset. For example, if each pixel corresponds to 0.03 meters on the ground, then 20.1 pixels correspond to 0.603 meters. The preset screening threshold simultaneously sets a displacement distance threshold and a morphological change threshold. For example, the displacement distance threshold is set to 0.30 meters, and the morphological change threshold is set to 0.25. When the displacement distance exceeds 0.30 meters or the morphological change exceeds 0.25, the candidate change region is identified as a non-fixed surface target region. For example, in the above example, the displacement distance is 0.603 meters and the morphological change is 0.375, both exceeding the corresponding preset screening threshold. Therefore, the candidate change region is identified as a non-fixed surface target region, and the region number, centroid position, region bounding rectangle, and cross-temporal calculation results are recorded as the basis for generating subsequent image culling masks.

[0068] The image culling mask is generated by mapping non-surface fixed target areas. The corresponding image areas are removed based on the culling mask, and interpolation filling and boundary smoothing are performed using neighboring static pixels. For each non-fixed target area on the ground surface, the set of contour points of the area is first mapped to the image pixel coordinates to generate a binary image culling mask. Pixels inside the mask are marked as 1, and pixels outside the mask are marked as 0. Then, the corresponding image area is removed based on the culling mask. The removal method is to set the value of the pixel inside the mask to invalid and exclude it from the subsequent processing index. On this basis, the neighboring static pixels are determined. The neighboring static pixels are selected as the set of pixels within the ring formed by extending 8 pixels outward from the mask boundary and not covered by the mask. It is required that these pixels are not identified as candidate change areas in both temporal phases. The interpolation filling adopts a weighted linear calculation. The weight is the reciprocal of the distance from the neighboring static pixels to the pixel to be filled and normalized. For example, 12 neighboring static pixels are selected to participate in the filling. After obtaining the filling value, it is written to the invalid mark position. Boundary smoothing processing performs value transition constraints in the boundary band with a width of 4 pixels inside and outside the mask boundary. The value difference between adjacent pixels in the boundary band does not exceed a preset amplitude threshold, for example, the amplitude threshold is set to 8, so as to avoid abrupt texture at the boundary of the filling area.

[0069] When performing interpolation filling, we can take a pixel to be filled within a removed region as an example. Assume the pixel to be filled is located at a point in the image. Twelve neighboring static pixels are selected as references outside the mask boundary. These neighboring static pixels maintain a static surface state in the dual-temporal image and are not marked as changed areas. The spatial distance between each neighboring static pixel and the pixel to be filled is calculated using the square root of the sum of the squares of the differences in planar coordinates. The reciprocal of each spatial distance is used as the initial weight of the corresponding neighboring static pixels, giving higher weight to neighboring static pixels closer to the pixel to be filled in the filling calculation. All initial weights are normalized so that the sum of all weights is one, thus eliminating the scale effect of different numbers or distributions of neighboring static pixels on the calculation results. The pixel value of each neighboring static pixel is multiplied by its corresponding normalized weight and summed to obtain the filling value of the pixel to be filled. This filling value is then written to the position originally marked as invalid. The fill value obtained by the above method is numerically within the range of neighboring static pixels and exhibits a smooth transition with changes in spatial position, thereby ensuring the continuity of brightness and texture between the filled area and the surrounding static surface area.

[0070] After image culling and masking of all non-surface fixed target areas, a static reliable image dataset is formed and its consistency is verified. Each image frame, after removal, interpolation filling, and boundary smoothing, is written into the static reliable image dataset, retaining the same temporal and spatial pose index relationships as the original dual-temporal image dataset, ensuring the static reliable image dataset still maintains temporal continuity. A verification process is also performed on the static reliable image dataset. This verification involves calculating the mean change in the gradient difference between the images before and after filling at the original non-surface fixed target area locations and comparing it to a threshold. For example, a mean change below 0.10 is considered acceptable, preventing the introduction of new high-gradient false boundaries during filling.

[0071] In one embodiment, to achieve multi-view 3D reconstruction using a static reliable image dataset, and to apply edge continuity constraints to areas of abrupt changes in surface height during the reconstruction process, so that the reconstructed 3D point cloud maintains continuous point cloud distribution and stable normal changes in the earthwork edge region, thereby generating a 3D point cloud model containing point-level reliability information, the following steps are performed:

[0072] A multi-view correspondence is established based on the spatial pose information of images from different perspectives in a static reliable image dataset. The pixel coordinate differences of corresponding feature points under each perspective are calculated, and the average value of the pixel coordinate differences of each perspective is used as the disparity consistency measure. Taking the same ground feature being identified as a corresponding feature point in images from three perspectives as an example, a corresponding feature point refers to a feature point that corresponds to the same physical location on the ground in images from different perspectives and can be identified as the same target point in terms of spatial position through image matching. The pixel coordinates of corresponding points in the three perspective images are read respectively, the pixel coordinate difference between any two perspectives is calculated, and the average of all differences is used as the disparity consistency measure; for example, if the differences are 1.6 pixels, 2.1 pixels, and 1.8 pixels, the disparity consistency measure is 1.83 pixels. The corresponding feature points with a disparity consistency measure less than the consistency threshold are retained. The retained corresponding feature points are used to complete 3D intersection with the spatial pose information to obtain an initial 3D point cloud. The initial 3D point cloud can contain 1,200,000 3D points, providing basic geometric input for subsequent height change localization.

[0073] In the initial 3D point cloud, the height difference between adjacent points is calculated, and regions with height differences greater than a preset height change threshold are identified as surface height abrupt change regions. Specifically, for each 3D point, a set of adjacent points is selected within a fixed radius neighborhood, and the height difference between the 3D point and each point in the adjacent point set is calculated. The maximum height difference is taken as the local height difference of the 3D point. For example, if the fixed radius neighborhood is 0.30 meters, the maximum height difference between a certain 3D point and its adjacent points in the neighborhood is 0.42 meters. When the preset height change threshold is 0.25 meters, the 3D point is marked as a height abrupt change point. Further, spatially adjacent height abrupt change points are aggregated to form regions. The several regions obtained after aggregation are the surface height abrupt change regions, thereby converging the scope of edge continuity constraints from the global to local high-risk areas, avoiding unnecessary geometric disturbances to flat surfaces. Aggregation refers to the process of merging multiple points or pixels that are spatially adjacent, have similar attribute characteristics, and meet the continuity condition into a whole region or object according to spatial adjacency.

[0074] Within the area of ​​abrupt change in surface height, calculate the spatial distance and normal angle between adjacent points. Points whose spatial distance exceeds a preset spacing threshold or whose normal angle exceeds a preset angle threshold are identified as outliers and removed. At the same time, new points are inserted between the retained points at fixed intervals to complete the local densification reconstruction. Specifically, for each 3D point within a region of abrupt change in surface height, several nearest neighboring points are selected, and the spatial distance between the 3D point and its neighboring points is calculated. The angle between the normal of the 3D point and the normal of its neighboring points is also calculated. For example, with a preset spacing threshold of 0.20 meters and a preset angle threshold of 25 degrees, when the spatial distance between a point and its nearest neighboring point is 0.31 meters or the angle between their normals is 38 degrees, the point is identified as an outlier and removed from the initial 3D point cloud. After removal, the remaining point cloud undergoes local densification and reconstruction. The densification method involves inserting new points at fixed intervals along the line connecting adjacent remaining points, for example, a fixed interval of 0.05 meters. The height of the inserted new point is obtained by linear interpolation of the heights of the two remaining points, and the normal of the inserted new point is obtained by normalized weighted average of the normals of the two remaining points. This results in a more uniform and continuous distribution of the point cloud in the earthwork edge region, and a smoother and more stable change in normal, meeting the target requirements of edge continuity constraints.

[0075] For each 3D point, calculate the average reprojection error at each viewpoint, and count the number of viewpoints at which the 3D point was successfully observed. The weighted result of the average reprojection error and the number of viewpoints is used as point-level reliability information and written into the 3D point cloud model. Specifically, each 3D point is projected back into the images from each viewpoint participating in the reconstruction using its spatial pose information. The predicted pixel position of the 3D point in each viewpoint image is obtained, and the pixel distance between the predicted pixel position and the pixel position of the corresponding feature point is calculated as the reprojection error. For example, if the reprojection errors of a 3D point in four viewpoints are 0.9 pixels, 1.2 pixels, 1.0 pixels, and 1.4 pixels, the average reprojection error is 1.125 pixels, and the number of viewpoints is 4. The average reprojection error and the number of viewpoints are weighted to synthesize point-level reliability information. For example, according to the principle that "the more viewpoints, the higher the reliability; the smaller the average reprojection error, the higher the reliability," the point-level reliability information is set as the combination result of the number of viewpoints and the average reprojection error and written into the point attribute field of the 3D point cloud model. Finally, a 3D point cloud model containing point-level reliability information is obtained, which provides a filterable and traceable geometric quality basis for subsequent joint analysis of soil surface morphology features and extraction of candidate areas for earthwork edges.

[0076] It should be noted that after generating a 3D point cloud model containing point-level reliability information, this information is used to reflect the geometric stability and observational consistency of each 3D point during multi-view reconstruction, rather than participating in subsequent calculations. When separating the point cloud based on soil surface morphology features, 3D points with point-level reliability information higher than a preset confidence threshold are prioritized for calculations of surface roughness, particle size characteristics, normal distribution, and height variation trends. This reduces the impact of unstable points introduced by reconstruction errors or insufficient viewpoints on soil similarity score calculations. Simultaneously, when extracting candidate earthwork edge regions, point-level reliability information is used to screen edge points, avoiding misclassification of points with low reliability as true edge points. This improves the continuity and geometric consistency of candidate earthwork edge regions, providing a quality-controlled point cloud foundation for subsequent stability determination of true earthwork boundaries and calculation of engineering earthwork volume.

[0077] In one embodiment, to achieve point cloud separation and extract candidate regions for earthwork edges, four types of features are calculated for each point using a fixed-radius neighborhood: surface roughness, particle scale features, normal distribution, and height variation trend. The fixed-radius neighborhood can be 0.30 meters. For any point in the 3D point cloud model, a set of neighborhood points within a 0.30-meter radius is collected. Surface roughness is obtained by performing plane fitting on the neighborhood point set, calculating the vertical distance from each neighborhood point to the fitted plane and averaging it. For example, if a point has 80 neighborhood points with an average distance of 0.012 meters, then the surface roughness of that point is 0.012. For particle scale features, spatial connectivity aggregation is first performed on the neighborhood point set, grouping points with a distance less than 0.05 meters into the same connected cluster. The connected cluster with the most points is selected, and the particle scale features are calculated... The volume envelope of the connected cluster in three-dimensional space is converted into an equivalent diameter. If the equivalent diameter is 0.18 meters, then the particle size feature is 0.18. The normal distribution is obtained by calculating the angle between the normal of each point in the neighborhood and the normal of the center point and calculating the variance. For example, if the variance of the angle is 36 square degrees, then the normal distribution is 36. The height change trend is obtained by determining the main direction of the neighborhood point set and calculating the average height difference along the main direction. For example, if the average height difference along the main direction is 0.22 meters, then the height change trend is 0.22. This provides four types of comparable feature inputs for subsequent soil similarity score calculation.

[0078] For each point, the four types of features are normalized and weighted to obtain a soil similarity score. Points with soil similarity scores below the separation threshold and spatially connected are identified as non-soil targets and removed from the 3D point cloud model. Normalization can be performed using preset upper limits, such as 0.05 meters for surface roughness, 0.50 meters for particle scale features, 100 square degrees for normal distribution, and 0.80 meters for height variation trend. The normalized values ​​are obtained by dividing each of the four features by these upper limits. The weighted summation can be set with weights of 0.30, 0.25, 0.25, and 0.20 to form the soil similarity score. For example, after normalization, the four features are 0.24, 0.36, 0.36, and 0.28, respectively, resulting in a soil similarity score. The score is 0.30×0.24+0.25×0.36+0.25×0.36+0.20×0.28=0.309; the separation threshold can be 0.35. When the soil similarity score is lower than 0.35, the point is marked as a non-soil candidate point. The non-soil candidate points are further spatially connected and grouped. If a connected group contains more than 300 points and the spatial range length exceeds 0.8 meters, the connected group is determined to be a non-soil target and removed from the 3D point cloud model so that the removed point set more concentratedly reflects the soil surface morphology.

[0079] Calculate the local point density gradient and normal mutation intensity for the removed soil point set, and mark the points whose point density gradient exceeds the point density candidate threshold or whose normal mutation intensity exceeds the normal candidate threshold as edge points. Local point density is obtained by dividing the number of points in a neighborhood with a fixed radius by the volume of the neighborhood. For example, if there are 90 points in a neighborhood with a radius of 0.30 meters, the local point density is 90 divided by the volume of a sphere with a radius of 0.30 meters, resulting in 796 points per cubic meter. The local point density gradient is obtained by comparing the average difference between the point density of the center point's neighborhood and the point density of its neighboring neighborhoods. For example, if the average difference is 210 points per cubic meter, the local point density gradient is 210. The normal abrupt change intensity is obtained by calculating the average angle between the normal of the center point and the normal of the neighboring points. For example, if the average angle is 28 degrees, the normal abrupt change intensity is 28. The candidate threshold can be set to 180 points per cubic meter for the point density gradient and 22 degrees for the normal abrupt change intensity. When the point density gradient exceeds 180 or the normal abrupt change intensity exceeds 22, the point is marked as an edge point, thus focusing edge detection on locations where the geometry changes rapidly.

[0080] Spatial connectivity closure processing is performed on edge points to output continuous earthwork edge candidate regions. This process involves first aggregating edge points based on spatial adjacency, merging edge points with a distance of less than 0.10 meters into edge connected segments, and then performing closure correction on these segments based on endpoint distance. When the distance between the first and last endpoints of an edge connected segment is less than 0.30 meters, the first and last endpoints are connected to form a closed boundary. For edge connected segments with breaks, the nearest edge connected segment is selected based on the consistency of the tangential direction at both ends of the break, thus forming a continuous boundary. After closure, the set of edge points enclosed by the closed boundary and its circumscribed range are output as continuous earthwork edge candidate regions. The edge point coordinate sequence is retained for subsequent construction of multiple profiles along the normal direction of the earthwork edge candidate regions, thereby providing stable input for determining the true earthwork boundary based on point cloud separation.

[0081] In one embodiment, to determine the true boundary of a closed earthwork, a profile sampling line is established using the normal direction of each edge point within the earthwork edge candidate region. The intersection sequence of dual-temporal 3D point clouds is then extracted from the profile sampling line at a fixed step size to form a profile surface. Each edge point in the earthwork edge candidate region is used as the profile center point. The normal direction of that edge point is read, and a profile sampling line passing through that edge point is constructed. The length of the profile sampling line can be 6.0 meters, with 3.0 meters along the positive normal direction and 3.0 meters along the negative normal direction. The fixed step size can be 0.05 meters, resulting in 121 sampling location points for each profile sampling line. For each sampling location point, the nearest point cloud point is searched in both the first and second phase 3D point clouds. The intersection point coordinate sequence is extracted based on the condition that "the projection of the sampling location point to the nearest point cloud point falls near the profile sampling line and the distance is less than 0.08 meters". For example, in a certain profile sampling line, 96 intersection points are obtained in the first phase and 94 intersection points are obtained in the second phase. The two intersection point coordinate sequences are paired according to the order of the sampling location points to form a profile surface, thereby transforming the edge determination problem into a two-phase profile geometric comparison problem, providing standardized input for subsequent stability analysis.

[0082] For each set of profiles, calculate the difference in normal projected distances and its mean square value at the intersection points of the edges. Intersection points with mean square values ​​below the stability threshold are considered stable edge locations. For each set of profiles, project the coordinates of the intersection points onto the normal direction in the profile sampling line coordinate system to obtain the first phase normal projected distance sequence and the second phase normal projected distance sequence. Subtract the two projected distances at the same sampling location to obtain the difference in normal projected distances. Calculate the mean square value of the difference in normal projected distances across the entire profile sampling line range, and take the average of the squared differences. For example, in a certain cross-section, if the mean square value of the difference in normal projection distance is 0.0025 square meters, and the corresponding stability threshold is 0.0064 square meters, then the cross-section meets the stability condition. Further, the edge intersection points within the cross-section can be located. The "location point with the largest gradient of normal projection distance" can be used as the candidate location point for the edge intersection point. The difference in the geometric position of this location point in the two time phases is also included in the mean square value statistics, so that the determination of edge intersection points and stability analysis are mutually constrained, reducing the impact of single-point noise on the edge position. Finally, the candidate location points of edge intersection points with a mean square value lower than the stability threshold are determined as stable edge positions, thereby filtering out edge positions that are greatly affected by dynamic disturbances or local terrain noise in the time dimension.

[0083] The median constraint is applied to the stable edge locations along the boundary tangent to suppress offsets caused by local undulations. After the set of stable edge locations is formed, the boundary tangent direction is determined along the boundary direction of the earthwork edge candidate area. For each stable edge location, five adjacent stable edge locations are selected along the boundary tangent direction to form a neighborhood sequence, and the offset of each stable edge location in the neighborhood sequence in the normal direction is calculated. These offsets are sorted, and the median value is taken as the normal offset correction value of the current stable edge location. The coordinates of the current stable edge location are then updated with the correction value. For example, the normal offsets of the neighborhood sequence of a stable edge position are 0.04 meters, 0.05 meters, 0.07 meters, 0.31 meters, 0.06 meters, 0.05 meters, 0.04 meters, 0.06 meters, 0.05 meters, 0.04 meters, and 0.05 meters, respectively. Among them, 0.31 meters is an abrupt shift caused by local fluctuations. The correction value obtained by median constraint is 0.05 meters, thereby suppressing the abrupt shift shift within the overall consistent level of the boundary tangential neighborhood. Through this processing, the stable edge position shows a continuous change in the boundary tangential direction, avoiding the formation of jagged boundaries and ensuring the geometric stability of subsequent closure verification and pruning.

[0084] The stable edge locations are connected according to spatial adjacency and a closure check is performed to obtain a closed true earthwork boundary. A clipped point cloud is then obtained by trimming the 3D point cloud model using this true earthwork boundary. Within the clipped point cloud, the cut-fill volume or volume difference is calculated with reference to the design datum or historical datum, and the earthwork volume is output. Specifically, the boundary points are first sorted according to the spatial coordinates of the stable edge locations, and then connected sequentially according to the adjacency rule of "the spatial distance between adjacent points is the smallest and less than 0.35 meters" to form boundary paths. After connection, the spatial distance between the first and last points of the boundary path is calculated, and the overall circumferential consistency of the boundary path is also calculated. When the spatial distance between the first and last points is less than 0.30 meters and the circumferential consistency meets the preset consistency condition, a closure check is performed, resulting in a closed true earthwork boundary. Subsequently, the 3D point cloud model is trimmed based on the closed earthwork boundary. The trimming method is to retain the points inside the boundary and discard the points outside the boundary to obtain the trimmed point cloud. A design reference plane or historical reference plane is introduced into the trimmed point cloud as a reference plane. Earthwork cutting and filling analysis or volume difference calculation is performed in the trimmed area to output the corresponding engineering earthwork volume results.

[0085] Closure verification involves first connecting the stable edge positions sequentially according to their relative spatial order after obtaining them, forming a continuous boundary path. Then, the spatial distance between the starting and ending positions is calculated along this boundary path, and the overall orientation of the boundary path is analyzed to determine if its circumferential direction remains consistent. If the spatial distance between the starting and ending positions is less than a preset distance threshold, and the boundary path exhibits no reversal or breakage in its spatial orientation, and its overall circumferential direction remains consistent, the boundary path is considered a closed boundary. If the spatial distance between the starting and ending positions exceeds the preset distance threshold, or the circumferential direction of the boundary path is inconsistent, the boundary path is considered not to have formed a valid closure, and therefore, the corresponding boundary path is rejected for subsequent 3D point cloud trimming and earthwork volume calculation.

[0086] The various thresholds in the embodiments are not set in isolation, but are determined comprehensively based on the overall goal of the present invention to stably and reliably extract the real earthwork boundary and accurately calculate the earthwork volume in dynamic and complex engineering scenarios, taking into account the image imaging conditions, the geometric accuracy of three-dimensional reconstruction, the spatial distribution characteristics of point cloud, and the actual morphological characteristics of the engineering soil. Specifically, the displacement distance threshold and morphological change threshold are used to distinguish non-surface fixed target areas that undergo significant spatial or morphological changes over time. These thresholds are set based on the imaging scale of the dual-temporal images and the movement amplitude of common dynamic targets at the construction site. The consistency threshold is used to filter spatially stable corresponding feature points from multiple perspectives. Its setting is based on image resolution and camera pose calculation accuracy. The height change threshold, spacing threshold, and angle threshold are used to identify areas of abrupt changes in surface height and constrain the continuity of the point cloud at the earthwork edge. These thresholds are set based on the height jump amplitude and normal change range exhibited by the earthwork edge in the real terrain. The separation threshold, along with the point density candidate threshold and normal candidate threshold, are used to distinguish between earthwork and non-earthwork targets and extract edge points with significant geometrical changes. These thresholds are set based on the statistical characteristics of relatively continuous earthwork surfaces and significant morphological differences in non-earthwork targets. The stability threshold and preset distance threshold are used to evaluate the consistency of the profile in the time dimension and the reliability of boundary path closure. These thresholds are set based on the multi-temporal reconstruction error level and the engineering survey requirements for boundary closure accuracy.

[0087] The weighted summation method used in this invention determines the weights based on the distinguishing ability and stability contribution of different features in engineering earthwork identification and boundary determination. Specifically, features that stably reflect the overall soil morphology and are less affected by noise and local disturbances are assigned relatively high weights, while features sensitive to environmental changes or easily affected by local anomalies are assigned relatively low weights. This ensures that the comprehensive result is both discriminative and robust. Furthermore, the specific values ​​of each weight are not fixed but are matched to image resolution, point cloud density, terrain undulation, and engineering measurement accuracy requirements. Without departing from the technical concept of this invention, the specific values ​​or acquisition methods of each weight and corresponding threshold are not limited to the embodiments described in this invention. Other settings in the prior art can also be used, as long as the setting method serves the same technical purpose, achieves the same or similar technical effect, and is applicable to the UAV engineering surveying and earthwork volume calculation environment in which this invention is used. Such settings should be considered to fall within the protection scope of this invention. The specific values ​​or determination methods can be adjusted according to different engineering site conditions, imaging accuracy requirements and measurement specifications, and can be reflected in engineering design documents, engineering implementation plans, engineering technical disclosure documents or engineering measurement operation documents during the engineering implementation process, without affecting the overall technical effect and application value of the present invention in engineering earthwork measurement.

[0088] Performing earthwork cut-fill analysis or volume difference calculation within the clipped area is an existing technique in the field of engineering surveying. It usually uses a design datum or historical datum as a reference and performs integration calculations on the surface elevation data within the defined area to obtain the fill volume, cut volume, or the volume difference between the two. However, in the existing technology, this type of volume calculation is highly dependent on the accuracy of the clipped area boundary and the geometric reliability of the input point cloud. If the boundary is drifted, adhered, or not closed, or if dynamic targets and non-soil interference are mixed in the point cloud, the volume calculation error will often be amplified. This invention does not improve upon the cut-and-fill analysis or volume difference calculation itself. Instead, it provides clear, geometrically stable, and reliable input conditions for existing volume calculations by performing temporal consistency analysis of dual-phase images, constructing a static reliable image dataset, generating a 3D point cloud model containing point-level reliability information, separating soil and non-soil targets, and stably determining the true boundary of the earthwork. This allows earthwork cut-and-fill analysis or volume difference calculation performed within the cropped area to more realistically reflect the actual soil range. Thus, it significantly improves the accuracy and engineering applicability of the earthwork volume results while using existing engineering calculation methods. Therefore, this invention will not elaborate on earthwork cut-and-fill analysis or volume difference calculation.

[0089] The above algorithms or formulas are all dimensionless and numerical calculations, and the results are obtained by software simulation based on a large amount of collected data to obtain the most recent real-world results. The preset parameters are set by those skilled in the art according to the actual situation.

[0090] It should be understood that in the various embodiments of this application, the order of the above-mentioned processes does not imply the order of execution. The execution order of each process should be determined by its function and internal logic, and should not constitute any limitation on the implementation process of the embodiments of this application.

[0091] Those skilled in the art will recognize that the units and algorithm steps of the various examples described in conjunction with the embodiments disclosed herein can be implemented in electronic hardware, or a combination of computer software and electronic hardware. Whether these functions are implemented in hardware or software depends on the specific application and design constraints of the technical solution. Those skilled in the art can use different methods to implement the described functions for each specific application, but such implementation should not be considered beyond the scope of this application.

[0092] Those skilled in the art will clearly understand that, for the sake of convenience and brevity, the specific working processes of the devices and units described above can be referred to the corresponding processes in the foregoing method embodiments, and will not be repeated here.

[0093] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application. Therefore, the scope of protection of this application should be determined by the scope of the claims.

Claims

1. A method for measuring engineering earthwork volume based on UAV machine vision, characterized in that, The method comprises the following steps: The UAV collects images of the target engineering area along the preset same flight route for at least two times in succession, and records corresponding time information and spatial pose information for each frame of collected image, thereby forming a dual-time-phase image data set with time continuity; Based on the dual-time-phase image data set, the pixel change conditions between different time images are analyzed for time sequence consistency, non-ground fixed target regions that have position changes or shape changes in the time dimension are removed from the dual-time-phase image data set, and a static reliable image data set reflecting only the static state of the ground is obtained; The static reliable image data set is used for multi-view three-dimensional reconstruction, and an edge continuity constraint is applied to the ground height mutation area in the reconstruction process, so that the three-dimensional point cloud obtained by the reconstruction keeps the point cloud distribution continuous and the normal change stable in the earthwork edge area, and a three-dimensional point cloud model containing point-level reliability information is generated; In the three-dimensional point cloud model, the point cloud is separated based on the surface shape characteristics of the soil body, the non-soil body target that is spatially attached to the soil body is separated from the soil body surface through joint analysis of the surface roughness, particle size characteristics, normal distribution and height change trend, and a continuous earthwork edge candidate area is extracted based on the separation result; A plurality of profiles are constructed along the normal direction of the earthwork edge candidate area, the profile change conditions corresponding to different times are analyzed for stability, the edge position that changes in position in the time dimension is retained, the edge shift caused by local undulating terrain is suppressed, the closed real earthwork boundary is determined, and the three-dimensional point cloud model is cropped based on the real earthwork boundary, and earthwork cutting and filling analysis or volume difference calculation is performed in the cropped area, and the corresponding engineering earthwork quantity result is output.

2. The method of claim 1, wherein, The camera pose parameters and imaging scale parameters are kept consistent during image collection. 3.The method of claim 2, wherein, The time sequence consistency analysis refers to: Geometric registration is performed on the dual-time-phase image data set according to the time information and spatial pose information, and a one-to-one correspondence relationship of pixels across time phases is established; The brightness difference, gradient difference and texture difference of the corresponding pixels are calculated, and the differences are normalized and fused to form a pixel change degree distribution; Based on the pixel change degree distribution, the adjacent pixels that meet the preset continuous distribution condition in change degree are merged as a constraint of the distance relationship between the adjacent pixels in space, and a plurality of spatially continuous change regions are formed.

4. The method of claim 3, wherein the method further comprises: The removal of non-ground fixed target regions refers to: The centroid position and region contour of the candidate change region in the dual-time-phase are matched, and the displacement distance and shape change amount across time phases are calculated; The candidate change region whose displacement distance or shape change amount exceeds the corresponding preset screening threshold is determined as a non-ground fixed target region; The non-ground fixed target region is mapped to generate an image removal mask, the corresponding image region is removed based on the removal mask, and the adjacent static pixels are used for interpolation filling and boundary smoothing processing to form a static reliable image data set.

5. The method of claim 4, wherein, The interpolation filling and boundary smoothing processing comprises the following steps: The spatial range of the removed region is determined, and the pixels that remain in a static state on the ground in different time images around the boundary of the region are selected as adjacent static pixels; According to the spatial distance relationship between the neighborhood static pixels and the pixels inside the culling region, the pixel values of the neighborhood static pixels are calculated by weighted linear calculation, and the calculation results are taken as the filling values of the pixels inside the culling region; In the junction zone of the filling region and the original image, the variation amplitude of the values of adjacent pixels is constrained, so that the pixel values are continuously transitioned along the spatial direction, and a smooth image boundary is formed.

6. The method of claim 5, wherein the method further comprises: The three-dimensional point cloud model generation logic containing point-level reliability information is as follows: Based on the spatial pose information of different view images in the static reliable image data set, the pixel coordinate difference values of the same feature points under different views are calculated, and the average value of the pixel coordinate difference values under different views is taken as the disparity consistency measure, and the feature points satisfying the consistency threshold are selected to generate an initial three-dimensional point cloud; In the initial three-dimensional point cloud, the height difference between adjacent points is calculated, and the region with a height difference greater than a preset height change threshold is determined as a ground surface height mutation region; In the ground surface height mutation region, the spatial distance and the normal angle between adjacent points are calculated, and the points with a spatial distance greater than a preset distance threshold or a normal angle greater than a preset angle threshold are determined as outlier points and removed, and meanwhile, new points are inserted between the remaining points at a fixed interval to complete local encryption reconstruction; For each three-dimensional point, the average value of the re-projection error under different views is calculated, and the number of views in which the three-dimensional point is successfully observed is counted, and the weighted result of the average value of the re-projection error and the number of views is taken as the point-level reliability information and written into the three-dimensional point cloud model.

7. The method of claim 6, wherein the method further comprises: The point cloud separation and extraction of the earthwork edge candidate region include: Four types of features, surface roughness, particle size feature, normal distribution, and height variation trend, are calculated for each point in a fixed radius neighborhood, the four types of features are normalized and weighted summed to obtain a soil similarity score, and the point set with a soil similarity score lower than a separation threshold and a connected distribution in space is determined as a non-soil object and removed from the three-dimensional point cloud model; The local point density gradient and the normal mutation intensity of the removed soil point set are calculated, and the points with a point density gradient greater than a point density candidate threshold or a normal mutation intensity greater than a normal candidate threshold are marked as edge points; The edge points are subjected to spatial connectedness closure processing, and a continuous earthwork edge candidate region is output.

8. The method of claim 7, wherein, The surface roughness takes the average distance of the neighborhood points to the fitted plane, the particle size feature takes the equivalent diameter of the spatially connected aggregated neighborhood point cloud, the normal distribution takes the variance of the neighborhood normal angle, and the height variation trend takes the average height difference in the main direction of the neighborhood. 9.The method of claim 8, wherein, Determining the real boundary of the earthwork refers to: A profile sampling line is established in the normal direction of each edge point in the earthwork edge candidate region, and a sequence of intersection points of the double-time-phase three-dimensional point cloud is extracted on the profile sampling line to form a profile pair; The normal projection distance difference and its mean square value of the edge intersection points are calculated for each profile pair, the intersection points with a mean square value lower than a stability threshold are taken as stable edge positions, and the stable edge positions are subjected to median constraint along the boundary tangent; The stable edge positions are connected according to the spatial adjacency relationship and subjected to closure verification to obtain a closed real boundary of the earthwork, and the three-dimensional point cloud model is cropped by the real boundary of the earthwork to obtain a cropped point cloud.

10. The method of claim 9, wherein, The execution of the closed check refers to: connecting the stable edge positions in spatial order to form a boundary path; calculating the spatial distance between the first and last positions along the boundary path and the consistency of the overall path in the ring direction; when the spatial distance between the first and last positions is less than a preset distance threshold and the path ring direction remains consistent, the boundary path is determined as a closed boundary; when the path does not meet the closed condition, the corresponding boundary path is rejected for use in three-dimensional point cloud clipping and volume calculation.

Images