A rock stratum automatic identification method and system based on multispectral remote sensing

Abstract

The present application relates to the field of artificial intelligence and data processing technology, in particular to a rock stratum automatic identification method and system based on multispectral remote sensing, specifically as follows: the adaptive normalization method based on continuum removal is used for the collected multispectral remote sensing data samples; the mixed pixel spectrum based on the physical mechanism of rock mass is used for decomposing the normalized samples, and the abundance coefficient of each end member is calculated; the first derivative value of each normalized sample in each wave band is calculated, and the fusion characteristic value of each dimension is calculated in combination with the abundance coefficient; the rock stratum identification model based on deep learning is constructed, the fusion characteristic value of each dimension is input for rock stratum category probability prediction; the rock stratum identification model is optimized and trained; the newly collected data is input into the trained model after processing, and the final identification result is obtained. The present application can improve the accuracy, stability and geological rationality of the rock stratum automatic identification result.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of artificial intelligence and data processing technology, and in particular to a rock stratum automatic identification method and system based on multispectral remote sensing. BACKGROUND

[0002] In geological exploration and mineral resources investigation, accurate identification of rock stratum distribution is of great significance for guiding field investigation, resource development and geological modeling. Multispectral remote sensing technology has become an important means for rock stratum identification due to its wide coverage and non-contact acquisition capability. However, existing rock stratum identification methods based on multispectral remote sensing still face many challenges in practical application. First, different rock strata have significant differences in background reflectivity, and traditional min-max normalization method cannot preserve local absorption characteristics, making it difficult to distinguish dark rocks and low reflectivity minerals. Second, remote sensing pixels often have mixed reflectance of multiple lithologies, and existing linear mixing models assume fixed endmembers, which cannot capture the nonlinear mixing effects caused by mineral composition variation. Third, conventional feature construction methods simply concatenate original spectra and abundance information without considering their physical sources and expression differences, resulting in insufficient model identification capability for key lithology characteristics. In addition, most deep learning models rely only on classification accuracy for training, ignoring the physical mechanism of mineral mixing that rock stratum identification should satisfy, leading to lack of model interpretability. Therefore, a new method that combines physical constraints and adaptive expression of spectral features is urgently needed to improve the accuracy, stability and geological reasonableness of rock stratum identification.

[0003] Chinese invention patent with publication number CN117830203A proposes a comprehensive remote sensing research field of basic geology, specifically relates to a rock stratum occurrence extraction method based on DEM and high-resolution optical remote sensing data. The main content includes: automatically extracting contour lines of corresponding accuracy according to DEM data of the study area using a computer; automatically extracting rock stratum boundary lines from high-resolution optical remote sensing data of the study area and performing manual correction; extracting intersection data of rock stratum boundary lines and contour lines; automatically extracting intersection data with more than or equal to 3 intersection points between a rock stratum surface and adjacent two contour lines; calculating different rock stratum surface inclination and dip angle parameters using the three-point coplanar method; labeling each point occurrence along the boundary line for each rock stratum surface in the study area and outputting in the form of a graph or table; efficiently and quickly obtaining rock stratum occurrence data in a large area; and outputting in the form of a graph or table for further analysis and utilization.

[0004] The prior art has the following objective defects: the conventional min-max normalization ignores local absorption characteristics, which easily leads to compression of dark rock characteristics and reduction of classification accuracy; the conventional mixed pixel processing relies on fixed end members and cannot cope with actual mineral composition fluctuations and nonlinear mixing effects; the conventional feature splicing method does not distinguish the physical differences between derivatives and abundances, which makes it difficult for the model to focus on key lithology characteristics; and the use of only cross-entropy loss ignores physical constraints, and the prediction result lacks physical consistency guarantee for rock layer mixing mechanism.

[0005] Therefore, the present application proposes a rock layer automatic identification method and system based on multispectral remote sensing to solve the above problems. SUMMARY

[0006] The present application is aimed at the deficiencies of the prior art and develops a rock layer automatic identification method based on multispectral remote sensing. The present application can improve the accuracy, stability and geological rationality of rock layer automatic identification results.

[0007] In one aspect, the technical scheme for solving the technical problem of the present application is a rock layer automatic identification method based on multispectral remote sensing, which is as follows:

[0008] An adaptive normalization method based on continuum removal is used for the collected multispectral remote sensing data samples to obtain normalized samples. The collected multispectral remote sensing data samples are derived from open source data sets or self-collected data sets, and the collected multispectral remote sensing data samples are all labeled with lithology class labels. The adaptive normalization method based on continuum removal is as follows:

[0009] For the original reflectivity of each normalized sample in each waveband, first, the influence of atmospheric scattering and illumination variation is eliminated by fitting a spectral envelope line. Specifically, a cubic polynomial fitting function is used, the center wavelength and the original reflectivity value of the input waveband are input, and the fitting value is obtained. Then, the maximum value of the original reflectivity value and the fitting value is taken as the continuum function value of each sample in each waveband, which is used to represent the spectral envelope line.

[0010] Next, the local absorption depth is normalized based on the difference between the continuum function value and the minimum reflectivity to obtain the reflectivity value of each normalized sample in each waveband, i.e. the normalized sample.

[0011] Based on the physical mechanism of the rock mass, the normalized sample is decomposed, each end member after decomposition is initialized, and the abundance coefficient of each end member is calculated. The abundance coefficient of each end member is calculated as follows:

[0012] A clustering method is used to generate sample clusters for each end member after decomposition, and a indicator function is used to determine whether the current waveband matches the known mineral spectral library based on the normalized samples, and the initial reflectance value of each end member at each waveband is calculated;

[0013] Then, the initial reflectance value of the fixed end member is calculated by minimizing the L2 norm square of the linear combination of the normalized reflectance value and each end member, and the abundance coefficient of each sample in each end member is calculated by combining the coefficient regularization term;

[0014] Wherein, the end member represents different rock types.

[0015] The first derivative value of each normalized sample at each waveband is calculated, and the multi-dimensional fusion feature value is obtained by combining the abundance coefficient and the characteristic fusion of the regularization attention weight;

[0016] For each waveband of each sample, the difference value of the normalized reflectance value in the adjacent waveband is calculated, and the first derivative value is obtained by combining the center wavelength difference value of the adjacent waveband and the derivative scaling factor;

[0017] Then, for each sample, the attention weight of the first derivative value in each dimension and the attention weight of the abundance coefficient are calculated, and then the weighted first derivative feature value and the weighted abundance coefficient are summed to obtain the fusion feature value in each dimension;

[0018] Wherein, the dimension includes the total number of wavebands and the total number of end members.

[0019] A rock layer identification model based on deep learning is constructed, and the fusion feature value in each dimension is predicted based on the spectral adaptive activation function and the multi-scale residual block architecture in the model. The operation in the rock layer identification model based on deep learning is as follows:

[0020] The fusion feature value in each dimension is input into the model, and the fusion feature value in each dimension is input into the linear layer, and the projection weight matrix of the input layer and the bias vector of the input layer are combined to obtain the initial activation value of each sample at each hidden layer node. The initial activation value is input into the multi-scale residual block architecture, and the activation function value of the previous layer residual block is combined with the activation function value of the first two layers of residual blocks multiplied by the cross-layer jump connection weight matrix to obtain the activation function value of the current layer residual block. After layer-by-layer calculation, the activation function value of the final layer residual block is obtained. The activation function value is combined with the output layer weight matrix and the output layer bias vector, and the prediction probability of the rock layer category is calculated by the Softmax function.

[0021] The lithological identification model is optimized by incorporating the physical constraint of spectral decomposition into the loss function, and the model parameters are updated by using a parameter grouping adaptive learning rate mechanism, to obtain the trained lithological identification model.

[0022] The physical constraint of spectral decomposition is incorporated into the loss function, and the total loss is calculated by combining the cross-entropy classification loss term and the physical regularization term.

[0023] The cross-entropy loss term is calculated by quantifying the deviation between the predicted probability of the lithological category calculated by the model and the true category.

[0024] For each sample, the L2 norm square of the linear combination of the normalized reflectance value, the predicted abundance coefficient and the endmember reflectance value is calculated, and the physical regularization term is obtained by combining the absolute value of the abundance coefficient.

[0025] The parameter grouping adaptive learning rate mechanism is used to group the trainable parameters according to their physical role in the model, and the dynamic adjustment of the cosine annealing learning rate is applied to the weight parameters that control the feature transformation to stabilize the convergence of the deep network, while the fixed small learning rate is applied to the bias parameters that control the offset to suppress overfitting.

[0026] The training process of the model is as follows:

[0027] The collected data is forward propagated, processed and input into the lithological identification model to obtain the predicted probability of the lithological category, and then the loss function is calculated, and the model parameters are updated through back propagation, and the model is iteratively optimized by the Adam optimizer until the iteration termination condition is met, the optimization is ended, and the trained lithological identification model is obtained.

[0028] The newly collected multi-spectral remote sensing data is processed and input into the trained model to obtain the final lithological category probability prediction result, and the lithological category corresponding to the maximum probability value is taken as the identification result.

[0029] On the other hand, the present application also provides a lithological automatic identification system based on multi-spectral remote sensing, which comprises modules for executing processing instructions of each step in the lithological automatic identification method based on multi-spectral remote sensing, comprising the following modules:

[0030] Data collection module: multi-spectral remote sensing data samples collected;

[0031] Data preprocessing module: normalizing, mixed pixel decomposition and special fusion operation on the collected samples;

[0032] Lithological identification module: inputting the preprocessed samples into the lithological identification module to predict the lithology of the samples;

[0033] An optimized training module: the stratum identification module parameters are optimized and updated by a model total loss function and an adaptive learning rate mechanism, so as to obtain a trained stratum identification module;

[0034] An identification result output module: after the sample to be identified is preprocessed, the sample is input into the trained stratum identification module, and a final identification result is output.

[0035] The effects provided in the summary are only the effects of the embodiments, not all the effects of the application, and the above technical solutions have the following advantages or beneficial effects:

[0036] The present application considers that the multispectral remote sensing data is in the form of one-dimensional spectral curve, is composed of reflectivity values of multiple bands, has the characteristics of high dimension, strong noise interference and inconsistent scale caused by light and atmospheric conditions, adopts a continuum removal normalization method, and the background reflection difference can be eliminated by continuum function calculation, that is, the influence of atmospheric scattering and light change is eliminated, the absorption feature integrity is retained, the key absorption feature of lithology is highlighted, and the identification ability of dark rock and low-reflectivity minerals is improved.

[0037] A physically constrained mixed pixel decomposition model is constructed, mixed reflection signals of multiple strata contained in a single pixel are processed, end members are initialized and combined with key absorption bands of a mineral spectral library, the geological interpretability of the decomposition result is enhanced, the problem that the end member spectrum dynamically changes due to mineral composition variation and causes insufficient decomposition accuracy and inability to capture nonlinear mixing effects can be avoided.

[0038] The present application adopts a first-order derivative and abundance coefficient double-path attention fusion mechanism, the first-order derivative of the spectrum is calculated to highlight the absorption edge features of the lithology, then the weighted derivative features and the weighted abundance coefficients are summed, so as to adaptively strengthen the high-discriminative features and avoid feature dimension redundancy and physical semantic confusion.

[0039] By jointly optimizing the classification loss term and the physical regularization term, the spectral reconstruction error and the abundance sparsity are included in the training, the prediction result of the model is forced to comply with the physical mechanism of rock and mineral mixing, and the physical consistency and classification robustness of the model are improved, so as to avoid the problem that the prediction result of the model is disconnected with the physical mechanism of lithology and causes decomposition and prediction fragmentation.

[0040] In summary, the present application can improve the accuracy, stability and geological rationality of the automatic stratum identification result. BRIEF DESCRIPTION OF DRAWINGS

[0041] The accompanying drawings are included to provide a further understanding of the present application, and constitute a part of the specification, together with the embodiments of the present application, to explain the present application, and do not constitute a limitation on the present application.

[0042] Figure 1This is a schematic diagram of the method flow of the present invention.

[0043] Figure 2 This is a schematic diagram of the reflectance of the original spectral lines and the continuous envelope.

[0044] Figure 3 This is a schematic diagram of the spectral reflectance after continuous normalization.

[0045] Figure 4 This is a schematic diagram illustrating the impact of the normalization method of this invention on the accuracy of rock strata identification compared to conventional normalization methods.

[0046] Figure 5 This is a schematic diagram of the spectral lines of the physically constrained endmember.

[0047] Figure 6 This is a schematic diagram of the abundance distribution of mixed pixels.

[0048] Figure 7 This is a schematic diagram of the attention weights of the spectral derivative.

[0049] Figure 8 This is a schematic diagram of the attention weights for abundance coefficients.

[0050] Figure 9 This is a schematic diagram of the distribution of rock strata in the feature space. Detailed Implementation

[0051] To clearly illustrate the technical features of this solution, the invention will be described in detail below through specific implementation methods and in conjunction with the accompanying drawings.

[0052] Example 1

[0053] like Figure 1 As shown, an automatic rock strata identification method based on multispectral remote sensing includes the following steps:

[0054] S1. Collected multispectral remote sensing data samples, and label the collected multispectral remote sensing data samples with lithological category labels;

[0055] S2. An adaptive normalization method based on continuum removal is used to process the collected multispectral remote sensing data samples to obtain normalized samples.

[0056] S3. Based on the physical mechanism of the rock mass, the mixed pixel spectrum is decomposed into normalized samples, each endmember of the decomposed sample is initialized, and the abundance coefficient of each endmember is calculated.

[0057] S4. Calculate the first derivative value of each normalized sample in each band, and combine it with the abundance coefficient memory attention weighted feature fusion to obtain multi-dimensional fusion feature values.

[0058] S5, construct a rock stratum identification model based on deep learning, and predict the rock stratum class probability of the fusion feature value of each dimension based on the spectral adaptive activation function and multi-scale residual block architecture in the model;

[0059] S6, optimize the rock stratum identification model by incorporating the physical constraints of spectral decomposition into the loss function, update the model parameters using the adaptive learning rate mechanism of parameter grouping, and obtain the trained rock stratum identification model;

[0060] S7, input the processed multi-spectral remote sensing data into the trained model to obtain the final rock stratum class probability prediction result, and take the lithology class corresponding to the maximum probability as the identification result.

[0061] In the specific implementation, S1 is specifically as follows:

[0062] The collected data can use an open-source standard dataset or a self-collected dataset;

[0063] The open-source standard dataset can use Landsat-8-OLI open-source dataset, which contains 7 visible to short-wave infrared bands, wavelength range 430-2300 nanometers, spatial resolution 30 meters, which meets the requirements of spectral resolution and spatial coverage for rock stratum identification; The open-source standard dataset already has labels, including labels of magmatic rocks such as granite and basalt, labels of sedimentary rocks such as carbonates and clastic rocks, labels of metamorphic rocks such as gneiss and marble, etc.

[0064] If the self-collected dataset, the standardized multi-spectral remote sensing image dataset covering the target area needs to be obtained through satellite or aerial platform, and the data source needs to meet the following core requirements:

[0065] 1) Band integrity: must contain key bands sensitive to lithology (such as Band 1-Band 7), especially emphasizing the short-wave infrared band, because it is sensitive to the absorption characteristics of silicate and carbonate minerals;

[0066] 2) Radiometric calibration and atmospheric correction: the original data needs to be converted to apparent reflectance through radiometric calibration to eliminate atmospheric scattering effects and ensure the accuracy of the physical meaning of reflectance;

[0067] 3) Regional representation: the data needs to cover typical bare rock stratum areas such as mountains and mining belts, avoiding invalid pixels with more than 20% vegetation coverage or more than 5% cloud pollution;

[0068] Self-collected datasets need to be labeled. The labeling method is manual labeling, specifically based on field-measured geological maps and mineral sample spectral libraries. Pure rock strata areas are delineated on the images as training samples. The labeling categories cover the main lithologies, such as granite, basalt, and limestone, with no less than 500 pixels labeled for each category.

[0069] In a specific implementation, S2 is as follows:

[0070] Since multispectral remote sensing data is presented in the form of one-dimensional spectral curves, consisting of reflectance values ​​of multiple bands, it has the characteristics of high dimensionality, strong noise interference, and scale inconsistency caused by illumination and atmospheric conditions. Conventional methods often use global min-max normalization to process such data, such as scaling the data to the 0 to 1 range. This method ignores the local absorption characteristics of the spectral curve and its physical meaning, and cannot eliminate the systematic bias caused by atmospheric scattering.

[0071] Therefore, this invention employs an adaptive normalization method based on continuum removal. First, it eliminates the influence of atmospheric scattering and illumination variations by fitting the spectral envelope. Then, it normalizes the local absorption depth based on the difference between the envelope and the minimum reflectance, thereby suppressing noise interference and scale differences while preserving key spectral features crucial for lithological identification. The specific steps are as follows:

[0072] S2.1 Calculation of Continuum Functions:

[0073] For each sample and each band's original reflectance value, a cubic polynomial fitting function is used. The center wavelength of the band and the original reflectance value are input to obtain the fitted value. Then, the maximum value between the original reflectance value and the fitted value is taken as the continuum function value, used to characterize the spectral envelope. This eliminates the influence of atmospheric scattering and illumination variations, and ensures that the envelope is always higher than the original spectrum to avoid underestimating the absorption characteristics. By fitting a cubic polynomial envelope higher than the original spectrum, atmospheric scattering and illumination variations are eliminated, preserving the integrity of the absorption characteristics. This is expressed as:

[0074] ,

[0075] In the formula, To indicate the first The sample at the th The continuum function value for each band, i.e. the calculated value of the spectral envelope in that band, is used to eliminate the effects of atmospheric scattering and changes in illumination. The sample index identifies different remote sensing pixels or observation points; For band indexing, identify different bands in multispectral data; Indicates the first The sample at the th The raw reflectance values of the individual wavebands are the multispectral input data directly observed; represents the central wavelength of the waveband, which is directly read from the remote sensing data metadata and determined by the sensor parameters; represents a polynomial fitting function, specifically a cubic polynomial, which is used to fit the given wavelength and the corresponding raw reflectance value ; represents the order of the polynomial fitting, i.e., the highest degree of the polynomial is 3; represents the max function;

[0076] It should be noted that spectral absorption features are manifested as reflectance dips, and the envelope line needs to cover the upper boundary of these dips to quantify the absorption depth. If the fitted curve is lower than the original data, it will lead to incorrect calculation of the absorption depth, The term ensures that the continuous function values of the fitting are always located at or above the original reflectance values , ensuring that the envelope line is above the original spectrum and avoiding underestimation of absorption features; the polynomial fitting is achieved through the least squares method, aiming to minimize the sum of squared residuals between the observed values and the fitted curve. Specifically, the input wavelength vector and the corresponding reflectance value are used to construct the design matrix . For example, assuming nm, , the design matrix obtained by fitting a cubic polynomial is , then the sparse vector , is the first element of the sparse vector , the second element of the sparse vector , the third element of the sparse vector , the fourth element of the sparse vector , and the output fitted value is , for example, set ; S2.2, absorption depth normalization:

[0077] Using the continuous function value and the minimum reflectance value of the sample in all wavebands, the difference between the original reflectance value and the continuous function value is calculated, and based on the difference between the continuous function value and the minimum reflectance value, the normalized reflectance value is obtained, highlighting the local absorption features and eliminating the global scale difference, represented as:

[0078]

[0079] , ​​

[0080] In the formula, The normalized first The sample at the th The reflectance values ​​of each band represent the spectral characteristics after eliminating global scale differences and highlighting local absorption features; To indicate the first Each sample in all bands The minimum reflectance value in the spectrum is used as the reference point for normalization, which usually corresponds to the deepest absorption feature in the spectrum of that sample. This is a band index, distinct from 'b', used to traverse all bands; To represent the continuum normalization scaling factor, the default value is set to 1. Its function is to control the range of the final normalized value, avoid the absolute value of the normalized result being too large due to the denominator being too small, and thus prevent the dynamic range of the feature from being over-compressed.

[0081] It should be noted that, unlike conventional normalization which uses global range or statistics, this method relies on local absorption characteristics rather than absolute reflectance intensity. The term represents the difference between the continuum value and the minimum reflectance, based on a sample-level minimum reflectance benchmark. Physically equivalent to mapping absorption depth to the scale of the current sample's potential maximum absorption capacity, it enables comparability of absorption intensity across samples. For example, the strong absorption of iron ore in the 700-900 nm range and the weak absorption of gypsum in the 1500 nm range are compressed to similar numerical ranges after normalization, preventing high-reflectance minerals from masking the key absorption of low-reflectance minerals. More importantly, it addresses the differences in background reflectance between different lithologies, such as dark basalt and light sandstone. Conventional fixed-range normalization would cause all band values ​​of dark rocks to approach 0, losing their separability. As a dynamic benchmark, it can automatically adapt to the differences in background reflection of different lithologies.

[0082] In a specific implementation, S3 is as follows:

[0083] The normalized spectral curve still faces the problem of mixed pixels, that is, a single pixel contains mixed reflection signals of multiple rock layers. Conventional methods usually use a linear mixing model to process it. This model assumes that the end-member spectrum is fixed. However, in practical applications, the end-member spectrum will change dynamically due to the variation of mineral composition, resulting in insufficient decomposition accuracy, inability to capture nonlinear mixing effects, and the decomposition results are inconsistent with the absorption characteristics of rocks and minerals.

[0084] Therefore, this invention combines lithological physical mechanisms and adopts a decomposition framework with endmember adaptive updating and abundance sparsity constraints. Through an iterative optimization process, it simultaneously learns endmember spectra and abundance coefficients to ensure that the decomposition results conform to the physical characteristics of rock and mineral absorption. The specific steps are as follows:

[0085] S3.1 Initialization of physical constraint endpoints:

[0086] Based on the normalized reflectance values, sample clusters are generated using clustering methods. For each endmember, its initial reflectance value is calculated as the average of the normalized reflectance values ​​across the key absorption band set. This average is used to characterize the spectral characteristics of pure rock layers, ensuring that the endmembers conform to the true mineral absorption characteristics and avoiding random initialization bias. This is expressed as:

[0087] ,

[0088] In the formula, Indicates the first The terminal in the first The initial reflectance values ​​of each band characterize the spectral features of the pure rock layer; Endmember indexes are used to identify different rock strata types; Indicates the first The sample clusters corresponding to each endmember are all clustered using k-means clustering. of Genesis, where the cluster center corresponds to a physically pure rock layer; This indicates an indicator function that outputs 1 when the condition is met and 0 otherwise. Indicates the first The key absorption band set of rock-like strata is predefined through a priori mineral spectral library to characterize physical constraints; To represent sample clusters The number of samples in the sample;

[0089] It should be noted that endmembers represent the spectral characteristics of pure rock types, indicating the reflectance curve of a single lithology under ideal conditions; key absorption band sets... The mineral spectral library is determined by matching characteristic absorption bands in a known mineral spectral library. This library can be the open-source ASTER (Advanced-Spaceborne-Thermal-Emission-and-Reflection-Radiometer). For example, if the... The rock strata are composed of hematite, and mineral spectral libraries show that they contain hematite. nm has strong absorption, then Includes band index satisfy nm, through physical constraints, ensures that the endmember initialization conforms to the true absorption characteristics of minerals, avoiding randomness; conventional clustering initialization applies equal weight to all bands, but the spectral characteristics of minerals are concentrated only in specific bands, such as epidote in the 2300 nm iron absorption band, through an indicator function. Limit the calculation to a predefined set of key absorption bands. A physical-driven feature filter is constructed to achieve adaptive suppression of noise and redundant bands. In the near-infrared band where there are no feature changes, the algorithm automatically ignores reflectance fluctuations in this region and extracts only spectral segments strongly correlated with mineral composition. Furthermore, because... By encoding domain knowledge into mathematical constraints from a prior knowledge base, endmember initialization is transformed from a purely data-driven approach to a physical-guided hybrid paradigm. When two minerals have similar reflections in non-characteristic bands, conventional methods may generate confusing endmembers, while band-selective averaging can force the separation of mineral-specific spectral morphology.

[0090] S3.2, Sparse Constraint Abundance Estimation

[0091] With fixed endmember reflectance values, the abundance coefficient is calculated for each sample. The abundance coefficient is obtained by minimizing the L2 norm squared of the linear combination of the normalized reflectance value and the endmembers, combined with a sparse regularization term. This reflects the spatial discontinuity of the rock strata and suppresses excessive mixing effects, and is expressed as:

[0092] ,

[0093] In the formula, Indicates the first In the nth sample The abundance coefficient of each endmember, characterizing the proportion of each rock layer in the mixed pixel, must satisfy the following conditions: To avoid negative values ​​leading to non-physical interpretations, and simultaneously satisfy the condition that the sum of the abundance coefficients of all endmembers is 1, i.e. ; This indicates that the function to be minimized is the function whose optimization objective is to minimize the subsequent expression. The total number of endmembers; Indicates the first The terminal in the first Reflectance values ​​for each band; Represents the L2 norm; This represents the sparsity regularization coefficient, used to enhance the sparsity of abundance coefficients, such as... =0.01;

[0094] It should be noted that, since pixels are usually dominated by a few rock layers, The term represents the abundance coefficient. Apply sparsity constraints to force The abundance coefficients of most endmembers are approximately 0 to reflect the spatial discontinuity of rock strata distribution, avoid over-mixing, and improve decomposition accuracy; initial reflectance value It is the initial endmember generated by clustering, and the reflectance value. This is the updated value after iterative optimization; because the abundance needs to satisfy... and Minimizing abundance forces non-dominant endmember coefficients to approach zero, while dominant endmembers are automatically assigned higher weights. The term, as a regularization term, represents the sparsity pressure applied to the abundance and sum. It works synergistically with the non-negative abundance constraint to effectively address the spatial sparsity of rock layer distribution and suppress excessive mixing effects.

[0095] Furthermore, the objective function is solved by optimizing the algorithm. To obtain the optimal abundance coefficient and reflectivity value The optimization algorithm can be gradient descent or non-negative least squares.

[0096] In a specific implementation, S4 is as follows:

[0097] While the abundance coefficients obtained from mixed pixel decomposition reduce the impact of the mixing effect, they lose detailed information from the original spectrum. Conventional methods often simply concatenate the original spectrum and abundance coefficients as features, which can easily lead to feature dimension explosion and redundancy, making it difficult for the model to focus on key bands that are crucial for lithological discrimination. Therefore, this invention adopts a multi-level feature enhancement mechanism. First, the spectral derivative is calculated to highlight the lithologically specific absorption edge features. Then, an attention mechanism is used to weightedly fuse the derivative features with the abundance coefficients, thereby adaptively enhancing lithological features with high discriminative power. The specific steps are as follows:

[0098] S4.1 Calculation of the first derivative of the spectrum:

[0099] For each sample and each band, the difference in normalized reflectance values ​​between adjacent bands is calculated. This difference is then combined with the difference in center wavelength between adjacent bands and the derivative scaling factor to obtain the first derivative value. This first derivative value is used to characterize the local slope of the spectral curve, amplify the rate of change in reflectance between bands, highlight lithologically specific absorption edge characteristics, and suppress smoothing noise. This is expressed as:

[0100] ,

[0101] In the formula, For the first The sample at the th The first derivative values ​​of each band characterize the local slope of the spectral curve and are used to identify lithology-specific absorption valleys. The normalized first The sample at the th The reflectance values ​​of each band represent the spectral characteristics after eliminating global scale differences and highlighting local absorption features; For the first The center wavelength of each band; For the first The center wavelength of each band is directly read from the remote sensing data metadata, determined by the sensor parameters, set , in nanometers; is a derivative scaling factor, used to amplify weak absorption signals while suppressing noise interference, set = 1.2;

[0102] It should be noted that the conventional derivative calculation directly uses the reflectance difference of adjacent bands, ignoring the slope distortion caused by uneven band spacing, The term represents the wavelength difference normalization operation, and the derivative scaling factor forms a coupling relationship, and the wavelength difference is used as the denominator to realize the first-order derivative comparable across bands, and the scaling factor can amplify the discriminative ability of weak absorption edges, and in essence, parameterizes the sensitivity of the model to the rate of spectral change, which can be dynamically adjusted according to the type of mineral, for example, the slow absorption valley of chlorite at 2250nm, the conventional derivative may be overwhelmed by noise, and by setting the hyperparameter, it is significantly highlighted through linear amplification, and because the derivative itself has a small value, the amplified value still remains within a reasonable range, which can form obvious peaks and valleys in the derivative space, and the difference with the surrounding rock is enlarged;

[0103] S4.2, attention weighted feature fusion:

[0104] For each sample, the attention weight of the derivative feature and the attention weight of the abundance coefficient are calculated, and then the weighted derivative feature and the weighted abundance coefficient are summed to obtain the fusion feature value, which is used to adaptively enhance the high discriminative lithology features, and is expressed as:

[0105] ,

[0106] In the formula, is the fusion feature value of the i-th sample in the j-th feature dimension, is the total dimension of the fusion feature, and the calculation method is represented as ; is the total number of bands; is the derivative feature attention weight of the i-th sample in the j-th band, and the calculation method is represented as ; is the natural exponential function; is the attention score function of the derivative feature, which is a learnable linear transformation; is the first-order derivative value of the i-th sample in the j-th band; ​​​​​​is the band index distinguishing from b, used to traverse all bands; is the abundance coefficient attention weight of the i-th endmember in the j-th band for the i-th sample, calculated as ; ; ; is the attention score function of the abundance coefficient, which is a learnable linear transformation; ; ; ; is the endmember index distinguishing from k, used to traverse all endmembers;

[0107] The attention score function of the first derivative feature and the attention score function of the abundance coefficient are both learnable linear transformation functions, taken as an example, the implementation is represented as , wherein, is the weight of the attention score function, is the bias of the attention score function, both are trainable parameters, and the calculated represents the attention score;

[0108] It should be noted that the calculation methods of the attention weight of the derivative feature and the abundance coefficient attention weight are both Softmax normalization calculation methods, which focus on high discriminative bands and strengthen important rock types while ensuring that the sum of weights is 1; conventional splicing or weighted averaging cannot distinguish the physical meaning difference between the derivative feature of local band change and the abundance coefficient of global lithology proportion, and in the calculation process of the feature value , the dual-path attention plays a synergistic technical effect, by calculating and through independent attention mechanisms, a feature type adaptive weight distributor is constructed to realize complementary strengthening of discriminative information, for example, when the derivative of a certain band is significantly caused by the 700nm absorption edge of hematite, its weight is automatically increased, and when the abundance of a certain endmember is prominent caused by a high proportion of granite, its contribution is improved, the dual-path design avoids feature crosstalk, and at the same time focuses on the derivative path of the altered mineral absorption edge and the abundance path of the host rock proportion.

[0109] In the specific implementation, S5 is specifically as follows:

[0110] Because the fused feature values ​​are high-dimensional and highly nonlinear, conventional fully connected networks use standard activation functions such as ReLU for processing. However, such functions are difficult to adapt to the continuous and smooth characteristics of the spectrum, easily leading to gradient vanishing and insufficient feature representation, thus making it difficult to capture the subtle spectral variation patterns required for lithology identification. Therefore, this invention adopts a spectral adaptive activation function and constructs a multi-scale residual block architecture, which can effectively preserve the continuity information of the spectrum in deep networks, while enhancing the feature representation capability of the model. The specific steps are as follows:

[0111] S5.1 Input Layer and Feature Projection:

[0112] The fused feature values ​​are input into the linear layer, and combined with the input projection weight matrix and the input layer bias vector to obtain the initial activation values. This achieves linear compression of the fused feature dimension, retains key lithological information, and reduces subsequent computational complexity, as shown below:

[0113] ,

[0114] In the formula, For the first The sample at the th The initial activation values ​​of the hidden layer nodes; The index of the hidden layer node identifies the different feature dimensions after dimensionality reduction. , The reduced-dimensionality feature representation provides input for subsequent residual blocks; The input projection weight matrix consists of trainable parameters that are learned through training to compress the feature mapping. The input layer bias vector is a trainable parameter used to adjust the feature offset. For the hidden layer dimension, for example, set to To balance expressive power and computational efficiency;

[0115] S5.2 Multi-scale residual block connection:

[0116] For each residual block, the activation value of the previous layer is applied with a spectral adaptive activation function, and then multiplied by the activation values ​​of the previous two layers by the cross-layer skip connection weight matrix to obtain the activation output value of the current layer. This maintains spectral continuity and alleviates gradient vanishing in deep networks, and is expressed as:

[0117] ,

[0118] In the formula, For the first The sample at the th Layer residual block The activation output value of each hidden layer node; For the layer index of the residual block, ; is the total number of residual blocks; is the weight matrix of the -th residual block, is a trainable parameter, performing feature transformation; is the weight matrix of the -th residual block, is a trainable parameter, performing feature transformation; is the output activation value of the -th hidden layer node of the -th residual block for the -th sample; is the output activation value of the -th hidden layer node of the -th residual block for the -th sample; is the bias vector of the -th residual block, is a trainable parameter; is the cross-layer skip connection weight matrix, directly transferring the features of the -th residual block to the -th residual block, is a trainable parameter; is the spectral adaptive activation function, set its input , the calculation method is represented as ; is the input value of the spectral adaptive activation function, corresponding to ; is the hyperbolic tangent function, the output range , providing basic nonlinearity; is the scaling coefficient, used to adjust the strength of the nonlinear term, set ;

[0119] It should be noted that the conventional ReLU activation function faces the contradiction between continuity preservation and mutation detection in spectral data processing. The hard truncation of the ReLU activation function destroys the continuity, and the Softplus activation function is too smooth, provides basic nonlinearity, and the output range is , represents the Gaussian decay term, and the decay speed is controlled by , makes the function approximately linear in the spectral smooth area and enhances the nonlinearity in the feature mutation area, which is consistent with the spectral continuous change characteristics. For example, when is small, it represents the spectral smooth area, , the spectral adaptive activation function is approximately linear, preserving the spectral continuity, and when is large, it represents the feature mutation area, , The term plays a dominant role, enhancing nonlinearity to capture absorption edges. Through adaptive smoothing-nonlinear switching, it improves sensitivity to subtle spectral changes. For example, in the quartz-feldspar mixed spectrum, the gradient of the weak absorption valley of feldspar may vanish under the ReLU activation function, while the spectral adaptive activation function maintains near-linear propagation with small inputs, enabling effective backpropagation of the gradient. Simultaneously, the parameters... and By controlling the nonlinear intensity and decay rate respectively, an adjustable spectral response surface can be formed;

[0120] S5.3, Probability Prediction of Rock Strata Types:

[0121] The activation values ​​of the final residual block are combined with the output layer weight matrix and the output layer bias vector, and the softmax function is applied to calculate the predicted probability of the rock strata category to support the classification decision, as shown below:

[0122] ,

[0123] In the formula, For the first The sample belongs to the first The predicted probabilities of each category, ; This represents the total number of rock strata categories. The softmax function converts the input vector into a probability distribution. is the output layer weight matrix, which consists of trainable parameters; The output layer bias vector is a trainable parameter.

[0124] In a specific implementation, S6 is as follows:

[0125] Conventional methods only use the cross-entropy loss function, focusing solely on classification accuracy and neglecting the physical consistency constraints that the spectral decomposition process should satisfy, such as the sparsity of abundance coefficients and the minimization of endmember reconstruction errors. This can easily lead to a disconnect between the model's prediction results and the physical mechanism of lithology, resulting in a separation between decomposition and prediction, and reducing the physical interpretability of rock strata identification. Therefore, this invention incorporates the physical constraints of spectral decomposition into the loss function, and forces the model's prediction results to conform to the physical mechanism of rock and mineral mixing by jointly optimizing the classification loss term and the physical regularization term.

[0126] S6.1.1 Constructing a framework for the joint loss function:

[0127] The total loss function is obtained by combining the cross-entropy classification loss term and the physical regularization term. The classification loss and physical regularization term are jointly optimized to balance model accuracy and the physical consistency of spectral decomposition, expressed as:

[0128] ,

[0129] where, is the total loss function, guiding the model end-to-end training; is the cross-entropy classification loss term, quantifying the deviation between predicted lithology class and true label; is the physical regularization term, constraining the network output to comply with the physical mechanism of spectral decomposition; is the physical regularization weight coefficient, controlling the strength of physical constraint, set as = 0.05;

[0130] S6.1.2, Cross-entropy classification loss calculation:

[0131] For all training samples and all classes, the average of the logarithmic product of true label and predicted probability is calculated to obtain the cross-entropy loss term, which quantifies the deviation between predicted probability and true label, drives the model to accurately identify the lithology class, and is expressed as:

[0132] ,

[0133] where, is the total number of training samples; is the total number of lithology classes; is the One-hot encoding of the th sample belonging to the th class, which is 1 when the sample belongs to the th class, otherwise 0;

[0134] S6.1.3, Physical regularization term construction:

[0135] For each sample, the L2 norm square of the linear combination of the normalized reflectance value and the predicted abundance coefficient and the endmember reflectance value is calculated, and the absolute value of the abundance coefficient is combined to obtain the physical regularization term, which constrains the network prediction to comply with the physical mechanism of mineral mixing, and is expressed as:

[0136] ,

[0137] where, is the abundance coefficient of the th endmember in the th sample predicted by the network through the auxiliary layer; is the absolute value function, imposing the sparsity constraint on the abundance coefficient;

[0138] The auxiliary layer can adopt a fully connected layer, which inputs the fused feature value and outputs the predicted abundance coefficient During training, the physical regularization term supervises the abundance coefficient ​The generation process is forced to approximate the true abundance distribution, and a Softmax function is applied after the auxiliary layer to satisfy... ;

[0139] It should be noted that, through reflectivity values... With abundance coefficient The hidden layer encoding of the forced classification network yields physically interpretable spectral decomposition results. The term characterizes the reconstruction error, requiring the network to predict abundance energy that can accurately reconstruct the input spectrum. Ensuring that the predicted abundance energy can accurately reconstruct the original spectrum forces physical consistency decomposition, which is equivalent to adding a self-supervised signal. The term represents the abundance sparsity. Since a few rock layers dominate in a pixel, the abundance sparsity can reflect the sparsity of rock layer distribution. During model training, the reflectance value... As a learnable parameter, the mineral mixing law is backpropagated through the physical regularization term, so that the hidden features of the network carry lithological physical interpretation. For example, when the classifier confuses gabbro and diorite, it can punish the reconstruction deviation of the endmember spectra of the two, driving the network to learn more discriminative mineral absorption feature expressions.

[0140] When training models with high-dimensional spectral data, it is easy to get stuck in local optima. Conventional Adam optimizers use a fixed learning rate update strategy for all parameters, which cannot distinguish the different roles of weight parameters and bias parameters in the model and their update requirements. This leads to unstable convergence of endmember-related parameters and reduces the model's ability to capture subtle differences in lithology. Therefore, this invention adopts a parameter grouping adaptive learning rate mechanism, which groups parameters according to their physical roles in the model. A dynamically adjusted cosine annealing learning rate is applied to the weight parameters that control feature transformation to stabilize the convergence of deep networks. At the same time, a fixed small learning rate is applied to the bias parameters that control the offset to suppress overfitting.

[0141] S6.2.1 Define parameter update rules:

[0142] Using the Adam optimizer, the trainable parameters are updated based on the bias correction terms from the first-order moment estimates, the bias correction terms from the second-order moment estimates, and the group-adaptive learning rate, as follows:

[0143] ,

[0144] In the formula, For the first The trainable parameters for the next iteration are in set form; For the first The trainable parameters for the next iteration are in set form; For the first The bias correction term for the first-order moment estimate in the next iteration is calculated as follows: ; the gradient first moment of the first iteration; the exponential decay rate of the gradient first moment; the exponential decay rate of the gradient first moment; the second moment of the first iteration; the bias correction term of the second moment estimate of the first iteration, calculated as ; the gradient second moment of the first iteration; the exponential decay rate of the gradient second moment; the exponential decay rate of the gradient second moment; the numerical stability constant to prevent denominator zero, set as ; the grouped adaptive learning rate; It should be noted that the exponential decay rate of the gradient first moment and the exponential decay rate of the gradient second moment

[0145] are hyperparameters of the Adam optimizer, used to eliminate the bias of the moment estimate in the early stage of optimization, and are usually set as , ;

[0146] S6.2.2, Grouped Learning Rate Setting:

[0147] According to the parameter physical effect, two groups are divided, and for the weight parameter group, the basic learning rate is multiplied by the cosine function value, and the cosine annealing is implemented to accelerate convergence. For the bias parameter group, the basic learning rate is multiplied by a fixed value, and the fixed small learning rate is used to suppress overfitting, which is expressed as:

[0148] ,

[0149] In the formula, is the trainable weight parameter group, which controls the feature transformation and information transmission; is the trainable bias parameter group, which controls the decision boundary offset; is the basic learning rate, which is set as =0.001; is the current training iteration number; is the maximum training iteration number; is the cosine function, the item implements the cosine annealing strategy of the weight parameter; is the trainable parameter, which is in set form; is the circular constant, calculated as ;

[0150] It should be noted that in the early stage of training, , , the weight convergence is accelerated; in the later stage of training,​​ , , fine-tuning parameters to stabilize the model; while the bias parameters are fixed small learning rate, , avoid the prediction bias caused by overfitting;

[0151] S6.3, the model training process is as follows:

[0152] The collected data is forward propagated, and after processing, it is input into the rock layer identification model to obtain the prediction probability of the rock layer category, then the loss function is calculated, and then the back propagation is performed, the model parameters are updated, and the model is iteratively optimized by the Adam optimizer, until the iteration termination condition is met, the optimization is ended, and the trained rock layer identification model is obtained;

[0153] The iteration termination condition is as follows: after each round of training, the validation set accuracy is evaluated, if the validation loss does not decrease for 10 consecutive rounds or the accuracy fluctuation is less than 0.5%, the training is terminated and the optimal model parameters are saved, to avoid overfitting, the maximum number of iterations is set to 500 rounds, and the basic learning rate is fixed at 0.001.

[0154] In the specific implementation, S7 is specifically as follows:

[0155] Input new multispectral remote sensing image, perform continuum normalization to eliminate illumination / atmospheric bias pixel by pixel, then perform physical constraint mixed pixel decomposition, output end member abundance coefficient, then perform attention weighted feature fusion, fuse derivative feature and abundance coefficient; the fused features are input into the trained multi-scale residual network, and linear projection, residual block nonlinear transformation are performed, and the lithology category probability distribution of each pixel is output, and the lithology category corresponding to the maximum probability value is taken as the recognition result of the pixel.

[0156] Embodiment 2

[0157] A rock layer automatic identification system based on multispectral remote sensing, comprising a module for executing processing instructions of each step in the rock layer automatic identification method based on multispectral remote sensing, comprising the following modules:

[0158] Data collection module: collect multispectral remote sensing data samples;

[0159] Data preprocessing module: normalize, mixed pixel decomposition and special fusion operation are performed on the collected samples;

[0160] Rock layer identification module: input the preprocessed samples into the rock layer identification module to predict the lithology of the samples;

[0161] Optimization training module: optimize and update the rock layer identification module parameters through the model total loss function and the adaptive learning rate mechanism, and obtain the trained rock layer identification module;

[0162] Recognition result output module: After the sample to be identified is preprocessed, it is input into the trained rock stratum recognition module and the final recognition result is output.

[0163] Example 3

[0164] like Figures 2-3 The figures shown are schematic diagrams of the reflectance of the original spectral curve and the envelope of the continuum, and the reflectance of the spectrum after the continuum is normalized. Figure 2 The original reflectance characteristics of different rock layers (granite, basalt, limestone, gneiss, and sandstone) in Landsat-8-OLI multispectral bands (seven bands from Band 1 to Band 7) are shown. The dashed lines represent the continuum envelope fitted by a cubic polynomial. This envelope is always above the original spectrum to accurately capture spectral absorption characteristics. This processing can effectively eliminate the influence of atmospheric scattering and illumination changes, and provide a physically accurate benchmark for subsequent normalization.

[0165] Figure 3 The spectral curves are presented after adaptive normalization, which is based on the difference between the continuum envelope and the minimum reflectance. This significantly highlights the key absorption characteristics of different rock layers, such as the strong absorption valley of basalt in Band 5. At the same time, it eliminates the scale inconsistency caused by the difference in background reflectance of different lithologies, making the characteristics of rock layers with large reflectance differences, such as dark basalt and light sandstone, comparable.

[0166] like Figure 4 As shown, a comparative analysis of normalization methods was conducted to verify the effect of continuum normalization on improving the accuracy of rock strata identification. The horizontal axis represents five typical rock strata: granite, basalt, limestone, gneiss, and sandstone, while the vertical axis represents the classification accuracy (dimensionless, numerical range 0-1). Figure 4 The medium gray bars represent the conventional min-max normalization method, while the blue bars represent the continuum normalization method of this invention. Experimental data show that the conventional min-max normalization method, by ignoring local absorption characteristics, results in significantly lower accuracy in identifying basalt (dark rocks) and limestone (low-reflectance minerals), proving its inability to adapt to the differences in background reflectance among different lithologies. The continuum normalization method of this invention achieves higher accuracy than the conventional method across all five types of rock formations. Because it employs a dynamic benchmark and spectral envelope fitting, it solves the problem of dark rock characteristics approaching zero, while preserving the key absorption characteristics of silicate and carbonate minerals. Figure 4 The uniform and stable height of the blue columns in the middle verifies that the continuum normalization achieves cross-lithological comparability through wavelength difference standardization.

[0167] Example 4

[0168] like Figures 5-6 As shown, Figure 5Three pure rock layer endmember spectra were demonstrated by clustering and key absorption band selection. These endmembers were calculated based on the initial reflectance of the predefined key absorption bands in the mineral spectral library (such as the silicate absorption characteristics of granite in Band 4-5), avoiding random initialization bias and ensuring that the endmembers conform to the true mineral absorption characteristics, thus providing physical mechanism constraints for the decomposition of mixed pixels.

[0169] Figure 6 The bar chart visually displays the composition ratio of different rock layers in a single pixel. The abundance coefficient is obtained through sparsity constraint optimization, reflecting the spatial discontinuity of rock layer distribution. The dominant rock layer (such as granite) has a high abundance value, while the abundance of non-dominant rock layers is close to zero, which can effectively suppress the excessive mixing effect.

[0170] Example 5

[0171] like Figures 7-8 As shown, Figure 7 The study revealed the degree of attention the model pays to the derivative features of different bands. Higher weights are concentrated in key bands that are sensitive to lithology, such as Band4-5 (640-880nm). The adaptive attention mechanism can enhance discriminative features such as the 700nm absorption edge of hematite and suppress noise interference in non-characteristic bands.

[0172] Figure 8 The model demonstrates its differentiated focus on the abundance characteristics of different rock strata, with limestone receiving the highest weight as the dominant stratum. This attention allocation allows the model to focus on spatially significant rock strata types and enhances the spectral contribution of the dominant strata during feature fusion.

[0173] like Figure 9 As shown, the feature enhancement mechanism of this invention is used to analyze the separation ability of different rock layers through two-dimensional feature space visualization. Figure 9 In the scatter plot shown, the horizontal axis represents "shortwave infrared characteristics (band 5-7 derivatives)" (dimensionless), and the vertical axis represents "abundance characteristic weighted values" (dimensionless). Figure 9 The five types of rock strata show obvious clustering distribution. Granite (red dot group) is concentrated in the upper right area, basalt (blue dot group) is clustered in the lower left area, limestone (green dot group) is distributed in the middle, gneiss (purple dot group) is located in the upper middle, and sandstone (cyan dot group) occupies the lower right.

[0174] The experimental results show that: 1) Similar rock strata clusters are closely clustered, proving that the feature representation has high consistency; 2) Dissimilar rock strata have clear boundaries, and the decision boundary can smoothly separate the boundaries of different categories; 3) The feature spatial structure conforms to geological laws, such as gneiss (metamorphic rock) being located in the transition zone between granite (igneous rock) and limestone (sedimentary rock).

[0175] The experimental results also show that about 10% of the granite samples intrude into gneiss area, which is consistent with the mixed rock phenomenon in the field measurement; and about 5% of the limestone samples overlap with sandstone, which reflects the similarity of the carbonate mineral composition of the two;

[0176] The experimental results verify the effectiveness of the attention-weighted feature fusion, and the discriminability of the original spectrum is improved by strengthening the derivative feature and the abundance feature.

[0177] The above describes the specific embodiments of the application in combination with the drawings, but is not a limitation on the protection scope of the application. Various modifications or variations made by those skilled in the art on the basis of the technical solutions of the application without creative labor are still within the protection scope of the application.

Claims

1. A method for automatically identifying rock strata based on multispectral remote sensing, characterized in that, The application comprises the following steps: An adaptive normalization method based on continuum removal is used on the collected multispectral remote sensing data samples to obtain normalized samples; Based on the physical mechanism of the rock mass, the normalized samples are decomposed, each endmember after decomposition is initialized, and the abundance coefficient of each endmember is calculated; The first derivative value of each normalized sample in each waveband is calculated, and the abundance coefficient is combined for attention weighted feature fusion to obtain multi-dimensional fusion feature values; a rock stratum recognition model based on deep learning is constructed, and the spectral adaptive activation function and multi-scale residual block architecture in the model are used to predict the rock stratum class probability of the fusion feature values in each dimension; The rock stratum recognition model is optimized by integrating the physical constraints of spectral decomposition into the loss function, and the model parameters are updated using the adaptive learning rate mechanism of parameter grouping to obtain the trained rock stratum recognition model; the newly collected multispectral remote sensing data is processed and input into the trained model to obtain the final rock stratum class probability prediction result, and the lithology class corresponding to the maximum probability value is taken as the recognition result.

2. The method for automatic lithological identification based on multispectral remote sensing according to claim 1, characterized in that, The adaptive normalization method based on continuum removal is as follows: For the original reflectance of each normalized sample in each waveband, first, the influence of atmospheric scattering and light change is eliminated by fitting the spectral envelope line, specifically, a cubic polynomial fitting function is used, the center wavelength and the original reflectance value of the input waveband are input, the fitting value is obtained, and then the maximum value of the original reflectance value and the fitting value is taken as the continuum function value of each sample in each waveband, which is used to represent the spectral envelope line; Then, the difference between the continuum function value and the minimum reflectance is used to normalize the local absorption depth to obtain the reflectance value of each normalized sample in each waveband, i.e. the normalized sample after normalization processing.

3. The method according to claim 2, characterized in that, The abundance coefficient calculation process of each endmember is as follows: A clustering method is used to generate sample clusters for each endmember after decomposition, based on the normalized samples, and an indicator function is used to determine whether the current waveband matches the known mineral spectral library, and the initial reflectance value of each endmember in each waveband is calculated; Then, the initial reflectance value of the endmember is fixed, and the abundance coefficient of each sample in each endmember is calculated by minimizing the L2 norm square of the linear combination of the normalized reflectance value and the coefficient regular term; Where endmember represents different rock stratum types.

4. The method according to claim 3, characterized in that, The calculation process of the multi-dimensional fusion feature value is as follows: For each waveband of each sample, the difference between the normalized reflectance value and the adjacent waveband is calculated, and the first derivative value is obtained by combining the center wavelength difference of the adjacent waveband and the derivative scaling factor; Then, for each sample, the attention weight of the first derivative value in each dimension and the attention weight of the abundance coefficient are calculated, and then the weighted first derivative feature value and the weighted abundance coefficient are summed to obtain the fusion feature value in each dimension; Where the dimensions include the total number of wavebands and the total number of endmembers.

5. The method according to claim 4, characterized in that, The operation in the rock stratum recognition model based on deep learning is as follows: The fusion feature values of each dimension are input into the model, the fusion feature values of each dimension are input into a linear layer, and the initial activation values of each sample at each hidden layer node are obtained by combining the projection weight matrix of the input layer and the bias vector of the input layer. The initial activation values are input into the multi-scale residual block architecture. Specifically, the activation function values of the previous layer residual block are combined with the activation function values of the first two layers of the residual block multiplied by the cross-layer skip connection weight matrix after passing through the spectral self-adaptive activation function, to obtain the activation function values of the current layer residual block. After layer-by-layer calculation, the activation function values of the final layer residual block are obtained. The activation function values are combined with the output layer weight matrix and the output layer bias vector, and the prediction probability of the rock layer class is obtained by calculating the Softmax function.

6. The method according to claim 5, characterized in that, The physical constraints of spectral decomposition are integrated into the loss function, and the total loss is calculated by combining the cross-entropy classification loss term and the physical regularization term. The cross-entropy loss term is calculated by quantifying the deviation between the prediction probability of the rock layer class calculated by the model and the true class. For each sample, the L2 norm square of the linear combination of the normalized reflectance value, the predicted abundance coefficient, and the endmember reflectance value is calculated, and the physical regularization term is obtained by combining the absolute value of the abundance coefficient.

7. The method according to claim 6, characterized in that, A parameter grouping adaptive learning rate mechanism is used to group the trainable parameters according to their physical effects in the model. A dynamically adjusted cosine annealing learning rate is applied to the weight parameters that control the feature transformation to stabilize the convergence of the deep network, and a fixed small learning rate is applied to the bias parameters that control the offset to suppress overfitting.

8. The method according to claim 7, characterized in that, The training process of the model is as follows: The collected data is forward propagated, processed, and input into the rock layer identification model to obtain the prediction probability of the rock layer class. Then the loss function is calculated, and the model parameters are updated through back propagation. The model is iteratively optimized by the Adam optimizer until the iteration termination condition is met, the optimization is ended, and the trained rock layer identification model is obtained.

9. The method according to claim 1, wherein: The collected multi-spectral remote sensing data samples are derived from open source data sets or self-collected data sets; The collected multi-spectral remote sensing data samples are labeled with lithology class labels.

10. A rock stratum automatic identification system based on multispectral remote sensing, characterized in that: The method comprises the following modules: Data collection module: collect multi-spectral remote sensing data samples; Data preprocessing module: normalize, mixed pixel decomposition, and special fusion operation on the collected samples; Rock layer identification module: input the preprocessed samples into the rock layer identification module to predict the lithology of the samples; Optimization training module: optimize and update the parameters of the rock layer identification module through the model total loss function and the adaptive learning rate mechanism to obtain the trained rock layer identification module; Identification result output module: input the preprocessed samples into the trained rock layer identification module to output the final identification result.

Images