A method for designing artificial forest thinning operation based on unmanned aerial vehicle image

Abstract

The application discloses a kind of based on unmanned aerial vehicle image artificial forest thinning operation design method, belong to artificial forest thinning operation design method field.The present application first collects target area unmanned aerial vehicle image and reconstructs image generation DOM and DSM data, then using computer image processing method is carried out single tree recognition and crown information extraction, then calculate unit grid forest attribute and carry out hot spot analysis, finally determine thinning area and thinning intensity.Adopt this method to carry out artificial forest thinning operation design, have three aspects of advantages:(1) make up the deficiency of traditional ground survey in the whole range single tree coordinate and stand spatial information;(2) grid subdivision is carried out to whole forest, the spatial heterogeneity in stand is fully considered, and the fine level of thinning design is improved;(3) provide the priority order determination scheme of tree thinning, reduce the subjectivity of thinning process, improve the scientific nature of object wood selection.

Description

TECHNICAL FIELD

[0001] The present application belongs to the field of artificial forest thinning operation design method, and particularly relates to a method for designing artificial forest thinning area and thinning intensity based on high-resolution images of unmanned aerial vehicles. BACKGROUND

[0002] Thinning has been widely used to restore different types of suppressed stands, as it can alleviate the competition pressure between trees, optimize the spatial structure of the stand, and thus achieve sustainable management and development of the stand. Generally, the forest thinning prescription is designed based on field plot investigation, and the forest practitioners subjectively select the trees to be cut according to the given prescription. However, the plot investigation method usually covers a relatively small spatial scale in practice, and due to the spatial heterogeneity of the stand, the accuracy is uncertain when extrapolated to a larger scale. Moreover, this method requires a large amount of manpower and time, especially in areas with large terrain undulations, poor transportation, and large stand density, the workload of investigation is greatly increased, and the accuracy is also linearly decreased. In addition, the current stand thinning prescription usually only specifies the total amount of trees to be cut or the number of trees of each size class or the percentage of cross-sectional area, resulting in a large number of spatial patterns of trees to be cut that meet the same thinning requirements of the stand. It is difficult for forest practitioners to evaluate the pros and cons of each selection mode, which may lead to large differences in the distribution of trees after thinning. Moreover, forest practitioners usually pursue "simple" tree extraction (i.e., close to the roadside or logging yard) as long as it meets the thinning requirements, without considering the overall regulation of the spatial structure of the stand.

[0003] The above limitations are mainly due to the lack of single-tree level information covering the entire stand that can be used to develop detailed thinning prescriptions. In recent years, the rapid development of unmanned aerial photography technology and image processing technology has created favorable conditions for capturing the spatial variability of trees within the stand and generating tree distribution maps and tree attributes. In order to solve the limitations of current thinning design practices, the present application fully considers the spatial heterogeneity of the stand, analyzes the spatial distribution pattern of trees and the compression of trees, and uses a computer to assist decision-making to develop a more detailed and accurate thinning operation plan. SUMMARY

[0004] The present application aims to overcome the defects in the prior art and provide an artificial forest thinning operation design method based on unmanned aerial vehicle images. The present application focuses on two key points of thinning operation design, thinning area and thinning intensity design, fully considers the spatial heterogeneity of the stand and the spatial distribution pattern of trees, solves the problems of where to cut and how much to cut, and provides technical support for optimizing thinning operation practices.

[0005] The specific technical solutions adopted by the present application are as follows:

[0006] The application provides a forest thinning operation design method based on unmanned aerial vehicle images, and specifically as follows.

[0007] S1, image collection is performed on a target thinning operation area;

[0008] S2, three-dimensional reconstruction is performed on the thinning operation area images collected in S1 to obtain an orthographic image and a digital surface model;

[0009] S3, based on the digital surface model obtained in S2, a local maximum algorithm is used to identify the positions of individual trees;

[0010] S4, according to the orthographic image obtained in S2, a grid computer is used to calculate VDVI pixel values; and based on the green light band pixel values and the VDVI pixel values, a threshold segmentation method is used to draw a forest canopy projection area;

[0011] S5, the thinning operation area images obtained in S1 are divided into a series of continuous grid blocks with spatial properties according to a user-defined grid size;

[0012] S6, according to the positions of individual trees obtained in S3 and the forest canopy projection area obtained in S4, the stock tree density SD of the unit grid in S5 is calculated j and the canopy density CD of the unit grid is calculated j .

[0013] S7, based on the results obtained in S6, a desired thinning threshold ET defined by a user is used to determine a thinning area;

[0014] S8, based on the results obtained in S6, a desired post-thinning density ED defined by a user is used to calculate a thinning intensity P.

[0015] As a preferred, in the step S1, the image collection is realized by means of a low-cost consumer unmanned aerial vehicle with an RGB camera.

[0016] As a preferred, in the step S4, the calculation formula of the VDVI pixel value is as follows:

[0017] VDVI = (2 × P green -P red -P blue ) / (2 × P green + P red + P blue )

[0018] Wherein, P green , P red , P blue represent the pixel values of the green light, red light and blue light bands respectively.

[0019] As a preferred, the step S7 is specifically as follows:

[0020] 1) The number density and canopy density of the strains obtained in S6 are analyzed for cold and hot spots, and the target stand is divided into three levels of high-density aggregation area, insignificant clustering area and low-density aggregation area according to the weight w value; wherein, w = 1 is the high-density aggregation area, w = 0 is the insignificant clustering area, and w = -1 is the low-density aggregation area;

[0021] 2) The two layers in step 1) are added according to the weight w to obtain a layer with w values of -2, -1, 0, 1, 2 five levels;

[0022] 3) Calculate the number density or canopy density of each level area, and compare the expected thinning threshold ET to determine the thinning area that exceeds ET.

[0023] As preferred, in the step S8, the calculation formula of the thinning intensity P is as follows:

[0024] P = (D 现 - ED) / D 现

[0025] Wherein, P represents the thinning intensity, D 现 and ED represent the number density or canopy density after the current and expected thinning, respectively.

[0026] Compared with the prior art, the present application has the following beneficial effects:

[0027] (1) It makes up for the deficiency of traditional ground survey in the whole range of single tree coordinates and stand spatial information;

[0028] (2) The whole forest is subdivided into grids, fully considering the spatial heterogeneity within the stand, and improving the fine level of thinning design;

[0029] (3) The priority order determination scheme of tree thinning is provided, which reduces the subjectivity of the thinning process and improves the scientificity of the object wood selection. BRIEF DESCRIPTION OF DRAWINGS

[0030] Figure 1 Flowchart of the present application;

[0031] Figure 2 (a) Single tree spatial distribution map; (b) canopy spatial distribution map;

[0032] Figure 3 Strain density and canopy density cold and hot spot analysis result map;

[0033] Figure 4 Five-level spatial distribution map of thinning area;

[0034] Figure 5 Thinning intensity map of each level area and grid area. DETAILED DESCRIPTION

[0035] The present application will be further illustrated and described with reference to the accompanying drawings and detailed description. The technical features of each embodiment of the present application can be combined accordingly without conflict.

[0036] As shown in the drawings, the present application provides a design method for artificial forest thinning operation based on unmanned aerial vehicle image, which specifically comprises the following steps: Figure 1

[0037] S1, image collection is performed on a target thinning operation area.

[0038] In actual application, image collection can be realized by means of a low-cost consumer drone equipped with an RGB camera.

[0039] S2, three-dimensional reconstruction is performed on the image of the thinning operation area collected in S1 to obtain a digital orthophoto map (DOM) and a digital surface model (DSM).

[0040] S3, based on the digital surface model obtained in S2, a local maximum algorithm is used to identify the position of a single tree.

[0041] S4, according to the digital orthophoto map obtained in S2, a raster computer is used to calculate the VDVI pixel value; and based on the green light band pixel value and the VDVI pixel value, a threshold segmentation method is used to draw the crown projection area.

[0042] Specifically, the calculation formula of the VDVI pixel value is as follows:

[0043] VDVI = (2 x P green -P red -P blue ) / (2 x P green + P red + P blue )

[0044] Wherein, P green , P red , P blue represent the pixel values of the green light, red light and blue light bands respectively.

[0045] S5, the image of the thinning operation area obtained in S1 is divided into a series of continuous grid blocks with spatial properties according to the user-defined grid size.

[0046] S6, according to the position of the single tree obtained in S3 and the crown projection area obtained in S4, the tree density SD j and the canopy density CD j of the unit grid in S5 are calculated respectively.

[0047] ​S7, based on the result of S6, determine the thinning area according to the user-defined expected thinning threshold ET (by plant density or canopy density).

[0048] In actual use, the step is specifically as follows:

[0049] 1) Perform cold-hot spot analysis on the plant density and canopy density obtained in S6 respectively, and divide the target stand into three levels of high-density aggregation area (w=1), clustering insignificant area (w=0), and low-density aggregation area (w=-1) according to the weight w value;

[0050] 2) Add the two layers in step 1) according to the weight w to obtain a layer with five levels of w values of -2, -1, 0, 1, and 2;

[0051] 3) Calculate the plant density or canopy density of each level area, compare the expected thinning threshold ET, and determine the area exceeding ET as the thinning area.

[0052] S8, based on the result of S6, calculate the thinning intensity P according to the user-defined expected post-thinning density ED (by plant density or canopy density).

[0053] Specifically, the calculation formula of the thinning intensity P is as follows:

[0054] P=(D 现 -ED) / D 现

[0055] Wherein, P represents the thinning intensity, D 现 and ED represent the plant density or canopy density before and after the expected thinning respectively.

[0056] Embodiment

[0057] The embodiment provides a method for artificial forest thinning operation design based on unmanned aerial vehicle image, and the specific steps are as follows:

[0058] S1 and S2, target area image acquisition and pretreatment: 30-year-old single-layer pure Chinese fir plantations in Qianchuan Town, Lin'an District are selected as test forests, the test forest covers 1.56 hectares, the stand density is about 1100 plants per hectare, and the average canopy density is about 0.8. In November 2018, a low-cost consumer unmanned aerial vehicle (model: DJI Phantom 4pro) was used to collect a total of 82 images required for the test, wherein the longitudinal and transverse overlap degrees of flight are 90% and 85% respectively. Image preprocessing is performed by agisoft photoscan (version 1.2.5.2594) through a standard workflow to generate a digital orthophoto map (DOM) and a digital surface model (DSM).

[0059] S3, Local maximum single tree recognition algorithm based on DSM: This algorithm uses a sliding window, which is approximately equal to the average crown size, to traverse the entire DSM from left to right and from top to bottom. The unique local maximum, i.e., the tree tip position, can be found by iteration from the sliding window at each position. This step realizes the automatic recognition of the position of each tree from the UAV DSM image, as shown in Fig. a, and the spatial distribution of single trees is accurately marked on the map in the form of points. Figure 2

[0060] S4, Threshold segmentation method based on G and VDVI: 1) Extract the green band (G) from the DOM, and calculate the visible difference vegetation index (VDVI) using a grid calculator, where VDVI = (2 x P green -P red -P blue ) / (2 x P green + P red + P blue ), P green , P red , and P blue represent the pixel values of the green, red, and blue bands, respectively; 2) The OTSU maximum between-class variance method is used to determine the segmentation threshold values of the green band and VDVI, respectively, and the non-crown area is extracted; 3) The non-crown parts of the two parts are combined to obtain the non-crown coverage area of the whole forest; 4) The crown coverage area is obtained by subtracting the non-crown coverage area from the whole forest area. This step realizes the extraction of the crown coverage area from the UAV DOM image, as shown in Fig. b, and the spatial distribution of the crown projection area is drawn on the map in the form of a plane. Figure 2

[0061] S5, Unit grid division: Based on the preliminary analysis of the forest stand, a 10m x 10m square is determined as the basic unit. In ArcGIS, the southwest corner of the image is taken as the starting point to create a 10m x 10m fishing net to refine the forest stand.

[0062] S6, Forest attribute calculation based on grid: According to the spatial distribution map of single trees obtained in S3, the tree density SD j of the unit grid is calculated, and according to the crown coverage area distribution map obtained in S4, the canopy density CD j of the unit grid is calculated.

[0063] S7, Thinning area determination: 1) The expected thinning threshold ET of 0.8 of the canopy density is user-defined in this test; 2) As shown in Fig. Figure 3 ​​As shown, the cold and hot spot analysis of the plant number density and canopy density respectively, are divided into three levels of high value aggregation area, clustering not significant area, low value aggregation area, and are assigned 1, 0, -1 respectively, and the distribution of cold and hot spots of the two is basically the same; 3) Add the two layers to get a layer with values of 2, 1, 0, -1, -2 in total 5 levels; 4) Calculate the average canopy density of each level area, compared with the expected threshold ET of thinning, if the canopy density is greater than 0.8, it is determined as the thinning area, otherwise not. As shown in Figure 4 The first and second level areas of the test area have a canopy density greater than 0.8, so they are determined as the thinning area.

[0064] S8, Thinning intensity determination: 1) The expected post-thinning density ED of 0.7 canopy density is defined by the user in this test; 2) According to the formula P = (D 现 -ED) / D 现 , calculate the thinning intensity, where P represents the thinning intensity, D 现 and ED represent the plant number density or canopy density before and after the expected thinning respectively. As shown in Figure 5 The green line frame defines the range of the final thinning required, and the thinning intensity of the first level area is 19.3%, and the thinning intensity of the second level area is 16.7%, and the thinning intensity of the unit grid size is also marked in the corresponding position.

[0065] The above-described embodiments are only a preferred scheme of the present application, and are not intended to limit the present application. Those skilled in the art can make various changes and modifications without departing from the spirit and scope of the present application. Therefore, any technical solutions obtained by equivalent replacement or equivalent transformation shall fall within the protection scope of the present application.

Claims

1. A method for designing thinning operations in artificial forests based on UAV imagery, characterized in that, Specifically as follows: S1, image collection is performed on a target intermediate cutting operation area; S2, three-dimensional reconstruction is performed on the image of the intermediate cutting operation area collected in S1 to obtain an orthographic image and a digital surface model; S3, based on the digital surface model obtained in S2, a local maximum value algorithm is used to identify the position of a single tree; S4. The orthographic image obtained in S2 is calculated using a raster computer VDVI pixel values; and based on the green band pixel values and VDVI pixel values, a threshold segmentation method is used to draw the canopy projection area; S5, the image of the intermediate cutting operation area obtained in S1 is divided into a series of continuous grid blocks with spatial properties according to a user-defined grid size; S6, calculate the tree density of unit grid in S5 according to the single tree position obtained in S3 and the crown projection area obtained in S4 SD j and the canopy density of unit grid CD j ; S7、based on the result obtained in S6, determining the gap cutting region according to the expected gap cutting threshold value defined by the user ET . S8, based on the result obtained in S6, calculating the intensity of the intermediate cutting P according to the expected intermediate post-cutting density defined by the user ED ; The step S7 is specifically as follows: 1) Perform hot and cold spot analysis on the plant density and canopy closure obtained from S6, and then analyze them according to their weights. w The value divides the target forest stand into three levels: high-density clustering areas, areas with insignificant clustering, and low-density clustering areas; among them, w = 1 indicates a high-density clustering area. w =0 indicates a region where clustering is not significant. w = -1 indicates a low-density clustering area; 2) add the two layers in step 1) by weight w to obtain w Layers with values of -2, -1, 0, 1, 2 3) Calculate the number of trees per hectare or canopy density for each rank area, compare with the expected thinning threshold ET areas that exceed ET the threshold are identified as thinning areas. 2.The method of designing thinning operation of a plantation based on UAV images according to claim 1, wherein, In the step S1, image collection is realized by means of a low-cost consumer drone with an RGB camera.

3. The method according to claim 1, wherein, In the step S4, VDVI The formula for calculating the pixel value is as follows: VDVI=(2×P green -P red - P blue ) / (2×P green +P red +P blue ) wherein, P green , P red , P blue represent the pixel values of the green, red and blue light bands, respectively.

4. The method according to claim 1, wherein, In the step S8, the intensity of thinning P The calculation formula is as follows: P =(D 现 - ED) / D 现 wherein, P represents the intensity of the intermediate cutting, D 现 and ED respectively represent the current and expected number density of trees or canopy density after intermediate cutting.

Images