Determination method for soil erosion of waste slag field based on three-dimensional model

Abstract

The invention provides a waste slag field soil erosion determination method based on a three-dimensional model, and relates to the technical field of engineering construction, and the method comprises the following steps: obtaining an initial terrain image; constructing a three-dimensional vector model of the waste slag field; a plurality of image control points are arranged on the periphery of the waste slag field; installing a plurality of fixed oblique photogrammetry devices with three-axis attitude control; the method comprises the following steps: regularly collecting multi-angle image data of a waste slag field by using an oblique photogrammetry device, obtaining image data of different visual angles, and carrying out encryption matching on the collected multi-angle image data until a preset coordinate resolution requirement is met; processing the matched image data to generate three-dimensional point cloud data, and utilizing the three-dimensional point cloud data to construct a three-dimensional irregular triangulation network model of the waste slag field; calculating elevation change in a preset time interval; adopting a numerical analysis method to extract change data of the waste slag body; and evaluating the stability of the waste slag field according to the change data. According to the invention, high-precision, omnibearing and dynamic monitoring of the waste slag field can be realized.

Description

Technical Field

[0001] The invention relates to the technical field of engineering construction, and in particular to a method for measuring soil erosion in a waste dump based on a three-dimensional model. Background Art

[0002] In the process of engineering construction, a large amount of earthwork excavation, concrete demolition or construction waste will inevitably be generated. These wastes are usually concentrated in a specific area to form the so-called waste dump. With the rapid development of social economy and the advancement of large-scale infrastructure construction, the number and scale of waste dumps are also growing rapidly. Waste dumps are not only related to construction safety and solid waste pollution prevention and control, but also an important source of soil erosion. Therefore, effective soil erosion monitoring and safety monitoring of waste dumps have become the focus of engineering monitoring.

[0003] At present, the monitoring of waste dumps usually includes multiple aspects such as the amount of stockpiling, the amount of erosion, the external displacement of the waste body, and the internal strain. The existing monitoring methods are mainly divided into two categories: one is the point observation method, that is, the contact equipment such as displacement meters, inclinometers, and runoff sediment monitoring equipment are deployed on the waste body, combined with manual measurement methods for monitoring; the other is the regular mapping method, which uses drone photogrammetry and on-site investigation for regular monitoring.

[0004] Although the point observation method can obtain accurate local data, it has problems such as high cost and difficulty in selecting control points. Especially during the construction of waste slag stacking, it is difficult to find effective control points that meet the quantity requirements. The regular surveying method can obtain more comprehensive data, but it is limited by the frequency of drone surveying and the accuracy of ground correction, making it difficult to accurately grasp the changes in the waste slag body and provide timely feedback. Summary of the invention

[0005] In view of this, the present invention proposes a method for measuring soil erosion in a waste dump based on a three-dimensional model. By combining unmanned aerial vehicle photogrammetry and fixed tilt photogrammetry devices, high-precision, all-round and dynamic monitoring of the waste dump can be achieved.

[0006] The technical solution of the present invention is implemented as follows: The present invention provides a method for measuring soil erosion in a waste dump based on a three-dimensional model, comprising:

[0007] S1. Use drones to collect terrain data of the waste dump and obtain initial terrain images;

[0008] S2. Construct a three-dimensional vector model of the waste dump based on the initial terrain image and engineering design documents;

[0009] S3, setting a plurality of image control points around the waste dump, including at least two image control points at the lower edge of the slag pile and at least one image control point upstream;

[0010] S4. Install multiple fixed tilt photogrammetry devices with three-axis attitude control at predetermined locations above the waste dump according to the three-dimensional vector model;

[0011] S5. Use the oblique photogrammetry device to regularly collect multi-angle image data of the waste dump, obtain image data of different viewing angles by adjusting the posture angle of the device, and encrypt and match the collected multi-angle image data until the preset coordinate resolution requirement is met;

[0012] S6, processing the matched image data to generate three-dimensional point cloud data, and using the three-dimensional point cloud data to construct a three-dimensional irregular triangulated network model of the waste dump;

[0013] S7, calculating the elevation change within a preset time interval based on the preset model equation and the image control point data;

[0014] S8. Using numerical analysis methods, based on the three-dimensional irregular triangulated network model and elevation changes, extract the change data of the spoil body;

[0015] S9. Evaluate the stability of the waste dump based on the change data.

[0016] Based on the above solution, preferably, step S1 includes:

[0017] S11. Plan the drone route to ensure coverage of the entire waste dump area and its surroundings;

[0018] S12, setting the flight altitude and heading overlap of the drone;

[0019] S13, perform automatic flight and collect high-resolution orthophotos of the waste dump;

[0020] S14, preprocessing the collected images, including image matching, color correction and geometric correction;

[0021] S15. Generate a digital orthophoto map DOM and a digital elevation model DEM of the waste dump.

[0022] Based on the above solution, preferably, step S2 includes:

[0023] S21. Import digital orthophoto DOM and digital elevation model DEM;

[0024] S22. Import the waste dump engineering design documents, including the waste dump range, design slope ratio and design dump height parameters;

[0025] S23. Use GIS software to overlay and analyze engineering design parameters and terrain data;

[0026] S24, constructing a three-dimensional vector boundary of the waste dump;

[0027] S25. Generate a three-dimensional surface model of the waste dump, including the original terrain surface and the designed waste dump surface.

[0028] Based on the above solution, preferably, step S5 includes:

[0029] S51. Use a fixed three-axis gimbal equipped with an oblique photogrammetry device to conduct field measurements of the waste slag body;

[0030] S52, obtaining the gimbal attitude data, including the pitch axis deflection angle p, the roll axis deflection angle r, and the yaw axis deflection angle y of the gimbal;

[0031] S53, calculating the ground coordinates corresponding to the center point of the image using the gimbal posture data;

[0032] S54, determining the coordinates of the image projection center in the three-dimensional space;

[0033] S55, setting target coordinate resolution requirements;

[0034] S56, calculating the coordinate resolution of the near point, center point and far point of the image;

[0035] S57, determine whether the image coordinate resolution meets the target coordinate resolution requirement, if so, output the matched image data, if not, adjust the gimbal posture and return to step S51 to re-measure.

[0036] Based on the above scheme, preferably, the target coordinate resolution requirement is set to 5 cm, and the calculation formulas for the near point, center point, and far point coordinate resolutions are as follows:

[0037]

[0038] Among them, GSD top is the near point coordinate resolution, GSD mid is the center point coordinate resolution, GSD bottom is the far-point coordinate resolution, δ is the sensor pixel size; h is the installation height of the oblique photogrammeter; f is the camera focal length, α y is the inclination angle, β y =arctan(b / f) is half the viewing angle.

[0039] Based on the above solution, preferably, step S6 includes:

[0040] S61, generating a point cloud for the matched image data, filtering the generated three-dimensional point cloud data, and removing noise points and non-ground points;

[0041] S62, using an improved convex hull algorithm to extract the point cloud boundary and determine the outer contour of the triangulated network;

[0042] S63, based on the point cloud data within the boundary, applying the Delaunay triangulation algorithm to generate an initial triangulated network;

[0043] S64. Use a local optimization algorithm to adjust the initial triangulated network, optimize the shape and distribution of triangles, and obtain an optimized three-dimensional irregular triangulated network model.

[0044] Based on the above solution, preferably, the improved convex hull algorithm includes:

[0045] A. Project the 3D point cloud data onto the XY plane to obtain a 2D point set P = {p1, p2, ..., p n};

[0046] B. Use Graham scanning method to preliminarily construct the convex hull:

[0047] Find the point p0 with the smallest Y coordinate as the starting point;

[0048] Calculate the polar angle θ of other points relative to p0 i =arctan((y i -y0) / (x i -x0)), (x0,y0) are the coordinates of point p0, (x i ,y i ) are the coordinates of any other point;

[0049] Sort by polar angle and get the ordered point set S = {s1, s2, ..., s n-1};

[0050] Use the stack structure to process the points in S in sequence, keeping the points in the stack in a counterclockwise direction;

[0051] Among them, the process of judging the counterclockwise direction is:

[0052] For three points p1(x1,y1),p2(x2,y2),p3(x3,y3);

[0053] cross product=(x2-x1)(y3-y1)-(y2-y1)(x3-x1);

[0054] If cross product>0, it is counterclockwise;

[0055] C. Apply the α-shape algorithm to refine the convex hull:

[0056] Define the α value, which represents the parameter that controls the degree of boundary detail;

[0057] For each edge e(p k,p l ), where p k and p l are two adjacent points on the convex hull: calculate the radius r of the circle with e as the chord: r = |e| / (2*sin(α)). If r is less than the preset threshold R, keep the edge; otherwise, delete the edge.

[0058] Reconnect the retained edges to form an α-shape;

[0059] D. Local density analysis to identify concave features:

[0060] Calculate the average neighbor distance d of the point cloud avg ;

[0061] For each point p on the boundary of the α-shape m :

[0062] Calculate the p m The center is k*d, and the radius is avg The point density ρ within the circle m ;

[0063] If m Greater than the threshold ρ threshold , then we think that p m is a concave feature point;

[0064] For the identified concave feature points, the RanSaC algorithm is used to fit the local boundaries;

[0065] E. Merge similar boundary segments and smooth the boundary curves;

[0066] F. Map the optimized two-dimensional boundary points back to three-dimensional space to obtain the final point cloud boundary.

[0067] Based on the above solution, preferably, step S64 includes:

[0068] S641, the initial triangulation network is {T i}, set the quality threshold Q threshold ;

[0069] S642, for each triangle T i Calculate its quality index Q i ; Among them, the quality index includes the minimum angle quality index and the aspect ratio quality index. If Q i threshold , then mark the triangle as a low-quality triangle, and traverse to obtain a set of low-quality triangles

[0070] S643, select a low quality triangle ​

[0071] S644, Inspection If there is an edge e and an adjacent triangle Share, then and Pre-flip using the shared edge e as the flip axis and calculate the two new triangles after pre-flipping and Quality indicators and like Then perform edge flipping and go to step S647, otherwise, execute step S645;

[0072] S645, Construction The circumcircle of On the circumcircle Mark multiple target locations in Insert a new vertex v at the current center position of Delete the three vertices of , and The three vertices of are connected to form three new triangles Calculate the quality index of the new triangle Then take the average value to get like Then keep three new triangles and use the vertices of the three new triangles as the original Otherwise, select another target position and insert the new vertex repeatedly. If the quality of all target positions is not improved, keep

[0073] S646, for Each vertex of Calculate the centroid g of its adjacent vertices, initialize the smoothing coefficient λ, and Move a small step in the direction of g: v'=v+λ(gv), calculate the average quality Q of the adjacent triangles before and after the move before and Q after , if Q before after , then accept this move, otherwise λ is halved and the attempt is repeated until λ<0.05 or an acceptable move is found;

[0074] S647, repeat steps S643-S646 until After traversal, update the triangulation network {T i};

[0075] S648, repeat steps S642-S647 until one of the stop conditions is met, the stop condition includes that the quality of all triangles is higher than the threshold Q threshold , reaching the preset maximum number of iterations N​max , the quality improvement of k consecutive iterations is less than the preset value ε;

[0076] S649. Perform global vertex smoothing and fine-tune the entire triangulation network.

[0077] Based on the above solution, preferably, in step S7, the preset model equation is as follows:

[0078] sin m sin n·x+sin m cos n·y+cos n·z=L;

[0079] Among them, x, y, z are spatial coordinates, m is the inclination angle of the spoil body, n is the inclination angle of the spoil body, and L is a constant. The determination method is as follows:

[0080] The distance from image control point 1 is L1, the distance from image control point 2 is L2, the distance from image control point 3 is L3, and so on. The distance from image control point n is L n , after spatial triangulation, determine the constant L;

[0081] Assume that the interval between two measurements of the oblique photogrammetry equipment is Δt; the constant after Δt is L Δt ; At the same time, the reference coordinates remain unchanged, so the x and y coordinates remain unchanged; the elevation changes are:

[0082]

[0083] Where z(t+Δt) is the elevation at time t+Δt, and z(t) is the elevation at time t.

[0084] Based on the above solution, preferably, step S8 includes:

[0085] S81, combining the optimized three-dimensional irregular triangulated network model with the original elevation data to generate a three-dimensional digital model including terrain features;

[0086] S82, discretize the three-dimensional scene into a regular grid, and for each grid point (i, j), obtain the calculated elevation change Δz(i, j), Δz(i, j) = z(i, j, t+Δt)-z(i, j, t);

[0087] S83. Analyze the changes of waste slag:

[0088] If Δz(i,j)=0, the grid point is considered to be stable;

[0089] If Δz(i,j)>0, it is considered that the grid point is piled up;

[0090] If Δz(i,j)<0, the grid point is considered to be eroded;

[0091] S84. Based on the analysis results, a spatial distribution map of the changes in the spoil body is generated, marking the accumulation, erosion and stability areas.

[0092] Compared with the prior art, the present invention has the following beneficial effects:

[0093] (1) By combining drone photogrammetry and fixed tilt photogrammetry devices, a three-dimensional vector model is constructed and multi-angle image data is collected to achieve high-precision, all-round, and dynamic monitoring of the waste dump. This method not only improves the continuity and accuracy of monitoring, but also can capture terrain changes in a timely manner and accurately calculate the amount of accumulation and erosion. By setting multiple image control points and using preset model equations, this method can effectively evaluate the stability of the waste dump and provide reliable data support for engineering management and safety decisions, while reducing labor costs and improving monitoring efficiency.

[0094] (2) By generating point clouds and filtering the matched image data, and applying the improved convex hull algorithm and Delaunay triangulation algorithm, this method can generate a high-quality three-dimensional irregular triangulated network model. This optimized model can more accurately characterize the terrain characteristics of the waste dump, improving the accuracy and reliability of subsequent change analysis.

[0095] (3) The improved convex hull algorithm is used to extract the point cloud boundary, which significantly improves the definition accuracy of the external contour of the waste dump. By combining the Graham scanning method and the α-shape algorithm, the algorithm can accurately capture complex terrain boundaries, including concave features. Through local density analysis, the algorithm can identify and retain important terrain details, avoiding the oversimplification that may be caused by the traditional convex hull method. This precise boundary extraction provides a reliable foundation for subsequent triangulation, ensuring that the generated three-dimensional model can truly reflect the actual shape of the waste dump.

[0096] (4) The quality of the irregular triangulated network model was greatly improved by applying the improved Delaunay triangulation algorithm and local optimization algorithm. The optimization process identified and improved low-quality triangles by setting quality thresholds, including adjusting the shape, size, and distribution of triangles. This method effectively reduced the appearance of narrow or deformed triangles and improved the geometric quality and expressiveness of the overall model. The optimized triangulated network model can more accurately characterize the subtle changes and complex terrain features of the waste dump surface, providing a more accurate data basis for subsequent soil erosion and accumulation calculations.

[0097] (5) Based on the preset model equations and the image control point data, this method can accurately quantify the vertical changes of the surface of the waste dump at different time points. By considering the inclination and dip angle of the waste dump, this method can more accurately reflect the changes in complex terrain. BRIEF DESCRIPTION OF THE DRAWINGS

[0098] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative work.

[0099] Figure 1 A flow chart of a method for measuring soil erosion in a waste dump based on a three-dimensional model provided in an embodiment of the present invention;

[0100] Figure 2 A technical implementation diagram of a method for measuring soil erosion in a waste dump based on a three-dimensional model provided by an embodiment of the present invention;

[0101] Figure 3 A schematic diagram of the layout provided for an embodiment of the present invention. DETAILED DESCRIPTION

[0102] The following will be combined with the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0103] like Figure 1 and Figure 2 As shown, the present invention provides a method for measuring soil erosion in a waste dump based on a three-dimensional model, comprising:

[0104] S1. Use drones to collect terrain data of the waste dump and obtain initial terrain images.

[0105] S2. Construct a three-dimensional vector model of the waste dump based on the initial terrain image and engineering design documents.

[0106] S3. Setting a plurality of image control points around the waste dump, including at least two image control points at the lower edge of the slag dump and at least one image control point upstream.

[0107] S4. According to the three-dimensional vector model, multiple fixed-tilt photogrammetry devices with three-axis attitude control are installed at predetermined positions above the waste dump.

[0108] S5. Use the oblique photogrammetry device to regularly collect multi-angle image data of the waste dump, obtain image data of different viewing angles by adjusting the posture angle of the device, and encrypt and match the collected multi-angle image data until the preset coordinate resolution requirement is met.

[0109] S6. Process the matched image data to generate three-dimensional point cloud data, and use the three-dimensional point cloud data to construct a three-dimensional irregular triangulated network model of the waste dump.

[0110] S7. Calculate the elevation change within a preset time interval based on the preset model equation and the image control point data.

[0111] S8. Use numerical analysis methods to extract the change data of the waste slag body based on the three-dimensional irregular triangulated network model and elevation changes.

[0112] S9. Evaluate the stability of the waste dump based on the change data.

[0113] Specifically, in one embodiment of the present invention, step S1 includes:

[0114] S11. Plan the drone route to ensure coverage of the entire waste dump area and its surroundings.

[0115] Design a comprehensive flight path to ensure that the drone can collect complete image data of the entire waste dump and its surrounding area. Route planning needs to take into account the shape, size and terrain characteristics of the waste dump, as well as the surrounding environment.

[0116] Using DJI GS Pro or similar flight planning software, set the spoil dump boundary as a polygonal area and extend 50 meters outside the boundary as a buffer zone. The flight direction is set to be perpendicular to the main slope of the spoil dump to ensure the best terrain capture effect. The flight spacing is set to 60 meters, which is calculated based on the camera parameters used and the desired overlap.

[0117] S12. Set the flight altitude and heading overlap of the drone.

[0118] The flight altitude affects the resolution and coverage of the image, while the heading overlap affects the image stitching quality and the accuracy of 3D reconstruction. For example, the flight altitude is set to 120 meters, the heading overlap is set to 80%, the lateral overlap is set to 70%, and the flight speed is set to 8m / s.

[0119] S13. Execute automatic flight route to collect high-resolution orthophotos of the waste dump.

[0120] Use the drone's built-in RTK positioning system to ensure centimeter-level positioning accuracy. Set the camera to shoot in equal time interval mode with an interval of 2 seconds to ensure sufficient front and back overlap. Monitor the drone status and image acquisition in real time during the flight.

[0121] S14. Preprocessing the collected images, including image matching, color correction and geometric correction.

[0122] First, image matching was performed, with the number of key points extracted set to 40,000 points / image and the matching threshold set to 0.7. Then, color correction was performed, using the histogram equalization method to ensure that the brightness and contrast of all images were consistent. Geometric correction was performed using the adaptive bundle adjustment method, with the reprojection error threshold set to 0.5 pixels.

[0123] S15. Generate a digital orthophoto map DOM and a digital elevation model DEM of the waste dump.

[0124] The ground sampling distance (GSD) was set to 5 cm to generate a digital orthophoto map (DOM). For the digital elevation model (DEM), the grid spacing was set to 10 cm and the inverse distance weighted interpolation method was used. Finally, the DOM and DEM were exported to GeoTIFF format to ensure that the geographic reference information was correct and the projection coordinate system was set to the commonly used local coordinate system.

[0125] Specifically, in one embodiment of the present invention, step S2 includes:

[0126] S21. Import digital orthophoto DOM and digital elevation model DEM.

[0127] Use ArcGIS Pro or QGIS and other GIS software to import the DOM and DEM files generated in the previous step. The DOM file format is GeoTIFF with a resolution of 5 cm / pixel; the DEM file is also in GeoTIFF format with a resolution of 10 cm / pixel.

[0128] S22. Import the waste dump engineering design documents, including the waste dump range, design slope ratio and design dump height parameters.

[0129] Import the engineering design files of the waste dump into the GIS software. These files include design drawings in CAD format, which need to be converted to match the coordinate system of DOM and DEM. The specific imported parameters include: waste dump range, design slope ratio, and design dump height.

[0130] S23. Use GIS software to overlay and analyze engineering design parameters and terrain data.

[0131] The analysis process is:

[0132] Superimpose the designed slag dump range with DOM to check the difference between the actual slag dump situation and the design.

[0133] Using DEM data, a slope raster layer is created to analyze the consistency between the actual slope and the designed slope ratio.

[0134] Perform elevation analysis to compare the actual elevation in the DEM with the designed pile height.

[0135] S24. Construct the three-dimensional vector boundary of the waste dump.

[0136] Based on the results of the overlay analysis, the 3D analysis tool of the GIS software was used to construct a three-dimensional vector boundary: the actual boundary line of the waste dump was extracted using DOM and DEM data. The boundary line was converted into a 3D feature and assigned an elevation value. The boundary was smoothed using the Bezier curve algorithm and the smoothing tolerance was set to 0.5 meters.

[0137] S25. Generate a three-dimensional surface model of the waste dump, including the original terrain surface and the designed waste dump surface.

[0138] Generate a TIN (Triangulated Irregular Network) model of the original terrain using DEM data. Create an ideal TIN model of the designed slag surface based on the design parameters. Merge the two TIN models to generate a comprehensive 3D surface model. Use the 3D analysis tools of the GIS software to calculate the volume difference between the actual terrain and the designed terrain.

[0139] Specifically, Figure 3 As shown, in one embodiment of the present invention, step S3 includes:

[0140] Set up no less than 3 image control points. Set up at least 2 image control points at the lower edge of the slag pile. Set up at least 1 image control point upstream of the waste dump.

[0141] Choose a stable, easily identifiable location around the dump. Ensure that the GCP will not be affected or buried by subsequent dumping activities. Choose a location that will remain stable for a long time to allow for long-term monitoring.

[0142] GCPs are marked with survey marks. High-precision GNSS RTK is used to measure the coordinates of GCPs. The number and coordinates (X, Y, Z) of each GCP are recorded in detail. GCP information is integrated into the GIS database for subsequent processing and analysis.

[0143] Specifically, in one embodiment of the present invention, step S4 includes:

[0144] Based on the three-dimensional vector model constructed in step S2, select a location on the original mountain that is higher than the designed dump height of the waste dump. Ensure that the selected location is not blocked by the mountain and can clearly observe the entire waste dump area. Select at least 2 different installation locations to ensure that the viewing angles complement each other and avoid monitoring blind spots.

[0145] Use a fixed-tilt photogrammetry device with a three-axis gimbal attitude control. The device should include a high-resolution camera, a three-axis gimbal, and an accurate attitude sensor. Establish a stable mounting platform at the selected location to ensure long-term stable operation of the device. When installing the device, it should face the waste dump with the lens at a depression or flat angle. Accurately adjust the device orientation to ensure that it covers the waste dump area and an area that is proportional to 20% of the area. Configure camera parameters, including focal length, aperture, ISO, etc., to adapt to different light conditions. Set the gimbal control parameters to ensure that the pitch axis deflection angle p, roll axis deflection angle r, and heading axis deflection angle y can be accurately adjusted.

[0146] Specifically, in one embodiment of the present invention, step S5 includes:

[0147] S51. Use a fixed three-axis gimbal equipped with an oblique photogrammetry device to conduct on-site measurement of the waste slag body.

[0148] S52, obtaining the gimbal attitude data, including the pitch axis deflection angle p, the roll axis deflection angle r, and the yaw axis deflection angle y of the gimbal.

[0149] S53. Calculate the ground coordinates corresponding to the center point of the image using the gimbal posture data.

[0150] S54, determining the coordinates of the image projection center in the three-dimensional space.

[0151] S55. Set the target coordinate resolution requirement.

[0152] S56. Calculate the coordinate resolution of the near point, center point, and far point of the image.

[0153] S57, determine whether the image coordinate resolution meets the target coordinate resolution requirement, if so, output the matched image data, if not, adjust the gimbal posture and return to step S51 to re-measure.

[0154] In this embodiment, the oblique photogrammetry device carried by the fixed three-axis gimbal is started to read the initial posture data of the gimbal, including: the pitch axis deflection angle p, the roll axis deflection angle r and the heading axis deflection angle y of the gimbal.

[0155] Calculate the ground coordinates corresponding to the center point of the image:

[0156]

[0157] Where [XYZ] T is the ground coordinate of the center point of the image, f is the focal length, (X s ,Y s ,Z s ) is the exposure point coordinate, R is the rotation matrix, and the gimbal attitude angle is calculated as follows:

[0158]

[0159] Combining the above two equations, the ground coordinates of the image center point can be simplified as follows:

[0160]

[0161] The coordinates of the image projection center in three-dimensional space are calculated using the GNSS antenna phase center coordinates and the gimbal attitude data.

[0162] The camera projection center, the gimbal GNSS antenna phase center, the three-dimensional coordinates of the ground object point and the image point coordinates meet the requirements.

[0163]

[0164] Where, [xy -f] T represents the coordinates of the image point in the image space coordinate system, λ is the denominator of the photographic scale; R is the rotation matrix; is the coordinate of the GNSS antenna phase center, is the coordinate of the GNSS antenna phase center in the image space coordinate system; [XYZ] T is the ground point coordinate corresponding to the image point.

[0165] When the adjustment converges, the image projection center can be calculated from the position of the GNSS antenna phase center. See the following formula for details.

[0166]

[0167] In the formula, [X s Y s Z s ] T are the coordinates of the image projection center in three-dimensional space.

[0168] Set the target ground sampling distance (GSD) to 5 cm.

[0169] Use the following formula to calculate the coordinate resolution of the near point, center point, and far point:

[0170]

[0171] Among them, GSD top is the near point coordinate resolution, GSD mid is the center point coordinate resolution, GSD bottom is the far-point coordinate resolution, δ is the sensor pixel size; h is the installation height of the oblique photogrammeter; f is the camera focal length, α y is the inclination angle, β y =arctan(b / f) is half the viewing angle.

[0172] Determine whether the calculated resolution meets the target requirement (5cm). If it does not meet the requirement, make the following adjustments:

[0173] a) Fine-tune the gimbal posture (adjust p, r, y angles).

[0174] b) If necessary, adjust camera parameters (such as focal length, aperture, etc.).

[0175] Repeat the above process until the target resolution requirement is reached.

[0176] On the basis of meeting the resolution requirements, control the gimbal to shoot at multiple angles: a) Set different pitch and heading angles. b) Shoot multiple images at each angle to ensure sufficient overlap. Record the exact shooting time of each image and the corresponding gimbal attitude data.

[0177] Perform preliminary processing on the collected multi-angle image data:

[0178] a) Image matching: Use feature point matching algorithms (such as SIFT or SURF). b) Coarse matching: Quickly identify overlapping areas. c) Fine matching: Accurately calculate the 3D coordinates of corresponding points.

[0179] Check the image quality in real time, including clarity, exposure, etc. Ensure that the image overlap meets the requirements (front-to-back overlap ≥ 80%, lateral overlap ≥ 60%).

[0180] Specifically, in one embodiment of the present invention, step S6 includes:

[0181] S61 , generating a point cloud for the matched image data, and filtering the generated three-dimensional point cloud data to remove noise points and non-ground points.

[0182] Use multi-view stereo vision algorithms (such as Structure from Motion, SfM) to process the matched image data.

[0183] Extract feature points and perform feature matching to establish correspondence between images.

[0184] Sparse reconstruction is performed to obtain preliminary 3D point cloud and camera pose.

[0185] Use dense matching algorithms (such as Semi-GlobalMatching, SGM) to generate high-density point clouds.

[0186] Filter the generated 3D point cloud data:

[0187] Apply statistical outlier filtering to remove noise points.

[0188] Use height threshold filtering to remove points that are obviously not on the ground.

[0189] Apply a segmentation algorithm to identify and remove non-ground points.

[0190] S62. Use the improved convex hull algorithm to extract the point cloud boundary and determine the outer contour of the triangulated network.

[0191] In this embodiment, the improved convex hull algorithm includes:

[0192] A. Project the 3D point cloud data onto the XY plane to obtain a 2D point set P = {p1, p2, ..., p n};

[0193] B. Use Graham scanning method to preliminarily construct the convex hull:

[0194] Find the point p0 with the smallest Y coordinate as the starting point;

[0195] Calculate the polar angle θ of other points relative to p0 i =arctan((y i -y0) / (x i -x0)), (x0,y0) are the coordinates of point p0, (x i ,y i ) are the coordinates of any other point;

[0196] Sort by polar angle and get the ordered point set S = {s1, s2, ..., s n-1};

[0197] Use the stack structure to process the points in S in sequence, keeping the points in the stack in a counterclockwise direction;

[0198] Among them, the process of judging the counterclockwise direction is:

[0199] For three points p1(x1,y1),p2(x2,y2),p3(x3,y3);

[0200] cross product=(x2-x1)(y3-y1)-(y2-y1)(x3-x1);

[0201] If cross product>0, it is counterclockwise;

[0202] C. Apply the α-shape algorithm to refine the convex hull:

[0203] Define the α value, which represents the parameter that controls the degree of boundary detail;

[0204] For each edge e(p k ,p l ), where p k and p lare two adjacent points on the convex hull: calculate the radius r of the circle with e as the chord: r = |e| / (2*sin(α)). If r is less than the preset threshold R, keep the edge; otherwise, delete the edge.

[0205] Reconnect the retained edges to form an α-shape;

[0206] D. Local density analysis to identify concave features:

[0207] Calculate the average neighbor distance d of the point cloud avg ;

[0208] For each point p on the boundary of the α-shape m :

[0209] Calculate the p m The center is k*d, and the radius is avg The point density ρ within the circle m ;

[0210] If m Greater than the threshold ρ threshold , then we think that p m is a concave feature point;

[0211] For the identified concave feature points, the RanSaC algorithm is used to fit the local boundaries;

[0212] E. Merge similar boundary segments and smooth the boundary curves;

[0213] F. Map the optimized two-dimensional boundary points back to three-dimensional space to obtain the final point cloud boundary.

[0214] The point density calculation formula is: Among them, n m is the number of points in the circle, and k is the density calculation radius coefficient. When smoothing the boundary curve, B-spline interpolation is used.

[0215] In this embodiment, considering that the traditional convex hull algorithm cannot handle depressions and complex boundaries in the terrain, which will cause the generated boundary to be too simple and unable to accurately reflect the terrain characteristics, the α-shape algorithm is introduced to refine the convex hull. By adjusting the parameter α, the detail level of the boundary can be controlled, thereby more accurately describing the true boundary of the terrain, and local density analysis is used to identify concave features, solving the problem that the traditional convex hull cannot identify the depressed area. The improved convex hull algorithm can generate more accurate terrain boundaries, improve the accuracy and reliability of the three-dimensional model, and is suitable for the feature description of the complex terrain of the waste dump.

[0216] S63 generates an initial triangulated network based on the point cloud data within the boundary using the Delaunay triangulation algorithm.

[0217] Project the 3D point cloud to the XY plane. Generate a 2D triangulated network using the standard Delaunay triangulation algorithm. Map the 2D triangulated network back to 3D space and restore the Z coordinate value.

[0218] S64. Use a local optimization algorithm to adjust the initial triangulated network, optimize the shape and distribution of triangles, and obtain an optimized three-dimensional irregular triangulated network model.

[0219] In this embodiment, step S64 includes:

[0220] S641, the initial triangulation network is {T i}, set the quality threshold Q threshold ;

[0221] S642, for each triangle T i Calculate its quality index Q i ; Among them, the quality index includes the minimum angle quality index and the aspect ratio quality index. If Q i threshold , then mark the triangle as a low-quality triangle, and traverse to obtain a set of low-quality triangles

[0222] S643, select a low quality triangle

[0223] S644, Inspection If there is an edge e and an adjacent triangle Share, then and Pre-flip using the shared edge e as the flip axis and calculate the two new triangles after pre-flipping and Quality indicators and like Then perform edge flipping and go to step S647, otherwise, execute step S645;

[0224] S645, Construction The circumcircle of On the circumcircle Mark multiple target locations in Insert a new vertex v at the current center position of Delete the three vertices of , and The three vertices of are connected to form three new triangles Calculate the quality index of the new triangle Then take the average value to get like ​Then keep three new triangles and use the vertices of the three new triangles as the original Otherwise, select another target position and insert the new vertex repeatedly. If the quality of all target positions is not improved, keep

[0225] S646, for Each vertex of Calculate the centroid g of its adjacent vertices, initialize the smoothing coefficient λ, and Move a small step in the direction of g: v'=v+λ(gv), calculate the average quality Q of the adjacent triangles before and after the move before and Q after , if Q before after , then accept this move, otherwise λ is halved and the attempt is repeated until λ<0.05 or an acceptable move is found;

[0226] S647, repeat steps S643-S646 until After traversal, update the triangulation network {T i};

[0227] S648, repeat steps S642-S647 until one of the stop conditions is met, the stop condition includes that the quality of all triangles is higher than the threshold Q threshold , reaching the preset maximum number of iterations N max , the quality improvement of k consecutive iterations is less than the preset value ε;

[0228] S649. Perform global vertex smoothing and fine-tune the entire triangulation network.

[0229] This embodiment uses the Delaunay triangulation algorithm to generate the initial triangulated network, but the initial triangulated network may contain low-quality triangles, such as slender or distorted triangles. Low-quality triangles will cause the model to be inaccurate in geometry and topology, affecting subsequent analysis and calculation. Therefore, local optimization processes such as edge flipping and vertex insertion are introduced to improve the quality of triangles, while performing global vertex smoothing to avoid falling into local optimality. This optimization process ensures the quality and stability of the model and is suitable for the dynamic monitoring needs of the waste dump.

[0230] Specifically, in one embodiment of the present invention, step S7 includes:

[0231] 1. The preset model equation is as follows:

[0232] sin m sin n·x+sin m cos n·y+cos n·z=L;

[0233] ​Among them, x, y, z are spatial coordinates, m is the inclination angle of the spoil body, n is the inclination angle of the spoil body, and L is a constant.

[0234] 2. Using the image control points arranged in step S3 (at least 3, including at least 2 at the lower edge of the slag pile and at least 1 upstream of the waste dump), accurately measure each image control point to obtain its three-dimensional coordinates (x, y, z).

[0235] 3. For each image control point, calculate its distance to the model plane: the distance to image control point 1 is L1, the distance to image control point 2 is L2, the distance to image control point 3 is L3, and so on. The distance to image control point n is L n , use these distance values ​​to perform aerial triangulation adjustment and determine the optimal constant L value.

[0236] 4. Use the determined L value and substitute it into the preset model equation to establish the model equation of the waste slag body at the initial moment.

[0237] 5. Set the measurement time interval: Determine the time interval Δt between two measurements (such as 1 day, 1 week or 1 month, depending on monitoring needs).

[0238] 6. Measurement after Δt time: After Δt time, repeat steps 2-4 to obtain a new constant value L Δt .

[0239] 7. If the reference coordinates remain unchanged, the x and y coordinates remain unchanged; the elevation changes are:

[0240]

[0241] Where z(t+Δt) is the elevation at time t+Δt, and z(t) is the elevation at time t.

[0242] Specifically, in one embodiment of the present invention, step S8 includes:

[0243] S81, combining the optimized three-dimensional irregular triangulated network model with the original elevation data to generate a three-dimensional digital model including terrain features;

[0244] S82, discretize the three-dimensional scene into a regular grid, and for each grid point (i, j), obtain the calculated elevation change Δz(i, j), Δz(i, j) = z(i, j, t+Δt)-z(i, j, t);

[0245] S83. Analyze the changes of waste slag:

[0246] If Δz(i,j)=0, the grid point is considered to be stable;

[0247] If Δz(i,j)>0, it is considered that the grid point is piled up;

[0248] If Δz(i,j)<0, the grid point is considered to be eroded;

[0249] S84. Based on the analysis results, a spatial distribution map of the changes in the spoil body is generated, marking the accumulation, erosion and stability areas.

[0250] In this embodiment, the three-dimensional irregular triangulated network model (TIN) is discretized into a regular grid (rasterization) to facilitate the calculation of elevation changes. The grid size is 5*5. For each grid point (i, j), the TIN model is used for spatial interpolation to obtain the elevation value z(i, j, t) at time t. The elevation value z(i, j, t+Δt) at time t+Δt is calculated using the elevation change data. After that, the elevation change is calculated: Δz(i, j) = z(i, j, t+Δt)-z(i, j, t).

[0251] After analyzing the changes in the spoil body, the change data can be quantified. The process is as follows:

[0252] Calculate the total accumulation: sum(Δz(i,j))for allΔz(i,j)>0;

[0253] Calculate the total erosion: sum(|Δz(i,j)|)for allΔz(i,j)<0;

[0254] Calculate the net change: sum(Δz(i,j))for all(i,j).

[0255] Spatial distribution analysis: Generate a spatial distribution map of the change type, mark the accumulation, erosion and stability areas, and calculate the area proportion of each type of area.

[0256] Hot spot analysis: Identify the areas with the most significant changes (places with the most severe accumulation or erosion) and calculate the rate of change and range of these hot spots. Rate of change DR = (Δz / Δt) * 100%.

[0257] Slope change analysis: Using the initial and final elevation data, the slope change is calculated for each grid point to identify areas with significant slope changes, which may indicate potential instability. Slope change SC = arctan((z2-z1) / d)-arctan((z1-z0) / d), where z0, z1, z2 are the elevation values ​​of three adjacent time points and d is the horizontal distance.

[0258] Volume change calculation: Using the elevation change data, calculate the volume change of the entire monitoring area, and calculate the volume change of the accumulation area and the erosion area respectively. Volume change VC = Σ(Δz*GridArea), where GridArea is the area of ​​each grid.

[0259] Statistical analysis: Calculate the statistical characteristics of elevation changes (mean, standard deviation, maximum, minimum, etc.) and conduct trend analysis.

[0260] Result visualization: Generate contour maps or heat maps of elevation changes and create 3D visualization models to intuitively display the changes.

[0261] Specifically, in one embodiment of the present invention, step S9 includes:

[0262] The spoil body change data extracted in step S8 are summarized, including elevation change, volume change, slope change and other data.

[0263] Stability index calculation:

[0264] Calculate the overall deformation rate: the average rate of change in elevation across the entire monitoring area. Calculate the local deformation rate: identify the area with the most significant change and calculate its rate of change. Calculate the slope change rate: pay special attention to areas with increased slope. Calculate the volume change rate: analyze the volume change ratio of accumulation and erosion.

[0265] Set warning thresholds for various indicators (based on engineering experience and geological conditions) and compare the calculated indicators with the preset thresholds.

[0266] Identify areas that exceed warning thresholds, paying special attention to areas with severe erosion at the foot of the slope or obvious accumulation at the top of the slope.

[0267] Combined with meteorological data (such as rainfall), the impact of environmental factors on stability is analyzed, and the impact of human activities (such as new slag piles, changes in drainage facilities, etc.) is considered.

[0268] Based on the above analysis, a qualitative assessment of the overall stability of the waste dump was conducted. The assessment results include: stable, basically stable, potentially unstable, unstable and other levels.

[0269] Conduct a detailed local stability analysis on the identified high-risk areas, and based on the assessment results, classify early warning levels for different areas of the waste dump, such as green (safe), yellow (need to pay attention), orange (need to be vigilant), red (high risk), etc.

[0270] Prepare a comprehensive report containing all analysis results, assessment conclusions and early warning information. The report should include data charts, spatial distribution maps and assessment basis.

[0271] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the present invention should be included in the protection scope of the present invention.

Claims

1. A method for measuring soil erosion in a waste dump based on a three-dimensional model, characterized in that: include: S1. Use drones to collect terrain data of the waste dump and obtain initial terrain images; S2. Construct a three-dimensional vector model of the waste dump based on the initial terrain image and engineering design documents; S3. Setting a number of image control points around the waste dump, including at least two image control points at the lower edge of the slag dump and at least one image control point upstream; S4. According to the three-dimensional vector model, several fixed tilt photogrammetry devices with three-axis attitude control are installed at predetermined positions above the waste dump; S5. Use the oblique photogrammetry device to regularly collect multi-angle image data of the waste dump, obtain image data of different viewing angles by adjusting the posture angle of the device, and encrypt and match the collected multi-angle image data until the preset coordinate resolution requirement is met; S6, processing the matched image data to generate three-dimensional point cloud data, and using the three-dimensional point cloud data to construct a three-dimensional irregular triangulated network model of the waste dump; S7, calculating the elevation change within a preset time interval based on the preset model equation and the image control point data; S8. Using numerical analysis methods, based on the three-dimensional irregular triangulated network model and elevation changes, extract the change data of the spoil body; S9. Evaluate the stability of the waste dump based on the change data.

2. A method for measuring soil erosion in a waste dump based on a three-dimensional model as claimed in claim 1, characterized in that: Step S1 includes: S11. Plan the drone route to ensure coverage of the entire waste dump area and its surroundings; S12, setting the flight altitude and heading overlap of the drone; S13, perform automatic flight and collect high-resolution orthophotos of the waste dump; S14, preprocessing the collected images, including image matching, color correction and geometric correction; S15. Generate a digital orthophoto map DOM and a digital elevation model DEM of the waste dump.

3. A method for measuring soil erosion in a waste dump based on a three-dimensional model as claimed in claim 2, characterized in that: Step S2 includes: S21. Import digital orthophoto DOM and digital elevation model DEM; S22. Import the waste dump engineering design documents, including the waste dump range, design slope ratio and design dump height parameters; S23. Use GIS software to overlay and analyze engineering design parameters and terrain data; S24, constructing a three-dimensional vector boundary of the waste dump; S25. Generate a three-dimensional surface model of the waste dump, including the original terrain surface and the designed waste dump surface.

4. The method for measuring soil erosion in a waste dump based on a three-dimensional model as claimed in claim 1, characterized in that: Step S5 includes: S51. Use a fixed three-axis gimbal equipped with an oblique photogrammetry device to conduct field measurements of the waste slag body; S52, obtaining the gimbal attitude data, including the pitch axis deflection angle p, the roll axis deflection angle r, and the yaw axis deflection angle y of the gimbal; S53, calculating the ground coordinates corresponding to the center point of the image using the gimbal posture data; S54, determining the coordinates of the image projection center in the three-dimensional space; S55, setting target coordinate resolution requirements; S56, calculating the coordinate resolution of the near point, center point and far point of the image; S57, determine whether the image coordinate resolution meets the target coordinate resolution requirement, if so, output the matched image data, if not, adjust the gimbal posture and return to step S51 to re-measure.

5. A method for measuring soil erosion in a waste dump based on a three-dimensional model as claimed in claim 4, characterized in that: The target coordinate resolution requirement is set to 5cm. The calculation formulas for the near point, center point, and far point coordinate resolutions are as follows: Among them, GSD top is the near point coordinate resolution, GSD mid is the center point coordinate resolution, GSD bottom is the far-point coordinate resolution, δ is the sensor pixel size; h is the installation height of the oblique photogrammeter; f is the camera focal length, α y is the inclination angle, β y =arctan(b / f) is half the viewing angle.

6. The method for measuring soil erosion in a waste dump based on a three-dimensional model as claimed in claim 2, characterized in that: Step S6 includes: S61, generating a point cloud for the matched image data, filtering the generated three-dimensional point cloud data, and removing noise points and non-ground points; S62, using an improved convex hull algorithm to extract the point cloud boundary and determine the outer contour of the triangulated network; S63, based on the point cloud data within the boundary, applying the Delaunay triangulation algorithm to generate an initial triangulated network; S64. Use a local optimization algorithm to adjust the initial triangulated network, optimize the shape and distribution of triangles, and obtain an optimized three-dimensional irregular triangulated network model.

7. A method for measuring soil erosion in a waste dump based on a three-dimensional model as claimed in claim 6, characterized in that: The improved convex hull algorithm includes: A. Project the 3D point cloud data onto the XY plane to obtain a 2D point set P = {p1, p2, ..., p n }; B. Use Graham scanning method to preliminarily construct the convex hull: Find the point p0 with the smallest Y coordinate as the starting point; Calculate the polar angle θ of other points relative to p0 i =arctan((y i -y0) / (x i -x0)), (x0,y0) are the coordinates of point p0, (x i ,y i ) are the coordinates of any other point; Sort by polar angle and get the ordered point set S = {s1, s2, ..., s n-1 }; Use the stack structure to process the points in S in sequence, keeping the points in the stack in a counterclockwise direction; Among them, the process of judging the counterclockwise direction is: For three points p1(x1,y1),p2(x2,y2),p3(x3,y3); cross product=(x2-x1)(y3-y1)-(y2-y1)(x3-x1); If the cross product is greater than 0, it is counterclockwise; C. Apply the α-shape algorithm to refine the convex hull: Define the α value, which represents the parameter that controls the degree of boundary detail; For each edge e(p k ,p l ), where p k and p l are two adjacent points on the convex hull: calculate the radius r of the circle with e as the chord: r = |e| / (2*sin(α)). If r is less than the preset threshold R, keep the edge; otherwise, delete the edge. Reconnect the retained edges to form an α-shape; D. Local density analysis to identify concave features: Calculate the average neighbor distance d of the point cloud avg ; For each point p on the boundary of the α-shape m : Calculate the p m The center is k*d, and the radius is avg The point density ρ within the circle m ; If m Greater than the threshold ρ threshold , then we think that p m is a concave feature point; For the identified concave feature points, the RanSaC algorithm is used to fit the local boundaries; E. Merge similar boundary segments and smooth the boundary curves; F. Map the optimized two-dimensional boundary points back to three-dimensional space to obtain the final point cloud boundary.

8. The method for measuring soil erosion in a waste dump based on a three-dimensional model as claimed in claim 6, characterized in that: Step S64 includes: S641, the initial triangulation network is {T i }, set the quality threshold Q threshold ; S642, for each triangle T i Calculate its quality index Q i ; Among them, the quality index includes the minimum angle quality index and the aspect ratio quality index. If Q i <Q threshold , then mark the triangle as a low-quality triangle, and traverse to obtain a set of low-quality triangles S643, select a low quality triangle S644, Inspection If there is an edge e and an adjacent triangle Share, then and Pre-flip using the shared edge e as the flip axis and calculate the two new triangles after pre-flipping and Quality indicators and like Then perform edge flipping and go to step S647, otherwise, execute step S645; S645, Construction The circumcircle of On the circumcircle Mark multiple target locations in Insert a new vertex v at the current center position of Delete the three vertices of , and The three vertices of are connected to form three new triangles Calculate the quality index of the new triangle Then take the average value to get like Then keep three new triangles and use the vertices of the three new triangles as the original Otherwise, select another target position and insert the new vertex repeatedly. If the quality of all target positions is not improved, keep S646, for Each vertex of Calculate the centroid g of its adjacent vertices, initialize the smoothing coefficient λ, and Move a small step in the direction of g: v'=v+λ(gv), calculate the average quality Q of the adjacent triangles before and after the move before and Q after , if Q before <Q after , then accept this move, otherwise λ is halved and the attempt is repeated until λ < 0.05 or an acceptable move is found; S647, repeat steps S643-S646 until After traversal, update the triangulation network {T i }; S648, repeat steps S642-S647 until one of the stop conditions is met, the stop condition includes that the quality of all triangles is higher than the threshold Q threshold , reaching the preset maximum number of iterations N max , the quality improvement of k consecutive iterations is less than the preset value ε; S649. Perform global vertex smoothing and fine-tune the entire triangulation network.

9. The method for measuring soil erosion in a waste dump based on a three-dimensional model as claimed in claim 6, characterized in that: In step S7, the preset model equation is as follows: sin m sin n·x+sin m cos n·y+cos n·z=L; Among them, x, y, z are spatial coordinates, m is the inclination angle of the spoil body, n is the inclination angle of the spoil body, and L is a constant. The determination method is as follows: The distance from image control point 1 is L1, the distance from image control point 2 is L2, the distance from image control point 3 is L3, and so on. The distance from image control point n is L n , after spatial triangulation, determine the constant L; Assume that the interval between two measurements of the oblique photogrammetry equipment is Δt; the constant after Δt is L Δt ; At the same time, the reference coordinates remain unchanged, so the x and y coordinates remain unchanged; the elevation changes are: Where z(t+Δt) is the elevation at time t+Δt, and z(t) is the elevation at time t.

10. A method for measuring soil erosion in a waste dump based on a three-dimensional model as claimed in claim 9, characterized in that: Step S8 includes: S81, combining the optimized three-dimensional irregular triangulated network model with the original elevation data to generate a three-dimensional digital model including terrain features; S82, discretize the three-dimensional scene into a regular grid, and for each grid point (i, j), obtain the calculated elevation change Δz(i, j), Δz(i, j)=z(i, j, t+Δt)-z(i, j, t); S83. Analyze the changes of waste slag: If Δz(i, j) = 0, the grid point is considered to remain stable; If Δz(i, j)>0, it is considered that the grid point is piled up; If Δz(i, j) < 0, the grid point is considered to be eroded; S84. Based on the analysis results, a spatial distribution map of the changes in the spoil body is generated, marking the accumulation, erosion and stability areas.

Images