Ecosystem monitoring method and system based on data analysis

Abstract

The application discloses an ecosystem monitoring method and system based on data analysis, and belongs to the technical field of remote sensing monitoring. The method comprises the following steps: obtaining multi-temporal remote sensing data of a target mining area and preprocessing the data to obtain a plurality of historical time sequence ecosystem type distribution maps and a current time ecosystem type distribution map; performing time sequence difference operation on the historical distribution maps, extracting boundary pixels of the ecosystem type transition and calculating the transition time to form a transition time grid map; positioning a to-be-checked block whose transition time is later than a preset reference date and cutting the local distribution map; using a change simulator constructed by a generative adversarial network to generate a simulated ecosystem distribution map of the to-be-checked block under the condition that the to-be-checked block is not disturbed by mining; calculating the structural similarity between the local distribution map and the simulated map to identify an abnormal evolution area; and superimposing the abnormal evolution area on an administrative planning base map to generate a warning map.

Description

Technical Field

[0001] This invention relates to the field of remote sensing monitoring technology, specifically to an ecosystem monitoring method and system based on data analysis. Background Technology

[0002] Monitoring of mining area ecosystems relies on comparative analysis of multi-temporal remote sensing images. Current techniques typically involve directly classifying and comparing two or a few periods of remote sensing images to generate change detection maps. This difference or ratio analysis method, oriented towards two periods of images, is prone to misinterpreting seasonal phenological fluctuations and short-term climate anomalies as genuine changes in ecosystem types, generating numerous false alarms in long-term series analysis. The gradual changes in natural ecosystems within mining areas and sudden mining disturbances both manifest as changes in type labels on images. Existing methods struggle to distinguish between natural succession and anomalous transitions caused by human development at the pixel level, thus failing to accurately pinpoint the extent of abnormal evolution caused by mining activities.

[0003] Solving these problems requires overcoming two obstacles. First, how to accurately characterize the time points of qualitative changes in ecosystem types at the pixel level within continuous time series, and eliminate temporary and reversible fluctuations caused by natural factors, thereby stably extracting transition events with irreversible characteristics. Second, how to construct a reliable reference benchmark to infer the ecosystem's distribution state at the current moment without mining activities, in the absence of a real, undisturbed control area, and use this as a comparison benchmark to separate deviations caused by mining development. Existing methods relying on historical averages or point interpolation fail to fully consider the constraints of complex topographic factors on the spatial distribution of ecosystems, resulting in overly smooth or ecologically unsound reference layers, leading to blurred boundaries and inaccurate positioning of anomaly identification results. Summary of the Invention

[0004] This invention provides an ecosystem monitoring method and system based on data analysis, which aims to perform pixel-by-pixel temporal difference and stability testing on multi-temporal remote sensing data, extract reliable transition moments and generate transition moment raster maps, and generate simulated ecosystem distribution maps under undisturbed conditions by combining a change simulator constructed with an adversarial generative network, and accurately identify abnormal evolution areas in mining areas through structural similarity comparison.

[0005] To achieve the above objectives, the present invention provides the following technical solution: The present invention provides an ecosystem monitoring method and system based on data analysis, which is applied to an ecosystem change monitoring platform. This method combines multi-temporal remote sensing data interpretation with intelligent simulation to accurately identify abnormal ecosystem evolution caused by mining development and generate early warning maps with spatial reference value.

[0006] As a technical solution of this invention, multi-temporal remote sensing data of the target mining area is acquired and preprocessed to obtain distribution maps of multiple ecosystem types in the target mining area over historical time series and distribution maps of ecosystem types at the current time. Temporal difference operations are performed on the multiple ecosystem type distribution maps to extract the set of boundary pixels where ecosystem types transition in the historical time series, and the transition time of each boundary pixel is calculated based on the boundary pixel set, forming a transition time raster map. Based on the transition time raster map, the blocks to be verified with transition times later than a preset baseline date are located, and a local distribution map corresponding to the blocks to be verified is cropped from the ecosystem type distribution map. The local distribution map is input into a change simulator constructed by a generative adversarial network to generate a simulated ecosystem distribution map of the blocks to be verified under conditions unaffected by mining development. The pixel-by-pixel structural similarity between the local distribution map and the simulated ecosystem distribution map is calculated, and abnormal evolution areas in the blocks to be verified caused by mining development are identified based on the pixel-by-pixel structural similarity. The abnormal evolution areas are overlaid on the administrative planning base map of the target mining area to generate an early warning map of abnormal ecosystem evolution. By introducing a raster map of transition times for screening, it is possible to effectively distinguish between natural succession and recent human disturbances. At the same time, by using a change simulator to construct a counterfactual comparison, the objectivity and accuracy of identifying abnormal evolution regions can be improved.

[0007] Preferably, the process of acquiring and preprocessing multi-temporal remote sensing data of the target mining area includes: acquiring historical and current remote sensing image sequences of the target mining area; performing radiometric calibration and atmospheric correction on the historical and current remote sensing image sequences respectively to generate a historical surface reflectance stack and a current surface reflectance layer; using a pre-trained deep semantic segmentation network to perform pixel-by-pixel classification on the historical surface reflectance stack and the current surface reflectance layer, outputting distribution maps of multiple ecosystem types in the target mining area over historical time series and distribution maps of ecosystem types at the current time. Through unified spectral preprocessing and deep semantic segmentation, consistency in category definition and spatial accuracy of the multi-temporal classification results is ensured.

[0008] Specifically, when performing pixel-by-pixel classification using a pre-trained deep semantic segmentation network (DSS), a DSS based on a dilated convolutional pyramid pooling module is constructed. The DSS is pre-trained using a publicly available ecosystem type annotation dataset, with the encoder parameters fixed. The single-layer reflectance map of each time point in the historical surface reflectance stack is input into the pre-trained DSS. After encoder downsampling and multi-scale dilated convolution feature extraction, the decoder upsampling restores the original spatial resolution, outputting the ecosystem type distribution map for each time point. The current surface reflectance layer is then input into the same pre-trained DSS, and the same pixel-by-pixel classification operation is performed, outputting the ecosystem type distribution map. The dilated convolutional pyramid pooling module captures multi-scale ground feature characteristics, improving the robustness of classification for complex mining area ecosystems.

[0009] As a preferred embodiment of the present invention, when forming the transition time raster map, multiple ecosystem type distribution maps are arranged in chronological order. A pixel-by-pixel type code comparison is performed on two ecosystem type distribution maps at adjacent times. If the type code at the same pixel location changes, the pixel is marked as a transition candidate pixel, and the next adjacent time is recorded as the initial transition time. For each transition candidate pixel, a stability check of the continuous type code is performed forward and backward along the time axis. Noise pixels whose type codes fluctuate repeatedly in a short period are removed, and the remaining transition candidate pixels that pass the stability check are formed into a boundary pixel set. For each boundary pixel in the boundary pixel set, the time when its type code first changes and remains stable for at least two consecutive times is determined as the transition time of that boundary pixel. The transition times of all boundary pixels are written into a blank raster of the same size as the ecosystem type distribution map, and empty values ​​are filled at the pixel locations where no transition has occurred, generating the transition time raster map. By filtering out spurious changes caused by classification noise through the stability check, the determination of the transition time is made more reliable.

[0010] Furthermore, the process of locating the blocks to be verified and cropping the local distribution map from the ecosystem type distribution map includes: traversing each pixel in the transition time raster map, reading the transition time value stored in that pixel, comparing the transition time value with a preset baseline date, and filtering out all pixels whose transition time values ​​are greater than the preset baseline date to form a pixel group to be verified; performing an octet aggregation operation on the pixel group to be verified, merging octet-connected pixels into a connected region, filtering out isolated noise regions with areas smaller than a preset area threshold, and treating each remaining connected region as a block to be verified, recording the minimum bounding rectangle boundary of each block to be verified; and performing a spatial cropping operation on the ecosystem type distribution map based on the minimum bounding rectangle boundary of each block to be verified, extracting the local distribution map corresponding to each block to be verified. By employing connected region aggregation and area filtering, fragmented noise can be eliminated, focusing on suspicious change blocks with spatial continuity.

[0011] As another preferred technical solution of the present invention, the process of generating a simulated ecosystem distribution map using a change simulator includes: obtaining the geographic coordinate range of the block to be verified, and selecting the closest image to the transition time from the historical remote sensing image sequence as a reference image; constructing an adversarial generative network (PGN) containing a generator and a discriminator, training the GPN using historical image sequences of areas within the target mining area that are not disturbed by mining development and their corresponding ecosystem type distribution maps, fixing the generator parameters after training, and using the trained generator as a change simulator; concatenating the ecosystem type distribution map corresponding to the reference image with the topographic feature encoding vector of the block to be verified and inputting it into the change simulator, which then performs gradual transformations through residual blocks and upsampling layers to output a simulated ecosystem distribution map. By utilizing the GPN to learn the natural succession laws of the mining area and combining topographic feature constraints, the simulation results better reflect the ecological distribution patterns under undisturbed conditions.

[0012] Preferably, when identifying anomalous evolution regions, a sliding window of size N x N pixels is set, and the sliding window is simultaneously slid across both the local distribution map and the simulated ecosystem distribution map. At each window position, the mean, variance, and covariance between the window sub-maps of the local distribution map and the simulated ecosystem distribution map are calculated. Based on the mean, variance, and covariance, combined with preset brightness and contrast constants, a structural similarity index is calculated for each window position, and this index is assigned to the center pixel of the window. After traversing all window positions, a structural similarity map is generated. The structural similarity index of each pixel in the structural similarity map is compared with a preset anomaly threshold. Pixels with a structural similarity index lower than the threshold are marked as anomalous pixels. Connectivity analysis is performed on all anomalous pixels, and isolated anomalous points are removed to form anomalous evolution regions. The structural similarity calculation based on the sliding window fully considers spatial neighborhood information, can more accurately characterize the structural differences between the real and simulated distributions, and suppress random noise interference.

[0013] Specifically, the process of generating an early warning map of abnormal ecosystem evolution includes: loading the administrative planning base map of the target mining area, which includes layers of mining area boundaries, mining right block boundaries, road networks, and residential distribution; extracting vector boundaries of the abnormal evolution areas to generate polygons of the abnormal evolution areas; calculating the geometric center of each polygon and assigning a unique identifier code and risk level label to each polygon; overlaying the polygons of the abnormal evolution areas onto the administrative planning base map with a semi-transparent fill style, and marking the identifier code and risk level label of each polygon at its geometric center; and exporting an early warning map of abnormal ecosystem evolution with a legend. Overlaying the abnormal evolution areas with planning layers such as mining rights and infrastructure can intuitively show the spatial relationship between abnormal evolution and mining development activities, assisting in regulatory decision-making.

[0014] As a refinement of the data acquisition process, when acquiring historical and current remote sensing image sequences of the target mining area, medium-resolution remote sensing satellite images passing over the target mining area are retrieved from publicly available satellite data distribution platforms. Images covering the target mining area with cloud cover below a preset cloud cover threshold are selected from the oldest to the most recent images to form a historical remote sensing image sequence. Images from the same satellite data distribution platform that are the most recent and have cloud cover below a preset cloud cover threshold are acquired as the current remote sensing image. The historical and current remote sensing image sequences are uniformly reprojected onto the same coordinate reference system, and spatial masking is performed according to the administrative boundaries of the target mining area. During the retrieval process, a geographic range query request carrying the set of vector boundary coordinates of the target mining area is sent to a publicly available satellite data distribution platform, specifying the image acquisition time span as a preset start year to the current year. The returned candidate image list is received, from which images with cloud cover greater than or equal to a preset cloud cover threshold, missing bands, or sensor anomalies are removed. The remaining candidate images are sorted from oldest to most recent acquisition time, and one image is extracted at preset time intervals. The extracted images are then arranged chronologically to form a historical remote sensing image sequence. Through rigorous cloud cover and image quality screening, combined with equally spaced sampling, a historical image sequence with uniform temporal distribution and reliable data quality can be obtained.

[0015] During radiometric calibration and atmospheric correction, metadata for each image in the historical remote sensing image sequence is read, radiometric calibration gain coefficients and bias coefficients are extracted, and a linear transformation is performed on the original pixel values ​​of each image to generate a corresponding radiance image. Aerosol optical thickness and water vapor data from the same period as each image are acquired, and atmospheric correction is performed on each radiance image using a lookup table to eliminate atmospheric scattering and absorption effects. Surface reflectance values ​​are output pixel by pixel, and the surface reflectance values ​​from all images in the historical remote sensing image sequence are stacked chronologically to form a historical surface reflectance stack. Radiometric calibration and atmospheric correction are repeated for the current remote sensing image to generate the current surface reflectance layer. This correction process effectively eliminates spectral distortion caused by sensor differences and changes in atmospheric conditions, providing standardized input data for subsequent classification and change analysis.

[0016] This invention also provides a data analysis-based ecosystem monitoring system, including a memory, a processor, and a computer program stored in the memory and running on the processor. When the processor executes the computer program, it implements the steps of the aforementioned data analysis-based ecosystem monitoring method. This system can automatically complete the entire process from data preprocessing, change detection, anomaly simulation to early warning map output, providing efficient technical support for ecological supervision in mining areas.

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

[0018] By arranging multiple ecosystem type distribution maps in chronological order and comparing pixel-by-pixel type codes, and performing temporal difference analysis on the type codes of the same pixel location at different times, the instantaneous locations of ecosystem type changes can be captured. After marking the changed locations as transition candidate pixels, the stability of the continuous type codes of each transition candidate pixel is checked along the time axis. Short-term, repeated jumps caused by natural disturbances or classification noise are eliminated, and only pixels whose type codes remain stable after a change are included in the boundary pixel set. Based on this, the moment when the type code first changes and remains unchanged for at least two consecutive moments is determined as the transition moment of that pixel, and the transition moment is written into the raster map to form a transition moment raster map. This processing flow performs multi-period verification of transition events in the temporal dimension, eliminating spurious change signals caused by seasonal fluctuations and occasional climate events, making the record of transition moments more time-accurate and ecologically significant. Subsequent blocks to be verified based on the transition moment raster map can more accurately reflect the locations of permanent or long-term ecosystem changes.

[0019] When constructing a change simulator using a generative adversarial network (GAN), the most recent historical remote sensing image of the area to be verified, preceding the geographic transition time, is used as a reference image. The ecosystem type distribution map corresponding to the reference image is concatenated with the topographic feature encoding vector of the area to be verified and then input into the generator. The residual blocks and upsampling layers within the generator, while preserving the ecological structure base of the reference period, combine topographic constraints to perform spatial refinement of the ecosystem type distribution within the area, outputting a simulated ecosystem distribution map under conditions unaffected by mining development. This simulated distribution map integrates the spatial constraints of topographic factors and the actual ecological configuration in the reference image, enabling a more precise reconstruction of the ecological pattern controlled by local topography compared to simple interpolation or historical mean synthesis. Using this simulated distribution map as a reference, when performing pixel-by-pixel structural similarity calculations on the local distribution map and the simulated distribution map, anomalous evolution areas show a significant decrease in structural similarity, thus clearly defining the spatial range of significant deviations from the natural evolution trajectory driven by mining development, reducing ambiguous transition zones, and improving the spatial accuracy of anomaly identification. Attached Figure Description

[0020] 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.

[0021] Figure 1This is a flowchart of a data analysis-based ecosystem monitoring method.

[0022] Figure 2 This is a flowchart of the process for generating ecosystem type distribution maps based on multi-temporal remote sensing imagery;

[0023] Figure 3 This is a flowchart of the raster map generation process for ecosystem type transition moments;

[0024] Figure 4 This is a flowchart of the process for extracting the blocks to be verified and clipping the local distribution map based on the raster map of the transition time;

[0025] Figure 5 This is a time-series map showing the change in the area ratio of ecosystem types in the target mining area;

[0026] Figure 6 This is an estimated graph of pixel frequency and kernel density at the moment of ecosystem type transition;

[0027] Figure 7 This is a map showing the area distribution of the smallest bounding rectangle of the block to be verified;

[0028] Figure 8 This is a schematic diagram of the structural similarity index distribution and anomaly threshold;

[0029] Figure 9 It is a scatter plot of the area and risk score of the region with abnormal evolution of the ecosystem. Detailed Implementation

[0030] 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.

[0031] See Figure 1This invention provides a data analysis-based ecosystem monitoring method applied to an ecosystem change monitoring platform. The method acquires and preprocesses multi-temporal remote sensing data of a target mining area to obtain distribution maps of multiple ecosystem types in historical time series and the current ecosystem type distribution map. Temporal difference operations are performed on the multiple ecosystem type distribution maps to extract the set of boundary pixels where ecosystem types transition, and the transition time of each boundary pixel is calculated to form a transition time raster map. Based on the transition time raster map, the blocks to be verified with transition times later than a preset baseline date are located, and corresponding local distribution maps are cropped from the ecosystem type distribution maps. The local distribution maps are input into a change simulator constructed using an adversarial generative network to generate a simulated ecosystem distribution map of the blocks to be verified under conditions unaffected by mining development. The pixel-by-pixel structural similarity between the local distribution map and the simulated ecosystem distribution map is calculated, and based on the pixel-by-pixel structural similarity, abnormal evolution areas in the blocks to be verified caused by mining development are identified. The abnormal evolution areas are overlaid onto the administrative planning base map of the target mining area to generate an early warning map of abnormal ecosystem evolution.

[0032] Example 1:

[0033] In specific implementation, please refer to Figure 2 The process involves acquiring historical and current remote sensing image sequences of the target mining area. Medium-resolution remote sensing satellite images passing over the target mining area are retrieved from publicly available satellite data distribution platforms. Images covering the target mining area with cloud cover below a preset cloud cover threshold are selected from the oldest to the most recent images to form a historical remote sensing image sequence.

[0034] In practice, a geographic range query request is sent to a publicly available satellite data distribution platform. This request includes the vector boundary coordinates of the target mining area and specifies the image acquisition time span as a preset start year to the current year. Upon receiving the query request, the publicly available satellite data distribution platform returns a candidate image list. The ecosystem change monitoring platform receives this list and removes images with cloud cover greater than or equal to a preset cloud cover threshold, as well as images with missing bands or sensor anomalies. The remaining candidate images are sorted by acquisition time from oldest to newest, and one image is extracted at preset time intervals. These extracted images are then arranged chronologically to form a historical remote sensing image sequence. The preset time interval is determined based on monitoring frequency requirements; for example, it can be set to one image per year to ensure the continuity and comparability of the historical sequence. The most recent image from the same satellite data distribution platform with cloud cover below the preset cloud cover threshold is used as the current remote sensing image.

[0035] In practice, historical and current remote sensing image sequences are uniformly reprojected onto the same coordinate reference system, and spatial masking is performed according to the administrative boundaries of the target mining area. The universal transverse Mercator projection can be used as the coordinate reference system to ensure the accuracy of area and distance calculations. Spatial masking uses the administrative boundary vector data of the target mining area to crop each reprojected image, obtaining a subset of images covering only the target mining area.

[0036] In practice, radiometric calibration and atmospheric correction are performed on historical and current remote sensing image sequences, respectively, generating a historical surface reflectance stack and a current surface reflectance layer. Metadata for each image in the historical remote sensing image sequence is read, and radiometric calibration gain and bias coefficients are extracted from the metadata. A linear transformation is performed on the original pixel values ​​of each image to generate a corresponding radiance image. The expression for the linear transformation is:

[0037]

[0038] in, The radiance value of a pixel in a radiance image is expressed in watts per square meter per steradian per micrometer. This represents the radiometric calibration gain coefficient, read from image metadata, used to convert raw pixel values ​​into radiance values. Its value is related to sensor band settings and imaging parameters. This represents the digital quantization value of a pixel in the original image; This represents the bias coefficient, read from image metadata, used to compensate for the sensor's zero-radiative response bias. Its value is determined by sensor laboratory calibration.

[0039] After generating radiance images, aerosol optical thickness and water vapor data contemporaneous with each image are acquired. Atmospheric correction is then performed on each radiance image using a lookup table. The lookup table is pre-generated by simulating the relationship between surface reflectance and entrance pupil radiance under different observation geometries and atmospheric conditions using an atmospheric radiative transfer model. For each pixel in the radiance image, the surface reflectance value after eliminating atmospheric scattering and absorption effects is obtained by quickly indexing or interpolating in the lookup table based on the corresponding observation zenith angle, solar zenith angle, relative azimuth angle, and the aforementioned aerosol optical thickness and water vapor data. Surface reflectance values ​​are output pixel by pixel, and the surface reflectance values ​​output from all images in the historical remote sensing image sequence are stacked in chronological order to form a historical surface reflectance stack. The dimensions of the historical surface reflectance stack are time, rows, columns, and bands, with the time dimension order consistent with the chronological order of the historical remote sensing image sequence.

[0040] Radiometric calibration and atmospheric correction are performed on the current remote sensing image repeatedly, using the same radiometric calibration coefficient acquisition method and atmospheric correction lookup table as historical images, to generate the current surface reflectance layer.

[0041] A pre-trained deep semantic segmentation network is used to perform pixel-by-pixel classification of historical surface reflectance stacks and current surface reflectance layers, outputting distribution maps of multiple ecosystem types in the target mining area over historical time series and at the current moment. In the specific implementation, a deep semantic segmentation network based on a dilated convolutional pyramid pooling module is constructed. The deep semantic segmentation network consists of three parts: an encoder, a dilated convolutional pyramid pooling module, and a decoder. The encoder consists of multiple convolutional modules and downsampling operations. The input is a single-layer surface reflectance map, and the output is a multi-scale feature map, where the highest-level feature map has the smallest spatial resolution and the largest receptive field. The encoder can adopt a deep residual network structure with fully connected layers removed. The dilated convolutional pyramid pooling module applies multiple dilated convolution operations with different dilation rates and a global average pooling operation in parallel to the highest-level feature map output by the encoder to capture multi-scale contextual information. After channel concatenation of the feature maps output by each branch, a 1×1 convolution is applied to compress the number of channels, outputting a feature map that integrates multi-scale information. The decoder receives the feature map output from the dilated convolution pyramid pooling module and gradually restores the spatial resolution through successive transposed convolutions or bilinear interpolation upsampling. After each upsampling, the decoder performs skip connections with the feature map of the corresponding resolution in the encoder to fuse shallow fine details and deep semantic information. The final output is a classification probability map with the same spatial size as the input surface reflectance map. For each pixel, a probability vector of length equal to the number of ecosystem categories is output, and the category code corresponding to the maximum probability is taken as the ecosystem type of that pixel.

[0042] A publicly available ecosystem type annotation dataset was used to pre-train a deep semantic segmentation network. This dataset contains multiple surface reflectance images and corresponding ground truth maps of ecosystem types. During pre-training, the surface reflectance images from the dataset were used as network input, and the corresponding ground truth maps of ecosystem types were used as supervision signals. The cross-entropy loss function was employed to calculate the difference between the predicted results and the ground truth. An adaptive moment estimation optimization algorithm was used to update all network parameters, with an initial learning rate of 0.0001 and a batch size of 16. Pre-training was completed after 100 iterations. After pre-training, all parameters of the encoder in the deep semantic segmentation network were fixed, and these parameters were not updated in subsequent pixel-by-pixel classification of the target mining area.

[0043] The single-layer reflectance map for each time point in the historical surface reflectance stack is input into a pre-trained deep semantic segmentation network. Each time point in the historical surface reflectance stack corresponds to a single layer of surface reflectance map, which has the same number of bands and spatial size as required by the network input. The input surface reflectance map is downsampled by the encoder and subjected to multi-scale dilated convolution feature extraction, and then upsampled by the decoder to restore the original spatial resolution, outputting the ecosystem type distribution map corresponding to that time point. The ecosystem type distribution map is a single-band raster, and the value of each pixel is an ecosystem type code, with the encoding rules consistent with the category definitions used during pre-training. The current surface reflectance layer is input into the same pre-trained deep semantic segmentation network, and the same pixel-by-pixel classification operation is performed to output the ecosystem type distribution map. The pixel values ​​of the ecosystem type distribution map use the same ecosystem type encoding system as the ecosystem type distribution map.

[0044] See Figure 5 The figure shows the time series trend of the area proportion of five ecosystem type codes (type codes 1 to 5) in the target mining area from 2000 to 2025. The horizontal axis represents time, spanning from 2000 to 2025, with a 4-year scale interval. The vertical axis represents the area proportion of each ecosystem type, in percentage, ranging from 0% to 100%. In the legend, different gray blocks represent type code 1 (light gray), type code 2, type code 3, and type codes 4 to 5 (black), respectively, and the width of each block corresponds to the area proportion of each type in the target mining area.

[0045] As can be observed from the graph, the area proportion of Type 1 gradually decreased from approximately 40% in 2000 to approximately 23% in 2025, showing a clear decreasing trend. The decline was mainly concentrated between 2000 and 2015, after which the trend leveled off. The area proportion of Type 2 showed a slow overall upward trend, increasing from approximately 20% in 2000 to approximately 28% in 2025, reflecting the gradual expansion of this type in the mining area ecosystem. The area proportion of Type 3 fluctuated slightly throughout the period, remaining between approximately 15% and 18%, with relatively small changes. The area proportion of Type 4 steadily increased from approximately 15% in 2000 to approximately 25% in 2025, showing a clear and relatively sustained growth trend. The area proportion of Type 5 was relatively stable, remaining in the range of approximately 10%, with the smallest fluctuation.

[0046] This figure reflects the classification results of the multi-temporal ecosystem type distribution map obtained from time-series remote sensing imagery in Example 1, demonstrating the dynamic changes in the spatial distribution of ecosystem types in the target mining area. The changing trends in area proportions lay the foundation for subsequent temporal difference operations and transition moment identification, helping to accurately locate the boundary pixel sets where ecosystem types transition. The overall trend indicates that ecosystem type 1 was affected and its area shrank during the monitoring period, while ecosystems type 2 and type 4 showed an increasing trend, revealing the evolutionary characteristics of the mining area's ecosystem structure.

[0047] Example 2:

[0048] In specific implementation, please refer to Figure 3 Multiple ecosystem type distribution maps are arranged chronologically. Each ecosystem type distribution map corresponds to a time node, and the order of these time nodes is determined by the acquisition time of the historical remote sensing image sequence. A pixel-by-pixel type encoding comparison is performed on two ecosystem type distribution maps at adjacent time points. Adjacent time points refer to two time nodes that are immediately adjacent in chronological order; the ecosystem type distribution map at the previous time point and the ecosystem type distribution map at the next time point form a comparison pair. Specifically, the pixel-by-pixel comparison is performed by traversing every pixel at the same spatial location in both ecosystem type distribution maps. The first type code of the ecosystem type distribution map at that pixel location is read from the previous time point, and the second type code of the ecosystem type distribution map at the same pixel location is read from the ecosystem type distribution map at the next time point. If the first type code and the second type code are inconsistent, the pixel is marked as a transition candidate pixel, and the time point following this adjacent time point is recorded as the initial time point of the transition. The initial time point of the transition identifies the preliminary time position of the change in ecosystem type.

[0049] In practice, for each transition candidate pixel, stability checks are performed on the continuous type encoding along the time axis, both forward and backward. The forward stability check involves tracing back from the initial transition time corresponding to the candidate pixel, reading the type encoding of that pixel's location in the ecosystem type distribution map at several consecutive earlier times, and checking whether the type encoding remains unchanged during the tracing period. The number of consecutive forward tracing times is set to two. If, within the two consecutive forward tracing times, the type encoding of the pixel location remains the same as the first type encoding in the ecosystem type distribution map at the previous time, then forward stability is considered satisfied. The backward stability check involves continuing backward from the initial transition time corresponding to the candidate pixel, reading the type encoding of that pixel's location in the ecosystem type distribution map at several consecutive later times, and checking whether the type encoding remains unchanged during the continued period. The number of consecutive backward tracing times is set to two. If, within the two consecutive backward tracing times, the type encoding of the pixel location remains the same as the second type encoding recorded at the initial transition time, then backward stability is considered satisfied. If a candidate pixel for a transition satisfies both forward and backward stability, it passes the stability test. If a candidate pixel fails to satisfy both forward and backward stability (i.e., its type encoding fluctuates repeatedly within a short period), it is considered a noise pixel and is removed. All remaining candidate pixels that pass the stability test are then grouped into a boundary pixel set.

[0050] In practice, for each boundary pixel in the boundary pixel set, the transition time of that boundary pixel is determined. The type code sequence of the boundary pixels is read from multiple ecosystem type distribution maps arranged chronologically, and the moment when the type code first changes is identified from the type code sequence. The moment of the first change is defined as the moment after which, from the earliest moment, the type code first appears different from the previous moment. A stability verification is performed on this moment of the first change. The condition for stability verification is that, after the moment of the first change, the type code of the boundary pixel remains consistent with the type code after the first change for two consecutive moments. If the stability verification passes, the moment of the first change is determined as the transition time of that boundary pixel. If the stability verification fails, the search continues for the next moment when the type code changes and satisfies the stability verification, until a moment that meets the condition is found as the transition time.

[0051] In practice, the transition times of all boundary pixels are written into a blank raster of the same size as the ecosystem type distribution map. The spatial range, pixel resolution, and number of rows and columns of the blank raster are completely consistent with the ecosystem type distribution map. When the blank raster is created, all pixel positions are filled with null values. For each boundary pixel in the boundary pixel set, its row and column coordinates in the raster, along with the transition time value obtained through the above determination process, are obtained and written into the corresponding pixel position in the blank raster. The transition time value is in date format, including year, month, and day sequence information. For pixel positions where no transition has occurred, the original null values ​​in the blank raster are left unchanged. After writing the transition times of all boundary pixels, the generated raster is the transition time raster map. In the transition time raster map, pixels with values ​​represent positions where ecosystem type transitions have occurred, and pixels with null values ​​represent positions where ecosystem type transitions have not occurred.

[0052] See Figure 6 In the figure, the horizontal axis represents the time of ecosystem type transition, in years, covering the period from 2000 to 2024; the vertical axis represents the pixel frequency at the corresponding transition time. The bar chart shows the distribution of the number of pixels at each transition time point where the ecosystem type transition occurred, and the broken line is a kernel density estimation curve based on the bar data, used to reflect the smooth trend of pixel frequency change at the transition time.

[0053] The bar chart and kernel density estimation curve show that the number of transition pixels was generally low between 2000 and 2009, mainly concentrated around 20 pixels, with slight fluctuations, reflecting relatively slow and stable ecosystem changes during this period. Around 2010, the pixel frequency increased significantly, with the number of transition pixels rapidly increasing to over 50 pixels, and remaining at a high level between 2010 and 2020, peaking at approximately 60 to 65 pixels. This indicates a significant increase in ecosystem type transition activity during this stage, with frequent and densely distributed transitions. After 2021, the number of transition pixels decreased somewhat, and the kernel density curve showed a downward trend, but still remained above 40 pixels, indicating that ecosystem type transition activity weakened but still existed.

[0054] Combining the temporal difference and stability verification steps described in Example 2, the image data reflects the distribution characteristics of the transition boundary pixel set extracted from historical time-series images of the target mining area's ecosystem type over time. Early stages with lower transition pixel frequencies correspond to historical states with less ecosystem type evolution, while stages with significantly increased transition pixel frequencies correspond to key time windows of significant ecosystem transitions, potentially related to ecological disturbances caused by mining development. This distribution characteristic of transition moments provides fundamental temporal data for subsequent identification of abnormal evolution areas, supporting the process in Example 3 of locating the blocks to be verified based on the transition moment raster map. The overall trend accurately reflects the temporal characteristics of the spatiotemporal evolution of ecosystem type transitions, meeting the technical requirements of the transition moment extraction and screening method in this embodiment.

[0055] Example 3:

[0056] In specific implementation, please refer to Figure 4 The system iterates through each pixel in the transition time raster, reading the transition time value stored in that pixel. The transition time raster is a two-dimensional raster data structure, with the same number of rows and columns as the ecosystem type distribution map. Each pixel position corresponds to a transition time value or a null value. The traversal operation is performed row by row and column by column. For the current pixel, the value stored in the current pixel is read. If the value stored in the current pixel is null, the current pixel is skipped, and the traversal continues to the next pixel. If the value stored in the current pixel is a valid transition time value, the transition time value read from the current pixel is compared with a preset baseline date. The preset baseline date is a date threshold pre-set by the ecosystem change monitoring platform to distinguish between natural evolution and ecosystem transitions suspected to be caused by mining development disturbances. The specific value of the preset baseline date is determined based on the time node when large-scale expansion of mining development activities in the target mining area begins, for example, a fixed date before the date of issuance of the mining license for the target mining area. The comparison operation determines whether the transition time value is greater than a preset reference date. If the transition time value is greater than the preset reference date, the current pixel is added to the pixel group to be checked. If the transition time value is less than or equal to the preset reference date, the current pixel is skipped, and the process continues to the next pixel. After the traversal is complete, all pixels in the transition time raster map whose transition time values ​​are greater than the preset reference date form the pixel group to be checked.

[0057] In practice, an 8-connected region aggregation operation is performed on the pixel group to be checked. The connectivity rule for an 8-connected region is defined as follows: for two pixels located at different row and column coordinates, if the two pixels are adjacent in the horizontal, vertical, or diagonal direction (i.e., the absolute value of the difference between their row coordinates is less than or equal to 1 and the absolute value of the difference between their column coordinates is less than or equal to 1), then these two pixels are considered 8-connected. All pixels in the pixel group to be checked are merged into a single connected region using a connected region labeling algorithm. The connected region labeling algorithm is implemented using a two-pass scanning method. The first pass assigns a temporary label to each pixel in the pixel group to be checked and records the equivalence relationships between different temporary labels. The second pass replaces the temporary labels with equivalent relationships with the smallest number, thus forming several connected regions. Each connected region consists of a group of spatially 8-connected pixels that all belong to the pixel group to be checked. Each connected region output by the connected region labeling algorithm corresponds to an initial selected region.

[0058] In practice, isolated noise regions with an area smaller than a preset area threshold are filtered out. For each initially selected region, its area is calculated by counting the total number of pixels it contains. The area of ​​the initially selected region is compared with the preset area threshold, which is set based on the spatial resolution of the remote sensing image and the minimum size of the feature of interest. If the area of ​​the initially selected region is smaller than the preset area threshold, it is identified as an isolated noise region and removed from the set of initially selected regions. If the area of ​​the initially selected region is greater than or equal to the preset area threshold, it is retained. Among all retained initially selected regions, each remaining connected region is treated as a block to be checked.

[0059] In practice, for each block to be verified, the minimum bounding rectangle boundary of the block is recorded. The minimum bounding rectangle boundary is calculated as follows: obtain the row and column coordinates of all pixels contained in the block; extract the minimum and maximum row coordinates from the row coordinates; extract the minimum and maximum column coordinates from the column coordinates; and use the minimum and minimum row and minimum column coordinates as the top-left corner point and the maximum row and maximum column coordinates as the bottom-right corner point to determine a rectangular boundary aligned with the raster row and column directions. This rectangular boundary is the minimum bounding rectangle boundary of the block to be verified. The minimum bounding rectangle boundary of each block to be verified is recorded, and the recorded information includes the minimum row coordinates, maximum row coordinates, minimum column coordinates, and maximum column coordinates.

[0060] In practice, a spatial clipping operation is performed on the ecosystem type distribution map based on the minimum bounding rectangle boundary of each block to be verified. The ecosystem type distribution map is raster data with the same geographic reference and spatial resolution as the raster map at the time of transition. The spatial clipping operation uses the four coordinates of the minimum bounding rectangle boundary as the clipping range, extracting all pixels within the rectangular area bounded by the minimum row coordinate to the maximum row coordinate and the minimum column coordinate to the maximum column coordinate from the ecosystem type distribution map, generating a sub-raster. The extracted sub-raster retains the ecosystem type encoding value of the ecosystem type distribution map. The extracted sub-raster is the local distribution map corresponding to the block to be verified. The above spatial clipping operation is performed on each block to be verified, generating a corresponding local distribution map for each block. The number of rows and columns of the local distribution map corresponding to different blocks to be verified is determined by the range of their respective minimum bounding rectangle boundaries.

[0061] See Figure 7 The chart uses a bar graph to show the frequency distribution of the minimum bounding rectangle area of ​​the blocks to be verified. The horizontal axis represents the size of the minimum bounding rectangle area in hectares, ranging from 0 hectares to approximately 170 hectares. The vertical axis represents the frequency of the blocks to be verified within the corresponding area range. The bar graph shows that there are significantly more blocks to be verified with smaller areas. The frequency is highest in the area between 0 and 10 hectares, and then gradually decreases as the area increases, showing a typical right-skewed distribution trend. The frequency of blocks with an area exceeding 50 hectares decreases significantly, and blocks with an area exceeding 100 hectares are almost rare.

[0062] The preset area threshold is marked by a vertical black dashed line in the graph, which is approximately 5 hectares. In the bar chart, the frequency of blocks below this threshold remains relatively high on the left, while the frequency of blocks exceeding the threshold gradually decreases. According to the description of the block screening in Example 3, this preset area threshold corresponds to the minimum area standard used to filter isolated noise areas. Blocks with an area below this threshold are determined to be isolated noise and are removed, while blocks with an area equal to or greater than the threshold are retained as valid blocks to be verified.

[0063] Overall, the data in the figure reflects the area distribution characteristics of the blocks to be verified after spatial aggregation. The vast majority of blocks are relatively small, with only a few being relatively large. By setting a preset area threshold, a large number of small, isolated noise areas are effectively eliminated, ensuring that subsequent verification of larger-scale ecosystem transition areas that may be affected by mining development is prioritized. This figure specifically demonstrates the effectiveness and rationality of the eight-connected domain aggregation and area threshold filtering steps in Example 3.

[0064] Example 4:

[0065] In practice, the geographic coordinate range of the area to be verified is obtained, and the closest image to the transition time is selected from the historical remote sensing image sequence as the reference base image. The geographic coordinate range of the area to be verified is determined by the actual geographic interval corresponding to the minimum bounding rectangle boundary of the area to be verified in the georeference information of the transition time raster map, including the minimum longitude coordinate, maximum longitude coordinate, minimum latitude coordinate, and maximum latitude coordinate. The specific value of the transition time of the area to be verified is obtained; the specific value of the transition time is the earliest occurrence among the transition times of all boundary pixels within the area to be verified. In the historical remote sensing image sequence, the acquisition time of each image is retrieved in chronological order from oldest to youngest, and all images acquired earlier than the specific value of the transition time are selected. From the selected images, the image whose acquisition time is closest to the specific value of the transition time is selected as the reference base image. The acquisition time of the reference base image represents the surface state of the area to be verified closest to the transition time before the ecosystem type transition.

[0066] In the specific implementation, an adversarial generative network (GAN) is constructed, comprising a generator and a discriminator. The GAN adopts a conditional generative adversarial network architecture. The generator consists of an encoder, a residual block, and an upsampling decoder. The encoder receives the input image and conditional information. The input image is a multi-channel feature tensor formed by concatenating the ecosystem type distribution map corresponding to the reference image with the topographic feature encoding vector of the area to be verified. The encoder consists of three downsampling convolutional modules. Each downsampling convolutional module contains a convolutional layer with a stride of 2, a batch normalization layer, and a leaky linear rectified activation function layer with a negative slope of 0.2. After passing through the three downsampling convolutional modules, the spatial resolution of the input multi-channel feature tensor decreases progressively while the number of channels increases progressively. The residual block section consists of six residual blocks cascaded together. Each residual block contains two convolutional layers, two batch normalization layers, and two rectified linear activation function (RCD) layers with leakage. The input of the residual block is element-wise added to the output of the second batch normalization layer within the residual block to obtain the output of the residual block. The convolutional kernel size of the residual block is 3x3 with a stride of 1 and padding of 1 pixel, without changing the spatial size and number of channels of the feature map. The upsampling decoder section consists of three upsampling modules. Each upsampling module contains a transposed convolutional layer with a stride of 2, a batch normalization layer, and a rectified linear activation function (RCD) layer with leakage. The negative slope of the activation function is set to 0.2. The transposed convolutional layer gradually restores the spatial resolution of the feature map to the same size as the input image. The generator's output layer is a convolutional layer with a kernel size of 7 rows and 7 columns, a stride of 1, and padding of 3 pixels. The number of output channels is equal to the total number of ecosystem types. The output layer uses a normalized exponential activation function to output the probability value of each pixel position belonging to each ecosystem type. The category with the highest probability is taken as the predicted type code for that pixel, forming a simulated ecosystem type distribution map.

[0067] The discriminator employs a Markov discriminator structure, consisting of five concatenated convolutional modules. The first four modules each contain a convolutional layer with a stride of 2, a batch normalization layer, and a leaky linear rectified activation function layer with a negative slope of 0.2 and a kernel size of 4x4. The fifth module contains a convolutional layer with a stride of 1, a kernel size of 4x4, and 1-pixel padding, outputting a single-channel discriminative feature map. The discriminator's input is a joint image formed by concatenating two images along the channel dimension. These two images represent either the input image and the corresponding ground truth ecosystem type distribution map, or the input image and a simulated ecosystem type distribution map generated by the generator. The discriminator outputs a 30x30 discriminative feature map, where the value of each pixel in the feature map represents the probability that the receptive field region of the input image corresponds to the true ecosystem type distribution map rather than the generated simulated ecosystem type distribution map.

[0068] In practice, historical image sequences of undisturbed areas within the target mining area, along with their corresponding ecosystem type distribution maps, are used to train the adversarial generative network (GDN). An area undisturbed by mining development is delineated from the target mining area as a training sample region, defined based on historical mining activity records and field survey results. Surface reflectance maps and corresponding ecosystem type distribution maps for each image in the historical remote sensing image sequence are obtained from the training sample region. The training samples are constructed by selecting any two image pairs from the time-series images of the training sample region. The image from the previous time moment serves as the input reference image, and the ecosystem type distribution map from the image from the next time moment serves as the ground truth distribution map. The ecosystem type distribution map corresponding to the input reference image is concatenated with the topographic feature encoding vector of the training sample region and used as the generator's input. The topographic feature encoding vector contains three topographic parameters for each pixel location in the training sample region: the standard deviation of elevation, the mean slope, and the mean profile curvature. Each parameter is normalized to the 0-1 interval and then concatenated by channel. The topographic parameters are extracted pixel-by-pixel from the digital elevation model (DEM) data. The generator generates a simulated ecosystem type distribution map based on the input. The discriminator is trained to output a discrimination value close to 1 for the joint image of the real ecosystem type distribution map and the input image, and a discrimination value close to 0 for the joint image of the simulated ecosystem type distribution map generated by the generator and the input image. The generator is trained to maximize the probability that the generated simulated ecosystem type distribution map is classified as real by the discriminator, while minimizing the pixel-wise cross-entropy loss between the generated simulated ecosystem type distribution map and the ground truth distribution map. The overall loss function of the adversarial generative network consists of a weighted combination of adversarial loss and pixel-wise class loss, with the weight coefficient of the adversarial loss set to 1.0 and the weight coefficient of the pixel-wise class loss set to 100.0. Training employs an adaptive moment estimation optimization algorithm. The initial learning rate for both the generator and the discriminator is set to 0.0002, the first-order moment decay coefficient is set to 0.5, the second-order moment decay coefficient is set to 0.999, the batch size is set to 1, and the training iterations are 200 epochs. After training, all parameters of the generator are fixed, and the trained generator is used as a change simulator.

[0069] In practice, the ecosystem type distribution map corresponding to the reference benchmark image is concatenated with the topographic feature encoding vector of the area to be verified and then input into the change simulator. The ecosystem type distribution map corresponding to the reference benchmark image is obtained from multiple ecosystem type distribution maps in historical time series. The method of acquisition is to locate the ecosystem type distribution map corresponding to the same time node as the acquisition time of the reference benchmark image. The topographic feature encoding vector of the area to be verified is extracted by calculating the elevation standard deviation, mean slope, and mean profile curvature pixel by pixel from the digital elevation model data of the area where the area to be verified is located. Each parameter is normalized to the interval between 0 and 1 and then concatenated by channel to form a three-channel topographic feature encoding vector consistent with the spatial range of the area to be verified. The type encoding of the ecosystem type distribution map corresponding to the reference benchmark image is converted into a one-hot encoding representation. The number of channels of the one-hot encoding is equal to the total number of ecosystem type categories. The one-hot encoded multi-channel feature map is concatenated with the topographic feature encoding vector in the channel dimension to form the input multi-channel feature tensor. The input multi-channel feature tensor is encoded by a downsampling convolution module in a transformation simulator, and then progressively transformed by residual blocks and upsampling layers to output a simulated ecosystem distribution map of the area to be verified under conditions unaffected by mining development. The simulated ecosystem distribution map is a single-band raster, and the value of each pixel is the ecosystem type code.

[0070] In the specific implementation, a sliding window of size N x N pixels is set, where N is 11. The width and height of the sliding window are 11 pixels, and the row and column coordinates of the window's center pixel determine the area covered by the window on the image. The sliding window is slid synchronously on both the local distribution map and the simulated ecosystem distribution map. Synchronous sliding means that the sliding window moves at the same step size on both maps, and the row and column coordinates of the window's center pixel remain the same. The sliding step size is set to 1 pixel, and the sliding method is to start from the upper left corner of the image and move the sliding window row by row and column by column until the window covers the lower right corner of the image. At each window position, a sub-image defined by the area covered by the sliding window is extracted from the local distribution map, denoted as sub-image A; a sub-image defined by the area covered by the sliding window with the same row and column coordinate range is extracted from the simulated ecosystem distribution map, denoted as sub-image B. At each window position, the mean value of all pixel values ​​in sub-image A is calculated and denoted as... Calculate the mean value of all pixel values ​​within subimage B, denoted as . Calculate the variance of all pixel values ​​within subimage A, denoted as . Calculate the variance of all pixel values ​​within subimage B, denoted as . Calculate the covariance between subgraph A and subgraph B, denoted as . .

[0071] In practice, the structural similarity index for each window location is calculated based on the mean, variance, and covariance, combined with preset brightness and contrast constants. The formula for calculating the structural similarity index is as follows:

[0072]

[0073] in, This represents the structural similarity index calculated at the current sliding window position, with a value ranging from -1.0 to 1.0. The closer the value is to 1.0, the more similar the structures of the two window subgraphs are. This represents the mean of sub-image A, which is a sub-image composed of pixel values ​​within the current sliding window coverage area on the local distribution map. The mean is obtained by summing all pixel values ​​in sub-image A and dividing by the total number of pixels. This represents the mean of subplot B, which is a subplot consisting of pixel values ​​within the current sliding window coverage area on the simulated ecosystem distribution map. The mean is obtained by summing all pixel values ​​in subplot B and dividing by the total number of pixels. The variance of subimage A is represented by the variance calculated by the ratio of each pixel value within subimage A to the variance of the subimage A. The mean of the sum of the squares of the differences is obtained. The variance of subimage B is represented by the variance of each pixel value within subimage B. The mean of the sum of the squares of the differences is obtained. This represents the covariance between subimage A and subimage B, calculated by plotting the covariance of each pixel value within subimage A with... The difference, the pixel value at the corresponding position in subimage B and The mean of the product of the differences is obtained. Represents the brightness constant, and its value is obtained through... Calculations show that Set to 0.01, This represents the dynamic range of the ecosystem type code, which is the difference between the maximum and minimum values ​​of the ecosystem type code. This represents the contrast constant, and its value is determined by... Calculations show that Set to 0.03, With the brightness constant The calculated structural similarity index is assigned to the center pixel of the sliding window. Pixels at the same row and column coordinates as the center pixel in both the local distribution map and the simulated ecosystem distribution map record this structural similarity index value. After the sliding window traverses all window positions, all assigned pixels constitute a structural similarity map.

[0074] In practice, the structural similarity index of each pixel in the structural similarity map is compared with a preset anomaly threshold. The preset anomaly threshold is a structural similarity boundary value set by the ecosystem change monitoring platform to distinguish between normal and abnormal evolution, determined based on the statistical distribution characteristics of structural similarity under simulated natural evolution conditions. If the structural similarity index of a pixel in the structural similarity map is lower than the preset anomaly threshold, that pixel is marked as an anomaly. Connectivity analysis is performed on all marked anomaly pixels, using an eight-connectivity rule. Each connected region output by the connected region marking algorithm is a set of spatially eight-connected pixels, all of which are marked as an anomaly. The total number of pixels contained in each connected region is calculated. If the total number of pixels is less than a preset minimum anomaly region area threshold, the connected region is determined as an isolated anomaly, and all pixel markings within that region are removed from the anomaly pixel markings. All connected regions remaining after removing isolated anomalies constitute an abnormal evolution region.

[0075] See Figure 8 The figure shows the distribution of the structural similarity index and the preset anomaly threshold. The horizontal axis represents the structural similarity index, ranging from 0.1 to 1.0, and the vertical axis represents the number of pixels within the corresponding structural similarity index range. The gray bar chart reflects the distribution of the number of pixels in each structural similarity index range. The bar chart exhibits a clear bimodal characteristic, with the left peak concentrated between 0.1 and 0.5, indicating that some pixels have low structural similarity and may have abnormal changes; the right peak is concentrated between 0.7 and 1.0, indicating that most pixels have high structural similarity and are consistent with normal evolution. The anomaly threshold, marked by a dashed line, is located at 0.55, serving as the boundary between normal and abnormal evolution. As can be seen from the distribution pattern, the number of structural similarity pixels below the anomaly threshold is significantly smaller, mainly concentrated in the low-value range on the left, indicating that these pixels may correspond to areas of abnormal evolution in the ecosystem. This distribution is consistent with the steps of identifying abnormally evolving pixels based on the structural similarity index in Example 4. The threshold setting is reasonable and effectively distinguishes between abnormal evolution caused by mining development interference and areas of natural evolution. The overall data trend reflects the spatial distribution characteristics of anomalous pixels in the structural similarity map obtained by sliding window calculation, providing a quantitative basis for subsequent connected component analysis and risk assessment of anomalous regions.

[0076] Example 5:

[0077] In implementation, an administrative planning base map of the target mining area is loaded. This base map is in vector data format and contains multiple layers. These layers include a mining area boundary layer, a mining right block boundary layer, a road network layer, and a settlement distribution layer. The mining area boundary layer stores the overall boundary line features of the target mining area, with each feature's attributes recording the mining area name and code. The mining right block boundary layer stores the boundary line features of each block within the target mining area that has been granted mining rights, with each feature's attributes recording the mining right license number, the name of the mining right holder, and the validity period. The road network layer stores the road lines within the target mining area, with each feature's attributes recording the road name and road grade. The settlement distribution layer stores the settlement area features within the target mining area, with each feature's attributes recording the settlement name and settlement type. All layers use the same coordinate reference system, consistent with the coordinate reference system used in the raster map at the transition time.

[0078] In practice, vector boundaries are extracted from the abnormal evolution regions. These regions are in raster format, each consisting of a set of connected pixels. Vector boundary extraction employs a boundary tracing algorithm. For each abnormal evolution region's raster data, starting from the top-left pixel, the pixel chain along the region's boundary is traced clockwise. The coordinates of the pixel center points of the traced pixel chain are converted to geographic coordinates, generating a closed polygon boundary. The polygon boundary is stored as a vector polygon feature; one polygon feature is generated for each abnormal evolution region, denoted as the abnormal evolution region polygon.

[0079] In practice, the geometric center of each polygon in the anomalous evolution region is calculated. The geometric center is calculated by extracting the coordinates of all vertices on the boundary of the polygon in the anomalous evolution region, taking the arithmetic mean of the longitude coordinates of all vertices as the longitude coordinate of the geometric center, and taking the arithmetic mean of the latitude coordinates of all vertices as the latitude coordinate of the geometric center.

[0080] In implementation, each polygon representing an abnormal evolution region is assigned a unique identifier code and risk level label. The identifier code consists of three parts: a fixed prefix, a generation date, and a serial number. The prefix is ​​the string "ECO_RISK," the generation date is the system date on which the abnormal evolution region was identified (in the format of eight digits "YYYYMMDD"), and the serial number is a three-digit number, starting from "001" and increasing sequentially in descending order of the abnormal evolution region's area. The risk level label is determined based on a combination of the area of ​​the abnormal evolution region and the minimum structural similarity within that region. The formula for calculating the risk score is defined as follows:

[0081]

[0082] in: The risk score represents the polygon of the abnormal evolution region, ranging from 0.0 to 1.0. The higher the value, the higher the risk level. This represents the area weighting coefficient, with a value set to 0.4, based on the importance of the area factor in risk assessment. The area of ​​the abnormal evolution region is represented in square meters and is obtained by multiplying the total number of pixels contained in the abnormal evolution region by the actual ground area corresponding to a single pixel. This represents the reference area constant, set to 10000.0 square meters, which serves as the benchmark for area normalization and is determined based on the area magnitude of typical mining development disturbance patches in the target mining area. This represents the structural similarity weighting coefficient, with a value set to 0.6, determined based on the importance of the degree of structural similarity reduction in risk assessment. This represents the minimum value among all pixel structural similarity indices within the abnormal evolution region, extracted from pixel values ​​within the corresponding spatial range of the structural similarity map. If... If the value is less than 0.3, the risk level label is set to "low risk"; if If the value is greater than or equal to 0.3 and less than 0.7, the risk level label is set to "medium risk"; if If the value is greater than or equal to 0.7, the risk level label is set to "high risk".

[0083] In practice, the polygons representing abnormally evolving areas are overlaid onto the administrative planning base map with a semi-transparent fill style. The semi-transparent fill style is set as follows: for each polygon representing an abnormally evolving area, the fill color is selected based on its risk level label: yellow for "low risk," orange for "medium risk," and red for "high risk." The transparency parameter of the fill color is set to 40%, i.e., an alpha channel value of 0.4. The boundary line color of the polygons representing abnormally evolving areas is the same as the fill color, the boundary line width is set to 1.5 pixels, and the boundary line transparency is set to completely opaque. The overlay process follows the layer stacking order, with all layers of the administrative planning base map drawn at the bottom layer and the polygon layer representing abnormally evolving areas drawn at the top layer.

[0084] In practice, the identification code and risk level label for each abnormal evolution region are marked at the geometric center of the polygon. The labeling style is set as follows: the identification code uses a 10-point font, black font color, and is displayed 10 pixels above the geometric center point; the risk level label uses an 8-point font, font color consistent with the fill color corresponding to the risk level, and is displayed below the identification code, with a vertical offset distance of the identification code text height plus 4 pixels. A white stroke with a stroke width of 2 pixels is added to the label text to ensure legibility under different background colors.

[0085] In practice, an early warning map of abnormal ecosystem evolution with a legend is exported. The legend includes symbols for mining area boundaries, mining right block boundaries, roads, settlements, and risk level symbols for abnormal evolution areas. The risk level symbols consist of a rectangular color block in the legend and corresponding risk level text. The fill color and transparency of the rectangular color block are consistent with the semi-transparent fill style of the polygons in the abnormal evolution area. The legend text indicates "Low-risk abnormal evolution area," "Medium-risk abnormal evolution area," and "High-risk abnormal evolution area." The early warning map of abnormal ecosystem evolution is exported in Portable Document format or geotagged image file format, with a resolution of 300 dots per inch, and the corresponding geographic coordinate reference information file is also exported.

[0086] See Figure 9 In the figure, the horizontal axis represents the area of ​​the abnormal evolution zone, in hectares, ranging from 0.00 to approximately 2.25 hectares; the vertical axis represents the risk score of the corresponding area, ranging from 0.0 to 1.0. The legend distinguishes three risk levels, with hollow circles representing low risk, triangles representing medium risk, and solid squares representing high risk.

[0087] As can be seen from the data distribution trend in the figure, there is a significant positive correlation between the area of ​​abnormal evolution zones and the risk score. Low-risk areas are mainly concentrated between 0.00 and 0.25 hectares, with risk scores ranging from 0.2 to 0.3, and the overall risk scores are low and dispersed. Medium-risk areas have a wider distribution, extending from approximately 0.00 hectares to 1.00 hectares, with risk scores roughly between 0.3 and 0.7, showing a trend of gradually increasing risk values ​​with increasing area, reflecting the risk level classification influenced by both area and structural similarity. High-risk areas are significantly concentrated in the range of areas larger than approximately 0.75 hectares and with risk scores greater than 0.7, with some areas reaching a maximum risk score of 1.0, indicating that the structural similarity of larger abnormal evolution zones is significantly reduced, and the risk level is extremely high.

[0088] Overall, the areas of abnormal evolution and the ranges of risk scores increase sequentially from low to high risk levels, verifying the rationality of the risk score calculation formula and level classification rules in Example 5. This figure helps to visually demonstrate the risk distribution characteristics of each abnormal evolution area and provides data support for risk level labeling in the ecosystem abnormal evolution early warning map.

[0089] 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 data-driven ecosystem monitoring method, characterized in that, The method, applied to an ecosystem change monitoring platform, includes: Acquire multi-temporal remote sensing data of the target mining area, and preprocess the multi-temporal remote sensing data to obtain distribution maps of multiple ecosystem types of the target mining area in historical time series and distribution maps of ecosystem types at the current time. Perform temporal difference operation on the distribution maps of the multiple ecosystem types to extract the set of boundary pixels where ecosystem types transition in the historical time series, and calculate the transition time of each boundary pixel based on the set of boundary pixels to form a raster map of transition times. Based on the aforementioned transition time raster map, the blocks to be verified that have transition times later than a preset baseline date are located, and a local distribution map corresponding to the blocks to be verified is cropped from the ecosystem type distribution map; The local distribution map is input into a change simulator constructed by an adversarial generative network to generate a simulated ecosystem distribution map of the block to be verified under conditions of no mining development interference; Calculate the pixel-by-pixel structural similarity between the local distribution map and the simulated ecosystem distribution map, and identify the abnormal evolution areas caused by mining development in the block to be verified based on the pixel-by-pixel structural similarity; The abnormal evolution area is overlaid onto the administrative planning base map of the target mining area to generate an early warning map of abnormal ecosystem evolution.

2. The ecosystem monitoring method based on data analysis according to claim 1, characterized in that, Acquire multi-temporal remote sensing data of the target mining area, and preprocess the multi-temporal remote sensing data to obtain distribution maps of multiple ecosystem types in the target mining area over historical time series and distribution maps of ecosystem types at the current time, including: The historical remote sensing image sequence and the current remote sensing image of the target mining area are acquired, and radiometric calibration and atmospheric correction are performed on the historical remote sensing image sequence and the current remote sensing image, respectively, to generate a historical surface reflectance stack and a current surface reflectance layer. A pre-trained deep semantic segmentation network is used to classify the historical surface reflectance stack and the current surface reflectance layer pixel by pixel, and output the distribution maps of multiple ecosystem types of the target mining area in the historical time series and the distribution map of ecosystem types at the current moment.

3. The data analysis-based ecosystem monitoring method according to claim 2, characterized in that, The pre-trained deep semantic segmentation network performs pixel-by-pixel classification on the historical surface reflectance stack and the current surface reflectance layer, outputting multiple ecosystem type distribution maps of the target mining area in historical time series and the ecosystem type distribution map at the current moment, including: A deep semantic segmentation network based on dilated convolution pyramid pooling modules is constructed, and the deep semantic segmentation network is pre-trained using a publicly available ecosystem type annotation dataset, with the parameters of the network encoder fixed. The single-layer reflectance map of each time node in the historical surface reflectance stack is input into the pre-trained deep semantic segmentation network. After encoder downsampling and multi-scale dilated convolution feature extraction, the decoder upsampling restores the original spatial resolution and outputs the ecosystem type distribution map corresponding to each time node. The current surface reflectance layer is input into the same pre-trained deep semantic segmentation network, and the same pixel-by-pixel classification operation is performed to output the ecosystem type distribution map.

4. The data analysis-based ecosystem monitoring method according to claim 3, characterized in that, The step of performing temporal difference operations on the distribution maps of the multiple ecosystem types to extract the set of boundary pixels where ecosystem types transitioned in the historical time series, and calculating the transition time of each boundary pixel based on the set of boundary pixels to form a raster map of transition times includes: The multiple ecosystem type distribution maps are arranged in chronological order. A pixel-by-pixel type code comparison is performed on two ecosystem type distribution maps at adjacent time points. If the type code at the same pixel position changes, the pixel is marked as a transition candidate pixel, and the time after the adjacent time point is recorded as the initial time of the transition. For each transition candidate pixel, the stability of the continuous type encoding is checked along the time axis forward and backward. Noise pixels whose type encoding fluctuates repeatedly in a short period of time are removed, and the remaining transition candidate pixels that pass the stability check are formed into the boundary pixel set. For each boundary pixel in the set of boundary pixels, the moment when its type code first changes and remains stable for at least two consecutive moments is determined as the transition moment of the boundary pixel. The transition moments of all boundary pixels are written into a blank raster of the same size as the ecosystem type distribution map, and the positions of pixels that have not undergone transitions are filled with null values ​​to generate the transition moment raster map.

5. The data analysis-based ecosystem monitoring method according to claim 4, characterized in that, The step of locating the blocks to be verified whose transition times are later than a preset baseline date based on the transition time raster map, and cropping the local distribution map corresponding to the blocks to be verified from the ecosystem type distribution map includes: Traverse each pixel in the transition time raster, read the transition time value stored in the pixel, compare the transition time value with the preset base date, and filter out all pixels whose transition time value is greater than the preset base date to form a pixel group to be checked; An octet aggregation operation is performed on the pixel group to be checked, merging pixels that are octet connected to each other into a connected region, filtering out isolated noise regions with an area smaller than a preset area threshold, and each remaining connected region is a block to be checked, and the minimum bounding rectangle boundary of each block to be checked is recorded. Based on the minimum bounding rectangle boundary of each block to be verified, a spatial clipping operation is performed on the ecosystem type distribution map to extract the local distribution map corresponding to each block to be verified.

6. The data analysis-based ecosystem monitoring method according to claim 5, characterized in that, The step of inputting the local distribution map into a change simulator constructed by an adversarial generative network to generate a simulated ecosystem distribution map of the block to be verified under conditions undisturbed by mining development includes: Obtain the geographic coordinate range of the block to be verified, and select the image closest to the transition time from the historical remote sensing image sequence as the reference baseline image. An adversarial generative network (PGN) containing a generator and a discriminator is constructed. The PGN is trained using historical image sequences of areas within the target mining area that are not disturbed by mining development and their corresponding ecosystem type distribution maps. After training, the parameters of the generator are fixed, and the trained generator is used as the change simulator. The ecosystem type distribution map corresponding to the reference base image is concatenated with the terrain feature encoding vector of the block to be verified and then input into the change simulator. After being transformed step by step by the residual block and upsampling layer in the change simulator, the simulated ecosystem distribution map is output.

7. The data analysis-based ecosystem monitoring method according to claim 6, characterized in that, The calculation of pixel-by-pixel structural similarity between the local distribution map and the simulated ecosystem distribution map, and the identification of anomalous evolution areas caused by mining development in the block to be verified based on the pixel-by-pixel structural similarity, includes: Set a sliding window of size N by N pixels, and slide the sliding window synchronously on the local distribution map and the simulated ecosystem distribution map respectively. At each window position, calculate the mean, variance and covariance between the window subplot of the local distribution map and the window subplot of the simulated ecosystem distribution map. Based on the mean, variance, and covariance, and combined with the preset brightness constant and contrast constant, the structural similarity index of each window position is calculated, and the structural similarity index is assigned to the center pixel of the window. After traversing all window positions, a structural similarity map is generated. The structural similarity index of each pixel in the structural similarity map is compared with a preset anomaly threshold. Pixels with a structural similarity index lower than the anomaly threshold are marked as anomaly pixels. Connectivity analysis is performed on all anomaly pixels, and isolated anomaly points are removed to form the anomaly evolution region.

8. The data analysis-based ecosystem monitoring method according to claim 7, characterized in that, The step of overlaying the abnormal evolution area onto the administrative planning base map of the target mining area to generate an early warning map of abnormal ecosystem evolution includes: Load the administrative planning base map of the target mining area, which includes layers of mining area boundaries, mining right block boundaries, road network, and residential distribution; Vector boundary extraction is performed on the abnormal evolution region to generate abnormal evolution region polygons. The geometric center of each abnormal evolution region polygon is calculated, and a unique identifier code and risk level label are assigned to each abnormal evolution region polygon. The polygons representing the abnormal evolution areas are overlaid onto the administrative planning base map with a semi-transparent fill style. The identification code and risk level label of each abnormal evolution area polygon are marked at its geometric center, and an early warning map of abnormal ecosystem evolution with a legend is exported.

9. The ecosystem monitoring method based on data analysis according to claim 2, characterized in that, The acquisition of historical and current remote sensing image sequences of the target mining area includes: Medium-resolution remote sensing satellite images passing over the target mining area are retrieved from publicly available satellite data distribution platforms. Images covering the target mining area with cloud cover below a preset cloud cover threshold are selected from the oldest to the newest time to form the historical remote sensing image sequence. The image that is the most recent one from the current date and has a cloud cover below a preset cloud cover threshold is obtained from the same satellite data distribution platform and is used as the current remote sensing image; The historical remote sensing image sequence and the current remote sensing image are uniformly reprojected to the same coordinate reference system, and spatial masking is performed according to the administrative boundary of the target mining area; The process of retrieving medium-resolution remote sensing satellite images passing over the target mining area from a publicly available satellite data distribution platform, and selecting multiple time-phase images covering the target mining area with cloud cover below a preset cloud cover threshold according to time from oldest to newest, constitutes the historical remote sensing image sequence, including: Send a geographic range query request to the publicly disclosed satellite data distribution platform. The geographic range query request carries the set of vector boundary coordinates of the target mining area and specifies the image acquisition time span as a preset start year to the current year. Receive the candidate image list returned by the publicly available satellite data distribution platform, remove images from the candidate image list whose cloud cover is greater than or equal to the preset cloud cover threshold, and remove images with missing stripes or sensor malfunctions. The remaining candidate images are sorted from farthest to most recent according to their acquisition time. One image is extracted at each preset time interval, and the extracted images are arranged in chronological order to form the historical remote sensing image sequence.

10. An ecosystem monitoring system based on data analysis, comprising a memory, a processor, and a computer program stored in the memory and running on the processor, characterized in that, When the processor executes the computer program, it implements the steps of the data analysis-based ecosystem monitoring method as described in any one of claims 1 to 9.

Images