A method for deploying forest permanent large sample plots using RTK

Abstract

The present invention discloses a method for arranging forest permanent plots using RTK, which relates to the technical field of forestry surveys. The aim is to solve the technical problems that the current methods of arranging points for forest permanent plots, such as manual point arrangement and tree species identification and marking for point arrangement, are time-consuming and laborious and cannot effectively improve efficiency and accuracy. The method includes the following steps: S1, point cloud data acquisition; using a drone equipped with lidar, airborne laser scanning is carried out above the original terrain, high-precision position coordinates are obtained through GNSS / INS combined positioning, and point cloud data is collected; S2, three-dimensional digital model establishment; the point cloud data is processed to filter out noise points and remove isolated points, and a digital elevation model is generated based on the processed point cloud data, and the center points of forest permanent plots are generated in the digital elevation model; S3, on-site plot planning; the layout area of the permanent plots is divided into several on-site plots. The present invention has the advantages of improving positioning accuracy and plot layout efficiency.

Description

Technical Field

[0001] The present invention relates to the technical field of forestry surveys, and more specifically, to a method for laying out fixed large forest plots using RTK. Background Art

[0002] With the in-depth study of the global forest ecosystem, there is an urgent need to establish a global forest monitoring system. Remote sensing technology plays an increasingly important role in the dynamic monitoring of forest resources and the research of forest ecosystems. The dynamic monitoring of forest resources mainly includes the real-time and regular monitoring of resource parameters such as forest volume, area, tree species distribution, and canopy.

[0003] In traditional forest ecology plot research, ground plots are established and corresponding geostatistical methods are applied to conduct small-scale forest resource surveys to obtain characteristic data of each tree species within the plots, and then models are established to simulate the dynamic changes of forest resources. The layout method of fixed forest plots mainly relies on traditional forest resource survey methods. This method requires a large amount of human and time costs and is not convenient for implementation within large ecosystems.

[0004] Using remote sensing image analysis methods to conduct dynamic monitoring of forest resources can save the human and time costs of forest surveys and improve efficiency. It is the main method at present. However, in the remote sensing analysis of forest resources, it is difficult to obtain accurate ground verification data using traditional methods, and there are also different opinions on methods such as image extraction. Especially when dealing with complex terrains such as hilly terrains, the analysis of remote sensing images cannot effectively obtain the true values of the images. Therefore, it is very necessary to apply the research of fixed forest ecology plots to strengthen the dynamic monitoring of forest resources.

[0005] At present, there are mainly two methods for laying out fixed forest ecology plots: manual point layout and tree species identification and marking layout. Manual point layout requires referring to previous data, actual terrain conditions, and the representativeness of sample points, and requires rich experience and professional knowledge, consuming a large amount of human and time; tree species identification and marking layout requires manual visual inspection of markers, which are special features that can be clearly identified from the air and require certain area, shape, etc. characteristics. It requires special investment in labor and time, and requires repeated field identification and marking, which is time-consuming and laborious and cannot effectively improve efficiency and accuracy.

[0006] In view of this, we propose a method for laying out fixed large forest plots using RTK. Summary of the Invention

[0007] The purpose of the present invention is to provide a method for laying out fixed large forest plots using RTK to solve the technical problems that the current layout of fixed large forest plots by manual layout and tree species identification and marking layout is time-consuming and laborious and cannot effectively improve efficiency and accuracy.

[0008] To solve the above technical problems, the present invention provides the following technical solutions: A method for deploying forest permanent plots using RTK, comprising the following steps:

[0009] S1. Point cloud data acquisition;

[0010] Using a drone equipped with a lidar, airborne laser scanning is carried out above the original terrain, high-precision position coordinates are obtained through GNSS / INS combined positioning, and point cloud data is collected;

[0011] S2. Three-dimensional digital model establishment;

[0012] Process the point cloud data, filter out noise points and remove isolated points, and generate a digital elevation model based on the processed point cloud data, and generate the center point of the forest permanent plot in the digital elevation model;

[0013] S3. Field plot planning;

[0014] Divide the deployment area of the permanent plot into several field plots, and determine the field plot range and the center coordinates of the field plot;

[0015] S4. Obtain the center of the drone's extended flight route;

[0016] Based on the center coordinates of the field plot, the center coordinates of the drone's extended flight route are calculated inversely, and the extended flight altitude is determined according to the area of the field plot and the height difference between the highest point and the lowest point in the permanent plot, and an extended flight route is generated;

[0017] S5. Obtain the coordinates and elevation of the extended sample points;

[0018] Through the drone, airborne laser scanning is carried out along the extended flight route to obtain the coordinates and elevation of the extended sample points;

[0019] S6. Deploy the permanent plot;

[0020] According to the coordinates and elevation of the extended sample points, the coordinates of the control sample points of the permanent plot are calculated inversely, and the permanent plot is deployed according to the coordinates of the control sample points.

[0021] Preferably, in S2, the drone is a fixed-wing drone, and the digital elevation model generation method is based on the construction of an irregular triangular network TIN, and the generation method is as follows: First, perform Delaunay triangulation on the processed point cloud data to construct an irregular triangular network. Let any triangle T abc in the triangular network, and its three vertex coordinates are A(x a , y a , z a ), B(x b , y b,z b ), C(x c , y c , z c ), for any point P(x, y) inside the triangle, its elevation z p is calculated by linear interpolation:

[0022]

[0023] By traversing the entire triangular network area, a digital elevation model is generated;

[0024] The processed point cloud data is mapped to the plane of the digital ground elevation model to establish a digital ground model. The original terrain is thinned in the three-dimensional digital model of the area, and the digital height of the digital ground model is added to the digital ground points in the thinned digital ground model.

[0025] Preferably, the method for planning field plots on the three-dimensional digital model of the original terrain and back-calculating the center of the UAV's expanded flight path through the center of the field plot includes:

[0026] a1. Based on the three-dimensional digital model, divide the layout area of the fixed large plot into several field plots;

[0027] a2. Determine the scope of the field plot, which includes the area and boundary of the field plot;

[0028] a3. Obtain the center coordinates of the field plot;

[0029] a4. Based on the center coordinates of the field plot, obtain the center coordinates of the UAV's expanded flight path.

[0030] Preferably, the method for back-calculating the center of the UAV's expanded flight path through the center of the field plot is based on the fact that within the scope of the field plot, the center of the fixed large plot is back-calculated through the center coordinates of the field plot. Let the center coordinates of the field plot be (X si , Y si ), where i = 1, 2,..., m, and m is the number of field plots. First, calculate the average value of all field plot center coordinates in the X direction and the average value in the Y direction The formulas are as follows:

[0031]

[0032] Then, taking as the initial estimate value, based on the coordinates of the expanded sample points (X spj , Y spj , H spj), j = 1, 2, …, n, where n is the number of extrapolated sample points, and weighted least squares adjustment optimization is performed;

[0033] Let the distance from the extrapolated sample point to the initial estimated center be Construct the error equation:

[0034]

[0035] where v Xj , v Yj are the observation residuals in the X and Y directions, δX and δY are the corrections of the center coordinates in the X and Y directions, and Δd j is the correction of the distance from the extrapolated sample point to the optimized center coordinates;

[0036] According to the least squares principle Construct the normal equation to solve for δX and δY, where p j is the weight of the extrapolated sample point, determined according to the accuracy of the point, and construct the normal equation to solve for δX and δY;

[0037] Finally, the calculation formula for fixing the center coordinates (X c , Y c ) of the large sample plot is:

[0038] Overlay the fixed center coordinates of the large sample plot with the three-dimensional digital model of the original terrain to obtain the center coordinates of the extrapolated flight path;

[0039] Let the center coordinates of the actual sample plot in the field be (X si , Y si ), the offset in the X direction of the fixed large sample plot center relative to the actual sample plot center is ΔX, and the offset in the Y direction is ΔY. Then, the calculation formula for the center coordinates (X o , Y o ) of the UAV extrapolated flight path is:

[0040] X o = X si + ΔX, Y o = Y si + ΔY.

[0041] Preferably, the method for determining the UAV extrapolated flight altitude through the center of the extrapolated flight path is to determine the height difference between the highest point and the lowest point in the fixed large sample plot based on the actual altitude of the sample plot, determine the distance from the lowest point to the center of the extrapolated flight path according to the height difference, and set the flight altitude as the elevation value from the lowest point in the fixed large sample plot to the center of the extrapolated flight path according to this distance;

[0042] Let the center coordinates of the actual sample plot in the field be Let the elevation of the highest point in the fixed large sample plot be H max , and the elevation of the lowest point be H min, if the horizontal distance from the lowest point in the sample plot to the center of the expanded flight line is D, then the expanded flight altitude H flight The calculation formula is: where ΔH = H max -H min .

[0043] Preferably, the expanded flight line expands outward from the center of the expanded flight line, and the scanning points are the expanded sample points. The method for calculating the coordinates and elevation of the expanded sample points from the expanded sample points includes:

[0044] b1. Superimpose the center coordinates of the fixed large sample plot with the three-dimensional digital model of the original terrain, align the center coordinates of the fixed large sample plot with the expanded flight line to obtain the expanded flight line;

[0045] b2. Use the unmanned aerial vehicle to fly at the elevation flight altitude according to the expanded flight line for expanded scanning to obtain the coordinates and elevation of the expanded sample points;

[0046] Let the center coordinates of the expanded flight line be (X o , Y o ). The offset distance of the expanded sample point relative to the center of the expanded flight line in the X direction is ΔX spj , and the offset distance in the Y direction is ΔY spj . The radius of the earth is R, and the angle between the line connecting the expanded sample point and the center of the expanded flight line and the due east direction is θ. Then the calculation formula for the coordinates (X spj , Y spj ) of the expanded sample point is:

[0047] X spj = X o + ΔX spj cosθ - ΔY spj sinθ

[0048] Y spj = Y o + ΔX spj sinθ + ΔY spj cosθ

[0049] where ΔX spj and ΔY spj are determined according to the set expanded distance interval and the serial number of the expanded sample point. Let the expanded step length be d, and the serial number of the f-th expanded sample point in this direction be f. Then ΔX spj = f × d × cosθ, ΔY spj = f × d × sinθ;

[0050] Let the expanded flight altitude be H flight , and the elevation of the center of the expanded flight line be H o, the vertical height deviation of the extrapolated sample point relative to the center of the extrapolated flight line is ΔH spi , then the elevation H of the extrapolated sample point spi The calculation formula is: H spi = H o + H flight + ΔH spi .

[0051] Preferably, the method for calculating the coordinates and elevation of the extrapolated sample point by extrapolating the extrapolated sample point is to obtain the distance deviation based on the extrapolated sample point and the center of the extrapolated flight line, determine the course distance between adjacent extrapolated sample points based on the distance deviation, and determine the longitude and latitude of the extrapolated sample point according to the course distance;

[0052] Let the coordinates of the center of the extrapolated flight line be (X o , Y o , H o ), the distance deviation between the extrapolated sample point and the center of the extrapolated flight line is Δd, and the course distance between adjacent extrapolated sample points is d h , and the radius of the earth is R, then the longitude and latitude of the extrapolated sample point are calculated as follows:

[0053] In the latitude direction: In the longitude λ direction:

[0054] where is the latitude of the extrapolated sample point, and λ sp is the longitude of the extrapolated sample point the latitude of the center of the extrapolated flight line, and λ o is the longitude of the center of the extrapolated flight line.

[0055] Preferably, the coordinates of the fixed large sample plot control sample point are calculated by inversely calculating the coordinates and elevation of the extrapolated sample point obtained by extrapolating the extrapolated sample point. The calculation method is as follows:

[0056] Let the set of coordinates of the extrapolated sample point be Based on the coordinates of the extrapolated sample point, perform spatial fitting, and use the least squares method to fit a plane equation Z = aX + bY + c, where a, b, and c are fitting parameters, and the fitting parameters are calculated by the following formula:

[0057]

[0058] Solve for a, b, and c;

[0059] Let the fixed large sample plot be a rectangle with a length of L and a width of W. Based on the coordinates of the center point of the large sample plot (X c , Y c , H c ) on the fitting plane, set m control sample points at equal intervals in the long direction i, and the interval Set k control sample points at equal intervals in the width direction j. Then the coordinates (X cpij , Y cpij , H cpij ) of the control sample point in the j-th row and i-th column are calculated as follows:

[0060] X cpij = X c + (i - 1)ΔX1, where i = 1, 2, …, m

[0061] Y cpij = Y c + (j - 1)ΔY1, where j = 1, 2, …, k

[0062] H cpij = aX cpij + bY cpij + c

[0063] Set up pole markers at the control sample points through a GNSS instrument, obtain the coordinates of the pole markers, and set the layout direction of the pole markers to be consistent with the fixed ground sample plot.

[0064] Preferably, it further includes the calculation of the accuracy of the pole marker layout coordinates. Calculate the accuracy value between the coordinates of the pole marker and the control sample point through the outer-expanded sample points and the outer-expanded flight line. The calculation method is as follows: Let the theoretical coordinates of the control sample point deduced from the outer-expanded sample points be (X ct , Y ct , H ct ), the actual coordinates of the pole marker be (X rb , Y rb , H rb ), and the accuracy threshold be ε. Then the formula for calculating the accuracy value P is:

[0065]

[0066] When P < ε, the pole marker layout meets the accuracy requirements; when P > ε, the pole marker layout does not meet the accuracy requirements.

[0067] Compared with the prior art, the beneficial effects of the present invention are:

[0068] 1. The present invention collects point cloud data by using a drone equipped with lidar, and successively calculates the center coordinates of the outer-expanded flight line of the drone, determines the outer-expanded flight altitude, obtains the coordinates and elevations of the outer-expanded sample points, and finally calculates the coordinates of the control sample points of the fixed large sample plot by inversion, greatly reducing the labor input and being beneficial to improving the positioning accuracy and the sample plot layout efficiency.

[0069] 2. The present invention can also effectively verify the accuracy of the layout of the rod markers by calculating the accuracy values of the rod marker coordinates and the control sample point coordinates and comparing them with the set accuracy threshold, thereby ensuring the reliability of the data for the layout of the entire sample plot.

[0070] 3. The present invention can also obtain rich topographic and ground object information through the point cloud data collected by the unmanned aerial vehicle. Compared with the limited field measurement data that may be relied on in the prior art or the remote sensing image data that is difficult to accurately interpret, the point cloud data can more comprehensively reflect the actual situation of the sample plot and provide more complete data support for forest resource investigation and research. Detailed implementation manners

[0071] To facilitate the understanding of the technical solution of the present invention by those skilled in the art, the technical solution of the present invention will be further described below.

[0072] A method for laying out forest permanent large sample plots using RTK includes the following steps:

[0073] S1. Point cloud data collection;

[0074] Use an unmanned aerial vehicle equipped with a lidar to perform airborne laser scanning above the original terrain, obtain high-precision position coordinates through GNSS / INS combined positioning, and collect point cloud data;

[0075] S2. Three-dimensional digital model establishment;

[0076] Process the point cloud data, filter out noise points and remove isolated points, and generate a digital elevation model based on the processed point cloud data, and generate the center point of the forest permanent large sample plot in the digital elevation model;

[0077] S3. Field sample plot planning;

[0078] Divide the layout area of the permanent large sample plot into several field sample plots, and determine the field sample plot range and the center coordinates of the field sample plot;

[0079] S4. Obtain the center of the outer expansion flight path of the unmanned aerial vehicle;

[0080] Based on the center coordinates of the field sample plot, calculate the center coordinates of the outer expansion flight path of the unmanned aerial vehicle in reverse, and determine the outer expansion flight altitude according to the area of the field sample plot and the height difference between the highest point and the lowest point in the permanent large sample plot, and generate an outer expansion flight path;

[0081] S5. Obtain the coordinates and elevations of the outer expansion sample points;

[0082] Perform airborne laser scanning along the outer expansion flight path by the unmanned aerial vehicle to obtain the coordinates and elevations of the outer expansion sample points;

[0083] S6. Lay out the permanent large sample plot;

[0084] According to the coordinates and elevations of the extended sample points, the coordinates of the control sample points of the fixed large sample plot are calculated inversely, and the fixed large sample plot is arranged according to the coordinates of the control sample points.

[0085] In the embodiment of the present invention, the method for filtering noise points from the point cloud data is based on a statistical filtering algorithm, and the algorithm formula is as follows: Let a certain point P i (x i , y i , z i ) in the point cloud data, calculate its average distance to k surrounding neighborhood points and the average distance mean μ d and the standard deviation σ d ;

[0086]

[0087] If then it is determined that point P i is a noise point and is filtered out, where n is a set threshold coefficient, taking 2-3;

[0088] The method for removing isolated points from the point cloud data is based on a clustering analysis algorithm. For any point P i in the point cloud data, calculate the number of neighborhood points N(P i , r) within a radius r. If N(P i , r) < m, then point P i is marked as an isolated point and removed, where m is the threshold of the minimum number of clustering points.

[0089] In the embodiment of the present invention, in S2, the unmanned aerial vehicle is a fixed-wing unmanned aerial vehicle, and the method for generating a digital elevation model is based on the construction of an irregular triangular network TIN, and the generation method is as follows: First, perform Delaunay triangulation on the processed point cloud data to construct an irregular triangular network. Let any triangle T abc in the triangular network, and the coordinates of its three vertices are A(x a , y a , z a ), B(x b , y b , z b ), C(x c , y c , z c ). For any point P(x, y) inside the triangle, its elevation z p is calculated by linear interpolation:

[0090]

[0091] By traversing the entire triangular network area, a digital elevation model is generated;

[0092] Set digital ground points and the plane of the digital ground elevation model on the point cloud data, map the digital ground points to the plane of the digital ground elevation model, establish a digital ground model, thin out the original terrain in the three-dimensional digital model of the area, and add the digital height of the digital ground model to the digital ground points in the thinned digital ground model.

[0093] The original terrain is low mountain and hilly terrain, and the flight course overlap degree of collecting lidar information above the original terrain is not less than 60%. The high overlap degree can effectively reduce the adjacent two flight course overlap areas.

[0094] In the embodiment of the present invention, the method for field sample plot planning for the three-dimensional digital model of the original terrain and calculating the center of the UAV's expanded flight route by the center of the field sample plot includes:

[0095] a1. Based on the three-dimensional digital model, divide the layout area of the fixed large sample plot into several field sample plots;

[0096] a2. Determine the scope of the field sample plot, and the scope of the field sample plot includes the area of the field sample plot and the boundary of the field sample plot;

[0097] a3. Obtain the center coordinates of the field sample plot;

[0098] a4. Based on the center coordinates of the field sample plot, obtain the center coordinates of the UAV's expanded flight route.

[0099] In the embodiment of the present invention, the area of the field sample plot is 20m×20m. The selection of the side length of the field sample plot is comprehensively considered according to the existing sample plot layout. If the sample plot area is small, it is easy to cause loss of point data. If the sample plot area is large, it is easy to cause the sample plot shape to be irregular. If the area of the field sample plot is too large, it will affect the sampling accuracy in the later field investigation. The area of the field sample plot is set within 20m×20m based on the actual investigation area and the terrain of the actual investigation area. The method for calculating the center of the UAV's expanded flight route by the center of the field sample plot is to calculate the center of the fixed large sample plot by the center coordinates within the scope of the field sample plot. Let the center coordinates of the field sample plot be (X si ,Y si ), where i = 1, 2, …, m, and m is the number of field sample plots. First, calculate the average value of all the center coordinates of the field sample plots in the X direction and the average value in the Y direction The formulas are as follows:

[0100]

[0101] Then, using as the initial estimate value, based on the coordinates of the expanded sample points (X spj ,Y spj ,H spj), j = 1, 2, …, n, where n is the number of extrapolated sample points, and weighted least squares adjustment optimization is performed;

[0102] Let the distance from the extrapolated sample point to the initial estimated center be Construct the error equation:

[0103]

[0104] where v Xj and v Yj are the observation residuals in the X and Y directions, δX and δY are the corrections of the center coordinates in the X and Y directions, and Δd j is the correction of the distance from the extrapolated sample point to the optimized center coordinates;

[0105] According to the least squares principle Construct the normal equation to solve for δX and δY, where p j is the weight of the extrapolated sample point, determined according to the accuracy of the point, and construct the normal equation to solve for δX and δY;

[0106] Finally, the calculation formula for fixing the center coordinates (X c , Y c ) of the large sample plot is:

[0107] Overlay the fixed center coordinates of the large sample plot with the three-dimensional digital model of the original terrain to obtain the center coordinates of the extrapolated flight path;

[0108] Let the center coordinates of the actual sample plot on the ground be (X si , Y si ), the offset in the X direction of the fixed large sample plot center relative to the actual sample plot center on the ground is ΔX, and the offset in the Y direction is ΔY. Then the calculation formula for the center coordinates (X o , Y o ) of the UAV extrapolated flight path is:

[0109] X o = X si + ΔX, Y o = Y si + ΔY.

[0110] In the embodiment of the present invention, the method for determining the UAV extrapolated flight altitude through the center of the extrapolated flight path is to determine the height difference between the highest point and the lowest point in the fixed large sample plot based on the actual altitude of the sample plot, determine the distance from the lowest point to the center of the extrapolated flight path according to the height difference, and set the flight altitude as the elevation value from the lowest point in the fixed large sample plot to the center of the extrapolated flight path;

[0111] Let the center coordinates of the actual sample plot on the ground be Let the elevation of the highest point in the fixed large sample plot be H max , and the elevation of the lowest point be H min, if the horizontal distance from the lowest point in the sample plot to the center of the expanded flight line is D, then the expanded flight altitude H flight The calculation formula is: where ΔH = H max - H min .

[0112] In the embodiment of the present invention, the method for calculating the coordinates and elevation of the expanded sample points by the expanded sample points includes:

[0113] b1. Superimpose the center coordinates of the fixed large sample plot with the three-dimensional digital model of the original terrain, align the center coordinates of the fixed large sample plot with the expanded flight line to obtain the expanded flight line;

[0114] b2. Use the unmanned aerial vehicle to fly at the elevation flight altitude according to the expanded flight line for expanded scanning to obtain the coordinates and elevation of the expanded sample points;

[0115] Let the center coordinates of the expanded flight line be (X o , Y o ), the offset distance of the expanded sample point relative to the center of the expanded flight line in the X direction is ΔX spj , and the offset distance in the Y direction is ΔY spj . The radius of the earth is R, and the angle between the line connecting the expanded sample point and the center of the expanded flight line and the due east direction is θ. Then the calculation formula for the coordinates (X spj , Y spj ) of the expanded sample point is:

[0116] X spj = X o + ΔX spj cosθ - ΔY spj sinθ

[0117] Y spj = Y o + ΔX spj sinθ + ΔY spj cosθ

[0118] where ΔX spj and ΔY spj are determined according to the set expanded distance interval and the serial number of the expanded sample point. Let the expanded step length be d, and the serial number of the f-th expanded sample point in this direction be f. Then ΔX spj = f × d × cosθ, ΔY spj = f × d × sinθ;

[0119] Let the expanded flight altitude be H flight , the elevation of the center of the expanded flight line be H o , and the vertical height deviation of the expanded sample point relative to the center of the expanded flight line be ΔHspi , then the elevation of the extended sample point H spi The calculation formula is: spi =H o +H flight +ΔH spi .

[0120] In an embodiment of the present invention, the outward expansion flight route is expanded outward from the center of the outward expansion flight route, the scanning point is the outward expansion sample point, the height difference between the outward expansion sample point and the center point of the outward expansion flight route is the set outward expansion flight height, the longitude and latitude of the outward expansion sample point are superimposed according to the height value of the outward expansion flight, and the outward expansion flight altitude is set to 50m. In order to avoid the elevation difference between the outward expansion route and the original terrain, 50m is the height difference between the outward expansion sample point and the center of the outward expansion route. The method for calculating the coordinates and elevation of the outward expansion sample point by the outward expansion sample point is to obtain the distance deviation based on the outward expansion sample point and the center of the outward expansion flight route, determine the heading distance of adjacent outward expansion sample points based on the distance deviation, and determine the longitude and latitude of the outward expansion sample point according to the heading distance;

[0121] Assume that the center coordinates of the extended flight route are (X o ,Y o ,H o ), the distance deviation between the outer expansion sample point and the center of the outer expansion flight route is Δd, and the heading distance between adjacent outer expansion sample points is d h , the radius of the earth is R, then the latitude and longitude of the expanded sample point are calculated as follows:

[0122] In latitude Direction: In the longitude λ direction:

[0123] in is the latitude of the extended sample point, λ sp Longitude of the extended sample point The center latitude of the extended flight route, λ o The center longitude of the outward flight route.

[0124] In an embodiment of the present invention, the coordinates and elevations of the expanded sample points are calculated by inferring the expanded sample points, and the coordinates of the control sample points of the fixed large sample are calculated by the following method:

[0125] Assume that the coordinate set of the external expansion sample points is Based on the coordinates of the expanded sample points, spatial fitting is performed and a plane equation Z=aX+bY+c is fitted using the least squares method, where a, b, and c are fitting parameters. The fitting parameters are calculated using the following formula:

[0126]

[0127] Solve for a, b, c;

[0128] Let the fixed large sample plot be rectangular, with length L and width W. Taking the coordinates (X c , Y c , H c ) of the center point of the large sample plot on the fitting plane as the reference, in the long direction i, m control sample points are evenly set at equal intervals, and the interval In the wide direction j, k control sample points are evenly set at equal intervals. Then the coordinates (X cpij , Y cpij , H cpij ) of the control sample point in the j-th row and i-th column are calculated as follows:

[0129] X cpij = X c + (i - 1)ΔX1, where i = 1, 2,..., m

[0130] Y cpij = Y c + (j - 1)ΔY1, where j = 1, 2,..., k

[0131] H cpij = aX cpij + bY cpij + c

[0132] During the on-site layout process, pole marks are set at the control sample points through GNSS instruments to obtain the coordinates of the pole marks. The number of pole marks set should be no less than 3 to be used for the layout of the sample plot. The layout direction of the pole marks is set to be the same as that of the fixed sample plot to ensure the layout efficiency. The number of pole marks set is no less than 3, which is convenient for determining the layout position of the fixed sample plot.

[0133] In the embodiment of the present invention, it also includes the calculation of the accuracy of the pole mark layout coordinates. The accuracy value between the pole mark coordinates and the control sample point coordinates is calculated through the extended sample points and the extended aerial flight routes. The calculation method is as follows: Let the theoretical coordinates of the control sample point deduced through the extended sample points be (X ct , Y ct , H ct ), the actual pole mark coordinates be (X rb , Y rb , H rb ), and the accuracy threshold be ε. Then the accuracy value P is calculated as follows:

[0134]

[0135] When P < ε, the pole mark layout meets the accuracy requirements; when P > ε, the pole mark layout does not meet the accuracy requirements.

[0136] The embodiments disclosed in the present invention are preferred embodiments, but are not limited thereto. Those of ordinary skill in the art can easily understand the spirit of the present invention based on the above embodiments and make different extensions and changes. However, as long as they do not depart from the spirit of the present invention, they are within the protection scope of the present invention.

Claims

1. A method for deploying forest permanent large sample plots using RTK, characterized in that, It includes the following steps: S1. Point cloud data acquisition; Using a drone equipped with a lidar, airborne laser scanning is carried out above the original terrain, high-precision position coordinates are obtained through GNSS / INS combined positioning, and point cloud data is collected; S2. Three-dimensional digital model establishment; Processing the point cloud data, filtering out noise points and removing isolated points, and generating a digital elevation model based on the processed point cloud data, and generating the center points of forest fixed large plots in the digital elevation model; S3. Field plot planning; Dividing the layout area of the fixed large plot into several field plots, and determining the field plot range and the center coordinates of the field plot; S4. Obtaining the center of the drone's extended flight route; Based on the center coordinates of the field plot, the center coordinates of the drone's extended flight route are inversely calculated, and the extended flight altitude is determined according to the area of the field plot and the height difference between the highest point and the lowest point in the fixed large plot, and an extended flight route is generated; S5. Obtaining the coordinates and elevations of the extended sample points; Through the drone, airborne laser scanning is carried out along the extended flight route to obtain the coordinates and elevations of the extended sample points; S6. Laying out the fixed large plot; According to the coordinates and elevations of the extended sample points, the coordinates of the control sample points of the fixed large plot are inversely calculated, and the fixed large plot is laid out according to the control sample point coordinates.

2. The method for deploying forest permanent plots using RTK according to claim 1, wherein: In S2, the drone is a fixed-wing drone, and the digital elevation model generation method is based on the construction of a Triangulated Irregular Network (TIN). The generation method is as follows: First, perform Delaunay triangulation on the processed point cloud data to construct an irregular triangular network. Let any triangle T in the triangular network abc , and the coordinates of its three vertices are A(x a , y a , z a ), B(x b , y b , z b ), C(x c , y c , z c ). For any point P(x, y) inside the triangle, its elevation z p is calculated by linear interpolation: By traversing the entire triangular network area, a digital elevation model is generated; Setting digital ground points and the plane of the digital ground elevation model on the point cloud data, corresponding the digital ground points to the plane of the digital ground elevation model, establishing a digital ground model, thinning the original terrain in the three-dimensional digital model in the area, and adding the digital height of the digital ground model to the digital ground points in the thinned digital ground model.

3. A method for laying out forest permanent large sample plots using RTK according to claim 2, characterized in that: The method for field plot planning of the three-dimensional digital model of the original terrain and inversely calculating the center of the drone's extended flight route through the center of the field plot includes: a1. Based on the three-dimensional digital model, dividing the layout area of the fixed large plot into several field plots; a2. Determining the field plot range of the field plot, and the field plot range includes the area of the field plot and the field plot boundary; a3. Obtaining the center coordinates of the field plot; a4. Based on the center coordinates of the field plot, obtaining the center coordinates of the drone's extended flight route.

4. A method for laying out forest permanent large sample plots using RTK according to claim 3, characterized in that: The method of calculating the center of the UAV's extended flight path by back-calculating from the center of the on-site sample plot is based on the on-site sample plot range. The center of the fixed large sample plot is obtained by back-calculating from the center coordinates of the on-site sample plot. Let the center coordinates of the on-site sample plot be (X si , Y si ), where i = 1, 2, …, m, and m is the number of on-site sample plots. First, calculate the average value of the center coordinates of all on-site sample plots in the X direction and the average value in the Y direction The formulas are as follows: Then, using as the initial estimate value, based on the coordinates of the extended sample points (X spj , Y spj , H spj ), where j = 1, 2, …, n and n is the number of extended sample points, weighted least squares adjustment optimization is performed; Let the distance from the externally extended sample point to the initial estimated center be Construct an error equation: where v Xj and v Yj are the observation residuals in the X and Y directions, δX and δY are the corrections of the center coordinates in the X and Y directions, and Δd j is the correction of the distance from the outer expansion sample point to the optimized center coordinate; According to the least squares principle Construct the normal equation to solve for δX and δY, where p j is the weight of the extrapolated sample points, determined according to the precision of the points, and construct the normal equation to solve for δX and δY; Final fixed large sample plot center coordinates (X c , Y c ) The calculation formula is: Overlaying the center coordinates of the fixed large plot with the three-dimensional digital model of the original terrain to obtain the center coordinates of the extended flight route; Let the central coordinates of the on-site plot be (X si , Y si ). If the offset of the center of the fixed large plot relative to the center of the on-site plot in the X direction is ΔX and the offset in the Y direction is ΔY, then the calculation formula for the central coordinates (X o , Y o ) of the UAV's expanded flight path is as follows: X o = X si + ΔX,Y o = Y si + ΔY.

5. A method for arranging forest permanent large sample plots using RTK according to claim 4, characterized in that, The method for determining the drone's extended flight altitude through the center of the extended flight route is to determine the height difference between the highest point and the lowest point in the fixed large plot based on the actual altitude of the plot, determine the distance from the lowest point to the center of the extended flight route according to the height difference, and set the flight altitude as the elevation value from the lowest point in the fixed large plot to the center of the extended flight route according to this distance; Let the central coordinates of the on-site plot be set, and the elevation of the highest point in the fixed large plot be H max , and the elevation of the lowest point be H min . The horizontal distance from the lowest point in the plot to the center of the expanded flight line is D, then the expanded flight height H flight The calculation formula is: where ΔH = H max - H min .

6. A method for laying out forest permanent large sample plots using RTK according to claim 5, characterized in that: The extended flight route extends outward from the center of the extended flight route, and the scanned points are the extended sample points. The method for obtaining the coordinates and elevations of the extended sample points includes: b1. Using the center coordinates of the fixed large plot to overlay with the three-dimensional digital model of the original terrain, aligning the center coordinates of the fixed large plot with the extended flight route to obtain the extended flight route; b2. Use a drone to fly at the elevation flight altitude according to the extended flight route for extended scanning to obtain the coordinates and elevation of the extended sample points. Let the central coordinates of the extended flight line be (X o , Y o ). The offset distance of the extended sample point relative to the center of the extended flight line in the X direction is ΔX spj , and the offset distance in the Y direction is ΔY spj . The radius of the earth is R, and the angle between the line connecting the extended sample point and the center of the extended flight line and the due east direction is θ. Then the calculation formula for the coordinates (X spj , Y spj ) of the extended sample point is as follows: X spj = X o + ΔX spj cosθ - ΔY spj sinθ Y spj = Y o + ΔX spj sinθ + ΔY spj cosθ where ΔX spj and ΔY spj are determined according to the set outer expansion distance interval and the outer expansion sample point number. Let the outer expansion step size be d, and the number of the f-th outer expansion sample point in this direction be f. Then ΔX spj = f × d × cosθ, and ΔY spj = f × d × sinθ; Let the external extended flight altitude be H flight , and the elevation of the center of the external extended flight route be H o . The vertical height deviation of the external extended sample point relative to the center of the external extended flight route is ΔH spi . Then the elevation H spi of the external extended sample point is calculated as follows: H spi = H o + H flight + ΔH spi .

7. A method for deploying forest permanent large sample plots using RTK according to claim 6, characterized in that: The method for obtaining the coordinates and elevation of the extended sample points includes: obtaining the distance deviation based on the extended sample points and the center of the extended flight route, determining the course distance between adjacent extended sample points based on the distance deviation, and determining the longitude and latitude of the extended sample points according to the course distance. Let the central coordinates of the extended flight line be (X o , Y o , H o ). Let the distance deviation between the extended sample point and the center of the extended flight line be Δd, and the course distance between adjacent extended sample points be d h . Given the radius of the Earth as R, the longitude and latitude of the extended sample point are calculated as follows: In the latitude direction: In the longitude λ direction: Among them is the latitude of the extrapolated sample point, λ sp is the longitude of the extrapolated sample point is the latitude of the center of the extrapolated flight line, λ o is the longitude of the center of the extrapolated flight line.

8. A method for deploying forest permanent large sample plots using RTK according to claim 7, characterized in that: Based on the coordinates and elevation of the extended sample points, the coordinates of the fixed large sample plot control sample points are calculated by inversion, and the calculation method is as follows: Let the set of coordinates of the extrapolated sample points be Perform spatial fitting based on the coordinates of the extrapolated sample points, and use the least squares method to fit a plane equation Z = aX + bY + c, where a, b, and c are fitting parameters, and the fitting parameters are calculated by the following formula: Solve for a, b, and c. Let the fixed large sample plot be rectangular, with length L and width W. Taking the coordinates (X c , Y c , H c ) of the center point of the large sample plot on the fitting plane as the reference, m control sample points are equally spaced in the length direction i, with an interval of In the width direction j, k control sample points are equally spaced, Then the coordinates (X cpij , Y cpij , H cpij ) of the control sample point in the j-th row and i-th column are calculated as follows: X cpij = X c + (i - 1)ΔX1, where i = 1, 2, …, m Y cpij = Y c + (j - 1)ΔY1, where j = 1, 2, …, k H cpij = aX cpij + bY cpij + c Set up a pole marker at the control sample point through a GNSS instrument to obtain the coordinates of the pole marker, and set the layout direction of the pole marker to be consistent with the fixed sample plot.

9. A method for deploying forest permanent plots using RTK according to claim 8, characterized in that: It also includes the calculation of the coordinate accuracy of the pole marker layout. The accuracy value between the coordinate of the pole marker and the coordinate of the control sample point is calculated through the expanded sample points and the expanded aerial flight routes. The calculation method is as follows: Assume that the theoretical coordinate of the control sample point deduced from the expanded sample points is (X ct , Y ct , H ct ), the actual coordinate of the pole marker is (X rb , Y rb , H rb ), and the accuracy threshold is ε. Then the calculation formula for the accuracy value P is: When P < ε, the layout of the pole marker meets the accuracy requirements; when P > ε, the layout of the pole marker does not meet the accuracy requirements.

Images