Mineral index calculation method and system based on global multispectral remote sensing image

Abstract

The invention discloses a mineral index calculation method and system based on a global multispectral remote sensing image, and belongs to the technical field of image data processing. The method comprises the following steps: acquiring original remote sensing image data and analyzing metadata; based on the metadata, screening the original remote sensing image data, and removing invalid data; re-sampling and slicing the remote sensing image data; data preprocessing is carried out, and data obtained after resampling and slicing are converted into data capable of directly calculating mineral indexes; performing image mosaic and image color homogenization, integrating different remote sensing image slices in the same geographic area into a single and continuous image, and reducing chromatic aberration among the slices to obtain a global image; a mineral index is calculated based on a global map. According to the method, the whole process from automatic data screening, data preprocessing, remote sensing embedding and color uniformizing to final mineral index calculation is achieved, and the pressure of storage, internal storage and the like faced by traditional software is avoided; and a lookup table mode is also used, so that the calculation time is greatly saved.

Description

Technical Field

[0001] The present invention relates to the technical field of image data processing, and in particular to a method and system for calculating mineral indices based on global multi-spectral remote sensing images Background Art

[0002] Altered rocks are rocks in which mineral composition, chemical composition, structure, texture, etc. have changed under the influence of hydrothermal processes. Since they are often found around hydrothermal deposits, they are called altered wall rocks, and altered wall rocks are an important ore prospecting indicator. Based on the spectral characteristics differences between mineralized altered rocks and wall rocks, image enhancement processing methods can be used to obtain image variables with enhanced mineralized alteration information, thereby ultimately achieving the purpose of extracting mineralized alteration information. Multi-spectral remote sensing data, with its rich spectral information and global coverage ability, has important application value in geological exploration and resource management. However, the processing and calculation of multi-spectral and multi-temporal remote sensing data still face many technical challenges globally. Remote sensing data is usually collected in a regional and time-segmented manner. Due to the imaging conditions and time differences in different regions, these data are inconsistent in radiation characteristics and geometric positions. At the same time, to achieve seamless stitching on a global scale, problems such as cloud removal, mosaicking matching, and color homogenization of cross-regional and multi-temporal data need to be solved, which poses extremely high requirements for algorithm design and data processing capabilities. In addition, geoscience researchers usually face the problem of limited computing resources and it is difficult to efficiently process large-scale data and complete complex preprocessing tasks

[0003] Traditional calculations of mineral alteration indices usually target small local areas, such as a certain metallogenic belt or ore point, screening and extracting corresponding remote sensing images, and then using remote sensing processing software for multiple steps of processing, which takes a relatively long time. This method is not suitable for the processing and calculation of large-scale remote sensing images, mainly for the following reasons. First, it is impossible to perform artificial data screening for large-scale remote sensing images. Second, ordinary remote sensing processing software cannot calculate large-scale remote sensing data and has many limitations in terms of time, efficiency, storage space, etc. Finally, ordinary software cannot perform fusion calculations on remote sensing images collected at different times and regions Summary of the Invention

[0004] Aiming at the deficiencies 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 method and system for calculating mineral indices based on global multi-spectral remote sensing images

[0005] The purpose of the present invention is achieved through the following technical solutions: A method for calculating mineral indices based on global multi-spectral remote sensing images, 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 short-wave 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%, the 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, the resampling and slicing of the filtered original remote sensing image data includes:

[0014] The original remote sensing image data is resampled, coordinate converted and sliced, specifically: set the target resolution, convert the coordinate systems of visible light data and short-wave infrared data into a unified coordinate system; set the target resolution, resample the short-wave infrared band to be consistent with the visible light band, so that it is spatially aligned with the visible light band; cut each scene image into image slices of 1024*1024 size, store them in tiff format, and write metadata files into tiff files 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 quantized value of the original data (DN), and convert it into the radiation brightness value (Radiance) through radiation correction; perform absolute calibration calculation on the data of each band according to the calibration coefficient to eliminate the influence of the sensor's own characteristics on the data, so that the data can accurately reflect the radiation characteristics of the ground object;

[0018] The atmospheric correction of remote sensing images is carried out by using the 6S model look-up table method to eliminate the influence of the atmosphere on the surface reflectivity.

[0019] Furthermore, the atmospheric correction of remote sensing images by using the 6S model look-up table method includes:

[0020] Construct an atmospheric correction coefficient look-up table for all parameters, discretize and sample the parameters input to the 6S model. The parameters include 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, call the 6S (Second Simulation of a Satellite Signal in the Solar Spectrum) model to calculate the correction coefficients in the look-up table, and store the correction coefficients in the look-up table (LOOK-UP TABLE) and index them according to the parameter combination.

[0021] For the image slices that need to be sliced, extract the solar zenith angle, solar azimuth angle, 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 look-up table, use the linear interpolation method to interpolate the parameter combinations input to the 6S model to obtain the correction coefficients. At the same time, use the linear interpolation method to correct the error caused by the input date.

[0023] Correct the image slices by the atmospheric correction coefficients.

[0024] Furthermore, the image mosaicking includes:

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

[0026] For the slice data on each grid, sort it according to the coverage and cloud content of the data on the slice.

[0027] Perform minimum slice synthesis based on shadow suppression. For a certain grid, fill the 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 as the shadow area, and the weighted synthesis is used for the shadow area. For the non-shadow area, directly fill the minimum value at this position in all slices.

[0028] Furthermore, for the slice data on each grid, the sorting is specifically 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 technology; the spectral features of the image fragments 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 coverage ratio is calculated; all slices are sorted according to the two dimensions of data coverage and cloud content.

[0030] Furthermore, the image color uniformity comprises:

[0031] Using MODIS surface reflectance data, after mosaicking and band normalization, a globally unified uniform color base map 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 the 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] Using the band ratio method, the value of each mineral index is calculated;

[0037] Acquire ground object mask data, and remove 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 slice module, used for resampling remote sensing images and converting slice formats;

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

[0044] A remote sensing mosaicking module for mosaicking multiple slices;

[0045] A color homogenization module for homogenizing remote sensing image slices;

[0046] A calculation module for calculating mineral indices;

[0047] A result storage module for data conversion and storage.

[0048] The beneficial effects of the present invention are as follows: The present invention proposes a calculation framework applicable to global remote sensing image processing, which transforms the traditional processing of remote sensing data in units of scenes into processing remote sensing data in units of slices, realizing the whole process from automatic data screening, remote sensing mosaicking and color homogenization to the final calculation of mineral indices. Moreover, this calculation framework divides the processing tasks into calculation sub-units, which can be placed in a parallel computing environment, avoiding the storage, memory and other pressures faced by traditional software processing; in the time-consuming step of atmospheric correction, the online calling algorithm is changed to a look-up table method, which can avoid repeatedly calling the atmospheric correction algorithm on large-scale slice data, greatly saving calculation time. Description of the Drawings

[0049] Figure 1 It is a flowchart of the method of the present invention;

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

[0051] Figure 3 It is a flowchart of the remote sensing mosaicking of the present invention;

[0052] Figure 4 It is a schematic diagram of the remote sensing image processing result of the present invention. Detailed Embodiments

[0053] Here, the exemplary embodiments will be described in detail, and the examples are shown in the drawings. When the following description refers to the drawings, unless otherwise indicated, the same numbers in different drawings represent the same or similar elements. The embodiments described in the following exemplary embodiments do not represent all embodiments consistent with the present invention. On the contrary, they are merely examples of devices and methods consistent with some 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 "said" used in this invention and the appended claims are also intended to include the plural forms unless the context clearly dictates otherwise. It should also be understood that the term "and / or" used herein refers to and encompasses any and all possible combinations of one or more of the associated listed items.

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

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

[0057] 1.1. Download ASTER multispectral remote sensing data as the original processed remote sensing image data. Select the first to the ninth bands of ASTER, including visible light and short-wave infrared data. Download the corresponding ASTER GDEM elevation data and Modis AOD data at the corresponding time as auxiliary data required for atmospheric correction processing. Download vegetation index and water body data as masking data.

[0058] 1.2. Parse the metadata in the xml of the ASTER remote sensing data to extract information such as time, cloud cover, geographical coordinates, and quality parameters.

[0059] Step 2. Data screening.

[0060] Conduct a preliminary screening of the downloaded original remote sensing data. Set screening rules according to cloud cover and image anomalies. When the quality score is 1 or the cloud cover is 40%, the data is excluded. Determine the availability of the data according to the extracted metadata information and mark the qualified data.

[0061] Step 3. Data resampling and slicing.

[0062] Due to the overlapping and intersecting areas in remotely sensed images collected at different times, resampling, coordinate transformation, and slicing are performed on the original remotely sensed images. The original visible light band has a resolution of 30 meters, and the short-wave infrared band has a resolution of 100 meters. The target resolution is set to 30 meters, and the coordinate systems of the visible light and short-wave infrared image data are converted to the unified EPSG: 3857 coordinate system. Then, with the target resolution set to 30 meters, the original visible light band (30-meter resolution) and the short-wave infrared band (100-meter resolution) are resampled to 30-meter resolution. Interpolation algorithms such as bilinear interpolation and cubic convolution interpolation can be used for resampling. Finally, each scene of the image is sliced into image slices of size 1024*1024 and stored in the tiff format. During the format conversion process, to ensure that key contents such as spectral information and spatial information of the data are not lost or distorted, metadata files will be written into the tiff files for subsequent data management and query.

[0063] Step 4. Data preprocessing.

[0064] Data preprocessing converts the original data quantization value (DN) collected by the remote sensing sensor into data that can directly calculate the mineral index after removing the influences of factors such as the sensor and atmospheric aerosol.

[0065] 4.1. Radiometric calibration. The original remote sensing digital quantization value DN is input and converted into radiance value (Radiance) through radiometric correction. According to the calibration coefficients provided by the USGS, precise absolute calibration calculations are performed on the data of each band to eliminate the influence of the sensor's own characteristics on the data, enabling the data to accurately reflect the radiation characteristics of the ground objects. The radiometric calibration calculation is shown in the following formula:

[0066] 。

[0067] Among them, m is the gain coefficient. The gain mode of this band is obtained from the information recorded in the image header data, and then the specific gain coefficient is obtained from the gain matrix in the USGS official 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 reflectance, realizes the precise correction of remote sensing data, and obtains more accurate surface reflectance information. In the present invention, the 6S model lookup table method is used to perform atmospheric correction on remotely sensed images. The traditional 6S model needs to query the aerosol optical thickness, elevation data corresponding to the remotely sensed image, and the metadata information of the remote sensing itself, input them into the model to obtain the correction coefficient, and perform transformation on the image. This process needs to be carried out for each slice, which is very time-consuming for the calculation of large-scale data. Therefore, the present invention uses the lookup table method to perform 6S model atmospheric correction.

[0069] 4.2.1. Construct a lookup table of atmospheric correction coefficients for all parameters, as shown in Figure 2 . For the solar zenith angle (SZA), solar azimuth angle (SAA), aerosol optical depth (AOD), digital elevation model (DEM), atmospheric profile, band, and date information input into the 6S model, discretize and sample the parameters as shown in Table 1 below. The distribution design of the sampling points takes into account the variation range of the parameters and the calculation efficiency. The sampling is dense in the low-value area and sparse in the high-value area to balance the accuracy of the lookup table and the 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 angle (SAA) Fixed at 90 degrees

[0071] After determining the sampling points, call the 6S (Second Simulation of a Satellite Signal in the Solar Spectrum) model to calculate the correction coefficients in the lookup table. The present invention uses the Py6S library encapsulated in Python to call the 6S model for calculation. Store the obtained correction coefficients in the lookup table (LOOK-UP TABLE) and index them according to the 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 coefficients (a, b, c). For the correction coefficient a, calculate it through the following formula:

[0074] .

[0075] Where , , , , , , , They respectively represent the calibration coefficient values at 8 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 respectively the solar zenith angle, aerosol optical depth, and elevation data of the image slice to be processed. , , and , , are respectively the values of two data points adjacent to their parameters.

[0078] Through the above interpolation method, the exact calibration coefficient a under any parameter combination can be calculated based on the discrete data points in the look-up table. Similarly, the calibration coefficients b and c can also be obtained by the same interpolation formula.

[0079] Through the above interpolation method, two sets of calibration coefficients on the dates of (2020 / 07 / 01, 2020 / 01 / 01) are calculated respectively, denoted as , Date effect correction. When the 6S model simulates the atmospheric radiative transfer, in addition to affecting the atmospheric profile, the input date also has other effects on the simulation process, and this effect is recorded as the date effect here. For the above date effect, further correction is required. The correction formula is as follows:

[0080] ;

[0081] ;

[0082] .

[0083] Among them, is an empirical value, and the specific calculation is .

[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 radiance value obtained after radiometric calibration, y is the surface reflectance after atmospheric correction, and ( , , ) is the atmospheric correction coefficient obtained from the previous step.

[0090] Step 5. Image Mosaic.

[0091] As Figure 3 shown, different remote sensing image slices within the same geographical area are integrated into a single, continuous image through methods such as data screening and spatial mosaicking, which specifically include 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. The sorting rule considers the coverage and cloud content of the data in the slice. Calculate the coverage of each slice on this grid, that is, the proportion of non-zero pixels. Perform cloud detection on each slice, specifically through spectral feature analysis and threshold segmentation techniques. First, calculate the spectral features of the image segment. Then, use the Otsu method to automatically threshold the spectral features to determine the boundary between cloud pixels and non-cloud pixels. Next, generate a cloud mask, where cloud pixels are marked as True and non-cloud pixels are marked as False. According to the cloud mask, calculate the cloud coverage ratio. Sort all slices according to the two dimensions of data coverage and cloud content.

[0094] Among them, the calculation steps for spectral feature extraction are as follows: Calculate the sum of the spectra of the visible and near-infrared bands (such as bands 0, 1, 2). Set the percentage range, calculate the spectral sum threshold, and store the pixels that meet the conditions. Calculate the mean value of each band for the selected pixels as the spectral features representing cloud pixels.

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

[0096] 。

[0097] For non-shadow areas, directly fill the minimum value at this position in all slices.

[0098] Step 6. Image color homogenization.

[0099] To solve the color difference between slices, the global data is color homogenized using the MODIS global reference base map, which specifically includes the following sub-steps:

[0100] 6.1. Construction of the MODIS global reference base map: Using the MODIS surface reflectance data (MOD09A1) with lower resolution but more consistent color brightness, mosaicking and band normalization are performed through the GEE software to output a globally unified color homogenized base map.

[0101] 6.2. Tasseled cap transformation or equal-weight transformation: The reference map and the data to be color homogenized are divided into overlapping small tiles, the tasseled cap transformation or equal-weight transformation coefficients are obtained, and the brightness characteristics are calculated.

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

[0103] 6.4. One-dimensional brightness histogram matching: Based on brightness, the spectrum is adjusted to ensure as much spectral consistency as possible.

[0104] Step 7. Mineral index calculation.

[0105] Using the processed global single map, mineral index indicators are calculated, which specifically includes the following sub-steps:

[0106] 7.1. Using the band ratio method, calculate the value of each mineral index. For example, according to the calculation method of the Ferric Iron index in the official USGS document, b2 / b1, divide the data of the second band by the value of the first band 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., in the calculated mineral index image, remove the data in these areas, that is, set the values to nodata.

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

[0109] The embodiment of the present invention also provides a mineral index calculation system based on global multi-spectral remote sensing images for implementing the above method, including:

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

[0111] A data screening module for automatically screening remote sensing data.

[0112] A resampling slicing module for resampling remote sensing images and converting the slice format.

[0113] A remote sensing image preprocessing module for preprocessing remote sensing data, including radiometric calibration and atmospheric correction processing.

[0114] A remote sensing mosaicking module for stitching multiple slices.

[0115] A color equalization module for equalizing remote sensing image slices.

[0116] A calculation module for calculating mineral indices.

[0117] A result storage module for data conversion and storage.

[0118] Those skilled in the art will readily conceive of other embodiments of the present application after considering the specification and practicing the content disclosed herein. The present application is intended to cover any variations, uses, or adaptations of the present application, which follow the general principles of the present application and include common general knowledge or conventional technical means in the technical field not disclosed in the present application. The specification and examples are only illustrative.

[0119] It should be understood that the present application is not limited to the precise structures described above and shown in the drawings, and various modifications and changes can be made without departing from its scope.

Claims

1. A method for calculating mineral index based on global multispectral remote sensing images, 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 short-wave 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; Perform image mosaicking to integrate different remote sensing image slices within the same geographic area into a single, continuous image; Perform image color grading to reduce color differences between slices and obtain a global image; 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: set the target resolution, convert the coordinate system of visible light data and short-wave infrared data into a unified coordinate system; set the target resolution, resample the short-wave infrared band to be consistent with the visible light band, so that it is spatially aligned with the visible light band; cut each scene image into image slices of 1024*1024 size, store them in tiff format, and write metadata files into tiff files to prevent the spectral information and spatial information of the data from being lost or distorted; wherein, the resampling selects any interpolation algorithm.

3. The method for calculating mineral index based on global multispectral remote sensing images according to claim 1, characterized in that: The original remote sensing image data uses ASTER multispectral remote sensing data.

4. The method for calculating mineral index based on global multispectral remote sensing images according to claim 3 is characterized in that: The data preprocessing includes radiation calibration and atmospheric correction, specifically: Input the quantized value of the original data and convert it into the radiation brightness value through radiation correction; perform absolute calibration calculation on the data of each band according to the calibration coefficient to eliminate the influence of the sensor's own characteristics on the data, so that the data can accurately reflect 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 by using the 6S model lookup table includes: Construct an atmospheric correction coefficient lookup table 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 point, call the 6s model to calculate the correction coefficient in the lookup table, and store the correction coefficient in the lookup table according to the parameter combination index; 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, a linear interpolation method is used to interpolate the parameter combination input by the 6s model to obtain the correction coefficient; at the same time, a 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 image mosaicking comprises: Preprocess each image slice of multiple bands to remove negative values ​​to ensure the validity and consistency of the data; For each slice data on the grid, the data is sorted according to its coverage and cloud content on the slice; Perform minimum slice synthesis based on shadow suppression. For a certain grid, fill data pixel by pixel, calculate the minimum and second minimum values ​​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.

7. The method for calculating mineral index based on global multispectral remote sensing images according to claim 6, characterized in that: The specific method for sorting the slice data on each grid according to the coverage and cloud content of the data on the slice is as follows: Cloud detection is performed on each slice through spectral feature analysis and threshold segmentation technology; the spectral features of the image fragments 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 coverage ratio is calculated; all slices are sorted according to the two dimensions of data coverage and cloud content.

8. The method for calculating mineral index based on global multispectral remote sensing images according to claim 1, characterized in that: The image color uniformity comprises: Using MODIS surface reflectance data, after mosaicking and band normalization, a globally unified uniform color base map 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 the overall uniform brightness deviation caused by cloud cover or cloud shadows; Adjust the spectrum based on brightness to ensure spectral consistency.

9. The method for calculating mineral index based on global multispectral remote sensing images according to claim 1, characterized in that: The mineral index indicators calculated based on the global map include: Using the band ratio method, the value of each mineral index is calculated; Acquire ground object mask data, and remove 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.

10. A mineral index calculation system based on global multispectral remote sensing images, 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 slice module, used for resampling remote sensing images and converting slice formats; Remote sensing image preprocessing module, used to preprocess remote sensing data, including radiation calibration and atmospheric correction processing; Remote sensing mosaic module, used to stitch multiple slices; Color homogenization module, used to homogenize remote sensing image slices; Calculation module, used for mineral index calculation; Result storage module, used for data conversion and storage.

Images