Water conservancy earthwork measurement method, device and equipment for radar point cloud real-time modeling and medium

Abstract

The invention relates to a radar point cloud real-time modeling water conservancy earthwork measurement method and device, equipment and a medium. According to the method, a vegetation root system interference area is scanned and recognized through a dual-frequency polarization radar, a coordinate set of the vegetation root system interference area is obtained, an elevation compensation model is constructed based on dielectric constant difference and fractal dimension to generate a compensation parameter set, and original point cloud is processed by combining a bare soil reference curved surface and adopting a dynamic constraint particle filtering algorithm to obtain real ground point cloud. Then, a corrected terrain curved surface is generated through parasitic energy field constraint optimization, and finally, compensation parameters are fused to realize high-precision earth volume calculation in a root system interference environment, so that point cloud distortion and a terrain curved surface systematic lifting effect caused by a vegetation root system are greatly reduced, and an earth measurement result of a complex vegetation coverage area meets an engineering acceptance standard.

Description

Technical Field

[0001] The present invention belongs to the technical field of earthwork measurement, and in particular relates to a water conservancy earthwork measurement method, device, equipment and medium for real-time modeling of radar point clouds. Background Art

[0002] In the fields of water conservancy project construction and slope management, accurate measurement of earthwork quantities is directly related to construction cost control and structural safety. While airborne LiDAR and photogrammetry, the current mainstream measurement methods, can quickly acquire three-dimensional surface point clouds, they struggle to penetrate dense vegetation cover, significantly obscuring the true slope topography. Particularly in areas with well-developed grass and shrub root systems, traditional point cloud filtering methods misidentify tangled roots as surface, resulting in false uplift in the reconstructed terrain surface and systematically inflating earthwork calculations. While recent research has attempted to enhance surface penetration using ground-based millimeter-wave radar, dielectric interference caused by vegetation roots distorts the echo signal, creating parasitic features in the point cloud data that resemble microscopic geological structures. This interference effect is nonlinearly amplified by changes in soil moisture, making it impossible for existing algorithms to distinguish between the true rock-soil interface and root artifacts.

[0003] In existing technologies, the engineering community usually uses manual sampling combined with empirical coefficient correction to compensate for vegetation interference. This method is not only time-consuming and labor-intensive, but the accuracy of compensation is also limited by the subjective judgment of technicians. Although the multi-band fusion algorithm proposed by some scholars can identify shallow roots, there is a significant deviation in the estimation of the penetration depth of deep roots. More importantly, there is a fundamental flaw in the current technology chain - the point cloud distortion caused by the root system has both signal delay offset caused by the sudden change in dielectric constant and multipath scattering effect caused by complex fractal structure. Both are often generally classified as "noise" in existing models and eliminated, but their systematic elevation effect on elevation is ignored. This technical blind spot directly leads to the actual filling volume often exceeding the design value during the acceptance of slope protection projects, triggering a large number of engineering claims disputes.

[0004] It's worth noting that errors in earthwork volume calculations are fundamentally due to systematic distortions in the terrain surface. Traditional solutions focus on either point cloud filtering (such as improved cloth filtering algorithms) or elevation correction in post-processing. However, no closed-loop compensation mechanism has been established, from dielectric property analysis to terrain reconstruction. This fragmented approach fails to fundamentally decouple the electromagnetic coupling between soil and roots, let alone quantitatively eliminate root interference.

[0005] Therefore, the water conservancy industry urgently needs a method that can solve root interference parameters in real time and dynamically embed compensation mechanisms into the entire terrain modeling process to end the industry dilemma of "the amount of earthwork with vegetation can never be accurately measured." Summary of the Invention

[0006] Based on this, it is necessary to provide a water conservancy earthwork measurement method, device, equipment and medium based on real-time modeling of radar point cloud to address the above technical problems.

[0007] In a first aspect, the present application provides a water conservancy earthwork measurement method based on real-time radar point cloud modeling, comprising:

[0008] S1. Scan the target slope vegetation area using a dual-frequency polarimetric radar to obtain raw point cloud data. Based on the dielectric response differences and fractal geometry characteristics of the point cloud clusters in the raw point cloud data, identify and mark the root interference area and obtain the root interference area coordinate set.

[0009] S2. Constructing an elevation compensation model and generating compensation parameters based on the root interference area coordinate set and the dielectric constant difference and fractal dimension of each root interference area;

[0010] S3. Based on the compensation parameters and the bare soil reference surface, the dynamic constrained particle filter algorithm is used to process the original point cloud data to obtain the ground point cloud data;

[0011] S4. Based on the ground point cloud data, a parasitic energy field constraint driven by compensation parameters is applied to the terrain surface initially generated in the root interference area. Through energy minimization optimization, a corrected terrain surface is output after eliminating the parasitic interference of the roots.

[0012] S5. Based on the corrected terrain surface and compensation parameters, the earthwork volume is calculated to obtain the earthwork engineering volume.

[0013] In a second aspect, the present application also provides a water conservancy earthwork measurement device for real-time modeling of radar point clouds, comprising:

[0014] The root interference area identification module is used to scan the target slope vegetation area using a dual-frequency polarization radar to obtain raw point cloud data. Based on the dielectric response differences and fractal geometry characteristics of the point cloud clusters in the raw point cloud data, the root interference area is identified and marked, and a set of root interference area coordinates is obtained.

[0015] The elevation compensation parameter generation module is used to construct an elevation compensation model and generate compensation parameters based on the coordinate set of the root interference area and the dielectric constant difference and fractal dimension of each root interference area;

[0016] The point cloud dynamic filtering module is used to process the original point cloud data using the dynamic constrained particle filtering algorithm based on the compensation parameters and the bare soil reference surface to obtain the ground point cloud data;

[0017] The terrain surface correction module is used to apply parasitic energy field constraints driven by compensation parameters to the terrain surface initially generated in the root interference area based on ground point cloud data. Through energy minimization optimization, it outputs the corrected terrain surface after eliminating the parasitic interference of the root system;

[0018] The earthwork volume real-time calculation module is used to calculate the earthwork volume based on the corrected terrain surface and compensation parameters to obtain the earthwork engineering volume.

[0019] In a third aspect, the present application also provides a computer device comprising a memory and a processor, wherein the memory stores a computer program, and when the processor executes the computer program, it implements a water conservancy earthwork measurement method based on real-time modeling of radar point clouds as in the first aspect.

[0020] In a fourth aspect, the present application further provides a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements a water conservancy earthwork measurement method based on real-time modeling of radar point clouds as in the first aspect.

[0021] The above-mentioned water conservancy earthwork measurement method, device, equipment and medium based on real-time modeling of radar point cloud identify the vegetation root interference area through dual-frequency polarization radar scanning and obtain its coordinate set, construct an elevation compensation model based on the dielectric constant difference and fractal dimension to generate a compensation parameter set, and use the dynamic constrained particle filter algorithm to process the original point cloud in combination with the bare soil reference surface to obtain the real ground point cloud, and then generate the corrected terrain surface through parasitic energy field constraint optimization, and finally integrate the compensation parameters to realize high-precision earthwork volume calculation in the root interference environment, thereby greatly reducing the point cloud distortion and the systematic uplift effect of the terrain surface caused by the vegetation root system, so that the earthwork measurement results in the complex vegetation coverage area meet the engineering acceptance standards. BRIEF DESCRIPTION OF THE DRAWINGS

[0022] In order to more clearly illustrate the technical solutions in the embodiments of the present application or related technologies, the following briefly introduces the drawings required for use in the embodiments or related technical descriptions. Obviously, the drawings described below are only some embodiments of the present application. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying any creative work.

[0023] Figure 1 A schematic flow chart of a water conservancy earthwork measurement method for real-time radar point cloud modeling provided by the present invention;

[0024] Figure 2 A schematic diagram of a process for generating ground point cloud data in an optional embodiment of the present invention;

[0025] Figure 3 This is a structural schematic diagram of a water conservancy earthwork measurement device based on real-time radar point cloud modeling provided by the present invention. DETAILED DESCRIPTION

[0026] In order to make the purpose, technical solutions and advantages of this application more clear, the following further describes this application in detail with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain this application and are not intended to limit this application.

[0027] refer to Figure 1 , which presents a flow chart of a water conservancy earthwork measurement method for real-time modeling of radar point clouds provided by this application, the method comprising the following steps:

[0028] S1. Scan the target slope vegetation area using a dual-frequency polarimetric radar to obtain raw point cloud data. Based on the dielectric response differences and fractal geometry characteristics of the point cloud clusters in the raw point cloud data, identify and mark the root interference areas and obtain the root interference area coordinate set.

[0029] Specifically, in water conservancy project construction and slope management scenarios, a dual-frequency polarization radar calibration device is used to scan the vegetation-covered area of ​​the target slope according to preset scanning trajectories and parameters. The two different frequency bands of electromagnetic waves emitted by the radar, after penetrating the vegetation canopy, undergo complex interactions with the vegetation roots and soil medium, including reflection, scattering, and refraction. For example, the low-frequency band is used to penetrate vegetation to detect deep soil structure, while the high-frequency band is used to obtain detailed surface and shallow features.

[0030] The acquisition of raw point cloud data involves time-domain sampling and spatial positioning of radar signals. After the radar receiving antenna captures the echo signal, the round-trip time is recorded by a clock synchronization system. Combined with the radar wave velocity (corrected to account for the dielectric constant of the medium), the spatial coordinates of the scattering point are calculated using geometric triangulation principles. Each point cloud data record also includes echo intensity information, which is related to factors such as the dielectric properties of the medium, surface roughness, and angle of incidence.

[0031] Based on the original point cloud data, the dielectric response difference characteristics and fractal geometry characteristics contained in the point cloud clusters are deeply explored. Different substances have different dielectric properties. There are significant differences in dielectric constants between vegetation roots and media such as bare soil, which leads to different time delays, attenuation and other responses when radar waves propagate inside them. By analyzing the echo signal characteristics of each point in the point cloud data, the dielectric response characteristic indicators related to the vegetation roots are extracted, which is used as one of the bases for identifying root interference areas. In terms of fractal geometry characteristics, the vegetation roots present a unique fractal structure in space, with specific parameter characteristics such as fractal dimension. With the help of fractal geometry theory and calculation methods, the spatial distribution morphology and complexity of point cloud clusters can be quantitatively analyzed, and the fractal dimension values ​​of each area can be calculated.

[0032] To identify root disturbance areas, a machine learning-based classification model can be constructed. First, multiple feature dimensions are extracted from the raw point cloud data: these include the amplitude decay rate of the echo signal (reflecting differences in dielectric response), the density gradient of the point cloud clusters (related to the fractal structure of the root system), and the spatial continuity of the point cloud (root systems typically exhibit a continuous fractal network structure). By training on a large number of labeled root disturbance samples and bare soil samples, the model learns the distribution patterns of different data categories in the feature space.

[0033] In the actual classification process, the raw point cloud data can be divided into several analysis windows. The above-mentioned feature vector is calculated for the point cloud within each window and input into the trained classification model for prediction. The model outputs the probability value of each analysis window belonging to the root disturbance area. When the probability value exceeds a set threshold (such as 0.8), the window is determined to be a root disturbance area and its boundary coordinates are recorded. In addition, to improve recognition accuracy, a spatial filtering algorithm can be introduced to perform spatial correlation analysis on the classification results of adjacent analysis windows, eliminating isolated misclassified areas and ultimately forming a coherent and accurate set of root disturbance area coordinates.

[0034] S2. Based on the coordinate set of the root interference area and the dielectric constant difference and fractal dimension of each root interference area, an elevation compensation model is constructed to generate compensation parameters.

[0035] Specifically, the core of building the elevation compensation model lies in establishing a quantitative relationship between the dielectric properties and geometric characteristics of the root disturbance area and the terrain elevation deviation. First, for each root disturbance area, the true dielectric constant is obtained using soil sample data collected on the spot. Based on the difference between the true dielectric constant and the equivalent dielectric constant measured by radar, the dielectric constant difference Δε is calculated. This difference is calculated by the formula Δε=ε measured -ε real Calculate. Among them, ε measured is the equivalent dielectric constant derived from radar inversion, ε real The actual soil dielectric constant measured in the laboratory.

[0036] The box dimension method can be used to calculate the fractal dimension D. On the point cloud data of the root disturbance area, a series of grid boxes with different side lengths are constructed, and the minimum number of boxes N required to cover the area is calculated. As the grid side length r decreases, the relationship between N and r follows the formula N~r -D , by linear fitting the logN-logr data points, we can get an estimate of the fractal dimension D.

[0037] The elevation compensation model can be constructed based on the energy conservation principle and electromagnetic wave propagation theory. The model assumes that the propagation delay Δt of the radar wave in the vegetation root medium is proportional to the dielectric constant difference Δε, and the path factor k related to the fractal dimension D is proportional to the propagation delay Δt of the radar wave in the vegetation root medium. D Related, that is: Δt=kε ·Δε+k D D. Where, k ε and k D is an empirical coefficient obtained through experimental calibration. The relationship between terrain elevation deviation ΔH and propagation delay Δt is: ΔH = 0.5·c·Δt, where c is the speed of the radar wave in a vacuum.

[0038] The compensation parameter generation process involves calculating the ΔH value for each root disturbance zone, combining it with its area weight (larger areas have a more significant impact on the overall terrain), and taking a weighted average of the ΔH values ​​for all root disturbance zones to obtain a global compensation parameter. Furthermore, to account for the dynamic influence of soil moisture on dielectric properties and root structure, the compensation model can also incorporate a moisture correction factor, derived by comparing and analyzing measurement data from different time periods, to further improve compensation accuracy.

[0039] S3. Based on the compensation parameters and the bare soil reference surface, the dynamic constrained particle filter algorithm is used to process the original point cloud data to obtain the ground point cloud data.

[0040] Specifically, during the initialization phase of the dynamic constrained particle filter algorithm, the initial state distribution of the particle filter is determined based on the root interference zone coordinate set and compensation parameters. Each particle in the particle swarm represents a possible ground point cloud classification hypothesis, containing a set of label vectors for the point cloud to be classified and a corresponding probability density function. The initial particle distribution can be based on prior knowledge of the root interference zone. Specifically, within the root interference zone, the probability of ground points is low, while the probability of non-ground points (such as vegetation and root artifacts) is high; the probability of ground points is high near the bare soil reference surface.

[0041] In the particle filter iteration process, a dynamic constraint mechanism is introduced. The constraints may include:

[0042] 1) Elevation constraint: Based on the compensation parameter, the elevation deviation range of the particle state is set. When the particle state deviates from the bare soil reference surface by more than the deviation allowed by the compensation parameter, it is given a lower weight.

[0043] 2) Spatial continuity constraint: The gradient information of the bare soil reference surface is used to constrain the movement direction of particles in space to avoid unreasonable classification results at places where the terrain changes suddenly.

[0044] 3) Fractal feature constraints: Near the root interference area, the spatial distribution of the particle state is constrained according to the fractal dimension characteristics, and the classification hypothesis that conforms to the root fractal structure is retained first.

[0045] The particle filter update process can be implemented using a sequential Monte Carlo method. At each iteration, a new particle swarm is generated through resampling based on the current particle state and weight. The probability of resampling is proportional to the particle weight, ensuring that high-weight particles (those that better align with the classification hypothesis of the actual terrain characteristics) account for a larger proportion of the new generation of particles.

[0046] To further improve filtering accuracy, the algorithm can also incorporate multi-resolution analysis techniques. In the initial iteration phase, the point cloud data is processed at a lower resolution to quickly remove obvious non-ground points. As the number of iterations increases, the resolution is gradually increased to finely classify detailed features. Ultimately, guided by global compensation parameters, after multiple rounds of iterative optimization, high-precision ground point cloud data is obtained.

[0047] S4. Based on the ground point cloud data, a parasitic energy field constraint driven by compensation parameters is applied to the terrain surface initially generated in the root interference area. Through energy minimization optimization, the corrected terrain surface after eliminating the parasitic interference of the roots is output.

[0048] Specifically, the optimization process of the corrected terrain surface can be based on the calculus of variations and partial differential equation theory. First, the initially generated terrain surface is represented as a discrete grid function z(x,y), where x and y are plane coordinates. The parasitic energy field U(x,y) is driven by the compensation parameters and can be defined as:

[0049] U(x,y)=α·(ΔH(x,y)) 2 +β·D(x,y);

[0050] Among them, α and β are weight coefficients, ΔH(x, y) is the elevation deviation of the point, and D(x, y) is the spatial distribution function of the fractal dimension.

[0051] The goal of energy minimization optimization can be to minimize the following energy functional:

[0052]

[0053] Among them, the first term is the smoothness constraint term of the terrain surface, which ensures that the optimized surface remains smooth in the non-interference area; μ is the parameter that controls the intensity of the parasitic energy field; the second term is the data fidelity term, which ensures that the optimized surface is consistent with the initially generated terrain surface z initial Maintain consistency overall; λ is a parameter that balances the weights of the two items.

[0054] Solving the Euler-Lagrange equation for this energy functional yields a second-order nonlinear partial differential equation. This partial differential equation is discretized using the finite difference method, constructing a system of linear equations at the grid points. This linear system can be solved using an iterative solver (such as the conjugate gradient method), gradually adjusting the elevation of the terrain surface to reduce the energy functional E[z].

[0055] During the optimization process, different iterative step-size strategies are employed for root-interference areas and non-interference areas. In root-interference areas, due to the higher parasitic energy field intensity, a smaller iterative step-size is used to fine-tune the elevation values; in non-interference areas, a larger iterative step-size is used for faster convergence. Furthermore, a multi-scale optimization approach can be employed, first optimizing on a coarse grid to obtain preliminary correction results, which are then interpolated onto a fine grid for further refinement, ultimately outputting a highly accurate corrected terrain surface.

[0056] S5. Based on the corrected terrain surface and compensation parameters, the earthwork volume is calculated to obtain the earthwork engineering volume.

[0057] Specifically, the earthwork volume calculation can adopt a hybrid calculation strategy that combines the volume integration method and the cross-section method. First, the corrected terrain surface is divided into multiple calculation sections. The section spacing is dynamically adjusted according to the complexity of the terrain changes, for example, the spacing is reduced in steep terrain areas. On each section, the terrain curve is fitted using the cubic spline interpolation method to calculate the volume between adjacent sections. The volume calculation formula is: V = ∫∫(z(x,y)-z design )dxdy. Among them, z(x,y) is the elevation of the corrected terrain surface, z design is the design elevation. The integration range covers the entire earthwork calculation area.

[0058] In order to improve the calculation accuracy, the compaction effect of soil and the influence of moisture content change on volume can be considered. The volume correction coefficient K is introduced, and its calculation formula is: K = 1 + γ (1-e -δ·H ). Where γ is the compaction coefficient related to soil type, δ is the moisture content sensitivity coefficient, and H is the average fill and cut height. The corrected earthwork volume is: V corrected =V·K.

[0059] Furthermore, the compensation parameters are combined to perform uncertainty analysis on the earthwork calculation results. Monte Carlo simulation can be used to propagate the uncertainty of the compensation parameters and evaluate the confidence interval of the earthwork calculation results, providing a reliable basis for engineering decision-making.

[0060] Throughout the calculation process, the algorithm divides the calculation area into grids, and the grid size can be dynamically adjusted based on the required terrain details and computational accuracy. When computing resources are limited, adaptive mesh refinement technology can be used to automatically refine the grid in areas with drastic terrain changes and appropriately relax the grid size in flat areas, ensuring a balance between computational efficiency and accuracy. Ultimately, the output earthwork volume data includes fill, excavation, and total earthwork volume, and displays the spatial distribution of earthwork in a 3D visualization, providing an intuitive quantitative basis for water conservancy project construction.

[0061] The above-mentioned water conservancy earthwork measurement method based on real-time modeling of radar point clouds identifies the vegetation root interference area through dual-frequency polarization radar scanning and obtains its coordinate set, constructs an elevation compensation model based on the dielectric constant difference and fractal dimension to generate a compensation parameter set, and uses a dynamic constrained particle filter algorithm to process the original point cloud in combination with the bare soil reference surface to obtain a real ground point cloud, and then generates a corrected terrain surface through parasitic energy field constraint optimization. Finally, the compensation parameters are integrated to achieve high-precision earthwork volume calculation in a root interference environment, thereby greatly reducing the point cloud distortion and the systematic uplift effect of the terrain surface caused by the vegetation root system, so that the earthwork measurement results in complex vegetation coverage areas meet the engineering acceptance standards.

[0062] In an optional embodiment, S1 includes the following steps:

[0063] S11. Acquire original point cloud data using a millimeter-wave radar operating at a first frequency and a second frequency; wherein the first frequency is lower than the second frequency.

[0064] Specifically, in the water conservancy earthwork measurement scenario, millimeter-wave radar equipment with first and second working frequencies is selected. The first frequency is relatively low, the wavelength of the electromagnetic wave is longer, and it has stronger vegetation penetration ability, which can effectively detect the deep soil structure under vegetation cover; the second frequency is relatively high, the wavelength is shorter, and it is mainly used to obtain fine feature information of the surface and shallow layer. The radar equipment scans the vegetation area of ​​the target slope according to the preset scanning path and parameters, emitting electromagnetic waves and receiving reflected echo signals. Through time domain sampling and processing of the echo signal, combined with the radar wave velocity (corrected wave velocity after considering the influence of the dielectric constant of the medium) and the geometric position and attitude information of the transmitting and receiving antennas, the principle of geometric triangulation is used to accurately calculate the spatial position coordinates of the scattering point, thereby obtaining the original point cloud data set containing spatial information and reflection characteristic information.

[0065] S12. For point cloud clusters in the original point cloud data whose diameters are smaller than the cluster size threshold, calculate the dual-frequency polarization intensity ratio. The calculation formula for the dual-frequency polarization intensity ratio is:

[0066]

[0067] Among them, R polar,i represents the dual-frequency polarization intensity ratio of the i-th point cloud cluster, I2 represents the second frequency, and I1 represents the first frequency.

[0068] Specifically, for the original point cloud data obtained, a cluster size threshold is first set to filter out point cloud clusters with diameters smaller than the threshold. These smaller point cloud clusters are more likely to be abnormal point cloud clusters affected by vegetation root interference. For each qualified point cloud cluster screened out, its dual-frequency polarization intensity ratio R is calculated. polar,i The specific calculation formula is: R polar,i =(I2-I1) / I1. Where I1 represents the echo intensity of the point cloud cluster at the first frequency, and I2 represents the echo intensity of the point cloud cluster at the second frequency. This intensity ratio can reflect the difference in echo characteristics of the point cloud cluster at different frequencies, providing a quantitative basis for the subsequent identification of potential interference areas. Generally speaking, interference sources such as vegetation roots will cause significant differences in echo intensity at different frequencies, making R polar,i The value is larger.

[0069] S13. When the dual-frequency polarization intensity ratio is greater than the dielectric threshold, the corresponding area is marked as a potential interference area, and the fractal dimension of the point cloud cluster in the potential interference area is calculated; the calculation formula of the fractal dimension is:

[0070]

[0071] Among them, FD i represents the fractal dimension of the i-th point cloud cluster, r i is the scale parameter corresponding to the i-th point cloud cluster, N(r i ) is expressed as i The number of point clouds contained in the lower grid cell.

[0072] Specifically, a dielectric threshold is set to determine whether the dual-frequency polarization intensity ratio of the point cloud cluster exceeds the normal range. polar,i When the value is greater than the dielectric threshold, it indicates that there is potential root interference in the area, and the corresponding area is marked as a potential interference area. For the point cloud clusters in these potential interference areas, the fractal dimension FD is further calculated. i , to quantitatively describe the complexity and fractal characteristics of its spatial distribution. The calculation formula of fractal dimension is: FD i =logN(r i ) / log(1 / r i ). Among them, r i is the scale parameter corresponding to the i-th point cloud cluster, which indicates the grid unit size used when analyzing the point cloud cluster; N(r i ) represents the scale r iThe number of point clouds contained in the grid cell below. Fractal dimension can reflect the spatial filling capacity and complexity of point cloud clusters at different scales. Vegetation roots usually have higher fractal dimension values ​​because they have a complex fractal structure in space.

[0073] S14. The potential interference areas with fractal dimensions greater than the fractal threshold are regarded as root interference areas, and a root interference area coordinate set of all root interference areas is output.

[0074] Specifically, a fractal threshold is set to further screen out the real root disturbance area. i Potential interference areas greater than this fractal threshold are identified as root interference areas. The point cloud clusters in these areas not only exhibit anomalies in the dual-frequency polarization intensity ratio but also possess high fractal dimensions consistent with the characteristics of vegetation roots, enabling accurate identification of areas affected by root interference. Finally, the coordinate information of all identified root interference areas is extracted and integrated, and a root interference area coordinate set is output. This coordinate set records in detail the spatial position of each root interference area within the target slope vegetation area, providing precise spatial positioning information for subsequent elevation compensation model construction and terrain modeling compensation, ensuring that subsequent processing steps can effectively compensate and correct these interference areas, thereby improving the accuracy of earthwork measurement.

[0075] In an optional embodiment, S2 includes the following steps:

[0076] S21. Select control points in the exposed soil area adjacent to the target slope vegetation area, calculate the median echo intensity of the control points, and calibrate the dielectric constant reference value of the soil in the target slope vegetation area based on the median echo intensity.

[0077] Specifically, multiple control points were selected within the exposed soil area adjacent to the target slope vegetation area to ensure that these control points were representative and could accurately reflect the dielectric properties of the soil in that area. Multiple radar scans were performed on each control point to obtain echo intensity data, and the median of this data was calculated to reduce the interference of random noise. Based on the calculated median echo intensity, the dielectric constant baseline value of the soil in that area was inferred using an electromagnetic wave inversion algorithm and a soil dielectric property model. This baseline value reflects the true dielectric properties of exposed soil in the absence of vegetation root interference, providing a reference standard for subsequent assessments of dielectric deviations in root interference areas.

[0078] S22, traverse each root interference cluster in the root interference area coordinate set, and inversely calculate the root dielectric constant based on the dual-frequency polarization intensity ratio of the root interference cluster; the calculation formula of the root dielectric constant is:

[0079] ε root,i =a×R polar,i +b;

[0080] Among them, ε root,i represents the root dielectric constant of the i-th root interference cluster, and a and b are the dielectric conversion coefficients.

[0081] Specifically, each root interference cluster in the root interference area coordinate set is traversed, and for each interference cluster, its dual-frequency polarization intensity ratio R polar,i , combined with the predetermined dielectric conversion coefficients a and b, the root dielectric constant ε is inversely calculated root,i , the calculation formula is: root,i =a×R polar,i +b. Here, a and b are dielectric conversion coefficients, obtained through experimental calibration or theoretical derivation, used to convert the dual-frequency polarization intensity ratio into the corresponding dielectric constant value. This step quantifies the root system dielectric constant, making the assessment of root disturbance more objective and operational.

[0082] S23. Calculate the dielectric difference corresponding to the root interference cluster based on the dielectric constant baseline value and the root system dielectric constant. The calculation formula for the dielectric difference is:

[0083] Δε i =ε root,i -ε soil ;

[0084] Among them, Δε i represents the dielectric difference corresponding to the i-th root interference cluster, ε soil is the reference value of dielectric constant.

[0085] Specifically, based on the dielectric constant reference value ε obtained in step S21 soil , calculate the dielectric difference Δε of each root interference cluster i The calculation formula is: Δε i =ε root,i -ε soil The dielectric difference reflects the difference in dielectric properties between the root disturbance area and the exposed soil area, and is a parameter for the subsequent construction of the elevation compensation model.

[0086] S24. Calculate the root depth corresponding to the root disturbance cluster based on the fractal dimension. The calculation formula for the root depth is:

[0087] δ root,i =c×(FD i )d;

[0088] Among them, δ root,i represents the root depth corresponding to the i-th root disturbance cluster, and c and d are the depth fitting coefficients.

[0089] Specifically, according to the fractal dimension FD of each root disturbance cluster i, combined with the depth fitting coefficients c and d, calculate the root depth δ root,i The calculation formula is: root,i =c×(FD i ) d Here, c and d are coefficients obtained by fitting experimental data and used to establish a quantitative relationship between fractal dimension and root depth. A higher fractal dimension indicates a more complex root structure, typically with a corresponding increase in root depth. This step converts the fractal dimension into root depth using a mathematical model, providing depth information for subsequent calculation of compensation parameters.

[0090] S25. Calculate the compensation coefficient of the corresponding root interference cluster based on the dielectric difference. The calculation formula of the compensation coefficient is:

[0091]

[0092] Among them, k i represents the compensation coefficient corresponding to the i-th root disturbance cluster, and h is the compensation fitting coefficient.

[0093] Specifically, the dielectric difference Δε is used i , combined with the compensation fitting coefficient h, calculate the compensation coefficient k of each root disturbance cluster i The calculation formula is: i =1-e -h×Δεi Where h is the compensation fitting coefficient determined through experiments or numerical simulations and is used to control the nonlinear relationship between compensation intensity and dielectric difference. The compensation coefficient reflects the degree of elevation compensation for the root disturbance area; the larger the dielectric difference, the larger the compensation coefficient.

[0094] S26. Calculate compensation parameters based on the compensation coefficient, dielectric difference, and root depth. The calculation formula for the compensation parameters is:

[0095] d comp,i =k i ×Δε i ×δ root,i ;

[0096] Among them, d comp,i represents the compensation parameter corresponding to the i-th root disturbance cluster.

[0097] Specifically, the comprehensive compensation coefficient k i , dielectric difference Δε i and root depth δ root,i , calculate the compensation parameter d comp,i The calculation formula is: comp,i =k i ×Δε i ×δ root,iThis compensation parameter comprehensively considers the influence of root disturbance intensity (dielectric difference), depth, and compensation coefficient, and can quantitatively characterize the impact of each root disturbance cluster on terrain elevation. By applying these compensation parameters to subsequent point cloud data processing and terrain modeling, root disturbance can be effectively compensated and corrected, improving the accuracy of earthwork measurements.

[0098] refer to Figure 2 In an optional embodiment, S3 includes the following steps:

[0099] S31. For the exposed soil area point cloud that does not belong to the root interference area coordinate set in the original point cloud data, a smooth and continuous reference elevation surface is generated using the Kriging spatial interpolation algorithm.

[0100] Specifically, in the original point cloud data, the point cloud data of the exposed soil area that does not belong to the root interference area coordinate set is screened out. These point cloud data are relatively less affected by vegetation root interference, have higher reliability, and can better reflect the elevation information of the real surface. The Kriging spatial interpolation algorithm is used to process these exposed soil area point clouds to generate a smooth and continuous reference elevation surface. The Kriging interpolation algorithm is an interpolation method based on spatial autocorrelation theory. It can comprehensively consider the spatial distribution characteristics and mutual relationships between data points. By constructing a variogram model and an estimation function, it can make the best unbiased estimate of the unknown elevation points, thereby generating a more accurate reference elevation surface, which provides a benchmark reference for subsequent point cloud data filtering and terrain modeling.

[0101] S32. Based on the cloth simulation filtering algorithm, set the dynamic quality for each cloth particle participating in the simulation, including:

[0102] If the particle p j The corresponding spatial position does not belong to any root interference cluster in the root interference area coordinate set, then set the particle p j The mass m j To preset standard quality;

[0103] If the particle p j The corresponding spatial position belongs to a root disturbance cluster, and the mass is increased based on the preset standard mass, as the particle p j The expression of mass increase processing is:

[0104] m j =m0+f×d comp,i ;

[0105] Where m0 is the preset standard mass and f is the mass adjustment coefficient.

[0106] Specifically, based on the cloth simulation filtering algorithm, dynamic mass setting is performed on each cloth particle involved in the simulation. j If its corresponding spatial position does not belong to any root interference cluster in the root interference area coordinate set, then its mass m is directly j Set to the preset standard mass m0, that is, m j =m0.

[0107] If the particle p j If the spatial position of belongs to a root disturbance cluster, the mass is increased based on the preset standard mass m0 to enhance the inertial effect of the particle during the simulation and better resist the influence of root disturbance. The expression of mass increase is: j =m0+f×d comp,i Among them, f is the quality adjustment coefficient, which is used to control the influence of compensation parameters on quality; d comp,i This is the compensation parameter corresponding to the root disturbance cluster, reflecting the intensity and characteristics of the root disturbance in that area. This dynamic quality setting allows the cloth particles to more realistically reflect the mechanical behavior of the terrain surface under root disturbance during the simulation, providing a more reasonable data foundation for subsequent filtering and terrain modeling.

[0108] S33. Calculate the additional force vector generated by the parasitic interference of the root structure on each particle belonging to the root interference cluster based on the reference elevation surface and the compensation parameter. The expression of the additional force vector is:

[0109]

[0110] Among them, F root,j represents the additional force vector, Represents particle p j The gradient vector of the reference elevation surface at the location. The negative sign indicates that the additional force points to the true shape of the surface. g is the force field adjustment coefficient.

[0111] Specifically, Represents particle p j The gradient vector of the reference elevation surface at the location reflects the degree and direction of the surface inclination at that location; the negative sign indicates that the additional force points to the true surface shape, that is, the direction of the force is opposite to the direction of terrain distortion caused by root interference, which is used to correct terrain distortion; g is the force field adjustment coefficient, which is used to control the intensity of the additional force; d comp,i is the compensation parameter corresponding to the root disturbance cluster, reflecting the degree of influence of the root disturbance on the particle. By calculating the additional force vector, we can quantitatively characterize the impact of the root disturbance on the particle, providing a basis for subsequent particle motion correction.

[0112] S34. In the motion simulation of the cloth particles, the additional force vector is introduced to correct the particle motion equation to obtain a corrected motion equation; the corrected motion equation is expressed as:

[0113]

[0114] Among them, r j Represents particle p j The position coordinate vector in three-dimensional space, G is the vertical downward gravity acceleration vector, T j For particle p j The elastic force vector between adjacent cloth particles, t is the time.

[0115] Specifically, G is the vertical downward gravitational acceleration vector, reflecting the gravitational force on the particle; T j For particle p j The elastic tension vector between adjacent cloth particles reflects the mutual constraint between them; t is the time variable. By introducing the additional force vector Froot,j into the equation of motion, particles are affected not only by gravity and elastic tension during the simulation, but also by additional forces related to root disturbance. This allows for a more realistic simulation of the mechanical behavior and morphological changes of the terrain surface under root disturbance, providing more accurate particle motion for subsequent filtering and terrain modeling.

[0116] S35. Perform iterative update of particle positions based on the corrected motion equation, and output ground point cloud data when a preset convergence condition is reached.

[0117] Specifically, based on the revised equation of motion, the iterative update of the particle position is performed. In each iteration, the position coordinates of the particle at the next moment are calculated based on the current particle velocity and acceleration, and the particle status information is updated. During the iterative update process, the changes in the particle position and the convergence state of the system are continuously monitored until the preset convergence conditions are reached. The convergence conditions include the threshold of the particle position change, the stability of the system energy, and the upper limit of the number of iterations. When the convergence conditions are met, the iterative update is stopped and the final ground point cloud data is obtained. After filtering and root interference correction, these ground point cloud data can more accurately reflect the morphological characteristics of the real surface, providing a high-quality data foundation for subsequent terrain modeling and earthwork calculation.

[0118] In an optional embodiment, S4 includes the following steps:

[0119] S41. Based on the ground point cloud data, a Poisson surface reconstruction algorithm is used to generate a preliminary terrain surface.

[0120] Specifically, the Poisson surface reconstruction algorithm is a 3D surface reconstruction method based on an indicator function. It treats point cloud data as sample points of the surface, constructs an implicit function by solving the Poisson equation, and then extracts the zero-isosurface as the reconstructed terrain surface. The specific steps are as follows:

[0121] 1) Construct an indicator function: Based on the ground point cloud data, construct an indicator function, which takes the value of 1 at the point cloud sample point and the value of 0 at other locations.

[0122] 2) Calculate the gradient field: The gradient field is obtained by calculating the indicator function. The gradient field reflects the direction and intensity information of the terrain surface.

[0123] 3) Solving the Poisson equation: The gradient field is used to solve the Poisson equation and obtain an implicit function that has a higher value near the sample points of the point cloud.

[0124] 4) Extracting isosurfaces: Extracting the zero isosurface from the implicit function is the preliminary generated terrain surface, which can better fit the ground point cloud data.

[0125] S42, within the root interference area coordinate set, for the root interference cluster C i The energy field function used to express the parasitic interference effect of the root system is constructed based on the surface area of ​​​​the root system. The expression of the energy field function is:

[0126]

[0127] Among them, S represents the surface to be optimized, dC i Indicates that in cluster C i In-region points, Represents the gradient vector of the preliminary terrain surface at the corresponding position, S ref Indicates the elevation value of the reference elevation surface at the corresponding position; i is the coupling coefficient, λ i =k i ×Δε i ×q, q is the energy field adjustment coefficient.

[0128] Specifically, the construction of the energy field function aims to quantify the degree of distortion of the terrain surface in the root disturbance area by introducing relevant parameters of root disturbance and provide an objective function for subsequent optimization. It represents the gradient vector of the preliminary terrain surface at the corresponding position, reflecting the terrain change trend of the preliminary surface; S ref Indicates the elevation value of the reference elevation surface at the corresponding position, providing the benchmark information of the real surface; FD i Root disturbance cluster C i The fractal dimension reflects the complexity of the root system structure; iis the coupling coefficient, and the calculation formula is: i =k i ×Δε i ×q. Among them, k i is the compensation coefficient, reflecting the intensity of root disturbance; Δε i The dielectric difference represents the difference in dielectric properties between the root system and the soil. q is the energy field adjustment coefficient, which is used to balance the weights of each term in the energy field function. This energy field function comprehensively considers the gradient difference between the terrain surface and the preliminary surface, the elevation difference from the reference elevation surface, and the characteristics of root interference. By minimizing this energy field function, the terrain surface can be effectively corrected in the root interference area.

[0129] S43. Using the reference elevation surface as the initial surface for iterative optimization, the conjugate gradient method is used to minimize the energy field function within the coordinate set of the root interference area.

[0130] Specifically, the conjugate gradient method is an iterative optimization algorithm suitable for large-scale unconstrained optimization problems. Its basic idea is to search for the step size that minimizes the objective function along the conjugate direction in each iteration, gradually approaching the optimal solution. During the optimization process, the gradient of the energy field function with respect to the surface S is calculated, and the surface elevation value is updated, so that the surface gradually converges to the actual surface shape in the root interference area. The specific steps are as follows:

[0131] 1) Initialization: Use the reference elevation surface as the initial surface S (0) , set the initial iteration step to 0, the maximum number of iterations and the convergence threshold and other parameters.

[0132] 2) Gradient calculation: Calculate the current surface S (n) The gradient of the energy field function E(S) at

[0133] 3) Determine the search direction: Determine the search direction for this iteration according to the formula of the conjugate gradient method.

[0134] 4) Step search: along the search direction, determine the step size that minimizes the energy field function through line search.

[0135] 5) Surface update: Update the surface S according to the step size and search direction (n) Get the new surface S (n+1) .

[0136] 6) Convergence judgment: Check whether the energy field function value or the surface update amplitude meets the convergence conditions. If so, stop the iteration; otherwise, return to step 2 to continue optimization.

[0137] Through the optimization of the conjugate gradient method, the distortion of the terrain surface in the root interference area is effectively corrected, and the terrain surface is closer to the actual surface morphology.

[0138] S44. When the iterative optimization process converges, the corrected terrain surface is output.

[0139] Specifically, after the iterative optimization process converges, a corrected terrain surface is output. This corrected terrain surface now takes into account the effects of root interference. By minimizing the energy field function, the effects of parasitic root interference on the terrain surface are eliminated, ensuring that the elevation and shape of the terrain surface in the root interference area more accurately reflect the actual surface conditions. This surface is used in subsequent earthwork calculations, providing a reliable basis for construction cost control and structural safety of water conservancy projects.

[0140] In an optional embodiment, S5 includes the following steps:

[0141] S51. Discretize the corrected terrain surface into a triangulated network composed of triangular facets; wherein the side length of each triangular facet is less than or equal to a preset value.

[0142] Specifically, the discretization process can use the Delaunay triangulation algorithm so that the side length of each triangular facet is less than or equal to a preset value. Delaunay triangulation can generate a triangular mesh that is relatively uniform in spatial distribution and has good geometric properties, avoiding triangles that are too long or too narrow, thereby ensuring the accuracy and stability of subsequent calculations. The preset side length threshold can be determined according to the complexity of the terrain and the required calculation accuracy, and can be set between a few centimeters and a few meters. In this way, the terrain surface is converted into a collection of a series of triangular facets, and the vertex coordinates and elevation information of each facet are derived from the corrected terrain surface.

[0143] S52, calculating the fill-cut height difference of each vertex in the triangulated network relative to the design elevation surface, and compensating for the fill-cut height difference based on whether the vertex is located within the root interference zone coordinate set, including:

[0144] When the vertex v k When the spatial position of is not located in the root interference cluster in the root interference area coordinate set, the vertex v k The fill and cut height is in, is the vertex v k The elevation value on the corrected terrain surface, H0 is the design elevation value;

[0145] When the vertex v k The spatial location belongs to the root disturbance cluster C i When the root system interference cluster C i The corresponding compensation parameter d comp,i , calculate vertex v k Fill height Among them, γ is the elevation compensation coefficient.

[0146] Specifically, when vertex v k When the spatial position of does not belong to any root disturbance cluster in the root disturbance area coordinate set, the calculation formula for the fill-cut height difference is: in, For vertex v k The elevation value on the corrected terrain surface, where H0 is the design elevation. This formula directly reflects the difference between the actual elevation at the vertex and the design elevation. A positive value indicates that fill is required, while a negative value indicates that cut is required.

[0147] When the vertex v k The spatial location belongs to the root disturbance cluster C i When considering the influence of root disturbance, the fill-cut height difference is compensated. The calculation formula of the compensated fill-cut height difference is:

[0148] Among them, d comp,i Root disturbance cluster C i The corresponding compensation parameter reflects the extent to which root disturbance affects elevation in that area; γ is the elevation compensation coefficient, used to adjust the compensation intensity. This formula corrects for elevation deviations caused by root disturbance by subtracting the compensation term, ensuring that the cut-and-fill height difference more accurately reflects actual construction requirements.

[0149] S53. Based on the area of ​​each triangular face and the fill-cut heights corresponding to the three vertices, the earthwork contribution corresponding to the triangular face is calculated using the truncated prism volume integral formula; the calculation formula for earthwork contribution is:

[0150]

[0151] Among them, V triangle is the earthwork contribution, S triangle is the area of ​​the triangular patch, Δh A , Δh B , Δh C They are the fill and cut heights corresponding to the three vertices of the triangular patch.

[0152] Specifically, S triangle is the area of ​​the triangular patch, which can be calculated from the vertex coordinates; Δh A , Δh B , Δh CThe height differences between the cut and fill heights corresponding to the three vertices of a triangular patch are calculated. This formula assumes that the elevation changes within the triangular patch are linearly distributed. By calculating the average of the cut and fill heights at the three vertices and multiplying it by the area of ​​the triangle, the earthwork contribution of the patch is obtained. This method effectively handles local elevation changes on curved terrain while maintaining accuracy, making it suitable for earthwork calculations on complex terrain.

[0153] S54. Accumulate the earthwork volume contributions of all triangular facets to obtain the total earthwork volume of the target slope vegetation area.

[0154] Specifically, the total earthwork volume is initialized first, and the initial value of the total earthwork volume Vtotal is set to 0. Then, the triangle face is traversed, and for each triangle face, the corresponding earthwork volume contribution V triangle Add to total earthwork volume V total After completing the traversal of all triangles, output the total earthwork volume V total , which is the total earthwork volume of the target slope vegetation area. By accumulating the earthwork contributions of all triangular patches, the earthwork volume of the entire area can be comprehensively and accurately calculated, providing reliable quantitative data support for construction planning and cost control of water conservancy projects.

[0155] The above-mentioned water conservancy earthwork measurement method based on real-time modeling of radar point clouds identifies the vegetation root interference area through dual-frequency polarization radar scanning and obtains its coordinate set, constructs an elevation compensation model based on the dielectric constant difference and fractal dimension to generate a compensation parameter set, and uses a dynamic constrained particle filter algorithm to process the original point cloud in combination with the bare soil reference surface to obtain a real ground point cloud, and then generates a corrected terrain surface through parasitic energy field constraint optimization. Finally, the compensation parameters are integrated to achieve high-precision earthwork volume calculation in a root interference environment, thereby greatly reducing the point cloud distortion and the systematic uplift effect of the terrain surface caused by the vegetation root system, so that the earthwork measurement results in complex vegetation coverage areas meet the engineering acceptance standards.

[0156] It should be understood that, although the various steps in the flowcharts involved in the various embodiments described above are displayed in sequence according to the instructions of the arrows, these steps are not necessarily executed in sequence in the order indicated by the arrows. Unless otherwise specified herein, there is no strict order restriction on the execution of these steps, and these steps can be executed in other orders. Moreover, at least a portion of the steps in the flowcharts involved in the various embodiments described above can include multiple steps or multiple stages, and these steps or stages are not necessarily executed and completed at the same time, but can be executed at different times, and the execution order of these steps or stages is not necessarily to be carried out in sequence, but can be executed in turn or alternately with other steps or at least a portion of steps or stages in other steps.

[0157] Based on the same inventive concept, embodiments of the present application also provide a device for implementing the aforementioned method for water conservancy earthwork measurement using real-time radar point cloud modeling. The solution provided by this device is similar to the solution described in the aforementioned method. Therefore, the specific limitations of the embodiments of one or more water conservancy earthwork measurement devices using real-time radar point cloud modeling provided below can be found in the limitations of the water conservancy earthwork measurement method using real-time radar point cloud modeling described above, and will not be further elaborated here.

[0158] In an exemplary embodiment, Figure 3 As shown, a water conservancy earthwork measurement device 30 for real-time modeling of radar point clouds is provided, comprising:

[0159] The root interference area identification module 31 is used to scan the target slope vegetation area through dual-frequency polarization radar to obtain raw point cloud data; based on the dielectric response differences and fractal geometric characteristics of the point cloud clusters in the raw point cloud data, the root interference area is identified and marked to obtain a set of root interference area coordinates.

[0160] The elevation compensation parameter generation module 32 is used to construct an elevation compensation model and generate compensation parameters based on the root interference area coordinate set and the dielectric constant difference and fractal dimension of each root interference area.

[0161] The point cloud dynamic filtering module 33 is used to process the original point cloud data using a dynamic constrained particle filtering algorithm based on the compensation parameters and the bare soil reference surface to obtain ground point cloud data.

[0162] The terrain surface correction module 34 is used to apply parasitic energy field constraints driven by compensation parameters to the terrain surface initially generated in the root interference area based on the ground point cloud data, and output the corrected terrain surface after eliminating the root parasitic interference through energy minimization optimization.

[0163] The earthwork volume real-time calculation module 35 is used to calculate the earthwork volume based on the corrected terrain surface and compensation parameters to obtain the earthwork engineering volume.

[0164] Optionally, the root disturbance zone identification module includes:

[0165] A dual-frequency point cloud data acquisition unit is used to obtain raw point cloud data through a millimeter wave radar having a first operating frequency and a second operating frequency; wherein the first frequency is lower than the second frequency.

[0166] The dual-frequency polarization intensity ratio calculation unit is used to calculate the dual-frequency polarization intensity ratio for point cloud clusters whose diameter is smaller than the cluster size threshold in the original point cloud data. The calculation formula of the dual-frequency polarization intensity ratio is:

[0167]

[0168] Among them, Rpolar,i represents the dual-frequency polarization intensity ratio of the i-th point cloud cluster, I2 represents the second frequency, and I1 represents the first frequency.

[0169] The potential interference area fractal feature analysis unit is used to mark the corresponding area as a potential interference area when the dual-frequency polarization intensity ratio is greater than the dielectric threshold, and calculate the fractal dimension of the point cloud cluster in the potential interference area; the calculation formula of the fractal dimension is:

[0170]

[0171] Among them, FD i represents the fractal dimension of the i-th point cloud cluster, r i is the scale parameter corresponding to the i-th point cloud cluster, N(r i ) is expressed as i The number of point clouds contained in the lower grid cell.

[0172] The root interference area determination and output unit is used to take the potential interference area with a fractal dimension greater than a fractal threshold as the root interference area, and output the root interference area coordinate set of all the root interference areas.

[0173] Optionally, the elevation compensation parameter generation module includes:

[0174] The soil dielectric reference calibration unit is used to select control points in the exposed soil area adjacent to the target slope vegetation area, calculate the median echo intensity of the control points, and calibrate the dielectric constant reference value of the soil in the target slope vegetation area based on the median echo intensity.

[0175] The root dielectric constant inversion unit is used to traverse each root interference cluster in the root interference area coordinate set and inversely calculate the root dielectric constant based on the dual-frequency polarization intensity ratio of the root interference cluster. The calculation formula of the root dielectric constant is:

[0176] ε root,i =a×R polar,i +b;

[0177] Among them, ε root,i represents the root dielectric constant of the i-th root interference cluster, and a and b are the dielectric conversion coefficients.

[0178] The dielectric difference calculation unit is used to calculate the dielectric difference corresponding to the root interference cluster based on the dielectric constant baseline value and the root system dielectric constant; the calculation formula of the dielectric difference is:

[0179] Δε i =v root,i -ε soil ;

[0180] Among them, Δε irepresents the dielectric difference corresponding to the i-th root interference cluster, ε soil is the reference value of dielectric constant;

[0181] The root depth calculation unit is used to calculate the root depth corresponding to the root disturbance cluster according to the fractal dimension; the calculation formula of the root depth is:

[0182] δ root,i =c×(FD i ) d ;

[0183] Among them, δ root,i represents the root depth corresponding to the i-th root disturbance cluster, and c and d are the depth fitting coefficients.

[0184] The compensation coefficient calculation unit is used to calculate the compensation coefficient of the corresponding root interference cluster according to the dielectric difference. The calculation formula of the compensation coefficient is:

[0185]

[0186] Among them, k i represents the compensation coefficient corresponding to the i-th root disturbance cluster, and h is the compensation fitting coefficient.

[0187] The compensation parameter generation unit is used to calculate the compensation parameter according to the compensation coefficient, dielectric difference and root depth. The calculation formula of the compensation parameter is:

[0188] d comp,i =k i ×Δε i ×δ root,i ;

[0189] Among them, d comp,i represents the compensation parameter corresponding to the i-th root disturbance cluster.

[0190] Optional point cloud dynamic filtering module includes:

[0191] The reference elevation surface generation unit is used to generate a smooth and continuous reference elevation surface using the Kriging spatial interpolation algorithm for the exposed soil area point cloud that does not belong to the root interference area coordinate set in the original point cloud data.

[0192] The particle dynamic quality setting unit is used to set the dynamic quality of each cloth particle participating in the simulation based on the cloth simulation filtering algorithm, including:

[0193] If the particle p j The corresponding spatial position does not belong to any root interference cluster in the root interference area coordinate set, then set the particle p j The mass m j To preset standard quality;

[0194] If the particle p j The corresponding spatial position belongs to a root disturbance cluster, and the mass is increased based on the preset standard mass, as the particle p j The expression of mass increase processing is:

[0195] m j =m0+f×d comp,i ;

[0196] Where m0 is the preset standard mass and f is the mass adjustment coefficient.

[0197] The additional force vector calculation unit is used to calculate the additional force vector generated by the parasitic interference of the root structure on each particle belonging to the root interference cluster based on the reference elevation surface and compensation parameters. The expression of the additional force vector is:

[0198]

[0199] Among them, F root,j represents the additional force vector, Represents particle p j The gradient vector of the reference elevation surface at the location. The negative sign indicates that the additional force points to the true shape of the surface. g is the force field adjustment coefficient.

[0200] The motion equation correction unit is used to introduce additional force vectors into the motion simulation of cloth particles to correct the particle motion equation and obtain the corrected motion equation. The corrected motion equation is expressed as:

[0201]

[0202] Among them, r j Represents particle p j The position coordinate vector in three-dimensional space, G is the vertical downward gravity acceleration vector, T j For particle p j The elastic force vector between adjacent cloth particles, t is the time.

[0203] The particle position iteration output unit is used to perform iterative updates of particle positions based on the modified motion equation and output ground point cloud data when the preset convergence conditions are reached.

[0204] Optional terrain surface correction module includes:

[0205] The preliminary terrain surface generation unit is used to generate a preliminary terrain surface based on ground point cloud data using a Poisson surface reconstruction algorithm.

[0206] The energy field function construction unit is used to generate the energy field function for the root disturbance cluster C in the root disturbance zone coordinate set. iThe energy field function used to express the parasitic interference effect of the root system is constructed based on the surface area of ​​​​the root system. The expression of the energy field function is:

[0207]

[0208] Among them, S represents the surface to be optimized, dC i Indicates that in cluster C i In-region points, Represents the gradient vector of the preliminary terrain surface at the corresponding position, S ref Indicates the elevation value of the reference elevation surface at the corresponding position; i is the coupling coefficient, λ i =k i ×Δε i ×q, q is the energy field adjustment coefficient.

[0209] The surface optimization unit is used to use the reference elevation surface as the initial surface for iterative optimization and to minimize the energy field function within the coordinate set of the root interference area using the conjugate gradient method.

[0210] The corrected terrain surface output unit is used to output the corrected terrain surface after the iterative optimization process converges.

[0211] Optional real-time earthwork calculation module includes:

[0212] The triangulated mesh generation unit is used to discretize the corrected terrain surface into a triangulated mesh composed of triangular facets; wherein the side length of each triangular facet is less than or equal to a preset value.

[0213] The cut-fill height compensation calculation unit is used to calculate the cut-fill height difference of each vertex in the triangulation network relative to the design elevation surface, and compensate for the cut-fill height difference based on whether the vertex is located in the root interference zone coordinate set, including:

[0214] When the vertex v k When the spatial position of is not located in the root interference cluster in the root interference area coordinate set, the vertex v k The fill and cut height is in, is the vertex v k The elevation value on the corrected terrain surface, H0 is the design elevation value;

[0215] When the vertex v k The spatial location belongs to the root disturbance cluster C i When the root system interference cluster C i The corresponding compensation parameter d comp,i , calculate vertex v k Fill height Among them, γ is the elevation compensation coefficient.

[0216] The total earthwork volume calculation unit is used to calculate the earthwork volume contribution corresponding to each triangular facet based on the area of ​​each triangular facet and the fill and cut heights corresponding to the three vertices using the truncated prism volume integral formula; the calculation formula for earthwork volume contribution is:

[0217]

[0218] Among them, V triangle is the earthwork contribution, S triangle is the area of ​​the triangular patch, Δh A , Δh B , Δh C They are the fill and cut heights corresponding to the three vertices of the triangular patch.

[0219] S54. Accumulate the earthwork volume contributions of all triangular facets to obtain the total earthwork volume of the target slope vegetation area.

[0220] An embodiment of the present application further provides a computer device, including a memory and a processor, wherein the memory stores a computer program, and the processor implements the steps in the above-mentioned method embodiments when executing the computer program.

[0221] An embodiment of the present application further provides a computer-readable storage medium on which a computer program is stored. When the computer program is executed by a processor, the steps in the above-mentioned method embodiments are implemented.

[0222] For the device embodiments, since they basically correspond to the method embodiments, the relevant parts can be referred to the partial description of the method embodiments. The device embodiments described above are merely illustrative, wherein the components described as separate parts may or may not be physically separated, and the parts displayed as units may or may not be physical units, that is, they may be located in one place, or they may be distributed on multiple network units. Some or all of the modules can be selected according to actual needs to achieve the purpose of the disclosed solution. A person of ordinary skill in the art can understand and implement it without expending creative work.

[0223] The above-described embodiments merely represent several implementation methods of the embodiments of the present application. While the descriptions are relatively specific and detailed, they should not be construed as limiting the scope of the patent application. It should be noted that a person skilled in the art may make various modifications and improvements without departing from the concept of the embodiments of the present application, and these modifications and improvements fall within the scope of protection of the embodiments of the present application.

Claims

1. A water conservancy earthwork measurement method based on real-time modeling of radar point clouds, characterized in that: The method comprises: S1. Scanning the vegetation area on the target slope using a dual-frequency polarization radar to obtain raw point cloud data; identifying and marking root interference areas based on the dielectric response differences and fractal geometry characteristics of point cloud clusters in the raw point cloud data to obtain a set of root interference area coordinates; S2. constructing an elevation compensation model and generating compensation parameters based on the root interference area coordinate set and the dielectric constant difference and fractal dimension of each root interference area; S3. Based on the compensation parameters and the bare soil reference surface, a dynamic constrained particle filter algorithm is used to process the original point cloud data to obtain ground point cloud data; S4. Based on the ground point cloud data, applying parasitic energy field constraints driven by the compensation parameters to the terrain surface initially generated in the root interference area, and outputting a corrected terrain surface after eliminating the root parasitic interference through energy minimization optimization; S5. Calculate the earthwork volume based on the corrected terrain surface and the compensation parameters to obtain the earthwork engineering volume.

2. The method according to claim 1, characterized in that Said S1 comprises: S11. Acquire the raw point cloud data using a millimeter-wave radar operating at a first frequency and a second frequency; wherein the first frequency is lower than the second frequency; S12. Calculate the dual-frequency polarization intensity ratio for the point cloud clusters in the original point cloud data whose diameters are smaller than the cluster size threshold. The calculation formula for the dual-frequency polarization intensity ratio is: Among them, R polar,i represents the dual-frequency polarization intensity ratio of the i-th point cloud cluster, I2 represents the second frequency, and I1 represents the first frequency; S13. When the dual-frequency polarization intensity ratio is greater than the dielectric threshold, mark the corresponding area as a potential interference area, and calculate the fractal dimension of the point cloud cluster in the potential interference area; the calculation formula of the fractal dimension is: Among them, FD i represents the fractal dimension of the i-th point cloud cluster, r i is the scale parameter corresponding to the i-th point cloud cluster, N(r i ) is expressed as i The number of point clouds contained in the lower grid cell; S14. The potential interference areas with fractal dimensions greater than a fractal threshold are taken as the root interference areas, and a root interference area coordinate set of all the root interference areas is output.

3. The method according to claim 2, characterized in that The S2 includes: S21, selecting control points in an exposed soil area adjacent to the target slope vegetation area, calculating median echo intensities of the control points, and calibrating a reference dielectric constant of the soil in the target slope vegetation area based on the median echo intensities; S22, traversing each root interference cluster in the root interference area coordinate set, and inversely calculating the root dielectric constant according to the dual-frequency polarization intensity ratio of the root interference cluster; the calculation formula of the root dielectric constant is: e root,i =a×R polar,i +b; Among them, ε root,i represents the root dielectric constant of the i-th root disturbance cluster, a and b are dielectric conversion coefficients; S23. Calculate the dielectric difference value corresponding to the root interference cluster according to the dielectric constant reference value and the root system dielectric constant; the calculation formula of the dielectric difference value is: No i =e root,i -e soil ; Among them, Δε i represents the dielectric difference corresponding to the i-th root interference cluster, ε soil is the dielectric constant reference value; S24. Calculate the root depth corresponding to the root disturbance cluster according to the fractal dimension; the calculation formula for the root depth is: d root,i =c×(FD i ) d ; Among them, δ root,i represents the root depth corresponding to the i-th root disturbance cluster, c and d are the depth fitting coefficients; S25. Calculate the compensation coefficient of the corresponding root interference cluster according to the dielectric difference; the calculation formula of the compensation coefficient is: Among them, k i represents the compensation coefficient corresponding to the i-th root disturbance cluster, and h is the compensation fitting coefficient; S26. Calculate the compensation parameter according to the compensation coefficient, the dielectric difference, and the root depth. The calculation formula of the compensation parameter is: d comp,i =k i ×Δε i ×δ root,i ; Among them, d comp,i represents the compensation parameter corresponding to the i-th root disturbance cluster.

4. The method according to claim 3, characterized in that The S3 includes: S31, generating a smooth and continuous reference elevation surface using a Kriging spatial interpolation algorithm for the exposed soil area point cloud that does not belong to the root interference area coordinate set in the original point cloud data; S32. Based on the cloth simulation filtering algorithm, set the dynamic quality for each cloth particle participating in the simulation, including: If the particle p j The corresponding spatial position does not belong to any root interference cluster in the root interference area coordinate set, then set the particle p j The mass m j To preset standard quality; If the particle p j The corresponding spatial position belongs to a certain root interference cluster, and the mass is increased based on the preset standard mass, as the particle p j The quality of the mass increase process is expressed as follows: m j =m0+f×d comp,i ; Wherein, m0 is the preset standard mass, and f is the mass adjustment coefficient; S33. Calculate, based on the reference elevation surface and the compensation parameter, an additional force vector generated by parasitic interference of the root structure on each particle belonging to the root interference cluster; the expression of the additional force vector is: Among them, F root,j represents the additional force vector, Represents particle p j The gradient vector of the reference elevation surface at the location, where the negative sign indicates that the additional force is pointing towards the true surface form, and g is the force field adjustment coefficient; S34. In the motion simulation of the cloth particles, the additional force vector is introduced to modify the particle motion equation to obtain a modified motion equation; the modified motion equation is expressed as: Among them, r j Represents particle p j The position coordinate vector in three-dimensional space, G is the vertical downward gravity acceleration vector, T j For particle p j The elastic tension vector between adjacent cloth particles, t is the time; S35. Perform iterative update of particle positions based on the corrected motion equation, and output the ground point cloud data when a preset convergence condition is reached.

5. The method according to claim 4, characterized in that The S4 includes: S41, generating a preliminary terrain surface using a Poisson surface reconstruction algorithm based on the ground point cloud data; S42: within the root interference area coordinate set, for the root interference cluster C i The energy field function used to express the parasitic interference effect of the root system is constructed based on the surface area of ​​​​the root system. The expression of the energy field function is: Among them, S represents the surface to be optimized, dC i Indicates that in cluster C i In-region points, Represents the gradient vector of the preliminary terrain surface at the corresponding position, S ref represents the elevation value of the reference elevation surface at the corresponding position; i is the coupling coefficient, λ i =k i ×Δε i ×q, q is the energy field adjustment coefficient; S43, using the reference elevation surface as the initial surface for iterative optimization, and minimizing the energy field function within the coordinate set of the root interference area using a conjugate gradient method; S44. After the iterative optimization process converges, the corrected terrain surface is output.

6. The method according to any one of claims 1 to 5, characterized in that The S5 includes: S51, discretizing the corrected terrain surface into a triangulated mesh composed of triangular facets; wherein the side length of each triangular facet is less than or equal to a preset value; S52, calculating the fill-cut height difference of each vertex in the triangulated network relative to the design elevation surface, and compensating the fill-cut height difference according to whether the vertex is located in the root interference area coordinate set, including: When the vertex v k When the spatial position of is not located in the root disturbance cluster in the root disturbance area coordinate set, the vertex v k The fill and cut height is in, is the vertex v k The elevation value on the corrected terrain surface, H0 is the design elevation value; When the vertex v k The spatial location belongs to the root disturbance cluster C i When the root system interference cluster C i The corresponding compensation parameter d comp,i , calculate vertex v k Fill height Among them, γ is the elevation compensation coefficient; S53. Based on the area of ​​each triangular facet and the fill-cut heights corresponding to the three vertices, calculate the earthwork contribution corresponding to the triangular facet using the truncated prism volume integral formula; the calculation formula for the earthwork contribution is: Among them, V triangle Contribution to the earthwork volume, S triangle is the area of ​​the triangular patch, Δh A , Δh B , Δh C are the fill-cut heights corresponding to the three vertices of the triangular patch respectively; S54: Accumulate the earthwork volume contributions of all the triangular facets to obtain the total earthwork volume of the target slope vegetation area.

7. A water conservancy earthwork measurement device based on real-time modeling of radar point clouds, characterized in that: The device comprises: A root interference area identification module is used to scan the target slope vegetation area using a dual-frequency polarization radar to obtain raw point cloud data; based on the dielectric response differences and fractal geometric characteristics of the point cloud clusters in the raw point cloud data, the root interference area is identified and marked to obtain a set of root interference area coordinates; An elevation compensation parameter generation module is used to construct an elevation compensation model and generate compensation parameters based on the root interference area coordinate set and the dielectric constant difference and fractal dimension of each root interference area; A point cloud dynamic filtering module is used to process the original point cloud data using a dynamic constrained particle filtering algorithm based on the compensation parameters and a bare soil reference surface to obtain ground point cloud data; a terrain surface correction module for applying parasitic energy field constraints driven by the compensation parameters to the terrain surface initially generated in the root interference area based on the ground point cloud data, and outputting a corrected terrain surface after eliminating the root parasitic interference through energy minimization optimization; The earthwork volume real-time calculation module is used to calculate the earthwork volume based on the corrected terrain surface and the compensation parameters to obtain the earthwork engineering volume.

8. A computer device comprising a memory and a processor, wherein the memory stores a computer program, wherein: When the processor executes the computer program, the method according to any one of claims 1 to 6 is implemented.

9. A computer-readable storage medium having a computer program stored thereon, characterized in that: When the computer program is executed by a processor, the method according to any one of claims 1 to 6 is implemented.

Images