Unmanned aerial vehicle remote sensing extraction method for plaque cultivated land boundary in flat terrain area

Abstract

The invention provides an unmanned aerial vehicle remote sensing extraction method for a plaque cultivated land boundary in a terrain flat area, which belongs to the technical field of remote sensing and comprises the following steps of: 1, processing an unmanned aerial vehicle image to obtain a digital surface model and an orthoimage of cultivated land; 2, obtaining a topographic relief degree image of the cultivated land by using a digital surface model, carrying out threshold segmentation binarization on the topographic relief degree image into a field parcel marginal region and a non-marginal region, and carrying out area threshold filtering and value negation operation in sequence to obtain a first initial cultivated land region; step 3, calculating a vegetation index based on the orthoimage and carrying out binaryzation to obtain a second initial cultivated land area; 4, obtaining a candidate cultivated land area by combining the first initial cultivated land area and the second initial cultivated land area, and obtaining a binary image of the candidate cultivated land area; and step 5, optimizing the pattern spot edge of the candidate cultivated land area obtained in the step 4 to obtain a patch-shaped cultivated land boundary extraction result. The method is simple and clear in steps, high in automation degree and beneficial to rapid extraction of plaque cultivated land boundaries in flat terrain areas.

Description

Technical Field

[0001] The invention belongs to the field of remote sensing technology, and particularly relates to a method for extracting patchy cultivated land boundaries using unmanned aerial vehicle (UAV) remote sensing. Background Art

[0002] Cultivated land is a natural resource on which human beings depend for survival. The "three-in-one" protection of cultivated land quantity, quality and ecology is the cultivated land protection policy that my country currently adheres to, and the survey and monitoring of cultivated land quantity play a fundamental role in this. With the advancement of high-standard basic farmland construction, patchy cultivated land is a common type of farming unit. The automatic and fine extraction of its boundaries is of great significance to cultivated land survey and cadastral mapping, efficient field management, etc. Patchy cultivated land in flat terrain areas is characterized by small plots, height differences between the field surface and its edge, and the coexistence of regular and irregular boundary shapes. It is inlaid with land features such as rivers, ditches, roads, and villages to form a typical agricultural landscape, which brings great challenges to its survey and mapping. Therefore, it is necessary to construct an automatic and fine extraction method for patchy cultivated land to achieve efficient mapping and dynamic monitoring of such farming units.

[0003] At present, the methods for mapping patchy cultivated land are mainly based on traditional instrument measurement, satellite remote sensing identification, and visual interpretation of drone orthophotos. Traditional field measurement methods are highly reliable and accurate, but they are costly and time-consuming, and are not suitable for large-scale cultivated land mapping. 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, but due to the limitation of spatial resolution (10 to 250 m), it is difficult to extract individual plots of fine cultivated land. The spatial resolution of high-resolution optical images can reach sub-meter level, such as China's Gaofen-2 satellite and the United States WorldView, which can reveal the details of cultivated land to a certain extent, but they are troubled by "different spectra for the same object and the same spectra for different objects", making it difficult to distinguish the edges of cultivated units; and due to the limitation of the revisit cycle, the reliability of image acquisition in cloudy, rainy and foggy areas is difficult to guarantee. In recent years, UAV remote sensing technology has been widely used in mapping due to its high accuracy, cost-efficiency, and flexibility. The mapping resolution can reach centimeters, which greatly facilitates the visual interpretation of cultivated land patches. However, it relies on manual experience and has low mapping efficiency. In addition, UAV image processing simultaneously obtains orthophotos representing vegetation coverage and digital surface models reflecting field elevation differences, which provides a good opportunity for mining images and elevation information of patchy cultivated land for automatic mapping of cultivated land. Therefore, there is an urgent need for an automatic and accurate extraction method of cultivated land such as patchy cultivated land that can coordinate the orthophotos and elevation features of UAV remote sensing. Summary of the invention

[0004] In view of this, the present invention provides a method for extracting patchy cultivated land boundaries in flat terrain areas by UAV remote sensing, which aims to combine digital image processing methods with cartography principles to achieve rapid and non-artificial extraction of cultivated land areas and realize patchy cultivated land boundary extraction with a high degree of automation, which is conducive to 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 extracting patchy cultivated land boundaries in a flat terrain area by using an unmanned aerial vehicle remote sensing method, characterized by comprising:

[0007] 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;

[0008] 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;

[0009] The step 2 specifically includes the following steps:

[0010] 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:

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

[0012] 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 value in the rectangular sliding window, h min is the minimum elevation value in the rectangular sliding window;

[0013] 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:

[0014]

[0015] 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;

[0016] 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:

[0017]

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

[0019] 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:

[0020]

[0021] 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).

[0022] Step 3: Calculate the vegetation index using the orthophoto obtained in step 1, and binarize it to obtain the second initial cultivated land area;

[0023] 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:

[0024]

[0025] 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;

[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 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):

[0027]

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

[0029] 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;

[0030] 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:

[0031] R1=F3×V2(7)

[0032] 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;

[0033] 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:

[0034]

[0035] 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;

[0036] 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:

[0037]

[0038] Among them, g2 represents the circular structural element, and R4 is the re-optimization of the cultivated land area;

[0039] 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:

[0040]

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

[0042] 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.

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

[0044] 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 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: 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. The shrinkage algorithm parameter value is 0.7, and then the area threshold filtering is used to obtain the final patchy cultivated land boundary extraction result. The threshold of the 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 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.

[0048] In summary, due to the adoption of the above technical solution, the beneficial effects of the present invention are:

[0049] 1. The present invention uses a digital surface model obtained by photogrammetry of a light and small unmanned aerial vehicle to calculate the terrain undulation, realizes the extraction of the field edge area and the initial cultivated land area grid of the cultivated land area, uses the obtained orthophoto to calculate the vegetation index, realizes the extraction of the initial cultivated land area grid, and obtains the intersection pixels of the two to obtain the candidate cultivated land area, and then uses digital image processing technology and cartography principles to realize the extraction of the final cultivated land area.

[0050] 2. The present invention is simple to operate and has high mapping efficiency, and realizes the drone remote sensing extraction of patchy cultivated land boundaries in flat terrain areas supported by drone photogrammetry products and subsequent fine cadastral management.

[0051] 3. The present invention has strong applicability. It is not only applicable to the specific terrain relief and visible light band difference vegetation index in the present invention, but can also be applied to other indicators that can characterize terrain elevation differences and crop coverage status. BRIEF DESCRIPTION OF THE DRAWINGS

[0052] The present invention will now be described by way of example with reference to the accompanying drawings, in which:

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

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

[0055] Figure 3 (a) and (b) are respectively the optimization results of regular-shaped plots, irregular-shaped plots and regular-shaped plots of cultivated land in the study area obtained by combining terrain information and visible light index according to an embodiment of the present invention.

[0056] Figure 4 (a) and (b) are comparison diagrams before and after optimization of irregular-shaped fields in the study area obtained by combining terrain information and visible light index according to an embodiment of the present invention.

[0057] Figure 5 This is a comparison diagram of the patchy farmland boundary extraction result obtained in an embodiment of the present invention and the reference value.

[0058] Figure 6 Schematic diagram of the process structure of this embodiment. DETAILED DESCRIPTION

[0059] In order to make the purpose, 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 drawings in the embodiments of the present invention. Obviously, the described embodiments are part of the embodiments of the present invention, not all of the embodiments. Generally, the components of the embodiments of the present invention described and shown in the 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 invention claimed for protection, but merely represents selected embodiments of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

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

[0062] It should be noted that similar reference numerals and letters denote similar items in the following drawings, and therefore, once an item is defined in one drawing, further definition and explanation thereof is not required in subsequent drawings.

[0063] In the present invention, unless otherwise clearly specified and limited, a first feature being "above" or "below" a second feature may include that the first and second features are in direct contact, or may include that the first and second features are not in direct contact but are in contact through another feature between them. Moreover, a first feature being "above", "above" and "above" a second feature includes that the first feature is directly above and obliquely above the second feature, or simply indicates that the first feature is higher in level than the second feature. A first feature being "below", "below" and "below" a second feature includes that the first feature is directly below and obliquely below the second feature, or simply indicates that the first feature is lower in level than the second feature.

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

[0065] Example 1

[0066] like Figure 6 As shown, an embodiment of the present invention discloses a method for extracting patchy cultivated land boundaries by UAV remote sensing in a flat terrain area, which is characterized by comprising:

[0067] 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;

[0068] 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;

[0069] The step 2 specifically includes the following steps:

[0070] 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:

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

[0072] 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 value in the rectangular sliding window, h min is the minimum elevation value in the rectangular sliding window;

[0073] 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:

[0074]

[0075] 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;

[0076] 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:

[0077]

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

[0079] 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:

[0080]

[0081] 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).

[0082] Step 3: Use the orthophoto obtained in step 1 to calculate the vegetation index to obtain the second initial cultivated land area;

[0083] 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:

[0084]

[0085] 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;

[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 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):

[0087]

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

[0089] 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;

[0090] 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:

[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 cultivated land area by combining F3 and V2;

[0093] 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:

[0094]

[0095] 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;

[0096] 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:

[0097]

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

[0099] 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:

[0100]

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

[0102] 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.

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

[0104] 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 fields and irregular-shaped fields; the area threshold set by the convex hull function is 105% of the pixel area of ​​each field.

[0105] 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. The shrinkage algorithm parameter value is 0.7, and then the area threshold filtering is used to obtain the final patchy cultivated land boundary extraction result. 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 candidate cultivated land areas, R6 represents the extraction result of patchy cultivated land boundary areas, and T3 is the filtering area threshold used.

[0108] Example 2

[0109] This embodiment is a specific implementation of embodiment 1. This embodiment is located about 50km west of Chengdu, Sichuan Province, in the Dujiangyan Irrigation District in the hinterland of the Chengdu Plain. The terrain is flat and the surface elevation is about 550m. The area is a typical two-crop agricultural planting structure, mainly planting wheat, rice and rapeseed; it has a subtropical humid monsoon climate with distinct four seasons, an average annual temperature of 16.1℃, and an average annual rainfall of about 1011mm. The farming conditions of this embodiment after land reclamation are good and convenient for mechanized agricultural production. It is mainly patchy cultivated land with a length of about 50 to 80m and a width of 35 to 60m. The width of the ridge line is about 0.9m, and the height difference relative to the field surface is about 0.6m.

[0110] A method for extracting patchy cultivated land boundaries in a flat terrain area by using a drone remote sensing method in this embodiment includes:

[0111] The multi-temporal UAV images of cultivated land were collected in the early and mature stages of rice growth. Based on the digital surface model and orthophoto obtained by three-dimensional reconstruction, the surface undulation was calculated using the digital surface model, and the threshold segmentation was adopted to binarize it into two types of objects: field edge area and non-field edge area. Then, the area threshold filtering was implemented to obtain a purer binary map of the field edge area. Then, the binary map value of the field edge area was inverted to obtain the initial cultivated land area obtained by using the terrain information. Then, the vegetation index was calculated using the orthophoto and the threshold was binarized to obtain the initial cultivated land area obtained by using the visible light index. Then, the cultivated land area shared by the visible light index and the terrain information was obtained through logical AND operation. Finally, the cultivated land area was divided into regular-shaped fields and irregular-shaped fields according to the area threshold optimized by the convex hull algorithm, and the edge optimization was performed respectively by the convex hull algorithm and the concave hull shrinkage algorithm. The specific implementation process is as follows:

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

[0113] 2) The terrain relief image in December 2023 calculated in step 1) is threshold segmented using the above formula (2), where the threshold is the absolute value of the difference between the mean of the terrain relief and half of the standard deviation, thereby extracting a binary map F1 of the edge area of ​​the candidate cultivated land plot;

[0114] 3) Calculate the area of ​​each patch of the candidate field edge area obtained in step 2), and then use the above formula (3) to perform area threshold filtering to purify the candidate field edge area obtained in step 2) to obtain a field edge area binary map F2;

[0115] 4) The filtered field edge area binary image obtained in step 3) is subjected to a value inversion operation using the above formula (4), thereby obtaining an initial cultivated land area binary image F3 extracted based on the digital surface model using terrain information;

[0116] 5) The orthophoto of August 2023 obtained in step 1) without the near-infrared band is used to calculate the visible-band difference vegetation index (VDVI) based on the vegetation spectral response characteristics and the visible light image band of the drone using the above formula (5) to obtain the visible light image vegetation index V1;

[0117] 6) According to the image histogram of the visible light band difference vegetation index obtained in step 5), determine the binarization adaptive threshold to maximize the variance between vegetation and non-vegetation, and use the above formula (6) to perform threshold segmentation to separate the foreground and background pixels of the image as much as possible, thereby obtaining 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 orthophoto in step 6), the above formula (7) is used to obtain the result R1 of the cultivated land area extracted by combining the terrain information and the visible light index. Based on the obtained R1, the image closing operation (8) is used to obtain the preliminary optimized cultivated land area R2. The holes of the obtained preliminary optimized cultivated land area are filled to obtain the image R3. Then the image opening operation (9) is performed to obtain the re-optimized cultivated land area R4. The area of ​​each patch in the cultivated land area obtained by the image opening operation is calculated. Finally, the candidate cultivated land area binary image R5 is obtained by formula (10);

[0119] 8) According to the regularity and irregularity of the shape of the field, the candidate cultivated land area binary map obtained in step 7) is divided into regular-shaped fields and irregular-shaped fields according to the area threshold set by the convex hull function. For the regular-shaped fields, the convex hull algorithm is used to set the area threshold for optimization. For the irregular-shaped fields, the concave hull shrinkage algorithm is used to set the shrinkage parameters of the field edge point coordinates for optimization. The results of the two optimizations are subjected to logical or operation, and then the area threshold is filtered through formula (11) to obtain the final patchy cultivated land boundary area R6.

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

[0121] Using orthophotos with a spatial resolution of 2.7 cm as the base map, the visual interpretation method was used to extract the cultivated land area, which was used as the verification data. The results showed that the accuracy of cultivated land field extraction based on the number of fields was 93.8%, and the accuracy of cultivated land area extraction based on pixels was 96.8%, corresponding to an overall accuracy of 95.0%, and a Kappa coefficient of 0.89, which proved the effectiveness of the present invention.

[0122] Although 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 the embodiments without departing from the principles and spirit of the present invention, and that 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. Needless to say, the content protected by the present invention does not involve improvements to software and methods.

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

[0125] The above description of the disclosed embodiments enables one skilled in the art to implement or use the present invention. Various modifications to these embodiments will be apparent to one skilled in the art, and the general principles defined herein may 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 widest 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; 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; 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 value 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 1 is characterized in that: 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 the pixel (x, y), and T3 is the segmentation 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: 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 binary map of candidate cultivated land area, and T3 is the filtering area threshold used.

5. The method for extracting patchy cultivated land boundaries by UAV remote sensing in flat terrain areas according to claim 4 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.

6. The method for extracting patchy cultivated land boundaries by UAV remote sensing in flat terrain areas according to claim 5 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.

7. The method for extracting patchy cultivated land boundaries by UAV remote sensing in flat terrain areas according to claim 5 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.

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

7.

9. The method for extracting patchy cultivated land boundaries by UAV remote sensing in flat terrain areas according to claim 5, 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