A mineral index calculation method and system based on global multispectral remote sensing imagery

Abstract

The present invention discloses a mineral index calculation method and system based on global multispectral remote sensing images, which belongs to the field of image data processing technology. The method includes: obtaining original remote sensing image data and parsing metadata; based on the metadata, screening the original remote sensing image data and eliminating invalid data; resampling and slicing the remote sensing image data; performing data preprocessing, converting the resampled and sliced data into data that can directly calculate the mineral index; performing image mosaicking and image color matching, integrating different remote sensing image slices in the same geographical area into a single, continuous image, and reducing the color difference between slices to obtain a global map; calculating the mineral index indicator based on the global map. The present invention realizes the whole process of automatic data screening, data preprocessing, remote sensing mosaicking and color matching to the final mineral index calculation, avoiding the storage, memory and other pressures faced by traditional software; and also uses a lookup table method to greatly save computing time.

Description

Technical Field

[0001] The present invention relates to the field of image data processing technology, and in particular to a mineral index calculation method and system based on global multispectral remote sensing images. Background Art

[0002] Altered rocks are rocks whose mineral composition, chemical composition, structure, and texture have changed under the influence of hydrothermal processes. Because they are often found around hydrothermal deposits, they are called altered host rocks and are an important indicator for mineral exploration. Based on the differences in spectral characteristics between mineralized altered rocks and host rocks, image enhancement methods can be used to obtain image variables that enhance the mineralization and alteration information, ultimately achieving the goal of extracting mineralization and alteration information. Multispectral remote sensing data, with its rich spectral information and global coverage, has important applications in geological exploration and resource management. However, the processing and computation of multispectral and multi-time series remote sensing data on a global scale still faces many technical challenges. Remote sensing data is typically collected in regions and time periods. Due to differences in imaging conditions and time periods across different regions, these data exhibit inconsistencies in radiometric properties and geometric positions. Furthermore, achieving seamless global stitching requires solving the problems of cloud removal, mosaic matching, and color uniformity for cross-regional and multi-time series data, placing extremely high demands on algorithm design and data processing capabilities. In addition, geoscience researchers often face the problem of limited computing resources, making it difficult to efficiently process large-scale data and complete complex preprocessing tasks.

[0003] Traditional calculations of mineral alteration indices typically focus on a small, localized area, such as a mineralization belt or mineral point. This involves screening and extracting corresponding remote sensing images, then using remote sensing processing software to perform multiple processing steps, which is time-consuming. This approach is unsuitable for processing and calculating large-scale remote sensing images for several reasons. First, manual data screening is not possible for large-scale remote sensing images. Second, conventional remote sensing processing software is incapable of calculating large-scale remote sensing data, due to limitations in time, efficiency, and storage space. Finally, conventional software cannot fuse and calculate remote sensing images collected over different time periods and regions. Summary of the Invention

[0004] In view of the shortcomings of the traditional method of processing remote sensing images using software on a scene-by-scene basis, the purpose of the present invention is to provide a mineral index calculation method and system based on global multispectral remote sensing images.

[0005] The object of the present invention is achieved through the following technical solution: a mineral index calculation method based on global multispectral remote sensing imagery, comprising the following steps:

[0006] Acquire original remote sensing image data and parse metadata; the original remote sensing image data includes visible light data and shortwave infrared data; the metadata includes extraction time, cloud content, geographic coordinates, and quality parameters;

[0007] Based on the metadata, the original remote sensing image data is screened to remove invalid data; for example, when the quality score is 1 or the cloud content is greater than 40%, data is removed;

[0008] Resample and slice the filtered remote sensing image data, store the image slices in tiff format, and write metadata into tiff files to prevent data information loss;

[0009] Perform data preprocessing to convert the resampled and sliced data into data that can directly calculate the mineral index;

[0010] Perform image mosaicking to integrate different remote sensing image slices within the same geographic area into a single, continuous image;

[0011] Perform image color grading to reduce color differences between slices and obtain a global image;

[0012] Based on the global map, mineral index indicators are calculated.

[0013] Furthermore, resampling and slicing the filtered original remote sensing image data includes:

[0014] The original remote sensing image data is resampled, coordinate converted and sliced, specifically: the target resolution is set, and the coordinate systems of the visible light data and shortwave infrared data are converted to a unified coordinate system; the target resolution is set, and the shortwave infrared band is resampled to be consistent with the visible light band so that it is spatially aligned with the visible light band; each scene image is cut into 1024*1024 image slices, stored in tiff format, and the metadata file is written into the tiff file to prevent the spectral information and spatial information of the data from being lost or distorted, so as to facilitate subsequent data management and query; wherein, the resampling uses interpolation algorithms such as bilinear interpolation and cubic convolution interpolation.

[0015] Furthermore, the original remote sensing image data is ASTER multispectral remote sensing data.

[0016] Furthermore, the data preprocessing includes radiation calibration and atmospheric correction, specifically:

[0017] Input the raw data quantization value (DN) and convert it into the radiance value (Radiance) through radiation correction. Based on the calibration coefficient, perform absolute calibration calculation on the data of each band to eliminate the influence of the sensor's own characteristics on the data, so that the data accurately reflects the radiation characteristics of the ground object.

[0018] The 6S model lookup table is used to perform atmospheric correction on remote sensing images to eliminate the influence of atmosphere on surface reflectivity.

[0019] Furthermore, the atmospheric correction of the remote sensing image using the 6S model lookup table includes:

[0020] A lookup table of atmospheric correction coefficients for all parameters is constructed. The parameters input by the 6s model are discretized and sampled. These parameters include the solar zenith angle (SZA), solar azimuth angle (SAA), aerosol optical depth (AOD), digital elevation model (DEM), atmospheric profile, band, and date information. After determining the sampling points, the 6s (Second Simulation of a Satellite Signal in the Solar Spectrum) model is called to calculate the correction coefficients in the lookup table. These correction coefficients are then stored in a lookup table indexed by parameter combination.

[0021] For image slices that need to be sliced, extract the solar zenith angle, solar azimuth, aerosol optical depth, digital elevation model, atmospheric profile, band and date information from their metadata;

[0022] Based on the discrete data points in the atmospheric correction coefficient lookup table, the parameter combination of the 6S model input is interpolated using a linear interpolation method to obtain the correction coefficient; at the same time, the linear interpolation method is used to correct the error caused by the input date;

[0023] The image slices are corrected using the atmospheric correction coefficients.

[0024] Furthermore, the image mosaicking comprises:

[0025] Preprocess each image slice of multiple bands and remove negative values to ensure the validity and consistency of the data;

[0026] For each slice data on the grid, the data is sorted according to its coverage and cloud content on the slice;

[0027] Perform minimum value slice synthesis based on shadow suppression; for a certain grid, fill data pixel by pixel, calculate the minimum value and the second minimum value of each pixel, if the difference between the two is greater than 0.2, it is determined to be a shadow area, and weighted synthesis is used for the shadow area; for non-shadow areas, directly fill in the minimum value of that position in all slices.

[0028] Furthermore, for each grid slice data, the sorting is performed according to the coverage and cloud content of the data on the slice:

[0029] Cloud detection is performed on each slice through spectral feature analysis and threshold segmentation techniques; the spectral features of the image segments are calculated; the spectral features are automatically threshold segmented using the Otsu method to determine the boundary between cloud pixels and non-cloud pixels; a cloud mask is generated, in which cloud pixels are marked as True and non-cloud pixels are marked as False; based on the cloud mask, the cloud cover ratio is calculated; and all slices are sorted according to the two dimensions of data coverage and cloud content.

[0030] Furthermore, performing image color uniformity includes:

[0031] Using MODIS surface reflectance data, after mosaicking and band normalization, a globally unified uniform color basemap is output;

[0032] Divide the reference image and the data to be uniformly colored into overlapping small blocks, obtain the tasseled cap transformation or equal-weight transformation coefficients, and calculate the brightness characteristics;

[0033] Perform local block segmentation to resolve overall uniform brightness deviation caused by cloud cover or cloud shadows;

[0034] Adjust the spectrum based on brightness to ensure spectral consistency.

[0035] Furthermore, the calculation of mineral index indicators based on the global map includes:

[0036] The value of each mineral index is calculated using the band ratio method;

[0037] Obtaining ground object mask data and removing the mask area data when calculating the mineral index image; the ground object mask data includes vegetation, water, ice and snow;

[0038] For each calculated index, the maximum and minimum values of the numerical part are calculated, the mineral index is normalized to [0, 255], and stored in PNG format for export.

[0039] The present invention also provides a mineral index calculation system based on global multispectral remote sensing images, comprising:

[0040] Data acquisition module, used for downloading remote sensing data and auxiliary data sets and extracting metadata;

[0041] Data screening module, used to automatically screen remote sensing data;

[0042] Resampling and slicing module, used for resampling remote sensing images and converting slice formats;

[0043] Remote sensing image preprocessing module, used for data preprocessing of remote sensing data, including radiometric calibration and atmospheric correction processing;

[0044] Remote sensing mosaic module, used to stitch multiple slices;

[0045] Color uniformity module, used to uniformize remote sensing image slices;

[0046] Calculation module, used for mineral index calculation;

[0047] Result storage module, used for data conversion and storage.

[0048] The beneficial effects of the present invention are as follows: It proposes a computational framework suitable for global remote sensing image processing, transforming the traditional scene-based remote sensing data processing into slice-based remote sensing data processing, thus achieving a complete process from automatic data screening, remote sensing mosaicking and color grading, to final mineral index calculation. Furthermore, this computational framework divides processing tasks into computational sub-units that can be placed in a parallel computing environment, avoiding the storage and memory pressures faced by traditional software processing. In the time-consuming atmospheric correction step, the online algorithm call is converted to a lookup table, which avoids repeated calls to the atmospheric correction algorithm on large-scale slice data and greatly saves computing time. BRIEF DESCRIPTION OF THE DRAWINGS

[0049] Figure 1 is a flow chart of the method of the present invention;

[0050] Figure 2 It is a schematic diagram of the atmospheric correction of the present invention;

[0051] Figure 3 It is the remote sensing mosaic flow chart of the present invention;

[0052] Figure 4 Schematic diagram of the remote sensing image processing results of the present invention. DETAILED DESCRIPTION

[0053] Exemplary embodiments will be described in detail herein, examples of which are illustrated in the accompanying drawings. In the following description, when referring to the drawings, like numbers in different figures represent like or similar elements unless otherwise indicated. The embodiments described in the following exemplary embodiments are not intended to represent all possible embodiments consistent with the present invention. Rather, they are merely examples of apparatus and methods consistent with certain aspects of the present invention, as detailed in the appended claims.

[0054] The terms used in this invention are for the purpose of describing specific embodiments only and are not intended to limit the invention. The singular forms "a," "the," and "the" used in this invention and the appended claims are also intended to include plural forms unless the context clearly indicates otherwise. It should also be understood that the term "and / or" as used herein refers to and includes any or all possible combinations of one or more of the associated listed items.

[0055] like Figure 1 As shown, an embodiment of the present invention provides a method for calculating mineral index based on global multispectral remote sensing images, comprising the following steps:

[0056] Step 1. Data download and metadata extraction.

[0057] 1.1. Download ASTER multispectral remote sensing data as the raw processed remote sensing image data. Select ASTER bands 1 to 9, including visible light and shortwave infrared data. Download the corresponding ASTER GDEM elevation data and the corresponding time-domain MODIS AOD data as auxiliary data for atmospheric correction. Download vegetation index and water body data as mask data.

[0058] 1.2. Parse the metadata in the XML of ASTER remote sensing data and extract information such as time, cloud content, geographic coordinates, and quality parameters.

[0059] Step 2. Data screening.

[0060] Perform a preliminary screening of the downloaded raw remote sensing data. Set screening rules based on cloud content and image anomalies. Eliminate data when the quality score is 1 or the cloud content is 40%. Based on the extracted metadata, determine data availability according to the screening rules and mark qualified data.

[0061] Step 3. Resample and slice the data.

[0062] Because remote sensing images collected at different times often overlap and intersect, the original remote sensing images were resampled, coordinate-converted, and sliced. The original visible light band had a resolution of 30 meters, and the shortwave infrared band had a resolution of 100 meters. With a target resolution of 30 meters, the coordinate systems of the visible and shortwave infrared image data were converted to a unified EPSG:3857 coordinate system. Then, with a target resolution of 30 meters, the original visible light band (30-meter resolution) and shortwave infrared band (100-meter resolution) were resampled to 30 meters. Interpolation algorithms such as bilinear interpolation and cubic convolution were used for resampling. Finally, each scene image was sliced into 1024*1024 image slices and stored in TIFF format. During the format conversion process, to ensure that key data elements such as spectral and spatial information were not lost or distorted, metadata files were written to the TIFF files for subsequent data management and querying.

[0063] Step 4. Data preprocessing.

[0064] Data preprocessing converts the raw data quantization value (DN) collected by remote sensing sensors into data that can be used to directly calculate the mineral index after removing the influence of factors such as sensors and atmospheric aerosols.

[0065] 4.1. Radiometric Calibration. The original remote sensing digital quantization value (DN) is input and converted to a radiometric brightness value (Radiance) through radiometric correction. Based on the calibration coefficients provided by the USGS, precise absolute calibration calculations are performed on the data for each band to eliminate the influence of the sensor's inherent characteristics on the data, ensuring that the data accurately reflects the radiometric characteristics of the ground object. The radiometric calibration calculation is shown in the following formula:

[0066] .

[0067] Where m is the gain coefficient. The gain mode of the band is obtained from the information recorded in the image header data, and the specific gain coefficient is obtained from the gain matrix in the official USGS document. The default value of n is 0.

[0068] 4.2. Atmospheric correction. Atmospheric correction simulates the scattering and absorption of light by the atmosphere, eliminates the influence of the atmosphere on the surface reflectivity, realizes the precise correction of remote sensing data, and obtains more accurate surface reflectivity information. The present invention uses a 6s model lookup table to perform atmospheric correction on remote sensing images. The traditional 6s model queries the aerosol optical depth, elevation data and metadata information of the remote sensing image corresponding to the remote sensing image, inputs the correction coefficient into the model, and transforms the image. This process must be performed on each slice, which is very time-consuming for large-scale data calculations. Therefore, the present invention uses a table lookup to perform 6s model atmospheric correction.

[0069] 4.2.1. Construct the atmospheric correction coefficient lookup table of all parameters, such as Figure 2 For the 6s model inputs, including the solar zenith angle (SZA), solar azimuth angle (SAA), aerosol optical depth (AOD), digital elevation model (DEM), atmospheric profile, band, and date information, parameters were discretized and sampled, as shown in Table 1. The sampling point distribution was designed based on the parameter range and computational efficiency, with dense sampling in low-value areas and sparse sampling in high-value areas to balance lookup table accuracy and storage overhead.

[0070] parameter Sampling point distribution Solar zenith angle (SZA) [0,2,4,6]+list(range(7,42,8))+list(range(43,60,4))+list(range(61,70,2))+list(range(71,81,1)) Aerosol Optical Depth (AOD) [2, 4, 6, 10, 20, 40, 70, 100, 150]+list(range(201,801,100))+list(range(801,2002,400)) Elevation (DEM) [0，200，1000，2000，4000，8000] Atmospheric profile type [‘1’，‘2’，‘3’，‘4’，‘5’] Band number [1，2，3，4，5，6，7，8，9] date [datetime(2020,1,1),datetime(2020,7,1)] Solar Azimuth (SAA) Fixed at 90 degrees

[0071] After determining the sampling points, the 6s (Second Simulation of a Satellite Signal in the Solar Spectrum) model is called to calculate the correction coefficients in the lookup table. This method uses the Python packaged Py6S library to call the 6s model for calculations. The resulting correction coefficients are stored in a lookup table, indexed by parameter combination.

[0072] 4.2.2. For the image slices to be processed, extract the solar zenith angle (SZA), solar azimuth angle (SAA), aerosol optical depth (AOD), digital elevation model (DEM), atmospheric profile, band, and date information from their metadata.

[0073] 4.2.3. In the atmospheric correction coefficient lookup table, use linear interpolation to interpolate the obtained parameter combination to obtain the correction coefficient (a, b, c). For the correction coefficient a, use the following formula to calculate:

[0074] .

[0075] in , , , , , , , Respectively represent the correction coefficient values at the eight vertices in the three-dimensional space. x, y, z are normalized interpolation parameters, and the calculation formula is:

[0076] , , .

[0077] Among them, solarz, aod, and dem are the solar zenith angle, aerosol optical depth, and elevation data of the image slice to be processed, respectively. 、 、 and 、 、 are the values of the two data points adjacent to their parameters.

[0078] Using this interpolation method, we can calculate the precise correction coefficient a for any parameter combination based on the discrete data points in the lookup table. Similarly, correction coefficients b and c can be calculated using the same interpolation formula.

[0079] By using the above interpolation method, two sets of correction coefficients are calculated for the dates $(2020 / 07 / 01,2020 / 01 / 01)$, which are recorded as , Date effect correction. When the 6S model simulates atmospheric radiation transfer, the input date not only affects the atmospheric profile but also has other effects on the simulation process. This effect is recorded here as the date effect. Further correction is required to address this date effect. The correction formula is as follows:

[0080] ;

[0081] ;

[0082] .

[0083] in, Is an experience value, specifically calculated as .

[0084] , .

[0085] After the date effect correction, the atmospheric correction parameters of the corresponding slice at its acquisition time can be obtained ( , , ).

[0086] 4.2.4. Use the atmospheric correction coefficient to correct the image slice. The formula is as follows:

[0087] ;

[0088] .

[0089] Among them, R is the radiation value obtained after radiation calibration, y is the surface reflectivity after atmospheric correction, ( , , ) is the atmospheric correction coefficient obtained in the previous step.

[0090] Step 5. Image mosaicking.

[0091] like Figure 3 As shown in the figure, different remote sensing image slices in the same geographical area are integrated into a single, continuous image through data screening, spatial stitching and other methods, which specifically includes the following sub-steps:

[0092] 5.1. Preprocess each image slice of multiple bands to remove negative values and ensure the validity and consistency of the data.

[0093] 5.2. Sort the slice data on each grid, and the sorting rules take into account the data coverage and cloud content of the slice. Calculate the coverage of each slice on the grid, that is, the proportion of non-zero pixels. Perform cloud detection on each slice, specifically through spectral feature analysis and threshold segmentation technology. First, calculate the spectral characteristics of the image segment. Then, use the Otsu method to automatically perform threshold segmentation on the spectral characteristics to determine the boundary between cloud pixels and non-cloud pixels. Next, generate a cloud mask, in which cloud pixels are marked as True and non-cloud pixels are marked as False. Based on the cloud mask, calculate the cloud coverage ratio. Sort all slices according to the two dimensions of data coverage and cloud content.

[0094] The spectral feature extraction calculation steps are as follows: Calculate the spectral sum of the visible and near-infrared bands (e.g., bands 0, 1, and 2). Set a percentage range, calculate the spectral sum threshold, and store the pixels that meet the criteria. For the selected pixels, calculate the mean of each band and use it as the spectral signature representing the cloud pixel.

[0095] 5.3. Perform shadow suppression based minimum value slice synthesis. For a grid, fill the data pixel by pixel, calculate the minimum value min_val and the second minimum value second_min_val of each pixel, and if the difference between the two is greater than 0.2, it is determined to be a shadow area. Then, a weighted synthesis method is used for the shadow area:

[0096] .

[0097] For non-shaded areas, directly fill in the minimum value of that position in all slices.

[0098] Step 6. Image color balancing.

[0099] In order to solve the color difference between slices, the global data is color-homed using the MODIS global reference basemap, which includes the following sub-steps:

[0100] 6.1. Construction of the MODIS Global Reference Basemap: Using the lower-resolution but more uniform color brightness MODIS surface reflectance data (MOD09A1), we mosaicked and normalized the bands using GEE software to output a globally unified uniform color basemap.

[0101] 6.2. Tasseled cap transform or equal-weight transform: Divide the reference image and the data to be uniformly colored into small overlapping blocks, obtain the tasseled cap transform or equal-weight transform coefficients, and calculate the brightness characteristics.

[0102] 6.3. Local block segmentation: Resolves the overall uniform color brightness deviation caused by cloud cover or cloud shadows.

[0103] 6.4. One-dimensional brightness histogram matching: Adjust the spectrum based on brightness to ensure the spectrum is consistent as much as possible.

[0104] Step 7. Mineral index calculation.

[0105] Using the processed global map, the mineral index is calculated, which includes the following sub-steps:

[0106] 7.1. Calculate the value of each mineral index using the band ratio method. For example, according to the calculation method b2 / b1 for the Ferric Iron index in official USGS documentation, divide the second band's value by the first band's value to obtain the Ferric Iron index.

[0107] 7.2. Remove masks such as water bodies, vegetation, ice and snow. Using the downloaded mask data of water bodies, vegetation, etc., remove the data in these areas in the calculated mineral index image, that is, set the values to nodata.

[0108] 7.3. Index normalization. For each index calculated, calculate the maximum and minimum values of the numerical part, then normalize the mineral index to [0, 255], save it in .png format and export it. The final result is as follows Figure 4 shown.

[0109] The embodiment of the present invention further provides a mineral index calculation system based on global multispectral remote sensing imagery, which is used to implement the above method, including:

[0110] Data acquisition module, used for downloading remote sensing data and auxiliary data sets and extracting metadata.

[0111] Data screening module is used to automatically screen remote sensing data.

[0112] The resampling and slicing module is used to resample remote sensing images and convert slicing formats.

[0113] The remote sensing image preprocessing module is used to preprocess remote sensing data, including radiation calibration and atmospheric correction.

[0114] Remote sensing mosaic module, used to stitch multiple slices.

[0115] The color uniformity module is used to uniformize remote sensing image slices.

[0116] Calculation module, used for mineral index calculation.

[0117] Result storage module, used for data conversion and storage.

[0118] Those skilled in the art will readily appreciate other embodiments of the present application after considering the specification and practicing the contents disclosed herein. This application is intended to cover any variations, uses, or adaptations of the present application that follow the general principles of this application and include common knowledge or customary techniques in the art not disclosed herein. The description and examples are to be considered as exemplary only.

[0119] It will be understood that the present application is not limited to the exact construction that has been described above and shown in the drawings, and that various modifications and changes may be made without departing from the scope thereof.

Claims

1. A mineral index calculation method based on global multispectral remote sensing imagery, characterized in that: The steps include: Acquire original remote sensing image data and parse metadata; the original remote sensing image data includes visible light data and shortwave infrared data; the metadata includes extraction time, cloud content, geographic coordinates, and quality parameters; Based on the metadata, the original remote sensing image data is screened to remove invalid data; Resample and slice the filtered remote sensing image data, store the image slices in tiff format, and write metadata into tiff files to prevent data information loss; Perform data preprocessing to convert the resampled and sliced data into data that can directly calculate the mineral index; Performing image mosaicking to integrate different remote sensing image slices within the same geographic area into a single, continuous image; performing image mosaicking includes: Preprocess each image slice of multiple bands and remove negative values to ensure the validity and consistency of the data; For each raster slice data, the data is sorted according to its coverage and cloud content. Specifically, cloud detection is performed on each slice using spectral feature analysis and threshold segmentation techniques; the spectral features of the image segment are calculated; the spectral features are automatically thresholded using the Otsu method to determine the boundary between cloud pixels and non-cloud pixels; a cloud mask is generated, in which cloud pixels are marked as True and non-cloud pixels are marked as False; based on the cloud mask, the cloud cover ratio is calculated; and all slices are sorted according to the two dimensions of data coverage and cloud content. Perform minimum value slice synthesis based on shadow suppression; for a certain grid, fill data pixel by pixel, calculate the minimum value and the second minimum value of each pixel, if the difference between the two is greater than 0.2, it is determined to be a shadow area, and weighted synthesis is used for the shadow area; for non-shadow areas, directly fill the minimum value of that position in all slices; Perform image color grading to reduce color differences between slices and obtain a global image. The image color grading includes: Using MODIS surface reflectance data, after mosaicking and band normalization, a globally unified uniform color basemap is output; Divide the reference image and the data to be uniformly colored into overlapping small blocks, obtain the tasseled cap transformation or equal-weight transformation coefficients, and calculate the brightness characteristics; Perform local block segmentation to resolve overall uniform brightness deviation caused by cloud cover or cloud shadows; Adjust the spectrum based on brightness to ensure spectral consistency; Based on the global map, mineral index indicators are calculated.

2. The method for calculating mineral index based on global multispectral remote sensing images according to claim 1, characterized in that: The resampling and slicing of the filtered original remote sensing image data includes: The original remote sensing image data is resampled, coordinate converted and sliced, specifically: setting the target resolution, converting the coordinate systems of visible light data and shortwave infrared data into a unified coordinate system; setting the target resolution, resampling the shortwave infrared band to be consistent with the visible light band, so that it is spatially aligned with the visible light band; cutting each scene image into 1024*1024 image slices, storing them in tiff format, and writing metadata files into tiff files to prevent the spectral information and spatial information of the data from being lost or distorted; wherein, the resampling uses any interpolation algorithm.

3. The mineral index calculation method based on global multispectral remote sensing imagery according to claim 1, characterized in that: The original remote sensing image data is ASTER multispectral remote sensing data.

4. The method for calculating mineral index based on global multispectral remote sensing images according to claim 3, characterized in that: The data preprocessing includes radiometric calibration and atmospheric correction, specifically: Input the quantized value of the original data and convert it into the radiation brightness value through radiation correction; according to the calibration coefficient, perform absolute calibration calculation on the data of each band to eliminate the influence of the sensor's own characteristics on the data, so that the data accurately reflects the radiation characteristics of the ground object; The 6S model lookup table is used to perform atmospheric correction on remote sensing images to eliminate the influence of atmosphere on surface reflectivity.

5. The method for calculating mineral index based on global multispectral remote sensing images according to claim 4, characterized in that: The atmospheric correction of the remote sensing image using the 6S model lookup table includes: Construct a lookup table of atmospheric correction coefficients for all parameters. Discretize and sample the parameters input by the 6S model, including solar zenith angle, solar azimuth, aerosol optical depth, digital elevation model, atmospheric profile, band, and date information. After determining the sampling points, call the 6S model to calculate the correction coefficients in the lookup table and store them in the lookup table, indexed by parameter combination. For image slices that need to be sliced, extract the solar zenith angle, solar azimuth, aerosol optical depth, digital elevation model, atmospheric profile, band and date information from their metadata; Based on the discrete data points in the atmospheric correction coefficient lookup table, the parameter combination of the 6S model input is interpolated using a linear interpolation method to obtain the correction coefficient; at the same time, the linear interpolation method is used to correct the error caused by the input date; The image slices are corrected using the atmospheric correction coefficients.

6. The method for calculating mineral index based on global multispectral remote sensing images according to claim 1, characterized in that: The calculation of mineral index indicators based on the global map includes: The value of each mineral index is calculated using the band ratio method; Obtaining ground object mask data and removing the mask area data when calculating the mineral index image; the ground object mask data includes vegetation, water, ice and snow; For each calculated index, the maximum and minimum values of the numerical part are calculated, the mineral index is normalized to [0, 255], and stored in PNG format for export.

7. A mineral index calculation system based on global multispectral remote sensing images, used to implement the method according to any one of claims 1 to 6, characterized in that: include: Data acquisition module, used for downloading remote sensing data and auxiliary data sets and extracting metadata; Data screening module, used to automatically screen remote sensing data; Resampling and slicing module, used for resampling remote sensing images and converting slice formats; Remote sensing image preprocessing module, used for data preprocessing of remote sensing data, including radiometric calibration and atmospheric correction processing; Remote sensing mosaic module, used to stitch multiple slices; Color uniformity module, used to uniformize remote sensing image slices; Calculation module, used for mineral index calculation; Result storage module, used for data conversion and storage.

Images