CN114612794B - A remote sensing identification method for land cover and planting structure in fragmented agricultural areas - Google Patents

Abstract

The invention belongs to the field of agricultural remote sensing, and relates to a remote sensing identification method of a ground covering and planting structure of a finely divided agricultural area. According to the invention, the agricultural area and the non-agricultural area are classified by the support vector machine SVM of the supervision algorithm, and then the crop type is further classified and judged by adopting the decision tree algorithm based on the obtained agricultural area classification raster data. The method does not depend on measured data, greatly reduces classification time and economic cost, and effectively avoids negative influence of the quality of the sample training set of the traditional machine algorithm on classification results. Meanwhile, the method meets the land utilization drawing precision required by agricultural remote sensing application on the premise of ensuring the classification high efficiency.

Description

A remote sensing identification method for land cover and planting structure in fragmented agricultural areas

Technical Field

The present invention belongs to the field of agricultural remote sensing, and relates to a remote sensing identification method for land cover and planting structure in fragmented agricultural areas. The method uses a machine learning algorithm to extract land use, does not rely on measured sample points, and obtains high-precision agricultural area land use rasterized data after automatic training.

Background technique

Accurate and timely acquisition of land cover information in agricultural areas is of great significance for monitoring agricultural conditions, adjusting agricultural structure, estimating production and formulating food policies. Traditional field survey methods are time-consuming and labor-intensive, and it is difficult to obtain large-scale regional land use information. Remote sensing technology can overcome the shortcomings of single-point observations with poor spatial representativeness under complex surface conditions, and provide multi-temporal, multi-spectral and multi-angle surface information for regional land use identification.

At present, the remote sensing inversion of land use is mainly divided into two methods: based on multi-source data phenology and polarization characteristics and based on measured sample point classification. The recognition algorithm based on sample points only needs to input the sample training set and use satellite remote sensing images to quickly realize large-scale regional land use inversion. Among them, algorithms such as support vector machine (SVM) and random forest are widely used in land use remote sensing identification because of their high training efficiency and accuracy that meets application requirements. However, this method has high requirements on sample accuracy and requires a large number of measured samples when applied, which is time and money costly. At the same time, for crops with similar spectral and texture characteristics, this type of algorithm will mix crop types during identification. The algorithm based on multi-source phenology and polarization characteristics only needs to establish specific classification rules, without the need for measured sample training sets, and can use the decision tree algorithm to realize land use inversion. This method is only for the identification of objects with very obvious differences in phenology/polarization characteristics, and it is necessary to process and analyze the original satellite image spectral data to calculate the phenology/polarization characteristics for classification. In summary, when the existing classification algorithms are applied to satellite remote sensing images, a single classification algorithm takes a long time to operate, has high requirements on computer configuration, requires a large number of measured samples, is time-consuming and labor-intensive, and therefore it is difficult to achieve large-scale land use identification.

In practical applications, various classifiers such as SVM, random forest, decision tree, etc. have their advantages and disadvantages. Therefore, the present invention combines the advantages of various algorithms and applies different algorithms according to the characteristics of different landforms in fragmented agricultural areas, taking advantage of their strengths and making up for their weaknesses, while improving the efficiency and accuracy of classification.

Summary of the invention

In view of the above technical problems, the purpose of the present invention is to provide a remote sensing identification method for land cover and planting structure in fragmented agricultural areas, which uses the spectral and texture characteristics of satellite remote sensing images, as well as the LSWI index, which is more effective in identifying agricultural and non-agricultural areas, to achieve the identification of agricultural and non-agricultural areas using a machine learning algorithm. On this basis, classification rules are established for the phenological characteristics of different crops in agricultural areas, and the phenological indicators calculated using NDVI time series data are used as input data. The further division of different crops in agricultural areas is achieved based on the decision tree algorithm, and finally high-precision land use remote sensing identification of fragmented agricultural areas is achieved.

In order to achieve the above object, the present invention provides the following technical solutions:

A remote sensing identification method for land cover and planting structure in fragmented agricultural areas comprises the following steps:

S1. Use the support vector machine (SVM) supervised learning algorithm to identify agricultural and non-agricultural areas in the original satellite images:

S1.1. Use bilinear interpolation to resample each band of multiple original satellite images of different periods and different meteorological conditions to a uniform spatial resolution; select an original satellite image of crops growing vigorously in early August, and select the visible light band that can better distinguish the ground objects and the red edge band, near infrared band or short-wave infrared band that can judge the growth status of plants in the image as the classification and identification bands; use the second-order probability statistics method to use a gray tone spatial correlation matrix, based on a 3×3 window and 64 gray quantization levels to calculate the texture value of the classification and identification band;

S1.2. Using the original satellite images under clear and cloudless conditions for multiple periods during the growing season, the land water index LSWI time series data set is obtained through formula 1. The land water index LSWI time series data set is resampled to a uniform time resolution using a linear interpolation method, and the land water index LSWI time series data in the data set are reassigned with a preset threshold. The pixels greater than the preset threshold are assigned a value of 0, and the pixels less than the preset threshold are assigned a value of 1. The reassigned land water index LSWI time series data are added and summed in sequence to obtain a land water index LSWI reclassification data;

Where: LSWI is the land water index; ρ _NIR represents the reflectivity of the near-infrared band; ρ _SWIR represents the reflectivity of the short-wave infrared band;

S1.3, select a land use sample set from Google Earth as training data, use the support vector machine (SVM) supervised classification algorithm to train and learn the classification identification bands and texture values obtained in step S1.1 and the LSWI reclassification data obtained in step S1.2, and calculate the land use classification results of agricultural and non-agricultural areas;

S2. Apply the decision tree algorithm to further divide the planting structure of the agricultural area and obtain the final high-precision land use data of the agricultural area:

S2.1, converting the land use classification result obtained in step S1 into raster data and re-assigning it, assigning the agricultural area to 1 and the non-agricultural area to NoData; converting the re-assigned raster data into vector data to obtain shp vector data containing only the agricultural area;

S2.2, based on the shp vector data containing only the agricultural area obtained in step S2.1, the original satellite images under clear and cloudless conditions for multiple periods during the growth period are cropped to obtain the original satellite images containing only the agricultural area, and the normalized vegetation index NDVI time series data set is calculated by formula 2, and the normalized vegetation index NDVI time series data set is resampled to a uniform time resolution by using a linear interpolation method;

Where: NDVI is the normalized vegetation index; ρ _NIR represents the reflectance of the near-infrared band; ρ _Red represents the reflectance of the red light band;

S2.3, using Savitzky-Golay filter to smooth the normalized vegetation index NDVI time series data set after resampling in step S2.2; establishing a phenological index judgment standard based on the smoothed NDVI growth period process curve, judging each pixel in the agricultural area, and finally obtaining phenological index raster data; the phenological index raster data includes three phenological indicators: the start time of the growth period, the end time of the growth period and the length of the growth period;

The process of establishing the phenological index judgment standard is as follows:

Determine the rising stage, falling stage and NDVI set value of the smoothed NDVI growth period process curve. When NDVI = set value, the moment corresponding to the rising stage is the start time of the growth period, the moment corresponding to the falling stage is the end time of the growth period, and the time length between the start time of the growth period and the end time of the growth period is the length of the growth period.

S2.4. Establish crop classification rules according to the value ranges of the three phenological indicators of different crops; based on the crop classification rules, use the decision tree algorithm to judge each grid in the phenological indicator grid data obtained in step S2.3, and obtain the agricultural planting structure classification results, so that each grid is finally divided into a specific crop type, thereby realizing crop identification in the agricultural area;

S3. Inlay the non-agricultural areas in the land use classification result obtained in step S1 with the agricultural planting structure classification result obtained in step S2 to obtain a land use and planting structure classification map of the agricultural planting area with classification label colors.

In step S1.2, the preset threshold is 0.2.

In the step S1.3, in the land use classification result, the agricultural area is farmland, and the non-agricultural area includes water bodies, sand dunes, residential areas and natural land.

In step S2.3, the setting value is 0.3.

The agricultural planting structure classification results obtained in step S2.4 include wheat, corn, sunflower, and fruits and vegetables.

In step S2.4, the value ranges of the three phenological indicators of different crops are determined based on local practical experience and field research.

In step S3, each pixel in the agricultural planting area land use and planting structure classification map has a classification, and different values and colors represent different classifications.

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

1. The present invention solves the problem that the number of field samples is limited, the quality requirements of measured samples are high, and they are difficult to obtain. The classification accuracy and stability are guaranteed without the need for measured samples.

2. Compared with traditional algorithms, the present invention reduces time and money costs and has strong operability. When applied to large-area high-precision remote sensing images, it can also use ordinary hardware equipment to achieve high-precision land use identification, making it possible to automatically obtain land use in large agricultural areas.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG1 is a flow chart of a remote sensing identification method for land cover and planting structure in fragmented agricultural areas according to the present invention;

FIG2 is an original satellite remote sensing image according to an embodiment of the present invention (Sentinel-2 historical images are used as an example in this description, and other satellite original images may also be used);

3 is a schematic diagram of a process for calculating texture values of an original satellite remote sensing image according to an embodiment of the present invention;

FIG4 is a diagram showing the calculation process of the supervised algorithm support vector machine (SVM) classification and the land use image obtained after classification according to an embodiment of the present invention;

FIG5 is a schematic diagram of phenological indicators selected by the present invention;

FIG6 is a decision tree classification rule for crop classification according to an embodiment of the present invention;

FIG7 is an image of an agricultural area after specific planting structure classification according to an embodiment of the present invention;

FIG8 is a land use and planting structure classification diagram finally obtained in an embodiment of the present invention;

FIG9 is a legend of each category in the classification results.

Detailed ways

The present invention will be further described below in conjunction with the accompanying drawings and embodiments.

As shown in FIG1 , a remote sensing identification method for land cover and planting structure in a fragmented agricultural area includes the following steps:

S1. Use the support vector machine (SVM) supervised learning algorithm to identify agricultural and non-agricultural areas in the original satellite images:

S1.1. Use bilinear interpolation to resample each band of multiple original satellite images of different periods and different meteorological conditions to a uniform spatial resolution; select an original satellite image of crops growing vigorously in early August, and select the visible light band that can better distinguish the ground objects and the red edge band, near infrared band or short-wave infrared band that can judge the growth status of plants in the image as the classification and identification bands; use the second-order probability statistics method to use a gray tone spatial correlation matrix, based on a 3×3 window and 64 gray quantization levels to calculate the texture value of the classification and identification band;

S1.2. Using the original satellite images under clear and cloudless conditions for multiple periods during the growing season, the land water index LSWI time series data set is obtained through formula 1. The land water index LSWI time series data set is resampled to a uniform time resolution using a linear interpolation method, and the land water index LSWI time series data in the data set are reassigned with a preset threshold. The pixels greater than the preset threshold are assigned a value of 0, and the pixels less than the preset threshold are assigned a value of 1. The reassigned land water index LSWI time series data are added and summed in sequence to obtain a land water index LSWI reclassification data;

Where: LSWI is the land water index; ρ _NIR represents the reflectivity of the near-infrared band; ρ _SWIR represents the reflectivity of the short-wave infrared band.

The preset threshold is 0.2.

S1.3, select a land use sample set from Google Earth as training data, use the support vector machine (SVM) supervised classification algorithm to train and learn the classification identification bands and texture values obtained in step S1.1 and the LSWI reclassification data obtained in step S1.2, and calculate the land use classification results of agricultural and non-agricultural areas;

In the land use classification results, agricultural areas are farmlands, and non-agricultural areas include water bodies, sand dunes, residential areas and natural land.

S2. Apply the decision tree algorithm to further divide the planting structure of the agricultural area and obtain the final high-precision land use data of the agricultural area:

S2.1, converting the land use classification result obtained in step S1 into raster data and re-assigning it, assigning the agricultural area to 1 and the non-agricultural area to NoData; converting the re-assigned raster data into vector data to obtain shp vector data containing only the agricultural area;

S2.2, based on the shp vector data containing only the agricultural area obtained in step S2.1, the original satellite images under clear and cloudless conditions for multiple periods during the growth period are cropped to obtain the original satellite images containing only the agricultural area, and the normalized vegetation index NDVI time series data set is calculated by formula 2, and the normalized vegetation index NDVI time series data set is resampled to a uniform time resolution by using a linear interpolation method;

Wherein: NDVI is the normalized difference vegetation index; ρ _NIR represents the reflectance of the near infrared band; ρ _Red represents the reflectance of the red light band.

S2.3, using a Savitzky-Golay (S-G) filter to smooth the normalized vegetation index NDVI time series data set resampled in step S2.2; establishing a phenological index judgment standard based on the smoothed NDVI growth period process curve, judging each pixel in the agricultural area, and finally obtaining phenological index raster data; the phenological index raster data includes three phenological indicators: the start time of the growth period, the end time of the growth period, and the length of the growth period;

The process of establishing the phenological index judgment standard is as follows:

Determine the rising stage, falling stage and NDVI set value of the smoothed NDVI growing period process curve. When NDVI = set value, the moment corresponding to the rising stage is the starting time of the growing period, and the moment corresponding to the falling stage is the ending time of the growing period. The length of time between the starting time and the ending time of the growing period is the length of the growing period.

The set value is 0.3.

S2.4. Establish crop classification rules according to the value range of the three phenological indicators of different crops; based on the crop classification rules, use the decision tree algorithm to judge each grid in the phenological indicator raster data obtained in step S2.3 to obtain the agricultural planting structure classification results, and finally each grid is divided into a specific crop type, thereby realizing crop identification in agricultural areas.

The agricultural planting structure classification results obtained in step S2 include wheat, corn, sunflower, and fruits and vegetables.

The value ranges of the three phenological indices for different crops were determined based on local practical experience and field surveys.

S3. Inlay the non-agricultural areas in the land use classification result obtained in step S1 with the agricultural planting structure classification result obtained in step S2 to obtain a land use and planting structure classification map of the agricultural planting area with classification label colors.

Each pixel in the agricultural planting area land use and planting structure classification map has a classification, and different values and colors represent different classifications.

Example

S1. Use the support vector machine (SVM) supervised learning algorithm to identify agricultural and non-agricultural areas in the original satellite images;

Each band of the Sentinel-2 original satellite image was resampled to a uniform spatial resolution of 10m using bilinear interpolation, as shown in Figure 2. An original satellite image of crops growing vigorously in early August was selected, and the visible light (red, green, and blue) bands that can better distinguish the ground objects in the image, as well as the red edge band, near infrared band, and short-wave infrared band that can judge the growth status of plants were selected as classification and identification bands. The texture value of the optimal band was calculated based on a 3×3 window and 64 grayscale quantization levels using a second-order probability statistical method using a gray tone spatial correlation matrix, and finally the texture value calculated in sequence based on eight texture filters (mean, variance, synergy, contrast, dissimilarity, information entropy, second-order moment, and correlation) was obtained, as shown in Figure 3.

Using Sentinel-2 original satellite images under clear and cloudless conditions during the growing season, the LSWI time series data set of the land water index was calculated by formula 1. The LSWI data set of the land water index was resampled to a uniform 5-day time resolution using the linear interpolation method, and re-assigned with a threshold of 0.2. Pixels greater than 0.2 were assigned a value of 0, and pixels less than 0.2 were assigned a value of 1. The re-assigned LSWI time series data of the land water index were added and summed in turn to obtain a LSWI reclassification data of the land water index.

A land use sample set (residential area, sand dune, wasteland, water body and farmland) was manually selected from Google Earth as training data. The support vector machine (SVM) supervised classification algorithm was used to train and learn the optimal band and texture value as well as the LSWI reclassification data. The classification results of land use in agricultural and non-agricultural areas were calculated, as shown in Figure 4.

S2, applying the decision tree algorithm to further divide the planting structure of the agricultural area and obtain the final high-precision land use data of the agricultural area;

The land use classification results are converted into raster data and re-assigned, with the agricultural area assigned to 1 and the non-agricultural area uniformly assigned to NoData. The re-assigned raster data are converted into vector data to obtain the shp vector data of the agricultural area.

Based on the vector data of agricultural areas, the original satellite images under clear and cloudless conditions during the growth period were cropped to obtain the original satellite images containing only the agricultural areas. The normalized difference vegetation index NDVI time series data set was calculated by formula 2, and the time series data set was resampled to a unified 5-day time resolution using the linear interpolation method.

The Savitzky-Golay (S-G) filter was used to smooth the NDVI time series data set during the growth period. As shown in Figure 5, based on the smoothed NDVI growth period process curve, the rising and falling stages of the curve were determined. When NDVI = 0.3, the rising and falling stages corresponded to the start time and end time of the growth period, respectively. The length of time between the start time and the end time of the growth period was the length of the growth period. Based on the above phenological index judgment standard, each pixel in the agricultural area was judged, and finally the phenological index grid data containing the three phenological indicators of the start time, end time and length of the growth period were obtained for further crop classification.

Based on local actual experience and field research, the value ranges of the three phenological indicators of different crops were determined and crop classification rules were established, as shown in Figure 6. Based on the classification rules, the decision tree algorithm was used to judge each grid in the phenological indicator grid data, and finally each grid was divided into a specific crop type, thereby realizing crop identification in agricultural areas and obtaining the classification results of agricultural area planting structure as shown in Figure 7.

The non-agricultural areas in the land use results are inlaid with the planting structure identification results, as shown in Figure 8, and then the high-precision land use classification results of the agricultural area are obtained. The classification result is an image with classification label colors, including water bodies, sand dunes, residential areas, natural non-agricultural land use types, and wheat, corn, sunflower and fruit and vegetable planting structure classification results. Examples of various land features and crops are shown in Figure 9.

Claims (7)

1. The remote sensing identification method for the ground cover and planting structure of the finely divided agricultural area is characterized by comprising the following steps:

S1, carrying out agricultural region and non-agricultural region identification on an original satellite image by using a support vector machine SVM supervised learning algorithm:

S1.1, resampling each wave band of a plurality of original satellite images with different meteorological conditions in different periods into uniform spatial resolution by adopting a bilinear interpolation method; selecting an original satellite image with vigorous growth of the early-stage August crops, and selecting a visible light wave band capable of distinguishing ground features and a red-edge wave band, a near-infrared wave band or a short-wave infrared wave band capable of judging plant growth conditions in the image as classification recognition wave bands; calculating texture values of the classification and identification bands based on a 3×3 window and 64 gray quantization levels by using a gray tone space correlation matrix by using a second-order probability statistical method;

S1.2, obtaining a land water index LSWI time sequence data set by utilizing an original satellite image under a cloudless condition of a multi-period clear sky in a growth period, resampling the land water index LSWI time sequence data set into uniform time resolution by adopting a linear interpolation method, reassigning the land water index LSWI time sequence data in the data set by using a preset threshold value, assigning pixels larger than the preset threshold value to 0, assigning pixels smaller than the preset threshold value to 1, sequentially adding and summing the reassigned land water index LSWI time sequence data to obtain land water index LSWI reclassification data;

wherein: LSWI is land water index; ρ _NIR represents the reflectivity in the near infrared band; ρ _SWIR represents the reflectivity of the short wave infrared band;

S1.3, selecting a land utilization sample set from Google Earth as training data, utilizing a support vector machine SVM supervised classification algorithm to train and learn the classification recognition wave band and texture value obtained in the step S1.1 and land water index LSWI reclassification data obtained in the step S1.2, and calculating to obtain land utilization classification results of agricultural areas and non-agricultural areas;

S2, carrying out further planting structure division on the agricultural area by applying a decision tree algorithm to obtain final agricultural area high-precision land utilization data:

S2.1, converting land utilization classification results obtained in the step S1 into raster data, reassigning the raster data, assigning an agricultural area to 1, and assigning a non-agricultural area to NoData; converting the reassigned raster data into vector data to obtain shp vector data only containing an agricultural area;

S2.2, cutting an original satellite image under a multi-period clear sky cloudless condition in a growing period based on shp vector data only comprising an agricultural area obtained in the step S2.1 to obtain an original satellite image only comprising the agricultural area, calculating through a formula II to obtain a normalized vegetation index NDVI time sequence data set, and resampling the normalized vegetation index NDVI time sequence data set to uniform time resolution by adopting a linear interpolation method;

Wherein: NDVI is normalized vegetation index; ρ _NIR represents the reflectivity in the near infrared band; ρ _Red represents the reflectivity of the red band;

S2.3, smoothing the normalized vegetation index NDVI time sequence data set resampled in the step S2.2 by utilizing a Savitzky-Golay filter; establishing a weather indicator judgment standard based on the smoothed NDVI growth period process curve, judging each pixel of the agricultural area, and finally obtaining weather indicator grid data; the weather indicator grid data comprises three weather indicators including a growth period starting time, a growth period ending time and a growth period length;

The establishment process of the weather indicator judgment standard is as follows:

Determining an ascending phase, a descending phase and an NDVI set value of the smoothed NDVI fertility phase process curve, wherein when NDVI=set value, the moment corresponding to the ascending phase is a fertility phase starting time, the moment corresponding to the descending phase is a fertility phase ending time, and the time length in the fertility phase starting time and the fertility phase ending time is a fertility phase length;

S2.4, establishing crop classification rules according to the value range of three weather indexes of different crops; based on crop classification rules, judging each grid in the weather indicator grid data obtained in the step S2.3 by utilizing a decision tree algorithm to obtain an agricultural planting structure classification result, and finally dividing each grid into a certain specific crop type, thereby realizing crop identification of an agricultural area;

and S3, inlaying the non-agricultural area in the land utilization classification result obtained in the step S1 and the agricultural planting structure classification result obtained in the step S2 together to obtain a land utilization and planting structure classification diagram of the agricultural planting area with classification label colors.

2. The method for remotely sensing and identifying land cover and planting structures in finely divided agricultural areas according to claim 1, wherein in the step S1.2, the preset threshold is 0.2.

3. The method for remotely identifying land cover and planting structures in finely divided agricultural areas according to claim 1, wherein in the step S1.3, the agricultural area is a farmland, and the non-agricultural area includes water, sand dunes, residential areas and natural lands.

4. The method for remotely sensing and identifying land cover and planting structures in finely divided agricultural areas according to claim 1, wherein in the step S2.3, the set value is 0.3.

5. The method for remotely identifying land cover and planting structures in finely divided agricultural areas according to claim 1, wherein the classification result of agricultural planting structures obtained in step S2.4 comprises wheat, corn, sunflower and fruits and vegetables.

6. The method for remotely sensing and identifying the ground cover and planting structure of a finely divided agricultural area according to claim 1, wherein in the step S2.4, the value ranges of three physical indexes of different crops are respectively determined based on local actual experience and field investigation.

7. The method for remote sensing identification of finely divided agricultural land cover and planting structures according to claim 1, wherein in said step S3, each pixel in the classification map of agricultural land utilization and planting structures has a classification, and different values and colors represent different classifications.