Dynamic Range Extraction and Analysis Method for Damaged Geomorphic Features in Mining Areas

Abstract

This invention relates to the field of image analysis technology, and more specifically, to a method for extracting and analyzing the dynamic range of images of damaged landforms in mining areas. The method includes: Step 1: Integrating the multispectral radiance within a set wavelength range, while subtracting two atmospheric scattering terms at each wavelength point to obtain the corrected reflectance; Step 2: Multiplying the Euclidean distance by the curvature amplitude in the near-infrared band to obtain a spectral collapse-curvature index sensitive to browning / chlorosis; Step 3: Calculating the peak elastic strain energy density per unit area of ​​overburden-coal mass; Step 4: Extracting all linear cracks from remote sensing images of the same time period and statistically analyzing their orientation histograms, then using the Shannon entropy formula to calculate the anisotropic entropy modulation of the cracks, thereby constructing a comprehensive dynamic range damage index to characterize the degree of damage to each pixel. This invention achieves accurate identification and intensity classification of surface damage areas under mining disturbance in mining areas.

Description

Technical Field

[0001] This invention belongs to the field of image analysis technology, specifically relating to a method for dynamic range extraction and analysis of images of damaged landform features in mining areas. Background Technology

[0002] In recent years, with the continuous intensification of mining activities, surface damage has become increasingly common, forming complex and varied subsidence basins, crack networks, exposed coal dust areas, and browning zones in the overlying strata. These geomorphological changes not only threaten operational safety and the geological environment but also directly affect the stability of the ecosystem and the long-term sustainable development of mining areas. Therefore, accurate monitoring and dynamic assessment of damaged landforms in mining areas have become an important technical requirement in the fields of geological disaster prevention and control, ecological restoration, and mine supervision. Traditional methods for monitoring landform damage mainly rely on ground surveys and point displacement monitoring, such as deploying dense monitoring networks using the Global Positioning System (GNSS) and installing inclinometers and crack gauges in key areas. Although these methods are highly accurate, they are limited by factors such as spatial sparsity, high deployment costs, small coverage areas, and long cycles, making it difficult to comprehensively perceive the damage pattern in large-scale, rapidly changing mining environments.

[0003] With the development of remote sensing technology, dynamic monitoring methods for mining areas based on optical satellite imagery or interferometric radar (InSAR) have begun to emerge. For example, optical imagery change detection identifies bare land expansion, while InSAR time-series analysis captures surface subsidence rates. These methods overcome the limitations of ground-based monitoring, enabling large-scale, periodic, millimeter-level observation of surface deformation, and have become the mainstream technology for current mining area environmental monitoring. However, existing technologies mainly focus on a single physical dimension; for example, InSAR emphasizes deformation, and optical change detection focuses on surface cover changes. Neither fully explores the multi-source coupling characteristics of mining damage processes. For instance, the surface reflectance characteristics of subsidence areas can change significantly due to vegetation degradation, coal dust accumulation, and brown rock exposure, but simple detection of changes in the National Displacement Vegetation Index (NDVI) often mistakenly interprets seasonal yellowing or farmland harvesting as mining-induced damage, leading to numerous false alarms. Similarly, a high deformation rate does not necessarily indicate catastrophic surface damage, as some plastic deformation processes, such as slow subsidence, may not be accompanied by crack propagation. Therefore, a single data source or a single physical indicator cannot fully characterize the essential features of the damaged landforms in mining areas. The lack of a comprehensive description of dynamic range, gray-scale variability and landform structure organization has become a major shortcoming of existing technologies. Summary of the Invention

[0004] The main objective of this invention is to provide a method for extracting and analyzing the dynamic range of damaged landform features in mining areas. By fusing the spectral collapse curvature features of multispectral images, deformation energy indices obtained from InSAR or GNSS, and anisotropic entropy modulation of crack directions in remote sensing images, a unified comprehensive dynamic range damage index is constructed. This enables accurate identification and intensity classification of surface damage areas under mining disturbance. This method requires no empirical parameters, possesses traceability of physical quantities, and cross-platform versatility. It significantly improves the responsiveness and discrimination accuracy of remote sensing images in complex geomorphic environments, providing highly reliable technical support for mine disaster monitoring, geological risk assessment, and ecological restoration planning.

[0005] To solve the above problems, the technical solution of the present invention is implemented as follows:

[0006] A method for extracting and analyzing the dynamic range of images showing damaged landforms in mining areas, the method comprising:

[0007] Step 1: Integrate the multispectral radiance within the set wavelength range, and subtract two atmospheric scattering terms at each wavelength point to obtain the corrected reflectance;

[0008] Step 2: Based on the corrected reflectance, extract the center reflectance of the green, red, near-infrared, and short-wave infrared bands; calculate the difference between the center reflectance of the green band and the center reflectance of the red band to obtain the first reflectance difference; calculate the difference between the center reflectance of the near-infrared band and the center reflectance of the short-wave infrared band to obtain the second reflectance difference; then calculate the Euclidean distance between the first and second reflectance differences to capture spectral collapse, and multiply the Euclidean distance by the curvature amplitude of the near-infrared band to obtain the spectral collapse-curvature index, which is sensitive to browning / fading.

[0009] Step 3: Based on the Laplace operator and first-order gradient mode of instantaneous subsidence, the coupling of subsidence amplitude and deformation curvature is obtained. Combined with the measured density of the overburden, the peak value of elastic strain energy density per unit area of ​​overburden-coal rock mass is calculated.

[0010] Step 4: Extract all linear cracks from remote sensing images of the same time period and calculate their orientation histograms. Normalize the orientation frequency of each degree into a probability, and then use the Shannon entropy formula to calculate the anisotropic entropy modulation of the cracks. In this way, construct a comprehensive dynamic range damage index to characterize the degree of damage of each pixel.

[0011] Furthermore, in step 1, the wavelength range is set to 0.45μm to 2.20μm.

[0012] Furthermore, the corrected reflectivity ρ for wavelength λ n (λ) is:

[0013]

[0014] Among them, L TOA (λ) represents the multispectral radiance at wavelength λ, in W / m². -2 sr -1 μm -1 L path (λ) represents atmospheric path-scattered radiation with wavelength λ, in W / m². -2 sr -1 μm -1 ;τ d (λ) represents the aerosol optical thickness at wavelength λ, which is dimensionless; L sky (λ) represents the backscattered radiation from the sky at wavelength λ, measured in W / m². -2 sr -1 μm -1 ; sec represents the secant function; θ s θ is the solar zenith angle with wavelength λ, in rad. v The zenith angle is the observed wavelength λ, expressed in rad.

[0015] Furthermore, the center wavelength G in the green band is 0.56 μm; the center wavelength R in the red band is 0.66 μm; the center wavelength NIR in the near-infrared band is 0.86 μm; and the center wavelength SW in the short-wave infrared band is... 1.6 It is 1.6 μm.

[0016] Furthermore, the spectral collapse-curvature index (SCCI) is:

[0017]

[0018] Where, ρ n (R) is the center reflectivity of the red band; ρ n (G) represents the center reflectivity of the green band; ρ n (NIR) is the center reflectivity in the near-infrared band; ρ n (SW 1.6 () represents the center reflectivity of the shortwave infrared band; This represents the curvature amplitude in the near-infrared band.

[0019] Furthermore, the instantaneous subsidence s is calculated using InSAR or GNSS points from the same period.

[0020] Furthermore, the peak elastic strain energy density SCSED of pixel (x, y) is:

[0021]

[0022] Where, ρ rock ρ is the average density of the overlying rock; g is the acceleration due to gravity; δ is the Laplace operator; δ is the pixel resolution of InSAR or GNSS. Let be the derivative of the instantaneous sinking s along the x-axis; Let be the derivative of the instantaneous sinking s in the y-axis direction; Let be the first-order gradient magnitude of the instantaneous subsidence s.

[0023] Furthermore, in step 4, the Canny-Hough combined algorithm is used to extract all linear cracks from the remote sensing image at the same time period t and to calculate their orientation histograms. The orientation frequency of each degree is normalized into a probability. Specifically, this includes: using the Canny edge extraction operator to extract the binary edge map of the remote sensing image; then performing Hough line detection on the binary edge map to obtain the detected lines and their corresponding orientation angles; dividing the 0 to 180° interval into 180 orientation boxes, and then calculating which orientation box each line's orientation angle belongs to, thereby constructing an orientation histogram; based on the orientation histogram, the orientation frequency is normalized into a orientation probability P with a sum of 1. t (θ); The crack anisotropic entropy modulation quantity CAEM is:

[0024]

[0025] Furthermore, in step 4, the Composite Dynamic Range Degradation Index (CDRDI(x,y)) for pixel (x,y) is:

[0026]

[0027] in, NDVI(x,y) represents the mean of the normalized vegetation index (NDVI) for all pixels in time period t; NDVI(x,y) represents the normalized vegetation index for pixel (x,y); σ NDVI，t denoted as the standard deviation of the normalized vegetation index for all pixels in time period t.

[0028] The method for dynamic range extraction and analysis of mining area damaged landform feature images of this invention has the following beneficial effects: it can achieve high-precision and high-reliability identification of mining-disturbed areas on remote sensing images. By constructing a spectral collapse and curvature coupling index, it comprehensively captures spectral morphological changes caused by vegetation retreat, coal dust deposition, and brown rock exposure, enhancing the image's response to browning and damaged areas; simultaneously, it introduces a composite index of instantaneous subsidence deformation, spatial gradient, and curvature, elevating traditional deformation monitoring to a dynamic assessment from an energy perspective, solving the problem that existing methods cannot quantify the intensity of landform damage; furthermore, it adopts a linear crack direction statistics and entropy modulation mechanism to incorporate the orderliness of image texture into the dynamic range adjustment model, realizing automatic control of the geometric complexity of gray-scale anomaly areas. The final constructed comprehensive index does not rely on any empirical weights or training samples, possesses consistency in physical quantity dimensions, cross-source data adaptability, and engineering scenario interpretability, and can effectively identify real high-risk areas and exclude pseudo-change areas in complex mining environments. This invention significantly improves the quantification, automation, and robustness of remote sensing disaster analysis in mining areas, providing reliable technical support for geological disaster prevention and control, safe mining, and ecological restoration. Attached Figure Description

[0029] Figure 1 This is a schematic diagram of the method flow for extracting and analyzing the dynamic range of damaged landform feature images in mining areas, as provided in an embodiment of the present invention. Detailed Implementation

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

[0031] refer to Figure 1 A method for dynamic range extraction and analysis of images of damaged landforms in mining areas, the method comprising:

[0032] Step 1: Integrate the multispectral radiance within the set wavelength range, and subtract two atmospheric scattering terms at each wavelength point to obtain the corrected reflectance;

[0033] Specifically, the process begins with the top-level atmospheric radiance recorded by the sensor. However, the original radiance is simultaneously affected by tropospheric molecular scattering and aerosol scattering. The former is mainly caused by Rayleigh scattering of incident light by atmospheric molecules, resulting in a background increase across the entire wavelength band. The latter originates from aerosol particles formed by dust, coal dust, and emissions from mining operations, which significantly contribute to sky backscattering. To eliminate these two types of interference, path scattering and dust-related scattering terms need to be subtracted at each discrete wavelength sampling point. Path scattering is usually obtained through dark pixel estimation or radiative transfer simulation, while dust scattering is obtained by multiplying the aerosol optical thickness retrieved from the on-site solar photometer with the simulated sky radiation field. After scattering subtraction, to avoid anisotropic errors caused by differences in the incident and observation angles in surface reflection, a cosine factor of the solar zenith angle and sensor viewing angle needs to be introduced into each pixel to correct for the effective incident energy dilution effect caused by oblique illumination along with the oblique path increment.

[0034] Subsequently, the corrected radiance at each wavelength is numerically integrated according to the wavelength step to obtain the integrated energy reflecting the overall radiative response of the pixel. This energy is then normalized using the solar spectral constant for the total incident energy within the same wavelength band, ensuring radiative comparability across different dates, solar altitudes, and sensor images. Since slope aspect variations in the mining area's subsidence basin alter local incident angles, the angle correction term in the formula promptly compensates for differences in shadows and bright spots caused by micro-topography. Meanwhile, aerosol scattering subtraction significantly reduces the false brightness caused by blasting and transportation dust, resulting in more accurate reflection characteristics for exposed coal seams, browned bare rock, and dust-laden areas after correction. Through this comprehensive integration and dual-scattering removal process, the final output spectral corrected reflectance retains the true physical reflection differences of the mining area's topography while maximally suppressing atmospheric and geometric noise. This provides a unified, stable, and consistent radiative benchmark for subsequent spectral collapse detection, strain energy field fusion, and crack direction entropy modulation. The entire processing chain can be automatically executed on large-scale, high-resolution, and multi-temporal remote sensing images of mining areas, ensuring that the grayscale dynamic range response of damaged landforms under different observation conditions remains comparable, thereby improving the accuracy and temporal consistency of the subsequent dynamic range damage index in identifying complex damages such as mining-induced subsidence, coal dust browning, and vegetation degradation.

[0035] Step 2: Based on the corrected reflectance, extract the center reflectance of the green, red, near-infrared, and short-wave infrared bands; calculate the difference between the center reflectance of the green band and the center reflectance of the red band to obtain the first reflectance difference; calculate the difference between the center reflectance of the near-infrared band and the center reflectance of the short-wave infrared band to obtain the second reflectance difference; then calculate the Euclidean distance between the first and second reflectance differences to capture spectral collapse, and multiply the Euclidean distance by the curvature amplitude of the near-infrared band to obtain the spectral collapse-curvature index, which is sensitive to browning / fading.

[0036] After radiometric correction, all pixels have been converted into cross-band comparable surface reflectance gratings. The core task of the second step is to analyze the spectral collapse characteristics caused by mining damage from these corrected reflectances and quantify them using a comparable single index. First, it is necessary to accurately locate four physically significant central bands in the corrected image: the green band records the high chlorophyll reflectance peak of healthy vegetation, the red band is located in the strong chlorophyll absorption trough, the near-infrared band corresponds to the specular high reflectance region of vegetation cell structure, and the 16-micron shortwave infrared band is most sensitive to changes in soil and dust moisture content. After mining, vegetation usually withers or is covered by dust, causing the green peak to weaken sharply, while the red trough rises with browning and oxidation, thus narrowing the difference between the two. Simultaneously, the high reflectance of near-infrared decreases sharply when vegetation turns green or coal dust accumulates, while shortwave infrared shows changes in different directions due to water evaporation or bare soil exposure. Therefore, the difference between near-infrared and shortwave infrared also shows a step-like convergence.

[0037] To weave these two difference values ​​into a highly sensitive grayscale scale for degradation and browning in mining areas, the first and second differences need to be treated as two orthogonal components on a Cartesian coordinate system. By calculating their Euclidean distance, the two differences are mapped to a single amplitude, which can be understood as the length of the "collapsed" vector in spectral space. A larger value indicates a more severe deviation of the pixel's spectral morphology from native vegetation or undisturbed surface. However, relying solely on the difference length is insufficient to capture the sharp bends in the browning spectrum of mining areas. Therefore, it is necessary to further utilize the second-order morphological features of the local spectral curve in the near-infrared band to enhance the response to minor degradation. In engineering implementation, a multispectral reflectance sequence can be extracted in wavelength order. For several equally spaced bands before and after the near-infrared center point, the second-order curvature amplitude can be calculated using the five-point central difference or Savitzky-Golay smoothing derivative method. A larger amplitude indicates a sharper bend in the near-infrared curve, typically corresponding to a precipitous drop in vegetation health or rapid coal dust coverage. Next, the curvature amplitude is directly multiplied and fused with the aforementioned Euclidean distance. This is equivalent to considering both the overall spectral collapse scale and the intensity of local bending in the grayscale space, which amplifies the response of the extreme browning region and suppresses the false collapse signal caused by a small amount of noise.

[0038] To ensure robust stability during large-scale batch processing, precise registration of the actual center wavelengths of the four center bands is required to avoid center value deviations caused by slight sensor drift. Furthermore, noise removal, cloud masking, and signal-to-noise weight adjustment must be performed on the images before difference and curvature calculations. If different satellite platforms are used for the mining area images, the discrete spectra of each platform can be converted to a unified center wavelength through linear spectral resampling or cubic spline interpolation. The entire process is executed in parallel at the pixel level, ultimately generating a single-valued spectral collapse-curvature index raster for each pixel. This raster, expressed as grayscale intensity, directly participates in the subsequent comprehensive calculation of the elastic strain energy field and crack direction entropy. Because this index fully integrates the spectral morphological changes caused by vegetation degreening, browning, coal dust deposition, and abrupt changes in water content, which are the most typical surface responses to mining disturbances, it can provide extremely high discrimination in dynamic range mapping. This makes key areas such as the collapse center, the toe of the spoil heap, and the dust channel significantly brighter or darker on the grayscale map, laying a solid spectral foundation for the subsequent quantitative classification of damage levels.

[0039] Step 3: Based on the Laplace operator and first-order gradient mode of instantaneous subsidence, the coupling of subsidence amplitude and deformation curvature is obtained. Combined with the measured density of the overburden, the peak value of elastic strain energy density per unit area of ​​overburden-coal rock mass is calculated.

[0040] After completing the collapse detection using spectral information, it is necessary to link the three-dimensional deformation field and radiation dynamic range of the surface caused by mining activities in the mining area. The key to step 3 is to convert the instantaneous subsidence into an energy scale that can be compared with the spectral index, so as to reveal the coupling characteristics of the collapsed surface in both gray-scale and mechanical spaces. This is achieved by first obtaining millimeter-level vertical displacements from time-series InSAR or high-frequency GNSS, and then unifying the displacement raster to an absolute vertical frame through multi-track registration, track error correction, and reference area constraints. Then, a first-order central difference is performed on the instantaneous subsidence in the pixel grid to obtain the spatial gradient fields along the east-west and north-south directions. The gradient field reflects the slope aspect and the degree of inclination at various locations in the subsidence depression. A larger gradient means that the same depth of subsidence is completed within a short distance, indicating that there is strong differential deformation in the supporting coal pillars or overburden beams below the goaf. Then, a Laplace convolution operation is performed on the subsidence grid to obtain the second-order spatial curvature distribution. The high curvature area corresponds to the position with the smallest bending radius, which is often where the overburden is about to undergo shear cracking or has already developed tensile cracks. Multiplying these two mechanical quantities can integrate the displacement amplitude information and deformation bending information into the same pixel index, thus simultaneously reflecting the subsidence depth and spatial abrupt changes at a single value level.

[0041] To give this index a clear physical energy meaning, the volumetric density of the overburden in the mining area needs to be introduced. This density is usually obtained through core experiments or wellbore acoustic logging. Multiplying its value by the constant gravitational acceleration of the region can elevate the displacement-curvature combination to an approximate scale of elastic potential energy per unit area. During calculation, to reduce discrete difference noise, the gradient and curvature are smoothed with Gaussian weighted kernels before and after convolution. The kernel scale is set based on the main fault spacing or working face length revealed by seismic exploration. This preserves the large-scale depression framework caused by underground mining while suppressing high-frequency interference. For scenarios involving segmented or multi-layered mining of thick coal, the same pixel may have superimposed displacements at different depths. Therefore, the multi-layer displacements need to be weighted and averaged according to the overburden burial depth before being included in the gradient and curvature calculations to avoid shallow micro-motions masking deep main mining signals. The processed elastic strain energy grid has the same spatial resolution as the spectral collapse index and is perfectly aligned with it through a unified geographic projection and step size, thus achieving a one-to-one correspondence at the pixel level when constructing the subsequent comprehensive dynamic range damage index.

[0042] High-value areas are often concentrated in the center of the working face, the intersection of the hanging wall and footwall of the failure zone, and the load edge of the spoil heap. These areas not only show strong abrupt changes in brightness and darkness in the radiometric images, but also accumulate the maximum potential energy in the mechanical field. Therefore, this energy peak index provides a mechanical constraint for disaster damage classification. To enhance dynamic monitoring capabilities, incremental difference analysis can be performed on the displacement field within different observation periods to obtain the instantaneous elastic energy growth rate over short periods. Then, a rate sequence is generated using time sliding window filtering to identify the location where secondary collapse or accelerated slope instability is about to occur. The entire process can be automatically and streamlined in the cloud, realizing the energy interpretation of large-scale, long-term, uninterrupted deformation in the mining area.

[0043] Step 4: Extract all linear cracks from remote sensing images of the same time period and calculate their orientation histograms. Normalize the orientation frequency of each degree into a probability, and then use the Shannon entropy formula to calculate the anisotropic entropy modulation of the cracks. In this way, construct a comprehensive dynamic range damage index to characterize the degree of damage of each pixel.

[0044] In the specific implementation, noise suppression and tonal equalization are first performed on high-resolution remote sensing images acquired at the same time stamp as the deformation field to ensure that the faint cracks under coal dust cover still maintain sufficient gradient contrast. Then, the Canny edge operator is used to establish a gradient amplitude and orientation grid for the entire pixel area. False edges generated by gravel and vegetation gaps are filtered out using a double-threshold connection, retaining only thin line pixels with consecutive lengths exceeding the threshold and consistent gradient directions. After obtaining the binary edge map, a probabilistic Hough transform is applied to accumulate straight-line features in polar coordinate space, outputting the pixel start and end coordinates of each crack. Mapping these coordinates back to geographic space yields the azimuth angle of each crack. The interval from 0 to 180 degrees is then divided into one-degree resolution orientation boxes, and all cracks are counted based on their azimuth angle boxes. Since different crack lengths contribute differently to the stretching of the visual dynamic range, the total crack pixel length is used instead of the number of cracks as the frequency benchmark during statistics to ensure that long-scale cracks occupy a higher weight in the probability distribution. After summing the lengths of each direction and dividing by the total length of the cracks across the entire field, a set of directional probabilities is obtained. Then, using the Shannon entropy calculation method in information theory, the logarithm of each directional probability is taken and multiplied by itself. Finally, the sum is obtained by integrating over a range of zero to 180 degrees and then normalizing to the maximum theoretical entropy.

[0045] The higher the anisotropic entropy modulation value of the cracks, the more uniform the distribution of crack lengths in different directions and the more chaotic the texture. Conversely, if the entropy value is close to the lower limit, it indicates that the crack orientation is concentrated and orderly, and the image grayscale is easily enhanced overall in terms of structure. To allow this texture complexity to play a suppressive role in the dynamic range, its square root is placed in the denominator of the comprehensive index, so that the high-entropy region produces the effect of denominator amplification and dynamic range reduction, while the low-entropy region maintains or enhances the dynamic response. To avoid the distortion of direction statistics by false crack linear objects such as cloud shadows and bare rocks, it is necessary to use the spectral thresholding method to remove high-reflectivity bare rocks and the DSM local slope aspect to remove steep cliff shadows before crack extraction, and then use the vegetation elevation texture mask to exclude canopy gaps. During the repeated blasting period of the working face, images from multiple time periods can be collected and the crack orientation probability can be averaged to weaken the orientation noise of temporary shallow cracks by time filtering, while the persistent and stable tension cracks contribute significantly to the orientation entropy, thus reflecting the long-term dynamic risk in the index.

[0046] Furthermore, to accommodate the spatial resolution differences between different satellites or UAV platforms, the crack pixel length needs to be normalized according to the actual ground size to ensure that the directional entropy output from images at different scales is comparable. Finally, the crack anisotropic entropy modulation amount, spectral collapse curvature index, and elastic strain energy density peak are superimposed at the pixel level. After denominator modulation and smoothing threshold processing, a comprehensive dynamic range damage index raster is generated. The value of the index directly corresponds to high-risk areas such as the collapse center of the mining area, the dense crack zone, and the unstable edge of the spoil heap. This achieves cross-domain coupling from local texture to global dynamics, making the entire extraction and analysis method form a complete closed loop with spectral, deformation, and geometric triple constraints.

[0047] Furthermore, in step 1, the wavelength range is set to 0.45μm to 2.20μm.

[0048] A wavelength window of 0.45 μm to 2.20 μm serves as the fixed boundary for integration and correction in step 1. This range covers visible light (blue-green-red), near-infrared, and two typical short-wave infrared bands. It includes the chlorophyll absorption-reflection region, which is most sensitive to vegetation, coal dust, and browning, as well as the 1.6 μm and 2.2 μm channels, which are most sensitive to dust moisture and ore oxidation. Therefore, it can completely capture all key spectral responses after mining activities in the mining area. At the same time, it excludes the long-wave infrared band above 2.2 μm, which has a low signal-to-noise ratio and is significantly affected by thermal radiation, to avoid irrelevant energy introducing grayscale distortion. In actual implementation, it should be ensured that scattering subtraction, angle correction, and integration accumulation are performed on each discrete center wavelength within this wavelength window for all input images. Missing channels can be resampled by interpolation of adjacent bands or by convolution using the sensor response function, thereby ensuring the integrity and comparability of the integrated energy within 0.45–2.20 μm and laying a unified benchmark for subsequent spectral collapse-curvature detection.

[0049] Furthermore, the corrected reflectivity ρ for wavelength λ n (λ) is:

[0050]

[0051] Among them, L TOA (λ) represents the multispectral radiance at wavelength λ, in W / m². -2 sr -1 μm -1 L path (λ) represents atmospheric path-scattered radiation with wavelength λ, in W / m². -2 sr -1 μm -1 ;τ d (λ) represents the aerosol optical thickness at wavelength λ, which is dimensionless; L sky (λ) represents the backscattered radiation from the sky at wavelength λ, measured in W / m². -2 sr-1 μm -1 ; sec represents the secant function; θ s θ is the solar zenith angle with wavelength λ, in rad. v The zenith angle is the observed wavelength λ, expressed in rad.

[0052] ρ n (λ) unifies multispectral observations into normalized reflectance, which can directly characterize surface energy exchange. The upper and lower limits of the equation are set between 0.45 μm and 2.20 μm. The physical motivation behind this is to encompass all key channels in the visible-near-infrared-shortwave infrared spectrum most sensitive to vegetation degradation, coal dust browning, bare rock oxidation, and water content changes, while simultaneously shielding unrelated mixtures of surface and atmospheric radiation in the longer-wavelength thermal infrared bands. The L term first appears in the molecule. TOA (λ) is the radiance of the top atmosphere received by the sensor. It includes not only the actual reflection information of the target pixel, but also Rayleigh and Mie scattering along the path from the sun to the surface and then to the sensor. If directly used for grayscale analysis in mining areas, it will misinterpret the brightness of non-surface sources as brightness variations caused by collapse or dust, thus requiring two-stage subtraction. The first stage subtracts L. path (λ), removing uniform atmospheric path scattering; the second stage utilizes aerosol optical thickness τ. d (λ) multiplied by the sky backscattered radiation L sky (λ) specifically removes wavelength-dependent brightness increases caused by aerosols such as dust from mining areas and blasting. This step is particularly important in mining scenarios because dust from mined-out areas often causes near-infrared low-reflectance areas to be raised by sky scattering; if not removed, it will mask vegetation decay signals. Following this is secθ. s (λ)+secθ v The (λ) term belongs to the first-order anisotropic correction. In the subsidence basin or steep slope of the spoil heap in the mining area, both the local incident angle and the observation angle change, causing the same material to exhibit drastically different brightness on different slope directions. By taking the secant of the solar zenith angle and the sensor viewing angle and then adding them together, it is equivalent to pressing the path increment caused by oblique illumination and oblique viewing angle back to the radiation reference of the vertical viewing angle in a geometric proportion, so as to achieve radiation comparability between pixels on different slope directions, while preserving the dynamic contribution of real shadow and highlight areas in grayscale.

[0053] The integral operation sums the wavelength point brightness after scattering subtraction and angle correction within the range of 0.45–2.20 μm to obtain the total energy actually returned by the pixel. To eliminate the influence of differences in imaging date, Earth-Sun distance, and solar altitude on the energy reference, the denominator is adjusted relative to the solar constant irradiance E0(λ)cosθ within the same band window. sIntegrating (λ) serves as a normalization scale: pixels with the same reflectance, regardless of whether they are in winter images at low solar altitudes or summer images at high solar altitudes, are mapped to a consistent 0–1 range, ensuring direct comparison of subsequent cross-temporal spectral collapse indices. Through this numerator-denominator structure, ρ... n (λ) effectively compresses the three dimensions of sensor observation space, solar radiation space, and atmospheric scattering space into a single surface reflectance dimension. This eliminates external noise from dust and path scattering while providing endogenous compensation for local topographic geometry, enabling the extraction of spectral differences representing the true surface properties even under multi-source disturbances in mining areas. Since this normalized reflectance maintains energy conservation at the pixel scale, it becomes the direct input for subsequent calculations of the first reflectance difference, the second reflectance difference, and the near-infrared curvature amplitude. This ensures that the response to browning and chlorosis during the spectral collapse detection stage originates entirely from the surface physical state rather than atmospheric or geometric artifacts, thus ensuring that the final integrated dynamic range deterioration index possesses both radiometric rigor and accuracy.

[0054] Furthermore, the center wavelength G in the green band is 0.56 μm; the center wavelength R in the red band is 0.66 μm; the center wavelength NIR in the near-infrared band is 0.86 μm; and the center wavelength SW in the short-wave infrared band is... 1.6 It is 1.6 μm.

[0055] In the dynamic range extraction process of images of damaged landforms in mining areas, 0.56μm, 0.66μm, 0.86μm, and 1.6μm were selected as the center wavelengths for green, red, near-infrared, and shortwave infrared, respectively, to maximize the capture of four complementary sensitive zones caused by mining disturbances on the surface spectral morphology. The green band at 0.56μm is located at the high reflectance peak of the chlorophyll absorption curve of healthy vegetation. When mining road excavation, spoil heaps, or dust cover cause vegetation deveining, this peak is the first to decay, thus it can sensitively record the brightness collapse caused by the decline in vegetation biochemical function. The adjacent 0.66μm red band is in the main chlorophyll absorption trough. Normal vegetation should show low reflectance in this band; however, coal dust adhesion or exposed loess caused by mining reduces the absorption depth of the red band, causing its brightness to rise. Therefore, the difference between the weakening green peak and the rising red trough narrows abruptly, forming a typical spectral "collapse" sign. Moving further into the longer wavelength range to 0.86 μm, the near-infrared band corresponds to the specular scattering region within the cell walls of vegetation. Healthy canopies exhibit extremely high reflectivity in this band. However, once vegetation is removed or coal dust covers the area, the near-infrared reflectivity drops drastically. Simultaneously, exposed overburden or water accumulation in collapsed basins creates additional curvature bends in the near-infrared spectrum. Therefore, using this band as a reference point for curvature calculation significantly amplifies sharp bends in the spectral curve, demonstrating extremely high sensitivity to browning and chlorosis. Finally, the 1.6 μm short-wave infrared center band is introduced because this band is particularly sensitive to surface moisture content, ore oxidation, and dust particle size distribution. Weathered rocks in spoil heaps and the floor of mined-out areas exhibit strong absorption differences in this band during alternating wet and dry conditions or oxidative browning, making the difference between this band and the high reflectivity of the near-infrared band a crucial criterion for determining the exposure of bare soil and rock in mining areas. At the same time, coal dust particles absorb less of the short-wave infrared than the near-infrared, further weakening the contrast between the two bands in areas with thick coal dust accumulation, exacerbating spectral collapse. By constructing two sets of differences from these four center wavelengths and calculating the Euclidean distance, the combined effects of vegetation degradation, coal dust cover, oxidative browning, and water content changes can be merged into a single amplitude in two-dimensional spectral space. Multiplying this amplitude by the near-infrared curvature couples the overall brightness attenuation with local curvature convexity, thereby generating the spectral collapse-curvature index, which is most discriminative of post-mining surface browning and devegetation. This band configuration avoids redundant calculations caused by using too many channels and ensures that all critical mining area damage signals are fully expressed within the selected wavelength window, providing the most physically comprehensive and spectrally redundant input basis for subsequent fusion with deformation energy and fracture entropy modulation.

[0056] Furthermore, the spectral collapse-curvature index (SCCI) is:

[0057]

[0058] Where, ρ n (R) is the center reflectivity of the red band; ρn (G) represents the center reflectivity of the green band;

[0059] ρ n (NIR) is the center reflectivity in the near-infrared band; ρ n (SW 1.6 () represents the center reflectivity of the shortwave infrared band; This represents the curvature amplitude in the near-infrared band.

[0060] In the method for extracting and analyzing the dynamic range of images of damaged landforms in mining areas, the Spectral Collapse-Curvature Index (SCCI) is designed to simultaneously measure the overall shrinkage and local bending effects of vegetation degradation, coal dust cover, bare rock oxidation, and changes in moisture in spoil heaps on the surface spectral morphology. Its core formula consists of the multiplication of two parts. Part One The four central reflectivities, after atmospheric scattering and BRDF correction, are organized into two-dimensional Euclidean distances: ρ n (G) captures the chlorophyll peak at 0.56 μm, which rapidly decays once vegetation is de-greened or covered by coal dust; ρ n (R) is located in the 0.66 μm chlorophyll absorption groove. When coal dust or bare soil covers the leaf surface, the reflection of the red groove will increase; the absolute contraction of the difference between the two directly quantifies the intensity of vegetation damage. At the same time, ρ n (NIR) recorded high reflectance of cellular structures at 0.86 μm, which decreased significantly in vegetation death or browning, while ρ n (SW 1.6 The reflectance at 1.6 μm is highly sensitive to moisture content and oxidation level. Changes in the dryness and wetness of the spoil heap or the browning of rocks will cause independent fluctuations in the reflectance of this band. The difference between the two is orthogonal to the previous difference, which can compensate for the moisture and oxidation information missed by a single difference. When the disturbance in the mining area causes the green-red and near-infrared-shortwave infrared differences to narrow simultaneously, the Euclidean distance becomes significantly smaller, which is equivalent to forming a "spectral collapse" vector in the four-dimensional gray space. When the difference widens or changes only in one direction, the distance reflects other environmental or background factors, ensuring that the index maintains separation among multiple spectral forms.

[0061] Part Two The second derivative amplitude of the near-infrared center band is used to quantitatively describe the curvature intensity of the spectral curve at the boundary between high vegetation reflectance and water absorption. When the vegetation is healthy, there is a smooth slope between the near-infrared and short-wave infrared bands. When browning or coal dust covers the area rapidly, the curve drops sharply from high reflectance to low reflectance, and the absolute value of the second derivative increases dramatically. This jump often precedes the green-red difference and can be detected in the early stages of deveining. In addition, when the soil moisture in the spoil heap evaporates or absorbs moisture rapidly in a short period of time, local peaks and valleys are introduced near this band, further amplifying the curvature. The design of multiplying the curvature amplitude by the Euclidean distance allows the SCCI to maintain a large Euclidean distance in healthy vegetation areas while remaining neutral overall due to the near-zero curvature. However, in areas such as the center of subsidence, thick coal dust accumulation, and spoil heap slopes, the Euclidean distance tends to converge while the curvature amplifies, resulting in a significant decrease or increase in the product output, creating a bidirectional stretching of the dynamic range. This two-factor coupling avoids the excessive amplification of shadows or sensor noise by a single difference and suppresses spurious peaks in curvature under noisy conditions. More importantly, SCCI, through two-dimensional differences with complementary physical meanings and second-order spectral information, synchronously compresses the radiation intensity jumps and curve bends caused by mining activities into a single-value index at the imaging pixel scale. This allows the dynamic range of the spectral dimension to be fully mapped into the grayscale space, providing a highly discriminative and numerically stable input for subsequent fusion with deformation energy fields and fracture direction entropy. Within the entire mining area, high SCCI values ​​typically correspond to edge areas where vegetation still exists but is affected by overburden stress, median values ​​correspond to undisturbed background surfaces, and low values ​​precisely fall in the core areas where coal dust coverage is most severe or subsidence and water accumulation are most significant. This gradient distribution spatially couples with the InSAR displacement gradient zone and the dense fracture zone, enabling the comprehensive dynamic range damage index to achieve pixel-level closed-loop correspondence in the multi-source information space.

[0062] Furthermore, the instantaneous subsidence s is calculated using InSAR or GNSS points from the same period.

[0063] Instantaneous subsidence *s* serves as a bridge connecting radiometric grayscale with the mechanical field, and its determination depends on InSAR or high-frequency GNSS points from the same period. For InSAR, the core mechanism utilizes the phase difference captured by repeated satellite orbit imaging; that is, the interference fringes obtained after interfering two complex radar images record subtle changes in the round-trip path from the satellite to the Earth's surface. Through phase unwrapping, the number of fringes can be converted into displacement along the line-of-sight. Vertical deformation caused by mining activities dominates, so given the radar incident angle, the line-of-sight displacement can be restored to an approximate vertical subsidence using simple geometric projection. To obtain millimeter-level accuracy and eliminate atmospheric delay and orbit control errors, a permanent or distributed scatterer temporal interferometry method is typically used. After applying stable reference area constraints to multiple images from the entire year or multiple months, the transient line-of-sight displacement for each period is jointly solved, and then differencing to a specified observation date to form a subsidence field grid for the same period. If using continuous GNSS stations, the principle is to deploy multiple high-precision receivers on the ground and use carrier phase observation equations to calculate the three-dimensional coordinates of the station center in real time. When the working face advances or the overlying rock becomes unstable, the vertical coordinates of the station will shift downwards by centimeters to decimeters. The instantaneous vertical displacement can be obtained by subtracting the coordinates of the station during the target period from the coordinates during the reference period. GNSS can provide an absolute displacement reference to correct potential vertical offsets in InSAR; while InSAR compensates for the sparseness of GNSS points with dense spatial sampling. By combining the two, the true GNSS values ​​can be embedded into the radar displacement plane through least squares adjustment to achieve a high-precision, high-resolution instantaneous subsidence field. To ensure a close correspondence with multispectral imagery, radar image orbital groups from the same day or week as the optical imagery should be selected, or coordinate solutions with a time difference of less than a few hours from the image acquisition time should be extracted from the continuous GNSS time series. This ensures that radiance and mechanical subsidence are synchronized in the time dimension, thus preventing time mismatch when multiplying pixel-level elastic strain energy and spectral collapse index in subsequent calculations. The subsidence surfaces in mining areas typically exhibit an elliptic parabolic distribution. The S-grid calculated using InSAR or GNSS not only contains the depression depth but also preserves subtle transitions in the boundary gradient. In the next step, the Laplace operator can be used to extract bending stress concentration zones, which, together with the first-order gradient mode, constitute the subsidence amplitude and curvature coupling quantity, providing reliable source data for elastic strain energy density. The entire deformation solution chain is based on radar phase physical delay or carrier phase integer deviation, allowing for traceable and quantifiable errors. This ensures that the geomorphic dynamics information obtained from the comprehensive dynamic range damage index is engineering-verifiable.

[0064] Furthermore, the peak elastic strain energy density SCSED of pixel (x, y) is:

[0065]

[0066] Where, ρ rock ρ is the average density of the overlying rock; g is the acceleration due to gravity; δ is the Laplace operator; δ is the pixel resolution of InSAR or GNSS. Let be the derivative of the instantaneous sinking s along the x-axis; Let be the derivative of the instantaneous sinking s in the y-axis direction; Let be the first-order gradient magnitude of the instantaneous subsidence s.

[0067] The peak elastic strain energy density, by coupling the amplitude, curvature, and spatial gradient terms of the instantaneous subsidence field at the same pixel scale, achieves a quantitative characterization of the potential elastic potential energy accumulation of the overburden-coal-rock mass. The formula's first factor ρ... rock g / 2 originates from the classical expression for the potential energy density of a deformable element, where ρ rock Obtained from core density measurements or wellbore acoustic inversion, it reflects the average unit volume weight of the overburden. g is a gravitational constant used to convert pure geometric displacement into gravitational potential energy, giving the index an energy dimension. The double coefficient in the denominator comes from the symmetry between elastic strain energy density and displacement Coulomb potential energy under the small strain assumption, ensuring dimensional closure of the amplitude. Two terms within the absolute value reflect the coupling between amplitude and curvature: s itself records the vertical settlement depth; the deeper the depth, the larger the overburden detachment space and the more severe the stress redistribution. Multiplying the displacement Laplace term by the square of the pixel resolution essentially converts the second-order spatial curvature of the pixel grid into an additional settlement with the same dimensions as the displacement amplitude. A higher absolute value of curvature indicates a smaller bending radius of the settlement surface and a more concentrated local tensile or compressive stress. Multiplying by δ... 2 This avoids artificially inflated curvature caused by grid refinement. By taking the absolute value of the outer layer, both the negative curvature depression in the center of the collapse basin and the positive curvature uplift of the spoil heap can be accumulated in the index in the form of positive energy, which is consistent with the fact that elastic energy can be stored in both the tension and compression zones in dynamics.

[0068] The gradient modulus term that follows This reflects the rate of spatial change of surface subsidence within a horizontal plane, specifically the steepness of the depression boundary. A steeper slope indicates a more significant difference in displacement between adjacent pixels, making it easier for shear or tensile stress to concentrate in the slope break zone. Therefore, it acts as a multiplicative amplifier in the index, resulting in a higher energy estimate at the collapse edge compared to the gently stressed central area. The entire expression projects the amplitude, curvature, and slope of the three-dimensional subsidence morphology onto a scalar quantity at the pixel level, with units of Jm. -2This index, along with the spectral collapse-curvature index of the same resolution, can be incorporated into the molecule of the comprehensive dynamic range damage index, injecting geomorphic mechanical intensity information into the brightness signal of the grayscale dynamic range. The formula structure also considers resolution portability: when using InSAR data from different orbits or GNSS data with different station densities, simply updating δ and maintaining the differential kernel scale ensures that the pseudo-gain caused by grid refinement on the curvature term is eliminated at the sub-pixel level, guaranteeing consistency in energy estimation across data sources. High-value areas output by this index typically correspond to the advanced arch curve in the working face advance direction, the shear ring at the collapse center, and the stress-induced area at the foot of the spoil heap in actual mining areas. These areas often simultaneously exhibit spectral collapse and dense crack textures in multispectral images, thus achieving a high degree of synergy between the mechanical and radiation fields when the comprehensive index is superimposed, improving the identification accuracy and dynamic range stretching effect of high-risk damaged units.

[0069] Furthermore, in step 4, the Canny-Hough combined algorithm is used to extract all linear cracks from the remote sensing image at the same time period t and to calculate their orientation histograms. The orientation frequency of each degree is normalized into a probability. Specifically, this includes: using the Canny edge extraction operator to extract the binary edge map of the remote sensing image; then performing Hough line detection on the binary edge map to obtain the detected lines and their corresponding orientation angles; dividing the 0 to 180° interval into 180 orientation boxes, and then calculating which orientation box each line's orientation angle belongs to, thereby constructing an orientation histogram; based on the orientation histogram, the orientation frequency is normalized into a orientation probability P with a sum of 1. t (θ); The crack anisotropic entropy modulation quantity CAEM is:

[0070]

[0071] On the one hand, roof fractures in goaf areas and instability in spoil heaps create tension cracks or shear cracks with predetermined orientations on the surface. These linear elements form brightness sharpening bands in high-resolution reflectance images. Once their directions converge, it's equivalent to superimposing high-amplitude, low-entropy strip signals on the spatial spectrum, causing the local grayscale distribution to shift from a Gaussian state to a peak state. On the other hand, when overburden enters late-stage waterlogging and collapse, or when spoil heap slopes undergo multiple rounds of sliding, cracks often expand randomly in a radial, fan-shaped, or grid-like pattern. At this point, the linear texture energy is diluted by the directional dimension, and the isotropic background still dominates. Even if image enhancement algorithms increase contrast, they will only result in mean shifts without regional brightness jumps. To allow the comprehensive dynamic range damage index to adaptively adjust between these two extremes, the entropy of the directional probability distribution must be used to constrain the final stretching amplitude of the grayscale dynamic range.

[0072] First, spectral equalization and median filtering are applied to the original image to eliminate local high-frequency noise. Then, the Canny edge operator is applied to extract significant gray-level gradient bands. Canny's dual-threshold connection mechanism can eliminate speckle pseudo-edges caused by uneven coal dust particle size while ensuring the continuity of crack pixels. The resulting binary edge map provides reliable pixel candidates for subsequent line recognition. Next, a probabilistic Hough transform is applied in polar coordinate space. The accumulator maps the possible line parameters of each pixel to a (ρ, θ) peak value. After thresholding, backtracking the pixel coordinates can separate all potential crack segments and obtain their orientation angle θ. Since the crack orientation does not distinguish between positive and negative, it is only necessary to take the modulus of the azimuth angle to the 0°~180° range to cover the entire orientation. To ensure that the contribution of long cracks to the orientation statistics is proportional to their amplification effect on the visual dynamic range, the pixel length is used to accumulate the orientation frequency instead of the number of cracks. The 0°~180° range is divided into 180 orientation boxes with a width of one degree. The length frequency f of all cracks is counted according to the orientation of each orientation box. t (θ), which is then normalized to the direction probability. The entropy measurement section introduces the squared weight θ of the direction angle. 2 This is because, in imaging geometry, vertical or large-angle cracks amplify the differences in local shadows and illumination more significantly, and have a greater visual impact on the grayscale dynamic range. If the crack direction is concentrated in the high-angle region, the upliftability of the geometric texture increases; conversely, if the distribution is uniform, the anisotropy is diluted, making it difficult to further enhance the contrast. Based on information theory, the higher the degree of randomness, the greater the entropy value. The classic form of directional entropy is -∑Pln P, which is multiplied by θ here. 2 Introducing the second moment of the angle, the final entropy is written as Considering that the maximum entropy occurs in a perfectly uniform distribution, the theoretical value is ln180°, using Normalization is performed to compress the entropy value into the range of 0 and 1, and a constant 1 is added to ensure that the index is always greater than 1, thus obtaining the anisotropic entropy modulation amount of the crack.

[0073] This index has three physical implications: First, if the fracture orientation is singular, such as a banded tensile fracture dominated by the main fault, then the majority of the probability is concentrated in a few directional boxes, the entropy term approaches zero, the CAEM approaches 1, and the comprehensive dynamic range damage index is close to zero in the denominator. The first form appears, thus minimizing the suppression of grayscale dynamic range, and further amplifying the ordered crack texture; secondly, if the crack direction is highly discrete, exhibiting a radial or checkerboard distribution, then P t (θ) is close to uniform, the entropy term approaches ln180°, CAEM approaches 2, and the increase in the denominator significantly weakens the dynamic response of grayscale, avoiding the false amplification of noise peaks caused by messy textures; thirdly, through θ 2After weighting, even with the same entropy value, high-angle cracks have a larger weight, thus vertical cracks generated on steep slopes at the collapse edge will have higher geometrically reinforced properties than low-angle cracks generated in the gentler areas of the spoil heap. It is worth noting that in real-world mining scenarios, cloud shadows, bare rock joints, or road edges are often misidentified as cracks. Therefore, a multispectral NDVI threshold and DSM slope are fused before Canny to filter out vegetation textures and abrupt changes in cliff elevation. A minimum length threshold is then set in the Hough results to remove short-distance pseudo-straight lines, ensuring that the directional probability truly represents the overlying tensile crack network. For multi-temporal images, to eliminate small crack noise induced by periodic road maintenance or temporary drainage ditches, P... t (θ) is used to perform a time-moving average, and then input into the entropy formula to obtain a robust directional chaos curve, thus reflecting the gradual collapse evolution rather than instantaneous noise in the comprehensive dynamic range damage index. The entire process is based entirely on the original image pixels and mathematical statistics, requiring no empirical parameters or training samples, and is applicable to satellite and UAV data with different spatial resolutions; cross-platform consistency can be guaranteed simply by updating the minimum line length of Hough detection and the θ sampling accuracy according to the resolution. Verification through numerous mining area cases has shown that CAEM is highly correlated with the anisotropic extensibility of the subsidence basin profile and the control effect of the main fault: In the high-tension working face of open-pit mines, early cracks develop rapidly in a main direction, and CAEM remains at around 1.05, with minimal suppression of grayscale enhancement, and the comprehensive index forms a high brightness wall at the subsidence boundary; as the working face expands, boundary cracks begin to derive in the vertical direction, CAEM rises to 1.4, and the brightness peak relatively declines, reflecting the weakening of texture enhancement by the tension-shear conversion; while in the subsidence area of ​​the spoil heap, cracks spread disorderly in a honeycomb pattern, CAEM rises to 1.8, the grayscale dynamic range is significantly compressed, and false indications of high disaster damage are avoided. It is evident that the crack anisotropic entropy modulation not only provides automatic weighting of the geometric texture dimension for the comprehensive dynamic range damage index, but also achieves adaptive geometric adjustment of the spectral collapse and deformation energy field of the mining area through a strict probability-entropy-angle second-order moment framework. This ensures that the gray-scale dynamic expression and the actual landform evolution maintain a high degree of consistency in time and space, thereby promoting a new engineering paradigm for monitoring surface disaster damage in mining areas from single-source radiation to multi-source coupling and high dynamic range quantification.

[0074] Furthermore, in step 4, the Composite Dynamic Range Degradation Index (CDRDI(x,y)) for pixel (x,y) is:

[0075]

[0076] in, NDVI(x,y) represents the mean of the normalized vegetation index (NDVI) for all pixels in time period t; NDVI(x,y) represents the normalized vegetation index for pixel (x,y); σ NDVI，t denoted as the standard deviation of the normalized vegetation index for all pixels in time period t.

[0077] The ln[1+SCSED] in the molecule first performs natural logarithmic compression on the peak energy density obtained from the coupling of instantaneous sinking amplitude, curvature, and gradient. This logarithmic mapping has two physical implications: First, elastic strain energy increases exponentially with sinking depth; directly using the original value would lead to energy dominating the synthesis in the extremely deep collapse zone, thus masking the moderate damage zone. Second, after taking the logarithm of the exponential change, the large and small energy regions are compressed back into a comparable range, achieving a "logarithmic stretching" of the dynamic range rather than nonlinear collapse. Meanwhile, [SCCI] 2 By squaring the product of the spectral collapse vector length and the near-infrared curvature amplitude, pixels with slight Euclidean distance contraction or curvature increase are geometrically magnified in the exponent, highlighting subtle changes in browning, fading, and coal dust coverage. Multiplying the two results in a numerator that encompasses both the severity of collapse at the mechanical potential level and dramatic jumps in spectral morphology at the radiation level, allowing truly high-risk pixels to exhibit a brightness potential far exceeding the background in the product domain. The denominator employs a two-term superimposed weighting structure: The square root of the crack entropy modulation is taken because an excessively high original entropy value would overly suppress spectral and energy information, while the square root can preserve the entropy's ability to distinguish the degree of texture chaos while avoiding denominator explosion. When the crack direction is highly dispersed, CAEM approaches the upper limit, and the increase in denominator weakens the collapse-curvature product, making it difficult to reflect the brightness enhancement of disordered texture regions. When the crack direction is uniform, CAEM approaches 1, and the denominator hardly suppresses the dynamic range, allowing the high grayscale potential of ordered crack bands to be fully released.

[0078] Second item Introducing the z-score square of NDVI has the physical meaning that when the vegetation index of a pixel is significantly higher than the mean for the same period, it is highly likely that it is a high reflectance area of ​​natural forest or farmland, rather than a brightness anomaly caused by mining damage. Adding the square of the deviation to the denominator will exert additional suppression on these high NDVI pixels, preventing the misidentification of naturally high-brightness areas as high-risk areas of damage against the backdrop of hilly woodlands or river valley oases. Conversely, if the NDVI of a pixel is much lower than the mean, it means that the vegetation has de-greened or is covered by coal dust; on the other hand, the z-score is negative, and its square is also positive, but the direction of deviation does not affect the magnitude of the suppression, thus maintaining the dynamic range response to degraded vegetation areas. The standard deviation σ of this term... NDVI，tBy adopting a time-based overall NDVI distribution, the differences in greenness across seasons are eliminated, allowing the NDVI deviation between winter deserts and summer oases to be assessed on the same scale. This also ensures that vegetation confidence degradation caused by mining activities within the mining area is consistently mapped over time. In summary, the algebraic construction of CDRDI follows three principles: First, the numerator uses multiplicative coupling to exponentially amplify high values ​​in both spectral and mechanical parameters, strengthening pixels representing true risks. Second, the denominator employs a dual weighting of interpretable entropy root terms and z-score square terms, automatically reducing the dynamic contribution of areas with chaotic textures or naturally high vegetation. Third, all input quantities are calculated from physically measurable indicators from previous steps, without any parameter tuning coefficients, ensuring the index remains comparable in absolute values ​​across satellite platforms, seasons, and even mining areas. Field measurements verified that in the subsidence area of ​​the inclined shaft, although the SCSED value of the large displacement area was high, the crack direction was concentrated, and the CAEM was as low as 1.1. In addition, the sparse vegetation and low NDVI resulted in a CDRDI as high as 18, which was significantly different from the background area of ​​2–4. In the open-pit spoil heap, although there was a high curvature curve, the crack direction was extremely chaotic, and the entropy value approached 1.9. In addition, the local herbaceous cover in the spoil heap resulted in a high NDVI, and the CDRDI was ultimately controlled at 6–8, avoiding misjudgment of large areas of medium and low risk areas. It can be seen that the CDRDI, which includes four transformations of logarithmic compression, square amplification, entropy root suppression and NDVI normalization, not only mathematically ensures the numerical stability of the gray-scale dynamic range, but also physically integrates the four major mechanisms of spectral collapse, subsidence energy, crack direction and vegetation background in the mining area. It provides a highly robust quantitative tool for large-scale mining area disaster monitoring that can be directly thresholded and consistent across sources.

[0079] We acquired local Sentinel-2 L1C multispectral imagery, the Sentinel-1 coherent permanent scatterer InSAR deformation sequence with the same orbital combination, and 0.15m resolution visible light aerial imagery acquired by UAV on the same day. To ensure cross-source consistency, all data were resampled to a δ=10m grid and projected onto WGS84 / UTM50N. The following demonstrates the calculation process of all indices using the pixel corresponding to the grid center coordinates (x0, y0). First, we performed band discrete integration on the radiometrically calibrated Sentinel-2 top atmospheric radiance within a band window from 0.45μm to 2.20μm, accumulating the term [L... TOA (λ)-L path (λ)-τ d (λ)L sky (λ)][secθ s (λ)+secθ v [(λ)] summed with a step size of 10 nm yields 1.92 × 10⁻⁶ 2 Wm -2 sr -1Simultaneously, using the ASTM E490 solar spectral constant and combining metadata, the solar zenith angle is used to represent ∫E0(λ)cosθ. s The integral of (λ)dλ yields 3.75 × 10⁻⁶. 2 Wm -2 Therefore, the normalized reflectance of this pixel is calculated as ρ. n =0.512.

[0080] Further readings were taken of the reflectivity of the four center bands: G (0.56μm) channel ρ n (G) = 0.42, R (0.66μm) channel ρ n (R) = 0.29, NIR (0.86μm) channel ρ n (NIR) = 0.34, SW 1.6 (1.6μm) channel ρ n (SW 1.6 = 0.27. Euclidean distance of the spectral collapse vector. To determine the curvature amplitude, five discrete reflectivities {0.29, 0.38, 0.34, 0.30, 0.28} are extracted in the range of 0.66 μm to 1.06 μm. These values ​​are then substituted into the second-order difference formula at the five-point center with a step size h = 0.1 μm. Therefore, the spectral collapse-curvature index SCCI of this pixel is 0.157(1+0.48) = 0.232. Subsequently, the line-of-sight displacement of -0.056m on the same day was extracted from the InSAR temporal difference field, and the instantaneous vertical subsidence s = -0.082m was obtained by projecting it through an incident angle of 38°.

[0081] A discrete Laplacian kernel [0, 1, 0; 1, -4, 1; 0, 1, 0] / δ is applied over a 3×3 pixel window. 2 Convolution yields Simultaneously, the horizontal gradient is obtained from the central difference. The mining geological exploration report gives the average density ρ of the limestone-mudstone overburden. rock =2.45×10 3 kg,m^{-3}. Substituting this into the formula for peak elastic strain energy density. in The resulting |-0.082+(-0.035)|=0.117m, gradient modulus Multiply by 2.45 × 10 3 Multiplying by 9.806 / 2 yields SCSED = 5.01 Jm -2 .

[0082] Next, using the 0.15m image, Canny extraction with σ=1.2 was performed in OpenCV to obtain 22261 edge lines. Probabilistic Hough detection, with a minimum line length of 50 pixels and a maximum discontinuity of 3 pixels, identified 143 cracks with a total length of 2348m. The cracks were concentrated in the directions of 85° (110m), 88° (950m), and 91° (670m), with the remaining length of 618m evenly distributed across the other 177 orientation boxes. Normalization yielded P... t (85°) = 0.047, P t (88°) = 0.405, P t (91°) = 0.285, and approximately 0.00148 for each of the other orientation boxes. Summation with second-order moment entropy. After normalization, CAEM = 1 + 2.86 / ln180 = 1.73. Vegetation index statistics were performed using NDVI raster data calculated on the same day using Sentinel-2B4 and B8 radii, with the average value over the time period. Standard deviation σ NDVI，t =0.12, target pixel NDVI 0.16, calculate z-score squared [(0.16-0.31) / 0.12] 2 =1.562.

[0083] Finally, the comprehensive dynamic range damage index is written as: Compared with the background pixels in the same image that are far from the mining area, the average CDRDI = 0.004 and the standard deviation of 0.006, the confidence level of this pixel is z = (0.0334-0.004) / 0.006 = 4.9, which is significantly higher than the 3σ threshold. It is automatically classified as the highest level of the three-level disaster damage. After visualization rendering, this pixel and its surrounding area form a brightness focused patch, which coincides with the settlement gradient zone and the unidirectional crack bundle. The crack width is 10-15cm and the vertical displacement of the slope is 9cm after on-site RTK re-measurement. The comprehensive index accurately reveals the high-risk slope section.

[0084] The 3σ threshold is a commonly used criterion derived from the normal distribution: if the value of a statistic falls outside the range of 3 standard deviations σ above and below the population mean μ, then that value is considered "extreme" or "abnormal". Expressed as a formula, the threshold range is μ ± 3σ. In the case of a perfectly normal distribution, approximately 68.27% of the sample is within μ ± 1σ, approximately 95.45% is within μ ± 2σ, and μ ± 3σ covers 99.73%. Therefore, the probability of falling outside this range is only about 0.27%, which is statistically considered a highly abnormal or significant event.

[0085] In the application of the comprehensive dynamic range deterioration index, the CDRDI of the background pixels of the entire mining area is regarded as a set of approximately normal reference distributions, and its mean is used. and standard deviation σ C Construct threshold If a pixel's CDRDI is greater than the upper limit, it means it is three standard deviations higher than the mining area average, falling at the extreme tail of 0.13%, and is very likely to correspond to high-risk damage areas such as collapse centers, densely fractured zones, or unstable edges of spoil heaps; conversely, if it is lower than the upper limit... In the context of damage analysis, such values ​​typically lack physical meaning and can be treated as abnormally low values ​​or even noise to be removed. By using a 3σ threshold, significant damage boundaries can be automatically drawn for the entire image without the need for manual empirical coefficients, achieving objective and repeatable level determination.

[0086] The above-described embodiments are only used to illustrate the technical solutions of the present invention, and are not intended to limit it. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features. Such modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for extracting and analyzing dynamic range of mine damage topographic feature images, characterized in that, The method comprises: Step 1: Integrate the brightness of multispectral radiation in a set wavelength range, and deduct two atmospheric scattering terms at each wavelength point to obtain a corrected reflectance; Step 2: Based on the corrected reflectance, extract the center reflectance of the green, red, near-infrared and short-wave infrared bands; calculate the difference between the center reflectance of the green band and the center reflectance of the red band to obtain a first reflectance difference; calculate the difference between the center reflectance of the near-infrared band and the center reflectance of the short-wave infrared band to obtain a second reflectance difference; then calculate the Euclidean distance of the first reflectance difference and the second reflectance difference to capture spectral collapse, and multiply the Euclidean distance by the curvature amplitude of the near-infrared band to obtain a spectral collapse-curvature index sensitive to browning / de-greening; Step 3: Based on the Laplacian operator and the first-order gradient module of the instantaneous subsidence, the coupling of the subsidence amplitude and the deformation curvature is obtained, and combined with the measured density of the overburden rock, the elastic strain energy density peak value of the unit area overburden-rock body is calculated; Step 4: On the remote sensing image of the same period, extract all linear cracks and count the direction histogram, normalize the direction frequency of each degree to probability, then use the Shannon entropy formula to calculate the crack anisotropy entropy modulation quantity, and build a comprehensive dynamic range damage index to represent the damage degree of each pixel. Wavelength of the corrected reflectivity is: ； where is the multi-spectral radiance at wavelength , in units of ; is the atmospheric path radiance at wavelength , in units of ; is the aerosol optical depth at wavelength , dimensionless; is the sky backscatter radiance at wavelength , in units of ; denotes the secant function; is the solar zenith angle at wavelength , in units of rad; is the observation zenith angle at wavelength , in units of rad; Pixel of the elastic strain energy density peak is: ； wherein, is the average density of the overburden; is the acceleration of gravity; is the Laplace operator; is the pixel resolution of InSAR or GNSS; is the instantaneous subsidence in the axial direction; is the instantaneous subsidence in the axial direction; is the instantaneous subsidence of the first gradient modulus; In step 4, the pixel of the integrated dynamic range impairment index is: ； wherein, is the mean of the normalized difference vegetation index of all pixels of the time period ; is the normalized difference vegetation index of the pixel ; is the standard deviation of the normalized difference vegetation index of all pixels of the time period .

2. The method of claim 1, wherein the method further comprises: In step 1, the set wavelength range is 0.45 µm to 2.20 µm.

3. The image dynamic range extraction and analysis method of claim 2, wherein, the center wavelength of the green band is 0.56 pm; the center wavelength of the red band is 0.66 pm; the center wavelength of the near infrared band is 0.86 pm; the center wavelength of the short-wave infrared band is 1.6 pm.

4. The image dynamic range extraction and analysis method of claim 3, wherein, Spectral collapse-curvature index is: ； wherein, Rred is the center reflectance for the red band; Rgreen is the center reflectance for the green band; Rnear-ir is the center reflectance for the near-infrared band; Rswir is the center reflectance for the short-wave infrared band; Cnear-ir represents the curvature amplitude for the near-infrared band.

5. The image dynamic range extraction and analysis method of claim 4, wherein, The instantaneous subsidence is calculated by InSAR point or GNSS point in the same period .

6. The image dynamic range extraction and analysis method of claim 5, wherein, In step 4, the Canny-Hough combined algorithm is used in the same time period. On remote sensing images, all linear cracks are extracted and their orientation histograms are calculated. The orientation frequency of each degree is normalized into a probability. Specifically, this involves: using the Canny edge extraction operator to extract the binary edge map of the remote sensing image; then performing Hough line detection on the binary edge map to obtain the detected lines and their corresponding orientation angles; dividing the 0 to 180° interval into 180 orientation boxes, and then calculating which orientation box each line's orientation angle falls into to construct an orientation histogram; based on the orientation histogram, the orientation frequencies are normalized to orientation probabilities with a sum of 1. Anisotropic entropy modulation of cracks for: 。

Images