Remote sensing identification method of irrigation area

Abstract

The present invention discloses a remote sensing identification method for irrigation area, including: S1, data preprocessing: collecting three sets of irrigation area data sets of the area to be identified and data sets related to the input features of the model, unifying the projection coordinate system and spatial resolution of each data set; S2, constructing irrigation possibility data based on mathematical uncertainty analysis; S3, constructing training samples based on statistical yearbooks and time filtering; S4, machine learning input feature data set processing; S5, training machine learning; S6. Irrigation area identification. The method enhances the ability to capture time changes, significantly improves the ability of machine learning in capturing time changes in irrigation areas by introducing continuous and discontinuous training samples, and can accurately identify in the case of sudden changes in irrigation areas; the method solves the problem of obtaining training samples; the method is applicable to any year, and improves the temporal continuity of irrigation area data sets and the applicability of irrigation area identification.

Description

Technical Field

[0001] The invention relates to the field of remote sensing technology, and in particular to a remote sensing identification method for irrigation areas. Background Art

[0002] Irrigation accounts for 70% of global freshwater consumption and is one of the most important human intervention activities in the terrestrial water cycle. It is also a key factor in ensuring global food security. In addition, irrigation also plays an important role in mitigating extreme events under climate change. Therefore, accurate simulation of irrigation activities is essential to achieve sustainable development of water resources and food security under the background of climate change.

[0003] However, the uncertainty of irrigation area data is one of the main constraints in irrigation simulation. The commonly used method is to use county-level statistical data and remote sensing land feature classification data (such as the Global Irrigation Area Dataset and the Global Irrigation Area Dataset) for downscaling. These datasets can provide reliable spatial distribution information of irrigation area, but due to the high cost of collecting statistical data, it is difficult to provide continuous data for many years, and thus it is impossible to provide information on the temporal changes of irrigation area.

[0004] Machine learning technology provides the possibility of pixel-level classification of irrigated and rain-fed farmland through remote sensing data, and can more efficiently identify long-term and large-scale irrigated area datasets. However, a major challenge in irrigated area identification based on machine learning is how to obtain accurate, reliable and sufficient training samples. Currently, the main means of obtaining high-quality sample data is manual labeling, which is time-consuming, labor-intensive and costly. To address this problem, some studies have proposed methods for automatically obtaining samples. These methods use irrigated area pixels that remain unchanged for many years as training samples, but the insufficient number of samples limits the classification ability of machine learning models.

[0005] Categorical triple collocation-based merging (CTCM) has been shown to improve the classification accuracy of binary classification datasets. Therefore, this method is expected to be used to generate reliable and sufficient training samples to support more accurate irrigated area identification. Summary of the invention

[0006] The purpose of the present invention is to provide an irrigation area remote sensing identification method that combines mathematical uncertainty analysis, automatic sampling and machine learning methods to provide reliable and sufficient training samples and solve the limitations of existing irrigation area datasets in capturing the temporal changes of irrigation area.

[0007] To this end, the present invention adopts the following technical solutions:

[0008] A remote sensing identification method for irrigation area comprises the following steps:

[0009] S1, data preprocessing:

[0010] Collect three sets of irrigation area datasets of the area to be identified and datasets related to the input features of the model, and unify the projection coordinate system and spatial resolution of each dataset;

[0011] S2, constructing irrigation possibility data based on mathematical uncertainty analysis, includes the following sub-steps:

[0012] S21, calculate whether the pixels in the irrigated area dataset are irrigated or not, there are:

[0013] x i =I+e i (1)

[0014] Where:

[0015] x i Indicates whether the pixel in the irrigation area dataset i is irrigated, i = 1 to 3;

[0016] I represents the true value of whether a pixel in the irrigated area dataset is irrigated, where the pixel irrigated or not is represented by -1 and +1, where -1 represents a non-irrigated pixel and +1 represents an irrigated pixel;

[0017] e i is the classification error;

[0018] S22, calculate the balance accuracy π i :

[0019] π i =0.5(A i +D i ) (2)

[0020] In the formula, A i represents the probability that a pixel in the irrigated area dataset is irrigated; D i represents the probability that a pixel in the irrigated area dataset is non-irrigated; the classification error e between different irrigated area datasets i are independent of each other, then the balance accuracy π and the covariance Q between any two sets of irrigation area data sets are positively correlated, as follows:

[0021] Q 12 =Cov(x1,x2)=f(I)(2π1-1)(2π2-1) (3)

[0022] Q 13 =Cov(x1,x3)=f(I)(2π1-1)(2π3-1) (4)

[0023] Q23 =Cov(x2,x3)=f(I)(2π2-1)(2π3-1) (5)

[0024] in:

[0025] f(I) is a function of whether a pixel in the irrigated area dataset is irrigated or not, which is unknown;

[0026] π i With v i Linear correlation, v i The covariance Q between the three irrigation area data sets preprocessed by S1 is calculated as follows:

[0027]

[0028] S23, obtain the weight w of the probability that the pixel in different irrigation area data sets is the true value i :

[0029]

[0030] Combined with the x u , and get the probability of pixel irrigation X p :

[0031] X p =w1x1+w2x2+w3x3 (10)

[0032] Where, X p The range is -1 to 1;

[0033] S3, construct training samples based on statistical yearbooks and time filtering:

[0034] Based on S2, combined with the regional irrigation area data in the statistical yearbook, we gradually checked from high to low and determined the regional sample (pixel) threshold X. T The irrigation area obtained is less than or equal to the yearbook data, based on the regional sample threshold X T Obtain potential training samples; then filter by time to select pixels that have been irrigated for many consecutive years as general continuous training samples, and pixels that have been irrigated only in a specified year as discontinuous training samples for the specified year;

[0035] S4, machine learning input feature dataset processing:

[0036] In the data sets related to each input feature after S1 preprocessing, the variables closely related to the irrigation area are selected as the input features of machine learning. The land remote sensing data set is used to perform cultivated land masking and screen out non-cultivated land objects to obtain the input data for machine learning.

[0037] S5, training machine learning:

[0038] The universal continuous training samples and the non-continuous training samples obtained in S3 and the input features obtained in S4 are input into the supervised learning model for training to generate a supervised learning model capable of identifying the irrigation area;

[0039] S6. Identification of irrigation area:

[0040] The supervised learning model obtained by S5 is used to obtain the spatiotemporally continuous irrigation area data through the input features of the area to be identified.

[0041] Preferably, in step S1, each data set is uniformly converted into the WGS84 projection coordinate system, and these data are resampled to a resolution of 500 meters.

[0042] Preferably, the supervised learning model used in step S3 is a random forest model.

[0043] The input features described in step S4 are normalized difference vegetation index NDVI, greenness index GI, precipitation P, potential evapotranspiration PET, daytime land surface temperature LST, agricultural suitability P-PET and digital elevation model DEM.

[0044] Compared with the prior art, the present invention has the following beneficial effects:

[0045] 1. The recognition method of the present invention enhances the ability to capture temporal changes: By introducing continuous and discontinuous training samples, the present invention significantly improves the ability of machine learning in capturing temporal changes in irrigation area, optimizes irrigation area recognition technology, and can provide more accurate recognition results. It can also perform accurate recognition in the case of sudden changes in irrigation area.

[0046] 2. The recognition method of the present invention solves the problem of obtaining training samples: Through the fusion algorithm based on classification triple correction, the present invention significantly improves the pain point of difficulty in obtaining machine learning training samples, reduces the difficulty of sample acquisition, and improves the quantity and quality of sample acquisition, thereby improving the remote sensing recognition and classification capabilities based on machine learning, and promoting the development of remote sensing recognition technology for irrigation areas.

[0047] 3. The identification method of the present invention is applicable to any year, improves the temporal continuity of the irrigation area data set, and can obtain irrigation area data sets for many consecutive years, which is conducive to the sustainable development of agriculture and land planning development.

[0048] 4. The identification method of the present invention improves the applicability of irrigation area identification and can be used to identify irrigation areas on a large scale from the whole country to the whole world, and finally obtains a large-scale, spatiotemporally continuous irrigation area data set. BRIEF DESCRIPTION OF THE DRAWINGS

[0049] Figure 1 It is a flow chart of the irrigation area remote sensing identification method of the present invention;

[0050] Figure 2 The process of constructing the universal continuous training samples and the non-continuous training samples of the specified year in the present invention;

[0051] Figure 3 This is the verification of the identification results of the present invention. DETAILED DESCRIPTION

[0052] In view of the current demand for irrigation area data in land surface models, earth system models and irrigation estimation, this paper selects three sets of common irrigation area data sets based on mathematical uncertainty analysis methods to construct continuous and discontinuous training samples. For variables closely related to irrigation area, the normalized difference vegetation index (NDVI), greenness index (GI), precipitation (P), potential evapotranspiration (PET), daytime land surface temperature (LST), agricultural suitability (P-PET) and digital elevation model (DEM) are selected as input features of the machine learning random forest model. By training the random forest model and performing post-processing, a large-scale, long-time series of irrigation area data sets are finally obtained.

[0053] The irrigation area remote sensing identification method of the present invention is described in detail below with reference to the accompanying drawings and implementation examples.

[0054] Example 1

[0055] refer to Figure 1 The irrigation area remote sensing identification method of the present invention comprises the following steps:

[0056] S1, data preprocessing:

[0057] The present invention needs to use three sets of irrigation area data sets to obtain training samples in the recognition process. These three sets of data sets come from existing data, such as Syn (https: / / figshare.com / articles / dataset / The_500-m_irrigated_cropland_maps_in_China_during_2000-2019_based_on_a_synergy_mapping_method / 19352501 / 1?file=34373714), CN (https: / / figshare. com / articles / dataset / IrriMap_CN_Improved_annual_irrigation_maps_across_China_in_2000_2019_based_on_satellite_imagery_and_machine-learning_method / 20363115) and GIA (https: / / doi.pangaea.de / 10.1594 / PANGAEA.884744) data, requiring that one of the datasets has at least 5 years of data. In this example, both Syn and CN datasets used 15 years of data.

[0058] At the same time, machine learning requires the use of some input feature data, mainly from remote sensing data and existing data, such as NDVI annual data with a spatial resolution of 500 meters from MOD13A1, GI annual data set with a spatial resolution of 500 meters from MOD09A1, monthly precipitation data with a spatial resolution of 1 km from existing data sets (https: / / zenodo.org / records / 3114194), potential evapotranspiration monthly data with a spatial resolution of 1 km (https: / / data.tpdc.ac.cn / zh-han / data / 8b11da09-1a40-4014-bd3d-2b86e6dccad4), annual daytime land surface temperature data with a spatial resolution of 1 km from MOD11A2, digital elevation model data with a spatial resolution of 250 meters from the Chinese Academy of Sciences (https: / / www.resdc.cn / data.aspx?DATAID=123), and AI The annual data of the land cover dataset (CLCD) with a spatial resolution of 30 meters on the Earth platform (https: / / engine-aiearth.aliyun.com / # / dataset / CHINA_CLCD_V01_WHU_IRSIP), the annual land use data with a spatial resolution of 30 meters from the National Geographic Information Resources Center (https: / / www.webmap.cn / commres.do?method=globeIndex), and the annual land use data with a spatial resolution of 1 km from the National Tibetan Plateau Data Center.

[0059] Since the present invention needs to perform pixel-by-pixel recognition, a unified projection coordinate system and spatial resolution are required.

[0060] In this example, the above data sets and input feature data are unified into the WGS84 projection coordinate system; the data resolutions used in this embodiment are concentrated at 30 meters, 250 meters, 500 meters and 1 kilometer, and the data are resampled to a resolution of 500 meters using the nearest value method to ensure the consistency and comparability of the data.

[0061] S2, constructing irrigation possibility data based on mathematical uncertainty analysis: The three sets of irrigation area datasets Syn, CN and GIA preprocessed by S1 are used as input data for the classification triple correction fusion method based on mathematical uncertainty analysis to obtain the possibility of irrigation pixel by pixel.

[0062] In scientific measurements and experiments, the simulation value can usually be expressed as the sum of the true value and the error. Specifically for the irrigation area concerned by the present invention, the simulation value formula is as follows:

[0063] x i =I+e i(1)

[0064] In the formula, x i Indicates whether the pixel i in the irrigation area dataset is irrigated, i = 1 to 3; I represents the true value (real value) of whether the pixel in the irrigation area dataset is irrigated, where irrigation or not is represented by -1 and +1, -1 represents non-irrigation pixels, and +1 represents irrigation pixels; e i is the classification error.

[0065] Analog value x i It is related to the balanced accuracy. The irrigation problem that this invention focuses on is a binary classification problem, so the balanced accuracy also evaluates the ability of the data to capture whether the pixel is irrigated or not. The balanced accuracy formula is as follows:

[0066] π i =0.5(A i +D i ) (2)

[0067] In the formula, A i represents the probability that a pixel in the irrigated area dataset is irrigated (+1); D i Represents the probability that a pixel in the irrigated area dataset is non-irrigated (-1). Classification error between different irrigated area datasets e i are independent of each other, then the balanced accuracy π is positively correlated with the covariance Q between any two sets of irrigation area data sets:

[0068] Q 12 =Cov(x1,x2)=f(I)(2π1-1)(2π2-1) (3)

[0069] Q 13 =Cov(x1,x3)=f(I)(2π1-1)(2π3-1) (4)

[0070] Q 23 =Cov(x2,x3)=f(I)(2π2-1)(2π3-1) (5)

[0071] where f(I) is a function of whether the pixel is irrigated or not, which is usually unknown. Therefore, π i It cannot be solved directly by formulas (3)-(5). Instead, the present invention solves i Linearly related v i To obtain indirectly,

[0072]

[0073] The covariance Q between the three input data sets can be used to obtain v i, and then combine the following formula to obtain the weight w of the probability that the pixel in the different irrigation area data sets is the true value i :

[0074]

[0075] Combined with the x i , and obtain the probability X of pixel irrigation in the irrigation area dataset p :

[0076] X p =w1x1+w2x2+w3x3 (10)

[0077] Where, X p The range is from -1 to 1, with values ​​closer to 1 indicating a higher probability of irrigation and vice versa.

[0078] S3, construct training samples based on statistical yearbooks and time filtering:

[0079] After S2, the irrigation possibility result of each pixel in the irrigation area dataset is obtained. Combined with the regional irrigation area data in the statistical yearbook, the regional sample threshold X is determined by stepping through the investigation from high to low. T The irrigation area obtained is less than or equal to the yearbook data, based on the regional sample threshold X T Potential training samples are obtained; then, through time filtering, pixels that have been irrigated for many consecutive years are selected as general continuous training samples, and pixels that have been irrigated only in a specified year are selected as discontinuous training samples in the specified year.

[0080] See also Figure 2 In this embodiment, the pixel-by-pixel irrigation probability data of region A from 2000 to 2002 is obtained based on the classification triple correction fusion method, and the sample threshold is determined in combination with the irrigation area results of region A from 2000 to 2002 in the statistical yearbook.

[0081] For example, the irrigated area of ​​region A in the statistical yearbook in 2000 was 1.3 km 2 , which is approximately equal to 5.2 pixels (after S1 processing, each pixel size is 500 meters × 500 meters). T =0.8, the number of samples is less than or equal to the data in the statistical yearbook, so the annual sample threshold is determined to be 0.8. Through time filtering, pixels that have remained unchanged for many years are selected. These pixels that have remained unchanged for many years are universal training samples and can participate in the training of other years. At the same time, non-continuous training samples that only appear in the specified year can also be obtained to capture mutation information in machine learning training.

[0082] The above method is used to obtain non-continuous and continuous training samples from 2000 to 2015, which provide a training basis for S5.

[0083] S4, machine learning input feature data set processing: Based on the pre-processed data obtained in S1, the variables closely related to the irrigation area are selected as the input features of machine learning, and the cultivated land mask is generated using the land cover data set to screen out non-cultivated land objects to improve the operation efficiency. At the same time, when constructing the relationship equation between the training samples and the input data, it is necessary to select variables closely related to the irrigation area to establish a more stable recognition relationship.

[0084] In this example, the normalized vegetation index (NDVI), greenness index (GI), precipitation (P), potential evapotranspiration (PET), daytime land surface temperature (LST), agricultural suitability (P-PET) and digital elevation model (DEM) are selected from the remote sensing data preprocessed by S1 and the existing data as input features. At the same time, based on S1, the cultivated land mask constructed by the classification of land objects is used to preliminarily screen out non-cultivated land areas to improve the calculation efficiency of the present invention.

[0085] S5, training machine learning: select the random forest model in machine learning as the specific implementation model of this example, and input the general continuous training samples and non-continuous training samples obtained in S3 and the input features obtained in S4 into the supervised learning model for training. In this embodiment, the random forest model is preferably used to finally generate a random forest model that can identify the irrigation area.

[0086] S6. Irrigation area identification: Using the random forest model obtained in S5, the spatiotemporally continuous irrigation area data can be obtained through the input features of the area to be identified.

[0087] The random forest model is applicable to any year.

[0088] In this embodiment, the irrigation area in the whole country from 2000 to 2022 is identified by processing the relevant data of S1 and S2 in combination with the random forest model obtained by S5.

[0089] The above results were verified using the national county-level statistical yearbook data from 2000 to 2019. Six counties were randomly selected across the country and time series were plotted. The results are shown in Figure 3 In the figure, ▲ is the measured result of the statistical yearbook, and * is the result obtained by the present invention. It can be seen from the figure that the method of the present invention can accurately identify the irrigation area information, and the method of the present invention can still accurately capture the A, B and E areas with obvious mutation phenomena.

[0090] Compared with the traditional method of obtaining training samples only by time filtering, the present invention can accurately and robustly expand the scale of training samples and reduce the spatial distribution deviation of irrigation area from 16.64% to 8.27%. Most importantly, the present invention provides a way to use continuous and non-continuous training samples at the same time, so as to capture the sudden change of irrigation area in time and improve the irrigation area time correlation from 0.15 to 0.54.

Claims

1. A remote sensing identification method for irrigation area, characterized in that The following steps are involved: S1, data preprocessing: Collect three sets of irrigation area datasets of the area to be identified and datasets related to the input features of the model, and unify the projection coordinate system and spatial resolution of each dataset; S2, constructing irrigation possibility data based on mathematical uncertainty analysis, includes the following sub-steps: S21, calculate whether the pixels in the irrigated area dataset are irrigated or not, there are: x i =I+e i (1) Where: x i Indicates whether the pixel in the irrigation area dataset i is irrigated, i = 1 to 3; I represents the true value of whether a pixel in the irrigated area dataset is irrigated, where the pixel irrigated or not is represented by -1 and +1, where -1 represents a non-irrigated pixel and +1 represents an irrigated pixel; e i is the classification error; S22, calculate the balance accuracy π i : π i =0.5(A i +D i ) (2) In the formula, A i represents the probability that a pixel in the irrigated area dataset is irrigated; D i represents the probability that a pixel in the irrigated area dataset is non-irrigated; the classification error e between different irrigated area datasets i are independent of each other, then the balance accuracy π and the covariance Q between any two sets of irrigation area data sets are positively correlated, as follows: Q 12 =Cov(x1,x2)=f(I)(2π1-1)(2π2-1) (3) Q 13 =Cov(x1,x3)=f(I)(2π1-1)(2π3-1) (4) Q 23 =Cov(x2,x3)=f(I)(2π2-1)(2π3-1) (5) in: f(I) is a function of whether a pixel in the irrigated area dataset is irrigated or not, which is unknown; π i With v i Linear correlation, v i The covariance Q between the three irrigation area data sets preprocessed by S1 is calculated as follows: S23, obtain the weight w of the probability that the pixel in different irrigation area data sets is the true value i : Combined with the x i , and obtain the probability X of pixel irrigation in the irrigation area dataset p : X p =w1x1+w2x2+w3x3 (10) Where, X p The range is -1 to 1; S3, construct training samples based on statistical yearbooks and time filtering: Based on S2, combined with the regional irrigation area data in the statistical yearbook, we gradually checked from high to low to determine the regional sample threshold X T The irrigation area obtained is less than or equal to the yearbook data, based on the regional sample threshold X T Obtain potential training samples; then filter by time to select pixels that have been irrigated for many consecutive years as general continuous training samples, and pixels that have been irrigated only in a specified year as discontinuous training samples for the specified year; S4, machine learning input feature dataset processing: In the data sets related to each input feature after S1 preprocessing, the variables closely related to the irrigation area are selected as the input features of machine learning. The land remote sensing data set is used to perform cultivated land masking and screen out non-cultivated land objects to obtain the input data for machine learning. S5, training machine learning: The universal continuous training samples and the non-continuous training samples obtained in S3 and the input features obtained in S4 are input into the supervised learning model for training to generate a supervised learning model capable of identifying the irrigation area; S6. Identification of irrigation area: The supervised learning model obtained by S5 is used to obtain the spatiotemporally continuous irrigation area data through the input features of the area to be identified.

2. The irrigation area remote sensing identification method according to claim 1, characterized in that: In S1, each dataset was uniformly converted to the WGS84 projection coordinate system, and these data were resampled to a resolution of 500 m.

3. The method for remote sensing identification of irrigation area according to claim 1, characterized in that: The supervised learning model used in S3 is the random forest model.

4. The method for remote sensing identification of irrigation area according to claim 1, characterized in that: In S4, the input features are normalized difference vegetation index NDVI, greenness index GI, precipitation P, potential evapotranspiration PET, daytime land surface temperature LST, agricultural suitability P-PET and digital elevation model DEM.

Images