Variational autoencoder based on spatial domain division and spatial clustering anomaly identification method

Abstract

The application provides a spatial domain division-based variational autoencoder and spatial clustering anomaly identification method, and belongs to the field of marine mineral resources exploration and geochemical anomaly identification. Multivariate features and spatial adjacency relationship are adopted to divide local background domains, a variational autoencoder is used to learn local background distribution in each spatial domain, a Mahalanobis distance is used to comprehensively form multi-dimensional reconstruction errors to form an anomaly score, and finally, discrete anomaly points are converted into continuous anomaly areas through density clustering and boundary processing, so as to significantly improve the anomaly identification precision under complex sea area background conditions, and to realize the purposes of reducing the influence of artificial interpretation and improving the accuracy of regionalization results of ore-prospecting target area delineation.

Description

Technical Field

[0001] This application relates to the field of marine mineral resource exploration and geochemical anomaly identification, specifically proposing a variational autoencoder based on spatial domain division and spatial clustering anomaly identification method applicable to the delineation of target areas for heavy mineral exploration in marine placer deposits. Background Technology

[0002] Marine placer deposits are important coastal mineral resources. The identification of anomalous deposits and the delineation of prospecting targets for such resources typically rely on information such as geochemical element content, sediment grain size, heavy mineral identification data, and seafloor geomorphology, provenance, and hydrodynamic conditions. Due to the complex marine sedimentary environment, sample points are often discretely distributed, and the elemental and grain size background is not spatially uniform. Single thresholds or manual interpretation methods are insufficient to reliably identify weak anomalies and combined anomalies.

[0003] Existing technologies typically fall into the following categories: First, statistical thresholding methods, such as mean plus standard deviation, cumulative frequency method, exploratory data analysis, and fractal methods. These methods are usually based on single-element or multi-element combinations, using thresholds to distinguish between background and anomalies. Due to their reliance on linear relationships, distribution assumptions, or manual thresholds, they are insufficient in identifying multi-element nonlinear combination anomalies and weak anomalies. Second, spatial statistical methods, such as spatial autocorrelation analysis, geographic weighted regression, and spatial weighted principal component analysis. These methods can consider certain spatial distribution characteristics, but usually still require strong statistical assumptions or a large number of manual parameters. Third, machine learning anomaly detection methods, such as isolated forests, Gaussian mixture models, principal component analysis combined with Mahalanobis distance, autoencoders, and convolutional autoencoders. These methods have strong data-driven capabilities and can be used for multivariate anomaly detection, but they do not systematically incorporate spatial continuity and local background differences into the anomaly identification process.

[0004] The aforementioned existing technologies generally treat the entire sample area as a unified background, or divide it into zones based solely on administrative regions, geological patches, or human experience. In marine placer deposit scenarios, different sea areas are significantly affected by sediment provenance, sedimentary dynamics, paleogeography, and surface sediment types, and the same element can have different background values ​​in different regions. If global modeling is directly adopted, it is easy to miss weak mineralization in high-background areas as background, and to misjudge normal fluctuations in low-background areas as anomalies.

[0005] Moreover, existing technologies typically express anomaly identification results as discrete points or interpolated graphs. While interpolated graphs are intuitive, they are prone to smoothing effects, causing anomaly boundaries to deviate from the actual distribution of point clusters. For practical exploration work, it is necessary to further transform anomaly points into anomaly areas with clear boundaries, making it difficult to directly form mineral exploration target areas with boundaries, and thus failing to adequately meet the needs of exploration deployment.

[0006] In view of the above, this application is hereby submitted. Summary of the Invention

[0007] This application proposes a variational autoencoder based on spatial domain partitioning and a spatial clustering anomaly identification method. To address the problems of existing technologies, it employs multivariate features and spatial adjacency relationships to partition local background domains. Within each spatial domain, a variational autoencoder learns the local background distribution, and Mahalanobis distance is used to synthesize multidimensional reconstruction errors to form anomaly scores. Finally, density clustering and boundary processing are used to transform discrete anomaly points into continuous anomaly regions. This aims to significantly improve the anomaly identification accuracy under complex marine background conditions, reduce the influence of manual interpretation, and improve the accuracy of regionalized results for mineral exploration target area delineation.

[0008] To achieve the above design objectives, this application proposes a variational autoencoder based on spatial domain partitioning and a spatial clustering anomaly identification method, establishing a system architecture including the following modules:

[0009] The spatial domain partitioning module is used to decompose the entire study area into multiple local regions with relatively uniform backgrounds, so as to facilitate subsequent modeling within local subdomains.

[0010] The domain-specific VAE background learning module is used to select the VAE model with the smallest reconstruction error in each domain through a multi-model optimization strategy, and reconstruct all samples in the domain to obtain the reconstruction error vector in the standardized space.

[0011] The anomaly scoring module uses Mahalanobis distance to convert the multidimensional reconstruction error of each sample into a comprehensive score, while eliminating the influence of the scale and correlation between features, so that the degree of anomaly of different samples can be directly compared and ranked.

[0012] The anomaly zone delineation module determines anomaly points based on thresholds and transforms discrete points into anomaly zones with clear geographical boundaries through spatial clustering and boundary generation.

[0013] Based on the above system structure, firstly, geochemical data and sediment grain size data are input into the spatial domain partitioning module to divide the study area into several spatially continuous subdomains with similar backgrounds and output domain labels. Then, the data of each domain are input into the VAE background learning module within the domain, and the reconstruction error vector is output. Secondly, the error vectors of all training samples are input into the anomaly scoring module, where the error between the sample input value and the reconstructed value represents the degree of deviation from the background. The inverse covariance matrix is ​​calculated to obtain the anomaly score for each sample. Finally, based on the threshold determined by the training set, samples with anomaly scores higher than the threshold are selected and sent to the spatial clustering anomaly delineation module. After cluster analysis and boundary extraction, the anomaly polygon is output.

[0014] The variational autoencoder based on spatial domain partitioning and spatial clustering anomaly identification method introduces known mining point data, associates it with test set samples through spatial nearest neighbor matching, and performs ROC curve analysis on the anomaly scores to verify its effectiveness.

[0015] Specifically, it includes the following steps:

[0016] Step (1), spatial domain division;

[0017] Step (1.1): Modeling different geochemical backgrounds;

[0018] The study area was divided into sub-regions with relatively uniform backgrounds, and background estimation and anomaly identification were performed in each sub-region. The global Moran's I and Geary's C indices were used to test spatial autocorrelation. The LISA index was introduced to decompose the global autocorrelation into each spatial unit, revealing the spatial correlation pattern and heterogeneity characteristics of TiO2 content at the local scale.

[0019] Step (1.2): Multivariate feature clustering based on K-means;

[0020] Based on the nature of the variables, the characteristics of the variables are divided into two categories: one is compositional variables, including major element oxides and trace elements; the other is physical parameter variables, namely particle size parameters.

[0021] For component-type variables, the percentage content data are first uniformly converted, and then processed using a centered logarithmic ratio transformation.

[0022] For the granularity parameter, the following Z-score standardization is used to eliminate the dimensional differences between different variables:

[0023]

[0024] in, These are the original particle size parameter values. and These are the mean and standard deviation of the corresponding variables, respectively. This is the standardized result;

[0025] Finally, the constituent variables processed by CLR are concatenated column-wise with the standardized granular parameters to form the input feature matrix for in-domain VAE training. , where n represents the number of samples in the subdomain and m represents the dimension of the input feature;

[0026] The K-means algorithm is used to perform initial clustering of the multivariate geochemical and granular characteristics of the samples; K-means is achieved by iteratively optimizing the following objective function:

[0027]

[0028] in, Indicates the first Clusters, It is the center of the cluster;

[0029] The input features of K-means are geochemical data after centering log-ratio transformation and granular data after Z-score normalization. , where m is the feature dimension, and the feature variables include geochemical indicators related to mineralization and grain size parameters representing sedimentary characteristics;

[0030] Step (1.3): Spatial connectivity component splitting and small domain fusion;

[0031] Based on K-means clustering, samples that are spatially disconnected within the same K-means cluster are split into different domains;

[0032] First, construct a spatial adjacency graph. Where the vertex set V corresponds to all samples, and the edge set E represents the sample...

[0033] Spatial adjacency relationship between samples; the condition for adjacency between samples i and j is the Euclidean distance d between them. ij Less than the given connection radius r conn ;

[0034] Then, after obtaining the spatial adjacency graph, for each K-means class c, consider the induced subgraph. ,in, Let c be the set of samples belonging to class c. for Edges between samples; extracted using breadth-first search or depth-first search. All connected components; samples within each connected component are spatially connected and have similar geochemical composition and sediment grain size structure characteristics, and are therefore defined as an independent candidate domain, thus splitting each initial K-means class into several spatially continuous candidate domains;

[0035] Secondly, the aforementioned smaller domains are merged into neighboring larger domains, and a minimum domain size threshold N is set. min When the number of samples in the candidate domain is less than the minimum domain size threshold N min The fusion operation is triggered at this time;

[0036] After the above fusion, each domain obtained is not only spatially continuous and similar in geochemical composition and grain size structure, but also has a sufficient sample size to support subsequent statistical modeling. All domains are renumbered consecutively starting from 0 to obtain the final domain label for each sample: domain∈{0,1,2,3……,D-1}, where D is the number of final spatial domains.

[0037] Step (2), variational autoencoder background learning;

[0038] Variational autoencoders are introduced in each subdomain to learn and reconstruct background features;

[0039] A VAE model is trained in each subdomain to represent the background distribution characteristics unique to that domain. The preprocessed feature samples are input into the trained VAE to obtain the corresponding reconstruction results. The reconstructed values ​​output by the model are regarded as the background estimates of the sample under the background constraints of the subdomain, and the subsequent anomaly score is calculated accordingly.

[0040] Step (3): Calculate the overall anomaly score;

[0041] Mahalanobis distance is introduced to comprehensively measure multidimensional error information. Given a multivariate dataset, let its mean vector be... The covariance matrix is For any sample point The Mahalanobis distance from the dataset center is defined as:

[0042]

[0043] The covariance matrix of the estimation error vector when calculating Mahalanobis distance This matrix reflects the distribution structure of the training set background samples in the multidimensional error space, including the variance of each error component and the covariance between the components.

[0044] Let the training set have a total of For each sample, the set of its reconstruction error vectors is: ,in ; Calculate the mean of the error vector:

[0045]

[0046] Since VAE reconstructs background samples unbiased, It is usually very close to the zero vector; the unbiased estimator of the covariance matrix is:

[0047]

[0048] The covariance matrix is A symmetric positive definite matrix whose diagonal elements Indicates the first The variance of each error component, off-diagonal elements Indicates the first The and the first The covariance of each error component;

[0049] After obtaining the regularized covariance matrix Then, the inverse matrix is ​​calculated using a matrix inversion algorithm. The inverse matrix will be used to calculate the Mahalanobis distance for all subsequent samples. The estimation of the covariance matrix is ​​based only on the training set samples, and the test set samples are not included in the calculation.

[0050] After obtaining the inverse covariance matrix, the outlier score of each sample is defined as the Mahalanobis distance of its error vector; for any sample, the reconstructed error vector is e, and the mean of the error vector is... The anomaly score S of this sample is defined as:

[0051]

[0052] The anomaly score S measures the degree to which the sample error vector deviates from the center of the background error distribution. The larger the value, the more likely the sample is to be anomaly. For background samples, their error vector e is usually located in the dense region of the background distribution, with a smaller Mahalanobis distance and a lower score. For anomaly samples, their error vector deviates from the center of the background distribution, resulting in a significantly larger Mahalanobis distance and a higher score. After obtaining the Mahalanobis distance scores of all samples in the training set, their distribution is analyzed to determine a reasonable anomaly threshold.

[0053] Step (4): Delineation of the anomaly zone space;

[0054] Step (4.1), drawing the anomaly diagram;

[0055] Spatial interpolation is performed based on the Mahalanobis distance scores of all samples to generate a continuously distributed comprehensive anomaly score map.

[0056] Spatial interpolation methods are used to extend discrete point data to the entire study area; inverse distance radial basis functions are used to spatially interpolate the anomaly scores of discrete sampling points to generate a continuous anomaly distribution map covering the entire area.

[0057] Step (4.2): Delineation of the abnormal area;

[0058] The density-based spatial clustering algorithm DBSCAN performs cluster analysis on outliers;

[0059] An adaptive method is used to estimate the neighborhood radius ε. For each outlier point in the set, the distance between each point and its k-th nearest neighbor is calculated, resulting in a distance array d. knn Take the median d of the array. medThe base radius is then multiplied by an adjustment factor. The final neighborhood radius ε is obtained as follows:

[0060]

[0061] in, Used to control the tightness of clustering, setting The minimum number of points is set to 5, and k is set to 4 to ensure that the cluster has a certain size and to avoid pseudo-clusters composed of a few points.

[0062] Input the coordinates of outliers into the DBSCAN algorithm to obtain the cluster label for each point; points with a cluster label of -1 are considered noise and do not participate in the subsequent generation of outlier regions; through clustering, spatially clustered outliers are merged into the same cluster, while isolated points are removed.

[0063] For all points within each anomaly cluster, with its planar coordinates as the center and radius... Generate a buffer polygon; the buffer radius should be set to balance the distribution range of the point group and the mapping requirements, as shown in the following formula:

[0064]

[0065] in, The neighborhood radius of DBSCAN. Let be the buffer coefficient, and take . It is 0.35;

[0066] Geometrically merge the buffer polygons of all points within the same cluster to obtain a preliminary covering polygon;

[0067] The merged polygons are subjected to morphological closing operations, and positive buffering is performed on the polygons. Then, negative buffering is performed on the positively buffered result. The generated polygons are then simplified using Douglas-Peucker simplification.

[0068] In step (1.2), values ​​less than or equal to 0 are replaced with a very small constant; subsequently, the following CLR transformation is performed on each sample:

[0069]

[0070] in, Let be the value of the i-th component variable in the sample, and D be the number of component variables; This is the result of the variable after CLR transformation.

[0071] Step (1.3) describes the adaptive determination of the connection radius r. conn For each sample, calculate its k-th... nn Find the distance to the nearest neighbor, and then take the median d of that distance for all samples.med And multiply by a coefficient , Let be the radius of the connection, then:

[0072]

[0073] Step (1.3) includes the following fusion strategy process: Graph adjacency-based fusion: First, using the constructed spatial adjacency graph G, the domain labels of other samples in the neighborhood of each sample in the small domain are counted; then, for the small domain D... small Count the non-D occurrences in the neighborhood of all its samples small The domain labels are selected, and the domain with the highest frequency of occurrence is chosen as the target domain. small All samples are merged into the target domain; this process is repeated until no new subdomains are generated or the maximum number of iterations is reached.

[0074] Forced fusion based on spatial nearest neighbors: If isolated small regions still exist after graph adjacency fusion, then the spatial nearest neighbor method is used for forced fusion; specifically, for each sample in the small region, its k-th neighbor is searched globally. nn For each nearest neighbor, the sample's domain label is changed to the largest domain label that appears most frequently among these nearest neighbors; if no large domain is found among the nearest neighbors, the domain label with the highest frequency is used. Considering both fusion stability and spatial local representativeness, the domain label is selected as follows: The value is 20. This parameter is used to determine the local spatial neighborhood range of isolated small-domain samples, using the labels of large-domain or majority-domain samples within the neighborhood as the basis for fusion.

[0075] Step (2) includes,

[0076] Step (2.1), Intra-domain VAE training;

[0077] A single-hidden-layer symmetric VAE structure is adopted. First, the input features are extracted through a fully connected hidden layer, and then the mean vectors of the latent variables are output separately. Sum of logarithmic variance vector The hidden layer uses the ReLU activation function to extract nonlinear information from the input features; for the input sample x, the encoding process can be represented as:

[0078]

[0079] in, For the input sample vector, and These are the weight matrix and bias vector from the input layer to the hidden layer, respectively. Output for hidden layer;

[0080] Based on this, the distribution parameters of the latent variables are further output:

[0081]

[0082]

[0083] in, This represents the mean vector of the latent variables. The logarithmic form of the variance of the latent variables. As a potential spatial dimension; and These are the weight matrix and bias vector of the corresponding output layer, respectively;

[0084] The decoder structure is basically symmetrical to the encoder, also containing a hidden layer and an output layer. The hidden layer uses the ReLU activation function. Since the input features have undergone CLR transformation and standardization, their numerical range is no longer limited to a fixed interval. Therefore, a linear form for the output layer is more in line with the numerical characteristics of the feature space. The reconstruction process is as follows:

[0085]

[0086] in, For the reconstruction result, W4 and b4 are the parameters from the latent layer to the decoded hidden layer, and W3 and b3 are the parameters from the decoded hidden layer to the output layer.

[0087] The network size is dynamically determined by the dimension of the input features; let the dimension of the input features be m, then the dimension of the hidden layers is set as:

[0088]

[0089] Where h1 represents the number of neurons in the hidden layer. This means taking the larger of the two values; this setting ensures that the network has sufficient expressive power on the one hand, and avoids the network being too narrow when the feature dimension is small on the other hand.

[0090] In the latent space, sampling is performed using a reparameterization technique. The dimension of the latent space is defined as:

[0091]

[0092] Where b represents the intermediate scale obtained by compressing the input dimension. Indicates rounding down;

[0093] The potential spatial dimension is set as follows:

[0094]

[0095] in, As a potential spatial dimension, This is the latent dimension scaling factor, used to control the degree of compression of the latent space;

[0096] Step (2.2), reconstruction error analysis;

[0097] Reconstruct all samples within the domain, including the training and test sets, and calculate the reconstruction error vector;

[0098] Let the input vector of the sample in the processed feature space be... The corresponding reconstruction result is ;in, Indicates the input feature dimension. These are the actual feature values ​​of the sample. The reconstructed value given by VAE;

[0099] The reconstruction error vector of this sample is defined as:

[0100]

[0101] in, To reconstruct the error vector, Indicates the first Reconstruction bias of each feature; if This indicates that the actual value of the feature is higher than the background reconstruction value given by the model; if If the value is less than the background reconstruction value, it means that the actual value of the feature is lower than the background reconstruction value; the larger the absolute value of the error vector, the more obvious the deviation of the sample from the background pattern on that feature.

[0102] In step (4.1), the sampling point location is set as follows: ( ), corresponding to the abnormal score Then the point to be estimated Interpolation results at It can be represented as a linear combination of radial basis functions:

[0103]

[0104] in, For Euclidean distance, the basis functions are taken in inverse distance form. , It is the distance between the point to be estimated and the sampling point, and the weighting coefficient. Through precise interpolation conditions The linear equations are solved to obtain the smoothing parameter, which is set to 0 to ensure that the interpolated surface strictly passes through all known sampling points.

[0105] This application also proposes an electronic device, including a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the computer program, it implements the above-mentioned variational autoencoder based on spatial domain partitioning and spatial clustering anomaly recognition method.

[0106] This application also proposes a computer-readable storage medium storing a computer program that, when executed, can implement the above-mentioned variational autoencoder based on spatial domain partitioning and spatial clustering anomaly identification method.

[0107] In summary, this application has the following advantages and beneficial effects compared with the prior art:

[0108] 1. This application can be specifically applied to the anomaly identification of heavy minerals in marine placer deposits such as ilmenite, rutile, zircon, magnetite, and monazite. It can also be extended to the identification of multi-element geochemical anomalies in coastal areas. It can adaptively handle the spatial variation of geochemical and sedimentary units under different geological backgrounds, fully explore the synergistic anomaly information among multiple characteristic variables, delineate prospecting target areas, and provide a scientific and interpretable decision-making basis for marine placer resource exploration.

[0109] 2. This application can quickly and accurately output the boundary of the anomaly area. It outputs closed anomaly area polygons through spatial clustering, buffer merging and morphological processing, which has high efficiency and economy in the deployment of mineral exploration target areas.

[0110] 3. This application has a high adaptability to spatial heterogeneity. It obtains sub-regions with different spatial backgrounds by dividing the spatial domain, instead of using a single global background. It models separately in local spatial domains with similar backgrounds to effectively reduce the risk of false alarms in high background areas and false alarms in low background areas.

[0111] 4. This application can comprehensively and efficiently extract multivariate nonlinear anomalies. The application of variational autoencoders can learn the complex relationship between geochemical elements and grain size parameters. Compared with linear dimensionality reduction or single-element thresholding methods, it is more suitable for complex sedimentary environments.

[0112] 5. The anomaly score of this application is more stable. By adopting Mahalanobis distance to incorporate the error covariance structure into the calculation, it can simultaneously handle feature correlation and scale difference, making the comprehensive anomaly score more objective. Attached Figure Description

[0113] The implementation process of the present application will now be further explained and illustrated with reference to the following figures.

[0114] Figure 1 This is a schematic diagram of the system structure for implementing the variational autoencoder based on spatial domain partitioning and the spatial clustering anomaly identification method described in this application;

[0115] Figure 2 This is a flowchart of the variational autoencoder and spatial clustering anomaly identification method based on spatial domain partitioning described in this application;

[0116] Figure 3 This is a schematic diagram of spatial domain division;

[0117] Figure 4 This is a schematic diagram of the spatial distribution of TiO2;

[0118] Figure 5 For LISA clustering diagrams;

[0119] Figure 6 The average cost function curves corresponding to different cluster numbers k;

[0120] Figure 7 This is a schematic diagram illustrating the splitting and merging of spatial domains.

[0121] Figure 8 This is a diagram of the VAE network architecture.

[0122] Figure 9 This is a schematic diagram of the VAE learning and training process;

[0123] Figure 10 This is a scatter plot showing the relationship between the measured TiO2 content and the predicted VAE background value.

[0124] Figure 11 A histogram showing the frequency distribution of Mahalanobis distance scores;

[0125] Figure 12 A comparison chart of ROC curves; Detailed Implementation

[0126] To better understand the above-mentioned objectives, features, and advantages of this application, the application will be further described below in conjunction with the accompanying drawings and embodiments. Many specific details are set forth in the following description to provide a thorough understanding of this application; however, this application may be implemented in other ways than those described herein, and therefore, this application is not limited to the specific embodiments disclosed below.

[0127] Example 1, such as Figure 1 As shown, this application proposes a variational autoencoder and spatial clustering anomaly identification method based on spatial domain partitioning. This method integrates spatial domain partitioning, variational autoencoder background learning, Mahalanobis distance anomaly scoring, and DBSCAN spatial clustering to achieve the identification and delineation of mineralization anomalies under complex background conditions. This provides methodological support for predicting target areas for heavy mineral exploration in marine placer deposits. To this end, a system structure including the following core modules is established:

[0128] The spatial domain partitioning module is used to decompose the entire study area into multiple local regions with relatively uniform backgrounds to facilitate subsequent modeling within local subdomains. This includes: addressing the spatial heterogeneity problem by performing initial K-means clustering based on geochemical multivariate characteristics and sediment grain size structure characteristics; then, combining spatial adjacency relationships to split spatially disconnected samples in the same K-means class into different connected components, and merging components that are too small in scale, ultimately obtaining several spatially continuous subdomains with similar backgrounds.

[0129] The domain-specific VAE background learning module is used to select the VAE model (Variational Autoencoder, or VAE for short, is an unsupervised deep learning model consisting of an encoder, latent variable sampling, and a decoder, used to learn the background distribution of samples and output the reconstruction result) with the minimum reconstruction error in each domain through a multi-model optimization strategy. This VAE is then used to reconstruct all samples within the domain, obtaining a reconstruction error vector in the standardized space. This vector reflects the degree to which each element of the sample deviates from the background pattern of the domain. Specifically, the module involves using VAEs to learn stable background patterns within each subdomain and then using the reconstruction error as a quantitative indicator of the degree of deviation. In each subdomain, a set of VAEs with identical network structures and parameters are trained independently. The VAE maps the input to the distribution parameters of latent variables through the encoder and then reconstructs the input through the decoder. Its loss function consists of the reconstruction error and the KL divergence.

[0130] The anomaly scoring module uses Mahalanobis distance to convert the multidimensional reconstruction error of each sample into a comprehensive score, while eliminating the influence of dimensions and correlations between features, allowing for direct comparison and ranking of the anomaly levels of different samples. Specifically, it aggregates the reconstruction error vectors of all samples, calculates the covariance matrix, adds a regularization term, and then inverts it to obtain a weight matrix representing the correlation between features. For any sample, its Mahalanobis distance is calculated using its reconstruction error vector and weight matrix and used as the anomaly score. Using Mahalanobis distance automatically eliminates the interference caused by dimensions and correlations between features, essentially compressing the multidimensional error into a value that can be compared horizontally. The higher the score, the further the error pattern of this sample deviates from the normal background, and the more likely it is to be an anomaly sample.

[0131] The anomaly zone delineation module identifies anomalies based on a threshold. Through spatial clustering and boundary generation, discrete points are merged into anomaly zones with clear geographical boundaries. Specifically, an anomaly threshold is determined based on the high quantile of the Mahalanobis distance in the training set (e.g., 95%), and samples with scores higher than the threshold are marked as anomalies. DBSCAN spatial clustering is performed on the anomalies to automatically identify spatially clustered anomalous points. Buffer analysis, merging, and morphological closing operations are performed on the points within each cluster to generate polygons that fit the shape of the point clusters. These polygons are then transformed and mapped to a spatial geographic coordinate system. Anomaly boundaries are superimposed on the anomaly map obtained through inverse distance radial basis function interpolation, ultimately resulting in anomaly zones with clear boundaries.

[0132] like Figures 2 to 12 As shown, based on the aforementioned system structure including core modules, the variational autoencoder and spatial clustering anomaly identification method based on spatial domain partitioning described in this application firstly inputs geochemical data and sediment grain size data into the spatial domain partitioning module, dividing the study area into several spatially continuous subdomains with similar backgrounds and outputting domain labels; then, the data of each domain is input into the VAE background learning module within the domain, outputting a reconstruction error vector; secondly, the error vectors of all training samples are input into the anomaly scoring module, using the error between the sample input value and the reconstructed value to represent the degree of deviation from the background, calculating the inverse covariance matrix, and obtaining the anomaly score for each sample; finally, based on the threshold determined by the training set, samples with anomaly scores higher than the threshold are selected from all samples and sent to the spatial clustering anomaly delineation module, outputting anomaly polygons after clustering analysis and boundary extraction.

[0133] Furthermore, this application can incorporate known mining site data, associate it with test set samples through spatial nearest neighbor matching, and perform ROC curve analysis on abnormal scores to verify its effectiveness.

[0134] Specifically, the method includes the following implementation steps:

[0135] Step (1), spatial domain division;

[0136] Spatial heterogeneity of geochemical background refers to the fact that the background element content is not uniformly distributed in space, but rather exhibits systematic differences with location. The causes of this phenomenon are usually closely related to regional geological backgrounds, including lithological differences, tectonic unit differentiation, surface geochemical processes, and anthropogenic disturbances. For example, in areas with basic rocks, the background Ti value is usually higher than in areas with acidic rocks; in alluvial plains and areas covered by residual and colluvial deposits, element content also varies due to differences in material sources and sedimentary sorting processes.

[0137] Step (1.1): Modeling different geochemical backgrounds;

[0138] This embodiment divides the study area into sub-regions with relatively uniform backgrounds, and performs background estimation and anomaly identification within each sub-region to effectively segment the area for different background patterns; for example... Figure 3 As shown;

[0139] To investigate the spatial distribution characteristics of TiO2 content, spatial autocorrelation tests were conducted using the global Moran's I and Geary's C indices. Considering the different focuses of Moran's I and Geary's C indices in judging spatial correlation, Moran's I primarily reflects the overall linear correlation between attribute values ​​in adjacent regions, while Geary's C focuses more on local differences between adjacent locations and is more sensitive to local variations. Using both in combination allows for mutual verification and avoids the bias that may arise from using a single indicator. Figure 4 The diagram shows the spatial distribution of TiO2 in a certain region of China. Moran's I value is 0.454 and Geary's C value is 0.560, both indicating significant positive spatial autocorrelation of TiO2 content. Furthermore, both p-values ​​are less than or equal to 0.001, indicating statistical significance. This suggests that high-value areas are adjacent to each other, and low-value areas are adjacent to each other, rather than being randomly distributed. It also indicates that the actual numerical differences between adjacent samples are small, and the spatial structure is relatively stable, suggesting that TiO2 content has a clear spatial aggregation structure. The spatial distribution of TiO2 is as follows: Figure 4 As shown, the existence of this spatial structure provides a statistical basis for using spatial domain division methods to address the differences in background values ​​across different geological units and sedimentary environments. However, global indicators can only determine whether there is an overall aggregation trend; they cannot reveal the specific spatial location of aggregation areas, nor can they identify whether there are outliers in local areas that contradict the overall trend.

[0140] To this end, this application further introduces the LISA (Local Spatial Autocorrelation) index, decomposing the global autocorrelation to each spatial unit, revealing the spatial correlation patterns and heterogeneity characteristics of TiO2 content at the local scale; such as Figure 5As shown in the LISA cluster diagram, the spatial distribution of TiO2 content exhibits a clear structural characteristic. The HH (high-high) cluster is mainly concentrated in the central-eastern part, exhibiting a blocky distribution and forming a continuous high-value zone. This area corresponds to the distribution of sandy sedimentary belts, areas with strong tidal alteration, and some ancient sedimentary bodies. The nearshore provenance is influenced by Precambrian metamorphic rocks and Mesozoic intrusive rocks, providing a source of stable heavy minerals such as ilmenite and zircon. At the same time, other nearshore areas in the study area have developed sandy sediments and strong hydrodynamic conditions, which are conducive to the transport or washing of fine-grained and light mineral components, resulting in the relative enrichment of high-density and stable titanium-bearing heavy minerals. Therefore, the HH area can be interpreted as a spatial unit with relatively favorable conditions for the enrichment of titanium-bearing heavy minerals. The LL (low-low) cluster is mainly distributed in the north, forming a clear spatial differentiation from the HH area. Considering the distribution of surface sediments and provenance conditions in the study area, this type of area is mostly associated with silty and muddy sediments or weak hydrodynamic environments. The accumulation of fine-grained sediments dilutes the content of heavy minerals, and the weak dynamic environment is unfavorable for the mechanical sorting and relative enrichment of heavy minerals, resulting in an overall low TiO2 content. The HH and LL zones each occupy relatively continuous spatial ranges, indicating significant differences in TiO2 content across different geological units or sedimentary environments, with relatively good internal homogeneity. Furthermore, a small number of HL (high-low) and LH (low-high) anomalies are present in the figure, reflecting certain local anomalies within the study area. These distribution characteristics indicate that TiO2 content is generally clustered spatially, but there are differences between different regions, with localized anomalies. These results validate the conclusions of global autocorrelation, clarify the spatial location and extent of the clusters, and provide spatial basis for subsequent spatial domain division, regional background modeling, and anomaly identification.

[0141] Step (1.2): Multivariate feature clustering based on K-means;

[0142] To ensure that different types of features can be effectively trained in the same model, variables are classified.

[0143] Specifically, based on the nature of the variables, the characteristics of the variables are divided into two categories: one is compositional variables, including major element oxides and trace elements; the other is physical parameter variables, namely particle size parameters.

[0144] For component variables, considering their closure effect and inconsistent dimensions, the percentage data are first uniformly converted, and then the centered logarithmic ratio transformation (CLR) is used for processing.

[0145] To avoid zero or negative values ​​in logarithmic operations, values ​​less than or equal to 0 are replaced with a very small constant. Subsequently, each sample undergoes the following CLR (Centered Logarithmic Ratio Transform, a logarithmic ratio transformation method used for component data processing to reduce the influence of summation constraints on elemental percentages or component contents) transformation:

[0146]

[0147] in, Let be the value of the i-th component variable in the sample, and D be the number of component variables. This is the result of the variable after CLR transformation;

[0148] The essence of the above CLR transformation is to compare each component variable with the geometric mean of all component variables within the sample, thereby mapping the component data that was originally constrained by the total amount to the real number space.

[0149] For the granularity parameter, the following Z-score standardization is used to eliminate the dimensional differences between different variables:

[0150]

[0151] in, These are the original particle size parameter values. and These are the mean and standard deviation of the corresponding variables, respectively. This is the standardized result;

[0152] Finally, the constituent variables processed by CLR are concatenated column-wise with the standardized granular parameters to form the input feature matrix for in-domain VAE training. , where n represents the number of samples in the subdomain and m represents the dimension of the input feature;

[0153] To initially define the spatial background, this embodiment uses the K-means algorithm to perform initial clustering of the samples based on their multivariate geochemical and granular characteristics. K-means is an unsupervised clustering algorithm whose goal is to divide n samples into k clusters such that the sum of the squared Euclidean distances from each sample to the center of its cluster is minimized. Mathematically, K-means achieves this by iteratively optimizing the following objective function:

[0154]

[0155] in, Indicates the first Clusters, It is the center of the cluster;

[0156] The input features of K-means are geochemical data after centering log-ratio transformation and granular data after Z-score normalization. Where m is the feature dimension, and the feature variables include geochemical indicators related to mineralization and grain size parameters representing sedimentary characteristics; number of clusters The choice of k needs to balance geological significance and computational efficiency. Based on the geological complexity of the study area and referring to the elbow rule, k can be chosen as 6. Figure 6 As shown;

[0157] As the number of clusters, k, gradually increases from 1, the average loss function value shows a decreasing trend. When k increases from 1 to 6, the decrease in the average loss function value is more significant; when k exceeds 6, the decrease in the average loss function value tends to level off, and the curve shows a clear inflection point. Based on this, the optimal number of clusters, k, is determined to be 6. The initial clustering results assign an initial category label to each sample. This label reflects the background similarity of samples at a global scale;

[0158] Step (1.3): Spatial connectivity component splitting and small domain fusion;

[0159] K-means clustering only considers the similarity of geochemical composition and grain size characteristics, but does not consider the spatial location of the samples. In actual geological environments, samples of the same lithology may be separated by tectonic fractures, resulting in samples of the same category being distributed in multiple spatially unconnected regions. These regions should be considered as independent background domains. Therefore, based on K-means clustering, it is necessary to further split spatially unconnected samples within the same K-means cluster into different domains according to their spatial proximity. The entire splitting and merging process is as follows: Figure 7 As shown;

[0160] Therefore, it is necessary to construct a spatial adjacency graph. Here, the vertex set V corresponds to all samples, and the edge set E represents the spatial adjacency relationship between samples; the condition for samples i and j to be adjacent is the Euclidean distance d between them. ij Less than the given connection radius r conn To adaptively determine r conn The following method can be used: For each sample, calculate its value up to the k-th digit. nn Find the distance to the nearest neighbor, and then take the median d of that distance for all samples. med And multiply by a coefficient In the experiment, take Let 1.5 be the radius of connection, and take k. nn If the value is 8, then:

[0161]

[0162] The adaptive method described above can automatically adjust the connection radius according to the density of the sample distribution, avoiding unreasonable connections or breaks caused by a globally fixed radius.

[0163] Then, after obtaining the spatial adjacency graph, for each K-means class c, consider the induced subgraph. ,in, Let c be the set of samples belonging to class c. for Edges between samples; extracted using breadth-first search or depth-first search. All connected components; samples within each connected component are spatially connected and have similar geochemical composition and sediment grain size structure characteristics. Therefore, they are defined as an independent candidate domain, thereby splitting each initial K-means class into several spatially continuous candidate domains.

[0164] After component partitioning, some candidate domains may be generated that are extremely small. These domains, due to their limited sample size, cannot support stable training of the subsequent variational autoencoder and lack statistical representativeness. These small domains need to be merged into neighboring larger domains; therefore, a minimum domain size threshold N can be set. min When the number of samples in the candidate domain is less than the minimum domain size threshold N min When this occurs, a fusion operation will be triggered, including the following fusion strategy process:

[0165] Graph adjacency-based fusion: First, using the constructed spatial adjacency graph G, we statistically analyze the domain labels of other samples within the neighborhood of each sample in a small domain; then, for the small domain D... small Count the non-D occurrences in the neighborhood of all its samples small The domain labels are selected, and the domain with the highest frequency of occurrence is chosen as the target domain. small All samples are merged into the target domain; this process is repeated until no new subdomains are generated or the maximum number of iterations is reached.

[0166] Forced fusion based on spatial nearest neighbors: If isolated small regions still exist after graph adjacency fusion, then the spatial nearest neighbor method is used for forced fusion; specifically, for each sample in the small region, its k-th neighbor is searched globally. nn For each nearest neighbor, the sample's domain label is changed to the largest domain label that appears most frequently among these nearest neighbors; if no large domain is found among the nearest neighbors, the domain label with the highest frequency is used. Considering both fusion stability and spatial local representativeness, the domain label is selected as follows: The value is 20. This parameter is used to determine the local spatial neighborhood range of isolated small-domain samples, and the labels of large domains or most domains within the neighborhood are used as the basis for fusion. If the value is too small, the fusion result may be affected by local noise and become unstable. If the value is too large, it may weaken the spatial locality.

[0167] After the above two steps of fusion, each domain obtained is not only spatially continuous and similar in geochemical composition and grain size structure, but also has a sufficient sample size to support subsequent statistical modeling. All domains are renumbered consecutively starting from 0 to obtain the final domain label for each sample: domain∈{0,1,2,3……,D-1}, where D is the number of final spatial domains.

[0168] Step (1.4), Analysis of Spatial Domain Division Results:

[0169] Based on the above spatial domain division method, the study area is divided into 5 spatially continuous subdomains with relatively uniform geochemical and sediment grain size characteristics, numbered from Domain 0 to Domain 4, as shown in Table 1 below.

[0170] Table 1. Statistical characteristics of TiO2 content in each spatial subdomain

[0171] 0 30 0.421 0.071 0.290 0.600 1 137 0.638 0.079 0.250 0.720 2 573 0.660 0.126 0.030 2.150 3 57 0.497 0.125 0.300 0.830 4 342 0.533 0.233 0.000 0.770

[0172] Therefore, Domain 2 not only has a high mean but also exhibits significant outliers, indicating the potential development of strong local anomalies within this domain. Domain 4 has a relatively wide box, suggesting a more complex internal background. Domains 0 and 1 have narrower boxes, indicating a more concentrated distribution of samples within their domains. Overall, the spatial domain division effectively differentiates the various geochemical and sediment grain size background types within the study area, providing a reasonable spatial framework for subsequent domain background learning and anomaly identification.

[0173] Step (2), variational autoencoder background learning;

[0174] Variational autoencoders (VAEs) are introduced in each subdomain to learn and reconstruct background features. VAEs are deep generative models that can better characterize the high-dimensional nonlinear relationships in geochemical and granular feature data. By regularizing the latent variable space, the robustness and reconstruction ability of the model are improved. VAEs can not only reconstruct input samples, but also learn the probability distribution characteristics of samples in the latent space, thus providing richer distribution information and statistical basis for subsequent multivariate anomaly identification.

[0175] Since samples within each subdomain have relatively consistent background conditions, a VAE model is trained separately within each subdomain to represent the background distribution characteristics unique to that domain. The preprocessed feature samples are then input into the trained VAE to obtain the corresponding reconstruction results. For background samples, the model typically performs well in reconstruction; however, for samples deviating from the domain's background pattern, the reconstruction error increases relatively. Therefore, the reconstructed value output by the model can be considered as the background estimate of the sample under the subdomain's background constraints, and this is used to calculate subsequent anomaly scores. Specifically,

[0176] Step (2.1), Intra-domain VAE training;

[0177] Employing a single-hidden-layer symmetric VAE structure, such as Figure 8 First, the input features are extracted through a fully connected hidden layer, and then the mean vectors of the latent variables are output. Sum of logarithmic variance vector The hidden layer uses the ReLU activation function to extract nonlinear information from the input features; for the input sample x, the encoding process can be represented as:

[0178]

[0179] in, For the input sample vector, and These are the weight matrix and bias vector from the input layer to the hidden layer, respectively. Output for hidden layer;

[0180] Based on this, the distribution parameters of the latent variables are further output:

[0181]

[0182]

[0183] in, This represents the mean vector of the latent variables. The logarithmic form of the variance of the latent variables. As a potential spatial dimension; and These are the weight matrix and bias vector of the corresponding output layer, respectively;

[0184] The decoder structure is basically symmetrical to the encoder, also containing a hidden layer and an output layer. The hidden layer uses the ReLU activation function. Since the input features have undergone CLR transformation and standardization, their numerical range is no longer limited to a fixed interval. Therefore, a linear form for the output layer is more in line with the numerical characteristics of the feature space. The reconstruction process is as follows:

[0185]

[0186] in, For the reconstruction result, W4 and b4 are the parameters from the latent layer to the decoded hidden layer, and W3 and b3 are the parameters from the decoded hidden layer to the output layer.

[0187] The network size is dynamically determined by the dimension of the input features. Let the dimension of the input features be m, then the dimension of the hidden layers is set as follows:

[0188]

[0189] Where h1 represents the number of neurons in the hidden layer. This means taking the larger of the two values; this setting ensures that the network has sufficient expressive power on the one hand, and avoids the network being too narrow when the feature dimension is small on the other hand.

[0190] In the latent space, sampling is performed using a reparameterization technique. The dimension of the latent space is defined as:

[0191]

[0192] Where b represents the intermediate scale obtained by compressing the input dimension. Indicates rounding down;

[0193] The potential spatial dimension is set as follows:

[0194]

[0195] in, As a potential spatial dimension, This is the latent dimension scaling factor, used to control the degree of compression of the latent space. The larger the factor, the more information is retained in the latent space; the smaller the factor, the stronger the compression.

[0196] In the learning process of a variational autoencoder (VAE), the input data is first mapped to the latent space by the encoder, and the mean and standard deviation of the latent variables are output. Then, through reparameterization techniques, latent variables are sampled from this distribution. This ensures the sampling process is differentiable. Latent variables. The data is passed to the decoder, which maps it to the original data space to obtain the reconstruction result. The VAE optimizes the model by calculating the reconstruction loss and KL divergence. The reconstruction loss measures the difference between the original data and the reconstructed data, and the KL divergence measures the difference between the latent distribution and the prior distribution. The total loss is the sum of the two. The gradient descent method is used to update the weights and biases to optimize the model parameters.

[0197] like Figure 9 The VAE learning and training process shown uses the Adam optimizer for model training, with a learning rate set to α = 1 × 10⁻⁶. −4α represents the parameter update step size, used to control the magnitude of weight adjustment in each iteration; the maximum number of iterations is set to 2000 to limit the number of training rounds for a single model. Samples are partitioned into spatial subdomains, with 20% of the samples randomly selected from each subdomain as the test set and the remaining samples as the training set. This approach maintains consistency in the sample proportions within each subdomain and facilitates subsequent validation of the model's performance on different background units.

[0198] To reduce the impact of random initialization and optimization processes on training results, five VAE models can be independently trained within each subdomain using the same network structure and training parameters. The model with the lowest final loss is then selected as the optimal model for that subdomain. This strategy can reduce the randomness of single training results and improve the stability of background learning results within the domain, while considering computational costs.

[0199] Step (2.2), reconstruction error analysis;

[0200] Reconstruct all samples within the domain, including the training and test sets, and calculate the reconstruction error vector, which reflects on which elements, to what extent, and in which direction the sample deviates from its normal background pattern within the domain.

[0201] Let the input vector of the sample in the processed feature space be... The corresponding reconstruction result is ;in, Indicates the input feature dimension. These are the actual feature values ​​of the sample. The reconstructed value given by VAE;

[0202] The reconstruction error vector of this sample is defined as:

[0203]

[0204] in, To reconstruct the error vector, Indicates the first Reconstruction bias of each feature. If This indicates that the actual value of the feature is higher than the background reconstruction value given by the model; if This indicates that the actual value of the feature is lower than the reconstructed value of the background. The larger the absolute value of the error vector, the more significantly the sample deviates from the background pattern on that feature.

[0205] In the specific calculations, the error matrices of the training set and the test set are denoted as follows:

[0206]

[0207]

[0208] in, and Let represent the input feature matrices of the training set and the test set, respectively. and These represent the reconstruction matrices of the training and test sets, respectively. and These represent the reconstruction error matrices of the training set and the test set, respectively; each row in the matrix corresponds to the error vector of a sample, and each column corresponds to an input feature;

[0209] Reconstruction error can reflect, to some extent, how much a sample deviates from the background distribution. For samples that conform to the background pattern within the domain, the model can usually provide relatively accurate reconstruction results, with a small overall error vector. However, for samples that deviate from the background pattern, the model's reconstruction ability will weaken accordingly, leading to an increase in some components of the error vector. When multiple feature components deviate significantly at the same time, the reconstruction error can reflect the comprehensive anomaly characteristics under the combined influence of geochemistry and sediment grain size information.

[0210] Generally, the error vectors corresponding to background samples exhibit relatively stable statistical characteristics: their mean is generally close to zero, indicating that the model does not have obvious systematic reconstruction bias; their variance is relatively small, indicating that most background samples are within the normal fluctuation range that the model can effectively represent; in addition, there is usually some correlation between the error components, reflecting that the coupling relationship between input features has not completely disappeared during the reconstruction process. Abnormal samples, on the other hand, often show significant shifts in local error components or an overall error structure that deviates from the background distribution corresponding to the training samples. The subsequent calculation of abnormal scores is based on the above-mentioned differences in error distribution.

[0211] After calculating the error vectors, the error vector set is stored, and these vectors are saved in matrix form, where rows correspond to samples and columns correspond to features. After obtaining the reconstruction error, this embodiment first utilizes the training set error matrix... The covariance structure of the background error distribution is estimated, and the distribution characteristics of the outlier scores in the training set are predicted accordingly. Then, under the same background constraints, a comprehensive outlier score is calculated for all samples. This allows the normal error range expressed by the training set to be used as a reference to measure the degree to which other samples deviate from the background pattern, facilitating subsequent performance evaluation of the model. Figure 10The correlation between the measured TiO2 content of the samples and the VAE background prediction values ​​is shown. Overall, most sample points are distributed around the 1:1 line, indicating that the model's reconstruction of the background samples is relatively close to the measured values ​​and has good background representation ability. In contrast, the abnormal sample points deviate from the 1:1 line more significantly, reflecting the difference between them and the background distribution. This shows that background learning based on in-domain VAE can provide a relatively reliable foundation for subsequent anomaly score calculation and anomaly sample identification.

[0212] Step (3): Calculate the overall anomaly score;

[0213] After completing the intra-domain reconstruction, the reconstruction error vector of each sample in the feature space can be obtained. Since the fluctuation range of different features is not consistent, and there is a certain correlation between the error components, it is difficult to directly measure the degree of sample anomaly by a single component or simple summation.

[0214] Therefore, Mahalanobis distance is introduced to comprehensively measure multidimensional error information, and based on this, outlier identification, spatial clustering, and outlier region delineation are further carried out. Mahalanobis distance, first proposed by Mahalanobis in 1936, is a commonly used distance metric in multivariate statistical analysis. Unlike Euclidean distance, which assumes that variables are independent and have consistent variances, Mahalanobis distance incorporates the covariance structure into the calculation process, simultaneously considering the correlation between variables and the scale differences between different variables. Based on this characteristic, it can more reasonably reflect the degree of deviation of sample points from the center of the population distribution, and has strong applicability in multivariate anomaly identification.

[0215] Specifically, given a multivariate dataset, let its mean vector be... The covariance matrix is For any sample point The Mahalanobis distance from the dataset center is defined as:

[0216]

[0217] The geometric significance of Mahalanobis distance lies in its ability to pass through... -1 A linear transformation is applied to the original space, rotating and scaling the correlated original coordinate axes into an isotropic standard space. Under the Mahalanobis distance metric, ellipses extending along the correlated directions in the original space are transformed into circles, thus making the distance calculation unaffected by variable scale and correlation. Mahalanobis distance automatically decorrelates correlated variables using the inverse of the covariance matrix. When two variables are highly correlated, their joint fluctuations in the original space are considered normal, while significant deviations of a single variable are amplified.

[0218] To calculate the Mahalanobis distance, it is necessary to estimate the covariance matrix of the error vector. This matrix reflects the distribution structure of the training set background samples in the multidimensional error space, including the variance of each error component and the covariance between the components.

[0219] Let the training set have a total of For each sample, the set of its reconstruction error vectors is: ,in ; Calculate the mean of the error vector:

[0220]

[0221] Since VAE reconstructs background samples unbiased, It is usually very close to the zero vector, but it is still retained in actual calculations to eliminate possible systematic biases.

[0222] The unbiased estimator of the covariance matrix is:

[0223]

[0224] The covariance matrix is A symmetric positive definite matrix whose diagonal elements Indicates the first The variance of each error component, off-diagonal elements Indicates the first The and the first The covariance of each error component;

[0225] In practical applications, if the number of training samples is limited or there is a high correlation between features, the covariance matrix may be close to singular, leading to an inverse matrix. Unstable computation is not possible. To address this issue, this embodiment introduces a regularization term to slightly perturb the covariance matrix:

[0226]

[0227] in, It is the identity matrix. Let be the regularization coefficient, and take . 10 -6 This is equivalent to adding a very small variance to all features, ensuring the positive definiteness of the matrix, while having almost no impact on the structure of the original covariance matrix;

[0228] After obtaining the regularized covariance matrix Then, the inverse matrix is ​​calculated using a matrix inversion algorithm. The inverse matrix will be used to calculate the Mahalanobis distance for all subsequent samples. The estimation of the covariance matrix is ​​based only on the training set samples, and the test set samples are not included in the calculation to ensure the unbiasedness of outlier scores.

[0229] After obtaining the inverse covariance matrix, the outlier score of each sample is defined as the Mahalanobis distance of its error vector; for any sample, the reconstructed error vector is e, and the mean of the error vector is... The anomaly score S of this sample is defined as:

[0230]

[0231] The anomaly score S measures the degree to which the sample error vector deviates from the center of the background error distribution; the larger the value, the more likely the sample is to be anomaly. For background samples, their error vector e is usually located in a densely distributed area of ​​the background, resulting in a smaller Mahalanobis distance and a lower score. For anomaly samples, their error vector deviates from the center of the background distribution, leading to a significantly larger Mahalanobis distance and a higher score. After obtaining the Mahalanobis distance scores of all samples in the training set, their distribution is analyzed to determine a reasonable anomaly threshold and provide a basis for subsequent anomaly identification.

[0232] like Figure 11 As shown, the frequency distribution characteristics of the Mahalanobis distance scores of the training set samples are illustrated. Overall, the sample scores are mainly concentrated in the lower value range, exhibiting a right-skewed distribution, indicating that most training samples deviate little from the center of the background distribution. Only a few samples are located in the high value range, with a small number of samples scoring above 5.91. This result indicates that the overall reconstruction error of the background samples is low, and the anomaly scores constructed based on Mahalanobis distance can effectively reflect the degree to which samples deviate from the background distribution, thus distinguishing between background samples and potentially anomalous samples.

[0233] The threshold is determined based on the sample frequency distribution and scores, with a threshold of 95%. An anomaly labeling for a point in a spatial location is determined by its Mahalanobis distance score. A point is considered an anomaly when its score is greater than or equal to the 95th percentile of the Mahalanobis distance in the training set; as shown in the following formula:

[0234]

[0235] in, Indicates position Mahalanobis distance score at T; 95 The training set uses the outlier scores corresponding to the 95th percentile of Mahalanobis distance. Under the assumption of a multivariate normal distribution, the squared Mahalanobis distance follows a chi-square distribution, and the significance level α corresponding to the 95th percentile is 0.05, which has a certain statistical basis. Furthermore, in subsequent spatial matching verification with known mining sites, the 95th percentile threshold achieved reasonable recall and precision on multiple algorithms.

[0236] Furthermore, this embodiment checks for missing values, outliers, and spatial outliers during the data preprocessing stage, but does not simply discard high-value samples. For example, for high-value data related to the enrichment of titanium-bearing heavy minerals such as TiO2, if their coordinates, test results, and spatial distribution with adjacent stations are reasonable, they are retained as potential mineralization anomaly information for modeling. Outliers that are clearly caused by input errors, unit errors, or missing key fields are corrected or removed. During the modeling stage, a comprehensive anomaly score is calculated using VAE reconstruction error and Mahalanobis distance, and the anomaly threshold is determined using the training set score quantile, thereby identifying high-anomaly samples. This processing method avoids the accidental deletion of high-value points with mineral exploration significance and reduces the impact of erroneous data on the model results.

[0237] Step (4): Delineation of the anomaly zone space;

[0238] Step (4.1), drawing the anomaly diagram;

[0239] To visually demonstrate the spatial distribution characteristics of background anomalies within the study area and to verify the rationality of the anomaly delineation, spatial interpolation was performed based on the Mahalanobis distance scores of all samples to generate a continuously distributed comprehensive anomaly score map. At the same time, the boundaries of the anomaly area and known enriched mineral points were overlaid, providing an intuitive visual basis for geological interpretation and target area delineation.

[0240] Since the original sampling points are discretely distributed, spatial interpolation methods are needed to extend the discrete point data to the entire study area. Considering the discrete spatial distribution of sampling points, this embodiment uses the inverse-distance radial basis function (RBF) to spatially interpolate the anomaly scores of discrete sampling points to generate a continuous anomaly distribution map covering the entire area. Specifically, let the sampling point location be... ( ), corresponding to the abnormal score Then the point to be estimated Interpolation results at It can be represented as a linear combination of radial basis functions:

[0241]

[0242] in, For Euclidean distance, the basis functions are taken in inverse distance form. , It is the distance between the point to be estimated and the sampling point, and the weighting coefficient. Through precise interpolation conditions Solving the linear equations yields the result; the smoothing parameter is set to 0 to ensure that the interpolation surface strictly passes through all known sampling points; this highlights the influence of neighboring sampling points on the interpolation result while maintaining the spatial continuity of abnormal distributions;

[0243] To generate an anomaly score spatial distribution map covering the entire study area, an interpolation grid needs to be constructed first. To reduce the influence of boundary regions on the interpolation results, the grid range should typically be slightly larger than the actual distribution range of the sampling points. Based on the minimum and maximum latitude and longitude of all sampling points, the grid boundaries are extended by 10% outside this range. In terms of grid division, both the longitude and latitude directions are set to 500 nodes, ultimately generating a uniformly distributed rectangular interpolation grid. To eliminate outliers during the interpolation process, the interpolation results are cropped, limiting them to the minimum and maximum values ​​of the original data.

[0244] Step (4.2): Delineation of the abnormal area;

[0245] Geochemical anomalies typically exhibit spatial clustering characteristics. Compared to isolated anomalies, spatial clusters of multiple anomalies often better reflect the true distribution pattern of anomalies and have greater significance for mineral exploration.

[0246] Considering that individual outliers may be affected by local noise or random perturbations, this embodiment introduces the density-based spatial clustering algorithm DBSCAN to perform cluster analysis on outliers. DBSCAN can identify clusters of arbitrary shapes based on the local density of sample points and mark samples in low-density areas as noise points, so it is suitable for identifying outlier clusters under complex spatial distribution conditions. The key parameters of this algorithm include the neighborhood radius and the minimum number of points, where the neighborhood radius is used to define the neighborhood range of sample points and the minimum number of points is used to determine the core points. The clustering results are gradually formed through the density reachability relationship and density connectivity relationship between core points.

[0247] To reduce manual intervention, an adaptive method is used to estimate the neighborhood radius ε. For the set of outliers, the distance between each point and its k-th nearest neighbor is calculated, resulting in a distance array d. knn Find the median d of the array. med The base radius is then multiplied by an adjustment factor. The final neighborhood radius ε is obtained as follows:

[0248]

[0249] in, Used to control the tightness of clustering, setting The minimum number of points is set to 5, and k is set to 4 to ensure that the cluster has a certain size and to avoid pseudo-clusters composed of a few points.

[0250] Input the coordinates of outliers into the DBSCAN algorithm to obtain the cluster label for each point; points with a cluster label of -1 are considered noise and do not participate in the subsequent generation of outlier regions; through clustering, spatially clustered outliers are merged into the same cluster, while isolated points are removed.

[0251] To ensure that the influence range of each point is related to its density and to avoid oversmoothing or underfitting due to a fixed radius, for all points within each anomaly cluster, a radius is set with its planar coordinates as the center. Generate a buffer polygon; the buffer radius setting needs to consider both the distribution range of the point group and the mapping requirements, as shown in the table below:

[0252]

[0253] in, The neighborhood radius of DBSCAN. Let be the buffer coefficient, and take . It is 0.35;

[0254] Geometrically merge the buffer polygons of all points within the same cluster to obtain a preliminary covering polygon;

[0255] Since the area covered by the point group may have internal gaps or uneven edges, a morphological closing operation needs to be performed on the merged polygon. A positive buffering operation is then performed on the polygon to connect adjacent separated regions and fill internal holes. A negative buffering operation is then performed on the positively buffered result to remove isolated small protrusions and smooth the boundaries. Finally, a complete covering polygon with no internal gaps and smooth edges is obtained.

[0256] The generated polygons were simplified using Douglas-Peucker simplification with a simplification tolerance of 200m. The simplified polygons still retained the main shape of the point group and removed unnecessary details, as shown in Table 2 below.

[0257] Table 2 Statistics of data in abnormal areas

[0258] 1 0.611 0.780 7.02 2 2 0.575 2.150 8.62 1

[0259] Note: 5.91 is used as the anomaly detection threshold, meaning that samples with a comprehensive anomaly score exceeding 5.91 are considered to have reached an anomaly level of over 95%. The number of enriched mineral deposits indicates the number of known enriched mineral deposits within each anomaly zone.

[0260] The table above shows the main statistical characteristics of the anomalous areas. The average outlier values ​​for the two anomalous areas are 7.02 and 8.62, respectively, both higher than the anomalous threshold of 5.91 set in this embodiment, indicating that the overall anomalous level has reached over 95%. Anomalous area 1 contains two enriched mineral deposits; anomalous area 2 contains one enriched mineral deposit. Its average and maximum TiO2 values ​​are both higher than those of anomalous area 1, and its average outlier value is the highest, indicating a stronger anomalous intensity and more pronounced local enrichment in this area. The anomalous areas show a good spatial correspondence with known mineralization points.

[0261] Step (5), Model Evaluation;

[0262] To systematically evaluate the performance of this application in mineral exploration anomaly identification, the following experimental evaluation scheme can be adopted.

[0263] Based on industrial grade sample information, a spatial nearest neighbor matching method is used to establish the correspondence between samples and model output anomaly scores, thereby obtaining true labels that can be used for quantitative evaluation. The model's recognition ability is evaluated from different dimensions, using metrics including ROC curves and AUC values ​​based on continuous anomaly scores, and accuracy (ACC) and recall based on binary classification results. Three representative anomaly detection algorithms are selected as comparison methods, and comparative experiments are conducted under a unified dataset and consistent evaluation process to improve the fairness and reproducibility of the experimental results.

[0264] Step (5.1): Spatial matching verification based on known mineral deposits;

[0265] In actual exploration, known mineral deposits are an important basis for evaluating the effectiveness of anomaly identification methods. However, mineral deposit data and sampling points do not correspond strictly; there may be spatial offsets between the two.

[0266] To establish the correlation between mining sites and model prediction results, this embodiment employs a spatial nearest neighbor matching method, associating each labeled sample with the nearest test set sample. Let the set of labeled samples be... , where N m Here, j represents the number of labeled samples, and j is the sample index. Each labeled sample contains geographic coordinates. and binary tags Where 1 indicates an anomaly at the mining site, and 0 represents a normal occurrence. The test set sample set is... , where N t `i` represents the number of samples in the test set, and `i` is the sample index. Each sample contains coordinates. and the anomaly scores calculated by each algorithm. To establish spatial relationships, the geographic coordinates of all points are projected onto the UTM coordinate system to obtain planar coordinates. For each labeled point Calculate the Euclidean distance between it and all test set samples, and find the closest point. .

[0267]

[0268] If the nearest distance is less than the set spatial threshold Then it is considered that the labeled sample is the same as the test set sample. Spatially matchable, and anomaly scores for the test set samples. As labeled sample points The predicted score is used; otherwise, it is considered that the sample point lacks corresponding sampling information and cannot be matched; in this embodiment, the predicted score is used. m. Through the above matching, a set of sample pairs with true labels and predicted scores is obtained. This is used for subsequent model evaluation. The matching and validation process is based solely on test set samples and does not use training set information, ensuring that the evaluation results are unbiased.

[0269] Known mineral deposit data is used only for performance evaluation after model training and does not participate in spatial domain partitioning, VAE background learning, Mahalanobis distance parameter estimation, or anomaly threshold determination. This embodiment is an unsupervised anomaly identification method. Spatial domain partitioning uses only the geochemical characteristics, sediment grain size characteristics, and spatial location of the sampling points themselves, without using known mineral deposit labels or constructing training positive or negative samples based on the neighborhood of known mineral deposits. Similarly, no mineral deposit labels are input during VAE model training; the model only learns the background distribution characteristics of samples within each spatial subdomain. In this embodiment, known mineral deposits are only used as external validation information to verify the spatial correspondence between the model's output anomaly score and actual mineralization information. Therefore, the model training process and the performance evaluation process are separated in terms of data usage, avoiding the circular reasoning problem caused by known mineral deposits participating in both training and evaluation.

[0270] Step (5.2), performance evaluation and parameter sensitivity analysis;

[0271] To verify the effectiveness of the DVAE-SC method, this embodiment selects three classic unsupervised anomaly detection algorithms for comparison, representing different anomaly detection paradigms.

[0272] Isolation Forest is an ensemble-based anomaly detection algorithm. Its core idea is that outliers are more easily isolated through random partitioning, resulting in shorter path lengths from the root node to the leaf node. The algorithm constructs multiple random trees, calculates the average path length for each sample, and normalizes it into an anomaly score. A higher score indicates a more likely anomaly. Because this method uses a random partitioning mechanism, it does not require strict assumptions about the data distribution, thus maintaining good applicability and computational efficiency even in high-dimensional and complex contexts.

[0273] Principal Component Analysis (PCA-Mahalanobis) is a linear dimensionality reduction method that projects the original features onto a new coordinate space through orthogonal transformation. The PCA-Mahalanobis method first performs PCA dimensionality reduction on the training set to obtain principal component scores; then, it reconstructs the dimensionality-reduced data and calculates the reconstruction error vector; finally, it uses Mahalanobis distance to synthesize the reconstruction error as anomaly scores. The principle of this method is similar to DVAE-SC, but it uses linear dimensionality reduction instead of the nonlinear modeling of VAE. Therefore, it is more suitable for data with significant linear relationships and has relatively limited ability to characterize complex nonlinear structures.

[0274] Gaussian Mixture Model (GMM) is a probabilistic generative model that assumes the data is a mixture of multiple Gaussian distributions. For each sample, GMM calculates the probability that it belongs to each component and outputs a log-likelihood value as an anomaly score; the lower the score, i.e. the smaller the log-likelihood, the less likely the sample is to be generated by the background distribution, and therefore the more likely it is to be an anomaly.

[0275] To measure the ability of each algorithm to distinguish between anomalies and background, this embodiment uses two types of evaluation metrics: ROC curves and AUC based on continuous scores, and traditional metrics based on binary classification, such as recall and precision.

[0276] The Receiver Operating Characteristic (ROC) curve is a classic tool for measuring the ranking ability of a classifier. It plots the false positive rate (FPR) on the horizontal axis and the true positive rate (TPR) on the vertical axis, drawing a curve that iterates through all possible classification thresholds. The Area Under the ROC Curve (AUC) is the area under the ROC curve, with values ​​ranging from 0.5 to 1. An AUC of 0.5 indicates that the model has no discriminative ability, while an AUC closer to 1 indicates a stronger ability to rank positive and negative samples. The advantage of AUC is that it is independent of a specific classification threshold and can objectively reflect the quality of the model's ranking of outlier scores. For the ROC curve, randomly selecting a threshold within the range of outlier scores can distinguish between outlier regions and normal background. As shown in the table below, using known mining sites as a reference, four indicators reflecting the outlier score results from different perspectives can be obtained for each sample's outlier score.

[0277] Table 3. Confusion Matrix of Classification Results

[0278] Including known mineral deposits Real example False counterexample No known mineral deposits False positives True counterexample

[0279] In this embodiment, for each algorithm, the matched sample pairs are used. Calculate the ROC curve and AUC value. A higher AUC value indicates that the algorithm is better able to rank genuine mining points before background points, meaning the algorithm has better discriminative power against outliers and a stronger ability to distinguish between mining points and background points. To ensure the fairness and reproducibility of the comparative experiments, all algorithms followed a unified experimental procedure and parameter settings. Parameter descriptions are shown in the table below:

[0280] Table 4 Parameter Description Table

[0281] DVAE-SC <![CDATA[For VAE training, the dimension of the hidden layer is 32, and the dimension of the latent variable is dynamically determined by latent_factor (default is 0.5). 5 models are trained for each domain and the optimal one is selected. The Adam optimizer is used, with a learning rate of 10 -4 , training is conducted for 2000 epochs, and the loss function is the sum of MSE reconstruction loss and KL divergence.]]> Isolation Forest contamination=0.15; random_state=42; the rest are default. PCA-Mahalanobis <![CDATA[PCA retains 95% variance; the Mahalanobis distance regularization coefficient is set to 10 consistent with DVAE-SC -6 > GMM n_components=3; covariance_type='full'; random_state=42

[0282] To further evaluate model performance, RBF interpolation was used to obtain anomaly maps for comparison of the four models. The anomaly maps obtained by different methods differed in anomaly area range, spatial continuity, and background discrimination. DVAE-SC identified a relatively concentrated distribution of anomalies, with good continuity in high-value areas, clearly reflecting the spatial aggregation characteristics of anomaly information and corresponding more closely to actual mining sites. While Isolated Forest could identify some anomaly locations, there was some diffusion in local areas, and the anomaly boundaries were relatively unclear. PCA-Mahalanobis could reflect the main anomaly centers, but the overall result was too smooth, with a relatively gentle transition between the anomaly area and the background. The anomaly response obtained by the GMM method was relatively scattered, with insufficient stability in local areas.

[0283] To simulate scenarios requiring explicit anomaly labels in real-world applications, the 95th percentile of the anomaly scores in the training set of each model was set as the classification threshold. Continuous scores were converted into binary prediction results, dividing the region into anomaly and normal background areas, and relevant evaluation metrics were obtained, such as... Figure 12 As shown.

[0284] As shown in the table below, this application achieved optimal results in all relevant evaluation metrics. Specifically, the AUC of 0.898 indicates that the method has a strong ability to distinguish between mineral deposit samples and background samples; the Recall of 0.889 indicates that it can effectively identify most known mineral deposits; and the ACC of 0.931 is the highest among all methods. These metrics demonstrate that the method has good accuracy and stability in the classification of mineral deposits and background points. DVAE-SC, by combining spatial domain partitioning with VAE nonlinear background modeling and using Mahalanobis distance to comprehensively express multidimensional error information, can effectively identify multivariate anomalies related to mineralization, with a relatively low false alarm level.

[0285] Table 5 Comparison of Performance Indicators

[0286] This application 0.898 0.889 0.931 IsolationForest 0.842 0.786 0.716 PCA-Mahalanobis 0.860 0.758 0.735 GMM 0.656 0.728 0.637

[0287] In summary, this embodiment details a variational autoencoder and spatial clustering anomaly identification method based on spatial domain partitioning. Through global spatial autocorrelation analysis and LISA local clustering feature analysis, it verifies the significant spatial heterogeneity of the study area background, demonstrating the necessity of spatial domain partitioning. Based on this, a spatial domain partitioning module is constructed. By combining multivariate K-means clustering, spatial connectivity component splitting, and small-domain fusion strategies, the study area is divided into five relatively uniform and spatially continuous subdomains. Subsequently, the variational autoencoder independently learns the background within each subdomain. The reconstruction results of most background samples are close to the measured values, while anomalous samples show more significant deviations. Mahalanobis distance is introduced to integrate multidimensional reconstruction errors, constructing a comprehensive anomaly score index. Mahalanobis distance, through the inverse mapping of the covariance matrix, automatically eliminates the interference of inter-element correlation and dimensions, making the quantification of anomaly degree more objective. Furthermore, DBSCAN spatial clustering and morphological closing operations are used to extract the boundaries of anomalous regions and spatially overlay them with the interpolated anomaly map, providing a reference for delineating target areas for marine placer mineral exploration.

[0288] The performance of this method can be intuitively and accurately verified through spatial matching and comparison experiments with known mineral deposits. The results show that this application significantly outperforms Isolation Forest, Principal Component Analysis-Mahalanobis (PCA-Mahalanobis), and Gaussian Mixture Model (GMM) in multiple indicators such as AUC (0.898), recall (0.889), and precision (0.931). Parameter sensitivity analysis further demonstrates that this application has strong robustness to fluctuations in the number of clusters and the dimension of latent variables, and can effectively solve the problem of anomaly identification and spatial delineation in complex backgrounds, providing a practical and feasible technical route for the delineation of target areas for marine placer deposit exploration.

[0289] This application also proposes a novel electronic device, which includes a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the program, it implements the above-mentioned variational autoencoder based on spatial domain partitioning and spatial clustering anomaly identification method.

[0290] This application also proposes a novel computer-readable storage medium storing a computer program that, when executed, enables the implementation of the aforementioned variational autoencoder based on spatial domain partitioning and spatial clustering anomaly identification method.

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

Claims

1. A variational autoencoder based on spatial domain partitioning and a spatial clustering anomaly identification method, characterized in that: Establish a system architecture that includes the following modules: The spatial domain partitioning module is used to decompose the entire study area into multiple local regions with relatively uniform backgrounds, so as to facilitate subsequent modeling within local subdomains. The domain-specific VAE background learning module is used to select the VAE model with the smallest reconstruction error in each domain through a multi-model optimization strategy, and reconstruct all samples in the domain to obtain the reconstruction error vector in the standardized space. The anomaly scoring module uses Mahalanobis distance to convert the multidimensional reconstruction error of each sample into a comprehensive score, while eliminating the influence of the scale and correlation between features, so that the degree of anomaly of different samples can be directly compared and ranked. The anomaly zone delineation module determines anomaly points based on thresholds and transforms discrete points into anomaly zones with clear geographical boundaries through spatial clustering and boundary generation. Based on the above system structure, firstly, geochemical data and sediment grain size data are input into the spatial domain partitioning module to divide the study area into several spatially continuous subdomains with similar backgrounds and output domain labels. Then, the data of each domain are input into the VAE background learning module within the domain, and the reconstruction error vector is output. Secondly, the error vectors of all training samples are input into the anomaly scoring module, where the error between the sample input value and the reconstructed value represents the degree of deviation from the background. The inverse covariance matrix is ​​calculated to obtain the anomaly score for each sample. Finally, based on the threshold determined by the training set, samples with anomaly scores higher than the threshold are selected and sent to the spatial clustering anomaly delineation module. After cluster analysis and boundary extraction, the anomaly polygon is output.

2. The variational autoencoder based on spatial domain partitioning and spatial clustering anomaly identification method according to claim 1, characterized in that: Known mining site data is introduced, and the results are correlated with test set samples through spatial nearest neighbor matching. ROC curve analysis is then performed on the outlier scores to verify the effectiveness.

3. The variational autoencoder based on spatial domain partitioning and spatial clustering anomaly identification method according to claim 1, characterized in that: Includes the following steps, Step (1), spatial domain division; Step (1.1): Modeling different geochemical backgrounds; The study area was divided into sub-regions with relatively uniform backgrounds, and background estimation and anomaly identification were performed in each sub-region. The global Moran's I and Geary's C indices were used to test spatial autocorrelation. The LISA index was introduced to decompose the global autocorrelation into each spatial unit, revealing the spatial correlation pattern and heterogeneity characteristics of TiO2 content at the local scale. Step (1.2): Multivariate feature clustering based on K-means; Based on the nature of the variables, the characteristics of the variables are divided into two categories: one is compositional variables, including major element oxides and trace elements; the other is physical parameter variables, namely particle size parameters. For component-type variables, the percentage content data are first uniformly converted, and then processed using a centered logarithmic ratio transformation. For the granularity parameter, the following Z-score standardization is used to eliminate the dimensional differences between different variables: in, These are the original particle size parameter values. and These are the mean and standard deviation of the corresponding variables, respectively. This is the standardized result; Finally, the constituent variables processed by CLR are concatenated column-wise with the standardized granular parameters to form the input feature matrix for in-domain VAE training. , where n represents the number of samples in the subdomain and m represents the dimension of the input feature; The K-means algorithm is used to perform initial clustering of the multivariate geochemical and granular characteristics of the samples; K-means is achieved by iteratively optimizing the following objective function: in, Indicates the first Clusters, It is the center of the cluster; The input features of K-means are geochemical data after centering log-ratio transformation and granular data after Z-score normalization. , where m is the feature dimension, and the feature variables include geochemical indicators related to mineralization and grain size parameters representing sedimentary characteristics; Step (1.3): Spatial connectivity component splitting and small domain fusion; Based on K-means clustering, samples that are spatially disconnected within the same K-means cluster are split into different domains; First, construct a spatial adjacency graph. Where the vertex set V corresponds to all samples, and the edge set E represents the sample... Spatial adjacency relationship between samples; the condition for adjacency between samples i and j is the Euclidean distance d between them. ij Less than the given connection radius r conn ; Then, after obtaining the spatial adjacency graph, for each K-means class c, consider the induced subgraph. ,in, Let c be the set of samples belonging to class c. for Edges between samples; extracted using breadth-first search or depth-first search. All connected components; samples within each connected component are spatially connected and have similar geochemical composition and sediment grain size structure characteristics, and are therefore defined as an independent candidate domain, thus splitting each initial K-means class into several spatially continuous candidate domains; Secondly, the aforementioned smaller domains are merged into neighboring larger domains, and a minimum domain size threshold N is set. min When the number of samples in the candidate domain is less than the minimum domain size threshold N min The fusion operation is triggered at this time; After the above fusion, each domain obtained is not only spatially continuous and similar in geochemical composition and grain size structure, but also has a sufficient sample size to support subsequent statistical modeling. All domains are renumbered consecutively starting from 0 to obtain the final domain label for each sample: domain∈{0,1,2,3……,D-1}, where D is the number of final spatial domains. Step (2), variational autoencoder background learning; Variational autoencoders are introduced in each subdomain to learn and reconstruct background features; A VAE model is trained in each subdomain to represent the background distribution characteristics unique to that domain. The preprocessed feature samples are input into the trained VAE to obtain the corresponding reconstruction results. The reconstructed values ​​output by the model are regarded as the background estimates of the sample under the background constraints of the subdomain, and the subsequent anomaly score is calculated accordingly. Step (3): Calculate the overall anomaly score; Mahalanobis distance is introduced to comprehensively measure multidimensional error information. Given a multivariate dataset, let its mean vector be... The covariance matrix is For any sample point The Mahalanobis distance from the dataset center is defined as: The covariance matrix of the estimation error vector when calculating Mahalanobis distance This matrix reflects the distribution structure of the training set background samples in the multidimensional error space, including the variance of each error component and the covariance between the components. Let the training set have a total of For each sample, the set of its reconstruction error vectors is: ,in Calculate the mean of the error vector: Since VAE reconstructs background samples unbiased, It is usually very close to the zero vector; the unbiased estimator of the covariance matrix is: The covariance matrix is A symmetric positive definite matrix whose diagonal elements Indicates the first The variance of each error component, off-diagonal elements Indicates the first The and the first The covariance of each error component; After obtaining the regularized covariance matrix Then, the inverse matrix is ​​calculated using a matrix inversion algorithm. The inverse matrix will be used to calculate the Mahalanobis distance for all subsequent samples. The estimation of the covariance matrix is ​​based only on the training set samples, and the test set samples are not included in the calculation. After obtaining the inverse covariance matrix, the outlier score of each sample is defined as the Mahalanobis distance of its error vector; for any sample, the reconstructed error vector is e, and the mean of the error vector is... The anomaly score S of this sample is defined as: The anomaly score S measures the degree to which the sample error vector deviates from the center of the background error distribution. The larger the value, the more likely the sample is to be anomaly. For background samples, their error vector e is usually located in the dense region of the background distribution, with a smaller Mahalanobis distance and a lower score. For anomaly samples, their error vector deviates from the center of the background distribution, resulting in a significantly larger Mahalanobis distance and a higher score. After obtaining the Mahalanobis distance scores of all samples in the training set, their distribution is analyzed to determine a reasonable anomaly threshold. Step (4): Delineation of the anomaly zone space; Step (4.1), drawing the anomaly diagram; Spatial interpolation is performed based on the Mahalanobis distance scores of all samples to generate a continuously distributed comprehensive anomaly score map. Spatial interpolation methods are used to extend discrete point data to the entire study area; inverse distance radial basis functions are used to spatially interpolate the anomaly scores of discrete sampling points to generate a continuous anomaly distribution map covering the entire area. Step (4.2): Delineation of the abnormal area; The density-based spatial clustering algorithm DBSCAN performs cluster analysis on outliers; An adaptive method is used to estimate the neighborhood radius. For each outlier in the set, calculate the distance between each point and its k-th nearest neighbor to obtain the distance array d. knn Take the median d of the array. med The baseline radius is then multiplied by an adjustment factor. The final neighborhood radius ε is obtained as follows: in, Used to control the tightness of clustering, setting The minimum number of points is set to 5, and k is set to 4 to ensure that the cluster has a certain size and to avoid pseudo-clusters composed of a few points. Input the coordinates of outliers into the DBSCAN algorithm to obtain the cluster label for each point; points with a cluster label of -1 are considered noise and do not participate in the subsequent generation of outlier regions; through clustering, spatially clustered outliers are merged into the same cluster, while isolated points are removed. For all points within each anomaly cluster, with its planar coordinates as the center and radius... Generate a buffer polygon; the buffer radius should be set to balance the distribution range of the point group and the mapping requirements, as shown in the following formula: in, The neighborhood radius of DBSCAN. Let be the buffer coefficient, and take . It is 0.35; Geometrically merge the buffer polygons of all points within the same cluster to obtain a preliminary covering polygon; The merged polygons are subjected to morphological closing operations, and positive buffering is performed on the polygons. Then, negative buffering is performed on the positively buffered result. The generated polygons are then simplified using Douglas-Peucker simplification.

4. The variational autoencoder based on spatial domain partitioning and spatial clustering anomaly identification method according to claim 3, characterized in that: In step (1.2), values ​​less than or equal to 0 are replaced with a very small constant; subsequently, the following CLR transformation is performed on each sample: in, Let be the value of the i-th component variable in the sample, and D be the number of component variables; This is the result of the variable after CLR transformation.

5. The variational autoencoder based on spatial domain partitioning and spatial clustering anomaly identification method according to claim 3, characterized in that: Step (1.3) describes the adaptive determination of the connection radius r. conn For each sample, calculate its k-th... nn Find the distance to the nearest neighbor, and then take the median d of that distance for all samples. med And multiply by a coefficient , Let be the radius of the connection, then: 。 6. The variational autoencoder based on spatial domain partitioning and spatial clustering anomaly identification method according to claim 3, characterized in that: Step (1.3) includes the following fusion strategy process, Graph adjacency-based fusion: First, using the constructed spatial adjacency graph G, we count the domain labels of other samples in the neighborhood of each sample in the small domain. Then, for the small domain D small Count the non-D occurrences in the neighborhood of all its samples small The domain labels are selected, and the domain with the highest frequency of occurrence is chosen as the target domain. small All samples are merged into the target domain; This process is repeated until no new small domains are generated or the maximum number of iterations is reached. Forced fusion based on spatial nearest neighbor: If isolated small regions still exist after graph adjacency fusion, then the spatial nearest neighbor method is used for forced fusion; specifically, for each sample in the small region, its k-th neighbor is searched globally. nn For each nearest neighbor, the domain label of the sample is changed to the largest domain label that appears most frequently among these nearest neighbors; if there is no large domain among the nearest neighbors, then the domain label that appears most frequently is taken; among these, considering both fusion stability and spatial local representativeness, the most frequent domain label is taken. The value is 20. This parameter is used to determine the local spatial neighborhood range of isolated small-domain samples, using the labels of large-domain or majority-domain samples within the neighborhood as the basis for fusion.

7. The variational autoencoder based on spatial domain partitioning and spatial clustering anomaly identification method according to claim 3, characterized in that: Step (2) includes, Step (2.1), Intra-domain VAE training; A single-hidden-layer symmetric VAE structure is adopted. First, the input features are extracted through a fully connected hidden layer, and then the mean vectors of the latent variables are output separately. Sum of logarithmic variance vector The hidden layer uses the ReLU activation function to extract nonlinear information from the input features; for the input sample x, the encoding process can be represented as: in, For the input sample vector, and These are the weight matrix and bias vector from the input layer to the hidden layer, respectively. Output for hidden layer; Based on this, the distribution parameters of the latent variables are further output: in, This represents the mean vector of the latent variables. The logarithmic form of the variance of the latent variables. As a potential spatial dimension; and These are the weight matrix and bias vector of the corresponding output layer, respectively; The decoder structure is basically symmetrical to the encoder, also containing a hidden layer and an output layer. The hidden layer uses the ReLU activation function. Since the input features have undergone CLR transformation and standardization, their numerical range is no longer limited to a fixed interval. Therefore, a linear form for the output layer is more in line with the numerical characteristics of the feature space. The reconstruction process is as follows: in, For the reconstruction result, W4 and b4 are the parameters from the latent layer to the decoded hidden layer, and W3 and b3 are the parameters from the decoded hidden layer to the output layer. The network size is dynamically determined by the dimension of the input features; let the dimension of the input features be m, then the dimension of the hidden layers is set as: Where h1 represents the number of neurons in the hidden layer. This means taking the larger of the two values; this setting ensures that the network has sufficient expressive power on the one hand, and avoids the network being too narrow when the feature dimension is small on the other hand. In the latent space, sampling is performed using a reparameterization technique. The dimension of the latent space is defined as: Where b represents the intermediate scale obtained by compressing the input dimension. Indicates rounding down; The potential spatial dimension is set as follows: in, As a potential spatial dimension, This is the latent dimension scaling factor, used to control the degree of compression of the latent space; Step (2.2), reconstruction error analysis; Reconstruct all samples within the domain, including the training and test sets, and calculate the reconstruction error vector; Let the input vector of the sample in the processed feature space be... The corresponding reconstruction result is ;in, Indicates the input feature dimension. These are the actual feature values ​​of the sample. The reconstructed value given by VAE; The reconstruction error vector of this sample is defined as: in, To reconstruct the error vector, Indicates the first Reconstruction bias of each feature; if This indicates that the actual value of the feature is higher than the background reconstruction value given by the model; if If the value is less than the background reconstruction value, it means that the actual value of the feature is lower than the background reconstruction value; the larger the absolute value of the error vector, the more obvious the deviation of the sample from the background pattern on that feature.

8. The variational autoencoder based on spatial domain partitioning and spatial clustering anomaly identification method according to claim 7, characterized in that: In step (4.1), the sampling point location is set as follows: ( ), corresponding to the abnormal score Then the point to be estimated Interpolation results at It can be represented as a linear combination of radial basis functions: in, For Euclidean distance, the basis functions are taken in inverse distance form. , It is the distance between the point to be estimated and the sampling point, and the weighting coefficient. Through precise interpolation conditions The linear equations are solved to obtain the smoothing parameter, which is set to 0 to ensure that the interpolated surface strictly passes through all known sampling points.

9. An electronic device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that: When the processor executes a computer program, it implements the variational autoencoder and spatial clustering anomaly identification method based on spatial domain partitioning as described in any one of claims 1 to 8.

10. A computer-readable storage medium storing a computer program, characterized in that: When the computer program is executed, it can implement the variational autoencoder and spatial clustering anomaly identification method based on spatial domain partitioning as described in any one of claims 1 to 8.

Images