Method for monitoring erosion and vegetation damage in construction zones based on stereoscopic pairs of remote sensing images

Abstract

This invention discloses a method for monitoring erosion and vegetation damage in engineering construction areas based on remote sensing stereo image pairs, belonging to the field of power engineering environmental and water conservation monitoring technology. The method includes: acquiring satellite-coordinated stereo image pairs and multispectral images at several time points before, during, and after construction; performing dense matching on the multi-temporal stereo image pairs to obtain three-dimensional point clouds, and then generating a digital elevation model (DEM) sequence based on the ground point clouds; performing difference operations between the DEMs of each temporal phase and the DEM before construction to calculate the net erosion and erosion modulus of the construction area; using multispectral images to invert the vegetation cover of the construction area at each temporal phase, and using the corresponding DEMs to perform terrain slope correction on the vegetation damage area to obtain the vegetation damage area of ​​the construction area. This invention employs three-dimensional quantitative monitoring technology to achieve accurate calculation of soil erosion and ecological damage in mountainous power transmission and transformation projects.

Description

Technical Field

[0001] This invention belongs to the field of environmental monitoring technology, specifically relating to a method for monitoring erosion and vegetation damage in engineering construction areas based on remote sensing stereo image pairs. Background Technology

[0002] The construction of power transmission and transformation projects, especially in mountainous and hilly areas, can drastically alter the landform and cover due to activities such as tower foundation excavation, construction road building, and temporary soil piling, leading to severe soil erosion and ecological damage. Currently, quantitative assessments of these impacts mainly rely on the following methods, which have significant limitations:

[0003] (1) Traditional engineering surveying method: using total station, GPS-RTK, etc. to measure a single cross section or pile body. This method is extremely inefficient, only measures discrete points, and cannot obtain information on the overall and continuous topographic changes in the construction area. In addition, the fieldwork is large and dangerous, and it is difficult to carry out periodic reviews.

[0004] (2) Single UAV aerial survey method: Periodic aerial surveys are conducted on local areas to generate a DEM and calculate the earthwork volume. This method is limited by flight endurance and airspace, making it difficult to economically and efficiently cover all scattered construction points along a route of tens of kilometers. Issues such as data splicing from different flights and differences in lighting conditions can also introduce errors, affecting the accuracy of time-series comparisons.

[0005] (3) Empirical formula and visual estimation method: Based on the disturbed area and soil type, the general soil loss equation is used for estimation, or the vegetation damage range is visually judged through satellite imagery. The former is difficult to localize parameters and has low accuracy; the latter cannot obtain three-dimensional change information, and the estimation of soil erosion is completely ineffective. In addition, the vegetation damage area is estimated in two-dimensional plane and does not take into account the actual surface area loss caused by topographic relief.

[0006] The aforementioned existing technologies typically separate topographic change monitoring from vegetation monitoring, using different data sources and time phases. This prevents the collaborative analysis of the complete causal chain of "topographic disturbance-soil loss-vegetation response" under the same spatiotemporal reference, resulting in assessment conclusions that lack systematicity and persuasiveness. There is an urgent need to develop a quantitative monitoring technology capable of synchronous, high-frequency, three-dimensional monitoring of the entire, dispersed construction disturbance areas along linear engineering projects. This would enable accurate accounting and scientific management of soil erosion and ecological damage caused by power transmission and transformation projects in mountainous areas. Summary of the Invention

[0007] This invention provides a method for monitoring erosion and vegetation damage in engineering construction areas based on remote sensing stereo image pairs. The three-dimensional quantitative monitoring technology enables accurate calculation of soil erosion and ecological damage caused by power transmission and transformation projects in mountainous areas.

[0008] To achieve the above technical objectives, the present invention adopts the following technical solution:

[0009] A method for monitoring erosion and vegetation damage in engineering construction areas based on remote sensing stereo image pairs includes:

[0010] Acquire satellite-in-orbit stereo image pairs and multispectral images at several time points before, during, and after construction;

[0011] Dense matching is performed on multi-temporal stereo image pairs to obtain three-dimensional point clouds, and then a digital elevation model sequence is generated based on the ground point clouds therein;

[0012] The difference calculation is performed between the digital elevation models of each time phase and the digital elevation model before construction to calculate the net erosion and erosion modulus of the construction area.

[0013] The vegetation coverage of the construction area at different time phases was retrieved using multispectral images, and the vegetation damage area was corrected for topographic slope using the digital elevation model of the corresponding time phase, thus obtaining the vegetation damage area of ​​the construction area.

[0014] Furthermore, before performing dense matching on the stereo image pairs, the acquired stereo image pairs are geometrically refined, specifically:

[0015] Establish a network of connection points between multiple images;

[0016] Regional network adjustment calculations were performed using ground control points to correct the coefficients of the satellite rational function model;

[0017] Using a rational function model with corrected coefficients, stereo correction is performed on the stereo image: the stereo image pairs are projected onto the same virtual imaging plane, so that the epipolar lines of the same object are aligned with the scan lines of the image.

[0018] Furthermore, dense matching of stereo image pairs specifically includes:

[0019] First, the stereo image pairs are radiometrically normalized, and a Gaussian pyramid is constructed.

[0020] Then, a matching cost calculation method based on Census transform is used to calculate the initial matching cost for each pixel in the left image within the disparity search range of its right image;

[0021] Then, a semi-global matching algorithm is applied to optimize the initial cost by minimizing the global energy, resulting in a whole-pixel disparity map. Left-right consistency checks and peak ratio filtering are then used to remove mismatched points. Finally, sub-pixel disparity optimization is achieved through quadratic curve fitting.

[0022] Finally, based on the principle of forward intersection, the three-dimensional point cloud is generated by using the optimized parallax point positions and the rational polynomial function model with corrected coefficients.

[0023] Furthermore, when calculating the initial matching cost using the Census transform-based matching cost calculation method, the template size is adaptively adjusted based on the image gradient. Specifically:

[0024] The horizontal gradient of each pixel in the left image is calculated using the Sobel operator. and vertical gradient Thus, the gradient magnitude is obtained. ;

[0025] Adaptive template determination: if the gradient magnitude of the pixel Select the preset large template; if the pixel gradient magnitude Select the preset medium template; if the pixel gradient magnitude Select a preset small template; among them, For weak texture threshold, Strong texture threshold;

[0026] Census transform coding: For pixels on the left image Using it as the center, iterate through and adaptively determine the remaining pixels within the template window. ; Set the grayscale value of the center pixel of the window grayscale values ​​of neighboring pixels Compare and generate a bit string: if The pixel is encoded as 1 if the encoding is 1, and 0 otherwise; all encoded results are concatenated to obtain the pixel. Census encoding ;

[0027] Matching cost calculation: for left image pixels Parallax of the right image is candidate points The Hamming distance encoded by Census is used as the initial matching cost between the two:

[0028] .

[0029] Furthermore, after generating a 3D point cloud by densely matching the stereo image pairs, the point cloud is further classified and filtered to obtain a ground point cloud, including:

[0030] First, the 3D point cloud generated by dense matching is preprocessed by denoising and regularization sampling.

[0031] Subsequently, a progressive triangulation filtering algorithm is adopted to iteratively separate ground points and non-ground points from the point cloud by dynamically adaptive distance and angle thresholds;

[0032] The method for adjusting the dynamic adaptive distance and angle threshold is as follows:

[0033] (1) Local terrain feature quantification: Before each iteration of encryption, a search neighborhood with radius R is set with the current point to be classified P as the center, and the set S of classified ground points in the neighborhood is counted;

[0034] (2) Dynamic distance threshold calculation: Based on the elevation values ​​of neighboring ground points, calculate the local elevation standard deviation. Therefore, a dynamic distance threshold is set. for:

[0035]

[0036] in, Based on the basic distance threshold, This is the adjustment coefficient for the dynamic distance threshold;

[0037] (3) Dynamic angle threshold calculation: Based on the neighboring ground point S, fit a local plane, calculate the angle between the normal vector of the plane and the horizontal plane, and obtain the local terrain slope. Therefore, a dynamic angle threshold is set. for:

[0038]

[0039] in, Based on the basic angle threshold, This is the adjustment coefficient for the dynamic angle threshold;

[0040] (4) Application of dynamic threshold: Apply the calculated dynamic threshold and Substitute it into the standard PTD judgment process for the classification decision of the current point to be classified, P;

[0041] Finally, the classification results are corrected through manual interactive inspection, and a digital elevation model is generated based on the classified ground point cloud through triangulation linear interpolation.

[0042] Furthermore, the net erosion and erosion modulus of the construction area are calculated, including:

[0043] S31, DEM difference operation:

[0044] ;

[0045] in, express Digital elevation model of the period express Digital elevation model of the period, express period relative to Elevation changes over a period of time;

[0046] S32. Volume Change Statistics: Based on the elevation change within each monitoring unit of the construction area. The volumes of all negative cells and all positive cells are summed separately. The volume of all negative cells is the total excavation volume, and the volume of all positive cells is the total fill volume. The net erosion is obtained by subtracting the total fill volume from the total excavation volume.

[0047] S33, Soil erosion modulus conversion:

[0048] ;

[0049] in, Erosion modulus, Net erosion amount For soil bulk density, The horizontal projected area of ​​the monitoring unit; This refers to the monitoring period.

[0050] Furthermore, the terrain slope of the vegetation damage area is corrected using a digital elevation model. The correction formula is as follows:

[0051]

[0052] in, For the first The horizontal projected area of ​​each damaged pixel For the first The slope of the location of each damaged pixel. This represents the corrected area of ​​vegetation damage.

[0053] Compared with existing technologies, the advantages of this invention are as follows: This invention provides a method for monitoring erosion and vegetation damage in engineering construction areas based on remote sensing stereo image pairs. Through time-series difference analysis, it performs integrated, high-precision quantitative calculation and assessment of changes in earthwork volume, soil erosion, and vegetation damage area caused during power transmission and transformation engineering construction, thereby achieving environmental monitoring. It has the following advantages:

[0054] (1) The three-dimensional volume change and the true surface area are used to replace the two-dimensional plane estimation, which improves the accuracy of the assessment of soil erosion (volume) and vegetation damage area (surface area).

[0055] (2) Physical quantities are directly measured by difference, avoiding empirical model errors. It eliminates the dependence on difficult-to-obtain local parameters and directly outputs verifiable earthwork volume, making the results more reliable.

[0056] (3) Under the same spatiotemporal reference, the changes in topography and vegetation are analyzed simultaneously. This provides a stable and reliable technical guarantee for the field of power transmission and transformation environmental and water conservation monitoring in mountainous areas. Attached Figure Description

[0057] Figure 1 This is a flowchart of the method for monitoring erosion and vegetation damage in engineering construction areas based on remote sensing stereo image pairs, as described in the embodiments of this application. Detailed Implementation

[0058] The embodiments of the present invention will be described in detail below. These embodiments are based on the technical solutions of the present invention and provide detailed implementation methods and specific operation processes to further explain the technical solutions of the present invention.

[0059] Example 1

[0060] This embodiment provides a method for monitoring erosion and vegetation damage in engineering construction areas based on remote sensing stereo image pairs, referencing... Figure 1 As shown, it includes:

[0061] Step 1: At several key time points before, during, and after construction, simultaneously acquire multiple high-resolution satellite co-orbit stereo image pairs and multispectral images covering the target construction area.

[0062] In this embodiment, the stereo image pair and the multispectral image are derived from the same observation on the same satellite platform to ensure the temporal synchronization and geometric consistency of the data. The purpose of this stage is to obtain high-quality, comparable raw data to lay the foundation for subsequent quantitative analysis. The specific embodiment also includes the following steps:

[0063] S11. Delineate Monitoring Units: In GIS software (such as ArcGIS or QGIS), based on the overall project layout plan, accurately delineate each independent construction disturbance area of ​​the power transmission and transformation project, such as tower foundation pits, construction access roads, material storage yards, and spoil disposal sites, using vector polygons. Each polygon should be accompanied by attributes such as "project station number" and "disturbance type" to form a monitoring unit vector file.

[0064] S12. Develop a monitoring plan: Based on the construction organization design, predict key time points such as peak earthwork periods and critical periods of vegetation damage. Plan at least four observation phases: T0 (baseline data before construction), T1 (mid-construction, active earthwork period), T2 (late-construction, significant vegetation damage period), and T3 (completion of construction or post-repair). Acquire satellite imagery data for each period. The core requirement is that each observation must simultaneously acquire "stereo image pairs" and "multispectral images." Stereo image pairs are used for 3D modeling, and multispectral images are used for vegetation analysis.

[0065] S13. Data Acceptance and Quality Inspection: Upon receiving the data, the primary task is to conduct a quality inspection. Check for cloud or snow cover in the imagery, especially in the monitoring unit area. Simultaneously, verify parameters such as "cloud cover percentage" and "side tilt angle" in the metadata. Generally, image data should have a cloud cover of less than 10% and a side tilt angle of less than 15° to ensure minimal imaging geometric distortion and maximum 3D modeling accuracy. Unqualified data requires replacement of the image data or a different time phase.

[0066] Step 2: For multi-temporal stereo image pairs, perform geometric refinement and dense matching to obtain three-dimensional point clouds, and then generate a multi-phase high-precision digital elevation model (DEM) sequence based on the ground point clouds.

[0067] When generating the DEM sequence in this step, it is necessary to filter the digital surface model by integrating concurrent multispectral data to remove the influence of non-surface attachments such as vegetation and temporary buildings, so as to obtain a high-precision DEM that reflects the real terrain.

[0068] Step 2.0, Engineering Data Import: Create a new project in professional photogrammetry processing software (such as ContextCapture, Pix4Dmatic, or Inpho). Import the forward-looking images, backward-looking images, multispectral images, and the satellite's built-in "rational polynomial coefficients" file from the same time phase. Set a unified output coordinate system for all time-phase data.

[0069] Step 2.1, Geometric Refinement: The regional network adjustment method with ground control points is used to perform joint geometric processing on multi-period stereo image pairs.

[0070] Step 2.1.1: Establish a network of connection points among multiple images by combining automatic matching with manual editing. In the photogrammetry software, activate the "Automatic Match Connection Points" function to automatically find a large number (usually hundreds to thousands per image) of image points with the same name as connection points in the overlapping areas of all images.

[0071] Step 2.1.2: Use high-precision ground control points to perform regional network adjustment calculations, correct the coefficients of the satellite rational function model, and complete the absolute positioning.

[0072] Acquire control points: On the reference image (usually the one with the best quality), accurately identify and click on ground feature points with known coordinates as control points.

[0073] For point transfer and adjustment, the software automatically predicts the image-side coordinates of these control points onto all other imagery. Subsequently, a regional network adjustment with control points is performed to solve for the optimal RPC model correction parameters.

[0074] Adjustment result verification and output: Use independent checkpoints that were not involved in the adjustment to evaluate absolute accuracy. View the "Checkpoint Mean Error" in the adjustment report.

[0075] Step 2.1.3 involves verifying the image through independent checkpoints and using a rational function model with corrected coefficients to perform stereo correction on the stereo image. This ensures that the multi-phase image products achieve sub-pixel level geometric alignment accuracy, providing a consistent geometric basis for subsequent temporal difference analysis.

[0076] The revised rational function model enables high-precision cross-calculation between image coordinates (rows and columns) and geographic coordinates (longitude, latitude, and altitude). In practice, the original RPC file can be replaced with the version with the corrected coefficients, and subsequent processing software (such as ENVI, ERDAS, PCI Geomatica, or ArcGIS) will automatically call this new file for coordinate calculations.

[0077] In this embodiment, a modified rational function model is used to project two overlapping images in a stereo pair onto the same virtual imaging plane, ensuring that the epipolar lines of corresponding ground features are strictly aligned with the scan lines of the images. This simplifies stereo matching by requiring only a one-dimensional search along the row direction, significantly improving efficiency and reliability. This geometric refinement process eliminates geometric distortions in the original images caused by pose and terrain, which is a prerequisite for achieving high-precision dense matching.

[0078] Step 2.2: Dense matching generates point clouds.

[0079] Step 2.2.1: Perform radiation normalization on the geometrically refined stereo image pairs and construct a Gaussian pyramid;

[0080] A Gaussian pyramid is a multi-scale image structure: by repeatedly applying Gaussian smoothing (blurring) and interleaved downsampling to the image, a series of image layers are generated, ranging from large to small and from clear to blurry. In the stereo matching of this invention, its core function is to achieve fast and robust matching "from coarse to fine": first, the approximate disparity is estimated on the top layer (smallest and most blurry) image, and then the resolution is refined layer by layer to the original resolution of the bottom layer, thereby significantly reducing the computational load and effectively handling large variations in ground feature disparity and suppressing noise interference. It is a key supporting technology for algorithms such as semi-global matching.

[0081] Step 2.2.2: The matching cost calculation method based on Census transform is used to calculate the initial matching cost for each pixel in the left image within the disparity search range of its right image;

[0082] The matching cost calculation method based on Census transform adaptively adjusts the template size based on image gradient when calculating the initial matching cost. Specifically:

[0083] The horizontal gradient of each pixel in the left image is calculated using the Sobel operator. and vertical gradient Thus, the gradient magnitude is obtained. ;

[0084] Adaptive template determination: if the gradient magnitude of the pixel Select a preset large template (e.g., 11×1111×11) to increase the uniqueness of the matching information; if the gradient magnitude of the pixel Select a preset medium template (e.g., 9×99×9); if the pixel gradient magnitude Select a preset small template (e.g., 7×77×7) to preserve edge details and prevent over-smoothing; among them, For weak texture threshold, Strong texture threshold;

[0085] Census transform coding: For pixels on the left image Using it as the center, iterate through every other pixel within the adaptive template window. ; Set the grayscale value of the center pixel of the window grayscale values ​​of neighboring pixels Compare and generate a bit string: if The pixel is encoded as 1 if the encoding is 1, and 0 otherwise; all encoded results are concatenated to obtain the pixel. Census encoding ;

[0086] Matching cost calculation: for left image pixels Parallax of the right image is candidate points The matching cost of its Census transform is the Hamming distance between the two Census codes, which is the number of corresponding bit values ​​that are different:

[0087] .

[0088] Step 2.2.3: Apply the semi-global matching (SGM) algorithm to optimize the global energy minimization of the initial cost to obtain an integer-pixel disparity map. Then, supplement it with left-right consistency checks and peak ratio filtering to remove mismatched points. Finally, sub-pixel-level disparity optimization is achieved through quadratic curve fitting.

[0089] In the SGM algorithm, the smoothness constraint is implemented through a penalty parameter in the energy function. This energy function defines the optimization objective of finding the optimal disparity map.

[0090]

[0091] The first term in the energy function is the data term: It is a pixel At parallax The matching cost at that time. It measures the points on the left image. How well does it match the corresponding point on the right image? The lower the cost, the better the match.

[0092] The second term in the energy function is the smoothing term: this part introduces the neighboring pixels. The constraint between them. Its logic is: the surface of objects in the real world is usually continuous, so the disparity (depth) of adjacent pixels should also change smoothly and should not jump arbitrarily. It is a penalty applied when the parallax change between adjacent pixels is 1 pixel. It is used to handle sloping or curved surfaces, allowing for small, continuous slope changes in parallax. It's a larger penalty applied when the disparity change between adjacent pixels is greater than one pixel. It's primarily used to handle object edges, where the depth truly changes abruptly.

[0093] Step 2.2.4: Based on the principle of forward intersection, the three-dimensional point cloud is generated by using the optimized disparity point position and the rational polynomial function model with coefficient correction.

[0094] The optimized disparity points refer to the sub-pixel level disparity map obtained in the previous step after left-right consistency checks, peak ratio filtering to eliminate mismatches, and quadratic curve fitting.

[0095] In particular, considering the characteristics of large terrain undulations and weak texture areas in the construction area, the smoothing constraint parameters in the SGM algorithm are adaptively adjusted, and multi-view geometric constraints can be introduced during the matching process. For example, left-right consistency checks can be used to eliminate suspicious matches to ensure the integrity and accuracy of the point cloud in complex construction scenarios.

[0096] Step 2.3, point cloud classification and filtering.

[0097] Step 2.3.1: Denoise and regular sampling preprocessing are performed on the 3D point cloud generated by dense matching;

[0098] Step 2.3.2: The progressive triangulation filtering algorithm is adopted to iteratively separate ground points and non-ground points from the point cloud by dynamically adaptive distance and angle thresholds.

[0099] The method for adjusting the dynamic adaptive distance and angle threshold is as follows:

[0100] (1) Local terrain feature quantification: Before each iteration of encryption, a search neighborhood with radius R is set with the current point to be classified P as the center, and the set S of classified ground points in the neighborhood is counted;

[0101] (2) Dynamic distance threshold calculation: Based on the elevation values ​​of neighboring ground points, calculate the local elevation standard deviation. Therefore, a dynamic distance threshold is set. for:

[0102]

[0103] in, Use the base distance threshold (e.g., 0.3 meters). This is an adjustment coefficient for the dynamic distance threshold (e.g., 0.5-1.5); when the local terrain has large undulations (i.e., elevation standard deviation). When the distance threshold is large, it is automatically widened to allow the preservation of real terrain features such as steep slopes and embankments.

[0104] (3) Dynamic angle threshold calculation: Based on the neighboring ground point S, fit a local plane, calculate the angle between the normal vector of the plane and the horizontal plane, and obtain the local terrain slope. Therefore, a dynamic angle threshold is set. for:

[0105]

[0106] in, The basic angle threshold (e.g., 5 degrees). This is an adjustment coefficient for the dynamic angle threshold (e.g., 0.2-0.8); when the local slope is large, the angle threshold is increased accordingly to adapt to changes in terrain.

[0107] (4) Application of dynamic threshold: Apply the calculated dynamic threshold and Substituting into the judgment process of the standard progressive triangulation encryption filtering algorithm, the classification decision for the current point P to be classified is: whether to add it to the ground point set.

[0108] The standard procedure for the standard progressive triangulation network encryption filtering algorithm includes: selecting seed points → constructing an initial TIN → calculating the distance and angle from the point to be determined to the TIN plane → comparing with a preset threshold → iterative encryption.

[0109] In step 2.3.2, during the separation of ground points and non-ground points, vegetation index information from concurrent multispectral images is fused to assist in the identification and removal of vegetation points; and for characteristic terrain features such as steep slopes and foundation pits formed during construction, a local optimization algorithm based on cross-section analysis is used to process them in order to preserve their true geometric shape.

[0110] The method of using vegetation index information from concurrent multispectral imagery to assist in vegetation point identification essentially involves fusing the spectral information of two-dimensional multispectral imagery with the geometric information of three-dimensional point clouds. Vegetation exhibits high reflectivity in the near-infrared band and strong absorption in the red band; therefore, indices such as NDVI can significantly distinguish between vegetation and non-vegetation. The specific implementation plan is as follows:

[0111] Step a1: Data registration and index creation.

[0112] 1. Spatial registration: Ensure that the point cloud data and the concurrent multispectral image have a completely consistent coordinate system and spatial range. This step is usually completed in the "geometric refinement" stage of step S2.

[0113] 2. Inverse Distance Weighted Mapping: For each 3D point (X,Y,Z) in the point cloud, based on its planar coordinates (X,Y), the corresponding pixel position is found on the multispectral image, and the multispectral band values ​​(such as red, green, blue, and near-infrared) of that pixel are assigned as attributes to the point. This is equivalent to adding "spectral information" to each point cloud.

[0114] Step a2: Vegetation Index Calculation. Multiple vegetation indices are calculated for each point cloud to enhance the distinction between vegetation and non-vegetation areas. Commonly used indices include: Normalized Difference Vegetation Index (NDVI): NDVI = (NIR - Red) / (NIR + Red); Excess Greenness Index (ExG): ExG = 2G - RB (suitable for point clouds with only RGB information); Green Leaf Index (LI): LI = (2G - RB) / (2G + R + B).

[0115] Step a3: Threshold segmentation or machine learning classification. Based on the distribution characteristics of vegetation indices, perform preliminary identification of vegetation points: set an empirical threshold (such as NDVI>0.3 or ExG>20), and mark points that meet the conditions as "suspected vegetation points".

[0116] Step a4: Integrate with geometric classification results.

[0117] If a point is labeled as a "ground point" by the PTD algorithm, but spectral indices (such as NDVI) indicate that it is vegetation, then the spectral information should be prioritized, and the point should be reclassified as a "vegetated point" and removed from the ground point set. The final output ground point cloud should consist of pure ground points verified by spectral analysis.

[0118] In addition, for the steep slopes, foundation pits, and other characteristic terrain features formed during construction, a local optimization algorithm based on cross-section analysis is used for processing, specifically including:

[0119] (1) Cross-section layout: Based on the local orientation and complexity of the terrain, a series of cross-section lines perpendicular to the terrain features are adaptively laid out, and a cross-section width threshold is set to extract cross-section point clouds.

[0120] (2) Feature recognition: Project the cross-sectional point cloud onto a two-dimensional plane, and identify abrupt changes in terrain by calculating indicators such as the slope change rate, elevation difference and curvature of adjacent points, and use a multi-section continuity verification mechanism to eliminate false changes;

[0121] (3) Local optimization: The identified terrain feature points are classified and processed, the point cloud at steep slopes, pit edges, etc. is densified and restored, and adaptive smoothing constraints are applied to the regions between feature points to ensure the overall continuity of the cross-section line while preserving the terrain change features.

[0122] (4) 3D fusion: The optimized point clouds of each section are fused back into the 3D space. Spatial consistency is checked by constructing a local triangulation network to ensure the continuity of feature points (such as steep slope lines) in the 3D space. Finally, the ground point cloud with feature fidelity is output.

[0123] Step 2.3.3 involves manually checking and correcting the classification results, and generating a digital elevation model based on the classified ground point cloud using triangular mesh linear interpolation.

[0124] Step 3 involves performing a difference calculation between the digital elevation models of the construction and restoration periods and the digital elevation model before construction. Based on the difference results, the net erosion and erosion modulus of the construction area are calculated. This step directly quantifies the earthwork migration and soil loss caused by construction by calculating the topographic changes.

[0125] Step 3.1, DEM difference calculation:

[0126] ;

[0127] in, express Digital elevation model of the period, express Digital elevation model of the period, express period relative to The change in elevation over a period of time.

[0128] Step 3.2, Volume Change Statistics: First, a change threshold is set (e.g., ±0.15 meters) to eliminate measurement noise from the DEM itself. Then, using the area statistics function, the elevation change within the vector polygon range of each monitoring unit is analyzed. Statistical analysis was conducted. The volumes of all negative value pixels (total excavation volume) and all positive value pixels (total fill volume) were summarized separately. Net erosion = total excavation volume - total fill volume, which directly reflects the net volume of soil loss.

[0129] Step 3.3, Soil erosion modulus conversion:

[0130] ;

[0131] in, Erosion modulus; Net erosion amount, in cubic meters; Soil bulk density can be obtained through on-site sampling experiments; To find the horizontal projected area of ​​the monitoring unit, simply calculate the area of ​​the vector polygon of the monitoring unit in GIS, and convert the unit to square kilometers. The monitoring period is calculated in years. The calculated erosion modulus can be compared with the "Classification and Grading Standard for Soil Erosion" to objectively evaluate the erosion intensity.

[0132] Step 4: Use multispectral images to invert the vegetation coverage of the construction area at each time phase, and use the digital elevation model of the corresponding time phase to correct the terrain slope, so as to obtain the vegetation damage area of ​​the construction area.

[0133] Step 4.1, Calculate vegetation cover: Calculate NDVI for each period of multispectral imagery, and invert vegetation cover FVC using a pixel-dichotomy model.

[0134] Step 4.2, detect the extent of damage and calculate. ; express Vegetation cover during the period express Vegetation cover during the period express period relative to The change in vegetation cover over a period of time; then, based on the threshold, it is determined whether the vegetation in the construction area has been damaged.

[0135] Step 4.3: Use the digital elevation model to correct the terrain slope of the vegetation damage area. The correction calculation formula is as follows:

[0136]

[0137] in, The first in the vegetation damage area The horizontal projected area of ​​each pixel The first in the vegetation damage area The slope of each pixel's location is calculated using a digital elevation model; The corrected area of ​​vegetation damage is the vegetation damage area.

[0138] Step 5: Spatial overlay and coupling analysis of the erosion / deposition distribution map and the vegetation damage distribution map, including but not limited to: statistically analyzing the proportion of vegetation damage area in severely eroded areas, analyzing the spatial correlation between vegetation damage areas and erosion / deposition areas, and generating a comprehensive quantitative assessment report.

[0139] This step involves integrating three-dimensional terrain change and two-dimensional vegetation change information, performing correlation analysis, and generating a comprehensive quantitative assessment report. Step 5 includes the following core steps:

[0140] Step 5.1: Convert the erosion intensity data obtained in Step 3 and the vegetation damage area obtained in Step 4 into vector surface layers.

[0141] Step 4, in the precise calculation of vegetation damage area based on terrain correction in GIS, actually involves statistically analyzing damaged pixels on the output vegetation damage data (raster data containing the spatial distribution of damaged areas, commonly referred to as a vegetation damage binary mask), rather than simply calculating the numerical value of the damaged area. Therefore, the vegetation damage binary mask raster generated in Step 4 can be converted into a vegetation damage vector surface layer using a raster vectorization tool, and attribute fields such as area and damage level can be added to each patch. Simultaneously, the erosion intensity raster generated in Step 3 (which has been reclassified according to national grading standards) can also be converted into an erosion intensity vector surface layer, which, together with the vegetation damage vector surface layer, serves as the basis for the spatial coupling analysis in Step 5.1.

[0142] S52. Spatial Coupling Analysis: In GIS, the "Erosion Zone Vector Layer" and the "Vegetation Damage Zone Vector Layer" are spatially overlaid (Intersect). The resulting layer will contain both the erosion level and damage level attributes of each polygon.

[0143] Example 2

[0144] This embodiment, based on Embodiment 1, uses the typical N87 tower foundation construction area (including the connecting access road) in a 500kV transmission line project in a mountainous area as the monitoring object. The aim is to fully demonstrate the method of this invention and achieve high-precision, quantitative environmental monitoring of the construction site throughout the entire construction cycle, from pre-construction to construction. This area is hilly with an average slope of approximately 15°, and some excavated slopes exceeding 30°. Traditional monitoring methods struggle to accurately assess the earthwork volume and ecological damage.

[0145] In step S11, when defining the monitoring unit, in the GIS software, a 100m×100m rectangular area is defined as the core monitoring unit with the center of the N87 tower base as the origin, and then extended outward to include a 200m construction access road, with a total monitoring area of ​​approximately 3ha.

[0146] In step S12, when developing the monitoring plan, three key time phases are determined based on the construction plan:

[0147] T0 (Local Period): May 10, 2025, construction machinery has not yet entered the site.

[0148] T1 (Peak Construction Period): On August 20, 2025, the foundation pit of the tower was excavated and the construction access road was completed, exposing a large amount of bare soil.

[0149] T2 (Initial Recovery Period): On November 15, 2025, the tower base was completed, and some areas were leveled.

[0150] In step S13, during data reception and quality control, satellite imagery data for each time period (T0, T1, T2) is acquired. The requirement is that "stereo image pairs" and "multispectral images" must be acquired simultaneously for each observation. Upon receiving the data, the primary task is quality control. Unqualified data must be replaced with new imagery or have a different timeframe. The following requirements apply to the data:

[0151]

[0152] Furthermore, in step 2.0, during the multi-phase data import and geometric processing, the professional photogrammetry software Pix4Dmatic was used. Six Gaofen-7 stereo image pairs (forward and backward views) from three time phases were imported together to create a new project. Four ground control points (GCPs) with known coordinates were selected at stable road intersections around the survey area, and these points were precisely marked on the imagery.

[0153] The report after adjustment in step 2.1 shows that the average reprojection error of the connection points is 0.42 pixels and the root mean square error of the control points is 0.31 meters, indicating that the geometric model has excellent accuracy.

[0154] In this embodiment, the semi-global matching (SGM) algorithm is selected in the Pix4Dmatic software for dense matching, generating dense point clouds for three time phases: T0, T1, and T2, with an average point density of approximately 16 points / m. 2 .

[0155] In this embodiment, when using the progressive triangulation filtering algorithm to separate ground points from non-ground points, the maximum building size parameter is set to 15m and the maximum terrain slope parameter is set to 70°, taking into account the characteristics of the construction area, so as to effectively preserve the steep slope characteristics formed by excavation.

[0156] In the point cloud classification process, this embodiment integrates vegetation index information from contemporaneous multispectral imagery to assist in the identification and removal of vegetation points: the NDVI values ​​calculated from the contemporaneous multispectral imagery are assigned to the point cloud to assist in the identification and removal of vegetation points (NDVI>0.3), ensuring that the generated DEM is a bare soil DEM. For characteristic terrain features such as steep slopes and foundation pits formed during construction, a local optimization algorithm based on cross-section analysis is used to process them, preserving their true geometric shape. Finally, the classification results are corrected through manual interactive inspection, and based on the classified ground point cloud, triangulation is performed on the classified ground points to generate a high-precision DEM with a grid spacing of 0.5 meters, named as follows: , , After verification at the checkpoints, the vertical accuracy error of the DEM was found to be 0.25 meters.

[0157] When performing DEM differential calculations, raster calculation tools are used in remote sensing or GIS software (such as ENVI or ArcGIS) to calculate elevation changes. In this example, the following operations are performed using the "Raster Calculator" in ArcGIS Pro:

[0158]

[0159]

[0160] In this example, the threshold for volume change statistics is set to ±0.15 meters (to exclude measurement noise). Areas with a median error less than -0.15m are classified as "cut areas," and areas with a median error greater than 0.15m are classified as "fill areas." (After rigorous adjustment and filtering, the actual effective median error estimate of the DEM is 0.08m, so twice the median error ≈ 0.16m. Taking 0.15m is a reasonable, rigorous, and conservative threshold.)

[0161] In this embodiment, statistics are compiled for the vector range of the N87 monitoring unit. The results are as follows: Total excavation volume =1850m 3 Total fill volume =1200m 3 Net erosion =650m 3 (The gap in the earthwork indicates that some earthwork has been removed or eroded.)

[0162] To align the results with soil and water conservation standards, the soil erosion modulus needs to be calculated. In this embodiment, the soil bulk density of the area was measured by on-site sampling: Monitoring period: =T1-T0=102 days≈0.28 years, Horizontal projected area of ​​monitoring unit: A=30000m² 2 =0.03km 2 Calculate the soil erosion modulus: = (650) 1.5) / (0.03) 0.28) = 116071t / (km) 2 ·a).

[0163] According to the "Classification and Grading Standards for Soil Erosion", the soil erosion intensity at this construction site has reached the "extremely severe erosion" level in the short term.

[0164] When extracting the extent of vegetation damage, the Normalized Difference Vegetation Index (NDVI) for time phases T0 and T1 is calculated, and the Free Variable Vegetation (FVC) is retrieved. This step also includes the following operations:

[0165] Step 1: Data preparation and preprocessing: including binary masking of vegetation damage, high-precision digital elevation model, and data preprocessing check.

[0166] Binary Vegetation Damage Mask: This data was obtained during the multispectral change detection process. The data is in raster format, with a pixel value of 1 representing a "damaged vegetation area" and a pixel value of 0 representing a "non-damaged area." This mask defines the spatial range within which the area to be calculated needs to be determined.

[0167] High-precision digital elevation model: A filtered DEM generated from stereo image pairs. Its spatial resolution (e.g., 0.5m) should match or exceed that of the vegetation damage mask, and the coordinate system must be completely consistent.

[0168] Data preprocessing check: In GIS software, the vegetation damage mask is overlaid with the DEM to ensure complete spatial registration. The spatial resolution of the DEM is used to calculate the planar projected area A_plane of a single pixel (e.g., for a 0.5m resolution DEM). = 0.25m 2 ).

[0169] Step 2: Calculate the Normalized Difference Vegetation Index (NDVI), including band determination, NDVI formula application, and quality check.

[0170] Band determination: The near-infrared band and red band are determined based on the image. In this example, the near-infrared band of the satellite image is B8, and the red band is B4.

[0171] NDVI Formula Application: Using the raster calculator, perform calculations pixel by pixel:

[0172]

[0173] This step outputs two floating-point grids: Theoretically, the range of pixel values ​​is [-1, 1].

[0174] Quality Check: Check the NDVI image. Densely vegetated areas should show high values ​​(0.6~0.8), bare soil / buildings around 0.2, and water bodies negative. If the results are abnormal, check if the input bands are correct. Statistical Value Check: A reasonable NDVI image should have a mean between 0.2 and 0.7. A large number of -1 or 1 values ​​may indicate a calculation error.

[0175] Step 3: Invert vegetation cover (FVC). The core of FVC inversion is the pixel-based binary model. Its principle is that the spectral information of a pixel is a linear mixture of its vegetation and bare soil components. The key to this operation is determining the NDVI values ​​of pure vegetation and pure bare soil. The model formula is:

[0176]

[0177] Where NDVI is the NDVI value of the current pixel. The NDVI value of a pure bare soil pixel. This represents the NDVI value of a pure vegetation pixel.

[0178] To be determined and Substitute the values ​​into the pixel binary model formula and calculate the FVC for each NDVI image period in the raster calculator. Confirm. and The method used in this implementation example is the empirical statistical method. The specific operation steps are as follows: 1. In 1. On the image, select areas of pure bare soil with high confidence (such as newly excavated surfaces or dried riverbeds). 2. Calculate the NDVI values ​​of all pixels within this area, and take the 5th percentile of the cumulative frequency as... 3. Select pure vegetation areas with high confidence levels (such as dense woodlands), calculate their NDVI values, and take the 95th percentile of the cumulative frequency as the baseline. Due to noise and model errors, the calculation results may contain outliers less than 0 or greater than 1. These outliers need to be forced to the valid range of [0, 1] using a conditional function (FVC<0 is set to 0, FVC>1 is set to 1). and Two grids, with cell values ​​representing the percentage of vegetation cover (0%~100%) at that location.

[0179] Step 4: Change detection and damage mask generation, including calculating the amount of change, setting the damage threshold, and generating a binary damage mask.

[0180] Calculate the change in vegetation cover:

[0181] Set a damage threshold: Determine a A negative change threshold is used to define "damage". This threshold needs to be set based on the actual field conditions. An empirical threshold is typically used. A coverage decrease of less than -0.3 (i.e., a decrease in coverage exceeding 30 percentage points) is considered significant damage. Statistical determination: In undisturbed, stable surrounding areas, the following parameters are calculated... The standard deviation is set, and the threshold is set to -2 standard deviations to exclude natural fluctuations.

[0182] Generate a binary damage mask: based on conditional judgment ( <Damage threshold), pixels that meet the condition are assigned a value of 1 (damaged area), otherwise 0 (undamaged area). This raster is the direct input for the next step of terrain correction area calculation.

[0183] In this embodiment, pixels with an FVC decrease exceeding 50% are extracted through change detection, and a binary mask of vegetation damage is generated. The planar projection shows the damaged area to be 5,200 m². 2 .

[0184] Step 4.3 calculates the actual terrain damage area. The purpose is to eliminate the influence of terrain undulations on the calculation of vegetation damage area, correcting the two-dimensional estimate based on planar imagery to a three-dimensional surface area that reflects the actual land surface conditions. The specific operation procedure is as follows:

[0185] Step 1: Data standardization and spatial alignment. This step is a prerequisite for ensuring that all subsequent calculations are performed on a unified geometric basis.

[0186] Data Preparation: Binary Vegetation Damage Mask Data: This data was obtained from previous vegetation change detection, with pixel values ​​of only 0 (representing undamaged areas) or 1 (representing damaged vegetation areas). High-Precision Digital Elevation Model: DEM data of the same region and phase as the vegetation damage mask.

[0187] Perform spatial consistency checks and processing: Load both raster layers simultaneously in the software. Through overlay display and comparison, ensure complete consistency in spatial coordinate system, geographic extent, cell size, and cell alignment. (If discrepancies exist, the "resample" or "register" tool must be used, with the DEM as the spatial reference, to resample and align the vegetation damage mask, ensuring each cell corresponds one-to-one.)

[0188] Step 2: DEM-based slope extraction. This step aims to obtain the terrain slope angle at each pixel location, providing geometric parameters for area correction.

[0189] To use the slope tool: In the software's "Spatial Analysis" or "Terrain Analysis" toolbox, select the "Slope" calculation tool.

[0190] Set the calculation parameters:

[0191] Input raster: Specify the high-precision DEM prepared in the previous step.

[0192] Output raster: Set the path and name of the output slope file (e.g., Slope_Degree).

[0193] Output unit: "degrees" must be selected as the unit of measurement for the slope angle. This is a fundamental requirement for substituting the slope into subsequent trigonometric function calculations.

[0194] Generate slope raster: Perform the calculation, and the software will output a new raster file. Each cell value in the file (denoted as θ) represents the tilt angle of the corresponding land surface location, with 0° for the horizontal plane and 90° for the vertical plane.

[0195] Step 3: Calculate the surface area correction coefficient. The physical principle behind this step is that the actual surface area represented by a pixel on a slope is its projected area on the ground plane. times.

[0196] Apply the correction formula: Use the "Raster Calculator" tool in your GIS software and enter the calculation formula: In calculation Beforehand, the slope value in degrees needs to be converted to radians. The complete expression is usually as follows: .

[0197] Addressing slope overflow: In the above calculation steps, when the slope θ approaches 90 degrees, cos(θ) approaches 0, causing the correction coefficient k to approach infinity, resulting in calculation overflow or distorted results. Therefore, a maximum effective slope threshold (such as 60 or 70 degrees) is set before calculation. A conditional function is used to ensure that pixels with slopes greater than this threshold are calculated according to this threshold. The formula is revised as follows:

[0198] ), ( ))

[0199] This step ensures stable calculations and conforms to engineering realities (vegetation is typically sparse on extremely steep slopes). In this embodiment, the average slope of the area is approximately 15°, and the average correction factor k≈1.04.

[0200] Step 4: Pixel-level true area calculation and accumulation. This step integrates damage location information, pixel size, and terrain coefficient to complete the accurate statistics of the true area.

[0201] 1. Pixel-level true area calculation: Continue using the "Raster Calculator" tool in the GIS software. The calculation formula is: True surface area A_real = Vegetation damage mask Correction coefficient k Pixel plane area A_plane. This operation is performed pixel-by-pixel: only pixels with a damage mask value of 1 are included in the calculation, and the result is the actual ground surface area corresponding to that pixel.

[0202] 2. Area Summarization: The raster data output from the above steps is statistically summed. Using the "Regional Statistics" or "Raster Attribute Statistics" tool, the sum of the A_real values ​​of all damaged pixels is obtained directly. This sum is the "true surface area of ​​vegetation damage after terrain correction." In this embodiment, after terrain correction, the true surface area of ​​vegetation damage is 5408 m². 2 .

[0203] Step 5: Result Verification and Comparative Analysis

[0204] 1. Calculate the planar projected area: As a reference, count the total number of pixels with a value of 1 in the vegetation damage mask, multiply by A_plane, and obtain the planar projected area using the traditional method.

[0205] 2. Quantifying the impact of terrain: Calculated area increase percentage = (Actual surface area - Planar projected area) / Planar projected area 100%.

[0206] In this embodiment, the traditional planar algorithm underestimated the vegetation damage area by approximately 4%. In slope areas with steeper gradients, the underestimation rate can reach over 15%.

[0207] Furthermore, step 5 includes the following steps:

[0208] Step 1: Spatial Coupling Analysis: This involves analyzing the "severely eroded zone" (| A spatial overlay analysis was performed on the "1-meter-wide" and "severely damaged vegetation zone (FVC decrease > 50%)" areas. In this embodiment, the spatial overlap between the two was as high as 88%, strongly demonstrating that the excavation and backfilling activities of construction machinery were the direct cause of vegetation destruction.

[0209] Step 2: Risk point identification: In this embodiment, the analysis found that there was a continuous erosion signal on the downstream slope of a temporary soil dump (fill area) and the vegetation was damaged. The system marked it as a "high risk point of soil erosion".

[0210] Step 3: Relevant data generation: In this embodiment, the following conclusions were finally obtained: "The N87 work site generated a net erosion of 650m³, with an erosion intensity of extremely strong; it caused actual vegetation damage of 5408m², mainly attributed to direct construction activities; one high-risk point for soil erosion was found, and it is recommended to reinforce the slope toe and set up an intercepting drainage ditch."

[0211] In this embodiment, based on the above findings, on-site management personnel accurately located the high-risk areas according to the report and took appropriate measures. During monitoring in phase T2, the erosion signal in this area weakened, achieving a closed-loop monitoring system.

[0212] This invention provides a quantitative assessment method for erosion and vegetation damage in power transmission and transformation engineering construction areas based on multi-phase remote sensing stereo image pairs. Compared with traditional methods, it offers higher accuracy and efficiency, improving earthwork volume accuracy to the cubic meter level. Vegetation area calculations correct for terrain errors, providing irrefutable evidence for ecological compensation. The entire process is based on remote sensing measurements and physical formulas, eliminating the subjectivity of human estimation. Through spatial coupling analysis, it reveals the inherent causal relationship between construction activities and environmental impacts, shifting management from post-event punishment to precise in-process prevention and control. It provides a feasible technical solution for the normalized and comprehensive environmental and water conservation supervision of long-distance linear engineering projects.

[0213] The above embodiments are preferred embodiments of this application. Those skilled in the art can make various changes or improvements based on them. Without departing from the overall concept of this application, these changes or improvements should fall within the scope of protection claimed in this application.

Claims

1. A method for monitoring erosion and vegetation damage in engineering construction areas based on remote sensing stereo image pairs, characterized in that, include: Acquire satellite-in-orbit stereo image pairs and multispectral images at several time points before, during, and after construction; Dense matching is performed on multi-temporal stereo image pairs to obtain three-dimensional point clouds, and then a digital elevation model sequence is generated based on the ground point clouds therein; The difference calculation is performed between the digital elevation models of each time phase and the digital elevation model before construction to calculate the net erosion and erosion modulus of the construction area. The vegetation coverage of the construction area at different time phases was retrieved using multispectral images, and the vegetation damage area was corrected for topographic slope using the digital elevation model of the corresponding time phase, thus obtaining the vegetation damage area of ​​the construction area.

2. The method for monitoring erosion and vegetation damage in engineering construction areas based on remote sensing stereo image pairs according to claim 1, characterized in that, Before performing dense matching on the stereo image pairs, the process also includes geometric refinement of the acquired stereo image pairs, specifically: Establish a network of connection points between multiple images; Regional network adjustment calculations were performed using ground control points to correct the coefficients of the satellite rational function model; Using a rational function model with corrected coefficients, stereo correction is performed on the stereo image: the stereo image pairs are projected onto the same virtual imaging plane, so that the epipolar lines of the same object are aligned with the scan lines of the image.

3. The method for monitoring erosion and vegetation damage in engineering construction areas based on remote sensing stereo image pairs according to claim 1, characterized in that, Dense matching of stereo image pairs specifically includes: First, the stereo image pairs are radiometrically normalized, and a Gaussian pyramid is constructed. Then, a matching cost calculation method based on Census transform is used to calculate the initial matching cost for each pixel in the left image within the disparity search range of its right image; Then, a semi-global matching algorithm is applied to optimize the initial cost by minimizing the global energy, resulting in a whole-pixel disparity map. Left-right consistency checks and peak ratio filtering are then used to remove mismatched points. Finally, sub-pixel disparity optimization is achieved through quadratic curve fitting. Finally, based on the principle of forward intersection, the three-dimensional point cloud is generated by using the optimized parallax point positions and the rational polynomial function model with corrected coefficients.

4. The method for monitoring erosion and vegetation damage in engineering construction areas based on remote sensing stereo image pairs according to claim 3, characterized in that, The matching cost calculation method based on Census transform adaptively adjusts the template size based on image gradient when calculating the initial matching cost. Specifically: The horizontal gradient of each pixel in the left image is calculated using the Sobel operator. and vertical gradient Thus, the gradient magnitude is obtained. ; Adaptive template determination: if the gradient magnitude of the pixel Select the preset large template; If the gradient magnitude of the pixel Select the preset medium template; If the gradient magnitude of the pixel Select a preset small template; among them, For weak texture threshold, Strong texture threshold; Census transform coding: For pixels on the left image Using it as the center, iterate through and adaptively determine the remaining pixels within the template window. ; Set the grayscale value of the center pixel of the window grayscale values ​​of neighboring pixels Compare and generate a bit string: if The pixel is encoded as 1 if the encoding is 1, and 0 otherwise; all encoding results are concatenated to obtain the pixel. Census encoding ; Matching cost calculation: for left image pixels Parallax of the right image is candidate points The Census-encoded Hamming distance is used as the initial matching cost between the two: 。 5. The method for monitoring erosion and vegetation damage in engineering construction areas based on remote sensing stereo image pairs according to claim 1, characterized in that, After generating a 3D point cloud by densely matching stereo image pairs, the point cloud is further classified and filtered to obtain a ground point cloud, including: First, the 3D point cloud generated by dense matching is preprocessed by denoising and regularization sampling. Subsequently, a progressive triangulation filtering algorithm is adopted to iteratively separate ground points and non-ground points from the point cloud by dynamically adaptive distance and angle thresholds; The method for adjusting the dynamic adaptive distance and angle threshold is as follows: (1) Local terrain feature quantification: Before each iteration of encryption, a search neighborhood with radius R is set with the current point to be classified P as the center, and the set S of classified ground points in the neighborhood is counted; (2) Dynamic distance threshold calculation: Based on the elevation values ​​of neighboring ground points, calculate the local elevation standard deviation. Therefore, a dynamic distance threshold is set. for: ； in, Based on the basic distance threshold, This is the adjustment coefficient for the dynamic distance threshold; (3) Dynamic angle threshold calculation: Based on the neighboring ground point S, fit a local plane, calculate the angle between the normal vector of the plane and the horizontal plane, and obtain the local terrain slope. Therefore, a dynamic angle threshold is set. for: ； in, Based on the basic angle threshold, This is the adjustment coefficient for the dynamic angle threshold; (4) Application of dynamic threshold: The calculated dynamic threshold is applied. and Substitute it into the standard PTD judgment process for the classification decision of the current point to be classified, P; Finally, the classification results are corrected through manual interactive inspection, and a digital elevation model is generated based on the classified ground point cloud through triangulation linear interpolation.

6. The method for monitoring erosion and vegetation damage in engineering construction areas based on remote sensing stereo image pairs according to claim 1, characterized in that, Calculate the net erosion and erosion modulus of the construction area, including: S31, DEM difference operation: ； in, express Digital elevation model of the period express Digital elevation model of the period express period relative to Elevation changes over a period of time; S32. Volume Change Statistics: Based on the elevation change within each monitoring unit of the construction area. The volumes of all negative cells and all positive cells are summed separately. The volume of all negative cells is the total excavation volume, and the volume of all positive cells is the total fill volume. The net erosion is obtained by subtracting the total fill volume from the total excavation volume. S33, Soil erosion modulus conversion: ； in, Erosion modulus, Net erosion amount For soil bulk density, The horizontal projected area of ​​the monitoring unit; This refers to the monitoring period.

7. The method for monitoring erosion and vegetation damage in engineering construction areas based on remote sensing stereo image pairs according to claim 1, characterized in that, The terrain slope is corrected for the vegetation damage area using a digital elevation model. The correction calculation formula is as follows: ； in, For the first The horizontal projected area of ​​each damaged pixel For the first The slope of the location of each damaged pixel. This represents the corrected area of ​​vegetation damage.

Images