A method for automatic identification of high-standard farmland plot boundaries based on satellite images

Abstract

The present invention discloses a method for automatically identifying the boundaries of high-standard farmland plots based on satellite images, which includes the following steps: A. preparing training data: downloading Google satellite images with a resolution of 1 meter of a representative area, using professional GIS software to draw the boundary vectors of high-standard farmland plots, making label data, rasterizing the drawn vectors and converting them into edge images. The present invention introduces the attention mechanism of transformer, uses high weights to focus on important features, and low weights to ignore irrelevant features, which can better capture the contextual information of boundary features. In complex terrain areas such as mountains and hills, it can retain local fine-grained information to the greatest extent, reduce feature loss, and use the LSD plot boundary closure detection method in the post-processing process to reduce the connection between different plots, further improve the accuracy of plot boundaries, and avoid artificial correction of the threshold of the binarization threshold segmentation step, simplifying the difficulty of use.

Description

Technical Field

[0001] The present invention relates to the technical field of remote sensing image processing, and in particular to a method for automatically identifying boundaries of high-standard farmland plots based on satellite images. Background Art

[0002] Current methods for obtaining high-standard farmland boundaries include manual drawing, traditional remote sensing image segmentation methods based on spectral, texture, shape and other features, and edge detection models based on deep learning:

[0003] First, manually drawn plot boundaries have the highest accuracy, but they require a lot of manpower and material resources, have a long cycle, and are slow to update, making them difficult to meet the requirements of large-scale high-standard farmland;

[0004] Second, traditional remote sensing image segmentation methods based on spectral, texture, shape and other features are subject to the mutual exclusivity between image spatial resolution and spectral resolution. High spatial resolution images lack spectral information, while high spectral resolution images have low spatial resolution. In addition, the spatial homogeneity of ground objects, such as "same spectrum, different objects" and "same object, different spectrum" phenomena, further increase the difficulty of accurate remote sensing image segmentation.

[0005] 3. With the rapid development of artificial intelligence technology, deep learning-based edge detection models (HED, RCF, etc.) have begun to be used to extract plot boundaries. This type of method first uses the edge detection model to calculate the plot boundary probability map, then performs threshold segmentation and binarization, and finally performs post-processing operations such as skeleton extraction and raster vectorization. Compared with traditional remote sensing image segmentation methods, although the accuracy and efficiency have been improved, when calculating the boundary probability map, the model performs downsampling to expand the receptive field, resulting in the loss of local fine-grained information. In addition, in the process of binarization threshold segmentation, there is a phenomenon of original field interconnection (different plots are connected as one). Due to these two reasons, the edge detection model based on deep learning performs poorly in complex terrains such as hills and mountains. Summary of the Invention

[0006] The purpose of the present invention is to provide a method for automatically identifying the boundaries of high-standard farmland plots based on satellite images to solve the problems raised in the above background technology.

[0007] To achieve the above objectives, the present invention provides the following technical solution: a method for automatically identifying the boundaries of high-standard farmland plots based on satellite images, the method comprising the following steps:

[0008] A. Prepare training data: Download 1-meter-resolution Google satellite images of representative areas, use professional GIS software to draw high-standard farmland boundary vectors, create label data, rasterize the drawn vectors and convert them into edge images, and create training data based on the original images and edge images.

[0009] B. Global context modeling: The input image is divided into 4*4 coarse-grained blocks. The global encoder calculates global attention based on the coarse-grained blocks. Then, a bidirectional multi-stage aggregate decoder is used to enhance the information flow in the encoder through two paths, top-down and bottom-up, to generate high-resolution features.

[0010] C. Local refinement: Use a non-overlapping sliding window to sample the entire image and divide it into multiple fine-grained block sequences. Then perform local encoding on each sequence in turn to capture short-range local information. All local information is integrated and input into a local bidirectional multi-stage aggregation decoder to further obtain pixel-level feature maps, i.e., fine-grained features.

[0011] D. Feature fusion: Using global context information as prior knowledge, adjust local fine-grained features to generate fused features that contain global context and fine-grained local details. The fused features are input into the local decision module to generate a boundary probability map.

[0012] E. Binarize the boundary probability map and extract the skeleton: Binarize the boundary probability map output in step D using the threshold segmentation method. Then use the skeleton extraction algorithm described in the paper "A fast parallel algorithm for thinning digital patterns" to simplify the binarized land contours into thin lines.

[0013] F. LSD plot closure detection and restoration: The extracted skeleton is scaled to reduce the jagged effect in the image. The gradients of the boundary pixels are then calculated and sorted by gradient magnitude. Pixels with small gradient values ​​correspond to smooth or slowly changing areas. Pixels with large gradient magnitudes are used as seeds. Points with deviations less than the tolerance are traversed within the eight-neighborhood. These points are added to the area and updated. A rectangular approximation is performed on the updated area. The alignment point density of the approximated rectangle is then calculated. The rectangle is cropped based on the alignment point density to obtain a completely closed plot.

[0014] G. Raster vectorization, hole removal, and smoothing to generate the final result: The results of the LSD plot closure detection output are vectorized to generate plot boundary lines, and linear features are converted into surface features. The area threshold method is used to remove small polygons and holes, and the smoothing function is used to remove boundary jagged edges to improve the visual effect and obtain the final high-standard farmland plot boundary data.

[0015] Preferably, step A is specifically as follows: ① Select a typical and representative high-standard farmland distribution area, obtain Google satellite images with a resolution of 1 meter, and use professional GIS software to outline the boundaries of the high-standard farmland in these areas to form a vector file of surface elements; ② Rasterize the outlined plot boundary vector file based on the Google image, set the boundary pixels to 1, and set the non-boundary pixels to 0 as label data; ③ Sliding crop the satellite image data and label data with a window size of 480*480 and a step size of 360 to generate a sample set, and divide the sample set into a training set, a validation set, and a test set in a ratio of 8:1:1.

[0016] Preferably, the step B is specifically as follows: ① normalize the training data to satisfy the normal distribution and normalize the data range to 0-1; ② input the normalized training data into the global encoder; ③ the global encoder uses the Vit method to divide the entire image into 4*4 coarse blocks, and passes through 24 Transformer blocks. It is the output of the first block to the 24th block. The output is divided into four groups, that is, the features are divided into four stages. Only the output of the last Transformer block in each group represents the output feature of this stage. Input into the global decoder for decoding, the formula is: Where, represents the output of the transformer block, G e Represents the global encoder; ④ Reshape the 2D features into 3D image features through 1*1 convolution, and then accumulate them in a top-down and bottom-up manner, multi-level aggregated features, and the 8 feature maps generated by the decoder are cascaded to obtain a coarse-grained edge map. The formula is: Where, f g represents the pixel-level global feature, G d Represents the global decoder.

[0017] Preferably, the step C is specifically as follows: ① Divide the picture into X 1 , X 2 , X 3 , X 4 Four blocks, each of which is further subdivided into 4*4 fine-grained blocks and fed into the shared Transformer local encoder to calculate the attention weight; ②X 1 Input into the transformer local encoder to generate X 1 of X 2 Input into the transformer local encoder to generate X 2 of X3 Input into the transformer local encoder to generate X 2 of X 4 Input into the transformer local encoder to generate X 4 of ③ Take the output features of the last fine-grained block in each group and convert X 1 , X 2 , X 3 , X 4 Four groups of output Connect in sequence, then reshape into a three-dimensional array, feed it into the local decoder, and output fine-grained features. The formula is: Where, f r Represents local features, R d Represents a local decoder.

[0018] Preferably, in step D, a spatial feature transformation block and two 3*3 convolutional layers are formed, followed by batch normalization and relu activation, and finally the fused features are fed into the local decision head to predict the boundary probability map. The formula is: E r =R h (FFM(f g , f r )), where R h It is a local decision head consisting of a 1*1 convolutional layer and a sigmoid operation.

[0019] Preferably, the step F is specifically as follows: ① reducing the skeleton data obtained in step E to 80% of the original image, then convolving it with a Gaussian kernel to smooth and eliminate the aliasing effect, and then downsampling; ② calculating the gradient of the image using a 2*2 convolution kernel, and the gradient calculation formula is:

[0020]

[0021]

[0022] Where f(x,y) is the pixel value, the amplitude G of the image gradient and the angle θ of the row and column lines are calculated as follows: ③ The pixel area with smaller gradient amplitude corresponds to the flat area in the image. Smaller gradient values ​​often appear in smooth areas. Pixels with gradient amplitude less than the threshold will be filtered out and not participate in the calculation of the support line area; ④ In the filtered pixel list, select unused pixels as seed points, and search for angles that meet the eight neighborhoods of unused points to form support line areas; ⑤ Perform rectangular approximation on each support line area and construct a specific rectangle containing all points in the area. Specifically, the entire area is first regarded as an entity. The value of each pixel in the area is the gradient size, so that the entire entity has a center of mass. The center of mass is used as the center point of the rectangle, and then the orientation angle of the rectangle is determined. The angle is set to the angle of the eigenvector related to the minimum eigenvalue of the matrix. After determining the center and orientation of the rectangle, the minimum rectangle containing all points is the approximated rectangle; ⑥ Calculate the alignment point density of the approximation rectangle, crop the approximation rectangle according to the alignment point density, and obtain a complete closed plot. The formula is: Where k is the number of alignment points, length is the length of the approximation rectangle, and width is the width of the approximation rectangle.

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

[0024] The present invention introduces the attention mechanism of transformer, which replaces the downsampling method of convolutional neural network. It uses high weight to focus on important features and low weight to ignore irrelevant features, which can better capture the contextual information of boundary features. In complex terrain areas such as mountains and hills, it retains local fine-grained information to the greatest extent, reduces feature loss, and improves the accuracy of broken and irregular plot boundaries. In the post-processing process, the LSD plot boundary closure detection method is used. Based on the gradient size and change direction of the existing boundary points, the plot boundary breaks and missing caused by threshold segmentation are repaired and closed to form a complete plot, reducing the connection between different plots, further improving the accuracy of the plot boundary, and avoiding the manual correction of the threshold of the binarization threshold segmentation step, simplifying the difficulty of use, and running automatically throughout the process, without the need for manual interaction, threshold adjustment and other operations, thereby improving production efficiency. BRIEF DESCRIPTION OF THE DRAWINGS

[0025] Figure 1 It is a flow chart of the present invention. DETAILED DESCRIPTION

[0026] The following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts are within the scope of protection of the present invention.

[0027] See also Figure 1 A method for automatically identifying boundaries of high-standard farmland plots based on satellite images comprises the following steps:

[0028] A. Prepare training data: Download 1-meter-resolution Google satellite images of representative areas, use professional GIS software to draw high-standard farmland boundary vectors, create label data, rasterize the drawn vectors and convert them into edge images, and create training data based on the original images and edge images.

[0029] B. Global Context Modeling: This enables coarse-grained feature extraction of image data. The input image is segmented into 4x4 coarse-grained blocks. The global encoder calculates global attention based on the coarse-grained blocks. A bidirectional multi-stage aggregate decoder then enhances the information flow in the encoder via both top-down and bottom-up paths, generating high-resolution features.

[0030] C. Local Refinement: This method extracts fine-grained features from images. It samples the entire image using a non-overlapping sliding window, dividing it into multiple fine-grained block sequences. Local encoding is then performed on each sequence to capture short-range local information. All local information is then integrated and input into a local bidirectional multi-stage aggregation decoder to obtain pixel-level feature maps, i.e., fine-grained features.

[0031] D. Feature fusion: Using global context information as prior knowledge, adjust local fine-grained features to generate fused features that contain global context and fine-grained local details. The fused features are input into the local decision module to generate a boundary probability map.

[0032] E. Binarize the boundary probability map and extract the skeleton: Binarize the boundary probability map output in step D using threshold segmentation (the threshold setting does not need to be adaptive for each image; the subsequent LSD plot closure detection will repair broken and missing boundaries). Then, use the skeleton extraction algorithm described in the paper "A fast parallel algorithm for thinning digital patterns" to simplify the binarized plot outlines into thin lines.

[0033] F. LSD plot closure detection and restoration: The extracted skeleton is scaled to reduce the jagged effect in the image. The gradient of the boundary pixels is then calculated and sorted by gradient magnitude. Pixels with small gradient values ​​correspond to smooth or slowly changing areas. These areas are filtered by setting a threshold so that they do not participate in subsequent region growing. Pixels with large gradient magnitudes are used as seeds. Points with deviations less than the tolerance are traversed within the eight-neighborhood. These points are added to the region and updated. A rectangular approximation is performed on the updated area. The alignment point density of the approximated rectangle is then calculated. The rectangle is cropped based on the alignment point density to obtain a completely closed plot.

[0034] G. Raster vectorization, hole removal, and smoothing to generate the final result: The results of the LSD plot closure detection output are vectorized to generate plot boundary lines, and linear features are converted into surface features. The area threshold method is used to remove small polygons and holes, and the smoothing function is used to remove boundary jagged edges to improve the visual effect and obtain the final high-standard farmland plot boundary data.

[0035] Step A is as follows: ① Select representative high-standard farmland distribution areas, obtain 1-meter resolution Google satellite images, and use professional GIS software to outline the boundaries of the high-standard farmland in these areas to form a vector file of surface features; ② Rasterize the outlined plot boundary vector file based on the Google image, set the boundary pixels to 1, and the non-boundary pixels to 0, as label data; ③ Sliding crop the satellite image data and label data with a window size of 480*480 and a step size of 360 to generate a sample set, and divide the sample set into a training set, a validation set, and a test set in a ratio of 8:1:1.

[0036] Step B is as follows: ① Standardize the training data to make it normal distribution and normalize the data range to 0-1; ② Input the standardized training data into the global encoder; ③ The global encoder uses the Vit method to divide the entire image into 4*4 coarse blocks, and passes through 24 Transformer blocks. It is the output of the first block to the 24th block. The output is divided into four groups, that is, the features are divided into four stages. Only the output of the last Transformer block in each group represents the output feature of this stage. Input into the global decoder for decoding, the formula is: Where, represents the output of the transformer block, G eRepresents the global encoder; ④ Reshape the 2D features into 3D image features through 1*1 convolution, and then accumulate them in a top-down and bottom-up manner to aggregate the features at multiple levels. The 8 feature maps generated by the decoder (already upsampled to the original image size) are cascaded to obtain a coarse-grained edge map. The formula is: Where, f g represents pixel-level global features (coarse-grained global context features), G d Represents the global decoder.

[0037] Step C is as follows: ① Divide the image into X 1 , X 2 , X 3 , X 4 Four blocks, each of which is further subdivided into 4*4 fine-grained blocks and fed into the shared Transformer local encoder to calculate the attention weight; ②X 1 Input into the transformer local encoder to generate X 1 of X 2 Input into the transformer local encoder to generate X 2 of X 3 Input into the transformer local encoder to generate X 3 of X 4 Input into the transformer local encoder to generate X 4 of ③ Take the output features of the last fine-grained block in each group and convert X 1 , X 2 , X 3 , X 4 Four groups of output Connect them in sequence, then reshape them into a three-dimensional array and feed them into the local decoder (the decoder structure is the same as the first stage, except that the convolution kernel is changed to 1*1), and output fine-grained features. The formula is: Where, f r Represents local features, R d Represents the local decoder. Unlike the global decoder, the 3*3 convolution in the local encoder is replaced by a 1*1 convolution layer to avoid artificial edges caused by padding.

[0038] In step D, it consists of a spatial feature transformation block and two 3*3 convolutional layers, followed by batch normalization and relu activation. Finally, the fused features are fed into the local decision head to predict the boundary probability map. The formula is: E r =R h (FFN(f g , fr )), where R h It is a local decision head consisting of a 1*1 convolutional layer and a sigmoid operation.

[0039] Step F is as follows: ① Reduce the skeleton data obtained in step E to 80% of the original image, then convolve it with a Gaussian kernel to smooth out the aliasing effect, and then downsample it to prevent other interference noise; ② Use a 2*2 convolution kernel to calculate the gradient of the image, ensuring the calculation speed while also ensuring the relative independence of the neighborhood direction distribution. The gradient calculation formula is:

[0040]

[0041]

[0042] Where f(x,y) is the pixel value, the amplitude G of the image gradient and the angle θ of the row and column lines are calculated as follows: ③ The pixel area with smaller gradient amplitude corresponds to the flat area in the image. Smaller gradient values ​​often appear in smooth areas. Pixels with gradient amplitude less than the threshold will be filtered out and not participate in the calculation of the support line area; ④ In the filtered pixel list, select unused pixels as seed points, and search for angles that meet the eight neighborhoods of unused points to form support line areas; ⑤ Perform rectangular approximation on each support line area and construct a specific rectangle containing all points in the area. Specifically, the entire area is first regarded as an entity. The value of each pixel in the area is the gradient size, so that the entire entity has a center of mass. The center of mass is used as the center point of the rectangle, and then the orientation angle of the rectangle is determined. The angle is set to the angle of the eigenvector related to the minimum eigenvalue of the matrix. After determining the center and orientation of the rectangle, the minimum rectangle containing all points is the approximated rectangle; ⑥ Calculate the alignment point density of the approximation rectangle, crop the approximation rectangle according to the alignment point density, and obtain a complete closed plot. The formula is: Where k is the number of alignment points, length is the length of the approximation rectangle, and width is the width of the approximation rectangle.

[0043] While embodiments of the present invention have been shown and described, it will be appreciated by those skilled in the art that various changes, modifications, substitutions, and variations may be made to these embodiments without departing from the principles and spirit of the invention, and that the scope of the invention is defined by the appended claims and their equivalents.

Claims

1. A method for automatically identifying boundaries of high-standard farmland plots based on satellite images, characterized by: The method comprises the following steps: A. Prepare training data: Download 1-meter-resolution Google satellite images of representative areas, use professional GIS software to draw high-standard farmland boundary vectors, create label data, rasterize the drawn vectors and convert them into edge images, and create training data based on the original images and edge images. B. Global context modeling: The input image is divided into 4*4 coarse-grained blocks. The global encoder calculates global attention based on the coarse-grained blocks. Then, a bidirectional multi-stage aggregate decoder is used to enhance the information flow in the encoder through two paths, top-down and bottom-up, to generate high-resolution features. C. Local refinement: Use a non-overlapping sliding window to sample the entire image and divide it into multiple fine-grained block sequences. Then perform local encoding on each sequence in turn to capture short-range local information. All local information is integrated and input into a local bidirectional multi-stage aggregation decoder to further obtain pixel-level feature maps, i.e., fine-grained features. D. Feature fusion: Using global context information as prior knowledge, adjust local fine-grained features to generate fused features that contain global context and fine-grained local details. The fused features are input into the local decision module to generate a boundary probability map. E. Binarization of boundary probability map and skeleton extraction: The boundary probability map output in step D is binarized using the threshold segmentation method, and then the skeleton extraction algorithm is used to simplify the binarized land contour into thin lines; F. LSD plot closure detection and restoration: The extracted skeleton is scaled to reduce the jagged effect in the image. The gradients of the boundary pixels are then calculated and sorted by gradient magnitude. Pixels with small gradient values ​​correspond to smooth or slowly changing areas. Pixels with large gradient magnitudes are used as seeds. Points with deviations less than the tolerance are traversed within the eight-neighborhood. These points are added to the area and updated. A rectangular approximation is performed on the updated area. The alignment point density of the approximated rectangle is then calculated. The rectangle is cropped based on the alignment point density to obtain a completely closed plot. G. Raster vectorization, hole removal, and smoothing to generate the final result: The results of the LSD plot closure detection output are vectorized to generate plot boundary lines, and linear features are converted into surface features. The area threshold method is used to remove small polygons and holes, and the smoothing function is used to remove boundary jagged edges to improve the visual effect and obtain the final high-standard farmland plot boundary data.

2. The method for automatically identifying boundaries of high-standard farmland plots based on satellite images according to claim 1, characterized in that: The step A is specifically as follows: ① select a typical high-standard farmland distribution area, obtain Google satellite images with a resolution of 1 meter, and use professional GIS software to outline the boundaries of the high-standard farmland in these areas to form a vector file of surface elements; ② rasterize the outlined plot boundary vector file based on the Google image, set the boundary pixels to 1, and set the non-boundary pixels to 0 as label data; ③ perform sliding cropping on the satellite image data and the label data with a window size of 480*480 and a step size of 360 to generate a sample set, and divide the sample set into a training set, a validation set, and a test set in a ratio of 8:1:

1.

3. The method for automatically identifying high-standard farmland plot boundaries based on satellite images according to claim 1, characterized in that: The specific steps of step B are as follows: ① normalize the training data to satisfy the normal distribution and normalize the data range to 0-1; ② input the normalized training data into the global encoder; ③ the global encoder uses the Vit method to divide the entire image into 4*4 coarse blocks, and passes through 24 Transformer blocks. It is the output of the first block to the 24th block. The output is divided into four groups, that is, the features are divided into four stages. Only the output of the last Transformer block in each group represents the output feature of this stage. Input into the global decoder for decoding, the formula is: , where represents the output of the transformer block, Represents the global encoder; ④ Reshape the 2D features into 3D image features through 1*1 convolution, and then accumulate them in a top-down and bottom-up manner, multi-level aggregated features, and the 8 feature maps generated by the decoder are cascaded to obtain a coarse-grained edge map. The formula is: , where represents pixel-level global features, Represents the global decoder.

4. The method for automatically identifying high-standard farmland plot boundaries based on satellite images according to claim 1, characterized in that: The specific step C is: ① Divide the picture into Four blocks, each of which is further subdivided into 4*4 fine-grained blocks and fed into the shared Transformer local encoder to calculate the attention weight;② Input into the transformer local encoder to generate of , Input into the transformer local encoder to generate of , Input into the transformer local encoder to generate of , Input into the transformer local encoder to generate of ; ③ Take the output features of the last fine-grained block in each group and Four groups of output Connect in sequence, then reshape into a three-dimensional array, feed it into the local decoder, and output fine-grained features. The formula is: , where Represents local features, Represents a local decoder.

5. The method for automatically identifying boundaries of high-standard farmland plots based on satellite images according to claim 1, characterized in that: In step D, it consists of a spatial feature transformation block and two 3*3 convolutional layers, followed by batch normalization and relu activation. Finally, the fused features are fed into the local decision head to predict the boundary probability map. The formula is: , where R h It is a local decision head consisting of a 1*1 convolutional layer and a sigmoid operation.

6. The method for automatically identifying boundaries of high-standard farmland plots based on satellite images according to claim 1, characterized in that: The specific steps of step F are as follows: ① reducing the skeleton data obtained in step E to 80% of the original image, then convolving it with a Gaussian kernel to smooth and eliminate the aliasing effect, and then downsampling it; ② calculating the gradient of the image using a 2*2 convolution kernel. The gradient calculation formula is: , , where f(x,y) is the pixel value, the amplitude G of the image gradient and the angle θ of the row and column lines are calculated as follows: , ;③ The pixel area with smaller gradient amplitude corresponds to the flat area in the image. Smaller gradient values ​​often appear in smooth areas. Pixels with gradient amplitude less than the threshold will be filtered out and will not participate in the calculation of the support line area;④ In the filtered pixel list, select unused pixels as seed points, and search for angles that meet the eight neighborhoods of unused points to form support line areas;⑤ Perform rectangular approximation on each support line area and construct a specific rectangle containing all points in the area. Specifically, the entire area is first regarded as an entity. The value of each pixel in the area is the gradient size, so that the entire entity has a center of mass. The center of mass is used as the center point of the rectangle, and then the orientation angle of the rectangle is determined. The angle is set to the angle of the eigenvector related to the minimum eigenvalue of the matrix. After determining the center and orientation of the rectangle, the minimum rectangle containing all points is the approximated rectangle;⑥ Calculate the alignment point density of the approximation rectangle, crop the approximation rectangle according to the alignment point density, and obtain a complete closed plot. The formula is: , where k is the number of alignment points, length is the length of the approximation rectangle, and width is the width of the approximation rectangle.

Images