Extraction method of cultivated land non-agriculturalization by integrating multi-dimensional and multi-scale remote sensing and terrain features

Abstract

The present invention discloses a method for intelligent monitoring of non-agriculturalization of cultivated land by integrating multi-dimensional multi-scale remote sensing and terrain features, including step 1, obtaining high-resolution images, hyperspectral images, DEM data, true vector data of cultivated land, and vector data of cultivated land that has been investigated and monitored; step 2, segmenting and clustering superpixel objects; step 3, extracting features within the range of each superpixel object object with superpixel objects as units; step 4, combining the features extracted from the training area into feature vectors of each superpixel object, and constructing superpixel cultivated land object training sample data; step 5, training a cultivated land object recognition model; step 6, extracting cultivated land in a test area; step 7, based on the extracted cultivated land in the test area and the cultivated land vector data that has been investigated and monitored, using spatial overlay analysis to automatically extract the spots of cultivated land becoming non-agricultural land, optimizing the spots, filtering and screening the non-agriculturalized spots of cultivated land. The present invention has the advantages of practicality, convenience, feasibility, etc.

Description

Technical Field

[0001] The present invention relates to the field of remote sensing image feature extraction, and in particular to a method for extracting cultivated land non-agriculturalization by integrating multi-dimensional and multi-scale remote sensing with terrain features. Background Art

[0002] my country has a large population and limited land, and arable land is one of the most precious resources. In recent years, we have severely punished the illegal occupation of arable land, strictly controlled the conversion of arable land to other agricultural land, and strengthened the control of arable land use. Therefore, using remote sensing technology to accurately interpret arable land and realize intelligent monitoring of arable land non-agriculturalization plays an important role in arable land protection and food security.

[0003] At present, the remote sensing monitoring methods of cultivated land non-agriculturalization are mainly divided into two types: visual interpretation and automated algorithms. Remote sensing visual interpretation of cultivated land non-agriculturalization monitoring is to find cultivated land non-agriculturalization spots based on the previous and later images and the manual visual whole-domain search. It requires certain prior knowledge and has huge manpower and time costs. In the era of remote sensing big data, facing high-frequency and large-scale cultivated land protection monitoring, visual interpretation obviously cannot meet the actual needs. The automated cultivated land non-agriculturalization monitoring algorithm uses remote sensing intelligent recognition technology to automatically identify cultivated land and cultivated land change areas from images, mainly including two methods: intelligent change detection of cultivated land non-agricultural areas and intelligent interpretation of cultivated land before finding non-agricultural areas. Intelligent change detection of cultivated land non-agricultural areas requires the production of a large number of high-quality cultivated land non-agricultural change samples, the construction of a change detection model, and finally the use of the trained model to automatically find cultivated land non-agricultural spots. This method relies on the collection and production of a large number of high-quality cultivated land non-agricultural change samples. The entire cycle is long and it is difficult to meet the urgent needs of actual business. Intelligent interpretation of cultivated land first and then finding non-agricultural areas is done by interpreting the early and late images to extract cultivated land, and then judging the changes in cultivated land and non-agriculturalization, or just interpreting the cultivated land in the late images, and then using the cultivated land vectors of historical results to obtain non-agricultural patches.

[0004] Especially in the southwestern mountainous areas, the terrain is undulating and the vertical differentiation of the surface is serious. The characteristics of cultivated land in the southwestern mountainous areas are changeable and small and scattered in remote sensing images. The characteristics of cultivated land are similar to those of forest and grassland, plantations, and ponds covered with aquatic plants. The problems of "same object different spectrum" and "different objects same spectrum" are serious, resulting in the low accuracy of the current mainstream automatic recognition method of cultivated land using only high-resolution images. From the current method and technology, it is difficult to accurately extract cultivated land in complex mountainous environments by using only the spatial characteristics of high-resolution images. Summary of the invention

[0005] In view of the shortcomings of the existing technology, the present invention provides a method for intelligent monitoring of non-agricultural conversion of cultivated land by integrating multi-dimensional and multi-scale remote sensing with terrain features, which solves the difficult problems of accurate identification of cultivated land in complex mountainous environments and intelligent monitoring of non-agricultural conversion of cultivated land. The monitoring results of this method are more accurate and are practical, convenient and feasible.

[0006] To achieve the above object, the technical solution adopted by the present invention is as follows:

[0007] An intelligent monitoring method for non-agriculturalization of cultivated land by integrating multi-dimensional and multi-scale remote sensing with terrain features comprises the following steps:

[0008] Step 1, divide the training area and the test area, obtain high-resolution images, hyperspectral images, digital elevation model DEM data, and true vector data of cultivated land in the training area, and obtain high-resolution images, hyperspectral images, digital elevation model DEM data, and vector data of cultivated land that has been investigated and monitored in the test area;

[0009] Step 2: firstly, the high-resolution image is segmented into superpixels using the simple non-iterative clustering (SNIC) algorithm embedded with edge detection Canny, and then the superpixel objects are clustered using the Gaussian mixture algorithm (GMM) in combination with the high-resolution image to obtain superpixel ground objects;

[0010] Step 3, taking the superpixel ground object obtained in step 2 as a unit, extracting features within the range of each superpixel ground object, wherein the features at least include spatial domain features of high-resolution images, frequency domain high-frequency features of high-resolution images, spectral domain features of hyperspectral images, and terrain features of DEM data;

[0011] Step 4: Combine the high-resolution image, high-resolution image, DEM data and true value vector data of the training area into a feature vector of each super-pixel object according to the features extracted in step 2 and step 3, determine the super-pixel object label according to the pixel area with the largest category label in the super-pixel object, and construct super-pixel cultivated land object training sample data;

[0012] Step 5: input the super-pixel cultivated land object training sample data constructed in step 4 into the random forest classifier, set the hyperparameters related to the random forest algorithm, and train the cultivated land object recognition model;

[0013] Step 6: Combine the high-resolution image, hyperspectral image and DEM data of the test area into a feature vector of each superpixel ground object according to the features extracted in step 2 and step 3, and then input the feature vector into the cultivated land object recognition model in step 5 for prediction, so as to obtain the category of each superpixel object in the test area and extract the cultivated land in the test area.

[0014] Step 7: Based on the cultivated land in the test area extracted in step 6 and the cultivated land vector data that has been investigated and counted, spatial overlay analysis is used to automatically extract the patches where cultivated land becomes non-cultivated land. First, the patches are deleting the fragmented and narrow patches through the patch area, aspect ratio, and compactness threshold. Then, the patches are filtered to eliminate holes, jagged edges, and sharp corners by using patch boundary simplification and rounding, hole filling, boundary index, and shape index, and candidate cultivated land non-agricultural patches are screened out. Finally, the area threshold is set according to the actual business to obtain the cultivated land non-agricultural patches.

[0015] Furthermore, in step 1, the training area and the test area are divided, and high-resolution images, hyperspectral images, digital elevation model DEM data, and true vector data of cultivated land in the training area are obtained; and high-resolution images, hyperspectral images, digital elevation model DEM data, and investigated and statistical cultivated land vector data in the test area are obtained, including:

[0016] An area with typical mountain cultivated land characteristics is selected to divide it into a training area and a test area, and high-resolution images and hyperspectral remote sensing images of the same area with similar phases are collected in the training area and the test area. The high-resolution images and hyperspectral remote sensing images are registered and geometrically precisely corrected, and DEM data of the corresponding area is collected. The cultivated land true value vector data is the cultivated land range vector extracted from the survey and monitoring results data of the corresponding phase to the high-resolution image, and the cultivated land range vector is the cultivated land true value vector; the cultivated land results vector data that have been investigated and monitored are the cultivated land results vector data that have been investigated and monitored before the current high-resolution image phase.

[0017] Furthermore, in step 2, the high-resolution image is first segmented into superpixels using the SNIC algorithm embedded in the edge detection Canny, and the specific process is as follows:

[0018] Step 2.1: For the high-resolution red, green, and blue bands, denoted as Red, Green, and Blue, respectively, use the Canny operator to extract the edges and obtain the edge detection results of the three bands, denoted as Red Canny 、Green Canny 、Blue Canny ;

[0019] Step 2.2: Perform multiplication and addition operations on the high-resolution red, green, and blue original bands and the edge detection results of the three bands to obtain the red, green, and blue bands with edge information, denoted as Red edge 、Green edge 、Blue edge , the calculation formula is as follows:

[0020]

[0021] Step 2.3, Red edge 、Green edge 、Blue edge The three-band image is converted to CIELAB color space, where l is the brightness channel in CIELAB color space, a and b are the two color channels in CIELAB color space, a is from magenta to green, and b is from yellow to blue;

[0022] Step 2.4: The l, a, and b channels in the CIELAB color space are used as the processing objects of the SNIC algorithm. The number of superpixels output by the SNIC algorithm is set to n according to the resolution of the high-resolution image, the minimum area of ​​cultivated land, and the area of ​​the area to be segmented. The weight for balancing color priority and space priority is c. According to the image resolution and the area size of the object to be interpreted, the values ​​of n and c are adjusted to generate compact and approximately uniform superpixels, so that the boundary of the superpixel is as close to that of the object as possible. The number of superpixels n output by the SNIC algorithm is calculated using the following formula:

[0023]

[0024] Where S is the area of ​​the high-resolution image to be segmented, in square meters; PNmin is the number of pixels set based on the image resolution and the minimum area of ​​cultivated land.

[0025] Furthermore, in step 3, the superpixel ground object obtained in step 2 is used as a unit to extract features within the range of each superpixel ground object, specifically including: taking the superpixel ground object as a unit,

[0026] A. Calculate the spatial domain features of the high-resolution image within each superpixel object using the gray-level co-occurrence matrix GLCM, wherein the spatial domain features include at least five spatial texture features: angular second moment ASM, contrast CON, correlation COR, entropy ENT, and homogeneity HOM, denoted as meanASM, meanCON, meanCOR, meanENT, and meanHOM;

[0027] B. Use Fourier transform to obtain the spectrum, filter the low-frequency information through Butterworth high-pass filter, and then use inverse Fourier transform to obtain the frequency domain high-frequency features of the high-resolution image within each super-pixel object, recorded as meanHighfrequency;

[0028] C. The spectral domain characteristics of the hyperspectral image within each super-pixel object range are obtained by using the first three principal components obtained by principal component analysis, the improved soil adjusted vegetation index MSAVI, and the sensitive response band characteristics of cultivated land, which are recorded as meanPCA1, meanPCA2, meanPCA3, meanMSAVI, meanBandtop1, meanBandtop2, and meanBandtop3;

[0029] D. Utilize ArcGIS software to analyze the raster surface and calculate the terrain features of the DEM data within the range of each superpixel object. The terrain features include at least slope, aspect, hillshade, undulation and surface roughness, which are recorded as meanSlope, meanAspect, meanHillshade, meanUndulation and meanRoughness.

[0030] Furthermore, the gray level co-occurrence matrix GLCM is used to calculate the spatial domain features of the high-resolution image within the range of each super-pixel object. The specific process is as follows:

[0031] First, the high-resolution image multi-band average is converted into a single-band grayscale image, and the grayscale level is compressed to L (L is generally 16, 32, 64). The w*w neighborhood of each pixel is used as a window. w is determined according to the image resolution and the area of ​​the object to be interpreted. The GLCM of four different angles, namely horizontal, vertical and two oblique directions, with a pixel interval is calculated. Four GLCMs are obtained for each neighborhood window. The GLCM is a matrix with L rows and L columns.

[0032] Secondly, for each texture feature, the four values ​​calculated based on the four GLCMs are averaged as the texture feature value of the central pixel of the neighborhood window, and the sliding window traverses all pixels to obtain the texture feature image of this kind, and finally five texture feature images are obtained;

[0033] The formulas for calculating the five texture features based on GLCM are:

[0034]

[0035] Among them, P(i,j) is the value of the i-th row and j-th column in GLCM, μ1 and μ2 are the means, σ1 and σ2 are the variances, and the calculation formula is as follows:

[0036]

[0037] Finally, for each superpixel object, the pixel average value of each texture feature image inside it is calculated, recorded as meanASM, meanCON, meanCOR, meanENT, and meanHOM, to obtain the five spatial texture feature values ​​of the superpixel object.

[0038] Furthermore, the frequency spectrum is obtained by Fourier transform, low-frequency information is filtered by Butterworth high-pass filter, and high-frequency features of the high-resolution image in the frequency domain within each super-pixel object are obtained by inverse Fourier transform. The specific process is as follows:

[0039] First, the high-resolution image multi-band average is converted into a single-band grayscale image f(x, y), and the spectrum F(u, v) is obtained by Fourier transform:

[0040]

[0041] Where x and y are the horizontal and vertical pixel positions of the single-band grayscale image, variables u and v are frequency variables in the spectrum graph, j is an imaginary unit, and M and N are the horizontal and vertical pixel sizes of the single-band grayscale image f(x,y);

[0042] Secondly, for the spectrum graph F(u,v), the low-frequency information is filtered through the Butterworth high-pass filter to obtain G(u,v):

[0043]

[0044] Where D0 is the radius threshold of the dot at the center of the spectrum graph, and D(u,v) is the distance from any position in the spectrum graph to the center;

[0045] Then, the frequency spectrum G(u,v) after Butterworth high-pass filtering is inversely transformed by Fourier transform to obtain the high-frequency information graph f(x,y) H :

[0046]

[0047] The variables u and v are frequency variables in the spectrum graph, j is an imaginary unit, and x and y are high-frequency information graphs f(x,y) H Horizontal and vertical pixel positions;

[0048] Finally, for each superpixel object, calculate its internal high-frequency information map f(x,y) H The average value of is recorded as meanHighfrequency, and the high-frequency feature value in the frequency domain of the super-pixel object is obtained.

[0049] Furthermore, the spectral domain characteristics of the hyperspectral image within each super-pixel object range are obtained by using the first three principal components obtained by principal component analysis, the modified soil adjusted vegetation index MSAVI, and the cultivated land sensitive response band characteristics, specifically including:

[0050] First, perform principal component analysis (PCA) on the hyperspectral image to obtain the first three principal components PCA1, PCA2, and PCA3;

[0051] Secondly, use the red band Red and near infrared band NIR of the hyperspectral image to calculate MSAVI;

[0052]

[0053] Then, pure pixel areas of cultivated land and background were selected on the hyperspectral image, and the pixel difference PD between cultivated land and background in each band was calculated. All bands were arranged in descending order according to PD values, and the top three bands with the highest PD values ​​were sensitive response bands of cultivated land.

[0054] PD=1-Cos(CL histogram ,BG histogram )

[0055] Among them, CL histogram and BG histogram is the pixel grayscale histogram of cultivated land and background in a certain band, and Cos is the cosine distance;

[0056] Finally, for each superpixel ground object, the internal PCA1, PCA2, PCA3, MSAVI, and the average values ​​of the three cultivated land sensitive response bands are calculated, recorded as meanPCA1, meanPCA2, meanPCA3, meanMSAVI, meanBandtop1, meanBandtop2, and meanBandtop3, to obtain the spectral feature values ​​of the superpixel ground object in the spectral domain.

[0057] Furthermore, in the step 4, the features extracted in the steps 2 and 3 are combined into a feature vector Features = [meanASM, meanCON, meanCOR, meanENT, meanHOM, meanHighfrequency, meanPCA1, meanPCA2, meanPCA3, mean MSAVI, meanBandtop1, meanBandtop2, meanBandtop3, meanSlope, meanAspect, meanHillshade, meanUndulation, meanRoughness] of each superpixel object.

[0058] Further, in step 4, determining the label of the superpixel ground object according to the largest proportion of the pixel area of the category label within the superpixel ground object specifically includes:

[0059] The label Label of the superpixel object is determined according to the largest proportion of the category label area within the object. 1 represents the cultivated land category, and 0 represents the background category. The formula is as follows:

[0060]

[0061] where S 耕地 and S 背景 are respectively the pixel areas of cultivated land and background within the object.

[0062] Further, step 6 also includes performing smoothing filtering on the categories of each superpixel object in the test area and extracting the cultivated land in the test area.

[0063] Further, for step 7, the specific process is as follows:

[0064] Step 7.1: Based on the surveyed and monitored cultivated land vector data results, perform spatial erase analysis on the cultivated land vector interpreted in the test area using the random forest model to obtain patches that are cultivated land in the surveyed and monitored cultivated land vector data results but not in the interpreted cultivated land vector, and obtain the preliminary suspected cultivated land non-agriculturalization patch vector shapedata1;

[0065] Step 7.2: For the preliminary suspected cultivated land non-agriculturalization patch vector shapedata1, calculate the area S, aspect ratio B, and compactness C of each patch. Set the thresholds S0, B0, and C0. If S < S0 or B > B0 or C < C0, then delete the fragmented and narrow patches, and then fill the small holes within the patches. Fill the patches with a hole area SK < SK0 to obtain the suspected cultivated land non-agriculturalization patch vector shapedata2;

[0066] The calculation formula for the aspect ratio B of the patch is:

[0067]

[0068] where len and wid are respectively the length and width of the minimum bounding rectangle of the patch;

[0069] The calculation formula for the compactness C of the patch is:

[0070]

[0071] where S and P are respectively the area and perimeter of the patch;

[0072] Step 7.3: For the suspected cultivated land non-agriculturalization patch vector shapedata2, use the Douglas–Peucker algorithm to simplify the patch, eliminate the sawtooth on the patch boundary, and then use the Bezier interpolation method to smooth the patch boundary and eliminate the sharp corners of the patch. Then calculate the boundary index BorderIndex and shape index ShapeIndex of each patch, and set the thresholds BI0 and SI0. If BorderIndex < BI0 and ShapeIndex < SI0, retain the patches that meet the conditions to obtain the candidate cultivated land non-agriculturalization patch vector shapedata3;

[0073] The calculation formulas for the boundary index BorderIndex and the shape index ShapeIndex are respectively:

[0074]

[0075] where S and P are the area and perimeter of the patch respectively, and extent P is the perimeter of the minimum bounding rectangle of the patch;

[0076] Step 7.4: According to the requirement that only the cultivated land reduction patches with a certain area in practice are cultivated land non-agriculturalization patches, for the candidate cultivated land non-agriculturalization patch vector shapedata3, select the threshold ST0 according to business requirements, and obtain the cultivated land non-agriculturalization patch result vector according to S > ST0.

[0077] The remarkable effect of the present invention is:

[0078] In view of the relevant methods of object-level feature extraction and fusion, the present invention first obtains superpixels by using simple non-iterative clustering SNIC that embeds edge detection Canny to segment high-resolution images, and then uses the Gaussian mixture algorithm GMM to cluster the superpixels to obtain ground objects, segment and cluster the high-resolution images, obtain accurate cultivated land objects, and compared with the existing object segmentation methods, the boundaries of the segmented objects are more consistent with the actual ground objects.

[0079] Then, taking the superpixel cultivated land object as the processing unit, taking into account the spatial domain features, frequency domain features, spectral domain features and terrain features, the spatial domain features and high-frequency features of the superpixel ground feature objects are extracted from the high-resolution image, the spectral features of the superpixel ground feature objects are extracted from the hyperspectral image, and the terrain features of the superpixel ground feature objects are extracted from the digital elevation DEM, and the obtained multi-dimensional and multi-scale features are combined into a feature vector; taking the random forest as the classifier, the random forest classification algorithm is used to construct the extraction model of the superpixel cultivated land object, and then accurately interpret the cultivated land, realizing the extraction of cultivated land and the extraction of cultivated land non-agricultural patches by integrating multi-dimensional and multi-scale remote sensing and terrain features, and improving the recognition accuracy of cultivated land in complex mountainous areas. This method fully considers the differences in spectral characteristics between cultivated land and other land types, coordinates hyperspectral images and high-resolution images, and makes full use of the spectral characteristics of hyperspectral images and the spatial structural characteristics of high-resolution images to improve the recognition accuracy of cultivated land. In addition, the cultivated land in the southwestern mountainous areas has strong terrain characteristics, such as sloping cultivated land and trough cultivated land with significant slope, slope direction, and mountain shadow characteristics. Therefore, terrain characteristics will play an important role in promoting the extraction of cultivated land and the monitoring of cultivated land non-agriculturalization in the southwestern mountainous areas. This method uses feature extraction and fusion at the superpixel object level to effectively solve the problem of accurate identification of cultivated land in complex mountainous environments and intelligent monitoring of cultivated land non-agriculturalization.

[0080] Finally, the cultivated land vector data results of the survey and monitoring statistics are used as the benchmark, such as the historical cultivated land results vector data of the third national land survey; the result vector of the random forest model interpreting cultivated land is erased and superimposed, and the map spots where cultivated land becomes non-cultivated land are automatically extracted to realize the non-agriculturalization monitoring of cultivated land in the southwestern mountainous areas, and the cultivated land non-agriculturalization map spots are filtered out by optimizing, screening, and area threshold, so that the data obtained is more accurate. The present invention has significant advantages such as practicality, convenience, and feasibility. BRIEF DESCRIPTION OF THE DRAWINGS

[0081] Figure 1 It is a schematic diagram of the structure of the present invention;

[0082] Figure 2 The high-resolution images, hyperspectral image data and DEM data of the present invention;

[0083] Figure 3 This is a graph showing the results of the method of the present invention for extracting cultivated land in the test area;

[0084] Figure 4 It is a map of the intelligent monitoring results of non-agriculturalization of cultivated land in the test area using the method of the present invention. DETAILED DESCRIPTION

[0085] The specific implementation manner and working principle of the present invention are further described in detail below with reference to the accompanying drawings.

[0086] See also Figure 1 and Figure 4 As shown in the figure, a method for extracting farmland non-agriculturalization by integrating multi-dimensional and multi-scale remote sensing with terrain features is presented. The specific steps are as follows:

[0087] Step 1: Divide the training area and the test area, obtain high-resolution images, hyperspectral images, digital elevation model (DEM) data, and true vector data of cultivated land in the training area, and obtain high-resolution images, hyperspectral images, DEM data, and investigated and statistical cultivated land vector data in the test area, including:

[0088] An area with typical mountain cultivated land characteristics is selected to divide the training area and the test area, and high-resolution images and hyperspectral remote sensing images of the same area with similar phases are collected in the training area and the test area. In this specific embodiment: high-resolution images and hyperspectral remote sensing images with a phase difference of less than one month are used to ensure the consistency of the collected data, and the hyperspectral remote sensing images are registered and geometrically precisely corrected, and the digital elevation model DEM data of the corresponding area is collected. The cultivated land true value vector data is the cultivated land range vector extracted from the survey and monitoring results data of the corresponding phase of the high-resolution image, that is, the cultivated land range vector extracted from the monitoring results data of the same year as the high-resolution image, and the cultivated land range vector is the cultivated land true value vector; the cultivated land results vector data that have been investigated and monitored are the cultivated land results vector data that have been investigated and monitored before the current high-resolution image phase. In this specific embodiment: the cultivated land vector data that have been investigated and monitored use the cultivated land range vector data of the third national land survey results.

[0089] like Figure 2 As shown, this embodiment uses worldview2's 0.5-meter high-resolution image, Zhuhai No. 1 10-meter hyperspectral image and 10-meter resolution DEM as experimental data. It is located in the jurisdiction of a city in the southwestern mountainous area. The imaging time of the 0.5-meter high-resolution image and the 10-meter hyperspectral image are both in September 2022. The production time of the 10-meter DEM data is 2020. The true value vector of cultivated land comes from the 2022 cultivated land monitoring results data of a city in the southwestern mountainous area. The training area and test area cover 420 square kilometers and 200 square kilometers respectively.

[0090] Step 2: firstly, the high-resolution image is segmented into superpixels using the simple non-iterative clustering (SNIC) algorithm embedded with edge detection Canny, and then the superpixel objects are clustered using the Gaussian mixture algorithm (GMM) in combination with the high-resolution image to obtain superpixel ground objects;

[0091] First, the high-resolution image is segmented into superpixels using the SNIC algorithm embedded in the edge detection Canny. The specific process is as follows:

[0092] Step 2.1: For the high-resolution red, green, and blue bands, denoted as Red, Green, and Blue, respectively, use the Canny operator to extract the edges and obtain the edge detection results of the three bands, denoted as Red Canny 、Green Canny 、Blue Canny ;

[0093] Step 2.2: Perform multiplication and addition operations on the high-resolution red, green, and blue original bands and the edge detection results of the three bands to obtain the red, green, and blue bands with edge information, denoted as Red edge 、Green edge 、Blue edge , the calculation formula is as follows:

[0094]

[0095] Step 2.3, Red edge 、Green edge 、Blue edge The three-band image is converted to CIELAB color space, where l is the brightness channel in CIELAB color space, with a value range of 0 (black) to 100 (white), a and b are the two color channels in CIELAB color space, a ranges from magenta to green, and b ranges from yellow to blue;

[0096] Step 2.4, the three channels l, a, and b in the CIELAB color space are used as the processing objects of the SNIC algorithm. According to the resolution of the high-resolution image, the minimum area of ​​cultivated land, and the area of ​​the area to be segmented, the number of superpixels output by the SNIC algorithm is set to n, and the weight for balancing color priority and space priority is c. According to the image resolution and the area size of the object to be interpreted, the values ​​of n and c are adjusted to generate compact and approximately uniform superpixels, so that the boundary of the superpixel is as close to the boundary of the object as possible. If possible, the weight c for balancing color priority and space priority in the SNIC algorithm is 10, so that the superpixel segmentation result is better and its boundary is consistent with the boundary of the object.

[0097] The number of superpixels n output by the SNIC algorithm is as follows:

[0098]

[0099] Where S is the area of ​​the high-resolution image to be segmented, in square meters; PNmin is the number of pixels set according to the resolution of the high-resolution image and the minimum area of ​​cultivated land.

[0100] In this specific embodiment, the area of ​​cultivated land in the southwestern mountainous area ranges from tens to thousands of square meters. To ensure that small cultivated land plots are segmented separately, the minimum area of ​​the cultivated land plot to be segmented is set to 50 square meters, which occupies 50 / (0.5*0.5)=200 pixels on the 0.5-meter image. The number of superpixels n output by the SNIC algorithm segmentation is determined according to the following formula:

[0101]

[0102] Therefore, the value of n in the training area is approximately 2100000, and the value of n in the test area is approximately 1000000.

[0103] In this specific embodiment, the use of the Gaussian mixture algorithm GMM in combination with high-resolution images to cluster superpixel objects specifically includes: using the GMM algorithm in combination with high-resolution images to cluster superpixel objects, with the number of clusters k=2, to obtain superpixel ground objects.

[0104] Step 3, taking the superpixel ground object obtained in step 2 as a unit, extracting features within the range of each superpixel ground object, wherein the features at least include spatial domain features of high-resolution images, frequency domain high-frequency features of high-resolution images, spectral domain features of hyperspectral images, and terrain features of DEM data;

[0105] Specifically include: taking super-pixel ground objects as units,

[0106] A. Calculate the spatial domain features of the high-resolution image within each superpixel object using the gray level co-occurrence matrix GLCM, wherein the spatial domain features include at least five spatial texture features: angular second moment ASM, contrast CON, correlation COR, entropy ENT, and homogeneity HOM;

[0107] In this specific embodiment, the spatial domain features include five texture features: angular second moment (ASM), contrast (CON), correlation (COR), entropy (ENT), and homogeneity (HOM). The specific process is as follows:

[0108] First, convert the high-resolution image into a single-band grayscale image and compress the grayscale to L. L is generally 16, 32, or 64. The w*w neighborhood of each pixel is used as a window. w is determined according to the image resolution and the area of ​​the object to be interpreted. The higher the image resolution and the larger the area of ​​the object, the larger w is required. Calculate the GLCM of four different angles, namely horizontal, vertical, and two oblique directions, with a one-pixel interval. Each neighborhood window obtains four GLCMs, and the GLCM is a matrix of L rows and L columns. In this specific example: compress the grayscale to L=64, use the 7*7 neighborhood of each pixel as a window, calculate the GLCM of four different angles, namely 0°, 45°, 90°, and 135°, with a one-pixel interval. Each neighborhood window obtains four GLCMs, and the GLCM is a matrix of 64 rows and 64 columns, which improves the redundancy of the data.

[0109] Then, for each texture feature, the four values ​​calculated based on the four GLCMs are averaged as the texture feature value of the central pixel of the neighborhood window, and the sliding window traverses all pixels to obtain the texture feature image, and finally five texture feature images can be obtained; the formula for calculating the five texture features based on GLCM is:

[0110]

[0111] Among them, P(i,j) is the value of the i-th row and j-th column in GLCM, μ1 and μ2 are the means, σ1 and σ2 are the variances, and the calculation formula is as follows:

[0112]

[0113] Finally, for each superpixel object, the pixel average value of each texture feature image inside it is calculated, recorded as meanASM, meanCON, meanCOR, meanENT, and meanHOM, to obtain the five spatial texture feature values ​​of the superpixel object.

[0114] B. Use Fourier transform to obtain the spectrum, filter the low-frequency information through Butterworth high-pass filter, and then use inverse Fourier transform to obtain the frequency domain high-frequency features of the high-resolution image within the range of each super-pixel object. The specific process is as follows:

[0115] First, the high-resolution image multi-band average is converted into a single-band grayscale image f(x, y), and the spectrum F(u, v) is obtained by Fourier transform:

[0116]

[0117] Where x and y are the horizontal and vertical pixel positions of the single-band grayscale image, variables u and v are frequency variables in the spectrum graph, j is an imaginary unit, and M and N are the horizontal and vertical pixel sizes of the single-band grayscale image f(x,y);

[0118] Secondly, for the spectrum graph F(u,v), the low-frequency information is filtered through the Butterworth high-pass filter to obtain G(u,v):

[0119]

[0120] Where D0 is the radius threshold of the dot at the center of the spectrum graph. In this embodiment, D0 is set to 0.5, and D(u, v) is the distance from any position in the spectrum graph to the center.

[0121] Then, the frequency spectrum G(u,v) after Butterworth high-pass filtering is inversely transformed by Fourier transform to obtain the high-frequency information graph f(x,y) H :

[0122]

[0123] The variables u and v are frequency variables in the spectrum graph, and x and y are high-frequency information graphs f(x,y) H The horizontal and vertical pixel positions, j is an imaginary unit.

[0124] Finally, for each superpixel object, calculate its internal high-frequency information map f(x,y) H The average value of is recorded as meanHighfrequency, and the high-frequency feature value in the frequency domain of the super-pixel object is obtained.

[0125] C. The spectral domain characteristics of the hyperspectral image within each super-pixel object are obtained by using the first three principal components obtained by principal component analysis, the modified soil adjusted vegetation index MSAVI, and the sensitive response band characteristics of cultivated land;

[0126] The specific process is as follows:

[0127] First, PCA principal component analysis is performed on the hyperspectral image to obtain the first three principal components PCA1, PCA2, and PCA3;

[0128] Secondly, the red band (Red) and near infrared band (NIR) of the hyperspectral image are used to calculate MSAVI; since the wavelength range of the Red band is 620-760 nanometers, and the wavelength range of the NIR band is 780-1100 nanometers, the Zhuhai No. 1 hyperspectral image used in this example has 32 bands, the Red band can be the 15th band (central wavelength is 686 nanometers) of the Zhuhai No. 1 hyperspectral image, and the NIR band can be the 32nd band (central wavelength is 940 nanometers) of the Zhuhai No. 1 hyperspectral image;

[0129]

[0130] Then, pure pixel areas of cultivated land and background are selected on the hyperspectral image, and the pixel difference PD between cultivated land and background in each band is calculated. In this example, the three bands with the highest PD values ​​are bands 8, 5, and 27, which are used as sensitive response bands for cultivated land.

[0131] PD=1-Cos(CL histogram ,BG histogram )

[0132] Among them, CL histogram and BG histogram is the pixel grayscale histogram of cultivated land and background in a certain band, and Cos is the cosine distance.

[0133] Finally, for each superpixel ground object, the average values ​​of its internal PCA1, PCA2, PCA3, MSAVI and the three cultivated land sensitive response bands are calculated, recorded as meanPCA1, meanPCA2, meanPCA3, meanMSAVI, meanBandtop1, meanBandtop2, and meanBandtop3, to obtain the spectral feature values ​​of the superpixel ground object in the spectral domain.

[0134] D. Calculate the terrain features of the DEM data within each superpixel object using ArcGIS software raster surface analysis. The terrain features include at least slope, aspect, hillshade, undulation and surface roughness. The specific process is as follows:

[0135] The raster surface analysis tool of ArcGIS software is used to calculate the terrain characteristics of DEM data within the range of each object, mainly including slope, aspect, hillshade, undulation and surface roughness characteristics. For each superpixel object, the average values ​​of Slope, Aspect, Hillshade, Undulation and Roughness are calculated, recorded as meanSlope, meanAspect, meanHillshade, meanUndulation and meanRoughness, to obtain the terrain characteristic value of the superpixel object.

[0136] In this specific embodiment, the Slope of DEM data is calculated using the Spatial Analyst Tools—>Surface—>Slope tool of ArcGIS software, the Aspect of DEM data is calculated using the Spatial Analyst Tools—>Surface—>Aspect tool of ArcGIS software, and the Hillshade of DEM data is calculated using the Spatial Analyst Tools—>Surface—>Hillshade tool of ArcGIS software; the Spatial Analyst Tools—>Neighborhood—>Focal Statistics tool of ArcGIS software is used to set the Statistic type to MAXIMUM, the size of the neighborhood is 5×5 to obtain the maximum value MAXvalue of the terrain, and then set the Statistic type to MINIMUM, the size of the neighborhood is 5×5 to obtain the minimum value MINvalue of the terrain, and finally the Spatial Analyst Tools—>MapAlgebra—>Raster Calculator tool of ArcGIS software is used to calculate MAXvalue-MINvalue to obtain Undulation; the Spatial Analyst Tools—>Map Algebra—>Raster Calculator tool of ArcGIS software can be used to calculate Roughness according to the following formula:

[0137]

[0138] Step 4: Combine the high-resolution image, high-resolution image, DEM data and true value vector data of the training area into a feature vector of each super-pixel object according to the features extracted in step 2 and step 3, determine the super-pixel object label according to the pixel area with the largest category label in the super-pixel object, and construct super-pixel cultivated land object training sample data;

[0139] The feature vector of the superpixel object is expressed as:

[0140] Features=[meanASM,meanCON,meanCOR,meanENT,meanHOM,meanHighfrequency,meanPCA1,meanPCA2,meanPCA3,meanMSAVI,meanBandtop1,meanBandtop2,meanBandtop3,meanSlope,meanAspect,meanHillshade,meanUndulation,meanRoughness]

[0141] Determining the superpixel ground object label according to the pixel area of ​​the category label in the superpixel ground object having the largest proportion specifically includes:

[0142] The label Label of the superpixel object is determined according to the maximum proportion of the category label area in the object, 1 represents the cultivated land class, and 0 represents the background class. The formula is as follows:

[0143]

[0144] Among them, S 耕地 and S 背景 are the pixel areas of cultivated land and background within the object, respectively.

[0145] Step 5: Input the superpixel cultivated land object training sample data constructed in step 4 into the random forest classifier, set the hyperparameters related to the random forest algorithm, and train the cultivated land object recognition model; in this specific embodiment: input the superpixel ground feature object sample data of the training area into the random forest classifier, set 20 trees in the random forest classifier, set the maximum depth to 10, and iterate the training of the cultivated land object recognition model until convergence.

[0146] Step 6: Combine the high-resolution image, hyperspectral image and DEM data of the test area into a feature vector of each superpixel object according to the features extracted in step 2 and step 3, and then input the feature vector into the cultivated land object recognition model in step 5 for prediction, obtain the category of each superpixel object in the test area, and extract the cultivated land in the test area. Figure 3 As shown in the figure, the method can achieve good results in extracting cultivated land in the test area, accurately distinguishing cultivated land from easily confused features such as gardens, woodlands, and grasslands, and the geometric morphology and boundary shapes of the cultivated land patches can be accurately superimposed with high-resolution images.

[0147] Possibly, the method further includes performing smoothing filtering on the category of each superpixel object in the test area to extract cultivated land in the test area, thereby improving the smoothness of the image.

[0148] Step 7: Based on the cultivated land in the test area extracted in Step 6 and the surveyed and statistically analyzed cultivated land vector data, use spatial overlay analysis to automatically extract the patches where cultivated land has become non-cultivated land. First, delete the fragmented and narrow patches from these patches according to the patch area, aspect ratio, and compactness threshold. Then, use patch boundary simplification and smoothing, hole filling, boundary index, and shape index to eliminate holes, filter out sawteeth and sharp corners, and screen out the candidate cultivated land non-agriculturalization patch vectors. Finally, set an area threshold according to the actual business to obtain the cultivated land non-agriculturalization patches.

[0149] The specific process is as follows:

[0150] Step 7.1: Based on the results of the cultivated land vector data of the third national land survey, perform spatial erase analysis on the cultivated land vector interpreted in the test area using the random forest model. In this example, use the Erase spatial erase tool in ArcGIS software to obtain the patches that are cultivated land in the results of the cultivated land vector data of the third national land survey but not cultivated land in the interpreted cultivated land vector, and obtain the preliminary suspected cultivated land non-agriculturalization patch vector shapedata1.

[0151] Step 7.2: For the preliminary suspected cultivated land non-agriculturalization patch vector shapedata1, calculate the area S, aspect ratio B, compactness C, etc. of each patch, set the thresholds S0, B0, C0, and delete the fragmented and narrow patches (S < S0 or B > B0 or C < C0). Then, fill the small holes inside the patches, and fill those with a hole area SK < SK0 inside each patch to obtain the suspected cultivated land non-agriculturalization patch vector shapedata2.

[0152] Optionally, for the area S, aspect ratio B, and compactness C of the patch, the calculation formula for the aspect ratio B of the patch is:

[0153]

[0154] where len and wid are the length and width of the minimum bounding rectangle of the patch respectively, and the calculation formula for the compactness C of the patch is:

[0155]

[0156] where S and P are the area and perimeter of the patch respectively.

[0157] In this example, considering the characteristics of small and scattered sloping cultivated land plots in the mountainous environment of the southwest, the values of S0, B0, C0, and SK0 are set to 50 square meters, 100, 0.0001, and 20 square meters respectively, and use the vector editing tool in ArcGIS software to fill the small holes inside the patches.

[0158] Step 7.3: For the suspected cultivated land non - agriculturalization polygon vector shapedata2, use the Douglas - Peucker algorithm to simplify the polygon, eliminate the sawteeth on the polygon boundary, and then use the Bezier interpolation method to smooth the polygon boundary and eliminate the sharp corners of the polygon. Then calculate the boundary index BorderIndex and shape index ShapeIndex of each polygon, set the thresholds BI0 and SI0, and retain the polygons that satisfy BorderIndex < BI0 and ShapeIndex < SI0 to obtain the candidate cultivated land non - agriculturalization polygon vector shapedata3;

[0159] The calculation formulas for the boundary index BorderIndex and shape index ShapeIndex are as follows:

[0160]

[0161] where S and P are the area and perimeter of the polygon respectively, and extent P is the perimeter of the minimum bounding rectangle of the polygon. Given the requirements of actual business for cultivated land non - agriculturalization polygons, BI0 and SI0 are set to 5 and 8 respectively.

[0162] In this specific embodiment, the Douglas - Peucker algorithm for polygon simplification uses the Cartography Tools—>Generalization—>Simplify Polygon tool in ArcGIS software, with the tolerance set to 5 meters. The Bezier interpolation method for polygon boundary smoothing uses the Cartography Tools—>Generalization—>Smooth Polygon tool in ArcGIS software.

[0163] Step 7.4: According to the requirement in actual business that only the cultivated land reduction polygons with a certain area are cultivated land non - agriculturalization polygons, for the candidate cultivated land non - agriculturalization polygon vector shapedata3, select the threshold ST0 according to business requirements. The value of ST0 is generally taken as 50, 100, or 200. If the value of ST0 is too small, it is meaningless in practice. Obtain the cultivated land non - agriculturalization polygon result vector according to S > ST0. In this specific example: the threshold ST0 = 200 square meters, and the cultivated land non - agriculturalization polygon result vector is obtained according to S > 200 square meters.

[0164] Such as Figure 4As shown in the figure, this example uses random forest to realize the arable land extraction and arable land non-agricultural monitoring by integrating multi-dimensional and multi-scale features of high-resolution images, hyperspectral images, and DEM data. The trained model can better identify the arable land in the test area. The method of the present invention uses a machine learning model to realize the arable land extraction by integrating multi-dimensional and multi-scale remote sensing and terrain features, and uses superpixel objects as processing units to improve the computing efficiency, automatically obtain arable land non-agricultural patches, and realize the arable land non-agricultural monitoring in the southwestern mountainous areas. Compared with the manual visual judgment of arable land changes in complex mountainous areas, it has significant advantages in practicality, convenience, feasibility, etc.

[0165] The technical solution provided by the present invention is described in detail above. Specific examples are used herein to illustrate the principles and implementation methods of the present invention. The description of the above embodiments is only used to help understand the method of the present invention and its core idea. It should be pointed out that for ordinary technicians in this technical field, without departing from the principle of the present invention, the present invention can also be improved and modified in a number of ways, and these improvements and modifications also fall within the scope of protection of the claims of the present invention.

Claims

1. An intelligent monitoring method for farmland non-agriculturalization integrating multi-dimensional and multi-scale remote sensing and terrain features, characterized in that: The steps include: Step 1: Divide the training area and the test area, obtain high-resolution images, hyperspectral images, digital elevation model DEM data, and true vector data of cultivated land in the training area, and obtain high-resolution images, hyperspectral images, digital elevation model DEM data, and vector data of cultivated land results that have been investigated and monitored in the test area; Step 2: firstly, the high-resolution image is segmented into superpixels using the simple non-iterative clustering (SNIC) algorithm embedded with edge detection Canny, and then the superpixel objects are clustered using the Gaussian mixture algorithm (GMM) in combination with the high-resolution image to obtain superpixel ground objects; Step 3, taking the superpixel ground object obtained in step 2 as a unit, extracting features within the range of each superpixel ground object, wherein the features at least include spatial domain features of high-resolution images, frequency domain high-frequency features of high-resolution images, spectral domain features of hyperspectral images, and terrain features of DEM data; Step 4: Combine the high-resolution image, high-resolution image, DEM data and true value vector data of the training area into a feature vector of each super-pixel object according to the features extracted in step 2 and step 3, determine the super-pixel object label according to the pixel area with the largest category label in the super-pixel object, and construct super-pixel cultivated land object training sample data; Step 5: input the super-pixel cultivated land object training sample data constructed in step 4 into the random forest classifier, set the hyperparameters related to the random forest algorithm, and train the cultivated land object recognition model; Step 6, combining the high-resolution image, hyperspectral image and DEM data of the test area into a feature vector of each super-pixel object according to the features extracted in step 2 and step 3, and then inputting the feature vector into the cultivated land object recognition model in step 5 for prediction, obtaining the category of each super-pixel object in the test area, and extracting the cultivated land in the test area; Step 7: Based on the cultivated land in the test area extracted in step 6 and the cultivated land vector data that has been investigated and monitored, spatial overlay analysis is used to automatically extract the patches where cultivated land becomes non-cultivated land. First, the patches are deleting the fragmented and narrow patches through the patch area, aspect ratio, and compactness threshold. Then, the patches are simplified and rounded, hole filling, boundary index, and shape index are used to filter out holes, jagged edges, and sharp corners of the patches, and the candidate cultivated land non-agricultural patches are screened out. Finally, the area threshold is set according to the actual business to obtain the cultivated land non-agricultural patches.

2. The intelligent monitoring method for non-agriculturalization of cultivated land by integrating multi-dimensional and multi-scale remote sensing with terrain features according to claim 1 is characterized by: The step 1: divide the training area and the test area, obtain high-resolution images, hyperspectral images, digital elevation model DEM data, and true vector data of cultivated land in the training area, and obtain high-resolution images, hyperspectral images, digital elevation model DEM data, and vector data of cultivated land results that have been investigated and monitored in the test area, including: An area with typical mountain cultivated land characteristics is selected to divide it into a training area and a test area, and high-resolution images and hyperspectral remote sensing images of the same area with similar phases are collected in the training area and the test area. The high-resolution images and hyperspectral remote sensing images are registered and geometrically precisely corrected, and DEM data of the corresponding area is collected. The cultivated land true value vector data is the cultivated land range vector extracted from the survey and monitoring results data of the corresponding phase to the high-resolution image, and the cultivated land range vector is the cultivated land true value vector; the cultivated land results vector data that have been investigated and monitored are the cultivated land results vector data that have been investigated and monitored before the current high-resolution image phase.

3. The intelligent monitoring method for non-agriculturalization of cultivated land by integrating multi-dimensional and multi-scale remote sensing with terrain features according to claim 1 is characterized by: In step 2, the high-resolution image is first segmented into superpixels using the SNIC algorithm embedded in the edge detection Canny, and the specific process is as follows: Step 2.1: For the high-resolution red, green, and blue bands, denoted as Red, Green, and Blue, respectively, use the Canny operator to extract the edges and obtain the edge detection results of the three bands, denoted as Red Canny 、Green Canny 、Blue Canny ; Step 2.2: Perform multiplication and addition operations on the high-resolution red, green, and blue original bands and the edge detection results of the three bands to obtain the red, green, and blue bands with edge information, denoted as Red edge 、Green edge 、Blue edge , the calculation formula is as follows: Step 2.3, Red edge 、Green edge 、Blue edge The three-band image is converted to CIELAB color space, where l is the brightness channel in CIELAB color space, a and b are the two color channels in CIELAB color space, a is from magenta to green, and b is from yellow to blue; Step 2.4: The l, a, and b channels in the CIELAB color space are used as the processing objects of the SNIC algorithm. The number of superpixels output by the SNIC algorithm is set to n according to the resolution of the high-resolution image, the minimum area of ​​cultivated land, and the area of ​​the area to be segmented. The weight for balancing color priority and space priority is c. According to the image resolution and the area size of the object to be interpreted, the values ​​of n and c are adjusted to generate compact and approximately uniform superpixels, so that the boundary of the superpixel is as close to that of the object as possible. The number of superpixels n output by the SNIC algorithm is calculated using the following formula: Where S is the area of ​​the high-resolution image to be segmented, in square meters; PNmin is the number of pixels set according to the resolution of the high-resolution image and the minimum area of ​​cultivated land.

4. The intelligent monitoring method for non-agriculturalization of cultivated land by integrating multi-dimensional and multi-scale remote sensing with terrain features according to claim 1 is characterized by: In step 3, the superpixel ground object obtained in step 2 is used as a unit to extract features within the range of each superpixel ground object, specifically including: taking the superpixel ground object as a unit, A. Calculate the spatial domain features of the high-resolution image within each superpixel object using the gray-level co-occurrence matrix GLCM, wherein the spatial domain features include at least five spatial texture features: angular second moment ASM, contrast CON, correlation COR, entropy ENT, and homogeneity HOM, denoted as meanASM, meanCON, meanCOR, meanENT, and meanHOM; B. Use Fourier transform to obtain the spectrum, filter the low-frequency information through Butterworth high-pass filter, and then use inverse Fourier transform to obtain the frequency domain high-frequency features of the high-resolution image within each super-pixel object, recorded as meanHighfrequency; C. The spectral domain characteristics of the hyperspectral image within each super-pixel object range are obtained by using the first three principal components obtained by principal component analysis, the improved soil adjusted vegetation index MSAVI, and the sensitive response band characteristics of cultivated land, which are recorded as meanPCA1, meanPCA2, meanPCA3, meanMSAVI, meanBandtop1, meanBandtop2, and meanBandtop3; D. Utilize ArcGIS software to analyze the raster surface and calculate the terrain features of the DEM data within the range of each superpixel object. The terrain features include at least slope, aspect, hillshade, undulation and surface roughness, which are recorded as meanSlope, meanAspect, meanHillshade, meanUndulation and meanRoughness.

5. The intelligent monitoring method for non-agriculturalization of cultivated land by integrating multi-dimensional and multi-scale remote sensing with terrain features according to claim 4 is characterized by: The gray level co-occurrence matrix GLCM is used to calculate the spatial domain features of the high-resolution image within each super-pixel object. The specific process is as follows: First, the high-resolution image multi-band average is converted into a single-band grayscale image, and the grayscale level is compressed to L. The w*w neighborhood of each pixel is used as a window, where w is determined according to the image resolution and the area of ​​the object to be interpreted. The GLCM of four different angles, namely horizontal, vertical and two oblique directions, with a pixel interval is calculated. Four GLCMs are obtained for each neighborhood window, and the GLCM is a matrix with L rows and L columns. Secondly, extract each texture feature, take the average of the four values ​​calculated based on the four GLCMs as the texture feature value of the central pixel of the neighborhood window, slide the window through all pixels to obtain the texture feature image, and finally obtain five texture feature images; The formulas for calculating five texture features based on the GLCM are: Among them, P(i,j) is the value of the i-th row and j-th column in GLCM, μ1 and μ2 are the means, σ1 and σ2 are the variances, and the calculation formula is as follows: Finally, for each superpixel object, the pixel average of each texture feature image inside it is calculated to obtain five spatial texture feature values ​​of the superpixel object, which are recorded as meanASM, meanCON, meanCOR, meanENT, and meanHOM.

6. The intelligent monitoring method for non-agriculturalization of cultivated land by integrating multi-dimensional and multi-scale remote sensing with terrain features according to claim 4 is characterized by: The frequency spectrum is obtained by Fourier transform, low-frequency information is filtered by Butterworth high-pass filter, and high-frequency features of the high-resolution image in the frequency domain within each super-pixel object are obtained by inverse Fourier transform. The specific process is as follows: First, the high-resolution image multi-band average is converted into a single-band grayscale image f(x, y), and the spectrum F(u, v) is obtained by Fourier transform: Where x and y are the horizontal and vertical pixel positions of the single-band grayscale image, variables u and v are frequency variables in the spectrum graph, j is an imaginary unit, and M and N are the horizontal and vertical pixel sizes of the single-band grayscale image f(x,y); Secondly, for the spectrum graph F(u,v), the low-frequency information is filtered through the Butterworth high-pass filter to obtain G(u,v): Where D0 is the radius threshold of the dot at the center of the spectrum graph, and D(u,v) is the distance from any position in the spectrum graph to the center; Then, the frequency spectrum G(u,v) after Butterworth high-pass filtering is inversely transformed by Fourier transform to obtain the high-frequency information graph f(x,y) H : The variables u and v are frequency variables in the spectrum graph, and x and y are high-frequency information graphs f(x,y) H The horizontal and vertical pixel positions, j is an imaginary unit; Finally, for each superpixel object, calculate its internal high-frequency information map f(x,y) H The average value of the frequency domain high-frequency feature value of the super-pixel object is obtained, which is recorded as meanHighfrequency.

7. The intelligent monitoring method for non-agriculturalization of cultivated land by integrating multi-dimensional and multi-scale remote sensing with terrain features according to claim 4 is characterized by: The spectral domain characteristics of the hyperspectral image within each super-pixel object range are obtained by using the first three principal components obtained by principal component analysis, the improved soil adjusted vegetation index MSAVI, and the sensitive response band characteristics of cultivated land. The specific process is as follows: First, perform principal component analysis (PCA) on the hyperspectral image to obtain the first three principal components, PCA1, PCA2, and PCA3; Second, calculate the modified soil adjusted vegetation index (MSAVI) using the red band (Red) and near-infrared band (NIR) of the hyperspectral image; Then, select pure pixel regions of cultivated land and background on the hyperspectral image, calculate the pixel difference (PD) between cultivated land and background for each band, and arrange all bands in descending order of PD values. The top three bands with the highest PD values are the cultivated land sensitive response bands; PD=1-Cos(CL histogram ,BG histogram ) Among them, CL histogram and BG histogram is the pixel grayscale histogram of cultivated land and background in a certain band, and Cos is the cosine distance; Finally, for each superpixel object, calculate the average values of PCA1, PCA2, PCA3, MSAVI, and the three cultivated land sensitive response bands within it to obtain the spectral domain spectral feature values of the superpixel object, denoted as meanPCA1, meanPCA2, meanPCA3, meanMSAVI, meanBandtop1, meanBandtop2, and meanBandtop3.

8. The intelligent monitoring method for non-agriculturalization of cultivated land by integrating multi-dimensional and multi-scale remote sensing with terrain features according to claim 4 is characterized by: In step 4, according to the features extracted in step 2 and step 3, combine them to form the feature vector of each superpixel object: Features = [meanASM, meanCON, meanCOR, meanENT, meanHOM, meanHighfrequency, meanPCA1, meanPCA2, meanPCA3, meanMSAVI, meanBandtop1, meanBandtop2, meanBandtop3, meanSlope, meanAspect, meanHillshade, meanUndulation, meanRoughness].

9. The intelligent monitoring method for non-agriculturalization of cultivated land by integrating multi-dimensional and multi-scale remote sensing with terrain features according to claim 1, characterized in that: In step 6, it also includes smoothing the class of each superpixel object in the test area and extracting the cultivated land in the test area.

10. The intelligent monitoring method for non-agriculturalization of cultivated land by integrating multi-dimensional and multi-scale remote sensing with terrain features according to claim 1, characterized in that: Step 7 specifically includes: Step 7.1: Based on the surveyed and monitored cultivated land vector data results, perform spatial erase analysis on the cultivated land result vector interpreted in the test area using the random forest model to obtain the patches that are cultivated land in the surveyed and monitored cultivated land vector data results but not in the interpreted result vector, and obtain the preliminary suspected cultivated land non-agriculturalization patch vector shapedata1; Step 7.2: For the preliminary suspected cultivated land non-agriculturalization patch vector shapedata1, calculate the area S, aspect ratio B, and compactness C of each patch, set thresholds S0, B0, and C0. If S < S0 or B > B0 or C < C0, then delete the fragmented and long and narrow patches, and then fill the small holes within the patches. Fill the patches with hole area SK < SK0 to obtain the suspected cultivated land non-agriculturalization patch vector shapedata2; The calculation formula for the aspect ratio B of the patch is: where len and wid are the length and width of the minimum bounding rectangle of the patch, respectively; The calculation formula for the compactness C of the patch is: where S and P are the area and perimeter of the patch, respectively; Step 7.3: For the suspected cultivated land non-agriculturalization patch vector shapedata2, use the Douglas–Peucker algorithm to simplify the patch, eliminate the serrations on the patch boundary, and then use the Bezier interpolation method to smooth the patch boundary and eliminate the sharp corners of the patch. Then calculate the boundary index BorderIndex and shape index ShapeIndex of each patch, and set the thresholds BI0 and SI0. If BorderIndex < BI0 and ShapeIndex < SI0, then retain the patches that meet the conditions to obtain the candidate cultivated land non-agriculturalization patch vector shapedata3; The calculation formulas for the boundary index BorderIndex and the shape index ShapeIndex are respectively: Among them, S and P are the area and perimeter of the patch respectively, and extent P is the perimeter of the minimum circumscribed rectangle of the patch; Step 7.4: According to the requirement that only the cultivated land reduction patches with a certain area in practice are cultivated land non-agriculturalization patches, for the candidate cultivated land non-agriculturalization patch vector shapedata3, select the threshold ST0 according to business requirements, and obtain the cultivated land non-agriculturalization patch result vector according to S > ST0.

Images