Method for extracting and analyzing dynamic range of damaged landform feature image of mining area

Abstract

The invention relates to the technical field of image analysis, in particular to a mining area damaged landform feature image dynamic range extraction and analysis method, which comprises the following steps of: 1, performing integral operation on multispectral radiation brightness in a set wavelength range, and deducting two atmospheric scattering terms from each wavelength point to obtain corrected reflectivity; 2, multiplying the Euclidean distance by the curvature amplitude of the near-infrared band to obtain a spectrum collapse-curvature index sensitive to browning / green fading; 3, calculating an elastic strain energy density peak value of the overlying strata-coal rock mass in unit area; and 4, on the remote sensing image in the same time period, extracting all linear cracks and counting the orientation histograms of the linear cracks, and then solving the anisotropic entropy modulation quantity of the cracks by using a Shannon entropy formula so as to construct a comprehensive dynamic range damage index to represent the damage degree of each pixel. According to the method, accurate identification and strength grading of the surface damaged area under mining disturbance of the mining area are realized.

Description

Technical Field

[0001] The present invention belongs to the technical field of image analysis, and in particular relates to a method for extracting and analyzing dynamic range of characteristic images of damaged landforms in mining areas. Background Art

[0002] In recent years, with the continued intensification of mining activities, surface damage has become increasingly common, forming complex and diverse collapse basins, fracture networks, exposed coal dust areas, and overburden browning zones. These landform changes not only threaten operational safety and the geological environment but also directly impact the stability of ecosystems and the long-term sustainable development of mining areas. Therefore, accurate monitoring and dynamic assessment of damaged landforms in mining areas have become a critical technical requirement for geological disaster prevention and control, ecological restoration, and mine supervision. Traditional landform damage monitoring methods rely primarily on ground surveys and point-based displacement monitoring, such as deploying dense monitoring networks using the Global Positioning System (GNSS) and installing inclinometers and crack meters in key areas. Although these methods offer high accuracy, they are limited by spatial sparsity, high deployment costs, limited coverage, and long deployment cycles, making it difficult to fully perceive damage patterns in large, rapidly changing mining environments.

[0003] With the advancement of remote sensing technology, dynamic mining monitoring methods based on optical satellite imagery or radar interferometry (InSAR) have been proposed. For example, optical image change detection can be used to identify bare ground expansion, while InSAR time-series analysis can be used to capture surface subsidence rates. These methods have overcome the limitations of ground-based monitoring, enabling large-scale, periodic, and millimeter-scale observations of surface deformation, becoming the mainstream technology for current mining environmental monitoring. However, existing technologies primarily focus on a single physical dimension: InSAR focuses on deformation, while optical change detection focuses on changes in surface cover. Neither fully exploits the multi-source coupling characteristics of mining damage processes. For example, the surface reflectance characteristics of subsidence areas can significantly change due to vegetation degradation, coal dust accumulation, and exposed brown rock. However, simple NDVI (Newton's Surface Vegetation Index) changes often misinterpret seasonal browning or harvesting as mining damage, resulting in a large number of false alarms. Similarly, high deformation rates do 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 characteristics of damaged landforms in mining areas. The lack of a comprehensive description of the dynamic range, grayscale variability and landform structure organization has become a major shortcoming of existing technologies. Summary of the Invention

[0004] The main purpose of this invention is to provide a method for extracting and analyzing the dynamic range of characteristic images of damaged landforms in mining areas. By fusing the spectral collapse curvature characteristics of multispectral imagery, deformation energy indicators obtained by InSAR or GNSS, and anisotropic entropy modulation of crack directions in remote sensing imagery, a unified comprehensive dynamic range damage index is constructed, enabling the precise identification and intensity classification of surface damaged areas disturbed by mining in mining areas. This method requires no empirical parameters, possesses physical quantity traceability, and is cross-platform compatible. It significantly improves the responsiveness and accuracy of remote sensing imagery for identifying actual damaged areas in complex landform environments, providing highly reliable technical support for mine disaster monitoring, geological risk assessment, and ecological restoration planning.

[0005] In order to solve the above problems, the technical solution of the present invention is achieved as follows:

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

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

[0008] Step 2: Based on the corrected reflectance, extract the central reflectance of the green, red, near-infrared, and short-wave infrared bands; calculate the difference between the central reflectance of the green band and the central reflectance of the red band to obtain a first reflectance difference; calculate the difference between the central reflectance of the near-infrared band and the central reflectance of the short-wave infrared band to obtain a second reflectance difference; then calculate the Euclidean distance between 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 / chlorosis;

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

[0010] Step 4: Extract all linear cracks from the remote sensing images of the same time period and calculate their directional histograms. Normalize the directional frequency of each degree into a probability. Then use the Shannon entropy formula to calculate the anisotropic entropy modulation of the cracks. This is used to 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 ρ at wavelength λ n (λ) is:

[0013]

[0014] Among them, L TOA (λ) is the multispectral radiance at wavelength λ, in Wm -2 sr -1 μm -1 ;L path (λ) is the atmospheric path scattered radiation with wavelength λ, in Wm -2 sr -1 μm -1 ; τ d (λ) is the aerosol optical depth at wavelength λ, dimensionless; L sky (λ) is the sky backscattered radiation of wavelength λ, in Wm -2 sr -1 μm -1 ; sec represents the secant function; θ s is the solar zenith angle of wavelength λ, in rad; θ v is the observation zenith angle of wavelength λ, in rad.

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

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

[0017]

[0018] Among them, ρ n (R) is the central reflectivity of the red band; ρ n (G) is the central reflectivity of the green band; ρ n (NIR) is the central reflectivity of the near-infrared band; ρ n (SW 1.6 ) is the central reflectivity of the shortwave infrared band; Represents the curvature amplitude in the near-infrared band.

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

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

[0021]

[0022] Among them, ρ rock is the average density of the overburden; g is the acceleration of gravity; is the Laplace operator; δ is the pixel resolution of InSAR or GNSS; is the derivative of the instantaneous sinking amount s in the x-axis direction; is the derivative of the instantaneous sinking amount s in the y-axis direction; is the first-order gradient norm of the instantaneous sinking amount s.

[0023] Furthermore, in step 4, the Canny-Hough combination algorithm is used to extract all linear cracks from the remote sensing image of the same time period t and to count their directional histograms, and the directional frequency of each degree is normalized into a probability. Specifically, the method 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 straight lines and corresponding directional angles; dividing the 0 to 180° interval into 180 directional bins, and then counting which directional bin the directional angle of each straight line falls in to construct a directional histogram; based on the directional histogram, the directional frequency is normalized to a directional probability P with a sum of 1. t (θ); the fracture anisotropy entropy modulation CAEM is:

[0024]

[0025] Furthermore, in step 4, the comprehensive dynamic range damage index CDRDI(x, y) of the pixel (x, y) is:

[0026]

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

[0028] The method for extracting and analyzing the dynamic range of characteristic images of damaged landforms in mining areas of the present invention has the following beneficial effects: it can realize high-precision and high-reliability identification of mining-disturbed areas on remote sensing images. By constructing a spectral collapse and curvature coupling index, the spectral morphological changes caused by vegetation greening, coal dust deposition and brown rock exposure are fully captured, and the image's response ability to browning and damaged areas is enhanced; at the same time, the composite index of instantaneous subsidence deformation, spatial gradient and curvature is introduced to elevate traditional deformation monitoring to a dynamic assessment from an energy perspective, solving the problem that existing methods cannot quantify the intensity of landform damage; further adopting linear crack direction statistics and entropy modulation mechanism, the orderliness of image texture is incorporated into the dynamic range adjustment model, realizing automatic control of the geometric structure complexity of grayscale anomaly areas. The final constructed comprehensive index does not rely on any empirical weights or training samples, has physical quantity dimension consistency, 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 loss analysis in mining areas, and provides reliable technical support for geological disaster prevention and control, safe mining and ecological restoration. BRIEF DESCRIPTION OF THE DRAWINGS

[0029] Figure 1 A schematic flow chart of a method for extracting and analyzing dynamic range of characteristic images of damaged landforms in mining areas provided by an embodiment of the present invention. DETAILED DESCRIPTION

[0030] In order to enable those skilled in the art to better understand the solutions of the present invention, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiments of the present invention. Obviously, the embodiments described are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts should fall within the scope of protection of the present invention.

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

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

[0033] Specifically, the process begins with the top-of-the-atmosphere radiance recorded by the sensor. However, the raw radiance is influenced by both tropospheric molecular scattering and aerosol scattering. The former is primarily due to Rayleigh scattering of incident light by atmospheric molecules, resulting in a global background elevation. The latter arises from aerosol particles formed by mining dust, coal dust, and operational emissions, which contribute significantly to sky backscatter. To eliminate these two types of interference, path scattering and dust-related scattering are subtracted separately at each discrete wavelength sampling point. Path scattering is typically obtained through dark pixel estimation or radiative transfer simulation, while dust scattering is derived by multiplying the aerosol optical depth retrieved from an on-site sunphotometer with a simulated sky radiation field. After scattering subtraction, to mitigate anisotropic errors in surface reflectance caused by differences in incident and observation angles, a cosine factor of the solar zenith angle and sensor viewing angle is introduced for each pixel to correct for the dilution of the effective incident energy due to oblique illumination and the oblique path increment.

[0034] Afterwards, the corrected radiance at each wavelength is numerically integrated according to the wavelength step size to obtain the integrated energy reflecting the overall radiation response of the pixel. The total incident energy within the same wavelength range is normalized using the solar spectrum constant, making the radiation comparable between images from different dates, different solar altitudes, and different sensors. Because the slope changes in the mining area collapse basin will change the local angle of incidence, the angle correction term in the formula can promptly compensate for the shadow and bright spot differences caused by microtopography; while the aerosol scattering subtraction significantly reduces the gray fog artifacts caused by blasting dust and transportation dust, making the exposed coal seams, browned bare rock, and dust deposition areas present more accurate reflectance characteristics after correction. Through this comprehensive integration and double scattering elimination process, the final output spectral correction reflectance not only retains the true physical reflectance differences of the mining area landforms, but also minimizes atmospheric and geometric noise, establishing a unified, stable, and cross-image consistent radiation benchmark for subsequent spectral collapse detection, strain energy field fusion, and fracture direction entropy modulation. The entire processing chain can be automatically executed on large-scale, high-resolution, multi-temporal remote sensing images of mining areas, ensuring the comparability of the grayscale dynamic range response of damaged landforms under different observation conditions, thereby improving the subsequent dynamic range damage index's recognition accuracy and temporal consistency for complex damage such as mining collapse, coal dust browning, and vegetation degradation.

[0035] Step 2: Based on the corrected reflectance, extract the central reflectance of the green, red, near-infrared, and short-wave infrared bands; calculate the difference between the central reflectance of the green band and the central reflectance of the red band to obtain a first reflectance difference; calculate the difference between the central reflectance of the near-infrared band and the central reflectance of the short-wave infrared band to obtain a second reflectance difference; then calculate the Euclidean distance between 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 / chlorosis;

[0036] After radiometric correction, all pixels are converted to a surface reflectance raster that is comparable across all bands. The second key task is to analyze the spectral collapse caused by mining damage from these corrected reflectances and quantify it using a single, comparable index. First, four physically distinct central bands must be precisely identified within the corrected imagery: the green band records the high chlorophyll reflectance peak of healthy vegetation; the red band lies within the strong chlorophyll absorption trough; the near-infrared band corresponds to the high specular reflectance region of vegetation cellular structure; and the 16-micron shortwave infrared band is most sensitive to changes in soil and dust moisture content. After mining operations, vegetation typically withers or is covered by dust, causing the green peak to weaken dramatically while the red trough rises with browning and oxidation, narrowing the difference between the two. Simultaneously, the high near-infrared reflectance decreases sharply when vegetation turns green or coal dust accumulates, while the shortwave infrared exhibits different changes due to water evaporation or exposure of bare soil, resulting in a step-like convergence in the difference between the near-infrared and shortwave infrared values.

[0037] To weave these two difference values into a grayscale scale highly sensitive to greening and browning in mining areas, the first and second differences must be treated as orthogonal components in Cartesian coordinates. Their Euclidean distance is then used to map the dual differences to a single amplitude. This amplitude can be interpreted as the length of the "collapse" vector in spectral space; larger values indicate a more severe deviation of the pixel's spectral morphology from native vegetation or undisturbed surface. Relying solely on the difference length is insufficient to capture the sharp bends in the browning spectrum in mining areas. Therefore, it is necessary to further exploit the second-order morphological characteristics of the local spectral curve within the near-infrared band to enhance the response to subtle greening. In engineering implementation, a multispectral reflectance sequence can be extracted by wavelength. The second-order curvature amplitude can be calculated using five-point central differences or the Savitzky–Golay smoothed derivative method for several equally spaced bands before and after the near-infrared center point. Larger amplitudes indicate sharper bends in the near-infrared curve, typically corresponding to a precipitous decline in vegetation health or rapid coal dust accumulation. Then, the curvature amplitude is directly multiplied and fused with the aforementioned Euclidean distance, which is equivalent to considering both the overall spectral collapse scale and the local bending intensity in the grayscale space. This not only amplifies the response of the extreme browning area, but also suppresses the false collapse signal caused by trace noise.

[0038] To ensure robustness and stability during batch processing of large scenes, the actual central wavelengths of the four central bands must be carefully aligned during implementation to avoid deviations in the central values due to slight sensor drift. Furthermore, the images must be subjected to noise removal, cloud masking, and signal-to-noise weighting adjustments before performing interpolation and curvature calculations. If different satellite platforms are used for mining area imagery, the discrete spectra of each platform can be converted to a unified central 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, in the form of grayscale intensity, directly contributes to the subsequent comprehensive calculations of the elastic strain energy field and fracture directional entropy. Because this index fully integrates the spectral morphological changes caused by vegetation degreening, browning, coal dust deposition and water content mutation, which are the most typical surface responses to mining disturbances in mining areas, it can provide extremely high discrimination in dynamic range mapping, so that key areas such as the collapse center, the foot of the spoil dump and the dust channel are significantly brightened or darkened 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 the first-order gradient norm of the instantaneous subsidence, the coupling of the subsidence amplitude and deformation curvature is obtained. Combined with the measured density of the overburden, the peak value of the elastic strain energy density per unit area of the overburden-coal rock mass is calculated;

[0040] After completing collapse detection using spectral information, it is necessary to link the three-dimensional deformation field of the mining area surface caused by mining with the radiation dynamic range. The key to step 3 is to convert the instantaneous subsidence into an energy scale that can be compared with the spectral index to reveal the coupled characteristics of the collapsed surface in both grayscale and mechanical space. This is achieved by starting with time-series InSAR or high-frequency GNSS to acquire millimeter-level vertical displacements. First, through multi-track registration, track error correction, and reference area constraints, the displacement grid is unified to the absolute vertical frame. Then, the instantaneous subsidence is firstly subjected to first-order central difference in the pixel grid to obtain the spatial gradient field along the east-west and north-south directions. The gradient field reflects the slope direction and steepness of the inclination at each location in the subsidence depression. The larger the gradient, the shorter the distance at which the same depth of subsidence is completed, indicating that the supporting coal pillars or overburden slab beams below the goaf have undergone strong differential deformation. A Laplace convolution operation is then performed on the subsidence grid to obtain a second-order spatial curvature distribution. The high-value curvature area corresponds to the location with the smallest bending radius, which is often the location where the overburden is about to undergo shear cracking or has already developed tension cracks. Multiplying these two mechanical quantities can integrate the displacement amplitude information and the deformation bending information into the same pixel indicator, thereby simultaneously reflecting the subsidence depth and spatial abrupt change at the single-value level.

[0041] To give this metric a clear physical energy meaning, it is necessary to incorporate the volume density of the overburden in the mining area. This density is typically obtained through core testing or borehole acoustic logging. Multiplying this density by the regional constant gravitational acceleration elevates the displacement-curvature combination to an approximate scale of elastic potential energy per unit area. To reduce discrete differential noise during calculation, the gradient and curvature are smoothed using a Gaussian weighted kernel before and after convolution. The kernel scale is set based on the spacing between major faults or the working face length revealed by seismic exploration. This preserves the large-scale depression framework caused by underground mining while suppressing high-frequency interference. In thick coal multi-stage mining or multi-layered mining scenarios, the same pixel may have overlapping displacements at different depths. The multi-layer displacements are weighted averaged according to the overburden depth before being used in the gradient and curvature calculations to prevent shallow micro-motion from masking the deeper main mining signal. 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, ensuring a one-to-one correspondence at the pixel level when constructing the integrated dynamic range damage index.

[0042] High-value areas are often concentrated in the center of the working face, at the intersection of the upper and lower walls of the damage zone, and at the edge of the spoil dump. These areas not only exhibit strong sudden changes in brightness and darkness in the radiation image but also accumulate the maximum potential energy in the mechanical field. Therefore, this energy peak indicator provides a mechanical constraint for damage classification. To enhance dynamic monitoring capabilities, incremental differences can be made in the displacement field over different observation periods to obtain the instantaneous elastic energy growth rate over a short period of time. A time sliding window filter is then used to generate a rate sequence, which is used to identify locations where secondary collapse or accelerated slope instability is imminent. The entire process can be automatically streamlined in the cloud, achieving energy translation for large-scale, long-term, and uninterrupted deformation in the mining area.

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

[0044] In the implementation, noise suppression and tonal equalization are first performed on high-resolution remote sensing images acquired at the same timestamp as the deformation field, ensuring sufficient gradient contrast for faint cracks obscured by coal dust. A Canny edge operator is then used to construct a gradient magnitude and direction grid for all pixels. A double-threshold concatenation is then used to filter out pseudo-edges caused by gravel and vegetation gaps, retaining only thin line pixels whose length exceeds the threshold and whose gradient direction is consistent. After obtaining a binary edge map, a probabilistic Hough transform is applied to accumulate linear features in polar coordinate space, outputting the pixel start and end coordinates of each crack. These coordinates are then mapped back to geographic space to determine the azimuth of each crack. The range from 0 to 180 degrees is then divided into azimuth bins with a resolution of one degree, and all cracks are counted based on the azimuth bins. Because different crack lengths contribute differently to the stretching of the visual dynamic range, the total length of the crack pixels is used as the frequency criterion, rather than the number of cracks, to ensure that long-scale tensile cracks receive a higher weight in the probability distribution. After adding up the lengths in each direction and dividing it by the total length of the crack in the entire field, we get a set of directional probabilities. Then, using the Shannon entropy calculation method in information theory, we take the logarithm of each directional probability and multiply it by itself. Finally, we integrate and sum it within the range of zero to one hundred and eighty degrees, and then normalize it to the maximum theoretical entropy.

[0045] Higher values of the fracture anisotropy modulation indicate a more uniform distribution of crack lengths in different directions and a more complex texture. Conversely, values close to the lower limit indicate concentrated and ordered crack orientations, facilitating overall structural enhancement of the image's grayscale. To allow this texture complexity to exert a suppressive effect on the dynamic range, its square root is placed in the denominator of the composite index. This amplifies the denominator and reduces the dynamic range in high-entropy regions, while maintaining or enhancing the dynamic response in low-entropy regions. To prevent directional statistics from being distorted by pseudo-crack linear objects such as cloud shadows and bare rock, spectral thresholding is used before crack extraction to remove highly reflective bare rock, local slope aspect analysis is used to remove cliff shadows, and vegetation elevation texture masks are used to exclude gaps in tree canopies. During repeated blasting at the working face, images collected over multiple time periods can be used to perform a sliding average of the crack orientation probabilities. This allows the azimuthal noise of temporary shallow cracks to be attenuated by temporal filtering, while persistent, stable tensile cracks contribute significantly to the directional entropy, thus reflecting long-term dynamic risk in the index.

[0046] Furthermore, to accommodate differences in spatial resolution across satellite or UAV platforms, the crack pixel lengths must be normalized to their actual ground dimensions to ensure comparability of directional entropy output from images of varying scales. Ultimately, the fracture anisotropic entropy modulation is superimposed at the pixel level with the spectral collapse curvature index and the peak elastic strain energy density. After denominator modulation and smoothing threshold processing, a comprehensive dynamic range damage index grid is generated. Its values directly correspond to high-risk areas such as mining subsidence centers, crack-dense zones, and unstable edges of spoil dumps, achieving cross-domain coupling from local texture to global dynamics, allowing the entire extraction and analysis method to form a complete closed loop with spectral, deformation, and geometric constraints.

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

[0048] The 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 the visible (blue-green-red), near-infrared, and two representative short-wave infrared bands. It encompasses the chlorophyll absorption-reflectance 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 fully captures all key spectral responses after mining activities. At the same time, it excludes the long-wave infrared band above 2.2 μm, which has a low signal-to-noise ratio and significant interference from thermal radiation, to avoid grayscale distortion introduced by irrelevant energy. In practical implementation, all input images should be scatter-subtracted, angle-corrected, and integrated at each discrete central wavelength within this band window. Missing channels can be resampled by interpolating adjacent bands or by convolution using the sensor response function. This ensures the integrity and comparability of the integrated energy within 0.45–2.20 μm, establishing a unified benchmark for subsequent spectral collapse and curvature detection.

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

[0050]

[0051] Among them, L TOA (λ) is the multispectral radiance at wavelength λ, in Wm -2 sr -1 μm -1 ;L path (λ) is the atmospheric path scattered radiation with wavelength λ, in Wm -2 sr -1 μm -1 ; τ d (λ) is the aerosol optical depth at wavelength λ, dimensionless; L sky (λ) is the sky backscattered radiation of wavelength λ, in Wm -2 sr-1 μm -1 ; sec represents the secant function; θ s is the solar zenith angle of wavelength λ, in rad; θ v is the observation zenith angle of wavelength λ, in rad.

[0052] ρ n (λ) Unifies multispectral observations into normalized reflectance that can directly represent surface energy exchange. The upper and lower limits of the formula are limited to 0.45μm to 2.20μm. Its physical motivation is to include all the key channels of visible-near infrared-shortwave infrared that are most sensitive to vegetation chlorosis, coal dust browning, bare rock oxidation, and water content changes, while screening out the irrelevant surface radiation and atmospheric radiation mixing terms in the longer-wavelength thermal infrared. The first appearance of L in the numerator TOA (λ) is the top atmospheric radiation brightness received by the sensor. It not only contains the actual reflection information of the target pixel, but also superimposes the Rayleigh scattering and Mie scattering on the path from the sun to the surface and then to the sensor. If it is directly used for grayscale analysis of mining areas, the brightness of non-surface sources will be mistakenly regarded as the brightness and darkness changes caused by collapse or dust, so it must be subtracted in two levels. The first level is to subtract L path (λ), removes uniform atmospheric path scattering; the second stage uses aerosol optical depth τ d (λ) multiplied by the sky backscattered radiation L sky (λ) specifically removes the wavelength-dependent brightness increase caused by aerosols such as mining dust and blasting dust. This step is particularly important in mining scenes because dust in goafs often causes the near-infrared low-reflection area to be scattered and raised by the sky. If not removed, it will mask the vegetation greening signal. s (λ)+secθ v The (λ) term belongs to the first-order anisotropy correction. In mining subsidence basins or steep slopes of spoil dumps, the local incident angle and observation angle both change, causing the same material to exhibit completely different brightness at different slope directions. Taking the secant of the solar zenith angle and the sensor viewing angle and adding them together is equivalent to geometrically compressing the path increments caused by oblique illumination and oblique viewing angle back to the radiation benchmark of the vertical viewing angle, achieving radiation comparability between pixels of different slope directions while retaining the dynamic contribution of true shadows and highlight areas in grayscale.

[0053] The integral operation accumulates the above-mentioned wavelength point brightness after scattering deduction and angle correction within 0.45–2.20 μm to obtain the total energy actually returned by the pixel; in order to eliminate the influence of the difference in imaging date, sun-earth distance and solar altitude on the energy benchmark, the denominator is the solar constant irradiance E0(λ)cosθ in the same band window. s(λ) is integrated, which acts as a normalization scale: pixels with the same reflectance are mapped to a consistent 0–1 range regardless of whether they are in winter images with low solar altitude or summer images with high solar altitude, ensuring that the subsequent cross-phase spectral collapse index can be directly compared. 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 not only eliminates external noise from dust and path scattering, but also intrinsically compensates for local terrain geometry, enabling the extraction of spectral differences representing true surface properties even in the multi-source disturbance environment of the mining area. Because this normalized reflectance maintains energy conservation at the pixel scale, it serves as a direct input for the subsequent calculation 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 phase is entirely due to the physical state of the surface, rather than atmospheric or geometric artifacts. This ensures that the final integrated dynamic range damage index has both radiometric rigor.

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

[0055] In the dynamic range extraction process for images of damaged landforms in mining areas, the central wavelengths of 0.56μm, 0.66μm, 0.86μm, and 1.6μm were selected as the green, red, near-infrared, and shortwave infrared wavelengths to maximize the capture of four complementary sensitive regions affected 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 vegetation loses green due to road excavation, spoil dumping, or dust cover in mining areas, this peak is the first to attenuate. Therefore, it can sensitively record the brightness collapse caused by the decline in vegetation biochemical function. The adjacent red band at 0.66μm lies in the main chlorophyll absorption trough, where normal vegetation should exhibit low reflectance. However, coal dust adhesion or exposed loess caused by mining reduces the absorption depth of the red band, increasing its brightness. As a result, the difference between the weakening green peak and the rising red trough suddenly narrows, forming a typical spectral "collapse." Moving further to the long wavelength of 0.86 μm, the near-infrared band corresponds to the specular scattering region within vegetation cell walls. Healthy canopies have extremely high reflectance in this region, but once vegetation is removed or covered by coal dust, the near-infrared reflectance drops dramatically. Furthermore, exposed overburden or water accumulation in collapsed basins creates additional curvature in the near-infrared. Therefore, using this band as a reference point for curvature calculations significantly amplifies the sharp inflection in the spectral curve, making it highly sensitive to both browning and chlorosis. Finally, the 1.6 μm short-wave infrared central band is introduced because it is particularly sensitive to surface moisture content, ore oxidation, and dust particle size distribution. Weathered rock in spoil dumps and goaf floors exhibit strong absorption differences in this band when exposed to alternating wet and dry conditions or oxidative browning. The difference between this band and the high near-infrared reflectance becomes a key indicator of exposed bare soil and rock in mining areas. Furthermore, coal dust particles absorb less in the short-wave infrared than in the near-infrared, further reducing the contrast between the two bands in areas with thick coal dust accumulation and exacerbating spectral collapse. By forming two sets of differences from these four central wavelengths and calculating the Euclidean distance, the combined effects of vegetation decline, coal dust cover, oxidative browning, and water content changes can be combined into a single amplitude in two-dimensional spectral space. Multiplying this amplitude with the near-infrared curvature couples the overall brightness decay with the local curvature convexity, generating a spectral collapse-curvature index that is most discriminative of post-mining surface browning and greening. This band configuration avoids redundant calculations caused by using too many channels while ensuring that all key mining damage signals are fully represented within the selected wavelength window. This provides the most comprehensive physically comprehensive input with minimal spectral redundancy for subsequent integration with deformation energy and crack entropy modulation.

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

[0057]

[0058] Among them, ρ n (R) is the central reflectivity of the red band; ρn (G) is the central reflectance of the green band;

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

[0060] In the dynamic range extraction and analysis method of characteristic images of damaged landforms in mining areas, the spectral collapse-curvature index (SCCI) is designed to simultaneously measure the overall shrinkage effect and local inflection effect of vegetation loss, coal dust coverage, bare rock oxidation, and changes in soil dump moisture on the surface spectral morphology. Its core formula consists of the multiplication of two parts. 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. Once mining causes vegetation to become green or covered by coal dust, the peak value decays rapidly; ρ n (R) is located at the 0.66μm chlorophyll absorption trough. When coal dust or bare soil covers the leaf surface, the red trough reflectance will rise; the absolute contraction of the difference between the two directly quantifies the intensity of vegetation damage. At the same time, ρ n (NIR) records high reflectance of cell structures at 0.86 μm, which decreases significantly when vegetation dies or browns, while ρ n (SW 1.6 ) is highly sensitive to moisture content and oxidation at 1.6μm. Changes in the wetness and dryness of the spoil dump or the browning of the rocks can cause independent fluctuations in the reflectance in this band. The difference between the two is orthogonal to the previous difference, compensating for moisture and oxidation information missed by a single difference. When mining disturbances cause both the green-red and near-infrared-shortwave infrared differences to narrow simultaneously, the Euclidean distance decreases significantly, equivalent to forming a "spectral collapse" vector in four-dimensional grayscale 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 across multiple spectral modalities.

[0061] Part 2 The second-order derivative amplitude of the central near-infrared band quantitatively describes the curvature of the spectral curve at the interface between high vegetation reflectance and water absorption. When vegetation is healthy, the curve shows a smooth, gentle slope between the near-infrared and short-wave infrared bands. However, when browning or coal dust accumulation occurs rapidly, the curve drops sharply from high to low reflectance, and the absolute value of the second-order derivative increases dramatically. This jump often precedes the green-red difference, making it detectable in the early stages of greening. Furthermore, when soil moisture in a dumpsite evaporates or absorbs moisture rapidly over a short period of time, local peaks and valleys can be introduced near this band, further amplifying the curvature. Multiplying the curvature amplitude by the Euclidean distance ensures that the SCCI retains a large Euclidean distance in healthy vegetation areas, while remaining neutral overall due to the near-zero curvature. In areas such as collapse centers, thick coal dust accumulation, and dumpsite slopes, the Euclidean distance converges while the curvature increases, causing the product output to significantly decrease or increase, resulting in a two-way stretch of the dynamic range. This dual-factor coupling not only avoids the over-amplification of shadows or sensor noise by a single difference, but also suppresses false peaks in curvature under noisy conditions. More importantly, SCCI uses two-dimensional differences and second-order spectral information, which have complementary physical meanings, to simultaneously compress the radiation intensity jumps and curve bends caused by mining activities into single-value indicators at the imaging element scale. This allows the dynamic range of the spectral dimension to be fully mapped into grayscale space, providing highly differentiated and numerically stable input for subsequent fusion with deformation energy fields and crack direction entropy. Within the entire mining area, high SCCI value areas generally correspond to marginal areas where vegetation still exists but has been affected by the stress of the overburden, the median value corresponds to the background undisturbed surface, and the low value falls precisely in the core areas with the most severe coal dust coverage or the most significant collapse and water accumulation. This gradient distribution forms a spatial coupling with the InSAR displacement gradient band and the crack-intensive zone, enabling the comprehensive dynamic range damage index to achieve pixel-level closed-loop response in the multi-source information space.

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

[0063] The instantaneous subsidence s is the bridge connecting the radiation grayscale and the mechanical field. Its calculation relies on InSAR points or high-frequency GNSS points in the same period. For InSAR, the core mechanism is to use the phase difference captured by the satellite's repeated orbital imaging. That is, the interference fringes obtained by interfering two radar complex images record the subtle changes in the round-trip path from the satellite to the surface. The number of fringes can be converted into displacement along the line of sight through phase unwrapping. The vertical deformation caused by mining in the mining area is dominant. Therefore, the line of sight displacement can be restored to approximate vertical settlement using simple geometric projection under the premise of known radar incident angle. In order to obtain millimeter-level accuracy and eliminate atmospheric delay and orbit control errors, permanent scatterer or distributed scatterer time-series interferometry methods are usually used. After applying stable reference area constraints to multiple images throughout the year or multiple months, the instantaneous line of sight displacement of each period is jointly solved, and then differentiated to the specified observation date to form a settlement field grid for the same period. If a GNSS continuous station is used, 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 overburden destabilizes, the vertical coordinates of the station will shift downward on the order of centimeters to decimeters. The instantaneous vertical displacement is obtained by subtracting the station coordinates during the target period from those during the reference period. GNSS provides an absolute displacement reference, which can be used to correct for potential vertical offsets in InSAR. InSAR, on the other hand, its dense spatial sampling compensates for the sparseness of GNSS points. By integrating the two, the true GNSS values can be embedded into the radar displacement plane through least-squares adjustment, resulting in a high-precision, high-resolution instantaneous subsidence field. To ensure close correspondence with multispectral imagery, radar image tracks from the same day or week as the optical imagery should be selected, or coordinate solutions within a GNSS continuous time series with a time difference of less than a few hours from the image acquisition time should be extracted. This ensures temporal synchronization between the radiance and mechanical subsidence, ensuring that there is no timing mismatch when the pixel-level elastic strain energy and the spectral collapse index are multiplied in parallel in subsequent calculations. Subsidence surfaces in mining areas typically exhibit an elliptical-parabolic distribution. The s-grid calculated using InSAR or GNSS not only captures the depth of the depression but also preserves the subtle transitions in the boundary gradient. This allows the Laplace operator to be used to extract the bending stress concentration zone. This, combined with the first-order gradient modulus, forms the coupling between the subsidence amplitude and curvature, providing reliable source data for elastic strain energy density. The entire deformation solution chain, based on the physical delay of the radar phase or the integer deviation of the carrier phase, allows for traceable and quantifiable errors, making the geomorphological dynamics information derived from the integrated dynamic range damage index engineering verifiable.

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

[0065]

[0066] Among them, ρ rock is the average density of the overburden; g is the acceleration of gravity; is the Laplace operator; δ is the pixel resolution of InSAR or GNSS; is the derivative of the instantaneous sinking amount s in the x-axis direction; is the derivative of the instantaneous sinking amount s in the y-axis direction; is the first-order gradient norm of the instantaneous sinking amount s.

[0067] The peak value of elastic strain energy density achieves a quantitative description of the potential elastic potential energy accumulation degree of overburden-coal rock mass by coupling the amplitude term, curvature term and spatial gradient term of the instantaneous subsidence field at the same pixel scale. rock g / 2 is derived from the classical expression for the potential energy density of a deformed solid unit, where ρ rock Derived from core density measurements or borehole acoustic inversion, it reflects the average unit volume weight of the overburden. g is the gravity constant, used to convert pure geometric displacement into gravitational potential energy, giving the indicator an energetic dimension. The denominator's double coefficient derives from the symmetry between elastic strain energy density and displacement Coulomb potential energy under the small strain assumption, ensuring dimensional closure of the amplitude. The two terms within the absolute value reflect the coupling between amplitude and curvature: s itself records the vertical settlement depth; deeper depth indicates greater space for overburden detachment and more intense stress redistribution. Multiplying the displacement Laplace term by the square of the pixel resolution actually converts the second-order spatial curvature on the pixel grid into an additional settlement of the same dimension as the displacement amplitude. The higher the absolute value of the curvature, the smaller the bending radius of the settlement surface and the more concentrated the local tensile stress or compressive stress. Multiplying by δ 2 This avoids the artificially high curvature caused by grid refinement. 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 in the spoil dump accumulate as positive energy in the index, which is consistent with the fact that both tension and compression zones can store elastic energy in dynamics.

[0068] The gradient modulus immediately following Reflects the spatial rate of change of surface subsidence in the horizontal plane, that is, the steepness of the slope of the depression boundary. The larger the slope, the more significant the difference in displacement between adjacent pixels, and the more likely the shear or tensile stress is to be concentrated in the slope break zone. Therefore, it acts as a multiplicative amplifier in the indicator, making the collapse edge obtain a higher energy estimate than the central stress-smooth area. The entire expression projects the amplitude, curvature, and slope of the three-dimensional settlement shape into a scalar at the pixel level, and its unit is Jm -2, can be incorporated into the numerator of the comprehensive dynamic range damage index along with the spectral collapse-curvature index of the same resolution, injecting geomorphological mechanical intensity information into the brightness signal of the grayscale dynamic range. The formula structure also takes into account resolution portability: when using InSAR data from different orbits or GNSS data with different station densities, simply updating δ while maintaining the adaptability of the differential kernel scale can eliminate the spurious gain in the curvature term caused by grid refinement at the sub-pixel level, ensuring consistent energy estimates across data sources. In actual mining areas, high-value areas output by this index typically correspond to the advance arch curve in the direction of working face advancement, the shear ring at the center of the collapse, and the stress lag zone at the foot of the spoil dump. These areas often exhibit both spectral collapse and dense crack textures in multispectral imagery. This allows for a high degree of synergy between the mechanical and radiation fields when the comprehensive index is overlaid, improving the accuracy of identifying high-risk damage units and the dynamic range stretching effect.

[0069] Furthermore, in step 4, the Canny-Hough combination algorithm is used to extract all linear cracks from the remote sensing image of the same time period t and to count their directional histograms, and the directional frequency of each degree is normalized into a probability. Specifically, the method 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 straight lines and corresponding directional angles; dividing the 0 to 180° interval into 180 directional bins, and then counting which directional bin the directional angle of each straight line falls in to construct a directional histogram; based on the directional histogram, the directional frequency is normalized to a directional probability P with a sum of 1. t (θ); the fracture anisotropy entropy modulation CAEM is:

[0070]

[0071] On the one hand, goaf roof fractures and dump instability can produce tensile or shear cracks with predetermined orientations on the surface. These linear features form brightness sharpening bands in high-resolution reflection images. Once their orientation is concentrated, they are equivalent to superimposing high-amplitude, low-entropy band 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 the overburden enters late-stage waterlogged collapse or the dump slope undergoes multiple rounds of sliding, cracks often expand disorderly in radial, fan-shaped, or grid-like patterns. At this time, the linear texture energy is diluted by the directional dimension, and the isotropic background still dominates. Even if the image enhancement algorithm increases the contrast, it will only produce a mean shift without regional bright-dark jumps. In order to ensure that the integrated dynamic range damage index can adaptively adjust between these two extremes, the entropy of the directional probability distribution must be used to constrain the final stretch of the grayscale dynamic range.

[0072] First, the original image is spectrally equalized and median filtered to eliminate local high-frequency noise. The Canny edge operator is then applied to extract significant grayscale gradient bands. Canny's dual-threshold connection mechanism ensures the continuity of crack pixels while eliminating speckle artifacts caused by uneven coal dust particle size. The resulting binary edge map provides reliable pixel candidates for subsequent line identification. A probabilistic Hough transform is then applied in polar coordinate space. An accumulator maps the parameters of the possible lines for each pixel to (ρ, θ) peaks. After thresholding, pixel coordinates are back-traced to isolate all potential crack segments and determine their azimuth angles θ. Since crack directions do not distinguish between positive and negative directions, simply modulo-scaling the azimuth angle to the range of 0° to 180° covers all directions. To ensure that the contribution of long cracks to directional statistics is proportional to their effect on the visual dynamic range, the directional frequency is accumulated using pixel length rather than the number of cracks. The 0° to 180° range is divided into 180 directional bins of equal width, and the length frequency f is calculated for all cracks within each bin. t (θ), and then normalized to the direction probability The entropy measurement part introduces the direction angle square weight θ 2 This is because in imaging geometry, vertical or large-angle cracks have a more significant amplification effect on local shadow and illumination differences, and visually have a greater impact on the grayscale dynamic range; if the crack direction is concentrated in the high-angle area, the geometric texture can be enhanced accordingly. Conversely, if the distribution is uniform, the anisotropy is diluted and it is difficult to further increase the contrast. Based on information theory, the higher the degree of randomness, the greater the entropy value. The classic form of directional entropy is -∑PlnP, which is obtained by multiplying by θ 2 Introducing the second-order moment of angle, the final entropy is written as Considering that the maximum entropy occurs in a completely uniform distribution, the theoretical value is ln180°, using Normalization is performed to compress the entropy value between 0 and 1, and a constant 1 is added to make the index always greater than 1 to obtain the anisotropic entropy modulation of the fracture.

[0073] This index has three physical meanings: First, if the crack direction is single, such as a strip-shaped tensile crack dominated by the main fault, most of the probability is concentrated in a few direction boxes, the entropy term approaches zero, the CAEM is close to 1, and the comprehensive dynamic range damage index is 1 in the denominator. The grayscale dynamic range is suppressed to the minimum, and the ordered crack texture is further amplified. Secondly, if the crack direction is highly discrete, radial or checkerboard, then P t (θ) is close to uniform, the entropy term approaches ln180°, and CAEM approaches 2. The increase in the denominator significantly weakens the grayscale dynamic response, avoiding the noise peak caused by texture clutter from being mistakenly enhanced. Third, through θ 2After weighting, even under the same entropy value, high-angle cracks have a greater weight, so vertical cracks generated by steep slopes at the edge of the collapse will have higher geometric enhanceability than low-angle cracks generated in the flat area of the spoil dump. It is worth noting that in real mining scenes, cloud shadows, bare rock joints or road edges are often mistakenly identified as cracks. For this reason, the multispectral NDVI threshold and DSM slope are fused before Canny to screen out plant texture and steep cliff elevation mutations, and then the minimum length threshold is set in the Hough result to eliminate short-distance pseudo-straight lines to ensure that the directional probability truly represents the overburden crack network. For multi-time series images, in order to eliminate small crack noise induced by periodic road maintenance or temporary drainage ditches, P t A time sliding average of (θ) is then applied to the entropy formula to generate a robust directional chaos curve, thereby reflecting a gradual collapse evolution rather than instantaneous noise in the integrated dynamic range loss index. The entire process is based entirely on native image pixels and mathematical statistics, requiring no empirical parameters or training samples. It is applicable to satellite and UAV data of varying spatial resolutions. Cross-platform consistency is ensured by simply updating the minimum line length of the Hough detector and the θ sampling accuracy based on resolution. Verified by a large number of mining case studies, CAEM is highly correlated with the anisotropic ductility of the collapse basin profile and the control effect of the main fault: in the strong tensile stress working face of the open-pit mine, early cracks develop rapidly in one main direction, and CAEM remains at around 1.05, with minimal inhibition on grayscale enhancement. The comprehensive index forms a high brightness wall at the collapse boundary; as the working face expands, boundary cracks begin to derive in the vertical direction, CAEM rises to 1.4, and the brightness peak drops relatively, reflecting the weakening of the tension-shear transformation on texture enhanceability; in the wet collapse area of the spoil dump, cracks spread disorderly in a honeycomb shape, CAEM rises to 1.8, and the grayscale dynamic range is significantly compressed, avoiding erroneous indications of high damage. It can be seen that the anisotropic entropy modulation of the fracture not only provides an automatic suppression weight of the geometric texture dimension for the comprehensive dynamic range damage index, but also realizes the 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, so that the grayscale dynamic expression and the real landform evolution maintain a high degree of consistency in time and space, thereby promoting the mining area surface disaster monitoring from a single radiation source to a new engineering paradigm of multi-source coupling and high dynamic range quantification.

[0074] Furthermore, in step 4, the comprehensive dynamic range damage index CDRDI(x, y) of the pixel (x, y) is:

[0075]

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

[0077] The ln[1+SCSED] in the numerator first applies a natural logarithmic compression to the peak energy density obtained by coupling the 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 will cause the energy in the deep collapse zone to dominate the composite energy, obscuring the moderately damaged area. Second, after taking the logarithm of the exponential change, the high-energy and low-energy areas are compressed back to comparable ranges, achieving a "logarithmic stretch" of the dynamic range rather than a linear collapse. At the same time, [SCCI] 2 By squaring the product of the spectral collapse vector length and the near-infrared curvature amplitude, pixels with slightly contracted Euclidean distances or slightly increased curvatures are geometrically amplified in the index, highlighting subtle changes in browning, chlorosis, and coal dust coverage. When the two are multiplied together, the numerator encompasses both the severity of the collapse at the mechanical potential energy level and the dramatic jumps in the spectral morphology at the radiation level, resulting in truly high-risk pixels exhibiting a brightness potential far higher than the background in the product domain. The denominator employs a two-term superimposed suppression weighting structure: The square root of the crack entropy modulation is taken because a high dimension of the original entropy value will over-suppress the spectral and energy information, while the square root can retain the entropy's ability to distinguish the texture chaos while avoiding denominator explosion. When the crack directions are highly dispersed, the CAEM approaches its upper limit, and the increase in the denominator weakens the collapse-curvature product, reflecting that it is difficult to enhance the brightness of the disordered texture area. When the crack direction is single, the CAEM approaches 1, and the denominator has almost no suppression on the dynamic range, so that the high grayscale potential of the ordered crack zone is fully released.

[0078] Item 2 The physical meaning of introducing the z-score square of NDVI is that when the vegetation index of a certain pixel is significantly higher than the mean value of the same period, it is very likely that it is a natural forest or farmland with high reflectivity, not a brightness anomaly caused by mining damage; adding the square of the deviation degree to the denominator will produce additional suppression on these high NDVI pixels, preventing natural highlights from being mistaken for high-risk damage areas against the background of hilly woodlands or river valley oases; on the contrary, if the pixel NDVI is far below the mean, on the one hand, it means that the vegetation is degreening or covered by coal dust, and on the other hand, the z-score is negative, and its square is also positive, but the deviation direction does not affect the size of the suppression weight, thereby maintaining the dynamic range response to the degraded vegetation area. The standard deviation of this item σ NDVI，tUsing the overall NDVI distribution over time eliminates seasonal differences in greenness, allowing NDVI deviations between winter deserts and summer oases to be assessed at the same scale. This also ensures that vegetation degradation within mining areas caused by mining activities is consistently mapped over time. Overall, the algebraic construction of the CDRDI follows three principles: First, the numerator uses multiplicative coupling, exponentially amplifying high spectral and mechanical values when they occur simultaneously, strengthening truly dangerous pixels. Second, the denominator uses a dual weighting scheme using interpretable entropy roots and z-score squares, allowing areas with chaotic textures or naturally high vegetation to automatically reduce their dynamic contributions. Third, all inputs are calculated from physically measurable indicators from previous steps, without any parameter adjustments. This ensures that the index remains comparable in absolute terms across satellite platforms, seasons, and even mining areas. Field measurements have verified that in dark inclined shaft collapse areas, while SCSED values are high in large displacement zones, crack directions are concentrated, CAEM is as low as 1.1, and combined with sparse vegetation and low NDVI, the CDRDI is as high as 18, significantly different from the background 2–4 areas. In open-pit dumps, despite high curvature curves, crack directions are extremely chaotic, with entropy values approaching 1.9. Combined with the high NDVI of local herbaceous cover in the dumps, the CDRDI is ultimately controlled at 6–8, avoiding misjudgment of large areas of medium- and low-risk areas. This demonstrates that the CDRDI, which incorporates a four-fold transformation consisting of logarithmic compression, square amplification, entropy root suppression, and NDVI normalization, not only mathematically ensures the numerical stability of the grayscale dynamic range but also physically fully incorporates the four major mechanisms of mining area spectral collapse, subsidence energy, crack direction, and vegetation background. This provides a highly robust quantitative tool for disaster loss monitoring in large-scale mining areas, with direct threshold classification and cross-source consistency.

[0079] We acquired local transiting Sentinel-2L1C multispectral images, a Sentinel-1 coherent permanent scatterer InSAR deformation sequence from the same orbit, and a 0.15m resolution visible light aerial survey image acquired by a drone on the same day. To ensure cross-source consistency, all data were resampled to a δ = 10m grid and projected to WGS84 / UTM50N. The following demonstrates the calculation process of all indicators using the pixel corresponding to the grid center coordinate (x0, y0). First, a band-discrete integration of the radiometrically calibrated Sentinel-2 top-of-atmosphere radiance was performed over the 0.45μm to 2.20μm band window, accumulating the term [L TOA (λ)-L path (λ)-τ d (λ)L sky (λ)][secθ s (λ)+secθ v (λ)] was summed over a 10 nm step to obtain 1.92×10 2 Wm -2 sr -1, and using the ASTM E490 solar spectrum constants and combined with the metadata solar zenith angle to ∫E0(λ)cosθ s (λ)dλ integral gives 3.75×10 2 Wm -2 , so the normalized reflectance of the pixel is calculated as ρ n =0.512.

[0080] Further read the reflectivity of the four central 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. Spectral collapse vector Euclidean distance To calculate the curvature amplitude, five discrete reflectivities {0.29, 0.38, 0.34, 0.30, 0.28} were extracted in the range of 0.66 μm to 1.06 μm, and the five-point center second-order difference formula was substituted with a step size of h = 0.1 μm to obtain The spectral collapse-curvature index (SCCI) for this pixel is 0.157(1+0.48) = 0.232. The line-of-sight displacement of -0.056 m on the same day was then extracted from the InSAR time-difference field and projected at an incident angle of 38° to obtain the vertical instantaneous subsidence s = -0.082 m.

[0081] The discrete Laplacian kernel [0, 1, 0; 1, -4, 1; 0, 1, 0] / δ is used on a 3×3 pixel window. 2 Convolution At the same time, the horizontal gradient is obtained by central difference The mining geological survey report gives the average density of the limestone-mudstone group overburden ρ rock =2.45×10 3 kg\,m^{-3}. Substitute into the formula for the peak value of elastic strain energy density in The obtained value is |-0.082+(-0.035)|=0.117m, and the gradient mode is Multiply by 2.45 × 10 3 ×9.806 / 2 to obtain SCSED = 5.01 Jm -2 .

[0082] Then, using the 0.15m image, we performed Canny extraction with σ=1.2 in OpenCV to obtain 22,261 edge lines. The probabilistic Hough detection set the minimum line length to 50 pixels and the maximum discontinuity to 3 pixels to identify 143 cracks with a total length of 2,348m. The directions were concentrated in 85° (110m), 88° (950m), and 91° (670m). The remaining length of 618m was evenly distributed in the remaining 177 direction boxes. Normalization yielded P t (85°)=0.047, P t (88°)=0.405、P t (91°) = 0.285, and the remaining direction boxes are about 0.00148 each. With the summation of the second-order angular moment entropy After normalization, CAEM=1+2.86 / ln180=1.73. Vegetation index statistics are calculated using Sentinel-2B4 and B8 on the same day to calculate NDVI grids, and the average value for each period is Standard deviation σ NDVI，t =0.12, target pixel NDVI 0.16, calculate z-score square [(0.16-0.31) / 0.12] 2 =1.562.

[0083] The final comprehensive dynamic range impairment index is written as Compared with the background pixels far away from the mining area in the same image with an average CDRDI of 0.004 and a standard deviation of 0.006, the confidence level of this pixel reached z = (0.0334-0.004) / 0.006 = 4.9, which was significantly higher than the 3σ threshold and was automatically classified as the highest level of the three-level disaster damage. After visual rendering, a brightness focus patch was formed in this pixel and its surroundings, and the position coincided with the settlement gradient zone and the unidirectional crack bundle. The on-site RTK re-measurement of the crack width of 10-15 cm and the vertical displacement of the slope of 9 cm verified that the comprehensive indicators accurately revealed the high-risk slope section.

[0084] The 3σ threshold is a common criterion derived from the normal distribution: if a statistic's value falls outside the interval of three standard deviations σ above and below the population mean μ, the value is considered "extreme" or "anomalous." Formulated as μ±3σ, the threshold interval is: μ±1σ. In a perfectly normal distribution, approximately 68.27% of the samples fall within μ±1σ, approximately 95.45% fall within μ±2σ, and 99.73% fall within μ±3σ. Therefore, the probability of falling outside this interval is only approximately 0.27%, which is generally considered a highly abnormal or significant event in statistics.

[0085] In the application of the comprehensive dynamic range damage index, the CDRDI of the background pixels in the whole mining area is regarded as a set of approximately normal reference distributions, and its mean is used as the and standard deviation σ C Construction threshold If the CDRDI of a pixel is greater than the upper limit, it means that it is three times higher than the average level of the mining area and falls in the extreme tail of 0.13%, which is likely to correspond to high-risk damage areas such as collapse centers, crack-intensive zones or unstable edges of spoil dumps. On the contrary, if it is lower than In the context of damage analysis, these values often have no physical meaning and can be eliminated as abnormally low values or even noise. Using a 3σ threshold, significant damage boundaries can be automatically delineated for the entire image without the need for manual empirical coefficients, enabling objective and repeatable grade determination.

[0086] As described above, the above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit the same. Although the present invention has been described in detail with reference to the above embodiments, those skilled in the art should understand that the technical solutions described in the above embodiments can still be modified, or some of the technical features thereof can be replaced by equivalents. However, these modifications or replacements do not deviate the essence of the corresponding technical solutions 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 characteristic images of damaged landforms in mining areas, characterized in that: The method comprises: Step 1: Integrate the multispectral radiance within the set wavelength range and deduct two atmospheric scattering terms at each wavelength point to obtain the corrected reflectance; Step 2: Based on the corrected reflectance, extract the central reflectance of the green, red, near-infrared, and short-wave infrared bands; calculate the difference between the central reflectance of the green band and the central reflectance of the red band to obtain a first reflectance difference; calculate the difference between the central reflectance of the near-infrared band and the central reflectance of the short-wave infrared band to obtain a second reflectance difference; then calculate the Euclidean distance between 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 / chlorosis; Step 3: Based on the Laplace operator and the first-order gradient norm of the instantaneous subsidence, the coupling of the subsidence amplitude and deformation curvature is obtained. Combined with the measured density of the overburden, the peak value of the elastic strain energy density per unit area of the overburden-coal rock mass is calculated; Step 4: Extract all linear cracks from the remote sensing images of the same time period and calculate their directional histograms. Normalize the directional frequency of each degree into a probability. Then use the Shannon entropy formula to calculate the anisotropic entropy modulation of the cracks. This is used to construct a comprehensive dynamic range damage index to characterize the degree of damage of each pixel.

2. A method for extracting and analyzing dynamic range of characteristic images of damaged landforms in mining areas, characterized in that: In step 1, the wavelength range is set to 0.45 μm to 2.20 μm.

3. The image dynamic range extraction and analysis method according to claim 2, wherein: Corrected reflectivity ρ at wavelength λ n (λ) is: Among them, L TOA (λ) is the multispectral radiance at wavelength λ, in Wm -2 sr -1 μm -1 ;L path (λ) is the atmospheric path scattered radiation with wavelength λ, in Wm -2 sr -1 μm -1 ; τ d (λ) is the aerosol optical depth at wavelength λ, dimensionless; L sky (λ) is the sky backscattered radiation of wavelength λ, in Wm -2 sr -1 μm -1 ; sec represents the secant function; θ s is the solar zenith angle of wavelength λ, in rad; θ v is the observation zenith angle of wavelength λ, in rad.

4. The image dynamic range extraction and analysis method according to claim 3, wherein: The central wavelength of the green band is G, which is 0.56 μm; the central wavelength of the red band is R, which is 0.66 μm; the central wavelength of the near-infrared band is NIR, which is 0.86 μm; the central wavelength of the short-wave infrared band is SW. 1.6 It is 1.6μm.

5. The image dynamic range extraction and analysis method according to claim 4, wherein: The spectral collapse-curvature index SCCI is: Among them, ρ n (R) is the central reflectivity of the red band; ρ n (G) is the central reflectivity of the green band; ρ n (NIR) is the central reflectivity of the near-infrared band; ρ n (SW 1.6 ) is the central reflectivity of the shortwave infrared band; Represents the curvature amplitude in the near-infrared band.

6. The image dynamic range extraction and analysis method according to claim 5, wherein: The instantaneous subsidence s is calculated using the InSAR or GNSS points during the same period.

7. The image dynamic range extraction and analysis method according to claim 6, wherein: The peak elastic strain energy density SCSED of the pixel (x, y) is: Among them, ρ rock is the average density of the overburden; g is the acceleration of gravity; is the Laplace operator; δ is the pixel resolution of InSAR or GNSS; is the derivative of the instantaneous sinking amount s in the x-axis direction; is the derivative of the instantaneous sinking amount s in the y-axis direction; is the first-order gradient norm of the instantaneous sinking amount s.

8. The image dynamic range extraction and analysis method according to claim 7, wherein: In step 4, the Canny-Hough combination algorithm is used to extract all linear cracks from the remote sensing image of the same time period t and calculate their directional histograms, and the directional frequency of each degree is normalized into a probability. Specifically, the following steps are performed: 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 corresponding directional angles; dividing the 0 to 180° interval into 180 directional bins, and then calculating which directional bin the directional angle of each line falls in to construct a directional histogram; based on the directional histogram, the directional frequency is normalized to a directional probability P with a sum of 1. t (θ); the fracture anisotropy entropy modulation CAEM is:

9. The image dynamic range extraction and analysis method according to claim 8, wherein: In step 4, the comprehensive dynamic range damage index CDRDI(x, y) of the pixel (x, y) is: in, is the mean of the normalized vegetation index of all pixels in time period t; NDVI(x, y) is the normalized vegetation index of pixel (x, y); σ NDVI，t is the standard deviation of the normalized vegetation index for all pixels in time period t.

Images