Mining area illegal expansion remote sensing change monitoring method based on sparse subspace clustering

Abstract

The invention discloses a mining area illegal expansion remote sensing change monitoring method based on sparse subspace clustering, and relates to the technical field of mining areas, and the method comprises the following steps: obtaining mining area image data, carrying out the preprocessing of the mining area image data, obtaining consistent image data, and extracting an image feature tensor; constructing a self-expression model with sparse constraint based on the image feature tensor, performing subspace division, and extracting a target change region set according to a subspace division result; performing spatial superposition comparison on the target change area set and pre-obtained mining area legal boundary data to generate a newly added change area set, and performing index extraction to obtain a change area quantitative index set; and judging an illegal expansion risk level of the newly added change area set, and generating a change monitoring report according to the illegal expansion risk level and the consistent image data. According to the invention, through illegal expansion of the risk marking layer, a supervision department is supported to carry out rapid positioning, grading management and graph and certificate output.

Description

Technical Field

[0001] The present invention relates to the technical field of mining areas, and in particular to a remote sensing change monitoring method for illegal expansion of mining areas based on sparse subspace clustering. Background Art

[0002] With the rapid development of remote sensing image acquisition and processing technology, time-series remote sensing images based on satellites or drones have been widely used in dynamic monitoring of land resources and ecological and environmental supervision. They play an irreplaceable and important role in the identification and evidence collection of illegal expansion of mining areas. Illegal expansion of mining areas refers to mining activities that are not within the legal boundaries and without obtaining legal mining licenses. They are highly concealed, sudden, and spread rapidly in space, posing a significant threat to land resources, the ecological environment, and regulatory order.

[0003] Currently, mainstream remote sensing change detection methods mainly include image change recognition methods based on threshold difference and land feature classification methods based on supervised learning. The former usually achieves preliminary extraction of change areas by setting fixed thresholds on the spectral band differences, vegetation index changes, or texture index changes of remote sensing images at different times. Although the method is simple and computationally efficient, it is prone to a large number of misjudgments and missed detections when faced with complex and mixed scenes of land features. The accuracy of detection results is greatly reduced in the presence of seasonal vegetation interference or inconsistent imaging conditions. The latter relies on a large number of labeled samples to construct a classifier and achieves the identification of change areas through supervised learning. However, in actual mining area remote sensing monitoring, samples are often scarce and the labeling cost is high. The change patches formed in the early stages of illegal expansion are small in scale, dispersed in distribution, and have fuzzy boundaries, which can easily lead to insufficient model training or poor generalization ability. In addition, although deep learning models have achieved certain results in image semantic segmentation tasks, they are extremely computationally expensive when applied to high-resolution and large-format images. They place strict requirements on GPU computing power and data preprocessing, making it difficult to meet the actual needs of rapid response and batch supervision in mining areas.

[0004] To sum up, there is an urgent need to propose innovative technical solutions in terms of algorithm mechanism, feature modeling and recognition logic.

[0005] Currently, no effective solutions have been proposed for the problems in related technologies. Summary of the Invention

[0006] In response to the problems in the related art, the present invention proposes a remote sensing change monitoring method for illegal expansion of mining areas based on sparse subspace clustering to overcome the above-mentioned technical problems existing in the existing related art.

[0007] To this end, the specific technical solutions adopted in the present invention are as follows:

[0008] A remote sensing change monitoring method for illegal expansion of mining areas based on sparse subspace clustering includes:

[0009] S1. Acquire mining area image data, pre-process the mining area image data to obtain consistent image data, and extract image feature tensors using the consistent image data;

[0010] S2. Construct a self-expression model with sparse constraints based on the image feature tensor, use the self-expression model to perform subspace division, and extract the target change region set based on the subspace division result;

[0011] S3. Spatially overlay and compare the target change area set with the pre-acquired legal boundary data of the mining area to generate a new change area set, and extract indicators from the new change area set to obtain a quantitative indicator set of the change area;

[0012] S4. Using the set of quantitative indicators of the changed areas, determine the illegal expansion risk level of the newly added set of changed areas, and generate a change monitoring report based on the illegal expansion risk level and the consistency image data.

[0013] Furthermore, the mining area image data is obtained and pre-processed to obtain consistent image data. The image feature tensor is extracted using the consistent image data, including:

[0014] S11, acquiring mining area image data, and preprocessing the mining area image data to obtain consistent image data;

[0015] S12, based on the spatial overlapping areas of the consistent image data at each time point, using a sliding window to perform sliding cropping on the consistent image data in a two-dimensional spatial coordinate domain to obtain an image block set;

[0016] S13, extracting image features based on the image block set, and performing cascade fusion on the extracted image features to obtain a unified feature representation vector, and stacking all the unified feature representation vectors to obtain an image feature tensor;

[0017] Image features include: spectral features, vegetation index features, water index features and synthetic aperture radar scattering features.

[0018] Furthermore, the spectral features are obtained by extracting the average pixel value of the image block corresponding to each band at each time point and concatenating them in chronological order;

[0019] The vegetation index feature is obtained by calculating the ratio of the difference and sum of the pixel values ​​of the near-infrared band and the red band of the image block at each time point;

[0020] The water index feature is obtained by calculating the ratio of the difference and sum of the pixel values ​​of the green band and the near-infrared band of the image block at each time point;

[0021] The synthetic aperture radar scattering characteristics are obtained by extracting the average backscattering intensity value of the image block under each polarization channel at each time point and concatenating them in series according to the time and channel sequence.

[0022] Furthermore, a self-expression model with sparse constraints is constructed based on the image feature tensor. The self-expression model is used to perform subspace division, and the target change region set is extracted according to the subspace division result, including:

[0023] S21. Setting a multi-scale window set based on the image feature tensor, and constructing a multi-scale image feature tensor for each scale in the multi-scale window set;

[0024] S22. Utilize the multi-scale image feature tensor to construct a sparse self-expression model under adaptive structural constraints, and obtain a multi-scale joint sparse coefficient matrix based on the sparse self-expression model;

[0025] S23. Perform subspace partitioning based on the multi-scale joint sparse coefficient matrix, and extract a target change region set according to the subspace partitioning result.

[0026] Furthermore, subspace partitioning is performed based on the multi-scale joint sparse coefficient matrix, and the target change region set is extracted according to the subspace partitioning result, including:

[0027] S231. Constructing a non-negative symmetric similarity matrix based on the multi-scale joint sparse coefficient matrix, and calculating the graph Laplacian matrix according to the non-negative symmetric similarity matrix;

[0028] S232, performing a spectral decomposition operation on the graph Laplacian matrix to obtain eigenvectors corresponding to the first K smallest eigenvalues, and concatenating the eigenvectors to form a spectral embedding matrix;

[0029] S233, performing a clustering operation on the spectral embedding matrix using K-means clustering to obtain a subspace category label set for each image block;

[0030] S234. Divide all subspaces in the subspace category label set into a stable background subspace set and a perturbation subspace set according to the compactness and stability of the sparse expression in each subspace;

[0031] S235 . Based on the image block index in the perturbation subspace set, extract the corresponding candidate change region set, and screen the candidate change region set in combination with the sparsity index and the spatial connectivity constraint to obtain the target change region set.

[0032] Furthermore, K-means clustering is used to perform clustering operations on the spectral embedding matrix, and the subspace category label set for each image block is obtained, including:

[0033] S2331, randomly selecting K data points from the spectral embedding matrix as initial cluster centers for K-means clustering according to a preset number of clusters K;

[0034] S2332, using Euclidean distance to calculate the distance between each data point in the spectral embedding matrix and the initial cluster center, and assigning each data point to the cluster to which the nearest cluster center belongs;

[0035] S2333, recalculating the average value of all points in each cluster to iteratively update the cluster center until the cluster center reaches a preset maximum number of iterations;

[0036] S2334. Assign a corresponding subspace category label to each image block in the spectrum embedding matrix to obtain a subspace category label set for each image block.

[0037] Furthermore, based on the image block index in the perturbation subspace set, the corresponding candidate change region set is extracted, and the candidate change region set is screened by combining the sparsity index and the spatial connectivity constraint. The target change region set includes:

[0038] S2351. Extracting a corresponding set of candidate changed regions based on the image block indexes in the perturbation subspace set;

[0039] S2352. Calculate a sparsity index for each image block in the candidate change region set, and based on a preset sparsity index threshold and spatial connectivity constraints, select image blocks that meet the conditions from the candidate change region set to obtain a target change region set.

[0040] Furthermore, the target change area set is spatially overlaid and compared with the pre-acquired legal boundary data of the mining area to generate a new change area set. The indicators of the new change area set are extracted to obtain a set of quantitative indicators of the change area, including:

[0041] S 31. Based on the two-dimensional polygon Boolean difference operation, the target change area set is spatially superimposed and compared with the pre-acquired legal boundary data of the mining area to generate a new change area set;

[0042] S32, extracting indicators from the newly added set of changed regions to obtain a set of quantitative indicators of the changed regions;

[0043] The set of quantitative indicators of the change area includes: geometric shape indicator, actual area indicator, expansion direction indicator and change intensity indicator.

[0044] Furthermore, the geometric shape index is obtained by taking the ratio of the actual area of ​​the newly added change region to the square of its boundary length and multiplying it by a constant factor;

[0045] The actual area index is obtained by counting the number of pixels contained in the newly added changed area and multiplying it by the spatial resolution of the image data in the horizontal direction and the spatial resolution in the vertical direction;

[0046] The expansion direction index is obtained by calculating the angular relationship between the central moments of the newly added change areas;

[0047] The change intensity index is obtained by calculating the difference between the unified feature representation vectors of the newly added change area at two time points before and after the change, and averaging the difference values ​​within the newly added change area.

[0048] Furthermore, the set of quantitative indicators of the changed areas is used to determine the illegal expansion risk level of the newly added set of changed areas, and a change monitoring report is generated based on the illegal expansion risk level and the consistent image data, including:

[0049] S41. Defining the illegal expansion risk level based on the joint distribution characteristics of the set of quantitative indicators of the change area;

[0050] S42. Determine the illegal expansion risk level of each newly changed area according to the defined illegal expansion risk level, and mark the risk level of each newly changed area using the illegal expansion risk level determination result to generate an illegal expansion risk marking layer.

[0051] S43. Perform three-dimensional visualization fusion of the illegal expansion risk marker layer and the consistent remote sensing image data to generate an interactive change map, and output a change monitoring report based on the interactive change map.

[0052] The beneficial effects of the present invention are:

[0053] 1. The present invention achieves collaborative modeling of the sparsity and local anomalies of illegal expansion patches through a sparse expression optimization strategy with local anomaly priors and spatial structure preservation constraints. An expression enhancement strategy is implemented for suspicious patches by introducing anomaly prior scores in the sparse expression optimization process. At the same time, a spatial structure preservation term is constructed through the structural similarity of adjacent image blocks. The expression sparsity and neighborhood consistency of local mutation areas are simultaneously constrained in the optimization function, thereby avoiding misjudging vegetation interference or terrain occlusion as illegal expansion changes. Compared with traditional sparse subspace clustering, the accuracy and stability of anomaly identification are improved, and it performs significantly in early illegal expansion areas with blurred boundaries or weak changes.

[0054] 2. The present invention performs Laplace spectral decomposition on the multi-scale sparse coefficient matrix and combines it with K-means clustering to realize subspace partitioning. The candidate areas extracted from the disturbance subspace are quantitatively analyzed based on expression sparsity, geometric shape, area, change intensity and directional consistency indicators, and a regularized illegal expansion risk grading model is constructed. In the newly added illegal areas outside the historical boundaries, small patches, high directional deviations or multiple disturbance features can be significantly identified, and an automated illegal expansion risk marker layer is generated to support regulatory authorities in rapid positioning, grading management and graphic output.

[0055] 3. This method extracts a composite feature vector consisting of spectral response, normalized difference vegetation index, normalized difference water index, and synthetic aperture radar scattering intensity. It constructs a remote sensing feature tensor using a multi-scale window. Based on this, it introduces a scale-adaptive sparse self-expression model to achieve sparse reconstruction of the heterogeneous features of land features at different spatial scales. Compared with traditional single-scale sparse representation methods, this method can more stably extract the heterogeneity of illegal expansion areas in mining scenarios with high interference and mixed features, improving sensitivity and robustness to regional structural changes. BRIEF DESCRIPTION OF THE DRAWINGS

[0056] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the following briefly introduces the drawings required for use in the embodiments. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying any creative work.

[0057] Figure 1 is a flow chart of a method for monitoring changes in illegal expansion of mining areas using remote sensing based on sparse subspace clustering according to an embodiment of the present invention;

[0058] Figure 2 It is a flow chart of the execution steps in the method for monitoring remote sensing changes in illegal expansion of mining areas based on sparse subspace clustering according to an embodiment of the present invention. DETAILED DESCRIPTION

[0059] To further illustrate each embodiment, the present invention provides drawings, which are part of the disclosure of the present invention. They are mainly used to illustrate the embodiments and can be used in conjunction with the relevant descriptions in the specification to explain the operating principles of the embodiments. By referring to these contents, ordinary technicians in this field should be able to understand other possible implementation methods and the advantages of the present invention.

[0060] According to an embodiment of the present invention, a method for monitoring changes in illegal expansion of mining areas through remote sensing based on sparse subspace clustering is provided.

[0061] The present invention will now be further described with reference to the accompanying drawings and specific embodiments. Figure 1-2 As shown, according to an embodiment of the present invention, a remote sensing change monitoring method for illegal expansion of mining areas based on sparse subspace clustering includes the following steps:

[0062] S1. Acquire mining area image data, pre-process the mining area image data to obtain consistent image data, and use the consistent image data to extract image feature tensors.

[0063] Specifically, the mining area image data is obtained and preprocessed to obtain consistent image data. The image feature tensor is extracted using the consistent image data, including:

[0064] S 11. Acquire mining area image data, and pre-process the mining area image data to obtain consistent image data;

[0065] S12. Based on the spatial overlapping areas of the consistent image data at each time point, use a sliding window to perform sliding cropping on the consistent image data in a two-dimensional spatial coordinate domain to obtain an image block set;

[0066] S 13. Extract image features based on the image block set, perform cascade fusion on the extracted image features to obtain a unified feature representation vector, and stack all the unified feature representation vectors to obtain an image feature tensor;

[0067] Image features include: spectral features, vegetation index features, water index features and synthetic aperture radar scattering features.

[0068] Specifically, the spectral features are obtained by extracting the average pixel value of the image block corresponding to each band at each time point and concatenating them in chronological order;

[0069] The vegetation index feature is obtained by calculating the ratio of the difference and sum of the pixel values ​​of the near-infrared band and the red band of the image block at each time point;

[0070] The water index feature is obtained by calculating the ratio of the difference and sum of the pixel values ​​of the green band and the near-infrared band of the image block at each time point;

[0071] The synthetic aperture radar scattering characteristics are obtained by extracting the average backscattering intensity value of the image block under each polarization channel at each time point and concatenating them in series according to the time and channel sequence.

[0072] Specifically, the mining area image data is multi-temporal remote sensing image data covering the target mining area and its surrounding areas. Radiometric normalization processing, terrain correction processing and joint registration processing are performed on the multi-temporal remote sensing image data to obtain consistent remote sensing image data with a unified coordinate reference and spectral response (i.e., consistent image data).

[0073] Specifically, based on the consistency of remote sensing image data I (t) (x, y, b) In the spatial overlapping area at each time point t, the consistent remote sensing image data is cut according to the sliding window on the two-dimensional spatial coordinate domain (x, y) to generate an image block set B. Each image block B in the image block set i Contains multi-band and multi-temporal information of fixed-size sliding window pixels; establishes an image block index matrix M, where each row in the image block index matrix records the corresponding image block B i The center two-dimensional space coordinate domain (x i ,y i ).

[0074] Specifically, the spectral feature vectors of each image block in the image block set at all time points and all bands are calculated. The spectral feature vectors are used to represent the spectral response change trend of the image block in the multi-temporal remote sensing image. The spectral feature vectors are obtained by extracting the average pixel value of the image block corresponding to each band at each time point and concatenating them in chronological order.

[0075] Specifically, the normalized difference vegetation index feature vector of each image block in the image block set at all time points is calculated. The normalized difference vegetation index feature vector is used to reflect the changes in surface vegetation coverage in the image block area. The normalized difference vegetation index feature vector is obtained by calculating the ratio of the difference and sum of the pixel values ​​of the near-infrared band and the red light band of the image block at each time point. The pixel values ​​of the near-infrared band and the red light band are both the average pixel values ​​of the image block in the band.

[0076] Specifically, the normalized difference water index feature vector of each image block in the image block set at all time points is calculated. The normalized difference water index feature vector is used to reflect the existence and change of water bodies in the image block area. The normalized difference water index feature vector is obtained by calculating the ratio of the difference and sum of the pixel values ​​of the green band and the near-infrared band of the image block at each time point, where the green band pixel value and the near-infrared band pixel value are both the average pixel values ​​of the image block in the corresponding bands.

[0077] Specifically, under the condition that the consistent remote sensing image data contains synthetic aperture radar information, the synthetic aperture radar scattering feature vector of each image block in the image block set at all time points and all polarization channels is calculated. The synthetic aperture radar scattering feature vector is used to enhance the response ability of the image block to the bare rock, stockpile and waste rock area in the mining area. The synthetic aperture radar scattering feature vector is obtained by extracting the average backscattering intensity value of the image block under each polarization channel at each time point, and concatenating them in series according to the time and channel order.

[0078] Specifically, the spectral feature vector f ispec , vegetation index characteristic vector f i ndvi , water body index characteristic vector f i ndwi With the radar characteristic vector f i sar Perform cascade fusion to construct a unified feature representation vector f i .

[0079] Specifically, for all image blocks B i The unified feature representation vector f i The images are stacked to form a remote sensing feature tensor F (i.e., image feature tensor).

[0080] Specifically, the overlapping areas of the consistent remote sensing image data are divided into image block sets according to a fixed size, and an image block index is established to maintain the spatial correspondence. The spectral features, vegetation index features, water index features and synthetic aperture radar scattering features of the image block set are extracted and fused to form a remote sensing feature tensor.

[0081] S2. Construct a self-expression model with sparse constraints based on the image feature tensor, use the self-expression model to perform subspace division, and extract the target change area set according to the subspace division result.

[0082] Specifically, a self-expression model with sparse constraints is constructed based on the image feature tensor, and subspace division is performed using the self-expression model. The target change region set is extracted based on the subspace division result, including:

[0083] S21. Setting a multi-scale window set based on the image feature tensor, and constructing a multi-scale image feature tensor for each scale in the multi-scale window set;

[0084] S22. Utilize the multi-scale image feature tensor to construct a sparse self-expression model under adaptive structural constraints, and obtain a multi-scale joint sparse coefficient matrix based on the sparse self-expression model;

[0085] S23. Perform subspace partitioning based on the multi-scale joint sparse coefficient matrix, and extract a target change region set according to the subspace partitioning result.

[0086] Specifically, subspace partitioning is performed based on the multi-scale joint sparse coefficient matrix, and the target change region set is extracted according to the subspace partitioning result, including:

[0087] S231. Constructing a non-negative symmetric similarity matrix based on the multi-scale joint sparse coefficient matrix, and calculating the graph Laplacian matrix according to the non-negative symmetric similarity matrix;

[0088] S232, performing a spectral decomposition operation on the graph Laplacian matrix to obtain eigenvectors corresponding to the first K smallest eigenvalues, and concatenating the eigenvectors to form a spectral embedding matrix;

[0089] S233. Perform a clustering operation on the spectral embedding matrix using K-means clustering to obtain a subspace category label set for each image block.

[0090] Specifically, the spectral embedding matrix is ​​clustered using K-means clustering to obtain the subspace category label set for each image block, including:

[0091] S2331, randomly selecting K data points from the spectral embedding matrix as initial cluster centers for K-means clustering according to a preset number of clusters K;

[0092] S2332, using Euclidean distance to calculate the distance between each data point in the spectral embedding matrix and the initial cluster center, and assigning each data point to the cluster to which the nearest cluster center belongs;

[0093] S2333, recalculating the average value of all points in each cluster to iteratively update the cluster center until the cluster center reaches a preset maximum number of iterations;

[0094] S2334. Assign a corresponding subspace category label to each image block in the spectrum embedding matrix to obtain a subspace category label set for each image block.

[0095] S234. Divide all subspaces in the subspace category label set into a stable background subspace set and a perturbation subspace set according to the compactness and stability of the sparse expression in each subspace;

[0096] S235 . Based on the image block index in the perturbation subspace set, extract the corresponding candidate change region set, and screen the candidate change region set in combination with the sparsity index and the spatial connectivity constraint to obtain the target change region set.

[0097] Specifically, based on the image block index in the perturbation subspace set, the corresponding candidate change region set is extracted, and the candidate change region set is screened by combining the sparsity index and the spatial connectivity constraint. The target change region set includes:

[0098] S2351. Extracting a corresponding set of candidate changed regions based on the image block indexes in the perturbation subspace set;

[0099] S2352. Calculate a sparsity index for each image block in the candidate change region set, and based on a preset sparsity index threshold and spatial connectivity constraints, select image blocks that meet the conditions from the candidate change region set to obtain a target change region set.

[0100] Specifically, a self-expression model with sparse constraints is constructed with the remote sensing feature tensor as input, and a sparse coefficient matrix is ​​obtained by solving it. The sparse coefficient matrix is ​​used to characterize the low-rank correlation structure between image blocks.

[0101] Specifically, based on the remote sensing feature tensor F, a multi-scale window set S is set, and for each scale s in the multi-scale window set k Construct the corresponding multi-scale remote sensing feature tensor Multi-scale remote sensing feature tensors are composed of image blocks at scale s k The unified feature representation vectors generated under the

[15] are stacked to characterize the structural distribution characteristics of illegal mining expansion at different spatial granularities.

[0102] Specifically, for each scale s k , based on multi-scale remote sensing feature tensor Construct a sparse self-expression model under adaptive structural constraints, and the objective function is:

[0103]

[0104] in, Indicates scale s k The lower sparse coefficient matrix; Indicates scale s k The sparse coefficient vector in the next i-th row also represents the linear expression contribution of any image block in the image block set to the other image blocks in the feature space; N represents the number of image blocks, Indicates scale s k The error matrix below also represents the residual term between the actual multi-scale remote sensing feature tensor and the sparse reconstruction, which is used to accommodate the unstructured disturbances caused by shadows, water reflections, and terrain occlusions in mining area remote sensing; Indicates scale s k The L1 norm of the sparse coefficient vector in the next i-th row is used to measure the image block B i In scale s k The sparse expression intensity under the condition of α is sparse, the smaller the value, the more concentrated the expression and the more significant the change; a i Represents image block B i The local anomaly prior score is obtained by calculating the difference between the NDVI, NDWI or SAR features and the phase average value, which is used to reflect the possibility of illegal changes in the image block; η represents the anomaly adjustment factor, which controls the prior anomaly score a i The sensitivity of the weighted sparse terms. The larger the value, the stronger the activation of the suspicious region. λ1 represents the robust error weight factor. represents the L2,1 norm of the error matrix, and also represents the overall strength of the error on each immediate image block, which encourages the error to concentrate in a specific area rather than diffuse globally, and helps to improve robustness. λ2 represents the structure preservation weight factor; N i Represents image block B i The spatial adjacent set of is defined as the image block B in geographic space. i The set of all image blocks that are adjacent or separated by a distance less than the set threshold; Indicates scale s k Lower image block B i With image block B j The adaptive spatial structure weights of With sparse coefficient vector The structural similarity of Represents image block B i and image block B j The square of the Euclidean distance between the expression coefficient vectors; Indicates that the diagonal of the sparse coefficient matrix is ​​zero, preventing the image block from expressing itself.

[0105] Specifically, the first term is the weighted L1 norm term (1+η·a i )||C i ||1, the weighting mechanism can strengthen the structural activation of illegal expansion areas during the optimization process, so that they obtain higher expression weights in the sparse subspace, thereby improving the model's sensitivity to weak change areas.

[0106] The second term is the reconstruction error term of the L2,1 norm It is used to enhance the robustness to local unstructured noise (terrain shadows, atmospheric disturbances, water reflections) and reduce its interference with the global expression structure.

[0107] The third term is the graph regularity constraint where N i Represents image block B i The set of spatially adjacent blocks of For its adjacent image block B j The regularization term is used to maintain the consistency of geographically adjacent areas in the sparse expression structure, suppress boundary fuzziness and local jump phenomena, and improve the spatial coherence of the partitioning results.

[0108] The sparse coefficient matrices corresponding to each scale in all multi-scale settings are weightedly fused to form a unified multi-scale joint sparse coefficient matrix. The multi-scale joint sparse coefficient matrix is ​​used to comprehensively express the sparse expression structure between image blocks at different spatial scales. Each element of the multi-scale joint sparse coefficient matrix is ​​obtained by the weighted sum of the corresponding sparse coefficient values ​​at all scales and its scale fusion weight. The scale fusion weight is a positive real number, and the sum of all scale fusion weights is 1, which is used to balance the contribution of different scales to the overall expression structure.

[0109] Multi-scale joint sparse coefficient matrix C fusion As input, a non-negative symmetric similarity matrix W is constructed. The non-negative symmetric similarity matrix is ​​obtained by symmetric normalization of the multi-scale joint sparse coefficient matrix:

[0110]

[0111] Among them, |C fusion | represents the multi-scale joint sparse coefficient matrix C fusion The absolute value of each element in is taken to ensure that the similarity is a non-negative real number. The similarity matrix W is used to quantify the mutual similarity between any two image blocks in the image block set in the sparse expression space. The corresponding graph Laplacian matrix L is calculated based on the similarity matrix W.

[0112] Perform spectral decomposition on the graph Laplace matrix L, calculate the eigenvectors corresponding to the first K smallest eigenvalues ​​of the graph Laplace matrix, and concatenate them to form the spectral embedding matrix Y∈R N×K , where K represents the number of preset subspaces, and each row of Y i ∈RK represents the image block B i Position in the low-dimensional spectral embedding space.

[0113] Perform K-means clustering (i.e., K-means clustering) on ​​the spectral embedding matrix Y to obtain the subspace category label set L for each image block. i ∈{1, 2, ..., K} represents the image block B i The subspace number to which it belongs is used, and all subspaces in the subspace category label set are divided into a stable background subspace set and a perturbation subspace set according to the compactness and stability of the sparse expression in each subspace.

[0114] According to the image block index in the perturbation subspace set, the corresponding candidate change region set is extracted. The candidate change region set is used to characterize the regional blocks showing abnormal structural changes in the remote sensing features.

[0115] For each image block B in the candidate change region set j , calculate the sparsity index ρj , the sparsity index is used to measure the concentration and abnormality of image blocks in the expression structure:

[0116]

[0117] Among them, ρ j represents the j-th row sparse coefficient vector in the multi-scale joint sparse coefficient matrix, ||·||1 represents the L1 norm; ||·||2 represents the L2 norm; ∈ represents a small value to prevent the denominator from being zero.

[0118] The sparsity index threshold τ ρ Based on the spatial connectivity constraint, the image blocks that meet the conditions are screened out from the candidate change area set to form the suspected change area set (i.e., the target change area set). Among them B j represents the jth suspected change image block obtained by sparse subspace partitioning and sparsity screening, N sus Represents the number of suspected change image blocks, where the sparsity index ρ j <τ ρ And image block B j It is adjacent to other low-sparse image blocks in terms of spatial connectivity constraints, and the adjacent blocks satisfy the Euclidean distance less than the threshold δ.

[0119] S3. Spatially overlay and compare the target change area set with the pre-acquired legal boundary data of the mining area to generate a new change area set, and extract indicators from the new change area set to obtain a quantitative indicator set of the change area.

[0120] Specifically, the target change area set is spatially superimposed and compared with the pre-acquired legal boundary data of the mining area to generate a new change area set. The indicators of the new change area set are extracted to obtain a set of quantitative indicators of the change area, including:

[0121] S 31. Based on the two-dimensional polygon Boolean difference operation, the target change area set is spatially superimposed and compared with the pre-acquired legal boundary data of the mining area to generate a new change area set;

[0122] S32, extracting indicators from the newly added set of changed regions to obtain a set of quantitative indicators of the changed regions;

[0123] The set of quantitative indicators of the change area includes: geometric shape indicator, actual area indicator, expansion direction indicator and change intensity indicator.

[0124] Specifically, the geometric shape index is obtained by taking the ratio of the actual area of ​​the newly added change region to the square of its boundary length and multiplying it by a constant factor;

[0125] The actual area index is obtained by counting the number of pixels contained in the newly added changed area and multiplying it by the spatial resolution of the image data in the horizontal direction and the spatial resolution in the vertical direction;

[0126] The expansion direction index is obtained by calculating the angular relationship between the central moments of the newly added change areas;

[0127] The change intensity index is obtained by calculating the difference between the unified feature representation vectors of the newly added change area at two time points before and after the change, and averaging the difference values ​​within the newly added change area.

[0128] Specifically, obtain the mining area legal boundary dataset B legal ; Set the suspected change region R sus and mining area legal boundary dataset B legal Perform spatial overlay analysis, which uses two-dimensional polygon Boolean difference operations to extract new areas outside any legal boundary, which are defined as the set of newly changed areas R new , the new change area set space logic expression is:

[0129] R new =R sus \(B plan ∪B permit ∪B history );

[0130] The set operator \ represents the spatial difference operation, and ∪ represents the spatial union operation, which is used to eliminate the changed areas within the legal range and only retain the expanded parts outside the illegal boundary.

[0131] The mining area legal boundary dataset includes: mining area planning red line data B plan , indicating the legal mining boundary approved by the corresponding department; mining license boundary data B permit , indicating the spatial scope of mining licenses actually mined by mining rights holders; historical mining boundary data B history , indicating the outline of the legal mining area confirmed in the historical remote sensing imagery.

[0132] Specifically, the geometric shape index of each newly added change area in the set of newly added change areas is calculated. The geometric shape index is used to measure the degree of regularity of the outline of the newly added change area. The geometric shape index is obtained by obtaining the ratio of the actual area of ​​the newly added change area to the square of its boundary length and multiplying it by a constant factor.

[0133] Specifically, the actual area index of each newly changed area in the set of newly changed areas is calculated. The actual area index is used to measure the physical expansion scale of the newly changed area in the remote sensing image. The actual area index is obtained by counting the number of pixels contained in the newly changed area and multiplying it by the spatial resolution in the horizontal direction and the spatial resolution in the vertical direction of the remote sensing image.

[0134] Specifically, the expansion direction index of each newly changed area in the set of newly changed areas is calculated. The expansion direction index is used to indicate whether the newly changed area has a significant directional extension trend in space. The expansion direction index is obtained by calculating the angular relationship between the central moments of the newly changed areas.

[0135] Specifically, the change intensity index of each newly added change area in the set of newly added change areas is calculated. The change intensity index is used to measure the overall change amplitude of the newly added change area in remote sensing characteristics. The change intensity index is obtained by calculating the difference between the unified feature representation vectors of the area at two time points before and after the change, and averaging the difference values ​​in the area.

[0136] Specifically, according to the shape index φ k , area index A k , expansion direction index θ k and the change intensity index Δ k Based on the joint distribution characteristics of the two regions, the illegal expansion risk levels are defined as low-risk area, medium-risk area and high-risk area.

[0137] Specifically, each newly added change area B k The risk levels are marked as low-risk areas, medium-risk areas and high-risk areas, and an illegal expansion risk marker layer L is constructed. risk ,The illegal expansion risk marking layer is used to mark the spatial location of each newly added change area and the corresponding risk level label.

[0138] Specifically, low-risk areas, medium-risk areas, and high-risk areas are defined as follows:

[0139] The low-risk area is the area that meets the shape index φ k ≥0.75, and area index A k ≤300m 2 , and the change intensity index Δ k ≤0.25·Δ max , and the expansion direction indicator The medium risk area is the area that meets the shape index 0.5≤φ k <0.75 or area index 300m 2 k ≤800mm 2 or change intensity index 0.25·Δ​max <Δ k ≤0.6·Δ max or expansion direction indicator The high-risk area is the area that meets the shape index φ k <0.5 or area index A k >800m 2 or change intensity index Δ k >0.6·Δ max or expansion direction indicator Among them, θ ref ∈[0,π) represents the reference value of the typical expansion direction of the historical mining area, which is used to determine whether the newly added area shows abnormal directional deviation.

[0140] Shape index φ k The threshold is 0.75 / 0.5, and its physical meaning is that the shape index is used to measure the regularity of the outline of the region. k =1 means a perfect circle, φ k →0 indicates a more complex shape, a more "straggling" or "sprawling" shape. The criteria for this are: a value above 0.75 indicates a compact area, typically a natural feature or engineering control boundary, and is considered low risk; a value between 0.5 and 0.75 indicates a complex edge with some ductility, potentially indicating temporary or unplanned construction; and a value below 0.5 indicates linear expansion, serpentine, or scattered patterns, often indicating informal encroachment, and is therefore considered high risk.

[0141] Area index A k The threshold is 300m 2 / 800m 2 , the physical meaning is the absolute coverage area of ​​the illegal expansion area. Setting basis (taking the remote sensing image resolution of 0.5-2.0m as an example): 300m 2 ≈30m×10m, equivalent to a newly opened small road or simple working zone; 800m 2 ≈New operating boundaries for a single block in a standard small or medium-sized open-pit mining area.

[0142] Change intensity index Δ k The threshold ratio is 0.25 / 0.6. Its physical meaning is to measure the degree of temporal variation in the remote sensing feature space of the area. Larger values ​​indicate more dramatic changes. The setting is based on the change value normalized to a range of 0-1 and then divided by quantiles. The top 25% is low-fluctuation, possibly due to seasonal changes; the 25-60% is intermediate-variability, indicating risk but requiring additional judgment; and values ​​above 60% indicate dramatic changes, typically due to human excavation or stockpiling.

[0143] Expansion direction indicator θ k The deviation angle threshold and The physical meaning is the degree of deviation between the main expansion direction of the newly added area and the main direction of the historical mining area. The setting basis is θ k and the historical mining area direction θ ref The greater the deviation, the more likely the expansion behavior will deviate from the approved working surface; π / 12≈15°, π / 6≈30°; according to the topographic engineering specifications and the mining area planning red line control line, the "error band" is allowed to be set. Generally, the error of human-controlled mining should not be greater than 15°. If it exceeds 30°, it is very likely that a new path has been opened or soil is taken in violation of regulations.

[0144] S4. Using the set of quantitative indicators of the changed areas, determine the illegal expansion risk level of the newly added set of changed areas, and generate a change monitoring report based on the illegal expansion risk level and the consistency image data.

[0145] Specifically, the illegal expansion risk level of the newly added change area set is determined using the change area quantitative indicator set, and a change monitoring report is generated based on the illegal expansion risk level and consistent image data, including:

[0146] S41. Defining the illegal expansion risk level based on the joint distribution characteristics of the set of quantitative indicators of the change area;

[0147] S42. Determine the illegal expansion risk level of each newly changed area according to the defined illegal expansion risk level, and mark the risk level of each newly changed area using the illegal expansion risk level determination result to generate an illegal expansion risk marking layer.

[0148] S43. Perform three-dimensional visualization fusion of the illegal expansion risk marker layer and the consistent remote sensing image data to generate an interactive change map, and output a change monitoring report based on the interactive change map.

[0149] Specifically, the change monitoring report includes time information, spatial location, magnitude of change, risk level and enforcement recommendations.

[0150] Example:

[0151] A natural resource monitoring platform received two sets of remote sensing image data covering an open-pit coal mine in a certain area: the first was an optical image acquired on August 25, 2024, with a resolution of 0.8 meters and four spectral bands; the second was a SAR radar image acquired on September 5, 2024, with a resolution of 10 meters and a VV polarization channel.

[0152] Without any manually annotated data, the platform used the proposed method to monitor and analyze whether there was illegal expansion in the mining area. This scenario involved a complex regulatory environment characterized by multi-temporal remote sensing data, mixed ground spectral data, and a small sample size for sudden changes.

[0153] First, the remote sensing data preprocessing module was invoked to perform radiometric normalization and terrain correction on the images from August 25 and September 5, ensuring consistent coordinates and spectral responses for the same features across the time series. A fixed window sliding partitioning was then performed within the approximately 12 square kilometers of overlapping area, generating a total of 18,200 image blocks, each measuring 256×256 pixels. Spectral signatures, NDVI, NDWI, and SAR backscatter intensity were extracted for each block and concatenated into a 124-dimensional remote sensing feature vector, forming a remote sensing feature tensor.

[0154] Without using any training samples, a multi-scale sparse expression model is constructed, and the sparse coefficient matrix is ​​optimized with structure preservation and anomaly excitation. During the sparse optimization process, 47 image blocks have anomaly scores of a in the NDWI or SARVV channels. i It is significantly higher than other areas (more than 3 times the global mean). The sparse expression weight of the image block numbered #8821 is mainly dominated by distant non-adjacent sub-blocks, with a sparse L1 norm of 0.021, an L2 norm of 0.523, and a sparsity index ρ. 8821 It is only 0.04, which is significantly lower than the average level of 0.57 in the whole area, so it is marked as a block with significant structural abnormality.

[0155] Subsequently, the joint sparse coefficient matrix was spectrally constructed and Laplacian matrix decomposition was performed to obtain the spectral embedding space Y. K-means clustering was used to divide all image blocks into five subspaces. The fourth subspace contained a large number of abnormal blocks with a sparsity index lower than 0.2, which were distributed in the southwest of the main mining area. It was designated as a disturbance subspace and marked as a candidate change area.

[0156] Furthermore, the candidate change areas were spatially superimposed with the legal boundary data of the mining area. The results showed that the four change area blocks numbered #8710, #8821, #8837, and #8842 were all located between the planning red line Bplan and the historical mining boundary B history In addition, a new set of change regions R is generated based on this new .

[0157] The newly added area was calculated and risk assessed. The area numbered #8821 is 1,130 square meters, with a shape compactness of φ = 0.36 and a change intensity index of Δ = 0.84·Δ max , the expansion direction θ = 1.98 rad, which deviates from the historical expansion direction by 0.78 rad. The area is marked as a high-risk illegal expansion area.

[0158] Finally, the illegal expansion risk mark layer L is generated riskA total of 6 high-risk areas, 12 medium-risk areas, and 21 low-risk areas were identified, and the high-risk area numbered #8821 was pushed to the regulatory end in real time.

[0159] In this area, the traditional method A (NDVI difference, threshold set to 0.12), method B (based on supervised random forest, 1200 training samples) and method C (based on deep FCN semantic segmentation model, training set same as B) were used simultaneously. The comparison results of the above traditional methods and the method D of the present invention for recognition are shown in Table 1 below.

[0160] Table 1 Identification comparison results of the method of the present invention

[0161]

[0162] The present invention accurately identifies all newly changed areas without using any training samples. Traditional method A has two false omissions, and method B has a serious omission due to insufficient training samples. Although method C performs well, its training time is 8.5 hours and the running time is 15 minutes. The method of the present invention can complete the entire process analysis on an ordinary workstation in just 5 minutes, which is significantly better than the traditional methods in terms of efficiency and responsiveness.

[0163] This Example 1 truly demonstrates the ability of the present invention to achieve high-sensitivity, low-false-alarm, and full-process automated identification of illegal expansion of mining areas under the background of limited regulatory resources, missing samples, and hidden and complex changing behaviors. It has extremely high engineering application value and deployment and promotion potential.

[0164] Based on multi-temporal remote sensing imagery with a resolution better than 2 meters, this method extracts composite feature vectors including spectral response, normalized difference vegetation index, normalized difference water index, and synthetic aperture radar scattering intensity. A remote sensing feature tensor is constructed using a multi-scale window. Based on this, a scale-adaptive sparse self-expression model is introduced to achieve sparse reconstruction of the heterogeneous features of land features at different spatial scales. Compared with traditional single-scale sparse representation methods, this method can more stably extract the feature heterogeneity of illegal expansion areas in mining scenarios with high interference and mixed features, improving sensitivity and robustness to regional structural changes.

[0165] The present invention proposes a sparse expression optimization strategy with local anomaly prior and spatial structure preservation constraints to achieve collaborative modeling of the sparsity and local anomaly of illegal expansion patches. By introducing anomaly prior scores constructed based on NDVI, NDWI or SAR indicators in the sparse expression optimization process, an expression enhancement strategy is implemented for suspicious patches. At the same time, a spatial structure preservation term is constructed through the structural similarity of adjacent image blocks. In the optimization function, the expression sparsity and neighborhood consistency of local mutation areas are constrained at the same time, thereby avoiding misjudging vegetation interference or terrain occlusion as illegal expansion changes. Compared with traditional sparse subspace clustering, the accuracy and stability of anomaly recognition are improved, and it performs significantly in early illegal expansion areas with blurred boundaries or weak changes.

[0166] The present invention performs graph Laplacian spectral decomposition on a multi-scale sparse coefficient matrix and combines it with K-means clustering to achieve subspace partitioning. The candidate areas extracted from the disturbance subspace are quantitatively analyzed based on expression sparsity, geometric shape, area, change intensity and directional consistency indicators, and a regularized illegal expansion risk grading model is constructed. In newly added illegal areas outside the historical boundaries, small patches, high directional deviations or multiple disturbance features can be significantly identified, and an automated illegal expansion risk marker layer can be generated to support regulatory authorities in rapid positioning, grading management and graphic output.

[0167] In summary, with the help of the above technical solutions of the present invention, the present invention realizes the collaborative modeling of the sparsity and local abnormality of illegal expansion patches through a sparse expression optimization strategy with local abnormality prior and spatial structure preservation constraints, implements an expression enhancement strategy for suspicious patches by introducing abnormality prior scores in the sparse expression optimization process, and constructs a spatial structure preservation term through the structural similarity of adjacent image blocks. In the optimization function, the expression sparsity and neighborhood consistency of local mutation areas are constrained at the same time, thereby avoiding misjudging vegetation interference or terrain occlusion as illegal expansion changes. Compared with traditional sparse subspace clustering, the accuracy and stability of anomaly recognition are improved, and it performs significantly in early illegal expansion areas with fuzzy boundaries or weak changes. The present invention performs graph Laplacian spectral decomposition on a multi-scale sparse coefficient matrix, and combines K-means clustering to perform ... The class implements subspace partitioning to quantitatively analyze the candidate areas extracted from the disturbance subspace based on expression sparsity, geometric shape, area, change intensity and directional consistency indicators, and constructs a regularized illegal expansion risk grading model. It can significantly identify small patches, high directional deviations or multiple disturbance features in newly added illegal areas outside the historical boundaries, and generate an automated illegal expansion risk marker layer to support regulatory authorities in rapid positioning, grading management and graphic output; the present invention extracts a composite feature vector including spectral response, normalized difference vegetation index, normalized difference water index and synthetic aperture radar scattering intensity, constructs a remote sensing feature tensor through a multi-scale window, and introduces a scale-adaptive sparse self-expression model on this basis to achieve sparse reconstruction expression of the heterogeneous features of land objects at different spatial scales. Compared with the traditional single-scale sparse representation method, it can more stably extract the feature heterogeneity of illegal expansion areas in mining scenes with many interferences and mixed land objects, and improve the sensitivity and robustness to regional change structures.

[0168] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the present invention should be included in the scope of protection of the present invention.

Claims

1. A remote sensing change monitoring method for illegal expansion of mining areas based on sparse subspace clustering, characterized by: The method includes: S1. Acquire mining area image data, pre-process the mining area image data to obtain consistent image data, and extract image feature tensors using the consistent image data; S2. Construct a self-expression model with sparse constraints based on the image feature tensor, use the self-expression model to perform subspace division, and extract the target change region set based on the subspace division result; S3. Spatially overlay and compare the target change area set with the pre-acquired legal boundary data of the mining area to generate a new change area set, and extract indicators from the new change area set to obtain a quantitative indicator set of the change area; S4. Using the set of quantitative indicators of the changed areas, determine the illegal expansion risk level of the newly added set of changed areas, and generate a change monitoring report based on the illegal expansion risk level and the consistency image data.

2. The method for monitoring illegal expansion of mining areas by remote sensing based on sparse subspace clustering according to claim 1 is characterized in that: The step of acquiring mining area image data, preprocessing the mining area image data to obtain consistent image data, and extracting image feature tensors using the consistent image data includes: S11, acquiring mining area image data, and preprocessing the mining area image data to obtain consistent image data; S12, based on the spatial overlapping areas of the consistent image data at each time point, using a sliding window to perform sliding cropping on the consistent image data in a two-dimensional spatial coordinate domain to obtain an image block set; S13, extracting image features based on the image block set, and performing cascade fusion on the extracted image features to obtain a unified feature representation vector, and stacking all the unified feature representation vectors to obtain an image feature tensor; The image features include: spectral features, vegetation index features, water body index features and synthetic aperture radar scattering features.

3. The method for monitoring illegal expansion of mining areas by remote sensing based on sparse subspace clustering according to claim 2 is characterized in that: The spectral features are obtained by extracting the average pixel value of the image block corresponding to each band at each time point and concatenating them in series in chronological order; The vegetation index feature is obtained by calculating the ratio of the difference and sum of the pixel values ​​of the near infrared band and the red light band of the image block at each time point; The water index feature is obtained by calculating the ratio of the difference and sum of the pixel values ​​of the green band and the near-infrared band of the image block at each time point; The synthetic aperture radar scattering feature is obtained by extracting the average backscattering intensity value of the image block under each polarization channel at each time point and concatenating them in series according to the time and channel sequence.

4. The method for monitoring illegal expansion of mining areas by remote sensing based on sparse subspace clustering according to claim 1, characterized in that: The self-expression model with sparse constraints is constructed based on the image feature tensor, the self-expression model is used to perform subspace division, and the target change region set is extracted according to the subspace division result, including: S21. Setting a multi-scale window set based on the image feature tensor, and constructing a multi-scale image feature tensor for each scale in the multi-scale window set; S22. Utilize the multi-scale image feature tensor to construct a sparse self-expression model under adaptive structural constraints, and obtain a multi-scale joint sparse coefficient matrix based on the sparse self-expression model; S23. Perform subspace partitioning based on the multi-scale joint sparse coefficient matrix, and extract a target change region set according to the subspace partitioning result.

5. The method for monitoring illegal expansion of mining areas by remote sensing based on sparse subspace clustering according to claim 4 is characterized in that: The subspace division based on the multi-scale joint sparse coefficient matrix and the extraction of the target change region set according to the subspace division result include: S231. Constructing a non-negative symmetric similarity matrix based on the multi-scale joint sparse coefficient matrix, and calculating the graph Laplacian matrix according to the non-negative symmetric similarity matrix; S232, performing a spectral decomposition operation on the graph Laplacian matrix to obtain eigenvectors corresponding to the first K smallest eigenvalues, and concatenating the eigenvectors to form a spectral embedding matrix; S233, performing a clustering operation on the spectral embedding matrix using K-means clustering to obtain a subspace category label set for each image block; S234. Divide all subspaces in the subspace category label set into a stable background subspace set and a perturbation subspace set according to the compactness and stability of the sparse expression in each subspace; S235 . Based on the image block index in the perturbation subspace set, extract the corresponding candidate change region set, and screen the candidate change region set in combination with the sparsity index and the spatial connectivity constraint to obtain the target change region set.

6. The method for monitoring illegal expansion of mining areas by remote sensing based on sparse subspace clustering according to claim 5 is characterized in that: The clustering operation of the spectral embedding matrix using K-means clustering is performed to obtain a subspace category label set for each image block, including: S2331, randomly selecting K data points from the spectral embedding matrix as initial cluster centers for K-means clustering according to a preset number of clusters K; S2332, using Euclidean distance to calculate the distance between each data point in the spectral embedding matrix and the initial cluster center, and assigning each data point to the cluster to which the nearest cluster center belongs; S2333, recalculating the average value of all points in each cluster to iteratively update the cluster center until the cluster center reaches a preset maximum number of iterations; S2334. Assign a corresponding subspace category label to each image block in the spectrum embedding matrix to obtain a subspace category label set for each image block.

7. The method for monitoring illegal expansion of mining areas by remote sensing based on sparse subspace clustering according to claim 5, characterized in that: Based on the image block index in the perturbation subspace set, the corresponding candidate change region set is extracted, and the candidate change region set is screened by combining the sparsity index and the spatial connectivity constraint to obtain the target change region set including: S2351. Extracting a corresponding set of candidate changed regions based on the image block indexes in the perturbation subspace set; S2352. Calculate a sparsity index for each image block in the candidate change region set, and based on a preset sparsity index threshold and spatial connectivity constraints, select image blocks that meet the conditions from the candidate change region set to obtain a target change region set.

8. The method for monitoring illegal expansion of mining areas by remote sensing based on sparse subspace clustering according to claim 1, characterized in that: The target change area set is spatially superimposed and compared with the pre-acquired legal boundary data of the mining area to generate a new change area set, and the indicators of the new change area set are extracted to obtain a quantitative indicator set of the change area, including: S31, based on the two-dimensional polygon Boolean difference operation, spatially overlay and compare the target changed area set with the pre-acquired legal boundary data of the mining area to generate a new changed area set; S32. Extract indicators from the newly added changed region set to obtain a set of quantitative indicators of the changed region; the set of quantitative indicators of the changed region includes: a geometric shape indicator, an actual area indicator, an expansion direction indicator, and a change intensity indicator.

9. The method for monitoring illegal expansion of mining areas by remote sensing based on sparse subspace clustering according to claim 8, characterized in that: The geometric shape index is obtained by taking the ratio of the actual area of ​​the newly added change area to the square of its boundary length and multiplying it by a constant factor; The actual area index is obtained by counting the number of pixels contained in the newly added changed area and multiplying it by the spatial resolution in the horizontal direction and the spatial resolution in the vertical direction of the image data; The expansion direction index is obtained by calculating the angular relationship between the central moments of the newly added change areas; The change intensity index is obtained by calculating the difference between the unified feature representation vectors of the newly added change area at two time points before and after the change, and averaging the difference values ​​within the newly added change area.

10. The method for monitoring illegal expansion of mining areas by remote sensing based on sparse subspace clustering according to claim 1, characterized in that: The method of determining the illegal expansion risk level of the newly added set of changed areas by using the set of quantitative indicators of the changed areas, and generating a change monitoring report based on the illegal expansion risk level and the consistent image data includes: S41. Defining the illegal expansion risk level based on the joint distribution characteristics of the set of quantitative indicators of the change area; S42. Determine the illegal expansion risk level of each newly changed area according to the defined illegal expansion risk level, and mark the risk level of each newly changed area using the illegal expansion risk level determination result to generate an illegal expansion risk marking layer. S43. Perform three-dimensional visualization fusion of the illegal expansion risk marker layer and the consistent remote sensing image data to generate an interactive change map, and output a change monitoring report based on the interactive change map.

Images