Automatic interpretation method for scope and degree of mine ecological destruction based on multi-source data fusion

Abstract

The present application belongs to the technical field of mines, and particularly relates to a method for automatically interpreting the scope and degree of mine ecological damage based on multi-source data fusion. The method comprises the following steps: Step 1: using a UAV platform equipped with an airborne hyperspectral imager and a thermal infrared imager, the position of each pixel in the multi-temporal dual-mode registration data set contains a time sequence composed of multi-dimensional spectral reflectance vectors and surface temperature values in chronological order; Step 2: marking the pixels with an ecological mutation confidence exceeding a preset mutation determination threshold as ecological damage candidate pixels; Step 3: performing fusion of spectral evidence and thermal infrared evidence based on evidence theory on the ecological damage candidate pixels, extracting independent damage areas and classifying the damage degree, and outputting a structured ecological damage evaluation report. The method fully integrates the complementary discrimination advantages of spectral and thermal infrared modalities, and realizes the automatic and quantitative interpretation of the scope and degree of mine ecological damage.

Description

Technical Field

[0001] This invention belongs to the field of mining technology and relates to the use of image recognition technology to judge the ecological damage of mines, specifically to an automatic judgment method for the scope and degree of ecological damage in mines based on multi-source data fusion. Background Technology

[0002] The development of remote sensing technology has provided a wide-ranging, multi-temporal observation method for monitoring mine ecological damage. Currently, commonly used remote sensing data sources include multispectral satellite imagery and high-resolution optical imagery. Multispectral satellite imagery (such as the Landsat and Sentinel-2 series) has a long historical archive and a high revisit frequency, and has been widely used for land use change detection and vegetation cover monitoring in mining areas. However, multispectral imagery has a limited number of bands (usually 4 to 13), and its spectral resolution is insufficient to distinguish land cover types with similar spectral characteristics in mining scenes. For example, the differences in the characteristic absorption positions of different mineral oxides are only on the nanometer scale, and wide-band multispectral imagery cannot capture these subtle differences. While high-resolution optical imagery is rich in spatial detail, it also suffers from insufficient spectral information and mainly relies on visual interpretation or automatic classification based on color and texture features in the visible light band. It cannot acquire information closely related to mine ecological damage, such as surface temperature anomalies and soil moisture changes, which are not reflected in the visible light band.

[0003] Hyperspectral remote sensing technology, with its hundreds of continuous narrow bands covering the entire spectral range from visible light to shortwave infrared, can acquire fine spectral characteristics of ground objects, providing a data foundation for the accurate identification of mining features. In recent years, with the popularization of UAV platforms, airborne hyperspectral imagers have been gradually applied to small-scale, detailed surveys in mining areas. Meanwhile, thermal infrared remote sensing technology can acquire surface temperature information. There are significant differences in thermal radiation characteristics between exposed areas in mining areas, spoil heaps, and tailings ponds and the surrounding natural vegetation; thermal infrared data can serve as an effective supplement to spectral data. However, in existing studies, hyperspectral and thermal infrared data are often used independently, lacking effective pixel-level fusion methods to fully combine the complementary advantages of the two information sources. Even in the few studies that have attempted joint analysis of multi-source remote sensing data, simple feature stitching or weighted voting strategies are often employed, failing to quantitatively model the differences in credibility and potential conflicts between different information sources at the evidentiary level. Summary of the Invention

[0004] The main objective of this invention is to provide an automatic method for judging the extent and degree of mine ecological damage based on multi-source data fusion. It fully integrates the complementary discrimination advantages of spectral and thermal infrared modes, and realizes the automated and quantitative judgment of the extent and degree of mine ecological damage.

[0005] To address the aforementioned technical problems, this invention provides an automatic method for determining the extent and degree of mine ecological damage based on multi-source data fusion, comprising the following steps: Step 1: Using an unmanned aerial vehicle (UAV) platform equipped with an airborne hyperspectral imager and a thermal infrared imager, the target mining area is repeatedly scanned and acquired at multiple different times. The acquired spectral images and thermal radiation images are preprocessed and registered to form a multi-temporal dual-mode registration dataset. Each pixel location in the multi-temporal dual-mode registration dataset contains a time series composed of a multi-dimensional spectral reflectance vector and surface temperature values ​​in chronological order. Step 2: For the spectral reflectance vector of each pixel location in each time phase, spectral features are encoded using a selective state-space model based on the Mamba architecture, outputting a spectral semantic representation vector; the spectral semantic representation vector is concatenated with the thermal feature vector to generate a spectral-thermal joint feature vector, forming a multi-temporal joint feature sequence; the multi-temporal joint feature sequence is input into the Bayesian online change point detection module to identify the moment of ecological state change, and pixels with ecological change confidence exceeding the preset change judgment threshold are marked as candidate pixels for ecological damage; Step 3: Based on evidence theory, perform fusion of spectral evidence and thermal infrared evidence on candidate pixels of ecological damage to determine the confidence distribution of damage type; perform spatial consistency optimization on the confidence distribution of damage type and then perform connected component analysis to extract independent damage regions and classify the degree of damage, and output a structured ecological damage assessment report.

[0006] Furthermore, in step 1, the airborne hyperspectral imager acquires continuous spectral images covering the 400 nm to 2500 nm band range, with no less than 256 spectral channels; the thermal infrared imager acquires thermal radiation images covering the 8 μm to 14 μm band range; the number of different time phases is no less than 3, and the interval between adjacent time phases is no less than 7 days.

[0007] Furthermore, in step 1, strip noise suppression, bad band removal, water vapor absorption band interpolation, and atmospheric radiation correction are sequentially performed on the continuous spectral images of each time phase to obtain a surface reflectance data cube; non-uniformity correction, radiometric calibration, and temperature inversion are performed on the thermal radiation images of each time phase to obtain a surface temperature distribution map.

[0008] Furthermore, in step 1, using the surface reflectance data cube of the first time phase as the reference image, the surface temperature distribution map is resampled by bilinear interpolation and then registered using a registration method based on maximizing mutual information to complete the pixel-level registration; for each subsequent time phase, cross-time phase registration with the reference image of the first time phase is added to ensure that pixels of the same geographical location in all time phases have a strict correspondence.

[0009] Furthermore, in step 2, the processing flow of the selective state-space model includes: setting the hidden state dimension, treating the spectral reflectance vector as a sequence of signals arranged along the band dimension and inputting them sequentially; for the input vector at each band position, generating a state transition modulation vector, an input gate vector, and an output gate vector through independent linear projection layers; the state transition modulation vector is processed by an activation function and then modulated element-wise with the state transition basis matrix to generate an adaptive state transition matrix; the hidden state vector is transformed by the adaptive state transition matrix and added to the input vector after element-wise multiplication by the input gate vector to obtain the updated hidden state vector; the updated hidden state vector is then multiplied element-wise with the output gate vector to obtain the spectral feature output vector; and the spectral feature output vector at the last band position is used as the spectral semantic representation vector.

[0010] Furthermore, in step 2, the execution flow of the Bayesian online change point detection module includes: initializing the running length to 0, setting the observation likelihood model to a multivariate Gaussian distribution, and setting the conjugate prior model to a normal inverse Wieshard distribution; for the spectral-thermal joint feature vector of each time phase arriving in chronological order, calculating the conditional likelihood value for the running segment continuation case and the marginal likelihood value for the change point occurrence case, and normalizing the conditional likelihood value and the marginal likelihood value to obtain the posterior probability of the change point occurrence in the current time phase; after traversing all time phases, selecting the maximum value of the posterior probability of the change point occurrence as the ecological mutation confidence, and recording the time phase position of the mutation and the Euclidean distance between the spectral-thermal joint feature vectors before and after the mutation as the feature deviation magnitude.

[0011] Furthermore, in step 3, fusion is performed based on Dempster-Shafer evidence theory, specifically including: constructing a spectral evidence body based on the spectral semantic representation vector of the post-mutation phase, and outputting basic probability allocation values ​​for five damage types: mining area, spoil heap, tailings dam, exposed surface, and damaged vegetation; calculating the temperature anomaly amplitude based on the deviation of the surface temperature value of the post-mutation phase from the median surface temperature of the non-ecological damage candidate pixels within the local window, and constructing a thermal infrared evidence body based on the temperature anomaly amplitude; and performing orthogonal fusion of the spectral evidence body and the thermal infrared evidence body using Dempster combination rules to obtain the damage type confidence distribution.

[0012] Furthermore, before performing orthogonalization and fusion, the conflict factor between the spectral evidence and the thermal infrared evidence is calculated. If the conflict factor is lower than the preset conflict threshold, the Dempster combination is performed directly. If the conflict factor is not lower than the preset conflict threshold, a discount coefficient that is negatively correlated with the conflict factor is applied to the thermal infrared evidence before the Dempster combination is performed.

[0013] Furthermore, in step 3, spatial consistency optimization is achieved using block-based conditional random field inference. The region consisting of candidate pixels of ecological damage is divided into overlapping sub-blocks of a preset size. A fully connected conditional random field model is constructed for each overlapping sub-block. The unary potential term of the fully connected conditional random field model is composed of the negative logarithm of the damage type confidence distribution, and the binary potential term of the fully connected conditional random field model is composed of the Gaussian kernel function of the spatial distance between pixels and the Euclidean distance of the spectral-thermal joint eigenvector. The mean field approximation method is used to iteratively infer each overlapping sub-block and then the mean value is taken in the overlapping region to output a raster map of the damage type determination result.

[0014] Furthermore, in step 3, the raster map of the damage type determination result is labeled with 8-neighborhood connected components, and the vector boundary polygon and area value of each independent damage area are extracted. Based on the average confidence value of ecological mutation and the average feature deviation of each pixel in each independent damage area, the degree of damage is divided into four levels: mild damage, moderate damage, severe damage and extremely severe damage, using a preset grading threshold or a grading threshold adaptively determined based on regional statistical distribution.

[0015] The automatic judgment method for the scope and degree of mine ecological damage based on multi-source data fusion of the present invention has the following beneficial effects: By using multi-temporal collaborative acquisition and pixel-level spatiotemporal registration of airborne hyperspectral imagers and thermal infrared imagers, spectral reflectance information and land surface temperature information are strictly aligned at the pixel scale and organized according to the time dimension, providing a high-quality multimodal data foundation for subsequent time series analysis. This overcomes the limitation of traditional methods that rely solely on single spectral or single-temporal remote sensing data, resulting in incomplete acquisition of ecological damage information.

[0016] In the feature encoding stage, a selective state-space model based on the Mamba architecture is used to recursively process the spectral reflectance vector along the band dimension. Through an adaptive state transition mechanism, the absorption features of iron oxides, vegetation red edges, and water-bearing minerals that are closely related to mine ecological damage are selectively enhanced in the latent state. Compared with traditional principal component analysis or fixed convolution kernel feature extraction methods, it can more effectively capture diagnostic spectral patterns that change continuously along the wavelength direction and improve the discriminative power of spectral semantic representation.

[0017] In the temporal mutation identification stage, Bayesian online change point detection performs phase-by-phase recursive analysis of multi-temporal joint feature sequences based on a rigorous probabilistic inference framework. It can automatically locate the moment when the ecological state undergoes a sudden change without pre-setting the change type or training a change detection classifier. At the same time, it outputs the ecological mutation confidence with clear probabilistic semantics, avoiding the problem of strong subjectivity and difficulty in quantifying uncertainty in the threshold selection of traditional change detection methods based on threshold comparison.

[0018] In the damage type identification stage, spectral evidence and thermal infrared evidence are orthogonally fused based on Dempster-Shafer evidence theory, fully utilizing the respective discriminative advantages and complementarity of the two sensor modes. Simultaneously, an adaptive discounting mechanism for conflict factors effectively suppresses the negative impact of evidence conflicts on the fusion results. Block-based conditional random field inference eliminates spatial noise generated by pixel-by-pixel independent discrimination while maintaining the true boundaries of the damaged area. This results in a structured ecological damage assessment report with high spatial consistency and quantitative reliability in terms of damage range vector boundaries, damage type, damage area, and damage severity level, directly serving mine ecological restoration planning and regulatory decisions. Attached Figure Description

[0019] Figure 1 This is a schematic diagram illustrating the spectral reflectance curves and diagnostic spectral feature annotations of typical mine features provided in an embodiment of the present invention. Figure 2 This is a schematic diagram of the adaptive state transition response and hidden state evolution of the selective state-space model provided in an embodiment of the present invention; Figure 3 This is a schematic diagram showing the comparison and classification of the degree of destruction of the spatial consistency optimization of the conditional random field provided in the embodiments of the present invention. Detailed Implementation

[0020] The method of the present invention will be further described in detail below with reference to the accompanying drawings and embodiments.

[0021] An automatic method for determining the extent and degree of mine ecological damage based on multi-source data fusion includes the following steps: Step 1: Using an unmanned aerial vehicle (UAV) platform equipped with an airborne hyperspectral imager and a thermal infrared imager, the target mining area is repeatedly scanned and acquired at multiple different times. The acquired spectral images and thermal radiation images are preprocessed and registered to form a multi-temporal dual-mode registration dataset. Each pixel location in the multi-temporal dual-mode registration dataset contains a time series composed of a multi-dimensional spectral reflectance vector and surface temperature values ​​in chronological order. Step 2: For the spectral reflectance vector of each pixel location in each time phase, spectral features are encoded using a selective state-space model based on the Mamba architecture, outputting a spectral semantic representation vector; the spectral semantic representation vector is concatenated with the thermal feature vector to generate a spectral-thermal joint feature vector, forming a multi-temporal joint feature sequence; the multi-temporal joint feature sequence is input into the Bayesian online change point detection module to identify the moment of ecological state change, and pixels with ecological change confidence exceeding the preset change judgment threshold are marked as candidate pixels for ecological damage; Step 3: Based on evidence theory, perform fusion of spectral evidence and thermal infrared evidence on candidate pixels of ecological damage to determine the confidence distribution of damage type; perform spatial consistency optimization on the confidence distribution of damage type and then perform connected component analysis to extract independent damage regions and classify the degree of damage, and output a structured ecological damage assessment report.

[0022] In one specific embodiment of the present invention, a drone platform is used as the flight carrier for remote sensing data acquisition of the target mining area. This drone platform carries two sensors: an airborne hyperspectral imager and a thermal infrared imager. The airborne hyperspectral imager operates using a pushbroom imaging method, acquiring continuous spectral images covering the 400 nm to 2500 nm wavelength range, with no fewer than 256 spectral channels and a spatial resolution better than 0.5 meters. Mine ecological damage involves complex land cover types. Iron oxides exhibit characteristic absorption around 850 nm to 900 nm, vegetation shows a red-edge jump between 680 nm and 750 nm, and hydrous minerals exhibit significant absorption valleys in the 1900 nm to 2200 nm wavelength range. Only an imager covering all of these wavelength ranges and with a sufficiently dense number of channels can ensure that these diagnostic spectral features are not smoothed out due to excessive channel spacing. The thermal infrared imager uses a focal plane array detector, acquiring thermal radiation images covering the 8 μm to 14 μm wavelength range, with a spatial resolution better than 2 meters. This band corresponds to the long-wave infrared region with the highest atmospheric transmittance. Surface radiation energy is minimally affected by atmospheric interference within this range, allowing for a relatively accurate reflection of the actual thermal state of the Earth's surface. The thermal response characteristics of exposed rocks in mining areas, spoil heaps, and tailings ponds under solar radiation differ significantly from the surrounding natural vegetation. Thermal infrared imagers can capture this temperature contrast information, providing auxiliary evidence independent of spectral information for subsequent identification of ecological damage types.

[0023] refer to Figure 1 , Figure 1This diagram illustrates the spectral reflectance curves and diagnostic spectral characteristics of typical landforms in a mine. The horizontal axis represents wavelength in nanometers, ranging from 400 nm to 2500 nm, covering the entire visible, near-infrared, and short-wave infrared bands. The vertical axis represents surface reflectance, ranging from 0 to 0.58, indicating the proportion of solar radiation reflected by various landforms at different wavelengths after atmospheric radiation correction. Five spectral reflectance curves are plotted, corresponding to five typical landform types in the mining area: healthy vegetation, damaged vegetation, bare rock (iron oxides) in the mining area, bare soil in the spoil heap, and tailings ponds (water-bearing minerals). The legend, located in the upper right corner, labels the landform name and line type for each curve. The healthy vegetation curve shows a significant trough near 680 nm caused by chlorophyll absorption, followed by a steep reflectance jump between 680 nm and 750 nm. This jump is the red-edge characteristic of vegetation. After the jump, a high reflectance plateau is maintained in the near-infrared band. The overall morphology of the damaged vegetation curve is similar to that of healthy vegetation, but the magnitude of the red edge rise is significantly weakened, and the reflectance level of the near-infrared plateau is significantly reduced. This difference directly reflects the physiological changes of damaged leaf cell structure and decreased chlorophyll content. The bare rock curve in the mining area exhibits a wide and shallow absorption valley in the 850-900 nm range, which is formed by iron oxide minerals. The crystal field transitions of ions are a key spectral indicator for identifying expansion areas in mining areas. The bare soil curve of the tailings dump shows a slow upward trend from the visible to the shortwave infrared band, lacking significant vegetation features, but exhibits absorption characteristics caused by hydrous minerals and clay minerals in the 1900-1950 nm range and near 2200 nm. The tailings dam curve has a low overall reflectance level, exhibiting the most significant deep absorption valley in the 1900-1950 nm range. This absorption valley is caused by hydrous minerals on the tailings surface. Key and The absorption caused by bond vibrations is a diagnostic spectral feature for identifying the wetted surface of tailings ponds. The figure marks three key diagnostic spectral feature intervals: the 680 nm to 750 nm vegetation red-edge interval, the 850 nm to 900 nm iron oxide absorption valley interval, and the 1900 nm to 1950 nm hydrous mineral absorption band interval. Each interval is marked as a vertical band in the figure and accompanied by textual descriptions. In addition, the figure also marks two water vapor absorption band intervals: 1350 nm to 1460 nm and 1790 nm to 1960 nm. In these two intervals, the strong absorption of solar radiation by atmospheric water vapor molecules makes the surface signal almost impossible for sensors to effectively receive; interpolation processing is required for these bands in actual data processing. Figure 1It can be intuitively seen that the reflectance differences among the five land cover types are the most significant in the aforementioned diagnostic band intervals, which provides a physical basis for the selective state-space model to adaptively enhance the coding response in these band intervals.

[0024] Data collection involves repeated scanning of the target mining area at multiple different time phases, with at least three different time phases and an interval of at least seven days between adjacent time phases. The minimum of three time phases is chosen because Bayesian online change point detection requires accumulating sufficient observation samples over time series to establish stable statistical characteristics for the operational segment. When the number of observation samples is less than three, the estimation of the sample covariance matrix becomes highly unstable, and the posterior probability distribution of the change point detection lacks statistical reliability. The consideration of an interval of at least seven days between adjacent time phases is that the ecological damage process in mining typically involves observable changes on a weekly scale, such as stope expansion, surface weathering of spoil heaps, and vegetation degradation. Repeated collection with too short an interval will generate highly redundant observation data, increasing the data processing burden and potentially causing the change point detection algorithm to misjudge noise fluctuations as ecological mutations. In an optional implementation, when the target mining area is in the rainy season or during periods of intensive mining activity, the time phase interval can be shortened to three to five days, and the number of time phases can be increased to more than twelve to capture rapid changes more precisely. In another alternative implementation, for mining areas with significant seasonal changes, data can be collected at monthly intervals within a complete annual cycle for 12 time phases to distinguish between natural seasonal vegetation changes and man-made mining ecological damage.

[0025] During each data collection flight, the airborne hyperspectral imager and thermal infrared imager operate synchronously. Flight routes are pre-planned based on the target mining area's extent and terrain, with a forward overlap rate set to 60% to 80% and a lateral overlap rate set to 30% to 50%. Flight altitude is determined based on the required spatial resolution and imager field of view, typically ranging from 200 to 500 meters. The inertial navigation system and satellite positioning system on the UAV platform continuously record flight trajectory and attitude parameters during flight, providing fundamental data for subsequent georeferencing.

[0026] For each time-phase acquired continuous spectral image, stripe noise suppression, bad band removal, water vapor absorption band interpolation, and atmospheric radiative correction are sequentially performed to obtain a surface reflectance data cube. Due to slight differences in the response characteristics of individual pixels in the pushbroom hyperspectral imager array, stripe noise perpendicular to the flight direction appears in images along the flight direction. Stripe noise suppression eliminates this systematic bias by statistically equalizing the radiative values ​​of each column of pixels. Bad bands refer to spectral channels with extremely low signal-to-noise ratios or where detector failure renders the data completely unusable. Typical bad bands appear in the detector response edge regions and within strong atmospheric absorption bands. Bad band removal marks the data from these channels as invalid and excludes them from subsequent processing. Water vapor absorption band interpolation targets spectral channels within the two main water vapor absorption ranges of 1350 nm to 1460 nm and 1790 nm to 1960 nm. Due to the strong absorption of water vapor molecules in the atmosphere in these bands, the surface reflection signal is almost completely blocked. Linear interpolation or spline interpolation is performed using the reflectance values ​​of the effective channels on both sides of the absorption range to restore the continuous spectral curve shape.

[0027] Atmospheric radiative correction is a crucial step in converting the radiance values ​​received by the sensor into surface reflectance values. The radiance actually received by the sensor is composed of three parts: surface reflected radiation, atmospheric path radiation, and radiation scattered by neighboring pixels. The core task of atmospheric radiative correction is to separate the reflectance, which is determined solely by the surface reflection characteristics. In one specific implementation, a method based on a radiative transfer model is used to perform radiative transfer calculations based on the atmospheric parameters of the day of data collection (including aerosol optical thickness, water vapor column content, and atmospheric profile) to obtain the atmospheric transmittance. Path radiation and spherical albedo ,in Indicates wavelength. Indicates wavelength Atmospheric transmittance at that location Indicates wavelength The radiation path that is scattered from the atmosphere towards the sensor. Indicates wavelength The hemispherical reflectance of radiation reflected from the Earth's surface by the atmosphere. Surface reflectance. The following relationship is obtained through inversion: ;in Indicates the sensor at wavelength The total radiance value received at the location, Indicates wavelength Solar irradiance at the location, This represents the solar zenith angle. After the above processing, the continuous spectral image of each time phase is converted into a surface reflectance data cube. The spatial dimension of this cube corresponds to the row and column positions of the pixels, and the spectral dimension corresponds to each band channel. The value of each voxel is the surface reflectance of the corresponding pixel at the corresponding band. In an optional implementation, when there is a ground reference target with known reflectance in the flight area, an empirical linear method can be used to directly establish a linear regression relationship between sensor radiance and surface reflectance, replacing the radiative transfer model method to achieve atmospheric correction.

[0028] For each time-phase acquired thermal radiation image, non-uniformity correction, radiometric calibration, and temperature inversion are performed sequentially to obtain a surface temperature distribution map. Non-uniformity correction addresses the inconsistencies in gain and bias among pixels of the focal plane array detector using a two-point correction method. This involves imaging a uniform high-temperature blackbody and a uniform low-temperature blackbody separately, and using the response values ​​at two known temperature points to calculate the gain coefficient and bias of each detector pixel, thus unifying the responses of all pixels to the same reference. Radiometric calibration converts the corrected digital quantization values ​​into physical radiance values. Temperature inversion is based on Planck's radiation law, using the radiance values ​​measured by the sensor to infer the surface temperature. For a thermal infrared imager covering a wide wavelength range of 8 to 14 micrometers, the surface temperature... The inversion relationship can be expressed as: ; in This represents the broadband radiance value measured by the thermal infrared imager. Indicates surface emissivity, This indicates upward radiation from the atmosphere. This indicates downward atmospheric radiation. Indicates atmospheric transmittance in the thermal infrared band. and These are the radiation constants obtained by integrating the Planck function in the 8-14 micrometer band, respectively. Surface emissivity. Emissivity can be estimated using the normalized difference in vegetation index (NDIE). The emissivity of exposed mining areas typically ranges from 0.92 to 0.96, while that of vegetated areas ranges from 0.97 to 0.99. After the above processing, the thermal radiation images for each time phase are converted into surface temperature distribution maps, where the value of each pixel represents the surface temperature in degrees Celsius or Kelvin.

[0029] The purpose of the registration process is to precisely align the two images acquired by the airborne hyperspectral imager and the thermal infrared imager to the pixel level in space. Since the two imagers have different spatial resolutions—better than 0.5 meters for the airborne hyperspectral imager and better than 2 meters for the thermal infrared imager—the surface reflectance data cube from the first time phase is used as the reference image. The surface temperature distribution map from that time phase is then resampled using bilinear interpolation to ensure its spatial resolution matches the reference image. Bilinear interpolation uses a weighted average of the temperature values ​​of the target pixel's four nearest-neighbor pixels in the original image. The weights are inversely proportional to the distance from the target pixel to each nearest neighbor pixel. This method improves resolution while maintaining the continuity of spatial temperature variations, avoiding the step artifacts that might be introduced by nearest-neighbor interpolation.

[0030] After resampling, a registration method based on maximizing mutual information was used to perform pixel-level precise registration between the land surface temperature distribution map and the land surface reflectance data cube. Mutual information measures the strength of the statistical dependency between two images, and its mathematical definition is: ;in and Let represent the random variables of grayscale values ​​of corresponding regions in the two images to be registered. The joint probability distribution is obtained by normalizing the joint gray-level histogram. and These represent their respective marginal probability distributions. When two images are spatially perfectly aligned, the grayscale values ​​of the same ground feature in both images exhibit the strongest statistical correlation, with mutual information reaching its maximum value. Mutual information is chosen as the registration similarity measure instead of the traditional cross-correlation coefficient because spectral and thermal images belong to different imaging modalities. The grayscale representation of the same ground feature in the two images may exhibit a non-linear relationship (e.g., exposed rocks have low reflectivity in the near-infrared band but high thermal radiation intensity). Mutual information can capture this non-linear statistical dependence, while the cross-correlation coefficient is only sensitive to linear relationships. The registration process maximizes the mutual information value by searching for the optimal translation and rotation transformation parameters, completing the pixel-level alignment of the two images within the first temporal phase and generating dual-mode registration data for that phase.

[0031] For each subsequent time phase, the dual-mode registration process described above is first repeated within that time phase, and then cross-time phase registration with the reference image of the first time phase is added. The necessity of cross-time phase registration lies in the fact that the flight trajectories, attitude angles, and positioning accuracy of UAVs in different flight batches cannot be completely consistent. If only the direct georeference of the airborne navigation system is relied upon, the pixel alignment error between different time phases may reach several pixels. This will cause the same pixel location to correspond to different ground locations in different time phases in subsequent time series analysis, resulting in serious deviations in the change point detection results. Cross-time phase registration also adopts the mutual information maximization method. Using the surface reflectance data cube of the first time phase as a fixed reference, affine transformation parameters are searched for on the surface reflectance data cubes of subsequent time phases to maximize the mutual information value. The obtained transformation parameters are synchronously applied to the surface temperature distribution map of the corresponding time phase to ensure a strict correspondence between pixels of the same geographical location in all time phases. In one alternative implementation, when the target mining area has large topographic relief, cross-temporal registration can use orthorectification based on digital elevation model as a preliminary step to eliminate geometric distortion caused by topography before performing mutual information registration, so as to further improve the alignment accuracy of multi-temporal pixels.

[0032] The dual-mode registration data, after undergoing the preprocessing and registration procedures described above, are organized chronologically to form a multi-temporal dual-mode registration dataset. Each pixel location in the multi-temporal dual-mode registration dataset contains one time series, and each temporal element of this time series consists of a 256-dimensional spectral reflectance vector and a 1-dimensional land surface temperature value. Taking three time phases as an example, the time series for a certain pixel location contains three elements: the first element is the 256-dimensional spectral reflectance vector and the 1-dimensional land surface temperature value measured at that location in the first time phase; the second and third elements correspond to the observation data at the same location in the subsequent two time phases, respectively. This multi-temporal dual-mode registration dataset forms the basis of the input data for subsequent spectral feature encoding and temporal abrupt change identification.

[0033] After obtaining the multi-temporal dual-mode registration dataset in step 1, it is necessary to perform spectral feature encoding and temporal abrupt change identification on the time series of each pixel location in the dataset, so as to screen out pixels with significant changes in ecological status as candidates for subsequent damage type discrimination and degree classification.

[0034] For the spectral reflectance vector of each pixel location at each time phase, a selective state-space model based on the Mamba architecture is used for spectral feature encoding. The reason for treating the spectral reflectance vector as a sequential signal rather than directly using dimensionality reduction methods such as fully connected networks or principal component analysis is that the diagnostic spectral characteristics of mine features do not exist in isolation in a single band, but rather exhibit a continuous variation pattern of reflectance values ​​along the wavelength direction within a specific band range. For example, the characteristic absorption of iron oxides is not merely a low point at 850 nm, but rather a complete curve shape in which reflectance gradually decreases from about 780 nm, forms a wide absorption valley in the 850-900 nm range, and then gradually rises after 950 nm; the red edge of vegetation shows a steep leap from the chlorophyll absorption valley near 680 nm to the high reflectance plateau at 750 nm. The identification of these spectral features depends on the contextual relationship between adjacent bands. The selective state-space model accumulates this contextual information gradually along the band dimension through recursive scanning. At the same time, its selective mechanism can dynamically adjust the retention and forgetting ratio of information according to the data content of the current band position, so that the diagnostic spectral range related to mine ecological damage is selectively enhanced in the latent state.

[0035] The specific structure and processing flow of the selective state-space model are as follows. The hidden state dimension is set to 64, which determines the capacity of feature representation within the model. Choosing 64 as the hidden state dimension strikes a balance between feature expressiveness and computational overhead—too low a dimension will fail to adequately represent the complex spectral mixing patterns in the mining scene, while too high a dimension will impose unnecessary computational burden when processing millions of pixels per pixel. In an optional implementation, when the target mining area involves extremely complex land cover types (e.g., the simultaneous presence of multiple mineral oxides, multiple vegetation types, and multiple artificial structures), the hidden state dimension can be increased to 128 to obtain stronger feature expressiveness.

[0036] The 256-dimensional spectral reflectance vector is treated as a sequence of signals arranged along the band dimension and input sequentially into the selective state-space model. Specifically, the first element of the spectral reflectance vector corresponds to the reflectance value of the shortest wavelength channel near 400 nm, and the 256th element corresponds to the reflectance value of the longest wavelength channel near 2500 nm. The processing order proceeds band by band from the short-wavelength end to the long-wavelength end. For each band position in the sequence, the reflectance value at that position is first projected into a 64-dimensional input vector through a 1D convolutional layer. The kernel width of this 1D convolutional layer is set to 4, which means that when generating the 64-dimensional input vector for the current band position, the reflectance value information of the current band and its four adjacent bands is actually utilized simultaneously. This design of the local receptive field ensures that the input vector contains small-scale spectral gradient information from the beginning, which is beneficial for the subsequent hidden state to capture subtle spectral changes such as the position of the absorption valley edge.

[0037] refer to Figure 2 The graph comprises three subplots (top, middle, and bottom) sharing the same horizontal axis range. The top subplot displays a spectral reflectance sequence of iron oxide bare rock input to the selective state-space model. The horizontal axis represents the band position number, ranging from 1 to 256, corresponding to 256 spectral channels from 400 nm to 2500 nm. A reference scale for the corresponding wavelength is appended at the top of the graph. The vertical axis represents the reflectance value. The curve shows a decrease in reflectance due to characteristic absorption of iron oxides in the band position range of approximately 52 to 62 (corresponding to 850 nm to 900 nm), and a dip in reflectance due to absorption by hydrous minerals in the band position range of approximately 184 to 190 (corresponding to 1900 nm to 1950 nm). The middle subplot shows the state transition modulation vector when the selective state-space model processes the above spectral sequence. At each band position The response intensity. The horizontal axis and the upper sub-axis. Figure 1 The vertical axis represents the state transition modulation intensity. The horizontal dashed line in the figure indicates the baseline modulation level, corresponding to the modulation response value in ordinary band intervals without diagnostic spectral features. In three intervals—the iron oxide absorption band (band positions approximately 52-62), the vegetation red edge (band positions approximately 34-42), and the hydrous mineral absorption band (band positions approximately 184-190)—the modulation intensity is significantly higher than the baseline level. This indicates that the selective state-space model actively identifies these diagnostic spectral intervals through an input-driven selective mechanism and generates larger state transition modulation values ​​to adjust the evolution trajectory of the hidden states. In the remaining band intervals, the modulation intensity remains near the baseline level, indicating that the model adopts a conservative default evolution strategy for non-diagnostic bands, not consuming additional representation capacity in these positions. The subfigure below shows four representative dimensions of the 64-dimensional hidden state vector (dimensions 8, 23, 41, and 57, denoted as...). , , , The evolution trajectory along the band position. The horizontal axis represents the band position number, and the vertical axis represents the latent state value. The four trajectories exhibit different evolution patterns: some dimensions show significant slope changes in the iron oxide absorption band, while others show amplitude jumps in the hydrous mineral absorption band, reflecting the division of labor mechanism where different latent state dimensions focus on encoding different spectral features. As the band position progresses from the shortwave end to the longwave end, the latent state values ​​of each dimension gradually accumulate contextual information, finally converging at the 256th band position into a 64-dimensional spectral semantic representation vector containing diagnostic spectral features across the entire band. In each of the three sub-figures, the positions of the iron oxide absorption band, the vegetation red edge, and the hydrous mineral absorption band are marked with vertical band-shaped regions, facilitating the observation of the correspondence between the selective enhancement effect and the diagnostic spectral intervals.

[0038] The selective state-space model maintains a 64-dimensional hidden state vector. subscript This represents the current band position index, ranging from 1 to 256. This hidden state vector is continuously updated during processing from the shortwave end to the longwave end to accumulate spectral context information. Before processing begins, Initialize to a vector of all zeros. For the current band position. 64-dimensional input vector State transition modulation vectors are generated through three independent linear projection layers. Input gating vector and output gating vector All three are 64-dimensional. Each of the three linear projection layers has its own independent weight matrix and bias vector, denoted as follows: , , Both are 64-row, 64-column matrices, and their corresponding bias vectors. , , Taking the state transition modulation vector as an example, its calculation method is as follows: .

[0039] State transition modulation vector After processing by the Softplus activation function, it is modulated element-wise with a preset state transition basis matrix to generate an adaptive state transition matrix for the current band position. The Softplus activation function is defined as follows: ,in Given a scalar input, this function maps any real number to a positive value, and in When the value is large, it approximates the identity mapping; When the value is small, the output approaches zero. Applying Softplus to each component of the state transition modulation vector yields... This ensures that all modulated values ​​are positive. The preset state transition basis matrix... The matrix is ​​a 64-dimensional square matrix, initialized using a diagonalized approximation of the HiPPO matrix. This initialization method ensures that the hidden states tend to retain long-range memories of historical inputs when not modulated by the input. The adaptive state transition matrix is ​​generated by... Each component and Multiplying corresponding elements on the diagonal, the resulting adaptive state transition matrix is ​​denoted as . The core value of this adaptive mechanism lies in the following: when the selective state-space model processes bands closely related to mine ecological damage (such as iron oxide absorption bands, vegetation red-edge bands, or water-bearing mineral absorption bands), the input signal features at these band locations will generate large state transition modulation values ​​through the linear projection layer, causing the adaptive state transition matrix to deviate significantly from the basis matrix. This allows the model to actively adjust the evolution trajectory of the hidden states to enhance the encoding of these diagnostic spectral features. In contrast, at band locations that are not significant for mine ecological discrimination, the state transition modulation values ​​are smaller, and the evolution of the hidden states is close to the default mode determined by the basis matrix. The model will not waste representation capacity on these bands.

[0040] The hidden state is updated as follows: the hidden state vector of the previous band position is updated... Through adaptive state transition matrix After transformation, it is compared with the input gate vector. Input vector after element-wise multiplication Adding them together yields the updated hidden state vector for the current band position, i.e. ,in This represents element-wise multiplication. The input gating vector controls the proportion of the input signal at the current band position that is written into the hidden state. A gating value close to zero means that the hidden state in that dimension mainly retains historical information and is not affected by the current input, while a larger gating value allows the current input signal to significantly update the hidden state. The updated hidden state vector is then... With output gate vector Element-wise multiplication yields the spectral feature output vector for the current spectral band position. The output gating vector determines which dimensions of information in the hidden state are released as output features for the current band position, and which dimensions of information continue to be transmitted only within the hidden state and are not reflected in the output.

[0041] After processing all 256 band positions sequentially according to the above procedure, the 64-dimensional spectral feature output vector of the last band position (i.e., the 256th band position, corresponding to around 2500 nm) is used as the 64-dimensional spectral semantic representation vector for that phase. The output of the last band position is selected as the summary representation of the entire spectral sequence because after a complete recursive scan from the short-wavelength end to the long-wavelength end, the hidden state at the final moment has accumulated contextual information from all 256 bands. Especially under the effect of the selective mechanism, key spectral modes such as iron oxide absorption features, vegetation red edge features, and water absorption features have been selectively enhanced in the hidden state and retained until the final moment. In an optional implementation, the spectral feature output vectors of all 256 band positions can also be globally averaged along the band dimension to obtain a 64-dimensional summary vector as the spectral semantic representation vector. This method has better robustness to situations where the band positions of diagnostic features in the spectral sequence are uncertain.

[0042] After obtaining the 64-dimensional spectral semantic representation vector for a specific time phase, the land surface temperature information for that time phase needs to be incorporated. The 1-dimensional land surface temperature value for that time phase is then expanded into a 16-dimensional thermal feature vector through a linear projection layer. This linear projection layer contains a 16-dimensional weight vector. and a 16-dimensional bias vector Regarding the surface temperature value The projection operation is ,in The resulting 16-dimensional thermal feature vector. Expanding the 1-dimensional scalar to a 16-dimensional vector aims to provide sufficient dimensionality in the feature space to encode the complex mapping between temperature and different damage types—for example, exposed areas in mining areas typically exhibit high-temperature anomalies, damaged vegetation areas show moderate temperature increases, and the moist surface of tailings ponds may exhibit relatively low temperatures. These different temperature patterns require sufficient dimensions for differentiation. The 64-dimensional spectral semantic representation vector and the 16-dimensional thermal feature vector are concatenated along the feature dimensions to generate an 80-dimensional spectral-thermal joint feature vector for this phase. The choice of 16 dimensions, rather than higher, to represent temperature information is based on the fact that the information entropy of temperature as a single physical quantity is far lower than the spectral reflectance of 256 channels; allocating too many feature dimensions to temperature would dilute the dominant role of spectral information in the joint representation.

[0043] The above encoding process is repeated for all temporal phases of the pixel location—that is, for each temporal phase, selective state-space model encoding from the spectral reflectance vector to the spectral semantic representation vector, linear projection from the temperature value to the thermal feature vector, and concatenation of the two are performed separately—resulting in a multi-temporal joint feature sequence composed of multiple 80-dimensional spectral-thermal joint feature vectors arranged in chronological order. Taking three temporal phases as an example, this sequence contains three 80-dimensional vectors, corresponding to the spectral-thermal joint features of the pixel location in the first, second, and third temporal phases, respectively. It is worth noting that each temporal phase shares the same set of selective state-space model parameters and temperature projection layer parameters, which ensures that the encoding results of different temporal phases are within the same feature space, making the subsequent change point detection on the time series comparable.

[0044] Multi-temporal joint feature sequences are input into Bayesian online change point detection to identify moments of abrupt ecological state changes. Bayesian online change point detection is a sequence analysis method strictly based on a probabilistic inference framework. Its core idea is to maintain a posterior probability distribution about the "run length" at each time point of the time series. The run length represents the number of temporal phases processed since the last change point. When the observed data at a certain time point deviates significantly from the statistical characteristics of the existing data in the current run segment, the posterior probability of zero run length rises sharply, indicating that a change point has most likely occurred at that time, that is, the ecological state jumps from one mode to another.

[0045] The specific execution flow of Bayesian online change point detection is as follows. The initial runtime is set to 0. The observation likelihood model is set to an 80-dimensional multivariate Gaussian distribution. This choice is based on the following considerations: there is correlation between the dimensions of the 80-dimensional spectral-thermal joint feature vector (for example, the spectral dimension related to vegetation and the temperature-related dimension tend to change co-evolved during vegetation degradation). The multivariate Gaussian distribution can characterize this inter-dimensional correlation structure through the covariance matrix. The conjugate prior model is set to a normal inverse Wissaud distribution. This distribution is the natural conjugate prior of the multivariate Gaussian distribution when both the mean vector and covariance matrix are unknown. Its advantage is that whenever a new observation arrives, the posterior distribution still maintains the form of a normal inverse Wissaud distribution. Bayesian inference can be completed simply by updating sufficient statistics, without the need for numerical integration or sampling. The normal inverse Wissaud distribution is controlled by four hyperparameters: the prior mean vector... (80-dimensional), mean prior strength scalar Scalar of the degrees of freedom of the inverse Wissaud distribution and the prior scale matrix (80 by 80 dimensions). In one specific implementation, Let it be the global mean of the joint feature vector of all pixels in the training sample. Set to 1, Set it to 81 (i.e., add 1 to the feature dimension, which is the minimum legal value of the degrees of freedom of the normal inverse Wissaud distribution to ensure the existence of the expected covariance matrix). Let the training sample covariance matrix be multiplied by The result.

[0046] Set the prior probability of the change point occurring. It is the reciprocal of the time phase interval. For example, when 12 time phases are collected, The intuitive meaning is that, in the absence of other information, it is assumed prior to be that the probability of a change point occurring between any adjacent time phases is 1. In one alternative implementation, if the mining production plan for the target mining area is known, the prior probability can be adjusted according to the expected frequency of mining activities in the plan.

[0047] For the 80-dimensional spectral-thermal joint feature vector of each phase arriving in chronological order in the multi-temporal joint feature sequence (in (Representing the phase number, counting from 1), perform the following recursive update. At each time step... running length The value range is from 0 to For the case where the running segment continues, i.e. It is necessary to calculate the current observations The conditional likelihood value given the accumulated data within the current running segment. Specifically, suppose the current running segment starts from time [time]. Start to Time Total of Based on these observations, the sufficient statistic of the normal inverse Wissaud distribution has been updated to the posterior mean. posterior mean intensity Posterior degrees of freedom and the posterior scale matrix Then the current observation The conditional likelihood (also known as the posterior predictive probability density) follows a multivariate order. distributed: ;in The mean is The scale matrix is Degrees of freedom are diversity The probability density function of the distribution. denoted as the dimension of the spectral-thermal joint eigenvector. When the accumulated observations within a running segment are few, the tail of the posterior prediction distribution is thick (i.e., the uncertainty is large), and it gradually approaches a multivariate Gaussian distribution as the observations accumulate.

[0048] For the situation where the change point occurs, i.e. The current observation is considered the first observation of a completely new operating segment, and its marginal likelihood value is calculated based on the prior model, also following a multivariate model. Distribution but with prior hyperparameter , , , Substitute into the formula above.

[0049] Combine the conditional likelihood of the run segment continuation scenario with the prior probability of run segment continuation. Multiplying these yields the joint probability of the continuation hypothesis; the marginal likelihood of the change point occurring is then combined with the prior probability of the change point occurring. Multiplying them yields the joint probability of the change point hypothesis. Normalizing the two joint probabilities—that is, dividing each by the sum of the two—gives the current time phase. Post-transformation probability If the current time phase is determined to be a continuation of the running segment (i.e., the change point has not occurred), then the current observations are used. Update sufficient statistics for the running segment: , , , superscript This represents the vector transpose. If the current phase is determined to be a change point, the run length is reset to 0 and sufficient statistics for the new run segment are initialized with prior hyperparameters.

[0050] After traversing all time phases, the maximum posterior probability of each phase transition point is selected as the confidence level of ecological abrupt change at that pixel location. This selection reflects the following logic: a pixel is sufficient to be included in the candidate range for further analysis if it experiences a high-confidence ecological state abrupt change at some point during the entire observation period. When the confidence level of ecological abrupt change exceeds a preset abrupt change judgment threshold, the pixel is marked as a candidate pixel for ecological damage. In one specific implementation, the abrupt change judgment threshold is set to 0.7. This value is set based on the following: in experimental verification of multiple mining areas, when the threshold is below 0.6, a large number of spurious change points caused by sensor noise or atmospheric correction residual errors are incorrectly included; when the threshold is above 0.85, some gradual ecological damage (such as slow vegetation degradation processes) is missed. The optimal balance between recall and precision is achieved around 0.7. In an optional implementation, the threshold can be adjusted within the range of 0.5 to 0.9 according to the specific characteristics of the target mining area, or the optimal operating point determined by the receiver operating characteristic curve on the verification dataset can be used.

[0051] Simultaneously, the temporal position of the mutation and the Euclidean distance between the spectral and thermal joint eigenvectors before and after the mutation are recorded as the characteristic deviation magnitude. The characteristic deviation magnitude is calculated by taking the 80-dimensional spectral and thermal joint eigenvector of the most recent temporal phase before the mutation point. 80-dimensional spectrothermal joint eigenvectors of the phase at which the change point occurs ( (where the time sequence number of the change point is used) to calculate the Euclidean distance between the two. ,in and Represent the first and second vectors respectively. Each component. The greater the deviation of the feature, the more drastic the change in the combined spectral and thermal characteristics of the pixel before and after the mutation, which will serve as an important basis for classifying the degree of damage in subsequent steps.

[0052] The above processing is performed on each pixel location in the multi-temporal dual-mode registration dataset one by one. Since the processing flow for each pixel is completely independent, it can be accelerated through parallel computing in practical engineering implementation. In an optional implementation, all pixels of the entire image can be organized into batches, and the parallel computing capabilities of the graphics processor can be used to process the spectral feature encoding and change point detection processes of thousands of pixels simultaneously.

[0053] The following describes the specific implementation process of the structured output, which defines the scope and degree of disruption under the constraints of multi-source data fusion and spatial consistency in step 3.

[0054] For all candidate pixels of ecological damage, the Dempster-Shafer evidence theory is used to fuse spectral and thermal infrared evidence to determine the damage type. Compared to traditional Bayesian classification methods, the Dempster-Shafer evidence theory has the unique advantage of allowing beliefs to be assigned to any subset of focal elements, not just a single category, thus enabling the explicit expression of "uncertainty" and "don't know" states. In the identification of mine ecological damage types, spectral and thermal infrared information each have types of damage they are good at identifying and types that are difficult to distinguish. The framework of the evidence theory allows each information source to assign probabilities to the union of multiple types when its own confidence is insufficient. After the two sources of evidence are fused, the intersection operation is used to focus on the specific type.

[0055] A spectral evidence body is constructed based on the spectral semantic representation vector of the post-mutation phase. Specifically, this 64-dimensional spectral semantic representation vector is input into a pre-trained 5-class fully connected layer, which contains a weight matrix. (5 rows and 64 columns) and bias vector (5-dimensional), calculation ,in It is a 64-dimensional spectral semantic representation vector. This is a 5-dimensional output vector. Then... Apply Softmax normalization, ,in The first in the spectral evidence body The basic probability assignment values ​​for each type of destruction. Indicates the first The focal element corresponding to each type of damage. The 5-dimensional output vector is the first Five components are used. The five damage types correspond to mining areas, spoil heaps, tailings ponds, exposed surface, and damaged vegetation, respectively. The parameters of this fully connected layer are trained on a pre-labeled dataset of mine ecological damage samples with the objective of minimizing cross-entropy loss. After training, the parameters remain fixed during the inference phase. The spectral semantic representation vector, after being encoded using a selective state-space model, has condensed key spectral information across 256 bands. The 5-class fully connected layer maps this information to the support for each damage type.

[0056] The temperature anomaly amplitude is calculated based on the deviation of the surface temperature value at the time of the abrupt change from the median surface temperature of non-ecologically damaged candidate pixels within the local window containing that pixel. The median, rather than the mean, of non-candidate pixels is chosen as the background temperature reference because the median is more robust to extreme values—there may be a few extremely hot heat sources (such as spontaneous combustion points) in local areas of the mining area. If the mean is used, the background temperature estimate will be inflated by these extreme values, leading to a systematic underestimation of the temperature anomaly amplitude. In one specific implementation, the size of the local window is set to a 51 x 51 pixel area centered on the current pixel, corresponding to a ground area of ​​approximately 25.5 meters by 25.5 meters at 0.5-meter resolution. The temperature anomaly amplitude is defined as... ,in This represents the surface temperature value of the current candidate pixel after the abrupt change. This represents the median of the surface temperature values ​​of non-candidate pixels within the local window.

[0057] Based on the magnitude of temperature anomalies Thermal infrared evidence is constructed using pre-defined segmented mapping rules, outputting basic probability allocation values ​​for five damage types. Different damage types exhibit regular differences in temperature anomalies: mining areas, due to large areas of exposed rock that is typically dark in color (stronger absorption of solar radiation), often show the largest temperature anomalies, with typical values ​​exceeding 15 degrees Celsius; spoil heaps and exposed surface areas lack the cooling effect of vegetation transpiration, resulting in moderate temperature anomalies, typically ranging from 8 to 15 degrees Celsius; damaged vegetation areas, despite reduced vegetation cover, still have some remaining vegetation, leading to relatively smaller temperature increases, typically ranging from 3 to 8 degrees Celsius; tailings pond surfaces, due to high water content, experience relatively low temperature anomalies and complex variation patterns due to evaporative cooling. Based on these thermal characteristics of the ground features, the design of the segmented mapping rules follows the following logic: when... When the temperature exceeds 15 degrees Celsius, the basic probability allocation value for the mining area type is set to the highest (0.6 in one specific implementation), and the remaining probability is equally distributed among the other four types; when When the temperature is between 8 and 15 degrees Celsius, the basic probability allocation value for exposed surface type is set to the highest (e.g., 0.5), followed by spoil heap type (e.g., 0.25), and the remaining probabilities are equally distributed among the other types; when When the temperature is between 3 and 8 degrees Celsius, the basic probability allocation value for the damaged vegetation type is set to the highest (e.g., 0.5), and the remaining probabilities are equally distributed among the other types. In an alternative implementation, the temperature threshold for segmented mapping can be adjusted according to the seasonal temperature background of the climate zone where the target mining area is located.

[0058] Before fusing spectral evidence and thermal infrared evidence, the conflict factor between the two pieces of evidence must first be calculated. Defined as: ;in For spectral evidence, the focal element The basic probability allocation value, For thermal infrared evidence, the focal element The basic probability distribution value is summed by traversing all focal elements, resulting in an empty set. The conflict factor reflects the degree of contradiction between two pieces of evidence; when the two pieces of evidence are highly consistent in their judgment of the type of destruction... Approaching zero, when the judgments of the two sets of evidence are severely contradictory. Approaching 1. When the conflict factor is lower than the preset conflict threshold (set to 0.3 in one specific implementation), the Dempster combination rule is directly applied to perform orthogonalization and fusion. The Dempster combination rule applies to each focal element in the fusion result. The basic probability allocation values ​​are calculated as follows: The numerator is the intersection of all focal elements in the two evidence bodies, which is exactly equal to... The sum of the products of the basic probability distribution values, denominator The normalization factor is used to ensure that the sum of the basic probability assignments for the fusion of all focal elements is 1. When the conflict factor is not lower than the preset conflict threshold, it indicates a significant discrepancy between spectral evidence and thermal infrared evidence. In this case, if the Dempster combination rule is applied directly, the normalization factor... A very small value can lead to an overemphasis on the probability of a few consistent focal elements in the fusion results, resulting in counterintuitive conclusions. To address this issue, a discount factor is applied to the thermal infrared evidence. Then, the Dempster combination is executed, and the discount factor is negatively correlated with the conflict factor. In one specific implementation, The thermal infrared evidence, after applying a discount, becomes... (For all non-universal focal elements) ,in This represents the focal element of the complete set containing all five destruction types. These are the probability values ​​assigned to the entire set in the original thermal infrared evidence. The effect of the discounting operation is to reduce the influence of thermal infrared evidence when there is significant conflict between the two pieces of evidence, making the fusion result rely more on spectral evidence. The reason for applying discounting to thermal infrared evidence rather than spectral evidence is that spectral information contains rich features across 256 bands and has undergone deep encoding by a selective state-space model, and its ability to discriminate damage types is generally superior to thermal infrared information that relies solely on the amplitude of a single temperature anomaly. After fusion, the damage type confidence distribution for each ecological damage candidate pixel is obtained, i.e., the fused basic probability assignment values ​​for each of the five damage types.

[0059] After obtaining the confidence distribution of damage types for all candidate ecological damage pixels, spatial consistency optimization is performed. The necessity of spatial consistency optimization lies in the fact that pixel-by-pixel change point detection and evidence fusion inevitably generate spatial "salt-and-pepper noise"—adjacent pixels, even if belonging to the same damage area on the ground, may be assigned to different damage types due to local spectral noise or temperature measurement fluctuations. Spatial consistency optimization is achieved using block-based conditional random field inference. The spatial region formed by candidate ecological damage pixels is divided into overlapping sub-blocks of a preset size. In one specific implementation, the sub-block size is set to 64 x 64 pixels, and the overlap width between adjacent sub-blocks is set to 16 pixels. A fully connected conditional random field model is constructed for each sub-block. The energy function of this model includes two parts: a univariate potential energy term and a binary potential energy term. The univariate potential energy term is composed of the negative logarithm of the damage type confidence distribution for each pixel. Assigned to type The univariate potential energy is ,in The first pixel after fusion The basic probability assignment values ​​for each type of destruction. The role of the univariate potential term is to maintain the original preference of each pixel for each type of destruction in the fusion result. The binary potential term is composed of a Gaussian kernel function of the spatial distance between pixels and the Euclidean distance of the spectral-thermal joint eigenvector. For each pixel... With pixels The binary potential energy when assigned to different types is: ; in and Each pixel and pixels Spatial coordinates (2D) and Each pixel and pixels The 80-dimensional spectral-thermal joint eigenvector, This refers to the spatial distance bandwidth parameter. For the feature distance bandwidth parameter, The spatial distance bandwidth parameter of the smoothing kernel. and These are the weighting coefficients of the two kernel functions. The first Gaussian kernel considers both spatial distance and feature similarity, causing spatially adjacent and feature-similar pixels to be assigned to the same type, while the constraint between spatially adjacent but feature-different pixels (such as pixels at the boundary between two disruptive types) is weaker, thus protecting the true type boundaries from over-smoothing. The second Gaussian kernel only considers spatial distance, providing a basic spatial smoothing effect. In one specific implementation, Set to 3 pixels, Set to 0.5 (based on the Euclidean distance scale after normalization of the joint feature vectors). Set to 1 pixel, Set to 5. Set it to 3.

[0060] refer to Figure 3The image contains three sub-images: left, center, and right. The horizontal and vertical axes of these sub-images represent the row and column numbers of pixels, indicating the spatial location of an 80x80 pixel local area within the mining area image. The left sub-image shows a raster map of the damage type determination results obtained from Dempster-Shafer evidence fusion performed independently pixel-by-pixel before performing Conditional Random Field spatial consistency optimization. The image includes six categories: background, mining area, spoil heap, tailings pond, exposed surface, and damaged vegetation. Each category is represented by a different grayscale level or pattern, and the legend labels each category. It can be seen that although the main outlines of the five damage types are generally discernible, there are numerous scattered pixels with inconsistent type labels within and at the boundaries of each damaged area. These pixels differ from the type determination results of their surrounding neighborhoods, creating significant spatial noise. This spatial noise originates from the local spectral noise of individual pixels during independent pixel-by-pixel processing, temperature measurement fluctuations, and random uncertainties in evidence fusion, failing to utilize the prior constraint that adjacent pixels within the same damaged area should have spatial consistency. The middle sub-image shows the raster map of the damage type determination results after block-based conditional random field spatial consistency optimization. Compared with the left sub-image, the scattered noise within each damage region is effectively eliminated, the region boundaries are smoother and more coherent, and pixels of the same damage type form continuous and complete patches in space. Meanwhile, the boundaries at the intersection of different damage types still maintain clear demarcation and are not excessively blurred due to spatial smoothing. This is attributed to the constraint effect of the Euclidean distance of the spectral-thermal joint eigenvector in the binary potential term of the conditional random field—the smoothing constraint imposed by the binary potential term on pixels on both sides of the true type boundary is weaker due to the large difference in joint features, thus protecting the sharpness of the boundary. The right sub-image shows the results of classifying each independent damage region into four levels of damage severity based on the mean confidence score of ecological mutation and the mean amplitude of feature deviation. The figure divides all damage regions into four levels: mild damage, moderate damage, severe damage, and extremely severe damage, represented by different gray levels. Spatially, the core area of ​​the mining area was identified as severely damaged, while the surrounding transition zone was severely damaged; the center of the spoil heap area was severely damaged, while the edges were moderately damaged; the tailings dam area was generally moderately damaged; and the exposed surface and damaged vegetation areas were mainly lightly to moderately damaged. This gradient distribution of damage from the core to the edge conforms to the general spatial pattern of mine ecological damage, that is, the intensity of damage decreases from the center of the disturbance source outward.

[0061] A mean-field approximation method is used to iteratively infer the fully connected conditional random field model within each sub-block. The mean-field approximation decomposes the joint probability distribution into the product of the edge distributions of each pixel, and iteratively updates the type assignment probability of each pixel by minimizing the KL divergence between this product distribution and the true joint distribution. In each iteration, for each pixel... Update its assigned type The probability is: ;in For the pixels in the previous iteration Assigned to type The probability of each pixel type is determined. The iteration count is set to 10. Experimental results show that the type assignment probability of each pixel converges after 10 iterations. After iterative inference for each sub-block, the average type assignment probability of each pixel in the overlapping area between sub-blocks is taken to eliminate the sub-block boundary effect. Finally, the type with the highest probability for each pixel is taken as the judgment result, and a spatially smoothed raster map of the damage type judgment result is output. In an optional implementation, when the spatial distribution of ecological damage candidate pixels is sparse and the area is small, a single fully connected conditional random field model can be directly constructed for inference of all candidate pixels without using a block-based strategy.

[0062] The raster map of the damage type determination results is labeled with 8-neighbor connected components. This labeling aggregates all spatially adjacent cells (horizontally, vertically, or diagonally) that are determined to have the same damage type into the same connected component, with each component corresponding to an independent damage region. In the implementation, the raster map is scanned row by row from the top left to the bottom right. When an unlabeled candidate cell is encountered, a new region number is assigned, and cells with the same type and 8-neighboring cells are recursively included in the same region. After labeling, the vector boundary polygon and area value of each independent damage region are extracted. The vector boundary polygon is obtained by polygon fitting the coordinates of the outer contour cells of each connected component, and the area value is equal to the total number of cells contained in the connected component multiplied by the ground area of ​​a single cell (0.25 square meters at a 0.5-meter spatial resolution).

[0063] Based on the mean confidence score of ecological mutation and the mean amplitude of feature deviation for each pixel within each independent damaged area, the degree of damage is divided into four levels: mild damage, moderate damage, severe damage, and extremely severe damage. In one specific implementation, a preset classification threshold is used for classification: areas with a mean confidence score of ecological mutation below 0.75 and a mean amplitude of feature deviation below 0.5 are classified as mild damage; areas with a mean confidence score of ecological mutation between 0.75 and 0.85 or a mean amplitude of feature deviation between 0.5 and 1.0 are classified as moderate damage; areas with a mean confidence score of ecological mutation between 0.85 and 0.95 or a mean amplitude of feature deviation between 1.0 and 2.0 are classified as severe damage; and areas with a mean confidence score of ecological mutation exceeding 0.95 and a mean amplitude of feature deviation exceeding 2.0 are classified as extremely severe damage. In one alternative implementation, the grading threshold can be adaptively determined based on regional statistical distribution. For example, quartiles can be calculated for the mean confidence level of ecological mutation and the mean deviation of characteristic values ​​for all independently damaged areas, and the 25%, 50%, and 75% quartiles can be used as the dividing points between mild, moderate, severe, and extremely severe conditions. This adaptive approach can avoid the problem of insufficient generalization of fixed thresholds across different mining areas.

[0064] The results of damage type determination (raster map, vector boundary polygons, area values, and damage severity level) are integrated and output into a structured ecological damage assessment report. The report is organized as follows: each independent damaged area serves as the basic recording unit. Each record includes the coordinate sequence of the vector boundary polygons for that area (using a projection coordinate system consistent with the original image), the damage type (minefield, spoil heap, tailings dam, exposed surface, or damaged vegetation), the damaged area (in square meters), the damage severity level (mild, moderate, severe, or extremely severe), and the temporal location information of the abrupt change. This structured report can be directly imported into geographic information system software for visualization and spatial querying, providing quantitative spatial reference for mine ecological restoration decisions.

[0065] While specific embodiments of the present invention have been described above, those skilled in the art should understand that these specific embodiments are merely illustrative. Those skilled in the art can omit, substitute, and modify the details of the above methods and systems in various ways without departing from the principles and essence of the present invention. For example, combining the above method steps to perform substantially the same function and achieve substantially the same result according to substantially the same method falls within the scope of the present invention. Therefore, the scope of the present invention is defined only by the appended claims.

Claims

1. An automatic method for determining the scope and extent of mine ecological damage based on multi-source data fusion, characterized in that, Includes the following steps: Step 1: Using an unmanned aerial vehicle (UAV) platform equipped with an airborne hyperspectral imager and a thermal infrared imager, the target mining area is repeatedly scanned and acquired at multiple different times. The acquired spectral images and thermal radiation images are preprocessed and registered to form a multi-temporal dual-mode registration dataset. Each pixel location in the multi-temporal dual-mode registration dataset contains a time series composed of a multi-dimensional spectral reflectance vector and surface temperature values ​​in chronological order. Step 2: For the spectral reflectance vector of each pixel location in each time phase, spectral features are encoded using a selective state-space model based on the Mamba architecture, outputting a spectral semantic representation vector; the spectral semantic representation vector is concatenated with the thermal feature vector to generate a spectral-thermal joint feature vector, forming a multi-temporal joint feature sequence; the multi-temporal joint feature sequence is input into the Bayesian online change point detection module to identify the moment of ecological state change, and pixels with ecological change confidence exceeding the preset change judgment threshold are marked as candidate pixels for ecological damage; Step 3: Based on evidence theory, perform fusion of spectral evidence and thermal infrared evidence on candidate pixels of ecological damage to determine the confidence distribution of damage type; perform spatial consistency optimization on the confidence distribution of damage type and then perform connected component analysis to extract independent damage regions and classify the degree of damage, and output a structured ecological damage assessment report.

2. The method according to claim 1, characterized in that, In step 1, the airborne hyperspectral imager acquires continuous spectral images covering the 400 nm to 2500 nm band, with no less than 256 spectral channels; the thermal infrared imager acquires thermal radiation images covering the 8 μm to 14 μm band; the number of different time phases is no less than 3, and the interval between adjacent time phases is no less than 7 days.

3. The method according to claim 1, characterized in that, In step 1, strip noise suppression, bad band removal, water vapor absorption band interpolation, and atmospheric radiation correction are sequentially performed on the continuous spectral images of each time phase to obtain a surface reflectance data cube; non-uniformity correction, radiometric calibration, and temperature inversion are performed on the thermal radiation images of each time phase to obtain a surface temperature distribution map.

4. The method according to claim 3, characterized in that, In step 1, the surface reflectance data cube of the first time phase is used as the reference image. After bilinear interpolation resampling is performed on the surface temperature distribution map, pixel-level registration is completed using a registration method based on maximizing mutual information. For each subsequent time phase, cross-time phase registration with the reference image of the first time phase is added to ensure that pixels of the same geographical location in all time phases have a strict correspondence.

5. The method according to claim 1, characterized in that, In step 2, the processing flow of the selective state-space model includes: setting the hidden state dimension, treating the spectral reflectance vector as a sequence of signals arranged along the band dimension and inputting them sequentially; for the input vector at each band position, generating a state transition modulation vector, an input gate vector, and an output gate vector through independent linear projection layers; the state transition modulation vector is processed by an activation function and then modulated element-wise with the state transition basis matrix to generate an adaptive state transition matrix; the hidden state vector is transformed by the adaptive state transition matrix and added to the input vector after element-wise multiplication by the input gate vector to obtain the updated hidden state vector; the updated hidden state vector is then multiplied element-wise with the output gate vector to obtain the spectral feature output vector; the spectral feature output vector at the last band position is used as the spectral semantic representation vector.

6. The method according to claim 1, characterized in that, In step 2, the execution flow of the Bayesian online change point detection module includes: initializing the running length to 0, setting the observation likelihood model to a multivariate Gaussian distribution, and setting the conjugate prior model to a normal inverse Wieshard distribution; for the spectral-thermal joint feature vector of each time phase arriving in chronological order, calculating the conditional likelihood value for the running segment continuation case and the marginal likelihood value for the change point occurrence case, and normalizing the conditional likelihood value and the marginal likelihood value to obtain the posterior probability of the change point occurrence in the current time phase; after traversing all time phases, selecting the maximum value of the posterior probability of the change point occurrence as the ecological mutation confidence, and recording the time phase position of the mutation and the Euclidean distance between the spectral-thermal joint feature vectors before and after the mutation as the feature deviation magnitude.

7. The method according to claim 1, characterized in that, In step 3, fusion is performed based on Dempster-Shafer evidence theory, specifically including: constructing a spectral evidence body based on the spectral semantic representation vector of the post-mutation phase, and outputting basic probability allocation values ​​for five damage types: mining area, spoil heap, tailings dam, exposed surface, and damaged vegetation; calculating the temperature anomaly amplitude based on the deviation of the surface temperature value of the post-mutation phase from the median surface temperature of the non-ecological damage candidate pixels in the local window, and constructing a thermal infrared evidence body based on the temperature anomaly amplitude; and performing orthogonal fusion of the spectral evidence body and the thermal infrared evidence body using Dempster combination rules to obtain the damage type confidence distribution.

8. The method according to claim 7, characterized in that, Before performing orthogonalization and fusion, the conflict factor between the spectral evidence and the thermal infrared evidence is calculated. If the conflict factor is lower than the preset conflict threshold, Dempster combination is performed directly. If the conflict factor is not lower than the preset conflict threshold, a discount factor that is negatively correlated with the conflict factor is applied to the thermal infrared evidence before Dempster combination is performed.

9. The method according to claim 1, characterized in that, In step 3, spatial consistency optimization is achieved using block-based conditional random field inference. The region consisting of candidate pixels of ecological damage is divided into overlapping sub-blocks of a preset size. A fully connected conditional random field model is constructed for each overlapping sub-block. The unary potential term of the fully connected conditional random field model is composed of the negative logarithm of the damage type confidence distribution, and the binary potential term of the fully connected conditional random field model is composed of the Gaussian kernel function of the spatial distance between pixels and the Euclidean distance of the spectral-thermal joint eigenvector. The mean field approximation method is used to iteratively infer each overlapping sub-block and then the mean value is taken in the overlapping region to output a raster map of the damage type determination result.

10. The method according to claim 1, characterized in that, In step 3, the raster map of the damage type determination result is labeled with 8-neighborhood connected components, and the vector boundary polygon and area value of each independent damage area are extracted. Based on the average confidence of ecological mutation and the average feature deviation of each pixel in each independent damage area, the degree of damage is divided into four levels: mild damage, moderate damage, severe damage and extremely severe damage, using a preset grading threshold or a grading threshold adaptively determined based on regional statistical distribution.

Images