Mine expansion remote sensing change monitoring method based on sparse subspace clustering

Abstract

The application discloses a mining area expansion remote sensing change monitoring method based on sparse subspace clustering, and relates to the technical field of mining areas. The method comprises the following steps: obtaining mining area image data, pre-processing the mining area image data to obtain consistent image data, and extracting an image feature tensor; constructing a self-expression model with sparse constraints based on the image feature tensor, performing subspace division, and extracting a target change area set according to the subspace division result; performing spatial superposition comparison between the target change area set and pre-obtained mining area legal boundary data to generate a newly added change area set, performing index extraction, and obtaining a change area quantitative index set; determining the expansion risk level of the newly added change area set, and generating a change monitoring report according to the expansion risk level and the consistent image data. The application supports rapid positioning, grading management and map output through an expansion risk marking layer.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of mining areas, in particular to a mining area expansion remote sensing change monitoring method based on sparse subspace clustering. BACKGROUND

[0002] With the rapid development of remote sensing image acquisition and processing technology, time-series remote sensing images based on satellites or unmanned aerial vehicles are widely used in dynamic monitoring of land resources and ecological environment, and play an irreplaceable important role in the identification and evidence collection of mining area expansion. Mining area expansion refers to the mining behavior that is not within the legal boundary and does not have a legal mining permit, which is characterized by strong concealment, high suddenness and fast spatial spread, and poses a significant threat to land resources, ecological environment and order.

[0003] At present, the mainstream remote sensing change detection methods mainly include image change identification method based on threshold difference and ground feature classification method based on supervised learning. The former usually sets a fixed threshold for the spectral band difference, vegetation index change or texture index change of different time remote sensing images to realize the preliminary extraction of the change area. Although the method is simple and efficient, it is easy to produce a large number of misjudgments and omissions in the face of complex ground feature mixed scenes, and the detection result accuracy decreases significantly in the presence of seasonal vegetation interference or inconsistent imaging conditions. The latter relies on a large number of labeled samples to build a classifier and realizes the discrimination of the change area through supervised learning, but the actual mining area remote sensing monitoring often lacks samples and has high labeling cost. The change patches formed in the early expansion stage have the characteristics of small size, scattered distribution and fuzzy boundary, which easily cause insufficient model training or poor generalization ability. In addition, although deep learning models have achieved certain results in image semantic segmentation tasks, they have high computational cost when applied to high-resolution large-format images, and require strict GPU computing power and data preprocessing, which is difficult to meet the actual demand of rapid response of mining areas.

[0004] In summary, there is an urgent need for innovative technical solutions in algorithm mechanisms, feature modeling and identification logic.

[0005] However, there is no effective solution to the problems in the related art. SUMMARY

[0006] To overcome the above technical problems existing in the prior art, the present application proposes a mining area expansion remote sensing change monitoring method based on sparse subspace clustering.

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

[0008] The mining area expansion remote sensing change monitoring method based on sparse subspace clustering comprises:

[0009] S1, obtain mine area image data, and pre-process the mine area image data to obtain consistent image data, and extract an image feature tensor from the consistent image data;

[0010] S2, construct a self-expression model with sparse constraints based on the image feature tensor, divide a subspace using the self-expression model, and extract a target change region set according to the subspace division result;

[0011] S3, perform spatial superposition comparison between the target change region set and pre-obtained mine area legal boundary data to generate a new change region set, and perform index extraction on the new change region set to obtain a change region quantitative index set;

[0012] S4, determine the expansion risk level of the new change region set using the change region quantitative index set, and generate a change monitoring report according to the expansion risk level and the consistent image data.

[0013] Further, obtaining mine area image data, and pre-processing the mine area image data to obtain consistent image data, and extracting an image feature tensor from the consistent image data comprises:

[0014] S11, obtaining mine area image data, and pre-processing the mine area image data to obtain consistent image data;

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

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

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

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

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

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

[0021] The synthetic aperture radar scattering feature is obtained by extracting the average backscattering intensity value of each image block in each polarization channel at each time point, and concatenating and splicing in time and channel order.

[0022] Further, based on the image feature tensor, a self-expression model with sparse constraints is constructed, the subspace is divided by using the self-expression model, and a target change region set is extracted according to the subspace division result, including:

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

[0024] S22, a sparse self-expression model under adaptive structure constraints is constructed by using the multi-scale image feature tensor, and a multi-scale joint sparse coefficient matrix is obtained by solving the sparse self-expression model;

[0025] S23, based on the multi-scale joint sparse coefficient matrix, the subspace is divided, and a target change region set is extracted according to the subspace division result.

[0026] Further, based on the multi-scale joint sparse coefficient matrix, the subspace is divided, and a target change region set is extracted according to the subspace division result, including:

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

[0028] S232, the spectral decomposition operation is performed on the graph Laplacian matrix, the eigenvectors corresponding to the first K smallest eigenvalues are obtained, and the eigenvectors are spliced to form a spectral embedding matrix;

[0029] S233, the spectral embedding matrix is clustered by using K-means clustering to obtain a subspace class label set of each image block;

[0030] S234, according to the compactness and stability of the sparse expression in each subspace, all subspaces in the subspace class label set are divided into a stable background subspace set and a disturbance subspace set;

[0031] S235, based on the image block index in the disturbance subspace set, a 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 a target change region set.

[0032] Further, the spectral embedding matrix is clustered by using K-means clustering to obtain a subspace class label set of each image block, including:

[0033] S2331、According to the preset cluster number K, randomly select K data points from the spectral embedding matrix as the initial cluster centers of K-means clustering;

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

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

[0036] S2334, assign each image block in the spectral embedding matrix to the corresponding subspace category label to obtain the subspace category label set of each image block.

[0037] Further, based on the image block index in the perturbation subspace set, a corresponding candidate change region set is extracted, and the candidate change region set is screened in combination with the sparsity index and the spatial connectivity constraint to obtain the target change region set, including:

[0038] S2351, based on the image block index in the perturbation subspace set, a corresponding candidate change region set is extracted;

[0039] S2352, calculate the sparsity index according to each image block in the candidate change region set, and screen the image blocks that meet the conditions from the candidate change region set according to the preset sparsity index threshold and the spatial connectivity constraint to obtain the target change region set.

[0040] Further, the target change region set is spatially superimposed and compared with the pre-acquired legal boundary data of the mining area to generate a new change region set, and the new change region set is subjected to index extraction to obtain a change region quantitative index set, including:

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

[0042] S32, index extraction is performed on the new change region set to obtain a change region quantitative index set;

[0043] The change region quantitative index set includes: a geometric shape index, an actual area index, an expansion direction index, and a change intensity index.

[0044] Further, the geometric shape index is obtained by obtaining the ratio of the actual area of the new 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 change region and multiplying the spatial resolution in the horizontal direction and the spatial resolution in the vertical direction of the image data;

[0046] The expansion direction index is obtained by calculating the angle relationship between the center moments of the newly added change region.

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

[0048] Further, the expansion risk level of the set of newly added change regions is determined by using the set of change region quantitative indexes, and a change monitoring report is generated according to the expansion risk level and the consistent image data, including:

[0049] S41, based on the joint distribution characteristics of the set of change region quantitative indexes, define the expansion risk level;

[0050] S42, determine the expansion risk level of each newly added change region according to the defined expansion risk level, and label the risk level of each newly added change region using the expansion risk level determination result, to generate an expansion risk marking layer;

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

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

[0053] 1、The present application realizes the collaborative modeling of the sparsity and local abnormality of the expansion patch by using the sparse expression optimization strategy with local abnormal prior and spatial structure preservation constraint, the expression enhancement strategy is implemented on the suspicious patch by introducing the abnormal prior score in the sparse expression optimization process, and the spatial structure preservation term is constructed by the structural similarity of adjacent image blocks, the expression sparsity and neighborhood consistency of the local mutation region are constrained in the optimization function, so that the vegetation interference or terrain obstruction is avoided to be misjudged as expansion change, the accuracy and stability of abnormal identification are improved compared with the traditional sparse subspace clustering, and the early expansion region with fuzzy boundary or weak change amplitude is significantly improved.

[0054] 2, The present application can realize the significant identification of small patches, high direction deviation or multiple disturbance characteristics in the region outside the historical boundary by quantitatively analyzing the candidate region extracted in the disturbance subspace according to the expression sparsity, geometric shape, area, change intensity and direction consistency index, and constructing a regularized expansion risk grading model, and generating an automatic expansion risk marking layer, thereby supporting rapid positioning, grading management and map output.

[0055] 3, The present application extracts a composite feature vector including spectral response, normalized difference vegetation index, normalized difference water body 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, thereby realizing the sparse reconstruction expression of the heterogeneity features of ground objects at different spatial scales. BRIEF DESCRIPTION OF DRAWINGS

[0056] In order to more clearly illustrate the technical solutions in the embodiments of the present application or the prior art, the drawings needed in the embodiments will be briefly introduced as follows. Obviously, the drawings in the following description are only some embodiments of the present application, and other drawings can be obtained by those skilled in the art without creative labor on the basis of these drawings.

[0057] Fig. 1 It is a flow chart of a mine expansion remote sensing change monitoring method based on sparse subspace clustering according to an embodiment of the present application.

[0058] Fig. 2 It is a flow chart of an execution step in a mine expansion remote sensing change monitoring method based on sparse subspace clustering according to an embodiment of the present application. DETAILED DESCRIPTION

[0059] To further illustrate the embodiments, the present application provides drawings, which are part of the disclosure of the present application, mainly used to illustrate the embodiments, and can explain the operating principle of the embodiments in conjunction with the related description of the specification. Those skilled in the art should understand other possible embodiments and advantages of the present application by referring to these contents.

[0060] According to an embodiment of the present application, a mine expansion remote sensing change monitoring method based on sparse subspace clustering is provided.

[0061] The present application will be further described in conjunction with the drawings and specific embodiments, such as Figs. 1-2As shown, according to the sparse subspace clustering-based mining area expansion remote sensing change monitoring method of the embodiment of the application, the method comprises the following steps:

[0062] S1, obtaining mining area image data, and pre-processing the mining area image data to obtain consistent image data, and extracting image feature tensors from the consistent image data.

[0063] Specifically, obtaining mining area image data, and pre-processing the mining area image data to obtain consistent image data, and extracting image feature tensors from the consistent image data comprises:

[0064] S11, obtaining mining area image data, and pre-processing the mining area image data to obtain consistent image data;

[0065] S12, based on the spatial overlap region of the consistent image data at each time point, using a sliding window to slide and cut the consistent image data in the two-dimensional spatial coordinate domain to obtain an image block set;

[0066] S13, image feature extraction is performed according to the image block set, and the extracted image features are cascaded and fused to obtain a unified feature representation vector, and all unified feature representation vectors are stacked to obtain an image feature tensor;

[0067] The 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 each image block corresponding to each band at each time point, and concatenating and splicing in time sequence;

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

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

[0071] The synthetic aperture radar scattering features are obtained by extracting the average backscattering intensity value of the image block under each polarization channel at each time point, and concatenating and splicing in 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 area, and the multi-temporal remote sensing image data is subjected to radiation normalization processing, terrain correction processing and joint registration processing to obtain consistent remote sensing image data (i.e. consistent image data) with a unified coordinate reference and spectral response.

[0073] Specifically, based on the consistent 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 in 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 a fixed size sliding window pixel; an image block index matrix M is established, each row in the image block index matrix records the center two-dimensional spatial coordinate domain (x i i , y i ) of the corresponding image block B.

[0074] Specifically, a spectral feature vector of each image block in the image block set at all time points and all bands is calculated, the spectral feature vector is used to represent the spectral response trend of the image block in the multi-temporal remote sensing image, and the spectral feature vector is obtained by extracting the average pixel value of the image block at each time point and each band, and concatenating and splicing in time sequence.

[0075] Specifically, a 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 change of ground vegetation coverage in the image block area, and the normalized difference vegetation index feature vector is obtained by calculating the ratio of the difference value to the sum value of the near-infrared band pixel value and the red band pixel value of the image block at each time point, wherein the near-infrared band pixel value and the red band pixel value are both the average pixel value of the image block at the band.

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

[0077] Specifically, under the condition that the consistent remote sensing image data contains synthetic aperture radar information, a 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 mine bare rock, stockyard and waste rock area, and the synthetic aperture radar scattering feature vector is obtained by extracting the average backscattering intensity value of the image block at each time point and each polarization channel, and concatenating and splicing in time and channel sequence.

[0078] Specifically, the spectral feature vector f i spec ​, vegetation index feature vector f i ndvi , water body index feature vector f i ndwi , and radar feature vector f i sar are cascaded and fused to construct a unified feature representation vector f i .

[0079] Specifically, the unified feature representation vectors f i of all image blocks B i are stacked to form a remote sensing feature tensor F (i.e., an image feature tensor).

[0080] Specifically, the image blocks are generated by cutting the consistent remote sensing image data in the overlapping area according to a fixed size, an image block index is established to maintain the spatial correspondence, spectral features, vegetation index features, water body index features, and synthetic aperture radar scattering features are extracted for the image block set, and are fused to form a remote sensing feature tensor.

[0081] S2, constructing a self-expression model with sparse constraints based on the image feature tensor, performing subspace division using the self-expression model, and extracting a target change region set according to the subspace division result.

[0082] Specifically, constructing a self-expression model with sparse constraints based on the image feature tensor, performing subspace division using the self-expression model, and extracting a target change region set according to the subspace division result include:

[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, constructing a sparse self-expression model under adaptive structure constraints using the multi-scale image feature tensor, and solving a multi-scale joint sparse coefficient matrix according to the sparse self-expression model;

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

[0086] Specifically, performing subspace division based on the multi-scale joint sparse coefficient matrix, and extracting a target change region set according to the subspace division result include:

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

[0088] S232, perform a spectral decomposition operation on the graph Laplacian matrix to obtain eigenvectors corresponding to the first K smallest eigenvalues, and concatenate 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 class label set for each image block.

[0090] Specifically, the clustering operation on the spectral embedding matrix using K-means clustering to obtain a subspace class label set for each image block includes:

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

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

[0093] S2333, iteratively update the cluster centers by recalculating the average of all points in each cluster until the cluster centers reach a preset maximum number of iterations;

[0094] S2334, assign each image block in the spectral embedding matrix to a corresponding subspace class label to obtain a subspace class label set for each image block.

[0095] S234, divide all subspaces in the subspace class label set into a stable background subspace set and a disturbance 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 disturbance subspace set, extract a corresponding candidate change region set, and screen the candidate change region set in combination with a sparsity index and a spatial connectivity constraint to obtain a target change region set.

[0097] Specifically, based on the image block index in the disturbance subspace set, extract a corresponding candidate change region set, and screen the candidate change region set in combination with a sparsity index and a spatial connectivity constraint to obtain a target change region set includes:

[0098] S2351, based on the image block index in the disturbance subspace set, extract a corresponding candidate change region set;

[0099] S2352, calculate the sparsity index for each image block in the candidate change region set, and screen the image blocks that meet the conditions from the candidate change region set based on a preset sparsity index threshold and a spatial connectivity constraint 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, which is used to represent the low-rank association structure between image blocks.

[0101] Specifically, based on the remote sensing feature tensor F, a set of multi-scale windows S is set, and for each scale sk in the set of multi-scale windows, a corresponding multi-scale remote sensing feature tensor The multi-scale remote sensing feature tensor is obtained by aggregating the image block set at scale s k The unified feature representation vector generated below is stacked to depict the structural distribution characteristics of the mine expansion at different spatial granularities.

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

[0103]

[0104] wherein, denotes the sparse coefficient matrix at scale sk; denotes the i-th row sparse coefficient vector at scale sk, which also represents the linear expression contribution of any image block in the image block set to the remaining image blocks in the feature space; N denotes the number of image blocks, denotes the error matrix at scale sk, which also represents the residual term between the actual multi-scale remote sensing feature tensor and the sparse reconstruction, and is used to accommodate the non-structural disturbances in the mine remote sensing caused by shadows, water reflection, and terrain obstruction; denotes the L1 norm of the i-th row sparse coefficient vector at scale sk, which is used to measure the sparse expression intensity of the image block B k at scale sk; the smaller the value, the more concentrated the expression, and the more significant the change; i denotes the local anomaly prior score of the image block B k , which is obtained by calculating the difference between its NDVI, NDWI, or SAR feature and the average value of the time phase, and is used to reflect the possibility of change of the image block; η denotes an anomaly adjustment factor, which controls the sensitivity of the sparse term weighting; λ1 denotes a robust error weight factor; i denotes the local anomaly prior score of the image block B i , which is obtained by calculating the difference between its NDVI, NDWI, or SAR feature and the average value of the time phase, and is used to reflect the possibility of change of the image block; η denotes an anomaly adjustment factor, which controls the sensitivity of the sparse term weighting; λ1 denotes a robust error weight factor; i ; the greater the value, the stronger the expression activation of the suspicious area; λ1 denotes a robust error weight factor; denotes the L2,1 norm of the error matrix, which also represents the overall intensity of the error on each image block, encourages the error to be concentrated in a specific area rather than globally spread, and helps to improve robustness; λ2 denotes a structure preservation weight factor; N i denotes the local anomaly prior score of the image block B ispatially adjacent set of image blocks B i all image blocks that are adjacent or less than a set distance threshold apart; denotes the scale s k lower image block B i adjacent to image block B j adaptive spatial structure weight, also denoted as sparse coefficient vector sparsity coefficient vector structure similarity; denotes the expression coefficient vector of image block B i and image block B j Euclidean distance square between the expression coefficient vectors of image block B i and image block B i ; denotes that the sparse coefficient matrix has zero diagonal to prevent image blocks from expressing themselves.

[0105] Specifically, the first term is a weighted L1 norm term (1+η·a i )||C i ||1, the weighting mechanism can strengthen the structural activation of the expanded area in the optimization process, so that it obtains higher expression weight in the sparse subspace, thereby improving the sensitivity of the model to weak and small change areas.

[0106] The second term is the reconstruction error term of L2,1 norm to enhance the robustness to local non-structural noise (terrain shadow, atmospheric disturbance, water reflection) and reduce its interference with the global expression structure.

[0107] The third term is a graph regularization constraint term where N i denotes the spatially adjacent block set of image block B i , is the structural similarity weight of 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 blur and local jump phenomenon, and improve the spatial coherence of the division result.

[0108] The sparse coefficient matrix corresponding to each scale in the entire multi-scale setting is weighted and 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 weighted summation between the corresponding sparse coefficient value at this position at all scales and the scale fusion weight, the scale fusion weight is a positive real number, the sum of all scale fusion weights is 1, and is used to balance the contribution degree of different scales to the overall expression structure.

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

[0110] where |C fusion | represents taking the absolute value of each element in the multi-scale joint sparse coefficient matrix C fusion , ensuring that the similarity is a non-negative real number, and the similarity matrix W is used to quantify the mutual similarity of any two image blocks in the sparse expression space in the image block set; the corresponding graph Laplacian matrix L is calculated according to the similarity matrix W.

[0111] Spectral decomposition is performed on the graph Laplacian matrix L to calculate the eigenvectors corresponding to the first K smallest eigenvalues of the graph Laplacian matrix, and the spectral embedding matrix Y ∈ R N×K is formed by splicing, where K represents a preset number of subspaces, and each row Y i ∈ RK represents the position of the image block B i in the low-dimensional spectral embedding space.

[0112] K-means clustering (i.e., K-means clustering) is performed on the spectral embedding matrix Y to obtain a subspace class label set L for each image block, where l i ∈ {1, 2,..., K} represents the subspace number to which the image block B i belongs, and all subspaces in the subspace class label set are divided into a stable background subspace set and a disturbance subspace set according to the compactness and stability of the sparse expression in each subspace.

[0113] According to the image block index in the disturbance subspace set, a corresponding candidate change region set is extracted, which is used to represent the region block that presents abnormal structural changes in the remote sensing feature.

[0114] For each image block B j in the candidate change region set, a sparsity index ρ j is calculated, which is used to measure the concentration and abnormality of the image block in the structural expression:

[0115]

[0116] where ρ j j represents the jth row sparse coefficient vector in the multi-scale joint sparse coefficient matrix, ||·||1 represents the L1 norm, ||·||2 represents the L2 norm, and ∈ represents a small numerical value to prevent the denominator from being zero.

[0117] A sparsity index threshold τ ρand the spatial connectivity constraint, the image blocks satisfying the condition are screened from the candidate change region set to form a suspected change region set (i.e. the target change region set) where B j represents the jth suspected change image block obtained by sparse subspace division and sparsity screening, N sus represents the number of suspected change image blocks, wherein the sparsity index ρ j <τ ρ and the image block B j is adjacent to other low sparsity image blocks in the spatial connectivity constraint, and the adjacency satisfies the condition that the Euclidean distance is less than a threshold δ.

[0118] S3, the target change region set is spatially superimposed and compared with the pre-acquired legal boundary data of the mining area to generate a new change region set, and the new change region set is subjected to index extraction to obtain a change region quantification index set.

[0119] Specifically, the target change region set is spatially superimposed and compared with the pre-acquired legal boundary data of the mining area to generate a new change region set, and the new change region set is subjected to index extraction to obtain a change region quantification index set, which includes:

[0120] S31, based on two-dimensional polygon Boolean difference operation, the target change region set is spatially superimposed and compared with the pre-acquired legal boundary data of the mining area to generate a new change region set;

[0121] S32, the new change region set is subjected to index extraction to obtain a change region quantification index set;

[0122] The change region quantification index set includes a geometric shape index, an actual area index, an expansion direction index and a change intensity index.

[0123] Specifically, the geometric shape index is obtained by obtaining the ratio of the actual area of the new change region to the square of the boundary length and multiplying it by a constant factor;

[0124] The actual area index is obtained by counting the number of pixels contained in the new change region and multiplying it by the spatial resolution in the horizontal direction and the spatial resolution in the vertical direction of the image data;

[0125] The expansion direction index is obtained by calculating the angle relationship between the central moments of the new change region;

[0126] The change intensity index is obtained by calculating the difference between the unified feature representation vectors of the new change region at two time points before and after the change, and then averaging the difference value within the new change region.

[0127] Specifically, the legal boundary data set Blegal ; a suspected change region set R sus and mining area legal boundary data set B legal Perform spatial overlay analysis, overlay analysis uses two-dimensional polygon Boolean difference operation, extract the new area outside any legal boundary, defined as a new change region set R new , the spatial logic expression of the new change region set is:

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

[0129] Wherein, the set operator \ represents the spatial difference set operation, ∪ represents the spatial union operation, used to remove the change region within the legal range, only keep the boundary expansion part.

[0130] The mining area legal boundary data set includes: mining area planning red line data B plan , indicating the legal mining boundary approved by the corresponding department; mining permit boundary data B permit , indicating the spatial permission range of the actual mining of the mining right holder; historical mining boundary data B history , indicating the legal mining area contour confirmed in the historical remote sensing image.

[0131] Specifically, calculate the geometric shape index of each new change region in the new change region set, the geometric shape index is used to measure the contour regularity of the new change region, and the geometric shape index is obtained by obtaining the ratio of the actual area of the new change region to the square of its boundary length and multiplying a constant factor.

[0132] Specifically, calculate the actual area index of each new change region in the new change region set, the actual area index is used to measure the physical expansion scale of the new change region in the remote sensing image, and the actual area index is obtained by counting the number of pixels contained in the new change region, and multiplying the spatial resolution in the horizontal direction and the spatial resolution in the vertical direction of the remote sensing image.

[0133] Specifically, calculate the expansion direction index of each new change region in the new change region set, the expansion direction index is used to represent whether there is a significant directional extension trend of the new change region in space, and the expansion direction index is obtained by calculating the angle relationship between the central moments of the new change region.

[0134] Specifically, the change intensity indicator of each new change region in the new change region set is calculated, and the change intensity indicator is used to measure the overall change amplitude of the new change region on the remote sensing feature. The change intensity indicator is obtained by calculating the difference between the unified feature representation vectors of the region at two time points before and after the change, and averaging the difference values in the region.

[0135] Specifically, according to the joint distribution characteristics of the shape indicator φ k , the area indicator A k , the expansion direction indicator θ k and the change intensity indicator Δ k , the expansion risk level is defined as a low risk region, a medium risk region and a high risk region.

[0136] Specifically, the risk level of each new change region B k is marked as a low risk region, a medium risk region and a high risk region, and an expansion risk marking layer L risk is constructed, which is used to mark the spatial position and corresponding risk level label of each new change region.

[0137] Specifically, the low risk region, the medium risk region and the high risk region are defined as follows:

[0138] The low risk region satisfies the shape indicator φ k ≥0.75, the area indicator A k ≤300m 2 , the change intensity indicator Δ k ≤0.25·Δ max , and the expansion direction indicator The medium risk region satisfies the shape indicator 0.5≤φ k <0.75, the area indicator 300m 2 <A k ≤800mm 2 , the change intensity indicator 0.25·Δ max <Δ k ≤0.6·Δ max , or the expansion direction indicator The high risk region satisfies the shape indicator φ k <0.5, the area indicator A k >800m 2 , the change intensity indicator Δ k >0.6·Δ max , or the expansion direction indicator Where θ ref ∈[0, π) represents a historical mining area typical expansion direction reference value, which is used to judge whether the new region presents abnormal directional deviation.

[0139] Shape index φ k The threshold value is 0.75 / 0.5, and the physical meaning is that the shape index is used to measure the contour regularity of the region, φ k = 1 represents a perfect circle, and φ k → 0 represents that the shape is more complex and more “tailing” or “spreading”. The setting basis is that a region shape is compact, usually a natural object or an engineering control boundary, which is considered to be low risk when the value is higher than 0.75; the region edge is complex and has a certain extension when the value is between 0.5 and 0.75, which may exist temporary or unplanned construction; and when the value is lower than 0.5, it is linear expansion, snake shape, scattered point shape, and is considered to be high risk.

[0140] Area index A k The threshold value is 300m 2 / 800m 2 , and the physical meaning is to represent the absolute coverage area of the expanded region. The setting basis (taking the remote sensing image resolution of 0.5-2.0m as an example) is that 300m 2 ≈ 30m x 10m, which is equivalent to a newly opened small road or a simple work belt; 800m 2 ≈ the standard small open-pit mining area single block new work boundary.

[0141] Change intensity index Δ k The threshold value is 0.25 / 0.6, and the physical meaning is to measure the time sequence change degree of the region in the remote sensing feature space, and the larger the value, the more intense the change. The setting basis is that after the change value is normalized to 0-1, the quantile point is divided, that is, the first 25% is a low fluctuation region, which may be due to seasonal changes; 25-60% is a medium change region, which has risks but needs auxiliary judgment; and higher than 60% is a dramatic change, which is usually human excavation and disturbance behavior.

[0142] Expansion direction index θ k The threshold value of the deviation angle is and The physical meaning is to represent the angular deviation degree between the main expansion direction of the new region and the main direction of the historical mining area. The setting basis is that the larger the deviation between θ k and the direction θ ref of the historical mining area, the more likely the expansion behavior deviates from the approved work surface; π / 12 ≈ 15°, π / 6 ≈ 30°; according to the terrain engineering specification and the red line control line allowed “error band” of the mine planning, the human control mining error should not be greater than 15°, and more than 30° is likely to be a newly opened path or illegal earth removal.

[0143] S4, using the change region quantification index set to determine the expansion risk level of the new change region set, and generating a change monitoring report according to the expansion risk level and the consistency image data.

[0144] Specifically, the expansion risk level of the newly added change region set is determined by using the change region quantification index set, and a change monitoring report is generated according to the expansion risk level and the consistent image data, which includes:

[0145] S41, defining an expansion risk level based on the joint distribution characteristics of the change region quantification index set;

[0146] S42, determining the expansion risk level of each newly added change region according to the defined expansion risk level, and labeling the risk level of each newly added change region using the expansion risk level determination result, to generate an expansion risk marking layer;

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

[0148] Specifically, the change monitoring report includes time information, spatial location, change amplitude, risk level and suggestion.

[0149] Embodiment:

[0150] A natural resource platform receives two periods of remote sensing image data covering a certain open-pit coal mine in a certain area: the first one is an optical image collected on August 25, 2024, with a resolution of 0.8 meters, containing 4 spectral bands; the second one is a SAR radar image obtained on September 5, 2024, with a resolution of 10 meters, providing a VV polarization channel.

[0151] Without obtaining any manually labeled data, the platform uses the method of the present application to monitor and analyze whether the mine area has expansion behavior. The scene belongs to the complex situation of "multi-temporal remote sensing data + mixed ground object spectrum + few samples of sudden change".

[0152] Firstly, the remote sensing data preprocessing module is called to complete the radiation normalization and terrain correction of the images on August 25 and September 5, to ensure that the same ground object maintains consistent coordinates and spectral response in the time series images. Then fixed window sliding segmentation is performed in the overlapping area of about 12 square kilometers, generating a total of 18200 image blocks, each with a size of 256x256 pixels. The spectral features, NDVI, NDWI and SAR backscatter intensity of each image block are extracted, and they are spliced into a remote sensing feature vector with a dimension of 124 to form a remote sensing feature tensor.

[0153] Without using any training samples, a multi-scale sparse representation model is constructed, and structure preserving and anomaly excitation optimization are performed on the sparse coefficient matrix. During the sparse optimization process, 47 image blocks have abnormal scores a isignificantly higher than other regions (more than 3 times the global mean), where the sparse expression weight of the image block numbered #8821 is mainly dominated by the distant non-adjacent sub-block, and its sparse L1 norm is 0.021, L2 norm is 0.523, and sparse degree index ρ 8821 is only 0.04, which is significantly lower than the average level of 0.57 in the whole region, so it is marked as a structural abnormal significant block.

[0154] Subsequently, the joint sparse coefficient matrix is constructed into a spectral graph, the graph Laplacian matrix decomposition is performed, the spectral embedding space Y is obtained, and all image blocks are divided into 5 subspaces using K-means clustering, wherein the 4th subspace contains a large number of abnormal blocks with a sparse degree index lower than 0.2, and is distributed in the southwest direction of the main mining area of the mine, which is classified as a disturbance subspace and marked as a candidate change area.

[0155] Further, the candidate change area is spatially superimposed with the legal boundary data of the mine area, and the results show that the four change area blocks numbered #8710, #8821, #8837 and #8842 are all located outside the planning red line Bplan and the historical mining boundary B history , and the newly added change area set R new is generated accordingly.

[0156] The index calculation and risk assessment are performed on the newly added area, the area numbered #8821 has an area of 1130 square meters, a shape compactness φ = 0.36, a change intensity index Δ = 0.84·Δ max , and an expansion direction θ = 1.98 rad, with a deviation of 0.78 rad from the historical expansion direction. The area is marked as a high-risk expansion area.

[0157] Finally, the expansion risk marking layer L risk is generated, and 6 high-risk areas, 12 medium-risk areas and 21 low-risk areas are identified. The high-risk area numbered #8821 is pushed to the terminal in real time.

[0158] The traditional method A (NDVI difference, threshold set to 0.12), method B (supervised random forest based, training sample 1200) and method C (semantic segmentation model based on deep FCN, training set same as B) are used in the area at the same time, and the comparison results of the above traditional methods and the method D of the present application are shown in Table 1.

[0159] Table 1 Comparison results of the method of the present application

[0160]

[0161] The application accurately identifies all new change regions without using training samples, while the traditional method A has two false negatives, the method B has a more serious missed detection under the condition of insufficient training samples, the method C has a good performance, but its training time is 8.5 hours, the running time is 15 minutes, and the whole process analysis of the application method can be completed on an ordinary workstation in only 5 minutes, which is significantly better than the traditional method in efficiency and responsiveness.

[0162] The embodiment 1 truly shows that the application has the ability of high sensitivity, low false alarm and full-process automatic identification for mine expansion under the condition of limited resources, missing samples and complex background of change behavior, and has high engineering application value and deployment and popularization potential.

[0163] The application extracts a composite feature vector including spectral response, normalized difference vegetation index, normalized difference water index and synthetic aperture radar scattering intensity based on multi-temporal remote sensing images with a resolution better than 2 meters, constructs a remote sensing feature tensor through a multi-scale window, and introduces a scale-adaptive sparse self-expression model on this basis to realize sparse reconstruction and expression of the heterogeneity features of ground objects under different spatial scales. Compared with the traditional single-scale sparse representation method, the features of the heterogeneity of the expansion region can be more stably extracted in the mine scene with many interferences and mixed ground objects, and the sensitivity and robustness of the regional change structure are improved.

[0164] The application proposes a sparse expression optimization strategy with local anomaly priori and spatial structure constraint, realizes the collaborative modeling of the sparsity and local anomaly of the expansion patch, implements the expression enhancement strategy on the suspicious patch by introducing the anomaly priori score constructed based on the NDVI, NDWI or SAR index in the sparse expression optimization process, and constructs a spatial structure preservation term through the structural similarity of adjacent image blocks, so that the expression sparsity and neighborhood consistency of the local mutation region are simultaneously constrained in the optimization function, thereby avoiding misjudgment of the vegetation interference or terrain obstruction as expansion change, improving the accuracy and stability of anomaly identification compared with the traditional sparse subspace clustering, and performing significantly in the early expansion region with fuzzy boundary or weak change amplitude.

[0165] The application realizes the subspace division by performing graph Laplacian spectral decomposition on the multi-scale sparse coefficient matrix and combining K-means clustering, quantitatively analyzes the candidate regions extracted in the disturbance subspace according to the expression sparsity, geometric shape, area, change intensity and direction consistency index, and constructs a regularized expansion risk grading model, so that the small patch, high direction deviation or multiple disturbance features in the region newly added outside the historical boundary can be significantly identified, and an automatic expansion risk marking layer is generated to support fast positioning, grading management and map output.

[0166] In summary, by means of the technical scheme of the present application, the sparse expression optimization strategy with local anomaly prior and spatial structure constraint is used to realize the collaborative modeling of the sparse expression and local anomaly of the expanding patch, the expression enhancement strategy is implemented on the suspicious patch by introducing the abnormal prior score in the sparse expression optimization process, and the spatial structure constraint term is constructed by the structural similarity of the adjacent image blocks to constrain the expression sparsity and neighborhood consistency of the local mutation region in the optimization function, so that the vegetation interference or terrain obstruction is avoided to be misjudged as the expanding change, the accuracy and stability of the anomaly recognition are improved compared with the traditional sparse subspace clustering, and the early expanding region with fuzzy boundary or weak change amplitude is significantly improved; the present application performs graph Laplacian spectral decomposition on the multi-scale sparse coefficient matrix, and combines K-means clustering to realize subspace division, the candidate region extracted in the disturbance subspace is quantitatively analyzed according to the expression sparsity, geometric shape, area, change intensity and direction consistency index, and a regularized expansion risk grading model is constructed, the small patch, high direction deviation or multiple disturbance features in the region outside the historical boundary can be significantly recognized, and an automatic expansion risk marking layer is generated to support rapid positioning, grading management and map output; the present application 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 by a multi-scale window, and introduces a scale-adaptive sparse self-expression model on this basis to realize the sparse reconstruction expression of the feature heterogeneity of the ground object at different spatial scales. Compared with the traditional single-scale sparse representation method, the feature heterogeneity of the expanding region can be more stably extracted in the mining area scene with interference and mixed ground objects, and the sensitivity and robustness of the regional change structure are improved.

[0167] The above merely describes the preferred embodiments of the present application and should not be used to limit the present application, and any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present application should be included in the protection scope of the present application.

Claims

1. A method for mining area expansion remote sensing change monitoring based on sparse subspace clustering, characterized in that, The method comprises: S1, obtaining mine area image data, and preprocessing the mine area image data to obtain consistent image data, and extracting an image feature tensor from the consistent image data; S2, constructing a self-expression model with sparse constraints based on the image feature tensor, performing subspace division using the self-expression model, and extracting a target change region set according to the subspace division result; S2 comprises: 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, constructing a sparse self-expression model under adaptive structure constraints using the multi-scale image feature tensor, and solving a multi-scale joint sparse coefficient matrix according to the sparse self-expression model; S23, performing subspace division based on the multi-scale joint sparse coefficient matrix, and extracting a target change region set according to the subspace division result; S23 comprises: S231, constructing a non-negative symmetric similarity matrix based on the multi-scale joint sparse coefficient matrix, and calculating a graph Laplacian matrix according to the non-negative symmetric similarity matrix; S232, performing spectral decomposition operation on the graph Laplacian matrix to obtain the eigenvectors corresponding to the first K smallest eigenvalues, and splicing the eigenvectors to form a spectral embedding matrix; S233, performing clustering operation on the spectral embedding matrix using K-means clustering to obtain a subspace class label set of each image block; S234, dividing all subspaces in the subspace class label set into a stable background subspace set and a disturbance subspace set according to the compactness and stability of sparse expression in each subspace; S235, extracting a corresponding candidate change region set based on the image block index in the disturbance subspace set, and screening the candidate change region set based on the sparsity index and spatial connectivity constraint to obtain a target change region set; S3, performing spatial superposition comparison between the target change region set and the pre-acquired mine area legal boundary data to generate a new change region set, extracting an index from the new change region set, and obtaining a change region quantitative index set; S4, determining the expansion risk level of the new change region set using the change region quantitative index set, and generating a change monitoring report according to the expansion risk level and the consistent image data.

2. The sparse subspace clustering based mining area expansion remote sensing change monitoring method according to claim 1, characterized in that, The method comprises: S11, obtaining mine area image data, and preprocessing the mine area image data to obtain consistent image data; S12, using a sliding window to slide and cut the consistent image data in a two-dimensional spatial coordinate domain based on the spatial overlap region of the consistent image data at each time point to obtain an image block set; S13, extracting image features according to the image block set, and performing cascade fusion on the extracted image features to obtain a unified feature representation vector, and stacking all unified feature representation vectors to obtain an image feature tensor; The image features include spectral features, vegetation index features, water index features, and synthetic aperture radar scattering features.

3. The sparse subspace clustering based mining area expansion remote sensing change monitoring method according to claim 2, characterized in that, The spectral feature is obtained by extracting the average pixel value of each image block corresponding to each wave band at each time point, and concatenating and splicing in time sequence; The vegetation index feature is obtained by calculating the ratio of the difference value and the sum value of the pixel values of the near-infrared wave band and the red wave band of the image block at each time point; The water body index feature is obtained by calculating the ratio of the difference value and the sum value of the pixel values of the green wave band and the near-infrared wave 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 and splicing in time and channel sequence.

4. The sparse subspace clustering based mining area expansion remote sensing change monitoring method according to claim 1, characterized in that, The clustering operation on the spectral embedding matrix by using K-means clustering to obtain the subspace category label set of each image block includes: S2331, according to the preset cluster number K, randomly selecting K data points from the spectral embedding matrix as the initial clustering centers of K-means clustering; S2332, calculating the distance of each data point in the spectral embedding matrix to the initial clustering center by using Euclidean distance, and assigning each data point to the cluster to which the nearest clustering center belongs; S2333, iteratively updating the clustering centers by recalculating the average value of all points in each cluster until the clustering centers reach the preset maximum iteration number; S2334, assigning each image block in the spectral embedding matrix with a corresponding subspace category label to obtain the subspace category label set of each image block.

5. The sparse subspace clustering based mining area expansion remote sensing change monitoring method according to claim 1, characterized in that, The candidate change region set corresponding to the image block index in the perturbation subspace set is extracted, and the candidate change region set is screened in combination with the sparsity index and the spatial connectivity constraint to obtain the target change region set, which includes: S2351, based on the image block index in the perturbation subspace set, extracting the corresponding candidate change region set; S2352, calculating the sparsity index according to each image block in the candidate change region set, and screening the image blocks that meet the conditions from the candidate change region set according to the preset sparsity index threshold and the spatial connectivity constraint to obtain the target change region set.

6. The sparse subspace clustering based mining area expansion remote sensing change monitoring method according to claim 1, characterized in that, The target change region set and the pre-acquired legal boundary data of the mining area are spatially superimposed and compared to generate a new change region set, and the change region quantization index set is obtained by extracting the index of the new change region set, which includes: S31, based on the two-dimensional polygon Boolean difference operation, the target change region set is spatially superimposed and compared with the pre-acquired legal boundary data of the mining area to generate a new change region set; S32, the index of the new change region set is extracted to obtain the change region quantization index set; The change region quantization index set includes: geometric shape index, actual area index, expansion direction index and change intensity index.

7. The sparse subspace clustering based mining area expansion remote sensing change monitoring method according to claim 6, characterized in that, The geometric shape index is obtained by obtaining the ratio of the actual area of the new change region to the square of the boundary length and multiplying a constant factor; The actual area index is obtained by counting the number of pixels contained in the new change region, and multiplying the spatial resolution in the horizontal direction and the spatial resolution in the vertical direction of the image data; The expansion direction indicator is obtained by calculating the angle relationship between the center moments of the newly added change region; The change intensity indicator is obtained by calculating the difference between the unified feature representation vectors of the newly added change region at two time points before and after the change, and averaging the difference values in the newly added change region.

8. The sparse subspace clustering based mining area expansion remote sensing change monitoring method according to claim 7, characterized in that, The expansion risk level of the newly added change region set is determined by using the change region quantification indicator set, and a change monitoring report is generated according to the expansion risk level and the consistent image data, which includes: S41, defining the expansion risk level based on the joint distribution characteristics of the change region quantification indicator set; S42, determining the expansion risk level of each newly added change region according to the defined expansion risk level, and labeling the risk level of each newly added change region using the expansion risk level determination result, to generate an expansion risk marking layer; S43, three-dimensional visualization fusion of the expansion risk marking layer and the consistent remote sensing image data is performed to generate an interactive change map, and a change monitoring report is output based on the interactive change map.

Images