Method for assessing the impact of railway tunnel works on the vegetation on the roof of the tunnel based on comprehensive data analysis

Abstract

The application discloses a railway tunnel engineering influence on vegetation on a hole top based on comprehensive data analysis, and concretely relates to the field of geographic information system and ecological remote sensing monitoring technology, and comprises the following steps: constructing a slope direction-slope gradient coupling surface orientation factor and a buffer elevation gradient offset index, introducing the normalized index into a dynamic buffer distance function to correct the topographic response of a traditional two-dimensional buffer zone; then, constructing a non-equidistant and direction-sensitive slope surface adaptive buffer zone based on digital elevation model data, and outputting a dynamic buffer layer in a vector form; superimposing multi-temporal remote sensing images, extracting vegetation index variation trends in different buffer level regions, and combining a benchmark control area differential correction and a spatial interpolation method to output a vegetation influence heat map and an error analysis map.

Description

Technical Field

[0001] This invention relates to the field of power system operation monitoring and risk prevention and control technology, specifically to an assessment method for the impact of railway tunnel engineering on the vegetation on the tunnel roof based on comprehensive data analysis. Background Technology

[0002] The impact of railway tunnel construction on roof vegetation refers to the damage or deterioration of growing conditions on the surface vegetation in the area above the tunnel (i.e., the tunnel roof) caused by engineering activities such as blasting, excavation, earthwork transportation, and changes in groundwater levels during tunnel construction. This manifests as decreased vegetation cover, damaged plant roots, loosened soil structure, and reduced water resources, thereby affecting the stability of the ecosystem and the natural recovery capacity of the vegetation. This impact is usually an important aspect of the ecological and environmental impact assessment of engineering projects.

[0003] In vegetation impact assessments of railway tunnels, traditional GIS buffer zone analyses are often based on two-dimensional planar construction, neglecting the influence of terrain elevation differences on actual distances, leading to the problem of "vertical projection error." Especially in mountainous or plateau regions, tunnel roofs often extend along steep slopes. If a buffer zone is set at a horizontal distance (e.g., 100m), the actual surface distance can increase significantly due to the slope (e.g., 141m at a 45° slope), causing the buffer zone to shift or "miss" into non-target areas, resulting in misjudgments or omissions of the affected vegetation area. This error is particularly severe in areas with dramatic terrain undulations, directly affecting the accuracy of environmental impact conclusions and the allocation of responsibility for subsequent ecological restoration. It represents a key shortcoming in existing technologies regarding insufficient three-dimensional topographic processing capabilities. Summary of the Invention

[0004] The purpose of this invention is to provide an assessment method for the impact of railway tunnel engineering on the vegetation above the tunnel roof based on comprehensive data analysis, so as to overcome the shortcomings of the prior art.

[0005] To achieve the above objectives, the present invention provides the following technical solution: a method for assessing the impact of railway tunnel engineering on tunnel roof vegetation based on comprehensive data analysis, comprising:

[0006] Obtain tunnel line vector data and digital elevation model data of the study area for railway tunnel engineering, and set the two-dimensional foundation buffer radius D0;

[0007] Based on digital elevation model data, the slope aspect and slope gradient within the study area are calculated, and the slope aspect-slope coupling orientation factor (SACDI) is further calculated.

[0008] Based on the surface elevation changes in the buffer boundary region, a buffer elevation gradient offset index (EGDF) is constructed to quantify the distance deviation between the region and the ideal plane.

[0009] After normalizing the SACDI and EGDF parameters, a dynamic buffer distance function is constructed, with the expression: D=D0·(1+α·SACDI+β·EGDF); where D represents the actual buffer distance of each pixel, and α and β are adjustment weight coefficients.

[0010] Based on the actual buffer distance, spatial expansion is performed with the tunnel line as the center to generate the buffer boundary, and the buffer result is vectorized to generate a dynamic buffer vector layer.

[0011] Multi-temporal remote sensing images are overlaid and vegetation index differences are calculated. Vegetation change trends at different distance levels within the buffer zone are extracted, and error analysis maps and vegetation impact heat maps are output.

[0012] Preferably, the process involves acquiring tunnel line vector data for the railway tunnel project and digital elevation model data of the study area, and setting a two-dimensional foundation buffer radius D0, specifically including:

[0013] Import the centerline or structural outline of the railway tunnel project into the GIS system, perform topology checks and simplification on the tunnel line vector data, and redefine the directional attribute from the tunnel exit point to the entrance point.

[0014] A digital elevation model with a resolution of no less than 10 meters was selected as the topographic data source for the study area, and the digital elevation model was reprojected to the same projection coordinate system as the tunnel line data using an equidistant projection method.

[0015] Based on the tunnel design depth, terrain undulation intensity, and regional vegetation cover, a two-dimensional basic buffer radius D0 is set by constructing a multi-factor linear weighted model.

[0016] Preferably, the calculation of the aspect-slope coupling orientation factor (SACDI) includes:

[0017] Before performing slope aspect calculations on the digital elevation model (DEM), the DEM is preprocessed using a terrain shadow filtering algorithm with directional weighted kernels.

[0018] A multi-scale sliding window joint calculation strategy for surface slope is adopted to extract the slope at different spatial scales, and the results of each scale are fused by local variance weighting to obtain the slope value.

[0019] The difference between the slope aspect value and the tunnel axis direction of each pixel is discretely classified and combined with its slope level to construct a two-dimensional directional response table. Each slope aspect-slope combination is assigned a directional level score, namely the slope aspect-slope coupling directional factor SACDI.

[0020] Preferably, the construction of the buffer elevation gradient offset index EGDF includes:

[0021] Extend several pixels outward from the boundary of the two-dimensional basic buffer along the normal direction, and extract a continuous sequence of elevation values ​​on each extension path;

[0022] For each buffer profile path, the ratio of continuous elevation difference to actual slope distance is calculated using a sliding window method;

[0023] Within the buffer zone, a stable terrain area is selected to construct a "reference elevation profile model" as the standard terrain gradient. The elevation gradient values ​​on each buffer path are compared with the reference profile model. The Z-score is used to quantify the degree of deviation, and the deviation values ​​are normalized to a dimensionless index between 0 and 1 to obtain the elevation gradient offset index EGDF.

[0024] Preferably, based on the actual buffer distance, spatial expansion is performed with the tunnel line as the center to generate a buffer boundary, and the buffering result is vectorized to generate a dynamic buffer vector layer, specifically including:

[0025] After obtaining the actual buffer distance value of each pixel, the buffer path is decomposed based on the slope direction, that is, the expansion direction of each pixel is limited to a fan-shaped range of ±45 degrees of its main slope direction.

[0026] The initially generated buffer raster layer is overlaid with the digital elevation model for analysis, and an elevation range control threshold is introduced, that is, the maximum acceptable elevation difference is set as the limit condition for the effective buffer area, and areas with significant elevation difference breaks within the same buffer distance are eliminated.

[0027] Apply a spatial clustering algorithm based on four or eight neighborhoods to the buffer raster layer to remove isolated pixel patches with an area smaller than a preset threshold, and perform boundary correction to generate buffer patch boundaries; convert the spatially corrected buffer raster layer into a vector polygon layer.

[0028] Preferably, the method involves overlaying multi-temporal remote sensing images and calculating vegetation index differences, extracting vegetation change trends at different distance levels within the buffer zone, and outputting error analysis maps and vegetation impact heat maps, specifically including:

[0029] The dynamic buffer is partitioned according to distance level, and multi-temporal NDVI values ​​are extracted for each segment. The NDVI change trend is calculated by combining the median difference method and the linear trend fitting method. The median difference reflects the abrupt change, and the linear fitting reflects the long-term trend.

[0030] The NDVI change values ​​calculated in each level of the buffer zone, together with the SACDI and EGDF topographic factors, are used as input weighting factors. A weighted interpolation method combining inverse distance weighting and variance-constrained Kriging is used to generate a vegetation impact heat map with spatial continuity. The color gradient in the map reflects the intensity, direction, and propagation trend of construction disturbance.

[0031] The technical effects and advantages provided by the present invention in the above technical solution are as follows:

[0032] 1. The proposed method for assessing the impact of railway tunnel engineering on roof vegetation effectively overcomes the technical bottlenecks of traditional two-dimensional buffer methods in complex terrain by introducing the aspect-slope coupled orientation factor (SACDI) and the buffer elevation gradient offset index (EGDF), combined with a dynamic buffer distance adjustment function. This method enables the buffer to adaptively expand according to the actual slope direction and terrain disturbance trend, significantly improving the spatial simulation accuracy of construction impact paths. It is particularly suitable for ecological assessment tasks in areas with drastic changes in landforms, such as mountains and plateaus.

[0033] 2. This invention achieves a quantitative expression of the intensity, extent, and direction of vegetation degradation by employing precise registration of multi-temporal remote sensing images, identifying hierarchical NDVI difference trends within buffer zones, introducing a baseline differential correction mechanism, and constructing a high-resolution vegetation impact heat map. This method not only possesses excellent automation and visualization capabilities but also provides a scientific basis for subsequent ecological restoration zone delineation, risk warning, and engineering management. Attached Figure Description

[0034] To more clearly illustrate the technical solutions in the embodiments of this application or the prior art, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments recorded in this invention. For those skilled in the art, other drawings can be obtained based on these drawings.

[0035] Figure 1 This is a flowchart of the method of the present invention. Detailed Implementation

[0036] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0037] For examples, please refer to Figure 1 As shown in this embodiment, the method for assessing the impact of railway tunnel engineering on the roof vegetation based on comprehensive data analysis includes:

[0038] Obtain tunnel line vector data and digital elevation model data of the study area for railway tunnel engineering, and set the two-dimensional foundation buffer radius D0;

[0039] Based on digital elevation model data, the slope aspect and slope gradient within the study area are calculated, and the slope aspect-slope coupling orientation factor (SACDI) is further calculated.

[0040] Based on the surface elevation changes in the buffer boundary region, a buffer elevation gradient offset index (EGDF) is constructed to quantify the distance deviation between the region and the ideal plane.

[0041] After normalizing the SACDI and EGDF parameters, a dynamic buffer distance function is constructed, with the expression: D=D0·(1+α·SACDI+β·EGDF); where D represents the actual buffer distance of each pixel, and α and β are adjustment weight coefficients.

[0042] Based on the actual buffer distance, spatial expansion is performed with the tunnel line as the center to generate the buffer boundary, and the buffer result is vectorized to generate a dynamic buffer vector layer.

[0043] Multi-temporal remote sensing images are overlaid and vegetation index differences are calculated. Vegetation change trends at different distance levels within the buffer zone are extracted, and error analysis maps and vegetation impact heat maps are output.

[0044] Railway tunnel engineering design data usually comes from CAD drawings, engineering survey results, or geological exploration databases. The data formats are diverse and there are problems such as inconsistent coordinate systems, non-topological connections of nodes, overlapping or offset line segments, and missing directional attributes. If they are not standardized, they will cause systematic errors in subsequent calculations based on directional parameters (such as aspect-slope coupling factor) or buffer paths.

[0045] This invention first imports the tunnel line layer (usually in polyline or polyline format) into a GIS platform (such as ArcGIS, QGIS, or MapGIS) and performs the following standardized operations:

[0046] Topology cleanup: Repairing breakpoints, eliminating duplicate nodes, and unifying line segment connections using topology tools ensures that the tunnel line has a continuous topology structure.

[0047] Geometric simplification: Merge line segments within the tolerance range for non-critical detail points (such as architectural annotation nodes) to reduce the burden of spatial data processing.

[0048] Directional attribute assignment: Extract the "inlet-exit" orientation information clearly defined in the tunnel design drawing, automatically calculate the azimuth angle of each line segment, and uniformly adjust the direction to ensure that the starting point orientation of all line segments is consistent with the tunnel axis, providing a basis for subsequent calculation of directional factors.

[0049] Through the above standardization process, the tunnel line layer has good geometric consistency and attribute standardization, which significantly improves the model accuracy and automated processing capabilities.

[0050] Digital elevation models (DEMs) are the core data for constructing topographic factors such as slope, aspect, and surface roughness. Since railway tunnel projects often traverse complex terrain areas such as mountains and hills, the choice of DEM resolution has a significant impact on the analysis results.

[0051] This invention recommends using high-precision DEM data with a spatial resolution of at least 10 meters, such as ASTER GDEM, ALOSPALSAR, and domestic elevation maps (e.g., DEM extracted from 1:5000 topographic maps). To ensure the accuracy and consistency of the spatial analysis, the implementation steps include:

[0052] DEM accuracy verification: Set up a comparison sample area around the tunnel to check the difference between the elevation represented by the DEM and the measured or topographic map data. The error should be within ±2 meters.

[0053] Unified projection system: The DEM data is reprojected to a projection coordinate system consistent with the tunnel line data. To ensure the accuracy of distance calculation, this invention preferably uses an equidistant projection system, such as UTM (Universal Transverse Mercator) or Gauss-Kruger projection, to avoid distance deviations caused by angular projection distortion.

[0054] DEM and Vector Line Alignment Correction: Use the "Align to Vector Boundary" tool to re-register the DEM raster, ensuring that the raster cell center falls precisely within the tunnel line and its buffer zone, reducing spatial overlay errors.

[0055] The above processing ensured the spatial accuracy matching between the elevation data and the tunnel engineering line, laying the foundation for subsequent buffer space construction, slope analysis, and terrain response modeling.

[0056] Traditional buffer zone analysis methods typically set a fixed radius (e.g., 100 meters, 200 meters) as the spatial analysis boundary. However, in areas with dramatic topographic relief or where tunnels are deeply buried and vary greatly, a fixed radius may lead to an underestimation or overestimation of the vegetation impact area, affecting the scientific validity of the assessment results.

[0057] To address this issue, this invention designs a basic buffer radius setting model based on a multi-factor combination, comprehensively considering the following three core variables:

[0058] Tunnel depth (H): This refers to the designed depth of the tunnel. The greater the designed depth of the tunnel, the more likely the terrain undulations above the tunnel roof will have an impact on vegetation, and the buffer radius needs to be expanded.

[0059] Terrain complexity (T): The standard deviation of slope within the buffer zone is calculated based on the DEM. Higher terrain complexity increases the risk of slope propagation path deviation.

[0060] Regional vegetation sensitivity (V): The average NDVI value is calculated from multiple remote sensing images to reflect the intensity of regional vegetation response to external disturbances (such as blasting and precipitation changes).

[0061] The initial buffer radius is equal to the baseline distance multiplied by a correction factor weighted by tunnel depth, terrain complexity, and vegetation sensitivity. The weights can be set by expert experience or obtained through fitting historical cases. For example, if the baseline buffer is 100 meters, the buffer radius can automatically expand to the 120–150 meter range when the tunnel is deeper, the terrain is rugged, or the vegetation sensitivity is high.

[0062] By dynamically setting the initial buffer distance, the adaptability of the method under different ecological zones and landform types can be greatly improved, avoiding over-reliance on human experience or fixed thresholds.

[0063] Digital elevation models (DEMs) are the primary data source for slope aspect calculations. However, due to factors such as high-frequency undulations in the terrain surface, data sampling errors, or disturbances caused by manual cutting and filling, local elevation noise often exists in the original DEM, leading to abrupt or irregular jumps in slope aspect results at valley, ridge, or steep slope boundaries. This error is particularly prominent when assessing tunnel impact paths, easily causing misjudgments.

[0064] Therefore, this invention performs terrain shadow direction filtering preprocessing on the DEM before slope aspect calculation. This method constructs a direction-weighted kernel based on the dominant illumination direction (e.g., the azimuth of the sun at typical observation times, set to 135 degrees). This kernel assigns higher weights to the dominant slope aspect direction (e.g., southeast) and lower weights to the opposite direction (e.g., northwest). By performing direction kernel convolution processing on the DEM data, the topographic structure in specific directions is "enhanced" while suppressing elevation disturbances in non-dominant directions.

[0065] Compared with traditional median or Gaussian filtering, this method effectively eliminates slope aspect jumps caused by terrain shading and gully interference while preserving the main slope features of the mountain, thus improving the spatial continuity of slope aspect extraction and the expressive power of mountain structure.

[0066] Slope is a key factor reflecting the degree of inclination of a slope in SACDI calculations. Generally, slope calculations use a fixed window (e.g., 3×3 pixels) to perform local elevation difference analysis on a DEM. However, in real terrain, features such as gentle slopes, steep slopes, and cliffs coexist, and single-scale calculations can easily lead to distorted results.

[0067] To improve the robustness of slope calculation, this invention employs a multi-scale sliding window fusion algorithm, the steps of which are as follows:

[0068] First, the slope is extracted using a 3×3 window to capture detailed terrain changes;

[0069] Then, a large-scale slope is extracted using a 7×7 window to reflect the overall slope direction;

[0070] The local elevation variances of the results at both scales are calculated separately and used as the basis for weighting.

[0071] Finally, the "fusion slope value" of each pixel is obtained by fusion of multi-scale results through weighted average.

[0072] This method overcomes the sensitivity of local noise to slope, making slope extraction for both gentle slopes and undulating terrains highly adaptable and expressive, providing a more stable foundation for subsequent directional computing.

[0073] After extracting stable slope aspect and slope, this invention does not directly calculate the SACDI value in the form of a continuous function, but instead adopts a more interpretable "response mapping table construction method".

[0074] This method is implemented through the following process:

[0075] Slope Aspect Difference Classification: The difference between the slope aspect of each pixel and the tunnel axis direction is used to calculate the "slope aspect angle". Based on the range of the angle, it is divided into 8 directional categories, including "positive" (0–22.5°), "with" (22.5–45°), "biased" (45–67.5°), "lateral" (67.5–112.5°), "biased" (112.5–135°), "opposite" (135–157.5°), and "opposite" (157.5–180°).

[0076] Slope classification: Based on the integrated slope value, the slope is divided into three levels: "gentle slope" (slope less than 8 degrees), "medium slope" (8-25 degrees), and "steep slope" (greater than 25 degrees);

[0077] Construct a two-dimensional mapping table: Create a two-dimensional cross table of direction category and slope level. Each combination (such as "forward + medium slope", "backward + steep slope", etc.) is assigned a directional response level value, ranging from 0.1 to 1.0, which represents the spatial response tendency of the slope to tunnel disturbance under the given direction and slope conditions.

[0078] Input the slope aspect category and slope grade of each cell into the response table, and obtain its corresponding grade score, which is the SACDI value of that cell.

[0079] This method has the following advantages:

[0080] It does not rely on complex function modeling, is easy to interpret and manually correct; it is flexibly extensible, supports expert experience and machine learning training optimization; and significantly enhances the ability to identify directions in heterogeneous terrain (such as the back slope of a mountain).

[0081] Conventional methods for calculating elevation gradients often employ a fixed neighborhood sampling method centered on a pixel (such as a 3×3 or 5×5 window) to calculate elevation differences along the main direction or any other direction. However, in buffer boundary analysis, the terrain expansion path should primarily follow the normal direction of the buffer boundary to simulate the natural propagation trend of disturbances on the terrain slope.

[0082] Therefore, this invention proposes a Layered Directional Expansion algorithm, which takes each cell on the boundary of the buffer zone as the starting point, calculates its spatial normal direction (i.e., the vertical direction from the boundary to the outside), and expands several layers (such as 3 to 5 cells) in this direction to extract continuous elevation values ​​and form a "topographic profile sequence" along the slope direction.

[0083] This method simulates the path propagating outward from the top of the tunnel while preserving the spatial properties of the actual slope trend. Compared to radial buffers or circular neighborhoods, this method can reveal local abrupt changes and trend variations in elevation along the slope direction more precisely, and is particularly suitable for complex terrain areas such as ridges, valleys, and cliffs.

[0084] After obtaining the directional terrain profile, this invention further processes it to extract an elevation gradient index representing the true slope characteristics of the path. Unlike the traditional simple slope calculation of dividing the elevation difference by the horizontal distance, this invention employs the following strategy:

[0085] Each profile path is treated as a one-dimensional high-order sequence;

[0086] By employing a sliding window mechanism (e.g., grouping 5 pixels together), the ratio of the maximum elevation difference within the window to the path length is calculated to obtain the local elevation change rate within the window.

[0087] The entire profile sequence was smoothed using the Savitzky–Golay derivative filter, a signal smoothing method based on local polynomial fitting. This filter can preserve the trend of data variation and remove random noise. It is widely used in spectral curves and biological signal processing, but it is still rare in DEM data slope profile processing.

[0088] Through the above processing, interference caused by abnormal elevation differences or cut and fill data can be effectively suppressed without sacrificing the fidelity of slope features, making the elevation gradient extraction more stable and accurate, especially suitable for areas with large slopes or mixed terrain.

[0089] To transform elevation gradient deviations from different paths into standardized and comparable EGDF indices, this invention introduces the concept of a "reference profile model." This model is constructed within a buffer zone characterized by consistent slope aspect, stable slope, and minimal topographic changes (such as stable slopes or platforms) to represent elevation variation patterns under ideal terrain conditions, serving as a benchmark.

[0090] The specific implementation is as follows:

[0091] Multiple profile paths are extracted from the candidate region, and their average elevation gradient values ​​are calculated. The measured elevation gradient of each buffer path is compared with the reference value. The degree of deviation is quantified using standardized deviation rate indicators, such as: using Z score (i.e., measured value minus reference mean, then divided by standard deviation); using percentile deviation rate (i.e., the position of the measured value in the reference distribution).

[0092] Finally, the offset results are converted into dimensionless values ​​between 0 and 1, forming the EGDF index. The closer the value is to 1, the more severe the elevation disturbance in the slope direction of the buffer path, the more obvious the terrain deviation, and the greater the risk of error in buffer construction.

[0093] EGDF, as a key parameter used in this invention to identify two-dimensional buffer errors, can be directly used as a weight term in the dynamic buffer adjustment function. For example, in the dynamic buffer distance calculation model, the EGDF value is used as a multiplication factor to affect the extent of local buffer radius expansion, realizing the terrain-responsive expansion logic of "the more complex the slope, the wider the buffer".

[0094] Furthermore, EGDF values ​​can also be used to: construct error heatmaps for visualizing error risk distribution; provide a spatial interpretation basis for areas with abnormal vegetation changes; and support the construction of ecological restoration priority ranking or risk warning models. Compared with traditional single slope factors, EGDF more comprehensively expresses the trend of elevation disturbance through profile modeling and directional extension, exhibiting stronger interpretability, sensitivity, and engineering adaptability. It is particularly suitable for the impact analysis of railway tunnels in special geomorphological areas such as steep slopes, plateaus, and fault zones.

[0095] This invention designs a dynamic adjustment model for buffer distance based on terrain response factors. The model uses the aspect-slope coupling orientation factor (SACDI) and the elevation gradient offset factor (EGDF) as driving parameters to construct a buffer distance function that can adapt to terrain changes and realize the "local" calculation of buffer boundaries.

[0096] Specifically, the model first normalizes SACDI and EGDF, unifying their values ​​to between 0 and 1 for easier participation in subsequent weighted calculations. Normalization methods can include linear scaling using maximum and minimum values, or Z-score standard deviation normalization. After normalization, the buffer distance function is defined as follows:

[0097] The final buffer distance for a given pixel is equal to the base buffer distance multiplied by a dynamic coefficient jointly controlled by SACDI and EGDF, specifically expressed as: D = D0·(1 + α·SACDI + β·αEGDF); where D represents the actual buffer distance for each pixel, and α and β are adjustment weight coefficients; the base buffer distance is the initial setting value (e.g., 100 meters), and α and β are adjustment coefficients used to control the influence weight of SACDI and EGDF on the result. The two weight coefficients can be fitted using historical engineering data, set using expert experience, or adaptively optimized based on machine learning. A larger coefficient indicates a stronger contribution of that factor to buffer expansion, while a coefficient of zero indicates that the factor does not participate in the adjustment.

[0098] This dynamic function has the following characteristics and advantages:

[0099] Directional responsiveness: The higher the SACDI value, the more the slope is oriented towards the tunnel axis, and the more direct the disturbance path is. Therefore, the buffer radius should be appropriately increased.

[0100] Elevation compensation: The higher the EGDF value, the greater the local terrain deviation and the more significant the two-dimensional buffer error, which needs to be corrected by stretching the buffer boundary;

[0101] Spatial adaptability: The buffer zone is no longer a regular circle or ring, but automatically deforms along the slope structure to form an "ecological response zone" that better fits the terrain;

[0102] Each extended result is controlled by well-defined parameters, facilitating analysis and debugging. The buffer layer generated after applying this function will be used for subsequent NDVI vegetation change analysis, risk heat map construction, and ecological restoration scope delineation, providing a more accurate spatial boundary reference for the ecological impact assessment of tunnel engineering.

[0103] Traditional GIS buffering operations default to an equidistant circular expansion method, which fails to reflect the directional orientation of actual terrain slopes. However, in mountainous tunnel scenarios, surface disturbances often propagate downwards or forwards along the slope, exhibiting a clear directional preference.

[0104] Therefore, this invention introduces a "slope direction-restricted expansion" strategy. Based on the previously calculated slope aspect value, the expansion direction of each pixel is restricted when constructing the buffer path, generating the buffer path only within an angle range of ±45 degrees of its main slope aspect, thus forming a "directional buffer sector". This method constructs a locally direction-sensitive expansion domain around each pixel, avoiding the problems of "reverse slope" or "slope crossing" during buffer expansion at steep slopes, faults, or slope inflections.

[0105] The generated buffer raster layer exhibits irregular geometric boundaries that conform to the slope morphology, which is more consistent with the actual ecological propagation path and provides a structural basis for subsequent geomorphological constraint processing.

[0106] Changes in slope morphology are often accompanied by abrupt changes in elevation. For example, there may be a difference of tens of meters in elevation between a valley and a ridge over a short distance. In this context, equidistant expansion can easily cross natural landform boundaries, mistakenly including areas of extreme elevation difference within the buffer zone, resulting in spatial "jumps" in vegetation impact assessment.

[0107] To overcome this problem, this invention proposes an elevation difference limiting mechanism, which sets a maximum acceptable elevation difference threshold (e.g., ±30 meters) to cull and trim the buffer area after grid expansion. The specific implementation is as follows:

[0108] The initially generated buffer raster is overlaid on the original DEM;

[0109] Calculate the elevation difference between each cell in the buffer and the nearest point on the tunnel line;

[0110] If the elevation difference exceeds the set threshold, it will be marked as a "terrain barrier area" and removed.

[0111] The remaining area constitutes a "continuous terrain buffer zone".

[0112] This treatment effectively prevents spatial errors such as buffer zones "crossing valleys" or "jumping to ridges," making the buffer zones more ecologically sound and feasible.

[0113] After generating a buffer raster based on slope orientation and elevation constraints, problems such as pixel discretization, jagged boundaries, or local patchiness may occur, especially when the slope is broken or there are abrupt change zones.

[0114] To improve the geometric continuity and patch stability of the buffer layer, this invention introduces a "neighbor clustering and boundary smoothing" strategy, including:

[0115] Pixel clustering analysis was performed using an eight-neighbor connection method to identify all buffer patches;

[0116] Remove isolated patches with an area smaller than a threshold (e.g., 20 pixels or 200 square meters);

[0117] Smooth the boundaries of the main buffer patches, such as by applying the Moran boundary filter or the Ramer–Douglas–Peucker simplification algorithm, to remove sharp corners and jagged edges and improve the regularity of the shape.

[0118] This process ensures the continuity, aesthetics, and practical usability of the final buffer boundary, facilitating the overlay of vegetation data for precise spatial analysis.

[0119] After completing the spatial construction of the buffer patch, in order to meet the requirements of subsequent analysis, mapping and modeling, the buffer layer needs to be converted from raster format to vector format and relevant attribute information needs to be integrated.

[0120] The specific process is as follows:

[0121] Perform a "Raster to Polygon" operation on the buffer raster layer to generate a closed polygon layer;

[0122] For each polygonal buffer patch, calculate its average SACDI value, EGDF value, mean slope, NDVI change value and other statistical indicators;

[0123] Write the above parameters into the vector layer attribute table in the form of fields;

[0124] If multiple buffer levels exist (e.g., 0–100 meters, 100–200 meters), they are labeled with the “Buffer Level” field, supporting tiered evaluation.

[0125] The final output vector layer not only provides complete spatial boundary information, but also contains multi-source data attributes, with extremely high data fusion. It can be directly used in multiple application scenarios such as risk heat map construction, vegetation change comparative analysis, and ecological restoration area priority division.

[0126] In vegetation change monitoring, spatial consistency of multi-temporal remote sensing images is a prerequisite for ensuring the accuracy of NDVI calculation. Especially in mountainous or complex terrain areas, there may be sub-pixel level offset errors between images, which will manifest as "false change bands" in NDVI calculation.

[0127] To address this problem, this invention employs a temporal image registration method that fuses sparse keypoint matching with dense optical flow analysis. The steps are as follows:

[0128] Extract topographic stability features (such as rocks, roads, and bare land) from images of different time phases as key points;

[0129] Initial registration was performed using the Scale Invariant Feature Transform (SIFT) algorithm;

[0130] The matching relationship in local regions is refined using dense optical flow algorithms such as Farneback or PWC-Net;

[0131] Register all images to a uniform reference (such as pre-construction images) to ensure strict alignment between pixels.

[0132] This method is more accurate than traditional affine or projection transformation methods, and is particularly suitable for registration of mountain remote sensing images with local deformation, laying a high-quality foundation for NDVI difference calculation.

[0133] To enable comparative analysis of vegetation changes within different distance zones, this invention divides the previously constructed dynamic buffer layer into several concentric hierarchical buffer zones, such as 0–100m, 100–200m, 200–300m, etc.

[0134] Within each range band, multi-temporal NDVI values ​​are extracted and processed as follows:

[0135] The NDVI median difference method (i.e., the median of the later period minus the median of the earlier period) is used to quickly identify areas of sudden disturbance (such as the blasting impact zone);

[0136] For images of three or more phases, the linear trend fitting method is used to calculate the temporal slope of NDVI and identify long-term degradation trends.

[0137] By combining the results of the two methods, a two-dimensional discrimination matrix of "intensity of change × direction of change" is constructed to classify pixels within the buffer zone (e.g., rapid degradation, stable recovery, large fluctuations, etc.).

[0138] It can output a statistical table of NDVI change trends, which can be summarized by distance level and vegetation change type.

[0139] This analytical method effectively integrates abrupt change detection with trend modeling, enhancing the ability to identify complex vegetation response mechanisms.

[0140] Changes in remote sensing indices not only reflect engineering disturbances but also include fluctuations caused by natural variations (such as seasonal, climatic, and radiation condition changes). Without correction, assessment results may confuse "natural degradation" with "human-caused damage."

[0141] Therefore, this invention proposes a differential correction method for a reference control region. Its implementation process is as follows:

[0142] Several "control areas" with similar ecological backgrounds and unaffected by construction were selected outside the study area;

[0143] The trend of NDVI variation in the control area over the same time period was calculated to obtain the "natural background fluctuation curve".

[0144] Using this curve as a baseline, differential correction is uniformly performed on the NDVI change values ​​in the buffer zone;

[0145] After removing background variation components, the remaining NDVI difference is the "net engineering impact" index.

[0146] This mechanism significantly improves the correlation between NDVI changes and actual disturbances, and solves the long-standing misjudgment problem of "change does not equal intervention" in ecological remote sensing monitoring.

[0147] To facilitate spatial representation and assist in the design of ecological restoration strategies, this invention integrates the NDVI difference results with various terrain and interference factors within the buffer zone in a weighted manner to construct a continuous vegetation impact heat map.

[0148] The specific process is as follows:

[0149] The NDVI difference is used as the main indicator, combined with SACDI (slope aspect response) and EGDF (elevation disturbance) as weighting parameters.

[0150] Inverse distance weighted interpolation (IDW) is used to perform preliminary spatial diffusion of NDVI values;

[0151] For regions with significant elevation differences, the variance-constrained Kriging method is introduced to correct the results, ensuring that the spatial statistical model is not skewed by outliers;

[0152] Output a raster-format heatmap, using hierarchical color coding to express the intensity of the impact (e.g., red indicates strong degradation, blue indicates stable recovery);

[0153] It can overlay vegetation type, land use or administrative boundary layers for the creation of ecological restoration zoning and risk warning maps.

[0154] This heat map not only provides a visual representation of the intensity of disturbance, but also supports the export of high-impact area patches by level, providing intuitive decision support for ecological engineering.

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

Claims

1. A method for assessing the impact of railway tunnel engineering on roof vegetation based on comprehensive data analysis, characterized by: include: Acquire tunnel line vector data and digital elevation model data of the study area for railway tunnel engineering, and set the two-dimensional foundation buffer radius. ; Based on digital elevation model data, the slope aspect and slope gradient within the study area are calculated. The slope aspect-slope coupling orientation factor (SACDI) is further calculated, including: Before performing slope aspect calculations on the digital elevation model (DEM), the DEM is preprocessed using a terrain shadow filtering algorithm with directional weighted kernels. A multi-scale sliding window joint calculation strategy for surface slope is adopted to extract the slope at different spatial scales, and the results of each scale are fused by local variance weighting to obtain the slope value. The difference between the slope aspect value and the tunnel axis direction of each pixel is discretized and classified. A two-dimensional directional response table is constructed by combining the slope level and assigning a directional level score to each slope aspect-slope combination, namely the slope aspect-slope coupling directional factor SACDI. Based on the surface elevation changes in the buffer zone boundary area, an elevation gradient offset index (EGDF) is constructed to quantify the distance deviation between this area and the ideal plane, including: Extend several pixels outward from the boundary of the two-dimensional basic buffer along the normal direction, and extract a continuous sequence of elevation values ​​on each extension path; For each buffer profile path, the ratio of continuous elevation difference to actual slope distance is calculated using a sliding window method; Within the buffer zone, a stable terrain area is selected to construct a reference elevation profile model as the standard terrain gradient. The elevation gradient values ​​on each buffer path are compared with the reference profile model. The Z-score is used to quantify the degree of deviation, and the deviation values ​​are normalized to a dimensionless index between 0 and 1 to obtain the elevation gradient offset index EGDF. After normalizing the SACDI and EGDF parameters, a dynamic buffer distance function is constructed, with the following expression: In the formula, This represents the actual buffer distance for each pixel. To adjust the weighting coefficients; Based on the actual buffer distance, spatial expansion is performed with the tunnel line as the center to generate a buffer boundary, and the buffering result is vectorized to generate a dynamic buffer vector layer, including: After obtaining the actual buffer distance value of each pixel, the buffer path is decomposed based on the slope direction, that is, the expansion direction of each pixel is limited to a fan-shaped range of ±45 degrees of its main slope direction. The initially generated buffer raster layer is overlaid with the digital elevation model for analysis, and an elevation range control threshold is introduced, that is, the maximum acceptable elevation difference is set as the limit condition for the effective buffer area, and areas with significant elevation difference breaks within the same buffer distance are eliminated. Apply a spatial clustering algorithm based on four or eight neighborhoods to the buffer raster layer to remove isolated pixel patches with an area smaller than a preset threshold, and perform boundary correction to generate buffer patch boundaries; convert the spatially corrected buffer raster layer into a vector polygon layer. Multi-temporal remote sensing images are overlaid and vegetation index differences are calculated. Vegetation change trends at different distance levels within the buffer zone are extracted, and error analysis maps and vegetation impact heat maps are output.

2. The method for assessing the impact of railway tunnel engineering on roof vegetation based on comprehensive data analysis according to claim 1, characterized in that: Acquire tunnel line vector data and digital elevation model data of the study area for railway tunnel engineering, and set the two-dimensional foundation buffer radius. Specifically, it includes: Import the centerline or structural outline of the railway tunnel project into the GIS system, perform topology checks and simplification on the tunnel line vector data, and redefine the directional attribute from the tunnel exit point to the entrance point. A digital elevation model with a resolution of no less than 10 meters was selected as the topographic data source for the study area, and the digital elevation model was reprojected to the same projection coordinate system as the tunnel line data using an equidistant projection method. Based on the tunnel design depth, terrain undulation intensity, and regional vegetation cover, a two-dimensional basic buffer radius is determined by constructing a multi-factor linear weighted model. .

3. The method for assessing the impact of railway tunnel engineering on roof vegetation based on comprehensive data analysis according to claim 1, characterized in that: By overlaying multi-temporal remote sensing images and calculating vegetation index differences, the vegetation change trends of different distance levels within the buffer zone are extracted, and error analysis maps and vegetation impact heat maps are output, specifically including: The dynamic buffer is partitioned according to distance level, and multi-temporal NDVI values ​​are extracted for each segment. The NDVI change trend is calculated by combining the median difference method and the linear trend fitting method. The median difference reflects the abrupt change, and the linear fitting reflects the long-term trend. The NDVI change values ​​calculated in each level of the buffer zone, together with the SACDI and EGDF topographic factors, are used as input weighting factors. A weighted interpolation method combining inverse distance weighting and variance-constrained Kriging is used to generate a vegetation impact heat map with spatial continuity. The color gradient in the map reflects the intensity, direction, and propagation trend of construction disturbance.

Images