An Unmanned Aerial Vehicle Remote Sensing Extraction Method for the Boundaries of Patchy Cultivated Land in Flat Terrain Areas

Abstract

The present invention proposes a method for remotely sensing and extracting the boundaries of patchy cultivated land in flat terrain areas by using an unmanned aerial vehicle (UAV), belonging to the field of remote sensing technology. The method includes the following steps: Step 1: Process the UAV images to obtain the digital surface model and orthophoto of the cultivated land; Step 2: Use the digital surface model to obtain the topographic undulation image of the cultivated land, perform threshold segmentation and binarization on it to obtain the edge area and non-edge area of the field block, and then successively apply area threshold filtering and inversion operation to obtain the first initial cultivated land area; Step 3: Calculate the vegetation index and perform binarization based on the orthophoto to obtain the second initial cultivated land area; Step 4: Combine the first initial cultivated land area and the second initial cultivated land area to obtain the candidate cultivated land area, and obtain the binary map of the candidate cultivated land area; Step 5: Optimize the patch edges of the candidate cultivated land area map obtained in Step 4 to obtain the extraction result of the patchy cultivated land boundary. The steps of the present invention are concise and clear, and have a high degree of automation, which is conducive to the rapid extraction of the boundaries of patchy cultivated land in flat terrain areas.

Description

Technical Field

[0001] The present invention belongs to the technical field of remote sensing, and particularly relates to a method for remotely sensing and extracting the boundaries of patchy cultivated land in flat terrain areas by using an unmanned aerial vehicle (UAV). Background Art

[0002] Cultivated land is a natural resource on which humans depend for survival. The "trinity" protection of the quantity, quality, and ecology of cultivated land is the current cultivated land protection policy in China, and the investigation and monitoring of the quantity of cultivated land play a fundamental role in it. With the advancement of the construction of high-standard basic farmland, patchy cultivated land is a common type of farming unit. The automatic and precise extraction of its boundaries is of great significance for cultivated land investigation, cadastral mapping, and efficient field management. Patchy cultivated land in flat terrain areas is characterized by small plots, a height difference between the field surface and its edge, and a coexistence of regular and irregular boundary shapes. It is mosaicked with ground features such as rivers, ditches, roads, and villages to form a typical agricultural landscape, which poses great challenges to its investigation and mapping. Therefore, it is necessary to construct an automatic and precise extraction method for patchy cultivated land to achieve efficient mapping and dynamic monitoring of such farming units.

[0003] At present, the mapping methods for patchy cultivated land mainly include traditional instrument field measurement, satellite remote sensing identification, and visual interpretation of UAV orthophotos. Traditional field measurement methods have high reliability and excellent accuracy, but they require high costs and long time-consuming, and are not suitable for large-scale cultivated land mapping needs. Commonly used medium- and low-resolution optical remote sensing images are helpful for the rapid identification of large-scale cultivated land plots, such as Sentinel-2, Landsat, and MODIS. However, due to the limitation of spatial resolution (10 - 250 m), it is difficult to extract individual plot patches of fragmented cultivated land. The spatial resolution of high-resolution optical images can reach the sub-meter level, such as China's Gaofen-2 satellite and the United States' WorldView. They can reveal the details of cultivated land to a certain extent, but are troubled by "same object, different spectra; different objects, same spectra" and it is difficult to distinguish the edges of farming units. Moreover, due to the limitation of the revisit cycle, the reliability of image acquisition in cloudy, rainy, and foggy areas cannot be guaranteed. In recent years, UAV remote sensing technology has been widely used in mapping due to its high precision, cost-effectiveness, and flexibility. The mapping resolution can reach the centimeter level, which provides great convenience for visual interpretation of cultivated land patches, but it relies on manual experience and has low mapping efficiency. In addition, UAV image processing simultaneously obtains orthophotos representing vegetation cover and digital surface models reflecting field height differences, providing a good opportunity for automatic mapping of cultivated land by mining the image and elevation information of patchy cultivated land. Therefore, there is an urgent need for an automatic and precise extraction method that can cooperate with the orthophoto and elevation features of UAV remote sensing to carry out patchy cultivated land such farming units. Summary of the Invention

[0004] In view of this, the present invention provides a method for remotely sensing and extracting the boundaries of patchy cultivated land in flat terrain areas by using an unmanned aerial vehicle (UAV), aiming to combine digital image processing methods with cartographic principles to quickly extract cultivated land areas without manual labor, achieve a relatively high degree of automation in extracting the boundaries of patchy cultivated land, and facilitate the identification of cultivated land boundaries and the investigation of cultivated land quantity.

[0005] The technical solution adopted by the present invention is as follows:

[0006] A method for remotely sensing and extracting the boundaries of patchy cultivated land in flat terrain areas by using an unmanned aerial vehicle (UAV), characterized by comprising:

[0007] Step 1: Perform three-dimensional reconstruction processing on the UAV images of cultivated land to respectively obtain the digital surface model and orthophoto of the cultivated land;

[0008] Step 2: Calculate the terrain undulation degree image of the cultivated land by using the digital surface model obtained in Step 1, perform threshold segmentation on the terrain undulation degree in the terrain undulation degree image, binaryize the terrain undulation degree image into two types of objects, namely the candidate field edge area and the non-edge area, and then sequentially apply area threshold filtering and value inversion operations to the candidate field edge area to obtain the first initial cultivated land area;

[0009] The specific steps of Step 2 include the following steps:

[0010] Step 2.1: Based on the digital surface model, mask and extract the images of the target cultivated land area around the main roads and branch roads, calculate the terrain undulation degree f of the masked and extracted digital surface model as shown in Equation (1), and obtain the terrain undulation degree image based on the terrain undulation degree f:

[0011] f = h max -h min (1)

[0012] In the formula, each pixel value of the terrain undulation degree f is the difference between the maximum value and the minimum value of the elevation in the rectangular sliding window, h max is the maximum value of the elevation in the rectangular sliding window, h min is the minimum value of the elevation in the rectangular sliding window;

[0013] Step 2.2: Perform threshold segmentation on the terrain undulation degree f by using the following formula to binaryize the terrain undulation degree image into two types of objects, namely the candidate field edge area and the non-edge area:

[0014]

[0015] where, F1 represents the filtered image of the candidate field edge area, f(x, y) is the pixel value at the pixel (x, y), and T1 is the adopted segmentation threshold;

[0016] Step 2.3: Calculate the area A1 of each patch in the candidate field edge area image processed in Step 2.2, and then set an area threshold to filter out small patches from the patches to complete the purification of the candidate field edge area obtained in Step 2.2, as shown in the following formula:

[0017]

[0018] In the formula, F2 represents the binary image of the filtered field edge area, and T2 is the filtering area threshold adopted;

[0019] Step 2.4: Perform an inversion operation on the binary image of the field edge area processed in Step 2.3 to obtain the first initial cultivated area, as shown in the following formula:

[0020]

[0021] Among them, F3 is the first initial cultivated area, and F2(x, y) is the binary value of the original image at the position (x, y).

[0022] Step 3: Use the orthophoto image obtained in Step 1 to calculate the vegetation index, and binarize it to obtain the second initial cultivated area;

[0023] Step 3.1: Based on the orthophoto image obtained in Step 1 and the vegetation spectral response characteristics and the UAV visible light image bands, calculate the visible light band difference vegetation index VDVI, as shown in the following formula:

[0024]

[0025] In the formula, V1 is the result of the visible light band difference vegetation index obtained from the orthophoto image, and its value range is [-1, 1], DN red 、DN green 、DN blue respectively represent the pixel values of the red, green, and blue 3 bands;

[0026] Step 3.2: According to the image histogram of the visible light band difference vegetation index VDVI obtained in Step 3.1, determine the binarization adaptive threshold to maximize the variance between vegetation and non-vegetation, and obtain the binarized image of the cultivated area by vegetation index threshold segmentation. The regional binarized image is the second initial cultivated area, as shown in Equation (6):

[0027]

[0028] Among them, V2 represents the second initial cultivated area obtained by binarization, V1(x, y) is the pixel value at the pixel (x, y), and T3 is the segmentation threshold adopted.

[0029] Step 4: Combine the first initial cultivated area and the second initial cultivated area to obtain a candidate cultivated area, and perform hole filling and area threshold filtering on the candidate cultivated area to obtain a binary image of the candidate cultivated area;

[0030] Step 4.1: Combine the first initial cultivated area and the second initial cultivated area to obtain a candidate cultivated area, and obtain the pixels where both the first initial cultivated area and the second initial cultivated area are cultivated land, as shown in the following formula:

[0031] R1 = F3 × V2 (7)

[0032] where F3 is the first initial cultivated area, V2 is the second initial cultivated area, and R1 is the result of extracting the cultivated area by combining F3 and V2;

[0033] Step 4.2: Based on R1 obtained in Step 4.1, perform an image closing operation to obtain a preliminarily optimized cultivated area, as shown in the following formula:

[0034]

[0035] where g1 represents a rectangular structuring element, R1 is the result of extracting the cultivated area by combining terrain information and visible light index, and R2 is the cultivated area preliminarily optimized by the image closing operation;

[0036] Step 4.3: Perform hole filling on the obtained preliminarily optimized cultivated area R2 to obtain an image R3, and then perform an image opening operation on R3 to obtain a cultivated area optimized again, as shown in the following formula:

[0037]

[0038] where g2 represents a circular structuring element, and R4 is the cultivated area optimized again;

[0039] Step 4.4: Calculate the area A2 of each patch in the cultivated area R4 optimized again, and complete the purification of the cultivated area R4 optimized again through area threshold filtering to obtain a binary image of the candidate cultivated area:

[0040]

[0041] where R5 represents the binary image of the candidate cultivated area, and T3 is the filtering area threshold used.

[0042] Step 5: Optimize the edges of the patches in the binary image of the candidate cultivated area obtained in Step 4 to obtain the final result of extracting the patchy cultivated land boundary.

[0043] Step 5.1: Perform median filtering on the obtained binary image of the candidate cultivated area; the window of the median filtering is the pixel size corresponding to the average width of the field ridges in the area;

[0044] Step 5.2: According to the area threshold set by the convex hull function, the binary image of the candidate cultivated land area is divided into regular-shaped fields and irregular-shaped fields; the area threshold set by the convex hull function is 105% of the pixel area of each field.

[0045] Step 5.3: Optimize the regular-shaped fields by setting the area threshold using the convex hull algorithm, and optimize the irregular-shaped fields by setting the contraction parameter for the coordinates of the field edge points using the concave hull contraction algorithm. The parameter value of the contraction algorithm is 0.7. Then, use the area threshold filtering to obtain the final extraction result of the patchy cultivated land boundary. The threshold for area filtering is selected as 95% of the pixel area of the smallest cultivated land patch in the study area:

[0046]

[0047] Among them, R5 represents the binary image of the candidate cultivated land area, R6 represents the extraction result of the patchy cultivated land boundary area, and T3 is the filtering area threshold adopted.

[0048] To sum up, due to the adoption of the above technical solutions, the beneficial effects of the present invention are as follows:

[0049] 1. The present invention calculates the terrain undulation degree using the digital surface model obtained by light and small unmanned aerial vehicle photogrammetry, realizes the extraction of the edge area of the fields in the cultivated land area and the grid of the initial cultivated land area, calculates the vegetation index using the obtained orthophoto image, realizes the extraction of the grid of the initial cultivated land area, and obtains the candidate cultivated land area by taking the intersection pixels of the two. Then, the digital image processing technology and the principles of cartography are used to realize the extraction of the final cultivated land area.

[0050] 2. The present invention is simple to operate and has high cartography efficiency, realizing the extraction of the patchy cultivated land boundary by unmanned aerial vehicle remote sensing in the flat terrain area supported by unmanned aerial vehicle photogrammetry products and subsequent fine cadastral management.

[0051] 3. The present invention has strong applicability. It is not only applicable to the specific terrain undulation degree and visible light band difference vegetation index in the present invention, but other indicators that can characterize the terrain elevation difference and crop coverage status can also be commonly used in this method. Description of the Drawings

[0052] The present invention will be described by way of examples with reference to the drawings, wherein:

[0053] Figure 1 In (a) and (b), they are respectively the comparison diagrams before and after the threshold segmentation of the terrain undulation degree of the study area obtained in the embodiment of the present invention.

[0054] Figure 2 In (a) and (b), they are respectively the comparison diagrams before and after calculating the vegetation index of the orthophoto image of the study area obtained in the embodiment of the present invention.

[0055] Figure 3 In (a) and (b) respectively are the optimized results of the regular-shaped field plots and irregular-shaped field plots of the study area obtained in the embodiments of the present invention after combining topographic information and visible light index.

[0056] Figure 4 In (a) and (b) respectively are the comparison diagrams of the irregular-shaped field plots before and after optimization of the study area obtained in the embodiments of the present invention after combining topographic information and visible light index.

[0057] Figure 5 Is the comparison diagram of the patchy cultivated land boundary extraction result obtained in the embodiment of the present invention and the reference value.

[0058] Figure 6 Is the schematic flow structure diagram of this embodiment. Detailed implementation manners

[0059] To make the objectives, technical solutions and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are part of the embodiments of the present invention, rather than all of the embodiments. Usually, the components of the embodiments of the present invention described and shown in the accompanying drawings here can be arranged and designed in various different configurations.

[0060] Therefore, the following detailed description of the embodiments of the present invention provided in the accompanying drawings is not intended to limit the scope of the claimed present invention, but merely represents the selected embodiments of the present invention. Based on the embodiments in the present invention, all other embodiments obtained by those of ordinary skill in the art without creative efforts belong to the scope of protection of the present invention.

[0061] It should be noted that, without conflict, the embodiments in the present invention and the features in the embodiments can be combined with each other.

[0062] It should be noted that: similar reference numerals and letters denote similar items in the following drawings. Therefore, once an item is defined in one drawing, it does not need to be further defined and explained in subsequent drawings.

[0063] In the present invention, unless otherwise clearly specified and defined, the first feature being "on" or "under" the second feature may include the direct contact between the first and second features, or may include the situation where the first and second features are not in direct contact but in contact through additional features therebetween. Moreover, the first feature being "above", "over" and "on top of" the second feature includes that the first feature is directly above and obliquely above the second feature, or merely indicates that the horizontal height of the first feature is higher than that of the second feature. The first feature being "under", "beneath" and "underneath" the second feature includes that the first feature is directly below and obliquely below the second feature, or merely indicates that the horizontal height of the first feature is lower than that of the second feature.

[0064] It should be noted that, without conflict, the embodiments in the present invention and the features in the embodiments may be combined with each other.

[0065] Embodiment 1

[0066] As Figure 6 shown, a method for remotely sensing and extracting the boundaries of patchy cultivated land in a flat terrain area by using an unmanned aerial vehicle (UAV) is disclosed in an embodiment of the present invention, which is characterized by comprising:

[0067] Step 1: Perform three-dimensional reconstruction processing on the UAV images of cultivated land to respectively obtain the digital surface model and the orthophoto of the cultivated land;

[0068] Step 2: Calculate the terrain undulation image of the cultivated land by using the digital surface model obtained in Step 1, perform threshold segmentation on the terrain undulation in the terrain undulation image, binarize the terrain undulation image into two types of objects, namely the candidate field edge area and the non-edge area, and then sequentially apply area threshold filtering and value inversion operations to the candidate field edge area to obtain the first initial cultivated land area;

[0069] The specific steps of Step 2 include the following steps:

[0070] Step 2.1: Based on the digital surface model, mask and extract the images of the target cultivated land area around the main roads and branch roads, calculate the terrain undulation f of the masked and extracted digital surface model as shown in Equation (1), and obtain the terrain undulation image based on the terrain undulation f:

[0071] f = h max -h min (1)

[0072] In the formula, each pixel value of the terrain undulation f is the difference between the maximum value and the minimum value of the elevation in a rectangular sliding window, h max is the maximum value of the elevation in the rectangular sliding window, and h min is the minimum value of the elevation in the rectangular sliding window;

[0073] Step 2.2: Perform threshold segmentation on the terrain undulation degree f using the following formula to binarize the terrain undulation degree image into two types of objects: the candidate farmland edge area and the non-edge area:

[0074]

[0075] where F1 represents the filtered candidate farmland edge area image, f(x, y) is the pixel value at pixel (x, y), and T1 is the segmentation threshold adopted;

[0076] Step 2.3: Calculate the area A1 of each patch in the candidate farmland edge area image processed in Step 2.2, and then set an area threshold to filter out small patches for the candidate farmland edge area obtained in Step 2.2 to complete purification, as shown in the following formula:

[0077]

[0078] In the formula, F2 represents the binary image of the filtered farmland edge area, and T2 is the filtering area threshold adopted;

[0079] Step 2.4: Perform an inversion operation on the binary image of the farmland edge area processed in Step 2.3 to obtain the first initial cultivated area, as shown in the following formula:

[0080]

[0081] where F3 is the first initial cultivated area, and F2(x, y) is the binary value of the original image at position (x, y).

[0082] Step 3: Calculate the vegetation index using the orthophoto image obtained in Step 1 to obtain the second initial cultivated area;

[0083] Step 3.1: Based on the vegetation spectral response characteristics and the UAV visible light image bands of the orthophoto image obtained in Step 1, calculate the visible light band difference vegetation index VDVI, as shown in the following formula:

[0084]

[0085] In the formula, V1 is the result of the visible light band difference vegetation index obtained from the orthophoto image, with a value range of [-1, 1], and DN red 、DN green 、DN blue represent the pixel values of the red, green, and blue 3 bands respectively;

[0086] Step 3.2: According to the image histogram of the visible light band difference vegetation index VDVI obtained in Step 3.1, determine the binary adaptive threshold to maximize the variance between vegetation and non-vegetation, and obtain the binary image of the cultivated land area segmented by the vegetation index threshold. The regional binary image is the second initial cultivated land area, as shown in Equation (6):

[0087]

[0088] Among them, V2 represents the second initial cultivated land area obtained by binarization, V1(x, y) is the pixel value at pixel (x, y), and T3 is the segmentation threshold adopted.

[0089] Step 4: Combine the first initial cultivated land area and the second initial cultivated land area to obtain a candidate cultivated land area, and perform hole filling and area threshold filtering on the candidate cultivated land area to obtain a binary map of the candidate cultivated land area;

[0090] Step 4.1: Combine the first initial cultivated land area and the second initial cultivated land area to obtain a candidate cultivated land area, and obtain the pixels where both the first initial cultivated land area and the second initial cultivated land area are cultivated land, as shown in the following formula:

[0091] R1 = F3 × V2 (7)

[0092] Among them, F3 is the first initial cultivated land area, V2 is the second initial cultivated land area, and R1 is the result of extracting the cultivated land area by combining F3 and V2;

[0093] Step 4.2: Based on R1 obtained in Step 4.1, use image closing operation to obtain a preliminary optimized cultivated land area, as shown in the following formula:

[0094]

[0095] Among them, g1 represents a rectangular structuring element, R1 is the result of extracting the cultivated land area by combining terrain information and visible light index, and R2 is the cultivated land area preliminarily optimized by image closing operation;

[0096] Step 4.3: Perform hole filling on the obtained preliminarily optimized cultivated land area R2 to obtain image R3, and then perform image opening operation on R3 to obtain a cultivated land area optimized again, as shown in the following formula:

[0097]

[0098] Among them, g2 represents a circular structuring element, and R4 is the cultivated land area optimized again;

[0099] Step 4.4: Calculate the area A2 of each patch in the cultivated land area R4 optimized again, and complete the purification of the cultivated land area R4 optimized again through area threshold filtering to obtain a binary map of the candidate cultivated land area:

[0100]

[0101] Among them, R5 represents the binary map of the candidate cultivated area, and T3 is the filtering area threshold adopted.

[0102] Step 5: Optimize the edge of the patches in the binary map of the candidate cultivated area obtained in Step 4 to obtain the final extraction result of the patch-shaped cultivated land boundary.

[0103] Step 5.1: Perform median filtering on the obtained binary map of the candidate cultivated area; the window of the median filtering is the pixel size corresponding to the average width of the field ridges in the area;

[0104] Step 5.2: Divide the binary map of the candidate cultivated area into regularly shaped field blocks and irregularly shaped field blocks according to the area threshold set by the convex hull function; the area threshold set by the convex hull function is 105% of the pixel area of each field block.

[0105] Step 5.3: Optimize the regularly shaped field blocks by setting the area threshold using the convex hull algorithm, and optimize the edge point coordinates of the irregularly shaped field blocks by setting the contraction parameter using the concave hull contraction algorithm. The value of the contraction algorithm parameter is 0.7, and then the final extraction result of the patch-shaped cultivated land boundary is obtained by area threshold filtering. The threshold of the area filtering is selected as 95% of the pixel area of the smallest cultivated land patch in the study area:

[0106]

[0107] Among them, R5 represents the binary map of the candidate cultivated area, R6 represents the extraction result of the patch-shaped cultivated land boundary area, and T3 is the filtering area threshold adopted.

[0108] Embodiment 2

[0109] This embodiment is a specific implementation manner of Embodiment 1. This embodiment is located about 50 km west of Chengdu, Sichuan Province, in the Dujiangyan Irrigation Area in the hinterland of the Chengdu Plain, with a flat terrain and a surface elevation of about 550 m. The area has a typical agricultural planting structure of two crops a year, mainly planting wheat, rice, rapeseed, etc.; it belongs to the subtropical humid monsoon climate, with distinct seasons, an average annual temperature of 16.1 °C, and an average annual rainfall of about 1011 mm. After land consolidation, the tillage conditions in this embodiment are good and convenient for mechanized agricultural production. It is mainly patch-shaped cultivated land, with a length of about 50 - 80 m, a width of 35 - 60 m, a width of the field ridge line of about 0.9 m, and a relative height difference from the field surface of about 0.6 m.

[0110] A method for remotely sensing the patch-shaped cultivated land boundary in a flat terrain area by using an unmanned aerial vehicle in this embodiment includes:

[0111] Select UAV images of cultivated land at multiple time phases during the early and mature growth stages of rice. Based on obtaining the digital surface model and orthophoto image through 3D reconstruction, calculate the surface undulation using the digital surface model. Adopt threshold segmentation to binarize it into two types of objects: the edge area of the field block and the non-edge area of the field block; then implement area threshold filtering to obtain a cleaner binary image of the edge area of the field block; then invert the obtained binary image value of the edge area of the field block to obtain the initial cultivated area obtained using terrain information; then calculate the vegetation index using the orthophoto image and perform threshold binarization to obtain the initial cultivated area obtained using the visible light index; then through logical AND operation, obtain the cultivated area common to both the visible light index and terrain information; finally, divide the cultivated area into regular-shaped field blocks and irregular-shaped field blocks according to the convex hull algorithm to optimize the area threshold, and use the convex hull algorithm and concave hull contraction algorithm to perform edge optimization respectively. The specific implementation process is as follows:

[0112] 1) Select UAV images of cultivated land in the study area for rice in August and December 2023 respectively. Add the images and ground control points in the ContextCapture software, and perform 3D reconstruction to obtain the Digital Surface Model (DSM) and Digital Orthophoto Map (DOM). The spatial resolution of both is 2.7 cm; then mask and extract the cultivated land image of the target area according to the study area boundary. Then, use the above formula (1) to calculate the terrain undulation image f with a pixel window of 33×33 corresponding to the average width of the ridge in this example (0.9 m) as the sliding analysis size for the digital surface model.

[0113] 2) Perform threshold segmentation on the terrain undulation image in December 2023 calculated in step 1) using the above formula (2). The threshold is the absolute value of the difference between half of the mean and standard deviation of the terrain undulation, so as to extract the binary image F1 of the candidate edge area of the cultivated land field block.

[0114] 3) Calculate the area of each patch for the candidate edge area of the field block obtained in step 2), and then perform area threshold filtering using the above formula (3) to purify the candidate edge area of the field block obtained in step 2), and obtain the binary image F2 of the edge area of the field block.

[0115] 4) Perform an inversion operation on the filtered binary image of the edge area of the field block obtained in step 3) using the above formula (4), so as to obtain the binary image F3 of the initial cultivated area extracted using terrain information based on the digital surface model.

[0116] 5) Using the orthoimage obtained in step 1) in August 2023 without the near-infrared band, based on the vegetation spectral response characteristics and the visible-light image bands of the unmanned aerial vehicle, and using the above formula (5) to calculate the visible-band difference vegetation index (VDVI), the vegetation index V1 of the visible-light image can be obtained;

[0117] 6) According to the image histogram of the visible-band difference vegetation index obtained in step 5), determine the binary adaptive threshold to maximize the variance between vegetation and non-vegetation, and use the above formula (6) for threshold segmentation to separate the foreground and background pixels of the image as much as possible, so as to obtain the initial cultivated land extraction area V2:

[0118] 7) According to the initial cultivated land extraction area F3 obtained from the digital surface model in step 4) and the initial cultivated land extraction area V2 obtained from the orthoimage in step 6), use the above formula (7) to obtain the result R1 of extracting the cultivated land area by combining terrain information and visible-light index. Based on the obtained R1, perform a closing operation (8) on the image to obtain the preliminary optimized cultivated land area R2. Fill the holes in the obtained preliminary optimized cultivated land area to obtain the image R3, and then perform an opening operation (9) on the image to obtain the re-optimized cultivated land area R4. Calculate the area of each patch in the cultivated land area obtained by the opening operation of the image, and finally obtain the binary image R5 of the candidate cultivated land area through formula (10);

[0119] 8) According to whether the shape of the field is regular or irregular, divide the binary image of the candidate cultivated land area obtained in step 7) into regular-shaped fields and irregular-shaped fields according to the area threshold set by the convex hull function. For regular-shaped fields, use the convex hull algorithm to set the area threshold for optimization. For irregular-shaped fields, use the concave hull contraction algorithm to set the contraction parameter for the coordinates of the field edge points for optimization. Perform a logical OR operation on the optimized results of the two, and then perform area threshold filtering through formula (11) to obtain the final patchy cultivated land boundary area R6.

[0120] Verification of the effectiveness of the present invention:

[0121] Taking the orthoimage with a spatial resolution of 2.7 cm as the base map, use the visual interpretation method to extract the cultivated land area, which is used as the verification data. The results show that the extraction accuracy of cultivated land fields based on the number of fields is 93.8%, the extraction accuracy of cultivated land area based on pixels is 96.8%, the corresponding overall accuracy is 95.0%, and the Kappa coefficient is 0.89, which proves the effectiveness of the present invention.

[0122] Although embodiments of the present invention have been shown and described, those of ordinary skill in the art will appreciate that various changes, modifications, substitutions, and variations can be made to these embodiments without departing from the principles and spirit of the present invention. The scope of the present invention is defined by the appended claims and their equivalents.

[0123] The circuits, electronic components, and modules involved are all prior art and can be fully implemented by those skilled in the art without further elaboration. The content protected by the present invention does not involve improvements to software and methods either.

[0124] In this specification, the various embodiments are described in a progressive manner. Each embodiment focuses on the differences from other embodiments, and the same or similar parts among the embodiments can be referred to each other.

[0125] The above description of the disclosed embodiments enables those skilled in the art to implement or use the present invention. Various modifications to these embodiments will be obvious to those skilled in the art, and the general principles defined herein can be implemented in other embodiments without departing from the spirit or scope of the present invention. Therefore, the present invention will not be limited to the embodiments shown herein, but rather to the broadest scope consistent with the principles and novel features disclosed herein.

Claims

1. A method for extracting patchy cultivated land boundaries from drone remote sensing in flat terrain areas, characterized in that: include: Step 1: Perform 3D reconstruction on the UAV images of cultivated land to obtain the digital surface model and orthophoto of the cultivated land respectively; Step 2: Using the digital surface model obtained in step 1, the terrain relief image of the cultivated land is calculated, and the terrain relief in the terrain relief image is segmented by threshold value, and the terrain relief image is binarized into two types of objects: the edge area of ​​the candidate field and the non-edge area. Then, the edge area of ​​the candidate field is sequentially subjected to area threshold filtering and value inversion operations to obtain the first initial cultivated land area; Step 3: Use the orthophoto obtained in step 1 to calculate the vegetation index to obtain the second initial cultivated land area; The step 3 specifically comprises the following steps: Step 3.1: Calculate the visible light band difference vegetation index VDVI based on the orthophoto obtained in step 1 according to the vegetation spectral response characteristics and the visible light image band of the drone, as shown in the following formula: Where V1 is the difference vegetation index result of the visible light band obtained by orthophoto, and its value range is [-1,1], and DN red DN green DN blue Respectively represent the pixel values ​​of the red, green and blue bands; Step 3.2: According to the image histogram of the visible light band difference vegetation index VDVI obtained in step 3.1, determine the binarization adaptive threshold so as to maximize the variance between vegetation and non-vegetation, and obtain the binary image of the cultivated land area segmented by the vegetation index threshold. The regional binary image is the second initial cultivated land area, as shown in formula (6): Among them, V2 represents the second initial cultivated land area obtained by binarization, V1(x, y) is the pixel value at pixel (x, y), and T3 is the segmentation threshold used; Step 4: Combining the first initial cultivated land area and the second initial cultivated land area to obtain a candidate cultivated land area, and performing hole filling and area threshold filtering on the candidate cultivated land area to obtain a binary map of the candidate cultivated land area; The step 4 specifically comprises the following steps: Step 4.1: Combine the first initial cultivated land area and the second initial cultivated land area to obtain candidate cultivated land areas, and obtain pixels where both the first initial cultivated land area and the second initial cultivated land area are cultivated land, as shown in the following formula: R1=F3×V2(7) Among them, F3 is the first initial cultivated land area, V2 is the second initial cultivated land area, and R1 is the result of extracting cultivated land area by combining F3 and V2; Step 4.2: Based on R1 obtained in step 4.1, the image closing operation is used to obtain the preliminary optimized cultivated land area, as shown in the following formula: Among them, g1 represents the rectangular structural element, R1 is the result of extracting the cultivated land area by combining terrain information and visible light index, and R2 is the cultivated land area initially optimized by image closing operation; Step 4.3: Fill the holes of the initially optimized cultivated land area R2 to obtain image R3, and then perform image opening operation on R3 to obtain the optimized cultivated land area again, as shown in the following formula: Among them, g2 represents the circular structural element, and R4 is the re-optimization of the cultivated land area; Step 4.4: Calculate the area A2 of each patch in the re-optimized cultivated land area R4, and purify the re-optimized cultivated land area R4 through area threshold filtering to obtain the candidate cultivated land area binary map: Among them, R5 represents the candidate cultivated land area binary map, T3 is the adopted filtering area threshold; Step 5: Optimize the patch edge of the candidate cultivated land area binary map obtained in step 4 to obtain the final patchy cultivated land boundary extraction result.

2. The method for extracting patchy cultivated land boundaries by UAV remote sensing in flat terrain areas according to claim 1 is characterized in that: The step 2 specifically includes the following steps: Step 2.1: Based on the digital surface model, perform mask extraction around the main roads and branches of the target cultivated land area. The terrain relief f is calculated from the mask-extracted digital surface model, as shown in formula (1). The terrain relief image is obtained based on the terrain relief f: f=h max -h min (1) In the formula, each pixel value of the terrain relief f is the difference between the maximum and minimum values ​​of the rectangular sliding window elevation, h max is the maximum elevation in the rectangular sliding window, h min is the minimum elevation value in the rectangular sliding window; Step 2.2: Use the following formula to perform threshold segmentation on the terrain relief f and binarize the terrain relief image into two types of objects: candidate field edge area and non-edge area: Among them, F1 represents the filtered candidate field edge area image, f(x,y) is the pixel value at pixel (x,y), and T1 is the adopted segmentation threshold; Step 2.3: Calculate the area A1 of each spot in the candidate field edge area image processed in step 2.2, and then set the area threshold to filter the spots and remove small spots, so as to purify the candidate field edge area obtained in step 2.2, as shown in the following formula: In the formula, F2 represents the binary image of the field edge area after filtering, and T2 is the filtering area threshold adopted; Step 2.4: Apply the value inversion operation to the binary image of the field edge area processed in step 2.3 to obtain the first initial cultivated land area, as shown in the following formula: Among them, F3 is the first initial cultivated land area, and F2(x, y) is the binary value of the original image at the position (x, y).

3. The method for extracting patchy cultivated land boundaries by UAV remote sensing in flat terrain areas according to claim 2 is characterized in that: The step 5 specifically comprises the following steps: Step 5.1: Perform median filtering on the obtained binary image of the candidate cultivated land area; Step 5.2: According to the area threshold set by the convex hull function, the candidate cultivated land area binary map is divided into regular-shaped plots and irregular-shaped plots; Step 5.3: For regular-shaped fields, the convex hull algorithm is used to set the area threshold for optimization. For irregular-shaped fields, the concave hull shrinkage algorithm is used to set the shrinkage parameters for the field edge point coordinates for optimization. Then, the area threshold filtering is used to obtain the final patchy cultivated land boundary extraction result: Among them, R5 represents the binary map of candidate cultivated land areas, R6 represents the extraction result of patchy cultivated land boundary areas, and T3 is the filtering area threshold used.

4. The method for extracting patchy cultivated land boundaries by UAV remote sensing in flat terrain areas according to claim 3 is characterized in that: In step 5.1, the window of the median filter is a pixel size corresponding to the average width of the ridges in the area.

5. The method for extracting patchy cultivated land boundaries by UAV remote sensing in flat terrain areas according to claim 3 is characterized in that: In step 5.2, the area threshold set by the convex hull function is 105% of the pixel area of ​​each field.

6. The method for extracting patchy cultivated land boundaries by UAV remote sensing in flat terrain areas according to claim 3 is characterized in that: In step 5.3, the shrinkage algorithm parameter is set to 0.

7.

7. The method for extracting patchy cultivated land boundaries by UAV remote sensing in flat terrain areas according to claim 3 is characterized in that: In step 5.3, the threshold of area filtering is selected as 95% of the pixel area of ​​the smallest cultivated patch in the study area.

Images