A method for calculating the area of forest encroached by ski resorts based on adaptive envelope reconstruction

Abstract

The application provides a ski field encroachment forest area measurement method based on adaptive envelope reconstruction, relates to the technical field of geographic spatial information processing and remote sensing monitoring, and the method is characterized in that: forest binary raster data and ski field boundary vector data are used for spatial superposition, a spatially constrained forest pixel set located in the ski field range is obtained, a spatial boundary is extracted and discrete point sampling is performed to obtain a unique boundary point set; according to the unique boundary point set and digital elevation data, local terrain slope and three-dimensional spatial distance are calculated to obtain an effective edge set; the effective edge set is used to construct a plane topological graph structure and perform closed path tracking to obtain an external envelope surface representing the continuous spatial distribution of the forest; based on the forest pixel set and the spatial resolution thereof, the actual area of the forest is calculated to obtain the ski field encroachment forest area. The application solves the problems of inaccurate forest extraction and lack of distinction between artificial and natural boundary features, thereby avoiding the problem of morphological reconstruction distortion.

Description

Technical Field

[0001] This invention relates to the fields of geospatial information processing and remote sensing monitoring technology, and in particular to a method for calculating the forest encroachment area of ​​ski resorts based on adaptive envelope reconstruction. Background Technology

[0002] In recent years, with the rapid development of the ice and snow industry, the construction scale of mountain ski resorts has been continuously expanding. The construction and operation of ski resorts often require excavation of mountains, construction of ski runs, and installation of cable cars, which inevitably leads to the felling and destruction of natural forests. In order to objectively assess the damage to the ecological environment and provide a scientific and quantitative technical basis for off-site compensation of forest land or forestry law enforcement, it is urgent to accurately calculate the area of ​​forest encroachment caused by the construction of ski resorts.

[0003] However, existing methods for calculating forest encroachment area mainly rely on manual field surveys or traditional basic remote sensing classification and area comparison techniques, which have significant shortcomings when dealing with complex mountain engineering projects such as ski resorts. First, during the remote sensing information extraction stage, ski resort areas often exhibit extreme conditions where high-albedo snow backgrounds intertwine with tree shadows, creating numerous mixed edges of half-trees and half-snow. Traditional optical remote sensing classification methods struggle to effectively suppress this strong background interference, easily leading to misclassification or omissions, resulting in low accuracy in basic forest boundary extraction.

[0004] Secondly, in the boundary morphology restoration stage, existing technologies often directly extract edge points using a uniform sampling interval, lacking an effective distinction between artificially logged boundaries and naturally grown boundaries. The edges of artificially logged snow runs are typically smooth and uniformly oriented, while the edges of natural forest belts exhibit high-frequency random fluctuations due to topographical influences. Traditional methods cannot adapt to these two characteristics, resulting in the reconstructed forest outline failing to accurately restore its undisturbed original state.

[0005] Finally, and most importantly, ski resorts are typically located in mountainous environments with significant elevation differences. Existing methods for reconstructing continuous forest surfaces mainly employ two-dimensional planar envelope algorithms (such as the conventional two-dimensional Alpha Shape algorithm). On steep mountain slopes, the actual physical distances are severely compressed and distorted when projected onto a two-dimensional plane; furthermore, conventional algorithms lack constraints on elevation-fractured terrain, easily generating erroneous lines crossing cliffs. This results in a significant discrepancy between the theoretical forest area calculated based on the two-dimensional projection surface and the actual mountain surface area, leading to a severely distorted final difference in the encroached area, making it unsuitable as a high-precision quantitative assessment basis. Summary of the Invention

[0006] To address the aforementioned shortcomings in existing technologies, this invention provides a method for calculating the forest encroachment area of ​​ski resorts based on adaptive envelope reconstruction. This method solves the problems in existing technologies, such as inaccurate forest extraction due to snow accumulation interference from ski resorts, morphological reconstruction distortion due to a lack of distinction between artificial and natural boundary features, and severe projection distortion and misconnection caused by conventional two-dimensional envelope algorithms in complex mountainous terrain.

[0007] To achieve the aforementioned objectives, the technical solution adopted by this invention is as follows: a method for calculating the forest encroachment area of ​​ski resorts based on adaptive envelope reconstruction, comprising: S1: Based on multi-resolution remote sensing image data and land use raster data, pixel feature extraction and adaptive threshold discrimination methods are used to obtain the forest binary raster data of the target area; S2: Using forest binary raster data and ski resort boundary vector data within the target area, spatial overlay processing is performed to extract a spatially constrained set of forest pixels located within the ski resort area; S3: Based on the forest pixel set, extract the spatial boundary and perform discrete point sampling. Use the boundary type discrimination and adaptive adjustment strategy based on local geometric consistency to obtain a unique boundary point set. S4: Based on the unique set of boundary points and digital elevation data, calculate the local terrain slope and three-dimensional spatial distance, and use the three-dimensional circle constraint method combined with slope correction to screen candidate points and construct an effective edge set; S5: Using the set of effective edges, construct a planar topological graph structure and perform closed path tracing to generate multiple candidate surface features. Spatial fusion and area filtering are performed on the candidate surface features, and topological repair processing is executed to obtain the outer envelope surface representing the continuous spatial distribution of the forest. S6: Based on the forest pixel set and its spatial resolution, calculate the actual forest area, and calculate the envelope area based on the outer envelope surface. Calculate the difference between the envelope area and the actual forest area to obtain the forest encroachment area of ​​the ski resort, thus completing the calculation of the forest encroachment area of ​​the ski resort.

[0008] Further, S1 includes: Acquire multi-resolution remote sensing image data and land use raster data; The land use raster data is classified and analyzed, and reflectance information of multiple bands is extracted from the multi-resolution remote sensing image data. Through calculation, the index features related to vegetation, water bodies and shortwave infrared suppression are obtained. Based on exponential features, a forest probability discrimination model is constructed, and an adaptive threshold is used to perform binary classification of each pixel to obtain the forest binary raster data of the target area.

[0009] Furthermore, the expression for the forest probability discriminant model is: ; ; ; ; ; in, This represents the normalized vegetation index at coordinates (x, y). This represents the water index at coordinates (x, y). This represents the shortwave infrared suppression index at coordinates (x, y). This represents the near-infrared reflectance extracted from multi-resolution remote sensing image data. Indicates the reflectivity in the red light band. Indicates the reflectivity in the green light band. Indicates the reflectivity in the shortwave infrared band. Let represent the forest discriminant response function at coordinates (x, y). , and These are the preset weighting coefficients corresponding to vegetation, water bodies, and shortwave infrared suppression features, respectively. Let represent the forest probability distribution function at coordinates (x, y). For the slope parameter, This is the bias coefficient.

[0010] Further, S3 includes: Based on the forest pixel set, the spatial boundary is extracted, and the discrete points of the spatial boundary are extracted as initial sampling points; Calculate the local curvature, rate of curvature change, and neighborhood orientation consistency at each initial sampling point. Dynamically adjust the basic sampling interval based on the local curvature. Construct a boundary discrimination model based on local geometric consistency to classify the spatial boundary into artificial boundaries and natural boundaries. For regions identified as artificial or natural boundaries, the corresponding sampling interval adjustment strategy is executed, and the generated sampling points are deduplicated to obtain a unique set of boundary points.

[0011] Furthermore, the expression for the boundary discrimination model is: ; in, This represents the artificial boundary confidence level of the i-th initial sampling point. This represents the local curvature at the i-th initial sampling point, calculated based on the spatial boundary. This represents the rate of change of curvature at the i-th initial sampling point, calculated based on the curvature difference between adjacent sampling points. This indicates the consistency of the neighborhood orientation at the i-th initial sampling point, calculated based on the angle between the orientations of adjacent sampling points. , and These are the corresponding preset weight coefficients.

[0012] Furthermore, the expression for the basic sampling interval is: ; in, This represents the dynamically adjusted sampling interval for the i-th initial sampling point. Indicates the preset base sampling interval. The preset curvature adjustment parameters, It represents the local curvature at the i-th initial sampling point.

[0013] Further, S4 includes: Obtain digital elevation data, and calculate the local terrain slope at each boundary point based on the unique set of boundary points and the digital elevation data; Based on the local terrain slope and the boundary points, an adaptive scale parameter with a slope correction factor and a three-dimensional spatial distance are constructed. The three-dimensional spatial distance and the adaptive scale parameter are used to filter out a set of three-dimensional candidate connections that meet the conditions. A three-dimensional circular constraint model based on the adaptive scale parameter is established on the local tangent plane, and a terrain consistency joint test is performed on the three-dimensional candidate connection set to construct an effective edge set.

[0014] Furthermore, the expression for the three-dimensional spatial distance and the adaptive scale parameter is as follows: ; ; ; ; in, Let v represent the three-dimensional spatial distance between the v-th unique boundary point and the u-th unique boundary point. and Let represent the plane coordinates and the elevation value obtained from digital elevation data for the v-th unique boundary point, respectively. and Let represent the plane coordinates and the elevation value obtained from digital elevation data for the u-th unique boundary point, respectively. This represents the slope influence factor at the v-th unique boundary point. This represents the local slope at the v-th unique boundary point calculated based on digital elevation data. The preset slope adjustment coefficient, This represents the adaptive scaling parameter at the v-th uniqueness boundary point. This indicates the preset basic scale of the artificial boundary. Let represent the combined scale of the candidate connections formed by the v-th unique boundary point and the u-th unique boundary point. This represents the function that takes the minimum value.

[0015] Furthermore, the three-dimensional circle constraint model is a candidate circle located on a local tangent plane, with its center... The expression is: ; ; The boundary equation of the candidate circle is: ; in, Let represent the planar distance between the v-th unique boundary point and the u-th unique boundary point. Let represent the coordinates of the midpoint of the line connecting the v-th unique boundary point and the u-th unique boundary point. Indicates the corresponding combination scale, Let represent the unit normal vector determined by the v-th unique boundary point and the u-th unique boundary point. This represents any point on the boundary of the candidate circle.

[0016] Further, S5 includes: A topological graph structure is constructed based on the set of valid edges. Unvisited edges are traversed, and path expansion is performed with the endpoint connection relationship of the edges and the spatial continuity of adjacent edges as constraints to generate a closed loop structure. The internal region of the closed loop structure is polygonized to generate multiple candidate surface features. Spatial fusion operation is performed on each candidate surface feature to eliminate local overlap, resulting in a fused surface set. Calculate the area of ​​each surface element in the fused surface set, retain the main surface element with the largest area and remove fragmented areas, and perform topological restoration processing on the main surface element, including self-intersection elimination, hole filling and fracture repair, to obtain the outer envelope surface representing the continuous spatial distribution of the forest.

[0017] The beneficial effects of this invention are as follows: This invention provides a method for calculating the forest area encroached upon by ski resorts based on adaptive envelope reconstruction. It acquires multi-resolution remote sensing images and land use raster data, uses pixel feature extraction and adaptive threshold discrimination methods to obtain binary forest data, and applies ski resort boundary constraints. This mechanism effectively overcomes the extreme background interference caused by the interplay of high-albedo snow cover and tree shadows at ski resorts, accurately locates the target study area, and greatly improves the anti-interference capability and reliability of basic forest information extraction.

[0018] During discrete point sampling, a boundary type discrimination and adaptive adjustment strategy based on local geometric consistency is employed. This method intelligently distinguishes between straight edges created by bulldozing and naturally growing, randomly undulating edges, and dynamically adjusts the sampling density accordingly. This ensures the accuracy of morphological representation of complex natural boundaries while reducing redundant data computation for artificially smoothed boundaries, significantly improving the accuracy and efficiency of subsequent reconstruction.

[0019] This paper creatively introduces digital elevation data to calculate local terrain slope and uses a three-dimensional circular constraint method combined with slope correction for candidate point selection. This mechanism completely breaks through the limitations of conventional two-dimensional algorithms in mountainous applications. It not only compensates for the physical distance and area compression errors caused by slope projection, but also effectively blocks erroneous connections across mountain cliffs through spatial domain detection, ensuring that the reconstructed skeleton network perfectly matches the real undulating terrain features.

[0020] By constructing a topological graph for closed-path tracing and employing spatial fusion, area maximization filtering, and multiple topological restoration techniques, a smooth, continuous outer envelope surface free of internal fragments can be automatically reconstructed from forest patches artificially fragmented by ski runs. Finally, by comparing the area of ​​this envelope with the actual remaining area, an objective and precise quantification of the forest area encroached upon by ski resorts is achieved, providing a highly robust measurement tool for ecological damage assessment. Attached Figure Description

[0021] This specification will be further described by way of exemplary embodiments, which will be described in detail with reference to the accompanying drawings. These embodiments are not limiting; in these embodiments, the same reference numerals denote the same structures, wherein: Figure 1 This is an exemplary flowchart illustrating a method for calculating the forest encroachment area of ​​ski resorts based on adaptive envelope reconstruction, according to some embodiments of this specification. Detailed Implementation

[0022] The specific embodiments of the present invention are described below to enable those skilled in the art to understand the present invention. However, it should be understood that the present invention is not limited to the scope of the specific embodiments. For those skilled in the art, various changes are obvious as long as they are within the spirit and scope of the present invention as defined and determined by the appended claims. All inventions utilizing the concept of the present invention are protected.

[0023] Example Figure 1 This is an exemplary flowchart illustrating a method for calculating forest encroachment area by ski resorts based on adaptive envelope reconstruction, according to some embodiments of this specification. Figure 1 As shown, the process includes the following steps. In some embodiments, the process may be executed by a processor.

[0024] S1: Based on multi-resolution remote sensing image data and land use raster data, pixel feature extraction and adaptive threshold discrimination methods are used to obtain the forest binary raster data of the target area.

[0025] Multi-resolution remote sensing image data refers to a dataset of digital images of the Earth acquired by sensors carried by satellites or aircraft, containing different spatial resolutions and multiple spectral bands. For example, multi-resolution remote sensing image data can include reflectance images in visible light bands, such as red, green, near-infrared, and short-wave infrared bands. These image arrays in different bands are used to reflect the photosynthetic characteristics of surface vegetation, water absorption characteristics, and the physical state of the Earth's surface.

[0026] In some embodiments, the processor can obtain multi-resolution remote sensing image data by receiving raw observation signals from remote sensing satellites. Specifically, after receiving the raw multispectral image file, the system first performs radiometric calibration to convert the digital quantization values ​​recorded by the sensor into radiance. Then, atmospheric correction is performed to eliminate the interference of atmospheric scattering and absorption on spectral reflectance, and geometric correction is performed to eliminate spatial distortion. Finally, multidimensional matrix data with consistent spatial location and spectral reflectance for each band is extracted.

[0027] Land use raster data refers to a two-dimensional spatial dataset that divides the land surface into regular grids, or pixels, and assigns a specific land cover type attribute value to each grid. For example, land use raster data may include digital codes representing different land cover classifications such as forests, grasslands, water bodies, built-up land, or bare land. These discrete digital codes are used to visually display the current spatial distribution and category boundaries of land cover within a target area.

[0028] In some embodiments, the processor can obtain land use raster data by parsing and resampling an existing basic geographic information database. Specifically, the system reads a standard format land cover vector or raster classification file, performs spatial reprojection and raster alignment on the original file according to the set target spatial resolution, and then extracts a pixel matrix containing classification attribute codes to obtain a basic land cover classification base map for subsequent classification extraction and masking processing.

[0029] Forest binary raster data refers to a simplified spatial distribution matrix containing only two types of pixel values, used to clearly distinguish between forest and non-forest areas within a target region. For example, forest binary raster data can include a set of pixels representing the forest category (represented by the value "1") and a set of pixels representing the non-forest category (such as artificial ski slopes, bare land, and buildings) (represented by the value "0"), thus forming a black-and-white distribution layer that retains only the target vegetation information.

[0030] In some embodiments, the processor can obtain forest binary raster data by thresholding the input land cover classification results or probability model results. Specifically, the system traverses each cell position of the input raster layer and determines whether the category attribute identifier or the calculated forest probability value of the cell meets the preset judgment criteria. If the judgment criteria are met, the cell value at that position is assigned a value of 1; if the criteria are not met, it is assigned a value of 0, thereby generating a uniformly expressed binary forest spatial distribution layer.

[0031] In some embodiments, the processor can acquire multi-resolution remote sensing image data and land use raster data; perform category analysis on the land use raster data; extract reflectance information of multiple bands from the multi-resolution remote sensing image data; and calculate the exponential features related to vegetation, water bodies, and shortwave infrared suppression. Based on the exponential features, a forest probability discrimination model is constructed, and an adaptive threshold is used to perform binarization classification on each pixel to obtain the forest binary raster data of the target area.

[0032] In some embodiments, the forest binary raster data of the target region is obtained by using the pixel feature extraction and adaptive threshold discrimination method, and its adaptive threshold classification expression is: ; in, Representing coordinates The forest binary results are used to construct forest binary raster data. Representing coordinates The probability distribution function of the forest at that location. The threshold is for adaptive determination.

[0033] In some embodiments, the spatial overlay processing based on the forest binary raster data and the ski resort boundary vector data is used to extract a spatially constrained set of forest pixels located within the ski resort area. The overlay and extraction expression is as follows: ; ; in, This represents a constrained cell layer after spatial overlay. This represents a unified forest distribution based on binary forest raster data. This indicates the mask range based on the ski resort boundary vector data; This represents the spatially constrained set of forest pixels extracted. In the set of forest pixels, the first... The spatial coordinates of each pixel.

[0034] Exponential features refer to a set of numerical indicators that utilize the mathematical calculations between reflectance in different spectral bands of remote sensing images to enhance information about specific ground features and suppress background interference. For example, exponential features may include the Normalized Difference Vegetation Index (NDVI) to highlight vegetation growth status and coverage, the Water Body Index to identify water bodies or humidity distribution, and the Shortwave Infrared Suppression Index specifically designed to suppress high-albedo backgrounds, such as snow cover interference. Together, they constitute a multidimensional feature vector describing the physical characteristics of a single pixel.

[0035] In some embodiments, the processor can obtain exponential features by performing band algebraic operations on the reflectance of each band of a multispectral remote sensing image pixel. Specifically, the system extracts the reflectance values ​​of the near-infrared, red, green, and short-wave infrared bands of the currently processed pixel. By calculating the ratio of the difference between the reflectances of different bands to the sum of the reflectance values ​​of the corresponding bands, the system calculates the vegetation index, water index, and short-wave infrared suppression index values ​​of the pixel, thereby converting the original reflectance into a feature vector for model discrimination.

[0036] A forest probability discrimination model is a mathematical evaluation system that maps the multidimensional feature vector of a pixel to the probability that the pixel belongs to a real forest category. For example, a forest probability discrimination model may include a continuous logistic distribution function whose output value is a continuous real number between 0 and 1. The closer the value is to 1, the higher the confidence that the pixel location is a natural forest, and the closer the value is to 0, the higher the probability that the location is an artificial ski slope or other non-forest land.

[0037] In some embodiments, the processor can obtain a forest probability discrimination model by combining multiple exponential features, performing a weighted summation, and applying a nonlinear function transformation. The system first obtains the values ​​of each exponential feature of the current pixel, assigns preset weight coefficients to each feature, and performs linear summation to calculate a comprehensive response value. Subsequently, the comprehensive response value is substituted into a smooth activation function containing slope and bias parameters for nonlinear mapping, and finally the probability distribution value of the pixel belonging to the forest category is calculated.

[0038] In some embodiments, the expression for the forest probability discriminant model is: ; ; ; ; ; in, This represents the normalized vegetation index at coordinates (x, y). This represents the water index at coordinates (x, y). This represents the shortwave infrared suppression index at coordinates (x, y). This represents the near-infrared reflectance extracted from multi-resolution remote sensing image data. Indicates the reflectivity in the red light band. Indicates the reflectivity in the green light band. Indicates the reflectivity in the shortwave infrared band. Let represent the forest discriminant response function at coordinates (x, y). , and These are the preset weighting coefficients corresponding to vegetation, water bodies, and shortwave infrared suppression features, respectively. Let represent the forest probability distribution function at coordinates (x, y). For the slope parameter, This is the bias coefficient.

[0039] In some preferred embodiments, for the extreme case of high-albedo snow background and tree shadows intertwined at ski resorts, the present invention specifically sets a specific weighting coefficient range when constructing the forest discriminant response function, that is, sets the weighting coefficient of the vegetation index. The range of values ​​is Weighting coefficients of water index The range of values ​​is Weighting coefficients of the shortwave infrared suppression index The range of values ​​is These parameter ranges are optimal anti-interference empirical thresholds determined through extensive cross-comparison experiments. They can effectively suppress misjudgments caused by snow reflection while ensuring the dominant position of vegetation signals, and accurately extract weak vegetation signals in complex and interwoven backgrounds.

[0040] S2: Using forest binary raster data and ski resort boundary vector data within the target area, spatial overlay processing is performed to extract a spatially constrained set of forest pixels located within the ski resort area.

[0041] Ski resort boundary vector data refers to polygonal geometric entities composed of a series of nodes with spatial coordinates and line segments connecting these nodes, used to precisely define the legal spatial scope for the construction and operation of a ski resort. For example, ski resort boundary vector data may include a series of latitude and longitude coordinates or projected coordinates of the ski resort's outer contour, as well as topological relationship records that form a closed polygon, used to define the absolute operating area for subsequent spatial analysis in a geographic information system.

[0042] In some embodiments, the processor can obtain ski resort boundary vector data by receiving an externally imported mapping boundary file and parsing its geometric node coordinates. Specifically, the system reads a standard format vector shape file, such as Shapefile format, parses the geometric attribute table in the file, extracts the ordered coordinate sequence that constitutes the outer ring of the polygon, and converts it into the unified spatial reference coordinate system of the current system, thereby obtaining the polygon vector boundary used for spatial range constraints.

[0043] A forest cell set refers to a list of coordinates of discrete spatial grid cells that simultaneously satisfy the conditions of being located within a specific spatial range and having their land cover attribute determined to be forest. For example, a forest cell set may include tens of thousands of two-dimensional planar coordinate pairs (x, y), each representing a specific spatial location located within the planned boundary of a ski resort and having a cell value of "1" in the binary raster data of the forest. This set constitutes the basic data source for subsequent boundary extraction.

[0044] In some embodiments, the processor can obtain a set of forest pixels by performing a logical intersection operation between the forest binary raster data and the spatial boundary. Specifically, the system converts the ski resort boundary vector data into a mask layer with the same resolution as the raster data, then spatially overlays and compares this mask layer with the forest binary raster data, filters out pixels that are simultaneously located within the effective area of ​​the mask and have a raster value of 1, extracts the center point coordinates of these pixels, and combines them to form a spatially constrained set of pixels.

[0045] S3: Based on the forest pixel set, extract the spatial boundary and perform discrete point sampling. Utilize the boundary type discrimination and adaptive adjustment strategy based on local geometric consistency to obtain a unique set of boundary points.

[0046] A unique set of boundary points refers to a series of discrete coordinate points that constitute the external geometric outline of a forest area and are not spatially repetitive or extremely clustered. For example, a unique set of boundary points may include a series of two-dimensional coordinate points distributed along the edge of a forest patch. After adjusting the sampling interval, these points maintain a reasonable spatial distance from each other, which can accurately delineate the undulating shape of the forest and eliminate overlapping coordinate data caused by computational redundancy.

[0047] In some embodiments, the processor can obtain a unique set of boundary points by performing distance deduplication and filtering operations on the extracted initial boundary sampling points. Specifically, the system first obtains a list of initial sampling points generated based on adaptive spacing, then calculates the spatial Euclidean distance between any two points in the list. When the distance between two points is found to be less than a preset minimum tolerance threshold, it is determined to be a duplicate point and one of them is deleted. After traversing the entire list to complete the removal, the final set of boundary point coordinates without duplicate positions is output.

[0048] In some embodiments, the processor can extract spatial boundaries based on the forest pixel set, extract discrete points of the spatial boundaries as initial sampling points; calculate the local curvature, rate of change of curvature, and neighborhood orientation consistency at each initial sampling point; dynamically adjust the basic sampling interval based on the local curvature; and construct a boundary discrimination model based on local geometric consistency to classify the spatial boundaries into artificial boundaries and natural boundaries; for regions determined to be artificial boundaries and natural boundaries, execute the corresponding sampling interval adjustment strategy, and perform deduplication processing on the generated sampling points to obtain a unique set of boundary points.

[0049] Initial sampling points refer to discrete spatial coordinate positions that are initially extracted along the extracted forest spatial boundary line according to certain initial rules to describe the boundary morphology. For example, initial sampling points may include a coordinate point extracted at fixed pixel intervals along the grid edge. These points are densely distributed along the boundary line between the forest area and the non-forest area, serving as the reference input data for subsequent curvature calculation and adaptive density adjustment.

[0050] In some embodiments, the processor can obtain initial sampling points by taking points along the edges of a binary image at fixed steps using an edge tracing algorithm. Specifically, the system performs contour tracing on a spatially constrained forest binary raster region to obtain the closed boundary pixel chain code; then, according to a preset basic sampling interval, for example, half the general width of a ski slope, the corresponding curve length is calculated along the pixel chain code, and the spatial coordinates of the current pixel are recorded every time an interval length is reached, thereby generating an initial discrete sampling point sequence.

[0051] Local curvature refers to a geometric property parameter that describes the degree of bending of a spatial boundary near a specific sampling point. For example, local curvature can include a specific numerical value. The larger the value, the sharper the bend and the more turbulent the boundary is at that point, which is commonly seen at the edges of naturally growing forests. The smaller the value, or even close to zero, the closer the boundary is to a straight line at that point, which is commonly seen at the edges of artificial ski slopes leveled by bulldozers.

[0052] In some embodiments, the processor can obtain the local curvature by performing circle fitting or difference calculation using the spatial coordinates of the current sampling point and several adjacent sampling points. Specifically, the system selects the current sampling point and two adjacent points to form a local point set, calculates the radius of the circle determined by these three points, and uses the reciprocal of this radius as the curvature value of the current point; or, it calculates the rate of change of the direction angle of adjacent line segments, thereby quantifying a numerical index describing the local bending characteristics of the boundary.

[0053] The rate of change of curvature is a second-order geometric parameter that describes the degree of curvature change of a spatial boundary as it extends along a specific direction; that is, the speed or abruptness of the change in curvature. For example, the rate of change of curvature can include a gradient value, which increases significantly when a boundary abruptly changes from a straight line to a sharp bend; at the edge of a natural forest belt, the value usually remains at a high level due to the high-frequency random fluctuations in curvature; and at artificially leveled edges, the value is relatively stable and low.

[0054] In some embodiments, the processor can obtain the rate of curvature change by calculating the gradient difference between the local curvatures of adjacent sampling points. Specifically, after calculating the local curvatures of all initial sampling points, the system extracts the curvature values ​​of the preceding and succeeding nodes for the currently processed sampling point, calculates the absolute value of the difference between the curvature of the current point and the curvature of the adjacent points, and divides it by the arc length or spacing between them to obtain the rate of change value reflecting the severity of curvature fluctuations.

[0055] Neighborhood orientation consistency refers to a structural characteristic index describing whether the local boundary line segment where a sampling point is located maintains uniformity or smoothness in its extension direction. For example, neighborhood orientation consistency can include a correlation coefficient or angle cosine value between 0 and 1. A higher value indicates that the boundary before and after the point is almost on the same straight line, showing strong artificial cutting marks; a lower value indicates that the boundary direction has undergone obvious deflection or irregular distortion at that point.

[0056] In some embodiments, the processor can obtain neighborhood direction consistency by calculating and normalizing the angle between adjacent boundary segments. Specifically, the system takes the current sampling point as the vertex and connects its predecessor and successor nodes to form two spatial vectors; then, it uses the vector dot product formula to calculate the spatial angle between these two vectors and calculates the cosine value or other normalization function value of the angle; the calculation result is used as an evaluation index to quantify whether the local boundary extension direction is consistent.

[0057] The base sampling interval refers to a default, uniform spatial distance span set before considering the specific shape and geometric features of the boundary for extracting discrete points along the forest boundary. For example, the base sampling interval can include a specific physical length value, such as 30 meters. This value is usually set proportionally based on the common characteristic width of artificial disturbances in the target area, such as ski slopes, and serves as the initial reference baseline for subsequent adaptive density adjustment algorithms.

[0058] In some embodiments, the processor can obtain the basic sampling interval by reading a preset configuration file or a manually input empirical standard value. Specifically, the system receives external parameter input, which is usually derived from the average width statistics of the ski slopes in the local ski resort; the system multiplies the obtained ski slope width value by a fixed scaling factor, such as half the width, and the result is used by the system as the fixed step size or basic distance when initially extracting boundary points.

[0059] A boundary discrimination model is a comprehensive classification and evaluation function used to distinguish whether each sampling point on a spatial boundary belongs to a naturally grown forest edge or an artificial edge formed by logging. For example, a boundary discrimination model may include a confidence calculation formula that incorporates local geometric features. When the output confidence score is higher than a set threshold, the model determines that the current region has low-frequency smoothness and strong directional consistency characteristics of logging; otherwise, it determines that it has natural fluctuation characteristics.

[0060] In some embodiments, the processor can obtain a boundary discrimination model by performing a linear combination operation on the local curvature, rate of change of curvature, and neighborhood orientation consistency of the sampling point. Specifically, the system obtains the curvature-related index and orientation consistency index of the current sampling point, multiplies these geometric indices or their transformed values ​​by pre-set weighting coefficients, and then sums up the product terms to calculate a comprehensive judgment score reflecting the probability of the artificial boundary attribute of the point.

[0061] In some embodiments, the boundary discrimination model based on local geometric consistency is based on a qualitative analysis of the physical process of ski resort construction: natural forest belts, constrained by terrain undulations and lighting conditions, inevitably exhibit high-frequency random fluctuations and frequent directional changes at their natural growth edges; while artificial ski slopes are the result of heavy bulldozing operations, and their excavated edges inevitably exhibit low-frequency smoothness and strong directional consistency. Therefore, this invention, by linearly combining three geometric parameters—local curvature, rate of curvature change, and neighborhood directional consistency—can accurately quantify and identify the smooth straight lines of artificial cutting and the tortuous edges of natural growth.

[0062] In some embodiments, the expression for the boundary discrimination model is: ; in, This represents the artificial boundary confidence level of the i-th initial sampling point. This represents the local curvature at the i-th initial sampling point, calculated based on the spatial boundary. This represents the rate of change of curvature at the i-th initial sampling point, calculated based on the curvature difference between adjacent sampling points. This indicates the consistency of the neighborhood orientation at the i-th initial sampling point, calculated based on the angle between the orientations of adjacent sampling points. , and These are the corresponding preset weight coefficients.

[0063] In some embodiments, the expression for the basic sampling interval is: ; in, This represents the dynamically adjusted sampling interval for the i-th initial sampling point. Indicates the preset base sampling interval. The preset curvature adjustment parameters, It represents the local curvature at the i-th initial sampling point.

[0064] Spatial boundaries refer to the geometric outlines that divide forested and non-forested areas in geographical space. Based on their formation, they can be divided into artificial boundaries and natural boundaries. For example, artificial boundaries can include those formed by bulldozing operations such as ski slope construction and road excavation, exhibiting a low-curvature, smooth, and uniformly oriented straight or regular arc-shaped outline. Natural boundaries, on the other hand, can include forest edges formed by natural conditions such as topography and sunlight, exhibiting high-frequency random fluctuations, meandering, and irregular shapes.

[0065] In some embodiments, the processor can obtain spatial boundaries by identifying locations where pixel attributes in the forest binary raster layer undergo abrupt changes and concatenating them. Specifically, the system traverses the range-constrained forest binary raster layer, using edge detection operators, such as differential checks based on eight or four neighborhoods, to find forest pixels surrounded by non-forest pixels. Subsequently, the outer edge segments of these boundary pixels are extracted, and the segments are connected end-to-end according to topological adjacency to construct a continuous geometric contour that distinguishes different land cover categories.

[0066] In some embodiments, the discrimination threshold used when classifying spatial boundaries as artificial or natural boundaries It can be dynamically determined based on the statistical distribution of boundary features, and its adaptive expression is: ; in, Indicates the discrimination threshold. This represents the mean of the artificial boundary confidence scores corresponding to all initial sampling points. Indicates the corresponding standard deviation. This indicates the preset threshold adjustment parameter.

[0067] S4: Based on the unique set of boundary points and digital elevation data, calculate the local terrain slope and three-dimensional spatial distance, and use the three-dimensional circle constraint method combined with slope correction to screen candidate points and construct an effective edge set.

[0068] Local terrain slope is a quantitative angular index of the degree of inclination of the ground surface at a specific sampling point, reflecting the rate of elevation change of the terrain at that point. For example, local terrain slope can include an angle value between 0 degrees and 90 degrees, where a value of 0 degrees represents that the terrain is completely flat, and a value of 45 degrees represents that the location is on an extremely steep mountain slope. This value is used to counteract the distortion of the two-dimensional planar projected area caused by mountainous terrain in subsequent calculations.

[0069] In some embodiments, the processor can obtain the local terrain slope by performing gradient calculations on the neighborhood elevation values ​​around the sampling point in the input digital elevation model (DEM) data. Specifically, the system extracts the elevation values ​​of the point and its surrounding adjacent grids (such as a 3×3 window) from the DEM data based on the planar coordinates of the boundary point; then, using a third-order inverse distance squared weighted difference or related terrain analysis algorithm, it calculates the elevation change rate of the point in the horizontal and vertical directions, and finally solves for the tilt angle corresponding to the maximum change rate as the slope value.

[0070] Three-dimensional spatial distance refers to the actual straight-line segment length between two spatial coordinate points in a three-dimensional Cartesian coordinate system, taking into account variations in terrain elevation. For example, three-dimensional spatial distance can include the comprehensive Euclidean distance calculated based on the coordinate differences of the three dimensions of longitude, latitude, and altitude between two boundary points. In steep mountain ski resort environments, this distance is used to replace the traditional two-dimensional planar projection distance, effectively preventing severe compression or distortion of boundary segment lengths due to terrain inclination, and ensuring the accuracy of the reconstructed scale.

[0071] In some embodiments, the processor can calculate the three-dimensional spatial distance by combining planar coordinates and digital elevation model data to perform vector magnitude calculation. Specifically, the system obtains the horizontal projected coordinates of two boundary points and extracts the corresponding elevation values ​​of these two points from the digital elevation model; then, it calculates the squared difference of the coordinates of the two points on the two orthogonal axes in the horizontal direction and the squared difference of the elevation in the vertical direction; after adding these three squared terms, it performs a square root operation to finally obtain the straight-line spatial distance that reflects the actual terrain undulations.

[0072] An effective edge set refers to a combination of reliable spatial line segments that, after a series of rigorous scale screenings, spatial domain detections, and terrain constraint determinations, are confirmed by the system to accurately reflect the continuous outline of the forest's exterior. For example, an effective edge set can include hundreds or thousands of three-dimensional vector line segments that can be connected end to end. These line segments exclude erroneous connections that cross cliffs and redundant connections that penetrate the existing forest interior, forming the most basic geometric element dataset used to subsequently stitch together and generate the skeleton of a closed polygon entity.

[0073] In some embodiments, the processor can obtain a set of valid edges by applying a three-dimensional circle constraint model to candidate connection point pairs and combining it with terrain consistency checks. Specifically, the system first places two boundary point pairs that satisfy the three-dimensional spatial distance threshold into a candidate list; then, it constructs a three-dimensional candidate circle on a local tangent plane and performs an empty circle detection to determine whether the circle contains other boundary points; for the connection that passes the empty circle detection, it further checks whether the slope gradient difference between its two endpoints is less than a set breakage compensation threshold; finally, it retains the connection segments that simultaneously satisfy both spatial and terrain constraints to form the set.

[0074] In some embodiments, the processor can acquire digital elevation data, calculate the local slope at each boundary point based on the unique set of boundary points and the digital elevation data; construct an adaptive scale parameter with a slope correction factor and a three-dimensional spatial distance based on the local slope and the boundary points; use the three-dimensional spatial distance and the adaptive scale parameter to filter out a set of three-dimensional candidate connections that meet the conditions; establish a three-dimensional circular constraint model based on the adaptive scale parameter on the local tangent plane; perform a joint terrain consistency check on the set of three-dimensional candidate connections to construct a set of valid edges.

[0075] Digital elevation data (DEM) refers to a spatial digital model matrix distributed in a regular array to record continuous topographic relief and elevation information of the Earth's surface. For example, DEM can include a raster file format containing the absolute sea-level elevation values ​​corresponding to the center points of each grid. Its accuracy is typically at the meter or sub-meter level. These elevation values ​​can intuitively reflect the three-dimensional topographic features of the mountains, such as ridges, valleys, and slopes, where ski resorts are located, and are the core basic data for slope analysis and three-dimensional distance correction.

[0076] In some embodiments, the processor can obtain digital elevation data by parsing point cloud results files generated by aerial photogrammetry or airborne lidar scanning. Specifically, the system receives a high-precision surface point cloud data source from an external source, uses morphological filtering algorithms to filter out non-terrain surface points such as vegetation canopies and buildings on the ground, and retains only the true bare ground elevation points; subsequently, it uses spatial interpolation algorithms such as inverse distance weighted or kriging interpolation to resample and reconstruct these discrete elevation points into a regular grid matrix, thereby outputting a standard digital elevation model layer.

[0077] Boundary points refer to the basic discrete spatial coordinate nodes that constitute the continuous spatial outline of the forest. They are the lowest-level vertex elements used in the envelope reconstruction algorithm for topological connection and polygon generation. For example, boundary points can include three-dimensional spatial nodes that have both two-dimensional planar projection coordinates (x, y) and corresponding terrain elevation attributes (z), which are retained after the initial deduplication process. These points are densely distributed on the boundary line between forest patches and external open areas or ski slopes, and their density and spatial location directly determine the shape of the final reconstructed envelope surface.

[0078] In some embodiments, the processor can obtain boundary points with three-dimensional information by performing elevation attribute assignment operations on a unique set of boundary points. Specifically, the system extracts the coordinate list of two-dimensional unique boundary points after adaptive sampling and duplicate removal; then, it performs spatial location matching of these two-dimensional coordinates with digital elevation data, and uses a bilinear interpolation algorithm to extract the terrain elevation of each point's location; finally, it appends the extracted elevation values ​​to the original two-dimensional coordinates to form node data containing complete three-dimensional coordinate information.

[0079] An adaptive scaling parameter is a variable spatial distance threshold used to dynamically control the search range of connection points and determine the size of the geometric constraint circle when performing a boundary envelope reconstruction algorithm. For example, the adaptive scaling parameter may include a dynamically corrected value that combines the general physical width of the local ski slope with the local terrain slope of the boundary point. When the point is in a flat area, this parameter is close to the basic set width; when the point is on a steep slope, the parameter will be adaptively reduced according to the slope influence factor to counteract the deformation error caused by the 3D spatial projection.

[0080] In some embodiments, the processor can obtain adaptive scale parameters by combining the artificial boundary base scale with the terrain slope extracted based on elevation data through a function correction. Specifically, the system first obtains a preset artificial boundary base scale value; simultaneously, it calculates a slope influence factor greater than or equal to one based on the local slope at the boundary point location; then, it divides the base scale value by the slope influence factor to calculate a specific scale value corrected for the current local terrain; for a candidate connection consisting of two boundary points, the smaller of the corrected scale values ​​of the two points is taken as the combined scale parameter of the connection.

[0081] In some embodiments, the expression for the three-dimensional spatial distance and the adaptive scale parameter is: ; ; ; ; in, Let v represent the three-dimensional spatial distance between the v-th unique boundary point and the u-th unique boundary point. and Let represent the plane coordinates and the elevation value obtained from digital elevation data for the v-th unique boundary point, respectively. and Let represent the plane coordinates and the elevation value obtained from digital elevation data for the u-th unique boundary point, respectively. This represents the slope influence factor at the v-th unique boundary point. This represents the local slope at the v-th unique boundary point calculated based on digital elevation data. The preset slope adjustment coefficient, This represents the adaptive scaling parameter at the v-th uniqueness boundary point. This indicates the preset basic scale of the artificial boundary. Let represent the combined scale of the candidate connections formed by the v-th unique boundary point and the u-th unique boundary point. This represents the function that takes the minimum value.

[0082] In some embodiments, given the extreme elevation differences in mountain ski resorts, the projected area of ​​forests on the slope would be severely compressed if a conventional two-dimensional envelope algorithm were used. Therefore, this invention employs a slope influence factor to counteract mountain distortion as a specifically designed terrain compensation mechanism. Furthermore, a fault compensation mechanism for extreme terrain is introduced: calculating the slope gradient difference between two adjacent candidate connection boundary points. ,when Greater than a certain threshold (preferably) When an error occurs, the system will identify the area as a region of abrupt elevation change or a cliff, and forcibly sever the candidate connections established based on two-dimensional Euclidean distance. This error prevention effect greatly improves the topological robustness of the 3D reconstruction model in complex mountainous terrain.

[0083] A 3D candidate connection set refers to a list of potentially valid boundary point pairs initially selected during the construction of the topology graph, where the 3D spatial distance between two points satisfies specific scale constraints. For example, the 3D candidate connection set may include a large number of unverified pairwise vertex indices and the line segments they form. The length of these line segments is less than or equal to twice the combined scale parameter dynamically calculated by the algorithm. They are temporarily stored as alternatives that may constitute the outer envelope contour, awaiting subsequent more stringent spatial exclusion checks.

[0084] In some embodiments, the processor can obtain a set of 3D candidate connections by traversing the list of boundary points and performing a 3D spatial distance comparison. Specifically, the system performs pairwise traversal of the set of boundary points containing 3D coordinate information, calculating the 3D spatial distance between any two boundary points; then, it calculates the adaptive combined scale parameter corresponding to these two points; it determines whether the obtained 3D spatial distance is less than twice the combined scale parameter; if the distance threshold condition is met, the two boundary points are recorded as a point pair and included in the candidate set.

[0085] The 3D circle constraint model is a spatial geometric criterion derived from the traditional 2D planar Alpha Shape algorithm, used to verify whether the line connecting two boundary points in 3D undulating terrain constitutes the true outer contour. For example, the 3D circle constraint model can include a virtual disk structure constructed by projecting along a local terrain tangent plane, the radius of which is determined by dynamically calculated adaptive combined scale parameters. The model determines whether the detected line segment is the true outermost boundary by detecting whether the interior of this virtual disk is completely empty, i.e., it does not contain other forest boundary points.

[0086] In some embodiments, the processor can construct a three-dimensional circle constraint model by combining spatial analytical geometry calculations and the spatial relationships of candidate point pairs. Specifically, the system extracts the three-dimensional coordinates of candidate point pairs, calculates the midpoint of the line connecting these two points and the unit normal vector of the line segment formed by these two points on the tangent plane; then, using a combination of scale parameters and the planar projection distance between the two points, the offset of the detection circle center in the direction of the normal vector is calculated using the Pythagorean theorem, thereby determining the precise three-dimensional coordinates of the virtual circle center; finally, a mathematical logic equation is established using this circle center and scale parameters to detect whether other points fall within the circle.

[0087] In some embodiments, the three-dimensional circle constraint model is a candidate circle located on a local tangent plane, with its center... The expression is: ; ; The boundary equation of the candidate circle is: ; in, Let represent the planar distance between the v-th unique boundary point and the u-th unique boundary point. Let represent the coordinates of the midpoint of the line connecting the v-th unique boundary point and the u-th unique boundary point. Indicates the corresponding combination scale, Let represent the unit normal vector determined by the v-th unique boundary point and the u-th unique boundary point. This represents any point on the boundary of the candidate circle.

[0088] In some embodiments, performing spatial domain detection based on the three-dimensional circle constraint model includes: determining whether the candidate circle contains other boundary points besides the v-th unique boundary point and the u-th unique boundary point; if not, retaining the candidate connection; and constructing a set of valid edges after combining terrain consistency filtering.

[0089] S5: Using the set of effective edges, construct a planar topological graph structure and perform closed path tracing to generate multiple candidate surface features. Spatial fusion and area filtering are performed on the candidate surface features, and topological repair processing is executed to obtain the outer envelope surface representing the continuous spatial distribution of the forest.

[0090] Candidate surface features refer to preliminary two-dimensional or three-dimensional polygonal spatial entities enclosed by a set of continuous, closed, and valid edge sequences after closed path tracing. For example, candidate surface features may include multiple closed vector graphic blocks of varying sizes and shapes. Among these initially generated blocks, some may represent large, contiguous main forest areas, while others may only be tiny fragmented surfaces composed of algorithmic noise or a very small number of scattered trees. Some blocks may even partially overlap or have tangent edges in space.

[0091] In some embodiments, the processor can obtain candidate face features by performing a polygonal filling operation on the internal spatial region of the generated closed loop structure. Specifically, the system receives a closed loop chain code or node sequence consisting of ordered valid edges; first, it uses a polygon generation algorithm, such as ray casting or parity checking, to define the internal two-dimensional planar region enclosed by the closed loop; then, it instantiates the internal planar region into a solid face feature with a defined geometric area, perimeter, and topological inclusion relationship; this operation is performed on all closed loops in sequence, outputting a set consisting of multiple initial polygons.

[0092] In some embodiments, the expression for generating multiple candidate surface features and performing spatial fusion and area filtering is as follows: ; ; ; in, Indicates by the first The first closed-loop structure generates the first Each candidate facet element This represents a function that seeks the internal spatial region enclosed by a closed loop. Represents the set of merged surfaces. This represents the face feature entities in the merged face set. Represents the union and fusion operations of spatial representations. This indicates the total number of candidate face features. This represents the initial state of the main surface feature with the largest area retained after screening, i.e., the outer envelope surface. A function that calculates the geometric area. This represents the variable function that takes the maximum area.

[0093] The outer envelope surface refers to the continuous and complete outer spatial distribution boundary polygon of the forest within the target area, assuming no human disturbance, such as the creation of ski slopes. For example, the outer envelope surface can include a single, huge vector surface structure that has undergone comprehensive topological restoration. This surface structure not only smoothly spans the internal gaps created by ski slope cutting, eliminating tiny holes and cracks that may appear during the algorithm reconstruction process, but also excludes minor area interference from non-subject elements. It is the core basis for ultimately calculating the theoretical forest area.

[0094] In some embodiments, the processor can obtain the outer envelope surface by performing area filtering and a series of spatial topology repair processes on the fused surface set. Specifically, after calculating the geometric area of ​​each fused surface element, the system sorts and compares the polygons, retaining only the one with the largest area value as the main contour surface, and removing all other fragmented surfaces with smaller areas. Subsequently, the system calls a geographic information topology repair tool to automatically identify and eliminate self-intersecting lines inside the main contour surface, fill in isolated non-data hole areas inside, and connect and repair minor boundary breaks, ultimately outputting a perfectly closed solid surface.

[0095] In some embodiments, the processor can construct a topological graph structure based on the set of valid edges, traverse unvisited edges, extend paths with the endpoint connections of edges and the spatial continuity of adjacent edges as constraints, and generate a closed loop structure; perform polygonization of the internal region of the closed loop structure to generate multiple candidate surface features, perform spatial fusion operation on each candidate surface feature to eliminate local overlap, and obtain a fused surface set; calculate the area of ​​each surface feature in the fused surface set, retain the main surface feature with the largest area and remove fragmented areas, and perform topological repair processing such as self-intersection elimination, hole filling and fracture repair on the main surface feature to obtain the outer envelope surface representing the continuous spatial distribution of the forest.

[0096] A topological graph structure is a mathematical graph theory model used to describe non-metric spatial relationships such as adjacency, connectivity, and containment between spatial data. For example, a topological graph structure can include an adjacency matrix or adjacency list structure, where the set of boundary points is the node in the graph, and the set of valid edges that have been verified and retained is the arc segment connecting adjacent nodes. This structure does not focus on the specific length and curvature of the line segments; it only records which point is connected to which point by which edge, providing an efficient underlying architecture for subsequent path search of closed polygon boundaries.

[0097] In some embodiments, the processor can obtain the topological graph structure by extracting the endpoint connections of geometric features and constructing a relational database table. Specifically, the system traverses all valid edge features, extracts the start and end coordinates of each valid edge; merges endpoints with the same coordinates into a single unique node; then, it establishes a node-edge association mapping table, recording which edges each node connects to, and which two nodes each edge connects to, thereby generating abstract network topological data in memory to support graph traversal algorithms and path connectivity analysis.

[0098] In some embodiments, the expression for constructing the topological graph structure and generating the closed loop structure based on the set of valid edges is: ; ; in, This represents the structure of the constructed topology graph. This represents the set of unique boundary points used to construct graph nodes. This represents the set of extracted valid edges. Indicates the generation of the tracking number A closed-loop structure The first edge that belongs to the set of valid edges and forms the closed loop is the... Edge sequence.

[0099] A closed loop structure refers to a chain of edges in a topological graph that starts from a specified starting node and extends forward along interconnected valid edges, with the final destination precisely returning to the starting node, and no other nodes are visited repeatedly except for the starting and ending points. For example, a closed loop structure can include a simple polygonal wireframe with completely connected ends. It serves as a crucial intermediate transition from discrete line segments to continuous surface features, representing a complete and independent closed spatial outline.

[0100] In some embodiments, the processor can obtain a closed loop structure by executing a depth-first or breadth-first closed path tracing algorithm in the topology graph. Specifically, the system selects an unvisited valid edge from the topology graph as the starting point; based on the endpoint node of the current edge, it searches the adjacency list for the next unvisited edge starting from that node and connects them; it continuously expands the path under the constraints of endpoint connection relationships and the spatial direction continuity of adjacent edges; when the latest endpoint node of the current path expansion is detected to coincide with the initial starting node, it is determined that one tracing is completed, and the ordered set of edges is recorded as a closed loop.

[0101] A merged surface set refers to a cleaner, more contiguous set of spatial polygons formed after spatial merging of all initially generated candidate surface features, eliminating local spatial overlaps or adjacent boundary contact issues. For example, a merged surface set may include merging two originally tangent semi-circular forest patches into a single complete elliptical forest patch, or combining a large main forest area and its slightly overlapping subsidiary forest areas into a single composite polygon entity, thereby ensuring that no duplicate calculations occur during spatial extent statistics.

[0102] In some embodiments, the processor can eliminate redundant structures and obtain a fused face set by performing a union operation on spatial geographic information. Specifically, the system uses all generated candidate face features as input layers; it uses the polygon Boolean union algorithm in computational geometry to determine whether there is a spatial intersection or a shared boundary between any two polygons; if so, the outer contours of the two polygons are re-extracted and stitched together, and redundant boundary segments that intersect or touch inside them are erased. After multiple iterations of merging, a final list of face features that are completely mutually exclusive in space is output.

[0103] S6: Based on the forest pixel set and its spatial resolution, calculate the actual forest area, and calculate the envelope area based on the outer envelope surface. Calculate the difference between the envelope area and the actual forest area to obtain the forest encroachment area of ​​the ski resort, thus completing the calculation of the forest encroachment area of ​​the ski resort.

[0104] Spatial resolution refers to the actual size of the smallest physical unit that a remote sensing image sensor or digital raster data can distinguish on the ground. For example, spatial resolution can include values ​​of 10 meters, 30 meters, or higher. A raster data with a spatial resolution of 10 meters means that each square pixel in its image represents a geographical area of ​​10 meters long and 10 meters wide in the real world, or 100 square meters. This metric is the only key multiplier standard for establishing the conversion relationship between the number of pixels in a digital image and the actual physical area of ​​the real world.

[0105] In some embodiments, the processor can obtain the spatial resolution by parsing the header metadata information of a remote sensing image file or a raster dataset. Specifically, when loading an input multispectral remote sensing image or land cover raster file, the system calls the underlying image reading library to parse the metadata tag block containing the geographic reference system; from the tag block, it directly extracts the physical unit values ​​of the pixel span in the east-west (X-axis) and north-south (Y-axis) directions of the image, and usually takes the absolute values ​​of the two as the ground pixel spatial scale parameter uniformly adopted for this layer.

[0106] The actual forest area refers to the total surface area of ​​forest land that currently exists within the target polygonal area bounded by the ski resort and has not been completely destroyed or logged by human construction activities. For example, the actual forest area can include a specific value in square meters or hectares, which reflects the actual remaining scale of the natural forests interspersed within the ski resort. It is one of the two core comparative benchmarks necessary for assessing the extent of forest damage and represents the current state.

[0107] In some embodiments, the processor can calculate the actual forest area by counting the total number of forest pixels and combining this with the physical area of ​​each pixel using scalar multiplication. Specifically, the system first performs a pixel counting operation on the spatially bounded binary raster data of the forest, accumulating the total number of pixels labeled as forest (pixel value 1); then, it obtains the spatial resolution value of the raster data and calculates the actual physical area of ​​a single square pixel, such as the square of the resolution; finally, it multiplies the total number of pixels by the physical area of ​​a single pixel to output the true forest spatial scale value.

[0108] The envelope area refers to the total area of ​​all internal two-dimensional spaces covered by the outer envelope polygon calculated using the adaptive reconstruction algorithm of this invention. For example, the envelope area can include a value representing the theoretically expected total size of forest land assuming the ski slope has never been excavated and the forest remains in a naturally continuous and contiguous state. This area not only includes the currently existing forest land but also fills in and covers the spaces of internal gaps, corridors, or fault zones caused by human interference, representing the original scale of the forest state.

[0109] In some embodiments, the processor can obtain the envelope area by performing geometric polygon area integration under equal-area projection on the finally generated outer envelope vector entity. Specifically, the system extracts the unique outer envelope vector graphic output after all topology repair processing; first, it confirms the spatial reference coordinate system of the graphic; if it is not an equal-area projection, it is forcibly transformed and reprojected to a high-precision equal-area projection coordinate system; then, using Green's formula or polygon vertex coordinate cross product accumulation algorithm based on the geographic information underlying layer, it directly calculates and outputs the theoretical plane area enclosed by the closed contour.

[0110] Forest encroachment by ski resorts refers to the specific area of ​​forest land lost due to the construction of ski resorts, such as leveling mountains to build ski slopes, erecting cable car facilities, and constructing visitor centers, resulting in the felling, destruction, and spatial encroachment of originally continuous natural forests. For example, the forest encroachment area can include a difference accurate to the square meter or hectare. This value is the final output target of the calculation method of this invention and is often used to provide scientific and quantitative technical basis for forestry law enforcement departments to define illegal forest occupation, assess ecological damage, or provide off-site compensation for forest land.

[0111] In some embodiments, the processor can obtain the area of ​​forest encroached upon by the ski resort by performing an arithmetic subtraction operation between the simulated envelope area and the statistically calculated actual forest area. Specifically, the system retrieves the envelope area numerical variable and the actual forest area numerical variable calculated and saved in the previous steps from memory; uses the processor to perform basic floating-point difference calculation, subtracting the actual value representing the current remaining undamaged forest area from the envelope area value representing the theoretical total area of ​​continuous forest; and outputs the calculated area difference result as the final measurement data representing the actual scale of forest encroachment and loss.

[0112] In some embodiments, the calculation of the actual forest area and its envelope area to obtain the forest area encroached upon by the ski resort is expressed as follows: ; ; ; in, This represents the calculated actual forest area. This represents the total number of pixels in a spatially constrained set of forest pixels. This indicates the spatial resolution of the raster data. This represents the calculated envelope area. This refers to the outer envelope surface after topology repair. The function for calculating area. This represents the area of ​​forest encroached upon by the ski resort, obtained by subtracting the two values.

[0113] In some preferred embodiments, to ensure the convergence of the algorithm model and its robustness under different environments, the preset constant parameters used in this system are set within the following empirical range: (1) When constructing the forest probability distribution function When, slope parameter The range of values ​​is set to Bias coefficient The range of values ​​is set to This combination ensures that the probability curve has the best classification steepness. (2) When performing adaptive threshold binarization classification, the adaptive decision threshold is... The range of values ​​is limited to To balance recall and false alarm rate. (3) When calculating the artificial boundary confidence, the geometric parameter weighting coefficients must meet the following constraints: (4) When using fixed empirical values ​​to determine the boundary type, the boundary discrimination threshold The range of values ​​is And is preferred in general application scenarios When the threshold is determined adaptively, the threshold adjustment parameter... The range of values ​​is (5) In three-dimensional scale reconstruction, the curvature adjustment parameter is used to control the magnitude of the influence of curvature on the sampling density. This is usually set as an empirical constant and can be adaptively adjusted according to terrain complexity; while the slope adjustment coefficient is used to calculate the slope influence factor. Its value range is set to To prevent overcompensation. Artificial boundary baseline scale. The width is usually set to the width of the ski slopes at the local ski resort. In the absence of prior statistical data, an empirical standard value range is used by default. The maximum value is used as the alternative benchmark.

Claims

1. A method for measuring the area of forest encroached by ski resorts based on adaptive envelope reconstruction, characterized in that, include: S1: Based on multi-resolution remote sensing image data and land use raster data, pixel feature extraction and adaptive threshold discrimination methods are used to obtain the forest binary raster data of the target area; S2: Using forest binary raster data and ski resort boundary vector data within the target area, spatial overlay processing is performed to extract a spatially constrained set of forest pixels located within the ski resort area; S3: Based on the forest pixel set, extract the spatial boundary and perform discrete point sampling. Use the boundary type discrimination and adaptive adjustment strategy based on local geometric consistency to obtain a unique boundary point set. S4: Based on the unique set of boundary points and digital elevation data, calculate the local terrain slope and three-dimensional spatial distance, and use the three-dimensional circle constraint method combined with slope correction to screen candidate points and construct an effective edge set; S5: Using the set of effective edges, construct a planar topological graph structure and perform closed path tracing to generate multiple candidate surface features. Spatial fusion and area filtering are performed on the candidate surface features, and topological repair processing is executed to obtain the outer envelope surface representing the continuous spatial distribution of the forest. S6: Based on the forest pixel set and its spatial resolution, calculate the actual forest area, and calculate the envelope area based on the outer envelope surface. Calculate the difference between the envelope area and the actual forest area to obtain the forest encroachment area of ​​the ski resort, thus completing the calculation of the forest encroachment area of ​​the ski resort.

2. The method for calculating the forest encroachment area of ​​ski resorts based on adaptive envelope reconstruction according to claim 1, characterized in that, S1 includes: Acquire multi-resolution remote sensing image data and land use raster data; The land use raster data is classified and analyzed, and reflectance information of multiple bands is extracted from the multi-resolution remote sensing image data. Through calculation, the index features related to vegetation, water bodies and shortwave infrared suppression are obtained. Based on exponential features, a forest probability discrimination model is constructed, and an adaptive threshold is used to perform binary classification of each pixel to obtain the forest binary raster data of the target area.

3. The method for calculating the forest encroachment area of ​​ski resorts based on adaptive envelope reconstruction according to claim 2, characterized in that, The expression for the forest probability discriminant model is: ； ； ； ； ； in, This represents the normalized vegetation index at coordinates (x, y). This represents the water index at coordinates (x, y). This represents the shortwave infrared suppression index at coordinates (x, y). This represents the near-infrared reflectance extracted from multi-resolution remote sensing image data. Indicates the reflectivity in the red light band. Indicates the reflectivity in the green light band. Indicates the reflectivity in the shortwave infrared band. Let represent the forest discriminant response function at coordinates (x, y). , and These are the preset weighting coefficients corresponding to vegetation, water bodies, and shortwave infrared suppression features, respectively. Let represent the forest probability distribution function at coordinates (x, y). For the slope parameter, This is the bias coefficient.

4. The method for calculating the forest encroachment area of ​​ski resorts based on adaptive envelope reconstruction according to claim 1, characterized in that, S3 includes: Based on the forest pixel set, the spatial boundary is extracted, and the discrete points of the spatial boundary are extracted as initial sampling points; Calculate the local curvature, rate of curvature change, and neighborhood orientation consistency at each initial sampling point. Dynamically adjust the basic sampling interval based on the local curvature. Construct a boundary discrimination model based on local geometric consistency to classify the spatial boundary into artificial boundaries and natural boundaries. For regions identified as artificial or natural boundaries, the corresponding sampling interval adjustment strategy is executed, and the generated sampling points are deduplicated to obtain a unique set of boundary points.

5. The method for calculating the forest encroachment area of ​​ski resorts based on adaptive envelope reconstruction according to claim 4, characterized in that, The expression for the boundary discrimination model is: ； in, This represents the artificial boundary confidence level of the i-th initial sampling point. This represents the local curvature at the i-th initial sampling point, calculated based on the spatial boundary. This represents the rate of change of curvature at the i-th initial sampling point, calculated based on the curvature difference between adjacent sampling points. This indicates the consistency of the neighborhood orientation at the i-th initial sampling point, calculated based on the angle between the orientations of adjacent sampling points. , and These are the corresponding preset weight coefficients.

6. The method for calculating the forest encroachment area of ​​ski resorts based on adaptive envelope reconstruction according to claim 4, characterized in that, The expression for the basic sampling interval is: ； in, This represents the dynamically adjusted sampling interval for the i-th initial sampling point. Indicates the preset base sampling interval. The preset curvature adjustment parameters, It represents the local curvature at the i-th initial sampling point.

7. The method for calculating the forest encroachment area of ​​ski resorts based on adaptive envelope reconstruction according to claim 1, characterized in that, S4 includes: Obtain digital elevation data, and calculate the local terrain slope at each boundary point based on the unique set of boundary points and the digital elevation data; Based on the local terrain slope and the boundary points, an adaptive scale parameter with a slope correction factor and a three-dimensional spatial distance are constructed. The three-dimensional spatial distance and the adaptive scale parameter are used to filter out a set of three-dimensional candidate connections that meet the conditions. A three-dimensional circular constraint model based on the adaptive scale parameter is established on the local tangent plane, and a terrain consistency joint test is performed on the three-dimensional candidate connection set to construct an effective edge set.

8. The method for calculating the forest encroachment area of ​​ski resorts based on adaptive envelope reconstruction according to claim 7, characterized in that, The expressions for the three-dimensional spatial distance and the adaptive scale parameter are as follows: ； ； ； ； in, Let v represent the three-dimensional spatial distance between the v-th unique boundary point and the u-th unique boundary point. and Let represent the plane coordinates and the elevation value obtained from digital elevation data for the v-th unique boundary point, respectively. and Let represent the plane coordinates and the elevation value obtained from digital elevation data for the u-th unique boundary point, respectively. This represents the slope influence factor at the v-th unique boundary point. This represents the local slope at the v-th unique boundary point calculated based on digital elevation data. The preset slope adjustment coefficient, This represents the adaptive scaling parameter at the v-th uniqueness boundary point. This indicates the preset basic scale of the artificial boundary. Let represent the combined scale of the candidate connections formed by the v-th unique boundary point and the u-th unique boundary point. This represents the function that takes the minimum value.

9. The method for calculating the forest encroachment area of ​​ski resorts based on adaptive envelope reconstruction according to claim 8, characterized in that, The three-dimensional circle constraint model is a candidate circle located on a local tangent plane, with its center... The expression is: ； ； The boundary equation of the candidate circle is: ； in, Let represent the planar distance between the v-th unique boundary point and the u-th unique boundary point. Let represent the coordinates of the midpoint of the line connecting the v-th unique boundary point and the u-th unique boundary point. Indicates the corresponding combination scale, Let represent the unit normal vector determined by the v-th unique boundary point and the u-th unique boundary point. This represents any point on the boundary of the candidate circle.

10. The method for calculating the forest encroachment area of ​​ski resorts based on adaptive envelope reconstruction according to claim 1, characterized in that, S5 includes: A topological graph structure is constructed based on the set of valid edges. Unvisited edges are traversed, and path expansion is performed with the endpoint connection relationship of the edges and the spatial continuity of adjacent edges as constraints to generate a closed loop structure. The internal region of the closed loop structure is polygonized to generate multiple candidate surface features. Spatial fusion operation is performed on each candidate surface feature to eliminate local overlap, resulting in a fused surface set. Calculate the area of ​​each surface element in the fused surface set, retain the main surface element with the largest area and remove fragmented areas, and perform topological restoration processing on the main surface element, including self-intersection elimination, hole filling and fracture repair, to obtain the outer envelope surface representing the continuous spatial distribution of the forest.

Images