A rapid reconstruction method of historical land cover based on GlobeLand30

Abstract

The present invention relates to a method for rapid reconstruction of historical surface cover based on GlobeLand30, comprising the following steps: 1) Extraction of effective patches from GlobeLand30: Using morphological operations to remove erroneous patches, and obtaining patches that can correctly express the spatial distribution of surface cover information. 2) Pseudo-sample selection: Calculating the spectral characteristics of the image, using the optimized patches as clustering units, and using an unsupervised algorithm to select a pseudo-sample set. 3) Global sample optimization: Constructing Gaussian mixture models for different types of land objects, and optimizing global samples through the process of solving the model. 4) Surface cover reconstruction: Obtaining a training sample set based on the results of the third step, determining whether change detection is required, inputting the sample set into a random forest, and reconstructing a surface cover classification map of the target historical phase. This algorithm takes into account the validity and uncertainty of the information in the surface cover product, and rapidly reconstructs the historical surface cover by effectively reducing uncertainty.

Description

Technical Field

[0001] The invention relates to a method for quickly reconstructing historical land cover based on GlobeLand30, and belongs to the technical field of remote sensing intelligent information extraction. Background Art

[0002] Land cover is a key indicator for describing the composition of the earth's surface and has a significant impact on the global land surface water, energy and material cycles. Supervised classification and interpretation of land cover by remote sensing images has become a mainstream technical method. Through visual interpretation of remote sensing images, training samples are manually selected, and multi-period land cover classification results are generated based on supervised classification algorithms. This type of method has been widely used in remote sensing land cover classification due to the effectiveness of supervised classification methods. The supervised classification results are highly dependent on training samples, and the acquisition of training samples is particularly important for the above methods. Sample selection is a process that relies on expert knowledge and takes a certain amount of time and manpower. This process requires high-resolution images as a reference or field surveys to obtain real surface data.

[0003] Therefore, for the reconstruction of historical image land cover, traditional remote sensing image classification methods based on machine learning have the following limitations: 1) The ideal supervised classification results rely on sufficient high-quality training samples, which requires a lot of time and manpower; 2) For historical images, high-resolution images used for marking references are very scarce, and it is difficult to cover the entire study area, so that the training sample set cannot well represent the distribution of objects in the entire study area. To address the above problems, a fast and unsupervised historical land cover reconstruction method is proposed, which makes full use of existing land cover products and realizes the automatic acquisition of historical image training samples. An unsupervised sample migration method is constructed and embedded in the land cover classification framework to quickly and accurately reconstruct historical land cover. Summary of the invention

[0004] The technical problem to be solved by the present invention is: to overcome the difficulties and high costs of obtaining training samples in historical land cover classification tasks, and to propose an unsupervised sample migration method. The method utilizes archived land cover product data, and fully reduces the spatiotemporal uncertainty of land cover products through the geometric constraints of product patches, the attribute constraints of image spectral characteristics, and the distribution constraints of land objects, automatically and quickly produces high-quality historical image training samples, and combines change detection and classification algorithms to quickly and accurately complete the reconstruction of historical land cover.

[0005] In order to solve the above technical problems, the present invention provides a method for rapid reconstruction of historical land cover based on GlobeLand30, comprising the following steps:

[0006] Step 1: Prepare the surface cover product GlobeLand30, remote sensing image Landsat data and DEM data, select the area of ​​interest and generate the vector file of the area of ​​interest, and crop GlobeLand30, DEM and remote sensing image Landsat data;

[0007] The second step is to optimize the land patches of the land cover product based on morphological opening operation to obtain patch units: the land cover product GlobeLand30 is divided into several binary layers according to categories, and the morphological opening operation with increasing window size is performed on each layer of binary image. According to the calculation results under different window sizes, when the current background pixel ratio is greater than the user-set threshold, the binary layers are merged to obtain the optimized patch units;

[0008] Step 3: Calculate the image normalized difference spectral vector NDSV: traverse all the pairwise combinations of the reflectance of the bands, calculate the normalized difference index for each pair of band combinations, and obtain the corresponding normalized difference spectral vector NDSV. All normalized difference spectral vectors NDSV constitute the image spectral feature set.

[0009] The fourth step is to use each patch unit as a clustering calculation unit, use the K-means algorithm for unsupervised clustering, and use the K-means++ algorithm to initialize the cluster center of each patch unit one by one: for any patch unit, randomly select a normalized difference spectral vector NDSV of a pixel as the cluster center of the first cluster. When there are k cluster centers, calculate the normalized difference spectral vector NDSV of other pixels in the current patch unit to the kth cluster center μ k The Euclidean distance is calculated, and the sampling weight of the corresponding pixel is set according to the Euclidean distance. The k+1th cluster center is randomly selected according to the sampling weight, and the process is repeated until the number of clusters meets the user's preset requirements;

[0010] Step 5: Use the variance ratio criterion VRC to optimize the number of clusters for each patch unit: define a set of cluster numbers, for any patch unit, input the normalized difference spectral vector NDSV of all pixels into the K-means algorithm for clustering, calculate the variance ratio criterion VRC result based on the clustering result, and obtain the cluster number K corresponding to the maximum variance ratio criterion VRC by traversing all elements in the cluster number set. j As the optimal number of clusters for clustering within the current patch unit;

[0011] Step 6. According to the methods in the fourth and fifth steps, the K-means method is performed on all patch units, and the global image pseudo sample set is obtained according to the clustering result of the optimal number of clusters: the spectral image feature set obtained in the third step is used as the K-means input feature, and all patch units obtained in the second step are traversed. The cluster center is obtained for all pixels in each patch unit using the method in the fourth step, and the optimal number of clusters is obtained using the method in the fifth step. On this basis, the clustering result is obtained, and the cluster with the largest number of pixels in the patch unit inherits the category information of the patch unit to obtain the pseudo sample set;

[0012] Step 7. Determine the land cover classification system based on the land cover product, select the corresponding key features for each land category, and construct multiple two-class Gaussian mixture models: calculate the improved normalized water index MNDWI of the global image as the key feature of the water body type, calculate the normalized vegetation index NDVI of the non-water area as the key features of the forest, grassland, and artificial surface types, and calculate the slope of the vegetation area as the key feature of the cultivated land type. The characteristic value distribution of each key feature is regarded as a mixture of two Gaussian distributions, and the following Gaussian mixture model is constructed:

[0013]

[0014] in is the probability density function, y represents the eigenvalue of the key feature, α 1 , α 2 are the mixing coefficients of the first and second Gaussian distributions, μ 1 , μ 2 are the mean vectors of the two distributions, ∑ 1 ,∑ 2 The distribution is the covariance vector of the two distributions, p(y|μ 1 ,∑ 1 ) and p(y|μ 2 ,∑ 2 ) are the first and second Gaussian distribution density functions respectively. j ,μ j ,∑ j )|j=1,2} are collectively referred to as mixing component parameters;

[0015] Step 8: Use the iterative expectation maximization EM algorithm to estimate the mixed component parameters in each Gaussian mixture model, complete the decomposition of the Gaussian mixture model, and thus complete the optimization of the pseudo sample set to obtain the training sample set;

[0016] Step 9. Use change detection and classification methods to complete the surface cover reconstruction of the target historical phase: determine whether the target phase and the surface cover product belong to the same year. If they belong to the same year, input the optimized training samples in step 8 into the random forest RF classifier to complete the training and obtain the surface cover classification result of the target phase; if they belong to different years, use the dual-phase image change detection method to obtain the invariant area, transfer the sample labels of the invariant area, combine the target phase image to form training samples, and input the training samples into the random forest RF classifier to complete the training and obtain the surface cover classification result of the target historical phase.

[0017] The historical land cover reconstruction method used in the present invention is an algorithm innovation. The core idea of ​​the algorithm is to use the information in the archived land cover products, while considering the validity and uncertainty of the product information, and preliminarily reduce the uncertainty of the land cover information of the product from the perspective of the land patch information in the land cover product. The present invention introduces an improved K-means algorithm, and quickly divides and obtains a pseudo sample set under the geometric constraints of the patch unit and the attribute constraints of the spectral characteristics. The present invention introduces a Gaussian mixture model, utilizes the distribution constraints of different land types in the key feature space, and optimizes the pseudo sample set through the mixed model decomposition to obtain a reliable training sample set. Based on two basic clustering algorithms, the present invention effectively reduces the uncertainty of product information from three perspectives: the geometric constraints of the land cover product patches, the attribute constraints of the image spectral characteristics, and the distribution constraints of the land type feature space, and automatically obtains a reliable historical phase training sample set. The algorithm proposed in the present invention is a fast and automatic historical land cover reconstruction method, which provides a solution to the problems of intelligent interpretation of remote sensing images and dynamic mapping of land cover. All steps are implemented by Matlab programming, and no manual participation is required except for setting the necessary parameters of the algorithm.

[0018] The main innovations of the algorithm are in the following three aspects:

[0019] 1. In the second step, a patch unit optimization method based on morphological opening operation is constructed. The purpose of removing small patches and erroneous patches in the original product is achieved through the opening operation of binary images category by category, so as to obtain a patch unit that better meets the algorithm requirements.

[0020] 2. In the sixth step, a pseudo-sample selection method based on an improved K-means algorithm is proposed. The algorithm uses the land patches optimized in the second step as clustering units and the NDSV features calculated in the third step as clustering features, thereby significantly reducing product uncertainty from local units. At the same time, the calculation units of the pseudo-sample selection algorithm are independent of each other, and parallel computing is used to ensure the efficiency of the algorithm operation.

[0021] 3. In the seventh step, a global sample optimization method based on a Gaussian mixture model is constructed. By selecting key features related to the surface cover type, the pseudo samples are globally optimized from the perspective of feature space data distribution to obtain a reliable training sample set.

[0022] The method of the present invention has good sample extraction effect and high classification accuracy in large-scale surface cover reconstruction of multiple historical years. It is more accurate than the original surface cover product and plays a good reference role in medium and high resolution historical surface cover reconstruction and surface cover update of new images. BRIEF DESCRIPTION OF THE DRAWINGS

[0023] The present invention will be further described in conjunction with the accompanying drawings.

[0024] Figure 1 Flowchart of the rapid reconstruction method of historical land cover based on GlobeLand30.

[0025] Figure 2 Flowchart of pseudo sample selection method.

[0026] Figure 3 Flowchart for global pseudo sample optimization based on Gaussian mixture models.

[0027] Figure 4(a) is the median composite image for the year 2000.

[0028] Figure 4(b) shows the land cover classification results in 2000.

[0029] Figure 5(a) is the median composite image for 1995.

[0030] Figure 5(b) shows the land cover classification results in 1995. DETAILED DESCRIPTION

[0031] The technical route and specific operation steps of the present invention are described in detail according to the accompanying drawings.

[0032] The example study area of ​​the present invention is the Taihu Lake Basin, the land cover product selected is the GlobeLand30 series, and the product data was acquired in 2000. For the convenience of description, the product will be referred to as GlobeLand30V20000 in the following, and the image acquisition time is 1995 and 2000 respectively.

[0033] This implementation case is based on the historical land cover rapid reconstruction method of GlobeLand30 ( Figure 1 ), including the following steps:

[0034] Step 1: Download GlobeLand30 V2000 data products (spatial resolution is 30m) and clip the data products according to the vector files of the test area. Obtain the Landsat median composite reflectance data (spatial resolution is 30m) through Google Earth Engine. The directory of this data product in GEE is "LANDSAT / LT05 / C02 / T1_L2". The Landsat data in this directory are accompanied by mask files of clouds, shadows, snow and other objects generated by the CFMask algorithm. Based on the mask files, obtain all clear sky observation data in 1995 and 2000 in GEE, output images by median synthesis, clip images according to the vector files of the test area, download the Digital Elevation Model (DEM) data, and clip according to the vector files of the test area.

[0035] Step 2: Use Matlab programming to complete the land patch optimization method. Read the raster format GlobeLand30V2000, decompose GlobeLand30 V2000 into multiple binary layers according to the category raster value, set the structural element shape to square, the domain size range to [3,21], and the step size to 2. Use the imopen function in the Matlab image processing toolbox to complete the morphological opening operation. After completing the morphological opening operation of each window value, merge the binary images and calculate the background pixel ratio, which is recorded as ξ. In order to ensure that the product information is retained as much as possible while optimizing the patch, the threshold of ξ is set to 5%. When the background pixel ratio is higher than 5% (user-set threshold), the optimization algorithm is stopped, and the merged result under the current window value is used as the optimized patch unit.

[0036] Step 3: Use Matlab programming to complete the extraction of the normalized difference spectral vector NDSV, read the Landsat 5 image synthesized by the median value in 2000, traverse the pairwise combination of all bands, calculate the normalized difference index for each pair of two bands, and obtain the corresponding normalized difference spectral vector NDSV. All normalized difference spectral vectors NDSV form a more discriminative image spectral feature set for each land type. The specific calculation formula of the normalized difference spectral vector NDSV is shown below:

[0037]

[0038]

[0039] where y d Represents the spectral feature vector of the dth pixel in the image, The bth pixel is the dth pixel. i The reflectivity of each band, The bth pixel is the dth pixel. j The reflectivity of each band, the value range of i and j is [1, B].

[0040] Step 4: Take each patch unit as a clustering calculation unit, use the K-means algorithm for unsupervised clustering, and use the K-means++ algorithm to initialize the cluster center of each patch unit one by one. j According to the third step NDSV calculation method, the corresponding spectral feature set is obtained from A normalized difference spectral vector is randomly selected as the cluster center of the first cluster, denoted by μ 1 ; then calculate The Euclidean distance from each normalized difference spectral vector to the cluster center in The next cluster center μ is randomly selected according to the sampling weight 2 , by normalizing the difference spectrum vector to μ 1 The Euclidean distance is used to set the corresponding sampling weight. For the nth feature vector y n In terms of sampling weight The calculation method is as follows:

[0041]

[0042] When the kth cluster center is obtained, Continue to randomly select the next cluster center according to the weight. The eigenvector y n The weight is calculated as follows:

[0043]

[0044] Repeat the process from k to k+1 until the number of clusters meets the preset requirements.

[0045] Step 5: Use Matlab programming to implement the variance ratio criterion VRC algorithm to automatically find the optimal number of clusters for each patch unit. j For example, K-means is used to transform the spectral feature set Clustering is K j clusters, and based on the clustering results, obtain the overall inter-cluster variance SS under the current clustering results B (overall between-cluster variance), the overall between-cluster variance SS under the current clustering result B The calculation method is as follows:

[0046]

[0047] where μ k It is a cluster The mean vector, K j is the number of clusters and μ is the overall mean vector of the data set. It is a cluster The total inter-cluster variance SS B The larger the value, the lower the similarity between clusters.

[0048] Calculate the overall intra-cluster variance SS under the current clustering result W (overall within-cluster variance), calculated using the following formula:

[0049]

[0050] where y n is the eigenvector, μ k It is a cluster The mean vector, K j is the number of clusters, and the smaller the overall intra-cluster variance is, the higher the similarity between the data within the cluster is.

[0051] The variance ratio criterion measures the clustering effect by calculating the ratio of the overall inter-cluster variance to the overall intra-cluster contrast, that is, ensuring the similarity of intra-cluster data and the difference of inter-cluster data. The specific measurement index is calculated as follows:

[0052]

[0053] Where N is the number of eigenvectors, K j is the number of clusters, and the optimal number of clusters is set to range from 1 to 10.

[0054] Step 6: Use Matlab programming to implement the pseudo sample selection algorithm ( Figure 2), the proposed pseudo sample selection algorithm requires three types of input data. First, the patch units optimized in the second step are sequentially encoded to obtain the encoded raster data of the patch units, which are read using Matlab; secondly, the feature image of the normalized difference spectral vector NDSV is read, and finally the surface cover product raster data is read. According to the order of the encoding information, individual patch units are obtained one by one as the geometric constraints of the pseudo sample selection algorithm, that is, the clustering units of the improved K-means method. At the same time, the NDSV feature image subset and surface cover product category information corresponding to the patch unit are obtained, and the NDSV features are used as attribute constraints, that is, the clustering features of the improved K-means method. For each independent patch unit, the improved K-means algorithm is executed. After completing the initialization and cluster number optimization, the kmeans function in the Matlab Statistics and Machine Learning Toolbox is used to select the input data as the NDSV features in the patch unit. The number of clusters is the result after optimization. The maximum number of iterations is set to 1000 to complete the K-means clustering. The algorithm defines the initial K j Cluster centers are calculated by Each normalized difference spectrum vector y in n To each cluster center Euclidean distance, y n Divide into clusters with the most similar cluster centers. Minimize the square error by continuously adjusting the cluster centers:

[0055]

[0056] Among them, μ k It is a cluster The mean vector of .

[0057] When the K-means algorithm converges, the current cluster division results are obtained, the cluster with the highest proportion is retained, and the surface cover product category information is assigned to the cluster. Through parallel computing, the pseudo-sample selection algorithm is completed for all patch units in the test area to obtain the pseudo-sample selection results for the entire test area.

[0058] Step 7: Select key features for each land type and choose the modified normalized difference water index (MNDWI) as the discriminant information between water bodies and other land types. MNDWI is obtained by the normalized difference between the reflectance of the green band and the reflectance of the mid-infrared band:

[0059]

[0060] where ρ MIR and ρ Green Represent the reflectivity of the mid-infrared band and the green light band respectively.

[0061] The Normalized Difference Vegetation Index (NDVI) is selected as the key indicator of forest land, grassland and artificial surface. NDVI is obtained by the normalized difference index between the near infrared band and the red light band:

[0062]

[0063] where ρ NIR and ρ Red Represent the reflectivity of the near-infrared band and the red light band respectively.

[0064] Slope is selected as the key indicator of cultivated land. The slope is calculated based on DEM data and obtained through ArcGIS software. First, select the spatial selection tool, open the "Surface" sub-bar, click the "Slope" tool, and enter the DEM calculation.

[0065] Using Matlab programming to build a multi-layer two-class Gaussian mixture model ( Figure 3 ), by considering the MNDWI of the overall image as a mixture of two types of Gaussian distributions, on the basis of solving the two-type Gaussian mixture model, the clusters are marked as water bodies and non-water bodies by comparing the cluster means, and the non-water body part of the global image continues to construct a Gaussian mixture model through NDVI and slope characteristics respectively. The Gaussian mixture model constructed by NDVI is used to optimize the pseudo samples of forest and grassland and artificial surface at the same time, and the pseudo samples of cultivated land are decomposed and optimized by the Gaussian mixture model constructed by slope. The construction of the two-type Gaussian mixture model is completed by programming using the fitgmdist function in the Matlab Statistics and Machine Learning Toolbox.

[0066] The Gaussian mixture model formula is as follows:

[0067]

[0068] in is the probability density function, y represents the eigenvalue of the key feature, α 1 , α 2 are the mixing coefficients of the first and second Gaussian distributions, μ 1 , μ 2 are the mean vectors of the two distributions, ∑ 1 ,∑ 2 The distribution is the covariance vector of the two distributions, p(y|μ 1 ,∑ 1 ) and p(y|μ 2 ,∑ 2 ) are the first and second Gaussian distribution density functions respectively. j ,μ j ,Σ j)|j=1,2} are collectively referred to as mixing component parameters.

[0069] Step 8: Use Matlab programming to complete the global pseudo sample optimization based on Gaussian mixture model decomposition to obtain a reliable training sample set. Use the EM algorithm to decompose the multiple second-class Gaussian mixture models constructed in step 7. First, use the Mahalanobis distance as the weight and use the weight sampling method to initialize the EM algorithm. Then, perform the E step of the EM algorithm and use the mixed component parameters of the current partition result to estimate the elements in the data set. The posterior probability of belonging to the first and second Gaussian mixture models is calculated as:

[0070]

[0071] in is the posterior probability of membership j, {(α j ,μ j ,∑ j )|j=1,2} is the mixed component parameter. On this basis, M steps are performed to adjust the mixed distribution parameters through maximum likelihood estimation. The log-likelihood function that needs to be maximized is:

[0072]

[0073] in is the current data set, N i is the number of elements in the dataset.

[0074] The stopping condition of the EM algorithm is set to the point where the likelihood function is no longer significantly changed. In this embodiment, the specific parameter is ΔLL≤10 -6 After the decomposition of multiple Gaussian mixture models is completed, the samples with wrong labels in the pseudo sample set are removed to obtain the final training sample results.

[0075] Step 9: Use Matlab programming to implement change detection and classification algorithms to obtain land cover reconstruction results. For the land cover reconstruction in 2000, the training samples obtained in step 8 and the image data in 2000 were directly input into the Random Forest classifier. The RF classifier was implemented using open source code (Classification and regression based on a forest of trees using random inputs, based on Breiman (2001)<DOI:10.1023 / A:1010933404324> .), where the number of base classifiers of RF is set to 500 and the dimension of the feature subset of RF is set to the square root of the input feature dimension. The land cover reconstruction results in 2000 were obtained using RF (Figure 4).

[0076] For the task of land cover reconstruction in 1995, considering the impact of land cover changes, the change vector analysis method combined with the Otsu threshold segmentation method was used to detect the unchanged areas between 1995 and 2000. The change detection algorithm was implemented by Matlab programming. First, the change intensity information was calculated:

[0077]

[0078] where ρ d is the change intensity information of the th pixel, and are the reflectance values ​​of the pixel in the first and second phases in the n-band, respectively, and N is the total number of image bands. Then, the Otsu method is used to binary the intensity image to obtain the change and unchanged areas. The Otsu method calculates the inter-class variance under different gray level thresholds by traversing all gray levels in the image. The threshold T corresponding to the maximum value in the set is taken as the optimal threshold. For an image with gray levels {0,1,2,3,…,l,…L-1}, the optimal threshold T * for:

[0079]

[0080] Where T is the threshold, L is the number of gray levels of the image, The intensity image is binarized according to the optimal threshold to obtain the invariant area, and the training sample labels of the invariant area are transferred from the 2000 image to the 1995 image to form the 1995 training samples, which are input into the RF classifier to complete the 1995 land cover reconstruction (Figure 5).

[0081] In addition to the above embodiments, the present invention may also have other implementation modes. Any technical solution formed by equivalent replacement or equivalent transformation falls within the protection scope required by the present invention.

Claims

1. A method for rapid reconstruction of historical land cover based on GlobeLand30, The following steps are involved: Step 1: Prepare the surface cover product GlobeLand30, remote sensing image Landsat data and DEM data, select the area of ​​interest and generate the vector file of the area of ​​interest, and crop GlobeLand30, DEM and remote sensing image Landsat data; The second step is to optimize the land patches of the land cover product based on morphological opening operation to obtain patch units: the land cover product GlobeLand30 is divided into several binary layers according to categories, and the morphological opening operation with increasing window size is performed on each layer of binary image. According to the calculation results under different window sizes, when the current background pixel ratio is greater than the user-set threshold, the binary layers are merged to obtain the optimized patch units; Step 3: Calculate the image normalized difference spectral vector NDSV: traverse all the pairwise combinations of the reflectance of the bands, calculate the normalized difference index for each pair of band combinations, and obtain the corresponding normalized difference spectral vector NDSV. All normalized difference spectral vectors NDSV constitute the image spectral feature set. The fourth step is to use each patch unit as a clustering calculation unit, use the K-means algorithm for unsupervised clustering, and use the K-means++ algorithm to initialize the cluster center of each patch unit one by one: for any patch unit, randomly select a normalized difference spectral vector NDSV of a pixel as the cluster center of the first cluster. When there are k cluster centers, calculate the normalized difference spectral vector NDSV of other pixels in the current patch unit to the kth cluster center μ k The Euclidean distance is calculated, and the sampling weight of the corresponding pixel is set according to the Euclidean distance. The k+1th cluster center is randomly selected according to the sampling weight, and the process is repeated until the number of clusters meets the user's preset requirements; Step 5: Use the variance ratio criterion VRC to optimize the number of clusters for each patch unit: define a set of cluster numbers, for any patch unit, input the normalized difference spectral vector NDSV of all pixels into the K-means algorithm for clustering, calculate the variance ratio criterion VRC result based on the clustering result, and obtain the cluster number K corresponding to the maximum variance ratio criterion VRC by traversing all elements in the cluster number set. j As the optimal number of clusters for clustering within the current patch unit; Step 6. According to the methods in the fourth and fifth steps, the K-means method is performed on all patch units, and the global image pseudo sample set is obtained according to the clustering result of the optimal number of clusters: the spectral image feature set obtained in the third step is used as the K-means input feature, and all patch units obtained in the second step are traversed. The cluster center is obtained for all pixels in each patch unit using the method in the fourth step, and the optimal number of clusters is obtained using the method in the fifth step. On this basis, the clustering result is obtained, and the cluster with the largest number of pixels in the patch unit inherits the category information of the patch unit to obtain the pseudo sample set; Step 7. Determine the land cover classification system based on the land cover product, select the corresponding key features for each land category, and construct multiple two-class Gaussian mixture models: calculate the improved normalized water index MNDWI of the global image as the key feature of the water body type, calculate the normalized vegetation index NDVI of the non-water area as the key features of the forest, grassland, and artificial surface types, and calculate the slope of the vegetation area as the key feature of the cultivated land type. The characteristic value distribution of each key feature is regarded as a mixture of two Gaussian distributions, and the following Gaussian mixture model is constructed: in is the probability density function, y represents the eigenvalue of the key feature, α 1 , α 2 are the mixing coefficients of the first and second Gaussian distributions, μ 1 , μ 2 are the mean vectors of the two distributions, ∑ 1 ,∑ 2 The distribution is the covariance vector of the two distributions, p(y|μ 1 ,∑ 1 ) and p(y|μ 2 ,∑ 2 ) are the first and second Gaussian distribution density functions respectively. j ,μ j ,∑ j )|j=1,2} are collectively referred to as mixing component parameters; Step 8: Use the iterative expectation maximization EM algorithm to estimate the mixed component parameters in each Gaussian mixture model, complete the decomposition of the Gaussian mixture model, and thus complete the optimization of the pseudo sample set to obtain the training sample set; Step 9. Use change detection and classification methods to complete the surface cover reconstruction of the target historical phase: determine whether the target phase and the surface cover product belong to the same year. If they belong to the same year, input the optimized training samples in step 8 into the random forest RF classifier to complete the training and obtain the surface cover classification result of the target phase; if they belong to different years, use the dual-phase image change detection method to obtain the invariant area, transfer the sample labels of the invariant area, combine the target phase image to form training samples, and input the training samples into the random forest RF classifier to complete the training and obtain the surface cover classification result of the target historical phase.

2. The method for rapid reconstruction of historical land cover based on GlobeLand30 according to claim 1, Features: In the land patch optimization method described in the second step, the shape of the morphological structure element is a square, the domain size range is [3,21], the step size is 2, and the background pixel ratio ξ is set to a threshold of 5%.

3. The method for rapid reconstruction of historical land cover based on GlobeLand30 according to claim 1, Features: The calculation formula of the normalized difference spectral vector NDSV feature described in the third step is as follows: where y d Represents the spectral feature vector of the dth pixel in the image, The bth pixel is the dth pixel. i The reflectivity of each band, The bth pixel is the dth pixel. j The reflectivity of each band, the value range of i and j is [1, B].

4. The method for rapid reconstruction of historical land cover based on GlobeLand30 according to claim 1, Features: The variance ratio criterion described in the third step sets the number of clusters to a natural number ranging from 1 to 10.

5. The method for rapid reconstruction of historical land cover based on GlobeLand30 according to claim 1, Features: In the sixth step, the clustering unit of the improved K-means algorithm is set to the patch unit in the second step, and the clustering feature is set to the normalized difference spectral vector NDSV feature calculated in the third step.

6. The method for rapid reconstruction of historical land cover based on GlobeLand30 according to claim 1, Features: In the eighth step, by setting a threshold for the change in the likelihood function between two times in the iterative expectation maximization EM algorithm, when the likelihood function no longer changes significantly, the iteration process is stopped, the decomposition of the Gaussian mixture model is completed, the pseudo sample set is removed from the wrong label, and the training sample set is obtained; the threshold for the change in the likelihood function is set to ΔLL≤10 -6 .

7. The method for rapid reconstruction of historical land cover based on GlobeLand30 according to claim 1, Features: In the ninth step, the parameters of the random forest RF classifier are set to the number of decision trees Ntree is set to 500, and the feature subset dimension Mtry of the base classifier is set to the square root of the feature set dimension of the input classifier.

8. The method for rapid reconstruction of historical land cover based on GlobeLand30 according to claim 1, Features: In the ninth step, the change detection algorithm is set to perform threshold segmentation using the Otsu method based on the change vector analysis intensity image results.

Images