## 1 Introduction
Among them, studies on using geostatistics to improve the accuracy of land cover classification products are limited. [PERSON] [25] used the sequential indicator simulation algorithm and the cokriging method to predict the area range of olive trees by using the classified image as the soft data and the interpretation sample from the aerial image as the hard index. [PERSON] et al. [14] used indicator Kriging to estimate the spatial variation in the accuracy of the source land cover products and use integration methods, which consider the local accuracy of each source product to obtain Africa's land cover products. [PERSON] et al. [26] improved the accuracy of land cover classification by using a direct sequential co-simulation algorithm combined with field observations and remote sensing images classified by maximum likelihood classification. [PERSON] et al. [27] utilized a multipoint geostatistic calculation method to take the maximum likelihood classification result as the training image, the maximum likelihood probability as the soft condition, and the training sample as the hard condition and improve the accuracy of classification products. [PERSON] and [PERSON] [28] proposed to use the Markov chain random field (MCRF) method to evaluate the uncertainty of the land cover classification. They carried out conditional random simulation at the sample pixels interpreted by experts referring to high-resolution images and other auxiliary information. Unlabeled pixels are regarded as uncertain regions, and the class is obtained by using the MCRF algorithm [29]. However, the disadvantage of this method is that a large number of expert-interpreted samples are needed to ensure the reliability of results; consequently, the method becomes time-consuming and laborious. [PERSON] et al. [30] further improved the algorithm by introducing a pre-classified image to the simulation process. The co-simulation algorithm of the MCRF (Co-MCSS) is used to improve the classification accuracy by integrating sample pixels and pre-classified images. In this study, the boundary of geostatistics is confined in each Geo-Eco zone, and the Co-MCSS method [30] is used to quantify the classification local accuracy by integrating the pre-classified image, the verification points, and the attribute data based on the Geo-Eco zoning. Three contributions are made in terms of determining, evaluating, and improving the accuracy of the land cover product. 1. Taking Geo-Eco zoning as the area of geostatistics. Generally the area of geostatistics is based on the selected image area, which often contains boundary effect, that is, the fact that a class may have statistically biased smaller frequencies of transitions if it has a higher chance of occurring at boundaries of the study area because boundary polygonsare incomplete and have no transition to other classes beyond the boundary [28]. Generally, the boundary of Geo-Eco zoning does not cross two different types of land cover to avoid the boundary effect. 2. In previous research, the large number of sample pixels were interpreted by experts to obtain a transiogram model [28], but this process is time-consuming and laborious. In this study, reliable verification data published by some websites or institutions related to land cover research are reused to reduce the work on visual interpretation. 3. The Co-MCSS method not only combines remote sensing pre-classified images but also takes various attributes of Geo-Eco zoning (such as DEM, slope, temperature, and humidity) as auxiliary data to participate in the simulation and calculation of cross field transition probability. With the combination of additional attribute information, the simulation algorithm becomes more robust.